Research and Verification of Trajectory Tracking Control of a Quadrotor Carrying a Load

Abstract

This paper assumes that considering the unknown and time-varying nature of different strong and weak wind field disturbances and considering the nonlinear, under-driven, strongly coupled quadrotor carrying, a load is disturbed by the complex and variable wind field and unmodeled parts when flying in the real external environment, which will reduce the control effect of the nonlinear controller and make the vehicle fail to affect the flight effect. In order to ensure that the quadrotor carrying a load can carry supplies in the harsh environment for stable trajectory tracking, a neural network adaptive control algorithm is introduced in the article. The neural network algorithm has the role of online dynamic approximation, the compensation of arbitrary external disturbance and the compensation of external disturbance. Its structure is simple and low computation. In the article, the Lyapunov method is used to design the adaptive weight and estimate the weight of the online neural network, and the stability of the system is proved. Finally, the comparison of three algorithms verified by simulation proves that the above interference problem can be solved effectively by the proposed algorithm.

1. Introduction

The quadrotor is the main type of multi-rotor UAV and is the most used type. Because it can perform rapid vertical takeoff and landing and free hovering in the air and other maneuvers that cannot be performed by ordinary aircraft, it has an increasingly wide range of applications. Initially, it was used in the military field, such as controlling aircraft to perform various tasks in dangerous environments, ensuring the safety of human life and property and even participating in rescue operations [1,2,3]. The use of quadrotors for transporting supplies in rescue activities is one of the important directions for quadrotor application research. Currently, most solutions involve transporting loads by using ropes to hang them. Not only does the quadrotor maintain its inherent flexibility during flight, but the quadrotor can hover stationary at a designated target location to facilitate the loading and unloading of loads [4,5,6,7]. Therefore, research on the anti-disturbance control of a quadrotor carrying a load has important theoretical significance and practical applications.

A quadrotor carrying a load is a typical underdriven system with four control input signals and eight degrees of freedom, including two swing angles. The control problem is more complex compared to the flight control of a quadrotor. Since a quadrotor carrying a load cannot be controlled directly by the drive, the motion of the load is driven by controlling the motion of the quadrotor, so the motion of the load depends on the quadrotor. At the same time, there are uncertain external disturbances, such as wind fields, during the hanging and transportation of the quadrotor. This problem greatly increases the difficulty of designing quadrotor stability control for hanging loads under external perturbations [8,9].

Some excellent teams at home and abroad have conducted much research on a quadrotor carrying a load and have achieved good results. A dynamic programming method for a quadrotor carrying a load has been proposed in the literature [10]. The method ensures that the load can follow the quadrotor through narrow boundaries without touching them. However, the literature does not consider the environmental interference of various complex interferences, such as the wind field magnetic field. At the same time, the swing angle of the system is not controlled, and the coupling degree is increased. The literature [11] presents a path planning algorithm for quadrotors that models the payload delivered by a hanging rope as a hybrid dynamical system. The experimental results show that this approach is effective in generating an aggressive hybrid and avoidance maneuver trajectory, including state transitions. However, the literature does not consider the suppression of the swing angle of the rope. In [12], a nonlinear adaptive integral backstepping control method is proposed to reduce the impact of load swing on the flight, aiming at the problem that the stability of quadrotor carrying a load system is easily affected by load swing caused by wind disturbance. However, the research on the influence of wind disturbance is relatively superficial, and the data are few. In [13], the authors propose a nonlinear stable quadrotor and its load controller. However, the literature only simulates the flight simulation of the quadrotor and its hanging load on a two-dimensional plane. In the actual three-dimensional space, the coupling degree of the nonlinear term of the system is deepened, and the design method suitable for a two-dimensional plane cannot be simply extended to a three-dimensional one. Therefore, the extension of the control system to a three-dimensional space while preserving the nonlinearity of the system has become the current research direction. In [14], the authors study the control algorithm of a quadrotor which uses rope to carry the load. By using Kane’s matrix method, the literature develops a nonlinear controller that stabilizes the system through the linear part of the quadrotor. Similarly, the literature does not consider the wind disturbances usually encountered during quadrotor flight. In [15], the authors use a downward camera to observe the load and use an extended Kalman filter method to estimate its state relative to the quadrotor for the first time to realize the flexible maneuverability of the lifting load and the closed-loop control of the maximum load angle. Regrettably, the literature does not verify the effectiveness of the algorithm under wind field disturbance. In [16], the authors linearize and decouple the quadrotor model at the hovering point to generate four subsystem models, which are linear. An improved adaptive L1 control method and a self-disturbance suppression linear controller are proposed. This method is very effective for interference suppression, but the suppression effect for the swing angle is not good enough. At the same time, it is too simplified for the model of the quadrotor carrying a load. In [17], two ways of transporting load with quadrotor, namely, grab and rope hanging, are studied, and different control methods and maneuvering strategies are analyzed. From the literature, there are two advantages of rope hanging. The first is that there is no requirement for the shape of the load, and there is no mismatch between the quadrotor and the load. The second is that the quadrotor can hover at a fixed point to facilitate the loading and unloading of the load. In [18], the authors use a new type of passive flexible grasping structure and servo positioning control, combined with detection algorithms, which can make the micro-air vehicle grasp and transport loads with different shapes and weights. However, the load cannot be too heavy during crawling, and the adjustment parameters of the algorithm are too complicated. In [19], the quadrotor robot transports the quadrotor robot with the load at an acceleration exceeding the ability of the previous collaborative method, but there is no clear communication between the collaborative robots. They rely entirely on visual cues and inertia obtained from airborne sensors. So, the airborne sensor will be interfered with and lose its function in the case of the mixed and complex interference of the wind field and the magnetic field. In [20], the authors show that the PID-controlled helicopter refuses to increase the load within the allowable load range and adjusts the period to simulate dominance. This means the load is only allowed to be placed once, and PID tuning is relatively slow.

The control of the trajectory tracking of a quadrotor carrying a load has been widely and deeply researched both at home and abroad. However, considering the uncertainty of the model parameters, the unknown time-varying wind field magnetic field perturbation under the actual working conditions and the need to combine specific engineering practices to achieve load swing suppression and accurate trajectory tracking, the research on quadrotors carrying load suspensions combining all the above points is relatively less. Therefore, on the basis of summarizing the research results of the predecessors on the quadrotor carrying a load [21,22,23,24,25,26,27,28,29] and the algorithms involved in the article [30,31], combined with our team’s non-linear control of the quadrotor [12,32,33], a neural network control algorithm based on adaptive integral backstepping trajectory tracking is proposed to solve the complex and variable wind field and the disturbance problem caused by unmodeled dynamics. The contributions of this paper are the following: In Section 1, the current status of research on a quadrotor carrying a load and its algorithm is described, and the previous research results and their shortcomings are summarized, which leads to the problem addressed in this paper. In Section 2, the quadrotor and its load are treated as a whole system for force analysis, and the mathematical model of the coupled quadrotor carrying a load with eight degrees of freedom is established. In Section 3, a neural network adaptive integral backstepping nonlinear controller is designed to estimate and compensate for the unknown external disturbances and uncertainties of the model parameters in response to the problem of complex external wind field disturbances. In this part, the roll angle in the attitude angle is analyzed in detail, and its stability is proved by the Lyapunov equation. In Section 4, the proposed algorithm is combined with the adaptive integral backstepping algorithm, the backstepping algorithm is investigated by the team for tracking control simulation to verify the spiral trajectory in 3D space and the superiority and efficiency of the method are demonstrated.

2. Establishment of a Quadrotor Carrying a Load

Building a suitable and complete model of a quadrotor carrying a load is the basis of the control algorithm. The quadrotor is a standard cross-shaped quadrotor with a general structure. The body coordinate system of the quadrotor is defined as $g=\left(x,y,z\right)$, the inertial coordinate system is defined as $b=\left(X,Y,Z\right)$ and $I={\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}^{\mathrm{T}}$ is the spatial position coordinate of the inertial coordinate system of the quadrotor. The front and rear direction of the quadrotor is the X-axis direction, and the backward direction is negative. The left and right direction is the Y-axis direction, and the right direction is positive. The up and down direction is the Z-axis direction, and the upward direction is positive. The moving speed of the quadrotor is expressed as $V={\left[\begin{array}{ccc}u& v& w\end{array}\right]}^{\mathrm{T}}$. $P={\left[\begin{array}{ccc}\varphi & \theta & \psi \end{array}\right]}^{\mathrm{T}}$ is the Euler angle of the space attitude. The roll angle $\varphi$ around the X-axis rolls to the right in the positive direction. At this time, the rotation speed of rotor 2 increases, and the rotation speed of rotor 4 decreases. When the nose is down, the pitch angle $\theta$ around the Y-axis is in the positive direction. In this case, the speed of rotor 1 increases, and the speed of rotor 3 decreases. The yaw angle $\psi$ around the Z-axis rotates counterclockwise in the positive direction. At this time, the speeds of rotor 1 and rotor 3 remain unchanged, and the speeds of rotor 2 and rotor 4 increase. The rotational angular velocity is $\omega ={\left[\begin{array}{ccc}p& q& r\end{array}\right]}^{\mathrm{T}}$. $\alpha$ is the projection angle of the hanging rope on the $xoz$ surface of the quadrotor, $\beta$ is the projection angle of the hanging rope on the $yoz$ surface of the quadrotor, L is the rope length, l is the half-length of the motor shaft, the weight of the quadrotor is M and the weight of the hanging load is m, so the resulting model of a quadrotor carrying a load is shown in Figure 1.

To simplify the modeling of a quadrotor carrying a load, the model should be simplified to some extent, so the following assumptions are made:

• The weight and moment of inertia of the quadrotor remain constant.
• The hanging rope is rigid and independent of the mass.
• The load always hangs from the center of the mass of the quadrotor.
• The payload is considered as a point-mass.
• The hanging load always hangs below the quadrotor, i.e., the swing angle of a quadrotor carrying a load is always between $\left(-\frac{\pi }{2},\frac{\pi }{2}\right)$.

Before modeling a quadrotor carrying a load, a mathematical model of a conventional quadrotor must be obtained first.

The conversion matrix between the airframe coordinate system and the geodesic coordinate system is given by the following formula:

$$C_b^g=\left[\begin{array}{ccc}cos\theta cos\psi & sin\varphi sin\theta cos\psi -cos\varphi sin\psi & cos\varphi sin\theta cos\psi +sin\varphi sin\psi \\ cos\theta sin\psi & sin\varphi sin\theta sin\psi +cos\varphi cos\psi & cos\varphi sin\theta sin\psi -sin\varphi cos\psi \\ -sin\theta & sin\varphi cos\theta & cos\varphi cos\theta \end{array}\right]$$

(1)

The total combined force U of the quadrotor is expressed in the inertial coordinate system as follows:

$$U=C_b^g{U}_1-\left[\begin{array}{c}0\\ 0\\ Mg\end{array}\right]=C_b^g\left[\begin{array}{c}0\\ 0\\ U_1\end{array}\right]-\left[\begin{array}{c}0\\ 0\\ Mg\end{array}\right]=\left[\begin{array}{c}{U}_1\left(cos\varphi sin\theta cos\psi +sin\varphi sin\psi \right)\\ U_1\left(cos\varphi sin\theta sin\psi -sin\varphi cos\psi \right)\\ U_1\left(cos\varphi cos\theta \right)-Mg\end{array}\right]$$

(2)

where $U_1$ represents the total lift generated by the four rotors, and g represents gravitational acceleration. After completion, the dynamic equation of the quadrotor’s inertial coordinate system is expressed as Equation (3):

$$\left\{\begin{array}{l}\ddot{\varphi }=\left[l{U}_2+\dot{\theta }\dot{\psi }\left(I_y-I_z\right)\right]/I_x+d_{\varphi }\\ \ddot{\theta }=\left[l{U}_3+\dot{\varphi }\dot{\psi }\left(I_z-I_x\right)\right]/I_y+d_{\theta }\\ \ddot{\psi }=\left[U_4+\dot{\varphi }\dot{\theta }\left(I_x-I_y\right)\right]/I_z+d_{\phi }\\ \ddot{X}=U_{1x}/M+d_x\\ \ddot{Y}=U_{1y}/M+d_y\\ \ddot{Z}=U_{1z}/M-g+d_z\end{array}\right.$$

(3)

where $U_{1x}$, $U_{1y}$ and $U_{1z}$ are denoted as follows:

$$\left\{\begin{array}{l}{U}_{1x}=\left(cos\varphi sin\theta cos\psi +sin\varphi sin\psi \right)U_1\\ U_{1y}=\left(cos\varphi sin\theta sin\psi -sin\varphi cos\psi \right)U_1\\ U_{1z}=\left(cos\varphi cos\theta \right)U_1\end{array}\right.$$

(4)

where $U_2$, $U_3$ and $U_4$ represent the roll force, pitch force and yaw moment. $I_x$, $I_y$ and $I_z$ represent the moment of inertia on the corresponding coordinate axes, and l is the half-length of the motor shaft; $d_i$ is the total disturbance term of model uncertainty, unknown external disturbance, air damping torque, propeller gyro torque and wind field disturbance.

Based on the quadrotor model established above, a hanging load is added for modeling.

The coordinates of the hanging loads in the inertial coordinate system are shown in Equation (5):

$$\left\{\begin{array}{l}{x}_1=x+Lsin\beta \\ y_1=y+Lsin\alpha cos\beta \\ z_1=z-Lcos\beta cos\alpha \end{array}\right.$$

(5)

where L is the length of the hanging rope, x, y, z are the position coordinates of the quadrotor and $x_1$, $y_1$, $z_1$ are the position coordinates of the hanging load.

Using Hamilton’s principle and Lagrange’s formula, the system is simplified as follows:

$$\frac{d}{d_{\tau }}\left(\frac{\partial \left(T-P\right)}{\partial \dot{n}}\right)-\left(\frac{\partial \left(T-P\right)}{\partial n}\right)=0$$

(6)

where T is the kinetic energy of the whole system, P is the potential energy of the whole system and n represents the $\alpha$ swing angle or the $\beta$ swing angle. Calculate T and P separately.

The first-order derivation of both sides of Equation (5) leads to Equation (7):

$$\left\{\begin{array}{l}{\dot{x}}_1=\dot{x}+Lcos\beta \dot{\beta }\\ {\dot{y}}_1=\dot{y}-L\left(sin\beta sin\alpha \dot{\beta }-cos\beta cos\alpha \dot{\alpha }\right)\\ {\dot{z}}_1=\dot{z}+L\left(sin\beta cos\alpha \dot{\beta }+cos\beta sin\alpha \dot{\alpha }\right)\end{array}\right.$$

(7)

The kinetic energy T of the entire system is shown in Equation (8):

$$T=\frac{1}{2}M\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)+\frac{1}{2}m\left({\dot{x}}_1^2+{\dot{y}}_1^2+{\dot{z}}_1^2\right)$$

(8)

where m represents the mass of the hanging load, which can be found from Equation (7) as T:

$$\begin{array}{l}T=\frac{1}{2}M\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)+\frac{1}{2}m\left[\begin{array}{l}{\left(\dot{x}+Lcos\beta \dot{\beta }\right)}^2\\ +{\left(\dot{y}-L\left(sin\beta sin\alpha \dot{\beta }-cos\beta cos\alpha \dot{\alpha }\right)\right)}^2\\ +{\left(\dot{z}+L\left(sin\beta cos\alpha \dot{\beta }+cos\beta sin\alpha \dot{\alpha }\right)\right)}^2\end{array}\right]\\ =\frac{1}{2}\left(M+m\right)\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)+\frac{1}{2}m{L}^2{\dot{\beta }}^2+\frac{1}{2}m{L}^2{cos}^2\beta {\dot{\alpha }}^2\\ +mLcos\beta \dot{\beta }\dot{x}+mL\left(cos\beta cos\alpha \dot{\alpha }-sin\beta sin\alpha \dot{\beta }\right)\dot{y}+mL\left(cos\beta sin\alpha \dot{\alpha }+sin\beta cos\alpha \dot{\beta }\right)\dot{z}\end{array}$$

(9)

When the potential energy P of the system is zero potential energy on the ground plane, the formula P can be written as follows:

$$P=Mgz+mg\left(z-Lcos\alpha cos\beta \right)$$

(10)

Substituting Equations (9) and (10) into Equation (6), when n takes $n=\alpha$ and $n=\beta$ and plugs into Hamilton’s formula, the relationship between the swing motion of the hangers and the quadrotor acceleration can be yielded as Equation (11):

$$\left\{\begin{array}{l}L\ddot{\beta }+cos\beta \ddot{x}-sin\beta sin\alpha \ddot{y}+sin\beta cos\alpha \ddot{z}+Lsin\beta cos\beta {\dot{\alpha }}^2+gsin\beta cos\alpha =0\\ cos\beta \ddot{\alpha }+cos\alpha \ddot{y}+sin\alpha \ddot{z}-2Lsin\beta \dot{\beta }\dot{\alpha }+gsin\alpha =0\end{array}\right.$$

(11)

Taking the second-order derivative of Equation (5), the acceleration of the load can be expressed in terms of the acceleration of the quadrotor, as shown in the following formula:

$$\left\{\begin{array}{l}{\ddot{x}}_1=\ddot{x}-L{\dot{\beta }}^{2}sin\beta +L\ddot{\beta }cos\beta \\ {\ddot{y}}_1=\ddot{y}-L{\dot{\alpha }}^{2}sin\alpha cos\beta -L{\dot{\beta }}^{2}sin\alpha cos\beta +L\ddot{\alpha }cos\alpha cos\beta -L\ddot{\beta }sin\alpha sin\beta -2L\dot{\alpha }\dot{\beta }cos\alpha sin\beta \\ {\ddot{z}}_1=\ddot{z}+L{\dot{\alpha }}^{2}cos\alpha cos\beta +L{\dot{\beta }}^{2}cos\alpha cos\beta +L\ddot{\alpha }sin\alpha cos\beta +L\ddot{\beta }cos\alpha sin\beta -2L\dot{\alpha }\dot{\beta }sin\alpha sin\beta \end{array}\right.$$

(12)

A force analysis was performed on the quadrotor. The tension F of the rope on the quadrotor is positive and downward:

$$\left\{\begin{array}{l}M\ddot{x}=U_{1x}+Fsin\beta \\ M\ddot{y}=U_{1y}+Fsin\alpha cos\beta \\ M\ddot{z}=U_{1z}-Fcos\alpha cos\beta -Mg\end{array}\right.$$

(13)

$$\left\{\begin{array}{l}m\ddot{x}=-Fsin\beta \\ m\ddot{y}=-Fsin\alpha cos\beta \\ m\ddot{z}=Fcos\alpha cos\beta -mg\end{array}\right.$$

(14)

Combining Equations (13) and (14), the following results are yielded:

$$\left\{\begin{array}{l}M\ddot{x}+m{\ddot{x}}_1=U_{1x}\\ M\ddot{y}+m{\ddot{y}}_1=U_{1y}\\ M\ddot{z}+m{\ddot{z}}_1=U_{1z}-\left(M+m\right)g\end{array}\right.$$

(15)

Substituting Equations (2) and (12) into the above Equation (15), the dynamics model of a quadrotor carrying a load can be yielded:

$$\left\{\begin{array}{l}\left(M+m\right)\ddot{x}-mL{\dot{\beta }}^{2}sin\beta +mL\ddot{\beta }cos\beta =U_{1x}\\ \left(M+m\right)\ddot{y}-mL{\dot{\alpha }}^{2}sin\alpha cos\beta -mL{\dot{\beta }}^{2}sin\alpha cos\beta +mL\ddot{\alpha }cos\alpha cos\beta \\ -mL\ddot{\beta }sin\alpha sin\beta -2mL\dot{\alpha }\dot{\beta }cos\alpha sin\beta =U_{1y}\\ \left(M+m\right)\ddot{z}+mL{\dot{\alpha }}^{2}cos\alpha cos\beta +mL{\dot{\beta }}^{2}cos\alpha cos\beta \\ +mL\ddot{\alpha }sin\alpha cos\beta +mL\ddot{\beta }cos\alpha sin\beta -2mL\dot{\alpha }\dot{\beta }sin\alpha sin\beta =U_{1z}-\left(M+m\right)g\end{array}\right.$$

(16)

By combing Equations (3), (11) and (16), the eight-degree-of-freedom mathematical model of a quadrotor carrying a load including attitude, position and swing angle can be yielded:

$$\left\{\begin{array}{l}\ddot{\varphi }=\left[l{U}_2+\dot{\theta }\dot{\psi }\left(I_y-I_z\right)\right]/I_x+d_{\varphi }\\ \ddot{\theta }=\left[l{U}_3+\dot{\varphi }\dot{\psi }\left(I_z-I_x\right)\right]/I_y+d_{\theta }\\ \ddot{\psi }=\left[U_4+\dot{\varphi }\dot{\theta }\left(I_x-I_y\right)\right]/I_z+d_{\phi }\\ \ddot{X}=\left[U_{1x}+mL\left({\dot{\beta }}^{2}sin\beta -\ddot{\beta }cos\beta \right)\right]/\left(M+m\right)+d_x\\ \ddot{Y}=\left[\begin{array}{ll}& U_{1y}+mL\ddot{\beta }sin\beta sin\alpha -mL\ddot{\alpha }cos\beta cos\alpha +mL{\dot{\beta }}^{2}cos\beta sin\alpha \\ & +mL{\dot{\alpha }}^{2}cos\beta sin\alpha +2mL\dot{\beta }\dot{\alpha }sin\beta cos\alpha \end{array}\right]/\left(M+m\right)+d_y\\ \ddot{Z}=\left[\begin{array}{l}{U}_{1z}-mL\ddot{\beta }sin\beta cos\alpha -mL\ddot{\alpha }cos\beta sin\alpha -mL{\dot{\beta }}^{2}cos\beta cos\alpha \\ -mL{\dot{\alpha }}^{2}cos\beta cos\alpha +2mL\dot{\beta }\dot{\alpha }sin\beta sin\alpha \end{array}\right]/\left(M+m\right)-g+d_z\\ \ddot{\alpha }=\left[2MLsin\beta \dot{\beta }\dot{\alpha }-cos\alpha U_{1y}-sin\alpha U_{1z}\right]/MLcos\beta +d_{\alpha }\\ \ddot{\beta }=\left[-MLsin\beta cos\beta {\dot{\alpha }}^2-cos\beta U_{1x}-sin\beta sin\alpha U_{1y}-sin\beta cos\alpha U_{1z}\right]/ML+d_{\beta }\end{array}\right.$$

(17)

According to the model Equation (17), the acceleration of the defined quadrotor is related to the swing angle of a quadrotor carrying a load, and each degree of freedom effects the external disturbances. Therefore, the designed controller should attenuate the swing angle oscillation in the environment of external time-varying disturbances while being able to stabilize quickly to track the given trajectory stably and improve the robustness of the system.

3. Controller Design

3.1. Attitude Controller Design

After modeling the eight-degree-of-freedom model and defining each quantity, it is shown in Figure 2.

In the figure, $Z_d$, $Y_d$, $X_d$ are the expected height, Y position and X position of the quadrotor, respectively. $u_{throttle}$ is the height control quantity output by the height controller to the motor, $u_x$ and $u_y$ are the output of the horizontal position controller, and combined with the output $u_{\alpha }$, $u_{\beta }$ of the hanging swing angle controller, the final horizontal direction control quantity $U_x$, $U_y$ can be obtained. Considering the coupling effect of the yaw angle and the expected roll angle, pitch angles ${\varphi }_d$ and ${\theta }_d$ are obtained after attitude calculation. Combined with the given expected yaw angle ${\phi }_d$ as the control variable of the attitude angle controller, the output to the motor’s horizontal and the control variables of the roll angle, pitch angle and yaw angle are $u_{roll}$, $u_{pitch}$, $u_{yaw}$, respectively. $u_1$, $u_2$, $u_3$, $u_4$ are obtained by allocating the four control variables of the motor, which is the PWM value of the motor, which defines the virtual control value $U_1$, $U_2$, $U_3$, $U_4$. The output of eight channels can be obtained from the system model $X/\dot{X}$, $Y/\dot{Y}$, $Z/\dot{Z}$, $\varphi /\dot{\varphi }$, $\theta /\dot{\theta }$, $\psi /\dot{\psi }$, $\alpha /\dot{\alpha }$, $\beta /\dot{\beta }$.

The relationship between the control variables of the four PWM motors and the control variables of the designed virtual channels is defined by the following Equation (18):

$$\left\{\begin{array}{l}{u}_1=u_{throttle}+u_{yaw}+u_{pitch}\\ u_2=u_{throttle}+u_{yaw}-u_{pitch}\\ u_3=u_{throttle}-u_{yaw}+u_{roll}\\ u_4=u_{throttle}-u_{yaw}-u_{roll}\end{array}\right.$$

(18)

As shown in Figure 2, $u_1$, $u_2$, $u_3$ and $u_4$ are the PWM value of the four rotor motors. $u_{throttle}$ represents the input signal of the throttle channel, and $u_{roll}$, $u_{pitch}$ and $u_{yaw}$ represent the control input values of the roll channel, pitch channel and yaw channel, respectively. The motors drive the rotor rotation of the quadrotor, and the generated lift and reaction moments are shown in Equation (19), as follows:

$$\left\{\begin{array}{l}{L}_i=M_t\frac{\omega }{s+\omega }{u}_i,i=1,2,3,4\\ Q_i=M_y\frac{\omega }{s+\omega }{u}_i,i=1,2,3,4\end{array}\right.$$

(19)

$L_i$ is the upward vertical lift generated by the rotor, $Q_i$ is the horizontal counter-torque, $M_t$ and $M_y$ are the lift and counter-torque coefficients, respectively, $\omega$ is the motor bandwidth and $u_i$ is the PWM control input signal to the motor. In calculating the relationship between the motor PWM signal and the force and the torque generated by the motor propeller, it is assumed that the motor is an ideal case with a sufficiently fast response, so the model can be simplified to the next relationship, see Equation (20):

$$\left\{\begin{array}{l}{L}_i=M_t{u}_i,i=1,2,3,4\\ Q_i=M_y{u}_i,i=1,2,3,4\end{array}\right.$$

(20)

The relationship between the virtual control quantity and the input signal is given by Equation (21):

$$\left\{\begin{array}{l}{U}_1=M_t(u_1+u_2+u_3+u_4)=M_t{\displaystyle \sum _{i=1}^4}{u}_i\\ U_2=M_t(u_3-u_4)\\ U_3=M_t(u_1-u_2)\\ U_4=M_y(u_1+u_2-u_3-u_4)\end{array}\right.$$

(21)

Assuming that a quadrotor carrying a load is in the environment of time-varying wind field disturbance, the attitude controller and swing angle controller are designed by the neural network adaptive method based on integral backstepping. The error integration term is introduced to reduce the steady-state error of the system and improves the steady-state performance of the system, and a neural network adaptive algorithm with strong adaptive and learning capabilities is incorporated to estimate online and compensate for the errors caused by model uncertainty and external disturbances.

Taking the roll angle channel of the quadrotor as an example, the roll angle part is extracted, the roll angle controller is designed and the roll angle dynamic model is organized as in Equation (22):

$$\ddot{\varphi }=U_2\frac{l}{I_x}+\frac{I_y-I_z}{I_x}\dot{\theta }\dot{\psi }+d_{\varphi }$$

(22)

Step 1: Based on the backstepping method, define the roll angle tracking error as $e_1$:

$$e_1={\varphi }_d-\varphi$$

(23)

where ${\varphi }_d$ is the expected rolling angle. The derivative of the rolling angle tracking error can be obtained by taking the derivative of Equation (23) as ${\dot{e}}_1$:

$${\dot{e}}_1={\dot{\varphi }}_d-\dot{\varphi }={\dot{\varphi }}_d-p$$

(24)

$p=\dot{\varphi }$ is the actual roll angular acceleration of the quadrotor. To stabilize $e_1$, the Lyapunov function of $e_1$ is chosen to be $V_1$:

$$V_1=\frac{1}{2}{e_1}^2\ge 0$$

(25)

The derivative of Lyapunov’s formula for Equation (25) yields Equation (26):

$${\dot{V}}_1=e_1{\dot{e}}_1=e_1\left({\dot{\varphi }}_d-p\right)$$

(26)

The actual roll angular velocity p is defined as the virtual control input, while the desired virtual control is defined as $p_{d1}$. To make ${\dot{V}}_1\le 0$, let

$$p_{d1}=k_1{e}_1+{\dot{\varphi }}_d$$

(27)

Step 2: Considering the need to eliminate the effect of model uncertainty, enhance robustness and reduce steady-state errors, an integration term ${\mu }_1{\gamma }_1$ is added to the virtual control variable:

$$p_{d1}=k_1{e}_1+{\dot{\varphi }}_d+{\mu }_1{\gamma }_1$$

(28)

where ${\gamma }_1={\int }_0^t{e}_1\left(\tau \right)d\tau$. $k_1$ and ${\mu }_1$ are constants and are greater than 0.

Step 3: Define the error of the virtual control input as $e_2$:

$$\begin{array}{l}{e}_2=p_{d1}-p=k_1{e}_1+{\dot{\varphi }}_d+{\mu }_1{\gamma }_1-p=k_1{e}_1+{\mu }_1{\gamma }_1+{\dot{e}}_1\end{array}$$

(29)

Substituting Equation (29) into Equation (24), obtain:

$$\begin{array}{l}{\dot{V}}_1=e_1{\dot{e}}_1=e_1\left(-k_1{e}_1-{\mu }_1{\gamma }_1+e_2\right)=-k_1{e_1}^2-e_1{\mu }_1{\gamma }_1+e_1{e}_2\end{array}$$

(30)

Choose a new Lyapunov function to simultaneously stabilize $e_1$, $e_2$ and define Lyapunov $V_2$ as:

$$V_2=V_1+\frac{1}{2}{e}_2^2+\frac{1}{2}{\mu }_1{\gamma }_1^2\ge 0$$

(31)

The derivative of Equation (31) yields:

$$\begin{array}{l}{\dot{V}}_2=e_1{\dot{e}}_1+e_2{\dot{e}}_2+{\mu }_1{\gamma }_1{\dot{\gamma }}_1\\ =-k_1{e_1}^2+e_1{e}_2+e_2\left[k_1{\dot{e}}_1+{\ddot{\varphi }}_d-\frac{\dot{\theta }\dot{\phi }\left(I_y-I_z\right)}{I_x}-\frac{l}{I_x}{U}_2-d_{\varphi }\right]\end{array}$$

(32)

Step 4: Since there are uncertain parameters and external wind disturbances $d_{\varphi }$ in Equation (32) and the upper limit of disturbance $d_{\varphi }$ is unknown, the approximation of the neural network control is designed to estimate and compensate for the uncertain parameters and external disturbances.

The algorithm of the neural network is as follows:

$$h_j=exp\left(∥x-{c_j∥}^2/\left(2b_j^2\right)\right)$$

(33)

$$d_{\varphi }=W_1^{*T}{h}_1\left(x\right)+{\zeta }_1$$

(34)

where x is the input to the neural network, $c_j$ is the center of the training data points, j is the number of nodes in the hidden layer and $b_j$ is a parameter representing the arrival rate or the rate at which the function drops to zero. The smaller the arrival rate, the narrower the image. $h={\left[h_j\right]}^T$ is the output of the Gaussian function of the neural network, W is the upper layer weight of the neural network and $\zeta$ is the approximation error of the neural network.

After designing the adaptive law of network weights, the RBF network input is $e_2$ and $e_2$ is the roll angular velocity tracking error. Then, the output of the network is:

$${\widehat{d}}_{\varphi }={\widehat{W}}_1^T{h}_1\left(x\right)$$

(35)

$$\begin{array}{l}{d}_{\varphi }-{\widehat{d}}_{\varphi }=W_1^{*T}{h}_1\left(x\right)+{\zeta }_1-{\widehat{W}}_1^T{h}_1\left(x\right)=-{\tilde{W}}_1^T{h}_1\left(x\right)+{\zeta }_1\end{array}$$

(36)

where $\tilde{W}=\widehat{W}-W^{*}$. Therefore, the rolling angle controller $U_2$ is designed as follows:

$$U_2=\frac{I_x}{l}\left[\begin{array}{l}\left(1-{k_1}^2+{\mu }_1\right)e_1+\left(k_1+k_2\right)e_2-k_1{\mu }_1{\gamma }_1+{\ddot{\varphi }}_d-\frac{\dot{\theta }\dot{\phi }\left(I_y-I_z\right)}{I_x}-{\widehat{d}}_{\varphi }\end{array}\right]$$

(37)

3.2. Stability Analysis

To demonstrate the stability of adaptive integral backstepping control of neural networks, a new Lyapunov function is designed, and the Lyapunov function $V_3$ is defined as:

$$V_3=V_2+\frac{1}{2\lambda }{\left(\widehat{W}-W^{*}\right)}^2\ge 0$$

(38)

The derivative of $V_3$ yields:

$$\begin{array}{l}{\dot{V}}_3={\dot{V}}_2+\frac{1}{\lambda }\left(\widehat{W}-W^{*}\right)\frac{d\widehat{W}}{dt}\\ =-k_1{e_1}^2+e_1{e}_2+e_2\left[\begin{array}{l}{k}_1{\dot{e}}_1+{\ddot{\varphi }}_d-\frac{\dot{\theta }\dot{\phi }\left(I_y-I_z\right)}{I_x}-\frac{l}{I_x}{U}_2-d_{\varphi }\end{array}\right]+\frac{1}{\lambda }\left(\widehat{W}-W^{*}\right)\frac{d\widehat{W}}{dt}\end{array}$$

(39)

Substituting Equations (36) and (37) into Equation (39), we obtain:

$$\begin{array}{l}{\dot{V}}_3=-k_1{e_1}^2-k_2{e_2}^2+e_2\left({\widehat{d}}_{\varphi }-d_{\varphi }\right)-\frac{1}{\lambda }\left(\widehat{W}-W^{*}\right)\dot{\widehat{W}}\\ =-k_1{e_1}^2-k_2{e_2}^2+e_2\left(\widehat{W}-W^{*}\right)h_1\left(x\right)-\frac{1}{\lambda }\left(\widehat{W}-W^{*}\right)\dot{\widehat{W}}\\ =-k_1{e_1}^2-k_2{e_2}^2+\left(\widehat{W}-W^{*}\right)\left(e_2{h}_1\left(x\right)-\frac{1}{\lambda }\dot{\widehat{W}}\right)\end{array}$$

(40)

At this point, the neural network is adjusted mainly by adjusting the weights of the neural network. Therefore, the neural network weight adaptive method is designed as $\dot{\widehat{W}}=\lambda e_2{h}_1\left(x\right)$. Then, we can obtain:

$${\dot{V}}_3=-k_1{e_1}^2-k_2{e_2}^2\le 0$$

(41)

Therefore, the designed roll angle control system is progressively stable.

In the same way, we can obtain the pitch angle controller and the yaw angle controller:

$$U_3=\frac{I_y}{l}\left[\begin{array}{l}\left(1-{k_3}^2+{\mu }_2\right)e_3+\left(k_3+k_4\right)e_4-k_3{\mu }_2{\gamma }_2+{\ddot{\theta }}_d-\frac{\dot{\varphi }\dot{\phi }\left(I_z-I_x\right)}{I_y}-{\widehat{d}}_{\theta }\end{array}\right]$$

(42)

$$U_4=I_z\left[\begin{array}{l}\left(1-{k_5}^2+{\mu }_3\right)e_5+\left(k_5+k_6\right)e_6-k_5{\mu }_3{\gamma }_3+{\ddot{\phi }}_d-\frac{\dot{\varphi }\dot{\theta }\left(I_x-I_y\right)}{I_z}-{\widehat{d}}_{\phi }\end{array}\right]$$

(43)

where $k_i$, ${\mu }_i$ are constants greater than 0, and

$$\left\{\begin{array}{l}{e}_3={\theta }_d-\theta ,{\gamma }_2={\int }_0^{\tau }{e}_3\left(\tau \right)d\tau \\ e_4=k_3{e}_3+{\dot{\theta }}_d+{\mu }_2{\gamma }_2-\dot{\theta }\\ e_5={\phi }_d-\phi ,{\gamma }_3={\int }_0^{\tau }{e}_5\left(\tau \right)d\tau \\ e_6=k_5{e}_5+{\dot{\phi }}_d+{\mu }_3{\gamma }_3-\dot{\phi }\end{array}\right.$$

(44)

3.3. Position Controller Design

The process of deriving the control rate of the position loop is similar to the design process of the attitude angle. The control rate of the position loop controller can be obtained in the following way.

The first step is to define the position error as:

$$\left\{\begin{array}{l}{e}_7=Z_d-Z\\ e_9=X_d-X\\ e_{11}=Y_d-Y\end{array}\right.$$

(45)

The expected intermediate virtual control amount is:

$$\left\{\begin{array}{l}{e}_8=k_7{e}_7+{\dot{Z}}_d+{\mu }_4{\gamma }_4-\dot{Z}\\ e_{10}=k_9{e}_9+{\dot{X}}_d+{\mu }_5{\gamma }_5-\dot{X}\\ e_{12}=k_{11}{e}_{11}+{\dot{Y}}_d+{\mu }_6{\gamma }_6-\dot{Y}\end{array}\right.$$

(46)

According to the derivation process of the attitude angle control law, the position control law is:

$$\left\{\begin{array}{l}{U}_1=\frac{\left[\begin{array}{l}\left(1-{k_7}^2+{\mu }_4\right)e_7+\left(k_7+k_8\right)e_8-k_7{\mu }_4{\gamma }_4+{\ddot{Z}}_d+g-{\widehat{d}}_z\\ +\frac{mL\left(\begin{array}{l}\ddot{\beta }sin\beta cos\alpha +\ddot{\alpha }cos\beta sin\alpha +{\dot{\beta }}^{2}cos\beta cos\alpha \\ +{\dot{\alpha }}^{2}cos\beta cos\alpha -2\dot{\beta }\dot{\alpha }sin\beta sin\alpha \end{array}\right)}{M+m}\end{array}\right]\left(M+m\right)}{cos\varphi cos\theta }\\ u_x=\frac{\left[\begin{array}{l}\left(1-{k_9}^2+{\mu }_5\right)e_9+\left(k_9+k_{10}\right)e_{10}-k_9{\mu }_5{\gamma }_5+{\ddot{X}}_d\\ -{\widehat{d}}_x-\frac{mL\left({\dot{\beta }}^{2}sin\beta -\ddot{\beta }cos\beta \right)}{M+m}\end{array}\right]\left(M+m\right)}{U_1}\\ u_y=\frac{\left[\begin{array}{c}\left(1-{k_{11}}^2+{\mu }_6\right)e_{11}+\left(k_{11}+k_{12}\right)e_{12}-k_{11}{\mu }_6{\gamma }_6+{\ddot{Y}}_d-{\widehat{d}}_y\\ -\frac{mL\left(\begin{array}{l}\ddot{\beta }sin\beta sin\alpha -\ddot{\alpha }cos\beta cos\alpha +{\dot{\beta }}^{2}cos\beta sin\alpha \\ +{\dot{\alpha }}^{2}cos\beta sin\alpha +2\dot{\beta }\dot{\alpha }sin\beta cos\alpha \end{array}\right)}{M+m}\end{array}\right]\left(M+m\right)}{U_1}\end{array}\right.$$

(47)

3.4. Swing Angle Controller Design

The swing angle of the hanging load cannot be directly controlled due to its under-driving characteristics. It needs to be converted into a position signal for indirect control by the controller. The swing angle controller is also a neural network adaptive integral backstepping controller. The design process is similar to the above, therefore, we can omit the detailed process and describe the abbreviated process as follows

The first step is to define the pendulum angle error as:

$$\left\{\begin{array}{l}{e}_{13}={\alpha }_d-\alpha \\ e_{15}={\beta }_d-\beta \end{array}\right.$$

(48)

The expected intermediate virtual control amount is:

$$\left\{\begin{array}{l}{e}_{14}=k_{13}{e}_{13}+{\dot{\alpha }}_d+{\mu }_7{\gamma }_7-\dot{\alpha }\\ e_{16}=k_{15}{e}_{15}+{\dot{\beta }}_d+{\mu }_8{\gamma }_8-\dot{\beta }\end{array}\right.$$

(49)

According to the above derivation process, the swing angle control law is:

$$\left\{\begin{array}{l}{u}_{\alpha }=\frac{\left[\left(1-{k_{13}}^2+{\mu }_7\right)e_{13}+\left(k_{13}+k_{14}\right)e_{14}-k_{13}{\mu }_7{\gamma }_7+{\ddot{\alpha }}_d-{\widehat{d}}_{\alpha }-\frac{2sin\beta \dot{\beta }\dot{\alpha }}{cos\beta }\right]MLcos\beta }{U_1}\\ u_{\beta }=\frac{\left[\left(1-{k_{15}}^2+{\mu }_8\right)e_{15}+\left(k_{15}+k_{16}\right)e_{16}-k_{15}{\mu }_8{\gamma }_8+{\ddot{\beta }}_d-{\widehat{d}}_{\beta }+sin\beta cos\beta {\dot{\alpha }}^2\right]ML}{U_1}\end{array}\right.$$

(50)

4. Simulation Verification

In this paper, the QBall2 quadrotor carrying a load is used as the simulation object. The experimental platform consists of three parts: the hanging flight system composed of the QBall2 quadrotor and the hanging load produced by Quanser, the OptiTrack positioning system and the flight control system of MATLAB/Simulink environment. The experimental platform is shown in Figure 3.

Table 1. Quadrotor parameters.

Parameter	Value	Parameter	Value
M	1.80 kg	$I_x$	0.03 kg $·{\mathrm{m}}^2$
l	0.20 m	$I_y$	0.03 kg $·{\mathrm{m}}^2$
m	0.2 kg	$I_z$	0.04 kg $·{\mathrm{m}}^2$
L	0.4 m	$M_t$	8.80 N
w	15 rad/s	$M_y$	0.40 N
g	9.81 m/(${\mathrm{s}}^2$)

below shows the parameters of each part of the quadrotor. The selected load is 0.2 kg, and the hanging rope is 0.4 m.

The actual outdoor wind disturbance nonlinear model is more complicated. Considering the convenience of the research, the outdoor wind disturbance is generally simplified into a simple and practical engineering wind field model, and the complex parts of the model are simplified. The simplified model includes the main influencing factors in the wind field. It can be regarded as a real wind field model. In order to facilitate the research, this paper superimposes the Dryden turbulent wind field model with variable wind direction and the average wind with variable wind speed. The Dryden model is selected for the turbulent wind field. The structure of the model is relatively simple, and the time-frequency spectrum function can be conjugate decomposed. The input signal’s Gaussian white noise is generated through the designed shaping filter to generate the numerical model of the atmospheric turbulent wind field to generate the turbulent wind field sequence. The transfer function $G\left(s\right)$ of the shaping filter obtained by splitting and simplifying is:

$$\left\{\begin{array}{l}{G}_u\left(s\right)=\frac{{\sigma }_u\sqrt{L_u/\pi V_u}}{\left(L_u/V_u\right)s+1}\\ G_v\left(s\right)=\frac{{\sigma }_v\sqrt{L_v/\pi V_v}}{\left(2L_v/\sqrt{3}{V}_v\right)s+1}\\ G_w\left(s\right)=\frac{{\sigma }_w\sqrt{L_w/\pi V_w}}{\left(2L_w/\sqrt{3}{V}_w\right)s+1}\end{array}\right.$$

(51)

Among them, ${\sigma }_i\left(i=u,v,w\right)$ represents the intensity of turbulence, $L_i\left(i=u,v,w\right)$ represents the scale of turbulence and $V_i\left(i=u,v,w\right)$ represents the flight speed of each axis.

The average wind speed is 4.4 m/s, the median value of a level 3 breeze, as weak interference and the median value of level 4 wind, 6.7 m/s, is strong interference. The superimposed wind field model can indicate the change in the wind direction and intensity. Differently, it can be expressed as a real wind field model.

During the experiment, the Qball2 quadrotor started to take off and began to ascend. At the fifth second, the quadrotor reaches an altitude of 0.5 m and starts to track according to: $\left[\begin{array}{ccc}X& Y& Z\end{array}\right]=\left[\begin{array}{cc}\begin{array}{cc}2cos\left(\frac{\pi }{10}t\right)& -2sin\left(\frac{\pi }{10}t\right)+2\end{array}& 0.1t\end{array}\right]$ orbit. At the 15th second, add the superimposed wind field as interference to the eight-degree-of-freedom mathematical model of a quadrotor carrying a load, and compare two different strengths of interference according to the strength of the wind field. The Simulink simulations are compared with the adaptive integral backstepping algorithm and the backstepping algorithm for these two strength disturbances. The simulation time is set to 35 s. The simulation plots are shown below, the left column is (a) and the right column is (b).

The simulation $Bks$ denotes the backstepping algorithm, $Adp$ denotes the adaptive integral backstepping algorithm and $Rbfibks$ denotes the neural network adaptive algorithm based on integral backstepping. Each channel of the controller is compared and simulated using the neural network adaptive integral backstepping algorithm and the existing adaptive integral backstepping algorithm. In Figure 4 and Figure 5, the neural network adaptive integral backstepping algorithm can track the desired curve quickly and stably, and the tracking effect is closer to the desired curve than the adaptive integral backstepping algorithm and the backstepping algorithm. Among them, the green dashed line is the trajectory curve of the hanging load, while a quadrotor and its load are given for one sampling moment during the simulation.

From Figure 6 and Figure 7, it can be seen that the control effect of the adaptive integral backstepping algorithm does not differ much from that of the algorithm proposed in the paper after being affected by weaker wind disturbances, but the calculation of the average error still shows that the algorithm proposed in this paper is slightly better. However, when the system is subjected to strong wind disturbance, it is obvious from the figure that the actual trajectory of the algorithm in this paper fits better with the expected trajectory, and the error fluctuation is smaller and more stable.

Table 2. X position average error.

Interference	Bks	Adp	Rbf
Weak	0.1062 m	0.0536 m	0.0489 m
Strong	0.1878 m	0.1156 m	0.0962 m

below shows the comparison of the average error of different algorithms on the X-axis under different intensity interference. It is obvious from the table that the average error of the algorithm proposed in this paper is the smallest, and its performance is significantly improved compared to the integral backstepping algorithm and the backstepping algorithm. In the table, $Bks$ represents the backstepping algorithm, $Adp$ represents the adaptive integral backstepping algorithm, and $Rbf$ represents the neural network adaptive integral backstepping algorithm.

In Figure 8 and Figure 9, the pink curve represents the backstepping algorithm, the red curve represents the adaptive integral backstepping algorithm, and the blue curve represents the algorithm proposed in this article. It can be seen from the figure, all three methods can be stabilized quickly, but from the perspective of local scaling, the adaptive integral backstepping method and backstepping method produce larger overshoot, generate larger fluctuations, and take a longer time. The neural network adaptive integral backstepping algorithm can be stabilized quickly and with less overshoot. From the error point of view, the error fluctuation of the adaptive integral backstepping method is larger, which indicates that the addition of a neural network can better ensure the stable tracking of the trajectory in the case of external disturbances of different intensities.

Table 3. Y-position average error.

Interference	Bks	Adp	Rbf
Weak	0.0668 m	0.0181 m	0.0168 m
Strong	0.3190 m	0.1355 m	0.0596 m

below shows the comparison of the average error of different algorithms on the Y-axis under different interference intensities. It is obvious from the table that the average error of the algorithm proposed in this paper is the smallest, and its performance is significantly improved compared to the integral backstepping algorithm and the backstepping algorithm.

From Figure 10 and Figure 11, it can be seen that the height control difference of the three algorithms is not significant, but from the height error policy graph, the error fluctuation of the neural network algorithm is significantly smaller than that of the adaptive algorithm, tracking slightly better.

Table 4. Z-position average error.

Interference	Bks	Adp	Rbf
Weak	0.0022 m	0.0010 m	0.0008 m
Strong	0.0047 m	0.0016 m	0.0011 m

below shows the comparison of the average error of different algorithms in height under different intensity interference. It is obvious from the table that the average error of the algorithm proposed in this paper is the smallest, and its performance is significantly improved compared to the integral backstepping algorithm and the backstepping algorithm.

Figure 12 shows the attitude angles. The figure shows that the attitude angle section can follow the attitude response closely. It shows that the neural network adaptive integral backstepping controller can control the wind disturbance with almost no influence on the input control signal of the system.

In Figure 13, $rbfibks$ denotes the neural network adaptive algorithm based on integral backstepping, the pendulum angle generated by the hanging load during the flight is shown. The black line indicates that the pendulum angle controller is not added, and the red line indicates that the pendulum angle controller is added. From the figures, it can be seen that the overshoot of the pendulum angle generated by adding the neural network pendulum controller is significantly smaller, and it is fast and stable with less oscillation, so that the pendulum angle of the hanging load is quickly stabilized at a smaller value, which proves that the proposed algorithm can quickly stabilize the pendulum angle while suppressing the disturbance.

A comprehensive analysis of the above simulation results in MATLAB shows that after adding the neural network algorithm, a variety of performance indicators of the original algorithm can be improved. The comparison of simulation plots shows that the algorithm proposed in this paper has a higher overlap rate with the desired trajectory in terms of curve fit, and the stronger the wind interference, the more obvious the comparison effect of the algorithm, which also shows that the inclusion of the neural network has a better suppression effect on strong interference, proving that the problem proposed in the introduction can be effectively solved by the algorithm, which provides a reasonable basis for the completion of the thesis.

5. Conclusions

In this paper, we address the problems of model uncertainty, modeling errors and external complex wind field perturbation simulation by taking the hanging system of a quadrotor as the control target. Our innovation is that, first, the load and quadrotor are analyzed as a whole. In the environment of model parameter uncertainty, modeling error and external complex wind field disturbance, an eight-degree-of-freedom mathematical model of the quadrotor load, including the pendulum angle, is obtained. Second, the neural network algorithm is added to the adaptive integral backstepping algorithm to achieve the fusion of multiple algorithms. Third, the new algorithm is applied to the eight-degree-of-freedom model of the quadrotor load, and the robustness and effectiveness of the proposed algorithm are demonstrated through simulation comparison with the adaptive integral backstepping algorithm and the backstepping algorithm. The quadrotor can be better stabilized in a jamming environment, enabling it to track the desired trajectory quickly and stably. The tracking control simulation of the 3D spiral trajectory is carried out in MATLAB/Simulink. The simulation results show that in an environment with external wind field interference, the quadrotor can reach the specified position quickly with less error, and the suspension swing angle is better than the adaptive integral backstepping algorithm and the backstepping algorithm by adding the neural network algorithm based on the original algorithm. It can be stabilized faster and reduce the influence of load swing, which effectively reflects the suppression of interference by the neural network and proves the effectiveness of the proposed algorithm. The safety of the quadrotor needs to be improved, so no practical experiments are conducted in this paper. This will also be the future research direction of the authors.