Main

Télécharger au format pdf ou txt
Télécharger au format pdf ou txt
Vous êtes sur la page 1sur 71

République Algérienne Démocratique et Populaire

Ministère de L’Enseignement Supérieur et de la Recherche Scientifique

Université Ferhat Abbas Sétif 1

Département d’Electrotechnique de la Faculté de Technologie

Polycopié de Cours
Diagnostic et Surveillance des Systèmes

Réalisé par : SID Mohamed Amine

Année Universitaire 2016-2017


Avant propos

E cours s’adresse aux étudiants de Master en Automatique, ainsi qu’aux élèves ingénieurs et
C plus généralement aux étudiants sollicités par des problèmes de diagnostic et détection de
défaut d’un système linéaire. Il devrait leur permettre de maitriser les concepts et les méthodes
de diagnostic afin de les appliquer à des domaines aussi variés que l’électricité, la mécanique, la
pétrochimie etc. Les pré-requis constitués des connaissances mathématiques dispensées durant les
deux premières années d’un cursus universitaire dans la filière "technologie" ou dans le cadre des
classes préparatoires aux écoles supérieures, pourront être complétés par les notions basiques du
traitement des signaux stochastique et de la notion de la robustesse en automatique.
Le contenu de ce cours, que j’assure son enseignement pour les étudiants de M2 au département
d’électrotechnique de l’université de Sétif1, a été rédigé dans le souci permanent d’éviter la pesante
et souvent l’inefficace présentation des notions mathématiques avancées à de rares exceptions des
démonstrations techniques relativement aux conditions d’applicabilité des méthodes proposées ;
en revanche certaines démonstrations abordables et ayant un intérêt pédagogique sont proposées
sous forme d’exercices. Cependant, à la fin de chaque chapitres, des problèmes de simulations sous
l’environnement MATLAB/SIMULINK sont proposés afin d’illustrer l’efficacité des méthodes
enseignées et d’aider le lecteur à bien assimiler les connaissances fournies à travers un outil
informatique interactif.

U NIVERSITÉ DE S ÉTIF .1

MOHAMED . SID @ UNIV- SETIF. DZ


Table des matières

I Notions et méthodes basiques du Diagnostic

1 Introduction au diagnostic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1 Intérêt pratique du diagnostic des défauts 7
1.2 Terminologie de base 10
1.2.1 Classification des défauts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Modélisation des défauts et des systèmes dynamiques 11
1.3.1 Modélisation du comportement nominal d’un système . . . . . . . . . . . . . . . . . . 12
1.3.2 Les systèmes perturbés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.3 Modélisation des défauts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4 Classification des méthodes de diagnostic 15
1.4.1 Méthodes à base de données (signal) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.4.2 Méthodes à base de modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.4.3 Performances d’un algorithme de diagnostic : . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.5 Exercices d’application 18

2 Espace de parité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1 Introduction 21
2.2 Sous-epaces vectoriels nuls 21
2.3 Redondance statique 23
2.4 Redondance dynamique 27
2.4.1 Auto-redondance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.2 Inter-redondance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5 Exercices d’application 30
3 Diagnostic à base d’observateurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.1 Introduction 33
3.2 Observateur de Luenberger 33
3.3 Exercises d’application 34

II Techniques de diagnostic robustes

4 Observateur à entrée inconnue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41


4.1 Introduction 41
4.2 Décoplage de l’entrée inconnue 41
4.3 Structuration de résidu 43
4.4 Exercices d’application 44

5 Placement des structures propres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47


5.1 Introduction 47
5.2 Synthèse de l’observateur 47
5.2.1 Structuration de résidu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.2.2 Placement des structures propres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3 Exercices d’application 50

6 Le filtre de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.1 Introduction 53
6.2 Estimation d’état d’un système stochastique 53
6.3 Le système linéaire stochastique 55
6.4 Le filtre de KALMAN 56
6.4.1 L’approche indirecte (Bayésienne) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.4.2 L’approche directe (EQMM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.5 Application au diagnostic 61
6.5.1 Conditions de convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.6 Exercices d’application 63

A Appendice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
A.1 La somme de deux Gaussiennes 69
A.2 Dérivée matricielle 70

Références bibliographiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Livres 71
Articles 71
I
Notions et méthodes
basiques du Diagnostic

1 Introduction au diagnostic . . . . . . . . . . . . 7
1.1 Intérêt pratique du diagnostic des défauts
1.2 Terminologie de base
1.3 Modélisation des défauts et des systèmes dyna-
miques
1.4 Classification des méthodes de diagnostic
1.5 Exercices d’application

2 Espace de parité . . . . . . . . . . . . . . . . . . . . 21
2.1 Introduction
2.2 Sous-epaces vectoriels nuls
2.3 Redondance statique
2.4 Redondance dynamique
2.5 Exercices d’application

3 Diagnostic à base d’observateurs . . . . 33


3.1 Introduction
3.2 Observateur de Luenberger
3.3 Exercises d’application
1. Introduction au diagnostic

1.1 Intérêt pratique du diagnostic des défauts


Le demande croissante d’augmentation des performances des systèmes techniques et d’amé-
lioration de la qualité de produit en tenant compte des coûts de production, a considérablement
augmenté le degré de complexité des procédés automatisés. Cependant, ce développement croissant
nécessite un renforcement des conditions de sécurité et de fiabilité des systèmes. Actuellement,
l’un des problèmes les plus critiques dans la conception des systèmes automatiques est comment
garantir leurs Sûreté de fonctionnement et leurs fiabilité. Une façon traditionnelle permettant d’as-
surer la fiabilité et la Sûreté d’un système est d’améliorer la qualité, la fiabilité et la robustesse
des composants (éléments) constituant le système tels que : les capteurs, les actionneurs et les
contrôleurs. Néanmoins, cette dernière méthode ne peut garantir un fonctionnement sans défaut
du système global. La surveillance de processus et le diagnostic de défauts sont donc de plus
en plus un composant indispensable d’un système automatique moderne et souvent cadré par
des lois réglementaires et des législations spécifiques. Initiée dans les années 1970, la technique
de diagnostic à base de modèle (sujet de ce cours) a remarquablement évolué en intégrant les
résultats de la théorie de commande (notamment l’estimation) dans la conception des algorithmes
de diagnostic des défauts. Son efficacité dans la détection de défauts dans un système dynamique
a été démontrée par un grand nombre d’applications réussies dans les processus industriels et les
systèmes de contrôle automatique. Aujourd’hui, les algorithmes de diagnostic à base de modèle sont
omniprésents intégrés dans les systèmes de contrôle de véhicules, robots, systèmes de transport, les
réseaux électriques, les procédés manufacturiers, etc.

Les informations fournies par l’algorithme du diagnostic sont indispensables pour lancer l’une des
actions correctives suivantes :
— Maintenance prévisionnelle : est une maintenance conditionnelle basée sur le franchis-
sement d’un seuil prédéfini qui permet de donner l’état de dégradation du bien avant sa
détérioration complète (défaillance). L’objectif de la maintenance préventive est de déter-
miner l’ensemble des actions à exercer sur le procédé afin de ne pas subir l’effet d’une
défaillance. On peut a cet effet distinguer deux approches possibles : la maintenance préven-
8 Chapitre 1. Introduction au diagnostic

tive systématique et la maintenance préventive conditionnelle.


Dans le premier cas, les activités de maintenance sont planifiées et on lieu selon un échéancier
basé sur le temps ou l’unité d’usage. Lors de ces interventions, les éléments sont remplacés
même s’ils ne sont pas défaillants (voir F IGURE 1.1).

défauts

u(t) y(t)
Procédé

Durée Planning
Unité d’usage d’intervention Activités de

maintenance

F IGURE 1.1 – Maintenance préventive systématique.

Dans le deuxième cas, les activités de maintenance sont déclenchés en fonction d’information
reflétant l’état de dégradation de l’équipement considéré. dans ce cas, les éléments sont
remplacés que si nécessaire (voir F IGURE 1.2).

défauts

u(t) y(t)
Procédé

Système de
diagnostic Activités de

maintenance

F IGURE 1.2 – Maintenance préventive conditionnelle.

La première approche qui est aussi la plus ancienne, a fait l’objet d’une attention considérable
de la part de la communauté scientifique. sa mise en oeuvre est maintenant bien maîtrisée et
l’on dispose pour cela d’une méthodologie générale. La deuxième approche, plus récente, est
d’une mise en oeuvre pratique plus délicate ; ceci s’explique par la multiplicité des approches
possibles ainsi que par la diversité des informations nécessaire à son élaboration.
On s’intéresse dans la suite à la maintenance conditionnelle, basé sur la surveillance en
continu de l’évolution du système considéré. Ce type d’approche nécessite la conception d’un
système de diagnostic permettant la détection précoce de déviations faibles par rapport à une
caractérisation du système en fonctionnement nominal, afin de prévenir un dysfonctionnement
avant qu’il n’arrive (voir F IGURE 1.3).
Les grandeurs surveillées peuvent être diverses : un courant électrique, une température, une
pression, un pourcentage de particule dans l’huile, un niveau de vibration et bien d’autres
encore. Le suivi d’une grandeur donnée s’effectue en mesurant sa valeur régulièrement dans
le temps. La grandeur surveillée subit, au cours du fonctionnement de l’installation deux
1.1 Intérêt pratique du diagnostic des défauts 9

Grandeur mesurée : y

yn + Δy

yn

yn - Δy

yn : valeur nominale de y

t1 t2 Temps

F IGURE 1.3 – Déviation d’une grandeur mesurée par rapport à sa valeur nominale.

types de variation : des fluctuations autour de sa valeur nominale qui sont dues aux diverses
perturbations agissant sur le système considéré, une dérive qui peut être due à un phénomène
d’usure, un phénomène de dégradation progressive ou à des conditions d’environnement
anormales. Le suivi de cette grandeur permet alors de vérifier qu’elle ne s’écarte pas, de
façon significative, de sa valeur nominale. Dans le cas d’une déviation jugée anormale (étape
de détection de défaut), il s’agit de localiser l’origine de cette anomalie (étape de diagnostic
du défaut), puis de prendre, en conséquence, les mesures qui s’imposent pour un retour à la
normale du fonctionnement du système (étape de prise de décision). Ces différentes étapes
seront réalisées à l’aide de ce que nous appellerons dans la suite un système de diagnostic.
Un système de diagnostic doit donc être en mesure de réaliser les trois étapes essentielles
suivantes : la détection d’un défaut, le diagnostic du défaut, la prise de décision pour un
retour à la normale. Avant de poursuivre, il convient de préciser les termes de détection
et de diagnostic des défauts. Le diagnostic a pour objectif de rechercher l’origine d’un
défaut constaté. Un défaut correspond à une déviation jugée anormale d’une grandeur
caractéristique du système. Toute méthode de diagnostic repose sur l’analyse d’un certain
nombre d’indicateurs, ou symptômes, permettant de caractériser l’état de fonctionnement (de
santé) du système. Le résultat de cette analyse débouche, dans le cadre de la maintenance
conditionnelle, sur une prise de décision destinée à un retour à la normale de l’installation
jugée défaillante.
La mise en oeuvre d’une telle approche (voir F IGURE 1.4) nécessite de générer des indicateurs
de défauts ou symptômes, puis à interpréter correctement ces symptômes afin de déterminer
l’origine de défaut, c’est à dire l’élément présentant un fonctionnement anormal et enfin à
prendre une décision pour un retour à un fonctionnement normal de l’installation.

Système de Diagnostic

Prise de Interprétation Génération


des d’indicateurs
décision
symptômes de défauts

Défaut Acquisition des grandeurs

Actions sur le système caractéristiques du système

entrées sorties
Procédé Industriel

F IGURE 1.4 – Structure générale d’un système de diagnostic


10 Chapitre 1. Introduction au diagnostic

L’élaboration d’une structure telle que celle présentée à la F IGURE 1.4 nécessite l’exploitation
de l’ensemble des connaissances disponibles sur l’installation. Ce dernier point fait l’objet
du paragraphe suivant.
— Le pronostic : Le pronostic est une procédure permettant d’évaluer la dégradation future du
système et d’en déduire sa durée de vie résiduelle (délai avant sa défaillance) en se basant sur
des informations fournies par le bloc du diagnostic.
— La commande tolérante aux défauts : une technique de commande permettant à un système
subissant un défaut de continuer à fonctionner, éventuellement de manière réduite (mode
dégradé), au lieu de tomber complètement en défaillance.

1.2 Terminologie de base


Rappelons tout d’abord quelques éléments de terminologie utilisés dans le domaine du diag-
nostic adoptées par le comité technique SAFEPROCESS de l’IFAC (International Federation of
Automatic contol). On définit :
Défaut (fault) : une déviation non-autorisée d’au moins une propriété caractéristique ou un
paramètre du système à partir de la condition standard / acceptable. Il peut être modélisé
comme une entrée externe ou sous forme d’écart de paramètre qui modifie les caractéristiques
du système. Les incertitudes et les perturbations sont modélisées de manière assez similaire
aux défauts comme une déviation paramétrique ou une entrée externe.
Défaillance (failure) : une interruption permanente de la capacité d’un système à effectuer les
fonctions attendues dans les conditions de fonctionnement nominales. Il est clair qu’une
défaillance implique l’apparition d’un défaut vue existence d’un écart entre la caractéristique
théorique et mesurée. Par contre, un défaut n’implique pas nécessairement une défaillance
puisque le dispositif peut très bien continuer à assurer sa fonction principale.
Panne : une panne est l’inaptitude d’un dispositif à accomplir une fonction requise. Une panne
résulte toujours d’une défaillance et donc d’un défaut :
Défaut −→ Défaillance −→ Panne
Dans le cadre de la maintenance préventive conditionnelle, le diagnostic doit permettre de
détecter et de localiser un défaut avant que celui-ci ne conduise à une défaillance ou à une
panne qui entrainerait l’arrêt du système.
Détection des défauts (faults detection) : Détermination des défauts présent dans un système et
de leur instants de détection.
Diagnostic des défauts (faults diagnosis) : localisation et l’identification des défauts.
Localisation des défauts (faults isolation) : Détermination du type, de la localisation et des
instants de détection des défauts.
Identification des Défauts (faults identification) : Détermination de l’amplitude et du comporte-
ment temporel des défauts.
Résidu Un signal indicateur de défauts, obtenu à partir de la différence entre les mesures et les
calculs.

1.2.1 Classification des défauts


Les défauts peuvent être classés selon les critères détaillés ci-après.
Selon leur rapidité de manifestation
Les défauts peuvent être classés selon leurs dynamiques d’évolution en : défauts brusques,
progressifs et intermittents. Les défauts brusques (biais) se produisent instantanément souvent à
cause de dommages matériels. Habituellement ils sont très graves car ils affectent les performances
et/ou la stabilité du système commandé. Les défauts progressifs (dérives) sont dus au changements
lents des valeurs paramétriques. Souvent dus au vieillissement, Ils sont plus difficiles à détecter
1.3 Modélisation des défauts et des systèmes dynamiques 11

en raison de leur dynamique lente, mais sont également moins graves. Les défauts intermittents
(valeurs aberrantes) sont des défauts qui apparaissent et disparaissent à plusieurs reprises, par
exemple à cause d’un câblage partiellement endommagé.

F IGURE 1.5 – Classification de défaut selon leurs rapidité d’évolution

Selon leurs emplacements


Dans les systèmes de commande, les défauts sont classés selon l’élément affecté par ce défaut
en trois catégories :
— Défauts actionneur.
— Défauts capteur.
— Défauts procédé (système).

F IGURE 1.6 – Classification des défauts selon leurs emplacements

Pour la détection, la localisation et l’identification des défauts, on fait appel à un “module du


diagnostic” qui est un calculateur exécutant un algorithme de diagnostic. Cet algorithme est conçu
pour la surveillance du système en utilisant les mesures en provenance des capteurs et les entrées
de commandes. Après la détection du défaut l’utilisateur peut choisir entre la reconfiguration de
la loi de commande (commande tolérante aux défauts) ou l’arrêt total du fonctionnement afin de
procéder à la maintenance corrective.

1.3 Modélisation des défauts et des systèmes dynamiques


Un système physique est constitué d’un ensemble de composants interconnectés entre eux
afin d’assurer une fonction requise. Cependant, pour pouvoir prédire, analyser et surveiller le
comportement dynamique du système, un modèle mathématique fiable doit être établi. Le diagnostic
peut ainsi faire appel à divers niveaux de modélisation permettant de comparer le comportement
observé et attendu afin de signaler toute incohérence liée probablement à l’apparition d’un ou
plusieurs défauts.
12 Chapitre 1. Introduction au diagnostic

1.3.1 Modélisation du comportement nominal d’un système


Dans le cadre de ce cours, nous considérons les systèmes linéaires sous leurs formes les plus
employées en automatique à savoir : la fonction de transfert et la représentation d’état.
La fonction de transfert
On appelle la fonction de transfert d’un système ou d’un procédé, le rapport de la transformée
de Laplace du signal de sortie y(t) à celui de l’entrée u(t) i.e :
Y (s)
Gyu (s) =
U(s)
où Y (s) et U(s) représentent les transformées de Laplace des signaux u(t) et y(t) respectivement.
En temps discret, la fonction de transfert est obtenue à partir de la transformée en Z des signaux
d’entrée et de sortie comme suit :
Y (z)
Gyu (z) = (1.1)
U(z)
Si les signaux y ∈ ℜny et u ∈ ℜnu sont multi-dimensionnels, alors Gyu est une matrice de ny lignes
et nu colonnes et dont chaque élément est une fonction de transfert.

F IGURE 1.7 – Vue externe d’un système

La représenation d’état
La forme standard de la représentation d’un système linéaire à temps continu dans l’espace
d’état est donnée par
ẋ(t) = Ax(t) + Bu(t) (1.2)
y(t) = Cx(t) + Du(t)
où x ∈ ℜnx ,u ∈ ℜnu et y ∈ ℜny représentent l’état, l’entrée et la sortie du système, respectivement.
Les matrices A, B, C et D sont des matrices réelles avec des dimensions appropriées. Cependant,
pour le système linéaire à temps discret nous utilisons la forme suivante
xk+1 = Axk + Buk (1.3)
yk = Cxk + Duk

R Une système linéaire possède une seule fonction de transfert et une infinité de représentations
d’état possibles.

Le passage de l’espace d’état vers la fonction de transfert est effectué comme suit :
Gyu (s) = C(sInx − A)−1 B + D (en temps continu)
Gyu (z) = C(zInx − A)−1 B + D (en temps discret)
Où Inx est une matrice identité de taille nx . Nous pouvons aussi écrire le système sous la forme
réduite suivante :
 
A B
Gyu (s) = (1.4)
C D
1.3 Modélisation des défauts et des systèmes dynamiques 13

1.3.2 Les systèmes perturbés


Les systèmes physiques subissent généralement l’effet néfaste des signaux exogènes à savoir :
les perturbations et les bruits de mesures. Ces signaux sont imprévisibles et nécessitent le recours
à des techniques robustes afin de diminuer leurs impact sur les performances de fonctionnement
requises. Selon une étude préalable, il est possible d’adopter l’un des trois modèles de signaux
exogènes :
— un processus stochastique ;
— une entrée inconnue ;
— une perturbation bornée ;
Le système stochastique
Un signal stochastique ou aléatoire est une famille de variables aléatoires indexée par le temps k.
De même, un système stochastique linéaire est un système linéaire perturbé par un signal aléatoire
et possède la représentation d’état suivante :
xk+1 = Axk + Buk + wk (1.5)
yk = Cxk + Duk + vk
où wk et vk représentent deux signaux aléatoires d’une distribution connue. Si wk et vk sont deux
signaux blancs, Gaussiens, non-corrélés et de covariances V  0 et W  0, alors le système (1.5) est
appelé un système de Gauss-Markov. Nous pouvons ainsi écrire : wk ∼ N (0,W ) et vk ∼ N (0,V )
pour indiquer la forme Gaussienne (Normalienne) de la distribution et E(vi v j ) = V δ (i − j) (δ
représente la fonction de Kronecker) pour montrer la blancheur du signal.
Le système à entrée inconnue
Les perturbations autour du processus considéré, les changements inattendus dans le processus
technique ainsi que les bruits de mesure sont souvent modélisés en tant qu’un vecteur à entrée
inconnue. Ce vecteur, noté dk , est intégré dans le modèle d’état (1.2) comme suit :
xk+1 = Axk + Buk + Ed dk (1.6)
yk = Cxk + Duk + Fd dk
Les matrices de distribution de l’entrée inconnue Ed et Fd indiquent la direction vers laquelle la
perturbation dk influe sur le système. L’équation ci-dessous donne une représentation du système
(1.6) dans le domaine fréquentiel.
Y (z) = Gyu (z)U(z) + Gyd d(z) (1.7)
avec
 
A B
Gyu (z) = = C(zInx − A)−1 B +C (1.8)
C D
et
 
A Ef
Gyd (z) = = C(zInx − A)−1 E f + Ff (1.9)
C Ff
La perturbation bornée
Il est possible d’établir une borne supérieure pour l’énergie de la perturbation dk dans le système
(1.6) comme suit :
! 12

k dk k2 = ∑ di> di < l où l ∈ ℜ† (1.10)
i=0
14 Chapitre 1. Introduction au diagnostic

La norme k · k2 , appelée la norme L2 , est liée directement à l’énergie du signal. Cependant, le


produit di> di correspond à la puissance instantanée du signal. La borne supérieure l peut être
obtenue empiriquement.

1.3.3 Modélisation des défauts


Il existe plusieurs façons de modélisation des défauts. L’extension du modèle (1.7) à

Y (z) = Gyu (z)U(z) + Gyd (z)d(z) + Gy f (z) f (z) (1.11)

est une forme largement adoptée au diagnostic, où f ∈ ℜn f est un vecteur inconnu représentant
tout les défauts possibles (nul dans le cas non-défaillant). Nous acceptons que Gy f (z) est matrice
de fonctions de transfert connues et f est une fonction déterministe. Une réalisation minimale de
(1.11) est donnée par

xk+1 = Axk + Buk + Ed dx + E f fk (1.12)


yk = Cxk + Duk + Fd dk + Ff fk

Les matrice E f ∈ ℜnx ×n f et Ff ∈ ℜny ×n f sont appelées les matrices de distribution de défaut. Il est
évident que E f , Ff indiquent l’endroit où un défaut se produit et son influence sur la dynamique du
système.
On a alors

Gy f (z) = C(zI − A)−1 E f + Ff (1.13)

Comme illustré à la F IGURE1.2.1, les défauts peuvent être classés en trois catégories :
— défaut capteur fs : ce sont les défauts qui agissent directement sur la mesure.
— défaut actionneur fa : Ces défauts provoquent des changements dans le fonctionnement des
actionneurs.
— défaut procédé f p : Ils sont utilisés pour indiquer des dysfonctionnements dans le système.
Si on a un défaut capteur, nous pouvons mettre Ff = I ce qui donne

yk = Cxk + Duk + Fd dk + fs (1.14)

Cependant, le défaut actionneur est modélisé en mettant E f = B et Ff = D, donc

xk+1 = Axk + B(uk + fA ) + Ed dk


yk = Cxk + D(uk + fA ) + Fd dk (1.15)

En fonction de leurs types et emplacements, les défauts de procédé peuvent prendre des différentes
des matrices du système i.e : E f = E p et Ff = Fp . Si un système subissant l’effet des défauts capteur,
actionneur et procédé alors
 
fA    
f =  fP  , E f = B E p 0 , Ff = D Fp I
fS

En raison de la manière dont ils affectent la dynamique du système, les défauts décrits par (1.19)
sont appelés défauts additifs. Il est très important de noter que l’apparition d’un défaut additif
n’affectera jamais la stabilité du système. Les défauts additifs typiques rencontrés dans la pratique
sont, par exemple, un offset dans les capteurs et les actionneurs ou une dérive dans les capteurs. Le
premier peut être décrit par une constante, tandis que le second par une rampe.
1.4 Classification des méthodes de diagnostic 15

En pratique, les défauts dans le procédé, les capteurs et actionneurs entraînent souvent des
modifications dans les paramètres du modèle. Ils sont appelés défauts multiplicatifs et modélisés
généralement en termes de changements de paramètres. On peut les décrire en étendant (1.6) au
modèle suivant

xk+1 = (A + ∆AF )xk + (B + ∆BF )uk + Ed dk (1.16)


yk = (C + ∆CF )xk + (D + ∆DF )uk + Fd dk (1.17)

où les matrices ∆AF , ∆BF , ∆CF , ∆DF représentent les défauts multiplicatifs dans le procédé, les
actionneurs et les capteurs, respectivement. Les défauts multiplicatifs sont caractérisés par leur
impact (possible) sur la stabilité du système.
Exercice 1.1 Considérons les défauts multiplicatifs suivants :

lA lB lC lD
∆AF = ∑ Ai θAi , ∆BF = ∑ Bi θBi , ∆CF = ∑ Ci θCi , ∆DF = ∑ Di θDi (1.18)
i=1 i=1 i=1 i=1

où Ai , i = 1, . . . , lA , Bi , i = 1, . . . , lB , Ci , i = 1, . . . , lC ,Di , i = 1, . . . , lD , sont des matrices connues


avec des dimensions appropriées.
1- Réécrire le modèle (1.16) sous la forme suivante :

xk+1 = Axk + Buk + Ed dx + Ē f fM (1.19)


yk = Cxk + Duk + Fd dk + F̄f fM

avec fM est une fonction inconnue qui dépend de xk et uk . 

1.4 Classification des méthodes de diagnostic


Dans les dernières décennies, plusieurs méthodes de diagnostic ont été développées par diffé-
rentes communautés de recherche et comme illustré à la Figure 1.8, ces méthodes sont généralement
classées en deux grands ensembles :
1. Méthodes à base de données (signal).
2. Méthodes à base de modèle.

1.4.1 Méthodes à base de données (signal)


L’idée centrale de ces méthode est d’extraire l’information de défaut à partir des signaux du
processus surveillé (consignes, mesures). Pour atteindre cet objectif, certaines propriétés du signal
(symptômes) peuvent être analysées dans le domaine temporel ou dans le domaine fréquentiel. Les
caractéristiques temporel des symptômes comprennent : l’amplitude, la moyenne (arithmétique
ou quadratique), les valeurs limites, les tendances, et les moments statistiques de la distribution
d’amplitude, etc.., tandis que les caractéristiques du domaine fréquentiel comprennent la densité
spectrale de puissance, les lignes spectrales de fréquence etc. Le diagnostic à base de signal est
utilisé pour surveiller les systèmes fonctionnant en régime permanent.

1.4.2 Méthodes à base de modèle


L’idée intuitive des méthodes à base de modèle consiste à remplacer la redondance matérielle par
un modèle de processus implémenté dans un calculateur. Le modèle implémenté est sollicité en ligne
par les mêmes entrées de commande que le processus réel. De cette façon, le comportement sain de
processus peut être reconstitué en ligne. Par analogie à la redondance matérielle, nous appelons
cette méthode “redondance logicielle” ou “redondance analytique” . Il est bien connu que les
16 Chapitre 1. Introduction au diagnostic

Méthodes de détection
et de diagnostic des
défauts

Méthodes à base de Modèle Méthodes à base de données

Méthodes quantitatives Méthodes qualitatives


Méthodes qualitatives
Méthodes quantitatives
Estimation Estimation Espace
Modèle causal
d’état Estimation État/paramètre de parité Systèmes
paramétrique Hiérarchie experts
d’abstraction
Graphe Réseau de Logique
Filtre de Espace
structurel Physique Neurones
MC/MCR d’état floue
Kalman qualitative
observateur étendu statistique
Analyse Filtre de Entrée/sortie Arbre de structurelle Analyse Des
de régression Kalman à défaillance fonctionnelle tendances
Filtre de Kalman deux étages qualitatives
ACP/MCP Classification Analyse de
Automatique
Fréquence
temps-fréquence

F IGURE 1.8 – Les méthodes de diagnostic de défaut

techniques à base de modèle sont plus puissantes que celles à base de signal grâce aux informations
supplémentaires fournies par la connaissance du comportement dynamique de système. Le modèle
du processus représentant le comportement qualitatif et quantitatif du procédé est obtenu en utilisant
la théorie bien établie d’identification et de modélisation des systèmes. Le modèle quantitatif
ou analytique du processus peut être représenté par un ensemble d’équations différentielles ou
des équations récurrentes alors que le modèle qualitatif est décrit par des fonctions qualitatives
centrées autour de différentes unités du processus. Comme illustré à la F IGURE 1.8, les méthodes
quantitatives à base de modèle sont généralement classées en trois ensembles : (i) méthode d’espace
de parité (ii) méthodes à base d’observateur (iii) méthodes à base d’estimation (identification)
paramétrique (iv) à base d’estimation d’état/paramètres. Dans ce cours nous nous consacrons sur
l’étude des méthodes à base de modèle, leurs fondements théoriques et les applications possibles
dans des divers domaines techniques.

Étapes du diagnostic
Comme illustré à la F IGURE 1.9, le diagnostic à base de modèle passe généralement par les
trois étapes suivantes :
1- Détection de défaut décider si un défaut est survenu. Cette étape détermine le temps au quel le
système est soumis à un défaut. Cette étape se compose de deux sous-étapes suivantes :
— Génération du résidu : le résidu est un signal qui reflète le défaut. La génération de résidus
se fait généralement en temps réel à partir des entrées et des sorties dy système.
— Évaluation du résidu : le résidu est analysé pour décider s’il y a ou non présence d’un
défaut. Cette prise de décision peut être réalisée en temps réel à l’aide d’un simple test de
dépassement de seuil sur les valeurs instantanées ou en analysant les propriétés statistiques
du signal résidu, ou bien encore par la décision floue. On génère aussi un symptôme.
1.4 Classification des méthodes de diagnostic 17

F IGURE 1.9 – Les étapes principales d’un algorithme de diagnostic

2- Localisation de défaut : trouver dans quelle composante du système le défaut s’est produit.
Cette étape détermine l’emplacement de la panne.
3- Identification de défaut : identifier le défaut et estimer son amplitude. Cet étape détermine le
type de défaut et sa gravité.
Notons que l’adoption d’une étape de plus parmi les étapes précédentes, permet d’accroitre
le niveau de sécurité de l’installation, et par conséquent, minimiser davantage le coût de dégâts
probables (voir F IGURE 1.10)

F IGURE 1.10 – Niveau de sécurité et étapes de diagnostic

1.4.3 Performances d’un algorithme de diagnostic :


Pour déterminer les performances d’un algorithme de diagnostic nous utilisons les critères
suivant :
— La sensibilité aux défauts : l’aptitude de l’algorithme de capter le défaut et d’éviter la
non-détection (missed alarm).
— La vitesse de détection : mesurée par le temps pris par l’algorithme pour signaler l’existence
d’un défaut après son apparition.
— La robustesse : l’aptitude de minimiser l’effet des signaux exogènes et des erreurs de
modélisation sur la décision du diagnostic.

R En implémentant un algorithme du diagnostic en temps réel, sa complexité temporelle (nombre


18 Chapitre 1. Introduction au diagnostic

des opérations de calcul) et spatiale (espace mémoire utilisé) s’avère d’une importance
primordiale dans l’évaluation de ces performances.

1.5 Exercices d’application

F IGURE 1.11 – Schéma d’un MCC

Exercice 1.2 La F IGURE 1.2 ci-dessus montre un schéma descriptif d’un moteur à courant continu.
Afin d’obtenir un modèle mathématique décrivant la dynamique du système, nous avons choisi
le courant IA et la vitesse du moteur Ω comme les éléments du vecteur d’état (variables d’état).
En appliquant les lois fondamentales de la mécanique et de l’électricité, nous aboutissons à la
représentation d’état suivante :

I˙A − RLAA − CΦ 1
        
LA IA LA 0
= + UA + ML
Ω̇ KM
J 0 Ω 0 − 1J

où J = 80.45 · 10−6 kg · m2 , CΦ = 6.27 · 10−3V /Rpm, KM = 0.06Nm/A, LA = 0.003H, RA =


3.13Ohm, KT = 2.5 · 10−3V /Rpm , TT = 5ms, ML représentent : l’inertie, constante de la ten-
sion, constante du moteur, inductance des armateurs, résistance, le gain, le constant du temps de
tachymètre et la charge entrainée par le moteur (inconnue), respectivement.
1. Quelles sont les variables d’entrée et de sortie.
2. Que représente charge Ml :
a) une entrée de commande ; b) une perturbation ; c) un bruit de mesure ; Choisissez la bonne
réponse.
3. Déduire la fonction de transfert du système ;
Comme illustré à la F IGURE 1.5 ci-dessous et afin de commander le moteur, deux boucles
de commande en cascade sont adoptées : la boucle de commande du courant et la boucle de
commande de vitesse.
Considérant le moteur + le contrôleur du courant IA − ctrl comme un seul système (processus)
commandé par le contrôleur de la vitesse de type PI Ω − ctrl.
4. Que représentent les variables : Ωref , Ω, Iref .
La dynamique du processus (IA ctrl+ moteur) est donnée par

y(s) = GyuU(s) + Gyd (s)d(s)

8.75 31.07
avec Gyu = (1+1.225s)(1+0.03s)(1+0.005s) et Gyd (s) = − s(1+0.005s)
1.5 Exercices d’application 19

D’autre part le contrôleur de la vitesse est donné par


1 + 1.225s
K(s) = 1.6
s
5. Quel est le type de ce contrôleur ? Donner l’expression de U(s).
6. Le modèle du système en boucle fermée est donné par

Y (s) = Gyw (s)W (s) + GBF


yd (s)d(s) (1.20)

7. Trouver Gyw (s) et GBF


yd (s). Dans la suite, nous considérons trois types de défauts :
— un défaut additif dans l’actionneur fA .
— un défaut additif dans le Tachymètre fs1 .
— un défaut multiplicatif dans le Tachymètre fs2 ∈ [0, −1].
Le modèle du système défaillant en boucle ouverte possède la forme suivante :

y(s) = Gyu (s)U(s) + Gyd (s)d(s) + Gy fA (s) fA (s) + Gy fS1 (s) fS1 (s) + ∆y(s) fS2 (s)

8. Trouver Gy fA (s), Gy fS1 (s), ∆y(s). D’autre part, nous pouvons écrire le modèle du système
défaillant en boucle fermée comme suit :

y(s) = Gyw (s)W (s) + GBF BF BF BF


yd (s)d(s) + Gy fA (s) f A (s) + Gy fS1 (s) f S1 (s) + ∆y (s) f S2 (s)

9. Trouver GBF BF BF BF
yd , Gy fA , Gy fS1 , ∆y .

F IGURE 1.12 – Commande en boucle fermée

Exercice 1.3 — Simulation MATLAB/SIMULINK. Considérons le système linéaire défini par la


représentation d’état discrète suivante :
   
(
xk+1 = Axk + Buk + F fk 0.5 1 0 0 1  
1 0 1
, A =  0 0.1 0  , B = 1 1  , C=
yk = Cxk 0 1 1
0 0 0.5 0 −1

1. A l’aide des blocs Simulink élémentaires (addition, intégration, multiplication), simuler


le comportement de ce système sollicité par une entrée de type échelon et subissant un
défaut égale à 0.6 à l’instant k = 25.
2. Refaite la même simulation en utilisant le bloc prédéfini "state-space".

2. Espace de parité

2.1 Introduction
L’idée de base de l’approche par espace de parité est de vérifier la cohérence entre les relations
mathématiques du système et les mesures (relations de redondance analytique). Supposons en effet,
qu’une mesure puisse s’exprimer en fonction des autres par une relation connue. La différence
entre la mesure et sa valeur calculée à l’aide du modèle est appelée résidu. Si le résidu est nul, les
mesures sont cohérentes par rapport au modèle et le système est déclaré sans défaut. Un résidu
non nul indique l’apparition d’un défaut. L’approche par espace de parité repose donc sur une
connaissance préalable du modèle mathématique du système. Il existe deux types de relations de
redondance analytique :
— La redondance statique : ensemble de relations algébriques entre les mesures fournies par les
différents capteurs.
— La redondance dynamique : ensemble d’équations différentielles ou récurrentes entre les
sorties des capteurs et les entrées du système.
Avant de développer ces relations, un petit rappel sur le calcul des noyaux associés à une matrice est
donné. Ceci permet de calculer la matrice de parité qui est d’un intérêt majeur dans l’établissement
du résidu.

2.2 Sous-epaces vectoriels nuls


Définition 2.2.1 Soit A ∈ ℜnl ×nc une matrice réelle quelconque où nc est le nombre des colonnes
et nl est le nombre de lignes. Le sous-espace vectoriel nul associé à A, aussi appelé le noyau de
A ou ker(A), est l’ensemble des vecteurs x solution du système d’équations homogène suivant
Ax = 0 ç-à-d :

ker(A) = {x ∈ ℜnc : Ax = 0} (2.1)

où 0 est un vecteur nul (composé de zéros).


 >  >
 Exemple 2.1 Les vecteurs x1 = −2 1 0 et x2 = −3 0 1 appartiennent au noyau
22 Chapitre 2. Espace de parité

F IGURE 2.1 – Sous epace vectoriel nul d’une matrice


 
1 2 3
de la matrice A =  1 2 3  

1 2 3

R La dimension de ker(A) est donnée par la relation suivante :


dim(ker(A)) = nc − rang(A) (2.2)

La dimension de l’espace vectoriel ker(A) correspond au nombre minimum de vecteurs qui com-
posent une base génératrice de cet
 espace.
 En effet, pour trouver la base de ker(A), il est possible de
A
former une matrice augmentée , puis on procède à échelonnement par l’élimination Gauss
Inc
 
M
pour aboutir finalement à la forme suivante . Une base du noyau de A peut être formée à
N
partir des colonnes non nulles de N qui se trouve directement au dessous des colonnes nulles de M.
 Exemple 2.2 Supposons que
 
1 0 3 0 2 −8
 0 1 5 0 −1 4 
A= 
 0 0 0 1 7 −9 
0 0 0 0 0 0
On a rang(A) = 3 ⇒ dim(ker(A)) = nc − rang(A) = 6 − 3 = 3 alors la base du noyau est formé de
trois vecteurs.
En augmentant la matrice A,
 
1 0 −3 0 2 −8
 0 1 5 0 −1 4 
 
 0 0 0 1 7 −9 
 
 0 0 0 0 0 0 
    
A  1 0 0 0 0 0 

= .
I6  0 1 0 0 0 0 
 0 0 1 0 0 0 
 
 
 0 0 0 1 0 0 
 
 0 0 0 0 1 0 
0 0 0 0 0 1
2.3 Redondance statique 23

Après l’échelonnement de la matrice augmentée, nous obtenons :


 
1 0 0 0 0 0
 0 1 0 0 0 0 
 
 0 0 1 0 0 0 
 
 0 0 0 0 0 0 
   
M  1 0 0 3 −2 8 
 
=
N  0 1 0 −5 1 −4 

 0 0 0 1 0 0 
 
 0 0 1 0 −7 9 
 
 
 0 0 0 0 1 0 
0 0 0 0 0 1
Les trois dernières colonnes de M sont nulles. Par conséquent, les trois derniers vecteurs de N,
−2
     
3 8
 −5   1   −4 
     
 1   0   0 
 ,  ,  
 0   −7   9 
     
 0   1   0 
0 0 1
Forment une base pour le noyau de A. 

Définition 2.2.2 Soit A ∈ ℜnl ×nc une matrice réelle quelconque. Le sous-espace vectoriel nul
à gauche associé à A, aussi appelé le conoyau de A ou co-ker(A), est l’ensemble des vecteurs
lignes x solution du système d’équations homogène suivant xA = 0 ç-à-d :

co-ker(A) = {x ∈ ℜ1×nl : xA = 0} (2.3)

où 0 est un vecteur ligne nul (composé de zéros).

R La dimension de co-ker(A) est donnée par la relation suivante :


dim(co-ker(A)) = nl − rang(A) (2.4)
Cette relation est utilisée pour le calcul de nombre de lignes de la matrice de parité, formée
principalement à partir d’une base génératrice du co- ker(A). Il est possible de trouver les base
de conoyau en appliquant la méthode d’élimination de Gauss pour la matrice A> . Ensuite, les
vecteurs colonnes trouvés doivent être transposés pour avoir les vecteurs lignes formant la
base du conoyau.

Exercice 2.1 Trouver une base génératrice du conoyau de A.


 
1 0 3 0 2 −8
 0 1 5 0 −1 4 
A= 
 0 0 0 1 7 −9 
0 0 0 0 0 0


2.3 Redondance statique


Dans un système physique, les variables mesurées sont généralement liées par un ensemble
de relations algébriques. L’objet de la redondance statique est de trouver les relations existantes
24 Chapitre 2. Espace de parité

Générateur de
Capteurs le résidu
Système Résidu
physique C V k

F IGURE 2.2 – La redondance statique

entre les mesures fournies par les différentes capteurs. Ceci est réalisée à l’aide d’un modèle
mathématique du système de mesure qui s’écrit, généralement, de la façon suivante :

yk = Cxk + fk (2.5)

où yk ∈ ℜny est le vecteur des mesures, C ∈ ℜny ×nx est la matrice d’observation, xk ∈ ℜnx est le
vecteur d’état et fk ∈ ℜny est le vecteur des défauts de capteurs. Les équations de redondance
sont obtenues par élimination du vecteur d’état xk dans (2.5). Donc, il est possible de trouver une
matrice C⊥ , dite de parité, orthogonale à C (c-à-d C⊥C = 0) et permettant de trouver des relations
indépendantes liant les mesures entre eux comme illustré à la F IGURE 2.2. La matrice C⊥ est
formée à partir d’une base génératrice de l’espace vectoriel nul à gauche co- ker(C). Ceci n’est
possible sauf si la dimension de co- ker(C) est non-nulle. A partir de (2.4), nous écrivons

dim(co- ker(C)) = ny − rang(C) > 0 ⇒ ny > rang(C) (2.6)

Sous cette dernière condition et en multipliant les deux côtés de l’équation (2.5) par la matrice
V ⊥ , nous obtenons :

rk = C ⊥ yk (2.7)

où rk est le vecteur résidu de dimension ny − rang(C). Cette relation est utilisée pour former les
équations de redondance statiques (résidus) en fonction des mesures en provenance des capteurs. Ces
relations sont implémentées dans un calculateur (générateur de résidu) afin d’assurer la détection
du défaut. Il est à noter que le résidu donne une image claire sur le vecteur défaut car implicitement
on a la relation suivante :

ny
rk = C⊥Cxk +C⊥ fk = C⊥ fk = ∑ Ci⊥ fi (2.8)
| {z } i=1
0

où Ci⊥ le ième vecteur colonne de la matrice C⊥ . Ceci permet de détecter l’apparition du défaut,
c-à-d de décider est ce que f = 0 ou f 6= 0. En cas d’absence de défaut, le vecteur résidu rk est
nul. Dans le cas de l’apparition d’un défaut fi (t), le vecteur résidu s’oriente vers la direction de Ci⊥
correspondant au vecteur défectueux.
 Exemple 2.3 soit l’équation de mesure suivante :
     
y1 1 f1
 y2  =  1  x +  f2 
y3 1 f3

La matrice C⊥ est choisie (par l’élimination de Gauss) comme suit


 
⊥ −1 1 0
V = (2.9)
−1 0 1
2.3 Redondance statique 25

le résidu obtenu est sous la forme


 
−y1 + y2
rk = (2.10)
−y1 + y3
sachant que (implicitement) on a la forme équivalente suivante :
     
−1 1 0
rk = f1 + f2 + f
−1 0 1 3
C1⊥ C2⊥ C3⊥

L’interprétation géométrique de ces relations est appelée résidu directionnel. Ce type de résidu

C3⊥

C2⊥
C1⊥
𝑟𝑘

F IGURE 2.3 – Le résidu directionel

donne une autre possibilité pour réaliser la tâche de localisation de défauts. Il est construit tel que,
en réponse à un défaut donné, le vecteur de résidus soit orienté suivant une direction bien précise
de l’espace des résidus. Dans notre exemple (F IGURE 2.3) le résidu est orienté vers le vecteur C2⊥ ,
ce qui indique l’apparition du défaut f2 . 

Supposons que la matrice C est de rang plein colonne i.e rang(C) = nx . Les lignes de C⊥
forment une base de l’espace de parité de dimension ny − nx . Une façon simple de déterminer une
matrice de parité C⊥ est de réarranger l’équation (2.5). En absence de défaut ( fk = 0), la relation
(2.5) peut être écrite sous la forme suivante :
   
ynx C1
ỹk = C̃xk , avec ỹk = , C̃ = (2.11)
yny −nx C2
La matrice C1 de la relation (2.11) est une matrice carrée composée de nx lignes indépendantes
de C,pour assurer l’existence de son inverse. La matrice C2 est obtenue à l’aide de ny − nx lignes
restantes de la matrice C. Nous obtenons alors :
yny −nx = C2 C1−1 ynx (2.12)
| {z }
xk

Cette relation est indépendante du vecteur d’état et permet de vérifier la cohérence des mesures.
La relation (2.12) peut encore s’écrire sous la forme suivante :
rk = −C2C1−1 Iny −nx ỹk = 0
 
(2.13)
La matrice de parité C̃⊥ = −C2C1−1 Iny −nx de dimension (ny − nx × nx ) est orthogonale à
 

C̃, ç-à-d (C̃⊥C̃ = 0).


26 Chapitre 2. Espace de parité

 Exemple 2.4 L’équation de mesure d’un système linéaire est donné par :
 
1 2
 1 0 
y=
 1
x
1 
2 0

Nous définissons les matrices C1 et C2 comme suit :


   
1 2 1 1
C1 = et C2 =
1 0 2 0

La matrice C ne subit pas un changement d’ordre des lignes, on a ainsi :

C̃ = C et ỹ = y

Les équations de redondance analytique sont données par


 
1 1 −2 0
y=0
0 4 0 −2

Donc
(
r1 = y1 + y2 − 2y3
r2 = 4y2 − 2y4


R La matrice C⊥ peut aussi être obtenue en sélectionnant ny − nx lignes linéairement indé-


pendantes de la matrice M = Iny −CC† . On a C† représente la matrice pseudo-inverse de la
matrice C.

Exercice 2.2 Soit le système linéaire suivant :


 
1 0 1

 1 2 1 

yk = Cxk + fk où C=
 2 0 2 ;

 1 0 1 
1 1 1

1. Calculer le rang de la matrice C.


2. Combien de relations de redondance statique existent elles ?
3. Construire la matrice de parité C⊥ .
4. Trouver les équations de redondance statique.


R La méthode de redondance statique possède deux inconvénients majeurs :


— Le grand nombre des capteurs utilisés ce qui augmente le coût de l’instrumentation du
système.
— La spécialisation dans la détection des défauts de capteurs limite son intérêt pratique.
Cependant, la méthode de redondance dynamique permet de surmonter ces difficultés en
employant une augmentation temporelle du modèle (lifting).
2.4 Redondance dynamique 27

2.4 Redondance dynamique


La redondance dynamique est une généralisation de la redondance statique dans le cas d’uti-
lisation d’un modèle dynamique du système étudié. L’objectif est de trouver des relations entre
les mesures fournies par les différents capteurs et les entrées du système à différents instants.
Considérons à cet effet un système supposé correctement représenté par le modèle d’état discret
suivant : On suppose que la dynamique du système illustré à la Figure 2.4 est décrite par
(
xk+1 = Axk + Buk + E fk
(2.14)
yk = Cxk + Duk + F fk

où xk ∈ ℜnx ,uk ∈ ℜnu ,yk ∈ ℜny , fk ∈ ℜn f sont les vecteurs : d’état, de commandes, de mesures, de
défauts, respectivement. Soit le vecteur ?s,k donné par

T
?T (k − s) ?T (k − s + 1) · · · ?T (k)

?s,k = (2.15)

où ? peut représenter u, y, f et s est la taille de la fenêtre d’observation.

F IGURE 2.4 – Génération de résidu via l’espace de parité

Proposition 2.4.1 L’évolution du système dans un horizon temporel de taille s peut être résumée
comme suit :

ys,k = Ho,s x(k − s) + Hu,s us,k + H f ,s fs,k (2.16)


 
D 0 ··· 0
 CB D ··· 0 
Hu,s =  (2.17)
 
.. .. .. .. 
 . . . . 
CAs−1 B CAs−2 B . . . D

La matrice H f ,s est obtenue en remplaçant la matrice B et D dans (2.17) par E et F, respectivement.


La matrice Ho,s , est appelée la matrice d’observabilité réduite et donnée par :
>
C> A>C> · · · (As )>C>

Ho,s = (2.18)

Démonstration. D’après (2.14), nous écrivons


28 Chapitre 2. Espace de parité

• à l’instant k − s nous obtenons


yk−s = Cxk−s + Duk−s (2.19)
• à l’instant k − s + 1, nous avons yk−s+1 = Cxk−s+1 + Duk−s+1 et xk−s+1 = Axk−s + Buk−s d’où

yk−s+1 = CAxk−s +CBuk−s + Duk−s+1 (2.20)


• à l’instant k − s + 2, nous avons yk−s+2 = Cxk−s+2 + Duk−s+2 et xk−s+2 = Axk−s+1 + Buk−s+1
d’où
yk−s+2 = CA2 xk−s +CABuk−s +CBuk−s+1 + Duk−s+2 (2.21)
..
.
• à l’instant k, nous avons yk = Cxk + Duk et xk = Axk−1 + Buk−1 d’où
yk−s+2 = CAs xk−s +CAs−1 Buk−s +CAs−2 Buk−s+1 + · · · +CBuk−1 + Duk (2.22)
En écrivant les relations (2.19)-(2.22) sous forme matricielle, nous obtenons la forme augmentée
(2.16).

Le résidu rs,k est donné par
rs,k = Vs⊥ (ys,k − Hu,s us,k ) (2.23)
La dynamique du résidu est décrite par l’équation suivante
rs,k = Vs⊥ H f ,s fs,k (2.24)
où Vs⊥ est la matrice de parité déterminée en résolvant les équations suivantes :
Vs⊥ Ho,s = 0 (2.25)

R Pour vérifier la condition de détectabilité, la matrice de parité doit en plus satisfaire la


condition suivante :
Vs⊥ H f ,s 6= 0 (2.26)

La taille de la fenêtre d’observation s est liée à l’indice d’observabilité du système ν expliqué dans
la définition suivante.
Définition 2.4.1 L’indice d’observabilité ν d’un système linéaire à temps discret (ou de la paire
(A,C)) est le plus petit entier naturel vérifiant la condition

rang(Ho,ν ) = rang(Ho,ν+1 ) (2.27)

Pour trouver les relations de redondance dynamique, la taille de la fenêtre d’observation doit vérifier
l’inégalité suivante :

ν
≤ s ≤ ν − rang(C) + 1 (2.28)
rang(C)
Si le système est observable et les lignes de la matrice C sont linéairement indépendantes, alors
l’inégalité (2.28) prend la forme :
nx
≤ s ≤ nx − ny + 1 (2.29)
ny
2.4 Redondance dynamique 29

Exercice 2.3 Démontrer la condition (2.28). 

2.4.1 Auto-redondance
Les relations d’auto-redondance sont obtenues en écrivant la relation (2.23) pour chaque
capteur. On ne conserve alors, pour un capteur donnée, que les relations indépendantes permettant
d’exprimer une partie de l’état. Par exemple, pour le capteur numéro j la relation, (2.23) s’écrit :

( j) ( j) ( j)
ys j ,k = Ho,s j x(k − s j ) + Hu,s j us j ,k (2.30)


 
D 0 ··· 0
( j) > ( j)
 C jB D ··· 0 
C> A>C> (As j )>C>

Ho,s j = ··· et Hu,s j = 
 
j j j .. .. .. .. 
 . . . . 
C j As j −1 B C j As j −2 B . . . D

où C j représente le vecteur ligne numéro j de la matrice C. Soit ν j l’indice d’observabilité de


la paire (A,C j ), alors la fenêtre temporelle pour chaque capteur j comporte ν j + 1 lignes ç-à-d
( j)
s j = ν j + 1. La matrice du rang maximum Ho,ν j est appelée matrice d’observabilité réduite.
L’équation d’auto-redondance est alors obtenue en éliminant l’état de la relation (2.30). Pour
cela, on cherche un vecteur ligne V ⊥( j) tel que :

( j)
V ⊥( j) Ho,s j = 0 (2.31)

L’équation d’auto-redondance s’écrit alors :


 
( j) ( j) ( j)
rk = V ⊥( j) ys j ,k − Hu,s j us j ,k (2.32)

( j)
En absence de défaut, la grandeur rk est réduite aux bruits des mesures et des structures. L’appari-
( j)
tion d’un défaut se traduit par une évolution de rk ; ce qui permettra la détection de l’anomalie.

2.4.2 Inter-redondance
Les relations d’inter-redondances permettent de relier les mesures provenant de plusieurs
capteurs. On les obtient en considérant les ν j ( j = 1 à q) relations indépendantes obtenues à partir
de (2.30) :
  
(1) (1)
  (1)  
yν1 ,k+ν1 Ho,ν 1
H u,ν 1 u ν ,k+ν

 .   1 1
 .   ..   .  . 
 ..   .. 
 .   . 
 ( j)   ( j)   
( j)  

y  = Ho,ν  xk + 

Hu,ν j   uν j ,k+ν j 
 (2.33)
 ν j ,k+ν j   j

 .  . 
 .   .   .  . 
 ..   ..    .  .
 
(q) (q) (q) u
yνq ,k+νq Ho,νq Hu,νq ν q ,k+νq

Le système (2.33) est composé de N relations indépendantes, avec N = ∑qj=1 ν j . Le vecteur d’état
étant de dimension nx , il existe N − nx relations d’inter-redondance indépendantes. Les équations
d’inter-redondance sont obtenues par élimination du vecteur d’état. Cela revient à trouver une
30 Chapitre 2. Espace de parité

matrice V ⊥ telle que :


 (1) 
Ho,ν
 . 1
 .. 
 
 ( j) 
V ⊥ Ho,ν =0 (2.34)
 . j
 . 
 . 
(q)
Ho,νq

Les relations d’inter-redondance s’écrivent alors :


   
(1) (1) 

yν1 ,k+ν1 
 .  Hu,ν1  uν1 ,k+ν1 
 .   ..   .. 
 .   .   . 

 ( j)   ( j)   
rk+ν ? = V  yν j ,k+ν j  −  Hu,ν j   uν j ,k+ν j  (2.35)
  
 .   .

 .  
 ..   .   .  
   .  . 

(q) (q) uνq ,k+νq
yνq ,k+νq Hu,νq

2.5 Exercices d’application


Exercice 2.4 Les relations de redondance statique
Soit le système de mesure suivant :
     
y1 (t) 1 8 1   f 1 (t)
y2 (t) 1 0 8 x1 (t)  f2 (t)
     
y3 (t) = 1 1 1 x2 (t) +  f3 (t)
     
y4 (t) 1 0 1 x3 (t)  f4 (t)
y5 (t) 8 0 8 f5 (t)

1. A partir de la matrice d’observation C, construire deux matrices C1 et C2 tels que : C1 carrée


et inversible, C2 représente les lignes restantes de C, nx le nombre des variables d’état et ny le
nombre des mesures.  
>  
C1 ynx
2. Déduire la matrice C̃ = et le vecteur mesure modifié ỹ = .
C2 yny −nx
3. En fonctionnement nominal, vérifier la relation : −C2C1−1 ynx + yny −nx = 0.
4. Déterminer la matrice de parité C⊥ tel que : C⊥C= 0.
5. Donner les équations de redondance statiques. Que ce passe t-il si le capteur de mesure y2 est
défectueux ?
Exercice 2.5 Les relations de redondance dynamique
Considérons le système défini par la représentation d’état discrète suivante :
   
(
x(k + 1) = Ax(k) + Bu(k) 0.5 1 0 0 1  
1 0 1
, A = 0 0.1 0 , B = 1 1 , C =
   
y(k) = Cx(k) 0 1 1
0 0 0.5 0 −1

avec x(k) = [x1 (k) x2 (k) x3 (k)]T , u(k) = [u1 (k) u2 (k)]T , y(k) = [y1 (k) y2 (k)]T .
I. Auto-redondance
1. Déterminer la matrice d’observabilité réduite par rapport à la sortie y1 . Modifier la matrice
précédente pour obtenir de la redondance. Trouver une relation entre x(k), u(k), u(k + 1) et
la sortie y1 .
2.5 Exercices d’application 31

2. Déduire la relation d’auto-redondance sur le premier capteur y1 , r1 (k).


3. Déterminer la matrice d’observabilité réduite par rapport à la sortie y2 . Modifier la matrice
précédente pour obtenir de la redondance. Trouver une relation entre x(k), u(k), u(k + 1) et
la sortie y2 .
4. Déduire la relation d’auto-redondance sur le second capteur y2 , r2 (k).
5. Déterminer le vecteur d’auto-redondance r(k).
6. Déterminer la table des signatures des défauts si on considère 2 défauts fy1 et fy2 liés aux
capteurs y1 et y2 , et 2 défauts fu1 et fu2 liés aux actionneurs d’entrées u1 et u2 . Interpréter les
résultats (surtout dans le cas des défauts actionneurs).
7. Donner la table d’orientation du vecteur de parité en fonction des défauts. Le problème
précédent est-il résolu ?
8. Tracer dans l’espace de parité (le plan (r1 , r2 )) les différentes directions, susceptibles d’être
occupées par le vecteur de parité, en fonction des défauts.
II. Inter-redondance
1. Sachant que les matrices d’observabilité réduite par rapport à y1 et y2 sont de rang 2, déduire
4 relations indépendantes et écrire une relation d’inter-redondance.
2. En éliminant le vecteur des états de l’équation précédente, déduire une relation d’inter-
redondance r3 en fonction de y1 , y2 et u1 . Interpréter ce résidu (surtout vis à vis l’actionneur
d’entrée u2 ).
3. Regrouper ce résidu avec les 2 autres (obtenus par les relations d’auto-redondance). Déduire
la table des signatures des défauts dans ce cas. Est-il possible de localiser maintenant l’origine
de chaque défaut ?
4. Donner la table d’orientation du vecteur de parité en fonction des défauts.

Exercice 2.6 Simulation MATLAB/SIMULINK


Considérons le système décrit dans l’exercice précédent. Nous disposons des fichiers de me-
sures reflétant l’état du système dans des différents modes de fonctionnement sur un intervalle de
temps = 50 secondes. Chaque fichier de mesure contient les informations suivantes enregistrées
dans l’ordre : k, u1 , u2 , y1 , y2 où k est le numéro de l’échantillon (la période d’échantillonnage
= 1 seconde), (u1 , u2 ) sont les composantes du vecteur de commande u(k) et (y1 , y2 ) sont les
composantes du vecteur de sortie y(k).

R Utilisez les lignes de commande MATLAB suivants pour lire k, u et y à partir d’un fichier
de mesures mesi pour i = 1 · · · 5 :
load mes5 ;
mes=mes5 ;
[n, m]=size(mes) ;
k=mes(1, :) ;
u( :,1)=mes( :,2) ;
u( :,2)=mes( :,3) ;
y( :,1)=mes( :,4) ;
y( :,2)=mes( :,5) ;

I. Génération de résidu
1. Par analogie avec l’exemple théorique précédent, construire un générateur de résidus en
exploitant les relations de redondance dynamique (auto-redondance et inter-redondance).
2. Tracer ces différentes résidus sur la même figure. (Utiliser subplot).
3. Tester votre algorithme avec des différentes mesures.
II. Evaluation de résidu
1. En fonction des résidus tracés, choisir un seuil T pour que les résidus soient insensibles
aux bruits.
32 Chapitre 2. Espace de parité

2. Déterminer la signature de chaque résidu.


3. A partir de la table des signatures, construire un système de localisation des défauts.
4. Tester votre algorithme avec des différentes mesures.

3. Diagnostic à base d’observateurs

3.1 Introduction
Le problème de l’estimation de l’état d’un système est d’une importance pratique considérable,
que se soit pour la mise en oeuvre d’une loi de commande ou pour l’élaboration d’une stratégie de
diagnostic. Le principe de base de la génération de résidus à l’aide d’observateurs est de réaliser
une estimation des sorties du système à partir des grandeurs accessibles au calculateur (les entrées
et les sorties). Le vecteur résidu est alors construit comme l’écart entre la sortie estimée notée ŷ
et la sortie mesurée y permettant d’établir l’erreur de l’estimation sur la sortie ey . Un observateur
d’état, appelé aussi reconstructeur d’état ou estimateur est un algorithme ayant comme entrées, les
commandes et les mesures du processus réel, et dont la sortie est une estimation de l’état du système
si l’on souhaite faire de la commande, ou le vecteur des résidus si l’on souhaite faire du diagnostic.
Le principe général d’un observateur est présenté à la F IGURE 3.1, il consiste essentiellement en un
modèle de simulation du processus bouclé par l’erreur d’estimation.

3.2 Observateur de Luenberger


Considérons le système à surveiller, supposé correctement décrit par la représentation d’état
suivante :
(
ẋ(t) = Ax(t) + Bu(t) + Fx f (t)
(3.1)
y(t) = Cx(t) + Fy f (t)

où Fx et Fy sont les matrices d’action des défauts f (t) à détecter. Les grandeurs et matrices
suivantes x(t), u(t), y(t), A, B, et C ont la signification habituelle. Si la paire (A,C) est observable, il
est possible de reconstruire l’état x(t) à l’aide d’un modèle du système nominal, corrigé par l’écart
entre la sortie mesurée et la sortie estimée.
L’observateur d’état de Luenberger, utilisé en tant que générateur de résidus ry (t), est donné
34 Chapitre 3. Diagnostic à base d’observateurs
f(t) d(t)
défauts perturbations

Système supposé correctement modélisé par la


u(t) représentation d’état : y(t)

Entrées ‫ݔ‬ሶ ሺ‫ݐ‬ሻ = ‫ݔܣ‬ሺ‫ݐ‬ሻ + ‫ݑܤ‬ሺ‫ݐ‬ሻ + ‫ܨ‬௫ ݂ሺ‫ݐ‬ሻ + ‫ܦ‬௫ ݀ሺ‫ݐ‬ሻ, ‫ݔ‬ሺ0ሻ = ‫ݔ‬଴  Sorties

‫ݕ‬ሺ‫ݐ‬ሻ = ‫ݔܥ‬ሺ‫ݐ‬ሻ + ‫ܨ‬௬ ݂ሺ‫ݐ‬ሻ

Reconstructeur d’état
Modèle de simulation du système nominal

‫ݔ‬ොሶ (‫)ݐ‬ න ‫ݔ‬ො(‫)ݐ‬ ‫ݕ‬ො(‫)ݐ‬ ry(t)


B C
Résidus

L A

݁௬ ሺ‫ݐ‬ሻ = ‫ݕ‬ሺ‫ݐ‬ሻ − ‫ݕ‬ො(‫) ݐ‬

F IGURE 3.1 – Principe de la génération de résidus à l’aide d’un observateur.

par :
˙

 x̂(t) = Ax̂(t) + Bu(t) + L(y(t) − ŷ(t)),
 x̂(0) = x̂0
ŷ(t) = Cx̂(t) (3.2)

ry (t) = ey (t) = y(t) − ŷ(t)

La matrice L (gain de l’observateur) est calculé de façon que l’estimation tende vers l’état x(t) du
système quand t tend vers l’infini, quels que soient les états initiaux x(0) et x̂(0). La dynamique de
l’erreur de l’estimation sur l’état ex (t) = x(t) − x̂(t), s’écrit :
˙ = Ax(t) + Bu(t) + Fx f (t) − x̂(t) − Bu(t) − Ly(t) + Lŷ(t)
ėx (t) = ẋ(t) − x̂(t)
avec, y(t) = Cx(t) + Fy f (t) et ŷ(t) = Cx̂(t), ce qui donne finalement :
ėx (t) = (A − LC)ex (t) + (Fx − LFy ) f (t) (3.3)
En l’absence de défaut ( f (t) = 0), l’erreur d’estimation devient ėx = (A − LC)ex . On veut que
limt→∞ ex (t) = 0, ceci est le cas en calculant L telle que la matrice (A − LC) soit de Hurwitz.

R En temps discret, l’observateur de Luenberger est donné par


˙

x̂k+1 = Ax̂k + Buk + L(yk − ŷk ), x̂(0) = x̂0

ŷk = Cx̂k (3.4)

ryk = eyk = yk − ŷk

Les pôles de (A − LC) doivent être placés à l’intérieur du cercle unité.

3.3 Exercises d’application


Exercice 3.1 Soit le système décrit par la représentation d’état suivante :
(
ẋ = Ax + Bu
(3.5)
y = Cx + F f
   
−10 −1 10
où A = ,B= , C = F = I2 .
2 −1 0
3.3 Exercises d’application 35

1. Notons L le gain de notre observateur, une matrice diagonale. Quel est le degré (l’ordre) du
système ? déduire le nombre des éléments de L.
2. Donner le polynôme caractéristique de cet observateur ρobs (s), en fonction des éléments de
L.
3. Utilisons la technique de placement de pole, ajuster les paramètres de L d’une manière à ce
que ρobs (s) admet une racine double (-20).
4. Donner la représentation d’état de cet observateur, en fonction de x̂ et ŷ.
5. Donner la matrice de transfert G f (s) entre le vecteur d’erreur ey (s) et le vecteur des défauts
f (s).
6. Déduire la table des signatures associée à cette matrice, G f . Es ce qu’on peut localiser
l’origine d’un défaut à l’aide de cette table des signatures ?
7. La structure obtenue en (6.) est dite non localisante. Afin de rendre cette table des signature
localisante, chercher une matrice Q(s) telle que Q(s)G f (s) soit diagonale.
8. Déduire l’équation de notre générateur de résidus, r(s), en fonction de ey (s) et Q(s).
Exercice 3.2 Observateur "deadbeat"    
−2 1 1
On considère un système discret xk+1 = Axk + Buk avec A = , B= ,C=
0 −1 1
 
1 3 .
1. Est-ce qu’un observateur trivial (Luenberger) fonctionne pour ce système ?
2. Déterminer un observateur ayant la propriété que l’état estimé x̂k correspond exactement à xk
après deux coups d’horloges, indépendamment des conditions initiales.
3. Vérifier le comportement
  avec un signal d’entrée u = [1, 2, −1, 3, · · · ] et avec la condition
2
initiale x0 = . La condition initiale de l’estimateur x̂0 est supposée nulle.
−1
Exercice 3.3 Étude d’un bras manipulateur à 1 ddl
Considérons le bras manipulateur d’un robot, supposé rigide et tourne autour de x, comme
illustré à la F IGURE 3.2 Où x1 est la position angulaire du bras, par rapport à la verticale z, et u est
la variable de commande de couple. En appliquant le principe fondamental de la dynamique par
rapport à la rotation, nous montrerons que l’équation différentielle de mouvement est de la forme :

J ẍ1 + bẋ1 + mgl sin(x1 ) = u

où m est la masse du bras, J son moment d’inertie par rapport à l’axe de rotation, l est la distance du
centre de gravité à l’axe de rotation, b est le coefficient de frottement visqueux et g est l’accélération
due à la gravité.
1. On pose x2 = ẋ1 , Montrer que le modèle d’état non linéaire est donné par :
(
ẋ1 = x2
ẋ2 = 1J u − bJ x2 − mgl
J sin(x1 )

2. On donne la linéarisation de la fonction sin(x) au point π4 est xcos( pi4 ), déduire le modèle
d’espace d’état linéarisé. √
3. On donne : J = 0.25, g = 10, l = 0.5,m = 7 202 , b = 0. Les unités sont en S.I. Calculer A et
B, sachant que :
 " #

X =
 x 1




 x2
Ẋ = AX + Bu




 y = CXh i

C = 1 0

36 Chapitre 3. Diagnostic à base d’observateurs

(P)

u
z
l

x1
mg

F IGURE 3.2 – Un bras manipulateur à 1ddl.

4. On souhaite asservir la sortie y à la grandeur de consigne r, au moyen d’une loi de commande


par retour d’état de la forme :
 
  x1
u = kr r − k1 k2
x2

Donner le polynôme caractéristique du système bouclé ρ(s), sans tenir compte de kr (on ne
tient compte de kr qu’en régime permanent).
5. Déterminer k1 et k2 qui donnent au système en boucle fermée une dynamique de second
ordre de la forme ρd (s) = s2 + 2ζ ω0 s + ω02 avec ζ = 0.5 et ω0 = 50rd/s.
6. Déterminer kr tel qu’en régime permanent, y = r.
7. Pour surveiller notre système, on va servir d’un observateur de Luenberger. Donner sa
représentation d’état.
 >
8. Supposons que le gain L de l’observateur est de forme L = l1 l2 . Donner le polynôme
caractéristique de cet observateur ρO (s).
9. Déduire l1 et l2 si les pôles de cet observateur
 sont placés à λ1 = λ2 = −10.
  0 0
10. On a Fy = 1 0 , Fx = . Calculer la matrice de transfert reliant les défauts aux
0 1
résidus.
11. Donner la table de signatures des défauts. Ces défauts sont-ils localisables ?

Exercice 3.4 Simulation MATLAB/SIMULINK


Soit le système décrit par la représentation d’état suivante :


ẋ(t) = Ax(t) + Bu(t)

y(t) = Cx(t) + Fy f (t)

x(t) = [x1 (t) x2 (t)]T , y(t) = [y1 (t) y2 (t)]T , f (t) = [ f1 (t) f2 (t)]T ,

3.3 Exercises d’application 37
       
−10 −1 10 1 0 1 0
A= , B= , C= , Fy =
2 −1 0 0 1 0 1
I. Construction de l’observateur
1. Réécrire la représentation d’état précédente sous la forme augmentée suivante :
(
ẋ = A1 x + B1U
Y = C1 x + D1U

Exprimer les matrices A1 , B1 , C1 , D1 en fonction de A, B, C, et Fy et les vecteurs U et Y ,


en fonction de u, y et f .
2. Notons L le gain de notre observateur. Calculer ce gain si l’observateur admet deux valeurs
propres (-5) et (-9).
N.B. Utiliser le code MATLAB :
K = PLACE(A,B,P) Calcule K tel que P est le vecteur propre de A-B*K.
3. Construire un modèle SIMULINK qui simule le comportement de cet observateur :
4. Donner des échelons aux entrées de commande et des zeros aux défauts, puis simuler le
fichier sur un intervalle de 10 secondes. Donner l’évolution des sorties et des erreurs.
II. Construction du générateur de résidus
1. Dans un premier cas, supposons qu’un biais de 0.2 est apparu sur l’entrée f1 à l’instant 4
secondes. Tracer l’évolution des sorties et des erreurs. Interpréter les résultats.
2. Dans un second cas, supposons qu’un biais de -0.3 est apparu sur l’entrée f2 à l’instant 6
secondes. Tracer l’évolution des sorties et des erreurs. Interpréter les résultats.
3. On veut localiser les défauts, chercher une matrice Q tels que les défauts soient localisables.
Tester de nouveau votre simulateur.

II
Techniques de diagnostic
robustes

4 Observateur à entrée inconnue . . . . . . 41


4.1 Introduction
4.2 Décoplage de l’entrée inconnue
4.3 Structuration de résidu
4.4 Exercices d’application

5 Placement des structures propres . . . . 47


5.1 Introduction
5.2 Synthèse de l’observateur
5.3 Exercices d’application

6 Le filtre de Kalman . . . . . . . . . . . . . . . . . . . 53
6.1 Introduction
6.2 Estimation d’état d’un système stochastique
6.3 Le système linéaire stochastique
6.4 Le filtre de KALMAN
6.5 Application au diagnostic
6.6 Exercices d’application

A Appendice . . . . . . . . . . . . . . . . . . . . . . . . . . 69
A.1 La somme de deux Gaussiennes
A.2 Dérivée matricielle

Références bibliographiques . . . . . . . . . 71
Livres
Articles
4. Observateur à entrée inconnue

4.1 Introduction
La génération robuste de résidus est la tâche la plus importante dans le diagnostic de défauts.
En fait, cette approche est généralement basée sur le découplage de l’effet des perturbations ou des
incertitudes de modélisation sur les performances d’estimation d’état. Ces contraintes (perturbation,
erreur de modélisation) sont considérées comme une entrée inconnue agissant sur la dynamique
du système linéaire. Bien que le vecteur d’entrée inconnue soit un signal inconnu, sa matrice de
distribution (sa direction) est supposée bien déterminée. Sur la base des informations fournies par
cette matrice de distribution, l’entrée inconnue (perturbation) peut être découplée du résidu afin de
garantir la robustesse du filtre.

4.2 Décoplage de l’entrée inconnue


Le principe de construction d’un observateur à entrées inconnues consiste à rendre l’erreur
d’estimation indépendante des perturbations non mesurables. Considérons le système à surveiller,
supposé correctement décrit par la représentation d’état suivante :
(
ẋ(t) = Ax(t) + Bu(t) + Fx f (t) + Dx d(t)
(4.1)
y(t) = Cx(t) + Fy f (t)
Dans le cas où le vecteur des entrées inconnues agit également sur le vecteur sortie, il est possible
moyennant une transformation linéaire, de se ramener à la forme ci-dessous. Comme illustrée à la
F IGURE 4.1, la structure de l’observateur généralement adoptée est la suivante :
(
ż(t) = Mz(t) + Nu(t) + Py(t)
(4.2)
x̂(t) = z(t) − Ly y(t)
où M, N, P et Ly sont des matrices de dimensions appropriées, qui vont être déterminées de façon
que l’estimé x̂(t) converge asymptotiquement vers l’état réel x(t) du système, malgré l’influence
des perturbations. L’erreur de construction d’état ex (t) = x̂(t) − x(t) s’écrit :
ex (t) = z(t) − Ly y(t) − x(t) = z(t) − (I + LyC)x(t) − Ly Fy f (t) (4.3)
42 Chapitre 4. Observateur à entrée inconnue

F IGURE 4.1 – Observateur à entrée inconnue.

en posant E = I + LyC, il vient :


ex (t) = z(t) − Ex(t) − Ly Fy f (t), (4.4)
Proposition 4.2.1 La dynamique de l’erreur d’estimation d’état s’écrit alors :

ėx (t) =Mex (t) + (ME + PC − EA)x(t) + (N − EB)u(t) + (MLy Fy + PFy


(4.5)
− EFx ) f (t) − Ly Fy f˙(t) − EDx d(t)
Démonstration. En dérivant l’équation (4.4), nous obtenons :
e˙x (t) = ż(t) − E ẋ(t) − Ly Fy f˙(t), (4.6)
Après la substitution de ż et ẋ par leurs expressions :
e˙x (t) = (Mz(t) + Nu(t) + Py(t)) − E(Ax(t) + Bu(t) + Fx f (t) + Dx d(t)) − Ly Fy f˙(t),
= Mz(t) − EDx d(t) + (N − EB)u(t) + (PC − EA)x(t) + (PFy − EFx ) f (t) − Ly Fy f˙(t),
Soit aussi z(t) = ex (t) + Ex(t) + Ly Fy f (t) alors
ėx (t) = M(ex (t) + Ex(t) + Ly Fy f (t)) − EDx d(t) + (N − EB)u(t) + (PC − EA)x(t)
+ (PFy − EFx ) f (t) − Ly Fy f˙(t),
= Mex (t) + (ME + PC − EA)x(t) + (N − EB)u(t) + (MLy Fy + PFy
− EFx ) f (t) − Ly Fy f˙(t) − EDx d(t)


Si les conditions suivantes sont remplies :



 M est Hurwitz




 ME + PC = EA


N = EB
(4.7)


 EDx = 0
MLy Fy + PFy − EFx 6= 0





Ly Fy 6= 0

4.3 Structuration de résidu 43

alors, la dynamique de l’erreur d’estimation devient indépendante de l’état, de l’entrée de commande


et de l’entrée inconnue, elle n’est sensible qu’aux défauts :

ėx (t) = Mex (t) + (MLy Fy + PFy − EFx ) f (t) − Ly Fy f˙(t) (4.8)

La résolution du système (4.7) consiste, en premier lieu, à assurer la condition de découplage


des entrées inconnues, c’est-à-dire à satisfaire EDx = 0. Or, E = I + LyC, d’où

(I + LyC)Dx = 0 ⇒ LyCDx = −Dx

Il s’agit donc de déterminer la matrice Ly telle que la relation ci-dessus soit satisfaite. Si l’inverse
généralisée de CDx , notée (CDx )† , existe, Ly peut être calculée à l’aide de la relation suivante :

Ly = −Dx (CDx )† , avec (CDx )† = [(CDx )> (CDx )]−1 (CDx )> (4.9)

La matrice Ly n’existe que si la matrice (CDx )> (CDx ) est inversible. Cette matrice étant de
dimension nd × nd , elle n’est inversible que si rang(CDx ) = nd , où nd représente le nombre d’entrée
inconnues. Finalement, le découplage n’est possible que si le rang de la matrice (CDx ) est égal au
nombre d’entrées à découpler. La synthèse de l’observateur peut être résumée comme suit :
1. Vérifier que le rang(CDx ) = nd , puis calculer Ly = −Dx [(CDx )> (CDx )]−1 (CDx )> .
2. A partir de Ly calculer E = I + LyC.
3. A partir de E calculer N = EB.
4. Imposer que M soit une matrice de Hurwitz. On peut à cet effet choisir M une matrice
diagonale faisant apparaître les valeurs propres désirées pour l’observateur.
5. Calculer la matrice P telle que : PC = EA − ME.
Exercice 4.1 En se basant sur la démarche précédente, concevoir un filtre à entrée inconnue
pour le système suivant (équation de mesure perturbée) :
(
ẋ(t) = Ax(t) + Bu(t) + Fx f (t) + Dx d(t)
(4.10)
y(t) = Cx(t) + Du u(t) + Fy f (t) + Dy d(t)


4.3 Structuration de résidu


Calculons maintenant la matrice de transfert reliant les défauts à l’erreur d’estimation en sortie.
Posons F = MLy Fy − EFx et F 0 = −Ly Fy . La transformation de Laplace de l’expression (4.8)
s’écrit alors :

sex (s) = Mex (s) + F f (s) + F 0 s f (s) → ex (s) = (sI − M)−1 (F + F 0 ) f (s) (4.11)

L’erreur d’estimation en sortie s’écrit :

ey (s) = ŷ(s) − y(s) = C(x̂(s) − x(s)) − Fy f (s) = Cex (s) − Fy f (s)

en remplaçant ex (s) par son expression (4.11), on obtient :


−1 0
ey (s) = [C(sI − M) (F + sF − Fy ] f (s)


F = MLy Fy + PFy − EFx (4.12)
F 0 = −Ly Fy


44 Chapitre 4. Observateur à entrée inconnue

Soit Q(s) une matrice de transfert stable et propre et générons un vecteur de résidus r(s) tel
que :
(
r(s) = Q(s)ey (s) = Q(s)G f (s) f (s)
(4.13)
G f (s) = C(sI − M)−1 (F + sF 0 ) − Fy
La matrice de paramétrisation Q(s) permet de structurer les résidus afin de faciliter la localisation
des défauts.

4.4 Exercices d’application


Exercice 4.2 Banc d’observateurs à entrée inconnue
Soit un système linéaire décrit par la représentation d’état suivante :
(
ẋ = Ax + Bu(t) + F f
y = Cx
Les matrices sont définies comme suit :
       
0 −1 0 1 1 0 0 1 0 0
A = 100 −10 −1 , B = 0 ,
   C = 0 1 0 , F = 0 0 1
0 2 −1 0 0 0 1 0 1 0
1. Donner les dimensions de x, y, u, f .
2. A partir de cette représentation, nous pouvons construire les trois modèles équivalents
suivants :
(
ẋ = Ax + Bu(t) + F̄i f̄i + Fi fi
pour i = 1, · · · , dim( f )
y = Cx

ou fi est le ième élément du vecteur f et f̄i est composé des éléments restants.
3. Évaluer F̄i , Fi .
4. Pour chacune de représentations précédentes, nous cherchons à concevoir un filtre à entrée
inconnue permettant d’éliminer l’effet de fi (découplage) et de garder l’effet de f̄i (sensibilité).
5. Vérifier la condition d’application du filtre à entrée inconnue pour chaque modèle.
(i)
6. Trouver les matrices Ly , Ei , Ni .
7. En choisissant Mi (Matrices de Hurwitz diagonale), déduire Pi .
8. Donner la représentation d’état de chaque filtre.
9. En implémentant ces observateurs simultanément, nous obtenons une architecture de sur-
veillance connue par le banc d’observateurs. Quelle est l’utilité de cette architecture ?

Exercice 4.3 Simulation MATLAB/SIMULINK


Soit le système décrit par la représentation d’état suivante :

     
 ẋ1 (t) x1 (t) f1 (t)
ẋ2 (t) = A x2 (t) + Bu(t) + Fx  f2 (t) + Dx d(t)





 ẋ3 (t)

x3 (t) f3 (t)
     
 y1 (t)
 x1 (t) f1 (t)
F





y 2 (t) = C  x2 (t) + y  f 2 (t)
y3 (t) x3 (t) f3 (t)

4.4 Exercices d’application 45

Les divers matrices sont les suivants :


       
0 −1 0 1 1 0 0 0
A = 100 −10 −1 , B = 0 ,
   C = 0 1 0 , Dx = 0 ,

0 2 −1 0 0 0 1 1
   
1 0 0 0 0 0
Fx = 0 0 0 , Fy = 0 1 0
0 0 0 0 0 1
Nous voulons construire un observateur à entrées inconnues qui possède les valeurs propres
suivantes : (λ1 = −5, λ2 = −6, λ3 = −7).
1. Recalculer Ly , E, N, M, et P.
2. Donner la représentation d’état de cet observateur (représentation augmentée).
3. Donner des échelons aux entrées de commande et des zeros aux défauts, puis simuler le
fichier sur un intervalle de 10 secondes. Donner l’évolution des sorties et des erreurs.
II. Construction du générateur de résidus
1. Par simulation, montrer que les perturbations sont découplées.
2. Simuler des différents défauts à des différents instants puis tracer l’évolution des sorties
et des erreurs. Interpréter les résultats. Localisation des défauts ?
3. On veut localiser les défauts, chercher une matrice Q qui permet de localiser les défauts.
Tester de nouveau votre simulateur.

5. Placement des structures propres

5.1 Introduction
Dans le chapitre précédent, une approche robuste basée sur l’observateur à entrée inconnue
a été introduite. Le principe de base de cette approche est de rendre l’erreur d’estimation d’état
indépendant des perturbations (entrées inconnues). Le résidu est défini comme étant l’erreur
d’estimation de la sortie (pondérée) qui est une transformation linéaire de l’erreur d’estimation
d’état. Le résidu généré par les observateurs à entrée inconnue est également indépendant des
perturbations, si le terme de perturbation n’apparaît pas dans l’équation de sortie ou si le terme de
perturbation dans l’équation de sortie a été annulé. Dans le diagnostic basé modèle, l’estimation
d’état n’est pas la mission principale de l’algorithme conçu, car l’information requise est le signal de
indicateur de défaut. Par conséquent, il n’est pas nécessaire de découpler l’erreur d’estimation d’état
des perturbations. Une approche directe de découplage de l’entrée inconnue est alors nécessaire.
Dans le contexte de placement des structures propres développé dans ce chapitre, le résidu lui-
même est découplé des perturbations, mais l’erreur d’estimation d’état peut ne pas l’être. On peut
s’attendre que les conditions d’application pour une telle approche directe est largement assouplies
par rapport à celles requises pour les observateur à entrée inconnue. L’approche directe la plus
importante pour concevoir des générateurs résiduels robustes (dans le sens de découplage des
perturbations) est de forcer certains vecteurs propres gauche de l’observateur à être orthogonaux
aux directions de distribution des perturbations. De cette manière, le résidu peut être rendu robuste
vis-à-vis ces perturbations.

5.2 Synthèse de l’observateur


Considérons le système à surveiller, supposé correctement décrit par la représentation d’état
suivante :
(
ẋ(t) = Ax(t) + Bu(t) + Fx f (t) + Dx d(t)
(5.1)
y(t) = Cx(t) + Fy f (t)
48 Chapitre 5. Placement des structures propres

où Dx représente la matrice d’action des perturbations d(t), Fx et Fy sont les matrices d’action
des défauts f (t) à détecter. Les grandeurs et matrices suivantes x(t), u(t), y(t), A, B, et C ont la
signification habituelle. Si la paire (A,C) est observable, il est possible de reconstruire l’état x(t)
à l’aide d’un modèle du système nominal, corrigé par l’écart entre la sortie mesurée et la sortie
estimée.
L’observateur d’état de Luenberger, utilisé en tant que générateur de résidus ry (t), est donné
par :
˙

 x̂(t) = Ax̂(t) + Bu(t) + L(y(t) − ŷ(t)), x̂(0) = x̂0

ŷ(t) = Cx̂(t) (5.2)

ry (t) = ey (t) = y(t) − ŷ(t)

La matrice L (gain de l’observateur) est calculé de façon que l’estimation tende vers l’état x(t) du
système quand t tend vers l’infini, quels que soient les états initiaux x(0) et x̂(0). La dynamique de
l’erreur de l’estimation sur l’état ex (t) = x(t) − x̂(t), s’écrit :
˙ = Ax(t) + Bu(t) + Fx f (t) + Dx d(t) − x̂(t) − Bu(t) − Ly(t) + Lŷ(t)
ėx (t) = ẋ(t) − x̂(t)
or, y(t) = Cx(t) + Fy f (t) et ŷ(t) = Cx̂(t), d’où finalement :
ėx (t) = (A − LC)ex (t) + (Fx − LFy ) f (t) + Dx d(t) (5.3)
En l’absence de défaut ( f (t) = 0) et en négligeant l’effet des entrées inconnues (d(t) = 0), l’erreur
d’estimation devient ėx = (A − LC)ex . On veut que limt→∞ ex (t) = 0, ceci est le cas en calculant L
telle que la matrice (A − LC) soit de Hurwitz.

5.2.1 Structuration de résidu

F IGURE 5.1 – Structuration de résidu à l’aide d’un post-filtre.

Intéressons nous maintenant à la matrice de transfert reliant les diverses entrées à l’erreur
d’estimation d’état. LA transformée de Laplace de l’expression (5.3) s’écrit (les conditions initiales
étant nulles par définition) :
sex (s) = (A − LC)ex (s) + (Fx − LFy ) f (s) + Dx d(s)
d’où l’expression de l’erreur d’estimation d’état :
ex (s) = [sI − (A − LC)]−1 (Fx − LFy ) f (s) + [sI − (A − LC)]−1 Dx d(s) (5.4)
5.2 Synthèse de l’observateur 49

Notons que cette erreur est sensibles aux défauts f (s) et pourrait donc être utilisée pour générer
un vecteur indicateur de défauts. Toutefois, l’état étant inconnu, cet écart ne peut être exploité
directement. Considérons alors l’erreur d’estimation en sortie ey (s), elle s’écrit :

ey (s) = y(s) − ŷ(s) = Cx(s) + Fy f (s) −Cx̂(s) = Cex (s) + Fy f (s)

En reportant dans ey (s) l’expression de ex (s) donnée par la relation (5.4), on obtient :

 ey (s) = G f (s) f (s) + Gd (s)d(s)

G f (s) = C[sI − (A − LC)]−1 (Fx − LFy ) + Fy (5.5)

Gd (s) = C[sI − (A − LC)]−1 Dx

On peut constater que cette erreur d’estimation est sensible aux défauts f (s), mais généralement
aux entrées inconnues d(s). Il est dès lors possible d’exploiter cet écart pour engendrer un ensemble
de signaux, sensibles aux défauts, appelés vecteur des résidus. Soit Q(s) une matrice de transfert
stable est propre et générons un vecteur de résidus tel que :

r(s) = Q(s)ey (s), soit aussi r(s) = Q(s)G f (s) f (s) + Q(s)Gd (s)d(s) (5.6)

Dans ces conditions, si l’on peut trouver une matrice de paramétrisation Q(s) telle que les deux
relations suivantes soient vérifiées :

Q(s)G f (s) 6= 0, Q(s)Gd (s) = 0 (5.7)

On obtiendra un vecteur des résidus r(s) insensibles aux perturbations d(s) et sensibles aux défauts
f (s). Dans le cas ou l’influence des perturbations est négligeable sur l’évolution des résidus, la
matrice Q(s) peut être entièrement utilisée pour faciliter la localisation des défauts, on parle alors
de structuration des résidus.

5.2.2 Placement des structures propres


Soit une matrice de paramétrisation Q(s) = Q f (s)Qd , où Q f est une matrice de transformation
stable et propre destinée à la structuration des résidus et Qd une matrice algébrique destinée au
découplage des entrées inconnues. D’après les relations (5.5) et (5.6), le vecteur des résidus s’écrit :


 r(s) = Q f (s)Qd G f (s) f (s) + Q f (s)Qd Gd (s)d(s)

G f (s) = C[sI − (A − LC)]−1 (Fx − LFy ) + Fy (5.8)

Gd (s) = C[sI − (A − LC)]−1 Dx

le problème est de déterminer Qd telle que Qd Gd (s) = 0, soit :

Qd C[sI − (A − LC)]−1 Dx = 0 (5.9)

Soit wTi le vecteur propre à gauche associé à la valeur propre λi de (A − LC). Ce vecteur propre
est toujours orthogonal aux vecteurs propres à droite v j associés aux nx − 1 valeurs propres λ j
de(A − LC), où λi 6= λ j . Nous écrivons pour a matrice (sI − Ac )−1 , avec Ac = A − LC, la relation
suivante :
v1 wT1 vn wTnx
(sI − Ac )−1 = +···+ (5.10)
s + λ1 s + λnx

où vi et wTi sont respectivement, les vecteurs propres à droite et à gauche de Ac associés à la


valeur propre λi . Montrons que si Qd CDx = 0 et que toutes les lignes de la matrice Wq = Qd C sont
50 Chapitre 5. Placement des structures propres

les vecteurs propres à gauche de Ac associés à q valeurs propres de Ac , alors la relation (5.9) est
satisfaite. en effet, si l’on construit Qd telle que :
 T
w1
wT 
 2
Qd CDx = 0 et Wq = Qd C =  .. 
 . 
wTq
on obtient Wq vi = 0 pour i = q + 1, . . . , n d’où :

v1 wT1 vn wTn
 
−1
Wq (sI − Ac ) Dx = Wq +···+ Dx
s + λ1 s + λn
Or, Qd CDx = Wq Dx = 0 que l’on peut aussi écrire wTi Dx = 0 pour i = 1, 2, . . . , q. Par conséquent, on
a bien Qd C(sI − Ac )−1 Dx = 0. La procédure de placement de structures propres peut être résumée
comme suit :
1. Calculer la matrice Qd telle Qd CDx = 0
2. Déterminer la structure propre de l’observateur. Les n valeurs propres de Ac sont choisies en
fonction de la dynamique de reconstruction souhaitée. Les q lignes de Wq = Qd C doivent être
q vecteurs propres à gauche associés à q valeurs propres de Ac = A − LC, les n − q vecteurs
propres restants peuvent êtres fixés arbitrairement.
3. Calculer la matrice gain L. Il s’agit de déterminer L de manière à imposer à l’observateur la
structure propre désirée. On peut à cet effet utiliser la méthode de placement de structure
propre en résolvant un problème de commande sur le système dual (voir annexe). Dans ces
conditions, le placement des vecteurs propres à gauche de l’observateur est transformé en un
placement des vecteurs propres à droite du système commandé dual.
Ce problème peut également être résolu à l’aide d’un placement de vecteurs propres à droite de
Ac . On montrerait en effet de la même façon que si QdCDx = 0 et que les colonnes de Dx sont les
vecteurs propres à droite associés à nc valeurs propres de Ac alors la relation (5.9) est satisfaite.
Bien entendu, le nombre d’entrées découplées diminue d’autant le nombre de résidus indé-
pendants qu’il est possible de générer. De façon plus précise, pour une matrice Dx donnée, le
découplage des perturbations peut être réalisé si et seulement si la condition suivante est satisfaite :

rang(Dx )  nx − q (5.11)

5.3 Exercices d’application


Exercice 5.1 Soit le système décrit par la représentation d’état suivante :

     
 ẋ1 (t)
 x1 (t) f1 (t)
ẋ (t) = A x2 (t) + Bu(t) + Fx  f2 (t) + Dx d(t)
 2



 ẋ3 (t) x3 (t) f3 (t)
     

 y1 (t) x1 (t) f1 (t)
y2 (t) = C x2 (t) + Fy  f2 (t)





y3 (t) x3 (t) f3 (t)

Les divers matrices sont les suivants :

       
0 −1 0 1 1 0 0 0
A = 100 −10 −1 , B = 0 , C = 0 1 0 , Dx = 0 ,
0 2 −1 0 0 0 1 1
5.3 Exercices d’application 51

   
1 0 0 0 0 0
Fx = 0 0 0 , Fy = 0 1 0
0 0 0 0 0 1

I. Construction de l’observateur
Nous voulons construire un observateur qui a la structure propre suivante (λ1 = −5, wT1 ),
(λ2 = −6, wT2 ) et (λ3 = −7, wT3 ), où les wti sont les vecteurs propres à gauches associées aux valeurs
propres λi
1. Quel est le rang de la matrice Dx ? déduire le nombre maximum de résidus indépendants que
l’on peut générer.  
q1 q2 q3
2. Soit Qd de la forme , chercher Qd telle que Qd CDx = 0. déduire wT1 et wT2 .
q4 q5 q6
 T
w1
3. Placer arbitrairement w3 tel que la matrice W = wT2  des vecteurs propres à gauches de
T 
wT3
l’observateur, soit de rang plein.
4. Si on veut réaliser le placement de cette structure propre conformément à la commande par
retour d’état sur le système dual :

(
ẋd = Ad xx + Bd ud
avec Ad = AT , Bd = CT
ud = −Kxd

où xd et ud représente l’état et la commande du système dual.


5. Calculer la matrice d’évolution du système bouclé par la loi de commande ud = −Kxd en
fonction de A,C, K.
6. Calculer K telle que la structure propre de ce système bouclé soit identique à notre observa-
teur.
7. Montrer que L = K T .
8. Déduire la représentation d’état de notre observateur.
II.Construction du générateur de résidus
1. Calculer les matrices de transfert Gd (s) et G f (s).
2. Déduire r(s) en fonction de Q f (s) et des défauts f1 (s), f2 (s) et f3 (s).
3. Afin de faciliter la localisation des défauts, nous souhaitons avoir une table des signatures de
la forme suivante :

f1 f2 f3
r1 1 0 1
r2 0 1 1

TABLE 5.1 – Table des signatures désirées

4. Donner la forme générale de Q f (s)Q


 d G f (s) 
q1 (s) q2 (s)
5. Supposons que Q f (s) est de forme .
q3 (s) q4 (s)
6. Déterminer un système à huit relations en fonction de q1 , q2 , q3 et q4 .
7. Donner une solution à ce système.
8. Déduire Q f (s).
9. Déduire l’équation du générateur de résidus, r(s), en fonction de f (s).
6. Le filtre de Kalman

6.1 Introduction
Après une étude statistique approfondie, il est possible d’établir un modèle probabiliste décrivant
le comportement moyen des signaux exogènes (perturbations et bruits). Ceci est réalisé en faisant
appel à la théorie des processus stochastiques ou aléatoires. En fait, l’état des systèmes dynamiques
subissant l’effet de ce type de signaux est lui même de nature stochastique et par conséquence
Ils sont aussi appelés systèmes stochastiques. Comme dans le cas déterministe, la génération des
résidus pour les systèmes stochastiques peut employer des algorithmes d’estimation d’état, et cela,
en affectant au résidu la même valeur que celle de l’erreur d’estimation de la sortie. L’estimation
dans le contexte stochastique consiste à trouver les propriétés statistiques de l’état en utilisant
l’ensemble des mesures fournies par les capteurs. Cette opération (estimation) est connue aussi
par le terme le filtrage d’état car elle permet de minimiser davantage l’effet néfaste des signaux
exogènes sur les performances de fonctionnement attendues. Pour un système linéaire stochastique
et Gaussien, le filtre de Kalman est l’algorithme le plus utilisé vue sa robustesse, sa simplicité et
ses hautes performances. Dans ce chapitre, nous mettons l’accent sur le filtre de Kalman et son
utilisation dans la surveillance des systèmes linéaires stochastiques.

6.2 Estimation d’état d’un système stochastique


Le filtrage consiste à trouver les propriétés statistiques de l’état en utilisant un ensemble de
mesures Ik en provenance des capteurs. Généralement, on cherche à trouver la fonction de densité
de probabilité (la PDF) de l’état conditionnée sur l’ensemble Ik ç-à-d :

p(xk | Ik ) (6.1)

Selon la composition de l’ensemble des mesures utilisées dans l’estimation d’état, nous distinguons
les trois cas particuliers suivants :
• Prédiction : Si Ik = Y0:k−n = {y0 , . . . , yk−n }
• Filtrage : Si Ik = Y0:k = {y0 , . . . , yk }
54 Chapitre 6. Le filtre de Kalman

uk
système

yk

Calcul de la p ( xk Y0:k ) x̂k état estimé


Critère
FDP
d’estimation
conditionnelle

Filtrage d’état

F IGURE 6.1 – Un schéma descriptif du filtrage d’état

• Lissage : Si Ik = Y0:k+m = {y0 , . . . , yk+m }

F IGURE 6.2 – L’estimation de la FDP conditionnelle de l’état

Une fois la fonction de densité de probabilité de l’état p(xk | Ik ) est déterminée, nous pouvons
calculer selon plusieurs critères considérés (F IGURE 6.3), une estimation optimale de la variable
d’état. Par exemple, l’estimation d’état à erreur quadratique moyenne minimale n’est que l’espérance
mathématique conditionnelle (EQMM) de l’état xk , ç-à-d :

x̂k = x̂kEQMM = arg min E (xk − x̃k )T (xk − x̃k )|Y0:K


 
(6.2)
x̃k

De même, l’estimation du maximum à posteriori (MAP), est obtenu à partir de la fonction de


6.3 Le système linéaire stochastique 55

densité comme suit


x̂k = x̂kMAP = arg max {p(xk | Y0:k )} (6.3)
xk

F IGURE 6.3 – Critères d’estimation d’état

6.3 Le système linéaire stochastique


Soit le système linéaire stochastique suivant (système de Gauss-Markov) :
(
xk+1 = Axk + Buk + wk
(6.4)
yk = Cxk + vk

où xk ∈ ℜnx est le vecteur d’état, uk ∈ ℜnu est la commande, yk ∈ ℜny est le vecteur de mesures.L’état
initial x0 , les perturbations wk et le bruit de mesure vk sont des processus aléatoires blancs non
corrélés avec
x0 ∼ Bb(x̄0 , P̄0 ), wk ∼ Bb(0,W ), vk ∼ Bb(0,V )

E(vi v j ) = V δ (i − j), E(wi w j ) = W δ (i − j), E(vi w j ) = 0


où P̄0 ,W  0 sont des matrices semi-définies positives et V  0 une matrice définie positive.

R Le bruit blanc Bb n’est pas nécessairement Gaussien, mais étant données la variance et la
moyenne des bruits, le filtre de Kalman reste le meilleur estimateur linéaire pour ce type de
système. Si le bruit est Gaussien, on remplace la notation Bb par N pour bien indiquer qu’il
s’agit d’une distribution gaussienne.

Proposition 6.3.1 Si le bruit et la perturbation affectant le système suivent une distribution Gaus-
sienne, alors la densité conditionnelle de l’état p(xk |Y0:k ) est aussi Gaussienne.

1 1
p(xk |Y0:k ) = N (xk , mk , Pk ) = det(2πPk )− 2 exp(− (xk − mk )Pk−1 (xk − mk )> ) (6.5)
2
Démonstration. La relation (6.5) est un résultat direct de la propriété de somme de deux Gaus-
siennes présentée dans l’appendice A.1. 
56 Chapitre 6. Le filtre de Kalman

6.4 Le filtre de KALMAN


Depuis sa publication en 1960, le filtre de Kalman est omniprésent dans une large gamme de
domaines technologiques (radar, vision électronique, communication ...). C’est un thème majeur de
l’automatique et du traitement du signal. Un exemple d’utilisation peut être la mise à disposition,
en continu, d’informations telles que la position ou la vitesse d’un objet à partir d’une série
d’observations relative à sa position, incluant éventuellement des erreurs de mesures. Le filtrage de
Kalman est aussi de plus en plus utilisé en dehors du domaine du traitement du signal, par exemple
en météorologie et en océanographie, pour l’assimilation de données dans un modèle numérique,
en finance ou en navigation et il est même utilisé dans l’estimation des états de trafic routier.
Cependant, pour le diagnostic des systèmes, le filtre de Kalman est utilisé comme un générateur de

F IGURE 6.4 – Le principe de fonctionnement de filtre de Kalman

résidu. En effet, l’erreur d’estimation de la sortie (innovation) est utilisée comme un indicateur de
défaut et qui est aussi un signal stochastique d’une moyenne nulle et une variance connue dans le
cas non-défaillant. La détection de changement dans la distribution de résidu permet de signaler
l’apparition d’un défaut (alarme). Comme illustrée à la F IGURE 6.4, le filtrage de Kalman consiste
à trouver à chaque instant k, la densité de probabilité conditionnelle de l’état xk . Cette dernière est
d’une distribution Gaussienne complètement décrite par sa moyenne mk et sa variance Pk (les seuls
paramètres inconnus).

F IGURE 6.5 – L’approche directe du filtrage de Kalman

Le filtre de Kalman est un estimateur récursif. Cela signifie que pour estimer la densité l’état
courant, seules l’estimation de l’état précédent et les mesures actuelles sont nécessaires. Donc,
6.4 Le filtre de KALMAN 57

pour trouver la densité de probabilité N (xk , mk , Pk ), il est possible de passer par les deux étape
distinctes : la Prédiction et la Correction. L’étape de prédiction utilise la densité conditionnelle
de l’état estimé de l’instant précédent, notée N (xk−1 , mk−1 , Pk−1 ) et le modèle du système (6.4)
afin de produire la densité conditionnelle de l’état prédit courant N (xk , m− −
k , Pk ). Dans l’étape de
correction, la mesure de l’instant courant yk est utilisée pour corriger la densité de l’état prédit
dans le but d’obtenir une densité conditionnelle plus précise, notée N (xk , mk , Pk ) (voir F IGURE
bayes). Après avoir explicité la densité de probabilité conditionnelle, l’estimation d’état choisie

F IGURE 6.6 – Calcul de la Gaussienne conditionnelle

dans l’algorithme de filtrage de Kalman est d’une erreur quadratique moyenne minimale (EQMM).
Dans la littérature, il existe deux approches différentes pour le développement de l’algorithme de
Kalman : l’approche directe (Bayésienne) et une approche indirecte (erreur quadratique moyenne
minimale). Ces deux approches feront l’objet des deux sous-sections suivantes.

6.4.1 L’approche indirecte (Bayésienne)


Dans la suite, nous développons la technique d’estimation Bayésienne d’une manière générale,
ensuite nous adaptons les résultats obtenus au système linéaire stochastique Gaussien (6.4). Suppo-
sons maintenant que p(xk−1 |Y0:k−1 ) la densité de l’état corrigé à l’instant k − 1 est connue. L’étape
de prédiction consiste à utiliser le modèle du système pour obtenir la densité prédite p(xk |Y0:k−1 )
via l’équation de Chapmann-Kolmogorov :
Z
p(xk |Y0:k−1 ) = p(xk |xk−1 )p(xk−1 |Y0:k−1 )dxk−1 (6.6)

A l’instant k et lorsque la mesure à yk devient disponible, nous procédons à la mise à jour de la


densité prédite en employant la règle de Bayes :
p(yk |xk , Y0:k−1 )p(xk |Y0:k−1 )
p(xk |Y0:k ) =
p(yk |Y0:k−1 )
p(xk |Y0:k−1 )p(yk |xk )
=
p(yk |Y0:k−1 )
p(xk |Y0:k−1 )p(yk |xk )
=R (6.7)
p(xk |Y0:k−1 )p(yk |xk )dxk
Donc, l’approche Bayésienne d’estimation de la densité conditionnelle de l’état du système se réside
dans l’évaluation des intégrales : ((6.6) pour la prédiction) et ((6.7) pour la correction). Considérons
maintenant le cas Gaussien ç-à-d p(xk |Y0:k−1 ) = N (xk , m− −
k , Pk ) et p(xk |Y0:k ) = N (xk , mk−1 , Pk ).
Le théorème suivant récapitule le résultat obtenu après l’évaluation des intégrales (6.6) et (6.7) pour
ce contexte particulier.
58 Chapitre 6. Le filtre de Kalman

Théorème 6.4.1 Pour le système (6.4), les relations récurrentes suivantes sont valables à chaque
instant k :

N (xk , m− − >
k , Pk ) = N (xk , Amk−1 + Buk−1 , APk−1 A +W ) (6.8)
N (xk , mk , Pk ) = N (xk , m− − −
k + Kk (yk −Cmk ), (In − KkC)Pk ) (6.9)
−1
où Kk = Pk−CT CPk−CT +V


Démonstration. Voir la page 56 de [1]. 


A partir du théorème 6.4.1, il est possible d’obtenir le filtre de Kalman (algorithme 1).

Algorithm 1 Filtre de Kalman


1. Initialisation :

m0 = x̄0
P0 = P̄0

2. Pour chaque pas d’échantillonnage k > 0 faire :


• Prédiction :

m−
k = Amk−1 + Buk−1
Pk− = APk−1 AT +W

• Corréction :
−1
Kk = Pk−CT CPk−CT +V


mk = m− −
k + Kk (yk −Cmk )
Pk = (In − KkC)Pk−

3. Fin Pour

Proposition 6.4.2 L’estimation d’état à erreur quadratique moyenne minimale (EQMM) n’est que
l’espérance mathématique conditionnelle de l’état xk ç-à-d :
Z
x̂k = x̂kEQMM = E(xk |Y0:k ) = xk p(xk | Y0:k )dxk (6.10)

Démonstration.
x̂kEQMM = arg min E (xk − x̃k )T (xk − x̃k )|Y0:K
 
x̃k

sachant que x̃k est une constante (et pas une variable aléatoire), montrons maintenant que
x̂kEQMM = arg min E (xk − x̃k )T (xk − x̃k )|Y0:K
 
x̃k

= arg min E(xk T xk |Y0:K ) + x̃kT x̃k − 2x̃kT x̂k



x̃k

= arg min{E(xk T xk |Y0:K ) + x̂kT x̂k + (x̂k − x̃k )T (x̂k − x̃k )}


x̃k | {z } | {z }
indépendant de x̃k terme>0

=⇒ x̂kEQMM = x̂k = E(xk |Y0:k )



6.4 Le filtre de KALMAN 59

D’après la proposition précédente, nous pouvons remplacer mk par x̂k et m− −


k par x̂k dans
l’algorithme 1 pour avoir la forme standard du filtre de Kalman.

Algorithm 2 Filtre de Kalman


1. Initialisation :

x̂0 = x̄0
P0 = P̄0

2. Pour chaque pas d’échantillonnage k > 0 faire :


• Prédiction :

x̂k− = Ax̂k−1 + Buk−1


Pk− = APk−1 AT +W

• Corréction :
−1
Kk = Pk−CT CPk−CT +V
x̂k = x̂k− + Kk (yk −Cx̂k− )
Pk = (In − KkC)Pk−

3. Fin Pour

6.4.2 L’approche directe (EQMM)


Dans l’approche directe, nous utilisons la proposition 6.4.2 pour déduire directement l’état
prédit x̂k− et estimé x̂k , et cela, sans avoir besoin de passer par la fonction de densité conditionnelle.
• La prédiction : D’après la proposition 6.4.2, nous pouvons définir l’état prédit (la moyenne de
l’état conditionnelle) ainsi que la matrice de covariance.

x̂k− = E (xk | Y0:k−1 )


Pk− = E (xk − x̂k− )(xk − x̂k− )T | Y0:k−1


∆ 

Supposons que l’état corrigé x̂k est connu, alors :



x̂k+1 = E(xk+1 | Y0:k )
= E(Axk + Buk + wk | Y0:k )
= E(Axk | Y0:k ) + E(Buk | Y0:k ) + E(wk | Y0:k )
| {z }
0
= Ax̂k + Buk

Dans la suite nous utilisons la notation suivante : soit une matrice quelconque M, le produit MM >
peut être écrit sous la forme : MM > = M(· · · )> .
La matrice de covariance conditionnelle Pk− ne dépend pas des observations alors

Pk− = E (xk − x̂k− )(· · · )T




= E ((Axk−1 + Buk−1 + wk−1 − (Ax̂k−1 + Buk−1 ))) (· · · )T


= E (A(xk−1 − x̂k−1 ) + wk−1 ) (· · · )T
60 Chapitre 6. Le filtre de Kalman

La perturbation wk n’est pas corrélée avec l’état, alors la matrice de covariance est donnée par
Pk− = APk−1 AT +W
• La correction : On corrige l’état prédit en utilisant la mesure yk comme suit :
x̂k = x̂k− + Kk (yk −Cx̂k− )
| {z }
innovation
et la covariance
Pk = (In − KkC)Pk− (In − KkC)T + KkV KkT (Forme de Joseph)
= Pk− − Pk−CT KkT − KkCPk− + Kk CPk−CT +V KkT

(6.11)
Le gain de Kalman Kk doit minimiser l’erreur quadratique moyenne e> >
k ek = (xk − x̂k ) (xk − x̂k ),
d’une manière équivalent, nous pouvons définir une fonction objective à minimiser J comme suit
nx
J = ∑ E e2i,k

i=1
D’autre part, on a l’expression de la matrice de covariance ci-dessous
E(e21,k )
 
E(e1,k e2,k ) · · · E(e1,k en,k )
 2 )
.. 
T
  E(e2,k e1,k ) E(e 2,k · · · . 
Pk = E ek ek =  .. .. .

 .. .. 
 . . . 
E(en,k e1,k ) ··· ··· 2
E(enx ,k )
La trace de la matrice de covariance n’est que la fonction objective J définie comme étant la somme
des erreurs quadratiques.
J = trace(Pk ) et Kk = arg min J
Kk

La fonction objective J est une fonction positive et convexe alors, dans la valeur minimale on a :

∂J
=0
∂ Kk
D’après la propriété A.2 de l’appendice, nous écrivons
∂J ∂trace(Pk )
= = 2(In − KkC)Pk− (−C> ) + 2KkV
∂ Kk ∂ Kk
∂J −1
= 0 =⇒ Kk = Pk−C> CPk−CT +V

∂ Kk
En remplaçant Kk dans l’équation de covariance nous obtenons :
Pk = (In − KkC)Pk−

R Le filtre de Kalman est un estimateur non biaisé c-à-d :


E(ek | Y0:k ) = 0 et E(e−
k | Y0:k−1 ) = 0
d’où
−T
Pk− = E (xk − x̂k− )(xk − x̂k− )T | Y0:k−1 = E(e− cov e−
 
k ek ) = k
| {z } | {z }
covariance de l’état covariance de l’erreur d’éstimation

Pour une application en temps réel des récurrences précédentes, nous pouvons considérer les mêmes
étapes décrites dans l’algorithme 2.
6.5 Application au diagnostic 61

6.5 Application au diagnostic


Considérons le système stochastique défaillant suivant
(
xk+1 = Axk + Buk + wk + Fx fk
yk = Cxk + vk + Fy fk
où fk ∈ ℜ f est le vecteur défaut. Pour surveiller le système, nous utilisons le processus d’innovation
comme un résidu indicateur de défaut :
rk = yk −Cx̂k− (6.12)
Le filtre de Kalman diagnostiquer est donné par l’algorithme 3. Le calcul du vecteur résidu est basé

Algorithm 3 Filtre de Kalman diagnostiquer


1. Initialisation :

x̂0 = x̄0
P0 = P̄0

2. Pour chaque pas d’échantillonnage k > 0 faire :


• Prédiction :

x̂k− = Ax̂k−1 + Buk−1


Pk− = APk−1 AT +W

• Génération de résidu :

rk = yk −Cx̂k−

• Corréction :
−1
Kk = Pk−CT CPk−CT +V
x̂k = x̂k− + Kk (yk −Cx̂k− )
Pk = (In − KkC)Pk−

3. Fin Pour

sur l’évaluation de l’état prédit x̂k− . Cependant, une version plus rapide du filtre de Kalman peut
être obtenue en remplaçant les équations de correction dans celle de correction. Ce filtre est appelé
le filtre prédiction car il ne fournie à chaque instant k que l’état prédit x̂k− et la covariance d’erreur
de prédiction Pk− (voir l’algorithme 4).

6.5.1 Conditions de convergence


Pour réduire la complexité temporelle du filtre de Kalman on peut utiliser une version station-
naire dont l’exécution est plus rapide. Pour cela, les conditions suivantes doivent être vérifiées avant
de passer à l’implémentation en temps réel.

Théorème 6.5.1 Les conditions de convergence du filtre de Kalman sont les suivantes :
1. (A,C) observable.
1 1 1T 1
2. (A,W 2 ) commandable (W 2 est la racine carrée de la matrice W telle que W 2 W 2 = W ).
Sous les conditions de convergence on a ( lorsque k → ∞) : Pk− → P∞− et Kkp → K∞ où K∞ et P∞−
62 Chapitre 6. Le filtre de Kalman

Algorithm 4 Filtre de prédiction


1. Initialisation :

x̂0− = x̄0
P0− = P̄0

2. Pour chaque pas d’échantillonnage k > 0 faire :



x̂k+1 = A x̂k− + Buk + Kkp (yk −Cx̂k− )
rk = yk − ŷk
−1
Kkp = APk−CT CPk−CT +V



Pk+1 = APk− AT − KkpCPk− AT +W

3. Fin Pour

sont deux matrices constantes.

Démonstration. Une preuve de cette condition peut être trouvée dans [3]. 

Si les conditions ci-dessus sont remplies, alors l’établissement du filtre stationnaire est possible.
Avant de déterminer le gain du filtre, il faut résoudre l’équation algébrique de Riccati suivante :

P∞− = (A − K∞pC)P∞− (A − K∞pC)T + K∞pV K∞p T +W (6.13)

Une fois la matrice de covariance P∞− est calculée, nous passons au calcul du gain constant du filtre.
−1
K∞p = AP∞−CT CP∞−CT +V

(6.14)

L’algorithme 5 récapitule les étapes d’application du filtre de Kalman stationnaire.

Algorithm 5 Filtre de Kalman stationnaire


1. Initialisation :

x̂0− = x̄0
P0− = P̄0

2. Pour chaque pas d’échantillonnage k > 0 faire :



x̂k+1 = A x̂k− + Buk + K ∞ (yk −Cx̂k− )
rk = yk −Cx̂k−

3. Fin Pour

Dans des conditions normales de fonctionnement (sans défaut), le résidu réduit rk suit une
distribution gaussienne rk ∼ N (rk , 0,CP∞−C> +V ). L’apparition d’un défaut modifie la moyenne,
la matrice de covariance ou rend le résidu non-blanc (coloré) ou non-gaussien. Ce changement
peut être détecté en utilisant des algorithmes de détection de changement basés sur l’analyse des
propriétés statistiques du signal résidu [2].
6.6 Exercices d’application 63

6.6 Exercices d’application


Exercice 6.1 Considérons le système illustré à la F IGURE 6.7 . La représentation d’état du système
s’écrit
(
xk+1 = Axk + Buk + wk
(6.15)
yk = Cxk + vk

où xk ∈ ℜnx est le vecteur d’état, uk est la commande. L’état initial x0 , les perturbations wk et
le bruit de mesure vk sont des processus aléatoires blancs et Gaussiens non corrélés donnés par,
x0 ∼ N (0, P0 ), wk ∼ N (0,W ) et vk ∼ N (0,V ) respectivement, où P0 ,W et V sont des matrices
définies positives.

F IGURE 6.7 – Estimation sous perte de paquets

Un défaut affectant le réseau reliant les capteurs et l’estimateur a causé des pertes d’informations
(mesures) permanentes. Ce phénomène peut être modélisé par un processus de Bernoulli θk qui
peut prendre la valeur 1 ou 0 avec Pr(θk = 1) = p, où la valeur "0" correspond à une perte de
paquet et "1" à la bonne réception du paquet. L’objectif du filtre du filtre conçu est de trouver la
meilleur estimation d’état x̂k|k à partir des mesures valables pour le filtre. L’estimation optimale
d’état est donnée par

x̂k| j = E [xk |yl , ∀l ≤ j telque θl = 1]

1. Modifier les équations du filtre de Kalman pour traiter ce contexte particulier. Noter que la
matrice de variance Pk est une variable aléatoire, donc, on peut s’intéresser à ces propriétés
statistiques et notamment son espérance mathématique E(Pk )
2. Nous cherchons à démontrer qu’il existe
— pc ∈]0, 1] tel que ∀p ≥ pc : ∃P0 ≥ 0 telque E(Pk ) soit non bornée.
— ∀p < pc et ∀P0 ≥ 0 E(Pk ) est uniformément borné (borne supérieur qui converge).
3. Pour simplifier l’étude on considère le cas scalaire (unidimentionel) régit par l’équation de
Riccati modifiée suivante

A2C2 Pk2
Pk+1 = A2 Pk +W − θ k
C2 Pk +V
4. Pour simplifier l’étude on considère le cas scalaire (unidimentionel) régit par l’équation de
Riccati modifiée suivante
A2C2 Pk2
Pk+1 = A2 Pk +W − θ k
C2 Pk +V
64 Chapitre 6. Le filtre de Kalman

a Donner une borne supérieure de E(Pk ) en modifiant l’équation de Riccati précédente et


en utilisant l’inégalité de Jensen.
b Trouver une condition (sur pc ) de convergence pour cette borne supérieure en utilisant
le théorème du point fixe de Banach (Pour P∞  0).
c Calculer pc pour le système suivant A = 2,C = 1,W = 1,V = 1.

F IGURE 6.8 – Le Rhin Supérieur

Exercice 6.2 Simulation MATLAB/SIMULINK


I- Estimation de la pollution au mercure avec le filtre de Kalman
L’objectif de cette partie est de montrer la capacité du filtre de Kalman dans l’estimation
de l’état d’un système stochastique. Le système considéré est une section d’une rivière, par
exemple 300Km du Rhin Supérieur (voir F IGURE 6.8), dans laquelle des déchets industriels
contenant du mercure sont renversés de temps en temps par un pollueur situé à Bâle (Suisse).
Imaginons que l’Université de Fribourg introduit un appareil de mesure à Breisach, au milieu du
Rhin supérieur, afin de détecter l’éventuelle pollution. Notre objectif consiste à utiliser cette
mesure unique pour estimer la concentration du mercure dans l’ensemble du Rhin supérieur de
Bâle à Bingen et notamment prédire la concentration du mercure à Strasbourg. Dans notre étude
nous considérons les données suivantes :
— La vitesse de l’eau : 8 Km/h
— La longueur totale du Rhin supérieur (de Bâle à Bingen) : 360 km.
— La longueur de la rivière de Bâle à Breisach : 64 km.
— La longueur de la rivière de Breisach à Strasbourg : 80 km.
Pour obtenir un modèle de la rivière, nous discrétisons le Rhin supérieur en plusieurs secteur et
on s’intéresse à l’évolution de l’état chaque 1 heure, c-à-d :
— Période d’échantillonnage : 1 h.
— La plupart de l’eau parcourt un secteur entier chaque heure. Mais 9D’autre part, 1secteur
suivant (plus bas).
— Nombre total de secteurs : 45.
6.6 Exercices d’application 65

— Secteur 1 : Bâle.
— Secteur 9 : Breisach.
— Secteur 19 : Strasbourg.
Dans notre problème, la variable d’état représente la quantité du mercure contenue dans chaque
secteur en [Kg], ç-à-d, nous avons 45 variables d’état X1 , X2 , . . . , X45 . L’entrée de ce système est
Kg
le taux du mercure jeté par le pollueur de Bâle . La seule mesure y est la concentration de
h
mercure à Breisach, qui est mise à l’échelle de telle sorte qu’elle correspond à la quantité totale
(en kg) contenue dans le volume d’eau du 9ème secteur.
Consignes :
1. Écrire le modèle du système sous la forme suivante :
(
x(k + 1) = Ax(k) + Bu(k)
y(k) = Cx(k)

   
   
   
   
   
   
A=

;B = 
 


   
   
   
   

et
 
C=

2. Simuler système pendant 72 heures (3 jours) avec un grand déversement de mercure au


niveau de Bâle à 5 heures du matin de la première journée. Pour ce faire, écrire une boucle
qui vous permet de simuler l’évolution des variables d’état.Vous pouvez introduire une
petite pause et tracer les variations des états pour obtenir une animation. Par exemple, on
peut écrire une boucle comme suit :
nx=45;
nu=1;
ny=1;
N=72;
u=zeros(nu,N);
x=zeros(nx,N);
y=zeros(ny,N);
u(1,5)=300;
x(:,1)=zeros(nx,1);
figure(1)
for i=1:N-1
plot(x(:,i));
pause(0.1);
x(:,i+1)=A*x(:,i) + B * u(:,i);
y(:,i)=C*x(:,i);
end
66 Chapitre 6. Le filtre de Kalman

3. Pour comparer votre simulation précédente avec les résultats obtenus par en utilisant
la commande lsim, écrire un script MATLAB définissant les matrice A, B, C, D et
le système à temps discret d’une période d’échantillonnage T=3600s en utilisant la
commande sys (A, B, C, D, 3600). Vous pouvez comparer votre propre simulation
avec le résultat de lsim comme suit :
figure(2);
subplot(2,1,1);
plot(y); subplot(2,1,2);
ylsim=lsim(sys,u);
plot(ylsim)
4. Supposons maintenant, que nous ne connaissons que les mesures y tandis que les états x
sont inconnus. En outre, comme dans la réalité, nous supposons que nous ne connaissons
pas l’entrée u. Pour définir un filtre de Kalman pour ce système, nous devons d’abord
penser à la manière par la quelle les perturbations et le bruit agissent sur le système.
Nous considérons l’entrée polluante u comme une perturbation du premier état (ç-a-d la
concentration de mercure à Bâle) et supposons que cette perturbation a un grande variance,
beaucoup plus grande que les autres perturbations. Cela nous permet de supposer l’entrée
u est nulle, donc nous pouvons l’abandonner complètement dans les équations, et nous
modélisons le système comme suit :
(
x(k + 1) = Ax(k) + w(k)
y(k) = Cx(k) + v(k)

où w(k) et v(k) représentent les perturbations et le bruit de mesure avec les matrice de
covariance diagonales W et V respectivement. Supposons que la variance de bruit de
mesure est très faible, par exemple [1kg2 ], et que les perturbations d’état pour tous les
secteurs en dehors de Bâle sont également faibles et de la même variance. Définir les
matrices W et V ,de dimensions adéquates, dans l’espace ci-dessous :
 
 
 
 
   
 
W =

 ;V =

 
 
 
 

5. Nous voulons maintenant implémenter les équations du filtre de Kalman

x̂k− = A x̂k−1
Pk− = APk−1 AT +W
−1
Kk = Pk−CT CPk−CT +V


x̂k = x̂k− + Kk (yk −Cx̂k− )


Pk = (In − KkC)Pk−

Programmer le filtre de Kalman et exécuter le, en utilisant les sorties précédemment


simulées ycomme les seules données d’entrée du filtre.
6.6 Exercices d’application 67

6. Tracer les états estimés x̂k sur la même figure avec les véritables états xk .
7. Les éléments diagonaux donnent la variance de l’erreur d’estimation de chaque variable
d’état.
8. Prenant la racine carrée, ajouter les marges d’erreur dans le tracé des états estimés.
9. Regarder l’animation et voir comment le filtre de Kalman détecte la pollution au mercure
et comment les barres d’erreur changent en fonction du temps et de l’espace.

R Pour l’initialisation du filtre de Kalman nous prenons x̂0 = 0 et P0− = In .

II- Détection d’une éventuelle fuite de mercure par le filtre de Kalman Supposons mainte-
nant que l’usine située à Bâle a subit un défaut f dans son système de traitement. Ce défaut a
causé une fuite importante de mercure dans le Rhin Supérieur.
1. Compléter le modèle de système en remplissant la matrice ci-dessous :

(
x(k + 1) = Ax(k) + w(k) + F fk
y(k) = Cx(k) + v(k)

 
 
 
 
 
 
F =



 
 
 
 

2. Peut-on détecter l’apparition de ce défaut à l’aide du filtre de Kalman.


3. Proposer un algorithme qui permet de générer un indicateur de défaut (résidu) rk .
4. Simuler le fonctionnement du filtre en choisissant un défaut constant avec une amplitude
égale à 5 fois la variance de la perturbation du première état.
5. Comment distinguer entre la pollution causée par le déversement quotidien du mercure et
la pollution causée par un défaut dans le système ? Utiliser une technique d’évaluation de
résidu adéquate.


Références
[1] Sarkka, Simo. Bayesian filtering and smoothing. Cambridge University Press, 2013.

[2] M. Basseville and I. Nikiforov. Detection of abrupt changes : Theory and application. Informa-
tion and System sciences, Prentice-Hall, 1993.

[3] T. Kailath, A. H. Sayed, and B. Hassibi. Linear estimation, volume 1. Prentice Hall New Jersey,
2000.
A. Appendice

A.1 La somme de deux Gaussiennes


Pour deux variables aléatoires indépendantes X ∼ pX et Y ∼ pY , la fonction de densité de probabilité
pZ de Z = X +Y est égale au produit de convolution de pX et pY , ç-à-d :

Z ∞
pZ (z) = pY (z − x)pX (x)dx (A.1)
−∞

Sachant que pX et pY sont deux fonctions de densité de probabilité Gaussiennes,

(x−µX ) 2
1 − 2
pX (x) = √ e 2σX (A.2)
2πσX

(y−µY ) 2
1 − 2
pY (y) = √ e 2σY (A.3)
2πσY

En remplaçant dans la convolution :

(z−x−µY ) 2 2
(x−µX )
1 1
Z ∞
− −
2σY2 2
pZ (z) = √ e √ e 2σX dx (A.4)
−∞ 2πσY 2πσX

 
 2 2
2
Z ∞
1

(z − (µX + µY ))2

1  x − σX (z−µ2 Y )+σY µX
σX +σY2

= exp − √ exp −  dx
 
−∞
√ q 2 2
2(σX + σY ) 2π √σX2σY 2
 2
2π σX2 + σY2 σX +σY

2 √ 2 2
σ X σY

σX +σY
(A.5)
70 Chapitre A. Appendice
 
 2 2
2
1

(z − (µX + µY ))2
Z ∞
1  x − σX (z−µ2 Y )+σY µX
σX +σY2

=q exp − √ exp −  dx
 
2(σ 2 + σ 2) √σX σY   2
2 2
2π(σX + σY ) X Y −∞ 2π  
σX2 +σY2 2 √ X2 Y 2
σ σ
σX +σY
(A.6)

L’expression dans l’intégrale est une fonction de densité de probabilité Gaussienne sur x, donc son
intégrale vaut 1. Finalement nous obtenons :

(z − (µX + µY ))2
 
1
pZ (z) = q exp − (A.7)
2π(σ 2 + σ 2 ) 2(σX2 + σY2 )
X Y

⇒ Z ∼ N (z, µs = µX + µY , Ps = PX + PY ) (A.8)

A.2 Dérivée matricielle


Pour deux matrices quelconques M1 et M2 on a :

∂ Trace(M2 M1> ) ∂ Trace(M1 M2> )


= = M2
∂ M1 ∂ M1

∂ Trace(M1> M2 M1 ) ∂ Trace(M1> M2 M1 )
= M2 M1 + M2> M1 ; = 2M2 M1 (si M2 est symétrique)
∂ M1 ∂ M1
Références bibliographiques

Livres
• Chen, J., & Patton, R. J. (2012). Robust model-based fault diagnosis for dynamic systems
(Vol. 3). Springer Science & Business Media.

• Ding, S. (2008). Model-based fault diagnosis techniques : design schemes, algorithms, and
tools. Springer Science & Business Media.

• Gertler, J. (1998). Fault detection and diagnosis in engineering systems. CRC press.

• Isermann, R. (2006). Fault-diagnosis systems : an introduction from fault detection to fault


tolerance. Springer Science & Business Media.

• Blanke, M., Kinnaert, M., Lunze, J., Staroswiecki, M., & Schröder, J. (2006). Diagnosis and
fault-tolerant control.Springer.

• Maquin, D. Ragot J. (2000) Diagnostic des systèmes linéaires, Collection Pédagogique


d’Automatique, Hermès Science Publications, Paris.

• Toscano, R. (2011) Commande et diagnostic des systèmes dynamiques - Modélisation,


analyse, commande par PID et par retour d’état, diagnostic, Ellipses, Paris.

Articles
• Isermann, R. (2005). Model-based fault detection and diagnosis : status and applications.
Annual Reviews in control, 29(1), 71-85.

• Willsky, A. S. (1976). A survey of design methods for failure detection in dynamic systems.
Automatica, 12(6), 601-611.

Vous aimerez peut-être aussi