Adaptive Feedforward Control for Gust-Induced Aeroelastic Vibrations

: This paper demonstrates the implementation of an adaptive feedforward controller to reduce structural vibrations on a wing typical section. The aeroelastic model includes a structural nonlinearity, which is modelled in a polynomial form. Aeroelastic vibrations are induced by several gusts and atmospheric turbulence, including the discrete “one-minus-cosine” and a notably good approximation in the time-domain to the von K á rm á n spectrum. The control strategy based on the adaptive feedforward controller has several advantages compared to the standard feedback controller. The controller gains, which are updated in real-time during the gust encounter, are found solving a minimization problem using the ﬁnite impulse responses as basis functions. To make progress with the application in aeroelasticity, a single-input single-output controller is designed measuring the wing torsional deformation. For both deterministic and random atmospheric shapes, the controller was found successful in alleviating the aeroelastic vibrations. The impact of the control action on the unmeasured structural modes was found minimal.


Introduction
In flight, aircraft regularly encounter atmospheric turbulence. Turbulence, lightning, hail and other phenomena can lead to injuries and discomfort on board and damage to the aircraft [1], resulting in huge cost to airlines. Poor weather detection and analysis can result in poor pilot decision making which could lead to otherwise completely avoidable danger to flights [2]. In addition, weather-related delays and cancellations cost airlines millions of dollars and cost countries' economies billions of dollars in lost productivity each year [3].
Models of atmospheric turbulence for aircraft gust load analysis (GLA) have been developed over time and are today required by certification authorities. The most common models of turbulence rely on a discrete and continuous representation of the gust. More details may be found in Ref. [4]. One of most serious forms of turbulence in flight is clear air turbulence (CAT). In 1966, a National Committee for CAT officially defined CAT as "all turbulence in the free atmosphere of interest in aerospace applications that is not in or adjacent to visible convective activity". Over time, less formal definitions of CAT have appeared but the most comprehensive definition is "turbulence encountered outside of convective clouds". CAT was recognized as a problem with the advent of high altitude jet operations in the 1950 s. CAT is particularly problematic because it is often encountered unexpectedly and frequently without visual clues to warn pilots of the hazard [5]. CAT is very difficult to predict accurately, due in part to the fact that CAT is spotty in both dimensions and time. Commonly, a turbulent area associated with a jet-stream extends in the order of 100 to 300 miles elongated in the direction of the wind, 50 to 100 miles wide and 5000 feet deep. These areas may persist from 30 min to a day. Despite the difficulty in forecasting CAT, there are certain rules that have been developed to identify those areas where CAT formation is likely.
The control for GLA emerged vigorously in 1964 due to a significant incident. During a routine mission, a B52-H bomber encountered severe turbulence with an estimated peak velocity of 35 m/s [4]. Nearly 80% of the ventral fin broke off in flight, resulting in an extensive investigation and several development programs on GLA. A wide range of feedback control strategies were used to minimize the adverse effects induced by gusts. The reader is invited to consult, as a sample of the vast literature on the topic, fundamental activities as those reported in [6,7]. In general, these techniques fed the structural response, sensed by accelerometers, to the controller that calculated the commanded control action. Today, strategies based on the linear quadratic regulator and derived from optimal control are considered mature in light of the extensive applications in aero-servo-elasticity. The drawback of a feedback-based controller is caused by the time delay between the sensed information, which must propagate through the entire plant before exhibiting the effect to be measured and the determination of the control action. As investigated in [8], this time delay is critical for the stability of the plant. This major disadvantage may be overcome by a feedforward control strategy when an a-priori knowledge of the disturbance is available. In principle, feedforward control can eliminate the effects of the measured disturbance on the process output. In presence of modelling errors, feedforward control can often reduce the effect of the measured disturbance on the output better than that achievable by feedback control alone [9]. An advantage of feedforward control is that there is no time delay between the disturbance and the control compensation and a corrective action can be taken before the output deviates from the set point [8,10].
The capability to measure the turbulence ahead of the aircraft is critical to the use of an on-board feedforward controller for GLA. Two types of sensors are available. The first is a light detection and ranging (LIDAR) sensor [11]. As an example, Honeywell offers a family of advanced weather radar systems (IntuVue ® , Primus ® , Honeywell International Inc., Morris Plains, NJ, USA) that scan the sky ahead of the aircraft at 17 different tilt angles. Several airliners have equipped their fleets with this equipment to ensure smoother, more comfortable and safer flights (This is a relevant subject as reported at: https://www.aviationtoday.com/2018/07/24/boeing-testing-use-autonomy-lidar-futureair-cargo-aircraft/ (accessed on 3 August 2018)). Using the reference system from a LIDAR sensor, Reference [12] developed an adaptive feedforward controller for the GLA of a linear aeroelastic model of the F/A-18 aircraft. The orthogonal finite impulse response (FIR) filter was used as the adaptive controller, which exhibited a good performance. The second type of sensor is referred to as an alpha probe [13]. Reference [14] presented a single-input single-output (SISO) FIR feedforward controller to suppress the gust-induced vibrations of a transport airplane. A similar strategy was extended to the multiple-input multiple-output case in [15]. Reference [16] investigated a hybrid strategy obtained by combining feedback and feedforward controls. It was found that the hybrid strategy significantly improved the reduction in wing bending moment when compared with the feedback strategy alone. Flight tests of an adaptive feedforward controller were carried out to alleviate wing bending vibrations [17]. Using a nose boom mounted flight log sensor, the feedforward compensation was found to provide a good reduction of dynamic loads and an improvement in ride comfort.
The already limited literature on the topic becomes scarce when nonlinearities exist in the plant to be controlled. Reference [18] investigated the impact that structural nonlinearities have on the effectiveness of a feedforward controller for a wing typical section. It was found that to suppress gust-induced vibrations of an intrinsically nonlinear plant, the control performance was degraded when using a linear representation for the internal model. The control strategy performed well when including all the nonlinearities in the model. Based on this forerunner knowledge, this paper aims at investigating the deployment of an adaptive feedforward control strategy for GLA of a nonlinear aeroelastic plant. The test case is for a model of a wing typical section that exhibited sub-critical limit cycle oscillations in wind tunnel tests. The work is built around three objectives. The first objective lies in the derivation of the aeroelastic equations in the time domain and in the numerical implementation of the adaptive feedforward controller. The second objective relates to the identification of the plant dynamics. The last objective assesses the effects of the proposed controller for GLA. It is worth noting that a follow-up study conducted flight tests on a high-altitude, long-endurance unmanned air vehicle with a wing span of 24 m [19].
The paper continues with a description of the aeroelastic testbed in Section 2. Models of atmospheric turbulence and gusts are given in Section 3. The control theory underlying the adaptive feedforward strategy is reviewed in Section 4. Then, numerical results are presented in Section 5. Concluding remarks are finally given in Section 6.

Aeroelastic Model
In this work, the aeroelastic test case is for a wing typical section, Figure 1. This is modelled as a two-dimensional symmetric airfoil section and the structural motion is characterized by two degrees of freedom: the plunge, h is measured at the elastic axis (e.a.) and is positive downward; and the pitch rotation, α, is the rotation of the airfoil section about the elastic axis, positive nose-up. As common in classic aeroelasticity, b is the airfoil semi-chord. The elastic axis is located at a distance a h b from the mid-chord, positive when the elastic axis is aft of the mid-chord. The center of gravity (c.g.) is located at x α b aft of the mid-chord. The aeroelastic model is restrained to the ground through a system of elastic springs, K α , K ξ and dampers, C α and C h (omitted in the figure for clarity), acting on both degrees of freedom. The nonlinearity in the stiffness of the elastic springs is modelled in a polynomial form including terms up to fifth order. lies in the derivation of the aeroelastic equations in the time domain and in the numerical implementation of the adaptive feedforward controller. The second objective relates to the identification of the plant dynamics. The last objective assesses the effects of the proposed controller for GLA. It is worth noting that a follow-up study conducted flight tests on a high-altitude, longendurance unmanned air vehicle with a wing span of 24 m [19]. The paper continues with a description of the aeroelastic testbed in Section 2. Models of atmospheric turbulence and gusts are given in Section 3. The control theory underlying the adaptive feedforward strategy is reviewed in Section 4. Then, numerical results are presented in Section 5. Concluding remarks are finally given in Section 6.

Aeroelastic Model
In this work, the aeroelastic test case is for a wing typical section, Figure 1. This is modelled as a two-dimensional symmetric airfoil section and the structural motion is characterized by two degrees of freedom: the plunge, ℎ is measured at the elastic axis ( . .) and is positive downward; and the pitch rotation, α, is the rotation of the airfoil section about the elastic axis, positive nose-up. As common in classic aeroelasticity, is the airfoil semi-chord. The elastic axis is located at a distance a b from the mid-chord, positive when the elastic axis is aft of the mid-chord. The center of gravity ( . .) is located at aft of the mid-chord. The aeroelastic model is restrained to the ground through a system of elastic springs, , and dampers, and (omitted in the figure for clarity), acting on both degrees of freedom. The nonlinearity in the stiffness of the elastic springs is modelled in a polynomial form including terms up to fifth order.
Herein, the interest is to analyze the dynamic response of the system induced by an external disturbance. This disturbance is representative of atmospheric turbulence and gusts, which are treated in more detail in Section 3. A flap is mounted at the trailing-edge of the airfoil, at a distance aft of the mid-chord, for GLA. The flap deflection, , is defined relative to the airfoil chord. It is assumed that the system has a horizontal equilibrium position at (ℎ = 0, = 0, = 0). The equations of motion governing the dynamics of the two degrees of freedom system are ℎ + + + ℎ + ℎ + ℎ = − ℎ + + + + + = (1) The coupled system of equations consists of the structural dynamics on the left-hand side and of aerodynamic loads on the right-hand side. The term (in kg) indicates the total mass that undergoes a vertical motion, (in kg·m) is the first moment of inertia about the elastic axis and (in kg·m 2 ) the second moment of inertia about the elastic axis. The lift, (in N/m), is defined positive upward as usual in aerodynamics, whereas the plunge deflection, ℎ , is positive downward as Figure 1. Schematic of an airfoil section with trailing-edge flap; the wind velocity is to the right and horizontal; e.a. and c.g. denote, respectively, the elastic axis and the center of gravity.
Herein, the interest is to analyze the dynamic response of the system induced by an external disturbance. This disturbance is representative of atmospheric turbulence and gusts, which are treated in more detail in Section 3. A flap is mounted at the trailing-edge of the airfoil, at a distance cb aft of the mid-chord, for GLA. The flap deflection, δ, is defined relative to the airfoil chord.
It is assumed that the system has a horizontal equilibrium position at (h = 0, α = 0, δ = 0). The equations of motion governing the dynamics of the two degrees of freedom system are The coupled system of equations consists of the structural dynamics on the left-hand side and of aerodynamic loads on the right-hand side. The term m (in kg) indicates the total mass that undergoes a vertical motion, S α (in kg·m) is the first moment of inertia about the elastic axis and I α (in kg·m 2 ) the second moment of inertia about the elastic axis. The lift, L (in N/m), is defined positive upward as usual in aerodynamics, whereas the plunge deflection, h, is positive downward as conventionally done in aeroelasticity. Hence, the negative sign in front of the lift. The aerodynamic moment, M (in N), is taken about the elastic axis.
The airfoil semi-chord, b and the uncoupled pitch natural frequency, ω α , are used to rewrite Equation (1) in terms of non-dimensional parameters. In non-dimensional form, the equations become It is worth noting that the derivative with respect to physical time, .
x = dx dt , in Equation (1) is replaced by the derivative with respect to non-dimensional time, The aerodynamics is for an incompressible two-dimensional flow [20]. The total aerodynamic loads consist of contributions arising from the airfoil motion, the flap deflection and the encounter with a time-varying gust The generalization of the aerodynamic loads to an arbitrary input time-history is obtained through convolution (or Duhamel integral). For a practical evaluation of the Duhamel integral, an exponential approximation is used for the Wagner and Küssner functions, which describe, respectively, the indicial build-up of the circulatory part of the lift and the lift build-up for the penetration into a sharp-edged gust [21]. This implies that the governing equations-Equation (2)-are a set of integro-differential equations (IDEs) when an expression for the aerodynamic loads is introduced.

State Space Form
The governing equations are a set of IDEs which are difficult to solve with state-of-the-art algorithms, primarily developed for ordinary differential equations (ODEs). The mathematical procedure to convert the set of IDEs into a set of ODEs is based on defining additional variables and equations describing the evolution of the aerodynamics. Following [22], the aeroelastic equations in state space form are rewritten as where the dependence on non-dimensional time, τ, is omitted for brevity. In a matrix-vector form the state space vector is The function f(x, τ) contains the (nonlinear) elastic contributions from the springs and is a nonlinear function of the state vector, x The term µ(τ) = δ indicates the control input (rotation of the trailing-edge flap) and d(τ) = w g denotes the exogenous external disturbance (atmospheric turbulence and gust). The elements of the state space vector w 1 through w g are aerodynamic states and a full derivation is given [22].

Atmospheric Turbulence and Gust Models
Models used for the prediction of the aircraft response need to accommodate those events that are perceived as discrete and usually described as gusts, as well as the phenomena described as continuous turbulence [4].

Discrete Deterministic Gust
The most common model of a discrete deterministic gust, which has evolved over the years from the isolated sharp-edge gust function in the earliest airworthiness requirements, is the "one-minus-cosine" function. Its formulation is where U ds is the design gust velocity and H is the gust gradient distance, that is, the distance over which the gust velocity increases to a maximum value. Airworthiness requirements prescribe a relationship between the design gust velocity and the gust gradient distance.

Continuous Atmospheric Turbulence
The engineering model of random turbulence at altitude has been developed over many years. It is now widely accepted that it is satisfactory to treat atmospheric turbulence as frozen, homogeneous and isotropic in relatively large patches. The frozen-field assumption, closely allied to Taylor s hypothesis, is that turbulent velocities do not change during the time of passage of the airplane. This is a valid assumption in most cases. Dryden and von Kármán models are considered adequate engineering models to predict the correlation and spectra, with the weight of experimental evidence favoring the latter.
A commonly used spectrum that matches experimental data is von Kármán model. The power spectral density (PSD, in m 2 /(s 2 Hz)) for the vertical direction, Φ z , according to the Military Specification MIL-F-8785C [23], is given by where Ω = ω/U ∞ is the scaled frequency (in rad/m), σ z the turbulence intensity (in m/s), U ∞ the freestream speed (in m/s) and L z the turbulence scale length (in m). The turbulence intensity is related to the wind speed at an altitude of 20 ft (6 m) by Values of wind speed are tabulated based on the turbulence severity: for light turbulence, w 20 = 15 kn; for moderate turbulence, w 20 = 30 kn; and for severe turbulence, w 20 = 45 kn.

Numerical Implementation
A time-domain representation of the external disturbance is required to simulate the non-linear dynamics of a system in the time-domain. Because of the irrational nature of von Kármán spectrum (in the frequency-domain, the power 11/6 is not an integer), there is no exact analytical formulation in the time-domain corresponding to Equation (8).
The approach used in this work consists of the following steps. First, the Fourier transform of a unit variance band-limited white noise signal, X(Ω), is calculated. This is passed through a filter, H z (Ω), which is defined as the square root of the PSD spectrum in Equation (8). Then, the output signal is calculated, given the filter and the input signal Finally, the inverse Fourier transform of the output signal, W g (Ω), is performed to obtain the continuous turbulence in the time-domain, w g (t), with statistical properties matching von Kármán spectrum.
The procedure requires calculating the Fourier transform and its inverse, requiring more computational resources as the time window increases but it was found to be more robust and accurate than other approaches. The method described herein is available in an open-source toolbox (both in MatLAB and Python), which is referred to as the Von Kármán Turbulence Generator (VKTG). The VKTG toolbox implements the mathematical representations of random turbulence defined in the Military Specification MIL-F-8785C and Military Handbook MIL-HDBK-1797, allowing for the dependence of the root mean square turbulent velocity and turbulence length scale on aircraft mission parameters and weather conditions. Figure 2 shows that, at higher frequencies, the PSD of the VKTG model achieves a better agreement with von Kármán spectrum of Equation (8) than the off-the-shelf MATLAB/SIMULINK model. For more information, the reader is invited to refer to Ref. [24] and the on-line documentation (http://daronchlab.com/ (Retrieved 10 August 2018)). The approach used in this work consists of the following steps. First, the Fourier transform of a unit variance band-limited white noise signal, (Ω), is calculated. This is passed through a filter, (Ω), which is defined as the square root of the PSD spectrum in Equation (8). Then, the output signal is calculated, given the filter and the input signal Finally, the inverse Fourier transform of the output signal, (Ω), is performed to obtain the continuous turbulence in the time-domain, ( ), with statistical properties matching von Kármán spectrum.
The procedure requires calculating the Fourier transform and its inverse, requiring more computational resources as the time window increases but it was found to be more robust and accurate than other approaches. The method described herein is available in an open-source toolbox (both in MatLAB and Python), which is referred to as the Von Kármán Turbulence Generator (VKTG). The VKTG toolbox implements the mathematical representations of random turbulence defined in the Military Specification MIL-F-8785C and Military Handbook MIL-HDBK-1797, allowing for the dependence of the root mean square turbulent velocity and turbulence length scale on aircraft mission parameters and weather conditions. Figure 2 shows that, at higher frequencies, the PSD of the VKTG model achieves a better agreement with von Kármán spectrum of Equation (8) than the off-the-shelf MATLAB/SIMULINK model. For more information, the reader is invited to refer to Ref. [24] and the on-line documentation (http://daronchlab.com/ (Retrieved 10 August 2018)).

Adaptive Feedforward Control
A typical block diagram of an adaptive feedforward control algorithm consists of two channels, Figure 3. The disturbance path includes the transfer function of the physical plant to be controlled, , between the external disturbance, and the system response, . The control path contains the transfer function of the adaptive feedforward controller, , to be designed. The sensed reference

Adaptive Feedforward Control
A typical block diagram of an adaptive feedforward control algorithm consists of two channels, Figure 3. The disturbance path includes the transfer function of the physical plant to be controlled, H, between the external disturbance, w g and the system response, x. The control path contains the transfer function of the adaptive feedforward controller, G c , to be designed. The sensed reference signal,ŵ g , is fed forward to the adaptive controller to generate the control command. The block G indicates the transfer function of the physical plant between the control input, u and the system response predicted by the model, y. The control input to the physical plant is commanded by the controller. Finally, the block with transfer functionĜ represents an approximation of −G. As an input to system G c , the filter output u drives the flap to produce a cancellation force and moment. The error signal, e, is the sum of the disturbance path output, x and the control path output, y. It is assumed that the wing typical section is equipped with an onboard LIDAR to provide the GLA system with a lead-time reference signal for the vertical gust. The objective of the GLA control is to adjust the weight coefficients of the controller to minimize the error signal . Furthermore, the dynamic characteristics of the actuator driving the rotation of the flap were not accounted. This decision was made based on the experimental observation that the dynamics of the actuator was fast enough, compared to the characteristic times of the aeroelastic plant, to be neglected [22]. Their inclusion is straightforward [25] and their effect would appear within the transfer function . When there is an exact knowledge of the system transfer functions, and , setting the error to be null yields the ideal feedforward controller = − As this is rarely the case in practice, the transfer functions of the plant are approximated using system identification (SI) methods. Generally, system identification methods seek to create a mathematical relationship between the output response and the input signal. These methods enjoy widespread utilization because the model creation is straightforward. A major drawback is the limited validity of predictions, restricted by the characteristics of the input signal, often referred to as the training signal, used to generate the model in first place [26]. Section 5.3 discusses the choice of the training signal based on a chirp input and evaluates the accuracy of the approximated transfer functions.
Once the system transfer functions and are identified, the block with transfer function is approximated as The output signal of the block is then calculated On the disturbance path, the system response ( ) is Substituting Equation (13) into Equation (14) and then substituting from the resulting expression into Equation (15) yields It is assumed that the wing typical section is equipped with an onboard LIDAR to provide the GLA system with a lead-time reference signal for the vertical gust. The objective of the GLA control is to adjust the weight coefficients of the controller to minimize the error signal e. Furthermore, the dynamic characteristics of the actuator driving the rotation of the flap were not accounted. This decision was made based on the experimental observation that the dynamics of the actuator was fast enough, compared to the characteristic times of the aeroelastic plant, to be neglected [22]. Their inclusion is straightforward [25] and their effect would appear within the transfer function G.
When there is an exact knowledge of the system transfer functions, G and H, setting the error to be null yields the ideal feedforward controller As this is rarely the case in practice, the transfer functions of the plant are approximated using system identification (SI) methods. Generally, system identification methods seek to create a mathematical relationship between the output response and the input signal. These methods enjoy widespread utilization because the model creation is straightforward. A major drawback is the limited validity of predictions, restricted by the characteristics of the input signal, often referred to as the training signal, used to generate the model in first place [26]. Section 5.3 discusses the choice of the training signal based on a chirp input and evaluates the accuracy of the approximated transfer functions. Once the system transfer functions G and H are identified, the block with transfer functionĜ is approximated asĜ = −G (13) The output signal of the blockĜ is then calculated On the disturbance path, the system response x(t) is Substituting Equation (13) into Equation (14) and then substituting w g from the resulting expression into Equation (15) yields This relation allows identifying the feedforward controller using a map between u a (t) and x(t) [27]. An adaptive filtering strategy to update the controller gains is adopted to increase the robustness of the weight coefficients of the controller against errors in the measurement of the external disturbance and nonlinear effects in the system response.

Control Model
The controller is considered as a discrete linear time-invariant system where G c (q) indicates the controller transfer operator and q is the forward shift operator, qw g (t) = w g (t + 1). The corresponding transfer function, G c (z) with z ∈ C, is formulated as A transfer function for a controller of order n is modelled as a sum of coefficients to be determined, L k , multiplied by an adequate choice of basis functions, B k (z). In this work, the basis functions are given by a FIR set B k (z) = z −k , k = 1, 2, . . . , n which captures any shape of response up to a frequency limit.

Exponentially Weighted Recursive Least-Square Algorithm
The calculation of the coefficients of the basis functions, L k , in Equation (18) is performed using an exponentially weighted recursive least-square algorithm. At the i-th time step, t i , the error between the desired response, e(i) and the FIR model output, r(i), is introduced where L(i) = (L 1 (i), L 2 (i), . . . , L n (i)) T is the unknown vector of coefficients (often referred to as tap weight vector) and Φ(i) = (u a 1 (i), u a 2 (i), . . . , u a n (i)) T is the vector containing the output of every basis function. A cost function is then defined as follows where N T indicates the total number of time steps and λ ∈ (0, 1] is the forgetting factor. The adaptive algorithm to calculate the coefficients of the basis functions follows the steps: • Initialization: parameters are initialized before entering the time iteration loop. L(0) = 0 is initialized to a column vector of zeros. P(0) = δ −1 I is the inverse correlation matrix (where δ = 0.1 in this case). • Time iteration: the recursive algorithm is repeated at each time step, t i , for the entire duration of the time vector. For i = 1, 2, . . . , N, calculate The algorithm is sensitive to the choice of the forgetting factor, λ. In this work, the value is set to 1. The reader is referred to Reference [28] for more details.

Results
The aeroelastic solver has previously been validated in many numerical and experimental studies. The reader is invited to refer to References [21,22,29] for further information. Based on previous findings, a non-dimensional time step of 0.05 was deemed sufficient to capture the dynamics of the system.
The aeroelastic parameters in this study model a wind tunnel configuration tested in low speed flow at the University of Liverpool, Figure 4. The experimental model consists of a rigid wing that spans the entire width of the wind tunnel cross-section and is suspended by a system of springs to the tunnel walls. The wing cross-section is based on the NACA0018 airfoil; the semi-chord is b = 0.175 m and the pitch natural frequency is ω α = 28.061 rad/s. The aeroelastic parameters are summarized in Table 1 and the procedure for the model updating is discussed in Reference [22]. The last two entries in the table are the cubic and quintic coefficients of the polynomial expansion of the stiffness in the plunge degree of freedom. The values indicate a typical hardening non-linearity.

Results
The aeroelastic solver has previously been validated in many numerical and experimental studies. The reader is invited to refer to References [21,22,29] for further information. Based on previous findings, a non-dimensional time step of 0.05 was deemed sufficient to capture the dynamics of the system.
The aeroelastic parameters in this study model a wind tunnel configuration tested in low speed flow at the University of Liverpool, Figure 4. The experimental model consists of a rigid wing that spans the entire width of the wind tunnel cross-section and is suspended by a system of springs to the tunnel walls. The wing cross-section is based on the NACA0018 airfoil; the semi-chord is = 0.175 m and the pitch natural frequency is = 28.061 rad/s. The aeroelastic parameters are summarized in Table 1 and the procedure for the model updating is discussed in Reference [22]. The last two entries in the table are the cubic and quintic coefficients of the polynomial expansion of the stiffness in the plunge degree of freedom. The values indicate a typical hardening non-linearity. The wind tunnel test rig exhibited an interesting non-linear dynamic behavior, characterized by the occurrence of sub-critical limit cycle oscillations (LCOs). From a control standpoint, it represents a challenging test case to demonstrate the control and suppression of gust-induced aeroelastic vibrations and LCOs.

Frequency-Domain Analysis
First, the aeroelastic stability is analyzed. Figure 5 shows the damped frequency and damping ratio for varying freestream speed. The solid line in figure is from numerical predictions that were  The wind tunnel test rig exhibited an interesting non-linear dynamic behavior, characterized by the occurrence of sub-critical limit cycle oscillations (LCOs). From a control standpoint, it represents a challenging test case to demonstrate the control and suppression of gust-induced aeroelastic vibrations and LCOs.

Frequency-Domain Analysis
First, the aeroelastic stability is analyzed. Figure 5 shows the damped frequency and damping ratio for varying freestream speed. The solid line in figure is from numerical predictions that were obtained solving an eigenvalue problem of the linearized aeroelastic model. Two sets of experimental data, labelled in figure by "WT Data", are available: one set resulting from control surface excitation and the other from shaker excitation of the plunge degree of freedom. At intermediate speeds, predictions of damped frequency agree well with measurements. A reason for the agreement degrading at lower speed and close to the flutter speed is the complication in conducting the experiments. At lower speeds, the shaker excitation is used to excite the modes, which in turn may affect the free vibrations of the test rig. Close to flutter, a difficulty is the coalescence of the pitch and plunge frequencies.
A reasonable agreement is observed for the aeroelastic damping ratio. Simulations predict a flutter speed of U f = 15.28 m/s.

Time-Domain Analysis
The hardening nonlinearity in the plunge degree of freedom, Table 1, was found to significantly affect the dynamics of the system. The aeroelastic system exhibits subcritical LCOs and the amplitude and frequency of LCOs were found in good agreement between measurements and predictions. One example is in Figure 6 where the response to an initial perturbation in angle of attack is simulated at two subcritical speeds (8 and 13 m/s). The onset of LCOs is at a speed of 12.871 m/s. The interest in this work is focused at testing the control strategy of Section 4 at a speed = 8 m/s below the onset point of LCOs.

Time-Domain Analysis
The hardening nonlinearity in the plunge degree of freedom, Table 1, was found to significantly affect the dynamics of the system. The aeroelastic system exhibits subcritical LCOs and the amplitude and frequency of LCOs were found in good agreement between measurements and predictions. One example is in Figure 6 where the response to an initial perturbation in angle of attack is simulated at two subcritical speeds (8 and 13 m/s). The onset of LCOs is at a speed of 12.871 m/s. The interest in this work is focused at testing the control strategy of Section 4 at a speed U ∞ = 8 m/s below the onset point of LCOs.

Time-Domain Analysis
The hardening nonlinearity in the plunge degree of freedom, Table 1, was found to significantly affect the dynamics of the system. The aeroelastic system exhibits subcritical LCOs and the amplitude and frequency of LCOs were found in good agreement between measurements and predictions. One example is in Figure 6 where the response to an initial perturbation in angle of attack is simulated at two subcritical speeds (8 and 13 m/s). The onset of LCOs is at a speed of 12.871 m/s. The interest in this work is focused at testing the control strategy of Section 4 at a speed = 8 m/s below the onset point of LCOs.

Model Identification
The system transfer function, G, is identified by post-processing the input-output relationship. The input, which describes the time evolution of the trailing-edge flap, is a chirp signal defined as where u 0 is a constant value and u A is the amplitude of the sinusoidal motion. The frequency of the sinusoidal motion, f , varies linearly in time and covers the frequency range [ f m , f M ] over the total simulation time, T. The frequency range is chosen to adequately excite the system over the desired frequency range. In this work, the parameters of the chirp signal were set to u 0 = 0.0 • and u A = 1.0 • . The aeroelastic model has two dominant natural frequencies at 2.646 Hz (ω h = 16.629 rad/s, plunge) and 4.466 Hz (ω α = 28.061 rad/s, pitch). For increasing dynamic pressure, the pitch and plunge frequencies tend to coalesce, Figure 5 and so the frequency range covered by the chirp signal was chosen between f m = 0.01 Hz and f M = 8.0 Hz. The input signal employed for the model identification is illustrated in Figure 7. The aeroelastic response induced by the control surface is shown in Figure 8. The resonant behavior in both degrees of freedom confirms the adequateness of the chosen frequency range.
A SISO model is considered in this exploratory work. The pitch degree of freedom is taken as the measured output, whereas the plunge is assumed unmeasurable. The transfer function, G, between the control effector and the pitch response is identified using the prediction error minimization algorithm available in MATLAB. A preliminary study was carried out to ensure independence from the number of zeros and poles used in constructing the transfer function. It was found that six zeros and seven poles, which are listed in Table 2, provide a good approximation of the aeroelastic system composed of 12 states in total. The real part of all poles is negative and the resulting transfer function is therefore stable. As confirmed in Figure 9, the identified transfer function (with six zeros and seven poles) is identical to the transfer function that relates the Laplace transform of the system output to the Laplace transform of the (chirp) input.
chosen to adequately excite the system over the desired frequency range.
In this work, the parameters of the chirp signal were set to = 0.0° and = 1.0°. The aeroelastic model has two dominant natural frequencies at 2.646 Hz ( = 16.629 rad/s, plunge) and 4.466 Hz ( = 28.061 rad/s, pitch). For increasing dynamic pressure, the pitch and plunge frequencies tend to coalesce, Figure 5 and so the frequency range covered by the chirp signal was chosen between = 0.01 Hz and = 8.0 Hz. The input signal employed for the model identification is illustrated in Figure 7. The aeroelastic response induced by the control surface is shown in Figure 8. The resonant behavior in both degrees of freedom confirms the adequateness of the chosen frequency range.  chosen to adequately excite the system over the desired frequency range. In this work, the parameters of the chirp signal were set to = 0.0° and = 1.0°. The aeroelastic model has two dominant natural frequencies at 2.646 Hz ( = 16.629 rad/s, plunge) and 4.466 Hz ( = 28.061 rad/s, pitch). For increasing dynamic pressure, the pitch and plunge frequencies tend to coalesce, Figure 5 and so the frequency range covered by the chirp signal was chosen between = 0.01 Hz and = 8.0 Hz. The input signal employed for the model identification is illustrated in Figure 7. The aeroelastic response induced by the control surface is shown in Figure 8. The resonant behavior in both degrees of freedom confirms the adequateness of the chosen frequency range.  A SISO model is considered in this exploratory work. The pitch degree of freedom is taken as the measured output, whereas the plunge is assumed unmeasurable. The transfer function, , between the control effector and the pitch response is identified using the prediction error minimization algorithm available in MATLAB. A preliminary study was carried out to ensure independence from the number of zeros and poles used in constructing the transfer function. It was found that six zeros and seven poles, which are listed in Table 2, provide a good approximation of the aeroelastic system composed of 12 states in total. The real part of all poles is negative and the resulting transfer function is therefore stable. As confirmed in Figure 9, the identified transfer function (with six zeros and seven poles) is identical to the transfer function that relates the Laplace transform of the system output to the Laplace transform of the (chirp) input.     First, the search for the worst-case-gust was carried out for the "one-minus-cosine" family considering a range of the gust gradient distance, H/b ∈ [10,200]. The search, which is an inexpensive process here for the small size of the numerical model, employed a standard Kriging interpolation as response surface method [30]. The worst-case-gust was identified to have a gust gradient distance H/b = 20. The gust intensity was set to a constant value, U ds /U ∞ = 0.1.
Then, the aeroelastic response caused by the worst-case-gust encounter was simulated with and without the adaptive feedforward controller. The comparison is reported in Figure 10. For the pitch, the response with control reveals a significant reduction in the peak-to-peak rotation during and immediately after the interaction with the gust. A reduction of the peak-to-peak rotation up to 50% is recorded. As the gust moves away from the surface of the wing typical section, the improvement in the loads alleviation vanishes because the effective reference input for the control path is null. For the plunge, there is no improvement in the response with control compared to the aeroelastic response without control. This is not unexpected because the controller is a SISO system, with the pitch rotation being the sole measurable quantity for the GLA.  Figure 11 shows the required control action for GLA and a sample of three weight coefficients of the controller. The trailing-edge flap is commanded during the gust encounter, whereas it is unactuated as the gust perturbation moves away from the wing typical section. As the gust propagates along the airfoil, the pitch rotation increases, as shown in Figure 10a. To counter the increasing pitch rotation, the control strategy commands the trailing-edge flap down, that is, a positive control input as shown in Figure 11. This control action has two effects. First, the lift increases, increasing in turn the vertical plunge deflection which is negative, see Figure 10b, according to the sign convention of Section 2. Second, deflecting the flap moves the center of pressure downstream aft of the elastic axis and the increased lift causes a pitch-down rotation. This highlights the capability of the controller to exploit structural flexibility effects to reduce gust-induced aeroelastic vibrations. It is also worth noting that the flap rotation and the angular speed meet fully physical limitations of the experimental test rig: the maximum allowed rotation is ±7° and the speed  Figure 11 shows the required control action for GLA and a sample of three weight coefficients of the controller. The trailing-edge flap is commanded during the gust encounter, whereas it is unactuated as the gust perturbation moves away from the wing typical section. As the gust propagates along the airfoil, the pitch rotation increases, as shown in Figure 10a. To counter the increasing pitch rotation, the control strategy commands the trailing-edge flap down, that is, a positive control input as shown in Figure 11. This control action has two effects. First, the lift increases, increasing in turn the vertical plunge deflection which is negative, see Figure 10b, according to the sign convention of Section 2. Second, deflecting the flap moves the center of pressure downstream aft of the elastic axis and the increased lift causes a pitch-down rotation. This highlights the capability of the controller to exploit structural flexibility effects to reduce gust-induced aeroelastic vibrations. It is also worth noting that the flap rotation and the angular speed meet fully physical limitations of the experimental test rig: the maximum allowed rotation is ±7 • and the speed is limited to 15 Hz. A preliminary study was performed to select the order of the controller. It was found that a controller of order 20 was a good compromise between cost and accuracy. The coefficients of the controller were initially trained for 10,000 time-steps using von Kármán turbulence model. Figure 11 shows the required control action for GLA and a sample of three weight coefficients of the controller. The trailing-edge flap is commanded during the gust encounter, whereas it is unactuated as the gust perturbation moves away from the wing typical section. As the gust propagates along the airfoil, the pitch rotation increases, as shown in Figure 10a. To counter the increasing pitch rotation, the control strategy commands the trailing-edge flap down, that is, a positive control input as shown in Figure 11. This control action has two effects. First, the lift increases, increasing in turn the vertical plunge deflection which is negative, see Figure 10b, according to the sign convention of Section 2. Second, deflecting the flap moves the center of pressure downstream aft of the elastic axis and the increased lift causes a pitch-down rotation. This highlights the capability of the controller to exploit structural flexibility effects to reduce gust-induced aeroelastic vibrations. It is also worth noting that the flap rotation and the angular speed meet fully physical limitations of the experimental test rig: the maximum allowed rotation is ±7° and the speed is limited to 15 Hz. A preliminary study was performed to select the order of the controller. It was found that a controller of order 20 was a good compromise between cost and accuracy. The coefficients of the controller were initially trained for 10,000 time-steps using von Kármán turbulence model.

Continuous Atmospheric Turbulence
For the random case, atmospheric turbulence was generated with the VKTG toolbox at an altitude h = 200 m with the intensity set to "moderate". The maximum instantaneous gust intensity was set to 10% of the freestream speed, U ∞ = 8 m/s. The time history of the vertical speed of the turbulence is shown in Figure 12. Figure 13 shows the aeroelastic response with and without control. Similar considerations to those already noted for the discrete gust hold true. The controller is effective at reducing the pitch vibrations throughout the entire duration of the gust encounter.

Continuous Atmospheric Turbulence
For the random case, atmospheric turbulence was generated with the VKTG toolbox at an altitude ℎ= 200 m with the intensity set to "moderate". The maximum instantaneous gust intensity was set to 10% of the freestream speed, = 8 m/s. The time history of the vertical speed of the turbulence is shown in Figure 12. Figure 13 shows the aeroelastic response with and without control. Similar considerations to those already noted for the discrete gust hold true. The controller is effective at reducing the pitch vibrations throughout the entire duration of the gust encounter.    To quantify the impact of the adaptive feedforward controller for GLA in the presence of atmospheric turbulence, the mean and standard deviation of the pitch and plunge time responses are listed in Table 3. When confronted with the uncontrolled aeroelastic response, the statistics of the pitch rotation are significantly improved by the control action: the mean value is reduced by approximately 80% and the standard deviation by approximately 35%. The increase in the statistics of the plunge is not unexpected, as it was assumed unmeasurable for the purpose of GLA. The tradeoff between the reduction in pitch statistics and the increase in plunge statistics is yet in favor of the former.
For the wing typical section, the pitch and plunge are structural modes. The input to the controller is a measure of the pitch degree of freedom, that is, one structural mode. Therefore, it is To quantify the impact of the adaptive feedforward controller for GLA in the presence of atmospheric turbulence, the mean and standard deviation of the pitch and plunge time responses are listed in Table 3. When confronted with the uncontrolled aeroelastic response, the statistics of the pitch rotation are significantly improved by the control action: the mean value is reduced by approximately 80% and the standard deviation by approximately 35%. The increase in the statistics of the plunge is not unexpected, as it was assumed unmeasurable for the purpose of GLA. The trade-off between the reduction in pitch statistics and the increase in plunge statistics is yet in favor of the former.
For the wing typical section, the pitch and plunge are structural modes. The input to the controller is a measure of the pitch degree of freedom, that is, one structural mode. Therefore, it is assumed that the sensor measures solely the torsional mode. An alternative approach to mitigate the effects of a SISO controller for a multi-modal system is to measure a combination of modes. This would be achieved, for the present test case, by placing the sensor at any location other than the elastic axis. Figure 14 shows the time evolution of the control input to achieve the load alleviation and a sample of three coefficients of the controller. It is worth mentioning that the flap rotation and angular speed meet fully physical limitations of the experimental apparatus (rotation up to ±7 • and speed up to 15 Hz). assumed that the sensor measures solely the torsional mode. An alternative approach to mitigate the effects of a SISO controller for a multi-modal system is to measure a combination of modes. This would be achieved, for the present test case, by placing the sensor at any location other than the elastic axis.  Figure 14 shows the time evolution of the control input to achieve the load alleviation and a sample of three coefficients of the controller. It is worth mentioning that the flap rotation and angular speed meet fully physical limitations of the experimental apparatus (rotation up to ±7° and speed up to 15 Hz).

Conclusions
The paper presents the implementation of an adaptive feedforward control strategy to alleviate gust-induced aeroelastic vibrations. The test case is a model problem, consisting of a typical wing section undergoing pitch and plunge motions. The aerodynamics are given by the classic theory of Theodorsen. The stiffness of the springs is modelled in a polynomial form, introducing a nonlinear element to a linear aeroelastic model. The aeroelastic parameters were set to match a wind tunnel test

Conclusions
The paper presents the implementation of an adaptive feedforward control strategy to alleviate gust-induced aeroelastic vibrations. The test case is a model problem, consisting of a typical wing section undergoing pitch and plunge motions. The aerodynamics are given by the classic theory of Theodorsen. The stiffness of the springs is modelled in a polynomial form, introducing a nonlinear element to a linear aeroelastic model. The aeroelastic parameters were set to match a wind tunnel test rig and a good comparison between predictions and measurements was found. The effects of structural nonlinearity, for this specific test case, give rise to the appearance of subcritical limit cycle oscillations. The demonstration of the control strategy is performed on a single point near but below the onset of limit cycle oscillations. The control strategy requires an approximate model of the aeroelastic plant and this was created using a system identification method retaining the nonlinear effects on the system response. Although the controller was designed to be a single-input single-output system, the response with control was found effective to reduce the peak-to-peak structural loads during the gust encounters. For a discrete "one-minus-cosine" gust, the peak-to-peak deformations were reduced by about 50% with a negligible effect on the unmeasured structural mode. For a random atmospheric turbulence encounter, which has statistical properties defined by von Kármán spectrum, the feedforward control significantly reduced the mean and amplitude of aeroelastic vibrations.
Directions for future work include: • the extension of the current work to a multi-input multi-output aeroelastic system; • an improved treatment of the nonlinearities in creating the approximate system transfer functions; • the implementation of the controller on an experimental aerial platform which will be flight tested. Initial flight tests have already been conducted to characterize the system dynamic response to forced excitations.

Acknowledgments:
Wang acknowledges the financial support received from the China Scholarship Council (CSC) to carry out the research presented in this work at the University of Southampton.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Aeroelastic parameters in dimensional and non-dimensional form and their definitions are reported below. airfoil semi-chord C L , C m lift and pitch moment coefficients C h , C α viscous damping in plunge and pitch, respectively C c h critical damping in plunge, 2 √ mK h C c α critical damping in pitch, 2 √ I α K α