STABILNOŚĆ LINIOWYCH UKŁADÓW AUTOMATYKI
1. Wprowadzenie do ćwiczenia
Problem stabilności każdego układu automatyki jest problemem podstawowym. Pojawienie się niestabilności w układach zazwyczaj pracujących stabilnie, może spowodować duże straty materialne, a nawet może zagrażać życiu ludzkiemu (np. sterowanie procesów o wielkiej mocy) i grozić uszkodzeniem urządzeń i maszyn.
Przez stabilność układu automatyki rozumiemy właściwość układu polegająca na tym, że układ wytrącony ze stanu ustalonego przez wymuszenie lub zakłócenie, po ustaniu wymuszenia lub zakłócenia, wraca w czasie skończonym do stanu ustalonego.
Warunki stabilności układu automatyki nie mogą zależeć od rodzaju pobudzania układu, to jest kształtu funkcji wielkości wejściowej lub zakłócenia, lecz tylko od struktury wewnętrznej układu.
Problem stabilności najprościej można wyjaśnić posługując się liniowym równaniem różniczkowym, opisującym przebieg wielkości wyjściowej y(t). Liniowy, jednowymiarowy układ dynamiczny opisany jest za pomocą równania różniczkowego (wzór 1).
u(t) y(t)
Rys. 1 Schemat blokowy układu dynamicznego
u b u b u b u b y a y a y a y an (n) + n 1 (n 1) +...+ 1 (1) + 0 = m (m) + m 1 (m−1) +...+ 1 (1) + 0 − − − , (1) czyli:
∑
∑
= = = n 0 i m 0 j j j j i i i dt ) t ( u d b dt ) t ( y d a , (2)gdzie: y(t) – wielkość wyjściowa,
u(t) – wielkość wejściowa,
i i i dt y d y() = dla (0,1,...,n), j j j dt u d u( )= dla (0,1, ...m), n≥m.
Rozwiązanie tego równania y(t), zgodnie z zasadą superpozycji, można przedstawić za
pomocą dwóch składowych:
- składowej swobodnej ys(t), zależnej tylko od parametrów układu,
- składowej wymuszonej yw(t), pochodzącej od wejścia wymuszającego.
Dla rozważanych układów liniowych można podać następującą definicję stabilności:
Układ jest stabilny, jeżeli dla dowolnych warunków początkowych składowa swobodna ys(t)
dąży do zera, gdy t dąży do nieskończoności.
Dla tak definiowanej stabilności, przy jej badaniu, wystarczy korzystać z równania
różniczkowego jednorodnego dynamiki układu (przy zerowym wymuszeniu):
0 ) ( ) ( ) ( 0 1 1 +⋅ ⋅⋅+ = + − − y t a y t a t y a n n n n (3) 2. Kryterium pierwiastkowe
Na podstawie definicji stabilności można określić warunki stabilności, korzystając z równania charakterystycznego układu automatyki (4), wynikającego bezpośrednio z równania (3):
0 0 1 1 1 +⋅ ⋅⋅+ + = + − − s a s a a s a n n n n . (4)
Kryterium to nazywane jest kryterium pierwiastkowym, ponieważ na podstawie pierwiastków równania charakterystycznego możemy określić stabilność układu. Pierwiastki równania charakterystycznego składają się z wyrażeń typu Ciesit, przy czym s
i=αi+jβi. W zależności
od parametrów α, β możemy rozpatrzyć następujące przypadki (tablica 1).
Tablica 1
Postać składowej przejściowej dla różnych pierwiastków równania charakterystycznego
Przypadek si=αi+jβi st
p i e t y ( )= 1. układ niestabilny (aperiodycznie) α >0; β=0 0 1 2 3 4 5 0 50 100 150 t y p (t )
2. układ stabilny (aperiodycznie) α <0; β=0 0 1 2 3 4 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t y p (t ) 3. układ niestabilny (oscylacyjnie) α >0; β≠0 0 1 2 3 4 5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 y p (t ) t 4. układ stabilny (oscylacyjnie) α<0; β≠0 0 1 2 3 4 5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 y p (t ) t 5. układ na granicy stabilności α =0; β≠0 0 1 2 3 4 5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 t y p (t )
Warunkiem koniecznym i wystarczającym stabilności układu regulacji jest, aby wszystkie pierwiastki równania charakterystycznego miały części rzeczywiste ujemne (przypadek 2, 4
w tablicy 1).
Re{si}=α<0 (5)
Przykład 1.1
Korzystając z kryterium pierwiastkowego wyznaczyć stabilność układu automatyki opisanego transmitancją Ko(s)
( )
26 s 27 s 9 s 3 1 s Ko 3 2 + + + =) s ( K 1 ) s ( K ) s ( K o o Z + =
Wstawiając Ko(s) do równania na Kz(s) otrzymujemy:
( )
27 s 27 s 9 s 3 1 s Kz 3 2 + + + =Równanie charakterystyczne układu jest następujące:
( )
s 3s 9s 27s 27 0Mz = 3+ 2 + + =
Do wyznaczenia pierwiastków równania charakterystycznego skorzystamy z polecenia roots
programu MATLAB.
W oknie poleceń MATLAB-a wpisujemy >> roots([3 9 27 27]) (wpisując współczynniki równania począwszy od najwyższego rzędu)
Równanie jest rzędu trzeciego otrzymujemy zatem trzy pierwiastki si. s1= -0.8363 + 2.4659i
s2= -0.8363 - 2.4659i s3= -1.3275
Zgodnie z warunkiem (5) sprawdzamy część rzeczywistą pierwiastków. Dla każdego z nich jest ona ujemna zatem układ ten jest stabilny.
Możemy także skorzystać z procedury pzmap. Na podstawie tego polecenia MATLAB-a otrzymamy wykres pierwiastków licznika i mianownika transmitancji, a zatem wykres zer i biegunów.
W oknie poleceń MATLAB-a wpisujemy >> pzmap([1],[3 9 27 27])
Rys. 2 Rozkład biegunów transmitancji Kz(s)
W praktyce stosowanie kryterium pierwiastkowego nie zawsze jest dogodne. W przypadku układów wyższego rzędu (tym samym wysokich rzędów równania charakterystycznego) pojawia się trudność w obliczeniu jego pierwiastków. Zostały opracowane inne metody pozwalające badać stabilność układu bez rozwiązywania równania charakterystycznego. Kryteria te dzielimy na analityczne i częstotliwościowe (analityczno-graficzne). Najbardziej znanym kryterium analitycznym jest kryterium Hurwitza. W wielu przypadkach nie jest znana postać transmitancji układu, natomiast znana jest (np. wyznaczona doświadczalnie)
charakterystyka częstotliwościowa układu. Najbardziej znanym kryterium
częstotliwościowym jest kryterium Nyquista (a także jego postać logarytmiczna –Bodego).
3. Kryterium Hurwitza
Jest to kryterium algebraiczne, które podaje warunki jakie powinny spełniać współczynniki równania charakterystycznego, aby pierwiastki tego równania miały części rzeczywiste ujemne (patrz kryterium pierwiastkowe).
Według kryterium Hurwitza muszą być spełnione dwa warunki nazywane warunkiem
koniecznym i wystarczającym:
1. wszystkie współczynniki równania charakterystycznego istnieją i są dodatnie ai>0 (warunek konieczny)
2. wyznacznik główny ∆n i wszystkie jego podwyznaczniki ∆i (i=1,2,...,n-1), utworzone z wyznacznika głównego, są większe od zera (warunek wystarczający).
Wyznacznik główny tworzymy wpisując na głównej przekątnej współczynniki równania charakterystycznego, począwszy od an-1 a skończywszy na a0. Następnie uzupełniamy
kolumny współczynnikami o kolejnych indeksach (względem współczynnika rozpatrywanego
na głównej przekątnej). Wpisujemy współczynniki o indeksach większych idąc w dół oraz o indeksach mniejszych idąc w górę, od współczynnika wpisanego na głównej przekątnej.
∆n n n n n n a a a a a a a a = − − − − 1 3 2 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 L L L M M L M M L L (6) Stąd podwyznacznik D1, ∆1 = an−1 (7) a podwyznacznik D2, ∆2 1 3 2 = − − − a a a a n n n n , (8)
analogicznie tworzymy pozostałe podwyznaczniki.
Przykład 1.2.
Dla transmitancji operatorowej Kz(s) z przykładu 1.1 wyznaczymy stabilność układu korzystając z kryterium Hurwitza
W jakim zakresie można zmienić k, aby spełnione były warunki stabilności układu ? Równanie charakterystyczne ma postać:
0 27 s 27 s 9 s 3 ) s ( Mz = 3 + 2 + + = .
1. Warunek konieczny jest spełniony (wszystkie współczynniki są większe od zera) 2. Warunek wystarczający dla Di >0
0 27 0 162 27 6 27 3 27 9 27 3 27 9 0 9 9 27 9 0 0 27 3 0 27 9 2 3 2 1 3 > ∆ ⋅ = ∆ > = ⋅ = ⋅ − ⋅ = = ∆ > = = ∆ = ∆
Wobec tego układ jest stabilny.
Wartość wyznacznika możemy obliczyć za pomocą procedury det.
W oknie MATLAB-a wpisujemy polecenia >> delta2=det([9 27;1 27]) oraz
delta3=det([9 27 0;1 27 0; 0 9 27]) Otrzymujemy wyniki: delta2 = ……..162 delta3 = 4374
A zatem układ jest stabilny.
4.1. Kryterium Nyquista
• Kryterium Nyquista pozwala badać stabilność zamkniętego układu regulacji na podstawie charakterystyki amplitudowo-fazowej układu otwartego.
• Metoda pozwala stwierdzić, na etapie projektu i budowy układu regulacji, czy po zamknięciu obwodu regulacyjnego układ będzie stabilny.
• Ważnym elementem kryterium Nyquista jest oparcie się na charakterystykach
częstotliwościowych, które mogą być wyznaczone doświadczalnie, a niekoniecznie na drodze analitycznej.
• Można określić stabilność układu także dla obiektów z opóźnieniem.
Eksperymentalne zdejmowanie charakterystyk amplitudowo-fazowych.
W celu określenia charakterystyki amplitudowo-fazowej zmieniamy częstotliwość sygnału sinusoidalnego wejściowego w1, w2 w3 ... i określamy parametry K( ωj ) oraz ϕ(ω).
K(jω)
x(t)=Xsinωt y(t)=XK(jω) sin(ωt+ϕ)
Im {K(jω)} Re {K(jω)} k ω=0 ω=∝ ω1 ω2 -ϕ(ω) K(jω) ω
Rys. 3. Eksperymentalne określanie charakterystyki częstotliwościowej
Jeżeli transmitancja widmowa układu otwartego jest równa Ko(jw), to transmitancja układu zamkniętego ma postać:
) ( 1 ) ( ) ( ω ω ω j K j K j K o o Z + = . (9)
Dla pewnej częstotliwości wr może zaistnieć przypadek, że będzie spełniony warunek
Ko(jwr)=-1.
0 ) (
1+Ko jω = (10)
Wówczas w równaniu (10) mianownik równa się zeru, co prowadzi do sytuacji, w której wyjściowy sygnał sinusoidalny może się sam podtrzymać, wywołując, poprzez sprzężenie zwrotne, na wejściu przebieg o tej samej amplitudzie i fazie (rys. 4). Sytuacja ta odpowiada układowi na granicy stabilności (rys.5a).
K(jωr)=-1 Y(ωr)sinωrt
x(t)=0 -Y(ωr)sinωrt
Rys. 4. Interpretacja kryterium Nyquista
Im {K(jω)} Re {K(jω)} ω=0 ω=ωr ω (-1,j0) a) Im {K(jω)} ω=0 ω (-1,j0) Re {K(jω)} b) Im {K(jω)} Re {K(jω)} ω=0 ω (-1,j0) c)
Rys. 5. Interpretacja wykresu charakterystyki częstotliwościowej Ko(jw) przy zastosowaniu
Według kryterium Nyquista w przypadku, gdy układ otwarty jest stabilny, będzie stabilny po zamknięciu (rys.5.c), jeżeli charakterystyka amplitudowo-fazowa układu otwartego nie obejmuje punktu (-1, j0).
W przypadku występowania bieguna zerowego w transmitancji Ko(s), (dotyczy to układów
zawierających człon całkujący) charakterystyka amplitudowo-fazowa Ko(jw) ma inny kształt
niż charakterystyki przedstawione na rys. 5. W tym przypadku, dla częstotliwości w0=0,
charakterystyka ( ) 1 *( ω )
ω
ω K j
j j
Ko = o rozpoczyna się na ujemnej półosi Im{Ko(jw)}
w nieskończoności (rys. 6). Nie ma to jednak wpływu na ocenę stabilności układu, gdyż także w tym przypadku kryterium Nyquista sprowadza się do nieobejmowania punktu (-1, j0) przez charakterystykę częstotliwościową Ko(jw). Im {K(jω)} Re {K(jω)} ω=0 ω=∝ (-1,j0)
Rys. 6. Przykład charakterystyki Ko(jw) z biegunem zerowym
Przykład 1.3.
Dla układu z przykładu 1.1 wyznaczyć stabilność korzystając z kryterium Nyguista.
W kroku pierwszym sprawdzamy czy układ otwarty jest stabilny. Korzystamy z kryterium Hurwitza. 0 26 0 26 3 27 9 27 3 26 9 0 9 9 26 9 0 0 27 3 0 26 9 2 3 2 1 3 > ∆ ⋅ = ∆ => ⋅ − ⋅ = = ∆ > = = ∆ = ∆
Układ otwarty jest stabilny. Wyznaczamy charakterystykę częstotliwościową Ko(jw).
Rys. 7. Charakterystyka Nyquista Ko(jw)
Układ będzie stabilny po zamknięcie jeżeli charakterystyka amplitudowo-fazowa układu otwartego nie będzie obejmowała punktu (-1;j0). Na podstawie rys. 7 widzimy, że układ po zamknięciu będzie stabilny. (Uwaga do rysunku. Polecenie MATLAB-a „nyquist” wyznacza charakterystykę Ko(jw) z przedziału częstotliwości (−∞;+∞)).
Ważną zaletą kryterium Nyquista jest łatwość określenia zapasu stabilności układu badanego. W pobliżu granicy stabilności stany nieustalone są oscylacyjne o tłumieniu coraz mniejszym, im bliżej granicy stabilności znajduje się układ. Dlatego istotnym jest, aby układ miał określony zapas stabilności. Zapas stabilności określa się za pomocą:
- zapasu modułu, - zapasu fazy.
Zapas modułu określa krotność o jaką musiałoby wzrosnąć wzmocnienie przy niezmienionym argumencie układu otwartego, aby układ zamknięty znalazł się na granicy stabilności.
Zapas fazy określa wartość zmiany argumentu transmitancji widmowej układu otwartego przy niezmiennym wzmocnieniu, która doprowadziłaby układ zamknięty do granicy stabilności. Pomocne w określeniu zapasu stabilności są analityczne warunki stabilności kryterium Nyquista, wynikające z równania (10):
− = < π ϕ ω) 1 ( j Ko (11a) lub
− > = 1 )} ( Re{ 0 )} ( Im{ ω ω j K j K o o . (11b)
4.2. Kryterium logarytmiczne Nyquista
Warunki stabilności wynikają bezpośrednio z warunku 11a i można je przedstawić następująco:
− = < . 0 ) ( log 20 π ϕ ω j Ko (12)
Po wyznaczeniu charakterystyk logarytmicznych sprawdzamy dwa warunki
częstotliwościowe:
- częstotliwość, przy której charakterystyka fazy przecina wartość −π (czyli –180o) i wtedy moduł powinien być ujemny,
- częstotliwość, przy której charakterystyka modułu przecina 0 dB i wtedy kąt przesunięcia fazowego powinien być większy niż −π (–180o).
Charakterystyka Nyquista przedstawiona w postaci logarytmicznej nazywana jest charakterystyką Bodego.
Aby wyznaczyć charakterystykę logarytmiczna amplitudy i fazy posłużymy się poleceniem
bode. W oknie MATLAB-a wpisujemy << bode([1],[3 9 27 26])
4.3. Zapas stabilności (modułu i fazy)
Zapas stabilności możemy wyznaczyć korzystając z polecenia margin. W oknie MATLAB-a wpisujemy << margin([1],[3 9 27 26])
Rys. 9. Zapas stabilności:
Poprawność otrzymanych wyników możemy sprawdzić analitycznie . Z charakterystyki przedstawionej na rys. 9 wynika, że dla pulsacji w=3 rad/s kąt przesunięcia fazowego jest równy -180 stopni. Wstawiamy do wyrażenia na Ko(jw) wartość w=3.
55 1 ) 81 81 ( j 81 26 1 26 3 j 27 ) 3 j ( 9 ) 3 j ( 3 1 ) 3 j ( Ko 3 2 =− − + − = + ⋅ + ⋅ + ⋅ =
Uzyskujemy tylko część rzeczywistą ujemną co odpowiada przesunięciu fazowemu o – 180 stopni. Moduł transmitancji widmowej jest równy
55 1 ) 3 j ( Ko =
Zgodnie z warunkiem (12) wyznaczamy
0 8 , 34 ) 3 j ( K log 20 o =− <
Zapas modułu wynosi 34,8 dB. Zapasu fazy nie musimy wyznaczać ponieważ moduł z transmitancji widmowej nie osiągnie wartości 1 (0 dB) dla żadnej częstotliwości.
5. Przebieg ćwiczenia laboratoryjnego
I. W programie MATLAB-Simulink utworzyć model układu regulacji jak na rysunku.
Rys. 10 Schemat blokowy układu regulacji
Wykorzystać bloki Simulinka: Transfer Fcn, PID Controller, Sum, Step oraz Scope. Struktura układu regulacji:
- obiekt opisany transmitancją drugiego rzędu - regulator PI.
Dla podanych, przez prowadzącego, parametrów obiektu przeprowadzić badania symulacyjne układu regulacji dla różnych nastaw parametrów regulatora PI.
Dobrać parametry regulatora, by uzyskać przebiegi aperiodyczne i oscylacyjne.
Dobrać parametry regulatora Kp oraz Ki tak aby uzyskać przebiegi dla układu: stabilnego, na granicy stabilności, niestabilnego.
II. dokonać oceny stabilności układu regulacji stosując różne kryteria stabilności. Skorzystać z poleceń \MATLAB-a omówionych we wprowadzeniu do ćwiczenia.
III. W sprawozdaniu należy zamieścić:
Wykresy charakterystycznych przebiegów regulacji uzyskanych w SIMULINK-u;
Obliczenia wykonane w MATLABie (roots, det, zapas modułu i zapas fazy).
Charakterystyki częstotliwościowe Nyquista i Bodego dla wybranych parametrów układu.
Obliczenia wykonane na podstawie wzorów, a dotyczące kryteriów Hurwitza, Nyquista, logarytmicznego Nyquista (Bodego).
Załącznik 1 do ćwiczenia - mplik MATLABA
% stabilnosc laboratorium 2011 % obiekt drugiego rzedu % kob(s)=1/(as^2+bs+c) % regulator PI % kr(s)=(kp+ki/s)=(kp*s+ki)/s % uklad otwarty % k0(s)=kob(s)*kr(s) clc; clear all;
% Przykladowe parametry obiektu k=1; a2=0.35; a1=2.4; a0=1; num1=k; den1=[a2 a1 a0]; kob=tf([num1],[den1]); % Parametry regulatora kp=1; ki=1; num2=[kp ki]; den2=[1 0]; kr=tf(num2,den2); ko=series(kob,kr); kz=feedback(ko,1); [numkz, denkz]=tfdata(kz,'v'); % kryterium pierwiastkowe p_rch=roots(denkz) n=length(denkz);
wyz_2=[denkz(2) denkz(4);denkz(1) denkz(3)];
wyz_3=[denkz(2) denkz(4) 0;denkz(1) denkz(3) 0;0 denkz(2) denkz(4)]; % kryterium Hurwitza k_hur2=det(wyz_2) k_hur3=det(wyz_3) subplot(2,2,1); % kryterium Nyquista nyquist(ko); subplot(2,2,2); % kryterium Bode bode(ko); grid; subplot(2,2,3); % zapas modulu i fazy margin(ko); subplot(2,2,4); % bieguny i zera pzmap(kz); figure; % przebiegi czasowe
% odpowiedz na wymuszenie sinusoidalne i wymuszenie stale t=0:0.01:20; y=step(kz,t); w=1; U=10; u=U*sin(w*t); [ysim x]=lsim(kz,u,t); subplot(2,1,1); plot(t,ysim,t,u); grid; subplot(2,1,2); plot(t,y); grid;
Załącznik 2 do ćwiczenia - mplik MATLABA
% stabilnosc laboratorium 2012 % obiekt drugiego rzedu % kob(s)=1/((T1s+1)(T2s+1)) % regulator PI % kr(s)=(kp+ki/s)=(kp*s+ki)/s % uklad otwarty % k0(s)=kob(s)*kr(s) %clc; clear all;
% Przykladowe parametry obiektu k=1; T1=10; T2=1; % Parametry regulatora kp=1; ki=1; liczkz=[kp*k k*ki]; miankz=[T1*T2 T1+T2 k*kp+1 k*ki]; % kryterium pierwiastkowe p_rch=roots(miankz) n=length(miankz);
wyz_2=[miankz(2) miankz(4);miankz(1) miankz(3)];
wyz_3=[miankz(2) miankz(4) 0;miankz(1) miankz(3) 0;0 miankz(2) miankz(4)]; % kryterium Hurwitza k_hur2=det(wyz_2) k_hur3=det(wyz_3) % kryterium Nyquista w=0.01:0.01:100; dlug=length(w); % charakterystyka A-F for i=1:dlug ko(i)=(k*ki+j*k*kp*w(i))/(-(T1+T2)*w(i)^2+j*(w(i)-T1*T2*w(i)^3)); end; % charakterystyka Bode mod_ko=abs(ko); M_ko=20*log10(mod_ko); for i=1:dlug if real(ko(i))<0 wsk_r(i)=1; else wsk_r(i)=0; end if imag(ko(i))>0 wsk_i(i)=1; else wsk_i(i)=0; end if (wsk_r(i)+wsk_i(i))==2 kat_ko(i)=-360+angle(ko(i))*180/pi; else kat_ko(i)=angle(ko(i))*180/pi; end end % wykresy figure; plot(ko); grid; figure; subplot(2,1,1); plot(log10(w),M_ko); ylabel('M[dB]'); xlabel('log10w'); grid;
subplot(2,1,2); plot(log10(w),kat_ko); grid;
ylabel('kat[%]'); xlabel('log10w');
