Dynamické systémy Filip Trojan 1. července 1999 2 Obsah 1 ÚVOD 5 2 SPOJITÉ DYNAMICKÉ SYSTÉMY 7 2.1 Teorie........................................ 7 2.2 Ekonomické aplikace................................ 9 2.2.1 Walrasův model tržní rovnováhy..................... 9 2.2.2 Model IS-LM................................ 10 2.2.3 Neoklasický model růstu.......................... 12 2.2.4 Lot kův- Volter r ův model lovce a oběti .................. 14 3 DISKRÉTNÍ DYNAMICKÉ SYSTÉMY 15 3.1 Teorie........................................ 15 3.2 Ekonomické aplikace................................ 19 3.2.1 Samuelsonův hospodářský cyklus (1939)................. 19 3.2.2 Hickův příspěvek k teorii obchodních cyklů (1950)........... 22 4 SIMULACE DYNAMICKÝCH SYSTÉMŮ 23 4.1 Teorie........................................ 23 4.2 Příklady....................................... 24 4.2.1 Pindyck-Rubinfeldův třírovnicový model................. 24 5 OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU 29 5.1 Optimální estimator................................ 29 5.2 Odhad stavu dynamického systému........................ 39 5.3 Kalmanův filtr................................... 41 5.4 Smoothing...................................... 45 5.5 Rozšířený Kalmanův filtr ............................. 46 5.6 Odhad nelineárních systémů metodou Monte Carlo............... 48 5.7 Shrnutí ....................................... 51 5.8 Příklady k procvičení................................ 51 6 IDENTIFIKACE DYNAMICKÝCH SYSTÉMŮ 53 6.1 Úvod......................................... 53 6.2 Klasický přístup .................................. 53 6.2.1 Metoda nejmenších čtverců (LS)..................... 54 6.2.2 Metoda maximální věrohodnosti..................... 56 6.3 Bayesův přístup................................... 56 3 4 7 TEORIE CHAOSU 61 A Přehled maticové algebry 63 A.l Blokové matice................................... 63 A.2 Maticový kalkulus ................................. 64 B Pravděpodobnost 65 B.l Základní pojmy................................... 65 B.2 Číselné charakteristiky (momenty)........................ 65 B.3 Normální rozdělení................................. 66 B.4 Podmíněná hustota pravděpodobnosti...................... 66 B.5 Podmíněné číselné charakteristiky......................... 67 B.6 Lineární regrese................................... 67 Kapitola 1 ÚVOD Slovo systém je jedním z nejčastéji používaných pojmů ve všech odvětvích vědy. Používá se pro označení pojmů abstraktních i konkrétních. Teorie systémů, která zaznamenala nebývalý rozvoj zejména v posledních padesáti letech, zkoumá objekty i jevy ve vzájemných souvislostech. Pojem systém je skutečně používán v nejrůznějších souvislostech. V tomto kurzu budeme pod pojmem dynamický systém chápat určitý časově neměnný vztah mezi okamžitými a minulými nebo budoucími hodnotami daných veličin. Přesná definice podle V. V. Němyckého zní takto: Nechť S je separabilní metrický prostor a T C R. Dynamický systém je pak množina transformací $:TxľxS-»S, které splňují: a) Pro všechna čo,či,Č2 S T, x G S platí $(Č2,čo,:c) = $(h,ti, $(ti,to,x)) b) Pro všechna č, to■ t a xn —>■ x platí $(tn,to,xn) —>■ $(t,to,x). Množinu S nazveme stavový prostor a libovolný bod x e S pak nazýváme stavem systému. Obraz $(t,to,x) představuje stav systému v čase č, byl-li systém v čase čo ve stavu x. Dynamické systémy se dělí podle několika kriterií: Podle typu množiny T Diskrétní systémy pro spočetnou množinu T (většinou T = Z) a spojité pro nespočetnou množinu T (většinou T = R). Podle typu stavové veličiny Deterministické systémy, jestliže S je množina reálných vektorů dimenze n a stochastické systémy, jestliže S je množina náhodných vektorů dimenze n. Podle vazby systému na okolí Rozlišujeme systémy neřízené (uzavřené), jestliže systém nemá žádné spojení s okolím, tj. množina transformací má jediný prvek, a systémy řízené (otevřené), kdy transformace F závisí na stavu okolí u G Rm, takzvaný řídící vektor. Podle časové závislosti vlastností systému Máme časově neměnné systémy (časově invariantní nebo též autonomní systémy), které splňují $(č, to,x) = $(č + Ač, čo + Ač, x), a časově proměnlivé systémy, jejichž vlastnosti se s časem mění. 5 6 KAPITOLA 1. ÚVOD Podle funkčních závislostí mezi veličinami Dělíme systémy na lineární, jestliže vnější popis systému je možný pomocí lineárních rovnic a nelineární, kde lineární popis je neadekvátní. Aplikace dynamických systémů je možno nalézt v těchto oborech: • elektrické systémy • mechanické systémy • hydraulické systémy • tepelné systémy • elektromechanické systémy • biologické systémy • ekonomické systémy Každý z uvedených oborů volí k dynamickým systémům poněkud odlišný přístup vzhledem k odlišnosti příslušných vědních oborů. Snad nejvíce jsou propracovány systémy elektrické, výsledky však nejsou vždy přenositelné. V tomto textu se budeme zabývat výhradně systémy ekonomickými, které jsou v mnoha ohledech komplikovanější a abstraktnější než klasické systémy technické. Přesto jejich výzkum zaznamenal v posledních desetiletích bouřlivého rozvoje a jejich aplikace se stále častěji objevují jako nástroje ekonomického rozhodování. Tato skutečnost ukazuje dynamické systémy jako disciplínu, která je hodna studia. Kapitola 2 SPOJITÉ DYNAMICKÉ SYSTÉMY 2.1 Teorie V této kapitole se omezíme pouze na deterministické systémy, takže stavový prostor S = M.n. Spojité dynamické systémy mohou být popsány mnoha různými způsoby: • diferenciální rovnicí • operátorovým přenosem • frekvenčním přenosem • frekvenční charakteristikou • impulsní odezvou • přechodovou charakteristikou • rozložením nul a pólů přenosu Všechny uvedené charakteristiky jsou vzájemně převoditelně (viz. např [?]). Věnujme se nyní representaci pomocí diferenciální rovnice. Dynamický systém je obecně popsán systémem obyčejných diferenciálních rovnic m-tého řádu ve tvaru x{m) = f(x^m-1\...,x,x) (2.1) kde / : Rn ->■ Rn je hladká funkce. Při splnění tzv. Cauchyho počáteční podmínky x(0) = & *(o) = 6 ... (2-2) ^(m-1)(o) = em-i je řešením systému jednoznačně určená funkce x : R —> W1. 7 8 KAPITOLA 2. SPOJITÉ DYNAMICKÉ SYSTÉMY Veta 2.1. Existuje soustava ODR 1. řádu taková, že systém m-tého řádu (2.1) je s touto soustavou ekvivalentní, tj. existuje bijekce mezi množinami všech řešení obou soustav. Důkaz. Položme y = (x,x,... ,x^m~l">) = (yo,yi, ■ ■ ■ ,ym-i)- Zřejmě platí y = (x,x,... ,x^m">). Přitom (2.3) x = 2/1 X = V2 x{m- _1) = ym-i X^ m) = ňvo, yi,---, Počáteční podmínky můžeme přepsat jako 2/0(0) = & 2/i(0) = 6 2/m-l(0) = 4"i-i což můžeme zkráceně napsat jako y = F{y) y(o) = í (2.4) (2.5) což je systém diferenciálních rovnic prvního řádu s Cauchyho počáteční podmínkou. Vztah mezi a; a y je bijekce, takže oba systémy jsou ekvivalentní. D Nadále se tedy budeme zabývat systémy prvního řádu x = f{x) (2.6) kde á = ( %S ..., ijf ) a / je funkce / : Rn -»■ Rn. T HF' • • •' ~3T__ Přiklad 2.2. Uvažujme skalární systém druhého řádu ve tvaru ý = ay + by + c s Cauchyho počáteční podmínkou 2/(0) = Vo 2/(0) = 2/1 Položme a; = (£1,^2) = (y, ý)- Potom platí: x\ = X2 (2.7a) ±2 = ax2 + 6a; 1 + c (2-7b) s počáteční podmínkou a;i(0) = yo »2(0) =2/i. Rovnice (2.7) definují spojitý dynamický systém prvního řádu s Cauchyho počáteční podmínkou. 2.2. EKONOMICKÉ APLIKACE 9 Systém diferenciálních rovnic (2.6) má i svou fyzikální interpretaci. Vektorová funkce / : Rra —> Rra definuje na W1 rychlostní pole. Vyskytne-li se v bodě x částice o jednotkové hmotnosti, udělí se jí rychlost f{x). Trajektorie částice je pak již jednoznačně určena funkcí /• Byl-li systém v čase to ve stavu xq, můžeme definovat zobrazení $t stavového prostoru do sebe takové, že ke každému stavu x(0) a času t přiřadí stav x(t). Toto zobrazení se nazývá operátor toku nebo také evoluční operátor. Toto zobrazení je totožné s transformací $ podle Němyckého a platí ®t(x) =®(t,to,x) Operátor toku je také řešením soustavy (2.6), neboť v každém bodě í je soustava (2.6) splněna a zároveň je možno na základě operátoru toku stanovit ke každému stavu x jeho budoucnost i minulost. Řešení soustavy (2.6) vyhovující počáteční podmínce x(to) = xq lze vyjádřit pomocí evolučního operátoru takto: X{t) = ^t-to(xo) Takovou křivku x (i) nazýváme trajektorie systému. Z definice evolučního operátoru vyplývají následující vztahy: <Š>s+t(x) = $s($t(a;)) (2.8) $t($-t(x)) = <í>o(x)=x (2.9) ^t"1 = $-t (2-10) Definice 2.3. Existuje-li x G S takové, že $t(%) = x pro všechna t, pak bod x nazveme kritickým bodem (pevným bodem) systému. Kritické body hrají při vyšetřování dynamických systémů významnou roli. Definice 2.4. Soustava parametrických křivek $t(%) prc> každé x e S tvoří tzv. fázový portrét systému. Zkoumáním velkého množství dynamických systémů se zjistilo, že tvary trajektorií systémů mohou nabývat množství různých podob. Některé z nich konvergují, monotónně nebo periodicky, k nějaké limitě, což může být bod a nebo množina bodů (kružnice). Mluvíme o stabilním systému. Některé další trajektorie se po vychýlení z rovnováhy pohybují směrem ven, aby se do ní opět časem vrátily (homoklinický orbit) nebo divergují do nevlastního bodu (pak mluvíme o nestabilním systému). Některé trajektorie jsou nejprve přitahovány do pevného bodu aby se, poté co dosáhnou jeho bezprostředního okolí, odchýlily do nekonečna (sedlový pevný bod). Některé trajektorie mají periodický charakter (těm pak říkáme orbity). Některé systémy jsou lokálně stabilní, tj. konvergují k limitě pouze v okolí pevného bodu, jiné jsou stabilní globálně, tj. konvergují z kteréhokoliv bodu. Existují dokonce trajektorie, které se v některém bodě rozdělují do několika větví. Takovému jevu říkáme bifurkace. Jsou známy systémy, které reagují nespojitým katastrofickým skokem na infinitesimální změnu některého parametru. Takové a ještě mnohé neprozkoumané podoby mohou mít dynamické systémy. 2.2 Ekonomické aplikace 2.2.1 Walrasův model tržní rovnováhy Walras ukazuje rovnováhu trhu jako výsledek přizpůsobovacího procesu typu aukce: poté, co licitátor vyhlásí cenu p, kupující vyhlásí poptávku D(j>) a prodávající vyhlásí nabídku S{p). 10 KAPITOLA 2. SPOJITÉ DYNAMICKÉ SYSTÉMY Rozhodující pro další postup je přebytek poptávky E (p) = D{p) — S (p). Je-li E (p) > 0, proběhne nové kolo při vyšší ceně. Je-li E (p) < 0, proběhne nové kolo při nižší ceně. Po dosažení rovnováhy E (p) = 0 se provede směna. Přizpůsobovací proces můžeme modelovat rovnicí p = k E (p) kde k je kladná konstanta, odrážející rychlost přizpůsobování. V našem příkladě uvažujeme lineární funkce poptávky D(p) = a + ßp a nabídky S(p) = 7 + öp. Počátek přizpůsobovacího procesu je p(0) = po- Upravíme stavovou rovnici: p = k (a + ßp — 7 — öp) = k{ß — ö)p + k (a — 7) Stavová rovnice je lineární diferenciální rovnice prvního řádu s konstantními koeficienty. Řešení hledáme ve tvaru p(ť) = a + bect Dosazením do stavové rovnice a počáteční podmínky vychází »= »-■£! <2-12> c = k(ß-ö) (2.13) takže řešením je křivka ^ÍB+^-iBr3-1" Systém je stabilní při splnění podmínky ß — ö < 0. Ale ß je sklon křivky nabídky a ö je sklon křivky poptávky. Je-li nabídka rostoucí a poptávka klesající, je podmínka splněna a trh je stabilní a konverguje k rovnovážné ceně 7-a p{cc)=pe = — při které je exces poptávky nulový, tzn. D(p(oo)) = S(p(oo)). 2.2.2 Model IS-LM Uvažujme jednoduchou uzavřenou ekonomiku, která se chová podle následujících pravidel: (i) Produkt Y roste s převisem agregátní poptávky D, tj. Y = h(D — S) (ii) Úroková míra r roste s převisem poptávky po penězích L, tj. ř = m(L — M), kde M je nabídka peněz určená centrální bankou. (iii) Poptávka po penězích je proporcionální produktu, tj. L = kY. Uvažuje se pouze transakční poptávka v duchu Fisherovy kvantitativní teorie peněz. (iv) Spotřeba C je proporcionální produktu, tj. C = cY. 2.2. EKONOMICKÉ APLIKACE 11 (v) Investice / klesají s rostoucím úrokem, tj. / = —ar. (vi) Agregátní poptávka je součtem spotřeby a investic, tj. D = C + I. (vii) Agregátní nabídka S je národní produkt, tj. S = Y. (viii) 0 < h,m,a, 0 < c < 1. Položíme-li mezní sklon k úsporám s = 1 — c, můžeme psát Ý = h(D -Y) = h(C + I -Y) = h{c - 1)Y - ahr = -hsY - ahr ř = mkY — rriM Derivujeme-li první rovnici a do vzniklého výrazu dosadíme rovnici druhou, dostaneme Y = —hsY — ahr = —hsY — ahmkY + ahmM neboli Y + hsY + ahmkY = ahmM což je lineární ODR druhého řádu. Napišme její charakteristickou rovnici A2 + hs\ + ahmk = 0 Zde se řešení rozpadá na tři části a) hs2 > Aamk V tomto případě existují 2 reálná řešení charakteristické rovnice - í — hs + \/h?s2 — iahmk j Al = 2 *» = 2 - í — hs — \Jh2s2 — iahmk j kterým odpovídá řešení DR ve tvaru Y(t) = ^- + CieXlt + C2eA2Í Toto řešení bude stabilní, jestliže bude platit Ai < 0, A2 < 0. Ukážeme, že tato podmínka je splněna Vzhledem k podmínce (viii) platí ahmk > 0. Pak také -Aahmk < 0 4^ h2s2 - Aahmk < h2s2 4^ \/h2s2 - Aahmk < \hs\ 4$ \/h2s2 - Aahmk < hs 4$ - (—hs + \/h2s2 - Aahmk) < 0 4$ Ai < 0 a podobně -\/h2s2 - Aahmk < hs 4$ - (—hs - \/h2s2 - Aahmkj < 0 44- A2 < 0 Z čehož vyplývá, že systém je stabilní a je přitahován k pevnému bodu M/k. 12 KAPITOLA 2. SPOJITÉ DYNAMICKÉ SYSTÉMY b) hs2 = 4amk Za této situace má charakteristická rovnice dvojnásobný reálný kořen 1 kterému odpovídá řešení Ai = A2 = A = --hs < 0 Y(t) = ^ + dext + C2text A2 = 2 Toto řešení je stabilní a konverguje k pevnému bodu M/k. c) hs2 < 4amk V tomto případě existují 2 komplexní kořeny charakteristické rovnice Ai = - (—hs + iyAahmk — h2s2) - (—hs — iyAahmk — h2s2 J Položíme-li a = — \hs & ß = \\l4ahmk — h2s2, pak můžeme řešení napsat ve tvaru Y(t) = ^- + eat(Ci cos ßt + C2 sin ßt) (2.14) Protože a < 0, je řešení stabilní a tlumenými sinusovými kmity se přibližuje k pevnému bodu M/k. Při vysokém sklonu k úsporám s je toto přibližování rychlejší, stejně tak jako při vysoké hodnotě koeficientu h, který zachycuje schopnost ekonomiky reagovat na zvýšenou poptávku růstem produktu. Periodické výkyvy produktu jsou obrazem hospodářských cyklů s periodou An/^Aahmk — h2s2. Reálná existence hospodářských cyklů v ekonomickém životě napovídá, že nejčastěji se vyskytuje případ c). Trajektorie skutečné ekonomiky však nikdy neodpovídá ideální křivce (2.14) vzhledem ke značně zjednodušeným předpokladům modelu, zejména izolaci od okolí a konstantní nabídce peněz. 2.2.3 Neoklasický model růstu Tento model, vypracovaný Swanem a Solowem (1956) je založen na následujících předpokladech: (i) Pracovní síla L roste konstantním tempem n, tj. L/L = n. (ii) Úspory jsou lineárně závislé na produktu, tj. S = sY. (iii) Všechny úspory jsou investovány do kapitálu K. (iv) Investice I = K + öK. (v) Produkce probíhá za podmínky konstantních výnosů z rozsahu, tj. Y = F(K,L) = LF(K/L, 1) = Lf(k), kde k = K/L je vybavenost práce kapitálem. (vi) M e (0,1). 2.2. EKONOMICKÉ APLIKACE 13 Spočtěme k: • KL-KĹ _ K K Ĺ I? ~Y~Yl z toho vyplývá, že k _ L -K Ĺ k~K ~K~L dosazením podmínky (i) dostáváme k K k = K-n z podmínek (ii) a (iii) vyplývá, že K = sY -SK dosazením dostáváme I = S| - (ó + n) = s^f{k) - (ó + n) = S-f(k) - (ó + n) Vynásobením k máme k = s f (k) -(6 + n)k = s f (k) - \k kde / je rostoucí konkávni funkce, splňující lanadovy podmínky lim^o f'(k) = oo a lim^oo f (k) = 0. V prípade Cobb-Douglasovy produkční funkce Y = KaLl~a dostáváme Y = KaL-a = ka = f (k) Li a diferenciální rovnice růstu je ve tvaru k = -Xk + ska což je Bernoulliho rovnice. Provedeme substitucí x = kl~a. Pak ka x = (1 — a)k~ak ==> k = x------- 1 — a k — hl~aha — rka Dosazením do DR máme x-------= -\xka + ska 1 — a odkud dostáváme nakonec x = —(1 — et)\x + (1 — a) s což je rovnice lineární. Jejím řešením je trajektorie x(t) = (x0 - j) e-^xt + i nebo ekvivalentně trajektorie kl-a = ^1-a _ j^ e-(l-a)Xt + * protože platí 1 — a> 0,X = ô + n > 0, je model stabilní s pevným bodem k* = i s \ 1~a ô + n 14 KAPITOLA 2. SPOJITÉ DYNAMICKÉ SYSTÉMY 2.2.4 Lotkův-Volterrův model lovce a oběti Tento model má svůj původ v letech po 1. světové válce. Lotka v roce 1925 a Volterra v letech 1931 a 1937 pozorovali výskyty dravých ryb (predators) a jejich kořistí (preys) v Jaderském Moři. Na základě tohoto pozorování zformulovali model, který vysvětluje fluktuace v počtech rybích populací. Tento model je vynikající ilustrací nelineárního dynamického systému v rovině. Populace kořistí x (i) roste přirozeným tempem a (tj. x = ax). Jejich počet je snižován o cxy díky přítomnosti dravců y (i), kteří loví v hejnech kořistí. Současně s tím přirozená populace dravců klesá tempem b (tj.y = —by), ale za přítomnosti obětí se tento pokles snižuje o dxy. Za předpokladu a,b,c,d> 0 pak máme systém x = ax — cxy (2-15) ý = dxy —by (2-16) Položíme-li x = 0, ý = 0, obdržíme dva kritické body systému, z = (0, 0) je tzv. rovnováha po vyhynutí (extinction equilibrium) a z = (|, ^) je tzv. rovnováha koexistence (coexistence equilibrium). Numerické řešení DR (2.15) je v Matlabu ve funkci odedemo. Kapitola 3 DISKRÉTNI DYNAMICKÉ SYSTÉMY 3.1 Teorie Jako diskrétní označujeme ty systémy, jejichž chování je definováno na diskrétní časové množině T = {kT :fceZ}. Proměnná k se nazývá krok. Délka kroku T bývá též označována jako vzorkovací interval nebo perioda vzorkování. Ze spojité funkce /(í) lze získat diskrétní funkci g{k) vzorkováním. Přitom platí /(í) = f(kT) = g(k). Pro vnější popis diskrétních systémů se běžně používá diferenční rovnice a operátorový přenos. Diferenční rovnice, popisující neřízený deterministický systém diskrétního typu má podobu y(k) = f(y(k-l),...,y(k-ly)) (3.1) kde y e W je výstup systému. Tato rovnice ukazuje závislost výstupu y (k) na minulých výstupech y (k — i). Rád systému ly pak udává paměť systému. Přidáme-li do systému prvek řízení, tj. řídící vektor u G Mm, dostaneme rovnici y{k) = f{y{k -l),...,y{k- ly),u(k -l),...,u(k- Q) (3.2) Takto napsaná rovnice naznačuje, že řízení u(k — í) přímo ovlivňuje výstup y{k). Účinek řízení je tedy zpožděn o jeden krok. Takto se chová převážná většina reálných systémů. Někdy je dokonce zpoždění větší. Pak se v rovnici (3.2) objevují členy až od u(k — nu). Sledujeme-li pouze závislost výstupu systému na minulých výstupech a vstupech, říkáme tomu vnější popis. Zavedeme-li stav systému jako veličinu, která plně určuje vývoj systému do budoucna, pak se systém dá napsat jako systém prvního řádu. Mluvíme pak o vnitřním nebo také stavovém popisu systému. Přesněji stav x(k) a řízení u(k) určují jednoznačně stavy x{k + l),x(k + 2),... a výstupy y(k),y(k + l),y(k + 2),.... Definice 3.1. Dynamický systém popsaný rovnicemi x(t + í) = f(x(t),u(t)) (3.3a) y(t) = g(x(t),u(t)) (3.3b) nazveme deterministickým dynamickým systémem ve stavovém tvaru 15 16 KAPITOLA 3. DISKRÉTNÍ DYNAMICKÉ SYSTÉMY Rovnici (3.3a) nazýváme stavovou rovnici (rovnice přechodu, state equation, plant equation, transition equation) a rovnice (3.3b) je zvána rovnicí výstupu (rovnice měření, observation equation, output equation, measurement equation). Systémy ve stavovém tvaru hrají v teorii dynamických systémů významnou roli, protože právě pro stavové systémy jsou vypracovány algoritmy optimální estimace, identifikace a optimálního řízení. Proto je vhodné zabývat se otázkou existence a jednoznačnosti stavového popisu systému, který je zadán svým vnějším popisem (3.2). Stavový popis vždy existuje, není však obecně jednoznačně určen. Jeden z možných převodů popisuje následující věta. Veta 3.2. Položíme-li x(k) = (y(k)T, ...,y(k-ly + 1)T, u(k - if, ...,u(k-lu + iff pak systém (3.2) lze napsat jako systém ve stavovém tvaru. Důkaz. Původní rovnici y(k) = f(y(k -l),...,y(k- ly),u(k -l),...,u(k- lu)) lze substitucí t + 1 = k napsat ve tvaru y{t +1) = f{y{t), ■ ■ ■, y{k - ly +1),u{k - i),...,u{k - lu +1)) = f{x{t)) Zřejmě platí x(t + 1) = (y(t + 1)T, ...,y(t-ly + 2)T, u(tf, ...,u(t-lu + 2)Tf takže první složka vektoru x(t + 1) už svou rovnici má. Najdeme stavové rovnice pro všechny zbývající složky. V těchto rovnicích mohou na pravých stranách vystupovat pouze složky vektorů x{ť) a u{ť). y(t + l) = f(x(t)) y{t) = AlX{t) Ayx(t) u(t) Bix(t) u(t -lu + l) = Blu_lX{ť) kde A\,... ,Aiy,B\,... ,Biu_i jsou permutační matice, vybírající z vektoru x{ť) vždy jen příslušnou část. Celkem lze tedy soustavu napsat jako x{t + l) = F{x{t),u{t)) y(t -ly + 2) u(t) u{t -1) což je popis systému ve stavovém tvaru D 3.1. TEORIE 17 Přiklad 3.3. Uvažujme lineární dynamický systém ve tvaru yt = aiyt-i + a2yt-2 + hut-i + b2ut-2 Posuneme čas o jeden krok dopředu a máme yt+i = am + a2yt-i + Mí + b2ut-\ Nyní položíme xt = (xit,x2t,x3t)T = (yt,yt-i,ut-i)T. Potom platí xt+i = (yt+i,yt,ut)T a můžeme psát xi,t+i = a\X\f, + a2x2tt + 62^3,* + hut x2>t+i = Xitt X3,t+1 = Ut yt = xi,t což je systém ve stavovém tvaru (3.3). Definice 3.4. Dynamický systém popsaný rovnicemi x(t + í) = f(x(t),u(t),v(ť)) (3.4a) y(t) = g(x(t),u(t),w(t)) (3.4b) kde v (ť) G Cnv, w(ť) G Cnw jsou náhodné vektory, splňující podmínku Ev(t) = 0,Ety(í) = 0, nazveme stochastickým dynamickým systémem ve stavovém tvaru. Náhodný vektor v(t) nazýváme šum procesu (process noise) a náhodný vektor w(t) nazýváme šum výstupu (output noise). V naprosté většině případů mají rovnice (3.4a) a (3.4b) speciální tvar x(t + l) = f(x(t),u(t))+v(t) y(t)=g(x(t),u(t)) + w(í) odtud je motivován požadavek na centrovanost obou šumů. Přechod od deterministického ke stochastickému systému je výrazným kvalitativním skokem, poněvadž původní stavový prostor reálných vektorů Rnx přešel ve stavový prostor Cnx náhodných vektorů. Stejně tak výstupy y(ť) jsou nyní náhodnými vektory. Odtud plyne, že ke zkoumání takových systémů je nutno použít teorii pravděpodobnosti. Definice 3.5. Nechť A G W^xn^,B G RniX™«,Ce Wh,xn*,D G Rnyxn^. Soustava x{t + 1) = Ax(t) + Bu{t) + v{t) (3.5) y(t) = Cx(t)+Du(t)+w(t) (3.6) popisuje lineární stochastický dynamický systém ve stavovém tvaru. Šumy v(t),w(ť) musí splňovat následující podmínky 18 KAPITOLA 3. DISKRÉTNÍ DYNAMICKÉ SYSTÉMY (i) Šumy vt a Wt jsou normální s nulovou střední hodnotou. (ii) Evtwf = 0 (iii) EvtvJ = Q (iv) Ewtwf = R Zabývejme se nyní stabilitou lineárního stochastického dynamického systému. Obecně je dynamický systém stabilní tehdy, když je pevný bod systému atraktorem neboli přitahujícím pevným bodem. Z libovolného počátečního stavu potom stav systému konverguje k tomuto pevnému bodu. Stochastický systém je stabilní právě tehdy, když k němu příslušný deterministický systém je stabilní. Pracujeme-li navíc se systémem řízeným, pak o konvergenci stavů je možno hovořit pouze tehdy, je-li vstup systému konstantní. Ve svém důsledku pak vlastně zkoumáme stabilitu deterministického neřízeného systému. Pevný bod lineárního systému s konstantním vstupem u určíme z rovnice x = Ax + Bu odkud dostáváme xe = (I-A)-lBu (3.7) Veta 3.6. Stochastický systém (3.5) je stabilní právě tehdy, když všechna vlastní čísla matice A leží v jednotkovém kruhu. Důkaz. Systém (3.5) je stabilní, právě když je stabilní systém x{t + 1) = Ax{t) + Bu (3.8) kde u je známý konstantní vstup systému. Pevný bod tohoto systému je (3.7). Zavedeme-li substituci y(ť) = x (i) — xe, pak stavová rovnice přejde v x(t + 1) — xe = Ax(ť) + Bu — xe <^ x(t + 1) — xe = Ax{ť) — Axe + Bu + Axe — xe <^ x(t + 1) - xe = A(x(ť) - xe) + Bu — (I - A)xe -^ y(t + 1) = Ay(t) +Bu-(I- A)(I - A)~lBu <& y(t + 1) = Ay(ť) Má-li x(t) konvergovat k xe, pak musí y{ť) konvergovat k nule. Vyjdeme-li z počáteční podmínky y(0) = j/o, pak stav y(t) se dá explicitně spočítat jako y(t) = Ao Je vidět, že má-li y{n) pro libovolné j/o konvergovat k nule, musí An konvergovat k nulové matici. Pro každou čtvercovou matici A existuje regulární matice P taková, že P~lAP = D kde D je diagonální matice, která má v diagonále vlastní čísla matice A D = dvág(Xi,..., XnJ 3.2. EKONOMICKÉ APLIKACE 19 Sloupce matice P pak jsou vlastní vektory matice A. Tato dekompozice je známá jako Schurova dekompozice matice A. Matici A pak můžeme napsat jako A = PDP'1, matici A2 zapíšeme jako A2 = PDP-lPDP~l = PD2P~l a tak dále až matici An zapíšeme jako An = pDnp-l Matice Dn je opět diagonální a obsahuje prvky Dra = diag(A?,...,A™J Posloupnost An konverguje k nule právě tehdy, když posloupnost Dn konverguje k nule a to je právě tehdy, když všechny posloupnosti A™ konvergují k nule. To je splněno, když pro všechna i e {1,..., nx} platí |Ai| < 1 (3.9) Protože vlastní čísla A» jsou obecně komplexní, neznamená (3.9) nic jiného, než že všechny vlastní čísla se v komplexní rovině nachází uvnitř jednotkového kruhu. D Chování stabilního stochastického systému je významně ovlivněno vlastnostmi stochastické složky. Bude-li šum slabý, bude systém po dostatečně dlouhém čase skutečně vykazovat stabilitu a bude jen nepatrně kolísat kolem pevného bodu. Naopak stabilní systém se silným šumem se může chovat naprosto chaoticky. Čím budou vlastní čísla blíže k nule, tím rychleji systém konverguje k pevnému bodu. Při hodnotách blízkých jedné už se konvergence značně zpomaluje a po překročení jednotkového kruhu se z přitahujícího pevného bodu stane odpuzující pevný bod a systém se stane nestabilním. 3.2 Ekonomické aplikace 3.2.1 Samuelsonuv hospodářský cyklus (1939) Samuelsonuv model je postaven na multiplikátoru a akcelerátoru. Struktura modeluje popsána těmito rovnicemi: Ct = Yt-x It = v(Ct - Ct-i) Gt = l Yt = Ct + It + G t Předpokládá se 0 < c < l,v > 0. 1/1 — c je multiplikátor a v je akcelerátor. Substitucí dostaneme typickou diferenční rovnici 2. řádu Yt-c{l+v)Yt_l + cvYt_2 = l Jednoduchým cvičením se dá dokázat, že platí: 20 KAPITOLA 3. DISKRÉTNÍ DYNAMICKÉ SYSTÉMY (a) Každý skalární násobek řešení homogenní rovnice je řešením homogenní rovnice. (b) Součet každých dvou řešení homogenní rovnice je řešením homogenní rovnice. (c) Součet libovolného řešení nehomogenní rovnice s libovolným řešením homogenní rovnice je řešením nehomogenní rovnice. Z uvedeného vyplývá, že množina všech řešení homogenní rovnice tvoří vektorový podprostor a množina všech řešení nehomogenní rovnice tvoří afinní podprostor v prostoru funkcí. Stačí tedy najít jedno řešení j/o nehomogenní rovnice (tzv. partikulární řešeni) a bázi yi,... ,yn všech řešení homegenní rovnice. Libovolné řešení nehomogenní rovnice pak dostaneme jako Lehce se ukáže, že funkce y = yo + aiyi + ... + anyn I/o (í) 1-c je partikulárním řešením naší diferenční rovnice. Hledejme nyní řešení homegenní rovnice ve tvaru Yt = A*. Pak musí platit Aí-c(1+v)A*-1 + ctAí-2 = 0 Vydělením rovnice výrazem A*-2 dostaneme X2 -c(l + v)X + cv = 0 Z této rovnice dostaneme A jako kořen kvadratické rovnice. Řešení závisí na znaménku diskriminantu A = c2(l + v) Acv a) c> Av/(l + v)2 V tomto případě existují 2 reálná řešení, která tvoří bázi řešení homogenní rovnice. A2 = kterým odpovídá obecné řešení ve tvaru Ai = - (c(l + v) + ^c2(l + v)2 - 4cv) - (c(l + v) - ^c2(l + v)2-W) Yt = 1 1-c + ai\\ + a2\\ Toto řešení bude stabilní, jestliže bude platit |Ai| < 1, | A21 < 1- |Ai| < 1 A |A2| < 1 4^ IAiA21 1 -Acv < 1 4^ CV < 1 3.2. EKONOMICKÉ APLIKACE 21 b) c = 4v/(l+v)2 Za této situace má kvadratická rovnice dvojnásobný reálný kořen Al = A2 = A = £(i±í) kterému odpovídá řešení Yt =------+aiA* + a2t\t 1 — c Toto řešení je stabilní za podmínky c(l + v) < 2 <^ v < 1. c) c < 4v/(l +v)2 V tomto případě existují 2 komplexní kořeny charakteristické rovnice Ai = A2 = - (c(l + v) + i^jAcv -c2(l + v)2) - (c(l + v) - ía/4ct-c2(1+v)2) Báze řešení homogenní rovnice je tvořena funkcemi u = \\ a v = A|, což jsou komplexní funkce. Ukážeme, že lze najít reálné funkce, které tvoří bázi řešení. Zaveďme goniometrické vyjádření komplexních kořenů Ai = r(cos9 + i s'm9) A2 = r(cos9 — is'm6) Pak platí ai\\ +a2\l = = air*(cos 6t + i sin 9ť) + e^r^cos 6t — i sin 9ť) = = r*((ai + CI2) cos oí + i{a\ — 02) sin 9ť) = = rt(bicos9t + b2sm9t) přitom r je absolutní hodnota společná komplexním číslům Ai, A2, která je rovna Systém bude stabilní za podmínky r < 1 <ř» cv < 1. Samuelsonuv ekonomický systém může tedy vykazovat 9 kvalitativně různých druhů chování, a to v závislosti na aktuální hodnotě dvou parametrů -vac. Každému typu chování přináleží v parametrickém prostoru (f,c) jedna množina. Těchto 9 množin pokrývá celý parametrický prostor, který jsme si hned v úvodu omezili na množinu (0, 00) x (0,1). 1. c > 4i)/(l + v)2 A cv < 1: monotónní konvergence - stabilní 2. c > 4v/(l + v)2 A cv = 1: neutrální stabilita 3. c > 4i)/(l + v)2 A cv > 1: monotónní divergence - nestabilní 4. c = 4i)/(l + v)2 A v < 1: jeden kriticky tlumený kmit - stabilní 22 KAPITOLA 3. DISKRÉTNÍ DYNAMICKÉ SYSTÉMY 5. c = 4t)/(l + v)2 A v = 1: lineární divergence - nestabilní 6. c = 4t)/(l + v)2 A f > 1: monotónní divergence - nestabilní 7. c < 4t»/(l + v)2 A cv < 1: periodická konvergence - stabilní 8. c < 4t»/(l + v)2 Acv = 1: konstantní periodická oscilace - stabilní 9. c < 4t»/(l + v)2 A cv > 1: periodická divergence - nestabilní 3.2.2 Hickův příspěvek k teorii obchodních cyklů (1950) K ilustraci předešlé věty použijeme Hicksův model hospodářského cyklu z roku 1950. Tento model vysvětluje hospodářský cyklus pomocí interakce mezi multiplikátorem a akcelerátorem. Model má následující strukturu: Ct = {l-s)Yt_l It = A0{l + g)t + v{Yt_l-Yt_2) Yt = Ct + It Spotřeba Ct je zpožděná lineární funkce Yt-\. Parametr s splňující podmínku 0 < s < 1 je mezní sklon k úsporám. Investice It mají dvě komponenty: autonomní investice Ao(l + g)ť rostoucí exponenciálně konstantním tempem g a indukované investice v{Yt-i — Yt-2), které jsou proporcionální růstu GNP v minulém období. Konstanta v > 0 je akcelerátor. Substitucí dostaneme rovnici Yt = (1 - s + v)Yt-1 + vYt_2 + A0(l + gf Autonomní investice můžeme chápat jako exogénni veličinu, kterou zde pro jednoduchost ztotožníme s řízením ut, tj. Yt = (l-s + v)Yt-i + vYt_2 + ut-i Posuneme čas o jeden krok dopředu a máme Yt+i = (1 - s + v)Yt + vYt_x + ut Nyní položíme xt = (xu,X2t) = (Yt,Yt-i) . Potom platí xt+i = (Yt+1,Yt)T odtud pak Xl,t+1 = (1 - s + v)xitt + vx2,t + ut {6.10) X2,t+1 = Xítt a také Yt = xht (3.11) rovnice (3.10) je stavovou rovnicí a rovnice (3.11) je výstupní rovnicí uvažovaného dynamického systému. Kapitola 4 SIMULACE DYNAMICKÝCH SYSTÉMŮ 4.1 Teorie Pojem simulace se úzce váže na použití výpočetní techniky jakožto nástroje pro její realizaci a výstup výsledků. Simulace dynamických systémů se významně liší podle toho, jestli se jedná o systém spojitý nebo diskrétní. Definice 4.1. Buď dán spojitý systém ve stavovém tvaru (2.6) s počáteční podmínkou. Dále nechť T = {tk : k £ {1,..., N}, N e N} je tzv. časový horizont. Simulovat daný systém na horizontu T pak znamená nalézt vzorky trajektorie systému v časech t k- Jinými slovy hledáme posloupnost {xo, ■ ■ ■ ,xn} takovou, že Xk = x(tk), kde funkce x (i) je řešením soustavy DR (2.6). Definice 4.2. Buď dán řízený stochastický diskrétní systém (3.3) s počáteční podmínkou. Dále nechť T = {1,..., N}, WeNje časový horizont. Buď {«(0),..., u(N)} známá posloupnost řídících vektorů—trajektorie řízeni a {v(0),... ,v(N)},{w(0),... ,w(N)} posloupnosti realizací náhodných šumů v,w. Simulovat systém (3.3) na horizontu T pak znamená nalézt trajektorii {x(0),... ,x(N)} stavů x a trajektorii {y(0),... ,y(N)} výstupů y tak, aby byly splněny rovnice (3.3). Pro simulaci spojitých systémů je potřeba numericky vyřešit soustavu diferenciálních rovnic a řešení vyčíslit v bodech tk- K tomu se využívá celá řada známých numerických metod, které lze najít např. v přehledu [?]. V systému Matlab jsou tyto procedury štandartne k dispozici ve funkcích ode23, případně ode45. Simulace diskrétních systémů je nesrovnatelně jednodušší. Vyřešit diferenční rovnici v podstatě znamená dosadit do její pravé strany. V případě stochastického systému je k tomu ještě potřeba vygenerovat realizaci příslušného šumu. Neurčitost stochastického systému tak způsobuje také neurčitost výsledku simulace, tj. trajektorií x a y. 23 24 KAPITOLA 4. SIMULACE DYNAMICKÝCH SYSTÉMŮ 4.2 Příklady 4.2.1 Pindyck-Rubinfeldův třírovnicový model Tento lineární model použil Pindyck v učebnici [?] při testování problému optimálního řízení. Model obsahuje 4 výstupní veličiny a 2 vstupní (řídící) proměnné. Je zadán svým vnějším popisem Ct = co + ciYt + c2Ct-i (4.1a) It = io + hYt + Í2{Yt-i-Yt-2)+hRt-A (4.1b) Rt = r0 + riFí + r2(Fí-Fí_1) + r3(Mí-Mí_1)+r4(iíí-i + iíí-2) (4.1c) Yt = Ct + It + Gt (4.1d) Poslední rovnice je ekonomická identita, která se nepočítá jako modelová rovnice a proto mluvíme o třírovnicovém modelu. Proměnné vystupující v modelu jsou: C soukromá spotřeba výstupní veličina / hrubé investice výstupní veličina R úroková míra výstupní veličina Y národní produkt výstupní veličina G vládní spotřeba řídící veličina M množství peněz řídící veličina Parametry modelu (4.1) byly odhadnuty metodou nejmenších čtverců (OLS). Jejich hodnoty jsou: c0 = -9.4540 ci = 0.05410 c2 = 0.9260 i0 = -66.1950 h = 0.16840 i2 = 0.2181 i3 = -11.2650 r0 = -0.5561 n = 0.00051 r2 = 0.0135 r3 = -0.0853 r4 = 0.4259 Pro každou rovnici byly vypočteny reziduálni rozptyly. Jejich hodnoty jsou: s\ = 11.2302 pro rovnici (4.1a) si = 23.8642 pro rovnici (4.1b) s2 = 0.8542 pro rovnici (4.1c) Model má dopravní zpoždění 0, což znamená, že vstup Ut přímo ovlivňuje výstup yt-Lineární soustavu zapíšeme maticově Vt = a0yt + aiyt-i + a2yt-2 + a3yt-3 + a4yt_4 + ß0ut (4.2) 4.2. PŘÍKLADY 25 kde jsme položili Vt OL\ Cr4 = cť It Rt Yt. C2 0 0 0 0 0 0 h 0 0 r4 -T2 0 0 0 0 0 0 0 0" 0 0 h 0 0 0 0 0 0 0 0 0 Ut a2 ßo = Mt - Mt-Gt 1 i 0 0 0 0 0 0 0 -ii 0 0 r4 0 0 0 0 0 0 0 co" 0 0 io r3 0 r0 0 1 0 a0 a3 0 0 0 ci ' 0 0 0 h 0 0 0 n + r2 110 o 0 0 0 0' 0 0 0 0 0 0 0 0 0 0 0 0 Poslední element vektoru Ut je jednička. Tato úprava formálně zbavuje soustavu úrovňových konstant Co,io,fo- Snadno nahlédneme, že model je konstruován dosud nezvyklým způsobem závislosti výstupních veličin. Současná hodnota yt je závislá na současné hodnotě yt- Takové dynamické systémy jsou v ekonomii celkem běžné a označují se jako interdependentní nebo také implicitní systémy. V případě lineárních systémů je lze převést na explicitní. Rovnici (4.2) upravíme následujícím způsobem: yt - a0yt = aiyt-i + a2yt-2 + a3yt-3 + a4yt_4 + ß0ut <í=> (/ - ao)yt = aiiyt-i + a2yt-2 + a3yt-3 + a4yt_4 + ß0ut <í=> yt = (I - a^aiyt-i + {I - a0)"1a2yt-2 + (I - a0)~1a3yt-3+ + (I - a0)~laAyt-A + {I - a0)~lßout <=^ yt = aiyt-i + a2yt-2 + a3yt-3 + a4|/í-4 + b0ut Tím jsme dostali model ARX, který zapíšeme blokově takto yt ' d\ a2 a3 a4 'yt-i "6o í-i í-2 = I 0 0 0 0 0 0 yt-2 yt-3 + 0 0 í-3. 0 0 I 0 JJt-4. 0 ut (4.3) 26 KAPITOLA 4. SIMULACE DYNAMICKÝCH SYSTÉMŮ kde a\ = (I — cto) ld\ = a2 = (/ - ao) 1a2 = a3 = (I - a0) la:i = a,4 = (I — ao) on = b0 = (I - a0) 1ß0 = 0.9904 0 0 0.0152 0.2006 0 0 0.2653 0.0167 0 0.4259 -0.0096 1.1910 0 0 0.2805 0 0 0 -0.0152 0 0 0 -0.2653 0 0 0.4259 -0.0039 0 0 0 -0.2805 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7838 0 0 0 -13.7049 0 0 0 -0.2030 0 0 0 -14.4887 0 0 0.0696 -14.7178" 0 0.2166 -82.5799 0.0853 0.0180 -1.9192 0 1.2862 -97.2977 což můžeme zkráceně napsat jako xt = Axt-i + But (4.4) Rovnice (4.4) je stavovou rovnicí našeho systému. Přitom x G M , u G R . Rozptyly s\,s\, S3 jsou nestrannými odhady rozptylů stochastické složky původního modelu (4.1). Tyto složky předpokládáme nezávislé a normální. Označme e šum původního modelu e ~ N 11.2302 0 0 0 23.8642 0 0 0 0.8542 Po restrukturalizaci se tento šum změnil na šum v, pro který platí Ev = (I-ao)-1Ee = 0 varv = (I — ao)_1vare((I — ao)-1) Matice varv nemusí být nutně diagonální a složky šumu v jsou tedy stochasticky závislé. Pro zkoumání stability systému je podle (3.9) nutné znát spektrum matice A. To bylo 4.2. PŘÍKLADY 27 spočteno programem Matlab: A = -0.4905 +0.4354í -0.4905 - 0.4354í 0.9878 0.9026 0.3937+ 0.5633Í 0.3937 - 0.5633« 0 0 0 0 0 0 0 0 0 0 Simulace systému (4.1) se rozpadá do čtyř relativně nezávislých částí • Příprava řídících strategií, nejlépe ve více variantách • Funkce z = pindyck(x,u) • Funkce [Y,X] = sim(xO,U) • Grafické zpracování výsledků, tj. prezentace matic U,X,Y. 28 KAPITOLA 4. SIMULACE DYNAMICKÝCH SYSTÉMŮ Kapitola 5 OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU 5.1 Optimální estimator Dynamické systémy, které jsme až doposud studovali, měly svůj vnitřní popis odvozený z popisu vnějšího a vektor stavů měl charakter umělé proměnné, která zprostředkovala zmíněný převod. Řada systémů je však často zadána přímo svým vnitřním popisem, ve kterém má vektor stavů svou konkrétní obsahovou interpretaci. Jako příklad je možno uvést celou řadu fyzikálních dějů, které probíhají v izolovaném prostředí. Pozorovateli je skryta přímá informace o hodnotách stavových proměnných. Namísto toho jsou do systému instalovány měřící přístroje, které podávají zprostředkovanou informaci o stavu systému. Vyvstává tak problém, jak z napozorovaných výstupů y určit co nejpřesněji stavy x. Protože se předpokládá, že jak samotný proces tak i měření jsou zatíženy náhodnými chybami, jedná se o poměrně komplikovaný problém na poli matematické statistiky. Tento problém byl zatím úspěšně vyřešen pro případ lineárních systémů a Gaussovských šumů. Příslušný algoritmus pro výpočet neznámých stavů je znám pod názvem Kalmanův filtr. Pro nelineární systémy jsou známy pouze přibližné metody. Předpokládejme neznámý náhodný vektor X G Cnx a vektor naměřených dat z G W12, který je realizací náhodného vektoru Z G Cnz. Předpokládejme dále, že známe simultánní hustotu pravděpodobnosti p(x,z). Předtím, než proběhne měření z, je naše informace o X obsažena pouze v hustotě p(x, y) a pravděpodobnostní rozložení náhodného vektoru X udává hustota «,<*) = /**,*)* Tato hustota se nazývá apriorní hustota náhodného vektoru X. Po provedení měření z je naše informace o X obohacena a pravděpodobnostní rozložení vektoru X dáno podmíněnou hustotou pravděpodobnosti pi(x) =p(x\y) kterou nazýváme apo ster iorní hustotou náhodného vektoru X (viz B.4). Často se namísto hustot určují pouze první a druhé momenty, tj. podmíněné střední hodnoty a podmíněné rozptyly, které mají některé významné vlastnosti. 29 30 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU Tyto úvahy jsou jádrem takzvaného Bayesovského přístupu k optimálnímu odhadu, který hojně používá podmíněné charakteristiky náhodných veličin. Z uvedeného se zdá, že jsme pro tentýž náhodný vektor zavedli dvě pravděpodobnostní míry, které si navzájem konkurují. Ve skutečnosti je jediným skutečným rozložením rozložení apriorní, avšak aposteriorní rozložení má oproti apriornímu některé výhodné vlastnosti právě s ohledem na co nejpřesnější odhad realizace vektoru X, nikoliv hustoty vektoru X. Přesnost odhadu se nejčastěji měří podle kriteria nejmenších čtverců. Definice 5.1. Buď E C W"x množina přípustných estimátorů náhodného vektoru X. Pro každý přípustný estimator y e 8 zaveďme střední kvadratickou chybu vzorcem J{y)=E{{X-y)T{X-y)) Estimator X G 8 nazveme optimálním podle kriteria nejmenších čtverců (least mean square), jestliže platí J(X)=minJ(y) X se také nazývá optimální mean square estimator (MS estimator). Veta 5.2 (optimální apriorní estimator). Střední hodnota EX je optimálním estimátorem na množině všech apriorních estimátorů. Důkaz. Napišme si střední kvadratickou chybu J(y)=E(X-y)T(X-y) = = E XTX - yTE X - E XTy + yTy = = E XTX - 2yTE X + yTy kde jsme použili E y = y protože y není náhodná veličina, ale reálný vektor. Vidíme, že J(y) je kvadratická funkce v proměnné y a globálního minima dosahuje v bodě X, kde diferenciál je roven nule. Zapsáno d.J — = -2EX + 2y = 0 dy nebo také X = EX (5.1) D Tímto jsme odvodili možná samozřejmý a určitě intuitivní výsledek, že při znalosti p(x) je nejlepší konstantou pro odhad X střední hodnota X. Zaveďme chybu estimace X e Cnx jako rozdíl X = X -X Pak platí EI = EI-í = EI-EI = 0, takže střední chyba estimace je rovna nule. Říkáme, že odhad X je nestranný (nevychýlený, unbiased). Ekvivalentně také platí EX = 0 ^ EÍ = EI 5.1. OPTIMÁLNÍ ESTIMÁ TOR 31 Jako míru přesnosti naší estimace spočtěme rozptyl chyby varX = E(X-EX)T(X-EX) = = E XTX = = E (X - xf{X -X) = = E (X - E X)T(X - E X) = = var X Bude-li rozptyl malý, je náš odhad dobrý. Inverzní matici (varX)-1 nazýváme informační matici. Veta 5.3 (optimální aposteriorní estimator). Podmíněná střední hodnota E(X|z) je optimálním estimátorem na množině všech aposteriorních estimátorů. Důkaz. Předpokládejme, že vektor z je neprázdný, tj. nz > 1, jinak by se jednalo o apriorní estimator. Přitom známe hustotu p(x,z). Optimální estimator bude deterministickou funkcí naměřeného vektoru z, tj. X : R""z —> R. Napišme funkcionál J J = / (x — X(z)) (x — X(z))p(x, z) dx dz kde integrujeme přes celý prostor W"x, případně Rraz. Podle řetězového pravidla (B.9) je p(x,z) = p{x\z)p{z) a funkcionál přejde do tvaru J = p(z) / (x — X (z)) (x — X(z))p(x\z) dx dz Protože p(z) a vnitřní integrál jsou nezáporné funkce, můžeme J minimalizovat tak, že minimalizujeme pro každé pevné z vnitřní integrál * = /<*-*M>T<*--*«>**>*» Integrál H můžeme napsat jako H = E((X-X(z))T(X-X(z))\z) = = E (XTX\z) - X(z)TE (X\z) - E (XT\z)X(z) + XT{z)X{z) = = E (XTX\z) - 2X(z)TE (X\z) + X{z)TX{z) Stejným způsobem jako minule dostáváme = -2E (XTX|Z) + 2X(Z) = 0 dX což dává X(z)=E(X\z) (5.2) D 32 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU To už není tak samozřejmý výsledek jako vzorec (5.1), ale je také intuitivní. Jestliže máme k dispozici realizaci z náhodného vektoru Z, pak nejlepším MS-estimátorem náhodného vektoru X je střední hodnota podmíněná Z = z. Rovnice (5.2) vypadá velmi jednoduše, avšak jednoduchost této rovnice skrývá komplikace při jejím konkrétním použití. Ve skutečnosti může být výpočet podmíněné střední hodnoty podstatným problémem. Navíc se předpokládá, že známe přesně simultánní rozložení p{x, z), což nemusí být pravda. V některých případech známe pouze první a druhé momenty. Optimální aposteriorní estimator je nestranný, což dokazuje výpočet EX = E(E(X|Z)) =E(E(X-E(X|Z)|Z)) = = E(X-E(X|Z)|Z) =E(X|Z) -E(X|Z) = 0 Míra spolehlivosti estimátoru je dána rozptylem chyby estimace, který se spočte jako varX = var[(X-X)|Z] = = E [{X - EX|Z)(X -EX\Z)T\Z] Výraz v hranatých závorkách je ale podle definice podmíněný rozptyl var (X\Z), takže rozptyl chyby můžeme psát jako varX = E(var(X|Z)) (5.3) Jestliže X a Z jsou nezávislé vektory, pak p(x, z) = p{x)p{z) a podle (B.13) platí p{x)p{z) E(X\z) = J x- dx = p(z) xp(x) dx = = EX Při nezávislosti tedy aposteriorní a apriorní estimator splývají, což můžeme interpretovat tak, že měření Z nám nepřináší žádnou novou informaci ohledně X. Veta 5.4 (gaussovské pozorování). Nechť platí N ßx ßz Pak aposteriorní náhodná veličina X\z má normální rozložení se střední hodnotou E{X\z)=ßx + TJXZTJ-\z-ßz) (5.4) a varianční maticí var (X | z) = E: Důkaz. Definujme náhodný vektor Y G Cnx + nz jako y V--1V ^XZ^z ^zx (5.5) Y = ßx ßz Hustota pravděpodobnosti náhodného vektoru Y je dána vzorcem p(y) ^{2n)n*+n* |Ej =e 2^ 5.1. OPTIMÁLNÍ ESTIMÁ TOR 33 kde Podle (A.5) je ßy = y-1 - ^y - ßx ßz Yjy — ^x ^xz Ľzx Ľ z D- U l^xz^z -2uz L^zxU 2uz -\- 2uz L^zxU l^xz^z kde Schurův komplement k T:x je D = Ex — TjXZľ~1TjZX. Podle (B.13) je p{y) p(x\z) = P(z) neboli p(x\z) = V^f ^{2n)n*+n* \Ľy 1 -.e 2 \[{y-ßy)T^y 1(y-ßy)-(z-ßz)T^z 1{z-ßz)\ =e 2 V(2tt)^|^|/|S, Napišme si exponent a upravujme ho (V - ßy)T^y\y - ßy) -{z- ßz)TZz\z - ßz) = (z-ßz)TTlz\z-ßz) \[{y-ßy)T^y 1(y-ßy)-(z-ßz)T^z1(z-ßz)] r T -, x - ßx ^ x - _ ßx z - - ßz y z - - ßz = (x-ßx)TD 1(x-ßx) - (x-ßx)TD 1yxz^z1(z- ßz)~ -(z- ßz)TYrzlyzxD-\x - ßx) + (z - ßz)T(Zzl + Ej'E^fl-^^Sj^z - ßz)- -(z-ßz)TYTz1(z- ßz) = = {x- ßx)TD~1(x - ßx) -{x- ßx)TD~1T1xZTlz1(z - ßz)- -(z- ß.fy-^.xD-^x - ßx) + (z- ßzfy-^xD-^xzZ-^z - ßz) = = {x- ßx)TD-\x - ßx) -{x- ßx)TD-l[JlxzL-l{z - ßz)]~ - [E^E;1^ - ßz)]TD~\x - ßx) + [ZxzZ-^z - ß^D-^xz^iz - ßz)] Z prvních dvou členů vzniklého výrazu vytkneme (x — ßx) D-1 a z dalších dvou členů vytkneme [Ľxz'Ľ~1(z — ßz)] D~l a dostaneme (x - ßxi)TD~l [{x - ßx) - F^X"1^ - ßz)]-\LxzTJ~1{z - ßz)] D~l [{x - ßx) - yXz^l{z - ßz)] Nyní můžeme vytknout zprava výraz [(x — ßx) — TJxz^1(z — ßz)], čímž dostaneme [(x - ßx) - yXz^l{z - ßz)] D~l [{x - ßx) - yXz^l{z - ßz)] = = [x- (ßx + yxz^l{z - ßz))] D~l [x - [ßx + yxz^l{z - ßz))] Podle (A.9) platí, že iv i — iv 11 ni \z-"y\ — \z-'z\ \J-J\ = \D\ 34 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU Položíme-li ßx\z = ßx + ^xz^z l{z — ßz), dostáváme nakonec, že V ' ' ^(2ir)n* \D\ tedy podmíněný náhodný vektor X\z má normálni rozložení se střední hodnotou ßx\z = ßx + ^xz^z \z ~ ßz) a rozptylem Jzx D — ^x ^xz^z ^-*z což bylo potřeba dokázat D Právě uvedený výsledek je velmi důležitý a bude všestranně užitečný pro další práci. Ze vzorce (5.4) vyplývá, že var(X|z) < varX. Jestliže nastane speciální případ var(X|z) < var X, pak jsme po provedení měření z více jistí správnou hodnotou X, než jsme byli předtím. To znamená, že aposteriorní náhodná veličina X\z je těsněji rozložena kolem své střední hodnoty E(X|z), než apriorní náhodná veličina X kolem své střední hodnoty EX. Rozptyl náhodné chyby odhadu je roven var((X-E(X|z))|z) = var (X\z) - var (E (X\z)) =var(X|z) (5.6) Důležitá skutečnost je ta, že zatímco E(X|z) závisí na změřené hodnotě z, varianční matice var(X|z) na změřené hodnotě z nezávisí, proto ji můžeme spočítat ještě předtím, než provedeme příslušné měření. Vraťme se nyní opět k obecnému vztahu mezi X a Z, danému simultánní hustotou pravděpodobnosti p(x, z). Ukázali jsme, že nejlepším estimátorem X při daném z je Xms = E(X\z). Výpočet tohoto estimatoru může být v některých případech velmi obtížný. Abychom tento výpočet výrazně zjednodušili, omezme množinu přípustných estimatoru jen na ty estimátory, které jsou lineární funkcí naměřeného vektoru z. Protože třída přípustných estimatoru je nyní omezena, lineární MS estimator nebude obecně tak dobrý jako nelineární MS estimator. Veta 5.5 (lineární MS estimator). Optimální MS estimator na přípustné množině lineárních estimatoru S = {y = Az + b : A e RUx XUz, b G RUx} je estimator Xlms = ßx + ^xz^ž1 (z - ßz) (5.7) kde ßx = EX, ßz = E Z, Y.XZ = cov(X, Z), Y,z = var Z Důkaz. Počítejme opět střední kvadratickou chybu J = E (X - X(Z))T(X - X{Z)) = E (X - AZ - bf {X -AZ-b) což můžeme vzhledem k (A.21) napsat jako J = E (tr (X - AZ - b){X - AZ - bf) = = trE (X - AZ - b){X - AZ - bf = = trE [(X - ßx) -{AZ + b- ßx)][(X - ßx) -{AZ + b- ßx)f = = tr [E (X - ßx){X - ßxf - E (X - ßx){AZ + b - ßxf-- E {AZ + b - ßx)(X - ßxf + E (AZ + b- ßx)(AZ + b- ßxf] 5.1. OPTIMÁLNÍ ESTIMÁ TOR 35 Označme Ex = varX = E (X — ßx)(X — ßx) . Předpokládejme navíc, že matice A je symetrická. Potom platí J = tr [Ex - T,XZAT - AT,ZX + E {AZ + b - ßx){AZ + b - ßxf] = = tr [Ex - 2AĽZX + E (A(Z - ßz) + (Aßz + 6 - jiíx))(A(Z - ßz) + (A^ + b - ßx))T] = = tr [Ex - 2AT,ZX + AY,ZAT + (A^ + 6 - /x;c)(A^ + b - ßx)T] = = tr [Ex + A{üz + ßzßzT)AT + {b- ßx)(b - ßx)T + 2Aßz(b - ßx)T - 2AY,ZX\ Podle (A.16) a (A.17) platí —tr (A(ĽZ + ßzßzT)AT)=2A(YJZ + ^/x/) 9 ——-tr {2Ae,-z p{z\x) = pz\x{z, x) = pv{z - ln -) = i x \0 pro x < e z sdružená hustota p(x, z) je dána vzorcem p(x, z) = p(z\x)p(x), neboli i \ I ie~z Pro x ^ e~z'i 0 < x < 1 p(x,z) = < I 0 jinak Marginalizací dostaneme hustotu p{z) = p z (z): Pz{z) = I p(x,z) dx = l 1 1 _* -e — z X dx = e ■z[lnx]l -z e~ ■*(()-ln e-') = ze ~z pro z >o Nyní už snadno spočteme aposteriorní hustotu p{x\z) jako p(x,z) Px\z{x,z) p{z) le-z ze~z {^ pro x > e~z, 0 < x < 1 0 jinak Aposteriorní estimator je střední hodnotou aposteriorní náhodné veličiny X\Z s hustotou Px\z, tedy Xms(z) = E(X\z) = / xp(x\z) dx = = /' 1 a;— da; = Je-Z a;z "i" — e~ *) = 1- e"2 38 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU Spočteme rozptyl aposteriorní chyby. Podle (5.3) máme varX = E(var(X|Z)) = r r (x — E (X\z))2p(x\z) dx dz = x2,p(x\z) dx — E (X\z)2 dz = ~ [l j-Lte-——-----p-^-Y---------p(xt \Dt-i)dxt = J p(vt\ut, Dt-i) _ p(xt+i,yt\ut,Dt-i) p(yt\ut, A-i) Tímto způsobem jsme ale ztratili aposteriorní hustoty, které jsou obecně lepší než hustoty apriorní. 5.3 Kalmanův filtr Kalmanův filtr je analogií výše uvedené nelineární rekurze pro speciální případ, kdy jak stavy x tak výstupy y jsou rozloženy normálně. Pak je podle 5.4 apriorní i aposteriorní hustota opět normální. Normálně rozdělenou náhodnou veličinu plně určuje její střední hodnota a rozptyl. Stačí tedy po celou dobu výpočtu sledovat pouze tyto dvě číselné charakteristiky, což v konečném důsledku umožňuje nahradit funkcionální rekurze rekurzemi algebraickými. Celý výpočet se tak snadno implementuje na počítači. Mají-li být stavy rozloženy normálně, musí být příslušný dynamický systém nutně lineární s gaussovskými šumy. To je největší slabinou Kalmanova filtru, kterou se snaží odstranit některé novější algoritmy, které jsou ovšem buď nepřesné nebo podstatně složitější. Veta 5.10 (Kalmanův filtr). Buď dán lineární diskrétní stochastický systém, popsaný rovnicemi xt+i = Axt + But + vt (5.20a) yt = Cxt + Dut + Wt (5.20b) splňující tyto podmínky: 42 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU (i) xo ~ N(/i0, E0) (ii) vt ~ -/V(0, E^) pro všechna í (iii) wt ~ -/V(0, Ew) pro všechna í (iv) ßutiüf = 0 (v) Exovf = 0 kde vektor /xo a matice Eo,Eí,,Ew jsou známé. Dále nechť jsou známa data Dt. Označme xt\k = xt\Dk- Pak apriorní estimator xt\t-\ a aposteriorní estimator xt\t jsou normální pro všechna t xt\t_x ~ A^(/xŕ|ŕ_x, Sŕ|ŕ_!) (5.21) xt\t ~ iV^it, Et|t) (5.22) kde střední hodnoty ßt\t-iilk\t a varianční matice Sr|É_1, Et|t se vypočtou podle vzorců A**|*—i + Kt(Vt - Cßt\t-\ - Dut) (5.23) £t|t-i-tftC£t|t-i (5.24) ^t|t + 5% (5.25) AEt|tAT + E„ (5.26) E^C^CE^C^ + E^)"1 (5.27) Všimněme si rekurentní charakter vztahů (5.23)-(5.27). Na základě počáteční podmínky i a dat Dt jsme schopni postupně spočítat odhady íCi|o5 ^i|i) ^2|i) ^2125 • • • ixt\t-iixt\t- Přitom ve výpočtu vystupují pouze první a druhé momenty těchto veličin. Rovnice (5.23),(5.24) tvoří datový krok algoritmu, což je přechod od xt\t-\ k xt\t. Naproti tomu rovnice (5.25),(5.26) tvoří predikční krok, tedy přechod od xt\t k xt+í\t. Dosadíme-li (5.27) do rovnice (5.24)(viz (5.29)), dostaneme diferenční rovnici pro Et|t, která je maticovou diskrétní variantou Riccatiho diferenciální rovnice tvaru x = ao + a\x + Důkaz. Předpokládejme, že platí (5.21). Označme yt\k = yt\Dj,. Pak Vt\t-i =Vt\Dt-i = = (Cxt + Dut + wt)\Dt_l = = Cxt\t_i +wt\Dt = = Cxt\t_x + wt Protože yt\t-i je lineární kombinací náhodných vektorů xt\t-\,wt, které jsou podle předpokladu normální, je yt\t-i podle (B.6) normální se střední hodnotou Cßt\t-\ + Dut a rozptylem var (yt\t-i) = var (Cxt\t-i + Dut + wt) = = Cvar (xt\t-i)CT + varwt = = C^í|í-iC + ^w ßt\t = sí|í = ßt+\\t = s*+i|í = Kt = 5.3. KALMANŮV FILTR 43 Spočítejme kovarianční matici cov^t-i^it-i) = cxw{xt\t_i,Cxt\t_i +Dut + wt) = var (xt\t_i)CT + cov(xt\t_i, wt) = sí|í-iC Spojený náhodný vektor (xt\t-i,yt\t-i) je normálně rozdělen xt\t-i Vt\t-i N ßt\t-i jí|í-i sí|í-iC T .CSt|t_i CStit.iC +s Jsou tedy splněny předpoklady věty 5.4 a pro podmíněnou náhodnou veličinu xt\t-i\Vt = xt\(yt,Dt-i) =xt\(yt,ut,Dt-i) = xt\Dt = xt\t platí, že je normální s parametry fJ>t\t = ßt\t-i + St|t_1CT(CEt|t_1CT + X»)"1^* - Cßt\t-i - Dut) ^t|í = ^t|í-i ~~ ^í|í-iC (C^í|í-iC + ^w) C^t|t-i Položíme-li pak můžeme psát (5.28) (5.29) Kt = £t|t_iC (CEti^iC + Sw) í"í|í = í"í|í-i + ^í(ž/í - Cju*|*-i - £>«*) sí|í = sí|í-i - ^íCSt|t_i což jsou rovnice (5.23) a (5.24), které tvoří datový krok algoritmu. Nyní odvoďme predikční krok. xt+1\t = (Axt + But + vt)\Dt = = A(xt\Dt) + But + (vt\Dt) = = Axt\t + But + vt Vidíme, že xt+\\t je lineární kombinací xt\t a Vt, které jsou normální. xt+\\t je tedy podle (B.6) také normální s parametry ßt+\\t = Aßt\t + But sí+i|í = AY t. Výsledkem je tedy podmíněný náhodný vektor xt\^. V porovnání s apriorním odhadem xt\t_i a aposteriorním odhadem xt\t máme pro odhad stavu k dispozici více informace a výsledný estimator tedy bude ještě lepší než oba uvedené typy. Naproti tomu smoothing nelze provádět v aktuálním čase t, zatímco rekurzivní odhady ano. Smoothing se tedy nejspíše uplatní tam, kde nás kromě aktuálního stavu xt zajímají také minulé stavy xq, ■ ■ ■ ,xt-i. Jak si ukážeme, lze odhady x0\t,... ,xt_i\t získat úpravou již existujících aposteriorních odhadů x0\0,... ,xt_i\t_i. Algoritmus přitom postupuje zpětně v čase. Veta 5.12. Nechť t < N. Pro podmíněnou hustotu p{xt\D^) platí: p(xt\DN) =p(xt\Dt) p(xt+i\xt,ut) p(xt+i\D N) p(xt+i\Dt dx í+i Důkaz. Počítejme p(xt\DN) = / p(xt,xt+i\DN)dxM p(xt,xt+i\DN) = p(xt\xt+i, DN)p(xt+\\DN) (5.35) (5.36) (5.37) Rozberme nyní podrobněji podmíněnou hustotu p{xt\xt+i,Dj^). Veškerá informace uložená v datech {yt+i,ut+i, ■ ■ ■ ,Vn,un} potřebná k predikci xt je koncentrovaná v hodnotě xt+i, což vyplývá z konstrukce modelových rovnic (5.12). Můžeme tedy prohlásit, že p(xt\xt+i,DN) =p(xt\xt+i,Dt) (5.38) 46 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU přitom p(xt\Dt) p(xt\xt+i, Dt) =p{xt+i\xt,Dt) = p{xt+i\xt,ut) p(xt+i\Dt) p(xt\Dt) p(xt+i\Dt) dosazením do (5.37) a (5.36) nakonec dostáváme p{xt\DN) = / p{xt+i\xt,ut)—— in Mxt+i\DN)dxt+i = J P{Xt+i\L>t) = p(xt\Dt) / p(xt+i\xt,ut) , t+1 | J^ dxt+i J p(xt+i\Dt) což jsme měli dokázat D Povšimněme si, že ve výpočtu nikde nevyužíváme naměřená data. Využíváme je pouze zprostředkovaně, a to prostřednictvím estimátorů xt+í\N a xt\t. Tyto dva estimátory jsou také postačující pro výpočet estimátorů xt\N- Rekurzivní vztah (5.35) se v případě gaussovského systému opět rozpadne na maticové vzorce pro střední hodnotu a rozptyl estimátorů. Tyto vzorce, známé pod názvem Rauch-Tung-Striebel Smoother, uvedeme bez důkazu. Veta 5.13 (Rauch-Tung-Striebel Smoother). Pro lineární systém 5.20a se odhad xt\N ~ N(/j,t\N, £t|jv) vypočte podle vzorců ßt\N = ßt\t + Ft(ßt+i\N - ßt+\\t) (5-39) St|jv = ^t\t - Ft{^t+i\t - ^t+i\N)FtT (5.40) Ft = St|tArE^1|t (5.41) 5.5 Rozšířený Kalmanův filtr Myšlenka rozšířeného Kalmanova filtru (Extended Kaiman Filter) spočívá v tom, že nelineární dynamický systém v každém kroku přibližně nahradíme lineárním dynamickým systémem a na tento pak aplikujeme obyčejný Kalmanův filtr. Výsledky věty 5.10 však nelze použít přímo, protože linearizace vede na poněkud jiný typ DS. Lemma 5.14. Buď dán lineární diskrétní stochastický systém, popsaný rovnicemi xt+i = Atxt + bt + vt (5.42) ílt = Ctxt + dt + wt (5.43) splňující podmínky (i)-(v) z věty 5.10. Pak pro apriorní estimator xt\t_i ~ N(/j,t\t_1,Ylt\t_1) a aposteriorní estimator xt\t ~ N(p,t\t, Et|t) platí ßt\t = ßt\t-i + Kt(yt - Ctßt\t-\ - dt) (5.44) St|t = Et|t-i - KtCt^t\t-i (5-45) ßt+i\t = Atßt\t + bt (5.46) Et+1|t = AtZt\tAtT + i:v (5.47) Kt = ^t\t-iCtT(CtXt\t-iCtT + E«,)-1 (5.48) 5.5. ROZŠÍŘENÝ KALMANŮV FILTR 47 Důkaz. Důkaz je pouze drobnou modifikací důkazu věty 5.10. D Vezměme obecný nelineární dynamický systém ve stavovém tvaru Xt+l = f(xt, Ut) + vt yt =g(xt,ut) + wt Linearizujeme-li funkce f(x, u) a g(x, u) proměnné x podle Taylorova rozvoje v bodě x = £, dostaneme f(x,u) = f&u) + |£(0(z -0 +0{x2) (5.49) g{x, u) = g(£, u) + ||(0(z - 0 + 0(a;2) (5.50) Aproximace bude tím lepší, čím menší bude rozdíl x — £. Zanedbáme-li členy 0(a;2), můžeme původní systém přibližně nahradit lineárním systémem df xt+i = m,ut)+At(xt-o+vt At = ú^ (5-51) Vt=g&ut)+Ct(xt-t)+wt Ct = ||(£) (5.52) Předpokládejme, že máme k dispozici estimator xt\t-\ prostřednictvím jeho momentů ßt\t-ii ^t\t-i-Linearizujme výstupní rovnici v okolí bodu ßt\t-\- dg yt = Ctxt + (g(nt\t-i,ut)-CtHt\t-i)+wt kde Ct =—(ßt\t-i) (5-53) Pak podle lemmatu 5.14 je datový krok dán rovnicí ßt\t = lh\t-i + Kt(yt-g(ßt\t-i,ut)) (5-54) sí|í = sí|í-i -KtCt?■ R je tzv. Diracova funkce , definovaná limitou ö(x) = \imöh(x), öh(x) = \hPrO0 0 tak, aby součet jednotlivých vah dával jedničku. Empirická hustota je potom n Vn{x) = ^2lwiö(x-xi) (5.70) í=l Rozložení náhodného vektoru x je tedy popsáno náhodným výběrem Sn a vektorem vah wn = (w\,... ,wn). Připomeňme obecné rekurzivní vztahy pro nelineární estimaci: p(xt\Dt) = fyfXt^ p(s*lA-i) (5.71) p{yt\ut, Dt-i) p{xt+i\Dt) = p(xt+i\xt,ut)p(xt\Dt)dxt (5.72) p{xt\DN) =p{xt\Dt) / p{xt+i\xt,ut) , t+1 | J^ dxt+i (5.73) J P(Xt+l\-L>t) Pracujeme-li s empirickými hustotami, pak aplikace těchto rekurzivních vztahů vyústí v algoritmus přepočtu vah Wi a vzorků Xi, který má název Weighted Bootstrap Algoritm. Vysvětlíme si ve stručnosti princip tohoto algoritmu: Buď dána apriorní hustota n Pn{xt\t-l) = ^ Wi(t\t - 1)5{X - Xí{t\t - 1)) í=l V datovém kroku hledáme aposteriorní hustotu pn(%t\t)- Položíme Xi(t\ť) = Xí{t\t — 1) a váhy Wi(t\ť) získáme pomocí (5.71). Clen p(yt\ut, -Dt-i) funguje ve vzorci (5.71) jako normalizační konstanta. Když tento člen vypustíme, obdržíme váhy Wi(t\ť) = p{y{ť)\xí{t\t - 1), u(t)) Wi(t\t - 1) (5.74) z nichž normalizací dostaneme výsledné váhy wdMt) = řS§m ("5) Tím je zakončen datový krok algoritmu. Při predikčním kroku nepoužijeme vztah (5.72), ale pro modelový krok využijeme vztahu xt+i\t = f(xt\t,u(t)) + v(t) který aplikujeme na vzorky Xí{t\ť). Dostaneme pak Xi(t + l|í) = f(xi(t\t), u(t)) + Vi(t) (5.76) kde Vi(t) jsou náhodně generované vzorky šumu procesu Vt- Váhy se při predikčním kroku nemění, tj. Wi(t + l|í) = Wi(t\t). Takto jsme dostali predikci xt+í\t. Konstrukce smoothingu má stejnou filozofii jako datový krok. Vzorky zůstávají stejné Xi(t\N) = Xí{t\ť) a upravují se jen váhy. Vypuštěním členu p{xtJri\Dt) ze vzorce (5.73) dostaneme zmenšené váhy Wi(t\N) podle vzorce n wi(t\N)=wi(t\t)^wj(t + l\N)p(xj(t + l\N)\xi(t\N),u(ť)) (5.77) i=i 50 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU a váhy Wi(t\N) dostaneme normalizací w'm) = nMm (5-78) V datovém kroku jsme potřebovali znalost výstupní rovnice skrze hustotu p(y\x, u), v predikčním kroku jsme zase použili funkci /. Při smoothingu se neobejdeme bez podmíněné hustoty p(xt+i\xt,ut), která je určená stavovou rovnici. Poslední nevyřešenou otázkou je zahájení algoritmu. Předpokládáme, že rozložení stavu Xq je známé. Pak v souladu s (5.19) položíme a?o|-i = xo- Potřebujeme vzorky Xi(0\ — 1) a váhy Wi(0\—1). Nejčastěji se vzorky získají generátorem náhodných čísel s požadovaným rozložením. Váhy jsou potom rovnoměrné, tedy Wi(0\ — 1) = l/n. Jiná možnost je zvolit rovnoměrnou síť bodů Xi(0\ — 1) ve vhodně zvolené oblasti a váhy volit Wi(0\ — 1) = po(xi(0\ — 1)) a potom je normalizovat. Dosažené výsledky shrneme v definici. Definice 5.16 (Weighted Bootstrap Algorithm). Buď dán dynamický systém Xt+l = f(xt,Ut,Vt) yt = g(xt,ut,wt) s počátečním stavem xq, který má empirickou hustotu n Pn(x) = Y, ™i(°l - 1)6(X - Xi(°\ - !)) í=l Buď dána data Dt (viz (5.13)). Pak odhady xt\t_i,xt\t,xt\N s empirickými hustotami n Pn{xt\t-l) = ^2 Wi(t\t - 1)0(X - Xi{t\t - 1)) í=l n Pn{xt\t) = ^2wi(t\t)Ö(x - Xi(t\t)) í=l n Pn(Xt\N) = J2Wi(t\N)6(x - Xi(t\N)) í=l jsou spočteny algoritmem Weighted Bootstrap, právě když platí Xi(t\ť) = Xi(t\t - 1) (5.79) Wi(t\t) = p{y{ť)\xí{t\t - 1), u(t)) Wi(t\t - 1) (5.80) a"(t|i) = EřfU ("81) Xi(t + l|í) = f(xi(t\t),u(t)) + Vi(t) (5.82) iuí(í + 1|í)=«;í(í|í) (5.83) Xi(t\N) = Xi(t\t) (5.84) n wi(t\N)=wi(t\t)^wj(t + l\N)p(xj(t + l\N)\xi(t\N),u(ť)) (5.85) i=i 5.7. SHRNUTÍ 51 5.7 Shrnutí Optimální estimator náhodného vektoru X na základě dat Z je podmíněná náhodná veličina X\Z případně podmíněná střední hodnota E(X|Z). Je-li sdružený vektor (X,Z) rozložen normálně, pak optimální estimator je lineární. Lineární estimator se dá použít i pro jiná rozdělení, je však obecně méně přesný než estimator X\Z. Při odhadu stavů dynamického systému (5.12) používáme apriorní odhad xt\t_i, aposte-riorní odhad xt\t a smoothing xt\N. První dva estimátory se vypočtou algoritmem dopředně rekurze pomocí vzorců (5.18) a (5.17). Smoothing probíhá zpětnou rekurzí s využitím vzorce (5.35). Pracujeme-li s lineárním gaussovským systémem (5.20a), pak odhady jsou normální a stačí spočítat jejich střední hodnoty a rozptyl. Dopředný algoritmus pro xt\t_i a xt\t se jmenuje Kalmanův filtr a zpětný algoritmus pro xt\^ se jmenuje Rauch-Tung-Striebel smoother. Pro nelineární systémy s diferencovatelnými funkcemi /, g se může použít rozšířený Kalmanův filtr. Ten v každém kroku lokálně linearizuje funkce /,j a na takto linearizovaný systém použije Kalmanových vzorců. Výsledné estimátory pak sice nejsou optimální, zato jejich výpočet je poměrně efektivní. Přesný nelineární estimator je dán vzorci (5.18), (5.17), (5.35). Jejich numerický výpočet je však se současnými prostředky nedosažitelný. Proto saháme k přibližným metodám. V Bootstrap algoritmu je rozvedena myšlenka reprezentace hustoty pomocí váženého náhodného výběru, tj. metoda Monte Carlo. 5.8 Příklady k procvičení Přiklad 5.17. Napište v systému Matlab funkci kf, která provede algoritmus Kalmanova filtru na časovém horizontu 0... N pro systém (5.20a). Funkce bude mít prototyp [Mp, Sp, Mf, Sf] = kf(A, B, C, D, Sv, Sw, mO, SO, Y, U), kde uvedené matice mají následující význam: Mp Matice typu nx x N + 1, obsahující ve svých sloupcích predikce /x01_i, • • •, Vn\n-i- Sp Matice typu nx x N + 1, obsahující ve sloupcích varianční matice Sŕ|ŕ_1. Mf Matice typu nxxN+l, obsahující ve svých sloupcích filtrované střední hodnoty /x0|o> • • • ,Vn\n- Sf Matice typu n2x x N + 1, obsahující ve sloupcích varianční matice Y,t\t. A Matice A typu ' ^X ' ^X i definující stavovou rovnici (5.20a). B Matice B typu nx x nu, definující stavovou rovnici (5.20a). C Matice C typu ny x nx, definující výstupní rovnici (5.20b). D Matice D typu ny x nu, definující výstupní rovnici (5.20b). Sv Varianční matice stavového šumu Y,v typu nx x nx, symetrická a pozitivně semidefinitní. Sw Varianční matice výstupního šumu Sw typu nyxny, symetrická a pozitivně semidefinitní. mO Střední hodnota počátečního stavu Xq typu nx x 1. 52 KAPITOLA 5. OPTIMÁLNÍ ODHAD STAVŮ SYSTÉMU SO Variance počátečního stavu xq typu nx x nx. Y Datová matice výstupů typu ny x N + 1, obsahující ve sloupcích t výstupy yt. U Datová matice vstupů typu nu x N + 1, obsahující ve sloupcích t vstupy ut. Kapitola 6 IDENTIFIKACE DYNAMICKÝCH SYSTÉMŮ 6.1 Úvod Uvažujme diskrétní dynamický systém závislý na parametru 9: xt+i = f(xt, ut, vt, 9) = fe(xt, ut, vt) (6.1a) yt = g(xt, ut, wt, 9) = g0(xt, ut, wt) (6.1b) Parametr 9 G 6 je veličina, na které jsou závislé vlastnosti systému (6.1). Množina 6 se pak nazývá parametrický prostor. Přiklad 6.1. Mějme dynamický systém, o němž víme, že je lineární: xt+i = Axt + But + vt yt = Cxt + Dut + wt, ale neznáme matice A, B,C, D. Tento systém je závislý na parametru 9 = [A B C D] a parametrický prostor 6 je potom 6 = Rn*xn* x Rra=xra« x M.nyxn* x M.nyxn^ s dimenzí Většinou se stává, že u dynamického systému neznáme konkrétně funkce fo,go, ale známe strukturu systému f,g. Struktura systému spolu s vektorem 9 jednoznačně určuje vlastnosti systému. Chybějí-li tedy informace o hodnotě 9, je žádoucí je doplnit. Jediným dostupným zdrojem informací přitom bývá chování systému v minulosti, tedy historie vstupů a výstupů. Techniky, které se zabývají odhadem parametrů systému na základě známé historie se nazývají techniky identifikace. Označíme-li 9 odhad parametru 9 a Dn = {uq, yo, ■ ■ ■, un, yn} historická data, pak identifikace systému (6.1) je úloha nalézt při známé struktuře f,g odhad 9 = F(Dn) parametru 9 tak, aby se co nejméně lišil od skutečného parametru 9. 6.2 Klasický přístup Při klasickém přístupu se neznámý parametr 9 uvažuje jako reálný vektor dimenze p. Parametrický prostor 0 je potom roven 0 C W. 53 54 KAPITOLA 6. IDENTIFIKACE DYNAMICKÝCH SYSTÉMŮ Při identifikaci je cílem dosáhnout maximální shody mezi odhadem 9 a skutečnou hodnotou parametru 9. Protože však rozdíl 9 — 9 není přímo měřitelný, používají se různá kriteria, která jsou vodítkem při určení 9. 6.2.1 Metoda nejmenších čtverců (LS) Metoda nejmenších čtverců (least squares method, LS) říká, že odhad parametru je tím lepší, čím lepší je shoda modelu se skutečností. Označíme-li výstup modelu jako ýt = ýt(9), pak shodu modelu se skutečností měříme podle kriteria nejmenších čtverců jako •/(*) = £ Iv*-fcWI5 (6.2) í=0 Namísto neměřitelného rozdílu 9 — 9 tak pracujeme s měřitelným rozdílem yt — ýt- Optimální odhad 9 bude takový vektor, pro který bude kriterium J minimální, tj. J (6) = minJ(0) 6>ee (6.3) Funkce J je zobrazení J : programování. L Tuto funkci minimalizujeme metodami matematického Přiklad 6.2 (Lineární systém). Pro jednoduchost předpokládejme nejprve lineární systém s pozorovatelnými stavy, neboli systém xt+i = Axt + But + vt yt = xt (6.4a) (6.4b) Máme identifikovat matice A, B na základě historie systému, zachycené v datech Dn = {uo,yo, ■■■ ,Un, yn} = {uo,xo, ■ ■ ■, un, xn} Kriteriální funkce J bude nyní vážit rozdíl mezi stavem xt+i a jeho 0-předpovědí Axt + Buf. TI—1 J(9) = J(A, B) = J2 \xt+i - Axt - But í=0 Označíme-li et = xt+i — Axt — But a matice et, xt+i,A, B rozepíšeme po řádcích et = ei,t e-"nx%t Xt+1 = Xl,t+1 xnx,t+l A = a\ B = pak platí TI—1 ti—1 nx nx ti—1 j(a b) = Y1 ie*i2 = Yl Yl 4 = Yl Y^+i - aixt - bíurf = Y J^a^ bí í=0 j = l j = l í=0 j = l í=0 6.2. KLASICKÝ PŘÍSTUP 55 kde jsme položili Jj(aj> bj) = (xj,t+i - ajXt - bjUtf Úloha minimalizovat J (A, B) se nám nyní rozpadla na nx úloh minimalizovat J j (a j, b j) pro j G {1,... ,nx}. Tyto dílčí úlohy jsou úlohami lineární regrese (viz B.6) ve tvaru Xj,t+1 = CLjXt + bjUt pri označeni xJ,i x 3,n. ßj = iaJ h\ u = T T Xo UO T T Xfi—1 Un—i (6.5) můžeme tyto úlohy zapsat maticově Zj = Ußj Podle (B.18) je optimálním odhadem vektoru ßj vektor ßj, daný vztahem ßj = {UTU)-lUTz3 Odhad 9 = [A B] složíme s odhadů ßj-. ~ßl (6.6) (6.7) [Á Ě] = Á Tu)-luT x\\ • • • xn \ ,T = [ßi ••• ßnf=[{UTU)-1UTz1 ••• (UTU)-^znx]T = ■ znx]TU(UTU)-1 = u(uTuyl =xu(uTuyl = {{utu)-Hjt[Zl ■■■ znx]}=[Zl T Celkem tedy máme kde jsme položili Xr, Trn-1 [Á B] = XU{UTU) (6.8) X = X\ xr. U = T T Xo1 «O T T %n—1 Un—i Přiklad 6.3 (Nelineární systém). Pro jednoduchost předpokládejme opět systém s pozorovatelnými stavy, s konstantními parametry a s aditivním šumem xt+i = f(xt, ut, 9) + vt yt = xt (6.9a) (6.9b) 56 KAPITOLA 6. IDENTIFIKACE DYNAMICKÝCH SYSTÉMŮ Na základě známých dat Dn máme odhadnout parametr 9 tak, aby součet čtverců TI—1 J(9) = Y,\%t+i-f(xt,ut,9)\2 (6.10) í=0 byl minimální na W. Jedná se o hledání volného extrému funkce p proměnných. Je-li funkce / třídy C1, pak z analýzy víme, že je-li bod 9 bodem minima funkce J, pak platí, že diferenciál v tomto bodě je nulový a tedy všechny parciální derivace jsou nulové. Spočteme parciální derivace: dJ ^j^^(xt,ut,6)(xjtt+i- fjíxt,^)) (6.11) r)9i ^ ^ 39, *=0 j=l Označíme-li M0)= [^(xt,ut,e) bt{9) = xjít+i - fj(xt, ut, 9) TI—1 F(9) = J2M0M9) í=0 pak můžeme podmínku |^ = 0 zapsat jako soustavu nelineárních rovnic řádu p F{9) = 0 (6.12) Tuto soustavu pak můžeme řešit některou ze známých numerických metod, jako je Newtonova metoda, Gauss-Seidelova metoda, apod. (viz [?]). V Matlabu je k dispozici funkce f solve, do které za funkci fun dosadíme funkci F. Jiný způsob je použití metod numerické optimalizace přímo pro minimalizaci kriteria (6.10). Známejšou metody nejrychlejšího spádu, metoda konjugovaných gradientů, metoda proměnné metriky, Newtonova metoda (viz [?]). V Matlabu minimalizaci provádí funkce fmins, do které jako argument F dosadíme funkci J podle (6.10). 6.2.2 Metoda maximální věrohodnosti Není-li kriterium nejmenších čtverců J{9) funkce kvadratická v 9, je mnohdy výhodné místo kriteria nejmenších čtverců uvažovat kriterium maximální věrohodnosti. 6.3 Bayesův přístup Bayesova metoda předpokládá, že parametr 9 je náhodný vektor, který se vyvíjí v čase. Jinými slovy posloupnost {9t}^Zo je náhodný proces. Takový proces je popsaný například sdruženou hustotou p(9o, ■ ■ ■, 9n-\). Tento náhodný proces je před identifikací již známý (je známa hustota pravděpodobnosti p). Podoba náhodného procesu {9t}™Zo patří mezi apriorní informace, které ovlivní konečný odhad. Při konstrukci modelu vývoje parametru 9 máme možnost ovlivnit závislost mezi parametry v různých obdobích i závislost mezi jednotlivými složkami parametru 9. Můžeme například předepsat, že parametry 9t a 9t+\ budou silně korelovány, zatímco složky #i;t,..., 9pj 6.3. BAYESŮV PŘÍSTUP 57 budou nezávislé a podobně. U Bayesovy metody se předpokládá, že náhodný proces {#t}™=0 je zadán jako dynamický systém (ve stavovém tvaru) et+i = h(et,et) (6.13) s počáteční podmínkou 9o = (, kde ( G Lp je náhodný vektor se známým rozložením a et je stochastická složka náhodného procesu (6.13) (bílý šum). Mezi nejčastější (a také nejjednodušší) modely vývoje parametrů patří model náhodné procházky (random walk) ve tvaru 6t+1 = 6t + et (6.14) Model náhodné procházky vykazuje některé charakteristické vlastnosti: (i) E(0t)=E(0o) proíe{0,...,n-l}. (ii) var (9t) = var (6>0) + S*=o var (ei) Pro t e {0,... ,n-1}. (iii) cov(0t, 6t+s) = var (9t) pro t G {0,..., n - 1}, s G {0,..., n - t - 1}. Předpokládejme, že je znám model vývoje parametrů (6.13) a je dán parametrický dynamický systém (6.1). Máme tedy systém ve tvaru xt+i = f(xt,ut,vt,9t) 6t+i = h(6t,et) yt = g(xt,ut,wt,9t) Definujme náhodný vektory Zt,et předpisem Zt et yrix+p následujícím způsobem: f(xt, ut, vt, 9t) h{9t,et) Definujme dále funkci F : Rn*+P x Rra- x Rra- F(zt,ut,et) = a nakonec funkci G : Rn*+P x RUu x RUw —> Rny předpisem G(zt, ut, wt) = g(xt, ut, wt, 9t) Potom dynamický systém (6.15) napíšeme pomoci (6.16), (6.17), (6.18) jako zt+i = F{zt, ut, ef) yt = G{zt,ut,wt) (6.15a) (6.15b) (6.15c) (6.16) (6.17) (6.18) (6.19a) (6.19b) Systém (6.19) je neparametrický dynamický systém ve stavovém tvaru. Úloha optimálního odhadu stavu zt systému (6.19) je ekvivalentní s úlohou optimálního odhadu stavu xt systému (6.15) a současně optimálního odhadu parametrů 9t systému (6.15). Tímto jsme převedli odhad parametrů na odhad stavů, což je Bayesova metoda odhadu parametrů dynamického systému. Algoritmicky se tato úloha dá realizovat například použitím Rozšířeného Kalmanova filtru. 58 KAPITOLA 6. IDENTIFIKACE DYNAMICKÝCH SYSTÉMŮ Přiklad 6.4. Vezměme nelineární parametrický dynamický systém s pozorovatelnými stavy a se známým modelem vývoje parametrů xt+i = f(xt,ut,vt,9t) 6t+i=h(6t,et) ílt = Xt Pro tento systém nechť známe data Dn = {uo,yo, ■ ■ ■ ,Un,Vn}- Odhadněme parametry 6t Bayesovou metodou. Zaveďme výstupní proměnnou rjt = Xt+i a stav zt = Qt- Dále definujme funkci g jako g(zt, ut, vt, t) = f(xt, ut, vt, zt) Funkce g je jednoznačně určena funkcí / a daty Dn. Nyní má systém podobu zt+i = h(zt,et) rjt = g(zt,ut,vt,ť) Jedná se tedy o neautonomní systém s nepozorovatelnými stavy zt s šumem procesu et a výstupním šumem vt, který byl původně šumem procesu. Pro odhad stavu Zt (parametru 6t) máme k dispozici vstupy uo,... ,un a výstupy í?o, • • •, f]n-i- Časový horizont je tedy zkrácený o jednu periodu. Přiklad 6.5. Buď dán skalární lineární systém x{t + 1) = a{t)x{t) + b{t)u{t) + v{t) y{t) = c{t)x{t) + d{t)u{t) + w{t) kde v(t) ~ N(0, a%),w(ť) ~ N(0, a^) s časově proměnnými parametry 0(t) = 'a(í) 6(í) c(t) d(í)_ Předpokládáme model vývoje parametrů ve tvaru náhodné procházky ôt+i = e{t) + e{t) kde e (i) ~ N(0, a"l) je bílý šum. Známe počáteční podmínky x0 ~ N{pt0,al) 9o~N{vo,tŮ) Sestrojme algoritmus pro optimální odhad parametru 6{ť) s použitím Bayesovy metody a Rozšířeného Kalmanova filtru. Definujeme si stavový vektor ~x{t)~ a{t) z{ť)= b{t) c(í) ß{t). 6.3. BAYESŮV PŘÍSTUP 59 a vektor šumu e(í) = ~v(ty £l(í) e2(í) £3(í) e4(í) Dále definujeme funkci F : vzorcem F(*(í),u(r),e(í)) = ^i(í)%(í)+2;3(í)«(í)+ei(í)' z2{ť) + e2{ť) %(í) + e3(í) z4(í) + e4(í) %(í) + e5(í) a funkci G : Dynamický systém vzorcem G(z(ť),u(ť),w(ť)) = Zi(ť)Z4(ť) + z5(t)u(t) + «;(*) z(í + l) = F(z(í),u(í),e(í)) y(t) = G(z(t),u(t),w(t)) je nyní přepisem původního lineárního dynamického systému. V rozepsané podobě je to systém 21 (i + 1) = zi(t)z2(t) + z3(ť)u(ť) + ei(í) z2{t + l)=z2{ť) + e2{ť) z3(t + 1) = z3(t) + e3(t) ^(í + 1) =zA{t) + eA{t) z5(t + l) =z5(ť) + e5(t) y(t) = zi(t)z4(t) + z5(t)u(t) + w(t), který zřejmě není lineární. Pro optimální odhad stavu z(ť) se tedy hodí Rozšířený Kalmanův Filtr. Matice A(t),C(t) jsou: A(í) z2(í) zi(í) u(í) 0 0' 1 0 0 0 0 0 1 0 0 0 0 0 10 0 0 0 0 10 0 0 0 0 1 C(t) = [zA{ť] 0 0 zi(t) u(t)] 60 KAPITOLA 6. IDENTIFIKACE DYNAMICKÝCH SYSTÉMŮ Kapitola 7 TEORIE CHAOSU 61 62 KAPITOLA 7. TEORIE CHAOSU Dodatek A Přehled maticové algebry A.l Blokové matice Necht matice A je zadána blokově A = Au A12 A21 A22 Pak definujeme Schurův komplement A22 jako ^22 = ^22 - A2iAn1A12 a Schurův komplement A\\ jako Du = Au - Al2A22lA2i Inverze k A pak může být psána ekvivalentně jedním ze tří vzorců A~l = A~l = Anl + AnlA12DŽ2lA2iAnl -AnlA12Drl 1 a Ä-i '22 -D22 A2iAn D 22 D 11 -Dn AÍ2A22 -A22 A2\Dll A22 + A22 A2\Dll AUA22 A~l = DÜ1 -A22 A2\Dll -Dn AÍ2A22 D 22 Porovnáním těchto tří tvarů se dá odvodit lemma o inverzní matici (Anl - A12A22A2i)-1 = An - AuA12(A2iAnA12 + A22l)-lA2iAu Dále platí \A\ = \Au\ \D22\ \A\ = \A22\\Du\ (A.l) (A.2) (A.3) (A.4) (A.5) (A.6) (A.7) (A.8) (A.9) 63 64 DODATEK A. PŘEHLED MATICOVÉ ALGEBRY A.2 Maticový kalkulus —> R funkce, xq g Rra je reálný vektor. Definujeme gradient funkce / v bodě xq Buď / : R" jako vektor Buď F : Rmxra matice funkce, Xq g Rmxra je reálná matice. Gradient funkce F v bodě Xq je dF dX (*>) = dF dX. (X0 ij 1,3 Některé užitečné gradienty: d_ dx d_ dx (yTx) = y (xTy) = y d_ dx d_ dx d_ dx (yTAx) = Ay (xTAy) = Ay (x1 Ax) = (A + A1 )x _d_ dÄ tľ(A) =1 4jtr (BAD) = BTDT dA v ' ^-tr (ABAT) = 2AB dA v ' d_ dÄ \BAD\ = \BAD\A-T (A.10) (A.ll) (A.12) (A.13) (A.14) (A.15) (A.16) (A.17) (A.18) kde jsme položili A T = (A lA Buď g : funkce, xq g Rra je reálný vektor. Jakobián funkce g v bodě xq je matice dg_ dx Některé užitečné Jakobiány: (xq) d dgi dx-. (xq) i,3 dx (Ax) = A (A.19) Stopa matice je součet jejích diagonálních prvků. Značíme tr A Platí tr (A + B) = tr A + tr B yTx = tr (yxT) (A.20) (A.21) Dodatek B Pravděpodobnost B.l Základní pojmy Označme (Q,A,V) jako pravděpodobnostní prostor, kde Q značí základní prostor, A jevové pole na Q a P je pravděpodobnostní míra (provděpodobnost) na A. Náhodná veličina (náhodný vektor) je zobrazení X : Q —> Rra. Distribuční funkce náhodného vektoru X G Rra je zobrazení F : Rra —> R definované vztahem F{x) = P{{uj G Q : X(w) < x}) Hustota pravděpodobnosti náhodného vektoru X G Rra je zobrazení F : Rra —> R a je definována jako 9raF V{X) = dXl---dxn{x) Hustota pravděpodobnosti má smysl pouze pro diferencovatelné distribuční funkce (spojité náhodné vektory). Tam, kde nemůže dojít k nedorozumění, budeme všechny distribuční funkce označovat symbolem F a všechny hustoty symbolem p. p{x) bude tedy hustota příslušná náhodné veličině X a p{y) bude hustota příslušná náhodné veličině Y. Spojení Z = (XT,YT) je náhodný vektor. Známe-li hustotu p(z), pak hustoty p{x) a p(y) získáme marginalizací p(x) = p(z) dy (B-l) p{y) = j p(z) dx (B-2) B.2 Číselné charakteristiky (momenty) Střední hodnota náhodného vektoru x (první moment) je vektor, definovaný jako E(x)= x dF{x) it» Druhý moment náhodného vektoru x je matice m2(x) = / xxT dF(x) Jw,n 65 66 DODATEK B. PRAVDĚPODOBNOST Rozptyl náhodného vektoru x je definován jako var (a;) = E (x — E (a;))(a; — E(x))T Platí, že var (a?) = m,2(x) — [E (x)]2. Zápisem X G Cn rozumíme, že X : Q —> Rra je náhodný vektor dimenze n s konečným druhým momentem. Takový vektor má i konečný první moment. {^il^l-oo Je spočetná posloupnost náhodných vektorů (náhodný proces). B.3 Normální rozdělení Náhodný vektor X G Cn má normální rozdělení, právě když hustota pravděpodobnosti je dána vzorcem p(x) = , l e|(x-/i)Ts-1(x-M) (B_3) kde ß G W1 je reálný vektor a S G Wixn je symetrická pozitivně definitní matice. Pak píšeme X~N(fi,Y,). Lze ukázat, že pro parametry /x, X platí EI = /i (B.4) varX = S (B.5) Normálně rozdělené náhodné vektory X &Y jsou nezávislé právě thedy, jsou-li nekorelované, tj. cov(X, Y) = 0. Mějme náhodný vektor X G Cn, X ~ N(/j,x,Y,x), vektor a G Mm a matici 13 G RmXí\ Pak náhodný vektor Y = a + 13 X má normální rozdělení a platí F ~ N (a + 5//x, BT,XBT). Zkráceně X - JVQíj., ^)^.a + M~ JV(a + Bßx, BZXBT) (B.6) B.4 Podmíněná hustota pravděpodobnosti Podmíněná hustota pravděpodobnosti p{x\y) je definována následujícím předpisem Pixly) = w ( 7) Na tuto funkci nahlížíme jako na funkci proměnné x při fixované hodnotě y. Je to hustota pravděpodobnosti náhodné veličiny x při dané hodnotě náhodné veličiny y, přičemž závislost obou náhodných veličin je vyjádřena simultánní hustotou pravděpodobnosti p{x,y). Jako hustota pravděpodobnosti splňuje podmínku p(x\y) dx = 1 (B-8) Z definice podměné hustoty se jednoduše odvodí řetězové pravidlo p(x, y) = p(x\y)p(y) = p(y\x)p(x) (B.9) B.5. PODMÍNĚNÉ ČÍSELNÉ CHARAKTERISTIKY 67 a také často používaný Bayesův vzorec p(y\x) =, Pixly)piy) p(x\y)p(y) p{x) Jp(x\y)p(y) dy Jsou-li náhodné vektory X a Y nezávislé, pak platí p(x\y) =p(x) (B.10) (B.ll) B.5 Podmíněné číselné charakteristiky Pro náhodnou veličinu X\z = (X\Z = z) zavádíme podmíněnou střední hodnotu E{X\z) jako E(X|z) xp{x\z) dx z pravidel pro počítání s podmíněnými hustotami vyplývají ekvivalentní vztahy E(X|Z)- /" ~.p(aľ'*) E(X|Z) E(X|Z) p{z) Jxp(x,z) dx J p(x,z) dx J xp{x\z)p{x) dx Jp{x\z)p{x) dx (B.12) (B.13) (B.14) (B.15) B.6 Lineární regrese Nechť j/,ee£l jsou náhodné veličiny, x\,..., xp jsou reálné proměnné, b\,..., bp jsou reálné konstanty (parametry, regresní koeficienty). Lineární regresní model je rovnice Y = b\X\ + ... + bpxp + e (B.16) popisující lineární závislost výstupu Y na regresorech Xj. Člen e je náhodná chyba, nebo tzv. šum měření s vlastnostmi E (e) = 0 var (e) = a2 kde a2 je známá konstanta. Lineární regresní model je přístupný pozorování. Provedeme n měření (pozorování) a dostaneme vektor výstupů y a matici plánu X: y = Ví Vr, X = xn %nl X\p X, np Při měření předpokládáme, že šumy e^ tvoří posloupnost stejně rozdělených a stochasticky nezávislých náhodných veličin. 68 DODATEK B. PRAVDĚPODOBNOST Na základě dat X, y máme odhadnout neznámé regresní koeficienty b = (b\,..., bp) tak aby výraz n p n ?(6) = D»* - E 6i*y)2 = Dw - y^2 (B-17) í=l J = l i=\ byl minimální. Funkce i - E hxij)2 = 5>i - y*)2 (B.19) " p i=i i=i i=\ se nazývá reziduálni rozptyl a je nestranným odhadem hodnoty