ˇ
´ vysoke
´ uc
ˇen´ı technicke
´ v Praze
Cesk
e
´
Fakulta elektrotechnicka
´ PRACA
´
DIPLOMOVA
Tvorba model˚
u a prediktivn´ı ˇ
r´ızen´ı budov
Praha, 2011
ˇ aˇ
Autor: Eva Z´
cekov´
a
Prohl´
aˇ
sen´ı
Prohlaˇsuji, ˇze jsem svou diplomovou pr´aci vypracovala samostatnˇe a pouˇzila jsem
pouze podklady (literaturu, projekty, SW atd.) uveden´e v pˇriloˇzen´em seznamu.
V Praze dne
podpis
i
Pod’akovanie
Na tomto mieste by som rada pod’akovala predovˇsetk´
ym ved´
ucemu svojej diplomovej
pr´ace Ing. Samuelovi Pr´ıvarovi za v´
yborne podmienky, ktor´e mi poskytol pri p´ısan´ı tejto
ˇ
pr´ace, za mnoˇzstvo cenn´
ych r´ad a pripomienok. Dalej
moja vd’aka patr´ı rodine za podporu
poˇcas cel´eho ˇst´
udia.
ii
Abstrakt
Pri obrovskom mnoˇzstve energie, ktor´e je vynakladan´e na vykurovanie budov, je
z hl’adiska zn´ıˇzenia n´akladov dˆoleˇzit´e zefekt´ıvnit’ proces vykurovania. Konvenˇcn´e met´ody
riadenia nie s´
u schopn´e zohl’adnit’ n´ahle zmeny poˇcasia, poˇziadavky na riadenie a pod.
Mnoho t´
ychto nedostatkov mˆoˇze odstr´anit’ pouˇzitie predikt´ıvneho regul´atora (MPC).
Ten m´a okrem mnoˇzstva v´
yhod, medzi ktor´e patri moˇznost’ zahrn´
ut’ obmedzenia
priamo do samotn´eho z´akona riadenia ˇci jednoduch´a aplikovatel’nost’ aj na MIMO syst´emy,
aj nev´
yhody, a to predovˇsetk´
ym fakt, ˇze k svojmu fungovaniu potrebuje matematick´
y
model syst´emu, ktor´
y ˇco najpresnejˇsie popisuje spr´avanie re´alneho syst´emu. V tejto
pr´aci bude preto zv´
yˇsen´a pozornost’ venovan´a hl’adaniu matematick´eho modelu syst´emu
ˇ eho
vhodn´eho pre pouˇzitie s predikt´ıvnych regul´atorom. Pri tvorbe modelu budovy Cesk´
ˇ
vysok´eho uˇcen´ı technick´eho (CVUT)
v Praze pouˇzijeme okrem klasick´
ych identifikaˇcn´
ych
met´od aj identifikaˇcn´e met´ody minimalizuj´
uce viackrokov´
u predikˇcn´
u chybu, ktor´e poskytuj´
u model s dobr´
ymi predikˇcn´
ymi vlastnost’ami, ktor´
y je vhodn´
y pre pouˇzitie s MPC.
Vybran´
y model bude n´asledne pouˇzit´
y s navrhnut´
ym predikt´ıvnym regul´atorom pre riaˇ
denie teploty v miestnostiach budovy CVUT.
iii
Abstract
Due to a huge amount of energy put on the heating in buildings, it is important
to make the building heating process more efficient in order to lower the costs. Conventional control methods are not able to take into account fast weather changes, control
requirements and others. Many of these drawbacks can be eliminated by making use
of predictive controller (MPC).
Besides plenty of advantages like an ability to include given constraints directly
into the control law or simple aplicability with MIMO systems, it has also some disadvantages, for examples it needs a mathematical model of system describing its behaviour as
accurately as posible in order to work properly. In this thesis, increased attention is paid
to searching for mathematical model of system suitable for use with predictive controller. For construnction of Czech Technical University (CTU) in Prague building model,
both classical methods and methods minimizing multistep ahead prediction error providing with a model with good prediction properties are used. The chosen model is then
employed with desinged predictive controller for temperature control in CTU building
rooms.
iv
Obsah
Zoznam obr´
azkov
viii
Zoznam tabuliek
x
´
1 Uvod
1
2 Popis riaden´
eho syst´
emu
3
3 Identifikaˇ
cn´
e met´
ody
7
3.1
Met´ody Grey box identifik´acie . . . . . . . . . . . . . . . . . . . . . . . .
9
3.1.1
Tvorba ˇstrukt´
ury modelu
. . . . . . . . . . . . . . . . . . . . . .
10
3.1.2
Odhad parametrov . . . . . . . . . . . . . . . . . . . . . . . . . .
13
Subspace identifikaˇcn´e met´ody . . . . . . . . . . . . . . . . . . . . . . . .
14
4 Identifik´
acia pre predikt´ıvne riadenie
´
4.1 Uvod
do MRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
3.2
4.2
4.3
4.4
MRI identifik´acia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
4.2.1
MRI u
´ˇcelov´a funkcia . . . . . . . . . . . . . . . . . . . . . . . . .
20
4.2.2
Dvojkrokov´
y algoritmus . . . . . . . . . . . . . . . . . . . . . . .
23
Jednokrokov´
y algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
4.3.1
Implement´acia (popis) algoritmu . . . . . . . . . . . . . . . . . .
30
MRI identifik´acia v kombin´acii s PLS . . . . . . . . . . . . . . . . . . . .
ˇ
4.4.1 Ciastoˇ
cn´a line´arna regresia . . . . . . . . . . . . . . . . . . . . . .
33
33
4.4.2
36
Popis algoritmu . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 Predikt´ıvne riadenie
5.1
19
38
N´avrh predikt´ıvneho regul´atora . . . . . . . . . . . . . . . . . . . . . . .
41
5.1.1
43
Predikt´ıvny regul´ator s obmedzen´ım . . . . . . . . . . . . . . . .
vi
5.1.1.1
Sledovanie konˇstantnej referencie . . . . . . . . . . . . .
Stabilita a rieˇsitel’nost’ . . . . . . . . . . . . . . . . . . .
46
MPC - formul´acia probl´emu . . . . . . . . . . . . . . . . . . . . . . . . .
49
5.1.1.2
5.2
6 Simul´
acie na re´
alnych datach
6.1
45
52
Identifik´acia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
6.1.1
Subspace identifikaˇcn´e met´ody . . . . . . . . . . . . . . . . . . . .
53
6.1.2
Grey box modelovanie . . . . . . . . . . . . . . . . . . . . . . . .
54
6.1.3
MRI identifikaˇcn´e met´ody . . . . . . . . . . . . . . . . . . . . . .
55
6.1.4
Grey box modelovanie v kombin´acii s MRI . . . . . . . . . . . . .
57
6.1.5
Kombin´acia MRI a PLS . . . . . . . . . . . . . . . . . . . . . . .
58
6.1.6
Porovnanie identifikovan´
ych modelov . . . . . . . . . . . . . . . .
59
6.2
Predikt´ıvny regul´ator . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
6.3
Zahrnutie vplyvu slnka . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
6.3.1
66
Simul´acie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7 Z´
aver
69
Literat´
ura
70
A Obsah priloˇ
zen´
eho CD
I
vii
Zoznam obr´
azkov
3
2.2
ˇ
Pohl’ad na budovu CVUT.
. . . . . . . . . . . . . . . . . . . . . . . . . .
ˇ
N´aˇcrt budovy CVUT.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3
Zjednoduˇsen´a sch´ema vykurovacieho syst´emu. . . . . . . . . . . . . . . .
4
2.4
Sch´ema riaden´eho objektu. . . . . . . . . . . . . . . . . . . . . . . . . . .
5
3.1
RC diagram pre rovnicu vedenia tepla. . . . . . . . . . . . . . . . . . . .
11
3.2
RC diagram pre blok s dvoma vykurovac´ımi okruhmi. . . . . . . . . . . .
11
4.1
Hodnota optim´alneho koeficientu ˇsumov´eho modelu v z´avislosti na P . . .
26
4.2
Priebehy RP v z´avislosti na P . . . . . . . . . . . . . . . . . . . . . . . . .
27
4.3
Hodnota optim´alneho koeficientu ˇsumov´eho modelu v z´avislosti na P . . .
28
4.4
Nyquistove charakteristiky pre rˆozne P . . . . . . . . . . . . . . . . . . . .
28
4.5
Priebeh pomeru Rp v z´avislosti na P . . . . . . . . . . . . . . . . . . . . .
32
4.6
Priebehy P -krokov´
ych predikci´ı v´
ystupov. . . . . . . . . . . . . . . . . .
32
5.1
ˇ
Strukt´
ura procesu riadenia. . . . . . . . . . . . . . . . . . . . . . . . . .
39
5.2
Blokov´a sch´ema predikt´ıvneho regul´atoru. . . . . . . . . . . . . . . . . .
39
5.3
Posuvn´
y horizont predikcie. . . . . . . . . . . . . . . . . . . . . . . . . .
40
5.4
Porovnanie klasickej kriteri´alnej funkcie ,,zone control“. . . . . . . . . . .
50
6.1
Predikcie teplˆot v referenˇcn´
ych miestnostiach (4SID). . . . . . . . . . . .
54
6.2
Odozvy na jednotkov´
y skok na vstupoch (4SID). . . . . . . . . . . . . . .
54
6.3
Predikcie teplˆot v referenˇcn´
ych miestnostiach (GB). . . . . . . . . . . . .
55
6.4
Odozvy na jednotkov´
y skok na vstupoch (GB).
. . . . . . . . . . . . . .
55
6.5
Predikcie teplˆot v referenˇcn´
ych miestnostiach (MRI). . . . . . . . . . . .
56
6.6
Odozvy na jednotkov´
y skok na vstupoch (MRI). . . . . . . . . . . . . . .
57
6.7
Predikcie teplˆot v referenˇcn´
ych miestnostiach (MRI+GB). . . . . . . . .
57
6.8
Odozvy na jednotkov´
y skok na vstupoch (MRI+GB). . . . . . . . . . . .
58
6.9
Predikcie teplˆot v referenˇcn´
ych miestnostiach (MRI+PLS). . . . . . . . .
58
2.1
viii
4
6.10 Odozvy na jednotkov´
y skok na vstupoch (MRI+PLS).
. . . . . . . . . .
59
6.11 Porovnanie viackrokov´
ych predikˇcn´
ych ch´
yb identifikovan´
ych modelov. .
61
6.12 Porovnanie fitfaktorov identifikovan´
ych modelov. . . . . . . . . . . . . . .
61
6.13 Porovnanie ekvitermnej regul´acie a MPC (juh). . . . . . . . . . . . . . .
63
6.14 Porovnanie ekvitermnej regul´acie a MPC (sever). . . . . . . . . . . . . .
63
6.15 Umiestnenie ˇcidla intenzity slneˇcn´eho ˇziarenia. . . . . . . . . . . . . . . .
64
6.16 Korel´acia medzi ϑrs/n a I. . . . . . . . . . . . . . . . . . . . . . . . . . .
65
6.17 Porovnanie rozptylu predikˇcn´
ych ch´
yb identifikovan´
ych modelov. . . . . .
67
6.18 Porovnanie fitfaktorov identifikovan´
ych modelov. . . . . . . . . . . . . . .
68
6.19 Porovnanie predikci´ı identifikovan´
ych modelov. . . . . . . . . . . . . . . .
68
ix
Zoznam tabuliek
2.1
Meran´e veliˇciny riaden´eho syst´emu. . . . . . . . . . . . . . . . . . . . . .
6
3.1
V´
yznam symbolov. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
6.1
Porovnanie u
´spor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
x
Kapitola 1
´
Uvod
Podl’a [33] predstavuje energia spotrebovan´a v budov´ach 20% − 40% celkovej spotreby
energie, priˇcom takmer 50% z tejto spotrebovanej energie tvoria n´aklady na vykurovanie,
ventil´aciu a klimatiz´aciu (HVAC - Heating Ventilation and Airconditioning), potom
nasleduj´
u n´aklady na ohrev u
´ˇzitkovej vody a svetlo. Aj napriek tomu, ˇze v posledn´
ych
rokoch sa v´
yrazne vylepˇsuje technick´
y stav prostriedkov HVAC, je v tejto oblasti st´ale
vel’k´
y priestor pre d’al’ˇs´ı v´
yskum a zlepˇsenie.
Najbeˇznejˇsie pouˇz´ıvanou strat´egiou pri regul´acii vykurovac´ıch syst´emov je tzv. ekvitermn´a regul´acia, kde regul´ator na z´aklade u
´daju o nameranej vonkajˇsej teplote podl’a
ekvitermnej krivky pre pr´ısluˇsn´
u poˇzadovan´
u teplotu v referenˇcnej miestnosti urˇc´ı teplotu vykurovacej vody (pr´ıpadne inej veliˇciny priamo u
´mernej mnoˇzstvu tepelnej energie
dodanej do s´
ustavy). Tak´eto riadenie m´a vˇsak niekol’ko nedostatkov, nespr´avnym naladen´ım ekvitermnej krivky mˆoˇze doch´adzat’ k nedodrˇzaniu teplotn´eho komfortu, pr´ıpadne
k zbytoˇcn´emu prek´
ureniu a zv´
yˇsenej spotrebe energie.
Ani v pr´ıpade spr´avneho naladenia vˇsak nemus´ı viest’ ekvitermn´a regul´acia k poˇzadovan´emu v´
ysledku, pretoˇze berie do u
´vahy len aktu´alnu vonkajˇsiu teplotu. Poveternostn´e podmienky sa mˆoˇzu v priebehu okamˇziku v´
yrazne zmenit’ a pri vel’k´
ych budov´ach
s ˇcasov´
ymi konˇstantami v r´adoch hod´ın mˆoˇze vel’mi jednoducho doch´adzat’ k prek´
ureniu
alebo nedok´
ureniu.
Tu sa naskyt´a pouˇzitie alternat´ıvneho pr´ıstupu, ktor´
y rieˇsi probl´em vykurovania budovy ako optimalizaˇcn´
u u
´lohu. Ide o predikt´ıvny regul´ator (MPC - Model Predictive
Control), ktor´
y pouˇz´ıva matematick´
y model syst´emu k predikcii bud´
uceho spr´avania
budovy ber´
uc do u
´vahy predpoved’ poˇcasia (pr´ıpadne predpoved’ slneˇcn´eho osvitu ˇci obsadenosti budovy). Predikcie v´
ystupov s´
u n´asledne pouˇzit´e pri hl’adan´ı minima u
´ˇcelovej
funkcie, ktor´a by mala byt’ volen´a tak, aby penalizovala na jednej strane nedostatoˇcn´
y
1
´
KAPITOLA 1. UVOD
2
tepeln´
y komfort a na strane druhej mnoˇzstvo energie vynaloˇzenej na vykurovanie.
Pre optim´alne fungovanie predikt´ıvneho regul´atora je dˆoleˇzit´a vhodn´a vol’ba u
´ˇcelovej
funkcie, ktor´a je podmienen´a dobrou znalost’ou riaden´eho procesu. Ovel’a zloˇzitejˇsou a
komplexnejˇsou u
´lohou je ale vo v¨aˇcˇsine pr´ıpadov n´ajst’ model procesu, ktor´
y dok´aˇze ˇco
najvernejˇsie predikovat’ bud´
uce spr´avanie (stavov´e, pr´ıpadne v´
ystupne veliˇciny) na dobu
cel´eho horizontu predikcie. Beˇzne pouˇz´ıvan´e identifikaˇcn´e met´ody [29] poskytuj´
u modely,
ktor´e s´
u optimalizovan´e len v zmysle jednokrokovej predikˇcnej chyby. Predikcie t´
ychto
modelov s´ıce odpovedaj´
u skutoˇcnosti na niekol’ko m´alo krokov dopredu, viackrokov´e predikcie vˇsak uˇz nie s´
u aˇz tak´e presn´e a ich pouˇzitie v kombin´acii s MPC vedie na suboptim´alne spr´avanie.
Rieˇsen´ım je pouˇzit’ met´ody, ktor´e minimalizuj´
u viackrokov´
u predikˇcn´
u chybu [16, 17,
25, 40, 41] a poskytuj´
u modely s uspokojiv´
ymi predikˇcn´
ymi vlastnost’ami aj pri vyˇsˇs´ıch
predikˇcn´
ych horizontoch a st´avaju sa tak vhodn´
ymi kandid´atmi pre nasadenie s predikt´ıvnym regul´atorom.
ˇ eho vysok´eho uˇcen´ı technick´eho
Ciel’om naˇsej pr´ace je navrhn´
ut’ pre budovu Cesk´
v Praze v Dejviciach predikt´ıvny regul´ator zaloˇzen´
y na modeli. Ked’ˇze sa vˇsak tak´
yto
regul´ator nezaob´ıde bez matematick´eho modelu syst´emu, ktor´
y mu poskytne ˇco najpresnejˇsie dlhodob´e predikcie bud´
uceho spr´avania, bude v¨aˇcˇsia ˇcast’ naˇsej pr´ace venovan´a
pr´ave hl’adaniu vhodn´eho modelu. V kapitole 2 sa obozn´amime s riaden´
ym syst´emom, uvedieme prehl’adn´
y popis princ´ıpu pouˇzit´eho vykurovacieho syst´emu a vysvetl´ıme v´
yznam
jednotliv´
ych procesn´
ych veliˇc´ın dˆoleˇzit´
ych pre identifik´aciu a n´asledn´e riadenie.
Kapitola 3 prin´aˇsa struˇcn´
y prehl’ad vybran´
ych identifikaˇcn´
ych met´od a je diskutovan´a
moˇznost’ pouˇzitia t´
ychto met´od pri identifik´acii modelu budovy. V tejto kapitole je d’alej
ˇ
odvoden´a fyzik´alna ˇstrukt´
ura vybran´eho bloku budovy CVUT.
Rozsiahla samostatn´a kapitola sa zaober´a identifikaˇcn´
ym metod´am minimalizuj´
ucim
viackrokov´
u predikˇcn´
u chybu, ktor´e poskytuj´
u modely vhodn´e pre MPC. Okrem podrobn´eho popisu vybr´anych met´od viackrokovej identifik´acie s´
u uveden´e ilustraˇcn´e pr´ıklady
demonˇstruj´
uce vlastnosti t´
ychto met´od.
ˇ ’ 5 venovan´a predikt´ıvnemu riadeniu obsahuje okrem teoretick´eho u
Cast
´vodu do probˇ
lematiky predikt´ıvneho riadenia aj n´avrh konkr´etneho regul´atoru pre budovu CVUT.
Dosiahnut´e v´
ysledky su prezentovan´e v kapitole 6. T´ato kapitola obsahuje s´
ubor simul´acii pre konkr´enty vybran´
y blok, ktor´
y zah´rn
ˇa porovnanie vybran´
ych identifikaˇcn´
ych
met´od, zhodnotenie u
´speˇsnosti navrhnut´eho MPC a anal´
yzu vplyvu intenzity slneˇcn´eho
ˇziarenia na kvalitu identifikovan´
ych modelov.
Kapitola 2
Popis riaden´
eho syst´
emu
V naˇsej pr´aci sa budeme zaoberat’ identifik´aciou a n´asledn´
ym predikt´ıvnym riaden´ım
ˇ eho vysok´eho
budov. Budova, ktor´a bola poskytnut´a pre u
´ˇcely tejto pr´ace, je budova Cesk´
uˇcen´ı technick´eho v Praze. Konkr´etne ide o osemposchodov´
u budovu Fakulty strojnej a
Fakulty elektrotechnickej1 , ktor´a je zobrazen´a na obr. 2.1.
ˇ
Obr. 2.1: Pohl’ad na budovu CVUT.
ˇ
Budova CVUT
pozost´ava z viacer´
ych samostatn´
ych blokov, z ktor´
ych ˇcast’ je zateplen´a
(na obr. 2.2 zn´azornen´a ˇcervenou farbou) a druh´a ˇcast’ nezateplen´a (zobrazen´a modrou
farbou). Objektom naˇsej pr´ace boli bloky oznaˇcen´e B1, B2, B3, ktor´e pozost´avaj´
u z
dvoch nez´avisl´
ych vykurovac´ıch okruhov.
1
Technick´
a 2, Praha 6
3
´
´
KAPITOLA 2. POPIS RIADENEHO
SYSTEMU
4
ˇ
Obr. 2.2: N´
aˇcrt budovy CVUT.
Kaˇzd´
y z blokov m´a totoˇzn´
y princ´ıp stropn´eho vykurovania, tzv. ”Crittall”. Ide o teplovodn´
y syst´em, ktor´
y bol patentovan´
y v roku 1907 britsk´
ym profesorom Arthurom H.
Bakerom [1], v roku 1909 bol odk´
upen´
y firmou ”Crittall” a v roku 1927 pantentovo
vylepˇsen´
y R. G. Crittallom a J. L. Musgravom [7]. Syst´em bol popul´arny uˇz v obdob´ı
prvej svetovej vojny, najv¨aˇcˇsieho rozmachu dosiahol po druhej svetovej vojne, najm¨a vo
ˇ
Vel’kej Brit´anii. V b´
yvalom Ceskoslovensku
bol tento syst´em hojne nasadzovan´
y najm¨a
ˇ
koncom ˇsest’desiatych rokoch vo vel’k´
ych budov´ach. V pr´ıpade budovy CVUT
ide o meandrovit´
y syst´em potrub´ı, ktor´
y je zabudovan´
y v bet´onovom monolitickom strope. Zjednoduˇsen´a sch´ema stropn´eho vykurovacieho syst´emu pre jeden vykurovac´ı okruh je zobrazen´a na obr. 2.3.
teplota vykurovacej
vody
vonkajšia teplota
stropný
vykurovací
systém
referenčná
miestnosť
teplota v miestnosti
výmenník
tepla
zmiešavač
+
rozdeľovač
teplota vratnej vody
Obr. 2.3: Zjednoduˇsen´a sch´ema vykurovacieho syst´emu.
´
´
KAPITOLA 2. POPIS RIADENEHO
SYSTEMU
5
K ohrevu sekund´arnej vody priv´adzanej do rozdel’ovaˇca doch´adza v parnom v´
ymenn´ıku, kde je sekund´arna voda ohrievan´a parou privedenou z tepl´arenskej pr´ıpojky na
konˇstantn´
u hodnotu. Z v´
ymenn´ıka je voda veden´a do rozdel’ovaˇca, odkial’ je s´
ustavou
teplovodn´
ych r´
ur distribuovan´a do jednotliv´
ych vykurovac´ıch okruhov. Pred vstupom
do dan´eho okruhu je voda zmieˇsan´a v trojcestnom zmieˇsavacom ventile, ktor´eho poloha
regulovan´a n´ızko´
urovˇ
nov´
ym regul´atorom (PID) ud´ava zmieˇsavac´ı pomer medzi vodou privedenou z rozdel’ovaˇca a vratnou vodou z vykurovacieho okruhu. Referenˇcn´e hodnoty pre
n´ızko´
urovˇ
nov´e regul´atory s´
u nastavovan´e regul´atorom vyˇsˇsej u
´rovne, ktor´
ym je v komˇ
plexe budov CVUT
MPC. Vratn´a ochladen´a voda z jednotliv´
ych vykurovac´ıch okruhoch
sa vracia do zberaˇca a d’alej do parn´eho v´
ymenn´ıku, kde op¨at’ doch´adza k ohriatiu vody
a cel´
y cyklus sa opakuje.
Kaˇzd´
y z blokov B1, B2, B3, ktor´
ymi sa budeme v tejto pr´aci zaoberat’, m´a dva samostatn´e vykurovacie okruhy (severn´
u a juˇzn´
u vetvu). Na kaˇzdom okruhu je meran´a teplota vykurovacej, vratnej vody a teplota v referenˇcnej miestnosti. Zjednoduˇsen´a sch´ema
jedn´eho bloku, z ktorej budeme vych´adzat’ pri identifik´acii a n´aslednom riaden´ı, je zobrazen´a na obr. 2.4.
Obr. 2.4: Sch´ema riaden´eho objektu.
.
V´
yznam znaˇcenia uveden´eho na obr. 2.4 moˇzno n´ajst’ v tabul’ke 2.1. Tri z meran´
ych
veliˇc´ın s´
u vstupn´e veliˇciny: ϑo , ϑrws , ϑrwn . Teploty vratnej vody v jednotliv´
ych vykurovac´ıch okruhoch s´
u riadiace veliˇciny (ich hodnota je nastavovan´a regul´atorom). Vonkajˇsia
teplota predstavuje meratel’n´
u poruchu. Je dˆoleˇzit´e poznamenat’, ˇze v pr´ıpade ϑo neide
o aktu´alnu meran´
u hodnotu ale o predpoved’ na niekol’ko hod´ın dopredu (dan´e predikˇcn´
ym
horizontom MPC), ktor´a je poskytovan´a agent´
urou NOAA (National Oceanic and Atmo-
´
´
KAPITOLA 2. POPIS RIADENEHO
SYSTEMU
6
Tabul’ka 2.1: Meran´e veliˇciny riaden´eho syst´emu.
Oznaˇcenie Popis
ϑrs
teplota v referenˇcnej miestnosti - juh
ϑrn
teplota v referenˇcnej miestnosti - sever
ϑrws
teplota vratnej vody - juh
ϑrwn
teplota vratnej vody - sever
ϑo
vonkajˇsia teplota (predpoved’)
ϑhws
teplota vykurovacej vody - juh
ϑhwn
teplota vykurovacej vody - sever
spheric Administration)2 . Tieto predpovede vyuˇz´ıva predikt´ıvny regul´ator k presnejˇs´ım
predikci´am v´
ystupov, ktor´
ymi s´
u teplota vratnej vody a teploty v referenˇcn´
ych miestnostiach.
ˇ
Dalej
budeme uvaˇzovat’ syst´em s nasleduj´
ucim usporiadan´ım vstupn´
ych a v´
ystupn´
ych
veliˇc´ın:
iT
h
y = ϑrs ϑrws ϑrn ϑrwn ,
iT
h
u = ϑo ϑhws ϑhwn .
(2.1)
Naˇsou u
´lohou je teda n´ajst’ model syst´emu, ktor´
y bude ˇco najvernejˇsie popisovat’ vzt’ahy
medzi vstupn´
ymi a v´
ystupn´
ymi veliˇcinami (2.1) a poskytne ˇco najpresnejˇsie predikcie
v´
ystupov.
2
http://www.noaa.gov
Kapitola 3
Identifikaˇ
cn´
e met´
ody
Ako uˇz bolo spomenut´e v u
´vode tejto pr´ace, znalost’ modelu riaden´eho syst´emu, ktor´
y
ˇco najvernejˇsie popisuje realitu, je nevyhnutn´a pre optim´alne fungovanie predikt´ıvneho
riadenia. Preto je potrebn´e venovat’ hl’adaniu modelu syst´emu rovnako vel’k´
u pozornost’
ako n´avrhu samotn´eho z´akona riadenia.
Existuje vel’k´e mnoˇzstvo rˆozn´
ych pr´ıstupov k identifik´acii dynamick´
ych syst´emov,
ktor´
ym je venovan´
y nemal´
y priestor v odbornej literat´
ure [29, 42, 48]. Prvou vol’bou,
ktor´
u mus´ıme pred samotnou identifik´aciou uˇcinit’, je v´
yber reprezent´acie identifikovan´eho
modelu. V pr´ıpade ˇcasovo invariantn´
ych line´arnych syst´emov, (budeme pre ne pouˇz´ıvat’
oznaˇcenie LTI - Linear Time Invariant), na ktor´e sa v tejto pr´aci obmedz´ıme, je pomerne
beˇzne pouˇzivan´
y Input - Output model. Najm¨a v pr´ıpade SISO syst´emov pon´
uka hned’
niekol’ko ˇstrukt´
ur od najjednoduchˇsieho FIR modelu cez ARX aˇz po zloˇzitejˇsie ˇstrukt´
ury
ako ARMAX ˇci Box Jenkins, ktor´e dovol’uj´
u modelovat’ nez´avisle nielen parametre derministickej ˇcasti syst´emu, ale aj ˇcast’ stochastick´
u [20].
V pr´ıpade identifik´acie syst´emov s vel’k´
ym poˇctom vstupov a v´
ystupov, medzi ktor´e
budovy urˇcite patria, je pouˇzitie tak´
ychto ˇstrukt´
ur znaˇcne obmedzen´e. Ovel’a v´
yhodnejˇsie
je v takomto pr´ıpade pouˇzit’ reprezent´aciu stavov´
ym popisom, ktor´
u v ˇsest’desiatych rokoch minul´eho storoˇcia formuloval Rudolf Kalman a je doteraz s obl’ubou pouˇz´ıvan´a v rade
ˇ
praktick´
ych aplik´acii. Dalej
sa budeme zaoberat’ len identifik´aciou line´arnych ˇcasovo invariantn´
ych diskr´etnych stochastick´
ych syst´emov s nasleduj´
ucim stavov´
ym popisom:
x(k + 1) = Ax(k) + Bu(k) + v(k)
y(k) = Cx(k) + Du(k) + e(k)
(3.1)
ˇ
kde u(k) je deterministick´
y vstup modelu. Sum
merania v(k) a ˇsum procesu e(k) s´
u
stacion´arne n´ahodn´e postupnosti s nulovou strednou hodnotou (biely ˇsum), ktor´e s´
u
7
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
v rˆoznych ˇcasov nez´avisl´e a pre ktor´e plat´ı:
("
#)
v(k)
E
= 0,
e(k)
"
("
# ("
#))
#
# "
#T  "

v(k)
v(k)
Q S
v(k)
v(k) 
cov
,
=E
=
.
,
 e(k)
e(k)
e(k)
ST R
e(k) 
8
(3.2)
Matice Q ≥ 0 a R ≥ 0 predstavuj´
u kovarianciu ˇsumu procesu a ˇsumu merania. Matica S vyjadruje vz´ajomn´
u kovarianciu t´
ychto ˇsumov. V d’alˇsom texte budeme uvaˇzovat’
vz´ajomne nekorelovan´e ˇsumy merania a procesu, t.j. S = 0.
Z hl’adiska vyuˇz´ıvania apri´ornej inform´acie a experiment´alnej znalosti d´at moˇzno rozdelit’ identifikaˇcn´e met´ody do troch z´akladn´
ych skup´ın:
• White box
Vych´adza len z apri´ornej znalosti popisovan´eho procesu a vˇseobecne platn´
ych fyzik´alnych z´akonov. V´
ysledkom je matematick´
y model, ktor´eho parametre s´
u z´ıskan´e
anal´
yzou termodynamick´
ych, chemick´
ych a mechanick´
ym vlastnost´ı identifikovan´eho
objektu bez pouˇzitia experiment´alnych d´at. Ked’ˇze prakticky je nemoˇzn´e zostrojit’ presn´
y model tak zloˇzit´eho syst´emu, ako je budova, len zo znalosti fyzik´alnochemick´
ych vlastnost´ı, je pouˇzitie tohto pr´ıstupu pri identifik´acii budov nie pr´ıliˇs
vhodnou vol’bou a d’alej sa t´
ymto pr´ıstupom nebudeme zaoberat’.
• Black box
Je pr´ıstupom predstavuj´
ucim opaˇcn´
y extr´em ako White box modelovanie. Pri tvorbe
modelu syst´emu s´
u pouˇzit´e len experiment´alne z´ıskan´e d´ata bez vyuˇzitia akejkol’vek
znalosti o ˇstrukt´
ure ˇci fyzik´alnej podstate riaden´eho syst´emu. Ked’ˇze teplota, ktor´a
je z´akladnou procesnou veliˇcinou termodynamick´eho spr´avania budovy, je pomerne
jednodnoducho meratel’n´a, existuje mnoˇzstvo presn´
ych a cenovo dostupn´
ych senzorov, a teda je pri identifik´acii v¨aˇcˇsinou dostupn´e dostaˇcuj´
uce mnoˇzstvo nameran´
ych
d´at. Z toho dˆovodu je pouˇzitie Black box pr´ıstupu jednou z moˇznost´ı ako n´ajst’
vhodn´
y model syst´emu. Z´akladnou nev´
yhodu tejto met´ody je fakt, ˇze tento pr´ıstup
ignoruje ak´
ukol’vek ˇstrukt´
uru procesu a z´ıskan´
y model s´ıce dobre fituje najm¨a testovacie d´ata, je vˇsak vo v¨aˇcˇsine pr´ıpadov nevhodn´
y pre riadenie.
Ako najvhodnejˇsie sa spomedzi Black box met´od jav´ı pouˇzitie pomerne nov´eho
pr´ıstupu - Subspace identifikaˇcn´
ych met´od, pre ktor´e sa v literat´
ure ˇcasto pouˇz´ıva
oznaˇcenie 4SID (Subspace State-Space IDentification). Tento pr´ıstup narozdiel
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
9
od beˇzn´
ych Black box met´od (odhaduj´
ucich parametre vstupno-v´
ystupn´
ych modelov ARX,ARMAX...) poskytuje model v tvare stavov´eho popisu, ˇc´ım sa st´ava
vhodn´
ym n´astrojom aj pre identifik´aciu MIMO syst´emov s vel’k´
ym poˇctom vstuˇ sou z v´
pov a v´
ystupov (ako napr´ıklad budovy). Dalˇ
yhod, ktor´a rob´ı z 4SID met´od
vhodn´eho kandid´ata pri hl’adan´ı parametrov modelu, ktor´
y bude n´asledne pouˇzit´
y
pre predikt´ıvne riadenie, je fakt, ˇze poskytnut´
y model je optim´alnym viackrokov´
ym prediktorom. Okrem obecnej nev´
yhody Black box modelovania, ktorou je
nereˇspektovanie fyzik´alnej podstaty modelovan´eho objektu, maj´
u Subspace met´ody
aj d’al’ˇsie nev´
yhody, medzi ktor´e patr´ı potreba vel’k´eho mnoˇzstva vstupno - v´
ystupn´
ych d´at a zloˇzit´a rekurzifik´acia. Vzhl’adom na vel’k´
y objem historick´
ych nameran´
ych d´at, ktor´e pri naˇsej pr´aci boli k dispoz´ıcii, n´as uveden´e skutoˇcnosti neobmedzovali pri pouˇzit´ı Subspace met´od. Podrobnejˇsie bud´
u met´ody Subspace identifik´acie pop´ısan´e v samostatnej podkapitole.
• Grey box
Je kombin´aciou oboch vyˇsˇsie uveden´
ych pr´ıstupov, ktor´a vyuˇz´ıva znalost’ fyzik´alnej
podstaty riaden´eho syst´emu pri tvorbe ˇstrukt´
ury modelu a nameran´e d´ata k odhadu konkr´etnych hodnˆot parametrov procesu. Vzhl’adom na spom´ınan´e mnoˇzstvo
dostupn´
ych d´at a na schopnost’ zjednoduˇsene pop´ısat’ tepeln´
u bilanciu budovy je
Grey box modelovanie d’alˇsou moˇznost’ou pri konˇstukcii modelu budovy. S ohl’adom
na uveden´e vlastnosti Grey box modelovacieho pr´ıstupu mu v tejto kapitole budeme
venovat’ zv´
yˇsen´
u pozornost’ a v samostatnej podkapitole najprv pop´ıˇseme postup,
ak´
ym moˇzno n´ajst’ ˇstrukt´
uru termodynamick´eho modelu budovy a n´asledne zo sekvencie vstupno-v´
ystupn´
ych d´at urˇcit’ aj konkr´etne hodnoty parametrov dan´eho modelu.
3.1
Met´
ody Grey box identifik´
acie
´
Ulohou
Grey box modelovania je n´ajst’ fyzik´alny popis (ˇstrukt´
uru) riaden´eho syst´emu a
n´asledne odhadn´
ut’ konkr´etne hodnoty parametrov z nameran´
ych vstupno-v´
ystupn´
ych
d´at. V prvej ˇcasti tejto podkapitoly najprv pop´ıˇseme spˆosob, ak´
ym moˇzno n´ajst’ fyzik´alny popis budovy, a tento spˆosob n´asledne demonˇstrujeme na pr´ıklade vybran´eho
ˇ
bloku budovy CVUT
v Prahe. V druhej ˇcasti bude uveden´
y popis odhadu parametrov
zo vstupno-v´
ystupn´
ych d´at.
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
3.1.1
10
Tvorba ˇ
strukt´
ury modelu
Z´akladn´
ym fyzik´alnym z´akonom, z ktor´eho budeme vych´adzat’ pri modelovan´ı termodynamick´eho spr´avania budovy, je rovnica vedenia tepla [10]:
dQ
= Kie (ϑe − ϑi )
dt
⇒
dQ dϑi
= Kie (ϑe − ϑi )
dϑi dt
|{z}
(3.3)
Ci
kde t predstavuje ˇcas, ϑi , ϑi s´
u teploty v dvoch bodoch (uzloch), medzi ktor´
ymi doch´adza
k vedeniu tepla, Q je tepeln´a energia, Ci je tepeln´a kapacita v bode i. Celkov´
y koeficient
prenosu tepla Ki,e mˆoˇze byt’ spoˇc´ıtan´
y podl’a nasleduj´
uceho vzt’ahu:
1
1
1
=
+
,
Kie
Ki
Ke
(3.4)
kde Ke a Ki s´
u s´
uˇcinitele tepelnej vodivosti pre materi´aly i a e, ktor´e z´avisia predovˇsetk´
ym
na vlastnostiach dan´eho materi´alu.
Ako si moˇzno vˇsimn´
ut’, z matematick´eho hl’adiska je rovnica (3.3) ekvivalentn´a Ohˇ
movmu z´akonu. Dalej
je zrejm´e, ˇze tepeln´a kapacita C je anal´ogiou k elektrickej kapacite
a tieˇz prevr´aten´a hodnota tepelnej vodivosti K (”tepeln´
y odpor”) v termodynamickej
s´
ustave pln´ı rovnak´
u funkciu ako odpor R v elektrickom obvode. Pr´ıˇcinou prenosu tepla Q
z jedn´eho bodu s´
ustavy do in´eho je rozdiel teplˆot v t´
ychto dvoch bodoch a ten ist´
y princ´ıp
plat´ı i pre s´
ustavy elektrick´
ych prvkov a prenos elektrick´eho n´aboja. Na popis termodynamick´
ych vlastnost´ı mˆoˇzeme teda nahliadat’ podobne ako na popis veliˇc´ın v elektrickom
obvode. Pre sch´ematick´e zn´azornenie vzt’ahov medzi prvkami v elektrick´
ych s´
ustav´ach sa
vd’aka ich prehl’adnosti tradiˇcne pouˇz´ıvaj´
u tzv. RC diagramy. Ked’ˇze sme uk´azali anal´ogie
medzi pojmami zaveden´
ymi v termodynamick´
ych a elektrick´
ych s´
ustav´ach, mˆoˇzeme modifikovan´
y RC diagram bez ujmy na vernosti popisu pouˇzit’ aj na zakreslenie rel´aci´ı medzi termodynamick´
ymi prvkami [18, 19, 39]. Na nasleduj´
ucom obr´azku 3.1 uv´adzame jednoduch´
y pr´ıklad RC siete, ktor´a odpoved´a diferenci´alnej rovnici vedenia tepla 3.3.
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
11
Obr. 3.1: RC diagram pre rovnicu vedenia tepla.
Pomocou vyˇsˇsie uveden´eho postupu sme schopn´ı zostrojit’ podobn´
u RC siet’ a k nej
prisl´
uchaj´
ucu s´
ustavu diferenci´alnych rovn´ıc, ktor´a bude popisovat’ termodynamick´e spr´aˇ
vanie jedn´eho z blokov nami uvaˇzovanej budovy CVUT
s dvoma vykurovac´ımi okruhmi
a dvoma rˆozne vel’k´
ymi referenˇcn´
ymi miestnost’ami:
1
1
1
−ϑ˙ rs =
(ϑrs − ϑo ) +
(ϑrs − ϑrn ) +
(ϑrs − ϑrws )
Crs Ros
Crs Rr
Crs Rrwrs
1
1
−ϑ˙ rws =
(ϑrws − ϑrs ) +
(ϑrws − ϑhws )
Crws Rrwrs
Crws Rws
1
1
1
−ϑ˙ rn =
(ϑrn − ϑo ) +
(ϑrn − ϑrs ) +
(ϑrn − ϑrwn )
Crn Ron
Crs Rr
Crs Rrwrn
1
1
(ϑrwn − ϑrn ) +
(ϑrwn − ϑhwn )
−ϑ˙ rwn =
Crwn Rrwrn
Crwn Rwn
Obr. 3.2: RC diagram pre blok s dvoma vykurovac´ımi okruhmi.
(3.5)
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
12
V´
yznam jednotliv´
ych symbolov pouˇzit´
ych v (3.5) a na obr´azku 3.2 je vysvetlen´
y v tabul’ke 6.1. Budeme uvaˇzovat’ vstupy a stavy (odpovedaj´
u priamo meran´
ym v´
ystupom)
Tabul’ka 3.1: V´
yznam symbolov.
Oznaˇcenie Popis
Crs
tepelen´a kapacita miestnosti - juh
Crn
tepelen´a kapacita miestnosti - sever
Crws
tepelen´a kapacita vratn´a voda - juh
Crwn
tepelen´a kapacita vratn´a voda - sever
Ros
tepeln´
y odpor vonkajˇsej steny - juh
Ron
Rrwrs
tepeln´
y odpor vonkajˇsej steny - sever
tepeln´
y odpor miestnost’ juh - miestnost’ sever
tepeln´
y odpor miestnost’ - vratn´a voda (juh)
Rrwrn
tepeln´
y odpor miestnost’ - vratn´a voda (sever)
Rws
tepeln´
y odpor vykurovacia - vratn´a voda (juh)
Rwn
tepeln´
y odpor vykurovacia - vratn´a voda (sever)
Rr
tak, ako bolo uveden´e v kapitole 2:
iT
h
x = ϑrs ϑrws ϑrn ϑrwn
h
iT
u = ϑo ϑhws ϑhwn .
(3.6)
S´
ustavu diferenci´alnych rovn´ıc (3.5) mˆoˇzeme potom prep´ısat’ na nasleduj´
uci tvar, ktor´
y
odpoved´a rovnici v´
yvoja stavu v popise line´arneho ˇcasovo invariantn´eho spojit´eho syst´emu:
(3.7)
x(k)
˙
= Ac x(k) + Bc u(k),
kde

− C 1R
rs os

1
− C R
− C 1R
rs rwrs
rs r
1
Crws Rrws
1
Crn Rr


Ac = 


1
Crs Rrwrs
− C 1R
− C 1R
rws rwrs
rws rws
0
0
1
C R
 rs os

0
Bc = 

1
 Crs Ron
0
0
0
0
1
Crws Rrws
0
0
0
0
1
Crwn Rrwn
1
Crs Rr
0
0
0
1
Crn Rrwrn
− C 1R
1
− C R
− C 1R
rn on
rn rwrn
rn r
1
Crwn Rrwrn
1
1
− C
rwn Rwn
rwn Rrwrn
−C









.


(3.8)
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
13
Z rovnice (3.7) mˆoˇzeme jednoducho pomocou diskretiz´acie z´ıskat’ matice A a B v stavovom popise (3.1)
A2 T 2
A = eAT = I + Ac Ts + c s + . . . ≈ I + Ac Ts ,
2
Z Ts
Z Ts
Idτ Bc = Ts Bc ,
eAc τ dτ ≈
B=
0
0
kde Ts predstavuje vzorkovaciu

a1

a2
A=
a
 3
0
3.1.2
(3.9)
peri´odu. Matice A a B bud´
u mat’ nasleduj´
ucu ˇstrukt´
uru:



b1 0 0
a4 a3 0





a5 0 0 
 , B =  0 b3 0  .
(3.10)
b 0 0 
0 a6 a8 

 2

0 0 b4
0 a7 a9
Odhad parametrov
K tomu, aby sme z´ıskali kompletn´
y model riaden´eho syst´emu, potrebujeme eˇste odhadn´
ut’ parametre zvolenej ˇstrukt´
ury z postupnosti vstupno-v´
ystupn´
ych dat. Pokial’ budeme uvaˇzovat’, ˇze vˇsetky stavy syst´emu s´
u priamo meratel’n´e a odpovedaj´
u v´
ystupom,
mˆoˇzeme zo vstupno-v´
ystupn´
ych d´at odhadovat’ priamo matice A, B stavov´eho popisu.
Pre u
´ˇcely identifik´acie je najvhodnejˇsie rovnicu v´
yvoja stavu zo stavov´eho popisu (3.1)
zap´ısat’ do rekurz´ıvneho maticov´eho tvaru:
#
"
h
i X N −1
0
X1N = AX0N −1 + BU0N −1 + E0N −1 = A B
+ E0N −1 ,
(3.11)
N −1
U0
kde N znaˇc´ı poˇcet vzorkov a
h
Xkk+N
i
= x(k) x(k + 1) . . . x(k + N ) ,
h
i
Ukk+N = u(k) u(k + 1) . . . u(k + N ) ,
h
i
Ekk+N = e(k) e(k + 1) . . . e(k + N ) .
(3.12)
Maticov´
u rovnicu (3.25) mˆoˇzeme prep´ısat’ do tvaru:
vecX1N
=
"
X0N −1
U0N −1
#
⊗ In
!T
h
i
vec A B + vecE0N −1 ,
(3.13)
kde In predstavuje jednotkov´
u maticu rozmerov n×n, kde n je r´ad syst´emu, ⊗ predstavuje
Kroneckerov produkt dvoch mat´ıc, vec je oper´ator vektoriz´acie.
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
14
V pr´ıpade, ˇze na parametre mat´ıc A, B nie s´
u kladen´e ˇziadne obmedzenia, je moˇzn´e
ˆ B
ˆ jednoducho tak, ˇze zavedieme:
z rovnice (3.13) z´ıskat’ odhady A,
Y = vecX1N ,
#
!T
"
X0N −1
,
⊗ In
X =
U0N −1
h
i
Θ= A B ,
(3.14)
ˆ dost´avame pouˇzit´ım najmenˇs´ıch ˇstvorcov:
a pre Θ
ˆ = YX T X X T
Θ
−1
.
(3.15)
ˆ nulov´e, je nutn´e vynechat’ pr´ısluˇsn´e
V pr´ıpade, ˇze poˇzadujeme niektor´e parametre Θ
riadky regresora X aj vektora Y a odhadovat’ len parametre, ktor´e s´
u nenulov´e.
3.2
Subspace identifikaˇ
cn´
e met´
ody
Subspace identifikaˇcn´e met´ody patria k pomerne nov´
ym pr´ıstupom k identifik´acii dynamick´
ych syst´emov. Ich z´aklady boli poloˇzen´e koncom minul´eho storoˇcia na Katol´ıckej
univerzite v belgickom Leuvene a pred p´ar rokmi boli u
´speˇsne pouˇzit´e pri identifik´acii
modelu budovy [5, 11, 35]. V tejto pr´aci sa preto nebudeme t´
ymito met´odami zaoberat’
podrobne, uvedieme len ich z´akladne princ´ıpy a v d’alˇs´ıch kapitol´ach modely z´ıskan´e 4SID
met´odami porovn´ame s modelmi, ktor´e boli z´ıskane ostatn´
ymi uveden´
ymi met´odami. Pre
hlbˇsie pochopenie sa mˆoˇze ˇcitatel’ obr´atit’ napr´ıklad na [8, 22, 43, 44].
Pre u
´ˇcely Subspace identifik´acie najprv zavedieme pojem stavov´
y model v inovaˇcnom
tvare, ktor´
y z´ıskame aplik´aciou Kalmanovho filtra na syst´em so stavov´
ym popisom v obvyklom tvare (3.1):
x(k + 1) = Ax(k) + Bu(k) + K e´(k),
y(k) = Cx(k) + Du(k) + e´(k) cov {´
e(k)} = Re
(3.16)
kde Re predstavuje kovarianciu inovov´aci´ı e´(k) a matica K predstavuje ust´alen´e Kalmanovo zosilnenie [20]:
K = AP C T (CP C T + R)−1 , Re = CP C T + R
(3.17)
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
15
ktor´e z´ıskame z ust´alen´eho rieˇsenia algebraickej Ricatiho rovnice:
P = AP AT − K(CP C T + R)−1 K T + Q.
(3.18)
Dˆovodom, pre ktor´
y zav´adzame stavov´
y popis v inovaˇcnom tvare, je skutoˇcnost’, ˇze
model s popisom (3.16) m´a v stochastickej ˇcasti menej stupˇ
nov vol’nosti ako (3.1) a tie
mˆoˇzu byt’ potom jednoznaˇcne identifikovan´e zo vstupno-v´
ystupn´
ych d´at. Tie s´
u usporiadan´e pre u
´ˇcely identifik´acie do Hankelovych mat´ıc:


u0
u1 · · ·
uj−1


 u1

u
·
·
·
u
2
j


..
..
 ..

...

.
.
.
  




Up
u
u
·
·
·
u
i−1
i
i+j−2 
 =

,
 ui
ui+1 · · · ui+j−1 
Uf




ui+j 
 ui+1 ui+2 · · ·
 .

..
..
...
 .

.
.
.


ui+h−1 ui+h · · · ui+h+j−2
(3.19)
i, h, j predstavuj´
u volitel’n´e parametre - poˇcet blokov´
ych riadkov resp. st´lpcov Hankelovej
matice. Parameter i urˇcuje rozdelenie d´at na minul´e (,,past“) Up a bud´
uce (,,future“) Uf .
Analogicky s´
u skonˇstruovan´e aj matice v´
ystupov Yp , Yf a stochastick´eho vstupu E´p , E´f .
Priestory bud´
ucich a minul´
ych d´at budeme znaˇcit’ Wp a Wf :
" #
" #
Up
Uf
Wp =
, Wf =
.
Yp
Yf
(3.20)
Stavov´e premenn´e s´
u podobne ako v´
ystupy a vstupy rozdelen´e na Xf a Xp a usporiadan´e
do nasleduj´
ucich mat´ıc:
i
h
∆
Xp = x0 . . . xi−1 xi . . . xj−1
h
i
∆
Xf = xi . . . xj−1 xj . . . xi+j−1
ˇ
Dalej
zavedieme pojem rozˇs´ırenej matice pozorovatel’nosti Γk nasledovne:


C


 CA 
∆ 

Γk =  . 
 .. 


CAk−1
(3.21)
(3.22)
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
16
a obdobne pojem sp¨atnej rozˇs´ırenej matice riaditel’nosti ∆dk a ∆sk deterministick´eho resp.
stochastick´eho vstupu:
∆
∆sk =
h
k−1
i
k−2
A B A B ... B ,
h
i
∆
∆dk = Ak−1 K Ak−2 K . . . K ,
kde k > n.
(3.23)
Po zaveden´ı Toeplitzovsk´
ych mat´ıc impulzn´
ych odoziev pre vstupy u a e´, Hkd a Hks :


D
0
0
··· 0


 CB
D
0
··· 0



d ∆
CAB
CB
D
·
·
·
0
Hk = 
,

..
..
..
. . . .. 

.
.
.
.


CAk−2 B CAk−3 B CAk−4 B · · ·

D
0
0
···
D
0

(3.24)


 CK

D
0
·
·
·
0



s ∆
CAK
CK
D
·
·
·
0
Hk = 
,

..
..
..
. . . .. 

.
.
.
.


k−2
k−3
k−4
CA K CA K CA K · · · D
mˆoˇzeme prep´ısat’ stavov´
y popis (3.16) do rekurz´ıvnej maticovej formy [44]:
Yp = Γi Xp + Hid Up + His E´p ,
Yf = Γh Xf + Hhd Uf + Hhs E´f ,
(3.25)
Xp = Ai Xp + ∆di Up + His E´p .
Z´akladnou myˇslienkou Subspace identifikaˇcn´
ych met´od (ako uˇz samotn´
y n´azov napoved´a)
je vyuˇzitie podpriestorov k odhadu mat´ıc stavov´eho popisu A, B, C, D, K. Postaˇcuj´
ucimi
podpriestormi k tomuto odhadu s´
u podpriestor rozˇs´ırenej matice pozorovatel’nosti Γh a
podpriestor stavov´eho priestoru bud´
ucich d´at Xf . Na z´ıskanie t´
ychto podpriestorov je
potrebn´e odhadn´
ut’ ˇclen
O h = Γ h Xf .
(3.26)
Jednou z moˇznost´ı, ako pomerne jednoducho odhadn´
ut’ Oh z druhej rovnice v (3.16), je
pouˇzitie ˇsikmej projekcie riadkov´eho priestoru bud´
ucich v´
ystupov Yf pozd´lˇz riadkov´eho
priestoru bud´
ucich vstupov do riadkov´eho priestoru minul´
ych d´at Wp , pre ktor´
u budeme
pouˇz´ıvat’ nasleduj´
uce oznaˇcenie:
Yf / Wp .
(3.27)
Uf
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
17
S vyuˇzit´ım nasleduj´
ucich troch vlastnost´ı ˇsikmej projekcie, ktor´
ych odvodenie moˇzno
n´ajst’ v [43]:
lim Xf / Wp = Xˆf ,
i,j→∞
Uf
lim Ef / Wp = 0,
i,j→∞
(3.28)
Uf
Uf / Wp = 0,
Uf
mˆoˇzeme v pr´ıpade dostaˇcuj´
uceho poˇctu d´at p´ısat’:
ˆf ,
Yf / Wp = Γ h X
(3.29)
Uf
ˇc´ım dost´avame poˇzadovan´
y odhad matice Oh . Poˇzadovan´e podpriestory bud´
ucich stavov
a rozˇs´ırenej matice pozorovatel’nosti z´ıskame z Oh pouˇzit´ım SVD rozkladu:
h
Oh = U SV T = Us Un
"
i S
s
0
0
Sn
#"
Vs
Vn
#T
,
(3.30)
ˆ f bude platit’:
kde pre odhadovan´e Γh a X
Γh = Us Ss1/2 ,
ˆ f = Ss1/2 VsT .
X
(3.31)
Nutno podotkn´
ut’, ˇze SVD rozkladom je z´ıskan´
y nielen odhad poˇzadovan´
ych podpriestorov, ktor´e bud´
u d’alej vyuˇzit´e k odhadu parametrov stavov´eho modelu, ale aj r´ad
odhadovan´eho modelu n, ktor´
y je rovn´
y poˇctu v´
yznamn´
ych singularn´
ych ˇc´ısel matice Oh ,
ktor´e tvoria diagon´alu matice Ss .
Teraz uˇz pozn´ame postupnost’ stavov Xf a k odhadu mat´ıc A, B, C, D je potrebn´e uˇz
len vyrieˇsit’ nasleduj´
ucu s´
ustavu rovn´ıc:
"
# "
#" #
ˆ i+1
ˆi
X
A B X
=
,
(3.32)
Yi
C D Ui
kde
h
i
h
ˆ i+1 = xˆi+1 · · · xˆi+j , X
ˆ i = xˆi · · ·
X
i
h
h
Yi = yi · · · yi+j−1 , Ui = ui · · ·
i
xˆi+j−1 ,
i
ui+j−1 .
(3.33)
Rieˇsenie s´
ustavy rovn´ıc (3.32) z´ıskame pomerne jednoducho pomocou najmenˇs´ıch ˇstvorcov.
Zavedieme podobn´e oznaˇcenie ako v predoˇslej podkapitole:
"
"
#
#
i
h
ˆ i+1
X
A B
ˆ i , Ui , Y =
Θ=
, X = X
.
Yi
C D
(3.34)
ˇ E
´ METODY
´
KAPITOLA 3. IDENTIFIKACN
ˆ B,
ˆ C,
ˆ D
ˆ pomocou LS dost´avame
Pre odhad syst´emov´
ych mat´ıc A,
"
#
ˆ
−1
Aˆ B
.
= YX T X X T
ˆ
Cˆ D
18
(3.35)
ˆ a Rˆe spoˇc´ıtame
Pre odhad stochastickej ˇcasti modelu so stavov´
ym popisom (3.16) K
rezidu´a Σ:
Σ=
ˆ
kde Y:
"
Σ11 Σ12
Σ21 Σ22
#
=
T
1
Y − Yˆ Y − Yˆ ,
j − (m + n)(n + l)
"
#
ˆ B
ˆ
A
Yˆ =
X.
ˆ
Cˆ D
(3.36)
(3.37)
ˆ potom dost´avame:
Pre odhad Rˆe a K
ˆ e = Σ22 ,
R
ˆ = Σ11 Σ−1
K
22 .
(3.38)
Kapitola 4
Identifik´
acia pre predikt´ıvne riadenie
Ako uˇz bolo neraz spomenut´e, pre spr´avne fungovanie predikt´ıvneho regul´atora je vel’mi
dˆoleˇzit´a dobr´a znalost’ matematick´eho modelu riaden´eho syst´emu. Krit´erium (´
uˇcelov´a
funkcia), ktor´e bude pouˇzit´e pri identifik´acii parametrov modelu, by malo byt’ volen´e
podl’a toho, pre ak´
y typ regul´atoru bude tento model n´asledne pouˇzit´
y. Predikt´ıvny regul´ator minimalizuje regulaˇcn´
u odch´
ylku poˇcas cel´eho predikˇcn´eho horizontu v z´avislosti
na predikcii v´
ystupov. Najˇcastejˇsie je pri n´avrhu MPC pouˇzit´a kvadratick´a u
´ˇcelov´a funkcia a regul´ator potom minimalizuje kvadr´at regulaˇcnej odch´
ylky. Model, ktor´
y bude
pouˇzit´
y pre predikt´ıvne riadenie, by mal byt’ teda predovˇsetk´
ym dobr´
ym viackrokov´
ym
prediktorom.
Beˇzne pouˇz´ıvan´e identifikaˇcn´e met´ody minimalizuj´
u vˇsak len jednokrokov´
u predikˇcn´
u
chybu (PEM - Prediction Error Methods, [29]). Preto sa pri tvorbe modelov pre MPC
(najm¨a pre vel’k´e predikˇcn´e horizonty) pouˇz´ıvaj´
u met´ody, ktor´e minimalizuj´
u viackrokov´
u predikˇcn´
u chybu. Pre tieto met´ody sa pouˇz´ıva s´
uhrnn´e oznaˇcenie MPC Relevant
Identification (MRI) a bud´
u podrobnejˇsie pop´ısan´e v nasleduj´
ucej podkapitole.
4.1
´
Uvod
do MRI
V posledn´
ych rokoch bolo publikovan´e pomerne vel’k´e mnoˇzstvo ˇcl´ankov, ktor´e sa zaoberaj´
u hl’adan´ım modelu, ktor´
y by mohol byt’ pouˇzit´
y pre predikt´ıvne riadenie.
Prv´
ymi pr´acami v tejto oblasti s´
u ˇcl´anky autorov Shook a Mohtadi publikovan´e
na zaˇciatku dev¨at’desiatych rokov [40, 41]. V t´
ychto ˇcl´ankoch poukazuj´
u na to, ˇze minimaliz´acia viackrokovej predikˇcnej chyby je ekvivalentn´a prefiltrovaniu vstupn´
ych a
19
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
20
v´
ystupn´
ych d´at filtrom, ktor´
y je z´avisl´
y na ˇsumovom modeli, a n´aslednej identifik´acii
minimalizuj´
ucej jednokrokov´
u predikˇcn´
u chybu. Pouˇzitie takejto met´ody je vˇsak podmienen´e znalost’ou skutoˇcn´eho ˇsumov´eho modelu, ktor´
y b´
yva vo v¨aˇcˇsine priemyseln´
ych
aplik´acii nezn´amy.
V rokoch 2002 aˇz 2004 vych´adzaj´
u ˇcl´anky [16,17], kde sa autori s neznalost’ou ˇsumov´eho
modelu vysporiadali rozdelen´ım MRI identifik´acie do dvoch krokov. V prvom kroku identifikuj´
u model deterministickej ˇcasti s´
ystemu, druh´
y krok algoritmu predstavuje odhad
koeficientov ˇsumov´eho modelu s vyuˇzit´ım znalosti o deterministickej ˇcasti syst´emu.
ˇ s´ım z moˇzn´
Dalˇ
ych pr´ıstupov k MRI identifik´acii je ,,multimodelov´
y pr´ıstup“, kde je
vytvoren´
y separ´atny model pre kaˇzd´e k ∈ {1, 2, 3, . . . , P }, P je horizont predikcie. Vˇsetky
z´ıskan´e modely s´
u s´
uˇcasne pouˇzit´e pre riadenie. Z´akladnou nev´
yhodou tohto pr´ıstupu je
vel’k´e mnoˇzstvo odhadovan´
ych parametrov najm¨a pre vel’k´e predikˇcn´e horizonty a syst´emy
s viacer´
ymi vstupmi a v´
ystupmi (MIMO). Ked’ˇze rozptyl odhadovan´
ych parametrov sa
zvyˇsuje s rast´
ucim poˇctom parametrov, touto met´odou sa nebudeme d’alej zaoberat’. Jej
podrobnejˇs´ı popis je moˇzn´e n´ajst’ v [38].
4.2
MRI identifik´
acia
4.2.1
MRI u
´ˇ
celov´
a funkcia
Je zn´ame, ˇze predikt´ıvny regul´ator minimalizuje odch´
ylku skutoˇcn´eho v´
ystupu syst´emu
od poˇzadovanej referenˇcnej hodnoty poˇcas cel´eho predikˇcn´eho horizontu. Vo v¨aˇcˇsine
pr´ıpadov je pouˇzit´e krit´erium kvadratick´e.
Definujme nasleduj´
ucu u
´ˇcelov´
u funkciu JM P C , ktor´a penalizuje kvadr´at odch´
ylky
od referencie:
JM P C
P
N
−P X
X
1
=
(yref (k + i) − y(k + i))2 .
(N − P ) P k=1 i=1
(4.1)
Ak pre v´
ystup pouˇzijeme y(k + i|k) = yˆ(k + i|k) + e(k + i|k), mˆoˇzeme krit´erium (4.1)
nap´ısat’ v tvare
JM P C
P
N
−P X
X
1
(yref (k + i) − yˆ(k + i|k) − (y(k + i) − yˆ(k + i|k)))2 . (4.2)
=
(N − P ) P k=1 i=1
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
21
ˇ s´ım rozpisovan´ım mˆoˇzeme vyjadrit’ krit´erium JM P C v tvare s´
Dalˇ
uˇctu troch v´
yrazov:
JM P C
P
N
−P X
X
1
(yref (k + i) − yˆ(k + i|k))2 +
=
(N − P ) P k=1 i=1
P
N
−P X
X
1
(y(k + i) − yˆ(k + i|k))2 −
(N − P ) P k=1 i=1
(4.3)
P
N
−P X
X
2
(yref (k + i) − yˆ(k + i|k)) × (y(k + i) − yˆ(k + i|k))
(N − P ) P k=1 i=1
Klasick´e MPC minimalizuje len prv´
y ˇclen rovnice (4.3) a to mˆoˇze viest’ na suboptim´alne rieˇsenie. Aby sme dospeli k optim´alnemu rieˇseniu, je potrebn´e minimalizovat’
aj dva zost´avaj´
uce ˇcleny. Tret´ı ˇclen predstavuje kr´ıˇzov´
u korel´aciu medzi chybou identifik´acie a chybou riadenia a jeho minimaliz´aciou sa zaober´a napr´ıklad [15] a v tejto pr´aci
mu nebude d’alej venovan´a pozornost’.
Mimimaliz´acia druh´eho ˇclenu je zabezpeˇcen´a pouˇzit´ım MRI identifikaˇcn´
ych algoritmov. Zaved’me si pre druh´
y ˇclen rovnice (4.3) nasleduj´
uce oznaˇcenie:
JNP =
P
P
N
−P X
N
−P X
X
X
1
1
(y(k + i) − yˆ(k + i|k))2 =
(ǫp (k))2 .
(N − P ) P k=1 i=1
(N − P ) P k=1 i=1 i
(4.4)
Na oznaˇcenie chyby predikcie na i krokov v ˇcase k, ǫpi (k), budeme d’alej pouˇz´ıvat’ znaˇcenie
ǫpk,i . V d’alˇsom texte sa budeme zaoberat’ uˇz len minimaliz´aciou vyˇsˇsie uvedenej u
´ˇcelovej
funkcie.
Dodajme, ˇze pre i-krokov´
u predikciu v´
ystupu budeme pouˇz´ıvat’ oznaˇcenie yˆ(k + i|k),
aby bola zachovan´a s´
uvislost’ medzi viackrokovou predikˇcnou chybou a u
´ˇcelovou funkciou
MPC 4.1.
Uvaˇzujme skutoˇcn´
y model syst´emu v tvare:
y(k) = Gp (q)u(k) + GL (q)e(k).
(4.5)
Gp (q) predstavuje prenos deterministickej ˇcasti syst´emu a GL (q) je tvarovac´ı filter ˇsumu.
Predpoklad´a sa, ˇze e(k) je norm´alne rozdelen´
y biely ˇsum s nulovou strednou hodnotou a
ˆ p (q, θ)
rozptylom σe2 . V t´
ychto v´
yrazoch q oznaˇcuje oper´ator posunu y(k + i) = qy(k). G
ˆ L (q, θ) predstavuj´
aG
u identifikovan´e modely, θ je nezn´amy odhadovan´
y parameter.
Po zaveden´ı tohoto oznaˇcenia mˆoˇzeme pre optim´alny viackrokov´
y prediktor yˆ(k + i|k)
p´ısat’ [29]:
ˆ iG
ˆ p (q)u(k + i) + 1 − W
ˆ i y(k + i).
yˆ(k + i|k) = W
(4.6)
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
22
ˆ p (q) p´ısat’ iba G
ˆ p . Pre W
ˆ i zo vzt’ahu (4.6)
Pre jednoduchost’ budeme d’alej namiesto G
plat´ı
ˆ i = Fˆi G
ˆ −1 ,
W
L
Fˆi =
i−1
X
gˆi q −1 ,
Fˆ1 = gˆ0 = 1,
(4.7)
i=0
kde koeficienty gˆi s´
u koeficienty impulznej odozvy odhadovan´eho tvarovacieho filtru ˇsumu.
V´
yrazom yˆ(k + i|k) rozumieme predikovan´e v´
ystupy v ˇcase k + i z d´at, ktor´e s´
u dostupn´e
v ˇcase k.
S pouˇzit´ım vzt’ahu pre optim´alny prediktor (4.6) mˆoˇzeme vyjadrit’ chybu predikcie ǫpi,k
(4.4) nasleduj´
ucim vzt’ahom:
ˆ iG
ˆ p (q)u(k + i)
ǫpk,i = Gp (q)u(k + i) + GL (q)e(k + i) − W
ˆ
ˆ
ˆ
− 1 − Wi Wi Gp (q)u(k + i) + GL e(k + j) .
Po d’aˇs´ıch u
´prav´ach dostaneme pre ǫpi,k :
p
ˆ
ˆ
ˆ i GL e(k + i).
ǫi,k = Wi Gp − Gp u(k + i) + W
(4.8)
(4.9)
Za predpokladu nekorelovanosti vstupov a ˇsumu dost´avame pre podmienen´
u stredn´
u
hodnotu µ v´
yraz:
ˆ i Gp − G
ˆ p u(k + i) = µi (k).
E(ǫpi,k |u) = W
(4.10)
Chyba predikcie je norm´alne rozdelen´a s urˇcitou strednou hodnotou µ a rozptylom σ 2 .
Stredn´a hodnota chyby ǫpi,k je z´avisl´a jednak na kvalite odhadu deterministickej ˇcasti
system´
u a aj na vlastnostiach ˇsumu.
Pre rozptyl chyby predikcie σi2 mˆoˇzeme nap´ısat’
2
2
ˆ
σi2 = W
G
i L σe .
(4.11)
2
Ako je vidno, rozptyl σi2 chyby predikcie je z´avisl´
y len vlastnostiach ˇsumu. S vyuˇzit´ım
(4.11) a (4.10) mˆoˇzeme pre MRI u
´ˇcelov´
u funkciu pouˇzit’ vzt’ah:
JNP
P
P
1 X 2 X 2
σ +
µi )
= (
P i=1 i
i=1
(4.12)
ˆ L ) sa probl´em minimaliz´acie viackrokoV pr´ıpade znalosti ˇsumov´eho modelu (GL = G
vej predikˇcnej chyby ǫpi,k redukuje len na minimaliz´aciu druh´eho ˇclenu (4.12). Prv´
y ˇclen
(4.12), ktor´
y predstavuje rozptyl chyby predikcie σ, je v tomto pr´ıpade zn´amy a nem´a
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
23
ˆ i pouˇzijeme
zmysel ho minimalizovat’. Jeho hodnotu mˆoˇzeme urˇcit’ nasledovne: pre W
vzt’ah (4.7), dosad´ıme do rovnice (4.11) a dost´avame
P
X
σi2 = P + (P − 1)g12 + · · · + gP2 −1 σe2 .
i=1
(4.13)
Jednou z moˇznost´ı minimaliz´acie (4.12) je vyuˇzit’ ekvivalenciu prefiltr´acie vstupn´
ych
a v´
ystupn´
ych d´at a minimaliz´acie viackrokovej predikˇcnej chyby, podobne ako je uveden´e
v [41] a [40]. Pokial’ vˇsak nepozn´ame vlastnosti ˇsumov´eho modelu, ˇco je v praktick´
ych
aplik´aciach vel’mi ˇcast´
y jav, je nutn´e mysliet’ aj na minimaliz´aciu rozptylu chyby predikcie
σi2 . Ako bolo spomenut´e v u
´vode tejto kapitoly, tento probl´em rieˇsia autori v ˇcl´ankoch
[16,17]. Podrobnejˇsi popis algoritmu, ktor´
y je oznaˇcovan´
y ako dvojkrokov´
y, bude uveden´
y
v nasleduj´
ucej podkapitole.
4.2.2
Dvojkrokov´
y algoritmus
Minimaliz´aciu MRI u
´ˇcelovej funkcie je moˇzn´e rozdelit’ do dvoch nez´avisl´
ych krokov:
ˆ L = 1),
1. Zafixovanie ˇ
sumov´
eho modelu - zafixujeme ˇsumov´
y model (napr. G
pre dan´
y ˇsumov´
y model z´ıskame model deterministickej ˇcasti syst´emu Gp .
2. Odhad modelu ˇ
sumu - s pouˇzit´ım modelu procesu, ktor´
y sme z´ıskali v prvom
ˆ L.
kroku algoritmu, odhadneme koeficienty tvarovacieho filtru ˇsumu G
Zafixovanie ˇ
sumov´
eho modelu
V prvom kroku dvojkrokov´eho algoritmu zafixujeme ˇsumov´
y model a hl’ad´ame model
ˆ p tak, aby bol minimalizovan´
procesu G
y druh´
y ˇclen krit´eria (4.12) pre zvolen´
y ˇsumov´
y
model. Rieˇsime teda nasleduj´
ucu optimalizaˇcn´
uu
´lohu
ˆ p = argmin
G
Gp
N
−P
X
P
X
µ2i
(4.14)
k=1 i=1
ˆ L = 1. Je
Ako najjednoduchˇsia vol’ba sa ukazuje zvolit’ ako fixn´
y ˇsumov´
y model G
ˆ L = 1 je MRI identifik´acia ekvivalentn´a klasickej
jednoduch´e uk´azat’, ˇze za predpokladu G
identifik´acii minimalizuj´
ucej jednokrokov´
u predikˇcn´
u chybu.
ˆ L = 1, potom pre koeficienty ipmulznej odozvy a filter W
ˆ i plat´ı
Ak si zvol´ıme G
gˆ0 = 1,
Fˆi = 1,
ˆ i = 1,
W
gˆi = 0, i = {1, . . . , P } ,
(4.15)
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
24
ˆ i = 1 vid´ıme, ˇze minimaliz´acia (4.9) pre l’ubovol’n´e P je ekvivaletn´a
Po dosaden´ı za W
minimaliz´acii jednokrokovej predikˇcnej chyby (P = 1).
K´
ym model procesu, ktor´
y vyhovuje minimaliz´acii jednokrokovej predikˇcnej chyby, je
moˇzn´e z´ıskat’ celkom jednoducho pouˇzit´ım line´arnej regresie, pri hl’adan´ı minima sumy
kvadr´atov viackrokov´
ych predikˇcn´
ych ch´
yb (4.9) pre P > 1 musia byt’ pouˇzit´e algoritmy
neline´arnej optimaliz´acie [49], ˇco je v´
ypoˇctovo omnoho n´aroˇcnejˇsie.
Odhad modelu ˇ
sumu
ˆ p , ktor´
Odhadnut´
y model procesu G
y sme z´ıskali v prvom kroku algoritmu, bude d’alej
pouˇzit´
y pre odhadovanie koeficientov tvarovacieho filtru ˇsumu GL .
Budeme uvaˇzovat’ model ˇsumu:
−1
−m
ˆ L = C(q) = 1 + c1 q · · · + cm q
G
,
D(q)
1 + d1 q −1 · · · + dn q −n
(4.16)
kde m ud´ava r´ad ˇcitatel’a modelu ˇsumu a n r´ad menovatel’a. Pripomeˇ
nme, ˇze koeficienty
impulznej odozvy ˇsumov´eho modelu gˆi s´
u funkciou koeficientov ci a di a moˇzno ich z´ıskat’
dlh´
ym delen´ım polyn´omov C/D.
Pre viackrokov´
u predikˇcn´
u chybu ǫk,i (θi ) mˆoˇzeme s pouˇzit´ım (4.16) nap´ısat’
ǫk,i (θi ) = y(k + i) − yˆ(k + i|θi )
i
D h
ˆ p u(k + i)
Fi y(k + i) − Fi G
ǫk,i (θi ) =
C
kde θi predstavuje vektor nezn´amych odhadovan´
ych parametrov:
h
i
θi = c 1 · · · c m 1 d 1 · · · d n
(4.17)
(4.18)
Objasnime eˇste v´
yznam doln´eho indexu i pri odhadovanom parametri θ: ako uˇz bolo
spomenut´e, koeficienty impulznej odozvy s´
u z´avisl´e na koeficientoch prenosu ci a di
Fˆi = f (c1 · · · ci , d1 · · · di ).
(4.19)
ˇ
Potom θi predstavuje vektor odhadovan´
ych parametrov, pre ktor´e plat´ı (4.19). Dalej
si
zavedieme nasleduj´
uce oznaˇcenie:
ˆ p u(k),
uf (k) = G
(4.20)
ˆ p u(k + i).
vk,i (θi ) = Fˆi y(k + i) − Fˆi G
(4.21)
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
25
Potom mˆoˇzeme vk,i (θi ) rozp´ısat’ do tvaru:
vk,i (θi ) = y(k + i) + gˆ1 y(k + i − 1) + · · · + gˆi−1 y(k + 1)
(4.22)
− uf (k + i) − gˆ1 uf (k + i − 1) − · · · − gˆi−1 uf (k + 1).
S pouˇzit´ım (4.22) mˆoˇzeme viackrokov´
u predikˇcn´
u chybu ǫk,i (θi ) vyjadrit’ ako
h
ǫk,i (θi ) = θi vk,i (θi ) vk−1,i (θi ) · · ·
vk−n,i (θi ) ǫk−1,i (θi ) · · ·
ǫk−m,i (θi )
iT
.
(4.23)
Teraz uˇz mˆoˇzeme druh´
y krok algoritmu formulovat’ ako optimalizaˇcn´
uu
´lohu
arg min
θ
N
−P
X
P
X
(ǫk,i )2
(4.24)
k=1 i=1
priˇcom pre ∀ i ∈ {1, . . . , P } mus´ı byt’ splnen´a obmedzuj´
uca podmienka (4.19).
Ide o klasick´
y probl´em neline´arnej optimaliz´acie, ktor´
y je moˇzn´e rieˇsit’ pouˇzit´ım neline´arnych najmenˇs´ıch ˇstvorcov. Na rieˇsenie tohto probl´emu je moˇzn´e pouˇzit’ napr´ıklad
solver lsqnonlin, ktor´
y je s´
uˇcast’ou Optimalizaˇcn´eho toolboxu v Matlabe [6].
Pr´ıklad 4.1 (MRI dvojkrokov´
a identifik´
acia): Uvaˇzujme model procesu ktor´
y je
dan´
y skutoˇcn´
ym modelom procesu Gp
Gp =
0.0077z 2 + 0.0212z + 0.0036
z 3 − 1.9031z 2 + 1.1514z − 0.2158
(4.25)
a modelom ˇsumu GL1 , GL2
GL1 =
1
,
z + 0.8
GL2 =
1
,
(z − 0.8)2
(4.26)
ktor´
y je buden´
y pseudon´ahodn´
ym bin´arnym sign´alom s amplit´
udou ±3 o d´lˇzke 500
vzorkov a zat’aˇzen´
y gaussovsk´
ym bielym ˇsumom s nulovou strednou hodnotou a rozptylom σe2 = 0.3. V prvom pr´ıpade GL1 predpokladajme, ˇze pozn´ame skutoˇcn´
u ˇstrukt´
uru
ˇsumov´eho modelu
1
(4.27)
z+d
a v pr´ıpade ˇsumov´eho modelu GL2 predpoklad´ame, ˇze nepozn´ame ˇstrukt´
uru ˇsumov´eho
ˆ L1 =
G
modelu
ˆ L2 =
G
1
z + db
(4.28)
Rieˇsenie: a, GL1 (zn´
ama ˇ
strukt´
ura ˇ
sumov´
eho modelu)
Na odhad koeficientu ˇsumov´eho modelu bol pouˇzit´
y dvojkrokov´
y algoritmus. V prvom
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
26
ˆ L = 1 a boli odhadnut´e koeficienty deterministickroku bol zafixovan´
y ˇsumov´
y model na G
kej ˇcasti. Tie boli v druhom kroku pouˇzit´e k odhadu nezn´ameho koeficientu d tvarovacieho
ˆ L . Na rieˇsenie optimalizaˇcn´eho probl´emu bol pouˇzit´
filtru ˇsumu G
y solver lsqnonlin.
Na obr. 5.1 je vykreslen´a optim´alna hodnota p´olu ˇsumov´eho modelu d v z´avislosti
na predikˇcnom horizonte P ∈ {1, 5, 10, . . . , 80}. Ukazuje sa, ˇze so zvyˇsuj´
ucim sa predikˇcn´
ym horizontom rastie aj kvalita odhadu a pri spr´avne zvolenej ˇstrukt´
ure ˇsumov´eho
modelu je dvojkrokov´
y algoritmus schopn´
y odhadn´
ut’ ˇsumov´
y model omnoho presnejˇsie
ako klasick´e jednokrokov´e identifikaˇcn´e met´ody. Na druhej strane pre vel’k´e predikˇcn´e
horizonty P sa odhadovan´
y parameter ust´ali na urˇcitej hodnote a d’alej sa men´ı len m´alo
a vplyv zvyˇsovania predikˇcn´eho horizontu na kvalitu odhadu je minim´alny.
To si moˇzno vysvetlit’ nasledovne: pri pohl’ade na vzt’ahy (4.22),(4.23) vid´ıme ˇze rozdiel medzi minimaliz´aciou u
´ˇcelovej funckie (4.24) pre rˆozne P je dan´
y koeficientami impulznej odzvy tvarovacieho filtra ˇsumu GL . Ked’ˇze tvarovac´ı filter ˇsumu m´a impulzn´
u
odozvu koneˇcnej d´lˇzky, koeficienty gi s´
u od urˇcit´eho i (v naˇsom pr´ıpade pre tvarovac´ı filter GL1 (4.28) cca 30) zanedbatel’ne mal´e a minimaliz´acia viackrokovej predikˇcnej chyby
je pre P > 30 nez´avisl´a na d’alˇsom zvyˇsovan´ı P .
0.8
0.7
0.6
odhadnuté d
skutočné d
0.5
d
0.4
0.3
0.2
0.1
0
−0.1
−0.2
10
20
30
40
P
50
60
70
80
Obr. 4.1: Hodnota optim´alneho koeficientu ˇsumov´eho modelu v z´avislosti
na P .
Ako mer´ıtko kvality odhadu sme pouˇzili pomer hodnˆot u
´ˇcelov´
ych funkci´ı:
PN −80 P80
ˆ 2
i=1 (ǫk,i (θP ))
.
RP = Pk=1
P
80
N −80
ˆ 2
i=1 (ǫk,i (θ1 ))
k=1
(4.29)
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
27
kde θˆP , P ∈ 1, 2, . . . , 80 predstavuje odhad nezn´amych parametrov z´ıskan´
y minimaliz´aciou P -krokovej predikˇcnej chyby. Na obr. 4.2, kde je zobrazen´
y priebeh RP v z´avislosti
na zvyˇsuj´
ucom sa predikˇcnom horizonte, je moˇzn´e pozorovat’ markantn´
y n´arast kvality
odhadu d so zvyˇsovan´ım P aˇz po urˇcit´
y predikˇcn´
y horizont, kde sa RP ust´ali na urˇcitej
hodnote a jeho zmeny s´
u zanedbatel’n´e.
1
0.98
0.96
RP
0.94
0.92
0.9
0.88
0.86
0.84
10
20
30
40
P
50
60
70
80
Obr. 4.2: Priebehy RP v z´avislosti na P .
b, GL2 (nezn´
ama ˇ
strukt´
ura ˇ
sumov´
eho modelu)
V druhom pr´ıpade zloˇzitejˇsieho ˇsumov´eho modelu GL2 sme podobne ako v prvom pr´ıpade
ˆ L2 = 1 a n´asledne v druhom kroku odhadli koefiodhadli koeficienty modelu procesu pre G
cient db . To, ako sa men´ı hodnota odhadnut´eho parametru db so zvyˇsuj´
ucim sa predikˇcn´
ym
horizontom P , je moˇzn´e vidiet’ na obr. 4.3. Ako mer´ıtko kvality odhadu bola tentokr´at
ˆ L2 .
pouˇzit´a Nyquistova frekvenˇcn´a charaketristika [13] odhadnut´eho ˇsumov´eho filtra G
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
28
−0.84
−0.86
−0.88
db
−0.9
−0.92
−0.94
−0.96
−0.98
10
20
30
40
P
50
60
70
80
Obr. 4.3: Hodnota optim´alneho koeficientu ˇsumov´eho modelu v z´avislosti
na P .
Na frekvenˇcn´
ych charakteristik´ach odhadnut´
ych ˇsumov´
ych modelov obr. 4.4 je vidno,
ˇze najbliˇzˇsie ku skutoˇcn´emu ˇsumov´emu modelu maj´
u modely z´ıskane pre predikˇcn´e horizonty okolo 10, d’alˇs´ım zvyˇsovan´ım predikˇcn´eho horizontu P sa kvalita odhadu nijak
nezlepˇsuje, naopak kles´a.
0
−5
Imaginary Axis
−10
Rastúce P
−15
Skutočný šumový model
−20
−25
P = 1
−30
0
10
20
30
Real Axis
40
50
60
Obr. 4.4: Nyquistove charakteristiky pre rˆozne P .
70
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
4.3
29
Jednokrokov´
y algoritmus
In´
y pr´ıstup k MRI identifik´acii prezentuj´
u v ˇcl´anku [25] autori Lauri a Mart´ınez a bude
mu venovan´a t´ato podkapitola.
Pre v´
ystup obecn´eho MIMO syst´emu mˆoˇzeme pouˇzit’ vzt’ah:
(4.30)
y(k) = Z(k)Θ + E(k),
kde
h
• y(k) = y1 (k) · · ·
i
yno (k) ,
h
• Z(k) = u(k − nk ) · · ·
i
u(k − nb ) y(k − 1) y(k − na ) .
Matica Z(k) je regresor, no znaˇc´ı poˇcet v´
ystupov, nb ud´ava poˇcet oneskoren´
ych
v´
ystupov v regresore (lagged outputs), na poˇcet oneskoren´
ych vstupov (lagged inputs)
v regresore, nk ud´ava oneskorenie prv´eho vstupu v regresore za aktu´alnym v´
ystupom
(nk = 0 znamen´a priamu v¨azbu vstup-v´
ystup) a Θ je vektor odhadovan´
ych parametrov
syst´emu
h
Θ = bn k · · ·
b n b a1 · · ·
an a
iT
.
(4.31)
Po zaveden´ı tejto konvencie mˆoˇzeme pre viackrokov´
y prediktor yˆ(k + i|k) nap´ısat’
ˆ T,
yˆ(k + i|k) = Z(k + i)Θ
i ∈ 1, 2, . . . , P ,
(4.32)
i
u(k + i − nb ) y(k + i − 1) y(k + i − na ) .
(4.33)
kde
h
Z(k + i) = u(k + i − nk ) · · ·
Treba si uvedomit’, ˇze kaˇzd´
y v´
ystup y(a) v regresnom vektore (4.33), kde a > k, nie je
d´atom, ktor´e by bolo dostupn´e v aktu´alnom ˇcase k, ale ide o predikciu v´
ystupu yˆ(a|k). Aby
sme z´ıskali predikciu yˆ(k + i|k), pouˇzijeme i-kr´at vzt’ah (4.32), t´
uto rekurz´ıvnu proced´
uru
odˇstartujeme z aktu´alneho v´
ystupu y(k).
Ako sme uˇz spomenuli v u
´vode podkapitoly 4.3, u
´lohou MRI identifik´acie je n´ajst’
nezn´ame parametre identifikovan´eho syst´emu tak, aby bol minimalizovan´
y druh´
y ˇclen
rovnosti (4.3), ktor´
y predstavuje sumu kvadr´atov viackrokov´
ych predikˇcn´
ych ch´
yb ǫk,i =
ǫ(k + i|k):
ǫk,i = y(k + i) − yˆ(k + i).
(4.34)
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
30
´ celov´
Uˇ
u funkciu, ktor´
u budeme znaˇcit’ podobne ako v [25] JLRP I (Long Range Predictive
Identification) a predstavuje spom´ınan´
y druh´
y ˇclen (4.3), mˆoˇzeme nap´ısat’ v nasleduj´
ucom
tvare:
JLRP I =
P X
N −i
X
(ǫ(k + i|k))2 .
(4.35)
i=1 k=1
K´
ym pre P = 1 ide o minimaliz´aciu klasickej jednokrokovej predikˇcnej chyby, pre P > 1
kvˆoli spom´ınanej rekurzii (4.32),(4.33) probl´em minimaliz´acie JLRP I op¨at’ predstavuje
probl´em neline´arnej optimaliz´acie a nemoˇzno pouˇzit’ line´arnu regresiu.
4.3.1
Implement´
acia (popis) algoritmu
Krit´erium 4.35 mˆoˇzeme nap´ısat’ aj v maticovej forme
JLRP I =
P X
N −i
X
(Ea )2 ,
(4.36)
i=1 k=1
kde

E a1
 . 
. 
Ea = 
 . ,
E aP



E ai = 

ǫ(1 + i|1)
..
.
ǫ((N − i) + i|N − i)


,

i ∈ {1, 2, . . . , P } .
(4.37)
Matica Eai predstavuje maticu i-krokov´
ych predikˇcn´
ych ch´
yb v kaˇzdom dikr´etnom ˇcase
k ∈ {1, 2, . . . , N − i}. Ide o vysok´
u maticu, ktor´a m´a poˇcet st´lpcov rovn´
y poˇctu v´
ystupov
syst´emu no . Poˇcet riadkov je dan´
y d´lˇzkou identifikaˇcn´
ych d´at. Ked’ˇze ide o i-krokov´e
predikˇcn´e chyby, pre kaˇzd´e i existuje pr´ave N − i predikˇcn´
ych ch´
yb. Pre Eai plat´ı:
Eai = Yai − Zai Θ,
(4.38)
kde Yai a Zai s´
u matice v´
ystupov resp. regresorov, ktor´e boli vytvoren´e podobn´
ym
spˆosobom ako matice predikˇcn´
ych ch´
yb:




X(1 + i)
Y1+i


 . 
..
.
..  , Zai = 
Yai = 
.




Z((N − i) + i)
YN
(4.39)
Rovnako, ako s´
u usporiadan´e matice vo vel’kej matici Ea , usporiadame aj matice
regresorov a v´
ystupov. Pre predikˇcn´
u chybu Ea teda dostaneme:
E a = Y a − Za Θ
(4.40)
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
31
a formulujeme optimalizaˇcn´
uu
´lohu:
arg min(Ya − Za (Θ)Θ)2 .
(4.41)
Θ
Takto formulovan´
y optimalizaˇcn´
y probl´em je op¨at’ moˇzn´e rieˇsit’ pouˇzit´ım solveru
pre neline´arne najmenˇsie ˇstvcorce lsqnonlin. Na rieˇsenie tejto ˇspecifickej u
´lohy je vˇsak
vhodnejˇs´ı solver lsqcurvef it hl’adaj´
uci koeficienty x, ktor´e rieˇsia probl´em:
arg min kF (x, xdata ) − ydata k22 .
(4.42)
x
Pri pohl’ade na (4.41) vid´ıme, ˇze ide o optimalizaˇcn´
uu
´lohu (4.42), kde nezn´amy parameter x je Θ, v´
ystupy Ya predstavuj´
u ydata a F je v naˇsom pr´ıpade s´
uˇcin Xa (Θ)Θ.
Pomocou pop´ısan´eho algoritmu (s pouˇzit´ım solveru lsqcurvef it) sme identifikovali
parametre MIMO syst´emu a v´
ysledku s´
u uveden´e v nasleduj´
ucom pr´ıklade.
Pr´ıklad 4.2 (MRI identifik´
acia MIMO syst´
emu): Uvaˇzujme model syst´emu destilaˇcn´eho st´lpca [47], ktor´
y je pop´ısan´
y maticou prenosu Gp pre vzorkovaciu peri´odu Ts =
60 s:
Gp =
"
0.06z −1 +0.63z −2 +0.072z −3
1−0.93z −1 −0.013z −2 −0.0023z 3
0.2z −1 −0.63z −2 +0.52z −3
1−2.18z −1 +1.58z −2 −0.39z 3
0.033z −1 +0.087z −2 −0.58z −3
1−1.57z −1 +0.72z −2 −0.13z −3
0.051z −1 +0.12z −2 −0.85z −3
1−1.55z −1 +0.71∗z −2 −0.13z 3
#
.
(4.43)
Vstupmi syst´emu s´
u prietoky refluxu (sp¨atn´
y tok) a var´aku. V´
ystupy predstavuj´
u
zloˇzenia produktov vystupuj´
ucich z horn´
ych a doln´
ych otvorov kol´ony. Pre u
´ˇcely identifik´acie sme budili syst´em pseudon´ahodn´
ym bin´arnym sign´alom s amplit´
udou ±2 a bol
zat’aˇzen´
y gaussovsk´
ym bielym ˇsumom s nulovou strednou hodnotou a rozptylom σe2 = 0.4.
Nezn´ame parametre syst´emu sme identifikovali z testovac´ıch d´at o d´lˇzke N = 500 a validovali na 400 vzorkoch. Pouˇzit´
y bol tvarovac´ı filter ˇsumu GL
"
−1
−1
GL =
0.1−0.3z
1−1.29z −1 +0.84z −2
0.3−0.95z −1
1−0.65z −1 +0.84z −2
0.37−0.11z
1−0.28z −1 +0.72z −2
0.52−0.14z −1
1−0.34z −1 +0.83z −2
#
.
(4.44)
Rieˇsenie: Ked’ˇze model z´ıskan´
y MRI identifik´aciou je pouˇz´ıvan´
y d’alej ako model pre MPC,
bude n´as zauj´ımat’, ako dobre dok´aˇze predikovat’ v´
ystupy na P krokov. Ako jedno z krit´erii, ktor´e pouˇzijeme na hodnotenie kvality odhadu, bude preto (podobne ako v pr´ıpade
dvojkrokov´eho algoritmu) pouˇzit´a hodnota u
´ˇcelovej funkcie 4.35.
Vykresl’ovan´
y je op¨at’ pomer hodnoty u
´ˇcelovej funkcie pre parametre z´ıskan´e i-krokovou
identifik´aciou (i ∈ {1, 2, . . . , P }) ku hodnote u
´ˇcelovej funkcie pre parametre z´ıskane jednokrokovou identifik´aciou. Tento pomer budeme oznaˇcovat’ RP . Priebeh na obr. 4.5
naznaˇcuje, ˇze podobne ako v pr´ıpade dvojkrokov´eho algoritmu kvalita odhadu rastie
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
32
so zvyˇsuj´
ucim sa predikˇcn´
ym horizontom aˇz do urˇcit´eho horizontu, kedy sa ust´ali na urˇcitej
hodnote a d’alej sa nemen´ı.
1
0.9
RP
0.8
0.7
0.6
0.5
0.4
10
20
30
P
40
50
60
Obr. 4.5: Priebeh pomeru Rp v z´avislosti na P .
ˇ
Dalej
porovn´ame 50-krokov´e predikcie v´
ystupu s re´alnymi v´
ystupmi syst´emu pre jednokrokov´
u a P -krokov´
u identifik´aciu. Aj toto porovnanie (obr. 4.6) potvrdzuje, ˇze model
z´ıskan´
y klasick´
ymi met´odami nie je ani zd’aleka tak´
ym dobr´
ym viackrokov´
ym prediktorom ako model z´ıskan´
y minimaliz´aciou (4.35).
50 − kroková predikcia
50 − kroková predikcia
reálne
15
reálne
20
k−krok
fit = 83.99
1−krok
k−krok
1−krok
fit = 75.81
15
10
10
5
2
y (−)
y1(−)
5
0
0
−5
fit = 33.96
−5
fit = 26.62
−10
−10
50
100
150
200
t(h)
250
300
350
50
100
150
200
t(h)
250
300
350
Obr. 4.6: Priebehy P -krokov´
ych predikci´ı v´
ystupov.
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
4.4
33
MRI identifik´
acia v kombin´
acii s PLS
Na predchadzaj´
ucich riadkoch boli pop´ısan´e met´ody MRI identifik´acie, ktor´e boli n´asledne
su
´spechom otestovn´e na jednoduch´
ych pr´ıkladoch. Data pouˇzit´e na identifik´aciu boli vˇsak
zat’aˇzen´e nekorelovan´
ym bielym ˇsumom, ˇco v pr´ıpade identifik´acie re´alnych riaden´
ych
syst´emov neb´
yva ˇcast´
ym javom, pretoˇze m´alokedy je k dispoz´ıcii dostatok tak kvalitn´
ych
d´at. Je takmer pravidlom, ˇze data ktor´e s´
u pri identifik´acii dostupn´e, s´
u zl´e podmienen´e
vd’aka pr´ıtomnosti kolinearity a ˇsumu, ˇco v´
yrazne zhorˇsuje spr´avanie beˇzne pouˇz´ıvan´
ych
met´od zaloˇzen´
ych na klasick´
ych najmenˇs´ıch ˇstvorcoch (OLS - Ordinary Least Squares),
ba niekedy mˆoˇze viest’ aj k nemoˇznosti pouˇzitia t´
ychto met´od. Tento probl´em mˆoˇze vyrieˇsit’ pouˇzitie ˇciastoˇcnej line´arnej regresie (PLS - Partial Least Squares), ktor´a je uˇz
niekol’ko rokov s obl’ubou pouˇz´ıvana v ekonometri, chemometrii a biometrii [4, 12, 27, 45].
Pouˇzitie algoritmu PLS v identifik´acii modelu riaden´eho syst´emu by s´ıce rieˇsilo probl´em s nekvalitn´
ymi identifikaˇcn´
ymi datami, ale tento model by nebol ekvivalentn´
y minimaliz´acii viackrokovej predikˇcnej chyby a nebol by optim´alny pre n´asledn´e MPC riadenie.
Met´odu, ktor´a kombinuje oba tieto pr´ıstupy - ˇciastoˇcn´
u line´arnu regresiu a MRI identifik´aciu - uviedli autori v ˇcl´anku [25]. V tejto ˇcasti najprv uvedieme struˇcn´
y popis met´ody
PLS a n´asledne ju pouˇzijeme v spojen´ı s MRI identifik´aciou.
4.4.1
ˇ
Ciastoˇ
cn´
a line´
arna regresia
Budeme uvaˇzovat’ line´arny regresn´
y model:
(4.45)
Y = ZΘ,
Y (N × no ), Z(N × m), Θ(no × m), kde m = na + nb je poˇcet parametrov v regresnom
vektore, N ud´ava poˇcet vzorkov, no je poˇcet z´avisl´
ych premenn´
ych (poˇcet v´
ystupov).
Met´oda najmenˇs´ıch ˇstvorcov vedie na nasleduj´
uce rieˇsenie
ˆ = Z TZ
Θ
−1
Z T Y.
(4.46)
V ide´alnom pr´ıpade s´
u riadky v matici Z line´arne nez´avisl´e. V pr´ıpade, ˇze tomu tak nie
je, matica (Z T Z)−1 je singul´arna a t´
ym p´adom jej inverzia neexistuje a rieˇsen´ım tohto
probl´emu moˇze byt’ pr´ave spom´ınan´a met´oda ˇciastoˇcnej lien´arnej regresie.
Met´oda PLS vyuˇz´ıva met´odu Principial Component Analysis (PCA) v spojen´ı s viacn´asobnou line´arnou regresiou. Met´oda PCA je podrobnejˇsie pop´ısana napr´ıklad v [21,44].
Tu uvedieme iba jej struˇcn´
y popis.
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
34
PCA
Na z´aklade PCA met´ody mˆoˇzeme maticu Z(N × m) rozloˇzit’ na nasleduj´
ucu sumu tak,
aby bolo zachovan´e ˇco najv¨aˇcˇsie mnoˇzstvo inform´acie obsiahnut´e v Z:
T
T
T
Z = t 1 pT
1 + t2 p2 + . . . ta pp + Ep = Tp Pp + Ep ,
(4.47)
kde th , h ∈ {1, 2, . . . , p} s´
u st´lpcov´e vektory d´lˇzky N , ktor´e s´
u navz´ajom kolm´e a naz´
yvaj´
u
sa sk´ore (score vectors), a ph , h ∈ {1, 2, . . . , p} s´
u st´lpcov´e vektory d´lˇzky m a naz´
yvaj´
u sa
z´at’aˇzov´e vektory (loadings vectors). V matici Ep (N × m) je uloˇzen´
y ˇsum a redundantn´a
inform´acia.
V jednoduchosti mˆoˇzeme povedat’, ˇze met´oda PCA sa snaˇz´ı presk´
umat’ ˇstrukt´
uru d´at
s dimenziou m a redukovat’ ich dimenziu na niˇzˇsiu dimenziu p t´
ym, ˇze je otoˇcen´a mnoˇzina
d´at a zobrazen´a do priestoru s niˇzˇsou dimenziou pozd´lˇz smerov, v ktor´
ych je obsiahnut´e
najv¨aˇcˇsie mnoˇzstvo inform´acie (smery s najv¨aˇcˇs´ım rozptylom). Z´at’aˇzov´e vektory ph reprezentuj´
u kos´ıny uhlov, ktor´e odpovedaj´
u smerom s najv¨aˇcˇs´ım mnoˇzstvom inform´acie,
sk´ore th zase predstavuj´
u projekcie jednotliv´
ych d´at pozd´lˇz vybran´
ych principi´alnych smerov. Pri hl’adan´ı vektorov th a ph sa pouˇz´ıva iterat´ıvny algoritmus parci´alnych najmenˇs´ıch
ˇstvorcov, tzv. Non-linear Iterative Partial Least Squares (d’alej len NIPALS) [14].
Jednou z moˇznost´ı, ako za pomoci PCA a NILAPS algoritmu odhadovat’ nezn´ame
parametre regresn´eho modelu (4.45), je pouˇzitie Principial Component Regression (PCR)
[21], kde je pre predikciu Y zo Z pouˇzit´e:
Y = Tp C,
(4.48)
kde C(p × no ) mˆoˇzeme z´ıskat’ pomocou OLS:
Cˆ = (T T T )−1 T T Y.
(4.49)
Ked’ˇze st´lpce matice T boli vyberan´e tak, aby boli na seba kolm´e, matica T T T je diagon´alna matica, ktorej inverzia urˇcite existuje.
Dodajme, ˇze pre Y mˆoˇzeme pouˇzit’ vzt’ah (4.49) namiesto Y = Tp PpT C, lebo st´lpce
matice T s´
u line´arnou kombin´aciou st´lpcov v Z a naˇs´ım ciel’om bolo n´ajst’ line´arnu komˆ
bin´aciu st´lpcov regresn´eho vektoru Z pre predikciu Y . Pre odhadovan´
y parameter Θ
z rovnice (4.45) potom plat´ı
ˆ = Pp C.
ˆ
Θ
(4.50)
T´ato met´oda m´a vˇsak t´
u nev´
yhodu, ˇze neberie do u
´vahy rel´aciu medzi Z a Y , ˇco mˆoˇze
mat’ za n´asledok fakt, ˇze komponenty bud´
u vybran´e tak, ˇze z d´at bude odstr´anen´a pr´ave
uˇzitoˇcn´a inform´acia namiesto ˇsumu.
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
35
Met´oda PLS, ktor´a bude pop´ısan´a na najbliˇzˇs´ıch riadkoch, tento probl´em rieˇsi, vyber´a
komponenty tak, aby rovnako dobre popisovali Z aj Y .
PLS Podobne ako v pr´ıpade regresnej matice Z, aj maticu v´
ystupov Y mˆoˇzeme nap´ısat’
ako nasleduj´
ucu sumu:
Z = r1 q1T + r2 q2T + . . . ra qpT + Fp = Rp QT
p + Fp ,
(4.51)
kde sk´ore th , h ∈ {1, 2, . . . , p} s´
u st´lpcov´e ortogon´alne vektory d´lˇzky N a z´at’aˇzov´e vektory
ph , h ∈ {1, 2, . . . , p} s´
u st´lpcov´e vektory d´lˇzky no . V matici Fp (N × no ) je uloˇzen´
y ˇsum
a redundantn´a inform´acia. Na zabezpeˇcenie korel´acie medzi Z a Y je pouˇzit´a v¨azba medzi score vectors rh a th . My sa obmedz´ıme na najjednoduchˇsiu moˇzn´
u line´arnu v¨azbu:
rh = b h t h .
(4.52)
Za pomoci rovn´ıc (4.51), (4.47), (4.52) mˆoˇzeme pre predikciu Yˆ nap´ısat’
ˆ T,
Yˆ = T BQ
(4.53)
ˆ
kde B(p×p)
je diagon´alna matica, ktor´a m´a na diagon´alach koeficienty bh , h ∈ {1, 2, . . . , p}.
Na n´asledn´
y PLS odhad regresnej matice Θ existuje viacero algoritmov. Prv´
y z nich
je v literat´
ure ˇcasto oznaˇcovan´
y aj ako PLS1 algoritmus a ide o pomerne jednoduch´
y
neiterat´ıvny algoritmus, ktor´
y vˇsak v pr´ıpade no > 1 odhaduje parametre regresn´eho
vektora Θ zvl´aˇst’ pre kaˇzd´
y v´
ystup, ˇco je nie vˇzdy ˇziad´
uce.
Pokial’ chceme odhadovat’ naraz regresn´e parametre pre vˇsetky v´
ystupy, naskyt´a sa
pouˇzitie zloˇzitejˇsieho iterativn´eho ale komplexnejˇsieho algoritmu, ktor´
y je zaloˇzen´
y na
NIPALS algoritme (jeho podrobnejˇs´ı popis je uveden´
y napr´ıklad v [46]).
V pr´ıpade pouˇzitia tohto algoritmu mˆoˇze byt’ pre predikciu Yˆ z X pouˇzit´
y nasleduj´
uci
vzt’ah
ˆ T = Z Θ.
ˆ
(4.54)
Yˆ = ZW (P T W )−1 BQ
Matica W (m × p) je ortonorm´alna matica, st´lpce matice wj su jednotkov´e vektory, ktor´e
s´
u na seba kolm´e. Tieto vektory s´
u vyberan´e tak, aby vektor wj popisoval smer v priestore
Z, ktor´
y vykazuje najv¨aˇcˇsiu kovariancu medzi Z a pr´ısluˇsn´
ym st´lpcom matice v´
ystupov
Y . To znamen´a, ˇze vel’k´a zmena v hodnote Z spˆosob´ı vel’k´
u zmenu v pr´ısluˇsnej hodnote Y .
Tento algoritmus je s´ıce omnoho komplexnejˇs´ı ako algoritmus oznaˇcovan´
y ako PLS1, ale je
iterat´ıvny a pomerne v´
ypoˇctovo n´aroˇcny. Tieto probl´emy rieˇsi algoritmus oznaˇcovan´
y ako
SIMPLS uveden´
y v [9]. Algoritmus SIMPLS nehl’ad´a z´at’aˇzov´e vektory qj iterat´ıvne ako
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
36
predchadzaj´
uci algoritmus, ale vyuˇz´ıva to, ˇze q1 predstavuje dominantn´
y vlastn´
y vektor
matice (X T Y )T (X T Y ), ktor´
y moˇzno n´ajst’ pomerne jednoducho prostredn´ıctvom SVD
rozkladu, ˇco znaˇcne zjednoduˇs´ı a zr´
ychli cel´
y algoritmus. Pre predikciu Y z X moˇzno
podobne ako v predoˇslom pr´ıpade pouˇzit’ vzt’ah (4.54).
4.4.2
Popis algoritmu
Postup pop´ısan´
y v predchadzaj´
ucom texte sa d´a v´
yhodne pouˇzit’ aj pri minimaliz´acii
u
´ˇcelovej funkcie (4.36). Ako uˇz bolo viackr´at spomenut´e, minimaliz´acia je u
´lohou neline´arnej viacrozmernej optimaliz´acie, ktorej met´ody s´
u podrobne pop´ısan´e v mnoˇzstve
kvalitnej literat´
ury [3, 32]. Tieto met´ody s´
u zaloˇzen´e na iterat´ıvnom hl’adan´ı optim´alneho
parametra Θ, pre ktor´
y v kaˇzdom kroku algoritmu plat´ı
Θk+1 = Θk + pk αk ,
(4.55)
kde pk predstavuje smer hl’adania vo viacerozmernom priestore s dimenziou rovnou poˇctu
odhadovan´
ych parametrov v regresnej matici Θ a αk predstavuje d´lˇzku kroku. Na urˇcenie
parametrov pk a αk sa pouˇz´ıvaj´
u rˆozne pr´ıstupy. V naˇsom pr´ıpade vyuˇzijeme pri hl’adan´ı
optimaln´eho smeru pk nasleduj´
uci postup, ktor´
y vyuˇz´ıva PLS a je pop´ısan´
y v [25].
Smer pk je hl’adan´
y tak, aby funkcia JLRP I dosahovala svojho minima v bode Θk + pk ,
ˇco dosiahneme prostredn´ıctvom nasleduj´
ucej podmienky:
∂JLRP I
= 0.
∂(Θk + pk )
(4.56)
∂JLRP I
T
T
= −2Za|k
Ya + 2Za|k
Za|k (Θk + pk ).
∂(Θk + pk )
(4.57)
Derivovan´ım dost´avame:
Deriv´aciu (5.8) poloˇz´ıme rovn´
u nule a pre pk dostaneme
T
T
pk = (Za|k
Za|k )−1 Za|k
Ya − Θ k ,
(4.58)
ˇco odpoved´a v´
yrazu pre (4.49) posunut´emu o θk . Pr´ave v tomto kroku mˆoˇze nastat’
T
probl´em s invertovan´ım matice Za|k
Za|k , ktor´
y vˇsak mˆoˇzeme jednoducho vyrieˇsit’ pr´ave
za pomoci PLS. Namiesto klasick´eho LS odhadu pre pk pouˇzijeme ˇciastoˇcn´
u line´arnu
regresiu. Potom pre pk dost´avame:
ˆ T − Θk ,
pk = W (P T W )−1 BQ
(4.59)
´
KAPITOLA 4. IDENTIFIKACIA
PRE PREDIKT´IVNE RIADENIE
37
ˆ z´ıskame aplikovan´ım PLS na Za|k a Ya . Ost´ava n´ajst’ uˇz len optim´alnu
kde W , P , Q, B,
d´lˇzku kroku αk tak, aby platilo
αk = min JLRP I (θk + αpk ).
α
(4.60)
Ako najjednoduchˇs´ı spˆosob minimaliz´acie (4.60) sa jav´ı aproxim´acia JLRP I (θk + αpk )
kvadratickou funkciou
JLRP I (θk + αpk ) ≈ a + bα + cα2 .
(4.61)
Aproxim´acia polyn´omom druh´eho stupˇ
na je volen´a kvˆoli faktu, ˇze hodnoty a, b, c a t´
ym aj
potrebn´
u aproxim´aciu sme schopn´ı urˇcit’ len zo znalosti JLRP I (θk + αpk ) v troch bodoch.
Pre minimum kvadratickej funkcie (4.62) existuje analytick´
y vzt’ah:
αk = −
b
.
2c
(4.62)
Dosaden´ım αk a pk (4.55) do rovnice dost´avame nov´
u hodnotu Θk+1 a algoritmus pokraˇcuje d’alˇsou iter´aciou za predpokladu:
JLRP I (Θk ) − JLRP I (Θk+1 ) > tol,
(4.63)
kde tol > 0 je volitel’n´
y parameter. Ak podmienka (4.63) nie je splnen´a, algoritmus skonˇcil
a hodnota Θk+1 je hl’adan´
ym optim´alnym parametrom.
Kapitola 5
Predikt´ıvne riadenie
Predikt´ıvne riadenie je uˇz niekol’ko rokov s obl’ubou pouˇz´ıvan´e pri riaden´ı rˆoznych priemyseln´
ych procesov. Ide o jednu z najˇcastejˇsie pouˇz´ıvan´
ych pokroˇcil´
ych met´od riadenia,
ktor´a nach´adza uplatnenie pri riaden´ı v najrozliˇcnejˇs´ıch oblastiach priemyslu.
Na obr. 5.1 je zobrazen´a ˇstrukt´
ura procesu riadenia. Najniˇzˇsou vrstvou je vrstva
inˇstrument´acie, ktor´a zah´rn
ˇa senzory a aktu´atory a tvor´ı rozhranie medzi samotn´
ym
technologick´
ym procesom a regul´atorom. Jadro riadiaceho syst´emu je tvoren´e vrstvou
z´akladn´eho riadenia. T´ato vrstva za vyuˇzitia z´akladn´
ych monitorovac´ıch a vizualizaˇcn´
ych
n´astrojov zaist’uje z´akladn´
u funkˇcnost’ a bezpeˇcnost’ a z´aroveˇ
n poskytuje podporu vrstve
pokroˇcil´eho riadenia. Hlavn´
ym ciel’om vrstvy pokroˇcil´eho riadenia je zabezpeˇcit’ optim´alne
spr´avanie riaden´eho procesu pri dodrˇzan´ı zadan´
ych podmienok a obmedzen´ı.
Ked’ˇze predikt´ıvny regul´ator dok´aˇze zabezpeˇcit’ optim´alne spr´avanie riadenej s´
ustavy,
ˇco je zabezpeˇcen´e vhodnou vol’bou u
´ˇcelovej funkcie, ale aj dodrˇzanie obmedzuj´
ucich podmienok, je vhodn´
ym kandid´atom na nasadenie vo vrstve pokroˇcil´eho riadenia. Najvyˇsˇsia
pl´anovacia vrstva definuje poˇziadavky pre riadenie ako napr´ıklad obmedzenia, optim´alne
podmienky a pod.
38
KAPITOLA 5. PREDIKT´IVNE RIADENIE
39
ˇ
Obr. 5.1: Strukt´
ura procesu riadenia.
Predikt´ıvne riadenie je ˇcasto aplikovan´e na mnohorozmern´e syst´emy ˇci syst´emy s dopravn´
ym oneskoren´ım. Narozdiel od in´
ych met´od poˇc´ıta MPC regul´ator akˇcn´
y z´asah
nielen pre nasleduj´
ucu vzorkovaciu peri´odu, ale na z´aklade predikcie bud´
uceho spr´avania
syst´emu n´ajde cel´
u optim´alnu postupnost’ na dobu horizontu predickcie P . Optim´alna
postupnost’ akˇcn´
ych z´asahov je urˇcen´a na z´aklade rieˇsenia optimalizaˇcnej u
´lohy [3, 49]
na koneˇcnom horizonte.
Obr. 5.2: Blokov´a sch´ema predikt´ıvneho regul´
atoru.
Z´akladn´
y princ´ıp predikt´ıvneho regul´atoru je na obr. 5.2. Na z´aklade modelu procesu
sa vypoˇc´ıtaj´
u predpokladan´e bud´
uce hodnoty jednotliv´
ych veliˇc´ın, spoˇc´ıtan´e predikovan´e
trajekt´orie stavov (v´
ystupov) sa pouˇzij´
u na rieˇsenie optimalizaˇcnej u
´lohy a stanov´ı sa
KAPITOLA 5. PREDIKT´IVNE RIADENIE
40
optim´alna postupnost’ bud´
ucich akˇcn´
ych z´asahov na dobu horizontu predikcie. T´ato postupnost’ je vypoˇc´ıtan´a tak, aby boli splnen´e vˇsetky zadan´e obmedzenia na riadiace ˇci
riaden´e veliˇciny procesu.
Pokial’ by sme aplikovali cel´
u vypoˇc´ıtan´
u postupnost’ akˇcn´
ych z´asahov na celom koneˇcnom horizonte predikcie, n´asledne spoˇc´ıtali nov´
u postupnost’ a aplikovali, riadili by
sme syst´em po prvom kroku v otvorenej sluˇcke. Tak´
yto postup vˇsak nie je pr´ıliˇs v´
yhodn´
y,
ked’ˇze ide o riadenie v otvorenej sluˇcke, nie je moˇzn´e reagovat’ na poruchy, ktor´e na re´alne
ˇ sou v´
syst´emy ˇcasto pˆosobia. Dalˇ
yraznou nev´
yhodou tohto pr´ıstupu je znaˇcn´a nerobustnost’.
Rieˇsen´ım tohto probl´emu je pouˇzitie tzv. k´lzav´eho horizontu predikcie (receding horizon) [24]. V kaˇzdom kroku sa spoˇc´ıta optim´alna riadiaca postupnost’ na dobu horizontu
predikcie P , ale aplikovan´
y je len prv´
y akˇcn´
y z´asah. Po zmeran´ı nov´eho stavu sa cel´
y
postup opakuje. T´
ymto je do riadenia zaveden´a sp¨atn´a v¨azba. Princ´ıp posuvn´eho horizontu predikcie vysvetl’uje obr. 5.3. Dodajme eˇste, ˇze aplik´acia posuvn´eho horizontu je
vlastne aproxim´aciou optimalizaˇcnej u
´lohy na koneˇcnom horizonte u
´lohou s nekoneˇcn´
ym
horizontom.
Jednou z nev´
yhod klasick´eho predikt´ıvneho regul´atora je fakt, ˇze cel´a optimalizaˇcn´a
u
´loha je rieˇsen´a on-line, pre kaˇzd´
y akˇcn´
y z´asah. Vel’k´a v´
ypoˇctov´a n´aroˇcnost’ ˇcasto zabraˇ
nuje pouˇzitiu predikt´ıvnych algoritmov riadenia v programovatel’n´
ych automatoch
alebo pri riaden´ı r´
ychlejˇs´ıch syst´emov s malou peri´odou vzorkovania.
Rieˇsen´ım mˆoˇze byt’ pouˇzitie explicitn´eho predikt´ıvneho regul´atoru, kde sa dan´
y probl´em vyrieˇsi naraz pre vˇsetky pr´ıpustn´e poˇciatoˇcn´e podmienky pomocou parametrick´eho
programovania [49]. V´
ysledn´
y z´akon riadenia je vo forme tabul’ky. Touto met´odou sa vˇsak
nebudeme bliˇzˇsie zaoberat’, jej podrobn´
y popis je uveden´
y napr´ıklad v [2, 34].
Obr. 5.3: Posuvn´
y horizont predikcie.
KAPITOLA 5. PREDIKT´IVNE RIADENIE
5.1
41
N´
avrh predikt´ıvneho regul´
atora
Pri rieˇsen´ı u
´lohy predikt´ıvnej regul´acie je dˆoleˇzit´a dobr´a znalost’ matematick´eho modelu syst´emu. Na z´aklade modelu sme schopn´ı predikovat’ bud´
uce stavy syst´emu a n´ajst’
optim´alnu postupnost’ akˇcn´
ych z´asahov vyrieˇsen´ım optimalizaˇcn´eho probl´emu tak, aby
boli splnen´e zadan´e obmedzenia. Budeme sa zaoberat’ n´avrhom predikt´ıvneho regul´atoru
pre line´arny dynamick´
y ˇcasovo invariantn´
y syst´em pop´ısan´
y nasleduj´
ucimi diferenˇcn´
ymi
stavov´
ymi rovnicami:
x(k + 1) = Ax(k) + Bu(k),
y(k) = Cx(k) + Du(k),
(5.1)
´
kde x(k) ∈ ℜn , u(k) ∈ ℜm , y(k) ∈ ℜp . Ulohou
predikt´ıvnej regul´acie je n´ajst’ optim´alnu
riadiacu postupnost’ minimalizuj´
ucu vhodne zvolen´e krit´erium na koneˇcnom horizonte
predikcie P . Vol’ba samotn´eho krit´eria z´avis´ı od konkr´etnej u
´lohy. Vol´ıme ho tak, aby
sme mohli penalizovat’ neˇziaduce spr´avanie ako napr. vel’k´
y akˇcn´
y z´asah, odch´
ylku od poˇzadovanej hodnoty alebo ust´alen´eho stavu.
ˇ
Castou
vol’bou je kvadratick´e krit´erium, ktor´e m´a n´azorn´
u fyzik´alnu interpret´aciu mnoˇzstvo energie a existuje pomerne vel’k´e mnoˇzstvo algoritmov pre rieˇsenie u
´loh kvadratickej optimaliz´acie. Najˇcastejˇsie pouˇz´ıvan´
ymi s´
u zobecnen´a Cholesky faktoriz´acia a LDU
(Biermanova) faktoriz´acia. Podrobnejˇs´ı popis vybran´
ych met´od je uveden´
y napr´ıklad
v [32, 49]. Pre n´aˇs pr´ıpad budeme volit’ kvadratick´e krit´erium tak, aby bola penalizovan´a odch´
ylka od ust´alen´eho stavu a vel’k´
y akˇcn´
y z´asah:
JP (x(k), u(k)) =
P
−1
X
x(k + i|k)T qx(k + i|k) + u(k + i|k)T ru(k + i|k).
(5.2)
i=0
Matice q = q T ≥ 0 a r = rT > 0 s´
u v´ahov´e matice, ktor´e s´
u spolu s horizontom predikcie
´
P hlavn´
ymi ladiacimi parametrami pri n´avrhu MPC. Ulohou
predikt´ıvnej regul´acie je
n´ajst’ minimum kvadratickej formy 5.2 vzhl’adom na premenn´
uu
u∗ = arg min JP (x(k), u(k)).
u
(5.3)
Obecne plat´ı, ˇze so zvyˇsovan´ım horizontu predikcie sa zvyˇsuje kvalita regul´acie, ale
z´aroveˇ
n aj v´
ypoˇctov´a zloˇzitost’. Je dˆoleˇzit´e zvolit’ horizont predikcie tak, aby bol poˇc´ıtaˇc
schopn´
y spoˇc´ıtat’ aktu´alny akˇcn´
y z´asah za dobu menˇsiu ako je vzorkovacia peri´oda. V´
yraz
x(k + i|k) predstavuje predikovan´
u hodnotu stavu v ˇcase k + i na z´aklade d´at v ˇcase k.
ˇ
Dalej
budeme pre zjednoduˇsenie pouˇz´ıvat’ znaˇcenie x(k + i|k) = xi .
KAPITOLA 5. PREDIKT´IVNE RIADENIE
42
Predikovan´e hodnoty stavov na horizonte predikcie P vyjadr´ıme ako
Matice Ap a Bp
Oznaˇcme
h
xT
1
xT
2
. . . xT
P −1
iT





Ap = 


A



,


A2
..
.
AP −1
iT
h
T
T
= Ap x o + B p u T
u1 . . . uP −1 .
0


B

 AB
B

Bp =  .
...
 ..

AP −2 B B
0



.


(5.4)
(5.5)
iT
h
T
T
U = uT
u1 . . . uP −1 .
0
iT
h
T
X = xT
x2 . . . xP −1 ,
1
Krit´erium 5.2 mˆoˇzeme potom p´ısat’ v maticovej forme:
JP = (Ap x0 + BpU )T Q (Ap x0 + BpU ) + U T RU.
h
i
h
(5.6)
i
V´ahov´e matice Q = IP −1 ⊗ q a R = IP −1 ⊗ r , kde ⊗ znaˇc´ı Kroneckerov produkt
mat´ıc a IP predstavuje jednotkov´
u maticu pr´ısluˇsnej dimenzie. Po rozn´asoben´ı a u
´prav´ach
dost´avame
T
T T
JP = U T (BpT QBp + R)U + xT
0 Ap QBp U + U Bp QAp x0 .
(5.7)
S´
uˇcast’ou krit´eria 5.7 by mal byt’ aj absol´
utny ˇclen, ktor´
y vˇsak pri minimaliz´acii krit´eria
nemus´ıme uvaˇzovat’.
V pr´ıpade, ˇze neuvaˇzujeme ˇziadne obmedzenia, mˆoˇzeme u
´lohu minimaliz´acie 5.7 rieˇsit’
’ yraz podl’a U a optim´alnu postupnost’
analyticky. Jednou z hmoˇznost´ı je zderivovat
iT v´
hl’adat’ podl’a
akˇcn´
ych z´asahov U = u∗T
u∗T
. . . u∗T
0
1
P −1
∂JP
= 0,
∂U
pre U = U ∗ .
(5.8)
Druh´
y moˇzn´
y spˆosob, ktor´
ym sa budeme zaoberat’ podrobnejˇsie, spoˇc´ıva v u
´prave na
u
´pln´
y ˇstvorec.
Zavedieme si oznaˇcenie
H = (BpT QBp + R),
T
f = xT
0 Ap QBp .
(5.9)
V´
yraz 5.7 uprav´ıme do nasleduj´
uceho tvaru:
JP = U T − U ∗ H (U − U ∗ ) − U ∗T HU ∗ .
(5.10)
KAPITOLA 5. PREDIKT´IVNE RIADENIE
43
Aby bolo 5.7 rovn´e 5.10, mus´ı platit’
−U ∗T HU = f U,
−U T HU ∗ = U T f T .
(5.11)
ˇ
Dalej
mˆoˇzeme p´ısat’
−U ∗T H = f,
−HU ∗ = f T .
(5.12)
Matica H je pozit´ıvne semidefinitn´a symetrick´a matica, H = H T , H ≥ 0. Po transponovan´ı prvej rovnice a sˇc´ıtan´ı oboch rovn´ıc 5.12 dost´avame
−HU ∗ = f T .
(5.13)
Potom pre optim´alnu postupnost’ akˇcn´
ych z´asahov U ∗ minimalizuj´
ucu kvadratick´
u formu
5.7 plat´ı
U ∗ = −(BpT QBp + R)−1 BpT QAp x0 .
5.1.1
(5.14)
Predikt´ıvny regul´
ator s obmedzen´ım
Takmer kaˇzd´a riaden´a s´
ustava obsahuje nejak´e obmedzenie, ˇci uˇz ide o fyzick´e obmedzenia senzorov a akˇcn´
ych ˇclenov, alebo o rˆozne technologick´e a ekonomick´e obmedzenia
riaden´
ych veliˇc´ın. Jednou z hlavn´
ych v´
yhod MPC regul´atora je moˇznost’ zahrnutia t´
ychto
obmedzen´ı.
Optimalizaˇcn´
uu
´lohu uˇz nemˆoˇzeme rieˇsit’ analyticky ako v pr´ıpade bez obmedzen´ı, ale
´
vedie na kvadratick´e programovanie. Ulohou
kvadratick´eho programovania je n´ajst’ tak´e
z, ktor´e minimalizuje funkciu:
J(z) =
1 T
z Hz + f T z
2
(5.15)
tak, aby boli splnen´e zadan´e obmedzenia, ktor´e s´
u vyjadren´e line´arnymi maticov´
ymi
nerovnost’ami:
Ainq z ≤ binq ,
Aeq z = beq .
(5.16)
Na minimaliz´aciu krit´eria 5.7 sa mˆoˇzeme pozerat’ ako na kvadratick´
y program, ktor´eho
optimalizaˇcnou premennou je vektor bud´
ucich akˇcn´
ych z´asahov U :
min U T HU + f T U,
U
Ainq U ≤ binq ,
Aeq U = beq .
(5.17)
Pomocou kvadratick´eho programu 5.17 n´ajdeme optim´alnu postupnost’ akˇcn´
ych z´asahov
na horizonte predikcie P tak, aby boli splnen´e obmedzuj´
uce podmienky. Na rieˇsenie u
´loh
KAPITOLA 5. PREDIKT´IVNE RIADENIE
44
kvadratick´eho programovania existuje mnoˇzstvo solverov, napr´ıklad v programovom prostred´ı Scilab1 je implementovan´a funkcia qpsolve, ktor´a m´a nasleduj´
uci z´apis
[U ] = qpsolve(H, f, Ainq , binq , me ),
(5.18)
kde me ud´ava index prvku vektora b, od ktor´eho bud´
u obmedzuj´
uce podmienky urˇcen´e
v tvare rovnosti. Ako uˇz bolo spomenut´e, mnohokr´at je potrebn´e obmedzovat’ nielen akˇcn´
y
z´asah, ktor´
y je optimalizaˇcnou premennou samotn´eho kvadratick´eho programu, ale i riaden´e veliˇciny (stavy) ˇci samotn´e v´
ystupy. Vyjadr´ıme si predikovan´e hodnoty stavov z 5.4
a obmedzuj´
uce podmienky pre stavov´e premenn´e bud´
u v tvare nasleduj´
ucich maticov´
ych
nerovnost´ı (rovnost´ı)
Ainq Bp U ≤ binq − Ainq Ap x0 ,
Aeq Bp U = beq − Aeq Ap x0 .
(5.19)
ˇ s´ım z ˇcast´
Dalˇ
ych obmedzen´ı je obmedzenie na zmenu riadiacej veliˇciny veliˇciny ∆u(k) =
u(k) − u(k − 1)
Ainq ∆U ≤ binq ,
Aeq ∆U = beq .
(5.20)
Pre zmenu akˇcnej veliˇciny ∆u(k) plat´ı
∆U = Di U −
h
uT
0
0 ···
0
iT
.
(5.21)
Tvar matice Di z´avis´ı na poˇcte vstupov syst´emu, pre syst´em s jedn´
ym vstupom plat´ı


1


−1 1



D=
(5.22)
.
.
.
.. ..




−1 1
Pomocou vzt’ahu 5.22 sme schopn´ı prep´ısat’ obmedzenia 5.20 do poˇzadovan´eho tvaru, kde
vystupuje len optimalizaˇcn´a premenn´a U .
Obmedzenia, ktor´e sme doposial’ spom´ınali, s´
u takzvan´e tvrd´e obmedzenia, nemˆoˇzu
byt’ za ˇziadnych okolnost´ı prekroˇcen´e. Zav´adzanie tvrd´
ych omedzen´ı je vhodn´e pokial’ ide
o limity samotn´eho procesu (fyzick´e obmedzenia senzorov a aktu´atorov).
Niekedy je vˇsak vhodn´e tak´eto obmedzenia nejak´
ym spˆosobom zm¨akˇcit’, a to tak, ˇze
na urˇcit´
y ˇcas povol´ıme prekroˇcenie dan´eho limitu za cenu urˇcitej penaliz´acie.
Deje sa tak zaveden´ım novej premennej ǫ:
AU ≤ b + ǫ
1
http://www.scilab.org/
(5.23)
KAPITOLA 5. PREDIKT´IVNE RIADENIE
45
Do kvadratick´eho krit´eria prid´ame ˇclen ǫT Eǫ, ktor´
ym budeme penalizovat’ pr´ıpadn´e prekroˇcenie dan´eho obmedzenia. V´ahov´a matica E je symetrick´a pozit´ıvne definitn´a matica
pr´ısluˇsn´
ych rozmerov, volen´a spravidla o niekol’ko r´adov vyˇsˇsia ako v´ahy na stavy a akˇcn´e
z´asahy.
Tak´eto obmedzenia sa pouˇz´ıvaj´
u v¨aˇcˇsinou na obmedzenia v´
ystupov alebo stavov pokial’ chceme penalizovat’ urˇcit´e neˇziaduce spr´avanie alebo zn´ıˇzenie kvality procesu. M¨akk´e
obmedzenia su dˆoleˇzit´e aj z hl’adiska rieˇsitel’nosti optimalizaˇcnej u
´lohy, o ktorej sa zmienime neskˆor.
ˇ
Dalej
je moˇzn´e prid´avat’ rozliˇcn´e obmedzenia, zmenit’ krit´erium (5.7) prid´avan´ım v´ah.
ˇ
Casto
je to obmedzenie na zmenu v´
ystupu alebo stavu ∆x(k) = x(k) − x(k − 1), pr´ıpadne
sa v krit´eriu nev´aˇzi samotn´
y akˇcn´
y z´asah, ale jeho zmena a vtedy m´a regul´ator integraˇcn´
y
charakter.
Treba vˇsak brat’ do u
´vahy fakt, ˇze optimaliz´acia prebieha spravidla on-line a pridan´ım
kaˇzd´eho d’alˇsieho obmedzenia ˇci v´ahy sa znaˇcne zv´
yˇsi vypoˇctov´a zloˇzitost’. Preto nie je
vhodn´e prid´avat’ ˇziadne zbytoˇcn´e obmedzenia. Nem´a napr´ıklad zmysel v´aˇzit’ alebo obmedzovat’ poˇciatoˇcn´
y stav x0 .
5.1.1.1
Sledovanie konˇ
stantnej referencie
Doteraz sme sa zaoberali u
´lohou ako n´ajst’ optim´alnu postupnost’ akˇcn´
ych z´asahov, ktor´a
prevedie syst´em z l’ubovol’n´eho poˇciatoˇcn´eho stavu x0 do poˇciatku s´
uradnicovej s´
ustavy.
V praxi vˇsak ˇcasto poˇzadujeme, aby syst´em sledoval zadan´
u referenˇcn´
u trajekt´oriu. Aby
sme mohli u
´lohu n´avrhu predikt´ıvneho regul´atoru modifikovat’ pre sledovanie po ˇcastiach
konˇstantnej referencie yref (k), mus´ıme predikovat’ jej hodnoty na dobu horizontu predikcie
N . Prvou moˇznost’ou je, ˇze pozn´ame referenˇcn´
u trajekt´oriu na celom horizonte predikcie,
ˇco bude aj n´aˇs pr´ıpad (programov´e riadenie) [20]. Ked’ hodnota referencie nie je zn´ama,
mˆoˇzeme ju modelovat’ ako n´ahodn´
u prech´adzku:
yref (k + 1) = yref (k) + w(k),
(5.24)
kde w(k) je biely ˇsum a stredn´
u hodnotu yˆref odhadujeme ako konˇstantu yˆref (k + 1) =
yˆref (k). Op¨at’ budeme uvaˇzovat’ model syst´emu pop´ısan´
y stavov´
ymi rovnicami (5.1). Definujeme si kvadratick´e krit´erium tak, aby bola penalizovan´a odch´
ylka od referenˇcn´eho
sign´alu e(k) = y(k) − yref (k):
JN =
N
−1
X
i=0
˜ui ,
˜ei + uT
eT
i r
i q
(5.25)
KAPITOLA 5. PREDIKT´IVNE RIADENIE
46
T
T
matice
u v´ahov´e matice. Predikovan´e hodnoty v´
ystupov
h q˜ = q˜ ≥ 0 a ir˜T = r˜ > 0 s´
T
na horizonte predikcie N vyjadr´ıme ako:
Y = y0T y1T . . . yN
−1
h
iT
T
T
Y = A˜p xo + B˜p uT
,
u
.
.
.
u
0
1
N −1
kde matice A˜p a B˜p s´
u v tvare:


C


 CA



˜
Ap =  .
,
 ..



CAP −1




˜
Bp = 


(5.26)

D
CB
..
.
D
CAP −2 B
...
CB D



.


(5.27)
Analogicky ako v pr´ıpade riadenia do poˇciatku nap´ıˇseme krit´erium v maticovej forme:
T ˜p U − Yref Q
˜ A˜p x0 + B
˜p U − Yref ) + U T RU,
˜
JP = A˜p x0 + B
(5.28)
kde Yref je vektor predikovan´
ych hodnˆot referenˇcnej trajekt´orie na dobu horizontu predikcie.
Po rozn´asoben´ı a u
´prav´ach 5.28 mˆoˇzeme optim´alnu postupnost’ riadiacich veliˇc´ın,
ktor´a zaist´ı sledovanie referenˇcn´eho sign´alu Yref , hl’adat’ pomocou kvadratick´eho programu
1 T˜
U HU + f˜T U, Ainq U ≤ binq ,
U 2
˜ v tomto pr´ıpade plat´ı
Pre v´ahove matice f˜ a H
min
Aeq U = beq .
˜ =B
˜pT Q
˜B
˜p + R,
˜
H
T
˜
˜
˜B
˜p .
f = Ap x0 − Yref Q
5.1.1.2
(5.29)
(5.30)
Stabilita a rieˇ
sitel’nost’
Riadenie met´odou k´lzav´eho horizontu eˇste nezaruˇcuje stabilitu, dokonca sa mˆoˇze stat’,
ˇze pre niektor´e poˇciatoˇcn´e podmienky nebude mat’ optimalizaˇcn´a u
´loha vˆobec rieˇsenie.
V dev¨at’desiatych rokoch minul´eho storoˇcia boli publikovan´e prv´e ˇcl´anky [26, 37], ktor´e
sa zaoberali garanciou stability pri pouˇzit´ı k´lzav´eho horizontu. Stabilitu moˇzno zaruˇcit’:
• pridan´ım v´ahy na stav na konci horizontu predikcie xP (terminal weight),
• pridan´ım obmedzenia na stav na konci horizontu predikcie xP (terminal constraint
set).
KAPITOLA 5. PREDIKT´IVNE RIADENIE
47
Tieto met´ody s´
u podrobnejˇsie vr´atane dˆokazov pop´ısane napr´ıklad v [31]. Ako bolo spomenut´e, s rast´
ucim horizontom predikcie P rastie aj kvalita regul´acie. Definujme si kvadratick´e krit´erium na nekoneˇcnom horizonte:
JP (x(k), u(k)) =
N
−1
X
i=0
x(k + i|k)T qx(k + i|k) + u(k + i|k)T ru(k + i|k) + ψP∞ .
(5.31)
Krit´erium sme rozdelili na dve ˇcasti. Prv´a ˇcast’ predstavuje optimaliz´aciu na koneˇcnom
horizonte predikcie a budeme ju op¨at’ rieˇsit’ pomocou kvadratick´eho programovania tak,
aby boli splnen´e vˇsetky poˇzadovan´e obmedzenia. Druh´a ˇcast’ ψP∞ predstavuje optimaliz´aciu na nekoneˇcnom horizonte a plat´ı:
ψP∞ = xT
P P xP ,
ui = −Kxi pre i ≥ P,
(5.32)
kde P je ust´alen´
ym rieˇsenim diferenˇcnej Riccatiho rovnice (ARE)
P = AT P A − AT P B(Rr + B T P B)−1 B T P A + Qr ,
(5.33)
a K je pr´ısluˇsn´e ust´alen´e zosilnenie (Kalmanovo zosilnenie)
K = (Rr + B T P B)−1 B T P A.
(5.34)
Predt´
ym, neˇz struˇcne pop´ıˇseme niektor´e z moˇznost´ı dosiahnutia zaruˇcenej stability pre k´lzav´
y horizont predikcie, uvedieme defin´ıcie pozit´ıvne invariantnej pr´ıpustnej mnoˇziny Ω
a mnoˇziny rieˇsitel’nosti Xf .
Defin´ıcia 5.1 (mnoˇ
zina rieˇ
sitel’nosti): Xf je tak´a mnoˇzina vˇsetk´
ych stavov xi ∈ ℜn ,
pre ktor´e je dan´
y kvadratick´
y program rieˇsitel’n´
y, to znamen´a, ˇze existuj´e tak´a optim´alna
postupnost’ akˇcn´
ych z´asahov, ktor´a sp´lˇ
na vˇsetky zadan´e obmedzenia.
◮
Defin´ıcia 5.2 (pozit´ıvne invariantn´
a pr´ıpustn´
a mnoˇ
zina): Ω pre dan´e riadenie ui =
h(xi ) je oblast’ stavov´eho priestoru, pre ktor´
u plat´ı:
• xi ∈ Ω ⇒ xi+j ∈ xi pre ∀j > 0,
• xi ∈ Ω ⇒ xmin ≤ xi ≤ xmax pre ∀i ≥ 0.
◮
Defin´ıcia 5.2 hovor´ı, ˇze pokial’ sa stav dostane do mnoˇziny Ω a bude zachovan´
y dan´
y
z´akon riadenia ui = h(xi ), uˇz v tejto mnoˇzine zostane a z´aroveˇ
n bud´
u zachovan´e vˇsetky
obmedzuj´
uce podmienky. Jedn´
ym z moˇzn´
ych spˆosobov ako zabezpeˇcit’ stabilitu je vol’ba
KAPITOLA 5. PREDIKT´IVNE RIADENIE
48
vhodnej v´ahy koncov´eho stavu xP . Budeme uvaˇzovat’ krit´erium 5.31. Na koneˇcnom horizonte N m´ame zaruˇcen´
u optimalitu a splnenie obmedzuj´
ucich podmienok, pokial’ budeme
koncov´
y stav xP v´aˇzit’ v´ahovou maticou P , ktor´a je ust´alen´
ym rieˇsen´ım diferenˇcnej Riccatiho rovnice, m´ame zaruˇcen´
u stabilitu. Pokial’ vˇsak nijak´
ym spˆosobom neobmedz´ıme
stav xP a bude platit’:
xP ∈ ℜn ,
(5.35)
tento spˆosob vedie na suboptim´alne rieˇsenie, lebo pre i ≥ P nem´ame zaruˇcen´e splnenie obmedzuj´
ucich podmienok. Ak vˇsak pre vˇsetky poˇciatoˇcn´e stavy v kaˇzdom kroku
optimaliz´acie plat´ı
x0 ∈ X f ,
(5.36)
tak pre dostatoˇcne vel’k´
y horizont predikcie bude zaruˇcen´a nielen stabilita ale aj splnenie vˇsetk´
ych obmedzen´ı. N´ajst’ tak´
uto d´lˇzku je vˇsak pomerne zloˇzit´e a zaoberaj´
u sa
t´
ym napr´ıklad autori ˇcl´anku [2]. Jednoduchˇs´ı spˆosob pre n´ajdenie optim´alne rieˇsenia je
pridanie obmedzenia na koncov´
y stav.
Jednou z moˇznost´ı ako obmedzit’ koncov´
y stav je striktn´e obmedzenie xP = 0 (terminal equality constraints). T´ato moˇznost’ je s´ıce najjednoduchˇsia, znaˇcne t´
ym vˇsak z´
uˇzime
mnoˇzinu rieˇsitel’nosti, preto ju vo v¨aˇcˇsine pr´ıpadov nie je vhodn´e pouˇzit’. Ovel’a vhodnejˇsie je volit’ menej striktn´e obmedzenie na koncov´
y stav, a to xP ∈ ΩM , kde ΩM je
maxim´alna pozit´ıvne invariantn´a mnoˇzina pre u = h(xi ) = −Kxi . Z druh´eho bodu definicie pozit´ıvne invariantnej mnoˇziny 5.2 je jasn´e, ˇze pre takto volen´e obmedzenie na stav
na konci horizontu predikcie bud´
u splnen´e vˇsetky obmedzenia aj pre i ≥ N .
Pre syst´emy, ktor´e s´
u stabiln´e v otvorenej sluˇcke, sa naskyt´a d’alˇsia moˇznost’ ako
∞
zabezpeˇcit’ stabilitu. Op¨at’ budeme uvaˇzovat’ krit´erium 5.31, ψN
vˇsak tentokr´at zvol´ıme
ako
∞
ψN
= xT
N Plyap xN ,
ui = 0 pre i ≥ N.
(5.37)
Plyap z´ıskame rieˇsen´ım Lyapunovej rovnice
AT Plyap A − Plyap + Qr = 0.
(5.38)
Ked’ˇze pre riadenie h = 0 plat´ı, ˇze ΩM = ℜn bez toho, aby sme museli prid´avat’ zvl´aˇstne
obmedzenia na koncov´
y stav, m´ame zabezpeˇcen´
u stabilitu aj splnenie obmedzen´ı pre i ≥
N.
KAPITOLA 5. PREDIKT´IVNE RIADENIE
5.2
49
MPC - formul´
acia probl´
emu
Pri n´avrhu predikt´ıvneho regul´atora, ktor´
y riadi teplotu v referenˇcn´
ych miestnostiach
ˇ
na jednotliv´
ych blokoch budovy CVUT,
sme vych´adzali jednak z teoretick´
ych poznatkov
pop´ısan´
ych v tejto kapitole, ale predovˇsetk´
ym z popisu samotn´eho riaden´eho syst´emu,
ktor´
y bol uveden´
y v kapitole 2. Naˇsou u
´lohou je navrhn´
ut’ tak´
y predikt´ıvny regul´ator,
ktor´
y bude v prvom rade udrˇziavat’ teplotn´
y komfort v referenˇcn´
ych miestnostiach a
z´aroveˇ
n bude minimalizovat’ mnoˇzstvo energie, ktor´e bude vynaloˇzen´e na vykurovanie.
Kladen´e s´
u nasleduj´
uce konkr´etne poˇziadavky:
• teplota v referenˇcn´
ych miestnostiach v¨aˇcˇsia ako 22◦ C vtedy, ked’ je budova obsaden´a
(pracovn´e dni od 8 : 00 do 18 : 00),
• teplota v referenˇcn´
ych miestnostiach v¨aˇcˇsia ako 19◦ C mimo obdobia, ked’ je budova
obsaden´a (v´ıkendy, sviatky, v noci),
• minimaliz´acia vynaloˇzenej energie (v pr´ıpade vykurovacieho syst´emu pop´ısan´eho v
2 predstavuje vynaloˇzen´
u energiu rozdiel ϑhw − ϑrw ),
• obmedzenie na teplotu topnej vody: 20◦ C ≤ ϑhw ≤ 55◦ C (dan´e maxim´alnou a
minim´alnou teplotou vody, ktor´a je dod´avana v´
ymenn´ıkom),
• obmedzenie na r´
ychlost’ zmeny teploty topnej vody: ∆ϑhw ≤ ±1◦ C/min.
Na to aby, sme mohli dosiahnut’ poˇzadovan´eho v´
ysledku a naplnili vyˇsˇsie uveden´e poˇziadavky, je nutn´e n´ajst’ spr´avny tvar u
´ˇcelovej funkcie vr´atane obmedzen´ı a spr´avneho
nastavenia v´ahov´
ych mat´ıc.
Klasick´
y tvar krit´eria 5.25, kde je v´aˇzeny rozdiel yref − y, je v takomto pr´ıpade nepouˇzitel’n´
y, ked’ˇze ciel’om nie je riadit’ na presne stanoven´
u referenˇcn´
u hodnotu, ale udrˇzat’
teplotu v poˇzadovan´
ych medziach (set range riadenie). V´aˇzenie odch´
ylky od referenˇcnej
hodnoty (rovnako v´aˇzen´e je nedok´
urenie aj prek´
urenie) vedie na situ´aciu, ktor´a je zobrazen´a na obr´azku 5.4 kr´
uˇzkovan´
ym priebehom. Ovel’a ˇziaducejˇs´ı je priebeh, ktor´
y je vyobrazen´
y na obr´azku s kr´ıˇzikmi, ked’ budova poˇcas doby, ked’ nie je obsaden´a (poˇzaduje sa
niˇzˇsia teplota v referenˇcnej miestnosti), len vel’mi pozvol’na vychl´ada a nenast´ava probl´em
s nedok´
uren´ım na zaˇciatku obsadenostn´
ych intervalov.
KAPITOLA 5. PREDIKT´IVNE RIADENIE
50
žiadna penalizácia
vysoká penalizácia
Referencia
Klasické kritérium
Zone control
čas
Obr. 5.4: Porovnanie klasickej kriteri´alnej funkcie ,,zone control“.
Najjednoduchˇs´ı spˆosob, ako tak´eto spr´avanie dosiahnut’, je zaviest’ premenliv´e v´aˇzenie odch´
ylky od referenˇcnej hodnoty - vyˇsˇsia v´aha v dobe, ked’ je budova obsaden´a,
niˇzˇsia v´aha v obdob´ı mimo doby, ked’ je budova obsaden´a. Takto s´ıce uspokojivo obmedz´ıme nedok´
urenie a zbytoˇcne vel’k´
u energiu vynaloˇzen´
u na udrˇzanie teplotn´eho komfortu na zaˇciatku nov´eho obsadenostn´eho intervalu v pr´ıpade, ˇze skok v referencii je len
jednor´azov´
y. Ak vˇsak poˇzadujeme viacn´asobn´
u zmenu referencie a t´
ym aj v´ahovej matice,
ˇ sia pon´
tak´
yto pr´ıstup zlyh´ava [36]. Dalˇ
ukan´a moˇznost’, ako zabezpeˇcit’ udrˇzanie vn´
utornej
teploty v poˇzadovan´
ych medziach, je zaviest’ tvrd´e obmedzenia, a to tak, aby tieto obmedzenia zdola obmedzovali teplotu v miestnosti hodnotou, ktor´a je dan´a referenˇcnou
hodnotou v z´avislosti na dennej dobe (19◦ C resp. 22◦ C). Tento pr´ıstup vˇsak mˆoˇze viest’
na nerieˇsitel’nost’ dan´eho probl´emu. Moˇzn´
ym rieˇsen´ım je zjemnenie t´
ychto obmedzen´ı,
pouˇzitie tzv. ”zone control”(niekedy sa taktieˇz pouˇz´ıva pojem funnel MPC) [30]. Ide o pr´ıstup, ked’ je penalizovan´e to, ˇze trajekt´oria optimalizovanej veliˇciny opust´ı poˇzadovan´
u
oblast’. V naˇsom pr´ıpade je ˇziaduce penalizovat’ pokles teploty v referenˇcnej miestnosti
pod 22◦ C resp.19◦ C. Tento pr´ıstup moˇzno pouˇzit’ aj pri minimaliz´acii vstupnej energie,
ked’ bude penalizovan´
y ak´
ykol’vek kladn´
y rozdiel ϑhw − ϑrw . Na z´aklade vyˇsˇsie uveden´
ych
poˇziadavok a n´asledn´
ych u
´vah budeme formulovat’ nasleduj´
uce krit´erium:
J=
min
a(k),b(k),u(k)
P
−1
X
kQ1 a(k)k22 + kQ2 b(k)k22
(5.39)
k=0
s ohl’adom na:
yref (k) − y1 (k) − a(k) ≤ 0 a ≥ 0,
y2 (k) − b(k) ≤ 0 b ≥ 0,
umin C ≤ u(k) ≤ umax ,
|u(k) − u(k − 1)| ≤ ∆umax ,
(5.40)
KAPITOLA 5. PREDIKT´IVNE RIADENIE
51
kde y1 (k) predstavuje riaden´e veliˇciny (teploty v referenˇcn´
ych miestnostiach ϑrs , ϑrn ),
y2 (k) s´
u rozdielov´e pomocn´e v´
ystupy, ktor´e predstavuj´
u v´aˇzeny rozdiel ϑhw − ϑrw , Q1
a Q2 s´
u v´ahov´e matice prisluˇsn´
ych rozmerov. ∆umax , umax , umin s´
u obmedzenia, ktor´e
odpovedaj´
u poˇziadavk´am uveden´
ym vyˇsˇsie.
Dvojica a(k), b(k) s´
u tzv.”slack variables” [23], ktor´e su nenulov´e iba za predpokladu,
ˇze je poruˇsen´e zadan´e obmedzenie (v pr´ıpade a(k) ked’ klesne teplota v miestnosti pod poˇzadovan´
u, v pr´ıpade b(k) ak je rozdiel ϑhw − ϑrw kladn´
y).
Pre v´
ystupy y1 (k) a y2 (k) plat´ı:
y1 (k) =A
k−1
x(0) +
y2 (k) =Ak−1 x(0) +
k−1
X
i=0
k−1
X
CA(k − 1 − i)(Bu(i)) + Du(k),
(5.41)
C2 A(k − 1 − i)(Bu(i)) + D2 u(k).
i=0
Tak´eto krit´erium moˇzno prep´ısat’ do formy u
´lohy kvadratick´eho programovania 5.15 a
n´asledne na rieˇsenie tejto u
´lohy pouˇzit’ jeden z mnoˇzstva dostupn´
ych solverov pre u
´lohy
kvadratick´eho programovania.
Vol’ba v´ahov´
ych mat´ıc Q1 a Q2 sa poˇcas pouˇz´ıvania uk´azala ako kl’u
´ˇcov´a pre poˇzadovan´e optim´alne spr´avanie MPC. V´
yber v´ahov´
ych mat´ıc bol uˇcinen´
y individu´alne
pre kaˇzd´
y blok, ked’ˇze pre kaˇzd´
y blok je pre udrˇzanie teplotn´eho konfortu potrebn´e in´e
mnoˇzstvo dodanej energie. Hodnoty mat´ıc Q1 a Q2 boli laden´e experimetn´alne na z´aklade
online vykresl’ovania oboch ˇclenov kriteri´alnej funkcie 5.39. Pri vyberan´ı konkr´etnych
hodnˆot v´ahov´
ych mat´ıc bolo dban´e na to, aby bol zachovan´
y ,,rozumn´
y“ pomer medzi
ˇclenmi kriteri´alnej funkcie. Poslednou vol’bou, ktor´
u bolo nutn´e pri nasaden´ı MPC uˇcinit’,
bol v´
yber predikˇcn´eho horizontu, ktor´
y bol napokon zvolen´
y na P = 160, ˇco pri vzorkovacej peri´ode 18 min´
ut predstavuje 2 dni.
Kapitola 6
Simul´
acie na re´
alnych datach
V tejto kapitole uvedieme sadu simul´aci´ı, na ktor´
ych uk´aˇzeme u
´speˇsnost’ jednotliv´
ych
identifikaˇcn´
ych met´od pop´ısan´
ych v predchadzaj´
ucom texte. Identifikaˇcn´e met´ody boli
ˇ
testovan´e na re´alnych d´atach z budovy CVUT
v Prahe. Zhodnotenie v´
ysledkov a simul´aci´ı
v tejto pr´aci uv´adzame len pre blok B1, podobne vˇsak boli identifikovan´e aj ostatn´e bloky
budovy Fakulty elektrotechnickej a Fakulty strojnej a vytvoren´e modely s´
u pouˇz´ıvan´e
pre MPC, ktor´eho u
´speˇsnost’ a porovnanie s konvenˇcn´
ymi met´odami riadenia (ekvitermn´a
regul´acia) bude diskutovan´e v druhej ˇcasti tejto kapitoly.
6.1
Identifik´
acia
Pri identifik´acii bola dostupn´a datab´aza historick´
ych d´at, ktor´a obsahovala u
´daje o teplote v miestnosti, teplote vykurovacej vody, teplote vratnej vody a arch´ıvne predpovede
vonkajˇsej teploty. Ked’ˇze data v datab´aze boli uloˇzen´e v nevhodnom form´ate, pred t´
ym,
neˇz boli data pouˇzit´e pri identifik´acii, museli byt’ upraven´e sadou predspracovac´ıch skriptov, ktor´e boli implementovan´e v programovom prostred´ı Scilab1 . D´ata boli n´asledne
prevzorkovan´e na vzorkovaciu peri´odu Ts = 18 min. Vzorkovacia peri´oda bola vybrat´a
s reˇspektovan´ım dynamiky budovy. Ked’ˇze zvolen´e v´
ystupy (2.1) maj´
u kaˇzd´
y rozliˇcne
r´
ychlu dynamiku, vol’ba 18 min´
utovej vzorkovacej peri´ody predstavuje kompromis medzi
vzorkovan´ım pomal´eho v´
ystupy teploty v miestnosti a relat´ıvne r´
ychlej teploty vratnej
vody [48]. Na identifik´aciu a testovanie vˇsetk´
ych uveden´
ych identifikaˇcn´
ych met´od sme
pouˇzili rovnak´
y u
´sek d´at z janu´ara 2011 a vˇsetky z´ıskan´e modely sme testovali na va1
http://www.scilab.org/
52
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
53
lidaˇcn´
ych d´atach z febru´ara 2011. Teraz uvedieme zoznam testovan´
ych identifikaˇcn´
ych
met´od, vr´atane znaˇcenia, ktor´e pre ne budeme v d’al’ˇsom texte pouˇz´ıvat’:
• 4SID- Subspace identifikaˇcn´e met´ody, ktor´
ych struˇcn´
y princ´ıp sme uviedli v 3.2,
• GB - identifik´acia zaloˇzen´a na Grey box modelovan´ı pop´ısan´a v podkapitole 3.1,
• MRI - identifikaˇcn´a met´oda zaloˇzen´a na minimaliz´acii viackrokovej predikˇcnej
chyby, pouˇzit´
y bol algoritmus 4.3,
• MRI+GB- identifikaˇcn´a met´oda zaloˇzen´a na minimaliz´acii viackrokovej predikˇcnej
chyby, ktor´a reˇspektuje fyzik´alnu ˇstrukt´
uru identifikovan´eho syst´emu,
• MRI+PLS - kombin´acia ˇciastoˇcnej line´arnej regresie a MRI identifik´acie podl’a
popisu uveden´eho v 4.4.
Dodajme eˇste, ˇze hl’ad´ame model, ktor´
y bude pouˇzit´
y pre predikt´ıvny regul´ator, preto
n´as bude zauj´ımat’ predovˇsetk´
ym to, ako dobre dok´aˇze tento model predikovat’ teploty
v referenˇcn´
ych miestnostiach poˇcas dvojdˇ
nov´eho predikˇcn´eho horizontu. Preto nebudeme
pri simul´aciach vykresl’ovat’ (ako je zvykom) simulovan´e v´
ystupy modelov, ale k-krokov´e
predikcie teplˆot v referenˇcn´
ych miestnostiach, konkr´etne predikcie na dva dni dopredu,
ˇco odpoved´a 160 krokov´emu predikˇcn´emu horizontu. Pre porovnavanie kvality modelov
z´ıskan´
ych jednotliv´
ymi met´odami nebud´
u pouˇzit´e len predikˇcn´e vlastnosti jednotliv´
ych
modelov, ale aj skutoˇcnost’, ˇci je dan´
y model vhodn´
y pre riadenie, ˇco moˇzno otestovat’
odozvami t´
ychto modelov na niektor´e typick´e vstupn´e sign´aly.
6.1.1
Subspace identifikaˇ
cn´
e met´
ody
Pri tvorbe modelov za pouˇzitia 4SID sme pouˇzili implement´aciu, ktor´a je s´
uˇcast’ou Identification toolbox pre Matlab [28]. Na z´aklade anal´
yzy singul´arnych ˇc´ısel matice, ktor´a
vznikla projekciou priestoru bud´
ucich v´
ystupov pozd´lˇz priestoru bud´
ucich vstupov do riadkov´eho priestoru minul´
ych d´at, bol odhadnut´
y r´ad syst´emu stanoven´
y na 4. Na obr´azku
6.1 vid´ıme predikcie teploty v referenˇcn´
ych miestnostiach na dve hodiny a dva dni. K´
ym
na p´ar hod´ın dok´aˇze model z´ıskan´
y Subspace identifikaˇcn´
ymi met´odami predikovat’ teplotu v miestnosti vel’mi presne, na dva dni s´
u uˇz predikcie tohto modelu omnoho nepresnejˇsie.
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
160 − kroková predikcia (juh)
160 − kroková predikcia (sever)
24
24
namerané data
predikcia
namerané data
predikcia
23.5
23.5
23
23
22.5
22.5
21.5
rn
rs
22
ϑ (°C)
22
ϑ (°C)
54
21.5
21
21
20.5
20.5
20
20
19.5
19.5
19
14.1. 0:00 16.1. 0:00 18.1. 0:00 20.1. 0:00 22.1. 0:00
19
14.1 0:00
16.1 0:00
18.1 0:00
20.1 0:00
22.1 0:00
Obr. 6.1: Predikcie teplˆot v referenˇcn´
ych miestnostiach (4SID).
ϑ(°C)
0.6
ϑo→ϑrs
0.4
ϑo→ϑrws
0.2
ϑ →ϑ
o
−0.2
rn
ϑo→ϑrwn
0
0
48
96
144
192
t(h)
ϑ(°C)
2
ϑhws→ϑrs
ϑhws→ϑrws
0
ϑhws→ϑrn
ϑhws→ϑrwn
−2
−4
0
48
96
144
192
t(h)
ϑ(°C)
4
ϑhwn→ϑrs
ϑ
2
→ϑ
hwn
rws
ϑhwn→ϑrn
0
ϑ
→ϑ
hwn
−2
0
48
96
144
rwn
192
t(h)
Obr. 6.2: Odozvy na jednotkov´
y skok na vstupoch (4SID).
6.1.2
Grey box modelovanie
ˇ s´ı z identifikaˇcn´
Dalˇ
ych pr´ıstupov, ktor´
y sme pouˇzili pri identifik´acii, bol postup, ktor´
y
sme podrobne pop´ısali v podkapitole 3.1. Na odhad parametrov mat´ıc stavov´eho popisu
A, B s poˇzadovanou ˇstrukt´
urou 3.10 sme pouˇzili skript pre Grey box identifik´aciu, ktor´
y
sme implementovali v programovom prostred´ı Scilab2 .
2
http://www.scilab.org/
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
160 − kroková predikcia (juh)
160 − kroková predikcia (sever)
24
24
namerané data
predikcia
namerané data
predikcia
23.5
23.5
23
23
22.5
22.5
21.5
rn
rs
22
ϑ (°C)
22
ϑ (°C)
55
21.5
21
21
20.5
20.5
20
20
19.5
19.5
19
14.1.0:00
16.1.0:00
18.1.0:00
20.1.0:00
19
14.1 0:00
22.1.0:00
16.1 0:00
18.1 0:00
20.1 0:00
22.1 0:00
Obr. 6.3: Predikcie teplˆot v referenˇcn´
ych miestnostiach (GB).
0.4
ϑ →ϑ
ϑ(°C)
o
rs
ϑ →ϑ
0.2
o
rws
ϑ →ϑ
o
−0.2
rn
ϑ →ϑ
0
o
0
48
96
144
rwn
192
t(h)
1
ϑ
→ϑ
ϑ
→ϑ
ϑ(°C)
hws
0.5
rs
hws
rws
ϑhws→ϑrn
ϑhws→ϑrwn
0
−0.5
0
48
96
144
192
t(h)
1
ϑ
→ϑ
ϑ
→ϑ
ϑ(°C)
hwn
0.5
hwn
rs
rws
ϑhwn→ϑrn
ϑhwn→ϑrwn
0
−0.5
0
48
96
144
192
t(h)
Obr. 6.4: Odozvy na jednotkov´
y skok na vstupoch (GB).
Model z´ıskan´
y Grey box identifik´aciou vykazuje podobne ako model z´ıskan´
y Subspace
identifik´aciami dobr´e predikˇcn´e spr´avanie na p´ar hod´ın dopredu, na dva dni dopredu vˇsak
uˇz nedok´aˇze odhadn´
ut’ teplotu v miestnosti presne. Narozdiel od modelu z´ıskan´eho 4SID
met´odami m´a vˇsak jednu v´
yhodu. Vd’aka tomu, ˇze zohl’adˇ
nuje ˇstrukt´
uru, skokov´e odozvy
6.4 si zachov´avaj´
u fyzik´alny zmysel.
6.1.3
MRI identifikaˇ
cn´
e met´
ody
Z dvojice uveden´
ych met´od minimalizuj´
ucich viackrokov´
u predikˇcn´
u chybu sme si vybrali
met´odu, ktor´
u sme oznaˇcili ako jednokrokov´
y algoritmus a bola prezentovan´a D. Laur´ım
a M. Mart´ınezom v [25]. Testovanie dvojkrokov´eho algoritmu na pr´ıklade 4.1 uk´azalo, ˇze
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
56
kvalita z´ıskan´eho odhadu je znaˇcne podmienen´a znalost’ou ˇstrukt´
ury modelu ˇsumu, ktor´
u
v pr´ıpade re´alneho syst´emu nem´ame k dispoz´ıcii. Algoritmus je d’alej vel’mi citliv´
y na vybudenost’. Ked’ˇze d´ata, ktor´e sme mali k dispoz´ıcii, poch´adzaju z uzavretej sluˇcky, boli
zat’aˇzen´e znaˇcnou korel´aciou (aktu´alny vstup - minul´e v´
ystupy), ˇco znemoˇznilo pouˇzitie
tejto met´ody, a preto bol pouˇzit´
y len algoritmus pop´ısany v 4.3. K implement´acii boli
pouˇzit´e solvery na rieˇsenie u
´loh neline´arnej optimaliz´acie, ktor´e s´
u s´
uˇcast’ou Optimization
toolbox pre Matlab [6]. Pri identifik´acii boli zvolen´e parametre ˇstrukt´
ury modelu na = 1,
nb = 1, nk = 1. V pr´ıpade, ˇze stavov´e premenn´e odpovedaj´
u priamo meran´
ym v´
ystupom
(ˇco je aj n´aˇs pr´ıpad), je jednoduch´e n´ajst’ s´
uvislost’ medzi ˇstrukt´
urou Θ (4.41) z´ıskan´
u
MRI identifik´aciou a maticami stavov´eho popisu A, B
" #
BT
Θ=
.
AT
(6.1)
D´lˇzka predikˇcn´eho horizontu bola zvolen´a P = 100, hoci predikˇcn´
y horizont MPC,
pre ktor´
y sa model pouˇzije, je P = 160. Ako sme vˇsak uˇz uk´azali v kapitole o MRI
identifikaˇcn´
ych met´odach, kvalita viackrokov´
ych predikci´ı rastie zo zvyˇsuj´
ucim horizontom len do urˇcitej hodnoty. Pre vyˇsˇs´ı horizont ako P = 100 sa kvalita identifikovan´eho
modelu uˇz nijak nezvyˇsuje, iba zbytoˇcne rastie v´
ypoˇctov´a n´aroˇcnost’.
160 − kroková predikcia (juh)
160 − kroková predikcia (sever)
24
24
namerané data
predikcia
namerané data
predikcia
23.5
23
23
22.5
22.5
22
22
ϑ (°C)
21.5
rn
rs
ϑ (°C)
23.5
21.5
21
21
20.5
20.5
20
20
19.5
19.5
19
14.1. 0:00 16.1. 0:00 18.1. 0:00
20. 0:00
22.1. 0:00
19
14.1 0:00
16.1 0:00
18.1 0:00
20.1 0:00
22.1 0:00
Obr. 6.5: Predikcie teplˆot v referenˇcn´
ych miestnostiach (MRI).
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
0.6
57
ϑ →ϑ
ϑ(°C)
o
rs
0.4
ϑ →ϑ
0.2
ϑ →ϑ
o
rws
o
rn
ϑ →ϑ
o
0
−0.2
0
48
96
144
rwn
192
t(h)
ϑ(°C)
1
ϑhws→ϑrs
ϑ
→ϑ
ϑ
→ϑ
ϑ
→ϑ
hws
0.5
rws
hws
0
−0.5
rn
hws
0
48
96
144
rwn
192
t(h)
ϑ(°C)
1
ϑhwn→ϑrs
ϑhwn→ϑrws
0.5
ϑhwn→ϑrn
ϑ
0
−0.5
→ϑ
hwn
0
48
96
144
rwn
192
t(h)
Obr. 6.6: Odozvy na jednotkov´
y skok na vstupoch (MRI).
6.1.4
Grey box modelovanie v kombin´
acii s MRI
T´ato met´oda taktieˇz pouˇziv´a k odhadu mat´ıc stavov´eho popisu MRI identifikaˇcn´
y algoritmus, ale narozdiel od predoˇslej met´ody dodrˇzuje fyzik´alnu podstatu syst´emu a identifikuje matice v ˇziadanej ˇstrukt´
ure 3.10, na ˇco je pouˇzit´
y solver na rieˇsenie neline´arnych
optimalizaˇcn´
ych u
´loh s obmedzen´ım.
Parametre sme volili podobne ako v pr´ıpade klasickej MRI identifik´acie: na = 1,
nb = 1, nk = 1, P = 100.
160 − kroková predikcia (juh)
160 − kroková predikcia (sever)
24
24
namerané data
predikcia
namerané data
predikcia
23.5
23.5
23
23
22.5
22.5
22
ϑ (°C)
rn
rs
ϑ (°C)
22
21.5
21.5
21
21
20.5
20.5
20
20
19.5
19.5
19
14.1.0:00
16.1.0:00
18.1.0:00
20.1.0:00
22.1.0:00
19
14.1.0:00
16.1.0:00 18.1. 0:00 20.1 0:00 22.1. 0:00
Obr. 6.7: Predikcie teplˆot v referenˇcn´
ych miestnostiach (MRI+GB).
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
0.6
58
ϑ →ϑ
o
rs
0.4
ϑ →ϑ
0.2
ϑ →ϑ
ϑ(°C)
o
rws
o
rn
ϑ →ϑ
o
0
−0.2
0
48
96
144
rwn
192
t(h)
1
ϑ
→ϑ
ϑ
→ϑ
ϑ
→ϑ
ϑ
→ϑ
ϑ(°C)
hws
0.5
hws
hws
0
−0.5
hws
0
48
96
144
rs
rws
rn
rwn
192
t(h)
ϑ(°C)
1
ϑhwn→ϑrs
ϑhwn→ϑrws
0.5
ϑhwn→ϑrn
0
−0.5
ϑhwn→ϑrwn
0
48
96
144
192
t(h)
Obr. 6.8: Odozvy na jednotkov´
y skok na vstupoch (MRI+GB).
6.1.5
Kombin´
acia MRI a PLS
Posledn´a uveden´a identifikaˇcn´a met´oda vyuˇz´ıva kombin´aciu MRI a PLS tak, ako bola
pop´ısan´a v [25]. Algoritmus bol implementovan´
y v programe Matlab a pri identifik´acii
boli op¨at’ zvolen´e parametre na = 1, nb = 1, nk = 1, P = 100. Tento algoritmus vyˇzaduje
okrem t´
ychto parametrov aj poˇcet komponent, ktor´e s´
u pouˇzit´e pri identifik´acii. V naˇsom
pr´ıpade bol poˇcet parametrov zvolen´
y na 6 zo 7 na z´aklade anal´
yzy vz´ajomnej z´avislosti
v´
ystupov a regresora.
160 − kroková predikcia (juh)
160 − kroková predikcia (sever)
24
24
namerané data
predikcia
namerané data
predikcia
23.5
23
23
22.5
22.5
22
22
ϑ(°C)
ϑ(°C)
23.5
21.5
21.5
21
21
20.5
20.5
20
20
19.5
19.5
19
14.1.0:00
16.1.0:00
18.1.0:00
20.1.0:00
22.1.0:00
19
14.1.0:00
16.1.0:00
18.1.0:00
20.1.0:00
22.1.0:00
Obr. 6.9: Predikcie teplˆot v referenˇcn´
ych miestnostiach (MRI+PLS).
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
59
0.6
ϑ →ϑ
ϑ(°C)
o
rs
0.4
ϑ →ϑ
0.2
ϑ →ϑ
o
rws
o
rn
ϑ →ϑ
o
0
−0.2
0
48
96
144
rwn
192
t(h)
1
ϑ
→ϑ
ϑ
→ϑ
ϑ
→ϑ
ϑ
→ϑ
ϑ(°C)
hws
0
rs
hws
rws
hws
−1
−2
rn
hws
0
48
96
144
rwn
192
t(h)
2
ϑ
→ϑ
ϑ
→ϑ
ϑ
→ϑ
ϑ
→ϑ
ϑ(°C)
hwn
hwn
1
hwn
0
−1
hwn
0
48
96
144
rs
rws
rn
rwn
192
t(h)
Obr. 6.10: Odozvy na jednotkov´
y skok na vstupoch (MRI+PLS).
6.1.6
Porovnanie identifikovan´
ych modelov
V tejto ˇcasti uvedieme porovanie predikˇcn´eho spr´avania modelov z´ıskan´
ych uvaˇzovan´
ymi
identifikaˇcn´
ymi met´odami. Pre porovnanie t´
ychto modelov sme pouˇzili rozptyl viackrokov´
ych predikˇcn´
ych ch´
yb vars/n (i) a nasleduj´
uci fit faktor:


ˆ
ϑ
(k
+
i)
−
ϑ
(k
+
i|k)
rs/n
rs/n
2  100%,
f its/n (i) = 1 − ϑrs/n (k + i) − E(ϑrs/n (k + i))
2
(6.2)
kde i ∈ {1, 2, . . . , 160} a ϑrs/n predstavuje teplotu v juˇznej/severnej referenˇcnej miestnosti. Na nasleduj´
ucich obr´azkoch 6.12, 6.11 moˇzno vidiet’ vykreslen´e rozptyly viackrokov´
ych predikˇcn´
ych ch´
yb a fit faktory z´ıskan´
ych modelov v z´avislosti na meniacom sa
poˇcte predikˇcn´
ych krokov. V pr´ıpade predikcie teploty na p´ar krokov vid´ıme, ˇze vˇsetky
z´ıskan´e modely najm¨a v pr´ıpade severnej miestnosti predikuj´
u pribliˇzne rovnako presne.
Pre vyˇsˇsie predikˇcn´e horizonty vˇsak viditel’ne lepˇsie spr´avanie vykazuj´
u modely z´ıskan´e
met´odami, ktor´e minimalizuj´
u viackrokov´
u predikˇcn´
u chybu. Najlepˇsie v´
ysledky vykazuje kombin´acia MRI+PLS, ktor´a vd’aka postupu, ktor´
y sme pop´ısali v ˇcasti 4.4, dok´aˇze
pri identifik´acii vyuˇz´ıvat’ len d´ata, ktor´e v sebe nes´
u releventn´
u inform´aciu, potlaˇcit’ vplyv
ˇsumu a t´
ym poskytn´
ut’ model s uspokojov´
ymi predikˇcn´
ymi vlastnost’ami. Len o nieˇco nepresnejˇsie predikcie (v pr´ıpade severnej miestnosti dokonca niekedy aj o nieˇco presnejˇsie)
poskytuje model z´ıskan´
y klasickou MRI identifik´aciou.
Oba tieto modely vˇsak maj´
u jeden z´asadny nedostatok, ktor´
y moˇzno pozorovat’ na
obr´azkoch 6.5 a 6.10, kde s´
u vykreslen´e odozvy na jednotkov´
y skok t´
ychto modelov.
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
60
Zar´aˇzaj´
uce s´
u najm¨a odozvy na vstup vykurovacej vody v juˇznej vetve. Oba modely
totiˇz predpovedaj´
u, ˇze zmena vykurovacej vody na vstupe z 0 na 1◦ C spˆosob´ı pokles
teploty nielen vratnej vody ale aj teploty v oboch miestnostiach. Tak´ato situ´acia samozrejme neodpoved´a skutoˇcnosti tak, ako bola pop´ısana v kapitole 2, ˇco je spˆosoben´e
predovˇsetk´
ym t´
ym, ˇze oba modely s´
u z´ıskan´e pr´ıstupom zaloˇzen´
ym na Black box modelovan´ı a nereˇspektuj´
u teda fyzik´alnu podstatu syst´emu. Pouˇzitie tak´ehoto modelu v kombin´acii s MPC by mohlo viest’ na nepr´ıjemn´e situ´acie, ako napr´ıklad snaha zv´
yˇsit’ teplotu
v miestnosti zn´ıˇzen´ım teploty vykurovacej vody a podobne.
Model z´ıskan´
y minimaliz´aciou viackrokovej predikˇcnej chyby, ktor´
y reˇspektuje ˇstrukt´
uru riaden´eho procesu (MRI+GB), predpoved´a s´ıce nepresnejˇsie ako predoˇsl´e dva uveden´e modely, ale pri pohl’ade na jeho skokov´e odozvy 6.8 vid´ıme, ˇze odpoved´a pribliˇzne
tomu, ˇco by sme oˇcakav´ali na z´aklade popisu riaden´eho syst´emu uveden´eho v kapitole 2. Schopnost’ predikcie zostavaj´
ucich dvoch modelov (4SID, GB) je znaˇcne horˇsia
ako u met´od minimalizuj´
ucich viackrokov´
u chybu, ˇco moˇzno vidiet’ aj na 6.1 a 6.3. Model
z´ıskan´
y 4SID nanajv´
yˇs d´ava odozvy 6.2, ktor´e nie v s´
ulade s dynamikou budovy.
Na z´aklade simul´aci´ı a u
´vah, ktor´e sme uviedli vyˇsˇsie, sa ako najvhodnejˇs´ı pre n´asledn´e
pouˇzitie v kombin´acii s MPC jav´ı model z´ıskan´
y met´odou, ktor´
u sme oznaˇcili MRI+GB.
Poskytuje uspokojiv´e predikcie aj na v¨aˇcˇsie horizonty a podl’a skokov´
ych odoziev odpoved´a dynamike riadenej budovy. Modely z´ıskan´e touto met´odou boli pouˇzit´e aj pri identifik´aci ostatn´
ych blokov a n´asledne boli rovnako ako na nami uvaˇzovanom bloku B1
ˇ
pouˇzit´e pre MPC, ktor´e pln´ı v budove CVUT
v Dejviciach u
´lohu regul´atora na vyˇsˇsej
u
´rovni riadenia.
Vˇsetky uvaˇzovan´e met´ody, ktor´e nebrali do u
´vahy ˇstrukt´
uru riaden´eho syst´emu, sa
javili pri s´
uˇcasnom poˇcte meran´
ych vstupov a v´
ystupov ako nepouˇzitel’n´e pre n´asledn´e
pouˇzitie s MPC. Za hlavn´
u pr´ıˇcinu tak´ehoto spr´avania mˆoˇzeme povaˇzovat’ fakt, ˇze na budovu pˆosobia okrem vstupov, ktor´e uvaˇzujeme v naˇsom modeli, aj in´e nemeratel’n´e vstupy.
Vel’k´
y vplyv na teplotu (predovˇsetk´
ym v juˇzne orientovan´
ym miestnostiach) m´a slneˇcn´e
ˇziarenie, ˇco je badatel’n´e aj z faktu, ˇze predikcie teploty v severne orientovan´
ych miestnostiach s´
u v´
yrazne presnejˇsie.
Nezahrnut´
ym vplyvom intenzity slneˇcn´eho ˇziarenia si moˇzno vel’mi jednoducho vysvetlit’ aj nezmyseln´
y priebeh skokov´
ych odoziev 6.2, 6.10, 6.6 u Black box modelov,
podl’a ktor´
ych by n´asledkom zv´
yˇsenia teploty vykurovacej vody v juˇznom vykurovacom
okruhu mal byt’ pokles teplˆot v jednotliv´
ych miestnostiach. N´azornejˇsie sa tento paradox
d´a pop´ısat’ tak, ˇze pokles teploty vykurovacej vody spˆosob´ı n´arast teploty v miestnostiach. Skutoˇcnost’ je vˇsak in´a. D´ata, ktor´e boli pouˇzit´e pri identifik´acii, s´
u znaˇcne korelo-
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
61
van´e, ˇco je v´
ysledok sp¨atnej v¨azby od regul´atora. Z toho dˆovodu je teplota vykurovacej
vody nepriamo u
´merne z´avisl´a na intenzite slneˇcn´eho ˇziarenia (sp¨atnov¨azbov´
y regul´ator
na zv´
yˇsenie intenzity reaguje zn´ıˇzen´ım teploty vykurovacej vody). Pokles teploty vykurovacej vody v skutoˇcnosti len odr´aˇza zv´
yˇsenie intenzity slneˇcn´eho ˇziarenia, ˇco samozrejme
spˆosob´ı zv´
yˇsenie vn´
utornej teploty.
ˇ s´ım vstupom, ktor´
Dalˇ
y nemoˇzno len tak zanedbat’, je mnoˇzstvo l’ud´ı, ktor´e sa v dan´
u
dobu v budove (miestnosti) nach´adza (tzv. occupancy). V pr´ıpade rozˇs´ırenia modelu
o tieto vstupy mˆoˇzeme oˇcak´avat’, ˇze aj modely z´ıskan´e met´odami zaloˇzen´
ymi na Black
box modelovan´ı sa stan´
u vhodn´
ymi pre nasadenie s MPC. Podrobnejˇsie budeme vplyv
ˇ
slneˇcn´eho osvitu na spr´avanie budovy CVUT
a moˇznosti d’alˇsieho vylepˇsenia modelu
diskutovat’ v nasleduj´
ucej podkapitole.
juh
0.4
s
var (i)
0.3
MRI+GB
MRI+PLS
MRI
GB
4SID
0.2
0.1
0
20
40
60
80
i
100
120
140
160
100
120
140
160
sever
0.4
n
var (i)
0.3
MRI+GB
MRI+PLS
MRI
GB
4SID
0.2
0.1
0
20
40
60
80
i
Obr. 6.11: Porovnanie viackrokov´
ych predikˇcn´
ych ch´
yb identifikovan´
ych
modelov.
juh
100
MRI+GB
MRI+PLS
MRI
GB
4SID
60
s
fit (i)(%)
80
40
20
20
40
60
80
i
100
120
140
160
100
120
140
160
sever
100
MRI+GB
MRI+PLS
MRI
GB
4SID
60
n
fit (i)(%)
80
40
20
20
40
60
80
i
Obr. 6.12: Porovnanie fitfaktorov identifikovan´
ych modelov.
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
6.2
62
Predikt´ıvny regul´
ator
Predikt´ıvny regul´ator, ktor´eho n´avrh sme pop´ısali podrobne v 5.2, je pouˇz´ıvan´
y pri riaden´ı
teploty v referenˇcn´
ych miestnostiach na vˇsetk´
ych blokoch budovy Fakulty elektrotechnickej a Fakulty strojnej. V tejto kapitole uvedieme porovnie spom´ınan´eho predikt´ıvneho
regul´atora a klasickej ekvitermnej regul´acie, ktor´a bola pouˇz´ıvan´a na nami uvaˇzovanej
ˇ
budove CVUT
pred nasaden´ım predikt´ıvnej regul´acie.
K tomu, aby sme mohli porovnat’ oba spˆosoby regul´acie aj napriek tomu, ˇze oba regul´atory pracovali za rˆoznych podmienok (predovˇsetk´
ym in´e poveternostn´e podmienky),
zavedieme pojmy ECM (Energy Consumption) a HDD(Heating Degree Day) nasledovne
HDD =
T
end
X
yref (i) − ϑo (i)
(6.3)
(ϑhws (i) − ϑrws (i)) + (ϑhwn (i) − ϑrwn (i)) ,
(6.4)
i=Tstart
a
ECM =
T
end
X
i=Tstart
kde Tstart predstavuje ˇcas zaˇciatku experimentu a Tend jeho koniec. Pre porovnanie kvality
pri rˆoznych poveternostn´
ych podmienkach pouˇzijeme pomer
ECM
HDD
vyjadruj´
uci energiu
spotrebovan´
u na vykurovanie severnej aj juˇznej vetvy 6.4 vztiahnut´
u na rozdiel medzi
poˇzadovanou referenˇcnou teplotou a vonkajˇsou teplotou 6.3.
Najprv sme porovnali regul´aciu ekvitermou, ktor´a bola nasaden´a na bloku B1 na zaˇciatku roku 2010, s predikt´ıvnym regul´atorom zo zaˇciatku roku 2011. Na obr´azkoch
6.13,6.14 vid´ıme vykreslen´e integr´aly rozdielu ϑhw − ϑrw , ktor´e predstavuj´
u energiu spotrebovan´
u na vykurovanie pre ekvitermn´
u regul´aciu aj MPC. Uˇz zo samotn´
ych obr´azkov je
badatel’n´e, ˇze integr´al dan´eho rozdielu m´a urˇcite menˇsiu plochu v pr´ıpade predikt´ıvneho
regul´atora ako v pr´ıpade klasickej ekvitermnej regul´acie, ˇco potvrdzuje fakt, ˇze MPC
moˇzno v pr´ıpade znalosti matematick´eho modelu a vhodne zvolenej u
´ˇcelovej funkcie s v´
yhodou pouˇzit’ pri riaden´ı vel’k´
ych budov a uˇsetrit’ tak znaˇcn´e mnoˇzstvo energie.
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
63
ekvitermná regulácia
50
← ϑhws
ϑ(°C)
40
30
←ϑ
rws
20
5.1.2010
9.1.2010
13.1.2010
17.1.2010
21.1.2010
25.1.2010
29.1.2010
MPC
50
ϑ(°C)
40
←ϑ
hws
30
←ϑ
20
5.1.2011
rws
9.1.2011
13.1.2011
17.1.2011
21.1.2011
25.1.2011
29.1.2011
Obr. 6.13: Porovnanie ekvitermnej regul´
acie a MPC (juh).
ekvitermná regulácia
55
50
45
← ϑ
hwn
ϑ(°C)
40
35
30
25
20
15
5.1.2010
← ϑrwn
9.1.2010
13.1.2010
17.1.2010
21.1.2010
25.1.2010
29.1.2010
MPC
55
50
45
ϑ(°C)
40
←ϑ
hwn
35
30
25
← ϑ
rwn
20
15
5.1.2011
9.1.2011
13.1.2011
17.1.2011
21.1.2011
25.1.2011
29.1.2011
Obr. 6.14: Porovnanie ekvitermnej regul´
acie a MPC (sever).
Ked’ˇze vˇsak pri porovnan´avani pouˇzit´
ych strat´egii neboli rovnak´e poveternostn´e podmienky, na potvrdenie u
´spor, ktor´e boli dosiahnut´e pomocou MPC, uv´adzame aj tabul’ku
6.1, ktor´a taktieˇz potvrdzuje energetick´
u v´
yhodnost’ MPC a jeho relat´ıvne u
´spory oproti
konvenˇcn´emu ekvitermn´emu riadeniu dokonca aˇz o 27%.
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
64
Tabul’ka 6.1: Porovnanie u
´spor.
6.3
ECM /HDD
poˇcet dn´ı relat´ıvne u
´spory
MPC
0.6012
59
EKVITERMA
0.8244
84
27.1%
Zahrnutie vplyvu slnka
Ani jeden z modelov, ktor´
y sme doteraz uvaˇzovali, nemal v sebe zahrnut´
y vplyv slneˇcn´eho
ˇziarenia, ktor´e m´a vˇsak urˇcite nezanedbatel’n´
y vplyv na spr´avanie budovy, a to najm¨a
v pr´ıpade juˇzne orientovan´
ych miestnosti, ˇco potvrdili aj simul´acie v sekcii 6.1, kde bolo
badatel’n´e, ˇze vˇsetky z´ıskan´e modely predikovali jednoznaˇcne presnejˇsie teplotu v juˇzne
orientovanej miestnosti.
Na to, aby sme vˇsak mohli zahrn´
ut’ do modelu pre MPC vplyv slneˇcn´eho ˇziarenia, by
sme potrebovali jeho dostatoˇcne presn´e predpovede na dobu horizontu predikcie. Ked’ˇze
vˇsak tak´eto predpovede neboli dostupn´e, uv´adzame v naˇsej pr´aci len anal´
yzu vplyvu
slneˇcn´eho ˇziarenia.
Ako vstup sme pre tieto anal´
yzy pouˇzili merania slneˇcn´eho osvitu nameran´e fotobunkou Siemens QLS60, ktor´a bola umiestnen´a na juˇznej strane budovy 6.15.
Obr. 6.15: Umiestnenie ˇcidla intenzity slneˇcn´eho ˇziarenia.
Zahnrut´ım vplyvu slneˇcn´eho ˇziarenia do modelu sa pˆovodn´
y vektor vstupov rozˇs´ıril
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
na:
65
h
iT
u = ϑo ϑhws ϑhwn I .
(6.5)
ˇ
kde I(Wm−2 ) predstavuje intenzitu slneˇcn´eho ˇziarenia. Strukt´
ura syst´emov´
ych mat´ıc A,
B sa z pˆovodne ˇstrukt´
ury (3.10) rozˇs´ırila na:




b1 0 0 b5
a1 a4 a3 0




 0 b3 0 0 
 a2 a5 0 0 




A=
 , B = b 0 0 b  .
a
0
a
a
 2
 3
6
6
8
0 0 b4 0
0 0 a7 a9
(6.6)
Pri samotnej identifik´acii vplyvu slneˇcn´eho ˇziarenia sme sa museli vysporiadat’ s jedn´
ym
z´asadn´
ym probl´em, ktor´
y bol spomenut´
y aj vyˇsˇsie. D´ata, ktor´e boli dostupn´e pri identifik´acii, poch´adzali bud’ z obdobia mimo vykurovaciu sez´onu, ked’ sa netopilo, a k identifik´acii boli nepouˇzitel’n´e kvˆoli nedostatoˇcnej vybudenosti. Druhou skupinou boli d´ata,
ktor´e boli zozbieran´e poˇcas vykurovacej sez´ony. V tomto obdob´ı vˇsak je teplota v miestnosti riaden´a sp¨atnov¨azbov´
ym regul´atorom, ktor´
y u
´speˇsne potl´aˇca vplyv poruchov´
ych
veliˇcin, medzi ktor´e patr´ı intenzita slneˇcn´eho ˇziarenia a vonkajˇsia teplota, a znemoˇznuje
spr´avne identifikovat’ vplyv t´
ychto vstupov na spr´avanie budovy. To potvrdzuje aj obr´azok
korel´acie medzi teplotou v miestnosti a intenzitou slneˇcn´eho ˇziarenia 6.16.
24
ϑ
rs
ϑ
23.5
rn
23
ϑrs/n(k+1)(°C)
22.5
22
21.5
21
20.5
20
19.5
0
100
200
300
400
500
600
700
800
900
−2
I(k)(Wm )
Obr. 6.16: Korel´acia medzi ϑrs/n a I.
Na rieˇsenie tohto probl´emu sme pouˇzili alternat´ıvny spˆosob identifik´acie pozostavaj´
uci
z dvoch krokov.
1. Identifik´
acia z nekorelovan´
ych d´
at
V prvom kroku algoritmu sme pouˇzili k identifik´acii data, ktor´e boli zozbieran´e
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
66
mimo vykurovaciu sez´onu, t.j. v obdob´ı, kedy sa netopilo a vplyv intenzity slneˇcn´eho
ˇziarenia na teplotu v miestnosti nebol potl´aˇcan´
y regul´atorom. Ked’ˇze teplota vykurovacej vody je v tomto obdob´ı konˇstantn´a, moˇzno identifikovat’ syst´em v nasleduj´
ucom tvare

b1 b5


0 0

x(k + 1) = Ax(k) + 
b b  d(k) + aff,
 2 6
0
h
iT

(6.7)
0
. Aby sme mohli teraz identifikovat’ maticu A a koeficienty
b1 , b2 , b5 , b6 , potrebujeme sa zbavit’ ˇclenu aff predstavuj´
uceho vplyv vykurovacej
vody, ktor´a je po cel´
y ˇcas konˇstant´a. Vd’aka linearite 6.8 toho moˇzno dosiahnut’
kde d =
ϑo I
identifik´aciou z rozdielov´
ych d´at xr a ur :
xr = xa − xb ,
ur = d a − d b ,
(6.8)
kde xa , da a xb , db predstavuj´
u s´
uvisl´e u
´seky d´at nameran´e mimo vykurovacej sez´ony.
2. Identifik´
acia s obmedzen´ım
V druhom kroku sme uskutoˇcnili klasick´
u identifik´aciu mat´ıc so ˇstrukt´
urou (6.6)
z d´at, ktor´e boli z´ıskan´e behom vykurovacej sez´ony. Aby vˇsak bola zabezpeˇcen´a
presn´a identifik´acia parametrov, ktor´e odpovedaj´
u vplyvu poruchov´
ych veliˇc´ın (intenzity slneˇcn´eho ˇziarenia a vonkajˇsej teploty), parametre b1 , b2 , b5 , b6 z´ıskan´e v prvom kroku identifik´acie posl´
uˇzili ako obmedzenia v druhom kroku, tieto parametre
boli hl’adan´e len na intervale, ktor´
y odpovedal ±5% od hodnoty pr´ısluˇsn´eho parametru z´ıskanej v prvom kroku.
Tak´
ymto spˆosobom je moˇzn´e identifikovat’ aj vplyv poruchov´
ych veliˇc´ın, ktor´e nie je
moˇzn´e spr´avne identifikovat’ z d´at ktor´e s´
u zat’aˇzen´e korel´aciou vd’aka pr´ıtomnosti regul´atoru, ˇco uk´aˇzeme aj na simul´aciach uveden´
ych v nasleduj´
ucej ˇcasti.
6.3.1
Simul´
acie
V nasleduj´
ucej ˇcasti na simul´aciach uk´aˇzeme, ako v´
yrazne mˆoˇze zahrnutie vplyvu intenzity slneˇcn´eho ˇziarenia vylepˇsit’ predikˇcn´e spr´avanie identifikovan´eho modelu, a z´aroveˇ
n
uk´aˇzeme, ak´
y v´
yznam m´a v tomto pr´ıpade dvojkrokov´a identifikaˇcn´a proced´
ura, ktor´
u
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
67
sme pop´ısali v 6.3. Budeme porovn´avat’ tri rˆozne modely: model1 predstavuje model, ktor´
y
zah´rn
ˇa aj vplyv intenzity slneˇcn´eho ˇziarenia, ktor´
y bol identifikovan´
y klasicky v jednom
kroku z d´at zozbieran´
ych poˇcas vykurovacej sez´ony, model2 nezah´rn
ˇa vplyv intenzity
slnˇcn´eho ˇziarenia, model3 bol z´ıskan´
y identifikaˇcnou proced´
urou pop´ısanou v 6.3.
Ako identifikaˇcn´e d´ata sme v tom pr´ıpade pouˇzili d´ata zo zaˇciatku decembra 2010,
priˇcom sme ich vyberali tak, aby v obdob´ı, z ktor´eho poch´adzaju, bolo ˇco najviac slneˇcn´
ych dn´ı. Ako validaˇcn´e d´ata sme pouˇzili d´ata z janu´ara 2011, rovnako sme volili obdobie,
kedy bolo relat´ıvne slneˇcno. Pri identifik´acii sme pouˇzili met´odu, ktor´
u sme oznaˇcili v u
´vode tejto kapitoly ako MRI+GB, a parametre sme volili podobne ako pri identifik´acii
predoˇsl´
ych modelov: P = 100, na = 1, nb = 1, nk = 1.
model1 v sebe s´ıce zah´rn
ˇa vplyv slnka, ale kvˆoli tomu, ˇze bol identifikovan´
y z d´at
poˇcas vykurovacej sez´ony, nebolo moˇzn´e identifikovat’ vplyv slnka presne, a na validaˇcn´
ych
dat´ach u
´plne str´aca predikˇcn´
u schopnost’ na viac ako p´ar hod´ın dopredu. Jeho predikcie
s´
u najm¨a na severnej strane ovel’a nepresnejˇsie ako u model2 , ktor´
y vˆobec nezah´rn
ˇa vplyv
slnka.
Naopak model3 , ktor´
y bol z´ıskan´
y dvojkrokovou identifikaˇcnou proced´
urou, predikuje jednoznaˇcne najpresnejˇsie. Ako jedin´
y dok´aˇze predpovedat’ v´
yrazne zv´
yˇsenie teploty
okolo obeda, kedy teplota v miesnosti znaˇcne st´
upne, a rovnako i v´
yrazn´
y pokles teploty
poˇcas noci. Zauj´ımav´a je situ´acia, ktor´
u je vidiet’ na obr´azku 6.19, kde predikcie v obdob´ı od 16.1. do 18.1. s´
u pre model1 a model3 v´
yrazne vych´
ylen´e, no model2 predikuje
vel’mi presne, ˇco je spˆosoben´e t´
ym, ˇze pr´ave tieto dva dni boli za dan´e obdobie zd’aleka
najslneˇcnejˇsie.
juh
0.4
model
1
model2
model
3
varrs(%)
0.3
0.2
0.1
0
20
40
60
80
i
100
120
140
160
100
120
140
160
sever
0.4
model1
model
2
varrn
0.3
model3
0.2
0.1
0
20
40
60
80
i
Obr. 6.17: Porovnanie rozptylu predikˇcn´
ych ch´
yb identifikovan´
ych modelov.
´
´
KAPITOLA 6. SIMULACIE
NA REALNYCH
DATACH
68
juh
100
model
1
90
model2
80
model
3
fitrs(%)
70
60
50
40
30
20
20
40
60
80
i
100
120
140
160
sever
100
model1
90
model2
80
model
3
fitrn(%)
70
60
50
40
30
20
20
40
60
80
i
100
120
140
160
Obr. 6.18: Porovnanie fitfaktorov identifikovan´
ych modelov.
160 − kroková predikcia (juh)
160 − kroková predikcia(sever)
24.5
24.5
24
24
23.5
23.5
23
23
22.5
22.5
22
ϑ(°C)
ϑ(°C)
22
21.5
21.5
21
21
20.5
20.5
20
namerané data
model
20
model
2
19.5
3
19
model2
model3
model
16.1.0:00 18.1.0:00 20.1.0:00 22.1.0:00 24.1.0:00
namerané data
model
1
1
19.5
19
16.1.0:00 18.1.0:00 20.1.0:00 22.1.0:00 24.1.0:00
Obr. 6.19: Porovnanie predikci´ı identifikovan´
ych modelov.
Kapitola 7
Z´
aver
V r´amci tejto pr´ace bola spracovan´a ˇst´
udia identifikaˇcn´
ych met´od vhodn´
ych pre predikt´ıvne riadenie. Vel’k´a pozornost’ bola venovan´a identifkaˇcn´
ym net´odam, ktor´e s´
u zaloˇzen´e
na minimaliz´acii viackrokovej predikˇcnej chyby a poskytuj´
u model, ktor´
y je priamo urˇcen´
y
pre pouˇzitie s MPC. Naˇsa ˇst´
udia potvrdila d’aleko lepˇsie predikˇcn´e vlastnosti t´
ychto
modelov aj pri testovan´ı na vzorov´
ych pr´ıkladoch, ale predovˇsetk´
ym, ˇco je hlavn´e, modely
z´ıskan´e MRI identifikaˇcn´
ymi met´odami potvrdili svoj potenci´al aj pri identifik´acii modelu
ˇ
konkr´etnych blokov budovy CVUT
v Praze. Teplotu v miestnosti dok´azali predikovat’
omnoho presnejˇsie ako doteraz pouˇz´ıvan´e modely z´ıskan´e met´odami Subspace identifik´acie
a Grey box modelovac´ım pr´ıstupom.
Hlavne vd’aka t´
ymto vlastnostiam boli MRI identifikaˇcn´e met´ody vybrat´e z mnoˇzstva
dostupn´
ych identifikaˇcn´
ych met´od a prostredn´ıctvom nich boli vytvoren´e modely vˇsetk´
ych
ˇ
blokov budovy CVUT
v Dejviciach, ktor´e vyuˇz´ıva predikt´ıvny regul´ator riadiaci teplotu
v jednotliv´
ych miestnostiach pre predikcie.
V d’alˇsej ˇcasti naˇsej pr´ace bol vykonan´
y n´avrh predikt´ıvneho regul´atora tak, aby
zabezpeˇcil nielen teplotn´
y komfort vo vybran´
ych miestnostiach, ale aby ˇco najviach zn´ıˇzil
n´aklady na vykurovanie. Podl’a porovnania, ktor´e sme uviedli v ˇcasti 6.2, predikt´ıvny
regul´ator vykazoval oproti klasickej ekvitermnej regul´acii relat´ıvne u
´spory okolo 27%.
Aj napriek znaˇcn´emu zn´ıˇzeniu n´akladov na vykurovanie, ktor´e MPC prinieslo, je v tejto oblasti eˇste vel’k´
y priestor na vylepˇsovanie, predovˇsetk´
ym v oblasti identifik´acie. Znaˇcn´e
vylepˇsenie vlastnost´ı identifikovan´
ych modelov mˆoˇze priniest’ zahrnutie vplyvu intenzity
slneˇcn´eho ˇziarenia. Ako bolo uk´azan´e na simul´aciach v 6.3, modely rozˇs´ıren´e o vplyv
intenzity slneˇcn´eho ˇziarenia dok´aˇzu hlavne poˇcas slneˇcn´
ych dn´ı presnejˇsie predikovat’
teplotu v miestnostiach. Pouˇzitie tak´ehoto modelu v kombin´acii s MPC by mohlo zabr´anit’
miernemu pret´apaniu najm¨a poˇcas obeda, ked’ je intenzita slneˇcn´eho ˇziarenia najvyˇsˇsia,
69
´
KAPITOLA 7. ZAVER
70
a zn´ıˇzit’ t´
ym n´aklady na vykurovanie.
Priestor pre d’al’ˇs´ı v´
yskum predstavuj´
u aj samotn´e MRI identifikaˇcn´e met´ody. Probl´emom, s ktor´
ym sme sa pot´
ykali aj pri naˇsej pr´aci, bol nedostatok dobre vybuden´
ych
d´at. Tento jav je beˇzn´
y pri identifik´acii takmer vˇsetk´
ych priemyseln´
ych syst´emov, ked’ˇze
pracuj´
u v uzavretej sluˇcke vd’aka pr´ıtomnosti regul´atora. Zohl’adenie tohto faktu pri identifik´acii zaloˇzenej na minimaliz´acii viackrokovej predikˇcnej chyby by mohlo znaˇcne zv´
yˇsit’
pouˇzitel’nost’ t´
ychto met´od pri identifik´acii pre predikt´ıvne riadenie.
Literat´
ura
[1] T.N. Adlam. Radiant heating. The Industrial Press, 1949.
[2] A. Bemporad, M. Morari, V. Dua, and E.N. Pistikopoulos. The explicit linear quadratic regulator for constrained systems* 1. Automatica, 38(1):3–20, 2002.
[3] S.P. Boyd and L. Vandenberghe. Convex optimization. Cambridge Univ Pr, 2004.
[4] W.W. Chin. The partial least squares approach to structural equation modeling.
1998), S, pages 295–336, 1998.
[5] J. Cigler and S. Pr´ıvara. Subspace identification and model predictive control for
buildings. In The 11th International Conference on Control, Automation, Robotics
and Vision – ICARCV2010, pages 750–755.
[6] T. Coleman, M.A. Branch, and A. Grace. Optimization Toolbox For Use with MATLAB User’s Guide Version 2. 1999.
[7] Richard Godfrey Crittall and Joseph Leslie Musgrave. Heating and cooling of buildings. GB Patent No. 210880, April 1927.
[8] K. De Cock and B. De Moor. Subspace identification methods. Contribution to
section, 5:933–979, 2003.
[9] S. de Jong. SIMPLS: an alternative approach to partial least squares regression.
Chemometrics and Intelligent Laboratory Systems, 18(3):251–263, 1993.
[10] E.R.G. Eckert. Introduction to the Transfer of Heat and Mass. McGraw-Hill, 1950.
ˇ
[11] Luk´aˇs Ferkl and Jan Sirok´
y. Ceiling radiant cooling: Comparison of ARMAX and
subspace identification modelling methods. Building and Enviroment, 45(1, Sp. Iss.
SI):205–212, 2010.
71
´
LITERATURA
72
[12] C. Fornell and F.L. Bookstein. Two structural equation models: LISREL and PLS
applied to consumer exit-voice theory. Journal of Marketing research, 19(4):440–452,
1982.
[13] J.D. Franklin, G. Powell. Feedback Control of Dynamic Systems. Prentice Hall, 2009.
[14] P. Geladi and B.R. Kowalski. Partial least-squares regression: a tutorial. Analytica
Chimica Acta, 185:1–17, 1986.
[15] M. Gevers. A decade of progress in iterative process control design: from theory to
practice. Journal of process control, 12(4):519–531, 2002.
[16] R.B. Gopaluni, R.S. Patwardhan, and S.L. Shah. Bias distribution in MPC relevant
identification. In Proceedings of IFAC world congress, Barcelona, pages 2196–2201.
Citeseer, 2002.
[17] R.B. Gopaluni, R.S. Patwardhan, and S.L. Shah. MPC relevant identification–tuning
the noise model. Journal of Process Control, 14(6):699–714, 2004.
[18] M. Gwerder, B. Lehmann, J. Todtli, V. Dorer, and F. Renggli. Control of thermallyactivated building systems (TABS). Applied energy, 85(7):565–581, 2008.
[19] D. Gyalistras and M. Gwerder. Use of weather and occupancy forecasts for optimal
building climate control (opticontrol): Two years progress report. Technical report,
Terrestrial Systems Ecology ETH Zurich, Switzerland and Building Technologies
Division, Siemens Switzerland Ltd., Zug, Switzerland, 2010.
ˇ
ˇ
[20] V. Havlena and J. Stecha.
Modern´ı teorie ˇr´ızen´ı. Praha: Vydavatelstv´ı CVUT,
2000.
[21] I. Jolliffe. Principal component analysis. 2002.
[22] T. Katayama. Subspace Methods for System Identification. Springer, 2005.
[23] E.C. Kerrigan and J.M. Maciejowski. Soft constraints and exact penalty functions
in model predictive control. In Control 2000 Conference, Cambridge, 2000.
[24] W.H. Kwon, S. Han, and S. Han. Receding horizon control: model predictive control
for state models. Springer Verlag, 2005.
[25] D. Laur´ı, M. Mart´ınez, J.V. Salcedo, and J. Sanchis. PLS-based model predictive
control relevant identification: PLS-PH algorithm. Chemometrics and Intelligent
Laboratory Systems, 100(2):118–126, 2010.
´
LITERATURA
73
[26] J.W. Lee, W. HYUN KWON, and J. Choi. On stability of constrained receding
horizon control with finite terminal weighting matrix. Automatica, 34(12):1607–
1612, 1998.
[27] F. Lindgren, P. Geladi, and S. Wold. The kernel algorithm for PLS. Journal of
Chemometrics, 7(1):45–59, 1993.
[28] L. Ljung. System identification toolbox. The MathWorks Inc, 1999.
[29] Lennart Ljung. System Identification: Theory for user. Prentice-Hall, 1999.
[30] J.M. Maciejowski. Predictive Control with Constraints. Prentice-Hall, 2002.
[31] D.Q. Mayne, J.B. Rawlings, C.V. Rao, and PO Scokaert. Constrained model predictive control: Stability and optimality. AUTOMATICA-OXFORD-, 36:789–814,
2000.
[32] J. Nocedal and S.J. Wright. Numerical optimization. Springer verlag, 1999.
[33] L. Perez-Lombard, J. Ortiz, and C. Pout. A review on buildings energy consumption
information. Energy and Buildings, 40(3):394–398, 2008.
[34] E.N. Pistikopoulos. Multi-parametric model-based control: theory and applications.
Vch Verlagsgesellschaft Mbh, 2007.
ˇ
[35] S. Pr´ıvara, J. Cigler, Z. V´an
ˇa, L. Ferkl, and M. Sebek.
Subspace identification of
poorly excited industrial systems. In Proceedings of the 49th IEEE Conference on
Decision and Control, pages 4405–4410, 2010.
ˇ
[36] S. Pr´ıvara, J. Sirok´
y, L. Ferkl, and J. Cigler. Model predictive control of a building
heating system: The first experience. Energy and Buildings, 2010.
[37] J. Rawlings and K. Muske. Stability of constrained receding horizon. IEEE Transaction on Automatic Control, 38, 1993.
[38] J.A. Rossiter and B. Kouvaritakis. Modelling and implicit modelling for predictive
control. International Journal of Control, 74(11):1085–1095, 2001.
[39] D. Schmidt. Modelling of hybrid building components with r-c networks in macro
elements.
´
LITERATURA
74
[40] D.S. Shook, C. Mohtadi, and S.L. Shah. A control-relevant identification strategy
for GPC. Automatic Control, IEEE Transactions on, 37(7):975–980, 2002.
[41] D.S. Shook, C. Mohtadi, and S.L. Shah. Identification for long-range predictive
control. In Control Theory and Applications, IEE Proceedings D, volume 138, pages
75–84. IET, 2002.
[42] T. S¨oderstr¨om and P. Stoica. System identification. 1989.
[43] Pavel Trnka. Subspace Identification Methods. PhD thesis, Czech Technical University in Prague, 2007.
[44] Peter Van Overschee and Bart De Moor. Subspace Identification for Linear Systems.
Kluwer Academic Publishers, 1999.
[45] S. Wold, M. Sjostrom, and L. Eriksson. PLS-regression: a basic tool of chemometrics.
Chemometrics and intelligent laboratory systems, 58(2):109–130, 2001.
[46] S. Wold, J. Trygg, A. Berglund, and H. Antti. Some recent developments in PLS
modeling. Chemometrics and Intelligent Laboratory Systems, 58(2):131–150, 2001.
[47] R.K. Wood and M.W. Berry. Terminal composition control of a binary distillation
column. Chemical Engineering Science, 28(9):1707–1717, 1973.
[48] Y. Zhu. Multivariable system identification for process control. Elsevier, 2001.
ˇ
ˇ
[49] J. Stecha.
Optim´
alni rozhodovan´ı a ˇr´ızen´ı. Praha: Vydavatelstv´ı CVUT,
1999.
Pr´ıloha A
Obsah priloˇ
zen´
eho CD
S´
uˇcast’ou tejto pr´ace je aj CD, ktor´e obsahuje okrem elektronickej verzie tejto pr´ace aj
naimplementovan´e zdrojov´e k´ody. CD m´a nasleduj´
ucu ˇstrukt´
uru:
• DP obsahuje diplomova praca.pdf,
• Scilab obsahuje zdrojov´e k´ody naimplementovan´e v programovom prostred´ı Scilab,
• Matlab obsahuje zdrojov´e k´ody naimplementovan´e v programovom prostred´ı Matlab.
I
Download

DIPLOMOV´A PR´ACA