PANNON EGYETEM Mérnöki Kar JÁRMŰRENDSZERTECHNIKAI LABORATÓRIUM
DIGITÁLIS SZABÁLYOZÁSI ALGORITMUSOK GYAKORLATI TERVEZÉSE
dr. Vass József
Veszprém, 2012.
Köszönetnyilvánítás TÁMOP-4.2.1/B-09/1/KONV-2010-0003 Mobilitás és környezet: Járműipari, energetikai és környezeti kutatások a Közép- és Nyugat-Dunántúli Régióban. A projekt a Magyar Állam és az Európai Unió támogatásával, az Európai Szociális Alap társfinanszírozásával valósul meg.
Tartalomjegyzék
_Toc321649288 1. A MINTAVÉTELES RENDSZEREK ALAPVETŐ TULAJDONSÁGAI ..........................6 1.1. A mintavételezés ..........................................................................................................6 1.2. A mintavételezett szabályozási rendszer .......................................................................6 1.3.Tartószervek ............................................................................................................... 10 1.4. A Z- transzformáció ................................................................................................... 11 1.4.1. A Z- transzformáció néhány alapvető tulajdonsága: .............................................12 1.4.2. A legfontosabb időfüggvények Laplace- és Z- transzformáltjai (1.1. táblázat): ....13 1.4.3. Az inverz Z- transzformáció ................................................................................ 15 1.5. Az impulzus-átviteli függvény ................................................................................... 16 1.5.1. Zárt digitális szabályozási körök impulzus-átviteli függvényei: ........................... 19 1.6. A mintavételezett rendszerek leírása differencia egyenletekkel: .................................22 2. DINAMIKUS RENDSZEREK INDENTIFIKÁCIÓJA..................................................... 24 2.1. Tipikus diszkrét idejű rendszermodellek ....................................................................25 2.1.1. OE modell ...........................................................................................................26 2.1.2. ARX modell ........................................................................................................27 2.1.3. ARMAX modell..................................................................................................27 2.1.4. Box-Jenkins (BJ) modell ..................................................................................... 29 2.2. Az identifikáció gyakorlati alkalmazásai ....................................................................31 2.2.1. Általános feladatok.............................................................................................. 31 2.2.2. MATLAB System IdentificationToolbox és SIMULINK alkalmazása ................. 35 2.2.3. Nemlineáris rendszerek identifikációja ................................................................ 49 Irodalomjegyzék ............................................................................................................... 51 3. DIGITÁLIS PID SZABÁLYOZÁSI ALGORITMUSOK ................................................. 52 3.1. Folytonos PID algoritmusokból származtatott digitális PID algoritmusok .................. 52 3.2. Általános digitális szabályozási algoritmusok impulzus-átviteli függvényei: ..............54 3.3. Első és másodrendű szabályozási algoritmusok .......................................................... 55 3.4. A diszkrét PID algoritmusok modifikációi .................................................................60 3.4.1. Egyéb P- I - D struktúrák..................................................................................... 61 3.5. Az I-tag telítődésének problémái (antiwindup) ........................................................... 63 3.6. Digitális PID algoritmusok hangolása - beállítási szabályai ........................................66
3.7. PID algoritmusok paraméterillesztése kritériumok alapján .........................................72 4. DAHLIN ALGORITMUS ................................................................................................ 77 4.1. A klasszikus Dahlin algoritmus elve ..........................................................................77 4.2. Példák ........................................................................................................................ 79 4.3. MATLAB program analitikus tervezéshez .................................................................86 4.4. Esettanulmány ...........................................................................................................89 4.4.1. A kísérleti rendszer felépítése ..............................................................................90 4.4.2. Az identifikáció ................................................................................................... 92 4.4.3. Diszkrét algoritmus analitikus tervezése .............................................................. 95 4.4.4. A zártköri viselkedés vizsgálata...........................................................................99 Irodalomjegyzék ............................................................................................................. 101 5. VÉGES BEÁLLÁSI IDEJŰ (DEADBEAT) SZABÁLYOZÁSI ALGORITMUSOK ..... 102 5.1. MATLAB program és mintapélda DEADBEAT algoritmus tervezéséhez ................ 109 Irodalomjegyzék ............................................................................................................. 113 6. BELSŐ MODELLT ALKALMAZÓ (IMC) DIGITÁLIS SZABÁLYOZÁSI ALGORITMUSOK ............................................................................................................ 114 6.1. Az IMC algoritmusok alapelve................................................................................. 114 6.2. MATLAB program digitális IMC algoritmus tervezéshez ........................................ 116 6.3.Esettanulmány .......................................................................................................... 120 Irodalomjegyzék ............................................................................................................. 123 7. NAGY HOLTIDŐVEL RENDELKEZŐ RENDSZEREK DIGITÁLIS SZABÁLYOZÓI ........................................................................................................................................... 124 7.1. A holtidős rendszerek tulajdonságai ......................................................................... 124 7.2. A SMITH prediktor ................................................................................................. 125 Irodalomjegyzék ............................................................................................................. 130 8. MINIMÁLIS SZÓRÁSÚ SZABÁLYOZÁSI ALGORITMUSOK .................................. 131 8.1. Minimális szórású szabályozók holtidő nélküli rendszerekhez.................................. 132 8.1.1. Stabilitás ........................................................................................................... 135 8.2. Minimális szórású algoritmusok holtidős rendszerek esetén ..................................... 140 Irodalomjegyzék ............................................................................................................. 144 9. ELŐRECSATOLÁSOK TERVEZÉSE ........................................................................... 146 9.1. Dinamikus előrecsatolás algoritmusa ....................................................................... 147 9.2. Előrecsatolási algoritmusok és visszacsatolások együttes használata ........................ 148 10. TÖBBVÁLTOZÓS RENDSZEREK DIGITÁLIS SZABÁLYOZÁSA ......................... 153
10.1. A többváltozós rendszerek struktúrái...................................................................... 153 10.2. MIMO rendszerek belső kereszthatásainak problémái visszacsatolt rendszerek esetén ....................................................................................................................................... 156 10.3. Az RGA (RelativeGainArray) ................................................................................ 159 10.4. Statikus kompenzátor tervezése MIMO rendszerekhez ........................................... 160 10.5. Dinamikus kompenzátor tervezése MIMO rendszerekhez ...................................... 164 10.6. Digitális MIMO PI szabályozási algoritmusok ....................................................... 168 Irodalomjegyzék ............................................................................................................. 173 11. DISZKRÉT ÁLLAPOTSZABÁLYOZÓK TERVEZÉSE ............................................. 174 11.1. Diszkrét állapotmodellek MATLAB/SIMULINK környezetben ............................. 174 11.2. Diszkrét állapotszabályozók ................................................................................... 178 11.3. Diszkrét állapot megfigyelők ................................................................................. 181 11.4. Diszkrét állapotszabályozó LUENBERGER megfigyelővel ................................... 184 11.5. Esettanulmány ....................................................................................................... 186 Irodalomjegyzék ............................................................................................................. 194
A mintavételezett rendszerek alapvető tulajdonságai
1. A MINTAVÉTELES RENDSZEREK ALAPVETŐ TULAJDONSÁGAI A számítógépes folyamatirányítás során meg kell határoznunk a rendszer matematikai modelljét, valamint meg kell ismernünk a szabályozási folyamat alatt használt jelek matematikai leírását. Mivel a szabályozott rendszer általában folytonos működésű, a rendszerbe történő beavatkozás pedig digitális számítógépekkel történik, első lépésben azokat a folyamatokat kell megismernünk, amelyek összehangolják, illesztik a folytonos működésű, szabályozni kívánt objektum és a szabályozást végző digitális számítógép működését.
1.1. A mintavételezés A fent említett probléma kiküszöbölése úgy történik, hogy a számítógép által feldolgozandó folytonos jelekből bizonyos időközönként mintát veszünk. Ezáltal a folytonos jelből egy időben diszkrét jelsorozatot kapunk. Ezt a folyamatot nevezzük mintavételezésnek. Az ún. egyszerű vagy közönséges mintavételezéskor a folytonos jelből állandó T időközönként mintát veszünk. Ezt a T időt mintavételezési időnek nevezzük és T0-al jelöljük. Ez a mintavételezési eljárás a leggyakoribb és a legtöbbet használt és a továbbiakban mi is ezt a fajta mintavételezést fogjuk használni.
1.2. A mintavételezett szabályozási rendszer Az 1.1. ábrán a mintavételezett szabályozási rendszer egyszerűsített hatásvázlata látható.
Gc(z) r(k)
e(k)
In1
Go(s)
Out1
y(t)
u(k)
Alapjel
In1 DA konverter
SZABÁLYOZÓ ALGORITMUS
Nulladrendű tartószerv
Out1
Kimenet SZABÁLYOZOTT SZAKASZ
y(k) AD konverter
Mintavevő
1.1. ábra. A mintavételezett (diszkrét) szabályozási rendszer struktúrája
6
A mintavételezett rendszerek alapvető tulajdonságai
ahol a:
szabályozott szakasz (folytonos működésű, tartalmazza a végrehajtó szervet illetve az érzékelőt és távadót, lineáris esetben folytonos átviteli függvénnyel jellemezhető)
szabályozó algoritmus (diszkrét működésű, impulzus/diszkrét átviteli függvénnyel jellemezhető lineáris tag)
nulladrendű tartószerv (a szakaszos jelet folyamatossá alakítja)
mintavevőszerv (T0 mintavételezési idővel mintát vesz a távadó kimenetéből)
A/D konverter (az feszültségimpulzusokat digitális jellé alakítja)
D/A konverter (digitális jelet feszültségimpulzusokká alakítja),
a.
b.
7
A mintavételezett rendszerek alapvető tulajdonságai
c.
... d. 1.2. ábra. A digitális szabályozási kör jelei (a-y(t),b-y(k),c-u(k),d-u(t))
a rendszer jelei (1.2. ábra):
y(t) – a szabályozott szakasz folytonos kimenete (pl. a távadó mA-es jeltartományú kimenete)
y(k) – a szabályozott szakasz mintavételezett kimenőjele
u(k) – a digitális szabályozási algoritmus kimenete (digitális végrehajtójel)
r(k) – digitális referenciajel/alapjel/parancsolt érték
e(k) – digitális hibajel( e(k)= r(k)-y(k) )
A fenti rendszerben a folyamatirányító számítógép az e(k) hibajel képzést illetve az u(k) digitális végrehajtójel számítását végzi el. Az A/D, D/A konverterek, a mintavevő szerv és a nulladrendű tartószerv elektronikája a számítógép be/kimeneti (I/O) perifériájában helyezkedik el. Ezen utóbbiak a számítógép számára a szabályozott szakaszhoz tartoznak és szabályozott objektumként jelöljük az
8
A mintavételezett rendszerek alapvető tulajdonságai
algoritmus tervezési módszerek bemutatásakor. Így a folyamatirányító számítógép és a folyamatperifériákkal kibővített zártköri modell a 1.3. ábra szerinti struktúrával bír.
Go(s)
Gc(z) r(k)
e(k)
In1
u(k)
Out1
y(t) In1
Out1
Alapjel
Kimenet SZABÁLYOZÓ ALGORITMUS
OBJEKTUM y(k)
1.3. ábra. A digitális szabályozó kör A mintavételezési eljárás, mint azt a 1.4. ábra is mutatja, egység-impulzussorozat amplitúdó modulációjaként is felfogható.
1.4. ábra. A mintavételezés elve
A bemenő jel f(t) folytonos időfüggvény, a moduláló jel i*(t)= (t nT0 ) egységimpulzusn 0
sorozat, a kimenőjel az f*(t)=f(t)i* (t)=f(t) (t nT0 ) a modulált impulzussorozat. n 0
9
A mintavételezett rendszerek alapvető tulajdonságai
1.3.Tartószervek A mintavételezésnél láttuk, ahhoz, hogy a számítógép a bemenetére érkező jeleket fel tudja dolgozni, azokat a jeleket mintavételezni és a számítógép által felhasználható jelsorozattá kell alakítani. A szabályozást végző számítógép kimenete is digitális kódolású diszkrét idejű jel, impulzus sorozat. Ahhoz, hogy ez a jel a szabályozott objektum számára feldolgozható legyen, folytonos formában kell megjelennie, tehát valamilyen jelátalakítás szükséges. Ezt az átalakítást egy D/A (digitális/analóg) konverter végzi el. Ez a művelet két jól elkülöníthető részre bontható, a dekódolásra és a tartásra. A dekódolás a digitális diszkrét idejű jelből analóg kódolású impulzussorozatot állít elő, például feszültségimpulzusokat. A tartás során a tartószerv ebből a feszültség impulzussorozatból folytonos idejű jelet állít elő. Ez a művelet nem egyértelmű, mivel két mintavételezési pont között a folytonos jel alakja tetszőleges lehet. A jelalak képzésének a függvényében háromféle tartószervet különböztetünk meg. A legegyszerűbb a nulladrendű tartószerv, amely az kT0és (k+1)T0 a jelnek az kT0időpontban felvett értékét tartja, így a jelet lépcsős görbeként állítja elő (1.5. ábra).
1.5. ábra. A nulladrendű tartószerv által szolgáltatott jel A nulladrendű tartószerv az impulzusokat T0 ideig integrálja. A kT0 és (k+1)T0 határok közötti jelátvitel, két egymástól T0 idővel eltolt kezdőpontú integrál különbségeként írható le, ezért a zérusrendű tartószerv folytonos átviteli függvénye (1.1): G H (s)
1 esT0 s
10
(1.1)
A mintavételezett rendszerek alapvető tulajdonságai
Elsőrendű tartószerv esetén a tartószerv a kT0 és az (k+1)T0 időpontok közé az x(kT0) és x((k-1)T0) értékekből egyenest, míg másodrendű x(kT0),x((k-1)T0) és x((k-2)T0) értékekből parabolát extrapolál. A magasabbrendű tartásnak a célja az x(t) görbealak minél tökéletesebb visszaállítása. Az extrapoláció viszont a múltbeli viselkedés előrevetítése, amely akkor hatásos, ha a görbe a továbbiakban is az extrapolációhoz használt pontokban megnyilvánuló tendenciát követi. Ezekből a pontokból nem jelezhető előre a görbealak váratlan megváltozása, amely a rekonstrukció által visszaadott görbe alakhűségét ronthatja.
1.4. A Z- transzformáció Az egyszerű vagy közönséges mintavételezési forma tulajdonképpen nem más, mint egy kapcsoló, amely T időközönként t ideig zárva van. A T időt mintavételezési időnek nevezzük és T0-lal jelöljük (1.6. ábra).
1.6. ábra. Mintavételezés Az időtartományban a mintavételezett jelet Dirac-impulzusokból álló sorozat írja le, a továbbiakban jelöljük * indexszel. Így egy diszkrét jelsorozat a következő formába írható (1.2): x* (t) = x(0)(t) + x(T0)(t-T0) + x(2T0)(t-2T0) …
(1.2)
A frekvencia tartományban a Laplace-transzformáció eltolási tétele szerint, T0 idővel való késleltetés e-sTo–lal történő szorzással egyenértékű. Ezáltal a (1.2) Laplacetranszformáltja (1.3): x*(s) = x(0) + x(T0) e-sTo+ x(2T0) e-2sTo+ …
11
(1.3)
A mintavételezett rendszerek alapvető tulajdonságai
Ezután, már csak a z = esTo vagy a z-1 = e-sTo helyettesítéssel élve, az alábbi függvényhez jutunk (1.4): x*(z) = x(0) + x(T0) z-1 + x(2T0) z-2+ …
(1.4)
Ezt a függvényt az x*(t) függvény Z- transzformációjának nevezzük. A képletet zárt formába írva (1.5): x*(z) =
x(k T0) z-k
(1.5)
k 0
1.4.1. A Z- transzformáció néhány alapvető tulajdonsága: Linearitási (szuperpozíció) tétel (1.6): Ha Z[x1(t)] = F1(z) és Z[x2(t)] = F2(z), akkor Z[x1(t)+ x2(t)] = F1(z) + F2(z) és Z[x1(t) - x2(t)] = F1(z) - F2(z)
(1.6)
Azaz az összegfüggvény Z- transzformáltja az összetevők Z- transzformáltjának az összege, valamint (1.7): Z[ax(t)] = aF(z)
(1.7)
Ha Z[x(t)] = F(z), akkor Z[x(t-dT0)] = F(z)z-d
(1.8)
Eltolási tétel:
Ahol d a diszkrét holtidő, mely a folytonos Th holtidő és a T0 mintavételezési idő hányadosa, azaz d=Th/T0, ahol d=0,1,2,…., aminek a jelentése az, hogy időtartományban d mintavételezési idővel történő késleltetés a Z- tartományban z-d tényezővel való szorzást jelent (1.8). Szorzás t-vel (1.9): d t * f (t ) = T0 z F ( z ) dz
(1.9)
Végértéktételek (1.10, 1.11):
lim f (t ) = lim t
z 1
z 1 F ( z) z
(1.10)
lim f (t ) = lim F ( z ) z
t 0
12
(1.11)
A mintavételezett rendszerek alapvető tulajdonságai
1.4.2. A legfontosabb időfüggvények Laplace- és Z- transzformáltjai (1.1. táblázat):
y(t)
Y(s)
Y(z)
1
1 s
z z 1
T
1 s2
T0 z ( z 1) 2
t2
2 s3
T0 z(z 1) ( z 1)3
e-at
1 sa
z z e aT0
te-at
1 ( s a) 2
T0ze aT0 ( z e aT0 )2
1-e-at
a s( s a)
1 - e -aT0 ( z 1)( z e aT0 )
sint
1 s 2
z sinT0 z 2 z cos T0 1
s s 2
z (z - cosT0 ) z 2 2 cos T0 1
2
2
cost
2
2
e-at sint
z e -aT0 sinT0 z 2 2 ze -aT0 cos T0 e -2aT0
( s a) 2 2
1.1. táblázat: Időfüggvények Laplace- és Z- transzformáltjai Most nézzünk néhány példát függvények Z- transzformáltjának a meghatározására. A példák megoldása során segítségünkre lehet a 1.1. táblázat. Példa: Határozzuk meg az egyik legegyszerűbb, ám a gyakorlatban mégis nagy jelentőséggel bíró függvény, az egységugrás függvény (1.7. ábra) Z- transzformáltját.
13
A mintavételezett rendszerek alapvető tulajdonságai
1.7. ábra. Az egységugrás függvény mintavételezése
x(t) = 1(t) , ahol x(t)=1, ha t>0 x(t)=0,ha t<0. x*(t) =
x(kT )(t kT ) , ahol n 0
0
x(kT0) = 1 és T0 a mintavételezési idő.
0
A (1.5) egyenlet szerint az X(z) a következő formába írható (1.12):
X(z) =
x(kT )z n 0
k
0
= 1 + z-1 +z-2 + …
(1.12)
A kifejezés (1.12) tehát egy geometriai sor, amely a következő összegképlettel kifejezhető (1.13): X ( z)
1 z 1 z 1 1 z
(1.13)
Példa: Következő lépésben nézzük meg egy bonyolultabb függvénynek az e-at függvénynek a Ztranszformáltját. Most az x*(t) a következőképpen néz ki:
x*(t) =
e n 0
akT0
( t kT0 ) ,ahol T0 a mintavételezési idő és a egy konstans.
Ismét használva a (1.5) képletet a következő kifejezéshez jutunk (1.14):
X(z) =
x(kT0 )z k = X(z) = n 0
e
akT0
z n = 1 e aT z 1 e 2aT z 2 ... 0
k 0
Meghatározva az összegfüggvényt az alábbi eredményre jutunk (1.15):
14
0
(1.14)
A mintavételezett rendszerek alapvető tulajdonságai
X ( z)
1 1 e
aT0
z
1
z z e aT0
(1.15)
Példa: A következő példában nem egy időfüggvénynek, hanem egy véletlen számsorozatnak a Ztranszformáltját határozzuk meg. A sorozat a következő: [x(nT0)] = 3, 1, 0, -2 Most újra a (1.5.) függvényre lesz szükségünk. Ennek alapján írhatjuk a következő formát:
X(z) =
x(nT0 )z n = 3 1z 1 0z 2 (2) z 3 = 3 + n 0
3z 3 z 2 2 1 2 = z z3 z3
1.4.3. Az inverz Z- transzformáció A Z- transzformáció inverz művelete az X(z) transzformált függvényből határozza meg az x*(t) impulzussorozatot. A gyakorlatban általában X(z) z-nek racionális törtfüggvénye, így az inverz Z- transzformáció a legegyszerűbben a következő két módszerrel határozható meg (1.16, 1.17): 1. A törtnek z negatív hatványai szerinti sorbafejtése, amely a számlálónak a nevezővel való osztásával oldható meg, z-1 hatványainak együtthatói közvetlenül megadják a mintavételezett jel értékeit a mintavételi pontokban.
(b0 z n b1 z n1 ... bn ) : ( z n a1 z n1 ... an ) c0 c1 z 1 c2 z 2 ...
(1.16)
2. Részlettörtekre bontással: X(z) a következő formába írható: X(z) =
Az Bz Cz Dz 2 ... 2 ( z 1) ( z 1) ( z a ) ( z cz d )
(1.17)
Ezek után az egyes törtrészeknek megfelelő időfüggvényeket visszakeressük a Ztranszformációs táblázatból, és megkapjuk a megfelelő időfüggvényt, ami jelen esetben a következő lenne (1.18): x(t) = A
Bt Ce at Deat sin t ... , ahol t = kT0 T0 15
(1.18)
A mintavételezett rendszerek alapvető tulajdonságai
és T0 a mintavételezési idő. (Természetesen a MATLAB lehetőséget ad arra, hogy az inverz transzformációkat könnyedén elvégezzük.)
1.5. Az impulzus-átviteli függvény Az irányítási rendszerek analízisében és szintézisében nagy jelentőséggel bírnak azok a függvények, amelyek a be- és kimenő jelek közötti kapcsolatot írják le operátor tartományban. Folytonos esetben a kimenő és a bemenő jelek Laplace- transzformáltjának a hányadosát nevezzük átviteli függvénynek. Diszkrét esetben is létezik olyan függvény, amely leírja az adott rendszer bemenete és kimenete közti kapcsolatot. Ez a függvény nem más, mint a rendszer impulzus-átviteli függvénye, vagy diszkrét átviteli függvénye, amely a kimenő- és bemenő jel Ztranszformáltjának a hányadosaként írható fel (1.8. ábra).
1.8. ábra. A mintavételezett rendszer modellje
A (1.8. ábra) modellben (1.19):
y(z) b 0 z n b1z n 1 ... b n G 0 (z) m u (z) z a 1 z m1 ... a m
(1.19)
vagy (filter alak) (1.20):
G 0 (z 1 )
y(z 1 ) b 0 b1z 1 ... b n z n d B(z 1 ) d z z u (z 1 ) 1 a 1z 1 ... a m z m A(z 1 )
(1.20)
Ahol B és A polinomok együtthatói állandó számok, az n és m pedig pozitív egész számok lehetnek. Tehát a tervezés folyamán az irányított objektumok matematikai modelljét
16
A mintavételezett rendszerek alapvető tulajdonságai
az impulzus-átviteli függvény segítségével adhatjuk meg, amely a mintavételezési időpontban írja le a modell dinamikai viselkedését. A modellben szereplő d- a diszkrét holtidő. A fenti modellekben valós fizikai rendszerek esetén b0=0, ami azt jelenti, hogy a bemenő jelváltozáshoz képest a kimenő jelváltozás egy mintavételezési idővel késik. Erre egyszerű példa szolgálhat. Példa: Legyen a folytonos rendszer differenciál egyenlete (1.21) dy u, ahol u=1(t), y(0)=0 dt
(1.21)
Legyen a mintavételezési idő T0=1s. Diszkretizálva a fenti differenciál egyenletet, a következő differencia egyenlethez jutunk (1.22):
y(k ) y(k 1) u (k ) T0
(1.22)
Figyelembe véve azt, hogy T0=1s, a differencia egyenletet felírhatjuk a következő formában (1.23): y(k)=y(k-1)+u(k)
(1.23)
a rendszer impulzus-átviteli függvénye (1.24): G(z 1 )
y( z ) 1 u (z) 1 z 1
A differencia egyenlet és a differenciálegyenlet megoldása különbözik (1.9. ábra).
17
(1.24)
A mintavételezett rendszerek alapvető tulajdonságai
1.9. ábra. A differenciál és differenciaegyenletek megoldásának összehasonlítása
A fenti esetben b0=1. Nézzük azt az esetet, amikor b0=0. Ebben az esetben a differencia egyenlet a következő lesz (1.25): y(k)=y(k-1)+u(k-1)
(1.25)
az impulzus-átviteli függvény (1.26):
G(z 1 )
y( z ) z 1 u (z) 1 z 1
(1.26)
és a differenciál és differencia egyenletek megoldásai egybeesnek a mintavételezési időpontokban (1.10. ábra).
18
A mintavételezett rendszerek alapvető tulajdonságai
1.10. ábra. A differenciál és differencia egyenletek megoldásának összehasonlítása
Az impulzus-átviteli függvényből meghatározható a vizsgált szakasz erősítési tényezője. Vegyünk példának egy objektumot, amely az alábbi impulzus-átviteli függvénnyel írható le:
0.8z 1 G p (z ) 1 0,4z 1 1
Az erősítési tényező a következő összefüggés alapján számítható (1.27): n
K0
b i 0 m
i
a i 0
(1.27) i
Tehát az erősítési tényező:
K0
b 0 b1 0.8 0,8 1,333.... a 0 a 1 1 0.4 0,6
1.5.1. Zárt digitális szabályozási körök impulzus-átviteli függvényei: A digitális szabályozási rendszerek általában a 1.11. ábra szerinti hatásvázlat alapján épülnek fel. 19
A mintavételezett rendszerek alapvető tulajdonságai
1.11. ábra. Digitális szabályozási kör felépítése
A zárt rendszer folytonos és szakaszos komponensekből tevődik össze. Az A/D átalakítót, a D/A átalakítót, a mintavételező elemet, a nulladrendű tartószervet tartalmazó rendszert a G 0 (z 1 ) impulzus átviteli függvénnyel tudjuk leírni (1.28):
G 0 (z 1 )
b1z 1 ... b n z n B(z 1 ) d z z d 1 1 m A( z ) 1 a 1z ... a m z
(1.28)
A G0 objektum bemenete az u(k) impulzussorozat kimenete pedig az y(k) impulzussorozat. A z-d komponens pedig a rendszer diszkrét holdideje. A szabályozási algoritmust a digitális számítógépen helyezzük el, amely a
G C (z 1 )
u (z) Q(z 1 ) q 0 q1z 1 ... q n z e(z) P(z 1 ) 1 p1z 1 ... p m z
(1.29)
impulzus-átviteli függvénnyel jellemezhető (1.29), melynek a bemenete r(k) – y(k) = e(k) hibajel sorozat. A zárt mintavételezett szabályozási kör átviteli függvénye a következőképpen írható (1.30):
Gr ( z 1 )
y( z 1 ) G 0( z 1 )Gc( z 1 ) r ( z 1 ) 1 G0 ( z 1 )GC ( z 1 ) (1.30)
vagy behelyettesítve a megfelelő polinomokat (1.31):
G r (z 1 )
Q(z 1 )B(z 1 )z d P(z 1 )A(z 1 ) Q(z 1 )B(z 1 )z d
20
(1.31)
A mintavételezett rendszerek alapvető tulajdonságai
Mivel a valóságban különböző zavaró jelek érhetik a rendszer, a tervezésnél ezeket a zavaró jeleket is figyelembe kell venni. A fenti szabályozási kört ki kell bővíteni a kimenetre szuperponált zavaró jel impulzus-átviteli függvényével, valamint az objektum bemenetére ható zavaró jelek figyelembevételével is. A zavarok figyelembevételével a zárt szabályozási körünk hatásvázlata a következőképpen alakul (1.12. ábra):
1.12. ábra. Zajcsatornát tartalmazó egyhurkos digitális szabályozó kör
Abban az esetben, ha az irányított objektum bemenetére uv(k) zavarójel szuperponálódik, a zárt rendszer impulzus-átviteli függvénye r(k)=0, v(k)=0 esetén (1.32):
G u (z 1 )
G o (z 1 ) y(z 1 ) P(z 1 )B(z 1 )z d u v (z 1 ) 1 G o (z 1 )G C (z 1 ) P(z 1 )A(z 1 ) Q(z 1 )B(z 1 )z d (1.32)
ha pedig a szakasz kimenetére n(k) jel szuperponálódik és feltesszük, hogy r(k)=0 és uv(k)=0 bármely k>0-ra (1.33, 1.34), akkor:
y(z 1 ) 1 P(z 1 )A(z 1 ) G n (z ) n (z 1 ) 1 G o (z 1 )G C (z 1 ) P(z 1 )A(z 1 ) Q(z 1 )B(z 1 )z d 1
(1.33)
és
y(z 1 ) A(z 1 )P(z 1 )D(z 1 ) G v (z ) v(z 1 ) A(z 1 )P(z 1 )C(z 1 ) Q(z 1 )C(z 1 )B(z 1 )z d 1
(1.34)
ahol a zavarójel csatorna impulzus-átviteli függvénye pl. ARMAX modell esetén (ld. 2.1.3 fejezet) (1.35):
21
A mintavételezett rendszerek alapvető tulajdonságai
D(z 1 ) 1 d1z 1 ... d m z m G ov (z ) C(z 1 ) 1 c1z 1 ... c m z m 1
(1.35)
1.6. A mintavételezett rendszerek leírása differencia egyenletekkel: A
lineáris,
koncentrált
paraméterű
digitális
rendszerek
dinamikai
viselkedésének
vizsgálatához az időtartományban differencia egyenleteket alkalmazunk, melyek az y(k) kimenőjel és az u(k) bemenőjel közötti kapcsolatot írják le. Ahhoz, hogy felírhassuk egy rendszer differenciaegyenletét, definiáljunk egy ún. eltolási operátort, melynek jele: E és a következő formában használható (1.36):
E n [ y(k )] y(k n)
(1.36)
E n [ y(k )] y(k n)
(1.37)
vagy az inverzét képezve (1.37):
ahol n az eltolás mértékét jelzi, értéke bármilyen pozitív egész szám lehet. Egy impulzus-átviteli függvényt könnyedén felírhatunk differencia egyenlet formájában hiszen az impulzus-átviteli függvényben szereplő z operátor kitevője is eltolást jelent. Vegyük példának a következő impulzus-átviteli függvényt (1.38): (
)
(
)
(
)
(1.38)
Ezt a függvényt differenciaegyenlet alakba írva az
y(k ) a1 y(k 1) ... y(k n) b0 u(k ) b1u(k 1) ... bn u(k n) egyenletet kapjuk (1.39). Példa: Legyen az impulzus-átviteli függvényünk a következő: (T0=4s) G(z 1 )
0,1387z 1 0,0889z 2 1 1,036z 1 0,2636z 2
mely az alábbi folytonos rendszer diszkrét alakja:
22
(1.39)
A mintavételezett rendszerek alapvető tulajdonságai
G(s)
1 37,5s 12,5s 1
A fenti rendszer differenciál egyenlete: 37,5
dy 2 dy 12,5 1 u, 2 dt dt
differencia egyenlete: y(k)=1,036y(k-1)-0,2636y(k-1)+0,1387u(k-1)+0,0889u(k-2). A differenciál illetve a differencia egyenlet megoldását a 1.13. ábra szemléleti.
1.13. ábra. A differenciál ill. differencia egyenletek megoldásai
23
Dinamikus rendszerek identifikációja
2. DINAMIKUS RENDSZEREK INDENTIFIKÁCIÓJA
Szabályozási rendszerek tervezésének megvalósításához szükségünk van a szabályozott objektum dinamikai tulajdonságainak ismeretére. Valamilyen leírási mód segítségével rendelkeznünk kell az objektum matematikai modelljével, hisz a tervezőrendszerünk kérni fogja a szabályozandó objektum matematikai modelljét. Szerencsés esetben az objektum matematikai modellje – fizikai, kémiai, biológiai paraméterekkel ismert, de a legtöbb esetben nem, így valamilyen közelítő munkaponti modellekkel kell a tervezést megkezdeni. Egy teljesen új technológiai rendszer felépítésénél a beruházási terv természetesen tartalmazza a gépészeti, technológiai, műszerezési terveket, azonban a szabályozási rendszerek tervezését nem, illetve csak a szabályozási körök elképzelt struktúráját adják meg a tervezők. A számítógépben elhelyezendő szabályozási algoritmusok struktúrája és paraméterei ezekben a kiviteli tervekben csak közelítőleg szerepelnek, hisz a tervező a visszacsatolt rendszerek pontos dinamikai tulajdonságait nem ismerik, ezeket a megépített rendszeren végzett vizsgálatokkal lehet csak pontosítani. A technológiai rendszer felépítésével, a műszerezés megvalósításával egy próbaüzemeltetés kapcsán elvégezhető a technológiai rendszer objektumainak ún. identifikációja, azaz meghatározhatjuk azt a matematikai modellt, amely az adott munkapontban visszatükrözi a rendszer időbeli viselkedését. Ezen modell bázisán a tervező képes lesz a visszacsatolás, előrecsatolás algoritmusainak megtervezésére. A technológiai rendszerek irányítását, szabályozását, vezérlését folyamatirányító számítógépek végzik. A számítógép számára tervezendő szabályozási algoritmusokat a folyamatirányító számítógép differencia egyenlet formájában kéri, valamilyen tervezett T0 mintavételezési idő mellett, azaz az algoritmus egy diszkrét algoritmus, így a számítógép is egy diszkrét rendszert „szeretne látni” szabályozott objektumként, amely a valóságban természetesen időben folytonos rendszer. Az időben folytonos rendszerek identifikációjával a szakirodalom több évtizede foglalkozik [1] és letisztultak azok az identifikációs eljárások, amelyeket a gyakorlat is közvetlenül alkalmazhat. Kialakultak azok a diszkrét rendszermodellek, melyek a gyakorlati problémák legtöbbjének leírására megfelelnek, melyek felhasználásával olyan matematikai modelleket kapunk, amelyek közvetlenül felhasználhatók digitális szabályozási algoritmusok
24
Dinamikus rendszerek identifikációja
tervezéséhez, adaptív rendszerek megalkotásához (pld. identifikációt alkalmazó adaptív rendszerek [2]).
2.1. Tipikus diszkrét idejű rendszermodellek
A tipikus diszkrét idejű modelleket tárgyalva a jelölésrendszerünkkel figyelembe vesszük a MATLAB/System IdentificationToolbox jelölésrendszerét, amelyet a következő fejezetekben másként fogunk használni. A két jelölésrendszer megfelelése a következő (2.1):
Polinomok:
A( z 1 ) 1 a1 z 1 .... am z m helyett A( z 1 ) 1 a1 z 1 .... ana z na B( z 1 ) b1 z 1 ... bm z m helyett B( z 1 ) b1 b2 z 1 ... bnb z nb 1 C ( z 1 ) 1 c1 z 1 ... cm z m helyett
(2.1)
C ( z 1 ) 1 c1 z 1 ... cnc z nc D( z 1 ) 1 d1 z 1 ... d m z m helyett D( z 1 ) 1 d1 z 1 ... d nd z z d
Diszkrét holtidő d helyett nk-t fogunk használni a MATLAB reprezentációkban. Még egy
fontos
kitétel
az,
hogy
a
zajcsatorna
tárgyalásánál
a
MATLAB
System
IndentificationToolbox szerzője [www.mathworks.com ] a
Gov ( z )
n( z ) C ( z 1 ) v( z ) D( z 1 )
(2.2)
átviteli függvényt (2.2) használja, az algoritmus tervezési fejezetükben azonban Gov(z) csatorna átviteli függvényét a
Gov ( z )
n( z ) D( z 1 ) v( z ) C ( z 1 )
átviteli függvény (2.3) formában fogjuk használni.
25
(2.3)
Dinamikus rendszerek identifikációja
2.1.1. OE modell Az OE (Output Error) modell e(k) mérési zajt tartalmazó modell a következő hatásvázlattal rendelkezik (2.1. ábra):
2.1. ábra. Az OE modell struktúrája
melynek impulzus-átviteli függvénye (2.4)
y( z)
B( z 1 ) d z u ( z ) e( z ), A( z 1 )
(2.4)
ahol e(z) a mérési zaj z-transzformáltja. Az impulzus-átviteli függvényt sok esetben a következő formális felírási móddal is írják (2.5): y (k )
B( z ) u (k d ) e(k ), vagy A( z )
(2.5)
a MATLAB szerzői pedig a következő formát használják (2.6): y (k )
B( z ) u (k nk ) e(k ), azaz A( z )
(2.6)
a d diszkrét holtidő jelölése nk. Ezen modell olyan stacionárius folyamatok leírására szolgál, melyek külső zavaró jelcsatornákat nem tartalmaznak, azonban az e(t) mérési zaj jelen van. A modell differencia egyenlete MATLAB jelölés szerint a következő (2.7):
y(k ) a1 y (k 1) ... ana y(k na) b1u (k nk ) b2u (k nk 1) ... bnbu (k nk nb 1) e(k )
26
(2.7)
Dinamikus rendszerek identifikációja
2.1.2. ARX modell Az ARX modell (AutoRegressive with eXternal input) a kimeneten színes zajt feltételez, aminek a következő modell felel meg (2.2. ábra). A(z)y(z) = B(z) z d + e(z)
(2.8)
2.2. ábra. Az ARX modell struktúrája 2.1.3. ARMAX modell Az ARMAX modell (AutoRegressive Moving Average modell with eXogenous signal) struktúrája a következő (2.3. ábra):
2.3. ábra. Az ARMAX modell struktúrája
ahol a zajcsatorna és az alapcsatorna impulzus-átviteli függvénye (2.9)
Gov ( z )
C ( z 1 ) , D( z 1 )
B( z 1 ) d G0 ( z ) z A( z 1 ) a megfelelő polinomok (2.10): 27
(2.9)
Dinamikus rendszerek identifikációja
C ( z 1 ) 1 c1 z 1 ... cuc z nc D( z 1 ) A( z 1 ) 1 a1 z 1 ... ana z na B( z 1 ) b1 b2 z 1 ... bnb z nb 1
(2.10)
Az ARMAX modell olyan esetben alkalmazható, ha a zajmodell – Gov – és az alapcsatorna karakterisztikus egyenlete megegyezik. Ezt az egyezést a rendszer részletes analízisével tudjuk meghatározni. Nézzük a következő példát. Legyen az identifikáló rendszer egy hőcserélő berendezés (2.4. ábra).
2.4. ábra. A hőcserélő – az identifikáció objektuma A gőzzel fűtött hőcserélő feladata az, hogy a kilépő folyadék Tki hőmérsékletét adott értékre állítsa egy számítógépes digitális szabályozási algoritmus segítségével. A TT hőmérséklet távadó mA-es kimenete y(t). A gőzágat több fogyasztó is terhelheti, ezért a P(t) gőznyomást stabilizáljuk. A gőzágba építik be a szabályozó szelepet, amely a mi esetünkben pneumatikus segédenergiával működik. A számítógép folyamatperifériájánál kapott mA-es u(t)
végrehajtó
jel
egy
elektropneumatikus
jelátalakítón
keresztül
működteti
a
végrehajtószervet. Tételezzük fel, hogy a belépő folyadék mennyisége és hőmérséklete állandó és az y(t)-re csak a mérési zaj szuperponálódik. Ebben az esetben az objektum jó közelítéssel modellezhető egy OE struktúrával (2.1. ábra). 28
Dinamikus rendszerek identifikációja
Ha gőz gerincvezetékének nyomását nem tudjuk stabilizálni és az egyéb fogyasztók véletlenszerűen terhelik a gőzrendszert, akkor a P(t) gőznyomás ingadozik, amely a kilépő y(t) hőmérséklet véletlenszerű változását vonja maga után. A hőmérséklet aktuális értékét a szelepre eső nyomásesés (bevitt hőenergia) határozza meg, mely nyomásesést vagy a gőznyomásváltozás, vagy a szabályozó szelep pozíciója befolyásolják. Így feltételezhető, hogy a zajcsatorna (gőznyomás-hőmérséklet) és az alapcsatorna (szeleppozíció-hőmérséklet) karakterisztikus egyenletei megegyeznek, azaz A(z-1)=D(z-1). Ilyen esetben az objektum munkapontban egy ARMAX modellel közelíthető. A harmadik esetként tételezzük fel, hogy a P(t) gőznyomás stabilizált, de a hőcserélőn átmenő Q(t) mennyiség sztochasztikusan változik, amely a kilépő y(t) hőmérsékletet befolyásolja. Ebben az esetben az alapcsatorna és a zavarójel csatorna (belépő folyadékáram hőmérséklet) karakterisztikus egyenletei nem fognak megegyezni, így a fent említett modellek helyett Box-Jenkins modellt ajánlott használni.
2.1.4. Box-Jenkins (BJ) modell A BJ modell a következő hatásvázlattal rendelkezik (2.5. ábra):
2.5. ábra. A BJ modell struktúrája
A BJ struktúra impulzus-átviteli függvénye (2.11): ( )
(
)
(
)
( )
(
)
(
)
( )
(2.11)
B( z 1 ) C ( z 1 ) u ( k nk ) e(k ), A( z 1 ) D( z 1 )
(2.12)
vagy MATLAB reprezentációja (2.12):
y (k ) ahol:
29
Dinamikus rendszerek identifikációja
D( z 1 ) 1 d1 z 1 ... dnd q nd
(2.13)
A bemutatott OE, ARX, ARMAX, BJ modellek azok, amelyeket a gyakorlatban legtöbbször alkalmazni kell, azonban a MATLAB kezelni tudja az AR, PEM modelleket (ld. MATLAB kézi könyv) is. A MATLAB konvenciója szerint a polinomokat együtthatóikkal kell megadni a z hatványainak csökkenő sorrendjében, tehát A(z) esetén az A 1 a1 a2 ... ana
sorvektor formájában. Speciális megkötés, hogy nk és B egyszerre kerül megadásra oly módon, hogy B elején nkdarab vezető nullát szúrunk be. Pl. mintavételezett kéttárolós tag és ARX modell esetén a diszkrét idejű rendszermodell (2.14): A( z) y(k ) B( z)u(k ) e(k ),
(2.14)
melyet az A 1 a1 a2 B 0 b1 b2
sorvektorokkal jellemezhetünk [3]. A rendszermodell speciális „th” adatstuktúrával jellemezhető. Ha a rendszer ismert (a polinomok együtthatói és a fehér zaj szórásnégyzete numerikusan ismertek, mert például ismert rendszeren akarunk jeleket előállítani szimulációs vizsgálatokhoz, pl. az identifikációs módszerek ellenőrzése céljából egy etalon rendszeren), akkor az IDENT a th poly 2th( A, B, C, F , , T )
függvényhívás hatására felépíti és th-ban tárolja a rendszert jellemző adatstruktúrát. A hátul álló C, D, F , , T paraméterek elhagyhatók híváskor, ilyenkor értékük 1 lesz. Szimulációs vizsgálatokhoz ezután zajos vagy zajmentes kimenetet generálhatunk ismert u és e bemenő jel és zaj sorozatok esetén az IDENT
y idsim(u e , th) y idsim(u, th) függvényhívásaival. A bemenő jelek sorozatait oszlopvektorok formájában kell megadni, a kimenő jel oszlopvektorokban keletkezik. Bemenő jeleket generálhatunk például a MATLAB rand és sign függvényei felhasználásával.
30
Dinamikus rendszerek identifikációja
Ismeretlen rendszer identifikációja esetén mérjük az u(t) bemenet és a zajos y(t) kimenet értékeit, ezeket összegyűjtjük az u és y oszlopvektorokban, felépítjük a megfigyelésekből álló és két oszloppal jelentkező z=[y u] mátrixot, megválasztjuk a rendszermodellt és az abban szereplő polinomok fokszámait, megválasztjuk a holtidő értékét és elvégezzük az identifikációt az IDENT szolgáltatásaival. Az identifikációs módszerek feladata, hogy meghatározza a polinomok valamilyen értelemben optimális paramétereit, miközben az e(t) zaj értékei ismeretlenek. Erre a célra a legkisebb négyzetek módszere (LS, leastsquaresmethod) vagy ennél általánosabb paraméterbecslési technikák, a zaj fehérítését megcélzó segédváltozós (I.V., instrumentalvariables) módszer, numerikus optimum-keresés, vagy esetleg ezek kombinációja alkalmazható. Az IDENT a rendszermodell típusától függően választja meg az identifikációhoz használt numerikus módszert. A segédváltozós módszer (IV) csak ARX modell esetén alkalmazható [3] Megjegyezzük,
hogy
több
itt
szereplő
függvénynek
létezik
általánosabb
paraméterezése is IDENT-ben, továbbá többváltozós (MIMO) rendszerek identifikációja is lehetséges. Az IDENT toolbox a továbbiakban algoritmus szinten is részletesen tárgyalt diszkrétidejű paraméter identifikációs módszereken kívül lehetővé teszi még folytonos idejű lineáris állapotegyenletek paramétereinek becslését is speciális zajstruktúrák esetén. Szolgáltatásai között szerepelnek a paraméterbecslés rekurzív realizációi is , amelyek adaptív irányításoknál jelentősek. Ezeken túlmenően lehetőség van nemparaméteres identifikációra, valamint a korrelációs függvények és a spektrumok számítására is. Ez utóbbi számítások algoritmusai gyors Fourier-transzformáción (FFT, Fast-Fourier-Transformation) és az FFT periodicitása miatt alkalmasan választott ablakozási technikán alapulnak [4]
2.2. Az identifikáció gyakorlati alkalmazásai Az identifikáció előkészítése nagyon körültekintő munkát igényel. A következőkben megnézzük a legfontosabb – identifikáció előtti – feladatokat, melyeket el kell végezni. 2.2.1. Általános feladatok a) Az állandósult állapot vizsgálata
31
Dinamikus rendszerek identifikációja
Nagyon gyakran a technológiai rendszerek különböző munkapontokon dolgoznak. Ezek a technológiai rendszerek nemlineárisak, pld. az erősítési tényezőjük jelentősen változhat a különböző működési tartományokban. Az erősítési tényező mellett változhatnak az időállandók, változhatnak egyéb modell paraméterek is, vagy a modell struktúra is más lesz. A vizsgált identifikációs eljárások eredménye egy impulzus átviteli függvény állandó paraméterekkel és struktúrával. Ezért fontos, hogy elérjünk egy állandósult állapotot (2.6. ábra) – és a bemenet-kimenet pont párokat innen gyűjtsük össze.
2.6. ábra. Az állandósult állapot vizsgálata b)Trend figyelés A gyakorlatban előfordulhat az átlagérték lassú kúszása pld. a távadók lassú driftje miatt, az objektum elemeinek kopása, vagy az aktuális hőcserélő vízkövesedése folytán, esetleg műszakváltás, időközi leállások miatt. Ebben az esetben a kúszást (2.7. ábra) trendfigyeléssel meg kell szüntetni, hisz a nemlineáris objektum tulajdonságok az identifikáció eredményét meghamisíthatják.
32
Dinamikus rendszerek identifikációja
2.7. ábra. Trend figyelés c) Kiugró értékek figyelése Az
állandósult
állapot
megfigyelésekor
találkozhatunk
szokatlan
kiugró
értékek
megjelenésével a mért jeleken (2.8. ábra). Ennek több oka lehet: a távadó vagy a kábelezés érzékeny
az
elektromágneses
tér
változásaira,
amikor
be/kikapcsolnak
villamos
berendezéseket, villámlás történt, stb. Az ilyen adatok nem alkalmasak az identifikációhoz.
2.8. ábra. Kiugró értékek figyelése
33
Dinamikus rendszerek identifikációja
d) Az érzékelő kikapcsol (2.9. ábra)
2.9. ábra. Az érzékelő kikapcsol e) A mintavételezési idő meghatározása A To mintavételezési idő meghatározása abban az esetben, ha ismerjük a vizsgált SISO objektum átmeneti függvényét az 2.10. ábra alapján történik. Az inflexiós pontban húzott érintő segítségével meghatározható a To mintavételezési idő. Abban az esetben, ha egy működő technológia (2.11. ábra) jeleit akarjuk mintázni, érdemes a mért jeleket kinagyítani (2.12. ábra) és kiválasztani azt a legnagyobb frekvenciát, mely még befolyásolja a modellalkotást. A To mintavételezési idő a legnagyobb frekvenciájú összetevő periódusidejének a húszad része legyen.
2.10. ábra. A mintavételezési idő meghatározása 34
Dinamikus rendszerek identifikációja
2.11. ábra. A legnagyobb frekvenciájú rész kiválasztása
2.12. ábra. A mintavételezési idő meghatározása 2.2.2. MATLAB System IdentificationToolbox és SIMULINK alkalmazása Az identifikációs technikák alkalmazásának tanulmányozására a MATLAB/SIMULINK modellező rendszert fogjuk használni és az identifikálandó objektumot is szimulációval állítjuk elő, természetesen folytonos rendszerként modellezzük azt. A 2.13. ábra segítségével egy másodrendű rendszer átmeneti függvénye alapján meghatározhatjuk a diszkrét modell mintavételezési idejét. A T0 mintavételezési időt lengésmentes önbeálló rendszer esetén kiválaszthatjuk úgy, hogy megkeressük az egy időállandós, holtidős közelítő modelljét és ezen modell T időállandójának tizedét fogadjuk el mintavételezési időnek. Legyen a folytonos rendszer modellje: G0 ( s)
2, 7 , 140s 45s 1 2
35
Dinamikus rendszerek identifikációja
közelítő egy időállandós, holtidős modellje: 2, 7 e2 s Gok ( s) 51s 1
ahol az időállandó T=50s, így a mintavételezési időt választhatjuk T0=T/10=50/10=5 s-ra (2.13. ábra).
2.13. ábra. A folytonos rendszer átmeneti függvénye és közelítése egyidőállandós, holtidős rendszerrel A 2.14. ábra azt a SIMULINK felületet mutatja be, ahol az identifikálandó objektum u(t) bemenetét állítjuk elő, ahol u(t)=us(t)+u*(t). A munkaponti identifikáció első fontos feladata a vizsgálat munkapontjának meghatározása, majd beállítása. Gyakori eset, hogy a folytonos technológiai rendszer különböző terheléssel működik és a rendszer nemlineáris volta miatt több munkaponti identifikáció szükséges. Ezen munkapontok meghatározásánál a technológusok segítségét, szakértelmét kell felhasználni. Az identifikáció tesztjelének megválasztása is fontos kérdés. A vizsgáló tesztjel frekvencia spektrumát és amplitúdóját úgy kell megváltoztatni, hogy az u(t) tesztjel hatása az y(t) kimeneten mérhető és értékelhető legyen. Off-line identifikáció esetén a leggyakoribb tesztjel forma a sávkorlátozott fehér zaj (Band-Limited White Noise a SIMULINK könyvtárból), vagy az álvéletlen kétállapotú jelsorozat (Psendo Random BinareSignal, PRBS). Ezek a jelalakok tartalmazzák a legtöbb 36
Dinamikus rendszerek identifikációja
frekvenciát, a frekvencia spektrumuk a legszélesebb a mesterségesen előállított jelsorozatok között. A PRBS egyik előállítási módját látjuk a 2.15. ábrán, ahol a mérési zaj modellezésére használjuk. Az us(t) folytonos jel segítségével állítjuk be a rendszert a kívánt munkapontba, az u*(t) jel pedig az identifikáció bemeneti tesztjele.
2.14. ábra. Az identifikáció előkészítése SIMULINK felületen Az objektum folytonos modellje a kimenetre szuperponálódó mérési zajjal a 2.15. ábrán látható:
2.15. ábra. Az objektum folytonos modellje 37
Dinamikus rendszerek identifikációja
A 2.14. ábrán az u(t) bemenőjel és az y(t) kimenőjel adatgyűjtését T0=5s-os mintavételezési idővel a ToWorkspace blokkal végzik, mely paraméterezését mutatja a 2.16. ábra, s ennek segítségével mátrix formában mentjük el a gyűjtött adatot.
2.16. ábra. Az adatgyűjtő blokk Állítsuk elő a bemenet-kimenet párok vektorát, és vizsgáljuk meg mérési eredményeinket a következő parancsok segítségével (2.17. ábra): z= [ y u ] idplot( z )
2.17. ábra. A mintavételezett adatpárok
38
Dinamikus rendszerek identifikációja
Az identifikációs algoritmusok nemlineáris rendszerek esetében csak akkor hoznak helyes eredményeket, ha a bemenőjel és a kimenőjel átlagértéke állandó, ezért a mi esetünkben az első 50 mintavételezett adatpárt eltávolítjuk a z1=[ y(50:400)
u(50:400)]
parancs segítségével (2.18.ábra).
2.18. ábra. A stacionárius be- kimeneti adatpárok Az identifikációs algoritmusok másik megkötése lehet az, hogy a bemeneti-kimeneti adatpárok nulla átlagértékűek legyenek, ezt pedig a z2 = dtrend(z1) paranccsal végezhetjük el. (2.19. ábra).
2.19. ábra. Az identifikációs algoritmus bemenetének végleges adatpárjai 39
Dinamikus rendszerek identifikációja
A véglegesnek tekinthető adatpárok meghatározása után az identifikáció számára meg kell határozni az illesztendő modell (pl. impulzus-átviteli függvény) számlálójának (nb) és nevezőjének (na) fokszámait, valamint a d diszkrét holtidőt, melyet nk-val jelöl a használt modellező rendszer. Ha az adatpárokból ezeket a paramétereket nem vagyunk képesek meghatározni, ezért egy egyszerű paranccsal ir = cra(z2) , azaz korrelációs analízis segítségével meghatározhatjuk a rendszerközelítő súlyfüggvényét, melyet a 2.20. ábrán mutatunk be.
2.20. ábra. A vizsgált rendszer közelítő súlyfüggvénye A közelítő súlyfüggvény, mint a vizsgált rendszer minimum másodrendű modellel jellemezhető, diszkrét holtideje d=0, mely a MATLAB realizációban nk=1-et jelent. A következő parancsokkal megvizsgálhatjuk a vizsgált rendszerközelítő átmeneti függvényét is (2.21. ábra): stepr = cumsum (ir) plot (stepr)
40
Dinamikus rendszerek identifikációja
2.21. ábra. A vizsgált rendszerközelítő átviteli függvénye A korrelációs analízis eredménye alapján másodrendű, holtidős rendszert jósolhatunk, amely számára állítsuk elő a th (téta) vektort, mely tartalmazni fogja az identifikálandó OE modell B(z-1),A(z-1) polinomjainak együtthatóit:
=[ a0 a1a2 …ana b1 b2 … bnb] a
th = oe (z2, [2 2 1])
parancs segítségével.
A present(th) parancs hatására megkapjuk a B(z-1) és A(z-1) polinomok együtthatóit: B=
0
A = 1.0000
0.1488
0.0877
-1.1129
0.2005
majd elvégezhetjük a modellünk és a tényleges mérési adatsorunk által szolgáltatott eredményeket ugyanarra a bemeneti jelre: compare (z2, th). A vizsgált rendszer esetén a mérési zaj e(k)=0 amplitúdóval rendelkezett és a modell kimenete ugyanazon vizsgáló jel esetén megegyezik a fizikai objektum kimenetével (2.22. ábra).
41
Dinamikus rendszerek identifikációja
2.22. ábra. A modell és a fizikai rendszer kimenetének összehasonlítása (e(k)=0)
Az OE modell impulzus-átviteli függvénye:
G0 ( z )
0.1488 z 1 0, 0877 z 2 y( z) 1 2 1 1,1129 z 0, 2005 z u( z)
melynek pólusait és zérusait mutatja a 2.23. ábra. Ezek az e(k)=0 esetre megegyeznek a tényleges
értékekkel.
A fizikai
rendszer
és
a modell
összehasonlítását mutatjuk be a 2.24. ábrán.
2.23. ábra. A modell pólusai és zérusai
42
frekvencia függvényének
Dinamikus rendszerek identifikációja
2.24. ábra. A folytonos rendszer és az identifikált modell frekvencia függvénye
A modellező rendszer SIMULINK felülete is rendelkezik identifikációs számítási blokkokkal (2.25. ábra), melyet a 2.26. ábra szerint helyezhetünk be az identifikációs rendszerbe, és a 2.27. ábrán bemutatott paraméter beviteli felülettel rendelkezik. Az így beillesztett ARX vagy más blokk az identifikációs eredményeket a MATLAB ablakban jeleníti meg a következő alakban: Transferfunction: num / den
0,1488 z 0, 0877 z 1,11289 z 0, 20046 2
2.25. ábra. A SIMULINK identifikációs blokkjai 43
Dinamikus rendszerek identifikációja
2.26. ábra. Az ARX blokk behelyezése
2.27. ábra. Az ARX blokk paraméterezése A 2.28. és 2.29. ábrákon mutatjuk be az identifikáció összehasonlító eredményeit e(k)=0 és 0,2 amplitúdójú álvéletlen (PRBS) mérési zaj esetén. Az ábrák alapján levonhatjuk azt a következtetést, hogy az e(k) mérési zaj amplitúdója az identifikáció pontosságát nagyban befolyásolja, tehát törekednie kell a mérési zajok amplitúdójának csökkentésére pl. úgy, hogy az y(t) kimenetet, vagy a már mintázott y(k) kimenetet megszűrjük. Természetesen a szűrő dinamikája is az identifikált modellünkbe integrálódik.
44
Dinamikus rendszerek identifikációja
2.28. ábra. Az identifikáció minősége e(k)=0 esetén
2.29. ábra. Az identifikáció minősége mérési zaj esetén A 2.30. ábrán egy ARMAX modell identifikációját illusztráljuk. Az alapcsatorna folytonos átviteli függvénye:
4(8s 5)e6 s G0 ( s) (20s 1)(32s 1)(11s 1) a zajcsatorna folytonos átviteli függvénye: Gov ( s)
8 (20s 1)(32s 1)(11s 1)
45
Dinamikus rendszerek identifikációja
A mintavételezési időt, To=6 s-ra választottuk.
2.30.ábra. Az ARMAX modell SIMULINK felületen A SIMULINK ARMAX blokkjának paraméterezését láthatjuk a 2.31. ábrán.
2.31. ábra. Az identifikátor paraméterezése
Az identifikáció minőségét mutatjuk be a 2.32. ábrán.
46
Dinamikus rendszerek identifikációja
2.32. Az identifikáció minősége A MATLAB felületen hozzáférhetünk az identifikált alapcsatorna és a zajcsatorna impulzus-átviteli függvényeihez, melyek a következők lettek:
G0 ( z )
Gov ( z )
0,86633z 2 0,37593z 0, 23062 z 4 1,844 z 3 1, 0601z 2 0,18238 z
z 3 0,83757 z 2 0,34081z 0, 085343 z 3 1,844 z 2 1, 0601z 0,18238
A Box-Jenkins modell bemutatására a következő folytonos zajcsatorna modellt használtuk: Gov ( s)
8 , (10s 1)(24s 1)(8s 1)
az alapcsatorna G0(s) modellje megegyezik az ARMAX modellével. A SIMULINK felületet láthatjuk a 2.33. ábrán, valamint a BJ-modell paraméterezését a 2.34. ábrán.
47
Dinamikus rendszerek identifikációja
2.33. ábra. Box-Jenkins modell SIMULINK felületen
2.34. ábra. Az identifikátor paraméterezése A Box-Jenkins modell alapján történt identifikáció minőségét mutatja a 2.35. ábra.
48
Dinamikus rendszerek identifikációja
2.35. ábra. Az identifikáció minősége
2.2.3. Nemlineáris rendszerek identifikációja A technológiai rendszerünk egy gőzös hőcserélő (2.36. ábra), gőzágba épített nemlineáris karakterisztikájú szabályozó szeleppel. Feladatunk megtervezni a visszacsatolásba építendő szabályozási algoritmust.
2.36. ábra. Gőzös hőcserélő
49
Dinamikus rendszerek identifikációja
A technológiai rendszer nemlineáris viselkedését láthatjuk a 2.37. ábrán. A szelep bemenetét változtatva különböző munkapontokon vizsgáltuk a hőcserélőből kilépő folyadék hőmérséklet változásait. A mintavételezési időt To=2 s-ra választottuk és az u(k) végrehajtó jelre +/- 1 mA PRBS jelet szuperponáltunk.
2.37. ábra. Gőzös hőcserélő nemlineáris viselkedése Vizsgáljunk meg egy tartományt, melynél az u(k) bemenet 10 mA-ről 14 mA-re változik. Ez a k=300-400 minta közötti rész minta párjai. OE modellt használva a következő átviteli függvény polinomokat kapjuk: B(q) = 0.06827 q-1 + 0.1446 q-2 - 0.2139 q-3 F(q) = 1 - 2.306 q-1 + 1.775 q-2 - 0.4706 q-3 Az identifikációt elvégezve a k=200-300 intervallumban az átviteli függvény polinomjai: B(q) = 0.02024 q-1 + 0.06032 q-2 - 0.06966 q-3 F(q) = 1 - 2.179 q-1 + 1.623 q-2 - 0.4214 q-3 Ez azt jelenti, hogy több munkapont környékén el kell végezni az identifikációt és a más-más munkapont környékén a szabályozó algoritmus paramétereit újra kell hangolni. Az alábbi ábra a nemlineáris rendszer szabályozását szemlélteti állandó paraméterű és struktúrájú szabályozó esetén (2.38. ábra).
50
Dinamikus rendszerek identifikációja
2.38. ábra. Nemlineáris rendszer állandó paraméterű és struktúrájú szabályozó esetén (A MATLAB >>ident parancsára nem térünk ki, lásd a System Identification Toolbox-ot)
Irodalomjegyzék [1] Eykhoff P.: System identification. John Wiley and SonsLtd, 1974 [2] Isermann R.:Digital Control Systems. Springer-Verlag, 1989 [3] Lantos B.: Irányítási rendszerek elmélete és tervezése, Akadémiai Kiadó, 2001. [4] Lenart Ljung: System Identification. Theory for theUser.Prentice Hall 1987
51
Digitális PID szabályozási algoritmusok
3. DIGITÁLIS PID SZABÁLYOZÁSI ALGORITMUSOK
3.1. Folytonos PID algoritmusokból származtatott digitális PID algoritmusok A digitális PID algoritmusokat az analóg PID algoritmusokból is származtatjuk, az analóg PID szabályozók differenciálegyenletének differenciaegyenletté alakításával. Egy ideális folytonos PID szabályozó átviteli függvénye legyen (3.1): G c (s)
u (s) 1 K(1 TD s) e(s) sTI
(3.1)
ez alapján felírva a differenciál egyenletet a következő kifejezést kapjuk (3.2): ( )
( )
∫
( )
( )
(3.2)
ahol az egyes paraméterek a következőket jelentik: K erősítési tényező TI integrálási időállandó TD differenciálási időállandó Abban az esetben, ha a mintavételezési idő T0 kicsi, akkor a differenciálegyenlet differencia egyenletté alakítható. A differenciáló tag egy egyszerű kivonással, míg az integráló komponens összegzéssel helyettesíthető. Ezek figyelembevételével fel tudjuk írni a diszkretizált egyenletet (3.3):
u (k ) K [e(k )
T0 TI
k 1
TD
i 0
0
e(i) T
(e(k ) e(k 1))]
(3.3)
Ezt az algoritmust pozíció algoritmusnak nevezzük, mivel a kimenete a beavatkozó szerv mindenkori állapotának megfelelő érték.
Nézzük meg a sebességalgoritmus levezetését, amely a beavatkozó szerv mindenkori állapotának kívánt, előjelhelyes megváltozását adja. A sebességalgoritmus levezetéséhez meg kell néznünk a beavatkozó jel egy mintavételezési idővel korábbi értékét u(k-1)-et (3.4):
u (k 1) K [e(k 1)
T0 TI
k 2
TD
i 0
0
e(i) T 52
(e(k 1) e(k 2))]
(3.4)
Digitális PID szabályozási algoritmusok
A levezetés folytatásában vegyük az u(k) és u(k-1) értékek különbségét és végezzük el a lehetséges összevonásokat (3.5):
u (k ) u (k 1) K (1
T TD T T )e(k ) (1 2 D 0 )e(k 1) K D e(k 2) T0 T0 TI T0
(3.5)
A hibajelek együtthatóit q-val jelölve a sebesség algoritmus az (3.6) formába írható:
u(k ) u(k ) u(k 1) q0 e(k ) q1e(k 1) q2 e(k 2)
(3.6)
ahol az együtthatók értékei (3.7):
q0 (1
q1 (1 2
TD ) T0
TD T0 ) T0 TI
q2 K
TD T0
(3.7)
A sebességalgoritmus impulzus átviteli függvénye (3.8): G c (z 1 )
u ( z) q 0 q1z 1 q 2 z 2 e(z)
(3.8)
Sebességalgoritmusokat akkor alkalmazunk, ha a rendszer végrehajtó szerve integráló jellegű tulajdonsággal bír. Példaként nézzünk meg egy PID pozíció (3.1. ábra) és egy sebességalgoritmus (3.2. ábra) viselkedését (átmeneti függvényét) e(k)=1 esetén.
3.1. ábra. Pozíció algoritmus átmeneti függvénye 53
Digitális PID szabályozási algoritmusok
3.2. ábra. PID sebességalgoritmus átmeneti függvénye
3.2. Általános digitális szabályozási algoritmusok impulzus-átviteli függvényei: Nagy mintavételezési idő esetében a folytonos algoritmusok diszkretizálása már nem végezhető el, mert a közelítések nagyon pontatlanok lesznek. Természetesen nagyon sok esetben tudunk megfelelően kicsi mintavételezési időt választani, azonban előfordulhatnak olyan esetek is amikor ez már nem lehetséges. Ilyenkor a diszkretizálás helyett más módszerhez kell folyamodnunk. Nézzük meg az ismert egyhurkú szabályozási kört (3.3. ábra).
3.3. ábra. A digitális szabályozó kör hatásvázlata
Ebben a szabályozási körben a szabályozott objektum impulzus-átviteli függvénye (3.9):
54
Digitális PID szabályozási algoritmusok
b1z 1 ... b m z m d y(z 1 ) B(z 1 ) G 0 (z) z u (z 1 ) A(z 1 ) 1 a 1z 1 ... a m z m
(3.9)
ahol d a rendszer diszkrét holdideje, a mintavételezési idő pedig T0. Az általános lineáris digitális szabályozó GC impulzus-átviteli függvénye (3.10): G C ( z)
u (z) Q(z 1 ) q 0 q1z 1 ... q z e(z) P(z 1 ) p 0 p1z 1 ... p z
(3.10)
ahol q0 0 , p0 1 . A szabályozó algoritmus polinomjainak a fokszáma a szabályozott objektum polinomjának a fokszámától függ, általában abból meg is határozható és a gyakorlatban általában kisebb az-objektum fokszámánál, tehát m és m d . Önbeálló rendszerek esetén pozíció algoritmust alkalmazva a GC(z) átviteli függvényben z=1 egyszeres pólus szükséges ahhoz, hogy a zárt digitális szabályozási körben ne legyen maradó szabályozási eltérés (3.11):
G C ( z)
Q(z 1 ) P * (z 1 )(z 1)
(3.11)
A legegyszerűbb lineáris digitális szabályozási algoritmus z = 1 pólussal (3.12):
Q(z 1 ) q 0 q1 z 1 ... q z G C (z) P(z 1 ) 1 z 1
(3.12)
Attól függően, hogy a számláló hányadfokú, beszélhetünk különböző algoritmusokról. Abban az esetben, ha υ = 1 PI típusú, ha υ = 2 PID típusú algoritmusokhoz jutunk. A (3.12) pozíció algoritmusból származtatott differencia egyenlet (3.13):
u(k ) u(k 1) q0 e(k ) q1e(k 1) ... q e(k )
(3.13)
3.3. Első és másodrendű szabályozási algoritmusok Másodrendű szabályozási algoritmus impulzus-átviteli függvénye (3.14): q 0 q1 z 1 q 2 z 2 G C ( z) 1 z 1
55
(3.14)
Digitális PID szabályozási algoritmusok
differencia egyenlete (3.15):
u(k ) u(k 1) q0 e(k ) q1e(k 1) q2 e(k 2)
(3.15)
A másodrendű algoritmus átmeneti függvényét ábrázolva felfedezhetők a PID viselkedés jelei (3.4. ábra).
3.4. ábra. A diszkrét PID algoritmus átmeneti függvénye és a q paraméterek kapcsolata Abban az esetben, ha az u(1)
u(k-1) akkor a diszkrét idejű szabályozó nagyban hasonlít az általában használt folytonos idejű PID szabályozóra. Vezessük be a következő jelöléseket (3.16):
K q0 q 2
erősítési tényező
c D q2 / K
differenciálási időállandó
c I (q0 q1 q2 ) / K
integrálási időállandó
Ezek figyelembevételével az előbbi 3.4. ábra a következőképpen alakul át (3.5. ábra):
3.5. ábra. A diszkrét PID algoritmus és a K, cD, cI paraméterek kapcsolata
56
(3.16)
Digitális PID szabályozási algoritmusok
Ezek az együtthatók kis mintavételezési időnél a folyamatos idejű PID szabályozó paramétereiből (K’, TD’, TI’) származtathatók a diszkretizálást követően (3.17) egyenlet szerint:
K=K’ , c D
T TD , cI 0 T0 TI
(3.17)
Az általános PID algoritmus esetében a következő összefüggések állnak fenn (3.18): c D 0; c I 0;
(3.18)
c I c D;
Ha ezeket az együtthatók segítségével Z- transzformált alakba írjuk fel a szabályozó impulzus-átviteli függvényét, akkor a következő impulzus-átviteli függvényhez jutunk (3.19):
G C (z)
u (z) K[(1 c D ) (c I 2c D 1)z 1 c D z 2 ] z 1 K 1 c c D (1 z 1 ) I 1 1 e(z) 1 z 1 z (3.19)
Így a folytonos PID – szabályozó analógiájára a P- a D- és az I - csatorna már szétválasztható és külön - külön is megvalósítható. Differencia egyenlet alakba írva az egyes PID komponenseket leíró függvények (3.20): u p (k ) Ke(k )
(P- tag)
u I (k ) u I (k 1) Kc I e(k 1)
(I- tag)
u D (k ) Kc D e(k ) Kc D e(k 1)
(D- tag),
(3.20)
melyek a folyamatirányító számítógépen közvetlenül programozhatóak. Az eredő differencia egyenlet ennek a három egyenletnek az összege (3.21):
u (k ) u P (k ) u I (k ) u D (k )
57
(3.21)
Digitális PID szabályozási algoritmusok
Ennek alapján fel tudjuk rajzolni a PID algoritmus blokkvázlatát is (3.6. ábra).
3.6. ábra. A diszkrét PID algoritmus blokkvázlata
A diszkrét idejű PID pozíciójú algoritmusnak az impulzus-átviteli függvénye q paraméterekkel (3.22): q 0 q1 z 1 q 2 z 2 G C ( z) 1 z 1
(3.22)
Ahol a q paraméterekkel felírt differencia egyenletrendszer a következő (3.23): u p (k ) (q0 q2 )e(k )
u I (k ) u I (k 1) (q0 q1 q2 )e(k 1)
(3.23)
u D (k ) q2 (e(k ) e(k 1)) u(k ) u p (k ) ui (k ) ud (k )
A másodrendű szabályozási algoritmusból kaphatjuk meg az elsőrendű szabályozási algoritmust ha q2 = 0. Így a szabályozó impulzusátviteli függvénye a következőképpen alakul (3.24): G C ( z)
q 0 q1 z 1 1 z 1
(3.24)
Ez alapján a differencia egyenlet (3.25):
u(k ) u(k 1) q0 e(k ) q1e(k 1) 58
(3.25)
Digitális PID szabályozási algoritmusok
mely egy digitális PI algoritmus egyenlete. Az elsőrendű algoritmus átmeneti függvényéből látható az algoritmus PI jellege (3.7. ábra):
3.7. ábra. Diszkrét PI algoritmus átmeneti függvénye és kapcsolata a q paraméterekkel
Ha u(1)>u(0), akkor az elsőrendű szabályozási algoritmus hasonlít folytonos idejű PI algoritmusra, ahol:
K q0
erősítési tényező
c I (q0 q1 ) / K
integrálási időállandó
A PI viselkedés csak pozitív integrálási együttható esetén teljesül. Tehát:
cI 0 Az előbbiek figyelembevételével a rendszer impulzus-átviteli függvénye (3.26): K[1 (c I 1)z 1 ] G C ( z) 1 z 1
(3.26)
A fenti algoritmusok paramétereinek megváltoztatásával újabb algoritmusokhoz juthatunk. Ha például a (3.24) egyenletben a q0= 0 helyettesítéssel élünk, akkor egy digitális integrátort kapunk, melynek impulzus-átviteli függvénye (3.27): q1 z 1 G C (z) 1 z 1
(3.27)
és ez alapján a differencia egyenlete (3.28):
u(k ) u(k 1) q1e(k 1) Ez a szabályozó egyenértékű egy digitális integrátorral.
59
(3.28)
Digitális PID szabályozási algoritmusok
Egy másik speciális esetben a (3.14) impulzus-átviteli függvényből az integráló tag elhagyásával egy PD algoritmushoz juthatunk, mely sebességalgoritmus (3.29): G C (z) q 0 q 2 z 1
(3.29)
Ezt az impulzus-átviteli függvényt differencia egyenlet alakba is írhatjuk (3.30):
u(k ) q0 e(k ) q2 e(k 1)
(3.30)
Még egyszerűbb esetben az I- és D-tag egyidejű elhagyásával a P szabályzóhoz jutunk, melynek jellemzésére a (3.31, 3.32) függvényeket használhatjuk:
G C ( z) q 0
(3.31)
u(k ) q0 e(k )
(3.32)
3.4. A diszkrét PID algoritmusok modifikációi Manapság több modifikációja ismert a diszkrét idejű PID szabályozóknak. Ezek a modifikációk folytonos PID szabályozók diszkretizálásával is előállíthatóak. Az egyik modifikáció szerint az
T T u (k ) u (k 1) K e(k ) e(k 1) 0 e(k 1) D (e(k ) 2e(k 1) e(k 2)) (3.33) TI T0 sebesség algoritmusban a derivált képzését nem a hibajelből, hanem az ellenőrző jelből képezzük, tehát a 3.33 egyenlet helyett írható az
T T u (k ) u (k 1) K e(k ) e(k 1) 0 e(k 1) D ( y(k ) 2 y(k 1) y(k 2)) TI T0 (3.34) egyenlet (3.34). Ahol e(k)=r(k)-y(k). Másik lehetséges modifikáció, amikor az e(k) referencia bemenetet csak a szabályozó integráló komponensében szerepeltetjük, így a következő sebesség algoritmushoz jutunk (3.35):
60
Digitális PID szabályozási algoritmusok
T T u (k ) u (k 1) K y(k ) y(k 1) 0 e(k 1) D ( y(k ) 2 y(k 1) y(k 2)) TI T0 (3.35) Ennek a sebességalgoritmusnak az ismeretében a következő pozíció algoritmust írhatjuk fel (3.36):
T u ( k ) K y( k ) 0 TI
k
TD
(r(i) y(i)) T i 0
0
( y(k ) y(k 1))
(3.36)
Így a rendszer struktúráját tekintve az alábbi formához jutunk (3.8. ábra):
3.8. ábra. Diszkrét PID algoritmus modifikáció hatásvázlata
A fenti PID modifikáció abban az esetben előnyös, ha a szabályozási rendszernek főleg az n(k) zavarójel hatását kell kompenzálnia, azaz a rendszert folyamatosan zavarójelek érik.
3.4.1. Egyéb P- I - D struktúrák A P, I, D tagokból összeállított diszkrét szabályozási algoritmusok fontos szerepet játszanak az ipari rendszerekben. Becslésünk szerint a különböző PID struktúrák modifikációi a visszacsatolt rendszerek több mint 90 %-ában megtalálhatók, alkalmazzák és egyszerűségük és hatékonyságuk életképessé teszi azokat. Mivel magyarázható ezen algoritmusok
61
Digitális PID szabályozási algoritmusok
hatékonysága, elterjedése? Milyen esetekben javasolható beépítésük a visszacsatolásba? Erre a válasz a következő: Abban az esetben, ha a visszacsatolandó rendszer fokszáma 2, azaz másodrendű objektum szabályozásáról van szó, a PID struktúrával minden minőségi követelményt ki tudunk elégíteni, mivel a PID struktúrával meg tudjuk valósítani másodrendű rendszer esetén az optimális állapot visszacsatolást. A szabályozó P és D tagja másodrendű rendszer esetén biztosítja a két állapotváltozó visszacsatolását, az I tag pedig a maradó szabályozási eltérést szünteti meg. Magasabb rendű objektumok esetén a PID struktúra csak közelítőleg tervezhető, mivel az állapotváltozók egynél magasabb deriváltját is be kéne tervezni a visszacsatolásba. A mintavételezési idő megfelelő kiválasztásával a gyakorlatban nagyon sok folyamat közelíthető másodrendű dinamikával, így a PID algoritmusok hatásosan alkalmazhatóak. A digitális P-ID alkalmazott struktúrák sokfélék. A 3.9. ábrán egy PI struktúrát mutatunk be, ahol az integráló tag egy pozitív visszacsatolással működik (3.37):
Gi ( z )
(1 e
T0 Ti
T0 Ti
1 e
) z 1
(3.37)
z 1
3.9. ábra. Digitális PI szabályozó modifikációja A 3.10. ábrán látható PID struktúra differencia egyenlet rendszere (3.38): ( ) ( )
(
( ) ( )
( ) )
( )
( ) ( )
( ) 62
(
) ( )
(3.38)
Digitális PID szabályozási algoritmusok
3.10. ábra. Egymásrahatásos PID struktúra A 3.11. ábrán bemutatott struktúra differencia egyenlet rendszere a fenti módszerrel szintén felírható.
3.11. ábra. Diszkrét PID modifikáció
3.5. Az I-tag telítődésének problémái (antiwindup)
A szabályozó algoritmust a beavatkozó szerv követi, amely rendszerint telítéssel, végpozíciókkal rendelkezik. Nagy hibajeleknél a szabályozó kimenete a telítéses jelleggörbe vízszintes szakaszára kerül, és ha az integrátor tovább integrál, akkor a telítéses jelleggörbe bemenő jele tovább nő és a hiba előjelének megfordulásakor még hosszú időnek kell eltelnie, amíg a bemenőjel ismét visszakerül a jelleggörbe lineáris szakaszára. Emiatt a szabályozási tranziensek ideje szükségtelenül megnő, és a szabályozási rendszer minőségi jellemzői leromlanak, nagy amplitúdójú lengések keletkezhetnek.
63
Digitális PID szabályozási algoritmusok
Kedvező megoldáshoz jutunk, ha a munkapont a lineáris és telítéses szakasz találkozásának helyére kerül, mert akkor a hiba előjelének megfordulásakor a munkapont azonnal visszatér a lineáris szakaszra. Ehhez a beavatkozó szerv telítéses jelleggörbéjét a szabályozási algoritmusba be kell építeni, képezni kell a jelleggörbe kimenő jelének és bemenő jelének különbségét, és megfelelő értékű Ts időállandó reciprokával megszorozva (felerősítve) visszavezetni az integrátor bemenetére. Az integrátor bemenetét folyamatosan korrigáljuk a telítéses szakaszon a visszavezetett jellel, amíg csak el nem tűnik a különbségi jel, tehát a két jelleggörbe szakasz találkozási pontjába nem kerül a munkapont. Az 3.12 ábrán mutatjuk be egy antiwindup-ot tartalmazó diszkrét PID algoritmus struktúráját, ahol a differenciáló tag az objektum y(k) kimenetét dolgozza fel, a P és I elemnek hibajel bemenete van.
3.12. ábra. Modifikált PID algoritmus az integráló tag antiwindup kiegészítésével
A 3.13. ábra egy másik modifikáció hatásvázlatát mutatja be, melynek átmeneti függvényét a 3.14. ábrán láthatjuk.
64
Digitális PID szabályozási algoritmusok
3.13. ábra. Egymásra hatásmentes PID struktúra
3.14. ábra. A korlátozóval ellátott szabályozó átmeneti függvénye (To=1; Kc=1; 1/Ts=0.1; Td/To=2; To/Ti=0.2) Az 3.15. ábrán egy olyan PID modifikált struktúrát mutatunk be, ahol a szabályozó D tagja sorba van kapcsolva egy egyidőállandós, egységnyi erősítésű taggal, s az így ernyedő tipusú differenciáló elemhez jutunk, mely az e(k) nagysebességű változásait csillapítja. Az algoritmus To=1s mintavételezési időre terveztük, az ernyedési időállandó, Ter = 2 s, a telítés határértékei: 0 és 5.
65
Digitális PID szabályozási algoritmusok
3.15. ábra. Ernyedési időállandót tartalmazó D-taggal rendelkező PID algoritmus és átmeneti függvénye
3.6. Digitális PID algoritmusok hangolása - beállítási szabályai A folyamatos működésű P, PI, PID szabályozók beállítási szabályai a szakirodalomból ismertek. Ezek a szabályozók átvihetők a digitális paraméteres szabályozók tervezéséhez is, feltéve, hogy a mintavételezési idő a rendszerben a szabályozott objektum időállandójához képest kicsi. A következőben felsorolt beállítási módszerek akkor fogják teljesíteni várakozásainkat, ha a To mintavételezési idő és az objektum domináns T időállandója közt az eltérés kb. egy nagyságrend (3.39): T 10 T0
(3.39)
és a szabályozott szakasz aperiodikus lengésmentes viselkedésű. A digitális PID algoritmus paramétereit a szabályozott szakasz átmeneti függvénye alapján lehet meghatározni. Az egyik módszer alapján az objektum dinamikáját a G( s)
K 0 Th s e Ts 1
66
(3.40)
Digitális PID szabályozási algoritmusok
egyidőállandós, holtidős taggal helyettesítjük (3.40). Az átviteli függvényben a K0 a közelítő egyidőállandós objektum erősítési tényezője, T az időállandója, a Th pedig a tiszta és a látszólagos holtidő összege. A Th és T idők az átmeneti függvény inflexiós pontjához húzott érintő segítségével határozhatók meg.
3.16. ábra. Az önbeálló szakasz átmeneti függvénye és a közelítő paraméterek A 3.16. ábrán egy rendszer paramétereinek meghatározását mutatjuk be, s ezek a paraméterek a fenti példában pld. a következők: becsült időállandó :
T 340 s
becsült látszólagos holtidő :
Th 70 s
becsült erősítési tényező :
K0 125
Ezen közelítő paraméterek mellett, feltételezve, hogy T/T010, megbecsülhetők egy GC ( z ) R
q0 q1 z 1 q2 z 2 1 z 1
(3.41)
alakú digitális PID algoritmus paraméterei (3.41). A módszerrel a következő P-I-D szabályozó paramétereket kapjuk (3.42, 3.43, 3.44): P szabályozó:
R 0,3
67
T K 0Th
(3.42)
Digitális PID szabályozási algoritmusok
q0 1 q1 1 q2 0 R 0,35
PI szabályozó:
T K 0Th
(3.43)
q0 1, 08 q1 1 q2 0 R 0, 6
PID szabályozó:
T K 0Th
(3.44)
Th T T q1 (1 10 h ) T T q2 5 h T q0 1,1 5
Az 3.16. ábra szerint viselkedő objektumhoz a (3.42, 3.43, 3.44) képletek alapján meghatározhatjuk a digitális P-I-D algoritmus paramétereit. Az 3.17 ábra ezen illesztési eljárás eredményét illusztrálja. Látható, hogy a módszer stabilis zártköri viselkedést eredményez, de a szabályozó hangolását javítani kell. R=0.0233 qo=2.1294 q1=-3.0588 q2=1.0294 R r(k)=10
Gain
qo.+q1.z-1+q2z-2 1-z-1
125 5 4 1.728e6s +2.261e6s +777600s3+25044s2+276s+1
Gc(z)
Go(z)
To=30 s
68
Scope
Digitális PID szabályozási algoritmusok
3.17. ábra. SIMULINK diagramm és a zárt kör viselkedése Az ilyen típusú paraméterillesztést „ökölszabály” szerinti paraméter beállításnak nevezzük, melyre az jellemző, hogy végeredményként egy nagyon közelítő illesztést eredményez, így a paramétereket tovább kell finomítani. Nagy előnye ezen módszereknek az, hogy van egy kiindulási paraméterbázis, melyről elindítható a paraméterek pontosítása. Beállítási szabályokat a szakirodalomban találhatunk a különböző PID algoritmusok modifikációihoz is. Ilyen módszert ismertetnek Y. TAKAHASHI és szerzőtársai, mely a következő sebesség algoritmus számára adja meg a beállítási paramétereket (3.45):
T T u (k ) K y (k ) y (k 1) 0 e(k ) D y (k ) 2 y (k 1) y (k 2) TI T0
(3.45)
A paraméterek számításának kiinduló adatait ugyanúgy, mint az első beállítási módszernél a szabályozandó objektum átmeneti függvénye alapján kell meghatározni. A kiindulási adatok tehát az objektum időállandója (T) látszólagos holtideje (Th) és erősítése (K0) valamint a T0 mintavételezési idő. Ezen adatok alapján az algoritmus számára meghatározhatók a K, T0/TI és TD/T0 paraméterek, melyek beállításával egy gyors túllendülő viselkedést kapunk ugrásfüggvény alapjel változás esetén.
69
Digitális PID szabályozási algoritmusok
E módszer alapján a P-I-D tagokból álló digitális szabályozási algoritmusok paramétereit a következőképpen számíthatjuk (3.46, 3.47, 3.48): P szabályozó esetén: KC
T (Th T0 ) K 0
(3.46)
PI szabályozó esetén:
K
0,135TT0 0,9T ; (Th T0 ) K 0 (Th T0 / 2)2 K 0
Tp TI
0, 27TT0 K (Th T0 / 2) 2
(3.47)
PID szabályozó esetén: K
0,3TT0 1, 2T (Th T0 ) K 0 (Th T0 / 2) ? K 0
T0 0, 6TT0 TI K (Th T0 / 2) 2
(3.48)
TD 0,5T T0 KT0
Az egyenletek csak akkor alkalmasak az algoritmus paramétereinek számítására, ha a Th/T0 hányadosa 1. Amint látható a legelső beállítási módszernél a paraméterszámító képletekben nem szerepelt a mintavételezési idő, mivel feltételezzük, hogy az egy nagyságrenddel kisebb az időállandónál; az utóbbi esetben azonban a számítási képletek a T0 mintavételezési időt is figyelembe veszik. A folytonos rendszerek P, PI és PID szabályozóinak beállítására a Zeigler-Nichols módszert is alkalmazzák, s a módszert átültették digitális P, PI és PID algoritmusok paramétereinek számítására is. A Zeigler-Nichols módszer olyan szabályozási körök tervezésénél alkalmazható, amikor a szabályozott objektum minimálfázisú. A módszer alapelve az, hogy a szabályozandó objektumot egy arányos P szabályozóval visszacsatoljuk és vagy kísérlettel, vagy számítógépes modellezéssel megkeressük a körben lévő P szabályozónak azt a kritikus erősítését (K*), amelynél a zárt kör
70
Digitális PID szabályozási algoritmusok
állandósult lengési állapotba kerül. A K* kritikus erősítés alatt a zárt rendszer lecsillapodik, fölötte pedig instabil lesz. Az 3.18. SIMULINK ábrán bemutatott rendszer esetében K*= ~0,0432 kritikus erősítést kapunk, az állandósult lengések periódusideje (T*) pedig ~ 400 s.
3.18. ábra. SIMULINK diagram és a zárt kör viselkedése
A kritikus K* erősítés és T* periódusidő segítségével megtervezhető azaz illeszthetők a q0 q1 z 1 q2 z 2 GC ( z ) R 1 z 1
(3.49)
szabályozó (3.49) R, q0, q1 q2 paraméterei a következő összefüggések alapján (3.50, 3.51, 3.52): P szabályozó esetén:
R 0,5 q0 1 q1 1 q2 0 PI szabályozó esetén: 71
K* K0 (3.50)
Digitális PID szabályozási algoritmusok
R 0, 45
K* K0
q0 1, 05
(3.51)
q1 1 q2 0 PID szabályozónál:
R 0, 6
K* K0
q0 3,5
(3.52)
q1 5,8 q2 2, 4 A mintavételezési idő ennél a tervezési (illesztési) eljárásnál kötött érték és meghatározását a T* periódusidő alapján végezhetjük el, ahol a mintavételezési idő (3.53): T0
T* 20
(3.53)
Takahashi és szerzőtársai a módosított digitális PID algoritmusok számára is kidolgozták a Zeigler-Nichols beállítási módszert, melynél a T0 mintavételezési idő szintén szabadon megválasztható néhány megkötéssel.
3.7. PID algoritmusok paraméterillesztése kritériumok alapján A digitális szabályozó algoritmusok paramétereinek tervezését a folytonos rendszerek tervezésénél használt minőségi integrálkritériumok segítségével is elvégezhetjük, azonban a kritériumot a szaggatott működés miatt diszkretizált alakban kell megadni. Ezek az összeg kritériumok a következők lehetnek (3.54):
összeg kritérium
négyzetes kritérium
72
Digitális PID szabályozási algoritmusok
abszolút érték kritérium
idővel súlyozott abszolút érték kritérium
idővel súlyozott négyzetes kritérium
(3.54)
Az I1 összegkritériumot olyan szabályozási rendszerek esetén lehet alkalmazni, melyek e(k) hibajele előjelet nem vált. Az egyik leggyakrabban használt kritérium az I2 négyzetes összegkritérium és az I3 abszolút érték kritérium. A négyzetes kritérium használata kis szabályozási eltérések esetén többször oszcilláló viselkedést eredményez, mely csökkenthető az I3 és I4 és I5 kritériumok használatával. Az I2 kritériumot legtöbbször az analitikus kezelhetősége miatt használják, matematikai előnyei miatt. A fenti összeg kritériumok használatával gyakran a szabályozott jellemző értéken tartása, illetve átvitele az egyik stacionárius pontból a másikba a gyakorlatban kivitelezhetetlen nagyságú végrehajtójelekkel történik, vagy olyan nagy beavatkozásokkal, melyet
az
irányított
rendszernél
nem
engedhetünk
meg.
Ilyen
esetekben
az
integrálkritériumokban a végrehajtó jelet is figyelembe vesszük az 3.55 egyenlet alapján:
I 6 e2 (k ) r u 2 (k ) k 0
(3.55)
ahol az r a végrehajtójel büntető faktora a u(k) a következő összefüggés alapján számítható az u stacionárius végrehajtójel segítségével (3.56): u(k ) u(k ) u
(3.56)
Az I6 összetett kritériumban az r büntető faktor nagyságára a szakirodalom ad felvilágosítást, mely szerint meghatározható az ugrásfüggvény alapjel váltásra tervezett rendszereknél az objektum K0 erősítési tényezője segítségével. 73
Digitális PID szabályozási algoritmusok
A konkrét optimalizálási feladat során mindenekelőtt ki kell választani azt a matematikai módszert, amely a végső eredményekhez egyrészt a legkevesebb számítási munka felhasználásával vezet, másrészt a felhasználó számára a legtöbb információt szolgáltatja. A szakirodalomban nagyon sok módszer terjedt el a paraméteroptimalizálás tématerületén belül, mi azonban egy szimulációval összekapcsolt eljárást ismertetünk. A szimulációval összekapcsolt paraméter-optimalizálási eljárás a tervező mérnökök számára, illetve a digitális szabályozástechnikát megismerni kívánók számára célszerű eljárásnak bizonyul. Ennél a módszernél a célunk nem a számítási idő lehető legnagyobb csökkentése, hanem az, hogy a tervező belelásson az optimalizáló eljárásokba, azt követni tudja, saját elképzeléseit tudja megvalósítani az eljárás segítségével. A szimulációval összekapcsolt optimalizáló eljárás általános célfüggvényei legyenek a következők (3.57, 3.58): I opt .1
1 M k N ( w(k ) y(k ) 2 r u 2 (k ) M 1 K 0
(3.57)
1 M k N (w(k ) y(k ) ru(k ) M 1 K 0
(3.58)
I opt .2
amelyek teljesen lefedik az I2-I6 célfüggvényeket. A MATLAB egyik tervező modulja (Nonlinear Control Design) a 3.57, 3.58 célfüggvények alapján képes megkeresni a célfüggvény minimumát optimalizációs algoritmusok alapján nemlineáris feltételek között is. Ezt illusztráljuk a következő példával (3.19. ábra). Legyen a digitális szabályozási rendszerünk struktúrája a következő:
74
Digitális PID szabályozási algoritmusok
3.19. ábra. NCD block a SIMULINK felületen ahol a nemlineáris objektum struktúrája tartalmazza a végrehajtó szerv telítési jelleggörbéjét (3.20. ábra):
3.20. ábra. Nemlineáris objektum A visszacsatolásban elhelyezett digitális PID algoritmus hatásvázlata (3.21. ábra):
3.21. ábra. Digitális PID algoritmus 75
Digitális PID szabályozási algoritmusok
A digitális PID algoritmus kezdeti paraméterei legyenek: Kc = 0.5000 Ki = To/Ti = 0.1 Kd = Td/To = 0 Az optimalizálást végző NCD algoritmus a fenti paraméterekről indulva az 3.22. ábra szerinti zártköri viselkedést mutatja (alsó görbe), majd az iterációk után a megadott 10 %-os túllendülés, 8 s-os felfutási idő, 20 s-os beállási idő és 5 %-os beállási pontosság megadásával megmutatja a tranziens viselkedést és megadja a digitális PID algoritmus paramétereit:
3.22. ábra. Az optimalizálás végeredménye Kp = 10.3795 Ki = To/Ti = 0.4473 Kd = Td/To = 30.1721 Az NCD MATLAB/SIMULINK szoftver blokk nagyon hatékony eszköz szabályozási paraméterek optimalizálásához.
76
DAHLIN algoritmus
4. DAHLIN ALGORITMUS
4.1. A klasszikus Dahlin algoritmus elve A számítógépes folyamatirányítás területén az ipari rendszerekben olyan folyamatokkal találkozunk, melynek tranziens viselkedését a zárt szabályozási körökben egyszerű átmenetekre kell megterveznünk. Az ilyen egyszerű viselkedésű zárt köröknek a folytonos átviteli függvénye holtidőt feltételezve legyen (4.1):
y ( s) e sTh Gr ( s ) w( s) TE s 1
(4.1)
ahol: TE – a zárt szabályozási kör időállandója TH - holtidő Tehát szeretnénk megvalósítani egy olyan szabályozási kört, amelynek a viselkedése egy elsőrendű holtidős, vagy holtidő nélküli rendszert képvisel. A Dahlin [6] algoritmus levezetéséhez induljunk ki a zárt digitális szabályozási kör összefüggéseiből. A zárt szabályozási kör átviteli függvénye (4.2):
Gr ( z )
GC ( z )G0 ( z ) 1 GC ( z )G0 ( z )
ahol: G0(z) a szabályozott objektum diszkrét modellje a tartószervvel együtt. A zárt rendszer felépítése (4.1. ábra):
4.1. ábra. A zárt szabályozási kör 77
(4.2)
DAHLIN algoritmus
Fejezzük ki a (4.2) egyenletből a GC(z) impulzus-átviteli függvényt. Így a digitális szabályozó impulzus-átviteli függvényét kapjuk (4.3):
GC ( z )
Gr ( z ) 1 G0 ( z ) 1 G r ( z )
(4.3)
A (4.1.) egyenlet folytonos rendszert ír le, melynek meghatározhatjuk a Z- átviteli függvényét (4.4): GW ( z )
T 1 (1 e T0 ) z 1 d z ;d h ; T0 1 T0 TE 1 e z
(4.4)
Ahhoz, hogy meg tudjuk határozni a GC(z) digitális szabályozó algoritmust, szükségünk van a G0(s) objektum Z- átviteli függvényére, amely ha elsőrendű a következő lesz (4.5):
G0 ( z )
b1 z 1 d z 1 a1 z 1
(4.5)
ahol az együtthatók értékei (4.6):
b1 K 0 (1 e a1 e
d
T0 T1
);
(K0 - az objektum erősítési tényezője és T1 - időállandója)
T0 T1
Th T0
(4.6)
Ha a (4.4) és a (4.5) kifejezést behelyettesítjük a (4.3) szabályozó egyenletébe, akkor a
(1 a1 z 1 )(1 e T0 ) Gc ( z ) b1 (1 e T0 z 1 (1 e T0 ) z d 1 )
(4.7)
digitális szabályozó átviteli függvényét kapjuk. A (4.7) szabályozási algoritmusban a szabályozónak csak egy hangolási paramétere van, ez pedig a λ = 1/TE. A gyakorlatban nagyon gyakran másodrendű holtidős rendszerekkel találkozunk, melyek folytonos átviteli függvénye (4.8):
K 0 e sTh G0 ( s ) (T1 s 1)(T2 s 1) 78
(4.8)
DAHLIN algoritmus
impulzus-átviteli függvénye (4.9):
G0 ( z )
K 0 (b1 z 1 b2 z 2 ) d z 1 a1 z 1 a2 z 2
(4.9)
ahol a folytonos rendszerből származtatható paraméterek a következők (4.10): b1 1
1 ( L2 e L1T0 L1e L2T0 ) L1 L2
b2 e ( L1 L2 )T0
1 ( L2 e L2T0 L1e L1T0 ) L1 L2
a1 e
a2 e
L1
T0 L1
T§ L1
e
T0 L2
e T0 L2
1 1 ; L2 T1 T2
(4.10)
Az elsőrendű esethez hasonlóan egyszerű behelyettesítéssel a (4.8) rendszerhez tartozó DAHLIN algoritmust megtervezhetjük, melynek impulzus-átviteli függvénye (4.11): ( )
( (
)(
) (
)(
)
)
(4.11)
melynek egyetlen hangolási paramétere a .
4.2. Példák A következő példákon keresztül bemutatjuk az algoritmus tulajdonságait tipikus rendszerekhez. 1. példa. Elsőrendű objektum Legyen az objektum folytonos átviteli függvénye: G 0 (s)
y(s) 0,6 , u (s) 67s 1
és szeretnénk, ha a zárt kör viselkedése a következő legyen:
79
DAHLIN algoritmus
G r (s)
y(s) 1 . r (s) 20s 1
Tervezzünk Dahlin algoritmust, mely biztosítja a TE=20 s időállandóval rendelkező elsőrendű átmenetet alapjelváltásra T0=2 s mintavételezési idővel. A folytonos G0(z) rendszer impulzus-átviteli függvénye:
y( z ) 0.01823z 1 G 0 (z) u (z) 1 0,4706z 1 a zárt kör impulzus-átviteli függvénye:
G r (z)
0,09516z 1 . 1 0,9048z 1
A (4.3) egyenlet alapján a Gc(z) szabályozó algoritmus átviteli függvénye:
G c (z)
u (z) 0,09516 0,1785z 1 0,08357z 2 e(z) 0,001823 0,03473z 1 0,0165z 2
az átmeneti függvénye pedig a 4.2. ábrán látható.
4.2. ábra. A megtervezett Dahlin algoritmus átmeneti függvénye Beillesztve a Gc(z) algoritmust a szabályozási körbe a rendszer viselkedését láthatjuk a 4.3. ábrán.
80
DAHLIN algoritmus
4.3. ábra. A zárt kör viselkedése Dahlin algoritmussal A Dahlin algoritmus pólusai és zérusait ábrázoljuk az 4.4. ábrán, ahol a 0,9048 póluszérus helyek egymást kompenzálják. (Megjegyzés: a MATLAB az egyszerűsítéseket nem végzi el!)
4.4. ábra. A Dahlin algoritmus pólusai és zérusai Pólusok: Zérusok:
1.0000 0.9048 0.9706 0.9048
81
DAHLIN algoritmus
A megtervezett algoritmust megvizsgálva láthatjuk, hogy ez egy PI algoritmus viselkedést mutat, hisz egy pólus és zérus azonos, így kiejtik egymást.
2. példa: Másodrendű objektum Legyen a szabályozott objektum egy másodrendű lengő rendszer, melynek átviteli függvénye: G 0 (z)
0,5 66s 12s 1 2
a mintavételezési idő: T0=4. A folytonos zárt rendszer időállandója legyen TE=12 s. A kiszámított Dahlin algoritmus:
G c (z)
5,9903 7,8719z 1 2,8945z 2 , 1 0,2161z 1 0,7839z 2
melynek átmeneti függvénye (4.5. ábra):
4.5. ábra. Másodrendű objektumhoz tervezett Dahlin algoritmus átmeneti függvénye Pólusait és zérusait pedig a 4.6. ábrán láthatjuk:
82
DAHLIN algoritmus
4.6. ábra. Másodrendű objektumhoz tervezett Dahlin algoritmus pólusai és zérusai Pólusok: Zérusok:
1.0000 0.7165
0.7165 0.6572+.2265i
-0.7839 0.6572-.2265i
Valamint az időtartománybeli viselkedést tanulmányozhatjuk az 4.7. ábrán:
4.7. ábra. A zárt kör viselkedése Dahlin algoritmussal
83
DAHLIN algoritmus
Megfigyelhetjük az ún. ringing jelenséget, melyet a -0,7839 helyen lévő negatív valós pólus okoz. A RINGING jelenséget ezen negatív valós pólus eltávolításával meg lehet szüntetni, így elkülöníthetők a stacionárius végrehajtójel körüli lengések, illetve a mintavételezési idők közötti kimenőjel hullámzás, lengés, amit az 4.8. ábra szemléltet.
4.8. ábra. A RINGING jelenség hatása a zárt körben A RINGING jelenséget megszüntető eljárás természetesen a zárt kör viselkedését megváltoztatja, de az u(k) végrehajtójel lengése és a mintavételezések közti y(t) lengések megszűnnek. A példában másodrendű rendszer szerepel, így a negatív valós gyök megszüntetésével egy PID algoritmust kapunk. Felhívjuk a figyelmet arra, hogy ilyen esetekben a módszer PID algoritmusok paraméter illesztésére is alkalmas. Írjuk fel a Dahlin algoritmusunk átviteli függvényének nevezőjét gyöktényezős alakba:
5,9903 7,8719z 1 2,8945z 2 G c (z) (1 z 1 )(1 0,7835z 1 ) és vezessük be a következő határértéket (z tart az 1-hez):
5,9903 7,8719z 1 2,8945z 2 G c (z) (1 z 1 ) lim(1 0,7835z 1 ) A határérték kiszámításával a Gc(z) algoritmus a következő PID tulajdonságú algoritmussá alakul: G
PID c
3,3578 4,4122z 1 1,6225z 2 . ( z) 1 z 1
Ezt az algoritmust helyettesítve a Dahlin algoritmus helyére a mintavételezések közti hullámzás/lengés megszűnik, de a tervezett zártköri viselkedés is módosul, amint a 4.9. ábra 84
DAHLIN algoritmus
bizonyítja:
4.9. ábra. A zárt kör viselkedése Dahlin algoritmusból származtatott PID algoritmussal A két algoritmus átmeneti függvényét összehasonlítva láthatjuk, hogy az integrálási tulajdonságaik megegyeznek.
4.10. ábra. A Dahlin és a származtatott PID algoritmus átmeneti függvénye 3. példa. Integráló viselkedésű objektum Legyen a szabályozandó objektum integráló jelleggel bíró rendszer, melynek folytonos átviteli függvénye: G 0 (s)
1 . s(33s 1)
Tervezzünk Dahlin algoritmust, ahol a zárt kör legyen szintén egyidőállandós, melynek 85
DAHLIN algoritmus
időállandója TE=10 s. A Dahlin algoritmus átmeneti függvénye látható a 4.11. ábrán:
4.11. ábra. Dahlin algoritmus átmeneti függvénye integráló jellegű objektum esetén mely alapján látható, hogy egy PD jellegű algoritmust kapunk, ugyanis az objektum rendelkezik integráló jelleggel. Ha megnézzük a szabályozó pólus-zérus elhelyezkedését láthatjuk, hogy az algoritmus egy negatív valós pólussal és egy pozitív zérushellyel rendelkezik. Pólusok: Zérusok:
1.0000 1.0000
0.6703 0.6703
-0.6204 0.2335
Ebből felírható a szabályozó átviteli függvénye.
4.3. MATLAB program analitikus tervezéshez A következő részben egy MATLAB m-file demonstrálja a tervezés egyszerűségét. Tételezzük fel, hogy To=4 s mintavételezéssel identifikáltuk az irányítandó objektumot, mely modelljét a programban Go-val jelölünk. Látható, hogy az identifikáló algoritmus d=2 diszkrét holtidőt azonosított, ami 8 s tiszta holtidőnek felel meg. A zárt köri viselkedést leíró Gr átviteli függvény szintén d=2 diszkrét holtidőt kell, hogy tartalmazzon. Az identifikált objektum B, A polinomjait, illetve a megtervezett Gc algoritmus Q,P polinomjait használja fel SIMULINK modellünk az időtartománybeli vizsgálathoz. 86
DAHLIN algoritmus
% DAHLIN ALGORITMUS TERVEZÉSE % (Dr. Vass József, Pannon Egyetem) % To- mintavételezési idő To=4 % Az identifikált objektum impulzus-átviteli függvénye (d=2) To- mintavételezési idő Go=tf([0.08232 0.0258],[1 -1.064 0.2881 0 0],To) % Az objektum B és A polinomjai [B,A]=tfdata(Go,'v') % A kívánt zártköri veselkedés folytonos átviteli függvénye holtidő nélkül gr=tf(1,[20 1]) % A kívánt zártköri veselkedés diszkrét átviteli függvénye d=1 diszkrét holtidővel Gr=tf(0.1813,[1 -0.8187 0 0],4) % A DAHLIN algoritmus Gc=1/Go*(Gr/(1-Gr)) % Az algoritmus átmeneti függvénye step(Gc) grid % Az algoritmus Q és P polinomjai [Q,P]=tfdata(Gc,'v') pause % Az algoritmus zérus-pólus elrendezése grid zero(Gc),pole(Gc) zpk(Gc)
4.12. ábra. A tervezőprogram szolgáltatásai Az identifikált Go objektum B,A polinomjait, illetve a megtervezett Gc algoritmus Q,P polinomjait használja fel SIMULINK modellünk az időtartománybeli vizsgálathoz.
87
DAHLIN algoritmus
4.13. ábra. Az algoritmus tesztelése SIMULINK segítségével Az algoritmus tesztelésekor az 4.14. ábra alapján azt mondhatjuk, hogy tervezésünk sikeres volt.
4.14. ábra. A zárt rendszer viselkedése holtidős rendszer esetében A zártköri tranziensek - u(k), y(k) - láthatók a 4.15. ábrán.
88
DAHLIN algoritmus
4.15. ábra. A végrehajtó és ellenőrző zártköri jelek
4.4. Esettanulmány Bevezetés
A Pannon Egyetemen üzembe helyeztünk egy DCMCT (DC Motor Control Trainer) rendszert. A rendszer folyamatirányítási felülete PC bázisú, mely USB felületen kapcsolódik egy PIC18F4550 mikrokontroller alapú egységhez, az pedig egy MAXON DC motorhoz (Graphite Brushless, 18 Watt, 26 mm). A DCMCT dinamikus viselkedését egy rúgó-tömeg mechanikus egység segítségével megváltoztattuk és a beépített PID algoritmus helyett diszkrét Dahlin algoritmust használunk. Az analitikusan megtervezett Dahlin algoritmus tartalmaz egy notch szűrőt/szabályozót, ill. egy PI algoritmust.
89
DAHLIN algoritmus
4.4.1. A kísérleti rendszer felépítése Az QUANSER rendszer [1] rendszertechnikai felépítése az 4.16. ábrán látható.
4.16. ábra. A QANSER DC motor Kit felépítése A rendszer USB vonalon keresztül kapcsolódik a számítógéphez, melyen fut egy realtime folyamatkezelő szoftver. A mintavételezési idő To= 0,01 s. A gyakorló eszköz a rugótömeg kiegészítéssel a 4.17.-4.18. ábrákon látható.
4.17. ábra. A QUANSER DC motor rugó-tömeg kiegészítő mechanikával
90
DAHLIN algoritmus
4.18. ábra. A rugó-tömeg ráépítés A QUANSER USB QICii szoftver kezelői felülete (4.19. ábra) pozíció és fordulatszám szabályozást enged meg, valamint rendelkezik egy egyszerű modellezési lehetőséggel
is.
Ez utóbbit
használtuk fel
a mechanikailag átalakított
rendszer
identifikációjánál az adatgyűjtéshez. Az 4.19. ábrán látható, hogy a rendszer keskenysávú lengéseket végez az átmeneti állapotokban. E viselkedés notch-szűrő alkalmazását indokolhatja.
4.19. ábra. Az adatgyűjtő QICii szoftver kezelői felülete 91
DAHLIN algoritmus
4.4.2. Az identifikáció Szabályozási rendszerek tervezésének megvalósításához szükségünk van a szabályozott objektum dinamikai tulajdonságainak ismeretére. Valamilyen leírási mód segítségével rendelkeznünk kell az objektum matematikai modelljével, hisz a tervezőrendszerünk kérni fogja a szabályozandó objektum matematikai modelljét. Szerencsés esetben az objektum matematikai modellje – fizikai, kémiai, biológiai paraméterekkel ismert, de a legtöbb esetben nem, így valamilyen közelítő munkaponti modellekkel kell a tervezést megkezdeni. A technológiai rendszerek irányítását, szabályozását, vezérlését számítógépek végzik. A számítógép számára tervezendő szabályozási algoritmusokat a folyamatirányító számítógép differencia egyenlet formájában kéri, valamilyen tervezett T0 mintavételezési idő mellett, azaz az algoritmus egy diszkrét algoritmus, így a számítógép is egy diszkrét rendszert „szeretne látni” szabályozott objektumként, amely a valóságban természetesen időben folytonos rendszer. Az időben folytonos rendszerek identifikációjával a szakirodalom több évtizede foglalkozik [2] és letisztultak azok az identifikációs eljárások, amelyeket a gyakorlat is közvetlenül alkalmazhat. Kialakultak azok a diszkrét rendszermodellek, melyek a gyakorlati problémák legtöbbjének leírására megfelelnek, melyek felhasználásával olyan matematikai modelleket kapunk, amelyek közvetlenül felhasználhatók digitális szabályozási algoritmusok tervezéséhez, adaptív rendszerek megalkotásához [3]. Az QICii rendszerből származtatott mérési adatok formátuma olyan, hogy a MATLAB System Identification Toolbox-a ezeket értelmezni tudja. Az 4.20. ábra mutatja a feldolgozandó be-kimeneti jelsorozatot.
92
DAHLIN algoritmus
4.20. ábra. Az identifikáció adatsora (To=0.01 s) A rendszerhez illesztett OE (Output Error) modell, mely mérési zajt tartalmazó modell - a következő hatásvázlattal jellemezhető (4.21. ábra):
4.21. ábra. Az OE modell struktúrája A modell impulzus-átviteli függvénye (4.12):
y(z)=
B(z -1 ) -d z u(z)+e(z), A(z -1 )
ahol e(z) a mérési zaj z-transzformáltja.
93
(4.12)
DAHLIN algoritmus
A zárt szabályozási kör mintavételezési idejének megválasztása A szabályozandó objektum átmeneti függvényének viselkedése (4.22. ábra) alapján a szabályozó kör mintavételezési idejét a lengési periódusidő tized részére állítjuk be, azaz legyen To= 0.01 s. Ezzel a mintavételezési idővel fogjuk üzemeltetni a megtervezendő szabályozási kört.
4.22. ábra. Az objektum átmeneti függvénye
Identifikáció MATLAB/System IdentificationToolbox segítségével A kísérleti rendszer adatsorozatát egy z=[y u] vektorban tároltuk el, majd felhasználva a Toolbox parancsát megkapjuk az identifikált rendszer OE modelljének B(z), A(z) polinomjait. Minkét polinom fokszáma három, a diszkrét holtidő pedig egy (d=1). th=oe(z,[3 3 2]) Discrete-time IDPOLY model: y(t) = [B(q)/A(q)]u(t) + e(t) B(q) = 2.009 q-2 - 3.575 q-3 + 1.897 q-4 A(q) = 1 - 2.603 q-1 + 2.489 q-2 - 0.8666 q-3 A compare(z,th) paranccsal az identifikáció minőségére kapunk választ (4.23.ábra):
94
DAHLIN algoritmus
4.23. ábra. A mérési adatok és az identifikált modell összehasonlítása időtartományban
4.4.3. Diszkrét algoritmus analitikus tervezése A számítógépes folyamatirányítás területén, az ipari rendszerekben olyan folyamatokkal találkozunk, melynek tranziens viselkedését a zárt szabályozási körökben egyszerű átmenetekre kell megterveznünk. Az ilyen egyszerű viselkedésű zárt köröknek a folytonos átviteli függvénye - holtidőt feltételezve (4.13):
G r (s)
y(s) e sTh w (s) TE s 1
(4.13)
ahol: TE – a zárt szabályozási kör időállandója, s TH – holtidő, s Szeretnénk megvalósítani egy olyan szabályozási kört, amelynek a viselkedése egy elsőrendű holtidős (vagy holtidő nélküli) rendszert képvisel. A Dahlin [4,5] algoritmus levezetéséhez induljunk ki a zárt digitális szabályozási kör összefüggéseiből. A zárt szabályozási kör átviteli függvénye (4.14):
G r (z)
G C ( z )G 0 ( z ) 1 G C ( z )G 0 ( z )
ahol: G0(z)- a szabályozott objektum diszkrét modellje a tartószervvel együtt
95
(4.14)
DAHLIN algoritmus
Gc(z) - a tervezendő algoritmus diszkrét átviteli függvénye Gr(z) - a zárt szabályozási kör viselkedése. A zárt rendszer szabályozási kör hatásvázlata látható a 4.24. ábrán.
4.24. ábra. A zárt szabályozási kör hatásvázlata Fejezzük ki a GC(z) impulzus-átviteli függvényt, és megkapjuk a digitális szabályozó algoritmus impulzus-átviteli függvényét (4.15):
G C (z)
G r (z) 1 G 0 ( z) 1 G r (z)
(4.15)
A zártköri viselkedés tervezése A zártköri viselkedés megadásakor figyelembe kell venni a diszkrét holtidőt, mely esetünkben egy mintavételezési idő, így a zártköri folytonos modell 0,3 s-os időállandóval: Grf(s) =
(
)
A c2d(Grf,0.01) parancs segítségével To=0.01s mintavételezési idő mellett megkapjuk a Grd(z) zártköri viselkedés impulzus-átviteli függvényét. Az ehhez tartozó átmeneti függvényt ábrázolja a 4.25. ábra. Grd(z) =
96
DAHLIN algoritmus
4.25. ábra. A zárt kör tervezett viselkedése
A szabályozási algoritmus Az összefüggést felhasználva a Gc=1/Go * (Grd/(1-Grd) MATLAB paranccsal megkapjuk a diszkrét algoritmus impulzus-átviteli függvényét: Gc(z) = és a 4.26.ábra mutatja a diszkrét algoritmus átmeneti függvényét.
4.26. ábra. A diszkrét DAHLIN szabályozó algoritmus átmeneti függvénye Abban az esetben ha a zpk(Gc) paranccsal megvizsgáljuk a megtervezett algoritmus zérus és pólus helyeit: 97
DAHLIN algoritmus
Gc(z) =
( (
)( )(
)( )(
) )
azt tapasztaljuk, hogy a megtervezett algoritmus tartalmaz egy diszkrét PI algoritmust és egy sorbakötött Notch-szűrőt (4.29. ábra), mely a keskenysávú lengéseket fogja kikompenzálni [5]
A nyitott köri Bode- (4.27. ábra) és a zérus-pólus (4.28. ábra) diagram is ezt reprezentálja.
4.27. ábra. Az objektum (Go), a Dahlin algoritmus (Gc) és a nyitott kör (Gc*Go) Bode diagrammjai
4.28. ábra. Az objektum és a diszkrét Dahlin algoritmus pólus-zérus helyei
98
DAHLIN algoritmus
4.4.4. A zártköri viselkedés vizsgálata A zárt szabályozási kör SIMULINK diagramjában a Dahlin algoritmust felbontottuk egy diszkrét PI és egy diszkrét notch-szűrő/szabályozó részre (4.29. ábra).
4.29. ábra. A zárt kör SIMULINK diagramja A 4.30. ábrán az y(k) ellenőrző jel és az u(k) végrehajtó jel viselkedését láthatjuk. A zárt kör lengésmentes, a zárt kör időállandója 0,3 s. A sztochasztikus zavarásnak kitett rendszer viselkedését a 4.31. ábra, változó alapjel esetében pedig a 4.32. ábra mutatja be a rendszer viselkedéseit.
4.30. ábra. Zártköri viselkedés ugrásfüggvény alapjel váltás esetén
99
DAHLIN algoritmus
4.31. ábra. Zártköri viselkedés sztochasztikus zavarójel esetében
4.32. ábra. A zárt rendszer viselkedése változó alapjel esetén Összefoglalás A Pannon Egyetemen üzembe helyezett és mechanikailag átalakított QUANSE DC motorhoz QUANSER USB QICii szoftver és a MATLAB/SIMULINK/IdentificationToolbox/Control System Toolbox bázisán diszkrét Dahlin algoritmust terveztünk. A megtervezett algoritmus tartalmaz egy diszkrét PI és egy notch-szűrő/szabályozó komponenst. Az esettanulmány bemutatja egy tervezés minden lépését, így felhasználható az irányításelmélet gyakorlati oktatásában, tervezési segédeszköze lehet a gyakorló mérnöki munkának.
100
DAHLIN algoritmus
Irodalomjegyzék [1] K. J. Aström, J. Apkarian: DC Motor Trainer (DCMCT), Instructor Workbook, QUANSEREngineering, Canada, 2008. [2] L. Ljung: System Identification, TheoryfortheUser, Prentice Hall, 1987. [3] R. Isermann: Digital Control Systems, Springer-Verlag, 1989. [4] K. Dutton,S.Thompson:The Art of ControlEngineering,Addison-Wesley, 1997. [5] B. C. Kuo: Digital Control Systems, Sounders College Publishing, 1992. [6] E. B. Dahlin: Designing and Tuning Digital Control Systems. Instr. Cont.Syst.,42, (1968)
101
DEADBEAT algoritmus
5. VÉGES BEÁLLÁSI IDEJŰ (DEADBEAT) SZABÁLYOZÁSI ALGORITMUSOK A véges beállási idejű (DEADBEAT) szabályozási algoritmusok diszkrét algoritmusok. Alaptulajdonságuk az, hogy előírt változású alapjel vagy zavaró jellemző esetén a szabályozott jellemző véges számú mintavételezési idő után eléri az állandósult állapotot, azaz egy m-ed rendű objektum esetén az m-edik mintavételezési időnél felveszi az r(k) parancsolt értéket. Előírt alapjel váltásra tervezik ezeket a szabályozókat (ugrásfüggvény, sebességugrás), amikor az alsó hierarchia szinten lévő szabályozó körök legfőbb feladata az, hogy a felsőszint utasításait minimális idő alatt teljesítsék. Legyen egy zárt szabályozási kör (5.1. ábra), melynek impulzus-átviteli függvénye (5.1) Gr ( z )
GC ( z )G0 ( z ) y( z ) w( z ) 1 GC ( z )G0 ( z )
(5.1)
melyből a szabályozó impulzus-átviteli függvénye (5.2): GC ( z )
Gr ( z ) u( z) 1 e( z ) Go ( z ) 1 Gr ( z )
(5.2)
5.1. ábra. A digitális szabályozókör hatásvázlata A rendszer bemenete az r(k) alapjel, melyek impulzus-átviteli függvénye (5.3): r ( z)
1 1 z 1
a zárt kör kimenő y(k) jele pedig viselkedjen az 5.2. ábra szerint.
102
(5.3)
DEADBEAT algoritmus
5.2. ábra. A véges beállási idejű szabályozóval visszacsatolt rendszer viselkedése
Az ábra alapján az egyes mintavételezési időhöz rendelt kimenőjelek transzformáltja felírható (5.4): y( z ) y(1) z 1 y(2) z 2 ... 1 z m z ( m1) ...
(5.4)
Ezt az y(k) kimeneti jelsorozatot, csak olyan u(k) jelsorozat hozhatja létre, mely m mintavételezési ciklusidő alatt szintén beáll a stacionárius végrehajtójel értékre (5.5): u( z ) u(o) u(1) z 1 ... u (m) z m z ( m1) ...
(5.5)
Ilyen kikötések alapján vizsgáljuk meg az y(k) kimeneti és u(k) szabályozó kimeneti jelsorozat z-transzformáltjának és a r(k) alapjel z-transzformáltjának hányadosait (5.6). y( z ) 1 z 1 y(1) z 1 y(2) z 2 ... z m y(1) z 1 y(2) y(1) z 2 ... 1 y(m 1) z m r ( z) (5.6)
ahol az egyes együtthatókat pi-vel jelölve a következőket kapjuk (5.7): y( z) P( z ) p1 z 1 p2 z 2 ... pm z m r ( z) Ugyanezt az eljárást elvégezve a szabályozó kimenőjelével is a
103
(5.7)
DEADBEAT algoritmus
u( z) Q( z ) q0 q1 z 1 ... qm z m r ( z) polinomokat kapjuk (5.8), ahol qm=u(m)-u(m-1).
(5.8)
A 5.7 és 5.8 polinomok együtthatóira a következő összefüggések igazak (5.9, 5.10): m
pi 1
(5.9)
i 0
valamint m
ai
m
qi u (m) 1/ K 0 i m0 bi
i 0
(5.10)
i 0
ahol a K0 a szabályozandó objektum erősítési tényezője (természetes, hogy az u(m) stacionárius végrehajtójel nagysága az objektum erősítésének reciprokával egyenlő egységugrás alapjel váltás esetén). Összevetve az 5.1 és 5.7 egyenleteket a zárt kör viselkedését a P(z) polinom határozza meg (5.11): Gr(z)=P(z)
(5.11)
illetve a szabályozott objektum impulzus-átviteli függvényét kifejezhetjük az 5.7 és 5.8 egyenletek alapján (5.12) G0 ( z )
y ( z ) u ( z ) y ( z ) P( z ) / r ( z ) r ( z ) u ( z ) Q( z )
(5.12)
Felhasználva a szabályozóra felírt 5.2 impulzusfüggvényt, melybe behelyettesítjük az 5.12 függvényt, a keresett véges beállási idejű (DEADBEAT) szabályozó algoritmust kapjuk (5.13): GC ( z )
Q( z ) P( z ) Q( z ) P( z ) 1 P( z ) 1 P( z )
(5.13)
q0 q1 z 1 ... qm z m 1 p1 z 1 ... pm z m
(5.14)
amelyet a polinomokat kifejtve a
GC ( z )
104
DEADBEAT algoritmus
impulzus-átviteli függvénnyel jellemezhetjük (5.14). Az 5.14. szabályozó paramétereit a szabályozott objektum impulzus-átviteli függvényének paramétereiből egyszerű úton származtathatjuk, mégpedig az egyes qi, pi együtthatók számíthatóak. Az 5.14. egyenletből kiindulva a véges beállási idejű szabályozók általános átviteli függvénye holtidő nélküli rendszerek esetén (5.15):
GC ( z )
qo A( z 1 ) 1 qo B( z 1 )
(5.15)
holtidővel rendelkező rendszerek esetén (d=Th/To) (5.16):
GC ( z )
qo A( z 1 ) 1 qo B( z 1 ) z d
(5.16)
ahol: q0
1 1 m b1 b2 ... bm b i
(5.17)
i 1
Ugrásfüggvény alapjel váltás esetén, mivel a szabályozási eltérés a k=0 időpillanatban a legnagyobb, a szabályozó algoritmus kimenete maximális lesz, melyet az 1/bi paraméter határoz meg. Ez az u(o) beavatkozójel több gyakorlati esetben a rendszer fizikai határait túllépi, nem realizálható. Van lehetőségünk arra, hogy e legnagyobb első beavatkozást megadjuk, előírjuk. Ebben az esetben a véges beállási idő egy mintavételezési idővel eltolódik tehát m+1 mintavételezés lesz, s a qiés pi paraméterek a következőképpen számíthatóak: -a szabályozó számlálójának együtthatói (5.18): qo=u(o) – általunk megszabott első beavatkozójel m
q1 q0 (a1 1) 1/ bi i 1
m
q2 q0 (a2 a1 ) a1 / bi i 1
(5.18)
m
qm qo (am am 1 ) am 1 / bi i 1
m
qm 1 am (q0 1) / bi i 1
105
DEADBEAT algoritmus
-a nevező polinomjának együtthatói (5.19): p1 q0bi
(5.19) m
pm q0 (bm bm 1 ) bm 1 / bi i 1
m
pm 1 bm (q0 1/ bi ) i 1
és a DEADBEAT szabályozó algoritmus fokszáma eggyel megnő (5.20):
Gc ( z )
q q z 1 ... qm1 z ( m1) Q( z ) 0 1 1 1 P( z ) 1 p1 z ... pm1 z ( m1)
(5.20)
Az általunk előre megadott u(0) első beavatkozójel nagyságát nagyon kicsire nem választhatjuk, mert akkor az u(1) érték ezt meghaladja. A szakirodalom szerint jó eredményeket szolgálhat a szabályozó kör üzemeltetésekor, ha u(1)u(o). A véges beállási idejű (DEADBEAT) algoritmusok alkalmazásánál a mintavételezési idő megválasztása is fontos kérdés, mivel az u(0) első beavatkozójel nagysága a mintavételezési idő csökkenésével növekedni fog, ellenkező esetben pedig csökken ugyanazon objektum esetén. Irodalmi adatok alapján [1] a mintavételezési időt:
T0 0,18T95
(5.21)
értékben kell megválasztani (5.21), ahol a T95 a szabályozott objektum 95 %-os beállási ideje. A fenti eset az első beavatkozás korlátozatlan esetére vonatkozik. Előírt u(0) értékre a mintavételezési idő a következőkből becsülhető (5.22):
T0 0,11T95
Példa Legyen az irányítandó objektum folytonos átviteli függvénye: Go ( s)
y( s) 90s 10 4 u ( s) 117304s 81848s 3 6140s 2 143s 1
106
(5.22)
DEADBEAT algoritmus
A mintavételezési idő kiválasztásánál használjuk fel az (5.21) összefüggést és az 5.3. ábra szerint kapott adatok alapján T0=50 s-os mintavételezési idő mellett:
5.3. ábra. A mintavételezési idő meghatározása az (5.21) ajánlás alapján az objektum impulzus-átviteli függvénye:
G0 ( z )
B( z 1 ) 1, 44 z 1 1, 472 z 2 0,04527 z 3 0,000051z 4 . A( z 1 ) 1 0,9353z 1 0, 2482 z 2 0,01674 z 3 0,00001z 4
Felhasználva az (5.15.) összefüggést a minimális beállási idejű algoritmus átviteli függvénye (5.4. ábra):
5.4. ábra. DEADBEAT algoritmus SIMULINK felületen A megtervezett algoritmus átmeneti függvényét láthatjuk az 5.5. ábrán. 107
DEADBEAT algoritmus
5.5. ábra. A DEADBEAT algoritmus átmenetei függvénye A zárt rendszer viselkedését SIMULINK felület segítségével végeztük el. Megvizsgáltuk a rendszer viselkedését egységugrás alapjel változásra és zavarójel hatására. A zárt köri viselkedést bemutató ábra is bizonyítja (5.6. ábra), hogy az ún. deadbeat viselkedés nem teljesül zavarójel kompenzáció esetén, ha az nem ugrásfüggvény vagy annak integráljai (5.7. ábra).
5.6. ábra. Rendszermodell SIMULINK felületen
108
DEADBEAT algoritmus
5.7. ábra. DEADBEAT algoritmus viselkedése alapjel váltásra és zavarójelre
5.1. MATLAB program és mintapélda DEADBEAT algoritmus tervezéséhez % VÉGES BEÁLLÁSI IDEJŰ ALGORITMUS (DEADBEAT) tervezése % (Dr. Vass József, Pannon Egyetem) % A folytonos rendszer modellje Gof=tf(1,[37.5 12.5 1]) % R.Isermann [1] step(Gof) grid pause % T95 beállási idő [s] meghatározása T95=45 % A To mintavételezési idő meghatározása (javaslat) % A T95 a beállási idő 98%-a To=T95*0.18 % Diszkrét modell generálása God=c2d(Gof,To) 109
DEADBEAT algoritmus
% A B(z), A(z) polinomok [B,A]=tfdata(God,'v') % A qo számítása: qo=1/(B(1)+B(2)+B(3)) qo=1/(B(1)+B(2)+B(3)) % A DEADBEAT algoritmus, Gc=R(z)/S(z) R=qo*A S=[1 0 0]-qo*B Gc=tf(R,S,To) step(Gc) A program a következő szolgáltatásokat és eredményeket biztosítja: A folytonos átmeneti függvény alapján a T95 beállási idő ~40 s. Felhasználva a T0 0,18T95 összefüggést a mintavételezési időt válasszuk To=8 s-ra (5.8. ábra). y(t)
5.8. ábra. A mintavételezési idő meghatározása az objektum átmeneti függvénye alapján A program által megtervezett DEADBEAT algoritmus átviteli függvénye: ( ) átmeneti függvénye (5.9. ábra):
110
DEADBEAT algoritmus
5.9. ábra. A DEADBEAT algoritmus átmeneti függvénye A tervezés eredményének elemzését SIMULINK környezetben végezzük. Az a) esetben a zavarójel egységugrás, b) esetben pedig egy elsőrendű dinamikával rendelkező jel (5.10. ábra).
5.10. ábra. Modellezés SIMULINK felületen
111
DEADBEAT algoritmus
5.11. ábra. A zárt kör viselkedése egységugrás bemenetekre Amint látható az 5.11. ábrán úgy az alapjel váltásra, mint a zavarásra a DEADBEAT jelleg megmarad, az 5.12. ábra pedig azt szemlélteti, milyen viselkedést mutat a rendszer nem egységugrás zavarójel esetén.
5.12. ábra. A zárt kör viselkedése nem egységugrás zavaró jel esetén Amennyiben az irányítandó objektum holtidővel rendelkezik (mi esetünkben T h= 8 s, tehát a diszkrét holtidő d=1 a programban a B és A polinomokat a következőképpen egészítjük ki: B=[0 0 0.3713 0.1521] A=[1 -0.5461 0.0695 0] valamint a szabályozó R,S polinomjai is megváltoznak: qo=1/(B(1)+B(2)+B(3)+B(4)) 112
DEADBEAT algoritmus
R=qo*A S=[1 0 00 ]-qo*B és a szabályozó átviteli függvénye: ( ) a zár kör viselkedésén pedig látható a d=1 holtidő hatása (5.13. ábra):
5.13. ábra. A zárt kör viselkedése d=1 holtidő esetén, valamint az algoritmus átmeneti függvénye
Irodalomjegyzék [1] R. Isermann: Digital Control Systems, Springer-Verlag, 1989. [2] K. Dutton,S.Thompson:The Art of ControlEngineering,Addison-Wesley, 1997. [3] Lantos B.: Irányítási rendszerek elmélete és tervezése. Akadémiai Kiadó, Budapest, 2001.
113
IMC algoritmus
6. BELSŐ MODELLT ALKALMAZÓ (IMC) DIGITÁLIS SZABÁLYOZÁSI ALGORITMUSOK (Internal Model Control, IMC)
6.1. Az IMC algoritmusok alapelve A belső modellt alkalmazó szabályozási algoritmust (IMC) tárgyaljuk a 6.1. ábra alapján (az a struktúra átalakításával kapjuk a b struktúrát).
OBJEKTUM
OBJEKTUM
6.1. ábra. A megegyező a és b IMC struktúrák Folytonos esetben a szabályozó kimenő jelét a következő átviteli függvény alapján számíthatjuk ki (6.1):
u (s)
Gc e(s) 1 G c G om
A rendszer átviteli függvénye a zavarójelet is figyelembe véve (6.2): 114
(6.1)
IMC algoritmus
y(s)
G ov (1 G omG c ) GoGc v(s) r (s) 1 G c (G o G om ) 1 G c (G o G om )
(6.2)
Ha nem rendelkezünk a Go(s) objektum Gom(s) modelljével, azaz Gom(s)=0 akkor a 6.2. kifejezés egyszerűsödik (6.3):
y(s)
G ov GoGc v(s) r (s) 1 GcG o 1 GoGc
(6.3)
ahol a G c (s) szabályozási algoritmus lehet például egy PID algoritmus is. Egy másik esetben a Gom(s) átviteli függvény egyezzen meg a Go(s) átviteli függvénnyel, azaz a szabályozott objektum modelljével (6.4):
G om (s) G o (s)
(6.4)
Ezt az esetet feltételezve a rendszer kimenete (6.5):
y(s) G ov (1 G o G c )v(s) G o G c r(s)
(6.5)
Legyen a szabályozási rendszerünk zavarójel mentes (v(t)=0) és tételezzük fel, hogy az r(t) parancsolt érték változást az y(t) kimenet tökéletesen követi. Ebben az esetben
y(s) G om G c r(s)
(6.6)
(6.6) egyenletből kifejezzük a szabályozási algoritmust: Gc
1 , Go
(6.7)
amely nem lesz más, mint az objektum modelljének inverze (6.7). A (6.7.) egyenlettel megadott szabályozási algoritmusnak több problémája van. Az első az, hogy a Go objektum legtöbbször nemlineáris, általában nem pontosan modellezhető, illetve az 1/Go nem realizálható. Ezen utóbbira nézzünk meg egy egyszerű esetet mintavételezett rendszerek esetén. Legyen egy elsőrendű objektumunk d=0 diszkrét holtidővel. Ebben az esetben az objektum Go(z) átviteli függvénye (6.8):
b1 z 1 G o (z) 1 a 1 z 1
115
(6.8)
IMC algoritmus
A G c (z) a (6.7) egyenlet alapján a G o (z) reciproka lesz (6.9)
G c ( z)
1 a 1z 1 u (z) 1 G o (z) e(z) b1z 1
(6.9)
melyből a Gc(z) algoritmus differencia egyenlete nem írható fel, nem realizálható az algoritmus. A realizálhatóság érdekében a (6.9.) Gc(z) algoritmussal sorba kötünk egy egységnyi erősítési tényezőjű elsőrendű digitális szűrőt, melynek átviteli függvénye (6.10): G f (z)
(1 )z 1 1 z 1
(6.10)
így a kibővített szabályozási algoritmus (6.11):
G c ( z)
1 G f ( z) G o ( z)
(6.11)
lesz, mely már realizálható, hisz példánk szerint az elsőrendű rendszer esetén a szabályozási algoritmus (6.12):
1 a 1z 1 (1 )z 1 (1 ) a 1 (1 )z 1 G c (z) , b1 z 1 1 z 1 b1 b1z 1
(6.12)
tehát a 6.1. ábrán látható Gc szabályozási algoritmus nem egyszerűen a realizálhatatlan objektummodell inverze, hanem a Gf szűrővel bővített változat.
A Gf(z) szűrő fokszáma m=n esetén mindig egy lesz, ha a diszkrét holtidő d=0. A Gf szűrő fokszáma általános esetben d+1.
6.2. MATLAB program digitális IMC algoritmus tervezéshez % DIGITÁLIS IMC ALGORITMUS TERVEZÉSE % másodrendű holtidő nélküli objektumhoz % (Dr. Vass József, Pannon Egyetem) % To- mintavételezési idő, s 116
IMC algoritmus
To=4 % Az identifikált objektum impulzus-átviteli függvénye; To- mintavételezési idő Gom=tf([0.08232 0.0258],[1 -1.064 0.2881],To) % Az objektum B és A polinomjai [B,A]=tfdata(Gom,'v') step(Gom) pause % Szűrő folytonos modellje gf=tf(1,[10 1]) % Az elsőrendű szűrő diszkrét átviteli függvénye. Mintavételezési idő: 4 s. Gf=c2d(gf,4) % A szűrő polinomjai [r,s]=tfdata(Gf,'v') % Az objektum inverze Gomi=1/Gom % Az objektum inverze a szűrővel Gomif=Gomi*Gf % A Gomif polinomjai [R,S]=tfdata(Gomif,'v') step(Gomif) pause % IMC algoritmus Gimc=Gomif/(1-Gomif*Gom) % IMC algoritmus polinomjai [Q,P]=tfdata(Gimc,'v') % IMC algoritmus átmeneti függvénye 117
IMC algoritmus
step(Gimc) zpk(Gimc)
A fenti m-file a következő vektorokat szolgáltatja a 6.3. ábra a) struktúrájához: R = [0.3297 -0.3508 0.0950] S = [0.0823 -0.0294 -0.0173] A 6.3. ábra b.) struktúrájában a Q(z) és P(z) polinomok együtthatói a következők lesznek: Q =[0.0271 -0.0674 0.0613 -0.0211 -0.0015 0.0027 -0.0005] P =[0.0068 -0.0143 0.0076 0.0017 -0.0019 0.0000 0.0001] Természetesen a szabályozó algoritmus fokszáma nem hat lesz, ugyanis a megegyező zérusok- pulzusok kiejtik egymást. Ez látható a következő alakból: Zero/pole/gain: (
( )(
)( )(
)( ) (
) )
A step(Gimc) parancs megmutatja a megtervezett IMC algoritmust átmeneti függvényét (6.2. ábra):
6.2. ábra. A megtervezett IMC algoritmus átmeneti függvénye
118
IMC algoritmus
a 6.3. ábra bemutatja a SIMULINK-es modellt.
6.3. ábra. Az IMC algoritmust tartalmazó szabályozókör jelei
A 6.3. ábra SIMULINK modellje által szolgáltatott tervezés eredményét láthatjuk a Scope1 és Scope3 ábrákon. Az előbbi alapján látható, hogy a zárt kör viselkedése ideális esetben megegyezik a tervezésnél felhasznált szűrő időtartománybeli viselkedésével (6.4. ábra). y(t), y(k)
u(k)
6.4. ábra. A zárt kör viselkedése és a végrehajtó jel IMC algoritmusnál 119
IMC algoritmus
6.3.Esettanulmány Az 4. fejezethez tartozó esettanulmányban elvégeztük a QUANSER DC motor identifikációját. Az identifikáció eredménye a következő modell volt (To=0.01 s, d=1) th=oe(z,[3 3 2]) Discrete-time IDPOLY model: y(t) = [B(q)/A(q)]u(t) + e(t) B(q) = 2.009 q-2 - 3.575 q-3 + 1.897 q-4 A(q) = 1 - 2.603 q-1 + 2.489 q-2 - 0.8666 q-3 Az IMC algoritmus tervező programban lecseréltük a Gom(z) átviteli függvényt, illetve a Gf(z) szűrő fokszámát d+1-re választottuk, vagyis két elsőrendű szűrőt sorba kapcsoltunk (6.5. ábra). A folytonos elsőrendű szűrök időállandója 0.1 s.
6.5. ábra. A Gf(z) szűrő átmeneti függvénye % Az identifikált objektum impulzus-átviteli függvénye; To- mintavételezési idő Gom=tf([2.009 -3.575 1.897],[1 -2.603 2.489 -0.8666 0],To) % Az objektum B és A polinomjai [B,A]=tfdata(Gom,'v') step(Gom) pause % Szűrő folytonos modellje (másodrendű!!!)
120
IMC algoritmus
gf=tf(1,[0.1 1]) % Az elsőrendű szűrő diszkrét átviteli függvénye. Mintavételezési idő: 0.1 s. Gf=c2d(gf,To) % A másodrendű szűrő átviteli függvénye Gf=Gf*Gf % A szűrő polinomjai [r,s]=tfdata(Gf,'v') A rendszer frekvencia függvényeiből látható, hogy a megtervezett IMC algoritmus tartalmaz egy Noch-szűrőt is. A Bode-diagramot a következő paranccsal kapjuk meg: bode(Gom,Gimc,Gimc*Gom)
6.6. ábra. Az objektum, az IMC algoritmus és a nyitott kör frekvencia függvényei
A Noch-szűrőt is tartalmazó IMC algoritmus átmeneti függvényét az 6.7. ábra mutatja:
121
IMC algoritmus
6.7. ábra. Noch-szűrőt is tartalmazó IMC algoritmus átmeneti függvénye A zárt kör SIMULINK modellje, ill. r(k)=100 ford/perc alapjel esetén a tranziens folyamatok az 6.8. – 6.9. ábrákon láthatók.
6.8. ábra. A zárt kör SIMULINK modellje
6.9. ábra. A zárt kör átmeneti függvénye és végrehajtó jele 122
IMC algoritmus
Irodalomjegyzék [1] K. Dutton,S.Thompson:The Art of ControlEngineering,Addison-Wesley, 1997. [2] Lantos Béla: Irányítási rendszerek elmélete és tervezése. Akadémiai Kiadó, Budapest 2001.
123
SMITH-prediktor
7. NAGY HOLTIDŐVEL RENDELKEZŐ RENDSZEREK DIGITÁLIS SZABÁLYOZÓI
7.1. A holtidős rendszerek tulajdonságai Tekintsünk egy általános visszacsatolt rendszert. Tételezzük fel, hogy ebben a körben, az objektum jelentős holtidővel rendelkezik, domináns időállandója nagyságrendjébe esik. A nagy holtidő okai lehetnek például:
A folyamat lehet egy folyadékot, szilárd anyagot hosszú úton szállító rendszer, vagy előfordulhat a hosszú lappangási idő jelensége, vagy rossz gépészeti tervezés, vagy például egy hőmérséklet érzékelő kerámia szigetelése.
A mérőeszköznek hosszú időre lehet szüksége, hogy elvégezze a mintavételezést és a kimenet mért értékeinek az analízisét például az anyagösszetétel megállapítását. Ilyenek a folyadék kromatográfok, a szakaszosan működő összetétel analizátorok stb.
A fent említett esetekben egy egyszerű negatív visszacsatolás nem vezet eredményre, mert:
Az objektumot érő bemeneti zavaró jelek csak hosszú idő elteltével észlelhetők.
A beavatkozás, ami a legutolsó mérési eredményen alapszik, nem lesz megfelelő, mert olyan eseményekre állítja elő a beavatkozó jelet, amelyek a régmúltban történtek, így soha nem az aktuális kimenet kerül kiértékelésre a szabályozási algoritmusban. A fent említett okok miatt a nagy mértékű holtidő az egyik fő oka lehet a visszacsatolt
rendszer instabilitásának, a lengési tulajdonságok felerősödésének, melyet ezen holtidős rendszerek nyitott köri frekvencia függvényei is bizonyítanak, azaz a Nyquist, illetve Bode stabilitási kritériumok. Az 7.1. ábrán egy folytonos átviteli függvénnyel rendelkező rendszer:
G 0 (s)
y(s) 5e Th s u (s) 0,5s 1
124
SMITH-prediktor
Bode diagramját láthatjuk különböző holtidők esetén. Természetesen a mintavételezett rendszerek esetén is érvényesek az ábrán vázolt holtidő hatások. A G0(s) rendszer Bodediagramját Th=0.01, 0.1, 0.2 s holtidő esetén ábrázoltuk. Az amplitúdó jelleggörbe mindhárom esetben megegyezik, azonban a fázismenetek különbözőek. Bizonyítható az, hogy a kritikus erősítési tényező a Th holtidő növelésével drasztikusan csökken. A stabilitás és fázistartalék biztosításához a szabályozó erősítési tényezőjét csökkenteni kell, azonban ezzel nagy beállási idejű, kúszó, lomha szabályozási rendszert kapunk (7.1. ábra).
7.1. ábra. A holtidő és a stabilitás kapcsolata (nyitott kör Bode diagramjai)
7.2. A SMITH prediktor A holtidő kompenzáció egyik megoldása a SMITH prediktor (holtidő kompenzátor) alkalmazása (7.2. ábra).
125
SMITH-prediktor
7.2. ábra. A SMITH prediktor elhelyezése a szabályozási körben Legyen az objektum modellje (7.1)
y(z) B(z 1 ) d G 0 (z) z u (z) A(z 1 )
(7.1)
ahol d a diszkrét holtidő ( d = Th / T0), valamint a szabályozási algoritmus impulzus-átviteli függvénye Gc(z). Az u(k) végrehajtójelet csatoljuk negatívan vissza egy G 0 (z) G 0 z d taggal az összehasonlító szerven keresztül (7.2, 7.3). Itt
G 0 (z)
G 0 (z)z d
y(z) B(z 1 ) u (z) A(z 1 )
(7.2)
y(z) B(z 1 ) d z u (z) A(z 1 )
(7.3)
Ebben az esetben az u(k) és r(k) közti impulzus-átviteli függvény a következő lesz (7.4). u(z) G c (z) r(z) 1 G cG 0 G cG 0z d
(7.4)
A SMITH prediktort tartalmazó nyitott kör átviteli függvénye (7.5):
y(t) G c (z)G 0 (z)z d r(z) 1 G c (z)G 0 (z) G c (z)G 0 (z)z d
126
(7.5)
SMITH-prediktor
a zárt kör eredő átviteli függvénye a következő lesz (7.6):
G r (z)
G c G 0 d y( z ) z r ( z) 1 G c G 0
(7.6)
melyből az következik, hogy a SMITH prediktort alkalmazva zárt kör stabilitását a d holtidő nem befolyásolja és a visszacsatolás megtervezhető a holtidő nélküli objektum alapján (7.3. ábra).
7.3. ábra. A prediktort tartalmazó rendszer eredő hatásvázlata
Példa Legyen egy folytonos rendszer munkaponti átviteli függvénye: G 0 (s)
0,25e 30s , 270s 2 40s 1
melynek átmeneti függvényét az 7.4. ábra szemlélteti:
7.4. ábra. A mintavételezési idő meghatározása 127
SMITH-prediktor
Válasszuk a mintavételezési T0 időt a domináns időállandó egytizedére (T0=T/10), mely T0=5s lesz. A SMITH prediktor tervezéséhez szükségünk van a holtidő nélküli objektum diszkrét modelljére, mely:
G 0 (z)
0,008819z 1 0,006643z 2 , 1 1,365z 1 0,4266z 2
illetve a holtidővel rendelkező modellre, ahol d =Th /T0 = 30/5=6
G 0 (z)z d
0,008819z 7 0,006643z 8 . 1 1,365z 1 0,4266z 2
Megterveztük a holtidő nélküli rendszer számára a digitális PI algoritmust, mely a következő: G c ( z)
5 4,3z 1 . 1 z 1
A holtidő nélküli rendszert visszacsatolva a fenti PI algoritmussal a zárt kör tranziens viselkedése (7.5. ábra):
7.5. ábra. A holtidő nélküli zárt kör viselkedése Érdekesség kedvéért nézzük meg a zárt szabályozási kör viselkedését a holtidővel rendelkező objektum visszacsatolásakor. Látható az 7.6. ábrán, hogy a holtidő a fázistartalékot erősen csökkentette.
7.6. ábra. A holtidővel rendelkező visszacsatolt rendszer viselkedése SMITH prediktor nélkül 128
SMITH-prediktor
Most hozzuk létre a 7.2. ábra szerinti SMITH prediktort tartalmazó visszacsatolt rendszert, és vizsgáljuk meg viselkedését változó holtidők esetén. A vizsgált rendszer MATLAB/SIMULINK felületét mutatja a 7.7. ábra:
7.7. ábra. A SMITH prediktor elhelyezése SIMULINK felületen A szimulációs vizsgálatok is bizonyítják, hogy SMITH prediktort alkalmazva holtidőtől független visszacsatolás tervezést végezhetünk. Az alábbi ábrán To=30 s és To=50 s holtidő esetén vizsgáljuk a zárt rendszer viselkedését. Láthatjuk, hogy a zárt szabályozási körök átmeneti függvényei csak holtidőben különböznek egymástól (7.8. ábra).
7.8. ábra. A zárt kör viselkedése különböző holtidők esetében A SMITH prediktorral tervezett nagy holtidővel rendelkező szabályozási rendszerek megfelelő pontosságú identifikáció esetén robosztus tulajdonsággal rendelkeznek.
129
SMITH-prediktor
Irodalomjegyzék [1] Lantos Béla: Irányítási rendszerek elmélete és tervezése. Akadémiai Kiadó, Budapest 2001.
130
MVC algoritmus
8. MINIMÁLIS SZÓRÁSÚ SZABÁLYOZÁSI ALGORITMUSOK A folytonos technológiai rendszerek általában huzamos ideig valamilyen állandó alapjellel üzemelnek, a szabályozó algoritmusok fő feladata a rendszert érő zavarások hatásának csökkentése a kimeneten. Nagyszerű példa erre egy hajó GPS-sel és ún. aktuátorral ellátott iránytartó navigációs rendszere, amely számára megadhatjuk a cél pontos koordinátáit (alapjel). Abban az esetben, ha nincs szél, nincs áramlás, hullámzás, a szabályozó rendszer nélkül is pontosan elérjük a kívánt végpozíciót, az utazási célt. Ha a fent említett zavarások léteznek, akkor csak visszacsatolt rendszer biztosíthatja a megfelelő irányítást. Természetesen ebben az esetben az a jó szabályozási stratégia, ha a kívánt/parancsolt útvonalat a kiindulási és a végcél pontja között minimális eltérésekkel, minimális szórással tudjuk megtenni. A zavaró jelek véletlenszerűek, nem mérjük őket. Ilyen és ehhez hasonló feltételek esetén javasolható visszacsatolásként az ún. minimális szórású szabályozó algoritmusok használata. A következőkben nézzünk egy egyszerű technológiai folyamatot, egy hőcserélő hőmérséklet szabályozását (8.1. ábra). A hőcserélő fűtőenergiáját egy több felhasználót ellátó gőzvezeték biztosítja, ahol a gőzfelhasználók véletlenszerűen terhelik a gőzellátó rendszert, melynek nyomása n(t) sztochasztikus jelnek tekinthető és nem mérjük. A melegítendő folyadék mennyisége n2(t) vagy esetleg belépő hőmérséklete n3(t) szintén véletlen – nem mért jelnek tekinthető. A hőmérsékletszabályozó TC legyen digitális, melynek bemenete y(k) jelsorozat, alapjele az r(k), kimenete u(k) jelsorozat. Tervezzünk
visszacsatolásként
ún.
minimális
algoritmust.
131
szórású
digitális
szabályozási
MVC algoritmus
8.1. ábra. Hőcserélő műszerezése és kapcsolata a számítógéppel
8.1. Minimális szórású szabályozók holtidő nélküli rendszerekhez A minimális szórású szabályozási algoritmusok állandósult állapotban /stacionárius esetben/ a szabályozási hiba négyzetének várható értékét minimumra csökkentik (8.1):
E y 2 (k) min
(8.1)
Általános esetben az a célfüggvény, amelynek minimumát keressük a minimális szórású irányítási stratégiánál w(k) = 0 és e(k) = -y(k) esetre a következő lesz (8.2):
I(k 1) E y 2 (k 1) ru 2 (k)
(8.2)
Holtidőt nem tartalmazó rendszerek esetén, mivel az u(k) bemenőjel b0 = 0 esetén (amely a valós rendszereknél általában fennáll) az y(k+1)/ kimenőjelet generálja (8.2. ábra).
132
MVC algoritmus
8.2. ábra. A zárt szabályozó kör minimális szórású szabályozási (MVC) algoritmussal
Az irányítási algoritmus megtervezése feltételezi, hogy a szabályozandó folyamat identifikált, tehát struktúrája és paraméterei ismertek és a következő átviteli függvényekkel jellemezhető (8.3):
y(z) b1 z 1 b 2 z 2 ..... b m z m G o (z) u (z) 1 a 1 z 1 .... a m z m
(8.3)
a zajcsatorna bemenete, v(k) nem mérhető, de a csatorna identifikálható és az identifikáció eredményeképpen a következő modellel jellemezhető (8.4):
G ov (z)
n(z) 1 d1z 1 ..... d m z m , v(z) 1 c1z 1 .... cmz m
(8.4)
ahol aav(k) zavarójel amplitúdója. A fenti ábra szerinti elrendezésben a v(k) zérus várható értékű, független jel, fehérzaj, melyre igazak a következő kikötések (8.5):
1, ha 0 E v(k).v(k ) 0, ha 0
Ev(k) 0
(8.5)
A teljes objektummodellt a következőképpen írhatjuk fel (8.6) [1]: ( )
(
)
(
)
( ) λ
(
)
(
)
( )
(8.6)
A célfüggvényben szereplő y(k+1) kimenőjel előre meghatározható az y(k), y(k-1, ... és u(k),u(k-1) értékekből, ahol ezt az un. predikciót a z transzformáció tulajdonságait felhasználva megadhatjuk. Az egyenletet átrendezve az 133
MVC algoritmus
A(z 1 )C(z 1 ) y(z)z B(z 1 )C(z 1 )u(z)z A(z 1 )D(z 1 )v(z)z
(8.7)
transzformációs egyenletet (8.7) kapjuk, ahol az egyes nagybetűvel jelölt polinomokat tagonként behelyettesítjük és a z-vel történő szorzást figyelembe vesszük, a következő differenciaegyenlethez jutunk (8.8):
y(k 1) (a 1 c1 ) y(k ) ... a m c m y(k 2m 1) b1 u (k ) (b 2 b1c1 )u (k 1) ... b m c m u (k 2m 1)
(8.8)
v(k 1) (a 1 d1 ) v(k ) .... a m d m v(k 2m 1) Az előző egyenletből kifejezve a célfüggvény számára szükséges y(k+1)predikciót és behelyettesítve a célfüggvény egyenletébe, majd a I(k 1) 0 u (k )
(8.9)
feltételt megvizsgálva (8.9) a minimális szórású szabályozási algoritmus általános egyenletét kapjuk (8.10) [1]: ( )
( )
(
)
( )
(
)
( (
)
(
) (
) )
( (
) ) (
)
(8.10)
ahol az MVC angol rövidítés (Minimum Variance Controller), b1 pedig a B polinom első együtthatója. A fenti egyenlet az ismert paraméterű rendszerek általános minimális szórású algoritmusa, amely nagymértékben egyszerűsödik úgy a zajmodell egyszerűsítése, mint a célfüggvény megváltoztatása miatt. Legyen az I(k+1) célfüggvényben a végrehajtójel amplitúdóját korlátozó faktor r = 0. Akkor a minimális szórású algoritmus (8.11): ( )
(
)
(
)
(
(
) (
) )
(8.11)
Abban az esetben, ha az objektum és a zajmodell karakterisztikus egyenletei megegyeznek, tehát C(z-1) = A(z-1) és r nem egyenlő 0-val az algoritmus (8.12): ( )
(
)
(
)
(
)
(
) (
)
(8.12)
és végül a C(z-1) = A(z-1) és r =0 esetre (8.13): ( )
(
134
(
) )
(8.13)
MVC algoritmus
8.1.1. Stabilitás Az irányított folyamat, illetve annak modelljének egyes tulajdonságai a minimális szórású algoritmusok alkalmazhatóságát illetve a zárt szabályozási kör stabilitását befolyásolják. A zárt kör stabilitására a több kikötést kell tennünk [1]. A zárt kör karakterisztikus egyenlete (8.14):
1 G c (z 1 )G 0 (z 1 ) 0
(8.14)
r A(z) A(z) zB(z) D(z) 0 b1
(8.15)
r A(z) zB(z) D(z) 0 b1
(8.16)
MV1-re (8.15):
ésMV3-ra (8.16):
A C(z) polinomot mindkét esetben elhagyjuk. Ebből következik, hogy a zárt kör stabilitása: MV1 és MV3 ( r 0 ) algoritmusok esetén: A zajcsatorna zérusainak D(z)=0 az egységkörön belül kell lennie a z tartományban. Az
r A(z) zB(z) 0 b1
(8.17)
zérusainak az egységkörön belül kell lenniük (8.17). MV1-re a rendszer pólusainak az egységkörön belül kell lenniük. MV2 és MV4 (r=0) algoritmusok esetén: A zárt kör karakterisztikus egyenlete r=0-ra (8.18): MV2 : zA(z)B(z)D(z) 0 MV4 : zB(z)D(z) 0
(8.18)
Ezért a rendszer B(z)=0 és a zajcsatorna D(z)=0 zérusainak az egységkörön belül kell lenniük. A zajcsatorna pólusai C(z)=0 (MV2-re) nem befolyásolják a karakterisztikus egyenletet, így akárhol lehetnek.
135
MVC algoritmus
Példa Egy mintapélda segítségével (Függelék 1) bemutatjuk a holtidő nélküli rendszer viselkedését MVC algoritmust alkalmazva. Az objektum átmeneti függvényei (8.3.-8.4. ábra): step(Go)
8.3. ábra. Az alapcsatorna átmeneti függvénye step(Gov)
8.4. ábra. A zajcsatorna átmeneti függvénye A minimális szórású algoritmus:
136
MVC algoritmus
Az algoritmus átmeneti függvénye (8.5. ábra):
8.5. ábra. Az MVC algoritmus átmeneti függvénye Az algoritmus zérusai, pólusai (8.6. ábra):
8.6. ábra. Az MVC algoritmus zérusai, pólusai
A zárt kör SIMULINK modellje (8.7. ábra):
137
MVC algoritmus
8.7. ábra. A zárt kör SIMULINK modellje és viselkedése r(k)=0 estben a következő lesz (8.8. ábra).
8.8. ábra. A zárt kör kimenete és a végrehajtó jel Mivel a zárt kör integráló tulajdonságú elemet nem tartalmaz, a tervezett rendszer maradó szabályozási eltéréssel képes csak üzemelni a parancsolt érték nem nulla (pl. r(k)=2) esetben (8.9. ábra).
8.9. ábra. A zárt kör maradó szabályozási eltérése és végrehajtó jele 138
MVC algoritmus
A maradó szabályozási eltérés megszüntetéséhez szükségünk van egy +1 helyen elhelyezkedő valós pólusra, melyet egy digitális integrátor beépítésével valósíthatunk meg (8.10. ábra). A zárt rendszer hatásvázlata így megváltozik:
8.10. ábra. Integrátor beépítése Természetesen az integrátor beépítésével a minimális szórás megnő, azonban az átlagérték a parancsolt érték lesz (8.11. ábra):
8.11. ábra. A beépített integrátor hatása
139
MVC algoritmus
8.2. Minimális szórású algoritmusok holtidős rendszerek esetén Holtidős rendszerek esetén az egy bemenetű és egy kimenetű rendszer impulzus-átviteli függvénye a holtidővel módosul (8.19):
G 0 (z)
b1z 1 ..... b m z m d y(z) B(z 1 ) d z z u (z) A(z 1 ) 1 a 1z 1 .... a m z m
(8.19)
a rendszer hatásvázlata pedig az ábrának megfelelő marad. Az irányított rendszer holtideje miatt az u(k) beavatkozójel a kimenetre k+d+1 időpillanatban hat, tehát a célfüggvény felépítése is módosul (8.20):
I(k 1) E y 2 (k d 1) ru 2 (k)
(8.20)
a rendszermodell a (8.21) lesz:
y( z )
B(z 1 ) d D(z 1 ) z u ( z ) v( z ) A(z 1 ) C(z 1 )
(8.21)
A célfüggvény felépítéséhez szükséges y(k+d+1) tag kifejezésére felírhatjuk (8.22):
z(d1) y(z)
B(z 1 ) D(z 1 ) (d1) zu(z) z v(z) A(z 1 ) C(z 1 )
(8.22)
Az egyenletből következik, hogy a v(k), v(k-1) megfigyelések pontosan számíthatók a v(k+1),v(k+2)...stb azonban nem, ezért célszerű ezeket a megfigyeléseket szétválasztani a következő módon (8.23):
z
( d 1)
1 ( d 1) L(z 1 ) B(z 1 ) y( z ) zu(z) F(z )z v( z ) A(z 1 ) C(z 1 )
(8.23)
ahol az egyes polinomok az együtthatók összehasonlítási módszerével egyszerűen számolhatók a
D(z 1 ) F(z 1 )C(z 1 ) z (d1) L(z 1 )
140
(8.24)
MVC algoritmus
8.12. ábra. A zajcsatorna helyettesítése holtidős objektum esetén összefüggésből (8.24), ahol az F és L polinomok a következő alakúak (8.25) [1]: F(z 1 ) 1 f1z 1 f 2 z 2 ... f d z d L(z 1 ) 10 11 z 1 ...1m1 z ( m1)
A
(8.25)
egyenletet
időtartományra
áttranszformálva
(8.25) és
behelyettesítve
az
integrálkritériumba, majd kiszámítva a I(k 1) 0 u (k )
(8.26)
szélsőértéket megkapjuk a holtidős rendszerek esetére (8.26) a minimális szórású szabályozási algoritmus alakjait. Általános esetben, mikor A(z-1) nem egyenlő C(z-1) és r nem egyenlő 0, az algoritmus (8.27) [1]. d G 1MVC ( z)
A(z 1 )L(z 1 ) r zB(z 1 )C(z 1 )F(z 1 ) A(z 1 )D(z 1 ) b1
(8.27)
ha pedig A(z-1) nem egyenlő C(z-1) és r=0, akkor (8.28): d G 2MVC (z)
A(z 1 )L(z 1 ) zB(z 1 )C(z 1 )F(z 1 ) 141
(8.28)
MVC algoritmus
valamint A(z-1) = C(z-1) és r nem egyenlő 0 (8.29): d G 3MVC (z)
L(z 1 ) r zB(z )F(z ) D(z 1 ) b1 1
(8.29)
1
és r=0 esetben (8.30): d G 4MVC ( z)
L(z 1 ) zB(z 1 )F(z 1 )
(8.30)
A zárt szabályozási kör stabilitására ugyanazok a szabályok érvényesek, mint a holtidő nélküli rendszerek esetében. Példa Tervezzünk minimális szórású algoritmust holtidős rendszer számára. Legyen az alapcsatorna folytonos átviteli függvénye
G0 ( s)
y( s) e4 s , u ( s) 37,5s 2 12,5s 1
és a mintavételezési időt válasszuk T0=4s-ra [1]. A mintavételezett alapcsatorna impulzus átviteli függvénye tehát:
y( z ) B( z 1 ) 0,1387 z 1 0, 0889 z 2 1 G0 ( z ) z. u ( z ) A( z 1 ) 1 1, 036 z 1 0, 2636 z 2 A zajcsatorna átviteli függvénye:
Gov ( z )
n( z ) D( z 1 ) 1 0,5 z 1 0, 25 z 2 z1 . 1 1 2 v( z ) C ( z ) 1 1, 036 z 0, 2636 z
Láthatjuk, hogy a Go(z) és Gov(z)-ből származtatható karakterisztikus egyenletek megegyeznek: A(z-1)=C(z-1) és legyen a végrehajtó jel büntető faktora: r=0, az objektum fokszáma: m=2. Ebben az esetben a (8.30) egyenletet alkalmazva a
Gc ( z ) G
4d MVC
u( z) L( z 1 ) ( z) e( z ) zB( z 1 )( Fz 1 )
minimális szórású algoritmust kell használni. Az L(z-1), F(z-1) polinomok számítására felhasználjuk a (8.22) egyenletet
D( z 1 ) F ( z 1 )C ( z 1 ) L( z 1 ) z ( d 1) 142
MVC algoritmus
összefüggést, ahol m=2 esetén:
F ( z 1 ) 1 f1 z 1 L( z 1 ) l0 l1 z 1
.
A fenti konkrét rendszert tekintve behelyettesítjük: 1 0,5z 1 0, 25z 2 (1 f1 z 1 )(1 1,036 z 1 0, 2636 z 2 ) (l0 l1 z 1 ) z 2
majd kiszámolhatóak a F és L polinom együtthatói:
f1 0,5 1, 036 1,536 l0 0, 25 0, 2636 1, 036 f1 1,5777 l1 0, 2636 1,536 0, 4049 azaz a keresett polinomok
F ( z 1 ) 1 1,536 z 1 L( z 1 ) 1,5777 0, 4049 z 1 a szabályozási algoritmus pedig
1,5777 0, 4049 z 1 Gc ( z ) 0,1387 0,3019 z 1 0,13655 z 2 lesz. A 8. 13. ábrán a visszacsatolt rendszer SIMULINK modelljét mutatjuk be,
8.13. ábra. Holtidős rendszer MVC algoritmussal majd a szabályozó kör kimenetének viselkedését (8.14. ábra):
143
MVC algoritmus
y(k)
y(t)
8.14. ábra. Az MVC algoritmust tartalmazó zárt kör viselkedése
melyhez az alábbi végrehajtójel tartozik (8.15. ábra): u(k)
8.15. ábra. Az MVC algoritmus kimenete
Megjegyezzük, hogy a holtidővel rendelkező rendszerek esetén is maradó szabályozási eltérés keletkezik az integráló tag hiánya miatt. Példánkban a parancsolt érték: r(k)=2 volt. A maradó szabályozási eltérés kiküszöbölésére az 1. példa szerint járunk el.
Irodalomjegyzék [1] R. Isermann: Digital Control Systems, Springer-Verlag, 1989. [2] Vass J.:Programcsomag digitális szabályozási rendszerek tervezéséhez. Mérés és Automatika, 36. évf., 1988. 1. szám 144
MVC algoritmus
Függelék 1. MATLAB program G 2MVC (z) tervezéséhez % MINIMÁLIS SZÓRÁSÚ SZABÁLYOZÁSI algoritmus tervezése % PANNON EGYETEM % Dr.Vass József % Másodrendű Box-Jenkins modell, d=0, r=0, Gmvc2. % Az identifikált modell polinom együtthatói; To=4 s. B=[0 .1387 .0889] A=[1 -1.036 .2636] C=[1 -1.138 .3536] D=[1 0.5 0.25] Go=tf(B,A,4) Gov=tf(D,C,4) % A 8.11 képletet felhasználva Goi=1/Go G=Gov-1 Gmvc=Goi*G [Q,P]=tfdata(Gmvc,'v') step(Gmvc)
145
Előrecsatolások tervezése
9. ELŐRECSATOLÁSOK TERVEZÉSE Az eddigi fejezetek szabályozási algoritmusai olyan esetekre vonatkoztak, amikor az ellenőrző jel mérése után a számítógépes algoritmus negatív visszacsatolásként a szabályozott objektum bemenetére hatott, tehát visszacsatolt rendszerekkel foglalkoztunk. A visszacsatolt rendszerek mellett az előrecsatolásos algoritmusok is előkelő helyet foglalhatnak el egyes szabályozástechnikai problémák megoldásánál. Tételezzük fel, hogy valamely objektum kimenetére zavaró jellemző szuperponálódik. Ez a zavaró jellemző természetesen lehet determinisztikus, vagy sztochasztikus jellegű és nagyon sok esetben közvetlenül mérhető is. Pl. egy gőzös hőcserélőnél a szabályozott jellemző lehet a kilépő folyadékáram hőmérséklete, zavaró jellemzője pedig a gőz nyomása, amely a különböző fogyasztók terhelései miatt erősen ingadozik és zavaró jellemzőként lép fel (ld. 2. fejezetben). Ilyen esetekben a rendszer struktúráját a 9.1. ábra alapján képzelhetjük el, ahol a Ge(z) előrecsatolt szabályozó algoritmus a v(k)-ban fellépő ingadozások alapján Gc(z) visszacsatolt szabályozó kimenetét megváltoztatja, nem várva meg azt, hogy v(k) változásai az y(k) kimenőjelben érvényesüljenek. Ideális esetben olyan Ge(z) előrecsatolást is meg lehet valósítani, hogy a v(k) bármilyen viselkedését az y(k)-ban kompenzálni lehet. A gyakorlatban determinisztikus jelek esetében paraméteroptimalizált előrecsatolásokat vagy sztochasztikus bemenetek esetén minimális szórású előrecsatolt algoritmusokat is alkalmaznak.
9.1. ábra. Az előrecsatolt szabályozó algoritmusok elhelyezkedése a hatásvázlatban 146
Előrecsatolások tervezése
A 9.1. ábra szerinti hatásvázlat egyes elemeit diszkrét esetben a következő impulzusátviteli függvényekkel jellemezhetjük (9.1, 9.2):
Go ( z )
y ( z ) B( z 1 ) b1z 1 ... bm z m d z u ( z ) A( z 1 ) 1 a1z 1 ... am z m
(9.1)
n( z ) D( z 1 ) d0 d1 z 1 ... dv z v v( z ) C ( z 1 ) 1 c1 z 1 ... cv z v
(9.2)
Gov ( z )
melyeket a további tervezési eljárásoknál ebben az alakban fogunk használni. A zajcsatorna identifikálása során olyan modellt kell választani, ahol d0=0 (pld. OE modellt)
9.1. Dinamikus előrecsatolás algoritmusa
Az ideális előrecsatolási algoritmus esetén a Gov(z) zajcsatorna szuperponáltn(k) jelét az előrecsatolt Ge(z) szabályozónak és a Go(z) objektumnak teljesen kompenzálnia kellene, amely azt jelenti, hogy ezek átvitele megegyezik, tehát érvényes a következő összefüggés (9.3):
Ge ( z )Go ( z ) Gov ( z ) 0
(9.3)
amelyből kifejezve a Ge(z) előrecsatolást a következő algoritmust kapjuk (9.4):
uv ( z ) Gov ( z ) A( z 1 ) D( z 1 ) K ( z) Ge ( z ) 1 1 d v( z ) Go ( z ) C ( z ) B( z ) z H ( z)
(9.4)
Ez az algoritmus a v(k) zavaró jellemzőt ideálisan kompenzálja, s egyszerűbb esetben, amikor a zajcsatorna impulzus-átviteli függvényében C(z-1)=A(z-1) a következőképpen egyszerűsödik (9.5): ( )
( (
) )
(9.5)
Az előrecsatolt szabályozó számlálójának fokszáma a 9.4. egyenlet alapján (m+v), nevezőjéé pedig (m+v+d). A 9.1. ábra alapján üzemelő rendszer eredő átviteli függvénye r(k)=0 alapjel esetén egyszerű összefüggéssel leírható (9.6): 147
Előrecsatolások tervezése
y ( z ) Ge ( z 1 )Go ( z 1 ) Gov ( z 1 ) Gve ( z ) v( z ) 1 GC ( z 1 )Go ( z 1 )
(9.6)
alakban, mely azt mutatja, hogy az előrecsatolt szabályozó algoritmus a szabályozott rendszer stabilitását nem befolyásolja, hisz karakterisztikus egyenletében (1+GCG0=0) az előrecsatolás nem játszik szerepet. A (9.4.) egyenlet általános alakja legyen (9.7):
Ge ( z )
hd z d
k0 k1 z 1 ... kmv z ( mv ) K ( z 1 ) h1 d z (1 d ) ... hmv d z ( mv d ) H ( z 1 )
(9.7)
mely alapján belátható, hogy d0 esetén az ideális előrecsatolás a gyakorlatban nem realizálható, valamint a következő kikötések érvényesek: a H(z) polinom gyökeinek ki kell elégíteni a (zi) < 1 feltételt, azaz a B és C polinomok gyökeinek az egységsugarú körön belül kell lenni a folyamatmodell és a zajcsatorna zérus helyeinek az egységsugarú körön belül kell lenni.
9.2. Előrecsatolási algoritmusok és visszacsatolások együttes használata Úgy a folytonos rendszereknél, mint a digitális irányítási rendszereknél az előrecsatolásokat önmagukban ritkábban alkalmazzák, felhasználásuk olyan szabályozási rendszereknél gyakoribb, amikor a rendszerben visszacsatolást is alkalmazunk (9.1. ábra). Természetesen a tervezésnek mindkét részletre ki kell terjednie: a visszacsatolt algoritmusra és az előrecsatolásra is. A visszacsatolt algoritmus (GC) az előző fejezetekben bemutatott módszerekkel megtervezhető, figyelembe véve azt a követelményt, hogy a tervezést alapjel váltásra, esetleg zavaró jellemzőre kell-e elvégeznünk. Példa Példaként vegyük a
y *( z ) 0, 087 z 1 0, 073z 2 B( z 1 ) Go ( z ) u( z) 1 1, 4 z 1 0, 48 z 2 A( z 1 )
148
Előrecsatolások tervezése
impulzus-átviteli függvénnyel leírható objektumot, melynél a mintavételezési idő T0=4 s, erősítési tényezője pedig: K0 2 . Az objektum kimenetére egy n(k) zavaró jellemző szuperponálódik, mely a következő zajcsatornán keresztül hat a kimenetre:
n( z ) 0,1 0,108 z 1 0, 032 z 2 D( z 1 ) Gov ( z ) v( z ) 1 1, 4 z 1 0, 48 z 2 A( z 1 ) A zajcsatorna erősítési tényezője pedig K0 3 . Előrecsatolást nem alkalmazva a GC szabályozóval visszacsatolt rendszer eredő impulzus-átviteli függvénye:
y( z ) D( z 1 ) P( z 1 ) Gv ( z ) v( z ) A( z 1 ) P( z 1 ) Q( z 1 ) B( z 1 ) ahol a visszacsatolt szabályozó polinomjainak meghatározását – az idővel súlyozott abszolút érték kritériumot választva – elvégeztük, s a PI szabályozóalgoritmus a következő lett:
uo ( z ) Q( z 1 ) 0, 41988 0,329733z 1 GC ( z ) e( z ) P( z 1 ) 1 z 1 Ezzel a szabályozó algoritmussal üzemeltetve a rendszert egységugrás zavarójel bemenet v(200)=1 és r(k)=1 alapjel váltás esetében a tranziens viselkedést a 9.2. ábra mutatja, ahol a zajhatás mértéke a kimeneten 1,5 egység.
9.2. ábra. Előrecsatolás nélküli visszacsatolt rendszer viselkedése zavaró jel esetén Alkalmazzuk a paraméteroptimalizált előrecsatolás legegyszerűbb esetét, a statikus előrecsatolást, amikor:
149
Előrecsatolások tervezése
( )
(9.8)
A (9.8) algoritmusra érvényes az, hogy erősítési tényezőjét (9.9) a következőképpen kell meghatározni: (9.9) ahol a Kov a zajcsatorna-, a K0 pedig az alapcsatorna erősítési tényezője. Ez a legegyszerűbb előrecsatolási algoritmus és statikus előrecsatolásnak nevezzük. A statikus előrecsatolás hatása akkor érzékelhető dominánsan, ha a zajcsatorna és az alapcsatorna dinamikus tulajdonságai hasonlóak és ezen csatornák holtideje megegyezik. A konkrét példánkban az előrecsatolás statikus algoritmusa tehát: Ge ( z )
3 1,5 . 2
Ezzel a statikus előrecsatolással kiegészítve a zárt szabályozási kört, a következő tranziens viselkedést kapjuk (9.3. ábra), ahol a v(200)=1 zavaró jel hatása sokkal kisebb amplitúdó változást okoz.
9.3. ábra. A zárt kör viselkedése statikus előrecsatolással Példa Ebben a feladatban a dinamikus előrecsatolás hatását mutatjuk be. Az objektum modelljeiből 9.4 egyenlet alapján számított dinamikus előrecsatolás beillesztését mutatja az 9.4. ábra:
150
Előrecsatolások tervezése
9.4. ábra. A dinamikus előrecsatolás beillesztése valamint a rendszer jeleit láthatjuk időtartományban. Az n(k) jelsorozatot ábrázoljuk a felső 9.5 részábrán, az alsó részábra pedig az y(k) kimenetet mutatja. A v(k) egy sztochasztikus zajforrás jele. (A 9.5. ábrán az y(k)ordinátája 0.0001-el szorzódik):
9.5. ábra. A dinamikus előrecsatolás hatása Ugyanez a rendszer statikus előrecsatolással a következő viselkedést mutatja (9.6. ábra):
9.6. ábra. A statikus előrecsatolás hatása
151
Előrecsatolások tervezése
A dinamikus előrecsatolással és a diszkrét PI algoritmussal történt visszacsatolással ellátott rendszer SIMULINK modelljét megalkotva v(k)=1 alapjel esetén a rendszer:
9.7. ábra. Előrecsatolás és visszacsatolás együttes alkalmazása viselkedését láthatjuk (9.7. ábra), ahol az n(k)és az y(k) jeleket ábrázoltuk (9.8. ábra).
9.8. ábra. A zárt kör veselkedése dinamikus előrecsatolással
152
MIMO tervezés
10. TÖBBVÁLTOZÓS RENDSZEREK DIGITÁLIS SZABÁLYOZÁSA (Multiple Input Multiple Output – MIMO – rendszerek)
10.1. A többváltozós rendszerek struktúrái A legtöbb gyártási folyamat az olajiparban, vegyiparban, un. többváltozós - MIMO - rendszer. A bonyolult repüléstechnikai eszközök, vagy az autók is a MIMO rendszerek struktúrájával rendelkeznek (10.1. ábra)
10.1. ábra. Autó, mint MIMO rendszer Egy autó esetében a rendszer bemenetei lehetnek például a gázpedál-, a fékpedál- és a kormánykerék pozíciók, kimenetek a megtett út, a kocsi sebessége, távolság egy balra/jobbra ívelő kanyarig. Zavaró jellemzőként tekinthető az úttest felületének minősége, stb. Úgy az un bemenetek, mint a vr zavarójelek hathatnak az ymkimenő jelekre. Az ilyen rendszerek általános felépítését látjuk a 10.2. ábrán.
10.2. ábra. A MIMO rendszerek általános felépítése
153
MIMO tervezés
Másik példaként vehetjük a legegyszerűbb desztilláló oszlop struktúráját (10.3. ábra). Bemeneti változók ebben az esetben a gőzzel bevitt energia mennyisége és a visszavezetés (reflux) mennyisége. Mindkettő befolyással van úgy a desztillátum, mint a fenéktermék koncentrációjára.
10.3. ábra. Desztilláló oszlop Egy desztilláló oszlop munkaponti 2x2-es MIMO modelljét [1] ábrázolja a 10.4. ábra SIMULINK felületen.
10.4. ábra. A desztilláló oszlop munkaponti SIMULINK modellje
154
MIMO tervezés
Az ilyen struktúrával rendelkező 2x2-es MIMO rendszert P kanonikus struktúrájúnak nevezzük, általános hatásvázlata a 10.5. ábrán látható.
10.5. ábra. Egy 2x2-es MIMO objektum P kanonikus struktúrával Vannak olyan technológiai rendszerek, ahol a MIMO struktúrában a kereszthatások más struktúrában jelentkeznek. A 10.6. ábrán gáz-folyadék elválasztó nyomás és szint kimenete hat egymásra egy egész más kölcsönhatásban. Ezt a struktúrát V kanonikus struktúrának (10.7. ábra) nevezzük. Természetesen a struktúrák keveredhetnek, de az megnehezíti a tervezési folyamatot.
10.6. ábra. Gáz-folyadék elválasztó (szeparátor)
155
MIMO tervezés
Ebben a fejezetben részletesen egy P struktúrájú MIMO objektum digitális szabályozásának
tervezésével
foglalkozunk,
végigvezetve
a
tervezés
lépéseit
egy
mintapéldán. Nem ez a többváltozós rendszerek szabályozásának egyetlen módszere, számos más módszert ismertet a szakirodalom [2,3].
10.7. ábra. Egy 2x2-es MIMO objektum V kanonikus struktúrával
10.2. MIMO rendszerek belső kereszthatásainak problémái visszacsatolt rendszerek esetén Egy 2x2-es többváltozós rendszer mintapéldáján keresztül [4] bemutatjuk a kereszthatások okozta problémákat, melyek bizonyos esetekben az egyszerű visszacsatolásokat sem engedik meg. A 10.8. ábrán bemutatott MIMO rendszer elemeinek átmeneti függvényei (10.9. ábra) alapján meghatározhatjuk a digitális szabályozási rendszerhez a To mintavételezési időt. Ebben a példában a To meghatározásához a legkisebb időállandó a G11 elemhez tartozik, mely értéke ~10 s, tehát To 1 s lesz.
10.8. ábra. MIMO rendszer SIMULINK felületen 156
MIMO tervezés
Step Response From: In(1)
From: In(2)
2
To: Out(1)
1.5 1
Amplitude
0.5 0 5
To: Out(2)
4 3 2 1 0
0
20
40
60 0
20
40
60
Time (sec)
10.9. ábra. A MIMO rendszer egyes elemeinek átmeneti függvényei A szabályozási rendszer vizsgálatánál megnéztük azt az esetet, amikor csak egy szabályozó kör üzemel. Ebben az esetben a rendszer stabilis viselkedésű, a szabályozási algoritmus hangolható. (A 10.10. ábra időfüggvényei az algoritmus másik ágba áthelyezésével változtak):
1-0.795z-1 1-z-1 r(k)=5
u1(t)
y1(t)
PI algoritmus
To=1 s y1(k),y2(k) u2(t)
y2(t)
u2(100)=0.5
To=1 s 2x2 MIMO rendszer
10.10. ábra. MIMO visszacsatolt rendszer viselkedése egy szabályozó kör esetén 157
MIMO tervezés
Abban az esetben, ha mindkét visszacsatolást megvalósítjuk (10.11. ábra) a kereszthatások miatt a szabályozási körök instabilak lesznek, a rendszer stabilizálása nem lehetséges. Ez kimutatható abban az esetben, ha megvizsgáljuk a MIMO rendszer relatív erősítési tényezőinek tömbjét, azaz az RGA-t
1-0.795z-1 1-z-1 r1(k)=5
y1(t)
PI algoritmus 0.3-0.279z-1 1-z-1
r2(k)=5
u1(t)
To=1 s y1(k),y2(k) u2(t)
y2(t)
PI algoritmus1
To=1 s 2x2 MIMO rendszer
10.11. ábra. MIMO visszacsatolt rendszer viselkedése két szabályozóval kompenzátor nélkül
158
MIMO tervezés
10.3. Az RGA (RelativeGainArray)
A 10.11. ábra alapján felépített visszacsatolások esetén e konkrét szabályozott rendszer nem stabilizálható, mely kimutatható a Bristol módszer [6] alapján. A Bristol módszer segít annak eldöntésében, hogy szükséges-e több bemenetű, több kimenetű szabályozó algoritmus, vagy valamilyen egyéb eszköz (kompenzátor) a rendszer stabilis működéséhez, vagy a szabályozási feladat megoldható SISO rendszereknél alkalmazott algoritmusokkal. A fenti módszer egy un. RGA (RelativeGainArray) mátrix segítségével segít megoldani a problémát. Az RGA mátrix a MIMO rendszer K erősítési tényezőinek mátrixából számítható az: RGA= inv(K’).*K
(10.1)
egyenlet alapján (10.1). A MATLAB újabb verziói ismerik az rga(G) parancsot, de az alábbi programot használva hozzájuthatunk az RGA mátrixhoz. % RGA számítása K=[1 2;4 5] RGA=inv(K').*K% ahol K’ a K mátrix pseudoinverse. A mi konkrét esetünkben az RGA 2x2-es mátrixa a következő lesz: RGA = -1.6667 2.6667 2.6667 -1.6667 Az RGA mátrixból elemeinek vizsgálatából a következő megállapításokhoz juthatunk:
RGA(i=j)=1 A szabályozó körök közt nincs egymásra hatás, i=j esetén
RGA ij=0
Az i-dik bemenet nincs hatással a j-dik kimenetre
159
MIMO tervezés
RGAij=0.5 Nagyfokú kölcsönhatás van;vagy másik szabályozó kör akkora hatással van a j-dik kimenetre, mint az i-dik bemenet
0.5
RGAij<0 A negatív átlós elemek azt jelzik, hogy a visszacsatolt rendszer bármely struktúrában instabil lesz. Ezt az utóbbi meghatározást mintapéldánk is mutatja (10.11. ábra), azaz a visszacsatolt
MIMO rendszer nem lehet stabilis.
10.4. Statikus kompenzátor tervezése MIMO rendszerekhez Egy m bemenetű, m kimenetű rendszer esetében az RGA mátrix elemeinek vizsgálatakor kiderülhet, hogy egyszerű visszacsatolásokkal a szabályozás nem oldható meg, szét kell választanunk a csatornákat, hogy a kereszthatások csökkenjenek, vagy teljesen megszűnjenek. Statikus kompenzátor tervezésére jól használható a Hawkins által javasolt módszer [6], mely frekvencia tartományban segít megoldani a feladatot. Hasonló eredményhez jutunk, ha az állandósult állapotban követelünk meg teljes kompenzációt, szétválasztást. Legyen a statikus kompenzátor struktúrája a MIMO objektummal megegyező (10.12.ábra). A k kompenzátor kij elemei és a MIMO objektum K erősítési tényezői állandók. Az állandósult állapotra (nulla frekvencián) felírható a következő követelmény: az u1*(t) bemenet valamilyen erősítési tényezővel (pld. K11-el) csak az y1(t) kimenetre legyen hatással, az y2(t) kimenetre a tranziensek lecsengése után, tehát stacionárius esetben a kereszthatás szűnjön meg. Ez legyen igaz arra az esetre is, ha csak az u2*(t) bemenet változik, az y1(t) stacionárius esetben ne változzon. Ehhez a következő algebrai egyenletrendszert kell megoldani egy 2x2-es MIMO rendszer esetén (10.2): k11*K11 + k21*K12 = K1 k11*K21 + k21*K22 = 0 160
MIMO tervezés
k22*K22 + k12*K21 = K2 k12 *K11 + k22*K12 = 0,
(10.2)
ahol K1 és K2 a k statikus kompenzátorral szétválasztott csatornák erősítési tényezője. 2x2-es STATIKUS KOMPENZÁTOR
2x2-es MIMO OBJEKTUM K11
1 u1*(t)
k11
u1(t)
1 20s2 +9s+1 G11(s)
1 y1(t) K21
4 80s2 +18s+1 G21(s)
k21
K12
2 18s2 +9s+1 G12(s)
k12
K22 5 2 u2*(t)
k22
110s2 +21s+1
u2(t)
2 y2(t)
G22(s)
10.12. ábra. A statikus kompenzátor és a MIMO objektum kapcsolata A 10.2. egyenlet természetesen kiterjeszthető m x m méretű rendszerre is. Példánkban az legyen K1=1 és K2=1, a MIMO objektum erősítési mátrixa K=[1 2;4 5], így kiszámíthatók a statikus kompenzátor elemeinek értékei az alábbi MATLAB parancsokkal: % A statikus kompenzátor tervezése B=[1 0;0 1] k=B/K A statikus kompenzátor tehát: k= -1.6667 0.6667 1.3333 -0.3333 , mely SIMULINK felületen ez a 10.13.ábrának felel meg.
161
MIMO tervezés
10.13. ábra. A 2x2-es statikus kompenzátor elemeinek értékei és struktúrája
Összecsatolva a megtervezett statikus kompenzátort és a MIMO objektumot (10.14. ábra) megvizsgálhatjuk az egymásrahatást. Az egymásra hatás vizsgálatánál az u1*(t) ugrás függvény bemenőjel amplitúdója 5, az u2*(t) jel a 100-dik másodpercben indul, amplitúdója 3. A 10.15. ábra mutatja be a statikus kompenzátor hatását. Az ábra alapján látható, hogy csak nulla frekvencián tökéletes a kompenzáció és az u1*(t) – y2(t) csatornán a kereszthatás nagy, az u2*(t) – y1(t) csatorna
u1*(t)
u1(t)
u1(t)
y1(t)
u1*(t)=5
To=1 s y1(k), y2(k) u2*(t)
u2(t)
u2(t)
y2(t)
u2*(100)=3
To=1 s Statikus kompenzátor
2x2 MIMO rendszer
10.14. ábra. A statikus kompenzátor illesztése MIMO objektumhoz közti kompenzáció nagyon sikeres. Ezzel már elérhetjük azt, hogy a visszacsatolt MIMO rendszer stabilis legyen, hiszen az u1* – y1 csatorna az u2* bemenettől független lett (10.15. ábra).
162
MIMO tervezés
10.15. ábra. A statikus kompenzátor hatása nyitott körben A statikus kompenzátorral visszacsatolt MIMO rendszer (10.16. ábra) viselkedése azt mutatja, hogy a rendszer stabilis, a statikus kompenzátor ezt biztosítja, de a szabályozás minőségét, az egymásra hatás nagyon befolyásolja (10.17. ábra). Ez a zárt körben létrejövő kereszthatás javítható, megváltoztatható, ha a statikus kompenzátor tervezésénél a K1 és K2 paramétereket más értékekre választjuk, illetve ha a digitális szabályozási algoritmus paramétereit megváltoztatjuk.
1-0.798z-1 1-z-1 r1(k)=5
u1(t)
u1(t)
y1(t)
To=1 s
PI1 algoritmus 0.5-0.4503z-1 1-z-1
r2(100)=5
u1*(t)
y1(k), y2(k) u2*(t)
u2(t)
u2(t)
y2(t)
To=1 s
PI2 algoritmus Statikus kompenzátor
2x2 MIMO rendszer
10.16. ábra. Statikus kompenzátorral kiegészített visszacsatolt MIMO rendszer
163
MIMO tervezés
10.17. ábra. Visszacsatolt MIMO rendszer viselkedése statikus kompenzátorral
10.5. Dinamikus kompenzátor tervezése MIMO rendszerekhez
Amint láthattuk a 10.4. fejezetben a statikus kompenzátorok a csatornák egymásra hatását csökkenthetik, de az egyes csatornák kereszthatásai nagyok lehetnek, a technológiai rendszer működését ez károsan érintheti, a szabályozás minősége egyes csatornák esetében nem lesz elfogadható. Amint említettük ez annak a következménye, hogy a kompenzáció nullánál nagyobb frekvenciákon nem működik. Az előrecsatolás tárgyalásánál bemutattuk, hogy dinamikus előrecsatolást alkalmazva a zavarójel csatorna hatását elvileg teljesen kompenzálni tudjuk, ha egyéb feltételek (pl. holtidők) is megfelelnek. A 10.18. ábra egy MIMO rendszer y2 kimeneti csatornájának ’zavarójel’ kompenzációját szemlélteti, ahol a ’zavarójel’ az u1 bemeneti jel. Természetesen a Ge1 ’előrecsatolás’ diszkrét algoritmus.
164
MIMO tervezés
2x2-es DISZKRÉT MIMO OBJEKTUM 0.02155z+0.01855 z2-1.598z+0.6376
u1(k) u1=5
Gd11(z) ELÕRECSATOLÁS t
0.0232z+0.02153 z2-1.787z+0.7985
v
Gd21(z)
Ge1(z)
Scope
0.04714z+0.0399 z2-1.563z+0.6065 Gd12(z)
u2(k)
0.02139z+0.02014 z2-1.825z+0.8338 Gd22(z)
u2(100)=4
10.18. ábra. MIMO rendszer előrecsatolás alkalmazásával Ezt alkalmazva a 10.19. ábrán látható, hogy az u1 hatása az y2 csatornára nem számottevő.
10.19. ábra. Az egy csatornán szétcsatolt MIMO rendszer viselkedése
Felírhatjuk tehát a 2x2-es MIMO rendszer számára mindkét dinamikus kompenzátor impulzus átviteli függvényét a következő egyenletekből kiindulva (10.3): Gd21(z) = - Ge1(z)*Gd22(z) Gd12(z)= - Ge2(z)*Gd11(z)
(10.3)
A dinamikus kompenzátorok impulzus átviteli függvényei a következők lesznek (10.4) [7]: Ge1(z) = - Gd21(z) / Gd22(z)
165
MIMO tervezés
Ge2(z)= - Gd12(z)/Gd11(z)
(10.4)
A MATLAB programmal megtervezett Ge1 ill. Ge2 dinamikus kompenzátorok számlálóinak és nevezőinek együtthatói az r, s, ill. a t, v vektorokban helyezkednek el (10.20. ábra).
u1*(t)
u1(t)
u1(t)
y1(t)
u1*(t)=5
To=1 s y1(k), y2(k) u2*(t)
u2(t)
u2(t)
y2(t)
u2*(100)=4
To=1 s 2x2 MIMO rendszer
Dinamikus kompenzátor
10.20. ábra. A dinamikus kompenzátor struktúrája
Megvizsgálva a fenti nyitott rendszert látható, hogy dinamikus kompenzátor esetén a 2x2-es MIMO rendszer két SISO rendszerrel jellemezhető, tehát a visszacsatolások tervezése leegyszerűsödik SISO rendszerek tervezésére.
10.21. ábra. A dinamikus kompenzátor hatása 166
MIMO tervezés
A 10.22. ábrán bemutatott digitális dinamikus kompenzátorral kiegészített digitális szabályozási rendszer tervezésénél a PI1 és PI2 algoritmusok helyett természetesen egyéb ismert digitális algoritmusok is alkalmazhatók, de tervezésnél a MIMO objektumot mindig ki kell egészíteni a kompenzátorral. Ez a statikus kompenzátorok alkalmazása esetén is igaz.
2-1.798z-1 1-z-1
u1*(t)
u1(t)
u1(t)
y1(t)
To=1 s
PI1 algoritmus
r1(k)=5
0.2-0.186z-1 1-z-1
r2(100)=5
PI2 algoritmus
y1(k), y2(k) u2*(t)
u2(t)
u2(t)
y2(t)
To=1 s Dinamikus kompenzátor
2x2 MIMO rendszer
10.22. ábra. Digitális szabályozási rendszer digitális dinamikus kompenzátorral
A
10.23.
ábrán
a
tranziens
folyamatokat
láthatjuk.
Digitális
dinamikus
kompenzátorokat alkalmazva a mintavételezés ill. a számábrázolás pontossága, az identifikáció bizonyos pontatlansága miatt nem kapunk teljes szétválasztást.
10.23. ábra. A zárt körök viselkedése
167
MIMO tervezés
10.6. Digitális MIMO PI szabályozási algoritmusok Egy SISO digitális PI algoritmus felírható a következő differencia egyenletrendszerrel (10.5):
x u (k 1) a x u (k ) b e(k ) (10.5)
u (k ) c x u (k ) d e(k ), ahol: e(k) – a hibajel, u(k) – a szabályozó kimenő jele, xu(k) – „segéd” változó és a=1, b=1, c=To/Ti (integrálási tényező), d= Kc (erősítési tényező), felépítése pedig az 10.24. ábrán látható.
10.24. ábra. SISO PI algoritmus modellje Adott paraméterek mellet (Kc=2, To/Ti=0.5, To=1) a fenti PI algoritmus átmeneti függvényét mutatja a 10.25. ábra. u(k)
Idő, s 10.25. ábra. A PI algoritmus átmeneti függvénye 168
MIMO tervezés
A digitális MIMO PI algoritmusok általános állapotegyenlete a következő (10.6): x u (k 1) A c x u (k ) B c e(k ) u(k ) C c x u (k ) D c e(k ),
(10.6)
ahol : xu(k)
a szabályozó állapotváltozóinak vektora,
e(k)
a hibajel vektora,
u(k)
a végrehajtó jel vektora,
és az Ac és Bc egységmátrixok, Dc tartalmazza a Kc arányossági tényezőket, Cc pedig a To/Ti állandók értékeit. Két bemenetű és kimenetű rendszerek esetén ez a következő (10.7):
Dc
kc11 kc 21 kc12 kc 22 (10.7)
Cc
cc11 cc 21 cc12 cc 22
A szabályozó struktúráját mutatja az 10.26. ábra.
10.26. ábra. Digitális MIMO PI szabályozó struktúrája
169
MIMO tervezés
Egységugrás hibajelek esetén egy diszkrét MIMO PI algoritmus átmeneti függvénye például a következő (lásd a majdani mintapéldát) (10.27. ábra):
10.27. ábra. MIMO PI szabályozó átmeneti függvénye Példa Legyen egy konkrét 2x2-es folytonos dinamikus rendszer, melynek átviteli mátrixa:
G 0 ( s)
G11( s) G21( s) G12 ( s) G22 ( s)
és SIMULINK programja az 10. 28. ábra [4]:
10.28. ábra. 2x2 –es MIMO rendszer 170
MIMO tervezés
A zárt szabályozó kör hangolását a MATLAB NCD BlockSet segítségével végeztük el. A MIMO diszkrét PI algoritmussal visszacsatolt fenti rendszer viselkedését To=1 s mintavételezési idő mellett néhány ábrával demonstráljuk. A megtervezett diszkrét szabályozó mátrixai a következők voltak: |
|
|
|
|
|
|
|
A szabályozó rendszer SIMULINK modellje az NCD felhasználásával (10.29. ábra):
10.29. ábra. MIMO PI algoritmus tervezése NCD segítségével
A megtervezett MIMO PI algoritmussal visszacsatolt rendszer viselkedését mutatja a 10.30. ábra. Az r1(k) alapjelváltás a 20., az r2(k) pedig a 120. időegységnél változott. A folytonos rendszer kimenetei: y1(k) és y2(k).
171
MIMO tervezés
10.30. ábra. Szabályozás MIMO PI algoritmussal Függelék 1: MATLAB program MIMO rendszerek szabályozásának tervezéséhez: % MIMO rendszerek tervezése % Pannon Egyetem, Dr.Vass József % A 2x2 -es vizsgált rendszer elemei [4] G11=tf(1,[20 9 1]) G21=tf(4,[80 18 1]) G12=tf(2,[18 9 1]) G22=tf(5,[110 20 1]) % A vizsgált rendszer átviteli mátrixa G=[G11 G12;G21 G22] % RGA számítása K=[1 2;4 5] RGA=inv(K').*K % A statikus kompenzátor tervezése B=[1 0;0 1] %K1=1, K2=1
172
MIMO tervezés
k=B/K % Dinamikus digitális kompenzátor tervezése % Diszkrét modell előállítása: G(s) >> G(z); To= 1 s Gd11=c2d(G11,1) Gd21=c2d(G21,1) Gd12=c2d(G12,1) Gd22=c2d(G22,1) % Dinamikus kompenzátor tervezése Ge1=-Gd12/Gd11 Ge2=-Gd21/Gd22 % A dinamikus kompenzátor vektoros formája [r,s]=tfdata(Ge1,'v') [t,v]=tfdata(Ge2,'v')
Irodalomjegyzék [1]
Housam Binous: Control of the Wood and Berry Distillation Column Using Simulink,2006.(http://www.mathworks.com/matlabcentral/fileexchange/12132control-of-the-wood-and-berry-distillation-column-using-simulink)
[2]
K. Dutton,S. Thompson:The Art of ControlEngineering,Addison-Wesley, 1997.
[3]
B. C. Kuo: Digital Control Systems, Sounders College Publishing, 1992.
[4]
Nagy L., Vass J.: MIMO szabályozási rendszerek tervezése. Mérés és Automatika 36.évf.,1988. 1. szám
[5]
Bristol, E. H., “On a newmeasure of interactionformultivariableprocesscontrol,” IEEE TransactionsonAutomaticControl, vol. 11, no. 1, pp. 133-134.
[6]
Hawkins, D. J.:Pseudodiagonalisation and inverseNyquistarray. Proc. IEE 116,103 (1969)
[7]
R. Isermann: Digital Control Systems, Springer-Verlag, 1989. 173
Diszkrét Állapotszabályozók Tervezése
11. DISZKRÉT ÁLLAPOTSZABÁLYOZÓK TERVEZÉSE
11.1. Diszkrét állapotmodellek MATLAB/SIMULINK környezetben Egy szakirodalomból vett tesztfeladaton keresztül bemutatjuk be az általunk felhasznált struktúrát, a megfigyelhetőségi kanonikus formát. A dinamikus objektum harmadrendű, folytonos modellje legyen [1]: 4 s
Go( s)
y( s) (1 2 s)e u( s) (1 10 s)(1 7 s)(1 3s) ,
impulzus átviteli függvénye To=4 s-os mintavételezési időnél, d=1 diszkrét holtidővel : T0=input('Mintavételezési idő (s)=') % A harmadrendű d=1 diszkrét holtidővel rendelkező SISO objektum % paraméterei[Isermann R. Digital Control Systems, Springer-Verlag,1989]: b1=.06525;b2=0.04793;b3=-.00750 a1=-1.49863;a2=.70409;a3=-.09978 Go=tf([b1 b2 b3],[1 a1 a2 a3 0],4) A mintavételezett objektum Go(z) impulzus-átviteli függvénye tehát: Transferfunction>>Go:
A lineáris diszkrét állapotegyenlet általános alakja: xk 1 A xk B uk y k C xk
(11.1)
A 11.1. diszkrét állapotegyenlet A,B,C,D mátrixainak felépítése a felhasználás céljaitól függően különbözők lehetnek [1]. Lehet diagonális kanonikus forma, vagy egyéb kanonikus struktúrák. Kiválasztásuk attól függ, milyen tervezési feladatot szeretnénk megoldani.
174
Diszkrét Állapotszabályozók Tervezése
Ha az általunk választott megfigyelhetőségi kanonikus struktúrát választjuk, akkor a harmadrendű holtidős dinamikus rendszerünk diszkrét állapotmodelljének mátrixai/vektorai az (11.2) struktúrával rendelkeznek: x1 k 1 0 x k 1 1 2 x3 k 1 0 x4 k 1 0 vagy
b3 x1 k 1 0 b2 x2 k 1 0 u k , b1 x3 k 1 0 0 x4 k 1 1
0 a3 0 a2 1
a1
0
0
x d k 1 A d x d k b d uk , x1 k x k y k 0 0 1 0 2 , x3 k x4 k vagy y k c Td x d k .
(11.2)
(Megj.: A 11.2. állapotegyenletben a d index holtidős rendszerre utal) A holtidő az állapotmodellben a 11.1. ábra alapján értelmezhető. (Az ábrán a B*u, A*u és C*u jelölések MATLAB sajátosságok)
u(k)
1
u(k-1)
z
1
1
z
z
x(k+1)
u(k-d) B* u
x(k)
1
y(k) C* u
z B
Unit Delay
C
A* u A
11.1. ábra. Bemeneti holtidővel rendelkező diszkrét állapotmodell struktúrája
A SIMULINK könyvtára tartalmazza a (11.1) állapotegyenletnek megfelelő felületet is.
A
mátrixok
megadhatók
MATLAB
paraméterfelületen.
175
felületen,
vagy
vektoros
formában
a
Diszkrét Állapotszabályozók Tervezése
y(n)=Cx(n)+Du(n) x(n+1)=Ax(n)+Bu(n) u(k)= 5
Scope Discrete State-Space
11.2. ábra. Diszkrét állapotmodell SIMULINK felületen A 11.2 ábra alapján programozott modellnek csak az y(k) kimenetéhez férhetünk hozzá, az egyes állapotváltozók viselkedését nem tanulmányozhatjuk. A 11.3. ábrán bemutatott modell azonban megengedi a hozzáférést az összes állapotváltozóhoz.
x(k) 1
B* u u(k)= 5
x(k) C* u
z B
Unit Delay
C
y(k)
A* u A
11.3. ábra. Az x(k) állapotváltozók hozzáférését biztosító modell. 176
Diszkrét Állapotszabályozók Tervezése
A modellt felépíthetjük egy olyan hatásvázlattal is, amiből könnyen következtethetünk az objektum differencia egyenlet rendszerére (11.4. ábra) [1]. Ez a differencia egyenletrendszer struktúra az állapotváltozókhoz a hozzáférést szintén biztosítja. 1 u(k)= 5
x4(k)
z b2
b3 b3
b1
x1(k+1)
z
To=4 s 1
1
1
-a3
b1
b2
z
x1(k)
-a2 -a3
x2(k)
z
x3(k)
y(k)
-a1 -a2
-a1
11.4. ábra. Az állapotmodell hatásvázlata A differencia egyenletrendszer a következő lesz (11.3):
x1 k a 3 x 3 k 1 b3 x 4 k 1
x 2 k x1 k 1 a 2 x 3 k 1 b2 x 4 k 1 x 3 k x 2 k 1 a1 x 3 k 1 b1 x 4 k 1 x 4 k uk 1
(11.3) y k x3 k
A d=1 diszkrét holtidő miatt a harmadrendű rendszer négy állapotváltozóval rendelkezik. Általánosan kijelenthető, hogy a diszkrét lineáris rendszer állapotváltozóinak száma m+d, ahol az m a rendszer fokszáma.
177
Diszkrét Állapotszabályozók Tervezése
11.2. Diszkrét állapotszabályozók A lineáris diszkrét dinamikus rendszer állapotegyenlete (11.4): x( k + 1 ) = A x( k ) + B u( k )
(11.4)
ahol az A, B konstans paramétermátrixok és a kezdeti feltétel x(0) adottak. Tegyük fel továbbá, hogy az összes állapotváltozó x(k) tökéletesen mérhető (11.5. ábra). x(0) u(k)
x(k+1)
Iz-1
B
x(k)
y(k)
C
A
11.5. ábra. Lineáris diszkrét állapotmodell Ekkor egy olyan szabályzót kell tervezni, ami az x(k) állapotváltozó vektorból olyan u(k) állapotváltozó vektort generál, ami az egész rendszert x(N)=0 végállapotba juttatja, és a négyzetes teljesítmény kritérium minimális lesz (11.5). N 1
I x T ( N ) Sx( N ) [ x T (k )Qx(k ) uT (k ) Ru(k )] k 0
(11.5)
ahol: S
pozitív szemidefinit, szimmetrikus mátrix
Q
pozitívszemidefinit, szimmetrikus mátrix
R
pozitivdefinit, szimmetrikus mátrix
azaz xTSx 0, xTQx és uTRu> 0. Meg kell jegyezni azt is, hogy a referencia változók és a külső zavar hatásától eltekintünk, és a kimeneti változó (11.6): y( k ) = C x( k )
(11.6)
nincs visszacsatolva. Tervezzük meg a K állapotszabályozót megadva a Q és R mátrixokat, vagyis oldjuk meg az ismert Riccati mátrix differencia egyenletet [1] a MATLAB segítségével. A Q mátrix a mintapélda alapján egy 4x4-es mátrix és csak a szabályozott változót súlyozzuk (11.7). 178
Diszkrét Állapotszabályozók Tervezése
0 0 Q 0 0
0 0 0 0 0 0 0 1 0 0 0 0
(11.7)
Mivel a vizsgált rendszer egy bemenetű – egy kimenetű így a (11.5) célfüggvényben az R mátrix skalár lesz. A mi esetünkben Az u(k) végrehajtójel súlyozásalegyen r=0.043 [1]. A Riccati egyenlet megoldása a következő program segítségével megadja az állapotszabályozó vektorát: % Az u(k) végrehajtójel súlyozási faktora a célfüggvényben r=input('r=') % A célfüggvény Q mátrixa Q=[0 000;0 0 00;0 0 1 0;0 0 00]; % A digitális állapotszabályozó K mátrixának számítása (SISO esetben ez % vektor) [K,S,e]=dlqr(A,B,Q,r) mely a következő lesz: K = [ 4.8191 5.0162 4.4578 0.5332] Az optimális állapotszabályozóval visszacsatolt rendszer hatásvázlata (11.6. ábra) alapján:
11.6. ábra. Lineáris állapotmodell K állapotszabályozóval [1] felépíthetjük a mintapéldánkban szereplő rendszert SIMULINK felületen (11.7.ábra).
179
Diszkrét Állapotszabályozók Tervezése
x(k) y(k)
1
B* u
B
C* u
z
x(k)
Unit Delay
C
A* u u(k)
A K K
Demux
K = [ 4.8191 5.0162 4.4578 0.5332]
11.7. ábra. Lineáris állapotmodell K állapotszabályozóval SIMULINK felületen Abban az esetben, ha az x(k) állapotváltozóknak egységnyi kezdeti értéket adunk megfigyelhetjük a visszacsatolt rendszer viselkedését (11.8.ábra).
11.8. ábra. Az állapotváltozók és a végrehajtó jel viselkedése a zárt rendszerben Az állapotváltozók viselkedésén látható, hogy a zárt rendszer komplex gyökpárral rendelkezik. Ezt a [K,S,e]=dlqr(A,B,Q,r) programsor hatására láthatjuk a futtatási eredményben is: e=
0.3660 + 0.3087i 180
Diszkrét Állapotszabályozók Tervezése
0.3660 - 0.3087i 0.2335 -0.0000 A Q és Rmátrixokkal tudjuk hangolni a szabályozó kört, de ezen mátrixok elemeinek megválasztására nincs szabályrendszer, több iterációt kell elvégezni.
11.3. Diszkrét állapot megfigyelők Az állapotszabályozó feltételezi, hogy az irányítandó folyamat összes állapotváltozója ismert, mérhető. Ez a feltétel az esetek legnagyobb részében nem így van, legtöbbször méréssel csak az y(k) kimenőjelekhez és az u(k) bemenőjelekhez férünk hozzá. Legyen a diszkrét dinamikus rendszer állapotegyenlete (11.8): xk 1 A xk B uk y k C xk
(11.8)
Kapcsoljunk a folyamathoz párhuzamosan egy ugyanilyen modellt a 11.9. ábrán bemutatott struktúrában és képezzük a hibajelet (11.9):
e k y k y k
(11.9)
Ezt egy H mátrix-xal súlyozzuk és visszacsatoljuk az x k 1 állapotváltozón keresztül.
11.9. ábra. A dinamikus rendszer és a megfigyelő kapcsolata [1] Ezt a megfigyelő struktúrát Lauenberger állapot-megfigyelőnek nevezzük. A H konstans értékű visszacsatoló mátrixot olyannak kell választanunk, hogy az ̂(k+1) 181
Diszkrét Állapotszabályozók Tervezése
aszimptotikusan közelítse az x(k+1)-et k esetén. A 11.9. ábra szerinti struktúra a következő állapotegyenlethez vezet (11.10):
x k 1 A x k B u k H e k
A x k B u k H y k C x k
(11.10)
A tényleges és a becsült állapotváltozók közti hiba (11.11, 11.12):
x~ k 1 x k 1 x k 1
(11.11)
x~ k 1 A HC x~ k
(11.12)
vagy
A becslő pólusait a H mátrix megfelelő megválasztásával tudjuk befolyásolni. A megfigyelő H mátrixának számításához meg kell változtatni a tervezési egyenletek mátrixait [1] a következőképpen:
Ao AT ;
Bo CT ;
K HT ;
illetve a [H,So,eo]=dlqr(Ao,Bo,Qo,ro) parancshoz elő kell állítani a Qo és ro mátrixokat:
1 0 Qo 0 0 0
0 0 0 1 0 0 0 1 0 0 0 1 0 0 0
0 0 0 ; ro 5 0 25
A [H,So,eo]=dlqr(Ao,Bo,Qo,ro) parancs a megfigyelő paramétereit szolgáltatja:
h T 0,061 0,418 0,984 1,217 1,217 és megadja a megfigyelő pólusait: eo =
0.2814 + 0.3290i 0.2814 - 0.3290i 0.2136 0.7383 0.0000
182
Diszkrét Állapotszabályozók Tervezése
A dinamikus rendszer és az állapot megfigyelő SIMULINK felületen a következőképpen programozható (11.10. ábra):
1 In1
x(k)
1
B* u
1 Out1
C* u
z B
C
Unit Delay
2 Out2
A* u A
1 In1
u *H H
2 Out2 2 In2
Bm* u
1
K*u
z B 0
Unit Delay1
1 Out1
C0
K*u AB 01
11.10. ábra. Az diszkrét állapotegyenlet és a megfigyelő SIMULINK felületen A teljes rendszer struktúráját SIMULINK felületen a 11.11. ábra mutatja meg. A dinamikus rendszer állapotváltozóinak kezdeti értéke x(0)=1. A tényleges x(k) és becsült ~x(k)állapotváltozók, valamint a tényleges és becsült kimenetek viselkedése a 11.12. ábrán látható.
183
Diszkrét Állapotszabályozók Tervezése
DINAMIKUS RENDSZER Out1
u(k)
y(k)
In1 Out2
u(k)=0
x(k)
delta e(k) ÁLLAPOT MEGFIGYELÕ In1
Out1
In2
Out2
~y(k) y(k),~y(k)
~x(k) x(k),~x(k)
11.11. ábra. Dinamikus rendszer diszkrét állapot megfigyelővel
11.12. ábra. A valós x(k), y(k)) és a becsült ~x(k), ~y(k) változók
11.4. Diszkrét állapotszabályozó LUENBERGER megfigyelővel Abban az esetben, ha az objektumot külső zavarások érik és w(k) alapjelet kap a zárt szabályozó kör - a hibajelünk (11.13): e (k ) w(k ) y(k ) n(k ) Cν(k ) x(k ) η(k ) Cε(k ).
lesz és megváltozik a rendszer struktúrája a 11.13. ábra szerint [1].
184
(11.13)
Diszkrét Állapotszabályozók Tervezése
11.13. ábra. Az állapotmegfigyelővel kiegészített szabályozó kör struktúrája [1] melynek SIMULINK programja a 11.14. ábra:
11.14. ábra. Az állapotmegfigyelővel kiegészített szabályozó kör SIMULINK felületen Az a differencia egyenletrendszer, melyet minden mintavételezéskor meg kell oldani az alábbi [1]: -
a megfigyelő kimeneti hibájának számítása (11.14):
e k 1 y k 1 x 3 k 1 185
(11.14)
Diszkrét Állapotszabályozók Tervezése
-
állapotváltozók becslése (11.15): xˆ 1 k a 3 xˆ 3 k 1 b3 xˆ 4 k 1 h1 e k
xˆ 2 k xˆ 1 k 1 a 2 xˆ 3 k 1 b2 xˆ 4 k 1 h2 e k xˆ 3 k xˆ 2 k 1 a1 xˆ 3 k 1 b1 xˆ 4 k 1 h e k xˆ 4 k xˆ 5 k 1 uk 1 h4 e k 1 xˆ 5 k xˆ 5 k 1 h5 e k 1
(11.15) -
az u(k) végrehajtójel számítása (11.16):
uk k1 xˆ 1 k k 2 xˆ 2 k k3 xˆ 3 k k 4 xˆ 4 k k5 xˆ 5 k
(11.16)
11.5. Esettanulmány Állapotszabályozó tervezése DC motorhoz A 4. fejezetben tanulmányozott QUANSER DCMCT rendszer rugó-tömeg kiegészítéssel (11.15. ábra) keskenysávú lengéseket végez átmeneti állapotban, amit mérésekkel igazoltunk. Ha ilyen rendszerekhez állapotszabályozót kívánunk tervezni, figyelembe kell venni ezt a tulajdonságot. A diszkrét állapotmodell mintavételezési idejét:
11.15. ábra. A QUANSER DC motor rugó-tömeg kiegészítő mechanikával
186
Diszkrét Állapotszabályozók Tervezése
11.16. ábra. A DC motor dinamikus viselkedése csökkenteni kell. Az 11.17. ábra alapján a lengések periódus ideje ~0.1 s, válasszuk a mintavételezési időt ezen periódusidő huszadára, azaz legyen To=0.005 s.
11.17. ábra. Az identifikált DC motor átmeneti függvénye Így a DC motor modelljét (ld.4. fejezet) B(q) = 2.009 q-2 - 3.575 q-3 + 1.897 q-4 A(q) = 1 - 2.603 q-1 + 2.489 q-2 - 0.8666 q-3
187
Diszkrét Állapotszabályozók Tervezése
a d2d parancs segítségével átalakítjuk To=0.005 s mintavételezési időhöz (ld. a program részletet) % Impulzus-átviteli függvény együtthatói a d2d parancsra b1=1.039;b2=-2.006;b3=1.01; a1=-2.863;a2=2.796;a3=-0.9309 T0=0.005 % DC motor modellje Gdc=tf([b1 b2 b3],[1 a1 a2 a3],T0) % Állapottér modell A=[0 0 -a3 b3;1 0 -a2 b2; 0 1 -a1 b1;0 000] B=[0;0;0;1]; C=[0 0 1 0]; D=0 Az állapotszabályozó tervezéséhez szükséges paramétereket megadva: % Az állapotszabályozó tervezése r=input('r=') % r=0.1 Q=[0 000;0 0 00;0 0 1 0;0 0 00]; megtervezzük a szabályozó K vektorát: [K,S,e]=dlqr(A,B,Q,r), mely a következő lesz: K = [0.9635 2.6420 4.9117 0.8575] A megfigyelő tervezésénél az egyes mátrixok kiegészítése és transzformációk után megkapjuk a megfigyelő paramétereit: %Megfigyelő tervezése: %A mátrixok kiegészítése: Am=[0 0 -a3 b3 0;1 0 -a2 b2 0; 0 1 -a1 b1 0;0 0 00 1;0 000 1] Cm=[0 0 1 0 0] % Transzformációk :A---A' és B--C' Ao=Am' Bo=Cm' 188
Diszkrét Állapotszabályozók Tervezése
Co=[0 0 1 0 0] Do=0 ro=input('ro=') % ro=5 Qo=[1 0 000;0 1 0 00;0 0 1 0 0;0 0 0 1 0;0 0 00 20] [H,So,eo]=dlqr(Ao,Bo,Qo,ro) Ahol a h vektor : h= [1.4327 -2.9820 1.8095 0.5723 0.5723] lesz, majd a 11.12. ábra szerint kiegészítjük a struktúrát a következő formában: Am=[A B;0 1] Bm=[-B;0] Km=[K 1] Cm=[C 0] A DC motor állapotszabályozó-megfigyelő rendszere az 11.18. ábrán látható. r(k)=100 0.09516 z-0.9048 u(k)
SZÛRÕ, To=0.005 DC MOTOR ÁLLAPOT MODELL y(k) In1
e(k)
Out1
u(k)
y(k)
MEGFIGYELÕ ee(k)
In1
Out1
In2
Out2
e^(k) x^(k)
u(k) x^(k) u(k) Km* u KI
K = 0.9635 2.6420 4.9117 0.8575 1.000 H = 1.4327 -2.9820 1.8095 0.5723 0.5723
189
Diszkrét Állapotszabályozók Tervezése
11.18. ábra. A DC motor szabályozásának SIMULINK felülete és tranziensei
190
Diszkrét Állapotszabályozók Tervezése
Függelék 1. Állapotszabályozók tervezése %Állapot szabályozók tervezése (Tesztfeladat [1]) %Pannon Egyetem, Dr.Vass József T0=input('Mintavételezési idő (s)=') % A harmadrendű d=1 diszkrét holtidővel rendelkező SISO objektum % paraméterei[Isermann R. Digital Control Systems, Springer,1989]: b1=.06525;b2=0.04793;b3=-.00750 a1=-1.49863;a2=.70409;a3=-.09978 Go=tf([b1 b2 b3],[1 a1 a2 a3 0],4) % Az objektum diszkrét állapottér modelljének mátrixai A=[0 0 -a3 b3;1 0 -a2 b2; 0 1 -a1 b1;0 000] B=[0;0;0;1] C=[0 0 1 0] D=0 Gss=ss(A,B,C,D,T0) [A,B,C,D]=ssdata(Gss) % Az u(k) végrehajtójel súlyozási faktora a célfüggvényben r=input('r=') % A célfüggvény Q mátrixa Q=[0 000;0 0 00;0 0 1 0;0 0 00]; % A digitális állapotszabályozó K mátrixának számítása (SISO esetben ez % vektor) [K,S,e]=dlqr(A,B,Q,r) %Megfigyelő tervezése %A mátrixok kiegészítése: 191
Diszkrét Állapotszabályozók Tervezése
Am=[0 0 -a3 b3 0;1 0 -a2 b2 0; 0 1 -a1 b1 0;0 0 00 1;0 000 1] Cm=[0 0 1 0 0] % A---A' és B--C' Ao=Am' Bo=Cm'; Co=[0 0 1 0 0] Do=0 ro=input('ro=') Qo=[1 0 000;0 1 0 00;0 0 1 0 0;0 0 0 1 0;0 0 00 25] [H,So,eo]=dlqr(Ao,Bo,Qo,ro) Goo=ss(Ao,Bo,Co,Do,T0); O=[0 000]; Am=[A B;O 1] Bm=[-B;0] Km=[K 1] Cm=[C 0]
Függelék 2. Állapotszabályozó-megfigyelő tervezése DC motorhoz % Állapotszabályozó-megfigyelő tervezése DC motorhoz % DC motor modellje % Impulzus-átviteli függvény b1=1.039;b2=-2.006;b3=1.01; a1=-2.863;a2=2.796;a3=-0.9309 T0=0.005 Gdc=tf([b1 b2 b3],[1 a1 a2 a3],T0) 192
Diszkrét Állapotszabályozók Tervezése
% Állapottér modell %[A,B,C,D]=ssdata(Gss) A=[0 0 -a3 b3;1 0 -a2 b2; 0 1 -a1 b1;0 000] B=[0;0;0;1] C=[0 0 1 0] D=0 % Az állapotszabályozó tervezése r=input('r=') % r=0.1 Q=[0 000;0 0 00;0 0 1 0;0 0 00]; [K,S,e]=dlqr(A,B,Q,r) %Megfigyelő tervezése: %A mátrixok kiegészítése: Am=[0 0 -a3 b3 0;1 0 -a2 b2 0; 0 1 -a1 b1 0;0 0 00 1;0 000 1] Cm=[0 0 1 0 0] % A---A' és B--C' Ao=Am' Bo=Cm' Co=[0 0 1 0 0] Do=0 ro=input('ro=') % ro=5 Qo=[1 0 000;0 1 0 00;0 0 1 0 0;0 0 0 1 0;0 0 00 20] [H,So,eo]=dlqr(Ao,Bo,Qo,ro) Goo=ss(Ao,Bo,Co,Do,T0); O=[0 000]; Am=[A B;O 1] Bm=[-B;0] 193
Diszkrét Állapotszabályozók Tervezése
Km=[K 1] Cm=[C 0]
Irodalomjegyzék [1] R. Isermann: Digital Control Systems, Springer-Verlag, 1989. [2] K. Dutton,S.Thompson:The Art of Control Engineering, Addison-Wesley, 1997. [3] K. J. Aström, J. Apkarian: DC Motor Trainer (DCMCT), Instructor Workbook, QUANSER Engineering, Canada, 2008.
194