1 

WYDZIAŁ TRANSPORTU PW                                        Zakład Sterowania Ruchem 

LABORATORIUM   PODSTAW  AUTOMATYKI III 

 Rok akad. 

2009/2010

  Data wykonania ćwiczenia 

2010r.  Uzyskane 

punkty 

za:

 

Specjalizacja 

Imiona  i  Nazwiska  studentów 

przygotowanie 

i realizację  

zaliczenie 

poprawę 

zaliczenia 

. . . . . . . . . .  

1. 

 

 

 

Zespół 

2. 

 

 

 

 

3. 

 

 

 

 

4. 

 

 

 

 

5. 

 

 

 

 

Ćwiczenie  nr  1 

Temat ćwiczenia: Badanie obiektów dynamicznych z zastosowaniem programu 

MATLAB 

 

1 

Wprowadzenie 

 
Dany jest schemat obwodu elektrycznego wzbudzenia przekaźnika o indukcyjności L i oporze 

czynnym RP (rysunek 1). W szereg z przekaźnikiem dołączony jest kondensator o pojemności C. 
Przewody  łączące  tak  powstały  dwójnik  ze  źródłem  zasilania  mają  opór  R

d

.  Zakładamy,  że 

indukcyjność  cewki  i  rdzenia  pozostaje  stała  w  czasie  (normalnie  ulega  zmianie  po  rozpoczęciu 
ruchu kotwicy). 

L

R

d

R

L

i

1

(t)

C

u(t)

i(t)

i

2

(t)

 

Rysunek 1. Obwód przekaźnika. 

 

Korzystając  z  równań  różniczkowych opisujących poszczególne elementy obwodu, oraz praw 

Ohm'a i Kirchoff'a, ułożyć można równanie różniczkowe zwyczajne drugiego rzędu (1) opisujące 
zadany dwójnik i stanowiące podstawę do dalszych rozważań. 

 

 

 

 

LC

R

t

u

t

i

LC

R

R

R

dt

t

di

C

R

L

R

dt

t

i

d

d

d

d

L

d

L

























1

2

2

 

(1) 

Przyjmując  za  sygnał  wejściowy  napięcie  u(t)  przyłożone  do  wyprowadzeń  dwójnika,  a  za 

sygnał  wyjściowy  prąd  płynący  przez  cewkę,  uzyskujemy  liniowy  układ  dynamiczny  o  jednym 
wejściu u(t) i jednym wyjściu i(t). Traktując obie strony równania (1) odwrotnym przekształceniem 
Laplace'a,  oraz  przyjmując  zerowe  warunki  początkowe  (kondensator  w  chwili  t  =  0

+

  jest 

rozładowany  i w żadnej  z  gałęzi  nie  płynie  prąd)  uzyskujemy  równanie  transformat  (2) 
rozpatrywanego układu. 

 

2 

 

 

 

 

LC

R

s

u

s

i

LC

R

R

R

s

si

C

R

L

R

s

i

s

d

d

d

L

d

L

























1

2

 

(2) 

Odpowiednie podzielenie stronami równania (2) sprowadza je do postaci (3) będącej zarazem

 

transmitancją operatorową układu. 

 

 

 





































LC

R

R

R

s

C

R

L

R

s

LC

R

s

u

s

i

s

G

d

d

L

d

L

d

1

1

2

 

 

 





d

L

d

L

d

R

R

s

L

C

R

R

LCs

R

s

G











2

1

 

(3) 

W zależności od wartości współczynników wielomianu mianownika wyrażenia (3), układ można 

utożsamiać  z  członem  inercyjnym  drugiego  rzędu,  bądź  ze  stabilnym  członem  oscylacyjnym 
(człon  niestabilny  posiada  dwa  ostatnie  współczynniki  o  wartości  ujemnej,  co  w  rozpatrywanym 
obwodzie  nie  zachodzi  –  wartości  rezystancji,  pojemności  i  indukcyjności  są  rzeczywistymi 
liczbami dodatnimi). 
W przypadku rzeczywistych pierwiastków równania charakterystycznego (4) 





0

2











d

L

d

L

d

R

R

s

L

C

R

R

LCs

R

 

(4) 

układ  jest  członem  inercyjnym  drugiego  rzędu,  w  przypadku  zespolonych  bądź  urojonych  – 
stabilnym członem oscylacyjnym. 
Wykonanie  ćwiczenia  polega  na  wyznaczeniu  określonych  charakterystyk  opisywanego  układu 
dla  zadanych  przez  prowadzącego  parametrów  za  pomocą  narzędzi  zaimplementowanych 
w środowisku MATLAB. 
 
2 

Zasady współpracy z programem MATLAB (niezbędne do wykonania ćwiczenia) 

 

Wykonanie  ćwiczenia  polega  na uruchomieniu programu MATLAB (wersja 5.3., dydaktyczna) 

przez  kliknięcie  na  ikonę  tego  programu  na  pulpicie  WINDOWS.  Po  zgłoszeniu  się  programu 
następuje  otwarcie  okna  poleceń  MATLAB  Command  Window  (gotowość  sygnalizowana 
kursorem:    >>    ).  Należy  wpisać  następnie  polecenie 

clear;clc 

i  zatwierdzić  je  naciskając 

klawisz Enter. Spowoduje to wyczyszczenie pola poleceń, zapamiętanych zmiennych i ustawienie 
kursora w lewym górnym rogu pola. Z menu File wybiera się polecenie Open a następnie należy 
wybrać katalog MATLABzad w którym wskazuje się na określony plik tekstowy typu skryptowego 
(tzn.  z rozszerzeniem  *.m,  znamiennym  dla  MATLAB’a).  Zbiór  tych  plików  umieszczonych 
w katalogu MATLABzad zawiera następujące pliki typu skryptowego: 

 

Każdy  z  tych  plików  zawiera  ciągi  poleceń  wykonujących  algorytm  obliczeń  i  ew.  kreślenia 

charakterystyk  lub  funkcji  napisanych  w  języku  programowania  wysokiego  poziomu 
akceptowanym  przez  środowisko  MATLAB’a.  Kliknięcie  na  określony  plik  i  akceptacja  otwarcia 
tego  pliku  spowodują,  że  automatycznie  zgłosi  się  program  MATLAB  Editor/Debugger,  który 

 

3 

otworzy  swoje  własne  okno  i  wyświetli  zawartość  wskazanego  pliku.  W  trybie  MATLAB 
Editor/Debugger dokonuje się edycji, poprawek i modyfikacji programów standardu MATLAB. 

Aby wykonać program zapisany w wybranym pliku należy (jeden z wielu sposobów): 

1.  Z menu wybrać opcje Edit a następnie  Select all i skopiować do pamięci poleceniem Copy 

(lub klikając w pasku poleceń na ikonę kopiowania), 

2.  Otworzyć  okno  MATLAB  Command  Window  i  poleceniem  Paste  z  menu  Edit  wprowadzić 

ciąg  poleceń  z  wybranego  pliku  tekstowego  do  obszaru  poleceń  MATLAB  Command 
Window, 

3.  Uruchomić wpisany w obszar poleceń program - przez akceptację przycisku Enter , 
4.  Jeśli wynikiem programu jest wykreślenie funkcji lub charakterystyk, następuje otwarcie okna 

graficznego Figure No. 1 z wykreślonymi charakterystykami lub funkcjami, jeśli na tle rysunku 
pojawią się dwie prostopadłe osie, oznacza to, że przeciągając mysz do wybranego miejsca 
na rysunku (charakterystyce), należy kliknąć w tym miejscu aby wstawić opis odnoszący się 
do rysunku (charakterystyki), 

5.  Jeśli  wynikiem  są  obliczenia,  następuje  wyświetlenie  wyników  w  formacie  zadanym  przez 

program  z  pliku  skryptowego  na  ekranie  monitora  bezpośrednio  po  sekwencji  poleceń 
odnoszących się do tego programu, 

6.  Po  uzyskaniu  wyników  wybranego  programu  można  uruchomić  program  następnego  pliku 

skryptowego otwierając (powrót do okna) okno MATLAB Editor/Debugger i zamykając otwarty 
uprzednio  w  tym  oknie  plik  poleceniem  File/Close  oraz  okno  rysunku  Figure  No  1,  po 
zamknięciu  okna  rysunku  i  zamknięciu  pliku  otworzyć  żądany  plik  poleceniem  File/Open 
i powtórzyć czynności od p.1., 

7.  Wyniki a w tym i rysunki należy wyprowadzić na drukarkę, wydruki opisać, 
8.  Aby zakończyć działanie MATLAB’a zamyka się kolejno okna Figure No. 1 (jeśli było otwarte), 

MATLAB Editor/Debugger i na końcu MATLAB Command Window w wyniku czego następuje 
powrót do systemu operacyjnego Windows. 

 

3 

Wykonanie ćwiczenia  

 
Realizacja  ćwiczenia  wymaga  kolejnego  wykonania  wszystkich  programów  zapisanych 

w plikach  od  P3_1_char_skokowe.m  do  P3_6_wsp_wielom_char.m.  Poniżej  zamieszczono  opis 
realizowanych  funkcji,  poparty  przykładami  m-plików  zrealizowanymi  dla  wybranych  parametrów 
elektrycznych  obwodu.  Zamieszczone  wyniki  powinny  zostać  wykorzystane  podczas odpowiedzi 
na postawione pytania. 
 
3.1  Wyznaczenie charakterystyk skokowych. 

 
Jedną z podstawowych charakterystyk opisujących układ dynamiczny jest jego odpowiedź na 

skok jednostkowy 1(t) zadany na wejście u(t). 

Skoro  u(t)  =  U·1(t),  toteż  transformatą  tego  sygnału  będzie  wymuszenie  u(s)  =  U·s

-1

. 

Współczynnik  U  będzie  wartością  napięcia,  jaką  podano na dwójnik. Korzystając z transmitancji 
operatorowej  (3)  wyznaczamy  odpowiedź  układu  w  dziedzinie  zmiennej  zespolonej  s  ze  wzoru 
(4). 

 

   









s

R

R

s

L

C

R

R

LCs

R

U

s

u

s

G

s

i

d

L

d

L

d















2

3

 

(4) 

Przekształcenie  operatorowe  Laplace’a  powyższego  wyrażenia  stanowi  wartość  natężenia 

prądu w funkcji czasu od chwili podania napięcia na zaciski wejściowe.  

Do  wyznaczenia  charakterystyki  skokowej  wykorzystamy  gotową  funkcję  step  znajdującą  się 

w dodatkowym pakiecie control narzędzi programu MATLAB. Postać funkcji jest następująca: 
[y, x, t]= step(L, M) 
gdzie: 
t – wektor kolejnych chwil czasu dla których obliczane są macierze wartości wektorów x i y, 

 

4 

y – macierz wartości sygnału wyjściowego w kolejnych chwilach symulacji, 
x – macierz wartości wektora stanu dla kolejnych chwil czasu, 
L – wektor współczynników wielomianu licznika transmitancji G(s), 
M - wektor współczynników wielomianu mianownika transmitancji G(s). 
 
Przykładowy  m.plik  powodujący  wykreślenie  charakterystyki  skokowej  przyjmie  zatem  postać 
(przyjęto U

1

 = L

1

 = C

1

 = R

d1

 = 1, R

L1 

= 4, U

2

 = L

2

 = 1, C

2

 = 0.2, R

d2

 = 5, R

L2 

= 1): 

 

%Przykład m-pliku; 
%JCK; 
%Charakterystyka skokowa obiektu o transmitancji; 
%           1; 
%G(s)=--------------; 
%      T2s^2+T1s+T0; 
%=================================================; 
U1=1; 
L1=1; 
C1=1; 
Rd1=1; 
Rl1=4; 
T21=Rd1*L1*C1; 
T11=Rl1*Rd1*C1+L1; 
T01=Rl1+Rd1; 
%wyznaczenie współczynników licznika transmitancji; 
l1=[0,0,U1]; 
%wyznaczenie współczynników mianownika transmitancji; 
m1=[T21,T11,T01]; 
% wyznaczenie wektora odpowiedzi; 
[y1,x1,t1]=step(l1,m1); 
%=================================================; 
U2=1; 
L2=1; 
C2=0.2; 
Rd2=5; 
Rl2=1; 
T22=Rd2*L2*C2; 
T12=Rl2*Rd2*C2+L2; 
T02=Rl2+Rd2; 
%wyznaczenie współczynników licznika transmitancji; 
l2=[0,0,U2]; 
%wyznaczenie współczynników mianownika transmitancji; 
m2=[T22,T12,T02]; 
% wyznaczenie wektora odpowiedzi; 
[y2,x2,t2]=step(l2,m2); 
%=================================================; 
% wykres obu odpowiedzi skokowych odpowiednio na jednym rysunku; 
plot(t1,y1,'g',t2,y2,'m') 
title('Charakterystyki czasowe (skokowe) obiektu'); 
ylabel('Amplituda odpowiedzi na wymuszenie jednostkowe 1(t)'); 
xlabel('czas [s]'); 
gtext('U1=1;L1=1;C1=1;Rd1=1;Rl1=4'); 
gtext('U2=1;L2=1;C2=0.2;Rd2=5;Rl2=1'); 
grid;pause;close;clear;clc; 

 
Uruchomienie powyższego m-pliku pozwala uzyskać następującą charakterystykę 

 

5 

 

 
3.2  Wyznaczenie charakterystyk amplitudowo-fazowych. 

 
Do  wykreślenia  charakterystyk  częstotliwościowych  układu  konieczna  jest  znajomość 

transmitancji  G(s)  liniowego  układu  ciągłego.  Na  podstawie  transmitancji  widmowej  G(jω) 
wyznacza  się  odpowiednio  wartości  (w  postaci  wektorów) części rzeczywistych P(ω) i urojonych 
Q(ω). Podstawmy do wzoru (3) s = jω. W efekcie otrzymamy 

 

 





d

L

d

L

d

R

R

j

L

C

R

R

j

LC

R

j

G

















2

1

 

 











L

C

R

R

j

LC

R

R

R

j

G

d

L

d

d

L











2

1

 

(5a) 

 





















2

2

2

2

2

1

1











L

C

R

R

LC

R

R

L

C

R

R

j

LC

R

R

R

j

G

d

L

d

L

d

L

d

d

L



















 

(5) 

Stąd część rzeczywista P(ω) i urojona Q(ω) przyjmują następujące postaci 

 

















2

2

2

2

2

1

1









L

C

R

R

LC

R

R

LC

R

R

R

j

P

d

L

d

L

d

d

L















 

(6) 

 

















2

2

2

2

1









L

C

R

R

LC

R

R

L

C

R

R

j

Q

d

L

d

L

d

L















 

(7) 

Dla wartości pulsacji ω (w poniższym programie zmienna oznaczająca pulsację występuje jako 

zmienna  w)  zmieniającej  się  od  0  do  100  rd/s  co  0,01.  Do  wykreślenia  charakterystyk 
amplitudowo-fazowych  zastosowano  standardową  funkcję  plot.  Funkcja  ta  ma  następującą 
postać:  

 

6 

plot(x, y, ‘typ linii’) 
gdzie: 
x – zbiór wartości na osi X (wartości zmiennej niezależnej) , 
y – zbiór wartości na osi Y (wartości funkcji), 
typ linii – znaki oznaczające kolor i rodzaj linii. 

 
%Przykład m-pliku; 
%JCK; 
%Charakterystyki częstotliwościowe - charakterystyki; 
%amplitudowo-fazowe obwodu przekaźnika; 
%           1; 
%G(s)=--------------; 
%      T2s^2+T1s+T0; 
%Wykreślenie ch-k polega na obliczeniu części wspólnej ws; 
%a następnie części rzeczywistej P i urojonej Q odpowiednio dla; 
%zadanych wartości parametrów elektrycznych obwodu a następnie; 
%dokonaniu wykreślenia przy pomocy funkcji plot; 
w=0:0.01:100; 
roz=size(w); 
%======================================================; 
L1=1; 
C1=1; 
Rd1=1; 
Rl1=4; 
ws1=1./(((ones(roz)-
(w.^2)*L1*C1)*Rd1+Rl1).^2+((Rl1*Rd1*C1+L1)^2).*(w.^2)); 
P1=((ones(roz)-(w.^2)*L1*C1)*Rd1+Rl1).*ws1; 
Q1=-((Rl1*Rd1*C1+L1)^2).*(w.^2).*ws1; 
%======================================================; 
L2=1; 
C2=0.2; 
Rd2=5; 
Rl2=1; 
ws2=1./(((ones(roz)-
(w.^2)*L2*C2)*Rd2+Rl2).^2+((Rl2*Rd2*C2+L2)^2).*(w.^2)); 
P2=((ones(roz)-(w.^2)*L2*C2)*Rd2+Rl2).*ws2; 
Q2=-((Rl2*Rd2*C2+L2)^2).*(w.^2).*ws2; 
%======================================================; 
%wykreślenie ch-k amplitudowo-fazowych 
plot(P1,Q1,'r',P2,Q2,'g'); 
 

%comet(P1,Q1,0.05);pause; 

 

%comet(P2,Q2,0.05);pause; 

title('Charakterystyki amplitudowo-fazowe obiektu dynamicznego'); 
xlabel('oś wartości rzeczywistych P1, P2'); 
ylabel('oś wartości urojonych Q1, Q2'); 
gtext('L1=1;C1=1;Rd1=1;Rl1=4'); 
gtext('L2=1;C2=0.2;Rd2=5;Rl2=1');grid;pause; 
%wyznaczenie współczynników licznika transmitancji; 
l1=[0,0,1]; 
l2=[0,0,1]; 
%wyznaczenie współczynników mianownika transmitancji; 
T21=Rd1*L1*C1; 
T11=Rl1*Rd1*C1+L1; 
T01=Rl1+Rd1; 
T22=Rd2*L2*C2; 
T12=Rl2*Rd2*C2+L2; 
T02=Rl2+Rd2; 

 

7 

m1=[T21,T11,T01]; 
m2=[T22,T12,T02]; 
nyquist(l1,m1,w);gtext('L1=1;C1=1;Rd1=1;Rl1=4');grid;hold on; 
nyquist(l2,m2,w);gtext('L2=1;C2=0.2;Rd2=5;Rl2=1');grid;pause; 
close;clear;clc; 

 
Alternatywną metodą do zastosowania funkcji plot jest zastosowanie funkcji comet. Funkcja ta 

pozwala  zaobserwować  animację  rysowania  wykresu  charakterystyki  amplitudowo-fazowej 
w miarę wzrostu pulsacji, i ma następującą postać: 
comet(x,y,dł) 
gdzie: 
x – odcięta trajektorii komety ruchomego wykresu imitującego jej lot, 
y – rzędna trajektorii komety ruchomego wykresu imitującego jej lot, 
dł – długość „ogona” komety w przedziale (0,1). 

Aby  uruchomić  tę  funkcję  w  m-pliku  należy  po  wklejeniu  usunąć  znaczki 

%

  sprzed  instrukcji 

comet, a umieścić jeden przed instrukcją plot. 

Wyznaczenie charakterystyk amplitudowo-fazowych (Nyquista) w zakresie pulsacji ω od 0 do 

+∞ i od -∞ do 0 odbywa się za pomocą standardowej funkcji nyquist o postaci: 
nyquist(L, M, w) 
gdzie: 
L  i  M  to  wektory  współczynników  wielomianów  licznika  –  L  i mianownika  –  M  transmitancji  G(s) 
o postaci 

G s

L

M

( ) 

 

(8) 

w - zakres pulsacji dla którego będzie kreślona charakterystyka. 
 
W wyniku uruchomienia opisywanego m-pliku uzyskuje się następujące charakterystyki 
 

 

 

8 

 

 

3.3  Charakterystyka amplitudowa i charakterystyka fazowa 

 
Amplitudą A(ω) nazywamy moduł transmitancji widmowej (5) o następującej postaci 

 













2

2

2

2

1

1







L

C

R

R

LC

R

R

R

A

d

L

d

d

L











 

(9) 

Faza  transmitancji  widmowej  stanowi  argument  tej  funkcji.  Wyliczamy  ją  na  podstawie  prostej 
zależności (10). 

 

 





 

 









2

1

arg















LC

R

R

R

L

C

R

R

arctg

P

Q

arctg

G

d

d

L

d

L













 

(10) 

Do wykreślenia charakterystyk amplitudowo-fazowych w funkcji pulsacji zastosujemy poznaną 

już funkcję plot.  

Do  wykreślenia  charakterystyk  logarytmicznych  modułu  i  fazy  zostanie  zastosowana 

standardowa funkcja bode o postaci: 
[moduł,faza,pulsacja]=bode(L, M) 
gdzie: 
L  i  M  to  wektory  współczynników  wielomianów  licznika  –  L  i  mianownika  –  M  transmitancji  G(s) 
o postaci (8). 

Pulsacja  ω  (w  poniższym  programie  zmienna  oznaczająca  pulsację  występuje  także  jako 

zmienna w) i zmienia się od 0 do 100 [rd/s], co 0,01 [rd/s]. 
 

%Przykład m-pliku; 
%JCK; 
% Charakterystyki częstotliwościowe - amplitudowe i fazowe; 
%           1; 
%G(s)=--------------; 
%      T2s^2+T1s+T0; 

 

9 

%Wykreślenie ch-k polega na obliczeniu części wspólnej ws; 
%a następnie części rzeczywistej P i urojonej Q odpowiednio; 
%dla zadanych wartości parametrów elektrycznych obwodu; 
%a następnie dokonaniu wykreślenia przy pomocy funkcji plot; 
w=0:0.01:100; 
roz=size(w); 
%======================================================; 
L1=1; 
C1=1; 
Rd1=1; 
Rl1=4; 
ws1=1./(((ones(roz)-
(w.^2)*L1*C1)*Rd1+Rl1).^2+((Rl1*Rd1*C1+L1)^2).*(w.^2)); 
P1=((ones(roz)-(w.^2)*L1*C1)*Rd1+Rl1).*ws1; 
Q1=-((Rl1*Rd1*C1+L1)^2).*(w.^2).*ws1; 
%======================================================; 
L2=1; 
C2=0.2; 
Rd2=5; 
Rl2=1; 
ws2=1./(((ones(roz)-
(w.^2)*L2*C2)*Rd2+Rl2).^2+((Rl2*Rd2*C2+L2)^2).*(w.^2)); 
P2=((ones(roz)-(w.^2)*L2*C2)*Rd2+Rl2).*ws2; 
Q2=-((Rl2*Rd2*C2+L2)^2).*(w.^2).*ws2; 
%======================================================; 
%wykreślenie ch-k amplitudowych i fazowych; 
plot(w,P1,'r',w,Q1,'r',w,P2,'g',w,Q2,'g'); 
ylabel('Q1, Q2, P1, P2'); 
xlabel('pulsacja'); 
title('Charakterystyki amplitudowe - części Re i Im w funkcji 
pulsacji'); 
gtext('U1=1;L1=1;C1=1;Rd1=1;Rl1=4'); 
gtext('U2=1;L2=1;C2=0.2;Rd2=5;Rl2=1'); 
grid;pause; 
%======================================================; 
T21=Rd1*L1*C1; 
T11=Rl1*Rd1*C1+L1; 
T01=Rl1+Rd1; 
%wyznaczenie współczynników licznika transmitancji; 
l1=[0,0,1]; 
%wyznaczenie współczynników mianownika transmitancji; 
m1=[T21,T11,T01]; 
T22=Rd2*L2*C2; 
T12=Rl2*Rd2*C2+L2; 
T02=Rl2+Rd2; 
%wyznaczenie współczynników licznika transmitancji; 
l2=[0,0,1]; 
%wyznaczenie współczynników mianownika transmitancji; 
m2=[T22,T12,T02]; 
bode(l1,m1);grid;hold on; 
bode(l2,m2);grid;hold on; 
title('Charakterystyki: amplitudowa (logarytmiczna) i fazowa w 
funkcji częstotliwości'); 
gtext('U1=1;L1=1;C1=1;Rd1=1;Rl1=4');grid; 
gtext('U2=1;L2=1;C2=0.2;Rd2=5;Rl2=1'); 
grid;pause; 
close;clear;clc; 

 

10 

 

W wyniku uruchomienia opisywanego m-pliku uzyskuje się następujące charakterystyki 

 

 

 

 

 

3.4  Graficzne  kryterium  stabilności  Michajłowa  na  podstawie  hodografu  oraz  przy  użyciu 

charakterystyk częstotliwościowych 

 

 

11 

Zgodnie z kryterium Michajłowa układ n-tego rzędu, określony równaniem charakterystycznym, 

będącym mianownikiem M(s) transmitancji operatorowej G(s) jest stabilny wtedy i tylko wtedy, gdy 
zmiana  argumentu  wektora  M(jω)  przy  zmianie  parametru  ω  od  0  do  +∞  wyniesie  n(/2),  czyli 
wykres  (hodograf)  M(jω)  przejdzie  przez  n  kwadrantów  (ćwiartek  układu  współrzędnych).  Do 
wykreślenia hodografu zastosujemy funkcję plot i obliczywszy wcześniej część rzeczywistą P(ω) 
i urojoną Q(ω) wektora M(jω). 

Poszukiwany wektor M(jω) jest odwrotnością wyrażenia (5a). 
Uruchomienie poniższego m-pliku pozwala uzyskać następującą charakterystykę 
 

 

 

%Przykład m-pliku; 
%JCK; 
%Wyznaczenie  stabilności  układu  dynamicznego  za pomoca graficznego 
kryterium Michajlowa; 
%           1; 
%G(s)=--------------; 
%      T2s^2+T1s+T0; 
w=0:0.01:10; 
roz=size(w); 
%======================================================; 
L1=1; 
C1=1; 
Rd1=1; 
Rl1=4; 
P1=Rl1+Rd1-(w.^2)*L1*C1*Rd1; 
Q1=(Rl1*Rd1*C1+L1).*w; 
L2=1; 
C2=0.2; 
Rd2=5; 
Rl2=1; 
P2=Rl2+Rd2-(w.^2)*L2*C2*Rd2; 

 

12 

Q2=(Rl2*Rd2*C2+L2).*w; 
plot(P1,Q1,'r',P2,Q2,'g'); 
xlabel('oś wartości rzeczywistych P1, P2'); 
ylabel('oś wartości urojonych Q1, Q2'); 
gtext('L1=1;C1=1;Rd1=1;Rl1=4'); 
gtext('L2=1;C2=0.2;Rd2=5;Rl2=1');grid;pause; 
close;clear;clc; 
 

3.5  Określenie sterowalności i obserwowalności układu 
 

Korzystając  z  równania  różniczkowego  (1)  będącego  podstawą  obliczeń  oraz  dwie  zmienne 

stanu x

1

(t) i x

2

(t) = x

1

’(t) można ułożyć następujący układ równań obiektu dynamicznego (11). 

 

 

 

 

 

 

 

 











































t

x

t

y

t

u

LC

R

t

x

C

R

L

R

t

x

LC

R

R

R

t

x

t

x

t

x

d

d

L

d

L

d

1

2

1

2

2

1

1

1

'

 

(11) 

Wielowymiarowy  liniowy  układ  dynamiczny  można  opisać  równaniami  stanu  x’(t)  =  Ax(t)  + 

Bu(t) i równaniem wyjścia y(t) = Cx(t) + Du(t). W rozpatrywanym przykładzie mamy do czynienia 
z układem o jednym sygnale wejściowym u(t) i jednym wyjściowym i(t), zatem wektory u(t) = [u(t)] 
i y(t) = [y(t)]. Wektor zmiennych stanu przyjmuje postać 

 

 















t

x

t

x

t

x

2

1

)

(

 

a macierze stanu A, wejścia B, wyjścia C i transmisyjna D postaci 









































C

R

L

R

LC

R

R

R

A

d

L

d

L

d

1

1

0

 





















LC

R

B

d

1

0

 





0

1



C

 

 

0



D

 

Układ jest całkowicie sterowalny, gdy macierz S 

S

b  Ab  A B .



2

..  = n

 

czyli jest rzędu n lub gdy wyznacznik macierzy S 

det

det

S

b  .





 b  Ab  A

..  

 0

2

. 

Układ jest całkowicie obserwowalny, gdy rząd macierzy O 





O

c





 





  c    A c     A

  ...  

2

 

jest równy n. n jest liczbą kolumn macierzy S i O. W rozpatrywanym przykładzie n = 2. 

Do  wyznaczenia  sterowalności  w  MATLAB,  stosuje  się  funkcję  rank(macierz),  która  oblicza 

rząd  macierzy,  w  naszym  przypadku:  S  lub  O.  Parametry  funkcji  rank  określa  się  na  podstawie 
macierzy typu A, B lub C. 

 
%Przykład m-pliku; 
%JCK; 
disp('Sprawdzenie sterowalności i obserwowalności'); 
disp('============================================================
======') 

 

13 

disp('Sprawdzenie  sterowalności  i  obserwowalności  dla  L1=1;  C1=1; 
Rd1=1; Rl1=4'); 
L1=1; 
C1=1; 
Rd1=1; 
Rl1=4; 
A1=[0 1;(-(Rd1+Rl1)/Rd1/L1/C1) (-Rl1/L1-1/Rd1/C1)]; 
B1=[0; (-1/Rd1/L1/C1)]; 
C1=[1 0]; 
disp('obliczenie rzędu macierzy S1 i wyznacznika S1'); 
disp('macierz stanu A1');disp(A1); 
disp('macierz wejścia B1');disp(B1); 
rzd=rank([B1 A1*B1]); 
disp('rząd macierzy S1'); disp(rzd); 
wyz=det([B1 A1*B1]); 
disp('wyznacznik macierzy S1'); disp(wyz); 
disp('obliczenie rzędu macierzy O1 i wyznacznika O1'); 
disp('macierz stanu A1');disp(A1); 
disp('macierz wyjścia C1');disp(C1); 
rzd=rank([C1' A1'*C1']); 
wyz=det([C1' A1'*C1']); 
disp('rząd macierzy O1'); disp (rzd); 
disp('wyznacznik macierzy O1'); disp(wyz); 
disp('============================================================
======') 
disp('Sprawdzenie 

sterowalnosci 

i 

obserwowalnosci 

dla 

L2=1; 

C2=0.2; Rd2=5; Rl2=1'); 
L2=1; 
C2=0.2; 
Rd2=5; 
Rl2=1; 
A2=[0 1;(-(Rd2+Rl2)/Rd2/L2/C2) (-Rl2/L2-1/Rd2/C2)]; 
B2=[0; (-1/Rd2/L2/C2)]; 
C2=[1 0]; 
disp('obliczenie rzędu macierzy S1 i wyznacznika S1'); 
disp('macierz stanu A2');disp(A2); 
disp('macierz wejścia B2');disp(B2); 
rzd=rank([B2 A2*B2]); 
disp('rząd macierzy S2'); disp(rzd); 
wyz=det([B2 A2*B2]); 
disp('wyznacznik macierzy S2'); disp(wyz); 
disp('obliczenie rzędu macierzy O2 i wyznacznika O2'); 
disp('macierz stanu A2');disp(A2); 
disp('macierz wyjścia C2');disp(C2); 
rzd=rank([C2' A2'*C2']); 
wyz=det([C2' A2'*C2']); 
disp('rząd macierzy O2'); disp (rzd); 
disp('wyznacznik macierzy O2'); disp(wyz); 
disp('Nacisnij dowolny klawisz by zakonczyc.'); 
pause;clear;clc;

 

 

Po uruchomieniu powyższego m-pliku uzyskujemy następujące wyniki:

 

 
Sprawdzenie sterowalności i obserwowalności 
================================================================== 
Sprawdzenie sterowalności i obserwowalności dla L1=1; C1=1; Rd1=1; 
Rl1=4 

 

14 

obliczenie rzędu macierzy S1 i wyznacznika S1 
macierz stanu A1 
     0     1 
    -5    -5 
 
macierz wejścia B1 
     0 
    -1 
 
rząd macierzy S1 
     2 
 
wyznacznik macierzy S1 
    -1 
 
obliczenie rzędu macierzy O1 i wyznacznika O1 
macierz stanu A1 
     0     1 
    -5    -5 
 
macierz wyjścia C1 
     1     0 
 
rząd macierzy O1 
     2 
 
wyznacznik macierzy O1 
     1 
 
================================================================== 
Sprawdzenie  sterowalnosci  i  obserwowalnosci  dla  L2=1;  C2=0.2; 
Rd2=5; Rl2=1 
obliczenie rzędu macierzy S1 i wyznacznika S1 
macierz stanu A2 
         0    1.0000 
   -6.0000   -2.0000 
 
macierz wejścia B2 
     0 
    -1 
 
rząd macierzy S2 
     2 
 
wyznacznik macierzy S2 
    -1 
 
obliczenie rzędu macierzy O2 i wyznacznika O2 
macierz stanu A2 
         0    1.0000 
   -6.0000   -2.0000 
 
macierz wyjścia C2 
     1     0 
 
rząd macierzy O2 
     2 

 

15 

 
wyznacznik macierzy O2 
     1 
 
Nacisnij dowolny klawisz by zakonczyc. 
 

3.6  Wyznaczenie współczynników i pierwiastków wielomianu charakterystycznego układu.  

 
Współczynniki  i  pierwiastki  wielomianu  charakterystycznego  wyznacza  się  na  podstawie 

macierzy stanu A za pomocą standardowej funkcji poly. 
poly(A) – dla macierzy kwadratowej A stopnia n tworzy n+1 elementowy wektor współczynników 
wielomianu  charakterystycznego  det(sI-A).  Współczynniki  uporządkowane  są  w  kierunku 
malejących potęg. 
Pierwiastki wielomianu charakterystycznego zostają wyznaczone przy pomocy funkcji roots (W). 
 

%Przykład m-pliku; 
%JCK; 
%Wyznaczenie 

współczynników 

i 

pierwiastków 

wielomianu 

%charakterystycznego; 
disp('Wyznaczenie 

współczynników 

i 

pierwiastków 

wielomianu 

charakterystycznego'); 
disp('************************************************************
************'); 
disp('dla A1[0 1;(-(Rd1+Rl1)/Rd1/L1/C1) (–Rl1/L-1/Rd1/C1)]'); 
L1=1; 
C1=1; 
Rd1=1; 
Rl1=4; 
L2=1; 
C2=0.2; 
Rd2=5; 
Rl2=1; 
a1=-(Rd1+Rl1)/(Rd1*L1*C1); 
b1=-(Rl1/L1)-1/(Rd1*C1); 
A1=[0 1;a1 b1]; 
disp('wspolczynniki=');disp(poly(A1)); 
W1=poly(A1) 
disp('pierwiastki=');disp(roots(W1)); 
disp('dla A2[0 1;(-(Rd2+Rl2)/Rd2/L2/C2) (–Rl2/L2-1/Rd2/C2)]'); 
a2=-(Rd2+Rl2)/(Rd2*L2*C2); 
b2=-(Rl2/L2)-1/(Rd2*C2); 
A2=[0 1;a2 b2]; 
disp('wspolczynniki=');disp(poly(A2)); 
W2=poly(A2); 
disp('pierwiastki=');disp(roots(W2)); 
disp('Nacisnij dowolny klawisz by zakonczyc.'); 
pause;clear;clc;

 

 
Po uruchomieniu powyższego m-pliku uzyskujemy następujące wyniki: 

 
Wyznaczenie 

współczynników 

i 

pierwiastków 

wielomianu 

charakterystycznego 
******************************************************************
****** 
dla A1[0 1;(-(Rd1+Rl1)/Rd1/L1/C1) (–Rl1/L-1/Rd1/C1)] 
wspolczynniki= 

 

16 

     1     5     5 
 
W1 = 
 
     1     5     5 
 
pierwiastki= 
   -3.6180 
   -1.3820 
 
dla A2[0 1;(-(Rd2+Rl2)/Rd2/L2/C2) (–Rl2/L2-1/Rd2/C2)] 
wspolczynniki= 
    1.0000    2.0000    6.0000 
 
pierwiastki= 
  -1.0000 + 2.2361i 
  -1.0000 - 2.2361i

 

 
Nacisnij dowolny klawisz by zakonczyc. 

 
3.7  Polecenia i pytania do plików skryptowych. 
 

 

Nazwa pliku 

Punkt 

instrukcji 

Działanie 

programu 

Zadania do wykonania w czasie zajęć 

1. 

P3_1_char_skokow
e.m 

3.1 

Wyznaczenie 
charakterystyk 
skokowych 
(czasowych) 
członu inercyjnego 
drugiego rzędu i 
stabilnego członu 
oscylacynego. 

  Wydrukować ch-ki odpowiedzi skokowych. 
  Uzasadnić na podstawie wykreślonych 

charakterystyk, jaki wpływ na ch-kę skokową mają 
wartości parametrów obiektu inercyjnego drugiego 
rzędu: T

1

, T

2

- stałe czasowe i k - współczynnik 

proporcjonalności (wzmocnienia). Jaki jest związek 
pomiędzy tymi parametrami a parametrami 
elektrycznymi charakteryzującymi rozpatrywany 
obwód? 

  Uzasadnić na podstawie wykreślonych 

charakterystyk, jaki wpływ na ch-kę skokową mają 
wartości parametrów obiektu oscylacyjnego: T

0

- 

stała czasowa, ξ

0

- współczynnik tłumienia i k - 

współczynnik proporcjonalności (wzmocnienia). 
Jaki jest związek pomiędzy tymi parametrami a 
parametrami elektrycznymi charakteryzującymi 
rozpatrywany obwód? 

  Podać interpretacje stałych czasowych obiektu 

inercyjnego drugiego rzędu. 

  Czy na podstawie uzyskanych ch-k skokowych 

można ocenić liniowość względnie nieliniowość 
układu dynamicznego? 

2. 

P3_2_char_ampfaz.
m 

3.2 

Wyznaczenie 
charakterystyk 
amplitudowo-
fazowych członu 
inercyjnego 
drugiego rzędu i 
stabilnego członu 
oscylacynego. 

  Wydrukować ch-ki amplitudowo-fazowe. 
  Od których parametrów i jak zależy 

charakterystyka amplitudowo-fazowa? Zaznaczyć 
strzałką na ch-kach kierunek wzrostu pulsacji. 

  Dla wybranej ch-ki amplitudowo-fazowej dla 

wartości pulsacji ω=0 i ω=10 wrysować wektory 
modułu transmitancji i odczytać z wykresu dla tych 
modułów wartości części rzeczywistych (P) i 

 

17 

 

Nazwa pliku 

Punkt 

instrukcji 

Działanie 

programu 

Zadania do wykonania w czasie zajęć 

urojonych (Q), wskazać kąt przesunięcia fazowego i 
obliczyć wartość tgφ tego przesunięcia. 

  Naszkicować lub wydrukować ch-ki amplitudowo-

fazowe Nyquista. 

  Czym różnią się ch-ki Nyquista od zwykłych 

charakterystyk amplitudowo-fazowych? 

  Porównać ch-kę Nyquista obiektu inercyjnego 

drugiego rzędu z ch-ką Nyquista stabilnego członu 
oscylacyjnego. 

  W jakim celu stosuje się ch-ki Nyquista? 
  Co można powiedzieć o stabilności układu 

regulacji na podstawie uzyskanych ch-k Nyquista? 

3. 

P3_3_char_amp_i_f
az.m 

3.3 

Wyznaczenie 
charakterystyk 
amplitudowych w 
funkcji pulsacji i 
logarytmicznych 
charakterystyk 
amplitudowo-
fazowych (tzw. 
(Bode Diagrams) 
członu inercyjnego 
drugiego rzędu i 
stabilnego członu 
oscylacynego. 

  Uzyskać odpowiednio 2 typy wykresów. 
  Wydrukować ch-ki amplitudowe i charakterystyki 

logarytmiczne amplitudowo-fazowe. 

  Określić różnice między charakterystykami 

amplitudowymi i logarytmicznymi amplitudowymi. 

  Jaki wpływ na przebieg charakterystyk mają 

wartości parametrów obiektu inercyjnego? Jaki 
wpływ mają parametry elektryczne obwodu? 

  Odczytać z wykresów ch-k amplitudowych 

wartości ekstremów części urojonych i 
odpowiadające im wartości pulsacji ω. W którym 
kierunku przesuwają się wartości ekstremów 
urojonych ze wzrostem pulsacji? 

  Na wykresach ch-k logarytmicznych 

zidentyfikować i opisać ch-ki fazowe. 

  Odczytać w [dB] „podbicie” ch-k amplitudowych 

a następnie wartości pulsacji dla których występują 
te maksima, obliczyć jakie to częstotliwości w [Hz]. 

4. 

P3_4_stab_Michajlo
w.m 

3.4 

Wyznaczenie 
stabiln ości 
graficznym 
kryterium 
Michajłowa członu 
inercyjnego 
drugiego rzędu i 
stabilnego członu 
oscylacynego. 

  Czy człony, dla parametrów dla których zostały 

wykreślone ch-ki, są obiekteami stabilnym? 

  Jakie są inne kryteria stwierdzania stabilności 

układów automatyki? 

5. 

P3_5_ster_obs.m 

3.5 

Określenie 
sterowalności i

 

obserwowalności 
członu inercyjnego 
drugiego rzędu i 
stabilnego członu 
oscylacynego. 

  Uzyskać obliczenia rzędu macierzy i wartości 

wyznaczników, tak, jak w przykładzie. 

  Dokonać interpretacji otrzymanych wyników 

oceniając sterowalność i obserwowalność. 

6. 

P3_6_wsp_wielom_ch
ar.m 

Przykład 

6.3 

Wyznaczenie 
współczynników i 
pierwiastków 
wielomianu 
charakterystyczneg

  Uzyskać obliczenia współczynników i 

pierwiastków wielomianu charakterystycznego - jak 
w przykładzie. 

  Napisać wielomiany charakterystyczne i porównać 

z wielomianem mianownika transmitancji układu. 

 

18 

 

Nazwa pliku 

Punkt 

instrukcji 

Działanie 

programu 

Zadania do wykonania w czasie zajęć 

o członu 
inercyjnego 
drugiego rzędu i 
stabilnego członu 
oscylacynego.

 

  Dokonać interpretacji otrzymanych wyników. 

 
1.1  Wnioski.  
 
W  sprawozdaniu  należy  udzielić  odpowiedzi  na  wszystkie  pytania  zawarte  w  instrukcji  oraz 
zalecone  przez  prowadzącego.  Sprawozdanie  powinno  zawierać  opisane  oryginalne  wydruki 
wyników  uzyskanych  podczas  realizacji  ćwiczenia.  We  wnioskach  należy  uzupełnić 
i uporządkować  odpowiedzi  na  postawione  pytania  lub  zagadnienia  należące  do  zadań 
wykonanych w czasie zajęć. 
 
2 

Niezbędny zakres tematyczny obowiązujący do przygotowania się do zajęć: 

 
Pojęcia:  obiekt  dynamiczny  ciągły  i  liniowy,  równanie  różniczkowe  opisujące  taki  obiekt, 
identyfikacja  obiektu  dynamicznego,  transmitancja  operatorowa  i  widmowa,  moduł  transmitancji, 
kąt  przesunięcia  fazowego,  podstawowe  człony dynamiczne (układy  – proporcjonalny, inercyjny, 
różniczkujący  rzeczywisty,  całkujący,  inercyjny  drugiego  rzędu  i  oscylacyjny),  transmitancje  tych 
członów  i  podstawowe  charakterystyki  członów,  obserwowalność, sterowalność, stabilność (oraz 
kryteria  i  wyznaczanie  obserwowalności,  sterowalności  i  stabilności)  układu  dynamicznego  i 
układu  regulacji,  równanie  charakterystyczne  i  wielomian  charakterystyczny,  stała  czasowa 
w podstawowych członach dynamicznych, 
Charakterystyki  układów  dynamicznych  –  charakterystyki  skokowe  i  częstotliwościowe, 
interpretacja 

charakterystyk, 

sposoby 

wyznaczania 

charakterystyk, 

obliczanie 

modułu 

transmitancji i przesunięcia fazowego, 
Równania  wektorowo-macierzowe  układu  dynamicznego.  Macierze:  stanu,  wejścia  i  wyjścia. 
Zastosowanie macierzy równań stanu. 
 
3 

Literatura  pomocna  do  wykonania  ćwiczenia,  przygotowania  sprawozdania  i  przygotowania 
się do zaliczenia: 

 
1.  Brzózka J., Ćwiczenia z automatyki w MATLABIE i SIMULINKU, Mikom, Warszawa, 1997. 
2.  Dębowski A.: Automatyka. Podstawy teorii, WNT, Warszawa 2008. 
3.  Frelek B. i inni, Laboratorium podstaw automatyki, W PW, 1984. 
4.  Mrozek B+Z, Matlab 5.x i Simulink 2.x, Wydawnictwo PLJ, Warszawa 1998. 
5.  Nowacki  P.,  Szklarski  L.,  Górecki  H.:  Podstawy  teorii  układu  regulacji  autonomicznej.  Tom  I, 

PWN, Warszawa 1970. 

6.  Pełczewski W., Teoria sterowania, WNT, Warszawa 1980. 
7.  Żelazny M., Podstawy automatyki, PWN, Warszawa 1976. 
8.  Badania  obiektów  dynamicznych  (z  zastosowaniem  programu  MATLAB)  -  Przykłady 

programów  do  wykonania  w  czasie  zajęć  laboratoryjnych.  Opracowanie  WTPW,  Warszawa 
2002.