\section{Základní úloha stochastického řízení} Tento oddíl se zabývá obecnou formulací úlohy stochastického řízení a pojmy s tím spojenými. \subsection{Systém a jeho popis} Ústředním pojmem v teorii řízení je systém. Systém je část světa, kterou chceme poznat či řídit. Ovlivňování systému, ať už za účelem jeho lepšího poznání, či za účelem řízení, provádíme pomocí vstupů (řídících zásahů). Ve většině případů je řešení úlohy stochastického řízení prováděno numericky, je proto účelné pracovat s diskrétním časem. Budeme-li proto uvažovat diskrétní povahu času, stav systému v časovém okamžiku $t$ podél konečného řídícího horizontu délky $N$ popisuje soustava rovnic \begin{equation} \label{sys} x_{t+1}=f_t(x_t,u_t,w_t), \qquad t=0,1,\ldots,N-1. \end{equation} Zde $x_t$ představuje stav systému v čase $t$, $u_t$ řídící zásah v čase $t$ a $w_t$ náhodnou veličinu reprezentující přítomnost šumu. Předpokládáme, že tvar rovnic $f_t$ je nám znám, například z fyzikálního rozboru úlohy, či ze znalosti konstrukce stroje, který popisujeme. Dále předpokládáme, že stav systému můžeme přímo pozorovat. Případ neúplného pozorování bude probrán později. \subsection{Ztrátová funkce a optimální řízení} Cílem řízení je pro systém popsaný soustavou \eqref{sys} navrhnout regulátor (posloupnost řídících zásahů), který bude stav systému udržovat co nejblíže požadované hodnotě. Pro tyto účely máme v úloze řízení k dispozici předepsanou ztrátovou (resp. účelovou) funkci \begin{equation} g(x_{1:N},u_{0:N-1}), \end{equation} která určuje, nakolik jsme vytyčených cílů dosáhli. Označme $U(x_t)$ neprázdnou množinu přípustných řídících zásahů pro systém nacházející se ve stavu $x_t$. Přípustnou řídící strategií $\pi=\mu_{0:N-1}$ budeme rozumět posloupnost zobrazení \begin{equation} \label{con} \mu_t(x_t)=u_t, \qquad t=0,1,\ldots,N-1, \end{equation} kde $\mu_t(x_t)=u_t \in U(x_t)$ je přípustný řídící zásah. Neprázdná množina $\Pi$ pak bude značit množinu všech přípustných řídících strategií. Pro danou řídící strategii označme očekávanou ztrátu jako \begin{equation} \label{los} J_\pi(x_0)=\E_{w_{0:N-1}}\left\{g(x_{1:N},\mu_{0:N-1}(x_{0:N-1}))\right\}. \end{equation} Úlohou řízení je potom najít takovou $\pi^*$, pro kterou platí \begin{equation} J_{\pi^*}(x_0)=\min_{\pi \in \Pi}J_\pi(x_0). \end{equation} Celkově se tedy jedná o optimalizační úlohu nalézt takovou posloupnost funkcí \eqref{con}, která minimalizuje očekávanou ztrátovou \eqref{los} za podmínek \eqref{sys}. \section{Úloha stochastického řízení s aditivní ztrátou} Úlohu stochastického řízení tak, jak byla definována v předchozí části, nelze obecně řešit. Je tedy potřeba úlohu nějak blíže specifikovat. \subsection{Aditivní ztrátová funkce} Jako vhodné se ukazuje omezit se na nějaký speciální tvar ztrátové funkce \eqref{los}. Budeme proto dále uvažovat tzv. aditivní tvar ztrátové funkce, tedy že existují funkce $g_t$ takové, že můžeme psát \begin{equation} \label{adi} g(x_{1:N},u_{0:N-1})=\sum_{t=0}^{N-1}g_t(x_{t+1},u_t). \end{equation} Očekávanou ztrátu \eqref{los} potom můžeme přepsat do tvaru \begin{equation} \label{ex} J_\pi(x_0)=\E_{w_{0:N-1}}\left\{\sum_{t=0}^{N-1}g_t(x_{t+1},\mu_t(x_t))\right\}. \end{equation} \subsection{Dynamické programování} Takto specifikovaná úloha stochastického řízení se dá řešit použitím dynamického programování \cite{bellman1957dynamic}. Dynamické programování je přístup k řešení optimalizačních úloh, na které se můžeme dívat jako na posloupnost rozhodnutí, pro které platí tzv. princip optimality. Ten říká, že optimální posloupnost rozhodnutí má tu vlastnost, že pro libovolný počáteční stav a rozhodnutí musí být všechna následující rozhodnutí optimální vzhledem k výsledkům rozhodnutí prvního. Platnost principu optimality pro očekávanou ztrátu tvaru \eqref{ex} je intuitivně snadno pochopitelná. Pokud by totiž nějaký úsek řídící strategie nebyl optimální, pak očekávanou ztrátu snížíme přechodem ke strategii, ve které onu neoptimální část nahradíme optimálním řešením podproblému na tomtéž úseku. Přesný důkaz platnosti principu optimality pro očekávanou ztrátu tvaru \eqref{ex} lze nalézt například v \cite{bertsekas1995dynamic}. \subsection{Použití dynamického programování při řešení úlohy stochastického řízení s aditivní ztrátou} Při řešení úlohy stochastického řízení s aditivní ztrátou je možné postupovat, jak je u úloh řešených pomocí dynamického programování zvykem. Za tímto účelem označme $J_t(x_t)$ minimální hodnotu střední ztráty od okamžiku $t$ do $N$ v závislosti na $x_t$. Dle \eqref{ex} pro ni můžeme psát \begin{gather} J_N(x_N)=0,\\ J_t(x_t)=\min_{u_t \in U(x_t)}\E_{w_t}\left\{g_k(x_{t+1},u_t)+J_{t+1}(x_{t+1})\right\}, \qquad t=0,\ldots,N-1. \end{gather} Při konstrukci optimální řídící strategie budeme postupovat od konce řídícího horizontu a postupně hledat $J_t(x_t)$. Pro výpočet $x_{t+1}$ se použije rovnice \eqref{sys}. Libovolná řídící strategie $\pi=\{\mu_0,\ldots,\mu_{N-1}\}$, která splňuje systém rovnic \begin{equation} \label{impl} J_t(x_t)=\E_{w_t}\left\{g_k(x_t,\mu_t(x_t),w_t)+J_{t+1}(f_t(x_t,\mu_t(x_t),w_t))\right\}, \qquad t=0,\ldots,N-1, \end{equation} pak bude optimální posloupností rozhodnutí. \section{Úloha stochastického řízení s nepřesnými daty} Při aplikaci matematického modelování na řešení nějaké konkrétní úlohy se obvykle potýkáme s problémem jak určit konstanty, které daný model určují. Zkoumáme-li například nějaký fyzikální systém, z rozboru fyzikálních zákonitostí obvykle známe tvar rovnic, které určují jeho vývoj v čase. Nicméně počáteční podmínky či parametry, které v rovnicích vystupují a jsou pro daný systém charakteristické, můžeme získat pouze nepřímo, obvykle měřením vhodných veličin. Tento oddíl se zabývá modifikací úlohy stochastického řízení pro případ přítomnosti neznámých parametrů. \subsection{Výstup systému a informační vektor} Informace o stavu systému $x_t$ v čase $t$ získáváme pomocí výstupu $y_t$, který je dán jako \begin{equation} \label{poz} y_0=h_0(x_0,v_0),\qquad y_{t+1}=h_{t+1}(x_{t+1},u_t,v_{t+1}), \qquad t=1,\ldots,N-1, \end{equation} kde $v_t$ je náhodná veličina charakterizující chybu měření. Předpokládáme znalost funkcí $h_t$. Počáteční stav $x_0$ je dán rozdělením pravděpodobnosti $P^{x_0}$ a další vývoj systému určuje soustava \eqref{sys}. Informace, které jsou v průběhu řízení k dispozici je zvykem psát ve formě tzv. informačního vektoru, který má tvar \begin{equation} I_0=y_0,\qquad I_{t+1}=(y_{0:t+1},u_{0:t}), \qquad t=1,\ldots,N-1. \end{equation} \subsection{Optimální řízení pro úlohu s nepřesnými daty} Řídící zásah nyní nemůže explicitně záviset na stavu systému, protože máme k dispozici pouze informační vektor. Podobně jako v předešlé kapitole proto zavádíme neprázdnou množinu $U(I_t)$ všech přípustných řídících zásahů za informace $I_t$. Přípustnou řídící strategií $\pi=\mu_{0:N-1}$ bude posloupnost \begin{equation} \label{icon} \mu_t(I_t)=u_t, \qquad t=0,1,\ldots,N-1, \end{equation} kde $\mu_t(I_t)=u_t \in U(I_t)$ je přípustný řídící zásah. Úkolem je najít přípustnou strategii, která by minimalizovala očekávanou ztrátu \begin{equation} \label{ilos} J_\pi=\E_{\substack{x_0,\ w_{0:N-1},\\ v_{0:N}}}\left\{\sum_{t=0}^{N-1}g_t(x_{t+1},\mu_t(x_t))\right\}, \end{equation} za podmínek \eqref{sys} a \eqref{poz}. \subsection{Převod na úlohu s úplnými daty} Protože v čase $t$ nemáme k dispozici přímo stav systému $x_t$, ale pouze informační vektor $I_t$, nemůžeme použít postup z předchozí kapitoly. Před tím je potřeba úlohu vhodně transformovat. Za tímto účelem zapíšeme informační vektor ve tvaru \begin{equation} \label{nep} I_0=y_0,\qquad I_{t+1}=(I_t,u_t,y_{t+1}), \qquad t=1,\ldots,N-1. \end{equation} Na tuto rovnost můžeme pohlížet jako na rovnice systému \eqref{sys}. Stav v čase $t$ je nyní $I_t$, vstup $u_t$ a $y_{t+1}$ náhodná veličina podmíněná $I_t$ a $u_t$ přes \eqref{poz}. Dále přejdeme k nové ztrátové funkci, kterou definujeme jako \begin{equation} \tilde{g}_t(I_{t+1},u_t)=\E_{x_{t+1}}\left\{g_t(x_{t+1},u_t)|I_t,u_t\right\}, \qquad t=1,\ldots,N-1, \end{equation} kde $x_{t+1}$ se počítá dle \eqref{sys} a $x_t$ se považuje za náhodnou veličinu podmíněnou informačním vektorem $I_t$. Očekávanou ztrátu podproblému od času $t$ do $N$ nyní můžeme psát ve tvaru \begin{gather} J_N(I_N)=0,\\ J_t(I_t)=\min_{u_t \in U_t}\E_{w_t,y_{t+1}}\left\{\tilde{g}_t(I_{t+1},u_t)+J_{t+1}(I_{t+1})|I_t,u_t\right\}, \qquad t=0,\ldots,N-1. \end{gather} Tato úloha již může být řešena pomocí dynamického programování. Pří řešení budeme postupovat od konce řídícího horizontu a postupně hledat $J_t(I_t)$. Potom libovolná $\pi=\{\mu_0,\ldots,\mu_{N-1}\}$, která nabývá minimální očekávané ztráty $J_0(y_0)$ je optimální řídící strategie. \section{Úloha řízení systému s neznámými parametry} Pokud chceme řídit systém, jehož výstup závisí na nějakém neznámém konstantním parametru $\theta$, můžeme využít znalosti řešení problému s neúplným pozorováním. Parametr $\theta$ bude reprezentovat stav systému $x_t$, který se nyní v čase nemění. \subsection{Systém s neznámými parametry, hyperstav} V této úloze máme výstupy systému $y_t$ popsány jako \begin{equation} \label{poz2} y_0=h_0(\theta,v_0),\qquad y_{t+1}=h_t(I_t^{(d)},\theta,u_t,v_{t+1}), \qquad t=0,\ldots,N-1, \end{equation} kde $I_t^{(d)}=(y_{t:t-d},u_{t-1:t-d})$ a číslo $d$ se nazývá řád modelu. Označme $T_t$ dostatečnou statistiku pro parametr $\theta$ založenou na informacích dostupných v čase $t$. Pokud dostatečná statistika neexistuje, pak bude $T_t$ označovat nějakou její vhodnou aproximaci. Označme dále $H_t=(I_t^{(d)},T_t)$ tzv. hyperstav systému. Předpokládejme dále, že o parametru $\theta$ máme nějakou apriorní informaci v podobě hustoty pravděpodobnosti $f(\theta|T_0)$. Aposteriorní hustotu $f(\theta|T_{t+1})$ získáme pomocí Bayesova vzorce \begin{equation} \label{bay} f(\theta|T_{t+1}) = \frac{f(y_{t+1} | \theta, I_t^{(d)},u_t) f(\theta| T_t)} {\int f(y_{t+1} | \theta, I_t^{(d)},u_t) f(\theta| T_t)\mathrm{d}\theta}. \end{equation} Rekurzivní použití vzorce \eqref{bay} pro odhad parametru $\theta$ se nazývá postup Bayesovského učení \cite{peterka1981bayesian}. Pro vývoj hyperstavu $H_t$ v čase můžeme na základě \eqref{bay} psát \begin{equation} \label{the} H_{t+1}=f_t(H_t,u_t,y_{t+1}), \qquad t=1,\ldots,N-1. \end{equation} Rovnici \eqref{the} můžeme podobně jako \eqref{nep} považovat za rovnici systému \eqref{sys} pro stav $H_t$ a vstup $u_t$ s šumem $y_{t+1}$. \subsection{Převod na úlohu s nepřesnými daty} Ztrátová funkce je nyní \begin{equation} \label{los2} g(y_{1:N},u_{0:N-1})=\sum_{t=0}^{N-1}g_t(y_{t+1},u_t). \end{equation} Úlohou je nalezení řídící strategie $\pi=\mu_{0:N-1}$, která by minimalizovala očekávanou ztrátu \begin{equation} \label{ilos2} J_\pi=\E_{\theta_0,v_{0:N-1}}\left\{\sum_{t=0}^{N-1}g_t(y_{t+1},\mu_t(H_t))\right\}, \end{equation} za apriorní informace $f(\theta|T_0)$, známého rozdělení šumu $v_t$ a podmínek \eqref{the} a \eqref{poz2}. Rovnice \eqref{the}, \eqref{poz2} a \eqref{los2} potom představují úlohu stochastického řízení s nepřesnými daty. Úlohu opět řešíme pomocí dynamického programování, tedy postupnou minimalizací očekávané ztráty od konce řídícího horizontu \begin{gather} J_N(H_N)=0,\\ J_t(H_t)=\min_{u_t \in U_t}\E_{y_{t+1}}\left\{g_t(y_{t+1},u_t)+J_{t+1}(H_{t+1})|H_t,u_t\right\}, \qquad t=0,\ldots,N-1, \end{gather} kde $H_{t+1}$ se počítá dle \eqref{the}. Střední hodnota vzhledem k $y_{t+1}$ se počítá pomocí \eqref{poz2} a $f(\theta|T_t)$ jakožto aktuálního odhadu na parametr $\theta$. \subsection{Kalmanův filtr} Pokud v rovnicích \eqref{poz2} popisujících výstup systému vystupuje aditivní gaussovký šum a neznámý parametr je separován jako lineární člen, můžeme vypočítat konkrétní tvar rovnice \eqref{the}, tzv. Kalmanův filtr \cite{kalman1960new}. Dle předpokladu má výstup v čase $t$ tvar \begin{equation} \label{sys2} y_{t+1}=\tilde{h}_t(I_t,u_t)+A_t(I_t,u_t)\theta+v_{t+1}, \qquad t=0,\ldots,N-1. \end{equation} kde $\tilde{h}_t(I_t,u_t)$, resp. $A_t(I_t,u_t)$ je známá funkce, resp. matice závisící na informačním vektoru a aktuální vstupu. Dále předpokládáme gaussovské rozložení šumu $v_{t+1}$ se známým rozptylem \begin{equation} v_{t+1}\sim N(0,Q_{t+1}), \end{equation} gaussovské rozložení odhadu neznámého parametru $\theta_t$ a jejich nekorelovanost, tedy \begin{gather} \theta_t\sim N(\hat{\theta}_t,P_t),\\ \cov(v_{t+1},\theta_t)=0. \end{gather} Dosazením do \eqref{bay} se odvodí, že aposteriorní hustota pravděpodobnosti $f(\theta|T_{t+1})$ je rovněž gaussovská a její parametry $(\hat{\theta}_{t+1}, P_{t+1})$ splňují rovnice \begin{gather} K_t=P_tA_t(A_t^TP_tA_t+Q_t)^{-1},\\ \hat{\theta}_{t+1}=\hat{\theta}_t+K_t(y_{t+1}-\tilde{h}_t(I_t,u_t)-A_t\hat{\theta}_t),\\ P_{t+1}=(I-K_tA_t)P_t. \end{gather} Odvození lze nalézt v \cite{peterka1981bayesian}. Alternativní odvození bez požadavku gaussovského šumu je možné provést za předpokladu, že odhadovací proceduru střední hodnoty $\hat{\theta}_{t+1}$ neznámého parametru $\theta$ budeme hledat ve tvaru lineární opravy střední hodnoty $\hat{\theta}_t$ úměrné neurčitosti v systému. Tedy že \begin{equation} \label{opr} \hat{\theta}_{t+1}=\hat{\theta}_t+K_t(y_{t+1}-\E_{\theta,v_t}y_{t+1}), \end{equation} kde $K_t$ je neznámá matice, kterou určíme z požadavku minimalizace výsledné matice rozptylu $P_{t+1}$. Pro šum $v_t$ budeme požadovat nulovou střední hodnotu a existenci druhého momentu. Matici rozptylu označíme opět $Q_t$. Pro matici $P_{t+1}$ jako funkci $K_t$ můžeme psát \begin{equation} P_{t+1}(K_t)=\E[(\theta-\hat{\theta}_{t+1})(\theta-\hat{\theta}_{t+1})^T]. \end{equation} Dosazením za $\hat{\theta}_{t+1}$ z \eqref{opr} a za $y_t$ ze \eqref{sys2} a úpravou dostaneme (pro libovolnou matici $B$ budeme pro lepší čitelnost namísto $BB^T$ psát zkráceně $B^2$) \begin{align*} P_{t+1}(K_t)&=\E_{\theta,v_t}\left\{(\theta-\hat{\theta}_t-K_t(y_{t+1}-\tilde{h}_t(I_t,u_t)-A_t\hat{\theta}_t))^2\right\} \\ &=\E_{\theta,v_t}\left\{((I-K_tA_t)(\theta-\hat{\theta}_t)-K_tv_t)^2\right\} \\ &=(I-K_tA_t)\E\left\{(\theta-\hat{\theta_t})^2\right\}(I-K_tA_t)^T-(I-K_tA_t)\cov(\theta,v_t)K_t^T-\\ &-K_t\cov(\theta,v_t)(I-K_tA_t)^T+K_t\E\left\{v_t^2\right\}K_t^T. \end{align*} Použitím definice $P_t$, $Q_t$ a předpokladu $\cov(\theta,v_t)=0$ máme \begin{equation} \label{Pt+1} P_{t+1}(K_t)=(I-K_tA_t)P_t(I-K_tA_t)^T+K_tQ_tK_t^T. \end{equation} Protože požadujeme minimální rozptyl odhadu $\hat{\theta}_{t+1}$, určíme $K_t$ z rovnice \begin{equation} \frac{\partial \tr( P_t)}{\partial K_t}=0. \end{equation} K provedení derivace použijeme vzorce \begin{gather} \frac{\partial\tr(MXN)}{\partial X}=M^TN^T,\\ \frac{\partial\tr(MXNX^TO)}{\partial X}=M^TO^TXN+OMXN, \end{gather} kde $M,N$ a $O$ jsou konstantní matice. Tím získáme lineární rovnici pro $K_t$ tvaru \begin{equation} -P_t^TA_t-P_tA_t+K_tA_tP_tK_t+K_tA_t^TP_tK_t+2Q_tK_t=0, \end{equation} která má řešení \begin{equation} \label{Kt} K_t=P_tA_t(A_t^TP_tA_t+Q_t)^{-1}. \end{equation} Dosazením \eqref{Kt} do \eqref{Pt+1} po úpravě dostaneme \begin{equation} \label{Pt+12} P_{t+1}=(I-K_tA_t)P_t \end{equation} Rovnice \eqref{opr}, \eqref{Kt} a \eqref{Pt+12} představují rovnice Kalmanova filtru.