# Chapitre 3

## Diagnosticabilité fonctionnelle et
algébrique

Depuis 2012, je m'intéresse aux problèmes de diagnostic suite au projet ANR MAGIC-SPS [102, 198, 202]. Le diagnostic, terme souvent entendu dans le milieu médical, consiste à identifier la cause (l'origine) d'une défaillance, d'un problème, à partir de symptômes<fnref n="1" /> relevés par des observations ou des tests du système physique. Sur la plupart des systèmes physiques, le diagnostic peut se faire à partir d'un modèle mathématique. Pour cela, nous faisons les hypothèses suivantes :

— le comportement (ou une partie) du système physique peut être modélisé par un modèle mathématique,

— les mesures des capteurs correspondent à des sorties du modèle,

— les fautes du système physique peuvent être représentées par des variables dans le modèle mathématique. Nous entendons par (multi) fautes un changement d'une (ou de plusieurs valeurs de) paramètres impliquant des dysfonctionnements du système. Nous supposons également qu'il n'y a pas de restriction sur le type de faute. Elles peuvent agir multiplicativement ou additivement sur les paramètres.

Le diagnostic d'un système est défini comme la détection et l'isolation des fautes (ou la localisation et l'identification). Étant donné le système et le modèle mathématique associé [84], la détection de fautes consiste à détecter un mauvais fonctionnement du système et l'isolation de fautes à isoler la faute responsable du dysfonctionnement.

Diverses méthodes peuvent être trouvées dans la littérature pour faire du diagnostic. Elles dépendent essentiellement du type de modèle choisi. Dans le cas de modèles non linéaires, les méthodes classiques pour faire le diagnostic sont basées sur les observateurs non linéaires ([180] par exemple) et/ou sur le relations analytiques redondantes (en anglais Analytical Redundancy Relations, ARRs) ([183, 187]). Ces dernières sont des relations liant les entrées, les sorties, leurs dérivées, les paramètres et les fautes du modèle (voir [53, 54, 183, 198, 225] pour les simples fautes et [46, 58, 113] pour des fautes multiples). Ainsi, l'étude du diagnostic se fait à partir de signatures de fautes qui, classiquement, sont des fonctions qui à un ensemble de fautes associe un ensemble d'indicateurs ou de relations appelés résidus.

Dans notre travail, nous avons supposé que la représentation mathématique du système physique se faisait par un modèle paramétrique non linéaire et que l'on disposait de signaux d'entrée et de sortie. Le choix du modèle amène naturellement à se poser la question :

*est-ce que le modèle mathématique que nous avons choisi peut discriminer d'éventuelles fautes à partir des données collectées sur le système physique ?*

Cette étude a priori du modèle consiste en une étude de diagnosticabilité. En analysant la diagnosticabilité, il est possible d'anticiper le pouvoir discriminant d'un appareil de diagnostic au moment de sa mise en marche et de proposer des solutions à d'autres problèmes importants comme celui du nombre minimal de capteurs et de leur place pour garantir la discriminabilité d'un ensemble anticipé de défauts.

Au cours de ce travail de recherche, nous avons abordé ces études de diagnosticabilité et de diagnostic en exploitant de nouvelles ARRs obtenues à partir de l'algèbre différentielle. L'originalité a été la façon d'exploiter ces relations ce qui nous a permis de développer de nouveaux algorithmes certifiés.

<footnote n="1">1. Les symptômes peuvent être définis comme les écarts constatés entre des valeurs mesurées ou estimées à partir d'un processus et des valeurs présupposées.</footnote>
## 3.1 Le modèle

Le modèle considéré dans ce chapitre est le suivant :

$$ egin{cases} 
\dot{x}(t,p,f) = g(x(t,p),u(t),f,p), \\ 
y(t,p,f) = h(x(t,p),u(t),f,\varepsilon(t),p), \\ 
x(t_0,p,f) = x_0, \\ 
t_0 \le t \le T. 
\end{cases} 
\quad (3.1) $$

où

- $x(t,p,f) \in \mathbb{R}^n$ et $y(t,p,f) \in \mathbb{R}^m$ représentent les variables d'état et les sorties respectivement,

- les fonctions $g$ et $h$ sont réelles et analytiques sur $M$, $M$ étant un ensemble ouvert de $\mathbb{R}^n$ tel que $x(t,p,f) \in M$ pour chaque $t \in [t_0, T]$,

- $u(t) \in \mathbb{R}^r$ est le vecteur de contrôle,

- $f \in \mathbb{R}^e$ est le vecteur de fautes,

- le vecteur de paramètres $p$ appartient à $U_P$, où $U_P \subseteq \mathbb{R}^q$ est l'ensemble des paramètres admissibles à priori connu,

- les conditions initiales $x_0$ sont supposées appartenir à un ensemble borné $X_0$, et être différents d'un point d'équilibre du système.

$f=0$ signifie qu'il n'y a pas de faute. Dans le cas d'un modèle non contrôlé, $u=0$.

Les auteurs de [53] considèrent qu'un système est diagnosticable si $f$ est algébriquement observable en $u$ et $y$. Si $f_i$ est la ième composante du vecteur de faute $f$, cela signifie que chaque composante $f_i$ peut s'écrire comme la solution d'une équation polynomiale en $f_i$ et des entrées $u$, des sorties $y$ et de leurs dérivées. Cette définition est équivalente à la définition d'identifiabilité proposée par [125]. Rappelons que dans [125], les polynômes sont obtenus à partir de l'algorithme de Ritt qui cherche à exprimer chaque paramètre/faute en fonction des entrées, sorties et de leurs dérivées. Dans [225], les auteurs ont proposé une méthode de détection et d'isolation de faute pour des systèmes non linéaires. Ils ont basé leur travail de diagnostic sur les relations obtenues par [125] conscients de la difficulté de les obtenir. Ces relations sont utilisées pour construire des résidus évalués grâce à une méthode de détection de changements statistiques. Toutefois, ils ne considèrent pas l'analyse de diagnosticabilité ni les conditions pour que les fautes soient discriminables.

Dans ce travail, nous sommes partis des relations entrées-sorties et les avons exploitées pour proposer une étude de diagnosticabilité, de détection et d'isolation de fautes. Elles ont été exploitées de deux façons différentes : i) en les considérant comme relations fonctionnelles; ii) en utilisant leurs coefficients en tant que relations algébriques.

Dans un premier temps, les relations polynomiales ES sont ré-écrites en reprenant le polynôme (2.4). Cette ré-écriture a permis dans un second temps de définir des signatures qui ont servi aux études de diagnosticabilité et de détectabilité.

Si on note $\theta$ le vecteur $(f,p)^T$, en reprenant (2.4), on obtient le polynôme :

$$ P_i(y, u, f, p) = m_{0,i}(y, u, p) - \sum_{k=1}^{s_i} \gamma_k^i(f, p) m_{k,i}(y, u) = 0. \quad (3.2) $$

La première partie du polynôme $\rho_i = m_{0,i}(y, u, p)$ est appelée résidu. Celle-ci ne contient pas de composantes de $f$. Sans faute ($f=0$), cette relation qui est calculable en les points de mesure est nulle. La seconde partie du polynôme (3.2), $w_{1,i}(y, u, f, p) := \sum_{k=1}^{s_i} \gamma_k^i(f, p) m_{k,i}(y, u)$, est appelée la forme résiduelle interne.

**Exemple 3.1.1.** Considérons une masse $m=1$ attachée à un ressort de force $k$. Notons $u$ une force extérieure au système non identiquement nulle et $d \ge 1$ une constante. Le mouvement de la masse est décrit par l'équation suivante

$$ \ddot{x} + kx - du = 0. \quad (3.3) $$

Supposons que deux fautes $f_1 \in [0,2)$ et $f_2 \in [0,2)$ impactent respectivement le ressort $k$ et le paramètre $d$ et que le modèle (3.3) prenne la forme :

$$ \ddot{x} + k(f_1 - 1)^2 x - (d + f_2)u = 0. \quad (3.4) $$

Supposons que $y := x$ est une sortie mesurée du système. L'équation (3.4) représente par définition une ARR du modèle (3.3). On peut la ré-écrire :

$$ P(y, u, f, p) = \ddot{y} + k(f_1 - 1)^2 y - (d + f_2)u = 0. \quad (3.5) $$
Remarquons que l'équation (3.4) peut se mettre sous la forme du modèle (3.1), ce qui donne :

$$ egin{cases} \dot{x}_1 = x_2, \\ \dot{x}_2 = -k(f_1 - 1)^2 x_1 + (d + f_2)u \end{cases} \tag{3.6} $$

et qu'on peut retrouver (3.5) en utilisant l'algorithme de Rosenfeld-Groebner.

Ainsi, l'ARR (3.5) peut se ré-écrire

$$ P(y, u, f, p) = m_0(y, u, p) - w_1(y, u, f, p) = 0 \tag{3.7} $$

$$ \text{avec } m_0(y, u, p) = \ddot{y} + k y - du \text{ et } w_1(y, u, p) = -k(f_1^2 - 2f_1)y + f_2 u. $$

La première définition de signature qui est apparue est le vecteur de booléen associant la valeur 1 si la faute apparaît dans l'ARR et 0 sinon. Dans cet exemple, cette signature est définie par $Sig(f_1) = Sig(f_2) = 1$ car les 2 fautes apparaissent dans la même ARR. Cette signature ne permet donc pas de discriminer quelle faute agit sur le système. Nous avons eu l'idée de définir d'autres signatures à partir de cette ARR en exploitant ses propriétés.

À partir de la forme interne du résidu $w_1$, nous avons proposé les deux définitions suivantes. La première est une signature fonctionnelle définie dans cet exemple par $FSig(f) = (-k(f_1^2 - 2f_1)y, f_2u)$. La première composante, obtenue en annulant $f_2$, caractérise la faute $f_1$ et la deuxième composante la faute $f_2$. Si les deux composantes sont distinctes alors il est possible de discriminer les deux fautes. La seconde signature est une signature algébrique définie à partir des coefficients de $w_1$. La signature algébrique prend la forme $ASig(f) = (k(f_1 - 1)^2, -d - f_2)$.

Nous allons voir dans la suite comment nous avons défini des critères pour discriminer les fautes à partir des signatures.

## 3.2 Diagnosticabilité fonctionnelle

### 3.2.1 Définition de diagnosticabilité basée sur la signature fonctionnelle

Soit $f_j$ la jème composante du vecteur de faute $f$. Notons $f_{[j]}$ le vecteur de faute dont toutes les composantes sont nulles sauf la jème égale à $f_j$. $f_j$ comme $f_{[j]}$ fait référence à la simple faute.

Supposons que $m$ résidus peuvent être obtenus. La définition classique de la signature de la faute $f_j$ est la suivante [187] :

**Définition 16.** La signature de $f_j$ est le vecteur $Sig(f_j)$ de dimension $m$ dont la ième composante est égale à 1 si le ième résidu $w_{1,i}(y,u,f,p)$ contient la faute $f_j$ et est égale à 0 sinon.

À partir de cette définition, les auteurs de l'article [187] ont proposé les définitions de discriminalité et de diagnosticabilité suivantes.

**Définition 17.** Deux fautes $f_i$ et $f_j$ sont discriminables si leurs signatures sont différentes. Quand toutes les fautes sont discriminables, le modèle est dit diagnosticable.

D'après cette définition, si deux fautes ont les mêmes résidus, le modèle n'est pas diagnosticable comme dans l'exemple 3.1.1.

Cette définition est loin d'exploiter les propriétés des résidus. Aussi, nous avons introduit la notion de signature fonctionnelle.

La signature fonctionnelle de la faute $f_j$ est donnée par la forme interne des $m$ résidus en considérant que seule la faute $f_j$ agit sur le système. En d'autres termes, la signature fonctionnelle est définie par le vecteur dont la ième composante est $w_{1,i}(\bar{y}, u, f_{[j]}, p)$, c'est à dire le polynôme obtenu à partir de $w_{1,i}(y, u, f, p)$ en considérant que toutes les composantes de $f$ sont nulles sauf la jème égale à $f_j$. Formellement, la définition est la suivante :

**Définition 18.** La signature fonctionnelle est la fonction $FSig$ qui, à la faute $f_j$ associe le vecteur $(w_{1,i}(\bar{y}, \bar{u}, f_{[j]}, p))_{i=1,...,m}$.

L'avantage de cette définition est qu'elle prend en compte explicitement l'effet de l'amplitude du défaut et de la valeur de l'entrée. Elle est donc plus informative que la définition classique.

Soit $FSig^{(i)}(f_j) = w_{1,i}(y, u, f_{[j]}, p)$ la ième composante de $FSig(f_j)$. La question qui se pose est maintenant de définir la notion de discriminalité par rapport à cette définition fonctionnelle.
**Exemple 3.2.1.** Considérons l'équation de Bernoulli :

$$ \dot{y}(t) = \beta_1 y(t) + \beta_2 y(t)^2, \quad \text{for } t \in [0,5], y(0) = -1. \tag{3.8} $$

Sa solution est $y(t) = \frac{\beta_1 e^{\beta_1 t}}{\beta_2 - \beta_1 - \beta_2 e^{\beta_1 t}}$.

Nous obtenons directement $w_{0,1}(y,p) = \dot{y}(t) - \beta_1 y(t) - \beta_2 y(t)^2$.

Supposons que deux simples fautes positives $f_1$ et $f_2$ impactent additivement les deux paramètres $\beta_1$ et $\beta_2$ respectivement. Alors $w_{1,1}(y,f,p) = f_1 y(t,p,f) + f_2 y(t,p,f)^2$ où $p = (\beta_1, \beta_2)^T$, $f = (f_1, f_2)^T$, et $T$ représente la transposée du vecteur considéré. On en déduit que $\text{Sig}(f_1) = \text{Sig}(f_2) = 1$ et $\text{FSig}(f_1) = f_1 y(\cdot, p, f_{[1]})$, $\text{FSig}(f_2) = f_2 y(\cdot, p, f_{[2]})^2$ où $y(\cdot, p, f_{[i]})$ représente la sortie impactée par la faute $f_i$ seulement. Là encore, les deux fautes ne sont pas discriminables au sens classique car leurs signatures sont égales quelque soit les valeurs de $\beta_1$ et $\beta_2$. Cependant, elles peuvent être discriminables au sens de la définition de la signature fonctionnelle pour $t \in [0,5]$. Par exemple, avec les valeurs de paramètres $\beta_1 = 1, \beta_2 = 2$, on peut montrer que $y(t,p,f_{[1]})$ est négatif sur $[0,5]$ donc $\text{FSig}(f_1)$ également alors que $\text{FSig}(f_2)$ est positif.

Suite à cet exemple, nous avons proposé les définitions suivantes pour lier les notions de discriminalité et de diagnosticabilité fonctionnelle. La première est vraie pour toute entrée, la seconde est vérifiée seulement pour une seule entrée.

**Définition 19.** Deux fautes $f_j$ et $f_l$ sont fortement fonctionnellement discriminables si pour toute entrée $u$, il existe au moins un indice $i^*$ et un temps fini $t_1 \in [t_0, T]$ tel que pour tout $t \in [t_0, t_1]$, $\text{FSig}^{(i^*)}(f_j)(t) \neq \text{FSig}^{(i^*)}(f_l)(t)$. Quand toutes les fautes sont fortement fonctionnellement discriminables, le modèle est dit fortement diagnosticable.

**Définition 20.** Deux fautes $f_j$ et $f_l$ sont faiblement fonctionnellement discriminables si il existe une entrée $u$, il existe au moins un indice $i^*$ et un temps fini $t_1 \in [t_0, T]$ tel que pour tout $t \in [t_0, t_1]$, $\text{FSig}^{(i^*)}(f_j)(t) \neq \text{FSig}^{(i^*)}(f_l)(t)$. Quand toutes les fautes sont faiblement fonctionnellement discriminables, le modèle est dit faiblement diagnosticable.

La définition peut-être étendue à des modèles non contrôlés.

La détectabilité est un cas particulier de discriminabilité qui nécessite que la situation sans défaut soit discriminable quelque soit la faute. La détectabilité fonctionnelle peut donc être définie sur la base des signatures fonctionnelles comme suit.

**Définition 21.** La faute $f_i$ est fonctionnellement détectable si la signature fonctionnelle $\text{FSig}(f_i)$ n'est pas égale au vecteur nul.

**Remarque 21.** Comme nous le verrons dans l'exemple 3.2.3, contrairement à la signature classique, la signature fonctionnelle apporte des informations sur la sensibilité des résidus.

Les signatures fonctionnelles de fautes peuvent être regroupées dans une matrice dont la composante $(i,j)$ contient le ième résidu en la faute $f_j$, c'est à dire $w_{1,i}(\bar{y}, \bar{u}, f_{[j]}, p)$.

**Exemple 3.2.2.** Considérons l'exemple suivant :

$$ \left\{ \begin{array}{l} \dot{x}_1 = (p_1 + f_1)(p_2 + f_2)x_1^2 + x_1x_2, \\ \dot{x}_2 = (p_2 + f_2)(p_3 + f_3)x_2^2 + x_2x_3, \\ \dot{x}_3 = (p_1 + f_1)(p_3 + f_3)x_3^2 + x_1x_3, \\ y_1 = x_1, y_2 = x_2, y_3 = x_3. \end{array} \right. \tag{3.9} $$

Il est facile de vérifier que :

$$ \begin{align} w_{0,1}(\bar{y}, \bar{u}, p) &= \dot{y}_1 - y_1 y_2 - p_1 p_2 y_1^2, \notag \\ w_{0,2}(\bar{y}, \bar{u}, p) &= \dot{y}_2 - y_2 y_3 - p_2 p_3 y_2^2, \tag{3.10} \\ w_{0,3}(\bar{y}, \bar{u}, p) &= \dot{y}_3 - y_1 y_3 - p_1 p_3 y_2^2. \notag \end{align} $$

et :

$$ \begin{align} w_{1,1}(\bar{y}, \bar{u}, f, p) &= (p_1 f_2 + p_2 f_1 + f_1 f_2)y_1^2, \notag \\ w_{1,2}(\bar{y}, \bar{u}, f, p) &= (p_2 f_3 + p_3 f_2 + f_2 f_3)y_2^2, \tag{3.11} \\ w_{1,3}(\bar{y}, \bar{u}, f, p) &= (p_1 f_3 + p_3 f_1 + f_1 f_3)y_3^2. \notag \end{align} $$
Les signatures fonctionnelles sont :

$$
\begin{align*}
FSig(f_1) &= (p_2 f_1 y_1^2, 0, p_3 f_1 y_3^2)^T, \\
FSig(f_2) &= (p_1 f_2 y_1^2, p_3 f_2 y_2^2, 0)^T, \\
FSig(f_3) &= (0, p_2 f_3 y_2^2, p_1 f_3 y_3^2)^T.
\end{align*}
$$

La diagnosticabilité peut alors être analysée grâce à la table suivante :

<table>
  <thead>
    <tr>
      <th>FSig(f)</th>
      <th>f<sub>1</sub></th>
      <th>f<sub>2</sub></th>
      <th>f<sub>3</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>FSig<sup>(1)</sup>(f)</th>
      <td>p<sub>2</sub>f<sub>1</sub>y<sub>1</sub><sup>2</sup></td>
      <td>p<sub>1</sub>f<sub>2</sub>y<sub>1</sub><sup>2</sup></td>
      <td>0</td>
    </tr>
    <tr>
      <th>FSig<sup>(2)</sup>(f)</th>
      <td>0</td>
      <td>p<sub>3</sub>f<sub>2</sub>y<sub>2</sub><sup>2</sup></td>
      <td>p<sub>2</sub>f<sub>3</sub>y<sub>2</sub><sup>2</sup></td>
    </tr>
    <tr>
      <th>FSig<sup>(3)</sup>(f)</th>
      <td>p<sub>3</sub>f<sub>1</sub>y<sub>3</sub><sup>2</sup></td>
      <td>0</td>
      <td>p<sub>1</sub>f<sub>3</sub>y<sub>3</sub><sup>2</sup></td>
    </tr>
  </tbody>
</table>

Clairement, pour $(j,l) \in \{1,2,3\}^2$, $j \neq l$, $FSig(f_j) \neq FSig(f_1)$ et le modèle est fonctionnellement diagnosticable.
Le système (3.9) est également diagnosticable au sens classique.

3.2.2 Lien entre identifiabilité et diagnosticabilité fonctionnelle

Rappelons que

$$
FSig(f_j) = \left( \sum_{k=1}^{s_i} \gamma_k^i(f_j, p) m_{k,i}(y, u) \right)_{i=1,\dots,m} .
$$

**Proposition 20.** Supposons qu'il existe un ensemble d'indices $I_0 \subseteq \{1, \dots, m\}$ tel que pour tout $i \in I_0$, les déterminants fonctionnels $\Delta w_{1,i}(y, u) = \det(m_{k,i}(y, u), k = 1, \dots, s_i)$ sont non identiquement nulles. Si, avec les polynômes $w_i(y, u, f, p)$, $i \in I_0$, le modèle est identifiable en les fautes alors le modèle est faiblement fonctionnellement diagnosticable. La réciproque n'est pas vraie.

De la proposition 1, on en déduit que l'injectivité de la fonction $\phi(\cdot)$ est une condition suffisante pour que le modèle soit diagnosticable mais elle n'est pas nécessaire. En effet, considérons la simple équation différentielle $\dot{y}(t,p) = (p_1 - p_2)y(t)$, $y(0) = 1$, $p = (p_1, p_2)^T$, dont la solution est $y(t,p) = e^{(p_1-p_2)t}$, $t \in \mathbb{R}$. Supposons que deux simples fautes additives et strictement positives $f_1$ et $f_2$ agissent sur les deux paramètres $p_1$ et $p_2$ respectivement. On déduit directement de l'équation que $w_1(\bar{y}, f, p) = (f_1 - f_2)y(t, p, f)$. Ce modèle n'est pas identifiable en les fautes car la fonction $\phi: (f_1, f_2) \mapsto f_1 - f_2$ n'est pas injective. Cependant, la signature fonctionnelle $FSig(f_1) = f_1y(t,p, f_{[1]}) = f_1e^{(p_1+f_1-p_2)t}$ est strictement positive et la signature fonctionnelle $FSig(f_2) = -f_2y(t,p, f_{[2]}) = -f_2e^{(p_1-p_2-f_2)t}$ est strictement négative. Les deux signatures sont donc distinctes sur $\mathbb{R}$ et le modèle est fonctionnellement diagnosticable.

Les signatures fonctionnelles dépendent directement des fautes a priori inconnues. Toutefois, il est possible
de ré-écrire ces expressions en fonction du résidu. En dehors de toute faute, le polynôme P (voir (3.2) dans le
cas où on a un seul polynôme car alors l'indice i est omis) est égal au résidu ρ = m₀(y, u, p). Lorsqu'une faute
agit sur le système, ce polynôme n'est plus nul et est égal à la signature fonctionnelle définie par la faute. Ainsi,
selon la valeur du résidu et des sorties mesurées, il est possible de détecter et d'isoler la faute.
Reprenons l'exemple 3.2.1 sur l'équation de Bernoulli.

**Exemple 3.2.3.** Pour l'exemple 3.2.1, le résidu vaut $\rho = w_{0,1}(\bar{y},p) = \dot{y}(t) - \beta_1 y(t) - \beta_2 y(t)^2$. Comme précédem-
ment, nous supposons que deux simples fautes dans $[0,1]$ agissent additivement sur les deux paramètres $\beta_1$ et $\beta_2$
respectivement. Rappelons que $w_{1,1}(\bar{y},f,p) = f_1 y(t,p,f) + f_2 y(t,p,f)^2$ et $FSig(f_1) = f_1 y(\cdot,p, f_{[1]})$, $FSig(f_2) =$
$f_2 y(\cdot,p, f_{[2]})^2$. Le déterminant fonctionnel est égal à $\Delta w_{1,1}(\bar{y}) = \dot{y}^2 \neq 0$ et la fonction $\phi: (f_1, f_2) \mapsto (f_1, f_2)$ est
injective. D'après les propositions 1 et 20, le modèle est globalement identifiable en les fautes et donc fonction-
nellement identifiable sur $[0,5]$. Ainsi, les signatures fonctionnelles sont distinctes lorsqu'une des fautes agit.
Numériquement, supposons que $\beta_1 = 1$, $\beta_2 = 2$, et qu'une des deux fautes agit de façon permanente à partir de
t = 0.5s. Supposons que la sortie soit perturbée par un bruit gaussien telle que l'erreur relative a une valeur
maximale de 0.1. Les figures 3.12 et 3.13 représentent les résidus pour différents scénarios consistant à faire
varier la valeur des fautes entre 0.1 et 1. On peut remarquer que l'évaluation des signatures à travers les résidus
permet de détecter les fautes agissant sur le système car leurs comportements sont différents.

**Remarque 22.** Remarquons que la signature fonctionnelle donne un moyen de détecter une faute et de l'isoler.
Lorsqu'il n'y a pas de fautes agissant sur le système, la signature fonctionnelle est nulle. Cependant, la présence
de bruit de mesures corrompt le signal mesuré conduisant à une signature fonctionnelle non nulle. Nous avons
donc introduit la notion de $\varepsilon$-fonctionnellement détectable. La faute est dite $\varepsilon$-fonctionnellement détectable s'il
<figure><img src="image_24.png" /><figcaption>FIGURE 3.1 – Valeur des signatures obtenues à partir de l’équation de Bernoulli quand une faute de magnitude comprise entre 0.1 et 1 agit sur $\beta_1$ à partir de t=0,5s.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 3.2 – Valeur des signatures obtenues à partir de l’équation de Bernoulli quand une faute de magnitude comprise entre 0.1 et 1 agit sur $\beta_2$ à partir de t=0,5s.</figcaption></figure>
existe un intervalle de temps $[t_1, t_2] \subseteq [t_0, T]$ tel que pour tout $t \in [t_1, t_2]$, $d(FSig^{(j)}(f_i), 0) > \varepsilon$, avec d'une
distance sur $\mathbb{R}$. Cette condition assure qu'au moins un des résidus est au dessus du seuil $\varepsilon$.
Pour l'exemple de Bernoulli ci-dessus, on constate sur la figure 3.13 que les fautes sont 0.02-fonctionnellement
détectables.

3.2.3 Un exemple complet : celui du réservoir

Considérons l'exemple de deux réservoirs couplés

$$
\left\{
\begin{array}{l}
\dot{x}_1(t,p) = a_1 u(t) - a_2 \sqrt{x_1(t,p)}, x_1(0) = 1, \\
\dot{x}_2(t,p) = a_3 \sqrt{x_1(t,p)} - a_4 \sqrt{x_2(t,p)}, x_2(0) = 0.6, \\
y_1(t,p) = a_5 \sqrt{x_1(t,p)}, \\
y_2(t,p) = a_6 \sqrt{x_2(t,p)},
\end{array}
\right.
\qquad (3.12)
$$

où $p = (a_i)_{i=1,...,6}, a_i \neq 0$, est le vecteur de paramètres du modèle, $x = (x_1, x_2)^T$ représente le vecteur d'état et correspond au niveau de l'eau dans chacun des réservoirs, $u \neq 0$ est l'entrée. Supposons que le niveau d'eau dans les réservoirs peut varier entre 0 et 10.

Soient $f_1$ une faute inconnue agissant additivement sur le signal d'entrée, $f_2$ et $f_3$ des fautes additives agissant sur les deux capteurs en sortie de chacun des réservoirs d'eau, et $f_4$ un défaut de colmatage. $f_4 = 1$ représente un tuyau complètement bouché et $0 < f_4 < 1$ un tuyau bouché partiellement. Pour le scénario avec faute, nous supposons qu'elles sont introduites au temps $t = 20s$.

Pour utiliser l'algorithme de Rosenfeld-Groebner implémenté dans Maple 16, les variables auxiliaires $z_1(t,p) = \sqrt{x_1(t,p)}$ et $z_2(t,p) = \sqrt{x_2(t,p)}$ sont introduites dans le modèle ainsi que les 4 fautes. Le modèle se réécrit :

$$
\left\{
\begin{array}{l}
\dot{x}_1(t,p) = a_1 (u(t) + f_1) - a_2 (1 - f_4) z_1(t,p), \\
\dot{x}_2(t,p) = a_3 (1 - f_4) z_1(t) - a_4 z_2(t,p), \\
z_1(t,p)^2 = x_1(t,p), \\
z_2(t,p)^2 = x_2(t,p), \\
y_1(t,p) = a_5 (1 - f_4) z_1(t,p) + f_2, \\
y_2(t,p) = a_6 z_2(t,p) + f_3, \\
\dot{f}_i = 0.
\end{array}
\right.
\qquad (3.13)
$$

L'algorithme de Rosenfeld-Groebner fournit les deux ARRs :

$$
\begin{align}
w_1(y, u, f, p) &= w_{0,1}(y, u, p) - w_{1,1}(y, u, f, p) \tag{3.14} \\
w_2(y, u, f, p) &= w_{0,2}(y, u, p) - w_{1,2}(y, u, f, p) \nonumber
\end{align}
$$

où

$$
\begin{align*}
w_{0,1} &= -ua_1 a_5^2 + (a_2 a_5 + 2 \dot{y}_1) y_1, \\
w_{1,1} &= 2 \dot{y}_1 f_2 - (f_4^2 a_2 a_5 - 2 f_4 a_2 a_5) y_1 - (-f_4^2 a_1 a_5^2 + 2 f_4 a_1 a_5^2) u \\
&\quad + f_4^2 f_2 a_2 a_5 + f_4^2 f_1 a_1 a_5^2 - 2 f_4 f_2 a_2 a_5 - 2 f_4 f_1 a_1 a_5^2 + f_2 a_2 a_5 + f_1 a_1 a_5^2, \tag{3.15} \\
w_{0,2} &= 2 a_5 \dot{y}_2 y_2 - a_3 a_6^2 y_1 + a_4 a_5 a_6 y_2, \\
w_{1,2} &= 2 \dot{y}_2 f_3 a_5 + f_3 a_4 a_5 a_6 - f_2 a_3 a_6^2.
\end{align*}
$$

D'où,

$$
\begin{align*}
FSig(f_1) &= (-f_1 a_1 a_5^2, 0)^T \\
FSig(f_2) &= (-2 \dot{y}_1 f_2 - f_2 a_2 a_5, f_2 a_3 a_6^2)^T \\
FSig(f_3) &= (0, -2 \dot{y}_2 f_3 a_5 - f_3 a_4 a_5 a_6)^T \\
FSig(f_4) &= (y_1 f_4^2 a_2 a_5 - 2 y_1 f_4 a_2 a_5 - u f_4^2 a_1 a_5^2 + 2 u f_4 a_1 a_5^2, 0)^T.
\end{align*}
\tag{3.16}
$$
Les déterminants fonctionnelles $\Delta w_{1,1}(\bar{y}, \bar{u}) = \det(1, u, y_1, \dot{y}_1) = \dot{u}(\ddot{y}_1 y_1^{(4)} - y_1^{(3)^2}) - \ddot{u}(\dot{y}_1 y_1^{(4)} - y_1^{(3)} \dot{y}_1) + u^{(3)}(\dot{y}_1 y_1^{(3)} - \dot{y}_1^2)$ et $\Delta w_{1,2}(\bar{y}, \bar{u}) = \det(\dot{y}_2, 1) = -\dot{y}_2$ ne sont pas identiquement nuls tant que $u$ n'est pas constant. La fonction

$$ \phi(f) = (-2f_2, f_4^2 a_2 a_5 - 2f_4 a_2 a_5, -f_4^2 a_1 a_5^2 + 2f_4 a_1 a_5^2, -f_4^2 f_2 a_2 a_5 - f_4^2 f_1 a_1 a_5^2 + 2f_4 f_2 a_2 a_5 \\ +2f_4 f_1 a_1 a_5^2 - f_2 a_2 a_5 - f_1 a_1 a_5^2, -2f_3 a_5, -f_3 a_4 a_5 a_6 + f_2 a_3 a_6^2) \quad (3.17) $$

est clairement injective pour $f_4 \in ]0, 1[$. On en déduit que le modèle est globalement identifiable en $f$ lorsque $f \in ]0, 1[$. D'après la proposition 20, le modèle est faiblement fonctionnellement diagnosticable.

Dans les simulations, un simple contrôleur est utilisé pour contrôler le niveau de l'eau dans le réservoir supérieur. Il suit un signal carré de référence. Les deux capteurs sont supposés être perturbés par un bruit Gaussien de telle sorte que l'erreur relative a une valeur maximale de 0.1. Les paramètres du modèle sont égaux à $a_1 = a_2 = a_3 = a_4 = 0.3$, $a_5 = a_6 = 1$. La figure 3.3 montre les niveaux d'eau dans les deux réservoirs dans le cas sans faute mais avec les données bruitées. Les résidus déduits de (3.14) sont $\rho_1 = -ua_1a_5^2 + (a_1a_5 + 2\dot{y}_1)y_1$

<figure><img src="image_26.png" /><figcaption>FIGURE 3.3 – Niveau d'eau dans le réservoir supérieur, $y_1$, et inférieur $y_2$ lors des simulations sans faute.</figcaption></figure>

et $\rho_2 = 2a_5\dot{y}_2y_2 - a_3a_6^2y_1 + a_4a_5a_6y_2$. Ils sont utilisés pour détecter et discriminer les fautes, celles-ci étant introduites à $t = 20s$. Les dérivées sont estimées en utilisant les différenciateurs HOSM ([123], [124]). Les figures 3.4-3.11 représentent les résidus $\rho_1$ et $\rho_2$ quand les différentes fautes agissent sur le système.

<figure><img src="image_27.png" /><figcaption>FIGURE 3.4 – Résidu $\rho_1$ quand la faute $f_1$ est introduite au temps $t = 20s$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIGURE 3.5 – Résidu $\rho_2$ quand la faute $f_1$ est introduite au temps $t = 20s$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>FIGURE 3.6 – Résidu $\rho_1$ quand la faute $f_2$ est introduite au temps $t = 20s$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE 3.7 – Résidu $\rho_2$ quand la faute $f_2$ est introduite au temps $t = 20s$.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>FIGURE 3.8 – Résidu $ρ_1$ quand la faute $f_3$ est introduite au temps $t = 20s$.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>FIGURE 3.9 – Résidu $ρ_2$ quand la faute $f_3$ est introduite au temps $t = 20s$.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>FIGURE 3.10 – Résidu $ρ_1$ quand la faute $f_4$ est introduite au temps $t = 20s$.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>FIGURE 3.11 – Résidu $ρ_2$ quand la faute $f_4$ est introduite au temps $t = 20s$.</figcaption></figure>

D'après les graphes des résidus, chaque faute est 0.03-détectable et il est possible de déduire la faute qui agit.

Remarquons que le résidu de la faute $f_1$ est positif contrairement à la faute $f_2$ qui est négatif et que le cas a priori ambigu est celui des fautes $f_1$ et $f_4$ car elles ont la même signature classique : la première composante est non nulle alors que la deuxième est nulle. Toutefois, le comportement des graphes est complètement différent. Par exemple, sur l'intervalle de temps [21, 29], elles sont 0.04-fonctionnellement discriminable, ce qui assure une certaine sécurité sur le résultat quant au bruit.

Cet exemple montre une certaine robustesse des résultats. Toutefois, ces résultats sont dépendants de l'identifiabilité des fautes et le cas multi-fautes lorsqu'il n'y a pas identifiabilité n'a pu être résolu pour le moment avec les signatures fonctionnelles. L'approche algébrique présentée à la section suivante nous a permis d'apporter une solution.

## 3.3 Diagnosticabilité algébrique

Dans [202], nous avons proposé une nouvelle approche pour étudier la diagnosticabilité et faire du diagnostic. Celle-ci consiste à exploiter les coefficients des ARRs qui sont vus comme des relations algébriques. Ce changement de point de vue et l'utilisation des outils semi-algébriques a permis tout d'abord d'intégrer des contraintes sur les paramètres et les fautes, comme les inégalités satisfaites par les paramètres et les contraintes déduites des conditions initiales. Comme nous le verrons par la suite, ces contraintes peuvent jouer un rôle fondamental dans l'analyse de la diagnosticabilité et le diagnostic. Ensuite, cette approche ne nécessite pas d'hypothèse forte sur le modèle comme i) son identifiabilité ii) la valeur de certains paramètres dans des cas particuliers iii) l'action additive des (multi)-fautes.

Cette méthode se résume ainsi :

1. calculer les ARRs en utilisant l'algorithme de Rosenfeld-Groebner sur le modèle;

2. utiliser les bases de Groebner pour obtenir une application algébrique que nous appellerons *signature algébrique*. Chacune de ses composantes dépend seulement des paramètres et des coefficients des ARRs. Elle peut être estimée numériquement à partir des entrées et des sorties du système. Par construction, chaque composante de la signature algébrique s'annule quand au moins une des (multi)-faute agit sur le système;

3. utiliser pour chaque (multi-)faute possible des outils semi-agébriques pour vérifier que des composantes de la signature s'annulent ou pas. Cette étape permet de certifier le résultat. Les résultats peuvent être résumés dans un tableau pré-calculé;
4. utiliser le tableau comme entrée d'un traitement numérique : à partir des mesures du système, il retourne les estimations des signatures algébriques. Leur comparaison avec les valeurs nominales permet de détecter et d'isoler les (multi-)faute(s).

Comme nous le verrons, nous avons automatisé toutes les étapes de cette procédure.

### 3.3.1 Préliminaires

Supposons que les fautes $f_i$ appartiennent à des ensembles connexes de $\mathbb{R}$ et notons $\mathcal{F}$ le produit cartésien de ces ensembles.

On suppose par la suite que les contraintes sur $p \in \mathcal{P}$ et $f \in \mathcal{F}$ et éventuellement les contraintes liant les fautes et les paramètres peuvent être formulées à l'aide d'équations et/ou d'inégalités algébriques. Cette hypothèse permet de travailler avec des ensembles semi-algébriques pour lesquels des outils de calculs algébriques ont été développés (voir [26, 73, 218] par exemple).

Rappelons que $\mathbb{R}[X_1, \dots, X_n]$ est l'ensemble des polynômes à coefficients réels et tels que $X_1, \dots, X_n$ sont les $n$ indéterminées. De plus, un ensemble de solutions réelles d'un ensemble fini d'équations polynomiales à plusieurs variables (de la forme $P = 0$) et/ou d'inégalités polynomiales (de la forme $Q \neq 0$) de $\mathbb{R}[X_1, \dots, X_n]$ est appelée ensemble semi-algébrique.

Soit $C_{p,f}$ l'ensemble de toutes les équations et inégalités vérifiées par les composantes des vecteurs de paramètres et de fautes du modèle et $C_{p,f}$ l'ensemble semi-algébrique défini par $C_{p,f}$. Pour prendre en compte les conditions initiales, les relations algébriques induites par ces conditions peuvent être rajouter à l'ensemble $C_{p,f}$.

**Exemple 3.3.1.** Reprenons l'exemple de la masse (Exemple 3.6) :

$$ \begin{cases} \dot{x}_1 = x_2, \\ \dot{x}_2 = -k(f_1 - 1)^2 x_1 + (d + f_2)u. \end{cases} \tag{3.18} $$

Supposons que les paramètres $k, d$ et les fautes, $f_1$ et $f_2$ vérifient les contraintes suivantes :

$$ \left\{ \begin{array}{l} 0 < k < 4 \\ 1 \le d \\ 0 \le f_1 < 2, \\ 0 \le f_2 < 2. \end{array} \right. \tag{3.19} $$

Ces contraintes algébriques peuvent être vues comme un ensemble d'équations et d'inégalités polynomiales dont les indéterminées sont $k, d, f_1$ et $f_2$. Ces inégalités algébriques définissent $C_{p,f}$ et $C_{p,f}$ est l'ensemble des valeurs admissibles des paramètres et des fautes correspondant.

Pour définir le vecteur de multi-fautes, nous introduisons les notations suivantes. Soient $\mathcal{N}$ un sous-ensemble de $\{1, \dots, e\}$ et $f_{\mathcal{N}}$ le vecteur de multi-fautes dont les composantes $f_i$ sont non nulles si $i \in \mathcal{N}$ et égales à 0 autrement. Naturellement, $f_{\mathcal{N}}$ appartient à $\mathcal{F}_{\mathcal{N}} = \{f \in \mathcal{F} | f_i \neq 0 \text{ si } i \in \mathcal{N} \text{ et } f_i = 0 \text{ si } i \notin \mathcal{N}\}$ et $\mathcal{F}_{\mathcal{N}}$ est un ensemble semi-algébrique par construction. Quand seule une composante de $f$ est non nulle, on retrouve le cas simple faute.

Comme la diagnosticabilité et la détectabilité du modèle (3.1) peut dépendre de $p = (p_1, \dots, p_m)$, nous considérons par la suite l'ensemble $R = \mathbb{R}[p_1, \dots, p_m]$ des polynômes en les indéterminées $p_i$ à coefficients réels.

### 3.3.2 Construction de la signature algébrique à partir du résumé exhaustif

La définition suivante introduit la notion de signature algébrique et permet de caractériser les multi-fautes. Elle diffère de la définition de diagnosticabilité fonctionnelle car elle est composée d'expressions algébriques. Par la suite, nous l'avons définie à partir du résumé exhaustif, la suite formée par les coefficients des polynômes ES.

**Définition 22.** Soit $ASig = (ASig_1, \dots, ASig_l)$ un vecteur d'expressions algébriques admettant $f_1, \dots, f_e$ comme indéterminées et à coefficients dans $R$. Une signature algébrique est une fonction $ASig$ définie par :

$$ \begin{array}{rcl} ASig: & \mathcal{F} & \to R^l \\ & f & \mapsto (ASig_1(f), \dots, ASig_l(f)). \end{array} $$

**Exemple 3.3.2.** En reprenant l'exemple 3.3.1, la fonction $\phi$ définie par
$$ \phi(f) = (k(f_1 - 1)^2, -d - f_2) $$

$$ \text{donne la signature algébrique } ASig(f) = (k(f_1 - 1)^2, -d - f_2). $$

La comparaison des images de ASig pour deux multi-fautes différentes donnera un moyen de les discriminer.

Si le modèle est contrôlé, nous proposons deux définitions : la forte diagnosticabilité algébrique et la faible diagnosticabilité algébrique, la première étant vraie pour toute entrée, la seconde pour un entrée.

**Définition 23.** Soit $\mathcal{N}$ et $\mathcal{N}'$ deux sous-ensembles distincts de $\{1, ..., e\}$. Les multi-fautes de $\mathcal{F}_\mathcal{N}$ et de $\mathcal{F}_{\mathcal{N}'}$ sont dites fortement algébriquement discriminables (resp. faiblement algébriquement discriminables) si il existe une signature algébrique telle que pour toute entrée u (resp. une entrée),

$$ ASig(\mathcal{F}_\mathcal{N}) \cap ASig(\mathcal{F}_{\mathcal{N}'}) = \emptyset. \qquad (3.20) $$

Cette égalité est satisfaite en particulier quand il existe un indice i tel que $ASig_i(\mathcal{F}_\mathcal{N}) \cap ASig_i(\mathcal{F}_{\mathcal{N}'}) = \emptyset$.

Si, pour deux sous ensembles distincts quelconques $\mathcal{N}$ et $\mathcal{N}'$ de $\{1, ..., e\}$, les multi-fautes de $\mathcal{F}_\mathcal{N}$ et $\mathcal{F}_{\mathcal{N}'}$ sont fortement algébriquement discriminables (resp. faiblement algébriquement discriminables), le modèle est dit fortement algébriquement diagnosticable (resp. faiblement algébriquement diagnosticable).

Dans le cas de modèles non contrôlés, la définition consiste à omettre la notion d'entrée dans les définitions précédentes.

La notion de détectabilité de fautes consiste à détecter la présence de fautes ou multi-fautes agissant sur le système. Nous avons proposé la définition suivante dans le cas de la signature algébrique.

**Définition 24.** Un ensemble de vecteurs multi-fautes $\mathcal{F}_\mathcal{N}$ est **algébriquement détectable si**

$$ ASig(\mathcal{F}_\mathcal{N}) \cap ASig(\mathcal{F}_\emptyset) = \emptyset, $$

$ASig(\mathcal{F}_\emptyset)$ étant la signature algébrique évaluée lorsque aucune faute agit sur le système.

La détectabilité consiste donc à comparer les valeurs de la signature algébrique avec celle obtenue quand l'ensemble $\mathcal{N}$ est vide, c'est à dire quand aucune faute agit sur le système.

**Exemple 3.3.3.** Reprenons l'exemple 3.1.1 (et 3.3.1). Les signatures possibles pour les multi-fautes sont $ASig(f_\emptyset) = (k, -d)$, $ASig(f_{\{1\}}) = (k(f_1 - 1)^2, -d)$, $ASig(f_{\{2\}}) = (k, -d - f_2)$ et $ASig(f_{\{1,2\}}) = (k(f_1 - 1)^2, -d - f_2)$. Les contraintes sur les paramètres et les fautes impliquent que les images des signatures algébriques ne s'intersectent pas. Par exemple, comme pour toute faute $f_2 \in (0,2)$, $-d \neq -d - f_2$, on a $ASig(f_{\{1\}}) \cap ASig(f_{\{1,2\}}) = \emptyset$. Ainsi, les multi-fautes peuvent être discriminées et le modèle est algébriquement diagnosticable.

Dans cet exemple, la fonction $\phi$ n'est pas injective : les valeurs $1/2$ et $3/2$ de $f_1$ donneront la même valeur de $\phi(f)$. Rappelons que le modèle est identifiable si et seulement si $\phi$ est injective [62]. Donc l'identifiabilité des fautes implique la diagnosticabilité algébrique car les images par $\phi$ des fautes auront une valeur distincte. L'exemple du ressort montre que la diagnosticabilité algébrique n'implique pas l'identifiabilité. En effet, si $\phi$ n'est pas injective, la discrimination des multi-fautes est possible.

La définition de la signature algébrique à partir du résumé exaustif nous a paru insuffisante dans le sens que deux fautes distinctes agissant sur les mêmes composantes ne peuvent être discriminées. De plus, elles ne sont pas exploitables numériquement car elles dépendent des fautes non connues a priori. Une approche naturelle consiste à obtenir une expression explicite des composantes des fautes en fonction des paramètres du modèle et des composantes $\phi_1, ..., \phi_N$ de $\phi(f)$. Cette approche revient à un problème d'inversion d'un système algébrique et n'aboutit pas en général. Nous avons donc développé une méthode pour obtenir des expressions algébriques ne dépendant pas des fautes et qui caractérisent leur présence. De telles expressions peuvent être calculées à l'aide d'une procédure automatique basée sur les bases de Groebner ([51, 76]) et utilisées pour définir une signature algébrique.

Pour alléger la présentation, nous supposons que $(\gamma_k^i)_{1 \le k \le n_i}$ sont des polynômes dans $R[f_1, ..., f_e]$. Lorsque ces expressions sont des fractions rationnelles, elles peuvent s'exprimer comme polynômes en introduisant de nouvelles variables égales aux inverses des dénominateurs. Des conditions sur la non nullité des dénominateurs peuvent être ajoutées à $C_{p,f}$.

Étant donné une multi-faute $f \in \mathcal{F}_\mathcal{N}$ ($\mathcal{N} \subset \{1, ..., e\}$), soit $E_\mathcal{N}$ un ensemble de polynômes

$$ E_\mathcal{N} = \{\gamma_1^1(p, f) - \phi_1, ..., \gamma_s^{n_s}(p, f) - \phi_N\} \cup \{v_i f_i - 1 | i \in \mathcal{N}\} \cup \{f_i | i \notin \mathcal{N}\} $$

où $v_i$ sont de nouvelles indéterminées.

Dans la définition de $E_\mathcal{N}$, les ensembles $\{v_i f_i - 1 | i \in \mathcal{N}\}$ et $\{f_i | i \notin \mathcal{N}\}$ caractérisent les multi-fautes de $\mathcal{F}_\mathcal{N}$. Le
premier signifie que la multi-faute $f \in \mathcal{F}_N$ a ses composantes $f_i$, $i \in \mathcal{N}$ non nulles et le deuxième que ses autres composantes $i \notin \mathcal{N}$ sont nulles.

Considérons l'idéal polynômial $I_N$ généré par $E_N$, c'est à dire l'ensemble de toutes les combinaisons linéaires des éléments de $E_N$ dans $R[v_1, \dots, v_e, f_1, \dots, f_e, \phi_1, \dots, \phi_N]$.

Une base de Groebner de l'idéal $I_N$ est calculée avec l'ordre d'élimination consistant à éliminer tout d'abord les indéterminées $v_i$ et $f_i$. L'intersection $G_N$ de cette base de Groebner et de $R[\phi_1, \dots, \phi_N]$ génère l'idéal $J_N = I_N \cap R[\phi_1, \dots, \phi_N]$ (Voir [51]). Clairement, tout polynôme de $G_N$ s'annule lorsqu'une multi-faute $f \in \mathcal{F}_N$ agit sur le système.

Pour toutes les multi-fautes possibles $f_N$, les ensembles $G_N$ sont calculés. Les polynômes de $\cup_{N \subset \{1, \dots, m\}} G_N$ s'annulant pour toutes les multi-fautes, i.e. les polynômes de $\cap_{N \subset \{1, \dots, m\}} I_N$, sont retirés de cet ensemble. Les polynômes restant sont utilisés pour définir les composantes de la signature algébrique.

L'algorithme retournant la signature algébrique est résumé ci-dessous.

Signature\_algébrique

1. Pour chaque sous ensemble $\mathcal{N}$ de $\{1, \dots, e\}$, nous considérons une multiple faute générique $f_N$ et nous appliquons les étapes suivantes à cette multiple faute.

(a) Calcul de la base de Groebner de l'idéal $I_N$ engendré par $E_N$ avec l'ordre lexicographique $v_{i_1} \succ \dots \succ v_{i_l} \succ f_1 \succ \dots \succ f_m \succ \phi_1 \succ \dots \succ \phi_N \succ p_1 \succ \dots \succ p_m$.

(b) Détermination de l'intersection, $G_N$, de cette base de Groebner et de $R[\phi_1, \dots, \phi_N]$.

2. Élimination dans $\cup_{N \subset \{1, \dots, m\}} G_N$ des polynômes s'annulant pour toutes les multi-fautes, i.e. les polynômes de l'idéal $\cap_{N \subset \{1, \dots, m\}} I_N$.

3. Classement de façon arbitraire de tous les polynômes restant pour former la suite $ASig = (ASig_1, \dots, ASig_l)$.

4. Renvoi de $ASig$.

Cet algorithme permet d'obtenir une signature algébrique de toutes les multi-fautes telle que :

$$ 
\begin{array}{rcl}
ASig: & \mathbb{R}^e & \to (R[\phi_1, \dots, \phi_N])^l \\
& f & \mapsto (ASig_1(\phi), \dots, ASig_l(\phi)).
\end{array}
\qquad (3.21) $$

**Exemple 3.3.4.** Considérons le résumé exhaustif $\phi(f_1, f_2) = (k(f_1 - 1)^2, -d - f_2)$ donné dans l'exemple 3.3.2. L'algorithme Signature\_algébrique retourne la signature $ASig$ définie par $ASig(f_1, f_2) = (\phi_1 - k, \phi_2 + d)$ dont chaque composante s'annule pour au moins une multi-faute.

Par construction, la signature $ASig(f)$ ne dépend pas explicitement de $f$. Cependant, la présence de (multi-)fautes apparaît dans les valeurs numériques de $(\phi_1, \dots, \phi_N)$ et donc dans $ASig$.

En comparant une estimation de $ASig(f)$ et les composantes attendues nulles de la liste $ASig(f_N)$, des (multi-)fautes peuvent être écartées. Toutefois, une telle comparaison peut ne pas être suffisante pour discriminer certaines multi-fautes. En effet, des polynômes de $G_N$ apparaissant dans $ASig(f_N)$ s'annuleront quand la faute $f_N$ a lieu mais les autres composantes de la signature $ASig(f_N)$ peuvent s'annuler également pour certaines valeurs particulières de paramètres et de fautes. C'est pourquoi, nous avons ajouté des critères supplémentaires pour améliorer la discrimination des multi-fautes.

### 3.3.3 Critères pour différentier les signatures

Pour élaborer ces critères, nous avons utilisé des outils semi-algébriques, en particulier, ceux permettant d'étudier les solutions réelles d'équations et d'inéquations polynomiales [10]. Cette approche permet de prendre en compte l'ensemble $C_{p,f}$ qui est l'ensemble des contraintes sur les paramètres et les fautes. Comme nous le verrons dans l'exemple 3.3.5, celles-ci peuvent jouer un rôle fondamental pour discriminer les signatures de (multi-)fautes.

Les trois résultats suivants sont basés sur le test d'ensemble vide d'ensembles semi-algébriques. Ce test peut-être fait grâce à des outils de calcul algébrique ([73, 218]).

Le premier critère (resp. le second) consiste à déterminer si la $k$-ème composante de $ASig(f_N)$ s'annule pour au moins une valeur réelle de la multi-faute $f \in \mathcal{F}_N$ (resp. ne s'annule jamais).

Pour tout $\mathcal{N} \subset \{1, \dots, m\}$, considérons l'ensemble $S_N$ des équations et inéquations polynomiales définies par

$$ S_N := \{\gamma_1^1(p, f) = \phi_1, \dots, \gamma_s^{n_s}(p, f) = \phi_N\} \cup C_{p,f} \cup \{v_i f_i = 1 \mid i \in \mathcal{N}\} \cup \{f_i = 0 \mid i \notin \mathcal{N}\} $$

où $v_i$ sont de nouvelles indéterminées.
**Critère 1.** Si l'ensemble semi-algébrique défini par $S_N \cup \{ASig_k(f_N) = 0\}$ est vide alors la kème composante de $ASig(f_N)$ ne s'annule jamais.

**Critère 2.** Si l'ensemble semi-algébrique défini par $S_N \cup \{v_k ASig_k(f_N) - 1 = 0\}$ est vide alors la kème composante de $ASig(f_N)$ est égale à 0.

Pour certains systèmes particuliers, une composante de la signature qui s'annule caractérise les multi-fautes dont la ième composante est nulle.

**Critère 3.** Soit $S$ l'ensemble semi-algébrique défini par $S := \{\gamma_1^1(p, f) = \phi_1, \dots, \gamma_s^{n_s}(p, f) = \phi_N\} \cup C_{p,f}$. Si l'ensemble des solutions réelles $S \cup \{ASig_j(f) = 0, v_i f_i - 1 = 0\}$ et $S \cup \{v_j ASig_j(f) - 1 = 0, f_i = 0\}$ est vide alors $ASig_j(f) = 0$ est équivalent à $f_i = 0$.

Si le critère 3 est satisfait pour toutes les composantes $f_i$ de $f = (f_1, \dots, f_m)$, il n'est clairement pas utile d'appliquer les critères 1 et 2 sur toutes les $m!$ possibles multi-fautes.

En utilisant les trois critères ci-dessus, les valeurs attendues de $ASig(f)$ quand une (multi-)faute $f$ agit peuvent être regroupées dans un tableau. Par la suite, nous prendrons les conventions suivantes : pour toute multi-faute $f$,

- une cellule contenant $\emptyset$ signifie que les critères 1 ou 3 assurent que la composante de $ASig(f)$ ne s'annulera jamais quand la multi-faute agit sur le système;

- un 0 dans une cellule signifie que les critères 2 ou 3 assurent que $ASig_i(f)$ est nécessairement égal à 0 quand la multi-faute agit sur le système;

- une cellule vide indique que la composante de la signature s'annule pour certaines valeurs de $(p, f)$ et ne s'annule pas pour d'autres valeurs de $(p, f)$.

**Exemple 3.3.5.** Continuons l'exemple 3.3.1 et reprenons la signature algébrique $ASig(f_1, f_2) = (\phi_1 - k, \phi_2 + d)$ retournée par notre algorithme dans l'exemple 3.3.4.

Si l'ensemble des contraintes $C_{p,f} = \{0 < k < 4, 1 \le d, 0 \le f_1 < 2, 0 \le f_2 < 2\}$ est pris en compte, les deux premiers critères fournissent les caractéristiques de la signature algébrique des multi-fautes possibles. Elles sont résumées dans le tableau suivant :

<table><thead><tr><th>f</th><th>ASig<sub>1</sub>(f)</th><th>ASig<sub>2</sub>(f)</th></tr></thead><tbody><tr><td>f<sub>{}</sub></td><td>0</td><td>0</td></tr><tr><td>f<sub>{1}</sub></td><td>Ø</td><td>0</td></tr><tr><td>f<sub>{2}</sub></td><td>0</td><td>Ø</td></tr><tr><td>f<sub>{1,2}</sub></td><td>Ø</td><td>Ø</td></tr></tbody></table>

Clairement, les valeurs de $ASig(f)$ et, plus précisément, les valeurs de $ASig_1(f)$ et de $ASig_2(f)$ sont suffisamment discriminantes pour toutes les multi-fautes possibles. Ce résultat peut être également obtenu en appliquant le critère 3 aux deux composantes de la signature : ce critère permet de montrer l'équivalence entre $f_1 = 0$ (resp. $f_2 = 0$) et $ASig_1 = 0$ (resp. $ASig_2 = 0$).

Si on ne considère pas les contraintes sur les paramètres et les fautes, le tableau suivant est obtenu :

<table><thead><tr><th>f</th><th>ASig<sub>1</sub>(f)</th><th>ASig<sub>2</sub>(f)</th></tr></thead><tbody><tr><td>f<sub>{}</sub></td><td>0</td><td>0</td></tr><tr><td>f<sub>{1}</sub></td><td></td><td>0</td></tr><tr><td>f<sub>{2}</sub></td><td>0</td><td>Ø</td></tr><tr><td>f<sub>{1,2}</sub></td><td></td><td>Ø</td></tr></tbody></table>

Ce tableau montre que les contraintes jouent un rôle important dans l'étude de la diagnosticabilité. Plus précisément, les outils semi-algébriques assurent que, pour des valeurs particulières de paramètres, la faute $f_{\{1\}}$ ne peut pas être détectée d'après les signatures algébriques de $f_{\{\}}$ et $f_{\{1\}}$. On peut faire la même remarque pour la discrimination des multi-fautes $f_{\{2\}}$ et $f_{\{1,2\}}$.

**Remarque 23.** D'autres critères algébriques utilisant les calculs de bases de Groebner peuvent être développés pour obtenir de l'information sur les valeurs possibles de $ASig(f_N)$ quand une multi-faute $f_N$ agit. Par exemple,

1. si la base de Groebner $E_N \cup \{ASig_k(f_N)\}$ est égale à $\{1\}$ alors la kème composante de $ASig(f_N)$ ne s'annule jamais. En effet, dans ce cas, les fonctions polynomiales associées aux polynômes de $E_N \cup \{ASig_k(f_N)\}$ n'ont pas de zéro complexe en commun ($[51]$) et par conséquent pas de solution réelle;
2. par construction de $I_N$, nous avons que si $ASig_k(f_N)$ appartient à $I_N$ alors $ASig_k(f_N) = 0$.

Même si les contraintes sur les paramètres impliquant les inégalités ne peuvent pas être prises en compte, ces critères peuvent être testés plus rapidement dans la pratique que les critères 1 ou 2.

### 3.3.4 Application

Dans cette section, nous allons reprendre l'exemple du réservoir 3.2.3 utilisé dans le cadre de la signature fonctionnelle. Nous allons montrer comment à partir de la signature algébrique et des critères que nous avons définis dans la section précédente on peut traiter la diagnosticabilité et le diagnostic d'un système dans le cas (multi-)fautes. Nous avons implémenté les algorithmes pour la construction de la signature algébrique et les critères 1, 2 et 3 en Maple 18. Cette première étape nous a permis de faire une étude de diagnosticabilité du système en nous fournissant la table donnant les valeurs attendues des signatures en fonction des multi-fautes possibles. Cette table a ensuite constitué l'entrée d'un programme Scilab permettant d'estimer numériquement les signatures algébriques des multi-fautes à partir des sorties simulées bruitées. La comparaison entre les valeurs numériques et les valeurs attendues de la signature ont permis de détecter et d'isoler les multi-fautes.

Reprenons l'exemple 3.2.3 du réservoir mais avec une seule observation².

$$ egin{cases} \dot{x}_1(t,p) = p_1 u(t) - p_2 \sqrt{x_1(t,p)}, & x_1(0) = 0.3, \\ \dot{x}_2(t,p) = p_3 \sqrt{x_1(t,p)} - p_4 \sqrt{x_2(t,p)}, & x_2(0) = 0.6, \\ y(t,p) = p_5 \sqrt{x_1(t,p)}, & \end{cases} \qquad (3.22) $$

où $p_1, \dots, p_5$ sont des paramètres positifs.

Les fautes sont $f_1$ une faute additive sur le signal d'entrée, $f_2$ une faute additive agissant sur le capteur à la sortie du premier réservoir et $f_3 \in ]0; 1[$ un défaut de colmatage.

Utilisant les changements de variables $z_1(t,p) = \sqrt{x_1(t,p)}$, $z_2(t,p) = \sqrt{x_2(t,p)}$ et l'algorithme de Rosenfeld Groebner, l'ARR suivante a été obtenue :

$$ 2y\dot{y} - p_5(f_3 - 1)^2(p_1p_5f_1 + p_2f_2) - p_1p_5^2(f_3 - 1)^2u + p_2p_5(f_3 - 1)^2y - 2f_2\dot{y} = 0 $$

d'où on déduit le résumé exhaustif :

$$ \phi(f_1, f_2, f_3) = (-p_5(f_3 - 1)^2(p_1p_5f_1 + p_2f_2), -p_1p_5^2(f_3 - 1)^2, p_2p_5(f_3 - 1)^2, -2f_2). $$

L'application de l'algorithme *Signature\_Algebrique* nous a donné la signature algébrique :

$$ ASig(f) = (\phi_1, \phi_4, p_1 p_5^2 + \phi_2, -p_2 p_5 + \phi_3, -\phi_3 \phi_4 + 2\phi_1, -p_2 p_5 \phi_4 + 2\phi_1). $$

L'étape suivante consiste à calculer les valeurs attendues de $ASig(f)$. L'ensemble des contraintes considérées est :

$$ C_{p,f} = \{0 < p_1, \dots, 0 < p_5, 0 \le f_3 < 1\}. $$

qui correspondent à la signification physique des paramètres et des fautes du modèle. Les critères 1 et 2 permettent de déduire le tableau 3.1 :

Ce tableau montre que les composantes $ASig_2(f)$, $ASig_4(f)$ et $ASig_5(f)$ permettent de discriminer toutes les multi-fautes pour n'importe quelle entrée $u$. En effet, $ASig_2(f)$, $ASig_4(f)$ et $ASig_5(f)$ ne dépendent pas de la composante $\phi_2$ qui est le seul coefficient contenant l'entrée $u$. Le modèle est donc fortement algébriquement diagnosticable.

**Remarque 24.** 1. Ce résultat est valable même si la valeur des paramètres $p_1, p_3$ et $p_4$ n'est pas connue. Ainsi, la connaissance de toutes les valeurs des paramètres internes n'est pas nécessaire pour détecter et discriminer toutes les multi-fautes possibles.

2. La conclusion que le modèle est fortement algébriquement diagnosticable peut être obtenue en appliquant le critère 3 à $ASig_2(f)$, $ASig_4(f)$ et $ASig_5(f)$.

Dans les simulations, un simple contrôleur est utilisé pour contrôler le niveau de l'eau dans le réservoir supérieur. On a supposé qu'il suivait un signal carré de référence. Les paramètres du modèle sont $p_1 = p_2 = p_3 = p_4 = 0.3$, $p_5 = 1$. La sortie simulée est perturbée par un bruit gaussien $\eta$ tel que $\eta(t) \in [-0.01; 0.01]$. Ainsi $y(t) = \bar{y}(t) + \eta(t)$ où $\bar{y}$ est la valeur exacte de la sortie correspondant à la valeur exacte des paramètres. Les observations sont supposées être faites aux temps discrets $(t_i)_{i=1,...,M}$ sur l'intervalle $[0,50]$, le pas de temps étant $h = 0.5$.

Dans le scénario avec faute, nous supposons que les fautes sont introduites au temps $t = 22s$. Les figures 3.12 et 3.13 représentent différents cas de simples et multi-fautes agissant sur le système.
<table><caption>TABLE 3.1 – Valeurs numériques attendues des signatures algébriques</caption><thead><tr><th></th><th>ASig1(f)</th><th>ASig2(f)</th><th>ASig3(f)</th><th>ASig4(f)</th><th>ASig5(f)</th><th>ASig6(f)</th></tr></thead><tbody><tr><th>f{}</th><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>f{1}</th><td>Ø</td><td>0</td><td>0</td><td>0</td><td>Ø</td><td>Ø</td></tr><tr><th>f{2}</th><td>Ø</td><td>Ø</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>f{3}</th><td>0</td><td>0</td><td>Ø</td><td>Ø</td><td>0</td><td>0</td></tr><tr><th>f{1,2}</th><td></td><td>Ø</td><td>0</td><td>0</td><td>Ø</td><td>Ø</td></tr><tr><th>f{1,3}</th><td>Ø</td><td>0</td><td>Ø</td><td>Ø</td><td>Ø</td><td>Ø</td></tr><tr><th>f{2,3}</th><td>Ø</td><td>Ø</td><td>Ø</td><td>Ø</td><td>0</td><td>Ø</td></tr><tr><th>f{1,2,3}</th><td></td><td>Ø</td><td>Ø</td><td>Ø</td><td>Ø</td><td></td></tr></tbody></table>

<figure><img src="image_35.png" /><figcaption>FIGURE 3.12 – Niveau d'eau dans le réservoir supérieur, y, durant la simulation sans faute et durant les simulations avec de simples fautes introduites à t = 22s.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>FIGURE 3.13 – Simulations du niveau d'eau dans le réservoir supérieur (valeurs de y) pour le cas sans faute (courbe noire) et les cas multi-fautes introduites à t = 22s (courbes bleues, vertes et rouges).</figcaption></figure>

Les dérivées sont estimées en utilisant une méthode basée sur les B-splines [97]. Pour estimer $\phi$, nous nous sommes ramenés à la résolution d'un système linéaire sur-dimensionné. Pour cela, nous avons réécrit l'ARR à chaque instant discret $t_i$, d'où un système linéaire de dimension $M$ dont les inconnues sont les composantes de $\phi$. Notons $y_p(t_i)$ l'estimée de $\dot{y}(t_i)$. Le système a la forme suivante :

— Dans le cas avec faute,

$$A^f X_f = b \quad (3.23)$$

$$\text{avec } X_f = (-p_5(f_3-1)^2(p_1p_5f_1+p_2f_2), -p_1p_5^2(f_3-1)^2, p_2p_5(f_3-1)^2, -2f_2)(=\phi(f)), \\ A_i^f = (1, u(t_i), y(t_i), y_p(t_i)) \text{ et } b_i = -2y(t_i)y_p(t_i).$$

— Dans le cas sans faute,

$$AX_0 = b \quad (3.24)$$

$$\text{avec } X_0 = (p_1p_5^2; p_2p_5), A_i = (u(t_i), y(t_i)) \text{ et } b_i = -2y(t_i)y_p(t_i).$$

Ces systèmes sont résolus en utilisant la factorisation QR.

Le système (3.24) est utilisé pour détecter le temps $t_d$ à partir duquel la multi-faute agit. À partir des 10 premiers pas de temps, la matrice $A$ et le vecteur $b$ sont construits. Ensuite, à chaque pas de temps, ils sont complétés et le système (3.24) est résolu avec cette nouvelle matrice $A$ et ce nouveau vecteur $b$. L'estimée de $X_0$ est comparée avec la valeur nominale. Si leur différence en norme 2 est plus grande que $\varepsilon = 10^{-3}$, une multi-faute est détectée.

Une fois une multi-faute détectée, le système (3.23) sert à discriminer la (multiple) faute. Au moins quatre valeurs de temps après la détection de la multiple faute en $t_d$ sont nécessaires puisque $X_f$ est un vecteur de dimension 4. Remarquons que les multi-fautes peuvent être détectées et discriminées tous les 0.5 correspondant aux mesures de sortie. La signature $ASig$ est alors estimée et le tableau 2.1 est utilisé pour discriminer la multi-faute qui agit. Les résultats sont résumés dans le tableau 3.2.

Même si nous avons fait le simple choix de définir un seuil constant, toutes les multi-fautes ont été détectées et discriminées. Ces résultats pourraient sans doute être améliorés en utilisant une autre méthode pour définir un seuil en fonction du bruit.

### 3.3.5 La détection de faute dans le cadre ensembliste

Cette partie reprend une partie de l'article [102] et fait le lien entre les définitions d'identifiabilité ensemblistes donnés à la section 2.5 et la détection.

Dans le cadre du projet ANR MAGIC-SPS, nous nous sommes intéressées aux méthodes de détection et d'identification des défauts pour les modèles non linéaires à incertitudes bornées. La détection de faute dans le cadre ensembliste peut se faire à partir de l'estimation des états du système ou des paramètres. Lorsque le modèle est

<footnote n="2">2. Ce choix est fait pour que le modèle ne soit pas identifiable presque partout.</footnote>
<table><caption>TABLE 3.2 – Temps de détection et de discrimination. Un temps de détection égal à 0 signifie que la multi-faute est détectée à 22s dès que la multi-faute agit. Un temps de discrimination de 1.5s signifie qu'il aura fallu 4 pas de temps (22; 22.5; 23; 23.5) pour discriminer la faute.</caption><thead><tr><th>(Multi-)fautes f</th><th>Temps de Détection</th><th>Temps de Discrimination</th></tr></thead><tbody><tr><td>f<sub>{1}</sub> = (0.5, 0, 0)</td><td>0</td><td>3</td></tr><tr><td>f<sub>{2}</sub> = (0, 0.5, 0)</td><td>0.5</td><td>1.5</td></tr><tr><td>f<sub>{3}</sub> = (0, 0, 0.5)</td><td>0.5</td><td>2</td></tr><tr><td>f<sub>{1,3}</sub> = (0.5, 0, 0.1)</td><td>0</td><td>1.5</td></tr><tr><td>f<sub>{1,3}</sub> = (0.5, 0, 0.7)</td><td>0.5</td><td>11*</td></tr><tr><td>f<sub>{1,2}</sub> = (0.5, 0.5, 0)</td><td>0</td><td>1.5</td></tr><tr><td>f<sub>{2,3}</sub> = (0, 0.5, 0.1)</td><td>0</td><td>1.5</td></tr><tr><td>f<sub>{2,3}</sub> = (0, 0.5, 0.7)</td><td>0</td><td>1.5</td></tr></tbody></table>

non linéaire, les ensembles à caractériser peuvent être non convexes et même se composer de plusieurs ensembles disjoints. Dans ce dernier cas, l'analyse par intervalles donne l'enveloppe connexe de l'union de ces ensembles entraînant le non déclenchement d'alarme. La difficulté est alors d'obtenir des vecteurs de solution intervalles plus conservateurs. Des méthodes récentes telles que [2] ou celle que nous avons proposées dans [102] fournissent des améliorations dans ce sens.

Dans [102], nous avons repris le travail d'identifiabilité et d'estimation de paramètres dans le cadre ensembliste et présenté aux sections 2.5.4 et 2.5.7. Ce travail tire parti de la méthode basée sur l'algèbre différentielle pour vérifier l'identifiabilité ensembliste et son homologue opérationnel, la μ-identifiabilité ensembliste. Ces notations fournissent un moyen d'étudier différents aspects de l'identifiabilité pour les systèmes à erreurs bornées, en particulier les systèmes qui représentent une famille infinie de systèmes non linéaires. En construisant le schéma d'estimation des paramètres sur l'analyse de l'identifiabilité, nous garantissons que l'ensemble de solutions se réduit à un ensemble connexe, évitant ainsi le pessimisme des méthodes ensemblistes. L'identifiabilité est étroitement liée à la diagnosticabilité car elle fournit la garantie que deux situations correspondant à des valeurs de paramètres différents sont distinguables. La méthode proposée a été appliquée à l'exemple décrivant la capacité d'un récepteur du macrophage à recevoir des macromolécule solubles spécifiques et présenté à la Section 2.1. Dans [103], le modèle a été donné dans un cadre stochastique. L'incertitude a été prise en compte au moyen d'hypothèses appropriées sur les distributions de probabilité de l'erreur du bruit et du modèle. Dans [102], la sortie y a été perturbée par un bruit additif borné $\eta$, $\eta(t) \in [\eta(t)]$ et nous nous sommes ramenées à la résolution d'un système $0 \in [A][x] - [b]$ avec l'algorithme SIVIA (voir Annexe 6).

Différents scénarios sans faute et avec faute ont été considérés. Le scénario sans faute a donné les intervalles nominaux de référence pour chaque paramètre. Pour chaque scénario, les paramètres ont été estimés correctement avec une précision raisonnable et les intervalles obtenus pour les paramètres avec fautes avaient une intersection disjointe avec les intervalles nominaux.

## 3.4 Conclusion

Dans le cadre de la méthode basée sur les ARRs, nous avons proposé deux définitions de diagnosticabilité : la diagnosticabilité fonctionnelle et la diagnosticabilité algébrique pour les modèles dynamiques non linéaires. Contrairement aux définitions classiques, l'étude de ces définitions met en évidence certaines des propriétés analytiques des relations de redondance liées à la faute agissant sur le système. Dans le premier cas, la propriété d'identifiabilité permet d'assurer qu'une ou des fautes pourront être discriminées en comparant et évaluant la faute en fonction de l'ensemble des trajectoires du résidu. Dans le second cas, une table obtenue à partir d'outils semi-algébriques certifie que des (multi)-fautes peuvent ou non être discriminées en évaluant des relations algébriques. Les méthodes développées s'appuient sur plusieurs outils d'algèbre (différentielle et algébrique) permettant de garantir les résultats théoriques. De celles-ci, nous en avons déduit des méthodes numériques pour détecter, discriminer et isoler la ou les fautes.

<footnote>* $f_3 \neq 0$ est tout d'abord détecté à $t = 22.5s$ et la multi-faute $f_{\{1,3\}}$ est discriminée à $t = 33s$.</footnote>
<figure><img src="image_37.png" /></figure>
# Chapitre 4

## Réseaux

### 4.1 Introduction

Je me suis également intéressée aux réseaux d'oscillateurs [30, 31, 50, 154]. Ceux-ci sont très utilisés pour décrire des processus dynamiques qui se produisent dans de nombreux domaines de la science comme la biologie, l'étude des populations. La caractéristique principale de ces réseaux d'interaction complexes est l'émergence de nouvelles propriétés comme la synchronisation entre les différents éléments.

J'ai tout d'abord étudié les conditions de synchronisation de réseaux dont les noeuds sont constitués de systèmes dynamiques non linéaires, non nécessairement identiques et dont les termes de synchronisation peuvent être non linéaires [154]. La résolution de ce problème nécessite la prise en compte des dynamiques individuelles des oscillateurs ainsi que la topologie du graphe, ce qui rend ce problème difficile à résoudre dans certains cas. La littérature sur ce thème est abondante et de nombreux résultats portant sur des conditions suffisantes pour obtenir la synchronisation complète de systèmes sont donnés (voir [159]). Classiquement, deux cycles limites couplés sont synchronisés dès que leur évolution en temps est périodique avec la même période et parfois avec la même phase [175]. Une question importante sur la synchronisation d'un réseau d'oscillateurs est de déterminer la stabilité de l'état de synchronisation. Il faut alors considérer les propriétés du réseau et les vecteurs d'état des oscillateurs (voir par exemple [16, 214, 215, 217, 226]). Deux méthodes ont été proposées dans la littérature. La première ("master stability function" en anglais), est basée sur le calcul de l'exposant de Lyapunov et les valeurs propres de la matrice de connectivité [156]. Cependant, le calcul des valeurs propres peut être difficile pour de grands réseaux d'oscillateurs. Une seconde méthode est la méthode de stabilité de connexion de graphe ("connection graph stability method" en anglais) [16]. Elle lie l'étude d'une fonction de Lyapunov à la topologie du graphe. Après avoir été développée sur des réseaux dont les termes de synchronisation sont linéaires, elle a ensuite été étendue aux graphes non orientés [13, 14, 15].

Nous avons étendu un des résultats existant dans le cas où les fonctions de couplage sont non linéaires. Les conditions suffisantes de synchronisation que nous avons proposées sont basées sur l'étude d'une fonction de Lyapunov et l'utilisation d'une pseudo-métrique. Ensuite, nous avons donné deux résultats pour obtenir la synchronisation complète. Le premier donne une condition sur la force de synchronisation. Le second est une version locale de la première, c'est à dire quand les oscillateurs sont proches de la variété de synchronisation. Nous avons également donné une condition d'existence des trajectoires des oscillateurs.

Ensuite, dans les articles [30, 31, 50], je me suis intéressée aux liens entre la synchronisation et les paramètres de bifurcations. Dans [30, 31], présentés dans le prochain chapitre, nous avons étudié la synchronisation de modèles décrivant des comportements humains en situation de catastrophe dans le cadre du contrôle. À un paramètre de bifurcation, nous avons associé un type de comportements. Selon sa valeur, la zone géographique représente une zone de sécurité ou une zone de danger.

Dans [50], nous avons étudié les connexions entre une bifurcation de Hopf d'un neurone modélisé par un système de Hindmarsh-Rose (HR) et une bifurcation d'un réseau de neurones acyclique direct formé par des systèmes HR distincts. En particulier, nous nous sommes demandées comment évoluaient les valeurs d'une bifurcation dans un réseau en fonction des valeurs des paramètres de couplage.
## 4.2 Synchronisation de systèmes non orientés

### 4.2.1 Extension de la méthode de stabilité de connexion de graphes

Soient $G$ un graphe connexe non orienté et $n$ son nombre d'arêtes, $\mathcal{E}$ l'ensemble de ses arêtes. L'arête de $G$, non orientée, du noeud $i$ au noeud $j$ sera notée $(i,j)$.

Les systèmes dynamiques considérés dans cette partie sont de la forme :

$$ egin{cases} 
\dot{X}_1 = F_1(X_1, t) - \epsilon \sum_{(1,j) \in \mathcal{E}} h(X_1, X_j), \\ 
\vdots \\ 
\dot{X}_n = F_n(X_n, t) - \epsilon \sum_{(n,j) \in \mathcal{E}} h(X_n, X_j), 
\end{cases} 
\quad (4.1) $$

où

— $X_i = (X_i^1, \dots, X_i^d)^T$ est le vecteur composé des $d$ coordonnées du $i$-ième oscillateur,

— $F_i = (F_i^1, \dots, F_i^d)^T$ un vecteur de fonctions,

— $h = (h^1, \dots, h^d)^T$ est la fonction de synchronisation qui définit le vecteur de couplage entre les oscillateurs,

— le paramètre réel $\epsilon$ est la force de synchronisation appliquée au réseau d'oscillateurs.

Le système (4.1) de conditions initiales $(X_1(0), \dots, X_2(0), \dots, X_n(0))^T$, est dit complètement synchronisé si, pour tout $(i,j) \in [1,n]$,

$$ \|X_i(t) - X_j(t)\| \xrightarrow[t \to +\infty]{} 0. $$

Cela signifie que le vecteur $(X_1, \dots, X_n)$ approche la variété de synchronisation définie par

$$ X_1(t) = X_2(t) = \dots = X_n(t). $$

En particulier, ceci implique que les oscillateurs ont le même comportement asymptotique (comme des trajets-toires chaotiques, des solutions stables et périodiques). La synchronisation complète du système (4.1) est dite *globale* si elle a lieu pour n'importe quel état initial des oscillateurs. Si le système (4.1) synchronise complètement dans un voisinage de la variété de synchronisation, la synchronisation est dite *locale*.

Par la suite, nous noterons les différences $\Delta_{i,j} = X_i^T - X_j^T$ et le vecteur constitué de l'ensemble des différences

$$ \Delta = (\Delta_{1,2}, \dots, \Delta_{1,n}, \Delta_{2,3}, \dots, \Delta_{2,n}, \dots, \Delta_{n-1,n})^T. $$

Ainsi, prouver la synchronisation complète du système (4.1) est équivalent à prouver que $\|\Delta(t)\| \xrightarrow[t \to +\infty]{} 0$.

Les résultats de synchronisation vont dépendre de la notion de pseudo-métrique.

Soit $D$ un ensemble non vide. Une pseudo-métrique est une application $\varphi : D \times D \to \mathbb{R}^+$ qui vérifie les axiomes suivants :

— $\varphi(z_1, z_1) = 0$,

— $\varphi(z_1, z_2) = \varphi(z_2, z_1)$ (propriété de symétrie),

— $\exists \rho > 0, \varphi(z_1, z_3) \le \rho(\varphi(z_1, z_2) + \varphi(z_2, z_3))$

On peut noter que toute métrique est une pseudo-métrique avec $\rho = 1$.

Soit $\varphi$ une pseudo-métrique sur l'ensemble $D$. Posons, pour tout $m \in \mathbb{N}^*$, $\rho(m)$ le plus petit réel tel que

$$ \varphi(z_1, z_{m+1}) \le \rho(m) [\varphi(z_1, z_2) + \dots + \varphi(z_m, z_{m+1})]. \quad (4.2) $$

Remarquons que $\rho(1) = 1$.

Dans les exemples suivants, les expressions de $\rho(m)$ sont les conséquences directes de la convexité des fonctions $x \to (x^2)^\alpha$ et $x \to x \tanh(x)$.

**Exemple 4.2.1.** 1. L'application $\varphi_\alpha : \mathbb{R}^2 \times \mathbb{R}^2 \to \mathbb{R}^+$ définie par

$$ \varphi_\alpha \left( \begin{pmatrix} x_1 \\ y_1 \end{pmatrix}, \begin{pmatrix} x_2 \\ y_2 \end{pmatrix} \right) = \left( (x_1 - x_2)^2 \right)^\alpha $$

avec $\alpha \ge 1/2$ est une pseudo-métrique pour laquelle $\rho(m) = m^{2\alpha-1}$.
2. Soit D la boule fermée de centre 0 et de rayon 1. L'application $\varphi: D \times D \rightarrow \mathbb{R}^+$ définie par

$$ \varphi\left(\begin{pmatrix} x_1 \\ y_1 \\ z_1 \end{pmatrix}, \begin{pmatrix} x_2 \\ y_2 \\ z_2 \end{pmatrix}\right) = (x_1 - x_2)\tanh(x_1 - x_2) $$

est une pseudo-métrique pour laquelle $\rho(m) = m$.

Nous avons les propriétés suivantes.

**Proposition 21.** ✓ La suite de réelles $(\rho(m))_{m \ge 1}$ est croissante.

✓ Pour tout $m \in \mathbb{N}^*$, $\rho(m) \le \rho(2)^{m-1}$ (voir [219]).

✓ Soit $\varphi_1$ et $\varphi_2$ deux pseudo-métriques sur D tels que $\rho_1(m)$ et $\rho_2(m)$ sont les plus petits réels vérifiant (4.2). Pour tout $\alpha > 0$ et $\beta > 0$, l'application $\alpha\varphi_1 + \beta\varphi_2$ est une pseudo-métrique sur D vérifiant $\rho(m) = \text{Max}\{\rho_1(m), \rho_2(m)\}$.

La définition d'une pseudo-métrique permet de définir une pseudo-longueur entre deux sommets de G. Considérons $\varphi$ une pseudo-métrique définie sur l'ensemble des variables des oscillateurs. $z_i$ étant la variable d'état du ième sommet de G et $\varphi(z_i, z_j)$ la pseudo-longueur entre les sommets $i$ et $j$ alors $\varphi(z_{i_1}, z_{i_2}) + \dots + \varphi(z_{i_{m-1}}, z_{i_m})$ est la pseudo-longueur de n'importe quel chemin $P_{i,j} = (i = i_1, i_2, \dots, i_m = j)$ du sommet $i$ au sommet $j$.

Les résultats de synchronisation du système (4.1) que nous avons établis (théorèmes 3 et 4 de la section 4.2.2) s'appuient sur un résultat de majoration de la somme des pseudo-longueurs entre deux oscillateurs quelconques. Nous avons montré que cette somme est bornée par la somme des pseudo-longueurs des chemins joignant deux oscillateurs quelconques. Cette majoration est définie à une constante multiplicative près $C(G)$. La constante $C(G)$ joue un rôle important dans les théorèmes 3 et 4 puisque la force de synchronisation $\epsilon$ est proportionnelle à cette constante. Le résultat de majoration est le suivant :

**Proposition 22.** Soit G un graphe connexe, $\mathcal{E}$ l'ensemble de ses sommets et $\varphi$ une pseudo-métrique sur l'ensemble D. Pour toute sommet $i$, soit $z_i \in D$ un vecteur associé au sommet $i$. Il existe une constante $C$ dépendant seulement de G telle que

$$ \sum_{i,j} \varphi(z_i, z_j) \le C \sum_{(i,j) \in \mathcal{E}} \varphi(z_i, z_j). \quad (4.3) $$

De plus, le plus petit réel C satisfaisant (4.3) et noté $C(G)$ est borné par

$$ \frac{n(n-1)}{2} \delta(G) \rho(\delta(G)), \quad (4.4) $$

où $\delta(G)$ est le diamètre de G.

**Remarque 25.** 1. La borne (4.4) de $C(G)$ peut ne pas conduire à une bonne estimation de $C(G)$ pour un graphe particulier. Néanmoins, cette borne est valable pour tout graphe avec n sommets.

2. Dans le cas d'une pseudo-métrique satisfaisant l'inégalité triangulaire, i.e. quand $\rho(n) = n$ pour tout $n \in \mathbb{N}^*$, une méthode considérant G comme entrée et donnant une borne de $C(G)$ est proposée dans [15]. Les deux principales étapes de la méthode proposée sont :

(a) pour tout $(i,j)$ avec $i > j$, choisir un chemin $P_{i,j}$ ; celui-ci est généralement choisi de longueur minimum (nombre d'arêtes du chemin) ;

(b) pour chaque arête e du graphe, déterminer la somme $B(e)$ des longueurs de tous les chemins $P_{i,j}$ contenant l'arête e. Une borne pour $C(G)$ est alors $\text{Max}\{B(e) : e \in \mathcal{E}\}$.

Pour chaque choix de chemin, ces deux étapes retournent une borne de $C(G)$. Clairement, le nombre de chemins possibles est important ces choix ne sont pas optimaux en général. Toutefois, les calculs des bornes de $C(G)$ restent possibles. En modifiant la première étape de la méthode, elle peut-être appliquée à notre cas : elle consiste à considérer, pour tout chemin $P_{i,j}$, la pseudo-longueur $\rho(|P_{i,j}|)$ plutôt que la longueur $|P_{i,j}|$.

3. Dans le cas de pseudo-métriques $\varphi$ satisfaisant $\rho(m) = m$, la méthode proposée dans [15, 16] et les bornes explicites de $C(G)$ calculées avec cette méthode peuvent être directement utilisées. C'est le cas de la seconde fonction dans l'exemple 4.2.1.
4.2.2 Synchronisations complètes

Nous avons traité deux cas. Le premier était la synchronisation complète pour laquelle les oscillateurs $X_1, \dots, X_n$ sont dans $D = \mathbb{R}^d$, le second étant la synchronisation complète pour laquelle les oscillateurs sont dans un voisinage de $D$ pour la variété $X_1 = X_2 = \dots = X_n$.

Les hypothèses sur le système (4.1) sont les suivantes :

— Pour tout $(i,j) \in \mathcal{E}$, il existe des réels positifs $a_1, \dots, a_d$ tels que pour tout $(X_i, X_j) \in D$,

$$ \varphi(X_i, X_j) = \sum_{k=1}^{d} a_k (X_i^k - X_j^k) h^k(X_i, X_j) \quad (4.5) $$

sont des pseudo-métriques où $h = (h^1, \dots, h^d)^T$ est la fonction de synchronisation.

— Pour tout $(i,j) \in [1,n]^2$ et, pour tout $t \ge t_0$ où $t_0 \in \mathbb{R}$, si $(X_i, X_j) \in D$, alors nous avons

$$ \sum_{k=1}^{d} a_k (X_i^k - X_j^k) (F_i^k(X_i, t) - F_j^k(X_j, t)) \le \varphi(X_i, X_j). \quad (4.6) $$

— Pour tout $(i,j) \in [1,n]^2$ et tout $(X_i, X_j) \in D$, si

$$ \sum_{k=1}^{d} a_k (X_i^k - X_j^k) (F_i^k(X_i, t) - F_j^k(X_j, t)) = 0 \quad \text{et/ou } \varphi(X_i, X_j) = 0 \text{ alors } (X_i = X_j). \quad (4.7) $$

**Remarque 26.** 1. L'hypothèse (4.5) implique que,

$$ h(X_i, X_j) = -h(X_j, X_i) \text{ (anti-symétrique).} $$

2. L'hypothèse (4.7) est nécessaire pour prouver la synchronisation du système (4.1) dans les théorèmes 3 et 4. La condition $\varphi(X_i, X_j) = 0$ dans cette hypothèse n'est pas toujours suffisante car elle n'implique pas les égalités de toutes les composantes des oscillateurs. Dans ce cas, la première condition est nécessaire pour prouver la synchronisation complète.

Dans les cas pratiques, une première difficulté est de prouver l'existence des trajectoires du système (4.1) pour $t$ suffisamment grand. La proposition suivante lie l'existence des trajectoires des systèmes couplés et non couplés.

**Proposition 23.** Pour tout $(i,j) \in [1,n]^2$, supposons que les hypothèses (4.5), (4.6) et (4.7) sont satisfaites et que, pour tout $t \ge t_0$, $X_i^T F_i(X_i, t) \le \Psi(||X_i||)$ où $\Psi$ satisfait les conditions

$$ \int_{t=t_0}^{+\infty} \frac{ds}{\Psi(t)} = +\infty \text{ et } \Psi(t) > 0 \text{ pour tout } t \ge t_0 \ge 0. $$

Alors, le problème de Cauchy défini par le système (4.1) et la condition initiale $$ \begin{pmatrix} X_1(t_0) \\ \vdots \\ X_n(t_0) \end{pmatrix} \in \mathbb{R}^{nd} \text{ a une solution définie sur } [t_0; +\infty). $$

Nous avons établi le théorème suivant qui permet d'obtenir la synchronisation complète du système (4.1). Sa preuve est basée sur l'application de la seconde méthode de Lyapunov et sur les propriétés de la pseudo-métrique $\varphi$.

**Théorème 3.** Supposons que les hypothèses (4.5), (4.6) et (4.7) soient satisfaites pour $D = (\mathbb{R}^d)^2$. Si $\epsilon > \frac{C(G)}{2n}$, où $C(G)$ est la borne optimale telle que l'inégalité (4.3) est vraie, alors le système (4.1) synchronise complètement.

Nous avons également établi un résultat de synchronisation locale. Pour cela, nous avons eu besoin des

notations suivantes. Soit $H$ la matrice diagonale $\text{Diag}(a_1, \dots, a_d)$ et $\mathcal{H} = \begin{pmatrix} H & 0 & \cdots & 0 \\ 0 & H & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & H \end{pmatrix}$ la matrice composée de $\frac{n(n-1)}{2}$ matrices $H$. L'application
$$ \| \cdot \|_V : \begin{array}{rcl} \mathbb{R}^{\frac{n(n-1)}{2} d} & \rightarrow & \mathbb{R}^+ \\ X & \rightarrow & \sqrt{\frac{1}{2} X^T H X} \end{array} \qquad (4.8) $$

est une norme puisque $a_1, \dots, a_d$ sont positifs. Considérons

$$ V(t) = \|\Delta(t)\|_V^2 = \frac{1}{2} \sum_{k=1}^{d} \sum_{i<j \le n} a_k (X_i^k(t) - X_j^k(t))^2. $$

**Théorème 4.** *Soit $\mathcal{B}$ la boule fermée*

$$ \{ X \in \mathbb{R}^{\frac{n(n-1)}{2} d} \mid \|X\|_V \le r \} $$

où $r$ est un réel positif. Supposons que les hypothèses (4.5), (4.6) et (4.7) sont vérifiées quand $\Delta$ appartient à $\mathcal{B}$ et supposons qu'à un instant $t_0$, $\Delta(t_0)$ appartient à l'intérieur $\mathring{\mathcal{B}}$ de $\mathcal{B}$.

Si $\epsilon > \frac{C(G)}{2n}$, $C(G)$ étant la borne optimale telle que l'inegalité (4.3) a lieu, alors le système (4.1) se synchronise.

### 4.2.3 Application

Le théorème 4 est appliqué à un réseau d'oscillateurs de Chua [42, 43, 131]. Si $X = (x, y, z)^T$, les équations pour un simple oscillateur sont données par $\dot{X} = F(X)$ où

$$ F(x, y, z) = \begin{pmatrix} a[y - x - f(x)] \\ x - y + z \\ -by - cz \end{pmatrix}, $$

$a > 0, b > 0, c > 0$ et $f$ est une fonction par morceaux définie par

$$ f(x) = dx + 1/2(d - e)(|x + 1| - |x - 1|). $$

Considérons une configuration en étoile pour laquelle le vecteur d'état du i-ème oscillateur est donné par

$$ \begin{pmatrix} \dot{x}_i \\ \dot{y}_i \\ \dot{z}_i \end{pmatrix} = \begin{pmatrix} a[y_i - x_i - f(x_i)] \\ x_i - y_i + z_i \\ -by_i - cz_i \end{pmatrix} + \epsilon \sum_{j \mid (i,j) \in \mathcal{E}} \begin{pmatrix} \tanh(x_i - x_j) \\ 0 \\ 0 \end{pmatrix}, \qquad (4.9) $$

et $\mathcal{E} = \{(1;2), (1;3), \dots, (1;n)\}$.

**Remarque 27.** 1. La fonction de synchronisation tanh modélise la situation commune dans laquelle l'influence des oscillateurs sur les autres satures quand elle est trop grande. Ce type de fonction apparaît naturellement dans les systèmes dynamiques [182].

2. Le fait que les solutions du système sont définis sur $\mathbb{R}$ est une conséquence directe de la proposition 23.

Soit $\delta$ un réel positif tels que $2d - e - 1 < \delta$. Remarquons que $2d - e$ borne l'ensemble des pentes

$$ \left\{ \frac{f(x) - f(y)}{x - y} \mid 0 < |x - y| \le 1 \right\}. $$

Le système (4.9) est équivalent à

$$ \begin{pmatrix} \dot{x}_i \\ \dot{y}_i \\ \dot{z}_i \end{pmatrix} = \begin{pmatrix} a[y_i - x_i - f(x_i)] \\ x_i - y_i + z_i \\ -by_i - cz_i \end{pmatrix} $$
$$ +\epsilon' \sum_{j|(i,j) \in \mathcal{E}} \left( \begin{array}{c} \delta \tanh(x_i - x_j) \\ 0 \\ 0 \end{array} \right), \qquad (4.10) $$

où $\epsilon' = \epsilon/\delta^1$. Dans le but d'appliquer le théorème 4, vérifions les hypothèses (4.5), (4.6) et (4.7). Pour cela, considérons la boule $\mathcal{B} = \{X \in \mathbb{R}^{\frac{3n(n-1)}{2}} | \|X\|_V \le 1\}$ et $a_1 = 1/a$, $a_2 = 1$, $a_3 = 1/b$. Comme $a > 0$ et $b > 0$, l'hypothèse (4.5) est évidente.

Vérifions l'hypothèse (4.6). Nous avons

$$ \begin{align} & \sum_{k=1}^{3} a_k (X_i^k - X_j^k) (F_i^k(X_i) - F_j^k(X_j)) \tag{4.11} \\ & = (x_i - x_j)(f(x_i) - f(x_j)) - (x_i - x_j)^2 \nonumber \\ & \qquad - (y_i - y_j)^2 - c/b(z_i - z_j)^2 \nonumber \\ & \le (2d - e - 1)(x_i - x_j)^2. \nonumber \end{align} $$

De la définition de $\delta$, nous en déduisons un voisinage fermé $\mathcal{B}$ de la variété de synchronisation,

$$ \begin{align} & \sum_{k=1}^{3} a_k (X_i^k - X_j^k) (F_i^k(X_i) - F_j^k(X_j)) \tag{4.12} \\ & \le \delta(x_i - x_j) \tanh(x_i - x_j) \nonumber \end{align} $$

Par conséquent, l'hypothèse (4.6) est vérifiée.

Montrons l'hypothèse (4.7). Supposons que $\varphi(x_i, x_j) = 0$, nous avons $x_i = x_j$ car tanh est décroissante. L'égalité (4.11) implique que $X_i = X_j$.

L'application $\varphi$ est une pseudo-métrique dans un voisinage de la variété de synchronisation satisfaisant, pour tout $m \in \mathbb{N}^*$, $\rho(m) = m$ (voir par exemple 4.2.1) et la borne $C(G)$ est donnée explicitement par $2n - 3$ (voir la remarque 25 et [16]).

Le théorème 4 peut être appliqué : si, pour un instant $t_0$, le système (4.9) est suffisamment proche de la variété de synchronisation, c'est à dire si $\Delta(t_0) \in \mathcal{B}$, et si $\epsilon > \frac{\delta C(G)}{2n} = \frac{\delta (2n - 3)}{2n}$ alors le système (4.9) synchronise complètement.

**Remarque 28.** L'approche que nous avons utilisée dans cet exemple peut être facilement étendue à d'autres fonctions f et à d'autres termes de synchronisation [109].

## 4.3 Liens entre synchronisation et bifurcations

Le modèle HR décrit en partie le fonctionnement d'un neurone et résulte de la simplification et la généralisation du modèle de Hodgkin-Huxley [95, 96]. Le modèle HR est un système autonome de trois équations différentielles ordinaires pouvant reproduire plusieurs dynamiques du neurone comme les spikes ou les comportements chaotiques. Il s'écrit

$$ \left\{ \begin{aligned} \dot{x} &= y + ax^2 - x^3 - z + I \\ \dot{y} &= 1 - dx^2 - y \\ \dot{z} &= \epsilon(b(x - c_x) - z) \end{aligned} \right. \qquad (4.13) $$

Les paramètres $a, b$ et $d$ sont déterminés expérimentalement. Le paramètre $I$ correspond au courant appliqué tandis que $c_x$ est la solution $x$ du système de dimension 2 défini par les deux premières équations de (4.13) quand $I = 0$ et $z = 0$. Finalement, le paramètre $\epsilon$ représente le rapport d'échelle de temps entre les flux rapides et lents qui ont lieu à travers la membrane d'un neurone. Les deux premières équations contrôlent la dynamique rapide du système tandis que la troisième la dynamique lente.

Les comportements périodiques sont des phénomènes fréquents dans beaucoup de systèmes naturels et sont souvent liés à des paramètres de bifurcation du modèle mathématique. Or, étant donné la difficulté des calculs,

<footnote n="1">1. Le réel $\delta$ est introduit pour borner la somme $\sum_{k=1}^{3} a_k (X_i^k - X_j^k) (F_i^k(X_i) - F_j^k(X_j))$ par la pseudo-métrique $\varphi(X_i, X_j) = \delta(x_i - x_j) \tanh(x_i - x_j)$ comme dans l'hypothèse (4.6)</footnote>
les analyses de bifurcation du modèle HR sont essentiellement numériques [59, 86, 185, 212] sauf dans [28] mais où une partie des paramètres est fixée.

Le paramètre $\epsilon$ joue un rôle important pour décrire l'activité neuronale et c'est ce paramètre que nous avons choisi comme paramètre de bifurcation, comme dans [49, 86]. Une bifurcation de Hopf réfère à un changement radical de la dynamique de sortie du modèle dû au changement de la valeur d'un des paramètres. Près d'un point critique, selon la valeur du paramètre de bifurcation, la sortie du modèle tend soit vers un point d'équilibre stable soit vers une solution périodique [90]. À la figure 4.1 est représenté le portrait de phase dans $(x,y,z)$ et les séries temporelles du système (4.13) pour deux valeurs différentes de $\epsilon$, avant et après la bifurcation de Hopf. Les valeurs des autres paramètres sont $a = 3, b = 4, D = 5$ et $I = 3,25$. Lorsque la valeur de $\epsilon$ est plus petite que $\epsilon_c = 0,125912$, le système admet une solution périodique. Lorsque sa valeur croît, cette solution disparaît, puis le système admet un point d'équilibre stable quand sa valeur est plus grande que $\epsilon_c$.

<figure><img src="image_38.png" /><figcaption>FIGURE 4.1 – (a) Portrait de phase dans $(x,y,z)$ et séries temporelles du système (4.13) avec les paramètres $a = 3, b = 4, D = 5$ et $I = 3,25$ et $\epsilon = 0.12 < \epsilon_c$. Ici, a solution asymptotique est un cycle limite stable. (b) Portrait de phase dans $(x,y,z)$ et séries temporelles du système (4.13) avec les paramètres $a = 3, b = 4, D = 5$ et $I = 3,25$ et $\epsilon = 0.13 > \epsilon_c$. Pour cette valeur de $\epsilon$ la solution asymptotique est un noeud stable.</figcaption></figure>

Dans [49, 86], les auteurs ont fait une étude de bifurcation sur le simple modèle HR par rapport à $\epsilon$ et en donnant une valeur numérique aux autres paramètres. Dans [50], nous avons proposé une étude de bifurcation par rapport à $\epsilon$ quelque soit les valeurs des autres paramètres. De plus nous avons étendu l'étude à un système couplé de N neurones. Un de nos défis était d'étudier la dynamique d'un réseau de neurones, comme l'émergence de comportements coordonnés du réseau, à partir des dynamiques d'un seul neurone.

En ce qui concerne les fonctions de couplage, nous avons choisi les fonctions linéaires. En effet, dans le cerveau, les neurones interagissent grâce aux synapses [87] qui peuvent être soient électriques soient chimiques. Dans le premier cas, les interactions sont modélisées par des fonctions de couplage linéaires, dans le second, par des fonctions de couplage non linéaires. Notre choix s'est porté sur les synapses électriques. Plus généralement, les réseaux que nous avons considérés sont des réseaux de neurones acycliques, orientés de modèles HR couplés linéairement. Notre choix de topologie inclut des topologies classiques si le réseau a un noeud source et pas de cycles (voir la figure 4.2).

Au niveau du réseau, nous avons travaillé sur un point d'équilibre particulier du réseau. Si les premiers systèmes découplés des autres commencent à osciller dues à la valeur critique de bifurcation, les autres sous-systèmes développeront également des cycles limites avec une amplitude dépendant de la force de couplage et le réseau présentera une bifurcation de Hopf. Nous avons été plus loin dans le raisonnement et avons montré qu'il est suffisant qu'un seul sous-système oscille pour avoir une bifurcation de Hopf du système complet. Dans le cas d'un réseau acyclique direct, nous avons également établi :

- une expression analytique du paramètre de bifurcation d'un sous-système en fonction des coordonnées du point d'équilibre et des paramètres de ce sous-système et de la somme des forces de couplage dépendant des systèmes précédents,

- les formules donnant la direction, la stabilité et la période des solutions périodiques de toute bifurcation de Hopf dans le réseau.
<figure><img src="image_39.png" /><figcaption>FIGURE 4.2 – Chaîne de réseau (si $k_{ij} = 0 \forall j \neq i - 1$); fermeture transitive de la chaîne (et tout état intermédiaire de cette fermeture; réseau en étoile (si $k_{ij} = 0 \forall j \neq 1$); tout graphe direct acyclique avec un noeud source.</figcaption></figure>

Ce travail théorique a permis d'étudier l'influence du couplage sur les bifurcations de Hopf dans le réseau
dans deux cas particuliers.

4.4 Conclusion

Les réseaux d'oscillateurs sont depuis toujours un des paradigmes les plus utilisés pour décrire les processus
dynamiques qui se produisent dans de nombreux domaines de la science. Étant donné la complexité des systèmes
étudiés, la plupart des résultats sont numériques. Nous avons pu toutefois généraliser des résultats théoriques
de synchronisation et proposer des expressions analytiques pour étudier des bifurcations de Hopf du réseau à
partir d'une bifurcation (de Hopf) survenant dans un neurone. Dans le prochain chapitre, je traite d'exemples
de réseaux modélisant une zone géographique impactée par une catastrophe et la propagation des différents
comportements adoptés par la population.
# Chapitre 5

## Croisement entre les mathématiques et les sciences humaines

Cette partie présente un travail qui a commencé en 2013 avec le projet PEPS HuMaIn (Humanité-Mathématiques-sciences de l'Information) : *Le comportement de foule en situation de catastrophe : les comprendre pour mieux les prévenir*. Celui-ci visait à décrire et modéliser les différents comportements adoptés lors d'une catastrophe en s'appuyant sur les connaissances en neurosciences. Il a conduit au développement du modèle Panique-Control-Réflexe (PCR) et à la thèse de Guillaume Cantin. Les différents intervenants que nous avons rencontrés au cours de ces années comme les géographes, les professionnels de gestion de crise et les psychologues ont fait évoluer notre compréhension sur les types de comportements et leurs transitions au cours d'un événement brutal. Ce travail a abouti au projet ANR Com2SiCa (https://www.com2sica.cnrs.fr/). Ainsi, le modèle PCR a été enrichi et est devenu le modèle Alerte-Panique-Contrôle (APC). Ce dernier est basé sur deux variables clés utilisées en psychologie pour identifier et caractériser les différents comportements humains : l'intensité de la charge émotionnelle et sa régulation. Cette caractérisation a conduit à définir trois états comportementaux adoptés lors d'une catastrophe : alerté, contrôlé et paniqué ainsi que les transitions entre ces états.

Un effort particulier de réflexion a été fait pour simplifier la représentation des différents comportements tout en gardant les caractéristiques principales des dynamiques individuelles et collectives dans ce qui est sans doute un domaine assez complexe de l'expérience humaine, celui du comportement humain en situation d'urgence. Ainsi, les déterminants psychologiques capturés sur la base des discussions avec des psychologues, d'enquêtes de terrain et de simulations ont été traduits en systèmes relativement simples d'équations différentielles ordinaires non linéaires. Cet effort s'est déroulé en deux étapes successives concordant au développement de deux modèles : le modèle PCR puis le modèle APC. Tous les deux ont été étudiés sous leur forme adimensionnelle. Pour ces deux modèles, nous nous sommes essentiellement concentrés sur un événement unique même si dans les premiers travaux sur le modèle PCR, nous avons intégré un processus d'enchaînement de catastrophes dans les transitions entre les trois états paniqués, contrôlés et réflexes.

La difficulté dans ces projets a été de valider les modèles mathématiques étant donné le peu de données de terrain. Toutefois, pour le modèle APC, des connaissances des thématiciens, nous avons pu donner un ordre de grandeur des paramètres. Puis, nous avons travaillé sur une expérience de réalité virtuelle réalisée en lieu clos et contrôlé au Laboratoire de Psychologie des Pays de la Loire (Nantes) en 2020. Elle s'appuyait sur un scénario, celui de l'arrivée d'une vague de tsunami d'origine sismique sur une plage de Nice [19]. Celui-ci est désormais intégré en France dans les outils réglementaires nationaux. L'étude portait plus spécifiquement sur l'analyse des réactions psychologiques, comportementales et physiologiques (mesures physiologiques du stress) des participants volontaires. Le scénario intégrait également des personnages virtuels (avatars) afin d'observer les réactions dans un contexte d'interaction sociale (entre la personne volontaire et les personnages virtuels). Le calage des données topographiques de la plage (largeur de la plage, présence d'escaliers, de galets, rambardes d'accès etc...) a été réalisé par des mesures précises réalisées durant les enquêtes de terrain à Nice [205]. Actuellement, nous ne disposons pas de données pour confronter et apprécier la différence entre réalité et réalité virtuelle. Cependant, dans les expériences de réalité virtuelle (RV) faites dans d'autres domaines des sciences humaines et sociales, il y a une bonne validité écologique, c'est-à-dire que les réactions humaines sont similaires entre la réalité et les expériences de RV comme expliqué dans [118]. Ce modèle ouvre donc de nouvelles perspectives dans la compréhension et la modélisation des comportements en situation de catastrophe comme nous le verrons dans la suite.
Modèles dynamiques

<figure><img src="image_40.png" /><figcaption>FIGURE 5.1 – Analyse et comparaison des deux approches de modélisation comportementale basée, soit sur des descriptions individuelles conduisant à des systèmes multi-agents (bas de la figure), soit sur des description de populations conduisant à des systèmes dynamiques (haut de la figure). L’approche à base d’agents nécessite une génération d’individus avec une part d’aléatoire puisqu’on ne peut se référer à des connaissances ou données qu’à l’échelle des populations (C. Bertelle, 2021)</figcaption></figure>

5.1 Pourquoi avons nous choisi la modélisation macroscopique à base d'EDOs ?

Deux échelles de modèles sont envisageables pour représenter les dynamiques et les interactions comportementales : au niveau individuel (approche micro utilisant généralement des modèles à base d'agents) et au niveau collectif (approche macro utilisant généralement des modèles à base de systèmes dynamiques) (voir Figure 5.1).

Dans la première phase de projet, deux points essentiels nous ont amenés à choisir la modélisation à base de systèmes dynamiques au détriment des modèles micro ou agent. Tout d'abord, dans les modèles micro ou agent, il est nécessaire d'introduire une part d'aléatoire dans le comportements des individus ; l'attribution individuelle des comportements est parfois difficile à justifier scientifiquement. Ces modèles, qui donnent l'apparence d'une grande précision, conduisent à un très grand nombre de degrés de liberté empêchant tout contrôle des dynamiques afin de comprendre les processus qui les génèrent. Ensuite, l'un des objectifs du projet était le développement d'une plateforme web cartographique de simulation des réactions humaines destinée à la communauté scientifique mais aussi aux professionnels en charge de la sécurité et sûreté des citoyens. Or, les politiques et notamment des politiques du risque souhaitent gérer et réguler les situations au niveau collectif et non de contrôler chaque individu. Ils agissent à une échelle macroscopique, celle de la population dans son intégralité. En raison de ces limites, nous n'avons pas retenu la modélisation à base d'agents mais une modélisation à base de systèmes dynamiques. Ces derniers permettent d'intégrer les connaissances que l'on a sur les propriétés collectives d'une population et sur le(s) processus(s). De plus, le paramétrage est bien adapté à la gestion des dynamiques comportementales collectives en situation de catastrophe. Par exemple, il suffit de modifier de façon continue un paramètre jugé important, comme celui de la culture du risque, pour en déduire comment varie la dynamique des comportements humains collectifs. Ces modèles sont donc intégrateurs de connaissances. Ils peuvent donner une aide à la formation dans une situation d'urgence et permettre de mieux comprendre le comportement émotionnel d'un groupe d'individus conduisant entre autres à des effets de panique.
## 5.2 Le modèle Panic-Control-Reflex (PCR)

Cette partie reprise des travaux [30, 32, 33, 164, 193, 199] est le fruit d'un travail en commun avec Guillaume Cantin, des géographes (Damienne Provitolo, Edwige Dubos-Paillard et Anne Tricot), des mathématiciens (Valentina Lanza, Moulay Aziz Alaoui) et des informaticiens (Rodolphe Charrier, Cyrille Bertelle).

### 5.2.1 D'une typologie des comportements au modèle mathématique

Dans ces projets, nous avons fait abstraction des dynamiques de déplacements collectifs [35, 47, 93, 130, 186, 223] et nous nous sommes focalisés sur des aspects souvent moins détaillés dans les modèles présents dans la littérature. Ces aspects correspondent à la modélisation des états comportementaux, leurs évolutions dans une situation de catastrophe et les phénomènes de propagation sociale de ces états. En effet, la plupart des modèles mathématiques font référence seulement aux paniques collectives [70, 106, 155] car ces réactions sont plus facilement observables et retranscrites tant lors de catastrophes naturelles, technologiques que sanitaires [52, 85, 94, 163, 165]. À une époque où les politiques et responsables territoriaux sont juridiquement responsables des incidents que peuvent provoquer la gestion de foules et des grands rassemblements, ces réactions de panique sont d'autant plus craintes par les acteurs en charge de la gestion des catastrophes. En effet, ils ne peuvent pas facilement les contrôler car elles se déclenchent et se diffusent comme un "feu de paille". Ces comportements collectifs de panique ne sont pas uniquement une réaction à une situation critique mais ils peuvent être l'origine de la dangerosité d'une situation qui aurait pu être contrôlée avec un comportement adapté comme lors de la panique au stade Heysel en 1985 ou lors des rassemblements à la Mecque, notamment sur le pont Djamarat en 2015. Pour décrire de manière plus précise ces comportements collectifs, les seules comportements de panique ne sont pas suffisants.

Des comportements contrôlés sont également relatés [68, 133, 153] : certains illustrent des réactions statiques comme l'immobilité ; d'autres des réactions dynamiques comme la fuite réfléchie, l'éloignement ou encore l'entraide [69, 78, 166]. Enfin, les recherches en psychologie et en neurosciences ont pu enseigner que l'être humain restait rarement figé dans un même état comportemental mais exprimait un enchaînement de réactions du début à la fin d'une catastrophe.

Provitolo et al. [164] ont proposé une typologie des comportements humains dans les zones impactées par un événement catastrophique. Cette typologie est basée sur les zones cérébrales impliquées dans la réponse comportementale mais également durant les phases temporelles de l'événement et ses alertes. Dans une zone impactée, deux catégories de comportements ont été décrites selon la zone cérébrale impliquée :

- la zone reptilienne du cerveau, le centre de nos instincts et la base de nos émotions. Il est la source des comportements instinctifs qui induisent des comportements impulsifs et d'urgence;

- la zone préfontale du cerveau qui conduit à des comportements contrôlés. Il adapte de manière plus réfléchie les réactions aux stimuli extérieurs.

La première catégorie regroupe les comportements d'évasion et de panique, mais aussi les comportements comme la sidération et les comportements d'automate [52]. Ce mécanisme permet aux gens de réagir rapidement soit en s'enfuyant le plus vite possible, soit en étant sidérés et en étant physiquement incapables de se déplacer dans l'espace. La deuxième catégorie concerne les comportements calmes et de maîtrise de soi.

Dans ce premier modèle mathématique, nous avons décidé de distinguer trois types de comportements différents en situation de catastrophe soudaine : deux faisant partie de la première catégorie et une faisant partie de la seconde. Malgré la diversité des comportements réflexes, nous avons décidé de les diviser en deux types.

- les comportements réflexes hors panique (sidération et comportements d'automates).

- les réactions de panique puisqu'elles ont un statut particulier dans les comportements réflexes.

Le comportement de panique a été différentié car il est le plus redouté : il peut provoquer des situations dangeruses dans une foule, dues au piétinement et à l'écrasement et peut provoquer des décès. De plus, il est difficile de stopper ce mécanisme une fois enclenché puisque il est plus lié à des dynamiques internes qu'à l'éloignement du danger [52]. On peut noter que ce comportement n'est pas toujours adopté comme, par exemple, lors d'un tremblement de terre au Japon, un pays où la population est formée pour réagir de façon adéquate en cas de catastrophe.

Enfin, puisque le cortex préfrontal s'empare du cerveau reptilien, le troisième type de comportement envisagé regroupe l'ensemble des réactions contrôlées, intelligentes et raisonnées. Ils peuvent prendre différentes formes lors d'une catastrophe, comme par exemple l'évacuation ou la fuite, le confinement, la mise à l'abri, la recherche de secours, le pillage, le vol, etc... Malgré leur variété, nous avons décidé de les globaliser.

Il faut bien noter qu'automaticité ne veut pas dire irrationalité. Le fait qu'une réaction soit automatique, instinctive, n'indique pas qu'elle soit irrationnelle, inadaptée ou erronée. Le même constat s'applique pour les comportements contrôlés, ces derniers n'étant pas systématiquement adaptés à la situation.

Les trois comportements précédents ne se produisent pas tous en même temps et respectent un certain ordre. En effet, le premier comportement d'un individu face au danger est un comportement réflexe, suivi dans un
second temps, d'un comportement contrôlé ou paniqué [83].

Cette catégorisation de ces comportements et leurs interactions a conduit à la création du modèle PCR (Panic-Control-Reflex) qui est décrit de façon plus précise à la section suivante.

## 5.2.2 Le modèle mathématique PCR et ses applications

Soit la classe $Q$ composée d'individus dans un comportement du quotidien. On suppose qu'il n'y a pas de naissance ou de décès au cours de la catastrophe. La population est donc constante et est composée de $N$ individus. Ainsi, on note

- $Q(t)$ le nombre d'individus ayant un comportement du quotidien avant que la catastrophe commence, donc $Q(0) = N$,

- $B(t)$ le nombre d'individus revenus à un comportement du quotidien durant ou après la catastrophe. Nous supposons qu'à la fin de la catastrophe, tous les individus sont dans ce comportement, donc $B(t_{end}) = N$.

D'après la section précédente, la population est décomposée en 3 sous-populations lors de la catastrophe, d'où les variables :

- $R(t)$ = nombre de personnes dans un comportement réflexe,

- $C(t)$ = nombre de personnes dans un comportement contrôlé,

- $P(t)$ = nombre de personnes dans un comportement paniqué.

Puisque nous supposons être en présence d'un événement soudain et imprévisible, toute la population concernée aura d'abord un comportement réflexe, correspondant à des comportements instinctifs. Ainsi, les comportements routiniers, représentés ici avec la variable $Q(t)$, ne peuvent être transformés qu'en comportements réflexes. Par la suite, les comportements réflexes peuvent devenir des comportements contrôlés ou paniqués. Puisque $B(t)$ représente le nombre d'individus qui reviennent à un mode de vie normal, il ne peut être alimenté que par les comportements contrôlés $C(t)$. En effet, un individu a besoin de reprendre le contrôle de soi afin de retrouver la routine quotidienne. De plus, nous supposons qu'une fois revenu à la normale, il reste dans leurs comportements habituels. Ainsi, les individus dans $B$ ne peuvent pas passer dans $Q$ et rentrer dans la boucle.

Selon les psychologues et les géographes [24, 52, 91, 164], lors d'événements catastrophiques, les individus ne restent pas dans le même état comportemental. Des transitions, récurrentes à chaque catastrophe, ont lieu entre ces sous-populations : les processus de causalité ($B_1, B_2, C_1, C_2$) et les processus d'imitation et de contagion ($\alpha_1, \alpha_2, \beta_1, \beta_2, \gamma_1, \gamma_2$) qui seront expliqués dans la suite. Ces états comportementaux et ces processus sont représentés à la figure 5.2. L'événement extérieur, c'est à dire la catastrophe est représentée par la fonction $\gamma$

<figure><img src="image_41.png" /><figcaption>FIGURE 5.2 – Graphe modélisant les trois types de comportements humains dans un contexte de catastrophe qui peut être une fonction discrète ou continue selon le type de catastrophe. Par exemple, un événement tels</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>FIGURE 5.3 – Représentation de la fonction η</figcaption></figure>

qu’un tsunami peut être modélisé par une fonction discrète alors qu’une inondation peut-être modélisée par une fonction continue puisqu’elle peut-être annoncée quelques heures avant son début. La fonction $\varphi$ modélise le retour au comportement du quotidien des personnes en comportement de contrôle.

$$
\left\{
\begin{aligned}
\frac{dR}{dt} &= \gamma(t)Q\left(1 - \frac{R}{N}\right) - (B_1 + B_2)R + \tilde{F}(R, C)\frac{R}{N}C + \tilde{G}(R, P)\frac{R}{N}P, \\
\frac{dC}{dt} &= B_1R + C_1P - C_2C - \tilde{F}(R, C)\frac{R}{N}C + \tilde{H}(C, P)\frac{C}{N}P - \varphi(t)C\left(1 - \frac{B}{N}\right), \\
\frac{dP}{dt} &= B_2R - C_1P + C_2C - \tilde{G}(R, P)\frac{R}{N}P - \tilde{H}(C, P)\frac{C}{N}P, \\
\frac{dQ}{dt} &= -\gamma(t)Q\left(1 - \frac{R}{N}\right), \quad \frac{dB}{dt} = \varphi(t)C\left(1 - \frac{B}{N}\right).
\end{aligned}
\right.
\qquad (5.1)
$$

La population est supposée constante, c'est à dire : $N = R + C + P + Q + B$.

Les fonctions d'imitation sont définies de la façon suivante :

$$
\left\{
\begin{aligned}
\tilde{F}(R, C) &= -\alpha_1\eta\left(\frac{R}{C+\varepsilon}\right) + \alpha_2\eta\left(\frac{C}{R+\varepsilon}\right) \\
\tilde{G}(R, P) &= -\beta_1\eta\left(\frac{R}{P+\varepsilon}\right) + \beta_2\eta\left(\frac{P}{R+\varepsilon}\right) \\
\tilde{H}(C, P) &= -\gamma_1\eta\left(\frac{C}{P+\varepsilon}\right) + \gamma_2\eta\left(\frac{P}{C+\varepsilon}\right).
\end{aligned}
\right.
\qquad (5.2)
$$

où $\alpha_i > 0, \beta_i > 0, \gamma_i > 0, i = 1, 2, 0 < \varepsilon \ll 1$ et $\eta(s) = \frac{s^2}{1+s^2}$.

Les fonctions décrivant le processus d'imitation entre deux comportements distincts dépendent de la proportion d'individus dans chacun des comportements. En particulier, la fonction d'imitation $\tilde{H}$ décrit une imitation qui peut se faire entre les $P$ et $C$ mais également entre les $C$ et $P$. Selon les valeurs des paramètres $\gamma_1$ et $\gamma_2$, elle peut décrire la loi du plus fort. En effet, lorsque la valeur de $s$ est petite, celle de $\eta(s)$ l'est également (voir la figure 5.3). Ainsi, s'il y a peu de contrôlés et beaucoup de paniqués, le ratio $C/(P+e)$ est petit et le ratio $P/(C+e)$ est grand. On en déduit d'après la forme de la courbe $\eta$ que le terme $\eta(C/(P+e))$ est petit par rapport à $\eta(P/(C+e))$. Cela signifie que si $\gamma_1$ et $\gamma_2$ sont égaux, peu de personnes en état de panique imiteront les comportements contrôlés alors que beaucoup de personnes en état de contrôle imiteront les paniqués et ceci dans une proportion $\gamma_1$.

Le paramètre $\epsilon$ est un paramètre qui modélise le début du processus d'interaction dû à l'imitation qui ne se produit que pour une certaine valeur des sous-populations $R$ et $C$.
Par la suite, nous avons a-dimensionné le système qui s’écrit alors :

$$
\left\{
\begin{aligned}
\frac{dr}{dt} &= \gamma(t)q(1-r) - (B_1+B_2)r + F(r,c)rc + G(r,p)rp, \\
\frac{dc}{dt} &= B_1r + C_1p - C_2c - F(r,c)rc + H(c,p)cp - \varphi(t)c(1-b), \\
\frac{dp}{dt} &= B_2r - C_1p + C_2c - G(r,p)rp - H(c,p)cp, \\
\frac{dq}{dt} &= -\gamma(t)q(1-r), \quad \frac{db}{dt} = \varphi(t)c(1-q)
\end{aligned}
\right.
\qquad (5.3)
$$

avec

— $q(t) = \text{densité de personnes ayant un comportement du quotidien avant que la catastrophe commence,}$
  donc $q(0) = 1$

— $b(t) = \text{densité de personnes qui sont revenus à un comportement du quotidien durant ou après la}$
  catastrophe

— $r(t) = \text{densité de personnes dans un comportement réflexe}$

— $c(t) = \text{densité de personnes dans un comportement contrôlé}$

— $p(t) = \text{densité de personnes dans un comportement paniqué}$.

et

$$
\left\{
\begin{aligned}
F(r, c) &= -\alpha_1 \eta \left( \frac{r}{c+\epsilon} \right) + \alpha_2 \eta \left( \frac{c}{r+\epsilon} \right) \\
G(r, p) &= -\beta_1 \eta \left( \frac{r}{p+\epsilon} \right) + \beta_2 \eta \left( \frac{p}{r+\epsilon} \right) \\
H(c, p) &= -\gamma_1 \eta \left( \frac{c}{p+\epsilon} \right) + \gamma_2 \eta \left( \frac{p}{c+\epsilon} \right).
\end{aligned}
\right.
\qquad (5.4)
$$

La population étant supposée constante, le système se ramène à l'étude d'un système de dimension 4 dont les variables sont $(r, c, p, q)$. En posant $u = (r, c, p, q)^T$ et en introduisant l'opérateur $\phi(u,t) = (\phi_i(u,t))_{1\le i\le 4}$, le système (5.3) se ré-écrit

$$
\frac{du}{dt} = \phi(u,t), \quad t \ge t_0. \tag{5.5}
$$

À ce système est ajouté la condition initiale $u_0 = (0,0,0,1)$.

Les paramètres du système (5.3) sont regroupés dans le vecteur $\lambda = (B_1, B_2, C_1, C_2, \alpha_1, \alpha_2, \delta_1, \delta_2, \mu_1, \mu_2)$ dont le domaine est donné par $\Lambda = (\mathbb{R}^{+*})^2 \times (\mathbb{R}^{+})^{10}$.

Pour toute valeur du paramètre $\lambda \in \Lambda$, et pour toute condition initiale $u_0 \in (\mathbb{R}^+)^5$, il existe $T > t_0$ tel que le problème (5.3) possède une unique solution locale $u$ définie et continûment dérivable sur $[t_0, T]$ dont les composantes sont positives. De plus, l'ensemble compact $\kappa = \{r, c, p, q, b) \in (\mathbb{R}^+)^4; r+c+p+q+b \le 1\}$ est positivement invariant pour le flot induit par le problème (5.5) et les solutions sont globales.

**Exemple de fonctions $\gamma$ et $\varphi$ prises dans les simulations :**

Dans le cas d'une catastrophe soudaine, comme un tsunami local, on a supposé que la population commençait
à être informée rapidement c'est-à-dire après 1 minute, et que toute la population concernée était informée
dans les 3 minutes suivantes, d'où la forme de la fonction $\gamma$ présentée à la figure 5.4. Le retour au quotidien
correspondant à la fonction $\varphi$, ne peut pas être immédiat. Nous avons supposé qu'il se faisait au bout de 5
minutes après le début de la catastrophe. Ce retour est supposé se faire très lentement, ce qui nous a amené à
considérer la fonction $\varphi$ comme à la figure 5.4. La forme des courbes peut être modulée selon le type d'événement
catastrophique (selon qu'il soit annoncé ou non).

**Remarque 29.** Il a été difficile de calibrer ce modèle car peu de données ont été trouvées dans la littérature. S. Boyd [24] estime que dans la plupart des catastrophes : "15% des individus manifestent des réactions pathologiques patentes, 15% conservent leur sang froid et 70% manifesteraient un comportement apparemment calme mais répondant en fait à un certain degré de sidération émotionnelle et de perte d'initiative qui relève du registre pathologique". Un rapport de l'ingénierie de la sécurité incendie [1] estime que 10 à 15% de personnes ont des comportements anormaux ou présentant des états confusionnels, de 12 à 25% réagissent correctement et 75% des personnes sont pétrifiées, inhibées ou hébétées. Enfin, J.S. Tyhurst [190] fait état, en situation exceptionnelle, de 10 à 25% d'individus qui présentent des comportements inadaptés, 12 à 25% qui réagissent de façon adaptée,
<figure><img src="image_43.png" /><figcaption>FIGURE 5.4 – Fonctions $\varphi$ et $\gamma$</figcaption></figure>

et 75% qui ont une réaction apparemment normale mais dans la réalité non réfléchie. Ces différentes sources révèlent une diversité des comportements. Nous avons opté pour la fourchette [164, 194]

— $r(t)=50$ à 75% de la population

— $c(t)=12$ à 25% de la population

— $p(t)=12$ à 25% de la population.

Ces pourcentages ont été une des justification pour notre choix de travailler avec le système a-dimensionné.

Ces différentes réactions sont de durée variable [206], avec :

— $r(t)$ : la durée du comportement réflexe varie de quelques minutes à 1 heure. Le plus souvent ce type de comportement ne dépasse pas 15 minutes. Mais il peut durer plus longtemps notamment lorsqu'il correspond au délai d'évacuation d'une région sinistrée. Dans ce cas on voit alors réapparaitre progressivement un retour aux comportements d'entraide et de recherche des proches et des sinistrés [52].

— $p(t)$ : la durée du comportement de panique varie de quelques minutes à 1 heure. Le plus souvent ce type de comportement ne dépasse pas 15 minutes. La panique collective se résout généralement spontanément. Cependant, parfois une intervention énergique extérieure permet à la population paniquée $p(t)$ de retrouver un état calme, de prostration, ce qui explique que cette dernière repasse par un comportement d'automate $r(t)$ avant d'adopter un comportement contrôlé $c(t)$.

— $r(t)+p(t)$ : la durée du comportement émotionnel non contrôlé ne dure pas plus d'1h30 en général. Dans le cadre de ce modèle, nous partons en effet de l'hypothèse qu'un individu ne peut pas rester 1 heure dans un comportement réflexe et encore 1 heure dans un état de panique.

— $c(t)$ : la durée du comportement contrôlé varie de quelques minutes à plusieurs heures, selon la rapidité et l'efficacité des acteurs institutionnels.

Lors des simulations, les paramètres du modèle ont été choisis pour que ces conditions soient vérifiées.

Le paramètre $C_1$ a un rôle particulier dans le système. La linéarisation du système PCR conduit à une valeur propre nulle. Cette stabilité est fortement liée à une bifurcation qui a lieu quand $C_1$ s'annule. La proposition suivante résume les résultats sur les points d'équilibre et leur nature [32].

**Proposition 24.** — Si $C_1 > 0$, alors le système (5.3) admet $O = (0,0,0,0)$ comme unique point d'équilibre.

— Si $C_1 = 0$ alors le système admet tout point $(0,0,\bar{p},0)$ avec $\bar{p} > 0$ comme point d'équilibre, en plus de $O$. Par la suite, cette valeur critique sera appelée persistance de panique.

— Le point d'équilibre $O$ est localement stable pour toute valeur de $\lambda \in \Lambda$ mais pas globalement stable.

Le dernier point a été montré en étudiant la fonction de Lyapunov correspondant à l'énergie liée au potentiel démographique de la population affectée par l'événement catastrophique étudié. Ainsi, le système admet deux classes de solutions, selon si la valeur du paramètre d'évolution $C_1$ du comportement de panique au comportement de contrôle est positive ou nulle. Dans le premier cas, c'est-à-dire $C_1 > 0$, la solution converge vers l'équilibre trivial, ce qui correspond à la situation attendue lorsque tous les individus reprennent un comportement du quotidien. Dans le deuxième cas, lorsque $C_1 = 0$, le système présente une persistance de panique, qui est la situation à éviter, au risque d'aggraver l'impact de l'événement catastrophique. Dans ce dernier cas, nous avons montré que les termes d'imitation quadratique peuvent inhiber ou exacerber la persistance de la panique, selon la direction du processus d'imitation de la panique au contrôle ou du contrôle à la panique, respectivement. Dans [194], nous avons traité l'exemple d'un tremblement de terre au Japon. Puisqu'au Japon, la culture du risque est bien établie, la population est entraînée à réagir rapidement. Ainsi, les processus de causalité de l'état
de réflexe et de l'état de panique à l'état de contrôle sont plus importants par rapport aux autres processus, même pour une zone définie à l'épicentre du tremblement de terre.

La prise en compte du relief géographique de la zone impactée par la catastrophe conduit naturellement à l'étude des réseaux. Des questions spécifiques se posent dans l'étude des systèmes complexes définis par de tels réseaux couplés, parmi lesquels la "synchronisation" ou la relation entre la topologie du réseau et la dynamique du système. Dans l'article de conférence [33], un exemple concernant le risque particulier de tsunami sur le littoral méditerranéen a été étudié. Le découpage géographique a été fait avec Damienne Provitolo et Edwige Dubos-Paillard, le scénario choisi étant celui d'un tsunami à Nice. À la figure 5.5 est représenté le découpage géographique de la plage de Nice et la promenade des anglais. Le couplage entre les différents noeuds est donné à la figure 5.6. Mathématiquement, les zones géographiques ont été définies selon la valeur du paramètre de bifurcation $C_1$. Une zone de danger est modélisée avec la valeur de paramètre $C_1 = 0$ et une zone de sécurité par $C_1 > 0$.

<figure><img src="image_44.png" /><figcaption>FIGURE 5.5 – Simulation d'un tsunami sur la côte méditerranéenne. En rouge, la zone de danger, en jaune la promenade des anglais, en vert, l'intérieur de la ville qui correspond à une zone de sécurité. Les arêtes entre les noeuds représentent les déplacements possibles des individus.</figcaption></figure>

Nous avons montré dans ce dernier article que l'évacuation des zones à risques correspondant aux plages, vers les zones refuges situées en centre-ville, joue un rôle déterminant. Par exemple, un couloir bouché peut provoquer une persistance de la panique. À l'opposé, un chemin d'évacuation supplémentaire de la plage vers le centre-ville peut aider les individus à reprendre leur comportement du quotidien. Dans [29], ces modèles ont été généralisés à des réseaux quelconques.

## 5.3 Le modèle Alert-Panic-Control (APC)

Cette partie est la suite du travail présenté à la section 5.2 qui a conduit au projet ANR Com2SiCa (https://www.com2sica.cnrs.fr/). Il reprend les articles [121, 201]. L'enjeu du projet était

1. en s'appuyant sur les connaissances des psychologues, avoir une meilleure connaissance de la diversité des comportements, la dynamique temporelle propre à chacun d'entre eux mais aussi les changements d'états comportementaux, c'est-à-dire le fait d'enchaîner plusieurs réactions pour faire face à un évènement catastrophique

2. de développer un modèle mathématique incluant ces connaissances

3. de développer une plate-forme pour l'analyse et la simulation comportementale des populations urbaines en situation de catastrophe.

Le modèle PCR (Panic-Control-Reflex) est basé sur les sciences cognitives et la modélisation mathématique des maladies infectieuses. Il décrit la dynamique des différents comportements lors d'événements catastrophiques. Ce modèle ne se focalise pas uniquement sur la panique, mais présente une analyse et une classification précise des comportements en fonction des différentes régions cérébrales activées et des différentes transitions comportementale. Toutefois, cette approche est basée sur les neurosciences et peut être insuffisante pour décrire les interactions sociales intervenant dans les évolutions comportementales. C'est dans cette optique que des psychologues ont été intégrés au projet. Leur contribution a permis d'intégrer dans la modélisation des populations
<figure><img src="image_45.png" /></figure>

FIGURE 5.6 – Représentation graphique du couplage linéaire. Chaque personne reste dans son état comportemental d'un noeud à un autre. Les forces de couplage ε traduisent l'intensité du flux et leurs valeurs dépendent du comportement.

en situation de catastrophe la cognition et la contagion sociale<fnref n="1" /> [77, 79] qui ont été appréhendées sous l'angle de la charge et la régulation émotionnelle [143]. Cette approche a permis de formaliser le modèle mathématique que l'on notera par la suite APC (Alert-Panic-Control). Celui-ci a été construit et enrichi suite à des

- Analyses de vidéos,

- Focus groupes avec des acteurs institutionnels et opérationnels intervenant sur les scènes de catastrophes (attentats, guerres, catastrophes naturelles, industrielles, etc.)<fnref n="2" />

- Entretiens immersifs auprès de la population sur des lieux réels, à Nice au mois de septembre 2018 et au Havre au mois de mai 2019 (45 minutes d'entretien individuel, totalisant 60 personnes enquêtées sur les 2 terrains). Ces protocoles ont permis d'immerger les participants dans une ambiance sonore et visuelle (un scénario de tsunami d'origine sismique à Nice, un scénario d'accident technologique au Havre), d'observer, enregistrer et identifier sur les deux terrains d'enquête leurs réactions comportementales, de cartographier les dynamiques de déplacements empruntées par les personnes pour sortir de la situation imaginée ou projetée, de mesurer sur une échelle qualitative la charge émotionnelle de chaque enquêtés [116].

L'ensemble des traitements réalisés sur ces différents corpus a permis de produire des données qualitatives, quantitative et spatialisées utiles à la modélisation mathématique et notamment sur les types de comportements, les changements d'états comportementaux et les cartographies des trajectoires de déplacement [120, 119]. Ces données ont pu être exploitées pour préciser le modèle APC.

<footnote n="1">1. La cognition sociale désigne l'ensemble des processus cognitifs impliqués dans les interactions sociales. La contagion sociale signifie qu'un grand nombre de personnes et pas seulement un individu imitent le comportement dominant par comparaison sociale. Le modèle dominant de comportement va ainsi se propager d'une personne à une autre, un peu comme s'il s'agissait d'un virus.</footnote>

<footnote n="2">2. Le Service départemental d'Incendie et de Secours des Alpes-Maritimes, la Préfecture des Alpes-Maritimes, la Direction Départementale des Territoires et de la Mer des Alpes - Maritimes, la Brigade des Sapeurs-Pompiers de Paris, la Police Nationale, la Police Municipale de Nice, l'Agence Régionale de la Santé, la Métropole Nice Côte d'Azur et la ville de Nice pour la Direction des Risques Majeurs, la Communauté Urbaine le Havre Seine Métropole pour la Direction des Risques Majeurs, l'Union Nationale des Centres Communaux d'Action Sociale</footnote>
### 5.3.1 Catégories comportementales

Pour qualifier et catégoriser la diversité des comportements humains<fnref n="3" /> deux variables importantes en psychologie des émotions ont été utilisées : la charge ou l'intensité émotionnelle (niveau de stress, nervosité, peur, etc.) et la régulation émotionnelle (la capacité des personnes à contrôler cette charge).

Ce croisement de variables a permis de qualifier les réponses ou les comportements dans des situations de catastrophe dans un éventail plus large que le seul comportement de panique (Figure 5.7). Ainsi, trois états comportementaux ont été dégagés :

- Les états d'alerte (A) désignent un ensemble de micro-comportements observables d'un point de vue moteur (sursaut, mouvement oculaire rapide, interrogation du regard/visuelle ou verbale des personnes proches). L'état d'alerte correspond à une phase d'assimilation, il marque une rupture avec le comportement habituel : un moment où la personne va rechercher dans un temps très court de l'information sur le scénario qu'elle est en train de vivre. En état d'alerte, la personne est en état d'incertitude, la charge émotionnelle est faible même s'il peut y avoir des signaux corporels marquant cette rupture avec le comportement quotidien, comme par exemple l'augmentation du rythme cardiaque.

- Les états ou comportements de contrôle (C) sont des comportements réfléchis qui désignent un ensemble de réactions dont la charge émotionnelle, plus ou moins forte, est régulée. Les comportements de contrôle sont divers : ils peuvent désigner des comportements pro-sociaux (entraide, protection, soin d'autrui) mais également des comportements anti-sociaux (vols, pillages, voyeurisme etc.). En d'autres termes, la notion de contrôle désigne le fait que la personne est capable de réguler ses émotions, d'agir et adapter son comportement au contexte de crise. Toutefois, cela ne veut pas dire que le comportement assure la sécurité de la personne ou sa survie, ni qu'il est vertueux, social ou exemplaire au regard d'une norme sociale.

- L'ensemble des comportements de panique (P) désigne des comportements non contrôlés où les émotions liées à la peur ont pris le dessus. Un comportement de panique implique une charge émotionnelle forte et une régulation faible, peu efficace pour retrouver un état contrôlé. Différents comportements peuvent signifier un état de panique : fuite panique, stupeur, agitation désordonnée [52].

L'intensité de la charge et de la régulation émotionnelle a donc permis de qualifier ces états, selon l'approche suivante :

- États d'alerte : une charge émotionnelle faible et une régulation émotionnelle faible.

- États de contrôle : une charge émotionnelle faible ou forte et une régulation émotionnelle forte.

- États de panique : une charge émotionnelle forte et une régulation émotionnelle faible.

Ces trois ensembles de comportements Alertés, Paniqués et Contrôlés donnent lieu à différentes intensités de charges émotionnelles et sont indépendantes du type de catastrophe. Il est à noter qu'ils peuvent être différentiés par des mesures physiologiques (voir [3] pour la détection du stress dans le cas d'un travail informatique). Par exemple, la panique provoque une augmentation significative du rythme cardiaque et de la transpiration contrairement aux comportements d'alerte.

<footnote n="3">3. Plus de 60 comportements humains ont été identifiés[176].</footnote>
<figure><img src="image_46.png" /><figcaption>FIGURE 5.7 – Classification des différents comportements durant un événement catastrophique selon la charge émotionnelle et sa régulation (N. Verdière).</figcaption></figure>

Il est connu qu'au cours d'une catastrophe, les personnes n'adoptent pas le même comportement et que des successions de réactions différentes peuvent être observées. Par exemple, les personnes adoptent un état d'alerte dès que la catastrophe est perçue, puis un comportement de panique ou de contrôle selon leur expérience et leur culture du risque.

Sur la base de la catégorisation proposée par les géographes et les psychologues du projet, les hypothèses thématiques suivantes ont été posées :

- Quelque soit la source d'information (aléas, cris), on passe en état d'alerte.

- L'état d'alerte peut être très rapide, il se situe le plus souvent tout au début de l'événement, cependant selon les caractéristiques de ce dernier, par exemple s'il revêt des répétitions comme dans les situations de séismes, l'état d'alerte peut revenir dans un deuxième temps, notamment après un état de contrôle, mais cela est plus rarement observé.

- C'est seulement en passant par un état de contrôle qu'on pourra retrouver un comportement pseudo-normal du quotidien, comportement qui ne peut pas être retrouvé tant qu'on est dans un état de panique.

- Même pour la panique, le pic émotionnel redescend. La panique peut donc évoluer vers le comportement d'alerte ou le comportement de panique.

- Le retour a un état d'alerte est un processus intrinsèque. Il ne peut pas se faire par imitation, et cela que l'on soit en état de panique ou en état de contrôle.

- Nous considérons les transitions comportementales selon deux processus distincts : l'agrégation d'une part, l'imitation et la propagation d'autre part. Ces processus présentés ci-après nous permettent d'étudier la dynamique comportementale :

✓ les transitions intrinsèques liées aux caractéristiques individuelles (âge, expérience, culture du risque...). Les comportements vont devenir collectifs par agrégation quand il y a une synchronisation de comportements individuels au même moment sans qu'il y ait nécessairement une action de coordination volontaire ou d'interaction du collectif.

✓ les transitions dues à l'imitation et la propagation qui nous permettent de comprendre l'émergence de comportements collectifs liés à l'interaction au sein d'une population ou d'une foule. Seuls les comportements de type alerte ne sont pas imitables car, de part leur nature et leur brièveté, ils ne sont pas imités par les populations en état de panique et de contrôle. Ces processus de propagation et d'imitation permettent d'étudier la diffusion des comportements dans le temps et dans l'espace.

- Malgré les entretiens et les focus group menés auprès d'une diversité d'acteurs, la connaissance des changements de comportements qui seraient liés à la présence de personnes décédées sur une scène de catastrophes est encore incertaine. C'est pourquoi, nous posons comme hypothèse qu'il n'y a pas de boucle de rétroaction de l'impact de la mortalité sur les comportements.
Dans la section suivante, je présente le modèle APC construit à partir de ces différentes hypothèses thématiques.

## 5.3.2 Présentation du modèle APC (Alert-Panic-Control)

Rappelons que $A$ décrit l'ensemble des personnes en état d'alerte, $C$ en état de contrôle et $P$ en état de panique. De plus, les comportements du quotidien avant la catastrophe sont pris en compte ($Q$ par la suite) ainsi que les retours au quotidien (notés $B$).

Le diagramme décrivant les interactions entre les différents comportements est donné à la Figure 5.8. Les transitions sont notées :

- $B_i, i = 1, \dots, 4$, et $C_i, i = 1, 2$ pour les transitions naturelles qui sont spécifiques à chaque individu selon son expérience et son passé;

- $\tilde{F}, \tilde{G}$ et $\tilde{H}$ pour les transitions dues à l'imitation ou à la comparaison sociale.

Les fonctions $\gamma$ et $\varphi$ modélisent l'aléa et le retour au quotidien.

<figure><img src="image_47.png" /><figcaption>FIGURE 5.8 – Schéma du modèle APC.</figcaption></figure>

De la même façon que pour le modèle PCR, le modèle APC a été ré-écrit en terme de densité (la population est supposée constante) :

$$ 
\left\{
\begin{aligned}
\frac{da(t)}{dt} &= \gamma(t)q(t) - (B_1 + B_2 + D_a)a(t) + B_3c(t) + B_4p(t) - F(a(t),c(t))a(t)c(t) - G(a(t),p(t))a(t)p(t), \\
\frac{dc(t)}{dt} &= B_1a(t) + C_1p(t) - (B_3 + C_2 + D_c)c(t) + F(a(t),c(t))a(t)c(t) + H(c(t),p(t))c(t)p(t) - \varphi(t)c(t), \\
\frac{dp(t)}{dt} &= B_2a(t) + C_2c(t) - (B_4 + C_1 + D_p)p(t) + G(a(t),p(t))a(t)p(t) - H(c(t),p(t))c(t)p(t), \\
\frac{dq(t)}{dt} &= -\gamma(t)q(t), \quad \frac{db(t)}{dt} = \varphi(t)c(t).
\end{aligned}
\right.
\qquad (5.6) $$

Les fonctions d'imitation sont construites de la même façon que pour le modèle PCR, la différence provient des fonctions $F$ et $G$ qui sont unidirectionnelles :

$$ 
\left\{
\begin{aligned}
F(a, c) &= \alpha\eta\left(\frac{c}{a+\varepsilon}\right), \quad G(a, p) = \beta\eta\left(\frac{p}{a+\varepsilon}\right) \\
H(c, p) &= -\gamma_1\eta\left(\frac{p}{c+\varepsilon}\right) + \gamma_2\eta\left(\frac{c}{p+\varepsilon}\right),
\end{aligned}
\right.
\qquad (5.7) $$

où $\alpha > 0, \beta > 0, 0 < \varepsilon \ll 1$ et $\eta(s) = \frac{s^2}{1+s^2}$ qui est représenté à la figure 5.3.

## 5.3.3 Une première calibration à partir des connaissances des thématiciens

Les sciences humaines disposent de très peu de données quantitatives dans ce domaine. Toutefois, il est reconnu que les phénomènes d'imitation sont le plus souvent observés dans les cas de foules denses. Ainsi, pour
<table><caption>TABLE 5.1 – Principales notations pour le système APC en terme de densité.</caption><thead><tr><th>Fonctions</th><th>Notation</th></tr></thead><tbody><tr><td>Comportement du quotidien avant la catastrophe</td><td>q(t)</td></tr><tr><td>Comportement d'Alerte</td><td>a(t)</td></tr><tr><td>Comportement du Contrôle</td><td>c(t)</td></tr><tr><td>Comportement de Panique</td><td>p(t)</td></tr><tr><td>Comportement du quotidien après la catastrophe</td><td>b(t)</td></tr><tr><td>Commencement de la catastrophe</td><td>γ(t)</td></tr><tr><td>Retour à un comportement du quotidien</td><td>φ(t)</td></tr><tr><td>Imitation entre alerte et contrôle</td><td>F</td></tr><tr><td>Imitation entre alerte et panique</td><td>G</td></tr><tr><td>Imitation entre panique et contrôle</td><td>H</td></tr></tbody></table><table><thead><tr><th>Paramètres</th><th>Notation</th></tr></thead><tbody><tr><td>Évolution de alerte à contrôle</td><td>B<sub>1</sub></td></tr><tr><td>Évolution de alerte à panique</td><td>B<sub>2</sub></td></tr><tr><td>Évolution de contrôle à alerte</td><td>B<sub>3</sub></td></tr><tr><td>Évolution de panique à alerte</td><td>B<sub>4</sub></td></tr><tr><td>Évolution de panique à contrôle</td><td>C<sub>1</sub></td></tr><tr><td>Évolution de contrôle à panique</td><td>C<sub>2</sub></td></tr><tr><td>Taux de mortalité</td><td>D<sub>R</sub>, D<sub>P</sub>, D<sub>C</sub></td></tr></tbody></table>

ce type de population, l'accent doit être mis sur les interactions d'imitation et non sur les transitions intrinsèques. Les coefficients des fonctions d'imitation seront compris entre 0,5 et 1 et dans le cas d'une population clairsemée, ces paramètres seront inférieurs à 0,5. Dans l'hypothèse d'une population à faible culture de risque, les transitions intrinsèques vers la panique sont les plus importantes. On a donc supposé $B_2 > B_1$, c'est-à-dire les transitions intrinsèques de l'alerte à panique sont plus importantes que celles d'alerte à contrôle et $C_2 > C_1$ c'est-à-dire que les transitions intrinsèques de contrôle à panique sont plus importantes que celles de panique à contrôle. Une plage de paramètres peut déjà être définie à partir de ces remarques. Elle est donnée dans la table 5.2.

<table><caption>TABLE 5.2 – Différent scénarios selon les valeurs des paramètres</caption><thead><tr><th>culture du risque<br>population</th><th>faible</th><th>forte</th></tr></thead><tbody><tr><th>dense</th><td><u>Jeux 1</u> : B<sub>2</sub> > B<sub>1</sub>,<br>C<sub>2</sub> > C<sub>1</sub><br>α, β, γ<sub>c→p</sub>, γ<sub>p→c</sub> > 0.5</td><td><u>Jeux 3</u> : B<sub>2</sub> < B<sub>1</sub>,<br>C<sub>2</sub> < C<sub>1</sub><br>α, β, γ<sub>c→p</sub>, γ<sub>p→c</sub> > 0.5</td></tr><tr><th>éparpillée</th><td><u>Jeux 2</u> : B<sub>2</sub> > B<sub>1</sub>,<br>C<sub>2</sub> > C<sub>1</sub><br>α, β, γ<sub>c→p</sub>, γ<sub>p→c</sub> < 0.5</td><td><u>Jeux 4</u> : B<sub>2</sub> < B<sub>1</sub>,<br>C<sub>2</sub> < C<sub>1</sub><br>α, β, γ<sub>c→p</sub>, γ<sub>p→c</sub> < 0.5</td></tr></tbody></table>

Pour vérifier que le modèle reproduit les observations des géographes et psychologues, pour chaque ensemble de paramètres, nous avons comparer les flux correspondant à l'imitation et aux transitions intrinsèques. Pour cela, les fonctions suivantes ont été comparées : ($t \in [0, T]$) :
$$ Im_p(t) = \gamma_1 \xi \left( \frac{p(t)}{c(t) + \epsilon} \right) p(t)c(t) + \beta \xi \left( \frac{p(t)}{a(t) + \epsilon} \right) p(t)a(t), \quad (5.8) $$

$$ In_p(t) = B_2 a(t) + C_2 c(t). \quad (5.9) $$

$Im_p(t)$ représente la densité des individus qui à chaque instant $t$ adopte un comportement de panique par imitation, alors que $In_p(t)$ représente la densité de ceux qui adoptent un comportement de panique dû aux transitions intrinsèques. La comparaison de ces deux fonctions permet de vérifier quelle transition entre le processus d'imitation et le processus intrinsèque est la plus importante.

Par analogie, nous notons $(t \in [0, T])$ :

$$ Im_c(t) = \alpha \xi \left( \frac{c(t)}{a(t) + \epsilon} \right) c(t) a(t) + \gamma_2 \xi \left( \frac{c(t)}{p(t) + \epsilon} \right) p(t) c(t). \quad (5.10) $$

$$ In_c(t) = B_1 a(t) + C_1 p(t) \quad (5.11) $$

Ainsi, $Im_c(t)$ représente la densité des individus qui adoptent un comportement contrôlé par imitation à chaque instant $t$ alors que $In_c(t)$ représente la densité de population qui adopte un comportement contrôlé dû aux transitions intrinsèques.

En choisissant des jeux de paramètres définis par le tableau 5.2, nous avons trouvé que :

- **Scénario 1, Jeux 1 (exemple d'une foule dense avec une faible culture du risque)** : les personnes qui adoptent un comportement de panique l'adoptent par imitation. De plus, certaines valeurs de paramètres permettent de faire converger le système vers un comportement contrôlé.

- **Scénario 2, Jeux 2 (exemple d'une foule éparpillée avec une faible culture du risque)** : les transitions intrinsèques prévalent sur les imitations et le comportement de panique est le plus adopté. En comparant les scénarios 1 et 2, nous avons obtenu plus ou moins les mêmes résultats qualitatifs : une prédominance de la panique. Or ce comportement de panique est la conséquence de deux processus différents : dans le premier cas la prévalence de la panique est due à des processus d'imitation, dans le second cas à des transitions intrinsèques.

- **Scénario 3, Jeux 3 (une population dense avec une forte culture du risque)** : les personnes qui adoptent un comportement contrôlé l'adoptent surtout par un processus d'imitation.

- **Scénario 4, Jeux 4 (une population éparse avec une forte culture du risque)** : les transitions intrinsèques prévalent sur les transitions d'imitation.

Ces résultats sont en accord avec les observations des psychologues : dans le cas d'une population dense, les comportements sont surtout adoptés par imitation.

Une étude des points d'équilibre du système lorsque $\varphi \equiv 0$ a montré le rôle des paramètres sur la dynamique du système. Pour chaque scénario, l'espace des paramètres peut être partitionné en plusieurs régions selon le nombre de points d'équilibre et selon qu'ils conduisent à un état contrôlé ou paniqué. Ce travail théorique explique par exemple les résultats du scénario 1 par rapport aux connaissances des psychologues. Dans le cas d'une foule dense avec une faible culture du risque, il y aura surtout une prédominance de comportements paniqués mais qu'ils ne sont pas automatiques.

# 5.4 Une calibration à partir de la réalité virtuelle

## 5.4.1 Introduction

Pour mieux comprendre les différents comportements vécus lors d'un événement catastrophique et les étudier dans des conditions idéales, une expérience de réalité virtuelle a été menée au laboratoire de Psychologie de Nantes en 2020. Cette expérience consistait à immerger le volontaire dans une catastrophe soudaine et imprévue sans signe avant-coureur. Grâce à cette nouvelle approche, des données d'observation ont été obtenues et ont permis d'obtenir un ordre de grandeur des paramètres. Le modèle calibré peut alors être un outil intéressant pour capturer et comprendre des dynamiques comportementales non directement mesurables lors d'un événement catastrophique.

Cependant, trois problèmes se sont posés lors de l'estimation des paramètres :

- La présence de non-linéarités dans le modèle mathématique dues aux processus d'imitation;

- La difficulté de différencier les états comportementaux même si l'expérience a été menée en laboratoire;

- Le faible nombre de données.

Ce contexte réaliste mais difficile implique une analyse fine préalable de faisabilité de l'estimation des paramètres, c'est pourquoi, nous avons dû répondre à deux questions principales avant d'estimer les paramètres :
— Est-il possible d'estimer l'ensemble des paramètres du modèle à partir des mesures physiologiques fournies par la réalité virtuelle ?

— Quel est le nombre minimal de mesures nécessaires pour être sûr des résultats de la procédure d'estimation de paramètres ?

Les réponses ont été apportées par trois études théoriques. La première est une analyse de sensibilité, la seconde une étude d'identifiabilité. La troisième a déterminé le nombre minimal de mesures de temps nécessaires pour estimer les paramètres inconnus. Une fois ces trois études menées sur le modèle mathématique, une procédure d'estimation des paramètres a été mise en place. Cette dernière est basée sur la résolution d'un problème des moindres carrés non linéaire résolu avec l'algorithme génétique [92]. Un algorithme global a été préféré à un algorithme local car il n'y avait aucune estimation initiale des valeurs des paramètres.

Les études étant conditionnées par l'expérience de RV elle-même, je vais tout d'abord la présenter et faire le lien avec le modèle mathématique APC puis je présenterai comment des trois analyses nous avons mis en place une procédure d'estimation de paramètres.

## 5.4.2 L'expérimentation de réalité virtuelle (RV)

L'expérience de RV a été menée par des psychologues du laboratoire LPPL de Nantes (France) dans un environnement clos et contrôlé. Son objectif était de déterminer à partir de mesures physiologiques, les différents états comportementaux et leurs changements survenant lors d'un événement catastrophique en prenant en compte le contexte environnemental tel que le type d'événement, la configuration spatiale et la présence ou l'absence d'emois.

### 5.4.2.1 La procédure de RV et l'échantillonnage

Cette expérience de réalité virtuelle consistait à immerger des individus dans la simulation d'un tsunami sur le littoral méditerranéen de Nice (France), un scénario événement qui a depuis été intégré aux outils réglementaires nationaux. Le scénario de la simulation était basé sur des données de recherche (voir [18]). Pour construire le contexte spatial, des informations topographiques ont été recueillies à Nice comme la plage telles que sa largeur, la présence d'escaliers, de galets et garde-corps d'accès etc... (voir [205]). Des personnages virtuels (emois) ont également été inclus pour observer les réactions comportementales dans le contexte des interactions sociales entre les répondants à l'enquête et les emois.

Trois scénarios ont été construits à partir de trois conditions expérimentales conçues pour simuler différentes réactions, les comparer et mieux comprendre comment le contexte affecte le comportement. Par exemple, une situation de danger ne génère pas automatiquement un sentiment de menace si la personne est seule. La réaction dépend de l'expérience de la personne et des comportements de ceux autour desquels le comportement peut être imité. Les trois scénarios suivants ont été réalisés :

— Premier scénario (scénario neutre) : dans cette condition expérimentale, la personne est seule sans personne virtuelle.

— Deuxième scénario (scénario organisé) : dans cette condition expérimentale, la personne est en présence de personnes virtuelles au comportement adapté (échapper au tsunami en marchant rapidement dans la bonne direction et sans crier).

— Troisième scénario (scénario non organisé) : dans cette condition expérimentale, la personne est en présence de personnes virtuelles au comportement inadapté (fuite en courant dans tous les sens et en criant).

Dans le deuxième scénario, les réactions physiques des emis expriment un comportement de contrôle alors que dans le troisième ils expriment la panique.

Quatre-vingt-huit volontaires non rémunérés ont été recrutés dans la presse locale et ont été affectés au hasard à l'une des trois conditions. L'échantillon était composé de 43 femmes et 45 hommes. Leur âge variait de 18 à 56 ans avec une moyenne de 26,76 ans et un écart type de 9,41. Les participants ont été répartis dans les différentes conditions en respectant un équilibre entre hommes et femmes, et une diversité dans les âges. Par conséquent, 29 personnes ont été affectées au premier et au troisième scénario et 30 au deuxième. La répartition des participants selon leur sexe et le scénario est fournie dans le tableau 5.3.

Lors de la simulation, chaque participant était immergé dans un environnement de réalité virtuelle qui fournit un retour visuel et auditif (figure 5.9). Le sujet était directement placé sur la plage au début de la simulation et avait des mouvements libres. Pour capter et maintenir l'attention du participant, la chronologie proposée s'est basée sur la perception d'un tsunami sur la côte méditerranéenne de Nice plutôt que sur la réalité de sa dynamique. Comme décidé, après 60 secondes, le niveau de la mer baissait. Ensuite, la vague se formait à l'horizon et les bruits de l'eau devenaient plus forts. La vague atteignait la plage après 3 minutes. Si les participants restaient immobiles sur la plage, l'application s'arrêtait lorsqu'ils étaient heurtés par la vague. Si les participants avaient un comportement de fuite à l'intérieur de la ville, l'application s'arrêtait lorsqu'ils atteignaient une sortie de ville. Dans tous les autres cas, l'application s'arrêtait 2 minutes après que la vague
<table><caption>TABLE 5.3 – Distribution des enquêtés selon le sexe et le scénario.</caption><thead><tr><th rowspan="2">Condition</th><th colspan="2">Sexe</th><th rowspan="2">Total</th></tr><tr><th>Femmes</th><th>Hommes</th></tr></thead><tbody><tr><td>Premier scénario</td><td>14</td><td>15</td><td>29</td></tr><tr><td>Scénario organisé</td><td>15</td><td>15</td><td>30</td></tr><tr><td>Scénario non organisé</td><td>14</td><td>15</td><td>29</td></tr><tr><td>Total</td><td>43</td><td>45</td><td>88</td></tr></tbody></table>

<figure><img src="image_48.png" /><figcaption>FIGURE 5.9 – Photographies de l'expérience virtuelle au laboratoire de Nantes. Chaque participant a été immersé dans un environnement de réalité virtuelle fournissant un retour visuel et auditif.</figcaption></figure>
ait submergé la plage. Les étapes sont résumées à la figure 5.10.

<figure><img src="image_49.png" /><figcaption>FIGURE 5.10 – Déroulement du scénario d’un tsunami dans le temps : en vert, progression de la vague due au tsunami ; en violet, les sons entendus par les personnes participant à l’expérience de réalité virtuelle ; en orange, l’action des avatars ; en rouge, le début et la fin des actions.</figcaption></figure>

### 5.4.2.2 Mesures

Des mesures physiologiques ont été enregistrées pour chaque participant afin de déterminer leur charge émotionnelle à chaque instant de la simulation. L'appareil Nexus X MKII était connecté au participant par des électrodes. Il a enregistré pour chacun d'eux un électrocardiogramme (ECG) et une conductance cutanée afin d'identifier d'éventuels signaux de stress face au tsunami. L'ECG est un enregistrement qui fournit différentes variables pertinentes pour la détection d'un état de stress [115]. Il rend compte des polarisations et dépolarisations successives du coeur traduites en pics caractéristiques à ce type d'enregistrement. En particulier, il fournit la fréquence cardiaque (FC) exprimée en battements par minute, calculée à partir de la fréquence d'apparition d'un type de pic pendant une minute d'enregistrement. La conductance cutanée mesure la variation continue des caractéristiques électriques de la peau. Les données sont collectées en appliquant une tension faible, indétectable et constante sur la peau, puis en mesurant la variation de la conductance cutanée.

Dans notre étude, seule la fréquence cardiaque a été retenue en raison de sa fiabilité [17]. Même si l'expérience a été menée en laboratoire, les légers mouvements des participants ont faussé certaines des mesures de conductance.

Avant de démarrer chaque simulation, les expérimentateurs ont réalisé un test pour vérifier le bon fonctionnement des appareils et pour familiariser les participants avec la RV. Ainsi, tous les participants ont vécu la même situation de test, qui était une immersion dans un village viking. Ce test nous a également permis d'enregistrer une ligne de base de la FC au repos spécifique à chaque individu. Tout au long de l'expérience de RV, l'état de stress de chaque participant a été déterminé par rapport à la ligne de base : un état de stress a été identifié dès que la FC du participant dépassait la ligne de base et donc augmentait.

À partir des mesures de FC, deux catégories de personnes ont été distinguées : les personnes calmes et les personnes stressées correspondant respectivement à une FC faible et à une FC élevée. Par conséquent, deux séries de mesures ont été obtenues.

La section suivante montre les liens entre les mesures physiologiques et les variables du modèle mathématique APC.

### 5.4.3 Le lien entre le modèle mathématique et les données expérimentales

Selon [151] et la section 5.3.1, trois types d'états comportementaux surviennent lors d'un événement catastrophique et dépendent de la charge émotionnelle et de la capacité des individus à la réguler : les états d'alerte,
de contrôle et de panique. Ces trois états ont conduits au modèle mathématique APC. Dans la section suivante, nous déterminons quelles variables ou combinaisons (de variables) du modèle APC peuvent être associées à l'une des mesures expérimentales.

### 5.4.3.1 Variables observées du modèle mathématique

Basés sur des connaissances en psychologie et en physiologie, les comportements d'alerte et de contrôle correspondent à un état de base dans lequel les personnes ont une fréquence cardiaque basse contrairement au comportement de panique (voir section 5.3.1 et [17]). Nous avons donc considéré que le premier ensemble de mesures correspondait aux comportements du quotidien, d'alerte et de contrôle dans le modèle APC et le second ensemble aux comportements de panique. Les personnes revenant à la vie quotidienne dans le modèle mathématique peuvent être identifiées aux personnes se mettant en sécurité dans le *le vieux Nice* ou terminant la simulation. Leur pourcentage pendant la durée de l'expérimentation correspond à la différence entre le nombre total de répondants faisant le scénario et le nombre encore dans la simulation.

Les données fournies par les psychologues sont rassemblées dans le tableau 5.4 et représentent les pourcentages de répondants à l'enquête présentant un comportement d'alerte-contrôle-quotidien, un comportement de panique et un retour au comportement quotidien dans les trois scénarios (neutre, organisé et non organisé). Une représentation graphique de la distribution de ces états comportementaux au cours du temps est fournie à la figure 5.11 pour le scénario neutre.

<figure><img src="image_50.png" /><figcaption>FIGURE 5.11 – Scénario neutre : représentation de la distribution (%) au cours du temps des comportements du quotidien-alerte-contrôle, des comportements de panique, et des comportements de retour au quotidien.</figcaption></figure>

L'étape suivante consistait à convertir les données du tableau 5.4 en densités pour le traitement mathématique. On note $S_m$ la densité de personnes à FC faible correspondant à la première série de mesures, $p_m$ la densité de personnes à FC élevée correspondant à la deuxième série et $b_m := 1 - S_m - p_m$ la densité de personnes ayant terminé le scénario de RV.

Les variables correspondantes du modèle mathématique (5.6) encore appelées sorties du modèle mathématique sont notées $y_1 = q + a + c$, $y_2 = p$ et $y_3 = b$.

### 5.4.3.2 Modélisation du début et de la fin de la catastrophe

À partir de la figure 5.10 et du Tableau 5.4, une représentation simple des fonctions $\gamma$ et $\varphi$ peut être trouvée. Dans chaque scénario, les personnes interrogées sont censées être engagées dans un comportement du quotidien jusqu'à 60 secondes, moment auquel le début de la vague peut être entendue. Ainsi, $\gamma$ qui modélise la catastrophe est supposée égale à 0 jusqu'à 60 secondes. Puis, à 90 secondes, lorsque la sirène retentit, on suppose que toutes les personnes prennent conscience de l'arrivée imminente du danger. Par conséquent, 100% des répondants à l'enquête ont adopté l'un des trois comportements : comportement d'alerte, de contrôle ou de panique. Entre 60 et 90 secondes, nous supposons que le pourcentage de personnes conscientes de la catastrophe
<table><caption>TABLE 5.4 – Tableau donnant le pourcentage de personnes ayant un comportement Alerté-Contrôlé-Quotidien (ACQ), un comportement de panique (P) et un Retour au comportement Quotidien (RQ) pour les trois scénarios : neutre, organisé et non organisé.</caption><thead><tr><th></th><th>Temps</th><th>40</th><th>60</th><th>80</th><th>100</th><th>120</th><th>140</th><th>160</th><th>180</th><th>200</th><th>220</th><th>240</th><th>260</th><th>280</th></tr></thead><tbody><tr><td rowspan="3">Neutre</td><td>P (%)</td><td>0.00</td><td>47.83</td><td>56.52</td><td>56.52</td><td>43.48</td><td>30.43</td><td>8.70</td><td>13.04</td><td>13.04</td><td>8.70</td><td>4.35</td><td>4.35</td><td>0.00</td></tr><tr><td>ACQ (%)</td><td>100.00</td><td>52.17</td><td>43.48</td><td>43.48</td><td>26.09</td><td>30.43</td><td>30.43</td><td>8.70</td><td>4.35</td><td>4.35</td><td>4.35</td><td>4.35</td><td>0.00</td></tr><tr><td>RQ (%)</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>30.43</td><td>39.13</td><td>60.87</td><td>78.26</td><td>82.61</td><td>86.96</td><td>91.30</td><td>91.30</td><td>100.00</td></tr><tr><td rowspan="3">Organisé</td><td>P (%)</td><td>0.00</td><td>44.00</td><td>52.00</td><td>56.00</td><td>52.00</td><td>36.00</td><td>12.00</td><td>8.00</td><td>8.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>ACQ (%)</td><td>100.00</td><td>56.00</td><td>48.00</td><td>44.00</td><td>44.00</td><td>28.00</td><td>20.00</td><td>4.00</td><td>4.00</td><td>4.00</td><td>4.00</td><td>4.00</td><td>0.00</td></tr><tr><td>RQ (%)</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>4.00</td><td>36.00</td><td>68.00</td><td>88.00</td><td>88.00</td><td>96.00</td><td>96.00</td><td>96.00</td><td>100.00</td></tr><tr><td rowspan="3">Non Organisé</td><td>P (%)</td><td>0.00</td><td>53.85</td><td>61.54</td><td>50.00</td><td>65.38</td><td>30.77</td><td>11.54</td><td>7.69</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>ACQ (%)</td><td>100.00</td><td>46.15</td><td>38.46</td><td>50.00</td><td>26.92</td><td>26.92</td><td>30.77</td><td>7.69</td><td>3.85</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>RQ (%)</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>7.69</td><td>42.31</td><td>57.69</td><td>84.62</td><td>96.15</td><td>100.00</td><td>100.00</td><td>100.00</td><td>100.00</td></tr></tbody></table>
augmente régulièrement. La fonction :

$$ \phi(s) = \begin{cases} 0 & \text{if } s < \tau_0 \\ 1 & \text{if } s > \tau_1 \\ \frac{1}{2} - \frac{1}{2} \cos\left(\frac{s-a}{b-a}\pi\right) & \text{if } \tau_0 \le s \le \tau_1 \end{cases} \quad (5.12) $$

est reprise avec $\tau_0 = 60$, $\tau_1 = 90$ pour obtenir une approximation de $\gamma$ (voir figure 5.12).

<figure><img src="image_51.png" /><figcaption>FIGURE 5.12 – Représentation de la fonction $\gamma$ qui modélise la catastrophe. Pour cela, la fonction (5.12) est reprise avec $\tau_0 = 60$, $\tau_1 = 90$.</figcaption></figure>

De la même manière, la fonction $\varphi$ décrivant le retour à un comportement quotidien est construite à partir de $\phi$ avec $\tau_0$ (resp. $\tau_1$) le temps à partir duquel le premier enquêté (resp. le dernier) quitte la simulation. $\varphi$ doit être calculé pour chaque scénario puisque les personnes n'adoptent pas les mêmes comportements. Par exemple, le scénario organisé se termine à 280 secondes contrairement au scénario non organisé qui se termine à 200 secondes. Cette différence montre que les processus d'imitation étaient différents selon les scénarios.

## 5.5 Sensibilité des sorties, identifiabilité et nombre minimal de données d'observation

Cette section décrit la stratégie mise en oeuvre pour estimer les paramètres à partir de l'ensemble des mesures fournies dans le tableau 5.4.

### 5.5.1 Étude de sensibilité

Une étude de sensibilité permet de quantifier l'influence des valeurs de paramètres sur la variabilité de la réponse du modèle. Il est ainsi possible de savoir quels sont les paramètres qui vont contribuer le plus à sa variabilité, les paramètres les moins influents et ceux qui interagissent entre eux.

En étudiant comment la réponse du modèle réagit aux variations de ses paramètres, l'analyse de sensibilité permet de répondre à un certain nombre de questions :

1. Le modèle est-il bien fidèle au phénomène/processus modélisé ? En effet, si l'analyse exhibe une forte influence d'un paramètre habituellement connu comme non influent, il sera nécessaire de remettre en cause la qualité du modèle ou (et) la véracité de nos connaissances sur l'impact réel des variables d'entrée.

2. Quelles sont les paramètres qui contribuent le plus à la variabilité de la réponse du modèle ? Si cette variabilité est synonyme d'imprécision sur la valeur prédite de la sortie, il sera alors possible d'améliorer la qualité de la réponse du modèle en restreignant le domaine d'appartenance du paramètre par exemple. Toutefois, cela n'est pas toujours possible, notamment lorsque la variabilité d'un paramètre est intrinsèque à la nature du paramètre et non due à un manque d'information ou à des imprécisions de mesures.

3. Quelles sont au contraire les paramètres les moins influents ? Il sera possible de les considérer comme des paramètres déterministes, en les fixant par exemple à leur espérance, et obtenir ainsi un modèle plus simple avec moins de paramètres.
4. Quels paramètres ou quels groupes de paramètres interagissent avec quels autres ? L'analyse de sensibilité peut permettre de mieux appréhender et comprendre le phénomène modélisé, en éclairant les relations entre les paramètres.

Pour mener ce type d'étude, il existe deux approches : les méthodes locales et les méthodes globales [179].
Alors que pour l'analyse de sensibilité locale, un seul paramètre varie lorsque les autres paramètres sont fixés,
la sensibilité globale vise à étudier la réponse du modèle lorsque tous les paramètres varient.

Les auteurs de [179] distinguent trois classes de techniques d'analyse de sensibilité globale : les méthodes de régression, les méthodes de filtrage et celles basées sur la variance. J'ai opté pour la dernière classe car elle fournit des résultats quantitatifs en considérant tout l'espace des paramètres possibles (voir l'Annexe 2). Ceci est important pour notre modèle car nous aimerions considérer différents scénarios conduisant à des valeurs très différentes des paramètres.

J'ai repris le package Sensitivity [98] implémenté dans le logiciel R dédié au calcul statistique.

Les résultats ont montré une très grande sensibilité de a, p et c par rapport aux paramètres intrinsèques comparés aux paramètres d'imitation (partie 2 de l'Annexe 2). Cette étude nous conduit à dissocier l'estimation de ces deux groupes de paramètres et favoriser d'abord l'estimation des paramètres intrinsèques qui pourra être faite à partir du premier scénario. En supposant que le profil des participants est le même dans les trois scénarios, les paramètres d'imitation pourront alors être calculés à partir des scénarios organisé et non organisé (s'ils sont identifiables).

5.5.2 Identifiabilité

Pour tester l'identifiabilité du modèle APC, l'algorithme *IdentifiabilityTree* développé dans [204] et présenté dans le premier chapitre a été utilisé.

Soient les hypothèses suivantes $\mathcal{H}$ :

1. Les fonctions $\gamma$ et $\varphi$ sont connues

2. $y_1(t) = q(t) + a(t) + c(t)$, $y_2(t) = p(t)$, $y_3(t) = b(t)$.

Identifiabilité du modèle APC sans imitation

Le premier scénario où il n'y a pas d'imitation implique pour le modèle mathématique que $F = G = H = 0$ (et donc que $\alpha = \beta = \gamma_1 = \gamma_2 = 0$) et que le vecteur paramètre à identifier est $\tilde{\theta} := (B_1, B_2, B_3, B_4, C_1, C_2)$.

La proposition suivante donne le résultat d'identifiabilité lorsque $\alpha = \beta = \gamma_1 = \gamma_2 = 0$.

**Proposition 25.** S'il n'y a pas d'imitation, sous les hypothèses $\mathcal{H}$, le modèle (5.6) est identifiable en $\tilde{\theta}$.

La preuve est basée sur les polynômes ES (5.13) qui seront utilisés dans la partie 5.5.3 lorsqu'on déterminera
le nombre minimal de mesures pour estimer les paramètres.

$$
\begin{align*}
P_1 &:= (\ddot{y}_2 + \dot{y}_2 \gamma) \varphi + \left(-B_1 C_2 + B_1 C_2^2 / B_2 - B_2 B_3 - B_2 C_2 + B_3 C_2 + C_2^2\right) \dot{y}_3 \\
&\quad + C_2 \dot{y}_3 \varphi + (B_2 - C_2) \dot{y}_3 \gamma + (B_1 - B_1 C_2 / B_2 + B_2 + B_4 + C_1) \dot{y}_2 \varphi \\
&\quad + (B_1 B_4 + B_1 C_1 - B_1 B_4 C_2 / B_2 - B_1 C_1 C_2 / B_2 + B_2 C_1 - C_1 C_2) y_2 \varphi + (B_4 + C_1) y_2 \gamma \varphi - B_2 y_1 \gamma \varphi = 0,
\end{align*}
$$

$$
P_2 := -(\dot{y}_3 \varphi + \ddot{y}_3 \varphi + \dot{y}_3 \varphi^2) + (B_1 C_2 / B_2 + B_3 + C_2) \dot{y}_3 \varphi - B_1 / B_2 \dot{y}_2 \varphi^2 + (-B_1 B_4 / B_2 - B_1 C_1 / B_2 - C_1) y_2 \varphi^2 = 0. \tag{5.13}
$$

Identifiabilité des paramètres d'imitation

Maintenant que les paramètres de la partie linéaire sont identifiables, une étude d'identifiabilité peut être
réalisée sur l'ensemble des paramètres d'imitation, c'est à dire $\tilde{\theta} := (\alpha, \beta, \gamma_{c \to p}, \gamma_{p \to c})$. Nous avons la proposition
suivante :

**Proposition 26.** On suppose que les paramètres intrinsèques sont connus. Sous l'hypothèse $\mathcal{H}$, le modèle (5.6) est identifiable par rapport au $\tilde{\theta}$.

Étant donné les non linéarités du modèle, il n'a pas été possible de déterminer les polynômes ES. La preuve a été faite en travaillant directement avec les équations du modèle.

**Remarque 30.** On peut montrer en utilisant le principe de la démonstration de la proposition 26 que le modèle (5.6) est identifiable.
Comme nous le verrons dans la suite, il n'est pas possible d'estimer tous les paramètres à partir du système complet et des données expérimentales. En effet, nous ne disposons que d'un faible nombre de mesures. Pour le justifier, la partie suivante présente un résultat sur le calcul du nombre minimal de données pour avoir l'identifiabilité locale et donc pour pouvoir estimer les paramètres de façon unique.

### 5.5.3 Nombre minimal de données d'observation

Pour déterminer le nombre minimal de mesures nécessaires pour avoir l'identifiabilité locale, nous avons adapté une méthode présentée dans [141, 220]. Le résultat est basé sur une fonction $\rho = \rho(\Theta, y, \dot{y}, \dots, y^{(k)})$, où $k$ est un entier positif tel que $\rho$ a des dérivées partielles continues par rapport à la variable $\Theta$. L'identifiabilité locale vient alors de l'inversion locale de $\rho$ en utilisant le théorème des fonctions implicites dans un voisinage où la fonction s'annule. [141] a noté que certains polynômes différentiels obtenus par des méthodes d'algèbre différentielle peuvent générer cette fonction $\rho$. Ce travail est repris et généralisé aux polynômes ES définis en (2.4). Pour simplifier les notations, on suppose qu'il n'y a qu'une seule sortie et donc un seul polynôme ES et on note $\Theta = (\gamma_k(\theta))_{1\le k\le d}$ la suite de coefficients du polynôme ES.

**Proposition 27.** *Supposons que*

1. la fonction

$$ \rho: \mathbb{R}^d \times \mathbb{R}^m \to \mathbb{R}^d \\ (\Theta, y) \mapsto (P(\theta, y), P(\theta, y)^{(1)}, \dots, P(\theta, y)^{(d-1)}) \qquad (5.14) $$

est une fonction différentiable en la première variable

2. $|\frac{\partial \rho}{\partial \Theta}|$ n'est pas identiquement nulle

3. la fonction $\phi$ définie en (2.5) est injective.

Le modèle est localement identifiable si le nombre minimal de mesures de y est $s_y+d$, où $s_y$ est la plus haute dérivée de y dans P.

**Remarque 31.** ✓ $|\frac{\partial \rho}{\partial \Theta}|$ n'est autre que le wronskien $W(m_1, \dots, m_d)$ défini à la proposition 1. La deuxième hypothèse permet d'exprimer $\Theta$ en fonction de y et de ses dérivées.

✓ L'injectivité de la fonction $\phi$ permet de retrouver les paramètres $\theta_i$ à partir des $\Theta_j$.

**Exemple 5.5.1.** Des polynômes $P_1$ et $P_2$ donnés en (5.13), on peut en déduire le nombre minimal d'observations pour avoir l'identifiabilité locale.

$P_1$ est composé de 7 coefficients et les indéterminées sont $y_1, y_2, y_3$. Son leader étant $\ddot{y}_2$, $7+2=9$ points sont nécessaires.

Pour le polynômes $P_2$, il faut $3+2=5$ mesures.

En conclusion, le nombre minimal de mesures requis pour avoir l'identifiabilité locale est de 9.

**Remarque 32.** — À partir du tableau 5.4, il est possible de mettre en place une procédure d'estimation de paramètres pour estimer les paramètres intrinsèques.

— Estimer tous les paramètres à partir des scénarios 2 et 3 ne donnera pas de bons résultats. En effet, l'ordre de dérivation des polynômes obtenus à partir du modèle complet sera plus élevé que celui de $P_1$ et on dispose de moins de points de mesure pour ces deux scénarios.

En résumé, l'identification des paramètres a été effectuée en trois étapes correspondant aux trois conditions de l'expérience de RV. L'hypothèse principale est que les participants ont le même profil dans les trois scénarios. De chacun d'entre eux, on a estimé un ensemble de paramètres. Ainsi, les paramètres intrinsèques ont d'abord été estimés à partir du scénario neutre dans lequel il n'y a pas de processus d'imitation puis les paramètres décrivant les processus d'imitation de contrôle à partir du scénario organisé et enfin les paramètres décrivant les processus d'imitation de panique à partir du scénario non organisé. La Figure 5.13 résume les étapes de l'identification.

### 5.5.4 Estimation des paramètres

Soient $t_i$, $i=1, \dots, M$ les temps de mesures des données physiologiques, et les vecteurs $S_m = (S_m(t_i))_{1,...,M}$, $p_m = (p_m(t_i))_{1,...,M}$ et $b_m = (b_m(t_i))_{1,...,M}$ regroupant les densités de personnes dans un état contrôlé, paniqué et dans un retour à un état du quotidien à chaque pas de temps. Soient les matrices $y^{\text{mes}} = (S_m, p_m, b_m)^T$ et $y(\theta) = (y_1(t_i), y_2(t_i), y_3(t_i))_{i=1,...,M}^T$ obtenus à partir des sorties du modèle évaluées à chaque $t_i$ et en le vecteur de paramètre $\theta$. $y^{\text{mes}} - y(\theta)$ est une matrice de dimension $M \times 3$ regroupant toutes les différences entre les mesures et les sorties du modèle.
<figure><img src="image_52.png" /><figcaption>FIGURE 5.13 – Diagramme résumant les étapes pour estimer les paramètres du modèle APC. Les transitions naturelles sont tout d'abord estimées à partir du premier scénario, puis viennent les processus d'imitation à partir des deuxième et troisième scénarios</figcaption></figure>

La fonction coût $J(\theta)$ est alors

$$J(\theta) = \|y^{mes} - y(\theta)\|_{2}^{2} = \sum_{i=1}^{M} \sum_{j=1}^{3} (y_{i,j}^{mes} - y(\theta)_{i,j})^{2}$$

où $y_{i,j}^{mes}$ (resp. $y(\theta)_{i,j}$) est la composante de la i-ème ligne et la j-ème colonne de la matrice $y^{mes}$ (resp $y(\theta)$). La solution de l'algorithme d'optimisation est notée

$$\hat{\theta} = \operatorname*{argmin}_{\theta \in \mathcal{P}} J(\theta) \qquad (5.15)$$

Le problème (5.15) a été résolu avec l'algorithme génétique implémentée dans Scilab avec la fonction `optim_ga` [82]. En effet, un algorithme global a été préféré puisque nous ne disposons pas de valeurs initiales pour utiliser un algorithme local. Dans [201], une étude sur l'erreur relative a été faite.

### 5.5.5 Résultats

#### Scénario sans avatar

Dans ce scénario, 29 personnes ont été enquêtées mais seulement 23 ont été retenues en raison de la qualité de leurs mesures. Dans le problème d'optimisation (5.15), les variables et leurs contraintes sont $\theta = (B_1, B_2, B_3, B_4, C_1, C_2) \in \mathcal{P}$ où

$$\mathcal{P} = \{\theta \in [0.01, 0.7]^2 \times [0.01, 0.2]^2 \times [0.01, 0.7]^2 / B_2 + C_2 < 1\}.$$

La condition $B_2 + C_2 < 1$ a été ajoutée pour que Runge-Kutta converge. En effet, le modèle étant très sensible aux paramètres intrinsèques, de grandes valeurs rendent Runge-Kutta divergent.

L'algorithme génétique a retourné les valeurs suivantes pour les paramètres intrinsèques : $B_1 = 0.1915$, $B_2 = 0.4374$, $B_3 = 0.1009$, $B_4 = 0.1741$, $C_1 = 0.1448$, $C_2 = 0.4676$.

La figure de gauche 5.14 représente les trajectoires des données réelles $S_m, p_m, b_m$ et les trajectoires des sorties calculées avec les paramètres estimés. La figure de droite représente les simulations de $a, p, c, b, q$.

**Remarque 33.** Nous pouvons noter que la population a une faible culture du risque. Les inégalités $B_2 > B_1$ et $C_2 > C_1$ impliquent que les individus ont tendance à adopter un comportement paniqué de leur propre initiative.
<figure><img src="image_53.png" /><figcaption>FIGURE 5.14 – Figure de gauche : Représentation des données réelles $S_m, p_m, b_m$ et des trajectoires de sorties : densités de comportements contrôlés ($y_1$), paniqués ($y_2$) et de retour au quotidien après la catastrophe ($y_3$). Figure de droite : Simulations de $a, p, c, b, q$ avec les paramètres obtenus à partir de l'algorithme génétique.</figcaption></figure>

## Scénario organisé

Dans le deuxième scénario, 25 personnes ont été sélectionnées parmi 29. Les valeurs des paramètres intrinsèques trouvées dans le premier scénario ont été utilisées. Dans ce scénario, puisqu'il n'y a pas d'emojis en état de panique, nous avons supposé que $\beta = \gamma_{c \to p} = 0$.

Le problème d'optimisation est résolu avec $\theta = (\alpha, \gamma_{p \to c})$ et $P = [0, 9]^2$. Les résultats obtenus avec l'algorithme génétique sont $\hat{\theta} = (\alpha, \gamma_{p \to c}) = (0.8338, 0)$. La figure 5.15 donne la représentation des données réelles avec les trajectoires de sortie évaluées en $\hat{\theta}$ puis les simulations des variables d'état du modèle mathématique.

<figure><img src="image_54.png" /><figcaption>FIGURE 5.15 – Scénario organisé. Figure de gauche : Représentation des données réelles $S_m, p_m, b_m$ et des trajectoires de sorties : densités de comportements contrôlés ($y_1$), paniqués ($y_2$) et de retour au quotidien après la catastrophe ($y_3$). Figure de droite : Simulations de $a, p, c, b, q$ avec les paramètres obtenus à partir de l'algorithme génétique.</figcaption></figure>

## Scénario non-organisé

Pour le troisième scénario, 26 personnes ont été retenues. Comme dans le scénario précédent, les valeurs des paramètres intrinsèques trouvées dans le premier cas ont été utilisées et comme il n'y avait pas d'emojis en état de contrôle, nous avons supposé que $\alpha = \gamma_{p \to c} = 0$. Dans le problème d'optimisation (5.15), $\theta = (\beta, \gamma_{c \to p})$ et $P = \{p \in [0, 1]^2 | \beta + \gamma_{c \to p} < 0.8\}$. Nous avons contraint les valeurs des paramètres pour limiter les valeurs de $p$. Le résultat obtenu avec l'algorithme génétique est $\hat{\theta} = (0.6599, 0.0273)$.

Dans les trois cas, la dynamique des différents comportements a été retrouvée.

## 5.6 Un modèle mathématique, pour quoi faire ?

L'exploitation du modèle mathématique nous a permis

- de définir de façon plus précise les comportements, leurs enchaînements et les déplacement associés des individus à partir les entretiens immersifs (voir Figure 5.7);
<figure><img src="image_55.png" /><figcaption>FIGURE 5.16 – Scénario non organisé. Figure de gauche : Représentation des données réelles $S_m, p_m, b_m$ et des trajectoires de sorties : densités de comportements contrôlés ($y_1$), paniqués ($y_2$) et de retour au quotidien après la catastrophe ($y_3$). Figure de droite : Simulations de $a, p, c, b, q$ avec les paramètres obtenus à partir de l'algorithme génétique.</figcaption></figure>

— de donner un ordre de grandeur des paramètres du modèle à partir de la réalité virtuelle dont le scénario était basé sur des situations simulées [201];

— de définir et étudier des moyens pour limiter les comportements de panique [142];

— d'adopter une approche réseau complexe pour modéliser l'influence de l'espace dans les transitions d'un comportement à un autre dans le cas d'un tsunami à Nice [120, 119];

— d'étudier l'importance de la capacité des zones de refuge sur les comportements de la population [120].

Dans la suite, je présente deux exemples : un exemple montrant que le travail de calibration permet une meilleure compréhension des dynamiques non observables [201]. Le deuxième concerne un modèle contrôlé visant à trouver la meilleure stratégie pour limiter les comportements de panique [142].

Une meilleure compréhension des dynamiques non observables

À partir du modèle mathématique et des résultats sur les paramètres obtenus à la section 5.5.5, une étude plus fine des processus de transition dans les changements de comportement a été faite lors de l'expérimentation de réalité virtuelle.

En reprenant les fonctions $Im_p, Im_n, Im_c$ et $Im_c$ définies en (5.8) et (5.10), nous avons obtenu les figures 5.17 pour le scénario organisé et 5.18 pour le scénario non organisé. La première constatation est que les processus intrinsèques dominent les processus d'imitation dans les deux scénarios. Ensuite, les processus d'imitation sont plus importants dans le scénario non organisé ce qui tend à montrer que le comportement de panique est plus fréquemment imité.

<figure><img src="image_56.png" /><figcaption>FIGURE 5.17 – Flux dus aux processus intrinsèque et aux processus d'imitation vers trinsèque et aux processus d'imitation vers le contrôle dans le scénario organisé par rapport au temps.
FIGURE 5.18 – Flux dus aux processus intrinsèque et aux processus d'imitation vers le contrôle dans le scénario non-organisé par rapport au temps.</figcaption></figure>

<caption>FIGURE 5.19 – Flux dus aux processus intrinsèques et aux processus d'imitation.</caption>
En mathématiques, une fonction de contrôle est une fonction qui va agir sur le système pour modifier son comportement dynamique. Dans notre exemple, les fonctions de contrôle doivent agir sur la panique en limitant ou en diminuant le nombre de paniqués dans la foule et en réduisant le nombre de décès. Dans notre exemple, une fonction de contrôle doit correspondre à un effet provenant d'une action menée par les opérationnels et les institutionnels, que ce soit avant (p.e. la prévention) ou pendant l'événement (p.e. la gestion de catastrophe). Ainsi, les variables de contrôle sont liées à :

- La présence d'opérationnels sur le lieu de la catastrophe

- La formation/l'expérience des populations avant une catastrophe

- Les consignes transmises pendant l'événement à la population. Ces consignes favorisent les comportements contrôlés comme l'évacuation, le confinement ou d'autres comportements contrôlés.

On va supposer par la suite que toute l'organisation des actions opérationnelles, des services, fonctionne bien, même si ce n'est pas toujours le cas. Trois types de contrôle ont été envisagés :

1. $u_1$: ce contrôle prend en compte l'effet de la formation et de l'expérience de la population. Ils sont acquis avant la catastrophe. Cet effet peut varier au cours du temps. Il agit sur la transition intrinsèque entre le comportement en alerte et le comportement contrôlé, soit sur le paramètre $B_1$. On aura donc $B_1+u_1$ au lieu de $B_1$ dans le système (5.6). Comme $B_1+u_1$ est un taux de transition, il est compris entre 0 et 1.

2. $u_2$: il représente les actions des opérationnels qui sont présents et qui interviennent sur le terrain dès le début de l'événement. Ils ont un effet apaisant tout au long de l'événement et on le suppose constant. $u_2$ est une constante qui agit sur deux transitions de comportements :

- le passage de comportement d'Alerte à Contrôle

- le passage de comportement de Panique à Contrôle

On suppose que cette constante agit de la même façon (c'est à dire avec la même valeur lors de la simulation) sur les deux transitions. $u_2$ agit sur la fonction d'imitation $F$ et en particulier sur le paramètre $\alpha$ pour forcer le passage de $A$ à $C$. On aura donc $\alpha+u_2$ dans le système contrôlé.

$u_2$ agit également sur la fonction d'imitation $H$ et en particulier sur le paramètre $\gamma_2$ pour forcer les passages de $P$ à $C$. On aura donc $\gamma_2+u_2$ dans le système contrôlé.

3. $u_3$: ce contrôle traduit l'effet de l'information institutionnelle transmise pendant l'événement, comme la transmission d'un message ou le déclenchement d'une sirène. Cette fonction est une fonction adaptative qui dépend du temps. En particulier, c'est une fonction constante par morceaux. Par exemple, à la figure 5.20, on a supposé qu'une sirène a été déclenchée au temps $t_1$ et a été suivie par un message des opérationnels au temps $t_2$. Chaque saut correspond donc à une information supplémentaire apportée à la population et on suppose que les effets s'accumulent.

Les personnes n'agiront pas de la même façon à l'information transmise pendant l'événement $u_3$ selon leur expérience. Par exemple, les personnes ayant suivi une formation sur la façon de réagir à l'arrivée d'un tsunami, adopteront des réactions adaptées lorsqu'elles entendront le bruit spécifique de la sirène pour ce type d'événement. Ainsi, le contrôle $u_3$ agit sur la transition intrinsèque entre Alerte-Contrôle et on aura $B_1+u_1+\rho u_3 \in [0,1]$ où $\rho \in [0,1]$. Ensuite, $u_3$ favorise l'imitation dans le sens Alerte-Contrôlé (donc agit sur $\alpha$) et Paniqué-Contrôlé (agit sur $\gamma_2$).

Les fonctions de contrôles sont prises constantes par morceau, voir par exemple une représentation de $u_3$ à la figure 5.20. Le système contrôlé s'écrit alors :

$$
\left\{
\begin{aligned}
\frac{da(t)}{dt} &= \gamma(t)q - (B_1 + u_1 + \rho u_3 + B_2 + D_a)a(t) + B_3 c(t) + B_4 p(t) - \bar{F}(a(t), c(t), u_2, u_3)a(t)c(t) - G(a(t), p(t))a(t)p(t), \\
\frac{dc(t)}{dt} &= (B_1 + \rho u_3)a(t) + C_1 p(t) - (B_3 + C_2 + D_c)c(t) + \bar{F}(a(t), c(t), u_2, u_3)a(t)c(t) \\
&\quad + \bar{H}(c(t), p(t), u_2, u_3)c(t)p(t) - \varphi(t)c(t), \\
\frac{dp(t)}{dt} &= B_2 a(t) + C_2 c(t) - (B_4 + C_1 + D_p)p(t) + G(a(t), p(t))a(t)p(t) - \bar{H}(c(t), p(t), u_2, u_3)c(t)p(t), \\
\frac{dq(t)}{dt} &= -\gamma(t)q(t), \quad \frac{db(t)}{dt} = \varphi(t)c(t)
\end{aligned}
\right.
\qquad (5.16)
$$
<figure><img src="image_57.png" /><figcaption>FIGURE 5.20 – Exemple de représentation de la fonction de contrôle $u_3$.</figcaption></figure>

avec les fonctions d'imitation :

$$ egin{cases} \bar{F}(a, c) = \frac{\alpha + u_2 + u_3}{N} \eta \left( \frac{c}{a + \varepsilon} \right) \\ \\ \bar{H}(c, p) = - \frac{\gamma_1}{N} \eta \left( \frac{p}{c + \varepsilon} \right) + \frac{\gamma_2 + u_2 + u_3}{N} \eta \left( \frac{c}{p + \varepsilon} \right), \end{cases} \quad (5.17) $$

L'existence et l'unicité de la solution du système contrôlé est classique pour tout contrôle

$$ u \in \mathcal{U}_T = \{(u_1, u_2, u_3) | \ u_i \text{ fonction constante par morceaux, } a_i \le u_i(t) \le b_i, i = 1, 2, 3, \forall t \in [0, T]\}, $$

$a_i, b_i$ étant des constantes comprises entre 0 et 1, et tout vecteur de paramètres $\theta \in \Theta := (\mathbb{R}_+)^6 \times (\mathbb{R}_+)^7$, ceci quelque soit la condition initiale positive ou nulle.

Considérons la fonction objectif suivante :

$$ J(u_1, u_2, u_3) = \int_0^T (C_p p(t) + C_d d(t) + A_1 u_1^2(t) + A_2 u_2^2(t) + A_3 u_3^2(t)) dt. \quad (5.18) $$

Les deux premiers termes représentent le bénéfice des comportements de panique et de décès. Les constantes $A_1, A_2$ et $A_3$ sont positives et correspondent aux poids qui permettent de réguler les coûts du contrôle pour la formation, les actions des secouristes et les informations externes, respectivement. Notre objectif est donc de minimiser cette fonction coût pour trouver le ou les contrôles minimisant le phénomène de panique.

Le principe de Pontryagin [188] permet de conclure à l'existence d'un contrôle optimal $u^* = (u_1^*, u_2^*, u_3^*) \in \mathcal{U}_T$ répondant à la question c'est à dire tel que

$$ \min_{(u_1, u_2, u_3) \in \mathcal{U}_T} J(u_1, u_2, u_3) = J(u_1^*, u_2^*, u_3^*). \quad (5.19) $$

Nous avons cherché à répondre à la question suivante :

**Pour un scénario donné, quelle est la meilleure stratégie pour réduire la panique et limiter le nombre de victimes ?**

À chaque contrôle $u_1, u_2$ et $u_3$ correspond un coût noté $A_1, A_2$, et $A_3$, respectivement. Pour chaque scénario, une stratégie correspond à une combinaison de ces coûts (ex : $A_1 = 0.01, A_2 = 0.01, A_3 = 1$). Pour modéliser les différentes stratégies, le système binaire suivant est adopté : $A_i = 0.01$ si la stratégie vise surtout à augmenter les efforts du contrôle $u_i$, et $A_i = 1$ si la stratégie vise à diminuer les efforts du contrôle $u_i$ ($i = 1, 2, 3$).

On rappelle que $u_1$ est le contrôle associé à la formation de la population, $u_2$ celui associé aux actions des opérationnels sur les lieux du sinistre, et $u_3$ le contrôle associé à la communication. Ainsi, une stratégie de gestion correspondra au triplet suivant ($A_1 = 1, A_2 = 0.01, A_3 = 1$), une stratégie de prévention basée sur
l'éducation de la population correspondra au triplet ($A_1 = 0.01$, $A_2 = 1$, $A_3 = 1$) et une stratégie de prévention basée sur la communication correspondra au triplet ($A_1 = 1$, $A_2 = 1$, $A_3 = 0.01$).
Nous avons testé les stratégies définies dans le tableau 5.5.

<table><caption>TABLE 5.5 – Les différentes stratégies imaginées</caption><thead><tr><th>Stratégies</th><th>$A_1$</th><th>$A_2$</th><th>$A_3$</th></tr></thead><tbody><tr><td>1</td><td>$10^{-2}$</td><td>$10^{-2}$</td><td>$10^{-2}$</td></tr><tr><td>2</td><td>$10^{-2}$</td><td>$10^{-2}$</td><td>1</td></tr><tr><td>3</td><td>$10^{-2}$</td><td>1</td><td>$10^{-2}$</td></tr><tr><td>4</td><td>$10^{-2}$</td><td>1</td><td>1</td></tr><tr><td>5</td><td>1</td><td>$10^{-2}$</td><td>$10^{-2}$</td></tr><tr><td>6</td><td>1</td><td>$10^{-2}$</td><td>1</td></tr><tr><td>7</td><td>1</td><td>1</td><td>$10^{-2}$</td></tr><tr><td>8</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

La stratégie 1 est une stratégie de référence qui permet d'observer les contrôles lorsqu'ils ont le même coût. Elle ne sera donc pas considérée pour déterminer la meilleure stratégie.

Afin de déterminer quelle est la meilleure stratégie qui réduit la panique, nous examinons celle dont les contrôles nous donnent le moins de personnes en état de panique et telle que la valeur de la fonction objectif est la plus petite. Nous notons :

$$I_{\Delta p} = \int_{0}^{T} (p(t) - p_{c}(t))dt \qquad (5.20)$$

où $p_c$ et $p$ correspondent aux densités d'individus dans un état de panique avec et sans contrôle. $I_{\Delta p}$ va permettre de juger l'efficacité des contrôles : plus la valeur de $I_{\Delta p}$ est grande, plus les contrôles réduisent la panique.

De la même manière, la quantité qui fournit une mesure similaire de l'efficacité des contrôles par rapport aux décès est définie par :

$$I_{\Delta d} = \int_{0}^{T} (d(t) - d_{c}(t))dt \qquad (5.21)$$

Par conséquent, pour chaque scénario, la meilleure stratégie pour réduire la panique correspondra à la stratégie ayant l'indice $I_{\Delta p}$ le plus élevé et la valeur la plus faible de la fonction de coût $J$ définie en (5.18).

De la même façon, la meilleure stratégie pour réduire les décès correspondra à la stratégie pour laquelle l'indice $I_{\Delta d}$ sera le plus élevé et la fonction coût $J$ la plus basse.

Dans toutes nos simulations, nous avons considéré $C_p = C_d = 1$.

### Scénario 1 : foule dense avec une faible culture du risque

Les paramètres sont : $B_1 = 0.05$, $B_2 = 0.051$, $B_3 = B_4 = 10^{-3}$, $C_1 = 0.05$, $C_2 = 0.01$, $D_r = D_c = 10^{-4}$, $D_p = 1.5 \times 10^{-5}$, $\alpha_1 = \beta_1 = \gamma_1 = \gamma_2 = 0.7$, $\rho = 0.5$ $^4$.

Dans la figure 5.21, afin d'avoir tous les indices d'un même graphique, nous avons dû multiplier $I_{\Delta d}$ par cent.

En testant toutes les stratégies du tableau 5.5, nous constatons que la meilleure stratégie pour réduire la panique et la mort (voir Fig. 5.21) correspond à la stratégie 2 : $A_1 = 10^{-2}$, $A_2 = 10^{-2}$, $A_3 = 1$. Cependant, si l'on préfère une stratégie dont les efforts ne portent que sur un seul contrôle (stratégies 4, 6 et 7), il serait intéressant de choisir la stratégie 6.

### Scénario 2 : foule éparpillée avec un faible risque de la culture

On rappelle que les paramètres sont : $B_1 = 0.1$, $B_2 = 0.11$, $B_3 = 10^{-3}$, $B_4 = 0.001$, $C_1 = 0.1$, $C_2 = 0.2$, $D_r = D_c = 10^{-4}$, $D_p = 1.5 \times 10^{-5}$, $\alpha_1 = \beta_1 = \gamma_1 = \gamma_2 = 0.2$, $\rho = 0.5$.

Les simulations numériques de la série de stratégies du tableau 5.5 et leurs résultats représentés sur la figure 5.25 montrent que la meilleure stratégie pour réduire la panique et les décès lors de la catastrophe serait la 2 :

<footnote n="4">4. L'échelle de temps n'est pas la même qu'à la section 5.4.2 ce qui explique la différence de valeur des paramètres.</footnote>
<figure><img src="image_58.png" /><figcaption>FIGURE 5.21 – 100 × $I_{\Delta d}$ (blue), $I_{\Delta p}$ (green) et $J$ (red) pour les 8 stratégies de contrôle dans le cas d'une foule dense avec une faible culture du risque.</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>FIGURE 5.22 – Contrôles optimaux.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>FIGURE 5.23 – Courbe des densités du nombre de paniqués et du nombre de morts.</figcaption></figure>

FIGURE 5.24 – Solutions pour la meilleure stratégie dans le cas d'une foule dense avec une faible culture du risque : ($A_1 = 10^{-2}$, $A_2 = 10^{-2}$, $A_3 = 1$). (a) Gauche : Contrôles optimaux $u_1$, $u_2$ et $u_3$. (b) Droite : En haut : Densité de comportements paniqués avec (trait plein) et sans contrôle (trait en pointillé). En bas : Densité de morts avec (trait plein) et sans contrôle (trait en pointillé).
$(A_1 = 10^{-2}, A_2 = 10^{-2}, A_3 = 1)$. Il s'agit d'une stratégie mixte alliant prise en charge et prévention basée sur la formation de la population. La meilleure stratégie impliquant un seul contrôle est la stratégie 6.

<figure><img src="image_61.png" /><figcaption>FIGURE 5.25 – $100 \times I_{\Delta d}$ (blue), $I_{\Delta p}$ (green) et $J$ (red) pour les 8 stratégies de contrôle dans le cas d'une population éparpillée avec un faible risque de culture.</figcaption></figure>

À la figure 5.28 sont représentés les contrôles optimaux et l'évolution en temps des densités de paniqués et de morts, avec ou sans contrôle.

<figure><img src="image_62.png" /><figcaption>FIGURE 5.26 – Contrôles optimaux.</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>FIGURE 5.27 – Courbe des densités du nombre de paniqués et du nombre de morts.</figcaption></figure>

FIGURE 5.28 – Solutions pour la meilleure stratégie dans le cas d'une foule éparpillée avec une faible culture du risque ($A_1 = 10^{-2}, A_2 = 10^{-2}, A_3 = 1$). À gauche sont représentés les contrôles optimaux $u_1, u_2$ et $u_3$. En haut à droite est représenté la densité de comportements paniqués avec (trait plein) et sans contrôle (trait en pointillé); en bas, la densité de morts avec (trait plein) et sans contrôle (trait en pointillé).

Ainsi, dans le cadre de la théorie (statistique) de la décision, il serait possible d'ouvrir la voie à un outil logiciel d'aide à la décision générant des scénarios probabilistes à partir d'un espace (éventuellement continu) d'actions possibles, du type de population et de fonctions mérites multi-objectifs. De cette manière, des combinaisons optimales d'efforts déployés en matière d'éducation, d'intervention et d'information pourraient être identifiées.
# Chapitre 6

## Perspectives

### 6.1 Suite du modèle APC

Le projet Com2SiCa a produit beaucoup de donnés et ouvert de nouvelles perspectives. Ainsi, j'envisage à court et moyen terme plusieurs suites.

#### Exploitation des autres données issues des enquêtes de terrain

Toutes les données expérimentales n'ont pas encore été exploitées comme celles issues des enquêtes de terrain. Durant l'enquête menée au Havre, les participants étaient immergés dans une ambiance sonore qui était découpée en séquences (bruits courants dans la ville, explosion, cris). À partir de mesures physiologiques obtenues avec des montres connectées et un questionnaire pointu, D. Provitolo a défini à chaque séquence de l'enregistrement sonore l'état émotionnel de chaque participant et la trajectoire qu'il aurait pris. Je souhaiterais reprendre ces données pour affiner la connaissance des paramètres du modèle spatial et a-spatial.

#### D'autres modèles contrôlés envisagés

- Le modèle APC : Ce premier travail sur le contrôle prend en compte la formation de la population, l'effet de la présence des forces de l'ordre et la communication lors de l'événement. Toutefois, le deuxième contrôle ne nous satisfaisait pas. En effet, une intervention rapide avec le bon nombre de secouristes au bon endroit permet soit d'éviter un mouvement de panique soit, si ces mouvements se déclenchent, de les réguler, de les contrôler et d'éviter ainsi qu'ils se prolongent de façon excessive. Cette réflexion a conduit à modéliser les forces de l'ordre comme une nouvelle population. Nous avons ainsi initié avec les géographes l'étude des impacts de deux variables de contrôle sur les tendances comportementales : i) la présence de secouristes (et non simplement leur effet) ; ii) l'effet de l'information transmise au cours de l'événement (sirène, téléphone portable). Il serait également intéressant de comparer ce modèle contrôlé avec un modèle hybride construit à partir du modèle APC et d'un modèle à base d'agents [129]. Ce type de modélisation constitue une approche intermédiaire entre les modèles continus et les modèles à base d'agents de systèmes d'individus en interaction. Tout en gardant le modèle APC pour modéliser une population concernée par un événement catastrophique, les opérationnels seraient considérés en tant qu'individus ou agents suivant un ensemble de règles inhérent et pouvant communiquer entre eux. Ce type de modélisation est bien adapté dans le cas où il y a un petit nombre d'agents.

- Le modèle APC et les réseaux : Actuellement, l'environnement spatial est modélisé par un réseau, les noeuds représentant un lieu doté de propriétés regroupant un ensemble d'individus soumis au panel de comportements du modèle APC, les arêtes symbolisant les axes de communication et les interactions entre noeuds. Une méthodologie mixant la théorie des équations différentielles non linéaires et la théorie des réseaux complexes nous a permis de commencer à étudier l'impact de la topologie des réseaux et de leur connectivité sur la dynamique de propagation des différents comportements. Les résultats sont actuellement basés sur des simulations. Une méthode de contrôle comme proposée dans [216] permettrait d'obtenir des résultats théoriques et de mieux comprendre, par exemple, la meilleure localisation des secouristes pour réguler les effets de panique.

#### Un modèle à différentes échelles

Une autre perspective à se projet est de développer un modèle microscopique hybride reprenant le modèle APC et les modèles de force sociale. Dans ces derniers, chaque individu se déplace selon l'action de plusieurs forces externes et internes s'exerçant sur le piéton. Les forces externes proviennent des interactions physiques
dues à l'environnement tandis que les forces internes proviennent de l'objectif que s'est fixé le piéton. Les forces internes sont la réponse à une configuration ou à un état spécifique de l'environnement. Il est possible également de représenter la tendance à se tenir à l'écart du danger ou de considérer l'effet attractif d'une scène. L'accélération instantanée optimale de chaque individu est ainsi déterminée à partir des forces agissant sur le piéton, de son environnement et des obstacles. Dans le modèle proposé par [47], on retrouve deux processus que nous avons décrits dans le projet :

— le comportement provenant du stress interne pris en compte en liant le niveau de stress à la vitesse de l'individu

— le processus de contagion ou d'imitation.

Notre expertise dans la compréhension et la modélisation des comportements humains et leurs transitions
doit permettre d'enrichir ce type de modèle.

Pour les opérationnels, il est important de planifier et concevoir le mouvement collectif d'un grand nombre de personnes dans de nombreuses situations à risque telles que l'évacuation de bâtiments, de stades... Bien que les premières études scientifiques sur l'évacuation des piétons aient été menées dans les années 30, ce sujet reste d'actualité notamment dû à une connaissance partielle des comportements complexes des êtres humains lors d'un tel événement. Or, lors de ce projet, les différents comportements adoptés par les individus en situation de catastrophe ainsi que leurs transitions comme l'imitation ont été analysés et décrits. De même, les liens entre états comportementaux (alerte, contrôle et panique) ont été associés à des types de comportements physiques. Il serait désormais pertinent d'intégrer ces connaissances dans des modèles microscopiques qui pourraient proposer ainsi des simulations plus réalistes. Pour cela, j'envisage le développement d'un système multi-agent (SMA) qui permet de modéliser les effets d'interactions entre les individus. Chaque individu ou agent peut rassembler différents niveaux de fonctionnement (raisonnement, comportement réactif et proactif) et différents niveaux d'interaction (coordination) en les exploitant de manière autonome. Ces modèles permettent donc de considérer les interactions locales et physiques et la possibilité d'incorporer la prise de décision individuelle et les processus d'apprentissage. Ainsi, toutes les connaissances acquises au cours de ce projet peuvent être reprises pour enrichir le SMA. Le SMA permettrait d'explorer par les simulations le fonctionnement de l'ensemble en adoptant une approche ascendante, "bottom-up" pour faire émerger des propriétés plus complexes.

## 6.2 Suite sur l'identifiabilité et la diagnosticabilité

Au cours de mon travail de calibration à partir de données réelles, notamment avec le modèle APC, j'ai été confrontée à un problème d'estimation des paramètres même si le modèle était identifiable. Différentes définitions d'identifiabilité des paramètres existent dans la littérature et peuvent être classés en trois classes : i) l'identifiabilité a priori, ii) l'identifiabilité a posteriori, iii) l'identifiabilité pratique.

L'identifiabilité théorique ou structurelle des paramètres du modèle que j'ai étudié examine la question de l'existence et de l'unicité d'une solution au problème d'estimation des paramètres, dans un cadre idéalisé où

— le système et le modèle ont une structure identique (pas d'erreur de caractérisation);

— les données sont non bruitées;

— les signaux d'entrée et les temps de mesure peuvent être choisis à volonté.

Cependant, ce n'est qu'une condition nécessaire qui ne garantit pas l'estimation réussie des paramètres à partir de données réelles.

L'identifiabilité a posteriori considère les hypothèses suivantes :

— le système et le modèle ont une structure identique ;

— les données sont sans bruit.

L'identifiabilité a posteriori est un cas particulier de la distinguabilité de la sortie pour une collection finie d'observations sans bruit et un signal d'entrée [65, 66].

Enfin, l'identifiabilité pratique s'appuie sur la seule hypothèse que le système et le modèle ont une structure identique [67, 192]. Dans ce cas, cette notion tient compte du bruit et est généralement associée à un critère d'optimisation. Pour cette raison, cette classe d'identifiabilité est souvent liée à la théorie de l'optimisation en mathématiques.

Le problème d'estimation rencontré dans le cadre de mon travail est un problème d'identifiabilité a posteriori et peut se voir sur une étude de sensibilité des sorties du modèle par rapport aux paramètres. Par exemple, en utilisant les indices de Sobol, les sorties du modèle APC sont très peu sensibles à certains paramètres. Des travaux ont été menés pour étudier les liens entre l'identifiabilité et la sensibilité du modèle mais tous les liens n'ont pas encore été établis (voir [66]).

En ce qui concerne la diagnosticabilité, les signatures que nous avons proposées ont ouvert la voie à de nouvelles méthodes prometteuses. Toutefois, il reste encore une question fondamentale pour la diagnosticabilité
fonctionnelle : comment distinguer les trajectoires en fonction des fautes qui agissent sur le système et sur quel critère ? J'envisage pour cela deux suites possibles :

— étendre les définitions de signatures/diagnosticabilité au cadre ensembliste;

— travailler sur les polynômes entrées-sorties et sur leurs propriétés.

## 6.3 À plus long terme

En 2020, nous avons développé avec G. Cantin un modèle mathématique structuré en âge décrivant la dynamique des écosystèmes forestiers et prenant en compte le phénomène de la transpiration forestière [31]. Ce dernier est un élément important du cycle hydrologique car il exerce une influence sur la dynamique forestière (voir par exemple [128]). Or, deux éléments au moins sont susceptibles de perturber ce phénomène de transpiration forestière. Tout d'abord, le dérèglement climatique induit des augmentations ou des diminutions localisées des quantités de précipitations, ainsi qu'une augmentation moyenne de la température, elle-même susceptible d'accroître la fréquence ou l'intensité des événements de sécheresse. D'autre part, l'activité humaine et les actions de déforestation, en diminuant les surfaces de forêt, perturbent la contribution forestière au cycle hydrologique. Nous avons également développé dans [31] un réseau complexe afin d'appréhender la structure hétérogène de certaines étendues forestières, tout en intégrant les phénomènes d'interaction entre îlots forestiers par échange de la ressource en eau. Dans ces modèles de réseaux, nous avons considéré un graphe dont les sommets modélisent les îlots forestiers, et dont les arêtes correspondent aux interactions entre les différents îlots. Chaque sommet du graphe est alors couplé avec un système d'équations reproduisant la dynamique des îlots. L'ensemble qui résulte de cette modélisation décrit donc la dynamique spatio-temporelle complexe d'une zone forestière hétérogène et notamment l'état d'équilibre et de coexistence des diverses catégories d'arbres (cet état d'équilibre faisant l'objet d'une attention particulière puisqu'il correspond à l'état que l'on souhaite préserver ou restaurer), mais aussi celui de disparition de l'écosystème. Cette démarche de modélisation permet ainsi de reproduire les conditions favorables à l'équilibre de l'écosystème. Dans une seconde étape, nous avons intégré sous la forme de perturbations les multiples événements susceptibles de menacer cet équilibre. En particulier, nous avons considéré les actions humaines de déforestation et les changements climatiques, dont les variations de température ou de précipitations pour étudier la robustesse des équilibres de l'écosystème forestier. Pour poursuivre ce projet de modélisation, nous sommes en contact avec des géographes et des écologistes pour récupérer des données de la forêt Amazonienne, région du monde qui voit depuis plusieurs années une augmentation de la fréquence et de l'ampleur des événements extrêmes. Ils disposent de donner qualitatives (études des extrêmes climatiques tels que définis en climatologie, leurs réponses environnementales et leurs perceptions, en examinant les similitudes et les contrastes fournis par chacun de ces prismes) et quantitatives (évaporation, densité d'arbres au $km^2$). Afin d'étudier l'interconnexion entre climat et écosystèmes, d'évaluer les risques et la vulnérabilité liés à l'augmentation des extrêmes, nous avons proposé de concevoir avec eux un modèle mathématiques et de le calibrer. Ce modèle s'appuierait sur le modèle que nous avons développé et intégrerait les interactions entre l'écosystème forestier tropical et la dynamique climatique à partir des données forestières du réseau Tropical Managed Forest Observatory (TmFO) entre autre. Il simulerait les réponses éco-dynamiques des forêts aux extrêmes climatiques. Comme dans le projet Com2SiCa, nous leur avons proposé de créer avec eux un outil de support pour la gestion forestière et la prise de décision.
<figure><img src="image_64.png" /></figure>
Annexes
<figure><img src="image_65.png" /></figure>