Identyfikacja systemów dynamicznych na przykładzie prostego obwodu elektrycznego
dla początkujących

Marcin Maślanka
student V roku, spec. Automatyka i Metrologia, IMiR, AGH

Streszczenie - Na podstawie danej odpowiedzi skokowej identyfikowanego układu RLC oraz znanej pojemności C, wyznaczono pozostałe parametry obwodu korzystając z trzech różnych metod identyfikacji. Przedstawiono podstawy teoretyczne wykorzystanych metod, dołączono także kody źródłowe m-plików realizujących wszystkie obliczenia.
Wykorzystano pakiet Matlab 5.3, Control System Toolbox 4.2, System Identification Toolbox 4.0.5

 SPIS TREŚCI

1. Przygotowanie dostarczonych danych empirycznych
2. Opis modelu rozpatrywanego systemu dynamicznego
3. Identyfikacja modelu parametrycznego ARX
    3.1 Wstęp teoretyczny
    3.2 Dobór optymalnej struktury modelu
    3.3 Wyznaczenie szukanych parametrów układu
4. Identyfikacja na podstawie analizy odpowiedzi skokowej
    4.1 Określenie cech charakterystycznych odpowiedzi skokowej
    4.2 Wyznaczenie szukanych parametrów układu
    4.3 Test symulacyjny
5. Identyfikacja na podstawie analizy charakterystyki amplitudowo – fazowej
    5.1 Wykreślenie charakterystyki amplitudowo-fazowej i lokalizacja punktów charakterystycznych
    5.2 Wyznaczenie szukanych parametrów układu
    5.3 Test symulacyjny
6. Podsumowanie
    6.1 Zestawienie wyników
    6.2 Komentarz
Literatura

1. Przygotowanie dostarczonych danych empirycznych

    Dostarczone dane empiryczne (zbiór obserwacji) obejmują szeregi czasowe sygnału wymuszającego oraz odpowiedzi układu na sygnał wymuszający przy znanej częstotliwości próbkowania. Plik part0.m (2kB) odczytuje w środowisku Matlab dostarczone w formie MAT-pliku dane (rr17.mat - 48kB) i przedstawia je w formie graficznej.
Rys.1 (8kB) przedstawia wynik działania funkcji part0.m (2kB).
Dane empiryczne obejmują odpowiedź układu na wymuszenie sygnałem prostokątnym o jednostkowej amplitudzie, częstotliwości 2.5 [Hz] i czasie trwania 2.5 okresu. Częstotliwość próbkowania wynosi 20 [kHz], ilość próbek: 20000.
Do celów identyfikacji zostaje wyodrębniony jeden z przebiegów przejściowego stanu nieustalonego odpowiedzi. Wybrano przebieg stanu nieustalonego odpowiadający trzeciemu półokresowi sygnału wymuszenia. Przebieg ten odpowiada fragmentom dostarczonych szeregów czasowych o indeksach od 801 do 1200.
Zakłada się liniowość badanego układu. Jeżeli zatem parametry układu nie zależą od amplitudy sygnału wymuszającego można przeskalować badany przebieg stanu nieustalonego do przebiegu jak na rys.2 (8kB), dokonując operacji zapisanej w notacji Matlaba jako: y=(y-min(y))/2; gdzie: y- odpowiedź układu na wymuszenie.

Tak przygotowany zestaw danych wykorzystywany będzie w stosowanych procedurach identyfikacji rozpatrywanego systemu.
Przyjmuje się następujące oznaczenia:
u – szereg czasowy wymuszający (wg rys.2 - 8kB)
y – szereg czasowy odpowiedzi układu na sygnał wymuszający (wg rys.2 - 8kB)

Parametry sygnałów u, y: Fs= 20 [kHz], t₀= 0, ilość próbek: 400.
Znana pojemność układu wynosi: C= 0.6[mikro F].
 

2. Opis modelu rozpatrywanego systemu dynamicznego

    Rozpatrywane jest połączenie szeregowe elementów RLC jak na rysunku poniżej.

C= 0.6[mikro F]

    Za model tego układu przyjmuje się liniowy model SISO (wartość amplitudy sygnału wejściowego nie ma wpływu na parametry modelu, rozpatrywane jest jedno wejście i jedno wyjście), sygnałem wejściowym do układu jest napięcie U₁, sygnałem wyjściowym jest napięcie natomiast U₂.
Poniżej wyznaczona zostaje transmitancja operatorowa tego modelu.

Napięcie na rezystorze wynosi:

mierzone w tym układzie napięcie na kondensatorze wynosi:

napięcie na cewce:

zatem zapisując równanie Kirchhoffa dla tego obwodu otrzymuje się:

Po przekształceniu Laplace’a:

zatem transmitancja operatorowa układu:

Ogólna postać transmitancji członu oscylacyjnego ma postać:

czyli:


Przy znanej pojemności C wyznaczenie parametrów R i L sprowadza się zatem do wyznaczenia podanych powyżej współczynników mianownika transmitancji operatorowej rozpatrywanego układu.
 

    W równaniach różniczkowych, za pomocą których opisywane są zjawiska dynamiczne modelowanych układów rzeczywistych, lub w ich dyskretnych odpowiednikach – równaniach różnicowych występuje skończona liczba parametrów, modele takie nazywa się stąd modelami parametrycznymi.
W przyborniku Matlaba System Identification Toolbox różnice pomiędzy poszczególnymi modelami parametrycznymi wynikają ze sposobu uwzględnienia wpływu zakłóceń na funkcjonowanie układu. Jednym z nich jest model autoregresji z zewnętrznym wymuszeniem ARX (AutoRegressive with eXogenous input) przedstawiony równaniem [1]:


gdzie:
u(t) - wejście obserwowane w dyskretnej chwili t,
y(t) - wyjście obserwowane w dyskretnej chwili t,
v(t) - zastępcze zakłócenie obserwowane w dyskretnej chwili t (nie skorelowany sygnał losowy typu „biały szum” o rozkładzie normalnym, z zerową wartością średnią i stałą wariancją),
t – zmienna losowa o zerowej wartości średniej i stałej wariancji
na - liczba czynników związanych z sygnałem wyjściowym y,
nb - liczba czynników związanych z sygnałem wejściowym u,

Aktualna wartość sygnału wyjściowego zależy zatem od:
    -  pewnej liczby na poprzednich wartości sygnału wyjściowego y (autoregresja),
    -  pewnej liczby kolejnych wartości sygnału wejściowego u,
    -  d wartości sygnału zakłócającego (oznaczonego przez epsilon).

Pierwiastki mianownika modelu nazywają się biegunami modelu, a pierwiastki licznika zerami. Model ARX jest stabilny gdy wszystkie jego bieguny leżą wewnątrz koła jednostkowego na płaszczyźnie zmiennej zespolonej z i jest minimalnofazowy jeśli warunek ten spełniają również jego zera. Struktura modelu jest jednoznacznie określona za pomocą trzech parametrów na, nb, nk.
Najprostszą metodę parametryczną stanowi metoda najmniejszych kwadratów. Chcąc znaleźć model określający zależności między sygnałem wejściowym u a wyjściowym y na podstawie N próbek odpowiadającym kolejnym chwilom czasu kT_P, k=1,2,...,N, poszukujemy współczynników modelu określonego zapisanym powyżej równaniem.

Równanie to jest spełnione jedynie w przybliżeniu, kiedy sygnał wejściowy, wyjściowy i zakłócenie są ze sobą skorelowane, estymator Q dla metody najmniejszych kwadratów jest obciążony - jego wartość oczekiwana różni się od wartości rzeczywistej, a więc powtarzając identyfikację wielokrotnie i obliczając średnią uzyskanych wyników otrzymamy wartość przesuniętą w stosunku do wartości rzeczywistej parametrów. W praktyce określa się estymator jako nieobciążony gdy stosunek sygnału do szumu jest mniejszy niż 10%.

Funkcją realizującą identyfikację modelu o takiej postaci metodą MNK (czyli przybliżającą w sensie najmniejszych kwadratów parametry modelu ARX) jest funkcja arx , należąca do biblioteki Identification Toolbox pakietu Matlab.
Postać wywołania tej funkcji :

th= arx([y u],[na nb nk])

gdzie:
y, u - wektory kolumnowe we/wy,
[na nb nk] - współczynniki zapisanego powyżej równania opisującego model, 'na' określa liczbę biegunów (pierwiastków mianownika), 'nb-1' określa liczbę zer (pierwiastków licznika),  
th - macierz wyniku w formacie THETA, z której poleceniem th2tf wydobywa się informację o współczynnikach licznika i mianownika transmitancji dyskretnej otrzymanego modelu. 

Ponieważ wektor współczynników [na nb nk] określający rząd licznika i mianownika transmitancji oraz opóźnienie nie jest znany poszukuje się go według kryterium minimalizacji błędu określanego podczas weryfikacji modelu.

    Optymalnego doboru struktury modelu dokonuje się metodami iteracyjnymi według kryterium minimalizacji błędu średniokwadratowego (mean square fit).
Weryfikacji modelu otrzymanego dla przyjętej arbitralnie struktury dokonuje się m.in. przy użyciu funkcji compare:
Uruchomienie następujących poleceń:

arx210 = arx([y u],[2 1 0]);
yh = compare([y u],arx210);

powoduje wyświetlenie wyników weryfikacji oraz wartości błędu średniokwadratowego (otrzymany wektor yh to szereg czasowy odpowiedzi symulowanej estymowanego modelu). Rys.4 (7kB) przedstawia przykładowy wynik weryfikacji modelu arx210 wygenerowany przez funkcję compare.

W celu szybkiego dokonania analizy większej liczby struktur modeli zastosowano postępowanie automatycznie generujące najlepsze wyniki weryfikacji poszukiwane wśród wszystkich możliwych struktur o parametrach z zakresu <1,4>.
Macierz parametrów wygenerowano poleceniem:

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

nn3 = selstruc(v1,c)

gdzie:
c - (np.: 'aic', 'mdl') kryterium informacyjne nadparametryzacji, według którego określana jest optymalna struktura modelu (w funkcji błędu średniokwadratowego, liczby estymowanych współczynników poszukiwanej transmitancji i ilości próbek danych, na podstawie których dokonywana jest estymacja).

Optymalna struktura zwracana przez funkcję selstruc to: [3 2 1] o wartości błędu średniokwadratowego: Fit= 0.015504.
Czasowy rozkład błędu przedstawiony został na rys.5 (8kB).

A oto jak przedstawiają się wartości błędu średniokwadratowego w zależności od dobranej struktury modelu:

Do dalszej analizy wybrano wyznaczoną przez funkcją selstruc strukturę o parametrach [3 2 1], (pomimo tego, że struktura [2 1 0] powinna teoretycznie najlepiej odwzorowywać opisany w pkt.2 model idealny):

arx321 = arx([y u],[3 2 1]);

3.3 Wyznaczenie szukanych parametrów układu

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

gd2= tf(l,m,1/Fs);

gc= d2c(gd2,’tustin’);

stąd: 

zatem dla tej struktury modelu:

L= 0.402 [H]        R= 314.40 [Om]

Wszystkie obliczenia realizuje funkcja part2.m (3kB).

4.3 Test symulacyjny

    Test symulacyjny polega na porównaniu odpowiedzi układu na zadane wymuszenie z symulowaną odpowiedzią na to wymuszenie estymowanego modelu. Test stanowi najprostszą z metod weryfikacji estymowanego modelu.

W celu dokonania testu symulacyjnego przeprowadzono symulację otrzymanego modelu przy zadanym wymuszeniu u, której wyniki porównano z szeregiem czasowym y.

Transmitancja operatorowa rozpatrywanego modelu:

num= [1]; den= [T0^2 2*T0*dzeta 1];
gc= tf(num,den);

gdzie:
T0 – okres drgań nietłumionych,
dzeta – współczynnik tłumienia.

Odpowiedź estymowanego modelu na wymuszenie 'u':

yh= lsim(gc,u,t);

Obliczenie błędu średniokwadratowego:

Fit = norm(yh -y)/sqrt(length(y));

Fit= 0.0155
Graficzne przedstawienie wyników testu symulacyjnego zamieszczono na rys.7 (8kB) oraz rys.8 (9kB).

5. Identyfikacja na podstawie analizy charakterystyki amplitudowo – fazowej
5.1 Wykreślenie charakterystyki amplitudowo-fazowej i lokalizacja punktów charakterystycznych

    Charakterystykę amplitudowo – fazową realizuje się przy znanym przebiegu odpowiedzi skokowej jako wykres na płaszczyźnie zespolonej transformaty Fouriera z szeregu czasowego odpowiedzi impulsowej układu, którą wyznacza się według zależności:

imp= [y(2)-y(1)  y(3)-y(2) ... y(n)-y(n-1)]

gdzie:
y – szereg czasowy odpowiedzi na skok jednostkowy,
n – ilość próbek szeregu y,
imp – szereg czasowy odpowiedzi impulsowej o ilości próbek n-1.

W celu szybkiego obliczenia zdefiniowanego powyżej szeregu odpowiedzi impulsowej wykorzystano służącą do tego funkcję diff. Polecenie:

imp = diff(y);

F= fft(imp,65536);
F= F(1:floor(length(F)/2));

której wykres na płaszczyźnie fazowej stanowi charakterystykę amplitudowo-fazową (rys.10 - 6kB).

Na przedstawionym na rys.10 (6kB) wykresie zlokalizowano punkty charakterystyczne, na podstawie których zostaną obliczone poszukiwane parametry.
Pierwszy z punktów (P45) określony dla Re= |Im| wyznaczony został jako punkt przecięcia ch-ki z funkcją Im= - Re (co zostało zobrazowane na wykresie), drugi punkt (P90) określony jest dla Re = 0.

Znalezione punkty:

P45= 901
P90= 1076

odpowiadające im punkty charakterystyki to F(P45) i F(P90).

5.2 Wyznaczenie szukanych parametrów układu

num= [1];den= [T0^2 2*T0*dzeta 1];
gc= tf(num,den);

yh = lsim(gc,u,t);

Fit = norm(yh - y)/sqrt(length(y));

¹ model o strukturze teoretycznie najlepiej dopasowanej do identyfikowanego układu
² model o najlepszym dopasowaniu do zbioru obserwacji [y u] wg kryterium informacyjnego Akaike'a
³ błąd średniokwadratowy porównania odpowiedzi symulowanej estymowanego modelu i odpowiedzi danej w zbiorze obserwacji

    Kluczowym problemem identyfikacji jest weryfikacja modelu. W opracowaniu tym estymowane modele poddawano jedynie testowi symulacyjnemu obliczając błąd średniokwadratowy jako kryterium określające zgodność modelu z dostarczonym zbiorem obserwacji. Okazuje się jednak, że kryterium to nie jest wystarczające, a wnioskować o tym można porównując wartości błędu dla różnych struktur oraz odpowiadające im wartości szukanych parametrów. Zestawienie wartości błędu dla różnych struktur (tabela), sugeruje że model [2 1 1] jest bardziej odpowiedni dla danego zbioru obserwacji niż model [2 1 0], modele o tych strukturach posiadają jednak takie same współczynniki mianownika transmitancji operatorowej (tabela), a zatem dają takie same wartości parametrów.

Przy wyborze kryterium określającego ze zbioru estymowanych modeli pewien szczególny model uważany za najbardziej odpowiedni brano pod uwagę kryterium informacyjne Akaike’a (w przypadku modeli ARX), wartość błędu średniokwadratowego i wreszcie informację o układzie rzeczywistym, którego funkcja przejścia była określona w pkt.2

Za najbardziej odpowiedni przyjmuje się model otrzymany na podstawie analizy charakterystyki amplitudowo-fazowej oraz model ARX z dwoma biegunami (na=2), bez zer (poszukiwano w liczniku jednego parametru stąd: nb=1) i bez opóźnienia nk=0 (modele z opóźnieniem nk=1 lub nk=2 określane są mniejszym błędem średniokwadratowym lecz dają takie same wartości szukanych parametrów tzn. przy wzroście nk modyfikowany jest jedynie licznik transmitancji operatorowej, który powinien być w przyjętym modelu jak najniższego rzędu).

Wyznaczane wartości parametrów wynosiły: R= 300 [Om], L= 0.4 [H] (podana pojemność układu: C= 0.6 [mikro F])