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{Iterativní dynamické programování}

Iterativní dynamické programování \cite{luus2000iterative} je jednou z variant klasického přístupu k  nalezení optimální strategie, která minimalizuje očekávanou ztrátu \eqref{ilos2}. Standardní numerický přístup k dynamickému progamování lze shrnout následovně
\begin{enumerate}
\item prostor proměnných $H_t$ se diskretizuje do mřížky,
\item postupně se od konce horizontu napočítává minimální očekávaná ztráta $J_t(H_t)$ pro každý bod diskretizace $H_t$. K výpočtu se používají již napočtené minimální očekávané ztráty v následujících časech,
\item optimální strategie bude ta, na které bude nabyto minimální očekávané ztráty z počátečního stavu na konec řídícího horizontu.
\end{enumerate}
Tento postup je přímočarou aplikací principu dynamického programování. Bohužel je velmi citlivý na dimenzi stavového prostoru $H_t$, který je potřeba diskretizovat, neboť počet bodů potřebných k disretizaci roste exponenciálně s dimenzí prostoru. Tato skutečnost se v anglické literatuře označuje jako "curse of dimenzionality".

Oproti klasickému dynamickému programování iterativní dynamické programování problém řeší v sérii iterací. V každé iteraci se vychází ze strategie spočtené v předchozím běhu a prostřednictvím perturbací tohoto (suboptimálního) řešení se hledá strategie, pro kterou bude očekávaná ztráta nižší. Tato se použije v následující iteraci. Výhodnost iterativního přístupu spočívá ve snížení citlivosti na dimenzi úlohy.

\subsection{Diskretizace prostoru}
Při hledání optimální strategie $\mu_t(H_t)$ je pro přesné vyčíslení očekávané ztráty \eqref{mon} na úseku řídícího horizontu $t\!:\!N$ nutné znát její analytické vyjádření. To ale není obvykle možné. Je proto nutné přejít k nějaké aproximaci, například 
\begin{enumerate}
\item 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ě,
\item 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). 
\end{enumerate}

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 diskretizace příliš málo, bude výpočet nespolehlivý, naopak pro příliš jemnou diskretizaci bude počet bodů v diskretizaci hyperstavu rychle stoupat a časová náročnost výpočtu pak prakticky znemožní jeho řešení. 

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í perturbací strategie spočtené v předchozím kroku se určí část prostoru, která je pro bezprostřední výpočet podstatná. Díky tomu stačí k dostatečně jemné diskretizaci podstatně méně bodů. Konkrétní implementace bude probrána dále.

KONVERGENCE
\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}. 

\subsection{Použití metody Monte Carlo k výpočtu očekávané ztráty}
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 by nám však dala k dispozici pouze odhad očekávané ztráty. Použití těchto aproximací v dalším výpočtu by chybu výpočtu  navyšovalo.

Namísto $J_t(H_t)$ je proto vhodné pro další výpočet uchovávat $\mu_t(H_t)$. Očekávanou ztrátu v čase $t$ pak lze počítat jako průměr $n$ realizací náhodných veličin 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(g_j(y_{j+1}^i,\mu_j(H_j)+\sum_{j=t+1}^{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 $\theta_{j+1}^i$ se generuje až ve chvíli, kdy je známé $I_j^i$, $u_j^i$, postačující statistika $T_j^i$ a $y_{j+1}^i$ a tedy přes \eqref{the2} i postačující statistika $T_{j+1}^i$. 

Výpočet je při uchovávání $\mu_t(H_t)$ namísto $J_t(H_t)$ časově náročnější. Namísto přečtení $J_t(H_t)$ je totiž nutné přečíst hodnotu $\mu_t(H_t)$ a následně vygenerovat trajektorii od času $t$ do konce horizontu. To obnáší vygenerovat generovat náhodnou realizace šumu $v_t$ a neznámého parametru $\theta$ (pomocí $T_t$), aplikovat řídící zásah, tedy dle \eqref{poz4} vypočítat $y_{t+1}$ a následně \eqref{the2} dle vypočítat $T_{t+1}$. Tím bude určen bod v $H_{t+1}$. Zde pak pomocí interpolace (a extrapolace) určíme optimální zásah, který aplikujeme. Podobně jako prve tak určíme následující bod v $H_{t+2}$, až nakonec se dostaneme na konec řídícího horizontu. Při výpočtu postupně napočítáváme hodnotu ztrátové funkce.

\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. Při použití iterativního dynamického programování se uchýlíme k diskretizaci prostoru hyperstavů a budeme používat interpolaci (popřípadě extrapolaci) napočtených hodnot. Poznamenejme, že díky předpokladu gaussovského rozdělení odhadu neznámého parametru $\theta$, diskretizace vzhledem k ${T_t}$ znamená diskretizaci vzhledem k ${(\hat{\theta}_t,P_t)}$.

\subsection{Algoritmus SIDP}
V tomto odílu je popsán algoritmus SIDP, tak jak byl navržen v \cite{thompson2005stochastic}. Parametry algoritmu 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$, počet bodů v diskretizaci je tedy $|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}

\subsection{Detaily implementace}
Část prosotoru, která se bude v následující iteraci algoritmu diskretizovat se určí pomocí aktuálního suboptimálního řešení a náhodných realizací šumu $v_{0:N}$ a neznámého parametru $\theta_{0:N}$. Pomocí těchto realizací vygenerujeme trajektorie v $H_{0:N}$. V každé časové úrovni pak diskretizujeme jen tu část prostoru, kterou takto vygenerované trajektorie prochází. Schématicky je situace znázorněna v \ref{tra}

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{tra}
\caption{Trajektorie v hyperstavu -- jednotlivá realizace trajektorie je napočítávána pomocí realizací šumu a neznámého parametru}
\label{tra}
\end{figure}


V této práci se pro diskretizaci zasažené části prostoru volí jednoduchá metoda, ve které se ve směru souřadných os spočte nejmenší hyperkvádr obahující vygenerované body. Prostor se poté diskretizuje pouze v této oblasti. V práci \cite{thompson2005stochastic}, kde je metoda SIDP navržena, je pro diskretizaci prostoru použit hyperkvádr s obecnou orientací. Metodu k jeho určení převzali autoři z  \cite{bh-eamvb-01}. Tento postup by měl vést k ještě efektivnější diskretizaci prostoru. Nicméně metoda, která je v naší práci použita se ukázala jako postačující. Navíc je implementačně podstatně jednodušší a vyhledávání v tabulce s orientací ve směru souřadných je rychlejší. Výhodou je i to, že můžeme snadno zaručit požadavek na kladný rozptyl $P_t$ neznámého parametru $\theta$, viz \ref{box} . 

\begin{figure}
\centering
\includegraphics[width=0.25\textwidth]{box}
\caption{Oblast určená k diskretizaci $H_t$ -- ačkoliv je objem obecně orientovaného hyperkvádru menší, body v něm nesplňují požadavek na kladný rozptyl $P$}
\label{box}
\end{figure}

Máme-li diskretizovanou požadovanou část prostoru, je nutné na ni namapovat dosavadní napočtené optimální řízení. K tomu se použije interpolace, popřípadě extrapolace napočteného řešení. V této práci je interpolace/extrapolace realizována jednoduše pomocí nejbližšího již napočteného bodu. Možným vylepšením by byla například lineární projekce či vážený průměr nejbližších napočtených bodů.

Pro každý z bodů (nejprve pro ty na konci řídícího horizontu) se optimální řídící zásah hledá pomocí perturbace stávajícího suboptimálního řešení. Pro daný bod se proto vygeneruje $m$ kandidátů na optimální zásah, rovnoměrně kolem optimálního zásahu z předcházející iterace. Jako jeden z kandidátů na optimální řízení se vždy ponechá stávající suboptimální řešení z minulé iterace.

Kandidáti na řízení se nyní porovnají pomocí metody Monte Carlo. Jak již bylo popsáno výše, pro každého kandidáta se vygeneruje $n$ realizací ztráty, přes které se dle \eqref{mon} spočte průměr. 

Namísto jednoduchého porovnání pomocí průměru lze kandidáty na optimální řídící zásah porovnat nějakým sofistikovanějším víceúrovňovým algoritmem. Jedno z možných vylepšení je použito i v \cite{thompson2005stochastic}. Konkrétně se jedná o dvouúrovňový algoritmus poposaný v \cite{nelson2001simple}. V první úrovni tohoto algoritmu se nejprve pro každého kandidáta $u_t$ 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ší. Navíc umožňuje efektivní 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í. Pro účely této práce postačuje základní verze metody Monte Carlo a proto je v následující implementaci SIDP použita.