POLITECHNIKA OPOLSKA

Wydział Elektrotechniki, Automatyki i Informatyki

Instytut Układów Elektromechanicznych i Elektroniki Przemysłowej

STUDIA DRUGIEGO STOPNIA STACJONARNE

AUTOMATYKA I ROBOTYKA

IDENTYFIKACJA PROCESÓW TECHNOLOGICZNYCH

PROJEKT

Proces identyfikacji obiektu automatyki z zastosowaniem pakietu przybornika System Identification Tool środowiska Matlab/Simulink

Prowadzący: dr hab. inż. Rafał Stanisławski, prof. PO	Pracę wykonał: Kamil KOLASA

Opole, czerwiec 2015

1. WPROWADZENIE

Celem projektu jest identyfikacja obiektu automatyki określonego jako element pakietu Simulink w postaci bloku funkcyjnego BlackBoxA. Wygląd elementy przedstawiono na rysunku 1.1.

Rys. 1.1. Wygląd elementu BlackBoxA w środowisku Simulink

Celem procesu identyfikacji jest uzyskanie transmitancji danego obiektu. W procesie tym wykorzystanie zostanie przybornik System Identification Tool będący częścią środowiska Matlab. Do głównych cech tego przybornika zaliczyć można:

• otrzymywanie transmitancji, modelu procesu lub przestrzeni stanu (ang. state-space) w oparciu o dane w dziedzinie czasu lub częstotliwości,

• otrzymywanie modeli: autoregresywnych (ARX, ARMAX), Box-Jenkins i Output-Error stosując techniki: najlepszego dopasowania, minimalizacji prognozowanego poziomu błędu (ang. PEM - Prediction-Error Minimalization) i identyfikacji podprzestrzeni systemu,

• estymacja parametrów modelu w trybie online,

• modelowanie czasowe (ang. Time-series modeling) (AR, ARMA) i prognozowanie,

• otrzymywanie modeli nieliniowych ARX i modeli Hammerstein-Wiener, z nieliniowościami wejść i wyjść takimi jak saturacja lub strefa martwa,

• liniowa i nieliniowa identyfikacja obiektów grey-box w celu estymacji modelu,

• estymacja opóźnienia, usuwanie dryftu głównego (ang. detrend), filtrowanie, zmiana częstotliwości próbkowania i rekonstrukcja brakujących danych.

2. PRZEBIEG PROCESU IDENTYFIKACJI

2.1. Analiza odpowiedzi czasowych obiektu BlackBoxA

Z wstępnej analizy wiadomo, że jest to obiekt dyskretny. W celu uzyskania informacji dotyczących zachowania obiektu na wejście układu podano skok jednostkowy (Step) oraz impuls Diraca (Delta Diraca). W ten sposób otrzymano odpowiedzi układu na wymuszenie reprezentowane przez zmienne u1, y1 oraz u2, y2. Rysunek 2.1 przedstawia układ realizujący to zadanie natomiast kod zawarty jest w skrypcie a1.m.

Rys. 2.1. Układ do analizy obiektu BlackBoxA w środowisku Simulink

2.2. Estymacja transmitancyjna - Transfer Function Models

W oparciu o dane wejściowe i wyjściowe (wektory wymuszenia oraz odpowiedzi układu) utworzone zostały modele transmitancyjne (Transfer Function Models) wykorzystując System Identification Tool. Obliczenia przeprowadzone zostały w dziedzinie czasu dyskretnego (Ts = 1). W oknie głównym narzędzia, przedstawionym na rysunku 2.2 widoczne są miniaturki opracowywanych charakterystyk w zakresie od tf1 do tf12. Charakterystyki w powiększeniu zostały pokazane na rysunkach 2.3, 2.4, 2.5 i 2.6. Cyfra zawierająca się w nazwie charakterystyki informuje o liczbie zer i biegunów danego modelu.

Rys. 2.2. Okno System Identification Tool z danymi estymowanymi modelem Transfer Function

Rys. 2.3. Charakterystyki: tf1, tf2, tf3	Rys. 2.4. Charakterystyki: tf4, tf5, tf6
Rys. 2.5. Charakterystyki: tf7, tf8, tf9	Rys. 2.6. Charakterystyki: tf10, tf11, tf12

Rysunek 2.3 przedstawia charakterystykę porównawczą odpowiedzi skokowych wszystkich otrzymanych modeli Transfer Function wraz z wskaźnikiem stopnia dopasowania do wzorca.

Rys. 2.7. Charakterystyka porównawcza wszystkich otrzymanych modeli Transfer Function wraz z współczynnikiem dopasowania

Biorąc pod uwagę kryterium sztywności modelu oraz optymalny poziom dopasowania wybrano model tf4. Rysunek 2.8 przedstawia transmitancję dyskretną wybranego modelu tf4 uzyskaną za pomocą skryptu a2.m.

Rys. 2.8. Transmitancja dyskretna uzyskana na bazie odpowiedzi skokowej

Odpowiedzi skokową i impulsową wybranego modelu oraz obiektu BlackBoxA przedstawiono na rysunkach 2.9 i 2.10, które generuje skrypt a3.m.

Rys. 2.9. Odpowiedzi skokowe układu o transmitancji uzyskanej z odpowiedzi skokowej oraz obiektu BlackBoxA	Rys. 2.10. Odpowiedzi impulsowe układu o transmitancji uzyskanej z odpowiedzi skokowej oraz obiektu BlackBoxA

2.3. Estymacja wielomianowa (Polynomial Models) - ARX

W oparciu o dane wejściowe i wyjściowe (wektory wymuszenia oraz odpowiedzi układu) utworzone zostały modele wielomianowe (ang. Polynomial Models) ARX (ang. AutoRegressive with eXogenous input) wykorzystując System Identification Tool. W oknie głównym narzędzia, przedstawionym na rysunku 2.11 widoczne są miniaturki opracowywanych charakterystyk w zakresie od arx111 do arx12121. Charakterystyki w powiększeniu zostały pokazane na rysunkach 2.12, 2.13, 2.14 i 2.15. Cyfry zawierające się w nazwie charakterystyki informuje o stopniach wielomianów: A(q) - na, B(q) + 1 - nb, opóźnienia wejścia-wyjścia - nk, w postaci ARX:[na nb nk].

Rys. 2.11. Okno System Identification Tool z danymi estymowanymi modelem ARX

Rys. 2.12. Charakterystyki: arx111, arx221, arx331	Rys. 2.13. Charakterystyki: arx441, arx551, arx661
Rys. 2.14. Charakterystyki: arx771, arx881, arx991	Rys. 2.15. Charakterystyki: arx10101, arx11111, arx12121

Rysunek 2.16 przedstawia charakterystykę porównawczą odpowiedzi skokowych wszystkich otrzymanych modeli ARX wraz z wskaźnikiem stopnia dopasowania do wzorca.

Rys. 2.16. Charakterystyka porównawcza wszystkich otrzymanych modeli ARX wraz z współczynnikiem dopasowania

Biorąc pod uwagę kryterium sztywności modelu oraz optymalny poziom dopasowania wybrano model arx10101 (na = 10, nb = 10, nk = 1). Rysunek 2.17 przedstawia transmitancję dyskretną wybranego modelu arx10101 otrzymaną poprzez zastosowanie funkcji idpoly() środowiska Matlab.

Rys. 2.17. Transmitancja obiektu wybranego modelu ARX (arx10101)

Odpowiedzi skokową i impulsową wybranego modelu arx10101 przedstawiono na rysunkach 2.18 i 2.19.

Rys. 2.18. Odpowiedź skokowa (arx10101)	Rys. 2.19. Odpowiedź impulsowa (arx10101)

2.4. Estymacja wielomianowa (Polynomial Models) - ARMAX

W oparciu o dane wejściowe i wyjściowe (wektory wymuszenia oraz odpowiedzi układu) utworzone zostały modele wielomianowe (ang. Polynomial Models) ARMAX (ang. AutoRegressive Moving Average with eXogenous input) wykorzystując System Identification Tool. W oknie głównym narzędzia, przedstawionym na rysunku 2.20 widoczne są miniaturki opracowywanych charakterystyk w zakresie od amx111 do amx12121. Charakterystyki w powiększeniu zostały pokazane na rysunkach 2.21, 2.22, 2.23 i 2.24. Cyfry zawierające się w nazwie charakterystyki informuje o stopniach wielomianów: A(q) - na, B(q) + 1 - nb, C(q) - nc, opóźnienia wejścia-wyjścia - nk, w postaci ARMAX:[na nb nc nk].

Rys. 2.20. Okno System Identification Tool z danymi estymowanymi modelem ARMAX

Rys. 2.21. Charakterystyki: amx1111, amx2221, amx3331	Rys. 2.22. Charakterystyki: amx4441, amx5551, amx6661
Rys. 2.23. Charakterystyki: amx7771, amx8881, amx9991	Rys. 2.24. Charakterystyki: amx1010101, amx1111111, amx1212121

Rysunek 2.25 przedstawia charakterystykę porównawczą odpowiedzi skokowych wszystkich otrzymanych modeli ARMAX wraz z wskaźnikiem stopnia dopasowania do wzorca.

Rys. 2.25. Charakterystyka porównawcza wszystkich otrzymanych modeli ARMAX wraz z współczynnikiem dopasowania

Biorąc pod uwagę kryterium sztywności modelu oraz optymalny poziom dopasowania wybrano model amx4441 (na = 4, nb = 4, nc = 4, nk = 1). Rysunek 2.26 przedstawia transmitancję dyskretną wybranego modelu amx4441 otrzymaną poprzez zastosowanie funkcji idpoly() środowiska Matlab..

Rys. 2.26. Transmitancja obiektu wybranego modelu ARMAX (amx4441)

Odpowiedzi skokową i impulsową wybranego modelu amx10101 przedstawiono na rysunkach 2.27 i 2.28.

Rys. 2.27. Odpowiedź skokowa (amx4441)	Rys. 2.28. Odpowiedź impulsowa (amx4441)

2.5. Estymacja wielomianowa (Polynomial Models) - Output-Error

W oparciu o dane wejściowe i wyjściowe (wektory wymuszenia oraz odpowiedzi układu) utworzone zostały modele wielomianowe (ang. Polynomial Models) OE (ang. Output-Error) wykorzystując System Identification Tool. W oknie głównym narzędzia, przedstawionym na rysunku 2.29 widoczne są miniaturki opracowywanych charakterystyk w zakresie od oe111 do oe12121. Charakterystyki w powiększeniu zostały pokazane na rysunkach 2.30, 2.31, 2.32 i 2.33. Cyfry zawierające się w nazwie charakterystyki informuje o stopniach wielomianów: B + 1 - nb, F - nf opóźnienia wejścia wyrażone w liczbie próbek - nk, w postaci OE:[nb nf nk].

Rys. 2.29. Okno System Identification Tool z danymi estymowanymi modelem OE

Rys. 2.30. Charakterystyki: oe111, oe221, oe331	Rys. 2.31. Charakterystyki: oe441, oe551, oe661
Rys. 2.32. Charakterystyki: oe771, oe881, oe991	Rys. 2.33. oe10101, oe11111, oe12121

Rysunek 2.34 przedstawia charakterystykę porównawczą odpowiedzi skokowych wszystkich otrzymanych modeli OE wraz z wskaźnikiem stopnia dopasowania do wzorca.

Rys. 2.34. Charakterystyka porównawcza wszystkich otrzymanych modeli OE wraz z współczynnikiem dopasowania

Biorąc pod uwagę kryterium sztywności modelu oraz optymalny poziom dopasowania wybrano model oe441 (nb = 4, nf = 4, nk = 1). Rysunek 2.26 przedstawia transmitancję dyskretną wybranego modelu oe441 otrzymaną poprzez zastosowanie funkcji idpoly() środowiska Matlab..

Rys. 2.35. Transmitancja obiektu wybranego modelu OE (oe441)

Odpowiedzi skokową i impulsową wybranego modelu oe441 przedstawiono na rysunkach 2.36 i 2.37.

Rys. 2.36. Odpowiedź skokowa (oe441)	Rys. 2.37. Odpowiedź impulsowa (oe441)

2.6. Estymacja wielomianowa (Polynomial Models) - Box-Jenkins

W oparciu o dane wejściowe i wyjściowe (wektory wymuszenia oraz odpowiedzi układu) utworzone zostały modele wielomianowe (ang. Polynomial Models) BJ (ang. Box-Jenkins) wykorzystując System Identification Tool. W oknie głównym narzędzia, przedstawionym na rysunku 2.38 widoczne są miniaturki opracowywanych charakterystyk w zakresie od bj11111 do bj12121212121. Charakterystyki w powiększeniu zostały pokazane na rysunkach 2.39, 2.40, 2.41 i 2.42. Cyfry zawierające się w nazwie charakterystyki informuje o stopniach wielomianów: B + 1 - nb, C + 1 - nf, D + 1 - nd, F + 1 - nf oraz opóźnienia wejścia wyrażone w liczbie próbek - nk, w postaci BJ:[nb nc nf nd nk].

Rys. 2.38. Okno System Identification Tool z danymi estymowanymi modelem BJ

Rys. 2.39. Charakterystyki: bj11111, bj22221, bj33331	Rys. 2.40. Charakterystyki: bj44441, bj55551, bj66661
Rys. 2.41. Charakterystyki: bj77771, bj88881, bj99991	Rys. 2.42. Charakterystyki: 101010101, 111111111, 121212121

Rysunek 2.43 przedstawia charakterystykę porównawczą odpowiedzi skokowych wszystkich otrzymanych modeli BJ wraz z wskaźnikiem stopnia dopasowania do wzorca.

Rys. 2.43. Charakterystyka porównawcza wszystkich otrzymanych modeli BJ wraz z współczynnikiem dopasowania

Biorąc pod uwagę kryterium sztywności modelu oraz optymalny poziom dopasowania wybrano model bj11111 (nb = 1, nc = 1, nfd = 1, nd = 1, nk = 1). Rysunek 2.26 przedstawia transmitancję dyskretną wybranego modelu bj11111 otrzymaną poprzez zastosowanie funkcji idpoly() środowiska Matlab..

Rys. 2.44. Transmitancja obiektu wybranego modelu BJ (bj11111)

Odpowiedzi skokową i impulsową wybranego modelu bj11111 przedstawiono na rysunkach 2.45 i 2.46.

Rys. 2.45. Odpowiedź skokowa (bj11111)	Rys. 2.46. Odpowiedź impulsowa (bj11111)

3. WNIOSKI

Wykonane w ramach projektu zadanie identyfikacji obiektu BlackBoxA zostało wykonane wykorzystując modele: Transfer Function Models i Polynomial Models. Każdy z przedstawionych modeli pozwolił uzyskać transmitancję o zadawalającym współczynniku podobieństwa. Narzędzie jakim jest System Identification Tool oferuje szeroki wachlarz narzędzi pozwalających na szybką i łatwą estymację parametrów wybranego modelu i otrzymanie obiektu wyjściowego w postaci transmitancji. Wygodny interfejs pozwala na łatwą interakcję z narzędziem. Analizując otrzymane modele można stwierdzić, że otrzymanie modelu matematycznego obiektu jest stosunkowo proste. Kluczowym aspektem jest odpowiedni poziom dopasowania parametrów modelu w stosunku do obiektu rzeczywistego. Zbyt niskie jak i zbyt wysokie wartości stopni estymowanych wielomianów powodują zbyt duży poziom niedopasowania lub przesztywnienie modelu, co jest niekorzystne. Oznaką przesztywnienia modelu jest bardzo dobre odwzorowanie dla danych, na których obiekt był identyfikowany oraz nieoczekiwane wartości sygnału wyjściowego w przypadku podania innych danych testowych na wejście obiektu.

4. ZAŁĄCZNIKI

Listing skryptu a1.m
clear all; close all; sim('identM.mdl',[1,150]); rozmiar = size(tout(2:size(tout,1)),1); y1 = outS(1:rozmiar,1); u1 = outS(1:rozmiar,2); y2 = outI(1:rozmiar,1); u2 = outI(1:rozmiar,2);

Listing skryptu a2.m
G1 = tf(tf1.num,tf1.den,1,'variable','z^-1'); G2 = tf(tf2.num,tf2.den,1,'variable','z^-1'); evalc('G1') evalc('G2') %G = G1; G = G2;

Listing skryptu a3.m

format compact; [yS, tS] = step(G,50); [yI, tI] = impulse(G,50); [XX1, YY1] = stairs(tS, yS); [XX2, YY2] = stairs(tout(2:size(tS,1)+1), y1(2:size(tS,1)+1)); [XX3, YY3] = stairs(tI, yI); [XX4, YY4] = stairs(tout(2:size(tI,1)+1), y2(2:size(tS,1)+1)); figure(1) hold on; grid on; plot(XX1(1:90), YY1(1:90),'-r'); plot(XX2(1:90), YY2(1:90),'-b'); title('Odpowiedzi skokowe'); legend('Model','Obiekt'); xlabel('Time t[s]'); ylabel('Response'); figure(2) hold on; grid on; plot(XX3, 10*YY3,'-r'); plot(XX4, YY4,'-b'); title('Odpowiedzi impulsowe'); legend('Model','Obiekt'); xlabel('Time t[s]'); ylabel('Response'); evalc('G')