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