Gust Load Alleviation with Robust Control for a Flexible Wing

Traditional methods for gust alleviation of aircraft are mostly proposed based on a specific flight condition. In this paper, robust control laws are designed for a large flexible wing with uncertainty in Mach number and dynamic pressure. To accurately describe the aeroelastic model over a large flight envelope, a nonlinear parameter-varying model is developed which is a function of both Mach number and dynamic pressure.Then a linear fractional transformation is established accordingly and amodifiedmodel order reduction technique is applied to reduce the size of the uncertainty block.The developed model, in which the statistic nature of the gust is considered by using the Dryden power spectral density function, enables the use of μ-synthesis procedures for controller design. The simulations show that the μ controller can always effectively reduce the wing root shear force and bending moment at a given range of Mach number and dynamic pressure.


Introduction
Gust load alleviation (GLA) design is used to improve passenger comfort and reduce the dynamic load at given position.This technique involves closely coupled interactions of unsteady aerodynamics, dynamic qualities of flexible structures, and control action [1].A historical review of the research efforts on the problem is given by Capello et al. [2].Different approaches including linear quadratic regulator theory [3,4], optimal control algorithms [5], and  ∞ robust control [6] have been investigated in the last years.Recently, the robust controllers have been extensively studied to account for model uncertainties, including variations in the nominal flight parameters and errors caused by model order reduction.Idan and Shaviv [7] designed a  controller for an aeroservoelastic (ASE) model with uncertainty in mass and proposed an order reduction technique of the uncertainty block.Blue and Balas [8] developed a linear parametervarying (LPV) model to facilitate the synthesis and analysis of flutter suppression controllers.This LPV model has linear dependence on Mach number and dynamic pressure and proves to have enough accuracy when representing a rigid wing model through a large flight envelope, but this accuracy cannot be guaranteed when elastic modes are taken into consideration.Moulin [9] designed a flutter suppression  controller for an aeroelastic plant incorporating airspeed and air density variations.Qian et al. [10] presented a new scheme to model the uncertainty block with reduced order, which is convenient for the design of robust controller.Both Idan and Qian's models do not account for the effect of Mach number on generalized aerodynamic force, which is essential for a flexible vehicle flying at subsonic speed.
Some researches aiming at improving the traditional control methods also proved promising.For example, Fonte et al. [11] proposed a GLA method based on static output feedback and gave a scheduled solution to assure adequate margins over the whole flight envelope; Alam et al. [12] developed a mixed feedforward/feedback approach to improve the robust performance of the feedforward GLA system.Besides, an improved model-predictive control formulation for stabilization and gust alleviation was proposed by Haghighat et al. [13].Application of intelligent material systems such as piezoelectric actuators in GLA system is also of interest for researchers in recent years [3,14].
The objective of this paper is to present a new algorithmic scheme to design a robust GLA controller for a flexible multiple-actuated wing (MAW) with uncertainty in Mach number and dynamic pressure.A nonlinear parametervarying (NPV) model is developed first to accurately represent the dependence of the aeroelastic system on Mach number and dynamic pressure with much smaller approximation error than the LPV model.Based on the structure of the NPV model, a state-space linear fractional transformation (LFT) model is constructed in order to design the robust controller.To facilitate the  controller design procedure, a revised order reduction technique based on the work of Idan and Shaviv [7] and Wang et al. [15] is applied to the uncertainly block of the LFT model.

Linear Time-Invariant (LTI) Models of Aeroservoelastic (ASE) System
This section sums up some necessary fundamentals of LTI models of ASE system.Doublet lattice method is used to calculate the three-dimensional aerodynamic forces for each flight condition.
Constructing the open-loop LTI model at a given flight condition is the starting point of the controller design process, in which the unsteady aerodynamic matrices can be obtained through rational function approximation by the minimum state method, which transforms frequency domain generalized forces to time domain [16]: where  = /,  is the Laplace variable,  is half of the reference chord length,  is the flight velocity, Q 0 , Q 1 , Q 2 , D, and E are all real matrix, and R is a diagonal matrix composed of real aerodynamic lag roots.Q  (), Q  (), and Q  () are generalized aerodynamic force matrices corresponding to the structural modes, control surfaces, and gust, respectively.To eliminate the coefficients associated with the second-order derivative of the gust velocity, approximation constraints are applied to the gust column to yield Q 2 () = 0 in (2).It is important to note that the results of rational function approximation are only valid for the given Mach number at which the frequency domain generalized forces are calculated.
The open-loop aeroelastic system with gust excitation for a fixed flight condition can be written in the state-space form [17]: and it can be rewritten as where where  denotes the -dimensional generalized modal coordinates,   is the   -dimensional aerodynamic lag states, M  and M  are generalized mass matrices, K  is the generalized stiff matrix, C  is the generalized damping matrix,  signifies the deflections of control surface, and   is gust velocity and  is air density.
For an acceleration sensor, the output equation can be written in the state-space form and for a displacement sensor, the output equation can be written as where Φ is the modal matrix of the structural grid nodes.The state-space form of actuator can be expressed as where x ac is the state vector of actuator, u  is the input signal, A ac is the dynamic matrix of actuator, B ac is the input gain matrix, and C ac is the output gain matrix.These matrixes can be obtained from the transfer function of actuator.
The state-space realization of the gust model is where x  is the gust state vector,  is white noise, A  is the dynamic matrix, B  is the input gain matrix, C  is the output gain matrix, and D  is the direct transfer matrix, all of which are derived from the power spectral density (PSD) function of Dryden form [18].
To design a gust GLA controller, the gust load, which can be calculated by modal displacement method, must be included in the output of the ASE model.The load equation is given by where C  is load coefficients matrix.
The open-loop state-space equation of an elastic aircraft is obtained by introducing the state-space equation of actuator and gust or in a more compact form where The output equation should be in the form where matrixes C and D can be easily obtained from ( 5), ( 6), and ( 9) when the output parameters are determined.

Nonlinear Parameter-Varying (NPV) Model of the ASE System
In this section, LTI models at some selected flight conditions are used to develop a NPV model of the ASE system.Mach number (Ma) and dynamic pressure () are used as parameters of the NPV model.The dependence of each element in the LTI state-space matrixes A, B, C, and D on  and Ma had been extensively studied to determine a basis function for the NPV model that can capture this dependence.
It can be easily observed from ( 2), (3), and ( 10) that part of the elements in matrixes A and B are affected by Ma, , and , among which the influence of Ma is indirectly reflected by the generalized aerodynamic forces.But for the case of standard atmosphere, only two of the parameters Ma, , and  are independent.In this paper, Ma and  are adopted as the independent uncertain parameters for convenience, so (11) can be rewritten as where where Δ and ΔMa represent the maximal variation of  and Ma and  1 and  2 represent the normalized uncertainty of dynamic pressure and Mach number within the range of In (11), the dependence of matrix elements on Ma and  can be approximated by a LPV model [8] and similarly for C(Ma, ) and D(Ma, ).While the elements of C(Ma, ) and D(Ma, ) corresponding to output y do depend on  and Ma, the output y is still linear combination of state variable x and input variable u.Thus, C(Ma, ) and D(Ma, ) can be easily constructed from A(Ma, ) and B(Ma, ).The accuracy of the LPV model was verified for a rigid wing by Blue and Balas [8].But later we found out that large approximation errors were introduced for the case of the large flexible wing, and it is supposed that this large approximation error is caused by the elastic modes whose unsteady aerodynamics destroys the linear characteristics of the matrix elements.
Through polynomial regression analysis we found that most of the elements in A(Ma, ) and B(Ma, ) showed a quadratic relationship with  and a linear relationship with Ma, so a new NPV model is built up here, and its accuracy is verified in Section 5.
This NPV model, in which flight conditions within a given range are all considered, significantly simplifies control design and analysis.

Linear Fractional Transformation (LFT) Model and Order Reduction
In this section, a reduced order NPV model is written as LFT form on the parameters  and Ma, which is required to utilize robust control theory, such as  analysis and synthesis technique applied in this paper.After checking ( 2) and ( 10) in detail, it can be found that only the rows corresponding to ξ and ẋ  (i.e., rows with subscripts from  + 1 to 2 +   ) have elements affected by Ma and , so these rows can be rewritten as Now, define the following new outputs z and inputs v, and then it is obvious that where Substituting ( 19) into (18) gives and then substituting ( 22) into ( 11) and combining with (19) gives and details of the matrix in (23) can be easily got and are omitted here for simplicity.The output equation is given by if the output is node acceleration, the matrix in (24) can be obtained through substituting ( 5) into (23); if the output is node displacement, then D  and D V should be null matrix while C  = [Φ 0 0] in which Φ corresponds to the generalized modal coordinates ; if the output is gust load at given position, the matrix in (24) can be obtained through substituting ( 9) into (23).Actually in the robust control model, the node acceleration or displacement should be exported as controlled output and input of controller while the gust load should be exported as controlled output, so the output y in (24) is a combined output.Combining ( 23) and (24) would produce and then define P() as the Laplace transform of S, so and then, closing the loop from z to v with Δ, given by  and is shown in Figure 1.
For controller synthesis and analysis it is necessary to use normalized parameters that vary in the range ±1.This transformation can be easily performed using an LFT with a scaling matrix [19].After normalizing the parameters the new LFT is denoted as   (P 1 , Δ 1 ), in which As the number of structural modes increases, the order of the uncertainty block Δ will increase accordingly.Furthermore, because the NPV model has more fitting items than the LPV model, this will also magnify the order of Δ.To reduce the dimension of the uncertainty block the technique of Idan and Shaviv [7] is used here.The LFT is rewritten by considering x as an uncertainty input and ẋ as an uncertainty output.Then the uncertainty block  1 I corresponding to  is separated from Δ 1 , considering its uncertainty inputs as "pseudoderivatives" and its uncertainty outputs as "pseudostates."Modes that are uncontrollable or unobservable with respect to the transform variable  1 can be canceled, and the size of the corresponding block  1 I can be reduced accordingly.
To further reduce the size of the uncertainty block, the least controllable and observable states are also truncated similar to the balanced truncation procedure.But for the current model, the system matrix is a null matrix, which can be concluded from (19).While the balance truncation procedure needs the system matrix to have all eigenvalues with negative real part, this can be fixed by adding a negative definite matrix − 1 I to the system matrix with  1 greater than 0 but much smaller than 1.Considering ( 20) and (23), it is equivalent to where v 1 and z 1 correspond to the first part in (20) corresponding to  with order 3(+  ).Equation ( 30) is equivalent to because | 1 | < 1 and 0 <  1 ≪ 1; the introduction of − 1 I almost has no influence on the bound of Δ 1 .The previous procedure should be repeated for the uncertainty block corresponding to Ma.

𝜇 Controller Design of a Multiple-Actuated
Wing (MAW) Model The MAW unsteady aerodynamic model is shown in Figure 3 with the locations of accelerometers.The model contains two trailing control surfaces close to mid span of the wing, and both are modeled by panel boxes mapped to structural nodes through spline interpolation.It should be noticed that the root of the MAW structural model is fixed to study the effects of elastic modes only.

NPV Model and Robust Control
Model.Dryden gust model is applied when constructing the LTI model for each flight condition.The root-mean-square value of the vertical gust velocity is 1.5 m/s and the scale of turbulence is 760 m.The vertical displacement or acceleration of node 1 and node 2 shown in Figure 3 is set as sensor outputs.The transfer functions of the actuators are both chosen as [20]  ac = 3.302 × 10 5  3 + 127.2 2 + 8789 + 3.302 × 10 5 . (32) A total of 24 flight conditions listed in Table 1 were used to build the NPV model, and these flight conditions have The italic numbers mean that these flight conditions are unstable.The italic numbers mean that these results are unacceptable.
covered the flutter dynamic pressure at each considered Mach number.Open-loop matched flutter analysis at Mach number from 0.4 to 0.7 was conducted through both the frequencydomain -method [21] and time domain root locus method.Fifteen generalized aerodynamic force matrices at reduced frequency values between  = 0.0 and  = 2.6 were calculated at each Ma before rational function approximation.Eight low-frequency elastic modes were used in the analysis.Structural damping of 0.01 was defined for all frequencies.The LTI model of every flight condition has 29 states including sixteen structural states, four aerodynamic states, six actuator states, and three gust states.The flutter analysis results are shown in Table 2.
The accuracy of the LPV model and NPV model was verified by calculating the flutter characteristics at the Mach number given in Table 1 with the time domain root locus method.The results are also shown in Table 2.It is obvious that the LPV model failed to approximate the original LTI model at Mach numbers 0.4 and 0.5 while the NPV model reached high accuracy.The accuracy of the NPV model was also proved by comparing the bode plots from the inner control surface to acceleration of node 1 at given flight conditions.Two cases are shown in Figure 4 in which the flight conditions are chosen at random from Table 1.Thus, it was concluded that the single NPV model adequately described the dynamics of the MAW for the purpose of controller synthesis and analysis.
The objective of designing the  controller is to suppress the vibration and to alleviate the gust load for the openloop aeroelastic system.Additionally, the controller should be effective when Ma varies between 0.4 and 0.7, and  varies between 10 kPa and 40 kPa.As a result, the values of Ma 0 , ΔMa,  0 , and Δ in (15) are set to 0.55, 0.15, 25, and 15 separately.
Then it is necessary to build up the robust control model.Figure 5 shows the structure of the closed-loop robust control model, in which K is the robust controller to be designed.P 1 is the nominal plant model after normalizing the uncertainty block, and Δ 1 is the normalized uncertainty block representing the dynamic pressure uncertainty and Mach number uncertainty.W  is the output weighting function.For the current model, the output contains the displacement of the measurement nodes and the wing root gust load.The displacement of the measurement nodes is chosen as input to the robust controller for the purpose of suppressing the vibration of low order modes which have main contribution to the wing root load.Because the node displacement is difficult to measure directly, integration of acceleration signal will be applied in the simulation.
The output weighting function of the displacement of node 1 and node 2 is where  1 (with value 0.1) and  2 (with value 0.27) are the maximum value of the transfer functions (among all the flight conditions within the considered range) from gust input to displacement of node 1 and node 2, respectively.These weights correspond to asking for a decrease to 30% of the open-loop response at the corresponding flight condition.The weight on the wing root shear force and wing root bending moment is where  1 (with value 3.3) and  2 (with value 13700) are calculated in the same way as  1 and  2 .Then  1 is set to 50 to put more weight in the low-frequency range, which contains most of the gust energy.
The actuator has a maximum deflection of ±15 deg or /12 rad.High frequency control commands should also be avoided since the servo dynamics have a bandwidth of 22 rad/s.Thus, the weight on the control command is chosen as In addition, it is assumed that the displacement measurements include a white, zero mean measurement noise   . is set to 0.01 to normalize this input signal.
The multiplicative unstructured output uncertainty is represented by a 2 × 2 complex uncertainty block Δ out and the weighting function is and it specifies a relative error of 10% at 20 rad/s and 85% at 1000 rad/s.Ma = 0.45, q = 40 kPa Ma = 0.65, q = 30 kPa The model in Figure 5 can be transformed to a structure suitable for controller design by the command sysic in Matlab robust control toolbox [22].The order of Δ 1 is reduced from 60 to 16 by using the technique in Section 4. The design process converged after 2 - iterations resulting in a controller of 189th order, which was order-reduced through balanced truncation method to 22.The final value of  corresponding to the order-reduced controller is 0.954, which is smaller than 1.That means with the robust controller the robust stability and robust performance of the closed-loop system are guaranteed.

Closed-Loop Performance Analysis.
The closed-loop performance is analyzed both in frequency domain and in time domain.In Figures 6 and 7, the Bode plots from gust disturbance to wing root load (including wing root shear force and bending moment) are compared between the open-and closed-loop models for two representative flight conditions.Clearly, wing root load caused by low-frequency disturbance is effectively alleviated, and this conclusion holds for any flight condition in the considered range.
To verify the performance of the  controller in time domain simulation, the Mach number is set to linearly increase from 0.4 to 0.7 in 20 seconds and the dynamic pressure linearly increases from 10 kPa to 40 kPa simultaneously.This means that the flight condition changes a lot in this process.It is also obvious that both the open-and closed-loop system are stable in this simulation.The open-and closedloop response of wing root load to the Dryden gust is shown in Figure 8.The RMS values of the wing root shear force and bending moment decrease by 28 and 34 percent, respectively, after active control.Furthermore, the motion of the two control surfaces is shown in Figure 9.It can be seen that the deflection angle remained within ±10 ∘ in the simulation.

Conclusions
The technique for modeling ASE systems perturbed by significant variations in Mach number and dynamic pressure is developed in this study.The generated LFT model is convenient for robust stability and performance analysis as well as for design of robust controllers of these systems.Structured uncertainties were used to represent the perturbations of Mach number and dynamic pressure and the order of the structural uncertainty block was reduced effectively.The developed technique was applied for design of GLA controller for a MAW in a wide range of flight conditions.The simulation results demonstrated efficiency of the presented methodology and tools.

Figure 4 : 1 Figure 5 :
Figure 4: Bode plots from the inner control surface to acceleration of node 1 at two given flight conditions.

FrequencyFigure 6 :
Figure 6: Open-and closed-loop PSD functions from wind gust to wing root load at Ma = 0.4 and  = 10 kPa.

Figure 7 :
Figure 7: Open-and closed-loop PSD functions from wind gust to wing root load at Ma = 0.7 and  = 40 kPa.

Figure 8 :Figure 9 :
Figure 8: Time response of wing root shear force and bending moment to wind-gust input as Ma and  vary.

Table 1 :
Flight conditions of LTI models.