7. ODE a SIMULINK Jednou z často používaných aplikací v Matlabu je modelování a simulace dynamických systémů. V zásadě můžeme postupovat buď klasicky inženýrsky (popíšeme systém diferenciálními rovnicemi a ty pak vyřešíme), nebo pomocí grafické nadstavby Matlabu právě pro účely modelování dynamických systémů - SIMULINK.
7.1. Řešení diferenciálních rovnic Nejprve velmi jednoduchý příklad s numerických řešením. Řešme rovnici 5 x′ + 4 x = 11
K tomu je nutné napsat funkci, která vrací derivaci, tedy upravíme rovnici na
x′ =
1 (11 − 4 x ) 5
function y = JednoduchaDR(t,x) % jednoducha ilustracni funkce na ODE % 5x_dot + 4x = 11 % prevedu rovnici kde na leve strane je derivace % x_dot = 1/5 * (11 - 4x) y = [1/5 * (11 - 4*x(1))];
A tuto funkci využijeme při volání řešiče diferenciálních rovnic (obyčejných - ODE) clear all; cas = [0 10]; % v sekundach, odkud kam pocatecni_hodnoty = [10]; % pocatecni hodnota x [t,y] = ode45(@JednoduchaDR,cas,pocatecni_hodnoty); plot(t,y(:,1)); xlabel('cas'); ylabel('x');
Řešičů je v Matlabu celá řada, v příkladu je použitý často využívaný řešič ODE45. Příklad na DR druhého řádu
Představme si parašutistu, který skáče s padákem. Uvažujme nejprve jednoduchý případ, kdy letí kolmo dolů, zajímá nás tedy jen jedna osa. Na parašutistu působí gravitace a odporová síla, která je závislá na rychlosti (prozatím předpokládejme že se padák otevře hned). Zaveďme souřadnicový systém xy, působení gravitace bude proti směru osy y (dolů), působení odporové síly dejme prozatím ve směru osy y. Tedy: my′′ = − mg + cy′2
Tohle je diferenciální rovnice druhého řádu. Matlab umí řešit jen rovnice prvního řádu, nebo soustavy rovnic prvního řádu. Využijeme tedy malý trik, kdy rovnici druhého řádu převedeme na dvě rovnice prvního řádu. function y = Parasutista(t,x) % jednoducha ilustracni funkce na ODE % y_dotdot * m = -Mg + const*v^2 % prevedu na dve rovnice, x_dot je rychlost, x_dotdot je zrychleni % x = (x,xdot) % y = (xdot, xdotdot) g = 10; % m/s^2 m = 80; % kg c = 4/15; % konstanta u rychlosti y1 = x(2); y2 = -g + c*x(2).^2/m; y = [y1;y2];
Tato funkce má jako vstupy čas a vektor x, který má nyní dva prvky, první odpovídá poloze a druhý rychlosti. Výstupy jsou opět dva, jeden odpovídá rychlosti (první derivace) a druhý zrychlení (druhá derivace). Trik je na řádku y1 = x(2);
kdy rychlost (vstup dva) jednoduše vrátím jako první výstup (rychlost). Použití je obdobné jako v prvním případě, jen okrajové podmínky jsou již dvě: clear all; cas = [0 15]; % v sekundach, odkud kam pocatecni_hodnoty = [600 0]; % 600m, pocatecni rychlost je nula [t,y] = ode45(@Parasutista2,cas,pocatecni_hodnoty); plot(t,y(:,1),t,abs(y(:,2).*5)); xlabel('cas'); ylabel('poloha, rychlost * 5'); legend('poloha','rychlost');
A můžeme se podívat jak se mění poloha rychlost parašutisty v čase
600 poloha rychlost
500
poloha, rychlost * 5
400
300
200
100
0
-100
0
5
10
15
cas
Je zde pěkně vidět jak rychlost narůstá, se zvětšující se rychlostí roste i odpor a tak se rychlost limitně blíží k určité hranici. Je zde ale fyzikální problém - parašutista spokojeně pokračuje v pádu i pod zemí. Potřebujeme tedy simulaci zastavit tehdy, když dojde k dosažení určité hodnoty některé z proměnných. K tomu musíme v Matlabu definovat takzvanou událost: function [udalost stop direction] = para_udalost(t,x) % funkce definujici udalost ktera zastavi provadeni simulace udalost(1) = x(1); % udalost (kdy je dana velicina rovna nule) % kdyz bychom chteli aby se simulace zastavila metr pod zemi, bylo by to % udalost = x(1)+1 stop(1) = 1; % zastavit simulaci direction(1) = 0; % podle toho z ktere strany se blizime k nule (nastaveni na 0 znamena ze je to jedno)
A handle na tuto funkci musí být obsažen ve volbách řešiče: clear all; cas = [0 15]; % v sekundach, odkud kam pocatecni_hodnoty = [600 0]; % 600m, pocatecni rychlost je nula options = odeset('Events',@para_udalost); [t,y] = ode45(@Parasutista2,cas,pocatecni_hodnoty,options); plot(t,y(:,1),t,abs(y(:,2).*5)); xlabel('cas'); ylabel('poloha, rychlost * 5'); legend('poloha','rychlost');
Nyní již bude simulace ukončena korektně:
600 poloha rychlost
500
poloha, rychlost * 5
400
300
200
100
0
-100
0
5
10
15
cas
Parametry
Občas je užitečné do funkce, která vrací derivace vložit nějaký parametr, například u našeho parašutisty by to mohla být jeho hmotnost, odpor těla, odpor padáku a doba volného pádu. Uděláme to jednoduše: do funkce popisující derivace přidáme jeden parametr. Do funkce volající rešič také přidáme parametr (který předem naplníme). A pozor, parametr musíme přidat i do funkce události (i když jej nevyužijeme). function y = Parasutista_parametry(t,x,parametry) % jednoducha ilustracni funkce na ODE - vlozeni parametru % parametry budou ve forme struktury: % parametry.hmotnost - hmotnost parasutisty % parametry.volnypad - doba volneho padu ve vterinach % parametry.odportela - odpor samotneho tela % parametry.odporpadaku - odpor padaku g = 10; % m/s^2 m = parametry.hmotnost; % kg c1 = parametry.odportela; % konstanta u rychlosti pri volnem padu c2 = parametry.odporpadaku; % konstanta u rychlosti s padakem y1 = x(2); if (t < parametry.volnypad) y2 = -g + c1*x(2).^2/m; else y2 = -g + c2*x(2).^2/m; end y = [y1;y2];
a funkce volání: % test parasutisty clear all; cas = [0 15]; % v sekundach, odkud kam pocatecni_hodnoty = [600 0]; % 600m, pocatecni rychlost je nula options = odeset('Events',@para_udalost); parametry.hmotnost = 80; parametry.volnypad = 8; parametry.odportela = 1/15; parametry.odporpadaku = 5/15; [t,y] = ode45(@Parasutista_parametry,cas,pocatecni_hodnoty,options,parametry); plot(t,y(:,1),t,abs(y(:,2).*5)); xlabel('cas'); ylabel('poloha, rychlost * 5'); legend('poloha','rychlost');
A je to ☺
600 poloha rychlost
500
poloha, rychlost * 5
400
300
200
100
0
-100
0
2
4
6
8 cas
10
12
14
7.2. Simulink 7.2.1 Co to je a jak to funguje
Simulink je nadstavba Matlabu na modelování, simulaci a analýzu dynamických systémů (systémů proměnných v čase). Pracovat můžeme se systémy lineárními i nelineárními, modelovanými spojitě nebo diskrétně. Model vytváříme pomocí grafického rozhraní, spustit ho lze jak z rozhraní tak přímo z Matlabu (to je výhodné pro spouštění série simulací). Použití se skládá ze dvou základních kroků. Nejprve vytvoříme model. Model popisuje vztahy mezi vstupy, stavy a výstupy systému. Poté spustíme simulaci. Simulace je v zásadě použití numerického řešiče který pomocí numerické integrace počítá aktuální stavy systému. Model se skládá z bloků. Každý blok obsahuje množinu vstupů, stavů a výstupů. Výstup je funkcí vstupů a stavů. Tato funkce závisí na typu bloku. Stav je proměnná která ovlivňuje výstup bloku a její aktuální hodnota je funkcí předchozí hodnoty stavu a vstupu. Blok tudíž musí hodnotu stavu ukládat, a bloky se stavy se označují také jako bloky s pamětí.
u (vstupy)
x (stavy)
y (vstupy)
Některé bloky stavy nemají, například blok Gain (zesílení) pouze vezme vstupní signál a vynásobí ho konstantní hodnotou zesílení kterou dále pošle na výstup. Gain tudíž nemá žádný (vnitřní) stav. Tak je tomu vždy když výstupní hodnota závisí přímo na vstupní. Každý blok obsahuje několik funkcí, které určují vztahy mezi vstupy, stavy a výstupy v závislosti na čase. Vždy je to: • • •
Výstupní funkce (funkce času, vstupu a stavu) Derivace (spojité systémy) – vrací derivaci stavů, závisí na čase, vstupu a stavu Update funkce (diskrétní systémy) – vrací novou hodnotu stavů v závislosti na čase, vstupu a stavu
7.2.2. První jednoduchý pokus
Spustíme simulink a vytvoříme si nový, prozatím prázdný model (v Matlabu na příkazové řádce napíšeme simulink. Otevře se okno Simulink library browser. V něm zadáme v menu File / New / Model). Při tvorbě modelu přenášíme jednotlivé bloky z knihovny do našeho okna, různě je spojujeme, až je model hotov ☺. Nejprve zkusíme jednoduchý model s bloky Simulink/Sources/SinWave, Simulink/Continuous/Integrator a Simulink/Sinks/Scope. Jednotlivé bloky natáhneme z knihovny (držíme levé tlačítko myši). Spojíme je dle obrázku (každý blok má vstupní a výstupní šipku. Mezi nimi natáhneme spojnici levým tlačítkem myši):
Model uložíme. Nyní můžeme spustit simulaci, v menu Simulation/Run. Po ukončení simulace otevřeme Scope ve kterém vidíme výslednou křivku. Při otevření ostatních bloků můžeme měnit jejich parametry. Vyzkoušejte změnu amplitudy a frekvence sinusovky. Parametry simulace (např. délku) můžeme měnit příkazem Simulation/ConfigurationParameters. Zkusme pozorovat jak sinusovku tak výsledek její integrace. Možností je několik, „rozdvojení“ signálu provedeme pravým tlačítkem myši.
Připojení dalšího Scope bloku
Paralelizace signálu blokem Simulink/SignalRouting/Mux 7.2.3. Troška podrobností o blocích a řešičích
Bloky mohou být spojité nebo diskrétní. Spojité bloky odpovídají na změny ve vstupech spojitě, zatímco diskrétní bloky reagují na změny ve vstupech pouze v násobcích fixního časového intervalu zvaného délka vzorku (sample time). Meze jednotlivými vzorky jsou
výstupy diskrétních bloků konstantní. Každý diskrétní blok má parametr sample time, který určuje délku vzorku. Některé bloky mohou pracovat spojitě i diskrétně (např. Gain). Takový blok má implicitní vzorkování, to znamená že pokud je některý ze vstupů spojitý, pracuje i blok spojitě. Pokud jsou všechny vstupy diskrétní, nastaví se taková délka vzorku která odpovídá nejmenšímu ze vstupů. Simulink umožňuje tvorbu subsystémů, kdy každý z nich představuje vlastní schéma. Vnořování subsystémů není omezeno. Je možné vytvářet vlastní bloky, nebo celé jejich knihovny. Můžete to udělat buď graficky pomocí již existujících bloků, nebo přímo programově vytvořením M-souboru který obsahuje příslušné funkce bloku. Takový M-soubor se pak označuje jako S-funkce. Vytváření S-funkcí tvoří samostatnou kapitolu. Matlab obsahuje celou řadu numerických řešičů (nastavení v Simulation / ConfigurationParameters). Řešiče mohou být s fixním nebo proměnným krokem. Řešiče s proměnným krokem mění délku integračního kroku v závislosti na rychlosti změn stavů (při rychlé změně stavu pracují s jemnějším krokem). Řešiče jsou k dispozici jak spojité tak diskrétní, na model obsahující oba stavy je nutné použít spojitý řešič. Během simulace Simulink kontroluje v každém kroku nespojitosti stavových proměnných. Pokud je nespojitost detekována, je určen přesný čas kdy k ní došlo a výpočet je zpřesněn před a po tomto čase. Je tomu tak proto že nespojitost souvisí s podstatnou změnou chování dynamické soustavy, například poloha poskakující kuličky je nespojitá v okamžiku dopadu na podložku. Tato událost musí být proto určena přesně, jinak by výsledkem mohla být změna rychlosti ve vzduchu, nikoliv přímo na podložce. Algebraické smyčky …
1. Demo Bouncing Ball Model - poskakující kulička
Gumová kulička je vymrštěna vzhůru rychlostí 15 m/s ve výšce 10 m. Elastické vlastnosti způsobí že při dopadu je rychlost odrazu snížena o 20 %. Jedinou vnější silou je gravitace: mx&& = −mg && x = −g ´ poč . podm. x&[0] = 15m / s; x[0] = 10m
Model v Simulinku:
Gravitační zrychlení zadáno blokem Konstanta. Rychlost se získá integrací zrychlení, poloha integrací rychlosti. Počáteční podmínka u polohy je vložena přímo do bloku. Integrátor polohy je omezen směrem dolů (poloha nemůže být záporná, kulička nemůže vletět do podlahy). Integrátor rychlosti má externí počáteční podmínku, která je na počátku rovna 15, v průběhu simulace je rovna vždy stavu bloku při jeho resetu (náraz kuličky na zem) vynásobeno elastickou konstantou (snížení rychlosti) a -1 (změna orientace). Stav bloku odpovídá výstupu bloku v momentu resetu (na počáteční podmínku). Externí reset je přiveden z bloku polohy. Příklad používá Scope s více bloky – je možno upravovat (pravé tlačítko, Signal & Scope Manager).
2. Bungee jumping
Chceme zjistit jaké lano použít. Síly působící na skokana: Tíhová síla G = mg , odpor vzduchu O = k1v + k2 v 2 , kde v =
Fe =
kx
pro x > 0
0
pro x < 0
Newtonův 2. poh. zákon:
∑ F = ma
G − O − Fe = ma
&& x=
1 mg − k1 x& − k2 x& 2 − kx ) ( m
dx , síla od lana dt
