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

MODELY PROCESů NAP2

Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

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ů,...).

Matematický model aparátu nebo systému

Identifikace parametrů modelu na základě experimentu

Vyhodnocení experimentů – diagnostika aparátů, např. diagnostika kolon, EEG

Optimalizace operačních parametrů (optimalizace návrhu)

MODELY PROCESů NAP2

Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

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.

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

Dvg

Dt

2( )( ) ( )

2gD v

u q v QDt

0D

vDt

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

hustota tepelného

toku

výkon vnitřních

sil

Většina rovnic, které znáte, se z těchto rovnic dá odvodit: Bernoulliho, Eulerova, Fourierova Kirchhoffova Navier Stokesova.

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

1( ) :

2iW u e d

iie uFdupdufuW

)(

2

0

1( )

2

t

i eS W W v d dt

deformační energie vnitřních sil (součin tenzoru

deformace a napětí)

práce vnějších sil (objemové,

povrchové, osamělé)u

Wi

We

W(u)

deformace i napětí lze vyjádřit jako funkce

posunutí u

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í).

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ů.

Příklady modelů-sušení NAP2

Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

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ů.

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

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ů w i (koeficientů zesílení signálů

mezi propojenými neurony)

NEURONOVÁ SÍŤ - MOZEK Neuron má několik vstupů - dendrity

a jen jeden výstup - axon

NEURONOVÁ SÍŤ - MOZEK Neuron má několik vstupů - dendrity

a jen jeden výstup - axon

Řízení průtoku substrátu

Řízení průtoku substrátu

Sensory 7

6

5

4

3

2

1

Vstupní vrstva (zde 2 neurony)

X1,X2

Skrytá vnitřní vrstva (zde 4 neurony)

Výstupní vrstva (zde 1 neuron) Y

1 1

1

1

( ) ( )

1/ (1 exp( ( ))

sgn( )

N N

i i i ii i

N

i ii

N

i ii

y f w x tgh w x

w x

w x

Odezva neuronu y na N-vstupních hodnot xi

Nejčastěji používané aktivační funkce f (tangens hyperbolický, sigmoidální funkce,

znaménková funkce)

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)

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.

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 w ij 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 w ij 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ší.

)2

)(exp()(

2

2

c

cc

xxxf

Regresní modely NAP2

Delvaux

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.

Regresní modely

),(),....,,( 21 pxfpppxfy M

N

i i

ii pxfyp

1

22 )),(

()(

NAP2

2

2

( ( , ))1

( )i i

i

y f x pr

y y

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)

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

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)

Hodnota r=0 odpovídá situaci, kdy by bylo lepší nahradit regresní model

konstantou (průměrem naměřených dat )y

Regresní modely

M

mmm xgppxf

1

)(),(

1 1 1

1

( ) .. ( )

...

...

( ) ( )

M

N M N

g x g x

A

g x g x

[ ] [[ ]][ ]predikcey A p

NAP2

2 ([ ] [[ ]][ ]) ([ ] [[ ]][ ])Ts y A p y A p

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)

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

Regresní modely lineární regrese

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é.

NAP2

0]]][[[]][[2]])[[][2(][

2

pAAAyp

s TTT

[[ ]] [[ ]][ ] [[ ]] [ ]T TA A p A y

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

maticí [[A]] NxM je výsledkem čtvercová matice s rozměrem MxM

[[C]]MxMMatice soustavy [[C]] umožňuje i odhad intervalu spolehlivosti vypočtených parametrů p1,…,pM

Transponovaná matice [[A]]T má rozměr MxN a když je vynásobena vektorem [y] Nx1

je výsledkem vektor M-hodnot

2 1 2pk kk yC

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)

Regresní modely lineární regrese

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);

Regresní modely příklad - zvlhčování NAP2

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

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

0

exp( )ne

e

X Xkt

X X

Xe je rovnovážná, XO počáteční vlhkost

zrna

Model má 2 parametry p1=k, p2=n. Nezávisle

proměnná t je čas

a ještě častěji používaný Pelegův model (viz článek týkající se fazolí)

01 2

tX X

p p t

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

Parametr p2 charakterizuje rovnovážnou vlhkost zrna

Pelegův model se používá pro luštěniny, soju, ořechy, cizrnu, hrách…

1 20( )

tp p t

X t X

…a pak se použije lineární regrese pro stanovení p1 a p2

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ý.

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

Vermeer

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.

* ** * *

2 * * *0

( ) , ef m em

e

D X XX X rr X r

t R r r r X X R

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

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 * * *

1

( , ) ( ) ( )i i ii

X t r c F r G 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*)

2* 2

* * *

1( )mi i

imef i i

dG dFR dr

D G dt r F dr dr

Toto závisí pouze na čase t Toto může záviset

pouze na 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).

Xe je rovnovážná vlhkost, X0 je počáteční vlhkost

NAP2 Analytická řešení difúze

Ř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)

2

2*

**

sin( )( ) , ( )

iefD t

i Ri i

rF r G t e

r

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

sin( ) 0i i i

Bezrozměrný koncentrační profil má tedy tvar 2 2

2*

* **

1

sin( )( , )

efi Dt

Ri

i

i rX t r c e

r

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 *

*1

sin( )1 i

i

i rc

r

NAP2 Analytická řešení difúze

NAP2 Analytická řešení difúze

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í)

1*2 * *

0

( ) ( ) 0i jr F r F r dr 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

*2 2 *2* *

*2 2 *2* *

( ) 0

( ) 0

ij i i j

ji j i j

dFdF r r FFdr dr

dFdF r r FFdr dr

To jsou vlastně definiční diferenciální rovnice pro funkce Fi a Fj. Můžeme je

odečíst a integrovat.

1 1*2 *2 * 2 2 *2 *

* * * *0 0

[ ( ) ( )] ( )jij i j i i j

dFdFd dF r F r dr r FF drdr dr dr dr

1

*2 1 2 2 *2 *0* *

0

[ ( )] ( )jij i j i i j

dFdFr F F r FF dr

dr dr

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)

NAP2 Analytická řešení difúze

Aplikujme ortogonalitu na předchozí rovnici (vynásobenou r2Fj a integrovanou)

1 1* * * * * *

10 0

1* * *

01

2 * *

0

sin( ) sin( )sin( )

sin( )2( 1)

sin ( )

ii

j

j

r j r dr c i r j r dr

r j r dr

cj

j r dr

Prosím, zkuste ověřit, že to tak je (je to jen malé cvičení integrace per partes)

*

*1

sin( )1 i

i

i rc

r

2 2

2*

* **

1

2( 1) sin( )( , )

efi Di tR

i

i rX t r e

ir

2 2

2*2 2

1

6( )

efi Dt

R

i

X t ei

Výsledkem je tedy koncentrační profil vlhkosti nebo (integrací přes objem koule) celková vlhkost jako

funkce času

exp( )exp( )aef wa

ED D bX

RT

Heat and mass transfer during the coffee drying process  Journal of Food Engineering, Volume 99, Issue 4, August 2010, Pages 430-436 Katrin Burmester, Rudolf Eggers

Modeling and simulation of heat and mass transfer during drying of solids with hemispherical shell geometry  Computers & 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

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

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.

NAP2 Analytická řešení difúze

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.

Optimalizační metody derivační NAP2

Gaussova metoda hledání minima součtu vážených čtverců odchylek

i

iii wfys 22 )(

0)(22

ii

j

iii

j

wp

ffy

p

s

0)( 0

i

ij

i

kk

k

iii w

p

fp

p

ffy

Váhové koeficienty mohou být třeba převrácené hodnoty variancí (čtverců očekávaných chyb)

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

data model

Optimalizační metody derivační NAP2

jkk

jk BpC

i

ik

i

j

ijk w

p

f

p

fC

i

iiij

ij wfy

p

fB )( 0

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ší).

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

f2f3

f4

L1 L2=0.618L1 L3=0.618L2

V původním intervalu neurčitosti (délky L1) se vypočtou f1 f2 ve zlatých řezech. Porovnáním f1>f2 plyne, ž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!

Zlatý řez NAP2

A BC1 C2 a

p q

qp

p

p

q

1)1( p

q

p

q

618.02

51

p

qkvadratická rovnice pro poměr zlatého

řezu q/p

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?

0.618

log log 0.618

log

0.209

nLL

mm n

mn

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 profiles  Advanced 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).

Optimalizační metody bezderivační NAP2

Metody jednorozměrného hledání lze aplikovat i na problémy s více parametry 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ý:

Animovaný gif převzatý z wikipedia.org

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ší

MATLAB NAP2

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)

Optimalizační metody MATLAB NAP2

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 …]’;

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

0 1 2 3 4 5 6 7 8 90.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

data vykreslíme příkazem plot(xdata,ydata,’*’)

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

0 1 2 3 4 5 6 7 8 90.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Polynom s koeficienty p(1)=0.0001, p(2)=0.0040 …vykreslíme např. takto

Y=polyval(p,xdata)

hold onplot(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))

vektor vypočtených koeficientů polynomu (ve směru klesajících mocnin)

Optimalizační metody MATLAB NAP2

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)

2 2

2*2 2

1

6( )

efi Dn tR

i

X t A ei

Fyzikálně správnější model sušící křivky může být založen na difúzním modelu

Optimalizační metody MATLAB NAP2

2 2

2*2 2

1

6( )

efi Dn tR

i

X t A ei

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);endxval=A*xv;

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ě.

tento text musí být uložen jako M-file se jménem souboru stejným jako je

název funkce, tedy xmodel.m

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

pa=[0.01,10]

p = fminsearch(@(p) xdev(@xmodel,p,pa,xdata,ydata),[.5;0.0003])

R=0.01 a počet členů řady n=10

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.

druhý parametr funkce fminsearch je vektor

nástřelu

Optimalizační metody MATLAB NAP2

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ů odchylekestimates = 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)) for i=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 endend

[estimates, model] = fitcurve(xdata,ydata)

Přesně to samé lze naprogramovat i do jediného M-filu

Co je třeba si pamatovatNAP2

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

Co je třeba si pamatovatNAP2

Dvg

Dt

2( )( ) ( )

2gD v

u q v QDt

0

Dv

Dt

Triáda bilančních rovnic

Lagrangeův variační princip (minimum Wi+We)

1( ) :

2iW u e d

iie uFdupdufuW

)(

Kriterium shody mezi modelem a realitou (chi-kvadrat)

N

i i

ii pxfyp

1

22 )),(

()(

Princip regrese (parametry p minimalizující chi-kvadrat) [[ ]] [[ ]][ ] [[ ]] [ ]T TA A p A y

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)