ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
Titulni strana
1
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
zadání
2
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
3
Vysoké učení technické v Brně
Vysoké učení technické v Brně Fakulta elektrotechniky a komunikačních technologií Ústav automatizace a měřicí techniky
Prediktivní regulátory s principy umělé inteligence v prostředí MATLAB – B&R Diplomová práce
Studijní program:
Kybernetika, automatizace a měření
Student:
Bc. Libor Matys
Vedoucí:
prof. Ing. Petr Pivoňka, CSc.
Abstrakt: Diplomová práce se zabývá problematikou prediktivního řízení zvláště Model (Based) Predictive Control (MBPC nebo MPC). V první části jsou porovnány identifikační metody. Algoritmus rekurzivní metody nejmenších čtverců je porovnán s identifikačními metodami založenými na neuronových sítí. Další části jsou věnovány prediktivnímu řízení. Je zde popsáno vytvoření MPC se sumační složkou a adaptivního MPC. Pevně nastavený PSD regulátor je porovnán s MPC. Jsou porovnány odezvy na poruchu a na změnu dynamiky regulované soustavy. Porovnání je provedeno na simulačních modelech v prostředí MATLAB/Simulink a na fyzikálním modelu připojeného k PLC firmy B&R.
Klíčová slova: Prediktivní řízení, Model Predictive Control, adaptivní MPC, MPC se sumačním členem, PSD regulátor, odezva na poruchu, odezva na změnu dynamiky soustavy, identifikace metodou nejmenších čtverců, učení metodami založenými na neuronových
sítí,
Levenberg-Marquardt.
učení
metodou
back-propagation,
učení
metodou
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
Brno University of Technology The Faculty of Electrical Engineering and Communication Department of Control and Instrumentation
Predictive controllers with principles of artificial intelligence Thesis
Specialisation of study:
Cybernetics, automation and measurement
Student:
Bc. Libor Matys
Supervisor:
prof. Ing. Petr Pivoňka, CSc.
Abstract: Master’s thesis deals with problems of predictive control especially Model (Based) Predictive Control (MBPC or MPC). Identifications methods are compared in the first part. Recursive least mean squares algorithm is compared with identification methods based on neural networks. Next parts deal with predictive control. There is described creation MPC with summing element and adaptive MPC. There is also compared fixed setting PSD controller with MPC. Responses on disturbance and changes of parameters of controlled plant are compared. Comparing is made on simulation models in MATLAB/Simulink and on physical model connected to PLC B&R.
Keywords: Predictive control, Model Predictive Control, adaptive MPC, MPC with summing element, PSD regulator, responses on disturbance, responses on changes parameter of plant, recursive least mean squares identification, identification based on neural networks, back-propagation learning method, Levenberg-Marquardt learning method.
4
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
Bibliografická citace MATYS, Libor. Prediktivní regulátory s principy umělé inteligence v prostředí MATLAB – B&R Diplomová práce. Brno: Vysoké učení technické v Brně, Fakulta elektrotechniky a komunikačních technologií, 2008. s 57., příloh. Vedoucí práce: prof. Ing.Petr Pivoňka, CSc.
5
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
Prohlášení „Prohlašuji, že svou diplomovou práci na téma Prediktivní regulátory s principy umělé intelenge v prostředí MATLAB – B&R jsem vypracoval samostatně pod vedením vedoucího diplomové práce a s použitím odborné literatury a dalších informačních zdrojů, které jsou všechny citovány v práci a uvedeny v seznamu literatury na konci práce. Jako autor uvedené diplomové práce dále prohlašuji, že v souvislosti s vytvořením této diplomové práce jsem neporušil autorská práva třetích osob, zejména jsem nezasáhl nedovoleným způsobem do cizích autorských práv osobnostních a jsem si plně vědom následků porušení ustanovení § 11 a následujících autorského zákona č. 121/2000 Sb., včetně možných trestněprávních důsledků vyplývajících z ustanovení § 152 trestního zákona č. 140/1961 Sb.“
V Brně dne :
Podpis:
Poděkování
Děkuji tímto prof. Ing. Petru Pivoňkovi, CSc. za cenné připomínky a rady při vypracování diplomové práce.
V Brně dne :
Podpis:
6
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
OBSAH 1. ÚVOD ...............................................................................................................12 2. PREDIKTIVNÍ ŘÍZENÍ.................................................................................14 2.1 Historie............................................................................................................14 2.2 Princip .............................................................................................................15 3. IDENTIFIKACE .............................................................................................18 3.1 Identifikace metodou nejmenších čtverců s modelem ARX...........................18 3.2 Identifikace pomocí neuronových sítí.............................................................21 3.2.1 Úvod .............................................................................................................21 3.2.2 Historie .........................................................................................................21 3.2.3 Umělý neuron ...............................................................................................22 3.2.4 Neuronové sítě ..............................................................................................24 3.2.5 Učení neuronových sítí .................................................................................25 3.2.6 Učení metodou back – propagation ..............................................................27 3.2.7 Učení metodou Levenberg – Marquardt.......................................................27 3.3 Realizace .........................................................................................................28 4. MPC ZALOŽENÝ NA VNĚJŠÍM POPISU .................................................31 4.1 Další části prediktivních regulátorů ................................................................31 4.1.1 Účelová (objektivní) funkce .........................................................................31 4.1.2 Referenční trajektorie ...................................................................................33 4.1.3 Omezení........................................................................................................34 4.1.4 Optimalizace .................................................................................................35 4.2 Prediktivní regulátor (MPC) s omezeními......................................................36 4.3 Adaptivní MPC ...............................................................................................37 5. MPC ZALOŽENÝ NA STAVOVÉM POPISU............................................40 5.1 Odvození algoritmu ........................................................................................40 5.2 Pseudostavový popis .......................................................................................42 5.3 MPC s integračním charakterem.....................................................................44 5.3.1 Sumační člen v algoritmu MPC....................................................................44 5.3.2 Sumační člen na výstupu MPC.....................................................................46 6. POROVNÁNÍ MPC S PSD REGULÁTOREM ...........................................50
7
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
6.1 Porovnání MPC s PSD regulátorem při skokové změně poruchy .................50 6.2 Porovnání MPC s PSD regulátorem při změně dynamiky soustavy...............52 6.3 Vizualizační prostředí .....................................................................................56 7. POUŽITÍ MPC PRO ŘÍZENÍ FYZIKÁLNÍCH MODELŮ.......................57 7.1 Skoková změna poruchy na vstupu soustavy..................................................58 7.2 Změna dynamiky soustavy .............................................................................60 8. ZÁVĚR .............................................................................................................63 9. LITERATURA ................................................................................................64
8
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
SEZNAM SYMBOLŮ ZKRATEK Symboly: A, B, C , D
matice stavového popisu
b
vektor pro řešení optimalizační funkce
E
jedničkový sloupcový vektor
ER
energie
f
vektor hodnot volné odezvy
h
vstup neuronu
F ( s)
spojitý tvar přenosu
F
matice dynamiky stavů
F ( z)
diskrétní tvar přenosu
G
první parciální derivace modelu podle proměnných
H
druhá parciální derivace modelu podle proměnných
I
jednotková matice
J
kriteriální funkce
Jc
Jacobiho matice
k
krok
K
zesílení MPC
Kr
suma prvků prvního řádku matice S
KR
zesílení regulátoru
Kx
první řádek matice Q
M
horizont řízení
N
filtrační konstanta
N1 , N 2
minimum a maximum horizontu predikce
P
horizont predikce
Pk
kovarianční matice
Q
matice zesílení stavů
r
referenční trajektorie
S
matice zesílení referenční trajektorie
9
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
TD
derivační časová konstanta
TI
integrační časová konstanta
Tt
filtrační konstanta
u
akční zásah
∆u
změna akčního zásahu
t
čas
Ts
perioda vzorkování
v
měřená porucha
V
vektor chyb
w
váha neuronu
x
stav
X
trénovací množina, množina vzorů
y ) y
hodnota výstupu
Z
dolní trojúhelníková matice
α
penalizace akčního zásahu
β
penalizace odchylky
δ
faktor L-M metody
ϕ
vektor minulých hodnot vstupů a výstupů
ε
chyba
Φ
matice tvořená z vektorů ϕ
η
momentum
λ
koeficient zapomínání
µ
konstanta učení
Θ
práh neuronu
θ
vektor parametrů soustavy
ξ
vstupní potenciál neuronu
predikovaná hodnota výstupu
10
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
Zkratky: ART
Adaptive Resonance Theory
ARX
AutoRegressive with eXternal input
B-P
back-propagation
CRHPC
Constrained Receding Horizont Predictive Control
DMC
Dynamic Matrix Control
EHAC
Extended Horizont Adaptive Control
EPSAC
Extended Prediction Self-Adaptive Control
GMV
Generalized Minimum Variance control
GPC
Generalized Predictive Control
L-M
Levenberg-Marquardt
MAC
Model-based Algorithmic Control
M(B)PC
Model (Based) Predictive Control
MPHC
Model Predictive Heuristic Control
MURHAC
Multipredictor Receding Horizont Adaptive Control
MUS-MAR
Multistep Multivariable Adaptive Control
NN
Neural Network
PFC
Predictive Functional Control
PID
proporcionálně integračně derivační regulátor
PLC
Program Logic Controller
PSD
proporcionálně sumačně derivační regulátor
QDMC
Quadratic Dynamic Matrix Control
RHTC
Receding Horizont Tracking Control
RLS
rekurzivní metoda nejmenších čtverců
SIORHC
Sable Input/Output Receding Horizont Control
SGPC
Stable GPC
UI
umělá inteligence
UPC
Unified Predictive Control
XOR
eXclusive OR
11
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
1.
ÚVOD
Pojem Model (Based) Predictive Control (MBPC nebo MPC) se objevuje v druhé polovině sedmdesátých let 20.století. Model Predictive Control není určený pro jednu specifickou strategii řízení, ale představuje poměrně široký rozsah řídicích metod, které využívají při výpočtu optimálního akčního zásahu odhad budoucí hodnoty výstupu regulované soustavy. Různé MPC algoritmy se mezi sebou liší jen v modelu použitého k reprezentování řízené soustavy, šumu a v použité účelové funkci pro optimalizaci. V současné době již existuje mnoho aplikací pro úspěšné použití prediktivního řízení a to nejen ve zpracovatelském průmyslu, ale například i v sušicích věžích a pro řízení robotických ramen [1] a v mnoha dalších aplikací. Úspěšné použití v těchto aplikacích ukazuje, jak může MPC dosáhnout vysoké účinnosti řízení systémů. Mezi výhody prediktivního řízení se může zahrnout schopnost řízení velmi různorodých procesů. Od systémů s poměrně jednoduchou dynamikou až po ty s více složitou, včetně systémů s dlouhým časovým zpožděním, systémů s neminimální fází nebo nestabilních systémů. Pokud známe referenční trajektorii (známe dopředu průběh žádané hodnoty), lze na změnu žádané hodnoty reagovat s předstihem a tím i dosáhnout lepší sledování žádané hodnoty. Je samozřejmé, že prediktivní řízení má i jisté nevýhody a úskalí. Je nutný vyšší výpočetní výkon než jaký vyžadují klasické PID (PSD) regulátory. Ovšem se soudobým stavem již tento problém přestává být aktuální. Asi největší problém představuje nutnost znát adekvátní model řízeného procesu. V praxi MPC prokázal být přijatelným řešením pro průmyslové řízení, navzdory původnímu nedostatku teoretického řešení v některých klíčových bodech jako je stabilita a robustnost [1]. Cílem této práce je vytvoření MPC a jeho ověření na fyzikálním modelu soustavy. Nejdříve bude nutné nastudovat metodiku návrhu prediktivních regulátorů. A protože MPC ke své správné činnosti potřebuje znát model řízené soustavy, je potřeba se seznámit s klasickými identifikačními metodami a metodami založenými na neuronových sítí. Jako klasická identifikační metoda bude uvažována metoda nejmenších čtverců.
12
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
Po seznámení s identifikačními metodami bude přistoupeno k vytvoření jednotlivých identifikačních algoritmů v prostředí MATLAB. A následně budou vlastnosti identifikačních algoritmů porovnány. Pak již bude možné vytvořit algoritmus MPC. Protože se jedná o řízení v otevřené smyčce, bude třeba se zabývat zaručením nulové ustálené odchylky. K tomu přirozeně vybízí vytvoření adaptivního MPC, který parametry modelu soustavy upravuje i během regulace. Předpokladem je, že ani adaptivní MPC nebude tím pravým řešením a bude nutné hledat i jiné možnosti. Následovat bude porovnání MPC s pevně nastaveným PSD regulátorem zvláště při působení poruchových signálů a při změně dynamiky soustavy. Bude vytvořeno vizualizační prostředí, které bude umožňovat porovnání MPC s PSD regulátorem. Vizualizační prostředí je určeno hlavně pro potřeby výuky. Aby bylo možné komunikovat s fyzikálním modelem, bude třeba využít pro
řízení soustavy programovatelného automatu. Na fyzikálních modelech bude ověřena schopnost regulace MPC při působení poruch a při změně dynamiky fyzikálního modelu.
13
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
2.
PREDIKTIVNÍ ŘÍZENÍ
2.1
HISTORIE
V roce 1978 J. Richalet, A. Rault, J. L. Testud a J. Papon publikovali Model Predictive Heuristic Control: Application to Industrial Processes (MPHC) [2]. Pro určení modelu používali impulsní odezvu. V roce 1980 uveřejnili C. R. Cutler a D. L. Ramakter algoritmus využívající pro identifikaci procesu odezvu na jednotkový skok. Svoji metodu nazvali Dynamic Matrix Control (DMC) [3]. Tyto algoritmy představují první generaci MPC technologie. Dokázali se prosadit i do aplikací v průmyslovém řízení. Nicméně tyto metody neumožňují optimální průběh regulace, protože neberou v úvahu množství omezujících podmínek, které se objevují v praktickém nasazení. MPC se rychle stává populárním, zvláště v chemickém průmyslu, díky jednoduchosti algoritmu a použitím impulsní či přechodové charakteristiky. Druhou generaci MPC algoritmů představuje QDMC (Quadratic Dynamic Matrix Control) s níž přišel v roce 1983 C. R. Cutler [4]. Zde je úloha řešena numericky jako úloha kvadratického programování (QP) s omezujícími podmínkami ve tvaru lineárních nerovností. Jedna v současnosti z nejvíce populárních metod je metoda Generalized Predictive Control (GPC) [5]. Byla publikována v roce 1987, jejími autory jsou D. W. Clarke, C. Mohtadi, P.S. Tuffs. Tato metoda vychází z adaptivní koncepce
řízení a je použitelná na systémy s neminimální fází, nestabilní systémy v otevřené smyčce, systémy s proměnným nebo neznámým zpožděním a systémy neznámého
řádu. Tato metoda se řadí do třetí generace MPC algoritmů. Z GPC vychází další metody, například RHTC (Receding Horizont Tracking Control) [6] lišící se v použitém modelu soustavy (stavový model) a v tom, že tato metoda je rekurzivní [7]. V dnešní době již existuje nespočet metod MPC, které během vývoje vznikly. Jmenujme Model-based Algorhitmic Control – MAC (Richalet), Prediction Self-Adaptive Control (Peterka), Extended Horizont Adaptive Control – EHAC
14
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
(Ydstie), Extended Prediction Self-Adaptive Control – EPSAC (De Keyser), Generalized Minimum Variance control – GMV (Clarke, Gawthrop), Predictive Functional Control – PFC (Richalet), Constrained Receding Horizont Predictive Control - CRHPC (Clarke, Scattolini), Sable Input/Output Receding Horizont Control – SIORHC (Mosca), Stable GPC – SGPC (Kouvaritakis), Multistep Multivariable Adaptive Control – MUS-MAR (Grenco), Multipredictor Receding Horizon Adaptive Control – MURHAC (Lemos, Mosca), Unified Predictive Control - UPC (Soeterboek) [1, 7] a mnohé další.
2.2
PRINCIP
Metodika všech regulátorů patřících mezi MPC může být popsána takto (obrázek 2.1):
u (t+k|t)
u(t)
y (t+k|t) y (t) P
t-1
t
t+1
…
t+k
…
t+P
Obrázek 2.1: Princip činnosti MPC [1]
1. Budoucí hodnoty výstupu pro určitý horizont predikce P jsou počítány v každém ) diskrétním časovém okamžiku t. Tyto predikované hodnoty výstupu y (t + k | t ) pro k = 1…P závisí na známých hodnotách až do času t (minulé hodnoty vstupů a
15
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
výstupů) a na budoucím řídicím signálu (akční zásah) u (t + k | t ) , k = 1…P-1. ) ) Zápis y (t + k | t ) můžeme číst, hodnota výstupu y v čase t + k se vypočítá v čase t. 2. Hodnoty akčních zásahů jsou počítány tak, aby byl průběh optimální podle zvoleného kritéria a výstup co nejlépe sledoval referenční trajektorii r (t + k | t ) . Kritérium optimalizace bývá obvykle kvadratické a vyhodnocuje se z rozdílu predikovaného výstupního signálu a predikované referenční trajektorie. Objektivní (účelová) funkce definuje kritérium optimalizace. 3. Pro regulaci je použita hodnota akčního zásahu u (t | t ) , protože v dalším kroku již známe y (t + 1) , s touto hodnotou se opakuje celý předchozí postup. Je získáno
u (t + 1 | t + 1) (které se od u (t + 1 | t ) liší, protože výpočet u (t + 1 | t + 1) je získán z nových informací) za použití ustupujícího horizontu.
účelová funkce
referenční trajektorie
predikovaná odchylka
predikce výstupu model
omezení
optimalizace predikované hodnoty akčního zásahu
minulé hodnoty výstupu a akčního zásahu
Obrázek 2.2: Základní schéma prediktivního regulátoru (MPC) [7] Pro názornost si lze prediktivní řízení představit na řízení automobilu, kdy
řidič zná referenční trajektorii pro konečný horizont predikce. Pomocí plynového pedálu, brzd a zatáčení se snaží sledovat požadovanou trajektorii. Řidič jede kupředu, vidí před sebe a je schopen dopředu reagovat na situaci. Zatímco při použití klasického řízení (například PID regulátoru) je akční zásah prováděn až na základě znalosti minulé chyby. Jakoby tedy řidič pro řízení využíval jen zpětné zrcátko. Tato přirovnání není úplně přesné, ale pro názornost je plně dostačující.
16
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
Struktura prediktivního regulátoru je zobrazena na obrázku 2.2. Blok model má rozhodující roli, protože vybraný typ modelu musí být schopen bezpečně a dostatečně přesně aproximovat dynamiku řízeného procesu. Proto je velmi důležité správně volit typ modelu a identifikační metodu. Jednou z nejvíce populárních metod je model impulsní charakteristiky. Pro určení modelu stačí znát výstupní hodnoty z procesu při buzení impulsem. Tato metoda je značně oblíbena v průmyslovém použití, protože je velmi intuitivní a může být použita pro vícerozměrné procesy. Nevýhodou je vysoký počet potřebných parametrů a možnost popsat jen systémy stabilní v otevřené smyčce. Těsně související s tímto modelem je model přechodové charakteristiky. Model přenosové funkce je nejvíce rozšířený na akademické úrovni a je využíván pro nejvíce navržených metod řízení. Tato reprezentace modelu vyžaduje jen několik parametrů a je možné ji využít na jakýkoliv proces. Pro vícerozměrné procesy lze také využít model soustavy ve stavovém prostoru [1]. V určitých aplikacích je možné se setkat také s modelem využívající neuronové sítě.
17
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
3.
IDENTIFIKACE
Pro kvalitní prediktivní řízení je nezbytně nutný kvalitní model procesu. Bez kvalitního modelu není možné správně predikovat výstupní hodnoty ŷ a tím se samozřejmě zhorší i celková kvalita regulace. Cílem modelu je přiblížit se k chování procesu. K tomuto účelu slouží různé metody identifikace. Obecně můžeme model získat matematickou analýzou pochodů v procesu nebo analýzou naměřených dat. V případě modelů pro adaptivní řízení a tím i pro prediktivní řízení výrazně dominuje analýza naměřených dat. V adaptivním řízení se v převážné míře používá regresní model (ARX) s metodou nejmenších čtverců, i když v současné době se začínají dosti prosazovat algoritmy založené na bázi neuronových sítí. Identifikace pro adaptivní řízení vychází z následujících podmínek:
•
Vstupy do procesu jsou generovány regulátorem.
•
Regulátor by měl kompenzovat poruchy a stabilizovat proces. Přítomnost poruch zhoršuje identifikaci parametrů soustavy a parametry regulátoru mohou být pak špatně určeny.
•
Struktura identifikovaného modelu se obvykle nemění.
•
Identifikační algoritmus musí být numericky spolehlivý a dostatečně rychlý.
•
Na počátku identifikace má poměrně důležitou roli počáteční nastavení parametrů regulátoru (zabránění nekorektním akčním zásahům). Odhady parametrů by měly již na začátku identifikace dostatečně reprezentovat proces, nebo regulátor má na počátku identifikace parametry přednastaveny tak, aby byl regulační děj přijatelný (s přetlumenou dynamikou, nestačí pouze stabilní děj).
3.1
IDENTIFIKACE METODOU NEJMENŠÍCH ČTVERCŮ S MODELEM ARX
Použití modelu ARX vychází z jeho vhodnosti použití pro průběžnou identifikaci. Parametry modelu ARX se určují metodou nejmenších čtverců. Výhodou této metody je, že může být použita pro libovolné počáteční podmínky i
18
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
19
Vysoké učení technické v Brně
působící poruchy. Předpokladem je jeden vstup a výstup systému, pak závislost mezi vstupem a výstupem je
y = br l
(3.1)
V případě j měření a při rozdílu mezi měřenými hodnotami můžeme určit optimální velikost parametru br . Kritérium je pak dáno j 1 j ( yim − yiv )2 = 1 ∑ ( yim − br li )2 ∑ 2 i =1 2 i =1
J (br ) =
(3.2)
kde yim je měřená hodnota a y iv je vypočítaná hodnota. Hledané minimum je j
br =
∑l y i
i =1
im
(3.3)
j
∑l i =1
2 i
Po rozšíření na n vstupů či vnitřních stavů se lze model popsat experimentální regresní rovnicí
y (i ) = ϕ1 (i )θ1 + ϕ 2 (i )θ 2 + ϕ 3 (i )θ 3 + ... + ϕ n (i )θ n = ϕ T (i )θ + ε
θ = (θ1 θ 2 θ 3 ... θ n )
T
(3.4)
ϕ T (i ) = (ϕ1 (i ) ϕ 2 (i ) ϕ 3 (i ) ... ϕ n (i )) kde y je odhad výstupní veličiny modelu, θ je hledaný vektor neznámých parametrů,
ϕ je vektor známých měřených funkcí, i je krok výpočtu, ε je chyba v kroku výpočtu. Proměnné φi označujeme jako regresní proměnné a model je tedy nazýván regresním modelem. Po n měřeních lze učinit první odhad parametrů modelu. Za předpokladu, že parametry vektoru θ (ty je snahou identifikovat) se během trvání identifikace nemění, po n-tém kroku lze určit jejich hodnotu. Převedeno do maticového tvaru
y = Φθ + ε
⇒
ε = y − Φθ
(3.5)
Hledané minimum účelové funkce
J (Φ ) =
1 T 1 ε ε = ( y − Φθ )T ( y − Φθ ) 2 2
(3.6)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
20
Vysoké učení technické v Brně
jejíž minimum lze získat, jestliže derivace podle vektoru parametrů θ bude rovna 0. Po vyřešení platí
(
θ = Φ TΦ
(
kde Φ T Φ
)
−1
)
−1
ΦT y
(3.7)
se označuje jako Pk (i ) , což je kovarianční matice. Matice Φ by se bez
redukce po každém kroku rozšiřovala o další řádek. Tomu lze zabránit vytvořením rekurentního vzorce s konstantním rozměrem matice Φ . S využitím rekurentního vzorce lze odvodit průběžná metoda nejmenších čtverců. Z (3.7) pak pro krok i platí θ (i ) = Pk (i )Φ T (i ) y
(3.8)
S použitím lemmy o inverzi matice, je pak po úpravách určeno v kroku i + 1
[
]
Pk (i + 1) = Pk (i ) − Pk (i ) (i + 1) 1 + ϕ T (i + 1)Pk (i ) (i + 1) ϕ T (i + 1)Pk (i ) −1
(3.9)
Odvození vztahu pro odhad parametrů podle [8]
θ (i + 1) = Pk (i + 1)[Φ T (i ) y (i ) + ϕ (i + 1) y (i + 1)]
(3.10)
θ (i + 1) = θ (i ) + Pk (i )ϕ (i + 1)[1 + ϕ T (i + 1)Pk (i )ϕ (i + 1)] ( y (i + 1) − ϕ T (i + 1)θ (i )) −1
(3.11)
Pro sledování pomalých změn parametrů identifikovaného procesu je možné použít exponenciální zapomínání. Minimalizací modifikovaného kritéria se vztahy (3.9, 3.11) upraví na výsledné rovnice doplněné o faktor exponenciálního zapomínání λ
(0 < λ
2
)
≤1 .
[
]
K (i + 1) = Pk (i )ϕ (i + 1) λ2 I + ϕ T (i + 1)Pk (i )ϕ (i + 1)
(
)
Pk (i + 1) = Pk (i ) − K (i + 1)ϕ T (i + 1)Pk (i ) / λ 2
(
−1
(3.12) (3.13)
)
θ (i + 1) = θ (i ) + K (i + 1) y (i + 1) − ϕ (i + 1)θ (i ) T
(3.14)
Vztah (3.14) odhaduje parametry pro následující krok. K je váhový součinitel, který určuje s jakým vlivem se má tento rozdíl vzít v úvahu a tím urychlit či zpozdit aktualizaci parametrů. Hodnoty v kovarianční matici Pk se volí co největší, obvykle Pk = 10 4 I .
postačí θ = (a1
a2
... a m
přenosu (3.15).
b1
Výsledkem
b2
je
pak
vektor
... bn ) , jednotlivé prvky odpovídají koeficientům
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
21
Vysoké učení technické v Brně
b1 z −1 + b2 z −2 + ... + bn z − n F (z ) = 1 + a1 z −1 + a 2 z − 2 + ... + a m z − m
3.2
(3.15)
IDENTIFIKACE POMOCÍ NEURONOVÝCH SÍTÍ
3.2.1 Úvod Pojem umělá inteligence je velmi těsně spjata s přirozenou inteligencí jakožto s vlastností živého organismu. Po zjištění, že mozek je tvořen spoustou neuronů navzájem propojených a že tyto neurony spolu komunikují pomocí elektrických impulsů, které lze měřit, začal Američan McCulloch pracovat na sestavení modelu biologického neuronu. Na počátku padesátých let se tak objevují reálné úvahy o sestrojení stroje obdařeného inteligencí a vzniká nový vědní obor zvaný umělá inteligence (UI). Od umělé inteligence a tedy i od neuronových sítí je požadováno, aby uměly uložit znalosti (vědomosti), aplikovat znalosti pro řešení daného problému (uvažování) a získat v průběhu experimentu nové znalosti (učení). Výzkum umělé inteligence se s vývojem počítačové techniky rozdělil do dvou základních směrů. Syntetický přístup, u kterého počítače pracují se symboly, které reprezentují objekty reálného světa. Dále zastánci tohoto přístupu tvrdili, že existence symbolů je nutnou a postačující podmínkou pro existenci inteligentní činnosti. Výsledkem tohoto přístupu jsou tzv. expertní systémy, z kterých vzniká a rozvíjí se Fuzzy logika. Druhým směrem vývoje je analytický přístup navazující na práci McCullocha. Zde je snaha v počítačích věrně modelovat funkce mozku, a právě tento přístup dal vzniknout umělým neuronovým sítím [9].
3.2.2 Historie V roce 1943 Američané W. S. McCulloch a jeho student W. Pitts vytvořili matematický model neuronu. Číselné hodnoty parametrů takového modelu byly bipolární, tzn. z množiny {-1; 0; 1}. Takový model tvoří základ většiny umělých neuronových sítí dodnes. Rosenblatt F. v roce 1958 vytvořil první funkční perceptronovou síť. Perceptron je zobecnění McCullochova a Pittsova modelu pro reálný číselný obor parametrů. Tato nová perceptronová síť byla ovšem schopná řešit
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
pouze problémy, které byly lineárně separabilní. Na tuto skutečnost upozornili M. Minsky a S. Papert. Úspěšně se jim podařilo spolu s dalšími nešťastnými událostmi zdiskreditovat výzkum umělých neuronových sítí ve prospěch jiných oblastí výzkumu umělé inteligence. Pro svou argumentaci využili faktu, že jeden perceptron nemůže realizovat jednoduchou logickou funkci XOR. Tento argument je správný, ale již tři perceptrony jsou schopny takovou funkci realizovat. Nicméně v té době nebyl znám algoritmus pro učení vícevrstvých perceptronů. Obecně se předpokládalo, že ani takový algoritmus vzhledem ke své složitosti není ani možný. K renesanci dochází v osmdesátých letech minulého století. V roce 1986 se nezávisle na sobě podařilo D. Rumelhartovi a LeCunovi odvodit algoritmus pro učení vícevrstvých neuronových sítí, který je dosud nejpoužívanější metodou učení neuronových sítí. Jedná se o tzv. Back - Propagation Algorithm. Vznikají také samoorganizující sítě (např. Kohenenova,), které ke svému učení nepotřebují učitele [9]. V osmdesátých letech vznikají i další typy sítí jako Hopfieldova síť, Grossbergova ART síť atd. [10].
3.2.3 Umělý neuron Umělý neuron je zobrazen na obrázku 3.3 se vstupy h1, h2, h3, …,hn, váhami spojení ) w1, w2, w3, …,wn a výstupem y . Agregace vstupních signálů, jejich porovnání s prahovou hodnotou Θ a následně jejich nelineární zobrazení, představují model těla neuronu. Umělý neuron je tedy charakterizován vnitřním potenciálem, prahem a aktivitou svého výstupu.
22
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
23
Vysoké učení technické v Brně
vstupní signály
synaptické váhy
h1
w1
h2
w2
h3
w3
hn
wn
práh neuronu h0 Θ w0
aktivační funkce
výstup neuronu
y f(ξ)
Obrázek 3.3: Schéma umělého neuronu [8] Vstupy do neuronu mohou být obecně jak číselné hodnoty, tak i nečíselné parametry (např. obrazce). Většinou se vyskytují právě číselné údaje. V závislosti na poloze neuronu mohou být vstupy výstupy jiných neuronů, popřípadě vstupy reprezentují podněty z vnějšího okolí. Synaptické váhy ovlivňují jednotlivé vstupy do neuronů a tím i celou neuronovou síť. Hodnota synaptické váhy reprezentuje citlivost, s jakou příslušný vstup ovlivňuje výstup z neuronu. Váhy jsou většinou dány jistým reálným číslem, právě změna hodnot vah a jejich postupné ladění představuje podstatnou část učicích algoritmů neuronových sítí. Podle biologického neuronu musí vstupní signál překonat prahovou hodnotu , aby se mohl dále šířit sítí. Hodnota prahu tedy určuje, kdy je neuron aktivní a neaktivní. Pokud je hodnota vstupního signálu nižší, je na výstupu neuronu signál odpovídající pasivnímu stavu neuronu. Po překročení prahové hodnoty je na výstupu signál daný agregační funkcí (3.16). Agregační funkce neuronu slučuje jistým způsobem vstupní signály neuronu, v našem případě (obrázek 3.3) funkci agregace představuje sumace. Pak pro vstupní potenciál platí n
n
i =1
i =0
ξ = ∑ hi wi + Θ = ∑ hi wi
(3.16)
Pro zjednodušení může být zvolen práh Θ = h0 w0 = 1 ⋅ w0 . Účelem aktivační funkce je převést hodnotu vstupního potenciálu na výstupní hodnotu z neuronu. V principu je možné tyto funkce rozdělit na lineární a nelineární
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
24
Vysoké učení technické v Brně
nebo také na spojité či nespojité. Obrázek 3.4 ukazuje nejčastěji používané aktivační funkce. V případě identické výstupní funkce je výstup neuronu stejný jako výstup z aktivační funkce, v případě neidentické výstupní funkce je výsledek z aktivační funkce dán nějakou funkční závislostí [9].
f(ξ)
f(ξ)
1
1
k
0
ξ
0
f(ξ)
a) lineární funkce f(ξ) = k ξ; k > 0
Θ
0
ξ
b) binární funkce (perceptron) pro ξ > Θ, f(ξ) = 1 pro ξ < Θ, f(ξ) = 0 f(ξ)
f(ξ)
c) binární funkce f(ξ) = sign (ξ)
f(ξ)
1
1
k 0
0 -1
ξ 0
d) lineární funkce s omezením
ξ -1
ξ
e) sigmoida f(ξ) =1/(1+e-kξ); k > 0
ξ
f) posunutá sigmoida f(ξ) =2/(1+e-kξ)-1; k > 0
Obrázek 3.4: Aktivační funkce neuronu [8] 3.2.4 Neuronové sítě Neuronová síť (Neural Network - NN) slouží jako univerzální funkční aproximátor. To znamená, že navržená a podle určitých pravidel vytrénovaná síť je schopná s určitou mírou přesnosti simulovat výchozí proces. Obecně může mít neuronová síť strukturu, která je popsána libovolným orientovaným grafem s vrcholy a s orientovanými hranami [11]. Vrcholy takového grafu tvoří neurony a orientované hrany reprezentují propojení neuronů s jistými vahami. Výstup z jednoho neuronu bývá vstupem do několika dalších neuronů. Společnou vlastností neuronových sítí je
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
jejich vrstvená struktura. Vrstvy dělíme na vstupní vrstvu (vstup signálu z okolí), na vrstvy skryté (všechny vrstvy ležící mezi vstupní a výstupní vrstvou) a vrstvu výstupní (předává signál do okolí). Podle toku signálu lze neuronové sítě značit jako dopředné (signál se šíří jedním směrem) a rekurentní (mezi neurony resp. vrstvami existují zpětné vazby).
h1
y1
hn
ym
Obrázek 3.5: Vícevrstvá dopředná neuronová síť [9] Pro návrh topologie neuronové sítě bohužel neexistuje jednoznačný návod. V různých literaturách je možné narazit na různá doporučení, která jsou ale vždy určena pro řešení určité úlohy a nejsou tedy obecně platná. Ze zkušeností vyplývá, že v řadě případů postačuje pro aproximaci řízeného procesu dvouvrstvá neuronová síť [9]. Počet neuronů ve vstupní a výstupní vrstvě je dán počtem vstupů a výstupů sítě. Počet neuronů ve skryté vrstvě má být tak velký, aby byla síť schopna aproximovat daný průběh s námi požadovanou přesností [8].
3.2.5 Učení neuronových sítí Základní vlastností neuronové sítě je její schopnost učení. To znamená, že síť je schopna reagovat a učit se na změněné podmínky. Ve většině případů se algoritmus učení soustřeďuje na adaptaci synaptických vah. Nastavitelnými parametry také mohou být strmosti aktivačních funkcí nebo změna struktury neuronové sítě. V oblasti řízení se pro učení dopředné neuronové sítě nejčastěji používají algoritmy Levenberg - Marquardt (L-M) a back – propagation (B-P). Při aproximaci přenosové funkce je dostačující síť s jednou skrytou vrstvou. Pokud je aktivační funkce ve tvaru sigmoidy nebo hyperbolického tangensu, je pak taková síť
25
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
schopna aproximovat i nelineární funkci. Ve vstupní vrstvě je použita lineární aktivační funkce a ve výstupní vrstvě je posunutá sigmoida nebo lineární funkce. Jestliže je neuronový model procesu naučen ještě před jeho použitím, mluvíme o učení off-line. Neuronový model se učí z dříve získaných dat bez připojení k procesu a podle algoritmu učení (B-P, L-M) se upraví synaptické váhy sítě (nebo se změní topologie apod.). Takto naučenou síť již může být použita pro identifikaci. Má-li být regulátor adaptivní, pokračuje se v učení i během řízení procesu. Takové učení je nazýváno on-line. Při on-line učení je model připojen k systému a učí se chovat jako daný systém. Nastavení vah modelu se tedy v každém kroku modifikuje podle současného stavu. Z toho vyplývá, že stav vah modelu se určuje ze současného stavu a předchozí stavy se postupně zapomínají. Následkem toho model zapomíná předchozí stavy a pokud se nemění dynamika systému, model se naučí na působení poruchových veličin. Ve výsledku pak nastavení vah neodpovídá systému [8]. Na obrázku 3.6 je znázorněno učení regresního modelu třetího řádu ze vstupních a výstupních hodnot procesu.
u(t)
y(t)
proces u(t-1)
z-1
u(t-2)
z-1
z-1 u(t-3) z y(t-3) y(t-2)
neuronová síť
ŷ(t)
-1
z-1 z-1
y(t-1)
Obrázek 3.6: Učení neuronové sítě se zpožděnými vstupy a výstupy [8]
26
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
27
Vysoké učení technické v Brně
3.2.6 Učení metodou back – propagation Je jedním z nejčastěji používaným algoritmem pro učení neuronových sítí. Cílem algoritmu je nalezení takového nastavení neuronové sítě, aby byla minimalizována chybová funkce
ER =
1 n ( y i − d i )2 ∑ 2 r =1
(3.17)
kde ER je energie, n je počet neuronů výstupu sítě, yi je i-tý výstup, di je i-tý požadovaný výstup. Učící algoritmus je dán vztahem
∆w z (i, j ) = − µ
∂E R ∂wk (i, j )
(3.18)
kde wz (i, j ) je váha spoje mezi neuronem i v jedné vrstvě a neuronem j v následující vrstvě, k určuje, mezi kterými vrstvami je váha počítána (z = 1 – váha mezi vstupní a vnitřní vrstvou, z = 2 – váha mezi vnitřní vrstvou a vrstvou následující), ∆w z (i, j ) je aktuální změna váhy, µ je konstanta učení. Pro naše potřeby lze algoritmus přepsat do podoby [12]
(
)
θ (k + 1) = θ (k ) + η (θ (k ) − θ (k − 1)) + µ y (k + 1) − ϕ T (k + 1)θ (k ) ϕ (k )
(3.19)
Pokud zvolíme větší hodnotu parametru µ bude krok tak velký, že malá lokální minima překročíme. Jestliže bude parametr µ ale moc malý, budou se váhy adaptovat velmi pomalu. Dalším parametrem je η , který se označuje jako momentum. Ten dodává do postupu váhovým prostorem setrvačnost a umožňuje překonávat lokální minima.
3.2.7 Učení metodou Levenberg – Marquardt Jedná se iterativní algoritmus, který numerický hledá minimum sumy čtverců obecně nelineárních funkcích. V podstatě lze každý reálný systém považovat za nelineární, protože obsahuje jistá omezení například nelinearitu typu saturace, vliv A/D, D/A převodníků. Identifikace pomocí L-M algoritmu pracuje na principu hledání globálního minima chyb mezi minulými výstupy a výstupy modelu z trénovací množiny (množina vzorů).
X (k ) = [ϕ (k ) ϕ (k − 1) L ϕ (k − p )]
(3.20)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
28
Vysoké učení technické v Brně
kde p je délka vektoru X. Minimalizaci algoritmus opakuje i počtu iterací v každém kroku k identifikace.
[
θ (i | k + 1) = θ (i | k ) − J c (i | k )J (i | k ) + δ I T
]
−1
J c (i | k )V (i | k ) T
(3.21)
kde V (i | k ) je vektor chyb mezi výstupem z modelu a odhadovaným výstupem ze soustavy T (k ) [12]. Jestliže je zvoleno δ malé, řešení je podobné jako při použití Newtonovy metody. Pokud δ nabývá větších hodnot, je rychlost konvergence parametrů podobná jako u metody back – propagation. V (k ) = T T (k ) − X T (k )θ (k )
(3.22)
T (k ) = [ y (k ) y (k − 1) L y (k − p )]
(3.23)
Jacobiho matice J c (i | k ) představuje nejlepší lineární aproximaci diferencovatelného vektoru získaného v blízkosti daného bodu a vypočtenou v každé iteraci:
J c (k ) =
3.3
(
)
∂ε (k ) ∂ T T (k ) − X T (k )θ (k ) = = − X T (k ) ∂θ (k ) ∂θ (k )
(3.24)
REALIZACE
Jedna z nejdůležitějších částí adaptivního prediktivního regulátoru je model soustavy. Protože prediktivní regulátor postrádá klasickou zpětnou vazbu známou z běžných PID (PSD) regulátorů, představuje model soustavy jedinou informaci o stavu reálného systému. Je proto nezbytně nutné, aby byl model dostatečně věrný. Pro identifikaci je použita rekurzívní metoda nejmenších čtverců (RLS) a dvě metody využívající neuronovou síť, učení metodu back - propagation (B-P) a učení metodou Levenberg – Marquardt (L-M). Neuronová síť byla vytvořena jako jeden lineární neuron se šesti vstupy a váhami. Váhy představují jednotlivé koeficienty modelu. Ve všech následujících případech je použit ARX model soustavy třetího
řádu reprezentovaný vztahem F (z ) =
b1 z −1 + b2 z −2 + b3 z −3 1 + a1 z −1 + a 2 z − 2 + a3 z −3
Pro ověření správnosti algoritmů jsem nechal identifikovat soustavu
(3.25)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
29
Vysoké učení technické v Brně
F1 (s ) =
1
(3.26)
(2s + 1)(s + 1)2
U metody B-P bylo momentum nastaveno na η = 0,001 a parametr učení
µ = 0,0005 . Parametr δ byl u metody L-M nastaven na 0,001. Počet vzorů pro učení neuronových sítí je dobré nastavit tak, aby v nich byla zachycena celá dynamika přechodové charakteristiky aproximované soustavy. Bylo zvoleno 250 vzorů. Počet iterací byl u B-P nastaven na 40 a u L-M na hodnotu 20. Koeficient zapomínání u metody nejmenších čtverců byl nastaven λ = 0,995 . Při identifikaci byla použita perioda vzorkování 0,5 s. Diskrétní ekvivalent soustavy (3.26) je tedy
Fz ( z ) =
0,007654 z −1 + 0,02249 z −2 + 0,004097 z −3 1 − 1,992 z −1 + 1,313z − 2 − 0,2865 z −3
(3.27)
U identifikačních algoritmů hraje dosti významnou roli i počáteční nastavení. U RLS bylo počáteční nastavení algoritmu zvoleno: Pk = 10 6 I a θ = (1 0 L 0) . A pro neuronové sítě je výhodné volit vektor parametrů soustavy θ nenulový. Jako počáteční hodnoty byly zvoleny náhodné hodnoty z intervalu 0; 1 . Z obrázku 3.7 je vidět, že identifikace byla úspěšná a získané modely jsou velmi podobné soustavě (3.26). Zvláště výsledky získané L-M a metodou nejmenších
čtverců jsou velmi podobné. Obecně se dá říci, že metoda nejmenších čtverců konverguje ke správnému výsledku rychleji než metoda B-P. Proto je vhodné, zvláště u metod založených na neuronových sítí, znát alespoň velmi jednoduchý model a podle něj před začátkem identifikace inicializovat algoritmus. Nevýhodou metody učení B-P je její neschopnost v otevřené smyčce správně aproximovat některé soustavy. Například při identifikaci soustavy s astatismem přestaly parametry konvergovat ke správným hodnotám. Pokud ovšem probíhala identifikace v uzavřené smyčce byl algoritmus numericky stabilní. Algoritmus učení L-M zvládá aproximovat i soustavy s astatismem v otevřené smyčce, tato metoda ale vyžaduje větší výpočetní nároky než B-P. Větší výpočetní náročnost je způsobena nutností počítat inverzi matic. Při zkracování periody vzorkování může metoda nejmenších
čtverců selhávat. Pokud se v soustavě nic neděje (nedochází ke změně vstupního signálu do soustavy) kovarianční matice
Pk roste a až nabude určité velikosti,
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
30
Vysoké učení technické v Brně
algoritmus přestane konvergovat. A právě v případech velmi krátké periody vzorkování je výhodnější použít algoritmus identifikace založený na neuronové síti. Prechodove charakteristiky soustavy a modelu
1 y(t) y (t) RLS
0.8
y y
(t)
B-P
(t)
L-M
y
0.6
0.4
0.2
0
0
5
10
15
20
25
t (s)
Obrázek 3.7: Přechodové charakteristiky soustavy F1(s) a identifikovaných modelů
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
31
Vysoké učení technické v Brně
4.
MPC ZALOŽENÝ NA VNĚJŠÍM POPISU
4.1
DALŠÍ ČÁSTI PREDIKTIVNÍCH REGULÁTORŮ
Pro určení predikované hodnoty výstupu procesu se použije známý model procesu. Predikovaná hodnota výstupu se zpravidla určuje ze součtu nucené a volné odezvy modelu. Volná odezva modelu závisí pouze na minulých hodnotách výstupu systému a minulých hodnotách řízení. Předpokladem je konstantní budoucí
řízení - tedy změna akčního zásahu je nulová. Zatímco nucená odezva modelu uvažuje nenulové změny akčního zásahu. y (t ) = y fc (t ) + y fr (t )
(4.1)
kde y(t) je celková odezva, yfr(t) je volná odezva a yfc(t) je nucená odezva modelu.
4.1.1 Účelová (objektivní) funkce
) Hlavním cílem účelové funkce je, aby budoucí výstupy y (t + k | t ) na uvažovaném horizontu predikce co nejlépe sledovaly danou referenční trajektorii
r (t + k ) a také aby změny akčního zásahu ∆u (t + k ) byly penalizovány. Pro objektivní funkci je pak dán výraz
J (N1 , N 2 , M ) =
N2
M
2 2 ∑ β (k )[ y(t + k | t ) − r (t + k )] + ∑ α (k )[∆u(t + k − 1)]
j = N1
)
(4.2)
j =1
V některých přístupech se druhá část výrazu pro penalizaci přírůstku akčního zásahu ) vůbec neuvažuje. Výraz za první sumou y (t + k | t ) − r (t + k ) představuje budoucí hodnoty regulační odchylky v horizontu. Parametry N1 a N2 určují minimum a maximum horizontu predikce P, M je horizont řízení. N2 se volí s ohledem na dynamiku. Vhodné je, aby horizont predikce byl srovnatelný s globální časovou konstantou systému [7]. Pokud je hodnota N1 velká, je to proto, že není důležité, jaká je odchylka v prvních okamžicích. To má za následek hladkou odezvu systému. Jestliže trvá zpracování nějaký čas d, není důvod, aby N1 bylo menší než d, protože výstup se bude odvozovat až od okamžiku t + d. Parametry α (k ) a β (k ) mají konstantní nebo exponenciální závislost. Tyto parametry určují penalizaci odchylky
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
32
Vysoké učení technické v Brně
β (k ) a penalizaci přírůstku akčního zásahu α (k ) . Například trváme-li na exponenciální závislosti β (k ) v průběhu horizontu predikce, můžeme použít vztah
β (k ) = τ N
2 −k
(4.3)
Leží-li hodnota τ na intervalu od 0 do 1, je nejvzdálenější chyba od okamžiku t penalizována více než ty ležící blíže. Důsledkem je hladší odezva s menšími akčními zásahy. Pokud je ovšem τ > 1 jsou první odchylky více penalizovány. Dále je uvažován pouze jeden konstantní parametr α (k ) pro penalizaci přírůstku akčního zásahu, parametr β (k ) = 1 . Tato úprava může být použita, protože minimální hodnota kritéria (2.2) závisí na poměru mezi parametry α (k ) a β (k ) , ne tedy přímo na hodnotách penalizačních parametrů. Předchozí výraz (4.1) a výraz pro výpočet kritéria (4.2) může být přepsán do maticového zápisu, s kterými se bude dále pracovat.
y = Gu + f
(4.4)
J = ( y − r ) ( y − r ) + α (u T u ) T
(4.5)
Po dosazení vztahu (4.4) do (4.5) platí
J = (Gu + f − r ) (Gu + f − r ) + α (u T u ) T
(4.6)
kde f představuje vektor hodnot volné odezvy a součin Gu určuje hodnoty nucené odezvy. Vektor u obsahuje změny akčního zásahu, r je vektor budoucích hodnot referenční trajektorie v horizontu P (horizont predikce). Matice G je Jacobián modelu (parciální derivace modelu podle proměnných) a v případě, že je systém kauzální, je matice G trojúhelníková. V případě lineárního modelu se matice G upraví na
g0 g1 G = g2 M g P
0 g0 g1 M g P −1
0 0 g0 M g P−2
j
j
i =1
i =0
L 0 L 0 L 0 O M L g P −M
,
kde
g j = −∑ ai g j −1 + ∑ bi a g k = 0 ∀k < 0 .
M
je
horizont
řízení
a
kde
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
33
Vysoké učení technické v Brně
Koeficienty ai, bi jsou získány z přenosu modelu soustavy
F (z ) =
b1 z −1 + b2 z −2 + ... + bn z − n 1 + a1 z −1 + a 2 z − 2 + ... + a m z −m
(4.7)
Vztah (4.6) lze upravit do tvaru pro kvadratické programování
J=
1 T u Hu + b T u + f 0 2
(4.8)
kde
(
H = 2 G TG + α I b = 2( f − r ) G
)
T
T
f0 = ( f − r)
T
(f
(4.9)
− r)
Je patrné, že člen f0 není závislý na vektoru u, proto nemá na optimalizační úlohu vliv a může být z rovnice vyloučen. Optimalizační úloha je pak ve tvaru
J=
1 T u Hu + b T u 2
(4.10)
Bez uvažování omezení signálů je možné úlohu zjednodušit na hledání nulového gradientu hodnoty kritéria J. Pak pro vektor u platí
(
u = − H −1b = G T G + α I
)
−1
G T (r − f )
(4.11)
Přírůstek akčního zásahu, který se posílá do systému, je pouze první prvek vektoru u
∆u (t ) = K (r − f )
(4.12)
(
kde K je první řádek matice G T G + α I
)
−1
GT .
4.1.2 Referenční trajektorie Jednou z výhod prediktivního řízení je znalost budoucího průběhu žádané hodnoty, tedy znalost referenční trajektorie, čehož lze náležitě využít. Systém je schopen reagovat na změnu žádané hodnoty dříve, než se skutečně vyskytne na vstupu regulátoru. A může tedy eliminovat vliv zpoždění systému, regulovat na děj bez překmitu, na nejkratší přechodný děj apod. Budoucí vývoj referenční trajektorie je v mnoha případech známý předem, například v robotických systémech, servosystémech. Referenční trajektorie může být zadána jako průběh nebo jako posloupnost nesoucí informaci o tom, v jaký čas a na jakou hodnotu se má referenční trajektorie změnit. Nebo může být zadána jako funkce závislá na výstupu.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
34
Vysoké učení technické v Brně
4.1.3 Omezení Každý proces obsahuje jistá omezení, snímače mají svá omezení, akční členy mají omezený rozsah akčních zásahů, technologická omezení, omezení s ohledem na bezpečnost procesu a uživatelů apod. Všechna taková omezení vedou na požadavek optimálního řízení. Mnoho prediktivních algoritmů počítá s omezeními. Většinou se však jedná o analyticky řešenou optimalizační část bez omezení a až po tomto výpočtu se aplikuje omezení typu saturace. Tento způsob ale nezaručuje optimální
řízení podle kritérií a také je možné jej aplikovat jen na veličiny vystupující z optimalizačního algoritmu (akční zásah a změna akčního zásahu). Pokud se řeší optimalizační úloha s již danými omezeními, je možné omezit nejen veličiny vystupující z algoritmu, ale i výstup soustavy či jednotlivé vnitřní stavy systému (při použití stavového modelu).
Tvrdá omezení Tvrdá omezení jsou taková, jejíž hranice nelze za žádných podmínek porušit. Toto jsou typické příklady tvrdých omezujících podmínek.
u min ≤ u (t ) ≤ u max
∆u min ≤ u (t ) − u (t − 1) ≤ ∆u max y min ≤ y (t ) ≤ y max
pro ∀t ≥ 0
(4.13)
pro ∀t ≥ 0
(4.14)
pro ∀t ≥ 0
(4.15)
kde (4.13) je vztah pro omezení akčního zásahu, další (4.14) pro omezení změny akčního zásahu a (4.15) pro omezení výstupu. Tato omezení lze napsat maticovým zápisem
Ru ≤ c I PM − I PM ZM kde R = P M − ZP M GP − GM P
∆u max E P1 − ∆u min E P1 (u − u (t − 1))E 1 P , c = max (− u min − u (t − 1))E P1 y max E P1 − f − y E1 − f min P
(
(
)
(4.16)
)
kde I je jednotková matice, Z je dolní trojúhelníková matice s nenulovými prvky rovny 1, matice I, Z, G mají P řádků a M sloupců. E je sloupcový vektor o délce P se všemi prvky rovny 1.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
35
Vysoké učení technické v Brně
Měkká omezení Jsou to taková omezení, která jdou za určitých okolností porušit. Ta jsou zavedena z důvodu nevyhovujících tvrdých omezení. A to proto, že v rámci tvrdých omezení nelze nalézt optimální řešení a z technologického a bezpečnostního hlediska je možné určité meze překročit. Jinými slovy tedy měkká omezení umožňují posunutí mezí. K tomu je možné přistupovat různými způsoby [7]:
•
Pevná omezení se průběžně posunují podle předem daných pravidel.
•
Optimalizační část se dělí na dvě části, kde první řeší optimalizaci bez omezení na první části horizontu a druhá řeší optimalizaci na zbylé části horizontu.
•
K pevnému omezení se přidá jistá tolerance posunutí mezí, ale zároveň je velikost posunutí penalizována v kritériu, což je nejuniverzálnější přístup, který je označován jako uvolnění mezí (constrain relaxation).
4.1.4 Optimalizace K získání hodnoty u (t + k | t ) je nezbytné minimalizovat funkční hodnotu J z ) (4.2). Hodnoty predikovaného výstupu y (t + k | t ) jsou počítány jako funkce minulých hodnot vstupů a výstupů a budoucích signálů. K jejich získání se využívá modelu soustavy a hledání minimálních hodnot po dosazení do účelové funkce. Analytické řešení může být získáno z kvadratické účelové funkce (4.11), kdy je model soustavy lineární a nejsou žádné omezující podmínky. Jinak se pro optimalizaci používají iterační metody. Získání řešení jakoukoliv metodou je poměrně obtížné, protože máme N 2 − N 1 + 1 nezávislých proměnných. Při řešení optimalizace při přítomnosti omezujících podmínek je nutné použít numerické optimalizační metody. Ty v případě optimalizace kvadratické účelové funkce (4.10) s omezujícími podmínkami vedou na kvadratické programování [1]. Výsledkem optimalizační funkce (4.2) je přírůstek akčního zásahu. Ten je získán součtem minulé hodnoty akčního zásahu a aktuálního přírůstku.
u (t ) = ∆u (t ) + u (t − 1)
(4.17)
Jestliže je predikovaná regulační odchylka (r-f) nulová, je nulová i změna akčního zásahu. V jiném případě je přírůstek akčního zásahu úměrný k budoucím odchylkám. Zatímco u konvečního zpětnovazebního řízení se akční zásah určuje podle minulé
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
36
Vysoké učení technické v Brně
hodnoty regulační odchylky, u prediktivního regulátoru se určuje přírůstek akčního zásahu z odhadované budoucí regulační odchylky.
Řešení optimalizace kvadratické funkce je poměrně výpočetně náročné, proto je vhodné zvolit dva různé horizonty P (horizont predikce) a M (horizont řízení). Zpravidla se volí P ≥ M . Pokud je horizont řízení kratší než horizont predikce, adekvátně se zmenší i rozměry matic pro následný výpočet kritéria a nalezení řešení není tak výpočetně náročné. Často se pro řešení kvadratické funkce s omezeními používá metody aktivních množin. Pro řešení optimalizačních úloh je možné také použít neuronové sítě.
4.2
PREDIKTIVNÍ REGULÁTOR (MPC) S OMEZENÍMI
Výpočetně nejjednodušší algoritmus optimalizace je takový, v kterém nejsou zahrnuta žádná omezení. Ovšem každý reálný systém má jistá omezení, které musíme při návrhu jakéhokoliv regulátoru zohlednit. Je možné omezení aplikovat až na akční zásah vystupující z regulátoru, tím je ovšem znehodnocen optimalizační algoritmus, protože takto omezené akční zásahy již nejsou optimální. Proto je vhodné zahrnout omezení přímo do optimalizačního algoritmu, pak jsou akční zásahy optimální nejen z hlediska kriteriální funkce, ale i z hlediska navržených omezení. Je nezbytné zdůraznit, že pokud je omezení zahrnuto přímo do algoritmu optimalizace, není možné řešit optimalizaci analyticky jako v případě bez omezení. Úloha s omezeními vede na kvadratické programování. V programu MATLAB je možné problém kvadratické optimalizační úlohy s omezeními řešit pomocí funkce quadprog. Na
obrázku 4.1 jsou porovnány dva prediktivní regulátory, z nichž jeden
s průběhem yB(t) je s omezením − 10 ≤ u ≤ 10 na výstupu a druhý s průběhem yS(t) je se stejným omezením přímo v optimalizačním algoritmu. Pokud nebude uvedeno jinak jsou do všech následujících MPC implementována takto nastavená omezení:
− 10 ≤ u (t ) ≤ 10
− 3 ≤ ∆u (t ) ≤ 3
(4.18)
Omezení je možno rozšířit i na výstupní veličinu, může jít o omezení maximální hodnoty výstupu, popřípadě požadavek na monotónní průběh výstupní veličiny – pak je přechodný děj bez překmitu.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
37
Vysoké učení technické v Brně
Vystup ze soustavy 10 r(t) y B(t)
5
y r,
y S(t)
0 -5 -10 40
50
60
70
80 t (s) Akcni zasah
90
100
110
10 uB(t) u
0
-10 40
uS(t)
50
60
70
80 t (s)
90
100
110
Obrázek 4.1: Porovnání přechodových charakteristik MPC s omezením v optimalizačním algoritmu yS(t), uS(t) a s omezením na výstupu regulátoru yB(t), uB(t) (Ts = 0,5 s; P = 50, M = 25, α = 0,2, řízená soustava F1(s)) 4.3
ADAPTIVNÍ MPC
Výhodou a zároveň i úskalím MPC je snaha vyregulovat predikovanou odchylku danou rozdílem referenční trajektorie a predikovanými hodnotami výstupu. Ne tedy skutečnou (minulou) hodnotu odchylky, a proto nemůže být zaručena nulová ustálená odchylka. Ta lze potlačit adaptivní variantou MPC. V případě pevně nastaveného regulátoru známe model soustavy, který byl identifikován off-line, čili před začátkem vlastního řízení procesu. Během regulace se již parametry modelu nemění. V případě použití adaptivního regulátoru je model procesu získáván on-line,
čili průběžně během regulace. Během každého kroku se tedy parametry modelu procesu upravují. Adaptace modelu vytváří zpětnou vazbu od procesu. Na obrázku
4.2 je porovnání průběhů pevně nastaveného a adaptivního MPC s identifikací metodou nejmenších čtverců. Do t = 170 s probíhá regulace bez problémů, po příchodu skokové změny poruchy na vstup soustavy F1(s) pevně nastavený MPC nezareaguje, protože nemá informaci o přítomnosti poruchy – jedná se o řízení v otevřené smyčce. Adaptivní MPC na poruchu zareaguje a vyreguluje ji. Na všech
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
38
Vysoké učení technické v Brně
obrázcích s průběhy je měřen skutečný akční zásah, který vstupuje do soustavy. Je v něm tedy zahrnuta i porucha vstupující do soustavy. Vystup ze soustavy
4
r(t) v(t) y PN(t)
2
y AD(t)
6
y r,
0 -2 -4 120
140
160
180
200 220 t (s) Akcni zasah
240
260
280
uPN(t)
10
u
300
uAD(t)
0
-10
120
140
160
180
200 t (s)
220
240
260
280
300
Obrázek 4.2: Porovnání průběhů pevně nastaveného MPC yPN(t), uPN(t) s adaptivním MPC yAD(t), uAD(t) s identifikací RLS λ = 0,995 (Ts = 0,2 s; P = 50, M = 25, α = 0,2, řízená soustava F1(s)) Dále byly porovnány tři identifikační algoritmy ve spojení s adaptivním MPC. Dva založené na neuronových sítích a algoritmus založený na metodě nejmenších čtverců. Byl použit ARX model 3. řádu s přenosem podle (3.25) a byla
řízena soustava F1(s). Nastavení jednotlivých algoritmů bylo následující: U rekurzivní metody nejmenších čtverců (RLS) byl koeficient zapomínání λ = 0,995. U neuronové sítě s učením metodou back-propagation (B-P) byla trénovací množina o 150 prvcích, počet iterací byl nastaven na 20, koeficient učení µ = 0,0001 a momentum bylo η = 0,0005 . U metody učení Levenberg-Marquardt (L-M) byl parametr δ = 0,001 , počet iterací byl 20 a velikost trénovací množiny byla nastavena na 250 vzorů. Protože pomocí metody B-P konvergují parametry modelu pomaleji (obrázek 4.3), nebude dále tato metoda používána. U metody RLS a L-M (pro malé
δ ) konvergují parametry soustavy přibližně stejně rychle.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
39
Vysoké učení technické v Brně
Vystup ze soustavy 6
2
r(t) v(t) y RLS(t)
0
y B-P(t)
-2
y L-M(t)
4
y r,
-4 -6
400
450
500 t (s) Akcni zasah
550 uRLS(t)
10
uB-P(t) uL-M(t)
0 u -10
400
450
500
550
t (s)
Obrázek 4.3: Porovnání identifikačních algoritmů RLS (yRLS(t), uRLS(t)), B-P (yB-P(t), uB-P(t)), L-M (yL-M(t), uL-M(t)) při skokovém vzniku poruchy o amplitudě 2 v t = 420 s (Ts = 0,2 s; P = 50, M = 25, α = 0,2, řízená soustava je F1(s)) Aby byla zajištěna nulová ustálená odchylka, muselo by být zesílení modelu naprosto shodné se zesílením soustavy, tento požadavek nemůže adaptivní MPC vždy zaručit. Ovšem adaptivní MPC je vhodný, jestliže se v průběhu regulačního děje mění dynamika řízené soustavy.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
40
Vysoké učení technické v Brně
5.
MPC ZALOŽENÝ NA STAVOVÉM POPISU
5.1
ODVOZENÍ ALGORITMU
Co ale dokáže spolehlivě odstranit ustálenou odchylku, je sumační člen v regulátoru. K tomu, aby bylo možné do MPC přidat sumační člen, je výhodné vytvořit regulátor, který akční zásahy vypočítává na základě stavového popisu soustavy. Pro diskrétní model SISO soustavy
x(k + 1) = Ax(k ) + Bu (k )
(5.1)
y (k ) = Cx(k ) + Du (k )
je hledána řídicí posloupnost u (k ) na horizontu predikce P minimalizující kritérium dané účelovou funkcí
J (N1 , N 2 , M ) =
N2
M
2 2 ∑ β (k )[y(t + k | t ) − r (t + k )] + ∑ α (k )[u (t + k − 1)]
)
j = N1
(5.2)
j =1
Což je stejná účelová funkce jako (4.2) s tím rozdílem, že zde je penalizována velikost akčního zásahu namísto jeho přírůstku. Opět platí
y = Gu + f
(5.3)
J = ( y − r ) ( y − r ) + α (u T u ) T
(5.4)
Po dosazení vztahu (5.3) do (5.4) platí
J = (Gu + f − r ) (Gu + f − r ) + α (u T u )
(5.5)
f = Fx(k )
(5.6)
T
kde
0 0 L C D D 0 L CA CB 2 F = CA , G = CAB CB D L M M O M M CA P −1 CA P − 2 B CA P −3 B L CB
0 0 0 M D
(5.7)
Matice G má P řádků a M sloupců, hodnoty v této matici jsou totožné s maticí G použitou ve vnějším popisu. To souvisí s tím, že tato matice určuje dynamické
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
41
Vysoké učení technické v Brně
vlastnosti systému. V účelové funkci je možné vážit změnu akčního zásahu
∆u (k ) = u (k ) − u (k − 1) , pak je již účelová funkce v typickém tvaru J (N1 , N 2 , M ) =
N2
M
∑ β (k )[ y(t + k | t ) − r (t + k )] + ∑ α (k )[∆u(t + k − 1)] )
j = N1
2
2
(5.8)
j =1
kde ∆u = Di u − u~ ,
1 1 L −1 1 O Di = M O O −1 L −1
1 u (k − 1) M ~ 0 ,u = 1 M 0 1
(5.9)
Z výrazu (5.5) lze odvodit vztah pro optimální posloupnost řízení pro analytický prediktivní regulátor (bez omezení)
u = −Qx(k ) + Sr
(5.10)
kde jednotlivé matice jsou
( S = (G
) G +α I)
Q = G TG + α I T
−1
−1
GTF
GT
(5.11)
Pro ustupující horizont je potřeba spočítat pouze první hodnotu z vektoru u , pak lze prediktivní regulátor vyjádřit ve tvaru
u (k ) = − K x x(k ) + K r r (k )
(5.12)
kde matice v této rovnici jsou P −1
K x = q1T , K r = ∑ s1, j
(5.13)
j =0
kde q1T je první řádek matice Q a s1, j jsou prvky prvního řádku matice S. Odvození bylo provedeno podle [13]. Z (5.12) je vidět, že prediktivní regulátor v analytickém tvaru je podobný stavovému regulátoru. Pro prediktivní regulátor s omezeními je opět nutné řešit problém kvadratického programování ve tvaru
J=
1 T u Hu + b T u 2
pro který platí
(5.14)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
42
Vysoké učení technické v Brně
H = G TG + α I bT = ( f − r) G T
(5.15)
Jediný rozdíl mezi stavovým a vnějším popisem je ve vytváření matic postihující dynamiku řízeného systému. U vnitřního popisu jsou matice získány ze stavového popisu modelu a u vnějšího z vnějšího popisu modelu.
5.2
PSEUDOSTAVOVÝ POPIS
U stavového popisu se získávají matice A, B, C, D ze stavových rovnic. Ne vždy jsou ale k dispozici všechny stavy systému. Pokud některý stav nemůže být přímo měřen, je nutné jej rekonstruovat např. Kalmanovým filtrem. Existuje však další možnost jak pouze za pomocí informací o výstupu systému aplikovat stavový algoritmus MPC. Stavová formulace řízeného systému totiž nemusí být definována v klasickém smyslu. Podmínkou je ale to, že musí být zachována matematická forma stavového popisu. To umožňuje vytvořit nový stavový vektor skládající se ze zpožděných vstupů a výstupů řízené soustavy. Pro systém třetího řádu je tedy nový stavový vektor
u (k − 1) u (k − 2 ) u (k − 3) xio = y (k − 1) y (k − 2 ) y (k − 3)
(5.16)
kde u (k ) je vstup do řízeného systému, tedy akční zásah regulátoru a y (k ) je výstup systému. Pak nové matice stavového modelu systému jsou
xio (k + 1) = Aio xio (k ) + Bio u (k ) y (k ) = C io xio (k )
(5.17)
Matice Dio je uvažována nulová. Pro ARX model platí
y (k ) = b1u (k − 1) + b2 u (k − 2 ) + b3u (k − 3) − a1 y (k − 1) − a 2 y (k − 2 ) − a 3 y (k − 3) Pak stavový model nového systému je definován pomocí následujících matic
(5.18)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
43
Vysoké učení technické v Brně
u (k ) 0 u (k − 1) 1 u (k − 2 ) 0 = y (k ) b1 ( − 1) 0 yk y (k − 2 ) 0
y (k ) = (b1
b2
0
0
0
0
0
0
0
0
1
0
0
0
b2
b3
− a1
− a2
0
0
1
0
0
0
0
1
b3
− a1
− a2
0 u (k − 1) 1 0 u (k − 2 ) 0 0 u (k − 3) 0 + u (k ) − a 3 y (k − 1) 0 0 y (k − 2 ) 0 0 y (k − 3) 0
(5.19)
u (k − 1) u (k − 2 ) u (k − 3) − a 3 ) y (k − 1) y (k − 2 ) y (k − 3)
Výsledkem je tedy takzvaná pseudostavová reprezentace řízeného systému. Pseudostavem je vektor zpožděných vstupů a výstupů soustavy, přičemž je matematická reprezentace systému zachována. Proto můžou být matice A, B, C nahrazeny novými Aio , Bio , C io získanými z ARX modelu [14]. Pseudostav tedy neodpovídá skutečným stavům jako je rychlost, poloha apod. Na obrázku 5.1 je zobrazen průběh MPC získaného ze stavového popisu pro různé hodnoty penalizačního parametru α. Protože je použita účelová funkce s penalizací velikosti akčního zásahu (5.2), neustálí se výstup na žádané hodnotě. Zvláště patrné je to při nastavení α = 0,2 . Pro MPC byl použit off-line získaný model řízené soustavy F2(s)
F2 (s ) =
1
(10s + 1)(s + 1)2
(5.20)
A bylo realizováno pouze omezení pro velikost akčního zásahu − 10 ≤ u ≤ 10 . Na přítomnost poruchy není schopen vytvořený MPC správně zareagovat. Pro vyregulování poruchy je nutné stavový model systému upravit.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
44
Vysoké učení technické v Brně
Vystup ze soustavy 2 1 r(t) v(t) y 0.01(t)
0
y r,
-1
y 0.2(t) -2 100
150
200 t (s) Akcni zasah
250
300
10
0
u
u0.01(t) u0.2(t)
-10 100
150
200 t (s)
250
300
Obrázek 5.1: Průběhy MPC pro různé hodnoty penalizačního parametru α = 0,2 (y0,2(t), u0,2(t)) a pro α = 0,01 (y0,01(t), u0,01(t)) v t = 220 s přichází na vstup soustavy skoková porucha (Ts = 0,2 s; P = 50, M = 25, řízená soustava je F2(s)) 5.3
MPC S INTEGRAČNÍM CHARAKTEREM
5.3.1 Sumační člen v algoritmu MPC Je tedy vytvořen MPC založený na stavovém popisu. Aby bylo možné dosáhnout nulové ustálené odchylky je nutné stavový systém rozšířit o sumační člen, do kterého vstupuje rozdíl mezi referenční trajektorií a skutečným výstupem ze soustavy.
x se (k + 1) = x se (k ) + r (k ) − y (k ) = x se (k ) + r (k ) − Cx(k ) y se (k ) = x se (k )
Stavové rovnice rozšířené o tento sumátor mají následně tvar
(5.21)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
45
Vysoké učení technické v Brně
xio (k + 1) Aio 0 xio (k ) Bio = + u (k ) = A′x ′(k ) + B ′u (k ) x se (k + 1) − C io 1 x se (k ) 0 x (k ) y ′(k ) = (C io − 1) io = C ′x ′(k ) x se (k )
(5.22)
kde y ′(k ) = y (k ) − y se (k ) . S takto upravenými stavovými rovnicemi již nic nebrání ověření funkčnosti algoritmu. Vystup ze soustavy 3 r(t) v(t) y sr(t)
2 1 y r,
0
y br(t)
-1 -2 -3 120
140
160
180 t (s) Akcni zasah
200
220 usr(t)
10
u
ubr(t)
0
-10 120
140
160
180
200
220
t (s)
Obrázek 5.2: Průběhy MPC se sumační složkou v algoritmu se znalostí referenční trajektorie (ysr(t), usr(t)) a bez znalosti referenční trajektorie (ybr(t), ubr(t) v t = 180 s přichází na vstup soustavy skoková porucha (Ts = 0,2 s; P = 50, M = 25, α = 1, řízená soustava je F2(s)) Na obrázku 5.2 jsou zachyceny odezvy výstupu a akčních zásahů při řízení soustavy F2(s) (5.20). Pro MPC byl použit off-line získaný model řízené soustavy F2(s) a bylo použito omezení pro velikost akčního zásahu − 10 ≤ u ≤ 10 . Je zde porovnán MPC se sumační složkou v algoritmu se znalostí budoucího průběhu referenční trajektorie a bez znalosti průběhu referenční trajektorie. Je patrné, že MPC se znalostí referenční trajektorie se snaží reagovat na změnu referenční trajektorie s předstihem, ovšem této akci brání sumační složka regulátoru, do které vstupuje hodnota skutečné odchylky. Sumační složka se pak snaží udržet průběh výstupu tak,
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
aby co nejlépe sledovala aktuální hodnotu referenční trajektorie. Proto se začne výstup co nejrychleji přibližovat k referenční trajektorii až po změně hodnoty referenční trajektorie. A také proto jsou ve výsledku průběhy se znalostí i bez znalosti referenční trajektorie skoro totožné. Pro regulaci je dokonce výhodnější použít MPC, který nevyužívá znalosti referenční trajektorie, protože během přechodného děje nedochází k zakmitnutí. Nicméně pokud je neznámý budoucí průběh referenční trajektorie je možné MPC se sumační složkou použít, protože je schopen zaručit nulovou ustálenou odchylku. To je patrné opět z obrázku 5.2, kde v t = 180 s je na vstup soustavy přivedena skoková porucha. Jinou možností by bylo sumační složku během přechodného děje vypínat. Pokud by ale během přechodného děje působila porucha, byl by i pak průběh přechodného děje zřejmě nepřijatelný.
5.3.2 Sumační člen na výstupu MPC Na docílení nulové ustálené odchylky je možné se dívat i z jiného pohledu. Pro docílení nulové ustálené odchylky stačí, aby sumační (integrační) člen byl přítomen v otevřené smyčce. Při návrhu regulátoru se pak postupuje tak, aby si regulátor myslel, že sumační složka je v soustavě. Toho se dá docílit tím, že do modelu soustavy je přidána sumační složka (do modelu soustavy je přidán astatismus). Regulátor pak počítá akční zásahy, které předpokládají řízenou soustavu s astatismem. Protože strukturu soustavy nelze měnit, přidá se sumační složka na výstup regulátoru. Výsledkem je, že regulátor si myslí, že soustava obsahuje sumační složku, i když v ní ve skutečnosti není, protože je umístěn v regulátoru. Pro realizaci MPC se sumační složkou na výstupu regulátoru byl použit stavový přístup s účelovou funkcí s penalizací velikosti akčního zásahu (5.2). Protože jsou ale akční zásahy počítány z modelu s přidaným astatismem, není ve skutečnosti výstupem optimalizační funkce velikost akčního zásahu, ale jeho přírůstek. Tento přírůstek vstupuje do sumačního členu, na jehož výstupu je už skutečná velikost akčního zásahu. Výhodou je možnost použít omezení na výstupní veličinu a akční zásah podle (4.16). Na obrázku 5.3 je porovnání obou přístupů k přidání sumační složky do regulátoru. MPC jsou nastaveny tak, aby měli stejnou odezvu na poruchový signál.
46
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
47
Vysoké učení technické v Brně
MPC se sumačním členem na výstupu regulátoru dává odezvu s menším překmitem na změnu hodnoty referenční trajektorie. U obou regulátorů je vypnuta znalost budoucích hodnot referenční trajektorie a definováno omezení − 10 ≤ u ≤ 10 . Dále bude používán pouze MPC se sumačním členem na výstupu regulátoru. Na obrázku
5.4 je zobrazen průběh regulace MPC se sumačním členem na výstupu regulátoru při neznalosti a znalosti budoucího průběhu referenční trajektorie. Pro regulaci byla použita soustava s neminimální fází F3 (s ) =
1 − 5s
(2s + 1)3
(5.23)
Při nevhodně zvolené periodě vzorkování např. Ts = 0,2 s nebyl schopen iterační algoritmus optimalizace (s omezeními) určit akční zásahy, které by vedly ke správnému sledování referenční trajektorie. Analyticky řešená optimalizační úloha (bez omezení) vykazovala větší numerickou stabilitu. Po dalším zkracování periody vzorkování i analyticky řešený algoritmus optimalizace přestal konvergovat. Nakonec byla použita perioda vzorkování Ts = 0,5 s a omezení byla implementována podle (4.18). Je patrné, že zvláště MPC se znalostí budoucího průběhu referenční trajektorie je schopen tuto soustavu řídit velmi dobře. Z průběhů akčních zásahů je vidět, že MPC reaguje na změnu referenční trajektorie 25 s předem. To odpovídá součinu PTs = 50 ⋅ 0,5 = 25 s . V případě referenční trajektorie ve tvaru lineárně narůstající funkce je MPC se znalostí budoucího průběhu referenční trajektorie schopen velmi přesně sledovat referenční trajektorii (obrázek 5.5). Chová se, jako by obsahoval alespoň dva sumační členy. Protože z teorie řízení plyne, aby byl výstup schopen přesně sledovat lineárně narůstající žádanou hodnotu, musí být v otevřené smyčce přítomnost alespoň dvou sumačních členů. V MPC je implementováno omezení akčního zásahu − 10 ≤ u ≤ 10 , proto se také snaží MPC v blízkosti tohoto omezení (t = 20 až 22 s) předstihnout průběh referenční trajektorie, aby byl výstup déle schopen sledovat referenční trajektorii. MPC s neznalostí budoucích hodnot referenční trajektorie sleduje referenční trajektorii s určitou konstantní odchylkou. Stejně se chová i PSD regulátor. Při lineárně narůstající poruše (obrázek 5.6) se výstup ustálí s konstantní odchylkou od referenční trajektorie.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
48
Vysoké učení technické v Brně
Vystup ze soustavy 3 r(t) v(t) y ast(t)
2 1 y r,
0
y sum(t)
-1 -2 -3 120
140
160
180 t (s) Akcni zasah
200
220 uast(t)
10
usum(t) 0
u
-10 120
140
160
180
200
220
t (s)
Obrázek 5.3: Průběhy MPC se sumační složkou v algoritmu (ysum(t), usum(t)) a na výstupu regulátoru (yast(t), uast(t) v t = 180 s přichází na vstup soustavy skoková změna poruchy (Ts = 0,2 s; P = 50, M = 25, α = 1, řízená soustava F2(s)) Vystup ze soustavy 4 2 y r,
0 r(t) y s (t)
-2
y b(t) -4 80
100
120
140
160 180 t (s) Akcni zasah
200
220
240 us (t) ub(t)
5
u
0
-5 80
100
120
140
160 t (s)
180
200
220
240
Obrázek 5.4: Průběhy MPC se znalostí referenční trajektorie (ys(t), us(t)) a bez znalosti referenční trajektorie (yb(t), ub(t) při regulaci soustavy s neminimální fází F3(s) (Ts = 0,5 s; P = 50, M = 25, α = 1)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
49
Vysoké učení technické v Brně
Vystup ze soustavy 8 6 y r,
r(t) y (t) s
4
y (t) b
2 0
0
5
10
15
20 t (s) Akcni zasah
25
30
35
10 u (t) s
u
5
0
ub(t)
0
5
10
15
20
25
30
35
t (s)
Obrázek 5.5: Odezva MPC na lineárně narůstající referenční trajektorii se znalostí referenční trajektorie (ys(t), us(t)) a bez znalosti referenční trajektorie (yb(t), ub(t) (Ts = 0,2 s; P = 50, M = 25, α = 1, řízená soustava F2(s)) Vystup ze soustavy 0.08 0.06 0.04
y r,
0.02 0
0
5
10
15 t (s) Akcni zasah
20
25
30
0
5
10
15 t (s)
20
25
30
0.4
0.2 u 0
Obrázek 5.6: Odezva MPC na lineárně narůstající poruchu o směrnici k = 0,25 (Ts = 0,2 s; P = 50, M = 25, α = 1, řízená soustava F2(s))
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
50
Vysoké učení technické v Brně
6.
POROVNÁNÍ MPC S PSD REGULÁTOREM
Po úspěšném vložení sumační složky do MPC je vhodné jej porovnat s PSD regulátorem. Pro porovnání regulátorů byly použity soustavy: F2 (s ) = F4 (s ) = F5 (s ) = F6 (s ) =
1
(10s + 1)(s + 1)2 1
s (s + 1)
(6.1)
(6.2)
2
1,2
(3s + 1)(s + 1)2 1,4
(10s + 1)(2s + 1)2
(6.3)
(6.4)
Jako pevně nastavený regulátor byl použit diskrétní PSD regulátor s dynamickým omezením přebuzení a filtrací derivační složky daný rovnicí Ts z −1 1 − z −1 FR ( z ) = K R 1 + +L TN −1 − s T ( 1 − z ) I TD −1 1− e z
T L= D Ts
TN − s 1 − e TD
(6.5)
(6.6)
kde Ts je perioda vzorkování, KR je zesílení regulátoru, TI je integrační konstanta, TD je derivační konstanta, N je filtrační koeficient. Vztah (6.6) byl získán z [16].
6.1
POROVNÁNÍ MPC S PSD REGULÁTOREM PŘI SKOKOVÉ ZMĚNĚ PORUCHY
Aby bylo porovnání adekvátní, je u MPC vypnuta znalost budoucího průběhu referenční trajektorie. A je použito pouze omezení velikosti akčního zásahu − 10 ≤ u ≤ 10 tak, jako u použitého PSD regulátoru. Model řízené soustavy byl získán offline před začátkem regulace. Nejdříve byla použita soustava F2(s) (6.1). Parametry PSD regulátoru byly nejdříve hrubě nastaveny metodou Ziegler - Nichols pro omezení kmitavého průběhu [15] a pak doladěny metodou pokus omyl tak, aby byl přechodný děj co nejkratší. Pro soustavu F2(s) jsou parametry PSD regulátoru
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
51
Vysoké učení technické v Brně
následující: KR = 6,5, TI = 5 s, TD = 0,75 s, Tt = 10 s a N = 3 při Ts = 0,5 s. Pomocí doladění parametrů PSD regulátoru bylo docíleno o něco rychlejšího přechodného děje než u MPC (obrázek 6.1). Ovšem je nutné upozornit, že MPC v tomto případě nezná budoucí průběh referenční trajektorie. MPC je schopen vyregulovat poruchu o něco rychleji než PSD regulátor. PSD regulátor by šel nastavit i tak, aby poruchu vyreguloval rychleji než MPC, pak lze ovšem očekávat, že odezva na změnu žádané hodnoty bude pomalejší a více kmitavá. Vystup ze soustavy 4 3 r(t) v(t) y MPC(t)
2 1 y r,
0
y PSD(t)
-1 -2 -3 40
60
80
100
120 t (s) Akcni zasah
140
160
180
uMPC(t)
10
uPSD(t) 0 u -10 40
60
80
100
120
140
160
180
t (s)
Obrázek 6.1: Porovnání PSD regulátoru (yPSD(t), uPSD(t)) s MPC (Ts = 0,5 s; P = 50, M = 25, α = 0,2) (yMPC(t), uMPC(t)) bez znalosti referenční trajektorie při regulaci soustavy F2(s) Při regulaci soustavy s astatismem F4(s) (6.2) byl PSD nastaven podle metody Ziegler - Nichols pro omezení kmitavého průběhu: KR = 0,5, TI = 6,9 s, TD = 0,86 s, Tt = 13,8 s a N = 3 při Ts = 0,2 s. Protože se jedná o soustavu na mezi stability, je
řízení soustavy s astatismem náročnější. Průběhy jsou zobrazeny na obrázku 6.2. Zde již MPC dává mnohem lepší odezvy s menším překmitem než PSD regulátor.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
52
Vysoké učení technické v Brně
Vystup ze soustavy 6
r(t) v(t) y MPC(t)
4 2 y r,
y PSD(t)
0 -2 -4 40
60
80
100
120 t (s) Akcni zasah
140
160
180 uMPC(t)
10
uPSD(t) u
0
-10 40
60
80
100
120
140
160
180
t (s)
Obrázek 6.2: Porovnání PSD regulátoru (yPSD(t), uPSD(t)) s MPC (Ts = 0,2 s; P = 50, M = 25, α = 0,2) (yMPC(t), uMPC(t)) bez znalosti referenční trajektorie při regulaci soustavy s astatismem F4(s) 6.2
POROVNÁNÍ MPC S PSD REGULÁTOREM PŘI ZMĚNĚ DYNAMIKY SOUSTAVY
Změna dynamiky soustavy se provádí skokovou změnou jedné řízené soustavy na jinou v t = 640 s. Nejdříve byla vyzkoušena možnost přepnutí soustavy z rychlejší dynamikou F5(s) (6.3) na pomalejší F2(s) (6.1). Pevně nastavený PSD měl parametry: KR = 1,8, TI = 5 s, TD = 0,4 s, Tt = 10 s a N = 3 při Ts = 0,5 s. Tyto parametry byly opět získány metodou Ziegler - Nichols pro omezení kmitavého průběhu, a pak již jen doladěny na kratší přechodný děj. Z obrázku 6.3 je vidět, že jak PSD, tak MPC jsou schopny i po změně dynamiky soustavy ji dále úspěšně řídit. Odezva na změnu referenční trajektorie je již ale více kmitavá. MPC je schopen regulovat změněnou soustavu s nulovou ustálenou odchylkou díky tomu, že je v něm přítomen sumační člen, který rozdíl mezi výstupem modelu soustavy a výstupem
řízené soustavy kompenzuje. Také je možné si všimnout, že při změně dynamiky (v t = 640 s) dává MPC více kmitavé akční zásahy než se ustálí.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
53
Vysoké učení technické v Brně
Vystup ze soustavy 5
y r,
r(t) y
0
(t)
MPC
y PSD(t) -5 500
550
600
650 700 t (s) Akcni zasah
750
800
850
10
0
u
uMPC(t) uPSD(t)
-10 500
550
600
650 t (s)
700
750
800
850
Obrázek 6.3: Porovnání PSD regulátoru (yPSD(t), uPSD(t)) s MPC (Ts = 0,5 s; P = 50, M = 25, α = 0,2) (yMPC(t), uMPC(t)) bez znalosti referenční trajektorie při změně soustavy z F5(s) na F2(s) v t = 640 s Vystup ze soustavy 6 4 2 y r,
0 -2
r(t) y MPC(t)
-4
y PSD(t) 500
550
600
650 700 t (s) Akcni zasah
750
800
850 uMPC(t) uPSD(t)
10
0 u -10 500
550
600
650 t (s)
700
750
800
850
Obrázek 6.4: Porovnání PSD regulátoru (yPSD(t), uPSD(t)) s MPC (Ts = 0,5 s; P = 50, M = 25, α = 0,2) (yMPC(t), uMPC(t)) bez znalosti referenční trajektorie při změně soustavy z F2(s) na F5(s) v t = 640 s
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
V případě, že proběhne změna dynamiky v opačném směru, z pomalejší soustavy F2(s) (6.1) na rychlejší F5(s) (6.3), je situace mnohem zajímavější (obrázek
6.4). Parametry PSD regulátoru jsou opět nastaveny pro omezení kmitání a následně doladěny na: KR = 6,5, TI = 5 s, TD = 0,75 s, Tt = 10 s a N = 3 při Ts = 0,5 s. I když má PSD regulátor velmi kmitavou odezvu, je stále ještě schopen regulovat soustavu. Za to MPC již ztrácí kontrolu nad řízenou soustavou. Při dalším zrychlení řízené soustavy selhává při regulaci i PSD regulátor. Z těchto průběhů se dá usoudit, že MPC je více náchylný na změnu dynamiky soustavy než pevně nastavený PSD regulátor. To je nejvíce ovlivněno tím, že model soustavy v MPC neodpovídá řízené soustavě a vypočtené akční zásahy nejsou adekvátní k řízené soustavě. Kmitání akčního zásahu nejvíce ovlivňuje přítomnost sumačního členu v MPC. Se změnou dynamiky si snadněji než MPC s pevně nastaveným modelem poradí adaptivní MPC, který provádí identifikaci průběžně i během regulace. Na
obrázku 6.5 je zobrazena stejná situace jak v předešlém případě (změna dynamiky z pomalejší soustavy F2(s) (6.1) na rychlejší F5(s) (6.3)), jen je použit adaptivní MPC. Pro adaptivní MPC byla použita průběžná identifikace s metodou učení L-M s parametrem δ = 0,001 a s počtem vzorů o 250 prvcích a 20 iterací. Adaptivní MPC se po změně řízené soustavy postupně adaptuje na novou soustavu F4(s). Při změně referenční trajektorie v t = 700 s není MPC ještě plně adaptován. To je převážně ovlivněno tím, že paměť vzorů u identifikační metody ještě není aktualizována pro soustavu F5(s). Počet vzorů byl nastaven na hodnotu 250, tzn. že při vzorkování po 0,5 s je v paměti posledních 125 s. Proto je při změně žádané hodnoty v t = 800 s již MPC správně adaptován. PSD regulátor je nastaven se stejnými parametry jako v předešlém případě. Aby se i PSD regulátor přizpůsobil změněné soustavě, je vhodné použít adaptivní PSD regulátor. Na obrázku 6.6 je zobrazena změna soustavy z F2(s) (6.1) na F6(s) (6.4) s tím rozdílem, že místo pevně nastaveného PSD je použit adaptivní PSD regulátor s identifikačním algoritmem založeným na metodě L-M. Nastavení identifikační
části je u PSD stejné jako u MPC. Výpočet parametrů PSD regulátoru se provádí metodou Ziegler – Nichols pro omezení kmitavého průběhu. Z průběhů je patrné, že se oba regulátory postupně adaptují na novou soustavu.
54
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
55
Vysoké učení technické v Brně
Vystup ze soustavy 6 4 2 y r,
0 -2
r(t) y MPC(t)
-4
y PSD(t) 500
550
600
650 700 t (s) Akcni zasah
750
800
850 uMPC(t) uPSD(t)
10
0 u -10 500
550
600
650 t (s)
700
750
800
850
obrázek 6.5: Porovnání PSD regulátoru (yPSD(t), uPSD(t)) s adaptivním MPC (Ts = 0,5 s; P = 50, M = 25, α = 0,2) (yMPC(t), uMPC(t)) bez znalosti referenční trajektorie při změně soustavy z F2(s) na F5(s) v t = 640 s Vystup ze soustavy 6 4
y r,
2
r(t) y
0
y PSD(t)
(t)
MPC
-2 -4 400
500
600
700 t (s) Akcni zasah
800
900
1000 uMPC(t)
10
uPSD(t)
0 u -10 400
500
600
700 t (s)
800
900
1000
Obrázek 6.6: Porovnání adaptivního PSD regulátoru (yPSD(t), uPSD(t)) s adaptivním MPC (Ts = 0,5 s; P = 50, M = 25, α = 0,2) (yMPC(t), uMPC(t)) při změně soustavy z F2(s) na F6(s) v t = 640 s
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
6.3
VIZUALIZAČNÍ PROSTŘEDÍ
Bylo vytvořeno jednoduché vizualizační prostředí, které umožňuje porovnání průběhů MPC s PSD regulátorem. Po odsimulování jsou zobrazeny výsledné průběhy výstupu a akčního zásahu. Program umožňuje nastavit parametry regulované soustavy do třetího řádu. Dále je možné volit amplitudu a periodu obdélníkové referenční trajektorie, také je možné nastavit skokovou poruchu působící na vstup soustavy. PSD regulátor je realizován s dynamickým omezením přebuzení a s filtrací derivační složky, přenos je dán vztahem (6.5), (6.6). Je možné určit zda MPC bude znát budoucí průběh referenční trajektorie. Další možností je použití adaptivní varianty MPC. Zde je před začátkem regulace řízená soustava zidentifikována metodou RLS a s takto získanými parametry soustavy je pak započata regulace. Adaptivní MPC používá pro průběžnou identifikaci během regulace metodu RLS s λ = 0,995 . U MPC je možné nastavit horizont predikce i
řízení, penalizaci změny akčního zásahu a omezení týkající se změny a velikosti akčního zásahu. Perioda vzorkování je společná pro MPC i PSD regulátor.
Obrázek 6.7: Hlavní obrazovka vizualizačního prostředí
56
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
57
Vysoké učení technické v Brně
7.
POUŽITÍ MPC PRO ŘÍZENÍ FYZIKÁLNÍCH MODELŮ
Celý systém, na kterém probíhalo řízení fyzikálních modelů se skládal z PC, kde v programu MATLAB běžel algoritmus MPC. Dále z programovatelného automatu B&R a k němu připojenému fyzikálnímu modelu. Jako fyzikální model byla použita soustava s třemi integrátory tvořenými operačními zesilovači. Pomocí přepínačů na předním panelu lze nastavovat časové konstanty přenosové funkce soustavy.
Omezení
pro
akční
zásah
fyzikálního
modelu
jsou
v rozsahu
− 10 V ≤ u ≤ 10 V . Použitý programový automat (PLC) B&R byl řady 2005. Obsahoval procesorovou jednotku CP360 a kartu analogových a digitálních vstupů a výstupů UM161. Na vstupu analogové karty byl 14ti bitový převodník. Propojení mezi PC a PLC bylo pomocí Ethernetu . V případě, že řídicí algoritmus běží v PC a PLC zajišťuje snímání potřebných dat z procesu a vytváří akční zásahy pro proces, hovoříme o soft real-time řízení (obrázek 7.1).
Analogové vstupy a výstupy
Ethernet
PC
PLC
fyz. model
obrázek 7.1: Soft real-time řízení fyzikálního modelu Řízeným soustavám odpovídají přibližně přenosy F2 (s ) = F7 (s ) = F8 (s ) =
1
(10s + 1)(s + 1)2 1,05
(0,5s + 1)(s + 1)2 1,05
(3,5s + 1)(s + 1)2
(7.1)
(7.2)
(7.3)
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
7.1
SKOKOVÁ ZMĚNA PORUCHY NA VSTUPU SOUSTAVY
Byly změřeny odezvy na poruchu vstupující do soustavy F8 (s ) (7.3). MPC měl v sobě implementována omezení pro velikost akčního zásahu − 10 V ≤ u ≤ 10 V a jeho změnu − 3 V ≤ ∆u ≤ 3 V . V případě pevně nastaveného MPC regulátoru, byl model soustavy získán off-line před začátkem regulačního děje. V případě skokové poruchy působící na vstupu soustavy je porucha bez problémů vyregulována (obrázek 7.2). V případě, že je na místo pevně nastaveného modelu v MPC použit adaptivní model s identifikací metodami RLS, L-M, je porucha také vyregulována
obrázek 7.3, 7.4). Pouze po příchodu poruchy se akční zásah mírně rozkmitá. Výstup je schopen stále sledovat referenční trajektorii. Kmitání akčního zásahu je způsobeno tím, že po příchodu poruchy se identifikační metody snaží aproximovat soustavu s poruchou, i když se ve skutečnosti přenos soustavy nezměnil. Model soustavy je ve výsledku naučen na jiný než na skutečný přenos soustavy. Podmínkou správné identifikace je zajištění nepůsobení poruchových signálů na soustavu, popřípadě je alespoň zajištěna možnost poruchu působící na soustavu měřit. V této kapitole jsou identifikační metody u adaptivních MPC nastaveny takto: RLS s
λ = 0,995 ; L-M s δ = 0,001 , počet iterací je nastaven na 20 a trénovací množina obsahuje 250 vzorů.
58
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
59
Vysoké učení technické v Brně
Vystup ze soustavy 3
U (V)
2 1 0 -1
r(t) v(t) y(t)
-2 -3 0
50
100
150
0
50
100
150
200 t (s) Akcni zasah
250
300
350
400
250
300
350
400
U (V)
10
0 -10 200 t (s)
Obrázek 7.2: Regulace pevně nastaveným MPC (Ts = 0,2 s; P = 50, M = 25, α = 0,2) soustavy F8(s) při příchodu poruchy v t = 130 s o amplitudě 1,5 V
Vystup ze soustavy 3 2 U (V)
1 0 -1
r(t) v(t) y(t)
-2 -3 0
50
100
150
0
50
100
150
200 t (s) Akcni zasah
250
300
350
400
250
300
350
400
U (V)
10
0
-10 200 t (s)
Obrázek 7.3: Regulace adaptivním MPC (RLS) (Ts = 0,2 s; P = 50, M = 25, α = 0,2) soustavy F8(s) při příchodu poruchy v t = 130 s o amplitudě 1,5 V
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
60
Vysoké učení technické v Brně
Vystup ze soustavy 4
U (V)
2
0 r(t) v(t)
-2
y(t) 0
50
100
150
0
50
100
150
200 t (s) Akcni zasah
250
300
350
400
250
300
350
400
U (V)
10
0
-10 200 t (s)
Obrázek 7.4: Regulace adaptivním MPC (L-M) (Ts = 0,2 s; P = 50, M = 25, α = 0,2) soustavy F8(s) při příchodu poruchy v t = 130 s o amplitudě 1,5 V 7.2
ZMĚNA DYNAMIKY SOUSTAVY
Byla ověřena schopnost MPC regulovat soustavu, které se v přibližně v t = 130 s změní dynamika. Při pevně nastaveném MPC byla nejdříve změněna
dynamika ze soustavy F2 (s ) (7.1) na rychlejší F8 (s ) (7.3). Prediktivní regulátor je i přes více kmitavý akční zásah stále schopen regulovat soustavu (obrázek 7.5). V případě změny ze soustavy F2 (s ) (7.1) na rychlejší F7 (s ) (7.2) již prediktivní regulátor o tuto schopnost přichází a soustavu není možné úspěšně řídit (obrázek
7.6). Při použití adaptivní varianty MPC je regulace i po změně soustavy z F2 (s ) (7.1) na rychlejší F7 (s ) (7.2) úspěšná (obrázek 7.7, 7.8). Parametry modelu se po změně přizpůsobují nové dynamice soustavy. U adaptivních MPC probíhá prvních 20 kroků počáteční nastavení modelu, v této době je soustava řízena P regulátorem o zesílení rovnu 1. A až po 20 krocích je k soustavě připojen MPC, který převezme
řízení.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
61
Vysoké učení technické v Brně
Vystup ze soustavy 3 2 U (V)
1 0 -1 r(t) y(t)
-2 -3 0
50
100
150
200
250
150
200
250
t (s) Akcni zasah
U (V)
10
0
-10 0
50
100 t (s)
Obrázek 7.5: Regulace pevně nastaveným MPC (Ts = 0,2 s; P = 50, M = 25, α = 0,2) při změně soustavy z F2(s) na F8(s) v t = 130 s Vystup ze soustavy 3 2 U (V)
1 0 -1 r(t) y(t)
-2 -3 0
50
100
150
200
250
150
200
250
t (s) Akcni zasah
U (V)
10
0
-10 0
50
100 t (s)
Obrázek 7.6: Regulace pevně nastaveným MPC (Ts = 0,2 s; P = 50, M = 25, α = 0,2) při změně soustavy z F2(s) na F7(s) v t = 130 s
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
62
Vysoké učení technické v Brně
Vystup ze soustavy 4
U (V)
2 0 r(t) y(t)
-2 0
50
100
150
200
250
150
200
250
t (s) Akcni zasah
U (V)
10
0
-10 0
50
100 t (s)
obrázek 7.7: Regulace adaptivním MPC (RLS) (Ts = 0,2 s; P = 50, M = 25, α = 0,2) při změně soustavy z F2(s) na F7(s) v t = 130 s Vystup ze soustavy 3 2 U (V)
1 0 -1 r(t) y(t)
-2 -3 0
50
100
150
200
250
150
200
250
t (s) Akcni zasah
U (V)
10
0
-10 0
50
100 t (s)
obrázek 7.8: Regulace adaptivním MPC (L-M) (Ts = 0,2 s; P = 50, M = 25, α = 0,2) při změně soustavy z F2(s) na F7(s) v t = 130 s
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií
63
Vysoké učení technické v Brně
8.
ZÁVĚR
Po nastudování problematiky prediktivního řízení a seznámení se s metodikou identifikace byl vytvořen identifikační algoritmus založený na metodě nejmenších
čtverců (RLS). Po seznámení se s neuronovými sítěmi (kapitola 3) byly vytvořeny identifikační
algoritmy
s metodami
učení
back – propagation
(B-P)
a
Levenberg - Marquardt (L-M). Po vzájemném porovnání podle uvedeného nastavení nejlépe vychází metody RLS a L-M, které dávají velmi podobné výsledky vzhledem k rychlosti konvergence parametrů modelu soustavy. Při velmi krátké periodě vzorkování se stává algoritmus RLS numericky nestabilní. Nejdříve byl v prostředí MATLAB/Simulink otestován a implementován analyticky řešený algoritmus MPC bez omezení s pevně nastaveným modelem získaným off-line metodou nejmenších čtverců. Dále byl algoritmus optimalizace akčních zásahů řešen numerickou metodou (využití funkce quadprog), aby bylo možné optimálně vyřešit účelovou funkci s danými omezeními. Také byly vytvořeny adaptivní varianty MPC využívající pro identifikaci soustavy metodu RLS a L-M. Pro odstranění ustálené odchylky jsou vyzkoušeny dva různé přístupy (kapitola 5). Zvláště výhodný se jeví použití MPC se sumačním členem na výstupu regulátoru, který je možné použít i při znalosti budoucích hodnot referenční trajektorie. U obou přístupů pro odstranění ustálené odchylky je zaručeno, že akční zásahy jsou stále optimální podle daného kritéria. Pro rychlejší vyregulování ustálené odchylky je možné zkrátit periodu vzorkování. Ale pro zachování stejné dynamiky přechodného děje je nutné adekvátně prodloužit i horizont predikce. Tím ale vzrostou rozměry matic vstupujících do optimalizačního algoritmu a zvětší se tak výpočetní nároky prediktivního regulátoru. Je proto nezbytné hledat určitý kompromis mezi těmito protichůdnými faktory. V kapitole 6 je porovnání PSD regulátoru s MPC se sumačním členem na výstupu regulátoru. Z průběhů je vidět, jak si oba regulátory dokáží poradit se skokovou změnou poruchy. Na změnu dynamiky soustavy je MPC více náchylný než PSD regulátor. V další kapitole je pak prezentována schopnost MPC úspěšně řídit fyzikální modely.
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
9. [1]
LITERATURA
CAMACHO, E. F. – BORDOS, C.: Model Predictive Control. Springer, London 1999. ISBN 3-540-76241-8.
[2]
RICHALET, J.– RAULT, A.– TESTUD, J. L.– PAPON, J.: Model Predictive Heuristic Control: Applications to Industrial Processes. Automatica, 14,
s. 413–428, 1978. [3]
CUTLER, C. R. – RAMAKER, D. L.: Dynamic matrix control – a computer control algorithm. Joint Automatic Control Conference, San Francisco,
California, 1980. [4]
CUTLER, C. R.: Quadratic-Program Dynamic Matrix Control. AIChE National Meeting, Houston, 1983.
[5]
CLARKE, D. W. – MOHTADI, C. – TUFFS, P. S.: Generalized Predictive Control. Part I. The Basic Algorithm. Automatica, 23, s. 137–148, 1987.
[6]
KWON, W.H. – BYUN, L.: Receding Horizont Tracking Control. Automatica, 25, s. 149–163, 1989.
[7]
NEPEVNÝ, P.: Prediktivní regulátory s prvky umělé inteligence. Pojednání o disertační práci. VUT, Brno, duben 2006.
[8]
PIVOŇKA, P.: Optimalizace regulátorů. Skriptum. VUT, Brno, 2005. Elektronický učební text.
[9]
DRÁBEK, O. – SEIDL, P. – TAUFER, I.: Umělé neuronové sítě: Základy teorie a aplikace. CHEMagazín: Časopis pro chemicko-technologickou a
laboratorní praxi. 2005-2006, ročník 15-16. Dostupné i z URL:
[10] ZELINKA, I: Umělá inteligence 1: Neuronové sítě a genetické algoritmy. VUTIUM, Brno, 1998. ISBN 80-214-1163-5. [11] NOVÁK, M. a kol.:Umělé neuronové sítě teorie a aplikace. C. H. BECK, Praha, 1998. 382 s. ISBN 80-7179-132-6. [12] VELEBA, V., PIVOŇKA, P.: Adaptive Controller with Identification Based on Neural Network for Systems with Rapid Sampling Rates. WSEAS Transactions
64
ÚSTAV AUTOMATIZACE A MĚŘICÍ TECHNIKY Fakulta elektrotechniky a komunikačních technologií Vysoké učení technické v Brně
on Systems, 2005, roč. 4, č. 4, str. 385 - 388, ISSN 1109-2777. Dostupné i z URL: [13] HAVLENA, V.: Moderní teorie řízení. Doplňkové skriptum. Skripta ČVUT, FEL 2002. ISBN 80-01-02036-3. Dostupné i z URL: [14] LORENC, J.: LQ regulátor. Návod do cvičení předmětu MOPR [15] PIVOŇKA, P.: Číslicová řídicí technika. Skriptum. VUT, Brno, 2003. Elektronický učební text. [16] SCHMIDT, M.:Diskrétní realizace derivační složky podle požadavků na tvar přechodové charakteristiky. VUT, Brno
65