Identyfikacja Procesów Technologicznych, Identyfikacja parametryczna obiektów dynamicznych, Nr �wiczenia :

Laboratorium

identyfikacji

proces�w

technologicznych

Wykonali :

Kierunek :

Numer wiczenia :

Temat :

Identyfikacja parametryczna

obiekt�w dynamicznych.

Data wykonania :

Data oddania :

Ocena :

Celem wiczenia byo przeprowadzenie procesu identyfikacji danego obiektu dynamicznego metod parametryczn, tworzc model typu ARX. W tym celu naleao skorzysta z metody MNK oraz zmiennej instrumental nej.

1. Wstp teoretyczny.

Zbyt skomplikowane lub wrcz nieznane prawa rzdzce modelowanym ukadem dynamicznym powoduj potrzeb przeprowadzenia analizy zwizk�w midzy jego sygnaami wejciowymi a wyjciowymi, czyli identyfikacj. Stosuje si dwa rodzaje identyfikacji : biern - jestemy biernymi obserwatorami tego, co wpywa na ukad i co z niego wychodzi; oraz czynn - czynnie oddziaywujemy na ukad pobudzajc jego wejcia. Modele wykorzystywane przez metody identyfikacji s modelami dyskretnymi. Spowodowane to jest dostpnymi sposobami pomiaru sygna�w wejciowych i wyjciowych badanego ukadu. Dopiero od modelu dyskretnego moemy przej do modelu w postaci cigej. Na wynik obserwacji, jak i na przebieg samego zjawiska wpyw maj nie tylko nasze wiadome oddziaywania, ale take niezalene od nas zak�cenia. Wynikiem tej sytuacji jest bd, kt�rym s obarczone wszystkie dane uzyskane do cel�w identyfikacji. Metody identyfikacji musz zatem uwzgldnia wpyw zjawisk przypadkowych.

Biblioteka Identification Toolbox Matlaba zawiera metody identyfikacji modeli dynamicznych o r�nych postaciach. R�nice pomidzy poszczeg�lnymi modelami wynikaj ze sposobu uwzgldnienia wpywu zak�ce na funkcjonowanie ukadu. Jednym z nich jest model

ARX ( AutoRegressive with eXogenous input )

przedstawiony r�wnaniem :

(1)

gdzie :

- na - liczba czynnik�w zwizanych z sygnaem wyjciowym y,

- nb - liczba czynnik�w zwizanych z sygnaem wejciowym u,

- nk - liczba takt�w op�nienia wystpujcego w ukadzie.

Nazwa ARX pochodzi z jzyka angielskiego i oznacza, e jest to model, w kt�rym aktualna warto sygnau wyjciowego zaley od :

- pewnej liczby poprzednich wartoci sygnau wyjciowego y, co okrela si mianem autoregresji (autoregressive),

- pewnej liczby kolejnych wartoci sygnau wejciowego u (exogenous input),

- od wartoci sygnau zak�cajcego e (zwizanego z szumem biaym).

Istnieje wiele metod prowadzcych do znalezienia parametr�w wyej przedstawionego modelu. Dokonano ich podziau na dwie grupy : parametryczne i nieparametryczne. Metody parametryczne wymagaj podania parametr�w modelu, czyli w naszym przypadku okrelenia wsp�czynnik�w [na,nb,nk]. Metody nieparametryczne same okrelaj parametry identyfikowanego modelu.

Najprostsz metod parametryczn stanowi metoda najmniejszych kwadrat�w (kt�ra zostaa ju przedstawiona w w.1 w odniesieniu do identyfikacji charakterystyki statycznej). Chcc znale model okrelajcy zalenoci midzy sygnaem wejciowym u a wyjciowym y na podstawie N pr�bek odpowiadajcym kolejnym chwilom czasu kT_P, k=1,2,...,N, poszukujemy wsp�czynnik�w modelu okrelonego r�wnaniem (1). R�wnanie to jest spenione jedynie w przyblieniu, wic dodajemy czynnik uwzgldniajcy (reprezentujcy) bd e_k i zapisujemy w postaci wektorowej :

0x01 graphic
,

oznaczajc :

otrzymujemy :

Uwzgldniajc fakt, e dysponujemy N wartociami sygnau y i u, w wyniku czego otrzymujemy ukad N r�wna o analogicznej do powyszej postaci :

0x01 graphic

przy czym naley uwzgldni procesy przejciowe (pocztkowe), co spowoduje, e ilo element�w w macierzach ulegnie zmianie. Po oznaczeniu wektora wartoci sygnau wyjciowego y jako Y, macierzy wartoci wyj i sterowa jako ", oraz wektora bd�w przyblienia jako E otrzymujemy ukad r�wna :

w kt�rym niewiadom jest wektor parametr�w .

Tak przedstawione zadanie jest identyczne z zadaniem rozwizania ukadu r�wna om�wionego w w. 1. Mona je wic rozwiza za pomoc operatora \ lub /. Znacznie wygodniejsz jednak funkcj jest w tym przypadku funkcja arx realizujca wanie t operacj. Otrzymujemy w�wczas oszacowane wartoci wsp�czynnik�w r�wnania modelu (1), kt�re stanowi estymator parametr�w modelu . Wartoci bdu e_k mona potraktowa jako losowe zak�cenie. Kiedy sygna wejciowy, wyjciowy i zak�cenie s ze sob skorelowane, estymator dla metody najmniejszych kwadrat�w jest obciony - jego warto oczekiwana r�ni si od wartoci rzeczywistej, a wic powtarzajc identyfikacj wielokrotnie i obliczajc redni uzyskanych wynik�w otrzymamy warto przesunit w stosunku do wartoci rzeczywistej parametr�w. Wpyw tej korelacji mona zmniejszy metod zmiennej instrumentalnej, kt�rej istota polega na takim przekonstruowaniu estymatora, by wspomniana korelacja miaa jak najmniejszy wpyw na jego statystyczne wasnoci.

Procedur identyfikacji mona podzieli na kilka etap�w :

- gromadzenie i przygotowanie danych pomiarowych,

- przyjcie postaci (struktury) modelu i okrelenie jego parametr�w za pomoc okrelonej metody identyfikacji,

- weryfikacja modelu.

2. Przebieg wiczenia.

2.1. Zebranie danych pomiarowych.

W pierwszym punkcie wiczenia dokonujemy symulacji badanego ukadu dynamicznego za pomoc blok�w dostpnych w Simulinku, tworzc wasny podsystem (Subsystem).

0x01 graphic

Struktura wewntrzna tego podsystemu przedstawiona jest poniej :

0x01 graphic

Po zestawieniu analizowanego elementu, umieszczamy go w ukadzie pomiarowym, sucym do zebrania danych pomiarowych. Dokonane zostan dwie serie pomiarowe. Kada skadaa si bdzie z N=300 pomiar�w. W wyniku pierwszej serii otrzymamy dane, kt�re bd wykorzystane w dalszym etapie procesu identyfikacji do budowy modelu metod parametryczn. Dane z tej serii zostan zawarte w macierzy Z1, o strukturze :

0x01 graphic

Ukad pomiarowy w kt�rym wykonana zostanie pierwsza seria pomiar�w :

0x01 graphic

natomiast graficzne przedstawienie wynik�w otrzymanych z pomiar�w wyglda w spos�b nastpujcy: 0x01 graphic

Przebiegi te uzyskujemy wywoujc funkcj :

plot(dana)

w naszym przypadku argument `dana' stanowi macierz zawierajca wyniki pomiar�w. Linia o szarym odcieniu przedstawia przebieg sygnau wejciowego u (wymuszenia), natomiast linie o kolorze czarnym stanowi odpowied ukadu.

Druga seria pomiar�w zostanie wykonana w celu uzyskania danych w postaci macierzy Z2 o takiej samej strukturze jak poprzednia macierz Z1. Dane te natomiast wykorzystamy do testowania (weryfikacji) otrzymanego modelu, dlatego te przed wykonaniem pomiar�w zmieniamy ustawienia (Seed) bloku Band-Limited White Noise (generator szumu biaego), aby otrzyma inny zestaw danych ni w serii pierwszej. Wyniki pomiar�w przekazane zostan do przestrzeni roboczej Matlaba po zmianie argumentu bloku To Workspace. Tak, wic ukad pomiarowy w kt�rym wykonana zostanie druga seria pomiar�w :

0x01 graphic

Posta graficzna pomiar�w drugiej serii :

0x01 graphic

2.2. Przygotowanie danych pomiarowych.

Zadanie przygotowania danych otrzymanych w punkcie wyej (z pomiar�w) do ich dalszej obr�bki w procesie identyfikacji ma na celu usunicie z nich wartoci redniej. Wykonanie tego zadania umoliwia nam funkcja:

dana1 = dtrend(dana)

Argumentem tej funkcji (`dana') s w naszym przypadku macierze otrzymane w punkcie 2.1. Wynikiem zastosowania tej funkcji s macierze Z11 oraz Z22 o budowie identycznej jak Z1 i Z2. Std dane do tworzenia modelu (teraz ju Z11) przyjmuj posta:

0x01 graphic

a dane z drugiej serii, czyli Z22 :

0x01 graphic

2.3. Identyfikacja (okrelanie parametr�w) modelu ARX za pomoc metody MNK.

2.3.1. Dob�r struktury.

Majc przygotowane dane pomiarowe przechodzimy do nastpnego etapu identyfikacji, w kt�rym okrelamy rodzaj modelu oraz jego struktur. Posugujemy si metod parametryczn, a konkretnie metod najmniejszych kwadrat�w dla kt�rej wybieramy model typu ARX. Funkcj realizujc identyfikacj modelu o takiej postaci metod MNK (czyli przybliajc w sensie najmniejszych kwadrat�w parametry modelu ARX) jest funkcja arx , naleca do biblioteki Identification Toolbox pakietu Matlab. Posta wywoania tej funkcji :

th = arx(z,nn)

gdzie :

- z=[y u], gdzie y i u s kolumnowymi wektorami zawierajcymi zmierzone wartoci sygna�w wejciowych (u) i wyjciowych (y),

- nn=[na nb nk], gdzie na , nb i nk s wsp�czynnikami r�wnania (1) opisujcego model.

Poniewa nie znamy parametru `nn' modelu (licznika i mianownika transmitancji oraz op�nienia), wic zakadamy je wstpnie jak najnisze, czyli nn=[1 1 1]. Dla tak zaoonej struktury modelu oraz dla zestawu danych pomiarowych Z11 stosujemy funkcj arx. Funkcja ta zwraca wyniki identyfikacji w postaci tzw. macierzy rezultat�w THETA o specjalnej budowie, z kt�rej potrzebne informacje mona uzyska `rcznie' lub za pomoc odpowiednich funkcji. Dlatego uywamy funkcji :

[num, den] = th2tf(th)

kt�ra zwraca dwa wektory : num - reprezentujcy licznik i den - reprezentujcy mianownik transmitancji dyskretnej otrzymanego modelu, na podstawie informacji zawartych wanie w macierzy THETA. Stosujc nastpnie funkcj printsys uzyskujemy t transmitancj w postaci :

W ten spos�b - dla wybranej struktury - otrzymalimy model ARX analizowanego ukadu dynamicznego o postaci :

2.3.2. Weryfikacja modelu.

Aby zweryfikowa otrzymane w poprzednim punkcie wyniki identyfikacji, wykorzystamy metod symulacyjn. Polega ona na por�wnaniu graficznych postaci odpowiedzi modelu oraz obiektu, W tym celu wykorzystujemy zestaw danych uzyskanych w drugiej serii pomiarowej, czyli macierz Z22. Funkcja Matlaba wykonujca tak postawione zadanie ma posta :

yh = compare(Z22,th)

Wida, e r�wnie ta funkcja korzysta z danych zawartych w macierzy THETA. Wynikiem zastosowania funkcji compare jest :

0x01 graphic

Wsp�czynnik Fit okrela w jakim stopniu odpowied modelu jest zbiena z odpowiedzi badanego obiektu. W przypadku idealnego pokrywania si tych dw�ch przebieg�w przyjmuje on warto r�wn 0.

2.3.3. Wyb�r `najlepszej' struktury modelu.

Por�wnujc przebiegi z punktu 2.4. mona zauway, e odpowied modelu o strukturze (najprostszej) zaoonej w punkcie 2.3. znacznie odbiega od odpowiedzi rzeczywistego obiektu. Dlatego staramy si tak zmodyfikowa struktur modelu aby jak najbardziej przybliy jego wasnoci do wasnoci obiektu, kierujc si przy tym wartociami wsp�czynnika Fit. Std, powtarzamy punkty 2.4. i 2.5. dla r�nych struktur modeli. W wyniku takiego postpowania otrzymujemy nastpujce pary :

nn=[na nb nk]	Fit
[2 2 1]	0,14460
[1 2 2]	0,15980
[1 2 1]	0,19000
[3 2 1]	0,09428
[2 3 3]	0,12600
[3 3 3]	0,11980
[3 2 3]	0,11900
[3 2 2]	0,09370
[3 3 2]	0,09048
[2 2 3]	0,13300

Najmniejsz wartoci wsp�czynnika Fit - wr�d struktur zawartych w tabeli - charakteryzuje si model ARX o strukturze okrelonej wektorem nn=[3 3 2], kt�rego transmitancja dyskretna wynosi :

a posta :

0x01 graphic

natomiast wynik dziaania funkcji compare, w odniesieniu do modelu o takiej postaci przedstawia si nastpujco :

0x01 graphic

W pakiecie Matlaba istnieje jednak funkcja :

v = arxstruc(Z,Z,n)

samoczynnie wyszukujca najkorzystniejsz struktur modelu z wybranych struktur umieszczonych w tzw. macierzy kandydat�w na struktury n, wykorzystujc dane pomiarowe Z.W naszym wiczeniu utworzylimy nastpujc macierz n:

0x01 graphic

na nb nk

Wyniki tej funkcji zapisywane s z kolei w macierzy v, kt�rej pierwszy element kadego wiersza przedstawia odchyk (bd) charakteryzujc zdolno odtwarzania wasnoci obiektu przez model o strukturze okrelonej przez wektor nn stanowicy pozostae elementy kadej kolumny. Analizujc t macierz wnioskujemy, e struktur dla kt�rej pierwszy element kolumny n przybiera najmniejsz warto, jest podobnie jak w punkcie 2.5. nn=[3 3 2]. Traktujemy macierz v jako argument funkcji :

nn1 = selstruc(v,0)

kt�ra ma na celu wyszukanie wr�d struktur zawartych w niej, struktury o najmniejszej odchyce. Otrzymamy ponownie wektor nn1=[3 3 2].

Aby przeanalizowa w procesie identyfikacji wiksz liczb struktur modeli ni dotychczas byy rozpatrywane, poprzez `rczne' zamieszczenie ich przez nas w macierzy n, naley zastosowa funkcj :

nn2 = struc(1:4, 1:4, 1:4)

Wygeneruje ona macierz nn2 (o budowie podobnej jak macierz n), zwierajc wszystkie moliwe kombinacje struktur modeli, czyli wektor�w nn[na nb nk]. Poszczeg�lne elementy kadego wektora bd z zakresu 1..4. Wywoujc nastpnie funkcj :

v1 = arxstruc(z11,z11,nn2);

otrzymamy macierz v1 (o strukturze podobnej do omawianej wyej macierzy v). Funkcja :

nn3 = selstruc(v1,0)

zwraca teraz wektor nn3=[4 4 1]. Wyznaczamy parametry modelu ARX dla tej struktury :

th = arx(z11,nn3)

[num,den] = th2tf(th)

printsys(num,den,'z')

std transmitancja modelu :

R�wnanie (1) uzyskanego modelu :

0x01 graphic

Dokonujc weryfikacji symulacyjnej funkcj :

compare(z22,th)

otrzymamy przebiegi :

0x01 graphic

2.4. Identyfikacja modelu ARX za pomoc metody zmiennej instrumentalnej.

W tym punkcie wiczenia wykorzystamy procedur stosowan w przypadku MNK, w niekt�rych tylko funkcjach naley okreli wyb�r metody zmiennej instrumentalnej,tak wic :

2.4.1. Dob�r `najlepszej' struktury modelu.

nn4 = struc(1:4, 1:4, 1:4)

v1 = ivstruc(z11,z11,nn4)

nn5 = selstruc(v1,0)

Taka sekwencja funkcji prowadzi do uzyskania struktury nn5=[4 1 1]. Wyznaczamy parametry modelu ARX dla tej struktury :

th = iv4(z11,nn5)

Funkcja ta przyblia parametry modelu ARX , korzystajc z techniki zmiennej instrumentalnej. Jej dziaanie odpowiada funkcji arx dla metody MNK. Wynikiem wykonania powyszych funkcji jest odpowiednia macierz THETA, skd wyr�niamy interesujce nas informacje :

[num,den] = th2tf(th)

printsys(num,den,'z')

czyli :

std posta modelu ARX :

2.4.2. Weryfikacja modelu.

Ten etap identyfikacji wykonujemy podobnie jak wczeniej :

compare(z22,th)

0x01 graphic

2.5. Por�wnanie (graficznych odpowiedzi) model�w ARX otrzymanych metod MNK oraz metod

zmiennej instrumentalnej.

Por�wnania dokonujemy w ukadzie :

0x01 graphic

Przebiegi wejcia i wyjcia obiektu oraz odpowied modelu ARX otrzymanego metod MNK dla sygnau wejciowego o czstotliwoci f=1 i amplitudzie A=1.

0x01 graphic

Przebiegi wejcia i wyjcia obiektu oraz odpowied modelu ARX otrzymanego metod zmiennej instrumentalnej dla sygnau wejciowego o czstotliwoci f=1 i amplitudzie A=1.

0x01 graphic

Przebiegi wejcia i wyjcia obiektu oraz odpowied modelu ARX otrzymanego metod MNK dla sygnau wejciowego o czstotliwoci f=1 i amplitudzie A=2.

0x01 graphic

Przebiegi wejcia i wyjcia obiektu oraz odpowied modelu ARX otrzymanego metod zmiennej instrumentalnej dla sygnau wejciowego o czstotliwoci f=1 i amplitudzie A=2.

0x01 graphic

Przebiegi wejcia i wyjcia obiektu oraz odpowied modelu ARX otrzymanego metod MNK dla sygnau wejciowego o czstotliwoci f=5 i amplitudzie A=1.

0x01 graphic

Przebiegi wejcia i wyjcia obiektu oraz odpowied modelu ARX otrzymanego metod zmiennej instrumentalnej dla sygnau wejciowego o czstotliwoci f=5 i amplitudzie A=1.

Posta m-pliku.

z11=dtrend(z1);

z22=dtrend(z2);

plot(z22);

hold on;

plot(z11);

nn=[1 1 1]

nn=[2 2 1]

nn=[1 2 2]

nn=[1 2 1]

nn=[3 2 1]

nn=[2 3 3]

nn=[3 3 3]

nn=[3 2 3]

nn=[3 2 2]

nn=[3 3 2]

nn=[2 2 3]

th=arx(z11,nn);

yh=compare(z22,th);

[num,den]=th2tf(th)

printsys(num,den,'z')

n=[2 2 1;1 2 3;2 2 3; 3 3 2; 3 2 2; 3 2 3; 3 2 1];

v=arxstruc(z11,z11,n)

nn1=selstruc(v,0);

nn2=struc(1:4,1:4,1:4);

v1=arxstruc(z11,z11,nn2);

nn3=selstruc(v1,0)

th=arx(z11,nn3);

%compare(z22,th);

[num,den]=th2tf(th);

printsys(num,den,'z')

nn4=struc(1:4,1:4,1:4);

v1=ivstruc(z11,z11,nn4);

nn5=selstruc(v1,0)

th=iv4(z11,nn5);

%compare(z22,th);

[num1,den1]=th2tf(th);

printsys(num1,den1,'z')

3. Uwagi i wnioski.

Powysze wiczenie zostao wykonane zgodnie z zaoonym wczeniej schematem postpowania. I tak pierwszym krokiem byo zgromadzenie danych pomiarowych. Posuyy one p�niej do tworzenia modelu ARX oraz jego weryfikacji. Nastpnym etapem przeprowadzanego procesu identyfikacji parametrycznej obiektu dynamicznego by wyb�r struktury tworzonego modelu. Wykonywalimy to w spos�b `rczny' - metod pr�b i bd�w - ale wykorzystywalimy r�wnie w tym celu funkcje dostpne w bibliotekach Matlaba, kt�re najpierw generoway moliwe struktury, a nastpnie dokonyway wyboru najlepszej z nich. W ten spos�b mona przeanalizowa znacznie wikszej liczby struktur. Majc ju zaoon struktur przystpowalimy do wyznaczania parametr�w modelu. Wynik tej operacji otrzymywalimy w postaci transmitancji dyskretnej. Ostatnim punktem identyfikacji bya weryfikacja otrzymanego modelu. Wykorzystywalimy tutaj weryfikacj graficzn (symulacyjn). Jej wynik stanowi podstaw do przyjcia lub odrzucenia wyznaczonego modelu i ponownego przeprowadzenia procesu identyfikacji ale ju dla zmodyfikowanych zaoe. W efekcie kocowym wydzielilimy dwa modele badanego elementu. Pierwszy z nich by modelem wyznaczanym przy wykorzystaniu metody najmniejszych kwadrat�w. W tym przypadku najkorzystniejsz okazaa si struktura [4 4 1]. Struktura [4 1 1] okazaa si natomiast najodpowiedniejsz dla modelu identyfikowanego metod zmiennej instrumentalnej, kt�ry zwraca wyniki czsto znacznie dokadniejsze ni analogiczne otrzymane funkcj arx - a wic poprzez metod MNK. Staje si to szczeg�lnie odczuwalne gdy zak�cenia s silnie skorelowane z sygnaem wejciowym. Zostay take wyznaczone przykadowe przebiegi wyjciowe modelu oraz obiektu, dla zadanych wartoci sinusoidalnego sygnau wejciowego.