# Contrôle garanti de systèmes quantifiés par décomposition de l'espace des états d'un modèle différentiel réduit

Adrien Le Coënt, Florian de Vuyst, Christian Rey, Ludovic Chamoin, Laurent Fribourg

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

## To cite this version:

Adrien Le Coënt, Florian de Vuyst, Christian Rey, Ludovic Chamoin, Laurent Fribourg. Contrôle garanti de systèmes quantifiés par décomposition de l'espace des états d'un modèle différentiel réduit. 12e colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01400456

HAL Id: hal-01400456

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

Submitted on 22 Nov 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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 CC BY 4.0 - Attribution - International License
# Contrôle garanti de systèmes quantifiés par décomposition de l'espace des états d'un modèle différentiel réduit

A. Le Coënt<fnref n="1" />, F. De Vuyst<fnref n="1" />, C. Rey,<fnref n="2" /> L. Chamoin,<fnref n="2" /> L. Fribourg<fnref n="3" />

**Résumé** — Nous considérons le problème de synthèse du contrôle de systèmes linéaires à temps discret lorsque ce contrôle est quantifié (*quantized*), c'est-à-dire quand l'ensemble des valeurs que peut prendre le contrôle est *fini*. Nous montrons comment synthétiser un contrôle stabilisant en appliquant une méthode de décomposition de l'espace des états sur un modèle réduit du système. Nous calculons également une borne supérieure de l'erreur induite sur le modèle original. La procédure est illustrée sur plusieurs études de cas.

**Mots clés** — réduction d'ordre de modèle, contrôle garanti, synthèse de contrôle, contrôle stabilisant.

## 1 Introduction

Nous considérons un système d'équations différentielles, noté $Σ$, discret, linéaire et invariant dans le temps (DLTI) de la forme :

$$ \Sigma: \begin{cases} x(t + \tau) = A_d x(t) + B_d u(t), \\ y(t) = C_d x(t), \end{cases} \qquad (1) $$

où $x \in \mathbb{R}^n$ est la variable d'état, $y \in \mathbb{R}^m$ est la variable de sortie, $u(\cdot)$ est une fonction de contrôle à valeurs dans $\mathbb{R}^p$, $A_d \in \mathbb{R}^{n \times n}$, $B_d \in \mathbb{R}^{n \times p}$, $C_d \in \mathbb{R}^{m \times n}$, ces matrices sont indicées d'après elles définissent le système discret, et $\tau$ est une constante d'échantillonnage du temps. Ce système correspond à la discrétisation temporelle du système continu linéaire invariant dans le temps suivant :

$$ \begin{cases} \dot{x}(t) = Ax(t) + Bu(t), \\ y(t) = Cx(t), \end{cases} $$

avec

$$ A_d = e^{A\tau}, \quad B_d = \int_0^\tau e^{A(\tau-t)} Bdt, \quad C_d = C. $$

De façon essentielle ici, nous imposons que $u(t)$ prenne ses valeurs dans un ensemble *fini* de vecteurs de $\mathbb{R}^p$, noté $U$. Chaque élément de $U$ est appelé *mode*. Un tel contrôle est dit *quantifié* (voir [1]). Le problème du contrôle (discret) quantifié est de trouver une loi $u(t)$ qui à chaque instant d'échantillonnage $\tau, 2\tau, \dots$ trouve le bon élément de $U$ afin de réaliser un objectif donné, tel que la stabilisation autour d'un point cible $y_{obj}$. En général, la stabilisation de tels systèmes ne peut être parfaite [?], et on recherche seulement une stabilisation "pratique" : on recherche non pas un point d'équilibre $y_{obj}$, mais seulement un voisinage de $y_{obj}$ à l'intérieur duquel le système finit par être confiné.

**Entrées du problème :**

1. un sous-ensemble $R_x \subset \mathbb{R}^n$ de l'espace des états, appelé *ensemble d'intérêt*,

2. un sous-ensemble $R_y \subset \mathbb{R}^m$ de l'espace des sorties, appelé *ensemble cible*.

**Objectif** : trouver une loi de contrôle $u(\cdot)$ qui, pour tout point initial $x_0 \in R_x$, stabilise la sortie $y$ dans la zone $R_y$.

Dans [2], nous avons proposé une procédure (appelée "procédure de décomposition") dans le but d'atteindre un objectif analogue. En raison de son coût exponentiel en fonction de la dimension, la méthode de décomposition sera appliquée ici à un modèle réduit $\hat{\Sigma}$ de dimension $n_r < n$ afin de stabiliser

<footnote n="1">¹ CMLA, ENS Cachan, {lecoent,devuyst}@cmla.ens-cachan.fr</footnote>

<footnote n="2">² LMT, ENS Cachan, {rey, chamoin}@lmt.ens-cachan.fr</footnote>

<footnote n="3">³ LSV, ENS Cachan, fribourg@lsv.ens-cachan.fr</footnote>
sa sortie dans l'ensemble cible $R_y$. Nous montrerons que la loi de contrôle ainsi synthétisée stabilise également la sortie du système original $\Sigma$ dans $R_y$ avec une tolérance $\varepsilon$. Nous illustrons la méthode sur plusieurs exemples.

## 2 Procédure de synthèse de contrôle

Etant donné le système $\Sigma$, une zone d'intérêt $R_x$ et une zone cible $R_y$, la méthode consiste d'abord à réduire le système $\Sigma$ de dimension $n$ en $\hat{\Sigma}$ de dimension $n_r < n$. Nous appliquons ici la méthode classique de troncature équilibrée (balanced truncation) (voir [3, 4]). Le système $\hat{\Sigma}$ est lui-même un système DLTI défini par des matrices $\hat{A}_d, \hat{B}_d, \hat{C}_d$. Le système $\hat{\Sigma}$ est décrit sous la forme :

$$ \hat{\Sigma} : \begin{cases} \hat{x}(t + \tau) = \hat{A}_d \hat{x}(t) + \hat{B}_d u(t), \\ y_r(t) = \hat{C}_d \hat{x}(t), \end{cases} $$

avec

$$ \hat{A}_d = e^{\hat{A}\tau}, \quad \hat{B}_d = \int_0^\tau e^{\hat{A}(\tau-t)} \hat{B} dt, \quad \hat{C}_d = \hat{C}. $$

La méthode de balanced truncation est une méthode de réduction par projection, notée ici $\pi_r$. La projection de $R_x$ sera notée $\hat{R}_x$. On va synthétiser un contrôle sur $\hat{\Sigma}$ en lui appliquant une adaptation de la procédure de décomposition de [2]. Cette procédure prend en entrées :

1. $\hat{R}_x$ la projection de l'ensemble d'intérêt $R_x$,

2. $R_y$ l'ensemble cible,

3. $k$ un entier fixant le nombre maximal de modes successifs recherchés par la procédure.

En sortie, cette procédure trouve une décomposition $\hat{\Delta}$ de la forme $\{V_i, \sigma_i\}_{i \in I}$ telle que :

1. $I$ est un ensemble fini d'indices,

2. chaque $V_i$ ($i \in I$) est un produit d'intervalle de dimension $n_r$ tel que $\bigcup_{i \in I} V_i = \hat{R}_x$,

3. chaque $\sigma_i$ ($i \in I$) est une $k$-séquence, c'est-à-dire une suite d'au plus $k$ modes de $U$, telle que l'image de $V_i$ après application successive des modes de $\sigma_i$ soit contenue dans $\hat{R}_x$ et la sortie correspondante dans $R_y$.

La décomposition $\hat{\Delta}$ induit un contrôle $u_{\hat{\Delta}}$ sur $\hat{R}_x$. Considérons en effet un point initial quelconque $x_0$ de $R_x$. Soit $\hat{x}_0 = \pi_r x_0$ sa projection dans $\hat{R}_x$: $\hat{x}_0$ appartient à $V_{i_0}$ pour un certain $i_0 \in I$; donc après application de $\sigma_{i_0}$, on obtient un point $x_1$ de $\hat{R}_x$; donc, $x_1$ appartient lui-même à $V_{i_1}$ pour un certain $i_1 \in I$; et ainsi de suite, on construit, à partir de $\hat{x}_0$, une suite de points $\hat{x}_1, \hat{x}_2, \dots$ obtenus par application des $k$-séquences $\sigma_{i_0}, \sigma_{i_1}, \dots$ (étapes (1),(2) et (3) de la procédure hors ligne figure 1).

Ce contrôle $u_{\hat{\Delta}}$ défini sur le système réduit $\hat{\Sigma}$, peut être appliqué au système original $\Sigma$: on construit, à partir de $x_0$, une suite de points $x_1, x_2, \dots$ obtenus par application des $k$-séquences $\sigma_{i_0}, \sigma_{i_1}, \dots$ (étapes (4),(5) et (6) de la procédure hors ligne figure 1).

Pour tout $x_0 \in R_x$ et pour tout $t \ge 0$, une majoration $\varepsilon_y^\infty$ de l'erreur $\|y(x_0, u_{\hat{\Delta}}, t) - y_r(\pi_r x_0, u_{\hat{\Delta}}, t)\|$ peut être calculée en utilisant des inégalités algébriques classiques et le principe de superposition pour les systèmes linéaires (voir [5]). L'erreur $\|y(x_0, u_{\hat{\Delta}}, t) - y_r(\pi_r x_0, u_{\hat{\Delta}}, t)\|$ est elle-même bornée par $\varepsilon_y^j$ quand une $k$-séquence de longueur $j$ est appliquée, ce qui donne les résultats (2) et (3). $\varepsilon_y^\infty$ existe dans le cas où le module des valeurs propres des matrices $e^{\tau A}$ et $e^{\tau \hat{A}}$ est strictement inférieur à 1, ce que nous supposons.

$$ \varepsilon_y^j = \|u(\cdot)\|_{\infty}^{[0, j\tau]} \int_0^{j\tau} \| [C - \hat{C}] \begin{bmatrix} e^{tA} \\ e^{t\hat{A}} \end{bmatrix} \begin{bmatrix} B \\ \hat{B} \end{bmatrix} \| dt + \sup_{x_0 \in \hat{X}_0} \| [C - \hat{C}] \begin{bmatrix} e^{j\tau A} \\ e^{j\tau \hat{A}} \end{bmatrix} \begin{bmatrix} x_0 \\ \pi_R x_0 \end{bmatrix} \| \quad (2) $$

$$ \varepsilon_y^\infty = \sup_{j \in \mathbb{N}_>0} \varepsilon_y^j \quad (3) $$

**Notation** : Pour tout produit d'intervalle $R_y$ de la forme $[a_1, b_1] \times \dots \times [a_m, b_m]$, on note $R_y + \varepsilon$ le produit d'intervalle $[a_1 - \varepsilon, b_1 + \varepsilon] \times \dots \times [a_m - \varepsilon, b_m + \varepsilon]$.
**Proposition 1** *Considérons un système DLTI Σ, une zone d'intérêt R<sub>x</sub>, une zone cible R<sub>y</sub>. Soit Σ̂ la projection de Σ par troncature équilibrée. Soit Δ̂ une décomposition obtenue pour Σ̂ avec R̂<sub>x</sub>, R<sub>y</sub> et k ∈ N<sub>>0</sub> en entrées. Pour tout point x<sub>0</sub> ∈ R<sub>x</sub>, le contrôle u<sub>Δ</sub> induit appliqué à Σ en x<sub>0</sub> est tel que pour tout j ∈ N<sub>>0</sub> y(jτ) appartient à R<sub>y</sub> + ε<sub>y</sub><sup>∞</sup> où ε<sub>y</sub><sup>∞</sup> est défini en (3).*

3 Extension pour l'utilisation en ligne

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Schémas de procédures hors ligne (à gauche) et en ligne (à droite) pour des simulations de longueur 3.</figcaption></figure>

Jusqu'ici, la procédure de synthèse de contrôle consistait à calculer une suite de *k*-séquences sur le
modèle réduit $\hat{\Sigma}$, et à l'appliquer au modèle complet $\Sigma$. Cette loi de contrôle est donc calculée hors-ligne.
Cependant, la décomposition $\hat{\Delta}$ peut être utilisée en ligne de la façon suivante : soit $x_0$ un point initial
quelconque de $R_x$ et $\hat{x}_0 = \pi_r x_0$ sa projection appartenant à $\hat{R}_x$ (étape (1) de la procédure en ligne figure
1) ; $\hat{x}_0$ appartient à $V_{i_0}$ pour un certain $i_0 \in I$ ; on peut appliquer la *k*-séquence $\sigma_{i_0}$ au modèle complet $\Sigma$,
ce qui emmène $\Sigma$ à l'état $x_1$ (étape (2) de la procédure en ligne figure 1), et la sortie correspondante $y_1$ est
envoyée dans $R_y + \varepsilon_y^{|\sigma_{i_0}|}$. Afin de calculer une nouvelle *k*-séquence, il faut s'assurer que $\pi_r x_1$ appartient
bien à $\hat{R}_x$ (étape (3) de la procédure en ligne figure 1).

Notons $\hat{x}_1$ l'état dans lequel est emmené le système réduit $\hat{\Sigma}$ après application de la $k$-séquence $\sigma_{i_0}$.
La majoration de l'erreur $\|\pi_r x_1 - \hat{x}_1\|$ (erreur de trajectoire d'état ou state trajectory error), et un renfor-
cement de la procédure permettant de prendre en compte cette erreur vont nous permettre de garantir que
la projection de l'état du système complet $\pi_r x_i$ reste dans $\hat{R}_x$ à chaque étape. Notons $\varepsilon_x^j$ une majoration
de $\|\pi_r x_1 - \hat{x}_1\|$ quand une $k$-séquence de longueur $j$ est appliquée (donnée ci-dessous en 4). L'erreur
$\|\pi_r x_1 - \hat{x}_1\|$ est donc majorée par $\varepsilon_x^{|\sigma_{i_0}|}$.

$$
\epsilon_x^j = \|u(\cdot)\|_{\infty}^{[0, j\tau]} \int_0^{j\tau} \left\| \begin{bmatrix} \pi_R & -I_{n_r} \end{bmatrix} \begin{bmatrix} e^{tA} & \\ & e^{t\hat{A}} \end{bmatrix} \begin{bmatrix} B \\ \hat{B} \end{bmatrix} \right\| dt + \\
\sup_{x_0 \in R_x} \left\| \begin{bmatrix} \pi_R & -I_{n_r} \end{bmatrix} \begin{bmatrix} e^{j\tau A} & \\ & e^{j\tau \hat{A}} \end{bmatrix} \begin{bmatrix} x_0 \\ \pi_R x_0 \end{bmatrix} \right\| . \quad (4)
$$

Nous modifions maintenant la procédure de décomposition. Cette nouvelle procédure, que nous ap-
pelons ε-décomposition, prend en entrées :
1. $\hat{R}_x$ la projection de l'ensemble d'intérêt $R_x$,

2. $R_y$ l'ensemble cible,

3. $k$ un entier fixant le nombre maximal de modes successifs recherchés par la procédure.

4. $\{\epsilon_x^j\}_{1\le j\le k}$, un vecteur contenant les erreurs $\epsilon_x^1$ à $\epsilon_x^k$.

En sortie, cette procédure trouve une $\epsilon$-décomposition $\hat{\Delta}$ de la forme $\{V_i, \sigma_i\}_{i \in I}$ telle que :

1. $I$ est un ensemble fini d'indices,

2. chaque $V_i$ ($i \in I$) est un produit d'intervalles de dimension $n_r$ tel que $\bigcup_{i \in I} V_i = \hat{R}_x$,

3. chaque $\sigma_i$ ($i \in I$) est une $k$-séquence, c'est-à-dire une suite d'au plus $k$ modes de $U$, telle que l'image de $V_i$ après application successive des modes de $\sigma_i$ soit contenue dans $\hat{R}_x - \epsilon_x^{|\sigma_i|}$ et la sortie correspondante dans $R_y$.

En appliquant la procédure modifiée comme ci-dessus, nous pouvons effectuer un contrôle en ligne. La séquence $\sigma_{i_0}$ a été calculée de façon à envoyer $V_{i_0}$ dans $\hat{R}_x - \epsilon_x^{|\sigma_{i_0}|}$. Ainsi, $\hat{x}_1$ appartient à $\hat{R}_x - \epsilon_x^{|\sigma_{i_0}|}$. Donc $\pi_r x_1$ appartient $\hat{R}_x$, puisque $\epsilon_x^{|\sigma_{i_0}|}$ est une majoration de la distance maximale entre $\hat{x}_1$ et $\pi_r x_1$. Nous pouvons donc calculer la $k$-séquence $\sigma_{i_1}$ et itérer la procédure. Nous avons alors la proposition suivante :

**Proposition 2** *Considérons un système DLTI $\Sigma$, une zone d'intérêt $R_x$, une zone cible $R_y$. Soit $\hat{\Sigma}$ la projection de $\Sigma$ par troncature équilibrée. Soit $\hat{\Delta}$ une $\epsilon$-décomposition obtenue pour $\hat{\Sigma}$ avec $\hat{R}_x, R_y$, $k \in \mathbb{N}_{>0}$ et $\{\epsilon_x^j\}_{1\le j\le k}$ en entrées. Pour tout point $x_0 \in R_x$, le contrôle $u_\Delta$ induit appliqué à $\Sigma$ en $x_0$ génère une suite de points en sortie $y_1, y_2, \dots$ qui appartiennent à $R_y + \epsilon_y^{|\sigma_{i_0}|}, R_y + \epsilon_y^{|\sigma_{i_1}|}, \dots$ où $\sigma_{i_0}, \sigma_{i_1}, \dots$ sont les $k$-séquences successivement appliquées.*

# 4 Exemple d'illustration

## 4.1 Colonne de distillation

On considère un modèle linéarisé de dynamique chimique dans une colonne de distillation [6] permettant de se ramener à un système DLTI de la forme de l'équation (1). L'état $x = (x_1, x_2, \dots, x_{11})^\top$ est de dimension 11 : $x_1, x_2, \dots, x_{10}$ correspondent aux concentrations du mélange aux différents étages, et $x_{11}$ correspond à la pression en haut de la colonne. La perturbation $\omega$ en apport de produit est négligée. La variable de contrôle $u \in U = \{0, 1\}$, correspond à l'état allumé (1) ou éteint (0) du chauffage. La sortie $y$ est de dimension 1 et correspond à la concentration du composant le plus volatile dans le produit récupéré en bas de la colonne.

<figure><img src="image_4.png" /><figcaption>**FIGURE 2** – Modèle de la colonne de distillation de dimension 11 tiré de [6].</figcaption></figure>
Les matrices A, B et C sont les suivantes :

$$
A = 10^{-2} \times \begin{bmatrix} -1.4 & -0.43 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.95 & -1.38 & 0.46 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.05 \\ 0 & 0.95 & -1.41 & 0.63 & 0 & 0 & 0 & 0 & 0 & 0 & 0.02 \\ 0 & 0 & 0.95 & -1.58 & 1.1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.95 & -3.12 & 1.5 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2.02 & -3.52 & 2.2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 2.02 & -4.22 & 2.8 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 2.02 & -4.82 & 3.7 & 0 & 0.02 \\ 0.255 & 0 & 0 & 0 & 0 & 0 & 0 & 2.02 & -5.72 & 4.2 & 0.05 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2.55 & -1.85 \end{bmatrix},
$$

$$
B = [ 0 \ 0 \ 0 \ 0 \ 0.01 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0 ]^T,
$$

$$
C = [ 0 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0.01 \ 0 ]^T.
$$

Le système réduit par troncature équilibrée, d'ordre $n_r = 2$, est le suivant :

$$
\hat{A} = 10^{-3} \times \begin{bmatrix} -0.11 & -0.26 \\ 0.26 & -6.98 \end{bmatrix}, \quad \hat{B} = 10^{-3} \times \begin{bmatrix} 1.56 \\ -1.75 \end{bmatrix}, \quad \hat{C} = 10^{-3} \times \begin{bmatrix} 1.56 & 1.75 \end{bmatrix}.
$$

L'ensemble d'intérêt est $R_x = [0,1]^{11}$. L'ensemble cible est $R_y = [0.0015,0.0025]$. On calcule la borne $\epsilon = 8,25 \times 10^{-4}$. Les méthodes hors ligne et en ligne permettent d'obtenir le partitionnement en dix-neuf sous-ensembles Figure 3, les deux méthodes aboutissent à la même décomposition du fait du caractère contractif du modèle. Des simulations hors ligne et en ligne sont données Figure 4. Les marqueurs correspondent à la fin des séquences $\sigma_i$. La sortie $y_r$ du système réduit est stabilisée dans $R_y$. Dans le cas hors ligne, la sortie $y$ du système complet est stabilisée dans $R_y + \epsilon_y^\infty$. Dans le cas en ligne, la sortie $y$ est stabilisée dans $R_y + \epsilon_y^{\ell_i}$, où $\ell_i$ est la longueur de la $k$-séquence appliquée.

<figure><img src="image_5.png" /><figcaption>FIGURE 3 – Décomposition $\hat{\Delta}$ de $\hat{R}_x = \pi_r[0,1]^{11}$ dans le plan $(\hat{x}_1, \hat{x}_2)$.</figcaption></figure>

4.2 Autres exemples traités

Cette méthode a été testée sur d'autres exemples de la littérature, en particulier la régulation de la température d'un appartement quatre pièces venant de [7], et l'atteignabilité de la température d'une plaque chauffée [8]. Nous étudions actuellement des exemples de plus grande taille.

5 Remarques finales

Nous avons donné deux procédures permettant de synthétiser un contrôle qui atteint un objectif de stabilisation pratique dans une zone cible donnée dans le cas où le contrôle est quantifié. Le contrôle hors ligne permet, au prix d'une majoration $\epsilon_y^\infty$ de l'erreur relativement pessimiste, d'obtenir le même
<figure><img src="image_6.png" /><figcaption>FIGURE 4 – Simulation de $y(t) = Cx(t)$ et $y_r(t) = \hat{C}\hat{x}(t)$ à partir de $x_0 = (0.55)^{11}$. À gauche : contrôle hors ligne ; à droite : contrôle en ligne.</figcaption></figure>

comportement que le modèle réduit. Le contrôle en ligne conduit à des erreurs moins pessimistes mais implique un comportement qui peut être différent du modèle réduit, des décalages peuvent apparaître et les cycles limites (séquences appliquées en régime stationnaire) peuvent être différents. Le comportement du modèle complet est donc moins maîtrisé, mais son utilisation peut être effectuée en ligne.

Des questions restent ouvertes avec les méthodes proposées ici. En effet, lors d'une utilisation réelle, seule la sortie y du système complet $\Sigma$ est connue, ce qui implique de reconstruire l'état réduit $\hat{x} = \pi_r x$ en ligne à partir de y (soit en reconstruisant l'état complet puis en le projetant, ou en reconstruisant directement l'état projeté). Jusqu'ici, cette reconstruction était supposée exacte. Nos futurs travaux se dirigent vers cette étape de reconstruction en ligne, et notre principale piste est l'utilisation de filtres de Kalman étendus [9].

Références

[1] Bruno Picasso et Antonio Bicchi. On the stabilization of linear systems under assigned i/o quantization. *IEEE Trans. Automat. Contr.*, 52(10):1994–2000, 2007.

[2] Laurent Fribourg et Romain Soulat. Stability controllers for sampled switched systems. In Parosh Aziz Abdulla and Igor Potapov, editors, *Proceedings of the 7th Workshop on Reachability Problems in Computational Models (RP'13)*, volume 8169 of *Lecture Notes in Computer Science*, pages 135–145, Uppsala, Sweden, September 2013. Springer.

[3] Peter Benner, Jing-Rebecca Li, et Thilo Penzl. Numerical solution of large-scale lyapunov equations, riccati equations, and linear-quadratic optimal control problems. *Numerical Linear Algebra with Applications*, 15(9):755–777, 2008.

[4] Dan Sorensen et Athanasios Antoulas. Approximation of large-scale dynamical systems : an overview. *International Journal of Applied Mathematics and Computer Science*, 11(5):1093–1121, 2001.

[5] Zhi Han et Bruce H. Krogh. Reachability analysis of hybrid systems using reduced-order models. In *American Control Conference*, pages 1183–1189. IEEE, 2004.

[6] D. Tong, W. Zhou, A. Dai, H. Wang, X. Mou, et Y. Xu. $h_\infty$ model reduction for the distillation column linear system. *Circuits Syst Signal Process*, 2014.

[7] Alessandro Vittorio Papadopoulos et Maria Prandini. Model reduction of switched affine systems : a method based on balanced truncation and randomized optimization. In *HSCC'14*, pages 113–122, 2014.

[8] Zhi Han et Bruce H. Krogh. Reachability analysis of large-scale affine systems using low-dimensional polytopes. In JoÃ£oP. Hespanha and Ashish Tiwari, editors, *Hybrid Systems : Computation and Control*, volume 3927 of *Lecture Notes in Computer Science*, pages 287–301. Springer Berlin Heidelberg, 2006.

[9] Bruce P. Gibbs. *Advanced Kalman filtering, least-squares and modeling : a practical handbook*. John Wiley & Sons, 2011.