Optimization method to obtain appropriate spacing parameters for crop cultivation

Considering the time-consuming and tedious work of the current methods to control plant layout, which is mostly based on expert experience or field trials, we propose an algorithm to optimize and simulate a planting layout based on a virtual plant model and an optimization algorithm. A functional-structural plant model, which combines the structure and physiological function of plants, is used to construct a planting scene. The planting and row spacing are set as the genetic factors and the chromosomes of the genetic algorithm are encoded with a binary method. The photosynthetic yield of the unit planting area is denoted as the fitness value. By using this method, the intercropping of maize and soybean plants and the sole cropping of rice plants are studied. Experimental results show that the proposed method can obtain a high yield planting plan.


Introduction 
In China, many kinds of crops are grown. Different types of crops need different environments to grow well. For example, the distribution and layout of plants in a field affect the light utilization of the crop canopy, the transport of materials, the resilience, and the competition for soil nutrients among the plants. A reasonable arrangement of row spacing and field planting modes can take full advantage of the physiological characteristics of crops and their regulation ability and reduce competition among plants, thereby avoiding overgrowth in early stages, ensuring the full utilization of light, heat, and other resources, and improving field microcirculation to accumulate additional organic matter. Therefore, the scientific and effective control of plant layout is necessary to regulate the growth of crops and improve their yield and quality [1] .
The optimal allocation of plant layout has been valued by agricultural producers and researchers [2] . The determination of the layout of traditional agricultural crops is mostly based on expert experience [3,4] or field trials [2][3][4][5][6][7][8][9][10][11] . Although the empirical model of experts is simple and practical, it tends to target certain characteristics of an area such as climate, geography, soil type, and crop planted, and therefore has less general applicability. Traditional field trials also target a limited number of varieties. Testing should be repeated after plants are changed, thereby resulting in time-consuming and tedious work.
With the popularization of mechanical sowing and planting, the spacing of sowing is often uniformly determined for the same crop, and cannot reflect the characteristics of different plant types. This traditional approach has failed to meet the requirements of modern precision agriculture for the precise positioning, quantification, and timing of agronomic measures and the modern agriculture requirements for smart farming technology [12] and production processes.
Plant spacing, row spacing, and field planting patterns are key factors affecting the regulation of crop spatial layout and population structure. Different combinations of those factors will produce different yields. Therefore, problems related to spatial layout optimization in plants are complex. Many solutions are available for this kind of optimization problem in practical agricultural production, which has a large solution space and involves spatial information. In recent years, there are researchers who utilized functional-structural plant models (FSPMs) and optimization algorithms for this kind of problem and good results have been obtained. For example, Qi et al. [13] used particle swarm optimization (PSO) to automatically optimize the sink strength parameters in a virtual corn model and obtained the ideal plant type of corn to increase the weight of corn, leaves, and stalks. Quilot-Turiona et al. [14] used genetic algorithms to automatically optimize the six most influential parameters in a virtual peach tree model and obtained the ideal peach plant type with improved fruit quality. However, their optimization factors are abstract and focused on describing the physiological processes in crops. Providing intuitive reference information to quantify the plant type design is difficult. Drewry et al. [15] used the multi-objective optimization algorithm to change the five characteristic parameters in the canopy to find the optimal combination of canopy structural parameters, and simultaneously improve plant yield, water use efficiency, and photosynthesis. However, the parameter to be optimized is a macroscopic description of the canopy structure, which cannot describe the geometric morphological characteristics of individual plants.
In this article, we propose an optimization strategy for plant spatial distribution based on a virtual plant model and a genetic algorithm. By using this method, we will study the optimization of plant layouts of two kinds: 1) intercropping of maize and soybean plants; 2) sole cropping of rice plants. To achieve this goal, starting with the solution based on intelligent computing, we combine the geometrical characteristics of the virtual plant models with the spatial structure and physiological characteristics of the plant population, and also the physical and geometric characteristics of the leaves. Our study will provide a technical reference to determine the appropriate spacing parameters for crop cultivation.

Abstraction of the problem
In agricultural practice, two methods have been developed to grow crops in the same area simultaneously: 1) planting one type of plant, or sole cropping and 2) planting different crops at intervals of a certain number of rows, such as 1:2 or 2:3, called intercropping. Intercropping is frequently used between high and short crop plants, such as maize and soy.
Here we consider crop varieties, field configuration, plant type, geographical latitude, planting season, and many other factors and to obtain the optimum morphological traits of crop plants, the optimum plant spacing, and the optimum row spacing parameter characteristics to obtain the maximum photosynthetic yield based on the functional structure model of plants under the specific light conditions. Before designing the algorithm, we should analyze the evolution of plant morphological characteristics of different plant species in different growth stages with different row spacing and planting patterns and the spatial and temporal variations of light distribution. In a monoculture or intercropping mode and different planting modes, such as equal row spacing, wide-narrow row spacing, and staggered row, the range of row and plant spacing should be explicit. The characteristics of monoculture and intercropping methods should be compared through studies of plant spatial layout, and the types of parameters that must be optimized should be determined. For a monoculture, the parameters to be optimized are plant spacing, row spacing, and their relative positions in the field.
For two different plants, namely, A and B, individual plant spacing and row spacing, the allocation of plant and row spacings between the two crops, and the ratio between the number of the rows of two crops should also be optimized. Assuming that f(x) represents the mapping from the combination x of different layout factors in the photosynthetic carbon assimilation of a plant; we may initially derive an optimization model describing the planting layout of a plant as follows: where, x i refers to a vector composed of layout factors, x 1 represents plants planted as a monoculture; x 2 indicates two kinds of plants planted in an intercropping manner. We want to maximize M(x) by using the functional-structural plant models and an optimization algorithm.

Establishment of functional-structural plant models (FSPMs)
A functional-structural plant model is a model that accurately reflects the interaction between the three-dimensional morphological structure and physiological functions of a plant based on its specific morphology. It provides a new way to study the operation of the plant system under different environmental or internal factors. To establish an FSPM model, one needs to establish the plant structure in terms of basic units, the rules of morphological development and the models of the metabolic processes that drive plant growth [16] . We utilize a software called GroIMP [17] (Grammar-related Interactive Modeling Platform Growth) to build the FSPMs we need. It is a modeling platform, based on Java and the tailored modeling language XL (eXtended L-systems), and includes model construction, visualization, interactive modules, etc. XL is developed based on the Java programming language, and further integrated with the parallel rewriting rules of the L-system [18] .

Modeling of plant structures
To simulate the topological structures of the plants and their growth process, the spatial characteristic of maize, soybean and rice plants at each growth stage, which include the topology and the sequence of activities of various plant modules, are analyzed to deduce the parameters of an eXtended L-system, such as axiom and productions. We divide the architecture of those three plants into base, growth unit, fruit, leaf, etc., and suppose that the growth unit is made of an internode, an axillary meristem, and a leaf. On this basis, we can deduce the axiom and productions to describe the spatial structure.
For example, through analyzing of the topological structure of soybean, the initial graphics structure axiom ω of the soybean generating structure and production rule set {R 1 , R 2 , R 3 , R 4 } of an eXtended L-system are as follows: where, B is the base, m is the meristem; P is the growth unit; I is the internode; L k1 is the leaf; G is the pod; & (180) is used to control the angle between soybean pod and stem; RL(α) is used to control the angle of stem and branches; φ(60) controls the trifoliolate leaves; parameters i and j presents the tiller number and internode number; maxRank is the maximum number of single stems, and maxOrder is the maximum number of internodes.
We adopt a method of polyhedral joints to simulate the stems and develop a function with four parameters (the spatial position, the direction, the radius, and the length of cylinder) to conveniently draw a cylinder. Since the stem is cylindrical, an internode was used as a unit in the model, and the stem consists of multiple connective internodes. The cylinder function provided by GroIMP was used to simulate the internodes. Simulation of the leaves in this paper was achieved by controlling of feature points and forming triangle meshes using GroIMP polygon Mesh functions. The male and the female tassels of a maize plant, and the pods of a soybean plant, are simulated based on the analysis of their shape and characteristic by using the suitable functions provided by GroIMP. The parameters of different organs at different growth times are saved into a database, and then the morphological models of those three plants are established, as shown in Figure 1

The models of the metabolic processes
The light environment, the radiation model, the photosynthesis model, the carbon partitioning model, and the growth functions for organs, are considered to build the plants' functional models while neglecting the effects of water and mineral uptake on plants.
We establish a virtual light environment [19,20] , which has two kinds of radiations: 1) the solar radiation directly reaching the Earth's surface, which named the direct radiation; and 2) the solar radiation being scattered after passing through the atmosphere, which called the diffuse radiation. According to the solar constant, the time-varying Earth-Sun distance, the real Earth-Sun distance and the elevation angle of the sun, the direct radiation (W/m 2 ) and the diffuse radiation (W/m 2 ) can be calculated.
The amount of Photosynthetic Active Radiation (PAR), intercepted by the canopy of plants is calculated using a radiation model in the platform of GroIMP [21] . For example, we divide the maize canopy into two parts: 1) the upper part, which is above the top of the soybean plant; and 2) the lower part, which is the rest. With the leaf area index of the upper and lower parts of the maize canopy, together with the unit area of the soybean calculated, radiation intensity of the upper and lower parts of the maize and the soybean canopies is calculated using the radiation. We calculate the crops' extinction coefficients [22] respectively. Then the light interception amount of the intercropping soybean and maize can also be calculated [23] .
A photosynthesis model which changed the species-specific parameters and the environmental parameters accordingly [24] was established to estimate the current average leaf photosynthetic rate (mol/m 2 s). Through the calculation of light distribution in the canopy and single leaf photosynthesis efficiency, the available photosynthetic products from all leaves are collected to a Common Assimilate Pool (CAP) [25] . The accumulation of photosynthetic production PR(T) (unit: g), which is produced in each growth day of soybeans or maize's leaves with the area of x (unit: m 2 ), is dynamically added to CAP, and the assimilation substance Y CAP (T) in day T is [28] : Growth is based on the source-sink hypothesis with organ sink forcing the current assimilates to be distributed to various organs for growth and development. During the growing process, the sink strengths of all organs are added up. In the case of no material transport resistance, each growing organ O obtains the assimilation production K o (T) (unit: g) from the pool through competition, according to their percentage of its own sink strength to the total sink strength [26] . Thus, the biomass amount According to allometric relationships between the biomass increment of organ O in growth cycle T and plant's 3D geometrical sizes, the diameter variation D O (T) of organ O is: where, φ O is the conversion factor between biomass and organ size [26] . To different kinds of plants or different kinds of organs in a plant, the values of φ O are different. However, in order to reduce the complexity of the model, we suppose that the values of φ O are the same for one type of organ in a plant.
In this paper, a day was considered as the basic time unit for the growth cycle. Plants grow once in each simulation step, and their intrinsic function and three-dimensional structure change once. Therefore, the growth dynamics of the maize/soybean plant before the termination of the plant growth are simulated.
The functional-structural model of rice used in this study is based on our previously published paper [28] .

Plant spacing optimization based on the optimization algorithm
With the preceding analysis, an optimization algorithm of a planting layout based on an evolutionary algorithm is proposed. The optimum design of plant layout factors is achieved by optimizing the growth conditions. Several problems should be considered during the algorithm design: the coding strategy and initial population design of individual chromosomes, the design strategies of genetic operators, and the evaluation of individual fitness. The coding format of an individual chromosome is {c 0 , c 1 ,…, c nc-1 }, where c i denotes one chromosome, and nc represents the number of chromosomes contained in the genome. The plant spacing, the row spacing, the relative position, the planting method (monoculture or intercropping), and the ratio between the number of rows are considered as chromosomes. Furthermore, chromosomes comprise multiple genes, and ng refers to the total number of genes in a single chromosome. The genetic code string length, population size, termination algebra optimization, and other parameters strongly influence algorithm efficiency and optimization results, so the optimization objective, algorithm efficiency, and computer performance are combined to set the parameters of the optimization algorithm. Figure 2 shows the flowchart of the plant spacing optimization algorithm. In the algorithm, parents continuously generate individual offspring through selection, crossover, and mutation operations, and the resulting individual offspring compete with their parents to join the updated population. The process continuously cycles until it meets the programmed termination condition. After several generations of optimization occur, the elite individuals are obtained and the appropriate spacing parameters for crop cultivation are also determined.

Coding scheme
In the plant spacing optimization algorithm, the genetic factors are a combination of plant and row spacings. In this study, individuals are coded in terms of their plant and row spacing, and each information bit of the genetic factor is 0 or 1. The range of a parameter is [U min , U max ], and the genetic code length of the parameter is l. The smallest parameter code 0 corresponds to the parameter U min , the largest parameter code 2 l -1 corresponds to the parameter U max , and the difference between two consecutive parameters is the numerical accuracy. As the parameters are evenly distributed, it is easy to calculate the numerical accuracy: Binary genetic factors can be converted into real numbers through the binary code decoding function to obtain the real values of the plant spacing factor G Dspacing and the row spacing factor G Hspacing . If a binary genetic code string is X:a l a l-1 …a 2 a 1 , then the numerical value x is expressed as follows:

Initialization of the population and genetic manipulations
Determining the number of optimal solutions is often difficult when the range of feasible solutions is large. Therefore, we need to conduct random sampling uniformly in the solution space of a problem to obtain individuals [29] . An individual with a length of LR is binary coded, and each bit of information on the binary code is randomly selected from {0, 1}. For population initialization with n c individuals, at least LR×n c random decisions are required.
Genetic manipulation includes crossing, mutation, and selection, which can simulate the natural selection of species in nature. We use the roulette algorithm [30] to implement a selection operation based on the fitness of an individual. The mutation strategy and the cross-operation proposed by Li et al. [29] are used to generate new individuals.

Design of fitness function
The establishment of a comprehensive evaluation model of individual fitness with reasonably high efficiency and low cost is a key step to design an optimization algorithm. The yield usually reveals whether an agronomic measure is appropriate or optimal in agricultural production. In this study, the yield per unit field planting area of the plant population is considered during the design of the fitness evaluation function. The individual fitness evaluation model for high-yielding crops is designed through the observation and analysis of the configuration of different row spacings of high-yielding plant groups and their morphological changes in different growth stages. In addition, the model must be combined with the plant growth cycle so that the optimized parameters can withstand continuous adaptation tests throughout the growing season.
The correlation between crop yield and photosynthetic carbon assimilation is significant, that is, the more the carbon assimilation per unit area of land is, the higher the plant yield will be. As such, the fitness function of the optimization algorithm should be designed on the basis of the amount of assimilated photosynthetic carbon. To accelerate the algorithm's search for the optimal solution, we propose the addition of a dynamically changing penalty function to the fitness function. An appropriate initial value of the penalty factor is set so that the search process can be simplified in the early evolution and that the algorithm can effectively search the solution space. During the algorithm process, the penalty factor constantly changes as the evolution algebra increases to separate the feasible and infeasible solutions effectively.
Yield is one of the most important criteria for measuring the suitability of plant spacing and row spacing in agriculture. Fitness is calculated on the basis of the functional-structural model of the plants [31] .
The fitness M of an individual is the total photosynthesis field per unit planting area, and is calculated as: where, N is the number of plants; Y (i) is the yield of the i-th plant, and s is the planting area.
Assuming that the plant and row spacing are inherited factors, we determine the plant and row spacing in terms of the yield per unit area. Before using genetic algorithms for optimization, we should set relevant parameters, including crossover probability p c , mutation probability p m , initial population size n c , and maximum evolutionary generation t max . We also determine the range of plant and row spacing and generate P 0 .

Genetic algorithm parameters
In the genetic algorithm, L R , n c , p m , p c , and t max of the genetic code string strongly influence computational efficiency and optimization results. According to the range of parameters proposed by [32], the parameters of the genetic algorithm are set on the basis of the optimized object, the efficiency of the algorithm, and the performance of the computer used in the experiment.
(1) L R . In this study, the plant and row spacing of rice were used as genetic factors. The length of the coding string is related to the range of values and optimization precision. Therefore, the coding length of the genetic factors is set on the basis of the experience and plant type of traditional agriculture.
(2) n c . A large population often contains extra genotypes and can simultaneously deal with different individuals; therefore, obtaining the global optimal solution is easy. However, the computation time of each generation is extended as the population size increases. Considering the range of optimization objects and the optimization accuracy, we set the population size to 20.
(3) p c . The convergence rate initially improves as p c is increased from zero.
Considering that the selected elite individuals in this algorithm have retained most of the good genes and that the influence of p c on the good gene structure is reduced, we set p c to 0.9.
(4) p m . Variation operations can effectively generate new individuals and prevent the genetic algorithm from degrading to a random search [33] . In this study, p m is set to 0.08.
(5) t max . If the optimal solution to a problem is known, then finding the optimal solution is the corresponding termination condition. However, if the optimal solution to a problem is unknown, then an iteration limit constant should be set as the termination condition [32] . The algorithm takes a long time because the photosynthesis yield per unit plant area taken as the fitness value takes a long time to compute; therefore, we set the iteration limit constant to 100 generations.

Experiments and analysis
The experiment is based on a GroIMP platform and is programmed in Java. Our computer hardware is Intel (R) Core i5 4590 with 2.66 GHz CPU, 8GB memory, and NVIDIA GTX 970 graphics card. The rice variety used in the experiment is double haploid rice (DH) [26] .

Optimization of intercropping planting of maize and soybean plants
The inter-planting mode is a mode that intercropping crops are planted in the same field in a certain proportion. There has been a great deal of maize and soybean intercropping pattern research in field trials. Based on the establishment of a maize/soybean intercropping function and structure model, the focuses of this study are the inter-planting patterns with ratio 2:3 and 1:2 of the maize and the soybean plants, to study the intercropping and compact planting, as shown in Figure 3. The distance between maize and soybeans, soybean and maize is defined as inter distance D in , the distance between maize and maize, soybean and soybean is defined as space distance D MM and D SS . In this algorithm, more light interception is the optimization target. Through continuous iteration of the genetic algorithm to optimize the value of D in , D MM and D SS within their specified ranges, the best distance combination which maximizes light interception will be obtained.
In the experiments, the ratios 2:3 and 1:2 were used as the maize/soybean intercropping mode for planting distance optimization. Furthermore, the intercropping space distances, D in between maize and soybean population is within the range [5 cm, 60 cm], and the simulated intercropping population was with 10 columns. As for the ratio 2:3, after many minutes running for growth, within the domain of D in and D MM or D SS , the optimal planting distance combination for optimized light interception can be calculated simulating unit land area, shown in Table 1. The change of light interception quantity in 120 days' growth process of initial plants and optimal plants are shown in Figure 4. After 110 days, the light interception is basically stable.
The light interception before optimization was 7416 mol/m 2 · d, while after optimization it is 9424 mol/m 2 · d, which is about 1.27 times more. Figure 5 shows the visualization simulation results of the soybean/maize population. As for the ratio 1:2, after many times running for growth, within the domain of [5 cm, 60 cm], the results are shown in Table 1. The change of light interception quantity in 120 days' growth process of initial plants and optimal plants are shown in Figure 4. The light interception before optimization was 7639 mol/m 2 · d, while after optimization is 10429 mol/m 2 · d, which is about 1.36 times more. The results show that the combination of virtual plant model, the genetic optimization algorithm and the ability of the computer to process data at high speed leads to obtain the best combination of planting distances in the designated domain and under the designated inter-specific mode.  Figure 4 Leaves' light interception before and after optimization Figure 5 Population morphology of maize/soybean (2:3)

Optimization of sole cropping of rice plants
In traditional agriculture, the large distance between adjacent plants is generally defined as row spacing, and the small distance is referred to as plant spacing. Most of the field trials have focused on narrow row spacing (15-33 cm), while some studies have investigated broad row spacing between 36 and 42 cm [34] . In most modern mechanized agricultural machinery, row spacing and plant spacing are fixed at 30 cm and at 12-16 cm, respectively. In this experiment, the changes in the plant spacing (D spacing ) and the row spacing (H spacing ) are controlled within a certain range. For example, D spacing and H spacing belong to the range of [10 cm, 30 cm] and the range of [15 cm, 45 cm], respectively. The initial row spacing is set to 20cm, and initial plant spacing is set to 15 cm.
Rice plantings in 4×4 and 5×5 patterns are constructed on the basis of the GroIMP platform to generate initial populations for genetic algorithms. The construction of virtual scenes and the simulation of the rice yield are based on the functional structure model proposed by [26]. Each experiment is conducted three times. The optimization results of plant spacing, row spacing, and yield per unit area are shown in Table 2. The yield per unit area of the optimized plant population improves compared with the planting conditions of 20 cm row spacing and 15 cm plant spacing. As can be seen in Figures 6a and 6c, the non-optimized rice visualization scene shows a dense population with generally small panicles in the rice plants. The middle rice plants are affected by the light received because of the interlacing of the leaves of the surrounding plants, and their growth is weak. The optimized rice visualization scenes are shown in Figures 6b  and 6d. In the population, the plant spacing and the number of rice ears are large, the overall growth is excellent, and the yield per unit area of the optimized population is improved. a. 4× 4 non-optimized b. 4× 4 optimized c. 5× 5 non-optimized d. 5× 5 optimized Figure 6 Comparisons between different layouts of the rice plants

Discussion
In comparison with traditional plant spacing optimization, the method based on the virtual plant model and the optimization algorithm is a new effort towards modeling, which reduces the duration of an experiment, resolves the disadvantages of traditional field experiments in selecting plant strains, and eliminates space and environmental constraints.
This research achieves the integration of an optimization algorithm with the model systems in both temporal and spatial dimensions. It is an exploration of a new application area of plant functional-structural models [28] .
The virtual plant growth scene constructed in this study can reveal the changes in the morphological structure of plants at different planting distances. The virtual plant growth scene combines physiological processes, such as photosynthesis, assimilation, and respiration. Some data and modeling methods are derived from the previous studies [25,26,31,35] .
In the proposed plant spacing optimization, the photosynthesis yield per unit area is considered the optimization goal. Although a high yield is an important indicator of the reasonability of plant spacing, the quality of rice grains [36] and the lodging resistance capacity of plants are also affected by plant spacing and indirectly or directly influenced by the rice plant. These factors are related to many complex conditions, such as plant traits, air conditions, climate, water, soil nutrients, and inorganic salt, which should be subjected to numerous field experiments. Therefore, the proposed optimization method should be further improved.
In the aspect of virtual model construction, data collection is based on manual measurement with low accuracy, although data regarding the morphological structure of plants used in the modeling of this study were obtained through field experiments.
In the future, data should be collected through image extraction and 3D scanning to enhance the accuracy of plant modeling parameters. In system development, virtual plant construction, canopy leaf area calculation, plant spacing optimization and closure time have been completed.
In future research, the application of graphics hardware (GPU) should be considered to compute the shadows affecting photosynthesis.

Conclusions
A method of plant spacing optimization based on a genetic optimization algorithm was proposed in this study. First, the virtual plant model, which combines the structure and physiological function of crop plants, was used to construct a planting scene. The planting strategies with a higher yield were obtained for the intercropping of maize and soybean plants or the sole cropping of rice plants by the optimization algorithm respectively.
This study provides new ideas for researchers studying the rational close planting of crops. In view of the fact that the optimization algorithm of planting distance specifies the pattern and the planting distance change domain, we still need to explore a more suitable optimization algorithm for intercropping planting distance, especially the optimization algorithm in a specific planting pattern.