Akademia Górniczo-Hutnicza

im. Stanisława Staszica

Wydział Inżynierii Mechanicznej i Robotyki

Analiza Sygnałów i Identyfikacja Procesów

Projekt –zadanie 11

Wykonał:

Przemysław Basiaga Gr17

Zadany układ

gdzie: R1=200Ω

R2=33,3 Ω

t_p=0.0001s

1. Identyfikacja obiektu G2 na podstawie odpowiedzi skokowej układu.

Badany układ jest obwodem RLC wiec jego transmitancja będzie miała postać:

Do obliczeń wykorzystamy jednak trochę inna postać, mianowicie:

Rys.1 Odpowiedź skokowa układu

Z wykresu odczytujemy:

y1=0.3511 y2=0.123286 Tosc=0.0422s

Wzory użyte do obliczenia poszczególnych parametrów :

ω=2π/ Tosc

Po wyliczeniu otrzymujemy:

Transmitancja przyjmuje postać:

Wyliczone parametry układu:

Kod z programu MatLab:

[y1,t1]=max(s);%wyznaczenie max y2

[y2,Ti]=min(s(t1:end)); %wyznaczenie min y2

y1=y1-1

y2=1-y2

l=y1/y2;

tp=0.0001; %czas próbkowania

T=Ti*tp*2 %wyznaczenie stałej czasowej

omega=(2*pi)/T %pulsacja drgań tłumionych

tlu=log(l)/(sqrt(pi^2+(log(l))^2)) %tłumienie

wstlum=(omega*tlu)/(sqrt(1-tlu^2)) %współczynnik tłumienia

omega0=sqrt(omega^2+wstlum^2) %pulsacja drgań własnych

T0=1/omega0

L1=[1];

M1=[T0^2 2*tlu*T0 1];

G=tf(L1,M1)

[f]=step(G,t);

figure(1)

plot(t,s,'r');%odpowiedź skokowa obiektu toru zakłóceń

hold on

plot(t,x,'g')

hold on

plot(t,f,'b')%odpowiedź skokowa obiektu zidentyfikowanego

grid;

C2=(2*tlu*T0)/Rz % pojemność C2

L2=(T0^2)/C2%indykcyjność L

Rys.2 Porównanie odpowiedzi układu zadanego i zidentyfikowanego

1. Identyfikacja obiektu G2 na podstawie char. amplitudowo-częstotliwościowej

Po poddaniu odpowiedzi skokowej transformacji Fouriera i wygładzeniu otrzymujemy wykres:

Rys.3 Charakterystyka amplitudowo-częstotliwościowa układu

Wyznaczamy częstotliwości dla których cześć urojona równa się rzeczywistej oraz dla których cześć rzeczywista się zeruje(wykres przecina się z osią urojoną)

Otrzymujemy:

Z poniższych wzorów wyliczamy:

Cała transmitancja ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

z=zeros(1,18400);

i(1601:20000)=z;

y=fft(i);%przekształcenie Fouriera

figure(1)

plot(y)% wykres nyquista

hold on

plot((0:0.0001:1.2),(0:-0.0001:-1.2),'g'); %prosta pod katem -45 stopni%

grid;

Ns=20000;

df=1/(tp*Ns);

w1=37*df*2*pi; %omega1

w2=50*df*2*pi; %omega2

tlu1=0.5*(w2/w1)*(1-(w1^2/w2^2));%tłumienie

T01=1/w2;

L1=[1];

M1=[T01^2 2*tlu1*T01 1];

G1=tf(L1,M1)

[f1]=impulse(G1,t);

figure(2)

plot(t,i(1:1601))%odpowiedź impulsowa obiektu toru zakłóceń

hold on

plot(t,f1*tp,'r')%odpowiedź impulsowa obiektu toru zakłóceń zidentyfikowanego

grid;

wstlum2=(w1*tlu1)/(sqrt(1-tlu1^2)) ; % współczynnik tłumienia

C2=(2*tlu1*T01)/Rz % pojemność C2

L2=(T01^2)/C2%indukcyjność L

Rys.4 Porównanie odpowiedzi układu zadanego i zidentyfikowanego w pkt. 2

1. Identyfikacja obiektu G2 na podstawie gęstości widmowej mocy

Na wejście podany został sygnał z.

Po poddaniu wejścia i wyjścia układu transformacji Fouriera i wygładzeniu otrzymujemy wykres:

Rys.5 Charakterystyka amplitudowo-częstotliwościowa układu (po uwzględnieniu wejścia i wyjścia)

Podobnie jak w punkcie 2 wyznaczamy częstotliwości

Dalej otrzymujemy:

Cała transmitancja ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

wz=fft(z);%transformata Fouriera wejścia z

wy=fft(y);%transformata Fouriera wyjścia y

Szz=wz.*conj(wz);

Szy=conj(wz).*wy;

gw3=Szy./Szz; %pierwszy sposób uśrednienia

gws31=decimate(gw3,20);

Szz_1=decimate(Szz,5); %drugi sposób uśredniania

Szy_1=decimate(Szy,5);

gws3=Szy_1./Szz_1;

figure(3)

plot(gws3(1:50),'r') %wykres nyquista z drugiego sposobu

hold on;

plot((0:0.0001:1.5),(0:-0.0001:-1.5),'g');

grid;

title('2 sp')

figure(4)

plot(gws31(1:15),'b') %wykres nyquista z pierwszego sposobu

hold on;

plot((0:0.0001:1.5),(0:-0.0001:-1.5),'g');

title('1 sp')

grid;

%wyliczenie omegi1 i omegi2 z wykresu nyquista z drugiego sposobu

Ns=6001;%ilość próbek

df=1/(tp*Ns);

w3=11*df*2*pi%omega1

w4=15*df*2*pi%omega2

tlu3=0.5*(w4/w3)*(1-(w3^2/w4^2))%tłumienie

T03=1/w4

L3=[1];

M3=[T03^2 2*tlu3*T03 1];

G3=tf(L3,M3)

[f3]=step(G3,t);

figure(5)

plot(t,s)%odpowiedź układu toru zakłóceń

hold on

plot(t,f3,'r')%odpowiedź układu toru zakłóceń zidentyfikowanego

grid;

Rz=33,3;

L1=T03*Rz/(2*tlu3) % indukcyjność L przy Rz=33.3

C1=2*tlu3*T03/Rz % pojemność C

Rys.6 Porównanie odpowiedzi układu zadanego i zidentyfikowanego w pkt. 3

4. Identyfikacja układu(G1, G2) z zakłóceniem (metoda korelacji)

Postępujemy podobnie ja powyżej i otrzymujemy wykres:

Rys.7 Charakterystyka amplitudowo-częstotliwościowa układu G1

Obliczając transmitancje G1 wystarczy policzyć To gdyż jest to układ inercyjny I rzędu

Cała transmitancja ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

wz=fft(z,1000000);

wy=fft(y,1000000);

wu=fft(u,1000000);

Sx1x1=wu.*conj(wu);

Sx2x2=wz.*conj(wz);

Sx1y=conj(wu).*wy;

Sx2y=conj(wz).*wy;

Sx1x2=conj(wu).*wz;

Sx2x1=conj(wz).*wu;

f=110;

Sx1x1=decimate(Sx1x1,f);

Sx2x2=decimate(Sx2x2,f);

Sx1y=decimate(Sx1y,f);

Sx2y=decimate(Sx2y,f);

Sx1x2=decimate(Sx1x2,f);

Sx2x1=decimate(Sx2x1,f);

G1ii=(Sx2x2.*Sx1y-Sx1x2.*Sx2y)./(Sx1x1.*Sx2x2-Sx1x2.*Sx2x1); %transmitancje widmowe tu szukamy próbek dla obiektu

figure(3)

plot(G1ii(1:1200),'r')%wykres częstotliwościowy układu toru głównego

hold on;

plot((0:0.0001:1),(0:-0.0001:-1),'g'); %prosta pod kątem -45stopnie

grid;

xlabel('Re(G(jw))')

ylabel('Img(G(jw))')

Ns=9091;%ilość próbek

df=1/(tp*Ns);

w5=7*df*2*pi%omega1

T05=1/w5

L5=[1];

M5=[T05 1];

G5=tf(L5,M5)

[f5]=step(G5,t);

figure(6)

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego

grid;

Ru=200;

C4=T05/Ru %pojemność C1

Rys.8 Odpowiedź skokowa układu obiektu zidentyfikowanego w pkt. 4

Postępujemy podobnie ja powyżej i otrzymujemy wykres:

Rys.9 Charakterystyka amplitudowo-częstotliwościowa układu G2

Wyznaczamy:

Dalej otrzymujemy:

Cała transmitancja toru zakłócenia ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

wz=fft(z,1000000);

wy=fft(y,1000000);

wu=fft(u,1000000);

Sx1x1=wu.*conj(wu);

Sx2x2=wz.*conj(wz);

Sx1y=conj(wu).*wy;

Sx2y=conj(wz).*wy;

Sx1x2=conj(wu).*wz;

Sx2x1=conj(wz).*wu;

f=110;

Sx1x1=decimate(Sx1x1,f);

Sx2x2=decimate(Sx2x2,f);

Sx1y=decimate(Sx1y,f);

Sx2y=decimate(Sx2y,f);

Sx1x2=decimate(Sx1x2,f);

Sx2x1=decimate(Sx2x1,f);

G2ii=(Sx1x1.*Sx2y-Sx2x1.*Sx1y)./(Sx1x1.*Sx2x2-Sx1x2.*Sx2x1); % tu dla zakłóceń

plot(G2ii(1:800),'b')%wykres częstotliwościowy układu toru zakłóceń

hold on;

plot((0:0.0001:1.2),(0:-0.0001:-1.2),'g'); %prosta pod kątem -45stopnie

grid;

xlabel('Re(G(jw))')

ylabel('Img(G(jw))')

Ns=9091;%ilość próbek

df=1/(tp*Ns);

w3=17*df*2*pi;%omega 1

w4=23*df*2*pi;%omega 2

tlu4=0.5*(w4/w3)*(1-(w3^2/w4^2));%tłumienie

T04=1/w4;

L4=[1];

M4=[T04^2 2*tlu4*T04 1];

G4=tf(L4,M4);

[f4]=step(G4,t);

figure(5)

plot(t,s)%odpowiedź skokowa obiektu toru zakłóceń

hold on

plot(t,f4,'r')%odpowiedź skokowa obiektu toru zakłóceń zidentyfikowanego w wpk.4

grid;

Rz=33,3;

L4=T04*Rz/(2*tlu4); % indukcyjność L przy Rz=33.3

C4=2*tlu4*T04/Rz; % pojemność C2

Rys.10 Porównanie odpowiedzi układu toru zakłóceń i zidentyfikowanego w pkt. 4

5.Identyfikacja układu(G1, G2) z zakłóceniem (identyfikacja metodami parametrycznymi)

1. Dla układu SISO

Transmitancja G1otrzymana metodą ‘arx’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Sarx=d2c(arx221);

[l,m]=th2tf(Sarx);

Sarx=tf(l/m(:,3),m/m(:,3))

[f]=step(Sarx,t);

figure(1)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'arx' w uk.SISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Sarx-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.11Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘arx’ w uk SISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1otrzymana metodą ‘amx’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Samx=d2c(amx2221);

[l,m]=th2tf(Samx);

Samx=tf(l/m(:,3),m/m(:,3))

[f]=step(Samx,t);

figure(2)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'amx'w uk.SISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Samx-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.12 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘amx’ w uk SISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1otrzymana metodą ‘iv’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Siv=d2c(iv331);

[l,m]=th2tf(Siv);

Siv=tf(l/m(:,4),m/m(:,4))

[f]=step(Siv,t);

figure(3)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'iv' w uk.SISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Siv-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.13 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘iv’ w uk SISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1otrzymana metodą’ oe’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Soe=d2c(oe221);

[l,m]=th2tf(Soe);

Soe=tf(l/m(:,3),m/m(:,3))

[f]=step(Soe,t);

figure(4)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'oe' w uk.SISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Soe-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.14Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘oe’ w uk SISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1otrzymana metodą ‘bj’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Sbj=d2c(bj22221);

[l,m]=th2tf(Sbj);

Sbj=tf(l/m(:,3),m/m(:,3))

[f]=step(Sbj,t);

figure(5)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'bj' w uk.SISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Sbj-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.15 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘bj’ w uk SISO i układu toru głównego zidentyfikowanego w pkt. 4

1. Dla układu MISO

Transmitancja G1otrzymana metodą ‘arx’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Marx=d2c(arx2211);

[l,m]=th2tf(Marx,1);

Marx=tf(l/m(:,3),m/m(:,3))

[f]=step(Marx,t);

figure(6)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'arx' w uk.MISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Marx-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.16 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘arx’ w uk MISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1otrzymana metodą’ amx’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Mamx=d2c(amx22211);

[l,m]=th2tf(Mamx,1);

Mamx=tf(l/m(:,3),m/m(:,3))

[f]=step(Mamx,t);

figure(7)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'amx'w uk.MISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Mamx-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.17 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘amx’ w uk MISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1 otrzymana metodą ‘iv’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Miv=d2c(iv2211);

[l,m]=th2tf(Miv,1);

Miv=tf(l/m(:,3),m/m(:,3))

[f]=step(Miv,t);

figure(8)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'iv' w uk.MISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Miv-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.18 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘iv’ w uk MISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1 otrzymana metodą ‘oe’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Moe=d2c(oe2211);

[l,m]=th2tf(Moe,1);

Moe=tf(l/m(:,3),m/m(:,3))

[f]=step(Moe,t);

figure(9)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'oe' w uk.MISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Moe-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.19 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘oe’ w uk MISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G1 otrzymana metodą ‘bj’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Mbj=d2c(bj222211);

[l,m]=th2tf(Mbj,1);

Mbj=tf(l/m(:,3),m/m(:,3))

[f]=step(Mbj,t);

figure(10)

plot(t,f)%odpowiedź skokowa obiektu toru głównego zidentyfikowanego metoda 'bj' w uk.MISO

hold on;

plot(t,f5,'r')%odpowiedź skokowa obiektu toru głównego zidentyfikowanego w pkc.4

grid;

title('Mbj-G1')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.20 Porównanie odpowiedzi skokowej układu toru głównego zidentyfikowanego metoda ‘bj’ w uk MISO i układu toru głównego zidentyfikowanego w pkt. 4

Transmitancja G2 otrzymana metodą ‘arx’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Marx=d2c(arx2211);

[l,m]=th2tf(Marx,2);

Marx=tf(l/m(:,3),m/m(:,3))

[f]=step(Marx,t);

figure(11)

plot(t,f)%odpowiedź skokowa obiektu toru zakłóceń zidentyfikowanego metoda 'arx' w uk. MISO

hold on;

plot(t,s,'r')%odpowiedź skokowa obiektu toru zakłóceń

grid;

title('Marx-G2')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.21 Porównanie odpowiedzi skokowej układu toru zakłóceń zidentyfikowanego metodą ‘arx’ w uk MISO i układu toru zakłóceń

Transmitancja G2 otrzymana metodą ‘amx’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Mamx=d2c(amx22211);

[l,m]=th2tf(Mamx,2);

Mamx=tf(l/m(:,3),m/m(:,3))

[f]=step(Mamx,t);

figure(12)

plot(t,f)%odpowiedź skokowa obiektu toru zakłóceń zidentyfikowanego metoda 'amx' w uk. MISO

hold on;

plot(t,s,'r')%odpowiedź skokowa obiektu toru zakłóceń

grid;

title('Mamx-G2')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.22 Porównanie odpowiedzi skokowej układu toru zakłóceń zidentyfikowanego metodą ‘amx’ w uk MISO i układu toru zakłóceń

Transmitancja G2 otrzymana metodą ‘iv’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Miv=d2c(iv2211);

[l,m]=th2tf(Miv,2);

Miv=tf(l/m(:,3),m/m(:,3))

[f]=step(Miv,t);

figure(13)

plot(t,f)%odpowiedź skokowa obiektu toru zakłóceń zidentyfikowanego metoda 'iv' w uk. MISO

hold on;

plot(t,s,'r')%odpowiedź skokowa obiektu toru zakłóceń

grid;

title('Miv-G2')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.23 Porównanie odpowiedzi skokowej układu toru zakłóceń zidentyfikowanego metodą ‘iv’ w uk MISO i układu toru zakłóceń

Transmitancja G2 otrzymana metodą ‘oe’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Moe=d2c(oe2211);

[l,m]=th2tf(Moe,2);

Moe=tf(l/m(:,3),m/m(:,3))

[f]=step(Moe,t);

figure(14)

plot(t,f)%odpowiedź skokowa obiektu toru zakłóceń zidentyfikowanego metoda 'oe' w uk. MISO

hold on;

plot(t,s,'r')%odpowiedź skokowa obiektu toru zakłóceń

grid;

title('Moe-G2')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.24 Porównanie odpowiedzi skokowej układu toru zakłóceń zidentyfikowanego metodą ‘oe’ w uk MISO i układu toru zakłóceń

Transmitancja G2 otrzymana metodą ‘bj’ ma postać:

Wyliczone parametry układu:

Kod z programu MatLab:

Mbj=d2c(bj222211);

[l,m]=th2tf(Mbj,2);

Mbj=tf(l/m(:,3),m/m(:,3))

[f]=step(Mbj,t);

figure(15)

plot(t,f)%odpowiedź skokowa obiektu toru zakłóceń zidentyfikowanego metoda 'bj' w uk. MISO

hold on;

plot(t,s,'r')%odpowiedź skokowa obiektu toru zakłóceń

grid;

title('Mbj-G2')

xlabel('Time(sec)')

ylabel('Amplitude')

Rys.25 Porównanie odpowiedzi skokowej układu toru zakłóceń zidentyfikowanego metodą ‘bj’ w uk MISO i układu toru zakłóceń

Porównanie wyników otrzymanych w punkcie V

`

	Układ	Metoda	Transmitancja otrzymana w punkcie V	T0		L	C
1.	SISO-G1	ARX			-	-
2.	SISO-G1	AMX			--	-
3.	SISO-G1	IV			-	-
4.	SISO-G1	OE			-	--
5.	SISO-G1	BJ			-	-
6.	MISO-G1	ARX			-	-
7.	MISO -G1	AMX			-	-
8.	MISO-G1	IV			-	-
9.	MISO-G1	OE			-	-
10.	MISO-G1	BJ				-
11.	MISO-G2	ARX
12.	MISO-G2	AMX
13.	MISO-G2	IV
14.	MISO-G2	OE
15.	MISO-G2	BJ