Wykład
Komputerowe wspomaganie projektowania
CAD
Projektowanie układów automatyki – MATLAB
Robert Kazała
26 czerwiec 2004
ver. 1.01
Wykład
Komputerowe wspomaganie projektowania
CAD
Projektowanie układów automatyki – MATLAB
Robert Kazała
26 czerwiec 2004
ver. 1.01
Wstęp
Matlab jest programem przeznaczonym do wykonywania różnorodnych obliczeń
numerycznych. Na całość pakietu składają się następujące elementy:
−
interpretator języka programowania wraz z bibliotekami podstawowych działań i obliczeń
na macierzach,
−
standardowe biblioteki procedur napisanych w języku programu Matlab – realizują one m.
in. przekształcenia macierzy, obliczenie wartości funkcji elementarnych i specjalnych,
całkowanie numeryczne, rozwiązywanie układów równań różniczkowych zwyczajnych,
−
biblioteki dodatkowe „toolbox” zawierające procedury wspomagające obliczenia
numeryczne wykonywane w różnych zastosowaniach m. in. w projektowaniu układu
sterowania, cyfrowym przetwarzaniu sygnału, metodach identyfikacji i optymalizacji, w
sieciach neuronowych i wielu innych,
−
nakładki – są to dodatkowe programy napisane w języku programu Matlab ułatwiające
realizację obliczeń określonego rodzaju. Przykładem takiej nakładki jest Simulink
umożliwiający interakcyjne definiowanie struktury układu sterowania oraz wygodną jego
symulację.
Control System Toolbox jest biblioteką narzędzi przeznaczonych do modelowania,
analizy i syntezy układów liniowych dynamicznych ciągłych i dyskretnych niezmiennych w
czasie oraz mieszanych dyskretno-ciągłych.
Układy te mogą być opisane za pomocą:
•
transmitancji,
•
iloczynu zer, biegunów i wzmocnienia
•
jak również w postaci równań stanu.
Funkcje Control System Toolbox umożliwiają:
•
tworzenie obiektów LTI ciągłych i dyskretnych w dowolnej z zdefiniowanych postaci
•
dokonanie wzajemnej konwersji różnych postaci modeli ciągłych i dyskretnych,
•
upraszczanie modeli zbudowanych z różnych bloków poprzez łączenie ich w różny
sposób,
•
wyznaczanie charakterystyk czasowych
•
wyznaczanie charakterystyk częstotliwościowych
•
określanie parametrów dynamicznych modeli
•
przekształcanie modeli opisanych w przestrzeni stanów
•
projektowanie regulatorów metodą miejsc geometrycznych pierwiastków
•
projektowanie estymatorów i regulatorów optymalnych dla modeli opisanych w
przestrzeni stanów
1. Tworzenie modeli obiektów liniowych ciągłych i dyskretnych
Układem liniowym można nazwać układ, w którym spełniona została zasada
superpozycji – odpowiedź na sumę wymuszeń jest równa sumie odpowiedzi na każde
wymuszenie z osobna. Układ liniowy ze względu na charakter występujących w nim
sygnałów można podzielić na:
- układ ciągły, w którym sygnały przesyłane są w sposób ciągły a sam układ dynamiczny
można opisać równaniami różniczkowymi,
- układ dyskretny, w którym przynajmniej jeden sygnał ma charakter dyskretny a sam
układ dynamiczny można opisać równaniami różnicowymi.
Układ liniowy w systemie MATLAB rozumiany jest jako obiekt, niżej oznaczany jako SYS i
możemy wykorzystać różne funkcje do jego tworzenia.
1.1. Modele opisane transmitancją
Najpopularniejszym sposobem opisu układów liniowych jest opis za pomocą transmitancji
będącej ilorazem wielomianów licznika i mianownika. W przypadku obiektu
jednowymiarowego opisuje się go za pomocą jednego równania. W przypadku układów
wielowymiarowych konieczne jest utworzenie macierzy transmitancji opisującej zależności
pomiędzy każdym wejściem i każdym wyjściem.
1.1.1. Obiekty SISO z jednym wejściem i jednym wyjściem
Obiekty opisane transmitancją
)
(
)
(
)
(
s
M
s
L
s
G
=
Do utworzenia transmitancji należy wykorzystać komendę tf
•
h = tf(num,den) tworzy układ ciągły
•
h = tf(num,den,T) tworzy układ dyskretny z czasem próbkowania T.
Przykład Tworzenie układu o transmitancji
2
2
15
9
)
(
2
3
2
+
−
−
+
+
=
s
s
s
s
s
s
G
W tym celu należy wykonać następujące wywołania:
num=[1 9 15];
den=[1 -2 -1 2];
h=tf(num, den)
lub
h=tf([1 9 15], [1 -2 -1 2])
Obiekty można także tworzyć poprzez zdefiniowanie specjalnej zmiennej s (od wersji 5.3)
s = tf('s')
i podanie transmitancji postaci
H=s/(s^2+2*s+10)
Zmienna s jest tworzona tylko raz.
Obiekty MIMO z wieloma wejściem i wieloma wyjściami
h11=tf([-1 1],[1 1]);
h21=tf([1 2],[1 4 5]);
lub
s=tf('s')
h11=(s-1)/(s+1);
h21=(s+2)/(s^2+4*s+5);
macierz wynikowa ma postać
H=[h11; h21]
Transmitancje można też tworzyć poprzez podanie tablic współczynników
N = {[1 –1];[1 2]};
D = {[1 1];[1 4 5]};
H=tf(N,D)
1.2. Opis poprzez podanie wzmocnień, zer i biegunów
ZPK ( zero-pole-gain ) układ opisany przez iloczyn zer, biegunów i wzmocnień.
SYS = ZPK(Z,P,K) tworzy układ ciągły.
SYS = ZPK(Z,P,K,T) tworzy układ dyskretny z czasem próbkowania T.
h=zpk(0,[1-i 1+i 2],-2)
tworzy
-2s
---------------
(s-2)(s^2-2s+2)
Korzystając z operatora s
s=zpk('s')
H=-2s/((s-2)*(s^2+2*s+2));
Obiekty MIMO
Z = {[],-5;1-i 1+i] []};
P = {0, [-1 –1];[1 2 3],[]}
K = [-1 3;2 0];
H=zpk(Z,P,K)
tworzy macierz
−
−
−
+
−
+
+
−
=
0
)
3
(
2
)(
1
(
)
2
2
(
2
)
1
(
)
5
(
3
1
)
(
2
s
s
s
s
s
s
s
s
s
G
1.3. Opis równaniem stanu
SS ( state-space ) układ opisany za pomocą równań stanu.
SYS = SS(A,B,C,D) tworzy układ ciągły.
SYS = SS(A,B,C,D,T) tworzy układ dyskretny z czasem próbkowania T.
dla macierzy
−
−
=
2
5
1
0
A
,
=
3
0
B
,
[
]
1
0
=
C
,
0
=
D
sys = ss([0 5;-5 –2],[0;3],[0 1],0)
1.4. Opis w postaci filtru
FILT układ dyskretny opisany w formacie DSP :
( )
( )
( )
1
1
1
−
−
−
=
z
DEN
z
NUM
z
H
(4.1.1)
H = FILT(NUM,DEN) tworzy układ dyskretny bez uwzględnienia czasu próbkowania.
H = FILT(NUM,DEN,TS) tworzy układ dyskretny z uwzględnieniem czasu próbkowania
TS.
2
1
1
1
3
2
1
5
.
0
1
)
(
−
−
−
−
+
+
+
=
z
z
z
z
H
h=filt([1 0.5],[1 2 3])
1.5. Łączenie obiektów
W celu uproszczenia obliczeń oraz zwiększenia przejrzystości schematu blokowego
można z minimum dwóch układów złożyć jeden obiekt. Funkcje systemu MATLAB
przedstawione poniżej, pozwalają na złączenie układów połączonych szeregowo, równolegle,
pętlą sprzężenia zwrotnego lub w gwiazdę, przy czym można łączyć zarówno całe wektory
wejść i wyjść jak i ich wybrane elementy. Złożenie układów jest możliwe nawet wtedy, gdy
są one opisane różnymi metodami.
APPEND łączy dwa systemy w jeden większy przez odpowiednie łączenie wektorów wejść i
wyjść jako elementów o różnych składowych
SYSOUT = APPEND(SYS1,SYS2, ...) wywołanie to jest jednoznaczne z powiązaniem
wektorów wejść i wyjść układów SYS1, SYS2,... w następujący sposób:
=
SysN
Sys
Sys
SysOut
...
0
0
...
....
...
...
0
...
2
0
0
...
0
1
[A,B,C,D] = APPEND(A1,B1,C1,D1,A2,B2,C2,D2) wywołanie to dołącza do układu
opisanego równaniami stanu określonymi macierzami A1, B1, C1, D1 układ o równaniach
stanu wyznaczonych przez macierze A2, B2, C2, D2. Daje to łączne równanie stanu w
następującej postaci:
+
=
2
1
2
0
0
1
2
1
2
0
0
1
2
1
.
.
u
u
B
B
x
x
A
A
x
x
PARALLEL tworzy układ wypadkowy z równolegle połączonych dwóch układów
SYS = PARALLEL(SYS1,SYS2) wywołanie to tworzy system z dwóch układów
połączonych równolegle
SYS = PARALLEL(SYS1,SYS2,INP1,INP2,OUT1,OUT2) wywołanie to tworzy system
z dwóch układów SYS1 i SYS2 połączonych równolegle w taki sposób że: INP1 zawiera
indeksy wejść układu SYS1, które mają zostać podłączone do odpowiadających im wejść
układu SYS2, co określa w ten sam sposób wektor INP2 natomiast OUT1 zawiera indeksy
wyjść układu SYS1, które mają zostać zsumowane z odpowiadającymi im wyjściami
układu SYS2, określonymi przez wektor OUT2.
Rys. 4.1 Schemat blokowy równoległego połączenia dwóch układów.
SERIES tworzy układ wypadkowy z szeregowo połączonych dwóch układów
SYS = SERIES(SYS1,SYS2,OUT1,INP2) wywołanie to łączy dwa układy w szereg tak,
że wyjścia układu SYS1 określone przez OUT1 są połączone z wejściami układu SYS2
określonymi przez INP2 w sposób pokazany na rysunku.
Jeżeli OUT1 i INP2 są pominięte, SERIES łączy SYS1 i SYS2 w kaskadę i zwraca
SYS = SYS1 * SYS2 .
Rys. 4.2 Schemat blokowy połączenia szeregowego dwóch układów.
FEEDBACK tworzy układ wypadkowy z układów połączonych w układ sprzężenia
zwrotnego.
SYS = FEEDBACK(SYS1,SYS2,SIGN) wywołanie to łączy układy w sprzężeniu
zwrotnym pokazane na rysunku. Jeśli SIGN =1 to sprzężenie zwrotne jest dodatnie,
natomiast jeśli SIGN = -1 to sprzężenie zwrotne jest ujemne. Pominięcie SIGN oznacza
użycie ujemnego sprzężenia zwrotnego.
Rys. 4.3 Schemat blokowy dwóch układów w sprzężeniu zwrotnym.
SYS = FEEDBACK(SYS1,SYS2,FEEDIN,FEEDOUT,SIGN) wywołanie to łączy układy
poprzez odpowiednie połączenie niektórych z ich wejść i wyjść w sposób pokazany na
Rys. 4.4.
Wektor FEEDIN zawiera indeksy wejść u układu SYS1 które są wykorzystywane przy
połączeniu sprzężeniem zwrotnym, podobnie wektor FEEDOUT zawiera indeksy wyjść y
układu SYS1 które są wykorzystywane przy połączeniu sprzężeniem zwrotnym z układem
SYS2.
Rys. 4.4 Schemat blokowy połączenia dwóch układów.
STAR tworzy układ wypadkowy z dwóch układów połączonych w gwiazdę
SYS = STAR(SYS1,SYS2,NU,NY) wywołanie to łączy dwa układy SYS1 i SYS2 w
jeden system w sposób pokazany na rys. 4.5.
Rys. 4.5 Schemat blokowy połączenia dwóch układów.
To sprzężenie zwrotne łączy pierwsze wyjście układu SYS2 z ostatnim wejściem układu
SYS1 (sygnał u) zawartymi w wektorze NU a ostatnie wyjście układu SYS1 z pierwszym
wejściem układu SYS2 (sygnał y) zawartymi w wektorze NY. W rezultacie otrzymujemy
układ SYS z wektorem wejść [w1;w2] i wektorem wyjść [z1;z2].
1.6. Metody konwersji obiektów
MATLAB określa układ za pomocą równań stanu, transmitancji oraz biegunów i zer. W
bibliotekach programu są opisane poniżej funkcje pozwalające na przekształcenie układu na
dowolny system zapisu, zamianę układu z ciągłego na dyskretny i odwrotnie przy
zastosowaniu różnych metod dyskretyzacji.
SS, ZPK, TF umożliwiają przekształcenie zapisu układu do określonej przez nie postaci:
równań stanu, zer i biegunów czy transmitancji.
Przykład: dany układ opisany przez zera, bieguny i wzmocnienie
)
5
)(
3
(
)
1
)(
2
(
2
−
−
−
+
=
s
s
s
s
SYS
należy przekształcić do postaci transmitancji SYS1 oraz opisać go równaniami stanu
SYS2.
W tym celu należy wykonać następujące wywołania:
>>
SYS = ZPK([-2 1],[3 5],2);
>> SYS1 = TF(SYS)
>> SYS2 = SS(SYS)
W wyniku ich realizacji otrzymamy transmitancję takiego układu i macierze równań
stanu:
Transfer function:
2 s^2 + 2 s - 4
---------------
s^2 - 8 s + 15
a = x1 x2
x1 8.00000 -3.75000
x2 4.00000 0
b = u1
x1 5.65685
x2 0
c = x1 x2
y1 3.18198 -1.50260
d = u1
y1 2.00000
C2D umożliwia przekształcenie układu ciągłego na dyskretny.
SYSD = C2D(SYSC,TS,METHOD) wywołanie to przekształca układ ciągły na dyskretny
z czasem próbkowania TS. Pole METHOD pozwala na dobór odpowiedniej metody
dyskretyzacji :
'zoh'
ekstrapolator zerowego rzędu,
'foh'
aproksymacja trójkątna (modyfikacja ekstrapolatora pierwszego
rzędu),
'tustin'
aproksymacja bilinearna używa przekształcenia:
1
1
2
'
+
−
=
z
z
T
s
s
(4.5.1)
'prewarp'
aproksymacja Tustina z uprzednim przekształceniem
częstotliwości:
(
)
1
1
2
'
+
−
=
z
z
T
tan
s
s
ω
ω
.
(4.5.2)
Częstotliwość krytyczna WC jest podana jako ostatnia w
wywołaniu C2D(SYSC,TS,'prewarp',WC),
'matched'
dopasowany do metody zer i biegunów tylko dla układów SISO.
Jeżeli nie zostanie wybrana żadna metoda program zastosuje ekstrapolator zerowego
rzędu.
D2C umożliwia przekształcenie układu dyskretnego na ciągły.
SYSC = D2C(SYSD,METHOD) wywołanie to przekształca układ dyskretny na ciągły z
czasem próbkowania TS. Pole METHOD pozwala na dobór odpowiedniej metody
dyskretyzacji :
'zoh'
ekstrapolator zerowego rzędu.
'foh'
aproksymacja trójkątna (modyfikacja ekstrapolatora pierwszego
rzędu)
'tustin'
aproksymacja bilinearna używa przekształcenia:
2
/
1
2
1
'
s
s
sT
sT
z
−
+
=
(4.5.3)
'prewarp'
aproksymacja Tustina z uprzednim przekształceniem
częstotliwości. Częstotliwość krytyczna WC jest podana jako
ostatnia w wywołaniu C2D(SYSD,TS,'prewarp',WC).
'matched'
dopasowany do metody zer i biegunów tylko dla układów SISO
Jeżeli nie zostanie wybrana żadna metoda program zastosuje ekstrapolator zerowego
rzędu.
D2D umożliwia zmianę czasu próbkowania lub dodanie opóźnienia na wejście.
SYS = D2D(SYS,TS) wywołanie to przepróbkowywuje układ dyskretny w celu
utworzenia układu dyskretnego z nowym czasem próbkowania TS.
SS2SS zmienia współrzędne stanu układu przestrzeni stanu
SYS = SS2SS(SYS,T) wywołanie to wykonuje podstawienie
Tx
x
=
w równaniach stanu
w wyniku czego następuje odpowiednie ich przekształcenie:
Du
Cx
y
Bu
Ax
x
+
=
+
=
&
→
=
Tx
x
Du
x
CT
y
TBu
x
TAT
x
+
=
+
=
−
−
1
1
&
(4.5.4)
CANON realizuje postać kanoniczną przestrzeni stanu.
CSYS = CANON(SYS,TYPE) wywołanie to oblicza postać kanoniczną przestrzeni stanu.
Pole TYPE określa postać kanoniczną :
'modal' – modalna postać kanoniczna gdzie wartości własne układu ukazane są na
przekątnej, macierz stanu A musi być diagonalna
'companion' – pierwsza postać kanoniczna gdzie wartości charakterystyczne
wielomianu ukazane są w prawej kolumnie
1.7. Odczytywanie informacji o obiektach
Układy, na których wykonuje się obliczenia podlegają różnym przekształceniom, w związku z
tym zmieniają się ich parametry. Zależnie od sposobu zapisu układu można użyć
przedstawionych niżej funkcji aby odczytać równania układu, jego bieguny, zera i
wzmocnienia czy wartości własnych. MATLAB dzięki funkcjom GET i SET pozwala też na
odczyt i modyfikację poszczególnych parametrów bez wykonania dodatkowych obliczeń.
SSDATA umożliwia odtworzenie macierzy opisujących system
[A,B,C,D] = SSDATA(SYS) wywołanie to zwraca wartości z macierzy A, B, C, D. Jeśli
SYS nie jest modelem opisanym równaniami stanu to następuje jego konwersja do postaci
przestrzeni stanu.
[A,B,C,D,TS,TD] = SSDATA(SYS) wywołanie to zwraca również czas próbkowania TS i
opóźnienie TD na wejściu. Dla układów ciągłych TD jest wektorem wierszowym, w
którym element TD(j) określa ile sekund opóźnienia ma j-te wejście. Dla układów
dyskretnych TD jest macierzą pustą.
ZPKDATA umożliwia odtworzenie zer, biegunów i wzmocnienia układu.
[Z,P,K] = ZPKDATA(SYS) wywołanie to zwraca zera, bieguny i wzmocnienie z każdego
kanału I/O układu.
[Z,P,K,TS,TD] = ZPKDATA(SYS) wywołanie to zwraca również czas próbkowania TS i
opóźnienia na wejściu, dla układów ciągłych TD jest wektorem wierszowym, w którym
element TD(j) określa ile sekund opóźnienia ma j-te wejście a dla układów dyskretnych
TD jest macierzą pustą.
[Z,P,K] = ZPKDATA(SYS,'V') wywołanie to zwraca zera i bieguny jako wektory
kolumnowe z podaniem ich wartości.
TFDATA umożliwia odtworzenie licznika i mianownika transmitancji opisującej system.
[NUM,DEN] = TFDATA(SYS) wywołanie to zwraca licznik i mianownik transmitancji.
[NUM,DEN,TS,TD] = TFDATA(SYS) wywołanie to zwraca również czas próbkowania
TS i opóźnienia na wejściu, dla układów ciągłych TD jest wektorem wierszowym, w
którym element TD(j) określa ile sekund opóźnienia ma j-te wejście a dla układów
dyskretnych TD jest macierzą pustą.
[NUM,DEN] = TFDATA(SYS,'V') wywołanie to zwraca licznik i mianownik jako
wektory kolumnowe z podaniem ich wartości.
SET wyświetla i modyfikuje własności układu.
SET(SYS) wywołanie to wyświetla wszystkie własności układu
SET(SYS,'PROPERTY') wywołanie to wyświetla możliwe wartości tej własności układu
SET(SYS,'PROPERTY',VALUE) wywołanie to przypisuje nową wartość VALUE
własności układu podanej w polu PROPERTY.
SET(SYS,'PROPERTY1',VALUE1,'PROPERTY2',VALUE2,...) wywołanie to
przypisuje podane wartości odpowiednim własnościom układu w jednym wywołaniu.
GET odczytuje zadane własności układu.
VALUE = GET(SYS,'PROPERTY') wywołanie to podaje wartość wybranej własności
PROPERTY układu SYS.
GET(SYS) wywołanie to wyświetla wszystkie własności układu z podaniem ich wartości.
1.8. Generowanie przykładowych modeli
System MATLAB umożliwia utworzenie przypadkowych, stabilnych układów ciągłych
bądź dyskretnych jeżeli zachodzi potrzeba przetestowania metod obliczeniowych na układzie
określonego rzędu.
DRSS tworzy przypadkowy stabilny i dyskretny układ opisany równaniami stanu.
SYS = DRSS(N) wywołanie to tworzy N-rzędowy układ SISO.
SYS = DRSS(N,P) wywołanie to tworzy N-rzędowy jednowejściowy układ z wyjściami
P.
SYS = DRSS(N,P,M) wywołanie to tworzy N-rzędowy układ z wyjściami P i
wejściami M.
DRMODEL tworzy przypadkowy stabilny i dyskretny N-rzędowy układ testowy.
[NUM,DEN]=DRMODEL(N) wywołanie to tworzy N-rzędowy układ SISO opisany
transmitancją.
[NUM,DEN]=DRMODEL(N,P) wywołanie to tworzy N-rzędowy jednowejściowy
układ z wyjściami P opisany transmitancją.
[A,B,C,D]=DRMODEL(N) wywołanie to tworzy N-rzędowy układ SISO opisany
równaniami stanu.
[A,B,C,D]=DRMODEL(N,P,M) wywołanie to tworzy N-rzędowy układ z wyjściami P
i wejściami M opisany równaniami stanu.
RSS tworzy przypadkowy, stabilny i ciągły układ opisany równaniami stanu
SYS = RSS(N) wywołanie to tworzy N-rzędowy układ SISO.
SYS = RSS(N,P) wywołanie to tworzy N-rzędowy jednowejściowy układ z wyjściami
P.
SYS = RSS(N,P,M) wywołanie to tworzy N-rzędowy układ z wyjściami P i wejściami
M.
RMODEL przypadkowy stabilny i ciągły N-rzędowy układ testowy.
[NUM,DEN]=RMODEL(N) wywołanie to tworzy N-rzędowy układ SISO opisany
transmitancją.
[NUM,DEN]=RMODEL(N,P) wywołanie to tworzy N-rzędowy jednowejściowy
układ z wyjściami P opisany transmitancją.
[A,B,C,D]=RMODEL(N) wywołanie to tworzy N-rzędowy układ SISO opisany
równaniami stanu.
[A,B,C,D]=RMODEL(N,P,M) wywołanie to tworzy N-rzędowy układ z wyjściami P i
wejściami M opisany równaniami stanu.
ORD2 układ oscylacyjny 2-go rzędu
[A,B,C,D] = ORD2(WN, Z) wywołanie to zwraca A, B, C, D reprezentujące układ
oscylacyjny 2-go rzędu gdzie WN jest częstotliwością rezonansową, a Z
współczynnikiem tłumienia.
[NUM,DEN] = ORD2(WN,Z) wywołanie to zwraca wielomiany reprezentujące
licznik i mianownik układu oscylacyjnego 2-go rzędu.
2. Charakterystyki modeli
2.1. Określanie właściwości układu
EIG, POLE znajduje bieguny układu.
P = EIG(SYS) lub P = POLE(SYS) wywołanie to zwraca bieguny układu jako wektor
kolumnowy P.
TZERO znajduje zera układu.
Z = TZERO(SYS) wywołanie to zwraca zera układu jako wektor kolumnowy Z.
[Z,GAIN] = TZERO(SYS) wywołanie to dla układu SISO zwraca zera i podaje jego
wzmocnienie.
Z = TZERO(A,B,C,D) wywołanie to zwraca zera pracując bezpośrednio na
macierzach przestrzeni stanu.
PZMAP oblicza bieguny i zera dla układów liniowych
PZMAP(SYS) dla układów SISO oblicza bieguny i zera układu i kreśli na
płaszczyźnie zmiennej zespolonej bieguny oznaczone krzyżykami oraz zera
oznaczone kółkami natomiast w przypadku układu MIMO kreśli bieguny i przekazuje
zera układu.
[P,Z] = PZMAP(SYS) to wywołanie umieszcza bieguny i zera układu w wektorach
kolumnowych P i Z, nie wywołuje natomiast wykresu.
Funkcje SGRID dla układu ciągłego lub ZGRID dla układu dyskretnego mogą być użyte
do nakreślenia na wykresie linii stałego współczynnika tłumienia i stałej częstotliwości
własnej.
DCGAIN określa wzmocnienie układu.
K = DCGAIN(SYS) wywołanie to oblicza wzmocnienie układu w stanie ustalonym.
NORM oblicza normy wektorów i macierzy.
NORM(SYS) wywołanie to oblicza normę H
2
dla układu ciągłego lub dyskretnego,
norma ta jest nieskończona jeśli układ jest niestabilny lub układ jest ciągły i ma
niezerowe wzmocnienie przy nieskończonej częstotliwości.
NORM(SYS,2) wywołanie to jest równoważne NORM(SYS).
NORM(SYS,INF) wywołanie to oblicza normę L
∞
układu, norma ta jest nieskończona
jeśli układ ciągły ma bieguny w urojonej części płaszczyzny zespolone.
NORM(SYS,INF,TOL) wywołanie to oblicza normę L
∞
układu z żądaną względną
dokładnością.
[N,FPEAK] = NORM(SYS,INF) wywołanie to zwraca również częstotliwość FPEAK
dla wzmocnienie o szczytowej wartości.
DAMP oblicza współczynniki tłumienia i pulsację drgań własnych.
[WN,Z] = DAMP(SYS) wywołanie to zwraca wektor WN i Z zawierające
współczynniki tłumienia i pulsację drgań własnych układu. WN i Z są pustymi
wektorami jeśli czas próbkowania TS jest niezdefiniowany.
[WN,Z,P] = DAMP(SYS) wywołanie to zwraca również bieguny układu w wektorze P
ESORT sortuje malejąco zespolone wartości własne.
S = ESORT(P) wywołanie to sortuje malejąco części rzeczywiste zespolonych
wartości własnych umieszczonych w wektorze P, niestabilne wartości własne
wpisywane są jako pierwsze.
[S,NDX] = ESORT(P) wywołanie to zwraca również wektor NDX zawierający
indeksy sortowanych wartości własnych.
DSORT sortuje malejąco zespolone dyskretne wartości własne.
S = DSORT(P) wywołanie to sortuje malejąco ze względu na moduł zespolone
dyskretne wartości własne zawarte w wektorze P, niestabilne wartości własne
wpisywane są jako pierwsze.
[S,NDX] = DSORT(P) wywołanie to zwraca również wektor NDX zawierający
indeksy sortowanych wartości własnych.
2.2. Wyznaczanie odpowiedzi czasowych
Charakterystyki czasowe układu można określić jako charakterystyki w stanie
nieustalonym układu, będące odpowiedzią układu na zadany, stały w czasie, sygnał
wejściowy. Ze względu na rodzaj tego sygnału można rozróżnić odpowiedź skokową,
impulsową oraz odpowiedź na stan początkowy – przebiegi stanu i wyjścia jakie powstaną w
układzie przy braku pobudzenia z zewnątrz, dla danych warunków początkowych. Funkcje
przedstawione poniżej pozwalają wykreślić te charakterystyki jak również znaleźć odpowiedź
układu na utworzony przez użytkownika sygnał sinusoidalny, prostokątny bądź okresowo
impulsowy.
INITIAL wyznacza odpowiedź układu na zadane warunki początkowe.
INITIAL(SYS,X0) wywołanie to kreśli odpowiedź układu na zadany warunek
początkowy X0, odpowiedź jest opisana przez następujące równanie:
Układ ciągły:
Ax
x
=
&
,
Cx
y
=
,
( )
0
0 x
x
=
Układ dyskretny:
x[k+1] = A x[k], y[k] = C x[k], x[0] = x0
Przedział czasu i liczba punktów są dobierane automatycznie.
INITIAL(SYS,X0,TFINAL) wywołanie to symuluje odpowiedź dla podanego czasu,
dla układu dyskretnego TFINAL określa liczbę kroków.
INITIAL(SYS,X0,T) wywołanie to wektor T określa czas symulacji odpowiedzi, dla
układu dyskretnego T ma postać Ti:Ts:Tf gdzie Ts określa czas próbkowania, dla
układu ciągłego T ma postać Ti:dt:Tf gdzie dt określa odstęp między kolejnym
próbkowaniem.
INITIAL(SYS1,SYS2,...,X0,T) wywołanie to kreśli odpowiedź impulsową dla wielu
układów na jednym wykresie z opcjonalnie użytym wektorem T, można też określić
kolor oraz styl linii dla każdego z układów np. initial(sys1,'r',sys2,'y--',sys3,'gx').
[Y,T,X] = INITIAL(SYS,X0,...) wywołanie to dla układu przestrzeni stanu zwraca
również wektor stanu X z kolejnych chwil symulacji.
Przykład: dla danego układu SYS opisanego równaniami stanu należy wyznaczyć
czasową na zadany sygnał X(0):
−
−
=
2
1
2
1
0
7814
.
0
7814
.
0
5572
.
0
x
x
x
x
&
&
[
]
=
2
1
4493
.
6
9691
.
1
x
x
y
=
0
1
)
0
(
x
W tym celu należy wykonać następujące wywołania:
>> A = [-0.5572 –0.7814 ; 0.7814 0];
>> C = [1.9691 6.4493];
>> X0 = [1 ; 0];
>> SYS=SS(A, [], C, [] );
INITIAL(SYS,X0)
W wyniku ich realizacji otrzymamy odpowiedź układu:
Rys. 4.6 Odpowiedź układu na zadane warunki początkowe.
STEP wyznacza odpowiedź układu na skok jednostkowy.
STEP(SYS) wywołanie to kreśli odpowiedź układu na skok jednostkowy podany na
wejście, przedział czasu i liczba punktów dobierane są automatycznie.
STEP(SYS,TFINAL) wywołanie to symuluje odpowiedź skokową dla zadanego czasu
TFINAL, natomiast dla układu dyskretnego bez ustalonego czasu próbkowania
TFINAL jest interpretowany jako liczba próbkowań.
Time (sec.)
A
m
pl
itu
de
Initial Condition Results
0
5
10
15
20
-2
-1
0
1
2
3
4
5
STEP(SYS,T) wywołanie to wektor T określa czas symulacji odpowiedzi skokowej,
dla układu dyskretnego T ma postać Ti:Ts:Tf gdzie Ts określa czas próbkowania, dla
układu ciągłego T ma postać Ti:dt:Tf gdzie dt jest czasem próbkowania z dyskretnej
aproksymacji na układ ciągły.
STEP(SYS1,SYS2,...,T) wywołanie to kreśli odpowiedź skokową dla wielu układów
na jednym wykresie z opcjonalnie użytym wektorem T, można też określić kolor oraz
styl linii dla każdego z układów np. (sys1,'r',sys2,'y--',sys3,'gx').
[Y,T] = STEP(SYS,...) wywołanie to zwraca wyniki symulacji w postaci wektora
odpowiedzi skokowej Y oraz wektora czasu trwania symulacji T.
[Y,T,X] = STEP(SYS,...) wywołanie to dla układu przestrzeni stanu zwraca również
wektor stanu X z kolejnych chwil symulacji.
IMPULSE wyznacza odpowiedź układu na impuls jednostkowy.
IMPULSE(SYS) wywołanie to kreśli odpowiedź układu na impuls jednostkowy
podany na wejście, przedział czasu i liczba punktów dobierane są automatycznie.
IMPULSE(SYS,TFINAL) wywołanie to symuluje odpowiedź impulsową dla
zadanego czasu TFINAL, natomiast dla układu dyskretnego bez ustalonego czasu
próbkowania TFINAL jest interpretowany jako liczba próbkowań.
IMPULSE(SYS,T) wywołanie to wektor T określa czas symulacji odpowiedzi
impulsowej, dla układu dyskretnego T ma postać Ti:Ts:Tf gdzie Ts określa czas
próbkowania, dla układu ciągłego T ma postać Ti:dt:Tf gdzie dt jest czasem
próbkowania z dyskretnej aproksymacji na układ ciągły.
IMPULSE(SYS1,SYS2,...,T) wywołanie to kreśli odpowiedź impulsową dla wielu
układów na jednym wykresie z opcjonalnie użytym wektorem T, można też określić
kolor oraz styl linii dla każdego z układów np. (sys1,'r',sys2,'y--',sys3,'gx').
[Y,T] = IMPULSE(SYS,T) wywołanie to zwraca wyniki symulacji w postaci wektora
odpowiedzi skokowej Y oraz wektora czasu trwania symulacji T.
[Y,T,X] = IMPULSE(SYS, ...) wywołanie to dla układu przestrzeni stanu zwraca
również wektor stanu X z kolejnych chwil symulacji.
LSIM symuluje odpowiedź czasową układu na zadany sygnał wejściowy.
LSIM(SYS,U,T) wywołanie to kreśli odpowiedź czasową układu na sygnał wejściowy
określony przez U i T przy zerowych warunkach początkowych, parametr U powinien
zawierać wierszami wektory sterowań dla kolejnych chwil czasu określonych w
wektorze T, stąd liczba wierszy macierzy U musi być równa liczbie elementów
wektora czasu T np. T = 0:0.01:5; U = sin(t); lsim(sys,u,t) symuluje odpowiedź
układu na sygnał wejściowy U przez 5 sekund.
LSIM(SYS,U,T,X0) wywołanie to dla układu opisanego równaniami stanu można
określić warunki początkowe X0.
LSIM(SYS1,SYS2,...,U,T,X0) wywołanie to symuluje odpowiedź kilku układów na
jednym wykresie przy czym warunki początkowe mogą ale nie muszą być określone,
każdemu układowi można przypisać odpowiedni kolor i rodzaj linii np.
lsim(sys1,'r',sys2,'y--',sys3,'gx',u,t).
[Y,T] = LSIM(SYS,U,...) wywołanie to zapamiętuje przebieg symulacji zapisując
przebieg wyjścia Y dla danych chwil czasu zapisanych w wektorze T, przy czym
macierz Y ma tyle samo wierszy co liczba elementów wektora T oraz tyle kolumn ile
jest wyjść układu.
[Y,T,X] = LSIM(SYS,U,...) wywołanie to dla układu opisanego równaniami stanu
zwraca też macierz stanów X, która ma tyle rzędów co liczba elementów wektora T i
tyle kolumn ile jest stanów układu.
Przykład: zasymulować i wykreślić odpowiedź danego układu SYS na zadany
sygnał prostokątny o okresie równym 4, czasie trwania 10 sekund i o czasie
próbkowania 0.1 sekundy:
5
1
2
+
+
−
=
s
s
s
SYS
W tym celu należy wykonać następujące wywołania:
>> [U,T]=GENSIG('SQUARE',4,10,0.1);
>> SYS=TF([1 -1],[1 1 5])
>> LSIM(SYS,U,T)
W wyniku ich realizacji otrzymamy odpowiedź układu:
Rys. 4.7. Odpowiedź układu na zadany sygnał wejściowy.
GENSIG tworzy sygnały wejściowe dla LSIM.
[U,T] = GENSIG(TYPE,TAU) wywołanie to tworzy skalarny sygnał U z pola TYPE o
okresie TAU. Generowane sygnały mogą być funkcjami typu:
TYPE = 'sin'
--- sinusoidalna
Time (sec.)
A
m
pl
itu
de
Linear Sim ulation Results
0
2
4
6
8
10
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
TYPE = 'square'
--- prostokątna
TYPE = 'pulse'
--- okresowo-impulsowa
GENSIG zwraca wektor T z czasem próbkowania i wektor U z wartościami sygnału
dla tych próbkowań, wszystkie generowane sygnały mają amplitudę równą jeden.
[U,T] = GENSIG(TYPE,TAU,TF,TS) wywołanie to dodatkowo można określić
długość trwania sygnału TF oraz czas próbkowania TS.
Przykład: utworzyć sygnał sinusoidalny o okresie równym 4, czasie trwania 10 sekund i o
czasie próbkowania 0.1 sekundy oraz przedstawić go na wykresie.
W tym celu należy wykonać następujące wywołania:
>> [U,T]=GENSIG('SIN',4,10,0.1);
>> plot(T,U)
>> AXIS([0 10 –1 2])
W wyniku ich realizacji otrzymamy następujący sygnał:
Rys. 4.8 Wykres zadanego sygnału sinusoidalnego.
2.3. Charakterystyki częstotliwościowe
Pełną informację o dynamice układu zawierają charakterystyki częstotliwościowe, które
są w zasadzie przeznaczone dla układów liniowych. Ilustrują one zależność stosunku
amplitud sygnału wyjściowego i wejściowego oraz przesunięcia fazowego tych sygnałów od
częstotliwości sygnału wejściowego harmonicznego. Przedstawione poniżej funkcje
programu MATLAB pozwalają wykreślić charakterystyki: Bodego, Nyquista, kartę Nicholsa,
obliczyć zapas modułu i fazy oraz wyznaczyć odpowiedź częstotliwościową dla zadanych
przez użytkownika częstotliwości.
BODE wyznacza odpowiedź częstotliwościową Bodego.
BODE(SYS) wywołanie to kreśli charakterystykę Bodego, częstotliwość i liczba
punktów są dobierane automatycznie.
0
1
2
3
4
5
6
7
8
9
10
-1
-0.5
0
0.5
1
1.5
2
BODE(SYS,{WMIN,WMAX}) wywołanie to kreśli charakterystykę Bodego dla
częstotliwości pomiędzy WMIN a WMAX w radianach na sekundę.
BODE(SYS,W) wywołanie to umożliwia samodzielne określenie częstotliwości w
radianach na sekundę dla których charakterystyka jest kreślona podanych w wektorze
W, do którego generowania można użyć funkcji LOGSPACE.
BODE(SYS1,SYS2,...,W) wywołanie to kreśli, z możliwością wykorzystania wektora
W, charakterystyki Bodego dla wielu układów na jednym wykresie, można też
określić kolor oraz styl linii dla każdego z układów np. bode(sys1,'r',sys2,'y-,sys3,'gx').
[MAG,PHASE,W] = BODE(SYS,...) wywołanie to zwraca wektor częstotliwości W
oraz odpowiadającym poszczególnym częstotliwościom macierze MAG i PHASE
zawierające odpowiednio wzmocnienie w dB oraz fazę w stopniach. Wywołanie to nie
kreśli wykresu. Jeżeli pominięty zostanie wektor W to funkcja określi go
samodzielnie.
Przykład:
W tym celu należy wykonać następujące wywołania:
>> SYS=TF([1 -1],[1 1 5])
>> BODE(SYS)
W wyniku ich realizacji otrzymamy odpowiedź układu:
Rys. 4.9 Charakterystyka częstotliwościowa Bodego.
NYQUIST wyznacza odpowiedź częstotliwościową Nyquista.
NYQUIST(SYS) wywołanie to kreśli charakterystykę Nyquista dla układu,
częstotliwość i liczba punktów dobierane są automatycznie.
NYQUIST(SYS,{WMIN,WMAX}) wywołanie to kreśli charakterystykę Nyquista dla
częstotliwości pomiędzy WMIN a WMAX w radianach na sekundę.
NYQUIST(SYS1,SYS2,...,W) wywołanie to kreśli, z możliwością wykorzystania
wektora W, charakterystyki Nyquista dla wielu układów na jednym wykresie, można
Frequency (rad/sec)
P
ha
se
(d
eg
);
M
ag
ni
tu
de
(d
B
)
Bode Diagrams
-20
-10
0
10
10
-1
10
0
10
1
-100
0
100
200
też określić kolor oraz styl linii dla każdego z układów np. bode(sys1,'r',sys2,'y--
',sys3,'gx').
[RE,IM,W] = NYQUIST(SYS,...) wywołanie to zwraca wektor częstotliwości W i
macierze RE i IM z rzeczywistą i urojoną częścią płaszczyzny zespolonej, wywołanie
to nie kreśli wykresu.
Przykład: dla danego układu SYS opisanego transmitancją należy wyznaczyć
charakterystykę częstotliwościową Nyquista:
2
2
15
9
2
3
2
+
−
−
+
+
=
s
s
s
s
s
SYS
W tym celu należy wykonać następujące wywołania:
>> SYS=TF([1 9 15],[1 -2 -1 2])
>> NYQUIST(SYS);
W wyniku ich realizacji otrzymamy odpowiedź układu:
Rys. 4.10 Charakterystyka częstotliwościowa Nyquista.
NICHOLS wyznacza odpowiedź częstotliwościową wykreśloną na karcie Nicholsa.
NICHOLS(SYS) wywołanie to kreśli charakterystykę Nicholsa dla układu,
częstotliwość i liczba punktów dobierane są automatycznie.
NICHOLS(SYS,{WMIN,WMAX}) wywołanie to kreśli charakterystykę Nicholsa dla
częstotliwości pomiędzy WMIN a WMAX w radianach na sekundę.
NICHOLS(SYS,W) wywołanie to umożliwia samodzielne określenie częstotliwości w
radianach na sekundę dla których charakterystyka jest kreślona podanych w wektorze
W, do którego generowania można użyć funkcji LOGSPACE.
NICHOLS(SYS1,SYS2,...,W) wywołanie to kreśli, z możliwością wykorzystania
wektora W, charakterystyki Nyquista dla wielu układów na jednym wykresie, można
też określić kolor oraz styl linii dla każdego z układów np.
nichols(sys1,'r',sys2,'y--',sys3,'gx').
[MAG,PHASE,W] = NICHOLS(SYS,...) wywołanie to zwraca wektor częstotliwości
W oraz odpowiadającym poszczególnym częstotliwościom macierze MAG i PHASE
Real Axis
Im
ag
in
ar
y
A
xi
s
Nyquist Diagrams
-1
0
1
2
3
4
5
6
7
8
-4
-3
-2
-1
0
1
2
3
4
zawierające odpowiednio wzmocnienie w dB oraz fazę w stopniach. Wywołanie to nie
kreśli wykresu. Jeżeli pominięty zostanie wektor W to funkcja określi go
samodzielnie.
NGRID wywołanie to kreśli na karcie Nicholsa linie stałych wartości modułu i
argumentu.
Przykład: dla danego układu SYS opisanego transmitancją należy wyznaczyć
charakterystykę częstotliwościową Nicholsa:
2
2
15
9
2
3
2
+
−
−
+
+
=
s
s
s
s
s
SYS
W tym celu należy wykonać następujące wywołania:
>> SYS=TF([1 9 15],[1 -2 -1 2])
>> NICHOLS(SYS);
>> NGRID
W wyniku ich realizacji otrzymamy odpowiedź układu:
Rys. 4.11 Karta Nicholsa.
MARGIN oblicza zapas modułu i fazy.
[GM,PM,WCG,WCP] = MARGIN(SYS) wywołanie to oblicza margines
wzmocnienia GM i fazy PM dla układu ciągłego SYS oraz odpowiadające im
częstotliwości graniczne odpowiednio WCG i WCP. Margines amplitudy jest
zdefiniowany jako
)
180
(
1
−
=
K
Gm
określający ile razy trzeba zwiększyć wzmocnienie by układ znalazł się na granicy
stabilności natomiast margines fazy jest podany w dB i liczony z równania
( )
Gm
10
log
20
.
[GM,PM,WCG,WCP] = MARGIN(MAG,PHASE,W) wywołanie to określa zapas
fazy i wzmocnienia na podstawie zadanych wektorów wzmocnienia MAG, fazy
Open-Loop Phase (deg)
O
pe
n-
Lo
op
G
ai
n
(d
B
)
Nichols Charts
50
100
150
200
250
300
350
-80
-60
-40
-20
0
20
40
6 dB
3 dB
1 dB
0.5 dB
0.25 d B
0 dB
-1 dB
-3 dB
-6 dB
-12 dB
-20 dB
-40 dB
-60 dB
-80 dB
PHASE oraz odpowiadającej im częstotliwości W pochodzących z charakterystyki
Bodego zarówno dla układów ciągłych jak i dyskretnych.
MARGIN(SYS) wywołanie to kreśli charakterystyki Bodego z zaznaczonymi
marginesami amplitudy i fazy.
Przykład: dla danego układu SYS opisanego transmitancją należy wyznaczyć zapas
amplitudy i fazy na charakterystyce częstotliwościowej Bodego
2
2
15
9
2
3
2
+
−
−
+
+
=
s
s
s
s
s
SYS
W tym celu należy wykonać następujące wywołania:
>> SYS=TF([1 9 15],[1 -2 -1 2])
>> MARGIN(SYS);
W wyniku ich realizacji otrzymamy odpowiedź układu:
Rys. 4.12 Zapas modułu i fazy.
FREQRESP oblicza odpowiedź częstotliwościową układu dla zadanych częstotliwości.
H = FREQRESP(SYS,W) wywołanie to oblicza odpowiedź częstotliwościową H
układu dla określonych w wektorze W pulsacji podanych jako liczby rzeczywiste w
radianach na sekundę.
Dla układu ciągłego opisanego transmitancją funkcja korzysta z przekształcenia s = j
ω
a dla
układu opisanego równaniami stanu korzysta z równania:
( )
(
)
B
A
I
j
D
j
H
1
−
−
+
=
ω
ω
(4.9.1)
natomiast dla układu dyskretnego przekształceniu ulega z = e
j
ω
Ts
gdzie T
s
jest czasem
próbkowania. Wyniki umieszczane są w trójwymiarowej tablicy H o wymiarach
liczba_wyjść
×
liczba_wejść
×
ilość_ elementów_ wektora_W.
Frequency (rad/sec)
P
ha
se
(d
eg
);
M
ag
ni
tu
de
(d
B
)
Bode Diagrams
-40
-20
0
20
Gm=11.5 dB (W cg=5.7); Pm=-58.2 deg. (W cp=2.6)
10
-1
10
0
10
1
10
2
0
100
200
300
2.4. Interfejs graficzny LTIView
LTIVIEW udostępnia interakcyjny graficzny interfejs użytkownika dla porównań czasowych
i częstotliwościowych odpowiedzi układów.
LTIVIEW(PLOTTYPE) wywołanie to generuje przypadkowy układ dla którego kreśli
odpowiedź. PLOTTYPE określa rodzaj odpowiedzi.
W pole PLOTTYPE można wstawić:
1) 'step’
2) 'impulse'
3) 'bode'
4) 'nyquist'
5) 'nichols'
6) 'sigma'
7) 'lsim'
8) 'initial'
LTIVIEW(PLOTTYPE,SYS) wywołanie to kreśli zadany rodzaj odpowiedzi dla
zadanego układu.
LTIVIEW(PLOTTYPE,SYS1,SYS2,...SYSN) wywołanie to kreśli odpowiedzi wielu
układów na jednym wykresie.
LTIVIEW(PLOTTYPE,SYS1,PLOTSTR1,SYS2,PLOTSTR2,...SYSN,PLOTSTRN)
wywołanie to kreśli odpowiedzi wielu układów na jednym wykresie ale umożliwia
wybór koloru i rodzaju linii dla każdego z tych układów poprzez wpisanie
odpowiednich znaków w pole PLOTSTR.
LTIVIEW(PLOTTYPE,SYS,OPTIONS) wywołanie to w polu OPTIONS pozwala
zadać użytkownikowi różne warunki dla których ma być nakreślona odpowiedź
układu np. jeśli odpowiedź ma być skokowa to w polu OPTIONS można umieścić
TFINAL określający czas trwania odpowiedzi.
3. Projektowanie regulatorów
3.1. Metoda miejsc geometrycznych pierwiastków
Projektując układ sterowania w dziedzinie częstotliwości koryguje się charakterystyki
częstotliwościowe układu tak, by po zamknięciu sprzężenia zwrotnego układ był stabilny, z
marginesami amplitudy i fazy o odpowiednich wartościach. W ten sposób pośrednio wpływa
się na położenie zer i biegunów układu zamkniętego. System MATLAB oferuje metodę
polegającą na projektowaniu układów sterowania poprzez bezpośredni wpływ na położenie
zer i biegunów układu zamkniętego, składającego się z obiektu i kompensatora połączonych
pętlą sprzężenia zwrotnego. Dla danego obiektu i postaci kompensatora bieguny i zera będą
położone w różnych miejscach płaszczyzny zespolonej w zależności od wzmocnienia
kompensatora. Poprzez odpowiedni dobór wzmocnienia można wpływać na własności
układu. Analizę położenia zer i biegunów układu wspomagają opisane niżej funkcje.
RLOCUS kreśli położenie pierwiastków równania.
RLOCUS(SYS) wywołanie to kreśli dla układów SISO opisanych transmitancją gdzie
N(s) jest licznikiem a D(s) mianownikiem położenie pierwiastków równania :
H(s) = D(s) + KN(s) = 0
(4.10.1)
dla dodatnich wzmocnień K, które zostaną utworzone w celu nakreślenia wykresu.
RLOCUS(SYS,K) wywołanie to używa wektora wzmocnień K określonego przez
użytkownika.
[R,K] = RLOCUS(SYS) lub R = RLOCUS(SYS,K) wywołania te pozwalają
obliczyć położenie pierwiastków równania i zapamiętać je wraz z odpowiadającymi
im wzmocnieniami. K jest wektorem kolumnowym zawierającym pulsacje, dla
których określano położenia biegunów. Bieguny układu zapamiętywane są w
macierzy R o liczbie wierszy równej liczbie elementów wektora K i liczbie kolumn
równej stopniowi mianownika transmitancji. W każdym wierszu macierzy R
pamiętane są bieguny odpowiadające danej wartości wzmocnienia.
Przykład: dla danego układu SYS opisanego transmitancją należy wyznaczyć zera i
bieguny układu.
2
2
15
9
2
3
2
+
−
−
+
+
=
s
s
s
s
s
SYS
W tym celu należy wykonać następujące wywołania:
>> SYS=TF([1 9 15],[1 -2 -1 2])
>> RLOCUS(SYS);
W wyniku ich realizacji otrzymamy rozmieszczenie zer ,,o’’ i biegunów ,,x’’ na
płaszczyźnie zespolonej :
-20
-15
-10
-5
0
5
-10
-8
-6
-4
-2
0
2
4
6
8
10
Real A xis
Im
ag
A
xi
s
Rys. 4.13 Położenie zer i biegunów układu na płaszczyźnie zespolonej.
RLOCFIND znajduje położenie i wzmocnienie wybranych przez użytkownika biegunów.
[K,POLES] = RLOCFIND(SYS,P) wywołanie to umożliwia wybranie myszką
żądanego bieguna układu SISO opisanego transmitancją lub równaniami stanu jeśli w
aktywnym oknie znajduje się wykres zależności położenia biegunów od wzmocnienia
uzyskany funkcją RLOCUS. Możliwe jest wybranie myszką zadanej przez parametr P
liczby biegunów. Pobrane wartości są zapamiętywane w wektorze K i macierzy
POLES analogicznej do macierzy R przedstawionej w opisie funkcji RLOCUS.
Punkty należy wybierać na lub w pobliżu trajektorii biegunów.
ACKER projektuje rozmieszczenie biegunów dla układu SISO.
K = ACKER(A,B,P) wywołanie to oblicza używając formuły Ackermana macierz
wzmocnienia sprzężenia zwrotnego K dla układu z jednym wejściem opisanym
równaniem
Bu
Ax
x
+
=
ˆ
ze sprzężeniem zwrotnym u = -Kx i wartościami określonymi
w wektorze P = eig(A-B*K).
PLACE projektuje rozmieszczenie biegunów dla układu MIMO.
K = PLACE(A,B,P) wywołanie to oblicza macierz sprzężenia zwrotnego K dla
wartości własnych określonych w wektorze P =EIG(A-B*K).
ESTIM oszacowuje stan dla przyjętego wzmocnienia.
EST = ESTIM(SYS,L) wywołanie to tworzy oszacowanie stanu i wyjścia układu dla
układu opisanego równaniem stanu z oszacowanym wzmocnieniem L. Wszystkie
wejścia w układu są procesami stochastycznymi a wszystkie wyjścia są pomierzone.
Dla układu ciągłego opisanego równaniami stanu:
ω
B
Ax
x
+
=
&
(4.10.2)
ω
D
Cx
y
+
=
(4.10.3)
EST generuje oszacowanie stanu i wyjścia z równań :
(
)
x
C
y
L
x
A
x
ˆ
ˆ
ˆ
−
+
=
&
(4.10.4) oraz
x
I
C
x
y
ˆ
ˆ
ˆ
=
.
(4.10.5)
SGRID, ZGRID wykreśla linie stałego współczynnika wzmocnienia oraz stałej
częstotliwości własnej dla układu ciągłego i dyskretnego.
SGRID wywołanie to kreśli linie stałego współczynnika wzmocnienia oraz stałej
częstotliwości własnej dla układu ciągłego
SGRID(‘NEW’) wywołanie to powoduje oczyszczenie okna graficznego przed
nakreśleniem linii
SGRID(Z,WN) wywołanie to pozwala na podanie przez użytkownika współczynnika
tłumienia i częstotliwości, dla których linie te mają być nakreślone
Funkcja ZGRID ma analogiczne wywołania.
3.2. Regulatory optymalne
Układy optymalne zapewniają uzyskanie ekstremalnej wartości zadanego z góry
kryterium optymalności przy zachowaniu ograniczeń na sterowanie i składowe wektora stanu.
W celu wyznaczenia tego rozwiązania wykorzystuje się zagadnienia algebry liniowej
przedstawione w podrozdziale 4.12.
LQR projektuje liniowo-kwadratowy regulator dla układów ciągłych.
[K,S,E] = LQR(A,B,Q,R,N) wywołanie to oblicza macierz sprzężeń zwrotnych K,
macierz S będącą rozwiązaniem równania Ricattiego:
(
)
(
)
0
1
=
+
+
+
−
+
−
Q
N
S
B
R
N
SB
SA
S
A
T
T
T
(4.11.1)
w którym macierze A, B są odpowiednio macierzą stanu i sterowania, a Q i R są
macierzami wag w optymalizowanym kryterium jakości o postaci:
( )
(
)
∫
∞
+
+
=
0
2
dt
Nu
x
Ru
u
Qx
x
u
J
T
T
T
(4.11.2)
dla układu ciągłego u = -Kx opisanego równaniem stanu
Bu
Ax
x
+
=
&
.
Pominięcie N zakłada wykonanie obliczeń dla N = 0.
Zwracany jest również wektor E zawierający wartości własne zamkniętej pętli sprzężenia
zwrotnego obliczone z wywołania : E = EIG(A-B*K).
DLQR projektuje liniowo-kwadratowy regulator dla układów dyskretnych.
[K,S,E] = DLQR(A,B,Q,R,N) wywołanie to oblicza macierz sprzężeń zwrotnych K i
macierz S będącą rozwiązaniem równania Ricattiego:
(
)(
) (
)
0
1
=
+
+
+
+
−
−
−
Q
N
SA
B
R
XB
B
N
SB
A
S
SA
A
T
T
T
T
T
(4.11.3)
w którym macierze A, B są odpowiednio macierzą stanu i sterowania, a Q i R są
macierzami wag w optymalizowanym kryterium jakości o postaci
( )
[ ] [ ] [ ] [ ]
[ ] [ ]
(
)
∑
∞
=
+
+
=
1
2
n
T
T
T
n
Nu
n
x
n
Ru
n
u
n
Qx
n
x
u
J
(4.11.4)
dla dyskretnego układu u[n] = -Kx[n] opisanego równaniami stanu:
[ ]
[ ]
[ ]
n
Bu
n
Ax
n
x
+
=
+
1
.
(4.11.5)
Pominięcie N zakłada wykonanie obliczeń dla N = 0.
Zwracany jest również wektor E zawierający wartości własne zamkniętej pętli sprzężenia
zwrotnego obliczone z wywołania : E = EIG(A-B*K).
LQRY projektuje liniowo-kwadratowy regulator z wagowym wyjściem.
[K,S,E] = LQRY(SYS,Q,R,N) oblicza optymalną macierz sprzężenia zwrotnego K dla
układu :
- ciągłego u = -Kx opisanego równaniami stanu:
Bu
Ax
x
+
=
&
,
Du
Cx
y
+
=
(4.11.6)
dla wskaźnika jakości o postaci:
( )
(
)
∫
∞
+
+
=
0
2
dt
Nu
y
Ru
u
Qy
y
u
J
T
T
T
(4.11.7)
- dyskretnego u[n] = -Kx[n] opisanego równaniami stanu:
x[n+1] = Ax[n] + Bu[n], y[n] = Cx[n] + Du[n]
(4.11.8)
dla wskaźnika jakości o postaci:
( )
[ ] [ ] [ ] [ ]
[ ] [ ]
(
)
∑
∞
=
+
+
=
1
2
n
T
T
T
n
Nu
n
y
n
Ru
n
u
n
Qy
n
y
u
J
(4.11.9)
Pominięcie N zakłada wykonanie obliczeń dla N = 0. Zwracany jest również wektor S
z rozwiązaniem algebraicznego równania Riccatiego oraz wektor E zawierający
wartości własne zamkniętej pętli sprzężenia zwrotnego obliczone z wywołania : E =
EIG(A-B*K).
LQRD projektuje dyskretny liniowo-kwadratowy regulator dla dyskretnego wskaźnika
jakości.
[K,S,E] = LQRD(A,B,Q,R,TS) wywołanie to oblicza optymalną macierz wzmocnienia
K dla dyskretnego układu u[n] = -Kdx[n] opisanego równaniami stanu:
x[n+1] = Ad x[n] + Bd u[n],
(4.11.10)
gdzie Ad i Bd wyznacza się z wywołania [Ad,Bd] = C2D(A,B,TS), dla wskaźnika
jakości o postaci :
( )
(
)
∫
∞
+
=
0
dt
Ru
u
Qx
x
u
J
T
T
.
(4.11.11)
Zwraca również rozwiązanie dyskretnego równania Riccatiego S oraz wektor wartości
własnych zamkniętej pętli z wywołania E = EIG(Ad-Bd*K).
KALMAN tworzy ciągły lub dyskretny estymator Kalmana.
[KEST,L,P] = KALMAN(SYS,QN,RN,NN) wywołanie to projektuje estymator
Kalmana KEST dla ciągłego obiektu
SYS.
Mając dany układ ciągły opisany przez :
Gw
Bu
Ax
x
+
+
=
&
równanie stanu
y = Cx + Du + Hw + v
równanie pomiaru
ze znanymi : wejściami u, przebiegiem białego szumu w, zmierzonym szumem v, i
kowariancjami szumu E(ww
T
) = QN, E(vv
T
) = RN, E(wv
T
) = NN można oszacować
stan xˆ dla minimalnej kowariancji błędu P określonej wzorem:
)
}
ˆ
}{
ˆ
({
lim
T
t
x
x
x
x
E
P
−
−
=
∞
→
.
(4.11.12)
Optymalnym rozwiązaniem jest filtr Kalmana opisany równaniami :
)
ˆ
(
ˆ
ˆ
Du
x
C
y
L
Bu
x
A
x
u
−
−
+
+
=
&
(4.11.13)
u
D
x
I
C
x
y
+
=
0
ˆ
ˆ
ˆ
(4.11.14)
gdzie wzmocnienie filtru L jest zdeterminowane przez rozwiązanie algebraicznego
równania Riccatiego.
Ten estymator używa znanych wejść u i zmierzonych y
u
do generowania wyjścia i
oszacowania stanu xˆ oraz yˆ .
Rys. 4.14 Schemat blokowy estymatora Kalmana.
[KEST,L,P,M,Z] = KALMAN(SYS,QN,RN,NN) wywołanie to projektuje estymator
Kalmana KEST dla dyskretnego plant SYS.
Mając dany układ dyskretny opisany równaniami :
x[n+1] = Ax[n] + Bu[n] + Gw[n]
(4.11.15)
y
v
[n] = Cx[n] + Du[n] + Hw[n] + v[n]
(4.11.16)
oraz mając dane kowariancje szumu :
E(w[n]w[n]
T
) = QN,
(4.11.17)
E(v[n]v[n]
T
) = RN,
(4.11.18)
E(w[n]v[n]
T
) = NN,
(4.11.19)
estymator Kalmana można opisać równaniami :
[
]
[
]
[ ]
[ ]
[
]
[ ]
(
)
n
Du
n
n
x
C
n
y
L
n
Bu
n
n
x
A
n
n
x
v
−
−
−
+
+
−
=
+
1
|
ˆ
1
|
ˆ
|
1
ˆ
(4.11.20)
[ ]
[ ]
(
)
[
]
(
)
[ ]
[ ]
−
−
+
−
−
−
=
n
y
n
u
M
MD
CM
D
CM
I
n
n
x
MC
I
MC
I
C
n
n
x
n
n
y
v
1
|
ˆ
|
ˆ
|
ˆ
(4.11.21)
Macierze wzmocnienia L i M są wyprowadzone z rozwiązania dyskretnego równania
Riccatiego. Wzmocnienie M jest użyte do uaktualnienia predykcji
[
]
1
|
ˆ
−
n
n
x
użytej
do nowego pomiaru y
v
[n] :
[ ] [
]
[ ]
[
]
[ ]
(
)
n
Du
n
n
x
C
n
y
M
n
n
x
n
n
x
v
−
−
−
+
−
=
1
|
ˆ
1
|
ˆ
|
ˆ
(4.11.22)
KALMD tworzy dyskretny estymator Kalmana dla ciągłego obiektu.
[KEST,L,P,M,Z] = KALMD(SYS,QN,RN,TN) wywołanie to wyznacza dyskretny
estymator Kalmana KEST o czasem próbkowania Ts dla ciągłego układu opisanego
przez:
Gw
Bu
Ax
x
+
+
=
&
równanie stanu
(4.11.23)
y = Cx + Du + v
równanie pomiaru
(4.11.24)
ze znanym przebiegiem i pomiarem szumów oraz ich kowariancjami :
E(w) = E(v) = 0,
(4.11.25)
E(ww
T
) = QN,
(4.11.26)
E(vv
T
) = RN,
(4.11.27)
E(wv
T
) = 0.
(4.11.28)
Wywołanie to zwraca estymator Kalmana KEST, macierze wzmocnień L i M oraz
błędy kowariancji P i Z.
LQGREG tworzy liniowo-kwadratowy Gaussowski regulator.
RLQG = LQGREG(KEST,K) wywołanie to tworzy regulator LQG poprzez połączenie
estymatora Kalmana KEST zaprojektowanego przez polecenie KALMAN ze
wzmocnieniem sprzężenia zwrotnego K zaprojektowanego przez LQR, DLQR lub
LQRY :
Rys. 4.15 Schemat blokowy regulatora LQR.
Regulator RLQG jest określony przez poniższe równania :
(
)
[
]
v
Ly
x
K
LD
B
LC
A
x
+
−
−
−
=
ˆ
ˆ&
(4.11.29)
oraz
x
K
u
ˆ
−
=
gdzie y
v
jest wektorem z mierzonymi wyjściami układu.
Zagadnienia algebry liniowej optymalnych układów sterowania.
LYAP rozwiązuje równanie Lapunowa dla układów ciągłych.
X = LYAP(A,Q) wywołanie to rozwiązuje równanie macierzowe Lapunowa
przedstawione poniżej:
0
=
+
+
Q
XA
AX
T
(4.12.1)
X = LYAP(A,B,C) wywołanie to rozwiązuje ogólne równanie macierzowe Lapunowa,
inaczej zwane równaniem Sylvestra przedstawione poniżej:
C
XB
AX
−
=
+
(4.12.2)
DLYAP rozwiązuje dyskretne równanie Lapunowa.
X = DLYAP(A,Q) wywołanie to rozwiązuje dyskretne równanie macierzowe
Lapunowa przedstawione poniżej:
Q
X
XA
A
T
=
−
(4.12.3)
CARE rozwiązuje algebraiczne równanie Riccatiego dla układów ciągłych
[X,L,G,RR] = CARE(A,B,Q) wywołanie to oblicza rozwiązanie X ciągłego
algebraicznego równania Riccatiego:
( )
0
=
+
−
+
=
Q
XBB
XA
X
A
X
Ric
T
T
.
(4.12.4)
Macierz X jest macierzą symetryczną związaną ze stabilnym rozwiązaniem równania
Eic(X) = 0. Zwracane są też:
- wartości własne L z A-BB
T
X ,
- macierz wzmocnienia G = B
T
X,
wartości względne residuów definiowane jako
( )
F
F
X
X
Ric
RR
=
(4.12.5)
[X,L,G,RR] = CARE(A,B,Q,R,S,E) wywołanie to oblicza rozwiązanie X ciągłego
algebraicznego równania Riccatiego:
( )
(
) (
)
0
1
=
+
+
+
−
+
=
−
Q
S
XE
B
R
S
XB
E
XA
E
XE
A
X
Ric
T
T
T
T
T
. (4.12.6)
Jeżeli pominięte zostaną R, S i E to obliczenia zostaną wykonane dla R=I, S=0, E=I.
Macierz wzmocnienia G obliczona jest ze wzoru:
(
)
T
T
S
XE
B
R
G
+
=
−
1
(4.12.7)
a wartości własne sprzężenia zwrotnego L=EIG(A-B*G,E).
DARE rozwiązuje dyskretne równanie Riccatiego.
[X,L,G,RR] = DARE(A,B,Q,R) wywołanie to oblicza rozwiązanie X dyskretnego
równania Riccatiego:
( )
(
)
0
1
=
+
+
−
−
=
−
Q
XA
B
R
XB
B
XB
A
X
XA
A
X
Ric
T
T
T
T
. (4.12.8)
Macierz X jest macierzą symetryczną związaną ze stabilnym rozwiązaniem równania
Ric(X) = 0. Zwracane są też:
- wartości własne L z
(
)
XA
B
R
XB
B
B
A
A
T
T
cl
1
−
+
−
=
,
- macierz wzmocnienia G z
(
)
XA
B
R
XB
B
G
T
T
1
−
+
=
,
- wartości względne residuów definiowane jako
( )
F
F
X
X
Ric
RR
=
[X,L,G,RR] = DARE(A,B,Q,R,S,E) wywołanie to oblicza jedyne symetryczne stabilne
rozwiązanie X dyskretnego równania Riccatiego:
0
)
(
)
)(
(
)
(
1
=
+
+
+
+
−
−
=
−
Q
S
XA
B
R
XB
B
S
XB
A
XE
E
XA
A
X
Ric
T
T
T
T
T
T
(4.12.9)
pominięcie R, S i E powoduje wykonanie obliczeń dla R = I, S = 0 i E = I. Zwraca
również macierz wzmocnienia jest obliczana ze wzoru:
)
(
)
(
1
T
T
T
S
XA
B
R
XB
B
G
+
+
=
−
,
(4.12.10)
wektor L z wartościami własnymi pętli sprzężenia z wywołania L=EIG(A-B*G,E)
oraz normę Frobeniusa z macierzy residuów.
3.3. Przykłady