6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 90
6.
OBLICZANIE PRZEBIEGU REGULACJI
Znane są następujące metody obliczania przebiegu procesu regulacji:
Metoda analitycznego rozwiązywania róŜniczkowych równań stanu,
Metoda operatorowa ,
Metoda modelowania analogowego,
Metoda numerycznego całkowania róŜniczkowych równań stanu,
6.1. RóŜniczkowe równania stanu
W nowoczesnej teorii regulacji posługującej się techniką komputerową bardzo przydatny okazał się zapis równań róŜniczkowych w postaci układu równań róŜniczkowych pierwszego rzędu nazywanych równaniami stanu.
Równania róŜniczkowe stanu otrzymuje się bezpośrednio z podstawowych zasad bilansowania: bilansu sił działających na układ mechaniczny, bilansu masy w układach napełniania zbiorników, bilansu ładunku elektrycznego w układach elektrycznych, bilansu energii w układach cieplnych.
Zakładamy, Ŝe na układ dynamiczny przedstawiony na rys. 6.1 działają sygnały wejściowe
[u , u , L, u ]
1
2
r . Są to sygnały deterministyczne o określonym przebiegu czasowym. Stan układu jest opisany sygnałami [x , x , L, x ]
1
2
n zwanych zmiennymi stanu. W układach
mechanicznych zmiennymi stanu są: przyspieszenie, prędkość, współrzędne połoŜenia.
.
u(t)
x(t)
x(t)
y(t)
f(x,u)
∫
g(x,u)
Rys. 6.1. Równania róŜniczkowe stanu układu dynamicznego
Dynamikę układu opisuje zbiór równań róŜniczkowych o postaci:
dx (t)
i
f (x , x , L
=
, x ; u , u , L, u ; t) , i = 1, 2, ... , n
(6.1)
dt
i
1
2
n
1
2
r
Nie wszystkie zmienne stanu mogą być mierzalne. Przyjmuje się, Ŝe moŜna mierzyć inne sygnały zwane sygnałami wyjściowymi [y , y , L, y
]
1
2
m , które są związane z sygnałami
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 91
wejściowymi [u , u , L, u ]
L
1
2
r i zmiennymi stanu [x , x ,
, x ]
1
2
n za pomocą równań:
y
g (x , x , L
=
, x ; u , u , L, u ; t)
j
j
1
2
n
1
2
r
j = 1, 2, ... m., m ≤ n
(6.2)
JeŜeli ilość równań jest duŜa, to układy równań (6.1), (6.2) zapisuje się w postaci wektorowo macierzowej.
dx
= x& = f(x, ,
u t)
(6.3)
dt
y = g(x, ,
u t)
(6.4)
przy czym u, x, y, x& , f (x, , u t) , g(x, ,
u t) są macierzami jednokolumnowymi (wektorami),
kreślonymi wzorami (6.5), (6.6).
u
x
x&
y
1
1
1
1
&
u 2
x 2
y
&
x 2
2
u =
x =
x =
y =
(6.5)
M
M
M
M
u
&
r
x n
x n
ym
f (x, ,
u t)
g (x, ,
u t)
1
1
f (x, ,
u t)
2
g (x, ,
u t)
2
f (x, ,
u t) =
g(x, ,
u t) =
(6.6)
M
M
f (x, ,
u t)
x u
n
g ( , , t)
m
Wektor x nazywany jest wektorem stanu, wektor u – wektorem sterowania, y – wektorem wyjściowym.
Jeśli funkcje f oraz g są liniowe względem zmiennych x oraz u, wtedy układ równań róŜniczkowych (6.3) rzędu pierwszego moŜna zapisać w postaci układu równań (6.7): x&
=
+
+ L +
+
+
+ L +
1 (t)
a11 x1 a12 x 2
a1n x n
b11 u1 b12 u 2
b1r u r
x&
=
+
+ L +
+
+
+ L +
2 (t)
a 21 x1 a 22 x 2
a 2 n x n
b 21 u1 b22 u 2
b 2 r u r
L
L
L
L
(6.7)
x&
=
+
+ L +
+
+
+ L +
n (t)
a n1 x1 a n 2 x 2
a n n x n
b n1 u1 bn 2 u 2
b n r u r
a układ równań algebraicznych (6.4) wiąŜących sygnały wyjściowe ze zmiennymi stanu i sygnałami sterującymi moŜna zapisać w postaci układu równań (6.8):
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 92
y =
+
+ L +
+
+
+ L +
1
1
c 1 x1 c12 x 2
1
c n x n
d11 u1 d12 u 2
d1r u r
y
=
+
+ L +
+
+
+ L +
2
c 21 x1 c22 x 2
c 2 n x n
d 21 u1 d22 u 2
d 2 r u r
L
L
L
L
(6.8)
y
=
+
+ L +
+
+
+ L +
n
c n1 x1 cn 2 x 2
cn n x n
d n1 u1 dn 2 u 2
d n r u r
PowyŜsze równania (6.7), (6.8) w zapisie wektorowo – macierzowym są następujące: x& = A(t) x + B(t)u
(6.9)
y = C(t) x + D(t)u
(6.10)
przy czym
a
a
L
a
b
b
L
b
11
12
1 n
11
12
1 r
a
L
L
21
a 22
a 2 n
b21
b22
b2 r
A =
B =
M
M
M
M
M
M
a
L
L
n1
a n2
a n n
bn1 bn2
bn r
L
macierze współczynników (6.11)
c
c
L
c
d
d
L
d
11
12
1n
11
12
1r
c
c
L
c
d
d
L
d
21
22
2 n
21
22
2 r
C =
D =
M
M
M
M
M
M
c
c
L
c
L
d
d
d
1
m
m2
m n
1
m
m2
m r
D
.
u(t)
x(t)
x(t)
y(t)
∫
B
C
∫
A
Rys. 6.2. Równania róŜniczkowe stanu układu liniowego w postaci macierzowej
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 93
Elementy macierzy A zaleŜą od parametrów obiektu, a elementy B pokazują wpływ sterowań na zmienne stanu. Macierz C pokazuje, jak sygnały wyjściowe są powiązane ze zmiennymi stanu, a macierz D pokazuje bezpośredni wpływ sterowań na sygnały wyjściowe.
Zapis wektorowo – macierzowy układu równań, w którym x jest wektorem, nazywa się przedstawieniem układu równań w przestrzeni stanów
(rys. 6.2). Zapis wektorowo – macierzowy jest stosowany w pracach teoretycznych dotyczących teorii regulacji i sterowania. Przez zastosowanie zapisu macierzowego uzyskuje się znaczne skrócenie czasu potrzebnego do rozwiązywania równań róŜniczkowych. Istnieje wiele programów komputerowych przeznaczonych do obliczania przebiegu procesów nieustalonych w układach regulacji (na przykład program MATLAB).
6.2.
Metody numerycznego całkowania róŜniczkowych równań stanu Metody całkowania numerycznego (rozwiązywania równań róŜniczkowych) dają rozwiązania przybliŜone, ale mają charakter uniwersalny. Nadają się do rozwiązywania równań liniowych i nieliniowych. Warunki początkowe mogą być dowolne; zerowe i róŜne od zera. Sygnały wymuszające mogą mieć dowolny przebieg w czasie podany w postaci funkcji analitycznej lub w formie tabeli zawierającej wartości funkcji w wybranych chwilach czasu.
Do rozwiązywania równań metodą numeryczną bardziej dogodne są równania róŜniczkowe pierwszego rzędu nazywane równaniami róŜniczkowymi stanu. Zwyczajne równania róŜniczkowe wyŜszego rzędu łatwo moŜna zamienić na układ równań pierwszego rzędu wprowadzając jako dodatkowe zmienne stanu prędkości zmian wybranych sygnałów na przykład w równaniach mechaniki oprócz współrzędnych połoŜenia wprowadza się pędy ciał.
W metodach numerycznych obliczenia wykonujemy na przybliŜonych równaniach róŜnicowych, które otrzymujemy z równań róŜniczkowych przez zastąpienia przyrostów nieskończenie małych przyrostami skończonymi.
Zasady całkowania numerycznego zostaną przedstawione na przykładzie pojedynczego równania róŜniczkowego pierwszego rzędu.
Przyjmujemy, Ŝe dynamikę układu opisuje równanie róŜniczkowe stanu dx(t)
= f (x(t),u(t))
(6.12)
dt
Gdzie x(t) przedstawia zmienną stanu, natomiast u(t) – sygnał wejściowy. Obliczeniowy przedział czasu dzielimy na małe kroki całkowania, najlepiej o równej szerokości ∆t = h. Punkty podziału numerujemy liczbami całkowitymi
L
k = 0
,1
, 2
,
, i.t.d. Otrzymujemy siatkę węzłów
obliczeniowych o współrzędnych czasu
L
L
t , t , t ,
, t , t
,
= ⋅
0
1
2
k
k 1
+
; t
k h
k
.
Obliczamy przybliŜone przyrosty x
∆ (t )
k zmiennej stanu x(t) w węzłach t k . W tym celu z
równania róŜniczkowego (6.12) obliczamy przybliŜony przyrost ∆x(t) zmiennej stanu w skończonym przedziale czasu ∆t.
dx(t) = f (x(t), u(t))dt
(6.13)
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 94
Całkujemy róŜniczkę dx(t) w przedziale od t
+ =
+
k do t
t
h
k 1
k
; w rezultacie otrzymamy
następujący wzór na obliczanie przewidywanego przyrostu x
∆ (t ) = x(t
)
+ − x(t )
k
k 1
k zmiennej
stanu x(t).
t k+
t
1
k +h
∆x(tk) = x(tk+1) − x(tk) = ∫dx(t) = ∫f(x(t),u(t))dt
(6.14)
t
t
k
k
Przyrost funkcji moŜemy wyrazić w postaci iloczynu średniej wartości funkcji f (x, u) w
śr
rozwaŜanym przedziale i kroku całkowania h.
t +h
k
x
∆ (t ) = ∫ f
=
⋅
k
(x(t),u(t))dt f(x,u) h
(6.15)
śr
t k
PoniewaŜ ani nie znamy wartości obliczanej funkcji x(t) w następnym węźle o numerze (k+1), ani przebiegu tej funkcji między węzłami, dlatego nie moŜemy dokładnie obliczyć przyrostu funkcji x
∆ (t )
k , ani średniej wartości funkcji podcałkowej f (x(t), u(t))
w
śr
rozwaŜanym przedziale całkowania.
Obliczamy oszacowanie przewidywanego przyrostu funkcji według wybranej metody numerycznej z pośród wielu opracowanych przez matematyków. Znane są róŜne metody numerycznego całkowania. Najprostszą, ale mało dokładną jest metoda Eulera. Bardzo dokładną jest metoda Runge - Kutty, ale jest trudna w realizacji na komputerze.
6.2.1. Metoda Eulera (prostokątów) do całkowania równań róŜniczkowych Najprostszą i najłatwiejszą do programowania obliczeń na komputerze jest metoda Eulera.
Zakładamy, Ŝe dane jest równanie róŜniczkowe pierwszego rzędu
dx(t)
x'(t) =
= f (x(t),u(t))
(6.16)
dt
z warunkiem początkowym
x(0) = xo
W metodzie Eulera (nazywanej teŜ metodą prostokątów) przyjmuje się, Ŝe w całym kroku całkowania pochodna zmiennej stanu dx(t)/dt jest stała i równa wartości funkcji f(x(t),u(t)) przyjmowanej na początku kroku całkowania, czyli w chwili tk .
t
x(t) k 1
+
≅ f (x,u)
= f (x(t ),u(t ) =
k
k )
f (x , u )
k
k
(6.17)
dt
t k
t k
Wobec tego przyrost funkcji jest równy
t +h
k
x
∆ (t ) = ∫ f
≅
⋅ =
⋅
k
(x(t),u(t))dt f(x,u) h f(x ,u ) h
k
k
(6.18)
k
t k
a przewidywana wartość zmiennej stanu (sygnału x(t)) w następnej chwili jest równa
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 95
x(t
)
+ = x(t ) + x
∆ (t ) ≅ x(t ) + f (x ,u ) ⋅ h
k 1
k
k
k
k
k
(6.19)
x + ≅x +f (x , u ) ⋅ h
k 1
k
k
k
(6.20)
Do obliczenia przewidywanej wartości funkcji w następnej chwili czasu (t + t
∆ ) wystarczy
znać wartość pochodnej w chwili bieŜącej (t).
x(t + t
∆ ) = x(t) + x
∆ (t) ≅ x(t) + f (x(t),u(t)) t
∆
(6.21)
Metoda Eulera ma najprostszy algorytm całkowania, bardzo łatwy do realizacji na komputerze (pisania programu), niestety, ma teŜ następujące powaŜne wady: 1) mała dokładność
2) kolejne nakładanie się błędów.
2
Błąd metody Eulera w pojedynczym kroku jest rzędu h .
W innych metodach numerycznych przyjmuje się inne wzory do oceny średniej wartości funkcji f(x(t),u(t)) wewnątrz kroku całkowania. Wybór metody wiąŜe się z dokładnością otrzymanego rozwiązania oraz czasem obliczeń.
6.2.2. Metoda Runge’go - Kutty do całkowania równań róŜniczkowych Zakładamy, Ŝe dane jest równanie róŜniczkowe pierwszego rzędu
dx(t)
x'(t) =
= f (x(t),u(t))
(6.22)
dt
z warunkiem początkowym
x(t
=
o )
xo
Wybieramy krok całkowania ∆t = h . W kroku całkowania wybieramy trzy punkty, w których obliczamy wartości funkcji f(x, u). Wprowadzamy oznaczenia:
t = t + i ⋅ h
=
=
i
o
, u
u(t )
i
i , x
x(t )
i
i , (
L
i = 0
,1
, 2
,
)
Obliczamy liczby:
k(i) = h ⋅ f (x(t ), u(t )) ,
(6.23)
1
i
i
k(i) = h ⋅ f
+
+
,
(6.24)
2
([x(t ) 1k(i)
i
2 1 ), u(t
1 h
i
)
2
]
k(i) = h ⋅ f
+
+
,
(6.25)
3
([x(t ) 1k(i)
i
2
2 ), u(t
1 h
i
)
2
]
k(i) = h ⋅ f
+
+
,
(6.26)
4
([x(t ) k(i)
i
3 ), u t
h
i
]
(
)
Zgodnie ze zwykłą metodą Rungego – Kutty kolejne wartości x = x(t ) i
i obliczanej zmiennej
stanu x(t) znajdujemy ze wzoru
x + =
+ =
+ ∆
= +∆
i 1
x(ti 1) x(ti )
x(ti ) xi
xi
(6.27)
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 96
gdzie
1
x
∆ =
=
+
+
+
i
x(ti )
k
2k
2k
k
(
L
i = 0
,1
, 2
,
).
(6.28)
6 [ (i)
(i)
(i)
(i)
1
2
3
4 ]
Wyprowadzenie powyŜszych wzorów jest dość zawiłe.
Metoda Runge’go Kutty jest bardziej skomplikowana od metody Eulera ale,ale bardzo dokładna. Nie ma wady sumowania błędów. Błąd metody Rungego – Kutty dla kaŜdego kroku jest rzędu 5
h przy załoŜeniu, Ŝe funkcja jest ciągła i ma ciągłe pochodne do 5 – go rzędu.
6.2.3. Metoda numeryczna Runge’go – Kutty do całkowania układów równań
róŜniczkowych
Dany jest układ trzech równań róŜniczkowych
dx (t)
1
= f1(x (t),x (t),x (t),u (t)u (t)u (t)
1
2
3
1
2
3
)
dt
dx (t)
2
= f2(x (t),x (t),x (t),u (t)u (t)u (t)
1
2
3
1
2
3
)
dt
dx (t)
3
= f3(x (t),x (t),x (t),u (t)u (t)u (t)
1
2
3
1
2
3
)
dt
W zapisie wektorowym
dx(t) = f(x(t), (
u t))
dt
i warunki początkowe
=
1
x (to ) x o
1
x
=
2 (to )
x2o
x
=
3 (to )
x o
3
W zapisie wektorowym
x( to ) = xo
Ustalamy krok całkowania t
∆ = h i wprowadzamy oznaczenia
t = t + i ⋅ h
=
=
∆
=
i
o
, u
u (t )
,
j i
j
i , x
x (t )
,
j i
j
i ,
x
+ −
,
j i
x ,j i 1 x ,j i
dla (
L
i = 0
,1
, 2
, ,
3 ,
4
), (j =1, 2, 3)
W zapisie wektorowym
i
u = (
u t )
i
x = (t )
i
x i
x
∆
=
+ −
i
x i 1 x i
Następnie obliczamy parametry pierwszego kroku o numerze i = 0
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 97
(0)
k
= h ⋅ f( (
u t ), x(t )
1
o
o )
(0)
k
= h ⋅ f u
+
x
+ k
2
( (t 1h),[ (t ) 1 (0)
o
o
2
2 1
)]
(0)
k
= h ⋅ f u
+
x
+ k
3
( (t 1h),[ (t ) 1 (0)]
o
o
2
2
2
)
(0)
k
= h ⋅ f u
+
x
+ k
4
( (t h),[ (t ) (0)
o
o
3
)]
Zgodnie z metodą Rungego – Kutty znajdujemy przybliŜoną wartość przyrostu zmiennych stanu z wzoru
1
x
∆
= ⋅
+ ⋅
+ ⋅
+
o
k
2 k
2 k
k
6 ( (0)
(0)
(0)
(0)
1
2
3
4 )
Stąd
=
+ ∆
1
x
xo
xo
W dalszym ciągu przyjmując
1
x jako dane wyjściowe, postępując podobnie znajdujemy x2 .
Analogicznie obliczamy xi (
L
i = ,
3
,
4 ,
5
).
Do rozpoczęcia obliczeń naleŜy znać:
•
Równanie róŜniczkowe zawierające funkcję f(x, u) zadaną w formie wzoru lub tabeli.
Funkcja f(x, u) przedstawia zaleŜność pochodnej zmiennej stanu od tej zmiennej x(t) oraz od sygnału wymuszającego u(t).
•
Czasowy przebieg sygnału wymuszającego u(t).
•
Wartość początkowa zmiennej stanu x(0) (sygnału obliczanego).
Metoda ta zostanie przedstawiona na trzech przykładach: pierwszy – obliczanie charakterystyki czasowej elementu inercyjnego 1 – go rzędu; drugi – obliczanie przebiegu regulacji obiektu inercyjnego 1 – go rzędu z opóźnieniem czasowym z zastosowaniem regulatora P; trzeci – obliczanie przebiegu regulacji prędkości obrotowej silnika wysokopręŜnego z zastosowaniem regulatora PI z uwzględnieniem nieliniowości obiektu regulacji i ograniczeń regulatora.
6.3.
Przykłady obliczeń numerycznych
A)
Obliczanie charakterystyki czasowej elementu inercyjnego
1 – go rzędu metodą numeryczną Eulera
Dane do obliczeń:
y(s)
K
Transmitancja operatorowa: G s
( ) =
=
x s
( )
Ts +1
Przebieg sygnału wymuszającego: u(t) = x(t) = (
1 t)
Wartości parametrów transmitancji: K = 1, T = 0,5 [s]
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 98
Wartość początkowa sygnału wyjściowego: y(0) = 0
Sygnał wyjściowy y(t) jest równocześnie zmienną stanu, dlatego w tych obliczeniach przyjęto oznaczenia sygnałów przyjmowane tradycyjnie przy opisie elementów automatyki.
Przygotowanie równań do obliczeń
Odtwarzamy równanie róŜniczkowe odpowiadające podanej transmitancji operatorowej.
y(s)
k
=
,
x s
( )
Ts +1
(Ts + )
1 y s
( ) = k x s
( )
1
s y s
( ) =
[k⋅x s()− y s()]
T
Odtworzone równanie róŜniczkowe odpowiadające podanej transmitancji operatorowej jest następujące:
d
1
y(t) =
[k⋅x(t)− y(t)]
dt
T
Odtworzone równanie w postaci całkowej odpowiadające podanej transmitancji operatorowej jest następujące:
1 1
y s
( ) = ⋅
[k⋅x s()− y s()]
s T
t
y(t) = y(0 + ∫ 1
)
[k ⋅x(t) − y(t)]dt
T
0
Równanie róŜnicowe do obliczania wartości funkcji
1
y(t + ∆t) = y(t) + y
∆ (t) = y(t) + [k ⋅ x(t) − y(t)]∆t
T
Dla wymuszenia skokowego jednostkowego x(t) = (
1 t) otrzymamy następujące równanie
róŜnicowe:
y(k + )
1 = y(k) + [1− y(k ] h
) ⋅
T
gdzie:
L
k = 0
,1
, 2
,
Przyjmujemy: krok całkowania h= 0,2 [sek], stałą czasową T = 0,5 [sek]
y(k + )
1 = y(k) + [1− y(k)] ,
0 2 s
[ ek]
⋅
5
,
0
s
[ ek]
Wyniki obliczeń charakterystyki czasowej przedstawiono w tabeli 6.1 oraz na rys. 6.3.
6.
OBLICZANIE PRZEBIEGU REGULACJI
str. 99
Tabela 6.1.
Numeryczne obliczanie charakterystyki czasowej elementu inercyjnego 1 – go rzędu metodą Eulera w arkuszu kalkulacyjnym EXCEL’a.
A
B
C
D
E
F
G
1
2
∆t =
0,2
dy/dt = (x-y)/T
3
T =
0,5
∆y = (dy/dt)*∆t
4
y(k) = y(k-1)+∆y(k)
5
Nr
czas
sygn.
sygn.
6
[s]
wej.
wyj.
7
k
t
x
x-y
dy/dt
∆y
y(k)
8
1
0
0
0
9
2
0,2
1
1,00
2,00
0,40
0,40
10
3
0,4
1
0,60
1,20
0,24
0,64
11
4
0,6
1
0,360
0,720
0,144
0,784
12
5
0,8
1
0,216
0,432
0,086
0,870
13
6
1
1
0,1296
0,2592
0,0518
0,9222
14
7
1,2
1
0,0778
0,1555
0,0311
0,9533
15
8
1,4
1
0,0467
0,0933
0,0187
0,9720
16
9
1,6
1
0,0280
0,0560
0,0112
0,9832
17
10
1,8
1
0,0168
0,0336
0,0067
0,9899
18
11
2
1
0,0101
0,0202
0,0040
0,9940
19
12
2,2
1
0,0060
0,0121
0,0024
0,9964
20
13
2,4
1
0,0036
0,0073
0,0015
0,9978
21
14
2,6
1
0,0022
0,0044
0,0009
0,9987
y
rozwiązanie metodą Eulera, krok całkowania = 0,2 sek
1.0
0.8
0.6
rozwiązanie dokładne: y = 1 - exp(- t/T), T = 0,5 sek
0.4
0.2
czas t
0.0
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0
2.2
sekund
Rys. 6.3. Porównanie przebiegu charakterystyki czasowej obliczonej metodą numeryczną Eulera (prostokątów) z rozwiązaniem dokładnym.