Active Disturbance Rejection Generalized Predictive Control of a Quadrotor UAV via Quantitative Feedback Theory

This paper investigates the design of a robust controller for a quadrotor unmanned aerial vehicle (UAV) system subject to parameter uncertainties and external disturbances. A new double closed-loop active disturbance rejection generalized predictive control (ADRC-GPC) is proposed for the trajectory tracking problem of the UAV. Considering the measurement noise, the ADRC-GPC control strategy can provide better dynamic performance and robustness against uncertainties and external disturbances than active disturbance rejection control (ADRC). Furthermore, a controller parameter tuning method for ADRC-GPC is proposed. In the framework of the classical two-degree-of-freedom equivalent model, the bandwidth of ESO is determined according to the parameter selection criterion. The criterion is a trade-off between the estimation accuracy and the immunity of measurement noise. In addition, the parameters of the ADRC-GPC controller are tuned by quantitative feedback theory to achieve the expected performance specifications. Finally, the proposed method is compared with several traditional control methods in the simulation experiments. The simulation results show that the proposed method has better set-point tracking performance, disturbance rejection performance, and robustness.


I. INTRODUCTION
A quadrotor UAV is a type of unmanned helicopter with four rotors. It is controlled by changing the lift forces of four rotors. Due to its simple structure, low inertia, and low cost, the quadrotor UAV is widely used in civil and military fields such as surveillance, search, surveying, and reconnaissance [1]. Trajectory tracking control is an important guarantee for the Quadrotor UAV to realize the above applications. The quadrotor UAV has six degrees of freedom, including three position variables and three attitude angle variables. However, there are only four control inputs. So it is a typical underactuated mechanical system. In addition, the quadrotor UAV has the characteristics of nonlinearity, multi-input The associate editor coordinating the review of this manuscript and approving it for publication was Ming Xu . multi-output (MIMO), strong coupling, and model uncertainty [2], traditional control strategies developed for fully actuated systems are not directly applicable to the quadrotor flight control systems. Research on the quadrotor flight control is one of the most interesting topics in the robotics and control community. In the past two decades, several control strategies for trajectory tracking were proposed [3]- [5], including PID control, sliding mode control, adaptive control, robust control, neural network control, etc.. In [6], a controller design method based on feedback linearization and backstepping was proposed for a quadrotor dynamic model. In [7], a further input saturation problem was considered. A backstepping controller based on the Nussbaum function was designed, and the latent singularities in the attitude extraction process caused by saturation nonlinearities were avoided. In [8]- [12], sliding mode control schemes were proposed to address the problems of underactuation, model uncertainty, actuator failure, and external disturbance. The stability of the closed-loop system was guaranteed with the Lyapunov theory. In [13], [14], model uncertainty and external disturbance were estimated by neural networks and the controller parameters were tuned online. Robust control is also an effective strategy to solve the control problem with parameter uncertainty. In [15]- [17], sliding mode observers and interval matrix approach were used to realize the trajectory tracking control of UAVs. It can be seen from the above references that underactuation, model uncertainty, and external disturbances are the main problems that need to be considered. The control strategy which does not rely on the exact model and has good disturbance rejection ability is a focused area in the trajectory tracking control of quadrotor UAVs.
Active disturbance rejection control (ADRC), which inherits many merits of PID controller, is emerging as a new robust control technology [18], [19]. In the framework of ADRC, the model uncertainties, nonlinearity, and external disturbance are regarded as the total disturbance, which is estimated and compensated by the extended state observer (ESO) in real-time. The disturbance rejection capability of ADRC has been verified in the field of UAV flight control. In [20], ADRC and wind feedforward compensation were used to realize the accurate flight path tracking control of a powered parafoil. The experiment showed that the ADRC method achieved good tracking performance and robustness against the variable wind disturbance. In [21], ADRC was applied to the attitude loop and the position loop at the same time. The simulation experiment proved that the proposed ADRC control scheme had good disturbance rejection performance and robustness. On this basis, ADRC was combined with sliding mode control to further improve the dynamic performance of the trajectory tracking control [22], [23]. In [24], an attitude controller was developed within the framework of ADRC and embedded model control, and a multi-step test strategy was proposed to assess the performance of ADRC controller. According to the above references, ADRC is good at dealing with the problem of parameter uncertainty and external disturbance. On the other hand, it should be noted that ADRC is an observer-based control method. When there is measurement noise or high-frequency disturbance in the output, the bandwidth of the ESO will be limited. At this time, the estimation accuracy of the total disturbance will be reduced, and the dynamic performance of the system will be difficult to guarantee.
To solve the problem of the ESO bandwidth limitation, a practicable solution is to replace the PD controller in ADRC with a more robust controller. In [25], [26], an active disturbance rejection generalized predictive controller (ADRC-GPC) was proposed. The system was compensated into an approximate integrator chain, and a generalized predictive controller was designed for the compensated system. The simulation result showed that the ADRC-GPC method had good dynamic performance. However, ADRC-GPC involves comparatively more parameters than ADRC. The parameter tuning is challenging in ADRC-GPC with model uncertainties. Quantitative feedback theory (QFT) [27], [28] is a robust control design methodology for the system with uncertainties. The model uncertainties and the performance specifications can be transformed into the boundaries in the Nichols chart, and loop shaping is realized by a graphic method. QFT has been successfully applied to a wide variety of industrial fields [29]- [31]. It is a feasible way to use QFT to guide the parameter tuning of the ADRC-GPC controller.
Motivated by afore-mentioned discussions, this paper proposes a trajectory tracking control structure based on ADRC-GPC and a parameter tuning method based on QFT. The main contributions of this paper are summarized as follows.
(1) According to the dynamic model of a quadrotor UAV, a closed-loop trajectory tracking control structure based on ADRC-GPC is proposed. By introducing the virtual control input, the original underactuated multivariable system is transformed into six single-input single-output (SISO) systems. Coupling effects and external disturbances between loops are estimated and compensated by ESO in real time.
(2) By the equivalent transformation, the two-degreeof-freedom equivalent model of ADRC-GPC is obtained. An ESO bandwidth selection criterion is proposed according to the equivalent model. The estimation accuracy and noise influence are taken as constraints in the criterion, and a reasonable bandwidth range can be obtained.
(3) A parameter tuning method of ADRC-GPC based on QFT is proposed. The controller parameters are tuned to meet the performance specifications in the Nichols chart. The system has good dynamic performance and disturbance rejection ability after parameter tuning.

II. DYNAMICS OF A QUADROTOR UAV
The structure of a quadrotor UAV used in this work is depicted in Figure 1. The dynamic model is introduced in detail in [32], and it contains six degrees of freedom: three position variables (x, y, z) and three attitude variables (θ, φ, ψ), θ represents the pitch angle, φ represents the roll angle, and ψ represents the yaw angle.
The simplified dynamic model of the quadrotor UAV is established according to Euler-Lagrange formulation as follow: where the parameter values of the UAV are shown in Table 1.
T is shown as follow: where k a is the coefficient between the rotor voltage and the lift generated by the rotor. The study objective in this work is to control three attitude angles and three position variables by the inputs

III. DESIGN OF TRAJECTORY TRACKING CONTROLLER FOR THE QUADROTOR UAV A. TRAJECTORY TRACKING CONTROL SCHEME OF THE QUADROTOR UAV
According to the dynamic model of the quadrotor UAV, the system includes six channels and two control loops. The longitudinal motion channel (x), the lateral motion channel (y), and the height motion channel (z) belong to the position control loop. The pitch angle channel (θ), the roll angle channel (φ), and the yaw angle channel (ψ) belong to the attitude control loop. For the attitude loop, each channel has a corresponding control input. However, for the position loop, there is only one control input U 1 . So three virtual control inputs are designed for the position control loop. The virtual  control inputs U x , U y , U z are expressed as According to equations (1) and (3), the real control input U 1 can be calculated as follows: where θ and φ in equations (5) and (6) are defined as the desired trajectories of the pitch angle and the roll angle, respectively. The designed trajectory tracking control scheme of the quadrotor UAV is shown in Figure 2.

B. DESIGN OF ADRC-GPC CONTROLLER
For the control structure depicted in Figure 2, the coupling effect between loops can be estimated and compensated by ESO. With the compensation, the system can be divided into six single input single output (SISO) systems. The estimation accuracy of the disturbance depends on the value size of ESO bandwidth. Larger bandwidth can obtain higher estimation accuracy, but it makes the controller sensitive to measurement noise. Thus, an ADRC-GPC is designed to enhance the robustness of the system when the bandwidth of ESO is limited. The specific control structure of ADRC-GPC is shown in Figure 3. The total disturbance of the system is estimated and compensated through ESO, and the generalized predictive controller is designed for the integrator chain after compensation. Coupling effects and external disturbances are defined as the total disturbance f , and the gravity compensation is added into the controller of the z-axis. The six channels of the quadrotor UAV can be rewritten as second-order systems as follow:ÿ where the total disturbance , a 1 and a 2 are the model parameters, w(t) is the sum of external disturbances and coupling effects, b and b 0 are the actual value and estimated value of the control input gains. Assume that the total disturbance f is differentiable and the derivative of f is h, equation (7) can be rewritten as the following state space form: where z 1 , z 2 and z 3 are the estimations of x 1 , x 2 and x 3 . β 1 , β 2 and β 3 are the error feedback gains of ESO. According to equation (9), the characteristic equation of ESO is Set the characteristic equation as λ(s) where ω o is the ESO bandwidth. With the pole configuration, the parameters of ESO are reduced to one. By tuning ω o , z 1 , z 2 and z 3 can estimate the values of y,ẏ and f accurately. The control law can be designed as Substituting equation (12) into equation (7), the relationship between u 0 and the output y is In ADRC-GPC, GPC controller is designed for the system (13). The following CARMA model is used as the prediction model: where u 0 (k) and y(k) are the input and output signals at time For simplicity, assume that C(z −1 ) = 1. With a zero-order hold, the discrete transfer function of equation (13) is where Z (·) is the z-transformation operator, T is the sampling time. According to equations (14) and (16), the CARMA model parameters of G 0 (z −1 ) can be expressed as In order to obtain the predictive model at k + j time, the following Diophantine equations are considered, where j = 1, 2, . . . , N , N is the prediction horizon. The forms of E j , F j , G j and H j are described as follows: Substituting equations (17) into Diophantine equations (18), the solution of the Diophantine equations can be obtained as According to equations (14)- (20), the prediction model of the system at k + j time is expressed as The vector form of equation (21) can be derived as , and N u is the control horizon. Consider the following performance index: where λ(λ > 0) is the control weighting factor. v(k + j) is the softening value of the reference signal, so that the output y(k) can reach the reference y r (k) smoothly, the softening reference trajectory is designed as where The performance index of equation (23) can be rewritten as According to equations (22) - (24), the solution of (25) can be calculated explicitly, Take the first element of the control sequence U as the current control input u 0 (k): Let h T = 1 0 · · · 0 (G T G + λI) −1 G T , and equation (27) is expressed as

IV. PARAMETER TUNING BASED ON QFT
In the ADRC-GPC algorithm, the parameters that need to be tuned include: ESO bandwidth ω o , sampling time T , predictive horizon N , control horizon N u , softening factor α, and control weighting factor λ. Due to the large number of parameters, it is difficult to achieve the desired performance by manual tuning. The two-degree-of-freedom unit feedback model is taken as a standard control structure in QFT, and the loop shaping of the controller can be realized by a graphical way in the Nichols chart. In this work, QFT is applied to the parameter tuning of ADRC-GPC, the controller parameters are tuned to meet the required performance specifications. According to equations (9) and (11), the outputs of ESO can be expressed as By equations (29), the internal model structure of the ESO is shown in Figure 4, where n(t) represents the effect of noise on the ESO, G(s) represents the transfer function of one channel in the UAV model. Take the position loop as an example. After a gravity compensation is introduced to the z-axis loop, and the effects of interactions are regarded as disturbances, the plants of the x, y, and z channels can be expressed as where a =  (12), (29), and (30), the equivalent plant between the system output y and the virtual control input u 0 can be readily derived as According to equation (28), the control law of GPC can be rewritten as where Assume that the discrete form of the equivalent plant G eq with a zero-order hold is expressed as G d (z −1 ), and there is y(k) = G d (z −1 )u 0 (k). By equation (32), the closed loop system is described as follow: where T (z −1 ) As shown in Figure 5, the closed-loop control structure of GPC is obtained according to equation (33). Further combining with Figure 4, the closed-loop control structure of ADRC-GPC can be obtained as shown in Figure 6. By proper transfer function transformation, the closed-loop structure shown in Figure 6 can be transformed into the two-degreeof-freedom equivalent model shown in Figure 7, where G eq is the equivalent controlled plant G eq = G p (s)G 0 (s).
G p (s) is the transfer function of ω o . F p (z −1 ) and C(z −1 ) are the prefilter and the controller respectively. The following parameter tuning will be carried out based on the above equivalent transformations. The parameter tuning rule will be divided into two steps: (1) the appropriate ESO bandwidth parameter will be selected with the parameter perturbation of the UAV plant; (2) QFT will be used for the parameter tuning of the ADRC-GPC controller.

B. SELECTION CRITERIA OF THE ESO BANDWIDTH
The total disturbance is observed and compensated through ESO, and the estimation accuracy of ESO depends on the value of bandwidth ω o . Generally, the estimation accuracy is higher with a large bandwidth. For the control structure shown in Figure 7, the purpose of ESO bandwidth selection is to ensure that G p (jω) is as close to 0dB as possible. Thus, the dynamic characteristics of the equivalent model G eq are close to the nominal model G 0 (s). However, the rise of ESO bandwidth will increase the influence of measurement noise on the system. Therefore, it is necessary to consider the effects of estimation accuracy and measurement noise at the same time. According to Figure 4, the influence of measurement noise on the observer can be expressed as     Figure 8, it can be seen that G p (jω) is close to 0dB with the increase of ω o . When ω o is fixed, the log magnitude-frequency characteristic plot of the system presents a monotonous increasing characteristic in the low frequency band. It tends to be horizontal in the high frequency band. According to the above characteristics, the solution of constraint (36) can be obtained by only selecting the minimum frequency point ω = 0.01 to verify the constraint. For any G j , the value of ω o should satisfy the constraint (36). The frequency characteristic of G n (s) is shown in Figure 9. The log magnitude-frequency characteristic plot of G n (s) presents a monotonous decreasing characteristic in the high frequency band. Therefore, the solution of constraint (37) can be verified VOLUME 10, 2022 by selecting the minimum frequency point ω = 200 to verify the constraint. Finally, the interval of ω o in the position loop that satisfies the constraints (36) and (37) is [23.4, 66.5]. Considering the coupling effects as the input disturbances, the interval of ω o in the attitude loop can also be obtained according to the above method in the range of [37.5, 66.7].

C. PARAMETER TUNING BASED ON QFT
Based on the determined bandwidth, QFT is adopted to design C(z −1 ) and F p (z −1 ) in Figure 7. The design procedure is mainly divided into four steps: QFT templates definition, performance specifications design, loop shaping and prefilter design.

1) QFT TEMPLATES DEFINITION
The uncertainties of the generalized plant G eq can be expressed as A discrete number of plants are defined by gridding the above uncertain parameters between the minimum and the maximum values. The distribution between the minimum and the maximum is linear. The frequency points of interest are selected according to the reference [28], and the selected set is ω = 0.01 0.05 0.1 0.5 1 5 10 20 rad/s. The QFT templates of the frequency points are shown in Figure 10. Each circle in Figure 10 is an uncertain template. The distribution of the circles represents the uncertainty of the plant at a specific frequency point. The sign ''+'' in Figure 10 indicates the nominal plant at the selected frequency. The nominal plant is used to design the controller and prefilter. It can be seen from Figure 10 that the open-loop phases of all templates are distributed around −180 degrees. This shows that the frequency characteristics of the generalized plant G eq (s) with model uncertainties is still very close to the series integral structure 1 s 2 . The distribution of the templates is concentrated, and this is beneficial to the following controller design.

2) PERFORMANCE SPECIFICATIONS DESIGN
In QFT, the desired performance specifications of stability and disturbance rejection need to be converted into a series of constraint boundaries for the frequency response L(jω) = C(jω)G eq (jω) of the nominal model in the Nichols chart. By designing the controller C(z −1 ) and the prefilter F p (z −1 ), the open-loop frequency response curve of the nominal plant will be adjusted to meet the boundary conditions. For the stability and disturbance rejection problem of the trajectory tracking control, the performance specifications are considered as follows: (1) Robust stability specification The magnitude constraint of the closed-loop system is defined as  where W s is the closed-loop resonance peak, the relationship between W s , magnitude margin (GM), and phase margin (PM) is shown as follows The robust stability boundary of the position loop is selected as W s = 1.17, the magnitude margin of the corresponding system is 5.4dB, and the phase margin is 50 • .
(2) Output disturbance rejection specification The output disturbance rejection specification of the system is designed as where the output disturbance rejection specification of the position loop is selected as δ 2 (jω i )= (jω i ) (jω i )+2 . (3) Input disturbance rejection specification The generalized input disturbances of the system include: the estimation error of ESO, the interaction between channels, and the system parameter perturbations. The specification of the input disturbance rejection is designed as where the input disturbance rejection specification of the position loop is selected as δ 3 (jω i )= (jω i ) (jω i )+0.1 . In the constraints (38)-(41), ω i ∈ ω, G eq ∈ G eq . G eq can be expressed in its polar form as G eq (jω i ) = pe jβ , the controller is expressed as C(jω i ) = qe jγ . When the selected frequency point ω i is fixed, δ 2 (jω i ) = δ 2 , δ 3 (jω i ) = δ 3 . p, β, δ 2 , and δ 3 are the constants. Substituting the polar forms into equations (38)-(41), the above performance specifications can be transformed into the following quadratic inequality equations: The phase γ of the controller C(jω i ) is discretized in the bounded interval −360 • , 0 • , like γ ∈ −360 • :5 • :0 • . When γ is fixed, the unknown parameter of the inequalities (42) is the controller magnitude q. According to the reference [33], the quadratic inequalities of equations (38), (40), and (41) are solved and translated into a set of curves on the Nichols chart for each frequency of interest and type of specification.

3) LOOP SHAPING
By constructing the gain, zero, and pole of the controller, the open-loop frequency response L(jω) can meet the boundaries of the performance specifications. L(jω) needs to be located above the boundary at each corresponding frequency, and does not intersect the robust stability boundary at the high frequency. The controller is designed by loop shaping graphically. The order of the controller C(z −1 ) is determined by the ADRC-GPC equivalent structural (34). According to equation (34), the zero and the pole of C(z −1 ) are mainly determined by the parameters N and N u in ADRC-GPC. The gain of C(z −1 ) is mainly determined by the parameters λ and α. According to the above rules, the controller parameters are tuned to meet the specification boundary. After parameter tuning, the controller of the position loop is designed as The  Figure 11. The black dashed and solid lines in the figure are the frequency response curves before and after loop shaping. The color curves are the boundaries calculated according to the performance specifications at each frequency point. The marking points of each color on the frequency response curves should be distributed above the corresponding color solid line or below the color dotted line. It can be seen from Figure 11 that the curve of L(jω) after correction is above the colored solid line at the corresponding frequency point in the Nichols chat. Therefore, the designed controller can meet the above performance specifications when there are uncertainties in the plant.

4) PREFILTER DESIGN
The tracking performance of the system can be achieved by designing the prefilter in Figure 7. The tracking boundaries of the position loop are constrained by the upper and lower bounds as follows: where (jω) = G eq (jω)C(jω) 1 + G eq (jω)C(jω) According to the ADRC-GPC parameters determined in the loop shaping, the prefilter is designed as The tracking performance of the system is shown in Figure 12. The output frequency response of the corrected system is within the tracking boundaries. If the parameters do not meet the performance specifications, it is necessary to return to tune the controller parameters λ and α again. The controller parameters of the attitude loop can also be tuned VOLUME 10, 2022

V. SIMULATION
In this section, comparative experiments are designed to highlight the performance and robustness of the ADRC-GPC method based on QFT. In order to demonstrate its numerical performance, ADRC and PID methods [21] are compared and analyzed. In the trajectory tracking experiment, the initial position of the UAV is  Table 2. It can be seen that ADRC-GPC  scheme achieves the best tracking performance in the position and attitude loops. The IAE index of the ADRC-GPC scheme is better than the ADRC and the PID schemes in the x, y, and ψ loops. The IAE index of the PID scheme is small in the z loop, but the PID scheme has a large overshoot. The ADRC-GPC scheme has a short settling time and no overshoot in the z loop. So the dynamic performance of the ADRC-GPC scheme is better than the PID scheme in the z loop.

B. DISTURBANCE REJECTION
In this section, the model uncertainties, measurement noise, and external disturbances are considered in the simulation experiment. The model parameters in Table 1 The reference trajectory of the position is from (1,1,1)m to (0,0,0)m, and the set-point of the yaw angle is ψ d = 1. The simulation results are shown in Figures 15-16. The IAE index is presented in Table 3. The IAE index of the ADRC-GPC scheme is better than the ADRC and the PID schemes in    the x, y, z and ψ loops. It can be seen that the PID control scheme is greatly affected by the disturbances and noise. The tracking trajectory of the PID control scheme has oscillation and the convergence speed is slow. By contrast, the dynamic performance of the ADRC and ADRC-GPC control schemes are much better than PID. The tracking error of ADRC-GPC is significantly smaller than the other two schemes. It is shown that the ADRC-GPC scheme has the best stability and disturbance rejection ability.

C. CIRCULAR TRAJECTORY TRACKING
In order to further verify the performance of the proposed method, a circular trajectory tracking experiment is carried   out with the disturbances. The model uncertainties, measurement noise, and external disturbances are consistent with   the previous section. The starting point is (1,1,1)m, and the desired tracking trajectory is a unit circle centered on the origin. The set-point of the yaw angle is ψ d = 1. The simulation results are shown in Figures 17-19. The IAE index is presented in Table 4. The IAE index of the ADRC-GPC scheme is better than the ADRC and the PID schemes in the z and ψ loops. The IAE index of the PID scheme is small in the x and y loop. But in Figures 17 and 18, it can be seen that the PID and the ADRC control schemes have large position tracking error in the initial time. The overshoot of the yaw angle in ADRC-GPC scheme is smaller than the PID and ADRC. According to Figures 19, the tracking trajectory of ADRC-GPC can track the reference signal more accurately and quickly. It is worth noting that the amplitudes of the x and y loops in the ADRC-GPC scheme are larger than the PID and ADRC schemes. This is caused by the measurement noise in the attitude and position loop. A filter design should be considered when the noise power is large.

VI. CONCLUSION
This paper proposed a new robust control scheme and a corresponding tuning rule for the trajectory tracking problem of a quadrotor UAV. The ADRC-GPC control scheme was designed based on active disturbance rejection control and generalized predictive control. The proposed control scheme improved the robustness and dynamic performance of the system when the ESO bandwidth was limited. The tuning rule of ADRC-GPC was proposed based on a two-degree-offreedom equivalent model and quantitative feedback theory. The tuned parameters of ADRC-GPC achieved the expected stability and performance specifications with model uncertainties and external disturbances. Moreover, the discrete form of ADRC-GPC was derived, and it was convenient for practical applications. Finally, the simulation experiments demonstrated the efficacy of the proposed control scheme and the tuning rule. Some practical applications and an ESO based on the Kalman filter will be carried out in future research.