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émy: 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 řízení. Připomeňme, že úlohou ří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{Duální řízení}
K navržení kvalitního řízení systému s neznámými parametry je potřeba eliminovat neurčitost, která způsobuje nepředvídané chování systému. Hledané řízení by tedy mělo zároveň
\begin{itemize}
\item minimalizovat aktuální ztrátu,
\item získat o systému co nejvíce informací pro minimalizaci budoucích ztrát.
\end{itemize}
Tento postup, nazývaný duální řízení, navrhl Alexander Feldbaum v práci \cite{feldbaum1965optimal}. Název metody pochází z faktu, že požadavky, které na dobré duální řízení klademe, jsou obvykle v rozporu. Abychom se o systému něco nového dozvěděli, je vhodné ho vybudit do stavu, kde může vykazovat nepředvídané chování. Naopak pro minimalizaci aktuální ztráty potřebujeme systém uvést do stavu předem známého.

Jak konkrétně volit tvar optimálního řídícího $\mu_t(H_t)$, aby odpovídal principům duálního řízení, je obtížná  úloha. Často (například \cite{wittenmark2002adaptive}) se tak volí aproximace tvaru 
\begin{equation}
\label{dual}
\mu_t(H_t)=\mu_t^{(1)}(H_t)+\mu_t^{(2)}(H_t),
\end{equation}
kde $\mu_t^{(1)}$ je řídící zásah snažící se o minimalizaci aktuální ztráty a  $\mu_t^{(2)}$ je budící zásah, který slouží k excitaci systému za účelem identifikace neznámých parametrů.

Problém je nyní ve výrazu \eqref{dual} vhodně určit funkce $\mu_t^{(1)}$ a $\mu_t^{(2)}$. Vhodným výběrem na řídící složku $\mu_t^{(1)}$  duálního řízení je například řízení metodou opatrného řízení (v anglické literatuře "cautious control", zkráceně pak CC), nebo metodou certainty equivalent (CE). Obě metody budou probrány dále. Pro budící složku je obvykle obtížné určit nějaký funkcionální tvar. Jako vhodná aproximace se ukazuje použití bílého šumu, který se zapne, jen pokud neurčitost v systému překročí nějakou mezní hodnotu. Nevýhodou této přímočaré aproximace je nutnost určení amplitudy budícího signálu. K tomu obvykle nějakou systematickou metodu k dispozici nemáme.

\section{Certainty equivalent}
Při použití metody certainty equivalent (CE, například \cite{bertsekas1995dynamic}) se v rovnici pro očekávanou ztrátu \eqref{los3} náhodná veličina $y_{t+1}$ nahradí střední hodnotou $\hat{y}_{t+1}$. Ta se vypočítá z \eqref{poz4} pomocí známého rozdělení $v_t$ a postačující statistiky $T_t$ pro parametr $\theta$. Očekávaná ztráta \eqref{los3} 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}
kde $\hat{H}_{t+1}$ se vypočte pomocí \eqref{the2} použitím $\hat{y}_{t+1}$.

Tím odpadne výpočet střední hodnoty a zbývá pouze problém minimalizace vzhledem k $u_t$. Jedná se  tedy o jednodušší úlohu, kterou se již může povést vyřešit. Nutno podotknout, že optimální řízení získané touto metodou nemusí být optimální pro původní úlohu. Některé vlastnosti CE jsou probírány dále při aplikaci na řízení jednoduchého systému. Podrobnější pojednání lze nalézt v \cite{bertsekas1995dynamic}.

\section{Opatrné řízení}
Metoda opatrného řízení (anglicky cautious control, například \cite{wittenmark2002adaptive}) spočívá v optimalizaci očekávané ztráty \eqref{los3} na horizontu délky $N=1$. Minimalizujeme tedy
\begin{equation}
J_0(H_0)=\min_{u_0 \in U_0}\E_{y_1}\left\{g_0(y_1,u_0)|H_0,u_0\right\}.
\end{equation}
Stačí proto spočíst střední hodnotu $g_0(y_1,u_0)$ vzhledem k $y_1$ a výsledek minimalizovat vzhledem k $u_t$. Poznamenejme, že výsledné řízení nebude zcela jistě duální. To plyne z toho, že minimalizujeme očekávanou ztrátu pouze jeden krok dopředu a tedy se nemůže projevit výhoda identifikace parametrů pomocí vybuzení systému mimo požadovaný stav.

\section{Iterativní dynamické programování}
\label{idp}
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{los3}. Standardní numerický přístup k dynamickému programová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$. Při diskretizaci daného prostoru totiž počet použitých bodů roste exponenciálně s jeho dimenzí. Tato skutečnost, v anglické literatuře označovaná jako "curse of dimenzionality", pak prakticky znemožňuje řešení komplexnějších úloh.

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. Jak bude ještě diskutováno, výhodnost iterativního přístupu spočívá v jeho nižší 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 diskretizaci 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í. Při něm totiž 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čí ta část prostoru, která je pro bezprostřední výpočet podstatná. Díky tomu stačí k dostatečně jemné diskretizaci podstatně méně bodů. Iterativní dynamické programování tak problém nárůstu časové náročnosti se zvyšující se dimenzí hyperstavu neřeší, nicméně díky tomu, že stačí používat výrazně méně bodů k diskretizaci prostoru, ho alespoň částečně redukuje.

\subsection{Konvergence metody}
Lze ukázat, že za určitých, relativně obecných, předpokladů, algoritmus iterativního dynamického programování konverguje k optimálnímu řešení. Předpoklady a hlavní teorémy zaručující konvergenci  metody jsou shrnuty v tomto oddíle. Jejich znění je převzato z  \cite{thompson2005stochastic}. Důkazy je možno nalézt tamtéž.

Označme $U=\{U_1,U_2,\ldots,U_{N-1}\}$ množinu přípustných řídících zásahů. Dále buď $J_0(u)$ hodnota očekávané ztráty v závislosti na zvolené posloupnosti řídících zásahů $u=\{u_1,u_2,\ldots,u_{N-1}\}$ pro konkrétní hodnotu počátečního stavu. Počáteční rozsah pro hledání optimálního $u^*$ pomocí perturbací dosavadního napočteného řešení označíme $\beta^{in}$, faktor redukující počáteční rozsah v průběhu iterací IDP pak $\gamma$ ($\gamma<1$). Řešení napočtené v $j$-té iteraci algoritmu budeme značit $u^{*[j]}$.

\begin{predpoklady}
\label{predp}
Základní předpoklady pro použití metody IDP jsou
\begin{itemize} 
\item $U \in  \mathbb{R}^{N}$ je kompaktní,
\item $J_0(u)$ je ostře konvexní na $U$,
\item $J_0(u)$ je zdola omezená na $U$,
\item množina $L(u^{*[j]})=\{u \in U| J_0(u) \le  J_0(u^{*[j]}) \}$ je kompaktní $\forall u^{*[j]}$
\end{itemize}
\end{predpoklady}

\begin{veta}
Za předpokladů \ref{predp} existuje $u^*$ splňující $u^*=\argmin_{u \in U} J_0(u)$. Navíc $u^*$ je konečné a právě jedno.
\end{veta}

\begin{veta}
Množina bodů $u^* \in \{ \argmin_{u \in U}J_0(u)\}$ jsou pevné body iterační posloupnosti $\{u^{*[j]}\}_{j=1}^\infty$ generované pomocí IDP.
\end{veta}

\begin{veta}
Posloupnost iterací $\{u^{*[j]}\}_{j=1}^\infty$ generovaná pomocí IDP je vždy konvergentní. Navíc pro dostatečně velké $\gamma$ (tj. pro $\gamma \to 1 $) konverguje iterační posloupnost pro libovolné počáteční $u^{*[0]}$ k $u^*$ s pravděpodobností rovnou jedné.
\end{veta}

\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ři 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 z $n$ realizací náhodných veličin přes které je prováděna střední hodnota (jsou to $\theta_{t:N-1}$ a $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áší generovat náhodnou realizaci š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ě z \eqref{the2} 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 optimálního řízení. 

\subsection{Algoritmus SIDP}
\label{sidp1}
V tomto oddílu je popsán algoritmus SIDP, tak jak byl navržen v \cite{thompson2005stochastic}. Parametry algoritmu jsou uvedeny v tabulce \ref{par1}. Jak plyne z uvedeného schématu, časová složitost algoritmu SIDP vzhledem k jeho parametrům a délce řídícího horizontu 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{table}
\begin{tabular}{|l|l|} 
\hline
počet opakování algoritmu & $n_{pass}$ \\
počet iterací algoritmu & $n_{iter}$\\
počet bodů v diskretizaci každé dimenze $H_t$ & $n_g$\\
apriorní řídící strategie & $\mu_{0:N-1}$\\
počet kandidátů na změnu  řídícího zásahu & $m$\\
počáteční rozsah pro hledání optimálního řídícího zásahu & $\beta^{in}$\\
parametr pro redukci $\beta^{in}$ při opakování algoritmu &$\gamma$\\
parametr pro redukci $\beta^{in}$ při iterování algoritmu  & $\lambda$\\
počet realizací pro odhad metodou Monte Carlo & $n$\\\hline
\end{tabular}
\caption{Parametry algoritmu SIDP}
\label{par1}
\end{table}

\begin{algorithm}
\begin{algorithmic}
\STATE {\large\textbf{Schéma algoritmu SIDP}}
\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}$, k tomu použij napočtenou řídící strategii $\pi^*$, její interpolace/extrapolace a realizace šumu a 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íť kolem bodů v $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}
\label{alg}
\end{algorithm}

\subsection{Detaily algoritmu}
\label{sidp}
Část prostoru, která se bude v následující iteraci algoritmu diskretizovat, se určí pomocí aktuálního suboptimálního řešení (jeho interpolace a extrapolace) 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í. Schematicky je situace znázorněna na obrázku \ref{tra}. Jak danou část prostoru diskretizovat je otázka implementace algoritmu na konkrétní systém, bude tedy probrána později.

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

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. Řídící zásah, pro který bude průměr přes realizace ztráty nejnižší, se ponechá do další iterace.

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 popsaný v článku  \cite{nelson2001simple}. V první úrovni tohoto algoritmu se pro každého kandidáta $u_t$ vygeneruje $n_0$ realizací. Na jejich základě se vyberou ti, na kterých 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.

Výstupem algoritmu je tabulka optimálních řídících zásahů bodech diskretizace. Samotné řízení systému je pak prováděno aplikací předpočtených optimálních zásahů.