VLIV NEJISTOT MĚŘENÍ NA VELIKOST INTERVALU SPOLEHLIVOSTI PROVOZNÍHO OPTIMA Javůrek M., Taufer I. Univerzita Pardubice, Fakulta elektrotechniky a informatiky, Katedra řízení procesů, e-mail: milan.javurek,
[email protected] Abstrakt: V příspěvku je prezentován způsob simulovaného stanovení optimálního režimu reaktoru, ve kterém probíhá následná reakce A → B → C. Je popsán algoritmus optimalizační metody zpětného kroku a uvedeny výsledky experimentálního stanovení optima. Jsou popsány způsoby stanovení mezních hodnot intervalů spolehlivosti a z experimentálních dat vyhodnoceny intervaly spolehlivosti v závislosti na různých hodnotách nejistot „měřených“ hodnot. Nejistota v datech způsobuje jednak zvětšování šíře intervalu spolehlivosti, jednak větší migraci odečítaných hodnot, což způsobuje zvýšení charakteristik těsnosti proložení regresní závislosti. Je konstatována dobrá shoda experimentálních výsledků a teoretických předpokladů. Klíčová slova: technologický proces, optimalizace, nejistoty měření, interval spolehlivosti Úvod Základním požadavkem na vedení každého technologického procesu je podmínka jeho efektivního provozu, což spočívá v zajištění minimálních nebo maximálních hodnot vybraných parametrů y, např. minimum provozních nákladů, maximum zisku, maximum výtěžnosti apod. Takový režim, který vykazuje extremní hodnoty těchto parametrů, označujeme jako režim optimální a řídicí parametry u, který tento režim zajišťují, parametry optimální. Při tom však je třeba vzít v úvahu, že technologický proces je proces dynamický, kdy se po změně řídicího parametru příslušná veličina ustálí až po uplynutí určité doby. Dále mohou na proces působit vnější vlivy, resp. chyby měření, souhrnně označované jako poruchy v, které nastavení nebo vyhodnocení optimálního režimu zkreslují. Blokové schéma takového procesu je uvedeno na obr. 1. v y
u TECHNOLOGICKÝ PROCES
Obr. 1 – Schéma technologického procesu Optimalizace Pod pojmem optimalizace resp. optimální řízení technologického procesu budeme chápat cílevědomou činnost člověka (obsluhy, technologa), která v daných provozních podmínkách (omezujících parametrech) zajišťuje co nejlepší výsledky provozu tohoto procesu [Drábek, 1990]. Úloha optimalizace pak spočívá v tom, abychom při daných technologických parametrech (např. výkon zařízení, kvalita suroviny a pomocných materiálů aj.), poruchách v (změna okolního prostředí, změna kvality suroviny, chybách měření apod.) našli takové 1
hodnoty řídicích veličin u (množství suroviny, teplota, tlak atd.), při nichž výstupní řízené veličiny y (kvalitativní ukazatele výrobku, výkonnost zařízení, spotřeba energie apod.) splňují požadovaná hlediska optimálnosti. Pro možnost posouzení tohoto hlediska je proto nutné definovat jednoznačné kritérium, kritérium optimalizace, které toto optimum kvantifikuje. Vzhledem k tomu, že v reálném technologickém procesu nemáme k dispozici matematický popis, který by charakterizoval vztah mezi kritériem optimalizace a jeho řídicími veličinami, je třeba k nastavení a vyhodnocení optimálního režimu přistoupit experimentálně. K tomu slouží celá řada metod a algoritmů. Algoritmy experimentální optimalizace Tyto algoritmy spadají do oblasti tzv. nelineárního programování [Drábek, 1990, Hrubina, 2001]. Jde o numerické metody iteračního charakteru, kdy se nové hodnoty nezávisle proměnných uk v daném iteračním kroku k počítají z hodnot kroku předcházejícího připočtením nějakého přírůstku Δuk. Takže platí: (1)
u k = u k −1 + Δu k
Je zřejmé, že pokud se hledá maximum dané funkce y(u), potom úspěšný krok bude charakterizován nerovností (2)
y (uk ) > y (uk −1 )
V opačném případě je přechod do stavu uk nežádoucí a pro výpočet další hodnoty uk je třeba určit novou hodnotu Δuk a výpočet opakovat. A právě výpočet této hodnoty Δuk charakterizuje příslušnou metodu řešení. Výpočet se ukončí, jestliže chyba řešení, vyjádřená jako norma rozdílu dvou po sobě jdoucích vektorů řídicích proměnných uk a uk-1, klesne pod určenou hodnotu. Pak platí: (3)
uk − uk −1 < ε
Technologický proces Průběh optimalizace technologického procesu bude demonstrován na modelu průtočného, dokonale míchaného chemického reaktoru, ve kterém probíhá následná chemická reakce A → B → C [Taufer, 1993]. Tuto reakci lze popsat matematickým vztahem xB =
kde
τ k1
(1 + τ k1 )(1 + τ k 2 )
(4)
x A0
k1 , k 2 jsou reakční rychlosti podle Arheniova vztahu: −
E
k = Ae RT
(5)
V Q
(6)
τ je doba zádrže reakční směsi v reaktoru τ=
x A0 - vstupní koncentrace látky A
A - frekvenční faktor E - aktivační energie R - plynová konstanta T - reakční teplota V - objem reaktoru Q - průtok směsi reaktorem 2
Technologické schéma tohoto průtočného, dokonale míchaného chemického reaktoru je uvedeno na obr. 2. Do reaktoru vstupuje reakční směs o průtoku Q a dané koncentraci xA0. V reaktoru je udržován konstantní objem V, vyjádřený výškou hladiny h, regulací výtoku reakční směsi z reaktoru. Požadovaná teplota T reakční směsi v reaktoru se reguluje průtokem chladícího média v chladícím potrubí. Cílem optimálního řízení je zajistit změnou teploty v reaktoru maximální výtěžnost reakce, která je vyjádřená koncentrací složky B.
Obr. 2 – Technologické schéma reaktoru Výše uvedená rovnice (4) představuje jednoparametrickou nelineární závislost koncentrace komponenty B na teplotě reakční směsi v reaktoru. Pro přiblížení se k realitě jsou všechny parametry rovnice (4) (vstupní koncentrace složky A, průtok reakční směsi reaktorem, teplota reakční směsi a výstupní koncentrace složky B) zatíženy náhodnou chybou s normálním rozložením, čímž jsou vyjádřeny nestabilní hodnoty (poruchy) těchto parametrů a nejistoty jejich měření. Dále je simulována i dynamika procesu zavedením přenosu 1. řádu. Blokové schéma tohoto procesu je uvedeno na obr. 3. vT
vA0
vQ
vB xB
T CHEMICKÝ REAKTOR
Obr. 3 – Blokové schéma reaktoru Pro demonstraci procesu optimalizace jednoparametrického hledání - metodu zpětného kroku.
3
použijeme
známou
metodu
Metoda zpětného kroku Základní princip metody [Drábek, 1990, Hrubina, 2001, Taufer, 2009] spočívá v zajištění změny řídicí veličiny u s krokem Δu ve směru stoupající hodnoty kritéria optimalizace (účelové funkce) y(u) podle vztahu (7)
u k = u k −1 + Δu k
Výchozí krok Δuk se volí jako určitý díl intervalu, ve kterém se má stanovit maximum funkce y(u). Výchozí krok můžeme zvolit např. podle vztahu Δu1 = h (u max − u min )
kde
(8)
h < 1. Tento krok je konstantní až do okamžiku, kdy přestane platit nerovnost (9)
y (u k ) > y (u k −1 )
Poté se hodnota kroku zmenší a obrátí se směr postupu hledání. Nová hodnota kroku se určí ze vztahu Δu k = −
kde
Δu k −1 z
(10)
z je dělicí faktor kroku. Výpočet se ukončí, když absolutní hodnota kroku klesne pod přípustnou chybu řešení (11)
Δu k ≤ ε
Obr. 4 – Průběh hledání optima
Experiment a) b) c) d)
Pro realizaci experimentu byly zvoleny: rozsah experimentu: umin = Tmin = 100 °C; umax = Tmax = 180 °C; počáteční dělicí poměr h = 1/3; průběžný dělicí faktor z = 3; chyba řešení ε = 0,5 °C.
Výsledky experimentu – hledání optima metodou zpětného kroku jsou uvedeny v tab. 1., kde Tž je podle výše uvedeného algoritmu vypočítaná žádaná hodnota teploty, T je teplota „měřená“ a xB rovněž „měřená“ koncentrace složky B. Proces optimalizace byl ukončen na 19. kroku experimentu, kdy absolutní hodnota rozdílu dvou následných hodnot 4
teploty klesla pod hodnotu chyby řešení ε. Získané optimální hodnoty jsou Topt = 132,7 °C a odpovídající hodnota koncentrace složky B xB = 83,5 %. Grafický průběh hledání optima je uveden na obr. 4., kde jsou vyneseny postupně „měřené“ hodnoty koncentrace xB složky B v závislosti na „měřených“ hodnotách teploty T v reaktoru. Tab. 1 – Výsledky experimentu k
Tž, °C
T, °C
xB, %
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
100,00 180,00 153,33 126,67 100,00 108,89 117,78 126,67 135,56 144,44 141,48 138,52 135,56 132,59 129,63 130,62 131,60 132,59 132,26
99,6 180,4 155,4 124,6 98,3 110,9 117,1 125,6 135,4 145,1 141,4 138,8 137,2 134,0 128,6 133,0 130,9 131,3 132,7
71,6 63,5 79,4 82,9 70,8 77,3 80,7 82,5 83,6 82,0 82,5 83,1 83,2 83,6 82,6 83,4 83,4 82,9 83,5
Vyhodnocení experimentu Předpokládané nejistoty výsledků měření Jak bylo výše uvedeno, na výstupní sledovaný parametr technologického procesu působí řada vnějších vlivů daných nestabilitou tohoto procesu (náhodné změny parametrů surovin a energií) a chybami měření. Je proto třeba vymezit interval pravděpodobných hodnot, které by „změřená“ veličina, koncentrace látky B xB a teplota reakční směsi T mohla vlivem poruch vstupních veličin získat. Tuto standardní nejistotu, vyjádřenou jako směrodatnou odchylku [Chudý, 1999], lze odhadnout použitím metody dvoubodové aproximace [Meloun, 2004], pro kterou platí vztah (po úpravě) s( x B ) =
kde
m
∑
{G [ zi,opt + s( zi )] + G [ zi,opt − s( zi )]}2
i =1
G [ z i ,opt ± s ( z i )]
4m
(12)
je funkční závislost (4), jejíž hodnoty jsou počítány postupně pro
jednotlivé mezní hodnoty sledovaných vstupních veličin, které jsou dány připočtením nebo odečtením jejich standardních nejistot k optimální hodnotě m – počet sledovaných poruchových veličin Dále lze vypočítat rozšířenou nejistotu jako interval spolehlivosti pro danou hladinu významnosti [Chudý, 1999], která bude vyjadřovat celkový možný rozsah pravděpodobných hodnot výstupní koncentrace xB
5
(13)
s e ( x B ) = t1−α ,ν s ( x B )
kde
t1−α ,ν je kritická hodnota Studentova rozdělení na hladině významnosti α a se stupněm
volnosti ν. V daném modelu reaktoru byly tyto poruchy simulovány jako náhodné veličiny s normálním rozložením a definovaným rozsahem, působící na vybrané parametry. Byly simulovány poruchy s rozsahem nejistot: teploty směsi v reaktoru sT = 0,1 – 0,25 – 0,5 – 1,0 – 2,5 °C, průtoku reakční směsi sQ = 0,04 – 0,1 – 0,2 – 0,4 – 1,0 m3 h-1, vstupní koncentrace sxA0 a výstupní koncentrace sxB = 0,2 – 0,5 – 1,0 – 2,0 – 5,0 %. Pro tyto rozsahy nejistot byly pak vypočítány: standardní nejistota s(xB) , rozšířená nejistota na 5% hladině významnosti se(xB) a interval spolehlivosti xBmin a xBmax . Tomu pak podle vztahu (4) odpovídá interval spolehlivosti pro teplotu Tmin a Tmax . Výsledky pro jednotlivé volby nejistot vstupních veličin jsou uvedeny v tab. 2. Tab. 2 – Vypočtené charakteristiky provozního optima pro volené nejistoty měřených veličin Zadané nejistoty měřených veličin Rozměr Hodnota 3 -1 m .h sQ 0,04 0,10 0,20 °C sT 0,10 0,25 0,50 % sxB, A0 0,20 0,50 1,00 Vypočtené charakteristiky Veličina Rozměr Hodnota % s(xB) 0,1312 0,3256 0,6515 s(xA0) pro 5% hl. v. % 0,2797 0,6940 1,3887 % xBmin 82,98 82,59 81,87 % xBmax 83,54 83,98 84,65 °C Tmin 128,80 126,00 122,50 °C Tmax 139,75 142,60 146,17 °C sr(T) 0,0442 0,0676 0,0647 sr(T) pro 5% hl.v. °C 0,0936 0,1433 0,1371 Koefic. polynomu a2 -0,00976 -0,0098 -0,0096 Koefic. polynomu a1 2,64072 2,6521 2,6169 Koefic. polynomu a0 -95,31179 -96,2875 -94,2413 Topt (proložení) °C 135,3 135,40 135,60 Koef. korelace r 0,99948 0,9990 0,9982 Rezid.rozptyl sr(xB) % 0,1925 0,2646 0,3542 Int.sp. se(xB) % 0,4081 0,5610 0,7510 Veličina
0,40 1,00 2,00
1,00 1,00 2,50
1,3031 2,7774 80,48 86,03 117,70 151,10 0,0821 0,1740 -0,0101 2,7282 -101,3915 135,60 0,9952 0,5956 1,2626
3,2575 6,9432 76,29 90,17 108,35 161,04 0,2324 0,4926 -0,0111 2,9784 -116,5260 134,50 0,9876 1,0441 2,2133
Statistické vyhodnocení experimentálních výsledků Výše získané výsledky představují maximální možné rozsahy vycházející z mezních hodnot poruchových veličin. Pro zhodnocení reálných výsledků „měření“ a jejich porovnání s výše vypočítanými mezními hodnotami je proto třeba následně experimentální data vyhodnotit statisticky. K tomu použijeme klasickou regresní analýzu, pomocí které lze mj. vypočítat reziduální rozptyl jakožto míru rozptylu experimentálních hodnot kolem regresní funkce yˆ = f ( x ) podle obecného vztahu 6
N
∑ ( y i − yˆ i )2 sr =
kde
i =1
(14)
N −1
jsou experimentální hodnoty yˆ i odpovídající hodnoty vypočítané pro příslušné nezávisle proměnné x N – počet experimentů yi
Pro získání standardní a následně rozšířené nejistoty teploty reakční směsi budeme vycházet z předpokladu lineární závislosti mezi experimentálně zjištěnou a žádanou teplotou. Pak bude platit N
∑ (T ž ,i − Ti )2 s r (T ) =
i =1
(15)
N −1
Byla takto vypočítána standardní nejistota sr(T) a odpovídající rozšířená nejistota na 5% hladině významnosti se(T) , kterou lze užít pro definici intervalu spolehlivosti. Tyto veličiny charakterizují subjektivní faktor - míru shody naměřené veličiny s její předvolenou hodnotou, tj. „pečlivost“ obsluhy při odečítání hodnot měřených veličin. Výsledky jsou uvedeny v tab. 2. Závislost koncentrace výstupního produktu xB na teplotě T byla aproximovaná regresním polynomem druhého stupně. Byl získán vztah
x B = a 0 + a1 T − a 2 T 2
(16)
s koeficientem korelace r, který představuje shodu mezi experimentem a regresní závislostí. Odpovídající reziduální rozptyl vyjadřuje standardní nejistota sr(xB) a interval spolehlivosti na 5% hladině významnosti se(xB). Výsledné hodnoty jsou opět uvedeny v tab. 2. Všechny získané výsledky jsou souhrnně graficky uvedeny na obr. 5 - 9.
Obr. 5 – Interval spolehlivosti pro zvolené nejistoty sQ=0,04, sT=0,1 a sxB=0,2
7
Obr. 6 – Interval spolehlivosti pro zvolené nejistoty sQ=0,1, sT=0,25 a sxB=0,5
Obr. 7 – Interval spolehlivosti pro zvolené nejistoty sQ=0,2, sT=0,5 a sxB=1
Obr. 8 – Interval spolehlivosti pro zvolené nejistoty sQ=0,4, sT=1 a sxB=2
8
Obr. 9 – Interval spolehlivosti pro zvolené nejistoty sQ=1, sT=2,5 a sxB=5
Plným kroužkem a pořadovým číslem experimentu jsou vyznačeny „měřené“ body koncentrace sledované látky B, kroužkem s pořadovým číslem 19 je označena „změřená“ optimální hodnota koncentrace xB. Čerchované vodorovné čáry vyznačují rozšířenou nejistotu mezních hodnot výstupní koncentrace. Rozšířené mezní nejistoty pro teploty jsou mimo rozsah grafu a proto nejsou zobrazeny. Plnou čarou je zobrazena část regresní závislosti xˆ B = f (T ) v okolí optima a na ní čtvercovou značkou její optimální hodnota. Tečkované čáry vyznačují nejistotu experimentálních hodnot koncentrace, svislé dvojitě čerchované čáry pak nejistotu experimentálních hodnot teploty. K řešení výše uvedených úloh byl využit tabulkový procesor EXCEL a jeho nástroje Řešitel, Hledání řešení a Přidání spojnic trendů. Závěr Při hledání optimálních podmínek procesu v reaktoru není možné se omezit pouze na vyjádření sledovaných veličin jedinou hodnotou. Na sledované veličiny působí jednak poruchy (vnější vlivy i nestabilita surovin), jednak je měření zatíženo nejistotami měření, jejich ustalování je dynamické a odečítání hodnot měřených veličin je zatíženo i subjektivním faktorem. Tyto vlivy je třeba vzít v úvahu a optimum veličin stanovit formou intervalu, ve kterém lze režim reaktoru ještě považovat za optimální. Přijatelný výsledek bude jen v tom případě, když hodnoty veličin zajišťující optimální provozní režim a jejich nejistoty budou ležet uvnitř vymezené oblasti mezních hodnot nejistot. Pro určení tohoto intervalu byly v našem případě použity metoda dvoubodové aproximace a metoda určení reziduálního rozptylu okolo regresní křivky proložené naměřenými body. Na základě získaných výsledků lze konstatovat, že takto experimentálně stanovený optimální režim je v rozsahu předpokládaných intervalů spolehlivosti a experimentální výsledky nevybočují z mezí zjištěných nejistot měření. Bylo provedeno testování vlivu velikosti nejistot měřených veličin na šíři intervalu spolehlivosti, který se při zvyšující se nepřesnosti měření zvětšuje. Je možno podle požadované velikosti intervalu definovat potřebnou minimální přesnost měření jednotlivých vstupních veličin. Se zvyšující nejistotou se kromě šíře intervalu spolehlivosti zhoršují charakteristiky těsnosti proložení regresní závislosti.
9
Problematika je řešena v rámci výzkumného záměru MŠM 0021627505 „Řízení, optimalizace a diagnostika složitých systémů“. Abstract: The contribution presents a method of simulated determination of the optimum regime of a reactor in which a competitive consecutive system of reactions takes place: A → B → C. The algorithm has been described for the optimisation method of reverse step, and the results of experimental determination of the optimum are given. Also described are the ways of determination of limiting values of reliability intervals, and from experimental data the reliability intervals have been evaluated depending on various values of uncertainty of “measured” quantities. The data uncertainty broadens the reliability interval and also increases the migration of reading values, which causes an increase in the characteristic of goodness of fit of the regression dependence. The experimental results and theoretical presumptions have been found to agree well. Key words: Technological process, optimization, uncertainties of measurement, confidence interval. Literatura DRÁBEK, O.; TAUFER, I. Automatizované systémy řízení technologických procesů. [Skriptum]. Pardubice : VŠCHT Pardubice, 1990. 90 s. ISBN 80-85113-16-3. HRUBINA, K.; et al. Mathematical modelling of technical processes. Prešov (Slovak Republic) : Fakulty of manufacturing Technologies on the UT Košice, 2001. 216 p. ISBN 80-88941-18-0. CHUDÝ, V.; PALENČÁR, R.; KUREKOVÁ, E.; HALAJ, M. Meranie technických veličín. Bratislava : STU Bratislava, 1999. 688 s. ISBN 80-227-1275-2. Dostupné též na
. KUBANOVÁ, J. Statistické metody pro ekonomickou a technickou praxi. Bratislava : STATIS, 2004. 249 s. ISBN 80-85659-37-9. MELOUN, M.; MILITKÝ, J. Statistická analýza experimentálních dat. Praha: Academia, 2004, 953 s., ISBN 80-200-1254-0. OLEHLA, M.; aj. Řešení úloh matematické statistiky ve FORTRANU. Praha : NADAS, 1982. 368 s. OD-31-035-82 – 01-13. SCALES, L.E. Introduction to Non-linear Optimization. Sutton. Surey : MacMillan Publishers Ltd., 1985, 243 s., ISBN 0-333-32553-2. TAUFER, I.; KOTYK, J.; HRUBINA, K.; TAUFER, J. Algoritmy a algoritmizace. Vývojové diagramy. Sbírka řešených příkladů. [Skriptum]. Pardubice : Univerzita Pardubice, 2009. 92 s. ISBN 978-80-7395-182-5. TAUFER, I.; KREJČÍ, S.; JAVŮREK, M. Laboratorní cvičení OPTIMALIZACE. Metody nelineárního programování. [Skriptum], Pardubice: VŠCHT, 1993.
10
