Matematický model

2.2  Simulace

Obsah kapitoly:

2.2.1  Princip numerického řešení obyčejných diferenciálních rovnic 2.2.2  Přesnost numerického řešení 2.2.3  Stabilita numerického řešení 2.2.4  Stanovení hodnot parametrů modelu 2.2.5  Simulační programovací jazyky 2.2.6  Příklady 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ů.

2.2.1  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

Soubor:Obrázek1.png

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ů  ti , pro něž platí

ti+1 = ti + h

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

(2.16)

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

Soubor:Obrázek2.png

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

(2.17)

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

Soubor:Obrázek3.png

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ě

Soubor:Obrázek4.png

Soubor:Obrázek5.png

(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

Soubor:Obrázek6.png

(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  Yi  ,  i=1, 2,..., k  na  t.

2.2.2  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

Soubor:Obrázek7.png

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.

Soubor:Obrázek8.png

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:

#

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

1. nalezneme řešení s tímto zvoleným krokem,
2. krok zmenšíme, např. na polovinu,
3. nalezneme řešení se zmenšeným krokem,
4. 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ší,
5. 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.

2.2.3  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  

Soubor:Obrázek9.png

Soubor:Obrázek10.png

(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	Soubor:Obrázek11.png
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.

 

na začátek kapitoly

 

2.2.4  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 yE(t). Stejný signál u(t) přivedeme na vstup modelu a zaznamenáme výstupní signál  yM(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

Soubor:Obrázek12.png

(2.21)

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

Soubor:Obrázek13.png

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

Soubor:Obrázek14.png

(2.23)

nebo obecněji

Soubor:Obrázek15.png

(2.24)

kde v1 a v2  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ů.

2.2.5  Simulační programovací jazyky

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.

 

na začátek kapitoly

 

2.2.6  Příklady simulace

Příklad S-1

Úkol:

Řešit diferenciální rovnici (2.20), tj.

Soubor:Obrázek16.png

Soubor:Obrázek18.png

Eulerovou metodou s krokem 0,2 pomocí jazyka PSI a simulací získaný výsledek porovnat s analytickým řešením (rovněž naprogramovaným v jazyku PSI) výpočtem chyby. Totéž řešení provést pomocí metody Runge-Kutta 2.řádu a výsledky porovnat.

 

Výpis programu:	% Příklad S-1 (řešení dif. rovnice) y     =  INT(-2*TIME*y PAR:1);    % numerické řešení yan  = exp(-(TIME^2));                  % analytické řešení ch   =  y-yan;                                   % chyba
Ukázka výsledků: Eulerova metoda, krok 0.2,  max.chyba 7.5 . 10-2
Metoda Runge-Kutta 2.řádu, krok 0.2, max.chyba -5 . 10-3
Program v PSI:	soubory: Prikl_S-1.PSM , Prikl_S-1.MOD

 

na začátek kapitoly

 

Příklad S-2

Úkol:

Simulovat reakci výstupní teploty elektrického ohřívače vody (viz příklad M-1) na skokovou změnu teploty vstupující vody. Matematický model ohřívače je

Soubor:Obrázek21.png

Soubor:Obrázek22.png

Hodnota teploty T0 odpovídá ustálenému stavu.

Doba simulace je 60 minut, skok vstupní teploty je proveden v 5.minutě. Číselné hodnoty parametrů a jejich rozměry jsou patrné z výpisu programu. K numerickému řešení je použito metody Runge-Kutta 2.řádu, krok výpočtu je 0,001.

Výpis programu:

% Příklad S-1 (simulace reakce T2 na skokovou změnu T1) T2 = INT(Q*(T1-T2)/V+P*eta/(V*rocp) PAR:T0); % teplota vody na výstupu (°C) V = 5;                                         % objem ohřívače (l) Q = 0.5;                                      % průtok vody (l/min) T0 = 37.91;                                 % výst.teplota v ust.stavu (°C) rocp= 1*4187;                            % součin hustoty a měr.tepla vody (J/l/K) T1 = BNG(time PAR:15,18,5);    % skok teploty na vstupu (°C) P = 60000;                                  % příkon topení (J/min) eta = 0.8;                                     % účinnost topení (1)

Ukázka výsledků: Soubor:Obrázek23.png

Program v PSI:

soubory: Prikl_S-2.PSM , Prikl_S-2.MOD

Z výsledku je vidět, že po změně vstupu se teplota výstupu ustaluje na nové hodnotě - v čase 59.04 minut je to  40.9092 °C (viz graf).

 

na začátek kapitoly

 

Příklad S-3

Úkol:

Provést identifikaci modelu elektrického ohřívače vody (viz příklad M-1) za předpokladu, že není známa hodnota účinnosti topení  h. Pro identifikaci byl proveden experiment, kdy byl po zapnutí topení naměřen následující časový průběh teploty :

čas (min)	0	5	10	15	20	25	30	35	40	45	50	55	60
teplota (°C)	15	26	31	34	37	39	40	40	40	40	41	41	41

Simulační model ohřívače je popsán v příkladu S-2, program je ještě třeba doplnit o příkaz umožňující zpracovávat experimentální data a příkaz pro výpočet optimalizačního kritéria. Protože model dává výsledky prakticky spojité, použijeme kritéria podle rovnice (2.22), kde meze integrálu budou 0 a 60 min, místo yM bude teplota T2M  vypočtená podle modelu a místo yE  naměřená teplota T2E. Výsledky experimentu musí být v tom případě k dispozici rovněž jako spojité, což zajistí příkaz FNG lineární interpolací v zadané tabulce. Odhadnutá hodnota (nástřel) účinnosti bude např. 0.8. Pro určení její skutečné hodnoty použijeme optimalizační metody, která je součástí jazyka PSI. Výsledný program je následující:

Výpis programu:

% Příklad S-3 (identifikace modelu ohřívače) % výst.teplota - model(°C) T2M = INT(Q*(T1-T2M)/V+P*eta/(V*rocp) PAR:T1);  V = 5;                                            % objem ohřívače (l) Q = 0.5;                                         % průtok vody (l/min) rocp= 1*4187;                               % součin hust. a měr.tepla vody (J/l/K) T1 = 15;                                         % teplota vody na vstupu (°C) P = 60000;                                     % příkon topení (J/min) eta = PAR(0.8);                              % účinnost (optimalizační proměnná) T2E = FNG(time PAR:0,60,1,T2mer);  % výst. teplota - experiment (°C) krit= INT((T2M-T2E)^2 PAR:0);         % kritérium optimalizace

Nastavení parametrů optimalizace:

1. položka menu Optimization / Optimization Data

Criterion:	krit	(jméno bloku, kde se počítá hodnota kritéria)
Initial Stepsize:	0.2	(počáteční krok změny optimalizační proměnné při vyhledávání)
Stop Accuracy:	0.0001	(přesnost pro ukončení optimalizace)
Maximum No of Iterations:	50	(maximální počet vyhledávacích kroků)
Optimization Method:	Pattern Search + PT	(metoda vyhledávání s předčasným ukončením)

2. položka menu Optimization / Variables to Optimize 

Variable:	eta	(jméno optimalizační proměnné - blok PAR nebo CON)
Value:	0.8	(počáteční hodnota optimalizační proměnné  - "nástřel")
Min:	0.5	(dolní mez optimalizační proměnné)
Max:	1.0	(horní mez optimalizační proměnné)
Gain:	1	(násobná konstanta - sjednocení řádů při více optimalizačních proměnných)

Optimalizace se spouští položkou menu Optimization / Start Optimization, průběh se vypisuje v textovém okně.

Výsledek simulace před optimalizací

Výsledek simulace po optimalizaci

Program v PSI:

soubory: Prikl_S-3.PSM , Prikl_S-3.MOD

Hodnota kritéria se optimalizací změnila z 367.924 na 7.26312 , výsledná hodnota účinnosti je 0.89858, čili prakticky 0.9.