A, G, Gd State and input matrices

Ac, Bc State and input matrices of transformed system

A*c Matrix of the system zero dynamics

a∞, p∞, ρ∞ Sound speed ,the pressure and air density of the undisturbed flow respectively

B Non-linear restoring moment

C0, C1, C, C3 Constants used in the bounded neural network composite weight matrix

e, r Tracking error and filtered tracking error, respectively

F, G Positive definite diagonal gain matrixes for update laws of ？ and V

g1, g2 Auxiliary saturation gains

Kz, Kv, Γ , Zb, Kd Controller gains

La(t), Ma(t) lifting and aerodynamic moment

p(y, t) Unsteady pressure

V Dimensionless flight speed

vz(t) Downwash velocity normal to the airfoil surface

W, V Ideal neural network interconnection weight matrices

？,

Estimated neural network interconnection weight matrices

w(t) Transverse deflection

x, u System state and input, respectively

Z, ZB Ideal neural networks composite weight matrix and its bound

Estimated neural networks composite weight matrix and mismatch

Auxiliary control input

ξ Dimensionless plunging displacement ξ = h / b

γ Isentropic gas coefficient (γ = 1.4 for dryair)

τ Dimensionless time τ = Ut / b

λ Aerodynamic factor

η(t) Vector of system states for analysis of zero dynamics

tr{？ } Trace of a square matrix defined as the sum of the elements on the matrix main diagonal

||？||F Frobenius norm defined as ||A||F =

< A, B> Inner product of two matrix, defined as tr{B*A}

In recent years, aeroelastic control and flutter suppression of flexible wings have been extensively investigated by numerous researchers. There are two basic problems associated with the aeroelastic instability of lifting surfaces ? the determination of the flutter boundary and of its character, i.e., the identification of the presence of a stable or unstable Limit Cycle Oscillation (LCO) in the proximity of the flutter boundary. Classical flutter analysis is based on the linearized aeroelastic equations, while LCO analysis requires a nonlinear approach [1]. The goal of the control is to expand the flight envelope above the uncontrolled flutter instability speed without weight penalties and eventually convert the catastrophic nature of flutter, associated with an unstable LCO typical of a subcritical Hopf-bifurcation behavior, into benign flutter, which conversely is associated with a stable LCO typical of a supercritical Hopf-bifurcation. A great deal of research activity devoted to the aeroelastic active control and flutter suppression of flight vehicles has been accomplished, e.g., see [2]. The model nonlinearities can help to stabilize the LCO or be detrimental by destabilizing the LCO [3]. The nonlinearities to be included in the aeroelastic model can be structural [4] (i.e., arising from the kinematic equations); physical [5] (i.e., those involving the constitutive equations); or aerodynamic appearing in the unsteady aerodynamic equations [1][3][6]. This issue is discussed in the context of panel flutter in [1][6][7].

A plethora of techniques is available for dealing with the effect of non-linear structural stiffness in the context of subsonic flow; linear control theory, feedback linearizing techniques, adaptive, and robust control techniques have been employed to account for these nonlinearities, e.g., [8]-[19]. Recently, neural-network-based (i.e., model-free) control approaches have been proposed in [20] and [21] to stabilize a nonlinear aeroelastic wing section. However, there is very little work (e.g., [3],[22]) dealing with the aeroelastic vibration suppression for a supersonic wing section in the presence of both structural and aerodynamic nonlinearities.

Motivated by our previous work in [19][21] and [23]-[25], a novel neural network (NN) based robust controller has been designed to asymptotically stabilize a supersonic aeroelastic system with unstructured nonlinear uncertainties. The nonlinearity of the model depends on the plunging distance and pitching angle. If the nonlinearity is known and could be linearly parameterized, then adaptive control is often considered to be the method of choice. In this paper, we assume unstructured uncertainty in the sense that the structure of the system nonlinearity is considered to be unknown. In contrast to existing neural network-based controllers that only achieve practical stability, the novel continuous control design in this paper is able to achieve asymptotic stability of the origin. A three-layer neural network is implemented to approximate the unknown nonlinearity of the system. While adaptive control relies on linear parameterizability of the system nonlinearity and the determination of a regression matrix, the universal approximation property of the NN controller enables approximation of the unstructured nonlinear system in a more suitable way. To compensate for the inevitable NN functional approximation error, an integral of a sliding mode term is introduced. Through a Lyapunov analysis, global asymptotic stability can be obtained for the tracking error in the pitching degree of freedom. Then, based on the fact that the system is minimum phase, the asymptotic stability of the plunging degree of freedom is also guaranteed. Simulation results show that this NN-based robust continuous control design can rapidly suppress the flutter and limit cycle oscillations of the aeroelastic system.

The rest of the paper is organized as follows. In Section II, the aeroelastic system dynamics are introduced. In Section III, the control objective is stated explicitly while zero dynamics of the system is analyzed. The open-loop error system is developed in Section IV to facilitate the subsequent control design while the closed-loop error system is developed in Section V. In Section VI, Lyapunov-based analysis of the stability of the closed-loop system is presented while the simulation results are shown in Section VII. Appropriate conclusions are drawn in Section VIII.

The aeroelastic governing equations of a supersonic wing section with plunging and twisting degrees-of-freedom (graphically represented in Fig. 1), accounting for flap deflections, and constrained by a linear translational spring and a non-linear torsional spring, are given as follows

The dimensionless plunging distance (positive downward) is expressed as ξ (≡ h / b), while α is the pitch angle (positive nose up),

are derivatives with respect to dimensionless time τ = U_t / b, and V = U / bω_α is the dimensionless flight speed. The parameter B represents the non-linear restoring moment and is defined as the ratio between the linear and non-linear stiffness coefficients, thus it measures of the degree of nonlinearity of the system; B > 0 corresponds to hard structural nonlinearities, B = 0 corresponds to a linear model, while B < 0 corresponds to soft structural nonlinearities. In addition, l_α and m_α represent the dimensionless aerodynamic lift and moment with respect to the elastic axis.

In order to account for flap deflections, some modifications need to be made to the non-linear Piston Theory Aerodynamics (PTA) which is used here to produce the aerodynamic loads on the lifting surface. To keep the paper self-contained, a short description of the PTA modified

to account for the flap deflection is presented next. Within the PTA, the unsteady pressure can be defined as follows

where v_z(t) and α_∞ represent the downwash velocity normal to the airfoil surface and the undisturbed speed of sound respectively, and are defined as follows

In the definition of v_z(t),

denotes the upper and lower surfaces, respectively, while U_∞ denotes the air speed of the undisturbed flow. In the expression (3), α_∞, p_∞ and p_∞ denote the pressure and air density of the undisturbed flow, respectively, while γ is the isentropic gas coefficient (γ = 1.4 for dry-air). The transverse deflection w(t) in (3) can be expressed as [26]

where x₀ and x₁ denote the dimensionless location of the elastic axis and of the torsional spring of the flap from the leading edge respectively, while β(t) represents the flap displacement. In the binomial expansion of (PTA), the pressure formula for PTA in the third-order approximation can be obtained by retaining the terms up to and including ( v_z / α_∞ ) as follows [7], [27]-[29]

The aerodynamic correction factor,

is used to correct the PTA to better approximate the pressure at low supersonic flight speed regime. It is important to note that (2) and (5) are only applicable as long as the transformation through contraction and expansion can be consider isentropic, i.e., as long as the induced show losses are negligible (low-intensity waves). For more details, see [1] [5][30]. PTA provides results in excellent accordance with those based on the Euler solution and the CFL3D code [31]. Considering that flow takes place on both the upper and lower surfaces of the airfoil, U⁺_∞= U^?_∞ = U; from (3)-(5), the aerodynamic pressure

difference can be expressed as

Notice that δ_p also accounts for the deflection of the flap β. Here, M = U_∞ / α_∞ is the undisturbed flight Mach number, while q_∞ = ρ _∞U²_∞ / 2 is the undisturbed dynamic pressure as presented in [1] and [3]. The model can be simplified to account only for the nonlinearities associated with α and discarding those associated with β. Even though this is an approximation, the magnitude of the nonlinearities associated with β is much smaller than those associated with α and will thus be omitted in this paper. In addition, it is assumed in the following development that the nonlinear aerodynamic damping in (6), i.e., the terms w³_t, w²_t w_x, and w_t w²_t will be discarded and consequently, the cubic nonlinear aerodynamic term reduces to w³_t only. Although nonlinear damping can be included in the model, this paper only considers linear damping and thus conservative estimates of the flutter speed are expected.

Finally, the nonlinear aerodynamic lifting and moment can be obtained from the integration of the difference of pressure on the upper and lower surfaces of the airfoil

where δp+_{x< bx1} and δp+_{x< bx1} are the aerodynamic pressure difference on the clean airfoil and on the flap. In the governing EOM presented in (1), l_α and m_α denote the counterpart of (7) and (8), which are defined as

Here, μ represent a the dimensionless mass ratio defined as m / 4ρb². Given the definitions above, the governing EOM can be transformed into the following form

where

is a vector of systems states, β(t) is a flap deflection control input, while A, G(z), G_d(z), and Φ (y) are defined as follows

where the explicit definitions for the constants c_i, k_i, ∀ _i = 1, …, 4 as well as p₂ and p₄ are reported in the Appendix.

The explicit control objective of this paper is to design a model-free aeroelastic vibration suppression strategy to guarantee the asymptotic convergence of the pitch angle α using the flap deflection β as a control input. The secondary objective is to ensure that all system states remain bounded at all times during closed-loop operation. It is assumed that the measurable variables available for control implementation are the pitch angle α, pitch angle velocity

plunging displacement ξ and plunging displacement velocity

Since the proposed control strategy is predicated on the assumption that the system of (12) is minimum phase, the stability of the zero dynamics of the system needs to be assured. For that purpose, the system of (11) is transformed into the following state-space form

Where

= U²β is an auxiliary comtrol input,

is a new vector of system states, while

are explicitly defined as follows

where θ_i ∀ _i = 1,2,3,4 are constants that are explicitly defined in the Appendix. In (13) above,

denotes a nonlinearity that encodes the nonlinear structural stiffness. It is to be noted here Φ(0) = 0. The state-space system of (13) can be expanded into the following from

Here, the stability of the zero dynamics is studied for the case when the pitch displacement is

regulated to the origin. Mathematically, this implies that

which implies from the second equation of (15) that

Since Φ₁(0) = 0. The zero dynamics of the system then reduce to the reduce to the third order system given by

Substituting (16) into the above set of equations for

we obtain the linear system of equations

and A^*_c is given by

For the nominal system of (15), the eigenvalues of A^*_c lie in the left half plane which implies that the zero dynamics of the system are asymptotically stable, i.e., this is a minimum phase system. This implies that asymptotic convergence of the pitching variable α assures the asymptotic convergence of the plunging variable z.

Given the definitions of (13) and (14),

can be expressed as follows

The tracking error

is defined where

denotes the desired output vector which needs to be smooth in deference to the requirements of the subsequent control design. For the control objective, one can simply choose α_d to be zero all the time or use another desirable smooth time-varying trajectory α_d(t) along which the actual pitching variable α can be driven towards the origin. In order to facilitate the ensuing control design and stability analysis, we also define the tracking error

and the filtered tracking error signal

as follows

where λ₁, λ₂ are positive constants. By utilizing the definitions above, one can obtain

By substituting (17) for

in the above expression, the open-loop dynamics for r can be obtained as follows

After a convenient rearrangement of terms, the open-loop dynamics can be rewritten as follows

In order to design a model-free controller, we define an auxiliary nonlinear signal N(·) as follows

By utilizing the definition of (22) above, the open-loop dynamics of the system can be compactly

rewritten as follows

Since the structure of the model is assumed to be unknown in the control design, standard adaptive control cannot be applied. In its lieu, a neural network feedforward compensator

along with a robustifying term is proposed to compensate for the function N as defined above in (22). By the universal function approximation property [32], the nonlinear function of the system N can be approximated as a three-layer network target function as follows

as long N is a general smooth function from

to

and the set of inputs to the function is restricted to a compact set S of

. In (24),

denotes the augmented input vector, vector

is the ideal first layer interconnection weight matrix between input layer and hidden layer,

denotes the sigmoidal activation function, while

denotes the ideal second layer interconnection weight matrix. In this work, the weight matrixes W and V are assumed th be constant and bounded as ∥W∥_F ≤ W_B and ∥V∥_F ≤V_B, where W_B and V_B are positive constants. The approximation error is assumed to be bounded in compact set ∥ε∥ < ε_N where ε_N is an unknown positive constant related to the number of nodes in the hidden layer.

After substituting the approximation from (24) into (23), one can rewrite the open-loop dynamics as follows

where

Motivated by the open-loop dynamics and the ensuing stability analysis, the control law is designed as follows

where K_v, K_d > 0 are constant control gains,

is typical three-level neural network compensator for target function

defined as follows

v is a robustifying term which will be defined later while ？ is an adaptive estimate for g. The dynamic update law for ？ is designed as follows

where the parameter projection operator proj{·} is designed to bound ？ in a known compact set Ω such that sgn(g₃)？(t) ≥ ε > 0 for all time. The projection operator defined here is meaningul because the minimum-phase nature of the system ensures that sgn(g₃)g(t)=g₃^?1g₁ is always positive. In (25), ？ and

are estimates for the neural network interconnection weight matrices that are dynamically generated as follows

where

and

are postive definite diagonal gain matrixes, while k > 0 is a scalar design parameter. By substituting the expression for control law in (26) into the open-loop dynamics of (25) and conveniently rearranging the terms, one can obtalin the closed-loop system dynamics as follows

where

is a parameter estimation error. Also note that we can write

where the weight estimation errors are defined as

while w is defined as follows

To facilitate the subsequent analysis, one can also obtalin a compact form representation for ∥w∥ follows

where C₀, C₁ and C₂ are all positive constants while the ideal composite weight matrix Z, estimated composite weight matrix

and the composite weight mismatch matrix

are given as follows

Per the boundedness property for ∥W∥ _F and ∥ W∥ _F as described above, there exists a constant Z_B such that Z_B > ∥ Z∥ _F. Based on the definition Z_B, the robustifying term v can be designed as

Where K_z is a positive constant. Finally, it is noted that the functional reconstruction error

is assumed to be bounded. Thus, the closed-loop dynamics can be finally written as

In this section, we provide the stability analysis for the proposed model-free controller. We begin by defining a nonnegative Lyapunov function candidate V₂ as follows

After differentiating V₂ along the closed-loop dynamics of r(t) as well as (28), one can obtain the following expression for

After applying the neural network weight update laws designed in (29), canceling out the matched terms and utilizing the definitions of (31), (35) can be upperbounded as

By substituting (30) and (32) into (36), it is possible to further upperbound

as

where the following relation has been used to derive

Based on the fact that

one can choose

K_z > C₂ such that (37) can be cast as

By defining C₃ = Z_B + C₁ / k and conveniently rearranging the terms, (40) yields

By choosing K_d > [C₀ ? kC²₃ / 4], one can obtain the following upperbound on

From (34) and (40), it is easy to see that r ∈ L₂ ∩ L_∞ while ？,

The boundedness of r implies that α,

are bounded by virtue of the definitions of (18) and (19). Since the system is minimum phase and relative degree one, the boundedness of the output guarantees that any first order stable filtering of the input will remain bounded. This implies that all system states remain bounded in closedloop operation which further implies that

stays bounded. Since (26) defines a stable filter acting on a bounded input, it is easy to see that β and

stay bounded; furthermore, the flap deflection control input β is continuous at all times. The boundedness of β implies in turn that

by virtue of the closed-loop dynamics of r. Thus, using previous assertions, one can utilize Barbalat's Lemma [33] to conclude that r → 0 as t → ∞ which further implies that

From the asymptotic stability of the zero dynamics, we can further guarantee that x₃, x₄ → 0 as t → ∞ . Thus, both the pitching and plunging variables show asymptotic convergence to the origin.

In this section, simulation results are presented for an aeroelastic system controlled by the proposed continuous robust controller. The nonlinear aerodynamic model is simulated using the dynamics of (1), (7) and (10). The nominal model parameters are list as follows

and the controller parameters are listed in Table 1.

The desired trajectory variables α_d,

are simply selected as zero. The inital conditions for pitching displacement α(t) and plunging displacement ξ(t) are chosen as α(0) = 5.729deg(about 0.1 radians) and ξ(0) = 0 m, while all other state variables are initialized to zero. The initial parameter estimate ？(0) is set to be -1.20, which is a 10% shift from its nominal value. The flap deflection β(t) is constrained to vary between ± 15 deg.

The effect of structural nonlinearities on LCO amplitude was analyzed before applying any control. As shown in [22], increase in structural stiffness factor denoted by B led to decrease in LCO amplitude provided the flutter speed

remains constant. Furthermore, we also explored the effect of the location of the elastic axis from the leading edge. It was shown in [22] that a decrease in x₀ leads to decrease in

LCO amplitude while the flutter speed increases. It was also shown that increasing the damping ratios ζ_h and ζ_α resulted in decrease of the amplitude of the LCO.

Fig 2 shows the dynamics of open-loop pitching displacement α and plunging displacement ξ at pre-flutter speed. The simulation is carried out in subcritical flight speed regime, M = 2, below the flutter speed of M_flutter = 2.15.

Without the controller, it is obvious that the oscillation of pitching degree-of-freedom α will converge within 3[s] while the plunging displacement is lightly damped and it takes over 3[s] to converge. In Fig 3, it is shown that the proposed robust controller suppresses the oscillation of α in less than 1.5[s] while the plunging displacement ξ is suppressed in 2.5[s]. The parameter estimate g is seen to converge to a constant value within less than 0.5[s].

Another set of simulations is run for post-flutter speed. As shown in Fig 4, when M is set to be 3, the system dynamics show sustained limit cycle oscillations in open-loop operation. Such LCOs is experienced due to the non-linear pitch stiffness and the aerodynamic nonlinearities. After applying the control to the plant, from Fig 5, it is shown that when the control is turned on at t=0[s], the oscillation of α is suppressed within 1.5[s]. The dynamic oscillatory behavior of the plunging displacement ξ is suppressed within 2.5[s].The control performance is very satisfactory when it start to work at t=0. Next simulation is for delay open of control. In Fig 6, control was turned on at t=4[s] after the system had gone into an LCO. It is seen that the oscillations of the pitching degree α and plunging displacement ξ are suppressed respectively in 1.5[s] and 3[s], while in [13] convergence time of the two states are 1.5[s] and 4[s]. The parameter estimate ？ also converges to a constant in less than 0.5[s].

These simulation results show that the proposed novel robust controller can effectively suppress the oscillation of both pitching and plunging degrees-of-freedom of the airfoil in both pre-flutter and post-flutter flight speed regimes.

A modular model-free continuous robust controller was proposed to suppress the aeroelastic vibration characteristics (including flutter and limit cycle oscillations in pre- and post-flutter

condition) of a supersonic 2-DOF lifting surface with flap. Differently from traditional adaptive control strategies, which strictly require the linear parameterization of the system, no prior knowledge of the system model is required for the method presented in this paper. A Lyapunov method based analysis was provided to obtain the global asymptotic stability result. Finally, the simulation results showed that this control strategy can rapidly suppress any aeroelastic vibration

in pre- and post-flutter flight speed regimes.

The auxiliary constants as well as and that were introduced in the model description of the state-space model are explicitly defined follows

The elements introduced in (15) are explicitly defined as follows