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. Tato kapitola se zabývá modifikací úlohy stochastického řízení pro případ přítomnosti neznámých parametrů.

\section{Formulace úlohy stochastického řízení s nepřesnými daty}
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í. 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. \emph{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}

Ří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 množinu  $U(I_t)$ všech přípustných řídících zásahů za informace $I_t$ a přípustnou řídící strategií bude $\pi=\mu_{0:N-1}$ 
\begin{equation}
\label{icon}
\mu_t(I_t)=u_t \, \qquad t=0,1,\ldots,N-1,
\end{equation}
kde $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-1}}}\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}.

\section{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 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{Ří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 konstantí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í.

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}$. 

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 ř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\\
\label{los3}
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ání č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ím derivace použijeme vzorce*ODVOZENI BUDE ASI AZ V DODATKU*
\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.