Comprehensive Engineering Frequency Domain Analysis and Vibration Suppression of Flexible Aircraft Based on Active Disturbance Rejection Controller

The crash of an aircraft with an almost vertical attitude in Wuzhou, Guangxi, China, on 21 March 2022, has caused a robust discussion in the civil aviation community. We propose an active disturbance rejection controller (ADRC) for suppressing aeroelastic vibrations of a flexible aircraft at the simulation level. The ADRC has a relatively simple structure and it has been proved in several fields to provide better control than the classical proportional-integral-derivative (PID) control theory and is easier to translate from theory to practice compared with other modern control theories. In this paper, the vibration model of the flexible aircraft was built, based on the first elastic vibration mode of the aircraft. In addition, the principle of ADRC is explained in detail, a second-order ADRC was designed to control the vibration model, and the system’s closed-loop frequency domain characteristics, tracking effect and sensitivity were comprehensively analyzed. The estimation error of the extended state observer (ESO) and the anti-disturbance effect were analyzed, while the robustness of the closed-loop system was verified using the Monte Carlo method, which was used for the first time in this field. Simulation results showed that the ADRC suppressed aircraft elastic vibration better than PID controllers and that the closed-loop system was robust in the face of dynamic parameters.


Introduction
In recent years, civil aviation companies have tended to adopt aircraft designs with larger aspect ratios to reduce their induced drag, achieve speed breakthroughs and reduce fuel consumption. However, this design sacrifices the aircraft's stability and produces flight vibrations [1]. Compared to conventional designs, slender aircraft have larger wave resistance and are subject to increased reaction forces. Therefore, as civil aircraft are improved, more pronounced and sustained vibrations occur [2,3]. Corresponding to the fatigue theory of engineering mechanics, elastic vibrations in flexible aircraft can reduce their service life and may even damage the mechanical structure at certain specific frequencies, leading to more severe effects [3].
In recent years, scholars have been focusing on effectively suppressing the elastic vibration of aircraft to improve its performance. Chatter suppression techniques can be divided into two main categories: passive suppression and active suppression. Passive chatter suppression mainly seeks more stable operation at the expense of vehicle performance, such as increasing airframe weight and reducing flight distance [4]. Active suppression is mainly to effectively suppress the flutter phenomenon by optimizing the aircraft's control system. Active flutter suppression (AFS) technology can be traced back to the research content of the famous project, active flexible wing (AFW), in the late 20th century [5]. Currently, the more mainstream AFS techniques include classical PID control, linear quadratic Gaussian (LQG) control, H ∞ Control, sliding mode (SM) control, and fuzzy logic controller (FLC) , all of which play essential roles in modern control theory and The invention of the PID control algorithm is primarily credited to Ziegler and Nichols, who introduced the concept of PID to the world in 1942 [6]. PID is the most widely used control algorithm and dominates the controller design subject. About 70 years later, Wang et al. developed a vibration model of the aircraft takeoff and landing process to facilitate the performance of the landing gear system through PID control [7]. After nearly a century of development, the functionality of PID has been over-exploited to the point that breakthroughs have been challenging. Many scholars have tried combining PID with different reference-seeking models in recent years to improve the control effect [8][9][10].
Gupta first proposed LQG control in 1980 and compared it with conventional controllers [11]. Ahmad, S., A. and Na, S., et al. have tried combining LQG with different feedback control strategies [12,13]. After 10 years of development, Franciszek combined the H ∞ method with LQG to give a theoretical solution for nonlinear aeroelastic vibrations [14], the H ∞ controller. Since its introduction in 1993, the control theory has been continuously discussed by scholars [15]. The works of Navya and Alexandr V. conducted a simulation analysis of the H ∞ controller in the aerospace field [16,17]. The sliding model was first proposed in 1990 and applied to aerospace control at the beginning of the 21st century [18]. In the following 20 years, scholars have improved SMC strategies with different compensation methods [19][20][21]. In recent years, deep learning has brought epic changes in various fields. Han et al. proposed a controller combining SMC and reinforcement learning to achieve greater robustness and control accuracy [22]. FLC is also an essential player in modern control theory, and its performance has been widely discussed in the aerospace control field [23][24][25][26]. However, the researchers realized that there are several problems with the above algorithms, such as the PID algorithm seems too simple to adapt to industry application. The assumptions of modern control algorithms, such as LQG and SMC, are valid. However, the mathematical form is too complex and requires a high level of controller design, making it difficult to achieve a breakthrough from simulation to application. The invention of the PID control algorithm is primarily credited to Ziegler and Nichols, who introduced the concept of PID to the world in 1942 [6]. PID is the most widely used control algorithm and dominates the controller design subject. About 70 years later, Wang et al. developed a vibration model of the aircraft takeoff and landing process to facilitate the performance of the landing gear system through PID control [7]. After nearly a century of development, the functionality of PID has been over-exploited to the point that breakthroughs have been challenging. Many scholars have tried combining PID with different reference-seeking models in recent years to improve the control effect [8][9][10].
Gupta first proposed LQG control in 1980 and compared it with conventional controllers [11]. Ahmad, S., A. and Na, S., et al. have tried combining LQG with different feedback control strategies [12,13]. After 10 years of development, Franciszek combined the H ∞ method with LQG to give a theoretical solution for nonlinear aeroelastic vibrations [14], the H ∞ controller. Since its introduction in 1993, the control theory has been continuously discussed by scholars [15]. The works of Navya and Alexandr V. conducted a simulation analysis of the H ∞ controller in the aerospace field [16,17]. The sliding model was first proposed in 1990 and applied to aerospace control at the beginning of the 21st century [18]. In the following 20 years, scholars have improved SMC strategies with different compensation methods [19][20][21]. In recent years, deep learning has brought epic changes in various fields. Han et al. proposed a controller combining SMC and reinforcement learning to achieve greater robustness and control accuracy [22]. FLC is also an essential player in modern control theory, and its performance has been widely discussed in the aerospace control field [23][24][25][26]. However, the researchers realized that there are several problems with the above algorithms, such as the PID algorithm seems too simple to adapt to industry application. The assumptions of modern control algorithms, such as LQG and SMC, are valid. However, the mathematical form is too complex and requires a high level of controller design, making it difficult to achieve a breakthrough from simulation to application. The active disturbance rejection control (ADRC) algorithm, as an innovative product of PID theory, is expected by developers to achieve a breakthrough from theory to practical application of modern control algorithms. Jingqing Han formally proposed ADRC in 1998. The core of ADRC is the extended state observer (ESO), which estimates the system's state and compensates for perturbations using the control method, turning the controlled object into a purely integral object [27]. Han pointed out the problems of PID algorithms and tried to use ADRC to solve the dilemma that PID has not made breakthroughs in development for decades [27][28][29]. Gao simplified Han's work to obtain the linear active disturbance rejection control (LADRC) technique, which makes parameter tuning easier [30].
Regarding vibration suppression in aerospace, the control effectiveness of ADRC has been affirmed by researchers at the simulation level. Liu et al. combined an electromechanical actuator (EMA) with LADRC to propose a dynamic aircraft servo system [31]. Yang et al. used ADRC to suppress the elastic vibration of a high aspect ratio UAV fuselage [32]. Chen and Zhao proposed a Support Vector Machine (SVM) ADRC control strategy for UAVs to suppress the system chattering [33]. Wang et al. proposed an improved ADRC strategy to simulate and analyze hypersonic aerospace models [34]. Duan et al. modeled the helicopter load system to study the effect of load sway on the helicopter and confirmed the effectiveness of ADRC [35]. Wang et al. analyzed the limit cycle vibration suppression in hovering conditions. The results showed that ADRC has outstanding performance in suppressing limit cycle vibration and eliminating attitude control errors [36]. The work of Qiao and Zhong et al. brought innovation to ADRC control theory and applied it to the simulation analysis of aerospace control [37,38]. ADRC has been applied in some fields, such as financial [39] and industrial [40] fields, and some degree of progress has been made due to the application of ADRC. For example, in 2010, ADRC was first applied to a factory, reducing its energy consumption by 41% [40]. Therefore, ADRC algorithm's role in aircraft controller design is also highly anticipated and has a high research value.
On 21 March 2022, the near-vertical crash of China Eastern Airlines plane MU5735 became a worldwide sensation. It was the deadliest air crash in China in 28 years. Across the globe, air crashes occur every year. Although air crashes do not occur as frequently as traffic accidents, the damage caused is enormous. After every air crash, there seems to be a warning that controllers need to be reformed to suit the needs of the airlines. In this study, we take a comprehensive engineering frequency domain analysis of the aeroelastic vibration model, mainly for flexible aircraft. In addition, we combined the developed aircraft elastic vibration model with the designed second-order active disturbance suppression controller to analyze the closed-loop control effect and the tracking performance and sensitivity of the closed-loop system. Finally, the Monte Carlo shooting method verified the system's robustness.  [38]. The general mathematical model of flexible aircraft pitch angle difference and the first elastic mode is given by [38,42]:

Flexible Aircraft Aerodynamic Vibration Model
..
The two equations in Equation (1) represent the rigid body aerodynamic model channel of the aircraft pitch angle and the elastic vibration channel of the flexible aircraft, respectively. Figure 2 is a schematic diagram of the mathematical model in a flexible aircraft. Where θ(t) is the pitch angle difference, that is, the difference of the pitch angle at immediate posture (in time) and the initial posture, u(t) is the control signal, q 1 (t) is the first generalized elastic mode of vibration, and it simulates the elastic vibration of the aircraft during flight. ξ and ω 1 are the damping ratio and the natural frequency of the elastic mode, respectively, c 3 represents the control gain of the elastic mode q 1 (t), c 1 , c 2 , b 11 and b 21 are the coupling coefficients between the pitch angle and the first mode, α is the attack angle, d 1 (t) and d 2 (t) simulates the external disturbance of the system.
The two equations in Equation (1) represent the rigid body aerodynamic model channel of the aircraft pitch angle and the elastic vibration channel of the flexible aircraft, respectively. Figure 2 is a schematic diagram of the mathematical model in a flexible aircraft. Where ( ) is the pitch angle difference, that is, the difference of the pitch angle at immediate posture (in time) and the initial posture, ( ) is the control signal, 1 ( ) is the first generalized elastic mode of vibration, and it simulates the elastic vibration of the aircraft during flight. and 1 are the damping ratio and the natural frequency of the elastic mode, respectively, 3 represents the control gain of the elastic mode 1 ( ), 1 , 2 , 11 and 21 are the coupling coefficients between the pitch angle and the first mode, is the attack angle, 1 ( ) and 2 ( ) simulates the external disturbance of the system. The dynamic system's Laplace-transfer matrix of zero initial condition is obtained. It is worth noting that the study in this section focuses only on the control object itself, so the external disturbance 1 ( ) and 2 ( ) are temporarily ignored in this section, and the angle of attack = 0 is assumed.  (1)). Among them, pitch angled difference represents physical quantity θ(t), and pitch angle is equal to climb angle plus attack angle (α).
The dynamic system's Laplace-transfer matrix of zero initial condition is obtained. It is worth noting that the study in this section focuses only on the control object itself, so the external disturbance d 1 (t) and d 2 (t) are temporarily ignored in this section, and the angle of attack α = 0 is assumed. According to Gramer's Law, the transfer function matrix is calculated: Among them, [38]. Bode plots of the transfer function G 1 and G 2 are shown in Figure 3.
According to Gramer's Law, the transfer function matrix is calculated: Among them, . For the transfer function 1 , the bandwidth ( ) is 1 = 1.28 rad/s, the resonant peak ( ) and resonant frequency ( ) are 1 = 3.26 dB, 1 = 0 rad/s respectively, while for the transfer function 2 , 2 is equal to 43.38 rad/s, 2 is 30.60 dB, while 2 is 30.00 rad/s (equals to 1 ). In general, the system bandwidth is relatively large, so the system has a strong tracking capability for step signals. The resonant peaks and frequencies are relatively large, so the system is less stable in response to step signals. For the transfer function G 1 , the bandwidth (w) is w 1 = 1.28 rad/s, the resonant peak (Mr) and resonant frequency (R f ) are Mr 1 = 3.26 dB, R f 1 = 0 rad/s respectively, while for the transfer function G 2 , w 2 is equal to 43.38 rad/s, Mr 2 is 30.60 dB, while R f 2 is 30.00 rad/s (equals to ω 1 ). In general, the system bandwidth is relatively large, so the system has a strong tracking capability for step signals. The resonant peaks and frequencies are relatively large, so the system is less stable in response to step signals.

Active Disturbance Rejection Controller Design
The ADRC control system consists of two main components, namelyu the extended state observer (ESO) and the proportional difference (PD) control method. In this section, we designed a second-order ADRC from these two parts, summarized how to simplify the ADRC parameter tuning, and analyzed the closed-loop system and its sensitivity.

Extended State Observer
Any dynamic system can be expressed by y (n) + a 1 y (n−1) + · · · + a n−1 y (1) + a n y = k(u (m) Among them, b 0 = b m /a n−2 , f contains all the terms in Equation (5) except y (2) and u. Therefore, f is expressed by a linear combination of other terms, that is Term f is defined as the generalized disturbance, including all perturbations inside and outside the system, which can be estimated by ESO to obtainf . The equation below is a control law involved in ADRC: which u 0 will be discussed in PD controller design section. Rewrite Equation (6): . .
where f −f is the estimation error, and Equation (6) can be rewritten as the following state-space equations: .
Among them, ESO is a stuff used to estimate x 1 , x 2 and f . The estimated state-space equations are shown in the equation is the gain vector of the ESO, and directly affects the performance of the observer. Rewriting Equation (12) yields: When the determinant of the matrix (A − LC) is equal to zero, the system can converge.
The characteristic polynomial p (n) can be expressed as: where α 1 , α 2 , α 3 are the adjustable parameters. To simplify the work of parameter adjustment, the equation below has been proposed [30].

PD Control Law
Assuming that the estimated valuef has accurate estimation, then we obtain . .
Then, u 0 should be concerned, and the problem is solved by classical PD control law, By rewriting Equations (17) and (18), the transfer function of the PD control law can be obtained To simplify the PD law's parameters, the equation below has been proposed [30], where k d = 2ω c and k p = ω 2 c , the final value theorem tells us that the pole occurs at −ω c , and if ω c is greater than zero, the system will eventually stabilize [3].

Closed-loop Analysis
For the flexible aircraft vibration model in the previous section, it is suggested to use the second-order ADRC controller described above for its control. In order to simulate the natural data acquisition pattern during flight, it is assumed that there is a specific relationship equation between the first mode q 1 (t) and the pitch angle error θ(t) [38], where θ m (t) is the state value obtained directly by the aircraft observer and 0.2 is the measurement factor. Combining ESO with the vibration model and rewriting the state-space Equation (10) obtains Among them, In order to obtain the control signal u(t), it can be obtained by solving Equation (18), u 0 (t), which is then converted to u(t) from Equation (8). Simplifying the progress obtains: The aircraft model and the ADRC controller were successfully combined to form a feedback closed-loop (see Figure 4), where G p represents the controlled object and G c and F represents the active disturbance rejection controller (ADRC) [43].
ESO estimates the state vector and obtains the estimating equation: In order to obtain the control signal ( ), it can be obtained by solving Equation (18), 0 ( ), which is then converted to ( ) from Equation (8). Simplifying the progress obtains: The aircraft model and the ADRC controller were successfully combined to form a feedback closed-loop (see Figure 4), where represents the controlled object and and represents the active disturbance rejection controller (ADRC) [43].  The transfer function of the closed-loop system needed to be analyzed to predict the characteristics of the system. The transfer function of ADRC we designed can be given by [43]: The transfer function of the controlled object is given by Equation (3), and combined with Equation (27), the closed-loop transfer function G cl can be found:

Parameter Tuning
By simplifying the parameters, there were finally three parameters left, namely, b 0 , ω 0 and ω c . The designed ADRC was simulated with the model for closed-loop feedback control, and the parameter values of the ADRC were trained using a unit-step setup signal. The parameter tuning of the model was limited by the range of values of some parameters in Table 1. Please refer to the following steps for the parameter tuning method.

•
As shown in Table 1, the six parameters of the elastic vibration model take random values within the intervals (simulation part), and each parameter has two boundary values. We know that if the system converges under a combination of 64 boundary values, any random set of parameters in the interval converge under the closed-loop system [38].  Through continuous training of the model, the controller parameters were finally determined, ω 0 = 7.8 rad/s, ω c = 1.68 rad/s and b 0 = −12.6, respectively. In real life, almost every field involving control can find the shadow of the PID controller. It is no exaggeration to say that among all controllers, PID is the most widely used controller at present. In order to show more clearly the advantages of ADRC compared to PID, a closedloop containing a PID controller was considered. Considering the parameter tuning method mentioned above and the optimization parameters of the genetic algorithm, the robust characteristics of the PID controller were ensured [44], and the parameters of the designed controller were the proportional coefficient (K P = −0.7), the integral coefficient (K I = −0.2), the derivative coefficient (K D = −0.7) and the filter coefficient (N = 1), respectively. Figure 5 illustrates the tracking effect of ADRC and PID closed-loop results. The ADRC closed-loop control was far superior to the PID loop for pitch angle difference and control signal input. At the same time, the PID controller lost credibility in the face of the fourth-order instability object. Focusing on Figure 5a, for ADRC loop, the pitch angle error peaked at 1.2942 rad at 3.58 s with a maximum overshoot of 29.42%. The 1% setting-point was used as the upper and lower limits of the acceptable state. The setting time of the closed-loop was considered as the point where it first maintained the acceptable range, that is, θ m (t) is in [0.99,1.01]. Therefore, the setting time of the ADRC loop was 8.65 s. Analyzing the control signal (Figure 5b), the minimum signal was minus 0.22 rad at 1.00 s, and the maximum signal was 0.43 rad at 3.34 s. Such a smooth control signal was appreciated. Overall, the designed ADRC controller achieved a good control effect by parameter adjustment.
Moreover, the frequency domain of the closed-loop transfer function needed to be analyzed. Analyzing the rad solid-lines of Bode diagrams (Figure 6) of the closed-loop TFs Equation (28), for the transfer function G cl1 , the amplitude margin (Gm) was equal to 49.43 dB, and the phase margin (Pm) accounted for 60.34 • , the cutoff frequency (C f ) was 30.00 Hz, while for G cl2 , the Gm = 28.08 dB, the Pm = −89.39 • , and the C f = 19.04 Hz. The relatively high cutoff frequency of the closed-loop system reflected the more extraordinary ability of the system to cope with the dynamic response. The high amplitude margin indicated the higher inclusivity and stability of the closed-loop system. Equation (28), for the transfer function 1 , the amplitude margin ( ) was equal to 49.43 dB, and the phase margin ( ) accounted for 60.34 ∘ , the cutoff frequency ( ) was 30.00 Hz, while for 2 , the = 28.08 dB, the = −89.39 ∘ , and the = 19.04 Hz. The relatively high cutoff frequency of the closed-loop system reflected the more extraordinary ability of the system to cope with the dynamic response. The high amplitude margin indicated the higher inclusivity and stability of the closed-loop system.

Closed-loop Sensitivity
Sensitivity is an essential indicator of the robustness of a closed-loop control system. For the closed-loop TFs (Equation (28)), the sensitivity function is given by [45]: Given that the parameters of the controller transfer function are fixed, while the parameter set of the control object is variable, bringing Equations (3) and (28) into Equation (29), we obtain, by rewriting: The Bode diagrams of 1 and 2 have been shown by the blue dotted-lines (see Figure

Closed-Loop Sensitivity
Sensitivity is an essential indicator of the robustness of a closed-loop control system. For the closed-loop TFs (Equation (28)), the sensitivity function is given by [45]: Given that the parameters of the controller transfer function G c are fixed, while the parameter set of the control object G p is variable, bringing Equations (3) and (28) into Equation (29), we obtain, by rewriting: The Bode diagrams of S 1 and S 2 have been shown by the blue dotted-lines (see Figure 6). The Nyquist lines of F(jω)G c (jω)G 1 (jω) and 0.2F(jω)G c (jω)G 2 (jω) are planned to be plotted in an imaginary coordinate system, and the point with the smallest distance from the plotted line to the coordinate (−1,0) is found, while the maximum value of sensitivity can be obtained. Its significance is elaborated in Wang's book [46], where the maximum amplitude of the sensitivity M s is given by the following Equation: where r refers to the radius of inscribed circle with the center point (−1,0), which is tangent to the Nyquist graph line in the imaginary coordinate system (see Figure 7).

Simulation
In this section, the estimation error of ESO is first discussed, and then the parameters of the aircraft model are randomly selected according to the Gaussian distribution. Simulation results for each set of parameters are analyzed, and the system's robustness is discussed. Figure 8 shows the mathematical model used for the simulation, implemented by the Simulink of Matlab.  Figure 7 shows the results of the sensitivity analysis of the closed-loop system with two inscribed circles of radius r 1 and r 2 , which were 0.62 and 0.57, respectively. Therefore, the maximum amplitude of the sensitivity was Ms 1 and Ms 2 , which were 1.61 and 1.74, respectively. They were numerically equal to the resonance peaks analyzed from the Bode plot. The maximum amplitudes were less than 2, indicating that the system was robust in the face of parameter variations [46].

Simulation
In this section, the estimation error of ESO is first discussed, and then the parameters of the aircraft model are randomly selected according to the Gaussian distribution. Simulation results for each set of parameters are analyzed, and the system's robustness is discussed. Figure 8 shows the mathematical model used for the simulation, implemented by the Simulink of Matlab.

Simulation
In this section, the estimation error of ESO is first discussed, and then the parameters of the aircraft model are randomly selected according to the Gaussian distribution. Simulation results for each set of parameters are analyzed, and the system's robustness is discussed. Figure 8 shows the mathematical model used for the simulation, implemented by the Simulink of Matlab.

ESO Estimation Error Analysis
A superiority of ADRC over the classic PID control method is its ability to provide an observational estimate of the state. The control effectiveness of ADRC for pitch angle vibration suppression is closely related to the accuracy of the ESO estimate.
As shown in Figure 9a-c, the maximum estimation error of θ m (t) occurred at 1.19 s with an error of 0.02 rad, and the estimation of θ m (t) by ESO was very accurate because it largely coincided with the true value. The maximum estimation error of . θ m (t) occurred at 1.13 s with an estimation error of 0.64 rad/s. In general, the ESO's estimation to . θ m (t) was relatively accurate. For the total perturbation f (t) estimation, the maximum estimation error of 14.78 rad/s 2 occurred at 1.00 s. The effectiveness of ESO in estimating the total perturbation seemed less satisfactory, but such a result was acceptable because it included the total disturbance of the system, which had larger uncertainty. The estimation error of f (t) stabilized in the range of [-1,1] after 2.77 s and, eventually, became stable as the system converged. Overall, the ESO estimates of the state variables largely overlapped with the system's output response, demonstrating the superiority of ADRC.

ESO Estimation Error Analysis
A superiority of ADRC over the classic PID control method is its ability to provide an observational estimate of the state. The control effectiveness of ADRC for pitch angle vibration suppression is closely related to the accuracy of the ESO estimate.
As shown in Figure 9a-c, the maximum estimation error of θ m (t) occurred at 1.19 s with an error of 0.02 rad, and the estimation of θ m (t) by ESO was very accurate because it largely coincided with the true value. The maximum estimation error of . θ m (t) occurred at 1.13 s with an estimation error of 0.64 rad/s. In general, the ESO's estimation to . θ m (t) was relatively accurate. For the total perturbation f (t) estimation, the maximum estimation error of 14.78 rad/s 2 occurred at 1.00 s. The effectiveness of ESO in estimating the total perturbation seemed less satisfactory, but such a result was acceptable because it included the total disturbance of the system, which had larger uncertainty. The estimation error of f (t) stabilized in the range of [−1,1] after 2.77 s and, eventually, became stable as the system converged. Overall, the ESO estimates of the state variables largely overlapped with the system's output response, demonstrating the superiority of ADRC.
was relatively accurate. For the total perturbation ( ) estimation, the maximum estimation error of 14.78 rad/s 2 occurred at 1.00 s. The effectiveness of ESO in estimating the total perturbation seemed less satisfactory, but such a result was acceptable because it included the total disturbance of the system, which had larger uncertainty. The estimation error of ( ) stabilized in the range of [-1,1] after 2.77 s and, eventually, became stable as the system converged. Overall, the ESO estimates of the state variables largely overlapped with the system's output response, demonstrating the superiority of ADRC.

Disturbance Rejection Test
The system's immunity to interference should be mentioned, and it is recommended to consider the step disturbance terms d 1 (t) and d 2 (t) at 11 s, where d 1 (t) = 0.07 rad/s 2 and d 2 (t) = −1.00 m/s 2 [38]. As shown in Figure 8d, the pitch angle maximum occurred at 18.56s at 1.0029 rad with a maximum overshoot percentage of 0.29%, which indicated that the system was always within the acceptable interval. At approximately 18s, the control signal u(t) reached a maximum value of 0.294 rad and, finally, converged to a control signal of 0.293 rad. The subtle changes in pitch angle and control signal proved that the system had relatively good immunity to interference.

Robustness Analysis
In this section, we designed an innovative ADRC robustness testing method, the main idea of which is based on Monte Carlo statistical analysis using Gaussian distribution. Under the 3σ principle, it randomly draws the parameter values and composes 22,000 groups, guaranteeing a 99.73% confidence level in computer simulation. The random range of values for each parameter is given in Table 1. The median of each variable parameter was the mean value (u) of the Gaussian distribution, and one-third of the difference between the mean and the upper or lower bound was used as the standard deviation (σ). Figure 10 shows the Gaussian probability density function curves for the six parameters in Table 1. 2 at 18.56s at 1.0029 rad with a maximum overshoot percentage of 0.29%, which indicated that the system was always within the acceptable interval. At approximately 18s, the control signal ( ) reached a maximum value of 0.294 rad and, finally, converged to a control signal of 0.293 rad. The subtle changes in pitch angle and control signal proved that the system had relatively good immunity to interference.

Robustness Analysis
In this section, we designed an innovative ADRC robustness testing method, the main idea of which is based on Monte Carlo statistical analysis using Gaussian distribution. Under the 3σ principle, it randomly draws the parameter values and composes 22,000 groups, guaranteeing a 99.73% confidence level in computer simulation. The random range of values for each parameter is given in Table 1. The median of each variable parameter was the mean value ( ) of the Gaussian distribution, and one-third of the difference between the mean and the upper or lower bound was used as the standard deviation ( ). Figure 10 shows the Gaussian probability density function curves for the six parameters in Table 1. Simulations performed for each set of parameters, and the outputs were analyzed. Three indices were planned to evaluate the system's output: setting time (TS), maximum percentage overshoot (MPO), and integral absolute error (IAE). Of the 22,000 experiments perform, 364 had parameters outside the upper and lower bounds of random values (3σ principle). Therefore, only 21,636 sets of data were available. Figure 11a,b show the simulation results for a total of 18,000 simulations, of which 17,706 sets of parameters satisfied the 3σ principal, and Figure 11cd are the convergence curves saved at the 2000th and 8000th simulations, respectively, during the experimental process. From Figure 11a, there was a clear stratification of the setting time, which mainly occurred at 6.64 s, 8.92 s and 11.32 s. The second stratification occurred at 8.92 s, and the Simulations performed for each set of parameters, and the outputs were analyzed. Three indices were planned to evaluate the system's output: setting time (TS), maximum percentage overshoot (MPO), and integral absolute error (IAE). Of the 22,000 experiments perform, 364 had parameters outside the upper and lower bounds of random values (3σ principle). Therefore, only 21,636 sets of data were available. Figure 11a,b show the simulation results for a total of 18,000 simulations, of which 17,706 sets of parameters satisfied the 3σ principal, and Figure 11c,d are the convergence curves saved at the 2000th and 8000th simulations, respectively, during the experimental process. From Figure 11a, there was a clear stratification of the setting time, which mainly occurred at 6.64 s, 8.92 s and 11.32 s. The second stratification occurred at 8.92 s, and the simulations with the setting time of less than 8.92 s accounted for 79.49% of the valid simulation. In addition, the amounts of simulation with a maximum overshoot percentage of less than 35% accounted for 78.59%, and the number when the integral absolute error was less than 1200 was 99.28%. Overall, the control system showed strong robustness in response to parameter changes. Based on the experimental results, Table 2 counts three percentages under different simulation times. The first one was the proportion of valid simulations (parameters satisfying the 3σ principal). The second one was the ratio of the number of simulations with the setting time less than 8.92s to the number of valid simulations. As can be seen from Figure 11b, there was a strong linear relationship between the maximum percent Based on the experimental results, Table 2 counts three percentages under different simulation times. The first one was the proportion of valid simulations (parameters satisfying the 3σ principal). The second one was the ratio of the number of simulations with the setting time less than 8.92s to the number of valid simulations. As can be seen from Figure 11b, there was a strong linear relationship between the maximum percent overshoot and the integral absolute error, and we chose one of them as the third percentage in Table 2, which was the ratio of the number of simulations with the maximum percent overshoot less than 35% to the valid number of simulations. It can be seen from Table 2 that when the number of experiments reached more than 12,000 times, the changes of the three percentages were less than 0.1%. Therefore, for a 6-parameter Monte Carlo shooting test, the number of experiments was maintained at more than 12,000 times, and the number of experiments could be considered sufficient. The first column counts the total number of simulations, the second, third, and fourth columns count three percentages under different simulation times, valid simulation is the ratio of the amounts of valid simulations to the total simulation, and the third column counts the ratio of the number of simulations with the setting time less than 8.92s to the number of valid simulations. The fourth column counts the ratio of the number of simulations with the maximum percent overshoot less than 35% to the valid number of simulations. When the number of experiments is greater than 12,000, there are subtle changes in the three percentages (less than 0.1%).

Discussion
The constant vibration frequency was the most significant deficiency of this study. In actual flight missions, the vibration frequency is variable, due to unknown external and internal disturbances. Therefore, it is recommended that deep learning algorithms, such as critical Recurrent Neural Networks (RNN) or Long Short-Term Memory (LSTM) neural networks, be considered in future studies to monitor and predict the system frequency changes at each moment, enabling ESO to achieve more accurate disturbance prediction. However, the combination of ADRC with deep learning algorithms can lose its usefulness to some extent.

Conclusions
This paper established a pitch angle vibration model of the flexible aircraft using Matlab/Simulink. The frequency domain characteristics of the system were analyzed comprehensively through a series of engineering analysis methods, while the design process of second-order ADRC was discussed in detail. Subsequently, an empirical summary of the ADRC parameter tuning process was proposed. In the presence of disturbances and uncertainties, the aircraft vibration model was controlled by two different robust control techniques, classical PID and ADRC techniques, respectively. Furthermore, the Monte Carlo shooting method was proposed to verify the robustness of the closed-loop control system. It was clarified that when the system parameters change in significant ranges, the designed controller showed a robust system within a 99.73% confidence interval by considering three indicators (MPO, TS, and IAE). In the test of referring to six dynamic parameters, more than 12,000 Monte Carlo simulations could reflect the statistical characteristics of the results. The structure of ADRC is also relatively simple compared with other modern control theories, and its control effect is better than the PID algorithm. The characteristics of robust control give ADRC the potential to evolve from simulation to practical engineering application.