Ačkoliv použití dynamického programování přináší významný pokrok v řešení úlohy duální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{ilos2} J_\pi=\E_{\theta_0,v_{0:N-1}}\left\{g_N(y_N)+\sum_{t=0}^{N-1}g_t(y_t,\mu_t(I_t,\theta_t),v_t)\right\}, \end{equation} za apriorní informace $\theta_0$ a podmínek \begin{gather} \label{the2} \theta_{t+1}=f_t(I_t,\theta_t,u_t,y_{t+1}),\\ \label{poz4} y_0=h_0(\theta_0,v_0),\qquad y_{t+1}=h_t(I_t,\theta_t,u_t,v_{t+1}), \qquad t=0,\ldots,N-1.\\ v_{t+1}\sim N(0,Q_{t+1})\\ \theta_t\sim N(\hat{\theta},P_t),\\ \cov(v_{t+1},\theta)=0. \end{gather} Ú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(I_N,\theta_N)=\E_{\theta_N,v_N}\left\{g_N(y_N)\right\},\\ \label{los3} J_t(I_t,\theta_t)=\min_{u_t \in U_t}\E_{y_{t+1},v_t}\left\{g_t(y_t,u_t,v_t)+J_{t+1}((I_t, ,u_t,y_{t+1},\theta_{t+1}))|I_t,\theta_t,u_t\right\}, \\ \qquad t=0,\ldots,N-1, \end{gather} kde $\theta_{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) [ref] se v rovnici pro očekávanou ztrátu nahradí náhodné veličiny svými středními hodnotami. Očekávaná ztráta tak přejde v \begin{gather} J_N(I_N, \theta_N)=g_N(y_N),\\ J_t(I_t, \theta_t)=\min_{u_t \in U_t}\left\{g_t(y_t,u_t,\hat{v}_t) +J_{t+1}(I_t,\theta_{t+1},u_t,\hat{y}_{t+1}))|I_t,\theta_t,u_t\right\}, \\ \qquad t=0,\ldots,N-1, \end{gather} \section{Metoda separace} Při použití metody separace [ref] 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{SIDP} Metoda stochastického iterativního dynamického programování (SIDP) [ref] 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. \subsection{Metoda Monte Carlo} Metoda Monte Carlo je statistická simulační metoda, kterou navrhl ... [ref]. 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{ilos2}. Pří běžném použití dynamického programování máme při výpočtu $J_t(I_t,\theta_t)$ k dispozici předpis pro následující očekávanou ztrátu $J_{t+1}(I_{t+1},\theta_{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(I_t,\theta_t)$ a očekávaná ztráta v čase $t$ se pak počítá jako průměr přes $n$ realizací náhodné veličiny $(\theta_{t:N-1},v_{t:N})$, tedy \begin{equation} \label{mon} \frac{1}{n}\sum_{i=1}^n\left(g_N(y_N^i)+\sum_{j=t}^{N-1}g_j(y_j^i,\mu_j(I_j^i,\theta_j),v_j^i)\right), \end{equation} kde $y_{j+1}^i$ se počítá podle \eqref{poz3} jako \begin{equation} y_{j+1}^i=h_j( I_j^i,\theta_j^i,\mu(I_j^i, \theta_j),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$, rozdělení $\theta_k$ a $y_{k+1}$ a tedy přes \eqref{the2} i rozdělení $\theta_{k+1}$. Tento jednoduchý postup lze vylepšit víceúrovňovým porovnáním. Jedním z možných vylepšeních je dvouúrovňový algoritmus poposaným [ref]. V tomto algoritmu se nejprve pro každého kandidáta vygeneruje $n_0$ realizací. Na základě realizací se vyberou ti kandidá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 upraveny 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í. Pro účely této práce nicméně postačuje základní verze metody Monte Carlo a je proto v následující implementaci SIDP použita. \subsection{Iterativní dynamické programování} Iterativní dynamické programování [ref] 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(I_t,\theta_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 $(I_t,\theta_t)$ a počítat $\mu_t(I_t,\theta_t)$ jen v bodech diskretizace a jinde se uchýlit k interpolaci (popřípadě extrapolaci). V této práci volíme druhou zmíněnou metodu. Poznamenejme, že díky předpokladu gaussovského rozdělení parametru ${\theta_t}$, diskretizace vyhledem k ${\theta_t}$ znamená diskretizaci vzhledem k ${(\hat{\theta}_t,P_t)}$. 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 $(I,\theta)_{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 [ref]. \subsection{Algoritmus SIDP} V tomto odílu je 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^2mn_g^{\dim H_N})$ (časová náročnost metody Monte Carlo je úměrná vzdálenosti od konce horizontu). \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}