An Integrated Yield Prediction Model for Greenhouse Tomato

Abstract

The commonly used greenhouse crop yield prediction models today have their specific application scenarios, which may not ensure the accuracy of the results if the greenhouse environment changes. This greatly restricts their use in the greenhouse environment. To solve this problem, two widely used tomato growth models were compared in the study: TOMGRO and Vanthoor, and then an integrated model was obtained. Through the extended Fourier amplitude sensitivity test (EFAST), the model parameters were divided into three categories: optimized, fixed and ignored. In addition, Bayesian optimization was used as an optimization algorithm, through which the parameters applicable to the greenhouse can be optimized based on the greenhouse data. Compared with TOMGRO and Vanthoor, the output of the integrated model was more reasonable and universal, and the RMSE in the integrated model was 2.5974 while that in TOMGRO and Vanthoor both were over 17, reflecting the fact that the model output was closer to the actual value. According to the verification results of four-year greenhouse data, the model had high performance in predicting yield.

1. Introduction

Crop growth model is an essential part for the optimization of the cultivation management. Besides, it can also benefit the environment control. As we all know, the greenhouse environment is a complex system with multi-variables, nonlinearity, strong coupling and large inertia. The crop growth greatly affects the environment factors with the transpiration and photosynthesis. However, due to the complex mechanism and numerous parameters, the growth model has few applications in the greenhouse and the several widely used growth models are based on specific application scenarios, which is not applicable if the environment changes. The lack of growth models has led to insignificant increases in crop yields and unsatisfactory control effects in greenhouses microclimate. Therefore, there is an urgent need for a more versatile and applicable crop yield model.

Growth models mainly consist of two categories according to the principle: descriptive model and explanatory model [1,2]. The descriptive model based on the existing theoretical knowledge and practical experience, determines the correlation between the research factors through statistical regression analysis of mass crop data [3]. Based on the dynamics principle, the explanatory model describes the relationship between environmental factors, cultivation management and crop growth, morphological development and yield forming process [3].

The two types of model above have their own advantages and disadvantages, but the explanatory model is more practical to reflect the actual growth process of crops. Moreover, with the deep research on crop physiological mechanisms and the development of computer technology, the explanatory model will be a trend for crop model development. For many years, scholars have conducted extensive researches on the explanatory model of crop growth. Jones et al. [4] proposed the TOMGRO model for the growth and development of tomato, parameters of which were mainly based on crop stems. Heuvelink [5] proposed the TOMSIM model, which focused on the study of crop canopy light interception and referenced the SUSROS87 model for the accumulation of dry matter. Israel and other countries [6] developed the common-use growth model HORTISIM, which provided the effective strategies for greenhouse environmental control and management. Van Keulen [7,8] proposed the SUCROS model, which simulated the weight of dry matter in various organs of wheat using solar radiation as the main factor affecting crop growth. Ritchie et al. [9,10,11,12] proposed the CERES model, which simulated the process of crop growth, soil and water balance, as well as energy balance. This study selected two representative tomato explanatory models that can be used for yield prediction for comparison: TOMGRO and Vanthoor.

TOMGRO is a tomato growth model proposed by Jones et al., which focuses on the study of the relationship between the growth of tomato and greenhouse environmental factors such as solar radiation, temperature and $C{O}_2$ to scientifically manage and predict the whole growth process of tomato. TOMGRO describes the entire model by using seven basic crop physiological factors (leaves number, nodes number, fruit number, leaves dry weight, nodes dry weight, fruit dry weight and leaf area index), and it is the most representative model in current greenhouse tomato growth and development simulation research [13,14,15]. The principle of TOMGRO is based on the source-sink theory. The biggest weakness of TOMGRO is that it is difficult to obtain some key parameters such as the potential growth rate in the real greenhouse environment, which has restricted the scope of its application [16].

Vanthoor is a well-structured explanatory model and its simulation performance has been extensively validated in European and North American greenhouses [17]. Vanthoor is actually a greenhouse environmental system model which includes a greenhouse microclimate model and a greenhouse crop growth model. Its crop growth model uses the buffer theory and basically obtains two crop growth status results (leaf area index and dry matter) [18,19,20]. The buffer theory is more abstract than the source-sink theory, because the former is an artificially added concept. Notably, many of its variables are more affected by microclimate parameters such as temperature or average temperature rather than by crop parameters. This makes it difficult to intuitively reflect some of the physiological characteristics of crops, but relatively, the value of its variables is easy to be obtained in the greenhouse.

It can be seen from the above comparison that TOMGRO is more suitable for cultivation management of crops and Vanthoor is more suitable for environmental regulation of greenhouses. In order to get a greenhouse tomato growth model with a wider applicable scope, the above two model was compared and obtained an integrated model in this study.The model combined the strength of TOMGRO and Vanthoor, that is, it retained a large number of crop physiological parameters, and ensured the capturing of key parameters in the greenhouse as well, greatly improving the practicality of the model. However, the model alone cannot dynamically adjusted parameters according to different greenhouse environments, so combination with the sensitivity analysis and optimization algorithm should be a good way for the repeating use of the model in various greenhouses.

2. Materials and Methods

The greenhouse data was collected from the Experimental Station (318500 N, 1218330 E), of the National Facility Agricultural Engineering Technology Center of China, in Chongming County, Shanghai City, China. Tomato plants, Gorioso, were grown in a 680 m${}^2$, 5 m high Venlo-type glass greenhouse with numerous actuators such as the ventilation equipment (roof-windows, side-windows and fans), heating equipment (heating pumps), curtains (shading screens and energy screens), lighting equipment (high voltage sodium lamps and light-emitting diodes) and irrigation equipment. The greenhouse contains eighteen rows of tomato crop. Crop planting density was 2.5 plants/m${}^2$.

The data collection began on September 20 and ended on December 28 each year (2014–2018). During the period, the side branches of tomatoes were cut off regularly, only the main stem were kept. The harvested fruit was weighed manually and recorded its weight daily. The crop was grown in the organic substrate (Rice husk charcoal), and sufficient irrigation was performed daily with the nutrient solution through the drip. Greenhouse climate data such as temperature (T, ${}^{\circ }\mathrm{C}$), $C{O}_2$ concentration ($C{O}_2$, ppm), photosynthetically active radiation ($PAR$, $\mathsf{\mu }$mol/m${}^2$/s), humidity (H, %) and so on were recorded every hour by sensors.

2.1. Integrated Model

The main difference between TOMGRO and Vanthoor is that the division of the growth stages in TOMGRO is mainly based on the number of plant nodes, while that in Vanthoor is based on the temperature. In addition, the main calculation time scale in TOMGRO is days, while that in Vanthoor is seconds. Though based on different growth theories, TOMGRO and Vanthoor are similar in the structure of computing the net dry matter. By comparing the two models, an integrated model was obtained, which can be summarized as the following five steps:

(1)

Dividing crop organs into different age classes.

(2)

Calculating the amount of dry matter provided by the environment: SUPPLY.

(3)

Calculating the amount of dry matter required by crops: DEMAND.

(4)

Comparing SUPPLY and DEMAND to screen out two situations: oversupply and undersupply.

(5)

Obtaining the final change rate of net dry matter, according to different supply and demand situations.

Based on the five steps, the integrated model is described in detail as below:

Step 1: dividing age classes:

Consistent with TOMGRO and Vanthoor, the integrated model divided crop organs into different age classes that crop physiological characteristics in the same age group were considered the same. The physiological age of the crop took reference to the idea of TOMGRO, which is described by the current number of nodes, in other words, the number of crop organs in a certain age class is related to the number of crop nodes. Taking the age group i as an example, the rate of change in the number of stems (${\dot{N}}_S\left(i\right)$), leaves (${\dot{N}}_L\left(i\right)$) and fruits (${\dot{N}}_F\left(i\right)$) was as shown in Equations (1)–(3):

$${\dot{N}}_S\left(i\right)=\left\{\begin{array}{cc}GEN{R}_S·n_{Plants}-r_L·n_F·N_S\left(1\right)\hfill & i=1\hfill \\ \left(N_S(i-1)-N_S\left(i\right)\right)·r_L·n_F\hfill & 1<i<n_F\hfill \\ N_S(i-1)·r_L·n_F\hfill & i=n_F\hfill \end{array}\right.$$

(1)

$${\dot{N}}_L\left(i\right)=\left\{\begin{array}{cc}GEN{R}_L·n_{Plants}-r_L·n_F·N_L\left(1\right)\hfill & i=1\hfill \\ \left(N_L(i-1)-N_L\left(i\right)\right)·r_L·n_F\hfill & 1<i<n_F\hfill \\ N_L(i-1)·r_L·n_F-P_L\hfill & i=n_F\hfill \end{array}\right.$$

(2)

$${\dot{N}}_F\left(i\right)=\left\{\begin{array}{cc}GEN{R}_F·n_{Plants}·R_c-r_F·n_F·N_F\left(1\right)\hfill & i=1\hfill \\ \left(N_F(i-1)-N_F\left(i\right)\right)·r_F·n_F\hfill & 1<i<n_F\hfill \\ N_F(i-1)·r_F·n_F-P_F\hfill & i=n_F\hfill \end{array}\right.,$$

(3)

where i represented the i-th age group. The newest crop organs were in the first group and the oldest ones in the last group. $GEN{R}_S$, $GEN{R}_L$ and $GEN{R}_F$ are the appearance rate of new nodes, leaves and fruits respectively ($no./plant/d$), $n_{Plants}$ is the plant density ($no./m^2$), $R_c$ is the ratio of supply and demand, $r_L$ and $r_F$ are the development rate function of leaves ($nodes$) and fruit respectively ($1/d$), $n_F$ is the number of organ age classes, $N_S\left(i\right)$, $N_L\left(i\right)$ and $N_F\left(i\right)$ are the number of nodes, leaves and fruits in the i-th age class respectively ($no./m^2$), and $P_L$ and $P_F$ are the leaves and fruit mortality respectively ($no./m^2/d$).

According to Jones et al. (1991), the appearance rate of new nodes was $GEN{R}_S=GENRAT·F_N\left(T\right)·F\left(C{O}_2\right)$, in which $GENRAT$ is a constant depending on tomato varieties. However, Cooman et al. (2006) [21] stated that $GENRAT$ is a linear function related to the number of stems. Given that the data of nodes in the conditions of specific temperature and humidity cannot be obtained in the actual greenhouse, and relevant studies [22] found that $C{O}_2$ inhibition function in TOMGRO had little effect on the appearance rate of the nodes, the integrated model used the calculation equation of the nodes as shown in Equation (4). The appearance rate of fruit in Vanthoor: $GEN{R}_F=a+b·T_{Can}^{24}$ is a linear function related to the daily average temperature. Temperature is a parameter that indirectly reflects crop growth, so it was much more inaccurate than that in TOMGRO. Therefore, the equations of stems appearance rate ($GEN{R}_S$), leaves appearance rate ($GEN{R}_L$) and fruit appearance rate ($GEN{R}_F$) in the integrated model were finally shown in Equations (4)–(6):

$$GEN{R}_S=\left(a_g+b_g·N_S\right)·F_N\left(T\right)$$

(4)

$$GEN{R}_L=\frac{GEN{R}_F}{1+TPL}$$

(5)

$$GEN{R}_F=GEN{R}_S·FPN,$$

(6)

where $a_g$ and $b_g$ are function coefficients of organ appearance rate, $F_N\left(T\right)$ is the temperature inhibition function, $TPL$ is the ratio of new trusses to new leaves, and $FPN$ is the ratio of new fruit to new nodes.

A set of experimental results of leaves and fruit development rates at different temperatures were given in TOMGRO, and the development rate in Vanthoor and Cooman was the linear function related to daily average temperature. Therefore, it was concluded that the linear function of daily average temperature was more accurate when characterizing the development rate, so the development rate of leaves ($r_L$) and fruit ($r_F$) in the integrated model is shown in Equations (7) and (8):

$$r_L=a_{rL}+b_{rL}·T_d$$

(7)

$$r_F=a_{rF}+b_{rF}·T_d,$$

(8)

where $a_{rL}$ and $b_{rL}$ are function coefficients of leaf development rate, $a_{rF}$ and $b_{rF}$ are function coefficients of fruit development rate and $T_d$ is the daily average temperature (${}^{\circ }\mathrm{C}$).

When the leaf area index (LAI) exceeded the critical value (XLAIM) in TOMGRO, the mature leaves fell off (the last age class) and the pruning operation (all age groups) was carried out; when LAI exceeded XLAIM in Vanthoor, the extra leaf carbohydrates were totally subtracted. Therefore, a leaf mortality ($P_L$) equation of leaves was proposed as shown in Equation (9) by combining the leaf mortality of the above two models. It simulated the leaf death and pruning based on age groups and maximum LAI; and the older leaves fell off and were trimmed more possibly. Although TOMGRO proposed the concept of fruit mortality, no specific equation and description were given, and there was no concept of fruit mortality in Vanthoor. Therefore, the study adopted Cooman’s fruit mortality ($P_F$) equation, a linear function of $R_c$, as shown in Equation (10):

$$P_L\left(i\right)=\left\{\begin{array}{cc}{N}_L\left(i\right)·\left(1-\frac{XLAIM}{LAI·k·i}\right)\hfill & LAI\ge \frac{XLAIM}{k·i}\hfill \\ 0\hfill & LAI<\frac{XLAIM}{k·i}\hfill \end{array}\right.$$

(9)

$$P_F=a_{PF}-b_{PF}·R_c,$$

(10)

where $XLAIM$ was the critical LAI ($m^2/m^2$), $a_{PF}$ and $b_{PF}$ were the function coefficients of fruit mortality, and k was the regulatory factor.

Step 2: calculating supply:

Supply represented the amount of carbohydrates produced by crops after absorbing external energy. Accordingly, this step was mainly used to calculate the photosynthesis rate and respiration rate of crops. Since most of the results from the photosynthesis models were similar, the integrated model here was consistent with TOMGRO. Equations (11)–(13) showed the calculation formula of SUPPLY, photosynthesis rate (A) and respiration rate ($R_m$):

$$SUPPLY=E·\left(A-R_m\right)$$

(11)

$$A=D\frac{P_{max}}{K}ln\frac{(1-m)P_{max}+\epsilon ·K·PPFD}{(1-m)P_{max}+\epsilon ·K·PPFD·e^{-\left(K·LAI\right)}}$$

(12)

$$R_m=Q_{10}^{0.1(T-20)}\left(RES{P}_L·\left(W_L+W_S\right)+RES{P}_F·W_F\right),$$

(13)

where E is the conversion efficiency ($g\left\{DM\right\}/g\left\{C{H}_{2}O\right\}$), D was the unit conversion factor ($g\left\{CH2O\right\}/umol\left\{C{O}_2\right\}/3600$), $P_{\max }$ is the maximum photosynthesis rate ($\mu mol/m^2/s$), K is the light extinction coefficient, m is the light transmission coefficient, $\epsilon$ is the quantum efficiency ($mol\left\{C{O}_2\right\}/mol\left\{photons\right\}$), $PPFD$ is the photosynthetic photon flux density ($\mu mol/m^2/s$), $Q_{10}$ was the sensitivity to temperature, $RES{P}_L$ and $RES{P}_F$ are the relative respiration requirement of leaves and fruit respectively ($g\left\{CH2O\right\}/g\left\{DM\right\}/d$), and $W_L$, $W_S$ and $W_F$ were dry matter of the leaves, nodes and fruit respectively ($g/m^2$).

Step 3: Calculating demand:

Demand represented the amount of carbohydrates required for the growth of crop organs. TOMGRO was consistent with Vanthoor in the fruit dry matter demand. Because Vanthoor did not set the age groups of leaves and nodes, it had no formula to calculate the required volume of the leaves and nodes. Therefore, the integrated model was consistent with TOMGRO in the calculation of the leaves and nodes. In addition, the $C{O}_2$ inhibition function had little effect on the results in the step, and there was no $C{O}_2$ inhibition function in Vanthoor, so the function was removed in the integrated model. The dry matter demand of stems ($S_{dem}$), leaves ($L_{dem}$) and fruit ($F_{dem}$) and leaf area change rate (${\dot{A}}_{LP}$) were as shown in Equations (14)–(17):

$$S_{dem}\left(i\right)=\frac{L_{dem}\left(i\right)·FRSTM·N_S\left(i\right)}{N_L\left(i\right)}$$

(14)

$$L_{dem}\left(i\right)=\frac{(1+FRPET)·{\dot{A}}_{LP}\left(i\right)}{SLA}$$

(15)

$$F_{dem}\left(i\right)=POF\left(i\right)·F_N\left(T\right)·N_S\left(i\right)$$

(16)

$${\dot{A}}_{LP}\left(i\right)=N_L\left(i\right)·POL\left(i\right)·F_N\left(T\right),$$

(17)

where $FRSTM$ is the ratio of node growth rates to leaf growth rates, $FRPET$ is the ratio of petiole weight to blade weight, $SLA$ is the average specific leaf area ($m^2/g\left\{DM\right\}$), $POF$ ($g\left\{DM\right\}/fruit/d$) and $POL$ ($g\left\{DM\right\}/leaf/d$) are the potential growth rate of fruit and leaves respectively (Equations (18) and (19)).

The experimental data of the potential growth rate of fruit and leaves at different organ development stages were given in TOMGRO, and the empirical equation at the development stages was used in Vanthoor. They were more precise than the discrete data in TOMGRO, so using a combination of Vanthoor and other related researches [21,23], the potential growth rates of fruit ($POF$) and leaves ($POL$) in the integrated model were as shown in Equations (18) and (19):

$$POF\left(i\right)=0.0458·a_F·e^{-b_F\left(S_{FD}-c_F\right)}·b_F·e^{-b_F\left(S_{FD}-c_F\right)}$$

(18)

$$POL\left(i\right)=1.248·a_L·e^{-b_L\left(S_{LD}-c_L\right)}·b_L·e^{-b_L\left(S_{LD}-c_L\right)},$$

(19)

where $a_F$, $b_F$ and $c_F$ are the coefficients in the potential growth rate function of fruits, $S_{FD}$ is the fruit development stage, $a_L$, $b_L$ and $c_L$ are the coefficients in the potential growth rate function of leaves, and $S_{LD}$ is the leaf development stage.

Step 4 Comparing supply and demand:

In TOMGRO, the supply–demand ratio was compared daily to determine the $R_c$, so that the actual organ dry matter growth rate was obtained, while Vanthoor used the idea of the cumulative amount of difference between supply and demand to determine whether to distribute dry matter, that is, it only carried out dry matter distribution when the supply reached a certain level. In comparison, the processing in TOMGRO was finer, so the integrated model used the idea of TOMGRO to handle the relationship of supply and demand. The supply–demand ratio ($R_c$) and the actual dry matter growth rate of fruit ($g_F$), stems ($g_S$) and leaves ($g_L$) were shown in Equations (20)–(23):

$$R_c=\frac{SUPPLY}{DEMAND}$$

(20)

$$g_F\left(i\right)=\left\{\begin{array}{cc}{F}_{dem}\left(i\right)& R_c\ge 1\\ F_{dem}\left(i\right)·R_c& 0\le R_c<1\end{array}\right.$$

(21)

$$g_S\left(i\right)=\left\{\begin{array}{cc}{S}_{dem}\left(i\right)& R_c\ge 1\\ S_{dem}\left(i\right)·R_c& 0\le R_c<1\end{array}\right.$$

(22)

$$g_L\left(i\right)=\left\{\begin{array}{cc}{L}_{dem}\left(i\right)& R_c\ge 1\\ L_{dem}\left(i\right)·R_c& 0\le R_c<1\end{array}\right.$$

(23)

Step 5: Calculating dry matter:

TOMGRO and Vanthoor had similar structure in calculating the rate of change of dry matter: the dry matter change rate in a specific age group = the actual dry matter growth rate in the current age group + the dry matter coming the previous age group − the dry matter entering into the next age group. So the change rate of the dry matter of stems (${\dot{W}}_S$), leaves (${\dot{W}}_L$) and fruit (${\dot{W}}_F$) in the integrated model were shown in the following Equations (24)–(26):

$${\dot{W}}_S\left(i\right)=g_S\left(i\right)+\left(W_S(i-1)-W_S\left(i\right)\right)·r_L·n_F$$

(24)

$${\dot{W}}_L\left(i\right)=g_L\left(i\right)+\left(W_L(i-1)-W_L\left(i\right)\right)·r_L·n_F$$

(25)

$${\dot{W}}_F\left(i\right)=g_F\left(i\right)+\left(W_F(i-1)-W_F\left(i\right)\right)·r_F·n_F.$$

(26)

See from the above description of the integrated model, the model used the shared structure of the source-sink theory and the buffer theory. The idea of dividing the age classes according to the number of nodes in TOMGRO was mainly adopted. With the rapid development of technologies such as computer vision, TOMGRO used the more intuitive data of nodes to reflect the state of crop growth, providing an interface for subsequent application of computer vision in greenhouse control such as fruit picking and pest warning, while the temperature in Vanthoor can only reflect the growth of crops indirectly. In addition, the integrated model provided empirical equations for parameters that were hard to be obtained based on other relevant models and studies such as Vanthoor, and optimized some of the content in the growth model, such as removing the $C{O}_2$ inhibition function when calculating the appearance rate of the new nodes and optimizing the leaf area processing function, etc.

2.2. Sensitivity Analysis

The integrated model should first find the parameters that have a great impact on the output and change with the greenhouse environment before optimization in a specific greenhouse environment. The main principle of finding the parameters to be optimized is to obtain these parameters by changing the input parameters, and the analysis of the influence degree of the model parameters on the model output generally depends on sensitivity analysis.

The sensitivity analysis is divided into two major categories: local sensitivity analysis and global sensitivity analysis [24,25]. The local sensitivity analysis is mainly applicable to scenes with no influence between the parameters in the model, while global sensitivity analysis is applicable to scenarios where the parameters in the model affect each other, and the related study [26,27] shows that global sensitivity analysis has better effect on crop model parameters analysis than local sensitivity analysis. Most of the parameters in the integrated model are coupled to each other, so the global sensitivity analysis was used in the study.

The global sensitivity analysis is composed of Morris method, sobol method, FAST method, EFAST method [28,29,30,31], etc. Currently, sobol and EFAST method are more commonly used. Because the EFAST method combines the FAST and sobol methods, it has the advantages of stability, high precision and fast calculation compared with other methods [32,33]. Therefore, the EFAST method is used to analyze the sensitivity of the relevant parameters in the integrated model.

According to the data collection table (see Appendix A) of the integrated model in the study, 13 parameters were selected for sensitivity analysis (because the parameters POL and POF which determine the sink-strength were one-dimensional vectors of the crop age class, it was difficult to conduct sensitivity analysis in the study. However, according to related studies [23] and simulations, POF has a great influence on the yield, so it needs to be optimized by default). The parameter range was set according to the average value of the parameters provided by Cooman [21], and it took reference to the +/−10% principle [34,35,36]. In this range, the parameters were evenly distributed [37,38], as shown in

Table 1. Parameters for sensitivity analysis.

Parameters	Description	Ranges	Distribution
$GENRAT$	Maximum rate of node initiation	[0.45,0.55]	uniform
$FPN$	Fruit initiated per new node	[0.45,0.55]	uniform
K	Light extinction coefficient	[0.522,0.638]	uniform
m	Leaf light transmission coefficient	[0.09,0.11]	uniform
$ϵ$	Leaf quantum efficiency	[0.058,0.071]	uniform
$\tau$	Carbon dioxide use efficiency	[0.062,0.076]	uniform
$P_{root}$	Supply of photosynthesis for root growth	[0.063,0.077]	uniform
$S_{C{O}_2}$	Effect of $C{O}_2$ on new stems	[0.00027,0.00033]	uniform
E	Conversion efficiency	[0.675,0.825]	uniform
$RES{P}_L$	Relative respiration requirement for leaf	[0.0135,0.0165]	uniform
$RES{P}_F$	Relative respiration requirement for fruit	[0.009,0.011]	uniform
$r_L$	Leaf development rate	[0.009,0.011]	uniform
$r_F$	Fruit development rate	[0.018,0.022]	uniform

.

2.3. Optimization Algorithm

After the integrated model was determined and sensitivity analysis was conducted, parameters optimization was carried out. It covered three cases: (1) the input and output data corresponding to the parameters to be optimized can be obtained in reality.In this case, the actual data can be directly used for optimization; (2) the input and output data corresponding to the parameters to be optimized cannot be actually obtained, but can be obtained indirectly from other parameters.In this case, the actual data can also be used for optimization; (3) The input and output data of the parameter to be optimized cannot be actually obtained. This requires optimization by soft measurement, such as state-based estimation. The sensitive parameters in the integrated model were not the first case. Taking the $r_F$ parameter as an example, it actually indicated the transfer rate of the number of fruits between the age groups.In TOMGRO, it was obtained in a laboratory environment based on a large number of data such as the fruit dry matter, but it was impossible to obtain them in a real greenhouse. Therefore, the method of indirect parameter optimization was adopted. Take the parameter $a_{rF}$ and other parameters in the equation where $r_F$ is located (see below) as an example.

$$r_F=a_{rF}+b_{rF}·T_d$$

(27)

$${\dot{N}}_F\left(i\right)=\left\{\begin{array}{cc}GEN{R}_F·n_{Plants}·R_C-r_F·n_F·N_F\left(1\right)\hfill & i=1\hfill \\ \left(N_F(i-1)-N_F\left(i\right)\right)·r_F·n_F\hfill & 1<i<n_F\hfill \\ N_F(i-1)·r_F·n_F-P_F\hfill & i=n_F\hfill \end{array}\right.$$

(28)

Therefore, the above equation was calculated to obtain the relationship between the total number of fruits at the current moment and $r_F$.

At specific time t, the total number of fruits was the sum of the number of fruits in all age groups:

$$N_{Ft}=\sum _{i=1}^{n_F}{N}_{Ft}\left(i\right).$$

(29)

From Equations (28–(30) can be calculated:

$$N_{Ft}=N_{Ft-1}+GEN{R}_F·n_{plants}·R_c.$$

(30)

It was found from Equation (30) that $r_F$ disappeared during the calculation, that is, $r_F$ was only an intermediate variable, which had no effect on the total number of fruits at a certain moment. In fact, Equation (30) indicated that the current number of fruits was the number of fruits in the previous moment plus the number of new fruits. This physical meaning can also verify the reason why $r_F$ disappeared.

What $r_F$ really affected was the number of mature fruits.The number of mature fruits was $N_F$ (10) (assuming the total age groups are 10), and it was related to the number of fruits in the previous age group, while the number of fruits in the previous age group was actually not available in the real greenhouse. It can be found that it was difficult to optimize the parameters locally in the model, therefore, the study finally selected the model’s microclimate as input and the mature fruit yield as the output, and used the overall growth model for parameter optimization. There were two reasons for this choice: 1. Microclimate and mature fruit yield data can be obtained in the real greenhouse; 2. all the sensitive parameters for yield can be optimized at one time.

Therefore, the optimized objective function was set as the RMSE between the predicted yield and the actual yield in the whole growth cycle, as shown in Equation (31).

$$f_{BO}=\sqrt{\frac{1}{n}\sum _{i=1}^t{\left(W_M\left(i\right)-W_{MA}\left(i\right)\right)}^2},$$

(31)

where $W_M$ was the simulated yield, $W_{MA}$ was the actual yield, n was the simulated days.

Take the commonly used optimization algorithm particle swarm optimization (PSO) as an example, it took more than one minute for the integrated model to run once (it simulated crop growth for 100 days on a i7-7700HQ PC with 8 GB RAM), and the time multiplied the iterations and populations quantity to get a very long optimization time. This problem was fundamentally a costly optimization problem, that is, it took a long time to calculate the objective function. Usually, the solution to such problem in the field of evolutionary computing is to use a surrogate model, or Bayesian optimization. Bayesian optimization is often used to calculate the target value that is costly and time-consuming [39]. It is equivalent to introducing domain knowledge during the search process, so the search efficiency is much faster than that of grid and random search.

In this paper, the Bayesian optimization selected the LCB function as acquisition function (Equation (32)) and the iterations were set to 100.

$$\mathrm{LCB}=-\mu \left(x\right)+k\sigma \left(x\right),$$

(32)

where $\mu \left(x\right)$ is the mean value of a certain point, $\sigma \left(x\right)$ was the variance of a certain point, and k was the adjustment parameter.

In Equation (32), the mean represented the final expected effect of the point; the smaller the mean value, the better the final index of the model. The variance indicated the uncertainty of this point; the larger the variance was, the more the point should be explored. The algorithm steps for Bayesian optimization were shown in Algorithm 1:

Algorithm 1: Outline. Steps of Bayesian optimization algorithm.

Input: Microclimate data (T, $C{O}_2$, $PAR$), mature fruit yield ($W_{MA}$), and the range of variables. Output: The best vector (solution). Step 1: Randomly sample n points in the sample space, and calculate the posterior probability distribution of the first n points by Gaussian process regression to obtain the expected mean and variance of each hyperparameter at each value point. Step 2: Get the next sample point according to the acquisition function (Equation (32)), and calculate the objective function value of the sample point (Equation (32)). Step 3: Determine whether the accuracy requirement or the number of iterations is reached. If the condition is not met, add the sample points to the sample point set, repeat the above steps; and if the conditions are met, stop the iteration.

3. Results and Discussion

The equipment used in the relevant simulations and experiments is an Intel® $Cor{e}^{TM}$ i7-7700HQ, armed with a 2.80 GHz processor and 8.00 GB in RAM; the OS of the computer is $Window{s}^{TM}$ 10 Professional Edition, and the machine is equipped with the $MATLA{B}^{TM}$ R2017a version software.

3.1. Integrated Model Validation

For the output of the integrated model, such as the dry weight of mature fruit and leaves, the comparison among the integrated model, TOMGRO and Vanthoor was shown in Figure 1.

According to Figure 1, the dry weight of mature fruits of the integrated model was within the range of the output of TOMGRO and Vanthoor. The dry weight of leaves in the integrated model changed more gently compared to TOMGRO when the leaf area approached the designated maximum (the maximum leaf area index was consistent, both were set to $2m^2/m^2$), while it was different from the dry weight of leaves in Vanthoor which remained completely constant after a certain time, and was closer to the actual situation. Figure 1 indicated the rationality of the results of the integrated model.

3.2. Parameter Analysis Results

After the integrated model was obtained, the parameters to be optimized were sorted out through sensitivity analysis and then the optimization algorithm was carried out.In this paper, EFAST was selected as the sensitivity analysis algorithm. The table and figure of sensitivity analysis results were shown in

Table 2. Sensitivity analysis results.

Parameters	First Order Sensitivity Index	Total Sensitivity Index
$GENRAT$	0.4247	0.446388
$FPN$	3.98 × 10${}^{-5}$	0.019022
K	0.0549	0.075865
m	0.000366	0.017531
$ϵ$	0.0536	0.074421
$\tau$	0.0199	0.039272
$P_{root}$	0.000436	0.019621
$S_{C{O}_2}$	0.0064	0.027562
E	0.0966	0.116923
$RES{P}_L$	0.0011	0.039419
$RES{P}_F$	0.000514	0.020610
$r_L$	0.0067	0.065230
$r_F$	0.1290	0.148063

and Figure 2:

A parameter with sensitivity index exceeding 0.1 is generally considered as the sensitive parameter [40]. The above figure indicated that among the 13 parameters selected, $GENRAT$, E and $r_F$ had over 0.1 sensitivity index, suggesting these parameters had a significant impact on the growth model. $GENRAT$ represented the maximum appearance rate of new nodes, which had a direct effect on the number of nodes, and the relevant studies [4] indicated that its value was related to tomato varieties meaning that re-optimization on different tomato varieties should be arranged to improve the accuracy of the model. E represented the conversion coefficient of photosynthesis assimilation rate to dry matter accumulation rate, which had a direct impact on $SUPPLY$, but its value was fixed in TOMGRO, so no extra optimization was conducted. $r_F$ was the fruit development rate which had a direct impact on the number of fruits. Its value was related to tomato varieties, meaning that re-optimization on different tomato varieties should be arranged to reduce errors of the model like $GENRAT$. In order to reduce the model error, the above two sensitive parameters ($GENRAT$ and $r_F$) should be optimized. Among the remaining parameters, it was found that the sensitivity value of some parameters which had been expected to be optimized [4] was very low, such as $FPN$, m, $\tau$, $P_{root}$, etc. Although these parameters changed with the greenhouse structure, the type of crops, etc., all these changes had little influence on the final results of the model. Therefore, these parameters can be considered as fixed, and their values are set as recommended values or average values according to some studies [21]. The last remaining value $S_{C{O}_2}$ was a fixed value, and its impact on the final results of the model was very small, verifying the idea of ignoring this parameter in the previous section.

Therefore, according to the above results, for the goal of yield, the parameter processing results were obtained as below:

Finally, combining the results in

Table 3. Parameter processing result.

Types	Parameters	Approach
High sensitivity value, parameters to be optimized	$GENRAT$,$r_F$	need to be optimized
High sensitivity value, fixed parameters	E	fixed
Low sensitivity value, parameters to be optimized	$FPN$,K,m,$ϵ$,$\tau$,$RES{P}_L$,$RES{P}_F$,$r_L$	fixed
Low sensitivity, fixed parameters	$S_{C{O}_2}$,$P_{root}$	ignored

and the parameter that need to be optimized by default, there are three parameters should be optimized, namely $GENRAT$, $r_F$, and $POF$. According to Equations (4), (8) and (18), the following parameters for optimization was obtained, and the default recommended values were attached, as shown in

Table 4. Parameters for optimization.

Parameters	Default Recommended Value
$a_{rF}$	0.0009389
$b_{rF}$	0.000756
$a_g$	0.779
$b_g$	−0.000458
$a_F$	1.2653
$b_F$	0.04295
$c_F$	46.34

. The Appendix Data Collection Table showed the information of the remaining parameters in the model.

3.3. Results of Yield Prediction

The parameters obtained after the sensitivity analysis of the model can be optimized based on greenhouse data, and the results of the yield prediction can be obtained finally.

The microclimate data was selected from Chongming Greenhouse in 2014 (the data selected from September 20 to December 28, an actual growth cycle of Chongming Greenhouse). The range of parameters to be optimized was +/−20% from the default value (After comparative testing, 20% was a reasonable and better range). The optimization results of the Bayesian optimization algorithm were shown in

Table 5. Bayesian optimization algorithm results.

Parameters	Parameter Range	Optimization Result
$a_{rF}$	[0.00075112,0.00112668]	0.00079752
$b_{rF}$	[0.0006048,0.0009072]	0.00086825
$a_g$	[0.6232,0.9348]	0.63451
$b_g$	[−0.0005496,−0.0003664]	−0.00046959
$a_F$	[1.01224,1.51836]	1.1745
$b_F$	[0.03436,0.05154]	0.051295
$c_F$	[37.072,55.608]	37.284

.

The RMSE changed in line with the iterations in Bayesian optimization algorithm and PSO were seen as follows:

It can be seen from the results in the Figure 3 that by using Bayesian optimization to identify the growth model parameters, the RMSE of the yield were 2.5217, 2.5974, 2.2309, 2.5974 and 2.5096 (ran five times continuously to reduce the impact of randomness, and the worst result was selected), The RMSE of PSO were 2.9534, 4.0855, 2.0732, 2.7776 and 3.5322. It can be seen that most of the RMSE of PSO were higher than the worst result of Bayesian optimization (except 2.0732), and the results of PSO were not as stable as the results of Bayesian optimization (a large RMSE: 4.0855 appeared in the results). In addition, the runtime of the Bayesian optimization was only about 2 h, 1/10 the time of the PSO. The reason was that the information of all previous iterations was utilized for the next search in the Bayesian optimization algorithm, making it more efficient than the grid and random search.

The optimized parameters were re-substituted into the model to compare with the actual yield, TOMGRO and Vanthoor, as shown in Figure 4:

In Figure 4, the RMSE between the TOMGRO and the actual yield was 17.2154, the RMSE of Vanthoor was 17.6157, and the RMSE of the model proposed in this paper was 2.5974. It was not hard to find that compared with the output of TOMGRO and Vanthoor (both with calibrated parameters), the output of the model proposed in this paper was closer to the actual yield. This proved that the model had higher accuracy in predicting greenhouse crop yield.

The integrated yield prediction model proposed in this paper also had the ability to “memorize”, meaning that after one parameter optimization in a particular greenhouse, the model can run accurately in the greenhouse for a long time without the second-time optimization. After conducting the optimization algorithm with the data of Chongming Greenhouse in 2014, the study used other growth cycle data (Chongming Greenhouse from 2015 to 2018) as a verification set to validate the scheme. The results were shown in Figure 5.

It can be seen from Figure 5 that after using the optimized parameters (tuned with the data of Chongming Greenhouse in 2014),the model output can approximate the actual yield well when different climate data (2015 to 2018) was input. By comparing with TOMGRO and Vanthoor, the RMSE of the integrated model, TOMGRO and Vanthoor was 9.7508, 15.4559 and 27.0760 in 2015, 6.7849, 13.6584 and 14.5267 in 2016, 7.4368, 16.9145 and 49.4001 in 2017, and 14.5421, 18.7890 and 22.0839 in 2018. Therefore, this proved once again that the integrated model combined with the optimization algorithm had a high accuracy in predicting yield. The reason was that the model combined the advantages of TOMGRO and Vanthoor. Compared with Vanthoor, it had a more detailed simulation of the crop mechanism, so the model output will not deviate too much from the actual yield. Compared with TOMGRO, it had more accurate parameters, so its output can be closer to the actual yield after the optimization algorithm.

4. Conclusions

To solve the problem of poor reusability of the existing greenhouse crop yield prediction models, an integrated yield prediction model for greenhouse tomato was proposed in the study, which included an integrated model, sensitivity analysis and an optimization algorithm. The summary of this paper was as follows:

(1) By analyzing the two representative growth theories of source-sinkand buffer on the basis of the research of TOMGRO, Vanthoor and related growth models, an integrated model for greenhouse tomato was proposed. Through the comparison of the model output, it was found that the output of the integrated model is reasonable and compared with TOMGRO and Vanthoor, it was more general.

(2) Then, the sensitive parameters which have a great influence on the model output and vary with the greenhouse environment should be determined by sensitivity analysis, and they were marked as parameters to be optimized. Other parameters were classified into fixed and ignored.

(3) Bayesian optimization was then used as the optimization algorithm. The RMSE of the yield in the whole growth cycle was taken as the objective function, and the actual greenhouse data was used for parameter optimization. The comparison results of the model output showed that the predicted value of the integrated model after parameter optimization was closer to the actual yield than TOMGRO and Vanthoor, indicating that the integrated yield prediction model has good performance in predicting yield.

(4) One-year greenhouse data was input to conduct parameter optimization for the integrated model. The model output can well track the actual yield for the next four years, compared to TOMGRO and Vanthoor, which proved once again that the model had high accuracy in predicting the yield, and indicated that once the greenhouse parameter optimization was carried out, the model can accurately predict the yield in the greenhouse for a long period of time.