# Gust Load Alleviation with Robust Control for a Flexible Wing

Authors: Xiang Liu; Qin Sun
Journal: Shock and Vibration
 (2016)

---

## Abstract

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 a modified model 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.


---

## Body

## 1. 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 H
∞ 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 parameter-varying (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 parameter-varying (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.


## 2. 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]:(1)
Q
s
-
≜
Q
ξ
ξ
s
-
Q
ξ
δ
s
-
Q
g
s
-
=
Q
0
+
Q
1
s
-
+
Q
2
s
-
2
+
D
s
-
I
-
R
-
1
E
s
-
,where s
-
=
s
b
/
V, s is the Laplace variable, b is half of the reference chord length, V 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
ξ
ξ
(
s
-
), Q
ξ
δ
(
s
-
), and Q
g
(
s
-
) 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
g
2
(
s
-
)
=
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]:(2)
ξ
˙
ξ
¨
x
˙
a
=
0
I
0
-
M
-
ξ
-
1
K
-
ξ
-
M
-
ξ
-
1
C
-
ξ
-
1
2
ρ
V
2
M
-
ξ
-
1
D
0
E
ξ
V
b
R
ξ
ξ
˙
x
a
+
0
0
0
-
1
2
ρ
V
2
M
-
ξ
-
1
Q
ξ
δ
0
-
1
2
ρ
V
b
M
-
ξ
-
1
Q
ξ
δ
1
-
M
-
ξ
-
1
M
ξ
δ
+
1
2
ρ
b
2
Q
ξ
δ
2
0
E
δ
0
δ
δ
˙
δ
¨
+
0
0
-
1
2
ρ
V
M
-
ξ
-
1
Q
g
0
-
1
2
ρ
b
M
-
ξ
-
1
Q
g
1
0
E
g
V
w
g
w
˙
g
,and it can be rewritten as(3)
x
˙
a
e
=
A
a
e
x
a
e
+
B
a
e
δ
a
e
+
E
a
e
w
a
e
,where(4)
x
a
e
=
ξ
T
ξ
˙
T
x
a
T
T
,
δ
a
e
=
δ
T
δ
˙
T
δ
¨
T
T
,
w
a
e
=
w
g
T
w
˙
g
T
T
,
M
-
ξ
=
M
ξ
ξ
+
1
2
ρ
b
2
Q
ξ
ξ
2
,
C
-
ξ
=
C
ξ
ξ
+
1
2
ρ
V
b
Q
ξ
ξ
1
,
K
-
ξ
=
K
ξ
ξ
+
1
2
ρ
V
2
Q
ξ
ξ
0
,where ξ denotes the n-dimensional generalized modal coordinates, x
a is the n
a-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 w
g is gust velocity and ρ is air density.For an acceleration sensor, the output equation can be written in the state-space form(5)
y
a
=
Φ
ξ
¨
,and for a displacement sensor, the output equation can be written as(6)
y
a
e
=
Φ
ξ
,where Φ is the modal matrix of the structural grid nodes.The state-space form of actuator can be expressed as(7)
x
˙
a
c
=
A
a
c
x
a
c
+
B
a
c
u
c
,
δ
a
e
=
C
a
c
x
a
c
,where x
a
c is the state vector of actuator, u
c is the input signal, A
a
c is the dynamic matrix of actuator, B
a
c is the input gain matrix, and C
a
c 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(8)
x
˙
g
=
A
g
x
g
+
B
g
η
,
w
a
e
=
C
g
x
g
+
D
g
η
,where x
g is the gust state vector, η is white noise, A
g is the dynamic matrix, B
g is the input gain matrix, C
g is the output gain matrix, and D
g 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(9)
y
l
=
C
L
ξ
,where C
L 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(10)
x
˙
a
e
x
˙
a
c
x
˙
g
=
A
a
e
B
a
e
C
a
c
E
a
e
C
g
0
A
a
c
0
0
0
A
g
x
a
e
x
a
c
x
g
+
0
B
a
c
0
u
c
+
E
a
e
D
g
0
B
g
η
,or in a more compact form(11)
x
˙
=
A
x
+
B
u
,where(12)
B
=
0
E
a
e
D
g
B
a
c
0
0
B
g
,
u
=
u
c
η
.The output equation should be in the form(13)
y
=
C
x
+
D
u
,where matrixes C and D can be easily obtained from (5), (6), and (9) when the output parameters are determined.


## 3. 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 (M
a) and dynamic pressure (q) 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 q 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 M
a, V, 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 M
a, V, and ρ are independent. In this paper, M
a and q are adopted as the independent uncertain parameters for convenience, so (11) can be rewritten as(14)
x
˙
=
A
M
a
,
q
x
+
B
M
a
,
q
u
,where(15)
q
=
q
0
+
Δ
q
·
δ
1
,
M
a
=
M
a
0
+
Δ
M
a
·
δ
2
,where Δ
q and Δ
M
a represent the maximal variation of q and M
a and δ
1 and δ
2 represent the normalized uncertainty of dynamic pressure and Mach number within the range of [
-
1,1
].In (11), the dependence of matrix elements on M
a and q can be approximated by a LPV model [8](16)
A
M
a
,
q
=
A
0
+
A
q
·
q
+
A
M
a
·
M
a
+
A
q
M
a
·
q
·
M
a
,
B
M
a
,
q
=
B
0
+
B
q
·
q
+
B
M
a
·
M
a
+
B
q
M
a
·
q
·
M
a
,and similarly for C
(
M
a
,
q
) and D
(
M
a
,
q
). While the elements of C
(
M
a
,
q
) and D
(
M
a
,
q
) corresponding to output y do depend on q and M
a, the output y is still linear combination of state variable x and input variable u. Thus, C
(
M
a
,
q
) and D
(
M
a
,
q
) can be easily constructed from A
(
M
a
,
q
) and B
(
M
a
,
q
). 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 inA
(
M
a
,
q
) and B
(
M
a
,
q
) showed a quadratic relationship with q and a linear relationship with M
a, so a new NPV model is built up here,(17)
A
M
a
,
q
=
A
0
+
A
q
·
q
+
A
M
a
·
M
a
+
A
q
M
a
·
q
·
M
a
+
A
q
q
·
q
2
+
A
q
q
M
a
·
q
2
·
M
a
,
B
M
a
,
q
=
B
0
+
B
q
·
q
+
B
M
a
·
M
a
+
B
q
M
a
·
q
·
M
a
+
B
q
q
·
q
2
+
B
q
q
M
a
·
q
2
·
M
a
,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.


## 4. Linear Fractional Transformation (LFT) Model and Order Reduction

In this section, a reduced order NPV model is written as LFT form on the parametersq and M
a, 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 x
˙
a (i.e., rows with subscripts from n
+
1 to 2
n
+
n
a) have elements affected by M
a and q, so these rows can be rewritten as(18)
x
˙
i
=
A
0
,
i
x
+
B
0
,
i
u
+
A
q
,
i
x
+
B
q
,
i
u
q
+
A
M
a
,
i
x
+
B
M
a
,
i
u
M
a
+
A
q
M
a
,
i
x
+
B
q
M
a
,
i
u
q
M
a
+
A
q
q
,
i
x
+
B
q
q
,
i
u
q
2
+
A
q
q
M
a
,
i
x
+
B
q
q
M
a
,
i
u
q
2
M
a
i
=
n
+
1
,
n
+
2
,
…
,
2
n
+
n
a
.Now, define the following new outputsz and inputs v,(19)
z
3
i
-
2
=
A
q
q
,
n
+
i
x
+
B
q
q
,
n
+
i
u
+
v
3
n
+
n
a
+
2
i
-
1
,
v
3
i
-
2
=
z
3
i
-
2
·
q
,
z
3
i
-
1
=
A
q
,
n
+
i
x
+
B
q
,
n
+
i
u
+
v
3
i
-
2
,
v
3
i
-
1
=
z
3
i
-
1
·
q
,
z
3
i
=
A
q
M
a
,
n
+
i
x
+
B
q
M
a
,
n
+
i
u
,
v
3
i
=
z
3
i
·
q
,
z
3
n
+
n
a
+
2
i
-
1
=
A
q
q
M
a
,
n
+
i
x
+
B
q
q
M
a
,
n
+
i
u
,
v
3
n
+
n
a
+
2
i
-
1
=
z
3
n
+
n
a
+
2
i
-
1
·
M
a
,
z
3
n
+
n
a
+
2
i
=
A
M
a
,
n
+
i
x
+
B
M
a
,
n
+
i
u
+
v
3
i
,
v
3
n
+
n
a
+
2
i
=
z
3
n
+
n
a
+
2
i
·
M
a
,
i
=
1,2
,
…
,
n
+
n
a
,and then it is obvious that(20)
v
=
Δ
·
z
,where(21)
Δ
=
q
·
I
3
n
+
n
a
M
a
·
I
2
n
+
n
a
.Substituting (19) into (18) gives(22)
x
˙
n
+
i
=
A
0
,
n
+
i
x
+
B
0
,
n
+
i
u
+
v
3
i
-
1
+
v
3
n
+
n
a
+
2
i
i
=
1,2
,
…
,
n
+
n
a
,and then substituting (22) into (11) and combining with (19) gives(23)
x
˙
z
=
A
B
u
B
v
C
z
D
z
u
D
z
v
x
u
v
,and details of the matrix in (23) can be easily got and are omitted here for simplicity.The output equation is given by(24)
y
=
C
y
x
+
D
y
u
u
+
D
y
v
v
;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
y
u and D
y
v should be null matrix while C
y
=
Φ
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(25)
x
˙
y
z
=
A
B
u
B
v
C
y
D
y
u
D
y
v
C
z
D
z
u
D
z
v
x
u
v
≜
S
x
u
v
,and then define P
(
s
) as the Laplace transform of S, so(26)
y
z
=
P
s
u
v
,and then, closing the loop from z to v with Δ, given by(27)
y
=
F
l
P
,
Δ
u
,yields a model equivalent to the NPV model. In (27), F
l
(
P
,
Δ
) represents a lower linear fractional transformation defined as(28)
F
l
P
,
Δ
=
P
11
+
P
12
Δ
I
-
P
22
Δ
-
1
M
12and is shown in Figure 1.Figure 1
Structure ofF
l
(
P
,
Δ
).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 F
l
(
P
1
,
Δ
1
), in which(29)
Δ
1
=
δ
1
·
I
3
n
+
n
a
δ
2
·
I
2
n
+
n
a
.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 x
˙ as an uncertainty output. Then the uncertainty block δ
1
I corresponding to q 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 -
a
1
I to the system matrix with a
1 greater than 0 but much smaller than 1. Considering (20) and (23), it is equivalent to(30)
v
1
=
δ
1
I
z
1
-
a
1
v
1
,where v
1 and z
1 correspond to the first part in (20) corresponding to q with order 3
(
n
+
n
a
). Equation (30) is equivalent to (31)
v
1
=
δ
1
1
+
a
1
δ
1
I
·
z
1
,because δ
1
<
1 and 0
<
a
1
≪
1; the introduction of -
a
1
I almost has no influence on the bound of Δ
1.The previous procedure should be repeated for the uncertainty block corresponding toM
a.


## 5.μ Controller Design of a Multiple-Actuated Wing (MAW) Model

### 5.1. Test-Case MAW Model

The technique developed for modeling of ASE systems with Mach number and dynamic pressure uncertainty was applied to design the feedback control system for GLA. A model of a MAW of a generic unmanned aerial vehicle was used. The MSC/NASTRAN structural model is shown in Figure2. It contains 609 grid nodes and 779 elements. The weight of the model is 120 kg. The span of the wing is 8.3 meters and the chord length varies from 2.1 meters at the wing root to 0.5 meters at the wing tip. 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.Figure 2
MAW structural model.Figure 3
MAW aerodynamic model.


### 5.2. 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 Figure3 is set as sensor outputs. The transfer functions of the actuators are both chosen as [20](32)
G
a
c
=
3.302
×
1
0
5
s
3
+
127.2
s
2
+
8789
s
+
3.302
×
1
0
5
.A total of 24 flight conditions listed in Table1 were used to build the NPV model, and these flight conditions have 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 frequency-domain g-method [21] and time domain root locus method. Fifteen generalized aerodynamic force matrices at reduced frequency values between k
=
0.0 and k
=
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.Table 1
Flight conditions of LTI models.
Ma
Dynamic pressureq (kPa)
0.4
10
20
30
40
50
60
0.5
10
20
30
40
50
60
0.6
10
20
30
40
50
60
0.7
10
20
30
40
50
60
The italic numbers mean that these flight conditions are unstable.Table 2
Open-loop flutter parameters.
Ma
g-method
LTI model
LPV model
NPV model
V
f, m/s
f, Hz
V
f, m/s
f, Hz
V
f, m/s
f, Hz
V
f, m/s
f, Hz
0.4
154.3
30.22
154.8
31.69
121.7
100.67
152.3
31.52
0.5
183.3
28.52
185.6
28.82
152.1
100.51
185.5
28.85
0.6
211.0
26.39
214.5
26.12
214.7
26.49
214.5
26.13
0.7
238.1
23.93
242.4
23.97
247.9
24.51
241.8
23.15
The italic numbers mean that these results are unacceptable.The accuracy of the LPV model and NPV model was verified by calculating the flutter characteristics at the Mach number given in Table1 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.Figure 4
Bode plots from the inner control surface to acceleration of node 1 at two given flight conditions.The objective of designing theμ controller is to suppress the vibration and to alleviate the gust load for the open-loop aeroelastic system. Additionally, the controller should be effective when M
a varies between 0.4 and 0.7, and q varies between 10 kPa and 40 kPa. As a result, the values of M
a
0, Δ
M
a, q
0, and Δ
q in (15) are set to 0.55, 0.15, 25, and 15 separately.Then it is necessary to build up the robust control model. Figure5 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
y 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.Figure 5
Structure of the closed-loop robust control model.The output weighting function of the displacement of node 1 and node 2 is(33)
W
y
=
1
0.7
k
1
1
0.7
k
2
,where k
1 (with value 0.1) and k
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(34)
W
l
=
1
0.7
l
1
1
s
/
w
1
+
1
1
0.7
l
2
1
s
/
w
1
+
1
,where l
1 (with value 3.3) and l
2 (with value 13700) are calculated in the same way as k
1 and k
2. Then w
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(35)
W
u
=
12
π
s
/
20
+
1
s
/
50
+
1
12
π
s
/
20
+
1
s
/
50
+
1
.In addition, it is assumed that the displacement measurements include a white, zero mean measurement noisew
n. ε is set to 0.01 to normalize this input signal.The multiplicative unstructured output uncertainty is represented by a2
×
2 complex uncertainty block Δ
o
u
t and the weighting function is(36)
W
o
u
t
=
0.05
s
/
20
+
1
s
/
500
+
1
,and it specifies a relative error of 10% at 20 rad/s and 85% at 1000 rad/s.The model in Figure5 can be transformed to a structure suitable for controller design by the commandsysic 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 D-K 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.


### 5.3. Closed-Loop Performance Analysis

The closed-loop performance is analyzed both in frequency domain and in time domain. In Figures6 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.Figure 6
Open- and closed-loop PSD functions from wind gust to wing root load atM
a
=
0.4 and q
=
10 kPa.Figure 7
Open- and closed-loop PSD functions from wind gust to wing root load atM
a
=
0.7 and q
=
40 kPa.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 closed-loop response of wing root load to the Dryden gust is shown in Figure 8.Figure 8
Time response of wing root shear force and bending moment to wind-gust input as Ma andq vary.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 Figure9. It can be seen that the deflection angle remained within ±
10
° in the simulation.Figure 9
Time response of control surfaces deflection to wind-gust input as Ma andq vary.


## 5.1. Test-Case MAW Model

The technique developed for modeling of ASE systems with Mach number and dynamic pressure uncertainty was applied to design the feedback control system for GLA. A model of a MAW of a generic unmanned aerial vehicle was used. The MSC/NASTRAN structural model is shown in Figure2. It contains 609 grid nodes and 779 elements. The weight of the model is 120 kg. The span of the wing is 8.3 meters and the chord length varies from 2.1 meters at the wing root to 0.5 meters at the wing tip. 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.Figure 2
MAW structural model.Figure 3
MAW aerodynamic model.


## 5.2. 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 Figure3 is set as sensor outputs. The transfer functions of the actuators are both chosen as [20](32)
G
a
c
=
3.302
×
1
0
5
s
3
+
127.2
s
2
+
8789
s
+
3.302
×
1
0
5
.A total of 24 flight conditions listed in Table1 were used to build the NPV model, and these flight conditions have 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 frequency-domain g-method [21] and time domain root locus method. Fifteen generalized aerodynamic force matrices at reduced frequency values between k
=
0.0 and k
=
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.Table 1
Flight conditions of LTI models.
Ma
Dynamic pressureq (kPa)
0.4
10
20
30
40
50
60
0.5
10
20
30
40
50
60
0.6
10
20
30
40
50
60
0.7
10
20
30
40
50
60
The italic numbers mean that these flight conditions are unstable.Table 2
Open-loop flutter parameters.
Ma
g-method
LTI model
LPV model
NPV model
V
f, m/s
f, Hz
V
f, m/s
f, Hz
V
f, m/s
f, Hz
V
f, m/s
f, Hz
0.4
154.3
30.22
154.8
31.69
121.7
100.67
152.3
31.52
0.5
183.3
28.52
185.6
28.82
152.1
100.51
185.5
28.85
0.6
211.0
26.39
214.5
26.12
214.7
26.49
214.5
26.13
0.7
238.1
23.93
242.4
23.97
247.9
24.51
241.8
23.15
The italic numbers mean that these results are unacceptable.The accuracy of the LPV model and NPV model was verified by calculating the flutter characteristics at the Mach number given in Table1 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.Figure 4
Bode plots from the inner control surface to acceleration of node 1 at two given flight conditions.The objective of designing theμ controller is to suppress the vibration and to alleviate the gust load for the open-loop aeroelastic system. Additionally, the controller should be effective when M
a varies between 0.4 and 0.7, and q varies between 10 kPa and 40 kPa. As a result, the values of M
a
0, Δ
M
a, q
0, and Δ
q in (15) are set to 0.55, 0.15, 25, and 15 separately.Then it is necessary to build up the robust control model. Figure5 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
y 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.Figure 5
Structure of the closed-loop robust control model.The output weighting function of the displacement of node 1 and node 2 is(33)
W
y
=
1
0.7
k
1
1
0.7
k
2
,where k
1 (with value 0.1) and k
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(34)
W
l
=
1
0.7
l
1
1
s
/
w
1
+
1
1
0.7
l
2
1
s
/
w
1
+
1
,where l
1 (with value 3.3) and l
2 (with value 13700) are calculated in the same way as k
1 and k
2. Then w
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(35)
W
u
=
12
π
s
/
20
+
1
s
/
50
+
1
12
π
s
/
20
+
1
s
/
50
+
1
.In addition, it is assumed that the displacement measurements include a white, zero mean measurement noisew
n. ε is set to 0.01 to normalize this input signal.The multiplicative unstructured output uncertainty is represented by a2
×
2 complex uncertainty block Δ
o
u
t and the weighting function is(36)
W
o
u
t
=
0.05
s
/
20
+
1
s
/
500
+
1
,and it specifies a relative error of 10% at 20 rad/s and 85% at 1000 rad/s.The model in Figure5 can be transformed to a structure suitable for controller design by the commandsysic 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 D-K 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.


## 5.3. Closed-Loop Performance Analysis

The closed-loop performance is analyzed both in frequency domain and in time domain. In Figures6 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.Figure 6
Open- and closed-loop PSD functions from wind gust to wing root load atM
a
=
0.4 and q
=
10 kPa.Figure 7
Open- and closed-loop PSD functions from wind gust to wing root load atM
a
=
0.7 and q
=
40 kPa.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 closed-loop response of wing root load to the Dryden gust is shown in Figure 8.Figure 8
Time response of wing root shear force and bending moment to wind-gust input as Ma andq vary.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 Figure9. It can be seen that the deflection angle remained within ±
10
° in the simulation.Figure 9
Time response of control surfaces deflection to wind-gust input as Ma andq vary.


## 6. 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.



---

*Source: 1060574-2016-01-06.xml*