Gust load alleviation: identification, control and wind tunnel testing of a 2D aeroelastic airfoil C. Poussot-Vassal, F. Demourant, A. Lepage

and Dominique Le Bihan

Abstract-One important element in the progression of aircraft environmental impact reduction is to reduce their overall weight (without impacting other consumption-oriented performance index, such as drag). In addition to the numerous work conducted in material and structural engineering, from a control viewpoint, this challenge is strongly connected to the need of the development and assessment of dedicated load control strategies in response to gust disturbances. Indeed, the load factors due to gust are considered as a sizing criteria during the aircraft conception steps and require specific verification according to the certification process. To this end, a dedicated experimental research program based on Wind Tunnel (WT) campaigns has been carried out. More specifically, the paper contributions are twofold: (i) to identify the gust load effect using two different versatile frequency-domain techniques, namely the Loewner interpolation and a modified subspace approach, and (ii) to design and implement an active closed-loop control to alleviate the gust main effect. The entire procedure is validated in a wind tunnel set-up, involving a gust generator device and a 2D aeroelastic airfoil, for varying configuration travelling from sub to transonic airflow and varying angles of attack, emphasizing the effectiveness and robustness of the overall approach.

Index Terms-Gust load control, Rational interpolation, Subspace, Robust control, Wind tunnel test.

I. INTRODUCTION

A. Context and experimental set-up

T HE topic of this paper concerns the design and the experimental validation of an active closed-loop gust load alleviation controller applied on a 2D wing in wind tunnel testing. Real experiments were achieved and validated within the ONERA Wind Tunnel (WT) facility. Real experiments were achieved and validated within the ONERA S3Ch Wind Tunnel (WT) facility for various test conditions -Angles of Attack (AoA) and Mach numbers 1 . The contributions made available in this paper encapsulate many research fields, going from mechanical design to unsteady aerodynamic characterisation of gust, including signal processing and aeroelasticity studies. However, the present paper is mainly focussed on the control engineering side, and more specifically, the open-loop dynamical model identification and closed-loop control design steps. An objective is to provided a rationalization of the gust control design.

From this control viewpoint, it is worth mentioning that the gust load alleviation problem has driven a considerable attention in the literature to substantially enhance load aircraft ONERA -The French Aerospace Lab, F-31055 Toulouse, France. Corresponding author: C. Poussot-Vassal, charles.poussot-vassal@onera.fr 1 http://www.onera.fr/en/news/gust-effect-control-a-european-first-at-onera (project webpage). performances, see e.g. [START_REF] Demourant | A frequency domain identificationcontrol approach for a flexible aircraft[END_REF], [START_REF] Moulin | Gust load alleviation using special control surfaces[END_REF], [START_REF] Torralba | Self-scheduling multiobjective control law design for a flexible aircraft[END_REF], [START_REF] Aouf | Robust gust load alleviation for a flexible aircraft[END_REF], [START_REF] Haghighat | A model predictive gust load alleviation controller for a highly flexible aircraft[END_REF], [START_REF] Haghighat | A model predictive gust load alleviation controller for a highly flexible aircraft[END_REF], [START_REF] Giesseler | Gust load alleviation based on model predictive control[END_REF], [START_REF] Demourant | A linear parameter-varying multiobjective control law design based on youla parametrization for a flexible blended wing body aircraft[END_REF], [START_REF] Xu | Aircraft design with active load alleviation and natural laminar flow[END_REF], [START_REF] Alama | Active gust load alleviation system for flexible aircraft: Mixed feedforward/feedback approach[END_REF], [START_REF] Schirrer | Feedback Control Designs in Modeling and Control for a Blended Wing Body Aircraft[END_REF] 2 . However a few involving real wind tunnel tests, and, to the authors' best of knowledge, most of them refer to tests performed in the low speed regime [START_REF] Tang | Experiments and Analysis for a gust generator in a wind tunnel[END_REF], [START_REF] Ricci | Wind tunnel testing of an active controlled wing under gust excitation[END_REF]. The originality of this work relies on the WT investigations achieved up to the transonic regime covering from the generation to the active alleviation of gust loads, which makes the results of this paper particularly relevant and the experimental demonstrations quasi-unique [START_REF] Silva | Development of aeroservoelastic analytical models and gust load alleviation control laws of a sensor craft wind-tunnel model using measured data[END_REF] 3 . In this work, the wind tunnel investigations cover the transonic air speed, which is the operation cruise regime of most modern aircraft and the most complex aerodynamic range. In aerodynamics, transonic flow involves sub and supersonic local flow in the same flow field (airflows surrounding an aircraft are locally below, at, and above the speed of sound). The transonic condition refers to compressibility effects and is characterised by discontinuities and nonlinear effects (drag divergence, shock wave, flow separation...). This constitutes the main challenge addressed in this work. As illustrated on Figure 1, the experimental research mean involved in this work consists of two main parts implemented within the WT facility. The first one is dedicated to the generation of the gust field, composed of two airfoils, and is called the gust generator. The second one is the 2D wing model which aerodynamic and aeroelastic behaviours were studied and tailored with objective of assessing the gust load alleviation functions. A schematic overview of the experimental set-up is also shown in Figure 2 4 . Fig. 2. Experimental set-up schematic view: gust wave (varying coloured surface) generated by the gust generator (left horizontal wings) attacking the controlled 2D airfoil (right horizontal wing). The airfoil is equipped with a movable surface (green) and 20 sensors (three of them, #4, #12, #15, are materialized by pink dots, and will be used either as control measurements or performance objectives). v is the wind stream vector of the wind tunnel and w is the vertical gust vector caused by the gust generator.

The gust generator (or disturbance generator). From an air flow stream generated by the WT, the gust generator aims at superposing parametric vertical gust velocities. This is achieved by the adequate actuation of the two airfoils performed by four servo-hydraulic jacks which synchronously rotate the two airfoils (i.e. pitch dynamic motions about the quarter chord point). Each airfoil is driven at its two sides and the actuation device includes a security system that stops the actuation whenever the two roots of one airfoil are not at the same angular position, thus preventing any unwanted torsional deformation of the airfoils. The two airfoils are made of an optimized architecture that combines carbon reinforced skins and metallic parts to reduce inertia and allow high frequencies of actuation under heavy aerodynamic loads. In addition to inertial and strength constraints, the dynamic behaviour of the airfoils was tuned such that the frequency bandwidth of the gust generator is about 100Hz.

The controlled 2D airfoil (controlled system). It represents the considered system to be first identified and then controlled. The model consists in a 2D airfoil mounted between the walls of the WT sections with specific boundary conditions allowing two flexible degrees of freedom: pitch and plunge/heave motions. Is is equipped with a movable full span trailing edge control surface (green area on Figure 2) which is actuated on each side by a high torque -high speed hydraulic actuator. This wing is symmetrically actuated to avoid torsional effects. The actuator's bandwidth is of 100Hz, but in the closed-loop set-up, will be limited to 50Hz to focus on the first structural mode (for this first attempt of active control demonstration). The wing is also equipped with 20 acceleration sensors located on the wingspan (some stations are materialized by pink marks on Figure 2). Figure 3 illustrates typical open-loop frequency responses (transfer gain) of the vertical acceleration z#15 at station #15 obtained in WT facility at four different configurations/combinations of Mach number and AoA, in response to both gust disturbance w(t) (left frame) and movable control surface input u(t) (right frame). The correspondence between each frequency amplification and the associated physical motions of the wing (obtained during the dynamic characterisation of the wing in laboratory), are also shown.

Gust load alleviation control objective. Accordingly to the experimental device and to standard aeronautical requirements, the main problem addressed here is, for all the considered n s Mach numbers/AoA combinations, to reduce the gust effect w(t) on the wing vertical acceleration at station #15, i.e. z#15 (t) (e.g. reduce Figure 3 left frame heave peaks) by an adequate actuation of the movable control surface input u(t) 5 . More specifically, attention will be given in disturbance attenuation over the frequency range Ω = [15 50]2π rad/s. From the control designer point of view and in a linear framework, by denoting S (i) (i = 1, . . . , n s ) as the nominal open-loop model of system for the ith Mach number/AoA configuration, with associated transfer function H (i) (s), this objective can be recast as finding K (s) such that,

K (s) := arg min K(s)∈H2 F l H (i) (s), K(s) w→z H∞ , (1)
for all i = 1, . . . , n s , where K (s) is the linear time invariant optimal controller that minimizes the H ∞ -norm 6 of the lower fractional interconnection (F l H (i) (s), K (s) ) between the system H(s) and aforementioned active controller. Problem (1) means that one is looking for a strictly proper and stable controller K (s) that (stabilizes and) attenuates the first peak, where gust generally occurs. Therefore, throughout the paper, focus on vertical accelerations reduction will be addressed. In the control literature, this objective can also be viewed as an active damping problem. Remark 1 (WT experiments): It is worth mentioning that, before obtaining the open-loop data, more than hundreds of experiments have been performed during two months to calibrate and adjust the gust generator, wing profile and perform adequate measurements. Then, as rooted on the numericalbased open-loop models and robust control controller synthesis (described in this paper), the WT closed-loop experiments have been performed in two days. This consideration is very important since WT experiments are relatively expensive.

B. Contributions

The main contributions of this work are twofold: first (i) in successfully identifying the dynamical model through two advanced frequency-domain methods: one grounded on the Loewner interpolation framework (see [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF], [START_REF] Ionita | Lagrange rational interpolation and its applications to approximation of large-scale dynamical systems[END_REF], [START_REF] Antoulas | Model reduction and approximation theory and algorithms[END_REF]) and the other, based on an innovative subspace formulation including eigenvalues constraints based on LMI regions 7 . Second, (ii) 5 Indeed, in the context of this study, this attenuation will lead to wing root bending moment minimization which is immediately linked to the mass at the wing/fuselage junction. In a more general framework of the actual civilian aircraft research, one objective is to the increase of the laminar flow working zone which are central objectives in consumption optimization [START_REF] Burlion | Longitudinal Manoeuvre Load Control of a Flexible Large-Scale Aircraft[END_REF]. 6 The H∞-norm, with appropriated weighting functions can also minimizes the root mean square over the frequency range Ω (see Section III). 7 It is shown that the two approaches are complementary since one is involved for the synthesis, while the other for numerical validation phases. in synthesizing and implementing a gust load alleviation controller using a structured H ∞ optimization criteria approach (see [START_REF] Apkarian | Nonsmooth H∞ Synthesis[END_REF]). As a glimpse of the obtained results, Figure 4 illustrates the vertical acceleration attenuation obtained in WT, for constant AoA and for various Mach numbers, in response to varying frequency gust wave disturbances. Figure 4 emphasizes the benefit of the proposed identification and control procedure, leading to the first frequency peak attenuation, while keeping the rest of the frequency response unchanged (see details in the rest of the paper).

Remark 2 (Media material, at Mach number 0.73 and AoA 0 degree): in addition to the present paper, a video is also made available at the following address: http://www.onera.fr/en/news/gust-effect-controla-european-first-at-onera, illustrating the attenuation of the displacement amplifications in response to a gust disturbance at the characteristic frequency 25Hz (heave frequency), with and without control law, in a transonic configuration.

C. Paper structure

The rest of the paper is dedicated to the description and illustration of the applied methodology to achieve results of Figure 4. The paper is organized in a top down framework. First, in Section II, two frequency-domain methods are presented and applied to the considered system to identify the linear dynamical models that accurately captures the dynamics in response to gust disturbance and surface deflection. Then, Section III is devoted to the design of the dedicated robust active gust load controller, followed by some simulations and experimental results (illustrating both the consistency between numerical and experimental data, and effectiveness of the methodology). Conclusions are finally given in Section IV.

D. Mathematical notations

Let us consider a stable Multiple Input/Multiple Output (MIMO) linear dynamical system S, with n u (resp. n y ) ∈ N * input(s) (resp. output(s)), represented by its n-th order transfer function H ∈ C ny×nu , where H(s) = C(sE -A) -1 B + D. Its associated realisation reads S : (E, A, B, C, D),

S : E ẋ(t) = Ax(t) + Bu(t) , y(t) = Cx(t) + Du(t). (2)
The input vector used here gathers the gust disturbance and control deflection as u(t) = [w(t) u(t)]. The considered output vector gathers three acceleration sensors at stations #4, #12 and #15, i.e. y(t) = [z #4 (t) z#12 (t) z#15 (t)]. Then, the approximate model, or reduced order model of order r n, will be denoted as Ĥ(s) = Ĉ(s Ê-Â) -1 B+ D. Its associated realisation reads Ŝ : ( Ê, Â, B, Ĉ, D),

Ŝ : Ê ẋ(t) = Âx(t) + Bu(t) , ŷ(t) = Ĉx(t) + Du(t). (3)

II. OPEN-LOOP MODEL IDENTIFICATION AND VALIDATION

PHASE

As illustrated on Figures 3 and4, the ONERA experimental WT facility allows the generation of accurate gust disturbances at varying frequencies (and amplitudes). Consequently, it is possible to obtain the frequency-domain responses Φ i ∈ C ny×nu for different frequency samples ω i (i = 1, . . . , N ). Then, one can write:

H(ıω i ) = Φ i , ( 4
)
where H is the exact 8 transfer function of the system. Based on the couple (ω i , Φ i ) N i=1 , the main purpose of this section is to obtain a low-order rational LTI model of the form,

Ĥ(s) = Ĉ(s Ê -Â) -1 B + D ∈ H ny×nu ∞
, with realisation Ŝ : ( Ê, Â, B, Ĉ, D), defined as: [START_REF] Haghighat | A model predictive gust load alleviation controller for a highly flexible aircraft[END_REF] that well matches the frequency samples (ω i , Φ i ) N i=1 and hopefully reproduces the real transfer H (in other words, ŷ should accurately reproduce y). Let x(t) ∈ R r , u(t) ∈ R nu and ŷ(t) ∈ R ny be the states, inputs and outputs vectors, respectively. To this aim, two approaches have been experimented: the first one is based on the Loewner framework (Section II-A), and the second one, more innovative, on the subspace one with LMI constraints (Section II-B). Both approaches are used in this work for different purpose. While the second approach will be used for the control design step, the former one is involved for numerical validation purpose 9 .

Ŝ : Ê ẋ(t) = Âx(t) + Bu(t) , ŷ(t) = Ĉx(t) + Du(t),

A. First approach: frequency-domain matrix function interpolation and model approximation

The first approach is based on the interpolatory framework well defined in [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF], [START_REF] Antoulas | Model reduction and approximation theory and algorithms[END_REF], involving the Loewner matrices, and recalled hereafter. It is composed of two steps: first an exact rational model interpolation, followed by a reduction procedure 10 . To this aim, let us first partition the collected 8 By exact, one intends the unknown, real exact function that describes the frequency responses obtained in WT facility. Of course, any model reproducing the data collected in experiment will be called exact if perfect match is ensured. However, every model can still be amended or questioned, and the search for the perfect model is hopelessly unachievable (and is not necessary for control design). 9 For a more detailed qualitative comparison between the Loewner and subspace approach, reader is invited to refer to [START_REF] Demourant | A new frequency-domain subspace algorithm with restricted poles location through LMI regions and its application to a wind tunnel test[END_REF]. 10 Most of the results here described are available in [START_REF] Ionita | Lagrange rational interpolation and its applications to approximation of large-scale dynamical systems[END_REF]. A very complete survey is also available in [START_REF] Antoulas | Model reduction and approximation theory and algorithms[END_REF].

data (ω i , Φ i ) N i=1 in two disjoint sets as follows (N = q + k): ı[ω 1 , . . . , ω N ] = [µ 1 , . . . , µ q ] ∪ [λ 1 , . . . , λ k ] [Φ 1 , . . . , Φ N ] = [ṽ 1 , . . . , ṽq ] ∪ [ w1 , . . . , wk ]. (6)
Then, define l j ∈ C ny×1 (j = 1, . . . , q) and r i ∈ C nu×1 (i = 1, . . . , k) the q left and k right tangential directions. Using these tangential directions, let us define v * j = l * j ṽj ∈ C 1×nu and w i = wi r i ∈ C ny×1 the left and right tangential data directions, respectively. Then, considering the right or column data and the left or row data, and assuming that λ i and µ j , are distinct, the right interpolation data are organized as:

Λ = diag(λ 1 , . . . , λ k ) ∈ C k×k R = [r 1 . . . r k ] ∈ C nu×k W = [w 1 . . . w k ] ∈ C ny×k    , (7)
and the left interpolation data are given as,

M = diag(µ 1 , . . . , µ q ) ∈ C q×q L * = [l 1 . . . l q ] ∈ C ny×q V * = [v 1 . . . v q ] ∈ C nu×q    . (8)
Based on ( 7) and ( 8), the corresponding Loewner L ∈ C q×k and shifted Loewner L σ ∈ C q×k matrices, are constructed as follows, for i = 1, . . . , k and j = 1, . . . , q:

[L] j,i = v * j r i -l * j w i µ j -λ i = l * j H(µ j ) -H(λ i ) r i µ j -λ i [L σ ] j,i = µ j v * j r i -λ i l * j w i µ j -λ i = l * j µ j H(µ j ) -λ i H(λ i ) r i µ j -λ i , (9)
which satisfy the following Sylvester equations:

LΛ -ML = LW -VR, L σ Λ -ML σ = LWΛ -MVR. ( 10
)
Exact rational interpolation. Based on the above data and tangential directions, the interpolation problem can be stated as follows:

Problem 1 (Data-driven interpolation [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF]): Given left interpolation driving frequencies {µ i } q i=1 ∈ C with left output or tangential directions {l i } q i=1 ∈ C ny , producing the left responses {v i } q i=1 and right interpolation driving frequen-

cies {λ i } k i=1 ∈ C with right input or tangential directions {r i } k i=1 ∈ C nu , producing the right responses {w i } k i=1
, find a (low order) system Ŝ that such that the resulting transfer function Ĥ(s) is an (approximate) tangential interpolant of the data, i.e. satisfies the following left and right interpolation conditions:

l * j Ĥ(µ j ) = v * j for j = 1, . . . , q and Ĥ(λ i )r i = w i for i = 1, . . . , k . (11)
Problem 1 can be solved thanks to the following theorem, proposed by [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF].

Theorem 1 (Loewner framework [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF]): Given right and left interpolation data as in ( 7)-( 8), and assume that k = q and let (L σ , L) be a regular pencil where λ i or µ j are (distinct and) not eigenvalue. The rational transfer function H(s), with realisation S : (E, A, B, C, 0) constructed as

E = -L, A = -L σ , B = V and C = W, (12)
with transfer minimal descriptor function

H(s) = W(L σ -sL) -1 V, (13)
interpolates the left and right constraints [START_REF] Tang | Experiments and Analysis for a gust generator in a wind tunnel[END_REF]. Theorem 1 allows to obtain a model S : (E, A, B, C, 0), whose transfer function interpolates the left and right constraints as stated in Problem 1. Moreover, an important property of the Loewner matrix is that,

rank(L) = n ≤ N, (14)
encodes the McMillian degree n of the rational function interpolation. Therefore, the exact Loewner interpolation provides a realisation of a rational transfer function that exactly interpolates all the data (ω i , Φ i ) N i=1 , with the minimal order realisation thanks to the Loewner matrix rank.

Approximate rational interpolation. Once an exact model S : (E, A, B, C, 0) with transfer function H(s) with McMillian degree of n ≤ N , and which tangentially interpolates the data is available, to obtain a reduced order model Ŝ of order r ≤ n that well approximates S one simply apply an SVD as follows:

L = Y 1 Y 2 Σ 1 Σ 2 X * 1 X * 2 , (15)
where

Σ 1 = diag(σ 1 , . . . , σ r ) ∈ R r×r , Σ 2 = diag(σ r+1 , . . . , σ n ) ∈ R (n-r)×(n-r) (where σ 1 ≥ σ 2 ≥ • • • ≥ σ n ) and Y 1 , Y 2 , X 1 , X 2 are of appropriate dimensions.
Then the reduced order model Ĥ is simply obtained by the Petrov-Galerkin projection 11 :

Ŝ : ( Ê, Â, B, Ĉ, 0) = (-Y * 1 LX 1 , -Y * 1 L σ X 1 , Y * 1 V, WX 1 , 0). (16)
In the objective -here considered -of identifying a model that accurately reproduces the gust and actuator effects on the wing profile, the main interest in Theorem 1 is that it enables to both find the exact rational function H (and realisation S) that matches all the frequency-domain data (ω i , Φ i ) N i=1 and, allows, at a very low additional cost, to obtain a good approximated model Ĥ (and realisation Ŝ) 12 .

Remark 3 (About optimality conditions): The truncation performed when applying the SVD in [START_REF] Burlion | Longitudinal Manoeuvre Load Control of a Flexible Large-Scale Aircraft[END_REF] does not guarantee the so-called H 2 -optimality conditions of the projected model (see e.g. [START_REF] Gugercin | H 2 Model Reduction for Large Scale Linear Dynamical Systems[END_REF]). These conditions can be achieved through an adequate interpolation data (µ i and λ i ) selection, as explained in recent papers [START_REF] Drmac | Vector fitting for matrix-valued rational approximation[END_REF], [START_REF]Quadrature-Based Vector Fitting For Discretized H 2 Approximation[END_REF], but not tested here. Nevertheless, as illustrated in the forthcoming applicative section, very accurate results are obtained.

Application to the WT data. For the considered application, Theorem 1 is then used with N = 292 sampled data points (k = q = 146), n y = 3 outputs (acceleration at stations #4, #12 and #15) and n u = 2 inputs (wind gust and actuator) and r = 20. First, exact interpolation is achieved with a state-space model of order n = 292. The vertical acceleration response corresponding to station #15 is illustrated on Figure 5, with red lines, providing a model that perfectly matches the measurements (blue points). Then, thanks to the Loewner matrix properties, a rank truncation is performed to reduce the model to an order r = 20 (singular values decay of L is not provided for space limitations), illustrating the good matching (dash-dotted pink lines).

Figure 5 illustrates the very good matching of the proposed models for two distinct configurations (sub and transonic). However, at this stage, no property on the spectrum associated to the obtained exact and approximate models can be given (e.g. stability of the model cannot be guaranteed). Indeed, the models might result in unstable ones in both cases (exact and approximate), which can be an issue since, in practice, the system is stable. This limitation is tackled in the next sub-section though an approach that also uses the frequency domain data, but in a different framework.

B. Second approach: frequency-domain subspace with LMI constraints

Similarly to the Loewner framework, the subspace approaches are well appropriate to identify and approximate frequency-domain sampled data, by a state-space realisation as in (5) (here with E = I n and D non necessarily null). However, as in the Loewner case, in its basic version, no constraint on the obtained model eigenvalues location can be enforced (e.g. stability constraint). The main result in this section is to propose a solution to constraint the poles of the identified model in a user-defined region using Linear Matrix Inequalities (LMI). As made clearer in what follows, these LMI regions can simply be stability domain or more complex regions (horizontal strip, ellipse etc.). To the authors' best of knowledge, in the subspace context, the stability constraint is treated in [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF], [START_REF] Maciejowski | Guaranteed stability with subspace methods[END_REF], [START_REF] Chui | Realization of stable models with subspace methods[END_REF], [START_REF] Gestel | Identification of stable models in subspace identification by using regularization[END_REF], [START_REF] Lacy | Subspace identification with guaranteed stability using constrained optimization[END_REF]. In [START_REF] Lacy | Subspace identification with guaranteed stability using constrained optimization[END_REF] a specific formulation based on discrete Lyapunov function is involved to impose stability whereas in the approach presented here a more general formulation is proposed to take into account LMI regions poles location and not only the stability one. The approaches developed in [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF], [START_REF] Maciejowski | Guaranteed stability with subspace methods[END_REF] can distort the solution and techniques proposed in [START_REF] Chui | Realization of stable models with subspace methods[END_REF], [START_REF] Gestel | Identification of stable models in subspace identification by using regularization[END_REF] are not direct and based on iterative algorithm. Here the proposed solution is direct without any tuning parameters for a large class of poles location, where the stability preservation is a particular configuration. That is why we will talk about LMI regions rather than stability domain. This point is illustrated in the application section by using an elliptic constraint to impose not only stability for discrete models but also to skip irrelevant and unrealistic under-damped poles.

Preliminary results on subspace methods. Fundamental subspace-based results are briefly recalled here. However additional information can be found in [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF], [START_REF] Liu | Frequency domain structural system identification by observability range space extraction[END_REF] 13 . Let

(ω i , Φ i ) N i=1 .
us first consider the frequency domain discrete state-space representation:

e ıω X(ω) = AX(ω) + BU(ω) Y(ω) = CX(ω) + DU(ω) , (17)
where X(ω) ∈ C r , U(ω) ∈ C nu and Y(ω) ∈ C ny are the Fourier transform of the state x(t), the input u(t) and the output y(t) vectors, respectively. Let us then denote X m (ω i ) the state at frequency ω i for an input vector U(ω i ) with one in the m-th row and zero elsewhere (i = 1, . . . , N and m = 1, . . . , n u ), then [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF] similarly reads as

e ıω X(ω) = AX(ω) + BU(ω) G(ω) = CX c (ω) + D , (18)
where X c (ω) = [X 1 (ω), . . . , X nu (ω)] ∈ C r×nu . Moreover, if a discrete frequency domain data set (ω i , Φ i ) N i=1 is considered then, one has G(ω i ) = Φ i and the following relation holds [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF] (with r + q ≤ N and q ≥ r):

G = OX c + ΓW (19)
where,

G =      Φ1 . . . Φ N e ıω 1 Φ1 . . . e ıω N Φ N . . . . . . . . . e ı(q-1)ω 1 Φ1 . . . e ı(q-1)ω N Φ N      ∈ C ny q×nuN , X c = X c (ıω 1 ) . . . X c (ıω N ) ∈ C r×nuN and O =      C CA . . . CA q-1      ∈ R nyq×r (20)
has full rank r for all values q ≥ r and finally,

W =      In u . . . In u e ıω 1 In u . . . e ıω N In u . . . . . . . . . e ı(q-1)ω 1 In u . . . e ı(q-1)ω N In u      ∈ C nuq×nuN , Γ =      D 0 . . . 0 CB D . . . 0 . . . . . . . . . . . . CA q-2 B CA q-3 B . . . D      ∈ R ny q×nuq .
Relation [START_REF] Antoulas | Model reduction and approximation theory and algorithms[END_REF] emphasizes that the extended observability matrix O (20), which depends on A and C, can be written as a combination of inputs/outputs data G, X c , and W, only. Then, an elegant way to extract A and C matrices from O, is to first apply the orthogonal projection W ⊥ defined as

W ⊥ = I -W T (WW T ) -1 W, (21)
to cancel the ΓW term 14 . Indeed, by right multiplying [START_REF] Antoulas | Model reduction and approximation theory and algorithms[END_REF] with W ⊥ as defined in [START_REF] Demourant | A new frequency-domain subspace algorithm with restricted poles location through LMI regions and its application to a wind tunnel test[END_REF], the following is obtained:

GW ⊥ = OX c W ⊥ .

Then, by denoting

W G = R 11 0 R 21 R 22 Q T 1 Q T 2 = RQ, (22)
an effective way to extract A and C is to use a QR and an SVD matrix decomposition of GW ⊥ ∈ R nyq×nuq , and by noticing that:

GW ⊥ = R 22 Q T 2 = Û Σ VT Q T 2 = [ Ûs Ûo ] Σs 0 0 Σo VT s VT o Q T 2 .
where Σs = diag(σ 1 , . . . , σ r ) ∈ R r×r and Σo = diag(σ r+1 , . . . , σ nyq , 0) ∈ R (nyq-r)×(nyq-r) (where It is shown that if q ≥ r, then Ûs = Ô = OS (S ∈ R r×r ) where Ô is the extended observability matrix of the estimated r-th order state-space representation (I r , Â, B, Ĉ, D) is equivalent to (I n , A, B, C, D) with  = S -1 AS and Ĉ = CS. It is well-known that if the frequency-domain data are generated by an LTI system of order r, a clear drop appears at the (r + 1) th singular value of R 22 [START_REF] Gugercin | H 2 Model Reduction for Large Scale Linear Dynamical Systems[END_REF], then the choice of the order is obvious. If data are obtained from experimental tests, this cut off is not necessarily clear. However this decomposition can be seen as tuning parameter to obtain a satisfactory complexity / accuracy trade-off. In brief, as in the Loewner framework, thanks to the SVD decomposition the main dynamic behaviour of the system can be extracted using the singular value drop. Consequently the system order r can be determined by the highest singular values of the R 22 block represented by Σs . Thus,

σ 1 ≥ σ 2 ≥ • • • ≥ σ n ) represent
 = ϕ † b and Ĉ = J 3 Ûs , (23)
where ϕ = J 1 Ûs , b = J 2 Ûs , M † = (M T M ) -1 M T and J 1 = [I (q-1)ny 0 (q-1)ny×ny ] J 2 = [0 (q-1)ny×ny I (q-1)ny ]

J 3 = [I ny 0 ny×(q-1)ny ]
and where

J 1 Ûs =      Ĉ Ĉ Â . . . Ĉ Âq-2      and J 2 Ûs =      Ĉ Â Ĉ Â2 . . . Ĉ Âq-1      . (24)
Let us assume that  and Ĉ are now known, the transfer function became affine in B and D and can consequently be obtained by solving the following standard convex optimization problem:

{ B, D} = arg min B ∈ R n×nu D ∈ R ny×nu N i=1 Φ i -Ĥ(ω i , B, D) 2 F .
where Ĥ(ω i , B, D) = Ĉ(e ıωi I r -Â) -1 B + D. In addition it is shown that in the noise-free case, ĤN -H H∞ = 0, if H is a stable system of order r and N ≥ r + q and q ≥ r are satisfied. Besides, this approach is consistent in the noisy case, i.e. lim N →∞ | ĤN -H H∞ = 0 with q ≥ r, if H is a stable system of order r and the noise affecting the data is a zero mean complex random variable and has bounded fourth-order moments (details are given in [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF]). Preliminary results on LMI regions for LTI systems. Before exposing the proposed algorithm, let us now recall some preliminary results on LMI regions which has initially been used in the context of controller synthesis [START_REF] Chilali | H∞ design with pole placement constraints: an LMI approach[END_REF], [START_REF] Chilali | Robust pole placement in LMI regions[END_REF].

Definition 1 (LMI region): A subset D of the complex plane is called LMI region if there exists a symmetric matrix P = P T ∈ R n×n and a matrix Q ∈ R n×n such that

D = {z ∈ C|f D (z) = P + Qz + Q T z < 0}. ( 25
)
Let us notice that the characteristic function f D (z) = P+Qz+Q T z is a n×n Hermitian matrix. Moreover, the LMI regions are symmetric with respect to the real axis 15 . Without lack of generalities, inequalities ( 26)-( 30) provide a set of characteristic functions provided as LMI regions, namely left half complex plane (26), Ellipse of centre c 0 with half-length axis a and b [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF], Vertical strip between α 1 and α 2 [START_REF] Maciejowski | Guaranteed stability with subspace methods[END_REF], Horizontal strip between ±β (29), Conic sector of tip c 0 and angle θ [START_REF] Gestel | Identification of stable models in subspace identification by using regularization[END_REF].

Remark 4 (About continuous and discrete stability regions): The complex domain setting holds for both continuous and discrete systems. In the present context the discrete-domain is used, therefore the stability constraint is obtained by imposing |z| < 1, i.e. all eigenvalues within the unit circle (by applying a = b = 1 and c 0 = 0 to the ellipse characteristic function [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF]). The equivalent constraint for continuous systems is then obtained by applying the standard transformation z = e sTe , where s is the Laplace variable and T e , the sampling period.

Remark 5 (Convex set intersection): The convexity property in preserved by the intersection operation on convex sets. Consequently the final subset is convex too, and, given two LMI regions D 1 and D 2 described by the characteristic function f D1 and f D2 respectively, the characteristic function

f D of D = D 1 ∩ D 2 is simply recast as f D = diag(f D1 , f D2 ).
Consequently, the formulation of poles placement objectives within LMI regions, is given by the following Theorem 2.

Theorem 2 (LMI regions [START_REF] Chilali | H∞ design with pole placement constraints: an LMI approach[END_REF], [START_REF] Chilali | Robust pole placement in LMI regions[END_REF]): Let us consider a dynamical system as in [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF], P = P T ∈ R n×n , Q ∈ R n×n , then A has all its eigenvalues in the LMI regions D defined by the characteristic function f D (z):

D = {z ∈ C|f D (z) = P + Qz + Q T z < 0}, (31)
if and only if there exists a symmetric Ψ = Ψ T > 0 such that:

P ⊗ Ψ + Q ⊗ AΨ + Q T ⊗ (AΨ) T < 0. ( 32
)
Subspace approach under LMI regions constraints. As rooted on the above preliminary results, we are now ready to state the main result of this section, namely, the subspace approach with LMI constraints, allowing to impose the identified model poles placement.

Proposition 1 (Subspace with LMI regions constraints): Given P = P T ∈ R r×r and Q ∈ R r×r characterizing the LMI region D defined as

D = {z ∈ C|f D (z) = P + Qz + Q T z < 0}. ( 33
)
The eigenvalues of  = ϕ † b ∈ R r×r lies in D if and only if the following problem: min Ã, Ψ, β subject to (35), (36), (37)

β ( 34
)
P ⊗ Ψ + Q ⊗ Ã + Q T ⊗ ÃT < 0 (35) I r (ϕ Ã -bΨ) (ϕ Ã -bΨ) T βI r > 0 (36) Ψ > 0 (37) Re(z) < 0 ⇔ z + z < 0 (26) Re(z) -c 0 a 2 + Im(z) b 2 < 1 ⇔     -1 -c 0 /a + z(1/2a -1/2b)+ z(1/2a -1/2b) -c 0 /a + z(1/2a + 1/2b)+ z(1/2a -1/2b) -r     < 0 (27) α 1 < Re(z) < α 2 ⇔ z + z -2α 2 0 0 z + z + 2α 1 < 0 (28) |Im(z)| < β ⇔ -2β z -z z -z -2β < 0 (29) Re(z) -c 0 tan θ < -|Im(z)| ⇔ (-2c 0 + z + z) sin θ (z -z) cos θ (z -z) cos θ (-2c 0 + z + z) sin θ < 0 ( 30
)
where à = ÂΨ ∈ R r×r and Ψ = Ψ T ∈ R r×r and β ∈ R, has a solution.

Proof: Thanks to Theorem 2, eigenvalues of  are located in the LMI region D if Ψ = Ψ T > 0 exists such that:

P ⊗ Ψ + Q ⊗ ÂΨ + Q T ⊗ ( ÂΨ) T < 0 ⇔ P ⊗ Ψ + Q ⊗ à + Q T ⊗ ÃT < 0, where à = ÂΨ. Besides, let J 1 Ûs  = J 2 Ûs ⇔ ϕ  = b with ϕ = J 1 Ûs and b = J 2 Ûs . Since ϕ  = b then ϕ ÂΨ = bΨ ⇔ ϕ à = bΨ with à = ÃΨ. Finding à is equivalent to solve min ϕ à -bΨ 2 2 .
Consequently, one can write:

min ϕ Ã -bΨ 2 2 ⇔ min β∈R β > 0 (ϕ Ã -bΨ) T (ϕ Ã -bΨ) < βI r .
Then, a quadratic cost function has been transformed into a linear cost function under quadratic constraints. Consequently, by invoking the Schur lemma, one obtains:

(ϕ Ã -bΨ) T (ϕ Ã -bΨ) < βI r ⇔ I r (ϕ Ã -bΨ) (ϕ Ã -bΨ) T βI r > 0
where matrix inequalities clearly are affine in decision variables β, Ã and Ψ.

A numerical procedure to obtain a dynamical model satisfying the eigenvalue oriented LMI constraint is sketched in what follows (see Algorithm 1).

Remark 6 (Handling noisy data): In the noisy data context, to preserve consistency, a weighting matrix N is used based on the noise relative variance. This matrix must be a priori known. The evaluation is strictly equivalent and can be included in our optimization problem as follows:

 = (J 1 N Û1 ) † J 2 N Û1 ⇔  = ϕ † b Ĉ = J 3 N Û1
where ϕ = J 1 N Û1 and b = J 2 N Û1 .

Remark 7 (LMI constraints): If no constraints is considered, the LMI least-square solution is exactly equivalent to the baseline subspace approach. Moreover, one should note that in practice, the consistency of the problem is not affected by imposing poles placement regions if these constraints are relevant with the system to identify [START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF]. Application to the WT data. For the considered application, Algorithm 1 has been involved with N = 292 sampled data points, n y = 3 outputs (acceleration at stations #4, #12 and #15) and n u = 2 inputs (wind gust and actuator) and a sample frequency f s = 4096Hz 16 . When applying Algorithm 1, the considered LMI region is an ellipse defined by the constraint ( 27) with a = 1, b = 0.92 and c 0 = 0. This ellipse leads to impose the stability since this LMI region is included in the unitary disk. Furthermore, the half-length axis b = 0.92 avoids to under-estimate the modes damping and to delete possible weakly damped poles which can appear and which are not physical.

First, the vertical acceleration output corresponding to station #15 is illustrated on Figure 6, with r = 20 (on the same configurations as the one shown on Figure 5 for the Loewner approach). Comparing with Figure 5, this figure illustrates the very good performance of the subspace approach. Then, Figure 7 compares the maximal singular values of the MIMO models obtained with either the exact and approximate Loewner approach and the subspace with LMI constraint one, for r = 16 (order of the model used in the control design Section III). It emphasizes the fact that the subspace with LMI constraint obtained model (black solid thin) well matches the measurement points (blue points).

In addition, Figures 8 and9 present the maximal singular value responses and eigenvalues of both the baseline subspace approach and the LMI-based subspace one. On Figure 8, the identified models for both approaches are illustrated. It is notable that both methods lead to a good approximation of the frequency data. More interestingly, Figure 9 shows the corresponding eigenvalues in the discrete domain (left) and the equivalent in the continuous one (right). First, one should note that the subspace approach without LMI constraint might lead to models with unstable modes (here two unstable modes). Besides, the proposed approach leads to models where the modes are stable and are gathered in the defined LMI region constraints (materialized by the solid blue line).

III. ACTIVE CLOSED-LOOP LOAD ALLEVIATION CONTROL DESIGN AND DEMONSTRATION

A. Control set-up description

Once, at each AoA and Mach number configurations, a model Ĥ(i) (i = 1, . . . , n s , where n s is the number of models/configurations) have been identified, the second contribution of this work is to design an active closed-loop feedback controller aiming at attenuating the first amplification peak (heave structural motion illustrated on Figure 1) through the actuation of the movable surface. Before describing the load controller design approach, two main remarks should be kept in mind for the design: (i) first, the controller should be easy to implement on real-time device17 (Figure 10) (ii) and second, it should be robust to all Mach number and AoA variations 18 . With reference to Figure 2, the control signal u(t) is the control surface deflection and the measurement vector y(t) = [z #4 (t) z#15 (t)] used for the control is composed of the acceleration sensors at station #4 and #15.

Remark 8 (Considerations about the selected models): As mentioned in the introduction, for the control design step, the models obtained using the subspace with LMI constraint approach have been selected (with dimension r = 16). The models of high order (n = 292) obtained by the exact Loewner approach has been used for validation purpose. Indeed, while the latter catches all the dynamics of the frequency response (but may result in instability), which are well appropriated for frequency domain validation, the former are more tailored for controller design since they preserve model stability and constrain badly damped modes 19 . From a practical point of view it is interesting to notice that both models are used for different purposes 20 .

B. Robust controller design

According to the above preliminary objective, a controller is synthesized using a robust approach in order to obtain K , 19 Indeed, avoiding artificial unstable modes alleviates the control optimization procedure to stabilize irrelevant and unrealistic modes. 20 The subspace approach, using LMI is quite sensitive to the dimension of the seek model and models should be kept quite small (i.e. in practice, r < 20) but can guarantee stability. On the other side, the Loewner models are not sensitive to the approximation order r 100 but might result in unstable models. Since it is not the scope of the paper, no further analysis is given here.

an LTI controller that is both robust to the Mach number and AoA variations. To address this objective, the synthesis problem is recast as a multi-channel H ∞ -norm minimization one 21 . Note that the robust consideration justifies the use of an H ∞ approach instead of the H 2 one, which might appear more appropriate in view of the root mean square minimization. Mathematically, the objective is to find the optimal controller K such that,

K := arg min dim(K) = nc K ∈ K ⊆ H 1×2 2     T (1)
w→z (K) . . .

T (ns) w→z (K)     KW k H∞ (38)
where K describes the considered controller family set (see hereafter for more details), n c ∈ N * is the controller order, and T

(i) w→z (K) (for i = 1, . . . , n s ) represents the K dependent performance transfer function from w to z (see hereafter and Figure 11). Finally, W k is an additional weighting function used (i) first to enforce transfer KW k to be stable, and thus K to be stable (W k being stable), and (ii) second, to shape K roll-off behaviour. More specifically, we have

W k = G k s 2 /w 2 ro s 2 /(10 3 w ro ) 2 + s1.4/(10 3 w ro ) + 1 (39)
where G k = 10 -2 and w ro = 2π30 are the filter low frequency gain and roll-off cutting frequency, respectively. 1) Considerations about K: With reference to (38), K is the set of strictly stable rational transfer functions of order n c . Indeed K ∈ K ⊆ H 1×2 2 is a structure-defined set of rational transfer functions in C 1×2 (2 measurements, 1 control signal) whose structure is defined as

K(s) = C K (sI nc -A K ) -1 B K , (40)
where

A K ∈ R nc×nc , B K ∈ R nc×2 and C K ∈ R nc×1
are tuning parameters to be determined. More specifically, the optimal controller K ∈ K will have the structure presented in (40). Obviously, the main interest in structuring the control law as in (40) is to find a law that is both simple and performs appropriately to attenuate the targeted peak. Moreover, thanks to constraint (39), the obtained controller should have a roll-off behaviour after w ro . Designing such a controller can be easily achieved thanks to recent advances in dedicated numerical tools. In the considered case, the approach developed in [START_REF] Apkarian | Nonsmooth H∞ Synthesis[END_REF] has been invoked.

2) Considerations about T (i) w→z (K): With reference to (38), the transfer functions depending on K, describing the the performance channels are defined as (i = 1, . . . , n s ):

T (i) w→z (K) = W I F l ( Ĥ(i) , K)W O (41)
where,

W O (s) = W p (s) 0 0 W u (s) and W I (s) = 1 (42)
and where (see also Figure 11),

• The performance channel is characterized by W p (s), defined as

W p (s) = G ∞ (s/10w p + 1) s/w p + 1 I 3 , (43)
where G ∞ is the frequency-limited H ∞ -norm (evaluated over Ω) of the performance transfer from w to z, and w p = 2π30rad/s. This performance channel is applied on all the considered output, i.e. the accelerations at stations #4, #12 and #15, with objective to attenuate the considered dimensioning heave peak. • Then, the control input signal attenuation level, aiming at preserving eventual control saturation is given as

W u (s) = s/w act + 1 s/10w act + 1 , (44)
where w act = 2π30rad/s. With reference to Figure 11, Ĥ(i) denotes the identified model obtained using the subspace approach described in Section II-B with order r = 16 (see Figures 7 and8), describing the system behaviour for varying AoA and Mach number. Figure 11 illustrates the interconnection for the ith transfer. For sake of completeness, u(t) and y(t) signals denote the control surface deflection signal u(t) and the vertical accelerations at stations #4 and #15, respectively. Then, w(t) and z(t) are the gust disturbance input w(t) and vertical accelerations output at stations #4, #12 and #15 (to be attenuated), respectively. Finally, w(t) and z(t) are the filtered exogenous input and output, respectively, defining the generalized plant and input/output transfer functions to be minimized by the H ∞ -norm objective.

Ĥ(i) W O W I z(t) w(t) z(t) w(t) K y(t) u(t)

C. Numerical and experimental validation

In practice, solving (38) with n s = 5 and a LTI controller K of order n c = 8 has been obtained using the approach described in [START_REF] Apkarian | Nonsmooth H∞ Synthesis[END_REF] 22 . The obtained controller presents a frequency behaviour consistent with the constraints given by (38), namely, K ∈ K, with a low gain action in low frequencies (indeed, only the peak has to be attenuated) and roll-off in higher ones (satisfying W k as stated in (39)) 23 . The achieved attenuation level, solution of the optimization problem (38), leads to an optimal H ∞ -norm value of 0.43473(< 1), meaning that the weighting function constraints are attenuated (especially the performance criteria W p one).

Before moving to the implementation step, in order to assess the obtained controller performance, simulation are first performed by applying the optimal K on the complete model obtained with the exact Loewner approach (Section II-A) 24 .

The left frame of Figure 12 both provides open-loop and closed-loop frequency responses and performance attenuations for varying configurations from sub to transonic (varying AoA and Mach number configurations).

Once the attenuation level is numerically validated, similar experiments are performed in the WT facility and results are reported on the right frame of Figure 12. By comparing the RMS gain displayed on both (left and right) frames of Figure 12, one can observe that they are quite close to each other. The only difference between both results is concerned with the notch-like behaviour observed on Figure 12 around 40Hz at configurations of AoA of 2 degrees, which has been observed in late experiments. This discrepancy can be explained by the friction effects that slightly modified the system behaviour due to intensive use of the experimental set-up 25 . Nevertheless, this varying behaviour illustrates that the implemented control law is also robust with respect to the plant dynamical variations.

Remark 9 (Implementation specifications & tested configurations): The synthesized controllers were implemented on a real time device which was comprised of several processors and input/output boards interlinked for fast internal communication and data exchange. The I/O interface was composed of a maximum of 15 analog inputs and 18 analog outputs. A dedicated computer was used for creating, compiling and implementing the control laws in the processor boards and a real-time man/machine interface was developed to monitor the signals and change control/command parameters. Synthesized in the continuous time-domain, the controllers were converted into discrete state-space models to be directly implemented. With respect to the bandwidth of the electrohydraulic actuator, a 2kHz sampling frequency was selected. Regarding the considered application within the WT, all the configurations used in the robust design approach have been tested, while in between configurations have not been tested for time limitations. However, extension of this work might include other tests such as a more complete design, involving a parametric / uncertain model instead of local models, as stated in (38).

In addition to the frequency-domain data collected in WT and compared with simulation (Figure 12), a time-domain gust response have also been performed in simulation only 26 . This time-domain simulation, which consists of a discrete gust response simulated by 1-cosine disturbance wave, is given as follows:

w(t) = U ds 2 × 1 -cos πV T AS H t . (45)
This disturbance shape is similar to the one used by aircraft manufacturers to validate the load alleviation functions where V T AS = 248m/s is the true airspeed, U ds the amplitude which varies from 11.8 to 19m/s and H the scale which lies from 1.52 to 7.62m. varying gust length and amplitude, derived from certification specifications.

Figure 13 clearly illustrates that the proposed algorithm attenuates the gust effect for all configurations, especially the critical ones, e.g. those where the acceleration is the larger, and which thus is dimensioning for the aircraft design. Finally, the gust effect on the so-called worst case configuration, is given on Figure 14, where the load alleviation function clearly shows a great attenuation level on the vertical structural responses.

IV. CONCLUSION AND PERSPECTIVES

In this paper a complete illustration of the benefits brought by an active closed-loop gust load alleviation controller has been presented and demonstrated on both sub and transonic configurations. This demonstration involves many research fields, going from the manufacturing and aeroelasticity to the identification and control design. Finally, the real implementation in a research-oriented experimental set-up has been achieved. Moreover, from a methodological point of view, two main contributions are also presented. The first one consists of the development an innovative frequency-based identification procedure with eigenvalues constraints using a subspace approach 27 . The second one is the use and application of the very recently developed structured H ∞ -norm minimization oriented tools to adjust a control law that can be used -as it -in the closed-loop experimental set-up of the ONERA wind tunnel facility. The overall approach provides a baseline/rationale to the gust load attenuation problem, through a model-based control methodology. More specifically, on the considered wind tunnel set-up, a RMS attenuation level of 80% (from open to closed-loop) has been achieved in between 15 and 50Hz, with a simple linear time invariant control law. Then, using time-domain simulations, and more specifically the worst case configuration, an attenuation level of 35% was observed, which is considerable and relevant for aeronautical applications. In view of a future application on a real aircraft, the next step is to assess the overall approach on a 3D model to increase the acceptability and reach aeronautical readiness level.

Finally, a media material is made available at http://www.onera.fr/en/news/gust-effect-control-a-europeanfirst-at-onera. It illustrates the benefit achieved by the active gust load control at 25Hz (bending mode) obtained in WT experimental facility of ONERA S3Ch, in transonic configuration (Mach number 0.73), with an AoA of 0 degree. The video, recorded with a high speed camera, emphasizes the control benefit. Open-loop: from 0 to 18 seconds. Closed-loop: from 18 seconds to the end. Indeed, in the open-loop configuration, the vertical structural response of 27 Which is compared to the very powerful Loewner approach, which main drawback is the lack of control on the obtained eigenvalues. the wing are large, while, when the active control is activated, a clear diminution of the displacement level is observed.

Fig. 1 .

 1 Fig. 1. Experimental set-up (internal view of the WT). Controlled 2D airfoil and electromechanical mechanism structure (foreground), gust generator (background) and air flow stream (graphical arrow).

Fig. 3 .

 3 Fig. 3. Open-loop frequency responses of the vertical acceleration z#15 of the wing at station #15 obtained in WT experiments in response to the gust disturbance (left) and to the movable surface actuator (right). Responses for four Mach number/AoA configurations. Since it represents the transfer from the control actuator, right frame exhibit two transfer, only, one for each Mach number.

  responses -Open-loop vs. Closed-loopOpen-loop Closed-loop

Fig. 4 .

 4 Fig. 4. Wind tunnel experimental data results. Comparison of the frequency response of the singular value of the vertical accelerations at station #15, for varying Mach number, with AoA of 0 degree. Open loop: dashed blue. Closed-loop: solid red.

Fig. 5 .

 5 Fig.5. First approach: comparison of the frequency response (gain and phase) of the vertical acceleration of the wing at station #15 obtained in WT experiments (dotted blue) with models obtained via Loewner-based exact interpolation (solid red) and reduced with r = 20 (dash-dotted magenta) interpolation. Left: Mach number 0.30 and AoA 0 degree. Right: Mach number 0.73 and AoA 0 degree. Vertical lines materialize the limit of the considered data set (ω i , Φ i ) N i=1 .

Algorithm 1 2 : 2 F

 122 Subspace with LMI constraints Require: Frequency-domain data set (ω i , Φ i ) N i=1 and D a LMI-defined region. Ensure: A reduced-order model Ĥ(s) = Ĉ(sI r -Â) -1 B+ D where Λ( Â) ∈ D 1: Evaluation of à = ÂΨ by the resolution of the following LMI optimization problem: min Ã, Ψ, β subject to (35), (36), (37) β Compute  = ÃΨ -1 and Ĉ = J 3 Ûs 3: Compute B and D as{ B, D} = arg min B ∈ R n×nu D ∈ R ny×nu N i=1 Φ i -H(ω i , B, D)where Ĥ(ω i , B, D) = Ĉ(e ıωi I r -Â) -1 B + D. 4: Recover the continuous form using Tustin transform.

Fig. 7 .

 7 Fig. 7. Comparison of the maximal singular values response (gain) of the wing vertical accelerations at station #4, #12 and #15 obtained in WT experiments (dotted blue) with models obtained via Loewner approach (exact in red, approximate in pink dashed) and subspace-based with LMI constraints (solid black) with r = 16 (dash-dotted magenta) interpolation. Vertical lines materialize the limit of the considered data set (ω i , Φ i ) N i=1 .

Fig. 8 .Fig. 10 .

 810 Fig. 8. Second approach (comparison): maximal singular value response (gain) of the vertical accelerations of the wing at stations #4, #12 and #15 obtained in WT experiments (dotted blue), with model obtained via subspacebased without LMI constraints (dash-dotted black) and with LMI constraints (solid red), with r = 16. Vertical lines materialize the limit of the considered data set (ω i , Φ i ) N i=1 .

Fig. 6 .Fig. 9 .

 69 Fig.6. Second approach: comparison of the frequency response (gain and phase) of the vertical acceleration of the wing at station #15 obtained in WT experiments (dotted blue) with models obtained via subspace-based with LMI constraints (solid red) and reduced with r = 20 (dash-dotted magenta) interpolation. Left: Mach number 0.30 and AoA 0 degree. Right: Mach number 0.73 and AoA 0 degree. Vertical lines materialize the limit of the considered data set (ω i , Φ i ) N i=1 .

Fig. 11 .

 11 Fig. 11. Generalized control scheme of the ith transfer T (i) w→z (K).

Figure 13 2 )Fig. 13 . 2 )

 132132 Fig. 13. Time-domain responses in open-loop (dashed blue) and closed-loop (solid red) for all the the gust configurations obtained (H ∈ [1.52 7.62]m and U ds ∈ [11.8 19]m/s at AoA 2 degree and Mach number 0.73).

Fig. 12 .

 12 Fig.[START_REF] Ricci | Wind tunnel testing of an active controlled wing under gust excitation[END_REF]. Frequency responses of the wing vertical acceleration at station #15, in response to gust -obtained in simulation (left) and WT facility (right), without (dashed) and with (solid) gust load alleviation functions. Varying colours for varying AoA and Mach number conditions. The RMS gain is also computed over the spectrum[15 50]Hz.

Indeed, since it is a sizing criteria during the conception phase for every aircraft, gust load alleviation plays a crucial role in aeronautics and has motivated many academic and industrial researches.

Transonic case lead to turbulent and contaminated nonlinear flow, resulting in drag and a loss of efficiency. This phenomena is widely studied in flow dynamics (see e.g.[START_REF] Mettot | Quasi-laminar stability and sensitivity analyses for turbulent flows: Prediction of low-frequency unsteadiness and passive control[END_REF]).

The experimental set-up was designed and manufactured through an ONERA -Aviation Design collaboration, Aviation Design being a French S.M.E. winner of a Call for Proposal. Additional results on the facility and aeroeslasticity results can be found in[START_REF] Lepage | A complete experimental investigation of gust load: from generation to active control[END_REF].

Practical considerations about complex arithmetic are also available in[START_REF] Ionita | Lagrange rational interpolation and its applications to approximation of large-scale dynamical systems[END_REF],[START_REF] Antoulas | Model reduction and approximation theory and algorithms[END_REF].

Obviously, from the original model H, other approximation methods might be applied to obtain Ĥ. For instance, the IRKA[START_REF] Gugercin | H 2 Model Reduction for Large Scale Linear Dynamical Systems[END_REF], DARPO[START_REF] Vuillemin | Poles Residues Descent Algorithm for Optimal Frequency-Limited H 2 Model Approximation[END_REF], etc.. However, this is another issue not covered in the present paper (see e.g. the MORE toolbox developed at ONERA[START_REF] Poussot-Vassal | Introduction to MORE: a MOdel REduction Toolbox[END_REF] for numerical tools encoding these algorithms).

For clarity purpose, the noise-free case is considered, only. Reader is invited to refer to[START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF] for the noisy case.

Note that to avoid complex arithmetic, one should consider the realrepresentation of G = [Re(G) Im(G)], W = [Re(W) Im(W)] and X = [Re(X c ) Im(X c )].

The eigenvalues of a dynamical system ẋ(t) = Ax(t) lies in LMI regions defined by D if these poles satisfied the characteristic function f D .

As the sampling frequency (fs = 4096Hz) is widely superior to the main system dynamics, the continuous and discrete frequency domain responses are strictly equal for frequencies inferior to 100Hz.

Since wind tunnel experiments are time consuming and expensive, the ability in easily implementing and re-tune the control law is primordial.

Indeed, for the considered application, one is interested in a single robust linear time invariant controller instead of a parameter varying one, which would be addressed in a future step of the study.

Without lack of generalities, by designing adequate weighting filters to generate the generalized plant for transfer minimization, the H 2 objective can be -in some sense -recast as an H∞ one. As it will be made clearer in the following, this formulation also allows to handle additional constraints on the controller properties.

Five configurations of open-loop systems for Mach number going from 0.30 to 0.73, and AoA from 0 to 2 degree.

One should note that the low gain in low frequency is non standard in traditional control set-up where tracking and low frequency disturbances are usually treated. However, in the considered context, neither tracking nor step disturbances are considered, which justifies this low frequency attenuation behaviour.

Note that here the exact Loewner model is used since, even is nominally unstable, it perfectly captures the frequency data over all the considered frequency range.

Note that this dynamic is linked to the pitch behaviour (see Figure2) which is, after recurrent wind tunnel experimental tests, affected by viscous friction effect.

It was experimentally impossible to generate the same kind of test in WT due to the limited amplitude range of the gust generator.

ACKNOWLEDGEMENT

The research leading to these results has received funding from the European Union's Seventh Framework Program (FP7/2007-2013) for the Clean Sky Joint Technology Initiative under grant agreement CSJU-GAM-SFWA-2008-001.


Guest Editors Introduction: Special Issue on Inductive Logic Programming

Jesse Davis • Jan Ramon Received: date / Accepted: date This special issue focuses on the field of Inductive Logic Programming (ILP), which is a subfield of machine learning that uses logic as a uniform representation language for examples, background knowledge and hypotheses. From these roots, ILP's scope has grown to encompass different approaches that address learning from structured relational data. One notable example is the area of statistical relational learning which focuses on extending ILP to model uncertainty.

The special issue is also in conjunction with the 24th International Conference on Inductive Logic Programming (ILP), which was held from September 14th to 16th, 2014, in Nancy, France in co-location with ECML/PKDD-2014. To avoid the redundancy between the conference proceedings and the special issue, authors with an accepted paper at ILP were asked to either have their paper appear in the conference proceedings or submit an extended version of the paper to the special issue. While associated with the ILP conference, there was an open call for submissions to this special issue.

The special issue received six submissions of which three were originally submitted to the ILP conference. Ultimately, four were accepted to appear in the special issue. The papers offer a nice reflection on the strengths of ILP and relational learning and where the field is headed. Namely, the articles build off ILP's established track record of being particularly well suited to addressing important applications and the vibrant recent work that focuses on modeling uncertainty in relational data.

One of the first application areas where ILP gained significant traction was in analyzing molecules, in particular for the task of drug design. Drugs are small molecules that affect disease by binding to a target protein in the human body. Approaches to drug design depend on whether properties and/or structure of the drug or target are known. Previously, ILP has been successfully applied to identify properties of a drug molecule responsible for it binding to KU Leuven, Dept. of Computer Science, Celestijnenlaan 200A, 3001 Heverlee, Belgium E-mail: {Jesse.Davis, Jan.Ramon}@cs.kuleuven.be the target protein. In contrast, the article "ILP-assisted De Novo Drug Design" by Kaalia, Srinivasan, Kumar, and Ghosh addresses this problem from a new perspective. Instead of analyzing the properties of potential drug molecules, it attempts to infer information about the active site of a protein. It uses ILP to identify common interaction patterns between a series of probes (small molecules) and a set of proteins related to the target. Technically, this paper exploits several well-known strengths of ILP. First, it exploits the ability to incorporate domain-knowledge in order to guide and customize the search process. Second, ILP returns learned models that are easy for experts to interpret.

The second paper, "Transductive Hyperspectral Image Classification: Towards Integrating Spectral and Relational Features via an Iterative Ensemble System" by Appice, Guccione, and Malerba combines several ideas from (statistical) relational learning with other machine learning techniques to tackle the practical problem of labeling pixels in hyperspectral images. The paper exploits two important ideas from relational learning to capture the spatial correlations in the data. First, it makes use of domain-specific relational features to capture information about a pixel's spatial neighborhood. Second, it employs an iterative collective classification scheme to jointly label the pixels to exploit the fact that neighboring pixels are more likely to have the same label. The authors find that their proposed framework results in improved accuracy on the standard benchmark data sets for this problem. Furthermore, the authors carefully assess the contribution that each subcomponent of the learning system makes to its overall performance.

The third paper, "Fast Learning of Relational Dependency Networks" by Schulte, Qian, Kirkpatrick, Yin, and Sun, falls into the area of statistical relational learning. It tackles the problem of learning the structure of a relational dependency network (RDN) from data. The proposed approach first learns a Bayesian network (BN), which is then converted into a relational dependency network. The paper provides a closed-form method for deriving the RDN's parameters from the learned BN. Theoretically, the authors show sufficient and necessary conditions for this transformation to result in a consistent RDN. The main empirical advantage is that the proposed approach provides improved scalability.

The final paper, "Probabilistic Logic Programming for Hybrid Relational Domains" by Nitti, De Laet, and De Raedt, combines the move towards incorporating uncertainty into ILP and focus on applications. The paper studies probabilistic logic programming, which incorporates probabilistic constructs into logic programming. Technically, it shows how to model both discrete and continuous variables within probabilistic logical programs for both static and dynamic settings. The paper then proposes and establishes the theoretical correctness of a particle-filter based inference algorithm for the language. Among other empirical results, the paper shows how the framework can successfully model real-world tracking problems that are common in robotics and vision. We hope that the readers enjoy the papers in this issue. We want to thank all the authors and reviewers for the time spent in putting this special issue together. Finally, we want to thank Filip Zelezny for serving as the action editor for "Probabilistic Logic Programming for Hybrid Relational Domains."