NUMERICKÁ ANALÝZA PROCESů

1. NUMERICKÁ ANALÝZA PROCESů NAP2 Fyzikální modely (transportní rovnice, energetické a entropické principy). Empirické modely (neuronové sítě a regresní modely). Modely, které využívají analytická řešení (difuze). Identifikace modelů na základě co nejlepší shody s experimentem (optimalizace). Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. Identifikace parametrů modelu na základě experimentu Matematický model aparátu nebo systému Vyhodnocení experimentů – diagnostika aparátů, např. diagnostika kolon, EEG Optimalizace operačních parametrů (optimalizace návrhu) MODELY PROCESů NAP2 Matematický model= matematický popis (program) charakteristik systému (rozložení teplot, rychlostí, tlaků, koncentrací, výkonu,...), jako funkce času t, místa a operačních parametrů (geometrie systému, materiálových vlastností, průtoků,...). Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

3. MODELY PROCESů NAP2 Obecná klasifikace matematických modelů • Black box(černá krabička) – model založen výhradně na porovnání s experimentem. Tři základní typy jsou • Neuronové sítě (analogie neuronů v mozku – metoda umělé inteligence) • Regresní modely (korelace např. ve tvaru mocninových funkcí typu Nu=cRenPrm ) • Identifikace přenosu systému (t), ze změřeného časového průběhu vstupu x(t) a výstupu y(t) systému. • White box(bílá krabička) – model je založen výhradně na ověřených fyzikálních principech: • na bilancích hmoty, hybnosti, energie(Fourierova a Fickova rovnice, Navier Stokesovy rovnice) • na energetických principech, hledané řešení minimalizuje celkovou energii systému (pružnost), • na simulacích náhodných pohybů a interakcí fiktivních částic (metody Monte Carlo, Boltzmann lattice…). • - Grey box(šedá krabička) - modely, které se nachází mezi výše uvedenými extrémy. Příkladem jsou kompartment modely průtočných systémů, respektující fyzikální princip zachování hmoty, ale další principy (třeba bilance impulsu) jsou u nich nahrazeny empirickými korelacemi nebo daty. Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

4. Fyzikální principy NAP2 Modely kontinua (bílé krabičky) Inženýrské modely založené na fyzikálních principech často vycházejí ze tří bilancí (hmoty, hybnosti a energie). V diferenciálním tvaru: Zachování hmoty Zachování hybnosti Zachování energie Cauchyho rovnice rovnováhy. Platí pro pružnost stejně jako pro proudění, stlačitelné i nestlačitelné. Celková energie, součet vnitřní, kinetické a potenciální energie výkon vnitřních sil hustota tepelného toku Většina rovnic, které znáte, se z těchto rovnic dá odvodit: Bernoulliho, Eulerova, Fourierova Kirchhoffova Navier Stokesova.

5. Wi W(u) u We Fyzikální principy NAP2 Modely kontinua (bílé krabičky) Zvláště v mechanice tuhé fáze se dává přednost principům založených na energii. V nauce o pružnosti je to Lagrangeův princip minima celkové potenciální energie deformace i napětí lze vyjádřit jako funkce posunutí u deformační energie vnitřních sil (součin tenzoru deformace a napětí) práce vnějších sil (objemové, povrchové, osamělé) ve stavu rovnováhy má celková potenciální energie (Wi+We) minimum Požadavek nulové variace funkcionálu energie W/u=0 je základem konstrukce mnoha konečných prvků. V dynamice je tento přístup reprezentován Hamiltonovým principem požadavku nulové variace funkcionálu, do něhož je zahrnuta i kinetická energie „Energetické“ principy nejsou tak obecné: z každého variačního principu lze odvodit jemu ekvivalentní Eulerovy Lagrangeovy diferenciální rovnice, ale obráceně to neplatí. Především u proudění lze aplikovat variační principy jen v některých případech (třeba plouživého proudění).

6. Fyzikální principy NAP2 Typy modelů souvisí s typem problémů a s typem inženýrů, kteří se řešením těchto problémů zabývají: FUNKCIONALISTÉ Představují si, že každý systém může nabývat nekonečně mnoha podob, a každé podobě přiřadí nějaké číslo, funkcionál. Pak si tato čísla vynesou do grafu a hledají místa, která jsou něčím zvláštní, obvykle minima, maxima, inflexe. A těmto místům přiřadí zvláštní mystický význam, třeba to, že reprezentují stav rovnováhy, stav ztráty stability apod. Pohybují se zejména v oblasti pružných těles, solidů a jsou tudíž obvykle SOLIDNÍ. Snadno je poznáte podle toho, že za jedinou použitelnou numerickou metodu považují konečné prvky (se skřípěním zubů snad i Boundary elements) a používají slůvka jako Cauchy Green tenzor, Kirchhoff - Piola stresses, sedmý invariant apod. DIFERENCIÁLOVÉ Věří zákonu zachování hmoty, hybnosti a energie a domnívají se, že z těchto bilancí dokáží popsat stav systému v nekonečně blízkém okolí libovolného bodu diferenciálními rovnicemi. Z řešení těchto rovnic pak věští realitu celku. Pohybují se zejména v oblasti transportních jevů a proudění, jsou to PROUDNÍCi. Při výběru numerické metody nejsou příliš vybíraví, každá metoda, která vyřeší jejich rovnice je jim dost dobrá, ale většinou volí metodu konečných objemů. Často používají slůvka jako materiálová derivace. PARTIKULARISTÉ Nevěří v mechaniku kontinua. Vše lze odvodit z mechaniky hmotného bodu. Existují pouze diskrétní veličiny jako v číslicových počítačích. Složitost reality je způsobena velkým množstvím paralelně běžících jednoduchých procesů. Poznáte je podle toho, že se straní všech ostatních a jsou DISKRÉTNÍ. FATALISTÉ Nevěří, že je v lidských silách poznat zákony božské prozřetelnosti a proto se omezují na popis pozorovaných jevů, jsou to EMPIRICI. Milují expertní systémy, metody umělé inteligence, inženýrské korelace a statistiku, i když se pragmaticky nevyhýbají ani používání šedých krabiček. CHAOTI Připouští existenci nějakých principů, ale nevěří v jejich smysluplnou řešitelnost. Nesnáší hladké křivky a raději je kazí třeba generátory náhodných čísel. Milují katastrofy, podivnosti a atraktory, jsou ATRAKTIVNÍ. Zvláště pro ženy, které okouzlují cizími slovy a barevnými květinami fraktálů.

7. Příklady modelů-sušení NAP2 Téměř vždy lze jako matematický model použít kterýkoliv z předchozích typů. Jako příklad (který nás v různých obměnách bude provádět celým kurzem) uvedeme sušení nebo zvlhčování materiálu, třeba zrna (škrobového, obilného, rýžového, kávového…). Vždy bude jedním ze základních výsledků modelu doba sušení nebo zvlhčování (resp. závislosti měrné vlhkosti na čase a na teplotě). Pro tento cíl většinou stačí jednoduché regresní modely nebo neuronové sítě. Pro stanovení mikrobiální aktivity, zdravotní nezávadnosti je třeba znát i distribuci vlhkosti uvnitř zrna. Pak nastupují modely bílé nebo šedé krabičky, modely založené na transportních rovnicích difúze a tepla. Když je geometrie zrna jednoduchá (kulička), a když není příliš významný vliv vlhkosti na difúzní součinitel, dá se řešení vyjádřit ve tvaru nekonečné řady (analytické řešení Fickovy a FK rovnice) nebo se spokojíme třeba jen s prvním (dominujícím) členem řady. Když je jednoduchá geometrie, ale silně nelineární koeficienty difúze a přenosu tepla, použije se zpravidla některá 1D numerická metoda (metoda sítí nebo spektrální metoda). U tvarově komplikovaných zrn (rýže, kávová zrna) se používají 3D numerické metody konečných prvků nebo kontrolních objemů (zpravidla se využijí komerční programy ANSYS, COMSOL nebo FLUENT). V současnosti je předmětem zájmu i měnící se vnitřní struktura materiálu, např. vznik trhlin, které zpětně transport vlhkosti výrazně ovlivňují (volná voda snáze proniká mikrotrhlinami zrna). To ovšem znamená, že kromě transportních rovnic se musí řešit i pole deformací a rozložení tenzoru mechanických napětí (pružnostní problém). Je to o to složitější, že u biologických materiálů jde o velké deformace (nelze tedy použít teorii malých deformací a např.Hookův zákon jako u běžných ocelových konstrukcí). Modely jsou zpravidla založené na metodě konečných prvků. Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

8. Neuronové sítě NAP2 Umělé neuronové sítě jsou technikou brutální umělé inteligence. Její obliba vyplývá z toho, že čím dál tím větší počet lidí ví o reálných procesech a jejich podstatě čím dál méně (snad by se to dalo nazvat zákonem exponenciálního růstu ignorance). ANN (Artificial Neural Network) jsou určeny právě pro ty, kteří o procesech, které chtějí modelovat, neví vůbec nic a ani nechtějí o jejich principech příliš přemýšlet. Mají jen spoustu experimentálních dat v nichž nejsou schopni se orientovat. A mají k dispozici MATLAB (například). Schiele

9. Řízení průtoku substrátu 3 1 4 7 Sensory 5 2 6 NEURONOVÁ SÍŤ - MOZEK Neuron má několik vstupů - dendrity a jen jeden výstup - axon Neuronové sítě NAP2 Neuron=modul, který vypočte jednu výstupní hodnotu z několika vstupních. Chování konkrétní neuronové sítě je dáno hodnotami synaptických váhových koeficientů wi (koeficientů zesílení signálů mezi propojenými neurony) Skrytá vnitřní vrstva (zde 4 neurony) Vstupní vrstva (zde 2 neurony) X1,X2 Výstupní vrstva (zde 1 neuron) Y Odezva neuronu y na N-vstupních hodnot xi Wi=koeficienty synaptických vah stanovené speciálními algoritmy „učení“ sítě Nejčastěji používané aktivační funkce f (tangens hyperbolický, sigmoidální funkce, znaménková funkce sgn). Vše je implementováno např. v MATLABu. Nejčastěji používané aktivační funkce f (tangens hyperbolický, sigmoidální funkce, znaménková funkce) Nejčastěji používané aktivační funkce f (tangens hyperbolický, sigmoidální funkce, znaménková funkce)

10. Neuronové sítě NAP2 Modeling of wheat soaking using two artificial neural networks (MLP and RBF)Journal of Food Engineering, Volume 91, Issue 4, April 2009, Pages 602-607M. Kashaninejad, A.A. Dehghani, M. Kashiri V tomto článku se dočtete jak byly vyhodnoceny experimenty zvlhčování obilních zrn. V experimentech byla zrna zvlhčována v destilované vodě při teplotách 25, 35, 45 55 a 65 0C po dobu cca 15 hodin (v 15 minutových intervalech byly váženy vzorky), celkem bylo k dispozici 154 hodnot měrné vlhkosti zrn pro různé teploty a časy. Z těchto hodnot bylo 99 dat použito pro trénování sítě, která měla dva neurony ve vstupní vrstvě (čas, teplota), 26 neuronů ve skryté vrstvě a jediný výstupní neuron (vlhkost). Zbývajících 55 dat (vlhkostí) bylo použito k ověření toho, že „natrénovaná“ neuronová síť dává rozumné výsledky a jaká je asi chyba predikce. Použily se dva typy sítě MLP (Multi Layer Perceptron) a RBF (Radial Base Functions) lišící se aktivačními funkcemi neuronů (sigmoidální, resp. Gaussovská bázová funkce). Takto vypadá výsledek predikce neuronové sítě (ANN). MR-Moisture Ratio jako funkce času a teploty. MLP je klasická neuronová síť. Aktivační funkce neuronů (tangens hyperbolický,…viz předchozí folie) nemají žádné nastavitelné parametry, optimalizují se jen váhové koeficienty wij propojující neuron i s neuronem j (a používají se stejné metody jako u dále popisovaných regresních modelů). Radiální bázové funkce neuronů RBF mají své vlastní nastavitelné parametry – souřadnice centra neuronu (určuje „vzdálenost“ centra neuronu od jednotlivých neuronů předchozí vrstvy) a „šířka“ bázové funkce. RBF jsou Gaussovské funkce RBF sítě mají jen jedinou skrytou vrstvu a váhové koeficienty wij se nastavují jen mezi skrytou a výstupní vrstvou. Parametry „radiálních“ neuronů (xc, c) se volí „ad hoc“ dle povahy problému, využívají se přitom statistické strategie typu „cluster analysis“. Je to složitější než MLP a výsledek (alespoň u sušení) bývá horší.

11. Regresní modely NAP2 Delvaux

12. Regresní modely NAP2 Regresní model má podobu relativně jednoduché funkce nezávisle proměnné x, a parametrů p1, p2,…,pM, které je třeba určit tak, aby se hodnota funkce co nejlépe shodovala s experimentálními hodnotami y. Upřímně řečeno neuronové sítě jsou téměř to samé, hledané parametry p1, p2,…,pM jsou koeficienty synaptických vah propojení neuronů. Rozdíl je v tom, že typ modelové funkce je u ANN víceméně unifikovaný a těch váhových koeficientů bývá mnohem více než počet parametrů běžně používaných regresních funkcí. Regresní funkce se volí na základě zkušeností nebo zjednodušených představ o modelovaném ději (požaduje se třeba rozumné a logicky vysvětlitelné chování při velmi malých nebo velkých hodnotách nezávisle proměnných). Je rovněž žádoucí aby parametry p1, p2,…,pM měly jasně definovaný fyzikální smysl. Nicméně je pravdou, že když neznáme fyzikální podstatu procesu nebo je příliš složitá (tj. stejná situace jako u neuronových sítí) volí se regresní funkce jako polynom y= p1+p2x+…+pMxM-1 nebo nějaká jiná „neutrální“ funkce.

13. Regresní modely NAP2 Uvažujme regresní model s hledaným vektorem M parametrů p Tyto parametry musí být stanoveny tak, aby se predikce modelu co nejvíce shodovala s N naměřenými daty (xi,yi). Nejčastěji používaným kriteriem shody je kriterium chi-kvadrát (součet čtverců odchylek) Regrese má smysl jen když je počet dat N vyšší než počet hledaných parametrů modelu M (kdyby byl stejný šlo by o interpolaci) Veličina i je směrodatná odchylka naměřené veličiny yi (chyba měření). Problém regrese je tedy optimalizační úloha hledání minima funkce 2 v prostoru parametrů p1, p2,…,pM (ale používají se i jiná robustnější kriteria shody, např. součet absolutních hodnot odchylek). Kvalita zvoleného modelu se hodnotí tzv. korelačním indexem, který by měl být blízký jedničce Hodnota r=0 odpovídá situaci, kdy by bylo lepší nahradit regresní model konstantou (průměrem naměřených dat )

14. Regresní modelylineární regrese NAP2 Lineární regresní model je lineární kombinace zvolených bázových funkcí, což mohou být třeba polynomy gm(x)=xm-1 nebo goniometrické funkce gm(x)=sin(mx) Na první pohled by se mohlo zdát, že bázové funkce jsou analogií aktivačních funkcí neuronů a parametry pm odpovídají synaptickým vahám. Jenže aktivační funkce všech neuronů jsou stejné, zatímco bázové funkce musí být rozdílné a lineárně nezávislé. Predikci modelu můžeme vyjádřit v maticovém tvaru [[A]] je návrhová matice jejíž N řádků odpovídá N-datům, a M sloupců bázovým funkcím. N>M, případ čtvercové matice N=M není regrese, ale interpolace, koeficienty [p] se stanoví řešením soustavy M rovnic pro M neznámých. U regrese je víc rovnic (N) než neznámých (M). Pokud je směrodatná odchylka měřených dat u všech bodů stejná, je 2 úměrné součtu čtverců odchylek vektoru predikce [ypredikce] a naměřených dat [y] Maticově je to vlastně skalární součin vektoru odchylek

15. Regresní modelylineární regrese NAP2 Podmínkou minima funkce s2 jsou nulové první derivace vzhledem ke všem parametrům modelu což je soustava M lineárních algebraických rovnic pro vektor neznámých parametrů Transponovaná matice [[A]]T má rozměr MxN a kdyžje vynásobena vektorem [y] Nx1 je výsledkem vektor M-hodnot Transponovaná matice [[A]]T má rozměr MxN a kdyžje vynásobena maticí [[A]] NxM je výsledkem čtvercová matice s rozměrem MxM [[C]]MxM Matice soustavy [[C]] umožňuje i odhad intervalu spolehlivosti vypočtených parametrů p1,…,pM pk je směrodatná odchylka vypočteného parametru pk pro k=1,2,…,M y je směrodatná odchylka měřených dat (předpokládáme, že je ve všech bodech stejná) Matici [[C]] je ovšem třeba invertovat (říká se jí kovarianční matice parametrů modelu)

16. Regresní modely příklad – filtrace šumu NAP2 Lineární regresní modely (používající obyčejné polynomy nízkého stupně) se výborně hodí pro vyhlazování zašuměného signálu. Princip (filtr Savitzky Golay) je jednoduchý: každému bodu vstupních dat xi,yi se přiřadí Nw bodů vlevo a Nw bodů vpravo (okénko), a těmito 2Nw+1 body se lineární regresí proloží polynom (k-tého stupně, kde k<2Nw). Hodnotou tohoto polynomu v bodě xi se pak nahradí původní hodnota yi. Počet dat N=1024, šířka okénka 50, kvadratický regresní polynom. V MATLABu je tato regresní filtrace implementována jako funkce SGOLAYFILT(X,K,F) kde X je vstupní vektor zašuměných dat, K je stupeň regresního polynomu a F=2Nw+1 je šířka okna. Application in extruder processing hmed(k)=median(h(i1:i2)); ... dpdt(1:nd)=sgolayfilt(dpd(1:nd),1,101);

17. Regresní modely příklad - zvlhčování NAP2 Pro modelování závislosti absorbované vlhkosti zrn na čase se používá několik různých regresních modelů, např. exponenciální Pageův model Model má 2 parametry p1=k, p2=n. Nezávisle proměnná t je čas Xe je rovnovážná, XO počáteční vlhkost zrna a ještě častěji používaný Pelegův model (viz článek týkající se fazolí) Parametr p2 charakterizuje rovnovážnou vlhkost zrna The application of Peleg's equation to model water absorption during the soaking of red kidney beans (Phaseolus vulgaris L.)Journal of Food Engineering, Volume 32, Issue 4, June 1997, Pages 391-401Nissreen Abu-Ghannam, Brian McKenna Pelegův model se používá pro luštěniny, soju, ořechy, cizrnu, hrách… Ani Pageův ani Pelegův model není lineární a tak se trochu „fixluje“. Pageův model lze převést na lineární dvojím logaritmováním, Pelegův model ještě jednodušeji úpravou na tvar To je teď nová závisle proměnná y …a pak se použije lineární regrese pro stanovení p1 a p2

18. Regresní modely zvlhčování NAP2 Na předchozí folii bylo uvedeno, že linearizací modelu se trochu „fixluje“. Jde o to, že když použijeme lineární regresí stanovené parametry v původním nelineárním modelu, nebude dosaženo minima součtu čtverců odchylek od měřené vlhkosti (a našla by se asi lepší kombinace parametrů p1, p2). V praxi se tento rozdíl zanedbává i když je trochu komické, že i pak se povinně počítají kovarianční matice, korelační index a podobné statistiky, jejichž význam je poněkud sporný.

19. Analytická řešenímodely difúze a tepla NAP2 Vermeer

20. Analytická řešení NAP2 Analytická řešení hledejte především u lineárních rovnic Obyčejné diferenciální rovnice: Parciální diferenciální rovnice: existují dva základní přístupy všechny analytické funkce, které znáte (i neznáte), např. sin, cos, exp, tgh, Besselovy funkce,… jsou vlastně definovány jako řešení určitých typů diferenciálních rovnic (řešení hledejte v příručkách, např. Kamke E: Differential Gleichungen, Abramowitz M., Stegun I.: Handbook of Mathematical functions). Obecný postup spočívá ve vyjádření řešení ve tvaru mocninové řady, jejíž koeficienty získáte dosazením do diferenciální rovnice. Fourierova metoda separace proměnných. Řešení F(t,x,y,z) se hledá ve tvaru součinu funkcí, které závisí jen na jediné proměnné F=T(t)X(x)Y(y)Z(z). Dosazením do parciální diferenciální rovnice získáte obyčejné diferenciální rovnice. Na diferenciální rovnici aplikujte integrální transformaci: Fourierovu, Laplaceovu, Hankelovu. Výsledkem je algebraická rovnice, kterou je třeba vyřešit a pak použít zpětnou transformaci.

21. Analytická řešenídifúze NAP2 Rozložení vlhkosti X (kg vody/kg sušiny) popisuje pro případ konstantního součinitele difúze vody v materiálu Fickova rovnice Xe je rovnovážná vlhkost, X0 je počáteční vlhkost m=0,1,2 pro desku, válec, kouli. Pro všechny tři základní geometrie lze nalézt analytické řešení Fourierovou metodou separace proměnných. Řešení se hledá ve tvaru součinu funkcí polohy r a času t Dosazením do Fickovy parciální diferenciální rovnice získáme obyčejné diferenciální rovnice pro funkce Gi(t) i Fi(r*) Nemůže to záviset ani na r ani na t, takže to musí být konstanta. Konstanty  (tzv. vlastní hodnoty) se volí tak, aby řešení splňovalo okrajovou podmínku na povrchu (například nulovou hodnotu X pro r=R). Toto závisí pouze na čase t Toto může záviset pouze na r

22. Analytická řešenídifúze NAP2 Řešením rovnice pro Gi(t) je ve všech třech případech exponenciála, řešení rovnice pro Fi(r) je cos(r) pro desku (m=0), Besselova funkce J0(r) pro válec (m=1) a pro kouli (m=2, což nás u kulových zrnek zajímá nejvíc) Pokud je na povrchu koule (r*=1) udržována fázová rovnováha (X=Xe), tj. když je na povrchu zajištěn dostatečně intenzivní přenos hmoty, je X*(t,1)=0 a tuto podmínku musí splňovat každá funkce Fi. A to je podmínka pro vlastní hodnoty Bezrozměrný koncentrační profil má tedy tvar Zatím neurčené koeficienty ci musí být stanoveny tak, aby řešení splňovalo počáteční podmínku, tj. rozložení koncentrace v čase nula. Pro konstantní koncentraci X=X0 tedy musí pro každé r*platit X*=1

23. Analytická řešenídifúze NAP2 Pro stanovení koeficientů ci využijeme zajímavou vlastnost funkcí Fi(r*), které se říká ortogonalita. Dvě funkce jsou ortogonální, když jejich skalární součin je roven nule. Skalární součin je v případě funkcí definován jako integrál jejich součinu (násobený event. nějakou váhovou funkcí) Ortogonalita funkcí je výborná vlastnost, která se uplatní i v řadě jiných aplikací. Když např. použijeme při regresi jako bázové funkce gm(x) ortogonální polynomy všechno se zjednoduší a zpřesní (nemusí se řešit soustavy rovnic, zmizí problém zaokrouhlovacích chyb) Důkaz ortogonality vlastních funkcí Fi To jsou vlastně definiční diferenciální rovnice pro funkce Fi a Fj. Můžeme je odečíst a integrovat. Integrace per partes. A to je nula, protože obě vlastní funkce musí splňovat okrajové podmínky (nulová derivace.v centru koule a nulová hodnota na povrchu)

24. Analytická řešenídifúze NAP2 Aplikujme ortogonalitu na předchozí rovnici (vynásobenou r2Fj a integrovanou) Prosím, zkuste ověřit, že to tak je (je to jen malé cvičení integrace per partes) Výsledkem je tedy koncentrační profil vlhkosti nebo (integrací přes objem koule) celková vlhkost jako funkce času

25. Analytická řešenídifúze NAP2 Difúzní součinitel Def ovšem závisí na teplotě, dokonce i na vlhkosti a stejně tak i rovnovážná vlhkost Xe je funkcí teploty. Je tedy třeba vzít do hry ještě další korelace, např. I když je použito analytické řešení, je výsledkem silně nelineární regresní model se třemi parametry p1=Dwa, p2=Ea a p3=b. Tento model použila Katrin Burmester při hledání optimálního režimu pražení kávových zrn Heat and mass transfer during the coffee drying processJournal of Food Engineering, Volume 99, Issue 4, August 2010, Pages 430-436 Katrin Burmester, Rudolf Eggers Poznamenejme, že i tak je použití analytického řešení jen aproximací (protože i při izotermním experimentu je difúzní součinitel Def proměnný) a i předchozí článek uvádí kromě analytického řešení, numerické řešení metodou sítí uvažující i proměnný teplotní profil v praženém zrnu kávy. Tomu se budeme věnovat později. Modeling and simulation of heat and mass transfer during drying of solids with hemispherical shell geometryComputers & Chemical Engineering, Volume 35, Issue 2, 9 February 2011, Pages 191-199I.I. Ruiz-López, H. Ruiz-Espinosa, M.L. Luna-Guevara, M.A. García-Alvarado

26. Optimalizace NAP2 • Pro hledání parametrů pi silně nelineárních modelů, které minimalizují 2, se používají dva základní typy metod: • Bezderivační, kterým pro stanovení minima stačí jen možnost vyčíslování s2 (nebo jiného zvoleného kriteria shody predikce a experimentu) pro libovolné hodnoty parametrů modelu. Metody tohoto typu jsou nutné u velmi komplikovaných regresních modelů, například když se predikce počítá metodou konečných prvků. • Derivační, které využívají i hodnoty všech prvních (někdy i druhých) derivací regresního modelu vzhledem ke všem parametrům p1,p2,…,pM.

27. Optimalizační metody derivační NAP2 Gaussova metoda hledání minima součtu vážených čtverců odchylek Váhové koeficienty mohou být třeba převrácené hodnoty variancí (čtverců očekávaných chyb) model data Podmínkou minima jsou nulové první derivace vzhledem ke všem parametrům j=1,2,…,M Regresní funkci je možné linearizovat, aproximovat ji pouze lineárními členy Taylorova rozvoje p je jen vektor přírůstku parametrů modelu v jedné iteraci

28. Optimalizační metody derivační NAP2 V každé iteraci Gaussovy metody se tedy řeší soustava lineárních rovnic pro vektor přírůstku parametrů V současnosti asi nejčastěji používaná metoda Marquardt Levenberg je mírnou modifikací Gaussovy metody: když iterace nekonvergují, přidává se k diagonále matice [[C]] konstanta . Když je  velmi velké stává se matice [[C]] diagonální a Gaussova metoda přechází na metodu gradientní (pravá strana je vlastně jen vektor gradientu minimalizované funkce s2). Hodnota  se během iterací mění, když iterace konvergují,  se zmenšuje a preferuje se velice rychlá Gaussova metoda, když divergují, tak se  zvyšuje (gradientní metoda je sice pomalejší, ale spolehlivější).

29. Optimalizační metody bezderivační NAP2 Nejjednodušší je případ optimalizace jediného parametru (M=1). Nejprve je třeba lokalizovat optimální hodnotu parametru do nějakého intervalu ve kterém je jen jediné minimum (interval neurčitosti). Přesná poloha minima se pak určí buď metodou půlení intervalu (ale na rozhodnutí, zda řešení leží „vlevo“ nebo „vpravo“ jsou třeba dvě hodnoty regresní funkce) nebo metodou zlatého řezu (Golden Section, interval neurčitosti se v každém kroku zkracuje sice jen v poměru 0.618 a ne v poměru 0.5, ale stačí zase jen jedna hodnota regresní funkce). Viz algoritmus Golden section search a následující folie(definice zlatého řezu) f1 f4 f3 f2 L1 L3=0.618L2 L2=0.618L1 V původním intervalu neurčitosti (délky L1) se vypočtou f1 f2 ve zlatých řezech. Porovnáním f1>f2plyne, že minimum nemůže být vlevo. V novém intervalu neurčitosti L2 opět potřebujeme dvě hodnoty f pro porovnání, ale jedna z nich (f2) už je spočítaná, protože bod ve zlatém řezu starého intervalu je i ve zlatém řezu nového intervalu – to je právě ta základní vlastnost zlatého řezu a magie čísla 0.618. Stačí tedy výpočet jediné hodnoty f3!

30. a B C2 C1 A p q Zlatý řez NAP2 kvadratická rovnice pro poměr zlatého řezu q/p

31. Optimalizační metody bezderivační NAP2 Problém do testu: Kolik kroků metody zlatého řezu a kolik funkčních hodnot regresní funkce je třeba vyčíslit, když požadujeme zmenšení intervalu neurčitosti třeba tisíckrát? Pro tisícinásobné zpřesnění (m=1000) tedy stačí cca 14 kroků. An approach to determine diffusivity in hardening concrete based on measured humidity profilesAdvanced Cement Based Materials, Volume 2, Issue 4, July 1995, Pages 138-144 D. Xin, D. G. Zollinger, G. D. Allen V tomto článku se řeší difúzní rovnice (ale ne zrna, tuhnutí betonu) a pro stanovení difúzního součinitele je použita metoda zlatého řezu (regresní model je definován jako numerické řešení difúzní rovnice-kolokační metoda).

32. Optimalizační metody bezderivační NAP2 Metody jednorozměrného hledání lze aplikovat i na problémy s víceparametry p1,…,pM , prostě se opakuje jednorozměná optimalizace pro každý parametr zvlášť (Rosenbrock). Asi nejčastěji používaná je simplexová metoda Nelder Mead . Její princip je prostý: 1. Generuje se simplex tvořený náhodně zvolenými M+1 vrcholy (pro dva parametry p1 p2 je to trojúhelník). 2. V každém vrcholu simplexu se vyčíslí regresní funkce f a vybere se nejhorší vrchol. Ten se nahradí buď překlopením, kontrakcí nebo expanzí. Krok 2 se opakuje dokud se simplex dostatečně nezmenší Animovaný gif převzatý z wikipedia.org

33. MATLAB NAP2

34. Optimalizační metody MATLAB NAP2 V MATLABU je snadná lineární polynomická regrese p=polyfit(x,y,m) a minimalizace funkce více proměnných (bez omezení „constraints“, metoda Nelder Mead) p=fminsearch(fun,p0) kde fun je odkaz (“handle” začínající symbolem @) na uživatelem definovanou funkci, která pro hodnoty parametrů modelu předá hodnotu, která má být mimalizována (třeba součet čtverců odchylek experimentálních dat od matematického modelu). Vektor p0 jsou zadávané počáteční hodnoty – nástřel. Pro nelineární regresi (nelineární modely) jsou k dispozici p = nlinfit(x,y,modelfun,p0) ci = nlparci(p,resid,'covar',sigma)umožňuje statistické zpracování výsledků (kovarianční matice, intervaly spolehlivosti predikce i počítaných parametrů modelu)

35. Optimalizační metody MATLAB NAP2 Příklad: Byla naměřena sušící křivka, 10 bodů (čas a vlhkost) matice a vektory jde definovat výčtem v lomených závorkách []. Prvky, které tvoří řádek matice jsou odděleny mezerami, středník znamená přechod na nový řádek. Apostrof [vektor]’ znamená transpozici, místo řádkového bude výsledkem sloupcový vektor time X(moisture) 0 0.9406 1 0.7086 2 0.7196 3 0.5229 4 0.4657 5 0.3796 6 0.3023 7 0.1964 8 0.1545 9 0.1466 xdata=[0 1 2 3 4 5 6 7 8 9]’; ydata=[0.9406 0.7086 …]’; data vykreslíme příkazem plot(xdata,ydata,’*’)

36. Optimalizační metody MATLAB NAP2 Body sušící křivky můžeme proložit polynom třetího stupně p1x3+…+p4 p=polyfit(xdata,ydata,3) p = 0.0001 0.0040 -0.1322 0.9103 vektor vypočtených koeficientů polynomu (ve směru klesajících mocnin) Polynom s koeficienty p(1)=0.0001, p(2)=0.0040 …vykreslíme např. takto Y=polyval(p,xdata) hold on plot(xdata,ydata,'*') plot(xdata,Y) vyčíslení hodnot polynomu do vektoru Y hold on umožňuje kreslit více křivek do jednoho grafu to celé šlo udělat i bez použití vektoru predikce Y jediným příkazem plot(xdata,ydata,'*',xdata,polyval(p,xdata))

37. Optimalizační metody MATLAB NAP2 Fyzikálně správnější model sušící křivky může být založen na difúzním modelu • kde A je „škálovací“ koeficient a fyzikální význam má difúzní koeficient Def a poloměr charakteristické částice R. Def a R ovšem nejsou nezávislé parametry, protože v modelu figuruje jen jejich poměr Def/R2, převrácená hodnota časové konstanty sušení. Pokud tedy poloměr R zvolíme (např. R=0.01 m) je cílem optimalizace nalezení dvou regresních parametrů A, Def, které minimalizují zvolené kriterium shody modelu a naměřených bodů sušící křivky (součet čtverců odchylek, nebo „robustnější“ kriterium součtu absolutních hodnot odchylek, které není tak citlivé na chybně naměřené body, tzv. „outliers“). Jak to udělat? • Definuj funkci xmodel(t,p,pa) jako součet řady (vektor parametrů p1=A, p2=Def, pa1=R, pa2=n) • Definuj funkci, která spočítá součty odchylek xdev(xmodel,p,pa,xdata,ydata) • použij knihovní funkci pro hledání minima p= fminsearch(xdev,p0)

38. Optimalizační metody MATLAB NAP2 Takže nejprve definice modelu function xval = xmodel(t,p,pa) A=p(1); D=p(2); R=pa(1); ni=pa(2); xv=0; for i=1:ni pii=(pi*i)^2; xv=xv+6/pii*exp(-pii*D*t/R^2); end xval=A*xv; tento text musí být uložen jako M-file se jménem souboru stejným jako je název funkce, tedy xmodel.m Takto lze definovat jakýkoliv model sušení, třeba i ten předchozí kubický polynom, Pelegův či Pageův model nebo model popisovaný diferenciálními rovnicemi. Model založený na diferenciálních rovnicích by však musel byl integrován pro každou hodnotu času t znovu – pro takové případy by byla efektivnější jiná strategie, vyhodnocení predikce modelu ve všech časech t (v celém vektoru xdata) současně.

39. Optimalizační metody MATLAB NAP2 Funkce definující zvolené kriterium shody (např. součet čtverců odchylek) function sums = xdev(model,p,paux,xdat,ydat) sums=0; n=length(xdat); for i=1:n sums=sums+(model(xdat(i),p,paux)-ydat(i))^2; end tento text musí být uložen jako M-file se jménem souboru stejným jako je název funkce, tedy xdev.m Definice pomocných parametrů a hledání regresních parametrů funkcí fminsearch R=0.01 a počet členů řady n=10 pa=[0.01,10] p = fminsearch(@(p) xdev(@xmodel,p,pa,xdata,ydata),[.5;0.0003]) druhý parametr funkce fminsearch je vektor nástřelu první parametr fminsearch může být pouze funkce s jediným parametrem, vektorem p optimalizovaných parametrů. To, že pro vyčíslení odchylky potřebujeme zadat i další hodnoty se v MATLABu řeší použitím tzv. anonymní funkce @(p) expression.

40. Optimalizační metody MATLAB NAP2 Přesně to samé lze naprogramovat i do jediného M-filu function [estimates, model] = fitcurve(xdata, ydata) start_point = [1 0.00005]; Hledané parametry jsou dva: A,D. Počáteční odhad. model = @expfun; Expfun je predikce modelu a výpočet součtu čtverců odchylek estimates = fminsearch(model, start_point); volání optimalizační Nelder Mead metody function [sse, FittedCurve] = expfun(params) A = params(1); optimalizovaný škálovací parametr D = params(2); optimalizovaný difuzní koeficient R=0.01; poloměr částice (když ho chcete změnit musíte opravit funkci fitcurve.m) ni=10; počet členů řady(správně nekonečno, ale 10 obvykle stačí) ndata=length(xdata); (počet bodů naměřené křivky sušení) sse=0; výsledek expfun – součet čtverců odchylek for idata=1:ndata xv=0; tady je třeba naprogramovat konkrétní model sušení (Y(X,params)) fori=1:ni pii=(pi*i)^2; xv=xv+6/pii*exp(-pii*D*xdata(idata)/R^2); end FittedCurve(idata) = A* xv; ErrorVector(idata) = FittedCurve(idata) - ydata(idata); sse=sse+ErrorVector(idata)^2; end end end [estimates, model] = fitcurve(xdata,ydata)

41. Co je třeba si pamatovat NAP2 Minimum toho, co byste si z této přednášky věnované především minimalizaci pamatovat (alespoň ke zkoušce), je na následující folii

42. Co je třeba si pamatovat NAP2 Triáda bilančních rovnic Lagrangeův variační princip (minimum Wi+We) Kriterium shody mezi modelem a realitou (chi-kvadrat) Princip regrese (parametry p minimalizující chi-kvadrat) Jak minimalizovat chi-kvadrat, když je regresní model nelineární funkcí parametrů p Minimalizační metody derivační (Marquardt Levenberg) a bezderivační (zlatý řez, simplexová metoda Nelder Mead)