Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
TP2 :Méthode des Moindres Carrés Généralisés
I. But :
on va étudier dans ce TP l’influence d’un bruit blanc dans la premiére partie en utilisant la méthode
des moindres carrés récursives et du bruit coloré dans la deuxiéme partie en utilisant la méthode
des moindres carrés généralisés.
II. Méthode des moindres carrés généralisé :
En présence de bruit sur le systéme la méthode des moindres carrés donne un estimateur biaisé
c’est pour ce là on utilise une nouvelle méthode appelé a méthode des moindres carrés généralisés.
Cette méthode consiste à introduire un filtre sur le vecteur mesure de façon à blanchir le bruit ce
qui autorise l’utilisation de [Link] paramètres du filtre sont eux aussi à estimer, ce qui conduit à
utiliser à deux reprises les MCR, comme le montre les étapes suivantes :
1 ère étape : choisir le filtre a priori
F(z)=1+f1(z^(-1)) + f2(z^(-2)) +………………..+ fn(z^(-n))
2 éme étape : calcul des valeurs filtrées :
Uf(k)=f(z)*u(k)
Yf(k) = f(z)*y(k)
3 éme étape : déterminer thêta (k) par la méthode des MCR appliqué au modèle.
a(z)*yf(k)=b(z)*uf(k)+e(k)
4éme étape :calculer les termes résiduels :
epsilon(k)=A(z)*v(k)-B(z)*u(k)
5 éme étape :estimer les coefficients du filtre f(z) par la méthode des MCR appliqué au modèle :
f(z)*epsilon(k)=e(k)
6 éme étape : reprendre l’étape 2 avec les nouvelles valeurs calculées du filtre F(z)
III. Travaille demandé : >> T=tf(2,[200 1])
1. Identification: T=
2
discrétisation la fonction de transfert T ( p ) = 2
200 p+1
avec une période d’échantillonnage T e =20 s ---------
200 s + 1
en utilisant la fonction c2d :
Continuous-time transfer function.
2EME ÉLECTRIQUE G1 TP2 1
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
Umes(z ) 0,1903 0,19∗z −1 0,19∗z −1 >> Td=c2d(T,20)
T(z) = = ≈ =
Uc (z) z−0,9048 (z−0,905)∗z−1 1−0.905∗z−1
Td =
U mes ( z )∗( 1−0.905∗z ) =U c ( z ) z ∗0.19
−1 −1
0.1903
U mes ( z )−U mes ( z ) z −1∗0.905=U c ( z ) z −1∗0.19
----------
U mes ( z )=U mes ( z ) z −1∗0. 905+U c ( z ) z−1∗0.19
U mes ( k )=0.905∗U mes (k −1)+0.19∗U c (k −1) z - 0.9048
Soit y ( k )=U mes ( k ) et U ( k −1 )=U c (k −1) Sample time: 20 seconds
Discrete-time transfer function.
Par identification, on obtient:
a1 réelle = 0,9048
b1 réelle = 0,1903
généralisation du signal d’excitation U c centrée d’amplitude unitaire et de langueur N=100 en
utilisant la fonction idinput
clc; clear all; close all;
%initialisation
a=0.905; b=0.19;y(1)=0;
%génération du signal d’excitation u
u=idinput(100,'prbs',[0 1],[0 1]);
2. systéme affecté d’un bruit e(k) :
1ér cas : e(k) est un bruit blanc :
étude de la robustesse de l’algorithme MCR pour différents niveaux de SNR (Signal to Noise Radio)
pour le bruit cible. En utilisant la fonction awgn (Add White Gaussian Noise to signal)
Syntaxe : y = awgn (x, SNR)
Description : awgn (x, SNR) ajoute un bruit gaussien blanc au signal vectoriel x. Le scalaire
SNR spécifie le rapport signal sur bruit par échantillon, en dB. Si x est complexe, awgn ajoute
un bruit complexe. Cette syntaxe suppose que la puissance de x est égale à 0 dBW.
2EME ÉLECTRIQUE G1 TP2 2
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
clc; clear all; close all;
%initialisation
a=0.905; b=0.19;y(1)=0; theta=[0 0]';p=(10^6)*eye(2);
thetaf(:,1)=theta;ye(1)=y(1);
%génération de signald’excitation u
u=idinput(100,'prbs',[0 1],[0 1]);
figure(1); plot(u); grid; xlabel('itération'); title('signal d''excitation');
%génération de signal de sortie y
for k=2:length(u)
y(k)=a*y(k-1)+b*u(k-1);
end
ysb=y;
y=awgn(y,50 );
figure(2); plot(y); grid; xlabel('itération'); title('signal de sortie');
%algo mc récursif
for k=2:length(u)
h=[y(k-1) u(k-1)];
p=p-((p*h'*h*p)/(1+h*p*h'));
theta=theta+p*h'*(y(k)-h*theta);
thetaf(:,k)=theta;
e(k)=y(k)-h*theta;
ye(k)=theta(1)*ye(k-1)+theta(2)*u(k-1);
trp(k)=trace(p);
end
a1=a*ones(1,length(u));
b1=b*ones(1,length(u));
figure(3); plot(thetaf(1,:),'b');
hold on
plot(thetaf(2,:),'r');
hold on
plot(a1,'b');
hold on
plot(b1,'r');
grid; xlabel('itération'); title('paramétres estimés');
figure(4);plot(e); grid; xlabel('itération'); title('erreur de préduction');
figure(5); plot(trp); grid; xlabel('itération'); title('trace de la matrice de gain
d''adaptation');
figure(6); plot(y,'b');
hold on
plot (ysb,'r');
grid; xlabel('itération'); title('comparaison y réelle et y
bruité');legend('y','ysb');
figure(7),plot(y,'b');
hold on
plot(ye,'r');grid;xlabel('itération'); title('comparaison y réelle et y
estimé');legend('y','ye');
2EME ÉLECTRIQUE G1 TP2 3
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
On prond SNR=50 :
2EME ÉLECTRIQUE G1 TP2 4
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
2EME ÉLECTRIQUE G1 TP2 5
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
On prond SNR=30 :
2EME ÉLECTRIQUE G1 TP2 6
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
On prond SNR=10 :
2EME ÉLECTRIQUE G1 TP2 7
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
Interprétation :
plus que la valeur de SNR diminue plus que le systéme devient de plus en plus robuste ainsi les
valeurs estimés s’éloigne de l’intervalles de justesse et par suite cette algorithme perte son
justesse.
Par exemple pour une valeur de SNR=10 les valeurs estimés s’éloigne de l’intervalles de justesse
avec une valeur d’erreur de l’ordre de 10^-1 qui est valeur importante
lorsque la valeur de SNR diminue le signal de sortie réelle et celle bruité s’éloigne et le temps de
convergence de gain devient très important.
Donc en présence de bruit cette méthode n’est pas bonne.
Conclusion :
L’identification avec méthode de moindre carré récursif et non récursif ne peuvent pas donner des
valeurs estimés proche de valeur réelle si le système est bruité .donc on a besoin d’une autre
méthode plus général c’est la méthode de moindre carré récursif généralisé.
2EME ÉLECTRIQUE G1 TP2 8
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
2éme cas : e(k) est un bruit coloré
clc; clear all; close all;
% premiere parie :intialisation
y2=0;P=1000*eye(2);teta_mcr=[0 ;0];y2F=0;Pg=1000*eye(2);
teta_mcg=[0;0]; Pe=1000;Ke=[];f1=0;e2=0;b2=0 ;
rand('seed',0);
u=idinput(100,'prbs',[0 0.05],[0 1]);
for k=2:length(u)
%deusieme partie: processus
b1=b2;
b2=0.1*(-0.5+rand);
y1=y2;
y2=0.2*y1+2*u(k)+1*b1+10*b2;
h=[y1 u(k)]';
%troisième partie :moindres carre recursifs
P=P-(P*h*h'*P)/(1+h'*P*h);
K=P*h;
teta_mcr=teta_mcr+K*(y2-h'*teta_mcr);
%Quatrième partie:moindres carres généralisés
uF=u(k)+f1*u(k-1);
y1F=y2F;
y2F=y2+f1*y1;
hg=[y1F uF]';
Pg=Pg-(Pg*hg*hg'*Pg)/(1+hg'*Pg*hg);
Kg=Pg*hg;
teta_mcg=teta_mcg+Kg*(y2F-hg'*teta_mcg);
e1=e2;
e2=y2-h'*teta_mcg;
he=e1;
Pe=Pe-(Pe*he*he'*Pe)/(1+he'*Pe*he);
Ke=Pe*he;
f1=f1+Ke*(e2-he'*f1);
end;
%cinquième partie :comparaison des résultats
sum((teta_mcr-[0.2;2]).^2)
sum((teta_mcg-[0.2;2]).^2)
On obtient les valeurs suivantes:
a1=
0.0149
b1 =
0.0087
Conclusion :
2EME ÉLECTRIQUE G1 TP2 9
Université de Monastir
Ecole Nationale d’Ingénieur de Monastir
On remarque que cette méthode est plus efficace que les autres méthodes.
2EME ÉLECTRIQUE G1 TP2 10