ZÁKLADY ANALÝZY A SYNTÉZY VE STAVOVÉM PROSTORU

TECHNICKÁ UNIVERZITA V LIBERCI Hálkova 6, 461 17 Liberec 1, CZ Fakulta mechatroniky a mezioborových inženýrských studií

Teorie automatického řízení II.

ZÁKLADY ANALÝZY A SYNTÉZY VE STAVOVÉM PROSTORU Studijní materiály Doc. Ing. Osvald Modrlák, CSc. Srpen 2004

Katedra řídicí techniky

Teorie automatického řízení II

Základy stavového řízení

Obsah Úvod

2

1. Stavový popis 2 1.1 Základní pojmy stavového popisu ................................. 5 1.2 Řešení stavových rovnic pomocí L- a Z-transformace .................... 7 1.3 Vztah mezi vnějším a vnitřním popisem ............................ 9 1.3.1 Normální forma řiditelnosti .................................. 10 1.3.2 Normální forma rekonstruovatelnosti ........................... 14 1.4 Vlastnosti stavových modelů dynamické soustavy .......................17 1.4.1 Stavový prostor, vektor stavu, stavová trajektorie ................... 17 1.4.2 Dosažitelnost, řiditelnost, pozorovatelnost, rekonstruovatelnost .......... 18 1.5 Transformace stavových rovnic .................................. 20 1.6 Model poruch ve stavovém popisu ................................ 21 2. Estimace stavu 25 2.1 Deterministický estimátor stavu .................................. 25 2.2 Kalmanova filtrace .......................................... 30 3. Syntéza stavových regulátorů 36 3.1 Struktura stavových regulátorů ................................... 36 3.2 Návrh stavového regulátoru metodou zadání pólů ...................... 37 3.2.1 Návrh stavového regulátoru .................................. 38 3.2.2 Stavový regulátor pro vyrovnání trvalé poruchy .................... 40 3.2.3 Diskrétní stavový regulátor pro vyrovnání trvalé poruchy .............. 44 3.3 Návrh stavového regulátoru podle kvadratického kritéria .................. 56 3.3.1 Rozšířený stavový popis .................................... 56 3.3.2 Podmínky pro vyrovnání skoků žádaných hodnot ................... 57 3.3.3 Algoritmy a softwarová podpora .............................. 58 3.3.4 Vyrovnání trvalých poruch pro kvadraticky optimální řízení ............ 64 4. Softwarová podpora analýzy a syntézy ve stavovém prostoru 74 4.1 Stavový popis a transformace ................................... 74 4.2 Návrh estimátorů ............................................ 77 4.3 Návrh stavových regulátorů ..................................... 78 Literatura

82

Je mojí milou povinností poděkovat pánům Doc. Ing. Josefu Janečkovi CSc. a Radku Reifovi za pečlivé čtení rukopisu a cenné připomínky.

Doc. Ing. Modrlák Osvald, CSc.

1

27.8.2004



ÚVOD Cílem této části je seznámit studenty se základy analýzy a syntézy ve stavovém prostoru. Omezíme výklad pouze na lineární systémy časově invariantní a výklad bude převážně prováděn souběžně jak pro spojitý tak i diskrétní popis. Hlubší znalosti o metodách analýzy a syntézy ve stavovém prostoru mohou studenti získat v kurzu "Řízení ve stavovém prostoru" v pátém ročníku a samozřejmě v knižních a časopiseckých publikacích [1-5]. Z přednášek "Teorie řízení I." je známo, že k matematickému popisu dynamických systémů s koncentrovanými parametry lze použít vnější nebo vnitřní popis. 1) Vnější popis systému je vyjádření dynamických vlastností pomocí relací mezi vstupní a výstupní veličinou. Tyto relace mohou být vyjádřeny ve tvaru diferenciálních rovnic, pro lineární časově invariantní systémy obrazovým přenosem atd. 2) Vnitřní popis systému chápeme jako relaci mezi vstupní veličinou u(t) a stavem systému x(t) a výstupní veličinou y(t). Tato relace má formu soustavy diferenciálních rovnic prvního řádu. Hovoří se pak o stavových rovnicích systému.

1. STAVOVÝ POPIS Z matematicko-fyzikální analýzy dynamických systémů při aplikaci makroskopických bilancí hmoty a energie dostáváme zpravidla soustavu lineárních rovnic prvního řádu, tedy přímo stavový popis. Každá stavová veličina pak představuje konkrétní fyzikální veličinu a struktura stavových rovnic pak vypovídá o vzájemných vazbách mezi stavovými-fyzikálními veličinami. Ukážeme to na následujícím příkladu ohřevu vody v průtokovém ohřívači. Příklad 1

Matematicko-fyzikální analýza průtokového ohřívače vody.

Dynamickou soustavu na obr.1-1 tvoří průtokový ohřívač PO. Vstupní veličinou je výkon P topné spirály TS. Výstupní veličinou je teplota vody T měřená čidlem teploty ČT. ČT

TS

Obr.1-1 Dynamická soustava PO Pro účely matematicko-fyzikální bude uvažována zjednodušená soustava podle obr.12 a vychází se z podmínek: dokonalé tepelné izolace, konstantního přítoku a objemu vody v ohřívači, promíchávání vody, teplota vody Tout a teplota ohřívače Th nezávisí na prostorových souřadnicích (Tout = Tout(t), Th = Th(t)). Doc. Ing. Modrlák Osvald, CSc.

2

27.8.2004



V … objem vody v ohřívači ρ … měrná hustota vody M … množství protékající vody Vh … objem topného tělesa TIN … teplota vody na vstupu TOUT … teplota vody na výstupu S … teplosměnná plocha topného tělesa … měrné specifické teplo ch topného tělesa ρh … měrná hustota topného tělesa α … koeficient přestupu tepla c … měrné specifické teplo v.

Obr.1.2 Matematicko-fyzikální analýza

Řešení: 1) Za uvedených předpokladů vycházíme ze zákona zachování energie, podle kterého platí Časová změna akum. Rychlost přívodu  Rychlost odvodu   =  −  tepelné energie , tepelné energie t epelné energie      

což je možno vyjádřit rovnicí

dU (t ) = PIN − POUT , dt kde je: U(t) PIN POUT

vnitřní (tepelná) energie rychlost přívodu tepelné energie - přiváděný tepelný výkon rychlost odvádění tepelné energie - odváděný tepelný výkon

2) Makroskopickou bilanci provedeme pro hmotu ohřívače a vody. Tepelná kapacita vody v ohřívači je rovna U = U(t) = c . ρ . V . TOUT(t). a tepelná kapacita tělesa ohřívače je rovna U h = U h (t ) = c h ⋅ ρ ⋅ h Vh ⋅ Th (t ). Nyní můžeme aplikovat zákon zachování energie na tělese ohřívače d [ch ⋅ ρ h ⋅ Vh ⋅ Th (t )] = α ⋅ S ⋅ [TOUT (t ) − Th (t )] + P(t ), dt

(1)

Aplikací rovnice energie na objem vody v průtokovém ohřívači dostaneme d [c ⋅ ρ ⋅ V ⋅ TOUT (t )] = M ⋅ c ⋅ [TIN (t ) − TOUT (t )] + α ⋅ S ⋅ [Th (t ) − TOUT (t )], dt


3

(2)

27.8.2004



Jednoduchou úpravou rovnic (1) a (2) vytvoříme soustavu dvou diferenciálních rovnic prvního řádu ve tvaru d P(t ) αS αS Th (t ) = − Th (t ) + TOUT (t ) + dt c h ρ hVh c h ρ hVh c h ρ hVh d Mc αS αS + Mc TOUT (t ) = Th (t ) − TOUT (t ) + TIN (t ). dt cρV cρV cρV Soustava rovnic obsahuje dvě budící funkce, P(t) a TIN(t). Akční veličinou je elektrický příkon u(t) = P(t), který lze měnit pomocí akčního členu. Druhá vstupní veličina je teplota vstupující vody TIN(t), kterou ale nemohu cíleně ovlivňovat a není ani měřena. Proto ji budeme pokládat za neměřenou poruchovou veličinu d(t) = TIN(t). Zavedeme-li následující označení  d  dt Th (t )   x&1 (t )   x (t )  Th (t )   P(t )  u (t )  , x(t ) =  1  =  , u b (t ) =  x& =  =   = , TIN (t ) d (t )  d T (t )  x& 2 (t )  x 2 (t ) TOUT (t )  dt OUT  pak je možno nalezenou soustavu rovnic prvního řádu maticově zapsat do tvaru

αS αS    1  , , 0  −    c ρV c ρV u (t )   x (t )  c ρ V  x& (t )  . x& (t ) =  1  =  h h h h h h  ⋅  1  +  h h h .  x& 2 (t )  αS , − αS + Mc   x 2 (t )  0 , Mc  d (t )   cρV cρV  cρV    Je třeba připomenout, že stavová veličina x1(t) je teplota topného elementu (ohřívače) a druhá stavová veličina x2(t) je teplota vody na výstupu z průtokového ohřívače. Naším úkolem bylo nalézti výstupní teplotu z ohřívače, takže z vektoru stavu ji vypočteme následovně  x1 (t )  y (t ) = x 2 (t ) = [0,1] ⋅   , x& (t ) = Ax(t ) + Bu(t ) ,  x 2 (t ) kde jsou A matice koeficientů soustavy, B matice buzení

αS αS   − c ρ V , c ρ V  − a ; a  A =  h h h h h h  =  11 12  , αS + Mc   a 21 ;−a 22   αS  cρV , − cρV   

1   , 0 c ρ V   B =  h h h Mc  .  0 ,  c ρ V  

Chceme-li formálně oddělit účinek akční a poruchové veličiny, přepíšeme stavovou rovnici do tvaru αS αS   , −  1  0   c h ρ hVh c h ρ hVh   x1 (t )    x&1 (t )   ⋅ x& (t ) =  (4) = + c h ρ hVh ⋅ u (t ) +  M .d (t ) .     &  x 2 (t )  αS , − αS + Mc   x 2 (t )  0   ρV     cρV cρV   Zjednodušeným maticovým zápisem dostaneme Doc. Ing. Modrlák Osvald, CSc.

4

27.8.2004



x& = Ax(t ) + B U u (t ) + B d d (t ),

 1  0    kde BU = c h ρ hVh , B d =  M  , BU, Bd jsou matice buzení akční a poruchové veličiny.      0   ρV 

1.1 ZÁKLADNÍ POJMY STAVOVÉHO POPISU Pro lineární časově invariantní soustavu s jedním vstupem a jedním výstupem má stavová rovnice tvar a) pro spojitý popis b) pro diskrétní popis x& (t ) = Ax (t ) + Bu(t ) ; rovnice výstupu y(t) = Cx(t) + Du(t), kde je A B x(t) u(t) C D

x(k + 1) = Mx(k) + Nu(k)

(1 – 1,2)

y(k) = Cx(k) + Du(k),

(1 – 3,4)

matice soustavy [n, n], M matice [n, n] diskrétních koeficientů matice [n, p] buzení, N matice [n, p] buzení vektor [n, 1] stavových veličin, x(k) vektor [n, 1] disk. stavových veličin, vektor [p, 1] buzení, u(k) vektor [p, 1] diskrétního buzení matice [1, n] výstupu p počet budících veličin matice [1, 1] převodu y(t), y(k) výstupní veličina

Pro stavovou rovnici se kreslí blokové schéma pomocí integračních a diskrétních zpožďovacích bloků viz obr.1-1a, b.

D x′(t)

x(t)

B

C

∫ dt

u(t)

y (t )

A Obr.1a Bloková struktura stavových rovnic spojitého systému

D x(k +1)

N u(k)

x(k) 1 z

C

y (k )

M Obr.1b Bloková struktura stavových rovnic diskrétního popisu Doc. Ing. Modrlák Osvald, CSc.

5

27.8.2004


Příklad 1.1


Uvažujme dynamický systém x(k + 1) = Mx(k) + N u(k) 0,6;−0,04 − 0,5 M = ; N=   ; C = [1,0; 0,2]; D = 0. 1,0; 0,1   1,0 

Určete: 1) Stavové diskrétní rovnice a rovnici pro y(k). 2) Spočtěte vektor stavu x(k) a y(k) pro u(k) = {1, 1, 0, -1}, x(0) = 0; k = 0, 1, 2, 3. 3) Nakreslete blokové schéma stavové rovnice. Řešení: 1) Stavové rovnice jsou x1 (k + 1) = 0,6 x1 (k ) − 0,04 x 2 (k ) − 0,5u (k ) x 2 (k + 1) =

x1 (k ) + 0,1x 2 (k ) + u (k ) y(k) = x1(k) +0,2 x2(k).

2) Výpočet diskrétních hodnot vektoru stavu x(k) a výstupní veličiny y(k) k = 0;

y(0) = 0,

− 0,5 0,6;−0,04 0 − 0,5 0,6;−0,04  x1 (0)  − 0,5 ⋅1 =  + ⋅ u ( 0) =  ⋅ + x(1) =   ,        1,0  1,0; 0,1  0  1,0  1,0; 0,1   x 2 (0)  1,0 

k = 1;

y(1) = x1(1) + 0,2x2(1) = -0,5 + 0,2 = -0,3

− 0,84 0,6;−0,04 − 0,5  − 0,5 0,6;−0,04  x1 (1)  − 0,5 ⋅1 =  + ⋅ u (1) =  ⋅ + x ( 2) =   ,        0,6  1,0; 0,1   1,0  1,0  1,0; 0,1   x2 (1)  1,0 

k = 2;

y(2) = x1(1) + 0,2x2(2) = -0,84 + 0,12 = -0,72

− 0,528  0,6;−0,04 − 0,84  − 0,5 0,6;−0,04  x1 (2)  − 0,5 ⋅0 =  + ⋅ u ( 2) =  ⋅ + x (3) =   ,       − 0,780 1,0; 0,1   0,60  1,0  1,0; 0,1   x2 (2)  1,0 

k = 3;

y(3) = x1(3) + 0,2x2(3) = -0,528 + 0,2⋅ 0,78 = -0,684

+ 0,214 0,6;−0,04 − 0,528  − 0,5 0,6;−0,04  x1 (3)  − 0,5 ⋅ (−1) =  + ⋅ u (3) =  ⋅ + x ( 4) =   .       − 1,606  1,0; 0,1  − 0,780  1,0  1,0; 0,1   x2 (3)  1,0 


6

27.8.2004



3) Bloková struktura je na obr.2.

-0,5

Σ

0,6

x1(k+1)

1 z

x1(k) 1,0

1 -,04 1,0

Σ

y(k)

Σ

u(k)

x2(k+1)

1 z

x2(k)

0,2

0,1 Obr.2 Bloková struktura

1.2 ŘEŠENÍ STAVOVÝCH ROVNIC POMOCÍ L- a Z-TRANSFORMACE Při řešení stavových rovnic se omezíme pouze na metodu založenou na Laplaceově obrazu stavových veličin, respektive jejich Z-obrazu. Postup ukážeme na spojitém i diskrétním popisu současně. Uvažujme stavový popis dle (1 – 1,2) a (1 – 3,4) ve tvaru a) pro spojitý popis

b) pro diskrétní popis

x& (t ) = Ax (t ) + Bu(t ) ; rovnice výstupu

x(k + 1) = Mx(k) + Nu(k)

y(t) = Cx(t) + Du(t),

y(k) = Cx(k) + Du(k).

Aplikujeme-li L a Z–transformaci vždy na obě strany stavového popisu, dostaneme sX(s) – x(0) = A⋅X(s) + B⋅U(s)

zX(z) – x(0) = M⋅X(z) + N⋅U(z)

Separací X(s) a X(z) obdržíme (sI – A) ⋅X(s) = x(0) + B⋅U(s)

(zI – M) ⋅X(z) = x(0) + N⋅U(z)

kde je I… jednotková matice [n, n]. Obrazy X(s) a X(z) stavového vektoru jsou X(s) = (sI – A)-1 ⋅[x(0) + B⋅U(s)];

X(z) = (zI – M)-1 ⋅[x(0) + N⋅U(z)].

Připomeneme si, že inverzní matice jsou rovny


7

27.8.2004


( sI − A) −1 =


adj ( sI − A) adj ( zI − M ) ; ( zI − M ) −1 = . det( sI − A) det( zI − M )

Obrazy X(s) a X(z) stavového vektoru můžeme upravit do tvaru X( s ) =

adj ( sI − A) adj ( zI − M ) ⋅ [x(0) + B ⋅ U( s )]; X ( z ) = ⋅ [ x (0) + N ⋅ U ( z )] . (1 – 5,6) det( sI − A) det( zI − M )

det(sI – A) = a0 + a1s + ... + an-1sn-1 + sn; (1 – 7,8) det(zI – M) = a0 + a1z + ... + an-1zn-1 + zn jsou charakteristickými polynomy stavové rovnice a kořeny těchto polynomů jsou vlastní čísla matice A resp. M. Determinanty

Obrazy Y(s), Y(z) výstupu jsou rovny (dále buzení u uvažujeme jako skalár u) Y(s) = CX(s) + DU(s); Příklad 1.2.1

Y(s) = CX(z) + DU(z).

(1 – 9,10)

Dynamická soustava je popsána stavovou rovnicí

 x&1 (t )   0; 1   x1 (t )  0 x& (t ) =   +   ⋅ u (t ); = ⋅  x& 2 (t ) − 0,5;−1,5  x2 (t ) 1 

 x1 (t )  y (t ) = [0,75;0]⋅  .  x2 (t )

Určete stavový vektor analyticky pomocí Laplaceovy transformace pro budící funkci u(t) = 1(t) a nenulové počáteční podmínky: x1(0) = 0,5 ; x2(0) = -1. Řešení: Obraz stavového vektoru určíme podle (1 – 5) X (s) =

adj ( sI − A) ⋅ [x (0) + B ⋅ U ( s )]; det( sI − A)

Matice (sI - A) a její determinant je roven  s;0   0; 1   s; − 1   s; − 1  ( sI − A) =   −  = ; det    + 0,5; s + 1,5 = s ( s + 1,5) + 0,5. 0; s  − 0,5;−1,5 + 0,5; s + 1,5   Adjungovaná matice je rovna  s; − 1   s + 1,5; + 1 adj ( sI − A) = adj   = + 0,5; s + 1,5  − 0,5 ; s 

Dosazením do (1 – 5) dostaneme  s + 1,5; + 1 − 0,5 ; s   0,5 adj ( sI − A)    + 0U ( s ) . X ( s) = ⋅ [ x (0) + B ⋅ U ( s )] = 2   det( sI − A) s + 1,5s + 0,5 − 1  1  


8

27.8.2004



Rozepsáním získáme obraz stavového vektoru ve tvaru 1 s + 1,5    ( s + 1)( s + 0,5) ; ( s + 1)( s + 0,5)  0,5 0 1   ⋅   +   ⋅  = X ( s) =  s − 0,5   − 1  1  s   ( s + 1)( s + 0,5) ; ( s + 1)( s + 0,5)    1  0,5s − 0,25     ( s + 1)( s + 0,5)   s ( s + 1)( s + 0,5)  +  = X P ( s) + X U ( s) = 1  − ( s + 0,25)     ( s + 1)( s + 0,5)   ( s + 1)( s + 0,5)     

Stavový vektor X(s) je dán součtem vektoru XP(s) a XU(s). Stavový vektor XP(s) popisuje dynamické účinky vyvolané nenulovými počátečními podmínkami a je roven  0,5s − 0,25   ( s + 1)( s + 0,5)   x (t )  − 0,5 exp(−t ) + exp(−0.5t )   →  1P  =  X P ( s) =  .  − ( s + 0,25)   x 2 P (t ) − 1,5 exp(−t ) + 0,5 exp(−0.5t )  ( s + 1)( s + 0,5)    Stavový vektor XU(s) popisuje dynamické účinky akční veličiny a je roven 1    s ( s + 1)( s + 0,5)   x (t )  (2 + 2 exp(−t ) − 4 exp(−0,5t ) )  →  1U  =  X U ( s) =  . 1    x 2U (t ) (+ 2 exp(−t ) − 2 exp(−0,5t ) )   ( s + 1)( s + 0,5)    Obraz výstupu a jeho předmět y(t) je roven 1    s ( s + 1)( s + 0,5)  0,75 = Y ( s ) = [0,75; 0] ⋅  → y (t ) = (1,5 + 1,5 exp(−t ) − 3 exp(−0,5t ) )1(t ). 1   s ( s + 1)( s + 0,5)  ( s + 1)( s + 0,5)    Konec příkladu

1.3 VZTAH MEZI VNĚJŠÍM A VNITŘNÍM POPISEM Uvažujme soustavu s jedním vstupem a jedním výstupem. Je známo, že na základě měření a jeho vyhodnocení získáme obrazový nebo diskrétní přenos zvolené struktury. Požadavky na kvalitu regulačních pochodů nás však často nutí provést návrh regulátoru ve stavo-


9

27.8.2004



vém prostoru. Stavový popis totiž poskytuje informace o stavech systému a tyto informace je možno pochopitelně využít při řízení a regulaci. Je tedy třeba vysvětlit, jakým způsobem přejdeme z obrazového přenosu na stavový popis (hledáme matici soustavy, buzení a výstupní matici). Jistě si dokážeme představit vektor stavových veličin, který bude obsahovat derivace výstupního signálu y(t). Převod diferenciální rovnice "n"-tého řádu na soustavu n rovnic prvního řádu realizuje např. metoda snižování řádu derivace nebo metoda postupných integrací. Budeme tedy vycházet ze schémat řešení diferenciálních rovnic pomocí integrátorů, sumátorů a násobení signálů konstantami. Omezíme se pouze na dvě metody 1) Zobecněnou metodu snižování řádu derivace, jejíž aplikace vede na stavový popis, který se označuje jako normální forma řiditelnosti. 2) Modifikovaná metoda postupných integrací, která vede na strukturu stavového popisu, nazývanou normální forma rekonstruovatelnosti.

1.3.1 Normální forma řiditelnosti 1) Spojitý popis. Uvažujme obrazový přenos F (s) =

B( s ) bn −1 s n −1 + bn − 2 s n − 2 + ... + b1 s + b0 Y ( s ) . = = A( s ) U (s) s n + a n −1 s n −1 + ... + a1 s + a0

(1 – 11)

Ze základů spojitého řízení je známo, že metoda snižování řádu derivace se hodí zejména pro řešení diferenciálních rovnic, které neobsahují derivace na pravé straně. Výhodou je, že ve schématu řešení existují signály odpovídající derivacím výstupního signálu. Vzhledem k uvažovanému přenosu (1 – 11) je tedy nutno provést modifikaci této metody. Modifikace spočívá a) v zavedení pomocné proměnné x(t), která nemá derivace vstupního signálu na pravé straně a b) ve vyjádření výstupního signálu y(t) jako váženého součtu derivací pomocné veličiny x(t). a) Obraz výstupní veličiny Y(s) můžeme vyjádřit z rovnosti A( s )Y ( s ) = B( s )U ( s ) → Y ( s ) = B( s ) ⋅

U ( s) = B ( s ) ⋅ X ( s ), A( s )

(1 – 12)

kde X(s) je pomocná proměnná. Pro pomocnou proměnnou X(s) platí X ( s) =

U ( s) → A( s ) X ( s ) = U ( s ) . A( s )

(1 – 13)

Z rovnosti (1 – 13) můžeme získat diferenciální rovnici pro pomocnou proměnnou x(t) ve tvaru (1 – 14) x(n) + an – 1x(n – 1) + ... + a1x(1) + a0x = u . b) Obraz výstupu Y(s) pak určíme z rovnosti (1 – 12) Y(s) = B(s)X(s) = (bn – 1sn – 1 + bn – 2sn – 2 + ... + b1s + b0)X(s)


10

(1 – 15)

27.8.2004



Výstup y(t) v časové oblasti je pak roven váženému součtu y(t) = bn – 1x(n – 1) + bn – 2x(n – 2) + ... + b1x(1) + b0x .

(1 – 16)

Stavové veličiny volíme jako derivace pomocné funkce x(t) x1 (t ) = x(t )

x&1 (t ) = x2 (t ) x& 2 (t ) = x3 (t )

x 2 (t ) = x(t )' x3 (t ) = x(t )' ' . . x n (t ) = x(t ) ( n −1)

⇒

x&3 (t ) = x4 (t ) . . x& n (t ) = −an−1 xn − ⋅ ⋅ ⋅ − a1 x2 − a0 x1 + u.

(1 – 17)

Dosadíme-li do (1 – 16) za derivace pomocné funkce x(t) stavové veličiny z (1 – 17), získáme pro výstupní veličinu y(t) rovnost y(t) = b0x1 + b1x2 + ... + bn – 2xn – 1 + bn – 1xn .

(1 – 18)

Struktura vektoru stavu x, matice systému A, matice buzení B a matice výstupu C (D = 0) je  x1   0 x   0  2    x3   0 x=  ; A= M   M  xn−1   0    − a0  xn 

Poznámka

1

0

L

0

1

L

0

0

L

0 − a1

0

L

− a2 L

0  0  0  ;  1   − a n−1 

0  0    0  B= ; M  0    1  C = [b0, b1, ..., bn – 1]

(1 – 19)

Nastane-li případ, že polynom B(s) obrazového přenosu (1 – 11) je roven B(s) = b0 (koeficienty b1 = b2 = ...= bn=0), pak x1(t) = y(t) a platí: x1 (t ) = y (t )

x&1 (t ) = x2 (t )

x2 (t ) = y (1) (t )

x& 2 (t ) = x3 (t )

x3 (t ) = y ( 2 ) (t ) ⇒

x&3 (t ) = x4 (t )

M xn (t ) = y ( n−1) (t )

x& n (t ) = − an−1 xn − ⋅ ⋅ ⋅ − a1 x2 − a0 x1 + u

Stavový prostor se pak nazývá fázovým prostorem. Stavové veličiny mají fyzikální význam a představuje-li y(t) dráhu, pak další složky jsou rychlost, zrychlení atd.


11

27.8.2004



Programové schéma normální formy řiditelnosti, které je realizováno bloky integrátorů, sumátorů a násobení konstantou, je na obr.1.3

1,5 určete stavovou rovnici v normální formě ři2s + 3s + 1 ditelnosti a programové schéma realizujte v SIMULINKu. 0,75 Řešení: Přenos upravíme F ( s ) = 2 . Zřejmě a0 = 0,5; a1 = 1,5; b0 = 0,75. s + 1,5s + 0,5 Podle (1 – 17) a (1 – 19) platí  x1 (t )   x& (t )   0 ; 1   x1 (t )  0 x& (t ) =  1  =  ⋅ +   ⋅ u (t ); y (t ) = [0,75;0] ⋅   ;   x 2 (t )  x& 2 (t ) − 0,5;−1,5  x 2 (t ) 1 

Příklad 1.3.1

K přenosu F ( s ) =

2

Programové schéma je na obr.1.3.1-1. Přechodová funkce a průběhy x1(t), x2(t) je na obr.1.3.1-2. x2 y

x1

x1

y

x2

Obr.1.3.1-1 Programové schéma

Konec příkladu


Obr.1.3.1-2 Odezva na u(t) = 1(t) 12

27.8.2004



2) Diskrétní popis. Uvažujme diskrétní (impulsní) přenos

B( z ) bn −1 z n −1 + bn − 2 z n − 2 + ... + b1 z + b0 Y ( z ) G( z) = = = , A( z ) U ( z) z n + a n −1 z n −1 + ... + a1 z + a 0 kde je

(1 – 20)

A(z) je polynom stupně n B(z) je polynom stupně n - 1

Této rovnici odpovídá diferenční rovnice

y (k + n) = − a n −1 y (k + n − 1) − ... − a1 y (k + 1) − a0 y (k ) + + bn −1u (k + n − 1) + bn − 2 u (k + n − 2) + ... + b1u (k + 1) + b0 u (k ).

(1 – 21)

Schéma řešení této diferenční rovnice je na obr.1.4. Zpoždění o jeden krokovací interval se realizuje blokem z-1. Za stavové veličiny jsou voleny jednotlivé výstupy na zpožďovacích členech. Tvary matic M, N a vektorů jsou shodné s (1 – 19)  x1 (k )   0  x (k )   0  2    x3 (k )   0 x=  ; M = M   M  xn−1 (k )  0    − a0  xn (k ) 

1

0

L

0

1

L

0

0

L

0

0

L

− a1

− a2 L

0  0   0  0    0  0  ; N= ;  M  0  1     − a n−1  1  C = [b0, b1, ..., bn – 1];

(1 – 22)

Σ − bn−2

− bn−1

u (k )

xn ( k )

Σ

z

−1

z − an−1

− b0 x2 (k )

−1

y (k )

xn−1 (k ) − an−2

z

−1

x1 ( k )

− a0

Σ Obr.1.4 Programové schéma diskrétní formy řiditelnosti pro diskrétní popis


13

27.8.2004



0.07339 z + 0.05716 určete: z 2 − 1.305 z + 0.4724 a) stavovou rovnici v normální formě řiditelnosti, b) programové schéma realizujte v SIMULINKu.

Příklad 1.3.2

K diskrétnímu přenosu G ( z ) =

Řešení: a) Stavová diskrétní rovnice je rovna  x1 (k + 1)   0.0 ;1.0   x1 (k )  0 x (k + 1) =   +   ⋅ u (k ); = ⋅  x2 (k + 1) − 0.4724; 1.305  x2 (k ) 1 

 x1 (k )  y (k ) = [0.05716;0.0734]⋅  .  x2 (k )

Programové schéma je na obr.1.3.2-1. Diskrétní přechodová funkce a průběhy x1(k), x2(k) je na obr.1.3.2-2.

x2(k)

y(k)

x1(k)

x2(k)

x1(k)

y(k)

Obr.1.3.2-1 Programové schéma

Konec příkladu

Obr.1.3.2-2.Diskrétní přechodová funkce a průběhy x1, x2 u(t) = 1(t)

1.3.2 Normální forma rekonstruovatelnosti 1) Spojitý popis. Struktura matice koeficientů systému, kterou označujeme jako normální forma rekonstruovatelnosti, vychází z metody postupných integrací. Uvažujeme-li obrazový přenos ve tvaru (1 - 11) , pak diferenciální rovnice má tvar

y(n) + an – 1y(n – 1) + ... + a1y(1) + a0y = bn – 1u(n – 1) + ... + b1u(1) + b0u.

(1 – 23)

Provedeme-li n-násobnou integraci od 0 až t levé a pravé strany této diferenciální rovnice, pak dostaneme integrální rovnici ve tvaru t

t t

t

t

t

t

t

0

0 0

0

0

0

0

0

y = − a n −1 ∫ ydτ − a n − 2 ∫ ∫ ydτ 2 − ⋅ ⋅ ⋅ − a 0 ∫ ⋅ ⋅∫ ydτ n + bn −1 ∫ udτ + ⋅ ⋅ ⋅ + b0 ∫ ⋅ ⋅∫ udτ n (1 – 24)


14

27.8.2004



Programové schéma je na obr.1.5. Jestliže označíme výstup za prvním integrátorem jako první stavovou veličinu x1(t) = y(t), pak můžeme výstup za druhým integrátorem označit jako stavovou veličinu x 2 (t ) atd. Ze schématu na obr.1.5 je zřejmé, že derivace stavových veličin jsou rovny

x&1 (t ) = −a n −1 x1 (t ) + x 2 (t ) + bn −1u (t ) , x& 2 (t ) = −a n − 2 x1 (t ) + x3 (t ) + bn − 2 u (t ) , . . x& n −1 (t ) = − a1 x1 (t ) + x n (t ) + b1u (t ) , x& n (t ) = −a 0 x n (t ) + b0 u (t ) .

Obr.1.5 Programové schéma metody postupné integrace Struktura vektoru stavu x, matice systému A, matice buzení B a matice výstupu C (D = 0) je  x1   − a n−1 x  − a  2   n−2  x3   − a n −3 x= ; A= M   M  xn−1   − a1     − a0  xn 

1 0 L 0 bn−1  b   0 1 L 0  n−2  bn−3  0 0 L 0  ; C = [1,0,0,...,0] ; B= M   b  0 0 L 1  1   0 0 L 0 b0 

(1 – 25)

2) Diskrétní popis. Struktura programovacího schématu je shodná se spojitým popisem, pouze integrátory jsou nahrazeny členy realizující zpoždění z-1 viz obr.1.6.


15

27.8.2004



Obr.1.6 Struktura programového schéma normální formy rekonstruovatelnosti Je zřejmé, že struktura matice systému M, matice buzení N a matice výstupu C bude zachována.  x1 (k )   − a n−1  x (k )  − a   2  n−2  x3 ( k )   − a n −3 x( k ) =   ;M =   M  M  xn−1 (k )  − a1     − a0  xn (k ) 

Příklad 1.3.3

1 0 L 0 bn−1  b  0 1 L 0  n −2  bn−3   0 0 L 0  ; C = [1,0,0,...,0]; ;N =  M   b  0 0 L 1  1   0 0 L 0 b0 

(1 – 26)

0.6 z + 0.15 určete: z 2 − z + 0.25 a) stavovou rovnici v normální formě rekonstruovatelnosti, b) programové schéma realizujte v SIMULINKu. K diskrétnímu přenosu G ( z ) =

Řešení: a) Vektor x(k) a matice M, N, C volíme podle (1 – 26). Stavová rovnice je rovna  x1 (k + 1)   + 1.00 ; 1.0   x1 (k )  0.6  x (k + 1) =  = + ⋅  ⋅ u (k );  x2 (k + 1) − 0.25 ; 0   x2 (k ) 0,15


16

 x1 (k )  y (k ) = [1.0; 0.0]⋅  ;  x2 (k )

27.8.2004



Programové schéma je na obr.1.3.3-1. Diskrétní přechodová funkce a průběhy x1(k), x2(k) jsou na obr.1.3.3-2.

x2(k + 1)

y(k) = x1(k) x2 ( k)

x1(k + 1)

Obr.1.3.3-1 Programové schéma normální formy rekonstruovatelnosti


17

27.8.2004



x1(k) = y(k)

x2(k)

Obr.1.3.3-2 Diskrétní přechodová funkce a průběhy x1(k), x2(k) u(k) = 1(k)

Konec příkladu

1.4 VLASTNOSTI STAVOVÝCH MODELŮ DYNAMICKÉ SOUSTAVY Matematické modely založené na vnějším popise – relace vstup/výstup poskytují pouze informace o výstupním signálu y(t). Pro sledování např. první derivace výstupního signálu je třeba použít derivačního členu (např. derivační člen PID regulátoru). Proto také kriteria jakosti regulace zpravidla zahrnují regulační odchylku a přírůstek akční veličiny. Stavový model popisuje vnitřní stavy systému. Stavový model např. ohřevu vody v průtokovém ohřívači obsahuje, jak jsme si odvodili, dvě stavové veličiny: x1(t) - teplotu topného elementu a x2(t) - teplotu vody. Můžeme si položit otázku: je možno akční veličinou dosáhnout v daném časovém intervalu požadovaných hodnot x1(t1) a x2(t1), které jsou na sobě nezávislé? Odpovědi na tuto otázku a podobné další jsou náplní této kapitoly.

1.4.1 Stavový prostor, vektor stavu, stavová trajektorie Zavedením stavového popisu se převedl matematický model na soustavu n diferenciálních rovnic prvního řádu, kterou obecně zapíšeme do tvaru

[

]

dxi (t ) = f i x1 (t ), x 2 (t ),⋅ ⋅ ⋅, x n (t ), u (t ), d j (t ) , dt

kde

xi(t) u(t) dj(t) fi

(1 – 27)

jsou stavové veličiny, i = 1, 2, …, n, je akční veličina, jsou všechny poruchové veličiny, obecně nelineární funkce.


18

27.8.2004



Pohyb systému je možno zobrazit v n-rozměrném Euklidovském prostoru Εn, kterému říkáme stavový prostor. Souřadnice stavového prostoru tvoří stavové veličiny systému x1, x2, …, xn. Těmito souřadnicemi je určen vektor stavu x(t), který je proměnný v čase a zapíšeme ho ve tvaru x (t ) T = [x1 (t ), x 2 (t ),⋅ ⋅ ⋅, x n (t )] . V každém časovém okamžiku t určují souřadnice x1, x2, …, xn polohu koncového bodu vektoru stavu ve stavovém prostoru. Koncový bod vektoru stavu při pohybu systému opisuje stavovou trajektorii viz.obr.1.7.

x2

Stavová trajektorie x(t1)

x(t 0 ) ... počáteční stav

x(t1 ) ... koncový stav

x1

x(t0)

x3 Obr.1.7 Stavový prostor a stavová trajektorie 1.4.2 Dosažitelnost, řiditelnost, pozorovatelnost, rekonstruovatelnost

Stavový model lineární časově invariantního systému je jednoznačně popsán maticemi A, B, C resp. M, N, C. Pro potřeby analýzy a syntézy je možno definovat vlastnosti systému pomocí stavového modelu. Vlastnosti soustav je možno definovat dosažitelností, řiditelností, pozorovatelností a rekonstruovatelností a dalšími. Odvození těchto vlastností přesahuje rámec kurzu a proto se spokojíme pouze s jejich definicemi a existenčními podmínkami. Dosažitelnost: Stav x(t1) je dosažitelný, jestliže existuje takový vstup u(t) na konečném časovém intervalu t1 - t0 (t0 < t1 ), kterým se převede soustava z počátečního stavu x(t0) = 0 do žádaného stavu x(t1). Podmínky dosažitelnosti:

Libovolný stav lineární stacionární soustavy je dosažitelný tehdy a jen tehdy, jestliže hodnost matice dosažitelnosti D = [B, AB, A2B, … ,An-1B] je rovna rozměru stavového prostoru. hodnost D = hodnost [B, AB, A2B,…, An-1B] = n.

(1 – 28)

Řiditelnost:

Stav x(t0) je řiditelný, jestliže existuje takový vstup u(t) na konečném časovém intervalu t1 - t0 (t1 > t0), kterým se převede soustava z počátečního stavu x(t0) do stavu x(t1) = 0 (do počátku souřadnic).


19

27.8.2004



Podmínky řiditelnosti:

Libovolný stav lineární stacionární soustavy je řiditelný tehdy a jen tehdy, jestliže hodnost matice řiditelnosti R = [A-1B, A-2B, …, A-nB] je rovna rozměru stavového prostoru.

hodnost R = hodnost [A -1B, A-2B, …, A-nB] = n, (1 – 29) hodnost {An [A -1B, A-2B, …, A-nB]} = hodnost [An-1B, …, AB, B] = n (1 – 30)

Hodnost matice (1 – 29) a matice (1 – 30) je stejná, přičemž pracnost je u matice (1 – 30) výrazně nižší. Pozorovatelnost:

Stav x(t0) systému v čase t0 je pozorovatelný, jestliže ho lze určit z průběhu vstupní veličiny u(t) a výstupní veličiny y(t) v konečném časovém intervalu t1 - t0 (t1 > t0). Podmínky pozorovatelnosti:

Lineární stacionární soustava je pozorovatelná tehdy a jen tehdy, jestliže hodnost matice pozorovatelnosti P = [CT, ATCT, A2TCT, …, (AT)n-1CT] je rovna řádu systému. hodnost P =hodnost [CT, ATCT, A2TCT, …, (AT)n-1CT] = n.

(1 – 31)

Rekonstruovatelnost:

Stav x(t0) systému v čase t0 je rekonstruovatelný, jestliže ho lze určit z průběhu vstupní veličiny u(t) a výstupní veličiny y(t) v konečném časovém intervalu t0 - t1 (t0 > t1, tj. z hodnot minulých). Podmínky rekonstruovatelnosti:

Lineární stacionární soustava je rekonstruovatelná tehdy a jen tehdy, je-li det A ≠ 0 a je-li hodnost matice rekonstruovatelnosti RE = [(AT)-nCT, (AT)(-n+1)CT, …, CT] je rovna řádu systému. hodnost [(AT)-nCT, (AT)(-n+1)CT,…, CT] = n. (1 – 32) T T -n T T (-n+1) T T T T T n/1 T hodnost {A [(A ) C , (A ) C , …, C ]} = hodnost [C , A C ,…, A C ] = n. (1 – 33)

Opět jsou hodnosti matice rekonstruovatelnosti (1 – 32) a matice pozorovatelnosti (1 – 33) stejné. Pro praktické ověření rekonstruovatelnosti u konkrétního modelu použijeme pro výpočet rovnost (1 – 33). Příklad 1.4.1

Rozhodněte o dosažitelnosti a rekonstruovatelnosti systému, jehož stavový model je

 x (k + 1)   0.0 ;1.0   x1 (k )  0 x (k + 1) =  1  +   ⋅ u (k ); = ⋅  x2 (k + 1) − 0.4724; 1.305  x2 (k ) 1 


20

 x (k )  y (k ) = [0.05716;0.0734]⋅  1 .  x2 ( k ) 

27.8.2004



Řešení: Sestavíme matici dosažitelnosti pro diskrétní systém D = [N, MN] v MATLABu viz. Obr.1.4.1a, její hodnost je 2, systém je řiditelný.

Obr.1.4.1a Matice dosažitelnosti

Obr.1.4.1b Modifikovaná matice rekonstruovatelnosti

Matice rekonstruovatelnosti dle (1 – 33) a její hodnost je na obr.1.4.1b . Konec příkladu

1.5 TRANSFORMACE STAVOVÝCH ROVNIC Stavový vektor je možno transformovat z jedněch souřadnic do druhých pomocí transformačních matic T. Uvažujme stavové rovnice ve tvaru (1 – 1) x& = Ax (t ) + Bu (t ) ; rovnice výstupu y = Cx(t) + Du(t). Nechť platí

x = Tz

(1 – 34)

pak dosazením (1 – 34) do rovnice (1 – 1,3) dostaneme Tz& = ATz + Bu (t ) , y = CTz + Du(t).

(1 – 35)

Vynásobíme-li rovnici (1 – 35) zleva inverzní maticí T-1, dostaneme stavový popis v nových souřadnicích z& = T −1 ATz + T −1Bu (t ) (1 – 36) y = CTz + Du(t). (1 – 37) kde je AT = T-1AT transformovaná matice soustavy a BT = T-1B transformovaná matice buzení.


21

27.8.2004



Transformační matice T je čtvercová matice [n, n] a musí býti invertovatelná. Je zřejmé, že odvozená transformační pravidla platí v celém rozsahu i pro diskrétní stavový popis.

1.6 MODEL PORUCH VE STAVOVÉM POPISU Z vnějšího popisu dynamických soustav víme, že poruchové veličiny dělíme na deterministické a stochastické. Deterministické poruchy pak dělíme na měřené a neměřené. Struktura modelu, která je založena na vnějším popise a aproximuje pomocí přenosů vliv poruch, je na obr. 1.7. Y(s) = Fu(s)U(s) + Fd(s)D(s) + Fm(s)Dm(s) (1 – 38) D(s) Fd(s) Yd(s) Fu(s) … aproximuje dynamické účinky akční veličiny u(t) Y Dm(s) m(s) F (s) … aproximuje dynamické účinky d Fm(s) neměřené poruchy d(t) F (s) … aproximuje dynamické účinky m Yu(s) Y(s) U(s) měřené poruchy dm(t) Fu(s) Obr.1.7 Vnější popis s modelem poruch 1) Ve stavovém popisu je dynamický účinek determinovaných poruch (Fm(s) = 0) popsán rovnicí x& (t ) = Ax(t ) + B U u (t ) + B d d (t ), x(k + 1) = Mx(k) + NUu(k) + Ndd(k), (1 – 39,40) y(t) = Cx(t) + Du(t), kde jsou BU, NU Bd, Nd

y(k) = Cx(k) + Du(k),

matice buzení akční veličiny matice buzení poruchové veličiny.

Struktura stavového modelu s neměř. poruchou d(t) a s maticí D = 0 je na obr.1.8a.

Obr.1.8a Stavový model s neměřenou poruchou 2) Účinek stochastických poruch vyjádříme vektorem stochastických signálů ω(t), ω(k), který se přičte ke stavovému vektoru. K výstupu y(t), y(k) se přičítá náhodný signál v(t), v(k). Stochastický signál v a vektorový stochastický signál ω jsou charakterizovány


22

27.8.2004



{ }

{

}

E{ν } = 0, E{ω} = 0 , E νν T = R, E ωωT = C ωω . Stavové rovnice pak mají tvar x& (t ) = Ax(t ) + BU u (t ) + ω(t ),

x(k + 1) = Mx(k) + NUu(k) + ω(k), (1 – 41, 42)

y(t) = Cx(t) + ν(t),

y(k) = Cx(k) + ν(k).

(1 – 43, 44)

Struktura diskrétního stavového modelu s účinky stochastických poruch je na obr.1.8b.

ν(k)

ω(k ) x (k + 1) N u(k)

y(k)

x(k )

1 z

C

M Obr.1.8b Struktura diskrétního stavového modelu s účinky stochastických poruch 3) Ukažme ještě odvození modifikované struktury normální formy rekonstruovatelnosti pro neměřenou poruchu d(t). Uvažujeme-li obraz Y(s) ve tvaru (1 – 38) a dosazením za FU, Fd dostaneme B( s) C ( s) Y ( s ) = FU ( s )U ( s ) + Fd ( s ) D( s ) = U ( s) + D( s) . A( s ) A( s ) Diferenciální rovnice má pak tvar y ( n ) + a n −1 y ( n −1) + a n − 2 y ( n − 2 ) + ⋅ ⋅ ⋅ + a1 y (1) + a 0 y = = bn −1u ( n −1) + ⋅ ⋅ ⋅ + b1u (1) + b0 u + c0 d + c1 d (1) + ⋅ ⋅ ⋅ + c n −1 d ( n −1) Provedeme-li n-násobnou integraci od t = 0 až t = ∞ levé a pravé strany této diferenciální rovnice, pak dostaneme integrální rovnici ve tvaru t

t t

t

t

0

0

y = − a n −1 ∫ ydτ − a n − 2 ∫ ∫ ydτ 2 − ⋅ ⋅ ⋅ − a 0 ∫ ⋅ ⋅∫ ydτ n + 0

0 0

t

t

t

t

t t

t

t

0

0

0

0

0 0

0

0

+ bn −1 ∫ udτ + ⋅ ⋅ ⋅ + b0 ∫ ⋅ ⋅∫ udτ n + c n −1 ∫ d (τ )dτ + c n − 2 ∫ ∫ d (τ )dτ 2 + ⋅ ⋅ ⋅ + c0 ∫ ⋅ ⋅∫ d (τ )dτ n

Programové schéma je na obr.1.9. Zvolíme-li první stavovou veličinu x1(t) = y(t), pak ze schématu je zřejmé, že derivace stavových veličin jsou rovny


23

27.8.2004



x&1 (t ) = −a n −1 x1 (t ) + x 2 (t ) + bn −1u (t ) + c n −1 d (t ) x& 2 (t ) = −a n −2 x1 (t ) + x3 (t ) + bn − 2 u (t ) + c n − 2 d (t ) . . . x& n −1 (t ) = − a1 x1 (t ) + x n (t ) + b1u (t ) + c1 d (t ) x& n (t ) = −a0 x1 (t ) + b0 u (t ) + c0 d (t ) .

(1 – 41)

Obr.1.9 Programové schéma s neměřenou poruchou d(t) Vektor stavu x, matice systému A, matice buzení BU, Bd a matice výstupu C mají strukturu (D = 0) C = [1, 0, 0, ..., 0]  x1   − a n−1 x  − a  2   n−2  x3   − a n −3 x= ; A= M   M  xn−1   − a1     − a0  xn 

Příklad 1.9

1 0 L 0 bn−1  cn−1  b  c   0 1 L 0  n−2   n−2  bn−3   c n −3  0 0 L 0  ; Bd =  .  ; BU =  M M      b  c  0 0 L 1  1   1   0 0 L 0 c0  b0 

Dynamický systém je popsán obrazovými přenosy 3s + 1 0,75 . ; Fd ( s ) = 2 s + 1,5s + 0,5 s + 1,5s + 0,5 1) Stavový model s poruchami v normální formě rekonstruovatelnosti. 2) Ověřte chování obvodu při vstupu šumového signálu na vektor stavu. FU ( s ) =

Vytvořte:

(1 – 42)

2


24

27.8.2004



Řešení: 1) Podle (1 – 42) jsou matice rovny − 1,5; 1  0,00 3 ; BU =  A= ; B d =  ;   − 0,5; 0 0,75 1 

C = [1, 0]

Programové schéma je na obr.1.9.1.

x1(t) = y(t)

x2(t)

Obr.1.9.1 Programové schéma stavového modelu s determinovanou poruchou d(t) Výsledné průběhy odezev stavových veličin na vstupní signál u(t) = 1(t); d(t) = 0 a u(t) = 0; d(t) = 1(t) jsou na obr.1.9.2a, b. 2) Vstupní šumové signály jsou na obr.1.9.3a, stavové veličiny jsou na obr.1.9.3b.

x1(t) = y(t)

x2(t) x1(t) = y(t) x2(t)

Obr.1.9.2a u(t) = 1(t - 1), d(t) = 0

Obr.1.9.2b u(t) = 0, d (t) = 1(t)

x2(t)

x2(t)

x1(t) x1(t)

Obr.1.9.3a Vstupní šumový signál ω1(k), ω2(k) Obr.1.9.3b Stavové veličiny s účinkem šumu


25

27.8.2004



2. ESTIMACE STAVU Matematický model reálného procesu ve formě stavových rovnic obsahuje stavový vektor. Informace o stavových veličinách v daném časovém okamžiku tohoto modelu získáme simulací nebo řešením stavových rovnic. Informace o reálném procesu získáváme měřením, pozorováním atd. Aby bylo možno provést měření zvolené veličiny, musí býti známé místo měření, což již samo o sobě může představovat komplikovaný úkol. Dále je nutno zabudovat čidlo nebo senzor, který musí býti k disposici a realizovat celou řadu podpůrných konstrukcí, prvků a pomocných obvodů. Musíme si uvědomit, že pro měření některých veličin ani vhodná čidla neexistují. Jedná se o vyšší derivace signálů apod. Náklady spojené s měřením nejsou také zanedbatelné. To jsou všechno důvody, proč zpravidla není možné všechny stavové veličiny měřit i když jsou pro stavovou regulaci nezbytné. Cílem této kapitoly je vysvětlit základní metody a postupy, které umožňují provádět výpočet vektorů stavových veličin pomocí stavových modelů a na základě měření vstupů u(t) resp. u(kT) a výstupu y(t) resp. y(kT). Výpočet stavového vektoru za uvedených podmínek se označuje jako odhad – estimace stavového vektoru. Blokové schéma prvků, pomocí kterých se odhad realizuje, se nazývá estimátor. Podle toho, zda měřený výstup je nebo není zatížen aditivním parazitním šumem, můžeme konstatovat: 1) Jestliže měřený signál y(t) resp. y(kT) neobsahuje aditivní šumový signál, pak je možno předpokládat, že na stavové veličiny modelu nepůsobí šumové signály. Estimaci stavu za této podmínky nazýváme deterministickou estimací stavu. 2) Jestliže měřený signál y(t) resp. y(kT) obsahuje aditivní šumový signál, pak jeho účinek v modelu je třeba aproximovat působením vhodných šumových signálů na stavových veličinách přes matici buzení. Estimaci stavu pak nazýváme Kalmanovou estimací stavu. Vzhledem k softwarové podpoře syntézy ve stavovém prostoru v MATLABu se omezíme na diskrétní estimátory stavu.

2.1 DETERMINISTICKÝ ESTIMÁTOR STAVU Návrh estimátoru, který popíšeme, byl navržen D. G. Luenbergrem v /5/ a vychází z předpokladu, že je znám lineární stavový model a matice M, N, C. Dále se předpokládá, že je měřen výstup ze soustavy y(k). Stavové rovnice jsou tvaru x(k + 1) = Mx(k) + Nu(k) y(k) = Cx(k) Je logické pro odhad stavu využít v prvním přiblížení přímo struktury modelu jako estimátoru viz obr.2.1. u(k)

x(k )

Dynamický

C

y(k)

systém xˆ (k )

yˆ (k )

C

Estimátor

Obr.2.1 Estimátor bez korekce z výstupu dynamického systému


26

27.8.2004



Estimátor pak bude popsán rovnicí xˆ (k + 1) = M E xˆ (k ) + N E u (k ) , yˆ (k ) = Cxˆ (k ) ,

kde jsou

ME, NE …matice estimátoru, yˆ (k ), xˆ (k ) …odhadovaný výstup a vektor stavu.

Zavedeme chybu estimace ∆x (k + 1) = x (k + 1) − xˆ (k + 1) = Mx (k ) + Nu (k ) − [M E xˆ (k ) + N E u (k )] = = Mx(k ) − M E x(k ) + (N − N E )u (k ).

Aby chyba a tedy ani odhad nezávisely na akční veličině u(k), musí býti splněna rovnost N = NE. Podobně nechceme, aby chyba estimace závisela na stavovém vektoru x(k) a jeho odhadu xˆ (k ), proto položíme M = ME. Chyba estimace je pak rovna ∆x (k + 1) = M∆x (k ) .

(2 – 1)

Je zřejmé, že jestliže ∆x(0) = x(0) − xˆ (0) ≠ 0 , pak estimace je podle (2 – 1) dynamický proces, který závisí na matici systému M. Tento estimátor nevyhovuje běžným požadavkům. Je proto třeba dalších informaci, pomocí kterých pak můžeme provést korekci estimace. Přirozená informace o kvalitě odhadu je obsažena v rozdílu ∆y (k ) = y (k ) − yˆ (k ) , který použijeme pro požadovanou korekci. Zavedeme-li tento rozdíl do estimátoru pomocí matice L, pak rovnice estimátoru bude xˆ (k + 1) = M E xˆ (k ) + N E u (k ) + L[ y (k ) − yˆ (k )] .

(2 – 2)

Vyjádříme-li odhad výstupu estimátoru pomocí vektoru odhadu xˆ (k ) , dostaneme xˆ (k + 1) = M E xˆ (k ) + N E u (k ) + L[ y (k ) − Cxˆ (k )] .

Chyba odhadu je pak rovna ∆x(k + 1) = x(k + 1) − xˆ (k + 1) = Mx(k ) + Nu (k ) − {M E xˆ (k ) + N E u (k ) + L[ y (k ) − Cxˆ (k )]}

Provedeme-li následující úpravy ∆x (k + 1) = Mx (k ) − Mxˆ (k ) + Mxˆ (k ) − M E xˆ (k ) + ( N − N E )u (k ) − Ly (k ) + LCxˆ (k ) = = M ⋅ ∆x (k ) + ( M − M E ) xˆ (k ) + ( N − N E )u (k ) − LCx (k ) + LCxˆ (k )

a zavedeme-li rovnosti N = NE, M = ME, dostaneme pro chybu odhadu rovnici Doc. Ing. Modrlák Osvald, CSc.

27

27.8.2004



∆x(k + 1) = (M – LC) ∆x(k) ,

(2 – 3)

kde L je zatím neurčená matice estimátoru. Je zřejmé, že dynamika chyby odhadu závisí na volbě právě této matice a volí se tak, aby chyba odhadu konvergovala k nule. Návrh estimátoru ukážeme na příkladě. Struktura estimace stavového vektoru je na obr.2.2a,b. u (t )

Dynamický systém

y (t )

u (k )

y (k )

Dynamický systém

b)

Estimátor stavu

L xˆ (k )

xˆ (k + 1) xˆ (t )

a)

z −1

N

C

C

yˆ (k )

M yˆ (t )

xˆ (k + 1) = M E xˆ (k ) + N E u (k ) + L[ y (k ) − Cxˆ (k )]

Obr. 2.2 Struktura estimace odhadu vektoru stavu Dynamické vlastnosti chyby odhadu jsou dány vlastními čísly matice (M – LC). Vlastní čísla je možno určit z determinantu (2 – 4). Je zřejmé, že koeficienty charakteristického polynomu a tedy i vlastní čísla matice (M – LC) závisí na matici estimátoru L. Výpočet det{zI − (M − LC)} = z n + a n −1 z n −1 + a n − 2 z n − 2 + ⋅ ⋅ ⋅ + a1 z + a 0

(2 – 4)

determinantu se zjednoduší, jestliže matice M má speciální strukturu, např. normální formu rekonstruovatelnosti. Pak determinant uzavřeného obvodu je roven det{zI − ( M − LC )} = z n + (l1 + a n −1 ) z n −1 + (l 2 + a n − 2 ) z n − 2 + ⋅ ⋅ ⋅ + (l n −1 + a1 ) z + (l n + a0 ) Celý postup návrhu estimátoru ukážeme na následujícím příkladě.

(2 – 4a)

Příklad 2.1 Pro dynamickou soustavu popsanou stavovou rovnicí (1) a (2) nalezněte: 1) Estimátor stavu. Matici estimátoru volte tak, aby vlastní čísla charakteristického polynomu estimátoru byla rovna nule. 2) Určete chybu estimace pro k = 1, 2, 3, jestliže počáteční vektor stavu soustavy je  x (k + 1)   1,00 ; 1   x1 (k )   0,6   1  x (0) =  , x(k + 1) =  1 = +  ⋅  ⋅ u (k ), − 0,5  x 2 (k + 1) − 0,25 ; 0  x 2 (k ) 0,15


28

(1)

27.8.2004



 x1 (k )   1,0 ;1   0,6  . M= ;N = y (k ) = [1 ;0] ⋅    ; C = [ 1 , 0 ] . − 0,25 ; 0 0,15  x 2 (k )

(2)

Řešení: 1) Estimátor stavu. Rovnice estimátoru je xˆ (k + 1) = M E xˆ (k ) + N E u (k ) + L[ y (k ) − Cxˆ (k )] , ∆x(k + 1) = (M – LC) ∆x(k). Je nutno vyjádřit matice (M - LC)  1 − l1 E ;1,0   1,0; 1  l1 E  [ ] 1 , 0 − ⋅ = ( M − LC ) =  − (0,25 + l ) ;0    − 0,25;0 l 2 E  2E   Determinant soustavy je  z , 0   1 − l1E ; 1,0    z + (l1E − 1) ; − 1   −  AE ( z ) = det{zI − (M − LC )} = det   = det   =  ( l 2 E + 0,25) ; z   0 , z  − (0,25 + l 2 E ) ;0 

= z[z + (l1E – 1)] + (l2E + 0,25) = z2 + z(l1E – 1) + (l2E + 0,25). Jestliže položíme koeficienty matice estimátoru l1E = +1,0; l 2 E = −0,25

0 , 1 (M − LC) =  0 , 0

 , 

pak charakteristický polynom estimátoru je roven AE(z) = z2. Vlastní čísla z1 = z2 = 0 a chyba odhadu bude nejrychleji konvergovat k nule.

x1(k)

x2(k)

Obr.2.1.1. Průběh estimace stavu při buzení u(k) = 1(k) a počátečním stavu na soustavě x(0).


29

27.8.2004



2) Chyba estimace pro k = 1, 2, 3.  1  Chyba estimace pro k = 0 je rovna ∆x(0) = x(0) =  . − 0,5 Dynamiku chyby vyjádříme rovnicí 0; 1  ∆x1 (k )  ∆x(k + 1) = (M − LC)∆x(k ) =   ⋅  . 0; 0 ∆x2 (k ) Chyba estimace pro k = 1, 2, 3 je rovna 0; 1  ∆x1 (0)  0; 1   1  − 0,5 ∆x(1) =   ⋅   =  ⋅ =  0; 0 ∆x2 (0) 0; 0 − 0,5  0 

0  0; 1  ∆x1 (1)  0; 1  − 0,5 0 ∆x ( 2 ) =   ⋅  =  ⋅ =   ; ∆x(3) =   .   0  0; 0 ∆x2 (1) 0; 0  0  0 Průběh estimace při buzení u(k) = 1(k) a počátečním stavu na soustavě x(0) je na obr.2.1.1.

Estimátor

Obr.2.1.2 Programové schéma estimace v SIMULINKU


30

Konec příkladu

27.8.2004



2.2 KALMANOVA FILTRACE Jestliže výstupní signál je zatížen šumovým parazitním signálem, je třeba uvažovat ve stavovém modelu působení šumového signálu nejen na výstupní veličinu y(t), ale i na jednotlivé složky stavového vektoru x(t). Stavová rovnice pak bude míti tvar

x(k + 1) = Mx(k) + Nu(k) + ω(k),

(2 – 5)

y(k) = Cx(k) + ν(k),

(2 – 6)

kde ν(k), ω(k) je posloupnost náhodných nezávislých vektorů. Strukturu estimátoru ponecháme, ovšem matici estimátoru bude nutno určit na základě jiných požadavků a podmínek. xˆ (k + 1) = Mxˆ (k ) + Nu (k ) + L[ y (k ) − Cxˆ (k )]

Chyba odhadu bude ∆x (k + 1) = x (k + 1) − xˆ (k + 1) = Mx (k ) + Nu (k ) + ω(k ) − {Mxˆ (k ) + Nu (k ) + L[ y (k ) − Cxˆ (k )]}.

Jestliže vyjádříme y(k) = Cx(k) + ν(k), pak po jednoduché úpravě je možno chybu odhadu zapsat ve tvaru ∆x(k + 1) = (M – LC) ∆x(k) + ω(k) – Lν(k). (2 – 7) Co tvoří vstupy diferenční rovnice chyby odhadu? Je to chyba odhadu v čase k = 0, ∆x(0) = x(0) − xˆ (0) ≠ 0 a stochastické vektory. Chceme-li tedy analyzovat chybu odhadu, je třeba si uvědomit, že vektor reprezentující tuto chybu, ω(k) a ν(k) představuje vektorový náhodný proces. Tento náhodný proces, jak je z (2 – 7) zřejmé, je popsán maticovou diferenční rovnicí. Parametry M, C maticové diferenční rovnice jsou dány. Matici estimátoru L je třeba určit. V rámci našeho kurzu jsme se již jednou s podobnou analýzou setkali. Provedli jsme rozbor vlastností vektoru chyb odhadu parametrů ∆ = pˆ − p . Základní informace o chybě odhadu poskytuje kovariační matice, která na hlavní diagonále obsahuje informace o rozptylech. Kovariační matici chyby odhadu v čase k je rovna

{

}

P(k ) = Ε [∆x(k ) − Ε{∆x(k )}] ⋅ [∆x(k ) − Ε{∆x(k )}] , T

kde E{ } je operátor střední hodnoty. Protože systém je lineární a podle předpokladu jsou stochastické vektory ω(k) a ν(k) stacionární s nulovou střední hodnotou, musí střední hodnota chyby odhadu býti rovna E{∆x(k)} = 0. Takže kovariační matice bude rovna P(k) = E{∆x(k) ∆x(k)T}.

(2 – 8)

Kovariační matici můžeme vyjádřit také pro časový okamžik k + 1 P(k + 1) = E{∆x(k + 1) ∆x(k + 1)T}.


31

27.8.2004



Dosadíme-li za ∆x(k+1) z rovnice (2 – 7) a provedeme-li transposici hranaté závorky, dostaneme P(k + 1) = Ε [(M − LC)∆x(k ) + ω(k ) − Lv(k )] ⋅ ∆x(k ) T (M − LC) T + ω(k ) T − v(k ) T LT

[

{

]}

Aplikujeme-li operátor střední hodnoty na jednotlivé členy, je možno vzhledem k předpokladům zapsat

{

}

{

}

{

}

Ε ω(k ) ∆x (k ) T = 0; Ε ν (k ) ∆x (k ) T = 0; Ε ω(k )ν (k ) T = 0 a kovariační matice je rovna

P (k + 1) = ( M − LC )Ε{∆x (k )∆x (k ) T }( M − LC ) T + Ε{ω(k )ω(k ) T }+ LΕ{v(k )v(k ) T }LT .

{

}

{

}

{

}

Zavedeme-li označení P(k ) = Ε ∆x(k )∆x(k ) T ; Q = Ε ω(k )ω(k ) T ; R = Ε v(k )v(k ) T , pak kovariační matice chyby odhadu bude mít tvar

P(k + 1) = (M - LC)⋅P(k)⋅(M - LC)T + Q + LRLT, kde

Q R

(2 – 9)

je kovariační matice [n, n] vektoru šumu ω (k), je kovariační matice [1, 1] vektoru šumu ν(k).

Označíme-li dále symbolem R0 kovariační matici chyb prvního odhadu stavového vekP(0) = R0. toru ∆xˆ (0), pak podle (2 – 8) je Kovariační matice chyby odhadu podle (2 – 9) se dynamicky v čase vyvíjí. Jediným zatím neurčeným parametrem je matice estimátoru L. Pomocí této matice je možno minimalizovat středně kvadratickou chybu. Středně kvadratickou chybu skalárního součinu ∆ = aTx je možno vyjádřit ve tvaru

{ }

{(

Ε ∆2 = Ε a T ∆x

) }= Ε{(a 2

T

)(

)}

{

}

∆x ∆x T a = a T Ε ∆x∆x T a = a T P(k )a .

Naším úkolem nyní bude navrhnout matici estimátoru L tak, aby středně kvadratická chyba ∆2 byla minimální. Rovnost (2 – 9) je kvadratická forma vzhledem k matici estimátoru L a její minimalizaci ve smyslu středně kvadratické chyby provedeme rozkladem na úplné čtverce. Rovnost (2 – 9) upravíme do tvaru

[

]

P(k + 1) = L C ⋅ P(k ) ⋅ CT + R ⋅ LT − 2M ⋅ P(k ) ⋅ CT LT + MP (k )M T + Q . (2 – 10) Zavedeme-li substituci A = C ⋅ P(k ) ⋅ C T + R; h = M ⋅ P(k ) ⋅ C T ; α = MP (k)M T + Q ; x = L, je možno kovariační matici zapsat do tvaru P(k + 1) = xAxT − 2hx T + α Minimalizaci provedeme rozkladem na úplné čtverce viz Příloha 1 a dle (P – 8,9) dostaneme Doc. Ing. Modrlák Osvald, CSc.

32

27.8.2004



[

x opt = LEopt (k ) = hA −1 = M ⋅ P (k ) ⋅ C T C ⋅ P (k ) ⋅ C T + R

]

−1

(2 – 11)

Dosazením (2 – 11) do (2 –10) dostaneme

[

P(k + 1) = M ⋅ P(k ) ⋅ C T C ⋅ P(k ) ⋅ C T + R

] ⋅ (M ⋅ P(k ) ⋅ C ) −1

T T

+ MP (k )M T + Q

a úpravou dostaneme P (k + 1) = ( M − L(k )C ) ⋅ P (k ) ⋅ ( M − L(k )C ) T + Q + L(k ) RL(k ) T .

(2 – 12)

Rovnice (2 –12) a (2 – 11) dávají rekurzivní vztah pro výpočet kovariační matice P(k+1) a matice estimátoru L(k) v jednotlivých krocích z výchozího stavu P(0). Pro neomezený počet kroků lim L(k ) = L . k →∞

Protože matice estimátoru L(k) nezávisí na měřených hodnotách, je možno matici estimátoru spočítat před vlastním odhadem stavu. Struktura odhadu vektoru s Kalmanovým filtrem je na obr.2.3. y (k )

Kalmanův filtr

u (k )

xˆ (k )

Soustava

y (k )

ω( k ) v(k )

Obr.2.3 Struktura odhadu stavu s Kalmanovým estimátorem Příklad 2.2

Pro dynamickou soustavu popsanou maticemi D = 0,

nalezněte: 1) matici Kalmanova filtru a porovnejte s maticí deterministického estimátoru, 2) Kalmanův filtr realizujte v SIMULNKu, 3) zobrazte průběh estimace pro počáteční chybu odhadu ∆x(0) = [1; 1; 1]T. Řešení: 1) Diskrétní stavový model je definován maticemi MD, ND, CD, D = 0. Matice deterministického estimátoru splňující podmínku det(zI – (MD – LD1⋅CD)) = z3 podle (2 – 4a) je na vektoru LD1. Výpočet Kalmanova filtru pomocí funkce kalman (viz kap.4.2 Př.4.5) je na obr.2.4. Matice Q, R byly zvoleny rovny jedné.


33

27.8.2004



Obr.2.4 Příkazy a výsledky výpočtu matice Kalmanova filtru. 2) Realizační schéma estimátoru je na obr.2.5. Obsahuje přenos spojité soustavy, diskrétní model v bloku Discrete State-Space, jejímž výstupem je diskrétní stavový vektor soustavy1. Matici estimátoru LD1 realizuje MatrixGain 4, matice MD, ND, CD, C jsou implementovány v Matrix Gain 2, 5, 3, 6. Vstupní signál generuje blok Step.

1

Pro získání všech stavových proměnných je v Discrete State-Space parametr C jednotkovou maticí. Takto budeme postupovat i v dalších schématech, kde bude tento blok použit. Pro získání skutečného skalárního výstupu soustavy použijeme Matrix Gain, jehož náplní je vektor C ze stavového popisu. Dále do parametru A v Discrete State-Space vepíšeme MD a místo B použijeme ND.


34

27.8.2004



Obr.2.5 Realizační schéma estimátoru 3) Průběh estimace s počáteční chybou odhadu ∆x(0) = [1; 1; 1]T a skoku u(k) = 1(k – 15) pro Kalmanův estimátor je na obr.2.5.1,2.

yˆ (k ) y( k)

Obr.2.5.1 Průběh y(k) a yˆ (k )

Obr.2.5.2 Průběh x(k ), xˆ (k )

Průběh estimace s deterministickým estimátorem je na obr.2.5.3 a 2.5.4. y( k)

yˆ (k )

Obr.2.5.3 Průběh y(k) a yˆ (k ) determ. est.

Obr.2.5.4 Průběh x(k ), xˆ (k ) determ. est. Konec příkladu


35

27.8.2004



PŘÍLOHA P1 MINIMALIZACE KVADRATICKÉ FUNKCE Předpokládejme, že je dána kvadratická funkce y (x) = x T Ax + 2h T x + α ,

(P – 1)

kde h, x jsou vektory a A je obecně nesymetrická regulární matice, ∝ je konstanta. Minimalizace kvadratické funkce spočívá v nalezení takového x, aby y → MIN. Nalezení minima kvadratické funkce provedeme pro jednoduchost technikou převedení (P – 1) na úplné čtverce. Vytvoření úplných čtverců vysvětlíme na úloze hledání minima funkce jedné proměnné. Můžete se přesvědčit, že platí y = x 2 + 2hx + c = ( x + h ) ⋅ ( x + h ) − h 2 + c

(P – 2)

Máme-li najít minimum funkce y v závislosti na x, pak je zřejmé, že funkce y bude minimální, jestliže x = -h. (P – 3) min y = −h 2 + c, pro xopt = − h. x

Podobně lze odvodit h  h  h2  y = ax + 2hx + c =  x +  ⋅ a ⋅  x +  − + c a  a a  2

min y = − x

h2 + c, a

pro

xopt = −

h a

(P – 4)

(P – 5)

nebo též h  h  h2  y = (ax) 2 + 2hx + c. =  ax +  ⋅  ax +  − 2 + c, a  a a  2 h h min y = − 2 + c, pro x opt = − 2 . x a a

(P – 6) (P – 7)

Pro minimalizaci kvadratické funkce (P – 1) použijeme rovnosti (P – 4,5). Úplný čtverec kvadratické funkce má tvar y (x) = x T Ax + 2h T x + α = (x + A −1h) T A(x + A −1h) − h T A −1 ⋅ h + α

(P – 8)

přičemž minimum y a x vyhovující podmínce minima je rovno

Konec přílohy P1

min y = −h T A −1h + α ; x opt = − A −1h. x


36

(P – 9)

27.8.2004



3. SYNTÉZA STAVOVÝCH REGULÁTORŮ Základní rozdíl mezi klasickým PID a stavovým regulátorem spočívá v tom, že stavový regulátor využívá pro výpočet akční veličiny všechny stavové veličiny. Platí uR(t) = -Kx(t);

kde je K x(t), x(k)

uR(k) = -Kx(k),

(3 – 1,2)

matice regulátoru [1; n], vektory stavu (spojitý nebo diskrétní).

Hovoříme pak, že využívá úplné informace o systému. Pokud nejsou všechny stavové veličiny přístupné měřením, pak se tyto stavové veličiny odhadují pomocí estimátoru a místo vektorů x(t), x(k) jsou dosazovány vektory estimované xˆ (t ), xˆ (k ) . Při návrhu matice regulátoru a estimátoru platí princip separability, což znamená, že návrhy se vzájemně neovlivňují.

3.1 STRUKTURA STAVOVÝCH REGULÁTORŮ Struktura obvodu se stavovým regulátorem je na obr.3.1. Vychází z předpokladu, že všechny stavové veličiny jsou měřitelné na regulované soustavě. Není-li možno stavové

Obr.3.1 Soustava se stavovým regulátorem

Obr.3.2 Soustava se stavovým regulátorem a estimátorem stavu veličiny měřit, je nutno je odhadovat a použít estimátor stavu viz obr.3.2. Takto navržený a zapojený stavový regulátor převede sice pro nenulové počáteční podmínky systém do počát-


37

27.8.2004



ku stavových souřadnic, ale nedokáže odstranit účinek trvalých poruch, např. ve tvaru skokových změn poruchových veličin, ani zajistit změny žádaných hodnot. Aby toho bylo možno dosáhnout, je třeba zařadit do schématu integrátor viz obr 3.3.

Obr.3.3 Soustava se stavovým regulátorem, estimátorem stavu a integrátorem Výstup z integrátoru je další stavová veličina, takže počet stavových veličin se zvýší na (n + 1). Je proto přirozené, že se rozšíří také matice regulátoru. Výstup z regulátoru je pak roven uR(t) = -Kn + 1xn + 1(t); uR(k) = -Kn + 1xn + 1(k) (3 – 3,4) kde je Kn+1 xn + 1(t), xn + 1(k)

matice regulátoru [1; n + 1], vektory stavu [(n + 1), 1] (spojitý nebo diskrétní).

Schémata na obr.3.1, 3.2 a 3.3 platí také pro diskrétní popis s tím, že integrátor se nahradí blokem jednotkového zpoždění.

3.2 NÁVRH STAVOVÉHO REGULÁTORU METODOU ZADÁNÍ PÓLŮ Metoda návrhu stavového regulátoru zadáváním pólů charakteristické rovnice patří při výkladu tématu syntézy stavového regulátoru z hlediska pedagogického přístupu i náročnosti k základům výkladu. Používá se v případech, kdy jsou póly charakteristického polynomu zadány, nebo je požadován charakteristický polynom s předem danými koeficienty. Je třeba si uvědomit, že koeficienty takto navrženého stavového regulátoru nerespektují tvar poruchy ani místo jejího vstupu, tedy ani dynamické vlastnosti poruch.


38

27.8.2004



3.2.1 Návrh stavového regulátoru Vlastní návrh stavového regulátoru bude proveden pro obvod dle obr.3.1 a vychází z těchto předpokladů 1) Obvod je vychýlen z rovnovážného stavu x(0) ≠ 0 2) Všechny stavové veličiny jsou měřitelné 3) Stavový popis bude uvažován v normální formě řiditelnosti Chování uzavřeného spojitého i diskrétního obvodu můžeme popsat následujícím způsobem Spojitý popis

Diskrétní popis

x& (t ) = Ax (t ) + Bu (t )

x(k + 1) = Mx(k) + Nu(k) Dosazením za u(t) z (3 – 1,2)

x& (t ) = Ax (t ) + B[− Kx (t ) + dU (t )] =

kde

x(k + 1) = Mx(k) + N[-Kx(k) + dU(k)] =

= (A – BK)x(t) + B⋅ dU(t) =

= (M – NK)x(k) + N⋅ dU(k) =

= AK x(t) + B⋅ dU(t),

= MK x(k) + N⋅ dU(k),

AK = (A – BK), MK = (M – NK)

(3 – 5,6)

jsou matice koeficientů uzavřeného obvodu.

Matice soustavy A, M a vektory buzení B, N v normální formě řiditelnosti jsou  0  0   0 A=  M  0  − a0

1

0

L

0

1

L

0

0

L

0

0

L

− a1

− a2 L

0   0 0    0 0  0      0 0  0  ;B =   ; M =    M M   0 0  1      − an−1  1  − a0

1

0

L

0

1

L

0

0

L

0

0

L

− a1

− a2 L

0  0   0  0    0  0  ;N =    M  0  1     − an−1  1 

Charakteristický polynom uzavřeného obvodu je roven determinantu det[sI – AK] = det [sI – (A – BK)],

det[sI – MK] = det [sI – (M – NK)].

Jestliže matice stavového regulátoru je

K = [k1; k2; k3; ...; kn - 1; kn], pak matice (A – BK), (M – NK) je rovna


39

27.8.2004



0 ; 1 ; 0 ; ⋅⋅⋅ ; 0     0 ; 0 ; 1 ; ⋅⋅⋅ ; 0     0 ; 0 ; 0 ; ⋅⋅⋅ ; 0 A − BK = M − NK =  .  ⋅ ⋅ ⋅ ⋅   0 ; 0 ; 0 ; ⋅⋅⋅ ; 1   − (a 0 + k1 );−(a1 + k 2 );−(a 2 + k 3 ); − ⋅ ⋅ ⋅ − ; − (a n −1 + k n )

Je možno se přesvědčit, že charakteristické polynomy jsou pak rovny det[sI − ( A − BK )] = s n + (a n −1 + k n ) s n −1 + ⋅ ⋅ ⋅ + (a 2 + k 3 ) s 2 + (a1 + k 2 ) s + (a 0 + k1 ) = = (s – s1)⋅ (s – s2)⋅ (s – s3) ⋅⋅⋅(s – sn),

(3 – 7)

det[zI − ( M − NK )] = z n + (an−1 + k n ) z n−1 + ⋅ ⋅ ⋅ + (a2 + k 3 ) z 2 + (a1 + k 2 ) z + (a0 + k1 ) = (z – z1)⋅ (z – z2)⋅ (z – z3) ⋅⋅⋅(z – zn),

(3 – 8)

kde jsou s i , z j …póly charakteristického polynomu uzavřeného obvodu. Z rovnosti (3 - 7,8) je zřejmé, že volbou koeficientů stavového regulátoru je možno ovlivnit všechny koeficienty charakteristického polynomu. Označíme-li symbolem α i požadované hodnoty koeficientů charakteristické rovnice, pak koeficienty stavového regulátoru určíme z rovnosti ki + 1 = αi – ai, pro i = 0, 1, 2, ..., n – 1,

(3 – 9)

kde αi jsou zadané koeficienty charakteristického polynomu. Pokud stavový vektor není možno měřit, pak se do (3 – 1,2) za vektor x(t) resp. x(k) dosadí vektor estimovaný xˆ (t ) resp. xˆ (k ) . Praktický výpočet stavového regulátoru pro spojitou soustavu ukážeme na následujícím příkladě. Příklad 3.1

Uvažujme stavový model soustavy z Př.1.3.1.

 x&1 (t )   0 ; 1   x1 (t )  0 x& (t ) =  =  +   ⋅ u (t ); ⋅  x& 2 (t ) − 0,5;−1,5  x 2 (t ) 1 

 x1 (t )  y (t ) = [0,75;0] ⋅  .  x 2 (t )

Určete : a) Koeficienty stavového regulátoru metodou zadání požadovaných pólů charakteristické rovnice uzavřeného obvodu char. pol. = s2 + 4s + 4. 1 b) Vlastnosti obvodu pro vyrovnání nenulových počátečních podmínek x (0) =   . 1 c) Vlastnosti obvodu pro vyrovnání poruchy dU(t) = 1(t) a x(0) = 0. Řešení: a) Koeficienty stavového regulátoru vypočteme podle (3 – 9)


40

27.8.2004



i = 0; k1 = α0 – a0 = 4 – 0,5 = 3,5 i = 1; k2 = α1 – a1 = 4 – 1,5 = 2,5

b) Programové schéma obvodu v SIMULINKu je na obr.3.4. Základní informace o softwarové podpoře naleznete v kap.4. Průběh regulačních pochodů při vyrovnávání nenulové počáteču(t)

uR(t)

y( t )

x1 ( t )

x(t) x2 ( t )

Obr.3.4.1 Vyrovnání nenulové počáteční podmínky

Obr.3.4 Stavová regulace

1 0  ní podmínky x (0) =   je na obr.3.4.1. Zřejmě tedy platí lim x (t ) =   . t →∞ 1 0  c) Vyrovnání poruchy na akční veličině dU = 1(t) je na obr.3.4.2 a 3.4.3. x1 ( t ) u(t)

du = 1(t)

du = 1(t) y( t )

x2 ( t )

Obr.3.4.2 Průběh x1(t), x2(t) při x(0) = 0

Obr.3.4.3 Průběh veličin u(t) a y(t) při x(0) = 0

Z průběhů regulačních pochodů na obr.3.4.2 a 3.4.3 je vidět, že stavový regulátor navržený dle (3 – 9) nekompenzuje zcela účinek trvalé poruchy dU(t). Jak již bylo řečeno v úvodní části, je třeba zařadit do obvodu integrační (sumační) složku. Konec příkladu 3.2.2 Stavový regulátor pro vyrovnání trvalé poruchy Schéma na obr.3.2 si ještě zjednodušíme tím, že do schématu nezakreslíme estimátor viz obr.3.5. Budeme pouze předpokládat, že do stavového regulátoru vstupují buď měřené nebo odhadované stavové veličiny. Omezíme se pouze na sledování skokových změn žádané hodnoty.


41

27.8.2004



u(t)

du(t)

y(t)

Regulovaná soustava x(t )

xˆ (t )

-K u R (t )

k n +1

x n+1 (t )

∫

x& n+1 (t )

w(t)

Obr.3.5 Zjednodušené schéma stavová regulace s integrátorem Na obrázku 3.5 je vidět stavový vektor soustavy x(t ), xˆ (t ) a integrátor, jehož úkolem je zajistit vyrovnání trvalé poruchy du(t) nebo dosažení skokové žádané hodnoty w(t). Vstup do integrátoru e(t) označíme jako derivaci další (rozšiřující) stavové veličiny x& n+1 (t ) . Víme, že zařazení integrátoru zvyšuje řád soustavy o jedničku. Akční veličina vstupující do systému u(t) je rovna u(t) = uR(t) + du(t).

(3 – 10)

Stavový popis regulované soustavy vyjádříme ve tvaru

x& (t ) = Ax (t ) + Bu (t ) = Ax (t ) + B ⋅ [u R (t ) + d u (t )] .

(3 – 11)

Stavová veličina x n+1 (t ) , o kterou se rozšířil počet stavových veličin po zabudování integrátoru dle obr.3.5, splňuje rovnici x& n +1 (t ) = w(t ) − y (t ) = w(t ) − Cx(t ) , y(t) = Cx(t)

Zavedeme-li rozšířený vektor stavových veličin xn + 1(t), pak platí  x& (t )   A ; 0   x (t )  0   B B = ⋅ +   ⋅ w(t ) +   ⋅ u R (t ) +   ⋅ d u (t ) . x& n+1 (t ) =     0 0   x& n+1 (t ) − C ; 0   xn+1 (t ) 1 

(3 – 12)

Vstup uR(t) je roven  x(t )  u R (t ) = −Kx(t ) + k n+1 ⋅ xn+1 (t ) = −[ K ; − k n+1 ]⋅   = −K n+1x n+1 (t ) ,  xn+1 (t )

kde je

 x (t )  x n +1 (t ) =  ,  x n +1 (t )


(3 – 13)

K n+1 = [K ;−k n+1 ] ,

42

27.8.2004


a


xn+1(t) je rozšířený stavový vektor viz obr.3.6. Kn+1 je matice rozšířeného stavového regulátoru.

Dosazením za u(t) z rovnosti (3 – 13) do (3 – 12) dostaneme  x& (t )   A ; 0   x (t )   B   x (t )   0  B = ⋅ −   ⋅ [ K ; − k n+1 ]⋅  +   ⋅ w(t ) +   ⋅ d u (t ) x& n+1 (t ) =      0   x& n+1 (t ) − C ; 0   xn+1 (t )  0   xn+1 (t )  1 

Zavedeme-li matice An + 1, Bn + 1 můžeme pro uzavřený obvod s rozšířeným stavovým vektorem psát  0 x& n +1 (t ) = ( A n +1 − B n +1K n +1 ) ⋅ x n +1 (t ) +   ⋅ w(t ) + B n +1 ⋅ d u (t ), (3 – 14)  1  x (t )  y (t ) = [ C ; 0 ] ⋅  ,  x n +1 (t )

kde

B  B n +1 =   je vektor [(n + 1); 1]. 0 

 A ;0  An+1 =   je matice [n + 1; n + 1]; − C ; 0  Charakteristický polynom je roven

det[sI – (An+1 – Bn+1 Kn+1)] = (s – s1) (s – s2) (s – s3)... (s – sn+1). Struktura obvodu s rozšířeným vektorem stavového regulátoru je na obr.3.6. du(t)

u(t)

y(t)

Regulovaná soustava

u R (t )

x(t )

-Kn+1

xˆ (t )

x n+1 (t ) x n+1 (t )

∫

w(t)

x& n+1 (t )

Obr.3.6 Schéma stavové regulace s rozšířeným vektorem stavového regulátoru Obecně není možno dosáhnout pro výpočet koeficientů stavového regulátoru s rozšířeným vektorem stavu uzavřeného tvaru jako je (3 – 9), protože matice uzavřeného obvodu již není v normální formě rekonstruovatelnosti. Pro výpočet použijeme proto softwarové podpory MATLABu.


43

27.8.2004


Příklad 3.2


Pro regulační obvod z Př.3.1 zajistěte nulovou regulační odchylku při působení trvalých poruch. Požadované póly charakteristické rovnice uzavřeného obvodu nechť jsou s1 = s2 = s3 = -2. Stavový model má tvar  x&1 (t )   0 ; 1   x1 (t )  0 x& (t ) =  =  +   ⋅ u (t ); ⋅  x& 2 (t ) − 0,5;−1,5  x2 (t ) 1 

 x1 (t )  y (t ) = [0,75;0] ⋅  .  x2 (t )

Řešení: Podle (3 – 14) je matice A n+1 a vektor B n+1 roven A n+1

 0 ; 1 ;0  A;0   =  = − 0,5 ; − 1,5 ; 0 − C ; 0  − 0,75 ; 0 ; 0 

0    ; B = B  = 1  . n +1 0       0    

Pro výpočet koeficientů regulátoru stavu byla použita funkce acker viz obr.3.7.1. Programové schéma regulace je na obr.3.7.3. Průběhy y(t), u(t) na skok w(t) = 1(t) jsou na obr.3.7.2.

y (t ) y (t ) y(t) u(t)

Obr.3.7.1 Výpočet koeficientů matice K

Obr.3.7.2 Průběh y(t) a u(t) na skok w(t) = 1(t)

u(t)

Obr.3.7.3 Programové schéma stavové regulace s integrátorem Průběhy stavových veličin x1(t), x2(t) a y(t) při skoku žádané hodnoty jsou na obr.3.7.4. Průběhy stavových veličin x1(t), x2(t) při vstupu poruchové veličiny dU(t) = 1(t) jsou na obr.3.7.5.


44

27.8.2004



x2(t)

y(t)

x1(t)

x1(t) y(t) x2(t)

Obr.3.7.4 Průběhy x(t), y(t) při w(t) = 1(t)

Obr.3.7.5 Průběhy x(t), y(t) při dU(t) = 1(t)

Z průběhů regulačních pochodů je zřejmé, že podmínky zadání jsou splněny jak pro skok žádané hodnoty, tak i skok trvalé poruchy na akční veličině. V dalším příkladě bude proveden a diskutován návrh diskrétního stavového regulátoru. Konec příkladu 3.2.3 Diskrétní stavový regulátor pro vyrovnání trvalé poruchy Obecný postup návrhu, který je založen na hledání stavového popisu rozšířené soustavy, která obsahuje v rozšířeném vektoru soustavy i vývoj žádané hodnoty, přesahuje rámec našeho kurzu. Omezíme se pouze na změnu žádané hodnoty ve tvaru jednotkového skoku. Formálně je možno schéma na obr. 3.6 upravit tak (viz. obr.3.8), že do schématu zařadíme sumátor místo integrátoru. Výstup za sumátorem je roven součtu regulačních odchylek a platí S (k ) = e(k ) + e(k − 1) + e(k − 2) + ⋅ ⋅ ⋅ + e(1) + e(0) = S (k − 1) + e(k ), (3 – 15) S (k − 1) = e(k − 1) + e(k − 2) + ⋅ ⋅ ⋅ + e(1) + e(0). du(k )

u(k)

y(k)


x (k )

xˆ (k )

u R (k )

K S (k ) k n +1

∑ 1 z

e(k) w(k)

S (k − 1)

Obr.3.8 Zjednodušené schéma se sumátorem Dynamický vývoj výstupu sumátoru popisuje diferenční rovnice (3 – 15). Předpokládáme, že diskrétní stavový model soustavy je ve tvaru Doc. Ing. Modrlák Osvald, CSc.

45

27.8.2004



y(k) = Cx(k); xT(k) = [x1(k), x2(k), ..., xn(k)].

x(k + 1) = Mx(k) + Nu(k);

(3 – 16)

Tento model rozšíříme o dynamiku sumátoru S(k), jehož vývoj v okamžiku (k + 1) popisuje diferenční rovnice S(k + 1) = S(k) + e(k + 1) = S(k) + w(k + 1) – y(k + 1), kde e(k + 1) je regulační odchylka v čase (k + 1). Regulační odchylka e(k + 1) = w(k + 1) – y(k + 1). Předpokládáme-li, že žádaná hodnota w(k+1) je známá a y(k+1) je možno vypočítat z vektoru x(k+1) pomocí (3 – 16), pak platí S(k + 1) = S(k) + w(k + 1) – Cx(k + 1) = S(k) + w(k + 1) – CMx(k) - CNu(k). (3 – 17)

Společný stavový model pro vektor x(k) a S(k) = xn+1(k) vytvoříme zavedením rozšířeného  x( k )   x( k )  stavového vektoru x n +1 (k ) =   a použitím rovnic (3 – 16) a (3 – 17) =  S ( k )   x n +1 ( k )  0   x(k + 1)  x(k + 1)   M ; 0   x(k )   N  u (k ) +   w(k + 1) x n+1 (k + 1) =  + ⋅ = =      1   S (k + 1)   x n+1 (k + 1) − CM ; 1  x n+1 (k ) − CN  (3 – 18) Rozšířený diskrétní stavový model je možno zapsat ve tvaru

xn+1(k + 1) = Mn+1xn+1(k) + Nn+1u(k) + Nww(k + 1),

kde

(3 – 19)

 x( k )  0   N   M ; 0 ; M n +1 =  x n +1 ( k ) =  = N1; N w =   . = M1; N n +1 =     1  − CN  − CM ; 1  x n +1 ( k ) 

Akční veličina je pak rovna u(k) = –Kn+1xn+1(k) = –K1⋅xn+1(k); Struktura rozšířeného stavového modelu (rovnice (3 – 18,19)) je na obr.3.9. w(k )

Nw d u (k )

x n+1 (k )

x n +1 (k + 1)

N n+1

z −1

M n +1 u (k )

y (k )

C n +1 x n+1 (k )

− K n +1

Obr.3.9 Označení matic M1, N1, C1, K1 bude používáno při označování matic zesílení v SIMULINKu. Celý postup si ukážeme na následujícím příkladu.


46

27.8.2004



Příklad 3.3 Regulovaná soustava je popsána přenosem F =

2 . s + 3s + 2s + 1 3

2

Určete pro zadanou regulovanou soustavu diskrétní stavový regulátor, pracující s periodou vzorkování T = 1 sec., který zajistí: 1) vyrovnání nenulové počáteční podmínky x(0) a 2) vyrovnání účinků trvalé poruchy a zajištění skoku žádané hodnoty. Pro oba regulátory provedete: a) Výpočet koeficientů stavového regulátoru metodou zadání požadovaných pólů charakteristické rovnice uzavřeného obvodu pro char.pol.1 = z3 nebo char.pol.2 = z4 b) Ověření vlastností obvodu pro vyrovnání nenulových počátečních podmínek c) Ověření vlastností při vyrovnání poruchy na akční veličině dU(k) = 1(k); w(k) = 0 a průběh regulované veličiny při skoku žádané hodnoty dU(k) = 0; w(k) = 1(k) Řešení: Pro zadanou periodu vzorkování určíme diskrétní přenos a stavový popis v normální formě řiditelnosti. 1 0   0 0  0.1669 z 2 + 0,3396 z + 0,03808   G( z) = 3 ; M= 0 0 1 ; N = 0 ; 2 z − 1,305 z + 0,6271z − 0,0498 0,0498; − 0,627;1,305 1  C = [0,03808; 0,3396; 0,1669].

1) Pro vyrovnání nenulové počáteční podmínky xT = [1; 1; 1] a) Koeficienty stavového regulátoru vypočteme podle (3 – 9) i = 0; k1 = α0 – a0 = 0 – (-0,04979) = 0,04979 i = 1; k2 = α1 – a1 = 0 – 0,62712 = –0,62712 i = 2; k3 = α2 – a2 = 0 – (-1,305) = 1,305

b) Programové schéma regulačního obvodu je na obr.3.10. Blok Discret State-Space 1 představuje neregulovanou soustavu a Discret State-Space regulovanou. Pomocí matice C[1,3] se vypočítává diskrétní regulovaná veličina y(k). Koeficienty stavového regulátoru jsou reprezentovány blokem Stav. reg. Vlastnosti diskrétní soustavy jsou porovnávány s dynamikou spojité soustavy. U bloků Discret State-Space je volena matice C3x3 jako jednotková matice, čímž dosáhneme toho, že výstupem z těchto bloků jsou stavové vektory regulované a neregulované


47

27.8.2004



u(k)

x(k) y( k)

uR(k)

Obr.3.10. Programové schéma diskrétního regulačního obvodu se stavovým regulátorem soustavy. Časové průběhy regulovaného a neregulovaného x(k) pro nenulové počáteční podmínky jsou na obr.3.10.1a a jejich diskrétní hodnoty jsou na obr.3.10.1b (matice XN neregulovaný vektor, matice XR-regulovaný, první řádka v matici diskrétní čas). Průběhy regulačních pochodů pro nenulové počáteční podmínky x(0)T = [1, 1, 1] jsou na obr.3.10.1-3. c) Pro vstup trvalé poruchy dU(k) = 1(k) jsou průběhy regulačních pochodů zobrazeny na obr.3.10.4-5.

xnereg(k)

x(k) Obr.3.10.1a Regulované a nereg. x(k)

Obr.3.10.1b Diskr. hodnoty reg. a nereg. x(k) y( k)

x(k) u(k)

Obr.3.10.2 Průběh x(k) při x(0)T = [1,1,1] Obr.3.10.3 Průběhy y(k), u(k) při x(0)T = [1,1,1]


48

27.8.2004



y( t )

x(k) y( k)

u(t)

Obr.3.10.4 Průběhy y(k),y(t),u(k) při dU(k) = 1(k) a x(0) = 0

Obr.3.10.5 Průběhy x(k) při dU(k) = 1(k) a x(0) = 0

2) Vyrovnání účinků trvalé poruchy a zajištění dosažení žádané hodnoty a) Podle (3 – 19) sestavíme matice M1 = Mn+1, N1 = Nn+1 viz obr.3.11. Pro zadané póly charakteristického polynomu z1 = z2 = z3 = z4 = 0 pak pomocí funkce acker vypočteme koeficienty stavového regulátoru viz obr.3.11. Význam matice TYU bude vysvětlen.

Obr.3.11 Výpočet stavového regulátoru

y( k)

Obr. 3.12 Struktura rozšířeného diskrétního stavového modelu Př.3.3


49

27.8.2004



b) Ověření vlastností regulačního obvodu je realizováno pomocí programového schéma dle obr.3.12. Průběhy regulačních pochodů pro nenulové počáteční podmínky x(0)T = [1, 1, 1, 0] jsou na obr.3.13.1-2. Diskrétní hodnoty k, y(k), u(k) jsou na obr.3.11 v matici TYU. u(k)

x(k)

y( k)

Obr.3.13.1 Průběh x(k) při x(0)T=[1,1,1,0] Obr.3.13.2 Průběhy y(k), u(k) při x(0)T=[1,1,1,0] c) Ověření vlastností při vyrovnání skoku žádané hodnoty a trvalé poruchy. Časový průběh stavových veličin a regulované a akční veličiny při skoku žádané hodnoty v okamžiku k = 0; w(k) = 1(k); dU(k) = 0 je na obr.3.13.3 a 3.13.4 u(k) y( t ) y( t ) y( k)

Obr.3.13.3a Čas. průběh y(t), y(k),w(k)=1(k) Obr.3.13.4a Čas. průběh y(t), u(k) při w(k)=1(k) Na obr.3.13.4a je vidět, že akční veličina je o jeden krok zpožděna a tím je zpožděn celý regulační pochod. Je třeba uvést, že počáteční podmínky na všechny stavové veličiny byly položeny rovny nule, čímž došlo k uvedenému zpoždění. Pomocí vhodné volby počátečních podmínek můžeme tento jev odstranit. Zavedeme-li na čtvrtou stavovou veličinu počáteční podmínku rovnu jedné x T (0) = [0 0 0 1] , což je vzhledem k definici xn+1(k) logické, pak dosáhneme průběhu uvedeném na obr.3.13.4. Časový průběh vektoru x(k) je na

y( t )

x(k)

y( k)

Obr.3.13.3 Časový průběh x(k), w(k) = 1(k) Obr.3.13.4 Čas. průběh y(t), y(k) při w(k) = 1(k)


50

27.8.2004

uR(k)



obr.3.13.3. Průběh diskrétní akční veličiny u(k) a výstupu ze spojité soustavy je na obr.3.13.5. Na tomto místě bych rád upozornil, že námi zvolené diskrétní stavové veličiny nemají žádný fyzikální význam a není je možné tedy ani měřit. Časový průběh stavových veličin a regulované a akční veličiny při vstupu poruchy na akční veličině v okamžiku k = 0; 1(k); w(k) = 0 je na obr.3.13.6 až 3.13.8

y( t )

x(k)

Obr.3.13.6 Čas. průběh x(k) na dU(k) = 1(k)

Obr.3.13.5 Průběh u(k), y(t) na w(k) = 1(k) y( t )

y( t )

y( k)

u(k)

Obr.3.13.7 Průběh y(k), y(t) na dU(k) = 1(k)

Obr.3.13.8 Průběh y(t), u(k) na dU(k) = 1(k)

Z praktického hlediska je vhodnější simulační schéma odpovídající struktuře na obr.3.8 viz obr.3.14, protože odpovídá struktuře zpětnovazebního regulačního obvodu. u(k)

Obr.3.14 Simulační schéma se sumátorem Doc. Ing. Modrlák Osvald, CSc.

51

27.8.2004



Odezvy na skok žádané hodnoty jsou na obr.3.14.1 a 3.14.2.

y( t ) y( t ) y( k)

u(k)

Obr.3.14.1 Průběh y(k),y(t) na w(k) = 1(k);

Obr.3.14.2 Čas. průběh u(k),y(t) na w(k) = 1(k)

V tomto simulačním schématu se neprojevuje zpoždění na skok žádané hodnoty. Je třeba zdůraznit, že podrobná analýza výsledků uzavřeného obvodu s diskrétním stavovým regulátorem přesahuje rámec našeho kurzu a bude jí věnována pozornost v pátém ročníku v kurzu „Řízení ve stavovém prostoru“. Konec příkladu

3.2.4 Princip separability V kapitole 3.2.2-3 jsme se seznámili s návrhem spojitých a diskrétních stavových regulátorů metodou zadáním pólů uzavřeného obvodu. Na jednoduchých příkladech jsme demonstrovali vlastnosti regulačních obvodů s takto navrženými stavovými regulátory za předpokladu, že všechny stavové veličiny jsou měřitelné. Při praktických aplikacích stavového řízení tomu tak není a logicky se nabízí možnost použít pro odhady stavu estimátorů. Akční veličina (dU(t) = 0) je pak rovna u (t ) = −Kxˆ (t ) ,

(3 – 20)

kde xˆ (t ) … je odhadovaný vektor stavu, který je možno vyjádřit rovnicí xˆ (t ) = ( Α − LC ) xˆ (t ) + Bu (t ) + Ly (t ) .

(3 – 21)

Otázkou však zůstává, jakým způsobem použití estimátoru ovlivní dynamiku uzavřeného obvodu. Analýzu uzavřeného obvodu s estimátorem omezíme pouze pro spojitý uzavřený regulační obvod, protože závěry budou platit i pro diskrétní systém v plném rozsahu. Předpokládáme, že regulovaná soustava je popsána stavovou rovnicí x& (t ) = Ax(t ) + Bu (t ) ,

y(t) = Cx(t).

(3 – 22)

Chyba estimace spojitého estimátoru je podle (2 – 3) ∆x(t ) = x(t ) − xˆ (t ) , a její dynamika je popsána stavovou rovnicí ∆x& (t ) = ( A − LC )∆x (t ) .


52

(3 – 23) (3 – 24)

27.8.2004



Dosazením (3 – 20) do (3 – 22) a s využitím (3 – 23) dostaneme pro vektor stavu x(t) rovnici x& (t ) = (A − B ⋅ K )x(t ) + B ⋅ K∆x(t ). (3 – 25a) Chyba estimace má rovnici (3 – 24), takže stavovou rovnici uzavřeného obvodu s estimátorem tvoří soustava rovnic x& (t ) = (A − B ⋅ K )x(t ) + B ⋅ K∆x(t ) ∆x& (t ) = ( A − LC )∆x (t ) .

 x(t )  Zavedeme-li vektor stavu uzavřeného obvodu x u (t ) =   , pak je možno odvodit ∆x(t )  x& (t )  ( A − BK ); BK   x(t )  ⋅ , x& u (t ) = A K x u (t ) =  = 0 ; ( A − LC) ∆x(t ) ∆x& (t ) 

(3 – 25b)

( A − BK ); BK  kde A K =  . 0 ; ( A − LC)  Charakteristická rovnice uzavřeného obvodu je pak rovna det AK = det(sI – A + BK)⋅det(sI – A + LC) = 0.

(3 – 26)

Z rovnosti (3 – 26) je vidět, že dynamika uzavřeného regulačního obvodu s estimátorem závisí odděleně na vlastních číslech uzavřeného obvodu bez estimátoru det(sI – A + BK) = sn + αn-1sn-1+…+α1s + α0 = αR(s) a na vlastních číslech matice chyby odhadu (A – LC) det(sI – A + LC) = sn + γn-1sn-1+…+γ1s + γ0 = γE(s). Jinými slovy na základě rovnice (3 – 26) vyplývá, že návrh regulátoru i estimátoru může probíhat odděleně a to libovolně zvoleným postupem. Proto se hovoří o principu separability, tj. o odděleném návrhu stavového regulátoru a estimátoru. Struktura zpětnovazebního obvodu s estimátorem je na obr.3.15.


53

27.8.2004



u (t )

d u (t )

u R (t )

−K

Soustava

xˆ (t )

y (t )

Estimátor

Obr.3.15 Struktura zpětnovazebního regulačního obvodu s estimátorem Ukažme na následujícím příkladě metodiku návrhu zpětnovazebního obvodu se stavovým regulátorem a estimátorem.

Příklad 3.4

Uvažujte regulovanou soustavu z Př.3.3 a navrhněte 1) simulační schéma v SIMULINKu, 2) stavový regulátor pro vyrovnání trvalých poruch, 3) estimátor pracující v uzavřené regulační smyčce. 4) Ověřte vlastnosti regulačních pochodů na a) du(t) = 1(t); w(t) = 0 Řešení: b) w(t) = 1 (t); du(t) = 0 1) Simulační schéma je na obr.3.16. Vychází ze zjednodušeného schématu uzavřeného obvodu se stavovým regulátorem, které je rozšířeno o sumátor dle obr.3.8. Sumátor je realizován Z-přenosem. Výstupem sumátoru je stavová veličina S(k), kterou není nutno estimovat, což vede ke zjednodušení estimátoru v porovnání s Př.3.3. Rozšířený stavový vektor xn+1(k + 1) rozměru (4 × 1) je vytvořen ze stavového vektoru x(k) rozměru (3 × 1) a součtu S(k) za multiplexem.


54

27.8.2004



ESTIMÁTOR

Obr.3.16 Struktura zpětnovazebního obvodu se stavovým regulátorem a estimátorem 2) Stavový regulátor převezmeme z příkladu 3.3. Návrh regulátoru byl proveden zadáním pólů uzavřeného obvodu z1 = z2 = z3 = z4 = 0 dle (3 – 8). Protože stavový popis není v normální formě řiditelnosti, byla použita funkce acker viz obr.3.17a.

Obr.3.17a Výpočet stavového regulátoru 3) Návrh estimátoru provedeme zadáním pólů matice chyby estimace dle (2 – 4a). Stavový popis není v normální formě rekonstruovatelnosti, takže je nutno použít funkci acker a její vstupy vhodně modifikovat. Postup je na obr.3.17b.


55

27.8.2004



Poznámka. Výstupem funkce acker je stavový regulátor, což je vektor (1xn) a vstupem jsou matice soustavy A a vektor B. Chceme-li vypočítat matici estimátoru L, pak tato je rozměru (nx1). Vstupem funkce acker musí být transponovaná matice A a transponovaná matice C. Výstupem bude transponovaný vektor estimátoru, který je nutno před implementací transponovat.

Obr.3.17b Výpočet estimátoru pomocí funkce acker 4) Průběhy regulačních pochodů na poruchu du(t) = 1(t); w(t) = 0 jsou na obr.3.18a,b,c. a průběhy regulačních pochodů na žádanou hodnotu w(t) = 1 (t); du(t) = 0 jsou na obr.3.18d,e,f. y( t )

uR(k)

y( k)

u(k)

Obr.3.18b Odezva uR(k), u(k)

Obr.3.18a Odezva y(t), y(k)

y( t )

x(k)

y( k)

Obr.3.18c Odezva x(k)

Obr.3.18d Odezva y(t), y(k)

u(k)

x(k)

uR(k)

Obr.3.18f Odezva x(k)

Obr.3.18e Odezva uR(k), u(k) Konec příkladu


56

27.8.2004



3.3 NÁVRH STAVOVÉHO REGULÁTORU PODLE KVADRATICKÉHO KRITERIA V technické praxi nalézá široké uplatnění metoda nejmenších čtverců. Také při návrhu regulátorů jsme se setkali s kritériem jakosti regulace jako součtem kvadrátů regulačních odchylek a přírůstků akční veličiny. Při návrhu diskrétního stavového regulátoru podle kvadratického kriteria pro vyrovnání nenulového počátečního stavu pro x(0) ≠ 0, lim x(t ) = 0 , t →∞

se používá rozšířené kvadratické kriterium ve tvaru N

[

]

J (u (k )) = ∑ x T (k )Qx(k ) + u T (k )Ru(k ) + xT (N )Px( N ) ,

(3 – 27)

k =1

kde Q, R jsou pozitivně definitní symetrické matice váhových koeficientů. Matice Q, R se volí zpravidla diagonální. Pro systémy s jedním vstupem a jedním výstupem je matice R skalár. 3.3.1 Rozšířený stavový popis Chceme-li kvadratické kriterium použít i pro vyrovnání trvalých poruch a N → ∞, je možno kriterium (3 – 27) upravit do tvaru J (u k ) =

∑ [xT (k )Qx(k ) + ∆u T (k )R∆u(k )]+ xT ( N )Px( N ) , N

(3 – 28)

k =1

kde je ∆u(k) = u(k + 1) – u(k). Dynamické vlastnosti soustavy, která má jako vstupní veličinu přírůstek akčních zásahů, jsou popsány obrazovým přenosem G( z) =

Y ( z) z z bn −1 z n −1 + bn − 2 z n − 2 + ... + b1 z + b0 = ⋅ G0 ( z ) = ⋅ , ∆U ( z ) z − 1 z −1 z n + a n −1 z n −1 + ... + a1 z + a0

kde G0 ( z ) =

Y ( z ) bn −1 z n −1 + bn − 2 z n − 2 + ... + b1 z + b0 = . U ( z) z n + a n −1 z n −1 + ... + a1 z + a0

(3 – 29)

Provedeme-li roznásobení polynomů v čitateli a jmenovateli, dostaneme Z-přenos ve tvaru G( z) =

bn −1 z n + bn − 2 z n −1 + ... + b1 z 2 + b0 z Y ( z) = ∆U ( z ) z n +1 + (a n −1 − 1) z n + (a n − 2 − a n −1 ) z n −1 + ... + (a1 − a 2 ) z 2 + (a0 − a1 ) z − a0

Abychom zjednodušili výklad, bude uvažován stavový popis v normální formě rekonstruovatelnosti, který přináší řadu zjednodušení. Rozšířený vektor x(k), rozšířené matice MD a ND mají pak strukturu


57

27.8.2004



 x1 (k )   (1 − a n −1 ) ;  x (k )  ( a − a ) ; n−2  2  n −1   x3 ( k )  (a n − 2 − a n −3 ); x( k ) =  ; MD =  ⋅ ⋅ ⋅ ⋅ ⋅   x n (k )   (a − a ) ; 1 0    + a0 ;  x n +1 (k ) 

1; 0 ; ⋅ ⋅ ⋅

; 0 bn −1   b  ; 0   n−2  bn −3  ; 0   ; ND =   . (3 – 30)  ⋅   b0  ; 1    ; 0  0 

0 ; 1; ⋅ ⋅ ⋅ 0; 0; ⋅ ⋅ ⋅ 0; 0; ⋅ ⋅ ⋅ 0;

0;

⋅⋅⋅

Stavová rovnice tohoto modelu je

x(k + 1) = MD ⋅ x(k) + ND ⋅ ∆u(k),

y(k) = CD ⋅ x(k) .

Akční zásah (přírůstek akční veličiny) je roven

∆u(k) = - KD ⋅ [x(k) – wx(k)],

(3 – 31)

matice stavového regulátoru kde je KD x(k) stavový vektor wx(k) vektor žádaných hodnot jednotlivých stavových veličin. Tento vektor vyjádříme rovností wx(k) = KW ⋅ w(k),

kde je KW w(k)

(3 – 32)

vektor, který je třeba určit z ustálených hodnot žádaná hodnota regulované veličiny y(k).

3.3.2 Podmínky pro vyrovnání skoků žádaných hodnot

Omezíme se v dalším výkladu pouze na skok žádané hodnoty regulované veličiny, pro kterou platí y(k) = x1(k). Předpokládejme skokovou změnu žádané hodnoty, takže platí w(k) = 1(k) = 1, pro k ≥ 0.

Pro ustálený stav k → ∞ platí x(k + 1) = MD ⋅ x(k), ∆u(k) = 0. Jestliže má ustálená hodnota regulované veličiny býti rovna y(∞) = 1 = x1(∞) = x1 a označíme-li xi(∞) = xi, pak musí platit  an a  n −1 an − 2 x3 = x1 ⋅ (a n − 2 − a n − 3 ) + x 4  ⋅⋅⋅ ⇒ x (∞ ) =  ⋅  ⋅ ⋅⋅⋅   a1 x n = x1 ⋅ (a1 − a0 ) + x n +1  a x n +1 = x1 ⋅ a0  0

x1 = x1 ⋅ (1 − a n −1 ) + x 2 x 2 = x1 ⋅ (a n −1 − a n − 2 ) + x3


58

          

(3 – 33, 34)

27.8.2004



Je-li x1(∞) = x1 = 1, můžeme určit ustálené hodnoty zbývajících stavových veličin z rovnic (3–33). Výsledek je uložen na vektoru x(∞) viz rovnost (3 – 34). Protože ∆u(∞) = 0, musí platit x(∞) = wx(∞). Protože platí (3 – 32), musí vektor KW býti roven KW = x(∞). Z rovnosti (3 – 35) je vidět, že složky tohoto vektoru tvoří koeficienty polynomu A(z) obrazového přenosu G0(z). Struktura uvažovaného obvodu je na obr.3.19 an    a n −1  K W = x(∞ ) =  a n − 2  .   an − 3   ⋅   a0 

(3 – 35)

Vektor KW je nutno ve schématu při skokových změnách žádané hodnoty naplnit podle rovnice (3 – 35). Chci připomenout, že ve schématu není zabudován estimátor ani model spojité soustavy.

1 z

ND

∆u(k)

CD

y(t)

x(k ) MD

∆du(k) ∆u R (k )

w(k) KD

KW

Obr.3.19 Regulační obvod s přírůstkem akční veličiny 3.3.3 Algoritmy a softwarová podpora Odvození výpočetních vzorců diskrétního stavového regulátoru podle rozšířeného kvadratického kriteria přesahuje rámec našeho kurzu. Pro zájemce bude zjednodušený výklad uveden v příloze. Bude též podrobně probírán v pátém ročníku ve speciálním kurzu „Řízení ve stavovém prostoru“. Protože jsou však pro praktické aplikace významné, budou uvedeny výsledné vzorce pro návrh diskrétního stavového regulátoru podle kvadratického kriteria, což umožní v aplikacích využít softwarové podpory MATLABu. Uvažujme stavový model regulované soustavy ve tvaru

x(k + 1) = Ax(k) + Bu(k), uR(k) = –Kx(k), pak podmínce minima rozšířeného kvadratického kriteria


59

27.8.2004

Teorie automatického řízení II ∞


[

]

J (u k ) = ∑ x T (k )Qx (k ) + ∆u T (k ) R∆u (k ) + 2 x T (k ) Nu (k ) → MIN

(3 – 36)

k =1

vyhovuje stavový regulátor s maticí regulátoru K, který vypočteme ze vzorce

K = (BTXB + R)-1(BTSA + N)

(3 – 37)

a z Riccatiho rovnice

(

)

AT SA − S − AT SB + N ⋅ (BXB + R )

−1

(B

T

)

SA + N T + Q = 0

(3 – 38)

Softwarová podpora výpočtů stavových regulátorů Použití softwarových prostředků ukážeme na dvou následujících příkladech. Pro soustavu z Př.2.1 navrhněte diskrétní stavový regulátor podle kvadratického kriteria pro vyrovnání počátečního rozvážení vektoru stavu x(0). Matice Q je jednotková matice a matice R = 0,1. Získané výsledky zkontrolujte metodou zadání pólů uzavřeného obvodu. Řešení: Matice M, N, Q a R jsou

Příklad 3.5

Použitím funkce dlqr vypočítáme matici stavového regulátoru K. Póly uzavřeného obvodu jsou ve vektoru e. Dosazením do funkce acker obdržíme stavový regulátor metodou "pole placement". Je zřejmé, že výsledky v rámci přesnosti výpočtu jsou shodné. Konec příkladu


60

27.8.2004


Příklad 3.6


Uvažujte dynamickou soustavu z Př.3.3. Vypracujte:

1) Model regulačního obvodu v SIMULINKu umožňující stavovou regulaci, vstup žádané hodnoty a poruchy na akční veličině ∞

{

2) Navrhněte stavový regulátor podle kvadratického kriteria J = ∑ [ w(k ) − y (k )]2 + ∆u (k ) 2

}

k =0

3) Volte w(k) = 0 a ověřte dynamické vlastnosti obvodu při a) vstupu skokové změny žádané hodnoty w(t) = 1(t), d(t) = 0 b) vstupu poruchy w(t) = 0, d(t) = 1(t) ∞

{

}

4) Navrhněte stavový regulátor podle kvadratického kriteria J = ∑ [ w(k ) − y (k )]2 , k =0

w(k) = 0 a porovnejte dynamické vlastnosti s regulátorem dle 3). Řešení: 1) Model regulačního obvodu včetně spojité části je na obr. 3.20.1. Diskrétní stavový přírůstkový model určíme z diskrétního přenosu G(z) tím, že se vynásobí členem z/(z – 1) G( z) =

0.1669 z 2 + 0,3396 z + 0,03808 0.1669 z 3 + 0,3396 z 2 + 0,03808 z z . ⋅ = z 3 − 1,305 z 2 + 0,6271z − 0,0498 z − 1 z 4 − 2,305 z 3 + 1,9321z 2 − 0.6763 z + 0,0498

Obr.3.20 Model regulačního obvodu se stavovým regulátorem k Př.3.6 Diskrétní matice systému MD, matice buzení ND byly sestaveny dle (1 – 26) a jsou zobrazeny na obr.3.20.3. Dynamika diskrétního přírůstkového modelu je definována blokem Discrete State-Space, jehož vstupem je přírůstek akční veličiny ∆u(k). Volbou matice výstupu CD jako jednotkové matice dosáhneme toho, že na výstupu tohoto bloku je vektor stavových veličin x(k). Diskrétní výstup modelu y(k) se realizuje pomocí Matrix Gain3. Dynamika spojité soustavy je určena blokem Transfer Fcn, jehož vstup je spojitá akční veličina, která se vytvoří sumátorem Discrete Transfer Fcn. Odezva systému na puls, který byl vytvořen pomocí bloků step4, step5, je na obr.3.20.2a,b. Matice stavového regulátoru je realizována v bloku Matrix Gain4.


61

27.8.2004

∆u(k)

x(k) Teorie automatického řízení II


Obr.3.20.2a Odezva ∆u(k), y(k), y(t) na puls

Obr.3.20.2b Odezva vektoru x(k) na puls

2) Navrhněte stavový regulátor podle kvadratického kriteria pro koeficient tlumení κ = 1. Podle zadání byla sestavena matice Q a vektor NW a výpočet regulátoru byl proveden pomocí funkce dlqr viz obr.3.20.3. Stavový regulátor je uložen v matici KD.

Obr.3.20.3 Zobrazení matic MD, ND, NW, Q a příkaz výpočtu matice KD 3) Ověřte dynamické vlastnosti obvodu a) Odezva y(k), y(t), ∆u(k) na skok změny žádané hodnoty w(t) = 1(t), d(t) = 0 je na obr.3.20.4a, průběh ∆u(k) a u(k) je na obr.3.20.4b u(k)

y(t) y(k) ∆u(k)

∆u(k)

Obr.3.20.4a Časový průběh y(k), y(t), ∆u(k)

Obr.3.20.4b Časový průběh ∆u(k) a u(k)

b) Odezva y(k), y(t), ∆u(k) na vstup skokové poruchy w(t) = 0, d(t) = 1(t) je na obr.3.20.5a, průběh ∆u(k) a u(k) je na obr.3.20.5b


62

27.8.2004



y(t)

u(k)

y(k)

∆u(k)

∆u(k)

Obr.3.20.5a Odezva y(k),y(t),∆u(k) na d(t)=1(t) Obr.3.20.5b Odezva ∆u(k), u(k) na d(t) = 1(t) 4) Seřízení stavového regulátoru podle kvadratického kriteria pro κ = 0. Regulátor pro takto zadané parametry je navržen příkazem MATLABu dlqr a výsledky jsou uvedeny na obr.3.20.6.

Obr.3.20.6 Výpočet matice regulátoru pomocí příkazu dlqr s výsledky Průběhy regulačních pochodů jako odezev na skok žádané hodnoty jsou na obr.3.20.7a,b. Regulační pochody jako odezvy na skok poruchové veličiny jsou na obrázcích 3.20.8a,b. Z porovnání regulačních pochodů na obr.3.20.8a, 3.20.5a vyplývá,

y(t)

u(k)

∆u(k)

∆u(k)

Obr.3.20.7a Odezva ∆u(k),u(k) na w(k)=1(k)

y(k)

Obr.3.20.7b Odezva y(t),y(k),u(k) na w(k)=1(k)

že regulátor seřízený pro κ = 0 vyrovná poruchu v pěti akčních zásazích. Chování regulátoru odpovídá seřízení pro minimální počet kroků regulace.


63

27.8.2004



∆u(k)

x(k)

y(t) y(k)

Obr.3.20.8a Odezva ∆u(k), y(k), y(t) na puls

Obr.3.20.8b Odezva x(k) na puls Konec příkladu

3.3.4 Vyrovnávání trvalých poruch pro kvadraticky optimální řízení Vyrovnávání trvalých poruch nebo změn žádaných hodnot, tak jak bylo realizováno v kap.3.2.3, vedlo na rozšířenou soustavu a následně pak ke zvyšování řádu regulátoru. Do kriteria jakosti regulace (3 – 27) vstupuje vektor stavu rozšířený o výstup ze sumátoru a akční veličinu, čímž dochází ke změnám hodnoty kriteria a dynamiky regulačních pochodů v porovnání s regulátorem navrženým na vyrovnání rozvážení počátečního x(0). Tomuto je možno se vyhnout zvolíme-li jiný postup, než je ten, který byl uveden v kap.3.2.3. Výklad rozdělíme na nalezení: a) vhodného rozšířeného vektoru stavu uzavřeného obvodu, b) modifikované struktury a zesílení zpětnovazebního regulačního obvodu se sumátorem, pomocí kterých dosáhneme splnění požadavků na vyrovnání trvalých poruch a skoků žádané hodnoty, c) zesílení zpětné vazby, která zajistí stejnou dynamiku, jakou má zpětnovazební obvod pro vyrovnání počátečního rozvážení, definované vektorem x(0). a) Základní rozdíl spočívá ve volbě rozšířeného stavového popisu, do kterého zařadíme místo výstupu ze sumátoru S(k+1) akční veličinu u (k + 1) , pro kterou platí u(k + 1) = -Kx(k + 1) = –K[Mx(k) + Nu(k)],

(3 – 39)

kde K je vektor stavového regulátoru, který je navržen podle kriteria (3 – 27) a optimálně vyrovnává nenulový vektor počátečních rozvážení x(0). Rozšířený stavový popis uzavřeného obvodu pak je roven x(k + 1)  M ; N  x(k ) u (k + 1) = − KM;−KN  ⋅ u (k ).      

(3 – 40)

Existuje transformační matice  1 ; 0 T=  , − K ; 1


64

27.8.2004



pro kterou platí  M − NK ; N   M ; N  . T −1  T =  0 ; 0   − KM ;− KN 

(3 – 41)

Z rovnosti (3 – 41) je vidět, že vlastní čísla uzavřeného obvodu v rozšířeném tvaru jsou rovny vlastním číslům matice (M - NK) a nulovému vlastnímu číslu, odpovídajícímu počáteční podmínce u(0), která odezní v jednom kroku. b) V modifikované struktuře obvodu je uvažován účinek vektoru trvalých poruch. Dynamika soustavy s účinkem trvalých poruch v souladu s (1 – 40) má pak tvar

x(k + 1) = Mx(k) + NUu(k) + Ndd(k) = Mx(k) + Nu(k) + d,

(3 – 42)

kde N = NU, d = Ndd(k). Aby bylo možno zajistit změny žádaných hodnot, musí schéma regulačního obvodu obsahovat sumátor. Použijeme strukturu sumátoru dle obrázku 3.21, což představuje modifikovaný sumátor z kap. 3.2.3 dle obr.3.8. Dynamiku pak vyjádříme následující diferenční rovnicí S (k + 1) = S (k ) + K 2 ⋅ e(k ) = (3 – 43) = S (k ) + K 2 (w(k ) − y (k ) ) = S (k ) + K 2 w(k ) − K 2 Cx (k ), kde K2 je zatím neurčené zesílení. Protože akční zásah počítáme v okamžiku „k“, je nutno pro akční veličinu „uS(k)“, kterou generuje sumátor, zařadit za S(k+1) zpoždění o jeden krok, což je realizováno členem z-1. Takže platí uS(k) = S(k). Zavedením zpoždění vzniká při skokové změně žádané hodnoty zpoždění na akční veličině o jeden krok. Toto je možno ošetřit dopřednou vazbou přímo od vstupu žádané hodnoty w se zesílením K2 viz obr.3.21, pro kterou platí uw(k) = K2w(k). u(k)

xˆ (k )

u R (k )

+ + uw

y(k)


x (k )

K1 S (k )

1 z

S (k + 1)

e(k)

∑

w(k)

K2

u S (k ) K2

Obr.3.21 Upravené schéma se sumátorem a dopřednou vazbou se zesílením K. Doc. Ing. Modrlák Osvald, CSc.

65

27.8.2004



Stavový regulátor ve zpětné vazbě generuje signál uR(k) = –K1x(k), kde K1 je zatím neurčený vektor stavového regulátoru modifikovaného schématu. Akční veličinu vstupující do soustavy pak tvoří tři složky a je rovna u(k) = uR(k) + uS(k) + uw = –K1x(k) + S(k) +K2w(k).

(3 – 44)

c) Pro takto definovanou akční veličinu pro časový okamžik (k+1) platí u(k+1) = –K1x(k+1) + S(k+1) +K2w(k+1).

(3 – 44a)

Dosazením za S(k+1) a x(k+1) z (3 – 43) a (3 – 42) dostaneme u(k+1) = –K1Mx(k) – K1Nu(k) – K1d + S(k) + K2w(k) –K2Cx(k) +K2w(k+1). (3 – 44b) Vyloučení proměnné S(k) zajistíme dosazením za S(k) z rovnice (3 – 44) u(k+1) = (K1–K1M–K2C )x(k) + (1– K1N )u(k) – K1d+K2w(k+1).

(3 – 44c)

Rozšířená stavová rovnice uzavřeného obvodu je pak ve tvaru ; M N   x (k )  0   I   x (k + 1)  u (k + 1)  =  K − K M − K C ; 1 − K N  ⋅ u (k )  +  K  ⋅ w(k + 1) + − K  ⋅ d ,    1   2 1 2 1    1

kterou můžeme upravit do tvaru  I   I ; 0   M − I ; N   x (k )  0   x (k + 1)  I ( 1 ) ⋅ + + + ⋅ = + w k   − K  ⋅ d . (3 – 45)         u (k + 1)  1     − K 1 ;− K 2   C ; 0  u (k )   K 2  Porovnáme-li matici uzavřeného obvodu z (3 – 45) a (3 – 40), kterou upravíme do tvaru  I ; 0  M − I ; N   M ; N  = + I  − K ;− K  ⋅  C ; 0  , − KM;−KN     1 2  

(3 – 46)

je možno dosáhnout stejné dynamiky (stejných pólů charakteristické rovnice uzavřeného obvodu), jestliže platí M − I ; N [− K 1 ; − K 2 ] ⋅  (3 – 47)  = [− KM ; − I − KN ] .  C ;0 Z rovnosti (3 – 47) je možno vypočítat hledaný vektor stavového regulátoru K1, a zesílení K2 M − I ; N [ K 1 ; K 2 ] = [KM ; I + KN] ⋅    C ;0


66

−1

(3 – 48 )

27.8.2004



Rozklad stavového regulátoru K na modifikovaný stavový regulátor K1 a zesílení K2 zachovávající původní dynamické vlastnosti existuje tehdy, existuje-li inverze matice M − I ; N   C ;0  .   Podrobnější analýza diskutované problematiky přesahuje rámec základního kurzu a bude diskutována v kurzu „Řízení ve stavovém prostoru“ v devátém semestru. Aplikaci probíraného tématu ukážeme na následujících příkladech.

Příklad 3.7

Řešení:

Uvažujte regulovanou soustavu z Př.3.3 a navrhněte 1) simulační schéma modifikovaného zpětnovazebního obvodu v SIMULINKu, 2) modifikovaný stavový regulátor pro vyrovnání trvalých poruch a skoků žádaných hodnot, 3) estimátor pracující v uzavřené regulační smyčce. 4) Ověřte vlastnosti regulačních pochodů na a) w(t) = 1(t); du(t) = 0 b) du(t) = 1(t); w(t) = 0

1) Simulační schéma modifikovaného zpětnovazebního obvodu v SIMULINKu je na obr.3.22

ESTIMÁTOR

Obr.3.22 Modifikovaná struktura zpětnovazebního obvodu s estimátorem


67

27.8.2004



2) Návrh modifikovaného stavového regulátoru pro vyrovnání trvalých poruch a skoků žádaných hodnot. a) Výpočet stavového regulátoru na vyrovnání počátečního rozvážení vektoru x(0) pomocí funkce acker je na obr. 3.22.1. Póly zpětnovazebního obvodu byly voleny z1 = z2 = z3 = 0.

Obr.3.22.1 Výpočet stavového regulátoru b) Výpočet modifikovaného stavového regulátoru K1 a zesílení K2 podle (3 – 48) je na obr.3.22.2. Bylo zavedeno označení M − I ; N  M − I ; N  WU =  ; WUI =     C ;0   C ;0 

−1

; WK = [KM ; I + KN ];

Obr.3.22.2 Výpočet modifikovaného stavového regulátoru K1 a zesílení K2 3) Estimátor pracující v uzavřené regulační smyčce vypočteme pro nerozšířený stavový popis viz obr.3.22.3 pomocí funkce acker.


68

27.8.2004



Obr.3.22.3 Výpočet stavového estimátoru 4) Ověření vlastností regulačních pochodů pro a) w(t) = 1(t); du(t) = 0 je na obr 3.23.1,2,3. uS(k)

uR(k) y(t)

y(k)

uw(k) u(k)

Obr.3.23.1 Odezva y(t),y(k) na w(k)=1(k) Obr.3.23.2 Odezva uR(k), u(k) na w(k)=1(k) b) Odezvy na du(t) = 1(t); w(t) = 0 jsou na obr.3.23.4,5,6.

x(k)

x(k)

Obr.3.23.3 Odezva x(k) na w(k) = 1(k)

Obr.3.23.4 Odezva x(k) na du(k) = 1(k) uR(k)

u w (k ) uw(k)

u(k) Obr.3.23.5 Odezva y(t),y(k) na du(k)=1(k)

uS(k)

Obr.3.23.6 Odezva uR(k), u(k) na du(k)=1(k)

Konec příkladu


69

27.8.2004


Příklad 3.8


Uvažujte regulovanou soustavu z Př.3.3 s tím, že na soustavu působí ještě neměřená porucha „d“, jejíž dynamické účinky jsou vyjádřeny přenosem

Fd ( s ) =

2s + 1 s + 3s 2 + 2s + 1 3

Úkol: 1) Navrhněte stavový regulátor podle kvadratického kriteria ∞

[

]

J = ∑ y (k )T Qy (k ) + u (k )T Ru (k ) pro Q = 1, R = 1. 1

2) Nalezněte modifikovaný stavový regulátor K1 a zesílení K2 pro vyrovnání trvalých poruch a skoků žádaných hodnot. 3) Vytvořte a použijte estimátor pomocí funkce estim (jako SS objekt) pracující v uzavřené regulační smyčce a nakreslete simulační schéma. 4) Ověřte vlastnosti obvodu pro vstupní veličiny a) w(t) = 1(t); du(t) = 0; d(t) = 0 b) du(t) = 1(t); w(t) = 0; d(t) = 0 c) d(t) = 1(t); du(t) = 0; w(t) = 0.

Řešení: 1) Stavový regulátor podle kvadratického kriteria. Protože porucha „d“ je neměřená, nebude zahrnuta do stavového popisu a nebude také zaváděna do estimátoru. Matice M, N, C a vytvoření diskrétního stavového popisu systému je na obr. 3.24.1. Výpočet stavového regulátoru pomocí funkce lqry pro zadané parametry je na obr.3.24.2.

Obr.3.24.1 Vytvoření diskrétního stavového popisu systému ssd

Obr.3.24.2 Výpočet stavového regulátoru pomocí funkce lqry 2) Modifikovaný stavový regulátor pro vyrovnání trvalých poruch a skoků žádaných hodnot. Výpočet modifikovaného regulátoru je na obr. 3.24.3. Doc. Ing. Modrlák Osvald, CSc.

70

27.8.2004



Obr.3.24.3 Výpočet modifikovaného stavového regulátoru 3) Estimátor jako SS objekt pomocí funkce estim vytvoříme následujícím způsobem (viz. kap.4.2). Bude použita matice estimátoru L z Př.3.7. Matice estimátoru a aplikace příkazu estim je na obr. 3.24.4. Vstupem je jedna akční veličina u, které odpovídá proměnná known = 1, a jedna měřená veličina y(k), které odpovídá proměnná sensor = 1. Obr.3.24.4 Estimátor jako SS objekt Vytvoření matic a, b, c, d pro definování SS objektu je realizováno příkazem ssdata na obr.3.24.5. Simulační schéma v SIMULINKu je na obr.3.24.6

Obr.3.24.5 Vytvoření matic a, b, c, d


71

27.8.2004



ESTIMÁTOR

Obr.3.24.6 Schéma zpětnovazebního obvodu s estimátorem jako SS objektem 4) Regulační pochody pro w(t) = 1(t); du(t) = 0; d(t) = 0 jsou na obr.3.25.1-3 a pro du(t) = 1(t); w(t) = 0; d(t) = 0 na obr.3.25.4-6, pro d(t) = 1(t); w(t) = 0; du(t) = 0 jsou na obr.3.25.7,8. uR(t)

y(t)

uS(t) uw(t)

y(k) y(k)

u (t)

Obr.3.25.1 Odezva y(t), y(k) na w(k) = 1(k)

Obr.3.25.2 Odezva uR(k), u(k) na w(k) = 1(k)

x(k) x(k)

Obr.3.25.3 Odezva x(k) na w(k) = 1(k)


Obr.3.25.4 Odezva x(k) na du(k) = 1(k)

72

27.8.2004



uR(k) y(t)

y(k)

uw(k) u(k) uS(k)

Obr.3.25.5 Odezva y(t), y(k) na du(k) = 1(k)

Obr.3.25.6 Odezva uR(k), u(k) na du(k) = 1(k)

u (k )

u(k)

y(t) y(k)

Obr.3.25.7 Odezva y(t), y(k) na d(k) = 1(k)

Obr.3.25.8 Odezva u(k) na d(k) = 1(k)

Z odezev je vidět, že jak skok žádané hodnoty tak i trvalé poruchy d(k), du(k) jsou vyrovnány. Konec příkladu V závěru této kapitoly je možno konstatovat, že takto navržený stavový regulátor zajistí vyrovnání trvalých poruch a skoků žádané hodnoty. Regulace je robustní v tom smyslu, že vyrovná i neměřenou trvalou poruchu, která není zavedena do estimátoru. Jak již bylo řečeno, podrobnější analýza přesahuje rámec našeho předmětu. Protože však klademe důraz na praktické aplikace, shrneme návrh stavového regulátoru do následujícího postupu: 1) Provede se návrh kvadraticky optimálního regulátoru řešením Riccatiho rovnice (3 – 38) a rovnosti (3 – 37), jehož výsledkem je regulátor K, který zajistí pouze optimální přechod z daného počátečního stavu x(0) do počátku souřadnic stavového prostoru. 2) Je-li požadováno vyrovnání trvalých poruch a skoků žádané hodnoty, je nutno provést modifikaci struktury zpětné vazby dle obr.3.21. 3) Podle (3 – 48 ) určit zesílení K1 a K2 modifikované zpětné vazby, která zajistí stejnou dynamiku modifikované zpětné vazby, jakou má zpětnovazební obvod pro vyrovnání počátečního rozvážení, definované vektorem x(0)


73

27.8.2004



4. SOFTWAROVÁ PODPORA ANALÝZY A SYNTÉZY VE STAVOVÉM PROSTORU Prostředí MATLABu ve spojením s „Control System Toolboxem“ umožňuje v rámci cvičení základního kurzu aplikovat softwarovou podporu jak pro spojitý tak i diskrétní stavový popis. Je možné též realizovat pomocí „Real Time Toolboxu“ přímo řízení a regulaci na úlohách v laboratořích řízení. V této kapitole je uveden výběr jen základních příkazů a stručný popis s příklady. Další funkce a podrobnosti jsou k dispozici v "Helpu" .

4.1 STAVOVÝ POPIS A TRANSFORMACE Předpokládá se základní znalost funkcí tf, c2d, arx, present, th2arx, th2tf, th2poly, které doplníme o následující funkce: Funkce ss

a) Vytvoření spojitého nebo diskrétního stavového modelu ze zadaných matic stavového popisu

x (k + 1) = Mx (k ) + Nu(k ) y (k ) = Cx (k ) + Du(k )

x& (t ) = Ax (t ) + Bu(t ) y (t ) = Cx (t ) + Du(t ) Syntaxe příkazu sys = ss(A,B,C,D) sys = ss(M,N,C,D,Ts)

kde je A B C D

matice soustavy matice buzení matice výstupu matice převodu

M N Ts

matice diskrétní soustavy diskrétní matice buzení vzorkovací perioda

Př .4.1

Konec příkladu


74

27.8.2004



b) Vytvoření spojitého nebo diskrétního stavového modelu ze spojitého nebo diskrétního přenosu

Syntaxe příkazu

sysd = ss(sd)

kde sd je spojitý nebo diskrétní přenos. Použití této funkce je demonstrováno v Př.4.2. Př .4.2

Konec příkladu Funkce tfdata

Vypíše polynomy jmenovatele a čitatele systému, který je ve tvaru tf, ss, zpk

Syntaxe příkazu kde je A B sys 'v' Př .4.3

[B,A] = tfdata (sys,'v') polynom jmenovatele v kladných mocninách "z" polynom čitatele v kladných mocninách "z" model ve tvaru tf, ss, zpk požadavek na přímé vyjádření koeficientů polynomů

Ke stavovému modelu z Př.4.1 určete polynomy diskrétního přenosu B, A.

Konec příkladu


75

27.8.2004


Funkce ssdata


Vytvoří a vypíše matice A, B, C, D systému, který je ve tvaru ss, zpk, tf

Syntaxe příkazu [A,B,C,D] = ssdata (sys)

kde je A,B,C,D sys Př .4.4

matice spojitého nebo diskrétního systému model ve tvaru tf, ss, zpk

Ze stavového modelu ss1 vypište matice M, N, C, D pro další aplikace v MATLABu. Použijte funkci ssdata.

Konec příkladu

Funkce estim

Vytvoří stavový estimátor na základě modelu ve stavovém popisu sys a matice estimátoru L.

Syntaxe příkazu EST= estim (sys,L,sensors,known)

kde je EST sys L sensors known

jméno estimátoru model ve tvaru ss matice estimátoru počet měřených výstupů y(k) soustavy počet známých (měřených) vstupů u(k) do soustavy

Výstupem z tohoto estimátoru je vektor, který obsahuje odhadovaný vektor výstupu yˆ (k ) a odhadovaný vektor stavu xˆ (k ) . Při implementaci v SIMULINKu se chová jako SS-objekt, kde je nutno definovat matice a, b, c, d, viz Př.3.8. Příkaz estim je třeba doplnit příkazem: [a,b,c,d,T] = ssdata(EST)


76

27.8.2004



4.2 NÁVRH ESTIMÁTORŮ Softwarová podpora je zaměřena především na syntézu Kalmanova filtru pomocí funkcí kalman, kalmd pro spojitou i diskrétní verzi a na vytvoření estimátoru (jeho struktury pro vypočtenou matici estimátoru a zadaný model soustavy) funkcí estim. V našem popisu se omezíme pouze na funkci kalman.

Funkce kalman

Výpočet matice Kalmanova estimátoru pro soustavy s modelem

x& = Ax + Bu + Gω , x(k + 1) = Mx(k ) + Nu(k ) + Gω (k ) ,

Stavová rovnice:

(4 –1)

Měřený výstup: y v = Cx + Du + Hω + ν , y v (k ) = Cx(k ) + Du(k ) + Hω (k ) + ν (k ) , (4 – 2) kde u je měřený-známý vstup do soustavy. Charakteristiky šumových signálů jsou

{ }

{ }

{ }

Ε{ω } = 0, Ε{ν} = 0, Ε ωω T = Q, Ε νν T = R, Ε ωv T = N .

(4 – 3)

Funkce kalman minimalizuje střední kvadratickou chybu kovariační matice

{

}

P = lim Ε [x − xˆ ] ⋅ [x − xˆ ]T . t →∞

(4 – 4)

Pro diskrétní systém platí

{

}

P = lim Ε e[k / k − 1] ⋅ e[k / k − 1]T , t →∞

{

}

Z = lim Ε e[k / k ] ⋅ e[k / k ]T , t →∞

e[k / k − 1] = x[k ] − x[k / k − 1]

(4 – 5)

e[k / k ] = x[k ] − x[k / k ] .

(4 – 6)

Syntaxe příkazu [kest,L,P] = kalman (sys,Qn,Rn,Nn) [kest,L,P,Z] = kalman (sys,Qn,Rn,Nn)

kde je kest L P sys Qn Rn Nn

model Kalmanova estimátoru matice Kalmanova estimátoru kovariační matice chyby estimace stavový model soustavy ve spojitém nebo diskrétním tvaru rozptyl bílého šumu ω kovariační matice vektoru šumu ν, ( p , 1) počet řádků je dán počtem výstupů ze soustavy MIMO kovariační matice Ε ων T

Z

Z = lim Ε e[k / k ] ⋅ e[k / k ]T ,

Př .4.5.

t →∞

{

{ }

}

e[k / k ] = x[k ] − x[k / k ]

Vypočtěte ke stavovému modelu ssd matici L Kalmanova estimátoru. Soustava má jeden vstup a jeden výstup: Zvolíme: Qn = 1, Rn = 1, Nn nezadáme.


77

27.8.2004



Konec příkladu

4.3 NÁVRH STAVOVÝCH REGULÁTORŮ Syntéza stavových regulátorů metodou zadávání pólů je reprezentována funkcemi acker, place a metody výpočtu podle kvadratického kriteria funkcemi lqr, dlqr, lqry, lqrd. V této kapitole budou uvedeny jen funkce acker, place a dlqr.

Funkce acker

Výpočet matice spojitého nebo diskrétního stavového regulátoru soustavy s jedním vstupem a výstupem metodou zadání pólů charakteristické rovnice uzavřeného obvodu.

Syntaxe příkazu K = acker(A,B,P)


78

27.8.2004


kde je K A,B P


matice stavového regulátoru matice soustavy a řízení spojitého nebo diskrétního stavového popisu vektor obsahující požadované póly charakteristické rovnice uzavřeného obvodu.

Př .4.6

Konec příkladu

Funkce place

Výpočet matice spojitého nebo diskrétního stavového regulátoru soustavy s více vstupy a výstupy metodou zadání pólů charakteristické rovnice uzavřeného obvodu.

Syntaxe příkazu K = place (A,B,P) [K,prec,message] = place (A,B,P)

kde je K A,B P prec message

matice stavového regulátoru matice soustavy a řízení spojitého nebo diskrétního stavového popisu vektor obsahující požadované póly char. rovnice uzavřeného obvodu udává dosaženou přesnost vypočtených kořenů na desetinná místa obsahuje hlášení, že vypočtený kořen se liší o 10% od požadovaného

Př .4.7


79

27.8.2004



Konec příkladu

Funkce dlqr

Výpočet matice diskrétního stavového regulátoru u (k ) = − Kx (k ) podle kvadratického kriteria (3 –27) ∞

[

]

J = ∑ x (k ) T Qx (k ) + u(k ) T Ru(k ) + 2 x (k ) T Nu(k ) ,

(4 – 7)

1

řešení Riccatiho rovnice a výpočet vlastních čísel uzavřeného obvodu. Syntaxe příkazu [K,S,e] = dlqr (a,b,Q,R) [K,S,e] = lqrd (a,b,Q,R,N) kde je K S e a,b Q R N Př .4.8

optimální diskrétní stavový regulátor minimalizující kriterium (4 – 7) řešení Riccatiho rovnice vlastní čísla uzavřeného obvodu matice systému, matice buzení váhová matice stavových veličin v kvadratickém kriteriu (4 – 7) váhová matice (koeficient) u akční veličiny v kvadratickém kriteriu (4 – 7) váhová matice (koeficient) bilineárního členu v kvadratickém kriteriu (4 – 7)

Pro soustavu z Př. 4.9 vypočtěte stavový regulátor pro zadané matice Q, R

Konec příkladu


80

27.8.2004


Př .4.9


Pro spojitou soustavu vypočtěte stavový regulátor .

Funkce lqry

Konec příkladu

Výpočet matice spojitého nebo diskrétního stavového regulátoru u = −Kx , který minimalizuje kvadratické kriterium. ∞

Spojitý systém:

(

)

J (u ) = ∫ y T Qy + u T Ru + 2 y T Nu dt

(4 – 8)

0

∞

Diskrétní systém:

[

]

J = ∑ y (k ) T Qy (k ) + u(k ) T Ru(k ) + 2 y (k ) T Nu(k ) ,

(4 – 9)

1

Syntaxe příkazu [K,S,e] = lqry (sys,Q,R) [K,S,e] = lqry (sys,Q,R,N)

kde je K S e sys Q R N

Př .4.10

optimální spojitý nebo diskrétní stavový regulátor minimalizující kriterium (4 – 9) řešení Riccatiho rovnice vlastní čísla uzavřeného obvodu: e = eig ( A − B ⋅ K ) definuje dynamický systém ( pomocí matic A, B, C, D nebo M, N, C, D) funkce ss váhová matice stavových veličin v kvadratickém kriteriu (4 – 8) váhová matice (koeficient) u akční veličiny v kvadratickém kriteriu (4 – 8) váhová matice (koeficient) bilineárního členu v kvadratickém kriteriu (4 – 8)

Pro zadaný diskrétní dynamický systém vypočtěte stavový regulátor.


81

27.8.2004



Konec příkladu

V závěru části „Základy analýzy a syntézy ve stavovém prostoru“ mohu pouze konstatovat, že kurz umožní pouze základní orientaci v problematice stavového řízení. Teoretický výklad byl omezen na minimum a v rámci softwarové podpory probíhají cvičení v laboratořích řídící techniky. Prohloubení jak teoretických tak i praktických znalostí mohou zájemci získat studiem literatury a v kurzu „Řízení ve stavovém prostoru“ v devátém semestru.

LITERATURA [1] STREJC V.: Stavová teorie lineárního diskrétního řízení. Academie Praha,1978. [2] HANUŠ,B., BALÁTĚ,J., ŠVARC,I., ZIKEŠ,F.: Teorie automatického řízení I. I.část. Skripta Liberec, 1982 [3] HAVLENA,V., ŠTECHA, J.: Moderní teorie řízení. Vydavatelství ČVUT, Praha 1996. [4] LUENBERGER D.G.: Observers for multivariable systems. IEEE Trans.of Automatic Control,(1967),(190-197) [5] PHILLPS, CH., L., HARBOR,R.,D.: Feedback Control Systems. Prentice Hall, 2000. [6] BURL, J. B.: Linear Optimal Control. Addison Wesley Longman, 1999. [7] GRACE,A., LAUB,J.A., LITTLE,J.N., THOMPSON,C.M.: Control System Toolbox. For Use with MATLAB. User's Guide


82

27.8.2004

ZÁKLADY ANALÝZY A SYNTÉZY VE STAVOVÉM PROSTORU

Recommend Documents