Ačkoliv použití dynamického programování přináší významný pokrok v řešení úlohy stochastického řízení, analytické řešení obvykle není možné získat. V každém časovém kroku se totiž potýkáme se dvěma obecně obtížnými problémemy: 1) výpočet střední hodnoty a 2) minimalizace vzhledem k $u_t$. Oba problémy obecně nemají analytické řešení a bez další specifikace úlohy je proto třeba přejít k aproximačním metodám.

V této kapitole se předkládá popis několika možných přístupů k aproximativnímu řešení úlohy duálního řízení. Přípomeňme, že úlohou duálního řízení  je nalezení řídící strategie $\pi=\mu_{0:N-1}$, která by minimalizovala očekávanou ztrátu
\begin{equation}
\label{ilos3}
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 $\theta_0$ a podmínek
\begin{gather}
\label{the2}
H_{t+1}=f_t(H_t,u_t,y_{t+1}),\\
\label{poz4}
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{gather}
kde $H_t=(I_t^{(d)},T_t)$ je hyperstav systému a $T_t$ dostatečná statistika pro neznámý parametr $\theta$ v čase $t$.

Úlohu řešíme pomocí dynamického programování, tedy postupnou minimalizací očekávané ztráty od konce řídícího horizontu
\begin{equation}
\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{equation}
kde $T_{t+1}$ a $y_{t+1}$ se počítá dle \eqref{the2} a \eqref{poz4}.

\section{Certainty equivalent control}
Při použití metody Certainty equivalent control (CEC) se v rovnici pro očekávanou ztrátu nahradí náhodná veličina $y_{t+1}$ střední hodnotou $\hat{y}_{t+1}$. Ta se vypočítá z \eqref{poz4} pomocí známých rozdělení na $v_t$ a postačující statistiky $T_t$. Očekávaná ztráta \eqref{los} tak přejde v
\begin{gather}
\label{CE}
J_t(H_t)=\min_{u_t \in U_t}\left\{g_t(\hat{y}_t,u_t)+J_{t+1}(\hat{H}_{t+1})|H_t,u_t\right\}, \qquad t=0,\ldots,N-1,
\end{gather}

Podrobnější pojednání s diskuzí aspektů použití CEC lze nalézt v \cite{bertsekas1995dynamic}.

\section{Metoda separace}
Při použití metody separace je proces řízení rozdělen do dvou fází: 1) indentifikace neznámého parametru a 2) řízení za použití odhadu $\hat{\theta}$ z první fáze.

První fáze slouží k nezávislému sběru dat, která jsou následně použita k odhadu neznámého parametru. K odhadu můžeme použít například rovnici \eqref{the2}. V druhé fázi pak po zbytek řídícího horizontu použijeme  pro návrh řídící strategie odhad $\hat{\theta}$ z první fáze.

\section{Duální řízení}
Hledané řízení by mělo nejen minimalizovat aktuální ztrátu, ale rovněž získat o systému co nejvíce informací pro minimalizaci budoucích ztrát. Tento postup se nazývá duální řízení [ref]. ODKAZ NA FILDEBAUMA, POPIS PRINCIPU... (napr JEDNOKROKOVA OPTIMALIZACE S BUZENIM - FILATOV)

\section{Metoda Monte Carlo}
Metoda Monte Carlo \cite{hammersley1964monte} je statistická simulační metoda. Její princip spočívá ve vzorkování nějaké náhodné veličiny za účelem odhadu její hledané charakteristiky, např. střední hodnoty. V této práci je metoda Monte Carlo použita k výpočtu očekávané ztráty \eqref{los3}. 

Pří běžném použití dynamického programování máme při výpočtu $J_t(H_t)$ k dispozici předpis pro následující očekávanou ztrátu $J_{t+1}(H_{t+1})$. Metoda monte Carlo nám však dá k dispozici pouze odhad očekávané ztráty a použití těchto aproximací v dalším výpočtu by chybu výpočtu  navyšovalo. Namísto toho se pro další výpočet uchovávají $\mu_t(H_t)$ a očekávaná ztráta v čase $t$ se pak počítá jako průměr přes $n$ realizací náhodných veličiny přes které je prováděna střední hodnota $(\theta_{t:N-1},v_{t:N})$, tedy
\begin{equation}
\label{mon}
\frac{1}{n}\sum_{i=1}^n\left(\sum_{j=t}^{N-1}g_j(y_{j+1}^i,\mu_j(H_j^i)\right),
\end{equation}
kde $y_{j+1}^i$ se počítá podle \eqref{poz4} jako
\begin{equation}
y_{j+1}^i=h_j( I_j^i,\theta_j^i,\mu(H_j^i),v_{j+1}^i), \qquad j=t,\ldots,N-1, \qquad i=1\ldots,n,
\end{equation}
a index $i$ označuje $i$-tou realizaci dané veličiny. Realizace $\theta_{t:N-1}$ se generují podél trajektorie \eqref{poz4}. To znamená, že dané  $\theta_{k+1}$ se generuje až ve chvíli, kdy je známé $I_k$, $u_k$, postačující statistika $T_k$ a $y_{k+1}$ a tedy přes \eqref{the2} i hustota pravděpodobnosti $f(\theta_{k+1})$.

Tento jednoduchý postup lze vylepšit víceúrovňovým porovnáním kandidátů na optimální řízení. Jedním z možných vylepšeních je dvouúrovňový algoritmus poposaný v \cite{nelson2001simple}. V první fázi tohoto algoritmu se nejprve pro každého kandidáta vygeneruje $n_0$ realizací. Na jejich základě se vyberou ti, na který je nabyto minima s pravděpodobností větší než je daná mez $\alpha_0$. Pro tyto se v druhé fázi vygeneruje dostatečný počet realizací tak, aby bylo možné nejlepší rozhodnutí zvolit s pravděpodobností alespoň rovné zadané mezi $\alpha_1$. Takto upravený algoritmus metody Monte Carlo je robustnější a umožňuje porovnání většího množství kandidátů, neboť počet realizací v první fázi může být poměrně nízký, slouží pouze k odfiltrování zjevně horších kandidátů na řízení. 

\section{Iterativní dynamické programování}
Iterativní dynamické programování \cite{luus2000iterative} je jední z přístupů k nalezení optimální strategie, která minimalizuje očekávanou ztrátu \eqref{ilos2}. Oproti dynamickému programování se problém řeší iterativně. Na začátku se zvolí nějaká apriorní strategie. V každé iteraci se potom vychází ze strategie spočtené v předchozím kroku a prostřednictvím perturbací tohoto (suboptimálního) řešení se hledá strategie, pro kterou bude očkávaná ztráta nižší. Tato se použije v následující iteraci.

\subsection{Diskretizace prostoru}
Při hledání optimální strategie $\mu_t(H_t)$ bychom pro přesné vyčíslení očekávané ztráty \eqref{mon} na úseku řídícího horizontu $t\!:\!N$ potřebovali její analytické vyjádření. To ale není obvykle možné. Je proto nutné přejít k nějaké aproximaci, například 1) předpokládat nějaký tvar optimální strategie a při výpočtu určit pouze konstanty, které výslednou strategii určí jednoznačně, nebo 2) diskretizovat prostor $(H_t)$ a počítat $\mu_t(H_t)$ jen v bodech diskretizace a jinde se uchýlit k interpolaci (popřípadě extrapolaci). 

Jakým způsobem efektivně diskretizovat prostor nezávislých proměnných pro aproximativní výpočet očekávané ztráty \eqref{mon} je při použití dynamického programování obtížná otázka. Bude-li bodů v diskretizaci příliš málo, bude výpočet nespolehlivý, naopak pro příliš jemnou diskretizaci bude časová náročnost výpočtu rychle stoupat (o časové náročnosti SIDP viz dále). Zde se ukazuje výhodnost použití iterativního dynamického programování, neboť stačí diskretizovat jen tu část prostoru která bude potřebná v následující iteraci. Pomocí strategie spočtené v předchozím kroku a náhodných realizací šumu $v_{0:N}$ a neznámého parametru $\theta_{0:N}$ vygenerujeme trajektorie v $(H)_{0:N}$. V každé časové úrovni pak diskretizujeme jen tu část prostoru, která byla zasažena. 

V této práci je volena jednoduchá metoda v které se spočte nejmenší hyperkvádr kolem zasažené tak, že se vezme nejmenší hyperkvádr orientovaný ve směru souřadných os, do kterého se vygenerované body vejdou. Prostor se poté diskretizuje pouze v této oblasti. Metodu k určení hyperkvádru s obecnou orientací lze najít v \cite{bh-eamvb-01}. 

\section{SIDP}
Metoda stochastického iterativního dynamického programování (SIDP) \cite{thompson2005stochastic} spočívá v současném použití  metody Monte Carlo k získání aproximace pro očekávanou ztrátu  a iterativního dynamického programování k nalezení optimální strategie. Pro účely této postačuje základní verze metody Monte Carlo a je proto v následující implementaci SIDP použita. Při použití iterativního dynamického programování se uchýlíme k diskretizovat prostoru hyperstavů a budeme používat interpolaci (popřípadě extrapolaci) napočtených hodnot. Poznamenejme, že díky předpokladu gaussovského rozdělení parametru ${\theta_t}$, diskretizace vzhledem k ${T_t}$ znamená diskretizaci vzhledem k ${(\hat{\theta}_t,P_t)}$.

\subsection{Algoritmus SIDP}
V tomto odílu je schématicky popsán algoritmus SIDP. Jeho parametry jsou

\begin{itemize}
\item $n_{pass}, \, n_{iter}$– počet opakování a iterací algoritmu
\item $N$ -- řídící horizont
\item $n_g$ -- počet bodů v diskretizaci každé dimenzi $H_t$, tj. $|H_t|=n_g^{\dim H_t}$
\item $\pi^*=\mu_{0:N-1}(H_{0:N-1})$ -- apriorní řídící strategie
\item $m$ -- počet kadnidátů na změnu  řídícího zásahu v jedné iteraci IDP
\item $\beta^{in}$ -- počáteční rozsah pro hledání optimálního řídícího zásahu
\item $\gamma,\, \lambda$ -- parametry pro redukci $\beta^{in}$
\item $n$ -- počet realizací pro odhad metodou Monte Carlo
\end{itemize}

Jak plyne z následujícího popisu, časová složitost SIDP vzhledem k jeho parametrům je $O(n_{pass}n_{iter}N^2mnn_g^{\dim H_N})$ (časová náročnost metody Monte Carlo je úměrná vzdálenosti od konce horizontu, proto je časová složitost úměrná druhé mocnině $N$).

\begin{algorithm}
\begin{algorithmic}
\FOR{$i = 1$ to $n_{pass}$} 
\FOR{$j = 1$ to $n_{iter}$} 
\STATE $\beta_{i,j} := \gamma^{j-1}\lambda^{i-1}\beta^{in}$
\FOR{$k = 1 $ to $|H_t|$}  
\STATE spočti trajektorii $H_{0,k}$, použij aktuální $\pi^*$, její interpolace a extrapolace a realizace neznámého parametru $\theta_0,\ldots,\theta_{N-1}$ podél této trajektorie
\ENDFOR
\FOR{$t = N-1 $ to $0$} 
\STATE vytvoř $\tilde{H}_t$ jakožto rovnoměrnou síť v oblasti bodů $H_t$
\STATE interpoluj (extrapoluj) $\mu_t^*(H_t)$ na $\mu_t^*(\tilde{H}_t)$
\FOR{$k = 1 $ to $|H_t|$}  
\FOR{$m=-\left[\frac{m-1}{2}\right]$ to $\left[\frac{m}{2}\right]$}
\STATE pro $\tilde{H}_{t,k}$ vygeneruj kandidáta na řízení $\mu_t(\tilde{H}_{t,k}) =  \mu_t^*(\tilde{H}_{t,k})+m\beta_{i,j}$ 
\STATE pomocí metody Monte Carlo spočti očekávanou ztrátu
\ENDFOR
\STATE rozhodnutí s nejnižší  očekávanou ztrátou uchovej jako nové optimální rozhodnutí pro $\tilde{H}_{t,k}$.
\ENDFOR
\ENDFOR
\ENDFOR
\ENDFOR
\end{algorithmic}
\end{algorithm}