Ostatnia aktualizacja: 2010-01-04



M. Tomera

1

Akademia Morska w Gdyni
Katedra Automatyki Okrętowej

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Mirosław Tomera

Można zaprojektować układ sterowania ciągłego i zaimplementować go w układach sterowania
cyfrowego stosując metody aproksymacji równań różniczkowych. Układ taki będzie mógł spełniać
zadane wymagania jakościowe jeśli częstotliwość próbkowania będzie przynajmniej 30 razy większa
od szerokości pasma układu. Na rysunku 1(a) pokazany jest typowy układ ciągły. Obliczenie sygnału
uchybu e i kompensacja dynamiczna C(s) może zostać zrealizowana w komputerze tak jak pokazano
to na rysunku 1(b). Podstawowe różnice pomiędzy tymi dwoma implementacjami są takie, że układ
cyfrowy przetwarza próbki pomierzonego sygnału wyjściowego, a nie sygnał ciągły i dynamika
opisana przez C(s) jest implementowana przez równania algebraiczne, które nazywają się
równaniami różnicowymi.

C(s)

u(t)

G(s)

H(s)

y(t)

e(t)

y(t)

r(t)

Regulator ciągły

Obiekt

Czujnik

(a)

Zegar

C/A

i hold

u(kT)

u(t)

Równania
różnicowe

C/A

y(t)

e(kT)

y(kT)

r(t)

T

r(kT)

H(s)

Czujnik

G(s)

Obiekt

T

Regulator cyfrowy

Impulsator

y(t)

(b)

Rys. 1. Podstawowy schemat blokowy układu sterowania, (a) układ ciągły, (b) z systemem mikropcesorowym.

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

2

Pewnym szczególnym sposobem zrealizowania aproksymaty dla komputera cyfrowego w celu
rozwiązania równania różniczkowego jest metoda Eulera. Metoda ta wyprowadzona została
z następującej definicji różniczki

t

x

x

t

δ

δ

δ

0

.

lim

→

=

(1)

gdzie x

δ

jest zmianą zmiennej x w czasie t

δ

, nawet jeśli t

δ

nie jest całkiem równe zero to ta

zależność może być prawdziwa po zastosowaniu poniższych aproksymat. W metodzie Eulera
wyróżnia się dwie metody:

•

aproksymację prostokątną w przód (forward rectangular rule)

( )

T

k

x

k

x

k

x

)

(

)

1

(

.

−

+

≅

(2)

gdzie:

k

k

t

t

T

−

=

+

1

(przedział próbkowania w sekundach)

kT

t

k

=

(dla stałego przedziału próbkowania)

k liczba całkowita
x(k) wartość x w chwili t

k

x(k+1) wartość x w chwili t

k+1

•

aproksymację prostokątną wstecz (backward rectangular rule)

( )

T

k

x

k

x

k

x

)

1

(

)

(

.

−

−

≅

(3)

gdzie:

1

−

−

=

k

k

t

t

T

(przedział próbkowania w sekundach)

kT

t

k

=

(dla stałego przedziału próbkowania)

k liczba całkowita
x(k) wartość x w chwili t

k

x(k

−

1) wartość x w chwili t

k

−

1

Aproksymacje te mogą być używane w miejsce wszystkich różniczek, które pojawiają się
w równaniach różniczkowych regulatora po to aby uzyskać zbiór równań, które mogą być rozwiązane
przez komputer cyfrowy. Równania te nazywane są równaniami różnicowymi i są rozwiązywane
cyklicznie z krokiem czasu o długości T. Przykład pierwszy ilustruje sposób uzyskiwania równań
różnicowych na podstawie równań różniczkowych metodą dyskretyzacji Eulera aproksymacją wprzód.

Przykład 1

Korzystając z metody Eulera w przód, dokonaj dyskretyzacji z okresem próbkowania T = 0.1 [s]
następującego równania różniczkowego

2

2

)

(

dt

t

y

d

+ 2

dt

t

dy )

(

+ 5

)

(

t

y

= 10

)

(

1

t

⋅

(1.1)

z warunkami początkowymi

)

0

(

y

=

−

4

(1.2)

)

0

(

)

1

(

y

=

−

10

(1.3)

a) Na podstawie uzyskanego równania różnicowego wygeneruj w Matlabie kolejne próbki

sygnału dyskretnego.

b) Zaimplementuj w

Simulinku powyższe równanie różniczkowe wraz z warunkami

początkowymi.

c) Dokonaj porównania na wykresie wyników uzyskanych w punkcie (a) i (b).

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

3

Rozwiązanie: Chcąc uzyskać postać dyskretną równania różniczkowego (1.1) w pierwszej
kolejności należy zastąpić pochodne odpowiednimi różnicami czyli równanie różniczkowe
przechodzi w następujące równanie różnicowe

2

2

)

(

T

kT

y

∆

+ 2

T

kT

y

)

(

∆

+ 5

)

(

kT

y

= 10

)

(

1

kT

⋅

(1.4)

Przy użyciu metody prostokątnej Eulera w przód pierwsza różnica zapisywana jest następująco:

)

(

kT

y

∆

=

[

]

T

k

y

)

1

(

+

−

( )

kT

y

(1.5)

natomiast druga różnica

)

(

2

kT

y

∆

=

[

]

( )

kT

y

T

k

y

∆

−

+

∆

)

1

(

=

[

]

T

k

y

)

2

(

+

−

[

]

T

k

y

)

1

(

2

+

⋅

+

( )

kT

y

(1.6)

Podstawiając różnice opisane wzorami (1.5) oraz (1.6) do równania różnicowego (1.4)
i przekształcając je uzyskuje się następującą postać nadającą się do implementacji
komputerowej

[

]

T

k

y

)

2

(

+

=

(

) (

)

[

]

T

k

y

T

1

2

2

+

−

+

(

)

( )

kT

y

T

T

⋅

⋅

−

⋅

+

−

2

5

2

1

+

)

(

1

10

2

kT

T

⋅

⋅

(1.7)

Po podstawieniu do wzoru (1.7) okresu próbkowania T = 0.1 [s] uzyskuję się następujące
równanie różnicowe

)

2

(

+

k

y

−

)

1

(

8

.

1

+

⋅

k

y

+

)

(

85

.

0

k

y

⋅

=

)

(

1

1

.

0

k

⋅

(1.8)

Pozostaje do przeliczenia jeszcze drugi warunek początkowy

)

0

(

)

1

(

y

=

T

y

y

)

0

(

)

1

(

−

(1.9)

Stąd wyznaczona zostanie wartość dla k = 1.

)

1

(

y

= y(0) +

)

0

(

)

1

(

y

T

⋅

(1.10)

Równanie różnicowe (1.7) wraz z warunkami początkowymi (1.8) oraz (1.9) można rozwiązać
bezpośrednio w Matlabie przy użyciu odpowiednio zapisanego programu.

Aby zaimplementować w Simulinku równanie różniczkowe (1.1) wraz z warunkami

początkowymi to najpierw należy to równanie zapisać w postaci zmiennych dynamicznych

2

1

.

x

x

=

4

)

0

(

1

−

=

x

)

(

1

10

2

5

2

1

2

.

t

x

x

x

⋅

+

−

−

=

10

)

0

(

2

−

=

x

(1.11)

1

x

y

=

Równania (1.10) można zamodelować w postaci następującego modelu Simulinka (Rys. 1.1),
który zachowany został pod nazwą

ex_1_sim.mdl

.

Rys. 1.1. Schemat pozwalający na rozwiązanie w Simulinku równań dynamicznych (1.11).

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

4

Na rysunku 1.2 znajdują się wyniki rozwiązania równania różniczkowego i uzyskanego
równania różnicowego dla T = 0.1 [s].

0

1

2

3

4

5

6

-6

-4

-2

0

2

4

6

Przybliżone rozwiązanie równania różniczkowego

t [s]

y(

k)

Rys. 1.2. Uzyskane wyniki rozwiązania równania różniczkowego (linia ciągłą) i równania różnicowego

(próbki) uzyskanego z dyskretyzacji dla okresu próbkowania T = 0.1 [s].

Wyniki przedstawione na rysunku 1.2. uzyskane zostały przy użyciu następującego kodu
programu zapisanego w Matlabie.

close all

% Zamkni

ę

cie wszystkich okien graficznych

clear

% Wyczyszczenie pami

ę

ci roboczej Matlaba

tmax = 6;

% Odcinek czasu

open_system('ex_1_sim')

% Otwarcie przygotowanego modelu Simulinka

sim('ex_1_sim',tmax)

% Wykonanie symulacji w zadanym odcinku
% czasu

close_system

% Zamkni

ę

cie modelu Simulinka

tC = tout;

% Podstawienie wektora czasu

yC = wyniki_sym(:,2);

% Pobranie z przestrzeni roboczej Matlaba
% uzyskanego rozwi

ą

zania

% Rozwi

ą

zanie dyskretne

Tp = 0.1;

% Okres próbkowania

y0 = -4;

% Warunki pocz

ą

tkowe

Dy0 = -10;
y1 = y0 + Dy0*Tp;

% Przeliczenie warto

ś

ci pocz

ą

tkowych

% na warto

ść

drugiej próbki

tk(1) = 0*Tp;

% Współrz

ę

dne pierwszych dwóch próbek

yk(1) = y0;
tk(2) = 1*Tp;
yk(2) = y1;
for k = 1:(tmax/Tp)-1;

tk(k+2) = (k+1)*Tp;
yk(k+2) = (2-2*Tp)*yk(k+1) - (1-2*Tp+5*Tp^2)*yk(k) + 10*Tp^2;

end;
% Wykre

ś

lenie uzyskanych wyników w Matlabie

figure(1)

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

5

plot( tC, yC, 'k-', tk, yk, 'ko')
title('Rozwi

ą

zanie równania ró

ż

nicowego')

xlabel('t [s]')
ylabel('y(k)')
grid on

Przykład 2

Korzystając z metody prostokątów Eulera z reguły wstecznej, dokonaj dyskretyzacji równania
różniczkowego (1.1) wraz z warunkami początkowymi (1.2) i (1.3) znajdującego się
w przykładzie 1. Okres próbkowania T = 0.1 [s]. Na wykresie dokonaj porównania wyników
uzyskanych obydwoma metodami.

Rozwiązanie: Chcąc uzyskać postać dyskretną równania różniczkowego (1.1) metodą
prostokątów Eulera wstecz postępuje się w podobny sposób jak w przykładzie 1. Do równania
(1.5) dokonuje się następujących podstawień: pierwsza różnica zapisywana jest następująco:

)

(kT

y

∆

=

( )

kT

y

−

[

]

T

k

y

)

1

(

−

(2.1)

natomiast druga różnica

)

(

2

kT

y

∆

=

( )

[

]

T

k

y

kT

y

)

1

(

−

∆

−

∆

=

( )

kT

y

−

[

]

T

k

y

)

1

(

2

−

⋅

+

[

]

T

k

y

)

2

(

−

(2.2)

Podstawiając różnice opisane wzorami (2.1) oraz (2.2) do równania różnicowego (1.4)
i przekształcając je uzyskuje się w ten sposób kolejną postać równania nadającą się do
implementacji komputerowej

)

(kT

y

=

(

)

[

]

T

k

y

T

T

T

1

5

2

1

2

2

2

−

+

+

+

−

(

)

[

]

T

k

y

T

T

2

5

2

1

1

2

−

+

+

+

)

(

1

5

2

1

10

2

2

kT

T

T

T

⋅

+

+

(2.3)

Po podstawieniu do wzoru (2.3) okresu próbkowania T = 0.1 [s] uzyskuje się następujące
równanie różnicowe

0

1

2

3

4

5

6

-6

-4

-2

0

2

4

6

t [s]

y(

k)

m. Eulera w przód

m. Eulera wstecz

Przybliżone rozwiązania równania różniczkowego

Rys. 2.1. Porównanie rozwiązań równań różnicowych uzyskanych z aproksymacji równania

różniczkowego (linia ciągłą) metodą prostokątów Eulera przy użyciu reguł w przód i wstecz
dla okresu próbkowania T = 0.1 [s].

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

6

)

(k

y

−

)

1

(

76

.

1

−

⋅

k

y

+

)

2

(

8

.

0

−

⋅

k

y

=

)

(

1

08

.

0

k

⋅

(2.4)

Warunki początkowe przelicza się w identyczny sposób jak w metodzie Eulera w przód.
Z porównania uzyskanych równań (1.8) oraz (2.4) widać, że różnią się między sobą
wartościami współczynników. Na rysunku 2.1 dokonane zostało porównanie uzyskanych
wyników.

W przykładzie 2 przedstawione zostało porównanie aproksymacji rozwiązań równania różniczkowego
przy użyciu metody prostokątnej Eulera. Z rysunku 2.1 wynika, że jeśli do aproksymacji zastosuje się
regułę Eulera w przód wówczas uzyskuje się rozwiązania z nadmiarem, natomiast przy zastosowaniu
reguły Eulara wstecz rozwiązania z niedomiarem.

ĆWICZENIA W MATLABIE

M.1. Poniższe równania różniczkowe

zaimplementuj w postaci schematów
Simulinka, następnie korzystając
z

metody prostokątów Eulera w

przód

(forward rectangular rule) wyznacz
aproksymujące równania różnicowe. Okres
próbkowania T = 0.1 [s]. Sekwencyjnie
przy użyciu Matlaba wykreśl rozwiązania
równań różnicowych i porównaj je
z rozwiązaniami uzyskanymi w Simulinku.

a)

3

3

)

(

dt

t

y

d

+ 5

2

2

)

(

dt

t

y

d

+ 4

dt

t

dy )

(

= 2

)

(t

δ

)

0

(

y

= 1

)

0

(

)

1

(

y

=

−

1

)

0

(

)

2

(

y

= 1

b)

2

2

)

(

dt

t

y

d

+ 5

dt

t

dy )

(

+ 6

)

(t

y

= 2

t

e

−

)

0

(

y

= 0

)

0

(

)

1

(

y

=

−

2

c)

dt

t

dy )

(

+ 2

)

(t

y

= 6

t

2

cos

)

0

(

y

=

−

2

d)

3

3

)

(

dt

t

y

d

+ 4

2

2

)

(

dt

t

y

d

+ 3

dt

t

dy )

(

= 3

)

(

t

δ

)

0

(

y

=

−

1

)

0

(

)

1

(

y

=

−

1

)

0

(

)

2

(

y

= 1

e)

dt

t

dy )

(

+ 2

)

(

t

y

=

2

2

1

t

)

0

(

y

=

−

1

f)

dt

t

dy )

(

+ 3

)

(

t

y

= 4

t

4

sin

)

0

(

y

= 1

g)

2

2

)

(

dt

t

y

d

+ 4

dt

t

dy )

(

+ 5

)

(

t

y

=

)

(

1

t

)

0

(

y

=

−

3

)

0

(

)

1

(

y

= 2

h)

2

2

)

(

dt

t

y

d

+ 3

dt

t

dy )

(

+ 2

)

(

t

y

= 3

t

e

3

−

)

0

(

y

= 0

)

0

(

)

1

(

y

= 1

i)

2

2

)

(

dt

t

y

d

+ 2

dt

t

dy )

(

+ 2

)

(

t

y

= sin2t

( )

0

y

= 2

( )

0

)

1

(

y

=

−

3

j)

2

2

)

(

dt

t

y

d

+ 2

dt

t

dy )

(

+

)

(

t

y

=

)

(

1

5

t

⋅

( )

0

y

= 1

( )

0

)

1

(

y

=

−

2

k)

2

2

)

(

dt

t

y

d

+

dt

t

dy )

(

+ 4

)

(

t

y

= cos3t

( )

0

y

= 1

( )

0

)

1

(

y

= 2

l)

2

2

)

(

dt

t

y

d

+ 3

dt

t

dy )

(

+ 6

)

(

t

y

=

t

e

t

3

sin

2

−

( )

0

y

= 0

( )

0

)

1

(

y

= 3

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

7

M.2. P

rzy użyciu metody prostokątów Eulera

wstecz (backward rectangular rule)
wyznacz aproksymujące równania
różnicowe dla równań różniczkowych
znajdujących się w zadaniu M.1. Okres
próbkowania T = 0.1 [s]. Sekwencyjnie
przy użyciu Matlaba wykreśl rozwiązania
i porównaj je z rozwiązaniami uzyskanymi
w zadaniu M.1.

ODPOWIEDZI DO WYBRANYCH ĆWICZEŃ

M1.

a)

)

(

02

.

0

)

(

54

.

0

)

1

(

04

.

2

)

2

(

5

.

2

)

3

(

k

k

y

k

y

k

y

k

y

δ

⋅

=

−

+

+

+

−

+

1

)

0

(

=

y

9

.

0

)

1

(

=

y

81

.

0

)

2

(

=

y

b)

(

)

k

k

y

k

y

k

y

9048

.

0

02

.

0

)

(

56

.

0

)

1

(

5

.

1

)

2

(

⋅

=

+

+

−

+

0

)

0

(

=

y

2

.

0

)

1

(

−

=

y

c)

(

)

k

k

y

k

y

⋅

⋅

=

−

+

2

.

0

cos

6

.

0

)

(

8

.

0

)

1

(

2

)

0

(

−

=

y

d)

)

(

03

.

0

)

(

63

.

0

)

1

(

23

.

2

)

2

(

6

.

2

)

3

(

k

k

y

k

y

k

y

k

y

δ

⋅

=

−

+

+

+

−

+

1

)

0

(

−

=

y

1

.

1

)

1

(

−

=

y

19

.

1

)

2

(

−

=

y

e)

2

0005

.

0

)

(

8

.

0

)

1

(

k

k

y

k

y

⋅

=

−

+

1

)

0

(

−

=

y

f)

(

)

k

k

y

k

y

⋅

⋅

=

−

+

4

.

0

sin

4

.

0

)

(

7

.

0

)

1

(

1

)

0

(

=

y

g)

)

(

1

01

.

0

)

(

65

.

0

)

1

(

6

.

1

)

2

(

k

k

y

k

y

k

y

⋅

=

+

+

−

+

3

)

0

(

−

=

y

8

.

2

)

1

(

−

=

y

h)

(

)

k

k

y

k

y

k

y

7408

.

0

03

.

0

)

(

72

.

0

)

1

(

7

.

1

)

2

(

⋅

=

+

+

−

+

0

)

0

(

=

y

1

.

0

)

1

(

=

y

i)

(

)

k

k

y

k

y

k

y

⋅

⋅

=

+

+

−

+

2

.

0

sin

01

.

0

)

(

82

.

0

)

1

(

8

.

1

)

2

(

2

)

0

(

=

y

7

.

1

)

1

(

=

y

j)

)

(

1

05

.

0

)

(

81

.

0

)

1

(

8

.

1

)

2

(

k

k

y

k

y

k

y

⋅

=

+

+

−

+

1

)

0

(

=

y

8

.

0

)

1

(

=

y

k)

(

)

k

k

y

k

y

k

y

⋅

⋅

=

+

+

−

+

3

.

0

cos

01

.

0

)

(

94

.

0

)

1

(

9

.

1

)

2

(

1

)

0

(

=

y

2

.

1

)

1

(

=

y

l)

(

)

(

)

k

k

y

k

y

k

y

k

⋅

⋅

=

+

+

−

+

3

.

0

cos

8187

.

0

01

.

0

)

(

76

.

0

)

1

(

7

.

1

)

2

(

0

)

0

(

=

y

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

8

3

.

0

)

1

(

=

y

M2.

a)

)

(

013

.

0

)

3

(

6494

.

0

)

2

(

2727

.

2

)

1

(

6234

.

2

)

(

k

k

y

k

y

k

y

k

y

δ

⋅

=

−

−

−

+

−

−

1

)

0

(

=

y

9

.

0

)

1

(

=

y

81

.

0

)

2

(

=

y

b)

(

)

k

k

y

k

y

k

y

9048

.

0

0128

.

0

)

2

(

6410

.

0

)

1

(

6026

.

1

)

(

⋅

=

−

+

−

−

0

)

0

(

=

y

2

.

0

)

1

(

−

=

y

c)

(

)

k

k

y

k

y

⋅

⋅

=

−

−

2

.

0

cos

5

.

0

)

1

(

8333

.

0

)

(

2

)

0

(

−

=

y

d)

)

(

021

.

0

)

3

(

6993

.

0

)

2

(

3776

.

2

)

1

(

6783

.

2

)

(

k

k

y

k

y

k

y

k

y

δ

⋅

=

−

−

−

+

−

−

1

)

0

(

−

=

y

1

.

1

)

1

(

−

=

y

19

.

1

)

2

(

−

=

y

e)

2

00042

.

0

)

1

(

8333

0

)

(

k

k

y

.

k

y

⋅

=

−

−

1

)

0

(

−

=

y

f)

(

)

k

k

y

k

y

⋅

⋅

=

−

4

.

0

sin

3077

.

0

)

(

7692

.

0

)

(

1

)

0

(

=

y

g)

)

(

1

0069

.

0

)

2

(

6897

.

0

)

1

(

6552

.

1

)

(

k

k

y

k

y

k

y

⋅

=

−

+

−

−

3

)

0

(

−

=

y

8

.

2

)

1

(

−

=

y

h)

(

)

k

k

y

k

y

k

y

7408

.

0

0227

.

0

)

2

(

7576

.

0

)

1

(

7424

.

1

)

(

⋅

=

−

+

−

−

0

)

0

(

=

y

1

.

0

)

1

(

=

y

i)

(

)

k

k

y

k

y

k

y

⋅

⋅

=

+

−

−

2

.

0

sin

0082

.

0

)

(

8197

.

0

)

1

(

8033

.

1

)

(

2

)

0

(

=

y

7

.

1

)

1

(

=

y

j)

)

(

1

0413

.

0

)

2

(

8264

.

0

)

1

(

8182

.

1

)

(

k

k

y

k

y

k

y

⋅

=

−

+

−

−

1

)

0

(

=

y

8

.

0

)

1

(

=

y

k)

(

)

k

k

y

k

y

k

y

⋅

⋅

=

−

+

−

−

3

.

0

cos

0088

.

0

)

2

(

8772

.

0

)

1

(

8421

.

1

)

(

1

)

0

(

=

y

2

.

1

)

1

(

=

y

l)

(

)

(

)

k

k

y

k

y

k

y

k

⋅

⋅

=

−

+

−

−

3

.

0

cos

8187

.

0

0074

.

0

)

2

(

7353

.

0

)

1

(

6912

.

1

)

(

0

)

0

(

=

y

3

.

0

)

1

(

=

y

LITERATURA

1. Amborski K., Teoria sterowania. Podręcznik programowany, PWN, Warszawa, 1985.
2. Franklin G.F, Powell J.D., Emami-Naeini A. Digital Control of Dynamic Systems, 3rd ed.

Addison-Wesley Publishing Company, 1998.

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

9

% Rozwi

ą

zanie równania ró

ż

niczkowego

%
% y''(t) + 2*y'(t) + 5*y(t) = 10*1(t)
%
% ró

ż

nymi metodami

%
% Copyright (c) by Mirosław Tomera
% Date: 2007/12/11

clear
close
clc

tmax = 6;

y0 = -4;

% Warunki pocz

ą

tkowe

Dy0 = -10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rachunek operatorowy Laplace'a
numC = [-4 -18 10]
denC = [1 2 5 0]

[rC, pC, kC] = residue( numC, denC)

M = abs( rC(1))
fi = angle( rC(1))*180/pi
A = 2*M
sigma = real( pC(1))
w = imag( pC(1))

tC = [0:0.1:tmax]';
yC = A*exp(sigma*tC).*cos(w*tC + fi*pi/180) + rC(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rozwi

ą

zania przybli

ż

one - metodami Eulera

Tp = 0.1;

% Okres próbkowania

% Rozwi

ą

zanie dyskretne metod

ą

Eulera w przód

y0a = y0;

% Warunki pocz

ą

tkowe

Dy0a = Dy0;
y1a = y0a + Dy0a*Tp;

% Przeliczenie warto

ś

ci pocz

ą

tkowych

% na warto

ść

drugiej próbki

tka(1) = 0*Tp;

% Współrz

ę

dne pierwszych dwóch próbek

yka(1) = y0a;
tka(2) = 1*Tp;
yka(2) = y1a;
for k = 1:(tmax/Tp)-1;

tka(k+2) = (k+1)*Tp;
yka(k+2) = (2-2*Tp)*yka(k+1) - (1-2*Tp+5*Tp^2)*yka(k) + 10*Tp^2;

end;

% Rozwi

ą

zanie dyskretne metod

ą

Eulera wstecz

y0b = y0;

% Warunki pocz

ą

tkowe

Dy0b = Dy0;
y1b = y0b + Dy0b*Tp;

% Przeliczenie warto

ś

ci pocz

ą

tkowych

% na warto

ść

drugiej próbki

tkb(1) = 0*Tp;

% Warto

ś

ci pierwszych dwóch próbek

ykb(1) = y0b;

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

10

tkb(2) = 1*Tp;
ykb(2) = y1b;

MM = (1 + 2*Tp + 5*Tp^2);
% Mianownik współczynników równania ró

ż

nicowego

a_1 = (2+2*Tp)/MM
a_0 = 1/MM
b_0 = 10*Tp^2/MM

for k = 3:(tmax/Tp)-1;

tkb(k) = (k-1)*Tp;
ykb(k) = ((2+2*Tp)/MM)*ykb(k-1) - (1/MM)*ykb(k-2) + 10*Tp^2/MM;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Symulacja równania w Simulinku
open_system('ex_1_sim')

% Otwarcie modelu Simulinka

sim('ex_1_sim',tmax)

% Wykonanie symulacji w zadanym odcinku
% czasu

close_system

% Zamkni

ę

cie modelu Simulinka

tS = tout;

% Podstawienie wektora czasu

yS = wyniki_sym(:,2);

% Pobranie z przestrzeni roboczej Matlaba
% uzyskanego rozwi

ą

zania

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zastosowanie funkcji ODE45
tp = 0;

% pocz

ą

tek symulacji

tk = tmax;

% czas ko

ń

cowy symulacji

y0a = -4;

% Warunki pocz

ą

tkowe

Dy0a = -10;

x0 = [y0; Dy0]; % Wektor warunków pocz

ą

tkowych

% OPTIONS - ustala maksymalny krok całkowania
OPTIONS = odeset('MaxStep',0.01,'InitialStep',0.01);
[tC_ode,x] = ode45(@plikfun,[tp,tk], x0, OPTIONS);
yC_ode = x(:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zastosowanie funkcji LSIM
tn = [0:0.01:tmax];

% Wektor czasu

[m,n] = size(tn);
un(1:n) = 10;

% Sygnał steruj

ą

cy

A = [0 1; -5 -2];

% Macierze opisuj

ą

ce równania

B = [0;1];

% dynamiczne

C = [1 0];
D = 0;
sys = ss(A, B, C, D);

%

yn = lsim( sys, un, tn, x0);

% Wykre

ś

lenie uzyskanych wyników w Matlabie

figure(1)
plot( tC, yC, 'k-', tS, yS, 'b:', tka, yka, 'ko', tkb, ykb, 'k*',...

tC_ode, yC_ode, 'r-', tn, yn, 'y-')

title('Rozwiazania rownania ró

ż

niczkowego')

xlabel('t [s]')
ylabel('y(k)')
grid on

Teoria sterowania

Dyskretyzacja równań różniczkowych

−

Matlab

Ostatnia aktualizacja: 2010-01-04



M. Tomera

11

% plikfun.m
% Plik zawieraj

ą

cy prawe strony równa

ń

ró

ż

niczkowych

% zapisanych w postaci macierzowej
%
% Rozwi

ą

zane równanie:

% y''(t) + 2*y'(t) + 5*y(t) = 10*1(t)
%
% Copyright (c) by Mirosław Tomera
% Date: 2007/12/11

function dxdt = plikfun(t, x)
u = 10;

% u(t) = A*1(t), gdzie A = 10

dxdt = [0 1; -5 -2]*x + [0;1]*u;