# On-the-field simulation of fertilizer spreading: Part 3 – Control of disk inclination for uniform application on undulating fields

El Mehdi Abbou-Ou-Cherif, Emmanuel Piron, Alaa Chateauneuf, Denis Miclet, Sylvain Villette

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

El Mehdi Abbou-Ou-Cherif, Emmanuel Piron, Alaa Chateauneuf, Denis Miclet, Sylvain Villette. On-the-field simulation of fertilizer spreading: Part 3 – Control of disk inclination for uniform application on undulating fields. Computers and Electronics in Agriculture, 2019, 158, pp.150-158. 10.1016/j.compag.2019.01.050. hal-02067177

HAL Id: hal-02067177

https://institut-agro-dijon.hal.science/hal-02067177

Submitted on 30 Sep 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On-the-field simulation of fertilizer spreading: Part 3 – Control of disk inclination for uniform application on undulating fields

E-M. Abbou-Ou-Cherif<sup>a,b</sup>, E. Piron<sup>a</sup>, A. Chateauneuf<sup>b,c</sup>, D. Miclet<sup>a</sup>, S. Villette<sup>d,*</sup>

To improve the quality of fertiliser centrifugal spreading, several control devices have already been developed to manage some disruptions occurring on horizontal fields. However, controls which manage disruptions on non-flat fields due to changes in tractor attitude (pitch and roll) have not been developed. In this study, the design of a new control device is developed for a twin-disk spreader by considering two new degrees of freedom for each disk and controlling the longitudinal and lateral tilts of each disk. These tilt corrections are derived from solving a constrained optimization problem. The cost function is the weighted sum of squared differences between the travelled distances of the particles obtained in the considered topography and those that would have been obtained on a horizontal surface. The weighting coefficients are provided by the experimental angular mass flow distribution. In order to reduce the computation time and expect the use of the method for a real-time correction device, a simplified optimization problem is proposed. The method is assessed for three ground surface con-figurations: longitudinal slope break, side slope break and combined slopes. The simulation demonstrates that using the new control device the application errors are lower than 10% while they reach up to 40% without tilt correction. Moreover the study shows that the problem of optimization uniformity on non-flat fields can be solved computing the travelled distance for only a very few number of particles of mean diameter. This helps in reducing the computational time required to solve the optimization problem and in making possible the de-velopment of real-time control devices.

## 1. Introduction

In agriculture, centrifugal spreaders are widely used for fertilizer application, with up to 90% of the market share in Europe (Van Liedekerke et al., 2009). In the last decades, numerous sensors and control mechanisms have been developed and marketed to improve spreading accuracy and environmental compliance (Yule and Grafton, 2013). Various new technologies have been used such as on-board weighing, variable rate application device, border management device, variable width spreading device, headland management system, start and stop control device, or self-calibration system. For instance, modern centrifugal spreaders can adjust the fertilization rate using variable rate technologies (VRT), based on field information provided by geographical information systems (GIS). Drop point and flow rate control devices are also commonly used to adapt the shape of the spread pattern deposition to the two-dimensional geometry of the field. Using

simulations, Virin et al. (2006) already demonstrated that such control systems improve the accuracy of the spreading, keeping the applied rate in a range of +/- 10% of the local targeted application for flat fields. Nevertheless, when spreaders are used on non-horizontal surfaces, the accuracy of the spreading may be called into question.

Few works addressed the effect of spreader level on the fertilizer distribution uniformity. Parish (2003) studied the effect of the longitudinal tilt (front to rear angle) of the spreader and demonstrated that small changes such as 5 degrees caused significant pattern distortion. Yildirim (2008) studied the effect of the transverse tilt (side to side angle) and demonstrated that the fertilizer distribution pattern was skewed when the spreader was not horizontal. These studies highlighted that an automatic leveling system would be required to maintain an even fertilizer distribution when the ground surface was horizontal but they did not consider the spreading distortions in the case of non-flat fields. More recently, Grafton et al. (2015) pointed out the

<footnote><sup>a</sup>Irstea, UR TSCF, Centre de Clermont-Ferrand, Les Palaquins, F-03150 Montoldre, France</footnote>

<footnote><sup>b</sup>Clermont Université, Université Clermont Auvergne, Institut Pascal, BP 10448, F-63000 Clermont-Ferrand, France</footnote>

<footnote><sup>c</sup>CNRS, UMR 6602, Institut Pascal, F-63171 Aubière, France</footnote>

<footnote><sup>d</sup>Agroécologie, AgroSup Dijon, INRA, Univ. Bourgogne Franche-Comté, 26, bd Docteur Petitjean, F-21000 Dijon, France</footnote>

<footnote>* Corresponding author.
E-mail address: sylvain.villette@agrosupdijon.fr (S. Villette).</footnote>
negative effect of slope on the spreading quality because of a skewing of the spread pattern. Horrocks et al. (2015) also demonstrated that spreading fertilizer on humps and hollows affects the distribution patterns and the uniformity of fertilizer application depending on the bout width, and the driving position with respect to humps, slopes and hollows.

Considering the proportion of the arable land located on slopes and hilly regions, the proportion of field areas that could be affected by non-uniform fertilizer application cannot be ignored. For example, in the hilly region of Bend Subcarpathians, in Romania, Chandes et al. (2009) analysed that 52% of arable lands are situated on low slope lands and 18% flat slow lands.

Thus, in order to improve the quality of the spreading, a new active control device needs to be developed to ensure even fertilizer distribution on undulating fields regardless of their slope. To develop this new technical solution, a traditional approach would lead to costly or unrealistic experiments because of the diversity of slope situations and the difficulty of in-field measurements. Thus, virtual experiments based on numerical simulations were preferred to overcome the limits of in-field experiments. The use of hybrid approaches (Reumers et al., 2003; Villette et al., 2017) was of particular interest to combine some theoretical motion models with experimental data measured at some steps of the spreading process. In order to obtain realistic results, Villette et al. (2017) suggested the use of the Hybrid Centrifugal Spreading Model which combined a mechanistic approach based on the use of mechanical relationships and a stochastic approach based on the use of statistical distribution input parameters (derived from actual experiments). Nevertheless, the model described in Villette et al. (2017) was developed for flat ground conditions and did not take into account any disruptions due to the slope in the case of uneven or sloped ground (i.e. non-flat fields).

In order to study the effect of the slope on the quality of fertilizer spreading and design a correcting active control device, three main steps were required.

First, a new model was developed (Abbou-ou-cherif et al., 2017b) to simulate fertilizer spreading during the motion of a tractor on non-flat fields simulated by a digital elevation model. The coordinate system, the geometrical parametrization of tractor and spreader, and the ballistic flight model required for this study had been described by Abbou-ou-cherif et al. (2017b). In this first step, various experimental measurements were also carried out to provide the input data required by the model for ammonium nitrate. The particle size distribution was derived from a sieve analysis conducted in accordance with the standard ISO 8397 (1988). The horizontal outlet angle and vertical mass flow distribution was deduced from previous measurements (Villette et al., 2013). The horizontal angular mass flow distribution and the particle drag coefficient were deduced from 2D spread patterns measured with a rotating test bench (Piron et al., 2010) for a concave disk equipped with two radial vanes (radius: 395 mm). The model was then assessed and validated by comparing simulated and measured 2D spread patterns for three disk rotation speeds (600, 800 and 1000 rpm) and three disk tilt angles (0°, 5°, 10°).

Second, the model was used (Abbou-ou-cherif et al., 2017a) to compute the local application rate for various simulated 3D slope configurations, to quantify the resulting over- and under-application, and to understand the mechanisms of spread pattern distortions. The use of the hybrid model ensured the study of application uniformity in theoretical but realistic conditions. Various situations had been investigated such as regular longitudinal or side slope and irregular slope with longitudinal or side slope break (i.e. change of slope). Application rate maps were simulated and over- and under-application areas were characterized by computing the mean longitudinal or transverse application rates as follows:

$$Q_{\text{mean}_i} = \frac{\sum_{i-1}^{N} m_{ij}}{\frac{N}{N_T}} \times 100 \quad (1)$$

where i is the index of a virtual straight row of collecting trays of 1 × 1 m (in the longitudinal or transverse direction), $m_{ij}$ is the accumulated mass of particles in the jth collecting tray of the ith row, N is the constant number of collecting trays in every row, M is the total mass of particles in the considered area, $N_T$ is the total number of collecting trays in the same area and $Q_{\text{mean}_i}$ is the mean application rate in the ith row.

In the case of regular non-flat fields (with constant slope), regarding the mean transverse application rates, Abbou-ou-cherif et al. (2017a) computed application errors within ± 10% even at steep slopes. However, in the case of a longitudinal change of slope (for 10° positive slope) the deviation from the expected uniform rate reached a maximum of +45% and a minimum of -25% (deviation observed on the mean longitudinal application rate). These deviations from the targeted rate clearly demonstrated the need for correcting systems to keep the application variation within an acceptable range. Considering the level of the spreader, when working on a horizontal surface, Parish (2003) and Yildirim (2008) already suggested the design of an automatic levelling system to ensure an even fertilizer distribution. In the case of non-horizonal fields, the simulation results obtained by Abbou-ou-cherif et al. (2017a) not only help in understanding why the spread pattern is affected by the slope (in terms of tractor attitude and in terms of intersection between the ground surface and particle trajectories) but also provide help in designing appropriate correcting devices to prevent application errors.

Third, the model is now used to design a new correction system. Thus, the objective of this paper is to design an active control device for dual-disk centrifugal spreader to dynamically adapt the disk tilt to the site-specific slope and thus limit application errors when fertilizers are spread onto non-flat fields, regardless of the tractor path and ground slope configuration. The main contributions of the paper are: (1) the design of an active system to control the longitudinal and the lateral tilt of each spinning disk independently, (2) the design of a control strategy derived from an optimization problem which is adapted to reduce the computational time and consider the development of real-time control, (3) the assessment of the control strategy using numerical simulations.

## 2. Formulation of the optimization problem

Abbou-ou-cherif et al. (2017a) demonstrated that fertilizer application errors were particularly significant in the case of irregular non-flat fields (i.e. with longitudinal or transverse change of slope). This is due to the difference between the inclination of the tractor and the ground surface that leads to important deviations in the lengths of the particle ballistic flights. These deviations also take place on regular fields with longitudinal or side slope but to a less extent. For a given situation (e.g. a particular inclination of the tractor and of the ground surface) the value of these deviations is not the same in all the directions and depends on the angular direction of the particle around the disk (see Abbou-ou-cherif et al. (2017a) for details). These changes in the lengths of the particle ballistic flights directly affect the shape of the spread pattern deposition and, consequently, affect the application uniformity. Thus, it appears that the control of the length of the ballistic flights is required to achieve a return to a spread pattern deposition close to the one that would be obtained on a horizontal surface and thus to keep the application variation within an acceptable range. Therefore, a method of correcting application errors on non-flat field surfaces is developed in this paper by controlling the lengths of the ballistic flights of the particles during the spreading process.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Kinematic diagram showing the new revolute joints providing two additional degrees of freedom to each disk.</figcaption></figure>

## 2.1. Optimal control strategy

Changing particle flight lengths with current marketed centrifugal spreaders can be achieved through several means such as varying the disk velocity, the disk cone angle or the length of vanes. However, these parameters affect the flight lengths of all the particles in the whole spread pattern, and is not adapted to expect a differentiated correction depending on the angular direction of the particles around the spinning disk. As a consequence, a new control system has been developed in this study. It is based on the introduction of two new degrees of freedom for each disk (Fig. 1) so that its orientation can be varied in a three-dimensional space and specified in a spherical coordinate system. Thus, this orientation system can not only control the flight length of the particles, but also target the areas of application errors (i.e. a particular angular sector of spreading). Based on the system represented in Fig. 1, two control variables have been derived for each disk at each new position of the tractor on its trajectory: T (longitudinal tilt) and S (side or lateral tilt). For this purpose, a control strategy was implemented using the simulation model detailed in Abbou-ou-cherif et al. (2017a) and Abbou-ou-cherif et al. (2017b). This model computes the fertilizer rate applied on the field based on input data pertaining to the centrifugal spreader parameters, the particle characteristics and the digital elevation model (DEM) of the field. A classical control strategy would consist of inverting this direct model in order to derive the values of the control variables from the required application rate. However, this task is difficult to achieve because the intersections of the particle ballistic flights with the DEM can only be computed numerically. As a consequence, the control strategy chosen in this study consists in deriving the control variables by solving an optimization problem where the objective of the cost function is to minimize the sum of squared deviations of particle travelled distances compared to distances that would be obtained on a flat field. The control strategy developed in the following for one disk is applied separately for both disks.

The purpose of a combined action of longitudinal and side (i.e lateral) tilt is to correct the changes of the particle travelled distances due to ground topography, with respect to the distances that would be obtained on a horizontal surface. As illustrated in Fig. 2, the travelled distance is defined by the curvilinear length that separates the disk center from the landing position of a particle on the DEM (Abbou-ou-cherif et al., 2017a). Considering a particle with a travelled distance $p_{0l}$ in a flat configuration, the same particle has a range $p_i$ on a non-flat field such as there is a non-zero error defined by:

$$ \delta_i = p_i - p_{0l} \qquad (2) $$

Then, the goal of the optimisation problem is to minimize the squared error $\delta_i^2$ by controlling the disk longitudinal and side tilts on

non-flat fields.

Considering the results presented in Abbou-ou-cherif et al. (2017a) and their analysis, it has been assumed that the application errors could be reduced by keeping particle travelled distances equal to those obtained on a flat field, regardless of the tractor orientation. This objective can be mathematically expressed by a cost function which needs to be minimized at each point of the tractor trajectory. The cost function is as follows:

$$ F_c(T, S) = \sum_{k=1}^{N} \delta^2(\varphi_k, d_k, v_k, \Omega_k, T, S) \qquad (3) $$

where N is the number of particles ejected by the spinning disk at the considered tractor location, $\varphi_k$ is the angular direction of the kth particle measured with respect to the opposite of the travel direction, $d_k, v_k$ and $\Omega_k$ are respectively the diameter, the outlet velocity and the vertical outlet angle (with respect to the disk plane) and the of the kth particle. The optimization output is a vector $\mathbf{x}^* = (T^*, S^*)$ with the optimal values of T and S to be used for the disk inclination. In practice, the value of $\varphi_k$ is in the interval $[-120^\circ, 120^\circ]$ where a single spread pattern can be located. These limits have been deduced from numerous experimental measurements of single spread patterns captured with the rotating test bench CEMIB (Piron and Miclet, 2005) used in Irstea (the National Research Institute of Science and Technology for Environment and Agriculture, France) to design and certify centrifugal spreaders.

## 2.2. Spreading configuration

Typical values have been chosen for the spreading parameters. The simulations were computed for a tractor speed of 12 km/h, a working width of 28 m and a flow rate of 130 kg/ha. The spinning disks were equipped of radial vanes with a radius of 395 mm and a vertical angle of 13.25°. The disk rotation speed was 800 rpm. Considering particle characteristics, the density was 1750 kg m⁻³, the drag coefficient was 0.54 and the mean diameter was 3.1 mm. The vertical and horizontal angular distributions of the flow rate around the disk were derived from experimental measurements (Abbou-ou-cherif et al., 2017b). These distributions were measured when the disk axle was vertical. Thus, in this work, it was assumed that these distributions (measured with respect to the disk) are still the same irrespective of the tilt of the disk.

## 2.3. Spatial discretization step

The numerical resolution of the optimization problem requires discretization in space along the tractor trajectory. With centrifugal spreaders, the fertilizer particles are discontinuously thrown from the disk. At each vane revolution, a set of particles is taken up by the vane, then undergoes a centrifugal acceleration and finally is thrown into the field when the particles reach the extremity of the vane. Thus, particles are periodically ejected from the spinning disk at each vane revolution. In order to take into account this discontinuous process, the coordinates of the particle landing points are computed after each disk revolution, at each position of the tractor along its trajectory. Considering the distance L crossed by the tractor, the distance $\Delta L$ travelled by the tractor during the time required for one revolution of the spinning disk is as follows:

$$ \Delta L = 2\pi \frac{v_i}{\omega_d} \qquad (4) $$

where $v_i$ (expressed in m s⁻¹) is the tractor velocity and $\omega_d$ is the rotational speed of the disk (expressed in rad s⁻¹).

Considering a tractor velocity of 12 km/h and a disk rotation speed of 800 rpm, the distance $\Delta L$ travelled by the tractor during one disk revolution is 0.25 m. This value is used in this study to discretize the tractor trajectory and thus define the trajectory points where the spread pattern needs to be computed.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Comparison of two particle travelled distances in the case of a flat field ($p_{01}$ and $p_{02}$) and in the case of a non-flat field ($p_1$ and $p_2$). (O,x,y,z) is a Cartesian coordinate system.</figcaption></figure>

## 2.4. Simplification of the cost function for real time computation

Deriving the optimal vector $x^*$ for each disk has to be done in real time at each point of the trajectory. However, the distance step$\Delta L = 0.25$ m is travelled in less than 0.075 s which is a time span significantly lower than the time necessary to simulate a single spread pattern, with at least 1000 particles (considering the usual working width, flow rate and velocity speed). As a consequence, the cost function in Eq. (3) needs to be simplified. Furthermore, to take into account the mechanical and dynamic limits of the actuators, the optimization problem has to be constrained.

The distortion of a single spread pattern can be represented by the displacements between the landing positions of the particles when they are spread on a horizontal surface and their landing positions on the non-flat field. For example, these displacements are illustrated in Fig. 3 in the case of a side slope break. Considering this example, two observations can be made. First, the travelled distance variation depends on the angular sector of spreading. In the case illustrated in Fig. 3a, the magnitude of the displacement vectors is zero for the particles ejected on the left side or at the rear of the tractor but the displacement can

reach 4 m for the particles ejected on the right side. Moreover, the value of the displacement magnitude changes very progressively as a function of the angular direction (e.g. in Fig. 3, the displacement magnitude increases from the rear to the right side). Second, similar displacement fields are obtained whether by considering all the particles of the single spread pattern (Fig. 3a), or considering only a selection of particles (Fig. 3b) chosen with only one diameter (equals to the mean diameter: 3.1 mm) and thrown out with a constant outlet velocity (i.e. equals to the mean value: 44.2 m s$^{-1}$) and a constant vertical outlet angle with respect to the disk plane (i.e. equals to the mean value: 8.7°).

Considering that the whole angular spread sector can be decomposed into a set of adjacent angular sectors of same angular width$\Delta\varphi$, the cost function described by Eq. (3) can be rewritten as follows:

$$ F_c(T, S) = \sum_{k=1}^{N_1} \delta^2(\varphi_{1k}, d_{1k}, v_{1k}, \Omega_{1k}, T, S) \\ + \sum_{k=1}^{N_2} \delta^2(\varphi_{2k}, d_{2k}, v_{2k}, \Omega_{2k}, T, S) \\ + \dots + \sum_{k=1}^{N_n} \delta^2(\varphi_{nk}, d_{nk}, v_{nk}, \Omega_{nk}, T, S) \quad (5) $$

<figure><img src="image_5.png" /><figcaption>Fig. 3. Displacement fields of a single spread pattern in the case of side slope break. (a) Randomly sampled particles were used, (b) Mean diameter particles and average ejection parameter values were used. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
$$F_c(T, S) = \sum_{j=1}^{n} \left( \sum_{k=1}^{N_j} \delta^2(\varphi_{jk}, d_{jk}, v_{jk}, \Omega_{jk}, T, S) \right) \qquad (6)$$

where $n$ is the total number of angular sectors, $N_j$ is the number of particles to be considered in the $j$th sector.

Considering that the angular directions of all the particles are included in the interval $[-120, 120°]$ (see Section 2.1), the value of $\varphi_{jk}$ is bounded as follows:

$$-120 + (j - 1) \times \Delta\varphi < \varphi_{jk} \leq -120 + j \times \Delta\varphi \qquad (7)$$

which can also be written:

$$\varphi_j < \varphi_{jk} \leq \varphi_{j+1} \qquad (8)$$

with:

$$\varphi_j = -120 + (j - 1) \times \Delta\varphi \qquad (9)$$

with $j$ from 1 to $n$.

As previously analysed (Fig. 3), the displacement computed for the whole range of particle characteristics can be estimated to the one computed with the mean characteristics of the particle (i.e. mean diameter and mean ejection parameters), the cost function can be simplified as follows:

$$F_c(T, S) = \sum_{j=1}^{n} \left( \sum_{k=1}^{N_j} \delta^2(\varphi_{jk}, \bar{d}, \bar{v}, \bar{\Omega}, T, S) \right) \qquad (10)$$

where $\bar{d}$, $\bar{v}$ and $\bar{\Omega}$ are the mean values of the particle diameter, the outlet velocity and the vertical outlet angle (with respect to the disk plane) respectively.

Moreover, the previous analysis (Fig. 3) shows that the displacements can be considered as constant within one angular sector. Thus Eq. (10) yields:

$$F_c(T, S) = \sum_{j=1}^{n} (N_j \times \delta^2(\varphi_j, \bar{d}, \bar{v}, \bar{\Omega}, T, S)) \qquad (11)$$

In order to express the cost function with respect to the proportion of fertiliser amount in each angular sector, a new expression of the cost function is derived as follows:

$$F_{cN}(T, S) = \frac{F_c(T, S)}{N} = \sum_{j=1}^{n} \left( \frac{N_j}{N} \times \delta^2(\varphi_{jk}, \bar{d}, \bar{v}, \bar{\Omega}, T, S) \right) \qquad (12)$$

$$F_{cN}(\mathbf{x}) = \sum_{j=1}^{n} \omega_j \delta_j^2(\mathbf{x}) \qquad (13)$$

where $\mathbf{x}$ is the vector $(T, S)$, $\delta_j^2(\mathbf{x})$ is computed for one particle considering the mean values of the particle diameter, the outlet velocity and the vertical outlet angle:

$$\delta_j^2(\mathbf{x}) = \delta^2(\varphi_{jk}, \bar{d}, \bar{v}, \bar{\Omega}, T, S) \qquad (14)$$

and the weighted coefficient $\omega_j$ expresses the proportion of the fertilizer rate in the $j$th angular sector:

$$\omega_j = \frac{N_j}{N} \qquad (15)$$

Thus, in Eq. (13), the squared error obtained for the $j$th angular sector is weighted by $\omega_j$ according to the fertilizer rate in the $j$th angular sector. Therefore, the weight of each angular sector in the single spread pattern distortion is taken into account. Using the density function of the horizontal angular distribution of particles around the disk, the coefficient $\omega_j$ can be calculated as the probability for a particle having an angular direction $\varphi$ to be included in a given angular sector. Using the notation defined in Eq. (9), $\omega_j$ is defined as follows:

$$\left\{ \begin{array}{l} \omega_j = P(\{\varphi_l < \varphi \leq \varphi_{j+1}\}) = F_a(\varphi_{j+1}) - F_a(\varphi_j) \\ F_a(\varphi) = \int_{-120}^{\varphi} f_a(\zeta) d\zeta \end{array} \right. \qquad (16)$$

where $P(\{\varphi_l < \varphi \leq \varphi_{j+1}\})$ is the probability of having $\varphi$ in the interval $[\varphi_{bi}, \varphi_{i+1}]$, $F_a$ is the cumulative distribution function and $f_a$ the probability density function of the particle horizontal angle.

In practice, $\omega_j$ is directly derived from the experimental measurement of the horizontal angular distribution which is used as input data for the model.

Regarding the simplified cost function defined in Eq. (13), the higher the number $n$ of angular sectors is, the closer the optimal vector $\mathbf{x}^*$ becomes to the value that would be obtained using the cost function defined in Eq. (3). However, as illustrated in Fig. 3b, the trend of the displacements is a gradual and continuous change with respect to the angular direction. As a consequence, the choice of a small number of angular sectors and thus a small number of particles can be sufficient to carry out the optimization at each point of the tractor trajectory. In practice, the choice of 3 angular sectors is sufficient to provide good results as shown in Section 3.

## 2.5. Optimization constraints

To take into account the mechanical and dynamic limits of the actuators, several constraints have to be expressed as lower and upper bounds concerning the variables that have to be optimized (i.e. $T$ and $S$) and their time derivatives. The set $A$ of constraints is defined as follows:

$$\left\{ \begin{array}{l} L_{b1} \leq T \leq U_{b1} \\ L_{b2} \leq S \leq U_{b2} \\ L_{b3} \leq \frac{dT}{dt} \leq U_{b3} \\ L_{b4} \leq \frac{dS}{dt} \leq U_{b4} \end{array} \right. \qquad (17)$$

where $L_{bi}$ and $U_{bi}$ are respectively the lower and upper boundary defining the constraint $C_i$.

As the trajectory of the tractor is discretized at every distance step $\Delta L$ (Eq. (4)), the constraints $C3$ and $C4$ are rewritten to be discretized by taking into account the tractor velocity $v_t$. For example, the constraint $C3$ is modified as follows:

$$\frac{L_{b3}}{v_t} \leq \frac{dT}{dt} \cdot \frac{1}{v_t} \leq \frac{U_{b3}}{v_t} \qquad (18)$$

$$\frac{L_{b3}}{v_t} \leq \frac{dT}{dt} \cdot \frac{dt}{dl} \leq \frac{U_{b3}}{v_t} \qquad (19)$$

where $l$ and $t$ are respectively the tractor displacement and the time in the continuous domain.

Thus, in the discrete domain, the following relationships are established:

$$\frac{L_{b3}}{v_t} \leq \frac{T_k - T_{k-1}}{\Delta L} \leq \frac{U_{b3}}{v_t} \qquad (20)$$

$$\frac{L_{b3}}{v_t} \Delta L + T_{k-1} \leq T_k \leq \frac{U_{b3}}{v_t} \Delta L + T_{k-1} \qquad (21)$$

where $k$ is the index of the tractor position, $T_{k-1}$ and $T_k$ are two successive values of $T$ regarding the spatial discretisation.

Similarly, considering the constraint $C4$, bounds are established for the variable $S$ at the location $k$ as a function of its value found in the previous position $k-1$. The final formulation of the constrained optimization problem is given by:
$$ \mathbf{x}_k^* = \min_{x \in A} \left[ \sum_{j=1}^{n} \omega_j \delta_j^2(\mathbf{x}) \right] $$

$$ A = \left\{ \begin{array}{ll} L_{b1} \le T_k \le U_{b1} & (C1) \\ L_{b2} \le S_k \le U_{b2} & (C2) \\ \frac{L_{b3}}{v_t} \Delta L + T_{k-1} \le T_k \le \frac{U_{b3}}{v_t} \Delta L & (C3) \\ \frac{L_{b4}}{v_t} \Delta L + S_{k-1} \le S_k \le \frac{U_{b4}}{v_t} \Delta L & (C4) \\ + T_{k-1} & \\ + S_{k-1} & \end{array} \right\} \quad (22) $$

For numerical computation, the tilt angles of the actuators were limited at $\pm 30^{\circ}$ ($L_{b1} = L_{b2} = -30^{\circ}$, $U_{b1} = U_{b2} = +30^{\circ}$) and the rotational speeds were limited at $\pm 4^{\circ}/s$ ($L_{b3} = L_{b4} = -4^{\circ}/s$, $U_{b3} = U_{b4} = +4^{\circ}/s$).

## 2.6. Optimization problem solving and algorithm

The optimization problem was solved using the function `fmincon` of the optimization toolbox of the software Matlab (2005). This function allows finding a local optimum thanks to several evaluations of the cost function using the simulation model of the ballistic flights. In order to improve the computation time, an analytical expression of the cost function is derived at each new position of the trajectory. For this purpose, it is proposed to fit a third-order polynomial surface to a fixed number of points of the 2D cost function at the trajectory index $k$. These points are in the vicinity of the optimal solution $\mathbf{x}_{k-1}^*$ at the previous position of the trajectory. The polynomial surface coefficients are obtained by the method of least squares to solve the following linear system:

$$ \begin{bmatrix} 1 & S_p & T_p & S_p^2 & T_p^2 & S_p T_p & S_p T_p^2 & T_p S_p^2 & S_p^3 & T_p^3 \\ [a] & [b] & [c] & [d] & [e] & [f] & [g] & [h] & [i] & [j] \end{bmatrix} = [F_{CN}(T_p, S_p)] \quad (23) $$

where $p$ is the $p$th point of fit (with $p$ from 1 to $n_p$ the number of points used to locally fit the polynomial surface). In this work, $n_p$ was set at 10 points.

In order to calculate the cost function necessary to derive $\mathbf{x}^*(k+1)$, the points required to fit the cost function were then translated to the vicinity of the previously found solution $\mathbf{x}^*(k)$. Once the analytical expression of the polynomial surface was known, its gradient and Hessian matrix were computed and provided as additional arguments to the function `fmincon`. This approach ensured finding systematically a solution for the optimization problem at each point of the tractor trajectory and in a limited time frame. Moreover, the convexity of the cost function ensured that each solution was the optimal solution of the optimization problem.

## 3. Results and discussions

To assess the effectiveness of the control strategy based on the disk inclination, simulations were carried out for three ground surface configurations (i.e. longitudinal slope break, side slope break and undulating surface with combined slopes) using three discrete elevation models. These are presented in Fig. 4 (upper row). The investigated configurations are representative of the elevation changes of agricultural fields regarding the longitudinal and transverse slope with respect to the tractor travel direction. In each case, the fertilizer application maps were computed with (Fig. 4, lower row) and without (Fig. 4, middle row) the disk inclination control strategy. Considering longitudinal and side slope breaks, the digital elevation models used a horizontal surface and a $10^{\circ}$ tilted surface, connected by a surface with a radius of curvature of 100 m. The undulating elevation model (with combined slopes) simulated undulating hills and valleys with slopes set between $-6^{\circ}$ and $+6^{\circ}$. All these parameter values ensured the realistic nature of the simulations as explained in Abbou-ou-cherif et al. (2017a). When the control strategy is used, the application map presented in Fig. 4 (lower row) results from the use of only 3 angular sectors ($n=3$) in the optimization problem. For the three discrete elevation models, comparing the middle and the lower row of Fig. 4 shows that, over-application and under-applications are much lower when the disk inclination control system is used. Thus, the corrections of the disk tilts significantly improve the uniformity of the fertilizer application for the three topographic reliefs simulated.

Fig. 5 presents the longitudinal and transverse application rates computed using Eq. (1) considering the areas delineated in Fig. 4 for the three ground surface configurations. When the control strategy was used, the mean application rates were computed considering various number $n$ of angular sectors for the optimization problem. In all cases, Fig. 5 shows that the application error magnitudes are significantly reduced since they are within the interval $\pm 10\%$ when the control strategy is used. Furthermore, Fig. 5 shows that choosing $n=3$ provides good results and no uniformity enhancement is obtained by increasing the number of angular sectors above $n=6$.

The analysis of the case of longitudinal slope break helps in understanding how the control strategy acts on the disk orientation. In this case, the difference between the tractor and the ground surface inclination varies in only one direction. Fig. 6 shows the curves of the optimal longitudinal ($T$) and side ($S$) tilt signals computed for each disk along the tractor trajectory when 3 angular sectors are considered ($n=3$). In this slope situation, the longitudinal tilt signal is the most significant with a maximum value up to $6^{\circ}$. The longitudinal tilt signal is the same for both disks. The side tilt signals are limited to small values less than $1^{\circ}$. This is due to the fact that in the case of a longitudinal slope break, compensation required to keep the particle travelled distances unchanged are more important in the longitudinal direction than in the transverse direction. Furthermore, the side tilt correction curves of both disks are symmetrical because when the right disk is tilted to the left, the left disk is tilted to the right. The longitudinal tilt signal has a bell shape adjusting the orientation of the disk with regards to the inclination of the ground surface. This shape corresponds to the inverted profile of the mean travelled distance of the particles presented in Abbou-ou-cherif et al. (2017a). In the area delimited by $x \le 5$ m, no control is active ($T=0$ and $S=0$) because the field is a horizontal surface. In the area delimited by $x \ge 30$ m, the tractor and the ground surface have the same inclination. Consequently, the longitudinal tilt control is active but with an absolute value less than $1^{\circ}$. The small negative value observed in this situation corresponds to the reduction of the vertical orientation of the particle outlet velocity required to reduce the travelled distances of the particles ejected at the rear of the spreader.

The previous results show that the control strategy based on disks inclination is effective in reducing the application errors under acceptable limits in various slope situations. They provide an overall validation of the control method based on the correction of the particle ballistic flight lengths and of the simplifying assumptions made to solve the optimization problem and reduce the computational complexity. In
<figure><img src="image_6.png" /><figcaption>Fig. 4. Digital elevation models (upper row) and corresponding fertilizer application maps simulated without (middle row) or with disk inclination control (lower row) for a 10° longitudinal change of slope (left column), a 10° transverse change of slope (middle column) and an undulating surface with combined slopes set between -6° and +6° (right column). White arrows on the discrete elevation models represent the adjacent passes of the tractor. The rectangular edge drawn on the applications maps delineate the collection areas that are analyzed in Fig. 5. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

particular, the study demonstrates that adequate solutions are obtained by using the mean values of particle diameter and ejection parameters; considering only a reduced number of angular sectors; and fitting the 2D cost function through a limited number of points. This is of particular interest to reduce the computation time and thus expect the development of an actual control device working in real-time.

Nevertheless, using the Matlab programming language and a 3.4 GHz quad-core processor, the optimal orientation of the disks could only be derived in less than 0.3 s. This is approximately four times the computation time required for real-time use considering a tractor speed of 12 km/h and a trajectory sampling of 0.25 m. The computation time could be reduced by improving the hardware and using a low level
<figure><img src="image_7.png" /><figcaption>Fig. 5. Mean application rates, with and without disk inclination control, calculated inside the rectangular areas shown on the application maps of Fig. 4 for the 10° longitudinal change of slope (a), the 10° transverse change of slope (b) and an undulating surface with combined slopes set between – 6° and +6° (c). The mean application rates were computed with various number *n* of angular sectors used in the optimization problem. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

programming language. In order to reduce computation time, the increase of the sampling step of the tractor trajectory could also be considered for further studies.

In this study, the vertical and horizontal mass flow distribution of fertilizer around the disk (and measured experimentally with respect to the disk) was assumed to be unchanging irrespective of the tilt of the disk. Further research is needed to measure potential changes in the mass flow distribution around the disk due to potential changes of the drop point position on the disk or over potential changes of fertilizer motion in the vanes. If the mass flow distribution around the disk changed depending on the disk tilt, then some relationships would be established to describe these changes and introduced in the mechanical model to take into account this behavior.

## 4. Conclusion

To sustain fertilizer application uniformity despite disruptions due to the slope and its variations in agricultural field, a new concept of active control of the disk attitude was developed. The control system consisted in offsetting particle travelled distance deviations by changing the inclination of the spinning disk axle. Considering two new degrees of freedom for the disk, the optimal longitudinal and side (i.e. lateral) tilts were derived from a constrained optimization problem. The objective was to minimize the quadratic sum of all the particle travelled

distance deviations at each position of the tractor trajectory in the field (i.e. minimizing the variance of the differences between the travelled distances of the particles obtained in the considered topography and the distances that would have been obtained on a horizontal surface). This cost function was simplified by considering a reduced number of angular sectors in the single spread patterns where the flight distance of only one particle weighted by a coefficient can be taken into account. This simplification reduced the computation time and paved the way for the development of real-time control devices. Using a simulation model, three cases of longitudinal, lateral and combined slopes were investigated. Simulations showed that the tilt correcting control reduces effectively the magnitudes of fertilizer application errors, below the acceptable limit of ± 10%, while the application errors could reach 40% when the application maps were simulated without control strategy. Furthermore, the study showed that the choice of a number of angular sectors less than or equal to six is enough to reach relevant disk tilt corrections. Accordingly, simulation results demonstrated the overall performance of the method to achieve adequate uniformity on undulating fields. In practice, the development of control devices could benefit from the development of on-board sensors that are useful to provide model parameters such as the outlet velocity or the mass flow distribution around the disk, in order to reduce preliminary experimental measurements or control these parameters with a closed loop system.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Control signals (longitudinal tilt T in red and side tilt S in green) computed in the case of the longitudinal slope break of 10° (presented in Fig. 4.) (a) Right disk, (b) Left disk. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
## Acknowledgment

This work was supported by the French National Research Institute of Science and Technology for Environment and Agriculture (Irstea). We would also like to acknowledge the grant of the European Regional Development Fund (FEDER) through the Auvergne Region of France. Finally, we would like to thank the company RAUCH Landmaschinenfabrik GmbH for its support.

## References

Abbou-ou-cherif, E.-M., Piron, E., Chateauneuf, A., Miclet, D., Lenain, R., Koko, J., 2017a. On-the-field simulation of fertilizer spreading: Part 2 - Uniformity investigation. Comput. Electron. Agric. 141 (Part C), 118-130.

Abbou-ou-cherif, E.M., Piron, E., Chateauneuf, A., Miclet, D., Lenain, R., Koko, J., 2017b. On-the-field simulation of fertilizer spreading: Part 1 – Modeling. Comput. Electron. Agric. 142 (Part A), 235-247.

Chendes, V., Dumitru, S., Simota, C., 2009. Analyzing the landforms-agricultural land-use types relationship using a DTM-based indicator. Scientific Papers, USAMV, Series A Agronomy 52, 135-140.

Grafton, M.C., Izquierdo, D.A., Yule, I.J., Willis, L.A., Manning, M.J., 2015. The effect of field conditions on in-field spread patterns from twin disk spreaders. 2015 ASABE Annual International Meeting. ASABE, St. Joseph, MI.

Horrocks, A., Thomas, S., Tregurtha, C., Meenken, E., Horrell, R., 2015. Developing guidelines for fertiliser spreading on west coast humps and hollows. In: 28th Annual Fertilizer and Lime Research Centre Workshop. Massey University, Palmerston North, New Zealand.

ISO 8397, 1988. Solid Fertilizers and soil Conditioners - Test Sieving. International

Organization for Standardization, Genève.

Matlab, 2005. Version 7.1. The MathWorks, Inc., Natick, Massachusetts, United States.

Parish, R.L., 2003. Effect of impeller angle on pattern uniformity. Appl. Eng. Agric. 19, 531-533.

Piron, E., Miclet, D., 2005. Centrifugal fertiliser spreaders: a new method for their evaluation and testing. In: Conference of the International Fertiliser Society, York, United Kingdom. The International Fertiliser Society, pp. 22p.

Piron, E., Miclet, D., Villette, S., 2010. CEMIB: an innovative bench for spreader eco-design. AgEng 2010, International Conference on Agricultural Engineering, Clermont-Ferrand, France.

Reumers, J., Tijskens, E., Ramon, H., 2003. Experimental characterisation of the cylindrical distribution pattern of centrifugal fertiliser spreaders: towards an alternative for spreading hall measurements. Biosyst. Eng. 86, 431-439.

Van Liedekerke, P., Tijskens, E., Dintwa, E., Rioual, F., Vangeyte, J., Ramon, H., 2009. DEM simulations of the particle flow on a centrifugal fertilizer spreader. Powder Technol. 190, 348-360.

Villette, S., Piron, E., Martin, R., Miclet, D., Jones, G., Paoli, J., Gée, C., 2013. Estimation of two-dimensional fertiliser mass flow distributions by recording granule impacts. Biosyst. Eng. 115, 463-473.

Villette, S., Piron, E., Miclet, D., 2017. Hybrid centrifugal spreading model to study the fertiliser spatial distribution and its assessment using the transverse coefficient of variation. Comput. Electron. Agric. 137, 115-129.

Virin, T., Koko, J., Piron, E., Martinet, P., Berducat, M., 2006. Application of optimization techniques for an optimal fertilization by centrifugal spreading. IEEE Proceedings of the International Conference on Intelligent Robots and Systems, Beijing, China.

Yildirim, Y., 2008. Effect of side to side spreader angle on pattern uniformity in single- and twin-disc rotary fertilizer spreaders. Appl. Eng. Agric. 24, 173-179.

Yule, I., Grafton, M., 2013. New spreading technologies for improved accuracy and environmental compliance. In: Currie, L.D., Christensen, C.L. (Eds.), Accurate and Efficient use of Nutrients on Farms. Fertilizer and Lime Research Centre, Massey University, Palmerston North, New Zealand, pp. 1-13.