Forecasting of Winter Wheat Yield: A Mathematical Model and Field Experiments

Abstract

An increase in world population requires growth in food production. Wheat is one of the major food crops, covering 21% of global food needs. The food supply issue necessitates reliable mathematical methods for predicting wheat yields. Crop yield information is necessary for agricultural management and strategic planning. Our mathematical model was developed based on a three-year field experiment in a semi-arid climate zone. Wheat yields ranged from 4310 to 6020 kg/ha. The novelty of this model is the inclusion of some stochastic data (weather and technological). The proposed method for wheat yield modeling is based on the theory of random sequence analysis. The model does not impose any restrictions on the number of production parameters and environmental indicators. A significant advantage of the proposed model is the absence of limits on the yield function. Consideration of the stochastic features of wheat production (technological and weather parameters) allows researchers to achieve the best accuracy. The numerical experiment confirmed the high accuracy of the proposed mathematical model for the prediction of wheat yield. The mean relative error (for the third-order polynomial model) varied from 1.79% to 2.75% depending on the preceding crop.

1. Introduction

The world population is increasing and it is predicted to exceed 9 billion by 2050 [1]. Agriculture is required to significantly increase (up to 110%) food production to meet food demands [2,3]. Thus, population growth has a negative impact on food resources [4]. The food and agricultural organization (FAO) states that food security has become an urgent problem for a number of countries [5,6].

Wheat is ranked third among food crops [7]. This crop covers around 21% of the world’s food demand. Around 220 million hectares of arable land worldwide are used for wheat cultivation [8]. Total wheat production exceeds 700 million tons [9,10]. The European Union is ranked first in wheat production [11]. China is a world leader in wheat production. Among European countries, Germany, France, and Ukraine are in the top ten wheat producers [12].

In the last decade, there has been a decrease in wheat production [13,14,15,16]. Climate change and biofuel production are among the primary problems of production growth. Average air temperature is rising [17]. Agricultural production, including wheat production, is sensitive to climate change, which affects crop yields [18,19,20,21]. The climate becomes drought. Thereby, irrigated agriculture is expected to expand. This kind of agriculture currently utilizes around 20% of arable land and produces up to 40% of total food [22]. The production of biofuels has increased. Reducing greenhouse gas emissions, diversifying vehicle fuels, and promoting renewable energy are the main reasons for the above [23]. The European Union has determined main targets to mitigate climate change. They include a 20% increase in renewable energy consumption [24]. Biofuels (based on crop origin) are an alternative to fossil fuels [25,26]. Their use in transport facilities is a priority for reducing carbon dioxide emissions in many countries [27]. In 2017, the European Union used around 7% of the total arable land for industrial crops [28]. Thus, biofuel production reduces available land for food production.

Agriculture must provide enough food for the growing population. Therefore, agricultural management and policymakers should use forecasting methods for crop production [6,29]. Statistics, modeling, time series, learning machine, etc., are used for prediction [30]. Reliable wheat yield forecasting is imperative. These methods help farmers to monitor yield and identify threats (weather conditions, fertilizer management, etc.) [31]. Crop yield information is necessary for strategic planning [32].

Yield is an important indicator for characterizing grain production. Yield forecasting is an important task for any country. The accuracy of forecasting determines the solution of some problems, including organizing reserve funds of food, volumes of grain storage, etc. It affects the formation of an effective foreign trade policy (including the import/export plan and price, optimum management of growing crops). The yield indicator is also the basis for assessing the profitability of agricultural companies. Therefore, estimating yield is an important tool for effective management.

In Ukraine, since 1990, wheat yield has ranged from 1980 to 4160 kg/ha [33]. It is behind the world’s leaders; 7530 kg/ha in Germany [34]. Irrigated winter wheat yield is in the range of 3550 to 5290 kg/ha [35]. In the world, irrigated winter wheat yield is up to 7990 kg/ha in arid and semi-arid zones [36]. Thus, there are reserves for increasing yields. It is necessary to optimize the use of available energy and material resources to realize these reserves. It is important for the arid and semi-arid zones of both Ukraine and other countries.

Significant weather and economic instability dictate the importance of yield forecasting. Crop yield forecasting is difficult because crop formation is associated with factors such as agricultural practice, weather conditions, the characteristics of biological systems, etc.

Currently, various methodical approaches to yield forecasting have been developed and applied in practice:

• The analysis of trends and cyclicality in yield dynamics [37,38];
• The identification of the year-analog [39,40,41];
• The building of regression models based on a set of statistics obtained on the basis of remote and meteorological observations [38,42];
• Modeling [38];
• The analysis of synoptic processes [43].

The approaches of the first, second, and fifth groups are not accurate enough. Groups 3 and 4 are the most widely used approaches. In most cases, meteorological data are used to build regressions or to model plant growth. This type of forecast does not take into account the actual state of the soil, the use of fertilizers, and other chemicals. Dynamic models are most widely used in practice [5]. However, they do not take into account the entire history of changes in yields and the conditions of grain production, which significantly limit the accuracy of existing models.

The main feature of yield is its stochastic changes. In this regard, the theory of random functions and random sequences must be used to predict crop yields. Methods and algorithms of the above theory take into account various random factors (precipitation, air temperature, soil temperature, etc.) [44], as well as the values of deterministic parameters (soil structure, crop variety, tillage practices, dosage and composition of fertilizers, etc.) [45,46,47,48,49]. Crop yields have been studied using simulation models [50,51]. These models most accurately reveal the impact of agronomic factors on crop yields [52,53,54].

However, practitioners and authorities require mathematical models for different crops, climate zones, tillage practices, etc. The purpose of this article is to develop a mathematical model for predicting the winter wheat yield in a semi-arid zone. The novelty of this study is the development of a mathematical model based on stochastic data, such as precipitation, plant density, fertilizer, micro fertilizers, the effective temperature sum of the autumn vegetation, the amount of water used for irrigation, and preceding crops. This model has been developed for a semi-arid zone. Its modeling algorithm was built based on random non-stationary sequences of input variables. The main requirement for the method developed is the absence of any significant restrictions on the random process of grain crop yields. The maximum consideration of stochastic characteristics will allow us to achieve the best accuracy of the modeling problem.

This study is based on previous publications [55,56,57].

2. Materials and Methods

This study focuses on the forecasting of winter wheat yield and proposes the use of a methodology combining statistical analysis and performing field experiments. This methodology comprises the following steps: the collection of field experiment data; the modeling of yield as a function of selected parameters. Field experiments were performed in the Mykolaiv region (Ukraine).

2.1. Field Experiment

The experiments were carried out in 2019 and 2020 in the Mykolaiv province, Ukraine (46°58′06″ N; 31°42′39″ E). The area of the experimental field is equal to 10 ha. Our experiment had a randomized design with three replications. Rapeseed and corn were preceding crops for winter wheat. The soil had the following properties: pH—from 6.8 to 7.2; organic carbon—from 2.9 to 3.2 g·kg⁻¹; phosphorus—from 31 to 38 mg·kg⁻¹; potassium—from 332 to 525 mg·kg⁻¹; bulk density—1380 kg·m⁻³. Winter wheat was grown under practice that is conventional for southern Ukraine (

Table 1. Production of winter wheat.

Farming Operation	Description
Tillage	Skimming (6–8 cm)
	Cultivation (8–10 cm)
	Harrowing
	Pre-seeding cultivation (3–4 cm)
Sowing	25 September–10 October; 4.0, 4.5 and 5 million seeds per hectare
Fertilization	N—from 60 to 120 kg/ha; P—15 kg/ha; K—15 kg/ha
Micronutrient fertilizers	I scheme: Fitohelp (0.5 L/ha) + Liposam (0.2 L/ha) II scheme: Quantum-grain (1.0 L/ha) + Liposam (0.2 L/ha)
Irrigation	I scheme: pre-sowing irrigation—600 m3/ha, other irrigation—1000 m3/ha II scheme: pre-sowing irrigation—700 m3/ha, other irrigation—900 m3/ha
Weed control	Chemicals: 0.025 kg/ha (glyphosate); 4 kg/ha (PIK)
Harvesting	June–July

).

The experiments were carried out over two years on irrigated lands. A field experiment was established by the method of randomized split plots. All studies, observations, and samplings were performed in quadruplicate. We used a sequential arrangement of plots in one tier. They were located in relation to organizational and technical factors: the convenience of tillage, fertilization, sowing, harvesting, etc. The total number of plots was 32. We investigated six factors in the field experiments. Factor A was the preceding crops (rapeseed and corn); factor B was the plant density, million plants per hectare: 4.0, 4.5, and 5.0; factor C was the fertilizer type and dosage (N:P:K); factor D was the micronutrient fertilizer type and dosage; factor E was the irrigation scheme; factor F was the total effective temperature of the autumn growing season. The sum of average air temperatures in autumn is the sum of average daily temperatures above +5 °C. This indicator characterizes the amount of heat necessary for the plant development process. Wheat during the autumn vegetation should gain a sum of effective temperatures from 300 to 350 °C. During the three-year experiments, differences in precipitation were observed. This study considered the influence of annual precipitation on wheat yield.

A preceding crop leaves nutrients in the soil. Thus, preceding crops influence yield. Two preceding crops (maize and rapeseed) were the limitation of this study. We chose maize because it is a widespread crop, and its production is around 50% of gross national grain production. Rapeseed is one of the best preceding crops.

2.2. Measurement of Yield

The method of mechanized harvesting was used to determine winter wheat yield. Wheat grain was harvested by a Sampo 500 combine harvester. Harvesting was carried out from a selective typical plot of the field. The yield was calculated by the equation:

$$Yd=MG\cdot P{A}^{-1}, \mathrm{kg}/\mathrm{ha},$$

(1)

where MG is the mass of wheat grain from a plot, kg; PA is the area of a plot, ha.

2.3. Methodology for the Synthesis of Winter Wheat Yield Models

The formation and use of models of changes in winter wheat yield indicators involves the implementation of the following stages:

• Stage 1. The collection of statistical data on grain yields and cultivation conditions;
• Stage 2. The estimation of moment functions $M\left[X_{}^{\lambda }\left(\nu \right)X_{}^s\left(i\right)\right]$ based on obtained statistical data;
• Stage 3. The determination of the optimal order of stochastic connections of the random vector $\left\{X\right\}$;
• Stage 4. The calculation of the characteristics of the canonical distribution of the random vector $\left\{X\right\}$;
• Stage 5. The calculation of the parameters of the mathematical model;
• Stage 6. The calculation of productivity indicators based on the predictive model under different initial conditions of production;
• Stage 7. The assessment of yield modeling accuracy.

2.4. Verification of a Developed Mathematical Model

The developed mathematical model was verified using such indicators as the mean relative error, the standard deviation of error, and the coefficient of error variation. A mean relative error is

$$\delta =\frac{{\displaystyle \sum _{i=1}^n\left|\frac{Y{m}_i-Y{e}_i}{Y{e}_i}\right|\cdot 100\%}}{n},$$

(2)

where Ym_i is the ith yield calculated by the mathematical model, kg/ha; Ye_i is the ith experimental yield, kg/ha; n is the number of experimental yields.

To find the standard deviation, we used the following formula

$$\sigma =\sqrt{\frac{{{\displaystyle \sum _{i=1}^n\left({\delta }_i-\delta \right)}}^2}{n}},$$

(3)

where δ_i is the ith error, %.

Finally, the coefficient of error variation is equal to

$$CV=\frac{\sigma }{\delta }\cdot 100\%.$$

(4)

3. Results and Discussion

3.1. Field Experiment

Three-year field experiments on winter wheat growing were carried out on a farm of Mykolaiv National Agrarian University (the Mykolaiv region). The results are presented in

Table A1. Yields, kg/ha (preceding crop—maize).

#	Plant Density, mln/ha	Fertilizer, N:P:K, kg/ha	Microfertilizers		Sum of Autumn Effective Temperatures	Irrigation
			Fitohelp (0.5 L/ha) + Liposam (0.2 L/ha)	Quantum-Grain (1.0 L/ha) + Liposam м (0.2 L/ha)		600 m3/ha + 1000 m3/ha	700 m3/ha + 900 m3/ha
1.	4.0	60:15:15	I		369	4310	4390
2.		60:15:15		II	369	4430	4480
3.		75:15:15	I		369	4410	4430
4.		75:15:15		II	369	4490	4540
5.		90:15:15	I		369	4480	4550
6.		90:15:15		II	369	4550	4620
7.		105:15:15	I		369	4530	4600
8.		105:15:15		II	369	4630	4710
9.		120:15:15	I		369	4790	4890
10.		120:15:15		II	369	4890	5020
11.	4.5	60:15:15	I		369	4410	4510
12.		60:15:15		II	369	4500	4600
13.		75:15:15	I		369	4470	4570
14.		75:15:15		II	369	4510	4650
15.		90:15:15	I		369	4560	4630
16.		90:15:15		II	369	4610	4740
17.		105:15:15	I		369	4680	4810
18.		105:15:15		II	369	4810	4850
19.		120:15:15	I		369	4950	4990
20.		120:15:15		II	369	5000	5100
21.	5.0	60:15:15	I		369	4500	4570
22.		60:15:15		II	369	4570	4680
23.		75:15:15	I		369	4550	4610
24.		75:15:15		II	369	4610	4680
25.		90:15:15	I		369	4630	4720
26.		90:15:15		II	369	4700	4790
27.		105:15:15	I		369	4740	4910
28.		105:15:15		II	369	4810	4940
29.		120:15:15	I		369	4980	5030
30.		120:15:15		II	369	5100	5150
31.	4.0	60:15:15	I		355	4520	4550
32.		60:15:15		II	355	4670	4780
33.		75:15:15	I		355	4710	4680
34.		75:15:15		II	355	4750	4810
35.		90:15:15	I		355	4810	4860
36.		90:15:15		II	355	4900	4980
37.		105:15:15	I		355	4930	4910
38.		105:15:15		II	355	4950	5020
39.		120:15:15	I		355	4970	5050
40.		120:15:15		II	355	5000	5100
41.	4.5	60:15:15	I		355	4600	4670
42.		60:15:15		II	355	4690	4730
43.		75:15:15	I		355	4720	4810
44.		75:15:15		II	355	4810	4910
45.		90:15:15	I		355	4950	5010
46.		90:15:15		II	355	5050	5120
47.		105:15:15	I		355	5090	5180
48.		105:15:15		II	355	5120	5210
49.		120:15:15	I		355	5180	5240
50.		120:15:15		II	355	5290	5370
51.	5.0	60:15:15	I		355	4720	4790
52.		60:15:15		II	355	4850	4910
53.		75:15:15	I		355	4890	4930
54.		75:15:15		II	355	4910	4940
55.		90:15:15	I		355	4930	4990
56.		90:15:15		II	355	4990	5090
57.		105:15:15	I		355	5110	5210
58.		105:15:15		II	355	5160	5240
59.		120:15:15	I		355	5220	5280
60.		120:15:15		II	355	5360	5400
61.	4.0	60:15:15	I		345	4700	4730
62.		60:15:15		II	345	4830	4850
63.		75:15:15	I		345	4910	4940
64.		75:15:15		II	345	4930	4970
65.		90:15:15	I		345	4990	5020
66.		90:15:15		II	345	5150	5210
67.		105:15:15	I		345	5190	5240
68.		105:15:15		II	345	5210	5290
69.		120:15:15	I		345	5250	5320
70.		120:15:15		II	345	5380	5410
71.	4.5	60:15:15	I		345	4800	4840
72.		60:15:15		II	345	4950	4990
73.		75:15:15	I		345	5020	5100
74.		75:15:15		II	345	5080	5120
75.		90:15:15	I		345	5120	5150
76.		90:15:15		II	345	5280	5300
77.		105:15:15	I		345	5310	5390
78.		105:15:15		II	345	5330	5410
79.		120:15:15	I		345	5350	5470
80.		120:15:15		II	345	5410	5590
81.	5.0	60:15:15	I		345	4970	4990
82.		60:15:15		II	345	5050	5060
83.		75:15:15	I		345	5080	5180
84.		75:15:15		II	345	5210	5290
85.		90:15:15	I		345	5290	5330
86.		90:15:15		II	345	5450	5510
87.		105:15:15	I		345	5510	5570
88.		105:15:15		II	345	5570	5630
89.		120:15:15	I		345	5610	5680
90.		120:15:15		II	345	5750	5860

(preceding crop—maize) and

Table A2. Yields, kg/ha (preceding crop—rapeseed).

#	Plant Density, mln/ha	Fertilizer, N:P:K, kg/ha	Microfertilizers		Sum of Autumn Effective Temperatures	Irrigation
			Fitohelp (0.5 L/ha) + Liposam (0.2 L/ha)	Quantum-Grain (1.0 L/ha) + Liposam м (0.2 L/ha)		600 m3/ha + 1000 m3/ha	700 m3/ha + 900 m3/ha
91.	4.0	60:15:15	I		369	4440	4480
92.		60:15:15		II	369	4530	4520
93.		75:15:15	I		369	4460	4500
94.		75:15:15		II	369	4520	4530
95.		90:15:15	I		369	4480	4510
96.		90:15:15		II	369	4640	4690
97.		105:15:15	I		369	4550	4620
98.		105:15:15		II	369	4720	4970
99.		120:15:15	I		369	4970	5080
100.		120:15:15		II	369	4990	5180
101.	4.5	60:15:15	I		369	4690	4720
102.		60:15:15		II	369	4750	4750
103.		75:15:15	I		369	4740	4780
104.		75:15:15		II	369	4820	4940
105.		90:15:15	I		369	4980	4990
106.		90:15:15		II	369	5120	5050
107.		105:15:15	I		369	5010	5030
108.		105:15:15		II	369	5040	5080
109.		120:15:15	I		369	4710	4800
110.		120:15:15		II	369	5090	5140
111.	5.0	60:15:15	I		369	4790	4850
112.		60:15:15		II	369	4840	4970
113.		75:15:15	I		369	4820	4870
114.		75:15:15		II	369	4860	5020
115.		90:15:15	I		369	5090	5090
116.		90:15:15		II	369	5210	5150
117.		105:15:15	I		369	5110	5140
118.		105:15:15		II	369	5230	5260
119.		120:15:15	I		369	4800	4990
120.		120:15:15		II	369	4990	5280
121.	4.0	60:15:15	I		355	4900	4950
122.		60:15:15		II	355	4990	5110
123.		75:15:15	I		355	4920	4970
124.		75:15:15		II	355	5060	5140
125.		90:15:15	I		355	5070	5100
126.		90:15:15		II	355	5190	5250
127.		105:15:15	I		355	5090	5140
128.		105:15:15		II	355	5220	5280
129.		120:15:15	I		355	5180	5200
130.		120:15:15		II	355	5290	5370
131.	4.5	60:15:15	I		355	5100	5150
132.		60:15:15		II	355	5150	5200
133.		75:15:15	I		355	5120	5170
134.		75:15:15		II	355	5190	5240
135.		90:15:15	I		355	5290	5390
136.		90:15:15		II	355	5410	5520
137.		105:15:15	I		355	5310	5420
138.		105:15:15		II	355	5570	5640
139.		120:15:15	I		355	5530	5610
140.		120:15:15		II	355	5640	5780
141.	5.0	60:15:15	I		355	5100	5090
142.		60:15:15		II	355	5120	5230
143.		75:15:15	I		355	5140	5160
144.		75:15:15		II	355	5200	5300
145.		90:15:15	I		355	5260	5300
146.		90:15:15		II	355	5410	5510
147.		105:15:15	I		355	5280	5330
148.		105:15:15		II	355	5440	5550
149.		120:15:15	I		355	5440	5530
150.		120:15:15		II	355	5610	5720
151.	4.0	60:15:15	I		345	5030	5050
152.		60:15:15		II	345	5080	5100
153.		75:15:15	I		345	5050	5070
154.		75:15:15		II	345	5110	5140
155.		90:15:15	I		345	5180	5250
156.		90:15:15		II	345	5240	5380
157.		105:15:15	I		345	5200	5270
158.		105:15:15		II	345	5290	5420
159.		120:15:15	I		345	5310	5400
160.		120:15:15		II	345	5400	5580
161.	4.5	60:15:15	I		345	5160	5180
162.		60:15:15		II	345	5180	5220
163.		75:15:15	I		345	5200	5230
164.		75:15:15		II	345	5240	5260
165.		90:15:15	I		345	5280	5370
166.		90:15:15		II	345	5320	5450
167.		105:15:15	I		345	5300	5410
168.		105:15:15		II	345	5340	5440
169.		120:15:15	I		345	5380	5490
170.		120:15:15		II	345	5410	5580
171.	5.0	60:15:15	I		345	4950	5000
172.		60:15:15		II	345	5110	5120
173.		75:15:15	I		345	5010	5030
174.		75:15:15		II	345	5130	5210
175.		90:15:15	I		345	5360	5450
176.		90:15:15		II	345	5460	5580
177.		105:15:15	I		345	5410	5510
178.		105:15:15		II	345	5490	5610
179.		120:15:15	I		345	5680	5730
180.		120:15:15		II	345	5910	6020

(preceding crop—rapeseed). Wheat yields were in the range from 4310 kg/ha to 6020 kg/ha. We can see that rapeseed was a better preceding crop than maize. This predecessor provided higher yields (from 4440 kg/ha to 6020 kg/ha). It was 3% higher compared to maize. Experimental data were used to develop our mathematical model.

We analyzed the impact of changing each factor (variable) on the yield. The results of the analysis are presented in

Table 2. The weight of each variable.

Variable	Unit	Range		Average Minimum	Average Maximum	Relative Increase, %
		Minimum	Maximum
preceding crop	-	maize	rapeseed	4953.89	5151.67	3.99
microfertilizer	-	I	II	5000.17	5105.39	2.10
nitrogen	kg/ha	60	120	4830.42	5956.88	23.32
plant density	mln/ha	4	5	4921.33	5155.25	4.75
the effective temperature sum of the autumn vegetation	°C	345	369	5265.67	4774.42	−9.33
irrigation scheme	m3/ha	600 + 1000	700 + 900	5017.50	5088.06	1.41

. We assumed that all variables, except for the one being studied, are constants. As can be seen, nitrogen fertilizers had the strongest effect on yield. The irrigation scheme had the least influence. The sum of average air temperatures in autumn ranked second, with a value of −9.33%. Plant density and a preceding crop had approximately the same effect on the yield.

We analyzed the variables that most strongly affect yield: nitrogen, the effective temperature sum of the autumn vegetation, and plant density (Figure 1, Figure 2 and Figure 3). Linear functions were built to determine their slopes. We can see that rapeseed was a better preceding crop for winter wheat. Rapeseed allows farmers to obtain higher yields. Slopes of linear functions varied from −22.17 to 8.64 (

Table 3. The slope of a linear function.

Variable	Preceding Crop
	Maize	Rapeseed
Nitrogen	8.68	6.24
The effective temperature sum of the autumn vegetation	−22.17	−19.27
Plant density	0.23	0.24

). The only variable that had a negative slope was the effective temperature sum of the autumn vegetation. An increase in the effective temperature sum of the autumn vegetation decreased winter wheat yield. However, fertilizer management could offset this negative impact.

3.2. Canonical Decomposition of a Random Sequence of Yield Index and Characteristics of Production Conditions

To form realizable algorithms for modeling random sequences, certain restrictions are imposed on the properties of the sequence under study. For example, it is assumed that:

(a)

the sequence under study is normal (has a normal probability distribution law) or a stationary sequence generated by a normal one in nonlinear systems;

(b)

the sequence is non-stationary, but with stationary increments;

(c)

the sequence is Markovian, etc.

For these classes of random sequences, there are quite efficient modeling algorithms. Therefore, to obtain a random sequence with a given correlation matrix (without taking into account distribution densities), the method of linear transformations is successfully used [58,59]. One of the varieties of the method of linear transformations—the canonical decomposition of V.S. Pugachev [60]—allows us to form the values of a sequence of random variables that are dependent within the framework of linear relationships, taking into account their one-dimensional distribution densities. The Fourier series is widely used to model a stationary random sequence [61]; the apparatus for modeling stationary normal sequences is well developed [62] (there are two operators for generating values and several approaches have been developed [63,64] for determining their parameters); the simplest solution is the problem of modeling Markov sequences [65], which is reduced to the method of conditional distributions for the simplest case—the use of only a two-dimensional distribution density. However, the use of simplifying assumptions about the properties of a random sequence in the formation of a modeling algorithm naturally reduces the accuracy of the representation of a random sequence. The most universal from the point of view of the restrictions (linearity, Markov property, stationarity, monotonicity, scalarness, etc.), which are superimposed on the properties of sequences of random variables, is the method based on a non-linear canonical decomposition [66].

The subject of study is a random sequence $\left\{X\left(i\right)\right\},$$i=\overline{1,I},$ where $X\left(i\right), i=\overline{1,I-1}$—random values that determine the conditions of production of grain crops (temperature, amount of precipitation, number of sunny days, volume of mineral and organic fertilizers, etc.), $X_{}\left(I\right)$—indicator of grain crop yield.

The nonlinear canonical decomposition of the investigated vector $\left\{X\left(i\right)\right\},$ can be written as [67]:

$$X\left(i\right)=M\left[X\left(i\right)\right]+{\displaystyle \sum _{\nu =1}^i{\displaystyle \sum _{\lambda =1}^N{W}_{\nu }^{(\lambda )}}}{\vartheta }_{\nu }^{(\lambda ;1)}\left(i\right), i=\overline{1,I},$$

(5)

The random coefficients $W_{\nu }^{(\lambda )}$ and non-random coordinate functions ${\vartheta }_{\nu }^{(\lambda ;h)}\left(i\right)$ of the mathematical yield model (1) are determined by the recurrence relations:

$$\begin{array}{c}{W}_{\nu }^{(\lambda )}=X^{\lambda }\left(\nu \right)-M\left[X^{\lambda }\left(\nu \right)\right]-{\displaystyle \sum _{\mu =1}^{\nu -1}{\displaystyle \sum _{j=1}^N{W}_{\mu }^{(j)}}}{\vartheta }_{\mu }^{(j;\lambda )}\left(\nu \right)-\\ -{\displaystyle \sum _{j=1}^{\lambda -1}{W}_{\nu }^{(j)}}{\vartheta }_{\nu }^{(j;\lambda )}\left(\nu \right), \lambda =\overline{1,N},\nu =\overline{1,I};\end{array}$$

(6)

$$\begin{gathered} {\vartheta }_{\nu }^{(\lambda ;h)}\left(i\right)=\frac{M\left[W_{\nu }^{(\lambda )}\left(X^h\left(i\right)-M\left[X^h\left(i\right)\right]\right)\right]}{M\left[{\left\{W_{\nu }^{(\lambda )}\right\}}^2\right]}=\frac{1}{D_{\lambda }\left(\nu \right)}\{M\left[X^{\lambda }\left(\nu \right)X^h\left(i\right)\right]- \\ -M\left[X^{\lambda }\left(\nu \right)\right]M\left[X^h\left(i\right)\right]-{\displaystyle \sum _{\mu =1}^{\nu -1}{\displaystyle \sum _{j=1}^N{D}_j}}\left(\mu \right){\vartheta }_{\mu }^{(j;\lambda )}\left(\nu \right){\vartheta }_{\mu }^{(j;h)}\left(i\right)- \\ -{\displaystyle \sum _{j=1}^{\lambda -1}{D}_j}\left(\nu \right){\vartheta }_{\nu }^{(j;\lambda )}\left(\nu \right){\vartheta }_{\nu }^{(j;h)}\left(i\right)\}, \lambda =\overline{1,h}, \nu =\overline{1,i}, h=\overline{1,N}, i=\overline{1,I}. \end{gathered}$$

(7)

$$\begin{array}{c}{D}_{\lambda }\left(\nu \right)=M\left[{\left\{W_{\nu }^{(\lambda )}\right\}}^2\right]=M\left[X^{2\lambda }\left(\nu \right)\right]-M^2\left[X^{\lambda }\left(\nu \right)\right]-\\ -{\displaystyle \sum _{\mu =1}^{\nu -1}{\displaystyle \sum _{j=1}^N{D}_j\left(\mu \right)}}{\left\{{\vartheta }_{\mu }^{(j;\lambda )}\left(\nu \right)\right\}}^2-{\displaystyle \sum _{j=1}^{\lambda -1}{D}_j\left(\nu \right)}{\left\{{\vartheta }_{\nu }^{(j;\lambda )}\left(\nu \right)\right\}}^2, \lambda =\overline{1,N},\nu =\overline{1,I};\end{array}$$

(8)

where $M\left[ \right]$ is the mathematical expectation; $D_{\lambda }\left(\nu \right)$ are the variances of the random coefficient $W_{\nu }^{(\lambda )}, \lambda =\overline{1,N},\nu =\overline{1,I}$.

Coordinate functions ${\vartheta }_{\nu }^{(\lambda ;h)}\left(i\right), \nu =\overline{1,i}; \lambda ,h=\overline{1,N}; i=\overline{1,I}$ are characterized by relations

$${\vartheta }_{\nu }^{(\lambda ;h)}\left(i\right)=\left\{\begin{array}{l}1, \mathrm{for} \left(h=\lambda \right)\wedge \left(\nu =i\right);\\ 0, \mathrm{if} \left(\mathrm{i}<\nu \right)\vee \left(\left(h<\lambda \right)\wedge \left(\nu =i\right)\right).\end{array}\right.$$

The nonlinear model (1) of the random vector $\left\{X\right\}=X\left(i\right), \mathrm{i}=\overline{1,I}$ contains $N$ arrays $\left\{W^{(\lambda )}\right\}, \lambda =\overline{1,N}$ of uncorrelated centered random coefficients $W_i^{(\lambda )}, \lambda =\overline{1,N},i=\overline{1,I}.$ Each of these coefficients contains information about the corresponding value $X^{\lambda }(i), \lambda =\overline{1,N}, i=\overline{1,I},$ and the coordinate functions ${\vartheta }_{\nu }^{(\lambda ;h)}\left(i\right), \lambda ,h=\overline{1,N}, \nu ,i=\overline{1,I}$ describe the probabilistic relations of the order $\lambda +h$ between the point $\nu$ and $i$ $\left(\nu ,i=\overline{1,I}\right).$ Expression (5) provides an optimal description of the studied sequence $\left\{X\right\}$(where $X\left(i\right), \mathrm{i}=\overline{1,I-1}$ are the technological parameters; $X\left(I\right)$ is the yield) according to the criterion of the minimum mean square of the modeling error. Expression (1) is also true if some stochastic relations of the random vector $\left\{X\right\}=X\left(i\right), \mathrm{i}=\overline{1,I}$ are missing. In this case, the corresponding coordinate functions take the value zero and these relations are automatically excluded from the canonical decomposition.

The legitimacy of the approach used to form representation (5) is confirmed by the proposition that it is possible to construct a canonical decomposition of the sequence $\left\{f_1\left({\overline{Y}}_1\right),\dots ,f_n\left({\overline{Y}}_n\right)\right\}$, where ${\overline{Y}}_{\nu }, \nu =\overline{1,n}$ is a vector random variable and $f_{\nu }\left(.\right), \nu =\overline{1,n}$ is a nonlinear function.

3.3. Predictive Model of Changes in Yield Indicators Depending on the Initial Conditions of Production

The sequential fixation of known values $x^{\nu }\left(j\right)$ in the canonical decomposition (5) (the values of random coefficients become known $W_j^{(\nu )}$) using the mathematical expectation operation obtains an extrapolation algorithm [67]:

$$m_x^{(\mu ,l)}\left(h,i\right)=\left\{\begin{array}{c}M\left[X^h\left(i\right)\right] \mathrm{if} \mu =0;\\ m_x^{(\mu ,l-1)}\left(h,i\right)+\left(x^l\left(\mu \right)-m_x^{(\mu ,l-1)}\left(l,\mu \right)\right){\vartheta }_{\mu }^{(l;h)}\left(i\right) \mathrm{if} l\ne 1,\\ m_x^{(\mu -1,N)}\left(h,i\right)+\left(x^l\left(\mu \right)-m_x^{(\mu -1,N)}\left(l,\mu \right)\right){\vartheta }_{\mu }^{(1;h)}\left(i\right) \mathrm{if} l=1.\end{array}\right.$$

(9)

The expression $m_x^{(\mu ,l)}\left(h,i\right)=M\left[X^h\left(i\right)/x^{\nu }\left(j\right),\right.$$j=\overline{1,\mu -1}, \nu =\overline{1,N}; \left.x^{\nu }\left(\mu \right), \nu =\overline{1,l}\right]$ (conditional expectation) for $h=1, l=N, \mu =k$ is an optimal estimate $m_x^{(k,N)}\left(1,i\right)$ of the future value $x\left(I\right)$ of the yield indicator, provided that the values $x^{\nu }\left(j\right), \nu =\overline{1,N}, j=\overline{1,I-1}$ are used to calculate this estimate; i.e., $I-1$ indicators that characterize the conditions of production of grain crops are known.

The expression for estimation $m_x^{(k,N)}\left(1,i\right)$ can be written in the following explicit form [68]:

$$m_x^{(k,N)}\left(1,i\right)=M\left[X\left(i\right)\right]+{\displaystyle \sum _{j=1}^k{\displaystyle \sum _{\nu =1}^N\left(x^{\nu }\left(j\right)-M\left[X^{\nu }\left(j\right)\right]\right)S_{((j-1)N+\nu )}^{(kN)}\left(\left(i-1\right)N+1\right)}},$$

(10)

$$\mathrm{where} S_{\lambda }^{(\alpha )}\left(\xi \right)=\left\{\begin{array}{l}{S}_{\lambda }^{(\alpha -1)}\left(\xi \right)-S_{\lambda }^{(\alpha -1)}\left(\alpha \right){\gamma }_k\left(i\right), \mathrm{if} \lambda \le \alpha -1;\\ {\gamma }_{\alpha }\left(\xi \right), \mathrm{for} \lambda =\alpha ;\end{array}\right.$$

(11)

$${\gamma }_{\alpha }\left(\xi \right)=\left\{\begin{array}{l}{\beta }_{[\alpha /N]+1}^{({\mathrm{mod}}_N(\alpha );1)}\left(\left[\alpha /N\right]+1\right), \mathrm{for} \xi \le kN;\\ {\beta }_{[\alpha /N]+1}^{({\mathrm{mod}}_N(\alpha );1)}\left(i\right), \mathrm{if} \xi =\left(i-1\right)N+1.\end{array}\right.$$

(12)

where ${\mathrm{mod}}_N( )$ is the division modulo $N$.

3.4. Synthesis of Models of Changes in Winter Wheat Yield

Based on statistical yield data from the period 2019–2021, as a result of conducting experiments at the innovative training ground of the Ukrainian National Academy of Sciences, it was determined that the main factors affecting winter wheat yield are as follows [69,70,71,72]:

• average annual precipitation, (m);
• sowing rate, (million seeds/ha);
• mineral nutrition (N:P:K—nitrogen-phosphorus-potassium; kg/ha);
• microfertilizers (L/ha or kg/ha);
• sum of effective temperatures of autumn vegetation (°C);
• volume of water used for irrigation (m³/ha).

That is, the random vector takes the form $\left\{X_{}\left(i\right)\right\},$$i=\overline{1,7}$: $X_{}\left(1\right)$— amount of average annual precipitation; $X_{}\left(2\right)$—sowing rate; $X_{}\left(3\right)$—amount of mineral nutrition; $X_{}\left(4\right)$—amount of microfertilizer; $X_{}\left(5\right)$—the sum of the effective temperatures of autumn vegetation; $X_{}\left(6\right)$—volume of water used for irrigation; $X_{}\left(7\right)$—winter wheat yield.

It was determined that the most significant stochastic relations affecting the yield index of winter wheat are relations of the third order: $X^{\nu }\left(i\right)X^{\mu }\left(j\right)\ne 0, i,j=\overline{1,7}, \nu ,\mu \le 3$; $X^{\nu }\left(i\right)X^{\mu }\left(j\right)=0, i,j=\overline{1,7}, \nu ,\mu >3$ (the analysis of the ripening process was not performed, therefore, the aftereffect interval in this case degenerates into a point: the moment of harvesting).

Table 4. Coordinate functions ${\beta }_{1\nu }^{(1)}\left(i\right); \nu ,\mathrm{i}=\overline{1.7}$ (preceding crop—maize).

0	1	2	3	4	5	6	7
1	0	1	2	3	4	5	6
2	1	1.00	0.18	−0.28	0.47	−0.21	0.59
3	2	0	1.00	0.75	0.33	0.45	−0.1
4	3	0	0	1.00	0.34	0.77	0.27
5	4	0	0	0	1.00	0.55	0.21
6	5	0	0	0	0	1.00	0.23
7	6	0	0	0	0	0	1.00

,

Table 5. Coordinate functions ${\beta }_{1\nu }^{(2)}\left(i\right); \nu ,\mathrm{i}=\overline{1,7}$ (preceding crop—maize).

0	1	2	3	4	5	6	7
1	1.00	−0.16	−0.12	0.34	−0.11	0.55	0.12
2	0	1.00	0.7	−0.22	0.33	−0.07	0.21
3	0	0	1.00	0.34	0.83	0.23	0.31
4	0	0	0	1.00	−0.79	0.11	−0.10
5	0	0	0	0	1.00	0.12	−0.15
6	0	0	0	0	0	1.00	0.21
7	0	0	0	0	0	0	1.00

,

Table 6. Coordinate functions ${\beta }_{1\nu }^{(3)}\left(i\right); \nu ,\mathrm{i}=\overline{1,7}$ (preceding crop—maize).

0	1	2	3	4	5	6	7
1	1.00	−0.11	0.11	0.27	−0.09	0.05	0.12
2	0	1.00	0.34	0.11	0.27	−0.02	−0.01
3	0	0	1.00	0.22	0.33	0.09	0.05
4	0	0	0	1.00	−0.22	0.09	0.05
5	0	0	0	0	1.00	0.11	−0.12
6	0	0	0	0	0	1.00	0.15
7	0	0	0	0	0	0	1.00

,

Table 7. Coordinate functions ${\beta }_{1\nu }^{(1)}\left(i\right); \nu ,\mathrm{i}=\overline{1,7}$ (preceding crop—rapeseed).

0	1	2	3	4	5	6	7
1	1.00	−0.33	0.55	0.43	−0.41	0.75	0.23
2	0	1.00	0.90	−0.41	0.97	−0.28	0.31
3	0	0	1.00	0.39	0.8	0.44	0.22
4	0	0	0	1.00	0.33	0.28	0.15
5	0	0	0	0	1.00	0.22	0.17
6	0	0	0	0	0	1.00	0.31
7	0	0	0	0	0	0	1.00

,

Table 8. Coordinate functions ${\beta }_{1\nu }^{(2)}\left(i\right); \nu ,\mathrm{i}=\overline{1,7}$ (preceding crop—rapeseed).

0	1	2	3	4	5	6	7
1	1.00	0.18	−0.28	0.47	−0.21	0.59	0.22
2	0	1.00	0.75	0.33	0.45	−0.1	0.07
3	0	0	1.00	0.34	0.77	0.27	0.17
4	0	0	0	1.00	0.55	0.21	0.15
5	0	0	0	0	1.00	0.23	−0.12
6	0	0	0	0	0	1.00	0.37
7	0	0	0	0	0	0	1.00

and

Table 9. Coordinate functions ${\beta }_{1\nu }^{(3)}\left(i\right); \nu ,\mathrm{i}=\overline{1,7}$ (preceding crop—rapeseed).

0	1	2	3	4	5	6	7
1	1.00	−0.18	0.21	0.33	−0.11	0.11	0.07
2	0	1.00	0.44	0.21	0.37	−0.02	−0.11
3	0	0	1.00	0.24	0.39	0.19	0.21
4	0	0	0	1.00	0.44	0.19	0.07
5	0	0	0	0	1.00	0.17	0.09
6	0	0	0	0	0	1.00	0.21
7	0	0	0	0	0	0	1.00

show the values of the coordinate functions of the canonical decomposition.

Using the values of the coordinate functions, mathematical models of changing winter wheat yield (t/ha) are formed:

The preceding crop is maize

$$\begin{array}{cc}{X}_{}\left(7\right)=m_k^{\left(6.3\right)}\left(1.7\right)\hfill & =4.902-0.0152\left(X\left(1\right)-M\left[X\left(1\right)\right]\right)-0.0224\left(X\left(2\right)-M\left[X\left(2\right)\right]\right)-\hfill \\ & -0.0529\left(X\left(3\right)-M\left[X\left(3\right)\right]\right)+0.0131\left(X\left(4\right)-M\left[X\left(4\right)\right]\right)-\hfill \\ & -0.006\left(X\left(5\right)-M\left[X\left(5\right)\right]\right)-0.0608\left(X\left(6\right)-M\left[X\left(6\right)\right]\right)+\hfill \\ & +0.0954\left(X^2\left(1\right)-M\left[X^2\left(1\right)\right]\right)+0.1598\left(X^2\left(2\right)-M\left[X^2\left(2\right)\right]\right)+\hfill \\ & +0.0009\left(X^2\left(3\right)-M\left[X^2\left(3\right)\right]\right)+0.0395\left(X^2\left(4\right)-M\left[X^2\left(4\right)\right]\right)-\hfill \\ & -0.0065\left(X^2\left(5\right)-M\left[X^2\left(5\right)\right]\right)-6.02\times 10^{-5}\left(X^2\left(6\right)-M\left[X^2\left(6\right)\right]\right)-\hfill \\ & -0.0075\left(X^3\left(1\right)-M\left[X^3\left(1\right)\right]\right)-0.0195\left(X^3\left(2\right)-M\left[X^3\left(2\right)\right]\right)-\hfill \\ & -2.524\times 10^{-6}\left(X^3\left(3\right)-M\left[X^3\left(3\right)\right]\right)+0.0833\left(X^3\left(4\right)-M\left[X^3\left(4\right)\right]\right)+\hfill \\ & +0.0027\left(X^3\left(5\right)-M\left[X^3\left(5\right)\right]\right)+2.304\times 10^{-7}\left(X^3\left(6\right)-M\left[X^3\left(6\right)\right]\right);\hfill \end{array}$$

(13)

The preceding crop is rapeseed

$$\begin{array}{cc}{X}_p\left(7\right)=m_p^{\left(6.3\right)}\left(1.7\right)\hfill & =5.044+0.0047\left(X\left(1\right)-M\left[X\left(1\right)\right]\right)-0.0353\left(X\left(2\right)-M\left[X\left(2\right)\right]\right)+\hfill \\ & +0.0794\left(X\left(3\right)-M\left[X\left(3\right)\right]\right)+0.0105\left(X\left(4\right)-M\left[X\left(4\right)\right]\right)-\hfill \\ & -0.1734\left(X\left(5\right)-M\left[X\left(5\right)\right]\right)+0.1192X\left(X\left(6\right)-M\left[X\left(6\right)\right]\right)-\hfill \\ & -0.462\left(X^2\left(1\right)-M\left[X^2\left(1\right)\right]\right)+0.2187\left(X^2\left(2\right)-M\left[X^2\left(2\right)\right]\right)+\hfill \\ & -0.0008\left(X^2\left(3\right)-M\left[X^2\left(3\right)\right]\right)+0.0336\left(X^2\left(4\right)-M\left[X^2\left(4\right)\right]\right)-\hfill \\ & -0.0091\left(X^2\left(5\right)-M\left[X^2\left(5\right)\right]\right)+0.0002\left(X^2\left(6\right)-M\left[X^2\left(6\right)\right]\right)-\hfill \\ & -0.0926\left(X^3\left(1\right)-M\left[X^3\left(1\right)\right]\right)-0.0288\left(X^3\left(2\right)-M\left[X^3\left(2\right)\right]\right)\hfill \\ & +2.96\times 10^{-6}\left(X^3\left(3\right)-M\left[X^3\left(3\right)\right]\right)+0.0704\left(X^3\left(4\right)-M\left[X^3\left(4\right)\right]\right)-\hfill \\ & -0.0209\left(X^3\left(5\right)-M\left[X^3\left(5\right)\right]\right)-7.0703\times 10^{-7}\left(X^3\left(6\right)-M\left[X^3\left(6\right)\right]\right).\hfill \end{array}$$

(14)

In Formulas (13) and (14), the first, second, and third initial moments of random variables $X\left(i\right), i=\overline{1.6}$ have the following values:

$$\begin{array}{c}M\left[X\left(1\right)\right]=3.706;M\left[X^2\left(1\right)\right]=13.779;M\left[X^3\left(1\right)\right]=51.56;\hfill \\ M\left[X\left(2\right)\right]=4.5;M\left[X^2\left(2\right)\right]=\mathrm{20,417};M\left[X^3\left(2\right)\right]=93.76;\hfill \\ M\left[X\left(3\right)\right]=90;M\left[X^2\left(3\right)\right]=8700;M\left[X^3\left(3\right)\right]=\mathrm{891,000};\hfill \\ M\left[X\left(4\right)\right]=0.75;M\left[X^2\left(4\right)\right]=0.625;M\left[X^3\left(4\right)\right]=0.5625;\hfill \\ M\left[X\left(5\right)\right]=6.5;M\left[X^2\left(5\right)\right]=42.5;M\left[X^3\left(5\right)\right]=279.5;\hfill \\ M\left[X\left(6\right)\right]=356.33;M\left[X^2\left(6\right)\right]=\mathrm{127,070};M\left[X^3\left(6\right)\right]=\mathrm{45,348,636.33}\hfill \end{array}$$

.

3.5. Verification of the Developed Model

We determined the approximation error. Mean relative errors were analyzed for three approximations: a linear model, a second-order polynomial model, and a third-order polynomial model. The linear model had the highest mean relative error of 11.9942%. The second-order polynomial model showed a mean error of 4.8582%. Therefore, the above models are not recommended for use.

We revealed that the mean relative error depends on the preceding crop. When the preceding crop was maize, the mean relative error for a third-order polynomial model was 1.7884%. Its value varied from 0.0056% to 7.2168%. The standard deviation was 1.4691. The coefficient of variation was rather high, and it was equal to 82.15%. If the preceding crop was rapeseed, the mean relative error was higher, 2.7532%. The standard deviation was almost the same (1.4532). The coefficient of variation was 52.78%. Hence, the errors varied in a wide range. Although, the mean relative error was of low value. Thereby, the developed third-order polynomial model was proven to predict winter wheat yield.

The maximum relative error of 7.61% (the preceding crop was a rapeseed) was under the following conditions: plant density—5.0 million per hectare, the effective temperature sum of the autumn vegetation—345°C, micro-fertilizer—scheme II, irrigation scheme—600 + 1000, preceding crop—rapeseed, nitrogen—120 kg/ha. The minimum relative error of 0.021% was under the following conditions: plant density—4.5 million per hectare, the effective temperature sum of the autumn vegetation—369 °C, micro-fertilizer—scheme II, irrigation scheme—700 + 900, preceding crop—rapeseed, nitrogen—120 kg/ha (Figure 4). The use of another preceding crop (maize) changed relative errors. The maximum relative error of 7.22% was at the plant density of 4.5 million per hectare and the effective temperature sum of the autumn vegetation of 369 °C. The minimum relative error of 0.0056% was observed at the plant density—4.0 million per hectare, the effective temperature sum of the autumn vegetation—369 °C, and the second scheme of micro-fertilizer (Figure 5).

To further verify the prediction ability of the developed mathematical model, we applied it to five farms in the Mykolaiv region in 2021. Their actual yields ranged from 5025 to 5284 kg/ha. The relationship between the prediction and observation was found to validate the accuracy of our model. Results are presented in

Table 10. Predicting for farms.

farm	preceding crop	nitrogen, kg/ha	plant density, mln/ha	the effective temperature sum of the autumn vegetation, °C	micro-fertilizer	irrigation scheme	actual yield, kg/ha	forecast, kg/ha	error, %
1	maize	120	5.0	355	I	I	5244	5123	2.31
2	maize	115	5.0	355	I	I	5017	5108	1.81
3	maize	120	4.5	355	I	I	5180	5037	2.76
4	rapeseed	90	4.0	355	I	I	5015	4941	1.48
5	rapeseed	90	4.5	355	I	I	5284	5071	4.03

. The developed model produced an average accuracy of 2.48%. Errors varied from 1.48 to 4.03%.

4. Conclusions

Wheat is an important food crop. Climate change and an increase in world population has increased the importance of wheat yield forecasting. It is a principal problem for both farmers and authorities.

In this study, proposed a mathematical method by which to solve a significant practical problem of modeling winter wheat yields. The mathematical model was developed based on the results of a three-year field experiment. The mathematical forecasting model can use an arbitrary number of variables affecting the yield, preceding crops, fertilizers, plant density, an irrigation scheme, the effective temperature sum of the autumn vegetation, and micro-fertilizers. The structure and computational algorithm do not depend on the number of variables and the order of the nonlinear stochastic model.

The modeling method can use various functional dependences of yield on random factors (linearity, stationarity, Markovianity, monotonicity, etc.). The mathematical model uses weather and technological stochastic indicators. It allows us to achieve the best accuracy. The mean relative error does not exceed 2.75%; whereas a linear extrapolation gives the mean relative error of 9–12%. Therefore, the only nonlinear model is suitable for practical application.

Further studies are planned to develop a mathematical model for energy, environmental, and economic assessment of wheat cultivation as a function of agricultural practices and weather conditions. They will build on our previous papers [55,73,74].