Stability Analysis of a Set of Uncertain Large-Scale Dynamical Models With Saturations: Application to an Aircraft System

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 irrational 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.


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 applications are 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 1 , 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 [22], satellites, launchers, fluid flow mechanics), civilian structures, electronics (e.g., [11]), where control engineers have to cope with many practical problems, including lightly damped modes, nonlinear actuator(s), etc. 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 of a physical dynamical system, which smoothly depends on a parameter p θ ∈  . This model is assumed to be only known through its linearized models G i at some parametric points i θ ( 1, , s i n =  ). Let G i be asymptotically stable large-scale Linear Time Invariant (LTI) dynamical models given by the state-space realizations: 11 12 ( ) ( ) ( ) 2 lin. 21 22 : 1 Of course, given that every model can always be questioned or amended, the approach is valid only according to the considered dynamical models, and additional precautions should be considered when it is applied to the real system.

x t B w t B u t z t C x t D w t D u t y t C x t D w t D u t
where ( ) i n i x t ∈  , ( ) w n w t ∈  , ( ) u t ∈ , ( ) z n z t ∈  and ( ) y n y t ∈  are the states, exogenous input, single control input, performance output and measurement signals, respectively. Moreover, let be given a robust th k n order LTI controller = ( , , , ) − + , looped between ( ) y t and ( ) u t , that ensures some robustness and performance specification(s) for all of the s n models. Such a controller could, for instance, be obtained with robust optimization tools, such as [3]. For an example of synthesis, see [21] and the references therein.
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 the paper and pursuant to Figure 1 and Algorithm 1, a three-step methodology is proposed: (i) approximate the s n dynamical models and bound the mismatch error, (ii) perform (inexact) interpolation of the reduced-order models with interpolation error bounds and, finally, (iii) assess the stability of the closed-loop model over both parametric variations and control input saturation limitations 2 . 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.
In comparison to [22] and [26] contributions, the proposed approach is accompanied with both approximation (Step (i)) and interpolation (Step (ii)) errors. Hence, the µ (structured singular value) and Integral Quadratic Constraint (IQC) analysis (Step (iii)) respectively provide sufficient stability conditions for the entire 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 multiplierbased approach developed in Step (iii) is an extension of [6]. It is shown that no solution is obtained by means of a rational multiplier and only a frequency domain approach can be used here to assess the closed-loop stability.
The paper is organized following the schematic view of Figure 1. First, 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. Then we illustrate the proposed procedure on a complex large-scale aeroservoelastic business jet aircraft model for various flight configurations, looped with an anti-vibration controller. To end, Conclusions are given.

Main result: Stability guarantee of a set of large-scale models subjec t to input saturations
With reference to Figure 1, the proposed contribution, in three steps, are 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 (i), Section "Multi-LTI model approximation and error bound"), then the interpolation and transformation into a Linear Fractional Representation (LFR) structure is achieved (Step (ii), Section "Bounded-error reduced-order LFR model generation"), and finally, the stability of the overall uncertain, parameter-dependent model is firstly assessed thanks to a µ analysis, and then, when subject to control input saturation, through a novel IQC technique (Step (iii), Section "Stability assessment").
Algorithm 1 -Global procedure • Close the open-loop LFR model ( ) P s with input saturation to obtain an augmented nonlinear standard form ( ) ( , ) M s ϕ − diag ∆ and check the robust stability by means of an IQC-based analysis test.
return A stability proof of the input-saturated closed-loop large-scale models.

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 [2] for a general overview of model reduction and refer to Box 1 for an overview of the tool used here to perform the model approximation step) 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 H 2 -norm has often been considered and several methods are now available to address the corresponding optimal H 2 model approximation problem (see e.g., [8,10]). However, in many cases, considering a limited frequency interval only is more relevant since (i) 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 the approximation accuracy to be increased, where the initial model is accurate. Besides (ii), controllers are usually designed to act over a limited frequency interval (due to actuator bandwidth or to prevent them from disturbing non-modeled dynamics), which means that a precise knowledge of the dynamics over the entire frequency domain is not necessarily useful. From the authors' point of view, the optimal approximation over a bounded frequency interval enables these practical considerations to be translated elegantly and is therefore preferred here. It is addressed through the use of the frequency-limited H 2 -norm in Section "Optimal frequency-limited H 2 model approximation". However, it is worth noticing that the overall methodology summarized in Algorithm 1 model reduction toolbox

Box 1 -The MORE toolbox
The more toolbox gathers a set of tools aimed at alleviating the numerical burden induced by the complexity of dynamical models (e.g. for simulation, control, optimization, etc.).
More specifically, it contains several model approximation techniques designed to cope with several large-scale problems as depicted below.
More formally, the problems that can be adressed are the following: • Reduction from state-space: considering a LTI dynamical model H represented by a large-scale differential equation, where ( ) x t ∈ , ( ) u n u t ∈  and ( ) y n y t ∈  are the state, command inputs and outputs of the model, respectively. The objective is to find a smaller model Ĥ represented by x t r n ∈   and ˆ( ) y n y t ∈  such that the input-output behaviors of H and Ĥ are close.
In the toolbox, this closeness is generally considered through optimality considerations based on the H 2 -norm of the approximation error or its restriction to a bounded frequency interval (as used in this paper).
• Reduction from data: the initial model is only known through a set of frequency data The objective is then to find a low-complexity model such as Ĥ in equation (B1-2) that matches the frequency data.
• Reduction of infinite dimensional models: the initial model is known through its irrational transfer function ( ) y u n n H s × ∈  obtained for instance from a partial differential equation (PDE), from a delayed differential equation, etc. Again, the objective is to build a lowcomplexity model Ĥ as in (B1-2) such that the input-output behavior of H is well reproduced (for instance in the H 2 sense).
For further information, interested readers should refer to the site of the toolbox : www.onera.fr/more.
with r n 

(+) Simulation (+) Analysis (+) Control (+) Optimization
Infinite order equations (require meshing) does not depend on the approximation strategy, since the approximation error is bounded in Section "Bound on the approximation error".

Optimal frequency-limited model approximation
Using the 2, Ω H -norm, one can formulate the approximation over a bounded frequency interval as an optimization problem. More specifically, given an asymptotically stable n th order large-scale model G and a frequency interval Ω, the optimal 2, Ω H model approximation problem consists in finding a reduced-order model Ĝ of order r n  that minimizes the 2, Ω H -norm of the approximation error − G G, i.e., 2, , Here, Problem (9) is addressed using the method called Descent Algorithm for Residue and Pole Optimization (DARPO), proposed in [27]. It relies on the pole-residue formulation of the 2, Ω H -norm [28] and finds the poles and associated residues of the reduced-order model that satisfy the first-order optimality conditions associated with Problem (9). Note that, since this problem is not convex, the reducedorder model obtained this way is only a local minimum.
With reference to Algorithm 1 (Step (i)), the approximation algorithm is applied to each large-scale model i G , = 1, , s i n  resulting in n s small-scale models ˆi G minimizing the 2, Ω H -norm of the approximation error with i G , as stated in (2).
Note that both the approximation order r and the frequency-interval Ω are tuning parameters that depend on the considered application. However, as mentioned before, the frequency interval Ω can 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 Ω Ω   (see e.g., [9,Chap. 4]), which can be viewed as the Hankel singular values in the frequencylimited 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.

Bound on the approximation error
Then, the stability of the set of uncertain models can be used, since one can always exhibit The controller is included here to be consistent with the interpolation step of Section "Bounded-error reduced-order LFR model generation".
The design of ( ) W s then consists in a trade-off between complexity and conservatism. Indeed, one must find a ( ) W s that is both an accurate modeling of the worst approximation error and whose complexity (order) is reasonable. For instance, (3). However, it does not offer an accurate model of the approximation error and might, therefore, be too conservative for stability analysis. A direct approach to design ( ) W s satisfying (10) would consist in using non-smooth H ∞ optimization tools [3] to solve the following problem However, depending on the application, the errors i Σ might be too large for such an approach to be tractable. In those cases, a heuristic approach may then be preferable.

Bounded-error reduced-order LFR model generation
Consider the parametrically-dependent set i n F s  of reducedorder models obtained above; the objective is now to derive a limited-size LFR, such that µ and IQC-based analysis tools can then be applied. In the general case, involving several parameters ( p θ ∈  ), the n s equations (4) must be solved for a parametric structure, e.g., should be kept as small as possible. Efficient solutions, based on multivariate sparse polynomial or rational interpolation techniques, are detailed in [14,5,22].
In the case of a scalar parameter (θ ∈ ), a specific technique can be developed to compute low-order LFR models whose ∆-block will both include the parametric variations (Θ) and a normalized real-valued uncertain operator (∆ P ). The latter is introduced to "cover" the interpolation errors, as illustrated by Equation (4). The proposed technique, based on a polynomial state-space data interpolation approach, can be broken down into three steps, which are briefly presented next.

Step 1: model 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 thus a good choice, but unfortunately leads to ill-conditioned matrices as the system order increases. As is also proposed in [7], a rescaled companion form will then be used. Using the notation + the system is rewritten as: where the scaling variables

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 p th order polynomial, the size of = n I θ θ Θ will then be limited to = 2 n p Θ .

Polynomial interpolation with guaranteed error bounds
Let us denote by i Y the last two lines in Equation (12): and focus on the polynomial approximation of the finite set ): 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 s n p − gets too small. This issue and possible remedies are further discussed in the applicative part.

Proposition 1
From Inequalities (14), E-dependent "shaping" matrices U (E ) and V (E ) of appropriate dimensions and a bounded, real-valued, blockdiagonal uncertain operator P ∆ : ( )  (15) can be easily defined, such that the function: satisfies the following statement:

Proof
The above proposition is trivially satisfied with the following (nonminimal) choice: ( ) polynomially depends on θ and affinely depends on P ∆ , standard algorithms (see [15] for further details) can be applied to compute the interconnection matrix  , such that: Next, standard LFR object manipulations implemented in the LFR toolbox [15] yield the required open-loop LFR models depicted in (4) and (5). Once again, standard manipulations are used to "construct" the closed-loop ( ) M s − ∆ standard forms that will include or not the saturation-type nonlinearity and will be used to check the stability.

Stability assessment
At this point, a low-order uncertain LFR model ( ( ), ) u P s ∆  covering the initial set 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, which can be viewed as a LFR model validation test. In a second step, an input saturation is introduced and the IQC-based analysis is considered.

Stability analysis without saturation using µ tools
Without saturation, the uncertain closed-loop model under consideration assumes an LTI standard form 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: where ( ) M µ ∆ , for any complex-valued matrix M, denotes the structured singular value with respect to ∆ and provides the inverse of the size of the critical uncertainty beyond which stability is no longer guaranteed (see [17] for further details). Testing (19) raises two difficulties. The computation of µ is an NP-hard optimization problem, which, in addition, must be solved for an infinite set of frequencies.
However, as is emphasized in [23], recent implementations (used in this paper) of this µ test in [4,24] provide quite efficient tools even for high-order plants with numerous and repeated uncertainties (see also [13]).

Remark 2
The proposed µ test is clearly a necessary stability condition. If there exists 0 , then the accuracy of the model should be improved in order to minimize the effects of P ∆ and ( ) R s ∆ .

Stability analysis with saturation using IQC
IQC-based analysis techniques enable a wide range of problems to be studied, namely the robust stability and performance properties of the interconnection ( ) M s − ∆ of an LTI operator ( ) M s with a structured model uncertainty ∆ containing nonlinearities, LTI and/or linear time-varying (LTV) parameters, neglected dynamics, delays, specific nonlinearities such as friction, hysteresis, etc. (see, e.g., [20]).
Here, standard IQC descriptions are used for both LTI uncertainties, ∆ and sector nonlinearities, denoted by ϕ. 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 KYP (Kalman-Yakubovitch-Popov)-based LMI conditions [16]. Theses standard approaches are however intractable for high-order models, since the number of scalar optimization variables quadratically increases with the closed-loop order [6]. Moreover, this approach is not compatible with the use of irrational multipliers 4 .

IQC generalities
An IQC describes a relation between the 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, ) ∞ , i.e. : The problem consists in analyzing the closed-loop that corresponds to the interconnection by a positive feedback of ( ) M s with ∆, where ∆ can be nonlinear and non-stationary. Let us suppose that input and output signals of ∆ satisfy the IQC defined by Π. The following result gives the stability criterion [16].
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 where ϕ represents one sector slope-restricted nonlinearity (0,1). The global multiplier Π corresponding to ∆ is built as follows (see [12,16,18] for additional details): 4 This constraint renders it necessary to fix the poles of the multipliers a priori (via a time-consuming trial-and-error process), without any guarantee on the optimality of the selected poles. 5 Note that ∆ is the same uncertain block as in Section "Bounded-error reducedorder LFR model generation" (containing the neglected model reduction dynamics ( ) R s s ∆ , parametric variations Θ and interpolation errors ∆ P ), augmented with ϕ, the saturation nonlinearity.
γ ≥ and λ ∈ . Closed-loop stability is ensured if a solution of the following LMI can be found, ω + ∀ ∈ :

Proposed innovative method
In this paper, the optimization problem is solved directly from frequency domain inequalities through a grid-based approach. A similar approach is used in [1], but without guarantee of the solution validity over the entire frequency domain. Here, in order to guarantee that the solution is valid over the entire frequency domain, a specific technique based on [25,4] is adapted to our problem [6]. In addition, 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.

Remark 4
The bilinear transformation invertible allows a positivity condition to be transformed into a weak gain condition: In the iterative approach, proposed in Algorithm 2, the validation step is performed a priori and 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, in order 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 throughout 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ω to be piecewise continuous. More specifically, between each i Ω , these multipliers are discontinuous, consequently no state-space representation for these multipliers can exist. Involving a state-space representation in order to parameterize multipliers would necessarily lead to constraining the solution and increasing the conservatism. Of course, it is also possible to use rational multipliers with a frequency domain resolution, by using the factorized form of ( ) X s presented previously [6]. The auxiliary matrix P is still avoided, but without the advantage of using irrational multipliers.

Application to an aeroelastic aircraft system
The methodology described in Section "Main result: Stability guarantee of a set of large-scale models subject to input saturations" and summarized in Algorithm 2 is now applied to check the stability of a set of n s = 3 large-scale models ( i n ≈ 600) representing the local behavior of an industrial aircraft for different Mach numbers, looped with K, an anti-vibration controller (n K = 6) [22].

Approximation
The n s = 3 large-scale models G i of order i n ≈ 600, are approximated by ˆi G of order r = 16 over The frequency interval Ω is chosen to keep the low frequency behavior of the large-scale models, since it is known to be accurate, whereas the dynamics above r ω are less accurately known and are therefore discarded.
The approximation order r is then chosen experimentally to achieve a low approximation error over Ω. The relative approximation errors, i.e., , i = 1, 2, 3, are respectively equal to 2.86 %, 2.39 % and 2.49 %. Figure 2 illustrates these low errors through the largest singular value of G 1 and of 1 G . Figure 2 illustrates that the dynamics occurring at higher frequencies than r ω (gray zone) are indeed discarded during the approximation step. By doing so, one can obtain very accurate reduced-order models over [0, ] r ω Ω = , as shown by the relative errors, which are all 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 captured more easily. This point is particularly obvious when comparing the relative errors obtained here to that obtained by optimal H 2 approximation of the same aircraft model in [22]. In the latter case, with an approximation order r = 16, the H 2 approximation error is above 30 %.
if (28)  An alternative illustration of the relevance of the frequency-limited approach in comparison to the standard H 2 approach is presented in Figure 4. The time responses of the approximation errors between the first input-output transfers of H 1 and 1 G and an optimal H 2 reducedorder model of order 16 for a sinusoidal input of frequency below and above r ω are shown. One can see that the frequency-limited approach leads to a significantly lower error when the input signal acts below r ω (left plot in Figure 4), while the H 2 approach is more efficient outside of the frequency interval (right plot).

Approximation error modeling
The order of the approximation errors ( The filter ( ) W s obtained here has an order n w = 25 and is plotted in Figure 3. One can observe that its singular value upper bounds the worst approximation error. In particular, with this filter, Step 2: Interpolation and LFR modeling (II-C) At this stage, a Mach-dependent family =1 3 { ( )} i i F s  of 22 nd order LTI models is available, together with a common weighting function ( ) W s shaping the worst-case approximation errors induced by the reduction process.

Polynomial approximation with guaranteed bounds
The interpolation technique summarized by the linear constraints (14) is initially applied with p = 2 and n s = 3. The scalar parameter θ is normalized in such a way that θ = -1 corresponds to the lowest Mach number of interest, while θ = 1 corresponds to the highest value. Since n sp = 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 (12), 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 i Y of (6). A result of this interpolation is plotted in Figure 5 for one of the most varying coefficient, namely 2,19 ( ) plot, while the dashed red lines represent lower and upper bounds, including the 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 ( 1) 2 ( 1) K r n n = × + + = × + ) coefficients, which permits the size of P ∆ to be reduced drastically in (15). Here, one obtains

LFR modeling
As has already been clarified in Section "Main result: Stability guarantee of a set of large-scale models subject to input saturations", ( , ) P θ ∆ Y is readily rewritten in a LFR format with the help of existing software [15]. Next, exposed in Equation (2) and has a minimal size that remains largely compatible with the specific µ and IQC based analysis tools to be applied next.

Preliminary tests via µ analysis
As mentioned in Subsection "Stability assessment", the validity of the global LFR model is preliminarily checked without saturation. An uncertain LTI closed-loop model is then built and the µ analysis test (19) is performed. Since the complexity of our algorithm is not directly affected by the number of states, but mainly depends on the size and structure of ∆, the results are obtained in a few seconds on any standard computer.
A guaranteed upper-bound of µ as a function of frequency is displayed in Figure 6. The yellow stars corresponding to lower-bounds reveal a rather low conservatism of our test, which can be summarized by: The continuum of closed-loop models, for any admissible uncertainty, then clearly remains stable, which concludes the preliminary validation phase.

Stability assessment via IQC-based analysis
An input saturation -converted to a deadzone operator ϕ, is now inserted in the uncertain closed-loop whose ∆-block is then augmented: . The initial frequency grid is = {1,5,10, 20,100} i ω rad/s with = 1,...,5 i . To limit the number of decision variables and then the computation time, ( ) X jω Θ and ( ) Y jω Θ are chosen to be diagonal, which leads to 17 scalar decision variables for each frequency, even though it is possible to use the general form if no solution was obtained. In addition, 3 decision variables , , x λ γ 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 × 104 + 3 = 1771. The following remarks can be made: • The solution ( ) X jω is a positive, complex, constant and piecewise continuous 6 × 6 matrix. For example, at iteration 8, for • 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, since a different or more complex parameterization may have enabled a solution to be found. Both points highlight the methodological superiority of irrational multipliers, which can only be considered from a frequency domain point of view.
• Finally, the stability of the uncertain and nonlinear closed-loop is proved on the large-scale dynamical model.

Conclusion and perspectives
In this paper, a methodology enabling the stability of a set of controlled SIMO large-scale LTI dynamical models subject to input saturation to be assessed has been presented. Firstly, the large-scale models are reduced, interpolated and the associated errors are bounded. This leads to a small-scale LFR, which 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 the approach to be drastically improved. The methodology is successfully validated on an industrial set of controlled largescale aircraft models subject to saturation limitations. 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. Similarly, determining whether the methodology can easily be extended to a broader class of models (e.g., descriptor models) requires further studies 