# Etudes qualitatives et quantitatives de modèles paramétriques. De la théorie à la pratique avec un exemple en sciences humaines.

Nathalie Verdière

► To cite this version:

Nathalie Verdière. Etudes qualitatives et quantitatives de modèles paramétriques. De la théorie à la pratique avec un exemple en sciences humaines. Mathématiques [math]. Université du Havre (ULH), 2022. tel-03947466

HAL Id: tel-03947466

https://hal.science/tel-03947466

Submitted on 3 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Études qualitatives et quantitatives de modèles paramétriques. De la théorie à la pratique avec un exemple en sciences humaines.

MÉMOIRE PRÉSENTÉ EN VUE DE L'OBTENTION DU
DIPLÔME D'HABILITATION À DIRIGER DES RECHERCHES
EN MATHÉMATIQUES APPLIQUÉES PAR

NATHALIE VERDIÈRE

UNIVERSITÉ DU HAVRE

## Rapporteurs :

Denis Dochain
Marie Doumic
Dimitri Lefebvre
Professeur - Université Catholique de Louvain
Directrice de recherche (HDR) - laboratoire J-L. Lions, Paris
Professeur - Université Le Havre

## Composition du Jury :

François Boulier
Andrea De Gaetano
Antoine Perasso
Tarek Raïssi
Louise Travé-Massuyès
Lilianne Denis-Vidal
Professeur - École Polytechnique Universitaire, Lille
Directeur de recherche - CNR-IASI Biomatlab
Maître de conférences (HDR) - Université de Franche-Comté
Professeur - CNAM, Paris
Directrice de recherche - LAAS-CNRS, Toulouse
Maître de conférences (HDR) - Université de Technologie de Compiègne

## Référente :

Carine Jauberthie
Maître de Conférences (HDR) - Université Paul Sabatier, Toulouse III
21 juin 2022
<figure><img src="image_1.png" /></figure>
Remerciements

Je tiens tout d'abord à remercier vivement tous les membres du jury pour leur soutien et leur encouragement dans ma décision de soutenir mon habilitation à diriger le recherche. Marie Doumic, Denis Dochain et Dimi-tri Lefebvre qui ont accepté, sans hésitation, la lourde tâche d'être les rapporteurs de cette HDR. J'apprécie sincèrement l'intérêt qu'ils portent à mes travaux ainsi que leurs commentaires avisés. Je suis également très reconnaissante à Louise Travé-Massuyès, François Boulier, Andrea DeGaetano, Tarek Raïssi et Antoine Perasso de faire partie du jury. Je remercie Carine Jauberthie qui a acceptée d'être ma référente et pour ses conseils avisés sur la préparation du dossier.

Je voudrais encore une fois exprimer ma gratitude à ma directrice de thèse, Lilianne Denis Vidal, qui m'a suivie tout au long de ces années, m'a guidée et m'a toujours soutenue dans mon travail et ma carrière. Ce manuscrit n'existerait sûrement pas sans elle.

Ce travail n'existerait pas non plus sans les personnes avec qui j'ai eu le plaisir de travailler. Discuter, argu-
menter, (dé)construire, programmer sont des processus permettant d'aller jusqu'au bout des idées. Merci donc
(dans le désordre) Valentina, Nathalie, Damienne, Sébastien, Rodolphe, Guillaume, Shouzeng...

D'une façon plus intime, j'adresse mes remerciements à la personne qui a été présente depuis le début, qui
m'a aidée à gérer le quotidien et m'a soutenue, mon mari Dominique. Merci également à mes deux petits amours
pour leur bonne humeur et leur joie de vivre de tous les jours, Delphine et Adrien.
<figure><img src="image_2.png" /></figure>
# Table des matières

1 Introduction 2

2 Identifiabilité et estimation des paramètres : des systèmes EDO aux systèmes EDP 6

2.1 Modèle mathématique et quelques préliminaires indispensables 6

2.2 Premières définitions de l'identifiabilité 8

2.3 Identifiabilité relative 9

2.3.1 Ensembles semi-algébriques 10

2.3.2 L'algorithme IdentifiabilityTree 11

2.3.2.1 Concepts utiles 11

2.3.2.2 Algorithme IdentifiabilityTree 12

2.3.3 Croissance microbienne dans un réacteur 13

2.3.4 Le modèle de la transmission du chikungunya 14

2.4 Estimation de paramètres à partir des polynômes ES 15

2.4.1 Méthode d'estimation à partir du polynôme ES 16

2.4.2 Approche basée sur les distributions 16

2.4.2.1 Erreur due à la méthode 17

2.4.2.2 Erreur due à la méthode numérique 18

2.4.2.3 Erreur due au bruit 19

2.4.2.4 Application au modèle de transmission du chikungunya 20

2.5 Identifiabilité ensembliste 24

2.5.1 Concepts utiles 27

2.5.2 Définitions 27

2.5.3 Liens entre la SM-identifiabilité et d'autres concepts 28

2.5.3.1 Liens avec la ε-identifiabilité globale 28

2.5.3.2 Liens avec l'injectivité partielle 29

2.5.3.3 Partitionnement d'un domaine en sous-ensembles (μ-)SM-identifiables 30

2.5.4 Estimation des paramètres par des méthodes ensemblistes et propriétés 31

2.5.5 Robustesse 34

2.5.6 La ε-consistance 34

2.5.7 Estimation des paramètres dans le cadre ensembliste 37

2.5.7.1 Une enveloppe de la solution initiale 37

2.5.7.2 Application 38

2.6 Identifiabilité dans le cadre des équations aux dérivées partielles (EDP) 39

2.6.1 Notations 41

2.6.2 Le modèle considéré 41

2.6.3 Cas général 42

2.6.4 Résultats d'identifiabilité 43

2.6.5 Exemple 46

2.6.5.1 Résultats d'identifiabilité 48

2.6.5.2 Exemple de l'estimation de paramètres à partir des polynômes ES dans le cadre des EDP 48

2.7 Conclusion 50

3 Diagnosticabilité fonctionnelle et algébrique 52

3.1 Le modèle 53

3.2 Diagnosticabilité fonctionnelle 54

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

3.2.2 Lien entre identifiabilité et diagnosticabilité fonctionnelle 56

3.2.3 Un exemple complet : celui du réservoir 58
<table>
  <tbody>
    <tr>
      <td>3.3</td>
      <td>Diagnosticabilité algébrique</td>
      <td>60</td>
    </tr>
    <tr>
      <td>3.3.1</td>
      <td>Préliminaires</td>
      <td>61</td>
    </tr>
    <tr>
      <td>3.3.2</td>
      <td>Construction de la signature algébrique à partir du résumé exhaustif</td>
      <td>61</td>
    </tr>
    <tr>
      <td>3.3.3</td>
      <td>Critères pour différentier les signatures</td>
      <td>63</td>
    </tr>
    <tr>
      <td>3.3.4</td>
      <td>Application</td>
      <td>65</td>
    </tr>
    <tr>
      <td>3.3.5</td>
      <td>La détection de faute dans le cadre ensembliste</td>
      <td>67</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Conclusion</td>
      <td>68</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Réseaux</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Introduction</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Synchronisation de systèmes non orientés</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>Extension de la méthode de stabilité de connexion de graphes</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Synchronisations complètes</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Application</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Liens entre synchronisation et bifurcations</td>
      <td>75</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Conclusion</td>
      <td>77</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Croisement entre les mathématiques et les sciences humaines</td>
      <td>78</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Pourquoi avons nous choisi la modélisation macroscopique à base d'EDOs?</td>
      <td>79</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Le modèle Panic-Control-Reflex (PCR)</td>
      <td>80</td>
    </tr>
    <tr>
      <td>5.2.1</td>
      <td>D'une typologie des comportements au modèle mathématique</td>
      <td>80</td>
    </tr>
    <tr>
      <td>5.2.2</td>
      <td>Le modèle mathématique PCR et ses applications</td>
      <td>81</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Le modèle Alert-Panic-Control (APC)</td>
      <td>85</td>
    </tr>
    <tr>
      <td>5.3.1</td>
      <td>Catégories comportementales</td>
      <td>87</td>
    </tr>
    <tr>
      <td>5.3.2</td>
      <td>Présentation du modèle APC (Alert-Panic-Control)</td>
      <td>89</td>
    </tr>
    <tr>
      <td>5.3.3</td>
      <td>Une première calibration à partir des connaissances des thématiciens</td>
      <td>89</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Une calibration à partir de la réalité virtuelle</td>
      <td>91</td>
    </tr>
    <tr>
      <td>5.4.1</td>
      <td>Introduction</td>
      <td>91</td>
    </tr>
    <tr>
      <td>5.4.2</td>
      <td>L'expérimentation de réalité virtuelle (RV)</td>
      <td>92</td>
    </tr>
    <tr>
      <td>5.4.2.1</td>
      <td>La procédure de RV et l'échantillonnage</td>
      <td>92</td>
    </tr>
    <tr>
      <td>5.4.2.2</td>
      <td>Mesures</td>
      <td>94</td>
    </tr>
    <tr>
      <td>5.4.3</td>
      <td>Le lien entre le modèle mathématique et les données expérimentales</td>
      <td>94</td>
    </tr>
    <tr>
      <td>5.4.3.1</td>
      <td>Variables observées du modèle mathématique</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.4.3.2</td>
      <td>Modélisation du début et de la fin de la catastrophe</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Sensibilité des sorties, identifiabilité et nombre minimal de données d'observation</td>
      <td>97</td>
    </tr>
    <tr>
      <td>5.5.1</td>
      <td>Étude de sensibilité</td>
      <td>97</td>
    </tr>
    <tr>
      <td>5.5.2</td>
      <td>Identifiabilité</td>
      <td>98</td>
    </tr>
    <tr>
      <td>5.5.3</td>
      <td>Nombre minimal de données d'observation</td>
      <td>99</td>
    </tr>
    <tr>
      <td>5.5.4</td>
      <td>Estimation des paramètres</td>
      <td>99</td>
    </tr>
    <tr>
      <td>5.5.5</td>
      <td>Résultats</td>
      <td>100</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Un modèle mathématique, pour quoi faire?</td>
      <td>101</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Perspectives</td>
      <td>108</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Suite du modèle APC</td>
      <td>108</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Suite sur l'identifiabilité et la diagnosticabilité</td>
      <td>109</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>À plus long terme</td>
      <td>110</td>
    </tr>
    <tr>
      <td colspan="2">Annexes</td>
      <td>113</td>
    </tr>
    <tr>
      <td colspan="2">Annexe 1 : L'algorithme IdentifiabilityTree</td>
      <td>114</td>
    </tr>
    <tr>
      <td colspan="2">Annexe 2 – Outils basiques de l'analyse par intervalle</td>
      <td>118</td>
    </tr>
    <tr>
      <td>1</td>
      <td>Définitions basiques</td>
      <td>118</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Vecteurs et matrices d'intervalle</td>
      <td>118</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Pavage et sous-pavage</td>
      <td>119</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Opérations arithmétiques</td>
      <td>119</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Les fonctions d'inclusion</td>
      <td>120</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Image réciproque d'ensembles via l'analyse par intervalle</td>
      <td>121</td>
    </tr>
  </tbody>
</table>
Annexe 3 – Analyse de sensibilité 122

<table><tr><td>1</td><td>Rappels sur les indices de Sobol</td><td>122</td></tr><tr><td>1.1</td><td>Théorie</td><td>122</td></tr><tr><td>1.2</td><td>Principe des calculs</td><td>125</td></tr><tr><td>2</td><td>Résultats</td><td>126</td></tr></table>
<figure><img src="image_3.png" /></figure>
# Chapitre 1

# Introduction

Ce mémoire résume l'ensemble de mes travaux de recherches en mathématiques appliquées depuis ma thèse "Utilisation de l'algèbre différentielle pour analyser l'identifiabilité des paramètres de modèles décrits par des équations aux dérivées partielles", soutenue en 2005 (directrice de thèse : Lilianne Denis-Vidal, Ghislaine Joly-Blanchard). Ce mémoire est donc l'occasion de faire le bilan sur mon travail qui porte sur l'étude qualitative et quantitative des modèles dynamiques et d'esquisser quelques pistes de recherches amorcées ou futures.

Au cours de ma thèse au sein du laboratoire de mathématiques appliquées de l'UTC, j'ai travaillé sur des méthodes pour étudier l'identifiabilité a priori de modèles aux équations aux dérivées ordinaires (EDO) ainsi que sur le développement de méthodes numériques pour estimer les paramètres de ces modèles.

L'identifiabilité est une notion importante lorsqu'on construit un modèle mathématique et que l'on souhaite estimer des paramètres non directement accessibles par des mesures. En effet, elle assure qu'à une sortie du modèle mathématique correspond un seul jeu de paramètres et donc garantit les résultats numériques sur les paramètres. Un modèle non identifiable pourrait produire les mêmes sorties avec plusieurs, voir une infinité de jeux de paramètres, ce qui pose un problème lorsque ces derniers ont une signification physique comme les taux de contact infectieux dans les modèles épidémiologiques qui doivent être positifs. Pour étudier l'identifiabilité, un travail théorique est fait dans un cadre idéal sans considérer le bruit de mesure, en supposant avoir toutes les mesures nécessaires et où le modèle cherché correspond à celui du système étudié. Cette propriété, encore appelée identifiabilité structurelle est alors inhérente à la structure du modèle et constitue une condition nécessaire pour l'identifiabilité numérique. Cette propriété structurelle dépend des entrées et des sorties du modèle considéré, des conditions initiales comme l'ont démontré [61, 178, 196] et plus généralement des contraintes des paramètres [204].

Mon travail de thèse s'appuyait sur des relations différentielles obtenues à partir d'outils d'algèbre différentielle. Elles se sont révélées très riches et leur potentiel pas complètement exploité aussi bien au niveau théorique que numérique. Plusieurs questions restaient ouvertes, pouvait-on par exemple

- réécrire certains modèles mathématiques pour obtenir des relations plus simples à partir de l'algèbre différentielle ? En effet, le calcul formel peut produire des relations très compliquées en terme d'expressions et en terme d'ordre de dérivation;

- incorporer formellement les contraintes sur les paramètres et les conditions initiales pour obtenir un résultat certifié sur l'identifiabilité du modèle ? Les outils d'algèbre différentielle ne prennent pas en compte les inégalités qui peuvent changer les résultats d'identifiabilité;

- développer un algorithme donnant les paramètres clés pour obtenir l'identifiabilité d'une partie du modèle voire de tout le modèle;

- élargir cette notion à d'autres cadres comme le cadre ensembliste ou le cadre des équations aux dérivées partielles (EDP) et quelles seraient les limites ?

Après ma thèse, j'ai pu apporter quelques réponses à ces questions qui sont abordées dans le chapitre 2.

La diagnosticabilité de fautes est apparue comme une suite naturelle de l'étude d'identifiabilité en se plaçant dans le cas de fautes constantes. Cette propriété établit quelles fautes peuvent être discriminées étant donnés les capteurs placés dans le système. Cette étude, si elle est menée lors de la conception du système, peut permettre de mettre en place des méthodes de diagnostic, c'est à dire des méthodes de détection et d'isolation de fautes. Actuellement, ces propriétés sont très étudiées dans le cadre des systèmes autonomes, en particulier dans l'ارسpatial pour les envois de navettes autonomes à de longues distances [57]. Nous avons proposé de nouvelles pistes pour mener cette étude théorique et numérique à partir de relations différentielles non encore totalement exploitées dans ce domaine. Nous avons étendu certaines définitions et développé de nouvelles techniques d'analyse dans le cadre ensembliste mais également dans un cadre plus traditionnel. Nous avons également apporté quelques réponses sur le cas multi-faute, cas qui est rarement abordé dans la littérature. Cette thématique constituera le chapitre 3 de ce manuscrit.
Au cours de ces années, je me suis également intéressée aux réseaux et aux propriétés émergentes qui pouvaient résulter des interactions entre ces différents constituants. En particulier, j'ai travaillé sur

— l'extension d'un résultat de synchronisation dans un réseau dont les noeuds sont constitués de systèmes dynamiques non linéaires, non nécessairement identiques et dont les termes de couplage sont non linéaires ;

— le lien entre la synchronisation d'un système acyclique direct formé de systèmes de Hindmarsh-Rose et les bifurcation de(s) système(s).

Ces résultats sont présentés au chapitre 4.

Ce travail sur l'identifiabilité, la diagnosticabilité et la synchronisation ont répondu à des problématiques ponctuelles et concrètes. L'identifiabilité et la diagnosticabilité ont été au coeur de la tâche 2 du projet ANR MAGIC-SPS, tâche qui proposait d'étendre ces notions au cadre ensembliste dans le cas de modèles dynamiques incertains. Quant aux notions d'identifiabilité et de synchronisation, elles ont été utilisées dans le cadre d'un travail initié en 2013 avec le projet PEPS HuMaIn (Humanité-Mathématiques-sciences de l'Information) et qui a abouti au projet ANR Com2SiCa en 2018. Ce travail est présenté au chapitre 5. Ces projets consistaient à mieux comprendre les différents comportements et leurs enchaînements lors d'une catastrophe puis à développer un modèle mathématique les décrivant en accord avec les géographes (Damienne Provitolo, Edwide Dubos-Paillard) puis les psychologues (Oscar Navarro, Abdel-Halim Boudoukha). Cette thématique connaît un vif écho dans notre société depuis plusieurs années qui sont marquées par une augmentation significative du nombre de catastrophes qu'elles soient naturelles, sociétales ou technologiques. Une meilleure compréhension des réactions humaines est devenue un enjeux majeur dans nos sociétés. D'ailleurs, ces tendances ne devraient pas s'inverser dans les années à venir puisque les facteurs de risque restent nombreux : changements climatiques, tensions géopolitiques, les risques émergents liés aux avancées technologiques, la détérioration de l'environnement et la pression urbaine etc. Or le comportement des individus et des populations lors de catastrophes reste encore méconnu en de nombreux points par les opérationnels comme la complexité des comportements véritablement adoptés [52], leur enchaînement, leur dynamique et interdépendance [164].

En faisant une revue de la littérature, différents modèles mathématiques à différentes échelles ont été développés pour décrire une dynamique de foule dans un contexte de situations extrêmes (voir. [11]). On peut citer par exemple :

— les modèles de force sociale [93] : les piétons sont identifiés à des particules soumises i) à une force extérieure due aux obstacles (murs), ii) une force extérieure qui prend en compte les interactions avec les autres piétons (typiquement la tendance à préserver une certaine distance avec les voisins). L'équation décrit la vitesse de l'individu qui va évoluer dans son environnement pour atteindre une vitesse désirée. C'est un modèle essentiellement réactif, il n'y a pas de gestion simultanée des interactions de nature différente ;

— les modèles granulaires [132] : chaque individu est identifié à un disque d'une certaine taille et ne pouvant pas s'intercepter avec son voisin. L'approche s'appuie sur la vitesse désirée par chaque individu (la vitesse qu'ils auraient s'ils étaient seuls) et le champ de vitesse global est défini comme la plus proche vitesse faisable sans que les disques s'interceptent. Ce modèle décrit la position de chaque individu, on considère l'ensemble des vitesses désirées et le cône de vitesses faisables. Les outils sont basés sur l'analyse convexe et la théorie des inclusions différentielles. Le problème est que les individus ont tendance à suivre une stratégie purement égoïste sans considérer les autres ;

— les automates cellulaires [9] : les agents ou personnes sont identifiés à des particules. Il y a discrétisation en espace et en temps. Pour l'espace, les agents sont localisés dans des cellules d'une grille fixée et pour le temps, l'évolution consiste en une succession de pas. Les positions des particules sont mises à jour l'une après l'autre de manière stochastique.

— les modèles macroscopiques : ils décrivent l'évolution en temps et en espace de la densité et de la vitesse moyenne d'une foule et sont basés sur des systèmes d'équations aux dérivées partielles. Ces modèles sont basés sur la dynamique des fluides et prennent en compte la conservation de la masse.

Deux points importants viennent limiter la portée des modèles que nous avons trouvés dans la littérature. Le premier est qu'ils ont été construits à partir de modèles de physiques (fluide, particules) et s'intéressent essentiellement aux dynamiques de foule le plus souvent paniquée (voir [11, 93]). Or, d'après les psychologues, les géographes et les neuro-scientifiques, lors de catastrophes, plusieurs comportements concurrents sont observés produisant une dynamique pouvant se révéler complexe. Le deuxième point est que ces modèles ne prennent pas ou très peu en compte les changements émotionnels des individus au cours d'une catastrophe comme : i) les transitions naturelles spécifiques à chaque individu selon son expérience et son passé ; ii) les transition dues à l'imitation ou la comparaison sociale qui est un important facteur dans les processus de changement de comportements [77]. Pour ces raisons, nous avons fait le choix de proposer un nouveau modèle basé sur les échanges avec les géographes des laboratoires de Géoazur (Nice) et de Géographie-Cités (Paris) et psychologues du laboratoire de LLPL (Nantes). Nous avons décidé de reprendre la formalisation des modèles SIR (Susceptibles-Infected-Removed) car ces modèles compartimentaux permettent de décrire en temps l'évolution du nombre
d'individus dans chaque compartiment. Un compartiment représentant un état comportemental. Ces modèles permettent de prendre en compte différents types d'interactions entre ces différents comportements (dans notre cas les transitions intrinsèques et les transitions dues à l'imitation ou la contagion sociale). Il est également possible de spatialiser ce type de modèle soit en utilisant la théorie des équations aux dérivées partielles pour une diffusion locale, soit en utilisant la théorie des réseaux pour une diffusion globale. Toutefois, cette thématique portant sur l'étude des comportements en situation de catastrophe a été rarement quantifiée car toute tentative de mesure est difficile à mettre en oeuvre en pratique. Dans le chapitre 5 de ce manuscrit, je vais présenter nos avancées dans cette thématique et en particulier

— comment le modèle a été construit avec les géographes et les psychologues;

— comment nous avons exploité les données existant dans la littérature et les connaissances des thématiciens pour donner un sens et un ordre de grandeur aux paramètres;

— comment nous avons exploité les données provenant d'une expérimentation de réalité virtuelle pour affiner la connaissance des valeurs des paramètres. Cette expérimentation a été réalisée en lieu clos et contrôlé au Laboratoire de Psychologie des Pays de la Loire (Nantes) en 2020 et s'appuie sur un scénario, celui de l'arrivée d'une vague de tsunami d'origine sismique sur une plage de Nice [19];

— les informations et les connaissances que ce type de modèle peut apporter.

Je finirai par un dernier chapitre qui présente les perspectives que j'envisage à court et moyen terme.
<figure><img src="image_4.png" /></figure>
# Chapitre 2

## Identifiabilité et estimation des paramètres : des systèmes EDO aux systèmes EDP

Les contributions présentées dans cette partie ont été publiées dans [100, 101, 173, 200, 204].

L'identifiabilité est une notion importante lorsqu'on construit un modèle mathématique et que l'on souhaite estimer des paramètres non directement accessibles par des mesures. Cette notion structurelle connaît depuis plusieurs années un intérêt croissant dans différentes domaines scientifiques [160, 207]. La question se pose lorsqu'on cherche à calibrer un modèle mathématique à partir de données réelles bruitées. Toutefois, si la structure du modèle n'est pas identifiable, le résultat de la procédure numérique mise en place pour estimer les paramètres peut donner des valeurs erronées, c'est à dire des valeurs ne correspondant pas à la vraie valeur des paramètres, ce qui pose un problème lorsque les paramètres ont une signification physique. Par exemple, l'algorithme utilisé pourrait renvoyer un taux de transition négatif intervenant dans un modèle biologique alors qu'il devrait être positif. Pour apporter une réponse à cette problématique, l'identifiabilité du modèle assure, qu'à une sortie du modèle mathématique, correspond un seul jeu de paramètres. Pour étudier l'identifiabilité, un travail théorique est fait dans un cadre idéal. Cette propriété, encore appelée identifiabilité structurelle est alors inhérente à la structure du modèle est constitue une condition nécessaire pour l'identifiabilité numérique. Plusieurs méthodes pour étudier l'identifiabilité structurelle peuvent être trouvées dans la littérature [36, 37, 60, 61, 64, 80, 99, 125, 161, 184, 209]. La propriété d'identifiabilité dépend des entrées et des sorties que l'on peut considérer dans le modèle mais également des conditions initiales [61, 178, 196] et plus généralement des contraintes sur les paramètres.

Ce chapitre présente mon travail dans ce domaine. Tout d'abord, je vais présenter la notion d'identifiabilité relative qui généralise la notion d'identifiabilité et permet de déterminer les paramètres clés permettant d'avoir l'identifiabilité de certains paramètres voire du modèle complet. Puis, je montrerai comment nous avons généralisé ce concept au cadre ensembliste et aux équations aux dérivées partielles. Dans chaque partie, j'ai travaillé sur le développement de méthodes basées sur des outils d'algèbre pour vérifier ces propriétés et conduisant à des méthodes d'estimation de paramètres.

### 2.1 Modèle mathématique et quelques préliminaires indispensables

Dans cette première section est présenté le modèle paramétrique utilisé au cours de tout ce chapitre ainsi que quelques rappels d'algèbres différentielles indispensables pour la suite.

Le modèle que nous allons considérer dans ce chapitre est le suivant :

$$ \left\{ \begin{array}{l} \dot{x}(t, \theta) = f(x(t, \theta), u(t), \theta), \\ y(t, \theta) = h(x(t, \theta), \theta), \end{array} \right. \qquad (2.1) $$

où

- Le vecteur de paramètres $\theta = (\theta_1, \dots, \theta_q)$ appartient au sous-ensemble $\mathcal{D}$ de $\mathbb{R}^q$ où $\mathcal{D}$ est a priori un ensemble connu de paramètres admissibles.

- La fonction $f(\cdot, \cdot, \theta)$ est supposée réelle, rationnelle et analytique<fnref n="1" /> sur $M \times C^\infty(\mathbb{R}, \mathbb{R}^r)$, où $M$ est un

<footnote n="1">1. Cette hypothèse n'est pas restrictive puisque de nombreux modèles peuvent être réduits à un modèle rationnel et analytique par un changement de variables (voir [8])</footnote>
ensemble ouvert connexe de $\mathbb{R}^n$ et tel que $x(t, \theta) \in M$ pour chaque $\theta \in \mathcal{D} \subseteq \mathbb{R}^q$ et chaque $t \in [0, T]$,

— $h(., \theta)$ est supposée réelle, rationnelle et analytique sur $M$,

— $u$ est un vecteur de dimension $r$ de fonctions régulières (infiniment différentiables). Dans le cas d'un système non contrôlé, $u$ est égal à 0,

Soit $x_0 = (x_{0,i})_{i=1,...,n}$ le vecteur de conditions initiales pour le système (2.1) et dont certaines composantes peuvent dépendre de paramètres à estimer.

**Exemple 2.1.1.** Reprenons l'exemple de l'article [196] qui décrit la capacité d'un récepteur du macrophage à recevoir des macromolécules solubles spécifiques aux transporteurs d'enzymes. Il s'écrit :

$$ egin{cases} 
\dot{x}_1 = \alpha_1(x_2 - x_1) - \frac{V_m x_1}{1 + x_1}, \\ 
\dot{x}_2 = \alpha_2(x_1 - x_2), \\ 
x_1(0) \in [0.62, 0.63], x_2(0) = 0, 
\end{cases} 
\quad (2.2) $$

où $x_1$ (resp. $x_2$) est la concentration d'enzymes à l'extérieur (resp. à l'intérieur) du macrophage et $p = (\alpha_1, V_m, \alpha_2)$ est le vecteur de paramètres inconnus à identifier. Le paramètre $\alpha_1$ est le taux de transfert du compartiment 1 (ou compartiment central, en pratique le plasma) au compartiment 2 (ou compartiment périphérique qui représente la partie du liquide extracellulaire accessible) et $\alpha_2$ est le taux de transfert du compartiment 2 au compartiment 1.

Si on suppose observer la concentration d'enzymes à l'extérieur du macrophage, en dérivant la première équation et en substituant $x_2$ et $\dot{x}_2$ par leur expression dépendant de $x_1$, on obtient le polynôme différentiel ($y := x_1$ et en supposant que $\alpha_1 \neq 0$),

$$ P(y, u) = \ddot{y}(1+y)^2 + \gamma_1 \dot{y}(1+y)^2 + \gamma_2 y(1+y) + \gamma_3 \dot{y} $$

$$ \text{avec } \gamma_1 = \alpha_1 + \alpha_2, \gamma_2 = \alpha_2 V_m \text{ et } \gamma_3 = V_m. $$

Nous verrons que ce polynôme a des propriétés très intéressantes pour mener des études sur les paramètres. Toutefois, pour obtenir de telles relations, les manipulations à la main peuvent être difficiles. L'algorithme de Rosenfeld-Groebner implémenté dans le package DifferentialAlgebra de Maple permet de le retrouver (voir [20, 22]). Par exemple, notons $\theta = (\alpha_1, \alpha_2, V_m)^T$. Considérons l'ordre d'élimination $[\theta] \prec [y] \prec [x_1, x_2]$ ($\theta$ est considéré comme une variable de dérivées nulles)$^2$. L'ordre d'élimination choisi consiste à éliminer les variables $x_1$, $x_2$ et leurs dérivées au profit de $y$ (et de ses dérivées) puis de $\theta$. À partir de cet ordre d'élimination, l'algorithme de Rosenfeld-Groebner triangularise le système d'équations (2.2) en utilisant les opérations d'addition, de soustraction, de multiplication et de dérivation. Dans notre exemple, il renvoie en particulier les deux listes suivantes :

$$ L_1 := \{-\alpha_2 x_2 + \alpha_2 y + \dot{x}_2, x_1 - y, V_m y + y \dot{y} + \dot{y}, \alpha_1\} $$

et

$$ L_2 := \{x_1 - y, \alpha_1 x_2 y - \alpha_1 y^2 - V_m x_1 + \alpha_1 x_2 - \alpha_1 y + (-y-1)\dot{y}, \\ -V_m \alpha_1 \alpha_2 y^2 + \alpha_1^2 y^2 \dot{y} - \alpha_1 \alpha_2 y^2 \dot{y} - V_m \alpha_1 \alpha_2 y + 2 \alpha_1^2 y \dot{y} - 2 \alpha_1 \alpha_2 y \dot{y} + \alpha_1 y^2 \ddot{y} + V_m \alpha_1 \dot{y} \\ + \alpha_1^2 \dot{y} - \alpha_1 \alpha_2 \dot{y} + 2 \alpha_1 y \ddot{y} + \alpha_1 \ddot{y}\}. \quad (2.3) $$

La liste $L_1$ est le cas particulier où $\alpha_1 = 0$ tandis que la liste $L_2$ est le cas général. En factorisant la dernière relation de $L_2$, on retrouve le polynôme $P$.

De façon générale, étant donné l'ordre d'élimination $[\theta] \prec [y, u] \prec [x]$ appliqué au système (2.1), l'algorithme divise le système initial en plusieurs sous-systèmes formés d'équations et d'inéquations. L'union de ces systèmes est encore appelée chaîne régulière. Les solutions de chacun de ces sous-systèmes sont solutions du système (2.1) et seulement l'un d'entre eux est le cas général, les autres étant des solutions particulières. Ainsi, l'ensemble solution du système initial est l'union de l'ensemble des solutions de tous les sous-systèmes. Le cas général est de la forme [62]

$$ C(\theta) = \{\dot{\theta}_1, \dots, \dot{\theta}_q, P_1(y, u, \theta), \dots, P_m(y, u, \theta), Q_1(x, y, u, \theta), \dots, Q_n(x, y, u, \theta)\} $$

où, pour $i = 1, \dots, m$, le leader du polynôme $P_i$, c'est à dire la variable dont l'ordre de dérivation est le plus élevé, est $y_i$.

<footnote n="2">2. Un classement < désigne un ordre total sur l'ensemble des dérivées et le bloc $[v]$, où $v$ est une variable, signifie que $v$ et ses dérivées sont ordonnées par le classement < de la façon suivante $v < \dot{v} < \ddot{v} < \dots$. Le symbole < indique quant à lui un classement entre les blocs dans le sens "classement inférieur à" entre deux blocs.</footnote>
Sous certaines conditions techniques presque toujours vérifiées [62], le polynôme $P_i$ appelé **polynôme Entrée-Sortie (ES)** peut s'écrire sous la forme :

$$P_i(y, u, \theta) = m_{i,0}(y, u) + \sum_{j=1}^{d_i} \gamma_{i,j}(\theta) m_{i,j}(y, u) \quad \text{et tel que} \quad P_i(y, u, \theta) = 0. \qquad (2.4)$$

Les $\gamma_{i,j}(\theta)$ pour $j = 1, \dots, d_i$ sont rationnels en $\theta$ et tels que $\gamma_{i,a} \neq \gamma_{i,b}$ ($a \neq b$). $(m_{i,j}(y, u))_{1 \leq j \leq d_i}$ est une suite de polynômes différentiels en $y$ et $u$. Enfin, $m_{i,0}$ vérifie $m_{i,0} \neq 0$.

Dans [62], les auteurs ont montré qu'il y a autant de polynômes ES que de sorties du modèle. Par la suite, pour simplifier les notations, nous considérerons qu'il n'y a qu'une sortie pour le modèle et donc un seul polynôme ES. L'indice $i$ dans (2.4) sera omis.

**Remarque 1.** Contrairement à d'autres auteurs [53, 225, 125], nous avons choisi ces polynômes plutôt que ceux obtenus à partir de l'algorithme de Ritt (voir [125] pour plus de détails). En effet, ces derniers requièrent beaucoup plus de manipulations des équations du modèle. Ainsi, les relations obtenues sont plus complexes en terme de nombre d'expression et l'ordre de dérivation plus élevé.

## 2.2 Premières définitions de l'identifiabilité

De manière générale, j'ai adopté les définitions classiques suivantes [8, 136, 138, 196] :

**Définition 1.** 1. Soit $\theta \in \mathcal{D}$ un vecteur de paramètres. Si il existe un contrôle $u$ tel que le système d'équations $y(t, \theta) = y(t, \tilde{\theta})$, avec $\tilde{\theta} \in \mathcal{D}$, possède

- seulement une solution $\theta = \tilde{\theta}$, le vecteur de paramètres $\theta$ est dit **globalement identifiable**,

- un nombre fini de solutions distinctes alors $\theta$ est dit **localement identifiable**,

- une infinité de solutions alors $\theta$ est dit **non identifiable**.

2. Une composante $\theta_i$ de $\theta$ est **globalement identifiable** si pour tout $\tilde{\theta} \in \mathcal{D}$ tel que $y(t, \theta) = y(t, \tilde{\theta})$ implique que $\theta_i = \tilde{\theta}_i$.

3. Le modèle (2.1) est dit **globalement identifiable** si toutes les composantes du vecteur de paramètres $\theta$ sont globalement identifiables.

Pour tester cette propriété, il est possible de se ramener à une étude d'injectivité sur les coefficients du polynôme ES (2.4), le résultat est le suivant et se trouve dans [62] :

**Proposition 1.** Supposons que le déterminant fonctionnel $\Delta P(y, u) = \det(m_j(y, u), j = 1, \dots, d)$ n'est pas identiquement nul$^3$. Le modèle (2.4) est globalement identifiable si et seulement si la fonction

$$\begin{array}{rcl} \phi: & \mathcal{D} & \to \mathbb{R}^d \\ & \theta & \mapsto (\gamma_j(\theta))_{1 \le j \le d} \end{array} \qquad (2.5)$$

est injective.

**Remarque 2.** 1. Dans la première version de la méthode décrite dans [62], les conditions initiales n'étaient pas considérées. Elles ont été intégrées par la suite démontrant leur rôle sur les résultats d'identifiabilité ([61, 178, 196]).

2. Plusieurs approches ont été proposées pour résoudre le système réel $\phi(\theta) = \phi(\tilde{\theta})$ [12, 62, 177, 196] comme l'utilisation des bases de Groebner. Toutefois, celles-ci ne permettent pas de considérer d'éventuelles inégalités vérifiées par les paramètres.

3. Si il y a M observations ($M \ge 1$) donc M polynômes ES alors pour chaque polynôme différentiel dont le déterminant fonctionnel est non identiquement nul, ses coefficients peuvent être ajoutés à l'image de la fonction $\phi$.

4. La suite $(\gamma_j(\theta))_{1 \le j \le d}$ est également appelée **résumé exhaustif**.

Souhaitant répondre à des questions plus précises sur la structure du modèle, je me suis intéressée à la généralisation des définitions classiques d'identifiabilité ainsi qu'au développement des méthodes associées.

<footnote n="3">3. Cette hypothèse consiste à vérifier l'indépendance linéaire des $m_j(y, u)$, $j = 1, \dots, d$ en vérifiant que le déterminant fonctionnel calculé à partir du Wronskien est non identiquement nul. Dans le cadre de l'algèbre différentielle, cette condition consiste à vérifier qu'il n'est pas dans l'idéal obtenu après l'élimination des variables d'état. En pratique, cette condition peut-être vérifiée grâce au package Differential-Algebra de Maple.</footnote>
## 2.3 Identifiabilité relative

La définition d'identifiabilité relative que nous avons proposée est apparue naturellement en étudiant des systèmes dynamiques non identifiables.

**Exemple 2.3.1.** Considérons le modèle de croissance microbienne dans un réacteur étudié dans les articles [37, 74] auquel est ajouté le facteur de maintenance m. Celui-ci correspond à la masse de substrat qu'une unité de biomasse consomme en une seule unité de temps. Le modèle dynamique pour le processus de croissance est gouverné par les deux équations suivantes :

$$ egin{cases} 
\dot{x}(t) = \frac{\mu s(t)x(t)}{K_S + s(t)} - mYx(t), \\ 
\dot{s}(t) = \frac{-\mu s(t)x(t)}{Y(K_S + s(t))}, 
\end{cases} 
\qquad (2.6) $$

où x est la concentration de micro-organismes, s la concentration de substrat limitant la croissance, μ la vitesse maximale de réaction, K_S la constante de Michaelis-Menten et Y le coefficient de rendement. Étant donné la nature des paramètres, μ, K_S, m et Y sont supposés être des réels positifs. Comme il est plus facile de mesurer x(t) que s(t) à l'intérieur du réacteur, la sortie mesurée est supposée être y := x. Le modèle étudié dans [37] est un cas particulier de (2.6) en posant K_d = mY.

Ce modèle n'est pas identifiable à partir de la seule mesure de concentration de microorganismes. En effet, l'algorithme de Rosenfeld-Groebner appliqué à ce système avec l'ordre d'élimination $[θ] \prec [y] \prec [x, s]$, $θ := (μ, K_S, m, Y)$ retourne le polynôme ES suivant :

$$ P(y) := (Y^3 m^3 - 2Y^2 m^2 \mu + Y m \mu^2) y^3 + (3Y^2 m^2 - 4Y m \mu + \mu^2) y^2 \dot{y} \\ + \mu K_S Y (y \ddot{y} - \dot{y}^2) + (3Ym - 2\mu) y \dot{y}^2 + \dot{y}^3 = 0. $$

Le déterminant fonctionnel $det(P(y)) = det(y^3, y^3\dot{y}, y\ddot{y} - \dot{y}^2, y\dot{y}^2)$ est non nul et la fonction $\phi: θ = (μ, K_S, m, Y) \mapsto (Y^3m^3 - 2Y^2m^2μ + Ymμ^2, 3Y^2m^2 - 4Ymμ + μ^2, μK_SY, 3Ym - 2μ)$ n'est pas injective. Donc, d'après la proposition 1, le modèle n'est pas identifiable.

Sur ce simple modèle, plusieurs questions viennent naturellement :

- étant donné la non identifiabilité du modèle, quels paramètres doit-on déterminer pour obtenir l'identifiabilité de paramètres significatifs non directement mesurables, voir l'identifiabilité de tout le système ?

- Quels rôles peuvent avoir les contraintes sur les résultats d'identifiabilité des paramètres ?

- Comment intégrer naturellement ces contraintes ainsi que les conditions initiales dans l'étude d'identifiabilité ?

En ce qui concerne la première question, les modèles construits dans le cadre de la biologie peuvent contenir des paramètres significatifs et importants pour l'expérimentateur. Leur estimation expérimentale représente parfois un coût financier non négligeable et si le modèle n'est pas identifiable, leur valeur estimée à partir du modèle n'est pas garantie. Dans ce cas de figure, il est courant de chercher et d'utiliser des combinaisons identifiables de paramètres à travers des approches numériques/statistiques ([27, 44, 72, 134, 213]), des méthodes symboliques-numériques ([41, 45]) et des approches algébriques ([6, 38, 62, 74, 75, 136, 137, 138, 139, 140, 196]) pour re-paramétriser le modèle initial et ainsi obtenir un modèle identifiable. Toutefois, cette re-paramétrisation introduit des paramètres qui n'ont plus de sens physique ou biologique conduisant à un modèle dont les résultats sont plus difficilement interprétables. Or, l'estimation de certains paramètres ou l'ajout de contraintes peuvent permettre d'avoir l'identifiabilité du modèle initial. C'est dans la perspective de laisser à l'expérimentateur le meilleur choix de stratégie que nous avons proposé la notion d'identifiabilité relative. Cette définition a conduit au développement d'un algorithme dont le résultat est un arbre donnant les paramètres clés pour avoir l'identifiabilité de certains paramètres voir du modèle complet.

La première contribution de ce travail est théorique. Comme l'identifiabilité d'un paramètre peut dépendre de l'identifiabilité d'autres paramètres, nous avons étendu la définition classique d'identifiabilité et introduit la notion d'identifiabilité relative, c'est à dire l'identifiabilité d'un paramètre donné quand d'autres sont supposés connus. Nous avons ensuite généralisé un résultat donné dans [61]. Ce résultat lie les trajectoires de sorties et les solutions de systèmes algébriques construits à partir des polynômes ES. Plus précisément, nous avons établi l'équivalence entre l'identifiabilité relative d'un paramètre et l'ensemble solution d'un système d'équations et/ou d'inéquations qui doit être vide. Cet ensemble solution peut être testé grâce à des logiciels de calcul algébrique. Cette approche théorique a permis d'obtenir une approche automatisée pour étudier l'identifiabilité (relative) des paramètres d'un modèle.

La seconde contribution est un algorithme basé sur cette équivalence et dont le résultat est un arbre. Sa construction se fait à partir d'outils semi-algébriques. Ces derniers, implémentés dans des logiciels de calcul
formel, n'ont pas encore été utilisés jusqu'à présent pour certifier l'identifiabilité structurelle de modèles dynamiques. Il permet d'intégrer les conditions initiales et les contraintes sur les paramètres. L'algorithme retourne un ensemble de listes, chacune d'entre elle donnant un ensemble de paramètres clés pour obtenir l'identifiabilité d'un ou des paramètres du modèle. Notre algorithme peut prendre en compte le domaine de définition des paramètres et des relations algébriques provenant par exemple des conditions initiales. Le résultat théorique obtenu sur l'équivalence assure que l'ensemble des listes fournies par notre algorithme est exhaustive pour le modèle considéré.

Dans beaucoup de modèles, il est naturel de considérer un sous-ensemble *C* de l'ensemble admissible des paramètres *D*. En effet, due à la signification même des paramètres, des conditions initiales ou pour assurer l'existence de solutions, des contraintes supplémentaires doivent être considérées. Par exemple, dans les études de chaînes trophiques, le nombre de proies mangées par les prédateurs ne peut pas être plus grand que la population totale de proies et en mécanique, une masse ne peut pas être négative. Pour voir comment intégrer ces contraintes qui apparaissent naturellement quand on étudie un problème, je vais rappeler la notion d'ensemble semi-algébrique qui sera à la base du travail présenté dans cette partie.

## 2.3.1 Ensembles semi-algébriques

Par la suite, nous supposons que les contraintes sur le vecteur de paramètres $\theta$ peuvent être formulées à l'aide d'équations et/ou d'inégalités algébriques. La définition d'un ensemble semi-algébrique est la suivante (cf [10]):

**Définition 2.** *Un ensemble de solutions réelles d'un ensemble fini d'équations polynômiales et/ou d'inégalités polynômiales est appelé ensemble semi-algébrique.*

L'introduction d'ensembles semi-algébriques permet d'utiliser les librairies de calculs formels comme la librairie *Raglib* de Maple [73] ou la librairie *SemiAlgebraicSetTools* [26, 39, 218]. Une introduction sur la façon de résoudre des systèmes semi-algriques à l'aide Maple est disponible à l'adresse suivante https://www.csd.uwo.ca/mmorenom/Publications/cdlmpxxx.pdf.

Dans la suite de cette partie, on suppose que :

- $C(\theta)$ est l'ensemble des équations et inéquations algébriques vérifiées par les composantes du vecteur de paramètres $\theta \in D$ du modèle;

- l'ensemble $C$ est un ensemble semi-algébrique défini par $C(\theta)$.

**Exemple 2.3.2.** Ainsi, pour l'exemple 2.3.1, nous pouvons prendre $C(\theta) = \{\mu > 0, K_S > 0, m > 0, Y > 0\}$.

La définition de l'identifiabilité relative d'un paramètre particulier par rapport à d'autres paramètres identifiables (éventuellement aucun) est la suivante :

**Définition 3.** *1. Le paramètre $\theta_{a_{r+1}}$ est dit relativement identifiable par rapport à l'ensemble de paramètres $\{\theta_{a_1}, \dots, \theta_{a_r}\}$ si il existe un contrôle u tel que, pour tout $\tilde{\theta} = (\tilde{\theta}_1, \dots, \tilde{\theta}_q) \in D$ dont les composantes vérifient $\theta_{a_1} = \theta_{a_1}, \dots, \theta_{a_r} = \theta_{a_r}$, l'équation $y(t,\theta) = y(t,\tilde{\theta})$ implique $\theta_{a_{r+1}} = \theta_{a_{r+1}}$.*

*2. L'étude de l'identifiabilité relative du système (2.1) est la détermination de l'identifiabilité relative de chaque paramètre par rapport à chaque ensemble de paramètres construit à partir de l'ensemble $\{\theta_1, \dots, \theta_q\}$.*

Par définition, $\theta_{a_{r+1}}$ n'est pas relativement identifiable par rapport à l'ensemble de paramètres $\{\theta_{a_1}, \dots, \theta_{a_r}\}$ (resp. pas identifiable) si, pour tout contrôle $u$, il existe $\tilde{\theta} = (\tilde{\theta}_1, \dots, \tilde{\theta}_q) \in D$ dont les composantes satisfont $\tilde{\theta}_{a_1} = \theta_{a_1}, \dots, \tilde{\theta}_{a_r} = \theta_{a_r}, \theta_{a_{r+1}} \neq \tilde{\theta}_{a_{r+1}}$ et $y(t,\theta) = y(t,\tilde{\theta})$.

Cette définition apparaît naturellement quand on souhaite déterminer quels paramètres doivent être estimés pour obtenir des paramètres identifiables et éventuellement l'identifiabilité du modèle.

**Remarque 3.** *1. On retrouve la définition classique de l'identifiabilité structurelle lorsque le paramètre $\theta_{a_{r+1}}$ est relativement identifiable par rapport à l'ensemble vide des paramètres. Par conséquent, l'identifiabilité du système (2.1) est un cas particulier de l'identifiabilité relative par rapport à toutes les composantes de $\theta$.*

*2. Les définitions 1 et 3 peuvent être étendues à des définitions locales en imposant que $\tilde{\theta}$ appartient à un voisinage $\nu(\theta) \subset D$ de $\theta$ plutôt que $D$.*

La proposition suivante est une conséquence directe de la définition 3.

**Proposition 2.** *Si $\theta_i$ est relativement identifiable par rapport à l'ensemble éventuellement vide $\mathcal{P} \subset \{\theta_1, \dots, \theta_q\}$ alors $\theta_i$ est relativement identifiable par rapport à tout sous-ensemble de $\{\theta_1, \dots, \theta_q\}$ contenant $\mathcal{P}$.*
Nous allons désormais voir comment réduire l'étude de l'identifiabilité relative à un problème d'ensemble semi-algébrique.

Pour prendre en compte les conditions initiales et/ou les contraintes supplémentaires sur les paramètres, le domaine de définition de la fonction $\phi$ (voir (2.5)) est le sous-ensemble $\mathcal{C}$ de $\mathcal{D}$. La proposition suivante donne le lien entre l'identifiabilité relative et le résumé exhaustif.

**Proposition 3.** Le paramètre $\theta_{a_{r+1}}$ est relativement identifiable par rapport à l'ensemble éventuellement vide des paramètres $\{\theta_{a_1}, \dots, \theta_{a_r}\}$ si et seulement si, pour tout $\theta = (\theta_j)_{j=1,...,q} \in \mathbb{R}^q$ et $\tilde{\theta} = (\tilde{\theta}_j)_{j=1,...,q} \in \mathbb{R}^q$ :

$$ 
\left\{ 
\begin{array}{l}
\theta \in \mathcal{C} \\
\tilde{\theta} \in \mathcal{C} \\
\theta_{a_1} = \tilde{\theta}_{a_1} \\
\vdots \\
\theta_{a_r} = \tilde{\theta}_{a_r} \\
\phi(\theta) = \phi(\tilde{\theta})
\end{array}
\right.
\quad \Rightarrow \quad \theta_{a_{r+1}} = \tilde{\theta}_{a_{r+1}} \qquad (2.7)
$$

**Remarque 4.** La proposition 3 peut-être étendue à une version locale en supposant que $\tilde{\theta}$ appartient à un voisinage semi-algébrique $V_{\theta}$ de $\theta$.

Le but de ce qui va suivre est de proposer une procédure automatique basée sur les outils semi-algébriques pour prouver l'implication (2.7) incluant les contraintes sur les paramètres.

Du système composant le membre de gauche dans l'implication (2.7), on obtient un système composé d'équations polynômiales et d'inégalités définissant un ensemble semi-algébrique. En effet, $\mathcal{C}$ est supposé être un ensemble semi-algébrique et le système $c_{i,j}(\theta) = c_{i,j}(\tilde{\theta})$, $i = 1, \dots, s$, $j = 1, \dots, d_i$ est équivalent à un système algébrique dont les inéquations polynômiales correspondent aux dénominateurs du système. Ces inéquations polynômiales sont alors rajoutées à $C(\theta)$ et $C(\tilde{\theta})$.

Soit $\theta = (\theta_1, \dots, \theta_q) \in \mathbb{R}^q$, $\tilde{\theta} = (\tilde{\theta}_1, \dots, \tilde{\theta}_q) \in \mathbb{R}^q$ et $0 \le r \le q$. Notons :

$$ S_{\theta_{a_1}, \dots, \theta_{a_r}} = C(\theta) \cup C(\tilde{\theta}) \cup \{\theta_{a_i} = \tilde{\theta}_{a_i} \mid 1 \le i \le r\} \cup \{c_{i,j}(\theta) = c_{i,j}(\tilde{\theta}) \mid 1 \le i \le s, 1 \le j \le n_i\}. $$

Le corollaire suivant est basé sur un outil classique en calcul formel : l'astuce de Rabinowitsch ([10, 168]). Pour montrer que le système algébrique admet une solution $u \neq 0$ où $u$ est une des indéterminée, il suffit d'ajouter l'équation $uv - 1 = 0$ au système, où $v$ est une nouvelle indeterminée, et de prouver que le nouveau système a une solution. Il est ensuite possible de tester si les ensembles semi-algébriques sont vides.

**Corollaire 1.** Soit $(\theta_1, \dots, \theta_q)$ et $(\tilde{\theta}_1, \dots, \tilde{\theta}_q)$ dans $\mathcal{C}$. Les conditions suivantes sont équivalentes :

1. le paramètre $\theta_{a_{r+1}}$ du système (2.1) est relativement identifiable par rapport à un ensemble éventuellement vide de paramètres $\{\theta_{a_1}, \dots, \theta_{a_r}\}$;

2. le système $S_{\theta_{a_1}, \dots, \theta_{a_r}} \cup \{v(\theta_{a_{r+1}} - \tilde{\theta}_{a_{r+1}}) - 1 = 0\}$, où $v$ est une nouvelle variable n'a aucune solution réelle $(\theta_{a_1}, \dots, \theta_{a_q}, \tilde{\theta}_{a_1}, \dots, \tilde{\theta}_{a_q}, v) \in \mathbb{R}^{2q+1}$.

**Remarque 5.** De la même façon qu'à la proposition 3, le corollaire 1 peut-être étendu à une version locale en substituant $C(\tilde{\theta})$ par $C(\theta) \cup V_{\theta}$ dans le corollaire où $V_{\theta}$ a été défini dans la remarque 4.

## 2.3.2 L'algorithmme IdentifiabilityTree

### 2.3.2.1 Concepts utiles

**Notation 1.** $[\theta_{a_1}, \dots, \theta_{a_r}]$ est la liste<fnref n="4" /> de $r$ éléments distincts ($1 \le r \le q$) pris parmi $\{\theta_1, \dots, \theta_q\}$.

**Définition 4.** 1. Soit $[\theta_{a_1}, \dots, \theta_{a_r}]$ une liste de $r$ éléments distincts de l'ensemble $\{\theta_1, \dots, \theta_q\}$. La liste vide et toute liste $[\theta_{a_1}, \dots, \theta_{a_i}]$ ($1 \le i \le r-1$) sont appelées préfixe de $[\theta_{a_1}, \dots, \theta_{a_r}]$.

2. Soit $l_1$ et $l_2$ deux listes. La liste concaténée $l_1$ cat $l_2$ est la liste obtenue en regroupant les deux listes $l_1$ et $l_2$ l'une à la suite de l'autre.

Pour pouvoir distinguer les paramètres identifiables et les non identifiables dans la liste, la notation suivante a été introduite.

<footnote n="4">4. En calcul formel, une liste est une suite finie d'éléments.</footnote>
**Notation 2.** Soit $l = [\theta_{a_1}, \dots, \theta_{a_r}]$ une liste de $r$ éléments pris dans l'ensemble $\{\theta_1, \dots, \theta_q\}$.

Si le paramètre $\theta_{a_1}$ n'est pas identifiable, il sera noté $\theta_{a_1}$ dans la liste $l$.

Plus généralement, si $\theta_{a_i}$ ($2 \le i \le r$) n'est pas relativement identifiable par rapport à $\{\theta_{a_1}, \dots, \theta_{a_{i-1}}\}$, nous
le noterons $\emptyset_{a_i}$ dans $l$.

**Exemple 2.3.3.** Pour l'exemple 2.3.1 avec $C(\theta) = \{\mu > 0, K_S > 0, m > 0, Y > 0\}$, la liste $[\mu, K_S, Y, m]$ a été
obtenue. Cette liste signifie entre autre que

- $\mu$ est identifiable mais pas $K_S$,

- $K_S$ n'est pas relativement identifiable par rapport à l'ensemble $\{\mu\}$, c'est à dire que l'identifiabilité de $\mu$
  et la résolution de $\phi(\theta) = \phi(\tilde{\theta})$ n'implique pas l'identifiabilité de $K_S$,

- Y est relativement identifiable par rapport à l'ensemble $\{\mu, K_S\}$ lorsque $K_S$ est identifiable : l'identifiabilité de $\mu$ et $K_S$ implique l'identifiabilité de Y et m.

Notons également que la liste $[\mu, K_S]$ est un préfixe de la liste $[\mu, K_S, Y]$ par exemple.

Par la suite, pour alléger les notations, un paramètre sera dit relativement identifiable par rapport à une
liste de paramètres s'il est identifiable par rapport à l'ensemble correspondant de paramètres.

L'algorithme développé est un arbre de décision basé sur l'utilisation de listes de paramètres. L'avantage
de ce type d'algorithme est qu'il évite un grand nombre de tests sur les ensembles semi-algébriques manipulés.
Il construit toutes les listes possibles formées des $\theta_1, \dots, \theta_m$ en indiquant les paramètres non identifiables des
listes.

Définissons d'abord la notion d'arbre d'identifiabilité qui sera la sortie de l'algorithme.

**Définition 5.** L'arbre d'identifiabilité $\mathcal{T}$ de (2.1) est l'ensemble de toutes les listes possibles de paramètres pris parmi $\{\theta_1, \dots, \theta_q\}$ tels que, pour toute liste $b \in \mathcal{T}$ et tout préfixe $p$ de $b$, $p$ est suivi dans $b$ par un paramètre identifiable par rapport à $p$ si il en existe un.

Si les listes de $\mathcal{T}$ ont un préfixe en commun, les tests fait sur ce préfixe ne seront fait qu'une seule fois. De
plus pour éviter des tests inutiles, les deux remarques suivantes sont utilisées dans l'algorithme.

**Remarque A** (Conséquence de la proposition 2) Si un préfixe $p$ de la liste $\mathcal{T}$ est suivie par des paramètres identifiables, ces paramètres peuvent être permutés dans cette liste puisqu'ils réfèrent à la même information : chacun d'entre eux est relativement identifiable par rapport à $p$.

**Remarque B** Si des listes complétant un préfixe donné $p_1$ en éléments de $\mathcal{T}$ sont connues, ces listes complètent n'importe quelle permutation de $p_1$ en éléments de $\mathcal{T}$.

**Exemple 2.3.4.** Toujours avec l'exemple 2.3.1 avec $C(\theta) = \{\mu > 0, K_S > 0, m > 0, Y > 0\}$ et la liste $[\mu, K_S, Y, m]$,
la remarque A implique que comme dans la liste, $K_S$ est suivi des deux paramètres identifiables Y et m, l'algorithme ne testera qu'une fois l'identifiabilité relative de Y et m par rapport à la liste $[\mu, K_S]$.

En ce qui concerne la remarque B, $p_1 = [\mu, K_S]$ est un préfixe de la liste $[\mu, K_S, Y, m]$. Si la liste $[K_S, \mu]$ apparaît
au cours de l'algorithme, il ajoutera directement la liste $[K_S, \mu, Y, m]$ à $\mathcal{T}$.

Ainsi, l'identifiabilité relative d'un paramètre par rapport à un préfixe $p$ dépend seulement de l'ensemble des
paramètres apparaissant dans $p$.

Cette dernière remarque conduit naturellement à un arbre transversal pré-ordonné [89], ce qui permet de compléter des préfixes avec des listes de $\mathcal{T}$ déjà calculées. De cette façon, des tests d'ensemble vide peuvent être évités diminuant le coût global du calcul de $\mathcal{T}$.

### 2.3.2.2 Algorithme IdentifiabilityTree

Rappelons que les tests d'identifiabilité sont réalisés en appliquant le corollaire 1. Pour tester l'identifiabilité relative de $\theta_{a_{r+1}}$ par rapport à un ensemble $\{\theta_{a_1}, \dots, \theta_{a_r}\}$, on teste l'existence d'une solution non nulle d'un système composé d'équations polynômiales et d'inégalités construites à partir de la fonction $\phi$, des contraintes des paramètres et des conditions initiales.

L'algorithme *IdentifiabilityTree* est un algorithme récursif qui prend en entrée un préfixe $p$ d'une liste de $\mathcal{T}$
et calcule toutes les listes de $\mathcal{T}$ admettant $p$ comme préfixe, initialement $p = []$. Les listes de $\mathcal{T}$ produites par
l'algorithme sont stockées dans l'ensemble $\mathcal{T}'$ initialement vide.

Durant les appels récursifs, le préfixe $p$ est complété d'abord avec tous les paramètres identifiables par rapport
à $p$ (l'ordre n'a pas d'importance d'après la remarque A), ensuite par chaque paramètre non identifiable. Les
préfixes ainsi construits sont alors complétés par des appels récursifs de l'algorithme. Quand c'est possible, la
remarque B est utilisée pour compléter les préfixes de $\mathcal{T}'$.

Les étapes de l'algorithme sont résumées ci-dessous :

IdentifiabilityTree ($p$)

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

1. Tant qu'il existe des paramètres identifiables par rapport à $p$ faire

1.1 Compléter $p$ avec un paramètre identifiable $\theta_v$;

1.2 Vérifier si la remarque B peut être appliquée pour compléter $p$ cat $[\theta]$ en de nouvelles listes de l'arbre d'identifiabilité. Si c'est le cas, ajouter ces nouvelles listes à $\mathcal{T}'$, arrêter cette boucle Tant que et ne pas faire l'étape 2.

2. Si il existe des paramètres non identifiables par rapport à $p$, pour chacun de ces paramètres, faire

2.1 Compléter $p$ avec un paramètre non identifiable $\emptyset$;

2.2 Si la remarque B peut-être appliquée pour compléter $p$ cat $[\emptyset]$ avec de nouvelles listes de l'arbre d'identifiabilité, ajouter ces nouvelles listes à $\mathcal{T}'$. Sinon, faire un rappel récursif de l'algorithme avec $p$ cat $[\emptyset]$ comme entrée.

La proposition suivante donne la complexité de l'algorithme.

**Proposition 4 (Complexité).** Le nombre de tests effectués par l'algorithme pour vérifier que les ensembles semi-algébriques sont vides est borné par $(2q - v + 2)2^{v-1}$ où $m$ (resp. $v$) est le nombre de paramètres (resp. de paramètres non identifiables) du système (2.1).

### 2.3.3 Croissance microbienne dans un réacteur

Reprenons l'exemple (2.3.1) sur la croissance de micro-organismes dans un réacteur. L'algorithme *IdentifiabilityTree* retourne avec $C(\theta) = \{\mu > 0, K_s > 0, m > 0, Y > 0\}$ l'arbre complet donné à la figure 2.1.

<figure><img src="image_6.png" /><figcaption>FIGURE 2.1 – Arbre transversal retourné par l’Algorithm IdentifiabilityTree</figcaption></figure>

La sortie de l'algorithme est $\mathcal{T} = \{[\mu, K_S, Y, m], [\mu, \mathcal{Y}, K_S, m], [\mu, \mathcal{X}, K_S, Y]\}$. Ces listes montrent que $\mu$ est identifiable et qu'aucun des autres paramètres n'est identifiable. Remarquons que si un des paramètres parmi $K_S, m$ ou $Y$ est connu alors le système est identifiable.

**Remarque 6.** 1. Dans cet exemple, le rôle des contraintes imposées sur les paramètres joue un rôle crucial pour établir l'identifiabilité du système (2.3.1). En effet, si les contraintes $C(\theta)$ ne sont pas considérées, i.e. si $C(\theta) = \emptyset$, l'algorithme renvoie

$$ \mathcal{T} = \{[\mu, K_S, \mathcal{Y}, \mathcal{X}], [\mu, K_S, \mathcal{X}, \mathcal{Y}], [\mu, \mathcal{Y}, K_S, \mathcal{X}], [\mu, \mathcal{X}, K_S, \mathcal{Y}], [\mu, \mathcal{X}, \mathcal{Y}, K_S], [\mu, \mathcal{X}, \mathcal{Y}, K_S]\}. $$

Contrairement à la situation précédente, les trois paramètres doivent être identifiés pour que le modèle soit identifiable.

2. Pour évaluer le rôle crucial des conditions initiales, introduisons les nouveaux paramètres inconnus $x_0, s_0, xp_0, sp_0$ correspondant respectivement à $x(0), s(0), \dot{x}(0), \dot{s}(0)$ et ajoutons à $C(\theta)$ les équations obtenues à partir du système (2.6) évaluées en $t = 0$.

Le calcul de $\mathcal{T}$ montre que dans toutes les listes, le paramètre $\mu$ est identifiable et que

— soit deux paramètres parmi $x_0, s_0, xp_0$ ne sont pas identifiables et exactement un parmi $K_S, Y, m$ ne l'est pas,
— soit trois paramètres parmi $x_0, s_0, sp_0, xp_0$ ne sont pas identifiables.

La sortie de l'algorithme donne ainsi les choix possibles sur les paramètres du modèle pour obtenir l'identifiabilité de chacun des paramètres.

3. De la même façon, sans considérer les conditions initiales, mais en supposant que les sorties $x$ et $s$ sont observées, l'arbre d'identifiabilité est réduit à $\{[K_S, Y, m, \mu]\}$. Ce qui signifie que le modèle est identifiable.

4. Les auteurs des articles [37, 74] ont prouvé que dans le cas $K_d = mY$, les paramètres et les combinaisons de paramètres $x(0), K_d m Y, s(0)/Y, s(0)/K_s$ sont identifiables, alors que les paramètres pris individuellement $m, Y, s(0), K_s$ ne sont pas identifiables. De plus, en utilisant les combinaisons de paramètres, ils en ont déduits que si $s(0)$ ou $K_s$ sont connus alors tous les paramètres sont identifiables. L'algorithme *Identifiability Tree* nous donne directement le résultat de façon automatique.

### 2.3.4 Le modèle de la transmission du chikungunya

Le second exemple concerne un modèle épidémiologique décrivant la transmission du chikungunya à l'homme. Il est repris de l'article [145] et est intéressant de part ses non linéarités et son nombre de paramètres.

Le système est composé de deux sous-systèmes. Le premier décrit l'évolution des trois états biologiques du cycle de vie du moustique quant au deuxième sous-système, construit à partir d'un modèle de population, il décrit la transmission du virus du chikungunya entre la population du moustique et la population humaine. Un humain est infecté lorsqu'un moustique infecté le pique et un moustique est infecté lorsqu'il pique un humain infecté. Ce modèle suppose qu'il n'y a pas de transmission verticale pour les moustiques et les hommes, c'est à dire qu'il n'y a pas de transmission aux générations futures. Le modèle ré-écrit en terme de densité est le suivant :

$$ 
\left\{ 
\begin{array}{l} 
\left\{ 
\begin{aligned} 
E'(t) &= bA(t) \left(1 - \frac{E(t)}{K_E}\right) - (s+d)E(t) \\ 
L'(t) &= sE(t) \left(1 - \frac{L(t)}{K_L}\right) - (s_L+d_L)L(t) \\ 
A'(t) &= s_L L(t) - d_m A(t) 
\end{aligned} 
\right. \\ 
\\ 
\left\{ 
\begin{aligned} 
S'_H(t) &= -(b_H + \beta_H I_m(t)) S_H(t) + b_H \\ 
I'_H(t) &= \beta_H I_m(t) S_H(t) - (\gamma + b_H) I_H(t) \\ 
I'_m(t) &= - \left( s_L \frac{L(t)}{A(t)} + \beta_m I_H(t) \right) I_m(t) + \beta_m I_H(t) 
\end{aligned} 
\right. 
\end{array} 
\right. 
\quad (a) 
\\ 
\quad (b) 
\end{array} 
\right. 
\quad (2.8) 
$$

Dans le système (2.8)(a), $E$ est le nombre d'oeufs, $L$ le nombre de larves et $A$ est le nombre d'adultes femelles. En effet, seules les moustiques femelles transmettent la maladie. $S_H$ correspond au nombre d'humains susceptibles d'être infectés, $I_H$ au nombre d'humains infectés et $I_m$ au nombre de moustiques infectés.

Le vecteur de paramètres inconnus est $\theta = (k_E, k_L, b_H, \beta_H, \beta_m, d_L, d_m, \gamma, b, d, s, s_L)$ où

— $b$ est le taux de ponte intrinsèque des femelles, $s$ (resp. $s_L$) est le taux de transfert de $E$ à $L$ (resp. de $L$ à $A$);

— $K_E$ et $K_L$ sont les capacité d'accueil du gîte en oeufs et en larves respectivement;

— $d, d_L$ et $d_m$ sont les taux de mortalités des oeufs, des larves et des adultes respectivement;

— $b_H$ est le taux de natalité;

— $\gamma$ est le taux de transfert des infectés vers les immunisés;

— $\beta_H$ (resp. $\beta_m$) est le taux de contact infectieux entre humains susceptibles et vecteurs (resp. moustiques susceptibles et humains).

Naturellement, les paramètres réels du système (2.8) sont supposés positifs ce qui conduit à considérer l'ensemble $C(\theta) = \{b_H > 0, \beta_H > 0, \gamma > 0, s_L > 0, \beta_m > 0, d > 0, d_L > 0, d_M > 0, b > 0, s > 0, K_E > 0, K_L > 0, b_H > 0\}$. Par la suite, nous supposons qu'aucun des 12 paramètres est connu.

Le système est défini sur $\Delta \times \Omega$ où

$$ \Delta = \left\{ (E(t), L(t), A(t)) \in (\mathbb{R}^+)^3 \mid \forall t \in [0, +\infty), \begin{array}{l} 0 \le E(t) \le K_E \\ 0 \le L(t) \le K_L \\ 0 \le A(t) \le \frac{s_L}{d_m} K_L \end{array} \right\} \quad (2.9) $$

et

$$ \Omega = \left\{ (S_H(t), I_H(t), I_m(t)) \in (\mathbb{R}^+)^3 \mid \forall t \in [0, +\infty), \begin{array}{l} 0 \le S_H(t) + I_H(t) \le 1 \\ 0 \le I_m(t) \le 1 \end{array} \right\}. \quad (2.10) $$
<table><caption>TABLE 2.1 – Ensemble des listes fournies par l’algorithme IdentifiabilityTree</caption>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, k<sub>E</sub>, d, s, b, d<sub>L</sub>, s<sub>L</sub>],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, k<sub>E</sub>, d, s, d<sub>L</sub>, b, s<sub>L</sub>],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, k<sub>E</sub>, d, s, s<sub>L</sub>, b, d<sub>L</sub>],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, b, k<sub>E</sub>, d, d<sub>L</sub>, s, s<sub>L</sub>],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, b, d, K<sub>E</sub>, d<sub>L</sub>, s, s<sub>L</sub>],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, b, d<sub>L</sub>, K<sub>E</sub>, d, s, s<sub>L</sub>],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, b, s, K<sub>E</sub>, d, d<sub>L</sub>, s<sub>L</sub>],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, b, s<sub>L</sub>, K<sub>E</sub>, d, d<sub>L</sub>, s],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, d, K<sub>E</sub>, s, b, d<sub>L</sub>, s<sub>L</sub>],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, d, K<sub>E</sub>, s, d<sub>L</sub>, b, s<sub>L</sub>],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, d, K<sub>E</sub>, s, s<sub>L</sub>, b, d<sub>L</sub>],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, d<sub>L</sub>, sl, k<sub>E</sub>, b, d, s],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, d<sub>L</sub>, s<sub>L</sub>, b, K<sub>E</sub>, d, s],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, d<sub>L</sub>, s<sub>L</sub>, d, K<sub>E</sub>, b, s],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, d<sub>L</sub>, s<sub>L</sub>, s, K<sub>E</sub>, b, d],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, s, K<sub>E</sub>, d, b, d<sub>L</sub>, s<sub>L</sub>],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, s, K<sub>E</sub>, d, d<sub>L</sub>, b, s<sub>L</sub>],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, s, K<sub>E</sub>, d, s<sub>L</sub>, b, d<sub>L</sub>],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, s<sub>L</sub>, d<sub>L</sub>, k<sub>E</sub>, b, d, s],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, s<sub>L</sub>, d<sub>L</sub>, b, K<sub>E</sub>, d, s],</td>
  </tr>
  <tr>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, s<sub>L</sub>, d<sub>L</sub>, d, d, K<sub>E</sub>, b, s],</td>
    <td>[K<sub>L</sub>, b<sub>H</sub>, β<sub>H</sub>, β<sub>m</sub>, d<sub>m</sub>, γ, s<sub>L</sub>, d<sub>L</sub>, s, K<sub>E</sub>, b, d]</td>
  </tr>
</table>

L'analyse de stabilité est détaillée dans [145].

Les biologistes, en plaçant des pièges dans lesquels les femelles vont pondre peuvent avoir une estimation de l'évolution du nombre de larves de moustiques chaque semaine. Ensuite, lors d'épidémies importantes, les personnels de santé remontent le nombre d'infectés aux instances nationales. C'est pourquoi, nous avons supposé que les sorties du modèle sont le nombre de larves $y_1 := L$, le nombre d'humains susceptibles d'être infectés $y_2 := S_H$ et le nombre d'humains infectés $y_3 := I_H$.

Dans le but d'appliquer l'algorithme de Rosenfeld-Groebner au système (2.8), la dernière équation de (2.8)(b) est multipliée par la fonction A. Cet algorithme avec l'ordre d'élimination $[y_1, y_2, y_3] \prec [S_H, I_H, L, I_M, E, A]$ retourne trois équations différentielles liant les paramètres inconnus et les fonctions $y_1, y_2, y_3$. L'une de ces équations n'a pas de coefficient constant $c_{ij}$ et par conséquent, ses coefficients ne peuvent être estimées qu'à une constante multiplicative près à partir de données expérimentales (voir [196]). Cette équation est alors divisée par un de ses coefficients non nul, $K_E K_L$, pour obtenir un polynôme de la forme (2.4).

**Remarque 7.** Dans [196], il a été montré que sous certaines conditions techniques, nous retrouvions les mêmes polynômes ES avec l'ordre d'élimination $[y_1, y_2, y_3] \prec [S_H, I_H, L, I_M, E, A]$ et l'ordre d'élimination $[\theta] \prec [y_1, y_2, y_3] \prec [S_H, I_H, L, I_M, E, A]$. L'avantage du premier est qu'il est plus rapide et aboutit plus facilement lorsque le système présente plusieurs non linéarités. En effet, moins de variables sont manipulées.

À partir de ces polynômes ES écrits sous la forme (2.4), nous avons obtenu un résumé exhaustif contenant 74 polynômes algébriques dépendant des paramètres. Notre algorithme a été appliqué au système obtenu à partir de $C(\theta)$ et la fonction $\phi$. Malgré les coûts de calcul, la taille du résumé exhaustif et donc le nombre d'expressions algébriques, l'algorithme a fourni l'arbre d'identifiabilité. Il a retourné l'ensemble donné à la table 2.1$^5$.

Pour les sorties considérées, $y_1, y_2$ et $y_3$, l'ensemble des listes montre que $K_L, b_H, \beta_H, \beta_M, d_m, \gamma$ sont identifiables. De plus, il indique quelles autres paires de paramètres doivent être estimées pour obtenir un modèle identifiable. Beaucoup d'autres informations peuvent être extraites de ces listes. Par exemple, si un expérimenterneur veut estimer le paramètre non identifiable $K_E$ à partir des données expérimentales, il doit estimer soit $d$ soit $s$.

## 2.4 Estimation de paramètres à partir des polynômes ES

L'estimation des paramètres est généralement effectuée en optimisant une fonction objectif sur l'espace des paramètres [110, 210]. Or, la plupart des méthodes d'optimisation qui sont raisonnables en temps nécessitent une première estimation. Elle peut être obtenue soit par des tests dans le domaine industriel, soit par une procédure numérique basée sur des mesures [21, 63, 126, 181, 196]. La qualité de cette première estimation

<footnote n="5">5. Le temps de calcul pour obtenir l'arbre d'identifiabilité est de 1634 s avec un Intel(R) Core(TM) i7 2.5 GHz processor avec 8 GO de RAM.</footnote>
initiale garantit que le processus d'optimisation prendra un temps raisonnable et ne convergera pas vers une solution locale irréaliste. Il est possible d'obtenir cette première estimation à partir des relations polynômiales. Par exemple, les auteurs de [181] utilisent les intégrales itératives dans le domaine des fréquences alors que [126] exploite les fonctions de modulation mais sur des exemples de modèles très simples.

Toutes ces méthodes utilisent de façon avantageuse les propriétés de l'intégrale. En effet, dans la plupart des cas, leur utilisation permet de diminuer l'ordre des dérivées dans les relations différentielles. L'estimation des dérivées surtout à partir du degré 3 est un problème mal posé. Les intégrales atténuent également l'effet de certains bruits [81, 195]. Enfin, elles permettent d'intégrer les conditions initiales apportant ainsi de l'information au problème d'optimisation à résoudre.

Dans [203], nous avons proposé une méthode basée sur l'utilisation des distributions. Dans cet article, nous avons défini un cadre rigoureux en étudiant les erreurs dues à la méthode, au choix de la méthode d'intégration et au bruit. Les erreurs que nous avons obtenues dépendent du choix des fonctions tests et de leur support. Ce travail théorique donne un degré de liberté supplémentaire pour améliorer les résultats numériques.

## 2.4.1 Méthode d'estimation à partir du polynôme ES

Dans cette partie est présentée la méthode classique utilisant directement le polynôme ES. Pour simplifier les notations, la méthode est présentée sur un seul polynôme et l'indice *i* est omis. De plus, nous supposons que le modèle étudié est identifiable.

À partir d'un ensemble de mesures discrètes ($t_i$)$_{1 \le i \le M}$ ($M \ge q$), le polynôme ES est évalué en chaque $t_{n_0}, \dots, t_{n_1}$, $1 \le n_0 < n_1 \le M$ conduisant à la résolution d'un système linéaire rectangulaire

$$A\gamma = b \qquad (2.11)$$

avec

$$A = (A_{i,j})_{i=n_0, \dots, n_1 \atop j=1, \dots, q} = (m_j(y(t_i), u(t_i)))_{i=n_0, \dots, n_1 \atop j=1, \dots, q} \qquad (2.12)$$

$$b = (b_i)_{i=n_0, \dots, n_1} = -(m_0(t_i))_{i=n_0, \dots, n_1}, \gamma = (\gamma_j(\theta))_{j=1, \dots, q}$$

Classiquement, la solution de ce système est trouvée en résolvant le problème de moindres carrés suivant :

$$\min_{(\gamma_1(\theta), \dots, \gamma_q(\theta))} \sum_{i=n_0}^{n_1} \left( \sum_{j=1}^{q} A_{i,j} \gamma_j(\theta) - b_i \right)^2 . \qquad (2.13)$$

## 2.4.2 Approche basée sur les distributions

Les fonctions tests sont des fonctions indéfiniment dérivables définies sur un support compact. Par la suite, les fonctions tests seront notées $\psi_i$ et leur support $supp(\psi_i) = [t_i - \epsilon, t_i + \epsilon]$, $\epsilon > 0$. Supposons que les fonctions $u, y_j, x_i, j = 1, \dots, m, i = 1, \dots, n$ vérifient l'hypothèse suivante : (H) : elles sont $C^\infty$ dans $]0, T[$ et sont égales à zéro sur $]-\infty, 0[$, $]T, +\infty[$. Elles coïncident dans $]0, T[$ avec les fonctions qui sont $C^\infty$ dans $[0, T[$.

Après avoir multiplié le polynôme ES par la fonction test $\psi_i$ et avoir intégré par parties sur $[t_i - \epsilon, t_i + \epsilon]$, une nouvelle relation est obtenue. Elle peut se réécrire :

$$\int_{t_i - \epsilon}^{t_i + \epsilon} P(y, u, \theta)(s) \psi_i(s) ds = \int_{t_i - \epsilon}^{t_i + \epsilon} m_0(y, u)(s) \psi_i(s) ds + \sum_{j=1}^{q} \gamma_j(\theta) \int_{t_i - \epsilon}^{t_i + \epsilon} m_j(y, u)(s) \psi_i(s) ds = 0. \qquad (2.14)$$

À partir de la relation (2.14), on peut construire le système rectangulaire suivant ($1 \le n_0 < n_1 \le M$) :

$$A^{\psi} \gamma^{\psi} = b^{\psi} \qquad (2.15)$$

avec

$$A^{\psi} = (A_{i,j}^{\psi})_{i=n_0, \dots, n_1 \atop j=1, \dots, q} = \left( \int_{t_i - \epsilon}^{t_i + \epsilon} m_j(y, u)(s) \psi_i(s) ds \right)_{i=n_0, \dots, n_1 \atop j=1, \dots, q}, \qquad (2.16)$$

$$b^{\psi} = (b_i^{\psi})_{i=n_0, \dots, n_1} = \left( - \int_{t_i - \epsilon}^{t_i + \epsilon} m_0(y, u)(s) \psi_i(s) ds \right)_{i=n_0, \dots, n_1} .$$
Dans ce cas, le problème de moindres carrés à résoudre est le suivant :

$$ \min_{(\gamma_1(\theta)^\psi, \dots, \gamma_q(\theta)^\psi)} \sum_{i=n_0}^{n_1} \left( \sum_{j=1}^{q} A_{i,j}^\psi \gamma_j(\theta)^\psi - b_i^\psi \right)^2 \quad (2.17) $$

Pour simplifier les notations, nous posons $m_j = m_j(y, u)$ pour $j = 1, \dots, q$. De plus, la norme infinie, $\|\cdot\|_\infty$, induite par la norme matricielle sera utilisée [167].

Dans les applications numériques, la factorisation QR a été utilisée pour résoudre (2.11) et (2.15). Plutôt que de résoudre un problème de moindres carrés basé sur un système rectangulaire de dimension $M \times q$, la factorisation QR consiste à résoudre un problème équivalent basé sur l'inversion d'une matrice carrée inversible de dimension $q \times q$. C'est pourquoi, par la suite, les erreurs seront données dans le cas particulier de systèmes carrés. Un exemple de calcul d'erreur a été donné dans [203].

Cette remarque justifie les deux premières hypothèses suivantes. Les deux autres hypothèses concernent les fonctions tests.

**Hypothèses 1.** 1. Les matrices $A$ et $A^\psi$ sont carrées et de dimension $q \times q$.

2. $A$ et $A^\psi$ sont inversibles.

3. Les fonctions tests $\psi_i$, $i = 1, \dots, q$ sont supposées avoir la même valeur maximale notée $\psi_{\max}$, valeur atteinte en $t_i$.

4. Pour tout $i, j = 1, \dots, q$, $m_j \psi_i \in C^2([t_i - \epsilon, t_i + \epsilon])$.

Nous allons tout d'abord donner l'erreur due à l'approximation du système (2.11) par le système (2.15). La seconde erreur portera sur la méthode d'intégration utilisée pour estimer les intégrales dans (2.15) et finalement la troisième sur le bruit du signal.

### 2.4.2.1 Erreur due à la méthode

Sans bruit, le système (2.15) est un système approché du système (2.11). En effet, soit $\epsilon > 0$. Par la méthode du point milieu, chaque coefficient $A_{i,j}^\psi$ est égal à :

$$ \int_{t_i - \epsilon}^{t_i + \epsilon} m_j(s) \psi_i(s) ds = 2\epsilon m_j(t_i) \psi_i(t_i) + \frac{\epsilon^3}{3} (m_j \psi_i)''(\zeta_i), \quad \zeta_i \in [t_i - \epsilon, t_i + \epsilon]. \quad (2.18) $$

Soit $E_{ij} = (m_j \cdot \psi_i)''(\zeta_i)$, $j = 1, \dots, q$ and $F_i = (m_0 \cdot \psi_i)''(\zeta_i)$. Selon l'hypothèse 3, $\psi_i(t_i) = \psi_{\max}$ pour tout $i = 1, \dots, M$ et les relations suivantes sont obtenues :

$$ \frac{A^\psi}{2\epsilon} = \psi_{\max} A + \psi_{\max} \delta A, \quad \frac{b^\psi}{2\epsilon} = \psi_{\max} b + \psi_{\max} \delta b, \quad (2.19) $$

avec

$$ \delta A = \frac{\epsilon^2}{6 \psi_{\max}} E, \quad \delta b = \frac{\epsilon^2}{6 \psi_{\max}} F, \quad (2.20) $$

où $E = (E_{i,j})_{i=1,\dots,q \atop j=1,\dots,q}$ et $F = (F_i)_{i=1,\dots,q}^T$. La substitution de la relation (2.19) dans l'équation (2.15) donne

$$ (A + \delta A) \gamma^\psi = b + \delta b. \quad (2.21) $$

Ainsi, le système (2.15) est un système perturbé du système (2.11) si $\|\delta A\|_\infty$ et $\|\delta b\|_\infty$ sont suffisamment petits.

La proposition suivante donne une estimation de l'erreur relative entre $\gamma$ et $\gamma^\psi$.

**Proposition 5.** Si $\|\delta A\|_\infty$ et $\|\delta b\|_\infty$ sont suffisamment petits, la solution du système (2.15) est une solution approchée du système (2.11).

Soit $K(A)$ le conditionnement de la matrice $A$. Soit $M_2 = \max(\|E\|_\infty, \|F\|_\infty)$.

Si

$$ \frac{\epsilon^2 M_2}{\psi_{\max}} < \frac{6 \|A\|_{\infty}}{K(A)}, \quad (2.22) $$

alors l'erreur relative entre la solution réelle $\gamma$ de (2.11) et la solution approchée $\gamma^\psi$ de (2.15) satisfait :

$$ \frac{\|\gamma^\psi - \gamma\|_\infty}{\|\gamma\|_\infty} \le \frac{1}{\frac{6 \psi_{\max}}{\epsilon^2 M_2} \frac{\|A\|_\infty}{K(A)} - 1} \left( \frac{\|A\|_\infty}{\|b\|_\infty} + 1 \right). \quad (2.23) $$
**Remarque 8.** — L'erreur relative dépend du choix de la fonction test.
— Selon (2.23), plus la longueur du support des fonctions test est petite, meilleure sera l'estimation.

Dans (2.22) et (2.23), $\|E\|_{\infty}$ et $\|F\|_{\infty}$ dépendent de l'inconnue $\zeta_i$ donnée dans l'équation (2.18). Par conséquent, (2.22) et le membre de droite de (2.23) ne peuvent pas être directement calculés. Le corollaire suivant donne un moyen d'obtenir une expression calculable.

**Corollaire 2.** Notons $\tilde{M}_2 = \max_{1 \le i \le q} \sum_{j=0}^q \max_{t \in [0,T]} (|m_j \psi_k)''(t)|)$. Sous l'hypothèse de la proposition 5, si

$$ \frac{\epsilon^2 \tilde{M}_2}{\psi_{\max}} < \frac{6 \|A\|_{\infty}}{K(A)}, \qquad (2.24) $$

l'erreur relative entre la solution réelle $\gamma$ de (2.11) et la solution approchée $\gamma^\psi$ de (2.15) satisfait :

$$ \frac{\|\gamma^\psi - \gamma\|_{\infty}}{\|\gamma\|_{\infty}} \le \frac{1}{\frac{6\psi_{\max}}{\epsilon^2 \tilde{M}_2} \frac{\|A\|_{\infty}}{K(A)} - 1} \left( \frac{\|A\|_{\infty}}{\|b\|_{\infty}} + 1 \right) \qquad (2.25) $$

Dans les applications numériques, les intégrales dans le système (2.15) sont estimées par une méthode d'intégration numérique, dans notre cas, la méthode composite de Newton-Cotes. Dans la partie suivante, l'erreur due à ce choix est donnée.

### 2.4.2.2 Erreur due à la méthode numérique

Dans toute procédure d'intégration numérique se pose la question du pas de temps utilisé pour définir l'intervalle d'intégration.

Les matrices $A$ et $A^\psi$ définies en (2.12) et (2.16) sont construites à partir de la suite $(t_i)_{i=1,...,M}$ qui peut être une sous-suite de $(\tilde{t}_l)_{l=1,...,V}$, $V \ge M$ où $\tilde{t}_1 = 0$, $\tilde{t}_V = T$, $\tilde{t}_l = \tilde{t}_1 + (l-1)\tilde{h}$, $l \ge 2$ et $\tilde{h} = 1/(V-1)$. Ainsi, les intégrales peuvent être calculées avec plus de temps d'intégration que ceux utilisés pour construire les matrices $A$ et $A^\psi$ permettant d'améliorer l'estimation de paramètres. Nous allons donner dans la suite une condition sur $\tilde{h}$ permettant de déduire une borne supérieure à l'erreur relative entre le système (2.15) et le système obtenu en approchant les intégrales par une méthode composite de Newton-Cotes.

Le support de la fonction test $\psi_i$ est supposé être égal à $[\tilde{t}_{\tau-r}, \tilde{t}_{\tau+r}]$, $r > 0$. Utilisant l'interpolation de Lagrange pour approcher $m_j \psi_i$ dans les intégrales du système (2.15), ces dernières peuvent se ré-écrire :

$$ I^{(i,j)} = I_{\nu}^{(i,j)} + \tilde{E}_{\nu}^{i,j} $$

où $I^{(i,j)} = \int_{supp(\psi_i)} (m_j \psi_i)(s) ds$, $I_{\nu}^{(i,j)}$ est une formule de quadrature dans laquelle les fonctions $m_j \psi_i$ ont été remplacées par des polynômes d'interpolation de degré $\nu$, $\nu \ge 0$ sur le sous-intervalle $[\tilde{t}_{\tau-r}, \tilde{t}_{\tau+r}]$. Par conséquent,

$I_{\nu}^{(i,j)}$ a la forme suivante : $I_{\nu}^{(i,j)} = \sum_{k=\tau-r}^{\tau+r} \alpha_k^{(i,j)} m_j(\tilde{t}_k) \psi_i(\tilde{t}_k)$. Finalement, pour toute formule de Newton-Cotes, il existe $\zeta \in [0, T]$ tel que l'erreur $\tilde{E}_{\nu}^{i,j}$ peut-être ré-écrite ([167], Théorème 9.2) :

$$ \text{— si } \nu \text{ est pair et } (m_j \psi_i) \in C^{\nu+2}([0, T]), \quad \tilde{E}_{\nu}^{i,j} = \frac{C_{\nu}}{(\nu+2)!} \tilde{h}^{\nu+3} (m_j \psi_i)^{(\nu+2)}(\zeta) $$

$$ \text{— si } \nu \text{ est impair et } (m_j \psi_i) \in C^{\nu+1}([0, T]), \quad \tilde{E}_{\nu}^{i,j} = \frac{C_{\nu}}{(\nu+1)!} \tilde{h}^{\nu+2} (m_j \psi_i)^{(\nu+1)}(\zeta) $$

avec $C_\nu$ la constante définie par

$$ \left\{ \begin{array}{ll} \displaystyle \int_0^\nu t \pi_{\nu+1}(t) dt & \text{si } \nu \text{ est pair,} \\ \\ \displaystyle \int_0^\nu \pi_{\nu+1}(t) dt & \text{si } \nu \text{ est impair,} \end{array} \right. \quad \text{où } \pi_{\nu+1}(t) = \prod_{i=0}^\nu (t-i). \qquad (2.26) $$

Dans les applications numériques, le système (2.15) est substitué par

$$ A^{\psi,a} \gamma^{\psi,a} = b^{\psi,a} \qquad (2.27) $$

$$ \text{avec } A^{\psi,a} := (I_{\nu}^{(i,j)})_{i=1,...,q}, \quad b^{\psi,a} := (-I_{\nu}^{i,0})_{i=1,...,q}. $$

Notons
✓ $\delta A^{\psi,a} := (\tilde{E}_{\nu}^{i,j})_{i=1,...,q}, \delta b^{\psi,a} := (\tilde{E}_{\nu}^{i,0})_{i=1,...,q},$

✓ $K(A^{\psi})$ le conditionnement de la matrice $A^{\psi}$,

✓ $\beta = \nu + 2$ si $\nu$ is impair et $\beta = \nu + 3$ si $\nu$ est pair,

✓ $M_{\beta} = \max_{1 \le i \le q} \sum_{j=0}^{q} \max_{t \in [0,T]} (|(m_j \psi_k)^{(\beta-1)}(t)|)$.

La proposition suivante donne le lien entre $\gamma^{\psi}$, $\gamma^{\psi,a}$ et $\delta A^{\psi}$, $\delta b^{\psi}$.

**Proposition 6.** Si $\| \delta A^{\psi,a} \|_{\infty}$ et $\| \delta b^{\psi,a} \|_{\infty}$ sont suffisamment petits, une solution de (2.27) est une solution approchée de (2.15).

Si

$$ \tilde{h}^{\beta} < \frac{\| A^{\psi} \|_{\infty} (\beta - 1)!}{K(A^{\psi}) |C_{\nu}| M_{\beta}} \quad (2.28) $$

alors l'erreur relative entre la solution réelle $\gamma^{\psi}$ de (2.15) et la solution approchée $\gamma^{\psi,a}$ de (2.27) vérifie :

$$ \frac{\| \gamma^{\psi} - \gamma^{\psi,a} \|_{\infty}}{\| \gamma^{\psi} \|_{\infty}} \le \frac{1}{\frac{(\beta-1)!! \| A^{\psi} \|_{\infty}}{K(A^{\psi}) |C_{\nu}| \tilde{h}^{\beta} M_{\beta}} - 1} \left( 1 + \frac{\| A^{\psi} \|_{\infty}}{\| b^{\psi} \|_{\infty}} \right) \quad (2.29) $$

**Remarque 9.** 1. Dans le cas de

— la formule des trapèzes, l'inégalité (2.29) est ré-écrite ($\nu = 1$) :

$$ \frac{\| \gamma^{\psi} - \gamma^{\psi,a} \|_{\infty}}{\| \gamma^{\psi} \|_{\infty}} \le \frac{1}{\frac{12 \| A^{\psi} \|_{\infty}}{K(A^{\psi}) \hbar^3 M_3} - 1} \left( 1 + \frac{\| A^{\psi} \|_{\infty}}{\| \delta b^{\psi} \|_{\infty}} \right). \quad (2.30) $$

— la formule de Simpson, l'inégalité (2.29) est ré-écrite ($\nu = 2$) :

$$ \frac{\| \gamma^{\psi} - \gamma^{\psi,a} \|_{\infty}}{\| \gamma^{\psi} \|_{\infty}} \le \frac{1}{\frac{90 \| A^{\psi} \|_{\infty}}{K(A^{\psi}) \hbar^5 M_5} - 1} \left( 1 + \frac{\| A^{\psi} \|_{\infty}}{\| \delta b^{\psi} \|_{\infty}} \right). \quad (2.31) $$

2. Plus le pas de temps $\tilde{h}$ est petit, plus l'erreur relative (2.29) est petite.

### 2.4.2.3 Erreur due au bruit

Dans [81], les auteurs ont mis en évidence deux sortes de perturbations. Les premières comme les perturbations constantes d'amplitude inconnue sont dites structurées et sont solutions d'une équation différentielle linéaire homogène. Elles peuvent être annihilées par des opérateurs différentiels linéaires. C'est le cas, par exemple, d'un bruit aléatoire de moyenne constante et inconnue. Les secondes sont dites non structurées et correspondent aux perturbations hautes fréquences qui peuvent être réduites par des filtres passe bas comme les intégrales itérées. Nous avons montré que la méthode que nous proposons permet d'éliminer une partie du bruit structuré.

Tout bruit $\eta(t)$ peut s'écrire $\eta(t) = \eta_0(t) + \kappa$ où $\eta_0(t)$ est un bruit de moyenne nulle et $\gamma$ est une constante représentant la moyenne de $\eta(t)$. Soit une sortie mesurée s'écrivant $y_{\eta}(t) = y(t) + \eta(t)$ où $y(t)$ est la sortie au temps $t$ sans perturbation. Dans le système (2.15), $y$ est alors substitué par $y_{\eta}$.

Supposons que les intégrales sont de la forme $\int_{supp(\psi_i)} y_{\eta}^{(l)}(s)\psi_i(s)ds$, $l \ge 0$. La forme du bruit dans (2.15) est alors :

$\int_{supp(\psi_i)} \eta(s)\psi_i^{(l)}(s)ds$. Comme $\eta(t) = \eta_0(t) + \kappa$, l'intégrale s'écrit :

$$ \int_{supp(\psi_i)} \eta(s) \psi_i^{(l)}(s) ds = \int_{supp(\psi_i)} \eta_0(s) \psi_i^{(l)}(s) ds + \kappa \int_{supp(\psi_i)} \psi_i^{(l)}(s) ds $$

soit,

$$ \int_{t_0 - \epsilon}^{t_0 + \epsilon} \eta(t) \psi_i^{(l)}(t) dt = \int_{t_0 - \epsilon}^{t_0 + \epsilon} \eta_0(t) \psi_i^{(l)}(t) dt + \gamma \int_{t_0 - \epsilon}^{t_0 + \epsilon} \psi_i^{(l)}(t) dt. $$

Si $l \ge 1$, $\int_{t_0 - \epsilon}^{t_0 + \epsilon} \psi_i^{(l)}(t) = \psi_i^{(l-1)}(t_0 + \epsilon) - \psi_i^{(l-1)}(t_0 - \epsilon)$ est nul dès que $0^+ \notin [t_0 - \epsilon, t_0 + \epsilon]$. Ainsi, les perturbations constantes sont annihilées lorsque $0^+$ n'est pas dans le support de $\psi_i$.

Enfin, le choix de la fonction test peut permettre d'atténuer l'effet du bruit de moyenne nulle. En effet, on a

$$ \| \left( \int_{t_0 - \epsilon}^{t_0 + \epsilon} \eta_0(t) \psi_i^{(l)}(t) dt \right)_{0 \le i \le N} \|_1^2 \le \int_{t_0 - \epsilon}^{t_0 + \epsilon} (\eta_0(t))^2 dt \times \psi_{\epsilon, \varphi, N} \quad (2.32) $$
et

$$E\left(\int_{t_0-\epsilon}^{t_0+\epsilon} (\eta_0(t))^2 dt \times \psi_{\epsilon,\varphi,N}\right) \le \int_{t_0-\epsilon}^{t_0+\epsilon} \text{var}(\eta_0(t)) dt \times \psi_{\epsilon,\varphi,N} \quad (2.33)$$

où var est la variance, E(X) la valeur attendue de la variable aléatoire X et $\psi_{\epsilon,\varphi,N} = \sum_{i=0}^{N} \int_{t_0-\epsilon}^{t_0+\epsilon} |\psi_i^{(l)}(t)|^2 dt$. $\psi_{\epsilon,N,\varphi}$ donne un critère pour choisir une fonction test. Il permet également de déterminer une borne inférieure de la longueur du support des fonctions test permettant d'atténuer l'effet du bruit. Des exemples numériques sont donnés dans [195].

**Remarque 10.** Selon la Remarque 8, un compromis doit être fait entre la longueur du support des fonctions tests : en diminuant la longueur du support l'erreur relative donnée par (2.23) est diminuée et en augmentant la longueur du support, l'effet du bruit est atténué.

Dans ce qui suit, l'erreur relative entre la solution de (2.15) et celle obtenue avec (2.15) lorsqu'il y a du bruit est étudiée.

Quand la sortie est perturbée par le bruit, $m_j(y, u)$ dans (2.15) est substitué par $m_j(y_\eta, u)$ et c'est le système suivant qui est résolu

$$(A^\psi + \delta A^\psi) \tilde{\gamma}^\psi = b^\psi + \delta b^\psi \quad (2.34)$$

$$\text{avec } A^\psi + \delta A^\psi = \left( \int_{t_i-\epsilon}^{t_i+\epsilon} m_j(y_\eta, u)(s) \psi_i(s) ds \right)_{i=1,...,q; j=1,...,q}, \\ b^\psi + \delta b^\psi = \left( - \int_{t_i-\epsilon}^{t_i+\epsilon} m_0(y_\eta, u)(s) \psi_i(s) ds \right)_{i=1,...,q}.$$

Dans le but de mesurer l'erreur due au bruit, nous supposons que $m_j(y_\eta, u)(t)$ s'écrit $m_j(y_\eta, u)(t) = m_j(t) + f_{\eta,j}(t)$ où $f_{\eta,j}$ contient toutes les perturbations.

La proposition suivante donne l'estimation de l'erreur relative entre les solutions de (2.15) et (2.34).

**Proposition 7.** Soit $\gamma^\psi \in \mathbb{R}^q$ la solution de (2.15) et $\tilde{\gamma}^\psi$ la solution de (2.34). Si $\| \delta A^\psi \|_\infty$ et $\| \delta b^\psi \|_\infty$ sont suffisamment petits, une solution de (2.34) est une solution approchée de (2.15).

Supposons que $m_j(y_\eta, u)(t)$ s'écrit $m_j(y_\eta, u)(t) = m_j(t) + f_{\eta,j}(t)$. Si $M_\eta = \max_{1 \le i \le q, 0 \le j \le q} |f_{\eta,j}(t)\psi_i(t)|$ existe et

$$M_\eta < \frac{\| A^\psi \|_\infty}{2q\epsilon K(A^\psi)} \quad (2.35)$$

alors, l'erreur relative entre $\gamma^\psi$ et $\tilde{\gamma}^\psi$ est donnée par :

$$\frac{\| \gamma^\psi - \tilde{\gamma}^\psi \|_\infty}{\| \gamma^\psi \|_\infty} \le \frac{1}{2\epsilon q K(A^\psi) M_\eta} \left( 1 + \frac{\| A^\psi \|_\infty}{\| b^\psi \|_\infty} \right). \quad (2.36)$$

**Remarque 11.** La proposition précédente montre là encore l'importance du choix de la fonction test et de son support. Celui-ci peut-être fait en faisant une étude préalable du système et du nombre de mesures (voir [195]).

### 2.4.2.4 Application au modèle de transmission du chikungunya

Cette méthode d'estimation a été appliquée au modèle décrit à la section 2.3.4. Dans cette partie, nous cherchons à estimer les paramètres $s, s_L, K_E, K_L, \beta_H$ et $\beta_m$. Rappelons que les sorties du modèle sont $y_1 = L, y_2 = S_H$ et $y_3 = I_H$. En supposant $b, b_H, d, d_L, d_m$ connus, les paramètres $s, s_L, \beta, \beta_H, K_E, K_L$ sont identifiables d'après l'algorithme *Identifiability Tree*.

Dans [196], les polynômes ES sont utilisés pour faire une première estimation des paramètres sans connaissance a priori de leur valeur. Il serait donc tentant de reprendre les relations obtenues à la section 2.3.4. Or, l'ordre de dérivation des polynômes ES est égal à 3 rendant difficile l'exploitation de ces relations. En effet, estimer les dérivées à partir de données bruitées est un problème mal posé.

Une autre difficulté cette fois-ci numérique provient de la grande différence de sensibilités des sorties par rapport à certains paramètres. Par exemple, la sensibilité de $L$ par rapport à $s$ et $s_L$ a un ordre de grandeur de 150 par rapport à $K_E$ et $A_0$ de l'ordre de 0.3 (cf figure 2.2). Ces difficultés ont été surmontées en travaillant tout d'abord sur les sous-systèmes en non pas sur le système en entier, en utilisant des "astuces" d'intégration et en considérant la sensibilité des sorties par rapport aux paramètres.

Par la suite nous supposerons que $K_L = K_E/2$ comme dans [146].
<figure><img src="image_7.png" /><figcaption>FIGURE 2.2 – Sensibilité de L par rapport aux paramètres s, $s_L$, $K_E$, $A_0$.</figcaption></figure>

## Simulations

Les simulations sont faites à partir des conditions initiales

$$ (E(0), L(0), A(0), S_H(0), I_H(0), I_m(0)) = (500, 200, 100, 0.9, 0.1, 0.2) $$

et des valeurs de paramètres donnés dans la table 2.1. Les sorties simulées $L$, $S_H$, $I_H$ sont ensuite perturbées par un bruit Gaussien telle que l'erreur relative a une valeur maximale de 5% puis de 10%. Les mesures sont supposées être faites durant 20 jours en $t_i$, $i = 1, \dots, M = 201$.

L'estimation des dérivées est faite en reprenant la méthode proposée dans [195] et qui est basée sur une approche utilisant la théorie des distributions.

Comme des intégrales similaires sont utilisées pour résoudre le système (2.15) et estimer les dérivées, on utilisera les mêmes fonctions test avec la même longueur de support. Les fonctions test choisies $\psi_i$, $i = 1, \dots, M$ sont

des fonctions centrées en $t_i$ et construites à partir de la fonction : $\psi(t) = \begin{cases} e^{\frac{t^2}{2-\epsilon^2}} & \text{if } |t| < \epsilon \\ 0 & \text{if } |t| \ge \epsilon \end{cases}$ dont le support est égal à $[-\epsilon, \epsilon]$. Des translations et des dilatations permettent de définir des fonctions test sur n'importe quel intervalle. L'erreur due à la méthode a été calculée en utilisant la Proposition 5 et le Corollaire 2 et nous a permis de mieux définir la longueur des supports des fonctions tests. Dans les applications, la longueur du support de toutes les fonctions test a été choisi égale à 2.9 et les intégrales ont été estimées par la méthode des trapèzes.

## Estimation des paramètres du modèle (2.8)

Pour estimer les paramètres $s, s_L, K_E, \beta_H, \beta_m$ dans (2.8a) et (2.8b), les deux systèmes couplés peuvent être considérés comme un système unique dans lequel $L$, $S_H$ et $I_H$ sont les sorties. Cependant, l'algorithme de Rosenfeld-Groebner donne des polynômes non utilisables pour définir une procédure numérique. En effet, le polynôme différentiel liant $L$, $s, s_L$ et $K_E$ contient 36 expressions et induit une matrice mal conditionnée de dimension $35 \times M$. Il est possible d'obtenir des relations plus simples en travaillant système par système [146, 227]. Du premier modèle, nous avons obtenu un premier polynôme contenant les paramètres $s, s_L, K_E$. Une fois estimés, nous avons repris l'expression explicite de $A$ dépendant uniquement de $L$ et donnée par l'algorithme de Rosenfeld-Groebner. $A$ et $L$ (donc $A/L$) considérés comme connus ont constitué les entrées du deuxième modèle (2.8b). En utilisant l'algorithme de Rosenfeld-Groebner, un deuxième polynôme $E$ a été
<table><caption>TABLE 2.2 – Valeurs des paramètres pour le modèle du chikungunya</caption><thead><tr><th>Paramètres</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>s</td><td>taux de transfert de E à L</td><td>0.7</td></tr><tr><td>s<sub>L</sub></td><td>taux de transfert de L à A</td><td>0.5</td></tr><tr><td>K<sub>E</sub></td><td>capacité d'accueil des oeufs</td><td>1000</td></tr><tr><td>A<sub>0</sub></td><td>conditions initiales pour A (A(0))</td><td>100</td></tr><tr><td>b</td><td>taux de ponte intrinsèque des femelles</td><td>6</td></tr><tr><td>d</td><td>taux de mortalité des oeufs</td><td>0.2</td></tr><tr><td>d<sub>L</sub></td><td>taux de mortalité des larves</td><td>0.2</td></tr><tr><td>d<sub>m</sub></td><td>taux de mortalité des adultes femelle</td><td>0.25</td></tr><tr><td>b<sub>H</sub></td><td>taux de natalité (humains)</td><td>0.0000457</td></tr><tr><td>β<sub>H</sub></td><td>taux d'infection entre humains et vecteurs</td><td>0.75</td></tr><tr><td>β<sub>m</sub></td><td>taux d'infection entre vecteurs et humains</td><td>0.5</td></tr><tr><td>γ</td><td>taux de transfert des infectés vers les immunisés</td><td>0.1428</td></tr></tbody></table>

trouvé et qui a permis de déterminer $\beta_m$ et $\beta_H$.

Le système (2.8a)

Dans [227], pour obtenir des polynômes ES plus simples, un système équivalent au système (2.8a) a été obtenu en intégrant la troisième équation. Cette opération conduit à poser $w(t) = e^{-dmt}$, $v(t) = e^{-dmt} \int_0^t L(s)e^{dms}ds$ et $A_0$, la condition initiale de A. Le système (2.8a) se ré-écrit :

$$ egin{cases} E'(t) = bA(t)\left(1 - \frac{E(t)}{K_E}\right) - (s+d)E(t) \\ L'(t) = sE(t)\left(1 - \frac{2L(t)}{K_E}\right) - (s_L + d_L)L(t) \\ A(t) = A_0w(t) + s_Lv(t) \\ v'(t) = L(t) - d_mv(t) \\ w'(t) = -d_mw(t) \end{cases} \quad (2.37) $$

Les systèmes (2.8(a)) et (2.37) sont équivalents si on impose les conditions initiales $v(0) = 0$ et $w(0) = 1$. Dans ce cas, l'analyse de stabilité du modèle (2.37) est la même que celle du modèle (2.8).

Comme $L$ est une sortie du modèle, $w$ et $v$ peuvent être estimées. Nous pouvons alors ajouter les équations $y_1 = L$, $y_2 = v$ et $y_3 = w$ au système (2.37). Dans le système (2.37), nous avons introduit la constante $A_0$. Cette constante est considérée comme un nouveau paramètre à estimer d'où le nouveau vecteur de paramètres $(A_0, s, s_L, K_E)$.

En utilisant le package *DifferentialAlgebra* de Maple et avec l'ordre d'élimination $[A_0, K_E, s, s_L] \prec [y_1, y_2, y_3] \prec [A, E, L, v, w]$, nous avons obtenu le polynôme

$$ egin{aligned} P_1 = A_0 K_E^2 b s y_3 + K_E^2 b s s_L y_2 & - K_E (b d_l s_L + 4 b s s_L + b s_L^2) y_2 y_1 \\ & - K_E^2 (d d_l + d s_L + d_l s + s s_L) y_1 - K_E (A_0 b d_l + 4 A_0 b s + A_0 b s_L) y_3 y_1 \\ & - K_E^2 (d + d_l + s + s_L) \dot{y}_1 - K_E b s_L y_2 \dot{y}_1 - A_0 K_E b y_3 \dot{y}_1 - K_E^2 \ddot{y}_1 \\ & + 2 (b d_l s_L + 2 b s s_L + b s_L^2) y_2 y_1^2 + 2 (A_0 b d_l + 2 A_0 b s + A_0 b s_L) y_3 y_1^2 \\ & + 2 K_E (d d_l + d s_L + d_l s + s s_L) y_1^2 + 2 b s_L y_2 \dot{y}_1 y_1 + 2 K_E (d + s) \dot{y}_1 y_1 \\ & + 2 A_0 b y_3 \dot{y}_1 y_1 - 2 K_E (\dot{y}_1^2 - \ddot{y}_1 y_1) = 0, \end{aligned} \quad (2.38) $$

et les deux polynômes

$$ P_A = -A_0 y_3 - s_L y_2 + A = 0, \quad (2.39) $$

$$ P_E = -K_E s E + (K_E d_l + K_E s_L + 2sE) y_1 + K_E \dot{y}_1 = 0. \quad (2.40) $$
$P_A$ et $P_E$ permettent d'estimer $A$ et $E$ respectivement si les paramètres $A_0, K_E, s, s_L$ sont connus.

Pour obtenir un polynôme de la forme (2.4), $P_1$ est divisé par $K_E y_1$. Nous avons ainsi travaillé avec le polynôme

$$ \tilde{P}_1 = \frac{P_1}{y_1 K_E}. \text{ Soit le polynôme } P_1^\psi = \int_{\text{supp}(\psi)} \tilde{P}_1(s) \psi(s) ds. $$

**Remarque 12.** 1. Divisé $P_1$ par $y_1$ permet tout d'abord de reporter la dérivée seconde dans l'expression

$$ \int_{\text{supp}(\psi)} \left( \frac{\dot{y}_1(s)^2}{y_1(s)} \psi(s) - \ddot{y}_1(s) \psi(s) \right) ds \text{ de } P_1^\psi \text{ sur la fonction test } \psi. \text{ Ce terme correspond à } b^\psi, \text{ le second membre de (2.15). Ensuite, la division diminue de façon importante le conditionnement de la matrice } A^\psi. $$

2. Les sortie bruitées $y_2$ et $y_3$ ne dépendant que de $y_1$ sont estimées à partir de l'intégration des deux dernières équations de (2.37) en utilisant l'algorithme du point milieu.

3. La méthode que nous avons développée dans [196] consistant à intégrer une fois le polynôme entrée-sortie $\tilde{P}_1$ n'a pas donné de résultat.

À partir des égalités

$$ \left\{ \begin{array}{l} \dot{y}_2 = y_1 - d_m y_2, \quad \dot{y}_3 = -d_m y_3, \\ \\ \displaystyle \frac{\ddot{y}_1}{y_1} = \left(\frac{\dot{y}_1}{y_1}\right) + \frac{\dot{y}_1^2}{y_1^2}, \end{array} \right. \qquad (2.41) $$

nous en avons déduit les opérations suivantes qui ont été appliquées à $P_1^\psi$ pour décroître l'ordre de dérivation de $y_1$ :

$$ \left\{ \begin{array}{l} \displaystyle \int \dot{y}_1 y_3 \psi = [y_1 y_3 \psi] + \int y_1 y_3 (d_m \psi - \dot{\psi}), \\ \\ \displaystyle \int \dot{y}_1 \psi = [y_1 \psi] - \int y_1 \dot{\psi}, \quad \int \ddot{y}_1 \psi = [\dot{y}_1 \psi - \dot{\psi} y_1] + \int y_1 \ddot{\psi}, \\ \\ \displaystyle \int \dot{y}_1 y_2 \psi = [y_1 y_2 \psi] - \int (y_1(y_1 - d_m y_2)\psi + y_1 y_2 \dot{\psi}), \\ \\ \displaystyle \int \frac{\ddot{y}_1}{y_1} \psi = [\frac{\dot{y}_1}{y_1} \psi] - \int \left(\frac{\dot{y}_1}{y_1} \dot{\psi} - \frac{\dot{y}_1^2}{y_1^2} \psi\right). \end{array} \right. \qquad (2.42) $$

Ces opérations nous ont permis d'obtenir un nouveau polynôme $\tilde{P}_1^\psi$ et d'en déduire un système linéaire de la forme (2.15) avec $A^\psi$ de dimension $M \times 15$. En effet, $\tilde{P}_1^\psi$ contient 16 expressions dont l'une d'elle ne dépend pas des paramètres à estimer. Cependant le conditionnement de la matrice $A^\psi$ est de l'ordre de $10^{12}$ et la factorisation QR donne pour certains paramètres des valeurs négatives.

En reprenant les courbes de sensibilité (Fig. 2.2), nous pouvons remarquer que les sorties sont plus sensibles aux paramètres $s, s_L$ qu'aux paramètres $K_E$ et $A_0$.

Comme le système est très peu sensible aux paramètres $K_E$ et $A_0$, la stratégie consiste à se focaliser sur l'estimation de $s, s_L$. En supposant que les deux paramètres $K_E$ et $A_0$ sont connus, le polynôme $P_1^\psi$ se ré-écrit :

$$ \tilde{m}_0(y_1, y_2, y_3) + \tilde{m}_1(y_1, y_2, y_3)s + \tilde{m}_2(y_1, y_2, y_3)s_L + \tilde{m}_3(y_1, y_2, y_3)s s_L + \tilde{m}_4(y_1, y_2, y_3)s_L^2 = 0, \qquad (2.43) $$

$\tilde{m}_i$ dépendant de $K_E$ et $A_0$. Dans ce cas, $\theta = (s, s_L, s_L, s_L)^T$ et $A^\psi \in \mathbb{R}^{M \times 4}$.

L'avantage de cette ré-écriture est que la matrice $A^\psi$ est réduite à quatre colonnes et son conditionnement est de l'ordre de 100.

Pour les deux autres paramètres, nous supposons qu'une estimation de $K_E$ et $A_0$ peut-être obtenue expérimentalement avec une erreur de 10% et 50% respectivement. Toutefois, ces estimés doivent être améliorées pour estimer les paramètres $\beta_H, \beta_m$ du système (2.8(b)). En effet, l'expression (2.39) qui permet de calculer $A$ et par conséquent $\beta_H, \beta_m$ contient la constante $A_0$.

$K_E$ et $A_0$ variant respectivement de 900 à 1100 et de 50 à 150 avec un pas égal à 10, nous estimons $s$ et $s_L$ pour chaque couple de valeurs $(K_E, A_0)$. L'erreur entre la courbe $L$ reconstruite à partir du système (2.8(a)) avec la courbe de mesures permet de garder les meilleures valeurs estimées. Cette procédure est résumée dans l'algorithme suivant :

## Algorithme 1.

* /* Entrées : Le vecteur de mesures $y_\eta$

* Pour $K_E$ variant de 900 à 1100 avec un pas de 10

* Pour $A_0$ variant de 50 à 150 avec un pas de 10

* Estimation de $s, s_L$

* Estimation de la sortie reconstruite $y_{est}$ avec les paramètres $A_0, K_E, s, s_L$
Estimation de l'erreur relative entre y<sub>η</sub> et y<sub>est</sub>.

fin

fin

Avec cette procédure, nous avons obtenus les estimées suivantes (n<sub>0</sub> = 1, n<sub>1</sub> = 150 dans (2.15)) :

— Pour un bruit de 5% :

— Méthode des distributions : A<sub>0</sub><sup>est</sup> = 110, K<sub>E</sub><sup>est</sup> = 1000, s<sup>est</sup> = 0.6277, s<sub>L</sub><sup>est</sup> = 0.2956

— Algorithme de Levenberg-Marquardt pour améliorer l'estimation de s et s<sub>L</sub> : s = 0.7212, s<sub>L</sub> = 0.5226

— Erreur relative (er) entre la courbe de mesures L et la courbe reconstruite à partir de (A<sub>0</sub>, K<sub>E</sub>, s, s<sub>L</sub>) = (110, 1000, 0.7212, 0.5226) : er = 0.016.

— Pour un bruit de 10% :

— Méthode des distributions : A<sub>0</sub><sup>est</sup> = 110, K<sub>E</sub><sup>est</sup> = 1000, s<sup>est</sup> = 0.669, s<sub>L</sub><sup>est</sup> = 0.251

— Algorithme de Levenberg-Marquardt pour améliorer l'estimation de s et s<sub>L</sub> : s = 0.7133, s<sub>L</sub> = 0.5169

— Erreur relative (er) entre la courbe de mesures L et la courbe reconstruite à partir de (A<sub>0</sub>, K<sub>E</sub>, s, s<sub>L</sub>) = (110, 1000, 0.7133, 0.5169) : er = 0.0364.

Le système (2.8b)

La troisième équation de (2.8(b)) lie la population humaine à la population du moustique avec le terme L(t)/A(t). Les courbes de sensibilité par rapport aux paramètres β<sub>H</sub> et β<sub>m</sub> sont données à la Figure 2.3.
Selon la section précédente, A(t) et L(t) peuvent être évaluées à chaque pas de temps t<sub>k</sub>, k = 1, ..., M. Ainsi, nous allons considérer par la suite que L/A constitue une entrée connue u du système (2.8(b)). En procédant comme précédemment, en ajoutant y<sub>4</sub> = I<sub>H</sub>, y<sub>5</sub> = S<sub>H</sub> au système (2.8(b)) et en considérant l'ordre d'élimination [β<sub>H</sub>, β<sub>m</sub>] ≺ [y<sub>4</sub>, y<sub>5</sub>, u] ≺ [I<sub>H</sub>, S<sub>H</sub>, I<sub>m</sub>], le polynôme ES suivant a été obtenu :

$$
\begin{equation}
\begin{split}
P_2 = {}& -b_H s_L u y_5 + b_H \dot{y}_5 + b_H s_L u y_5^2 + s_L u \dot{y}_5 y_5 - \dot{y}_5^2 + \dot{y}_5 y_5 \\
& + (b_H \beta_m + \beta_H \beta_m) y_4 y_5^2 + \beta_m (y_4 \dot{y}_5 y_5 - b_H y_4 y_5).
\end{split}
\tag{2.44}
\end{equation}
$$

Le polynôme $P_2^\psi = \int_{supp(\psi)} \frac{P_2(s)}{y_5(s)} \psi(s)ds$ est utilisé pour estimer les paramètres $\beta_H$ et $\beta_m$ après avoir utilisé les opérations suivantes :

$$
\left\{
\begin{array}{l}
\displaystyle \int \dot{y}_5 \psi = [y_5 \psi] - \int y_5 \dot{\psi}, \\[1em]
\displaystyle \int \ddot{y}_5 \psi = [\dot{y}_5 \psi - y_5 \dot{\psi}] + \int y_5 \ddot{\psi}.
\end{array}
\right.
\qquad (2.45)
$$

Avec l'approche des distributions et les valeurs de A₀, s_L calculées précédemment, la valeur des paramètres
trouvée en résolvant le système (2.15) est (n₀ = 1, n₁ = 169) :

— Pour un bruit de 5% :

— Méthode des distributions : $\beta_H^{est} = 0.8120$ et $\beta_m^{est} = 0.5135460$

— Algorithme de Levenberg-Marquardt : $\beta_H = 0.7529$ et $\beta_m = 0.4971$

— Erreur relative entre les sorties et les courbes reconstruites avec (β<sub>H</sub>, β<sub>m</sub>) = (0.7529, 0.4971) : er = 0.034.

— Pour un bruit de 10% :

— Méthode des distributions : $\beta_H^{est} = 0.2235$ et $\beta_m^{est} = 0.2215$

— Algorithme de Levenberg-Marquardt : $\beta_H = 0.7494$ et $\beta_m = 0.5021$

— Erreur relative entre les sorties et les courbes reconstruites avec (β<sub>H</sub>, β<sub>m</sub>) = (0.7494, 0.5021) : er = 0.0643.

2.5 Identifiabilité ensembliste

Le travail présenté dans cette partie provient des articles [100, 101] écrits dans le cadre du projet ANR MagicSPS. L'un des objectifs du projet était de développer des méthodes d'estimation de paramètres dans un contexte d'erreurs bornés et caractériser l'ensemble de toutes les valeurs des paramètres qui sont consistant avec les données mesurées. En effet, les systèmes complexes sont soumis à des perturbations qui compliquent la modélisation du système étudié. Dans le cadre stochastique, l'incertitude est prise en compte en ajoutant des hypothèses appropriées sur le bruit et les distributions de probabilité des erreurs de modèle. Or la caractérisation du bruit et des perturbations peut être difficile à déterminer dans certains cas comme les tolérances de fabrication sur les valeurs de paramètres. Cette contrainte peut rendre difficile l'utilisation des méthodes probabilistes [7, 211].
<figure><img src="image_8.png" /><figcaption>FIGURE 2.3 – Courbes de sensibilité de $S_H$, $I_H$ par rapport aux paramètres $\beta_H$, $\beta_m$.</figcaption></figure>

Dans le cadre du projet, nous avons considéré une approche alternative en supposant que les valeurs incertaines des variables et des paramètres appartenaient à des ensembles bornés nous amenant à étendre la définition classique d'identifiabilité. Nous avons ensuite exploité les polynômes ES dans le cadre ensembliste et développé des méthodes certifiées d'estimation de paramètres, de détection et d'isolation de fautes. Lorsqu'on parle de méthode certifiée, on sous-entend que les ensembles solutions calculés contiennent l'ensemble des vecteurs de paramètres consistant avec les mesures [107]. Dans ce travail, nous parlerons d'ensembles d'appartenance, en abrégé SM pour Set-Membership en anglais. Ceux-ci peuvent avoir des formes différentes comme des ellipsoïdes [114], des boîtes [108], des parallélotopes [40], des zonotopes [4] ou d'autres polytopes.
L'exemple suivant illustre les questions que nous nous sommes posées.

**Exemple 2.5.1.** Soit l'équation différentielle $\dot{x}(t) = x(t) + t \cos(p)$ avec les conditions initiales $x(0) = 1$ associées. La solution est $x(t) = (1 - t - e^t) \cos(p)$. Supposons que le paramètre $p$ est à incertitude bornée et que son ensemble admissible est l'ensemble $U_p = [0, 2\pi]$. Nous pouvons remarquer qu'il existe deux valeurs de paramètres dans $P^* := [\frac{\pi}{2}, \frac{3\pi}{2}]$ qui vont générer les mêmes trajectoires (figure 2.4) mais que les paramètres dans $P^*$ et ceux de son complémentaire ne générerons pas les mêmes trajectoires (figure 2.5).

Cet exemple amène à se poser les deux questions suivantes : étant donné un intervalle borné solution trouvé par un algorithme d'optimisation,

1. comment peut-on être sûr que les valeurs de paramètres dans cet intervalle et dans son complémentaire vont générer des trajectoires différentes ?

2. comment peut-on distinguer les intervalles solutions sur lesquels il y a identifiabilité et ceux pour lesquels cette propriété n'est pas vérifiée ?

Sur la base de l'exemple suivant, nous avons étendu la notion d'identifiabilité classique à : i) la SM-identifiabilité, la $\mu$-SM-identifiabilité, et $\varepsilon$-SM-identifiabilité.

Nous nous sommes ensuite intéressées à la façon dont ces propriétés impactent le problème d'estimation des paramètres dans le cadre ensembliste. Ce problème est caractérisé par deux nouvelles propriétés. La première garantie que l'ensemble faisable de paramètre (FPS) est réduit à un seul ensemble connexe borné qui sera dit robuste (Soundness en anglais). La seconde est une propriété numérique que nous avons appelée la $\varepsilon$-consistance (ou $\varepsilon$-consistence en anglais). Celle-ci garantit que le FPS est composé d'un même nombre d'ensembles connexes disjoints que l'ensemble solution retourné par un algorithme d'estimation de paramètres de précision $\varepsilon$. Malgré
<figure><img src="image_9.png" /><figcaption>FIGURE 2.4 – Représentation de la fonction $w \mapsto \cos(w)$ sur l’intervalle $[0, 2\pi]$. Dans l’intervalle $P^* := [\frac{\pi}{2}, \frac{3\pi}{2}]$, deux valeurs de $w$ donnent la même valeur de $\cos(w)$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 2.5 – Représentation de la fonction $w \mapsto \cos(w)$ sur l’intervalle $[0, 2\pi]$. Pour deux valeurs de $w$, l’une dans $P^* := [\frac{\pi}{2}, \frac{3\pi}{2}]$, l’autre dans son complémentaire, les valeurs de $\cos(w)$ sont distinctes.</figcaption></figure>
une abondante littérature sur l'identification dans le cadre ensembliste [104, 108, 169, 208], nous n'avions pas trouvé de travaux sur ces relations.

Par la suite, nous supposons que les incertitudes du modèle portent sur les conditions initiales et les paramètres. Ainsi, nous supposons que $x_0 \in X_0$ où $X_0$ est un ensemble borné et $p \in \mathcal{P}$ est un ensemble connexe de $\mathcal{U}_P$ où $\mathcal{U}_P \subseteq \mathbb{R}^q$ est un ensemble a priori connu de paramètres admissibles.

Le modèle (2.46) se ré-écrit

$$ egin{cases} 
\dot{x}(t,p) = f(x(t,p),u(t),p), \\
y(t,p) = h(x(t,p),p), \\
x(t_0,p) = x_0 \in X_0, \\
p \in P \subset \mathcal{U}_P, t_0 \le t \le T,
\end{cases} 
\quad (2.46) $$

où $x(t,p) \in \mathbb{R}^n$ et $y(t,p) \in \mathbb{R}^m$.

Dans la suite, nous noterons $Y(P)$ l'ensemble des trajectoires solutions de (2.46) générées par les paramètres $p \in P$. $P^c$ désignera le complémentaire de $P$ dans $\mathcal{U}_P$.

## 2.5.1 Concepts utiles

Considérons un ensemble connexe non vide $\Pi$ de $\mathbb{R}^\alpha$, $\alpha \in \mathbb{N}^* \cup \{\infty\}$, $\|\cdot\|$ une norme de $\mathbb{R}^\alpha$ et $d$ la distance associée.

La distance<fnref n="6" /> entre les deux ensembles $\Pi_1$ et $\Pi_2$ de $\mathbb{R}^p$ est définie par :

$$ d(\Pi_1, \Pi_2) = \min_{\pi_1 \in \Pi_1, \pi_2 \in \Pi_2} d(\pi_1, \pi_2). $$

Définissons $\delta(\Pi)$ comme le diamètre de $\Pi$. $\delta(\Pi)$ est défini comme la plus petite borne supérieure de $\{d(\pi_1, \pi_2), \pi_1, \pi_2 \in \Pi\}$. Si $\Pi$ n'est pas borné, nous définissons la distance comme $\delta(\Pi) = +\infty$ [23]. Sur l'espace métrique $(\Pi, d)$, une contraction $\mu$ est une fonction continue de $\Pi$ dans $\Pi$ telle qu'il existe un nombre strictement positif $k < 1$ et que pour tout $\Pi_1, \Pi_2 \subseteq \Pi$, $d(\mu(\Pi_1), \mu(\Pi_2)) < kd(\Pi_1, \Pi_2)$. Dans la suite, $\|\cdot\|$ désignera la norme Euclidienne, $\|\cdot\|_\infty$ la norme infinie, et $\|\cdot\|_1$ la norme 1. Elles seront définies sur $\mathbb{R}^\alpha$, où $\alpha \in \{n, m, q\}$, selon le cas considéré.

## 2.5.2 Définitions

Les définitions suivantes présentent les concepts d'identifiabilité dans le cadre ensembliste.

**Définition 6.** Étant donné le système contrôlé (2.46), considérons un ensemble connexe non vide $P^* \subseteq \mathcal{U}_P$. $P^*$ est globalement SM-identifiable si il existe une entrée $u$ telle que $Y(P^*) \neq \emptyset$ et pour tout ensemble $\bar{P} \subseteq \mathcal{U}_P$, on a $Y(P^*) \cap Y(\bar{P}) \neq \emptyset \implies P^* \cap \bar{P} \neq \emptyset$.

Dans le cas d'un modèle non contrôlé, la définition est la suivante.

**Définition 7.** Étant donné le système (2.46) avec $u \equiv 0$, considérons un ensemble connexe non vide $P^* \subseteq \mathcal{U}_P$. $P^*$ est globalement SM-identifiable si $Y(P^*) \neq \emptyset$ et pour tout ensemble $\bar{P} \subseteq \mathcal{U}_P$, on a $Y(P^*) \cap Y(\bar{P}) \neq \emptyset \implies P^* \cap \bar{P} \neq \emptyset$.

Ainsi, on dira qu'un ensemble connexe $P^*$ est globalement SM-identifiable si les ensembles de sorties générées par le système (2.46) à partir de $P^*$ et tout autre ensemble connexe dans le complémentaire de $P^*$ n'ont pas de trajectoire en commun.

Il est clair que dans l'exemple 2.5.1, le système n'est pas globalement identifiable. Par exemple, toute paire $(p_1 = \pi - \alpha, p_2 = \pi + \alpha)$, avec $\alpha \in [0, \pi]$, produit la même trajectoire car $\cos(\pi - \alpha) = \cos(\pi + \alpha)$ pour $\alpha \in [0, \pi]$. Cependant, les trajectoires provenant de tout ensemble $P^* = [\pi - \alpha, \pi + \alpha]$ avec $\alpha \in [0, \pi]$ sont différentes de toutes trajectoires provenant de son complémentaire. $P^*$ est donc globalement SM-identifiable.

Pour assurer que l'ensemble $P^*$ peut-être contracté tout en gardant la propriété d'identifiabilité ensembliste, nous avons introduit la définition de $\mu$-SM-identifiabilité. Pour cela, une contraction $\mu$ est appliquée à $P^*$ et, d'après le théorème du point fixe, le diamètre de $\mu^n(P^*)$ tend vers zéro quand $n$ tend vers l'infini [148]. La définition est la suivante :

**Définition 8.** Un ensemble connexe non vide $P^* \subseteq \mathcal{U}_P$ est globalement $\mu$-SM-identifiable si $\mu(P^*)$ est globalement SM-identifiable pour toute contraction $\mu$ de $P^*$ à $P^*$.

<footnote n="6">6. Nous faisons délibérément un abus de langage en appelant $d(\Pi_1, \Pi_2)$ une distance entre les deux ensembles $\Pi_1$ et $\Pi_2$ de $\mathbb{R}^p$ même si certaines hypothèses comme l'inégalité triangulaire ne sont pas vérifiées.</footnote>
La proposition suivante découle directement des définitions :

**Proposition 8.** Si l'ensemble connexe non vide $P^* \subseteq \mathcal{U}_P$ est globalement $\mu$-SM-identifiable alors il est globalement SM-identifiable. La réciproque n'est pas vraie.

Le système (2.5.1) donne un contre-exemple de la réciproque. En effet, l'ensemble $P^* = [\pi - \alpha, \pi + \alpha]$ avec $\alpha \in [0, \pi]$ n'est pas $\mu$-SM-identifiable. En supposant $\alpha_1, \alpha_2 \in ]0, \pi[, \alpha_1 < \alpha_2$, tout ensemble $P_1^* = [\pi - \alpha_1, \pi - \alpha_2] \subseteq P^*$ a des trajectoires en commun avec son ensemble complémentaire $P_1^{*c}$.

Si le diamètre de $\mu(P^*)$, $\delta(\mu(P^*))$, ne peut pas être plus petit que $\varepsilon$ sans perdre la SM-identifiabilité, nous nous référerons à la $\varepsilon$-SM-identifiabilité [100].

**Définition 9.** *Considérons un ensemble $P^* \subseteq \mathcal{U}_P$ connexe, non vide et SM-identifiable. Alors $P^*$ est globalement $\varepsilon$-SM-identifiable si il existe un ensemble de contraction $\mu$ de $P^*$ dans $P^*$ tel que*

$$1. \delta(\mu(P^*)) = \varepsilon$$

2. $\mu(P^*)$ est globalement SM-identifiable

3. pour tout $\tilde{\mu}$ tel que $\tilde{\mu}(P^*) \subseteq \mu(P^*)$, $\tilde{\mu}(P^*)$ n'est pas globalement SM-identifiable.

Pour résumer, interpréter l'identifiabilité dans le cadre ensembliste nous a amené à proposer trois définitions. La première considère l'ensemble solution comme une partie agrégée (SM-identifiabilité) contrairement à la seconde qui prend en compte les propriétés de ses sous-ensembles ($\mu$-SM-identifiabilité). Ainsi, la $\mu$-SM-identifiabilité peut être vue comme une généralisation de l'identifiabilité classique : si $P^*$ est $\mu$-SM-identifiable, alors tout $p \in P^*$ est identifiable dans le sens classique [125]. La troisième définition, la $\varepsilon$-SM-identifiabilité peut être vue comme une propriété structurelle de la $\mu$-SM-identifiabilité puisque les sous ensembles SM-identifiables de diamètre $\varepsilon$ et non $\mu$-SM-identifiables sont acceptés.

Ces définitions peuvent-être étendues à des versions structurelles locales.

### 2.5.3 Liens entre la SM-identifiabilité et d'autres concepts

Les liens entre la ($\mu$)-SM-identifiabilité, l'identifiabilité classique au sens de la définition 1 et l'identifiabilité par intervalle ont été donnés dans [101, 102]. Dans cette section, nous allons nous intéresser au lien entre la ($\mu$)-SM-identifiabilité, la $\varepsilon$-identifiabilité [25] et l'injectivité partielle [117]. Ces liens nous ont permis de développer avec Laleh Hosseini une méthode pour partitionner un ensemble admissible de paramètres en sous-ensembles SM-identifiables et $\mu$-SM-identifiable.

#### 2.5.3.1 Liens avec la $\varepsilon$-identifiabilité globale

L'identifiabilité globale dans $P^* \subseteq \mathcal{U}_P$ (g.i.i.$P^*$) a été proposée par [25] comme un moyen de fournir une conclusion plus forte que l'identifiabilité structurelle. Elle garantit que des régions atypiques dans l'espace des paramètres pour lesquelles il n'y a pas identifiabilité n'existent pas. La définition est la suivante :

**Définition 10.** *Étant donné $(u, x_0) \in \mathbb{R}^T \times X_0$, le paramètre $p_i$ est globalement identifiable dans $P^*$ (g.i.i.$P^*$) si*

$$\forall (p, \bar{p}) \in P^*, y(., p) \equiv y(., \bar{p}) \Rightarrow p_i = \bar{p}_i, \qquad (2.47)$$

et le vecteur de paramètre $p$ est g.i.i.$P^*$ si toutes ses composantes sont g.i.i.$P^*$.

L'originalité du travail de [25] est de proposer un moyen pratique de formuler la condition de la Définition 10. La condition peut en effet se ré-écrire :

$$\mathfrak{H}(p, \bar{p}) \in P^* \times P^* \text{ tel que } y(., p) \equiv y(., \bar{p}), \quad \| \bar{p} - p \|_{\infty} > 0. \qquad (2.48)$$

Cette condition se ramène à un problème de satisfaction de contraintes (constraint satisfaction problem ou CSP en anglais) qui peut-être résolu de façon garantie par un algorithme de propagation de contraintes d'intervalle (interval constraint propagation ou ICP en anglais). En pratique, résoudre (2.48) revient à résoudre le problème suivant :

$$\mathfrak{H}(p, \bar{p}) \in P^* \times P^* \text{ tel que } y(., p) \equiv y(., \bar{p}), \quad \| \bar{p} - p \|_{\infty} > \varepsilon, \qquad (2.49)$$

qui définit la $\varepsilon$-g.i.i.$P^*$. Remarquons que cette dernière définition ne donne aucune information sur les sous-ensembles de $P^*$. Le lien entre cette définition et la $\mu$-SM-identifiable est donné dans la proposition suivante.

**Proposition 9.** 1. $P^*$ est $\mu$-SM-identifiable en $P^*$ (dans le sens que $\mathcal{U}_P$ est réduit à $P^*$) si et seulement si (2.48) est vérifié.

2. Si $P^*$ est $\varepsilon$-SM-identifiable en $P^*$, alors la condition (2.49) est vérifiée.

Pour le deuxième point, la réciproque est fausse car quand la condition (2.49) est vérifiée, cela ne fournit aucune information sur les sous-ensembles de $P^*$ de diamètre inférieur à $\varepsilon$.

La condition (2.49) ne permet pas de décider entre la $\mu$-SM-identifiabilité et la $\varepsilon$-SM-identifiabilité. Elle permet toutefois de considérer la $\mu$-SM-identifiabilité en acceptant une précision numérique de $\varepsilon$.
### 2.5.3.2 Liens avec l'injectivité partielle

La définition de l'injectivité partielle a été introduite dans [117]. Cette notion caractérise parfaitement la $\mu$-SM-identifiabilité.

**Définition 11.** Soit une fonction $f: A \to B$ et un ensemble $A_1 \subseteq A$. La fonction $f$ est une injection partielle de $A_1$ sur $A$, ou $(A_1, A)$-injective, si pour tout $a_1 \in A_1$ et tout $a \in A$,

$$a_1 \neq a \Rightarrow f(a_1) \neq f(a).$$

f est dite A-injective si elle est $(A, A)$-injective.

**Exemple 2.5.2.** D'après la figure 2.4, la fonction cosinus est $([0, \frac{\pi}{2}, [0, \frac{3\pi}{2}])$-injective puisque la fonction cosine est injective sur $[0, \frac{\pi}{2}]$ et l'image de $[0, \frac{\pi}{2}]$ par la fonction cosine a une intersection vide avec l'image de $[\frac{\pi}{2}, \frac{3\pi}{2}]$. Par contre, la fonction cosine n'est pas $([\frac{\pi}{2}, \frac{3\pi}{2}, [0, \frac{3\pi}{2}])$-injective car la fonction cosine n'est pas injective dans $[\frac{\pi}{2}, \frac{3\pi}{2}]$.

[117] a proposé un algorithme basé sur l'analyse par intervalle pour tester l'injectivité d'une fonction différentiable. Les auteurs ont également développé un algorithme appelé ITVIA (Injectivity Test Via Interval Analysis) implémenté en C++ qui, pour une fonction donnée, partitionne un ensemble en deux domaines : un domaine sur lequel la fonction est partiellement injective correspondant à une liste d'ensembles $\mu$-SM identifiables et un domaine indéterminé dans lequel la fonction peut être ou non injective. Notons que ITVIA ne donne pas d'information sur ces domaines indéterminés.

Pour caractériser la SM-identifiabilité, nous avons introduit la notion d'injectivité partielle restreinte (restricted-partial injectivity en anglais).

**Définition 12.** Soit une fonction $f: A \to B$ et un ensemble $A_1 \subseteq A$. La fonction $f$ est dite partiellement injective restreinte de $A_1$ à $A$, ou $(A_1, A)$-R-injective, si :

$$\forall a_1 \in A_1, \forall a \in A_1^c, f(a_1) \neq f(a).$$

**Exemple 2.5.3.** D'après la figure 2.4, la fonction cosinus est $([\frac{\pi}{2}, \frac{3\pi}{2}, [0, \frac{3\pi}{2}])$-R-injective car l'image de $[\frac{\pi}{2}, \frac{3\pi}{2}]$ par la fonction cosine a une intersection vide avec l'image de $[0, \frac{\pi}{2}]$. Comme on peut le voir sur cet exemple, il n'est plus imposé à la fonction cosine d'être injective sur $A_1 = [\frac{\pi}{2}, \frac{3\pi}{2}]$.

Tester numériquement la $(P^*, U_P)$-R-injectivité peut se faire en adaptant l'algorithme ITVIA proposé dans [117]. Toutefois, il ne permet pas de décider entre la $\mu$-SM-identifiabilité et la $\varepsilon$-SM-identifiabilité ou la SM-identifiabilité et la $\varepsilon$-SM-identifiabilité.

Dans la proposition suivante, l'injectivité partielle et l'injectivité partielle restreinte sont interprétées en terme de trajectoires ce qui permettra par la suite de faire le lien direct avec les définitions de SM-identifiabilité et $\mu$-SM-identifiabilité. Considérons l'ensemble des sorties $S_u$ résultant de $U_P$ pour une entrée donnée $u$.

**Proposition 10.** Étant donné le modèle (2.46), $P^*$ est SM-identifiable (resp. $\mu$-SM-identifiable) pour une entrée $u$ si et seulement si la fonction $\varphi: U_P \to S_u: p \to y(.,p)$ est $(P^*, U_P)$-R-injective (resp. $(P^*, U_P)$-injective).

Des propositions 9 et 10, nous en avons déduit le corollaire suivant :

**Corollaire 3.** Les propriétés suivantes sont équivalentes :

- $P^*$ est $\mu$-SM-identifiable,

- la fonction $\varphi: U_P \to S_u: p \to y(.,p)$ est $(P^*, U_P)$-injective,

- La condition (2.48) est satisfaite.

Enfin, nous avons l'implication :

**Corollaire 4.** $P^*$ est $\varepsilon$-SM-identifiable implique que $\varphi$ est $(\tilde{P}, U_P)$-R-injective, avec $\tilde{P} \subseteq P^*$ et $\delta(\tilde{P}) \geq \varepsilon$. La réciproque est fausse.
<figure><img src="image_11.png" /><figcaption>FIGURE 2.6 – Exemple de partitionnement d'un domaine en une liste de sous ensembles (μ)-SM-identifiables.</figcaption></figure>

### 2.5.3.3 Partitionnement d'un domaine en sous-ensembles (μ-)SM-identifiables

La proposition 10 montre qu'il est suffisant de trouver une fonction $\varphi : \mathcal{U}_P \to S_u : p \to y(.,p)$ définie de l'ensemble des paramètres à l'ensemble des trajectoires et ayant des propriétés d'injectivité pour vérifier la SM ou la μ-SM-identifiabilité.

La fonction $\varphi$ définie dans la proposition 1 est une candidate naturelle pour la fonction $\varphi$. La Proposition 1 peut ainsi être reformulée de la façon suivante.

**Proposition 11.** [101] Supposons que le déterminant fonctionnel $\Delta R(y, u)$ n'est pas identiquement nul. Soit $P^*$ un sous-ensemble connexe de $\mathcal{U}_P$. Si la fonction $\varphi : p = (p_1, \dots, p_q) \mapsto (\gamma_1(p), \dots, \gamma_d(p), y(t_0^+, p), \dots, y^{(l-1)}(t_0^+, p))$ est $(P^*, \mathcal{U}_P)$-R-injective (resp. $(P^*, \mathcal{U}_P)$-injective) alors $P^*$ est globalement SM-identifiable (resp. μ-SM-identifiable). De plus, si pour une contraction μ, $\mu(P^*)$ a un diamètre égal à $\varepsilon$ et $\varphi$ est $(\mu(P^*), \mathcal{U}_P)$-R-injective mais pas $(\mu(P^*), \mathcal{U}_P)$-injective alors $P^*$ est $\varepsilon$-SM-identifiable. Dans les deux cas, si le coefficient de $y^{(l)}$ in (2.4) n'est pas égal à $t_0$, alors la réciproque est vraie$^7$.

Les domaines sur lesquels la fonction $\varphi$ est partiellement injective restreinte mais pas partiellement injective seront les ensembles SM-identifiables (cf. corollaire 3) et les ensembles sur lesquels la fonction $\varphi$ est partiellement injective seront μ-SM-identifiables. On peut résumer les résultats comme suit :

<table><thead><tr><th>En terme d'identifiabilité</th><th></th><th>En terme d'injectivité</th></tr></thead><tbody><tr><td>P μ-SM identifiable</td><td>$\leftrightarrow$</td><td>$\phi$ partiellement injective de P sur $\mathcal{U}_P$,</td></tr><tr><td>P SM-identifiable</td><td>$\leftrightarrow$</td><td>$\phi$ partiellement injective restreinte de P sur $\mathcal{U}_P$,</td></tr><tr><td>P $\varepsilon$-SM-identifiable</td><td>$\leftrightarrow$</td><td>$\phi$ partiellement injective restreinte à $\mu(P)$ tels que $\delta(\mu(P)) > \varepsilon$.</td></tr></tbody></table>

La proposition 11 a été utilisée dans [173] pour développer un algorithme partitionnant un ensemble admis- sible de paramètres en sous-ensembles SM-identifiables et μ-SM-identifiables. Cet algorithme détermine tout d'abord les sous-ensembles pour lesquels on a la μ-SM-identifiabilité grâce à l'algorithme ITVIA puis, dans le complémentaire, les ensembles SM-identifiables. Par exemple, à la figure 2.6, le sous-ensemble $A_i$ est μ-SM iden- tifiable et les sous-ensembles $A_{i, i} = 0, 2, 3, 4, 5$ sont SM-identifiables mais non disjoints. Notre algorithme trouve d'abord $A_1$ en utilisant ITVIA puis il renvoie la liste $F = \{F0 = A0, F2 = A3\}$ qui est constituée d'ensembles disjoints SM-identifiables c'est à dire que les trajectoires générées par ces ensembles connexes sont distinctes de celles générées par leurs complémentaires.

Rappelons quelques notions en analyse par intervalle. Soit $f : \mathbb{R}^n \to \mathbb{R}^m$ ($f \in C^1$) une fonction vectorielle différentiable définie sur une boîte de dimension $n$, $[x] \in \mathbb{R}^n$, où $\mathbb{R}$ est l'ensemble de tous les intervalles et $[f]$ sa fonction d'inclusion ([104] et voir l'Annexe 6.3). Un sous-pavage de $[x]$ est un ensemble de boîtes non superposées incluses dans $[x]$ (voir [34]). Il peut être considéré soit comme une collection (liste) de boîtes

<footnote n="7">7. Quand les conditions initiales ne sont pas considérées, la fonction $\varphi$ est alors $\varphi : p = (p_1, \dots, p_q) \mapsto (\gamma_1(p), \dots, \gamma_d(p))$ et la réciproque du théorème n'est plus vraie.</footnote>
$\mathbb{K} = \{[x]^{(1)}, [x]^{(2)}, ... \}$ ou comme une union $\mathbb{K} = [x]^{(1)} \cup [x]^{(2)} \cup ...$. Par conséquent, un sous-pavage peut être considéré comme un sous-ensemble discret de $\mathbb{R}^n$ ou comme un sous-ensemble convexe de $\mathbb{R}^n$. Les sous-pavages permettent d'approximer des ensembles convexes avec une précision arbitraire. Un sous-pavage régulier est un sous-pavage généré par bissections successives et il peut être facilement représenté par un arbre binaire.

L'algorithme (inspiré de [104]) peut se résumer ainsi :

1. partitionnement du domaine en sous-pavés réguliers non-minimaux de telle façon que la longueur de chaque boîte est inférieure à $\varepsilon$

2. détermination de la liste $\mathcal{L}$ d'ensembles disjoints connexes ($\mu$)-SM-identifiables (pour chaque sous-pavés réguliers, et étant donné une fonction d'inclusion, les images de boîtes de longueur $\varepsilon$ sont calculées et comparées)

3. régularisation de la liste $\mathcal{L}$ en sous-pavés réguliers.

Dans l'étape 2, la détermination des ensembles SM-identifiables est basée sur deux sous algorithmes. Le premier permet de déterminer le second optimum d'une fonction dans un intervalle, ce qui permet de délimiter éventuellement un ensemble SM-identifiable. Le deuxième vérifie quels sont les ensembles connexes d'une liste. L'étape 3 est importante et a nécessité de modifier l'algorithme proposé dans [104] section 3.4.2. Cet algorithme peut s'appliquer à une liste de boîtes se chevauchant ou non. Or, comme on peut le voir sur les figures 2.7, la liste initiale de boîtes se chevauchant (figure de gauche) est incluse dans un sous-pavage régulier.

<figure><img src="image_12.png" /><figcaption>FIGURE 2.7 – À gauche est représenté une liste $\mathcal{L}$ de boîtes se chevauchant. À droite est représenté un sous-pavage régulier externe contenant $\mathcal{L}$.</figcaption></figure>

Les images de fonctions réelles rationnelles sont surestimées, parfois de façon importante dans le cadre de l'analyse par intervalle. Ainsi, des domaines qui ne sont pas SM-identifiables pourraient être inclus dans des domaines SM-identifiables. La dernière étape a donc été ajoutée pour limiter les effets de surestimation et obtenir un sous-pavage régulier. Un exemple est donné à la figure 2.8. La figure 2.8(a) représente la liste initiale de boîtes se chevauchant, la figure 2.8(b) le sous-pavage régulier obtenu à l'étape 2 de notre algorithme, la figure 2.8(c) le sous-pavage obtenu à l'étape 3 après régularisation.

**Exemple 2.5.4.** Considérant l'équation de Bernoulli dy(t) = -p₁²y(t) + p₂²y(t)², l'algorithme appliqué à l'équation

$\phi: \mathbb{R}^2 \rightarrow \mathbb{R}^2$, $\phi(p) = \begin{bmatrix} p_1^2 \\ -p_2^2 \end{bmatrix}$ avec une précision $\epsilon = 0.05$ a donné la figure 2.9. En rouge est représenté l'ensemble $\mu$-SM-identifiable et en gris et noir les ensembles SM-identifiables.

## 2.5.4 Estimation des paramètres par des méthodes ensemblistes et propriétés

Dans cette section, je vais présenter le problème d'estimation de paramètres dans le cadre ensembliste ainsi que deux propriétés que nous avons introduites : la propriété de robustesse et la propriété de consistance.

Dans le cadre de l'estimation paramétrique classique, on dispose d'une série de données bruitées $y_m(t_i), i = 0, ..., h$ ($y_m(\cdot) \in \mathbb{R}^m$) générées par un système physique sur l'intervalle $[0, T]$. Le problème consiste alors à trouver le vecteur de paramètre $p^*$ du modèle mathématique telles que ses sorties soient le plus proches possibles,
<figure><img src="image_13.png" /><figcaption>FIGURE 2.8 – (a) : liste de boîtes se chevauchant, (b) : sous-pavage régulier externe (c) : sous-pavage régulier interne</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 2.9 – En rouge est représenté l'ensemble µ-SM-identifiable et en gris et noir les ensembles SM-identifiable.</figcaption></figure>

en un certain sens, des mesures. La méthode des moindres carrés est une méthode très utilisée et se formule de la façon suivante :

$$p^* = \arg\min_{p \in \mathcal{P}} \sum_{t=t_0}^{t_h} \|y_m(t) - y(t, p)\|^2.$$

Dans le cadre ensembliste, les données sont supposées variées dans un intervalle borné traduisant les incertitudes que l'utilisateur peut avoir sur leurs valeurs. Les bornes des intervalles peuvent être définies à partir du bruit de mesure ou de la précision des appareils et on a $y_m(t_i) \in Y_m(t_i), i = 0, \dots, h$, les $Y_m(t_i)$ étant des ensembles connexes de $\mathbb{R}^m$. L'estimation des paramètres par des méthodes ensemblistes consiste à trouver l'ensemble des vecteurs de paramètres $\mathcal{P} \subseteq \mathbb{R}^p$ telles que les trajectoires générées par $\mathcal{P}$ coïncident avec des trajectoires des données expérimentales, plus précisément :

$$p^* \in \mathcal{P} \Leftrightarrow y(t_i, p^*) \in Y_m(t_i), \forall i = 0, \dots, h.$$
P est appelé l'ensemble des paramètres réalisables. Ces problèmes sont généralement résolus avec une méthode d'arborescence de type branch and bound. L'espace des paramètres étant partitionnés, l'ensemble des trajectoires générées par chacune de ces partitions est comparé aux trajectoires mesurées. Le sous-ensemble de paramètres est alors soit réalisable, non réalisable ou indéterminé. Les ensembles non réalisables sont exclus alors que les ensembles indéterminés sont partitionnés et testés pour déterminer s'ils sont réalisables. L'algorithme réitère le processus jusqu'à ce que le diamètre des solutions candidates soit plus petit ou égal à un seuil $\epsilon$ donné par l'utilisateur.

$\epsilon$ est appelé le seuil de précision de l'algorithme. L'algorithme SIVIA (Set Inversion Via Interval Analysis) [105] rappelé dans l'annexe 6.3, section 6, utilise ce principe dans le cadre de l'analyse par intervalle. Le nombre de bissections ou de partitions du domaine est en général très important. Des algorithmes utilisent des techniques de propagation de contraintes pour réduire la dimension des sous-ensembles ou boîtes à tester. Dans ce contexte, le modèle est interprété comme l'ensemble des contraintes d'un problème de satisfaction de contraintes (CSP). Pour le résoudre, différents types de contracteurs peuvent être utilisés [34].

Il est à noter que ces algorithmes, indépendamment du temps de calcul, fournissent des solutions garanties. La solution retournée est une sur-estimation de l'ensemble des paramètres réalisables et est l'union convexe d'un ou plusieurs ensembles candidats dits réalisables et indéterminés (cf. [104]).

Regardons sur un exemple les problèmes qui peuvent être rencontrés.

**Exemple 2.5.5.** Considérons le modèle :

$$ egin{cases} 
\dot{x}_1 = (p_1 + 2(1-p_2)\cos(p_1))x_1^2 + (1-p_2)x_2, \\ 
\dot{x}_2 = \sin(p_1)x_1, \\ 
y = x_1, 
\end{cases} 
\quad (2.50) $$

où $(p_1, p_2) \in [-1,4] \times [0,1/10] = U_P$.

En posant $c_1 = \sin(p_1)$, avec l'ordre d'élimination $\{c_1, p_2\} < \{y\} < \{x_1, x_2\}$, l'algorithme de Rosenfeld-Groebner donne le polynôme différentiel :

$$ R(y, u) = \ddot{y} - 2(p_1 + 2(1-p_2)\cos(p_1))\dot{y} - (1-p_2)\sin(p_1)y. \quad (2.51) $$

Le déterminant fonctionnel est réduit à $\Delta R(y) = \det(\dot{y}, y) = -y^2\ddot{y}$ et n'est pas identiquement nul.

En étudiant l'injectivité de la fonction $\phi: (p_1, p_2) \to (p_1 + 2(1-p_2)\cos(p_1), (1-p_2)\sin(p_1))$ et grâce à l'algorithme que nous avons développé et rappelé à la section 2.5.3.3, nous obtenons la figure 2.10 (à droite). Ainsi $U_P = [-1,4] \times [0,1/10]$ est partitionné en deux domaines : le domaine rouge sur lequel la fonction $\phi$ est partiellement injective et donc ce domaine est $\mu$-SM-identifiable et deux sous-ensembles blancs pour lesquels leur image par $\phi$ a une intersection non vide.

Or tout algorithme d'optimisation est défini selon une précision $\epsilon$ (ou un seuil) que ce soit la précision machine ou la précision décidée par l'utilisateur. Pour savoir s'il permet de conserver le nombre d'ensembles solutions, on a cherché à répondre aux deux questions suivantes. Si deux ensembles de paramètres sont distants de moins de $\epsilon$, quelles sont les conditions sur le système et la précision de mesure pour décider si

— P est composé d'un seul ensemble connexe ou non ?

— les trajectoires calculées sont distinguables au sens de la précision de mesure ?

La première question nous a conduit à définir la robustesse d'un problème d'estimation de paramètres dans le cadre ensembliste et la deuxième la $\epsilon$-consistance.

Comme dans le cadre de l'estimation paramétrique classique, on souhaite savoir si lors d'un problème d'estimation de paramètres dans le cadre ensembliste, P est réduit à un seul ensemble connexe, cette propriété indiquant que le problème est bien posé mathématiquement. Aussi, nous avons introduit la définition suivante de robustesse qui est représentée à la figure 2.11.

**Définition 13.** Un problème d'estimation de paramètres dans le cadre ensembliste est dit _robuste_ si $P \subseteq U_P$ est réduit à un seul ensemble connexe. Dans ce cas, P est dit également _robuste_.

Étant donné un algorithme d'estimation de paramètres dans le cadre ensembliste avec un seuil de précision $\epsilon$, on note $P_\epsilon$ l'ensemble solution. Nous allons regarder les propriétés de $P_\epsilon$ par rapport à P (voir la figure 2.12).

**Définition 14.** Supposons que P est égal à l'union de $\kappa \ge 1$ ensembles connexes deux à deux disjoints, alors l'ensemble solution $P_\epsilon$ est dit _$\epsilon$-consistant_ si $P_\epsilon$ est égal à l'union de $\kappa_\epsilon$ ensembles connexes deux à deux disjoints et $\kappa_\epsilon = \kappa$.

La précision fixée des algorithmes peut conduire à une surestimation de P. Dans ce cas, $\kappa_\epsilon < \kappa$, c'est à dire une au moins des composantes de $P_\epsilon$ contient plusieurs ensembles composant P.

La robustesse est analysée à la section 2.5.5 et la $\epsilon$-consistance à la section 2.5.6. Un exemple dans cette dernière section permettra de tester numériquement la $\epsilon$-consistance par rapport à la précision $\epsilon$ de l'algorithme.
<figure><img src="image_15.png" /></figure>

FIGURE 2.10 – A gauche, une boule de diamètre ε comprise entre P₁ et P₂ définit la limite à partir de laquelle l’algorithme ne peut plus les distinguer. À droite, le domaine de paramètres a été partitionné en deux sous-domaines : en rouge les sous-domaines µ-SM-identifiables ; en noir, les sous-domaines SM-identifiables ; en blanc, les domaines ayant la même image. Quelle doit être la valeur de ε sur la figure de gauche et les conditions sur le système pour que les trajectoires calculées soient distinguables au sens de la précision de mesure ?

### 2.5.5 Robustesse

La condition de robustesse est liée à la condition de SM-identifiabilité. En effet, par définition, si $\mathcal{P}$ est globalement SM-identifiable, les trajectoires de (2.46) générées à partir de $\mathcal{P}$ sont distinctes des trajectoires générées par le complémentaire $\mathcal{P}^c = U_P - \mathcal{P}$. De plus, $\mathcal{P}$ étant connexe, il est robuste. Réciproquement, si $\mathcal{P}$ est robuste, par définition il est globalement SM-identifiable. D'où la proposition suivante.

**Proposition 12.** *Considérons le système (2.46) et supposons que l'ensemble $\mathcal{P} \subseteq U_P$ est l'ensemble des paramètres réalisables d'un problème d'estimation de paramètres dans le cadre ensembliste pour (2.46), alors $\mathcal{P}$ est robuste si et seulement si $\mathcal{P}$ est globalement SM-identifiable.*

Si, en plus d'être SM-identifiable, $\mathcal{P}$ est μ-SM-identifiable alors $\mathcal{P}$ préserve la robustesse quand l'erreur sur les données diminue. Dans ce cas, $\mathcal{P}$ est dit *μ-robuste* pour (2.46).

**Proposition 13.** *Étant donné l'ensemble des données $Y_m(t_i), i = 1, \dots, h$, supposons que $\mathcal{P}$ est robuste. Alors, si $\mathcal{P}$ est μ-SM-identifiable pour (2.46), l'ensemble des paramètres réalisables du même problème avec des ensembles de données contractées $\mu_i(Y_m(t_i)), i = 1, \dots, h$, où les $\mu_i$ sont des contractions, est également robuste.*

Si un problème d'estimation de paramètres est formulé tel que l'ensemble des paramètres réalisables $\mathcal{P}$ est dans $U_P$, nous pouvons désormais décider si $\mathcal{P}$ est robuste ou non. En reprenant l'exemple 2.5.5, toute partie connexe rouge est μ-robuste et toute partie connexe noire est robuste.

### 2.5.6 La ε-consistance

La ε-consistance définie à la définition 14 est une propriété de l'ensemble solution $\mathcal{P}_\varepsilon$ retournée par un algorithme d'estimation de paramètres dans le cadre ensembliste avec une précision seuil $\varepsilon$. Parmi les problèmes qui impactent la ε-consistance, nous avons regardé deux problèmes en particulier :

- l'algorithme d'estimation de paramètres dans le cadre ensembliste peut ne pas être capable de séparer des ensembles disjoints connexes de $\mathcal{P}$ en testant des ensembles de solutions candidats topologiquement pertinents,

- les trajectoires issues de paramètres solution peuvent ne pas être distinguées des trajectoires issues de paramètres non-solution, compte tenu de la précision des capteurs.

Si $\mathcal{P}$ est robuste, il est réduit à un seul ensemble connexe. D'après le principe des algorithmes branch and bound, l'ensemble solution $\mathcal{P}_\varepsilon$ est donc réduit à un seul ensemble connexe même si c'est une surestimation de $\mathcal{P}$. Ainsi, l'ensemble solution $\mathcal{P}_\varepsilon$ est ε-consistant pour tout $\varepsilon$.
<figure><img src="image_16.png" /><figcaption>FIGURE 2.11 – Représentation d'un problème d'estimation de paramètres dans le cadre ensembliste. P est l'ensemble des paramètres réalisables donc consistants avec l'ensemble des données. Il est dit robuste s'il est réduit un seul ensemble connexe.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 2.12 – Représentation d'une problématique lorsqu'on résout un problème d'estimation de paramètres dans le cadre ensembliste. Ce type de problème peut renvoyer une sur-estimation de l'ensemble P. Ici, l'ensemble P des paramètres réalisables est constitué de 2 ensembles connexes disjoints. Si l'algorithme renvoie deux ensembles solutions, l'ensemble solution P<sub>ε</sub> est dit consistant.</figcaption></figure>
Supposons maintenant que $\mathcal{P}$ n'est pas robuste et consiste en $\kappa$ ensembles mutuellement dis joints $\mathcal{P}_i, i = 1, \dots, \kappa$. Le fait que l'algorithme est capable de séparer les $\mathcal{P}_i$ est un problème topologique impliquant la distance entre les $\mathcal{P}_i$ et le diamètre des plus petits ensembles solutions candidats donnés par l'algorithme branch and bound. Une condition nécessaire pour qu'un ensemble solution $\mathcal{P}_\varepsilon$ soit $\varepsilon$-consistant est que $d(\mathcal{P}_i, \mathcal{P}_j) > \varepsilon$, pour tout $i, j = 1, \dots, \kappa, i \neq j$.

Considérons maintenant le second problème concernant les ensembles de sorties de données $Y_m(t_i), i = 0, \dots, h$, fournies par des capteurs dont la précision est de $\lambda$, i.e. $v_{mes} = v \pm \lambda$ où $v$ est la vraie valeur et $v_{mes}$ la valeur mesurée. Dans ce cas, deux trajectoires $y(\cdot, p)$ et $y(\cdot, \bar{p})$ doivent être distantes de $\lambda$, i.e. être telles qu'il existe $t \in [t_0, T], \|y(t,p) - y(t, \bar{p})\|_{\infty} > \lambda$, pour être distinguables. Si on ne peut distinguer les trajectoires provenant de paramètres qui ne sont pas solutions et celles provenant de paramètres solutions $\mathcal{P}_i$ alors $\mathcal{P}_\varepsilon$, l'ensemble solution renvoyé par l'algorithme, peut ne pas être $\varepsilon$-consistant.

La proposition suivante dont la preuve est basée sur le lemme de Gronwall montre que, sous certaines conditions, $y$ est lipschitz continue par rapport au vecteur de paramètres. Elle donne explicitement la constante de Lipschitz $K_{y,p}$. Cette proposition permet également de déterminer la distance maximale entre deux vecteurs de paramètres $p$ et $\bar{p}$ pour que la distance entre les trajectoires $y(t,p)$ et $y(t, \bar{p})$ soit inférieure à $\lambda$.

Rappelons que, si une fonction $g$ est réelle et analytique sur $M$, un ensemble ouvert de $\mathbb{R}^n$, pour tout ensemble compact $K \subset M$, il existe une constante $K > 0$ telle que pour tout $x$ dans $K$, on a la majoration suivante :

$$ \| \frac{dg}{dx}(x) \|_{\infty} \leq K, $$

$f$ et $h$ définissant le modèle (2.46) sont supposées être réelles et analytiques sur $M$. Si $x \in K$, $K$ un compact, alors $f$ et $h$ sont Lipschitz continues selon $x$.

Considérons les hypothèses suivantes :

i) $f$ et $h$ définies sur $[t_0, T] \times \mathcal{U}_\mathcal{P}$ sont Lipschitz continues selon $p$ et telles que leurs constantes de Lipschitz sont notées $K_{f,p}$ et $K_{h,p}$ respectivement;

ii) la solution $x(t,p)$ de (2.46) est dans le compact $K$;

iii) si les conditions initiales dépendent de $p$, la fonction $p \mapsto x(t_0,p)$ est supposée être Lipschitz continue selon $p$ et sa constante de Lipschitz est notée $K_{x_0,p}$.

**Proposition 14.** Supposons que les conditions i), ii), et iii) sont vérifiées. Alors $y$ est Lipschitz continue selon $p$ et sa constante de Lipschitz $K_{y,p}$ est donnée par :

$$ K_{y,p} = K_{h,x}(K_{x_0,p} + K_{f,p}(T - t_0))e^{K_{f,x}(T - t_0)} + K_{h,p}. $$

Si les conditions initiales ne dépendent pas de $p$, $K_{y,p}$ est donnée par :

$$ K_{y,p} = K_{h,x} K_{f,p} (T - t_0) e^{K_{f,x}(T - t_0)} + K_{h,p}. $$

De plus, les deux trajectoires $y(\cdot, p)$ et $y(\cdot, \bar{p})$ vérifient :

$$ \left( \forall p, \bar{p} \in \mathcal{U}_{\mathcal{P}}, \|p - \bar{p}\| \le \frac{\lambda}{K_{y,p}} \right) \Rightarrow (\forall t \in [t_0, T], \|y(t,p) - y(t, \bar{p})\|_{\infty} \le \lambda). $$

Ce résultat signifie également que $\mathcal{P}_i$ formant une partie de $\mathcal{P}$ générera des trajectoires qui ne seront pas distinguables de celles générées par son voisinage.

La proposition suivante donne une condition pour avoir la non $\varepsilon$-consistance.

**Proposition 15.** Soit le système (2.46) et les hypothèses i), ii), iii). Supposons que $\mathcal{P}$ n'est pas robuste et est constitué de $\kappa$ ensembles connexes mutuellement dis joints $\mathcal{P}_i, i = 1, \dots, \kappa$. Soit $\lambda$ la distance minimale permettant de distinguer les trajectoires et $\varepsilon$ la précision de l'algorithme d'estimation de paramètres dans le cadre ensembliste.

Si $d(\mathcal{P}_i, \mathcal{P}_j) \le \varepsilon + \frac{2\lambda}{K_{y,p}}$, pour tout $i, j = 1, \dots, \kappa$ alors l'ensemble solution $\mathcal{P}_\varepsilon$ retourné par l'algorithme d'estimation de paramètres dans le cadre ensembliste n'est pas $\varepsilon$-consistant.

**Remarque 13.** La réciproque est vraie si la fonction d'inclusion $P \to [Y](P)$ utilisée dans l'algorithme d'estimation (d'estimation de paramètres dans le cadre ensembliste) pour prédire l'ensemble des trajectoires générées par un ensemble de paramètres candidat $P$ est telle que $[Y](P) = Y(P)$, ce qui est rarement le cas.

**Exemple 2.5.6.** Considérons le système suivant défini sur $[0, T]$ :

$$ \left\{ \begin{array}{l} \dot{x}_1(t,p) = \cos(p)x_2(t,p), x_1(0) = (\frac{1}{\sqrt{2}} + 1)/2, \\ \dot{x}_2(t,p) = -x_1(t,p)\cos(p) + (1 - x_1(t,p)^2 - 2x_2(t,p)^2)x_2(t,p), x_2(0) = 0, \\ y(t,p) = x_1(t,p). \end{array} \right. \qquad (2.52) $$
Les fonctions f et h sont définies par $f(x,p) = (\cos(p)x_2, -x_1 \cos(p) + (1-x_1^2-2x_2^2)x_2)^T$, et $h(x,p) = x_1$ où
$x = (x_1,x_2)^T$ et $T$ représente la transposée du vecteur considéré. La solution $(x_1(t),x_2(t))^T$ reste dans l'anneau
$\mathbb{R}$ défini par les deux cercles centrés en $(0,0)$ de rayons $\frac{1}{\sqrt{2}}$ et $1$ respectivement. Pour déterminer les constantes de
Lipschitz, considérons $z = (z_1,z_2)^T \in \mathbb{R}$. Clairement, $K_{h,p}=0$, $K_{h,x}=1$ et $K_{f,p}=1$ comme $\|x\|=\|(x_1,x_2)^T\|<$
$1$. Pour $K_{f,x}$, en ordonnant les termes et en ajoutant $x_1^2z_2 - x_1^2z_2$ à la ligne 3, nous obtenons :

$$
\begin{align}
\| f(x,p) - f(z,p) \|_{1} &\le |\cos(p)|(|x_2 - z_2| + |x_1 - z_1|) + |x_2 - z_2| + |(x_1^2 - 2x_2^2)x_2 - (z_1^2 - 2z_2^2)z_2| \nonumber \\
&\le (|x_1 - z_1| + 2|x_2 - z_2|) + |x_1^2 x_2 - x_1^2 z_2 + x_1^2 z_2 - z_1^2 z_2 - 2(x_2^3 - z_2^3)| \tag{2.53} \\
&\le (|x_1 - z_1| + 2|x_2 - z_2|) + x_1^2 |x_2 - z_2| + |z_2| |x_1^2 - z_1^2| + 2|x_2^3 - z_2^3|. \nonumber
\end{align}
$$

Comme $|x_1| < 1$, $|z_1| < 1$, nous en déduisons d'une part que $|x_1^2 - z_1^2| \le |x_1 - z_1|(|x_1| + |z_1|) \le 2|x_1 - z_1|$ et d'autre part que $|z_2^3 - x_2^3| = |x_2 - z_2||z_2^2 + 2x_2z_2 + x_2^2| \le 4|x_2 - z_2|$. Ainsi

$$
\| f(x,p) - f(z,p) \|_{1} \le (|x_1 - z_1| + 2|x_2 - z_2|) + |x_2 - z_2| + 2|x_1 - z_1| + 8|x_2 - z_2| \le 11 \|x - z\|_{1}. \quad (2.54)
$$

En utilisant l'équivalence de la norme 1 et de la norme infinie, il vient :

$$
\|f(x,p) - f(z,p)\|_{\infty} \le 22 \|x - z\|,
$$

et donc $K_{f,x} = 22$. D'après la proposition 14, la constante de Lipschitz $K_{y,p}$ est égale à $T e^{22T}$. Supposons
que $T = 1$, que la précision du capteur est $\lambda = 0,01$ et que le seuil de précision de l'algorithme d'estimation
de paramètres dans le cadre ensembliste est égal à $\varepsilon = 10^{-16}$. Alors d'après la Proposition 15, si pour tout
$i,j \in \{1,\dots,\kappa\}$, $i \neq j$, $d(\mathcal{P}_i, \mathcal{P}_j) \le \varepsilon + \frac{2\lambda}{K_{y,p}} \simeq 10^{-16} + 5,579.10^{-12} \simeq 5,579.10^{-12}$, l'ensemble solution n'est pas
$\varepsilon$-consistant. Ainsi, si cette distance est inférieure à $\varepsilon + \frac{2\lambda}{K_{y,p}}$, le nombre d'ensembles renvoyés par l'algorithme
sera inférieur au nombre attendu d'ensembles solutions.

2.5.7 Estimation des paramètres dans le cadre ensembliste

2.5.7.1 Une enveloppe de la solution initiale

En reprenant les notations de la partie 2.4, nous avons étendu les méthodes d'estimation de paramètres
au cadre de l'analyse par intervalle. Nous avons supposé que les bruits de mesure et les incertitudes sur les
paramètres p ∈ P, P étant un vecteur d'intervalle, étaient pris en compte dans le bruit de la sortie y. Ce bruit
est supposé additif et borné η(t) ∈ [η(t)]. Γ<sub>j</sub>(P) obtenu à partir de γ<sub>j</sub>(p) en substituant p par l'ensemble connexe
P est alors un ensemble connexe. Les problèmes d'estimation proposés dans la partie 2.4 se réécrivent de la
façon suivante.

1. Méthode basée sur les polynômes ES

Le système (2.11) peut s'interpréter dans le cadre ensembliste de la façon suivante

$$
[A][\gamma] = [b] \tag{2.55}
$$

où $[A]_k = ([m_j(y_k, u_k)])_{j=1,...,q}$ and $[b_k] = -[m_0(y_k, u_k)]$. L'ensemble des solutions du système (2.55)
défini par le vecteur intervalle $(\Gamma_j(P))_{1\le j\le q}$ peut être trouvé en résolvant le système

$$
0 \in [A][\gamma] - [b]. \tag{2.56}
$$

2. Méthode basée sur les polynômes intégro-differentiels

Le système (2.15) interprété dans le cadre ensembliste s'écrit

$$
[A^{\psi}][\gamma] = [b^{\psi}] \tag{2.57}
$$

où $[A^\psi]_i = \left( \left[ \int_{t_i-\epsilon}^{t_i+\epsilon} m_j(y,u)(s)\psi_i(s)ds \right]_{j=1,...,q} \right)$ et $[b]^\psi_i = -\left[ \int_{t_i-\epsilon}^{t_i+\epsilon} m_0(y,u)(s)\psi_i(s)ds \right]$. Les solutions de
ce système peuvent être obtenues en résolvant le système suivant ($\gamma = (\gamma_j)_{k=1,...,q}$)

$$
0 \in [A^{\psi}][\gamma] - [b^{\psi}]. \tag{2.58}
$$

La question qui vient naturellement est de savoir quelles sont les conditions pour que l'ensemble solution du système $[A^\psi][\gamma] = [b^\psi]$ contient l'ensemble solution du système $[A][\gamma] = [b]$ ? Le système (2.15) peut-être vu comme un système perturbé du système (2.11). Ainsi, une condition pour que les ensembles solutions du système (2.58) enveloppent les ensembles solutions de (2.56) est que $0 \in [\delta A][\gamma^\psi] - [\delta b]$ (voir [197]).
### 2.5.7.2 Application

Les fonctions tests sont celles de la partie 2.4.2.4 et l'enveloppe de $\int_a^b f(s)ds$ est obtenue en étendant la méthode classique des trapèzes au cadre ensembliste. Elle est notée $I_{[a,b]}(f)$. Enfin, l'enveloppe des dérivées est calculée en utilisant les différenciateurs HOSM [124, 172] avec les paramètres $\lambda_0 = 3$, $\lambda_1 = 0.2$ et $\lambda_2 = 0.1$.

L'exemple choisi est celui de deux réservoirs couplés l'un au dessus de l'autre par

$$ 
\left\{ 
\begin{aligned}
\dot{x}_1(t,p) &= a_1 u(t) - a_2 \sqrt{x_1(t,p)}, \\
\dot{x}_2(t,p) &= a_3 \sqrt{x_1(t,p)} - a_4 \sqrt{x_2(t,p)}, \\
y_1(t,p) &= \sqrt{x_1(t,p)}, \\
y_2(t,p) &= \sqrt{x_2(t,p)}
\end{aligned}
\right. 
\qquad (2.59) 
$$

où $p = (a_i)_{i=1,...,4}, a_i \neq 0$, est le vecteur de paramètre, $x = (x_1, x_2)^T$ représente le vecteur d'état et correspond au niveau d'eau dans chaque réservoir. $u \neq 0$ est le vecteur d'entrée. Le niveau d'eau dans les réservoirs peut varier entre 0 et 10.

Pour utiliser l'algorithme de Rosenfeld-Groebner implementé dans Maple 16, nous avons ajouté les variables auxiliaires $z_1(t,p) = \sqrt{x_1(t,p)}$ et $z_2(t,p) = \sqrt{x_2(t,p)}$ et le modèle contenant quatre fautes est ré-écrit

$$ 
\left\{ 
\begin{aligned}
\dot{x}_1(t,p) &= a_1 u(t) - a_2 z_1(t,p), \\
\dot{x}_2(t,p) &= a_3 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) &= z_1(t,p), y_2(t,p) = z_2(t,p),
\end{aligned}
\right. 
\qquad (2.60) 
$$

D'après l'algorithme de Rosenfeld-Groebner, les polynômes entrée-sortie sont :

$$ 
\begin{aligned}
R_1(y, u, p) &= -a_1 u + a_2 y_1 + 2 y_1 \dot{y}_1, \\
R_2(y, u, p) &= a_4 y_2 - a_3 y_1 + 2 y_2 \dot{y}_2
\end{aligned} 
\qquad (2.61) 
$$

Le déterminant fonctionnel $\Delta R_1 = u\dot{y}_1 - \dot{u}y_1$ n'est pas identiquement nul tant que $u \neq 0$. Donc $\Delta R_2 = y_2\dot{y}_1 - \dot{y}_2y_1$ n'est pas identiquement nul si $u \neq 0$ et la fonction $\phi : (a_1, a_2, a_3, a_4) \mapsto (a_1, a_2, a_4, a_3)$ étant clairement injective, tout sous-ensemble connexe $P^*$ inclus dans l'ensemble des paramètres admissibles est $\mu$-SM-identifiable.

1. Méthode basée sur les polynômes ES

En notant $y_{1,p}^{(1)}(t_k)$ (resp. $y_{2,p}^{(1)}(t_k)$) l'estimée de $\dot{y}_1(t_k)$ (resp. $\dot{y}_2(t_k)$), le premier système à résoudre est $[A_1][\gamma] = [b_1]$ où

$$ 
\begin{aligned}
[A_1]_k &= ([-u(t_k)], [y_1(t_k)]), \\
[b_1]_k &= [-2y_1(t_k)y_{1,p}^{(1)}(t_k)].
\end{aligned} 
\qquad (2.62) 
$$

Le second système est composé de la matrice intervalle $[A_2]$ et du vecteur intervalle $[b_2]$ tel que $[A_2]_i = ([-y_1(t_i)], [y_2(t_i)])$ et $[b_2]_i = [-y_2(t_i), y_{2,p}(t_i)]$.

2. Méthode basée sur les polynômes intégro-differentiels

Le premier système à résoudre est $[A_1^\psi][\gamma] = [b_1^\psi]$ où

$$ 
\begin{aligned}
[A_1^\psi]_i &= ([-I_{[t_i-\epsilon, t_i+\epsilon]}(u\psi_i)], [I_{[t_i-\epsilon, t_i+\epsilon]}(y_1\psi_i)]), \\
[b_1^\psi]_i &= [I_{[t_i-\epsilon, t_i+\epsilon]}(y_1^2\dot{\psi}_i)].
\end{aligned} 
\qquad (2.63) 
$$

Le deuxième est composé de $[A^\psi]$ et $[b^\psi]$ tel que

$$ 
\begin{aligned}
[A_2^\psi]_i &= ([-I_{[t_i-\epsilon, t_i+\epsilon]}(y_1 \psi_i)], [I_{[t_i-\epsilon, t_i+\epsilon]}(y_2 \psi_i)]), \\
[b_2^\psi]_i &= [I_{[t_i-\epsilon, t_i+\epsilon]}(y_2^2 \dot{\psi}_i)].
\end{aligned} 
\qquad (2.64) 
$$

Pour les tests de simulation, les paramètres sont $a_1 = 0.2$, $a_2 = 0.05$, $a_3 = 0.1$ et $a_4 = 0.1$. les bornes de l'intervalle pour $\eta(t)$ sont données par $[-0.1 \ 0.1]$. Les intervalles initiaux pour les paramètres sont $\mathbf{p}_0 = [0; 5]$ et
<figure><img src="image_18.png" /><figcaption>FIGURE 2.13 – À gauche est représenté le niveau d'eau dans le réservoir supérieur; à droite, le niveau d'eau dans le réservoir supérieur.</figcaption></figure>

le seuil de bissection pour l'algorithme SIVIA est de 0.025 pour chaque paramètre. Les figures 2.13 représentent les niveaux d'eau des réservoirs supérieurs et inférieurs.

Les méthodes d'estimation ont été appliquées à l'intervalle de temps [0, 20] avec une discrétisation $(t_i)_{i=1,...,M}$ de pas égal à 0.1.

Les résultats obtenus à partir de la première méthode sont donnés aux figures 2.14 et 2.15. Le temps de calcul est de 24.00 secondes. Ceux de la deuxième méthode sont représentés aux figures 2.16 et 2.17 et le temps de calcul est de 16.41 secondes.

Les boîtes rouges sont les boîtes qui ont été rejetées, en jaune, les boîtes indéterminées, c'est à dire les boîtes contenant des solutions mais dont la longueur est plus petite que la condition d'arrêt. Enfin, les boîtes vertes ne contiennent que des solutions.

Le boîtes vertes contiennent les valeurs exactes des paramètres en utilisant la première ou la seconde méthode. Cependant, la deuxième méthode fournit des solutions beaucoup moins conservatives, i.e. les régions vertes et jaunes sont plus petites.

Les tableaux 2.3 et 2.4 donnent le pourcentage de boîtes initiales éliminées par rapport aux domaines $\bar{S}$ (boîtes de solution et indéterminées) et $\bar{S}$ (boîtes de solution uniquement), respectivement définis dans l'équation (6.10).

Le pourcentage éliminé $\%\mathbf{p}$ est calculé par $\%\mathbf{p} = 1 - \frac{w([\mathbf{p}])}{w([\mathbf{p}_0])}$ où $w([\alpha])$ est la largeur de l'intervalle $[\alpha]$.

<table><caption>TABLE 2.3 – Pourcentage de boîtes éliminées pour les deux méthodes (boîtes indéterminées et solutions).</caption><thead><tr><th>Paramètres</th><th>%<b>p</b><sub>method1</sub></th><th>%<b>p</b><sub>method2</sub></th></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>86.00</td><td>95.00</td></tr><tr><td>a<sub>2</sub></td><td>90.00</td><td>97.00</td></tr><tr><td>a<sub>3</sub></td><td>82.00</td><td>89.00</td></tr><tr><td>a<sub>4</sub></td><td>70.00</td><td>91.00</td></tr></tbody></table>

## 2.6 Identifiabilité dans le cadre des équations aux dérivées partielles (EDP)

Suite à mes travaux sur l'identifiabilité dans les systèmes d'EDOs, il m'a paru intéressant de voir si il était possible de généraliser la méthode entrée-sortie aux systèmes d'EDPs et quels en seraient les limites. En effet, étant donné un système d'EDOs et/ou d'EDPs, il est toujours possible, du moins formellement, d'utiliser l'al-
<figure><img src="image_19.png" /><figcaption>FIGURE 2.14 – Représentation des solutions du système (2.62) par rapport aux paramètres admissibles $a_1, a_2$ ($a_1$ sur l'axe x, $a_2$ sur l'axe y). Les boîtes rouges sont les boîtes qui ont été rejetées, en jaune, les boîtes indéterminées, c'est à dire les boîtes contenant des solutions mais dont la longueur est plus petite que la condition d'arrêt. Enfin, les boîtes vertes ne contiennent que des solutions.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 2.15 – Représentation des solutions par rapport aux paramètres admissibles $a_3, a_4$ pour la méthode 1 ($a_3$ sur l'axe x, $a_4$ sur l'axe y).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 2.16 – Représentation des solutions du système (2.63) par rapport aux paramètres admissibles $a_1, a_2$ ($a_1$ sur l'axe x, $a_2$ sur l'axe y). Les boîtes rouges sont les boîtes qui ont été rejetées, en jaune, les boîtes indéterminées, c'est à dire les boîtes contenant des solutions mais dont la longueur est plus petite que la condition d'arrêt. Enfin, les boîtes vertes ne contiennent que des solutions.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIGURE 2.17 – Représentation des solutions par rapport aux paramètres admissibles $a_3, a_4$ pour la méthode 1 ($a_3$ sur l'axe x, $a_4$ sur l'axe y).</figcaption></figure>

La difficulté dans cette méthode consiste alors à obtenir un polynôme ES bien défini au sens de la régularité de la solution du modèle initial. En effet, il apparaît souvent que le polynôme contient des dérivées en la sortie d'ordre plus élevé que la régularité de la solution. Pour utiliser tout de même ces relations polynômiales aussi bien d'un point de vue théorique que numérique, nous avons été amené à construire de nouvelles relations adaptées au problème initial. L'exemple qui nous a guidé pour écrire les articles [200, 229, 228] et définir un cadre général pour les EDP est celui de la transmission du chikungunya à l'homme. En effet, le simple modèle EDO ne prend pas en compte un phéno-
<table><caption>TABLE 2.4 – Pourcentage de boîtes éliminées pour les deux méthodes (boîtes solutions).</caption><thead><tr><th>Paramètres</th><th>%<b>p</b><sub>method1</sub></th><th>%<b>p</b><sub>method2</sub></th></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>93</td><td>95.80</td></tr><tr><td>a<sub>2</sub></td><td>93</td><td>97.60</td></tr><tr><td>a<sub>3</sub></td><td>92.20</td><td>93.00</td></tr><tr><td>a<sub>4</sub></td><td>88.00</td><td>93.60</td></tr></tbody></table>

mène important dans les maladies épidémiologiques qui est la ré-émergence de la maladie du fait de la mobilité humaine. Pour spatialiser ce modèle, nous avons choisi d'ajouter un terme de diffusion comme dans [149] ou [127] et de considérer sa solution dans $L^2$.

Dans la littérature, l'identifiabilité dans les systèmes d'équations aux dérivées partielles porte essentiellement sur les problèmes sources ou le problème de Calderòn (voir [5, 56, 71, 88, 150, 191, 224]). Ceux-ci consistent à identifier un terme source ou à déterminer un paramètre variable en espace à partir de mesures au bord. Dans notre cas, nous avons supposé que la solution est connue sur tout le domaine ou sur une partie de celui-ci. Même si cette hypothèse semble moins considérée dans la littérature, on trouve quelques travaux. Par exemple, les auteurs de [48, 221] ont proposé une étude d'identifiabilité sur un modèle de réaction-diffusion en épidémiologique de type SIS (Susceptible - Infecté - Susceptible) en utilisant une méthode de contrôle optimale. Dans [157], les auteurs utilisent des relations entrée-sortie pour étudier l'identifiabilité d'un système d'équations intégrDIFFÉRENTIELLE non linéaires de type équation de transport. Une approche basée sur la théorie des semi-groupes pour identifier le taux de mortalité dans un modèle de population structuré en âge a été proposée par [158]. Dans le cas de systèmes couplés d'ODEs et d'EDPs, une étude d'identifiabilité a été faite dans [122] en utilisant les estimateurs de type Carleman. Dans [189], les auteurs ont utilisé une approche basée sur l'algèbre différentielle lorsque le modèle EDP est un modèle de population structuré en âge par rapport au début de l'infection.

## 2.6.1 Notations

Cette première partie introduit les notations utiles pour la suite de la section où on considérera

— $\Omega$ un domaine borné de $\mathbb{R}^n$ et $C^2$ au bord;

— $C_l([0, \infty); C(\bar{\Omega}))$ l'espace des fonctions continues par morceaux de $[0, \infty)$ sur $C(\bar{\Omega})$ et $C_l^1([0, \infty); C(\bar{\Omega}))$ l'espace des fonctions différentiables par morceaux de $[0, \infty)$ sur $C(\bar{\Omega})$;

— $T_F$ la distribution régulière de $F$ qui est une fonction localement intégrable;

— $H_N^2(\Omega) := \{v \in H^2(\Omega) \mid \frac{\partial v}{\partial \nu} = 0 \text{ on } \partial\Omega\}$, où $\nu$ est le vecteur normal unitaire sortant de $\partial\Omega$;

— $\langle f \rangle := \frac{1}{|\Omega|} \int_{\Omega} f(x)dx$ où $f \in L^1(\Omega)$.

## 2.6.2 Le modèle considéré

Le type de modèle que nous avons considéré est un modèle maître-esclave de deux espèces. Sa particularité est que l'espèce esclave se déplace contrairement à l'espèce maître. Ce type de modèle peut s'appliquer par exemple au modèle de transmission du Chikungunya. En effet, dans la majorité des cas, le vecteur de transmission, c'est à dire le moustique tigre ne se déplace que très peu par rapport à l'homme.

Le modèle maître est alors un système d'EDOs s'écrivant :

$$ \left\{ \begin{array}{l} \partial_t \zeta = F(t, \zeta), \quad \text{in } (0, \infty), \\ \zeta(0) = \zeta_0. \end{array} \right. \qquad (2.65) $$

où $\zeta_0 \in C(\bar{\Omega})$ et $F: \mathbb{R} \times \mathbb{R}^p \to \mathbb{R}^p$ est continue par morceaux en la première composante et lipschitzienne par rapport à la seconde composante. D'après le théorème de Cauchy-Lipschitz, il existe une unique solution différentiable par morceaux $\zeta$ telle que $\zeta \in C_l^1([0, \infty); C(\bar{\Omega}))$.

Le modèle esclave est quant à lui un système d'équations de réaction-diffusion défini sur l'espace de Banach $X := C(\Omega; \mathbb{R}^{q_0}) \times L^2(\Omega; \mathbb{R}^{q_1})$ et qui dépend de la solution $\zeta$ du modèle maître :

$$ \left\{ \begin{array}{l} \partial_t U + BU = G(t, \zeta, U), \quad \text{in } (0, \infty), \\ U(0) = U_0, \end{array} \right. \qquad (2.66) $$
où $U_0 \in X$ et l'opérateur diagonal $B := \text{diag}(B_i)$ sur $X$ satisfait :

- pour tout $i \in \{1, \dots, q_0\}$, $B_i := I$,

- pour tout $i \in \{q_0+1, \dots, q\}$ (où $q := q_0+q_1$), $B_i := d_i A + b_i$ avec $d_i, b_i > 0$ et $A$ est la réalisation du Laplacien dans $L^2(\Omega)$ avec des conditions homogènes de Neumann sur le bord $\partial\Omega$.

Ainsi, (2.66) est composé de $q_0$ EDOs, $q_1$ EDPs et qui sont couplées.

De plus, $G(t, \zeta, \cdot)$ est un opérateur non linéaire de $D(B^\eta)$ dans $X$ où $0 < \eta < 1$, qui satisfait la condition de type lipschitz suivante : il existe une fonction croissante $\phi$ telle que, pour tout intervalle compact $I$ sur lequel $\zeta$ est différentiable et pour tout couple $(t, U), (s, V) \in I \times D(B^\eta)$, $G$ vérifie

$$ \|G(t, \zeta, U) - G(s, \zeta, V)\|_{X} \le \phi(\|U\|_{X} + \|V\|_{X}) \\ \times [\|B^{\eta}(U-V)\|_{X} + (\|B^{\eta}U\|_{X} + \|B^{\eta}V\|_{X} + 1)(|t-s| + \|U-V\|_{X})], \quad (2.67) $$

où $\|\cdot\|_X$ est la norme de $X$.

Une simple adaptation du Théorème 4.4 de [222] (p. 188 et p. 199 pour le cas non autonome) implique qu'il existe $T_{U_0} > 0$ tel que (2.66) admet une unique solution locale

$$ U \in C_l([0, T_{U_0}]; D(B)) \cap C_l^1((0, T_{U_0}]; X). $$

Notons que $U$ est continue et différentiable par morceaux due à sa dépendance à $\zeta$.

De plus, s'il existe $C_{U_0} > 0$ telle que la solution locale $U$ vérifie

$$ \|U(t)\|_2 \le C_{U_0}, \quad (2.68) $$

sur tout intervalle compact $I$ où $U$ est continue par rapport au temps, alors $U$ est une solution globale c'est à dire :

$$ U \in C_l([0, \infty); D(B)) \cap C_l^1((0, \infty); X). $$

Soit $\theta := (\theta_a, \theta_b) \in L^\infty(\Omega; \mathcal{U}_p)$, où $\mathcal{U}_p$ est un sous-ensemble ouvert de $\mathbb{R}^p$, et le vecteur sortie. Les modèles (2.65) et (2.66) peuvent se ré-écrire :

$$ (\Gamma_a) \left\{ \begin{array}{ll} \partial_t \zeta(t, \theta_a) & = f(t, \zeta, \theta_a), \quad t \in [0, \infty), \\ y_a(t, \theta_a) & = h_a(\zeta), \quad t \in [0, \infty), \end{array} \right. \quad (2.69) $$

complété des conditions initiales $\zeta(t=0, \theta_a) = \zeta_0$, et

$$ (\Gamma_b) \left\{ \begin{array}{ll} \partial_t U(t, \theta_b) + B(\theta_b) U(t, \theta_b) & = g(\zeta, U, \theta_b), \quad t \in [0, \infty), \\ y_b(t, \theta_b) & = h_b(\zeta, U), \quad t \in [0, \infty), \end{array} \right. \quad (2.70) $$

complété des conditions $U(t=0, \theta_b) = U_0$ et telles que :

- $f := (f_1, \dots, f_p)$ où $f_i$, $1 \le i \le p$, est une fonction rationnelle en $\zeta$ et $\theta_a$,

- $g := (g_1, \dots, g_q)$ où $g_i$, $1 \le i \le q$, est une fonction rationnelle en $\zeta$, $U$ et $\theta_b$,

- $h_a := (h_{a,1}, \dots, h_{a,m_a})$ (resp. $h_b := (h_{b,1}, \dots, h_{b,m_b})$) où $h_{a,i}$, $1 \le i \le m_a$ (resp. $h_{b,i}$, $1 \le i \le m_b$), sont des fonctions rationnelles en $\zeta$ (resp. $\zeta$ et $U$).

Par la suite, $(\Gamma_\theta)$ désignera $(\Gamma_a)$ ou $(\Gamma_b)$.

### 2.6.3 Cas général

La définition classique de l'identifiabilité (1) a été adaptée à notre problème. En effet, les paramètres peuvent être des données spatiales. La définition est la suivante :

**Définition 15.** Le modèle $(\Gamma_\theta)$ est identifiable en $\theta$ si il existe un temps $t_1 > 0$ et un sous-ensemble $\bar{\Omega}$ de $\Omega$ tel que, pour tout $\theta, \tilde{\theta} \in L^\infty(\Omega; \mathcal{U}_p)$, nous avons

$$ \forall (t, \mathbf{x}) \in [0, t_1] \times \bar{\Omega}, \quad y(t, \theta)(\mathbf{x}) = y(t, \tilde{\theta})(\mathbf{x}) \implies \theta = \tilde{\theta}. $$

Ainsi, pour qu'un modèle soit identifiable, il doit exister un intervalle de temps et un domaine tel que deux jeux de paramètres distincts ne généreront pas les mêmes sorties. Cette définition peut-être facilement étendue à des sorties continues par morceaux en considérant l'identifiabilité du modèle sur chaque sous intervalle sur lequel les sorties sont définies.
Là encore, l'étude va se ramener à étudier les coefficients de polynômes en les sorties et les entrées du modèle. Ces derniers peuvent être obtenus grâce à un logiciel de calcul formel. Toutefois, dans ce cas, la procédure d'élimination est faite de façon formelle et ne nous garantit pas que la fonction polynômiale associée au polynôme entrée-sortie soit bien définie par rapport à la régularité de la solution initiale. Pour remédier à ce problème et gagner de la régularité sur la fonction polynômiale, nous avons proposé des méthodes permettant de montrer l'identifiabilité du modèle, résultat que nous n'avions pas en utilisant directement le polynôme de départ.

Rappelons que le polynôme entrée-sortie peut se ré-écrire :

$$ P(y_{\theta}, \theta, u) = m_0(y_{\theta}, u) + \sum_{k=1}^{n} \gamma_k(\theta) m_k(y_{\theta}, u) \quad \text{et} \quad P(y_{\theta}, \theta, u) = 0, \qquad (2.71) $$

où $(\gamma_k(\theta))_{1\le k\le n}$ ($\theta_\alpha \ne \theta_\beta$ si $\alpha \ne \beta$ et $\alpha, \beta \in \{1, \dots, n\}$) sont rationnels en $\theta$, $(m_k(y_\theta, u))_{0\le k\le n}$ sont des polynômes différentiables en $y_\theta$ et $u$, et $m_0 \ne 0$.

**Remarque 14.** Pour le système $(\Gamma_b)$, la fonction $\zeta$ peut-être considérée comme une entrée.

La proposition 3 ne peut pas être reprise directement pour montrer l'identifiabilité du modèle $\Gamma_\theta$. En effet, les ensembles semi-algébriques sont définis sur $K = \mathbb{R}$ ou $\mathbb{C}$ or ici nous travaillons avec des fonctions $L^\infty$. Toutefois, elle peut être adaptée. Supposons que la fonction polynômiale associée au polynôme $P$ est bien définie au sens de la régularité de la solution, alors on a la proposition suivante :

**Proposition 16.** Supposons que la famille de fonctions $m_k(y_\theta, u)$, $1 \le k \le d$ définie en (2.71) est linéairement indépendante. Si pour tout $(\theta, \tilde{\theta}) \in L^\infty(\Omega, \mathcal{U}_p)^2$, nous avons

$$ \forall k \in \{1, \dots, d\}, \quad \gamma_k(\theta) = \gamma_k(\tilde{\theta}) \implies \theta = \tilde{\theta}, $$

alors le modèle est identifiable en $\theta$.

Dans le cas où la dérivée de plus haut degré de $y$ par rapport au temps est trop grande par rapport la régularité de la solution, nous proposons de ré-écrire le polynôme au sens des distributions. Supposons que dans (2.71), les $m_k(y_\theta, u)$, $k = 0, \dots, d$, sont définis presque partout et sont Lebesgue intégrables. On peut alors remplacer $m_k(y_\theta, u)$ par sa distribution $T_{m_k(y_\theta, u)}$. Ainsi, la nouvelle relation obtenue est :

$$ T_{m_0(y_\theta, u)} + \sum_{k=1}^{n} \gamma_k(\theta) T_{m_k(y_\theta, u)} = 0. \qquad (2.72) $$

Similairement à la proposition 16, nous avons le résultat suivant au sens des distributions.

**Proposition 17.** Supposons que le famille de distributions $(T_{m_k(y_\theta, u)})_{1\le k\le n}$ provenant de (2.72) est linéairement indépendante. Si, pour tout $(\theta, \tilde{\theta}) \in L^\infty(\Omega, \mathcal{U}_p)^2$, nous avons

$$ \forall k \in \{1, \dots, n\}, \quad \gamma_k(\theta) = \gamma_k(\tilde{\theta}) \implies \theta = \tilde{\theta}, $$

alors le modèle considéré au sens des distributions est identifiable en $\theta$.

**Remarque 15.** Il est facile de montrer que si les $m_k(y_\theta, u)$, $k = 1, \dots, n$ sont linéairement indépendants alors les $T_{m_k(y_\theta, u)}$, $1 \le k \le n$ sont linéairement indépendants également. Cette remarque peut être importante dans les applications pour prouver l'hypothèse d'indépendance.

### 2.6.4 Résultats d'identifiabilité

Résultat d'identifiabilité pour $(\Gamma_a)$ quand $\theta_a \in L^\infty$

La sortie $y_a$ est en général différentiable par morceaux en temps (ou juste continue). Utilisant cette propriété, on peut obtenir une relation polynômiale plus informative et prenant en compte les points de discontinuité de $y_a$ (rappelons que seule la dérivation en temps est considérée). Cette nouvelle relation peut alors permettre de conclure sur l'identifiabilité du modèle. Regardons l'exemple suivant pour nous convaincre.

**Exemple 2.6.1.** Considérons le simple système d'EDO suivant :

$$ \begin{cases} \dot{x}_1 = u, \\ \dot{x}_2 = x_1 + k u, \end{cases} $$

où $k \in \mathbb{R}$ est le paramètre à identifier ($\theta = k$) et $u$ est une fonction entrée.
Supposons que $u$ est $C^1$ et que la sortie soit $y := x_2$. Alors, $y$ est $C^2$ et nous avons

$$ \ddot{y} = \dot{x}_1 + k \dot{u} = u + k \dot{u}. \qquad (2.73) $$

Nous en déduisons que le polynôme ES est :

$$ P(y, k, u) := \ddot{y} - u - k\dot{u}. $$

Si $u$ est constant, $k$ n'est pas identifiable.

Maintenant, supposons que $u$ est une fonction constante par morceaux et est définie par

$$ u(t) := \begin{cases} u_1 & \text{if } t \in (0, t_1) \\ u_2 & \text{if } t \in (t_1, T), \end{cases} $$

où $u_1, u_2 \in \mathbb{R}$, $u_1 \neq u_2$. L'équation (2.73) écrite dans le sens des distributions donne

$$ T_y'' = T_u + k(u_2 - u_1)\delta_{t_1}. $$

Si $u_1 \neq u_2$, le modèle est identifiable grâce à cette nouvelle relation.

Revenons au système $(\Gamma_a)$ pour lequel la solution $\zeta$ est différentiable par morceaux en temps. Soit $s$ la dérivée la plus élevée en temps de $y_a$ apparaissant dans (2.71). Si la fonction $y_a^{(s)}$ peut s'exprimer formellement dans (2.71) comme un terme de la forme $d(y_a, y_a', \dots, y_a^{s-1})y_a^{(s)}$ où $d$ est une fonction polynômiale différentiable en $y_a$ et qui ne s'annule pas, alors (2.71) peut se ré-écrire :

$$ y_a^{(s)} = \tilde{m}_0(y_a) + \sum_{k=1}^{d} \gamma_k(\theta_a) \tilde{m}_k(y_a), \qquad (2.74) $$

où, pour tout $k=0, \dots, d$, $\tilde{m}_k(y_a)$ est une fonction rationnelle supposée localement lebesgue intégrable en temps$^8$.

Supposons que $y_a$ est $s$ fois différentiable par morceaux par rapport au temps et soit $(t_\nu)_\nu$ la suite ordonnée de points correspondant aux points de discontinuité de $y_a$. Étant donné les discontinuités de $y_a$, la relation (2.74) peut être considérée au sens des distributions :

$$ T_{y_a}^{(s)} = T_{\tilde{m}_0(y_a)} + \sum_{k=1}^{n} \gamma_k(\theta_a) T_{\tilde{m}_k(y_a)}. \qquad (2.75) $$

Notons pour tout $j=0, \dots, s$, le saut en $t_\nu$ de $y_a^{(j)}$ est donné par

$$ \sigma_{y_a, \nu}^{j}(\theta_a) := y_a^{(j)}(t_{\nu}^{+}) - y_a^{(j)}(t_{\nu}^{-}). $$

Remarquons que le saut dépend de $\theta_a$ puisque $y_a$ dépend lui-même de $\theta_a$. Nous obtenons

$$ T_{y_a}^{(s)} = T_{y_a^{(s)}} + \sum_{\nu} \sum_{j=0}^{s-1} \sigma_{y_a, \nu}^{j} \delta_{t_{\nu}}^{(s-1-j)}, $$

où $\delta_{t_\nu}$ représente la distribution de Dirac au point $t_\nu$. Ainsi, l'équation (2.75) devient

$$ T_{y_a}^{(s)} = T_{\tilde{m}_0(y_a)} - \sum_{\nu} \sum_{j=0}^{s-1} \sigma_{y_a, \nu}^{j}(\theta_a) \delta_{t_\nu}^{(s-1-j)} + \sum_{k=1}^{d} \gamma_k(\theta_a) T_{\tilde{m}_k(y_a)}. \qquad (2.76) $$

### Remarque 16.

1. Ces relations sont plus informatives que (2.74) car les contraintes dues aux discontinuités et pouvant dépendre des paramètres à identifier sont prises en compte directement dans l'étude.

2. Cette opération peut permettre d'obtenir des relations bien définies par rapport à la régularité de la solution.

La proposition suivante donne une condition suffisante pour obtenir l'identifiabilité des paramètres.

<footnote n="8">8. La dépendance en temps a été omise pour simplifier les notations.</footnote>
**Proposition 18.** Supposons que la famille de distributions $(T_{\tilde{m}_k(y_a)})_{1\le k\le d}$ est linéairement indépendante. Si, pour tout $(\theta_a, \tilde{\theta}_a) \in L^\infty(\Omega, \mathcal{U}_p)^2$, nous avons

$$ \left\{ \begin{array}{l} \gamma_k(\theta_a) = \gamma_k(\tilde{\theta}_a) \\ \forall \nu, \quad \sum_{j=0}^{s-1} \sigma_{y_a, \nu}^j(\theta_a) = \sum_{j=0}^{s-1} \sigma_{y_a, \nu}^j(\tilde{\theta}_a) \quad \implies \theta_a = \tilde{\theta}_a, \end{array} \right. $$

alors le modèle est identifiable en $\theta_a$.

Le résultat vient de l'indépendance linéaire des familles $(\delta_{t_\nu}^{(j)})_{j=0,\dots,s}$ et $(\tilde{T}_{m_k(y_\theta)})_{k=1,\dots,n}$.

Ces relations différentielles (2.71), (2.72), (2.76) peuvent également être utilisées pour obtenir une première estimation des paramètres inconnus comme nous le verrons dans la section suivante.

### Résultat d'identifiabilité pour $(\Gamma_b)$

Soit le vecteur de paramètres $\theta_b := (\theta_c, \theta_d)$ où $\theta_c \in \mathbb{R}^{q_0}$ et $\theta_d := (d_i)_{i \in I}$ avec $I \subset \{q_0+1, \dots, q\}$. Soit $U := (U_i)_{i=1,\dots,q}$ la solution de $(\Gamma_b)$ et $y_b$ donné par (2.70). Nous avons obtenu le premier lemme suivant :

**Lemme 1.** Soit $y_b$ défini par $y_b := (y_c, y_d)$ où $y_d := (U_i)_{i \in I}$. Supposons qu'il existe $\nu \in \mathbb{N}$ tel que, pour tout $t \in (t_\nu, t_{\nu+1})$,

1. il existe une famille de réels positifs $(\lambda_i)_{i \in I}$ telle que

$$ \sum_{i \in I} \lambda_i (g_i(\zeta(t), U(t), \theta_d) - \partial_t U_i) $$

est indépendant de $\theta_d$,

2. la famille $(U_i(t, \cdot))_{i \in I}$ est linéairement indépendante à une constante additive près dans $\Omega$,

3. les $b_i$ définis dans $B_i = d_i A + b_i$ sont supposés être connus.

Alors $(\Gamma_b)$ est identifiable en $\theta_d$.

**Remarque 17.** Le lemme 1 donne une condition suffisante pour obtenir l'identifiabilité des coefficients de diffusion $(d_i)_{i \in I}$ quand les solutions associées $(U_i)_{i \in I}$ sont connues.

En effet, soit $\theta_d, \tilde{\theta}_d$ tel que pour tout $(t,x) \in (t_\nu, t_{\nu+1}) \times \Omega$, $y_b(t,x, \theta_d) = y_b(t,x, \tilde{\theta}_d)$. D'après l'hypothèse (1), nous avons

$$ \begin{align*} \sum_{i \in I} \lambda_i B_i(\theta_d) U_i &= \sum_{i \in I} \lambda_i (g_i(\zeta(t), U(t), \theta_d) - \partial_t U_i) \\ &= \sum_{i \in I} \lambda_i (g_i(\zeta(t), U(t), \tilde{\theta}_d) - \partial_t U_i) \\ &= \sum_{i \in I} \lambda_i B_i(\tilde{\theta}_d) U_i, \end{align*} $$

qui implique que

$$ \sum_{i \in I} \lambda_i (d_i - \tilde{d}_i) A U_i = 0. $$

Donc $V := \sum_{i \in I} \lambda_i (d_i - \tilde{d}_i) U_i$ est un vecteur propre de A associé à la valeur propre nulle. Ainsi, $V$ est constant. Comme pour tout $t \in (t_\nu, t_{\nu+1})$, $(U_i(t, \cdot))_{i \in I}$ est linéairement indépendant à une constante additive près dans $\Omega$, nous en déduisons que pour tout $i \in I$

$$ \lambda_i (d_i - \tilde{d}_i) = 0, $$

qui implique le résultat puisque $\lambda_i \neq 0$.

Pour obtenir l'identifiabilité de $\theta_c$, il est possible d'utiliser les propositions 16 et 17. Cependant, le calcul du wronskien pour prouver l'indépendance linéaire de la famille $(m_k(y_b))_{k=1,\dots,n}$ peut être difficile. Dans le lemme suivant, nous proposons de substituer ce test par un autre.

**Lemme 2.** Soit $\tilde{P} := \langle P \rangle$. Supposons que

1. $(\Gamma_b)$ est identifiable en $\theta_d$,
2. le polynôme P défini par (2.71) avec $\theta := \theta_c$ satisfait, pour tout $k = 1, \dots, n$,

$$ \gamma_k(\theta_c) = \gamma_k(\tilde{\theta}_c) \implies \theta_c = \tilde{\theta}_c, $$

3. la famille de fonctions ($\langle m_k(y_b) \rangle$)$_{k=1,...,n}$ est linéairement indépendante.
Alors ($\Gamma_b$) est identifiable en $\theta_c$.

La preuve est basée sur l'expression obtenue en faisant la différence entre $\tilde{P}(y_b, \theta_c)$ et $\tilde{P}(y_b, \tilde{\theta}_c)$ évaluée en $t_j$ où $y_b := y_{\theta_b}(\cdot, \theta_c)$ et $\tilde{y}_b := y_{\tilde{\theta}_b}(\cdot, \tilde{\theta}_c)$.

**Remarque 18.** 1. Si la relation $\tilde{P}$ obtenue en prenant la valeur moyenne de (2.71) n'est pas bien définie par rapport à la variable de temps, la relation peut être considérée dans le sens des distributions. Dans ce cas, pour l'hypothèse 3, il est suffisant de prouver l'indépendance linéaire de la famille $\langle m_k(y_b) \rangle$)$_{k=1,...,n}$. En effet, si la famille $\langle m_k(y_b) \rangle$)$_{k=1,...,n}$ est linéairement indépendante alors la famille $T_{\langle m_k(y_b) \rangle}$$_{k=1,...,n}$ est également linéairement indépendante.

2. Une autre façon de traiter l'hypothèse d'indépendance linéaire est de supposer qu'il existe n points distincts $t_1, \dots, t_n$ tel que le déterminant

$$ \left| \begin{array}{ccc} \langle m_1(y_b) \rangle (t_1) & \dots & \langle m_n(y_b) \rangle (t_1) \\ \langle m_1(y_b) \rangle (t_2) & \dots & \langle m_n(y_b) \rangle (t_2) \\ \vdots & & \vdots \\ \langle m_1(y_b) \rangle (t_n) & \dots & \langle m_n(y_b) \rangle (t_n) \end{array} \right| \qquad (2.77) $$

est non nul.

Dans les applications numériques, les n points distincts $t_1, \dots, t_n$ peuvent être choisis de façon aléatoire pour évaluer ce déterminant.

### 2.6.5 Exemple

Reprenons l'exemple 2.3.4. Contrairement au déplacement des moustiques qui est très souvent limité à quelques mètres, nous modélisons la mobilité humaine en ajoutant des termes de diffusion comme dans [149]. Le modèle se ré-écrit :

$$ \left\{ \begin{array}{l} \left( \Gamma_a \right) \left\{ \begin{array}{l} \partial_t E_m(t,x) = bA_m(t,x)\left(1 - \frac{E_m(t,x)}{K_E(t,x)}\right) - (s+d)E_m(t,x), \\ \\ \partial_t L_m(t,x) = sE_m(t,x)\left(1 - \frac{L_m(t,x)}{K_L(t,x)}\right) - (s_L+d_L)L_m(t,x), \\ \\ \partial_t A_m(t,x) = s_L L_m(t,x) - d_m A_m(t,x), \\ \\ \partial_t I_m(t,x) = -(s_L \frac{L_m(t,x)}{A_m(t,x)} + \beta_m I_H(t,x))I_m(t,x) + \beta_m I_H(t,x), \end{array} \right. \\ \\ \left( \Gamma_b \right) \left\{ \begin{array}{l} \partial_t S_H(t,x) = -(b_H + \beta_H I_m(t,x))S_H(t,x) + b_H + d_1 \Delta S_H(t,x), \quad (b) \\ \\ \partial_t I_H(t,x) = \beta_H I_m(t,x)S_H(t,x) - (\gamma + b_H)I_H(t,x) + d_2 \Delta I_H(t,x), \end{array} \right. \end{array} \right. \qquad (2.78) $$

où $d_1, d_2 > 0$ sont des coefficients de diffusion, $x$ est dans le domaine borné $\Omega \subset \mathbb{R}^n$ avec des conditions $C^2$ au bord et $t > 0$. Le modèle est complété avec les conditions initiales suivantes

$$ \begin{aligned} & (E_m(0,x), L_m(0,x), A_m(0,x)) = (E_0(x), L_0(x), A_0(x)), \\ & (I_m(0,x), S_H(0,x), I_H(0,x)) = (I_0^m(x), S_0^H(x), I_0^H(x)). \end{aligned} $$

De plus, nous supposons qu'il n'existe pas de population traversant la frontière $\partial\Omega$ de $\Omega$. Ainsi, le modèle est complété par les conditions de Neumann suivantes :

$$ \frac{\partial S_H}{\partial \nu} = \frac{\partial I_H}{\partial \nu} = 0 \text{ dans } \mathbb{R}_+ \times \partial\Omega, $$

où $\nu$ est la normale sortante de $\partial\Omega$. Nous supposons également que
1. $(E_0, L_0, A_0) \in C(\bar{\Omega})^3, I_0^m \in C(\Omega)$ et $(S_0^H, I_0^H) \in L^2(\Omega)^2$.

2. Les capacités d'accueil de gîte $K_E$ et $K_L$ satisfont

$$K_E, K_L \in L^\infty([0, \infty) \times \Omega) \text{ avec } K_L(t,x) = K_E(t,x)/2 \neq 0, \text{ pour presque tout } (t,x) \in \mathbb{R}_+ \times \Omega.$$

3. $K_E$ et $K_L$ sont des fonctions constantes par morceaux en temps et connues, ce qui peut correspondre à des changements brutaux de l'environnement.

En posant, pour $x$ fixé dans $\Omega$,

$$\forall t \in [0, \infty), \quad \zeta_1(t) := \frac{E_m(t, x)}{K_E(t, x)}, \quad \zeta_2(t) := \frac{L_m(t, x)}{K_E(t, x)}, \quad \zeta_3(t) := \frac{A_m(t, x)}{K_E(t, x)},$$

Comme $\tilde{\zeta}_2$ dépend seulement de $d_m$ et de la variable mesurée $\zeta_2$, nous la considérons comme variable mesurée. Les trois observations sont

$$y_{a,1} := \zeta_2, \quad y_{a,2} := \tilde{\zeta}_2 \text{ et } y_{a,3} := \lambda.$$

Le système (2.78)-(a) se ré-écrit :

$$\left\{
\begin{aligned}
\partial_t \zeta_1 &= b\zeta_3(1-\zeta_1)-(s+d)\zeta_1 \\
\partial_t \zeta_2 &= s\zeta_1(1-2\zeta_2)-(s_L+d_L)\zeta_2 \\
\zeta_3 &= \zeta_3(0)\lambda + s_L \tilde{\zeta}_2 \\
\partial_t \tilde{\zeta}_2 &= \zeta_2 - d_m \tilde{\zeta}_2 \\
\partial_t \lambda &= -d_m \lambda \\
y_{a,1} &= \zeta_2, \quad y_{a,2} = \tilde{\zeta}_2, \quad y_{a,3} = \lambda
\end{aligned}
\right.
\qquad (2.79)$$

avec les conditions initiales :

$$(\zeta_1(0,x), \zeta_2(0,x), \zeta_3(0,x), \tilde{\zeta}_2(0,x), \lambda(0)) = (\zeta_{1,0}(x), \zeta_{2,0}(x), \zeta_{3,0}(x), 0, 1).$$

Soit

$$\Sigma := \left\{ (E, L, A) \in C(\bar{\Omega})^3 \mid 0 \le E \le k, 0 \le L \le \frac{k}{2}, 0 \le A \le \frac{s_L k}{2d_m} \right\},$$

où nous avons posé $k := \|K_E\|_\infty > 0$.

Une simple adaptation des Lemmes 4.2 et 4.3 de [145] permet de passer du cas $K_E$ constant au cas où $K_E$ est une fonction dépendant de $x$ et constante par morceaux par rapport au temps, d'où :

**Proposition 19.** Soit $(E_0, L_0, A_0) \in \Sigma$ et

$$R := \frac{s_L}{d_m} \frac{b}{s+d} \frac{s}{s_L+d_L}.$$

alors, si $R > 1$, le système (2.78)-(a) admet une unique solution $(E_m, L_m, A_m)$ dans l'espace des fonctions différentiables par morceaux $C_l^1([0, \infty); \text{int}(\Sigma))^3$ où $\text{int}(\Sigma)$ est l'intérieur de $\Sigma$.

**Remarque 19.**

1. Les points de discontinuité de $E_m, L_m$ et $A_m$ sont les mêmes que ceux de $K_E$.

2. De [145], la condition de survie des oeufs, larves et des populations de moustiques adultes est donnée par $R > 1$. Nous supposerons par la suite que $R > 1$ c'est à dire que $E_m, L_m, A_m > 0$. Nous faisons en particulier l'hypothèse suivante :

$$(\mathcal{H}) \quad R > 1 \text{ et il existe } c > 0 \text{ tel que, pour tout } t \ge 0, \zeta_3(t) := A_m(t, \cdot) > c.$$

Comme $(\mathcal{H})$ est vérifiée,

$$r := \frac{L_m}{A_m} \in C_l^1([0, \infty); C(\bar{\Omega})).$$

En posant $U := (I_m, S_H, I_H)$, le système (2.78)-(b) est de la forme (2.66) où l'espace de Banach $X$ est

$$X := \{(u_1, u_2, u_3) \mid u_1 \in C(\bar{\Omega}), u_2, u_3 \in L^2(\Omega)\},$$

et l'opérateur $B$ est défini par

$$B := \text{diag}(I, -d_1 A + b_H, -d_2 A + \gamma + b_H).$$
A est la réalisation du Laplacien dans $L^2(\Omega)$ avec les conditions de Neumann homogènes sur le bord $\partial\Omega$, et $I$ est l'opérateur identité sur $C(\bar{\Omega})$. Par conséquent, l'opérateur non-linéaire $G(t, \cdot)$ est défini par

$$G(t, U) := -(s_L r(t) + \beta_m u_3 - 1)u_1 + \beta_m u_3, -\beta_H u_1 u_3 + b_H, \beta_H u_1 u_2).$$

Nous avons montré dans [200] que $G$ vérifie (2.67) et ainsi que (2.78)-(b) admet une unique solution locale

$$U \in C_l([0, T_{U_0}]; D(B)) \cap C_l^1((0, T_{U_0}]; X)$$

où $U_0 := (I_0^m, S_0^H, I_H^0)$.

De plus, si

$$X_+ := \{(u_1, u_2, u_3) \in X \mid u_1, u_2, u_3 \ge 0\}, \quad (2.80)$$

le théorème suivant donné dans [200] permet de conclure que la solution est globale.

**Théorème 1.** Sous l'hypothèse $(\mathcal{H})$, pour toute $(I_0^m, S_0^H, I_H^0) \in X_+$, (2.78)-(b) possède une unique solution globale positive telle que

$$\left\{ \begin{array}{l} I_m \in C_l^1([0, \infty); C(\bar{\Omega})), \\ S_H, I_H \in C_l([0, \infty); H_N^2(\Omega)) \cap C_l^1((0, \infty); L^2(\Omega)). \end{array} \right.$$

### 2.6.5.1 Résultats d'identifiabilité

Nous reprenons les hypothèses de la section 2.4.2.4, c'est à dire que les paramètres à identifier sont $s$, $s_L$ pour le système (2.78)-(a) et $\beta_H$, $\beta_m$, $d_1$, $d_2$ pour le système (2.78)-(b). Les sorties mesurées sont $L_m$, $S_H$, $I_H$.

Dans [229], l'identifiabilité du modèle (2.78) a été étudiée mais dans le cas de fonction régulières : pour le système (2.78)-(a), les solutions étaient $C^\infty$ et pour le système (2.78)-(b), elles étaient continues.

Nous avons démontré le résultat suivant dans [200] :

**Théorème 2.** Supposons que les constantes $b, d, d_L, d_m, b_H, \gamma$ sont connues et qu'il existe $\nu \in \mathbb{N}$ tel que

- pour tout $(t, x) \in (t_\nu, t_{\nu+1}) \times \Omega$, $L_m(t, x)$, $S_H(t, x)$ et $I_H(t, x)$ sont des sorties du modèle,

- pour chaque $t \in (t_\nu, t_{\nu+1})$, $S_H(t, \cdot)$ et $I_H(t, \cdot)$ sont indépendants linéairement à une constante additive près dans $\Omega$.

Nous rappelons que la notation $\langle \cdot, \cdot \rangle$ désigne la moyenne sur tout l'espace $\Omega$. Alors, le système (2.78) est identifiable en $s, s_L, A_0, d_1$ et $d_2$.

Soit $\tilde{I}_m := \beta_H I_m$. Sous les mêmes hypothèses, si il existe $\tilde{t} \in (t_\nu, t_{\nu+1})$ tel que $\langle I_H S_H^2 \rangle (\tilde{t}) \neq 0$ et $\langle \tilde{I}_m I_H S_H^2 \rangle (\tilde{t}) \neq 0$, alors le système (2.78) est identifiable en $\beta_m$ et $\beta_H$.

### Remarque 20.

1. Dans la démonstration, un résultat supplémentaire a été obtenu sur l'identifiabilité de la condition condition initiale $A_0 \in C(\bar{\Omega})$. Par conséquent, on peut mettre en place une procédure d'estimation de paramètres pour estimer de façon unique $A_0$ et les autres paramètres à partir de $L_m$, $S_H$ et $I_H$.

2. L'identifiabilité des paramètres $\beta_m$ et $\beta_H$ a été démontrée à partir d'un polynôme ES obtenu en appliquant l'algorithme de Rosenfeld-Groebner au système (2.78)-(b).

3. L'identifiabilité des paramètres $d_1$ et $d_2$ a été obtenue à partir d'un polynôme ES qui est repris à la section suivante pour estimer ces paramètres.

### 2.6.5.2 Exemple de l'estimation de paramètres à partir des polynômes ES dans le cadre des EDP

Dans cette section, nous présentons une méthode numérique basée sur les polynômes ES pour estimer les paramètres $d_1$ et $d_2$.

La relation ES est obtenue en faisant la somme de la seconde et la troisième équation de (2.78)-(b), ce qui nous donne

$$d_1 \Delta y_{b,1}(t, \cdot) + d_2 \Delta y_{b,2}(t, \cdot) = \partial_t y_{b,1}(t, \cdot) + \partial_t y_{b,2}(t, \cdot) + b_H (y_{b,1}(t, \cdot) - 1) + (\gamma + b_H) y_{b,2}(t, \cdot). \quad (2.81)$$

Pour obtenir une relation bien définie, la relation (2.81) est ré-écrite au sens des distributions :

$$d_1 \int_{\Omega} S_H(t, x) \Delta \psi(x) dx + d_2 \int_{\Omega} I_H(t, x) \Delta \psi(x) dx = \int_{\Omega} (\partial_t S_H(t, x) + \partial_t I_H(t, x)) \psi(x) dx \\ + (\gamma + b_H) \int_{\Omega} I_H(t, x) \psi(x) dx - b_H \int_{\Omega} (1 - S_H(t, x)) \psi(x) dx,$$
pour tout $\psi \in C_c(\Omega)$. En intégrant la relation sur $(0,t)$, nous obtenons

$$
\begin{align}
d_1 \int_0^t \left( \int_{\Omega} S_H(\tau, s) \Delta \psi(s) ds \right) d\tau + d_2 \int_0^t \left( \int_{\Omega} I_H(\tau, s) \Delta \psi(s) ds \right) d\tau \nonumber \\
&= \int_{\Omega} (S_H(t, s) + I_H(t, s)) \psi(s) ds - \int_{\Omega} (S_H(0, s) + I_H(0, s)) \psi(s) ds \tag{2.82} \\
&\quad + (\gamma + b_H) \int_0^t \left( \int_{\Omega} I_H(\tau, s) \psi(s) ds \right) d\tau - b_H \int_0^t \left( \int_{\Omega} (1 - S_H(\tau, s)) \psi(s) ds \right) d\tau, \nonumber
\end{align}
$$

pour tout $\psi \in C_c(\Omega)$ et tout $t > 0$.

Pour les fonctions test, nous avons choisi :

$$
\psi(x,y) := \begin{cases} e^{\frac{1}{1-x^2-y^2}} & \text{if } (x,y) \in \Omega, \\ 0 & \text{else,} \end{cases}
$$

et une discrétisation $0 = t_1 < t_2 < \dots < t_{M+1} = T$ de $[0, T]$. Pour $k = 1, 2, \dots, M$, nous posons

$$
A(k, 1) := \int_{0}^{t_{k+1}} \left( \int_{\Omega} S_{H}(\tau, s) \Delta \psi(s) \, ds \right) d\tau,
$$

$$
A(k, 2) := \int_{0}^{t_{k+1}} \left( \int_{\Omega} I_{H}(\tau, s) \Delta \psi(s) \, ds \right) d\tau,
$$

et

$$
\begin{align*}
B(k) &:= \int_{\Omega} (S_{H}(t_{k+1}, s) + I_{H}(t_{k+1}, s)) \psi(s) ds - \int_{\Omega} (S_{H}(0, s) + I_{H}(0, s)) \psi(s) ds \\
&\quad + (\gamma + b_{H}) \int_{0}^{t_{k+1}} \left( \int_{\Omega} I_{H}(\tau, s) \psi(s) ds \right) d\tau - b_{H} \int_{0}^{t_{k+1}} \left( \int_{\Omega} (1 - S_{H}(\tau, s)) \psi(s) ds \right) d\tau.
\end{align*}
$$

De (2.82), on déduit le système suivant :

$$
Ad = B, \text{ où } d := (d_1, d_2)^T \tag{2.83}
$$

qui est résolu en utilisant la méthode des moindres-carrés basée sur la factorisation QR. Les résultats sont résumés dans le tableau 2.5 où $\lambda$ est l'erreur relative entre $(S_H, I_H)$ et $(S_H^{meas}, I_H^{meas})$, tandis que $re$ est l'erreur relative entre $(d_1, d_2)$ et $(d_1^a, d_2^a)$, la solution de (2.83).

Choix des paramètres pour les simulations : Les paramètres utilisés pour les simulations sont données dans le tableau 1 de [146] (également dans [147]) et sont rappelés dans le tableau 2.2. Le domaine considéré est $\Omega := (0,1)^2$ et le temps final $T = 20$. Nous faisons également l'hypothèse que la population susceptible diffuse plus vite que la population infectée dans le domaine en prenant $d_1 = 0.1$ et $d_2 = 0.01$.

Au temps initial, les variables $\zeta_i$ sont constantes. La population infectée est supposée concentrée au centre du domaine, ce qui pourrait modéliser un hôpital par exemple et la population susceptible est localisée aux extrémités du domaine. Les conditions initiales sont alors :

$$
\left\{
\begin{array}{l}
\zeta_1(t=0, x, y) = 0.1, \quad \zeta_2(t=0, x, y) = 0.04, \quad \zeta_3(t=0, x, y) = 0.01, \\
I_m(t=0, x, y) = 0.2, \\
S_H(t=0, x, y) = 1 - z(x, y), \\
I_H(t=0, x, y) = z(x, y),
\end{array}
\right.
$$

où

$$
z(x,y) := \begin{cases} e^{-\frac{1}{4-(x-\frac{1}{2})^2-(y-\frac{1}{2})^2}}, & \text{if } (x-\frac{1}{2})^2+(y-\frac{1}{2})^2 < \frac{1}{4} \\ 0 & \text{sinon.} \end{cases}
$$

Les simulations de $S_H$ et $I_H$ sont obtenues en résolvant le système (2.78) par la méthode des différences finies explicites.
Les états mesurés $S_H^{meas}$ et $I_H^{meas}$ sont construits en utilisant la fonction `randn` de Matlab dont le coefficient $\sigma$ est choisi pour avoir une valeur maximale de l'erreur relative entre $S_H$ et $S_H^{meas}$ (et entre $I_H$ et $I_H^{meas}$) égale à 1%, 5%, 10% et 20%.

Résultats numériques :

<table><caption>TABLE 2.5 – Résultats de la résolution du système (2.83) par la méthode des moindres-carrés.</caption><thead><tr><th>λ</th><th>1%</th><th>5%</th><th>10%</th><th>20%</th></tr></thead><tbody><tr><td>d<sub>1</sub><sup>a</sup></td><td>0.1019</td><td>0.1011</td><td>0.0836</td><td>0.1153</td></tr><tr><td>d<sub>2</sub><sup>a</sup></td><td>0.0186</td><td>0.0159</td><td>0.0102</td><td>0.0191</td></tr><tr><td>re</td><td>0.0880</td><td>0.0595</td><td>0.1636</td><td>0.1773</td></tr></tbody></table>

Cette méthode donne de bonnes approximations de $d_1$ et $d_2$ et montre une stabilité à un bruit fort.

## 2.7 Conclusion

Une partie du travail qui a suivi ma thèse a consisté à élargir la notion d'identifiabilité à d'autres cadres mathématiques. Dans le même temps, il m'a paru important de concevoir et de développer des méthodes pour démontrer cette propriété.

Tout d'abord, la notion d'identifiabilité a été généralisée à l'identifiabilité relative qui étudie l'identifiabilité de certains paramètres par rapport à d'autres. Cette propriété permet de connaître les paramètres clés à déterminer pour obtenir l'identifiabilité de paramètres plus difficiles à estimer voir l'identifiabilité du système complet, le résultat étant certifié grâce aux outils semi-algébriques. Puis, ce concept a été généralisé aux systèmes dans le cadre d'erreurs bornées. Un algorithme a été développé pour partitionner le domaine des paramètres en parties présentant des propriétés différentes d'identifiabilité dues au contexte ensembliste. Enfin, l'identifiabilité ainsi que la méthode basée sur la théorie de l'élimination ont été généralisées au cadre EDP montrant l'intérêt qu'elles peuvent présenter pour certains types de modèle. Ce travail a été complété par l'étude et le développement de méthodes numériques pour estimer les paramètres dans chacun de ces trois cadres.

Dans le chapitre suivant, je présente une application de cette propriété à un autre domaine, celui du diagnostic et le lien existant entre ces deux notions. Une des conséquences est le développement de méthodes en diagnostic dans le cas simple et multi-fautes. L'identifiabilité sera vue également au chapitre 5 dans un cas pratique, celui de la modélisation des comportements en situation de catastrophe.
<figure><img src="image_23.png" /></figure>