VYSOKÉ UENÍ TECHNICKÉ V BRN FAKULTA STAVEBNÍ ALEŠ DRÁB HYDROINFORMATIKA I MODUL M03 ÚVOD DO MATLAB

VYSOKÉ U ENÍ TECHNICKÉ V BRN FAKULTA STAVEBNÍ

ALEŠ DRÁB

HYDROINFORMATIKA I MODUL M03 ÚVOD DO MATLAB

STUDIJNÍ OPORY PRO STUDIJNÍ PROGRAMY S KOMBINOVANOU FORMOU STUDIA

Hydroinformatika I · Modul 3

© Aleš Dráb, Brno 2005

- 2 (34) -

Obsah

OBSAH 1 Úvod ...............................................................................................................5 1.1 Cíle ........................................................................................................5 1.2 Požadované znalosti ..............................................................................5 1.3 Doba pot ebná ke studiu .......................................................................5 1.4 Klí ová slova.........................................................................................5 1.5 Metodický návod na práci s textem ......................................................6 2 Za ínáme s MATLAB ..................................................................................7 2.1 Základy programování v MATLAB .....................................................7 2.1.1 Definice prom nných..............................................................8 2.1.2 Priorita matematických operátor ...........................................9 2.2 Využití nápov dy v MATLAB .............................................................9 2.3 Práce s maticemi .................................................................................10 2.4 Import dat ............................................................................................11 2.5 Export dat a uložení MATLAB workspace ........................................13 2.6 Použití soubor typu M-file ................................................................14 2.6.1 Script M-file..........................................................................14 2.6.2 Function M-file .....................................................................16 2.6.3 Rozhodovací struktury a cykly .............................................17 3 ešení optimaliza ních úloh ......................................................................21 4 ešení oby ejných diferenciálních rovnic ................................................23 5 ešení soustav oby ejných diferenciálních rovnic ..................................27 6 Záv r ............................................................................................................33 6.1 Shrnutí.................................................................................................33 6.2 Studijní prameny .................................................................................33 6.2.1 Seznam použité literatury .....................................................33 6.2.2 Seznam dopl kové studijní literatury ...................................34 6.2.3 Odkazy na další studijní zdroje a prameny ...........................34

- 3 (34) -

Úvod

1

Úvod

MATLAB p edstavuje programové prost edí, které poskytuje inženýr m a v dc m snadno použitelné nástroje pro ešení širokého okruhu výpo tových úloh. Obecn se dá íci, že MATLAB je vhodný pro ešení úloh, které svým rozsahem i komplikovaností p esahují možností efektivního ešení v programu EXCEL. Výhody programu MATLAB oproti EXCEL jsou p edevším tyto: •

poskytuje kvalitn jší nástroje pro operace s maticemi,

•

nabízí v tší množství vestav ných funkcí a nástroj ,

•

je univerzáln jší p i použití podmín ných p íkaz a cykl ,

•

v n kterých p ípadech poskytuje v tší možnosti grafické prezentace.

1.1

Cíle

Obsah tohoto modulu je navržen tak, aby ho byl schopen absolvovat kdokoliv bez hlubších znalostí programování. Je založen na využití nápov dy, která je p ímo sou ástí MATLAB. P edkládaný text tvo í pr vodce vybranými stránkami nápov dy, dopl uje ji o p ipomínky a konkrétní p íklady aplikovatelné ve vodohospodá ské praxi. V tomto modulu se nau íte: •

Základy programování v MATLAB.

•

Využívat systému nápov dy MATLAB.

•

Pracovat s maticemi.

•

Importovat a exportovat data.

•

ešit optimaliza ní úlohy.

•

ešit oby ejné diferenciální rovnice a jejich soustavy.

1.2

Požadované znalosti

K úsp šnému zvládnutí tohoto modulu se p edpokládají základní znalosti numerických metod, hydrauliky a pasivní znalost anglického jazyka.

1.3

Doba pot ebná ke studiu

Doba pot ebná ke zvládnutí textu odpovídá 6 výukovým hodinám, tedy cca 6 x 50 min.

1.4

Klí ová slova

Numerické metody, hydraulika, matematický model, numerický model, lineární programování, nelineární programování. - 5 (34) -


1.5

Metodický návod na práci s textem

Text modulu je t eba procházet postupn (nep eskakovat mezi kapitolami) a pr b žn zpracovávat uvedené p íklady a úlohy. Spln ním úlohy je mnohdy podmín no pochopení dalšího textu, který navazuje na poznatky získané b hem jejího ešení. Student by m l p i pro ítání textu sou asn realizovat nazna ené postupy na po íta i v programu MATLAB.

- 6 (34) -

Za ínáme s MATLAB

2

Za ínáme s MATLAB

Po spušt ní aplikace Matlab bude vaše obrazovka odpovídat p ibližn obrázku 2.1. V tuto chvíli pracovní plocha obsahuje t i okna: Command Window, Current Direktory a Command History. Na pracovní ploše si ponechte pouze okno Command Window a ostatní uzav ete pomocí k ížku v pravém horním rohu. V okn Command Window budeme provád t na ítání dat, vkládat p íkazy, funkce, atd. P ed zahájením další práce si vytvo te pracovní adresá u ený pro tento modul a nastavte ho jako pracovní (Current directory) pomocí rozbalovací nabídky v horní ásti obrazovky.

Obr. 2.1- Pracovní plocha aplikace Matlab.

2.1

Základy programování v MATLAB

Program tvo í posloupnost p íkaz , funkcí a prom nných uložených v pam ti po íta e. Programování je procedura sestávající z následujících krok : 1. Definice problému, který má program ešit. 2. Návrh programu (rozhodnutí jaké p íkazy pot ebuji k ešení problému a kde je uložím). 3. Zápis a uložení p íkaz a funkcí. 4. Testování (lad ní) programu. 5. Dokumentace programu (poskytnutí instrukcí pro p ípadné další uživatele programu). - 7 (34) -


Zvládnutí základ programování nep edstavuje nic složitého a vyžaduje pouze pochopení n kolika základních princip a konvencí.

2.1.1

Definice prom nných

Definice prom nných, datové typy a priority používaných operátor p edstavují základní znalosti, které je t eba zvládnout. Do p íkazového okna MATLAB napište:

Tímto p íkazem jste práv p i adili prom nné hodnotu 5. Po íta vytvo il v pam ti prostor, který nazval a uložil do n ho hodnotu 5. V MATLAB jsou všechny prom nné uchovávány ve form matic. je p íkladem nejmenší možné matice s rozm ry 1 x 1. Pokud napíšete:

vytvo íte matici rozm r 2 x 2, nazvanou . Napište:

ímž jste definovali textový et zec, který je uložen matici 1 x 1. Pokud prom nné p i adíte novou hodnotu, pak p vodní je p epsána. Nap íklad:

p epíše p vodní hodnotu 5 na 6. Nyní zadejte:

Tímto znovu p epíšete stávající hodnotu 6 zp t na hodnotu 5. Poznámka B hem naší práce budeme využívat t i základní datové typy prom nných: • reálná ísla – REAL, • celá ísla – INTEGER, • textové et zce – STRING. Pokud v MATLAB definujete novou íselnou prom nnou, pak je jí automaticky p i azen datový typ REAL. Nap íklad:

je uloženo jako reálné íslo na p ibližn 16 desetinných míst. Pokud bychom cht li do této prom nné uložit striktn celé íslo nebo textový et zec, museli bychom to provést nap . takto:

nebo - 8 (34) -

Za ínáme s MATLAB

!"

2.1.2

#$!

Priorita matematických operátor

MATLAB používá konvence b žné z matematiky. Vybrané operátory m žeme set ídit podle priority od nejv tší po nejmenší takto: • ^ (mocnina). • * (násobení), / (d lení). • + (s ítání), - (od ítání).

2.2

Využití nápov dy v MATLAB

Pro vyvolání nápov dy stiskn te klávesu F1 nebo z hlavního menu zvolte položku Help. Na obrázku 2.2 jsou nazna eny t i hlavní cesty, jak je možné se dostat k požadované informaci: • Rejst ík a vyhledávání. • Rozbalovací nabídka obsahující již navštívené stránky nápov dy. • Šipka umož ující návrat o jednu stránku zp t.

Rejst ík a vyhledávání

O jednu stranu zp t

Obr. 2.2- Nápov da aplikace MATLAB.

- 9 (34) -

Rozbalovací nabídka


2.3

Práce s maticemi

Úkol 2.1 V tomto úkolu se samostatn nau íte základy práce s maticemi za pomocí nápov dy. V úvodní nabídce nápov dy, tak jak ji vidíte na obrázku 2.2, zvolte nejprve Getting Started, dále klikn te na Matrices and Arrays a nakonec Matrices and Magic Squares. i te se instrukcemi udávanými nápov dou, prostudujte jednotlivé strany a pokra ujte dokud nedosáhnete strany s nadpisem “The load function”. Tuto stranu již neprocházejte. K pohybu po jednotlivých stánkách nápov dy se výborn hodí tla ítka v pravé horní ásti obrazovky. Informace P i definování prom nných je t eba mít na pam ti, že MATLAB rozlišuje velká a malá písmena, tj. je r zné od . Pokud vložíte p íkaz do p íkazového okna, objeví se vám automaticky odpov . N kdy m že nastat situce, kdy toto zobrazení nepožadujete. V tom p ípad za napsaný p íkaz vložte znak ; a odpov bude uložena, nikoliv však zobrazena. Nápov du k p íkaz m a funkcím je možné vyvolat i p ímo z p íkazového okna zápisem % & '" + hledný p íkaz (nap . % & '" (& ). Po jednotlivých, dosud zadaných, p íkazech v p íkazovém okn je možné se pohybovat pomocí kurzorových šipek a . Dále je možné využít okno s historií zadaných p íkaz (Command History), které jsme uzav eli na za átku cvi ení. Op tovn je možné ho vyvolat z hlavního menu pomocí nabídky Desktop. Pokud jste splnili úkol 2.1 zadaný v této kapitole, ovládáte již tyto operace s maticemi: • Zadání hodnot matice z klávesnice do p íkazového okna. • S ítání ádk , sloupc a diagonál matice. • Výb r zvolené hodnoty z matice a její uložení, jako samostatné prom nné. • Volba sloupce nebo ádku matice a její uložení jako matice nové. • Vytvo ení nulové matice nebo matice vypln né náhodnými ísly.

- 10 (34) -

Za ínáme s MATLAB

P íklad 2.1 Vytvo te matici (12 ádk , 5 sloupc ) vypln nou náhodnými ísly nabývajícími hodnot z uzav eného intervalu od 50 do 100. Použijte p íkaz , data nevkládejte z klávesnice ru n ! P edpokládejme, že vámi vytvo ená data p edstavují syntetickou adu m sí ních srážkových úhrn v letech 2001-2005. Na p íklad takto vytvo ených dat z obrázku 2.3 (vaše hodnoty se samoz ejm mohou lišit) je srážkový úhrn, odpovídající lednu 2005, 58,81 mm.

Obr. 2.3- P íklad vytvo ené matice.

Se tete hodnoty v jednotlivých sloupcích tak, aby jste dostali srážkové úhrny za jednotlivé roky, a tyto hodnoty uložte do samostatné matice ) . Kopii hodnot ze druhého ádku matice uložte do zvláštní matice * , ímž dostanete srážkové úhrny za m síc únor v jednotlivých letech.

2.4

Import dat

Data, která budeme zpracovávat v MATLAB nebudeme obvykle vkládat p ímo do p íkazového okna, ale asto je obdržíme nap . v podob textového nebo binárního souboru. Mohou to být nap . data získaná m ením v terénu, hydraulickým výpo tem, analýzou v GIS apod. Úkol 2.2 Nyní se vra te do nápov dy MATLAB, kterou jste opustili na stran s nadpisem “The load function”. Pokra ujte v úkolech uvedených na této stran nápov dy. K vytvo ení požadovaného textového souboru použijte aplikaci NOTEPAD z prost edí Windows. Informace Pokud používáte p íkaz ' pro na tení dat, MATLAB vyhledává soubor s daty v nastaveném pracovním adresá i (Current Directory v hlavním menu). Pokud chcete na íst soubor z jiného umíst ní, je nutné zadat krom názvu i celou cestu (nap . +,- & $ - &$ -$ # . / ).

- 11 (34) -


Soubor s daty, která chcete na íst p íkazem ' nebo *.dat.

by m l mít koncovku *.txt

Pro výpis soubor v pracovním adresá i zadejte p íkaz

.

Pro vymazání obsahu p íkazového okna lze použít p íkaz +'+. Úkol 2.3 Pokud jste úsp šn splnili p edchozí zadaný úkol, pokra ujte v nápov d od strany s nadpisem “Concatenation” až do strany “Scalar Expansion”. P íklad 2.2 V tomto cvi ení využijete dva vzorové soubory Prutok.txt a Koncentr.txt. Otev ete oba soubory a podívejte se na jejich obsah. Na obrázku 2.4 je soubor Prutok.txt, obsahující seznam p ítok ze zdroj zne išt ní do eky Svratky, sledované b hem t í dn v roce 2005. Jednotlivé zdroje zne išt ní jsou ozna eny kódem. Druhý soubor Konc.txt obsahuje vyhodnocení vzork vody pro jednotlivé zdroje zne išt ní odebraných ve stejných dnech, kdy bylo provád no m ení pr toku. V souboru jsou uvedeny hodnoty BSK.

Obr. 2.4- Obsah souboru Prutok.txt.

Nyní na teme íselná data z obou soubor do MATLAB. Aby byly p i na ítání vynechány ádky s textem, použijeme na za átku znak 0 (viz obrázek 2.5). Tímto zp sobem je ádek ozna en jako komentá .

Obr. 2.5- Soubor Prutok.txt upravený pro na tení do MATLAB.

Import dat provedeme p íkazy: '

1 #

'

3

/2 +/ 2

Tyto p íkazy uloží na tená íselná data do matic nazvaných 1 #

- 12 (34) -

a3

+.

Za ínáme s MATLAB

Úkol 2.4 Vypo ítejte celkové zne išt ní vyjád ené v CHSK, vypoušt né do eky Svratky v jednotkách kg/den, v každém ze t í monitorovacích dn . Výsledek uložte do matice 4 5 6 ) . Limitní množství celkového vypušt ného zne išt ní v CHSK bylo stanoveno hranicí 1500 kg/den. V které dny nebylo toto množství p ekro eno?

2.5

Export dat a uložení MATLAB workspace

Použijte p íkaz 7 % $ k prohlednutí velikosti a typu všech prom nných, které jste dosud vytvo ili. Ke stejnému ú elu lze využít i okno Workspace (aktivace z hlavního menu položkou Desktop). V tuto chvíli chceme ponechat pouze prom nné Prutok, Konc, TDLS. Ostatní vymažeme p íkazem +'& . Nap íklad: +'& Data, která jsme vytvo ili v MATLAB je možné samoz ejm uložit v r zných formátech a využít pro další zpracování. Nap . import do výpo tových program , GIS, CAD atd. Zadejte: $ 8& 4 5 6) /2 4 5 6)

)

99

Tento p íkaz vezme vaši matici TDLS a uloží ji do pracovního adresá e jako textový soubor (ASCII) s názvem TDLS.txt. Otev ete vytvo ený soubor v textovém editoru a zjistíte, že jednotlivé hodnoty jsou na ádcích odd leny mezerami. V p ípad , že chcete použít jíný odd lova , nap . tabulátor: $ 8& 4 5 6) /2 4 5 6) : )

99 4

)

Informace K získání podrobn jších informací o p íkazu $ 8 & m žete použít nap . již zmi ovaný p íkaz % & '" $ 8 & / asto pot ebujete uložit všechny prom nné najednou, pak zadejte: $ 8&

(& 8 $ # .

#

Tímto se vytvo í soubor s názvem (& 8 $ # . # / v MATLAB – z hlavního menu File nebo p íkazem ' prost edí (workspace) pod názvem 8 +& ; : $ 8&

8 +& ;

- 13 (34) -

, který je možné otev ít . Uložte vaše pracovní


2.6

Použití soubor typu M-file

Po íta ový program se obvykle skládá z množství p íkaz . Krom velice jednoduchých program jsou tyto p íkazy zpravidla rozd leny do více soubor . Tento zp sob ukládání zaru uje, že je program snadn ji modifikovatelný a efektivn využívaný, tedy celkov “upravený“. V MATLAB jsou tyto soubory s p íkazy ozna ovány jako M-files. Rozlišujeme u nich dva typy: • script M-files, • function M-files.

2.6.1

Script M-file

Script M-files jsou soubory ádk s textovými et zci, které MATLAB interpretuje jako p íkazy. Nap íklad skript, který vytvo í a uloží stejný typ náhodných dat, který byl popsán v p íkladu 2.1 je ukázán na obrázku 2.6.

Obr. 2.6- P íklad script M-file.

Vytvo te v aplikaci NOTEPAD (Poznámkový blok) nový soubor a uložte ho pod názvem srazky.m do pracovního adresá e. Potom spus te tento skript napsáním jeho názvu do p íkazového okna: $

(

Všimn te si, že všechny vypo ítané hodnoty jsou zobrazeny v p íkazovém okn (pokud nevložíte na konec ádku ). P íklad 2.3 Na obrázku 2.6 je uveden p íklad souboru typu script M-file, který provádí výpo et z úkolu 2.4.

Obr. 2.7- P íklad script M-file.

- 14 (34) -

Za ínáme s MATLAB

Spus te uvedený M-file napsáním jeho názvu do p íkazového okna. Pokud se objeví ervené chybové hlášení v p íkazovém okn , ov te následující skute nosti: • Jsou matice Prutok a CHSK na teny v pracovním prost edí? • Jsou Prutok a CHSK pojmenovány p esn tak, jak ve uvedeno na obrázku 2.7? (Nezapome te, že MATLAB rozlišuje malá a velká písmena). • Mají Prutok a CHSK správný formát (matice o 7 sloupcích a 3 ádcích). • Je M-file s p íkazy uložen v pracovním adresá i. Poznámka Povšimn te si, že v p íkazovém okn se vám po spušt ní souboru M-file z obrázku 2.7 objeví pouze hodnota prom nné CelkZatiz. Je to dáno již zmi ovaným použitím znaku ; na konci ádku. MATLAB používá první dva ádky souboru M-file pro ú ely nápov dy. Zkuste do p íkazového okna zadat p íkaz: % & '" <

(= " )

Informace V datových souborech, souborech typu M-file a p íkazovém okn slouží znak % na za átku ádku k ozna ení komentá e. Tento text je p i zpracování programem MATLAB ignorován a není považován za p íkaz. Komentá e jsou velice d ležité, protože napomáhaní k pochopení obsahu souboru, a to jak autorovi, tak p ípadným dalším uživatel m. Úkol 2.5 Matice Zatiz, kterou jste vytvo ili v rámci p edchozího úkolu obsahuje celkem 21 hodnot. Napište script M-file, který : •

provede konverzi matice Zatiz do sloupcového vektoru (využijte stránku nápov dy s nadpisem Concatenation).

•

Stanovte maximální a minimální hodnotu ze všech 21 hodnot matice Zatiz (použijte funkce max a min).

•

Stanovte po et hodnot, které p ekra ují kritérium 150 kg/den (použijte funkce find a numel).

Informace •

Vytvo te nový M-file použitím položky z hlavního menu File | New | M-File a následn ho uložte pod názvem ZatizStat.m.

•

P edtím než za nete psát do souboru M-file, m žete si jednotlivé p íkazy ov it zadáním do p íkazového okna.

- 15 (34) -


•

2.6.2

Nápov du k jednotlivým funkcím m žete rychle získat zadáním nap / % & '" # & ' do p íkazového okna.

Function M-file

Function M-file se mírn liší od skript , protože vstupní anebo výstupní parametry funkce jsou definovány na prvním ádku souboru M-file. Otev ete pomocí hlavního menu MATLAB volbou File | Open soubor s názvem ZatizVypF.m. Jeho obsah je na obrázku 2.8. Napište do p íkazového okna MATLAB: <

(> & ' <

(

6

'+* 1 #

>

)3

MATLAB zašle vstupní parametry (Prutok, CHSK) do funkce LoadCalcF, provede výpo et a zašle zp t výstupní parametry [Zatiz, CelkZatiz], které se zobrazí v p íkazovém okn .

Obr. 2.8- Function M-file s názvem ZatizVypF.m.

Poznámka Použití funkcí oproti skript m m že být v mnoha ohledech výhodn jší. Nap íklad, pokud chcete provád t nad více sadami vstupních dat stejné operace. Prakticky si to m žeme ukázat na funkci ZatizVypF. P edpokládejme, že krom vstupních dat Prutok, CHSK máme ješt další sadu nazvanou Prutok2, CHSK2. Výsledky chceme uložit do prom nných Zatiz2, CelkZatiz2. Pokud bychom cht li použít skript, bylo by nutné zasáhnout p ímo do souboru ZatizVypF.m a upravit názvy prom nných. Pokud však použijeme funkci, je postup mnohem jednodušší. Sta í pouze zm nit názvy prom nných p i volání funkce. Úkol 2.6 Vytvo te nové prom nné Prutok2, CHSK2, které budou obsahovat dvojnásobné hodnoty z prom nných Prutok, CHSK (jedná se pouze o fiktivní p íklad, dvojnásobné hodnoty jsou zvoleny zcela náhodn , pouze za ú elem vytvo ení nové sady dat). Nyní zavolejte funkci ZatizVypF se zm n nými parametry: <

( > &' <

(

6

'+* 1 #

- 16 (34) -

>

)3

Za ínáme s MATLAB

2.6.3

Rozhodovací struktury a cykly

Sou ástí p evážné v tšiny po íta ových program jsou rozhodovací stuktury a cykly. Nap íklad na vodohospodá ském dispe inku shromaž ovány informace o hodnotách pr tok v m rných profilech na vodních tocích. Jestliže (rozhodovací struktura IF) tyto hodnoty p ekro í ur itou mez, je t eba vyhlásit povodový stupe . Mimo to kontrolní systém pr b žn kontroluje údaje z jednotlivých m rných profil (probíhá smy ka, která cykluje p es jednotlivé sníma e). MATLAB obsahuje rozhodovací struktury typu if, elseif a else a cykly typu while nebo for. Úkol 2.7 S použitím nápov dy se seznamte s výše uvedenými funkcemi. P íklad nápov dy je uveden na obrázku 2.9.

Obr. 2.9- Nápov da k rozhodovací struktu e if, elseif, else.

P íklad 2.4 Senzor_alarm.m je soubor M-file obsahující funkci, která prochází data v matici nazvané Senzor_Data. Tato matice obsahuje hodnoty m ení ze senzor snímajících koncentraci nebezpe ných plyn v ovzduší. Funkce kontroluje, zda není p ekro ena kritická mez udaná hodnotou Krit_Hodnota. Výsledky srovnání jsou ukládány do matice Poplach pomocí ozna ení (V) vysoké nebezpe í,(S) st ední a (N) nízké.

- 17 (34) -


Obr. 2.10- Obsah souboru Senzor_alarm.m.

Spus te funkci Senzor_alarm nap . s následujícími vstupními parametry: 1 ?& @) & ( 3

@

A/ A/ A/

) & ( @5

A/ A/ A/

Následn zavolejte funkci 1 " ' +%

)& ( @ '

) & ( @5

>3

@

> 1 +& @ ) & (

Výsledkem bude: 1B 16

=C )

Informace MATLAB je navržen tak, aby pot eba použití cykl byla co možná nejmenší. Pokud používáte n které další programovací jazyky, jako je nap . Visual Basic, pot ebujete cyly por práci s poli. Naproti tomu MATLAB má již tyto funkce p ímo vestav né. Pokud používáte cykly je t eba dávat pozor, abyste nevytvo ili zadáním nevhodných kritérií pro ukon ení tzv.nekone ný cyklus. MATLAB pak b ží stále dokola v nekone né smy ce. P erušení je v takovém p ípad možné stiskem kláves Ctrl+C.

- 18 (34) -

Za ínáme s MATLAB

Úkol 2.8 V rámci tohoto úkolu si vyzkoušíte použití cykl a rozhodovacích struktur. Dále se seznámíte s možnostmi tvorby graf v MATLAB. K vyjád ení pr toku vody v daném profilu na toku použijeme matematický model vyjád ený rovnicí: Q B (t k ) = aQ A (t k = m ) + bQ B (t k −1 ) ,

(2.1)

kde •

QB(tk) je pr tok v dolním profilu (sm rem po toku) v ase tk,

•

QA(tk-m) je pr tok v horním profilu (sm rem po toku) v ase tk-m,

•

m je po et asových krok za které se zm na pr toku v horním profilu projeví v profilu dolním,

•

a,b jsou parametry, které závisí na charakteru vodního toku.

Tento model m že být využit nap . k p edpov di velikosti pr toku QB níže na toku, pokud známe hodnoty pr toku QA v profilu horním (typickým p íkladem je povod ová situace, která vyžaduje p edpov , kdy a jak velká bude kulminace pr toku v profilu B). Abychom mohli uvedený model použít pot ebujeme znát: •

Hodnoty pr tok QA ve všech asech a po áte ní hodnotu pr toku QB.

•

První p edpov pr toku QB je možná nejd íve za m asových krok od provedení prvního m ení QA.

•

Hodnoty parametr a,b získáme kalibrací.

Postup je následující: 1. Vytvo te matici vstupních parametr a, b. Budeme p edpokládat, že a=1 a b=1. 1 . 2. Zapište funkci M-file nazvanou Predpoved.m, která: –

na te soubor s daty hydrog.dat (viz obrázek 2.11),

–

spustí model pro odhad pr toku QB jako funkce QA pro as t z intervalu 3 až 30 hod, s asovým krokem 1 hod,

–

vykreslí graf pro pr toky QA, QB(m ené), QB(model),

–

zobrazí výstrahu, jestliže pr tok QB(model) p ekro í hodnotu 100 m3/s.

Jediným vstupním parametrem funkce Predpoved bude matice vstupních parametr Parm. 3. Zkoušejte spoušt t funkci Predpoved s r znými vstupními parametry tak, aby se asový pr b h hodnoty QB(model) co nejvíce p iblížil hodnotám QB(m ené). 4. Nezapome te si uložit vaše pracovní prost edí p íkazem >>save ukol2_8. - 19 (34) -


Výsledný graf pro parametry a=1 a b=1ve uveden na obrázku 2.12.

Obr. 2.11- Soubor hydrog.dat.

Obr. 2.12- Výsledný graf pro parametry a=1 a b=1.

- 20 (34) -

ešení optimaliza ních úloh

3

ešení optimaliza ních úloh

MATLAB obsahuje dv základní optimaliza ní funkce nazvané fminbnd (optimalizuje jednu prom nnou) a fminsearch (optimalizuje více prom nných). Navíc MATLAB disponuje vestav ným souborem nástroj , který zahrnuje pokro ilejší algoritmy než zmi ované funkce. V tomto cvi ení použijeme funkci fminsearch k nalezení vstupních parametr a,b, které jsme používali v úkolu 2.8 p i odhadu pr toku QB. Úkol 3.1 Získejte podrobn jší informace o funkci fminsearch s využitím nápov dy MATLAB. Zjist te, jaký uvedená funkce používá algoritmus a popište jeho princip. Z nápov dy v úkolu 3.1 jste zjistili, že syntaxe funkce fminsearch má ve své nejjednodušší podob tento tvar: 2

$&

+% D

# >2A

kde 2 je matice s hodnotami ur enými k optimalizaci, # je název ú elové funkce a 2A je matice s po áte ním odhadem hodnot 2. P i zavolání této funkce hledá MATLAB ešení s maximální odchylkou 0,01%. Základní princip spo ívá v hledání minima ú elové funkce # pro vstupní hodnoty z matice 2. V našem p ípad se p i hledání hodnot parametr a,b budeme snažit minimalizovat odchylku mezi hodnotami QB(model) a QB(m ené). Použijeme principu metody nejmenších tverc a ú elová funkce bude vyjád ena jako suma kvadratických odchylek (QB(m ené)-QB(model))2. Zápis ú elové funkce E * bude mít v syntaxi MATLAB tvar: E * $#

F

,

F

,

/G

Úkol 3.2 K nalezení hodnot parametr a,b využijte ú elovou funkci definovanou v p edchozím odstavci. V této podob ji dopl te na konec souboru M-file Predpoved.m. Doporu ujeme však ješt p edtím provést tyto úpravy: •

Soubor Predpoved.m uložte pod jménem OptAB.m, aby se vám zachovala jeho p vodní verze.

•

V souboru OptAB.m dopl te na konce ádk s p íkazy znak “;“, odstra te vykreslování grafu a vypisování výstrah.

K nalezení parametr a,b spl ujících daná kritéria zadejte p íkaz: 1

$&

+% D B "

>

- 21 (34) -


Jaké jsou hodnoty nalezených parametr a, b? Ov te jejich správnost spušt ním funkce Predpoved s nov nalezenými hodnotami: 1 & " 8&

1

Informace Použitím znaku “;” na konci p íkazových ádk zabráníte MATLAB vypisování výsledk do p íkazového okna v každém itera ním kroku. Toto opat ení vám podstan zmenší celkovou dobu zpracování úlohy. Úkol 3.3 V p edchozím úkolu 3.2 použijte p íkaz optimset ke zm n nastavení optimalizace. Nap íklad: " 1

$

"

$& $&

H

29 & !> A

+% D B "

>

> "

$

Toto nastavení omezí maximální po et iterací na 10. Prozkoumejte další možnosti nastavení optimalizace. M žete doporu it n jaká další vylepšení funkce OptAB?

- 22 (34) -

ešení oby ejných diferenciálních rovnic

4


Postup numerického ešení oby ejných diferenciálních rovnic budeme demonstrovat na p íklad úloh prázdn ní a pln ní nádob, které vychází ze základní objemové rovnice: Qo − Q p dz , =− dt S

(4.1)

kde z(t) zna í polohu hladiny vzhledem k ose výtokového otvoru, Qp(t) je p ítok do nádrže, Qo(z) je odtok z nádrže, S(z) je plocha hladiny a t je as. Uvedená diferenciální rovnice má analytické ešení za p edpokladu Qp=konst. P íklad 4.1 Kónická nádoba výšky h=1,0 m se svislou osou má horní pr m r D1=0,80 m, pr m r dna D2=0,30 m. Ve dn je kruhový otvor o pr m ru Dv=0,30 mm (µv=0,62). Je t eba ur it dobu úplného vyprázdn ní nádoby, je-li napln na až po okraj vodou. P ítok do nádrže Qp=0 m3/s [6]. ešení provedeme jak analyticky, tak numericky. Nejprve je t eba vyjád it funk ní závislost plochy hladiny S=f(z): S=

kde A=

π 4

(D2 + Az )2 ,

D1 − D2 h

(4.2)

.

(4.3)

Dále vyjád íme odtok z nádrže Qo jako: Qo = µ v S v 2 gz

.

(4.4)

Vztahy 4.2 až 4.4 dosadíme do rovnice 4.1 a integrací dostaneme hledaný vztah pro dobu prázdn ní (analytické ešení): t=

2

µ v Dv2

2g

D22 h +

2 1 D2 Ah 3 / 2 + A 2 h 5 / 2 3 5

.

(4.5)

Dosazením p íslušných hodnot ze zadání do vztahu 4.5 obdržíme výsledek analytického ešení t=194,2 s. Numerické ešení oby ejné diferenciální rovnice 4.1 provedeme RungovouKuttovou metodou 4. až 5. ádu. K tomuto ú elu použijeme funkci MATALAB s názvem ode45. Její spušt ní se provádí zadáním p íkazu v následujícím tvaru:

Matice s výsledky obsahující as to až tf

>(

&

D C

Matice s výsledky obsahující hodnoty z v asechto až tf

Funkce dle použité numerické metody

.>

A

I

>

M-file obsahující ešenou oby . dif. rovnici

- 23 (34) -

Hodnoty to a tf

Po áte ní hodnota z v ase to


Obsah souboru M-file s názvem Nadob.m je patrný z obrázku 4.1.

Obr. 4.1- Obsah souboru Nadob.m

Po spušt ní výpo tu obdržíme výsledek numerického ešení t=193,88 s. Úkol 4.1 Prove te analytické ešení p edchozího p íkladu tak, aby jste získali hodnoty pro vynesení grafu funk ní závislosti z=f(t). Do grafu rovn ž dopl te výsledky numerického ešení. Pro lepší p ehlednost slu te získané výsledky do jedné matice Vysl_nadob, která bude obsahovat tyto t i sloupce: as, z(analytické ešení) a z(numerické ešení). Kontrolní otázky ím je zp soben vzniklý rozdíl mezi analytickým a numerickým ešením v úkolu 4.1? Jaký vliv bude mít použití jiné metody pro numerické ešení oby ejné parciální diferenciální rovnice (ov te prakticky v MATLAB)? Úkol 4.2 V této úloze provedete posouzení transforma ního ú inku nádrže p i pr chodu povod ové vlny. K ešení máte k dispozici následující vstupní data: •

ára zatopených ploch a objem pro ešenou nádrž (viz obrázek 4.2).

•

asovou závislost p ítoku do nádrže (viz obrázek 4.4).

•

Technické parametry výpustného za ízení a bezpe nostního p elivu.

ešená nádrž má zemní, sypanou hráz, jejíž sou ástí je objekt spodních výpusti a bo ní, nehrazený, bezpe nostní p eliv (viz obrázek 4.3). V okamžiku p íchodu povod ové vlny ( as t=0 s) je kóta hladiny v nádrži pod úrovni maximálního zásobního prostoru (315,50 m n.m.). Budeme p edpokládat, že od tohoto okamžiku za nou spodní výpusti vypoušt t z nádrže neškodný pr tok 22 m3/s. Tento pr tok z stane konstantní po celou dobu pr chodu povod ové vlny, ehož je dosaženo regulací pomocí uzáv r na spodních výpustech. Od doby, kdy velikost p ítoku do nádrže p esáhne vypoušt né množství za ne hladina v nádrži stoupat. V okamžiku, kdy dosáhne kóty ko-

- 24 (34) -


runy bezpe nostního p elivu (316,40 m n.m.), zapo ne sou asné prázdn ní nádrže spodními výpustmi a bezpe nostním p elivem. Velikost pr toku p epadajícího p es pevnou hranu bo ního p elivu je samoz ejm závislé na aktuální poloze hladiny v nádrži. Na p elivu budeme uvažovat po celou dobu dokonalý p epad. Ú inná délka p elivné hrany je 26 m a hodnotu p epadového sou initele zvolte m=0,51. Pr tok z bezpe nostního p elivu je skluzem odvád n do koryta pod hrází. Stanovte asový pr b h pr tok v koryt pod hrázi od chvíle p íhodu povod ové vlny až po okamžik, kdy hladina v nádrži poklesne zp t na úrove p elivné hrany bo ního p elivu (tj. hladina maximálního zásobního prostoru.) •

Vypo tený asový pr b h pr tok v koryt pod hrází vyneste do stejného grafu jako velikost p ítok do nádrže. Prove te srovnání a zd vodn te vzniklé rozdíly.

•

ešte stejnou úlohu ve variant , kdy hladina v nádrži je v okamžiku p íchodu povod ové vlny práv v úrovni maximální hladiny zásobního prostoru (315,50 m n.m.).

Obr. 4.2- áry zatopených ploch a objem dle [10]

- 25 (34) -


Obr. 4.3- Podélný profil nádrže dle [10]

asový pr b h p ítok do nádrže 120 108 100

Q [m3/s]

80

60 50 40 20

20

15

11 3

0 0

10

20

30

40

t [hod]

Obr. 4.4- asový pr b h p ítok do nádrže

- 26 (34) -

50

60

ešení soustav oby ejných diferenciálních rovnic

5


Postup ešení soustavy oby ejných diferenciálních rovnic si p iblížíme na p íklad ešení vyrovnávací komory na p ivad i vodní elektrárny. Poznámka Totožný p íklad je ešen v modulu 2 pomocí programu EXCEL. Tento postup byl zvolen zám rn , aby bylo možné srovnat zp soby ešení stejné úlohy pomocí odlišných program (EXCEL, MATLAB). Pro p ípad, že by text modulu 3 byl využíván nezávisle, bez modulu 2, uvádíme kompletní zadání v etn teoretického úvodu. Vyrovnávací komory tvo í ást tlakového p ivad e nebo odpadu vodní elektrárny. Ú el vyrovnávacích komor je dvojí: 1. zmírn ní vodního rázu, 2. vytvo ení nádrže, která v prvních okamžicích pozm n pracovního režimu pojme p ebyte né množství vody nebo dodává tlakovému potrubí chyb jící pr tok. Zmírn ní vodního rázu vyrovnávací komorou se projevuje omezením škodlivého p sobení vodního rázu na krátké tlakové potrubí, zatímco dlouhý p ivad z stává prakticky uchrán n. Dále pak zkrácením doby p sobení p ímého rázu v tlakovém potrubí, takže se sníží maximum tlakového p evýšení. Jakákoli zm na pracovního režimu vyvolá v soustav vodní nádrž – tlakový p ivad – vyrovnávací komora – tlakové potrubí – elektrárna neustálený pohyb vody, který se projeví jednak vodním rázem, jednak oscila ním pohybem vody v p ivad i a vyrovnávací komo e. Principem hydraulického ešení vyrovnávacích komor je hledání závislosti zm ny rychlosti proud ní (v) v p ivad i a polohy hladiny (z) ve vyrovnávací komo e na ase (t) pro známou asovou zm nu pr toku Q=f(t). Pro ur ení obou neznámých je nutno sestavit dv diferenciální rovnice, kterým se vzhledem k periodickým výkyv m hladiny ve vyrovnávací komo e íká oscila ní. První rovnice 5.1 je rovnicí kontinuity a vyjad uje rovnost pr toku p ivad em p ed vyrovnávací komorou, p ítokem do vyrovnávací komory a pr tokem turbínou: dz 1 = (Q − Sv) , dt S k

(5.1)

kde t je as, z je okamžitá poloha hladiny ve vyrovnávací komo e (orientace je kladná, když je hladina zaklesnutá pod hydrostatickou hladinou v nádrži a záporná v opa ném p ípad ), Sk zna í plochu hladiny ve vyrovnávací komo e, Q je okamžitá hodnota pr toku od vyrovnávací komory k turbín a v st ední rychlost v p ivad i mezi akumula ní nádrží a vyrovnávací komorou (kladná orientace je ve sm ru proud ní k vyrovnávací komo e).

- 27 (34) -


Druhá rovnice 5.2 (pohybová) je odvozena na základ Newtonových zákon a vyjad uje závislost zrychlení vodní hmoty na poloze hladiny ve vyrovnávací komo e a na velikosti ztrát t ením v p ivad i: dv g = (z dt l

Zt ) ,

(5.2)

kde g zna í gravita ní zrychlení, l délku p ivad e p ed vyrovnávací komorou a Zt souhrn tlakových ztrát v p ivad i mezi akumula ní nádrží a vyrovnávací komorou. Soustavu oby ejných diferenciálních rovnic budeme ešit numericky metodou Rungovou-Kuttovou s použitím MATLAB v rámci následujícího p íkladu. Informace Stru ný teoretický popis Rungeovy-Kuttovy metody je uveden v modulu 2 (kapitola 5). P íklad 5.1 Vy ešte asový pr b h výkyv hladiny ve vyrovnávací komo e válcového tvaru s volnou hladinou p i náhlém uzav ení p ivad e na elektrárnu za komorou (viz obrázek 5.1).

Obr. 5.1- Schéma vyrovnávací komory

Je dán pr tok p ed uzav ením Q=25,0 m3/s, pr m r kruhového p ivad e D=3,57 m, plocha pr ezu komory Sk=100 m2, rychlost v p ivad i p ed uzav ením vo=2,5 m/s, délka p ivad e mezi akumula ní nádrží a vyrovnávací komorou l=3000 m a stupe drsnosti p ivad e n=0,016. ešení zahájíme výpo tem ztrát v p ivad i za pomocí Chezyho rovnice: Zt = l

Q2 C 2S 2R

,

(5.3)

- 28 (34) -


kde C je rychlostní sou initel dla Manninga a R hydraulický polom r. Dosazením z C po úprav dostaneme vztah: Zt = v2

n 2l R

4/3

= v 2ζ

.

(5.4)

Nyní již máme p ipraveny všechny pot ebné vztahy a m žeme p istoupit k ešení v MATLAB. Postup ešení je obdobný jako v p ípad ešení oby ejných diferenciálních rovnic, s tím rozdílem, že nyní máme 2 neznámé veli iny v a z v podob vektoru x. Soubor M-file nazvaný Komora.m je uveden na obrázku 5.2. Spušt ní ešení pomocí funkce ode45 provedeme následujícím p íkazem: 4 ><

&

D 3

>A

AAA > /

/ I

Vektor po áte ních podmínek – rychlost (v) a poloha hladiny v komo e (z) v ase to

Hodnota po áte ní rychlosti v p ivad i (v) v ase to=0 s je dána zadáním. Poloha hladiny (z) ve vyrovnávací komo e v ase to=0 s se stanoví výpo tem ztráty (Zt) dle rovnice 2.10. Hodnota je kladná, protože hladina v komo e je pod úrovní hydrostatické hladiny v nádrži

Obr. 5.2- Soubor Komora.m.

- 29 (34) -


Po skon ení výpo tu provedeme vykreslení výsledk do grafu p íkazem: "'

4 > < ,> >4 >< ,>

Výsledek je uveden na obrázku 5.3. 12 z v

10 8

z [m], v [m/s]

6 4 2 0 -2 -4 -6 -8

0

200

400

600

800

1000 1200 Cas [s]

1400

1600

1800

2000

Obr. 5.3- Pr b h st edních rychlostí proud ní v p ivad i a kolísání hladiny ve vyrovnávací komo e v p ípad náhlého odstavení elektrárny.

Informace Znaménko “-“ je v p íkazu plot uvedeno pouze z d vodu v tší názornosti vykreslených výsledk . Jak již bylo uvedeno v p edchozích odstavcích, hodnota z nabývá kladných hodnot v p ípad , že je hladina v komo e pod úrovní hydrostatické hladiny v akumula ní nádrži. Úkol 5.1 Dopl te výsledky p edcházejícího p íkladu o následující informace: •

Popište slovn chování hladiny ve vyrovnávací komo e a pr b h rychlosti v p ivad i mezi akumula ní nádrží a vyrovnávací komorou. Jinými slovy, prove te interpretaci výsledk .

•

Jakých maximálních a minimálních hodnot dosahuje kóta hladiny v nádrži? P edpokládejte, že hladina v akumula ní nádrži je na kót 455,56 m n.m.

•

Zd vodn te na jaké kót se hladina ustálí.

•

Navrhn te úpravu rozm r vyrovnávací komory tak, hladina p i náhlém úplném uzav ení p ívodu vody na turbínu dosáhla maximáln kóty 463,56 m n.m.

•

Úpravou po áte ních podmínek prove te simulaci stavu, kdy na po átku je elektrárna mimo provoz a dojde k jejímu náhlému spušt ní. Odb r na turbíny iní Q=25 m3/s. Výsledky výpo tu by m ly odpovídat obrázku 5.4. Kontrolou je pro vás, že kóta ustálené hladiny v komo e a rych-

- 30 (34) -


lost proud ní vody v p ivad i by m la být shodná s hodnotami v p vodním zadání (tj. v=2,5 m/s a z=5,6 m). •

Jakým konstruk ním opat ením na vyrovnávací komo e, pop . zm nou manipulace na díle, by bylo možné zmenšit po áte ní náhlý pokles hladiny ve vyrovnávací komo e? Navrženou zm nu ov te výpo tem.

ešení 4 2 0

z [m], v [m/s]

-2 -4 -6 -8 -10 -12

z v

-14 -16

0

200

400

600

800

1000 1200 Cas [s]

1400

1600

1800

2000

Obr. 5.4- Pr b h st edních rychlostí proud ní v p ivad i a kolísání hladiny ve vyrovnávací komo e v p ípad náhlého spušt ní elektrárny.

- 31 (34) -

Záv r

6

Záv r

6.1

Shrnutí

V rámci tohoto modulu jste úsp šn zvládli ešení vybraných úloh z oboru vodního hospodá ství v programu MATLAB. Nyní jste schopni s pomocí toho programu importovat a exportovat data, pracovat s maticemi, ešit diferenciální rovnice, pop . jejich soustavy. Seznámili jste se s funkcemi programu MATLAB pro ešení optimaliza ních úloh. Rovn ž ovládáte použití rozhodovacích struktur a cykl b hem výpo t . Výsledky jste schopni prezentovat ve form graf . U ební text je pouze úvodem do programu MATLAB. Jeho možnosti využití jsou samoz ejm daleko širší. V tuto chvíli jste však schopni si své znalosti z ovládání programu MATLAB rozši ovat dle vašich aktuálních pot eb. Domníváme se, že z ešených p íklad je patrný rozdíl v použití program MATLAB a EXCEL na ešení r zných typ úloh, jak již bylo zmi ováno v úvodu u ebního textu modulu.

6.2

Studijní prameny

6.2.1

Seznam použité literatury

[1]

erná, J., Machalický, M., Vogel, J., Zlatník, . Základy numerické matematiky a programování. SNTL. Praha 1987.

[2]

Škrášek, J., Tichý, Z. Základy aplikované matematiky I. SNTL. Praha 1989.

[3]

Škrášek, J., Tichý, Z. Základy aplikované matematiky II. SNTL. Praha 1986.

[4]

Vitásek, E. Základy teorie numerických metod pro ešení diferenciálních rovnic. Academia. Praha 1994.

[5]

Kolá , V. a kol. Hydraulika. SNTL. Praha 1966.

[6]

Bém, J., Ji ínský, K. Hydraulika v p íkladech. VUT Praha 1975.

[7]

Jandora, J., Uhmannová, H. Základy hydrauliky a hydrologie. VUT Brno 1999.

[8]

Broža, V., Gabriel, P., 1990.

[9]

Stara, V., Veselý, J. Hydraulika – P íklady ke cvi ení. VUT. Brno 1988.

[10]

Pernica, M. a kol. Vodní dílo Slušovice. SZN. Praha 1981.

[11]

Škrášek, J., Tichý, Z. Základy aplikované matematiky III. SNTL. Praha 1990.

[12]

Kiely, G. Enviromental Engineering. McGraw-Hill International 1997.

ihák, F. Využití vodní energie.

- 33 (34) -

VUT. Praha


6.2.2 [13]

6.2.3 [14]

Seznam dopl kové studijní literatury íha, J. a kol. Matematické modelování hydrodynamických a disperzních jev . VUT Brno 1997.

Odkazy na další studijní zdroje a prameny Výrobce programu MATLAB - http://www.mathworks.com/

- 34 (34) -

VYSOKÉ UENÍ TECHNICKÉ V BRN FAKULTA STAVEBNÍ ALEŠ DRÁB HYDROINFORMATIKA I MODUL M03 ÚVOD DO MATLAB

Recommend Documents