WOJSKOWA AKADEMIA TECHNICZNA
im. Jarosława Dąbrowskiego
ZAKA
AD AWIONIKI I UZBROJENIA LOTNICZEGO
Przedmiot:
PODSTAWY AUTOMATYKI I AUTOMATYZACJI
ĆWICZENIE LABORATORYJNE Nr 2
MODELOWANIE UKA
ADU REGULACJI AUTOMATYCZNEJ Z
WYKORZYSTANIEM PAKIETU MATLAB/SIMULINK
Warszawa 2014
1
1.1. CEL ĆWICZENIA
Celem ćwiczenia jest zapoznanie studentów z zapisem modeli
dynamicznych projektowanych układów automatycznego sterowania w
środowisku MATLAB-SIMULINK oraz symulacji i analizy własności tych
układów.
1.2. PRZEDMIOT ĆWICZENIA
Przedmiotem ćwiczenia jest szczegółowa analiza układów
automatycznego sterowania z wykorzystaniem nowoczesnych narzędzi
modelowania i symulacji komputerowej. Studenci zapoznają się ze sposobem
zapisu i prezentacji układów automatyki, określą parametry statyczne
i dynamiczne układu na podstawie otrzymanych charakterystyk czasowych
i częstotliwościowych.
1.3. WIADOMOŚCI OGÓLNE
1.3.1. Wprowadzenie
MATLAB jest interakcyjnym środowiskiem do wykonywania naukowych
i inżynierskich obliczeń oraz do wizualizacji danych. Zakres zastosowań
obejmuje różne dziedziny nauki i techniki, w tym biologię, medycynę,
ekonomię, meteorologię i wiele innych.
Zasadnicze jego zalety to: możliwość szybkiego uzyskania rezultatów
złożonych obliczeń i przedstawienia ich w postaci wykresów dwu- lub
trójwymiarowych, a także jako mapy wielobarwnej.
MATLAB jest w zasadzie językiem programowania wysokiego poziomu.
Podstawowym typem danych, na których operuje jest macierz rzeczywista lub
zespolona. W programie tym nie deklaruje się zmiennych, są one
przechowywane w przestrzeni roboczej i są dostępne poprzez nazwę.
Najważniejszymi zaletami pakietu MATLAB są:
�� Otwarta architektura pakietu, na którą składają się:
� M-pliki, umożliwiające definiowanie poleceń i algorytmów
obliczeniowych;
� MEX-pliki i pliki ASCII służące do wymiany danych i wyników
obliczeń pomiędzy MATLAB em, a innymi programami;
� Grafika służy do wizualizacji danych i wyników obliczeń,
możliwa jest prosta animacja i efekty dz
więkowe;
� GUI (ang. Graphics User Interface)- interfejs graficzny dający
możliwość pracy interaktywnej za pomocą okienek
edycyjnych, przycisków, suwaków i menu.
� Usługi DDE (ang. Dynamic Data Exchange) realizują
statyczną lub dynamiczną wymianę danych tekstowych i
graficznych pomiędzy programami w środowisku MS
Windows. Za pomocą tej usługi można przesyłać obiekty
graficzne oraz teksty i dane liczbowe pomiędzy MATLAB em,
a innymi programami Windows.
�� Toolbox- jest to ponad 20 wyspecjalizowanych pakietów
oprogramowania, m.in. SIMULINK;
SIMULINK jest to interaktywny pakiet do modelowania i symulacji
układów dynamicznych. Umożliwia on tworzenie wielopoziomowych
schematów blokowych. Obiekty biblioteczne są umieszczone w okienkach jako
ikony. Program ten może współpracować z:
�� Rael Time Toolbox (RTW)- pakietem czasu rzeczywistego z
generatorem kodu C
i asemblera;
�� Z rzeczywistym obiektem regulacji poprzez układy sprzęgające;
�� SIMULINK accelerator- automatycznie kompiluje i linkuje modele
SIMULINK a, dając znaczną poprawę szybkości symulowania
modelowanego systemu.
1.3.2. Modele matematyczne układów dynamicznych
Układem nazywamy umownie wyodrębniony ze środowiska układ
fizyczny lub jego część. Wielkości charakteryzujące oddziaływanie środowiska
na układ nazywać będziemy wymuszeniami lub wielkościami wejściowymi.
Wymuszenia dzielimy na wielkości sterujące i wielkości zakłócające.
Wielkościami sterującymi u(t) nazywamy wymuszenia zmieniane celowo, a
wielkościami zakłócającymi z(t) wielkości podlegające zmianom
przypadkowym (losowym). Wielkości charakteryzujące oddziaływanie układu
na środowisko nazywać będziemy odpowiedziami lub wielkościami
wyjściowymi y(t).
wejście wyjście
u1(�t)� y1(�t)�
x1(�t)�
�� ł�
u2(�t)�
y2(�t)�
ę�x
(�t)�ś�
2
ę� ś�
x(�t)�=�
ę� ś�
M�
M� M�
un(�t)� yn(�t)�
ę�
(�t)�ś�
��xn ��
Rys.1. Schemat układu dynamicznego.
Przebieg procesu dynamicznego w czasie zależy nie tylko od wartości
wymuszeń w danej chwili, ale również od wartości tych wymuszeń w
przeszłości. Można więc powiedzieć, że proces dynamiczny ma pamięć, w
której są gromadzone skutki przeszłych oddziaływań. Dla prezentacji tej
pamięci wprowadza się pojęcie stanu procesu.
Stanem procesu nazywamy najmniejszy liczebnie zbiór wielkości x1(t),
x2(t), ..., xn(t) określających w pełni skutki przeszłych oddziaływań na układ,
który jest wystarczający do przewidywania przebiegu procesu w przeszłości.
Wielkości te nazywamy zmiennymi stanu. Znajomość stanu procesu w chwili t0
oraz wymuszeń w przedziale [t0,t1) wystarcza do określenia przebiegów
odpowiedzi i stanu procesu w przedziale [t0,t1).
1.3.2.1 Opis układów dynamicznych w przestrzeni stanu
Układy, których miarą zmiany procesu w czasie jest pochodna wektora
stanu x(t), a stan procesu dla t>t0 zależy tylko od stanu w chwili t0 oraz od
sygnału sterującego u(t) w przedziale [t0,t] można opisać równaniem:
&�
x(�t)� =� f [�x(�t)�,u(�t)�,t]� (1)
3
o danym warunku początkowym:
x(�t0)�=� x0 (2)
Dla pełnego opisu układu potrzebne jest równanie wiążące wektor
odpowiedzi y(t) z wektorem stanu x(t) i wektorem wymuszenia u(t) postaci:
y(�t)� =� g[�x(�t)�,u(�t)�,t]� (3)
Dla układów liniowych niestacjonarnych (o parametrach zależnych od
czasu) równania (1) i (2) przyjmują postać:
&�
x(�t)�=� A(�t)�x(�t)�+� B(�t)�u(�t)� (4)
y(�t)� =� C(�t)�x(�t)�+� D(�t)�u(�t)� (5)
gdzie: A(t)-macierz układu o wymiarze n x n, B(t)- macierz sterowania o wymiarze n x p,
C(t) - macierz odpowiedzi o wymiarze q x n, D(t)- macierz transmisyjna o wymiarze q x p.
Schemat blokowy takiego układu przedstawia rys.2.
D(t)
.
+
u(t) x(t) x(t)
+ +
y(t)
B(t) C(t)
��
+
A(t)
Rys.2. Schemat blokowy ciągłego układu liniowego niestacjonarnego.
Dla układów liniowych stacjonarnych (o parametrach niezależnych od
czasu) elementy macierzy A, B, C, D są stałe i nie zależą od czasu. W tym
przypadku równania (4) i (5) przyjmują postać:
&�
x(�t)� =� Ax(�t)�+� Bu(�t)� (6)
y(�t)�=� Cx(�t)�+� Du(�t)� (7)
1.3.2. Opis układow dynamicznych za pomocą transmitancji
operatorowe i widmowej
Transmitancją operatorową G(s) jednowymiarowego układu liniowego
stacjonarnego nazywamy stosunek transformaty odpowiedzi Y(s) do
transformaty wymuszenia U(s) tego układu przy zerowych warunkach
początkowych.
4
Y(�s)�
G(�s)�=� (8)
X(�s)�
Dla układów liniowych stacjonarnych opisanych równaniami
różniczkowymi liniowymi o stałych współczynnikach transmitancja
operatorowa jest ilorazem dwóch wielomianów o postaci:
L(�s)�
G(�s)� =� (9)
M(�s)�
gdzie:
m
L(�s)�=� bisi (10)
��
i=�0
n
j
M(�s)� =� s , n ł� m (11)
��a
j
j=�0
Pierwiastki równania L(s)=0 nazywamy zerami, a pierwiastki równania
M(s)=0 nazywamy biegunami transmitancji operatorowej opisanej
równaniem (9).
Transmitancją widmową G(jw�) liniowego układu stacjonarnego nazywać
będziemy wielkość określoną jako stosunek wartości zespolonej składowej
wymuszonej odpowiedzi Y(jw�) wywołanej wymuszeniem sinusoidalnym do
wartości zespolonej tego wymuszenia U(jw�):
Y(� jw�)�
G(� jw�)�=� (12)
U(� jw�)�
Transmitancja widmowa jest wielkością zespoloną zależną od
parametrów układu i pulsacji wymuszenia w�:
G(�jw�)�=� P(�w�)�+� jQ(�w�)� (13)
gdzie:
P(�w�)�=� Re[�G(�jw�)�]� (14)
Q(�w�)�=� Im[�G(�jw�)�]� (15)
Transmitancja widmowa G(jw�) z transmitancją operatorową G(s)
związana jest zależnością:
G(�jw�)�=� G(�s)� (16)
s=� jw�
Zależność tą otrzymamy z definicji transmitancji operatorowej (8)
i transmitancji widmowej (12) po uwzględnieniu, że przekształcenie Fouriera
5
jest szczególnym przypadkiem dla s=jw� przekształcenia Laplace a oraz, że
transformata Fouriera przebiegu sinusoidalnie zmiennego jest równa wartości
zespolonej tego przebiegu.
1.3.3. Kryteria jakości regulacji
Do najczęściej stosowanych kryteriów jakości regulacji należą:
�� Kryteria zapasu stabilności;
�� Kryteria rozkładu pierwiastków równania charakterystycznego;
�� Kryteria czasowe;
�� Kryteria częstotliwościowe;
�� Kryteria całkowe.
Miarą zapasu stabilności jest odległość charakterystyki amplitudowo-
fazowej układu otwartego od punktu (-1,j0) (rys.3.a). Zapas stabilności
amplitudy (modułu) jest określony następującymi wielkościami:
D�L1 =� 20log(�k1)�, D�L2 =� 20log(�k2)� (17)
gdzie k1 i k2 są określone w sposób podany na rys.3.
Zapas amplitudy jest tym większy, im większe są D�L1 i D�L2. W dobrze
tłumiących układach wielkości te wynoszą od 6 do 12 dB.
Rys.3. Określenia zapasu stabilności i fazy na podstawie charakterystyki
amplitudowo-fazowej układu otwartego.
W przypadku szczególnym, gdy charakterystyka amplitudowo-fazowa
układu otwartego ma punkty przecięcia z osią liczb rzeczywistych tylko na
prawo od punktu (-1,j0), zapas stabilności amplitudy określamy za pomocą
jednej wielkości (rys.3.b):
D�L =� 20log k (18)
Zapasem stabilności fazy D�j� nazywamy kąt zawarty między osią liczb
rzeczywistych, a półprostą łączącą początek układu współrzędnych 0 z
punktem B przecięcia okręgu o promieniu równym 1 i środku w początku
układu współrzędnych z charakterystyką amplitudowo-fazową układu
otwartego (rys.3.), czyli:
6
 D�j� =�180o� +�j� (19)
gdzie: j� jest argumentem transmitancji widmowej układu otwartego K(jw�), odpowiadającym
modułowi równemu 1.
Zapas stabilności fazy jest tym większy, im większy jest kąt D�j�, który
zwykle przyjmuje się w granicach od 30��do 60��. Zadając określoną wartość
zapasu stabilności amplitudy i fazy określamy obszar, przez który nie może
przechodzić charakterystyka amplitudowo-fazowa układu otwartego.
Rys.4. Określenia zapasu stabilności amplitudy i fazy na podstawie logarytmicznych
charakterystyk amplitudowej i fazowej układu otwartego
Dla logarytmicznych charakterystyk amplitudowych i fazowych układu
otwartego zapas stabilności zaznaczamy jak na rys.4, a definicje tych
wielkości są następujące:
�� Zapasem stabilności amplitudy nazywamy ujemne odchylenie
logarytmicznej charakterystyki amplitudowej układu otwartego od
wartości 0 dB dla pulsacji odcięcia fazy w�0, przy której
logarytmiczna charakterystyka fazowa przecina prostą -180��.
�� Zapasem stabilności fazy nazywamy dodatnie odchylenie
logarytmicznej charakterystyki fazowej układu otwartego od
wartości -180�� dla pulsacji odcięcia modułu w�m, przy której
logarytmiczna charakterystyka amplitudowa przecina oś odciętych
(wartość 0 dB).
�� Wymagania stawiane układom dotyczące zapasu stabilności i
szybkości działania można sformułować w postaci warunków, aby
pierwiastki równania charakterystycznego układu znajdowały się w
określonej części lewej półpłaszczyzny zmiennej s. W ten sposób
budujemy kryterium położenia biegunów i zer układu,
z którym związane są pojęcia stopnia stabilności h� oraz
oscylacyjności m�.
7
Stopniem stabilności h� nazywamy wartość bezwzględną części
rzeczywistej pierwiastka położonego w lewej półpłaszczyz
nie najbliżej osi
urojonej (rys.5). Natomiast oscylacyjnością m� nazywamy największą wartość
bezwzględną stosunku części urojonej b�k do części rzeczywistej a�k pierwiastka
zespolonego sk=-a�ką�jb�k spośród wszystkich pierwiastków równania
charakterystycznego układu:
Imsk b�k
m� =� max =� max (20)
sk
Re sk sk a�k
Rys.5. Określenia stopnia stabilności.
Założenie dla układu określonej oscylacyjności m� jest równoważne
wymaganiu, aby pierwiastki równania charakterystycznego układu leżały w
części lewej półpłaszczyzny zmiennej zespolonej s ograniczonej z prawej
strony dwiema półprostymi nachylonymi względem osi rzeczywistej pod kątem
nie mniejszym niż j�=arctg(m�).
Kryteria czasowe określają pewne parametry lub właściwości przebiegów
wielkości regulowanych w stanie nieustalonym wywołanych skokową zmianą
wielkości zadających lub zakłóceń. Skłonność układu do oscylacji można w
przybliżeniu ocenić za pomocą wielkości zwanej przeregulowaniem, którą
definiujemy jako maksymalną różnicę między wartością chwilową wielkości
regulowanej w stanie nieustalonym, a jej wartością ustaloną, odniesioną do tej
wartości ustalonej:
ymax -� yu
d� =� , yu ą� 0 (21)
yu
8
gdzie: ymax jest maksymalną wartością chwilową wielkości regulowanej, a yu jest wartością
ustaloną tej wielkości. (rys.6.). Przyjmuje się, że dopuszczalne jest przeregulowanie w
granicach od 0,1 do 0,3.
Kolejnym parametrem określającym szybkość działania układu jest czas
regulacji tr zdefiniowany jako czas (mierzony od chwilowej skokowej zmiany
wymuszenia), po upływie którego wielkość regulowana różni się od wartości
ustalonej trwale mniej niż założona wartość D�, przyjmowana zwykle od 2 do
5% wartości ustalonej.
Między czasem regulacji tr, a stopniem stabilności h� zachodzi zależność:
1 1
tr =� ln (22)
h� D�
Z kryterium częstotliwościowym oceny jakości układu związane są
następujące pojęcia: szczyt rezonansowy Mp, pulsacja rezonansowa w�r oraz
polsacja odcięcia w�o.
Szczytem rezonansowym Mp nazywamy największą wartość
charakterystyki amplitudowej układu zamkniętego. Im większy jest zapas
stabilności i im słabsze jest tłumienie oscylacji, tym większy jest szczyt
rezonansowy Mp. Dla układów dobrze tłumiących oscylacje wartość tego
parametru powinna być zawarta w granicach od 1,1 do 1,5.
Pulsacją rezonansową w�r nazywamy pulsację odpowiadającą
największej wartości charakterystyki amplitudowej układu zamkniętego.
Korzystając z charakterystyki amplitudowo-fazowej układu otwartego możemy
w przybliżeniu określić pulsację rezonansową w�r jako pulsację odpowiadającą
punktowi tej charakterystyki położonemu najbliżej punktu (-1,j0). Pulsacja
rezonansowa niewiele różni się od pulsacji drgań występujących w układzie w
stanie nieustalonym. Czas regulacji tr z pulsacją rezonansową związany jest
zależnością:
3p�
tr � (23)
w�r
Pulsacją odcięcia w�o nazywać będziemy pulsację odpowiadającą
punktowi przecięcia charakterystyki amplitudowo-fazowej układu otwartego
K(jw�) z okręgiem o promieniu jednostkowym i środku w początku układu
współrzędnych, |K(jw�)|=1. Pulsację odcięcia można określić również jako
pulsację, przy której logarytmiczna charakterystyka amplitudowa układu
otwartego przecina oś odciętych, L(w�o)=0.
Ostatnim kryterium określającym jakość układu regulacji jest kryterium
całkowe. Pozwala ono ocenić jakość regulacji w stanie ustalonym (dokładność
statyczna), jak i w stanie nieustalonym (zapas stabilności i szybkość działania
układu). Za całkowe kryteria jakości regulacji przyjmuje się funkcjonały typu:
tk
&� &�&�
I =� f (�t,e,e,e,...)�dt (24)
��
0
gdzie f jest funkcja, e=e(t)- uchybem regulacji, a tk=Ą�.
9
Najczęściej za całkowe kryteria jakości regulacji przyjmuje się
następujące przypadki szczególne funkcjonału (24):
tk
I20 =� e2 dt (25)
��
0
tk
p
I2 p =� e2 dt ( p -� dana liczba naturaIna) (26)
��t
0
tk
2
&�
I2 =� (�e2 +� T e2)�dt (T -� dany wspólczynnik) (27)
��
0
tk
I0 =� e dt (28)
��
0
tk
p
I2 p =� e dt ( p -� dana liczba naturaIna) (29)
��t
0
tk
p
I2 p =� e dt ( p -� dana liczba naturaIna) (30)
��t
0
1.4. ANALIZA UKA
ADÓW DYNAMICZNYCH W ŚRODOWISKU MATLAB-
SIMULINK
1.4.1. Postaci modeli dynamicznych wykorzystywanych w pakiecie
Control Toolbox środowiska MATLAB
Rys.6. Okno programu MATLAB
W pakiecie tym wykorzystywane są następujące postaci liniowych modeli
dynamicznych:
�� Równania stanu;
�� Macierze transmitancji dla układów SIMO (jedno wejście wiele
wyjść);
10
�� Macierze transmitancji (układy SIMO) w postaci iloczynu zer,
biegunów i wzmocnienia.
Transmitancja układu dynamicznego jest funkcją wymierną operatora s
(dla układów ciągłych) lub z (dla układów dyskretnych). W bibliotece Control
Toolbox transmitancję podaje się w postaci pary wektorów zawierających
współczynników licznika i mianownika, przy czym umieszcza się je tam wg.
malejących potęg operatorów s lub z, np.: L=[1 2]; M=[1 3 2];
s +� 2
Takiemu zapisowi odpowiada transmitancja: G(�s)�=� . Aby określić
s2 +� 3s +� 2
transmitancję układu SIMO należy:
�� Podać wektor współczynników mianownika transmitancji;
�� Podać macierz zawierającą w kolejnych wierszach współczynniki
liczników odpowiadających kolejnym wyjściom układu.
Na przykład, zapis L=[1 2; 3 1]; M=[1 3 2] odpowiada transmitancji:
s +� 2
�� ł�
2
ę�s +� 3s +� 2ś�
G(�s)� =� .
ę� ś�
3s +�1
ę� ś�
��s2 +� 3s +� 2��
W programie MATLAB posługujemy się wyłącznie transmitancjami
spełniającymi warunek mówiący, ze stopień mianownika powinien być większy
od stopnia licznika. Aby określić macierz transmitancji przez podanie zer,
biegunów i wzmocnienia należy podać:
�� Kolumnowy wektor biegunów transmitancji P (ang. poles);
�� Macierz Z (ang. zeros) zawierającą w kolejnych kolumnach zera
odpowiadające kolejnym wyjściom układu;
�� Podać wektor kolumnowy G (ang. gains) zawierający wzmocnienia
odpowiadające kolejnym wyjściom układu.
Opisując układ dynamiczny w programie MATLAB należy pamiętać, że
najkorzystniejsze własności ma przedstawienie układu w postaci równań
stanu. Wykorzystujące podczas obliczeń postać transmitancji operaatorowej
układu program staje się bardziej podatny na błędy. Dotyczy to zwłaszcza
transmitancji wysokich rzędów lub o bardzo dużych lub bardzo małych
współczynnikach.
W programie tym przy konstruowaniu modeli dynamicznych wykorzystuje się
następujące polecenia:
�� Zmiana postaci modelu:
Skrót Zapis Opis
ss2tf [N, D]=ss2tf(A,B,C,D,iu) Funkcja ta zamienia równania stanu na odpowiadającą mu
transmitancję liczoną względem wejścia o numerze iu.
Macierz N zawiera w kolejnych wierszach współczynniki
liczników transmitancji odpowiadające kolejnym wyjściom
układu. Wektor D zawiera współczynniki mianownika
transmitancji.
tf2ss [A, B, C, D]=tf2ss(N, D) Funkcja ta wykonuje operacje odwrotne do funkcji ss2tf.
ss2zp [Z,P,K]=ss2zp(A,B,C,D,i Funkcja ta zamienia równania stanu na odpowiadającą im
u) macierz transmitancji liczoną względem wejścia o numerze
iu i przedstawia ją w postaci zer Z, biegunów P i
wzmocnień K.
zp2ss [A, B, C, D]=zp2ss(Z, P, Funkcja ta wykonuje operacje odwrotne do funkcji ss2zp.
K)
tf2zp [Z,P,K]=tf2zp(N, D) Funkcja ta znajduje zera Z, bieguny P i wzmocnienia K dla
11
transmitancji lub macierzy transmitancji układu SIMO
określonej parametrami N- współczynniki licznika, D-
współczynniki mianownika.
zp2tf [N, D]=zp2tf(Z, P, K) Funkcja ta wykonuje operacje odwrotne do funkcji zp2tf.
�� Zmiana struktury równań:
Skrót Zapis Opis
appen [A,B,C,D]= Funkcja ta dołącza do układu opisanego
d append(A1,B1,C1 D1,A2,B2, C2, równaniami stanu określonymi macierzami A1, ...
D2) układ o równaniach stanu wyznaczonymi przez
macierze A2,... .
augst [Ab, Bb, Cb, Db]=augstate(A,B,C,D) Funkcja ta wykonuje operacje odwrotne do funkcji
ate ss2tf.
ss2zp [Z,P,K]=ss2zp(A,B,C,D,iu) Funkcja ta rozszerza wektor wyjść układu o
współrzędne stanu. Zmiana ta wpływa wyłącznie
na równanie wyjścia układu.
ssdel [A, B, C, D]= Funkcja ta usuwa zadane wejścia WEJ, wyjścia
ete ssdelete(A1,B1,C1,D1,WEJ,WYJ,S WYJ i stany STANY z równań stanu o macierzach
TANY) A1,... .
sssele [A, B, C, D]= Funkcja ta umożliwia utworzenie podukładu na
ct ssdelete(A1,B1,C1,D1,WEJ,WYJ,S podstawie układu zadanego przez macierze
TANY) A1,....Parametry WEJ, WYJ, STANY są wektorami
zawierającymi indeksy wejść, itd., które mają
zostać włączone do tworzenia podukładu.
�� Tworzenie modeli dla różnych struktur układu:
Skrót Zapis Opis
cloop [A,B,C,D]=cloop(A1,B1,C1,D1,ZNAK) Funkcja ta oblicza macierze stanu układu
[A,B,C,D]=cloop(A1,B1,C1,D1,WEJ,WYJ) zamkniętego sprzężeniem zwrotnym o
[L,M]=cloop(L1,M1, ZNAK) wzmocnieniu 1. ZNAK określa znak
sprzężenia zwrotnego
feedba [A,B,C,D]= Funkcja ta oblicza macierze stanu układu
ck feedback(A1,B1,C1,D1,A2,B2,C2,D2,ZN zamkniętego sprzężeniem zwrotnym
AK) spisanym macierzami A2,.. i macierzami w
[A,B,C,D]= torze głównym A1,... ZNAK określa znak
feedback(A1,B1,C1,D1,A2,B2,C2,D2,WE sprzężenia zwrotnego
J,WYJ)
[L,M]= feedback (L1,M1,L2,M2,ZNAK)
paralle [A,B,C,D]=parallel(A1,B1,C1,D1,A2,B2,C Funkcja ta oblicza macierze stanu dla
l 2,D2) obiektów połączonych równolegle.
[A,B,C,D]=
feedback(A1,B1,C1,D1,A2,B2,C2,D2,
WEJ1,WYJ1,WEJ2,WYJ2)
[L,M]= parallel(L1,M1,L2,M2)
series [A,B,C,D]= series Funkcja ta oblicza macierze stanu dla
(A1,B1,C1,D1,A2,B2,C2,D2) obiektów połączonych szeregowo.
[A,B,C,D]= series
(A1,B1,C1,D1,A2,B2,C2,D2,
WYJ1,WEJ2)
[L,M]= series (L1,M1,L2,M2)
12
�� Inne funkcje wspomagające konstruowanie modeli dynamicznych:
Skrót Zapis Opis
rmode [L,M]=rmodel(n) Funkcje te generują losowe modele liniowe i dyskretne.
l [L,M]=rmodel(n,p) n-rząd modelu, p- ilość wyjść modelu, m-ilość wejść modelu
dmod [A,B,C,D]=rmodel(n)
el [A,B,C,D]=rmodel(n,p,
m)
Ord2 [A,B,C,D]=ord2(wn,dze Funkcja ta generuje odel układu oscylacyjnego II rzędu.
ta)
wn- w�n- częstotliwośc rezonansowa;
[L,M]=ord2(wn,dzeta)
dzeta-x�- współczynnik tłumienia.
Pade [A,B,C,D]=pade(T,n) Funkcja ta zwraca macierze równań stanu lub transmitancję
[L,M]= pade(T,n) układu przybliżającego opóz
nienie o dany czas T.
Przybliżenie to jest realizowane w oparciu o rozwinięci funkcji
w szereg Taylora.
1.4.2. Charakterystyki czasowe i częstotliwościowe układów
W programie MATLAB możemy wygenerować odpowiedz
układu w dziedzinie
czasu na następujące wymuszenia:
�� Odpowiedz
na skok jednostkowy;
�� Odpowiedz
impulsowa;
�� Odpowiedz
układu na stan początkowy- przebiegi stanu i wyjścia
jakie powstaną w układzie przy braku pobudzenia z zewnątrz i
danych warunkach początkowych;
�� W dziedzinie częstotliwości możemy wykreślić następujące
charakterystyki układu dynamicznego:
�� Charakterystykę Bode;
�� Charakterystykę Nyquista;
Powyższe zadania realizują następujące funkcje:
Skrót Zapis Opis
impuls impulse(A,B,C,D,iu, Funkcje te wyznaczają odpowiedz
układu ciągłego (dimpulse-
e t) dyskretnego) na impuls jednostkowy. W przypadku układu
dimpul impulse(L,M,t) opisanego równaniami stanu kreślona jest odpowiedz
wszystkich
se wyjść na impuls pojawiający się na wejściu o numerze iu. Dla
układów ciągłych można dołączyć własny wektor czasu symulacji
odpowiedzi. W przypadku pominięcia parametru t wektor czasu
zostanie pominięty.
initial initial(A,B,C,D,X0,T Funkcja ta wyznacza odpowiedz
układu opisanego równaniami
dinitial ) stanu na warunek początkowy określony parametrem X0. Czas
symulacji i wektor chwil czasu jest określany automatycznie lub
przez użytkownika (parametr T).
step step(A,B,C,D,iu) Funkcja ta wyznacza odpowiedz
układu ciągłego na skok
dstep step(A,B,C,D,iu,t) jednostkowy. W przypadku układu opisanego równaniami stanu
step(L,M) kreślona jest odpowiedz
wszystkich wyjść na skok jednostkowy
pojawiający się na wejściu o numerze iu. Dla układów ciągłych
można dołączyć własny wektor czasu symulacji odpowiedzi. W
przypadku pominięcia parametru t wektor czasu zostanie
pominięty.
bode bode(A,B,C,D,iu,w) Funkcja ta wyznacza charakterystykę częstotliwościową Bode.
bode(L,M,w) Częstotliwości i liczba punktów dobierana jest automatycznie lub
zgodnie z wymaganiami użytkownika (parametr w). W przypadku
układu opisanego równaniami stanu kreślona jest odpowiedz

wszystkich wyjść na skok jednostkowy pojawiający się na wejściu
o numerze iu.
13
nyquist nyquist(A,B,C,D,iu, Funkcja ta wyznacza charakterystykę Nyquista. Częstotliwości i
w) liczba punktów dla których kreślona jest charakterystyka dobierane
nyquist(L,M,w) są automatycznie lub zgodnie z wymaganiami użytkownika
(parametr w). W przypadku układu opisanego równaniami stanu
kreślona jest odpowiedz
wszystkich wyjść na skok jednostkowy
pojawiający się na wejściu o numerze iu.
1.4.3. Analiza właściwości dynamicznych obiektu sterowania.
Analiza właściwości obiektu sterowania polega na określeniu stabilności
układu, czy jest on obserwowalny, czy też nie, obliczeniu zapasu amplitudy i
fazy obiektu jak również wykreśleniu dróg pierwiastków układu. Zadania te
realizują następujące funkcje:
Skrót Zapis Opis
margin [Gm, Pm, Wcg, Wcp] = Funkcja ta oblicza zapas amplitudy (Gm) i fazy
margin(A,B,C,D) (Pm) dla układu opisanego równaniami stanu
[Gm, Pm, Wcg, Wcp] = margin(L,M) lub transmitancją oraz odpowiadające im
częstotliwosci graniczne- odpowiednio Wcg i
Wcp.
imargi [Gm, Pm, Wcg, Wcp] = Funkcja ta działa podobnie jak funkcja margin z
n margin(ampl,faza,w) tym, że obliczenia opierają się o zadane
wektory amplitudy i fazy oraz odpowiadające im
częstotliwości, które można uzyskać funkcją
bode lub dbode, dzięki czemu funkcja ta działa
zarówno dla układów ciągłych, jak i
dyskretnych.
rlocus rlocus(A,B,C,D) Funkcja ta wyznacza drógi pierwiastków układu
rlocus(L,M) otwartego.
rlocus (L,M,K)
rlocus(A,B,C,D,K)
rlocfin [K,R]=rlocfind(A,B,C,D) Jeżeli w aktywnym oknie znajduje się wykres
d [K,R]=rlocfind (L,M) zależności położenia biegunów od wzmocnienia
[K,R]=rlocfind (L,M,n) uzyskany funkcją rlocus, to funkcja rlocfind
[K,R]=rlocfind (A,B,C,D,n) umożliwia wybranie myszką żądanego bieguna
układu SISO opisanego transmitancją lub
równaniami stanu.
1.5. SIMULINK
Wygodnym i skutecznym narzędziem symulacyjnym oferowanym przez
środowisko MATLAB jest program SIMULINK. Umożliwia on utworzenie w
oknie graficznym struktury układu sterowania zbudowanej z bloków różnych
typów reprezentujących obiekty dynamiczne, z
ródła sygnału i przyrządy
pomiarowe. Definiując obiekty możemy odwołać się do istniejących w pamięci
zmiennych, dostępnych normalnie z wiersza poleceń MATLAB a. Po
skonstruowaniu struktury i określeniu parametrów poszczególnych elementów
ustalamy parametry symulacji i wybieramy polecenie Simulation/Start. Wyniki
symulacji zostaną przedstawione w standardowych oknach graficznych
programu MATLAB.
14
Rys.7. Okno programu SIMULINK
W programie tym mamy do dyspozycji nastepujące grupy elementów
dynamicznych (rys.7.):
�� y
ródła (ang. sources)- np. skok jednostkowy, generator sygnałów
zmiennch;
�� Elementy wizualizacji wyników symulacji (ang. sinks)- np.
oscyloskop, itd.;
�� Elementy ciągłe (ang. continuous);
�� Elementy dyskretne (ang. discrete);
�� Elementy realizujące określone zadania matematyczne (ang.
math);
�� Funkcje i tabele (ang. functions & tables);
�� Elementy nieliniowe (ang. nonlinear);
�� Inne elementy niezbędne do zbudowania układu symulacji. Blok
Signals & Systems- np. multipleksery i demultipleksery, wejścia i
wyjścia układu, itd.
Rozwinięcie bloków elementów z rys.7. przedstawiono na rys.8. i rys.9.
Rys.8. Biblioteki: Sources, Signal & Systems, Math.
15
Program ten ściśle współpracuje z środowiskiem MATLAB. Dlatego też,
pozwala on na analizę układu zamodelowanego w programie SIMULINK za
pomocą funkcji linmod oraz dlinmod o nastepującym zapisie: [A, B, C,
D]=linmod( nazwa ). W miejsce  nazwa należy wpisać nazwę pliku
zawierającego odpowiednio przygotowany model analizowanego układu.
Właściwe przygotowanie modelu polega na wskazaniu w programie
SIMULINK wejścia i wyjścia układu (rys.9.).
Rys.9. Biblioteki progarmu SIMULINK: Continuous, Discrete, Nonlinear,
Functions & Tables, Sinks
16
1.6. PRZEBIEG ĆWICZENIA LABORATORYJNEGO
1.6.1. Opis stanowiska laboratoryjnego
Stanowisko laboratoryjne składa się komputera klasy PC wraz z
oprogramowaniem MATLAB- SIMULINK.
1.6.2. Wykonanie ćwiczenia
1) Badasz układ opisany transmitancją operatorową postaci:
k
G(�s)� =�
2
T s2 +� 2x�Ts +�1
dla danych współczynników:
a) k=1.3, T=1, � =1.9
b) k=1.2, T=1, � =2
c) k=1.4, T=1, � =1.8
d) k=1.1, T=1, � =1.7
e) k=1.4, T=1, � =2.1
f) k=1.2, T=1, � =1.9
g) k=1.3, T=1, � =1.7
h) k=1.5, T=1, � =1.9
2) Korzystając z danych współczynników obliczyć na kartce postać
transmitancji G(s)
3) W programie Matlab zdefiniować transmitancję badanego układu:
a) zdefiniować licznik np.: dla równania 3s2 +� 5 zapisać:
L = [3 0 5]
b) zdefiniować mianownik np.: dla równania 2s3 +� 4s2 +� 7 zapisać:
M = [2 4 0 7]
c) funkcją tf(licznik, mianownik) definiujemy transmitancję np.:
dla powyższych danych: G = tf (L, M)
4) Wyznaczamy charakterystykę skokową układu opisanego
transmitancją G. Zapisujemy: step(G).
a) Zapisujemy: step(G).
b) Na wykresie klikamy prawym klawiszem myszy
characteristics i wybieramy wartości charakterystyczne, które
chcemy wyświetlić. np.:
i) peak response  wartość maksymalna
ii) settling time  czas regulacji
iii) rise time  czas narastania
iv) steady state  stan ustalony
c) komentujemy odczytane wartości
5) Wyznaczamy charakterystykę impulsową układu opisanego
transmitancją G.
a) Zapisujemy: impulse(G).
17
b) Na wykresie klikamy prawym klawiszem myszy
characteristics i wybieramy wartości charakterystyczne, które
chcemy wyświetlić. np.:
i) peak response  wartość maksymalna
ii) settling time  czas regulacji
c) komentujemy odczytane wartości
6) Wyznaczamy charakterystyki amplitudową i fazową układu
opisanego transmitancją G (wykresy Bodego).
a) Zapisujemy: bode(G).
b) Na wykresie klikamy prawym klawiszem myszy
characteristics i wybieramy wartości charakterystyczne, które
chcemy wyświetlić. np.:
i) peak response  wartość maksymalna
ii) All stability margins  zapas amplitudy i fazy
c) komentujemy odczytane wartości
7) Wyznaczamy charakterystykę amplitudowo fazową układu
opisanego transmitancją G (wykres Nyquista).
a) Zapisujemy: nyquist(G).
b) Na wykresie klikamy prawym klawiszem myszy
characteristics i wybieramy wartości charakterystyczne, które
chcemy wyświetlić. np.:
i) peak response  wartość maksymalna
ii) All stability margins  zapas amplitudy i fazy
c) komentujemy odczytane wartości
8) Zamykamy układ pętlą ujemnego sprzężenia zwrotnego ze
wzmocnieniem:
wykorzystując funkcję:
[Lf, Mf] = feedback(L, M, L1, M1, znak sprzężenia),
gdzie:
- L  licznik układu w torze głównym
- M  mianownik układu w torze głównym
- L1  licznik układu w torze sprzężenia zwrotnego
- M1  mianownik układu w torze sprzężenia zwrotnego
- Lf  licznik układu zamkniętego
- Mf  mianownik układu zamkniętego
np.:
dla wzmocnienia k = 3 i ujemnego sprzężenia zwrotnego:
18
 [Lf, Mf] = feedback(L, M, [3], [1], -1)
9) Zdefiniować transmitancję układu zamkniętego zapisując:
10) F = tf (Lf, Mf)
11) Wyznaczyć charakterystyki czasowe układu zamkniętego F wg
punktów 4 i 5
12) Porównać wartości charakterystyczne dla układów otwartych i
zamkniętych.
13) Dobrać wzmocnienie w torze sprzężenia zwrotnego w celu
minimalizacji wybranej wartości charakterystycznej układu
zamkniętego np. czasu regulacji.
14) Ćwiczenie wykonać dla różnych współczynników transmitancji
układu w torze głównym wymienionych w pkt. 1 a) - h).
1.7. WYKONANIE SPRAWOZDANIA
W sprawozdaniu powinno zostać zawarte:
�� Model badanego układu dynamicznego;
�� Wydruk M-pliku;
�� Wydruk charakterystyk dynamicznych układu;
�� Wnioski do otrzymanych wyników.
1.8. LITERATURA:
1. B. Mrozek, Z. Mrozek,  MATLAB. Uniwersalne środowisko do obliczeń
naukowo-technicznych , Kraków 1995.
2. T. Kaczorek  Teoria sterowania , Warszawa 1977
3. A. Zalewski, R. Cegieła  Matlab- obliczenia numeryczne i ich
zastosowania , Poznań 1996.
19