Harmonic Extended State Observer Based Anti-Swing Attitude Control for Quadrotor with Slung Load

During the flight of the quadrotor, the existence of a slung load will exert a swing effect on the system and the motion of which will significantly change the dynamics of the quadrotor. The external torque caused by the slung load can be considered as a kind of disturbance and it is a threat to the attitude control stability of the system. In order to solve this problem, a high precision disturbance compensation method is presented in this paper, based on the harmonic extended state observer (HESO). Firstly, a generic mathematical model for the quadrotor-slung load system is obtained via the Lagrangian mechanics, and according to the analysis of the slung load motion, we obtain the disturbance as a form of periodic equation. Secondly, based on the dynamic model of the disturbance, we propose a HESO to achieve high precision disturbance estimation and its stability is proved by Lyapunov theory. Thirdly, we designed an attitude tracking controller based on backstepping method, and discussed the stability of the entire system. Finally, numerical simulations and real time experiments are carried out to evaluate the performance of the proposed method. Our results show that the robustness of the quadrotor subject to slung load has been improved.


Introduction
In recent years, the research on quadrotor vehicles has attracted great interest due to the wide range of civil and military applications, and many achievements have been made [1,2].As a new kind of unmanned aerial vehicle (UAV), the quadrotor is a small rotorcraft with four propellers driven by four direct current (DC) motors respectively [3].Compared with traditional helicopters, the structure of the quadrotor is simpler and more efficient, and has significant advantages in precise hovering, aggressive maneuver, vertical take-off and landing (VTOL) [4,5], etc.
Like traditional helicopters, quadrotor vehicles have many important applications in carrying slung load, such as deploying supplies in military operations, or delivering first-aid kits for personal assistance to the victims in disasters like floods, earthquakes, fires, industrial accidents, etc. [6,7], and the research work addressing quadrotor vehicles with slung load becomes an attractive topic.On this application, the major difficulty in modeling and control study is the coupled effects between the quadrotor vehicles and the slung load [8].The external slung load behaves like a pendulum and the motion of which will significantly change the dynamics of the quadrotor.Moreover, the quadrotor vehicle is a typical underactuated, strong coupled, nonlinear system [9], and inherently unstable without close-loop controller [10].If the pendulous motion of the load exceeds certain limits, the stability of control system will be broken because of the changes in dynamic characteristics of the plant.Therefore, improving the robustness of the controller subject to the oscillation of slung load is very necessary.
In the relevant literatures on the quadrotor, the problem of addressing quadrotor-slung load system has been discussed in some publications.The mathematical model of the entire system is derived by the Newton-Euler formulation in [8,11], however, there is no further analysis on the dynamic characteristics of the slung load motion in detail.In [12], the slung load is modeled as a point mass spherical pendulum, and a related adaptive control method is proposed to handle the additional forces and torques acting on the quadrotor.In [13], a generalized approach is presented using an iterative optimal control algorithm, and a series of complex tasks are thus solved without the need for manual manipulation of the system dynamics, heuristic simplifications, or manual trajectory generation.In [14], a nonlinear dynamic model is presented, and an interconnection and damping assignment-passivity based Control (IDA-PBC) methodology is used for precise payload's positioning with stabilization of the swing angles to the minimum of the desired energy function.
Alternatively, the effect of slung load on the quadrotor can be considered as a kind of disturbance [15].Therefore, in order to improve system reliability and achieve the requirements of high precision control, the disturbance and uncertainty estimation and attenuation (DUEA) method would be a potential solution.In recent years, this method has been widely used and achieved good results [16,17] The framework of DUEA can be divided into two parts, namely, a disturbance and uncertainty estimator (DUE) and a feedback controller (FC).In the first part, DUE is designed to estimate the disturbances so that they could be compensated in the feedforward loop.Then the FC in the second part is designed to guarantee fast convergence of the closed-loop system.In this framework the DUE plays an important role because the performance of the closed-loop system is largely determined by the estimation accuracy of it.Therefore, a series of observers have been proposed as the DUE so far to improve estimation accuracy under different conditions, such as disturbance observer (DO) [18,19], extended state observer (ESO) [15,20,21] and proportional integral observer (PIO) [22,23] etc.By the appropriate use of the observer, disturbance rejection performance and robustness of the existing control system could be significantly improved.
The extended state observer (ESO), known as the key module of active disturbance rejection control, can estimate both the states of system and the total disturbances with less dependence on model information [24,25].This method was first proposed by Han in 1990s and the basic idea behind ESO is to view disturbance as an extended state and utilize observer to estimate it [24].Traditionally, ESO approaches focus primarily on dealing with slowly changing disturbances.However, it's obvious that the disturbance caused by the slung load is periodic, which cannot be estimated by traditional ESO thoroughly [26].Therefore, an enhancement of ESO that can handle periodic disturbance is necessary in this field.In [27], a higher-order ESO is investigated, from the results, it can be seen that the higher-order ESO can improve the estimation accuracy of sinusoidal external disturbances more or less, while, there still exists a periodic estimation error that will in turn decrease the control accuracy of the closed-loop system.Furthermore, the higher level of the observer order will lead to a higher observer gain, which will in return excite the sensor noise and introduce them into the control loop.The internal model principle is applied for generalized ESO in [28] and harmonic disturbance observer (HDO) in [29] by embedding the disturbance dynamics into the observer, and the disturbance estimation performance is improved.
Motivated by these methods, a harmonic extended state observer based anti-swing attitude control method is proposed for a quadrotor subject to slung load in this article, where the periodic disturbance caused by slung load motion is compensated by the estimated state signals produced of the HESO.The main contributions of this paper are summarized as follows: (1) Build the generic mathematical model for the quadrotor-slung load system by the Lagrangian mechanics, and analyze the characteristics of the slung load motion in detail.(2) Propose a HESO based on the characteristics of the slung load motion.
The outline of this paper is as follows: The mathematical model and the control problems of quadrotor-slung load system are formulated in Section 2. A HESO is designed in Section 3. In Section 4, an attitude tracking controller is designed via backstepping method.Numerical simulation and real time experimental results are presented in Section 5, and the conclusions are summarized in Section 6.

Notations and Assumptions
Throughout this paper, the following notations will be used.R is the set of real numbers.Let • denote the 2-norm of a vector or a matrix.For a given vector and for a given matrix A ∈ R n×n , A = λ max A T A , where λ max (•) is the maximal eigenvalue of the matrix.In addition, the operator S(•) maps a vector T to a skew symmetric matrix as: In this section, the mathematical model of the quadrotor-slung load system is established.As shown in Figure 1, we consider this system consists of three parts: a quadrotor, a cable and a payload.And before our work, the following assumptions are made: (1) The slung load is considered as a particle and only swings in a plane.The length of the connecting cable is constant and known.(2) The inelastic cable is massless and always tight and no consideration of the energy loss caused by the friction force in the swing.(3) The aerodynamic effects on the load are neglected.

Quaternion Operations
In order to avoid the singularity problem of trigonometric functions, unit quaternion q = q 0 q T v T ∈ R 4 , q = 1 is used to represent rotation [30].Following are the operations we used.
The quaternion multiplication is: q 01 q 02 − q T v1 q v2 q 01 q v2 + q 02 q v1 − S(q v2 )q v1 (2) The relationship between rotation matrix R B A and q is calculated as: .
The derivative of a quaternion is given by the quaternion multiplication of the quaternion q and the angular velocity of the plant ω: The quaternion error q e is given as the quaternion multiplication of the conjugate of the desired quaternion q d and the actual quaternion q: q e = q * d ⊗ q = q 0e q ve = q 0 q 0d + q v T q vd q 0d q v − q 0 q vd + S(q v )q vd (6) where q * = q 0 −q v T is the conjugate of the quaternion q.
2.2.Mathmetical Model of Quadrotor-Slung Load system Firstly, the quadrotor can be considered as a rigid cross frame attached with four rotors, and the center of gravity coincides with the body-fixed frame origin.The simplified model of the quadrotor is presented in Figure 1, rotors R1 and R3 rotate counterclockwise, and rotors R2 and R4 rotate clockwise, each propeller rotates at the angular speed Ω i and produces a force F i (i = 1,2,3,4) along the negative z-direction relative to the body frame [30,31]: where k T > 0 denotes the aerodynamic coefficient which consists formed of the atmospheric density ρ, the radius of the propeller r, and the thrust coefficient c T .In addition, due to the spinning of the rotors, a reaction torque M i (i = 1,2,3,4) is generated on the quadrotor body by each rotor: where k D > 0 denotes the drag coefficient of the rotor, which depends on the same factors as k T .
The variation of the orientation is achieved by varying the angular speed of a specific rotor.The total force and torque acting on the quadrotor are defined as follows: where u represents the attitude control signal to be designed.
Electronics 2018, 7, x FOR PEER REVIEW 4 of 18 q q q q q q q q q q (6) where is the conjugate of the quaternion q .

Mathmetical Model of Quadrotor-Slung Load system
Firstly, the quadrotor can be considered as a rigid cross frame attached with four rotors, and the center of gravity coincides with the body-fixed frame origin.The simplified model of the quadrotor is presented in Figure 1, rotors R1 and R3 rotate counterclockwise, and rotors R2 and R4 rotate clockwise, each propeller rotates at the angular speed i Ω and produces a force Fi (i = 1,2,3,4) along the negative z-direction relative to the body frame [30,31]: where denotes the aerodynamic coefficient which consists formed of the atmospheric density ρ , the radius of the propeller r , and the thrust coefficient T c .In addition, due to the spinning of the rotors, a reaction torque Mi (i = 1,2,3,4) is generated on the quadrotor body by each rotor: where denotes the drag coefficient of the rotor, which depends on the same factors as T k .
The variation of the orientation is achieved by varying the angular speed of a specific rotor.The total force and torque acting on the quadrotor are defined as follows: where u represents the attitude control signal to be designed.In the mathematical model of quadrotor, three coordinate frames are considered: the non-moving inertial coordinate frame E I : {o I , x I , y I , z I }, the body-fixed coordinate frame E B : {o B , x B , y B , z B } and the desired frame E D : {o D , x D , y D , z D } to represent the actual attitude and desired attitude of quadrotor respectively.Note that NED coordinates are used to define all frames.In the inertial frame E I , the position of the quadrotor is X q = x q y q z q T ∈ R 3 , the attitude angle is Θ = φ θ ψ T ∈ R 3 and the quaternion expression of the attitude is q = [ q 0 q v ] T ∈ R 4 .In the body-fixed frame E B , the angular velocity of is ω = [ ω x ω y ω z ] T ∈ R 3 .Thus, the rotation matrix from E B to E I can be represented as: and the rotation velocity transfer matrix can be given as: Thus, the rotational kinematic equations with respect to the inertial frame E I can be expressed as: .
Secondly, we study the slung load.The position of the hook point is defined as c = 0 0 c T in body frame E B .According to previous assumptions, the relationship between the position of the slung load X l = x l y l z l T ∈ R 3 and X q are: where γ ∈ (−90 • , 90 • ) is the pendulum angle between the cable and the positive orientation of o I z I , β is the angle between the pendulum plane and the x I o I z I plane, and L c is the length of the cable.Moreover, according to the former assumptions, β and L c are constant.Finally, we use Lagrangian mechanics to summarize the dynamics of the quadrotor-slung load system.Compared with the Newtonian mechanics, this choice eliminates the need for the constraint forces to enter into the resultant system of equations, so that fewer equations are needed.As shown in Figure 1, the entire system has seven degrees of freedom (DOF): pendulum angle of the slung load γ, position of the quadrotor X q = x q y q z q T and attitude of the quadrotor . Therefore, we defined the generalized coordinate of the quadrotor-slung load system as η = x q y q z q φ θ ψ γ T .The expressions for the kinetic and potential energies will be presented in order to obtain the Lagrangian of the system.The total kinetic energy function of the quadrotor-slung load system, resulting from the translational and rotational motions can be portioned as the sum of the translational kinetic energy: And the rotational kinetic energy of the entire system: where the inertia matrix J = diag(J x , J y , J z ) is diagonal.
And the total potential energy function of the system results from the sum of the potential energies of the quadrotor and the slung load: Thus, the Lagrangian of the generalized system can be defined as: Notice that gravity G is the conservative force in this system, while total lift F = and torque u are non-conservative forces, the Euler Lagrangian equations of the second kind can be obtained as: where is the non-conservative generalized force.When the quadrotor is carrying a slung load, Large-scale maneuvers should be avoided, therefore, the coupling of angular velocity between different channels can be ignored.Then, we assume P ≈ I 3 when calculating the dynamic quadrotor-slung load system.According to Equations ( 14)-( 18), the dynamic model can be expressed as: ∑ R : ..

γ.
According to [30], we can obtain the attitude kinematic of the quadrotor as:

Analysis of Slung Load Motion
In this article, we consider the effect of slung load on the quadrotor as a kind of disturbance.In order to realize high precision estimation of the disturbance, the characteristics of the disturbance need to be analyzed.
The derived equations of motion of the load and quadrotor are highly nonlinear and strongly coupled, thus, they are difficult to be used for motion analysis.Therefore, we trimmed the mathematical near hovering or uniform linear flight.In this condition, we have .. X q ≈ 0. And in actual flight c L c , then we have cL −1 c ≈ 0. From Equation (20), we can see that the load movement can hardly affect ψ, therefore, disturbance in pitch and roll channels are considered in this condition.The trimmed results of the rotational dynamic are: where Then, analyze the motion of slung load.When γ is a small angle, we have sin γ ≈ γ and cos γ ≈ 1.Thus, we can rewrite Equation (23) as: ..
According to Equation ( 24), we can calculate that: where A 0 and χ 0 are the unknown amplitude and phase related to initial conditions.Substituting Equation (25) into Equation ( 23), we have: where ∆ i (i = φ, θ) are the unmodeled residual of the disturbance in each channel.
According to the above analysis, the disturbance torque caused by the slung load is a sum-of-sinusoids function which has two frequencies 1 = gL −1 c , 2 = 3 gL −1 c , while the amplitude A 1,i , A 2,i and phase χ 1,i , χ 2,i (i = φ, θ) are unknown.

Problem Formulation
In order to study the transient and steady-state characteristics of the quadrotor, the dynamics of attitude error are introduced.We use T and q d = [ q 0d q vd ] T to denote the desired angular velocities and attitude respectively, thus: where ω e = [ ω e,x ω e,y ω e,z ] T is the tracking error vector of the angular velocities.Then, we can obtain the dynamics of ω e according to Equations ( 4), ( 20) and ( 27): .
where J 0 = diag(J x + m l c 2 , J y + m l c 2 , J z ), R B D can be calculated according to Equations ( 3) and (6).And according to Equations ( 5), ( 6) and ( 27), we can obtain the kinematics of attitude tracking error .
The problem we try to tackle in this work is to design a continuous control law u using only the measurable system output ω and q such that the error of attitude ω e and q e converge to zero in presence of the slung load.In order to ensure the robustness of the controller, the DUEA strategy is necessary Figure 2 illustrates the control scheme that we designed.Based on the DUEA control methodology, the attitude control problem for quadrotor can be divided into two components: design the feedforward loop so that the periodic disturbance is estimated by HESO and compensated this way; and design the feedback loop that regulates the orientation to track the desired attitude produced by the commander timely.Therefore, the control signal u contains two parts as: where u N is the nominal control input vector and u E is the disturbances attenuation input vector.compensated this way; and design the feedback loop that regulates the orientation to track the desired attitude produced by the commander timely.Therefore, the control signal u contains two parts as: where N u is the nominal control input vector and E u is the disturbances attenuation input vector.

Design and Stability Analysis of HESO
In this section, the design of HESO, which provides the disturbance estimate for the controller, is described in detail.As for the DUEA control methodology, the control performance of closed loop system will be largely determined by the observation performance.However, the disturbances acting on quadrotor are periodic, which cannot be accurately estimated by traditional ESO thoroughly [28].Therefore, in order to enhance the performance of feedback controller, a HESO is designed to estimate the periodic disturbances, and the stability analysis is carried out afterwards.

Design of HESO
According to the analysis of the disturbances in previous section, we can conclude that the external disturbances caused by the slung load have the same characteristic in roll and pitch channel.
Without loss of generality, the design process in i ( ) Firstly, we can rewrite the time varying disturbance as: sin sin where ≤ .In order to establish the model of equivalent disturbance, we write Equation (31) as a form of state equation where

Design and Stability Analysis of HESO
In this section, the design of HESO, which provides the disturbance estimate for the controller, is described in detail.As for the DUEA control methodology, the control performance of closed loop system will be largely determined by the observation performance.However, the disturbances acting on quadrotor are periodic, which cannot be accurately estimated by traditional ESO thoroughly [28].Therefore, in order to enhance the performance of feedback controller, a HESO is designed to estimate the periodic disturbances, and the stability analysis is carried out afterwards.

Design of HESO
According to the analysis of the disturbances in previous section, we can conclude that the external disturbances caused by the slung load have the same characteristic in roll and pitch channel.Without loss of generality, the design process in i(i = φ, θ) channel is demonstrated in detail.Firstly, we can rewrite the time varying disturbance as: ∆ i = δ i is upper bounded δ i ≤ σ i .In order to establish the model of equivalent disturbance, we write Equation (31) as a form of state equation . where is the system state and δ d,i is the perturbation term of the disturbance.
Then, consider the SISO nonlinear Equation (23), we have Define as the new state vector, the reconstructed system can be written as: Note that in the reconstructed model of quadrotor, it can be verified that the pair (A, C) is observable [28].
From the reconstructed system model in Equation ( 34), the observer is designed as follows: .
where xi = is the stats of the observer, ŷi is the estimate of the output y i , and T is the observer gain to be designed.
By defining the estimation error of the observer as x i = x i − xi , we have the following equation: .
From Equation (36), we can find that the estimation accuracy is partly determined by the perturbation term δ i , thus in order to increase the estimation accuracy of the observer, δ i should be decreased.Compared with the traditional higher order ESO, the equivalent disturbance model is introduced in HESO which makes the model more precise and the perturbation δ d,i smaller [28].

Convergency Analysis of HESO
Frist, define the Lyapunov candidate function V 0,i (i = φ, θ) as: Take the time derivative of V 0,i , we have: .
Then, consider Equation (38).Under the condition that A − L e,i C is Hurwitz, for any k 0 > 0 there exists a positive definite symmetric matrix P 0 satisfying: Substituting Equation (39) into Equation (38), we can get: .
where δ d,i is bounded with δ d,i ≤ σ i , and λ max (P 0 ) is the maximum eigenvalue of P 0 .It is obvious that .
V 0,i < 0 whenever x d,i > 2σ i λ max (P 0 )k −1 0 .Therefore, the upper bound for estimation error x i will be constrained by the bounded ball B r = r| r ≤ 2σ i λ max (P 0 )k −1 0 .Moreover, the estimation error decreases with the increase of the model accuracy.

Design of Attitude Controller
In this section, the main procedures of attitude controller integrated with HESO are presented for effectively handling disturbance caused by slung load.The quadrotor UAV is an underactuated system with six DOF and four control inputs.In order to derive its model, backstepping method is used in the design of attitude controller.
Step 1. Design the control strategy to ensure that q e (t) converges to zero.According to the attitude error kinematics subsystem Equation ( 29), we select the candidate Lyapunov function as: V 1 = q e T q e + (1 Take the time derivative of V 1 , we have: .
. q 0e = (q T ve S(q ve ) + q 0e q T ve I 3 + (1 ve ω e (42) Then, we design a virtual control scheme as: where K 1 is the gain matrix of the controller, which is diagonal positive definite.If the angular velocity tracking error ω e is equal to the virtual control input ω ed , .
According to the Lyapunov stability theorem, we can conclude that q e converges to zero, under the condition that the virtual control ω e converges to −K 1 q e .
Step 2. Design the control signal u to ensure that ω e track the desired virtual control input ω ed .We define the error between ω e and ω ed as: In order to discuss the stability of the entire system including DUE and FC, we define the following candidate Lyapunov function V 2 as: Take the time derivative of V 2 , and substitute Equations ( 28), (45), and (47) into .
V 2 , then we have .
V 2 = q T ve ω e − q T ve K 1 q ve + ω T e J 0 . .
Define the nominal control input vector as: .
) and λ min (K 2 ) are the minimal eigenvalue of K 1 and K 2 respectively.Define z T = (q T e , ω T e , x T φ , x T θ ) as the uniformed vector of errors, and ξ = min λ min (K 1 ), > 0, then Equation (49) can be reduced to .
V 3 < 0 whenever z > 2σλ max (P 0 )ξ −1 .Notice that (q e , ω e ) is a linear diffeomorphism of (q e , ω e ), hence (q e , ω e ) can converge into a compact set.We can conclude that, the attitude error (q e , ω e ), virtual control input error ω e and the estimation error x φ and x θ are uniformly ultimately bounded and exponentially converges to the bounded ball In general, when we chose a larger ξ, the bounded ball B z will become smaller and consequently, z will also become smaller, so a larger ξ is preferred.However, larger ξ will lead to larger control gains which can excite the sensor noise and undesirable high frequency dynamics of the system.Thus, the tuning of controller parameters is a tradeoff between the demand of performance and the real conditions.Moreover, by improving the model accuracy of the disturbance, σ will be reduced, therefore, the accuracy of attitude control can be improved.

Simulation and Experimental Results
In order to evaluate the performance of the proposed HESO, simulation and real world experimental results are presented in this section.

Simulation Results
We present the numerical simulations of the proposed GESO based DUEA control strategy on a model generated by the online toolbox of Quan and Dai [32], and the values of the nominal model parameters are list in Table 1.Length of the cable 1 m The values of gain parameters used in our controller are given as K 1 = 7I 3 and K 2 = 2I 3 , and the observer gains of HESO are L e,i = [ 6 27 27 8 8 2 2 ] T .Traditional 2nd order ESO is used as the comparison, and the bandwidth is fix at 3 rad/s, the observer gain is L = 9 27 27 T .

Comparison in Attitude Stabilization Performance
This part involves attitude stabilization control in the presence of slung load.The initial condition of the slung load is γ 0 = 30 • and β = 20 • .The initial conditions used of the quadrotor in the simulation are zero, and we select the desired attitude signal as follows Three comparative simulations were conducted, and the performances of the proposed controller with different ESOs are compared in pitch and roll channel.Figures 3 and 4 show the curves of the vehicle's attitude response during its flight.We can see that although the proposed controller (without GESO) was able to ensure the stabilization of the attitude angles, the control accuracy was reduced under the influence of the slung load motion.When we introduced Traditional 2nd order ESO and HESO as the DUE, the performances of the proposed controller were improved, and through the comparison of attitude control results, we can conclude that HESO performed better than traditional ESO.Moreover, Figures 5 and 6 show the disturbance estimation result of HESO and Traditional ESO.Obviously, we can see that the estimation error decreases when we adopt HESO. .
Three comparative simulations were conducted, and the performances of the proposed controller with different ESOs are compared in pitch and roll channel.Figures 3 and 4 show the curves of the vehicle's attitude response during its flight.We can see that although the proposed controller (without GESO) was able to ensure the stabilization of the attitude angles, the control accuracy was reduced under the influence of the slung load motion.When we introduced Traditional 2nd order ESO and HESO as the DUE, the performances of the proposed controller were improved, and through the comparison of attitude control results, we can conclude that HESO performed better than traditional ESO.Moreover, Figures 5 and 6 show the disturbance estimation result of HESO and Traditional ESO.Obviously, we can see that the estimation error decreases when we adopt HESO.

Comparison in Attitude Tracking Performance
In this case, the numerical simulation demonstrates the effectiveness of the proposed control scheme for attitude tracking.The initial condition of the slung load is γ 0 = 0 • .We chose desired attitude signal as φ d = 0 • and The attitude tracking performances of the proposed controller with different ESOs are illustrated in Figures 7 and 8.We can see that the proposed controller alone was hard to ensure the tracking of the desired attitude angles, and the control accuracy was severely affected by the slung load.When we introduced Traditional 2nd ESO and HESO as the DUE, the performances of the proposed controller were improved, and the simulation results show that HESO performed better than traditional ESO.Moreover, as shown in Figure 9, we can see that the estimation error decreases when we adopt HESO.illustrated in Figures 7 and 8.We can see that the proposed controller alone was hard to ensure the tracking of the desired attitude angles, and the control accuracy was severely affected by the slung load.When we introduced Traditional 2nd ESO and HESO as the DUE, the performances of the proposed controller were improved, and the simulation results show that HESO performed better than traditional ESO.Moreover, as shown in Figure 9, we can see that the estimation error decreases when we adopt HESO.The attitude tracking performances of the proposed controller with different ESOs are illustrated in Figures 7 and 8.We can see that the proposed controller alone was hard to ensure the tracking of the desired attitude angles, and the control accuracy was severely affected by the slung load.When we introduced Traditional 2nd ESO and HESO as the DUE, the performances of the proposed controller were improved, and the simulation results show that HESO performed better than traditional ESO.Moreover, as shown in Figure 9, we can see that the estimation error decreases when we adopt HESO.

Experimental Results
We have also tested the proposed control scheme on a self-assembled GF360 quadrotor, and PIXHAWK [33,34] was used as the autopilot of the quadrotor.To evaluate the stability and robustness of the proposed control scheme, the experiments were carried out as follows.The weight of the slung load l m is 0.5 kg and cable length c L is 1 m, the distance between the hook of the

Experimental Results
We have also tested the proposed control scheme on a self-assembled GF360 quadrotor, and PIXHAWK [33,34] was used as the autopilot of the quadrotor.To evaluate the stability and robustness of the proposed control scheme, the experiments were carried out as follows.The weight of the slung load m l is 0.5 kg and cable length L c is 1 m, the distance between the hook of the cable and the mass center of quadrotor c is about 20 cm.
Subject to the sampling frequency and noise of MEMS gyro, the observer gain of HESO is chosen as Figure 10 presents the history of pitch angle tracking errors θ e = θ − θ d obtained by the proposed controller with different ESO in the presence of slung load.Figure 11 shows the experiment results of pitch rate ω y in these three cases.It can be observed that the control performance is improved by HESO compared with the one without ESO and the one with traditional 2nd order ESO.In addition.The root mean square (RMS) errors obtained by proposed controller with different ESO are given in Table 2.It is illustrated that HESO can achieve 50.11% reduction on the RMS error of pitch angle.

Experimental Results
We have also tested the proposed control scheme on a self-assembled GF360 quadrotor, and PIXHAWK [33,34] was used as the autopilot of the quadrotor.To evaluate the stability and robustness of the proposed control scheme, the experiments were carried out as follows.The weight of the slung load l m is 0.5 kg and cable length c L is 1 m, the distance between the hook of the cable and the mass center of quadrotor c is about 20 cm.Subject to the sampling frequency and noise of MEMS gyro, the observer gain of HESO is chosen as

Conclusions
In this paper, the problem of anti-swing attitude control for quadrotor in the presence of slung load is investigated.In the proposed approach, we developed a HESO based DUEA control scheme to accurately attenuate the disturbance torque caused by slung load.Through the analysis of the quadrotor-slung load system, the external disturbance torque can be described as a periodic function with two frequencies 1 = gL −1 c , 2 = 3 gL −1 c , while the amplitude and phase are unknown.Then the harmonic extended state observer is designed according to the model of the disturbance in propose of enhancing the robustness of feedback nonlinear controller.Subsequently, an tracking controller is designed based on the backstepping method.From the simulation and experimental results, we can conclude that compared to the traditional 2nd order ESO, the HESO can achieve a better performance in estimating the periodic disturbances, so that the robustness of the proposed controller can be further improved.

γ 2 sin
γ sin β and d θ = −m l cL c .. γ cos γ − .γ 2 sin γ cos β are the torque disturbance caused by the slung load in roll and pitch channel respectively.

Figure 2 .
Figure 2. Block diagram of the proposed control scheme.

Figure 2 .
Figure 2. Block diagram of the proposed control scheme.

5. 1 . 1 .
Comparison in Attitude Stabilization Performance This part involves attitude stabilization control in the presence of slung load.The initial condition of the slung load is 0 30 γ =  and 20 β =  .The initial conditions used of the quadrotor in the simulation are zero, and we select the desired attitude signal as follows [

Figure 3 .
Figure 3. Simulation curves of the attitude angle in attitude stabilization: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 3 .
Figure 3. Simulation curves of the attitude angle in attitude stabilization: (a) roll channel; (b) pitch channel; (c) yaw channel.Electronics 2018, 7, x FOR PEER REVIEW 13 of 18

Figure 4 .
Figure 4. Simulation curves of the attitude rate in attitude stabilization: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 4 .
Figure 4. Simulation curves of the attitude rate in attitude stabilization: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 4 .
Figure 4. Simulation curves of the attitude rate in attitude stabilization: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 5 .
Figure 5. Estimation results comparison of harmonic extended state observer (HESO) and traditional extended state observer (ESO) in attitude stabilization.

Figure 6 .
Figure 6.Estimation errors comparison of HESO and traditional ESO in attitude stabilization.

Figure 5 .
Figure 5. Estimation results comparison of harmonic extended state observer (HESO) and traditional extended state observer (ESO) in attitude stabilization.

Figure 4 .
Figure 4. Simulation curves of the attitude rate in attitude stabilization: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 5 .
Figure 5. Estimation results comparison of harmonic extended state observer (HESO) and traditional extended state observer (ESO) in attitude stabilization.

Figure 6 .
Figure 6.Estimation errors comparison of HESO and traditional ESO in attitude stabilization.Figure 6. Estimation errors comparison of HESO and traditional ESO in attitude stabilization.

Figure 6 .
Figure 6.Estimation errors comparison of HESO and traditional ESO in attitude stabilization.Figure 6. Estimation errors comparison of HESO and traditional ESO in attitude stabilization.

Figure 7 .
Figure 7. Simulation curves of the attitude angle in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 8 .
Figure 8. Simulation curves of the attitude rate in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 7 .
Figure 7. Simulation curves of the attitude angle in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 7 .
Figure 7. Simulation curves of the attitude angle in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 8 .
Figure 8. Simulation curves of the attitude rate in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.Figure 8. Simulation curves of the attitude rate in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.

Figure 8 . 18 Figure 9 .
Figure 8. Simulation curves of the attitude rate in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.Figure 8. Simulation curves of the attitude rate in attitude tracking: (a) roll channel; (b) pitch channel; (c) yaw channel.Electronics 2018, 7, x FOR PEER REVIEW 15 of 18

Figure 9 .
Figure 9. Estimation errors comparison of HESO and traditional ESO in attitude tracking.

12 8 12 8 2 2 ]
T .And subject to the limitations of computing power, HESO is only used in pitch channel, and the slung load only swing in longitudinal plane.The controller gains were chosen as K 1 = diag[ 7 7 2.8 ] and K 2 = diag[ 0.15 0.15 0.2 ].Three comparative experiments with different ESO are presented.

Figure 9 .
Figure 9. Estimation errors comparison of HESO and traditional ESO in attitude tracking.

1 diag
to the limitations of computing power, HESO is only used in pitch channel, and the slung load only swing in longitudinal plane.The controller gains were chosen as experiments with different ESO are presented.Figure 10 presents the history of pitch angle tracking errors e d θ θ θ = − obtained by the proposed controller with different ESO in the presence of slung load.Figure 11 shows the experiment results of pitch rate y ω in these three cases.It can be observed that the control performance is improved by HESO compared with the one without ESO and the one with traditional 2nd order ESO.In addition.The root mean square (RMS) errors obtained by proposed controller with different ESO are given in Table2.It is illustrated that HESO can achieve 50.11% reduction on the RMS error of pitch angle.

Figure 11 .
Figure 11.Experimental curves of ω y with different ESO.

Table 1 .
Quadrotor parameters used in simulation.

Table 2 .
Comparison of control performance in the tracking error of pitch angle with different ESO (RMS error).Without ESO Traditional ESO HESO θe (°) 2.4696 1.8890 1.2320

Table 2 .
Comparison of control performance in the tracking error of pitch angle with different ESO (RMS error).