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 Lotkův-Volterrů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í estimátor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 OBSAH 7 TEORIE CHAOSU 61 A Přehled maticové algebry 63 A.1 Blokové matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 A.2 Maticový kalkulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 B Pravděpodobnost 65 B.1 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 souvis- lostech. 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: Necht' S je separabilní metrický prostor a T R. Dynamický systém je pak množina transformací : T × T × S S, které splňují: a) Pro všechna t0, t1, t2 T, x S platí (t2, t0, x) = (t2, t1, (t1, t0, x)) b) Pro všechna t, t0 T a x S a pro všechny posloupnosti {tn}, {xn}, pro které platí tn t a xn x platí (tn, t0, xn) (t, t0, x). Množinu S nazveme stavový prostor a libovolný bod x S pak nazýváme stavem systému. Obraz (t, t0, x) představuje stav systému v čase t, byl-li systém v čase t0 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 vek- torů 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 Rm, takzvaný řídící vektor. Podle časové závislosti vlastností systému Máme časově neměnné systémy (časově in- variantní nebo též autonomní systémy), které splňují (t, t0, x) = (t + t, t0 + t, 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 = Rn. 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 f : Rn Rn je hladká funkce. Při splnění tzv. Cauchyho počáteční podmínky x(0) = 0 ˙x(0) = 1 x(m-1) (0) = m-1 (2.2) je řešením systému jednoznačně určená funkce x : R Rn. 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-1)) = (y0, y1, . . . , ym-1). Zřejmě platí ˙y = ( ˙x, ¨x, . . . , x(m)). Přitom ˙x = y1 ¨x = y2 x(m-1) = ym-1 x(m) = f(y0, y1, . . . , ym-1) (2.3) Počáteční podmínky můžeme přepsat jako y0(0) = 0 y1(0) = 1 ym-1(0) = m-1 (2.4) což můžeme zkráceně napsat jako ˙y = F(y) y(0) = (2.5) což je systém diferenciálních rovnic prvního řádu s Cauchyho počáteční podmínkou. Vztah mezi x a y je bijekce, takže oba systémy jsou ekvivalentní. Nadále se tedy budeme zabývat systémy prvního řádu ˙x = f(x) (2.6) kde ˙x = dx1 dt , . . . , dxn dt T a f je funkce f : Rn Rn. Priklad 2.2. Uvažujme skalární systém druhého řádu ve tvaru ¨y = a ˙y + by + c s Cauchyho počáteční podmínkou y(0) = y0 ˙y(0) = y1 Položme x = (x1, x2) = (y, ˙y). Potom platí: ˙x1 = x2 (2.7a) ˙x2 = ax2 + bx1 + c (2.7b) s počáteční podmínkou x1(0) = y0 x2(0) = y1. 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 f : Rn Rn definuje na Rn 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í f. Byl-li systém v čase t0 ve stavu x0, 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, t0, x) Operátor toku je také řešením soustavy (2.6), nebot' v každém bodě t 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(t0) = x0 lze vyjádřit pomocí evolučního operátoru takto: x(t) = t-t0 (x0) Takovou křivku x(t) nazýváme trajektorie systému. Z definice evolučního operátoru vyplývají následující vztahy: s+t(x) = s(t(x)) (2.8) t(-t(x)) = 0(x) = x (2.9) -1 t = -t (2.10) Definice 2.3. Existuje-li x S takové, že t(x) = 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(x) pro každé x 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í, monotonně nebo period- icky, 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(p) 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 = kE(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) = + p a nabídky S(p) = + p. Počátek přizpůsobovacího procesu je p(0) = p0. Upravíme stavovou rovnici: ˙p = k( + p - - p) = k( - )p + k( - ) Stavová rovnice je lineární diferenciální rovnice prvního řádu s konstantními koeficienty. Řešení hledáme ve tvaru p(t) = a + bect Dosazením do stavové rovnice a počáteční podmínky vychází a = - - (2.11) b = p0 - - - (2.12) c = k( - ) (2.13) takže řešením je křivka p(t) = - - + p0 - - - ek(-)t 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ě p() = pe = - - při které je exces poptávky nulový, tzn. D(p()) = S(p()). 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. ˙r = 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 I klesají s rostoucím úrokem, tj. I = -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 ˙Y = h(D - Y ) = h(C + I - Y ) = h(c - 1)Y - ahr = -hsY - ahr ˙r = mkY - m M Derivujeme-li první rovnici a do vzniklého výrazu dosadíme rovnici druhou, dostaneme ¨Y = -hs ˙Y - ah ˙r = -hs ˙Y - ahmkY + ahm M neboli ¨Y + hs ˙Y + ahmkY = ahm M což je lineární ODR druhého řádu. Napišme její charakteristickou rovnici 2 + hs + ahmk = 0 Zde se řešení rozpadá na tři části a) hs2 > 4amk V tomto případě existují 2 reálná řešení charakteristické rovnice 1 = 1 2 -hs + h2s2 - 4ahmk 2 = 1 2 -hs - h2s2 - 4ahmk kterým odpovídá řešení DR ve tvaru Y (t) = M k + C1e1t + C2e2t Toto řešení bude stabilní, jestliže bude platit 1 < 0, 2 < 0. Ukážeme, že tato podmínka je splněna Vzhledem k podmínce (viii) platí ahmk > 0. Pak také -4ahmk < 0 h2 s2 - 4ahmk < h2 s2 h2s2 - 4ahmk < |hs| h2s2 - 4ahmk < hs 1 2 -hs + h2s2 - 4ahmk < 0 1 < 0 a podobně - h2s2 - 4ahmk < hs 1 2 -hs - h2s2 - 4ahmk < 0 2 < 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 = 2 = = - 1 2 hs < 0 kterému odpovídá řešení Y (t) = M k + C1et + C2tet 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 1 = 1 2 -hs + i 4ahmk - h2s2 2 = 1 2 -hs - i 4ahmk - h2s2 Položíme-li = -1 2 hs a = 1 2 4ahmk - h2s2, pak můžeme řešení napsat ve tvaru Y (t) = M k + et (C1 cos t + C2 sin t) (2.14) Protože < 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 hos- podářských cyklů s periodou 4/ 4ahmk - 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) s, (0, 1). 2.2. EKONOMICKÉ APLIKACE 13 Spočtěme ˙k: ˙k = ˙KL - K ˙L L2 = ˙K L - K L ˙L L z toho vyplývá, že ˙k k = L K ˙k = ˙K K - ˙L 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 - K dosazením dostáváme ˙k k = s Y K - ( + n) = s L K f(k) - ( + n) = s k f(k) - ( + n) Vynásobením k máme ˙k = sf(k) - ( + n)k = sf(k) - k kde f je rostoucí konkávní funkce, splňující Ianadovy podmínky limk0 f (k) = a limk f (k) = 0. V případě Cobb-Douglasovy produkční funkce Y = K L1- dostáváme Y L = K L- = k = f(k) a diferenciální rovnice růstu je ve tvaru ˙k = -k + sk což je Bernoulliho rovnice. Provedeme substitucí x = k1-. Pak ˙x = (1 - )k- ˙k = ˙k = ˙x k 1 - k = k1- k = xk Dosazením do DR máme ˙x k 1 - = -xk + sk odkud dostáváme nakonec ˙x = -(1 - )x + (1 - )s což je rovnice lineární. Jejím řešením je trajektorie x(t) = x0 - s e-(1-)t + s nebo ekvivalentně trajektorie k1- = k1- 0 - s e-(1-)t + s protože platí 1 - > 0, = + n > 0, je model stabilní s pevným bodem k = s + n 1 1- 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 fluktu- ace 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(t) roste přirozeným tempem a (tj. ˙x = ax). Jejich počet je snižován o cxy díky přítomnosti dravců y(t), 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) ˙y = dxy - by (2.16) Položíme-li ˙x = 0, ˙y = 0, obdržíme dva kritické body systému. z = (0, 0)T je tzv. rovnováha po vyhynutí (extinction equilibrium) a z = (b d, a c ) T je tzv. rovnováha koexistence (coexistence equilibrium). Numerické řešení DR (2.15) je v Matlabu ve funkci odedemo. Kapitola 3 DISKRÉTNÍ 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 : k Z}. 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 f(t) lze získat diskrétní funkci g(k) vzorkováním. Přitom platí f(t) = 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 - 1), . . . , y(k - ly)) (3.1) kde y Rp je výstup systému. Tato rovnice ukazuje závislost výstupu y(k) na minulých výstupech y(k - i). Řád systému ly pak udává pamět' systému. Přidáme-li do systému prvek řízení, tj. řídící vektor u Rm, dostaneme rovnici y(k) = f(y(k - 1), . . . , y(k - ly), u(k - 1), . . . , u(k - lu)) (3.2) Takto napsaná rovnice naznačuje, že řízení u(k-1) 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 + 1), x(k + 2), . . . a výstupy y(k), y(k + 1), y(k + 2), . . .. Definice 3.1. Dynamický systém popsaný rovnicemi x(t + 1) = 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 rovnicí (rovnice přechodu, state equation, plant equa- tion, transition equation) a rovnice (3.3b) je zvána rovnicí výstupu (rovnice měření, observa- tion 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 - 1)T , . . . , u(k - lu + 1)T ) T pak systém (3.2) lze napsat jako systém ve stavovém tvaru. Důkaz. Původní rovnici y(k) = f(y(k - 1), . . . , y(k - ly), u(k - 1), . . . , u(k - lu)) lze substitucí t + 1 = k napsat ve tvaru y(t + 1) = f(y(t), . . . , y(k - ly + 1), u(k - 1), . . . , u(k - lu + 1)) = f(x(t)) Zřejmě platí x(t + 1) = (y(t + 1)T , . . . , y(t - ly + 2)T , u(t)T , . . . , u(t - lu + 2)T ) T 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(t) a u(t). y(t + 1) = f(x(t)) y(t) = A1x(t) ... ... y(t - ly + 2) = Aly x(t) u(t) = u(t) u(t - l) = B1x(t) ... ... u(t - lu + 1) = Blu-1x(t) kde A1, . . . , Aly , B1, . . . , Blu-1 jsou permutační matice, vybírající z vektoru x(t) vždy jen příslušnou část. Celkem lze tedy soustavu napsat jako x(t + 1) = F(x(t), u(t)) což je popis systému ve stavovém tvaru 3.1. TEORIE 17 Priklad 3.3. Uvažujme lineární dynamický systém ve tvaru yt = a1yt-1 + a2yt-2 + b1ut-1 + b2ut-2 Posuneme čas o jeden krok dopředu a máme yt+1 = a1yt + a2yt-1 + b1ut + b2ut-1 Nyní položíme xt = (x1t, x2t, x3t)T = (yt, yt-1, ut-1)T . Potom platí xt+1 = (yt+1, yt, ut)T a můžeme psát x1,t+1 = a1x1,t + a2x2,t + b2x3,t + b1ut x2,t+1 = x1,t x3,t+1 = ut yt = x1,t což je systém ve stavovém tvaru (3.3). Definice 3.4. Dynamický systém popsaný rovnicemi x(t + 1) = f(x(t), u(t), v(t)) (3.4a) y(t) = g(x(t), u(t), w(t)) (3.4b) kde v(t) Lnv, w(t) Lnw jsou náhodné vektory, splňující podmínku E v(t) = 0, E w(t) = 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 + 1) = f(x(t), u(t)) + v(t) y(t) = g(x(t), u(t)) + w(t) 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 Lnx náhodných vektorů. Stejně tak výstupy y(t) 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. Necht' A Rnx×nx , B Rnx×nu , C Rny×nx , D Rny×nu . 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(t) 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) EvtwT t = 0 (iii) EvtvT t = Q (iv) EwtwT t = 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ý determin- istický 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)-1 Bu (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(t) = x(t) - xe, pak stavová rovnice přejde v x(t + 1) - xe = Ax(t) + Bu - xe x(t + 1) - xe = Ax(t) - Axe + Bu + Axe - xe x(t + 1) - xe = A(x(t) - xe) + Bu - (I - A)xe y(t + 1) = Ay(t) + Bu - (I - A)(I - A)-1 Bu y(t + 1) = Ay(t) Má-li x(t) konvergovat k xe, pak musí y(t) konvergovat k nule. Vyjdeme-li z počáteční podmínky y(0) = y0, pak stav y(t) se dá explicitně spočítat jako y(t) = At y0 Je vidět, že má-li y(n) pro libovolné y0 konvergovat k nule, musí An konvergovat k nulové matici. Pro každou čtvercovou matici A existuje regulární matice P taková, že P-1 AP = D kde D je diagonální matice, která má v diagonále vlastní čísla matice A D = diag(1, . . . , nx ) 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-1 PDP-1 = PD2 P-1 a tak dále až matici An zapíšeme jako An = PDn P-1 Matice Dn je opět diagonální a obsahuje prvky Dn = diag(n 1 , . . . , n nx ) Posloupnost An konverguje k nule právě tehdy, když posloupnost Dn konverguje k nule a to je právě tehdy, když všechny posloupnosti n i konvergují k nule. To je splněno, když pro všechna i {1, . . . , nx} platí |i| < 1 (3.9) Protože vlastní čísla i 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. Chování stabilního stochastického systému je významně ovlivněno vlastnostmi stocha- stické složky. Bude-li šum slabý, bude systém po dostatečně dlouhém čase skutečně vyka- zovat 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 Samuelsonův hospodářský cyklus (1939) Samuelsonův model je postaven na multiplikátoru a akcelerátoru. Struktura modelu je popsána těmito rovnicemi: Ct = Yt-1 It = v(Ct - Ct-1) Gt = 1 Yt = Ct + It + Gt Předpokládá se 0 < c < 1, v > 0. 1/1 - c je multiplikátor a v je akcelerátor. Substitucí dostaneme typickou diferenční rovnici 2. řádu Yt - c(1 + v)Yt-1 + cvYt-2 = 1 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í y0 nehomogenní rovnice (tzv. partikulární řešení) a bázi y1, . . . , yn všech řešení homegenní rovnice. Libovolné řešení nehomogenní rovnice pak dostaneme jako y = y0 + a1y1 + . . . + anyn Lehce se ukáže, že funkce y0(t) = 1 1 - c je partikulárním řešením naší diferenční rovnice. Hledejme nyní řešení homegenní rovnice ve tvaru Yt = t. Pak musí platit t - c(1 + v)t-1 + cvt-2 = 0 Vydělením rovnice výrazem t-2 dostaneme 2 - c(1 + v) + cv = 0 Z této rovnice dostaneme jako kořen kvadratické rovnice. Řešení závisí na znaménku diskriminantu = c2 (1 + v)2 - 4cv a) c > 4v/(1 + v)2 V tomto případě existují 2 reálná řešení, která tvoří bázi řešení homogenní rovnice. 1 = 1 2 c(1 + v) + c2(1 + v)2 - 4cv 2 = 1 2 c(1 + v) - c2(1 + v)2 - 4cv kterým odpovídá obecné řešení ve tvaru Yt = 1 1 - c + a1t 1 + a2t 2 Toto řešení bude stabilní, jestliže bude platit |1| < 1, |2| < 1. |1| < 1 |2| < 1 |12| < 1 1 4 c(1 + v) + c2(1 + v)2 - 4cv c(1 + v) - c2(1 + v)2 - 4cv < 1 1 4 4cv < 1 cv < 1 3.2. EKONOMICKÉ APLIKACE 21 b) c = 4v/(1 + v)2 Za této situace má kvadratická rovnice dvojnásobný reálný kořen 1 = 2 = = c(1 + v) 2 kterému odpovídá řešení Yt = 1 1 - c + a1t + a2tt Toto řešení je stabilní za podmínky c(1 + v) < 2 v < 1. c) c < 4v/(1 + v)2 V tomto případě existují 2 komplexní kořeny charakteristické rovnice 1 = 1 2 c(1 + v) + i 4cv - c2(1 + v)2 2 = 1 2 c(1 + v) - i 4cv - c2(1 + v)2 Báze řešení homogenní rovnice je tvořena funkcemi u = t 1 a v = t 2, což jsou komplexní funkce. Ukážeme, že lze najít reálné funkce, které tvoří bázi řešení. Zaved'me goniometrické vyjádření komplexních kořenů 1 = r(cos + i sin ) 2 = r(cos - i sin ) Pak platí a1t 1 + a2t 2 = = a1rt (cos t + i sin t) + a2rt (cos t - i sin t) = = rt ((a1 + a2) cos t + i(a1 - a2) sin t) = = rt (b1 cos t + b2 sin t) přitom r je absolutní hodnota společná komplexním číslům 1, 2, která je rovna 1 1 = 12 = cv Systém bude stabilní za podmínky r < 1 cv < 1. Samuelsonův 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ů - v a c. Každému typu chování přináleží v parametrickém prostoru (v, 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, ) × (0, 1). 1. c > 4v/(1 + v)2 cv < 1: monotonní konvergence - stabilní 2. c > 4v/(1 + v)2 cv = 1: neutrální stabilita 3. c > 4v/(1 + v)2 cv > 1: monotonní divergence - nestabilní 4. c = 4v/(1 + v)2 v < 1: jeden kriticky tlumený kmit - stabilní 22 KAPITOLA 3. DISKRÉTNÍ DYNAMICKÉ SYSTÉMY 5. c = 4v/(1 + v)2 v = 1: lineární divergence - nestabilní 6. c = 4v/(1 + v)2 v > 1: monotonní divergence - nestabilní 7. c < 4v/(1 + v)2 cv < 1: periodická konvergence - stabilní 8. c < 4v/(1 + v)2 cv = 1: konstantní periodická oscilace - stabilní 9. c < 4v/(1 + v)2 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 = (1 - s)Yt-1 It = A0(1 + g)t + v(Yt-1 - Yt-2) Yt = Ct + It Spotřeba Ct je zpožděná lineární funkce Yt-1. Parametr s splňující podmínku 0 < s < 1 je mezní sklon k úsporám. Investice It mají dvě komponenty: autonomní investice A0(1 + g)t rostoucí exponenciálně konstantním tempem g a indukované investice v(Yt-1 - 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(1 + g)t Autonomní investice můžeme chápat jako exogenní veličinu, kterou zde pro jednoduchost ztotožníme s řízením ut, tj. Yt = (1 - s + v)Yt-1 + vYt-2 + ut-1 Posuneme čas o jeden krok dopředu a máme Yt+1 = (1 - s + v)Yt + vYt-1 + ut Nyní položíme xt = (x1t, x2t)T = (Yt, Yt-1)T . Potom platí xt+1 = (Yt+1, Yt)T odtud pak x1,t+1 = (1 - s + v)x1,t + vx2,t + ut x2,t+1 = x1,t (3.10) a také Yt = x1,t (3.11) rovnice (3.10) je stavovou rovnicí a rovnice (3.11) je výstupní rovnicí uvažovaného dynam- ické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. Bud' dán spojitý systém ve stavovém tvaru (2.6) s počáteční podmínkou. Dále necht' T = {tk : k {1, . . . , N}, N N} je tzv. časový horizont. Simulovat daný systém na horizontu T pak znamená nalézt vzorky trajektorie systému v časech tk. Jinými slovy hledáme posloupnost {x0, . . . , xN } takovou, že xk = x(tk), kde funkce x(t) je řešením soustavy DR (2.6). Definice 4.2. Bud' dán řízený stochastický diskrétní systém (3.3) s počáteční podmínkou. Dále necht' T = {1, . . . , N}, N N je časový horizont. Bud' {u(0), . . . , u(N)} známá posloup- nost řídících vektorů--trajektorie řízení 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 standartně 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 = c0 + c1Yt + c2Ct-1 (4.1a) It = i0 + i1Yt + i2(Yt-1 - Yt-2) + i3Rt-4 (4.1b) Rt = r0 + r1Yt + r2(Yt - Yt-1) + r3(Mt - Mt-1) + r4(Rt-1 + Rt-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 I 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 hod- noty jsou: c0 = -9.4540 c1 = 0.05410 c2 = 0.9260 i0 = -66.1950 i1 = 0.16840 i2 = 0.2181 i3 = -11.2650 r0 = -0.5561 r1 = 0.00051 r2 = 0.0135 r3 = -0.0853 r4 = 0.4259 Pro každou rovnici byly vypočteny reziduální rozptyly. Jejich hodnoty jsou: s2 1 = 11.2302 pro rovnici (4.1a) s2 2 = 23.8642 pro rovnici (4.1b) s2 3 = 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ě yt = 0yt + 1yt-1 + 2yt-2 + 3yt-3 + 4yt-4 + 0ut (4.2) 4.2. P ŘÍKLADY 25 kde jsme položili yt = Ct It Rt Yt ut = Mt - Mt-1 Gt 1 0 = 0 0 0 c1 0 0 0 i1 0 0 0 r1 + r2 1 1 0 0 1 = c2 0 0 0 0 0 0 i1 0 0 r4 -r2 0 0 0 0 2 = 0 0 0 0 0 0 0 -i1 0 0 r4 0 0 0 0 0 3 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 = 0 0 0 0 0 0 i3 0 0 0 0 0 0 0 0 0 0 = 0 0 c0 0 0 i0 r3 0 r0 0 1 0 Poslední element vektoru ut je jednička. Tato úprava formálně zbavuje soustavu úrovňových konstant c0, i0, r0. 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 - 0yt = 1yt-1 + 2yt-2 + 3yt-3 + 4yt-4 + 0ut (I - 0)yt = 1yt-1 + 2yt-2 + 3yt-3 + 4yt-4 + 0ut yt = (I - 0)-1 1yt-1 + (I - 0)-1 2yt-2 + (I - 0)-1 3yt-3+ + (I - 0)-1 4yt-4 + (I - 0)-1 0ut yt = a1yt-1 + a2yt-2 + a3yt-3 + a4yt-4 + b0ut Tím jsme dostali model ARX, který zapíšeme blokově takto yt yt-1 yt-2 yt-3 = a1 a2 a3 a4 I 0 0 0 0 I 0 0 0 0 I 0 yt-1 yt-2 yt-3 yt-4 + b0 0 0 0 ut (4.3) 26 KAPITOLA 4. SIMULACE DYNAMICK ÝCH SYSTÉMŮ kde a1 = (I - 0)-1 1 = 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 a2 = (I - 0)-1 2 = 0 0 0 -0.0152 0 0 0 -0.2653 0 0 0.4259 -0.0039 0 0 0 -0.2805 a3 = (I - 0)-1 3 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a4 = (I - 0)-1 4 = 0 0 -0.7838 0 0 0 -13.7049 0 0 0 -0.2030 0 0 0 -14.4887 0 b0 = (I - 0)-1 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-1 + But (4.4) Rovnice (4.4) je stavovou rovnicí našeho systému. Přitom x R16, u R3. Rozptyly s2 1, s2 2, s2 3 jsou nestrannými odhady rozptylů stochastické složky původního mod- elu (4.1). Tyto složky předpokládáme nezávislé a normální. Označme e šum původního modelu e N 0 0 0 , 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í E v = (I - 0)-1 E e = 0 var v = (I - 0)-1 var e (I - 0)-1 T Matice var v 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: = -0.4905 + 0.4354i -0.4905 - 0.4354i 0.9878 0.9026 0.3937 + 0.5633i 0.3937 - 0.5633i 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(x0,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í estimátor 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á infor- mace 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ě kompliko- vaný 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 Lnx a vektor naměřených dat z Rnz , který je realizací náhodného vektoru Z Lnz. 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 p0(x) = p(x, z) dz 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 p1(x) = p(x|y) kterou nazýváme aposteriorní 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. Bud' E Rnx množina přípustných estimátorů náhodného vektoru X. Pro každý přípustný estimátor y E zaved'me střední kvadratickou chybu vzorcem J(y) = E ((X - y)T (X - y)) Estimátor ^X E nazveme optimálním podle kriteria nejmenších čtverců (least mean square), jestliže platí J( ^X) = min yE J(y) ^X se také nazývá optimální mean square estimátor (MS estimátor). Veta 5.2 (optimální apriorní estimátor). Střední hodnota E X 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 XT X - yT E X - E XT y + yT y = = E XT X - 2yT E X + yT y 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 J y = -2E X + 2y = 0 nebo také ^X = E X (5.1) 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. Zaved'me chybu estimace ~X Lnx jako rozdíl ~X = X - ^X Pak platí E ~X = E X - ^X = E X - E X = 0, takže střední chyba estimace je rovna nule. Říkáme, že odhad ^X je nestranný (nevychýlený, unbiased). Ekvivalentně také platí E ~X = 0 E ^X = E X 5.1. OPTIM ÁLNÍ ESTIM ÁTOR 31 Jako míru přesnosti naší estimace spočtěme rozptyl chyby var ~X = E ( ~X - E ~X) T ( ~X - E ~X) = = E ~XT ~X = = E (X - ^X) T (X - ^X) = = E (X - E X)T (X - E X) = = var X Bude-li rozptyl malý, je náš odhad dobrý. Inverzní matici (var ~X)-1 nazýváme informační maticí. Veta 5.3 (optimální aposteriorní estimátor). 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í estimátor. Přitom známe hustotu p(x, z). Optimální estimátor bude deterministickou funkcí naměřeného vektoru z, tj. ^X : Rnz R. Napišme funkcionál J J = (x - ^X(z)) T (x - ^X(z))p(x, z) dx dz kde integrujeme přes celý prostor Rnx , případně Rnz . 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)) T (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 min- imalizujeme pro každé pevné z vnitřní integrál H = (x - ^X(z)) T (x - ^X(z))p(x|z) dx Integrál H můžeme napsat jako H = E ((X - ^X(z)) T (X - ^X(z))|z) = = E (XT X|z) - ^X(z) T E (X|z) - E (XT |z) ^X(z) + ^XT (z) ^X(z) = = E (XT X|z) - 2 ^X(z) T E (X|z) + ^X(z) T ^X(z) Stejným způsobem jako minule dostáváme H ^X = -2E (XT X|Z) + 2 ^X(Z) = 0 což dává ^X(z) = E (X|z) (5.2) 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í estimátor je nestranný, což dokazuje výpočet E ~X = 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 var ~X = var [(X - ^X)|Z] = = E [(X - E X|Z)(X - E X|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 var ~X = 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í E (X|z) = x p(x)p(z) p(z) dx = = xp(x) dx = = E X Při nezávislosti tedy aposteriorní a apriorní estimátor 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í). Necht' platí X Z N x z , x xz zx z Pak aposteriorní náhodná veličina X|z má normální rozložení se střední hodnotou E (X|z) = x + xz-1 z (z - z) (5.4) a varianční maticí var (X|z) = x - xz-1 z zx (5.5) Důkaz. Definujme náhodný vektor Y Lnx + nz jako Y = x z Hustota pravděpodobnosti náhodného vektoru Y je dána vzorcem p(y) = 1 (2)nx+nz |y| e- 1 2 (y-y)T -1 y (y-y) 5.1. OPTIM ÁLNÍ ESTIM ÁTOR 33 kde y = x z y = x xz zx z Podle (A.5) je -1 y = D-1 -D-1xz-1 z --1 z zxD-1 -1 z + -1 z zxD-1xz-1 z kde Schurův komplement k x je D = x - xz-1 z zx. Podle (B.13) je p(x|z) = p(y) p(z) neboli p(x|z) = (2)nz |z| (2)nx+nz |y| e- 1 2 [(y-y)T -1 y (y-y)-(z-z)T -1 z (z-z)] = = 1 (2)nx |y| / |z| e- 1 2 [(y-y)T -1 y (y-y)-(z-z)T -1 z (z-z)] Napišme si exponent a upravujme ho (y - y)T -1 y (y - y) - (z - z)T -1 z (z - z) = = x - x z - z T -1 y x - x z - z - (z - z)T -1 z (z - z) = = (x - x)T D-1 (x - x) - (x - x)T D-1 xz-1 z (z - z)- - (z - z)T -1 z zxD-1 (x - x) + (z - z)T (-1 z + -1 z zxD-1 xz-1 z )(z - z)- - (z - z)T -1 z (z - z) = = (x - x)T D-1 (x - x) - (x - x)T D-1 xz-1 z (z - z)- - (z - z)T -1 z zxD-1 (x - x) + (z - z)T -1 z zxD-1 xz-1 z (z - z) = = (x - x)T D-1 (x - x) - (x - x)T D-1 [xz-1 z (z - z)]- - [xz-1 z (z - z)] T D-1 (x - x) + [xz-1 z (z - z)] T D-1 [xz-1 z (z - z)] Z prvních dvou členů vzniklého výrazu vytkneme (x - x)T D-1 a z dalších dvou členů vytkneme [xz-1 z (z - z)] T D-1 a dostaneme (x - x)T D-1 (x - x) - xz-1 z (z - z) -[xz-1 z (z - z)] T D-1 (x - x) - xz-1 z (z - z) Nyní můžeme vytknout zprava výraz (x - x) - xz-1 z (z - z) , čímž dostaneme (x - x) - xz-1 z (z - z) T D-1 (x - x) - xz-1 z (z - z) = = x - (x + xz-1 z (z - z)) T D-1 x - (x + xz-1 z (z - z)) Podle (A.9) platí, že |y| = |z| |D| = |y| |z| = |D| 34 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU Položíme-li x|z = x + xz-1 z (z - z), dostáváme nakonec, že p(x|z) = 1 (2)nx |D| e- 1 2 (x-x|z)T D-1(x-x|z) tedy podmíněný náhodný vektor X|z má normální rozložení se střední hodnotou x|z = x + xz-1 z (z - z) a rozptylem D = x - xz-1 z zx což bylo potřeba dokázat 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) var X. 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 E X. 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í mat- ice 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í. Vrat'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 estimátoru může být v některých případech velmi obtížný. Aby- chom tento výpočet výrazně zjednodušili, omezme množinu přípustných estimátorů jen na ty estimátory, které jsou lineární funkcí naměřeného vektoru z. Protože třída přípustných estimátorů je nyní omezena, lineární MS estimátor nebude obecně tak dobrý jako nelineární MS estimátor. Veta 5.5 (lineární MS estimátor). Optimální MS estimátor na přípustné množině lineárních estimátorů E = {y = Az + b : A Rnx×nz , b Rnx } je estimátor ^XLMS = x + xz-1 z (z - z) (5.7) kde x = E X, z = E Z, xz = cov(X, Z), 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 - b)T (X - AZ - b) což můžeme vzhledem k (A.21) napsat jako J = E (tr (X - AZ - b)(X - AZ - b)T ) = = tr E (X - AZ - b)(X - AZ - b)T = = tr E [(X - x) - (AZ + b - x)][(X - x) - (AZ + b - x)]T = = tr [E (X - x)(X - x)T - E (X - x)(AZ + b - x)T - - E (AZ + b - x)(X - x)T + E (AZ + b - x)(AZ + b - x)T ] 5.1. OPTIM ÁLNÍ ESTIM ÁTOR 35 Označme x = var X = E (X - x)(X - x)T . Předpokládejme navíc, že matice A je symet- rická. Potom platí J = tr [x - xzAT - Azx + E (AZ + b - x)(AZ + b - x)T ] = = tr [x - 2Azx + E (A(Z - z) + (Az + b - x))(A(Z - z) + (Az + b - x))T ] = = tr [x - 2Azx + AzAT + (Az + b - x)(Az + b - x)T ] = = tr [x + A(z + zz T )AT + (b - x)(b - x)T + 2Az(b - x)T - 2Azx] Podle (A.16) a (A.17) platí A tr (A(z + zz T )AT )=2A(z + zz T ) A tr (2Azx) =2xz A tr (2Az(b - x)T ) =2(b - x)z T Podle (A.21) a (A.10) pak máme b tr ((b - x)(b - x)T )=2(b - x) b tr (2Az(b - x)T ) =2Az Můžeme tedy napsat parciální derivace J A =2A(z + zz T ) - 2xz + 2(b - x)z T = 0 (5.8) J b =2(b - x) + 2Az = 0 (5.9) Z rovnice (5.9) vypočteme b = x - Az a dosazením do (5.8) dostaneme Az - xz = 0 což dává A = xz-1 z Optimální lineární MS estimátor je potom dán vzorcem ^XLMS = x + xz-1 z (z - z) což je (5.7) Prozkoumejme strukturu vzorce (5.7). Definujeme-li výstupní reziduum jako ~Z = Z - E Z (5.10) pak vzorec (5.7) přejde v ^XLMS = x + xz-1 z ~z. Člen x představuje apriorní estimátor pro X a druhý člen je korekce, založená na výstupním reziduu. Bude-li závislost mezi X a Z 36 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU velká, bude apriorní odhad silněji korigován. Naopak bude-li veliká nejistota ohledně Z (velký rozptyl z), bude korekční člen potlačen. Konečně změříme-li z = E Z, pak se v podstatě náš apriorní odhad potvrdil a k žádné korekci nedojde. Srovnáme-li nelineární odhad s lineárním, pak vidíme, že zatímco ke konstrukci nelineárního odhadu jsme potřebovali znát dokonale simultánní hustotu p(x, z), v případě lineárního odhadu úplně postačí první a druhé momenty. To je výrazná úspora dat při reálném výpočtu. Veta 5.6. Estimátor ^XLMS je nestranný. Důkaz. E ( ^XLMS) = E x + xz-1 z E (Z - z) = = x + xz-1 z (z - z) = = E X Odtud plyne, že E ~X = 0 a můžeme psát var ~X = E ( ~X - E ~X)( ~X - E ~X) T = E ~X ~XT = E (X - ^X)(X - ^X) T = = E ((X - x) - xz-1 z (Z - z))((X - x) - xz-1 z (Z - z)) T = = x - xz-1 z zx - xz-1 z zx + xz-1 z z-1 z zx = = x - xz-1 z zx Celkem tedy var ~X = x - xz-1 z zx (5.11) Porovnáním vzorců (5.7),(5.11) se vzorci (5.4),(5.5) dospějeme k důležitému závěru: Je-li vektor (XT , ZT ) T rozložen normálně, pak estimátory ^XMS a ^XLMS splývají. Priklad 5.7. Necht' náhodná veličina X má rozložení X Rs(0, 1). Tuto náhodnou veličinu měříme pomocí náhodné veličiny Z. Vztah mezi X a Z je určen rovnicí měření Z = ln 1 X + V kde V Ex(1) je šum měření. Spočítáme optimální estimátory ^XMS a ^XLMS. Nejprve si napišme hustoty pravděpodobnosti p(x), p(v): p(x) = pX(x) = 1 pro 0 x 1 0 jinak p(v) = pV (x) = e-x pro x 0 0 pro x < 0 Hustota pX je apriorní hustota a její střední hodnota 1 2 je apriorním estimátorem náhodné veličiny X. Rozptyl chyby činí 1 12 . 5.1. OPTIM ÁLNÍ ESTIM ÁTOR 37 Ze znalosti výstupní rovnice a rozložení šumu můžeme ihned určit podmíněnou hustotu p(z|x): p(z|x) = pZ|X(z, x) = pV (z - ln 1 x ) = e-(z-ln 1 x ) = 1 x e-z pro x e-z 0 pro x < e-z sdružená hustota p(x, z) je dána vzorcem p(x, z) = p(z|x)p(x), neboli p(x, z) = 1 x e-z pro x e-z, 0 x 1 0 jinak Marginalizací dostaneme hustotu p(z) = pZ(z): pZ(z) = p(x, z) dx = = 1 e-z 1 x e-z dx = = e-z [ln x]1 e-z = = e-z (0 - ln e-z ) = = ze-z pro z 0 Nyní už snadno spočteme aposteriorní hustotu p(x|z) jako pX|Z(x, z) = p(x, z) p(z) = = 1 x e-z ze-z = = 1 xz pro x e-z, 0 x 1 0 jinak Aposteriorní estimátor 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 e-z x 1 xz dx = = 1 z (1 - e-z ) = = 1 - e-z z 38 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU Spočteme rozptyl aposteriorní chyby. Podle (5.3) máme var ~X = E (var (X|Z)) = = (x - E (X|z))2 p(x|z) dx dz = = x2 p(x|z) dx - E (X|z)2 dz = = 0 1 e-z x2 1 xz dx - (1 - e-z)2 z2 dz = = 0 1 z 1 e-z x dx - (1 - e-z)2 z2 dz = = 0 1 2z (1 - e-2z ) - (1 - e-z)2 z2 dz = = 0 1 - e-2z 2z - 1 - 2e-z + e-2z z2 dz = = 0 (z - 2) + 4e-z - (z - 2)e-2z 2z2 dz = Nyní se budeme věnovat lineárnímu estimátoru ^XLMS, který je určen vzorcem (5.7). Podle předchozího x = 0.5. Zbývá spočítat z, z, xz: z = zpZ(z) dz = = 0 z2 e-z dz = = 2 E XZ = xzp(x, z) dx dz = = 0 1 e-z xzze-z dx dz = = 1 2 0 z2 e-z (1 - e-2z ) dz = = 3 4 = 0.75 xz = E XZ - xz = = 3 4 - 2 2 = = - 1 4 = -0.25 5.2. ODHAD STAVU DYNAMICKÉHO SYSTÉMU 39 E Z2 = z2 pZ(z) dz = = 0 1 e-z z2 ze-z dx dz = = 0 z2 ze-z (1 - e-z dz = = 6 z = E Z2 - 2 z = = 6 - 22 = = 2 Podle (5.7) se lineární estimátor vypočte jako ^XLMS = x + xz-1 z (z - z) = = 1 2 - 1 4 1 2 (z - 2) = = 3 4 - 1 8 z = 0.75 - 0.125z 5.2 Odhad stavu dynamického systému Zkoumejme nyní diskrétní dynamický systém (3.4) z hlediska optimálního odhadu stavů. Dynamický systém necht' je definován rovnicemi xt+1 = f(xt, ut, vt) (5.12a) yt = g(xt, ut, wt) (5.12b) a počátečním rozložením p(x0) (5.12c) Budeme postupně konstruovat optimální odhady stavů x0, x1, . . . , xt rekurentním způsobem. Pro odhad stavů máme k dispozici pouze časovou řadu vstupů a časovou řadu výstupů. Rekurentní způsob výpočtu je optimální pro potřeby praxe. V realitě je předmětem zájmu aktuální stav xt, který se jednoduše vypočte na základě minulého odhadu xt-1 a nových dat ut, yt. Protože rozložení náhodného šumu vt pokládáme za známé, je rovnicí (5.12a) jednoznačně určena podmíněná hustota pravděpodobnosti p(xt+1|xt, ut). Analogicky je rovnicí (5.12b) jed- noznačně určena hustota p(yt|xt, ut). Označme symbolem Dt data, vzniklá pozorováním systému na horizontu 0, . . . , t Dt = {u0, y0, . . . , ut, yt} (5.13) Protože jsou obě rovnice (5.12a) a (5.12b) nezávislé na datech Dt-1, můžeme podle (B.11) psát p(xt+1|xt, ut, Dt-1) = p(xt+1|xt, ut) (5.14) p(yt|xt, ut, Dt-1) = p(yt|xt, ut) (5.15) 40 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU Vztah (5.15) říká, že stav systému xt obsahuje veškerou informaci obsaženou v datech Dt-1, která je potřebná k predikci výstupu yt. Tuto vlastnost můžeme pokládat za ekvivalentní definici stavu. Z rovnice (5.12a) dále vyplývá, že řízení ut neovlivňuje přímo stav xt, což můžeme zapsat jako p(xt|Dt-1, ut) = p(xt|Dt-1) (5.16) Tato vlastnost se někdy označuje jako přirozené podmínky řízení 1. Předpokládejme, že se v současnosti nacházíme v časovém okamžiku t N a máme k dispozici data Dt-1. Chtěli bychom nyní na základě změřených dat Dt-1 co nejlépe odhad- nout neznámou hodnotu stavu xt. Tuto hodnotu nejlépe aproximuje podmíněná hustota pravděpodobnosti p(xt|Dt-1). Nazveme ji apriorní hustota. Předpokládejme dále, že apri- orní hustota je známá funkce. Po doplnění datového souboru aktuálními hodnotami ut, yt máme k dispozici zaktualizovaná data Dt, na jejichž základě máme možnost zaktualizovat také náš odhad stavu xt. Hledáme tedy novou podmíněnou hustotu p(xt|Dt)--tzv. aposteri- orní hustotu. Veta 5.8 (aposteriorní hustota). Aposteriorní hustota p(xt|Dt) je dána vzorcem p(xt|Dt) = p(yt|xt, ut) p(yt|ut, Dt-1) p(xt|Dt-1) (5.17) Důkaz. Použitím Bayesova vzorce (B.10) a vztahů (5.14),(5.16) dostáváme p(xt|Dt) = p(xt|yt, ut, Dt-1) = = p(yt|xt, ut, Dt-1)p(xt|ut, Dt-1 p(yt|ut, Dt-1) = = p(yt|xt, ut) p(yt|ut, Dt-1) p(xt|Dt-1), což je vzorec (5.17) Jmenovatel p(yt|ut, Dt-1) je konstanta, protože yt, ut, Dt-1 jsou známé hodnoty. Toto číslo funguje jako normalizační konstanta. Přechod od apriorní hustoty k aposteriorní hustotě pomocí (5.17) se nazývá datový krok. Je dobré si všimnout, že jsme při něm použili pouze znalost výstupní rovnice (5.12b). Aby bylo možné počítat rekurentně aposteriorní hustoty, je potřeba odvodit přechod od aposteriorní hustoty p(xt|Dt) k apriorní hustotě p(xt+1|Dt). Tento přechod se někdy nazývá časový nebo též predikční krok odhadu stavu. Veta 5.9. Predikční krok odhadu stavu je dán vztahem p(xt+1|Dt) = p(xt+1|xt, ut)p(xt|Dt) dxt (5.18) Důkaz. Podle řetězového pravidla (B.9) platí p(xt+1, xt|Dt) = p(xt+1|xt, Dt)p(xt|Dt) 1 Poznamenejme ovšem, že zejména u ekonomických systémů nebývá často tato podmínka splněna - viz příklad 4.2.1. 5.3. KALMANŮV FILTR 41 Vztah (5.14) dává p(xt+1|xt, Dt) = p(xt+1|xt, ut, yt) stav xt+1 je však nezávislý na výstupu yt, takže můžeme psát p(xt+1|xt, Dt) = p(xt+1|xt, ut) marginalizací podle (B.1) pak dostaneme p(xt+1|Dt) = p(xt+1|xt, ut)p(xt|Dt) dxt což jsme měli dokázat Posledním problémem je zahájení rekurze. Je nutné znát apriorní hustotu p(x0|D-1) = p(x0) (5.19) Apriorní hustota p(x0) je bud'to předem známa, nebo se určuje zkusmo a podle zkušenosti, což zanáší do výpočtu subjektivní prvek. Datový a časový krok algoritmu lze spojit do jediného kroku p(xt+1|Dt) = p(xt+1|xt, ut, yt)p(yt|xt, ut) p(yt|ut, Dt-1) p(xt|Dt-1) dxt = = p(xt+1, yt|ut, Dt-1) p(yt|ut, Dt-1) 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 bud' nepřesné nebo podstatně složitější. Veta 5.10 (Kalmanův filtr). Bud' dán lineární diskrétní stochastický systém, popsaný rovnicemi xt+1 = 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) x0 N(0, 0) (ii) vt N(0, v) pro všechna t (iii) wt N(0, w) pro všechna t (iv) EvtwT t = 0 (v) Ex0vT t = 0 kde vektor 0 a matice 0, v, w jsou známé. Dále necht' jsou známa data Dt. Označme xt|k = xt|Dk. Pak apriorní estimátor xt|t-1 a aposteriorní estimátor xt|t jsou normální pro všechna t xt|t-1 N(t|t-1, t|t-1) (5.21) xt|t N(t|t, t|t) (5.22) kde střední hodnoty t|t-1, t|t a varianční matice t|t-1, t|t se vypočtou podle vzorců t|t = t|t-1 + Kt(yt - Ct|t-1 - Dut) (5.23) t|t = t|t-1 - KtCt|t-1 (5.24) t+1|t = At|t + But (5.25) t+1|t = At|tAT + v (5.26) Kt = t|t-1CT (Ct|t-1CT + w)-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 x1|0, x1|1, x2|1, x2|2, . . . , xt|t-1, xt|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-1 k xt|t. Naproti tomu rovnice (5.25),(5.26) tvoří predikční krok, tedy přechod od xt|t k xt+1|t. Dosadíme-li (5.27) do rovnice (5.24)(viz (5.29)), dostaneme diferenční rovnici pro t|t, která je maticovou diskrétní variantou Riccatiho diferenciální rovnice tvaru ˙x = a0 + a1x + a2x2. Důkaz. Předpokládejme, že platí (5.21). Označme yt|k = yt|Dk. Pak yt|t-1 = yt|Dt-1 = = (Cxt + Dut + wt)|Dt-1 = = Cxt|t-1 + wt|Dt = = Cxt|t-1 + wt Protože yt|t-1 je lineární kombinací náhodných vektorů xt|t-1, wt, které jsou podle předpokladu normální, je yt|t-1 podle (B.6) normální se střední hodnotou Ct|t-1 + Dut a rozptylem var (yt|t-1) = var (Cxt|t-1 + Dut + wt) = = Cvar (xt|t-1)CT + var wt = = Ct|t-1CT + w 5.3. KALMANŮV FILTR 43 Spočítejme kovarianční matici cov(xt|t-1, yt|t-1) = cov(xt|t-1, Cxt|t-1 + Dut + wt) = = var (xt|t-1)CT + cov(xt|t-1, wt) = = t|t-1CT Spojený náhodný vektor (xt|t-1, yt|t-1) je normálně rozdělen xt|t-1 yt|t-1 N t|t-1 Ct|t-1 + Dut , t|t-1 t|t-1CT Ct|t-1 Ct|t-1CT + w Jsou tedy splněny předpoklady věty 5.4 a pro podmíněnou náhodnou veličinu xt|t-1|yt = xt|(yt, Dt-1) = xt|(yt, ut, Dt-1) = xt|Dt = xt|t platí, že je normální xt|t N(t|t, t|t) s parametry t|t = t|t-1 + t|t-1CT (Ct|t-1CT + w)-1 (yt - Ct|t-1 - Dut) (5.28) t|t = t|t-1 - t|t-1CT (Ct|t-1CT + w)-1 Ct|t-1 (5.29) Položíme-li Kt = t|t-1CT (Ct|t-1CT + w)-1 pak můžeme psát t|t = t|t-1 + Kt(yt - Ct|t-1 - Dut) t|t = t|t-1 - KtCt|t-1 což jsou rovnice (5.23) a (5.24), které tvoří datový krok algoritmu. Nyní odvod'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+1|t je lineární kombinací xt|t a vt, které jsou normální. xt+1|t je tedy podle (B.6) také normální s parametry t+1|t = At|t + But t+1|t = At|tAT + v což jsou rovnice (5.25),(5.26). Tím je věta dokázána 44 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU Priklad 5.11. Lod' se pohybuje po rovníku východním směrem rychlostí 10 námořních mil za hodinu. Okamžitou rychost lodi ovšem ovlivňují náhodné poryvy větru a nárazy vln. Navigátor lodi má za úkol každou hodinu odhadnout zeměpisnou délku l v minutách a rychlost lodi s = dl/dt v mph. V čase t = 0 navigátor odhadnul polohu l0 = 0 a rychlost s0 = 10. Každou hodinu pak změřil sextantem zeměpisnou šířku a údaje zaznamenal do tabulky: Čas Poloha 1 9.0" 2 19.5" 3 29.0" Označíme-li lt, st polohu a rychlost lodi v čase t, pak úloha navigátora je úlohou optimálního odhadu veličin lt, st. Počáteční odhady navigátora můžeme matematicky modelovat jako nezávislé náhodné veličiny s normálním rozložením. Rozptyly navigátorových odhadů jsou dlouhodobě sledovány a jejich hodnoty jsou var l0 = 2, var s0 = 3. První krok při řešení problému je namodelovat dynamiku systému. Během hodiny t se lod' pohybuje rychlostí st mph, takže její poloha se změní na2 lt+1 = lt + st (5.30) Veličina st+1 by měla být rovna st, protože lodivod udržuje konstantní rychlost 10 mph. Protože ale je rychlost náhodně ovlivněna účinkem větru a vln, přidáme do rovnice bílý gaussovský šum e, jehož rozptyl byl dlouhodobým měřením odhadnut na číslo 1. st+1 = st + et (5.31) Definujeme-li stavový vektor xt jako xt = lt st , můžeme (5.30),(5.31) napsat jako jedinou stavovou rovnici systému xt+1 = 1 1 0 1 xt + vt, vt N( 0 0 , 0 0 0 1 ) (5.32) Šum procesu vt je singulární gaussovský šum. Počáteční podmínka je podle předchozího x0 N(0, 0) = N( 0 10 , 2 0 0 3 ) Tímto je kompletně popsána dynamika systému. Rovnici měření zapíšme jako yt = 1 0 xt + wt, wt N(0, 2) (5.33) kde w = 2 je rozptyl měření sextantu, udávaný výrobcem. Z rovnic (5.32),(5.33) je zřejmé, že se jedná o lineární gaussovský systém bez deterministického vstupu ut. Ted' je vše připraveno k 2 Jedna ujetá námořní míle odpovídá změně zeměpisné šířky o jednu úhlovou minutu. 5.4. SMOOTHING 45 optimálnímu odhadu stavů Kalmanovým filtrem. Použijeme obecné rekurentní vztahy (5.23)­ (5.27) s počáteční podmínkou 0|-1 = 0 10 , 0|-1 = 2 0 0 3 (5.34) Výsledky uvedeme přehledně v tabulce t yt t|t-1 t|t-1 t|t t|t 0 -- 0.000 10.000 2.000 0.000 0.000 3.000 0.000 10.000 2.000 0.000 0.000 3.000 1 9.0 10.000 10.000 5.000 3.000 3.000 4.000 9.286 9.571 1.429 0.857 0.857 2.714 2 19.5 18.857 9.571 5.857 3.571 3.571 3.714 19.336 9.864 1.491 0.909 0.909 2.091 3 29.0 29.200 9.864 5.400 3.000 3.000 3.091 29.054 9.783 1.460 0.811 0.811 1.875 Všimněme si, že variance v modelovém kroku vzroste v důsledku nejistoty stavové rovnice. Naproti tomu v datovém kroku variance klesá v důsledku zapracování dodatečné informace do naší předpovědi. 5.4 Smoothing Smoothing je zkráceně řečeno optimální odhad stavu xt dynamického systému (5.12), založený na datech DN , kde N > t. Výsledkem je tedy podmíněný náhodný vektor xt|N . V porovnání s apriorním odhadem xt|t-1 a aposteriorním odhadem xt|t máme pro odhad stavu k dispozici více informace a výsledný estimátor 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. Smooth- ing se tedy nejspíše uplatní tam, kde nás kromě aktuálního stavu xt zajímají také minulé stavy x0, . . . , xt-1. Jak si ukážeme, lze odhady x0|t, . . . , xt-1|t získat úpravou již existujících aposteriorních odhadů x0|0, . . . , xt-1|t-1. Algoritmus přitom postupuje zpětně v čase. Veta 5.12. Necht' t < N. Pro podmíněnou hustotu p(xt|DN ) platí: p(xt|DN ) = p(xt|Dt) p(xt+1|xt, ut) p(xt+1|DN ) p(xt+1|Dt) dxt+1 (5.35) Důkaz. p(xt|DN ) = p(xt, xt+1|DN ) dxt+1 (5.36) Počítejme p(xt, xt+1|DN ) = p(xt|xt+1, DN )p(xt+1|DN ) (5.37) Rozberme nyní podrobněji podmíněnou hustotu p(xt|xt+1, DN ). Veškerá informace uložená v datech {yt+1, ut+1, . . . , yN , uN } potřebná k predikci xt je koncentrovaná v hodnotě xt+1, což vyplývá z konstrukce modelových rovnic (5.12). Můžeme tedy prohlásit, že p(xt|xt+1, DN ) = p(xt|xt+1, Dt) (5.38) 46 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU přitom p(xt|xt+1, Dt) = p(xt+1|xt, Dt) p(xt|Dt) p(xt+1|Dt) = = p(xt+1|xt, ut) p(xt|Dt) p(xt+1|Dt) dosazením do (5.37) a (5.36) nakonec dostáváme p(xt|DN ) = p(xt+1|xt, ut) p(xt|Dt) p(xt+1|Dt) p(xt+1|DN ) dxt+1 = = p(xt|Dt) p(xt+1|xt, ut) p(xt+1|DN ) p(xt+1|Dt) dxt+1 což jsme měli dokázat 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+1|N a xt|t. Tyto dva estimátory jsou také postačující pro výpočet estimátoru 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átoru. 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(t|N , t|N ) vypočte podle vzorců t|N = t|t + Ft(t+1|N - t+1|t) (5.39) t|N = t|t - Ft(t+1|t - t+1|N )Ft T (5.40) Ft = t|tAT -1 t+1|t (5.41) 5.5 Rozšířený Kalmanův filtr Myšlenka rozšířeného Kalmanova filtru (Extended Kalman 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. Bud' dán lineární diskrétní stochastický systém, popsaný rovnicemi xt+1 = Atxt + bt + vt (5.42) yt = Ctxt + dt + wt (5.43) splňující podmínky (i)­(v) z věty 5.10. Pak pro apriorní estimátor xt|t-1 N(t|t-1, t|t-1) a aposteriorní estimátor xt|t N(t|t, t|t) platí t|t = t|t-1 + Kt(yt - Ctt|t-1 - dt) (5.44) t|t = t|t-1 - KtCtt|t-1 (5.45) t+1|t = Att|t + bt (5.46) t+1|t = Att|tAt T + v (5.47) Kt = t|t-1Ct T (Ctt|t-1Ct T + w)-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. Vezměme obecný nelineární dynamický systém ve stavovém tvaru xt+1 = 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) + f x ()(x - ) + O(x2 ) (5.49) g(x, u) = g(, u) + g x ()(x - ) + O(x2 ) (5.50) Aproximace bude tím lepší, čím menší bude rozdíl x - . Zanedbáme-li členy O(x2), můžeme původní systém přibližně nahradit lineárním systémem xt+1 = f(, ut) + At(xt - ) + vt At = f x () (5.51) yt = g(, ut) + Ct(xt - ) + wt Ct = g x () (5.52) Předpokládejme, že máme k dispozici estimátor xt|t-1 prostřednictvím jeho momentů t|t-1, t|t-1. Linearizujme výstupní rovnici v okolí bodu t|t-1: yt = Ctxt + (g(t|t-1, ut) - Ctt|t-1) + wt kde Ct = g x (t|t-1) (5.53) Pak podle lemmatu 5.14 je datový krok dán rovnicí t|t = t|t-1 + Kt(yt - g(t|t-1, ut)) (5.54) t|t = t|t-1 - KtCtt|t-1 (5.55) kde jsme použili Kalmanovo zesílení Kt = t|t-1Ct T (Ctt|t-1Ct T + w)-1 Nyní máme k dispozici aposteriorní odhad xt|t. Linearizujeme-li stavovou rovnici v okolí bodu t|t, obdržíme lineární stavovou rovnici xt+1 = Atxt + (f(t|t, ut) - Att|t) + vt kde At = f x (t|t) (5.56) Podle lemmatu 5.14 napíšeme predikční krok t+1|t = f(t|t, ut) (5.57) t+1|t = Att|tAt T + v (5.58) Někdy se odhady vylepšují ještě smoothingem podle věty 5.13, takže algoritmus rozšířeného Kalmanova filtru se pak rozpadne na dopředný běh (filtr) a zpětný běh (smoother). Dosažené výsledky přehledně zformulujeme v definici 48 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU Definice 5.15. Pro nelineární dynamický systém xt+1 = f(xt, ut) + vt yt = g(xt, ut) + wt s počáteční podmínkou x0 N(0, 0) definujeme Rozšířený Kalmanův filtr jako algoritmus výpočtu estimátorů xt|t-1 N(t|t-1, t|t-1) xt|t N(t|t, t|t) xt|N N(t|N , t|N ) podle vzorců t|t = t|t-1 + Kt(yt - g(t|t-1, ut)) (5.59) t|t = t|t-1 - KtCtt|t-1 (5.60) t+1|t = f(t|t, ut) (5.61) t+1|t = Att|tAt T + v (5.62) Kt = t|t-1Ct T (Ctt|t-1Ct T + w)-1 (5.63) t|N = t|t + Ft(t+1|N - t+1|t) (5.64) t|N = t|t - Ft(t+1|t - t+1|N )Ft T (5.65) Ft = t|tAT -1 t+1|t (5.66) kde matice At, Ct jsou Jakobiány At = f x (t|t) Ct = g x (t|t-1) (5.67) 5.6 Odhad nelineárních systémů metodou Monte Carlo Rozšířený Kalmanův filtr je průmyslovým standartem pro nelineární rekurzivní estimaci. V roce 1992 se v článku [?] objevil zcela nový algoritmus, založený na aplikaci známé metody Monte Carlo při bayesovské estimaci. Metoda Monte Carlo je založena na tom, že informaci o rozložení náhodné veličiny nese náhodný výběr z tohoto rozložení. Čím je náhodný výběr delší, tím je informace přesnější. V předchozích metodách informaci o pravděpodobnostním rozložení estimátoru obsahovala bud'to hustota pravděpodobnosti, nebo první dva momenty v gaussovském případě. Pro náhodný vektor x Lnx máme náhodný výběr Sn = (x1, . . . , xn) délky n, obsahující vzorky xi Rnx . Empirická hustota pravděpodobnosti, získaná z náhodného výběru, je aproximací skutečné hustoty pravděpodobnosti náhodného vektoru x. Empirickou hustotu pravděpodobosti můžeme napsat jako pn(x) = 1 n n i=1 (x - xi) (5.68) kde (x) : Rnx R je tzv. Diracova funkce , definovaná limitou (x) = lim h0 h(x), h(x) = 1 h pro 0 < x < h 0 jinak (5.69) 5.6. ODHAD NELINE ÁRNÍCH SYSTÉMŮ METODOU MONTE CARLO 49 V článku [?] je metoda Monte Carlo ještě vylepšena tím, že se jednotlivým vzorkům xi přiřadí váha wi 0 tak, aby součet jednotlivých vah dával jedničku. Empirická hustota je potom pn(x) = n i=1 wi (x - xi) (5.70) Rozložení náhodného vektoru x je tedy popsáno náhodným výběrem Sn a vektorem vah wn = (w1, . . . , wn). Připomeňme obecné rekurzivní vztahy pro nelineární estimaci: p(xt|Dt) = p(yt|xt, ut) p(yt|ut, Dt-1) p(xt|Dt-1) (5.71) p(xt+1|Dt) = p(xt+1|xt, ut)p(xt|Dt) dxt (5.72) p(xt|DN ) = p(xt|Dt) p(xt+1|xt, ut) p(xt+1|DN ) p(xt+1|Dt) dxt+1 (5.73) Pracujeme-li s empirickými hustotami, pak aplikace těchto rekurzivních vztahů vyústí v algo- ritmus 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: Bud' dána apriorní hustota pn(xt|t-1) = n i=1 wi(t|t - 1)(x - xi(t|t - 1)) V datovém kroku hledáme aposteriorní hustotu pn(xt|t). Položíme xi(t|t) = xi(t|t - 1) a váhy wi(t|t) získáme pomocí (5.71). Člen p(yt|ut, Dt-1) funguje ve vzorci (5.71) jako normalizační konstanta. Když tento člen vypustíme, obdržíme váhy wi(t|t) = p(y(t)|xi(t|t - 1), u(t)) wi(t|t - 1) (5.74) z nichž normalizací dostaneme výsledné váhy wi(t|t) = wi(t|t) n j=1 wj(t|t) (5.75) 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+1|t = f(xt|t, u(t)) + v(t) který aplikujeme na vzorky xi(t|t). Dostaneme pak xi(t + 1|t) = 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 + 1|t) = wi(t|t). Takto jsme dostali predikci xt+1|t. Konstrukce smoothingu má stejnou filozofii jako datový krok. Vzorky zůstávají stejné xi(t|N) = xi(t|t) a upravují se jen váhy. Vypuštěním členu p(xt+1|Dt) ze vzorce (5.73) dostaneme zmenšené váhy wi(t|N) podle vzorce wi(t|N) = wi(t|t) n j=1 wj(t + 1|N) p(xj(t + 1|N)|xi(t|N), u(t)) (5.77) 50 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU a váhy wi(t|N) dostaneme normalizací wi(t|N) = wi(t|t) n j=1 wj(t|t) (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 f. Při smoothingu se neobejdeme bez podmíněné hustoty p(xt+1|xt, ut), která je určená stavovou rovnicí. Poslední nevyřešenou otázkou je zahájení algoritmu. Předpokládáme, že rozložení stavu x0 je známé. Pak v souladu s (5.19) položíme x0|-1 = x0. 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) = 1/n. Jiná možnost je zvolit rovnoměrnou sít' bodů xi(0| - 1) ve vhodně zvolené oblasti a váhy volit wi(0| - 1) = p0(xi(0| - 1)) a potom je normalizovat. Dosažené výsledky shrneme v definici. Definice 5.16 (Weighted Bootstrap Algorithm). Bud' dán dynamický systém xt+1 = f(xt, ut, vt) yt = g(xt, ut, wt) s počátečním stavem x0, který má empirickou hustotu pn(x) = n i=1 wi(0| - 1)(x - xi(0| - 1)) Bud' dána data Dt (viz (5.13)). Pak odhady xt|t-1, xt|t, xt|N s empirickými hustotami pn(xt|t-1) = n i=1 wi(t|t - 1)(x - xi(t|t - 1)) pn(xt|t) = n i=1 wi(t|t)(x - xi(t|t)) pn(xt|N ) = n i=1 wi(t|N)(x - xi(t|N)) jsou spočteny algoritmem Weighted Bootstrap, právě když platí xi(t|t) = xi(t|t - 1) (5.79) wi(t|t) = p(y(t)|xi(t|t - 1), u(t)) wi(t|t - 1) (5.80) wi(t|t) = wi(t|t) n j=1 wj(t|t) (5.81) xi(t + 1|t) = f(xi(t|t), u(t)) + vi(t) (5.82) wi(t + 1|t) = wi(t|t) (5.83) xi(t|N) = xi(t|t) (5.84) wi(t|N) = wi(t|t) n j=1 wj(t + 1|N) p(xj(t + 1|N)|xi(t|N), u(t)) (5.85) wi(t|N) = wi(t|t) n j=1 wj(t|t) (5.86) 5.7. SHRNUTÍ 51 5.7 Shrnutí Optimální estimátor 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í estimátor je lineární. Lineární estimátor se dá použít i pro jiná rozdělení, je však obecně méně přesný než estimátor X|Z. Při odhadu stavů dynamického systému (5.12) používáme apriorní odhad xt|t-1, 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-1 a xt|t se jmenuje Kalmanův filtr a zpětný algoritmus pro xt|N se jmenuje Rauch-Tung-Striebel smoother. Pro nelineární systémy s diferencovatelnými funkcemi f, g se může použít rozšířený Kalmanův filtr. Ten v každém kroku lokálně linearizuje funkce f, g 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í estimátor 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í Priklad 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, m0, S0, Y, U), kde uvedené matice mají následující význam: Mp Matice typu nx × N + 1, obsahující ve svých sloupcích predikce 0|-1, . . . , N|N-1. Sp Matice typu n2 x × N + 1, obsahující ve sloupcích varianční matice t|t-1. Mf Matice typu nx×N+1, obsahující ve svých sloupcích filtrované střední hodnoty 0|0, . . . , N|N . Sf Matice typu n2 x × N + 1, obsahující ve sloupcích varianční matice t|t. A Matice A typu nx × nx, definující stavovou rovnici (5.20a). B Matice B typu nx × nu, definující stavovou rovnici (5.20a). C Matice C typu ny × nx, definující výstupní rovnici (5.20b). D Matice D typu ny × nu, definující výstupní rovnici (5.20b). Sv Varianční matice stavového šumu v typu nx × nx, symetrická a pozitivně semidefinitní. Sw Varianční matice výstupního šumu w typu ny ×ny, symetrická a pozitivně semidefinitní. m0 Střední hodnota počátečního stavu x0 typu nx × 1. 52 KAPITOLA 5. OPTIM ÁLNÍ ODHAD STAVŮ SYSTÉMU S0 Variance počátečního stavu x0 typu nx × nx. Y Datová matice výstupů typu ny × N + 1, obsahující ve sloupcích t výstupy yt. U Datová matice vstupů typu nu × 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 : xt+1 = f(xt, ut, vt, ) = f(xt, ut, vt) (6.1a) yt = g(xt, ut, wt, ) = g(xt, ut, wt) (6.1b) Parametr je veličina, na které jsou závislé vlastnosti systému (6.1). Množina se pak nazývá parametrický prostor. Priklad 6.1. Mějme dynamický systém, o němž víme, že je lineární: xt+1 = Axt + But + vt yt = Cxt + Dut + wt, ale neznáme matice A, B, C, D. Tento systém je závislý na parametru = A B C D a parametrický prostor je potom = Rnx×nx × Rnx×nu × Rny×nx × Rny×nu s dimenzí p = nxnx + nxnu + nynx + nynu. Většinou se stává, že u dynamického systému neznáme konkrétně funkce f, g, ale známe strukturu systému f, g. Struktura systému spolu s vektorem jednoznačně určuje vlastnosti systému. Chybějí-li tedy informace o hodnotě , 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 ^ odhad parametru a Dn = {u0, y0, . . . , un, yn} historická data, pak identifikace systému (6.1) je úloha nalézt při známé struktuře f, g odhad ^ = F(Dn) parametru tak, aby se co nejméně lišil od skutečného parametru . 6.2 Klasický přístup Při klasickém přístupu se neznámý parametr uvažuje jako reálný vektor dimenze p. Para- metrický prostor je potom roven Rp. 53 54 KAPITOLA 6. IDENTIFIKACE DYNAMICK ÝCH SYSTÉMŮ Při identifikaci je cílem dosáhnout maximální shody mezi odhadem ^ a skutečnou hodno- tou parametru . Protože však rozdíl - ^ není přímo měřitelný, používají se různá kriteria, která jsou vodítkem při určení ^. 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 ^yt = ^yt(), pak shodu modelu se skutečností měříme podle kriteria nejmenších čtverců jako J() = n t=0 |yt - ^yt()|2 (6.2) Namísto neměřitelného rozdílu - ^ tak pracujeme s měřitelným rozdílem yt - ^yt. Optimální odhad ^ bude takový vektor, pro který bude kriterium J minimální, tj. J(^) = min J() (6.3) Funkce J je zobrazení J : Rp R. Tuto funkci minimalizujeme metodami matematického programování. Priklad 6.2 (Lineární systém). Pro jednoduchost předpokládejme nejprve lineární systém s pozorovatelnými stavy, neboli systém xt+1 = Axt + But + vt (6.4a) yt = xt (6.4b) Máme identifikovat matice A, B na základě historie systému, zachycené v datech Dn = {u0, y0, . . . , un, yn} = {u0, x0, . . . , un, xn} Kriteriální funkce J bude nyní vážit rozdíl mezi stavem xt+1 a jeho -předpovědí Axt + But: J() = J(A, B) = n-1 t=0 |xt+1 - Axt - But|2 Označíme-li et = xt+1 - Axt - But a matice et, xt+1, A, B rozepíšeme po řádcích et = e1,t ... enx,t xt+1 = x1,t+1 ... xnx,t+1 A = a1 ... anx B = b1 ... bnx pak platí J(A, B) = n-1 t=0 |et|2 = n-1 t=0 nx j=1 e2 jt = nx j=1 n-1 t=0 (xj,t+1 - ajxt - bjut)2 = nx j=1 Jj(aj, bj) 6.2. KLASICK Ý P ŘÍSTUP 55 kde jsme položili Jj(aj, bj) = (xj,t+1 - ajxt - bjut)2 Úloha minimalizovat J(A, B) se nám nyní rozpadla na nx úloh minimalizovat Jj(aj, bj) pro j {1, . . . , nx}. Tyto dílčí úlohy jsou úlohami lineární regrese (viz B.6) ve tvaru xj,t+1 = ajxt + bjut (6.5) při označení zj = xj,1 ... xj,n j = aj bj T U = x0 T u0 T ... ... xn-1 T un-1 T můžeme tyto úlohy zapsat maticově zj = Uj (6.6) Podle (B.18) je optimálním odhadem vektoru j vektor ^j, daný vztahem ^j = (UT U)-1 UT zj (6.7) Odhad ^ = [ ^A ^B] složíme s odhadů ^j: ^A ^B = ^T 1 ... ^T nx = ^1 ^nx T = (UT U)-1UT z1 (UT U)-1UT znx T = = (UT U)-1 UT z1 znx T = z1 znx T U(UT U)-1 = = x1,1 xnx,1 ... ... x1,n xnx,n T U(UT U)-1 = XU(UT U)-1 Celkem tedy máme ^A ^B = XU(UT U)-1 , (6.8) kde jsme položili X = x1 T ... xn T U = x0 T u0 T ... ... xn-1 T un-1 T Priklad 6.3 (Nelineární systém). Pro jednoduchost předpokládejme opět systém s po- zorovatelnými stavy, s konstantními parametry a s aditivním šumem xt+1 = f(xt, ut, ) + vt (6.9a) yt = xt (6.9b) 56 KAPITOLA 6. IDENTIFIKACE DYNAMICK ÝCH SYSTÉMŮ Na základě známých dat Dn máme odhadnout parametr tak, aby součet čtverců J() = n-1 t=0 |xt+1 - f(xt, ut, )|2 (6.10) byl minimální na Rp. Jedná se o hledání volného extrému funkce p proměnných. Je-li funkce f třídy C1, pak z analýzy víme, že je-li bod ^ 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: J i = -2 n-1 t=0 nx j=1 fj i (xt, ut, )(xj,t+1 - fj(xt, ut, )) (6.11) Označíme-li At() = fj i (xt, ut, ) bt() = xj,t+1 - fj(xt, ut, ) F() = n-1 t=0 At()bt() pak můžeme podmínku J = 0 zapsat jako soustavu nelineárních rovnic řádu p F() = 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 fsolve, 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ámé jsou 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() funkce kvadratická v , 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 je náhodný vektor, který se vyvíjí v čase. Jinými slovy posloupnost {t}n-1 t=0 je náhodný proces. Takový proces je popsaný například sdruženou hustotou p(0, . . . , n-1). Tento náhodný proces je před identifikací již známý (je známa hustota pravděpodobnosti p). Podoba náhodného procesu {t}n-1 t=0 patří mezi apriorní informace, které ovlivní konečný odhad. Při konstrukci modelu vývoje parametru máme možnost ovlivnit závislost mezi parametry v různých obdobích i závislost mezi jednotlivými složkami parametru . Můžeme například předepsat, že parametry t a t+1 budou silně korelovány, zatímco složky 1,t, . . . , p,t 6.3. BAYESŮV P ŘÍSTUP 57 budou nezávislé a podobně. U Bayesovy metody se předpokládá, že náhodný proces {t}n-1 t=0 je zadán jako dynamický systém (ve stavovém tvaru) t+1 = h(t, t) (6.13) s počáteční podmínkou 0 = , kde Lp je náhodný vektor se známým rozložením a t 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 t+1 = t + t (6.14) Model náhodné procházky vykazuje některé charakteristické vlastnosti: (i) E (t) = E (0) pro t {0, . . . , n - 1}. (ii) var (t) = var (0) + t-1 j=0 var (j) pro t {0, . . . , n - 1}. (iii) cov(t, t+s) = var (t) pro t {0, . . . , n - 1}, s {0, . . . , n - t - 1}. Předpokládejme, že je znám model vývoje parametrů (6.13) a je dán parametrický dy- namický systém (6.1). Máme tedy systém ve tvaru xt+1 = f(xt, ut, vt, t) (6.15a) t+1 = h(t, t) (6.15b) yt = g(xt, ut, wt, t) (6.15c) Definujme náhodný vektory zt, et předpisem zt = xt t et = vt t (6.16) Definujme dále funkci F : Rnx+p × Rnu × Rnv Rnx+p následujícím způsobem: F(zt, ut, et) = f(xt, ut, vt, t) h(t, t) (6.17) a nakonec funkci G : Rnx+p × Rnu × Rnw Rny předpisem G(zt, ut, wt) = g(xt, ut, wt, t) (6.18) Potom dynamický systém (6.15) napíšeme pomocí (6.16), (6.17), (6.18) jako zt+1 = F(zt, ut, et) (6.19a) yt = G(zt, ut, wt) (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ů t 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Ů Priklad 6.4. Vezměme nelineární parametrický dynamický systém s pozorovatelnými stavy a se známým modelem vývoje parametrů xt+1 = f(xt, ut, vt, t) t+1 = h(t, t) yt = xt Pro tento systém necht' známe data Dn = {u0, y0, . . . , un, yn}. Odhadněme parametry t Bayesovou metodou. Zaved'me výstupní proměnnou t = xt+1 a stav zt = t. Dále definujme funkci g jako g(zt, ut, vt, t) = f(xt, ut, vt, zt) Funkce g je jednoznačně určena funkcí f a daty Dn. Nyní má systém podobu zt+1 = h(zt, t) t = g(zt, ut, vt, t) Jedná se tedy o neautonomní systém s nepozorovatelnými stavy zt s šumem procesu t a výstupním šumem vt, který byl původně šumem procesu. Pro odhad stavu zt (parametru t) máme k dispozici vstupy u0, . . . , un a výstupy 0, . . . , n-1. Časový horizont je tedy zkrácený o jednu periodu. Priklad 6.5. Bud' 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, 2 v), w(t) N(0, 2 w) s časově proměnnými parametry (t) = a(t) b(t) c(t) d(t) Předpokládáme model vývoje parametrů ve tvaru náhodné procházky t+1 = (t) + (t) kde (t) N(0, 2 ) je bílý šum. Známe počáteční podmínky x0 N(0, 2 0) 0 N(0, 2 0 ) Sestrojme algoritmus pro optimální odhad parametru (t) s použitím Bayesovy metody a Rozšířeného Kalmanova filtru. Definujeme si stavový vektor z(t) = x(t) a(t) b(t) c(t) d(t) 6.3. BAYESŮV P ŘÍSTUP 59 a vektor šumu e(t) = v(t) 1(t) 2(t) 3(t) 4(t) Dále definujeme funkci F : R5 × R × R5 vzorcem F(z(t), u(t), e(t)) = z1(t)z2(t) + z3(t)u(t) + e1(t) z2(t) + e2(t) z3(t) + e3(t) z4(t) + e4(t) z5(t) + e5(t) a funkci G : R5 × R × R vzorcem G(z(t), u(t), w(t)) = z1(t)z4(t) + z5(t)u(t) + w(t) Dynamický systém z(t + 1) = F(z(t), u(t), e(t)) 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 z1(t + 1) = z1(t)z2(t) + z3(t)u(t) + e1(t) z2(t + 1) = z2(t) + e2(t) z3(t + 1) = z3(t) + e3(t) z4(t + 1) = z4(t) + e4(t) z5(t + 1) = z5(t) + e5(t) y(t) = z1(t)z4(t) + z5(t)u(t) + w(t), který zřejmě není lineární. Pro optimální odhad stavu z(t) se tedy hodí Rozšířený Kalmanův Filtr. Matice A(t), C(t) jsou: A(t) = z2(t) z1(t) u(t) 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 C(t) = z4(t) 0 0 z1(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.1 Blokové matice Necht matice A je zadána blokově A = A11 A12 A21 A22 (A.1) Pak definujeme Schurův komplement A22 jako D22 = A22 - A21A-1 11 A12 (A.2) a Schurův komplement A11 jako D11 = A11 - A12A-1 22 A21 (A.3) Inverze k A pak může být psána ekvivalentně jedním ze tří vzorců A-1 = A-1 11 + A-1 11 A12D-1 22 A21A-1 11 -A-1 11 A12D-1 22 -D-1 22 A21A-1 11 D-1 22 (A.4) A-1 = D-1 11 -D-1 11 A12A-1 22 -A-1 22 A21D-1 11 A-1 22 + A-1 22 A21D-1 11 A12A-1 22 (A.5) A-1 = D-1 11 -D-1 11 A12A-1 22 -A-1 22 A21D-1 11 D-1 22 (A.6) Porovnáním těchto tří tvarů se dá odvodit lemma o inverzní matici (A-1 11 - A12A22A21)-1 = A11 - A11A12(A21A11A12 + A-1 22 )-1 A21A11 (A.7) Dále platí |A| = |A11| |D22| (A.8) |A| = |A22| |D11| (A.9) 63 64 DODATEK A. P ŘEHLED MATICOVÉ ALGEBRY A.2 Maticový kalkulus Bud' f : Rn R funkce, x0 Rn je reálný vektor. Definujeme gradient funkce f v bodě x0 jako vektor f x (x0) = ( f x1 (x0), . . . , f xn (x0)) T Bud' F : Rm×n R funkce, X0 Rm×n je reálná matice. Gradient funkce F v bodě X0 je matice F X (X0) = F Xij (X0) i,j Některé užitečné gradienty: x (yT x) = y (A.10) x (xT y) = y (A.11) x (yT Ax) = Ay (A.12) x (xT Ay) = Ay (A.13) x (xT Ax) = (A + AT )x (A.14) A tr (A) = I (A.15) A tr (BAD) = BT DT (A.16) A tr (ABAT ) = 2AB (A.17) A |BAD| = |BAD| A-T (A.18) kde jsme položili A-T = (A-1) T . Bud' g : Rn Rm funkce, x0 Rn je reálný vektor. Jakobián funkce g v bodě x0 je matice g x (x0) = gi xj (x0) i,j Některé užitečné Jakobiány: x (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 (A.20) yT x = tr (yxT ) (A.21) Dodatek B Pravděpodobnost B.1 Základní pojmy Označme (, A, P) jako pravděpodobnostní prostor, kde značí základní prostor, A jevové pole na a P je pravděpodobnostní míra (provděpodobnost) na A. Náhodná veličina (náhodný vektor) je zobrazení X : Rn. Distribuční funkce náhodného vektoru X Rn je zobrazení F : Rn R definované vztahem F(x) = P({ : X() x}) Hustota pravděpodobnosti náhodného vektoru X Rn je zobrazení F : Rn R a je definována jako p(x) = nF x1 xn (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 , Y T ) T 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.1) p(y) = 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) = Rn x dF(x) Druhý moment náhodného vektoru x je matice m2(x) = Rn xxT dF(x) 65 66 DODATEK B. PRAVDĚPODOBNOST Rozptyl náhodného vektoru x je definován jako var (x) = E (x - E (x))(x - E (x))T Platí, že var (x) = m2(x) - [E (x)]2. Zápisem X Ln rozumíme, že X : Rn je náhodný vektor dimenze n s konečným druhým momentem. Takový vektor má i konečný první moment. {xj} j=- je spočetná posloupnost náhodných vektorů (náhodný proces). B.3 Normální rozdělení Náhodný vektor X Ln má normální rozdělení, právě když hustota pravděpodobnosti je dána vzorcem p(x) = 1 (2)n || e 1 2 (x-)T -1(x-) (B.3) kde Rn je reálný vektor a Rn×n je symetrická pozitivně definitní matice. Pak píšeme X N(, ). Lze ukázat, že pro parametry , platí E X = (B.4) var X = (B.5) Normálně rozdělené náhodné vektory X a Y jsou nezávislé právě thedy, jsou-li nekorelo- vané, tj. cov(X, Y ) = 0. Mějme náhodný vektor X Ln, X N(x, x), vektor a Rm a matici B Rm×n. Pak náhodný vektor Y = a + BX má normální rozdělení a platí Y N(a + Bx, BxBT ). Zkráceně X N(x, x) = a + BX N(a + Bx, BxBT ) (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 p(x|y) = p(x, y) p(y) (B.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) = p(x|y)p(y) p(x) = p(x|y)p(y) p(x|y)p(y) dy (B.10) Jsou-li náhodné vektory X a Y nezávislé, pak platí p(x|y) = p(x) (B.11) 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 (B.12) z pravidel pro počítání s podmíněnými hustotami vyplývají ekvivalentní vztahy E (X|Z) = x p(x, z) p(z) dx (B.13) E (X|Z) = xp(x, z) dx p(x, z) dx (B.14) E (X|Z) = xp(x|z)p(x) dx p(x|z)p(x) dx (B.15) B.6 Lineární regrese Necht' y, e L1 jsou náhodné veličiny, x1, . . . , xp jsou reálné proměnné, b1, . . . , bp jsou reálné konstanty (parametry, regresní koeficienty). Lineární regresní model je rovnice Y = b1x1 + . . . + 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) = 2 kde 2 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 = y1 ... yn X = x11 x1p ... ... xn1 xnp Při měření předpokládáme, že šumy ei 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 = (b1, . . . , bp)T tak aby výraz q(b) = n i=1 (yi - p j=1 bjxij)2 = n i=1 (yi - ^yi)2 (B.17) byl minimální. Funkce q(b) se nazývá reziduální součet čtverců. Je to kvadratická forma v proměnných b a její globální minimum se nachází v bodě ^b, který se vypočte jako ^b = (XT X)-1 XT y (B.18) Je vidět, že ^b je lineární funkcí y. Výraz s2 = q(^b) = 1 n - p n i=1 (yi - p j=1 ^bjxij)2 = n i=1 (yi - ^yi)2 (B.19) se nazývá reziduální rozptyl a je nestranným odhadem hodnoty 2.