Incremental Nonlinear Control for Aeroelastic Wing Load Alleviation and Flutter Suppression

Abstract

This paper proposes an incremental nonlinear control method for an aeroelastic system’s gust load alleviation and active flutter suppression. These two control objectives can be achieved without modifying the control architecture or the control parameters. The proposed method has guaranteed stability in the Lyapunov sense and also has robustness against external disturbances and model mismatches. The effectiveness of this control method is validated by wind tunnel tests of an active aeroelastic parametric wing apparatus, which is a typical wing section containing heave, pitch, flap, and spoiler degrees of freedom. Wind tunnel experiment results show that the proposed nonlinear incremental control can reduce the maximum gust loads by up to 46.7% and the root mean square of gust loads by up to 72.9%, while expanding the flutter margin by up to 15.9%.

1. Introduction

Modern transport aircraft commonly feature high-aspect-ratio wings to increase aerodynamic efficiency. A disadvantage of increasing the wing aspect ratio is the increased susceptibility to gust and manoeuvre loads in addition to the onset of aeroelastic phenomena such as flutter and divergence. In the literature, there are typically two approaches to alleviate gust loads and suppress flutter: passive aeroelastic tailoring and active control [1]. The passive approach has a long history and is achieved by exploiting the anisotropic properties of composite material to steer dynamic and static aeroelastic behaviour. By contrast, the active approach designs feedforward and/or feedback controllers to actuate the leading-edge and/or trailing-edge control surfaces for load redistribution and closed-loop dynamic modification. This paper focuses on the active approach because typically, it has a better adaptability to variations in flight and load conditions than its passive counterpart [1].

The majority of active gust load alleviation (GLA) and flutter suppression control algorithms are designed based on a reduced-order, linear, time-invariant state-space aero(servo)elastic model. These algorithms include proportional–integral–derivative control [2], pole placement [3], linear quadratic regulator/Gaussian [4], eigensystem synthesis, $\mu$ analysis [5], and linear robust control (${\mathcal{H}}_2$ and ${\mathcal{H}}_{\infty }$ [6]). Although these linear control approaches have shown their effectiveness in practice, the resulting controllers only have a guaranteed stability and performance around the linearization point; thus, the additional and tedious gain-scheduling method [7] is required to expand these linear controllers to a wider flight envelope. Furthermore, it is challenging for linear controllers to passively tolerate some specific nonlinearities (i.e., free-play, backlash hysteresis [8], bifurcation [9]), sudden faults in actuators and/or sensors, and structure damage.

On the contrary, nonlinear control methods, especially those that have guaranteed stability in the Lyapunov sense, have shown great potential in solving nonlinear aeroservoelastic system control problems without requiring the gain-scheduling technique. An immersion and invariance controller was proposed in [10] for nonlinear flutter suppression and free-play compensation. Platanitis et al. [11] used the nonlinear dynamic inversion approach together with the model reference adaptive control technique for the limit cycle oscillation suppression of a typical aeroelastic wing section. Recurrent neural networks have been used for nonlinear model identification and active flutter suppression in [12]. However, these nonlinear control approaches have a relatively high model dependency, while offline and/or online model identification of a nonlinear aeroservoelastic system is a nontrivial task. When sudden faults occur during flight, the convergence rate of online model identification can be insufficient to guarantee stability.

Different from model-based nonlinear control methods in the literature, the incremental nonlinear dynamic inversion (INDI) control is a sensor-based approach [13]. It is derived from nonlinear dynamic inversion (NDI) or feedback linearization, which linearizes the input–output mapping of a nonlinear system via feedback, resulting in a chain of integrators that can be easily stabilized by a linear virtual control [13]. The INDI method inherits the merits of NDI. More importantly, it greatly reduces the model dependency of NDI via exploiting sensor measurement [13]. In conventional NDI, since a perfect model is never known and external disturbances always exist in reality, the ideal linearization never exists, leading to robustness issues. By contrast, INDI makes full use of sensing information and simultaneously reduces the model dependency and improves control robustness. In the literature, INDI has been applied to a free-flying flexible aircraft tracking problem [14] and a morphing-wing gust load alleviation problem [15]. However, the effectiveness of INDI on flutter suppression, especially in a real-world environment, remains to be proven.

In this paper, we propose to use the sensor-based incremental nonlinear dynamic inversion control to tackle the active gust load alleviation and flutter suppression problems of an aeroservoelastic system. The goal is to use one single controller to simultaneously achieve these two objectives without requiring a gain adjustment, control architecture variations, or gain scheduling. The reduced model-dependency of INDI also reduces its practical implementation effort. The performance, robustness, and implementation ease of the proposed control method are validated by wind tunnel experiments on our newly developed active, aeroelastic, parametric wing apparatus.

The rest of this paper is structured as follows. The main methods are detailed in Section 2. Section 3 explains the experimental setup. The results are presented and discussed in Section 4, followed by the conclusion in Section 5. Finally, an outlook on future research is given in Section 6.

2. Methods

2.1. Dynamic Model for a Typical Aeroelastic Wing Section

The equations of motion for a typical aeroelastic wing section can be written as [16]

$$\left[\begin{array}{ccc}m& S& S_{\beta }\\ S& I_{\theta }& I_{\theta \beta }\\ S_{\beta }& I_{\theta \beta }& I_{\beta }\end{array}\right]\left[\begin{array}{c}\ddot{h}\\ \ddot{\theta }\\ \ddot{\beta }\end{array}\right]+\left[\begin{array}{ccc}{K}_h& 0& 0\\ 0& K_{\theta }& 0\\ 0& 0& K_{\beta }\end{array}\right]\left[\begin{array}{c}h\\ \theta \\ \beta \end{array}\right]=\left[\begin{array}{c}-L_h\\ M_{\theta }\\ M_{\beta }\end{array}\right]$$

(1)

where h is the vertical displacement or plunge of the airfoil; $\theta$ is the pitch angle of the wing section; $\beta$ is the deflection angle of the control surface; m is the mass per unit length of the wing section; S is the static mass moment of the wing around the pitch axis; $I_{\theta }$ is its mass moment of inertial around the same axis; the static mass moment of the control surface around the hinge axis is denoted by $S_{\beta }$; $I_{\beta }$ is the moment of inertia of the control surface around the hinge axis; $I_{\theta \beta }$ is the product of inertia; $K_h,K_{\theta },K_{\beta }$ denote the extension spring’s stiffness, the torsional spring’s stiffness related to $\theta$, and the torsional spring’s stiffness related to $\beta$, respectively; $L_h$ is the lift; $M_{\theta }$ is the pitching moment of the wing section around the pitch axis; and $M_{\beta }$ is the pitching moment of the control surface around the hinge axis. We assume structural damping is negligible, which is a common assumption in the field. An illustration of a typical aeroelastic wing section is show in Figure 1.

Theodorsen’s unsteady aerodynamic theory was adopted [16]. We indicate the half-chord length as b. Using the exponential approximation of Wagner’s function $\mathsf{\Phi }\left(t\right)=1-{\varphi }_1{e}^{-{ϵ}_1\frac{Ut}{b}}-{\varphi }_2{e}^{-{ϵ}_2\frac{Ut}{b}}$ with ${\varphi }_1=0.165$, ${\varphi }_2=0.335$, ${ϵ}_1=0.0455$, and ${ϵ}_2=0.3$ [16], we obtain the complete equations of motion for a pitch–plunge airfoil with control surface as

$$\begin{array}{cc}\hfill & (\mathit{A}+\rho \mathit{B})\ddot{\mathit{y}}+(\mathit{C}+\rho U\mathit{D})\dot{\mathit{y}}+(\mathit{E}+\rho U^2\mathit{F})\mathit{y}+\rho U^3\mathit{W}\mathit{w}={\left[\begin{array}{ccc}-L_h& M_{\theta }& M_{\beta }\end{array}\right]}_{\mathrm{external}}^{\mathsf{T}}\hfill \\ & \dot{\mathit{w}}={\mathit{W}}_1\mathit{y}+U{\mathit{W}}_2\mathit{w}\hfill \end{array}$$

(2)

where $\mathit{y}={[h,\theta ,\beta ]}^{\mathsf{T}}$ and $\mathit{w}={[w_1,w_2,\cdots ,w_6]}^{\mathsf{T}}$ represents the aerodynamic lag states, $\rho$ is the density of the flow, U is the speed of the flow. The detailed expressions of the matrices can be derived from Equation (1) and can also be found in [16]. The external lift $L_h$ and moments $M_{\theta },M_{\beta }$ in Equation (2) are gusts- and control-surface-deflection-induced components of the total lift and moments in Equation (1).

The loads due to a gust are represented as ${[L_h^g,M_{\theta }^g,M_{\beta }^g]}^{\mathsf{T}}$, in which the lift is a purely circulatory event. The gust velocity $u_g$ is implemented in circulatory lift calculations with the Küssner function. The Küssner function is approximated with exponential form $\Psi \left(t\right)=1-{\psi }_1{e}^{-{\epsilon }_{k_1}\frac{Ut}{b}}-{\psi }_2{e}^{-{\epsilon }_{k_2}\frac{Ut}{b}}$ with ${\psi }_1=0.5$, ${\psi }_2=0.5$, ${\epsilon }_{k_1}=0.13$, and ${\epsilon }_{k_2}=1$. Denote the gust input velocity as $u_g$; denote $w_{g1}$ and $w_{g2}$ as the aerodynamic lag states due to gust; then, we have

$$\begin{array}{ccc}\hfill w_{g1}& ={\int }_0^t{e}^{-{\epsilon }_{k_1}(t-\tau )\frac{U}{b}}{u}_g\left(\tau \right)\mathrm{d}\tau ,w_{g2}\hfill & \hfill ={\int }_0^t{e}^{-{\epsilon }_{k_2}(t-\tau )\frac{U}{b}}{u}_g\left(\tau \right)\mathrm{d}\tau \end{array}$$

(3)

After applying the Leibnitz integration rule and using the Küssner function, the complete expressions for the aerodynamic loads due to a gust are defined as:

$$\begin{array}{c}{\dot{\mathit{w}}}_g=\left[\begin{array}{c}{\dot{w}}_{g1}\\ {\dot{w}}_{g2}\end{array}\right]=\left[\begin{array}{cc}-\frac{U{\epsilon }_{k_1}}{b}& 0\\ 0& -\frac{U{\epsilon }_{k_2}}{b}\end{array}\right]\left[\begin{array}{c}{w}_{g1}\\ w_{g2}\end{array}\right]+\left[\begin{array}{c}1\\ 1\end{array}\right]u_g\triangleq {\mathit{W}}_{g1}{u}_g+U{\mathit{W}}_{g2}{\mathit{w}}_g\\ {\mathit{y}}_g=\left[\begin{array}{c}-L_h^g\\ M_{\theta }^g\\ M_{\beta }^{\theta }\end{array}\right]=2\pi \rho U^2\left[\begin{array}{cc}-{\varphi }_1{\epsilon }_{k_1}& -{\varphi }_2{\epsilon }_{k_2}\\ -(a+0.5)b{\varphi }_1{\epsilon }_{k_1}& -(a+0.5)b{\varphi }_2{\epsilon }_{k_2}\\ \frac{b{T}_{12}}{2\pi }{\varphi }_1{\epsilon }_{k_1}& \frac{b{T}_{12}}{2\pi }{\varphi }_2{\epsilon }_{k_2}\end{array}\right]{\mathit{w}}_g\triangleq {\mathit{C}}_g{\mathit{w}}_g\end{array}$$

(4)

where the parameter a and $T_{12}$ are defined as in [16].

Assume that the inertial coupling force/moment for the control surface with regard to the rest of the wing section are negligible; then, the independent servo actuator dynamics is governed by

$$\ddot{\beta }=-2{\zeta }_s{\omega }_s\dot{\beta }-{\omega }_s^2\beta +{\omega }_s^2{\beta }_c$$

(5)

where ${\beta }_c$ is the control command for the aileron, and the coefficients ${\zeta }_s,{\omega }_s$ are obtained using the system identification of the servo.

Then, the dynamics for state $\overline{\mathit{y}}={[h,\theta ]}^T$ is governed by

$$\begin{array}{cc}\hfill & (\overline{\mathit{A}}+\rho \overline{\mathit{B}})\ddot{\overline{\mathit{y}}}+(\overline{\mathit{C}}+\rho U\overline{\mathit{D}})\dot{\overline{\mathit{y}}}+(\overline{\mathit{E}}+\rho U^2\overline{\mathit{F}})\overline{\mathit{y}}+\rho U^3\overline{\mathit{W}}\mathit{w}={\left[\begin{array}{cc}-L_h& M_{\theta }\end{array}\right]}_{\mathrm{external}}^{\mathsf{T}}\hfill \\ & {\dot{\mathit{w}}}_{1-4}={\mathit{W}}_{1,1-4}\overline{\mathit{y}}+U{\mathit{W}}_{2,1-4}{\mathit{w}}_{1-4}\hfill \\ & {\dot{\mathit{w}}}_{5-6}={\mathit{W}}_{1,5-6}\beta +U{\mathit{W}}_{2,5-6}{\mathit{w}}_{5-6}\hfill \end{array}$$

(6)

where $\overline{\mathit{A}},\overline{\mathit{B}},\overline{\mathit{C}},\overline{\mathit{D}},\overline{\mathit{E}},\overline{\mathit{F}}\in {\mathbb{R}}^{2\times 2}$ and $\overline{\mathit{W}}\in {\mathbb{R}}^{2\times 6}$ are the corresponding matrices for $h,\theta$ fragmented from Equation (2).

The resulting external force and moment due to the control surface deflection $\beta$ are modelled as

$$\left[\begin{array}{c}{L}_h^{\beta }\\ M_{\theta }^{\beta }\end{array}\right]=\left[\begin{array}{ccc}\rho U^2{\mathit{F}}_{1,3}& \rho U{\mathit{D}}_{1,3}& \rho {\mathit{B}}_{1,3}\\ -\rho U^2{\mathit{F}}_{2,3}& -\rho U{\mathit{D}}_{2,3}& -\rho {\mathit{B}}_{2,3}\end{array}\right]\left[\begin{array}{c}\beta \\ \dot{\beta }\\ \ddot{\beta }\end{array}\right]$$

(7)

where

$$\begin{array}{cc}\hfill & {\mathit{F}}_{1,3}=2b{T}_{10}\mathsf{\Phi }\left(0\right)+\Xi b^2{T}_{11},{\mathit{F}}_{2,3}=b^2(T_4+T_{10})-2b^2(a+0.5)T_{10}\mathsf{\Phi }\left(0\right)-\Xi b^3(a+0.5)T_{11}\hfill \\ \hfill & {\mathit{D}}_{1,3}=-b^2{T}_4+\mathsf{\Phi }\left(0\right)b^2{T}_{11},{\mathit{D}}_{2,3}=b^3(T_1-T_8-(c_h-a)T_4+0.5T_{11})-\mathsf{\Phi }\left(0\right)b^3(a+0.5)T_{11}\hfill \\ \hfill & {\mathit{B}}_{1,3}=-b^3{T}_1,{\mathit{B}}_{2,3}=-b^4(T_7+(c_h-a)T_1),\mathsf{\Phi }\left(0\right)=1-{\psi }_1-{\psi }_2,\Xi =\frac{{\psi }_1{ϵ}_1}{b}+\frac{{\psi }_2{ϵ}_2}{b}\hfill \end{array}$$

and $c_h$ and $T_1$ to $T_{14}$ are defined as in [16].

Ignoring the inertia coupling between the aileron actuator and the rest of the wing section and substituting Equation (5) into Equation (7), we have

$$\left[\begin{array}{c}-L_h^{\beta }\\ M_{\theta }^{\beta }\end{array}\right]=\left[\begin{array}{cc}-\rho U^2{\mathit{F}}_{1,3}+\rho {\mathit{B}}_{1,3}{\omega }_s^2& -\rho U{\mathit{D}}_{1,3}+2\rho {\mathit{B}}_{1,3}{\zeta }_s{\omega }_s\\ -\rho U^2{\mathit{F}}_{2,3}+\rho {\mathit{B}}_{2,3}{\omega }_s^2& -\rho U{\mathit{D}}_{2,3}+2\rho {\mathit{B}}_{2,3}{\zeta }_s{\omega }_s\end{array}\right]\left[\begin{array}{c}\beta \\ \dot{\beta }\end{array}\right]+\left[\begin{array}{c}-\rho {\mathit{B}}_{1,3}{\omega }_s^2\\ -\rho {\mathit{B}}_{2,3}{\omega }_s^2\end{array}\right]{\beta }_c& \triangleq {\mathit{C}}_{\beta }\left[\begin{array}{c}\beta \\ \dot{\beta }\end{array}\right]+{\mathit{D}}_{\beta }{\beta }_c$$

(8)

Now, choose state $\mathit{x}={[h,\theta ,\dot{h},\dot{\theta },{\mathit{w}}^{\mathsf{T}},\beta ,\dot{\beta },{\mathit{w}}_g^{\mathsf{T}}]}^{\mathsf{T}}$, control output $y_c=h$, and control input $u={\beta }_c$. The following control-oriented state-space model is obtained:

$$\dot{\mathit{x}}\left(t\right)={\mathcal{A}}_c\left(t\right)\mathit{x}\left(t\right)+{\mathcal{B}}_c\left(t\right)u\left(t\right)+{\mathcal{D}}_c{u}_g\left(t\right)$$

(9)

where

$$\begin{array}{cc}\hfill & {\mathcal{A}}_c=\left[\begin{array}{cccccc}\mathbf{0}& {\mathit{I}}_2& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ -{\overline{\mathit{M}}}_{ae}^{-1}{\overline{\mathit{K}}}_{ae}& -{\overline{\mathit{M}}}_{ae}^{-1}{\overline{\mathit{C}}}_{ae}& -{\overline{\mathit{M}}}_{ae}^{-1}\rho U^3{\overline{\mathit{W}}}_{1-4}& -{\overline{\mathit{M}}}_{ae}^{-1}\rho U^3{\overline{\mathit{W}}}_{5-6}& {\overline{\mathit{M}}}_{ae}^{-1}{\mathit{C}}_{\beta }& {\overline{\mathit{M}}}_{ae}^{-1}{\overline{\mathit{C}}}_g\\ {\mathit{I}}_4& \mathbf{0}& U{\mathit{W}}_{2,1-4}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& U{\mathit{W}}_{2,5-6}& {\mathit{I}}_2& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{A}}_{\beta }& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& U{\mathit{W}}_{g2}\end{array}\right]\hfill \\ \hfill & {\mathcal{B}}_c={\left[\begin{array}{ccccccccc}0& 0& {\overline{\mathit{M}}}_{ae}^{-1}{\mathit{D}}_{\beta }& 0& 0& 0& {\omega }_s^2& 0& 0\end{array}\right]}^{\mathsf{T}}\hfill \\ & {\mathcal{D}}_c={\left[\begin{array}{ccc}{\mathbf{0}}_{1\times 10}& 1& 1\end{array}\right]}^{\mathsf{T}}\hfill \end{array}$$

(10)

with

$$\begin{array}{cc}\hfill & {\overline{\mathit{M}}}_{ae}=\overline{\mathit{A}}+\rho \overline{\mathit{B}},{\overline{\mathit{K}}}_{ae}=\overline{\mathit{E}}+\rho U^2\overline{\mathit{F}},{\overline{\mathit{C}}}_{ae}=\overline{\mathit{C}}+\rho U\overline{\mathit{D}}\hfill \end{array}$$

(11)

$${\mathit{A}}_{\beta }=\left[\begin{array}{cc}0& 1\\ -{\omega }_s^2& -2{\zeta }_s{\omega }_s\end{array}\right],{\overline{\mathit{C}}}_g=2\pi \rho U^2\left[\begin{array}{cc}-{\varphi }_1{\epsilon }_{k_1}& -{\varphi }_2{\epsilon }_{k_2}\\ -(a+0.5)b{\varphi }_1{\epsilon }_{k_1}& -(a+0.5)b{\varphi }_2{\epsilon }_{k_2}\end{array}\right]$$

(12)

2.2. Incremental Nonlinear Dynamic Inversion

The incremental nonlinear dynamic inversion (INDI) method can control the following nonlinear system:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\right)+\mathit{G}\left(\mathit{x}\right)\mathit{u}+\mathit{d}\left(t\right),\mathit{y}=\mathit{h}\left(\mathit{x}\right)$$

(13)

where $\mathit{f}:{\mathbb{R}}^n\to {\mathbb{R}}^n$ and $\mathit{h}:{\mathbb{R}}^n\to {\mathbb{R}}^p$ are smooth vector fields. $\mathit{G}$ is a smooth function mapping ${\mathbb{R}}^n\to {\mathbb{R}}^{n\times m}$, whose columns are smooth vector fields. The external disturbance vector is $\mathit{d}\left(t\right)\in {\mathbb{R}}^n$, which is assumed to satisfy ${\parallel \mathit{d}\left(t\right)\parallel }_2\le \overline{d}$. The external disturbances in the real world can easily satisfy this boundedness assumption.

In Equation (13), $\mathit{y}\in {\mathbb{R}}^p$ denotes the controlled output vector. This paper considers the case where $p\le m$, which means the system is either fully actuated or overactuated. The vector relative degree [17] of the system is defined as $\mathit{\rho }={[{\rho }_1,{\rho }_2,...,{\rho }_p]}^{\mathsf{T}}$, which satisfies $\rho ={\parallel \mathit{\rho }\parallel }_1={\sum }_{i=1}^p{\rho }_i\le n$. From Equation (13), the input–output mapping of the nonlinear system is

$${\mathit{y}}^{\left(\mathit{\rho }\right)}=\mathit{\alpha }\left(\mathit{x}\right)+\mathcal{B}\left(\mathit{x}\right)\mathit{u}+{\mathit{d}}_y$$

(14)

In Equation (14), $\mathit{\alpha }\left(\mathit{x}\right)={[{\mathcal{L}}_f^{{\rho }_1}{h}_1,{\mathcal{L}}_f^{{\rho }_2}{h}_2,...,{\mathcal{L}}_f^{{\rho }_p}{h}_p]}^{\mathsf{T}}$, $\mathcal{B}\left(\mathit{x}\right)\in {\mathbb{R}}^{p\times m}$, and ${\mathcal{B}}_{ij}={\mathcal{L}}_{g_j}{\mathcal{L}}_f^{{\rho }_i-1}{h}_i$, where ${\mathcal{L}}_f^{{\rho }_i}{h}_i,{\mathcal{L}}_{g_j}{\mathcal{L}}_f^{{\rho }_i-1}{h}_i$ are the corresponding Lie derivatives [18]. When ${\rho }_i=1$ for all $i=1,...,p$, ${\mathit{d}}_y={[{\mathcal{L}}_d{h}_1,{\mathcal{L}}_d{h}_2,...,{\mathcal{L}}_d{h}_p]}^{\mathsf{T}}$. When $\rho =n$, the system given by Equation (13) is full-state feedback linearizable. Otherwise, $n-\rho$ internal dynamics exist.

INDI considers system variations in one sampling interval $\Delta t$. The incremental dynamic equation is derived by taking the first-order Taylor series expansion of Equation (14) around the condition at $t-\Delta t$ (denoted by the subscript 0) as:

$${\mathit{y}}^{\left(\mathit{\rho }\right)}={\mathit{y}}_0^{\left(\mathit{\rho }\right)}+\frac{\partial \left[\mathit{\alpha }\right(\mathit{x})+\mathcal{B}(\mathit{x}\left)\mathit{u}\right]}{\partial \mathit{x}}{|}_0\Delta \mathit{x}+\mathcal{B}\left({\mathit{x}}_0\right)\Delta \mathit{u}+\Delta {\mathit{d}}_y+{\mathit{R}}_1$$

(15)

in which $\Delta \mathit{x}$, $\Delta \mathit{u}$, and $\Delta {\mathit{d}}_y$ represent the state, control, and disturbance increments in one sampling time step $\Delta t$, respectively. ${\mathit{R}}_1$ is the expansion remainder. Define the internal state vector as $\mathit{\eta }=\mathit{\varphi }\left(\mathit{x}\right)$ and the external state vector as $\mathit{\xi }={[{\mathit{\xi }}_1^{\mathsf{T}},...,{\mathit{\xi }}_p^{\mathsf{T}}]}^{\mathsf{T}}$, ${\mathit{\xi }}_i={[h_i\left(\mathit{x}\right),...,{\mathcal{L}}_f^{{\rho }_i-1}{h}_i\left(\mathit{x}\right)]}^{\mathsf{T}}$, $i=1,...,p$. In a stabilization problem, the reference for the controlled output equals zero and the control increment is designed to satisfy

$$\overline{\mathcal{B}}\left({\mathit{x}}_0\right)\Delta {\mathit{u}}_{\mathrm{indi}}={\mathit{\nu }}_c-{\mathit{y}}_0^{\left(\mathit{\rho }\right)},{\mathit{\nu }}_c=-\mathit{K}\mathit{\xi }$$

(16)

where $\overline{\mathcal{B}}$ is an estimation of $\mathcal{B}$. The gain matrix $\mathit{K}=\mathrm{diag}\left\{{\mathit{K}}_i\right\},i=1,2,...,p$, and ${\mathit{K}}_i=[K_{i,0},...,K_{i,{\rho }_i-1}]$. ${\mathit{y}}_0^{\left(\mathit{\rho }\right)}$ can be directly measured or estimated. The total control command for the actuator is ${\mathit{u}}_{\mathrm{indi}}={\mathit{u}}_{\mathrm{indi},0}+\Delta {\mathit{u}}_{\mathrm{indi}}$, in which ${\mathit{u}}_{\mathrm{indi},0}$ can be measured or estimated using an actuator model. Considering the internal dynamics, the resulting closed-loop dynamics are:

$$\begin{array}{ccc}\hfill \dot{\mathit{\eta }}& =& {\mathit{f}}_{\eta }(\mathit{\eta },\mathit{\xi },\mathit{d})=\frac{\partial \mathit{\varphi }}{\partial \mathit{x}}(\mathit{f}\left(\mathit{x}\right)+\mathit{d}\left(t\right)){|}_{\mathit{x}={\mathit{T}}^{-1}\left(\mathit{z}\right)}\hfill \\ \hfill \dot{\mathit{\xi }}& =& ({\mathit{A}}_c-{\mathit{B}}_c\mathit{K})\mathit{\xi }+{\mathit{B}}_c[\delta (\mathit{x},\Delta t)+(\mathcal{B}\left({\mathit{x}}_0\right)-\overline{\mathcal{B}}\left({\mathit{x}}_0\right))\Delta {\mathit{u}}_{\mathrm{indi}}+\Delta {\mathit{d}}_y]\hfill \\ & \triangleq & ({\mathit{A}}_c-{\mathit{B}}_c\mathit{K})\mathit{\xi }+{\mathit{B}}_c{\mathit{\epsilon }}_{\mathrm{indi}}\hfill \end{array}$$

(17)

where $\mathit{z}=\mathit{T}\left(\mathit{x}\right)={[{\mathit{\eta }}^{\mathsf{T}},{\mathit{\xi }}^{\mathsf{T}}]}^{\mathsf{T}}$ is a diffeomorphism. The term $\delta (\mathit{x},\Delta t)$ is the closed-loop value of the variations and expansion reminder:

$$\mathit{\delta }(\mathit{x},\Delta t)=\left[\frac{\partial \left[\mathit{\alpha }\right(\mathit{x})+\mathcal{B}(\mathit{x}\left)\mathit{u}\right]}{\partial \mathit{x}}{|}_0\Delta \mathit{x}+{\mathit{R}}_1\right]{|}_{\mathit{u}={\mathit{u}}_{\mathrm{indi}}}$$

(18)

${\mathit{A}}_c=\mathrm{diag}\left\{{\mathit{A}}_0^i\right\}$, ${\mathit{B}}_c=\mathrm{diag}\left\{{\mathit{B}}_0^i\right\}$, ${\mathit{C}}_c=\mathrm{diag}\left\{{\mathit{C}}_0^i\right\}$, $i=1,2,...,p$, and $({\mathit{A}}_0^i,{\mathit{B}}_0^i,{\mathit{C}}_0^i)$ is the canonical form representation of a chain of ${\rho }_i$ integrators. The gain matrix $\mathit{K}$ is designed to guarantee that ${\mathit{A}}_c-{\mathit{B}}_c\mathit{K}$ is Hurwitz [13].

Theorem 1.

If $\parallel {\mathit{\epsilon }}_{\mathrm{indi}}{\parallel }_2\le \overline{\epsilon }$ is satisfied for all $\mathit{\xi }\in {\mathbb{R}}^{\rho }$, ${\mathit{f}}_{\eta }(\mathit{\eta },\mathit{\xi },\mathit{d})$ is continuously differentiable and globally Lipschitz in $(\mathit{\eta },\mathit{\xi },\mathit{d})$, and the origin of $\dot{\mathit{\eta }}={\mathit{f}}_{\eta }(\mathit{\eta },\mathbf{0},\mathbf{0})$ is globally exponentially stable, then the external state ξ in Equation (17) is globally ultimately bounded by a class $\mathcal{K}$ function of $\overline{\epsilon }$, while the internal state η in Equation (17) is globally ultimately bounded by a class $\mathcal{K}$ function of $\overline{\epsilon }$ and $\overline{d}$.

Proof .

This can be proved by applying Theorem 1 in Ref. [15] and setting the reference vector to zero. □

Theorem 2.

If $\parallel {\mathit{\epsilon }}_{\mathrm{indi}}{\parallel }_2\le \overline{\epsilon }$ is satisfied for all $\mathit{\xi }\in {\mathbb{R}}^{\rho }$, ${\mathit{f}}_{\eta }(\mathit{\eta },\mathit{\xi },\mathit{d})$ is continuously differentiable, and the origin of $\dot{\mathit{\eta }}={\mathit{f}}_{\eta }(\mathit{\eta },\mathbf{0},\mathbf{0})$ is exponentially stable, then there exists a neighborhood $D_z$ of $\mathit{z}={\mathbf{0}}^{\mathsf{T}}$ and ${\epsilon }^{*}>0$, such that for every initial state $\mathit{z}\left(0\right)\in D_z$ and $\overline{\epsilon }<{\epsilon }^{*}$, the external state ξ in Equation (17) is ultimately bounded by a class $\mathcal{K}$ function of $\overline{\epsilon }$, while the internal state η in Equation (17) is ultimately bounded by a class $\mathcal{K}$ function of $\overline{\epsilon }$ and $\overline{d}$.

Proof .

This can be proved by applying Theorem 2 in Ref. [15] and setting the reference vector to zero.

Theorems 1 and 2 prove that under the perturbation of bounded uncertainties and disturbances and with stable internal dynamics, the closed-loop system under INDI control is stable in the Lyapunov sense. □

2.3. INDI Design for an Aeroelastic Wing

In this research, we focused on stabilizing the heave degree of freedom of an aeroelastic wing using the trailing-edge flap. Consequently, in Equation (9), the controlled output is h, while the input $u={\beta }_c$. This leads to an input–output mapping with a relative degree equal to two. Applying Equation (16), the aileron’s input for stabilization is designed as:

$$\Delta u_{\mathrm{indi}}=({\nu }_c-{\ddot{h}}_0)/\overline{\mathcal{B}}\left({\mathit{x}}_0\right),{\nu }_c=-K_d\dot{h}-K_{p}h$$

(19)

where $K_d>0$ and $K_p>0$ are differential and proportional gains. The aileron’s control input command equals ${\beta }_c=u_{\mathrm{indi}}=u_{\mathrm{indi},0}+\Delta u_{\mathrm{indi}}$, in which $u_{\mathrm{indi},0}$ is the control input at the previous time point.

In Section 3, we validate the performance of the INDI method for gust load alleviation and active flutter suppression challenges. It is noteworthy that this single controller can solve both issues without changing the control architecture or the control gains. A control and block diagram is presented in Figure 2.

3. Experiment Setup

The proposed controller for active gust load alleviation and flutter suppression was validated using an aeroelastic wing apparatus by wind-tunnel testing. The self-designed active wing section [19] itself was mounted in an aeroelastic apparatus (AA) developed by Gjerek et al. [20]. The AA consisted of a rectangular, acrylic section mounted on the gust generator and providing heave and pitch degrees of freedom, with adjustable stiffness. In addition to the adjustable stiffness, weights can be added forward or rearward of the pitch axis, allowing both the mass distribution of the wing to be easily changed. The wing section was equipped with a movable full-span trailing-edge control surface (aileron) and a spoiler, which were actuated by high-bandwidth, electric servo actuators. The main focus of this experiment was to use the aileron for controlling the heave degree of freedom. The wing was also equipped with an MPU-9250 inertial measurement unit (IMU), as well as linear and rotation variable differential transformers (LVDT-Sentech 75DC-500/RVDT-Midori QP-2HC) for the measurements of acceleration, angular rates, heave, and pitch. All combined, this AA provided a setup that closely resembled the typical section, a two-dimensional wing section with an aileron and heave and pitch degrees of freedom (DOFs), in combination with the aerodynamic model developed by Theodorsen [21], which was presented in Section 2.1.

Wind tunnel testing was performed at the low-speed W-tunnel at Delft University of Technology, which is an open-circuit blow-down tunnel with a 0.4 m × 0.4 m test section, with low turbulence levels and a maximum attainable speed of 35 m/s. Attached to the wind tunnel is a gust generator capable of generating sinusoidal and $1-$ cosine gust excitations with gust frequencies ranging from 0.5 Hz to 12 Hz in 0.5 Hz increments [22]. The aeroelastic apparatus in the wind tunnel is shown in Figure 3.

The heave DOF was provided by two pairs of cantilever leaf springs, with one end of the springs clamped to the AA and the other end connected to a pitch assembly. The axles protruding from both sides of the wing were connected to bearings in the pitch assembly. The length of the leaf springs could be changed, providing a variable spring stiffness in the heave. Torsional stiffness was provided by a pair of axial springs connected to one of the axles by a pulley. The torsional stiffness could be varied by changing the diameter of the pulley or exchanging the axial springs. A top view of the AA is seen in Figure 4, showing the top half of the heave and pitch mechanism.

The specifications of the aeroelastic wing apparatus for GLA and flutter suppression are shown in

Table 1. Configuration parameters of the aeroelastic apparatus.

Parameter		Value	Unit
Heave stiffness	$K_h$	710	N/m
Pitch stiffness	$K_{\theta }$	3.14	Nm/rad
Elastic Axis	$EA$	$0.4c$	-

. When the blades of gust generator deflected at a certain frequency, gusts were generated in the wing section test field with a gust-induced angle of attack ${\alpha }_g=\frac{1}{2}{A}_{g0}{a}_m(1-cos\left(2\pi f_g(t-t_0)\right))$, where $t_0$ is the initial time, $a_m$ is the amplitude of the blades’ deflection angle, and $A_{g0}$ is the gust coefficient relevant to the velocity U and the gust generator frequency $f_g$. A second-order model for the aileron’s actuation mechanism was given by $G_{\mathrm{act}}\left(s\right)=\frac{2.6s+347.8}{s^2+34.7s+358.3}$. The minimal and maximal deflections were limited to $-{20}^{°}$ and ${20}^{°}$, and the maximum deflection rate was estimated to be 750${}^{\circ }$/s. The IMU used on the apparatus had a bandwidth of 200 Hz. The LVDT had a cut-off frequency of 200 Hz and was read by a 12-bit analogue-to-digital converter, resulting in a resolution of approximately 6 × 10${}^{-3}$ mm. The sampling interval was set as 0.002 s for capturing the high-frequency aeroelastic modes. The measured outputs used by the controller were the heave acceleration $\ddot{h}$ from the IMU, the heave displacement h from the LVDT, and the control surface deflection ${\beta }_0$.

4. Results and Discussion

In this section, the experimental results are discussed. First, the GLA results are presented in Section 4.1. Then, this is followed by the flutter suppression results in Section 4.2.

4.1. Gust Load Alleviation Results

The GLA results are presented in this subsection. In the wind tunnel test, $K_p=0.12$ and $K_d=0.009$ were chosen for the controller gains based on the desired eigenvalues of the closed-loop system. The nominal control effectiveness for 12 m/s was identified from a wind tunnel test as $\widehat{\overline{\mathcal{B}}}=4.0$. The open- and closed-loop experiments were performed for gust frequencies $f_g$ of 3, 3.5, 4, 4.5, and 5 Hz, and a gust amplitude of $a_m=15deg$. In addition to the directly measured outputs $\ddot{h}$ and h from the IMU and LVDT, a Luenberger observer with eigenvalues $[-150,-30]$ was applied to provide an estimation of $\widehat{\dot{h}}$ by using the measurements of the heave displacement and acceleration. $\widehat{\dot{h}}$ was required for the implementation of the control law as shown in Equation (19). An overview of the measured and estimated signals is shown in Figure 5.

For the GLA experiments, the wing was subjected to a series of $1-$ cosine gusts. First, the open-loop gust response was determined, after which the experiment was repeated with the controller enabled. Figure 6, Figure 7 and Figure 8 show the recorded GLA data for gust frequencies of 3, 4, and 5 Hz. The results show significant improvements in reducing the amplitude of h. The top plot of each subfigures shows the theoretical gust input in terms of the gust-induced angle of attack ${\alpha }_g$. At each frequency, the experiments were repeated four times to evaluate the coherence of the results.

The middle subplots show comparisons of the mean values of the open-loop (dashed line) and closed-loop responses (solid line). Also included in each middle subplot is the standard deviation of the heave response as a shaded region, which nearly coincides with the mean data, proving the coherence and repeatability of the experimental results. The transient, open-loop heave response in the middle plots can be subdivided into two parts. During the first half-period of the oscillation, the response is dominated by the gust input. For the remainder of the transient response, the wing oscillates with a frequency of approximately 4 Hz, close to the frequency of the first heave mode of the wing section. Whereas the open-loop results are highly underdamped, the closed-loop responses are mostly damped out after one full period.

The commanded and actual aileron deflections are detailed in the bottom subplots of Figure 6, Figure 7 and Figure 8. It can be observed that the aileron is settled on a steady-state deflection after the convergence of the heave response. Results from Schildkamp et al. [19] show the magnitude and phase responses of the servo actuator degrades for frequencies higher than 2 Hz. Since the servo actuator is now connected to the aileron, adding inertia and friction to the system, the frequency response is expected to be worse than what was presented in [19]. The commanded aileron deflection has a frequency similar to that of the heave response, explaining the differences in magnitude and phase between the commanded and actual aileron deflection. The larger variability seen in both aileron deflection signals is attributed to noise in the signals driving the controller. The excitation frequency of 5 Hz is close to the pitch mode natural frequency of 6.39 Hz. Since this study focused on the heave motion, the uncontrolled pitch motion inertially coupled with the control surface’s movement, leading to the high frequency and low magnitude oscillations in control deflections in Figure 8, although the steady-state oscillations in heave were damped out.

A summary of the GLA results in terms of the reduction of the absolute peak value and RMS heave values for all previously mentioned gust frequencies is shown in

Table 2. GLA controller reduction rate of the heave displacement h.

${\mathit{f}}_{\mathit{g}}$ (Hz)	Reduction Rate of $max\left(\right|\mathit{h}\left|\right)$ (%)			Reduction Rate of RMS $\left(\mathit{h}\right)$ (%)
	Mean	Min	Max	Mean	Min	Max
3	38.3	33.8	44.6	51.1	46.0	54.5
3.5	29.3	27.0	30.8	58.7	58.7	59.7
4	32.4	31.5	34.4	63.5	61.2	64.1
4.5	33.2	29.8	35.8	61.7	44.0	67.4
5	44.2	40.6	46.7	71.4	63.0	72.9

, where the absolute peak heave value relates to the peak load endured by the wing, and the RMS value gives a measure of the vibrational loads the wing endures, related to the fatigue life of a structure. The mean, minimum, and maximum reduction rates for both the absolute peak and RMS values are given. Overall, the proposed INDI control method, without adjusting control gains, provides an attenuation higher than $27\%$ for the vibration amplitude due to a gust disturbance, and an attenuation higher than $44\%$ for the RMS(h). The greatest reduction in absolute peak and RMS mean values are achieved for a gust frequency of 5 Hz, with values of 71.4 and 44.2%, respectively. Previous results show the overall lowest damping coefficient for the open-loop gust response at 5 Hz, giving the highest RMS value. The controller quickly damps out the gust-induced oscillations, leading to the second-lowest closed-loop RMS value, hence the greatest reduction in RMS. Similarly, this is also the case for the absolute peak value, with the highest open-loop peak value and the second lowest peak value.

In addition to the open- and closed-loop heave responses plotted in the time domain, Figure 9 shows the power spectral density (PSD) of these heave responses in the frequency domain for gust frequencies of 3, 4, and 5 Hz. Also indicated in this figure are the (1) first heave mode, (2) first pitch mode, and (3) first rocking mode at 3.55, 6.39, and 11.10 Hz, respectively, as identified using a ground vibration test (GVT) described in [19]. The open-loop response results in two distinct peaks, around 4 and 12 Hz, at frequencies slightly higher than the first heave and rocking modes. Changes to the hardware and wiring of the wing section were made after the GVT, likely affecting the identified frequencies. This will be verified during a future GVT.

Three observations can be made from the closed-loop PSDs. First, as expected, the energy near the heave and rocking modes is reduced as the heave motion is primarily influenced by the use of the controller. Secondly, an increase in energy can be seen near the first pitch mode. This result is also expected, as the implemented controller does not directly damp out the pitch mode, and it is well-known that if a trailing-edge control device is used for GLA, then the root bending moment is alleviated at the expense of the amplification in the root torsional moment [23]. Furthermore, deflecting the aileron not only induces a heave motion but also a pitching motion around the elastic axis. Finally, no differences can be observed for frequencies greater than 15 Hz owing to the finite actuator bandwidth.

4.2. Flutter Suppression Results

In this subsection, we evaluated the flutter suppression ability of the INDI controller designed in Equation (19). Both the control architecture and control gains ($K_p,K_d$) were kept the same as in the GLA experiment. It is noteworthy that the control effectiveness $\overline{\mathcal{B}}\left({\mathit{x}}_0\right)$ is a function of dynamic pressure. Therefore, it was scaled by $U^2$ during the experiment. The static flow velocity U was known and was converted from the tunnel’s RPM value. The open-loop flutter speed was determined to be $U_f=14.5\mathrm{m}/\mathrm{s}$ using the parametric flutter margin method [19]. Moreover, the control reversal speed was also determined for this configuration as $U_R=14.4$ m/s, which was slightly below the flutter speed.

To test the performance of the controller for flutter suppression, the flow velocity was gradually increased in steps past the flutter speed up to a velocity of 18.5 m/s. At each velocity step, the wing section was excited from its equilibrium position to trigger flutter. For consistent excitations of the wing section, it was subjected to a predefined $1-$ cosine gust. First, open-loop flutter was recorded. To prevent any damage to the wing section or test setup, the wing section was manually stopped and returned to its equilibrium position after flutter occurred. The manual stopping of the wing was visible in the plotted gust responses by clipping and sharp peaks in the heave responses seen in Figure 10 and Figure 11. After recording the open-loop flutter, the controller was activated, and the wing section was excited again to record the closed-loop flutter response.

Figure 10, Figure 11, Figure 12 and Figure 13 show a comparison of the open- and closed-loop heave response and the commanded and actual aileron deflection for the closed-loop response for increasing flow velocities past the open-loop flutter velocity. All open-loop responses to the excitation start with a 5 mm amplitude and expectedly oscillate with a frequency around 4 Hz, again close to the first heave mode. Figure 10 and Figure 11 show a clear diverging open-loop response. The closed-loop responses show the reduction of the heave response to the initial excitation compared to the open-loop response, after which the disturbance is damped out within approximately one period.

Similar to the GLA results, the actual aileron deflection shows a slight lag in time of 0.02 s. For the higher velocities, Figure 10 and Figure 11, a nonzero aileron deflection can be observed after the initial disturbance has been damped out. Since the test conditions are beyond the control reversal speed, the heave response, however, remains constant, indicating that the controller has found a new equilibrium position where the increase in lift due to the aileron deflection is offset by the decrease in lift due to the change in pitch of the wing section.

Overall, the flutter suppression tests have shown the implemented INDI controller is able to increase the closed-loop flutter speed to 16.8 m/s, an increase of 15.9%.

5. Conclusions

In this paper, we designed a single nonlinear controller for the gust load alleviation and active flutter suppression problems of a typical wing section. The effectiveness of the proposed incremental nonlinear dynamic inversion (INDI) controller was validated by wind tunnel tests.

The GLA performance of the INDI controller was tested by subjecting the wing section to gust frequencies of 3, 3.5, 4, 4.5, and 5 Hz. Time-domain analyses of the results showed a reduction in peak heave displacement of up to 46.7% and a reduction in heave RMS of up to 72.9%. Frequency-domain results showed a decrease in energy near the first heave and rocking modes compared to the open-loop response, while an increase in energy was observed near the first pitch mode. This was expected as the controller focused on damping out the heave motion and it is well-known that the use of a trailing-edge control device typically amplifies pitch motions.

The effectiveness of the proposed INDI controller on flutter suppression was also validated by wind tunnel experiments. Neither the control architecture nor the control parameters need to be changed. The open-loop flutter speed was found at $U_f=13.5$ m/s. Using the INDI controller, the flutter speed was increased to 16.8 m/s, achieving an increase of 15.9%.

6. Outlook

In this research, we focused on alleviating gust responses and suppressing flutter in the heave degrees of freedom. We have planned future wind tunnel tests for including the pitch degree of freedom in the feedback signal, and we will explore the harmonious usage of both aileron and spoiler. Furthermore, we plan to include an automatic adaptation to the INDI controller, which will enable it to self-adapt to variations in free-streaming velocities.