Na tej strani lahko vidite del protokola za laboratorijske vaje pri predmetu Identifikacije. Prikazano je besedilo naloge, potek naloge in odziv na eno vzbujanje. Prav tako je tukaj tudi program v programskem jeziku Matlab, s katerim lahko pridemo do odzivov pri drugih vzbujanjih.

Besedilo naloge:

Imamo proces 1., 2. in 3. reda. S pomočjo osebnega računalnika in A/D-D/A vmesnika posnemite odzive na različne testne signale. Nato s temi odzivi in testnimi signali identificirajte posamezne procese.

Parametri procesa:

• proces 1. reda: T₁=30s
• proces 2. reda: T₁=30s T₂=9s
• proces 3. reda: T₁=30s T₂=9s T₃=9s

V laboratoriju smo snemali odziv sistema, kateremu smo lahko spreminjali časovne konstante T₁, T₂ in T₃. Najprej smo posneli sistem kot da je sistem prvega reda, zato smo časovni konstanti T₂ in T₃ nastavili na nič. Sledilo je, da smo nastavili tudi časovno konstanto T₂ in tako dobi časovni odziv sistema drugega reda. Nato smo dodali še časovno konstanto T₃ in dobili sistem tretjega reda. Pri samem snemanju so se nam izmerjene vrednosti (vhod in izhod) vpisovale v datoteke.

Pri vsakem redu sistema smo snemali odzive na vzbujanje s stopnico in na vzbujanje s PRBS signalom. Tako smo izvedli šest meritev in dobili šest ločenih datotek. 

Po meritvi smo morali izmerjene sisteme tudi identificirati. To smo opravljali s programskim paketom Matlab. 

Razlaga programa v Matlabu:

Ker smo opravili merjenja na treh različnih sistemih z dvema različnima signaloma, sem po definiranju posameznih originalnih sistemov, za kasnejšo primerjavo, v programu dal možnost izbire kateri sistem želimo identificirati in s kakšnim signalom. Sama identifikacija je pri vsakem sistemu in signalu potekala enako.

Za začetek identifikacije smo prebrali izmerjene podatke iz datoteke vhodno-izhodnih vrednosti za tisti primer sistema in signala, katerega smo se odločili identificirati. Te prebrane vrednosti smo zapisali v matriko Z. Matrika je v prvem stolpcu vsebovala izhodne vrednosti, v drugem pa vhodne vrednosti. Stolpca sem posebej zapisal v vektorja Y (izhod) in U (vhod) in jih tudi izrisal.

Prava identifikacija sistema se je začela s klicem funkcije RARMAX (rekurzivni armax model). V samo funkcijo RARMAX smo morali vnesti vhod, red, zakasnitev armax modela, kateri algoritem faktorja pozabljanja bomo uporabili in sam faktor pozabljanja. Funkcija nam je kot rezultat v matriko THM vrnila ocenjene parametre a, b in c, ter v YHAT odziv modela sistema. Pri tem so morali parametri a in b v času konvergirati, zato smo jih izrisali. Iz poteka teh parametrov je bilo opaziti da proti koncu meritve ti parametri že konvergirajo, da pa smo napako zaradi kratkega časa meritve držali čim nižjo, smo kot parametre vzeli povprečje zadnjih petih vrednosti, kjer so bili že dovolj blizu končne vrednosti.

Identifikacija se je s tem končala. Ostalo nam je le še preveritev, kako dober (če je sploh dober) je naš dobljeni model. Zato smo morali preveriti ojačanje (statična verifikacija), primerjati parametre a in b in opazovati razlike, ter primerjati izhod sistema Y in izhod modela Y_MODELA (dinamična verifikacija).

• Statično verifikacijo smo računali po formuli: 

• Pri primerjanju parametrov poznanega sistema in parametrov modela smo opazili, da do odstopanj pride, vendar so le ta odstopanja majhna. To se pojavi zaradi diskretnega izhoda, ki je vzrok da dobimo lahko le boljšo ali slabšo aproksimacijo z nekim polinomom.

• Pri dinamični verifikaciji smo model našega sistema vzbujali z enakim vhodom kot poprej poznan sistem in opazovali izhod. Iz razlike med odzivoma obeh sistemov smo dobili občutek pogreška modela.

Identifikacijo sem opravil pri dejanskih izmerjenih rezultatih ter pri izmerjenih rezultatih, katerim srednje vrednosti sem nastavil na nič!

Sledijo posneti odzivi in program napisan v Matlabu:

Sistem 1.reda – vzbujanje s stopnico:

 

 

Dobljeni parametri:

Ojačanje K	0.9962
Parameter a1	-0.9775	Dejanski parameter a1	-0.9672	Razlika	0.0103
Parameter a2	0.0224	Dejanski parameter b1	0.0328	Razlika	0.0104

 

Program v jeziku Matlab:

% I D E N T I F I K A C I J E 

% Računalniške vaje 

% 

% Študent: Boštjan Štumberger Datum: 26.11.1997 

% Časovne konstante za mojo nalogo in vzorčni čas

% T1 = 0,5 min = 30 s

% T2 = 0,15 min = 9 s

% T3 = 0,15 min = 9 s

% Ts = 1 s

 

%-----------------------------------------------------------------

% Finkcija, v kateri izberemo kakšen sistem bomo identificirali

% (prvega, drugega ali tretjega reda, oz. i=1,2 ali 3

%-----------------------------------------------------------------

while(1) 

clc; 

disp(' IZBERI SISTEM, KI GA ŽELIŠ IDENTIFICIRATI:'); 

disp(' '); 

disp(' 1) sistem 1. reda'); 

disp(' 2) sistem 2. reda'); 

disp(' 3) sistem 3. reda'); 

disp(' '); 

disp(' 0) Izhod'); 

red=input('Vpisi željeno številko: ');

%-----------------------------------------------------------------

% Čas tipanja:

Ts=1;

%-----------------------------------------------------------------

% SISTEM 1. REDA % Zapis prenosne funkcije prvega reda

% 1

% G1(s) = ------

% (30s+1)

imen1=[30 1]; % Imenovalec funkcije prvega reda

%-----------------------------------------------------------------

% SISTEM 2. REDA % Zapis prenosne funkcije drugega reda

% 1

% G2(s) = --------------

% (30s+1)(9s+1)

imen2=conv(imen1,[9 1]); % Imenovalec funkcije drugega reda

%-----------------------------------------------------------------

% SISTEM 3. REDA

% 1

% G3(s) = ---------------------

% (30s+1)(9s+1)(9s+1)

imen3=conv(imen2,[9 1]); % Imenovalec funkcije tretjega reda

%-----------------------------------------------------------------

% Izpis sistema izbranega reda

% Imenovalcu izbranega sistema priredimo matriko 'imen'

%-----------------------------------------------------------------

clc; 

if (red==1)

imen=imen1;

disp('Prenosna funkcija sistema 1. reda: ');

elseif (red==2)

imen=imen2;

disp('Prenosna funkcija sistema 2. reda: ');

elseif (red==3)

imen=imen3;

disp('Prenosna funkcija sistema 3. reda: ');

elseif (red==0)

break;

end;

% Izbrano prenosno funkcijo izpišemo

disp(' ');

disp('Sistem v S prostoru:');

printsys(1,imen,'s');

%-----------------------------------------------------------------

% Z transformacija izbranega sistema z nictim redom

%-----------------------------------------------------------------

[std,imd]=c2dm(1,imen,Ts,'zoh');

disp(' ');

disp('Po Z transformaciji sistema dobimo: ');

% Izpišemo prenosno funkcijo v z prostoru

printsys(std,imd,'z');

disp(' ');

disp('Pritisni tipko za nadaljevanje');

pause;

clc

%-----------------------------------------------------------------

% Funkcija za izbiro vzbujalne funkcije

%-----------------------------------------------------------------

clc;

disp('IZBERI VZBUJANJE:');

disp(' ');

disp('1) stopnica');

disp('2) PRBS');

disp(' ');

disp('0) Izhod');

disp(' ');

% Spremenljivka vzbujanje zavzame izbiro

vzbujanje=input('Izberi vzbujanje: ');

%-----------------------------------------------------------------

% Čitanje matrike vhodno-izhodnih izmerjenih vrednosti

% Te vrednosti damo v matriko Z

%-----------------------------------------------------------------

if (red==1) % Če je sistem prvega reda

if (vzbujanje==1) % Kakšno vzbujanje smo izbrali

fid=fopen('1_red_s.dat','r');

Z=fscanf(fid,'%f%f',[2 inf]);

fclose(fid);

elseif (vzbujanje==2)

fid=fopen('1_red_p.dat','r');

Z=fscanf(fid,'%f%f',[2 inf]);

fclose(fid);

elseif (vzbujanje==0)

break;

end;

elseif (red==2) % Če je sistem drugega reda

if (vzbujanje==1) % Kakšno vzbujanje smo izbrali

fid=fopen('2_red_s.dat','r');

Z=fscanf(fid,'%f%f',[2 inf]);

fclose(fid);

elseif (vzbujanje==2)

fid=fopen('2_red_p.dat','r');

Z=fscanf(fid,'%f%f',[2 inf]);

fclose(fid);

elseif (vzbujanje==0)

break;

end;

elseif (red==3) % Če je sistem tretjega reda

if (vzbujanje==1) % Kakšno vzbujanje smo izbrali

fid=fopen('3_red_s.dat','r');

Z=fscanf(fid,'%f%f',[2 inf]);

fclose(fid);

elseif (vzbujanje==2)

fid=fopen('3_red_p.dat','r');

Z=fscanf(fid,'%f%f',[2 inf]);

fclose(fid);

elseif (vzbujanje==0)

break;

end;

end;

% Iz matrike prebranih vrednosti ločimo vhodne vrednosti od 

% izhodnih 

Z=Z';

U=Z(:,2); % Vzbujanje sistema - vhodni signal;

Y=Z(:,1); % Odziv sistema - izhodni signal;

%-----------------------------------------------------------------

% Vzbujanju //in odzivu// damo srednjo vrednost na nič

%-----------------------------------------------------------------

sestevek=0;

dolzina=0;

dolzina=length(Z);

zacetna=U(1);

if vzbujanje==1

disp('Vzbujanje = 1')

for i=1:dolzina

U(i)=U(i)-zacetna;

Y(i)=Y(i)-zacetna;

end;

end;

if vzbujanje==2

disp('Vzbujanje = 2')

for i=1:dolzina

sestevek=sestevek+U(i);

end;

srednja=sestevek/dolzina;

for i=1:dolzina

U(i)=U(i)-srednja;

% Y(i)=Y(i)-srednja;

end;

end;

zacetna=Y(1);

if zacetna<0

for i=1:dolzina

Y(i)=Y(i)+abs(zacetna);

end;

else

for i=1:dolzina

Y(i)=Y(i)-zacetna;

end;

end;

Z=[Y U];

%-----------------------------------------------------------------

% Izris vhodnih in izhodnih vrednosti sistema

figure(1);

hold off

plot(U,'b');

hold on

plot(Y,'r');

grid;

title('Vhod(modra) in izhod(rdeca) sistema');

disp(' ');

disp('Pritisni tipko za nadaljevanje');

pause;

clc

%-----------------------------------------------------------------

% Identifikacija sistema z funkcijo RARMAX

% Klicanje funkcije rarmax ki nam v matriko THM vpiše

% ocenjene parametre [a,b,c] in v YHAT izhod modela

% rarmax(vhod,redi in zakasnitve armax modela,algoritem faktorja 

% pozabljanja, faktor pozabljanja,začetne vrednosti parametrov,

% začetne vrednosti matrike P)

%-----------------------------------------------------------------

disp(' ');

% lambda=input(' Vpiši vrednost za lambdo (med 0,8 in 1) : '); 

alfa=input(' Vpiši vrednost za alfo (med 10 in 1e12) : '); 

lambda=1;

if (vzbujanje==10)

if (red==1)

alfa=1e9;

elseif (red==2)

alfa=1e9;

elseif (red==3)

alfa=1e12;

end;

end;

if (vzbujanje==20)

if (red==1)

alfa=500;

elseif (red==2)

alfa=0.1;

elseif (red==3)

alfa=10e9;

end;

end;

zak=1;

nn=[red red red zak];

spr=zeros(red*3);

TH0=spr(:,1); %Vpišem začetne vrednosti parametrov

P0 = alfa*eye(3*red); %Vpišem začetne vrednosti matrike P

[THM,YHAT] = rarmax(Z,nn,'ff',lambda,TH0,P0); 

disp(' ');

%-----------------------------------------------------------------

% Iz THM preberemo parametre a in b

%-----------------------------------------------------------------

if (red==1) % Ce je funkcija prvega reda

a1=THM(:,1); % a1 je v prvem stolpcu

b1=THM(:,2); % b1 je v drugem stolpcu

figure(2);

subplot(2,1,1),plot(a1,'b');title('parameter a_1');grid;

subplot(2,1,2),plot(b1,'b');title('parameter b_1');grid;

pause;

% Povprečimo parametre a in b (vzamem zadnjih deset vrednosti)

k1=size(a1); % Ugotovimo koliko je število vrednosti

k1=k1(:,1);

a1p=0;

for i=0:4 % Čitamo 5 zadnjih vrednosti

a1p=a1p+a1(k1-i);

end;

a1p=a1p/5;

l1=size(b1);

l1=l1(:,1);

b1p=0;

for i=0:4

b1p=b1p+b1(l1-i);

end;

b1p=b1p/5;

ap=[a1p];

bp=[b1p];

elseif (red==2) % Če je funkcija drugega reda

a1=THM(:,1); % Prvi stolpec od THM je a1 itd.

a2=THM(:,2);

b1=THM(:,3);

b2=THM(:,4);

figure(2);

subplot(2,2,1),plot(a1,'b');title('parameter a_1');grid;

subplot(2,2,3),plot(a2,'b');title('parameter a_2');grid;

subplot(2,2,2),plot(b1,'b');title('parameter b_1');grid;

subplot(2,2,4),plot(b2,'b');title('parameter b_2');grid;

pause;

% Povprečimo parametre a in b (vzamem zadnjih deset vrednosti)

k1=size(a1); % Ugotovimo koliko je število vrednosti

k1=k1(:,1);

a1p=0;

for i=0:4 % Čitamo 10 zadnjih vrednosti

a1p=a1p+a1(k1-i);

end;

a1p=a1p/5;

l1=size(b1);

l1=l1(:,1);

b1p=0;

for i=0:4

b1p=b1p+b1(l1-i);

end;

b1p=b1p/5;

k2=size(a2);

k2=k2(:,1);

a2p=0;

for i=0:4

a2p=a2p+a2(k2-i);

end;

a2p=a2p/5;

l2=size(b2);

l2=l2(:,1);

b2p=0;

for i=0:4

b2p=b2p+b2(l2-i);

end;

b2p=b2p/5;

ap=[a1p;a2p];

bp=[b1p;b2p];

elseif (red==3) % Če je funkcija trtjega reda

a1=THM(:,1); % Prvi stolpec od THM je a1

a2=THM(:,2);

a3=THM(:,3);

b1=THM(:,4);

b2=THM(:,5);

b3=THM(:,6);

figure(2);

subplot(3,2,1),plot(a1,'b');title('parameter a_1');grid;

subplot(3,2,3),plot(a2,'b');title('parameter a_2');grid;

subplot(3,2,5),plot(a3,'b');title('parameter a_3');grid;

subplot(3,2,2),plot(b1,'b');title('parameter b_1');grid;

subplot(3,2,4),plot(b2,'b');title('parameter b_2');grid;

subplot(3,2,6),plot(b3,'b');title('parameter b_3');grid;

pause;

% Povprečimo parametre a in b (vzamemo zadnjih deset vrednosti)

k1=size(a1);

k1=k1(:,1);

a1p=0;

for i=0:4

a1p=a1p+a1(k1-i);

end;

a1p=a1p/5;

l1=size(b1);

l1=l1(:,1);

b1p=0;

for i=0:4

b1p=b1p+b1(l1-i);

end;

b1p=b1p/5;

k2=size(a2);

k2=k2(:,1);

a2p=0;

for i=0:4

a2p=a2p+a2(k2-i);

end;

a2p=a2p/5;

l2=size(b2);

l2=l2(:,1);

b2p=0;

for i=0:4

b2p=b2p+b2(l2-i);

end;

b2p=b2p/5;

k3=size(a3);

k3=k3(:,1);

a3p=0;

for i=0:4

a3p=a3p+a3(k3-i);

end;

a3p=a3p/5;

l3=size(b3);

l3=l3(:,1);

b3p=0;

for i=0:4

b3p=b3p+b3(l3-i);

end;

b3p=b3p/5;

ap=[a1p;a2p;a3p];

bp=[b1p;b2p;b3p];

end

st_m=[0 bp'];

im_m=[1 ap'];

%-----------------------------------------------------------------

% Razlika med Y in YHAT ( dinamična verifikacija )

%-----------------------------------------------------------------

e1=Y-YHAT;

figure(3);

subplot(2,1,1);

plot(Y,'r'),title('Y(rdeca) in Y_H_A_T(modra)');

hold on,plot(YHAT,'b'),grid,hold off;

subplot(2,1,2),plot(e1,'r'); grid;

title('Razlika med Y in Y_H_A_T');

pause;

%-----------------------------------------------------------------

% Primerjava parametrov a in b funkcije H(z) in matrike THM

%-----------------------------------------------------------------

clc

disp('Razlika parametrov med pravim sistemom in modelom')

disp(' ')

disp('Povprečne vrednosti parametrov b:')

disp(bp);

disp('Povprečne vrednosti parametrov a:')

disp(ap);

disp(' ');

disp('Pritisni tipko za nadaljevanje');

pause;

clc

%-----------------------------------------------------------------

% Preverjanje ojačanja K (statična verifikacija) 

% (mora biti okrog 1)

%-----------------------------------------------------------------

disp('Ojačanje:')

if (red==1)

K=b1p/(1+a1p)

elseif (red==2)

K=(b1p+b2p)/(1+a1p+a2p)

elseif (red==3)

K=(b1p+b2p+b3p)/(1+a1p+a2p+a3p)

end;

disp(' ');

disp('Pritisni tipko za nadaljevanje');

pause;

clc;

%-----------------------------------------------------------------

% Vzbujamo identificirani sistem

%-----------------------------------------------------------------

[ymodela3,xmodela]=dlsim(st_m,im_m,U);

%if vzbujanje==1

% izbira=45;

% else

% izbira=5;

%end;

%i=1;

%clear ymodela2;

%clear ymodela3;

%ymodela2(1)=-5;

%ymodela2(2)=-5;

%Zac_U(2)=-5;

%Zac_Y(2)=0;

% while ymodela2<Y(1)

% Zac_U(i)=izbira;

% Zac_Y(i)=Y(1);

% vmesna=i;

% [ymodela2,xmodela]=dlsim(st_m,im_m,Zac_U');

% disp(ymodela2(i))

% i=i+1;

% end;

% VektorU=[Zac_U';U];

% VektorY=[Zac_Y';Y];

% [ymodela2,xmodela]=dlsim(st_m,im_m,VektorU);

%d1=size(ymodela2);

%d2=size(Y);

%d1=d1(:,1);

%d2=d2(:,1);

%for b=1:d2-1

% ymodela3(b)=ymodela2(d1-d2+b+1);

%end;

% ymodela3(d2)=ymodela2(d1);

% ymodela3=ymodela3';

%-----------------------------------------------------------------

% Primerjava Y_modela in Y (Odzivi)

%-----------------------------------------------------------------

e1=Y-ymodela3; 

figure(4); 

clf

subplot(2,1,1); 

plot(Y,'b');hold on;title('Y(rdeča) in Y_modela(modra)'); 

plot(ymodela3,'r');grid;

subplot(2,1,2),plot(e1,'b'); grid; title('Pogrešek'); 

pause;

%-----------------------------------------------------------------

% end; %Prvi while(1)

end; %Drugi while(1)