# LFT modelling and $\mu$-based robust performance analysis of hybrid multi-rate control systems

Jean-Marc Biannic, Clément Roos, Christelle Cumer

► To cite this version:

Jean-Marc Biannic, Clément Roos, Christelle Cumer. LFT modelling and $\mu$-based robust performance analysis of hybrid multi-rate control systems. 2024 IEEE 63rd Conference on Decision and Control (CDC), Dec 2024, Milan, Italy. pp.3662-3668, 10.1109/CDC56724.2024.10886203. hal-04631666

HAL Id: hal-04631666

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

Submitted on 2 Jul 2024

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
# LFT modelling and $\mu$-based robust performance analysis of hybrid multi-rate control systems

Jean-Marc Biannic*, Clément Roos and Christelle Cumer

**Abstract**— This paper focuses on robust stability and $H_\infty$ performance analyses of hybrid continuous/discrete time linear multi-rate control systems in the presence of parametric uncertainties. These affect the continuous-time plant in a rational way which is then modeled as a Linear Fractional Transformation (LFT). Based on a zero-order-hold (ZOH) LFT discretization process at the cost of bounded quantifiable approximations, and then using LFT-preserving down-sampling operations, a single-rate discrete-time closed-loop LFT model is derived. Interestingly, for any step inputs, and any admissible values of the uncertain parameters, the outputs of this model cover those of the initial hybrid multi-rate closed-loop system at every sampling time of the slowest control loop. Such an LFT model, which also captures the discretization errors, can then be used to evaluate both robust stability and guaranteed $H_\infty$ performance with a $\mu$-based approach. The proposed methodology is illustrated on a realistic and easily reproducible example inspired by the validation of multi-rate attitude control systems.

## I. INTRODUCTION

Although a few methods were already proposed over 30 years ago [1], [2], [3], [4] but also quite recently [5] to design multi-rate control systems taking into account the hybrid nature (simultaneously involving continuous and discrete time dynamics) of a closed-loop system, common engineering practice still consists in designing controllers in the continuous-time domain. This is easily explained by the high level of maturity and the numerous tools available in the field of robust continuous-time control, which has gained in popularity in recent years.

As a result, to shorten time-consuming simulation-based validation campaigns, there is a growing need for advanced and possibly cheap validation methodologies. Their main objective, which is also central in this paper, is to guarantee (at a reasonably low cost) closed-loop stability and performance after digital implementation of control laws often involving several loops operating at different rates. Unsurprisingly, after the early work of Kranc [6], modeling and analysis techniques for hybrid multi-rate control systems have then received a great deal of attention from the control community in recent decades [7], [8], [9], [10].

Many contributions in this field are based on time or frequency lifting techniques which, given a few assumptions on the sampling rates, allow the multi-rate system to be rewritten as an augmented single-rate model on which classical analysis techniques (such as gain or phase margins evaluations) become applicable [7], [8], [10]. More precisely, it can be observed that a multi-rate sampled-data system

can be first rewritten as a Linear Periodically Time Varying System (LPTVS) [11] and next transformed into an equivalent single-rate LTI model with expanded inputs and outputs. Such methods are not directly compatible, however, with the introduction of parametric uncertainties in a suitable way for $\mu$ or IQC based robustness analysis techniques, which rely on LFT modeling. An alternative approach is then developed in this paper to generate a single-rate full discrete time LFT model from the initial uncertain hybrid system. Given an uncertain continuous-time linear process in feedback with a discrete-time multi-rate controller, the proposed approach consists of the following steps:

1) zero-order-hold (ZOH) LFT-based discretization of the continuous-time model at the sampling rate of the fastest control loop,

2) single-rate closed-loop LFT generation by iterative closed-loop evaluations and LFT-preserving down-sampling operations from the fastest to the slowest control loop.

3) application of advanced $\mu$-analysis tools<fnref n="1" /> for robust stability and $H_\infty$ performance evaluation.

One of the key issues in the above procedure is the LFT-based discretization phase. As is clarified in the following sections, the exact ZOH discretization which must be considered in our context introduces non-rational terms (matrix exponentials) in the model. The LFT structure is then lost, and an important contribution of this paper is to show how a discrete-time LFT model can be obtained at the price of a quantified approximation error. It is also clarified that, under certain conditions, down-sampling operations preserve the LFT structure at the price of increased complexity.

The remainder of the paper is organized as follows. An overview of the modeling process is given in section II, which also clarifies the class of hybrid multi-rate systems under consideration. Next, the central results of the paper regarding LFT preserving discretization and down-sampling techniques are discussed in section III. Then, $\mu$-based robust performance is briefly detailed for the considered problem in section IV and an easily reproducible illustrative example is presented in section V. Final comments and future directions conclude the paper.

*Notation.* Given two operators $M$ and $N$ where $M$ is block-partitioned in a compatible way with the dimensions of

<footnote n="1">¹This evaluation can be performed either directly in the discrete-time domain for which $\mu$-based tools exist or can be adapted. Alternatively, a bilinear transformation (preserving both stability and the $H_\infty$ norm) can be applied to the discrete-time LFT so that more standard continuous-time $\mu$ tools become directly applicable.</footnote>
N, the lower-LFT (Linear Fractional Transformation), when it exists, is defined as $\mathcal{F}_l(M, N) = M_{11} + M_{12}N(I - NM_{22})^{-1}M_{21}$. Similarly, for a possibly different suitable partition, the upper-LFT is defined as $\mathcal{F}_u(M, N) = M_{22} + M_{21}N(I - NM_{11})^{-1}M_{12}$. Given any structure X whose elements are real/complex valued matrices or vectors, $\mathcal{B}_X$ denotes either the unit ball $\{X \in \mathbf{X}/\bar{\sigma}(X) \le 1\}$ or the hypercube $\{X \in \mathbf{X}/|X_i| \le 1\}$.

## II. OVERVIEW OF THE MODELING PROCESS

Consider the class of uncertain hybrid continuous/discrete time multi-rate systems illustrated by the closed-loop diagram in Fig. 1, where the state-space matrices of the $n^{th}$ order linear continuous-time uncertain process

$$G(s, \delta) := \begin{cases} \dot{x} = A(\delta)x + B(\delta)w \\ z = C(\delta)x + D(\delta)w \end{cases} \quad (1)$$

are assumed to depend rationally<fnref n="2" /> on the vector of normalized parametric uncertainties $\delta = [\delta_1, \dots, \delta_{n_\delta}]^T \in \mathcal{B}_\delta$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Uncertain hybrid continuous/discrete multi-rate closed-loop system</figcaption></figure>

Closed-loop stability and performance – evaluated through the $H_\infty$ norm of the transfer from the exogenous inputs $w_G$ to the exogenous outputs $z_G$ – are ensured by a multi-rate controller involving r loops. As shown in Fig. 1, each (possibly multivariable) loop operates at a specific frequency $1/T_i$ verifying Assumption 2.1, and is composed of the following three elements in cascade:

* $S_{T_i}(\cdot): n_{y_i} \times n_{y_i}$ sampler where $n_{y_i}$ represents the size of the vector $y_i$ of measurements updated every $T_i$ seconds. Note that $\{y_i\}_{i=1,...,r}$ form a partition of $y$. Then we have $\sum_{i=1}^r n_{y_i} = p$.

* $C_i(z):$ multivariable $n_{v_i} \times n_{y_i}$ discrete-time controller with the same period $T_i$ as above.

* $H_{T_i}(\cdot): n_{v_i} \times n_{v_i}$ zero-order-hold operator.

**Assumption 2.1:** The sampling periods $\{T_i\}_{i=1,...,r}$ are assumed to verify $T_{i+1} = q_i T_i$ with $q_i \in \mathbb{N}$.

Each loop finally delivers a continuous-time signal $v_i(t)$ from which the global control input $u(t)$ is obtained as the output of a linear static gain $L_\sigma \in \mathbb{R}^{m \times \sum n_{v_i}}$. Note that this matrix operator is designed to perform elementary operations (addition, subtraction) on the signals $v_i$ of compatible dimensions. It will then be essentially composed of 0, +1, or -1.

The main objective of the modeling process is to convert the uncertain hybrid closed-loop model of Fig. 1 into a single-rate discrete-time LFT as illustrated in Fig. 2, where

* $M(z)$ denotes an LTI discrete-time interconnection operating at the slowest rate $1/T_r$,

* $\Delta$ captures all the parametric uncertainties $\delta$ from the continuous-time process and also incorporates discretization errors such that, for any step input profile $w_G(t)$ with constant values $w_G(kT_r)$ on each interval $[kT_r, (k+1)T_r[$, the sampled trajectories $z_G(kT_r)$ of the hybrid multi-rate system are covered by those of the single-rate and full discrete-time version $z_M(kT_r)$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. LFT-based single-rate full discrete-time approximation</figcaption></figure>

To achieve this goal, the procedure described in Algorithm 1 is proposed. It involves two key technical steps (LFT preserving discretization and down-sampling), summarized by relations (2) and (4) further discussed in section III.

**Algorithm 1** Single-rate discretization of an LFT-based uncertain hybrid multi-rate system

**Require:** Uncertain hybrid closed-loop system of Fig. 1

**Ensure:** Single-rate discrete-time model of Fig. 2

**Initial step** Perform a ZOH discretization (see III-A) of the uncertain continuous-time model $G(s, \delta)$ with sampling period $T_1$ and generate a discrete-time LFT:

$$G(s, \delta) \rightarrow \mathcal{F}_u(H(z), \Delta_H) \quad (2)$$

Set $\Delta_1 = \Delta_H$ and close the fastest loop with sampling period $T_1$ to generate $M_1(z)$ by evaluating the lower-LFT:

$$M_1(z) = \mathcal{F}_l(H(z), C_1(z)) \quad (3)$$

**for** $i = 1, \dots, r-1$ **do**

**if** $q_i > 1$ **then**

Down-sample the upper-LFT from $T_i$ to $T_{i+1}=q_i T_i$:

$$\mathcal{F}_u(M_i(z), \Delta_i) \rightarrow \mathcal{F}_u(M_{i+1}(z), \Delta_{i+1}) \quad (4)$$

**else**

Set $M_{i+1}(z) = M_i(z)$ and $\Delta_{i+1} = \Delta_i$

**end if**

Update $M_{i+1}(z)$ by closing the loop with $C_{i+1}(z)$:

$$M_{i+1}(z) \leftarrow \mathcal{F}_l(M_{i+1}(z), C_{i+1}(z)) \quad (5)$$

**end for**

**Final step:** Set $M(z) = M_r(z)$ and $\Delta = \Delta_r$

<footnote n="2">²This assumption is required for LFT modeling, but is not so restrictive in practice since any continuous nonlinear function assumes rational approximations on a bounded set.</footnote>
### III. LFT PRESERVING ZOH DISCRETIZATION & DOWN-SAMPLING TECHNIQUES

#### A. ZOH discretization of a continuous-time LFT model

As detailed in [12], there exist many different techniques to discretize a continuous-time LFT model. Based on a bilinear transformation of the Laplace variable $s$, the well-known Tustin's method is very interesting as it fully preserves the LFT structure without any augmentation of the $\Delta$ block. Unfortunately, this approach is not suitable here, as the discrete-time output signals are not guaranteed to match the continuous-time signals at the sampling times ($z_k \neq z(kT)$). A possible alternative would be a "full ZOH" discretization of $M_G(s)$ where $G(s, \delta) = \mathcal{F}_u(M_G(s), \Delta_G(\delta))$. However, as observed in [12], this approach becomes inexact as soon as $\delta \neq 0$ and the only rigorous ZOH discretization of (1) reads

$$ \begin{bmatrix} x_{k+1} \\ z_k \end{bmatrix} = \mathcal{H}_0(\delta) \begin{bmatrix} x_k \\ w_k \end{bmatrix} \quad (6) $$

with

$$ \mathcal{H}_0(\delta) = \begin{bmatrix} e^{A(\delta)T} & \int_0^T e^{A(\delta)\tau} B(\delta) d\tau \\ C(\delta) & D(\delta) \end{bmatrix} \quad (7) $$

Interestingly, for any step input signal $w(t) = w_k$, with $t \in [kT, (k+1)T]$, the above expression ensures that the outputs of the continuous-time system (1) match exactly those of (6) at every sampling time: $\forall k \ge 0, z(kT) = z_k$. However, the presence of matrix exponentials unfortunately destroys the rational dependence on the uncertain parameters $\delta$. Consequently, an LFT model cannot be derived directly from (6), which must first be rewritten using polynomial or rational approximations of the exponentials.

1) **First-order rational approximation:** A commonly used rational approximation of the matrix exponential is based on the first-order Padé approximant:

$$ e^{A(\delta)T} \approx \left(I - \frac{T}{2} A(\delta)\right)^{-1} \left(I + \frac{T}{2} A(\delta)\right) \quad (8) $$

By introducing a quantification of the error in the above approximation and then substituting the exponentials in (7), the following lemma is obtained:

**Lemma 3.1:** $\forall T \ge 0$ there exists a bounded positive real $\mu_1(T)$ such that

$$ \forall \delta \in \mathcal{B}_{\delta}, \exists \Delta_{\epsilon} / \mathcal{H}_0(\delta) = \mathcal{F}_u(\mathcal{H}_1(\delta, \Delta_{\epsilon}), T.I_n) \quad (9) $$

with

$$ \mathcal{H}_1(\delta, \Delta_{\epsilon}) = \begin{bmatrix} \frac{1}{2} A(\delta) & (I + \Delta_{\epsilon}) A(\delta) & (I - \Delta_{\epsilon}) B(\delta) \\ I & I & 0 \\ 0 & C(\delta) & D(\delta) \end{bmatrix} \quad (10) $$

and

$$ \bar{\sigma}(\Delta_{\epsilon}) \le \mu_1(T) \quad (11) $$

**Proof:** Using the power series of the matrix exponential $e^X = \sum_{k=0}^{\infty} \frac{1}{k!} X^k$, it is readily checked after standard matrix manipulations that (9) holds with

$$ \Delta_{\epsilon} = -\frac{T^2 A(\delta)^2}{12} \underbrace{\sum_{k=0}^{\infty} \frac{T^k A(\delta)^k}{(1 + \frac{k}{2})(1 + \frac{k}{3})k!}}_{\underset{T \to 0}{\longrightarrow} I_n} \quad (12) $$

from which the bound $\mu_1(T)$ is obtained as

$$ \mu_1(T) = \frac{T^2}{12} \max_{\delta \in \mathcal{B}_{\delta}} \bar{\sigma} \left( \sum_{k=0}^{\infty} \frac{T^k A(\delta)^{k+2}}{(1 + \frac{k}{2})(1 + \frac{k}{3})k!} \right) \underset{T \to 0}{\longrightarrow} 0 \quad (13) $$

If $A(\delta)$ affinely depends on $\delta$ (often verified in practice), a cheap approximation of $\mu_1(T)$ for low values of $T$ reads

$$ \mu_1(T) \approx \frac{T^2}{12} \max_{\delta \in V(\mathcal{B}_{\delta})} \bar{\sigma} (A(\delta)^2) \quad (14) $$

where $V(\mathcal{B}_{\delta})$ denotes the set of vertices of $\mathcal{B}_{\delta}$.

**Corollary 3.1:** There exists a linear interconnection matrix $\mathcal{L}_H$ and a block-diagonal uncertain operator

$$ \Delta_{\delta} = \operatorname{diag}(\delta_1 I_{n_1}, \dots, \delta_{n_{\delta}} I_{n_{n_{\delta}}}) \quad (15) $$

such that

$$ \mathcal{H}_1(\delta, \Delta_{\epsilon}) = \mathcal{F}_u(\mathcal{L}_H, \operatorname{diag}(\Delta_{\delta}, \Delta_{\epsilon})) \quad (16) $$

**Proof:** As emphasized by equation (10), this straightforward consequence of Lemma 3.1 results from the fact that $\mathcal{H}_1(\delta, \Delta_{\epsilon})$ affinely depends on

(i) the approximation error $\Delta_{\epsilon}$

$$ \mathcal{H}_1(\delta, \Delta_{\epsilon}) = \mathcal{H}_1(\delta, \mathbf{0}) + \begin{bmatrix} \Delta_{\epsilon} \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix} [\mathbf{0} \ A(\delta) - B(\delta)] \quad (17) $$

(ii) the state-space matrices $A(\delta)$, $B(\delta)$, $C(\delta)$, $D(\delta)$ which, themselves, by assumption, rationally depend on $\delta$.

The construction of $\mathcal{L}_H$ associated with the block-diagonal uncertainty structure $\Delta_H = \operatorname{diag}(\Delta_{\delta}, \Delta_{\epsilon})$ is easily realized from (17) with the help of uncertainty modeling tools based on uss<fnref n="3" /> or gss objects [13].

**Remark 3.1:** Both Lemma 3.1 and its corollary are illustrated by Fig. 3, where the LFT-based discretization (2) introduced in Algorithm 1 clearly appears:

$$ z_k = \mathcal{F}_u \left( H(z), \underbrace{\operatorname{diag}(\Delta_{\delta}, \Delta_{\epsilon})}_{\Delta_H} \right) w_k \quad (18) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. A graphical illustration of Lemma 3.1 and its corollary</figcaption></figure>

<footnote n="3">³see: https://www.mathworks.com/help/robust/ref/uss.html</footnote>
2) **Second-order approximation:** The upper-bound on the approximation error (11) can be drastically reduced by considering higher-order rational approximations of the matrix exponential. Notably, the second-order approximation:

$$e^X \approx \left( I - \frac{X}{2} + \frac{X^2}{12} \right)^{-1} \left( I + \frac{X}{2} + \frac{X^2}{12} \right) \quad (19)$$

permits to update equation (9) of Lemma 3.1 by replacing $H_1(\delta, \Delta_\epsilon)$ with $H_2(\delta, \Delta_\epsilon)$ as follows

$$H_2(\delta, \Delta_\epsilon) = \begin{bmatrix} \left(I - \frac{T}{6} A(\delta)\right) \frac{A(\delta)}{2} & (I + \Delta_\epsilon)A(\delta) & (I - \Delta_\epsilon)B(\delta) \\ I & I & 0 \\ 0 & C(\delta) & D(\delta) \end{bmatrix} \quad (20)$$

where $\Delta_\epsilon$ becomes

$$\Delta_{\epsilon} = -\frac{T^4 A(\delta)^4}{720} \sum_{k=0}^{\infty} \frac{T^k A(\delta)^k}{(1 + \frac{k}{3})(1 + \frac{k}{4})(1 + \frac{k}{5})k!} \quad (21)$$

so that the approximated bound (14) for small values of $T$ and affine parametric dependency becomes

$$\mu_2(T) \approx \frac{T^4}{720} \max_{\delta \in V(B_\delta)} \bar{\sigma} (A(\delta)^4) \quad (22)$$

In this way, the second-order approximation saves two orders of magnitude on the norm of $\Delta_\epsilon$. The price to pay is a greater complexity of the LFT model (18) with an increase in the size of $\Delta_\delta$ in (15), since $H_2(\delta, \Delta_\epsilon)$ not only depends on $A(\delta)$ but also on $A(\delta)^2$.

3) **Extension to higher orders:** Finally, on specific cases involving few uncertainties, higher-order approximations can be considered to reach greater accuracy. Indeed, the upper-bound $\mu_n(T)$ on the $n^{th}$ order approximation error verifies $\mu_n(T) \propto \max_{\delta \in B_\delta} \bar{\sigma} ((T.A(\delta))^2^n)$.

**Remark 3.2:** Assuming the spectral radius of $T.A(\delta)$ is bounded by 1 throughout the uncertainty domain (always verified in practice), it is easily checked that $\bar{\sigma}(\Delta_\epsilon) \underset{n \to \infty}{\longrightarrow} 0$.

## B. Down-sampling: an LFT preserving operation

Let us now focus on the down-sampling operation summarized by equation (4) of Algorithm 1, and clarify that the LFT structure of the model is preserved. To do so, the following lemma is introduced.

**Lemma 3.2:** Consider $F_\tau(z) = C(zI - A)^{-1}B + D$, a discrete-time system with $n$ states, $m$ inputs and sampling period $\tau$. Then, for all $q \in \mathbb{N}$, the down-sampled system with period $q\tau$ is $F_{q\tau}(z) = C(zI - A_q)^{-1}B_q + D$ with

$$[A_q \quad B_q] = [I_n \quad 0] \begin{bmatrix} A & B \\ 0 & I_m \end{bmatrix}^q \quad (23)$$

*Proof:* This standard result is readily obtained by computing $x_{k+q} = A_q x_k + B_q w_k$ from $x_{k+1} = A x_k + B w_k$ and holding $w_k$ constant ($w_k = w_{k+1} = \dots = w_{k+q-1}$). ■ Now, following the notation of Algorithm 1, after closing the first $i$ control loops, the system $F_u(M_i(z), \Delta_i)$ to be down-sampled from period $T_i$ to $T_{i+1} = q_i T_i$ is an upper-LFT. Then, replacing the fixed matrices $[A \quad B]$ in Lemma 3.2 by $[A(\Delta_i) \quad B(\Delta_i)]$, equation (23) can be re-interpreted as the

product of $q = q_i$ LFTs. As a result, and further assuming without a significant loss of generality in practice that the $C$ and $D$ matrices are fixed, it is readily checked that $\Delta_{i+1}$ in the down-sampled model $F_u(M_{i+1}(z), \Delta_{i+1})$ verifies

$$\Delta_{i+1} = \operatorname{diag} \left( \underbrace{\Delta_i, \dots, \Delta_i}_{q_i \text{ times}} \right) \quad (24)$$

As a result

$$\Delta = \operatorname{diag} \left( \underbrace{\Delta_H, \dots, \Delta_H}_{N \text{ times}} \right) \quad \text{with} \quad N = \prod_{i=1}^{r-1} q_i \quad (25)$$

## IV. ON $\mu$-BASED ROBUST PERFORMANCE ANALYSIS

At this point, a standard discrete-time uncertain LFT model $z_M = F_u(M(z), \Delta)w_M$ is obtained, whose stability and $H_\infty$ performance can now be analyzed.

**Remark 4.1:** Since both stability and the $H_\infty$ norm are preserved under any bilinear transformation $s = k \frac{z-1}{z+1}$, which also preserves the LFT structure without altering the $\Delta$ block, the analysis can be carried out in the continuous-time domain on $F_u(P(s), \Delta)$ with $P(s) = M(-\frac{s+k}{s-k})$.

Let $\Delta$ be the set of all matrices with the same structure as $\Delta$ in (25). Recall that all uncertainties are normalized, so the considered uncertainty domain is simply $B_\Delta$. The following two quantities are evaluated using $\mu$-analysis:

* the robust stability margin:

$$k_r = \max \{k \ge 0 : P(s) - \Delta \text{ is stable } \forall \Delta \in kB_\Delta\}$$

* the worst-case $H_\infty$ performance level (if $k_r > 1$):

$$\gamma_{wc} = \max_{\Delta \in B_\Delta} \|F_u(P(s), \Delta)\|_{\infty}$$

The underlying theory is not presented here due to space limitations, but the interested reader can for example refer to [14], [15]. Only a few facts are briefly recalled to facilitate the understanding of section V-C. $\mu$-analysis basically consists of computing the peak value over the entire frequency range of the structured singular value $\mu_\Delta$. This computation being NP-hard in general, lower / upper bounds $\underline{\mu} / \bar{\mu}$ are usually determined instead of the exact value, from which bounds $\underline{k_r} / \bar{k_r}$ on $k_r$ and $\gamma_{wc} / \bar{\gamma}_{wc}$ on $\gamma_{wc}$ are derived. Much work has been done in the past decades to reduce the gap between these bounds, and (almost) exact values of $k_r$ and $\gamma_{wc}$ are now obtained in most cases with a reasonable computational time [16]. The main reason why the gap sometimes remains non-negligible and the computational time significant is the presence of uncertainties repeated many times in $\Delta$, which is precisely the case in this paper, see (25). A branch-and-bound algorithm can be used to overcome this issue. The uncertainty domain $B_\Delta$ is cut into smaller and smaller subsets until the relative gap between the highest lower bound and the highest upper bound on $k_r$ or $\gamma_{wc}$ computed on all subsets becomes less than a user-defined threshold. This algorithm is known to converge for uncertain systems with only real parametric uncertainties. However, it suffers from an exponential growth of computational complexity as a function of the number of
uncertainties. This can be alleviated using the μ-sensitivities,
which provide a very efficient way to detect the most critical
uncertainties, and therefore to decide in which directions to
cut the uncertainty domain to quickly reduce the gap. This
strategy is implemented in the routine *mubb* of the Matlab
SMART Library [17], which computes tight bounds on $k_r$
and $\gamma_{wc}$ with a very reasonable CPU time (see section V-C).

V. ILLUSTRATIVE EXAMPLE

The hybrid closed-loop model under consideration is illustrated by Fig. 4. The continuous plant $G(s, \delta)$ is a double integrator (pure inertia $J = 1$) with a positive feedback involving a poorly damped second order transfer function ($\alpha = 0.5, \xi = 0.001, \omega = 4$ rad/s). The actuator is represented by a first-order model with time-constant $\tau = 0.2$ s. Such a model is commonly used to represent the dynamics of a flexible satellite consisting of a central body and a solar panel. Normalized multiplicative uncertainties $\delta_J, \delta_\alpha, \delta_\omega$ and $\delta_\xi$ are then applied to the plant parameters to introduce variations of 10% on each of them. The gains $K_p = 1.65, K_i = 0.5$ and $K_d = 2.7$ of the PID controller are tuned by a standard continuous-time approach, so that the dominant mode of the nominal closed-loop system has a frequency of 0.5 rad/s. The nominal delay margin associated with this tuning is $\tau_d = 0.08$ s. This PID controller is then digitalized and implemented in a multi-rate setting as shown in Fig. 4. The derivative loop operates at sampling time $T_1$, while the proportional and integral (for which a forward-Euler approximation is used) loops share the same larger sampling time $T_2 = qT_1$.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Hybrid closed-loop model with a multi-rate PID controller</figcaption></figure>

A. Preliminary analysis: a counter-intuitive result?

A preliminary analysis is performed by applying a unit
step disturbance $w_p$ on the control input signal for 5 seconds.
The output signal $z = \dot{\theta}$ is first plotted (Fig. 5) in the
nominal parametric configuration ($\delta = 0$). Next, a critical
configuration (maximum values of $\alpha$ and $\omega$, and minimum
value of $J$) is displayed in Fig. 6. In both cases, the
continuous-time reference (black) is compared with two
hybrid implementations. The first one (visualized in red)
corresponds to a single rate implementation ($T_2 = T_1 = T$)
where $T \approx 2\tau_d$ is set at almost twice the delay margin (with
$\tau_{d0} \approx 0.08$ s in the nominal case and $\tau_d^* \approx 0.05$ s in the
critical case). The second one (magenta plots) implements a

multi-rate set-up with fixed sampling times $T_2 = 2T_1 = 0.2$ s
whatever the parametric configuration.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Nominal configuration: continuous-time vs hybrid single & multi-rate implementations</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Critical configuration: continuous-time vs hybrid single & multi-rate implementations</figcaption></figure>

In both configurations, the single-rate simulations (in red)
confirm the classical result, which predicts the stability of
the hybrid continuous/discrete system as long as the sam-
pling period remains less than twice the delay margin. But,
interestingly, as shown in Fig. 5, the proportional and integral
loops can accommodate a higher sampling period ($T_2 > 2\tau_d$)
by reducing that of the derivative loop ($T_1 < 2\tau_d$). With the
proposed tuning $T_2 = 2T_1 = 0.2$ s the oscillations on $\dot{\theta}$ are
strongly attenuated after 15 s.

More surprisingly (*a priori*), in the critical configuration,
the multi-rate implementation outperforms the single-rate
version, even though in both cases the derivative loop op-
erates at the same rate ($T_1 = 0.1$ s). Thus, by reducing the
rate (increasing $T_2$) of the proportional and integral "outer-
loops", the stability is improved, which actually makes sense,
but is not easily captured by analysis techniques in an
uncertain context. This simple example is then a very good
benchmark to evaluate the proposed methodology.

B. *LFT modeling of the hybrid multi-rate closed-loop model*

The first step is to run Algorithm 1. To do so, the open-loop plant $G(s, \delta) = F_u(M_G(s), \Delta_G)$ is first written in the LFT format. A fifth-order ($n_s = 5$) model is easily obtained with $\Delta_G = \text{diag}(\delta_J, \delta_\alpha, \delta_\omega I_2, \delta_\xi)$. Next, using the rational approximations discussed in subsection III-A, this model is discretized with $T_1 = 0.1$ s and the approximation errors are evaluated. Then, the derivative loop is closed, a down-sampling operation from $T_1$ to $T_2 = 2T_1$ is applied (see subsection III-B) and finally the proportional & integral
loops are closed. For comparison purposes, the algorithm is also applied after a Tustin discretization (LFT preserving but unfortunately not suited) or a full ZOH discretization (exact for $\delta = 0$, but inexact and even theoretically wrong for $\delta \neq 0$) of $M_G(s)$. The results are summarized in table I.

<table><caption>TABLE I
COMPLEXITY OF $F_u(M(z), \Delta)$ vs QUALITY OF THE APPROXIMATION</caption><thead><tr><th>Discretization method</th><th>structure of Δ<sub>δ</sub></th><th>size of Δ<sub>ε</sub></th><th>σ̄(Δ<sub>ε</sub>)</th></tr></thead><tbody><tr><td>Rational approx. (order 1)</td><td>I<sub>2</sub> ⊗ Δ<sub>G</sub></td><td>10 × 10</td><td>0.15</td></tr><tr><td>Rational approx. (order 2)</td><td>I<sub>4</sub> ⊗ Δ<sub>G</sub></td><td>20 × 20</td><td>0.0012</td></tr><tr><td>Full ZOH</td><td>I<sub>2</sub> ⊗ Δ<sub>G</sub></td><td>NA</td><td>NA</td></tr><tr><td>Tustin</td><td>I<sub>2</sub> ⊗ Δ<sub>G</sub></td><td>NA</td><td>NA</td></tr></tbody></table>

For each of the above four models, the simulations of subsection V-A are now replayed in the nominal case ($\delta = 0$) and for the same critical parametric configuration. As can be seen from the lower left subplot in Fig. 7, the method based on the full ZOH discretization exactly reproduces the output of the hybrid multi-rate system. This result was expected, since this discretization is exact without uncertainty. However, this is no longer the case in the critical configuration, where severe instability appears as revealed by the same subplot in Fig. 8.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Nominal configuration: hybrid multi-rate vs $F_u(M(z), 0)$</figcaption></figure>

Conversely, the Tustin-based approach (displayed on the lower-right subplots) introduces strong distortions whatever the configuration and fails to capture the oscillatory behavior (notably in the critical case).

Finally, for both configurations, the two closed-loop LFT models (see upper subplots in Fig. 7 and 8) generated from the rational approximations of subsection III.A produce good (even excellent with the second-order approximation) results.

**Remark 5.1:** Despite a higher complexity, the second-order approximation offers a major advantage in this application. As shown by the two upper-right subplots, and confirmed by the very small value of $\bar{\sigma}(\Delta_{\epsilon}) \le 1.2 \times 10^{-3}$, the accuracy of the associated LFT is indeed extremely high.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Critical configuration: hybrid multi-rate vs $F_u(M(z), \Delta^*)$</figcaption></figure>

### C. Robust stability & performance analysis

As discussed in section IV, a bilinear transformation is finally applied to the above discrete-time LFT models, where $\Delta_{\epsilon}$ is also introduced to obtain guaranteed results. In the general case, from equations (12) or (21), the simplest way would be to consider $\Delta_{\epsilon}$ as a real non-structured full block, for which specific scaling operators should be introduced in the $\mu$ upper-bound characterization. In this application however, due to the particular structure of the matrix $A(\delta)$, there are many zeros in $\Delta_{\epsilon}$. This property can be exploited to represent the approximation error by a limited number of scalar parametric uncertainties, which considerably reduces the conservatism.

1) **Robust stability:** Upper and lower bounds on the robust stability margin $k_r$ are first evaluated via a standard call to the *mubb* routine (see section IV). A first-order rational approximation (#1) without and with modeling errors is initially considered. Since the bound on $\Delta_{\epsilon}$ remains rather high in this case, robust stability cannot be ensured on the entire uncertainty domain ($\mu = 1/k_r = 2.31 > 1$). In the second-order case (#2), however, robust stability is easily ensured with $\bar{\mu} = 1/k_r = 0.86 < 1$. Moreover, this result is obtained with a very limited number of iterations of the branch-and-bound algorithm, which explains why only 5 s are required despite a more complex model. Note also that the conservatism introduced with the modeling error is very low ($\bar{\mu} = 0.86$ compared to 0.84 without the error), this time since the bound on $\Delta_{\epsilon}$ is very small.

<table><caption>TABLE II
MULTI-RATE ROBUST STABILITY ANALYSIS RESULTS</caption><thead><tr><th>Model type</th><th>μ̄</th><th>μ̅</th><th>CPU time<sup>4</sup></th></tr></thead><tbody><tr><td>Rational #1 (Δ<sub>ε</sub> = 0)</td><td>0.79</td><td>0.83</td><td>5 s</td></tr><tr><td>Rational #1 ((Δ<sub>ε</sub> ≠ 0)</td><td>2.31</td><td>2.43</td><td>27 s</td></tr><tr><td>Rational #2 (Δ<sub>ε</sub> = 0)</td><td>0.80</td><td>0.84</td><td>3 s</td></tr><tr><td><b>Rational #2</b> (Δ<sub>ε</sub> ≠ 0)</td><td><b>0.82</b></td><td><b>0.86</b></td><td><b>5 s</b></td></tr><tr><td>Full ZOH</td><td>2.41</td><td>2.53</td><td>0.2 s</td></tr><tr><td>Tustin</td><td>0.30</td><td>0.31</td><td>113 s</td></tr></tbody></table>
For comparison purposes, robust stability analysis is also applied to the models generated from the full ZOH and Tustin discretizations. In the first case, with $\bar{\mu} = 2.53$, stability is only proved for 4% variations of the parameters, while in the second case, with $\bar{\mu} = 0.31$, stability remains guaranteed much beyond the parametric domain, up to 33% variations. These results confirm the time-domain simulations in Fig. 8 and the invalidity of these two models.

2) **Robust performance:** Robust performance is now considered through the evolution of the $H_\infty$ norm of the transfer $T_{w_p \to z}(s)$ on the uncertainty domain. Upper and lower bounds on the worst-case $H_\infty$ performance level $\gamma_{wc}$ are computed for four different models generated from the second-order approximation without ($\Delta_\epsilon = 0$) and with ($\Delta_\epsilon \neq 0$) modeling error. The first two models (denoted MR) correspond to the multi-rate case ($T_2 = 2T_1 = 0.2$ s) while the last two (denoted SR-HF) correspond to a single-rate configuration ($T_2 = T_1 = 0.1$ s) at the highest frequency.

<table><caption>TABLE III
ROBUST $H_\infty$ PERFORMANCE ANALYSIS: MULTI-RATE vs SINGLE-RATE</caption><thead><tr><th>Model type</th><th>γ<sub>wc</sub></th><th>γ̃<sub>wc</sub></th><th>CPU time</th></tr></thead><tbody><tr><td>MR / Rational #2 (Δ<sub>ε</sub> = 0)</td><td>6.71</td><td>6.91</td><td>40 s</td></tr><tr><td><b>MR / Rational #2</b> (Δ<sub>ε</sub> ≠ 0)</td><td><b>7.12</b></td><td><b>8.55</b></td><td><b>100 s</b></td></tr><tr><td>SR-HF / Rational #2 (Δ<sub>ε</sub> = 0)</td><td>10.48</td><td>11.00</td><td>50 s</td></tr><tr><td>SR-HF / Rational #2 (Δ<sub>ε</sub> ≠ 0)</td><td>11.44</td><td>12.01</td><td>37 s</td></tr></tbody></table>

Very interestingly, a significant increase (beyond 40%) of the two bounds on the worst-case $H_\infty$ norm can be observed when the control loops all operate at the highest frequency. This result confirms the preliminary analysis and the fact that the flexible mode is more rapidly damped (whatever the parametric configuration) when the integral and proportional loops operate at a slower rate.

## VI. CONCLUSIONS AND FUTURE WORKS

A new approach has been proposed in this paper to describe an uncertain hybrid and multi-rate closed-loop system by a discrete-time and single-rate LFT model. Moreover, modeling errors have been quantified, allowing them to be integrated into a guaranteed $\mu$-based robustness analysis framework. Both the modeling process and the robustness analysis have been successfully evaluated on an easily reproducible and realistic benchmark.

Future work will be devoted to more specific adaptations of the robustness analysis tools to this new context of discrete-time and multi-rate systems. Notably, the case of the robust $H_2$ norm, whose value is not preserved by the bilinear transformation, should be carefully investigated. This metric is indeed of high practical interest to quantify pointing errors in space-oriented control applications.

## REFERENCES

[1] M. Berg, N. Amit, and J. Powell, “Multirate digital control system design,” *IEEE Transactions on Automatic Control*, vol. 33, no. 12, pp. 1139–1150, 1988.

[2] H. M. Al-Rahmani and G. F. Franklin, “Multirate control: A new approach,” *Automatica*, vol. 28, no. 1, pp. 35–44, 1992.

[3] P. G. Voulgaris and B. Bamieh, “Optimal $H_\infty$ and $H_2$ control of hybrid multirate systems,” *Systems and Control Letters*, vol. 20, no. 4, pp. 249–261, 1993.

[4] T. Chen and L. Qiu, “$H_\infty$ design of general multirate sampled-data control systems,” *Automatica*, vol. 30, no. 7, pp. 1139–1152, 1994.

[5] M. Cimino and P. R. Pagilla, “A design technique for multirate linear systems,” *IEEE Transactions on Control Systems Technology*, vol. 17, no. 6, pp. 1342–1349, 2009.

[6] G. Kranc, “Input-output analysis of multirate feedback systems,” *IRE Transactions on Automatic Control*, vol. 3, no. 1, pp. 21–28, 1957.

[7] I. S. Apostolakis and D. Jordan, “Multirate system performance evaluation using singular value analysis,” in *Proceedings of the American Control Conference*, 1990, pp. 1502–1509.

[8] H. Itop, H. Ohmori, and A. Sano, “Stability analysis of multirate sampled-data control systems,” *IMA Journal of Mathematical Control and Information*, vol. 11, no. 4, pp. 341–354, 12 1994.

[9] S. Longhi, “Structural properties of multirate sampled-data systems,” *IEEE Transactions on Automatic Control*, vol. 39, no. 3, pp. 692–696, 1994.

[10] M. van Haren, L. Blanken, and T. Oomen, “Frequency domain identification of multirate systems: A lifted local polynomial modeling approach,” in *Proceedings of the 61st IEEE Conference on Decision and Control*, 2022, pp. 2795–2800.

[11] S. Bittanti and P. Colaneri, *Periodic systems: Filtering and control*. London: Springer, 2009.

[12] R. Toth, M. Lovera, P. S. C. Heuberger, M. Corno, and P. M. J. Van den Hof, “On the discretization of Linear Fractional Representations of LPV systems,” *IEEE Transactions on Control Systems Technology*, vol. 20, no. 6, pp. 1473–1489, 2012.

[13] J.-M. Biannic and C. Roos, “Generalized State-Space objects: a new Matlab class to model uncertain and nonlinear systems as Linear Fractional Representations,” February 2016, available with the SMAC toolbox at http://w3.onera.fr/smac/gss.

[14] K. Zhou, J. Doyle, and K. Glover, *Robust and optimal control*. Prentice-Hall, 1996.

[15] G. Ferreres, *A practical approach to robustness analysis with aeronautical applications*. Kluwer Academic, 1999.

[16] C. Roos and J.-M. Biannic, “A detailed comparative analysis of all practical algorithms to compute lower bounds on the structured singular value,” *Control Engineering Practice*, vol. 44, pp. 219–230, 2015.

[17] C. Roos, “Systems Modeling, Analysis and Control (SMAC) toolbox: an insight into the robustness analysis library,” in *Proceedings of the IEEE International Symposium on Computer-Aided Control System Design*, 2013, pp. 176–181, available with the SMAC Toolbox at w3.onera.fr/smac/smart.