Tp1 : Methodes de la moindre carre récursifs

Objectif du TP
 Analyser la méthode des moindres carrés généralisés pour identifier sans biais des
modèles (système + perturbation).
 Appliquer cette méthode à un procédé de régulation de la température d’un four.
Introduction
L’identification consiste à déterminer les paramètres d’un modèle mathématique en fonction
d’un critère spécifique. Les paramètres sont obtenus en minimisant l’écart entre le signal mesuré
et le signal estimé, en utilisant des critères comme les moindres carrés ou l’erreur quadratique
moyenne. Nous nous concentrons ici sur la méthode des moindres carrés récursifs, qui permet
d'identifier en temps réel les paramètres d’un modèle d’un système ou d’un procédé industriel
(par exemple, les machines à courant continu ou les machines synchrones).
Cette méthode vise à minimiser l'erreur quadratique moyenne (EQM) entre la sortie du modèle et
celle du système observé (comme une machine asynchrone) pour ajuster les paramètres du
modèle. L’algorithme récursif permet de mettre à jour les paramètres sans devoir recommencer
les calculs à chaque nouvelle observation, ce qui le rend efficace pour des applications en temps
réel.
Nous appliquons cette méthode à un système modélisé par un ensemble d’équations linéaires, où
l'objectif est d'identifier les paramètres du modèle à partir de mesures d’observations.
L'algorithme des moindres carrés est utilisé pour résoudre ce problème d'identification de
manière optimale.
Travail demande :
Le but de ce travail est d'étudier la régulation de la température d'un four à chauffage en utilisant
l'algorithme des moindres carrés récursifs. Ce processus consiste à identifier en temps réel les
paramètres du modèle de régulation en minimisant l'écart entre la température mesurée et celle
estimée.
L'objectif est d'améliorer la précision et l'efficacité du système de commande malgré les
perturbations.

A/
Modélisation :
1. Calcul de la fonction de transfert :
Après avoir appliqué la transformée de Laplace, on obtient une équation reliant
Δ T ( p)
F (p)=
Uf ( p)
ΔT ( p ) K ( f )∗R
En simplifiant, on arrive à la fonction de transfert : F (p)= = ,
Uf ( p) 1+ RCP
Avec : R=0.1°C/W
C=2000J/°C

Cela correspond à un système de premier ordre

qe (p) = k f * U f(p)
 1
qs (p) = * ∆T(p)
R
1 1 1
∆T(p) = *(q e (p)- q s(p)) = (k f*Uf (p) - *∆T(p))
C.p C.p R
1 1
Donc ∆T(p) + ∆T(p) = k f*U f(p)
RC . p C.p
1
∗k
1 1 ΔT ( p) C.p f k f ∗R
∆T(p) *(1 + )= k f*U f(p)  F (p) = = =
RC . p C.p U ( p) 1 1+ RC . p
1+
RC . p
2. Parameters du système :
o Le gain est calculé comme kc=kf*R=40*0.1=4
o La constante de temps est τ=R⋅C=2000*0.1=200 secondes
3. Fonction de transfert pour le signal mesuré :
On relie ensuite le signal mesuré Umes(p)= U(t)*Kb) = k b * U t = k b * k t * ∆T(p)
= k b * k t * Uf (p) * F (p) avec Uf (p) = k a* UC
4
Umes (p) = k a * k b* k t * UC *
1+ 200 p
Finalement, on exprime la température T(p) en fonction de Uc(p)
Umes( p) 4 2
T (p) = = 10*1250*40*10-6 * =
Uc ( p) 1+ 200 p 1+ 200 p
 Cela montre que le système peut être modélisé comme un filtre de premier ordre, avec
une réponse caractéristique à une entrée Uc(p)
B/ Identification :
Umes( p) 2
1) T(p)= =
Uc ( p) 1+ 200 p
T ( p) 2
T(z) =Z {P(p)T(p)} = (1-z-1) Z { } = (1-z-1) Z { }
p p (1+200 p)
1/200 1−z
= 2(1-z-1)Z{ }=2 z ¿¿ =2 ¿ ¿
p (1/200+ p) z
b 0.1903
 Avec Te=20s on trouve T (z) = =
z−a z−0.9048
 −1
0. 4424 0. 4424∗z
T(z)= =
z−0. 7788 1−0. 7788∗z −1

Umes ( z )∗( 1−0. 7788∗z −1 ) =Uc ( z ) (0. 4424∗z−1)

Umes(k)- 0.7788*Umes(k-1)=0.4424*Uc(k-1) Umes(k)=0.7788*Umes(k-1)+0.4424*Uc(k-1)
Par identification, on conclut que :
Umes(k)=a*Umes(k-1)+b*Uc(k-1)
a=0.7788
b=0.4424

2) Le signal d'excitation est un SBPA (Signal Bipolaire à Période d'Amplitude Unitaire)

centré, d'amplitude unitaire et de longueur N=500. La génération automatique d'un SBPA
sous MATLAB est réalisée grâce à la commande « idinput ».
PROGRAMME

Type : 'sine'

signal d`excitation Uc signal de sortie Umes

1 1.5

0.8
1
0.6
0.5
0.4

0.2 0

0
-0.5
-0.2

-0.4 -1

-0.6
-1.5
-0.8

-1 -2
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Iteration Iteration
 Type : 'prbs'
signal d`excitation Uc 3) signal de sortie Umes
1 2

0.8
1.5

0.6
1
0.4
0.5
0.2

0 0

-0.2
-0.5

-0.4
-1
-0.6
-1.5
-0.8

-1 -2
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Iteration Iteration

La détermination du vecteur de paramètres peut être réalisée à l'aide de plusieurs méthodes.

Parmi celles-ci, on trouve la méthode des moindres carrés non récursive, qui s'effectue
selon l'algorithme suivant :
 trajectoires E/S
2

1.5

0.5

-0.5

-1

-1.5

-2
0 50 100 150 200 250 300 350 400 450 500
Iteration

4) On implémente une fonction sous Matlab, θ = mcr(Umes, Uc, P_0, θ_0), qui estime les
paramètres du modèle à partir des mesures (Umes, Uc) et des valeurs initiales (P_0, θ_0),
afin de fournir une estimation optimale du vecteur de paramètres θ.

PROGRAMME:
  MC_R

>>clc;clear all;close all;

% declaration des donneés
a=0.78;
b=0.44;
N=500;
Umes(1)=0;
p=10^6*eye(2); theta_e=[0 0]';Umes_e(1)=Umes(1);
theta_i=theta_e;pi=p;
% generation du signal d'excitation
Uc=idinput(N,'prbs',[0 0.1],[-1 1]);
figure(1);plot(Uc);grid;xlabel('Iteration');
% generation du signal de sortie
for k=2:length(Uc)
Umes(k)=a*Umes(k-1)+b*Uc(k-1);
end
hold on
plot(Umes);
title('trajectoires E/S');
% estimation
[theta_e,theta_ef,e,trp,Umes_e]=mcr(Uc,Umes,theta_e,p,Umes_e)
% visualisation
ar=a*ones(1,length(Uc));
kr=b*ones(1,length(Uc));
theta_ef(:,1)=theta_i;trp(1)=trace(pi);
figure(2);
plot(theta_ef(1,:),'b');
hold on
plot(ar,'b');grid;xlabel('k');
hold on
plot(theta_ef(2,:),'r');
hold on
plot(br,'r');
grid;xlabel('itérations');legend('a_e','b_e');
title('Evolution des parametres estimes');
figure(3);
plot(e);
grid;xlabel('iteration');
title('Erreur de prediction');
figure(4);
plot(trp);
grid;xlabel('iteration');
title('trace de la matrice du gain de correction');
figure(5);
plot(Umes,'b');
hold on
plot(Umes_e,'r');
grid;xlabel('iteration');
legend('Umes','Umes estime');
title('comparaison Usortie reelle et Usortie estime');

 mc_r

function [theta_e,theta_ef,e,trp,Umes_e]=mcr(Uc,Umes,theta_e,p,Umes_e)
for k=2:length(Uc)
% construction de la matrice d'observation
h=[Umes(k-1) Uc(k-1)];
% mise à jour de la matrice du gain d'adaptation
p=p-p*h'*h*p/(1+h*p*h')
% mise a jour du parametres estimes
theta_e=theta_e+p*h'*(Umes(k)-h*theta_e);
% enregistrement des parametres estimes
theta_ef(:,k)=theta_e;
% erreur d'estimation
e(k-1)=Umes(k)-h*theta_e;
% trace de la matrice du gain de correction
trp(k)=trace(p);
% sortie estime
Umes_e(k)=theta_e(1)*Umes_e(k-1)+theta_e(2)*Uc(k-1);
end
end

Interpretation:
Dans le cadre de l'estimation continue et de l'ajustement des paramètres du modèle, chaque étape
de l'algorithme joue un rôle important :

 Partie d'initialisation : Cette étape consiste à initialiser le vecteur de paramètres ainsi

que le gain d'adaptation P0.Ce dernier permet de corriger l'erreur de prédiction à chaque
instant.
 Partie génération de l'entrée : Cette étape est identique à celle de la méthode des
moindres carrés non récursifs, où l'entrée est générée pour le modèle.
 Partie génération de la sortie : Cette partie est également similaire à celle de la méthode
des moindres carrés non récursifs, où la sortie du modèle est calculée.
 corps de l'algorithme : C'est ici que, à chaque instant, les éléments suivants sont
déterminés :
o Le vecteur de paramètres
o Le gain d'adaptation P
o L'erreur de prédiction e
Analyse des courbes :
 1ere cas : θ initiale =[0,0]’ et p0=10^6*eye(2)
{‘avec eye(2) : c’est la matrice d’identité 2*2’} :
-Vitesse de correction est rapide
- On remarque une chute en 3éme valeur donc c’est la convergence.
-On remarque que les 2 courbes Umes réelle et estimé sont superposé (tant que u mes est
inférieur à 10-4 donc c’est un modèle parfait (choix des valeurs initiaux est parfait)

θ
 2éme cas : initiale =[100,100]

θ
D`après les courbes obtenu par initiale =[100,100], On remarque qu'en changeant le vecteur de
paramètres initial, les valeurs estimées de a et b se rapprochent des valeurs réelles, mais plus
rapidement que dans le premier cas. La tension de sortie estimée et réelle deviennent presque
identiques. Par contre, le gain d'adaptation descend vers zéro, mais plus lentement que dans le
premier cas.

θ
 3éme cas : initiale =[0,0]’ et p0=10^2*eye(2)
 -On remarque que dans ce cas, le choix de ce gain est incorrect, car les valeurs estimées ne
convergent pas vers les valeurs réelles.
-L'erreur augmente également.

 La méthode des moindres carrés non récursive donne un vecteur de paramètres proche du
réel, mais nécessite beaucoup d'itérations, ce qui rend les calculs lourds.

Donc pour simplifier, nous utiliserons la méthode des moindres carrés récursive, qui
est plus efficace.

Conclusion
En résumé, la méthode des moindres carrés non récursifs est un outil puissant pour l'estimation
des paramètres d'un modèle mathématique représentant un procédé industriel. Elle permet de
trouver les valeurs des paramètres qui minimisent l'écart entre les observations expérimentales et
les résultats prévus par le modèle, assurant ainsi une bonne approximation du comportement du
système étudié. Toutefois, cette méthode peut être limitée par le nombre élevé d'itérations
nécessaires, ce qui peut rendre les calculs longs et complexes, surtout pour des systèmes de
grande taille. Malgré ces contraintes, elle reste une méthode privilégiée dans les situations où la
précision de l'estimation est essentielle.
Pour rendre les calculs plus rapides, il est possible d'utiliser des méthodes alternatives, comme
les moindres carrés récursifs

Tp2: Application des moindres carres récursives sur une

maquette didactique DIDALAB
de régulation de débit de température d`air

Introduction

Dans ce TP, nous travaillons avec un modèle représenté par l'équation Y=Hθ+e tel que :
Y : est un vecteur de mesures,
Θ : est un vecteur de paramètres que nous devons estimer,
H : est une matrice connue et
e : représente les erreurs ou le bruit dans les mesures.
Si le vecteur e contient des éléments qui sont liés entre eux, cela signifie que le bruit est "coloré".
Dans ce cas, la méthode des moindres carrés ordinaires utilisée dans le TP1 ne donne pas de bons
résultats, car elle introduit des erreurs dans l'estimation des paramètres.
Pour éviter ce problème, il existe des méthodes plus adaptées. Ces méthodes supposent que le
bruit dans les mesures peut être vu comme la sortie d'un filtre linéaire inconnu, excité par du
bruit blanc non corrélé.
Deux méthodes populaires pour obtenir des estimations sans biais des paramètres sont les
moindres carrés généralisés (MCG) et la méthode de la variable instrumentale (VI).
L'objectif de ce TP est de comprendre ces méthodes. Nous commencerons par tester l'algorithme
des moindres carrés récursifs avec du bruit blanc, puis avec du bruit coloré. Ensuite, nous
étudierons les moindres carrés généralisés pour identifier des modèles sans biais. Enfin, nous
appliquerons ces méthodes pour estimer les paramètres d'un modèle de régulation de la
température d'un four.
Travail demande :
Nous allons ajuster le programme précédent pour y ajouter du bruit au signal. En effet, c'est la
commande awgn(y,k) qui permet d'introduire ce bruit.et nous minimiser le longueur ‘N=300’

PROGRAMME:

clc;clear all;close all;

% declaration des donneés
a=0.78;b=0.44;N=300;Umes(1)=0;
p=10^6*eye(2); theta_e=[0 0]';Umes_e(1)=Umes(1);
theta_i=theta_e;pi=p;
% generation du signal d'excitation
Uc=idinput(N,'prbs',[0 0.1],[-1 1]);
figure(1);plot(Uc);grid;xlabel('Iteration');
% generation du signal de sortie
for k=2:length(Uc)
Umes(k)=a*Umes(k-1)+b*Uc(k-1);
end
hold on
plot(Umes);title('trajectoires E/S');
Umes_sb=Umes; %sans bruit
Umes=awgn(Umes,30);
figure(2); plot(Umes_sb,'b');
hold on
plot(Umes,'r');grid; xlabel('Iteration');
legend('Umes sans bruit','Umes avec bruit');
title('comparaison de Umes avec bruit et sans bruit');
% estimation
[theta_e,theta_ef,e,trp,Umes_e]=mcr(Uc,Umes,theta_e,p,Umes_e)
% visualisation
ar=a*ones(1,length(Uc));
kr=b*ones(1,length(Uc));
theta_ef(:,1)=theta_i;trp(1)=trace(pi);
figure(3);
plot(theta_ef(1,:),'b');
hold on
plot(ar,'b');
hold on
plot(theta_ef(2,:),'r');
hold on
plot(kr,'r');
grid;xlabel('itérations');legend('a_e','b_e');
title('Evolution des parametres estimes');
figure(4);
plot(e);
grid;xlabel('iteration');
title('Erreur de prediction');
figure(5);
plot(trp);
grid;xlabel('iteration');
title('trace de la matrice du gain de correction');
figure(6);
plot(Umes,'b');
hold on
plot(Umes_sb,'r');
hold on
plot(Umes_e,'g');
grid;xlabel('iteration');
legend('Umes','Umes_sans bruit','Umes estime');
title('comparaison sortie bruitée , sorie non bruitée et sortie estimée');

Analyse des courbes :

1ere cas : SNR=30
o Nous obtenons les courbes suivantes :
On interprète :
 Erreur de prédiction oscille
 On troupe un modèle approximative
 Gain statique ks = 2
 Gain statique estimé ks estimé = 1.999
 Erreur acceptable

2ieme cas :SNR=20

o Nous obtenons les courbes suivantes :
On interprète :

 pas de convergence
 pas de superposition
 erreur est grand Déformation totale
 Gain statique ks = 2
 Gain statique estimé ks estimé = 1.986
 Erreur acceptable
3ieme cas: SNR=50
 Trace prend plus de temps pour la chute
 On trouve le modèle mais après un peut du temps  bonne identification

 L'algorithme montre des limites pour un SNR > 20, où les sorties estimée et réelle se
confondent, ce qui n'est pas fréquent en pratique. Pour un SNR = 10, bien que l'erreur soit
 plus importante, les résultats restent raisonnables. À un SNR = 30, les sorties estimée et
réelle sont presque identiques, et pour un SNR > 50, l'algorithme est très précis, mais un
tel SNR est rarement observé dans les situations réelles.

Conclusion
Les méthodes classiques des moindres carrés, qu'elles soient récursives ou non, échouent à
donner de bonnes estimations lorsque le système est exempt de bruit ou que le bruit est blanc.
Dans ce contexte, une approche plus avancée, comme la méthode des moindres carrés
récursifs étendus, est nécessaire pour améliorer les résultats.