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

HAL
open science

# Robust stability and performance analysis in the presence of LTI and LTV uncertainties with bounded rates of variation

Tommaso Casati, Clément Roos, Jean-Marc Biannic, Hélène Evain

► To cite this version:

Tommaso Casati, Clément Roos, Jean-Marc Biannic, Hélène Evain. Robust stability and performance analysis in the presence of LTI and LTV uncertainties with bounded rates of variation. European Journal of Control, 2025, 10.1016/j.ejcon.2025.101298. hal-05228616

HAL Id: hal-05228616

https://hal.science/hal-05228616v1

Submitted on 28 Aug 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Robust stability and performance analysis in the presence of LTI and LTV uncertainties with bounded rates of variation

Tommaso Casati<sup>a, ID, *</sup>, Clément Roos<sup>a, ID</sup>, Jean-Marc Biannic<sup>a, ID</sup>, Hélène Evain<sup>b</sup>

## ARTICLE INFO

Recommended by T. Parisini

**Keywords:**

Uncertain systems

IQCs

LTV uncertainties with bounded rates of variation

Robust stability

Robust $\ell_2$-gain performance

Space applications

## ABSTRACT

The assessment of robust stability and performance of systems affected by Linear Time-Invariant (LTI) and Linear Time-Varying (LTV) uncertainties with bounded rates of variation is still an open problem and yet of importance. To this purpose, Integral Quadratic Constraints (IQCs) can be used to derive sufficient conditions in terms of Frequency Domain Inequalities (FDIs) involving both dynamic and static multipliers. In this paper, a new fully frequency-based approach is proposed which offers two main advantages. Firstly, LTI uncertainties are characterized by the independent frequency-by-frequency multipliers of $\mu$-analysis, which are much less conservative than those used in IQC approaches, for which a basis must be chosen arbitrarily. Secondly, sub-optimal (but relevant) solutions can be obtained by iteratively optimizing the static and frequency-dependent multipliers. This lightens the computational load and makes the method more applicable to real systems.

## 1. Introduction

The need for increasingly autonomous space missions with reduced development times requires new Verification and Validation (V&V) techniques to assess performance in the presence of uncertainties and external disturbances for Attitude and Orbit Control Systems (AOCS) (Dennehy et al., 2022). Among the explored validation methods, those based on $\mu$-analysis have produced very promising results in the presence of many Linear Time-Invariant (LTI) parametric uncertainties (Biannic et al., 2021; Roos, 2018), but they allow only a partial validation of the control laws because of the required assumptions.

A general approach to treat systems including non-linearities, time variations and uncertainties has been developed by Megretski and Rantzer (1997). This method describes the characteristics of the perturbations by means of Integral Quadratic Constraints (IQCs) and formulates analysis tests in terms of Frequency Domain Inequalities (FDIs) which have to be satisfied by the plant. As the obtained FDIs require solving an infinite number of Semi-Definite Programs (SDPs), the Kalman-Yakubovich-Popov (KYP) lemma is generally used to transform them into state-space Linear Matrix Inequalities (LMIs). In this way, robust stability or performance can be assessed with only one SDP, but the dimension of the problem rapidly grows with the number of states. This fact is even more evident by considering that the initial state-space representation of the plant shall be augmented with the state-space representations of the dynamic multipliers, which are introduced to reduce the conservatism of the results. In this regard, it

is well-known from $\mu$-analysis that frequency-dependent DG scalings may be adopted when dealing with LTI uncertainties, whereas constant DG multipliers are used for arbitrarily-fast Linear Time-Varying (LTV) parameters. To fill the gap between these two extremes, peculiar multipliers have also been investigated in the literature to describe LTV uncertainties with bounded rates of variation (Helmersson, 1999; Jönsson & Rantzer, 1996). More recently, Köroğlu and Scherer proposed IQC tests based on a generalized version of the swapping lemma adopted by Helmersson (1999), to perform robust stability and $\ell_2$-gain performance analysis (Köroğlu & Scherer, 2007). In practice, the IQC-based theoretical results have been implemented in toolboxes such as IQC-$\beta$ (Kao, Megretski, Jönsson, & Rantzer, 2004), the Sedumi-based IQC solver of SMAC (Biannic et al., 2016) and the most recent IQClab (Veenman et al., 2021), which all use the KYP lemma to derive state-space LMIs from FDIs. An interesting approach was also proposed by Demourant, who used a frequency-domain IQC-based method to perform robust stability analysis of a fighter aircraft model (Demourant, 2017). The problem is first solved on a finite frequency grid, and a Hamiltonian-based technique inspired by Ferreres and Biannic (1998) is then adopted to validate the solution on the entire frequency range. In this case the number of decision variables depends only on the number of uncertainties, their occurrences in the model and the order of the selected scalings. The IQC parameterization for LTV uncertainties, however, is the one of Jönsson and Rantzer (1996), which was proved conservative by Helmersson (1999).

<footnote><sup>a</sup> DTIS, ONERA, Université de Toulouse, Toulouse, 31000, France</footnote>

<footnote><sup>b</sup> CNES - The French Space Agency, Toulouse, 31400, France</footnote>

<footnote>* Corresponding author.
E-mail address: tommaso.casati@onera.fr (T. Casati).</footnote>
The present paper advances the research on the topic by proposing a new approach to perform robust stability and performance analysis of uncertain systems characterized by both LTI and LTV uncertainties with bounded rates of variation. In particular, two main achievements with respect to existing literature are presented. First, following the approach of Demourant (2017), the stability and performance IQC-based tests of Megretski and Rantzer (1997) are assessed in the frequency domain on a finite grid and a Hamiltonian-based technique is then used to validate the results on a continuous frequency range. With respect to Demourant (2017), however, the multipliers of Köroğlu and Scherer (2007) are used to account for LTV uncertainties with bounded derivatives, which were proved less conservative than the ones of Jönsson and Rantzer (1996). Secondly, the multipliers used in the paper to describe LTI uncertainties correspond to the well-known DG scaling matrices of µ-analysis, which are independent frequency by frequency. Such scalings do not require a state-space representation and are therefore less conservative than the multipliers currently adopted in the context of IQC-based robustness analysis. They also do not imply the choice of arbitrary parameters such as the type, the length or the poles and zeros location of some basis functions, which makes solutions less sensitive to engineering judgment. Note that the adoption of frequency-by-frequency independent multipliers is made possible by the avoidance of the KYP lemma, which requires state-space realizations for all scalings. The main drawback of using such DG multipliers, however, is that the number of optimization parameters greatly increases with the size of the finite frequency grid over which the coupled FDIs shall be solved. In order to reduce the computational burden, therefore, a suboptimal optimization algorithm inspired by Ferreres and Roos (2007) is also implemented. The idea consists in iteratively optimizing with respect to constant and frequency-dependent variables. In this way, multiple small-size SDPs are solved one after the other until the algorithm converges to a suboptimal solution.

The paper is organized as follows. Section 2 discusses some preliminaries to robustness analysis in the presence of LTI and LTV uncertainties with bounded rates of variation. Sufficient conditions to assess stability and $\ell_2$-gain performance are then demonstrated in Section 3. The Hamiltonian-based technique used to validate the results on a continuous frequency range is detailed in Section 4. Eventually, the theoretical results are implemented in Section 5 and applied to test cases of increasing complexity in Section 6.

### Notation

The notation used in the paper is standard. Bold font is adopted for vectors and matrices, whereas roman characters indicate scalars. The supremum of a set or a function is denoted by $\sup$. Given a matrix $\mathbf{A}$, $\mathbf{A}^T$ is its transpose and $\mathbf{A}^*$ its conjugate transpose. $I$ stands for the identity matrix. The inequality relations between Hermitian matrices are expressed by $<$, $\leq_i$, $>$ and $\geq$. Given $n$ matrices $\{\mathbf{A}_i\}_{i=1,...,n}$, diag $(\mathbf{A}_1, ..., \mathbf{A}_n)$ denotes the standard block-diagonal augmentation. The transfer matrix of a system in terms of its state-space matrices is expressed as

$$ \begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{bmatrix} = C(s\mathbf{I} - \mathbf{A})^{-1} \mathbf{B} + \mathbf{D} $$

$\| \cdot \|$ indicates the $\ell_2$-gain of a linear bounded operator $T$: $\|T\| = \sup_{0 \neq x \in \ell_2^n} \frac{\|Tx\|_2}{\|x\|_2}$, where $x \in \ell_2^n$ means that $x$ is square-summable and $\|x\|_2$ is the $\ell_2$-norm of $x$. Open intervals are delimited by $]-[$, whereas closed intervals are bounded by $[ -]$. $\mathbb{R}$ represents the real space and $\mathbb{RH}$ denotes the space of all real-rational and proper transfer matrices whose poles have strictly negative real part. The Kronecker product is indicated with the standard symbol $\otimes$. Given two matrices of the form

$$ \Pi_i = \begin{bmatrix} \Pi_{i,11} & \Pi_{i,12} \\ \Pi_{i,12}^* & \Pi_{i,22} \end{bmatrix} $$

<figure><img src="image_4.png" /><figcaption>Fig. 1. Standard interconnections for robust stability (left) and robust performance (right).</figcaption></figure>

the duag operator represents the following relation:

$$ \text{duag}(\Pi_1, \Pi_2) = \begin{bmatrix} \Pi_{1,11} & \mathbf{0} & \mathbf{0} & \Pi_{1,12} & \mathbf{0} \\ \mathbf{0} & -\Pi_{1,12} & -\Pi_{2,11} & \mathbf{0} & -\Pi_{2,12} \\ \mathbf{0} & \Pi_{2,12}^* & \mathbf{0} & \Pi_{2,22} & \mathbf{0} \end{bmatrix} $$

## 2. Preliminaries

### 2.1. Definitions

Consider the standard interconnections represented in Fig. 1 for robust stability (left) and performance analysis (right). Let $\mathbf{M} = [\mathbf{M}_q, \mathbf{M}_u] = [\mathbf{M}_{qp}, \mathbf{M}_{qu}]$ be a continuous-time, stable and proper real-rational transfer matrix which represents the nominal closed-loop system. $\mathbf{\Delta} = \text{diag}(\Delta_I, \Delta_S)$ is a structured block-diagonal causal and bounded operator composed of $\Delta_I = \text{diag}(\Delta_{I_1}, ..., \Delta_{I_{n_I}})$, which gathers all the LTI uncertainties in the model, and $\Delta_S = \text{diag}(\Delta_{S_1}, ..., \Delta_{S_{n_S}})$, which contains all the time-varying parameters. The rate of variation of $\Delta_S$ at time $t$ is defined as in Köroğlu and Scherer (2007), Poolla and Tikku (1995):

$$ \bar{\sigma}(\dot{\Delta}_S(t)) = \lim_{dt \to 0} \bar{\sigma} \left( \frac{\Delta_S(t+dt) - \Delta_S(t)}{dt} \right) $$

where $\bar{\sigma}$ denotes the largest singular value of a matrix. Define also $\mathcal{V}_{\Delta_S}$ as the so-called variation operator such that $\mathcal{V}_{\Delta_S} = s\Delta_S - \Delta_S s$ where $s$ is the standard differentiation operator (Veenman, Scherer, & Köroğlu, 2016). Each $\Delta_I$ of $\Delta_I$ is a time-invariant diagonal matrix $\Delta_I = \delta_I \mathbf{I}$, where $\delta_I$ is a real or complex parametric uncertainty. Similarly, each $\Delta_I$ of $\Delta_S$ is a causal and bounded time-varying diagonal matrix $\Delta_I(t) = \delta_I(t) \mathbf{I}$ with bounded rate of variation. Let $\mathbf{\Delta}$ be the set of admissible uncertainties, i.e. the set of all uncertainties $\mathbf{\Delta}$ with the size and structure defined above. Let then

$$ \mathbf{B}_{\Delta} = \left\{ \mathbf{\Delta} \in \mathbf{\Delta} : \bar{\sigma}(\mathbf{\Delta}_I) \le \alpha, \bar{\sigma}(\mathbf{\Delta}_S(t)) \le \alpha, \bar{\sigma}(\dot{\mathbf{\Delta}}_S(t)) \le \beta \quad \forall t \in \mathbb{R} \right\} $$

where $\alpha$ and $\beta$ are scalar values which bound the size and the rate of variation of the uncertainties respectively. The input $p$ and output $q$ of $\mathbf{\Delta}$, assumed square-summable, are said to satisfy the IQC defined by $\Pi(j\omega)$ if

$$ \int_{-\infty}^{\infty} \left[ \begin{matrix} p(j\omega) \\ q(j\omega) \end{matrix} \right]^* \Pi(j\omega) \begin{bmatrix} p(j\omega) \\ q(j\omega) \end{bmatrix} d\omega \geq 0 \quad (1) $$

Thanks to IQCs it is possible to describe energy relations between two signals in a convenient form, capturing useful properties of uncertainties, time-varying parameters and non-linearities. In general, $\Pi(j\omega)$ can be any measurable Hermitian function, but in practice it is usually assumed to be rational.
## 2.2. Problem formulation

The present paper proposes a new approach to robust stability and performance analysis in the presence of LTI and LTV uncertainties with bounded rates of variation. The two problems can be formalized as follows:

**Problem 1 (Robust stability).** Determine if the interconnection of Fig. 1 (left) is stable $\forall \Delta \in B_\Delta$.

**Problem 2 (Robust $\ell_2$-gain performance).** Compute an upper bound $\gamma$ on the worst-case $\ell_2$-gain of the upper Linear Fractional Transformation (LFT) $\mathcal{F}_u(M, \Delta)$ between $u$ and $y$ in Fig. 1 (right):

$$ \sup_{\Delta \in B_\Delta} \| \mathcal{F}_u(M, \Delta) \| < \gamma $$

This performance indicator is particularly significant as it represents the well-known extension of the $H_\infty$ norm to LTV systems.

## 3. Proposed approach

Let $\mathcal{X} = \{X_I = X_I^* > 0 : \forall \Delta \in \Delta, X_I \Delta_I = \Delta_I X_I\}$ be the set of Hermitian, positive definite matrices $X_I$ which commute with $\Delta_I$ and correspond to the well-known $\mathcal{D}$ scalings of $\mu$-analysis. Let also $\mathcal{Y} = \{Y_I = -Y_I^* : \forall \Delta \in \Delta, Y_I \Delta_I = \Delta_I^* Y_I\}$ be the set of skew-Hermitian matrices which correspond to the $jG$ scalings of $\mu$-analysis. For each parametric LTV uncertainty $\delta_I(t)$ repeated $m_I$ times, define the following set of stable multipliers: $T_I(s) = C_{T_I}(sI - A_{T_I})^{-1} B_{T_I} + D_{T_I} \in \mathcal{R}H^{l_1 \times l_1}$, $V_I(s) = (sI - A_{T_I})^{-1} B_{T_I} \in \mathcal{R}H^{l_2 \times l_2}$ and $U_I(s) = C_{T_I}(sI - A_{T_I})^{-1} \in \mathcal{R}H^{l_1 \times l_2}$, where $l_1$ and $l_2$ can be chosen independently of each other, but to limit the number of optimization variables, the relation $v_I = l_1 - 1$ is imposed. The global multipliers $T(s)$, $V(s)$ and $U(s)$ have a block-diagonal structure $T(s) = \text{diag}(T_1(s) \otimes I_{m_1}, \dots, T_{n_S}(s) \otimes I_{m_{n_S}})$, $V(s) = \text{diag}(V_1(s) \otimes I_{m_1}, \dots, V_{n_S}(s) \otimes I_{m_{n_S}})$ and $U(s) = \text{diag}(U_1(s) \otimes I_{m_1}, \dots, U_{n_S}(s) \otimes I_{m_{n_S}})$. In addition, the following optimization matrices are defined for each parametric LTV uncertainty: $X_{S_I} = X_{S_I}^T > 0$, $X_{S_I} \in \mathbb{R}^{l_1 m_I \times l_1 m_I}$, $X_{V_{S_I}} = X_{V_{S_I}}^T > 0$, $X_{V_{S_I}} \in \mathbb{R}^{l_2 m_I \times l_2 m_I}$, $Y_{S_I} = -Y_{S_I}^T$, $Y_{S_I} \in \mathbb{R}^{l_1 m_I \times l_1 m_I}$, $Y_{V_{S_I}} = Y_{V_{S_I}}^T$, $Y_{V_{S_I}} \in \mathbb{R}^{l_2 m_I \times l_2 m_I}$. Analogously to the multipliers' case, the global optimization matrices can be defined with a block-diagonal structure as follows: $X_S = \text{diag}(X_{S_1}, \dots, X_{S_{n_S}})$, $X_{V_S} = \text{diag}(X_{V_{S_1}}, \dots, X_{V_{S_{n_S}}})$, $Y_S = \text{diag}(Y_{S_1}, \dots, Y_{S_{n_S}})$ and $Y_{V_S} = \text{diag}(Y_{V_{S_1}}, \dots, Y_{V_{S_{n_S}}})$.

Sections 3.1 and 3.2 present sufficient conditions to guarantee robust stability and $\ell_2$-gain performance of an uncertain linear system respectively. For LTI uncertainties new frequency-by-frequency independent scaling matrices are introduced, which in Section 3.3 are demonstrated to be less conservative than multipliers usually adopted in the literature.

### 3.1. Stability

The following theorem is a sufficient condition for stability of uncertain linear systems characterized by LTI and LTV uncertainties with bounded rates of variation. New frequency-by-frequency independent scaling matrices are introduced to describe the LTI uncertain parameters.

**Theorem 1.** Assume that the interconnection between $M_{qp}$ and $\tau\Delta$ (see Fig. 1 (left)) is well posed $\forall \Delta \in B_\Delta$ and $\forall \tau \in [0, 1]$. Let

$$ \Pi(j\omega) = \left[ - \frac{\Pi_{11}(j\omega)}{\Pi_{12}(j\omega)^*} - \frac{\Pi_{12}(j\omega)}{\Pi_{22}(j\omega)} \right] = \bar{T}(j\omega)^* P(\omega) \bar{T}(j\omega) $$

where

$$ \bar{T}(j\omega) = \begin{bmatrix} I & 0 & 0 & 0 & 0 \\ 0 & T(j\omega) & 0 & 0 & 0 \\ 0 & -\beta V(j\omega) & 0 & 0 & 0 \\ 0 & 0 & 0 & I & 0 \\ 0 & 0 & 0 & 0 & I \end{bmatrix} $$

and $P(\omega) = \text{duag}(P_I(\omega), P_S)$ with

$$ P_I(\omega) = \begin{bmatrix} - \frac{X_L(\omega)}{Y_I(\omega)^*} - \frac{Y_L(\omega)}{a^2 X_I(\omega)} \\ - \frac{1}{a^2} X_I(\omega) \end{bmatrix} $$

and

$$ P_S = \begin{bmatrix} X_S & 0 & 0 & 0 & 0 \\ 0 & X_{VS} & 0 & 0 & 0 \\ \frac{1}{Y_S^T} & -\frac{1}{0} & -\frac{0}{X_S} & -\frac{Y_{VS}}{0} & 0 \\ 0 & Y_{VS}^T & 0 & -X_{VS} & 0 \end{bmatrix} \quad (2) $$

If there exist $X_I(\omega) \in \mathcal{X}$, $Y_I(\omega) \in \mathcal{Y}$, $X_S$, $X_{VS}$, $Y_S$, $Y_{VS}$ defined as above such that

$$ \left[ \begin{matrix} - \frac{M_{qp}(j\omega)}{I} & 0 \\ 0 & I \end{matrix} \right]^* \Pi(j\omega) \left[ \begin{matrix} - \frac{M_{qp}(j\omega)}{I} & 0 \\ 0 & I \end{matrix} \right] < 0 \quad (3) $$

holds true $\forall \omega \in \mathbb{R}$, then the interconnection in Fig. 1 (left) is robustly stable $\forall \Delta \in B_\Delta$.

**Proof.** The objective of this proof is to demonstrate that $\Pi(j\omega)$ satisfies the IQC defined by Eq. (1). As $\Delta = \text{diag}(\Delta_I, \Delta_S)$, it is possible to prove the satisfaction of Eq. (1) separately for LTI and LTV uncertainties. For what concerns the LTV part, the proof is detailed in Kőroğlu and Scherer (2007) and is therefore omitted. Focusing on LTI uncertainties and calling $p_I(j\omega)$ and $q_I(j\omega)$ the input and output of $\Delta_I$ respectively, Eq. (1) is satisfied by $\Pi(j\omega)$ iff

$$ \int_{-\infty}^{\infty} \left[ \frac{p_I(j\omega)}{q_I(j\omega)} \right]^* P_I(\omega) \left[ \frac{p_I(j\omega)}{q_I(j\omega)} \right] d\omega \ge 0 \quad (4) $$

By developing the products, Eq. (4) becomes

$$ \int_{-\infty}^{\infty} p_I(j\omega)^* X_I(\omega) p_I(j\omega) d\omega + \int_{-\infty}^{\infty} q_I(j\omega)^* Y_I(\omega) p_I(j\omega) d\omega + \int_{-\infty}^{\infty} p_I(j\omega)^* Y_I(\omega) q_I(j\omega) d\omega - \frac{1}{a^2} \int_{-\infty}^{\infty} q_I(j\omega)^* X_I(\omega) q_I(j\omega) d\omega \ge 0 $$

Recalling from Fig. 1 (left) that $q_I(j\omega) = \Delta_I p_I(j\omega)$, one obtains

$$ \int_{-\infty}^{\infty} p_I(j\omega)^* X_I(\omega) p_I(j\omega) d\omega + \int_{-\infty}^{\infty} p_I(j\omega)^* \Delta_I^* Y_I(\omega)^* p_I(j\omega) d\omega + \int_{-\infty}^{\infty} p_I(j\omega)^* Y_I(\omega) \Delta_I p_I(j\omega) d\omega - \frac{1}{a^2} \int_{-\infty}^{\infty} p_I(j\omega)^* \Delta_I^* X_I(\omega) \Delta_I p_I(j\omega) d\omega \ge 0 $$

Furthermore, as the definition of $Y_I(\omega)$ yields $\Delta_I^* Y_I(\omega)^* = -Y_I(\omega) \Delta_I \forall \omega \in \mathbb{R}$, it holds that

$$ \int_{-\infty}^{\infty} p_I(j\omega)^* X_I(\omega) p_I(j\omega) d\omega - \frac{1}{a^2} \int_{-\infty}^{\infty} p_I(j\omega)^* \Delta_I^* X_I(\omega) \Delta_I p_I(j\omega) d\omega \ge 0 $$

or equivalently

$$ \int_{-\infty}^{\infty} X_I(\omega) d\omega - \frac{1}{a^2} \int_{-\infty}^{\infty} \Delta_I^* X_I(\omega) \Delta_I d\omega \ge 0 \quad (5) $$

Furthermore, as the definition of $X_I(\omega)$ implies that $X_I(\omega) \Delta_I = \Delta_I X_I(\omega) \forall \omega \in \mathbb{R}$, Eq. (5) becomes

$$ \int_{-\infty}^{\infty} X_I(\omega) d\omega - \frac{1}{a^2} \int_{-\infty}^{\infty} \Delta_I^* \Delta_I X_I(\omega) d\omega \ge 0 $$
If $\Delta \in B_{\Delta}$, moreover, $\bar{\sigma}(\Delta_I) \le \alpha$ and therefore

$$X_I(\omega) - \frac{1}{\alpha^2} \Delta_I^* \Delta_I X_I(\omega) \ge 0$$

holds true $\forall \omega \in \mathbb{R}$ and $\forall \Delta \in B_{\Delta}$, which verifies the IQC defined by Eq. (1). As the upper left corner $X_I(\omega)$ of $P_I(\omega)$ is positive definite and the lower right corner $-\frac{1}{\alpha^2}X_I(\omega)$ of $P_I(\omega)$ is negative definite, the IQC defined by Eq. (1) is verified by $\tau\Delta \forall \tau \in [0, 1]$ (Helmersson, 1999). Robust stability is eventually guaranteed by Megretski and Rantzer (1997). $\square$

## 3.2. $\ell_2$-gain performance

Similarly to what has been presented for stability, Theorem 2 is a sufficient condition for robust $\ell_2$-gain performance, where the multipliers associated to LTI uncertainties are frequency-by-frequency independent.

**Theorem 2.** Assume that the interconnection between $M$ and $\tau\Delta$ (see Fig. 1 (right)) is well posed $\forall \Delta \in B_{\Delta}$ and $\forall \tau \in [0, 1]$. Let

$$\Pi_{\gamma}(j\omega) = \begin{bmatrix} \Pi_{11}(j\omega) & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_{n_y} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I}_{n_{22}} \end{bmatrix} \begin{bmatrix} \Pi_{12}(j\omega) & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_{n_{22}} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I}_{n_{u}} \end{bmatrix}$$

where $n_u$ and $n_y$ are the size of $u$ and $y$ respectively, and $\Pi(j\omega) = \begin{bmatrix} \Pi_{11}(j\omega) & \mathbf{0} \\ \mathbf{0} & \Pi_{12}(j\omega) \end{bmatrix}^T$ is defined as in Theorem 1. If there exist $X_I(\omega) \in \mathbb{X}$, $Y_I(\omega) \in \mathbb{Y}$, $X_S$, $X_{VS}$, $Y_S$, $Y_{VS}$ defined as above, and $\gamma > 0$ such that

$$\begin{bmatrix} -\bar{M}(j\omega) \\ \bar{I} \end{bmatrix}^* \Pi_{\gamma}(j\omega) \begin{bmatrix} -\bar{M}(j\omega) \\ \bar{I} \end{bmatrix} < 0 \quad (6)$$

with $\bar{M}(j\omega) = [M_q(j\omega) \quad 0 \quad M_u(j\omega)]$, holds true $\forall \omega \in \mathbb{R}$, then the interconnection in Fig. 1 (right) is robustly stable $\forall \Delta \in B_{\Delta}$ and

$$\sup_{\Delta \in B_{\Delta}} ||\mathcal{F}_u(M, \Delta)|| < \gamma$$

**Proof.** If Eq. (6) is verified, the upper-left block of the LMI is definite negative (see Observation 7.1.2 of Horn & Johnson, 2012). By developing the upper-left block of Eq. (6), therefore, one obtains

$$\begin{bmatrix} M_{qp}(j\omega)^* \\ \mathbf{0} \end{bmatrix} \Pi(j\omega)_{11} \begin{bmatrix} M_{qp}(j\omega) & \mathbf{0} \end{bmatrix} + \begin{bmatrix} M_{qp}(j\omega)^* \\ \mathbf{0} \end{bmatrix} \Pi_{12}(j\omega) + \Pi_{12}(j\omega)^* \begin{bmatrix} M_{qp}(j\omega) & \mathbf{0} \end{bmatrix} + \Pi_{22}(j\omega) + \begin{bmatrix} M_{qy}(j\omega)^* \\ \mathbf{0} \end{bmatrix} \begin{bmatrix} M_{qy}(j\omega) & \mathbf{0} \end{bmatrix} < 0 \quad \forall \omega \in \mathbb{R}$$

or equivalently

$$\begin{bmatrix} -\bar{M}_{qp}(j\omega) & \mathbf{0} \\ \bar{I} & \mathbf{0} \\ \mathbf{0} & \bar{I} \end{bmatrix}^* \Pi(j\omega) \begin{bmatrix} -\bar{M}_{qp}(j\omega) & \mathbf{0} \\ \bar{I} & \mathbf{0} \\ \mathbf{0} & \bar{I} \end{bmatrix} + \begin{bmatrix} M_{qy}(j\omega)^* \\ \mathbf{0} \end{bmatrix} \begin{bmatrix} M_{qy}(j\omega) & \mathbf{0} \end{bmatrix} < 0 \quad \forall \omega \in \mathbb{R} \quad (7)$$

Knowing that $\begin{bmatrix} M_{qy}(j\omega)^* \\ \mathbf{0} \end{bmatrix} \begin{bmatrix} M_{qy}(j\omega) & \mathbf{0} \end{bmatrix} > 0 \quad \forall \omega \in \mathbb{R}$, Eq. (7) implies Eq. (3). The robust stability of the system is therefore guaranteed by Eq. (6).

Furthermore, by respectively pre- and post-multiplying both sides of Eq. (6) by $\begin{bmatrix} q(j\omega) \\ q_e(j\omega) \\ u(j\omega) \end{bmatrix}^*$ and $\begin{bmatrix} q(j\omega) \\ q_e(j\omega) \\ u(j\omega) \end{bmatrix}$, where $q_e(j\omega) = V_{\Delta_S}V_S(j\omega)p_S(j\omega)$

and $p_S(j\omega)$ is the input of $\Delta_S$, one obtains

$$\begin{bmatrix} q(j\omega) \\ q_e(j\omega) \\ u(j\omega) \end{bmatrix}^* \begin{bmatrix} -\bar{M}(j\omega) \\ \bar{I} \\ \bar{I} \end{bmatrix}^* \Pi_{\gamma}(j\omega) \begin{bmatrix} -\bar{M}(j\omega) \\ \bar{I} \\ \bar{I} \end{bmatrix} \begin{bmatrix} q(j\omega) \\ q_e(j\omega) \\ u(j\omega) \end{bmatrix} < 0 \quad \forall \omega \in \mathbb{R} \quad (8)$$

Knowing also that $p(j\omega) = M_{qp}(j\omega)q(j\omega) + M_{up}(j\omega)u(j\omega)$ and $y(j\omega) = M_{qy}(j\omega)q(j\omega) + M_{uy}(j\omega)u(j\omega)$, Eq. (8) becomes

$$\begin{bmatrix} p(j\omega) \\ y(j\omega) \\ q(j\omega) \\ q_e(j\omega) \\ u(j\omega) \end{bmatrix}^* \Pi_{\gamma}(j\omega) \begin{bmatrix} p(j\omega) \\ q(j\omega) \\ q_e(j\omega) \\ u(j\omega) \end{bmatrix} < 0 \quad \forall \omega \in \mathbb{R}$$

or equivalently

$$\begin{bmatrix} p(j\omega) \\ q(j\omega) \\ q_e(j\omega) \end{bmatrix}^* \Pi(j\omega) \begin{bmatrix} p(j\omega) \\ q(j\omega) \\ q_e(j\omega) \end{bmatrix} + y(j\omega)^* y(j\omega) - \gamma^2 u(j\omega)^* u(j\omega) < 0 \quad \forall \omega \in \mathbb{R} \quad (9)$$

Integrating Eq. (9) between $-\infty$ and $+\infty$, one obtains

$$\int_{-\infty}^{+\infty} \begin{bmatrix} p(j\omega) \\ q(j\omega) \\ q_e(j\omega) \end{bmatrix}^* \Pi(j\omega) \begin{bmatrix} p(j\omega) \\ q(j\omega) \\ q_e(j\omega) \end{bmatrix} d\omega + \int_{-\infty}^{+\infty} y(j\omega)^* y(j\omega) d\omega - \gamma^2 \int_{-\infty}^{+\infty} u(j\omega)^* u(j\omega) d\omega < 0 \quad (10)$$

Note that $p(j\omega)$, $q(j\omega)$ and $q_e(j\omega)$ satisfy the IQC defined $\Pi(j\omega)$ (see the proof of Theorem 1 for LTI uncertainties and Helmersson (1999), Köroğlu and Scherer (2007) for LTV uncertainties with bounded rates of variation). Eq. (10), therefore, implies that

$$\int_{-\infty}^{+\infty} y(j\omega)^* y(j\omega) d\omega - \gamma^2 \int_{-\infty}^{+\infty} u(j\omega)^* u(j\omega) d\omega < 0$$

which yields

$$\sup_{\Delta \in B_{\Delta}} \left( \frac{\int_{-\infty}^{+\infty} y(j\omega)^* y(j\omega) d\omega}{\int_{-\infty}^{+\infty} u(j\omega)^* u(j\omega) d\omega} \right)^{1/2} = \sup_{\Delta \in B_{\Delta}} ||\mathcal{F}_u(M, \Delta)|| < \gamma \quad \square$$

## 3.3. Comparison with existing IQC multipliers

Existing IQC solvers use the KYP lemma to transform the resulting FDIs into a single state-space LMI (Biannic et al., 2016; Kao et al., 2004; Veenman et al., 2021). The KYP lemma, however, requires a state-space representation of the FDIs and it cannot deal with the well-known frequency-by-frequency independent multipliers usually adopted in $\mu$-analysis to treat LTI uncertainties. To avoid this issue, Veenman et al. (2021) describes the LTI uncertainties with the same IQC multipliers used for LTV parameters with bounded rates of variation, namely $\Pi(j\omega) = \hat{T}(j\omega) P_S \hat{T}(j\omega)$, where

$$\hat{T}(j\omega) = \begin{bmatrix} T(j\omega) & \mathbf{0} & \mathbf{0} \\ \beta V(j\omega) & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & T(j\omega) & - \beta \bar{U}(j\omega) \end{bmatrix}$$

and $P_S$ is a constant matrix defined as in Eq. (2). Imposing $\beta = 0$ for LTI uncertainties and developing the stability condition expressed in Eq. (3) one obtains

$$M_{qp}(j\omega)^* T(j\omega)^* X_S T(j\omega) M_{qp}(j\omega) + T(j\omega)^* Y_S^T T(j\omega) M_{qp}(j\omega) + M_{qp}(j\omega)^* T(j\omega)^* Y_S T(j\omega) - \frac{1}{\alpha^2} T(j\omega)^* X_S T(j\omega) < 0$$

which coincides with the well-known stability condition of classical $\mu$-analysis (Young, Newlin, & Doyle, 1995):

$$T(j\omega)^* X_S T(j\omega) \rightarrow D \text{ scaling}$$

$$T(j\omega)^* Y_S^T T(j\omega) \rightarrow jG \text{ scaling}$$
$$ \frac{1}{\alpha^2} \rightarrow \beta^2 $$

The DG scalings of μ-analysis, however, are independent frequency-by-frequency, whereas $T(j\omega)^* X_S T(j\omega)$ and $T(j\omega)^* Y_S^T T(j\omega)$ imply a state-space representation. The multipliers proposed in the paper are, therefore, less conservative than the ones used by existing IQC solvers and do not require choosing arbitrary transfer matrices. This benefit is particularly compelling when the number of LTI uncertainties is large with respect to the number of LTV parameters, which is often the case in many applications, typically in the aerospace sector.

## 4. Hamiltonian-based method

The LMIs defined in Theorems 1 and 2 are frequency-dependent, which requires solving an infinite-dimensional problem. To make it finite, Eq. (3) (resp. Eq. (6)) is first solved on a finite frequency grid and a Hamiltonian-based technique is then used to validate the obtained solution on the entire frequency range. The method was initially formulated for continuous-time systems in the context of μ-analysis (Ferreres & Biannic, 1998) and it has then been applied by Parrilo (1999) to solve FDIs in case of strictly proper systems. A similar approach was also used by Demourant (2017) to perform robust stability analysis of time-varying systems. With respect to Demourant (2017), however, the present paper introduces frequency-by-frequency independent scalings for LTI uncertainties and adopts the multipliers proposed by Helmersson (1999) for LTV parameters with bounded rates of variation, which are less conservative than the ones of Jönsson and Rantzer (1996). Considering a solution of Eq. (3) (resp. Eq. (6)) composed of a matrix $P_S$ and N matrices $P_{I_i} = P_I(\omega_i)$ valid at N grid points $\omega_i$, the Hamiltonian-based algorithm is run N times, once for each $P_{I_i}$. At each iteration the objective is to identify the set $\Omega_i = \bigcup_j I_{I_j}(\omega) \subseteq [-\infty, +\infty[$ of all intervals $I_{I_j}$ for which $P_S$ and $P_{I_i}$ satisfy Eq. (3) (resp. Eq. (6)). The global validity region $\Omega_{valid}$ is the union of all $\Omega_i: \Omega_{valid} = \bigcup_j \Omega_i$. If $\Omega_{valid} \equiv [-\infty, +\infty[$, then the obtained solution is validated on the entire frequency range. Otherwise, some frequencies which do not belong to $\Omega_{valid}$ are added to the grid and the SDP defined by Theorem 1 (resp. Theorem 2) is solved again on the augmented grid, until a valid solution is obtained. The implemented Hamiltonian-based method is based on Theorem 3, inspired by Ferreres and Biannic (1998) and formulated in terms of robust stability.

**Theorem 3.** Assume that $P_{I_i}$ and $P_S$ verify Eq. (3) at $\omega_i$ and let $\psi_i(j\omega)$ be the left-hand side of Eq. (3), where $P(\omega)$ is replaced with $P_i = \text{duag}(P_{I_i}, P_S)$ in $\Pi(j\omega)$. Introduce a complex-valued matrix

$$ L = \begin{bmatrix} L_{11} & L_{12} \\ L_{21} & L_{22} \end{bmatrix} / \psi_i(j\omega) = F_u(L, j\delta_\omega I_{n\psi}) \quad (11) $$

where $n_{\psi}$ denotes the order of $\psi_i$ and $\delta_\omega = \omega - \omega_i$. Then $\psi_i(j\omega)$ becomes singular for all frequencies $\omega = \omega_i + \delta_\omega$ such that

$$ \det (I - j\delta_{\omega}H) = 0 \quad (12) $$

$$ \text{with } H = L_{11} - L_{12}L_{22}^{-1}L_{21}. $$

**Proof.** As $M_{ap}(s)$, $T(s)$, $V(s)$ and $U(s)$ are stable, it is readily checked that the LFT (11) is well-posed $\forall \omega \in \mathbb{R}$ as $L_{11}$ has no eigenvalue on the imaginary axis, thus implying that $I - j\delta_\omega L_{11}$ is non-singular. Observe then that $\psi_i(j\omega)$ is singular iff $\det(\psi_i(j\omega)) = \det(L_{22} + j\delta_\omega L_{21}(I - j\delta_\omega L_{11})^{-1}L_{12}) = 0$. By assumption, $L_{22} = \psi_i(j\omega_i) < 0$ is invertible and the above determinant equality can equivalently be rewritten as

$$ \det (I + j\delta_{\omega}L_{21}(I - j\delta_{\omega}L_{11})^{-1}L_{12}L_{22}^{-1}) = 0 $$

which also reads

$$ \det (I - j\delta_{\omega}L_{11} + j\delta_{\omega}L_{12}L_{22}^{-1}L_{21}) = 0 $$

or equivalently

$$ \det (I - j\delta_{\omega}(L_{11} - L_{12}L_{22}^{-1}L_{21})) = 0 $$

which is equivalent to Eq. (12). ☐

As a straightforward corollary of the above theorem, assuming that:

• $(\tilde{\omega}_k, \tilde{\omega}_{k+1})$ denotes a pair of two consecutive frequencies for which $\psi_i(j\omega)$ is singular,

• $\exists \hat{\omega} \in [\tilde{\omega}_k, \tilde{\omega}_{k+1}]$ [such that $\psi_i(j\hat{\omega}) < 0$]

then $\psi_i(j\omega) < 0 \quad \forall \omega \in [\tilde{\omega}_k, \tilde{\omega}_{k+1}]$. A similar condition can also be derived for the $\ell_2$-gain performance analysis. In this case, however, $P_{I_i}$ and $P_S$ shall satisfy Eq. (6) at $\omega_i$ for a given $\gamma > 0$ and $\psi_i(j\omega)$ is the left-hand side of Eq. (6), where $P(\omega)$ is replaced with $P_i$ in $\Pi_\gamma(j\omega)$. Note that in Theorem 3, $\delta_\omega$ is not arbitrarily chosen a priori, but is determined from the resolution of Eq. (12). As a consequence, this parameter does not affect either the numerical stability or the computational efficiency of the proposed method. Note also that the Hamiltonian-based technique guarantees that the solution obtained at a given frequency is valid on a continuous frequency range. There is therefore no need for a confidence level, as no approximation of the framework is carried out. Eventually, it is worth pointing out that LMIs and Hamiltonian-based method would still be applicable if the uncertainty block of the LFT contained a non-linear parameter. The only difference would be the introduction of an additional IQC to describe the non-linearity, which would increase the size of the SDP to be solved.

## 5. Implemented algorithms

The present section details two algorithms which have been implemented to assess robust stability and performance of linear systems characterized by LTI and LTV uncertainties with bounded rates of variation. Both of them are inspired by the method followed in Ferreres and Roos (2007) to perform robust feedforward design in the presence of LTI and arbitrarily-fast LTV uncertainties.

### 5.1. Optimal algorithm

The first algorithm consists in solving the LMI expressed in Eq. (3) (resp. Eq. (6) for $\ell_2$-gain performance analysis) on a finite frequency grid. The so-obtained solution is then validated on the entire frequency range through the Hamiltonian-based technique. Algorithm 1 presents the detailed scheme of the implemented approach, where $\psi_i(j\omega)$ indicates the left-hand side of Eq. (3) (resp. Eq. (6)) as in Theorem 3.

**Algorithm 1 Robust analysis: optimal algorithm.**

1: Select an initial (small size) frequency grid $\{\omega_i\}_{i=1,...,N}$

2: Solve the SDP of Theorem 1 (resp. Theorem 2 in case of $\ell_2$-gain performance analysis) on the grid to get one matrix $P_S$ and N matrices $P_{I_i}$ which satisfy Eq. (3) (resp. Eq. (6)) in $\omega_i$, $i = 1, ..., N$. If no feasible solution is found, the algorithm stops and the robust stability (resp. performance) is not guaranteed.

3: **while** The Hamiltonian-based method yields that Eq. (3) (resp. Eq. (6)) is not verified $\forall \omega \in \mathbb{R}$ **do**

4: Find all frequency intervals $\{\bar{I}_j\}_{j=1,...,N_p}$ for which $\nexists P_{I_i}$, $i = 1, ..., N$, such that Eq. (3) (resp. Eq. (6)) is satisfied

5: Update the current frequency grid by adding the middle frequency of each $\bar{I}_j$

6: Apply Step 2

7: **end while**

8: If the while loop ends, the system is robustly stable (resp. it has robust $\ell_2$-gain performance $\gamma$), otherwise stability (resp. performance) is not guaranteed
Note that when running the stability analysis, Algorithm 1 assesses the feasibility of Eq. (3) for given values of $\alpha$ and $\beta$ using the MATLAB built-in function `feasp`. The performance analysis, on the other hand, minimizes $\gamma$ for given values of $\alpha$ and $\beta$ with `mincx`. Note also that it is sometimes possible to speed up the convergence of the algorithm by imposing the LMIs at multiple frequencies inside each interval $\bar{I}_j$, $j = 1, \dots, N_p$ for which the obtained solution is still not valid. By doing so, the number of optimization variables does not change, but the LMIs shall hold on more than one frequency point inside a given interval (and not only in the middle frequency). Note eventually that Algorithm 1 does not require any arbitrary initial guess and it therefore converges to the optimal solution. Its main disadvantage, however, is that the LMIs to be solved cannot be decoupled frequency by frequency as $P_S$ is a constant matrix which shall satisfy Eq. (3) (resp. Eq. (6)) $\forall \omega \in \mathbb{R}$. The computational time can therefore become very high when the size of the frequency grid increases. This is due to the fact that the number of optimization matrices $P_{I_i}$ corresponds to the number of grid points. To overcome this issue, a new algorithm is proposed in Section 5.2, which allows solving multiple frequency-by-frequency decoupled SDPs, but does not guarantee the convergence to the optimal solution.

## 5.2. Suboptimal algorithm

The second implemented algorithm separately optimizes constant and frequency-by-frequency independent optimization parameters in an iterative way. By doing so, it is possible to reduce the number of optimization variables of each SDP, which allows decreasing the computational time with respect to Algorithm 1 in the presence of many uncertainties and large frequency grids. The algorithm, however, requires an initial guess for the constant (i.e. frequency independent) optimization variables of $P_S$ and it therefore converges to suboptimal solutions only. Note that a first guess for $P_S$ can be obtained by performing Step 2 of Algorithm 1 on a coarse frequency grid. The algorithm to perform suboptimal optimization is similar to Algorithm 1 except that Step 2 is replaced by the following sequences:

2a. Set a first guess of $P_S$.

2b. For the current values of $P_S$, solve the SDP defined by Theorem 1 (resp. Theorem 2 in case of $\ell_2$-gain performance analysis) $N$ times, once at each $\omega_i$ to determine $P_{I_i}$, $i = 1, \dots, N$.

2c. For the current values of $X_{I_i}$ and $Y_{I_i}$, $i = 1, \dots, N$, solve the SDP defined by Theorem 1 (resp. Theorem 2) on the finite frequency grid $\{\omega_i\}_{i=1,\dots,N}$ with respect to $P_S$.

2d. Repeat Steps 2b and 2c until a feasible solution is found. If a maximum number of iterations is reached, the algorithm stops and the robust stability (resp. performance) is not guaranteed.

Note that, analogously to Algorithm 1, the suboptimal algorithm uses the MATLAB built-in function `feasp` to assess the feasibility of Eq. (3) for given $\alpha$ and $\beta$ when running the stability analysis. If dealing with $\ell_2$-gain performance, however, the algorithm only assesses the feasibility of Eq. (6) for given values of $\alpha$, $\beta$ and $\gamma$ with `feasp`. This is due to the fact that alternatively minimizing $\gamma$ with respect to $P_{I_i}$, $i = 1, \dots, N$ and $P_S$ often yields unfeasible solutions. This could be avoided by modifying some parameters of the optimization (for example relaxing the required accuracy of the solutions), but this operation would be quite arbitrary and different models would require different settings. To make the algorithm numerically more stable, therefore, a threshold $\gamma$ is imposed before running the computations. Furthermore, the suboptimal algorithm can also be accelerated by imposing the LMIs at multiple frequencies inside each interval $\bar{I}_j$, $j = 1, \dots, N_p$ (and not only in the middle frequency). Note eventually that this algorithm may become ineffective when the number of LTI parameters is large with respect to the number of LTI uncertainties. In this scenario, in fact, it would not be possible to significantly reduce the size of the

SDPs to be solved by separately optimizing frequency dependent and constant variables. Classical IQC-based solvers may therefore outcome the performance of the proposed approach. In aerospace applications, however, time-varying parameters are usually less common than LTI uncertainties.

# 6. Applications

All results have been obtained using LMILAB on a Windows 10 laptop with a 13th Gen Intel(R) Core(TM) i5-1345U as CPU and 16 GB of RAM.

## 6.1. Validation example

Consider the two-inputs one-output interconnection proposed by Veenman et al. (2016) and represented in Fig. 2. $H = H(s) = \begin{bmatrix} A_w & B_w \\ C_w & D_w \end{bmatrix}$ is a stable linear uncertain system with

$$A_w = \begin{bmatrix} 0 & -1 + \delta_1 \\ 1 & -0.5 + \delta_2(t) \end{bmatrix}, \quad B_w = \begin{bmatrix} -2 & 0 \\ 1 & 0 \end{bmatrix},$$

$$C_w = \begin{bmatrix} 1 & 0 \\ -1 & 1 \end{bmatrix}, \quad D_w = \begin{bmatrix} 0 & 0 \\ 0 & 0.01 \end{bmatrix}$$

where $\delta_1$ is an LTI uncertainty and $\delta_2(t)$ is an LTV parameter such that $|\delta_2(t)| < 0.1 \forall t \in \mathbb{R}$. $E = E(s)$, furthermore, is a stable linear filter with state-space representation:

$$E(s) = \begin{bmatrix} -191 & 190 \\ 109.5 & -109 \\ 25.6 & 0 \end{bmatrix} \begin{bmatrix} -7.461 \\ 4.24 \\ 0 \end{bmatrix}$$

The objective of the filter is to process the measurement signal $y_2$ in order to estimate the signal $y_1$ so that the $\ell_2$-gain from the disturbance $u$ to the estimation error $y$ is minimized (Veenman et al., 2016). Suppose the following transfer matrix is selected as multiplier for the LTV uncertainty:

$$T(j\omega) = \begin{bmatrix} 1 & \frac{1}{j\omega - \rho} & \frac{1}{(j\omega - \rho)^2} & \dots & \frac{1}{(j\omega - \rho)^n} \end{bmatrix}^T \otimes I_m \quad (13)$$

where $\rho$ and $\nu$ are the pole location and the degree of the multiplier respectively and $m$ is the number of repetitions of the uncertainty, in this case equal to 1. $U(s)$ and $V(s)$ can be constructed by considering the minimal state-space representation of Eq. (13) proposed by Veenman et al. (2016):

$$T(s) = \begin{bmatrix} \rho & 0 & \dots & \dots & 0 & 1 \\ 1 & \ddots & \ddots & \ddots & \vdots & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \dots & 0 & 1 & \rho & 0 \\ \multicolumn{6}{c}{\begin{bmatrix} 0 \\ I_\nu \end{bmatrix}} \end{bmatrix} \otimes I_m$$

Choosing $\rho = -1$ rad/s and $\nu = 1$, Algorithm 1 is run starting from a coarse frequency grid composed of 0 rad/s, 1 rad/s and 2 rad/s. For $\alpha = 0.49$ the algorithm takes $\sim 0.35$ s to guarantee stability (note that for $\alpha = 0.5$ stability is no more guaranteed). This result is checked with IQClab, which requires defining a multiplier also for the LTI uncertainty. By choosing a scaling of the form of Eq. (13), with $\rho = -1$ rad/s and $\nu = 1$, stability is guaranteed only for $\alpha < 0.47$. However, by imposing $\rho = -0.1$ rad/s for the multiplier of the LTI uncertainty, IQClab takes $\sim 0.35$ s to guarantee stability for $\alpha = 0.49$ (when $\alpha = 0.5$ stability is no more guaranteed). The same sensitivity to the choice of $\rho$ is found when running the robust $\ell_2$-gain performance analysis. In particular, Fig. 3 shows the values of $\gamma$ which are obtained for different values of $\alpha$. The results of Algorithm 1 are compared with the outcome of IQClab, where the multiplier of the LTV parameter is kept as before, whereas different combinations of $\rho$ and $\nu$ are tried for the multiplier of the LTI uncertainty.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Validation example block diagram (Veenman et al., 2016).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Robust ℓ₂-gain performance analysis with respect to α.</figcaption></figure>

From the figure, it is possible to assess that the results of IQClab are sensitive to the choice of ρ and v, which is generally quite arbitrary (see Section 3.3). Algorithm 1, on the contrary, does not require choosing a multiplier for LTI uncertainties and its results overlap the best ones obtained with IQClab, namely for ρ = −0.1 rad/s and v = 2. Note that similar considerations hold when fixing a different multiplier for the LTV uncertainty.

## 6.2. Satellite with sloshing

The suboptimal algorithm is then tested on a more demanding satellite benchmark. Since coupling effects are often negligible, the satellite is represented as a single-axis model. Let θ be the attitude angle of the satellite and J its moment of inertia, which carries an uncertainty of 20% around its nominal value of 1000 kg m². The SISO dynamics of the model can be expressed as

$$J\ddot{\theta} = T_W + T_S + T_F + T_D$$

where T<sub>W</sub> is the control input torque, T<sub>S</sub> is induced by the propellant slosh effects, T<sub>F</sub> is produced by the flexible modes of the solar arrays and T<sub>D</sub> includes all remaining disturbance torques. As it is common in the literature, T<sub>S</sub> and T<sub>F</sub> are represented as poorly damped second-order linear models:

$$T_S + T_F = \sum_i \frac{L_i s^2}{s^2 + 2\xi_i \omega_i s + \omega_i^2} \ddot{\theta}$$

where L<sub>i</sub>, ω<sub>i</sub> and ξ<sub>i</sub> indicate respectively the magnitude, the frequency and the damping of the i<sup>th</sup> mode. Their nominal values are L<sub>1</sub> = 30 kg m², L<sub>2</sub> = 300 kg m², L<sub>3</sub> = L<sub>4</sub> = 100 kg m², L<sub>5</sub> = 50 kg m², ω<sub>1</sub> = 0.1 rad/s, ω<sub>2</sub> = 0.6 rad/s, ω<sub>3</sub> = 1 rad/s, ω<sub>4</sub> = 5 rad/s, ω<sub>5</sub> = 10 rad/s and ξ<sub>i</sub> = 10<sup>-3</sup>, i = 1, ..., 5. The first mode is due to sloshing and ω<sub>1</sub> is therefore considered as an LTV parameter. It is allowed to vary of 30% around its nominal value with a maximum rate of variation of 3 × 10<sup>-4</sup> rad/s². All L<sub>i</sub>, ω<sub>i</sub> and ξ<sub>i</sub> of the flexible modes, namely for i = 2, ..., 5, furthermore, are supposed to carry an LTI uncertainty of 20% around their nominal values. The control input torque T<sub>W</sub> is generated from the commanded torque T<sub>C</sub> by a reaction wheel, approximated as a

unitary first-order transfer function with time constant τ = 0.5 s. A 5-th order controller is designed to ensure closed-loop stability and reject as much as possible the disturbance torques T<sub>D</sub>. The transfer function between T<sub>D</sub> and θ is used to evaluate the performance of system.

A first analysis is performed on a satellite model including only the sloshing mode and the first flexible mode. The chosen multiplier for the LTV parameter is of the form of Eq. (13), with ρ = −0.1 rad/s and v = 1. Starting from a frequency grid composed by a single point at 0.1 rad/s, Algorithm 1 guarantees stability for α = 1.4 in ~ 1.9 s. Analogously, starting from an initial guess computed at 0.1 rad/s with g<sub>ep</sub>, the suboptimal algorithm guarantees stability for α = 1.4 in ~ 0.3 s. These results are then compared with the outcome of IQClab. To do so, the same multiplier is maintained for the LTV parameter, whereas for each LTI uncertainty new scalings are chosen of the form of Eq. (13). The analysis is repeated multiple times for different values of ρ (between −1 rad/s and −0.1 rad/s with steps of 0.1 rad/s) and v (between 0 and 3). Imposing α = 0.9, however, no feasible solution is identified by IQClab. This result does not mean that the system is unstable, because there could exist a combination of multipliers for the LTI uncertainties so that IQClab yields a feasible solution. In practice, however, it is not possible to verify all possible combinations of the multipliers' parameters (note that IQClab takes ~ 1.9 s and ~ 6 s to analyze the stability of the system when v is equal to 1 and 2 respectively). The same procedure is then conducted to analyze the robust ℓ₂-gain performance of the system. Starting from an initial frequency grid equal to [0.1, 0.6] rad/s, Algorithm 1 takes ~ 7.6 s to obtain γ = 0.012 for α = 1. The solution is confirmed by the suboptimal algorithm, which, starting from an initial guess computed with mincx at [0.1, 0.6] rad/s, takes ~ 5.0 s to assess the feasibility of γ = 0.012 when α = 1. On the other hand, IQClab returns feasible solutions where stability analysis failed. Imposing α = 1, the program yields γ ~ 0.016 in ~ 3 s when v = 1 and γ ~ 0.010 in ~ 10 s when v = 2.

The results so far highlighted become even more evident when the number of LTI uncertainties is largely increased. Suppose now to include in the satellite's model all five sloshing and flexible modes. The so-obtained dynamic system has 1 LTV uncertainty (namely ω<sub>1</sub>) and 13 LTI parameters. Considering that the uncertainties related to the natural frequencies are repeated twice, the size of Δ is equal to 19. All the following analyses are performed by choosing the multiplier of the LTV parameter as in Eq. (13) with v = 1 and ρ = −0.1 rad/s. Starting from an initial guess obtained with g<sub>ep</sub> at 0.1 rad/s, the suboptimal algorithm guarantees stability for α = 0.1 in ~ 24 s. The same analysis is then performed with IQClab, choosing multipliers as in Eq. (13) with v = 1 rad/s for all the uncertainties. The resulting SDP to be solved, however, is computationally burdensome, as it involves 1773 optimization parameters. By imposing ρ = −0.1 rad/s for all the multipliers, IQClab takes ~ 332 s to return infeasibility when α = 0.1. The same result is obtained if, for each multiplier of the LTI uncertainties, ρ is set to a value taken from a grid of 10 equispaced frequency points between −0.1 rad/s and −1 rad/s. A similar procedure is then conducted for the ℓ₂-gain performance analysis. Starting with an initial guess computed with mincx at 0.1 rad/s, the suboptimal algorithm takes ~ 3 min to assess that ||F(M, Δ)|| < γ = 0.5 for α = 0.1. The optimization of IQClab performed with mincx, instead, takes ~ 1.23 hr to yield infeasibility (ρ is set to −0.1 rad/s for all the multipliers). Note that similar considerations hold when fixing a different value of ρ for the scaling of the LTV parameter. The method proposed in the paper, therefore, results particularly efficient when the number of LTI uncertainties is considerable with respect to the number of LTV parameters (which is often the case in practical applications).

## 7. Conclusions and future work

A new method is proposed in this paper to assess robust stability and ℓ₂-gain performance in the presence of LTI and LTV uncertainties with
bounded rates of variation. Sufficient conditions are verified in the frequency domain on a discrete grid and a Hamiltonian-based technique is then used to validate the obtained solutions on a continuous frequency range. This approach is particularly interesting for two main reasons. Firstly, the LTI uncertainties are described with the well-known μ-analysis frequency-by-frequency independent scalings, which are less conservative than state-space multipliers and do not require choosing arbitrary transfer matrices. Secondly, suboptimal solutions can be obtained by separately optimizing constant and frequency-dependent decision variables in an iterative way, which reduces the size of the SDPs to be solved. Future work will focus on possible extensions of the proposed approach to discrete-time systems and saturating actuators. Additional advantages to compute hard bounds on the generalized $H_2$-norm of uncertain time-varying systems should also be investigated.

### CRediT authorship contribution statement

**Tommaso Casati:** Conceptualization, Formal analysis, Methodology, Writing – original draft, Data curation, Investigation, Project administration, Writing – review & editing. **Clément Roos:** Conceptualization, Investigation, Project administration, Supervision, Visualization, Formal analysis, Methodology, Resources, Validation, Writing – review & editing. **Jean-Marc Biannic:** Conceptualization, Investigation, Project administration, Supervision, Visualization, Formal analysis, Methodology, Resources, Validation, Writing – review & editing. **Hélène Evain:** Funding acquisition, Resources, Validation, Writing – review & editing, Project administration, Supervision, Visualization, Methodology.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### References

Biannic, J.-M., Burlion, L., Demourant, F., Ferreres, G., Hardier, G., Loquen, T., et al. (2016). The SMAC toolbox: a collection of libraries for systems modeling, analysis and control. Available: https://w3.onera.fr/smac/ssiqc.

Biannic, J.-M., Roos, C., Bennani, S., Boquet, F., Preda, V., & Girouart, B. (2021). Advanced probabilistic μ-analysis techniques for AOCS control laws. *European Journal of Control*, **62**, 120–129.

Demourant, F. (2017). Stability analysis by a new algorithmic approach based on integral quadratic constraints: Application to an aircraft benchmark. *AerospaceLab Journal*, (13).

Dennehy, C., Bennani, S., Shankar, U., Vandersteen, J., VanZwieten, T., Von der Porten, P., et al. (2022). Verification and validation (V&V) of guidance & control systems: results from the first inter-agency workshop on GNC V&V. In *Proceedings of the Annual American Astronautical Society Guidance, Navigation, and Control Conference*.

Ferreres, G., & Biannic, J.-M. (1998). A μ analysis technique without frequency gridding. In *Proceedings of the American Control Conference*.

Ferreres, G., & Roos, C. (2007). Robust feedforward design in the presence of LTI/LTV uncertainties. *International Journal of Robust and Nonlinear Control*, **17**(14), 1278–1293.

Helmersson, A. (1999). An IQC-based stability criterion for systems with slowly varying parameters. In *Proceedings of the IFAC World Congress*.

Horn, R. A., & Johnson, C. R. (2012). *Matrix Analysis* (2nd ed.). Cambridge University Press.

Jönsson, U., & Rantzer, A. (1996). Systems with uncertain parameters - time-variations with bounded derivatives. *International Journal of Robust and Nonlinear Control*, **6**(9–10), 969–982.

Kao, C.-Y., Megretski, A., Jönsson, U., & Rantzer, A. (2004). A MATLAB toolbox for robustness analysis. In *Proceedings of the IEEE International Conference on Robotics and Automation*.

Köroğlu, H., & Scherer, C. W. (2007). Robust performance analysis for structured linear time-varying perturbations with bounded rates-of-variation. *IEEE Transactions on Automatic Control*, **52**(2), 197–211.

Megretski, A., & Rantzer, A. (1997). System analysis via integral quadratic constraints. *IEEE Transactions on Automatic Control*, **42**(6).

Parrilo, P. A. (1999). On the numerical solution of LMIs derived from the KYP lemma. In *Proceedings of the Conference on Decision and Control*.

Poolla, K., & Tikku, A. (1995). Robust performance against time-varying structured perturbations. *IEEE Transactions on Automatic Control*, **40**(9), 1589–1602.

Roos, C. (2018). *Advanced control laws design and validation - a set of methods and tools to bridge the gap between theory and practice* (Habilitation à diriger des recherches), Université de Toulouse.

Veenman, J., Scherer, C. W., Ardura, C., Bennani, S., Preda, V., & Girouart, B. (2021). IQClab: A new IQC based toolbox for robustness analysis and control design. *IFAC-PapersOnLine*, **54**(8), 69–74.

Veenman, J., Scherer, C. W., & Köroğlu, H. (2016). Robust stability and performance analysis based on integral quadratic constraints. *European Journal of Control*, **31**, 1–32.

Young, P. M., Newlin, M. P., & Doyle, J. C. (1995). Computing bounds for the mixed μ problem. *International Journal of Robust and Nonlinear Control*, **5**(6), 573–590.