Attitude Control Method of Unmanned Helicopter Based on Asymmetric Tracking Differentiator and Fal-Extended State Observer

In order to meet the constraints of velocity and acceleration of displacement and attitude motion of an unmanned helicopter during an automatic carrier landing mission, an asymmetric tracking differentiator, which could set the speed and acceleration limits in two directions of tracked signal motion respectively, was derived based on the tracking differentiator in the active disturbance rejection control method. Based on the proposed asymmetric tracking differentiator, a fal-extended state observer based on the fal function was added to construct the attitude and angular velocity controller which is programmable during the transition process of unmanned helicopters. The mathematical simulation and result analysis show that the newly proposed attitude estimation algorithm effectively compensates for the deficiencies of the existing methods, improves the antijamming capability and the accuracy of attitude estimation in the maneuvering process, achieving the expected design purposes.


Introduction
In modern military actions, unmanned helicopters are widely used and play an increasingly important role. Carrier-borne helicopters of UAV cannot only take off and land on large aircraft carriers in multiple battle groups, but also take off and land on medium and small ships and fly over the ocean. It can fly singularly or in a group over the ocean to carry out escort, anti-submarine, anti-ship and beyond-visual-range (BVR) guidance and other combat missions. It can fly and hover, rapidly and slowly, highly and lowly, back and forth, left and right, and is easy to use. In addition, it demonstrates a faster speed than naval vessels at sea, a longer flight range than aircrafts on land, higher flexibility than other carrier-borne aircrafts. All of these combat advantages, are attracting more and more attention from navies of other countries around the globe. Attitude control is the basis of motion control of an unmanned helicopter. Higher attitude control accuracy and better transition process quality enables helicopters to accurately and rapidly execute the attitude command generated by the guidance system during flight and ensure the accuracy of aircraft movement. The flight control system of an unmanned helicopter is a strong-coupled nonlinear system with multiple inputs and outputs. The stability boundary of the control system changes with the flight conditions. The control technology of unmanned helicopters can be divided into linear control and nonlinear control methods. The linear control method is simple in design and easy to implement, but it has large limitations. It is generally used to ensure the stability of helicopters under a specific flight mode. The nonlinear control method is complicated and difficult to implement, and requires a certain precision of the nonlinear modeling, but it can make the helicopter effective in a wider range of flight modes. There are many kinds of UAV helicopter control methods, such as PID control, neural network control, fuzzy control; control methods have their own characteristics for different systems and applications. The dynamic inversion [1] method has been widely applied in flight control systems. The advantages of the dynamic inversion method are that it is simple, intuitive and easy to understand, and the system model targeted by dynamic inversion is not limited by affine nonlinear conditions. At the same time, the dynamic inversion method also has some shortcomings. It is suitable for nonlinear systems where the nonlinear part is solvable and it also depends on the model precision to a certain extent. Andrew Richard Conway [2] used four independent PD controllers to realize the control of the altitude, direction, forward flight and crabbing of the unmanned helicopter. David H. Shim et al. [3] adopted cascade split loop control to realize hovering control of an unmanned helicopter. Kyoto University in Japan [4] realized the hovering and autonomous flight of an unmanned helicopter by using a neural network controller. Tokyo Institute of Technology [5] designed a flight control system based on fuzzy control and successfully realized the autonomous takeoff and landing of unmanned helicopters.
Two-degree-of-freedom control based on disturbance observer, also known as disturbance observation and suppression technology [6], can effectively solve the conflict between disturbance rejection performance and tracking performance in single-degree-offreedom feedback control. Since the 1960s, a series of disturbance estimators have been proposed, among which the design of Extended State Observer (ESO) requires the least model information. It is the core part of Active Disturbance Rejection Control (ADRC) and Tracking Differentiator (TD) proposed by Han J.Q. in 1998 [7]. In the filtering characteristics of ESO, some good results have also been achieved in recent years. Although the designed ESO can suppress measurement noise to a certain extent, it still faces many difficulties in dealing with time-delay type systems, especially the phase lag of signals. On the basis of analyzing the essence of Smith's predictive control of time-delay systems, Han proposed an algorithm with strong noise suppression ability, which realized the functions of "phase lead" and "phase lag", thus solving the control problem of time-delay systems.TD is a nonlinear tracking algorithm which can extract differential signals from target signals. Compared with the pure difference calculation, it can overcome the noise amplification effect to some extent. TD is often used for reference transition process generation and signal filtering. In the literature [8], aiming at the phenomenon that the input noise signal of the TD is amplified and causes the oscillation of the tracking signal, a design method of the phase lead compensator is proposed: the differential signal generated by the TD is connected in series with the TD to eliminate the oscillation of the output signal. Aiming at the problems of the traditional differentiator with the flutter phenomenon, slow dynamic response, and poor filtering ability, the literature [9] proposed an improved nonlinear TD, which takes into account the requirements of rapidity and accuracy, and realizes the signal tracking. For the application research of TD, the most prominent ones are spacecraft [10], unmanned aerial vehicles [11] and maglev trains [12,13], all of which have achieved very good results. Reference [14] proposed a nonlinear posterior linearization filter method robust to outliers for the nonlinear system filtering problem with abnormal observations. This paper firstly introduces the tracking differentiators in active disturbance rejection control algorithms. Compared with the original tracking differentiator, a speed and acceleration limiting mechanism whose threshold in two directions can be adjusted separately is added in the new asymmetric tracking differentiator (ATD). This speed limiting mechanism makes the ATD suitable for applications including the aircraft attitude control system. In such systems, limitation on rotational speed is mandatory. The asymmetric acceleration adjustment capability enables the ATD to be used in applications including attitude control, where flexible load adjustment capability shows important value. The proposed method is also verified by mathematical simulations and compared with the PID attitude control method.

Asymmetric Tracking Differentiator
Tracking-differentiator (TD) is an important part of the system of active disturbance rejection control technology developed by Professor Han J.Q. [15][16][17][18][19] at the end of the last century. Tracking differentiator is a nonlinear tracking algorithm which can extract differential signals from target signals. Compared with the pure difference calculation, it can overcome the noise amplification effect to some extent. Tracking differentiators are often used for reference transition process generation and signal filtering [20,21]. Because of its ability to generate the transition process without overshoot pending on given maximum required acceleration, it is especially suitable for aircraft guidance mechanism with strict requirements on available overload [22,23]. In this section, the basic tracking differentiator is introduced concisely. Then, according to the characteristics and usual mission requirements of a small unmanned helicopter, the asymmetric tracking differentiator suitable for the guidance mechanism of an unmanned helicopter is proposed.

Basic Tracking Differentiator
A basic discrete tracking differentiator is shown in Equation (1): where k is sampling time, ∆t is sampling interval; r is the tracked signal; x is the output signal of TD; v and ζ are respectively the first-order difference (velocity) and second-order difference (acceleration) of the relative time of x; a lmt is the available acceleration of TD, h is the account step of the function f han, and f han is the comprehensive function of discrete time-optimal control. Assume that the reference signal is r(t) = 1. The initial conditions are x(0) = 0, v(0) = 0. The available acceleration is set as a lmt = 1, the step length of TD is set as h = 5 × 10 −3 s, and the computer simulation of which the length is 3 s and the simulation step is 2.5 × 10 −3 s is carried out. The simulation curve is obtained as shown in Figure 1. system. In such systems, limitation on rotational speed is mandatory. The asymmetric acceleration adjustment capability enables the ATD to be used in applications including attitude control, where flexible load adjustment capability shows important value. The proposed method is also verified by mathematical simulations and compared with the PID attitude control method.

Asymmetric Tracking Differentiator
Tracking-differentiator (TD) is an important part of the system of active disturbance rejection control technology developed by Professor Han J.Q. [15][16][17][18][19] at the end of the last century. Tracking differentiator is a nonlinear tracking algorithm which can extract differential signals from target signals. Compared with the pure difference calculation, it can overcome the noise amplification effect to some extent. Tracking differentiators are often used for reference transition process generation and signal filtering [20,21]. Because of its ability to generate the transition process without overshoot pending on given maximum required acceleration, it is especially suitable for aircraft guidance mechanism with strict requirements on available overload [22,23]. In this section, the basic tracking differentiator is introduced concisely. Then, according to the characteristics and usual mission requirements of a small unmanned helicopter, the asymmetric tracking differentiator suitable for the guidance mechanism of an unmanned helicopter is proposed.

Basic Tracking Differentiator
A basic discrete tracking differentiator is shown in Equation (1): where k is sampling time, t  is sampling interval; r is the tracked signal; x is the output signal of TD; v and  are respectively the first-order difference (velocity) and second-order difference (acceleration) of the relative time of x ; lmt a is the available acceleration of TD, h is the account step of the function fhan , and fhan is the comprehensive function of discrete time-optimal control.
Assume that the reference signal is The initial conditions are (0) 0, (0) 0 xv == . The available acceleration is set as 1 lmt a = , the step length of TD is set as  As can be seen from the figure, TD always accelerates or decelerates at the maximum available acceleration in the process of tracking the reference signals. TD is a signal tracking algorithm with excellent performance because the curve () xt generated by TD has As can be seen from the figure, TD always accelerates or decelerates at the maximum available acceleration in the process of tracking the reference signals. TD is a signal tracking algorithm with excellent performance because the curve x(t) generated by TD has second derivative and no overshoot for the step signal. When TD is applied to the attitude control system, the attitude command generated by the guidance system is the tracked signal to generate the reference trajectory including attitude Angle, attitude angular rate and attitude angular acceleration, which is used for error calculation and feedback control. When TD is applied to the guidance system, the distance between the helicopter and the target navigation point or the cross-track-error between the target routes is treated as the tracking variable, and 0 is used as the reference signal. As a result, a smooth reference transition process without overshoot can be generated, and reference speed and reference acceleration are available. This is exceedingly beneficial to improve the bandwidth of the aircraft attitude control system and guidance system and improve the accuracy of attitude control and route tracking. This helps to ensure safety during landing and flying when there are high obstacles (trees, buildings, ships, etc.) near landing platforms, target points, or flight paths.
Although TD can provide a tracking curve with outstanding performance, it can be found from algorithm description and simulation curve that TD has the following two defects: (1) TD has no mechanism to limit the change speed v of the tracking signal. Small unmanned helicopters are mostly equipped with low-cost MEMS gyroscopes with limited range. Meanwhile, the small unmanned helicopter has tiny body size, high control moment-rotational inertia ratio and strong attitude maneuverability. If TD is used directly in the attitude control system, the angular velocity of an aircraft may exceed the maximum range of on-board gyroscope. In the process of autonomous flight, there is usually a maximum velocity limit or the requirement of flying at a given flight speed. High maneuvers are often accompanied by a limit on the maximum rate of climb and rate of descent. (2) TD has only one parameter a lmt to limit acceleration. A small unmanned helicopter has strong attitude maneuverability. Sometimes the helicopter is required to start and brake with different accelerations on a certain trajectory segment, and the aircraft is also required to decelerate with a more moderate trajectory to reduce the possibility of collision and shock when approaching the ground or landing platforms. However, it is difficult to deal with the above different cases simultaneously with a single parameter.
TD is improved on these two points by adding a mechanism for setting speed and acceleration limits in two directions.

Asymmetrical Tracking Differentiator
In light of the shortcomings of the TD algorithm proposed above, the maximum available acceleration and maximum available deceleration limits are set on the TD, so that it has the ability to limit the maximum change rate of the tracking signal.
Firstly, the synthesis function of the asymmetric discrete time-optimal control is derived to calculate the required acceleration. One then introduces the maximum and minimum available accelerations: l aup > 0 and l alw < 0, respectively. The difference equation of the second-order linear discrete system of the object is: where the control quantity u(k) satisfies: The system state at the initial moment is denoted as [x 1 (0), x 2 (0)] T , then the solution of Equation (2) is: If [x 1 (k), x 2 (k)] T = [0, 0] T , the reverse control sequence u(k) ∈ l ahv, l aup , in Equation (4) can make the system (2) reach the origin [0, 0] T from the initial value [x 1 (0), x 2 (0)] T within k steps. The phase plane OX 1 X 2 is introduced, and the horizontal and vertical axes correspond to the system state variables x 1 and x 2 , respectively. When x 2 < 0 and the starting point that exactly reaches the origin of the phase plane (0,0) after continuously exerting control u = l aup through k step is denoted as a + k . Coordinate values of a + k satisfy: when x 2 < 0, at the first step u = l alw is taken, and during the following (k − 1) step u = l aup is taken, and the starting point that exactly reaches the origin after the total k steps is denoted as b + k . The coordinate value of b + k satisfies: when x 2 < 0, u = 0 is taken at the first step, and u = l aup is taken during the next (k − 1) step, and the starting point that exactly reaches the origin after the total k steps is denoted as c + k . Coordinate values of point c + k satisfies: By eliminating k from Equations (5)- (7), it can be known that a + k , b + k and c + k are respectively on the part of the curve x 2 < 0, which is described by Equations (8)-(10).
Relatively, when x 2 > 0, a − k is denoted as the starting point that exactly reaches the origin after applying control quantity u = l alw continuously through k steps. b − k is denoted as the starting point that exactly reaches the origin under the circumstances that the control quantity of the first step is u = l aup , and the control quantity of other (k − 1) steps is set as u = l alw . c − k is denoted as the starting point that exactly reaches the origin under the circumstances that the control quantity of the first step is u = 0 and the control quantity of other (k − 1) steps is set as u = l alw . Then, a − k , b − k and c − k are located on the part of the curve x 2 > 0 described by Equations (11)-(13), respectively. See Figure 2.
Firstly, take the case where the minimum number of steps is 2 or less into consideration. If it is true, the initial state value of the system [ 1 (0), 2 (0)] is within the quadrilateral (including the boundary) formed by points 2 + , 2 + , 2 − , 2 − . The kinetic equation of the system reaching the origin within two steps is: when h > 0, directly work out the control quantity u(0) by using [ 1 (0), 2 (0)] as a known condition: Remove the label "(0)", the obtained u = −( 1 + 2ℎ 2 )/ℎ 2 is the control quantity when the system state [ 1 , 2 ] is within the quadrilateral 2 − (including the boundary). To determine 2 − , the coordinates of four points 2 + , 2 + , 2 − , 2 − should be determined at first. According to Equations (8), (9), (11) and (12), through the substitution of k = 2, the coordinate of each point is obtained as follows:  Firstly, take the case where the minimum number of steps is 2 or less into consideration. If it is true, the initial state value of the system [x 1 (0), x 2 (0)] is within the quadrilateral (including the boundary) formed by points a + 2 , b + 2 , a − 2 , b − 2 . The kinetic equation of the system reaching the origin within two steps is: when h > 0, directly work out the control quantity u(0) by using [x 1 (0), x 2 (0)] T as a known condition: Remove the label "(0)", the obtained u = −(x 1 + 2hx 2 )/h 2 is the control quantity when the system state [x 1 , should be determined at first. According to Equations (8), (9), (11) and (12), through the substitution of k = 2, the coordinate of each point is obtained as follows: According to the coordinates of the four points, the equations of the four lines at which the four sides of the quadrilateral a + 2 b + 2 a − 2 b − 2 are respectively located can be solved. The equations of the lines where the four lines a (17) to (20), respectively. See Figure 3.
According to the coordinates of the four points, the equations of the four lines at which the four sides of the quadrilateral 2  (17) to (20), respectively. See Figure 3. According to the above analysis, when the system state variable [ 1 , 2 ] meets the two condition 1 + ℎ 2 ∈ [ℎ 2 , ℎ 2 ], 1 + 2ℎ 2 ∈ [−ℎ 2 , −ℎ 2 ] the corresponding control quantity is: Then consider the case where the minimum number of steps is larger than two. If the initial point is on Γ + , the state variable [ 1 , 2 ] can return to the origin at the fastest speed by applying control quantity = continuously. If the initial point is on Γ + , then the control quantity = applied in the first step causes the state variable According to the above analysis, when the system state variable [x 1 , x 2 ] T meets the two condition x 1 + hx 2 ∈ h 2 l alw , h 2 l aup , x 1 + 2hx 2 ∈ −h 2 l aup , −h 2 l alw the corresponding control quantity is: Then consider the case where the minimum number of steps is larger than two. If the initial point is on Γ a+ , the state variable [x 1 , x 2 ] T can return to the origin at the fastest speed by applying control quantity u = l aup continuously. If the initial point is on Γ b+ , then the control quantity u = l awl applied in the first step causes the state variable [x 1 , x 2 ] T to fall on Γ a+ after the first step. Then the control quantity u = l aup can be applied continuously to make the state variable [x 1 , x 2 ] T return to the origin during the following (k − 1) steps. For example, the initial point of the system state is b + 5 , then it falls to the point a + 4 after the first step and then returns to the origin along the Γ a+ . If the initial point is on Γ c+ , then the control quantity u=0 applied in the first step causes the state variable [x 1 , x 2 ] T to fall on Γ a+ after one step. Then the control quantity u = l aup can be applied continuously to make the state variable [x 1 , x 2 ] T return to the origin during the following (k − 1) steps. For example, the initial point of the system state is c + 5 , then it falls to the point a + 4 after the first step and then returns to the origin along the Γ a+ . For the case where the initial point is located at Γ (·)− , there are similar conclusions, which will not be discussed repeatedly. Obviously, if the initial point is located in the region below the polyline Γ a+ Γ b− , the control quantity of the first step should be u = l aup . If the initial point is located in the region above the polyline Γ a− Γ b+ , the control quantity of the first step should be u = l alw . See following (k − 1) steps. For example, the initial point of the system state is 5 + , then it falls to the point 4 + after the first step and then returns to the origin along the Γ + . For the case where the initial point is located at Γ (·)− , there are similar conclusions, which will not be discussed repeatedly. Obviously, if the initial point is located in the region below the polyline Γ + Γ − , the control quantity of the first step should be = . If the initial point is located in the region above the polyline Γ − Γ + , the control quantity of the first step should be = . See Figure 4. If the shortest distance from the point P to the polyline Γ is denoted as dist(p,Γ), when the initial point is located at the region between the two polylines Γ + Γ − , Γ − Γ + the control quantity should have the following characteristics.
when the state point is located in the region between the curve Γ + Γ − and Γ − Γ + (as shown in Figure 4, interpolation is performed according to the transverse distance between the state point [ 1 , 2 ] and the curve Γ − Γ + to calculate the control quantity u. ̃1(Γ, 2 ) is the value of the abscissa 1 at the corresponding ordinate 2 on the curve Γ, then the calculation formula of u is as follows: If the shortest distance from the point P to the polyline Γ is denoted as dist(p,Γ), when the initial point is located at the region between the two polylines Γ a+ Γ b− , Γ a− Γ b+ the control quantity should have the following characteristics.
when the state point is located in the region between the curve Γ a+ Γ b− and Γ a− Γ b+ (as shown in Figure 4, interpolation is performed according to the transverse distance between the state point [x 1 , x 2 ] T and the curve Γ c− Γ c+ to calculate the control quantity u.
is the value of the abscissa x 1 at the corresponding ordinate x 2 on the curve Γ, then the calculation formula of u is as follows: According to Equations (5)- (7), it can be seen that point c Based on the above content, the partition boundaries of Figure 4) and specific expressions can be obtained: A saturation function sat(x,δ) is introduced, as shown in Algorithm 1: Based on the above content, the improved asymmetric discrete time−optimal control synthesis function fm is built as shown in Algorithm 2. Introduce variables l vup and l vlw , which are the maximum and minimum speed limits, respectively. Then the improved asymmetric tracking-differentiator (ATD) with the function of limiting the speed of the tracked signal and respectively setting the maximum and minimum of the acceleration is shown in Algorithm 2.
The reference signal r(t) is defined as follows: At the initial moment, x(0) = 0, v(0) = 0 the upper limit of available acceleration is l aup = 3 and the lower limit of available acceleration is l alw = −1. The speed limits are l vup = −0.5 and l vlw = −1, respectively. Computer simulation with a length of 6 s and simulation step of 2.5 × 10 −3 s is carried out, and a simulation curve was obtained as shown in Figure 5.
At the initial moment, x(0) = 0, v(0) = 0 the upper limit of available acceleration is = 3 and the lower limit of available acceleration is = −1. The speed limits are = −0.5 and = −1, respectively. Computer simulation with a length of 6 s and simulation step of 2.5 × 10 −3 s is carried out, and a simulation curve was obtained as shown in Figure 5. Based on the algorithm and simulation curves ( Figures 1 and 5), it can be seen that ATD can limit the change rate v and acceleration ζ of the tracked signal with different limits under the condition that the tracked signal x has the second derivative and there is no overshoot in the step response. ATD is introduced into the attitude control system, and are set according to the range of on−board gyroscope and , are set according to the attitude maneuverability of the body. When the ATD is introduced into the helicopter guidance system to generate the reference trajectory of plane motion, can be set according to the available overload of the helicopter, can be set according to the requirements of the deceleration process. Set the value of and according to the flight speed requirements, realizing the control and limitation of the flight speed. When ATD is used to generate the vertical reference track, set according to the limit condition of the maximum climb rate and set according to the limit condition of the maximum descent rate to realize the limit of the motion speed in the height direction. Based on the data set above and Based on the algorithm and simulation curves ( Figures 1 and 5), it can be seen that ATD can limit the change rate v and acceleration ζ of the tracked signal with different limits under the condition that the tracked signal x has the second derivative and there is no overshoot in the step response. ATD is introduced into the attitude control system, l vlw and l vup are set according to the range of on-board gyroscope and l aup , l alw are set according to the attitude maneuverability of the body. When the ATD is introduced into the helicopter guidance system to generate the reference trajectory of plane motion, l aup can be set according to the available overload of the helicopter, l alw can be set according to the requirements of the deceleration process. Set the value of l vup and l vlw according to the flight speed requirements, realizing the control and limitation of the flight speed. When ATD is used to generate the vertical reference track, set l vup according to the limit condition of the maximum climb rate and set l vlw according to the limit condition of the maximum descent rate to realize the limit of the motion speed in the height direction. Based on the data set above and

Attitude Control System ATD
The helicopter attitude control system generates main rotor loop screw pitch commands and tail rotor collective screw pitch commands according to the attitude commands given by the guidance system and attitude estimation value given by the attitude estimation system. According to the configuration characteristics, the attitude control system of the unmanned helicopter with a single main rotor is divided into two parts. One is the balance controller, which realizes the Angle/angular speed control of the roll and pitch direction by adjusting the loop screw pitch input of the main rotor. The other one is the tail-lock controller, also known as direction controller, which realizes the angular speed control along the yaw direction by adjusting the collective screw pitch input of the tail-rotor.
The attitude control is divided into two layers depending on different processed signals. As shown in the Figure 6 above, the first layer is a reference attitude angular speed generator containing roll and pitch directions. According to the vehicle attitude kinematics, reference attitude trajectory and reference angular speed commands are generated from reference roll Angle and reference pitch Angle. The triaxial angular speed controller on the second layer generates the control commands of three channels according to the reference attitude angular speed command and the output of on-board control gyroscope, which are then converted into the actuating gear commands by the hybrid controller. For small unmanned helicopters, the actuators are mainly driven by servo motor (steering gear), and the output of the attitude control system is rudder deflection commands.

Attitude Control System ATD
The helicopter attitude control system generates main rotor loop screw pitch commands and tail rotor collective screw pitch commands according to the attitude commands given by the guidance system and attitude estimation value given by the attitude estimation system. According to the configuration characteristics, the attitude control system of the unmanned helicopter with a single main rotor is divided into two parts. One is the balance controller, which realizes the Angle/angular speed control of the roll and pitch direction by adjusting the loop screw pitch input of the main rotor. The other one is the tail-lock controller, also known as direction controller, which realizes the angular speed control along the yaw direction by adjusting the collective screw pitch input of the tailrotor.
The attitude control is divided into two layers depending on different processed signals. As shown in the Figure 6 above, the first layer is a reference attitude angular speed generator containing roll and pitch directions. According to the vehicle attitude kinematics, reference attitude trajectory and reference angular speed commands are generated from reference roll Angle and reference pitch Angle. The triaxial angular speed controller on the second layer generates the control commands of three channels according to the reference attitude angular speed command and the output of on-board control gyroscope, which are then converted into the actuating gear commands by the hybrid controller. For small unmanned helicopters, the actuators are mainly driven by servo motor (steering gear), and the output of the attitude control system is rudder deflection commands.

ATD-Based Attitude Transition Process Design
Roll and pitch channels are discussed firstly. Two asymmetric tracking differentiators ATD ϕ and ATD θ are used to track the command attitude angles and from the guidance system ϕ gd and θ gd respectively, to generate roll and pitch reference transition trajectory including reference attitude Angle ϕ r , θ r , reference attitude angular velocity   Take r ATD ϕ (k) = ϕ r (k) as the input to the ATD ϕ , and r ATD θ (k) = θ r (k) as the input to ATD. The output of ATD (·) can be used to generate reference trajectory (·) T , .  As is shown in Equations (30) to (32): . ..
The attitude Angle change rate of fuselage is calculated according to the attitude value ϕ, θ, ψ given by the attitude estimation system, and the angular velocity of the body measured by the airborne gyroscope ω x , ω y , ω z . According to the attitude kinematics equations, the relationship between the angular velocity of fuselage and the attitude Angle change rate is shown in Equation (33): The calculation formula of attitude angular velocity can be obtained through inversion of the conversion matrix on the right side of Equation (33): The equation is singular at the pitch Angle θ = ±π/2. It is rare during the usual helicopter flight, and generally the pitch Angle |θ| < π/4 during unmanned helicopter autonomous flight. When θ = ±π/2 occurs, cos(θ) is substituted with a value whose absolute value is tiny for the calculation. The deductions in the following passage are based on assumptions θ = ±π/2.
In order to realize attitude control, PI feedback is used for roll Angle error and pitch .  .. ..
For the yaw channel, the guidance system transmits the yaw angular speed instruction . ψ gd . Due to the limitation of GNSS update frequency and the computing power of the airborne computer, the operational frequency of the guidance system is generally much lower than that of the attitude control system. In order to match signal sampling rate, a third differentiator ATD ψ can be used for tracking According to Equation (33), calculate the angular velocity of the reference fuselage: For occasions where the signal quality is good or there are higher requirements on the dynamic characteristics of the aircraft, the reference angular acceleration can be calculated according to the above formula: . .
Equations (39)-(42) do not contain ψ, which ensures that the reference angular velocity generator is free from the disturbance from the yaw Angle estimation error. When the attitude estimation algorithm only uses magnetometer for directional correction, it is affected by magnetic field interference or the dual-antenna GNSS measurement system is out of lock, even when the directional estimation value has a large error, the reference attitude trajectory generator itself will not be affected.
Then the influence of roll Angle and pitch Angle estimation errors generated by the attitude estimator on the reference body's angular velocity commands are discussed. Assume the real roll Angle and pitch Angle of the helicopter are ϕ and θ. The corresponding attitude estimation errors are ε ϕ and ε θ . The estimated value of the roll Angle and pitch Angle can be expressed as ϕ = ϕ + ε ϕ , θ = θ + ε θ .
The trim pitch Angle of the unmanned helicopter with the adjusted center of mass is almost zero when hovering, and the trim roll Angle is generally less than 5 • . In the hovering state with a small attitude estimation error, Equation (39) can be approximately expressed as follows: It can be seen that under such conditions, the pitch estimation error will cause a disturbance term which is proportionate to command yaw rate in the command angular rate of O b X b channel. The roll Angle estimation error will lead to a command angular velocity of O b Y b axis coupling with that of O b Z b axis. When an unmanned helicopter performs a fixed-point hover mission, generally . ϕ m = 0, . θ m = 0, but it is also possible to make direction adjustments (e.g., adjust the load vector) to ensure . ψ m = 0. In hovering state, the disturbance generated by attitude estimation error on the reference angular velocity command can be reduced by lowering the speed limit of ATD ψ .
When flying forward, the helicopter generally has a large pitch Angle, but the roll Angle remains a little more or less than the trim pitch Angle. When roll estimation error ε ϕ is small, Equation (39) can be approximately regarded as: For the attitude estimation method described above, the maneuver overload corresponding to the helicopter speed change is the main source of error. According to the above analysis, the coupling between the reference yaw rate and the reference angular velocity of O b X b axis will be strengthened during the starting acceleration and braking deceleration, and ε θ is in the same sign with θ. For the application of mobile platform autonomous takeoff and landing, direction adjustment should be mainly completed in the approaching stage for the sake of tracking accuracy. There should not be large yaw rate instruction in the following stage.

ATD-falESO Angular Rate Controller
According to the mathematical models in chapter 2, the relationship between the control torque generated by the main rotor and tail rotor and the control input is complex, and it is obviously affected by the flight state. However, the simplified force and torque expressions show that the transmission paths of the helicopter's basic control action in three channels of the corresponding axes of fuselage are similar. The angular speed control models upwards along three axes, namely O b X b , O b Y b and O b Z b , can be expressed as: .
where . ω is angular acceleration. b(t) is the control gain coefficient determined by rotor structure parameters, rotation rate and other working conditions. u(t) is the control input. The two channels O b X b and O b Y b are the main rotor loop pitch input; the channel of axis O b Z b is the collective pitch input of the tail rotor. The d(t) is the sum of interference effects, including airflow, geometric parameters simplification, external interference, and other factors.
In order to estimate these disturbances quickly and eliminate them by adjusting the control quantity u(t), a disturbance observer is installed in supplement to the error feedback mechanism. In order to adapt to the research content of this paper, the computation quantity of the selected disturbance observer, so as to realize the control algorithm on the airborne MCU with low computing power. In addition, the single-rotor helicopter with a tail rotor has many movable parts, a complex transmission mechanism and large vibration during operation. The selected disturbance observer should have strong robustness to the measurement noise. At the same time, the control actuator of a small unmanned helicopter is mostly driven by a servo motor. Because the output of the disturbance observer is directly used to calculate the control signal, in order to ensure the servo motors normally work for a long time, there should not be serious shake and vibration phenomenon when using the selected disturbance observer. The mathematical model of helicopter is complex and the disturbance sources are diverse, so the selected disturbance observer should not have strict requirements on the boundary and character of the disturbance value.
Through research and comparisons of different kinds of disturbance observers which are suitable for first-order uncertainty system, the results show that there are two kinds of disturbance controllers suitable for the research object in this paper: high power disturbance observer HPDOB and falESO. Because of the linear structure of the fal function near zero point, the observation results are less affected by measurement noise than HPDOB, hence more suitable for the research content of this paper. So falESO is chosen as the disturbance observer of the attitude control system.
The first-order system model is shown in Equation (46).
where x is the system state; f is the known system model; d is the uncertain disturbance including the unknown system model and external disturbance, and it is sectionally continuous and has the first derivative versus time; u is the control effect. Assume that the state measurement value containing noise v(t) is x = x + v. For this system, a basic falESO is as follows: where fal is a fractional power function: Design parameters α and δ meet α ∈ (0, 1), δ > 0.
The convergence of the falESO of this configuration is proved and the tracking accuracy is estimated through a literature review. According to the attitude angular velocity control model described by Equation (45), Equation (47) is rewritten into a discrete form which is suitable for digital computer execution, as follows: In view of the limited computing power of airborne MCU, in order to run the attitude control programs at a higher frequency, the calculation amount should be considered when setting falESO parameters. Selection of α = 1 2 can avoid calling POW function and it also takes computation and observer performance into consideration. By using fast square root algorithm or DSP instruction set can further reduce the amount of computation and improve the calculation speed.
In falESO, Z 2 is the interference estimation term. Add it into the PI error feedback mechanism to form disturbance rejection angular velocity controller ATD-falESO: The angular acceleration can be introduced into Equation (51) to further improve the response speed of the angular velocity controller when the reference angular acceleration signal is good. In order to avoid amplifying the noise caused by gyro output signal difference, only reference angular acceleration information is used to form the direct feedback term. As is shown in Equation (51).
After introducing the falESO, the value of control gain coefficient b will have a crucial influence on interference estimation and controller operation. The math models in chapter 2 show that the control force and torque produced by the rotor are directly related to the its rotational speed. With other conditions being certain, the resulting control force and torque are approximately proportional to the square of the rotational speed.
The previous study shows that the amplification gain coefficient of a small unmanned helicopter can be adjusted according to the rotational speed of the rotors. Assume that the nominal rotational speed is Ω 0 , and the nominal control gain coefficient calculated by experiment or simulation under nominal conditions is b 0 , then the current control gain coefficient can be estimated according to the following formula: Based on the previous research work, this paper uses Equation (53) as the mapping number to adjust the angular velocity control gain of an unmanned helicopter in real time during flight. In addition, K 0p and K 0i are respectively the nominal scale and integral gain coefficients measured through experiments or calculated through simulations under nominal rotational speed.
The forward feedback coefficient K f can be directly estimated according to the simplified mathematical model in Chapter 2. For the O b X b and O b Y b channels, the output of angular velocity controller is the main rotor loop pitch control quantity, and the forward feedback coefficient is estimated by Equations (54) and (55), respectively. For the O b Z b channel, the output of the angular velocity controller is the collective pitch control quantity of the tail rotor, and the forward feedback coefficient is estimated by formula (56), and the signs are determined according to the installation direction of the tail rotor. Where I XX , I yy and I zz are the main moments of inertia on the three axes of the fuselage respectively, and K AB are the ratio of the main rotor loop pitch input and the periodic flapping Angle along corresponding direction, which can be determined by simulation.

Fuselage Parameter Calculation of ATD Attitude Controller
Small unmanned helicopters generally have a higher load/weight ratio [24][25][26]. The payload quantity maximum of helicopters of 25 kg level and below exceeds 50% their takeoff weight maximum. If there is a load dropping or spreading task during the flight, the inertial parameters of the fuselage will have obvious changes as tasks progress. Many small unmanned helicopters are powered by electricity. The voltage and discharge capacity (the ability to provide current) of high−power battery change obviously during discharge [27,28]. Temperature variation also has a significant effect on the discharge capacity of the battery. In some tasks, even if there is a motor speed control device, the variation of the main rotor of the electric helicopter can reach 15% of the nominal speed.
The missions of an unmanned helicopter, including the autonomous take−off and landing process all have relatively high requirements on the precision of position control, attitude control and dynamic performance [29]. Angular rate control is the basis of attitude control and position control, and it is often necessary to achieve better dynamic performance indicator by the fine tuning of control parameters. In ATD attitude control system, the inertial parameters of the fuselage affect the forward feedback coefficient k f in Equation (51), while the control gain coefficient b directly affects the observation and elimination of disturbance (as shown in Equations (49) and (51)). If a mechanism can be established to estimate the inertial parameters and control gain coefficients b in each channel in real time during flight, it will be exceedingly helpful to improve the performance of the whole control system.
Based on the input-output relationship of the system and CAD data, this paper presents a method to estimate the control gain coefficient and inertial characteristics of the system, and to adjust the control parameters accordingly.

Estimation of Control Gain Coefficient
According to the first-order kinetics difference equation of the system, the control gain coefficient of the system is estimated by measuring the control input and the system output, and applying forgetting factor least squares estimation. Take the angular velocity control model of a single axis into consideration. Discretize Equation (45) and ignore the interference term to obtain the discretized angular velocity control model of a single axis: where ∆t is the sampling time interval, u is the control quantity, and ω is the angular velocity component of the helicopter on an axis of the system. Take O b Z b axis as an example, u is the output rudder deviation of the tail rotor, ω is the angular velocity on the fuselage O b X b shaft. According to the general form of first-order discrete system, the estimation model is defined as follows: Among them y(k) = ω(k), β 0 = ∆t·b. Defines the estimated parameter matrix θ = α 1 ,β 0 T , whereα 1 ,β 0 are the estimates of α 1 and β 0 respectively. The matrix of measurement data is ∅(k) = [y(k − 1), u(k − 1)] T . Forgetting factor least squares method was used to obtain the following recursive estimation formula: In the initial moment, k = 0, initialize P with an initial value of e i I, where e i is a positive number much less than one. Take the estimation of control gain coefficient of a small unmanned helicopter O b Z b shaft. Angular rate signal and steering engine output signal measured in a test are shown in Figure 7, and the sampling frequency is 80 Hz.
where ∆ is the sampling time interval, u is the control quantity, and ω is the angular velocity component of the helicopter on an axis of the system. Take axis as an example, u is the output rudder deviation of the tail rotor, ω is the angular velocity on the fuselage shaft. According to the general form of first-order discrete system, the estimation model is defined as follows: Among them y(k)= ( ), 0 = ∆ • . Defines the estimated parameter matrix ̂= [̂1,̂0] , where ̂1, ̂0 are the estimates of 1 and 0 respectively. The matrix of measurement data is ∅( ) = [ ( − 1), ( − 1)] . Forgetting factor least squares method was used to obtain the following recursive estimation formula: In the initial moment, k = 0, initialize P with an initial value of , where is a positive number much less than one. Take the estimation of control gain coefficient of a small unmanned helicopter shaft. Angular rate signal and steering engine output signal measured in a test are shown in Figure 7, and the sampling frequency is 80 Hz.  Calculated by Equation (59),β 0 and the main rotor speed Ω is shown in Figure 8. Calculated by Equation (59), ̂0 and the main rotor speed Ω is shown in Figure 8. The control gain coefficient value b and the main rotor speed value are calculated according to the value of ̂0 in the stable state, as is shown in the Table 1.
According to the mathematical model of the helicopter, the directional control moment along produced by the tail rotor is approximately proportional to the square of the tail rotor speed. In addition, because the main rotor and tail rotor are generally connected by transmission machinery with fixed reduction ratio, the least square fitting of the estimated value of b and the square of the main rotor speed value are carried out. The equation of the fitting line is shown in Equation (60). Where the units of b and Ω are respectively s −2 and rad/s. The control gain coefficient value b and the main rotor speed value are calculated according to the value ofβ 0 in the stable state, as is shown in the Table 1. According to the mathematical model of the helicopter, the directional control moment along O b Z b produced by the tail rotor is approximately proportional to the square of the tail rotor speed. In addition, because the main rotor and tail rotor are generally connected by transmission machinery with fixed reduction ratio, the least square fitting of the estimated value of b and the square of the main rotor speed value are carried out. The equation of the fitting line is shown in Equation (60). Where the units of b and Ω are respectively s −2 and rad/s.
Data points and fitting lines are shown in Figure 9. The control gain coefficient value b and the main rotor speed value are calculated according to the value of ̂0 in the stable state, as is shown in the Table 1.
According to the mathematical model of the helicopter, the directional control moment along produced by the tail rotor is approximately proportional to the square of the tail rotor speed. In addition, because the main rotor and tail rotor are generally connected by transmission machinery with fixed reduction ratio, the least square fitting of the estimated value of b and the square of the main rotor speed value are carried out.  For all the data points shown in Table 1, the maximum deviation between the fitted value of b and its estimated value output by the estimation algorithm is about 22% of the fitted value. For all the data points shown in Table 1, the maximum deviation between the fitted value of b and its estimated value output by the estimation algorithm is about 22% of the fitted value.
The main rotor speed was selected as 1200 rpm (125.66 rad/s) and 2400 rpm (251.33 rad/s), respectively, to check the estimation formula (60). Under two working conditions, the variation of angular velocity output of the fuselage ω versus time is shown in Figure 10. The maximum reverse rudder deviation u = −1.18 rad corresponds to the rapid maneuvering process starting at 0.5 s. According to the measured data lines, the maximum acceleration of this part under two working conditions is respectively −20.21 rad/s 2 and −66.25 rad/s 2 . The corresponding control gain coefficient b under these two working conditions are obtained as 17.12 s −2 and 56.14 s −2 . The estimated values obtained from Equation (60) are respectively 15.16 s −2 and 60.64 s −2 . The deviation between the estimated value and the measured value is 13% and 7% of the estimated value respectively, and the match conditions are decent. celeration of this part under two working conditions is respectively −20.21 rad/s 2 and −66.25 rad/s 2 . The corresponding control gain coefficient b under these two working conditions are obtained as 17.12 s −2 and 56.14 s −2 . The estimated values obtained from Equation (60) are respectively 15.16 s −2 and 60.64 s −2 . The deviation between the estimated value and the measured value is 13% and 7% of the estimated value respectively, and the match conditions are decent.

Estimation of Inertial Parameters
Almost all modern small unmanned helicopters are designed and produced by computer aided design software and computer aided machining software. By virtue of the rich functions of CAD software and the powerful processing capacity of modern microcomputers, ordinary engineers can easily estimate the body inertial characteristics of unmanned helicopters, including mass, location of center of mass, moment of inertia and other parameters through electronic drawings. See Figure 11. Based on well−precision modeling and following industry standard databases of standard parts and materials, the mass, main moment of inertia and geometric feature estimation errors of small unmanned helicopters can easily be controlled within 10%.

Estimation of Inertial Parameters
Almost all modern small unmanned helicopters are designed and produced by computer aided design software and computer aided machining software. By virtue of the rich functions of CAD software and the powerful processing capacity of modern microcomputers, ordinary engineers can easily estimate the body inertial characteristics of unmanned helicopters, including mass, location of center of mass, moment of inertia and other parameters through electronic drawings. See Figure 11. Based on well−precision modeling and following industry standard databases of standard parts and materials, the mass, main moment of inertia and geometric feature estimation errors of small unmanned helicopters can easily be controlled within 10%. The change of inertial characteristics of an unmanned helicopter is mainly caused by load dropping, spreading, picking up and fuel consumption. Either of these conditions will be accompanied by changes in mass. In order to reduce the change of the location of center of mass, when designing the helicopters, usually the payload and fuel tank is placed near the location of center of mass or near the axis. Most of the flight time unmanned helicopters are required to maintain a given altitude. At this point, the balance between the lift force generated by the main rotor and the gravity has the following approximate relationship. Where ℎ is the mass of the body, and g is the gravitational acceleration.
It can be seen from the second chapter that, with other parameters fixed, the tension of the main rotor in the wind shaft system is approximately proportional to the collective pitch input 0 and the square of the main rotor speed Ω 2 . These relationships are the theoretical basis for estimating the change of collective moment of inertia from the The change of inertial characteristics of an unmanned helicopter is mainly caused by load dropping, spreading, picking up and fuel consumption. Either of these conditions will be accompanied by changes in mass. In order to reduce the change of the location of center of mass, when designing the helicopters, usually the payload and fuel tank is placed near the location of center of mass or near the O b Z b axis.
Most of the flight time unmanned helicopters are required to maintain a given altitude. At this point, the balance between the lift force generated by the main rotor and the gravity has the following approximate relationship. Where m veh is the mass of the body, and g is the gravitational acceleration.
It can be seen from the second chapter that, with other parameters fixed, the tension of the main rotor in the wind shaft system F MR w T is approximately proportional to the collective pitch input θ 0 and the square of the main rotor speed Ω 2 . These relationships are the theoretical basis for estimating the change of collective moment of inertia from the collective pitch input of the main rotor. Taking a 25 kg unmanned helicopter as an example, the relationship between the fuselage mass and the moment of inertia of each axis in the process of load consumption is firstly discussed, and then the relationship between the collective pitch input of the main rotor and the moment of inertia is established.
The maximum load capacity of the 25 kg unmanned helicopter is 10 kg, which is mainly liquid load. The load is placed in the load barrels on both sides of the fuselage and it is consumed with the spreading during the mission. The assembly drawing of the fuselage (including the load barrel, but excluding the main rotor and tail rotor) is shown in Figure 12.
placed near the location of center of mass or near the axis. Most of the flight time unmanned helicopters are required to maintain a given altitude. At this point, the balance between the lift force generated by the main rotor and the gravity has the following approximate relationship. Where ℎ is the mass of the body, and g is the gravitational acceleration.
It can be seen from the second chapter that, with other parameters fixed, the tension of the main rotor in the wind shaft system is approximately proportional to the collective pitch input 0 and the square of the main rotor speed Ω 2 . These relationships are the theoretical basis for estimating the change of collective moment of inertia from the collective pitch input of the main rotor. Taking a 25 kg unmanned helicopter as an example, the relationship between the fuselage mass and the moment of inertia of each axis in the process of load consumption is firstly discussed, and then the relationship between the collective pitch input of the main rotor and the moment of inertia is established.
The maximum load capacity of the 25 kg unmanned helicopter is 10 kg, which is mainly liquid load. The load is placed in the load barrels on both sides of the fuselage and it is consumed with the spreading during the mission. The assembly drawing of the fuselage (including the load barrel, but excluding the main rotor and tail rotor) is shown in Figure 12. Loads vary from 0% to 100% at an interval of 10%. For each load amount, the body mass m veh and the moments of inertia of the three axes, namely I XX , I yy and I zz , are given by CAD software, as shown in Table 2. It can be seen from Table 2 that I XX , I yy and I zz are in an approximately linear relationship with mass m veh . The least square fitting is carried out, and the equation of the fitting line obtained is shown as Equations (62)-(64). Where the unit of I XX , I yy and I zz is kg·m 2 , and the unit of m veh is kg. Data points and fitting lines are shown in Figure 13.
I yy = 6.332 × 10 −3 m veh + 1.890, (63) It can be seen from Table 2 that IXX, Iyy and Izz are in an approximately linear relationship with mass m veh . The least square fitting is carried out, and the equation of the fitting line obtained is shown as Equations (62)-(64). Where the unit of IXX, Iyy and Izz is kg·m 2 , and the unit of m veh is kg. Data points and fitting lines are shown in Figure 13. As can be seen from the figures above, for this type of helicopter, the effect of load changes are relatively obvious on IXX and Izz, and the most apparent on IXX . While the influence of load changes on Iyy is hardly seen, which is related to the general configuration of helicopters. When the load is hung on either side of the fuselage, the distance from the axis will be further than that of most parts and will account for a large part of the calculation of IXX. When calculating Iyy, the mass of the body tail and head account for the largest weight, especially the mass of the tail that expands backwards. The center of mass of the load is roughly located in the plane , and when the load changes, basically only the change of its own moment of inertia has an effect on Iyy. The fuselage Izz is roughly in the same situation as Iyy, basically the same order of magnitude. However, when the As can be seen from the figures above, for this type of helicopter, the effect of load changes are relatively obvious on I XX and I zz , and the most apparent on I XX . While the influence of load changes on I yy is hardly seen, which is related to the general configuration of helicopters. When the load is hung on either side of the fuselage, the distance from the O b X b axis will be further than that of most parts and will account for a large part of the calculation of I XX . When calculating I yy , the mass of the body tail and head account for the largest weight, especially the mass of the tail that expands backwards. The center of mass of the load is roughly located in the plane O b Y b Z b , and when the load changes, basically only the change of its own moment of inertia has an effect on I yy . The fuselage I zz is roughly in the same situation as I yy , basically the same order of magnitude. However, when the load is on both sides of the body, the center of mass is also in some distance from the axis O b Z b , roughly equal to the distance from the axis O b X b . In this case, the changes of I XX and I zz caused by load changes are roughly equivalent. However, compared with the other two axis components, I XX of the helicopter body with general configuration is the smallest, so load changes have the most significant impact on I XX .
A low−pass filter is introduced. When the target altitude is unchanged, the helicopter mass is estimated according to Equation (65) with the estimated attitude angles of the main rotor collective pitch, main rotor speed, estimated values of roll and pitch attitude Angles being input. The estimated quality of helicopter ism veh , and the calculation formula of the kth execution cycle is: where η ∈ (0, 1) is used to adjust the filter bandwidth, the coefficient K t can be measured by experiment or calculated by simulation. For the above example, when the main rotor speed is 141.37 rad/s (1350 rpm), the simulation results show the relationships and linear fitting lines between the body massm veh and the collective pitch input θ 0 under various load conditions when hovering, as shown in Figure 14, and K t = 9.81 × 10 −2 mkg·rad −3 is obtained through calculation.
whereη∈ (0,1) is used to adjust the filter bandwidth, the coefficient can be measured by experiment or calculated by simulation. For the above example, when the main rotor speed is 141.37 rad/s (1350 rpm), the simulation results show the relationships and linear fitting lines between the body mass ̂ℎ and the collective pitch input 0 under various load conditions when hovering, as shown in Figure 14, and = 9.81 × 10 −2 mkg·rad −3 is obtained through calculation. Combined with Equations (62)-(65), the estimated recursive formula of the rotational inertia of this type of 25 kg unmanned helicopter in hovering or stable flight is obtained as follows. Where the unit of rotational inertia is kg·m 2 , and the unit of each angular quantity is rad.

Simulation System Configuration
The simulation system is shown in Figure 15. As static variables, helicopter parameters, atmospheric and wind data, planetary model parameters are called by the simulation programs during run−time. The helicopter parameters include the mass, size, main rotor, tail rotor installation position and geometric parameters, rotational inertia matrix and so on. Based on the research content of Section 2, the helicopter dynamics module takes the output of the main rotor and tail rotor control actuators and the motion state of the helicopter as input, calls the helicopter parameter table and atmospheric data, and calculates the flapping motion of the rotor according to its flapping kinetics. The aerodynamic force on the helicopter is calculated according to the force and moment calculation formula. Its internal structure is shown in Figure 16. Firstly, coordinate transformation is carried out according to the body attitude parameters and rotor installation data, and the component matrix of the relative airflow velocity of main rotor and tail rotor in the hub coordinate system is obtained. Then, according to the flapping motion equation of the main rotor, Equation (8), take the flow velocity and control pitch in the hub coordinate system as input, and work out the flapping motion of the main rotor. According to Equations (19)- (24) and (37)-(40), the forces and torques generated by the main rotor and tail rotor are calculated and output after synthesis. According to the rigid body dynamics and kinematics equations, the "real" motion state of the helicopter relative to the reference coordinate system O I X I Y I Z I was calculated by taking the force and moment output from the helicopter dynamics module as input. According to the measurement of "real" motion state by sensor and GNSS module, the attitude estimation module uses the CGD attitude estimation algorithm in Section 3 to estimate the attitude of helicopter. The autopilot takes GNSS module measurement, pre−stored waypoint information and target state as input, and calculates the reference attitude Angle and angular velocity by using the guidance algorithm described in Section 4. The attitude and angular rate controller take the measurement of airborne motion sensor, the attitude estimation value of the fuselage output by the attitude estimation module and the reference attitude Angle and angular velocity output by the autopilot as input, and generates the actuator instructions according to the attitude control method discussed in Section 5. The actuator generates the control pitch of the main rotor and tail rotor according to the commands, and outputs to the helicopter dynamics module to form a closed loop.
Equation (8), take the flow velocity and control pitch in the hub coordinate system as in-put, and work out the flapping motion of the main rotor. According to Equations (19)- (24) and (37)-(40), the forces and torques generated by the main rotor and tail rotor are calculated and output after synthesis. According to the rigid body dynamics and kinematics equations, the "real" motion state of the helicopter relative to the reference coordinate system was calculated by taking the force and moment output from the helicopter dynamics module as input. According to the measurement of "real" motion state by sensor and GNSS module, the attitude estimation module uses the CGD attitude estimation algorithm in Section 3 to estimate the attitude of helicopter. The autopilot takes GNSS module measurement, pre−stored waypoint information and target state as input, and calculates the reference attitude Angle and angular velocity by using the guidance algorithm described in Section 4. The attitude and angular rate controller take the measurement of airborne motion sensor, the attitude estimation value of the fuselage output by the attitude estimation module and the reference attitude Angle and angular velocity output by the autopilot as input, and generates the actuator instructions according to the attitude control method discussed in Section 5. The actuator generates the control pitch of the main rotor and tail rotor according to the commands, and outputs to the helicopter dynamics module to form a closed loop.  All algorithms are programmed with C language, and the same code is used for mathematical simulation and computer control system. For mathematical simulation, C-MEX files are generated by a VSC/C++ compiler and run in a Matlab/Simulink. The same code and the Arm embedded development tool chain are used to generate Flash images which are suitable for MCU in the airborne controller for actual flight or HIL simulation. See Figure 17. All algorithms are programmed with C language, and the same code is used for mathematical simulation and computer control system. For mathematical simulation, C-MEX files are generated by a VSC/C++ compiler and run in a Matlab/Simulink. The same code and the Arm embedded development tool chain are used to generate Flash images which are suitable for MCU in the airborne controller for actual flight or HIL simulation. See Figure 17.
All algorithms are programmed with C language, and the same code is used for mathematical simulation and computer control system. For mathematical simulation, C-MEX files are generated by a VSC/C++ compiler and run in a Matlab/Simulink. The same code and the Arm embedded development tool chain are used to generate Flash images which are suitable for MCU in the airborne controller for actual flight or HIL simulation. See Figure 17.  Table 3 shows the main rotor parameters for simulation, and the tail rotor parameters are shown in Table 4. The nominal inertia parameters of fuselage are shown in Table 5.    Table 3 shows the main rotor parameters for simulation, and the tail rotor parameters are shown in Table 4. The nominal inertia parameters of fuselage are shown in Table 5. Position array of hub center of main rotor in fuselage coordinate system is P MR b = [0, 0, −0.238] T (m), and installation Angle is i s = 0. Position array of hub center of tail rotor in fuselage coordinate system is P MR b = [−1.365, 0, 0.070] T (m).

Roll Control Simulation
For the control simulation of roll Angle, set the control parameters K Pϕ = 6.0, K 1ϕ = 0, K p = 0.3, K i = 0.8. The longitudinal loop pitch limit is ±0.2 rad. Reference roll Angle ϕ r is shown as follows: The curves of attitude Angle ϕ, angular velocity ω x and loop pitch control quantity A 1s which are corresponding to ordinary PID, vs. time controller are shown in Figure 18.

Roll Control Simulation
For the control simulation of roll Angle, set the control parameters = 6.0, 1 = 0, = 0.3, = 0.8. The longitudinal loop pitch limit is ±0.2rad. Reference roll Angl is shown as follows: The curves of attitude Angle , angular velocity and loop pitch control quantit 1 which are corresponding to ordinary PID, vs. time controller are shown in Figure 18 (a) (b) (c) By comparison with Figure 18a,c, it can be seen that after the step change of referenc attitude Angle , the output of the control quantity 1 maintains the maximum ±0.2ra in a short time because the feedback gain coefficient is large. When approaching th reference value 0.3 rad/s, then decelerate with 1 = −0.2 rad Due to the saturation lim By comparison with Figure 18a,c, it can be seen that after the step change of reference attitude Angle ϕ r , the output of the control quantity A 1s maintains the maximum ±0.2 rad in a short time because the feedback gain coefficient is large. When ϕ approaching the reference value 0.3 rad/s, then decelerate with A 1s = −0.2 rad Due to the saturation limit of actuator, the energy accumulated during acceleration motion of helicopter attitude cannot be dissipated in a short time, resulting in obvious overshoot. This phenomenon is more obvious when t = 2 s, the step amplitude of the signal is doubled. As can be seen from Figure 18b, the control quantity jumps between the maximum and minimum values, resulting in significant changes in the angular velocity along O b X b axis of the body. The absolute value of the maximum angular velocity is greater than 6 rad/s, which requires a large range of airborne gyroscope to realize a safe control process.
ATD-falESO attitude controller is constructed by mapping ATD reference trajectory generator and falESO disturbance observer with constant feedback gain. For ATD, set the speed limit as l vup = 3 rad/s, l vlw = −3 rad/s, and acceleration limit is l aup = 10 rad/s 2 , l alw = −10 rad/s 2 . The control gain coefficient of falESO is set as b = 170 according to the simulation results. The generated change process of roll Angle ϕ, angular velocity ω x and control quantity A 1s are shown in Figure 19.
ATD-falESO attitude controller is constructed by mapping ATD reference trajectory generator and falESO disturbance observer with constant feedback gain. For ATD, set the speed limit as = 3 / , = −3 / , and acceleration limit is = 10 / 2 , = −10 / 2 . The control gain coefficient of falESO is set as b=170 according to the simulation results. The generated change process of roll Angle , angular velocity and control quantity 1 are shown in Figure 19. Due to the function of ATD, the underlying angular velocity controller can obtain rich reference trajectory information including reference angular velocity and angular acceleration, which increases the bandwidth of the system. It can be seen from Figure 19a,b that the roll Angle has definite acceleration and deceleration stages during the transition process. The transition process of attitude Angle is smooth and the overshoot is small. Although its rising speed is smaller than that of PID controller, due to the small overshoot, there is almost no oscillation, and the actual transition process time is about 50% of that of PID controller. By contrast with Figure 19b, it can be seen that the angular velocity is controlled within ±3rad/s and there is little redundant adjustment. Compared with PID controller, gyroscope with smaller range can be used to achieve the control purpose. Figure 19c shows that the output control quantity of ATD-falESO controller 1 does not reach saturation limit in the whole simulation process, which shows an obvious improvement in efficiency when compared with PID controller.

Pitch Control Simulation
For the pitching attitude control, the control parameters are = 6.0, = 0, = 0.8 and = 2.0. The lateral loop pitch limit is ±0.2rad. And the equation for reference pitch Angle is the same with (69). Due to the large rotational inertia of the pitching axis, the response of the pitching axis is "slower" when the available loop pitch limit is the same. During the step reference pitch Angle tracking, the phenomenon of overshoot and oscillation is obviously weakened. However, the error accumulated in the integral in the ascending process is difficult to eliminate in a short time, resulting in the extension of the transition process. See Figure   Figure 19. Due to the function of ATD, the underlying angular velocity controller can obtain rich reference trajectory information including reference angular velocity and angular acceleration, which increases the bandwidth of the system. It can be seen from Figure 19a,b that the roll Angle ϕ has definite acceleration and deceleration stages during the transition process. The transition process of attitude Angle is smooth and the overshoot is small. Although its rising speed is smaller than that of PID controller, due to the small overshoot, there is almost no oscillation, and the actual transition process time is about 50% of that of PID controller. By contrast with Figure 19b, it can be seen that the angular velocity is controlled within ±3 rad/s and there is little redundant adjustment. Compared with PID controller, gyroscope with smaller range can be used to achieve the control purpose. Figure 19c shows that the output control quantity of ATD-falESO controller A 1s does not reach saturation limit in the whole simulation process, which shows an obvious improvement in efficiency when compared with PID controller.

Pitch Control Simulation
For the pitching attitude control, the control parameters are K P θ = 6.0, K I θ = 0, K P = 0.8 and K i = 2.0. The lateral loop pitch limit is ±0.2 rad. And the equation for reference pitch Angle θ r is the same with (69).
Due to the large rotational inertia of the pitching axis, the response of the pitching axis is "slower" when the available loop pitch limit is the same. During the step reference pitch Angle θ r tracking, the phenomenon of overshoot and oscillation is obviously weakened. However, the error accumulated in the integral in the ascending process is difficult to eliminate in a short time, resulting in the extension of the transition process. See ADT reference trajectory generator and falESO disturbance observer are introduced under the condition that the feedback gain is fixed. Setting parameters of ATD are the ADT reference trajectory generator and falESO disturbance observer are introduced under the condition that the feedback gain is fixed. Setting parameters of ATD are the same as roll axis, and the control gain coefficient of falESO as b = 50, the pitch Angle θ, angular velocity of pitch axis ω y and lateral loop pitch B 1s are shown in Figure 21. ADT reference trajectory generator and falESO disturbance observer are introduced under the condition that the feedback gain is fixed. Setting parameters of ATD are the same as roll axis, and the control gain coefficient of falESO as b = 50, the pitch Angle θ, angular velocity of pitch axis and lateral loop pitch 1 are shown in Figure 21.
(a) (b) (c) According to Figures 20a and 21a, similar to the rolling direction, the control response is significantly improved because the reference transition process includes high−frequency information. In comparison to Figure 19a, it can be seen that the reference transition trajectory information generated under the same circumstances is also the same due to the identical ATD parameters. Therefore, the ATD-falESO controller produces very similar transition processes in the pitch and roll channels. In addition, because falESO disturbance observer has a faster response speed, the influence of accumulated errors on the transition process is obviously weakened after it is introduced into the controller. As can be seen from Figure 21c, although the output of the controller is saturated with the increase of the moment of inertia, since the reference trajectory already contains definite acceleration and deceleration stages, the controller saturation does not cause an obvious According to Figures 20a and 21a, similar to the rolling direction, the control response is significantly improved because the reference transition process includes high−frequency information. In comparison to Figure 19a, it can be seen that the reference transition trajectory information generated under the same circumstances is also the same due to the identical ATD parameters. Therefore, the ATD-falESO controller produces very similar transition processes in the pitch and roll channels. In addition, because falESO disturbance observer has a faster response speed, the influence of accumulated errors on the transition process is obviously weakened after it is introduced into the controller. As can be seen from Figure 21c, although the output of the controller is saturated with the increase of the moment of inertia, since the reference trajectory already contains definite acceleration and deceleration stages, the controller saturation does not cause an obvious increase in overshoot quantity. Lateral loop pitch has no obvious back and forth adjustment process, which means higher efficiency than PID controller.

Taillock Control Simulation
Taillock controller is an angular velocity tracking controller, which takes reference direction angular velocity . ψ r as input. Set K p = 1.0, K i = 3.0. The collective pitch limit of tail rotor output θ 0 TR by the controller is ±0.2 rad. Reference direction angular velocity signal is shown in Equation (70): Yaw rate . ψ corresponding to PID angular velocity controller and tail rotor collective pitch instruction θ 0 TR output by controller are shown in Figure 22. It can be seen that when there is no maneuver, the controller also has a collective pitch instruction of tail rotor of about 0.11 rad, which is used to balance the anti-torque generated by the main rotor rotation. Due to the anti-torque, the motion states along the two directions are obviously different. When maneuvering along the positive direction of the axis O b Z b the directional control moment is the sum of the anti-torque of the main rotor and the yaw moment generated by the tail rotor, and the consequent acceleration as well as the overshoot during the transition process is relatively large. When maneuvering in the opposite direction along the axis O b Z b , the directional control moment is the difference between the yaw moment generated by the tail rotor and the anti-torque. Under the same collective pitch limit, both the resulting acceleration and the overshoot in the transition process are small. rotor of about 0.11 rad, which is used to balance the anti-torque generated by the main rotor rotation. Due to the anti-torque, the motion states along the two directions are obviously different. When maneuvering along the positive direction of the axis the directional control moment is the sum of the anti-torque of the main rotor and the yaw moment generated by the tail rotor, and the consequent acceleration as well as the overshoot during the transition process is relatively large. When maneuvering in the opposite direction along the axis , the directional control moment is the difference between the yaw moment generated by the tail rotor and the anti-torque. Under the same collective pitch limit, both the resulting acceleration and the overshoot in the transition process are small. With the feedback gain fixed, introduce ATD reference trajectory generator and fa-lESO disturbance observer to constitute ATD-falESO angular speed controller. For the ATD reference trajectory generator, set = =̇, = 5.0 rad/s 2 , = −5.0 rad/s 2 . The curves of corresponding ̇a nd 0 are shown in Figure 23. Due to the acceleration limitation capability of ATD, the generated reference trajectory of yaw angular velocity has the same change rate in both positive and negative directions of . The yaw motion of the fuselage has the same performance in both directions even when subject to anti-torque of main rotor.

Conclusions
In this paper, according to the characteristics and application requirements of small unmanned helicopters, an asymmetrical tracking differentiator (ATD) is proposed. A speed and acceleration limiting mechanism whose threshold in two directions can be adjusted separately is added in the new ATD. This speed limiting mechanism makes the ATD suitable for applications including aircraft attitude control system. Based on the new ATD and fal−ESO disturbance observer, a programmable high−precision helicopter attitude control algorithm during transition is established. Aiming at the structure and task characteristics of the attitude control algorithm of an unmanned helicopter, a parameter estimation method for estimating the control gain coefficient and the main inertial characteristics of the fuselage is proposed. The simulation results show that, under the same conditions, the newly proposed ATD−based attitude control method has better control accuracy, faster transition process, smaller overshoot and higher control efficiency than the existing PID controllers. In the future, the application of ATD in trajectory control will be expanded. Due to the acceleration limitation capability of ATD, the generated reference trajectory of yaw angular velocity has the same change rate in both positive and negative directions of O b Z b . The yaw motion of the fuselage has the same performance in both directions even when subject to anti-torque of main rotor.

Conclusions
In this paper, according to the characteristics and application requirements of small unmanned helicopters, an asymmetrical tracking differentiator (ATD) is proposed. A speed and acceleration limiting mechanism whose threshold in two directions can be adjusted separately is added in the new ATD. This speed limiting mechanism makes the ATD suitable for applications including aircraft attitude control system. Based on the new ATD and fal-ESO disturbance observer, a programmable high-precision helicopter attitude control algorithm during transition is established. Aiming at the structure and task characteristics of the attitude control algorithm of an unmanned helicopter, a parameter estimation method for estimating the control gain coefficient and the main inertial characteristics of the fuselage is proposed. The simulation results show that, under the same conditions, the newly proposed ATD-based attitude control method has better control accuracy, faster transition process, smaller overshoot and higher control efficiency than the existing PID controllers. In the future, the application of ATD in trajectory control will be expanded.
Author Contributions: Data curation, C.Z. and Y.Y.; Funding acquisition, C.Z.; Methodology, C.C., Z.W. and C.Z.; Project administration, C.Z. All authors have read and agreed to the published version of the manuscript.