Modelovani a simulace

Modelování a simulace

V praxi často potřebujeme zjišťovat chování daného procesu nebo zařízení za různých podmínek, ale nemůžeme s ním z technologických nebo ekonomických důvodů provádět složitější experimenty. Efektivním, levným a z technologického hlediska bezrizikovým postupem je v takovém případě vytvoření matematického modelu procesu a studium jeho chování pomocí metod simulace.

Matematické modely jsou však vhodné nejen jako náhrada experimentů, ale také pro řízení technologie a pro předvídání průběhu procesu a jeho výsledku. Každý technolog by proto měl znát základní principy odvozování matematických modelů, neboť ho v této činnosti nemůže nikdo nahradit. Tato práce totiž vyžaduje dokonalou znalost technické a technologické podstaty modelovaných dějů, kterou má především technolog.

Vytváření matematických modelů

Stav procesu je charakterizován hodnotami veličin, které jej popisují (tzv. stavové veličiny). Jsou to teplota, tlak, množství či koncentrace složek apod. Matematickým modelem procesu nebo zařízení rozumíme matematické vztahy mezi těmito veličinami, nejčastěji v závislosti na čase. Hovoříme pak o tzv. dynamice procesu. Matematickými vztahy jsou nejčastěji diferenciální rovnice (obyčejné i parciální) nebo jejich soustavy, ale také nelineární rovnice a jejich soustavy, soustavy lineárních rovnic apod. Pro účely automatického řízení se často užívá matematického popisu ve tvaru přenosové funkce, což má některé výhody, jak se dozvíme později.

Matematické modely můžeme podle způsobu vytváření rozdělit na dvě základní kategorie: modely stochastické a modely deterministické (nebo také fyzikální).

Stochastické modely předpokládají, že modelovaný systém je „černá skříňka“ s definovanými vstupy a výstupy. Vytváříme je zhruba následujícím postupem :

• na základě předpokládaného chování systému zvolíme vhodný tvar matematického popisu (často to bývá přenosová funkce), který obsahuje parametry neznámých hodnot,

• na vstup systému přivádíme proměnlivý signál (může být skutečný nebo uměle generovaný) a v pravidelných (krátkých) časových intervalech měříme hodnoty vstupního a výstupního signálu tak dlouho, až získáme dostatečně dlouhé řady odpovídajících dvojic hodnot,

• výsledky měření zpracujeme tak, abychom získali hodnoty parametrů zvoleného matematického popisu systému.

Tohoto postupu používáme tehdy, když nejsme schopni odvodit deterministický model, nebo při aplikaci některých speciálních algoritmů automatického číslicového řízení.

Deterministické modely jsou odvozovány na základě analýzy a matematického popisu skutečné podstaty děje. Těmito modely se budeme podrobněji zabývat a shrneme hlavní zásady a postupy jejich vytváření.

Obecný postup vytváření deterministických matematických modelů

Proces (systém), který chceme matematicky popsat, se řídí podle určitých fyzikálních, fyzikálně-chemických a chemických zákonů, které mají své matematické vyjádření. Z tohoto vyjádření zákonů při sestavování deterministických modelů vycházíme. Někdy takto můžeme popsat všechny zadané podmínky a vztahy modelovaného procesu úplně a získat tak přesný matematický model. Takový matematický model ale bývá tak složitý, že jeho řešení je prakticky neproveditelné. Navíc většinou ani úplný popis získat nemůžeme, protože průběh děje do potřebných podrobností neznáme. V praxi ale obvykle nepotřebujeme naprosto přesné výsledky, stačí, když model vystihuje podstatné vlastnosti a chování procesu. Můžeme si tedy dovolit některé méně podstatné vlivy a vztahy zanedbat nebo zjednodušit. Konečným kritériem kvality a použitelnosti modelu je vždy jeho souhlas s realitou v souladu s účelem ke kterému byl vytvořen. Na obr. 2.1  jsou schématicky znázorněny jednotlivé fáze vytváření deterministického modelu.

Obr.2.1 Postup vytváření matematického modelu

První fází je důkladná analýza procesu, při níž musíme provést následující kroky:

1. a)specifikujeme děje, které v procesu probíhají a určíme jejich podstatu,
2. b)vymezíme vlivy, které na proces a jeho průběh působí,
3. c) určíme veličiny, které proces popisují,
4. d) o každém z výše uvedených prvků rozhodneme, zda jej budeme považovat za podstatný a do modelu jej zahrneme, či nikoliv,
5. e)o každém z výše uvedených prvků zahrnutých do modelu rozhodneme, zda jej budeme brát v úvahu v celé šíři a ve všech souvislostech, nebo zda použijeme nějaká zjednodušení.

Výsledkem popsaného postupu je tzv. teoretický model. Ten sice nepopisuje zcela přesně skutečnost, ale jeho výhoda spočívá v jednoduchosti, přehlednosti a následně i ve snazším řešení výsledných rovnic.

V bodu e) se hovoří o zjednodušování některých vztahů při tvorbě teoretického modelu. Je to je velmi důležitý krok, který musíme vždy provádět a navíc má podstatný vliv na kvalitu výsledného modelu. Uvádíme proto výčet nejpoužívanějších zjednodušujících předpokladů:

• rozdělujeme složitý systém na jednodušší subsystémy a zpracováváme je odděleně, ovšem se zřetelem na jejich vzájemné souvislosti,
• zavádíme neexistující formy, např. ideální plyn, ideální kapalina,
• předpokládáme nezávislost látkových vlastností na teplotě,
• předpokládáme, že materiál je homogenní a isotropní, 
• u systémů, kde probíhají současně pomalé a rychlé děje, předpokládáme, že rychlý děj již dosáhl rovnovážného stavu,
• zanedbáváme ztráty, 
• nelineární závislosti linearizujeme,
• zavádíme různé korekční (obvykle empirické) koeficienty,
• používáme empiricky zjištěných vztahů a závislostí mezi veličinami,
• zjednodušujeme geometrické proporce a volíme vhodné souřadnicové soustavy, abychom mohli co nejjednodušeji popsat geometrii systému,
• zaváděním některých předpokladů odstraňujeme závislost sledovaných veličin na souřadnicích - získáváme tím tzv. modely se soustředěnými parametry.

Ve druhé fázi postupu sestavujeme na základě teoretického modelu matematický popis procesu. Znamená to provést následující kroky:

1. a)výběr matematického popisu zákonitostí použitých v teoretickém modelu, 
2. b) vytvoření modelových rovnic, tj. doplnění vybraných vztahů o zjednodušující předpoklady a potřebné matematické úpravy,
3. c)určení podmínek řešení (obvykle počáteční a okrajové podmínky pro řešení diferenciálních rovnic).

V této fázi obvykle využíváme matematických rovnic vyjadřujících známé zákony a vztahy z fyziky, fyzikální chemie, chemie atd. Výsledkem tohoto postupu je obecný matematický model procesu. Může se však také stát, že náš teoretický model matematicky popsat nedokážeme; potom musíme buď změnit teoretický model (zavést jiné nebo další předpoklady), nebo v krajním případě musíme použít vhodných experimentálních postupů a jejich vyhodnocením získat náhradní popis děje (např. stochastický model).

Třetí fází celého postupu je řešení modelu, to znamená vytvoření simulačního programu. Představuje to provést tyto kroky:

1. a) volbu metody řešení modelových rovnic,
2. b) zpracování modelových rovnic pro řešení na počítači, tj. nalezení vhodného algoritmu řešení,
3. c) sestavení výpočetního programu, nejlépe s využitím vhodného simulačního jazyka.

Čtvrtou fází je vytvoření tzv. simulačního modelu. Jejím výsledkem je počítačový program vhodný pro používání v praxi. Tato fáze zahrnuje následující kroky:

1. a) identifikace modelu, tj. nalezení neznámých hodnot parametrů modelu (např. porovnáním získaných výsledků řešení s údaji z literatury, s experimentálními hodnotami apod.), 
2. b)verifikace modelu, tj. řešení kontrolních úloh a analýza jejich výsledků za účelem ověření správnosti modelu v celé předpokládané oblasti použití, posouzení přesnosti a vhodnosti modelu pro daný účel. 

Pokud model vyhovuje, práci jsme úspěšně dokončili a můžeme modelu používat. Jestliže nevyhovuje, musíme zjistit, v čem jsme se dopustili chyby, který předpoklad byl neoprávněný, které zjednodušení nepřípustné, který vliv jsme opominuli nebo neoprávněně zanedbali. Musíme zformulovat nový teoretický model a celý postup znovu zopakovat.

Je zřejmé, že tvorba matematického modelu není nijak jednoduchou záležitostí. Má však jednu velkou přednost: přinutí nás důkladně si ujasnit podstatu a mechanismus studovaného děje.

Na závěr bychom rádi zdůraznili jednu velmi důležitou zásadu: je třeba začínat vždy od teoretických modelů maximálně jednoduchých a komplikovat je teprve tehdy, když výsledky nevyhovují našim představám a požadavkům nebo praktickým zkušenostem.

Vytváření matematických modelů na základě bilancí

Bilancemi rozumíme aplikaci zákonů o zachování na daný systém. Jejich využití při vytváření matematických modelů je velice časté - odhaduje se, že je to 70 až 80 % případů. Bilancovat lze jen tzv. extenzivní veličiny. V praxi se jedná nejčastěji o zákony o zachování hmoty a energie - hovoříme pak o látkových a energetických bilancích.

Chceme-li se zabývat bilancováním, musíme si nejprve vyjasnit některé základní pojmy. Je to především bilancovaný systém. Rozumíme jím nějakou (nejčastěji prostorovou) oblast, která má definované hranice (rozhraní) s okolím. Okolí je tedy to, co systém obklopuje. Další pojem je bilanční časový interval, což je časový úsek, během kterého bilancované veličiny sledujeme. Základní bilanční rovnice má tvar

VSTUP + ZDROJ = VÝSTUP + AKUMULACE

(2.1)

VSTUP  je množství bilancované veličiny, které za bilanční časový interval vstoupí z okolí přes rozhraní do systému.

VÝSTUP  je množství bilancované veličiny, které za bilanční časový interval vystoupí ze systému přes rozhraní do okolí.

ZDROJ  je množství bilancované veličiny, které za bilanční časový interval uvnitř systému vznikne nebo zanikne (rozumí se vznik z jiné, nebilancované veličiny, nebo změna na jinou, nebilancovanou veličinu).

AKUMULACE  je změna množství (zádrže) bilancované veličiny uvnitř systému za bilanční časový interval.

Při bilancování musíme nejprve zvolit velikost bilančního časového intervalu. Volba je dána jednoduchým pravidlem:

• pokud se bilance týká ustáleného stavu systému (tj. hodnoty veličin se s časem nemění), pak můžeme volit časový interval zcela libovolný, nejčastěji jednotkový,
• pokud chceme modelovat systém jako proměnný v čase (tzv. dynamický model), pak musíme volit bilanční časový interval infinitesimálně malý, tj. dt.

Dalším důležitým krokem je volba hranic systému. Vycházíme při ní jednak ze skutečného (konstrukčního) tvaru systému, jednak z toho, zda chceme vytvořit model systému se soustředěnými parametry (v matematickém popisu nefigurují souřadnice jako nezávisle proměnné), nebo s rozloženými parametry (alespoň jedna souřadnice je v matematickém popisu použita jako nezávisle proměnná). Pravidlo pro volbu velikosti bilancovaného systému je opět (podobně jako u časového intervalu) jednoduché:

• jedná-li se o model se soustředěnými parametry, volíme rozměry bilancovaného systému libovolně (např. podle konstrukčního tvaru aparátu),
• jedná-li se o model s rozloženými parametry, pak rozměr bilancovaného systému ve směru té souřadnice, na níž problém závisí (figuruje v popisu jako nezávisle proměnná), musí být infinitesimálně malý (dx apod.).

S volbou rozhraní úzce souvisí volba vhodné souřadnicové soustavy. Tu se snažíme volit vždy tak, aby získaný matematický model byl co nejjednodušší. Základní souřadnicové soustavy jsou kartézská (pravoúhlá), cylindrická (válcová) a sférická (kulová). Jejich popis najdete v matematické literatuře. Při volbě se řídíme následujícími zásadami:

• nejčastěji používáme kartézské souřadnice (jako nejobecnější),
• pro problémy souměrné kolem osy používáme cylindrické souřadnice,
• pro středově souměrné problémy používáme sférické souřadnice.

Velice důležitá je volba kladných směrů ve vybrané souřadnicové soustavě. Jestliže ji důsledně neprovádíme a nedodržujeme, vznikají neřešitelné znaménkové problémy vedoucí k nesprávným výsledkům. Pokud však jsme v jejich dodržování důslední, nemusíme se o znaménka vstupu, výstupu a akumulace vůbec starat - vycházejí automaticky. Fyzikální zákony už totiž ve svém matematickém vyjádření správné znaménko zahrnují. Znaménko musíme úvahou určit jedině pro zdrojový člen: jedná-li se o vznik veličiny, je jeho znaménko kladné, jedná-li se o zánik, je záporné.

Je třeba si uvědomit, že všechny odvozené rovnice musí být rozměrově homogenní. Doporučujeme tuto skutečnost vždycky zkontrolovat.

Počáteční a okrajové podmínky

U modelů popsaných diferenciálními rovnicemi musíme popis doplnit příslušným počtem okrajových a počátečních podmínek. Pro každou nezávisle proměnnou potřebujeme tolik vzájemně nezávislých podmínek, jaký je nejvyšší v rovnicích se vyskytující řád derivace podle této proměnné. Formulace počátečních a okrajových podmínek je nedílnou součástí vytváření matematického modelu. Některé podmínky vyplývají zcela jednoduše ze zadání úlohy (např. na počátku je teplota ve všech bodech stejná a rovná určité hodnotě), jiné musíme odvodit stejnými postupy jako matematický model (např. na základě bilance). Jako kontrola správnosti jejich odvození nám může sloužit skutečnost, že obecně co do matematického tvaru existuje jen několik druhů podmínek, a tedy v konkrétním případě musíme dosáhnout shody s jedním z nich.

Pro popis obecného matematického tvaru druhů počátečních a okrajových podmínek použijeme následující označení veličin:

u ... závisle proměnná,              t ... čas,            x, y, z ... souřadnice,  

f ... funkční předpis, jehož tvar známe a hodnotu funkce můžeme kdykoliv vypočíst.

Počáteční podmínka je obvykle jedna a definuje situaci na počátku řešení. Obecně ji můžeme zapsat ve tvaru

pro  t = t0:               u = f (x,y,z)

(2.2)

čili na počátku (v čase t₀) je závisle proměnná u známou funkcí souřadnic x, y, z.

Okrajové podmínky se vyskytují v případech, kde jako nezávisle proměnné vystupují souřadnice. Rozeznáváme tři základní typy okrajových podmínek :

1. a)okrajová podmínka 1.druhu (Dirichletova) :

pro x = x0 :                      u = f (y,z,t)

(2.3)

1. hodnota závisle proměnné v místěx₀ je známou funkcí ostatních souřadnic a času.
2. b)okrajová podmínka 2.druhu (Neumannova) :

pro x = x0 :

(2.4)

čili hodnota derivace závisle proměnné podle jedné souřadnice (např. podle x  v bodě x₀)  je známou funkcí ostatních souřadnic a času. Často se setkáváme s okrajovou podmínkou 2. druhu ve tvaru

pro x = x0 :

(2.5)

1. derivace závisle proměnné podle normály k nějaké ploše je nulová. 
2. c)okrajová podmínka 3.druhu (Newtonova) :   

pro x = x0 :

(2.6)

1. hodnota lineární kombinace hodnoty závisle proměnnéuv bodě x₀  a její derivace podle x v místě x₀ je známou funkcí ostatních souřadnic a času; konstanty a, b jsou koeficienty lineární kombinace. 

Ukázkové příklady na vytváření matematických modelů

Příklad M-1

Úkol: Vytvořit matematický model elektrického průtokového ohřívače vody za účelem sledování časového průběhu výstupní teploty.

Známe:

objem ohřívače V, průtok vody Q, vstupní teplotu vody T₁, příkon topení P, účinnost topení h , 
hustotu  vody r , měrné teplo vody  c_p

Předpoklady (tj. teoretický model):

Model budeme odvozovat pro systém se soustředěnými parametry (předpoklad dokonalého promíchávání obsahu ohřívače, tj. všude uvnitř je teplota stejná jako na výstupu), hustota a měrné teplo jsou nezávislé na teplotě.

Bilance:

Jako bilancovaný systém volíme celý ohřívač (nezávislost na souřadnici), jako bilanční časový interval volíme dt (hledáme model pro neustálený stav). Bilancovaná veličina bude energie (resp. entalpie). Jednotlivé členy bilanční rovnice jsou:

	VSTUP :	QrcpT1.dt
	VÝSTUP :	QrcpT2.dt
	ZDROJ :	Ph.dt
	AKUMULACE :	Vrcp.dT2

Dosazením do základní bilanční rovnice (2.1) a úpravou dostaneme výslednou diferenciální rovnici pro teplotu na výstupu ohřívače

(2.7)

Počáteční podmínku zformulujeme např. z předpokladu, že teplota na výstupu se na počátku rovná teplotě na vstupu, tedy

pro t = 0 :

T2(0) = T1

(2.8)

Rovnice (2.7) a (2.8) tvoří hledaný matematický model.

Příklad M-2

Úkol: Vytvořit matematický model pro stanovení časových změn teplotního průběhu podél teploměrné jímky vyčnívající ze stěny zařízení do proudícího média.

Známe:

rozměry jímky (délku L, průřez S, obvod B), teplotu proudícího média T_M, teplotu stěny T_S, 
tepelnou vodivost materiálu jímky l, hustotu materiálu jímky r, měrné teplo materiálu jímky c_p, 
koeficient přestupu tepla mezi povrchem jímky a proudícím médiem a

Předpoklady (tj. teoretický model):

Teplota je funkcí pouze souřadnice x a času t (tj. zanedbáváme rozdělení teploty na průřezu jímkou), teplo se vyměňuje s proudícím médiem na celém povrchu rovnoměrně, koncem jímky (x=L) se teplo nepředává, hodnoty fyzikálních parametrů jímky jsou v uvažovaném teplotním rozmezí konstantní.

Bilance:

Zanedbáváme rozložení teploty na průřezu jímky, takže vystačíme s jedinou souřadnicí x. Kladný směr osy x je zleva doprava. Jako bilancovaný systém volíme element („plátek“) vyříznutý z jímky kolmo k její délce, jeho délka je dx  (protože nás zajímá závislost na souřadnici), vzdálenost tohoto elementu od počátku je obecně x. Jako bilanční časový interval volíme dt (protože nás zajímá i závislost na čase). Bilancujeme teplo a jednotlivé členy bilanční rovnice jsou:

	VSTUP :
	VÝSTUP :	vedením                                  konvekcí
	ZDROJ :	0
	AKUMULACE :

První člen ve výrazu pro výstup (vedením) upravíme provedením naznačené derivace a po zanedbání derivace  dx podle x (diference 2.řádu) dostaneme

Pak dosadíme do základní bilanční rovnice  (2.1)  a upravíme. Výsledkem je hledaná parciální diferenciální rovnice pro teplotu T jako funkci souřadnice x a času t :

(2.9)

K této rovnici potřebujeme jednu počáteční a dvě okrajové podmínky. Počáteční podmínka definuje stav na počátku, tj. teplota jímky ve všech bodech se rovná teplotě stěny:

pro t = 0 :

T(x,0) = TS

(2.10)

Okrajová podmínka pro x=0 vychází ze zadání, že teplota jímky se trvale rovná teplotě stěny, tedy

pro x = 0 :

T(0,t) = TS

(2.11)

a okrajová podmínka pro konec jímky x=L se odvodí na základě předpokladu, že se koncem jímky teplo nepředává, tedy gradient teploty je nulový:

pro x = L :

(2.12)

Hledaný matematický model je tedy tvořen rovnicemi (2.9 - 2.12). Jen pro úplnost uvádíme, že rovnice (2.11) je okrajová podmínka 1.druhu (viz rov. 2.3) a rovnice (2.12) je okrajová podmínka 2.druhu (viz rov. 2.4).

 Příklad M-3

Úkol: Vytvořit matematický model zásobníku kapaliny určený pro sledování výšky hladiny kapaliny v něm v závislosti na změnách přítoku a odtoku. Na odtoku z nádrže je instalován ventil, jehož míru otevření je možné měnit, přítok do nádrže může kolísat.

Známe:

tvar zásobníku (stojatý válec), plochu dna A, přítok kapaliny Q₁ a charakteristiku ventilu (vztah vyjadřující závislost průtoku Q na míře otevření x, na tlakovém spádu na ventilu Dp a hustotě kapaliny r; K_V je konstanta ventilu)

Předpoklady (tj. teoretický model):

Pro náš případ je tlakový spád na ventilu dán hydrostatickým tlakem, tj. Dp = rgh, zanedbávají se tlakové ztráty prouděním kapaliny, hustota kapaliny je konstantní.

Bilance:

Jako bilancovaný systém volíme celou nádrž, bilancovanou veličinou je množství kapaliny, což lze v našem případě vyjádřit pouze objemem. Jednotlivé členy bilanční rovnice jsou:

	VSTUP :	Q1.dt
	VÝSTUP :
	ZDROJ :	0
	AKUMULACE :	A.dh

Dosazením do základní bilanční rovnice (2.1) a úpravou dostaneme výslednou diferenciální rovnici pro hladinu v nádrži

(2.13)

Počáteční podmínku zformulujeme např. z předpokladu, že výška hladiny na počátku odpovídá ustálenému stavu (přítok = odtok)

pro t = 0 :

h(0) = h0

(2.14)

Rovnice (2.13) a (2.14) tvoří hledaný matematický model.

Simulace

Simulací (pro účely matematického modelování) rozumíme řešení rovnic matematického modelu, přičemž měníme podle potřeby hodnoty vstupních signálů a parametrů a sledujeme jejich vliv na hodnoty výstupních signálů. U modelů dynamických systémů to znamená řešit diferenciální rovnice nebo jejich soustavy. Postupů analytického řešení známých z matematiky zde lze také použít, ale spíše v jednodušších případech. Většinou volíme řešení numerické. Navíc jsou k dispozici simulační jazyky, jejichž používání je většinou snadné a nevyžaduje příliš hluboké znalosti z programování a práce s počítačem. Abychom však dovedli výsledky správně interpretovat, je rozumné znát základní principy, na kterých je numerické řešení diferenciálních rovnic založeno, a je nutné vědět něco o přesnosti výsledků.

Princip numerického řešení obyčejných diferenciálních rovnic

Budeme hovořit o numerickém řešení obyčejných diferenciálních rovnic (ODR) 1.řádu s počáteční podmínkou, které jsou lineární (tj. derivace je v nich vždy v první mocnině). Princip řešení si ukážeme na nejjednodušší metodě Eulerově, ale je společný všem tzv. vícekrokovým metodám. Mějme diferenciální rovnici obecně zapsanou ve tvaru

s počáteční podmínkou   y(t0) = y0

(2.15)

(t ... čas, obecně nezávisle proměnná, y ... závisle proměnná, g ... funkční předpis, jehož tvar známe). Každou diferenciální rovnici 1.řádu na tento tvar můžeme převést osamostatněním derivace. Princip řešení spočívá v tom, že spojitý interval nezávisle proměnné na kterém hledáme řešení diskretizujeme, tj. rozdělíme jej na stejně velké díly o velikosti h, takže dostaneme řadu bodů  t'_i , pro něž platí

ti+1 = ti + h

pro  i =  0, 1, 2, ...

(2.16)

Hodnoty závisle proměnné v bodech t_i  pak vypočteme z Eulerovy formule

pro  i =  0, 1, 2, ...

(2.17)

[]

kde hodnota Y_i+1  je přibližně rovna hledanému řešení  y(t_i+1), přičemž na počátku Y₀ = y₀ .

Pro názornost je užitečný grafický výklad principu Eulerovy metody, který je vidět na obr.2.2. Metoda vychází z náhrady průběhu analytického řešení na úseku mezi dvěma sousedními body ti , ti+1 přímkou. V počátečním bodě (t0, y0) jsme schopni vypočíst hodnotu pravé strany rovnice (2.13) g(t0, y0), tedy vlastně hodnotu derivace dy/dt (směrnice tečny). Vzorec (2.17) pak udává, že bod (t1, Y1) získáme jako průsečík tečny v počátečním bodě a pořadnice t1. Takto získaný bod (t1, Y1) můžeme nyní považovat za počáteční a postup opakovat. Můžeme tak postupně získat body (t2, Y2) , (t3, Y3) atd. až do konce intervalu t  na němž hledáme řešení. Výsledkem numerického řešení diferenciální rovnice není tedy vzorec jako u řešení analytického, ale tabulka bodů (ti, Yi) ,   i = 0, 1, 2,...

Obr.2.2  Princip Eulerovy metody

Postup řešení (kreslení tohoto obrázku) je možné si prohlédnout 

zde

.

Na stejném principu je založeno řešení soustavy ODR 1.řádu s počátečními podmínkami. Soustava obecně k diferenciálních rovnic má  k  závisle proměnných a jednu společnou nezávisle proměnnou. Rovnici (2.15) si můžeme představit zapsanou vektorově

(2.18)

kde je y vektor závisle proměnných, y' vektor derivací závisle proměnných avektor funkcí (pravých stran po osamostatnění derivací). Stejně tak můžeme vektorově zapsat Eulerovu formuli

(2.19)

takže vlastně dostaneme samostatný Eulerův vzorec pro každou závisle proměnnou a v každém kroku řešení vypočteme pro danou hodnotu nezávisle proměnné  t  hodnoty všech závisle proměnných podle těchto vzorců. Výsledkem je tabulka vyjadřující závislosti proměnných  Y_i  ,  i=1, 2,..., k  na  t.

Přesnost numerického řešení

Numerické řešení diferenciální rovnice závisí nejen na nezávisle proměnné t , ale i na velikosti kroku h . Aby přibližné řešení Y odpovídalo analytickému (přesnému) řešení y , musí platit podmínka konvergence

tedy přesnost numerického řešení je tím větší, čím menší je krok h . U krokových metod vyjadřujeme přesnost resp. chybu metody pomocí mocniny h a hovoříme o metodě příslušného řádu. V případě Eulerovy metody je to přesnost 1.řádu h. Tento údaj nám slouží pouze k hrubému odhadu chyby, tedy např. je-li  h = 0.1 , dá se očekávat, že přesné bude první desetinné místo výsledku.

O odhad se jedná především proto, že chyba závisí nejen na velikosti kroku, ale i na tvaru funkce v okolí právě počítaného bodu. Z toho logicky vyplývá, že chyba numerického řešení se může v jeho jednotlivých částech (pro různá t) lišit. Na obr. 2.2 jsou chyby numerického řešení v bodech t1 a t2 označeny jako D1 a D2. Tato chyba je chyba metody a protože souvisí s diskretizací nezávisle proměnné, říká se jí chyba diskretizační. Při numerickém řešení se však do výpočtů vnáší ještě další chyba související se skutečností, že výpočetní prostředek (počítač) pracuje pouze na omezený počet platných míst. Tato chyba souvisí také s velikostí kroku h, protože v rovnici (2.17) může docházet při násobení hodnoty funkce  g  krokem  h  ke značným řádovým posunům. Je to chyba zaokrouhlovací a její velikost je tím větší, čím menší je krok h.  Celková chyba je pak součtem obou, jak je znázorněno na obr. 2.3.

Obr.2.3  Chyba Eulerovy metody

Diskretizační chyba u Eulerovy metody závisí na kroku lineárně, ale obecně tomu tak není. Zmíněné dvě chyby působí proti sobě. Abychom zjistili skutečnou chybu numerického řešení konkrétní rovnice (a to v praxi obvykle potřebujeme), musíme postupovat empiricky následujícím způsobem:

1. podle požadované přesnosti a řádu přesnosti zvolené metody numerického řešení ODR odhadneme potřebnou velikost kroku h,

2. nalezneme řešení s tímto zvoleným krokem,

3. krok zmenšíme, např. na polovinu,

4. nalezneme řešení se zmenšeným krokem,

5. hodnoty závisle proměnné z obou řešení ve stejných bodech nezávisle proměnné porovnáme a předpokládáme, že správné jsou ty dekadické řády výsledků, které se u obou řešení neliší,

6. pokud nejsme se získanou přesností spokojeni, krok dále zmenšíme a porovnáváme vždy výsledky dvou posledních kroků; to děláme tak dlouho, až nám přesnost vyhovuje.

[]

V praxi se používají metody s vyšší přesností, které pracují na podobném principu jako metoda Eulerova, ale mají poněkud složitější vzorce. K nejznámějším patří metody Runge-Kutta, a to 2.řádu a 4.řádu. Zájemci najdou jejich popis v každé příručce numerické matematiky, vzorce uvádíme 

zde

Stabilita numerického řešení

Vedle problémů s chybou numerického řešení diferenciálních rovnic mohou ještě vyvstat problémy s jeho stabilitou. Je-li numerické řešení stabilní, tak se výše popsané chyby vnášené do výpočtu v každém kroku vhodně utlumují a nerostou. Naopak, při nestabilitě se tyto chyby kumulují a mohou narůst tak, že výsledky jsou zcela nepoužitelné. Na stabilitu řešení má opět vliv velikost kroku h, a to tak, že s rostoucím krokem se pravděpodobnost nestability zvyšuje. Pro každou diferenciální rovnici existuje krok při kterém nestabilita začíná - hovoříme o mezním kroku stability řešení. Problém nestability je ilustrován na obr. 2.4., kde jsou tabulka a graf numerického řešení rovnice  

(2.20)

získané Eulerovou metodou s různým krokem. Při kroku 0.2 je řešení stabilní a blíží se řešení analytickému, při kroku 0.5 je sice řešení ještě stabilní, ale informace o průběhu funkce je velmi hrubá a při kroku 1.0 je řešení nestabilní, kmitá a výsledky jsou zcela nesmyslné - z tvaru analytického řešení jednoznačně vyplývá, že výsledky musí být z intervalu (0;1ñ .

t''' i  Y''' i   Y''' i Y''' i    y''' anal (h=0.2) (h=0.5) (h=1.0) 0 1 1 1  1  0.2 1.000 0.961 0.4 0.920 0.852 1.5 1.000 0.779 0.6 0.773 0.698 0.8 0.587 0.527 1.0 0.399 0.500 1.000 0.368 1.2 0.240 0.237 1.4 0.125 0.141 1.5 0.000 0.105 1.6 0.055 0.077 1.8 0.020 0.039 2.0 0.006 0.000 -1.000 0.018 :  : :  :  :  3.0 0.000 0.000 3.000 0.000
Obr. 2.4.  Ukázka k problému stability numerického řešení ODR

Pokud je řešení nestabilní, pomůže ve většině případů zmenšení kroku.

Z toho, co bylo řečeno o chybách a stabilitě, vyplývá, že nemůžeme numerické řešení matematických modelů provádět mechanicky, ale musíme nad ním přemýšlet a vypracovat si vždy způsob jak výsledky kontrolovat. Protože lze předpokládat, že budeme většinou řešit problémy z praxe, je nám velkou pomocí právě možnost kontroly výsledků srovnáním s realitou, s fyzikální podstatou veličin a s praktickými omezeními jejich hodnot.

Stanovení hodnot parametrů modelu

Při odvozování matematického modelu nebo i jen při vyjadřování některých dílčích závislostí veličin v modelu často používáme parametrů (konstant), jejichž hodnoty předem neznáme a musíme je nějakým způsobem stanovit. Hovoříme o tzv. identifikaci modelu. Problematika identifikace představuje dost rozsáhlou oblast, my se zde budeme zabývat pouze nejjednodušší metodikou, kterou bychom mohli běžně použít. V matematice se probírá postup aproximace empirických dat metodou nejmenších čtverců (prokládání přímky či jiné funkce naměřenými daty apod.). Tato metodika je sice obecně použitelná, ale nedá se využít tehdy, kdy je potřeba dosáhnout souhlasu výsledků řešení diferenciálních rovnic s nějakými empirickými daty. V tomto případě je vhodnější použít metod optimalizace.

Optimalizace  je matematický resp. početní postup, kterým hledáme extrém nějaké funkce. Této funkci říkáme účelová nebo kriterium optimalizace a její tvar musíme předem určit tak, abychom optimalizací dosáhli požadovaného cíle. Hodnotu kritéria optimalizace můžeme ovlivnit změnou zvolených parametrů matematického popisu optimalizační úlohy, kterým říkáme optimalizační proměnné. Extrém (buď maximum nebo minimum) hledáme na nějaké množině, z níž hodnoty optimalizačních proměnných vybíráme. Na této množině může být více extrémů, kterým říkáme extrémy lokální (jsou na podmnožinách dané množiny) a jeden extrém globální, který právě hledáme. Při většině metod numerické optimalizace musíme zadat nějakou počáteční hodnotu pro každou optimalizační proměnnou, tzv. nástřel. Lokální extrémy jsou  nepříjemné tím, že v nich optimalizační postupy někdy „uvíznou“ a pak nemohou dosáhnout extrému globálního. Pomoci si můžeme tak, že zvolíme jiný nástřel, ale musíme být schopni (podle znalosti řešeného problému a charakteru výsledků) rozpoznat, zda je nalezený extrém lokální či globální, což nebývá vždy jednoduché a vyžaduje to určitou zkušenost. Pokud se týče optimalizačních proměnných, může být jejich výběr co do rozsahu hodnot neomezený, tj. z intervalu (-¥,+¥), pak mluvíme o hledání volného extrému, nebo omezený, např. pro každou optimalizační proměnnou máme určen uzavřený interval z něhož mohou být její hodnoty vybírány, a pak mluvíme o hledání vázaného extrému. Ve většině praktických optimalizačních úloh hledáme vázaný extrém, protože vždy existují nějaká technická nebo technologická omezení.

  Základním problémem optimalizace je volba vhodného optimalizačního kritéria. V případech identifikace matematických modelů by mělo vyjadřovat shodu modelu se skutečností. Při identifikaci hledáme takové hodnoty parametrů modelu, aby odezva modelu na změnu vstupu byla co nejbližší odezvě reálného systému na stejnou změnu vstupu. Předpokládejme tedy, že jsme na vstup reálného systému přivedli definovaný (tj. se známým časovým průběhem) signál u(t) a na výstupu jsme naměřili časově proměnný signál y_E(t). Stejný signál u(t) přivedeme na vstup modelu a zaznamenáme výstupní signál  y_M(t). Pokud má model odpovídat skutečnosti, neměly by se časové průběhy těchto signálů podstatně lišit. Optimalizační kritérium vytvoříme tak, aby kvantitativně vyjadřovalo míru shody obou průběhů. Prakticky vždy se používá osvědčené kritérium součet čtverců odchylek, které má pro případ porovnávání obou průběhů v n bodech tvar

(2.21)

nebo v případě hodnocení shody na spojitém intervalu <a, b>  tvar

(2.22)

Tato kritéria nabývají při naprosté shodě obou průběhů hodnoty 0, jinak se nule blíží. Jestliže jich použijeme jako účelové funkce, pak optimalizací hledáme jejich minimum. Jestliže bychom potřebovali porovnávat průběh více výstupů současně, např.  y(t) a z(t), pak jednoduše použijeme účelové funkce ve tvaru (např. pro spojitý případ)

(2.23)

nebo obecněji

(2.24)

kde v₁ a v₂  jsou váhy, pomocí kterých můžeme podle potřeby upravit vliv výsledku porovnání jednotlivých průběhů na hodnotu účelové funkce. Podle stejného principu postupujeme i v případě porovnávání tří a více průběhů.

Simulační programovací jazyk

Pro řešení matematických modelů existuje celá řada simulačních programovacích jazyků, které umožňují efektivní práci bez hlubších znalostí programování počítačů. Jejich podstata je různá, např. existují jazyky rovnicově orientované (programují se zápisem rovnic modelu), blokově orientované (programují se propojováním funkčních bloků, které realizují základní operace, např. integraci, různé matematické funkce atd.), modulově orientované (obsahují moduly představující již hotové modely zařízení nebo procesů a ty se pak propojují pomocí signálů reprezentujících toky hmot a energií) apod. Mají obvykle zabudovány různé numerické metody pro řešení modelových rovnic, dobrou grafickou prezentaci výsledků simulace a možnost práce se vstupními i výstupními datovými soubory. Jsou velmi často oborově zaměřené.

My pro naši potřebu budeme používat jednoduchý, poměrně univerzální, blokově orientovaný jazyk PSI, který je k dispozici na síti v učebnách VŠCHT. Je to blokově orientovaný jazyk a byl vyvinut na univerzitě v holandském Delftu. Jeho předností je to, že je speciálně zaměřen na řešení soustav obyčejných diferenciálních rovnic. Z hlediska uživatele je velice příznivé, že je v základě zvládnutelný zhruba během dvou hodin. Součástí této učebnice je i jeho volně šiřitelná demoverze, která umožňuje řešit úlohy malého rozsahu. Popis jazyka PSI a možností demoverze je v příloze této učebnice.

Dalším simulačním nástrojem, který je k dispozici na síti VŠCHT, je rovněž blokově orientovaný SIMULINK, který je součástí programového balíku MATLAB. Jeho určitou nevýhodou je, že k jeho používání je nutné nejprve zvládnout základy MATLABu, což vyžaduje poměrně dost času.

Použité zdroje

Kadlec K., Kmínek M., Učebnice měřicí a řídicí techniky, dostupné online https://uprt.vscht.cz/ucebnice/mrt/

 (původní texty a obrázky jsou použity se souhlasem autorů)

O projektu

Kybernetika wiki byla vytvořena na 

Ústavu matematiky, informatiky a kybernetiky

 VŠCHT zejména pro studenty oboru 

Chemická kybernetika

.