Modélisation et Identification des systèmes (Master) 2018/2019
Chapitre III
Estimation paramétrique des systèmes dynamiques
par la méthode de l’erreur de prédiction
3.1 Introduction
Ce chapitre est réservé à l’estimation paramétrique des systèmes dynamiques en utilisant la méthode de
l’erreur de prédiction et les techniques des moindres carrés tout en se basant sur la connaissance de
plusieurs mesures (couples entrée-sortie) issues de ce système. Cette méthode d’estimation consiste en la
minimisation d’un critère quadratique portant sur l’erreur de prédiction (erreur entre la sortie du système
dynamique et la sortie prédite du modèle ajustable) en vue d’aboutir aux meilleurs paramètres estimés du
modèle mathématique à élaborer.
On considère dans ce chapitre les systèmes dynamiques non commandés opérant dans un environnement
déterministe, qui peuvent être décrits par la classe des modèles mathématiques de type entrée-sortie,
linéaires, monovariable, déterministes (ou légèrement bruités) et à paramètres inconnus, mais constants.
Le problème consiste donc en le développement d’une structure algorithmique permettant d’estimer au
mieux les paramètres d’un système dynamique pouvant être décrit par un modèle mathématique
appartenant à cette classe de modèles mathématiques.
La formulation du problème d’estimation paramétrique posé sera menée en utilisant la méthode des
moindres carrés ordinaires (Gauss 1795) qui est la base de la formulation de la plupart des méthodes
d’estimation paramétrique et structurelle, sui sont présentées dans la littérature.
3.2 Hypothèses retenues et formulation du problème
Ce paragraphe est consacré à la présentation des hypothèses retenues et à la formulation du problème
d’estimation paramétrique des systèmes dynamiques déterministes ou légèrement bruités.
Considérons un système dynamique décrit par un modèle mathématique de type entrée-sortie DARMA
légèrement bruité, tel que :
𝐴(𝑞 −1 )𝑦(𝑘) = 𝐵(𝑞 −1 )𝑢(𝑘) + 𝑒(𝑘) (3.1)
Où 𝑢(𝑘) et 𝑦(𝑘) représentent respectivement l’entrée et la sortie du système à l’instant discret 𝑘,
𝑒(𝑘) désigne le bruit qui agit sur cette sortie , et 𝐴(𝑞 −1 ) et 𝐵(𝑞 −1 ) sont des polynômes donnés par :
𝐴(𝑞 −1 ) = 1 + 𝑎1 𝑞 −1 + ⋯ + 𝑎𝑛 𝑞 −𝑛 (3.2)
𝐵(𝑞 −1 ) = 𝑏1 𝑞 −1 + ⋯ + 𝑏𝑚 𝑞 −𝑚 (3.3)
Où 𝑎𝑖 et 𝑏𝑗 , 𝑖 = 1 … 𝑛 , 𝑗 = 1 … 𝑚 (𝑛 > 𝑚) sont des paramètres inconnus mais supposés constants.
Remarque :
Le niveau du bruit 𝑒(𝑘) agissant sur le système est supposé faible (i.e., sa variance étant assez faible) de
façon qu’il n’ait pas beaucoup d’influence sur la qualité d’estimation et que l’utilisation d’un algorithme
d’estimation déterministe permet de donner des résultats satisfaisants.
3.2.1 Hypothèses retenues
Nous retenons les hypothèses suivantes qui sont généralement admises lors de la formulation d’un
problème d’estimation paramétrique d’un système dynamique pouvant être décrit par un modèle
mathématique entrée-sortie :
15
Modélisation et Identification des systèmes (Master) 2018/2019
1.Le système est en boucle ouverte,
2.Le système est stable,
3.Le système est observable,
4.La classe et la structure du modèle mathématique considéré sont bien connues. Ce modèle
mathématique est supposé ‘’équivalent ‘’au système,
5. L’entrée appliquée au système est supposée stationnaire et suffisamment ‘’riche ‘’,
6. Le nombre de mesure expérimentales 𝑀 (entrée-sortie) issues du système à identifier est
suffisamment grand ,
7. Le bruit agissant sur le système est supposé stationnaire, c’est-à-dire son comportement
dynamique est invariant dans le temps dans tout changement de l’origine du temps. En outre, ce
bruit est supposé ergodique, c’est-à-dire que ses éléments caractéristiques, notamment sa moyenne
et sa variance, peuvent être déterminés à partir d’un seul essai expérimental, en considérant un
intervalle de temps suffisamment grand. Nous supposons que { 𝑒(𝑘), 𝑘 = 1 … 𝑀} est une séquence
de variables aléatoires indépendantes de moyenne nulle et de variance 𝜎 2 .
3.2.2 Formulation du problème
La sortie 𝑦(𝑘) du système dynamique considéré peut s’écrire sous la forme suivante :
𝑦(𝑘) = −𝑎1 𝑦(𝑘 − 1) − ⋯ − 𝑎𝑛 𝑦(𝑘 − 𝑛) + 𝑏1 𝑢(𝑘 − 1) + ⋯ 𝑏𝑚 𝑢(𝑘 − 𝑚) + 𝑒(𝑘) (3.4)
Ou de manière équivalente, en optant l’écriture matricielle
𝑦(𝑘) = 𝜃 𝑇 𝜓(𝑘) + 𝑒(𝑘) (3.5)
Où 𝜃 𝑇 = [𝑎1 … 𝑎𝑛 𝑏1 … 𝑏𝑚 ] (3.6)
𝜓𝑇 (𝑘) = [−𝑦(𝑘 − 1) … − 𝑦(𝑘 − 𝑛) 𝑢(𝑘 − 1) … 𝑢(𝑘 − 𝑚)] (3.7)
Nous définissons la sortie prédite a priori 𝑦̂(𝑘) du système considéré, qui correspond à la sortie du
modèle ajustable, par l’expression suivante :
𝑦̂(𝑘) = −𝑎̂1 (𝑘 − 1)𝑦(𝑘 − 1) − ⋯ − 𝑎̂𝑛 (𝑘 − 1)𝑦(𝑘 − 𝑛) + 𝑏̂1 (𝑘 − 1)𝑢(𝑘 − 1) + ⋯ 𝑏̂𝑚 (𝑘 − 1)𝑢(𝑘 −
𝑚) (3.8)
Où 𝑎̂1 (𝑘 − 1) désigne l’estimé du paramètre 𝑎1 à l’instant discret 𝑘 − 1, etc.
Remarque :
La sortie prédite 𝑦̂(𝑘) est calculée à l’instant 𝑘, tout en se basant sur les valeurs des paramètres estimés à
l’instant discret 𝑘 − 1 , et à partir de la connaissance des informations accessibles à cet instant et aux
instants antérieurs (i.e. 𝑦(𝑘 − 1), 𝑢(𝑘 − 1), 𝑒𝑡𝑐. )
Nous pouvons récrire la sortie prédite a priori 𝑦̂(𝑘) sous la forme compacte suivante :
𝑦̂(𝑘) = 𝜃̂ 𝑇 (𝑘 − 1)𝜓(𝑘) (3.9)
Où 𝜃̂ 𝑇 (𝑘 − 1) = [𝑎̂1 (𝑘 − 1) … 𝑎̂𝑛 (𝑘 − 1) 𝑏̂1 (𝑘 − 1) … 𝑏̂𝑚 (𝑘 − 1)] (3.10)
Dans un contexte idéal, les mesures effectuées sur un système dynamique sont non bruitées et les
hypothèses de travail sont supposées parfaites (modèle mathématique et système étant équivalents). Dans
ce cas, nous pouvons avoir 𝑦(𝑘) = 𝑦̂(𝑘). Cependant, dans un contexte réel (cas d’une application
pratique), ces mesures sont toujours entachées de perturbations aléatoires. Il existe donc un certain écart,
appelé erreur de prédiction, entre la sortie du système et la sortie prédite du modèle mathématique
ajustable.
L’erreur de prédiction ϵ(𝑘) est définie a priori par :
𝜖(𝑘) = 𝑦(𝑘) − 𝑦̂(𝑘) (3.11)
16
Modélisation et Identification des systèmes (Master) 2018/2019
qui peut se récrire comme suit :
𝜖(𝑘) = 𝑦(𝑘) − 𝜃̂ 𝑇 (𝑘 − 1)𝜓(𝑘) (3.12)
A partir des expressions (3.5) et (3.9), nous pouvons récrire l’erreur de prédiction a priori 𝜖(𝑘) sous la
forme matricielle suivante :
𝜖(𝑘) = 𝜓𝑇 (𝑘) [ 𝜃 − 𝜃̂ (𝑘 − 1)] + 𝑒(𝑘) (3.13)
Il est donc bien évident que, si le vecteur 𝜃̂ (𝑘 − 1) présente la meilleure estimation du vecteur des
paramètres 𝜃, alors l’erreur de prédiction 𝜖(𝑘) correspond à la meilleur estimation du bruit 𝑒(𝑘).
Autrement dit, la valeur du bruit 𝑒(𝑘), qui n’est pas mesurable à l’instant discret 𝑘, peut étre estimée à
partir du calcul de la valeur de l’erreur de prédiction 𝜖(𝑘) à cet instant.
Nous pouvons donc écrire les deux expressions suivantes :
1
lim ∑𝑀 ̂
𝑘=1‖𝜃 (𝑘)‖ → ‖𝜃‖ (3.14)
𝑀→∞ 𝑀
1
lim ∑𝑀 2
𝑘=1 𝜖 (𝑘) → 𝜎
2
(3.15)
𝑀→∞ 𝑀
Où 𝜎 2 est la valeur de la variance du bruit 𝑒(𝑘).
La formulation du problème d’estimation des paramètres intervenant dans le vecteur 𝜃, qui est défini par
(3.6) peut être menée en se basant sur la minimisation d’un critère portant sur l’erreur de prédiction 𝜖(𝑘).
Ce critère doit avoir certaines propriétés pour pouvoir aboutir à une solution optimale unique (notamment
un critère quadratique).
Ainsi, nous proposons d’opter pour la minimisation d’un critère portant sur la somme des carrés des écarts
entre les valeurs mesurées de la sortie du système et celle prédites du modèle ajustable, et ce, sur un
horizon de 𝑘 mesures. Ainsi, nous pouvons assurer une estimation paramétrique à minimum de variance
qui se ramène à un problème d’optimisation.
La méthode d’estimation paramétrique qui en résulte est appelée méthode de l’erreur de prédiction, dans
laquelle nous utilisons aussi la méthode des moindres carrés ordinaires.
L’erreur de prédiction 𝜖(𝑘), telle que définie par (3.12), peut s’écrire sous la forme développée
suivante :
𝜖(𝑘) = 𝑦(𝑘) − [−𝑎̂1 (𝑘 − 1)𝑦(𝑘 − 1) − ⋯ − 𝑎̂𝑛 (𝑘 − 1)𝑦(𝑘 − 𝑛)
+𝑏̂1 (𝑘 − 1)𝑢(𝑘 − 1) + ⋯ +𝑏̂𝑚 (𝑘 − 1)𝑢(𝑘 − 𝑚)] (3.16)
Les valeurs des séquences d’entrée {𝑢(𝑘)} et de la sortie {𝑦(𝑘)} sont mesurées aux instants discrets
𝑘 = 1, … , 𝑀. D’après l’expression (3.16) de l’erreur de prédiction 𝜖(𝑘), nous constatons que le calcul de
la séquence {𝜖(𝑘)} se fait à partir de l’instant discret 𝑘 = 𝑛 + 1 jusqu’à 𝑘 = 𝑀, et ceci en raison de la
présence des valeurs non disponibles dans le système d’acquisition et de traitement des données. En fait,
les valeurs 𝑦(0), 𝑢(0), 𝑦(−1), 𝑢(−1), … 𝑦(1 − 𝑛)𝑒𝑡 𝑢(1 − 𝑚) ne sont pas retenues ici.
Un estimateur 𝜃̂(𝑘) optimal du vecteur de paramètres 𝜃 (3.6) a été développé et ce à partir de la
minimisation du critère quadratique 𝐽(𝑘) suivant :
2
𝐽(𝑘) = ∑𝑘𝑖=𝑛+1[𝑦(𝑖) − 𝜃̂ 𝑇 (𝑘)𝜓(𝑖)] (3.17)
Où 𝜃̂(𝑘) est le vecteur des paramètres estimés à l’instant discret 𝑘, qui est défini par :
𝜃̂ 𝑇 (𝑘) = [𝑎̂1 (𝑘) … 𝑎̂𝑛 (𝑘) 𝑏̂1 (𝑘) … 𝑏̂𝑚 (𝑘)] (3.18)
17
Modélisation et Identification des systèmes (Master) 2018/2019
Notons que le critère quadratique 𝐽(𝑘) est une fonction linéaire par rapport aux paramètres estimés
intervenant dans le vecteur 𝜃̂(𝑘). Il en résulte donc que la minimisation du carré de cette erreur de
prédiction peut être faite en utilisant des techniques d’optimisation simples et faciles à mettre en œuvre.
Le terme 𝜃̂ 𝑇 (𝑘)𝜓(𝑖) intervenant dans l’expression du critère 𝐽(𝑘) peut être donné par :
𝜃̂ 𝑇 (𝑘)𝜓(𝑖) = −𝑎̂1 (𝑘)𝑦(𝑖 − 1) − ⋯ − 𝑎̂𝑛 (𝑘)𝑦(𝑖 − 𝑛) + 𝑏̂1 (𝑘)𝑢(𝑖 − 1) + ⋯ 𝑏̂𝑚 (𝑘)𝑢(𝑖 − 𝑚) (3.19)
Nous pouvons mettre l’expression (3.19) sous la forme suivante :
𝜃̂ 𝑇 (𝑘)𝜓(𝑖) = 𝑦̂ (𝑖 | 𝜃̂ (𝑘)) (3.20)
Le terme 𝜃̂ 𝑇 (𝑘)𝜓(𝑖) correspond donc à la sortie prédite à l’instant discret 𝑖, calculée à partir de la
connaissance des paramètres estimés à l’instant discret 𝑘, 𝑖 ≤ 𝑘.
3.3 Méthode d’estimation non récursive des moindres carrés ordinaires
On s’intéresse dans cette partie du cours à la formulation de l’estimateur non récursif des moindres carrés
ordinaires permettant d’estimer les paramètres intervenant dans le modèle mathématique (3.4).
Le critère quadratique 𝐽(𝑘) donné par (3.17) peut être réécrit sous la forme matricielle suivante :
𝐽(𝑘) = 𝐸 𝑇 (𝑘) 𝐸(𝐾) (3.21)
Où 𝐸(𝑘) est le vecteur d’erreur de prédiction tel que :
𝐸(𝑘) = 𝑌(𝑘) − Ψ(𝑘)𝜃̂(𝑘) (3.22)
Dans lequel, nous avons
𝐸 𝑇 (𝑘) = [𝜖(𝑛 + 1) … 𝜖(𝑘)] (3.23)
𝑌 𝑇 (𝑘) = [𝑦(𝑛 + 1) … 𝑦(𝑘)] (3.24)
et
−𝑦(𝑛) −𝑦(𝑛 − 1) … −𝑦(1) 𝑢(𝑛) 𝑢(𝑛 − 1) ⋯ 𝑢(𝑛 + 1 − 𝑚)
−𝑦(𝑛 + 1) −𝑦(𝑛) … −𝑦(2) 𝑢(𝑛 + 1) 𝑢(𝑛) ⋯ 𝑢(𝑛 + 2 − 𝑚)
Ψ(𝑘) = [ ]
⋮ ⋮ … ⋮ ⋮ ⋮ ⋯ ⋮
−𝑦(𝑘 − 1) −𝑦(𝑘 − 2) ⋯ −𝑦(𝑘 − 𝑛) 𝑢(𝑘 − 1) 𝑢(𝑘 − 2) ⋯ 𝑢(𝑘 − 𝑚)
(3.25)
La matrice Ψ(𝑘) peut être écrite sous la forme matricielle suivante :
𝜓𝑇 (𝑛 + 1)
𝑇
Ψ(𝑘) = 𝜓 (𝑛 + 2) (3.26)
⋮
𝑇
[ 𝜓 (𝑘) ]
3.3.1 Estimateur non récursif des moindres carrés ordinaires
L’objectif consiste à déterminer l’estimé optimal 𝜃̂(𝑘) du vecteur des paramètres 𝜃 , tel que donné par
(3.6), et ce en cherchant le minimum critère quadratique 𝐽(𝑘), ainsi, l’estimé optimal de 𝜃̂(𝑘) correspond
à la solution minimale de 𝐽(𝑘). La détermination de cette solution revient à annuler la dérivée du critère
quadratique 𝐽(𝑘) par rapport au vecteur 𝜃̂ (𝑘).
Pour résoudre ce problème d’optimisation, on propose d’exprimer 𝐽(𝑘), donné par (3.21), sous la forme
développée suivante:
18
Modélisation et Identification des systèmes (Master) 2018/2019
𝐽(𝑘) = 𝑌 𝑇 (𝑘)𝑌(𝑘) − 2𝑌 𝑇 (𝑘)Ψ(𝑘)𝜃̂(𝑘) + 𝜃̂ 𝑇 (𝑘)ΨT (𝑘)Ψ(𝑘)𝜃̂(𝑘) (3 .27)
La dérivée de ce critère quadratique 𝐽(𝑘) par rapport au vecteur 𝜃̂(𝑘) donne :
𝜕𝐽(𝑘)
̂ (𝑘) = 2[ΨT (𝑘)Ψ(𝑘)𝜃̂(𝑘) − ΨT (𝑘)𝑌(𝑘)] (3.28)
𝜕𝜃
L’estimé optimal 𝜃̂(𝑘) du vecteur de paramètres 𝜃, tel que donné par (3.6), s’obtient donc en annulant
l’expression (3.28), ce qui donne :
[ΨT (𝑘)Ψ(𝑘)]𝜃̂(𝑘) = ΨT (𝑘)𝑌(𝑘) (3.29)
En supposant que la matrice 𝜓 T (𝑘)𝜓(𝑘) est inversible, on aura :
𝜃̂(𝑘) = [ΨT (𝑘)Ψ(𝑘)]−1 ΨT (𝑘)𝑌(𝑘) (3.30)
Remarques :
Dans le cas où la matrice Ψ(𝑘) est carrée, il est facile de montrer que l’estimateur non récursif des
moindres carrés ordinaires peut se mettre sous la forme suivante:
𝜃̂(𝑘) = Ψ−1 (𝑘)𝑌(𝑘) (3.31)
La condition d’unicité de l’estimateur optimal 𝜃̂(𝑘) peut être vérifiée en calculant la dérivée
seconde du critère quadratique 𝐽(𝑘) par rapport à 𝜃̂(𝑘) tel que :
𝜕2 𝐽(𝑘) T (𝑘)Ψ(𝑘)
̂ (𝑘) = 2Ψ
𝜕𝜃
(3.32)
Ψ T (𝑘)Ψ(𝑘)
est définie positive , alors le critère 𝐽(𝑘) admet in estimé optimal unique 𝜃̂ (𝑘) qui
correspond à un minimum absolu.
3.3.2 Propriétés asymptotiques de l’estimateur non récursif des moindres carrés ordinaires
Etant donné que l’erreur de prédiction 𝜖(𝑘) correspond au meilleur estimé du bruit 𝑒(𝑘), il convient donc
d’étudier les propriétés asymptotiques de l’estimateur non récursif des moindres carrés ordinaires donné
par (3.30). Une telle étude permet de déterminer les conditions qui permettent d’avoir des paramètres
estimés sans biais : c’est l’écart entre les valeurs des paramètres du système 𝜃 et celles des paramètres
estimés 𝜃̂(𝑘). Dans ce cas, on peut dire que l’estimateur non récursif des moindres carrés ordinaires ne
converge pas vers les paramètres du système 𝜃, mais il converge vers des paramètres estimés avec un
certain biais (écart).
Depuis l’équation (3.5), nous pouvons exprimer 𝑒(𝑘) à partir de l’instant 𝑛 + 1 jusqu’à l’instant discret
𝑘 sous la forme matricielle suivante :
𝐸(𝑘) = 𝑌(𝑘) − Ψ(𝑘)𝜃 (3.33)
Où le vecteur de sortie 𝑌(𝑘) et la matrice d’observation Ψ(𝑘) sont définis respectivement par (3.24) et
(3.25) . 𝐸(𝑘) est le vecteur bruit donné par :
𝐸 𝑇 (𝑘) = [𝑒(𝑛 + 1) … 𝑒(𝑘)] (3.34)
Nous pouvons écrire le vecteur de sortie 𝑌(𝑘) comme suit :
𝑌(𝑘) = Ψ(𝑘)𝜃 + 𝐸(𝑘) (3.35)
En tenant compte de (3.35), l’estimateur non récursif des moindres carrés peut être récrit ainsi :
𝜃̂(𝑘) = [ΨT (𝑘)Ψ(𝑘)]−1 ΨT (𝑘)𝑌(𝑘)
= 𝜃 + [ΨT (𝑘)Ψ(𝑘)]−1 ΨT (𝑘)𝐸(𝑘) (3.36)
19
Modélisation et Identification des systèmes (Master) 2018/2019
Nous remarquons que l’expression de 𝜃̂(𝑘), dépend du vecteur 𝐸(𝑘) qui est aléatoire (non mesurable),
d’où nous pouvons considérer l’espérance mathématique (moyenne statique) du vecteur des paramètres
estimés tel que :
𝐸{𝜃̂(𝑘)} = 𝜃 + E{ [ΨT (𝑘)Ψ(𝑘)]−1 ΨT (𝑘)𝐸(𝑘)} (3.37)
Le premier terme de cette expression correspond au vecteur de paramètres du système. Le deuxième terme
représente le biais de l’estimateur non récursif des moindres carrés ordinaires.
Nous définissons ainsi l’erreur d’estimation paramétrique par :
𝜃̃(𝑘) = 𝜃̂(𝑘) − 𝜃 (3.38)
Ou d’une manière équivalente :
𝜃̃(𝑘) = [ΨT (𝑘)Ψ(𝑘)]−1 ΨT (𝑘)𝐸(𝑘) (3.39)
Nous notons que le vecteur du biais d’estimation 𝐵𝑒 (𝑘) est donné par
𝐵𝑒 (𝑘) = 𝐸{𝜃̃(𝑘)} (3.40)
T (𝑘)Ψ(𝑘)]−1
= 𝐸{[Ψ ΨT (𝑘)𝐸(𝑘)} (3.41)
Etant donné que les éléments de la matrices d’observation Ψ(𝑘) sont bien connus à l’instant 𝑘, le calcul
du vecteur du biais d’estimation peut être fait comme suit :
𝐵𝑒 (𝑘) = [ΨT (𝑘)𝜓(𝑘)]−1 ΨT (𝑘)𝐸{𝐸(𝑘)} (3.42)
= [ΨT (𝑘)Ψ(𝑘)]−1 𝜓T (𝑘)𝑚
̅ 𝐸 (3.43)
Où 𝑚
̅ 𝐸 désigne la moyenne statique telle que :
1
̅ 𝐸 = 𝑘−𝑛 ∑𝑘𝑖=𝑛+1 𝐸(𝑖)
𝑚 (3.44)
Remarque :
Depuis les expressions (3.42) et (3 .43), on peut déterminer deux propriétés asymptotiques relatives à
l’estimateur non récursif des moindres carrés ordinaires (3.30)
Estimation biaisée : c’est-à-dire que le terme 𝐸{𝐸(𝑘)} ≠ 0 ou également 𝑚̅ 𝐸 ≠ 0. On distingue
deux cas qui donnent une estimation biaisée, qui sont :
1. Le vecteur 𝐸(𝑘) est de moyenne non nulle , c’est-à-dire 𝐸{𝐸(𝑘)} ≠ 0 (i.e. 𝑚̅ 𝐸 ≠ 0)
2. Le vecteur 𝐸(𝑘) est corrélé avec la matrice d’observations Ψ(𝑘), c’est-à-dire
𝐸{Ψ(𝑘)𝐸(𝑘)} ≠ 0 ou également Ψ(𝑘)𝑚 ̅ 𝐸 ≠ 0.
Estimation non biaisée : c’est-à-dire que le terme 𝐸{𝐸(𝑘)} = 0 ou également 𝑚 ̅ 𝐸 = 0. Dans ce
cas, le vecteur 𝐸(𝑘) est de moyenne non nulle et qu’il non corrélé avec la matrice d’observations
Ψ(𝑘). Ceci nous permet d’écrire les deux expressions suivantes :
lim 𝐸{𝜃̂(𝑘)} → 𝜃 (3.45)
𝑘→∞
et
𝑇
lim 𝐸{[𝜃̂(𝑘) − 𝜃][𝜃̂(𝑘) − 𝜃] } → 0 (3.46)
𝑘→∞
Remarque :
Suite à une analyse détaillée des erreurs d’estimation paramétrique de l’estimateur non récursif des
moindres carrés ordinaires, on peut définir la matrice de variance-covariance de vecteur des paramètres
estimés 𝜃̂(𝑘), comme suit ;
20
Modélisation et Identification des systèmes (Master) 2018/2019
𝑇
𝑣𝑎𝑟[𝜃̂(𝑘)] = 𝐸 {[𝜃̂(𝑘) − 𝜃][𝜃̂(𝑘) − 𝜃] } (3.47)
= 𝜎 2 [ΨT (𝑘)Ψ(𝑘)]−1 (3.48)
𝜎 2 est la variance du bruit 𝑒(𝑘) qui peut être exprimée par l’équation suivante :
1
𝜎̂ 2 (𝑘) = 𝑁−(𝑛+𝑚) [𝑌(𝑘) − Ψ(𝑘)𝜃̂(𝑘)]𝑇 [𝑌(𝑘) − Ψ(𝑘)𝜃̂(𝑘)] (3.48)
Où 𝑁 est le nombre de lignes de la matrice d’observations Ψ(𝑘) tel que 𝑁 > (𝑛 + 𝑚).
3.3.3 Mise en œuvre pratique de l’estimateur non récursif des moindres carrés ordinaires
Nous pouvons souligner que la mise en œuvre pratique de l’estimateur non récursif des moindres carrés
ordinaires présente plusieurs inconvénients. Nous citons ici plus particulièrement :
Détermination non récursive des paramètres estimés 𝜃̂(𝑘), qui est ennuyeux dans une application
en temps réel ;
Nécessité d’une inversion matricielle, qui peut poser des problèmes numériques, notamment la
perte de positivité de la matrice Ψ𝑇 (𝑘) Ψ(𝑘) pendant une certaine étape de calcul ;
Capacité importante de stockage et de traitement des données, qui se traduisent par une
augmentation des dimensions de la matrice d’observation Ψ(𝑘) et du vecteur de sortie 𝑌(𝑘) au
fur et à mesure de l’introduction de nouvelle données (couples entrée-sortie). En outre, ceci
permet d’augmenter le temps de calcul de 𝜃̂ (𝑘) d’une étape à une autre.
3.4 Méthode d’estimation récursive des moindres carrés ordinaires
Ce paragraphe est réservé à l’élaboration d’une version récursive de l’estimateur non récursif des
moindres carrés ordinaires (3.30). Dans ce contexte, nous allons développer une structure algorithmique
permettant d’estimer récursivement les paramètres intervenant dans un vecteur 𝜃, tel que donné par (3.6).
Dans une telle structure d’estimation récursive, les nouveaux paramètres estimés sont calculés à partir de
la connaissance des estimés précédents plus un terme de correction ; ceci peut être illustré comme suit :
‘’Paramètres estimés à l’instant en cours’’=’’Paramètres estimés à l’instant antérieur’’ +’’correction’’
Ceci peut se traduire par l’expression suivante :
𝜃̂(𝑘) = 𝜃̂(𝑘 − 1) + 𝐶(𝑘) (3.49)
Dans laquelle 𝜃̂(𝑘) et 𝜃̂ (𝑘 − 1) représentent les vecteurs des paramètres estimés aux deux instants
discrets 𝑘 et 𝑘 − 1, respectivement, et 𝐶(𝑘) désigne le terme de correction.
3.4.1 Algorithme d’estimation récursif des moindres carrés ordinaires RLS
Pour le développement d’une structure algorithmique permettant d’estimer récursivement les paramètres
intervenant dans le vecteur 𝜃, nous allons reprendre le problème de minimisation du critère quadratique
𝐽(𝑘) donné par (3.17). Le vecteur de paramètres estimés optimal 𝜃̂(𝑘), qui correspond au minimum de
ce critère quadratique, s’obtient en résolvant la condition suivante :
𝜕𝐽(𝑘)
̂ (𝑘) = −2 ∑𝑘𝑖=𝑛+1[𝑦(𝑖) − 𝜃̂ 𝑇 (𝑘)𝜓(𝑖)] 𝜓(𝑖) = 0 (3.50)
𝜕𝜃
Nous pouvons donc écrire :
∑𝑘𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) = [∑𝑘𝑖=𝑛+1 𝜓(𝑖)𝜓(𝑖)𝑇 ]𝜃̂ (𝑘) (3.51)
21
Modélisation et Identification des systèmes (Master) 2018/2019
Ou encore, en supposant que la somme des matrices 𝜓(𝑖)𝜓(𝑖)𝑇 existe,
−1
𝜃̂ (𝑘) = [∑𝑘𝑖=𝑛+1 𝜓(𝑖)𝜓(𝑖)𝑇 ] ∑𝑘𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) (3.52)
Nous avons supposé, jusqu’à maintenant, que l’estimation des paramètres intervenant dans le vecteur
𝜃, se fait en se basant sur des valeurs mesurées (couples entrée-sortie) du système à identifier, et ce, à
partir de l’instant 𝑛 + 1 jusqu’à l’instant discret 𝑘 + 1 ; L’adjonction d’un nouveau couple d’entrée-
sortie permet d’estimer ces paramètres à l’instant discret 𝑘 + 1. Ainsi, en considérant l’équation (3.52),
nous pouvons écrire :
𝜃̂ (𝑘 + 1) = 𝑃(𝑘 + 1) ∑𝑘+1
𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) (3.53)
Où la matrice 𝑃(𝑘 + 1), appelée matrice du gain d’adaptation, est définie par :
−1
𝑃(𝑘 + 1) = [∑𝑘+1 𝑇
𝑖=𝑛+1 𝜓(𝑖)𝜓 (𝑖) ] (3.54)
Nous pouvons écrire l’expression suivante :
∑𝑘+1 𝑘
𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) = ∑𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) + 𝑦(𝑘 + 1)𝜓(𝑘 + 1) (3.55)
qui se réécrit en considérant l’équation (3.51)
∑𝑘+1 𝑘 𝑇 ̂
𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) = [∑𝑖=𝑛+1 𝜓(𝑖)𝜓 (𝑖) ]𝜃 (𝑘) + 𝑦(𝑘 + 1)𝜓(𝑘 + 1) (3.56)
En ajoutant et en retranchant la quantité 𝜓(𝑘 + 1)𝜓𝑇 (𝑘 + 1)𝜃̂(𝑘) dans le membre droit de l’équation
précédant, nous obtenons :
∑𝑘+1 𝑘+1 𝑇 ̂ ̂𝑇
𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) = [∑𝑖=𝑛+1 𝜓(𝑖)𝜓 (𝑖) ]𝜃 (𝑘) + 𝜓(𝑘 + 1)[ 𝑦(𝑘 + 1) − 𝜃 (𝑘) 𝜓 (𝑘 + 1)] (3.57)
Ainsi, on déduit que :
∑𝑘+1 −1 ̂
𝑖=𝑛+1 𝑦(𝑖)𝜓(𝑖) = 𝑃 (𝑘 + 1) 𝜃 (𝑘) + 𝜓(𝑘 + 1)𝜖(𝑘 + 1) (3.58)
Ainsi, compte tenu des deux expressions (3.53) et (3.58), nous pouvons exprimer le vecteur des
paramètres estimés optimal 𝜃̂(𝑘) sous la forme récursive suivante :
𝜃̂(𝑘 + 1) = 𝜃̂(𝑘) + 𝑃(𝑘 + 1)𝜓(𝑘 + 1)𝜖(𝑘 + 1) (3.59)
De l’expression (3.59), nous remarquons que le calcul du vecteur des paramètres estimés 𝜃̂ (𝑘 + 1) se
base sur la connaissance du vecteur des paramètres estimés 𝜃̂(𝑘) . Cependant, le calcul dépend de la
matrice du gain d’adaptation paramétrique 𝑃(𝑘 + 1) qui se détermine de façon non récursive à partir de
l’expression (3.54). Ainsi, pour avoir un vecteur des paramètres estimés optimal 𝜃̂ (𝑘) sous une forme
totalement récursive, nous devons chercher une forme récursive pour la matrice du gain d’adaptation
paramétrique 𝑃(𝑘 + 1). Notons toutefois que cette matrice 𝑃(𝑘 + 1) est symétrique définie positive.
A partir de l’expression (3.54), nous pouvons écrire :
𝑃−1 (𝑘 + 1) = ∑𝑘𝑖=𝑛+1 𝜓(𝑖)𝜓𝑇 (𝑖) + 𝜓(𝑘 + 1)𝜓𝑇 (𝑘 + 1) (3.60)
Ou encore
𝑃(𝑘 + 1) = [𝑃 −1 (𝑘) + 𝜓(𝑘 + 1)𝜓𝑇 (𝑘 + 1)]−1 (3.61)
Où la matrice 𝑃(𝑘) est définie par :
−1
𝑃(𝑘) = [∑𝑘𝑖=𝑛+1 𝜓(𝑖)𝜓𝑇 (𝑖) ] (3.62)
22
Modélisation et Identification des systèmes (Master) 2018/2019
Lemme 1 :
Soient la matrice de gain d’adaptation 𝑃(𝑘) et le vecteur d’observation 𝜓(𝑘), nous pouvons donc écrire
la propriété d’inversion matricielle suivante :
[𝑃−1 (𝑘) + 𝜓(𝑘 + 1)𝜓 𝑇 (𝑘 + 1)]−1 = 𝑃(𝑘) − 𝑃(𝑘)𝜓(𝑘 + 1)[1 + 𝜓 𝑇 (𝑘 + 1)𝑃(𝑘)𝜓(𝑘 + 1)]−1 𝜓 𝑇 (𝑘 + 1)𝑃(𝑘)
(3.63)
En appliquant l’inversion matricielle (3.63) à l’équation (3.61), nous pouvons exprimer la matrice de gain
d’adaptation paramétrique 𝑃(𝑘) sous la version récursive suivante :
𝑃(𝑘 − 1)𝜓(𝑘)𝜓𝑇 (𝑘)𝑃(𝑘 − 1)
𝑃(𝑘) = 𝑃(𝑘 − 1) − (3.64)
1 + 𝜓 𝑇 (𝑘)𝑃(𝑘 − 1)𝜓(𝑘)
Finalement, nous pouvons définir l’algorithme d’estimation récursif des moindres carrés ordinaires
𝑅𝐿𝑆 (Recursive Least Squares), qui permet d’estimer les paramètres intervenant dans le vecteur 𝜃 des
paramètres, par le système d’équations suivant :
𝜃̂(𝑘) = 𝜃̂(𝑘 − 1) + 𝑃(𝑘)𝜓(𝑘)𝜖(𝑘)
𝑃(𝑘 − 1)𝜓(𝑘)𝜓𝑇 (𝑘)𝑃(𝑘 − 1)
𝑃(𝑘) = 𝑃(𝑘 − 1) − (3.65)
1 + 𝜓𝑇 (𝑘)𝑃(𝑘 − 1)𝜓(𝑘)
𝜖(𝑘) = 𝑦(𝑘) − 𝜃̂ 𝑇 (𝑘 − 1)𝜓(𝑘)
Soulignons que le terme de correction 𝐶(𝑘) intervenant dans l’expression (3.49) est tel que :
𝐶(𝑘) = 𝑃(𝑘)𝜓(𝑘)𝜖(𝑘) ; Il représente donc le vecteur de correction dans l’algorithme d’estimation
récursif des moindres carrés ordinaires RLS (3.65).
Remarque 1 :
L’algorithme d’estimation récursif RLS permet le traitement séquentiel des valeurs mesurées (entrée,
sortie) issues du système à identifier, qui sont disponibles à chaque pas de calcul. Dans ce contexte, ces
valeurs mesurées peuvent être évaluées au fur et à mesure que le système évolue. Ainsi, l’estimation
paramétrique du système considéré, en se basant sur une structure algorithmique récursive, offre plusieurs
avantages par rapport à la structure non récursive tels que :
Capacité réduite de stockage et de traitement des données expérimentales. En fait, le calcul récursif
des paramètres estimés à un instant discret 𝑘 donné se base uniquement sur quelques informations
antérieures, et ce, en fonction de l’ordre du modèle mathématique considéré(e.g., pour un modèle
mathématique d’ordre 2, nous nous limitons uniquement à des informations disponibles aux
instants 𝑘 − 1 et 𝑘 − 2) ;
Temps de calcul relativement faible. Notons toutefois que la valeur de ce temps ne varie pas d’une
itération à une autre ;
Mise en œuvre pratique simple et flexible, notamment dans une structure numérique informatisée.
Ajoutons que cette mise en œuvre pratique ne nécessite pas de techniques particulières ;
Implémentation facile et efficace en temps réel, et plus particulièrement dans un schéma de
commande adaptative. En outre, l’utilisation des systèmes informatiques (cartes d’interfaçage,
DSP, etc.) dans des schémas d’identification ou de commande numérique donne plus de souplesse
quant à l’implémentation de l’algorithme d’estimation récursif RLS.
23
Modélisation et Identification des systèmes (Master) 2018/2019
Remarque 2 :
Comme toute forme récursive, l’implémentation de l’algorithme d’estimation récursif des moindres
carrés ordinaires RLS nécessite son initialisation. Dans ce cas, nous devons :
1. Initialiser le vecteur des paramètres estimés 𝜃̂ (0).
2. Initialiser la matrice de gain d’adaptation paramétrique 𝑃(0). En fait, le choix de 𝑃(0) est lié à la
distance paramétrique entre les paramètres du système à identifier et les paramètres du vecteur
𝜃̂(0). C’est-à-dire, si la distance paramétrique est faible, nous n’avons pas besoin d’un gain
d’adaptation important, et inversement. En général, si on n’a aucune information sur le système,
nous optons pour le choix suivant : 𝑃(0) = 𝛼𝐼 où 𝛼 est un paramètre positif et 𝐼 est une matrice
unité.
Il est très important de souligner que le choix de 𝜃̂(0) et 𝑃(0) ne conditionne pas la convergence de
l’algorithme RLS , par contre, il influe sur sa vitesse de convergence
Pour tester les performances et les efficacités de l’algorithme d’estimation récursif des moindres carrés
ordinaires RLS lors de son application, en simulation numérique, à un système dynamique pouvant être
décrit par un modèle mathématique DARMA légèrement bruité, et pour évaluer la qualité d’estimation
paramétrique obtenue, nous pouvons considérer la distance paramétrique 𝑑(𝑘) définie par :
𝑁 𝑚 2 0,5
2
𝑎𝑖 − 𝑎̂𝑖 (𝑘) 𝑏𝑗 − 𝑏̂𝑗 (𝑘)
𝑑(𝑘) = [∑ [ ] +∑ [ ] ] (3.66)
𝑎𝑖 𝑏𝑖
𝑖=1 𝑗=1
Nous pouvons facilement montrer que l’algorithme d’estimation récursif des moindres carrés
ordinaires RLS (3.65) peut se mettre sous la forme suivante :
𝜃̂(𝑘) = 𝜃̂(𝑘 − 1) + 𝐾(𝑘)𝜖(𝑘)
𝑃(𝑘 − 1)𝜓(𝑘)
𝐾(𝑘) =
1 + 𝜓𝑇 (𝑘)𝑃(𝑘 − 1)𝜓(𝑘)
(3.67)
𝑃(𝑘 − 1)𝜓(𝑘)𝜓𝑇 (𝑘)𝑃(𝑘 − 1)
𝑃(𝑘) = 𝑃(𝑘 − 1) −
1 + 𝜓 𝑇 (𝑘)𝑃(𝑘 − 1)𝜓(𝑘)
𝜖(𝑘) = 𝑦(𝑘) − 𝜃̂ 𝑇 (𝑘 − 1)𝜓(𝑘)
Où 𝐾(𝑘) représente le vecteur du gain d’adaptation paramétrique.
Nous remarquons que le calcul de la sortie prédite a priori 𝑦̂(𝑘), telle que donnée par l’équation
(3.9), se fait à partir de la connaissance des paramètres estimés à l’instant discret 𝑘 − 1 . Toutefois,
il est possible de calculer la sortie prédite a posteriori, 𝑦̂0 (𝑘), et ce en se basant sur les valeurs des
paramètres estimés à l’instant discret 𝑘. Ainsi, nous pouvons définir cette sortie prédite a
posteriori, 𝑦̂0 (𝑘) par l’expression suivante :
𝑦̂0 (𝑘) = 𝜃̂ 𝑇 (𝑘)𝜓(𝑘) (3.68)
Dans laquelle 𝜃̂(𝑘) représente le vecteur des paramètres estimés à l’instant discret 𝑘.
Ainsi, nous définissons l’erreur de prédiction a posteriori 𝜖0 (𝑘), qui représente l’écart entre la
sortie 𝑦(𝑘) du système et celle prédite a posteriori 𝑦̂0 (𝑘) du modèle ajustable, par :
24
Modélisation et Identification des systèmes (Master) 2018/2019
𝜖0 (𝑘) = 𝑦(𝑘) − 𝜃̂ 𝑇 (𝑘)𝜓(𝑘) (3.69)
D’où il est facile à montrer, depuis les expressions (3.67) et (3.69) que :
𝜖(𝑘)
𝜖0 (𝑘) = 𝑇
(3.70)
1 + 𝜓 (𝑘)𝑃(𝑘 − 1)𝜓(𝑘)
Compte tenu de cette expression, nous pouvons affirmer que le comportement asymptotique de l’erreur
de prédiction a priori 𝜖(𝑘) peut s’écrire sous la forme suivante, et ce , en fonction de l’erreur a
posteriori 𝜖0 (𝑘), c’est-à-dire :
lim 𝜖(𝑘) = 𝜖0 (𝑘) (3.71)
𝑘→∞
si les conditions suivantes sont respectées :
1. Le gain d’adaptation paramétrique utilisé dans l’algorithme d’estimation récursif des moindres
carrés ordinaires RLS est décroissant, c’est-à-dire les paramètres de la matrice 𝑃(𝑘) décroissent
asymptotiquement vers zéro.
2. Les composantes de la matrice d’observation 𝜓(𝑘) sont finies
Nous pouvons donc définir l’algorithme d’estimation récursif RLS, utilisant une estimation a posteriori,
par :
𝜃̂(𝑘) = 𝜃(𝑘 − 1) + 𝑃(𝑘)𝜓(𝑘)𝜖0 (𝑘)
𝑃(𝑘 − 1)𝜓(𝑘)𝜓𝑇 (𝑘)𝑃(𝑘 − 1)
𝑃(𝑘) = 𝑃(𝑘 − 1) − (3.72)
1 + 𝜓𝑇 (𝑘)𝑃(𝑘 − 1)𝜓(𝑘)
𝑦(𝑘) − 𝜃̂ 𝑇 (𝑘 − 1)𝜓(𝑘)
𝜖0 (𝑘) =
1 + 𝜓𝑇 (𝑘)𝑃(𝑘 − 1)𝜓(𝑘)
Remarque 3 :
Nous pouvons énoncer les deux remarques suivantes concernant l’algorithme d’estimation récursif des
moindres carrés ordinaires RLS
1. Il n’est pas robuste vis-à-vis d’une certaine variation des paramètres du système considéré.
2. Il n’est pas robuste vis-à-vis d’un certain niveau de bruit 𝑒(𝑘), c’est-à-dire à partir d’une certaine
valeur de la variance 𝜎 2 de ce bruit, la qualité d’estimation est médiocre.
3.4.2 Analyse de convergence de l’algorithme d’estimation récursif RLS
Toute analyse de convergence vise à satisfaire principalement les deux propriétés suivantes :
lim 𝐸{𝜃̂(𝑘)} → 𝜃 (3.73)
𝑘→∞
et
𝑇
lim 𝐸{[𝜃̂(𝑘) − 𝜃][𝜃̂(𝑘) − 𝜃] } → 0 (3.74)
𝑘→∞
25
