Modélisation et Identification des Processus

Ben Hamouda KAMEL

October 4, 2021
Faculté des Sciences de Bizerte
Identification des Modèles
Approche Physique

• Méthode basée sur une description physique du système.

• Identification d’un modèle continu.
• Une série d’expériences ciblées permet de mesurer indépendamment
chaque paramètre du modèle.
• Cette approche peut être envisagée sous certaines conditions :
• Les équations du système sont parfaitement connues.
• Les sous-systèmes peuvent être isolés de manière à mesurer
indépendamment chaque paramètre.
• Parfois, nécessité d’utiliser des instruments de mesure additionnels

1
Exemple : Identification des paramètres d’un MCC

• Identification des paramètres électriques :

• Blocage du rotor
• Application d’un échelon indicielle de tension d’induit
di
u = Ri + L ⇒ U (s) = RI (s) + LsI (s)
dt
U (s) U0
⇒ I (s) = avec U (s) =
R + Ls s
U0  L

⇒ i (t) = 1 − e− R t
R
U0 R
i (∞) = et τ =
R L

2
Exemple : Identification des paramètres d’un MCC

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
0 1 2 3 4 5 6 7 8 9 10

3
Exemple : Identification des paramètres d’un MCC

• Identification de la constante électro-mécanique :

• Entraı̂nement externe par un moteur à une vitesse ω0 connue
• Mesure à l’aide d’un voltmètre de la tension à vide aux bornes de l’induit
u0 (moteur fonctionnant en génératrice)
u0
u0 = Kω0 ⇒ K =
ω0

4
Exemple : Identification des paramètres d’un MCC

• Identification des paramètres mécaniques :

• On active l’asservissement de courant
• On considère que la boucle de courant est infiniment rapide, donc sa
fonction de transfert est égale à 1
• On applique un échelon indicielle de courant d’induit :
dω
Γm − Γr = J
dt
dω
Ki − f ω = J ⇒ KI (s) = (f + Js) Ω (s)
dt
K I0
Ω (s) = I (s) avecI (s) =
f + Js s
KI0  J

ω (t) = 1 − e− f t
f

5
Méthodes Graphiques

• Méthodes simples utilisant la réponse indicielle

• Identification d’un modèle continu
• Limitations :
• Système d’ordre peu élevé
• ne permet pas d’identifier les zéros
• Système peu bruité
• Paramètres identifiés grossièrement (convient pour la commande mais pas
pour la prédiction).

6
Systèmes du premier ordre
K E
G (s) = 1+τ s avec U (s) = s

Réponse Indicielle
1

0.9 X: 4
Y: 0.9502
0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
0 1 2 3 4 5 6 7 8 9 10

7
Systèmes du second ordre sous-amortis

• Forme canonique

Kωn2
G (s) =
s2 + 2ξωn s + ωn2

• Caractéristiques de la réponse indicielle d’un tel système :

• La tangente à l’origine est de pente nulle
• Réponse oscillatoire (Dépassement)
• Pseudo-pulsation

8
Systèmes du second ordre sous-amortis

• D2 ' 0

√−πξ
ξ / D1 = Ke 1−ξ2
1.8
Réponse indicielle système du second ordre sous-amorti et
ωn / t1 = √π 2
1.6

1.4

1.2
ωn 1−ξ
1

0.8

0.6

0.4

0.2

0
0 2 4 6 8 10 12 14 16 18 20

• D2  0
1 D1
2π ln D2
ξ= r
1
1+ D1
4π 2 ln2
D2

et
2π
√
ωn =
(t2 −t1 ) 1−ξ 2 9
Systèmes du second ordre sous-amortis

• D2 ' 0

√−πξ
ξ / D1 = Ke 1−ξ2
1.8
Réponse indicielle système du second ordre sous-amorti et
ωn / t1 = √π 2
1.6

1.4

1.2
ωn 1−ξ
1

0.8

0.6

0.4

0.2

0
0 2 4 6 8 10 12 14 16 18 20

• D2  0
1 D1
2π ln D2
ξ= r
1
1+ D1
4π 2 ln2
D2

et
2π
√
ωn =
(t2 −t1 ) 1−ξ 2 9
Systèmes amortis d’ordre supérieur : Méthode de Strejc

Soit le système amorti d’ordre supérieur

Ke−τ s
G (s) = n
(1 + Ts)

Y(t)

Dy = K Du
Point d’inflexion

0
T1 T2 t 10
Systèmes amortis d’ordre supérieur : Méthode de Strejc

1. On trace la tangente au point d’inflexion I pour déterminer T1 et T2

2. Relever T1 et T2 en déduire l’ordre n en utilisant le tableau ci-joint.
T2
3. Déterminer la constante du temps T à partir de la colonne T .
4. Déterminer le retard τ quand il existe à partir de la différence entre la
valeur de T1 mesurée et celle donnée par la colonne TT1 .

NB : Entre deux lignes du tableau, on choisit la valeur de n la plus faible.

11
Systèmes amortis d’ordre supérieur : Méthode de Strejc

n T1 T2 T1
T T T2
1 0 1 0
2 0.28 2.72 0.1
3 0.8 3.7 0.22
4 1.42 4.46 0.32
5 2.10 5.12 0.41
6 2.81 5.70 0.49

12
Systèmes avec un intégrateur
−τ s
e
Soit le système G (s) = K s(1+Ts)n

- le coefficient d’intégration K * L’ordre est donné par la

- le retard τ courbe suivante :
- la constante du temps T
- l’odre n n

* La réponse indicielle étant la 3

suivante :
2

0.1 0.15 0.2 0.25 0.3 0.35 AB/AC

Y(t) On mesure le rapport : AB/AC

D2 D2 // D1

Asymptote D1

C
B
0
t
A’ A

13
Systèmes avec un intégrateur

A0 A
• Si n est entier, calculer T = n et le temps mort τ est nul
AB
• Si n n’est pas entier , déterminer le nouveau rapport AC correspondant à
la partie entière de n. Pour cela déplacer D2 parallèlement à D1 vers D1
pour obtenir ce nouveau rapport. Le temps mort τ est égale à la
translation effectuée par D2 . Calculer T à partir de A0 A = τ + nT
( AAC0 A )
• Calculer le coefficient d’intégration K = ∆u

14
Simplification de modèle

• Si le système contient plusieurs pôles réels différents, on peut garder le

pôle le plus lent, c’est-à-dire celui qui est le plus proche de 0.
• Si le système contient un ensemble de pôles réels et complexes
conjugués de partie réelle différente, on garde le pôle ou la paire de
pôles complexes conjugués dont la partie réelle est la plus proche de 0.

15
Méthodes Harmoniques

Méthodes basées sur la réponse fréquentielle du système

• Le signal d’entrée est une sinusoı̈de dont la fréquence balaie la plage de

fonctionnement du système
• Pour un système linéaire, la sortie est une sinusoı̈de ayant subit un
déphasage et une modification d’amplitude par un gain
• Le relevé du gain et du déphasage pour plusieurs fréquences permet de
tracer le diagramme de Bode du système puis d’en déduire sa fonction
de transfert en relevant les points de cassure

NB : Méthodes adaptées à l’identification de modèles linéaire continus,

stables, à retard, avec intégrateurs et zéros.

16
Méthodes par corrélation

• Méthode basée sur l’extraction du signal utile en sortie en le corrélant

avec une sinusoı̈de
• Méthode adaptée aux cas où la mesure est très bruitée
• nécessite la numérisation des données et l’utilisation d’un ordinateur
pour leur traitement
• Il existe des dispositifs commerciaux qui implémentent ce type
d’identification

17
Démarche expérimentale

• Utilisation d’un GBF pour • Utilisation d’un ordinateur

générer le signal u et un équipé d’un CnA et d’un CAn
oscilloscope numérique pour pour respectivement générer
acquérir u et y. le signal d’entrée u et
Le modèle échantillonné est le échantillonner le signal de
suivant : sortie y.
Le modèle est le suivant :

18
Méthodes par corrélation

• Utilisation d’un GBF et d’un oscilloscope :

Z {ye (t)} Y (z)

G (z) = Z {G (s)} = =
Z {ue (t)} U (z)

• Utilisation d’un ordinateur :

 
−1

G (s) Y (z)
G (z) = 1 − z Z =
s U (z)


• Réponse harmonique : G z = ejωTe pour 0 < ω < ∞
• Si Te → 0 la réponse harmonique de G (z) tend vers celle de G (s) dans
sa bande passante.
• L’utilisation d’un ordinateur est préférable car G (z) intègre le BOZ.

19
Méthodes par corrélation

Soit u [k] et y [k] les n valeurs échantillonnées de u (t) et y (t) à la période Te

avec u (t) = E cos (ωt).

u [k] = Re EejωkTe


g [l] est la réponse impulsionnelle du système

∞
∞ 
P  jω(k−l)Te
P jω(k−l)Te
y [k] = g [l] Re Ee = Re g [l] Ee
l=0 l=0
∞
 
g [l] e−jωlTe = [Link] ejωkTe G ejωTe


= [Link] ejωkTe
P

l=0


= E G ejωTe cos ωkTe + arg G ejωTe

20
Méthodes par corrélation

Soit les grandeurs suivantes :

n n
1X 1X
Ic = y [k] cos (ωkTe ) Is = y [k] sin (ωkTe )
n n
k=0 k=0

Supposons que la sortie du système en régime permanent est bruitée, la

sortie peut s’écrire comme suit :

y [k] = E G ejωTe cos ωkTe + arg G ejωTe

+ v [k]

 
n
1
E G ejωTe cos ωkTe + arg G ejωTe  cos (ωkTe )
P


Ic =

n
k=1 | {z } | {z }
Gω ϕω
n
+ n1
P
v [k] cos (ωkTe )
k=1
21
Méthodes par corrélation

n
11
P
Ic = n 2 EGω [cos (2ωkTe + ϕω ) + cos (ϕω )]
k=1

n
+ 1n
P
v [k] cos (ωkTe )
k=1

On obtient alors

n
EGω EGω P
Ic = 2 cos (ϕω ) + 2n cos (2ωkTe + ϕω )
k=1

n
+ 1n
P
v [k] cos (ωkTe )
k=1

22
Méthodes par corrélation

n
EGω P
* 2n cos (2ωkTe + ϕω ) → 0 quand n → ∞ car la valeur
k=1
moyenne d’un cosinus et nulle.
n
* 1n
P
v [k] cos (ωkTe ) → 0 quand n → ∞ puisque le bruit est
k=1
supposé non coloré (bruit blanc), ne disposant donc pas de composante
fréquentielle à la pulsation ω.
EGω
* On obtient finalement Iclim = 2 cos (ϕω )
* On démontre aussi de la même façon que Islim = − EG2 ω sin (ϕω )
* d’où
q
2 Ic2lim + Is2lim
Gω = et ϕω = arctan 2 (−Islim , Iclim )
E

23
Prédiction
Objectif

Identification de modèles paramétriques à partir de données expérimentales

24
Définitions

• Système (ou processus) : la partie de l’univers que nous avons décidé

de modéliser et avec laquelle nous interagissons
• Modèle (M) : règles permettant de calculer le plus précisément possible
les grandeurs du système (ses sorties notamment) en fonction de
grandeurs connues ou mesurées (entrées appliquées et mesures passées)
• Entrées (u) : grandeurs connues agissant sur le système dont nous
maı̂trisons l’évolution
• Sorties (y) : grandeurs caractéristiques du système que nous observons
et mesurons
• Perturbations et bruits (e) : ce sont des grandeurs mal connues et non
maı̂trisables qui agissent sur le système ou sur ses sorties
• États (x) : ce sont les grandeurs qui permettent de décrire
complètement l’état du système à un instant donné

25
La procédure d’identification

L’identification est une procédure en 5 étapes

• Acquisition des données expérimentale :
• Utilisation d’un ordinateur équipé d’une carte d’entrées/sorties permet
d’avoir des données échantillonnées
• Conditionnement des données
• Choix de la structure de modèle :
• Plusieurs candidats possible
• Physiquement, il n’existe pas de ”vrai” modèle
• Choix du critère :
• Choisir un critère J à minimiser permettant de comparer les modèles entre
eux
• Calcul du jeu de paramètres θ∗ optimal
• Critique des résultats et validation :
• Soumettre le modèle obtenu à un ensemble de tests
• Décider si le modèle est acceptable. Sinon, recommencer la procédure

26
Rappel

• Définition de la transformée en z :
∞
X
Z {f (t)} = f (kTe ) z−k = F (z)
k=0

• notations :
q−1 : Opérateur de retard d’une période
d’où q−1 f (t) = f (t − Te )
• sans perte de généralité on prend Te = 1
q−1 f (t) = f (t − 1) avec t = kTe = k
ainsi f [k] = f (kTe ) = f (t)

27
Rappel

Soit {h (t)} = {h0 , h1 , h2 , h3 ....} réponse impulsionnelle du système.

Z {h (t)} = H {z} sa transmittance
La réponse y (t) du système à entrée f (t) est :
∞
X
y (t) = h (k) u (t − k)
k=0

Par définition :
∞
X
y (t) = h (k) u (t − k) = H (q) u (t)
k=0

Remarque :
Si H (z) = h0 + h1 z−1 + h2 z−2 + ... + hn z−n
alors Y (z) = H (z) U (z)
et y (t) = H (q) u (t) = h0 u (t) + h1 u (t − 1) + h2 u (t − 2) + ... + hn u (t − n)

28
Rappel

• Soit e une variable aléatoire. La probabilité P que e soit comprise entre

a et b est donnée par :
Zb
P (a ≤ e ≤ b) = fe (x) dx
a
2


• fe (x) = N m, σ est la fonction densité de probabilité FDP de la
variable e.

29
Rappel

• Espérance mathématique ou moyenne

Z
E (e) = xfe (x) dx
R

• Densité de probabilité jointe fy,z (xy , xz )

probabilité d’un événement xy et xz (y et z variables non
indépendantes).

30
Rappel

• Densité de probabilité conditionnelle

fy,z (xy , xz )
fz|y (xz |xy ) =
fy (xy )
•
Zb
P (a ≤ z ≤ b|y = xy ) = fz|y (xz |xy ) dxz
a

• Si y et z sont deux variables indépendantes :

fy,z (xy , xz ) = fy (xy ) fz (xz )

fz|y (xz |xy ) = fz (xz )

31
Prédiction à un pas

• Soit le modèle générique suivant :

y (t) = G (q) u (t) + H (q) e (t) (1)

• e (t) est un bruit blanc de valeur moyenne nulle.
• H (q) est monique tel que h (0) = h0 = 1

32
Prédiction à un pas

• Le bruit filtré v (t) est tel que :

∞
X ∞
X
v (t) = H (q) e (t) = h (k) e (t − k) = e (t)+ h (k) e (t − k) (2)
k=0 k=1
| {z }
m(t−1)

• considérons m (t − 1) la partie du bruit v (t) connue jusqu’à l’instant

(t − 1)
∞
X
m (t − 1) = h (k) e (t − k) (3)
k=1

33
Prédiction à un pas

• Soit fe (x) la fonction densité de probabilité de la variable e (t)

P (x ≤ e (t) ≤ x + dx) = fe (x) dx (4)

• Soit fv (x) la fonction densité de probabilité de la variable v (t)

fv (x) dx = P (x ≤ v (t) ≤ x + dx|v (−∞) , ..., v (t − 1))

= P (x ≤ e (t) + m (t − 1) ≤ x + dx)
= P (x − m (t − 1) ≤ e (t) ≤ x − m (t − 1) + dx) = fe (x − m (t − 1)) dx

fv (x) = fe (x − m (t − 1))
(5)

Remarque : La connaissance des réalisation antérieures

{v (−∞) , ..., v (t − 1)} implique la connaissance de m (t − 1) ce qui permet
le passage de la ligne 1 à la ligne 2 de l’équation précédente.
34
Prédiction à un pas

• e (t) est un bruit blanc, sa FDP est représentée par une loi NORMALE
de moyenne nulle, fe (x) admet un maximum pour x = 0.
fv (x) = fe (x − m (t − 1)) admet un maximum à m (t − 1).
• Soit v̂ (t|t − 1) la valeur la plus probable de v (t) à l’instant t,
connaissant les réalisations antérieures jusqu’à l’instant t − 1

v̂ (t|t − 1) = m (t − 1) (6)
• Remarque : Il s’agit d’une estimation au sens de Maximum A
Posteriori (MAP) ou Maximum Likelihood Estimator (MLE)

35
Prédiction à un pas

∞ ∞
!
X X
−k
v̂ (t|t − 1) = h (k) e (t − k) = h (k) q e (t) = (H (q) − 1) e (t)
k=1 k=1

v (t)
= (H (q) − 1)
H (q)

= 1 − H −1 (q) v (t)


(7)

• La prédiction à un pas devient :

ŷ (t|t − 1) = G (q) u (t) + v̂ (t|t − 1)

= G (q) u (t) + 1 − H −1 (q) v (t)

= G (q) u (t) + 1 − H −1 (q) (y (t) − G (q) u (t))

(8)
36
Prédiction à un pas

• On obtient finalement :

ŷ (t|t − 1) = H −1 (q) G (q) u (t) + 1 − H −1 (q) y (t)

(9)

ou encore

H (q) ŷ (t|t − 1) = G (q) u (t) + (H (q) − 1) y (t)

• L’erreur de prédiction est alors :

y (t) − ŷ (t|t − 1) = H −1 (q) (y (t) − G (q) u (t)) = e (t) (10)

e (t) est appelé l’innovation à l’instant t

37
Principaux modèles
ARX : AutoRegressive eXogenous Model

• Considérons l’équation récursive suivante :

y (t) + a1 y (t − 1) + ... + ana y (t − na )

| {z }
Partie Autoregressive
(11)
= b1 u (t − 1) + ... + bnb u (t − nb ) + e (t)
| {z } |{z}
Partie Exogene Bruit

• Le vecteur des paramètres θ est comme suit :

T
θ = (a1 a2 ... ana b1 b2 ... bnb )

38
ARX : AutoRegressive eXogenous Model

• Considérons

A (q) = 1 + a1 q−1 + a2 q−2 + ... + ana q−na

B (q) = 1 + b1 q−1 + b2 q−2 + ... + bna q−nb

B(q) 1
G (q) = A(q) et H (q) = A(q)

• Le système devient donc :

A (q) y (t) = B (q) u (t) + e (t)

(12)
⇒ y (t) = G (q) u (t) + H (q) e (t) = y (t, θ)

39
ARX : AutoRegressive eXogenous Model

Modèle ARX

40
ARX : AutoRegressive eXogenous Model

• Prédicteur à un pas :
- D’après la relation (9)
ŷ (t|t − 1, θ) = B (q) u (t) + (1 − A (q)) y (t) (13)
• Définissant le régresseur :

T
ϕ (t) = (−y (t − 1) ... − y (t − na ) u (t − 1) ... u (t − nb ))
• L’équation (13) peut se réécrire comme suit :

ŷ (t|t − 1, θ) = ϕT (t) θ (14)

C’est une régression linéaire.

41
ARMAX : AutoRegressive Moving Average with eXogeneous inputs

• La relation récursive d’entrées-sorties étant la suivante :

y (t) + a1 y (t − 1) + ... + ana y (t − na )

| {z }
Autoregressive Part

= b1 u (t − 1) + ... + bnb u (t − nb )
| {z } (15)
Exogeneous Part

+ e (t) + c1 e (t − 1) + ... + cnc u (t − nc )

| {z }
Moving Average

• Le vecteur des paramètres θ de ce modèle est comme suit :

T
θ = (a1 a2 ... ana b1 b2 ... bnb c1 c2 ... cnc )

42
ARMAX : AutoRegressive Moving Average with eXogeneous inputs

• Considérons

A (q) = 1 + a1 q−1 + a2 q−2 + ... + ana q−na

B (q) = 1 + b1 q−1 + b2 q−2 + ... + bna q−nb
C (q) = 1 + c1 q−1 + c2 q−2 + ... + cnc q−nc

B(q) C(q)
G (q) = A(q) et H (q) = A(q)

• Le système devient donc :

A (q) y (t) = B (q) u (t) + C (q) e (t)

(16)
⇒ y (t) = G (q) u (t) + H (q) e (t) = y (t, θ)

43
ARMAX : AutoRegressive Moving Average with eXogeneous inputs

Modèle ARMAX

44
ARMAX : AutoRegressive Moving Average with eXogeneous inputs

• Prédicteur à un pas
• D’après la relation (9) on a :

   
AB A B A
ŷ (t|t − 1, θ) = u (t) + 1 − y (t) = u (t) + 1 − y (t)
CA C C C
(17)

• NB : pour alléger les notations on utilisera le raccourci X au lieu de

X (q) pour un polynôme en q−1 .

45
ARMAX : AutoRegressive Moving Average with eXogeneous inputs

• L’équation de prédiction peut être réécrite comme suit :

Cŷ (t|t − 1, θ) = Bu (t) + (C − A) y (t) (18)

• En isolant la valeur actuelle de la prédiction de ses valeurs antérieures

on obtient :

ŷ (t|t − 1, θ) = c1 q−1 ŷ (t|t − 1, θ) + ... + cnc q−nc ŷ (t|t − 1, θ)

+ Bu (t) + (C − A) y (t) (19)

= (1 − C) ŷ (t|t − 1, θ) + Bu (t) + (C − A) y (t)

46
ARMAX : AutoRegressive Moving Average with eXogeneous inputs

• Or (C − A) = (C − 1) − (A − 1)
• L’équation (19) devient :

ŷ (t|t − 1, θ) = Bu (t) + (1 − A) y (t) + (C − 1) [y (t) − ŷ (t|t − 1, θ)]

(20)
• Soit ε (t) l’erreur de prédiction telle que :

ε (t) = y (t) − ŷ (t|t − 1, θ) (21)

• Le régresseur est tel que :

ϕ (t, θ) = [−y (t − 1) ... − y (t − na ) u (t − 1) ... u (t − nb )

ε (t − 1) ... ε (t − nc ) ]T

Remarque : Le régresseur fait intervenir l’erreur de prédiction, il

dépend alors du vecteur paramètres θ.

47
ARMAX : AutoRegressive Moving Average with eXogeneous inputs

• L’équation (20) devient :

ŷ (t|t − 1, θ) = ϕT (t, θ) θ (22)

• La relation décrite par l’équation (22) est non linéaire, il s’agit d’une
régression pseudo-linéaire.

48
OE : Output Error

• La relation récursive d’entrées-sorties étant la suivante :

y (t) + f1 y (t − 1) + ... + fnf y (t − nf )
| {z }
Autoregressive Part

= b1 u (t − 1) + ... + bnb u (t − nb )
| {z } (23)
Exogeneous Part

+ e (t) + f1 e (t − 1) + ... + fnf e (t − nf )

| {z }
noise

• Le vecteur des paramètres θ de ce modèle est comme suit :

T
θ = b1 b2 ... bnb f1 f2 ... fnf
• On définit F tel que :

F (q) = 1 + f1 q−1 + ... + fnf qnf

49
OE : Output Error

• Ce modèle considère une erreur de sortie

B (q)
G (q) = et H (q) = 1 (24)
F (q)
• Le système devient donc :

y (t) = FB u (t) + e (t)

(25)
⇒ y (t) = G (q) u (t) + e (t) = y (t, θ)

50
OE : Output Error

Modèle OE

51
OE : Output Error

• Prédicteur à un pas
• D’après la relation (9) on a :

ŷ (t|t − 1, θ) = 1 FB u (t) + (1 − 1)y (t)

(26)
= FB u (t)

52
OE : Output Error

• L’équation de prédiction peut être réécrite comme suit :

Fŷ (t|t − 1, θ) = Bu (t) (27)

• En isolant la valeur actuelle de la prédiction de ses valeurs antérieures

on obtient :

ŷ (t|t − 1, θ) = −f1 q−1 ŷ (t|t − 1, θ) − ... − fnf q−nf ŷ (t|t − 1, θ) + Bu (t)

= Bu (t) + (1 − F) ŷ (t|t − 1, θ)
(28)

53
OE : Output Error

• Le régresseur est tel que :

ϕ (t, θ) = [u (t − 1) ...u (t − nb )

ŷ (t − 1|t − 2, θ) ...ŷ (t − nf |t − nf − 1, θ) ]T

• L’équation (28) devient :

ŷ (t|t − 1, θ) = ϕT (t, θ) θ (29)

• C’est une régression pseudo-linéaire.

54
BJ : Box and Jenkins

• On définit

B (q) = 1 + b1 q−1 + ... + bnb q−nb

C (q) = 1 + c1 q−1 + ... + bnc q−nc
D (q) = 1 + d1 q−1 + ... + dnd q−nd
F (q) = 1 + f1 q−1 + ... + fnf q−nf
• L’équation récursive est donnée par la relation suivante :

y (t) = FB u (t) + DC e (t)

B C
G (q) = F et H (q) = D

y (t) = G (q) u (t) + H (q) e (t) = y (t, θ)

• Le vecteur des paramètres θ de ce modèle est comme suit :

T
θ = b1 ...bnb c1 ...cnc d1 ...dnd f1 ...fnf
55
BJ : Box and Jenkins

Modèle BJ

56
BJ : Box and Jenkins

• Prédicteur à un pas
• D’après la relation (9) on a :
 
DB D
ŷ (t|t − 1, θ) = u (t) + 1 − y (t) (30)
CF C
• Soient :
B C
v0 (t) = u (t) et v (t) = e (t) = y (t) − v0 (t)
F D
• On montre que :
ŷ (t|t − 1, θ) = (C − 1) e (t) + Bu (t) + (1 − D) v (t) + (1 − F) v0 (t)
(31)
• Le régresseur est tel que :
ϕ (t, θ) = [u (t − 1) ...u (t − nb ) e (t − 1) ...e (t − nc )

− v (t − 1) ... − v (t − nd ) − v0 (t − 1) ... − v0 (t − nd ) ]T

ŷ (t|t − 1, θ) = ϕT (t, θ) θ
57
Modèle d’état

• Modèle de représentation d’état :

(
x (t + 1) = Ax (t) + Bu (t) + v (t)
y (t) = Cx (t) + w (t)

• v (t) et w (t), bruits blancs à valeur moyenne nulle.

58
Principe de parcimonie

• Principe lié au choix de la complexité du modèle

• Par exemple :

b1 q−1 b1 q−1
M1 (θ1 ) = et M2 (θ2 ) =
1 + a1 q−1 1 + a1 q−1 + a2 q−2

• L’ensemble des comportements que peut générer M1 (.) est inclus dans
l’ensemble des comportements que peut générer M2 (.).
• Pour un jeu de données identique, M2 (.) modélisera le comportement
toujours plus finement que M1 (.).

59
Principe de parcimonie

• Si la structure réelle du processus est M1 (.) mais que les mesures sont
bruitées, M2 (.) modélisera toujours mieux le comportement que M1 (.).
Par contre M2 (.) aura un moins bon comportement prédictif.

• Principe de parcimonie : choisir la structure de modèle la moins

complexe possible pour conserver un bon pouvoir prédictif.

60
Méthodes Quadratiques
Régression linéaire

• La régression linéaire consiste à approcher la réponse temporelle du

système par un modèle linéaire.
• Un modèle linéaire par rapport à ses paramètres a une expression du
type :
x (t) = a1 f1 (t) + a2 f2 (t) + ... + ak fk (t) (32)
où fi (t) sont les formants du modèle.
• Exemple : Approximation polynomiale

x (t) = a0 + a1 t + ... + ar tr (33)

61
Formulation du problème

• Le modèle linéaire définit par l’équation (32) peut se mettre sous une
forme matricielle :

x (ti ) = ϕT (ti ) θ
(34)
T T
ϕ (ti ) = [f1 (ti ) f2 (ti ) ... fk (ti )] et θ = [a1 a2 ... ak ]
• ainsi :

X = Φθ

T  T
X = [x (t1 ) x (t2 ) ... x (tn )] et Φ = ϕT (t1 ) ϕT (t2 ) ... ϕT (tn )
(35)
• notons :
• t1 , t2 , ..., tn les instants des n mesures.
• yi = y (ti ) les sorties mesurées aux instants ti .
• xi = x (ti ) les expressions du modèle aux instants de mesure.
62
Formulation du problème

• Soit k la dimension de θ. Si n = k, alors Φ est une matrice carrée.

θ = Φ−1 X
• Cependant, si les mesures sont bruitées, il peut être intéressant d’avoir
des mesures supplémentaires pour assurer une redondance de
l’information.
• En introduisant des mesures supplémentaires, la matrice Φ n’est plus
carrée, et l’équation Y = Φθ n’admet pas de solution dans Im (Φ).

Im (Φ) : Image du sous-espace engendré par les vecteurs de la matrice Φ.

63
Formulation du problème

• Soit εi l’écart entre la mesure de la sortie (yi ) et l’expression du modèle

(xi ) à l’instant ti .

εi = ε (ti ) = y (ti ) − ϕT (ti ) θ

Y = X + ε = Φθ + ε (36)

T
ε = [ε1 ε2 ...εn ]
Les ε (ti ) sont souvent appelés résidus.
• En vue d’éviter la compensation des erreurs entre elles, on considère le
critère J tel que :
n n
X X
2
J = εT ε = ε2i = y (ti ) − ϕT (ti ) θ (37)
i=1 i=1

64
Formulation du problème

• L’objectif est de trouver un estimé θ̂ de θ, à partir des n mesures

permettant la minimisation du critère de l’erreur quadratique J.

65
Recherche du minimum de J

• On suppose qu’il y a indépendance des paramètres, le critère de l’erreur

quadratique J admet un minimum pour :

∂J ∂2J
=0 et >0 i = 1, ..., k (38)
∂ai ∂a2i
• Le minimum est obtenu pour un jeu de valeurs dites optimales pour le
critère considéré (critère des moindres-carrés).

∂J Resolution des k equations

=0 ⇒ â1 , â2 , ..., âk (39)
∂ai
• Le modèle le plus proche des mesures au sens du critère J est :

x̂ (t) = â1 f1 (t) + â2 f2 (t) + ... + âk fk (t) = ϕT (t) θ̂ (40)

66
Recherche du minimum de J

• On montre que :

−1
θ̂ = ΦT Φ ΦT Y (41)

minimise le critère de l’erreur quadratique J.

−1
• Le terme ΦT Φ ΦT est appelé pseudo-inverse de Φ.

67
Recherche du minimum de J : Démonstration (Méthode I)

• Reprenons l’équation (37)

n n
X X
2
J = εT ε = ε2i = y (ti ) − ϕT (ti ) θ
i=1 i=1

∂J
• θ̂ est solution ∂θ = 0, ainsi :
n  
∂J
ϕ (ti ) y (ti ) − ϕT (ti ) θ̂ = 0
P
∂θ = −2
i=1

n n
ϕ (ti ) ϕT (ti ) θ̂
P P
ϕ (ti )y (ti ) =
i=1 i=1

 n
−1 n
ϕ (ti ) ϕT (ti )
P P
θ̂ = ϕ (ti )y (ti )
i=1 i=1

68
Recherche du minimum de J : Démonstration (Méthode I)

• Or
 
n
X h i y (t1 )
T
ϕ (ti )y (ti ) = ϕ (t1 ) ... ϕ (tn )  ...  = Φ Y
 
i=1 y (tn )

• et
 
n
X h i ϕ (t1 )
ϕ (ti ) ϕT (ti ) = ϕ (t1 ) ... ϕ (tn ) T
 ...  = Φ Φ
 
i=1 ϕ (tn )

• il vient donc :

−1
θ̂ = ΦT Φ ΦT Y

69
Rappel : Dérivation matricielle

T
• Soit X = [x1 x2 ... xn ] et f (X) = f (x1 , x2 , ... , xn ), une fonction à
valeurs réelles.

−−→
• Notons gradf (X) le gradient de la fonction f (X) tel que :
 ∂f (X) 
∂x1
−−→ →
− ∂f (X) 
gradf (X) = ∇f (X) = = ...

∂X

∂f (X)
∂xn

70
Rappel : Dérivation matricielle

• Soit f (X) = AT X
∂AT X
=A
∂X
• Soit f (X) = X T X
∂X T X
= 2X
∂X
• Soit A une matrice carrée symétrique

∂X T AX
= 2AX
∂X
• Soit A(n,p) et Y(p,1)
∂X T AY
= AY
∂X

71
Recherche du minimum de J : Démonstration (Méthode II)

• En se référant aux équations (36) et (37)

T
J = εT ε = (Y − Φθ) (Y − Φθ)

T T
= Y T Y − Y T Φθ − (Φθ) Y + (Φθ) Φθ

= Y T Y − 2θT ΦT Y + θT ΦT Φθ

∂J
• La valeur approchée θ̂ est telle que ∂θ =0

∂J
∂θ = −2ΦT Y + 2ΦT Φθ̂ = 0

−1
θ̂ = ΦT Φ ΦT Y

72
Pondération du critère J

• On définit la matrice de pondération Q telle que :

 
q1 . . . 0
 . .
. . ... 

Q=  .. 
0 · · · qn

• La matrice de pondération Q est définie positive (toutes ses valeurs

propres sont strictement positive)
• Q est symétrique
• Le critère de l’erreur quadratique devient :

J = εT Qε (42)

73
Pondération du critère J

• Minimisation du critère pondéré

h i
∂J ∂ T
θ̂ / ∂θ = ∂θ (Y − Φθ) Q (Y − Φθ) = 0

   T
= −ΦT Q Y − Φθ̂ − Y − Φθ̂ QΦ = 0
 
= −2ΦT Q Y − Φθ̂ = 0 (43)

−1
⇒ θ̂ = ΦT QΦ ΦT Q Y
| {z }
pseudo−inverse de Φ

• Si Q = I nous retrouvons l’équation (41).

74
Identification des modèles
Identification d’un modèle déterministe

• Soit un système donné par sa transmittance G (z) :

b1 z−1 + b2 z−2 + ... + bm z−m

G (z) = (44)
1 + a1 z−1 + a2 z−2 + ... + an z−n
• l’équation récursive étant la suivante :

ŷ (t) = −a1 y (t − 1) − ... − an y (t − n) + b1 u (t − 1) + ... + bm u (t − m)

(45)
• qui peut être mise sous forme matricielle

ŷ (t) = ϕT (t) θ

T
ϕ (t) = [−y (t − 1) ... − y (t − n) u (t − 1) ... u (t − m)] (46)

T
θ = [a1 ... an b1 ... bm ]

75
Détermination d’un modèle ARX par la Méthode des Moindres Carrés
Ordinaires (MCO)

• notons :

T
ϕ (t) = [−y (t − 1) ... − y (t − na ) u (t − 1) ...u (t − nb )]
(47)
T
θ = [a1 a2 ... an b1 b2 ... bm ]

• La somme des erreurs quadratiques de prédiction à un pas est donnée

par :
n n
X 2
X
2
J= (y (t) − ŷ (t|t − 1, θ)) = y (t) − ϕT (t) θ (48)
t=1 t=1

• On obtient :

−1
θ̂ = ΦT Φ ΦT Y (49)

76
Biais de l’estimation des paramètres

• Soit le modèle parfait X tel que :

X = Φθ̃ (50)

Où θ̃ est le vecteur des paramètres théoriques.

• Les mesures sont bruitées par e :

Y = X + e = Φθ̃ + e (51)

• La différence entre la valeur estimée et la valeur théorique (vraie

valeur), est donnée par :

−1
θ̂ − θ̃ = ΦT Φ T
 Y − θ̃
Φ

−1

= ΦT Φ ΦT Y − ΦT Φ θ̃

−1   (52)
= ΦT Φ ΦT Y − Φθ̃

−1
= Φ T Φ ΦT e

77
Biais de l’estimation des paramètres

• Biais de l’estimation des paramètres : Le biais d’un estimateur est

l’espérance de l’écart entre l’estimation et les vraies valeurs.
• NB : L’estimateur non-biaisé fournit une valeur non-décalée par
rapport à la vraie valeur
h i h
−1 i
E θ̂ − θ̃ = E ΦT Φ ΦT e (53)
• Si le bruit e est non corrélé à Φ, l’espérance devient :
h i h
−1 i
E θ̂ − θ̃ = E ΦT Φ ΦT E [e]
(54)
T

−1 T
= Φ Φ Φ E [e]

• Rappel : E [XY] = E [X] E [Y] si X et Y sont indépendants.

78
Biais de l’estimation des paramètres

• L’estimateur des moindres-carrés est non-biaisé, càd, que la valeur

estimée des paramètres est proche de la vraie valeur en moyenne si :

E [e] = 0 (55)

• La précédente égalité est vérifiée si le bruit a une valeur moyenne nulle

(bruit blanc).

Si le bruit est coloré, cet estimateur est biaisé

79
Variance de l’estimation des paramètres

• Variance d’une variable aléatoire

n
1X 2
h
2
i
Var (α) = σα2 = (αi − ᾱ) = E (α − E (α))
n
i=1
• Considérons la grandeur vectorielle A telle que :

T
A = [α1 α2 ... αn ]
• La matrice de variance/covariance d’une grandeur vectorielle est
comme suit :
h

T i
ΓA = A − Ā A − Ā


2



α1 − ᾱ1 α1 − ᾱ1 α2 − ᾱ2 ... α1 − ᾱ1 (αn − ᾱn )
| {z } | {z }
Terme de

Variance Terme de Covariance
ΓA = E 


2


α2 − ᾱ2 α1 − ᾱ1 α2 − ᾱ2 ... α2 − ᾱ2 (αn − ᾱn ) 
... ... ... ...
(αn − ᾱn )2



(αn − ᾱn ) α1 − ᾱ1 (αn − ᾱn ) α2 − ᾱ2 ...

80
Variance de l’estimation des paramètres

• Variance de l’estimateur des moindres-carrés (pour un bruit centré) :

  T 
Γθ = E θ̂ − θ̃ θ̂ − θ̃
 
T

−1 T   T
−1 T T
Γθ = E Φ Φ Φ e Φ Φ Φ e
 
T

−1 T T 
T

−1 T
Γθ = E Φ Φ Φ ee Φ Φ Φ

• Pour un bruit non-corrélé à la matrice Φ , la matrice de

variance/covariance de l’estimation est donnée par :
h
−1 i   
−1 T

Γθ = ΦT Φ ΦT E eeT Φ ΦT Φ

81
Variance de l’estimation des paramètres

• Cas de bruit blanc sur les mesures

• chaque mesure est supposée entachée d’un bruit blanc additif
• il n’y a pas de corrélation du bruit de mesure entre 2 mesures

e2 (t1 ) σ12
   
0 ... 0 0 ... 0
 T  0 e2 (t2 ) ... 0  =  0
  σ22 ... 0 
E ee = E 
 ...

... ... ...   ... ... ... ... 
0 0 ... e2 (tn ) 0 0 ... σn2

• Pour un bruit constant sur chaque mesure, de variance σY2 , le bruit sur
l’estimation des paramètres du modèle est :
h
−1 iT
Γθ = σY2 ΦT Φ

82
Estimateur du maximum de vraisemblance

• On prend comme matrice de pondération l’inverse de la matrice de

variance/covariance du bruit, ce qui a pour effet de minimiser le poids
des mesures à bruit élevé.

Γbruit = E eeT
 

J = εT Qε = εT Γ−1 bruit ε

−1
θ̂ = ΦT Γ−1 bruit Φ ΦT Γ−1 bruit Y

83
Rejet des mesures aberrantes

• Le calcul des paramètres du modèle permet de calculer X̂ tel que :

X̂ = Φθ̂
• La variance de l’estimation est :
h i h i h i
σX̂2 = E X̂ X̂ T = E Φθ̂θ̂T ΦT = ΦE θ̂θ̂T ΦT = ΦΓθ ΦT
• Pour des mesures bruitées par un bruit blanc
h i h
−1 iT T
σX̂2 = E X̂ X̂ T = σY2 Φ ΦT Φ Φ

• Rejet des mesures y (ti ) dont l’écart avec le modèle x̂ (ti ) est supérieur à
l’écart type σX̂ avec k fixé par l’utilisateur ( de l’ordre de 2 à 3).

84
Méthode des Variables Instrumentales (IV)

• Supposons que les données répondent à la loi suivante :

y (ti ) = ϕT (ti ) θ̃ + e (ti ) (56)

• Soit z (ti ) vecteur des variables instrumentales, un vecteur colonne de
même dimension que ϕ (ti ) et non corrélé à e (ti ) (z (ti ) et e (ti ) sont
orthogonaux).

n n
1X 1X h i
lim z (ti ) e (ti ) = lim z (ti ) y (ti ) − ϕT (ti ) θ̃ = 0 (57)
n→∞ n n→∞ n
i=1 i=1
" #−1 " n #
 Xn X 
lim z (ti ) ϕT (ti ) z (ti ) y (ti ) = θ̃ (58)
n→∞  
i=1 i=1

85
Méthode des Variables Instrumentales (IV)

• Choix de z (ti ) :
• L’expression la plus simple de z (ti ) est la suivante :

z (t) = [u (t − 1) ... u (t − na − nb )]
• A moins d’avoir une contre-réaction de la sortie vers l’entrée, ce choix
permet de garantir une totale indépendance de z (ti ) par rapport au bruit
e (ti ).
• Cependant il est important de mentionner que ce choix présente un
problème de stabilité numérique.

86
Méthode des Variables Instrumentales (IV)

• Choix pratique de z (ti ) :

• On génère ym la sortie simulée du système. Pour cela, on utilise la valeur
de θlsm donnée par un algorithme des moindres carrés :

pour i = 1...n,

z (ti ) = [−ym (ti − 1) ... − ym (ti − na ) u (ti − 1) ... u (ti − nb )]T

ym = zT (ti ) θlsm

fin
• θiv l’estimation de θ par la méthode IV est donnée par :

−1
θiv = Z T Φ ZT Y
(59)
Z = zT (ti ) ... zT (tn )
 

87
Moindre-Carrés Récursifs
Position du problème

• Calcul sur n mesures :

• On dispose de n observations auxquelles correspondent n valeurs du
modèle

   
y (t1 ) f1 (t1 ) f2 (t1 ) ... fk (t1 )
 y (t2 )   f1 (t2 ) f2 (t2 ) ... fk (t2 ) 
Yn =  ; Xn = Φn θ =  θ
   
 ...   ... ... ... ... 
y (tn ) f1 (tn ) f2 (tn ) ... fk (tn )
(60)
• la solution optimale calculée sur les n mesures est :
 −1
θ̂n = ΦTn Φn ΦTn Yn (61)
• elle minimise le critère calculé sur ces n mesures :

Jn = εTn εn = (Yn − Xn )T (Yn − Xn ) (62)

88
Position du problème

• Inconvénients de la méthode classique :

• La dimension (n x k) de la matrice Φ croı̂t avec le nombre de mesures n .
• L’ajout d’une (n + 1)ème mesure impose de recommencer la totalité du
calcul.

• Avantages de la méthode récursive :

• Permet de traiter les informations en ligne au fur et à mesure.
• Permet de répartir dans le temps la charge CPU en allégeant le calcul de la
matrice pseudo-inverse.

89
Mise en évidence de la solution

• Utiliser la technique de calcul

θ̂n = R−1
n Qn
(63)
avec Rn(k x k) = ΦTn Φn et Qn(k x 1) = ΦTn Yn
ème
• Ajout de la (n + 1) mesure :

y (t1 ) f1 (t1 ) f2 (t1 ) ... fk (t1 )

   

 y (t2 ) 


 f1 (t2 ) f2 (t2 ) ... fk (t2 ) 

Yn+1 =
 ...  ; Xn+1 = Φn+1 θ = 
  ... ... ... ... θ

 y (tn )   f1 (tn ) f2 (tn ) ... fk (tn ) 
y (tn+1 ) f1 (tn+1 ) f2 (tn+1 ) ... fk (tn+1 )
(64)
ème
• L’ajout d’une (n + 1) mesure impose de recommencer la totalité du
calcul.

90
Mise en évidence de la solution

• Soit :
h i
ϕT (n + 1) = f1 (tn+1 ) f2 (tn+1 ) ... fk (tn+1 ) (65)
• On dispose de n + 1 observations auxquelles correspondent n + 1
valeurs du modèle :

! ! !
Yn Xn Φn
Yn+1 = ; Xn+1 = = θ
y (tn+1 ) x (tn+1 ) ϕT (n + 1)
(66)
• Posons :

−1
θ̂n+1 = ΦTn+1 Φn+1 ΦTn+1 Yn+1 (67)

91
Mise en évidence de la solution

• Calculons ΦTn+1 Φn+1

!
  Φn
ΦTn+1 Φn+1 = ΦTn ϕ (n + 1) T
ϕ (n + 1)
(68)

= ΦTn Φn + ϕ (n + 1) ϕT (n + 1)
• Calculons ΦTn+1 Yn+1
!
  Yn
ΦTn+1 Yn+1 = ΦTn ϕ (n + 1)
y (tn+1 )
(69)

= ΦTn Yn + ϕ (n + 1) y (tn+1 )

92
Mise en évidence de la solution

• Expression récursive de base

−1
θ̂n+1 = ΦTn Φn + ϕ (n + 1) ϕT (n + 1) ΦTn Yn + ϕ (n + 1) y (tn+1 )

(70)
• Posons

−1
Pn = R−1 T
n = Φ n Φn

−1
Pn+1 = R−1 T T
n+1 = Φn Φn + ϕ (n + 1) ϕ (n + 1)

93
Détermination de la relation entre Pn+1 et Pn

• Afin d’éviter de calculer à chaque pas l’inverse d’une matrice, on a

recours au lemme d’inversion matricielle suivant :

−1 −1
= A−1 − A−1 B DA−1 B + C−1 DA−1

[A + BCD] (71)
• Posons

A = P−1 T
n = Φn Φn

B = D−1 = ϕ (n + 1)

C=1

94
Détermination de la relation entre Pn+1 et Pn

• L’application du lemme d’inversion sur Pn+1 permet d’écrire :

h i−1
Pn+1 = Pn − Pn ϕ (n + 1) I + ϕT (n + 1) Pn ϕ (n + 1) ϕT (n + 1) Pn
| {z }
Kn+1
(72)
• On définit le gain de Kalman Kn+1 tel que :
−1
Kn+1 = Pn ϕ (n + 1) I + ϕT (n + 1) Pn ϕ (n + 1)

(73)
d’où

Pn+1 = Pn − Kn+1 ϕT (n + 1) Pn = I − Kn+1 ϕT (n + 1) Pn

 
(74)

95
Calcul récursif évolué

• En se référant à l’équation (70), on a :

θ̂n+1 = Pn+1 ΦTn Yn + ϕ (n + 1) y (tn + 1)

= Pn − Kn+1 ϕT (n + 1) Pn ΦTn Yn + ϕ (n + 1) y (tn + 1)

 

= Pn ΦTn Yn + Pn ϕ (n + 1) y (tn + 1) − Kn+1 ϕT (n + 1) Pn ΦTn Yn

− Kn+1 ϕT (n + 1) Pn ϕ (n + 1) y (tn + 1)
 

= θ̂n + Kn+1 y (tn + 1) − ϕT (n + 1) θ̂n 

 
| {z }
ŷ(tn +1)
| {z }
erreur de prédiction

96
Algorithme des Moindres Carrés Récursifs (Recursive Least Squares)

1. Calcul du gain de correction

−1
Kn+1 = Pn ϕ (n + 1) I + Pn ϕ (n + 1) ϕT (n + 1)


2. Actualisation de P

Pn+1 = Pn − Kn+1 ϕT (n + 1) Pn = I − Kn+1 ϕT (n + 1) Pn

 

3. Prédiction de la sortie

ŷ (tn + 1) = ϕT (n + 1) θ̂n
4. Actualisation des paramètres

θ̂n+1 = θ̂n + Kn+1 (y (tn + 1) − ŷ (tn + 1))

97
Initialisation

• La mise en œuvre pratique de la récurrence évoluée nécessite

l’initialisation de la matrice P, et le vecteur des paramètres θ̂.
• Désignant par :
• P00 matrice carrée de dimension k x k
• θ̂00 vecteur colonne k x 1

−1
• Initialisation de P00 : Sachant que P00 = ΦT00 Φ00 , initialement Φ00
est inexistante (absence de lignes)

98
Initialisation

• Première Méthode : Acquisition de k mesures, de sorte que

R00 = ΦT00 Φ00 soit inversible. Initialiser avec :

−1 T

−1
 P00 = RK = ΦK ΦK



θ̂00 = P00 ΦTK YK


• Deuxième méthode : P étant de dimension constante k x k , on propose

le choix d’initialisation suivant :
(
P00 = GIK avec G = 10 ou 100
θ̂00 : choix arbitraire, souvent on prend θ̂00 = 0

99
Initialisation

• Première Itération
 
K1 = P00 ϕ (1) 1 + ϕT (1) P00 ϕ (1)
 
P1 = I − K1 ϕT (1) P00
ŷ (1) = ϕT (1) θ̂00
θ̂1 = θ̂00 + K1 [y (1) − ŷ (1)]

• Deuxième Itération
 
K2 = P1 ϕ (2) 1 + ϕT (2) P1 ϕ (2)
 
P2 = I − K2 ϕT (2) P1
ŷ (2) = ϕT (2) θ̂1
θ̂2 = θ̂1 + K2 [y (2) − ŷ (2)]

• Récurrence générale

100
Méthode des moindres Carrés
récursifs avec facteur d’oubli
Méthode des moindres Carrés récursifs avec facteur d’oubli

• Objectif :
Diminuer les ”poids” des anciennes mesures au profit
des plus récentes.

La méthode des moindres carrés récursifs avec facteur d’oubli exponentiel

est utilisée si les paramètres du modèle évoluent de manière lente.

101
Méthode des moindres Carrés récursifs avec facteur d’oubli

• Le critère quadratique pondéré J étant le suivant :

n+1
X
λn−i+1 y (ti ) − ϕT (i) θ


J= (75)
i=1

• Le vecteur des paramètres est obtenu en minimisante J tel que :

" n+1 #
X
λn−i+1 y (ti ) − ϕT (i) θ


θ̂ (n + 1) / min(J) = min (76)
i=1

102
Méthode des moindres Carrés récursifs avec facteur d’oubli

• on définit de la même façon que précédemment :

n+1 n
λn−i+1 ϕ (i) ϕT (i) = λn−i+1 ϕ (i) ϕT (i) +ϕ (n + 1) ϕT (n + 1)
P
P  
Rn+1 =
i=1 i=1

= λRn + ϕ (n + 1) ϕT (n + 1)
(77)
et
n+1 n
λn−i+1 (ϕ (i) y (i)) = λn−i+1 [ϕ (i) y (i)] +ϕ (n + 1) y (n + 1)
P P
Qn+1 =
i=1 i=1

= λQn + ϕ (n + 1) y (n + 1)
(78)
Rn(kxk) = ΦTn Φn et Qn(kx1) = ΦTn Yn

103
Méthode des moindres Carrés récursifs avec facteur d’oubli

θ̂n+1 = R−1
n+1 Qn+1 = Pn+1 Qn+1

−1
= λRn + ϕ (n + 1) ϕT (n + 1) (λQn + ϕ (n + 1) y (tn+1 ))
(79)

−1
= λP−1 T
n + ϕ (n + 1) ϕ (n + 1) (λQn + ϕ (n + 1) y (tn+1 ))

−1
= λP−1 T
n + ϕ (n + 1) ϕ (n + 1) (λQn + ϕ (n + 1) y (tn+1 ))

• Soit

−1
Pn+1 = λP−1 T
n + ϕ (n + 1) ϕ (n + 1)

• Posons
A = λRn = λP−1
n

B = DT = ϕ (n + 1)

C=1
104
Méthode des moindres Carrés récursifs avec facteur d’oubli

Pn ϕ(n+1)
Kn+1 = λ+ϕT (n+1)Pn ϕ(n+1)

Pn+1 = λ−1 Pn − Kn+1 ϕT (n + 1) Pn

(80)

 
θ̂n+1 = θ̂n + Kn+1 y (tn+1 ) − ϕT (n + 1) θ̂n

105
Détermination d’un modèle
ARMAX par la Méthode des
Moindres Carrés Étendus
Introduction

• La structure du modèle ARMAX étant la suivante :

- u (t) : entrée du système

- y (t) : sortie du système
- e (t) : bruit blanc
- v (t) : bruit filtré

106
Introduction

• La structure du modèle ARMAX étant la suivante :

- u (t) : entrée du système

- y (t) : sortie du système
- e (t) : bruit blanc
- v (t) : bruit filtré

107
Travaux Dirigés
TD Prédiction

Exercice 1 :
Soit le système suivant :

1. Montrer qu’on peut mettre le système sous forme ARMAX ?

Définir A, B et C. Peut-on aussie le mettre sous forme ARX ?
2. Ecrire l’équation récursive liant y (t) à y (t − 1), u (t − 1) et e (t).
3. Déterminer la meilleure prédiction à 1 pas pour ce système.
4. Retrouver ce résultat à l’aide de la formule générale.

108
TD Prédiction

Exercice 2 :
Soit le système suivant :

1. Trouver la structure de modèle correspondante.

2. Déterminer la prédiction à un pas.

109
TD Prédiction

Exercice 3 :
Soit le système suivant :

1. Trouver la structure de modèle correspondante.

2. Déterminer la prédiction à un pas.

110
TD Prédiction

Exercice 4 :
Soit le système suivant :

q−1
y (t) = u (t) + e (t)
1 − 0.5q−1
Avec :
e (0) = 0.1; e (1) = −0.05; e (2) = −0.03
e (3) = 0.08; e (4) = −0.1; e (5) = −0.01 et u (t) est un échelon .

1. En supposant les conditions initiales nulles, déterminer les six premiers

échantillons de la sortie y (0) ... y (5).
2. Calculer les valeurs simulées (en absence de bruit) des mêmes échantillons.
3. Vérifier que l’écart entres valeurs réelles et celles simulées est égal à e (t).
4. Calculer la prédiction à 1 pas ŷ (t|t − 1) pour t = 1...5.
5. Comparer les résultats de la simulation et de la prédiction en utilisant les
mesures bruitées.
6. Recalculer les prédictions à 1 pas à partie de t = 2 111