# Stability analysis of a set of uncertain large-scale dynamical models with saturations: application to an aircraft system

Pierre Vuillemin, Fabrice Demourant, Jean-Marc Biannic, Charles

Poussot-Vassal

► To cite this version:

Pierre Vuillemin, Fabrice Demourant, Jean-Marc Biannic, Charles Poussot-Vassal. Stability analysis of a set of uncertain large-scale dynamical models with saturations: application to an aircraft system. IEEE Transactions on Control Systems Technology, 2016, 25 (2), pp.661-668. 10.1109/TCST.2016.2558542. hal-01398012

HAL Id: hal-01398012

https://hal.science/hal-01398012

Submitted on 4 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stability Analysis of a Set of Uncertain Large-Scale Dynamical Models With Saturations: Application to an Aircraft System

P. Vuillemin, F. Demourant, J.-M. Biannic, and C. Poussot-Vassal

**Abstract**—From a sparse set of large-scale linear time-invariant dynamical models, a methodology to generate a low-order parameter-dependent and uncertain model, with guaranteed bounds on the approximation error, is first obtained using advanced approximation and interpolation techniques. Second, the stability of the aforementioned model, represented as a linear fractional representation and subject to actuator saturation and dynamical uncertainties, is addressed through the use of an irational multiplier-based integral quadratic constraint approach. The effectiveness of the approach is assessed on a complex set of aeroservoelastic aircraft models used in an industrial framework for control design and validation purposes.

**Index Terms**—Aeroelastic aircraft stability, integral quadratic constraint (IQC), model approximation, model interpolation.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Global process of the proposed approach (Algorithm 1).</figcaption></figure>

## I. INTRODUCTION

Many techniques have been developed to model, control, and assess the stability and performance of dynamical systems. When complex systems are considered, dedicated numerical software is usually used to accurately reproduce their dynamical behavior. The obtained models then result in large-scale ones equipped with a prohibitively high number of variables. Although complex models have a high degree of likeness with reality,¹ in practice, due to finite machine precision and computational burden, they are problematic to manipulate. This is the case in many engineering fields such as aerospace (e.g., aircraft [1], satellites, launchers, and fluid flow mechanics), civilian structures, and electronics (e.g., [2]), where control engineers have to cope with many practical problems including lightly damped modes and nonlinear actuator(s). Moreover, parametric uncertainties usually affect such models, accounting for variabilities and uncertainties. In most cases, the parametric dependency is not a priori known and local linear models, representing the system at frozen configurations, are often considered.

Let us consider a model $G(\theta)$ of a physical dynamical system, which smoothly depends on a parameter $\theta \in \mathbb{R}^p$. This model is assumed to be known only through its linearized models $G_i$ at some parametric points $\theta_i$ ($i = 1, \dots, n_s$).

Let $G_i$ be asymptotically stable large-scale linear time-invariant (LTI) dynamical models given by the state-space realizations

$$G(\theta_i) \stackrel{\text{in.}}{=} G_i: \begin{cases} \dot{x}_i(t) = A^{(i)}x_i(t) + B_1^{(i)}w(t) + B_2^{(i)}u(t) \\ z(t) = C_1^{(i)}x_i(t) + D_{11}^{(i)}w(t) + D_{12}^{(i)}u(t) \\ y(t) = C_2^{(i)}x_i(t) + D_{21}^{(i)}w(t) + D_{22}^{(i)}u(t) \end{cases} \quad (1)$$

where $x_i(t) \in \mathbb{R}^{n_i}$, $w(t) \in \mathbb{R}^{n_w}$, $u(t) \in \mathbb{R}$, $z(t) \in \mathbb{R}^{n_z}$, and $y(t) \in \mathbb{R}^{n_y}$ are the states, exogenous input, single control input, performance output, and measurement signals, respectively. Moreover, let be given a robust $n_k$th-order LTI controller $K = (A_K, B_K, C_K, D_K)$ with transfer $K(s) = C_K(sI_{n_k} - A_K)^{-1}B_K + D_K$, looped between $y(t)$ and $u(t)$, that ensures some robustness and performance specification(s) for all the $n_s$ models.

The problem of assessing the stability of such a high-dimensional controlled system over the continuum of parametric variations, when the single control input $u(t)$ is subject to saturations, is addressed here. To this aim, as clarified in the rest of this paper and following Fig. 1 and Algorithm 1, a three-step methodology is proposed: 1) approximate the $n_s$ dynamical models and bound the mismatch error; 2) perform (inexact) interpolation of the reduced-order models with interpolation error bounds; and 3) assess the stability of the closed-loop model over both parametric variations and control input saturation limitations.<fnref n="2" />

In comparison with the contributions of [1] and [3], the proposed approach is accompanied with both approximation (step 1) and interpolation (step 2) errors. Hence, the $\mu$ (structured singular value) and integral

<footnote>The authors are with ONERA–The French Aerospace Lab, Toulouse 31055, France (e-mail: contact@pierre-vuillemin.fr; fabrice.demourant@onera.fr; jean-marc.biannic@onera.fr; charles.poussot-vassal@onera.fr).</footnote>

<footnote n="2">²Note that in practice, people usually reduce and perform the analysis in a trial-and-error way, which is of course tedious and time consuming.</footnote>
quadratic constraint (IQC) analysis (step 3), respectively, provide sufficient stability conditions for the whole set of closed-loop models without and with saturation. This represents the main contribution of this paper. It is also worth mentioning that the irrational multipliers-based approach developed in step 3 is an extension of [4]. It is shown that no solution is obtained by rational multiplier, and only a frequency-domain approach can be used here to assess the closed-loop stability.

This paper is organized following the schematic of Fig. 1. In Section II, the main result, i.e., the procedure to assess the stability of a set of large-scale models looped with a control law subject to saturations, is described. Section III illustrates the proposed procedure on a complex large-scale aeroservoelastic business jet aircraft model for various flight configurations, looped with an antivibration controller. Conclusions are given in Section IV.

**Notations:** Given three operators $P(\cdot)$, $M(\cdot)$, and $\Delta(\cdot)$ of compatible dimensions, the lower and upper linear fractional transformations (LFTs) are, respectively, defined (for appropriate partitions of $P$ and $M$) by $\mathcal{F}_l(P, \Delta) = P_{11} + P_{12}\Delta(I - P_{22})^{-1}P_{21}$ and $\mathcal{F}_u(M, \Delta) = M_{22} + M_{21}\Delta(I - M_{11})^{-1}M_{12}$. The star product $\star$ of $P$ and $M$ is defined by

$$P \star M = \begin{bmatrix} \mathcal{F}_l(P, M_{11}) & P_{12}(I - M_{11}P_{22})^{-1}M_{12} \\ M_{21}(I - P_{22}M_{11})^{-1}P_{21} & \mathcal{F}_u(M, P_{22}) \end{bmatrix}. \quad (2)$$

Given a matrix $M \in \mathbb{C}^{p \times m}$, $M_{j,k} = M(j,k)$ (with $1 \le j \le p$ and $1 \le k \le m$) denotes the scalar coefficient in the jth row and kth column of $M$, $M^*$ denotes the conjugate transpose of $M$, and $\bar{\sigma}(M)$ is its largest singular value. The frequency-limited $\mathcal{H}_2$-norm, denoted by $\mathcal{H}_{2,\Omega}$-norm, is defined as the restriction of the $\mathcal{H}_2$-norm over the interval $\Omega = [0, \omega]$ with $\omega \in \mathbb{R}_+$, where $\mathbb{R}_+$ denotes the set of positive real numbers. Given an asymptotically stable LTI model realization $\mathbf{H}$ with transfer function $H(s)$, $\|H\|_{\mathcal{H}_{2,\Omega}} := (1/\pi \int_{\Omega} \|H(j\omega)\|_F^2 d\nu)^{1/2}$ [5], [6].

## II. MAIN RESULT: STABILITY GUARANTEE OF A SET OF LARGE-SCALE MODELS SUBJECT TO INPUT SATURATIONS

With reference to Fig. 1, the proposed contribution, in three steps, is summarized in Algorithm 1. More specifically, an optimal frequency-limited approximation algorithm is first applied followed by the creation of a frequency-dependent mismatch bound (step 1, Section II-A), then the interpolation and transformation into a linear fractional representation (LFR) structure is achieved (step 2, Section II-B), and finally, the stability of the overall uncertain parameter-dependent model is first assessed thanks to a $\mu$ analysis, and then, when subject to control input saturation, through a novel IQC technique (step 3, Section II-C).

### A. Multi-LTI Model Approximation and Error Bound

Generally speaking, the main objective of the approximation step is to capture, with a stable low-order model, the initial large-scale model most relevant dynamics. Various approaches exist for the approximation of large-scale LTI models (see [7]

**Algorithm 1 Global Procedure**

**Data:** $G_i$ ($i = 1, \dots, n_s$) describing a system at various frozen parameters combinations values $\theta_i \in \mathbb{R}^p$ and a robust LTI controller $K$.

**Result:** Stability assessment

**begin step (i) (Section II-A)**

* Compute $\hat{G}_i(s)$ ($i = 1, \dots, n_s$) such that

$$\hat{G}_i := \arg \min_{H \in \mathcal{H}_\infty, \text{rank}(H)=r} \|G_i - H\|_{\mathcal{H}_{2,\Omega}} \quad (3)$$

* Determine a low-order weighting function $W(s)$ s.t. $\forall i = 1 \dots n_s, \exists \Delta_{R_i} \in \mathcal{H}_\infty, \|\Delta_{R_i}\|_{\mathcal{H}_\infty} \le 1$ and:

$$F_i(s) = \hat{F}_i(s) + W(s)\Delta_{R_i}(s) \quad (4)$$

with $F_i(s) = K(s)G_i(s)$ and $\hat{F}_i(s) = K(s)\hat{G}_i(s)$.

**return** A set of reduced-order approximations $\hat{F}_i(s)$.

**begin step (ii) (Section II-B)**

* Compute a parameter-dependent LFR approximation $\hat{P}(s)$ associated with the normalized and lowest-size block-diagonal parametric structure $\Theta(\theta)$ such that, for each parametric configuration $\Theta_i = \Theta(\theta_i)$, there exists a real-valued norm-bounded structured uncertainty $\Delta_P$ capturing the interpolation errors such that:

$$\hat{F}_i(s) = \mathcal{F}_u(\hat{P}(s), \text{diag}(\Theta_i, \Delta_P)) \quad (5)$$

* Combine (4) and (5), construct $P(s)$ including all errors, where $\Delta = \text{diag}(\Theta_i, \Delta_P, \Delta_R(s))$, such that,

$$F_i(s) = \mathcal{F}_u(P(s), \Delta) \text{ with } \|\Delta\|_{\mathcal{H}_\infty} \le 1 \quad (6)$$

**return** A low-order uncertain LFR model $\mathcal{F}_u(P(s), \Delta)$ covering the initial set $\{F_i(s)\}_{i=1 \dots n_s}$.

**begin step (iii) (Section II-C)**

* Close the open-loop LFR model $P(s)$ without input saturation, build the standard form $M(s) - \Delta$ and check robust stability by $\mu$ test:

$$\forall \omega \ge 0, \mu_{\Delta}(M(j\omega)) \le 1. \quad (7)$$

* Close the open-loop LFR model $P(s)$ with input saturation to get an augmented nonlinear standard form $M(s) - \text{diag}(\varphi, \Delta)$ and check robust stability by an IQC-based analysis test.

$$\forall \omega \in \mathbb{R} \ [M(j\omega)^* I] \Pi(j\omega) [M(j\omega)^* I]^* < 0 \quad (8)$$

**return** A stability proof of the input-saturated closed-loop large-scale models.

for an overview), and one of them consists in formalizing the model approximation problem as an optimization one. The problem then consists in finding a reduced-order model that minimizes a given norm of the approximation error.

In the literature, the $\mathcal{H}_2$-norm has often been considered and several methods are now available to address the corresponding optimal $\mathcal{H}_2$ model approximation problem [8], [9]. Yet, in many cases, considering a limited
frequency interval only is more relevant since the following conditions hold.

1) The system dynamics might not be perfectly known over the whole frequency domain, meaning that the model is inaccurate in some frequency intervals. Discarding these areas enables us to increase the approximation accuracy where the initial model is accurate.

2) Controllers are usually designed to act over a limited frequency interval (due to actuator's bandwidth or to prevent them from disturbing nonmodeled dynamics), which means that a precise knowledge of the dynamics over the whole frequency domain is not necessarily useful.

From the authors' point of view, the optimal approximation over a bounded frequency interval enables us to elegantly translate these practical considerations and is therefore preferred here. It is addressed through the use of the frequency-limited $H_2$-norm in Section II-A1. Yet it is worth noting that the overall methodology summarized in Algorithm 1 does not depend on the approximation strategy as the approximation error is bounded in Section II-A2.

1) Optimal Frequency-Limited $H_2$ Model Approximation: Using the $H_{2,\Omega}$-norm, one can formulate the approximation over a bounded frequency interval as an optimization problem. More specifically, given an asymptotically stable nth-order large-scale model $G$ and a frequency interval $\Omega$, the optimal $H_{2,\Omega}$ model approximation problem consists in finding a reduced-order model $\hat{G}$ of order $r \ll n$, which minimizes the $H_{2,\Omega}$-norm of the approximation error $G - \hat{G}$

$$ \hat{G} = \arg \min_{H \in \mathcal{H}_{\infty}, \text{rank}(H)=r} \|G - H\|_{\mathcal{H}_{2,\Omega}}. \quad (9) $$

Here, Problem (9) is addressed using the method called descent algorithm for residues and poles optimization (DARPO), proposed in [6]. It relies on the poles-residues formulation of the $H_{2,\Omega}$-norm [10] and finds the poles and associated residues of the reduced-order model that satisfy the first-order optimality conditions associated with (9).

With reference to Algorithm 1 (step 1), the approximation algorithm is applied to each large-scale model $G_i$, $i = 1, \dots, n_s$ resulting in $n_s$ small-scale models $\hat{G}_i$ minimizing the $H_{2,\Omega}$-norm of the approximation error with $G_i$, as stated in (3).

Note that both the approximation order $r$ and the frequency interval $\Omega$ are tuning parameters that depend on the considered application. Yet, as mentioned before, the frequency interval $\Omega$ may be chosen as the interval that contains the most relevant dynamics of the physical systems. Observing the decay of the eigenvalues of the product of the frequency-limited gramians $P_{\Omega}Q_{\Omega}$ (see [11, Ch. 4]), which may be viewed as the Hankel singular values in the frequency-limited case, can give an idea of the adequate approximation order $r$.

The stability analysis must take into account the error induced by the approximation step. For that purpose, a low-complexity model upper bounding the worst approximation error is built in the next section.

2) Bound on the Approximation Error: Let us denote by $F_i(s) = K(s)G_i(s)$ and $\hat{F}_i(s) = K(s)\hat{G}_i(s)$ the open loops

from the inputs of the large- and small-scale models to the output of the controller<fnref n="3" /> $K$. Let us denote the order of $\hat{F}_i(s)$ as $n = r + n_K$. The objective of this section is to model the approximation error $\Sigma_i(s) = F_i(s) - \hat{F}_i(s)$ ($i = 1, \dots, n_s$) as a low-order additive output uncertainty. More specifically, one seeks a low-order filter $W(s)$ such that $\forall i = 1, \dots, n_s, \exists \Delta_{R_i} \in \mathcal{H}_{\infty}$ with $\|\Delta_{R_i}\|_{\mathcal{H}_{\infty}} \le 1$ and $F_i(s) = \hat{F}_i(s) + W(s)\Delta_{R_i}(s)$.

Then the stability of the set of uncertain models $\{\hat{F}_i(s) + W(s)\Delta_R(s), \|\Delta_R\|_{\mathcal{H}_{\infty}} \le 1\}$ implies the stability of the finite set of models $\{F_i(s)\}_{i=1,\dots,n_s}$. Note that any invertible filter $W(s)$ such that

$$ \max_{i=1,\dots,n_s} \|W^{-1}\Sigma_i\|_{\mathcal{H}_{\infty}} \le 1 \quad (10) $$

can be used as one can always exhibit $\Delta_{R_i}(s) = W^{-1}(s)\Sigma_i(s)$ such that $F_i(s) = \hat{F}_i(s) + W(s)\Delta_{R_i}(s)$.

The design of $W(s)$ then consists in a tradeoff between complexity and conservatism. Indeed, one must find a $W(s)$ which is both an accurate modeling of the worst approximation error and which complexity (order) is reasonable. For instance, $W = \max_{i=1,\dots,n_s} \|\Sigma_i\|_{\mathcal{H}_{\infty}}$ obviously satisfies (10). Yet, it does not offer an accurate modeling of the approximation error and might therefore be too conservative for stability analysis. A direct approach to design $W(s)$ satisfying (10) would consists in using nonsmooth $\mathcal{H}_{\infty}$ optimization tools to solve $\min \|W\|_{\mathcal{H}_{\infty}}$ such that $\|W^{-1}\Sigma_i\|_{\mathcal{H}_{\infty}} \le 1$ (for $i = 1, \dots, n_s$). Yet depending on the application, the errors $\Sigma_i$ might be too large for such an approach to be tractable. In those cases, a heuristic approach may then be preferable.

## B. Bounded-Error Reduced-Order LFR Model Generation

Consider the parametrically dependent set $\{\hat{F}_i(s)\}_{i=1,\dots,n_s}$ of reduced-order models obtained above, the objective is now to derive a limited-size LFR such that $\mu$ and IQC-based analysis tools can then be applied. In the general case, involving several parameters ($\theta \in \mathbb{R}^p$), the $n_s$ equations (5) must be solved for a parametric structure, e.g., $\Theta_i = \text{diag}(\theta_{1i}, I_{n_{\theta_1}}, \dots, \theta_{pi}, I_{n_{\theta_p}})$, whose size $n_{\Theta} = \sum_{k=1}^p n_{\theta_k}$ should be kept as small as possible. Efficient solutions, based on multivariate sparse polynomial or rational interpolation techniques, are detailed in [1], [12], and [13].

In the case of a scalar parameter ($\theta \in \mathbb{R}$), a specific technique can be developed to compute low-order LFR models whose $\Delta$-block will both include the parametric variations ($\Theta$) and a normalized real-valued uncertain operator ($\Delta_P$). The latter is introduced to cover the interpolation errors as illustrated by (5). The proposed technique based on a polynomial state-space data interpolation approach can be decomposed into three steps, which are briefly presented next.

1) Step 1: Models Rewriting in a Rescaled Companion Form: Reduced-size LFR models are easier to obtain when all varying data appear in a limited number of rows (or columns) of each state-space representation. A companion form is then a good choice but unfortunately leads to ill-conditioned matrices as the system order increases. As also proposed in [14], a rescaled companion form will then be used. Using the

<footnote n="3">³The controller is included here to be coherent with the interpolation step of Section II-B.</footnote>
notation $\hat{F}_i(s) = C_i(sI_n - A_i)^{-1}B_i + D_i$, the system is rewritten as

$$ \left( \frac{A_i}{C_i} \middle| B_i \middle| \frac{D_i}{C_i} \right) = \begin{pmatrix} 0 & \lambda_1 & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & \lambda_{n-1} \\ \frac{a_1^{(i)}}{c_1^{(i)}} & \frac{a_2^{(i)}}{c_2^{(i)}} & \dots & \frac{a_n^{(i)}}{c_n^{(i)}} \\ \hline \frac{b^{(i)}}{d^{(i)}} \end{pmatrix} \quad (11) $$

where the scaling variables {$\lambda_k$}$_{k=1...n-1}$, with the help of standard numerical balancing techniques, are tuned to optimize the average condition number of each matrix $A_i$. Note that the standard companion form is recovered for $\lambda_k = 1$.

**Remark 1:** In the context of LFR modeling, the above description is of high interest since the varying state-space data all appear in the last two rows. Assuming that every coefficient is approximated by a pth-order polynomial, the size of $\Theta = \theta I_{n\theta}$ will then be limited to $n_{\Theta} = 2p$.

2) Polynomial Interpolation With Guaranteed Error Bounds: Let us denote by $Y_i$ the last two lines in (11)

$$ Y_i = \begin{pmatrix} a_1^{(i)} & a_2^{(i)} & \dots & a_n^{(i)} & b^{(i)} \\ c_1^{(i)} & c_2^{(i)} & \dots & c_n^{(i)} & d^{(i)} \end{pmatrix} \in \mathbb{R}^{2 \times (n+1)} \quad (12) $$

and focus on the polynomial approximation of the finite set {$Y_i$}$_{i=1...n_s}$ with guaranteed and minimized error bounds. Given $p$, the order of the polynomial, the problem reduces to the determination of an error matrix $E \in \mathbb{R}_{+}^{2 \times (n+1)}$ and a set of matrices {$X_q$}$_{q=0...p}$ such that the nonnegative entries of $E$ are minimized under the following linear constraints (with $j=1, 2$ and $k=1...n+1$):

$$ \left| \left[ X_0 + \sum_{q=1}^{p} \theta_i^q X_q - Y_i \right]_{j,k} \right| \le E_{j,k}, \quad i = 1 \dots n_s. \quad (13) $$

The above optimization problem is easily solved by any standard linear programming solver. However, the order $p$ of the polynomial should be carefully chosen. Low orders will indeed result in rough approximations yielding conservative models with large entries in $E$. Conversely, high-order polynomials will improve the accuracy at the interpolation points. Moreover, critical oscillations are likely to appear between the interpolation points when the difference $n_s - p$ gets too small. This issue and possible remedies are further discussed in Section III.

3) LFR Modeling:

**Proposition 1:** From inequalities (13), E-dependent shaping matrices $U(E)$ and $V(E)$ of appropriate dimensions and a bounded real-valued block-diagonal uncertain operator $\Delta_P$

$$ \Delta_P = \mathbf{diag}(\delta_{p_1} I_{n_{\delta p_1}}, \dots, \delta_{p_r} I_{n_{\delta p_r}}) \quad (14) $$

can be easily defined such that the function

$$ \mathcal{Y}(\theta, \Delta_P) = X_0 + \sum_{q=1}^{p} \theta^q X_q + U(E) \cdot \Delta_P \cdot V(E) \quad (15) $$

satisfies the following statement:

$$ \forall i = 1, \dots, n_s, \exists \Delta_P / |\delta_{p_k}| \le 1 \text{ and } \mathcal{Y}(\theta_i, \Delta_P) = Y_i. \quad (16) $$

*Proof:* The above proposition is trivially satisfied with the following (nonminimal) choice:

$$ \begin{align*} \Delta_P &= \mathbf{diag}(\delta_{p_1}, \dots, \delta_{p_{2n+2}}) \in \mathbb{R}^{(2n+2) \times (2n+2)} \\ U &= \begin{pmatrix} 1 & \dots & 1 & 0 & \dots & 0 \\ 0 & \dots & 0 & 1 & \dots & 1 \end{pmatrix} \in \mathbb{R}^{2 \times (2n+2)} \end{align*} $$

and $V(E) = \mathbf{diag}(E_{1,1}, \dots, E_{1,n+1}, E_{2,1}, \dots, E_{2,n+1})$. $\square$

Remarking that $\mathcal{Y}(\theta, \Delta_P)$ polynomially depends on $\theta$ and affinely depends on $\Delta_P$, standard algorithms (see [15] for further details) can be applied to compute the interconnection matrix $\mathcal{X}$ such that

$$ \mathcal{Y}(\theta, \Delta_P) = \mathcal{F}_u(\mathcal{X}, \Theta, \Delta_P) = \mathcal{F}_u(\mathcal{X}, \mathbf{diag}(\theta I_{2p}, \Delta_P)). \quad (17) $$

Next, standard LFR object manipulations implemented in the LFR toolbox [15] yield the required open-loop LFR models depicted in (5) and (6). Once again, standard manipulations are used to construct the closed-loop $M(s) - \Delta$ standard forms that will include or not the saturation-type nonlinearity and will be used to check the stability.

### C. Stability Assessment

At this point, a low-order uncertain LFR model $\mathcal{F}_u(P(s), \Delta)$ covering the initial set {$F_i(s)$}$_{i=1...n_s}$ is available. The objective of this section is to prove the stability of the closed-loop LFR model $P(s)$ both with and without input saturation. As summarized in Algorithm 1, the proposed analysis method consists of two steps. No saturation is considered in the first one, which can be viewed as an LFR model validation test. In a second step, an input saturation is introduced and the IQC-based analysis is considered.

1) Stability Analysis Without Saturation Using $\mu$ Tools: Without saturation, the uncertain closed-loop model under consideration assumes that an LTI standard form $M(s) - \Delta$, where $\Delta = \mathbf{diag}(\Theta, \Delta_P, \Delta_R(s))$ is a normalized LTI structured uncertainty block. As a result, the stability of the continuum (covering the initial set of full-order plants) of closed-loop models obtained for any admissible uncertainty inside the unit ball is guaranteed if and only if

$$ \forall \omega \ge 0, \quad \mu_{\Delta}(M(j\omega)) \le 1 \quad (18) $$

where $\mu_{\Delta}(M)$, for any complex-valued matrix $M$, denotes the structured singular value with respect to $\Delta$ and provides the inverse of the size of the critical uncertainty beyond which stability is no longer guaranteed (see [16] for further details). Testing (18) raises two difficulties. The computation of $\mu$ is an NP-hard optimization problem that, in addition, must be solved for an infinite set of frequencies. Yet, as emphasized in [17], recent implementations (used in this paper) of this $\mu$ test in [18] and [19] provide quite efficient tools even for high-order plants with numerous and repeated uncertainties (see also [20]).

**Remark 2:** The proposed $\mu$ test is clearly a necessary stability condition. If there exists $\omega^* \ge 0$ such that $\mu_{\Delta}(M(j\omega^*)) > 1$, then the accuracy of the model should be improved in order to minimize the effects of $\Delta_P$ and $\Delta_R(s)$.
2) *Stability Analysis With Saturation Using IQC*: IQC-based analysis techniques enable us to study a wide range of problems, namely, the robust stability and performance properties of the interconnection $M(s) - \Delta$ of an LTI operator $M(s)$ with a structured model uncertainty $\Delta$ containing nonlinearities, LTI and/or linear time-varying parameters, neglected dynamics, delays, and specific nonlinearities such as friction and hysteresis [21].

Here, standard IQC descriptions are used for both LTI uncertainties, $\Delta$, and sector nonlinearities, denoted by $\varphi$. The originality of our approach resides in the specific algorithm that has been developed to reduce the computational burden. Indeed, standard IQC-oriented analysis methods consist in solving Kalman-Yakubovitch-Popov-based LMI conditions [22]. These standard approaches are, however, intractable for high-order models since the number of scalar optimization variables quadratically increases with the closed-loop order [4]. Moreover, this approach is not compatible with the use of irrational multipliers.<fnref n="4" />

a) **IQC generalities:** An IQC describes a relation between input and output signals of an operator. Since these two formulations are completely equivalent, these constraints can be defined either in the time or the frequency domain. Nevertheless, frequency-domain constraints are often preferred since they lead to simpler stability conditions. The definition of an IQC is given in the frequency domain.

**Definition 1:** Two signals, respectively, of dimension $m$ and $p$, square integrable on $[0, \infty)$, i.e., $v \in L_2^m[0, \infty)$, $w \in L_2^p[0, \infty)$, satisfy the IQC defined by $\Pi: j\mathbb{R} \to \mathbb{C}^{(m+p)\times(m+p)}$, and Hermitian-valued function, iff

$$ \int_{-\infty}^{\infty} \left[ \begin{array}{c} \tilde{v}(j\omega) \\ \tilde{w}(j\omega) \end{array} \right]^* \Pi(j\omega) \left[ \begin{array}{c} \tilde{v}(j\omega) \\ \tilde{w}(j\omega) \end{array} \right] d\omega \geq 0 \quad (19) $$

where $\tilde{v}(j\omega)$ and $\tilde{w}(j\omega)$, respectively, correspond to Fourier transforms of $v$ and $w$, such as $w = \Delta v$.

The problem consists in analyzing the closed loop that corresponds to the interconnection by a positive feedback of $M(s)$ with $\Delta$, where $\Delta$ can be nonlinear and non stationary. Let us suppose that the input and output signals of $\Delta$ satisfy the IQC defined by $\Pi$. The following result gives the stability criterion [22].

**Theorem 1:** Let us suppose that $M(s)$ is stable and $\Delta$ is a causal and bounded operator, if the following conditions hold.

1) Interconnection $M-\tau\Delta$ is well posed for any $\tau \in [0, 1]$.

2) $\tau\Delta$ satisfies the IQCs defined by $\Pi$, $\forall \tau \in [0, 1]$.

3) It exists $\epsilon > 0$ such as

$$ \forall \omega \in \mathbb{R} \underbrace{ \begin{bmatrix} M(j\omega) \\ I \end{bmatrix}^* \Pi(j\omega) \begin{bmatrix} M(j\omega) \\ I \end{bmatrix} }_{Z(j\omega)} \leq -\epsilon I \quad (20) $$

then the closed-loop system is stable.

Let us consider a stable $M(s)$, forming the constant block of the LFR and an augmented block $\Delta \leftarrow \text{diag}(\varphi, \Delta)^5$ where $\varphi$

represents one sector slope-restricted nonlinearity $(0, 1)$. The global multiplier $\Pi$ corresponding to $\Delta$ is built as follows (see [22]–[24] for additional details):

$$ \begin{align} \Pi(j\omega) &= \begin{bmatrix} 0 & 0 & x+j\omega\lambda + \omega^2\gamma & 0 \\ 0 & X(j\omega) & 0 & Y(j\omega) \\ x-j\omega\lambda + \omega^2\gamma & 0 & -2x-2\omega^2\gamma & 0 \\ 0 & Y^*(j\omega) & 0 & -X(j\omega) \end{bmatrix} \nonumber \\ X(j\omega) &= \text{diag}(X_{\Theta}(j\omega), X_P(j\omega), x_{\Delta_R}) \nonumber \\ Y(j\omega) &= \text{diag}(Y_{\Theta}(j\omega), Y_P(j\omega), 0) \tag{21} \end{align} $$

where $X_{\Theta}(j\omega) = X_{\Theta}^*(j\omega) \geq 0 \in \mathbb{C}^{n_{\Theta} \times n_{\Theta}}$, $X_P(j\omega) = X_P^*(j\omega) \geq 0 \in \mathbb{C}^{(2n+2) \times (2n+2)}$, $Y_{\Theta}(j\omega) = -Y_{\Theta}^*(j\omega) \in \mathbb{C}^{n_{\Theta} \times n_{\Theta}}$, $Y_P(j\omega) = -Y_P^*(j\omega) \in \mathbb{C}^{(2n+2) \times (2n+2)}$, $x \geq 0$, $x_{\Delta_R} \geq 0$, $\gamma \geq 0$, and $\lambda \in \mathbb{R}$. Closed-loop stability is ensured if a solution of the following LMI can be found, $\forall \omega \in \mathbb{R}_+$:

$$ \begin{bmatrix} M(j\omega) \\ I \end{bmatrix}^* \Pi(x, \lambda, \gamma, X(j\omega), Y(j\omega)) \begin{bmatrix} M(j\omega) \\ I \end{bmatrix} < 0. \quad (22) $$

b) **Proposed Innovative Method:** In this paper, the optimization problem is directly solved from frequency-domain inequalities through a grid-based approach. A similar approach is used in [25] but without guarantee of the solution validity on the whole frequency domain. Here to guarantee that the solution is valid on the whole frequency domain, a specific technique based on [18] and [26] is adapted to our problem [4]. Besides another advantage is to limit the number of LMI constraints since only active constraints are added in the LMI optimization problem. Here the main result is presented.

Let $\Xi = (A_\Xi, B_\Xi, C_\Xi, D_\Xi)$ be the realization of $\Xi(s)$ (of order $m$), with $\Xi(j\omega) = (I - Z(j\omega))(I + Z(j\omega))^{-1}$ $[(I+Z)$ is invertible], where $Z(j\omega) = Z^*(j\omega)$ is the stability criterion (20), and $\Xi(j(\omega_0 + \delta\omega)) = F_l(S(\omega_0), \delta\omega I_m)$, with $\forall \delta\omega \geq -\omega_0$, i.e., $S(\omega_0)$ is interconnected to $\delta\omega$ as a lower LFT, where $\delta\omega$ is a real parameter. $S(\omega_0)$ reads

$$ S(\omega_0) = \begin{pmatrix} D_\Xi & \frac{C_\Xi}{\sqrt{j}} \\ \frac{B_\Xi}{\sqrt{j}} & -jA_\Xi \end{pmatrix} \star \left( \frac{1}{\omega_0} \begin{pmatrix} I & I \\ -I & -I \end{pmatrix} \right). \quad (23) $$

**Proposition 2:** If $\bar{\sigma}(\Xi(\omega_0)) < 1$, then $\bar{\sigma}(F_l(S(\omega_0), \delta\omega I_m)) < 1$ holds true for $\omega_0 + \delta\omega \in [\underline{\omega}, \bar{\omega}]$, where $\underline{\omega}$ and $\bar{\omega}$ are computed as $\underline{\omega} = \omega_0 + (1/\eta_n)$ and $\bar{\omega} = \omega_0 + (1/\eta_p)$, where $\eta_n$ and $\eta_p$ are the maximal magnitude real negative and positive eigenvalues of $T$, respectively, defined as

$$ T = \begin{bmatrix} S_{22} & 0 \\ 0 & S_{22}^* \end{bmatrix} - \begin{bmatrix} 0 & S_{21} \\ S_{12}^* & 0 \end{bmatrix} X^{-1} \begin{bmatrix} S_{12} & 0 \\ 0 & S_{21}^* \end{bmatrix} \quad (24) $$

where

$$ S(\omega_0) = \begin{bmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{bmatrix} \quad \text{and} \quad X = \begin{bmatrix} I & S_{11} \\ S_{11}^* & I \end{bmatrix}. \quad (25) $$

**Remark 3:** When $\bar{\sigma}(\Xi(+\infty)) = 1$ and $\bar{\omega} = +\infty \Leftrightarrow \eta_p = 0$, a null eigenvalue is obtained, which means that $\bar{\sigma}(\Xi(\omega))$ crosses the 0 dB axis for $\omega = +\infty$. But the intersection of the stability criterion with the 0 dB axis has no physical meaning.

**Remark 4:** The bilinear transformation $\Xi(j\omega) = (I - Z(j\omega))(I + Z(j\omega))^{-1}$ with $(I + Z)$ invertible

<footnote n="4">⁴This constraint imposes to fix a priori (via a time-consuming trial-and-error process) the poles of the multipliers without any guarantee on the optimality of the selected poles.</footnote>

<footnote n="5">⁵Note that $\Delta$ is the same uncertain block as in Section II-B (containing the neglected model reduction dynamics $\Delta_R(s)$, parametric variations $\Theta$, and interpolation errors $\Delta_p$), augmented with $\varphi$, the saturation nonlinearity.</footnote>
**Algorithm 2 Iterative IQC Resolution**

**Data:** $M(j\omega)$ the stable fixed block of the LFR, multiplier $\Pi(\omega_i)$ and $\omega_i \in \mathbb{R}_+$, $i = 1, \dots, n_f$.

**Result:** A stability proof of the LFR model including nonlinear sector saturations.

while stability not checked do
For $i = 1, \dots, n_f$, check the stability criterion
$$ \begin{bmatrix} M(j\omega_i) \\ I \end{bmatrix}^* \Pi(\omega_i) \begin{bmatrix} M(j\omega_i) \\ I \end{bmatrix} < 0. \quad (27) $$

if (27) has solutions then

* Set $\Pi_i \leftarrow \Pi(\omega_i)$ be the solution obtained at $\omega_i$.

* Set $\omega_0 \leftarrow w_i$ and apply Proposition 2.

* For each solution $\Pi_i$, a frequency-domain

$\Omega_i = [\omega_i, \bar{\omega}_i]$ is obtained. $\Omega_{\text{VALID}} = \bigcup_{i=1, \dots, n_f} \Omega_i$.

if $\Omega_{\text{VALID}} = [0 + \infty)$ then

The solution composed by the set of $\Pi_i$ is validated on the whole frequency domain.

Stability is proved, stop.

else

* Determine the complementary set
$\Omega_{\text{NOVALID}} = \mathbb{C}_{[0, +\infty)} \Omega_{\text{VALID}}$.

* Select one or several frequencies in $\Omega_{\text{NOVALID}}$ and update the grid.

else

Stability cannot be proved, stop.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Largest singular value of $G_1$ and of the 16th-order reduced-order model $\hat{G}_1$ obtained with **DARPO** with $\Omega = [0, \omega_r]$. The gray area denotes the discarded frequencies (i.e., above $\omega_r$).</figcaption></figure>

$n_s = 3$ large-scale models ($n_i \approx 600$) representing the local behavior of an industrial aircraft for different Mach numbers, looped with **K**, an antivibration controller ($n_K = 6$) [1].

A. Step 1: LTI Approximation and Error Bound (Section II-B)

1) **Approximation:** The $n_s = 3$ large-scale models $G_i$ of order $n_i \approx 600$ are approximated by $\hat{G}_i$ of order $r = 16$ over $\Omega = [0, \omega_r]$. The frequency interval $\Omega$ is chosen to keep the low-frequency behavior of the large-scale models as it is known to be accurate, whereas the dynamics above $\omega_r$ are less accurately known and are therefore discarded. The approximation order $r$ is then chosen experimentally to achieve a low approximation error over $\Omega$. The relative approximation errors, i.e., $e_i = \|G_i - \hat{G}_i\|_{H_{2,\Omega}} / \|G_i\|_{H_{2,\Omega}}$, $i = 1, 2, 3$, are, respectively, equal to 2.86%, 2.39%, and 2.49%. Fig. 2 illustrates these low errors through the largest singular value of $G_1$ and $\hat{G}_1$.

Fig. 2 illustrates that the dynamics occurring at higher frequencies than $\omega_r$ (gray zone) are indeed discarded during the approximation step. By doing so, one can obtain very accurate reduced-order models over $\Omega = [0, \omega_r]$ as shown by the relative errors that all are below 3%.

The high-frequency dynamics require a complex model to be accurately captured, while the low-frequency ones, which contain the rigid behavior and the first flexible modes of the aircraft, can be caught more easily. This point is particularly obvious when comparing the relative errors obtained here to the one obtained by optimal $H_2$ approximation of the same aircraft model in [1]. In the latter case, with an approximation order $r = 16$, the $H_2$ approximation error is above 30%.

2) **Approximation Error Modeling:** The order of the approximation errors $\Sigma_i(s) = F_i(s) - \hat{F}_i(s) = K(G_i(s) - \hat{G}_i(s))$ prevents from using optimization tools to design the filter $W(s)$ efficiently. That is why it is built here in a heuristic manner. More specifically, $W(s)$ is designed as a product of simple first-order filters $W(s) = k \Pi_{i=1}^{n_w} (s - z_i / s - p_i)$, where the poles $p_i$, zeros $z_i$, and gain $k$ are adjusted for $W(s)$ to be as close as possible of the approximation errors while still ensuring that $\max_{i=1, \dots, n_s} \|W^{-1}\Sigma_i\|_{H_\infty} \le 1$. The filter $W(s)$ obtained here has an order $n_W = 25$ and is plotted in Fig. 3. One can observe that its singular value upper bounds the worst

allows us to transform a positivity condition into a weak gain condition

$$ \bar{\sigma}(\Xi) \le 1 \Leftrightarrow Z + Z^* \ge 0. \quad (27) $$

In the iterative approach, proposed in Algorithm 2, the validation step is a priori performed during the LMI optimization problem resolution. The choice of the initial grid has no influence on the feasibility problem. It is possible to choose a singleton at the first iteration. However, to limit the number of iterations, and consequently the calculation time, without any a priori knowledge, it is recommended to take some frequencies roughly spread on the frequency domain. It is possible, when first solutions are obtained, to tune this initial frequency grid to decrease the number of iterations.

This approach allows the frequency-domain irrational multipliers $X(j\omega)$ to be piecewise continuous. More specifically, between each $\Omega_i$, these multipliers are discontinuous, and consequently, no state-space representation for these multipliers can exist. Involving a state-space representation in order to parameterize multipliers would necessarily lead to constrain the solution and increase the conservatism. Of course it is also possible to use rational multipliers with a frequency-domain resolution using the factorized form of $X(s)$ presented previously [4]. The auxiliary matrix $P$ is still avoided, but without the advantage to use irrational multipliers.

III. APPLICATION TO AN AEROELASTIC AIRCRAFT SYSTEM

The methodology described in Section II and summarized in Algorithm 1 is now applied to check the stability of a set of
<figure><img src="image_4.png" /><figcaption>Fig. 3. Singular values of $W(s)$ and $F_i(s) - \hat{F}_i(s)$ ($i = 1, \dots, n_s$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Illustration of a second-order polynomial interpolation result with minimized guaranteed error bound for the coefficient $c_{19} = \mathcal{Y}_{2,19}$.</figcaption></figure>

approximation error. In particular, with this filter, one obtains $\max_{i=1, \dots, n_s} \|W^{-1}\Sigma_i\|_{\mathcal{H}_{\infty}} = 0.99 < 1$.

## B. Step 2: Interpolation and LFR modeling (Section II-C)

At this stage, a Mach-dependent family $\{\hat{F}_i(s)\}_{i=1...3}$ of 22nd-order LTI models is available together with a common weighting function $W(s)$ shaping the worst case approximation errors induced by the reduction process.

1) Polynomial Approximation With Guaranteed Bounds: The interpolation technique summarized by the linear constraints (13) is initially applied with $p = 2$ and $n_s = 3$. The scalar parameter $\theta$ is normalized in such a way that $\theta = -1$ corresponds to the lowest Mach number of interest, while $\theta = 1$ corresponds to the highest value. Since $n_s - p = 1$, this first trial yields an exact approximation at each of the three interpolation points, but a poor behavior is observed elsewhere. Reducing the order $p$ to 1 would yield a rough and unacceptable approximation. The only remaining option then consists in adding fictitious models for intermediate Mach numbers. This is achieved here by generating additional coefficients in (11) with a standard linear interpolation technique. Two models are then generated for Mach 0.825 and 0.875, and a new interpolation is thus realized with $n_s = 5$ for each of the 46 coefficients contained in the matrices $Y_i$ of (12). A result of this interpolation is plotted in Fig. 4 for one of the most varying coefficients, namely, $\mathcal{Y}_{2,19}(\theta)$. The blue solid line corresponds to the nominal plot, while the red dashed lines visualize lower and upper bounds including the

<figure><img src="image_6.png" /><figcaption>Fig. 5. Visualization of $\mu$ upper and lower bounds for the evaluation of robust stability margins: stability proved for $\|\Delta\|_{\infty} \le 0.43^{-1} = 2.32$.</figcaption></figure>

five interpolation points. Note that the three coefficients from the initial set of models are all located on the same bound (the upper bound for this coefficient). Quite interestingly, this property holds true for the 46 [$=2 \times (r + n_K + 1) = 2 \times (n + 1)$] coefficients, which permits to drastically reduce the size of $\Delta_P$ in (14). One here obtains $\Delta_P = \delta_P I_2$, and (15) boils down to

$$ \mathcal{Y}(\theta, \Delta_P) = X_0 + \theta X_1 + \theta^2 X_2 + \text{diag}(\delta_P, \delta_P) \cdot V(E). \quad (28) $$

2) LFR Modeling: As already clarified in Section II, $\mathcal{Y}(\theta, \Delta_P)$ is readily rewritten in an LFR format with the help of existing software [15]. Next, exposed in (6), a global 47th-order ($=r+n_W+n_K$) dynamic LFR model encompassing the whole initial set of full-order open-loop plants is obtained. The structure of its $7 \times 7$ $\Delta$-block reads

$$ \Delta = \text{diag}(\theta I_4, \delta_P I_2, \Delta_R(s)) \quad (29) $$

and has a **minimal** size that remains largely compatible with the specific $\mu$ and IQC-based analysis tools to be applied next.

## C. Step 3: Stability Analysis (Section II-D)

1) Preliminary Tests via $\mu$ Analysis: As mentioned in Section II-C, the validity of the global LFR model is preliminarily checked without saturation. An uncertain LTI closed-loop model is then built, and the $\mu$ analysis test (18) is performed. Since the complexity of our algorithm is not directly impacted by the number of states, but mainly depends on the size and structure of $\Delta$, the results are obtained in a few seconds on any standard computer. A guaranteed upper bound of $\mu$ as a function of frequency is displayed in Fig. 5. The yellow stars corresponding to lower bounds reveal a rather low conservatism of our test, which can be summarized by

$$ \sup_{\omega \ge 0} \mu_{\Delta}(\mathbf{M}(j\omega)) = 0.43 \ll 1. \quad (30) $$

The continuum of closed-loop models for any admissible uncertainty then clearly remains stable, which concludes the preliminary validation phase.

2) Stability Assessment via IQC-Based Analysis: An input saturation converted to a deadzone operator $\phi$ is now inserted in the uncertain closed-loop whose $\Delta$-block is then augmented: $\Delta \leftarrow \text{diag}(\phi, \Delta)$. The initial frequency grid is $\omega_i = \{1, 5, 10, 20, 100\}$ rad/s with $i = 1, \dots, 5$. To limit the number of decision variables and then the computation time, $X_{\Theta}(j\omega)$ and $Y_{\Theta}(j\omega)$ are chosen diagonal, which leads
to 17 scalar decision variables for each frequency, even if it is
possible to use the general form if no solution was obtained.
In addition, three decision variables x, λ, and γ come from
the multiplier, which corresponds to the static nonlinearity φ.
A solution has been obtained in 8 iterations and
104 frequencies. The total number of decision variables
is $17 \times 104 + 3 = 1771$. The following remarks can be made.

1) The solution $X(j\omega)$ is a positive, complex, constant, and piecewise continuous 6 × 6 matrix. For example, at iteration 8, for $\omega_3 = 10$ rad/s, the solution $\Pi_3(j\omega)$ is valid on the frequency domain $\Omega_3 = [9.72, 32.82]$ rad/s. Finally, after eight iterations, $\Omega_{\text{VALID}} = \bigcup_{i=1,...,104} \Omega_i = [0, \infty)$, and consequently, the solution is validated on the whole frequency domain.

2) An a priori trial-and-error approach to determine the parameterization for multipliers is not required here. Furthermore, with rational multipliers, if no solution is obtained with a specific parameterization, it is still impossible to conclude on the feasibility problem as a different or more complex parameterization may have enabled to find a solution. Both points highlight the methodological superiority of irrational multipliers that can be considered only from a frequency domain point of view.

3) Finally, the stability of the uncertain and nonlinear closed loop is proved on the large-scale dynamical model.

IV. CONCLUSION AND PERSPECTIVES

In this paper, a methodology that enables us to assess the
stability of a set of controlled SIMO large-scale LTI dynamical
models subject to input saturation has been presented. First,
the large-scale models are reduced and interpolated, and the
associated errors are bounded. This leads to a small-scale LFR
that represents both the parametric variation of the initial set
of models and the errors induced during the reduction and
interpolation steps. The stability analysis is then achieved with
an innovative algorithmic approach based on IQC techniques.
Unlike standard methods that require a possibly conservative
parameterization of the multiplier, here, no parameterization
is required. This decrease in the conservatism enables us
to drastically improve the approach. The methodology is
successfully validated on an industrial set of controlled large-
scale aircraft models subject to saturation limitations.<fnref n="6" />

The extension of the methodology to MIMO models is
conditioned by the use of an interpolation technique with
guaranteed error bounds. The development of such a technique
is still under investigation.

REFERENCES

[1] C. Poussot-Vassal, C. Roos, T. Loquen, P. Vuillemin, O. Cantinaud, and J.-P. Lacoste, "Control-oriented aeroelastic BizJet low-order LFT modeling," in *Control-Oriented Modelling and Identification: Theory and Practice*, M. Lovera, Ed. Edison, NJ, USA: IET, 2014, ch. 11, pp. 241–268.

[2] R. Ionutiu, J. Rommes, and W. H. A. Schilders, "SparseRC: Sparsity preserving model reduction for RC circuits with many terminals," *IEEE Trans. Comput.-Aided Design Integr. Circuits Syst.*, vol. 30, no. 12, pp. 1828-1841, Dec. 2011.

[3] P. Vuillemin, F. Demourant, J.-M. Biannic, and C. Poussot-Vassal, "Global stability validation of an uncertain large-scale aircraft model," in *Proc. IEEE Conf. Control Appl.*, Oct. 2014, pp. 152-157.

[4] F. Demourant, "New algorithmic approach based on integral quadratic constraints for stability analysis of high order models," in *Proc. Eur. Control Conf.*, Jul. 2013, pp. 359-364.

[5] D. Petersson and J. Löfberg, "Model reduction using a frequency-limited H₂-cost," *Syst. & Control Lett.*, vol. 67, pp. 32-39, May 2014.

[6] P. Vuillemin, C. Poussot-Vassal, and D. Alazard, "Poles residues descent algorithm for optimal frequency-limited H₂ model approximation," in *Proc. Eur. Control Conf.*, 2014, pp. 1080-1085.

[7] A. C. Antoulas, *Approximation of Large-Scale Dynamical Systems* (Advances in Design and Control). Philadelphia, PA, USA: SIAM, 2005.

[8] P. Fulcheri and M. Olivi, "Matrix rational $H_2$ approximation: A gradient algorithm based on Schur analysis," *SIAM J. Control Optim.*, vol. 36, no. 6, pp. 2103-2127, 1998.

[9] S. Gugercin, A. C. Antoulas, and C. Beattie, "$H_2$ model reduction for large-scale linear dynamical systems," *SIAM J. Matrix Anal. Appl.*, vol. 30, no. 2, pp. 609-638, 2008.

[10] P. Vuillemin, C. Poussot-Vassal, and D. Alazard, "Spectral expression for the frequency-limited $H_2$-norm of LTI dynamical systems with high order poles," in *Proc. Eur. Control Conf.*, Jun. 2014, pp. 55-60.

[11] W. K. Gawronski, *Advanced Structural Dynamics and Active Control of Structures*. Germany: Springer Mechanical Engineering, 2004.

[12] M. Lovera, C. Novara, P. L. dos Santos, and D. Rivera, "Guest editorial special issue on applied LPV modeling and identification," *IEEE Trans. Control Syst. Technol.*, vol. 19, no. 1, pp. 1-4, Jan. 2011.

[13] J. De Caigny, R. Pintelon, J. F. Camino, and J. Swevers, "Interpolated modeling of LPV systems," *IEEE Trans. Control Syst. Technol.*, vol. 22, no. 6, pp. 2232-2246, Nov. 2014.

[14] G. Ferreres, "Computation of a flexible aircraft LPV/LFT model using interpolation," *IEEE Trans. Control Syst. Technol.*, vol. 19, no. 1, pp. 132-139, Jan. 2011.

[15] J.-F. Magni. (2006). "Linear fractional representation toolbox for use with MATLAB," French Aerosp. Lab, ONERA, Châtillon, France, Tech. Rep. TR 5/10403.01F DCSD, Oct. 2005. [Online]. Available: http://w3.onera.fr/smac/

[16] A. Packard and J. Doyle, "The complex structured singular value," *Automatica*, vol. 29, no. 1, pp. 71-109, 1993.

[17] C. Roos and J.-M. Biannic, "A detailed comparative analysis of all practical algorithms to compute lower bounds on the structured singular value," *Control Eng. Pract.*, vol. 44, pp. 219-230, Nov. 2015.

[18] J.-M. Biannic and G. Ferreres, "Efficient computation of a guaranteed robustness margin," in *Proc. IFAC World Congr.*, 2005, pp. 1-6.

[19] C. Roos, F. Lescher, J.-M. Biannic, C. Döll, and G. Ferreres, "A set of μ-analysis based tools to evaluate the robustness properties of high-dimensional uncertain systems," in *Proc. IEEE Int. Symp. Comput.-Aided Control Syst. Design (CACSD)*, Sep. 2011, pp. 644-649.

[20] C. T. Lawrence, A. L. Tits, and P. Van Dooren, "A fast algorithm for the computation of an upper bound on the μ-norm," *Automatica*, vol. 36, no. 3, pp. 449-456, Mar. 2000.

[21] H. Pfeifer and P. Seiler, "Robustness analysis of linear parameter varying systems using integral quadratic constraints," *Int. J. Robust Nonlinear Control*, vol. 25, no. 15, pp. 2843-2864, 2015.

[22] A. Megretski and A. Rantzer, "System analysis via integral quadratic constraints," *IEEE Trans. Autom. Control*, vol. 42, no. 6, pp. 819-830, Jun. 1997.

[23] U. T. Jonsson and A. Rantzer, "A unifying format for multiplier optimization," in *Proc. Amer. Control Conf.*, vol. 5. 1995, pp. 3859-3860.

[24] P. Park, "Stability criteria of sector- and slope-restricted Lur'e systems," *IEEE Trans. Autom. Control*, vol. 47, no. 2, pp. 308-313, Feb. 2002.

[25] M. S. Andersen, S. K. Pakazad, A. Hansson, and A. Rantzer, "Robust stability analysis of sparsely interconnected uncertain systems," *IEEE Trans. Autom. Control*, vol. 59, no. 8, pp. 2151-2156, Aug. 2014.

[26] A. Sideris and R. S. S. Pena, "Robustness margin calculation with dynamic and real parametric uncertainty," *IEEE Trans. Autom. Control*, vol. 35, no. 8, pp. 970-974, Aug. 1990.

<footnote n="6">⁶Note that the reduction, interpolation, and stability analysis tools used in this paper are gathered in the MORE Toolbox and in the SMAC Toolbox available from w3.onera.fr/more and from w3.onera.fr/smac, respectively.</footnote>