Development and Evaluation of Drone Based Spraying System for Precision Agriculture Application

4.1 Translation motion

$\begin{gathered}{\left[\begin{array}{c}\dot{x} \\ \dot{y} \\ \dot{z}\end{array}\right]=\left[R_B^I\right]\left[\begin{array}{c}u \\ v \\ w\end{array}\right]} =\left[\begin{array}{ccc}c \psi c \theta & s \theta c \psi s \varphi-s \psi c \varphi & c \theta s \psi c \varphi-s \psi s \varphi \\ s \psi c \theta & s \theta s \psi s \varphi+c \psi c \varphi & s \theta s \psi c \varphi-c \psi s \varphi \\ -s \theta & c \theta s \varphi & c \theta s \varphi\end{array}\right]\left[\begin{array}{c}u \\ v \\ w\end{array}\right]\end{gathered}$     (5)

The linear velocity equations in the inertial frame are:

$\begin{gathered}\dot{x}=(c \psi c \theta) u+(s \theta c \psi s \varphi-s \psi c \varphi) \\ +(c \theta s \psi c \varphi-s \psi s \varphi)\end{gathered}$     (6)

$\begin{gathered}\dot{y}=(s \psi c \theta) u+(s \theta s \psi s \varphi+c \psi c \varphi) v \\ +(s \theta s \psi c \varphi-c \psi s \varphi) w\end{gathered}$     (7)

$\dot{z}=-(s \theta) u+(c \theta s \varphi) v+(c \theta s \varphi) w$     (8)

4.2 Rotational motion

$\begin{gathered}\omega_B=T_E^B \times \omega_E \\ {\left[\begin{array}{l}p \\ q \\ r\end{array}\right]=\left[\begin{array}{ccc}1 & 0 & -\sin \theta \\ 0 & \cos \psi & \sin \psi \cos \theta \\ 0 & -\sin \psi & \cos \psi \cos \theta\end{array}\right] \times\left[\begin{array}{c}\dot{\varphi} \\ \dot{\theta} \\ \dot{\psi}\end{array}\right]}\end{gathered}$     (9)

Euler angle rates in body frame can be expressed as $(p q r)$ and inertial frame can be expressed as $\{\dot{\varphi}, \dot{\theta}, \dot{\psi}\}$. The set of angular equations of motion in the inertial frame can be written as:

$\left[\begin{array}{c}\dot{\varphi} \\ \dot{\theta} \\ \dot{\psi}\end{array}\right]=\left[\begin{array}{c}p+q(\sin \varphi \tan \theta)+r(\cos \varphi \tan \theta) \\ q(\cos \varphi-r \times \sin \varphi) \\ q(\sin \varphi / \cos \theta+r * \cos \varphi / \cos \theta\end{array}\right]$

4.3 Newton – Euler equations

The Newton-Euler equations are mathematically represented in the body frame as follows:

$m \dot{V}_B+m \omega_B \times V_B=\sum F_B$     (10)

$I \dot{\omega}_B+\omega_B \times I \omega_B=\sum M_B$     (11)

$m$: The mass of rigid body.

$I$: The inertia matrix.

$\sum F_B$: The external Forces.

$\sum M_B$: The external moments.

$\sum \mathrm{F}_{\mathrm{B}}=\left[\begin{array}{c}F_x \\ F_y \\ F_z\end{array}\right], \sum M_B=\left[\begin{array}{l}\mathrm{M}_{\mathrm{x}} \\ \mathrm{M}_{\mathrm{y}} \\ \mathrm{M}_{\mathrm{z}}\end{array}\right]$

4.4 Forces equations

The forces are described with the help of Eq. (10).

$\left[\begin{array}{c}F_x \\ F_y \\ F_z\end{array}\right]=m \times\left[\begin{array}{c}\dot{u} \\ \dot{v} \\ \dot{w}\end{array}\right]+\left[\begin{array}{c}p \\ q \\ r\end{array}\right] \times m\left[\begin{array}{c}u \\ v \\ w\end{array}\right]$

The linear motion equations can be expressed in terms of linear accelerations in the body frame by rearranging the dynamical equations, focusing on the accelerations along the x, y, and z axes and their relation to forces, resulting in Eq. (12).

$\left[\begin{array}{c}\dot{u} \\ \dot{v} \\ \dot{w}\end{array}\right]=\left[\begin{array}{c}\frac{F_x}{m}+r v-q w \\ \frac{F_y}{m}+(p w-r u) \\ \frac{F_z}{m}+(q u-p u)\end{array}\right]$     (12)

4.5 Moment equations

Simplified:

$\left[\begin{array}{l}M_x \\ M_y \\ M_z\end{array}\right]=\left[\begin{array}{c}L \\ M \\ N\end{array}\right]=\left[\begin{array}{l}I_{x x} \dot{p}+\left(I_{z z}-I_{y y}\right) q r \\ I_{y y} \dot{q}+\left(I_{x x}-I_{z z}\right) r p \\ I_{z z} \dot{r}+\left(I_{y y}-I_{x x}\right) p q\end{array}\right]$     (13)

To express the equations with a focus on angular accelerations, it is essential to isolate the angular accelerations on the left side of the equation. This approach yields the equations, as illustrated in the equations provided below:

$\dot{p}=\frac{1}{I_{x x}}\left[\left(I_{y y}-I_{z z}\right) q r+M_x\right]$     (14)

$\dot{q}=\frac{1}{I_{y y}}\left[\left(I_{z z}-I_{x x}\right) r p+M_y\right]$     (15)

$\dot{r}=\frac{1}{I_{z z}}\left[\left(I_{x x}-I_{y y}\right) q p+M_z\right]$     (16)

4.6 External forces

Various external forces impact the flight dynamics of a quadrotor. The collective effect of these forces determines the overall net force exerted on the quadrotor. Eq. (17) outlines the components that make up the total external forces acting on the quadrotor.

$\sum F_B=F_M+F_G+F_w$     (17)

where, $F_G$ is the gravitational force, and $F_M$ represents the force generated by the motors, and $F_w$ denotes the force exerted by the wind.

$F_M$ represents:

$F_M=\left[\begin{array}{c}0 \\ 0 \\ f_m\end{array}\right]$     (18)

where, the $f_m$ is the sum of generated by the four motors:

$f_m=f_{m 1}+f_{m 2}+f_{m 3}+f_{m 4}$     (19)

And $F_G, F_w$ are represented by the matrix below:

$F_G=R_I^B\left[\begin{array}{c}0 \\ 0 \\ m g\end{array}\right]$     (20)

$F_w=\left[\begin{array}{l}f_{w x} \\ f_{w y} \\ f_{w z}\end{array}\right]$     (21)

The external forces can be simplified by removing wind force and applying the rotation matrix to the gravity vector, as shown in Eq. (22).

$\left[\begin{array}{c}F_x \\ F_y \\ F_z\end{array}\right]=\left[\begin{array}{c}-m g(\sin \theta) \\ m g(\sin \varphi \cos \theta) \\ m(\cos \varphi \cos \theta)-f_m\end{array}\right]$     (22)

4.7 State equations

State-space representation is a robust method for describing dynamical systems, particularly beneficial for multiple input-multiple output (MIMO) systems. It consists of two main equations [22, 23]:

The state-space representation is composed of two primary equations.

$\dot{x}=A x+B u$     (23)

$y=C x+D u$     (24)

where, $u$: Control input vector, $x$: State victor, $y$: Output victor, $A$: System matrix, $B$: Input matrix, $C$: Output matrix, $D$: Feed forward matrix.

Typically, a quadcopter system encompasses twelve state variables and four control inputs. both the state variables and the control inputs can be described using the following representation.

$x$ vector and $u$ are the input vector can be written as:

$x=[\dot{x} \dot{y} \dot{z} \ddot{x} \ddot{y} \ddot{z} \dot{\psi} \dot{\theta} \dot{\varphi} \dot{p} \dot{q} \dot{r}]^T, U=\left[\begin{array}{c}U_1 \\ U_2 \\ U_3 \\ U_4\end{array}\right]$

4.8 Newton's equations of motion

Newton's equation of motion, originally formulated in the body frame, can be transposed to the inertial frame:

$m \dot{V}_I=T_B^I \sum F_I$     (25)

Forces include thrust $F_M$ and gravity $F_G$

$m \dot{V}_I=T_B^I\left[\begin{array}{c}0 \\ 0 \\ -\mathrm{U}_1\end{array}\right]+T_B^I T_I^B\left[\begin{array}{c}0 \\ 0 \\ m g\end{array}\right]=T_B^I\left[\begin{array}{c}0 \\ 0 \\ -\mathrm{U}_1\end{array}\right]+\left[\begin{array}{c}0 \\ 0 \\ m g\end{array}\right]$

4.9 Linear and angular velocities

The derivatives of the linear and angular velocities can be simplified:

$\dot{V}_I=\left[\begin{array}{c}\ddot{x} \\ \ddot{y} \\ \ddot{z}\end{array}\right]=\left[\begin{array}{c}\frac{-U_1}{m}(\cos \varphi \sin \theta \cos \psi+\sin \varphi \sin \psi) \\ \frac{-U_1}{m}(\cos \varphi \sin \theta \sin \psi-\sin \varphi \cos \psi) \\ g-\frac{U_1}{m}(\cos \varphi \cos \theta)\end{array}\right]$     (26)

Assuming small tilt angles close to zero: 

$\omega_B=\omega_E$

Taking the derivative of both sides’ yields:

$\begin{aligned} & \dot{\omega}_{\mathrm{E}}=\dot{\omega}_{\mathrm{B}} \\ & {\left[\begin{array}{c}\ddot{\varphi} \\ \ddot{\theta} \\ \ddot{\psi}\end{array}\right]=\left[\begin{array}{c}\dot{\mathrm{p}} \\ \dot{\mathrm{q}} \\ \dot{\mathrm{r}}\end{array}\right]}\end{aligned}$

We can construct the final form of the derivative of the Euler rates vector, $\dot{\omega}_E$.

$\left[\begin{array}{c}\ddot{\varphi} \\ \ddot{\theta} \\ \ddot{\psi}\end{array}\right]=\left[\begin{array}{l}\frac{1}{I_{x x}}\left[\left(I_{x x}-I_{z z}\right) q r+M_x\right] \\ \frac{1}{I_{y y}}\left[\left(I_{z z}-I_{x x}\right) r p+M_y\right] \\ \frac{1}{I_{z z}}\left[\left(I_{x x}-I_{y y}\right) q p+M_z\right]\end{array}\right]$     (27)

The equations derived from the modelling of a quadrotor include twelve state variables:

$\begin{gathered}\dot{x}=(c \psi c \theta) u+(s \theta c \psi s \varphi-s \psi c \varphi) v+(c \theta s \psi c \varphi-s \psi s \varphi) w \\ \dot{y}=(c \theta) u+(s \theta s \psi s \varphi+c \psi c \varphi) v+(s \theta s \psi c \varphi-c \psi s \varphi) w \\ \dot{y}=(c \theta) u+(s \theta s \psi s \varphi+c \psi c \varphi) v+(s \theta s \psi c \varphi-c \psi s \varphi) w\end{gathered}$

$\begin{gathered}\ddot{x}=\frac{-U_1}{m}(c \varphi s \theta c \psi+s \varphi s \psi) \\ \ddot{y}=\frac{-U_1}{m}(c \varphi s \theta s \psi-s \varphi c \psi) \\ \ddot{z}=g-\frac{U_1}{m}(c \varphi c \theta) \\ \dot{\psi}=p+t \theta(q s \varphi+r c \varphi) \\ \dot{\theta}=q c \varphi-r s \varphi \\ \dot{\varphi}=1 / c \theta(q s \varphi+r c \varphi) \\ \dot{p}=\frac{1}{I_{x x}}\left[\left(I_{x x}-I_{z z}\right) q r+M_x\right] \\ \dot{q}=\frac{1}{I_{y y}}\left[\left(I_{z z}-I_{x x}\right) r p+M_y\right] \\ \dot{r}=\frac{1}{I_{z z}}\left[\left(I_{x x}-I_{y y}\right) q p+M_z\right]\end{gathered}$

4.10 Linearization

Linearization simplifies system dynamics for control design. Using small angle approximations:

$\begin{gathered}\sin (\varphi, \theta, \Psi) \cong(\emptyset, \theta, \Psi), \cos (\varphi, \theta, \Psi) \cong(1,1,1) \\ \varphi \psi \cong 0, \theta \psi \cong 0\end{gathered}$

This assumption will replace the nonlinear expressed probability in the Eq. (12) and give the result that these equations are linear

Resulting linear equations:

$\begin{gathered}\dot{x}=u, \quad \dot{y}=v, \quad \dot{z}=w \\ \ddot{x}=-g \theta, \quad \ddot{y}=g \theta, \ddot{z}=-\frac{U_1}{m} \\ \dot{\psi}=p, \quad \dot{\theta}=q, \quad \dot{\varphi}=r \\ \dot{p}=\frac{1}{I_{x x}} U_2, \quad \dot{q}=\frac{1}{I_{y y}} U_3, \quad \dot{r}=\frac{1}{I_{z z}} U_4\end{gathered}$

State space model

$x=\left[\begin{array}{lll}x \ y \ z \ u \ v \ u \ p \ q \ r \ \varphi \ \theta \ \Psi\end{array}\right]^T$

And input vector can be defined as:

$u=\left[F M_x M_y M_z\right]^T$

And finally, from previous equations, A and B matrices can be constructed as in equations, the state space matrices for A and B can be states as:

$A=\left.\frac{\partial f(x, u)}{\partial x}\right|_{\frac{x=x}{u=\bar{u}}}=\left[\begin{array}{llll}0001000 & 0 & 0000 \\ 0000100 & 0 & 0000 \\ 0000010 & 0 & 0000 \\ 0000000 & -g & 0000 \\ 000000 g & 0 & 0000 \\ 0000000 & 0 & 0000 \\ 0000000 & 0 & 0100 \\ 0000000 & 0 & 0010 \\ 0000000 & 0 & 0001 \\ 0000000 & 0 & 0000 \\ 0000000 & 0 & 0000 \\ 0000000 & 0 & 0000\end{array}\right]$     (28)

$B=\left.\frac{\partial f(x, u)}{\partial u}\right|_{\frac{x=x}{u=\bar{u}}}=\left[\begin{array}{cccc}0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -1 / m & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 1 / I_{x x} & 0 & 0 \\ 0 & 0 & 1 / I_{y y} & 0 \\ 0 & 0 & 0 & 1 / I_{z z}\end{array}\right]$     (29)

The linearized quadcopter dynamics are expressed in state-space form Eq. (23).

$\dot{x}=A x+B u$

where, A and B are system and input matrices, respectively.

4.11 Controller design

The proposed model uses six PID controllers to enhance flight control, ensuring accurate regulation and stable flight path [24]. These controllers, including altitude, orientation, and position, are central to the model, allowing for easy implementation and tuning, making them an advantageous choice for quadrotor simulations.

$u(t)=K_p e(t)+K_I \int_0^t e(\tau) d \tau+K_D \frac{d}{d t} e(t)$     (30)

$K_p$: Proportional gain. $K_I$: Integral gain. $K_D$: Derivative gain. $e(\mathrm{t})$: The error term. $u(t)$: The control input.

Simulink diagram

First, define the diagram of the quadcopter designed in Simulink (Table 2), providing an overview of each cascade in this diagram.

Table 2. Parameters of the quadcopter designed

The Simulink diagram in Figure 11 illustrates a quadrotor design with six PID controllers for stability and responsiveness. These controllers manage position and altitude, ensuring the quadrotor follows a predefined trajectory and maintains a steady height. The position controller handles lateral movement, while the altitude controller adjusts vertical position. This setup allows for detailed analysis and optimization of the quadrotor's performance.

11.png

Figure 11. Schema of Simulink diagram

Trajectory

The trajectory is essential for simulating a drone's performance in agricultural missions. It involves the drone ascending to 5 meters and following a predefined path, moving 30 meters along the x-axis and 100 meters along the y-axis. This path provides clear guidance for evaluating the drone's real-world effectiveness as shown in Figure 12.

12.png

Figure 12. Trajectory given to follow

Position controllers

The position controller in a drone is crucial for maintaining precise positioning in three dimensions (Z Y X). It receives reference inputs for desired positions and calculates errors to adjust the drone's actuators, ensuring it maneuvers effectively to achieve the desired positions as shown in Figure 13.

13.png

Figure 13. Schema position controller

Altitude controller

To calculate errors between desired positions and actual measurements, adjusting actuators to achieve desired altitude and ensuring stability during flight as shown in Figures 14.

14.png

Figure 14. Schema altitude controller

4.12 Motor mixing

Motor mixing is a crucial process in multi-rotor drone control, determining power distribution among motors to achieve desired flight behavior as shown in Figure 15. It ensures stability by adjusting motor speeds based on control inputs like roll, pitch, and yaw angles.

15.png

Figure 15. Schema motor mixing

4.13 Quadcopter dynamics

In this section, the dynamics of a quadcopter are examined, with particular attention paid to input processing, output processing, and the six degrees of freedom (6DOF). It draws attention to the relationship that exists between the quadcopter's dynamics and the mathematical model and explains the importance that relationship is to understand and controlling flight behavior.

4.14 Result of simulation

The initial section presents trajectory tracking results for the quadrotor along the X, Y, and Z axes. This analysis compares the actual and reference trajectories to evaluate the control system's accuracy and effectiveness. Figure 16 is a visual representation of the system's behavior in each dimension.

This examination evaluates the quadrotor system's performance in following a prescribed trajectory and maintaining stable flight. The chosen trajectory, shown in Figures 16 and 17, was simulated using MATLAB. Figure 16 compares the desired path (in blue) with the actual path followed by the drone (in red).

161.png

162.png

Figure 16. Quadcopter dynamics

17.png

Figure 17. The drone's trajectory along the X, Y, and Z axes and trajectory followed by drone

Figure 17 illustrates the comparison between the reference (in red) and actual (in blue) positions of X, Y, and Z, as well as the yaw pitch and roll angle Figure 16, when employing a PID controller.

4.15 PID controller response analysis

This section assesses the effectiveness of PID controllers in quadrotor system trajectory tracking, focusing on their contribution to accuracy and maintaining the desired path. Visual demonstrations and simulation results are provided in Table 3 and Figure 18.

Table 3. PID controller gain values

181.png

182.png

183.png

(a)

185.png

186.png

187.png

(b)

Figure 18. Angle step input results and, position step input results and position step input results

Tuning the six parameters—X, Y, and Z for position control, and yaw, pitch, and roll for orientation control—ensures drone stability during flight. By adjusting these parameters within the PID controllers, the system effectively regulates the drone's position and orientation, minimizing oscillations and deviations from the desired trajectory.

Finally, the six PID controllers stabilize the drone, as evidenced by its ability to follow the prescribed path, simulating real-world conditions for spraying a square field. The 3D and 2D trajectory plots demonstrate how the drone accurately follows the given trajectory around the X, Y, and Z axes. The tuning of the PID controllers effectively regulates the drone's position and orientation, minimizing oscillations and deviations during flight.