Design of anti-load perturbation flight trajectory stability controller for agricultural UAV

The unmanned aerial vehicle(UAV) used in agricultural fields usually has a load, which affects stability of flight trajectory. This problem is a great technical challenge and critical issue which has fascinated many researchers in this area. This research paper presents a flight dynamics model for agricultural UAV under load condition. The robust T-S fuzzy control method is applied in the attitude angle control part, and the traditional PID controller is used in the position loop control part. Results of simulation depict that the flight trajectory of agricultural UAV can reach certain stability when the load is relatively small.

1 Introduction

UAVs are widely used in pesticide spraying, sowing, pollination and other agricultural fields. Compared with traditional agricultural machinery, UAV has the advantage of high precision and efficiency (Hu et al., 2022). However UAV is a nonlinear, strongly coupled, under-driven system (Wang and Zhou, 2021). In operating processes, load of UAV changes in real time, which affects speed and accuracy of its flight attitude and trajectory stability (Luo et al., 2020). To improve the flight trajectory stability of UAV against load perturbations for different mission requirements has been a research focus for researchers (Qi and Ping, 2018).

A variety of algorithms are designed to control the flight trajectory stability of UAV,such as PID (Li et al., 2016), linear quadratic regulator (LQR) (Dong et al., 2015), H-infinity (Emam and F Ak.Harian, 2016), sliding mode control (Zhou et al., 2016), back-stepping (Shi and Lin, 2018), and adaptive control (Fan D et al., 2018) and so on. Specifically, A PID intelligent fuzzy control system is designed for spraying operation of agricultural UAV (Zheng and Chen, 2021). An on-line trajectory planning method based on the reinforcement learning is designed for UAV trajectory stability control against swing motion of the slung-payload (Xian et al., 2021). A double closed-loop adaptive control strategy is designed to solve the parameter uncertainty and external disturbance in trajectory tracking process of quad-rotor UAV problem (Gao et al., 2021). The adaptive sliding mode control strategy combined with self-anti-perturbation technique is used, considering the effects on flight performance due to sudden load changes caused by loading/unloading and external perturbation in short-distance transportation (An et al., 2018).

In summary, researches on trajectory control of UAV against load perturbation are remarkable, but the previous researches mostly focus on controlling the suspension load swing perturbation, and some of the research still has improving space of trajectory tracking considering the sudden change of load. Therefore, we designed a different controller for agricultural UAV flight trajectory control against load perturbation. In the attitude angle part, robust T-S fuzzy control method is used. And in the position loop, PID controller method is used. The proposed control method is verified by simulation tests. Load change is simulated by mass change of UAV. Results are showed by outputting UAV flight trajectory tracking image and the changes of roll, pitch and yaw angle.

2 Dynamics modelling

To describe the position and attitude of UAV system, two reference coordinate systems are established: inertial coordinate system{E}, and body coordinate system{B} as shown in Figure 1.

FIGURE 1

Figure 1 Inertial coordinate system and body coordinate system.

1) inertial coordinate system{E},(O_E,X_E,Y_E,Z_E). It is fixed to the Earth and follows the right-hand rule, forming a right-handed, right-angle coordinate system in the order of “North, West and Sky”, i.e. the axis X_E points to the North, the axis Y_E points to the West, and the axis Z_E is perpendicular to the ground and back from the center of the Earth. O_E is the origin.

2) body coordinate system{B},(O_B,X_B,Y_B,Z_B). It is a right-angle coordinate system solidly connected to UAV, following right-hand rule. O_B is the origin, located at the center of UAV gravity. The axis X_B in symmetry plane of UAV, parallels to UAV design axis and directs to the nose. The axis Y_B is perpendicular to symmetry plane of UAV and points to the left of the fuselage. The axis Z_B is in symmetry plane of UAV, perpendicular to plane X_B-Y_B and points above the fuselage.

The Euler angle $\text{Θ}={[\begin{array}{ccc}\varphi & \text{θ}& \text{ψ}\end{array}]}^{\text{T}}$ is used to represent the attitude angle vector of UAV in the inertial system{E}, which is defined as follows.

ϕ : Roll angle means rotation angle of UAV body around axis X_B, with positive rotation to the right.

θ : Pitch angle means rotation angle of UAV body around axis Y_B, with positive rotation upward.

ψ : Yaw angle refers to rotation angle of UAV body around axis Z_B, with positive rotation to the right.

In the beginning, the above two coordinate systems coincide. Each coordinate vector of {B} can be converted to corresponding coordinate vector of {E} by coordinate transformation and vice versa. Rotation matrix is described by the parameters in Euler angles $\text{Θ}={[\begin{array}{ccc}\varphi & \text{θ}& \text{ψ}\end{array}]}^{\text{T}}$. The specific derivation process is as follows.

When rotating around Z-axis in inertial coordinate system, as shown in Figure 2A, the transformation equation is

FIGURE 2

Figure 2 Coordinate transformation. (A) Rotating around Z-axis. (B) Rotating around Y-axis. (C) Rotating around X-axis.

$\begin{array}{ll}\{\begin{array}{c}x=\cos (\text{ψ})x+\sin (\text{ψ})y\\ y=-\sin (\text{ψ})x+\cos (\text{ψ})y\\ z=z\end{array}& (1)\end{array}$

The conversion matrix is

$\begin{array}{ll}\text{R}(\text{z},\text{ψ})=\left[\begin{array}{ccc}\cos \text{ψ}& \sin \text{ψ}& 0\\ -\sin \text{ψ}& \cos \text{ψ}& 0\\ 0& 0& 1\end{array}\right]& (2)\end{array}$

When rotating around Y-axis, as shown in Figure 2B, the transformation equation is

$\begin{array}{ll}\{\begin{array}{c}x=\cos (\text{θ})x+\sin (\text{θ})y\\ y=y\\ z=\sin (\text{θ})x+\cos (\text{θ})z\end{array}& (3)\end{array}$

The conversion matrix is

$\begin{array}{ll}\text{R}(\text{y},\text{θ})=\left[\begin{array}{ccc}\cos \text{θ}& 0& -\sin \text{θ}\\ 0& 1& 0\\ \sin \text{θ}& 0& \cos \text{θ}\end{array}\right]& (4)\end{array}$

When rotating around X-axis, as shown in Figure 2C, the transformation equation is

$\begin{array}{ll}\{\begin{array}{c}x=x\\ y=cos(\varphi )y+\sin (\varphi )z\\ z=-\sin (\varphi )y+cos(\varphi )z\end{array}& (5)\end{array}$

The conversion matrix is

$\begin{array}{ll}\text{R}(x,\varphi )=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \varphi & \sin \varphi \\ 0& -\sin \varphi & \cos \varphi \end{array}\right]& (6)\end{array}$

Transformation from {E} to{B} can be achieved by following yaw angle, then pitch angle, and then roll angle. The transformation matrix is equal to chain-multiplication of the matrices obtained from sequential rotations, i.e.

$\begin{array}{ll}{\text{R}}_{\text{E}}^{\text{B}}=\text{R}(x,\varphi )\text{R}(y,\text{θ})\text{R}(z,\text{ψ})& (7)\end{array}$

Substituting equations (2), (4), and (6), into equation (7), rotation matrix from {E} to {B} is obtained: (8); similarly, we can get rotation matrix from {B} to {E}, which is an orthogonal array: (9)

$\begin{array}{ll}{\text{R}}_{\text{E}}^{\text{B}}=\left[\begin{array}{ccc}\cos \text{θ}\cos \text{ψ}& \cos \text{θ}\sin \text{ψ}& 0\\ \sin \varphi \sin \text{θ}\cos \text{ψ}-\cos \varphi \sin \text{ψ}& \sin \varphi \sin \text{θ}\sin \text{ψ}+\cos \varphi \cos \text{ψ}& \sin \varphi \cos \text{θ}\\ \cos \varphi \sin \text{θ}\cos \text{ψ}+\sin \text{θ}\sin \text{ψ}& \cos \varphi \sin \text{θ}\sin \text{ψ}-\sin \varphi \cos \text{ψ}& \cos \varphi \cos \text{θ}\end{array}\right]& (8)\end{array}$

$\begin{array}{ll}{\text{R}}_{\text{B}}^{\text{E}}=\left[\begin{array}{ccc}\cos \text{θ}\cos \text{ψ}& \sin \varphi \sin \text{θ}\cos \text{ψ}-\cos \varphi \sin \text{ψ}& \cos \varphi \sin \text{θ}\cos \text{ψ}+\sin \text{θ}\sin \text{ψ}\\ \cos \text{θ}\sin \text{ψ}& \sin \varphi \sin \text{θ}\sin \text{ψ}+\cos \varphi \cos \text{ψ}& \cos \varphi \sin \text{θ}\sin \text{ψ}-\sin \varphi \cos \text{ψ}\\ -\sin \text{θ}& \sin \varphi \cos \text{θ}& \cos \varphi \cos \text{θ}\end{array}\right]& (9)\end{array}$

According to ${\text{ν}}_{\text{E}}={\text{R}}_{\text{B}}^{\text{E}}{\text{V}}_{\text{B}}$, where velocity under body coordinate system ${\text{V}}_{\text{B}}=({\dot{X}}_{\text{B}},{\dot{Y}}_{\text{B}},{\dot{Z}}_{\text{B}})$, velocity under inertial coordinate system ${\text{V}}_{\text{E}}=({\dot{X}}_{\text{E}},{\dot{Y}}_{\text{E}},{\dot{Z}}_{\text{E}})$. Position state equations can be solved.

Similarly, angular velocity transformation relation from {B} to {E} can be obtained, taking X-axis as an example.

$\begin{array}{ll}{\omega }^{\text{E}}=\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]=\text{R}(x,\varphi )\text{R}(y,\text{θ})\left[\begin{array}{c}0\\ 0\\ \dot{\psi }\end{array}\right]+\text{R}(x,\varphi )\left[\begin{array}{c}0\\ \dot{\theta }\\ 0\end{array}\right]+\left[\begin{array}{c}\dot{\varphi }\\ 0\\ 0\end{array}\right]=\left[\begin{array}{ccc}1& 0& -\sin \theta \\ 0& \cos \varphi & \cos \text{θ}\sin \varphi \\ 0& -\sin \varphi & \cos \varphi \cos \text{θ}\end{array}\right]\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]& (10)\end{array}$

Transforming the formula(10) gives angular velocity transformation matrix from {E} to {B}.

$\begin{array}{ll}\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]=\left[\begin{array}{ccc}1& \tan \theta \sin \varphi & \tan \theta \cos \varphi \\ 0& \cos \varphi & -\sin \varphi \\ 0& \sin \varphi /\cos \text{θ}& \cos \varphi /\cos \text{θ}\end{array}\right]& (11)\end{array}$

The purpose of establishing UAV dynamics model is to analyze the variations of UAV trajectory under external forces and moments. The input of kinetics model is external forces and moments provided by the propeller. The output of kinetics model is velocity and angular velocity, which is input of the kinematics model. Then output of the kinematic model is position and attitude angle. The relationships are as shown in Figure 3.

FIGURE 3

Figure 3 A rigid body model of UAV flight control.

2.1 Kinetic model

According to Newton-Euler Equations, motion of a rigid body equals translation of centroid plus rotation around centroid. Translation of centroid is described by Newton’s second law, i.e. $F=m\frac{dv}{dt}$, where F is resultant external force; v is velocity. Rotation round centroid is described by Euler equation, i.e. $M=J\dot{\omega }+\omega \times J\omega$, where M is resultant external moment; J is 3×3 inertia matrix, and ω is angular velocity.

Then position dynamics model is ${\dot{\upsilon }}^E=\frac{f^B}{m}-g^E$, where the letter E indicates vectors in {E} and the letter B indicates vectors in {B}. For convenience, force f is converted to Earth coordinate system by left multiplying rotation matrix ${\text{R}}_{\text{B}}^{\text{E}}:{\dot{\upsilon }}^e={\text{R}}_{\text{B}}^{\text{E}}\frac{f^B}{m}-g^E$.

φ, θ, ψ respectively indicate roll angle of X-axis rotation, pitch angle of Y-axis rotation, and yaw angle of Z-axis rotation, i.e., the Euler angle, as shown in Figure 4.

FIGURE 4

Figure 4 Attitude angle of agricultural UAV.

Expanding the equation and substituting ${\text{R}}_{\text{B}}^{\text{E}}$ into it can be obtained as

$\begin{array}{ll}\{\begin{array}{c}{\dot{\upsilon }}_x=\frac{f}{m}(\cos \text{ψ}\sin \text{θ}\cos \varphi +\sin \text{ψ}\sin \varphi )\\ {\dot{\upsilon }}_y=\frac{f}{m}(\sin \text{ψ}\sin \text{θ}\cos \varphi -\sin \varphi \cos \text{ψ})\\ {\dot{\upsilon }}_z=\frac{f}{m}(\cos \text{θ}\cos \text{ψ})-g\end{array}& (12)\end{array}$

This yields the equations regarding combined external force and velocity——the position kinetic model. Then we talk about the equations regarding combined moment and angular velocity——the attitude kinetic model.

The angular velocity equation is that the body makes a rotational motion around the axis under the combined external moment, so that the UAV attitude Euler angle $\text{Θ}={[\begin{array}{ccc}\varphi & \text{θ}& \text{ψ}\end{array}]}^{\text{T}}$changes. Moments action that UAV is subjected to when flying mainly includes: aerodynamic effects, inertial counter-torsional moments and gyroscopic moments.

For moment analysis of UAV, the angular motion equation can be obtained according to Euler’s equation as $J{\dot{\omega }}^B+{\omega }^B\times J{\omega }^B=G_a+\tau$, where, ω^B denotes angular velocity in {B}; G_a denotes the gyroscopic moment; τ denotes the counter-torsional moment generated by propeller on the body axis, including roll moment τ_x around axis O_bx_b, pitch moment τ_y around axis O_by_b, and the yaw moment τ_z around axis O_bz_b.

Regarding ω^B , for convenience, the three components ω_xb , ω_yb , ω_zb are denoted by p, q, and r respectively, i.e.

$\begin{array}{ll}{\omega }^{\text{B}}=\left[\begin{array}{c}{\omega }_{xb}\\ {\omega }_{yb}\\ {\omega }_{zb}\end{array}\right]=\left[\begin{array}{c}p\\ q\\ r\end{array}\right]& (13)\end{array}$

Inertia matrix of UAV in the body coordinate system consists of nine components, denoted as

$\begin{array}{ll}\text{J}=\left[\begin{array}{ccc}{\text{I}}_{xx}& -{\text{I}}_{xy}& -{\text{I}}_{xz}\\ -{\text{I}}_{yx}& {\text{I}}_{\text{yy}}& -{\text{I}}_{yz}\\ -{\text{I}}_{zx}& -{\text{I}}_{zy}& {\text{I}}_{zz}\end{array}\right]& (14)\end{array}$

where I_xx , I_yy , I_zz are moments of inertia, while I_xy , I_xz , I_yz , etc. are products of inertia.

$\begin{array}{ll}\{\begin{array}{c}{\text{I}}_{xx}={\displaystyle \sum m_i(y_i^2+z_i^2)}={\displaystyle \int (y^2+z^2)dm}\\ {\text{I}}_{yy}={\displaystyle \sum m_i(z_i^2+x_i^2)}={\displaystyle \int (z^2+x^2)dm}\\ {\text{I}}_{\text{zz}}={\displaystyle \sum m_i(x_i^2+y_i^2)}={\displaystyle \int (x^2+y^2)dm}\end{array}& (15)\end{array}$

$\begin{array}{ll}\{\begin{array}{c}{\text{I}}_{xy}={\text{I}}_{yx}={\displaystyle \sum m_i{x}_i{y}_i}={\displaystyle \int (xy)dm}\\ {\text{I}}_{xz}={\text{I}}_{zx}={\displaystyle \sum m_i{x}_i{z}_i}={\displaystyle \int (xz)dm}\\ {\text{I}}_{\text{yz}}={\text{I}}_{\text{zy}}={\displaystyle \sum m_i{y}_i{z}_i}={\displaystyle \int (yz)dm}\end{array}& (16)\end{array}$

Since the UAV structure is uniformly symmetric and the origin of body coordinate system coincides with the centroid. The integral of an odd function in the symmetry region is zero, thus the inertia matrix is

$\begin{array}{ll}\text{J}=\left[\begin{array}{ccc}{\text{I}}_{xx}& 0& 0\\ 0& {\text{I}}_{\text{yy}}& 0\\ 0& 0& {\text{I}}_{\text{zz}}\end{array}\right]& (17)\end{array}$

The second term on the left side of Euler’s equation can be written as

$\begin{array}{ll}{\omega }^B+J{\omega }^B=\left[\begin{array}{ccc}i& j& k\\ p& q& r\\ {\text{I}}_{xx}\text{p}& {\text{I}}_{\text{yy}}\text{q}& {\text{I}}_{\text{zz}}\text{r}\end{array}\right]=\left[\begin{array}{c}qr({\text{I}}_{\text{zz}}-{\text{I}}_{\text{yy}})\\ pr({\text{I}}_{\text{xx}}-{\text{I}}_{\text{zz}})\\ pq({\text{I}}_{\text{xx}}-{\text{I}}_{\text{yy}})\end{array}\right]& (18)\end{array}$

Gyroscopic torque G_a : When a motor rotates at high speed, it is equivalent to a gyroscope. A gyroscope spinning at high speed is a very stable individual with the ability to keep its axial direction unchanged. Therefore if an external force tries to change the direction of the gyroscopic axis, a gyroscopic moment will be generated to resist this change.

The equation for the gyroscopic moment is given by G_a=J_RPΩ×ϖ , where Ω is the propeller angular velocity vector, ϖ is the rotational angular velocity of motor rotor shaft, and we consider the rotor shaft angular velocity to be equal to the body angular velocity. Gyroscopic moment of the rotor assembly k (with motor and propeller) is G_a,k=J_RPΩ_kb₃×ϖ , where b₃ is direction of the propeller angular velocity. Because the propeller only has angular velocity in Z-axis direction in {B}, $b_3={[\begin{array}{ccc}0& 0& 1\end{array}]}^T$. J_RP is total rotational inertia of the motor rotor and the propeller. Ω_k is the kth propeller angular velocity, according to the right-hand rule, the clockwise rotation around z-axis is positive. When k=1, the propeller speed is negative, Ω_k=−ϖ₁ , the speed of the remaining propellers and so on.

$\begin{array}{ll}{b}_3\times \varpi =-(\varpi \times )b_3=\left[\begin{array}{ccc}0& -r& q\\ r& 0& -p\\ -q& p& 0\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]=\left[\begin{array}{c}q\\ -p\\ 0\end{array}\right]& (19)\end{array}$

“× “ denotes the cross operator.

Formula for G_a can be obtained as

$\begin{array}{ll}{G}_a=\left[\begin{array}{c}{G}_{a,\varphi }\\ G_{a,\theta }\\ G_{a,\psi }\end{array}\right]=-{\text{J}}_{\text{RP}}{\omega }_k\left[\begin{array}{c}q\\ -p\\ 0\end{array}\right]=\left[\begin{array}{c}-q{\text{J}}_{\text{RP}}({\varpi }_1-{\varpi }_2+{\varpi }_3-{\varpi }_4)\\ p{\text{J}}_{\text{RP}}(-{\varpi }_1+{\varpi }_2-{\varpi }_3+{\varpi }_4)\\ 0\end{array}\right]& (19)\end{array}$

Where ϖ₁,ϖ₂,ϖ₃,ϖ₄ indicates speed of propellers 1, 2, 3, and 4.

Organizing the above equations, we have

$\begin{array}{ll}\left[\begin{array}{c}{I}_{xx}\dot{p}\\ I_{yy}\dot{q}\\ I_{zz}\dot{r}\end{array}\right]=\left[\begin{array}{c}qr({\text{I}}_{\text{yy}}-{\text{I}}_{\text{zz}})\\ pr({\text{I}}_{\text{zz}}-{\text{I}}_{\text{xx}})\\ pq({\text{I}}_{\text{xx}}-{\text{I}}_{\text{yy}})\end{array}\right]=\left[\begin{array}{c}-q{\text{J}}_{\text{RP}}({\varpi }_1-{\varpi }_2+{\varpi }_3-{\varpi }_4)\\ p{\text{J}}_{\text{RP}}(-{\varpi }_1+{\varpi }_2-{\varpi }_3+{\varpi }_4)\\ 0\end{array}\right]+\left[\begin{array}{c}{\tau }_x\\ {\tau }_y\\ {\tau }_z\end{array}\right]& (20)\end{array}$

2.2 Kinematic model

The position equation for UAV is ${\dot{p}}^E=v^E$, where ${\dot{p}}^E={[\begin{array}{ccc}x& y& z\end{array}]}^T$ is coordinate position in {E}, and the expansion yields: ${[\begin{array}{ccc}\dot{x}& \dot{y}& \dot{z}\end{array}]}^T={[\begin{array}{ccc}{v}_x& v_y& v_z\end{array}]}^T$.

Then we use Euler’s angle to express the attitude. The change rate of attitude angle is related to body rotational angular velocity as follows.

$\begin{array}{ll}\dot{\Theta }=\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]=\text{W}·{\omega }^{\text{B}}=\left[\begin{array}{ccc}1& \tan \theta \sin \varphi & \tan \theta \cos \varphi \\ 0& \cos \varphi & -\sin \varphi \\ 0& \sin \varphi /\cos \text{θ}& \cos \varphi /\cos \text{θ}\end{array}\right]\left[\begin{array}{c}p\\ q\\ r\end{array}\right]& (21)\end{array}$

Where $\text{Θ}={[\begin{array}{ccc}\dot{\varphi }& \dot{\theta }& \dot{\psi }\end{array}]}^{\text{T}}$, ${\omega }^{\text{B}}={[\begin{array}{ccc}{\omega }_{xb}& {\omega }_{yb}& {\omega }_{zb}\end{array}]}^{\text{T}}={[\begin{array}{ccc}p& q& r\end{array}]}^{\text{T}}$When small perturbations happen, i.e., provided that variation of each angle is small, and change rate of attitude angle is approximately equal to the body rotational angular velocity, then we have

$\begin{array}{ll}\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]=\left[\begin{array}{c}p\\ q\\ r\end{array}\right]& (22)\end{array}$

The rigid-body flight control model of UAV consists of the kinetic model and the kinematic model. According to the above equations, there are

$\begin{array}{ll}\{\begin{array}{c}{\ddot{x}}_x=\frac{f(c\varphi s\theta c\psi +s\varphi s\psi )}{m}\\ \ddot{y}=\frac{c\varphi s\theta c\psi -s\varphi s\psi }{m}{\text{U}}_1\\ \ddot{z}=\frac{c\varphi c\theta }{m}{\text{U}}_1-g\end{array}& (23)\end{array}$

$\begin{array}{ll}\{\begin{array}{c}\ddot{\varphi }=\frac{I{\text{U}}_2+\dot{\theta }\dot{\psi }({\text{I}}_{\text{yy}}-{\text{I}}_{\text{zz}})}{{\text{I}}_{\text{xx}}}\\ \ddot{\theta }=\frac{I{\text{U}}_3+\dot{\varphi }\dot{\psi }({\text{I}}_{\text{zz}}-{\text{I}}_{\text{xx}})}{{\text{I}}_{\text{yy}}}\\ \ddot{\psi }=\frac{I{\text{U}}_4+\dot{\varphi }\dot{\psi }({\text{I}}_{\text{zz}}-{\text{I}}_{\text{xx}})}{{\text{I}}_{\text{zz}}}\end{array}& (24)\end{array}$

The following equations can be obtained.

$\begin{array}{ll}\{\begin{array}{c}{\dot{p}}^E=v^E\\ {\dot{v}}^E=g{\text{e}}_3-\frac{f}{m}R{\text{e}}_3\\ \dot{\Theta }=W·{\omega }^{\text{B}}\\ J·{\omega }^{\text{B}}=-{\omega }^{\text{B}}\times (J·{\omega }^{\text{B}})+{\text{G}}_{\text{a}}+\tau \end{array}& (25)\end{array}$

3 Controller design

3.1 Attitude angle control

The robust T-S fuzzy control method is proposed in attitude angle control part, which demands a linear expression of the nonlinear model. Therefore, we need to build the robust T-S fuzzy control model. A robust T-S fuzzy control system is a linear time-varying system whose matrix depends on a time-varying parameter vector that varies within a known interval.

In order to represent the dynamics of UAV with a set of local linear models, a nonlinear representation of the sector based on the weight and affiliation function is applied to attitude angle. Assuming that the aerodynamic effects including blade grinding, rotor inertia or drag force, are negligible, the attitude angle can be rewritten as a multiple linearized model based on the operating points. The operating points are mainly determined by the roll rate, pitch rate and yaw rate (i.e., the first order derivatives of the corresponding angles), i.e. $\text{P}=[\begin{array}{ccc}\dot{\varphi }& \dot{\theta }& \dot{\psi }\end{array}]$. The dynamics linearized model of the attitude angle is

$\begin{array}{ll}\dot{X}={\text{A}}_{\text{i}}{\text{X}}_{\text{E}}+\text{BU}& (26)\end{array}$

Where $X_E={[e_{\varphi } \dot{e_{\varphi }} e_{\theta } \dot{e_{\theta }} e_{\psi } \dot{e_{\Psi }} q_1 q_2 q_3]}^T$, $\text{U}={[\begin{array}{ccc}{\text{U}}_2& {\text{U}}_3& {\text{U}}_4\end{array}]}^{\text{T}}$, $\dot{q}={[\begin{array}{ccc}{\dot{q}}_1& {\dot{q}}_2& {\dot{q}}_2\end{array}]}^{\text{T}}={[\begin{array}{ccc}{\text{e}}_{\varphi }& {\text{e}}_{\text{θ}}& {\text{e}}_{\text{ψ}}\end{array}]}^{\text{T}}$, θ^ref , ψ^ref are given by inverse solution of the position control loop, ψ^ref =0;i=1,2,.,r, r denotes the number of control rates.

The concept of sector nonlinearity is applied to this T-S fuzzy model. If the jth operating point is restricted to α_max and α_min, the weighting function is expressed as follows

$\begin{array}{ll}{\omega }_1^{\text{j}}=1-{\omega }_0^{\text{j}}=1-\frac{{\text{α}}_{\text{max}}-{\text{P}}_{\text{j}}}{{\text{α}}_{\text{max}}-{\text{α}}_{\text{min}}}& (27)\end{array}$

Then an affiliation function based on the weighting function is defined as

$\begin{array}{ll}{h}_i=\prod _{j=1}^p{\omega }_{i_j}^j& (28)\end{array}$

where p is the number of operator points and i=1,······,2p , p=3, then r=2p=8. Defining h_i>0, P_i is restricted to α_max and α_min. P_i∈[−2,2],rad/s . According to the above equation, we can obtain ${\sum }_{\text{i}=1}^{\text{r}}{h}_i=1$. Using this affiliation function, the linearized model of attitude angle dynamics can be expressed as

$\begin{array}{ll}\dot{X}={\sum }_{\text{i}=1}^{\text{r}}{h}_i{K}_i{X}_E+\text{BU}& (29)\end{array}$

where K_i is the controller gain calculated by the robust T-S fuzzy system, i.e., every possible upper and lower bound arrangement of the components in P.

The eight rules of the robust T-S fuzzy system are shown in Table 1.

TABLE 1

Table 1 T-S control rate.

To design robust T-S fuzzy system controller gain, the robust T-S fuzzy system optimal control problem is expressed as an optimization problem with LMI constraints.

Given the LMI parameters Q=QT≥0,R=RT>0, the optimal performance bound μ>0 (such that J(x, u)<μ), and the decay rate η>0. Then, controller gain K_i calculated by the robust T-S fuzzy system is obtained by solving P∈S^9×9 , W_i∈S^3×9 , such that μ is minimized while meeting the following LMIs:

$\begin{array}{ll}\{\begin{array}{c}({\text{A}}_{\text{i}}P+B{\text{W}}_{\text{i}})+{({\text{A}}_{\text{i}}\text{P}+{\text{BW}}_{\text{i}})}^{\text{T}}+2\eta P\le 0\\ trace({\text{Q}}^{\frac{1}{2}}P{({\text{Q}}^{\frac{1}{2}})}^{\text{T}})+trace(Y)\le \mu \\ \left|\begin{array}{cc}-\text{Y}& {\text{R}}^{\frac{1}{2}}{\text{W}}_{\text{i}}\\ {({\text{R}}^{\frac{1}{2}}{\text{W}}_{\text{i}})}^{\text{T}}& -\text{P}\end{array}\right|\le 0\end{array}& (30)\end{array}$

Then K_i=−W_iP⁻¹ , it guarantees global stability while providing the required transient behavior.

Where

$$\text{R}=\left[\begin{array}{ccc}2& 0& 0\\ 0& 0.5& 0\\ 0& 0& 0.5\end{array}\right],\text{Q}={\text{I}}_9.$$

Based on the above LMI pole configuration method and T-S fuzzy system, the attitude angle robust T-S fuzzy controller is designed, and we get the attitude angle control law as

$\begin{array}{ll}\text{U}=\left[\begin{array}{c}{\text{U}}_2\\ {\text{U}}_3\\ {\text{U}}_4\end{array}\right]={\sum }_{\text{i}=1}^{\text{r}}{h}_i{K}_i{X}_E& (31)\end{array}$

3.2 Position loop control

For the position loop control part, a PID controller is used, with the error defined as

$\begin{array}{ll}{e}_i=X_d=X_i-X_i^{ref}& (32)\end{array}$

Where $X=[\begin{array}{ccc}x& y& z\end{array}]$, denotes the real-time coordinates and attitude angle values of UAV obtained by the sensors, and

$X^{ref}=[\begin{array}{ccc}{x}^{ref}& y^{ref}& z^{ref}\end{array}]$ denotes the desired coordinates of UAV. The desired coordinates of UAV are given by its desired trajectory. The desired attitude angle θ^ref , ψ^ref are given by the following equations.

$\begin{array}{ll}\{\begin{array}{c}{\varphi }^{ref}=arcsin(\frac{m(U_{x}sin(\Psi )-U_{y}sin(\Psi ))}{U_1})\\ {\theta }^{ref}=arcsin(\frac{m{U}_x-U_{1}sin(\varphi )cos(\Psi ))}{U_{1}cos(\varphi )cos(\Psi )})\\ U_x=cos(\varphi )sin(\theta )cos(\psi )+sin(\varphi )sin(\Psi )\\ U_y=cos(\varphi )sin(\theta )cos(\psi )-sin(\varphi )sin(\Psi )\end{array}& (33)\end{array}$

The error matrix of the position loop is defined as $e=[\begin{array}{ccc}{x}_d& y_d& z_d\end{array}]$. According to the PID control law, controller for error handling can be written as

$\begin{array}{ll}u=k_{p}e+k_i{\displaystyle \int edt+k_d\frac{de}{dt}}& (34)\end{array}$

Then the control law of position loop is $U_1=\sqrt{U_x^2+U_y^2+{(U_z-g)}^2}$, U_x , U_y , U_z is given by

$\begin{array}{ll}\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]=K_{PID}{[x_d \dot{x_d} y_d \dot{y_d} z_d {\dot{z}}_d {\displaystyle \int x_d}{\displaystyle \int y_d}{\displaystyle \int z_d}]}^{\text{T}}& (34)\end{array}$

K_PID is the control parameter.

4 Simulation verification

The simulation experiments are based on UAV parameters as shown in Table 2. Considering flight trajectory stability of UAV under load perturbation, load perturbation is simulated by mass change of UAV.

TABLE 2

Table 2 UAV parameters (Shen et al., 2021; Lotufo et al., 2020).

4.1 No change in UAV mass

In this scenario, the UAV mass is 1.121 kg. The flight trajectory of UAV is shown in Figures 5A–D. From these diagrams, we can see that UAV can track the expected flight trajectory well under the designed control algorithm, which reflects good control performance. There is no obvious deviation in the whole, and the real-time tracking of the flight trajectory is realized.

FIGURE 5

Figure 5 UAV flight trajectory with no mass change. (A) Three-dimensional trajectory of UAV, (B) UAV X-axis trajectory, (C) UAV Y-axis trajectory, (D) UAV Z-axis trajectory.

To further verify tracking performance of UAV under the designed control algorithm, angle variation diagrams are obtained. From Figures 6A–C, it can be seen that in this case, only at the start, roll angle and pitch angle have an angular change. After stabilization change of roll angle and pitch angle is less than ±0.1°. The yaw angle has been maintained at 0°. The designed control algorithm achieves the stability of UAV flight trajectory.

FIGURE 6

Figure 6 Angle variation of UAV under no-load condition. (A) Variation of roll angle, (B) Variation of pitch angle, (C) variation of yaw angle.

4.2 Small increase in UAV mass

Under this simulation scenario, UAV mass is initially 1.121kg, and changes abruptly to 1.5kg at 10s. The mass change diagram is shown in Figure 7.

FIGURE 7

Figure 7 Mass change of UAV.

From Figures 8B–D, it can be seen that UAV tracks the desired trajectory well in X and Y axis with small mass increase, consistent with constant UAV mass. While in the vertical direction, a tracking error is generated, but eventually remains smooth. The angle variation diagram is obtained from Figure 8.

FIGURE 8

Figure 8 UAV flight trajectory with small mass change. (A) Three-dimensional trajectory of UAV, (B) UAV X-axis trajectory, (C) UAV Y-axis trajectory, (D) UAV Z-axis trajectory.

It should be noted that in Figure 8D, when UAV mass changes at small scale at 10s, z-axis trajectory deviates from the desired trajectory at 10-15s. And it does not return to the desired trajectory after the deviation. Mass variation has an influence on z-axis trajectory. So the control method is affected by change of load.

From Figures 9A–C, it can be seen that when small increase of UAV mass happens, the roll angle and pitch angle have an angular change only at start-up. After stabilization the change is less than ±0.1°. The yaw angle has been maintained at 0°. The algorithm can realize control of UAV flight trajectory stability.

FIGURE 9

Figure 9 Angle variation of agricultural UAV under small load. (A) Variation of roll angle, (B) Variation of pitch angle, (C) Variation of yaw angle.

4.3 Significant increase in UAV mass

Considering a significant increase in UAV mass, the mass of UAV is also initially set to 1.121kg, which makes a significant increase to 1.5kg at 10s, indicating that a perturbation of the load occurs. The mass change diagram is shown in Figure 10.

FIGURE 10

Figure 10 Significant change of UAV mass.

The simulation gives flight trajectory as in Figure 11A. As can be seen in Figure 11B–D, when a significant increase in UAV mass happens, the UAV shows tracking errors in X and Y axis, which cannot track the desired trajectory well comparing to cases of constant UAV mass and small increase in UAV mass. In the vertical direction, the tracking error is further enlarged in comparison with the case of small UAV mass increase, and overshoot is more obvious.

FIGURE 11

Figure 11 UAV flight trajectory with significant mass change. (A) Three-dimensional trajectory of UAV, (B) UAV X-axis trajectory, (C) UAV Y-axis trajectory, (D) UAV Z-axis trajectory.

From Figure 12A–C, it can be seen that roll angle and pitch angle have an angular change situation only at start-up when UAV mass increases significantly. After stabilization, the change in roll angle and pitch angle is less than ±0.1°, and the yaw angle has been maintained at 0°. The designed controller can realize the flight trajectory stability under significant load.

FIGURE 12

Figure 12 Angle variation of UAV under significant load.

4.4 UAV mass change with time

To better reflect effect of load perturbation, UAV mass was designed to change with time. In this scenario, UAV mass is initially 1.121kg, then it increases to 3kg from 10s to 20s. And this weight is maintained until 30s, after which the mass starts to decrease and reaches the initial mass at 40s as shown in Figure 13.

FIGURE 13

Figure 13 Mass change of UAV.

This case better represents the effect of load perturbation on the stability of UAV flight trajectory, the perturbation shows dynamic changes with time, and this is closer to the real situation. The trajectory in this case is shown in Figure 14A.

FIGURE 14

Figure 14 UAV flight trajectory with mass change. (A) Three-dimensional trajectory of UAV, (B) UAV X-axis trajectory, (C) UAV Y-axis trajectory, (D) UAV Z-axis trajectory.

From Figures 14B–D, it can be seen that under the condition of time-varying UAV mass, the trajectory of the UAV in X and Y axis can track the desired trajectory well when the mass does not vary much. The tracking error appears when the mass increases more. In the vertical direction, the tracking error of the UAV varies with the UAV mass change.

From Figures 15A–C, it can be seen that roll angle and pitch angle have an angular change only at start-up. After stabilization, the change in roll angle and pitch angle is less than ±0.1°. The yaw angle has been maintained at 0°. The designed controller can realize the flight trajectory stabilization.

FIGURE 15

Figure 15 Angular variation for UAV under load variation. (A) Variation of roll angle, (B) Variation of pitch angle, (C) Variation of yaw angle.

According to the simulation results, when the UAV has a significant change in mass, i.e., a large load, the UAV will have a significant error in tracking the desired trajectory, especially in the vertical axis direction.

5 Conclusion

In the scenario of small increase in mass, UAV can track the desired trajectory well in X and Y axis, and a tracking error is generated in the vertical direction, but it eventually remains steady. In the scenario of a significant increase in mass, UAV tracked the desired trajectory tracking in X and Y axis, with a bigger tracking error in the vertical direction. In the scenario of time-varying mass, the UAV tracked the desired trajectory well in X and Y axis when the mass did not change much, and tracking errors appeared when the mass increased more. In vertical direction, the tracking error varies with the mass of the UAV.

In summary, taking agricultural UAV as research object, we adopt robust T-S fuzzy control method in attitude angle control part and PID controller in position loop control part to make UAV flight trajectory achieve stability in the simulation environment of load variation. Due to the sudden change of mass, the trajectory of UAV will deviate from the expected path, but its tracking error is small, so that the designed controller can effectively solve the trajectory tracking problem of small variable load against perturbation.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Author contributions

XW: Conceptualization, Methodology, Supervision. WX: Software, Validation. LJ: Writing - Original Draft, Data Curation. YY: Formal analysis, Writing - Review & Editing. All authors contributed to the article and approved the submitted version.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

An, Sh, Yuan, S., Li, H. D. (2018). Design of trajectory tracking controller for a quad-rotor UAV with varied payload. Electron. Optics Control 25(3), 59–63. doi: 10.3969/j.issn.1671-637X.2018.03.013

CrossRef Full Text | Google Scholar

Dong, Y., Fu, J., Yu, B., Zhang, Y., Ai, J. (2015). “Position and heading angle control of an unmanned quadrotor helicopter using LQR method”. In Chinese Control Conference, CCC (Vol. 2015-September, pp. 5566–5571) (HangZhou, China: IEEE Computer Society). doi: 10.1109/ChiCC.2015.7260508

CrossRef Full Text | Google Scholar

Emam, M., F Ak.Harian, A. (2016). “Attitude tracking of quadrotor UAV via mixed H₂/H_∞, controller: An LMI based approach,” in proceeding of 2016 24th Mediterranean Conference on Control and Automation (MED), (Athens, Greece: IEEE), 390–395. doi: 10.1109/MED.2016.7535919

CrossRef Full Text | Google Scholar

Fan, D., Xing, X., Zhao, Y., Yu, H., University, N. P. (2018). Research on UAV 4D trajectory control based on adaptive fuzzy PID. Comput. Measurement &Control 1), 112–114+119. doi: 10.16526/j.cnki:11-4762/tp.2018.01.027

CrossRef Full Text | Google Scholar

Gao, J., Duan, L. Y., Deng, L. W. (2021). Anti-interference trajectory tracking control of quadrotor UAV. Control Decision 36 (2), 8. doi: 10.13195/j.kzyjc.2019.0875

CrossRef Full Text | Google Scholar

Hu, W. Y., Xi, Y. Q., Yin, X. M. (2022). Research progress on pesticide application technology of agricultural plant protection UAV in China. Modern Agric. Sci. Technol. 10), 5. doi: 10.3969/j.issn.1007-5739.2022.10.027

CrossRef Full Text | Google Scholar

Li, F. L., Li, T. Y., Wang, Y. L. (2016). Flight attitude stability control optimization for quad rotor unmanned air vehicle. Comput. Simulation 33 (10), 43–47. doi: 10.3969/j.issn.1006-9348.2016.10.010

CrossRef Full Text | Google Scholar

Lotufo, M. A., Colangelo, L., Novara, C. (2020). Control design for UAV quadrotors via embedded model control. IEEE Trans. Control Syst. Technol. 28 (5), 1741–1756. doi: 10.1109/TCST.2019.2918750

CrossRef Full Text | Google Scholar

Luo, D. X., Jiang, J., Xu, D. W. (2020). Dual closed-loop trajectory tracking control of quadrotor load UAV. Small Special Electrical Machines 48 (9), 39–44. doi: 10.3969/j.issn.1004-7018.2020.09.011

CrossRef Full Text | Google Scholar

Qi, J. T., Ping, Y. (2018). Survey on flight control technology for hanging load UAV. Unmanned Syst. Technol. 1(1), 83–90. doi: 10.19942/j.issn.2096-5915.2018.01.008

CrossRef Full Text | Google Scholar

Ma, D. L., Xia, Y. Q., Shen, G. H., Jiang, H. R., Hao, C. X. (2021). Practical fixed-time disturbance rejection control for quadrotor attitude tracking. IEEE Trans. Ind. Electron. 68 (8), 7274–7283. doi: 10.1109/TIE.2020.3001800

CrossRef Full Text | Google Scholar

Shi, C., Lin, D. (2018). Adaptive integral backstepping control of quadrotor. Appl. Res. Comput. 35 (11), 3338–3342. doi: 10.3969/j.issn.1001-3695.2018.11.033

CrossRef Full Text | Google Scholar

Wang, H., Zhou, L. (2021). Design of a backstepping integral adaptive controller for quadrotor UAV. Acta Armamentarii 42 (6), 1283–1289. doi: 10.3969/j.issn.1000-1093.2021.06.019

CrossRef Full Text | Google Scholar

Xian, B., Zhang, S. J., Han, X. W., Cai, J. M., Wang, L. (2021). Trajectory planning for unmanned aerial vehicle slung-payload aerial transportation system based on reinforcement learning. J. Jilin University:Engineering Technol. Edition 51 (6), 2259–2267. doi: 10.13229/j.cnki.jdxbgxb20200577

CrossRef Full Text | Google Scholar

Zheng, H., Chen, J. (2021). Optimization analysis of intelligent control system for spraying operation of agricultural plant protection machine. J. Agric. Mechanization Res. 43 (12), 203–206. doi: 10.3969/j.issn.1003-188X.2021.12.036

CrossRef Full Text | Google Scholar

Zhou, X., Liu, R., Zhang, J., Zhang, X. “Stabilization of a Quadrotor With Uncertain Suspended Load Using Sliding Mode Control.” Proceedings of the ASME 2016 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Volume 5A: 40th Mechanisms and Robotics Conference. Charlotte, North Carolina, USA: IEEE ASME. doi: 10.1115/DETC2016-60060

Google Scholar