AKADEMIA GÓRNICZO-HUTNICZA

im. Stanisława Staszica Wydział Inżynierii Mechanicznej i Robotyki

Katedra Automatyzacji Procesów

Identyfikacja Procesów Technologicznych.

Marcin Baran

Rok studiów IV Semestr 7

Kierunek studiów: Automatyka i Robotyka

Specjalność: Automatyka i Metrologia

Projekt Zaliczeniowy.

Konsultował: dr inż. Andrzej Podsiadło.

Ocena: ......................................................................

Kraków, 2008

Schemat zadanego układu:

DANE:

K1 = 10 000 [Nm/s]

K2 = 0,005 [Nm/s]

SZUKANE:

m1 = ?

m2 = ?

D1 = ?

D2 = ?

Układ rozpatrujemy jako dwa oddzielne układy w postaci:

Mamy do wykonania:

1. Wykresy z odpowiedzią A1(t), wyznaczyć m i D

2. Wyznaczyć S_xx S_xy S_yy S_yz S_xz

3. Identyfikacja (gęstość widmowa).

1) Metoda 1

a)układ 1

y1=11.38;

y2=3.468;

te1=1046;

te2=1092;

k1=200000

Dt=(te2-te1)*0.005

T=2*Dt

w=2*pi/T

ksi=log(y1/y2)/sqrt(pi^2+(log(y1/y2))^2)

w0=sqrt(w^2-ksi^2)

T0=1/w0

M1=T0^2*k1

D1=2*ksi*T0*k1

M1=1072 kg

D1=10364

b) układ 2

k2=10000;

N=50000;

y0=y(1:5000);

y0(5001:N)=0;

z0=z(1:5000);

z0(5001:N)=0;

Sy0y0=fft(y0).*conj(fft(y0));

Sy0z0=conj(fft(y0)).*fft(z0);

Sy0y0=Sy0y0(1:1000);

Sy0z0=Sy0z0(1:1000);

R=25;

Sy0y0=decimate(Sy0y0,R);

Sy0z0=decimate(Sy0z0,R);

g1i=Sy0z0./Sy0y0;

figure(1)

plot(g1i);grid

deltaf=1/(0.005*N)*R;

f1=deltaf*9;

w1=2*pi*f1;

f2=deltaf*13;

w2=2*pi*f2;

T0=1/w2;

M2 = T0^2*k2

ksi=(1/2*w2/w1*(1-w1^2/w2^2));

D2=2*ksi*k2*T0

M2 =149.8834

D2 =920.8176

2) Metoda 2

a) układ 1

N=50000;

xi=x1(1:10000);

xi(10001:N)=0;

yi=y1(1:10000);

yi(10001:N)=0;

Sxixi=fft(xi).*conj(fft(xi));

Sxiyi=conj(fft(xi)).*fft(yi);

Sxixi=Sxixi(1:1150);

Sxiyi=Sxiyi(1:1150);

R=25;

Sxixi=decimate(Sxixi,R);

Sxiyi=decimate(Sxiyi,R);

g1i=Sxiyi./Sxixi;

figure(1)

plot(g1i);grid

deltaf=1/(0.005*N)*R;

f1=deltaf*13;

w1=2*pi*f1;

f2=deltaf*18;

w2=2*pi*f2;

T0=1/w2;

M1 = T0^2*k1

ksi=(1/2 * w2/w1 * (1 - w1^2/w2^2));

D1=2*ksi*k1*T0

M1=1563 kg

D1= 11714

b) układ 2

N=50000;

we=y1(1:10000);

we(10001:N)=0;

wy=z1(1:10000);

wy(10001:N)=0;

Swewe=fft(we).*conj(fft(we));

Swewy=conj(fft(we)).*fft(wy);

Swewe=Swewe(1:1150);

Swewy=Swewy(1:1150);

R=25;

Swewe=decimate(Swewe,R);

Swewy=decimate(Swewy,R);

g1i=Swewy./Swewe;

figure(1)

plot(g1i);grid

deltaf=1/(0.005*N)*R;

f1=deltaf*9;

w1=2*pi*f1;

f2=deltaf*13;

w2=2*pi*f2;

T0=1/w2;

M2 = T0^2*k2

ksi=(1/2 * w2/w1 * (1 - w1^2/w2^2));

D2=2*ksi*k2*T0

M2=149.88 kg

D2=920

3) Metoda 3

a) uklad 1

Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t)

A(q) = 1 - 1.97 q^-1 + 0.9753 q^-2

B(q) = 0.002478 q^-1 + 0.002458 q^-2

Estimated using ARX with focus from data set m1

Loss function 1.10619e-022 and FPE 1.10707e-022

Sampling interval: 0.005

b =

0 0.0025 0.0025

a =

1.0000 -1.9704 0.9753

Transfer function:

0.002478 z + 0.002458

---------------------

z^2 - 1.97 z + 0.9753

Sampling time: 0.005

Transfer function:

-1.837e-006 s + 200

-------------------

s^2 + 9 s + 200

M1= 1000 kg

D1=9000

b) układ 2

Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t)

A(q) = 1 - 1.949 q^-1 + 0.9512 q^-2

B(q) = 0.002458 q^-1 - 2.065e-005 q^-2

Estimated using ARX with focus from data set m2

Loss function 1.91823e-008 and FPE 1.91977e-008

Sampling interval: 0.005

b =

0 0.0025 -0.0000

a =

1.0000 -1.9488 0.9512

Transfer function:

0.002458 z - 2.065e-005

-----------------------

z^2 - 1.949 z + 0.9512

Sampling time: 0.005

Transfer function:

0.2521 s + 99.96

---------------------

s^2 + 9.998 s + 99.97

M2=100.04 kg

D2=1000

Tabela:

	1 metoda	2 metoda	3 metoda
M1	1072	1563	1000
D1	10364	11714	9000
M2	149	149.88	100.04
D2	920	920	1000

Wyniki obserwacji

Wyniki jakie otrzymaliśmy w trzech metodach są bardzo zbliżone do siebie.

Różnice w wartościach wyników są spowodowane przybliżeniami jakie zostały wprowadzone w pierwszej metodzie i małą liczbą próbek jakie posiadaliśmy w drugiej metodzie.

Najdokladniesze wyniki wyszły w metodzie parametrycznej .

---