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.