Introduction au filtrage de Kalman
Introduction au filtrage de Kalman
Sciences de la Mesure
Positionnement dynamique
Jérôme Verdun
* Photo réalisée par KuoShen Choong, chercheur au département d'ingénierie aéronautique de l'université de Sydney (Australie).
3
Avant-propos
Le ltrage de Kalman est une méthode d'estimation de paramètres dont l'un des principaux atouts
1
est de pouvoir prendre en compte, en temps réel, toute nouvelle observation. Cette méthode est lar-
gement répandue dans toutes les applications qui nécessitent la fusion de données en provenance de
divers capteurs - telle la robotique - et constitue l'outil de base de tous les systèmes de navigation
modernes qui équipent les satellites articiels, les missiles balistiques, les avions, les drones, les navires,
les submersibles habités et autonomes ainsi que les voitures haut de gamme .
Le principe de base du ltrage de Kalman repose sur des idées simples, trop souvent occultées par la
lourdeur du formalisme mathématique. Ainsi, le présent document a-t-il pour objectif premier d'illus-
trer ce principe sur un exemple simple an que le lecteur en acquière une vision intuitive. Puis, les
équations du ltrage de Kalman sont établies et justiées à partir de la méthode des moindres car-
rés selon une approche développée par [Sage]. Parmi les nombreuses voies qui conduisent au ltre de
Kalman, cette approche progressive privilégie l'interprétation physique des équations qui conduit plus
aisément à la mise en ÷uvre du ltrage de Kalman dans les applications. La lecture assidue du présent
document devrait donc vous permettre de comprendre les utilisations du ltre de Kalman notamment
dans le domaine du positionnement dynamique, puis de développer assez rapidement vos propres ap-
plications.
Bonne lecture !
1. Rudolf Emil Kalman, ingénieur hongrois né en 1930 qui a travaillé aux États-Unis et en Suisse.
Table des matières
1 Navigation astronomique 7
1.1 Posons le problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Ajoutons une mesure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 De la dynamique en plus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4 Cas limites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Bibliographie 31
5
Chapitre 1
Navigation astronomique
Cet exemple est librement inspiré de l'ouvrage de [Maybeck].
Figure 1.1 Perdu en pleine mer de nuit, comment retrouver votre route ?
7
8 CHAPITRE 1. NAVIGATION ASTRONOMIQUE
Cette densité présente un maximum en x = z1 ; la probabilité que la valeur de x(t1 ) se situe dans
l'intervalle [z1 − σz1 , z1 + σz1 ] (intervalle à 1σ ) est égale à 68, 3%. Vous pouvez donc admettre que la
meilleure estimation x̂(t1 ) de x(t1 ) compte tenu des informations disponibles à l'instant t1 (une seule
mesure) est donnée par
x̂(t1 ) = z1 , (1.1)
avec un écart-type associé σx (t1 ) donné par
px(t1)|z1 (x)
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05 σz1
0.00
−2 0 2 4 z1 6 8 10 12
x
Figure 1.2 Représentation de la densité de probabilité conditionnelle px(t1 )|z1 (x) en fonction des
valeurs possibles x de la position. L'aire hachurée correspond à la probabilité que la variable aléatoire
x(t1 ) se situe dans l'intervalle [z1 − σz1 , z1 + σz1 ], soit 68, 3% pour une distribution normale.
px(t2)|z2 (x)
0.45
0.40
0.35
0.30
σz2
0.25
0.20
0.15
0.10
0.05
0.00
−2 0 2 4 6 z2 8 10 12 x
Figure 1.3 Représentation de la densité de probabilité conditionnelle px(t2 )|z2 (x) en fonction des
valeurs possibles x de la position. La mesure z2 étant plus précise que z1 , la courbe se resserre autour
du maximum x = z2 (trait continu) par rapport à la courbe précédente (trait pointillé).
La question cruciale à ce niveau est de savoir comment combiner les deux mesures z1 et z2 disponibles à
l'instant t2 pour déterminer une meilleure estimation de x(t2 ) ≃ x(t1 ) que celle obtenue avec seulement
l'une des deux mesures. D'après les résultats théoriques que nous établirons plus loin (cf chapitre 3),
l'estimation optimale de x(t2 ) compte tenu des mesures z1 et z2 disponibles correspond à la simple
moyenne pondérée µ suivante :
σz22 σz12
x̂(t2 ) = µ = z1 + z2 , (1.3)
σz12 + σz22 σz12 + σz22
Si les deux mesures sont de même qualité alors σz1 = σz2 et les deux poids des mesures z1 et z2 sont
égaux à 12 dans l'équation 1.3 :
σz22 σz12 1
= = .
σz12 + σz22 σz12 + σz22 2
10 CHAPITRE 1. NAVIGATION ASTRONOMIQUE
px(t2)|z1,z2 (x)
0.45
0.40
0.35
0.30
σ
0.25
0.20
0.15
0.10
0.05
0.00
−2 0 2 4 6 µ 8 10 12 x
Figure 1.4 Représentation de la densité de probabilité conditionnelle px(t2 )|z1 ,z2 (x) en fonction des
valeurs possibles x de la position. Cette courbe est plus étroite autour du maximum x = µ que ne le
sont celles (trait pointillé) représentant la même distribution avec seulement une mesure disponible
(z1 ou z2 ).
À partir de l'expression de x̂(t2 ), il est possible de donner une autre interprétation de la combinaison
1.3. Il vient très simplement :
σz12
x̂(t2 ) = z1 + (z2 − z1 ),
σz12 + σz22
1.3. DE LA DYNAMIQUE EN PLUS 11
De même, l'évolution de la variance sur l'estimation peut être obtenue en exprimant la variance σx2 (t2 )
en fonction de σx2 (t1 ) à partir des relations 1.4 et 1.6 ; il vient successivement :
1 σz21 + σz22 1
2
= 2 2
=
σ σz1 σz2 K(t2 )σz2
1 1 1
= − ,
K(t2 ) σ 2 σz21
d'où, il résulte :
1 1 1
−1 2
= ,
K(t2 ) σ K(t2 )σz21
soit
σ 2 = (1 − K(t2 ))σz21 .
Sachant que σx2 (t2 ) = σ 2 et σx2 (t1 ) = σz21 , il vient nalement :
σx2 (t2 ) = σx2 (t1 ) − K(t2 )σx2 (t1 ). (1.8)
Puisque 1 − K(t2 ) ≤ 1, on en déduit que σx2 (t2 ) ≤ σx2 (t1 ) ⇒ σx (t2 ) ≤ σx (t1 ), et ce d'autant plus que le
gain K(t2 ) est proche de 1. Il est intéressant de noter à ce stade que toute l'information que nous avons
sur notre position à l'instant t2 est donnée par les valeurs de l'estimation x̂(t2 ) et la variance σx2 (t2 ) ;
elles permettent en eet de xer complètement la densité de probabilité conditionnelle px(t2 )|z1 ,z2 (x)
représentée sur la gure 1.4.
La gure 1.5 montre l'évolution de la densité de probabilité conditionnelle px(t)|z1 ,z2 (x) d'après le
modèle d'évolution temporelle 1.9, depuis l'instant t2 de la dernière mesure jusqu'à l'instant t− 3 qui
précède l'acquisition d'une troisième mesure à l'instant t3 . L'estimation de la position est conforme
à un déplacement du navire à la vitesse constante u ; la relation entre les estimations x̂(t−
3 ) et x̂(t2 )
s'écrit donc :
x̂(t−
3 ) = x̂(t2 ) + u(t3 − t2 ). (1.10)
Dans le même temps, l'incertitude sur cette estimation croît de par la présence du bruit sur le modèle
et puisqu'aucune nouvelle information autre que ce modèle d'évolution ne vient la contraindre. La
3 ) et σx (t2 ) s'écrit alors :
relation qui lie les variances σx2 (t− 2
σx2 (t− 2 2 2
3 ) = σx (t2 ) + σw (t3 − t2 ) . (1.11)
px(t)|z1,z2 (x)
0.40
0.35
0.30
0.25
0.20
t2
0.15
t > t2
σ(t2) t−
3 > t
0.10
0.05 σ(t)
σ(t−
3)
0.00
−10 0 10 20 30 40 50 60 70 80 x
x̂(t2) x̂(t)
x̂(t−
3)
Figure 1.5 Représentation de l'évolution de la densité de probabilité conditionnelle px(t)|z1 ,z2 (x) en
fonction des valeurs possibles x de la position entre les instants t2 et t−
3 . La courbe s'élargit progressi-
vement autour du maximum x = x̂(t) ce qui indique que l'incertitude sur l'estimation de la précision
s'accroît de par la présence du bruit de modèle w sur la vitesse u.
1.4. CAS LIMITES 13
Ainsi, avant la mesure prévue à l'instant t3 , l'estimation de la position est-elle x̂(t3 ) avec l'écart-type
σx (t3 ).
Soit z3 la nouvelle mesure disponible à l'instant t3 . En appliquant rigoureusement les relations 1.5, 1.8
et 1.7, l'estimation de la position à l'instant t3 peut s'écrire :
x̂(t3 ) = x̂(t− −
3 ) + K(t3 )(z3 − x̂(t3 )), (1.12)
L'équation 1.12 rend compte de la modication à apporter à l'estimation x̂(t− 3 ) donnée par le modèle
d'évolution lors de l'apparition d'une nouvelle mesure ; elle permet donc d'apporter une correction
3 ) à l'estimation prédite σx (t3 ). L'eet de l'apport d'une nouvelle mesure est donc une
K(t3 )σx2 (t− 2 −
diminution de la variance proportionnelle à σx2 (t− 3 ), égale à K(t3 )σx (t3 ). Le rôle fondamental du gain de
2 −
Kalman apparaît donc clairement puisque ce dernier intervient à la fois sur l'amplitude de la correction
et celle de la diminution de la variance (et donc l'augmentation de la précision) de l'estimation.
Supposons tout d'abord que la mesure z3 soit de mauvaise qualité par rapport à l'estimation dis-
ponible x̂(t−
3 ), ce qui se traduit par σz3 ≫ σx (t3 ). La relation 1.14 indique alors que K(t3 ) ≪ 1. Dans
−
La nouvelle mesure n'a alors pas été prise en compte dans l'estimation x̂(t3 ) qui reste égale à celle
prédite x̂(t−
3 ). Le modèle d'évolution est ici prépondérant pour la détermination de l'estimation x̂(t3 ).
Supposons à présent que ce soit le modèle d'évolution vitesse constante qui soit à son tour peu
able, ce qui se traduit par σw2 ≫ σx (t2 ). Dans ce cas, σx2 (t−
3 ) ≫ 1 (relation 1.10) et K(t3 ) ≫ 1. À la
limite, il vient
σw2 → ∞ ⇒ σx2 (t−
3 ) → ∞ ⇒ K(t3 ) → 1.
Le mesure z3 est donc ici prépondérante puisque c'est elle qui xe la valeur de l'estimation x̂(t3 ) avec
une incertitude nulle σx (t3 ) = 0 (équation 1.13). Le modèle d'évolution est donc ici complètement
occulté.
Enn, si une estimation de la position devient quasiment certaine juste avant la mesure z3 , autre-
ment dit σx (t−
3 ) = 0, alors le gain de Kalman est nul K(t3 ) = 0 et plus aucune mesure supplémentaire
ne pourra changer la valeur de l'estimation x̂(t3 ) = x̂(t−
3 ).
Bien que nous n'ayons encore rien démontré, l'examen de ce petit exemple permet d'illustrer les aspects
principaux du ltrage de Kalman : une méthode d'estimation qui permet d'estimer un paramètre qui
14 CHAPITRE 1. NAVIGATION ASTRONOMIQUE
varie dans le temps en combinant de façon optimale des mesures disponibles et un modèle d'évolution
temporelle tous deux aectés de leurs propres incertitudes. Les équations 1.10 à 1.14 respectent la
structure des phases dites de prédiction et de correction de l'algorithme du ltre de Kalman qui sont
présentes dans toutes les formulations de ce ltre.
Fin du chapitre 1
Chapitre 2
Petit retour sur les moindres carrés
Les développements mathématiques présentés dans la suite sont tirés de [Rogers].
An de comprendre la notion de combinaison optimale évoquée dans l'exemple du chapitre 1, il nous
faut revenir sur la méthode d'estimation par moindres carrés. Dans ce qui suit, nous supposerons que
les paramètres estimés par a méthode ne varient pas dans le temps.
z = Hx + v, (2.1)
où H est une matrice de dimension p × n et v ∈ Rp le vecteur des résidus dont les composantes sont
les diérences entre les observations z et leur modèle linéaire Hx.
Nous supposerons, en outre, que les mesures sont redondantes, c'est-à-dire p ≫ n et que la matrice H
est de rang égal à n.
La meilleure estimation x̂ du vecteur des paramètres x au sens des moindres carrés est celle qui, si elle
existe, minimise la fonction coût JLS dénie par :
JLS (x) = (z − Hx)T (z − Hx) = ∥z − Hx∥22 , (2.2)
où (z − Hx)T désigne la transposée de la matrice z − Hx et ∥ . . . ∥2 la norme 2 (norme euclidienne)
sur Rp .
Le taux de variation de cette fonction coût entre les points x et x + ∆x s'exprime par :
JLS (x + ∆x) − JLS (x) = (z − H(x + ∆x))T (z − H(x + ∆x)) − (z − Hx)T (z − Hx)
= (z − Hx)T − (H∆x)T ((z − Hx) − (H∆x)) − (z − Hx)T (z − Hx)
15
16 CHAPITRE 2. PETIT RETOUR SUR LES MOINDRES CARRÉS
Les deuxième et troisième terme de cette relation correspondent donc respectivement aux diérentielles
première et seconde de la fonction coût. Les extremums de la fonction coût satisfont donc l'équation
matricielle (z − Hx)T H = 0TRn où 0Rn désigne le vecteur nul de Rn . En prenant la transposée membre
à membre de cette dernière égalité, on obtient le système normal d'équations :
H T Hx = H T z. (2.4)
Si le rang de H est égal à n, alors la matrice carrée de dimension n, H T H , est inversible. Le système
normal d'équations admet alors une solution unique x̂, extremum de la fonction coût, qui s'exprime
par :
x̂ = (H T H)−1 H T z. (2.5)
De plus, la matrice H H est également symétrique, dénie et positive. Ainsi, pour tout ∆x ∈ Rn ,
T
∆xT H T H∆x ≥ 0 et ∆xT H T H∆x = 0 uniquement lorsque ∆x = 0Rn . On en déduit que x̂ correspond
au minimum de la fonction coût et constitue l'estimation du vecteur des paramètres x par moindres
carrés.
Exercice 1 (Théorème fondamental) : Montrer que si H est une matrice de dimension p × n de rang
n alors H T H est une matrice inversible de dimension n, symétrique, dénie et positive.
La diagonale de cette matrice contient toutes les variances des paramètres estimés σx̂2i , i = 1, 2, . . . , n.
Les termes non diagonaux correspondent aux covariances entre les paramètres estimés cov(x̂i , x̂j ) pour
(i, j) ∈ {1, 2, . . . , n}2 , i ̸= j .
Il vient successivement d'après 2.6 :
Σε = E ((H T H)−1 H T v)(v T H(H T H)−1 )
Il vient alors :
Σv̂ = (Ip − H(H T H)−1 H T )E(vv T )(Ip − H(H T H)−1 H T )T ,
d'où nalement,
Σv̂ = (Ip − H(H T H)−1 H T )Σv (Ip − H(H T H)−1 H T )T . (2.9)
Exercice 2 (Solution des moindres carrés pondérés) : Démontrer que l'estimation des moindres carrés
pondérés est donnée par :
x̂ = (H T W H)−1 H T W z, (2.11)
et que la matrice de covariance de l'erreur d'estimation Σε s'exprime par :
Σε = (H T W H)−1 H T W Σv W H(H T W H)−1 .
Dans le cas très fréquent en pratique où la matrice de poids est prise égale à l'inverse de la matrice de
covariance des résidus W = Σ−1v , il vient :
Σε = (H T W H)−1 H T W
| W
−1 T −1 T −1 T T −1
{z } W H(H W H) = |(H W H){z H W H}(H W H) ,
Ip In
d'où il résulte :
Σε = (H T W H)−1 . (2.12)
Exercice 3 (Formule de propagation de la variance) : Soit X un vecteur aléatoire et ΣX sa matrice
de covariance. On suppose que le vecteur Y est lié au vecteur X par la relation matricielle Y = AX
où A est une matrice complètement déterministe.
h i
Exprimer la matrice de covariance ΣY de Y à partir de la relation ΣY = E (Y − E(Y )) (Y − E(Y ))T
en fonction de ΣX . La relation obtenue est appelée formule de propagation de la variance.
Exercice 4 : Démontrer les relations 2.8 et 2.9 en utilisant la formule de propagation de la variance.
Fin du chapitre 2
Chapitre 3
Mise à jour des paramètres
Nous allons à présent étudier le problème de la mise à jour de l'estimation de paramètres par moindres
carrés lorsque de nouvelles observations sont disponibles. Nous envisagerons successivement deux cas,
le premier dans lequel les paramètres seront supposés indépendants du temps (cas statique), le second
dans lequel les paramètres seront susceptibles de varier dans le temps en suivant une loi d'évolution
connue (cas dynamique). Pour alléger les notations, nous noterons P la matrice de covariance de
l'erreur d'estimation Σε et R la matrice de covariance des résidus a priori Σv .
Supposons à présent qu'une nouvelle (m + 1)ième mesure zm+1 soit disponible de sorte qu'il faille
modier la relation z = Hx + v en
z H v
= x+ ,
zm+1 hT vm+1
On se propose de déterminer la nouvelle estimation par moindres carrés x̂m+1 sous la forme :
19
20 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES
En notant temporairement
z H
Z= et H = ,
zm+1 hT
et en utilisant la relation 2.11, il vient :
soit nalement :
x̂m+1 = Pm+1 H T R−1 z + hr−1 zm+1 , (3.3)
A−1 = B −1 + C T D−1 C
Ce sont là les équations fondamentales qui permettent de calculer, à partir d'une estimation établie
pour un ensemble de mesures, la nouvelle estimation lorsqu'une nouvelle mesure est disponible.
Remarque 1 : Le vecteur ∆x dénie par la relation 3.1 est nalement donné par :
Remarque 2 : Le vecteur des paramètres estimés x̂ est appelé vecteur d'état. Les quantités x̂m et
x̂m+1 correspondent donc à deux états qui dièrent par leur nombre de mesures (une de plus pour
x̂m+1 ).
sur la nouvelle mesure par l'intermédiaire de la variance r présente dans le terme (hT Pm h + r)−1 .
Le terme restant hT Pm h correspond à la matrice de covariance du vecteur hT x̂m selon la formule de
propagation de la variance. Ce vecteur fournit une estimation de la nouvelle observation zm+1 à partir
des paramètres x̂m estimés avec les m mesures zi , i = 1, 2, . . . , m ; la matrice hT Pm h correspond donc
à l'incertitude sur cette observation.
Remarque 4 : La matrice de covariance de l'erreur d'estimation Pm voit parfois sa dénomination
raccourcie en covariance de l'erreur . Il ne faut néanmoins moins jamais oublier qu'il s'agit bien d'une
matrice.
Exercice 6 : À l'aide des relations 3.6, retrouver les résultats de l'exemple de la navigation astrono-
mique traitée au chapitre 1.
Exercice 7 : Vous disposez de m mesures zi , i = 1, 2, . . . , m d'une grandeur x réputée constante.
Quelle est la meilleure estimation de x au sens des moindres carrés en supposant que toutes les me-
sures sont décorrélés, de variance respectives σi2 , i = 1, 2, . . . , m ? Comment est modiée cette estima-
tion lorsque l'on dispose d'une mesure zm+1 supplémentaire, décorrélée des autres et de variances σm+1 2
?
La démarche présentée ici peut aisément se généraliser au cas où tout un vecteur de nouvelles mesures
est disponible.
Exercice 9 : On suppose à présent qu'un nouvel ensemble de mesures est disponible représenté par
le vecteur zs de dimension s. La fonction coût de ce problème d'estimation peut alors s'écrire sous la
forme matricielle suivante :
T −1
z H R 0m,s z H
JRW LS (x) = − x − x ,
zs Hs 0Tm,s Rs−1 zs Hs
où Hs est une matrice de dimension s × n et Rs la matrice de covariance des résidus a priori sur les s
nouvelles mesures.
On appelle respectivement x̂ et x̂s les estimations obtenues avant et après l'ajout des nouvelles mesures
et P et Ps les matrices de covariance de l'erreur associées à ces estimations. En utillisant la même
démarche que celle employée dans le paragraphe 3.1.1, démontrer les relations suivantes :
x̂s = x̂ + ks (zs − Hs x̂), (3.7)
avec ks = P HsT (Hs P HsT + Rs )−1 ,
et Ps = (In − ks Hs )P.
À partir des relations 3.6, en multipliant par hT à gauche l'équation donnant x̂m+1 , il vient :
hT k →1
hT x̂m+1 = hT x̂m + hT km (zm+1 − hT x̂m ) = hT x̂m + hT km zm+1 − hT km hT x̂m −→
m
zm+1 ,
Ce résultat indique que la mesure zm+1 est obtenue sans erreur par hT x̂m+1 avec une variance hT Pm+1 h
nulle. L'estimation x̂m+1 est en quelque sorte contrainte de restituer sans erreur la mesure réputée de
très grande qualité.
et
Pm+1 → (In − r−1 Pm hhT )Pm ≈ In Pm = Pm .
Ce résultat indique qu'une nouvelle mesure zm+1 de mauvaise qualité n'est que peu prise en compte
dans l'estimation x̂m+1 puisque à la limite, cette dernière égale x̂m , estimation obtenue sans la mesure
zm+1 , avec la même incertitude (Pm+1 = Pm ).
Ainsi, les deux cas limites de prépondérance de la mesure et du modèle d'évolution illustrés dans
le chapitre 1, trouvent-ils là leur justication. Le réglage d'un ltre de Kalman consiste à mettre en
balance les contributions de la mesure et du modèle d'évolution à un instant donné en jouant sur
les niveaux de bruit accordés respectivement aux mesures (bruit de mesure ) et au modèle d'évolution
(bruit de modèle ).
Exercice 10 : Retrouver les résultats de l'étude des cas limites de l'exemple du navire (cf chapitre 1)
à partir des résultats généraux précédemment établis.
Un modèle d'évolution temporelle de x(t) est alors constitué par un système d'équations diérentielles
impliquant les paramètres à estimer.
24 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES
Dans l'exemple du chapitre 1, la trajectoire du navire est supposée parcourue à vitesse constante
ce qui a été traduit par le modèle d'évolution :
dx
= u,
dt
où u représente la vitesse du navire.
Il s'agit là d'une équation diérentielle linéaire du premier ordre, qui est ici complètement déter-
ministe ; en eet, la connaissance d'une condition initiale x(t0 ) = x0 où t0 est l'instant où la position
du navire est x0 permet de déterminer, de façon exacte, l'évolution de la position du navire à un instant
ultérieur ; en eet, la postion à l'instant t est donnée par :
x(t) = x0 + u(t − t0 ).
L'incertitude sur la pertinence de ce modèle a été prise en compte en introduisant le bruit w sous la
forme :
dx
= u + w,
dt
ce qui rend alors le modèle d'évolution non plus déterministe, mais stochastique.
Dans la suite, nous supposerons que le vecteur des paramètres x(t) à l'instant t satisfait le système
diérentiel du premier ordre suivant :
dx
= −Ax(t) + Bu(t) + Cw(t), (3.8)
dt
où A est une matrice carrée de dimension n, B , C deux matrices de dimensions respectives n × q et
n × l, u(t), un fonction vectorielle déterministe de R → Rq appelée commande, et w(t) un vecteur de
bruit blanc à l composantes.
Remarque 5 : Le choix d'un système diérentiel linéaire du premier ordre n'est pas réducteur
puisque :
tout système non linéaire peut se ramener à un système linéaire au voisinage d'une solution
approchée ;
toute équation diérentielle linéaire d'ordre n peut se ramener à la résolution d'un système
diérentiel de n équations à n inconnues.
Exercice 11 : Vérier la première assertion de la remarque précédente en donnant le système dié-
rentiel correspondant à l'équation diérentielle suivante :
d2⃗r
= −k⃗r,
dt2
où k est une constante positive et ⃗r le rayon vecteur de coordonnées cartésiennes ⃗r = [x(t), y(t)]T .
dx dy
On utilisera le vecteur d'état X = [x(t), ẋ(t), y(t), ẏ(t)]T où ẋ(t) = et ẏ(t) = .
dt dt
Exercice 12 : Reprendre l'exercice précédent avec l'équation diérentielle suivante :
d2⃗r
= −k∥⃗r∥2⃗r.
dt2
3.2. CAS DYNAMIQUE 25
Remarque 6 : Le bruit blanc w(t) est un processus aléatoire qui satisfait les deux propriétés suivantes :
∀t ∈ R, E[w(t)] = 0 ;
∀(t, t ) ∈ R , E[w(t)wT (t′ )] = Q(t)δ(t − t′ ),
′ 2
Remarque 7 : Dans l'expression e−(t−t0 )A x(t0 ), le terme e−(t−t0 )A représente l'exponentielle de matrice
dénie par :
1 1
e−(t−t0 )A = In + (t − t0 )A + ((t − t0 )A)2 + . . . + ((t − t0 )A)i + . . . , i ∈ N.
2! i!
Il s'agit donc d'une matrice carrée de même dimension que A. Il conviendra de ne pas inverser e−(t−t0 )A
et x(t0 ) dans le produit e−(t−t0 )A x(t0 ).
Exercice 13 : Dans le cas où les matrices A, B et C qui interviennent dans la relation 3.8 dépendent
explicitement du temps, montrer que la solution du système diérentiel peut s'écrire :
Rt Z t Rt
Z t Rt
−
x(t) = e t0 A(τ )dτ
x(t0 ) + e− τ A(τ ′ )dτ ′
B(τ )u(τ )dτ + e− τ A(τ ′ )dτ ′
C(τ )w(τ )dτ. (3.10)
t0 t0
An de préciser la meilleure estimation de x(t) selon ce modèle d'évolution temporelle, il convient de
poser l'hypothèse suivante :
(3.11)
∀τ ∈ R, E (x(t0 ) − E[x(t0 )])wT (τ ) = 0n,l ,
donc :
∀τ ∈ R, E x(t0 )wT (τ ) = 0n,l .
Cette hypothèse traduit la complète décorrélation entre le bruit w(t) avec la condition initiale x(t0 ) et
donc avec tout les vecteurs de la forme e−(t−t0 )A x(t0 ) puisque pour tout τ ∈ R :
De plus, le bruit w(t) est naturellement décorrélé avec le terme de commande de l'équation 3.9 puisque
pour tout t′ ∈ R :
Z t Z t
−(t−τ )A ′
E e Bu(τ )dτ T
w (t ) = e−(t−τ )A Bu(τ ) E(wT (t′ )) dτ = 0n,1 .
t0 t0 | {z }
0T
l,1
Enn, l'étude de la corrélation avec le terme de bruit de l'équation 3.9 conduit, pour tout t′ ∈ R, à :
Z t Z t
−(t−τ )A ′
E e Cw(τ )dτ T
w (t ) = e−(t−τ )A CE(w(τ )wT (t′ ))dτ
t0 t
Z 0t
= e−(t−τ )A CQ(t′ )δ(t′ − τ )dτ.
t0
An de garantir également la décorrélation du bruit de modèle w(t) et le bruit de mesure représenté
par le vecteur v(t) dans l'équation d'observation z(t) = Hx(t) + v(t) où z(t) est le vecteur des m
observations disponibles à l'instant t et H une matrice de dimension m × n, la condition suivante est
également imposée :
∀(t, t′ ) ∈ R2 , E w(t)v T (t′ ) = 0. (3.12)
À ce stade, on peut utiliser l'équation d'évolution 3.9 pour déterminer le vecteur des paramètres à
l'instant tk+1 , x(tk+1 ), en fonction de fonction de celui obtenu à l'instant tk , x(tk ) où tk et tk+1 , k ∈ N,
sont deux instants consécutifs de la suite croissante d'instants (tn )n∈N . Nous verrons que cette suite
correspond aux instants de disponibilité des mesures. Il vient alors :
Z tk+1 Z tk+1
−(tk+1 −tk )A −(tk+1 −τ )A
x(tk+1 ) = e x(tk ) + e Bu(τ )dτ + e−(tk+1 −τ )A Cw(τ )dτ.
tk tk
Posons xk+1 = x(tk+1 ) et xk = x(tk ) vecteurs d'état respectivement aux instants tk et tk+1 . On dénit
également la matrice de tansition Φk:k+1 entre les états k et k + 1 par la relation suivante :
Φk:k+1 = e−(tk+1 −tk )A . (3.13)
Enn, étant donné la proximité des instants tk et tk+1 en pratique, il est légitime de réaliser les
approximations ci-après :
Z tk+1
e−(tk+1 −τ )A Bu(τ )dτ ≈ (tk+1 − tk )Bu(tk+1 ) = Bk+1 uk+1 ;
tk
Z tk+1
e−(tk+1 −τ )A Cw(τ )dτ ≈ (tk+1 − tk )Cw(tk+1 ) = Ck+1 wk+1 ,
tk
où l'on a posé :
Bk+1 = (tk+1 − tk )B et Ck+1 = (tk+1 − tk )C, (3.14)
puis,
uk+1 = u(tk+1 ) et wk+1 = w(tk+1 ). (3.15)
En tenant compte des équations 3.13, 3.14 et 3.15, l'équation d'évolution peut s'écrire sous la forme :
xk+1 = Φk:k+1 xk + Bk+1 uk+1 + Ck+1 wk+1 . (3.16)
Cette dernière équation correspond à une forme discrétisée de l'équation d'évolution 3.9. Cette mo-
délisation discrète est absolument nécessaire pour tenir compte du fait que les mesures ne sont pas
continûment disponibles dans le temps, mais seulement à des instants discrets.
3.2. CAS DYNAMIQUE 27
En appliquant l'opérateur linéaire espérance uniquement sur les termes aléatoires, il vient :
x̂k+1 = E (Φk:k+1 xk + Bk+1 uk+1 + Ck+1 wk+1 ) = Φk:k+1 E(xk ) + Bk+1 uk+1 + Ck+1 E(wk+1 ),
et puisque E(xk ) = x̂k et E(wk+1 ) = 0, il vient nalement :
x̂k+1 = Φk:k+1 x̂k + Bk+1 uk+1 . (3.17)
Déterminons à présent la matrice de covariance de l'erreur d'estimation Pk+1 dénie par :
Pk+1 = E (xk+1 − x̂k+1 )(xk+1 − x̂k+1 )T .
Avec les équations 3.16 et 3.17, l'erreur d'estimation dans l'état k + 1, εk+1 est donnée par :
εk+1 = xk+1 − x̂k+1 = Φk:k+1 xk + Bk+1 uk+1 + Ck+1 wk+1 − (Φk:k+1 x̂k + Bk+1 uk+1 )
= Φk:k+1 (xk − x̂k ) + Ck+1 wk+1 = Φk:k+1 εk + Ck+1 wk+1 ,
où l'on a posé εk = xk − x̂k .
La détermination de la matrice Pk+1 peut donc se poursuivre selon les calculs suivants :
Pk+1 = E(εk+1 εTk+1 ) = E (Φk:k+1 εk + Ck+1 wk+1 )(Φk:k+1 εk + Ck+1 wk+1 )T
Puisque le bruit w est décorrélé de x(t) à tout instant, il vient : E(xk wk+1
T
) = 0n,l où 0n,l est la matrice
nulle de dimension n × l ; on en déduit nalement que :
T
) = 0n,l et E wk+1 εTk = 0Tn,l .
E(εk wk+1
Finalement, l'équation de propagation de la variance entre les états k et k + 1 est donnée par :
Pk+1 = Φk:k+1 Pk ΦTk:k+1 + Ck+1 Qk+1 Ck+1
T
. (3.18)
Les équations 3.17 et 3.18 constituent les équations de propagation de l'estimation des paramètres au
cours du temps, sachant que ces derniers obéissent à un modèle stochastique d'évolution temporelle.
Exercice 14 : En utilisant les équations 3.17 et 3.18, démontrer les relations ci-après utilisées dans
l'exemple du navire au chapitre 1 :
x̂(t−
3 ) = x̂(t2 ) + u(t3 − t2 ),
avec
σx2 (t− 2 2 2
3 ) = σx (t2 ) + σw (t3 − t2 ) .
28 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES
x̂−
k+1 = Φk:k+1 x̂k + Bk+1 uk+1 ,
−
Pk+1 = Φk:k+1 Pk ΦTk:k+1 + Ck+1 Qk+1 Ck+1
T
.
de matrice de covariance Rk+1 ; ce dernier est supposé lié au vecteur des paramètres xk+1 par :
zk+1 = Hk+1 xk+1 + vk ,
où Hk+1 est une matrice de dimension m × n et vk le vecteur des m résidus a priori à l'instant tk+1 .
−
Pk+1 = (In − Kk+1 Hk+1 )Pk+1 .
Supposons que l'on cherche, à chaque instant k , l'estimateur x̂k de sorte à minimiser l'espérance
E eTk Sk ek où Sk est une matrice symétrique, dénie, positive. Remarquons que si Sk est une matrice
diagonale, l'expression de E eTk Sk ek s'écrit :
E eTk Sk ek = Sk (1)E e2k (1) + Sk (2)E e2k (2) + . . . + Sk (n)E e2k (n) ,
3.3. RÉCAPITULATION DES ÉQUATIONS DU FILTRE DE KALMAN 29
où Sk (1), Sk (2), . . . , Sk (n) sont les éléments de la diagonale de Sk et ek (1), ek (2), . . . , ek (n) sont les
composantes de ek .
comme un algorithme qui blanchit les mesures et donc en extrait le maximum d'informations.
L'analyse statistique des innovations permet de vérier la conformité de la loi de probabilité empirique
de l'innovation avec un bruit blanc centré avec la bonne matrice de variance-covariance. Cette analyse
permet en particulier de vérier le bon fonctionnement d'un ltre de Kalman.
3.3.1 Formulation
Hypothèses fondamentales
Les bruits de modèle wk et de mesure vk sont des bruits blancs décorrélés entre eux :
Initialisation
Le premier vecteur d'état est xé par la condition initiale : x̂0 = E(x(t0 )) ; la matrice de covariance de
l'erreur d'estimation est P0 .
Prédiction
x̂−
k+1 = Φk:k+1 x̂k + Bk+1 uk+1 ,
−
Pk+1 = Φk:k+1 Pk ΦTk:k+1 + Ck+1 Qk+1 Ck+1
T
.
30 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES
Correction
x̂k+1 = x̂− −
k+1 + Kk+1 zk+1 − Hk+1 x̂k+1 ,
− −
Kk+1 = Pk+1 T
Hk+1 (Hk+1 Pk+1 T
Hk+1 + Rk+1 )−1 ,
−
Pk+1 = (In − Kk+1 Hk+1 )Pk+1 .
bien que très simple à programmer, ne garantit pas que la matrice obtenue soit symétrique ou dénie,
positive. Aussi Peter Joseph a proposé dans les années 60 une formulation dite stabilisée qui s'exprime
par [Simon] :
−
Pk+1 = (In − Kk+1 Hk+1 )Pk+1 (In − Kk+1 Hk+1 )T + Kk+1 Rk+1 Kk+1
T
. (3.20)
La formulation donnée par l'équation (3.20) est celle à privilégier dans toute programmation de l'al-
gorithme du ltre de Kalman.
Fin du chapitre 3
Bibliographie
[Sage] ,
Sage, A. P. et Melsa, J. L. Estimation Theory with Applications to Communications and
Control, Éd. MacGraw-Hill, New York, 1971.
[Maybeck] Maybeck Peter S., Stochastic Models, Estimation, and Control, Mathematics in
Science and Engineering, Volume 141-1, Éd. Richard Bellman, Université de Californie Sud,
1979 ; republié par Navtech Book and Software Store, 1994.
Premier chapitre disponible à l'adresse : [Link]
[Rogers] Rogers Robert M., Applied Mathematics in Integrated Navigation Systems, Éd. AIAA
Education Series, Joseph A. Schetz, Reston, Virginie, 2003.
[Simon] ,
Simon Dan Optimal State Estimation , Éd. John Wiley & Sons, 2006.
**********
31