A new frequency-domain subspace algorithm with restricted poles location through LMI regions and its application to a wind tunnel test

ABSTRACT In this paper, an innovative method is presented to identify models with a modified frequency-domain subspace method. This new approach allows to introduce in the subspace resolution constraints on the identified model poles location. Here, a very general formulation is proposed to take into account regions in continuous/discrete map. This formulation is based on an LMI (linear matrix inequalities) description where the stability domain represents a particular case. These LMI constraints are combined with the frequency-domain subspace resolution to obtain identified models whose poles are situated in the specified LMI regions. This approach is benchmarked with the Loewner one, which belongs to the class of frequency-domain input–output model identification and approximation methods. Besides the fact that they both belong to the data-driven model approximation class, they result to have slightly different objectives and show complementary performances. This discussion is illustrated in practice with experimentations that have been performed for the identification and control of the gust disturbance over a 2D wing span, from sub to transonic, in a wind tunnel facility.


Forewords and problem description
The analysis, control and optimisation of complex dynamical systems are very challenging tasks when designing and enhancing the behaviours of advanced systems. In practice, most of the associated techniques are grounded on a mathematical description of the system by means of a dynamical model. Indeed, the main existing control-oriented techniques require an accurate model whose primary function is to accurately catch the main dynamical behaviours. Utilising a low-order dynamical model as a surrogate in order to reproduce the inputoutput behaviour of the original model or system, is, for many practical reasons (e.g. intensive simulation), a key element in any model-based control design, analysis and optimisation procedures. This statement holds true in most of the engineering fields such as aerospace (Poussot-Vassal & Roos, 2012), micro-electronics, building structure, fluid flow (Poussot-Vassal & Sipp, 2015), and biological systems.
Within the many model approximation methods (Antoulas, Lefteriu, & Ionita, 2016), when the internal dynamics are not accessible, e.g. when a dynamical model is not a priori known or easily accessible, the so-called CONTACT F. Demourant fabrice.demourant@onera.fr data-driven input-output techniques are the most appropriated methods to approximate the system responses. Among these approximation methods, two broad families can be considered: the time-domain and the frequencydomain ones. The former one has driven a lot of attention providing solutions close to the system identification domain and is not necessarily appropriate to identify flexible structure and more generally control-oriented models. Usually it is necessary to catch the system response in a time-domain window, which can represent a difficult task for weakly damped models with the drawback to have finally a very large number of data. On the other hand, the latter is well appropriated for our problem. It is rather straightforward to take into account control objectives with a frequency-domain approach. Indeed, it is possible to take into account bandwidth and to use a more dense frequency-domain points set to identify frequency-domain peaks which correspond to flexible modes. Besides, this kind of approach can be involved with a very limited frequency-domain points. All these points represent, from both methodological and practical points, strong advantages. From a more formal aspect, let us consider that we are given a set of N ∈ N * frequencies and responses as where i ∈ C n y ×n u and ω i ∈ R + . (1) The objective is to construct a low-order rational ODE (or DAE) LTI model with realisationŜ : (Â,B,Ĉ,D) of the formŜ : wherex(t ) ∈ R r :=X is the internal variables, u(t ) ∈ R n u := U the input function and t ) ∈ R n y :=Ŷ the output function withÂ ∈ R r×r ,B ∈ R r×n u ,Ĉ ∈ R n y ×r and D ∈ R n y ×n u as constant matrices and which corresponding transfer function iŝ such thatĤ(s) (orŜ) well matches the frequency samples {ω i , i } N i=1 , and hopefully reproduces the real transfer, denoted by H (in other words,ŷ should accurately reproduce y, the output of H).

Remark 1.1 (DescriptorÊ):
In the general formulation of the Loewner approach, a descriptorÊ is considered and the objective is to obtain a low-order model with reali-sationŜ : (Ê,Â,B,Ĉ,D). As this descriptor is not completely used in our case, it is assumed thatÊ is invertible and thus one can useÊ ← I r ,Â ←Ê −1Â andB ← E −1B . Remark 1.2 (Model identification/approximation connections): In this discussion, we do not intend using the presented techniques (interpolation and subspace) for a model reduction purpose since we consider that no realisation are available (when a realisation is available, reader should refer to, e.g. Gugercin, Antoulas, & Beattie (2008) or to Poussot-Vassal & Vuillemin, 2012). We rather consider the proposed technique for a combined identification and model reduction objective. This objective makes sense when a black box simulation code or an experimental set-up is available to provide outputs, given some forced inputs. Therefore, we believe that this framework is at the boundary between identification and model approximation. Remark 1.3 (About the frequency-data considerations and collection): In this paper, we focus on the frequency data, only. Basically, this representation is rather more appropriate for control-oriented identification to take into account more explicitly and easily performance and robustness specifications such as bandwidth, neglected dynamics, and high accuracy level for the frequency-domain peaks identification, thanks to a specific frequency-domain gridding (Demourant & Ferreres, 2002). Besides, in the context of the considered electro-mechanical application, frequency-domain responses were naturally used to evaluate the load alleviation obtained between the closed and open loops. This kind of representation is very classical and useful for active control of flexible modes. Remark 1.4 (Considerations about linear descriptor systems): In this note, we will restrict the model to be linear (descriptor) only. This limitation is also a differentiating point with advanced model identification techniques that are able to encapsulate nonlinear phenomena (see e.g. Janot, Young, & Gautier, 2016). However, as made clearer throughout the paper, this restriction is balanced by some strong and formal results, accompanying the proposed methodology (e.g. interpolation, pole location, optimality properties). Therefore, once again, we believe that this should be of interest in some applications.

Outlines
In what follows, given a set of input-output frequencydomain data as in (1), we propose an innovative technique to constraint poles of the identified models in LMI regions. These constraints are specifically implemented with a frequency-domain subspace algorithm. We will show that the relation which allows to extract the estimated matrixÂ can be rewritten as an optimisation problem and can be recast, thanks to a specific variable substitution as an LMI optimisation problem when LMI poles constraints are added. Of course, this technique treats stability problem of identified models, but more generally allows to constraint poles of the identified models in specific region of the s-plane or z-plane, where stability appears as a particular case. This technique is particularly adapted for control-oriented models. This new approach will be compared with another frequencydomain approximation, namely the Loewner approach. We propose a discussion of these two complementary methods for combined MIMO model identification and reduction methods. More specifically, Section 2 focuses on the model interpolation and approximation in the socalled Loewner framework, while Section 3 focuses on the identification in the subspace framework, using an original and new LMI-based poles placement methodology to enforce some model spectral properties. Then, both methods are compared from a methodological point of view on numerical test beds available in the literature in Section 4. Then, Section 5 presents a more complete comparison of the aforementioned methods on an experimental wind tunnel set-up, where the objective is to identify (and control) the effect of a gust disturbance on a 2D wing profile (for varying configurations). Conclusive remarks are drawn in Section 6.

Interpolatory methods for data-driven model approximation
In the first step, let us introduce the Loewner framework and some of its associated properties. The exposition in this section follows Mayo and Antoulas (2007) and Ionita (2013). It is also inspired from a recent and very didactic survey of Antoulas et al. (2016). These papers should be consulted for an overview of interpolatory reduction methods. The interpolation problem is stated as follows.

Problem 2.1 (Data-driven interpolation): Given left interpolation driving frequencies {μ
, find a (low order) systemŜ of the form (2), such that the resulting transfer functionĤ(s) of the form (3) is an (approximate) tangential interpolant of the data, i.e. satisfies the following left and right interpolation conditions: Besides the very general appearance of this problem and the large variety of solutions that can be derived from it, interestingly, the Loewner framework results to be a very powerful tool to address it. In the following section, we derive the main procedure to find such aĤ.

Loewner framework for linear systems data-driven interpolation
As introduced in Section 1, the main ingredient to solve Problem 2.1 is the Loewner matrix, which was developed in a series of papers (see e.g. Mayo & Antoulas, 2007). In the sequel, we provide an overview of the Loewner framework, in the MIMO case. In this section, we formulate the main results for the above introduced data-driven interpolation problem. This will involve the more general tangential interpolation problem defined in (4). By considering λ i and μ j , to be distinct, the right and left interpolation data can be gathered as follows: One should note that the pulsation ω i are the so-called driving frequencies and are directly linked with and M matrices. The associated 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: Interestingly, the L and L σ are the solutions of the following Sylvester equations: As a major result, if the data (5) are sampled from the fictitious n u inputs n y outputs transfer function H δ (s) = C δ (sE δ − A δ ) −1 B δ , with E δ , A δ ∈ R n×n , and by denoting the left interpolation data v * j = l * j H(μ j ) and right interpolation data w i = H(λ i )r i , it follows that: Moreover, if one defines O q ∈ C q×n and R k ∈ C n×k matrices as it follows that the Loewner pencil constructed from tangential data has a system theoretic interpretation in terms of the tangential controllability and observability matrices. Indeed, following (9), the following holds true: A very powerful result in the Loewner framework is then stated in the following Lemma.

Lemma 2.1 (McMillan degree and rank):
Given (tangential) samples of a rational function defined in terms of a minimal descriptor realisation S δ : (E δ , A δ , B δ , C δ ), construct the associated Loewner and shifted Loewner matrices L, L σ . Assuming that we have enough samples, and that the left, right tangential directions l j and r i are chosen so that O q and R k have full rank, the following holds: where n is the McMillan degree of the underlying rational function. Moreover, The above Lemma states that the Loewner pencil encodes the minimal McMillan degree of a system. Now, we are ready to recall the main result concerning the construction of interpolants using the Loewner pencil. The following theorem summarises the interpolatory model construction from data. Theorem 2.1 (Loewner framework; Mayo & Antoulas, 2007): Given right and left interpolation data as in (5), and assume that k = q and let (L σ , L) be a regular pencil where λ i or μ j are not an eigenvalue. The rational trans- (13) is a minimal descriptor realisation of an interpolant of the data, i.e. H δ (s) = W(L σ − sL) −1 V interpolates the left and right constraints (4).
Invoking the above theorem, one can construct a model S δ : (E δ , A δ , B δ , C δ ) with transfer function H δ (s) that has a McMillan degree of n ࣘ k + q, and which tangentially interpolates the data (4). To obtain a reducedorder modelŜ of order r ࣘ n that well approximates S δ , it is possible to apply an SVD as follows: where 1 ∈ R r×r , 2 ∈ R (n−r)×(n−r) and Y 1 , Y 2 , X 1 , X 2 are of appropriate dimensions. Then, the reduced-order modelĤ is simply obtained by the classical Petrov-Galerkin projection:

Properties and discussions
Let us now consider the case where more data are used due to oversampling, which is realistic for applications where one has a poor a priori knowledge of the system. Note that this case is not exactly the same as the noisy data one, which should be treated in another manner. In this case, the problem has a solution provided that, for all Then, a minimal realisation of an interpolant of the data is given by the system: where X ∈ C k×n * and Y ∈ C k×n * are the orthogonal factors of the short SVD, αL − L σ = Y X * , where ∈ C n * ×n * . Note that in this case, as in the Gramian-based approach (and the balanced realisation), the SVD aims at revealing the non-relevant state variables, which appears as zero entry singular value matrix . If one then selects r < n * as the singular value matrix order, the projected quadruple (E, A, B, C), given by is a descriptor realisation of an (approximate) interpolant of the data with McMillan degree r < n * = rank L.

Remark 2.1 (Noisy data in the Loewner framework):
The noise-free case can be encountered when the model to be identified/approximated is obtained from numerical simulations, coming from finite element, fine volume software, etc. (see Section 4). Obviously, in the case where data are subject to noise, which appears in the applicative context when sensors are included in the acquisition process (see Section 5), so far, to the best of authors' knowledge, no dedicated solution can be stated in the Loewner framework. However, an a priori filtering and selection of the data {ω i , i } N i=1 (e.g. using linear phase filters, bandpass filters, etc.) or using an iterative Loewner procedure can be used.
The above remark is one of the justifications for studying another approach based on subspaces, described in what follows, which objective is no longer to interpolate the data but rather to minimise a mismatch between the data and the rational functionĤ. This shade in the approach will lead to other interesting properties, especially in view of control design.

Main contribution: approximation methods in the subspace framework with LMI constraints
Similarly to the Loewner framework, the subspace approaches are well tailored to identify and approximate frequency-domain sampled data, by a state-space realisation as defined in (2). 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 by eigenvalue real part negativity). In this section, we propose a solution to constraint the poles of the identified model in a user-defined region using 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 knowledge, in the subspace context, the stability constraint is treated in Chui and Maciejowski (1996), Lacy and Bernstein (2003), Maciejowski (1996), McKelvey et al. (1996), and Van Gestel, Suykens, Van Dooren, and De Moor (2001. In Lacy and Bernstein (2003), 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 stability region. Approaches in McKelvey et al. (1996), Maciejowski (1996) can distort the solution and techniques proposed in Chui and Maciejowski (1996), Van Gestel et al. (2001) 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 in this article and not only stability domain. This point is illustrated in the application section by using an ellipsoid constraint to impose not only stability for discrete models but also to skip irrelevant and unrealistic under-damped poles. After some reminder about subspace methods and LMI regions, a subspace approach with LMI region constraint is proposed, and an algorithm is derived.

Preliminary results on subspace methods
Fundamental subspace-based results are briefly recalled here, but additional information can be found in Liu, Jacques, and Miller (1996) and McKelvey et al. (1996).. Let us first consider the frequency-domain discrete statespace representation: where X(ω) ∈ C r , U(ω) ∈ C n u and Y(ω) ∈ C n y 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 (19) similarly reads as where is considered, then, one has G(ω i ) = i and the following relation holds (McKelvey et al., 1996) (with r + q ࣘ N and q ࣙ r): where has full rank r for all values q ࣙ r and, finally, Relation (21) emphasises that the extended observability matrix O (22), 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 an orthogonal projection W ⊥ defined as to cancel the W term. Indeed, by right multiplying (21) with W ⊥ as defined in (23), the following is obtained: Then, by denoting an effective way to extract A and C is to use a QR and an SVD matrix decomposition of GW ⊥ ∈ R n y q×n u q , and by noticing that whereˆ 1 ∈ R r×r andˆ 2 ∈ R (n y q−r)×(n y q−r) represent the singular values, and whereÛ 1 , It is well known that if the frequency-domain data are generated by an LTI system of order r, the choice of the identified system order, and, equivalently, the choice of the highest singular values are obvious since in the singular value decomposition of R 22 (24), a clear drop appears at the (r + 1)th singular value. When frequency-domain data are obtained from experimental system, this cutoff is not necessarily clear, but this decomposition can be seen as the 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ˆ 1 . Thus, and where Let us now assume thatÂ andĈ are now known, the transfer function became affine inB andD and can consequently be obtained by solving the following standard convex optimisation problem: 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 which is a zero mean complex random variable and has bounded fourth-order moments (details are given in McKelvey et al., 1996;Pintelon, 2002).

Preliminary results on LMI regions for LTI systems
Now, before exposing the proposed algorithm, let us recall some preliminary results on LMI regions which have initially been used in the context of controller synthesis (Chilali & Gahinet, 1996;Chilali, Gahinet, & Apkarian, 1999).

Definition 3.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 (27) Let us notice that the characteristic function f D (z) = P + Qz + Q T z is an n × n Hermitian matrix. Moreover, the LMI regions are symmetric with respect to the real axis. Without lack of generalities, inequalities (28)-(32) provide a set of characteristic functions provided as LMI regions, namely left half complex plane (28), ellipse of centre c 0 with half-length axis a and b (29), vertical strip between α 1 and α 2 (30), horizontal strip between ±β (31), conic sector of tip c 0 and angle θ (32).

Re(z)
Remark 3.1 (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 (29)). The equivalent constraint for continuous systems is then obtained by applying the standard transformation s = 1 T e ln(z), where s is the Laplace variable and T s the sampling period.

Remark 3.2 (Convex sets intersection):
The convexity property is 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 functions f D 1 and f D 2 , respectively, the characteristic function Consequently, the formulation of poles placement objectives within LMI regions is given by Theorem 3.1. ;Chilali & Gahinet, 1996;Chilali et al., 1999): Let us consider a dynamical system as in (19), P = P T ∈ R r×r , Q ∈ R r×r , then A has all its eigenvalues in the LMI regions D defined by the characteristic function f D (z):

Theorem 3.1 (LMI regions
if and only if there exists a symmetric = T > 0 such that

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 poles placement of the identified model. Proposition 3.1 (Subspace with LMI regions constraints): Given P = P T ∈ R r×r and Q ∈ R r×r characterising the LMI region D defined as The eigenvalues ofÂ = ϕ † b ∈ R r×r lies in D if and only if the following problem: whereÃ =Â ∈ R r×r and = T ∈ R r×r and β ∈ R, has a solution.
Proof: Thanks to Theorem 3.1, eigenvalues ofÂ are located in the LMI region D if = T > 0 exists such that Then, a quadratic cost function has been transformed into a linear cost function under quadratic constraints. Consequently, by invoking the Schur lemma, one obtains where matrix inequalities clearly are affine in decision variables β,Ã and .
A numerical procedure to obtain the dynamical model satisfying the eigenvalue LMI constraint is sketched in what follows (see Algorithm 1).

Algorithm 1 Subspace with LMI constraints
Evaluation ofÃ =Â by the resolution of the following LMI optimisation problem: Recover the continuous form using Tustin transform.

Remark 3.3 (Handling noisy data):
In the noisy data context, to preserve consistency, a weighting matrix N is derived from the noise relative variance, which must be known a priori, for each frequency sample. The evaluation is strictly equivalent and can be included in our optimisation problem as follows:

Remark 3.4 (LMI constraints):
If no constraint 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 (McKelvey et al., 1996).
Other extensions can be investigated as in the nonlinear context (see e.g. Lacy & Bernstein, 2005;Marchesiello & Garibaldi, 2008;Noel & Kerschen, 2013). Generally speaking, the condition to satisfy in order to involve poles location is to write the problem under the following form ϕÂ = b .

Numerical applications
Now two frequency-based model identification and reduction techniques have been described, and we are here interested in comparing their performances in a qualitative manner. Indeed, authors stress that it is not the objective of claiming that an approach is better than another, but rather to point out some strengths and weaknesses, in a numerical and reproducible set-up.

Benchmarks description
To do so, the following standard academical numerical benchmarks, available in the COMPl e ib (Leibfritz, 2003) library, are used. Each one has some interesting characteristics: r The Los Angeles Hospital model (LAH) is a strictly proper and stable SISO model of order n = 48. This model has several significant dynamics which require a relatively high reduced order to be accurately reproduced. It also presents a derivative action in low frequency.
r The clamped beam model (CBM) is a strictly proper and stable SISO model of order n = 348. This model has a high gain dynamics at low frequency and several other small dynamics at higher frequencies.
Catching the low-frequency dynamic directly yields a low approximation error.
r The CD player model (CDP) is a strictly proper and stable MIMO (n y = n u = 2) model of order n = 120. This model has very different dynamics, going from lightly to highly damped modes.
r The International Space Station model (ISS1) is a strictly proper and stable MIMO (n y = n u = 3) model of order n = 270. Aside from being MIMO, this model has high gain dynamics at several frequencies which makes it particularly well suited to illustrate flexible effects and thus a high-order model is required to catch its behaviour.
r The mass damper model (CM6) is a strictly proper and stable SIMO (n y = 2, n u = 1) model of order n = 960. This model has a lot of very lightly damped modes resulting in very high and narrow peaks.
For the subspace approach, a frequency sample F s = 500 Hz is chosen for all benchmarks except with the CDP one where the frequency sample is increased to 20,000 Hz. The bandwidth considered is [0 100] rad/s for ISS1, CBM, LAH and CM6 benchmarks, and [10 1e5] rad/s for CDP one.

Noise-free use-cases
Let us first analyse the performances of both Loewner and subspace approaches in the noise-free context. The different results obtained have been gathered in Table 1, and some properties and comments are given in what follows.
The following notation is adopted: r In the Loewner approach, it is not necessary to choose the identified model order since, as described above, the algorithm is able to find the 'optimal' order which allows to exactly fit the frequency-domain data generated by the LTI model. However, to compare the two approaches, the same order has been chosen. The results in Table 1 show clearly that the residual error represented by the column relative H 2 error, obtained with a non-optimal order by the Loewner approach, is negligible; r It would be possible to choose a higher order for the identified model, and eventually the same order as the original benchmarks, but it is obviously not necessarily relevant. In the context of this article, the first objective is to obtain control-oriented models, i.e. low-order models which contain the main dynamics. For example, low frequencies and high frequencies model behaviour are not critical in the control context. Indeed, in this case, a high gain strategy for the controller in the low frequencies allows to erase the static error and a weak gain strategy for high frequencies allows to guarantee stability face to neglected dynamics. In the same way, a slight identification error on the flexible modes is perfectly and easily handled by the controller with an active control strategy (Demourant & Ferreres, 2002). Second, the relative H 2 error in Table 1 shows that the chosen order is sufficient enough to obtain a model which accurately approximates the frequency-domain response. This point is illustrated in Figures 1 and 2. Since the LAH benchmark order is relatively low (n = 48), the same order r for the identified model has been chosen; r The H 2 reduction error obtained by the Loewner approach is lower to the one obtained with the subspace approach. However, as pointed in Figures 1  and 2, even if the Loewner approach is more accurate, the identification error obtained with the subspace approach is negligible. It is almost impossible to distinguish the difference between the frequencydomain responses obtained by the Loewner and subspace approaches; r Regarding the computation time, it can be surprising to have such a difference. For the subspace approach, to identify the different models which contain a large number of flexible modes, it is necessary to increase the value of q to obtain accurate identified models. For example, for model CBM, a value of q = 4000 is necessary to catch the main dynamics. With a value of q = 3000, the result is not satisfactory and several flexible modes are missed. Of course, it is possible to tune more precisely q to improve the trade-off between computational burden and accuracy, but the former will not strongly vary. And as N > q + n, the matrix size might significantly increase. For the benchmark CM6, it is necessary to impose q = 2000 to obtain a good matching between data and the identified model. With q = 1500, some flexible modes are badly identified. However, it is important to keep in mind that the computational time in favour of the Loewner approach has been obtained on specific benchmarks which contain a large number of flexible modes. Other kinds of benchmarks could lead to a calculation time in favour of the subspace approach.
r No stability constraint has been added in the noise-free case for the subspace approach since obtained models were all stable. This result derives from the good properties of Loewner and subspace approaches, i.e. the correctness property for the subspace approach and the exact interpolation for the Loewner approach. Of course, these properties are true for the 'optimal' order, i.e. when the identified model order is equal to the reference LTI model order. Even if no formal proof can guarantee the stability of a reduced-order model, the consequence in our case is that the stability is preserved. Let us remind that this result is obtained in the noise-free case and with frequency-domain data generated by an LTI model of finite dimension.
As the first conclusion, the Loewner approach is particularly dedicated to provide a reduced-order LTI models from frequency data, which are here generated by a full-order model in the noise-free context. Indeed, the Loewner approach is an interpolation-based method which is exact when the order is not imposed. This property allows to obtain residual error which is, in general case, weaker than with the subspace approach. More generally, the Loewner approach is an efficient method to interpolate frequency-domain responses obtained by a numerical black box simulator.

Remark 4.1:
It is noticed on the CDP benchmark ( Figure 1) that high frequencies (>1e4 rad/s) are not correctly identified. It is possible to improve the result with the Loewner approach and to exactly fit the frequencydomain response with a reduced model of order r = 80. It is numerically difficult to obtain a similar result with the subspace approach since, in this case, it is necessary to impose a very high sample frequency (>1e6 rad/s), which can lead, in practice, to numerical issues.

Noisy use cases
In this section, the Loewner and subspace approaches will be compared in the noisy context. Two different noises have been considered: r An additive noise which acts on low gain frequency response and which can represent a modeling error; Finally, the noisy frequency-domain response of different benchmarks is disturbed as follows: where ε m and ε a are tuned to affect the frequencydomain responses, η m and η a are generated by the Matlab function rand.m with a standard uniform distribution on the open interval (−0.5, 0.5) and (0, 1), respectively. Then, ε m = 0.2 is chosen and ε a is tuned to act on the weak gain frequency-domain response only. The aim in this section is to evaluate the different algorithms with data subject to noise, without any a priori knowledge of it. Figure 3 illustrates the noise influence on the frequency-domain response of the benchmark LAH. Results are reported in Table 2. This table presents additional columns: r Column No physical peaks corresponds to an analysis of the peaks that appear when the algorithm tends    For the sake of completeness, in Table 3, the LMI constraints involved to identify models with the subspace approach are also provided. These constraints are applied in the discrete plane since the subspace approach provides discrete-time models. This kind of constraint leads to impose the stability since this LMI region is included in the unitary disk. Furthermore, the half-length axis b ࣘ 1 avoids to underestimate the modes damping and to delete possible weakly damped poles which can appear and which are not physical. It is important to notice that these constraints strongly depend on the sampling period (T s ). For example, if the LMI constraints between benchmarks LAH and CDP are compared, one should not consider that constraints on CDP benchmark are stronger. Indeed, it suffices to represent them in the continuous time by the transformation s = 1 T s ln(z) (see Figure 6) and to notice that, as the sample period T s is different between the two benchmarks, constraints applied on CDP benchmark to identify it are weakest.
From results gathered in Table 2, the following remarks can be made: r Concerning the subspace approach, the computational burden remains limited with the LMI resolution. More specifically, when the LMI poles placement constraints are not active, typically when the stability is imposed on an identified model which is already stable in spite of noise, the resolution time is comparable to the resolution time without LMI constraints (benchmarks ISS, CDP, CBM). In the case where LMI constraints are active, it is clear that the computation time is increased (benchmarks LAH and CM6) while still remaining acceptable. As in the previous case (noise-free case), the resolution time obtained by the Loewner approach remains largely inferior to the one obtained by the subspace approach (with or without LMI); r As expected, the Loewner approach sometimes leads to instability or weakly damped modes, which do not correspond to any structural or physical flexible modes. Except for the benchmark ISS1, identified models by the Loewner approach are unstable and/or contain parasitical dynamics. It is particularly true for the CM6, CDP and CBM benchmarks. The baseline subspace approach has a more robust behaviour with respect to the noise since only LAH and CM6 benchmarks are unstable; r The LMI resolution involved in the subspace approach allows to improve the quality of the identified models, since it is possible without any a priori knowledge of the noise, to obtain relevant and accurate identified models. Figure 7 shows the main advantage to constraint dynamics in LMI regions. In Figure 7, the identified model obtained by the Loewner approach is slightly unstable. Furthermore, this kind of results is obtained with a frequencydomain response which efficiently approximates data (see Figure 5); r More specifically, it is to be noticed that the H 2 relative error can be significant for the Loewner and subspace approaches. Still, it is important to distinguish these two approaches for the results interpretation: • Concerning the Loewner approach, the H 2 relative error essentially comes from parasitical peaks which do not correspond to any physical behaviours (see Figure 4). This kind of problem can be very pejorative for the controller design process since this one will take into account these (very) flexible modes by involving either a weak gain strategy or an active control strategy. However, these two strategies have a cost in terms of results and specification feasibility, all the more when these parasitical dynamics do not exist on the true system; • Concerning the subspace approach, the H 2 relative error essentially comes from notches in the frequency-domain responses which are not correctly fitted in high frequencies (benchmarks LAH, CBM and ISS1). Nevertheless, these points are easily handled in a closed-loop set-up. A notch   in the frequency-domain response of the system is not viewed as a constraint and do not lead to specific strategies for the controller design. And error in high frequencies domain is naturally treated by a strictly proper controller, which is a suited and classical constraint to be robust to neglected dynamics and to limit control signal effort in presence of measure noise.
r The instability obtained on benchmarks LAH, CBM and CM6 has not a negative influence on the accuracy of the frequency-domain response since it still very accurately reproduces the frequency-domain data with an unstable model. However, this property is unacceptable when it is necessary to perform time-domain simulations to compare open-loop or closed-loop responses. This is particularly the case in the aeronautical domain for gust load alleviation functions (see Section 5). Furthermore, specific analysis with time-domain simulations can be performed to evaluate the closed-loop behaviour face to saturation or dead-zone on actuators, and of course, results are radically different if the open loop is stable or unstable; r It is important to notice that the relative H 2 norm is quasi-identical for the baseline and LMI-based subspace approaches. It is perfectly coherent with Remark 3.4, where it is claimed that poles constraints do not affect consistency if these constraints are coherent with model dynamics (McKelvey et al., 1996). This is particularly true for benchmarks LAH and CM6, where LMI constraints are active. Let us notice that even with parasitical peaks, the H 2 relative errors obtained by subspace approach with and without LMI constraints are very close due to the fact that their energetic contribution is very weak; r Globally, the cumulated relative H 2 error for the three approaches in the noisy context are close, and quasi-identical between the baseline and LMI-based subspace approaches.
r The instability and/or the presence of parasitical peaks depend on models and noise level. For different noise levels, results would be different in terms of instability and weakly damped poles. But the objective here is not to evaluate quantitatively the influence of the load level on the modelling approaches but to show qualitatively the possible consequences of noisy data on the behaviour of different algorithms. And it is particularly interesting to notice that similar results in Section 5 have been obtained with experimental data which are weakly noisy.
The Loewner approach is a very efficient approach to obtain models to validate closed-loop frequency-domain response. First, because the residual error is in general case very weak, even in the noisy context. And, second, the calculation time in our applications were particularly limited in comparison with the subspace approach with or without LMI. That is why the Loewner approach is involved to obtain an LTI model to evaluate the frequency-domain load alleviation. We will see in Section 5.4 that a validation LTI model of optimal order n = 292 is obtained in a few seconds, with an exact interpolation by the Loewner approach. The subspace approach with LMI constraints is particularly dedicated to obtain relevant synthesis models of low-order models which is specifically well adapted for the controller design step. Controlled D airfoil and electromechanical mechanism structure (foreground), gust generator (background) and air flow stream (graphical arrow). Right: experimental set-up schematic view. Gust wave generated by the gust generator (left horizontal wings) attacking the controlled D airfoil (right horizontal wing). The airfoil is equipped with a movable surface and  sensors (three of them, #, #, #, are materialised by dots, and will be used either as control measurements or performance objective). v is the wind stream vector of the wind tunnel and w is the vertical gust vector caused by the gust generator.

Experimental applications
Now that the first numerical comparison has been done, allowing to draw the main advantages and drawbacks of each approach, it is interesting to perform another analysis on a real experimental test measurement set. This is shortly described in this section, where an identification, followed by a control design, has been done in a wind tunnel (WT) facility. Briefly, the purpose of this experiment was to identify the effect of a gust disturbance and movable control surface on a 2D wing profile (i.e. its vertical accelerations), for varying angles of attack and Mach numbers, going from sub to transonic velocities. Additional details are given in Lepage et al. (2015) and in Figure 8. The bandwidth [15 85] rad/s is considered here.

Application to the WT data by the Loewner approach
For the considered application, Theorem 2.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 in Figure 9, with solid lines, providing a model that perfectly matches all measurement points. Then, thanks to the Loewner matrices 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 lines).  Figure 9 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. Among others, no stability guarantee can be enforced. 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 drawback is handled in the next subsection using the subspace approach.

Application to the WT data by the subspace approach
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 = 4096 Hz. For the resolution, q = 100 is chosen. When applying Algorithm 1, the considered LMI region is an ellipse defined by constraint (29) 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 underestimate 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 in Figure 10, with r = 20 (on the same configurations as the one shown in Figure 9 for the Loewner approach). Comparing with Figure subspace approach. Then, Figure 11 compares the maximal singular values of the MIMO models obtained with 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 step, and actually used in the experimental context -see Section 5.3). It emphasises the fact that the subspace with LMI constraint obtained model (black solid thin) well matches the measurement points (blue points).  In addition, Figure 12 presents the eigenvalues of both the baseline subspace approach and the LMI-based subspace one. Interestingly, Figure 12 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 (materialised by the solid line).

Robust controller design
To complete this section, and illustrate the previous remarks concerning the stability issues and noise properties that the sub-space approach provides, let us now briefly describe the controller design procedure that has been used to control the considered phenomena. This presentation is done from a very high-level point of view (additional details are available in Poussot-Vassal, Demourant, Lepage, & Le Bihan, 2016). Regarding the identification and reduction step, it is important to note that we have obtained, after many wind tunnel tests, n s = 5 dynamical models, at varying Angles of Attack (AoA) and Mach numbers (going from sub to transonic). As the consider system is nominally stable, the models used in the control design step are the one obtained with the subspace approach with stability constraint. Indeed, the Loewner framework will be used for frequency validation purpose only.
Using these models, a controller is synthesised using a robust approach in order to obtain K , an LTI controller that is both robust to the Mach number and AoA. To address this objective, the synthesis problem is recast as a multi-channel H ∞ -norm minimisation problem. Mathematically, the objective is to find the optimal controller K such that (see also Apkarian & Noll, 2006), where K describes the considered controller family set, 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 fromw toz, i.e. from the gust to the weighted vertical accelerations at stations #4, #12 and #15. Finally, W k is an additional weighting function used: (1) first enforce transfer KW k to be stable, and thus K to be stable, and (2) second, to shape K roll-off behaviour. Since it is not the topic of this paper, interested reader should refer to Poussot-Vassal et al. (2016) for additional details.
For the sake of completeness, let just notice that u(t ) and y(t ) signals denote the control signal deflection 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 generalised plant and input/output transfer to be minimised by the H ∞ -norm objective.

Numerical and experimental validation
In practice, solving (41) with n s = 5 and a LTI controller K of order n c = 8 has been obtained using the approach described in Apkarian and Noll (2006) using the models obtained with the subspace with LMI constraints approach. The achieved attenuation level, solution of the optimisation problem (41), leads to an optimal H ∞ -norm value of 0.43473( < 1), meaning that the weighting function constraints have well been attenuated. Before moving to the implementation step, in order to assess the obtained controller performance, simulations are first done by applying the optimal K on the complete model obtained with the exact Loewner approach (Section 2). The left frame of Figure 13 provides both 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 13. By comparing the RMS gain displayed on both (left and right) frames of Figure 13, one can observe that they are quite close to each other. So it is for the frequency shape that well reduces the peak at 25 Hz (heave peak and main objective of the control). The only difference between both results is concerned with the notch-like behaviour observed in Figure 13 around 40 Hz 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. Nevertheless, this varying behaviour illustrates that the implemented control law is also robust to dynamical variations. Remark 5.1 (Implementation specifications): The synthesised controllers were implemented on a real time device which was composed of several processors and input/output boards interlinked for fast internal communication and data exchange. The input-output 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. Synthesised in the continuous time-domain, the controllers were converted into discrete state-space models model to be directly implemented. With respect to the bandwidth of the electrohydraulic actuator, a 2 kHz sampling frequency was selected.
In addition to the frequency-domain data collected in WT and compared in simulation ( Figure 13), a timedomain gust response have also been performed in simulation only using the models obtained with the subspace with LMI constraints approach. This time-domain simulation, which consists of a discrete gust response simulated by 1-cosine disturbance wave, is given as follows: This disturbance is similar to the one used by aircraft manufacturer to validate the load alleviation functions, where V TAS = 248 m/s is the true airspeed, U ds the amplitude which varies from 11.8 to 19 m/s and H the scale which lies from 1.52 to 7.62 m. Figure 14 gathers all the open-loop (dashed line) and closed-loop (solid line) responses of the vertical acceleration at station #15 to the wind gust disturbance (42) for varying gust length and amplitude, derived from certification organisms. Figure 14 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 conception. Finally, the gust effect on the so-called worst case configuration (the one which is dimensioning the aircraft) is given in Figure 15, where the load alleviation function clearly shows a great attenuation.

Conclusions
In this paper, a discussion on two combined model identification and approximation techniques, using frequencydomain data, is done. More specifically, the interpolation in the Loewner framework and the subspace approaches are discussed. Working with frequency-domain data allows identifying the considered system behaviour over a limited frequency support, avoiding overparametrisation of the model, and focusing on the purpose of the model, e.g. control, validation. We show, through both numerical and experimental benchmarks, that the subspace and Loewner approaches should not be opposed but rather combined when analysing a phenomena, since they present complementary performances. While the Loewner approach provides an excellent framework when data come from numerical noise-free simulators (indeed, it provides an accurate model in a short time), the subspace approach combined with the LMI constraints (slower and less accurate in the aforementioned case) is the perfect alternative when data come from experimental set-up, handling the noise and allowing to ensure some spectral properties of the identified model. These complementary properties have been successfully highlighted in the paper through reproducible numerical simulations. In addition, wind tunnel experiment, combining the two approaches and a successful controller design and implementation, allowing to attenuate the gust disturbance effect on a 2D wing profile, from sub to transonic configurations, provides another important result of this paper.