Main
Main
Main
Polycopié de Cours
Diagnostic et Surveillance des Systèmes
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
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
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
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
défauts
u(t) y(t)
Procédé
Durée Planning
Unité d’usage d’intervention Activités de
maintenance
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
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
entrées sorties
Procédé Industriel
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.
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é.
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
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
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
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
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
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
Méthodes de détection
et de diagnostic des
défauts
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
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)
des opérations de calcul) et spatiale (espace mémoire utilisé) s’avère d’une importance
primordiale dans l’évaluation de ces performances.
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
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
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 :
9. Trouver GBF BF BF BF
yd , Gy fA , Gy fS1 , ∆y .
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.
1 2 3
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
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 :
Générateur de
Capteurs le résidu
Système Résidu
physique C V k
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
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
L’interprétation géométrique de ces relations est appelée résidu directionnel. Ce type de résidu
C3⊥
C2⊥
C1⊥
𝑟𝑘
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 à
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
C̃ = C et ỹ = y
Donc
(
r1 = y1 + y2 − 2y3
r2 = 4y2 − 2y4
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)
Proposition 2.4.1 L’évolution du système dans un horizon temporel de taille s peut être résumée
comme suit :
où
D 0 ··· 0
CB D ··· 0
Hu,s = (2.17)
.. .. .. ..
. . . .
CAs−1 B CAs−2 B . . . D
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
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
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)
où
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
( j)
V ⊥( j) Ho,s j = 0 (2.31)
( 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é
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
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é
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.
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
Entrées ݔሶ ሺݐሻ = ݔܣሺݐሻ + ݑܤሺݐሻ + ܨ௫ ݂ሺݐሻ + ܦ௫ ݀ሺݐሻ, ݔሺ0ሻ = ݔ Sorties
൜
ݕሺݐሻ = ݔܥሺݐሻ + ܨ௬ ݂ሺݐሻ
Reconstructeur d’état
Modèle de simulation du système nominal
L A
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.
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 :
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
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 ?
ẋ(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
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.
ėx (t) = Mex (t) + (MLy Fy + PFy − EFx ) f (t) − Ly Fy f˙(t) (4.8)
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)
sex (s) = Mex (s) + F f (s) + F 0 s f (s) → ex (s) = (sI − M)−1 (F + F 0 ) f (s) (4.11)
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.
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 ?
ẋ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
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.
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.
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 :
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 :
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.
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
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
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)
ẋ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)
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
f1 f2 f3
r1 1 0 1
r2 0 1 1
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.
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
Filtrage d’é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 :
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 )
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
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).
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
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.
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
m0 = x̄0
P0 = P̄0
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
x̂0 = x̄0
P0 = P̄0
• Corréction :
−1
Kk = Pk−CT CPk−CT +V
x̂k = x̂k− + Kk (yk −Cx̂k− )
Pk = (In − KkC)Pk−
3. Fin Pour
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
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−
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
x̂0 = x̄0
P0 = P̄0
• 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).
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
x̂0− = x̄0
P0− = P̄0
−
Pk+1 = APk− AT − KkpCPk− AT +W
3. Fin Pour
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 :
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)
x̂0− = x̄0
P0− = P̄0
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
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.
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
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
— 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=
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 =
x̂k− = A x̂k−1
Pk− = APk−1 AT +W
−1
Kk = Pk−CT CPk−CT +V
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.
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 =
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
Z ∞
pZ (z) = pY (z − x)pX (x)dx (A.1)
−∞
(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
(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)
∂ 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.
• Blanke, M., Kinnaert, M., Lunze, J., Staroswiecki, M., & Schröder, J. (2006). Diagnosis and
fault-tolerant control.Springer.
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.