Akademia Górniczo-Hutnicza

im. Stanisława Staszica w Krakowie

Wydział Inżynierii Mechanicznej i Robotyki

Katedra Automatyzacji Procesów

Analiza sygnałów i identyfikacja

Projekt

Opracował

Jakub Kleszcz

Opis zadania

Tematem zadania jest identyfikacja obiektu dwuwejściowego przedstawionego na poniższym schemacie. Transmitancje Gu(s) i Gz(s) opisują obiekty RLC lub RC.

Otrzymałem zestaw danych numer 36. Zawierają się w nim następujące dane:

u - wejście w postaci szumu o rozkładzie równomiernym

y - wyjście obiektu

dt - czas próbkowania

ys - wyjście obiektu odizolowanego od zakłóceń

Dodatkowe dane do obiektów RLC i RC:

Ru = 90 [ohm]

Rz = 1 [kohm]

Układ RC jest elementem inercyjnym o transmitancji

G(s)=

Układ RLC jest elementem oscylacyjnym o transmitancji

G(s)=

1) Identyfikacja parametrów na podstawie odpowiedzi skokowej

Na podstawie odpowiedzi skokowej możemy rozróżnić, że nasz obiekt jest oscylacyjny. Będziemy posługiwać się wzorem

Następnie z wykresu odczytujemy wartości:

y1 = 0,7274 (pierwsze największe wychylenie)

y2 = 0,5287 (pierwsze najmniejsze wychylenie)

Tosc = 0,0058 (okres oscylacji)

Szukane parametry wyliczamy z następujących wzorów:

Obliczeń dokonujemy za pomocą programu Matlab. Wyniki przedstawiają się nastepująco:

$$\omega = 1083,3\frac{\text{rad}}{s}$$

ζ = 0.101

δ = 108, 891

$$\omega_{o} = 1088,8\frac{\text{rad}}{s}$$

T_o = 0.00091847s

Transmitancja Gu(s) przyjmuje postać:

$$\mathbf{G}_{\mathbf{u}}\mathbf{(s) =}\frac{\mathbf{1}}{\mathbf{8.436}\mathbf{*}\mathbf{10}^{\mathbf{-}\mathbf{7}}\mathbf{s}^{\mathbf{2}}\mathbf{+ 0.0001856s + 1}}$$

Możemy teraz wyliczyć C1 oraz L

C1 = 2.062*10⁻⁶F

L = 0.4091H

Aby upewnić się o poprawności otrzymanych przez nas wartości porównajmy odpowiedź skokową zadaną z tą przez nas zinterpretowaną

2) Identyfikacja parametrów modelu na podstawie odpowiedzi impulsowej

Teraz należy zidentyfikować model na podstawie odpowiedzi impulsowej, która będzie wyznaczona z odpowiedzi skokowej.

Odpowiedź impulsową uzyskujemy poprzez uzupełnienie odpowiedzi skokowej zerami i zróżniczkowanie.

Następnym krokiem jest wyznaczenie krzywej Nyquista, którą otrzymujemy za pomocą transformaty Fourierra. Na podstawie punktów przecięcia krzywej z osiami układu współrzędnych wyznaczamy parametry transmitancji.

Kolejny krok to znalezienie dwóch punktów na wykresie. Aby łatwiej było je znaleźć nakładamy na jednym wykresie część rzeczywistą i odbicie względem osi OX części urojonej. Pierwszy z punktów to miejsce przecięcia się wykresów. Drugi z kolei to miejsce, w którym Re = 0. Dla znalezionych punktów potrzebny nam jest jedynie numer próbki, czyli współrzędna X.

Znaczne przybliżenie wykresu pozwala nam bardzo dokładnie odczytać numer próbki. Odczytane numery to:

P1 = 1626 (miejsce przecięcia się wykresów)

P2 = 1790 (Re = 0)

Dzięki nim obliczyć możemy parametry występujące w szukanej transmitancji.

ω_{1 = 2 * π * P_Od1}

ω_{2 = 2 * π * P_Od2}

Parametry obliczone za pomocą programu Matlab:

$$\mathbf{G}_{\mathbf{u}}\mathbf{(s) =}\frac{\mathbf{1}}{\mathbf{7.906}\mathbf{*}\mathbf{10}^{\mathbf{-}\mathbf{7}}\mathbf{s}^{\mathbf{2}}\mathbf{+ 0.0001711s + 1}}$$

C1 = 1.9016*10⁻⁶ F

L = 0.4157 H

Dla pewności, czy znalezione parametry są poprawne dokonajmy porównania odpowiedzi skokowych układu zadanego i zinterpretowanego.

3) Identyfikacja parametrów modelu obiektu SISO

Identyfikacja tą metodą polega na wykorzystaniu transformat Fouriera sygnałów wyjściowego i wejściowego, wyznaczeniu gęstości widmowej mocy sygnału wejściowego oraz gęstości widmowej mocy wzajemnej wejścia i wyjścia. Dzięki temu wyznaczymy charakterystyki amplitudowo-fazowej. Podobnie jak w poprzednim punkcie znaleźć musimy charakterystyczne numery próbek.

Otrzymujemy następującą charakterystykę Nyquista:

W tym momencie powtarzamy proces z poprzedniego zadania: znajdujemy miejsce przecięcia się wykresów i próbki dla której Re = 0.

Odczytane punkty:

P1 = 162 (miejsce przecięcia się wykresów)

P2 = 178 (Re = 0)

Jesteśmy gotowi do obliczenia transmitancji oraz parametrów C i L. Używamy tych samych wzorów co w poprzednim podpunkcie.

$$\mathbf{G}_{\mathbf{u}}\mathbf{(s) =}\frac{\mathbf{1}}{\mathbf{7.995}\mathbf{*}\mathbf{10}^{\mathbf{-}\mathbf{7}}\mathbf{s}^{\mathbf{2}}\mathbf{+ 0.0001687s + 1}}$$

C1 = 1.8742*10⁻⁶ F

L = 0.4266 H

I po raz kolejny dla sprawdzenia porównujemy odpowiedzi skokowe układów.

4) Identyfikacja parametrów modelu obiektu MISO

Identyfikacja metodą MISO również polega na wykorzystaniu transformat Fouriera sygnałów wyjściowego i wejściowego. Następnie wyznaczamy gęstości widmowe mocy sygnałów wejściowego i zakłócającego oraz gęstości widmowe mocy wzajemnej osobno wejść (do obiektu i zakłócającego) i wyjścia. Potem wyznaczamy charakterystykę amplitudowo-fazową modelu obiektu oraz zakłóceń i odczytujemy z niej charakterystyczne wartości częstotliwości: ω2 dla Re = 0 oraz ω1 dla Re = - Im

u - wejście obiektu (sygnał losowy o rozkładzie równomiernym)

z - wejście zakłócające (sygnał losowy o rozkładzie normalnym)

y - wyjście obiektu.

Za pomocą funkcji gęstości widmowych mocy własnej i wzajemnej wyznaczam współczynniki transmitancji podobnie jak w po przednim punkcie. Model zakłócenia jest obiektem pierwszego rzędu odczytywany będzie jedynie punkt przecięcia dla części rzeczywistej równej części urojonej.

Po raz kolejny, podobnie jak w poprzednich punktach, odczytujemy charakterystyczne punkty z wykresów.

P1 = 32

P2 = 36

$$\mathbf{G}_{\mathbf{u}}\mathbf{(s) =}\frac{\mathbf{1}}{\mathbf{7.818}\mathbf{*}\mathbf{10}^{\mathbf{-}\mathbf{7}}\mathbf{s}^{\mathbf{2}}\mathbf{+ 0.0002088s + 1}}$$

C1 = 2.3196*10⁻⁶ F

L = 0.337 H

Wykres dla zakłóceń wygląda następująco:

P3 = 27

Obiekt zakłócający to obiekt RC dla którego transmitancja ma postać:

$$G_{z} = \frac{1}{R_{\text{z\ \ }}C_{z}*s + 1}$$

Za pomocą Matlaba obliczamy transmitancję obiektu zakłócającego

$$\mathbf{G}_{\mathbf{u}}\mathbf{(s) =}\frac{\mathbf{1}}{\mathbf{0.0001179s + 1}}$$

C1 = 1.1789*10⁻⁶ F

Dla sprawdzenia poprawności obliczeń:

5) Identyfikacja metodami parametrycznymi

IDENTYFIKACJA SISO

Za pomocą ToolBoxa wyznaczam transmitancje obiektów w oparciu o modele parametryczne :

- ARX (AutoRegressive with eXogenous input)

- IV (ARX with Instrumental Variables)

- ARMAX (AutoRegressive Moving Average with eXogenous input)

- OE (Output Error)

- BJ (Box-Jenkins)

Na początku ustawiamy w Toolboxie wszystkie potrzebne okienka, tzn: Time domain data (Input, Output, Data Name, Starting time, Sampling interval), następnie dodaje Polynominal Models.

Zanim zaimportujemy modele do Workspace i porównamy z odpowiedzią wzorcową za pomocą 'Model output' możemy dowiedzieć się w jakim stopniu się one pokrywają. Jeżeli chcemy wyznaczyć transmitancję z dużą dokładnością interesują nas te najbardziej pokrywające się modele. Pokażę jednak każdy z nich z osobna niezależnie od stopnia pokrycia.

Tak przedstawiają się wygenerowane odpowiedzi skokowe w porównaniu z odpowiedzią wzorcową:

Metoda BJ Metoda OE

Metoda IV Metoda ARX

Metoda ARMAX

Od razu widać, że najbardziej z odpowiedzią wzorcową pokrywają się metody BJ, OE oraz IV. Ich stopnie pokrycia są tak zbliżone, że patrząc na wykresy nie jesteśmy w stanie określić, która z nich jest najlepsza. Dzięki wcześniej zaprezentowanemu 'Model output' wiemy jednak, że najlepsza okazała się metoda BJ. Metod ARX i ARMAX w ogóle nie bierzemy pod uwagę - ich odpowiedzi najmniej pokrywają się z wzorcową.

Skoro najlepszy wynik uzyskała metoda BJ, zajmijmy się nią i wyliczmy transmitancję w sposób identyczny jak w punkcie pierwszym.

Z wykresu odpowiedzi skokowej odczytujemy potrzebne wartości:

ymax = 0.757

ymin = 0.5729

Tosc = 0.00589

I tak przedstawiają się interesujące nas wyniki obliczeń:

$$\omega = 1066,8\frac{\text{rad}}{s}$$

ζ = 0.0879

δ = 93, 39

$$\omega_{o} = 1070,8\frac{\text{rad}}{s}$$

T_o = 0.0009339s

$$\mathbf{G}_{\mathbf{u}}\mathbf{(s) =}\frac{\mathbf{1}}{\mathbf{8.721}\mathbf{*}\mathbf{10}^{\mathbf{-}\mathbf{7}}\mathbf{s}^{\mathbf{2}}\mathbf{+ 0.0001641s + 1}}$$

C1 = 1.8237*10⁻⁶F

L = 0.4782H

IDENTYFIKACJA MISO

Identyfikacja przebiega niemalże identycznie. Tym razem jednak bierzemy pod uwagę wpływ zakłóceń. W ustawieniach pojawi się więc dodatkowe wejście 'z'.

Wstępne wyniki pokazują nam o wiele lepsze pokrycie z odpowiedzią wzorcową. Dla metod BJ i OE wynosi ona aż 99,99%. Pozostałe metody pokrywają się w nie mniej niż 70 procentach, co jest lepszym wynikiem od najlepszego wyniku przy identyfikacji SISO.

Metoda BJ Metoda OE

Metoda IV Metoda ARX

Metoda ARMAX

Metody BJ oraz OE mają niemalże identyczne pokrycie z odpowiedzią wzorcową. W poprzednim punkcie badaliśmy wynik metody BJ, więc teraz zajmijmy się metodą OE.

Z wykresu odczytujemy następujące dane:

ymax = 0,727

ymin = 0,53

Tosc = 00565

I tak przedstawiają się interesujące nas wyniki obliczeń:

$$\omega = 1112,1\frac{\text{rad}}{s}$$

ζ = 0.1001

δ = 110, 76

$$\omega_{o} = 1117,6\frac{\text{rad}}{s}$$

T_o = 0.00089487s

Transmitancja Gu(s) przyjmuje postać:

$$\mathbf{G}_{\mathbf{u}}\mathbf{(s) =}\frac{\mathbf{1}}{\mathbf{8.007}\mathbf{*}\mathbf{10}^{\mathbf{-}\mathbf{7}}\mathbf{s}^{\mathbf{2}}\mathbf{+ 0.0001791s + 1}}$$

Możemy teraz wyliczyć C1 oraz L

C1 = 1.9904*10⁻⁶F

L = 0.4023H

Podsumowanie:

Naszym zadaniem była identyfikacja obiektu, jego transmitancji i parametrów. Wykorzystaliśmy do tego aż pięć metod. Teraz czas na podsumowanie.

Najprostszą metodą była identyfikacja na podstawie odpowiedzi skokowej. Dzięki możliwości bardzo dokładnego odczytania z wykresu interesujących nas wartości okazała się być dokładna (co obala mit o tym, że najprostsze metody są najgorsze i najmniej dokładne). Metoda jednak ma też swoją wadę - nie uwzględnia wpływu zakłóceń na obiekt.

Posługując się wcześniejszą metodą poszliśmy dalej i dokonaliśmy identyfikacji obiektu na podstawie odpowiedzi impulsowej. Efekt był bardzo podobny. Porównując wykresy odpowiedzi wzorcowej i odpowiedzi zidentyfikowanych widzimy bardzo podobny wynik (wizualnie możemy jednak stwierdzić, że druga metoda okazała się minimalnie lepsza).

Identyfikacja obiektu SISO przy wykorzystaniu funkcji gęstości widmowej mocy jest mniej dokładna. Tutaj również nie było uwzględnione zakłócenie.

Identyfikacja obiektu MISO przy wykorzystaniu funkcji gęstości widmowej mocy pozwala zidentyfikować zarówno tor główny jak i zakłócenie. Wynik tej identyfikacji uznajemy za zadowalający.

Identyfikacja modelami parametrycznymi za pomocą System Identification Toolbox dała nam możliwość przeprowadzenia weryfikacji modelu z wykorzystaniem zaawansowanych algorytmów. Narzędzie to pozwoliło nam również sprawdzić stopień dopasowania odpowiedzi wzorcowej i zidentyfikowanej. Zarówno w przypadku SISO jak i MISO najbardziej pasującym modelem okazał się BJ. W przypadku MISO równie dobrze wypadł model OE, gdzie pokrycie z odpowiedzią wzorcową wyliczone zostało na 99%.