\documentclass{article}
\parskip=7pt   % mezera mezi odstavci

\usepackage[czech]{babel}      % česky psaná práce
\usepackage[T1]{fontenc}
\usepackage[cp1250]{inputenc}  % vstupní znaková sada: Windows 1250

\usepackage{amsmath}
\usepackage{epsfig}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{dcolumn}
\usepackage{amsfonts}

\DeclareMathOperator*{\E}{E}
\DeclareMathOperator*{\cov}{Cov}
\DeclareMathOperator*{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}

\begin{document}

\section{Popis úlohy}
Chceme řídit systém s výstupem $y_t$ popsaným jako 
\begin{equation}
\label{poz2}
y_0=h_0(\theta,v_0),\qquad y_{t+1}=h_t(I_t,\theta,u_t,v_{t+1}), \qquad t=0,\ldots,N-1,
\end{equation}
kde $I_t=(y_{0:t},u_{0:t-1})$.

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,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,u_t) f(\theta| T_t)}
{\int f(y_{t+1} | \theta, I_t,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í.

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=0,\ldots,N-1.
\end{equation} 

Ztrátová funkce je 
\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}.

Úloha je řešitelná 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,\\
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,
\label{los3}
\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$.

Pro konkrétní systém však vlastní výpočet obvykle není možné provést. Důvodem je obtížný výpočet parametrického integrálu (parametr je $u_t$) a jeho následná minimalizace vzhledem k $u_t$. V další sekci je proto představen aproximativní algoritmus k nalezení suboptimálního řízení.

\section{Algortimus SIDP se separovanou budící složkou (SIDPS)}
\label{sidp'}
Algoritmus je založený na SIDP, tj. na aproximaci dynamického programování pomocí iterování suboptimálního řízení a metody Monte Carlo.

Oproti klasickému SIDP specializujeme tvar hledaného řízení. To hledáme v tvaru
\begin{equation}
\label{rozklad}
u_t(H_t)=u^{(1)}_t(H_t)+u^{(2)}_t(\tilde{H}_t), \qquad  t=0,\ldots,N-1.
\end{equation}
kde $u^{(1)}_t$ je zásah starající se o minimalizaci ztráty a $u^{(2)}_t$ budí systém za účelem jeho lepší identifikace. Tvar $u^{(1)}_t$ předpokládáme za známý, předmětem optimalizace je $u^{(2)}_t$, který závisí pouze na $\tilde{H}_t \subset H_t$. Část $H_t$, na které $u^{(2)}_t$ nezávisí označíme $\tilde{H}^{(c)}_t$. Tedy je $H_t=\{{\tilde{H}_t,\tilde{H}^{(c)}_t}\}$.

Další změnou oproti SIDP je minimalizace ztrátové funkce na ubíhajícím horizontu délky  $\tilde{N}$ (obvykle volíme $\tilde{N}<<N$). Namísto iterování řízení na mřížkách $H_{0:N-1}$ za účelem nalezení optimálního $u_t$ pro každý bod $H_t$ (případ SIDP) tak iterujeme pouze mřížku obsahující $\tilde{H}_0$ a v každém bodě optimalizujeme jen $u^{(2)}_0$. Optimální zásah pro $j$-tý bod diskretizovaného $H_0$ ($=H_{0,j}$) označíme $u_{0,j}$.  Připoměňme, že optimalizace probíhá na ubíhajícím horizontu délky  $\tilde{N}$. Vzhledem k tomu, že v každém časovém okamžiku používáme řízení napočtené na jediné mřižce $\tilde{H}_0$, ztrácí se výhodnost lokální diskretizace hyperstavu (viz článek o SIDP). To vyplývá z předpokladu, že v průběhu řízení se bude neurčitost v systému vlivem procesu učení  měnit (i výrazně). Bude tedy třeba šírší spektrum hodnot $\tilde{H}_0$ pro které se bude hledat optimální řázení. V případě jediné mřížky bychom tak měli diskretizovat $\tilde{H}_0$ v celém rozsahu předpokládaných hodnot, pro které chceme systém řídit.

Schematickou podobu SIDPS zachycuje \ref{alg}.

\begin{algorithm}
\begin{algorithmic}
\STATE diskretizuj $\tilde{H}_0$
\STATE v bodech diskretizace $\tilde{H}_0$ (je jich $|\tilde{H}_0|$) inicializuj $\pi_1$ (např. $\equiv 0$)
\STATE zvol libovolně $\tilde{H}^{(c)}_0$ (dle předpokladu na této volbě $u^{(2)}_0$ nezávisí)
\FOR[iteruj $\pi_i=\{u_{0,1:|\tilde{H}_0|}\}$]{$i= 1$ to $n_{iter}$}
\FOR[pro $\tilde{H}_{0,j}$ najdi (lokálně) nejlepší řízení]{$j = 1 $ to $|\tilde{H}_0|$}  
\STATE perturbací $u^{(2)}_{0,j}$ vygeneruj $m$ kandidátů na řízení (ozn. $u_{1:m}$).

\FOR[spočti ztrátu $J^*$ pro $\tilde{\pi}_i=\{\tilde{u}_{0,1:k-1},u_k,u_{k+1:|\tilde{H}_0|}\}$]{$k = 1 $ to $m$}  
\FOR[$J^*$ aproximuj jako průměr z $n_0$ realizací  $J^{(l)}$]{$l = 1 $ to $n_0$}
\STATE $J_0^{(l)}=0$
\FOR[generuj realizaci ztráty $J_l$]{$t = 0 $ to $\tilde{N}-1$}
\STATE generuj $\theta$ a  $v_t$ (použijeme $f(\theta,T_t)$ a známé rozdělení šumu)
\STATE interpolací/extrapolací $\tilde{\pi}_i$ urči zásah $u^{(2)}_t$ pro $\tilde{H}_t$
\STATE spočti $y_{t+1}=h_t(I_t,u_t,v_{t+1},\theta)$ 
\STATE přičti aditivní část ztráty, tj. $J_{t+1}^{(l)}=J_{t}^{(l)}+g_t(y_{t+1},u_t)$
\STATE dopočítej $H_{t+1}=f_t(H_t,u_t,y_{t+1})$
\ENDFOR
\STATE $J^{(l)}=J_{\tilde{N}}^{(l)}$
\ENDFOR
\STATE $\bar{J}=\frac{1}{n_0}\sum_{l=1}^{n_0}J_l$
\ENDFOR
\STATE $u_k$ s nejnižší hodnotou $\bar{J}$ (ozn. $\tilde{u}_{0,j}$) uchovej jako řízení pro $\tilde{H}_{0,j}$.
\STATE $\pi_i:=\{\tilde{u}_{0,1:j-1},\tilde{u}_{0,j},u_{0,k+1:|\tilde{H}_0|}\}$
\ENDFOR
\STATE $\pi_{i+1}:=\{\tilde{u}_{0,1:\tilde{H}_0|}\}$
\ENDFOR
\end{algorithmic}
\caption{Schéma algoritmu SIDPS}
\label{alg}
\end{algorithm}

\section{Integrátor s neznámým ziskem}
Výstup systému je popsán jako
\begin{gather}
\label{simple}
y_{t+1}=y_t+\theta u_t+v_{t+1},	\qquad t=0,\ldots,N-1,\\
v_{t+1}\sim N(0,\sigma^2),
\end{gather}
kde $\theta\neq0$ je neznámý parametr a rozptyl šumu $\sigma^2$. Počáteční hodnota výstupu je nastavena na $y_0=1$ a optimální řízení je takové, pro které $y_t=0$ . Systém je možné v případě potřeby vhodně posunou či přeškálovat, není tedy potřeba uvažovat jiné hodnoty počátečního výstupu a referenčního signálu .

O neznámém parametru $\theta$ máme v čase $t$ informaci v podobě dostatečné statistiky $T_t=(\hat{\theta}_t,P_t)$ tvořené střední hodnotou a rozptylem. Předpokládáme nekorelovanost odhadu $\theta_t$ s šumem, tedy že
\begin{equation}
\cov(v_{t+1},\theta_t)=0.
\end{equation}

Ztrátová funkce je kvadratická v $y_{t+1}$, čili
\begin{equation}
\label{glos}
g(y_{1:N},u_{0:N-1})=\sum_{t=0}^{N-1}y_{t+1}^2.
\end{equation}

Odhadovací procedurou pro parametr $\theta$ je Kalmanův filtr. Pro systém \eqref{simple}  má tvar
\begin{gather}
K_t=\frac{u_tP_t}{u_t^2P_t+\sigma^2},\\
\label{kal}
\hat{\theta}_{t+1}=\hat{\theta}_t+K_t(y_{t+1}-u_t\hat{\theta}_t),\\
P_{t+1}=(1-K_tu_t)P_t. 
\end{gather}

Hyperstav systému $H_t$ tvoří vektor $(y_t,\hat{\theta}_t,P_t)$. Očekávaná ztráta je
\begin{equation}
J_t(H_t)=\min_{u_t \in U_t}\E_{y_{t+1}}\left\{y_{t+1}^2+J_{t+1}(H_{t+1})|H_t,u_t\right\}, \qquad t=0,\ldots,N-1.
\end{equation}

Ta po dosazení z \eqref{simple} a částečném provedení střední hodnoty přejde na tvar
\begin{gather}
\label{dos}
J_t(H_t)=\min_{u_t \in U_t}\left\{(y_t+\hat{\theta}_tu_t)^2+u_t^2P_t+\sigma^2+\E_{y_{t+1}}(J_{t+1}(H_{t+1}))|y_t,\hat{\theta}_t, P_t,u_t\right\}.
\end{gather}

\subsection{Transformace rovnic systému}
Od proměnných $(y_t,\hat{\theta}_t,P_t,u_t)$ přejdeme k  $(\eta_t,\beta_t,\zeta_t,\nu_t)$ dle vztahů
\begin{gather}
\eta_t=\frac{y_t}{\sigma}, \\
\label{trab}
\beta_t=\frac{\hat{\theta}_t}{\sqrt{P_t}}, \\
\zeta_t=\frac{1}{\sqrt{P_t}}, \\
\nu_t=\frac{u_t\sqrt{P_t}}{\sigma}.
\end{gather}

Současně můžeme neurčitost ve výstupu \eqref{simple} reprezentovat jedinou normalizovanou náhodnou veličinou
\begin{equation}
s_t=\frac{y_{t+1}-y_t+\hat{\theta}_tu_t}{\sqrt{u_t^2P_t+\sigma^2}} \sim N(0,1).
\end{equation}
 
Rovnice pro výstup \eqref{simple} a následující odhad neznámého parametru \eqref{kal} tak přejde v
\begin{gather}
\label{eta}
\eta_{t+1}=\eta_t+\beta_t\nu_t+\sqrt{1+\nu_t^2}s_t,\\
\label{beta}
\beta_{t+1}=\sqrt{1+\nu_t^2}\beta_t+\nu_ts_t.
\end{gather}

Přejdeme-li k vhodně upravené očekávané ztrátě, dostaneme 
\begin{align}
\label{bel}
V_t(\eta_t,\beta_t,\zeta_t)&=\frac{J_t(y_t,\hat{\theta}_t,P_t)}{\sigma^2}\\
&=\min_{\nu_t }\left\{(\eta_t+\beta_t\nu_t)^2+\nu_t^2+1+\E_{s_t}(V_{t+1}(\eta_{t+1},\beta_{t+1},\zeta_{t+1}))\right\}.
\end{align}

Nyní spočteme očekávanou ztrátu pro $N-1$ jako
\begin{equation}
V_{N-1}(\eta_{N-1},\beta_{N-1},\zeta_{N-1})=\min_{\nu_{N-1}}\left\{(\eta_{N-1}+\beta_{N-1}\nu_{N-1})^2+\nu_{N-1}^2+1\right\}.
\end{equation}

Pomocí diferenciálního počtu pak získáme optimální zásah ve tvaru
\begin{equation}
\label{optcon}
\nu_{N-1}=-\frac{\eta_{N-1}\beta_{N-1}}{1+\beta_{N-1}^2}
\end{equation}
a očekávanou ztrátu rovnu
\begin{equation}
V_{N-1}(\eta_{N-1},\beta_{N-1},\zeta_{N-1})= \frac{\eta_{N-1}^2+1}{\beta_{N-1}^2+1}.
\end{equation}

Protože optimální zásah $\nu_{N-1}$ ani očekávaná ztráta $V_{N-1}$ nezávisí na $\zeta_{N-1}$, díky tvaru $V_t$ nebude rovněž optimální zásah $\nu_t$ a očekávaná ztráta $V_t$ záviset na $\zeta_t$. K nalezení optimálního řízení tedy stačí v každém čase $t$ uvažovat pouze dvourozměrný hyperstav $H_t=(\eta_t,\beta_t)$. Navíc můžeme bez újmy na obecnosti počítat optimální zásah pouze pro kladné hodnoty $\eta_t$ a $\beta_t$. Optimální řízení tedy napočteme pouze pro kladné hodnoty $y_t$ a $\theta_t$. Pro ostatní možnosti pak díky tvaru \eqref{simple} dostaneme požadovaný řídící zásah vhodnou volbou znaménka napočteného zásahu. 

\subsection{Opatrné řízení}
Metoda opatrného řízení 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.

Optimální řídící zásah je pro metodu opatrného řízení dán dle \eqref{optcon} jako 
\begin{equation}
\label{optcon2}
\nu_t=-\frac{\eta_t\beta_t}{1+\beta_t^2}.
\end{equation}

Pro malé hodnoty $\beta_t$ (tedy pro velké hodnoty neurčitosti v identifikaci parametru $\theta$, viz \eqref{trab}) bude řízení velmi konzervativní (odtud název). Pro $\beta_t=0$ (to odpovídá $\hat{\theta}_t=0$) pak bude $\nu_t=0$ a regulátor tedy nebude poskytovat žádné aktivní řízení. Opatrné řízení je vhodným kandidátem na $u^{(1)}_t$ v \ref{rozklad}. 

\subsection{Klasický přístup k dynamickému programování}
Při použití klasického numerického přístupu k dynamickému programování je prostor hyperstavů diskretizován do mřížky (zde 64x64). Pro každý bod hyperstavu se napočítá očekávaná ztráta, mimo body mřížky se použije interpolace (zde kubická). K numerické integraci se použije klasická Simpsonova metoda a k nalezení minima pak jednoduchá interpolační metoda (zde se každými třemi body na mřížce proložila parabola, nalezlo její minimum a to se pak testovalo, zdali je globálním minimem očekávané ztráty).

Optimální řízení na výsledné mřížce bylo nakonec parametrizováno analytickou formulí tvaru
\begin{equation}
\label{ast}
\nu_t=-\frac{0,56+\beta_t}{2,2+0,08\beta_t+\beta_t^2}\eta_t-\frac{1,9}{1,7+\beta_t^4}.
\end{equation}
Analytická aproximace způsobí zvýšení ve ztrátové funkci o méně než 1 \%.

První člen v \eqref{ast} můžeme interpretovat jako modifikované opatrné řízení, druhý člen pak jako budící složku řízení. 

\subsection{Implementace SIDPS pro jednouchý systém}
\label{implSIDPS}
Budící složka v \eqref{ast} závisí pouze na $\beta$. I to nás vede k tomu, že řízení budeme hledat ve tvaru
\begin{equation}
\nu_t=-\frac{\eta_t\beta_t}{1+\beta_t^2} -\nu^{(2)}_t(\beta_t),
\end{equation}
Volíme tedy $\tilde{H}_t=\beta_t$ a optimalizujeme $\nu^{(2)}$.

Dle popisu (a následné transformace) systému \eqref{simple} je pro výpočet optimálního $\nu^{(2)}_t$ třeba diskretizovat (kladnou) část $1$-dimenzionálního prostoru proměnné $\tilde{H}_0$.

Implementace všech částí algoritmu byla provedena v souladu s oddílem \ref{sidp'}. Konkrétní nastavení parametrů algoritmu zachycuje tabulka \ref{par2}. Mimo body diskretizace byla použita lineární interpolace. Pro porovnání kandidátů na řízení byl navíc implementován RSSS algoritmus, což vedlo k drobnému zlepšení výsledného řízení (snížení ztráty o cca 1\% a stabilnější kovergence $\pi_i$).

\begin{table}
\begin{tabular}{|l|ll|} 
\hline
délka ubíhajícího horizontu & $\tilde{N}$ & 10\\
počet iterací algoritmu & $n_{iter}$ & 4\\
počet bodů v diskretizaci každé dimenze $H_t$ & $n_g$ & 10\\
počet kadnidátů na změnu  řídícího zásahu & $m$ & 7\\
počáteční rozsah pro hledání optimálního řídícího zásahu & $\beta_0$ & 5\\
parametr pro redukci $\beta_0$ při opakování algoritmu &$\gamma$ & 0,2\\
počet realizací pro odhad metodou Monte Carlo & $n$ & 20\\\hline
\end{tabular}
\caption{Konkrétní volba parametrů algoritmu SIDPS}
\label{par2}
\end{table}

Otázkou zůstává, jak volit $\tilde{H}_0^{(c)}=\eta_0$. Ačkoliv předpokládáme, že $\nu^{(2)}_t$ na $\tilde{H}_0^{(c)}$ nezávisí, skutečnost může být jiná. V tom případě je zřejmě vhodné volit $\tilde{H}_0^{(c)}$ tak, aby co nejvíce odpovídala stavu, ve kterém se bude systém skutečně nacházet. Jako vhodný se proto jeví návrh $\eta_0=0$. Opodstatnění volby spočívá v tom, že v případě úspěšného řízení bude systém v okolí optimální (tedy nulové) hodnoty.

Optimální zásah $\nu^{(2)}_0$ byl počítán pro hodnoty
\begin{equation}
\tilde{H}_0=\beta_0=\{0.0001+ 0.75\times j|j=0,\ldots,9\}.
\end{equation}
Výsledné řízení se pak dalo s úspěchem použít pro širokou škálu hodnot $(b_0,P_0,y_0,\sigma)$ bez nutnosti opakovat vlastní výpočet. Veškeré experimenty tak byly prováděny s řízením získaným při tomto jediném výpočtu. Poznamenejme, že v případě klasického SIDP by to z důvodu lokální povahy algoritmu nebylo možné. 

Výpočet za daných parametrů trval v řádech jednotek vteřin (cca 5). Charakter konvergence $\pi_i$ je diskutován dále.

\section{Srovnání jednotlivých přístupů}
\subsection{Kvantitativní srovnání}
\label{kvansrov}
Řízení bylo testováno vzhledem k různým hodnotám počátečního odhadu $\theta$, tedy k různým hodnotám $(\hat{\theta}_0,P_0)$. Řídící horizont byl vždy $N=100$, rozptyl šumu $\sigma=1$. Každá simulace byla opakována 1000x, uvedené hodnoty celkové ztráty jsou průměrem přes jednotlivé realizace. 

Pro každé jednotlivé opakování simulace byla skutečná hodnota parametru $\theta$ pro počáteční hodnoty $(\hat{\theta}_0,P_0)$ zvolena jako první nenulová realizace náhodné veličiny s rozdělením $N(\hat{\theta}_0,P_0)$. Tedy apriorní informace odpovídá skutečnosti. Pro takto vygenerovanou hodnotu $\theta$ byly postupně aplikovány jednotlivé řídící algoritmy. Pro snížení vlivu náhodnosti při porovnání kvality řízení byly všechny realizace šumu (v rámci jednoho opakování simulace) v průběhu řízení jednotlivými algoritmy voleny stejně.

Na obrázku \ref{kvan} jsou zachyceny výsledky pro různé hodnoty apriorní informace $(\hat{\theta}_0,P_0)$. Hodnoty ztráty jsou uvedeny relativně k průměrné ztrátě, které bylo ve stejné situaci dosaženo pomocí DP regulátoru. Dle obdržených výsledků tak v průměru nový algoritmus úspěšně konkuruje DP řízení, jehož výpočet je mnohonásobně složitější (dle informací v článku o SIDP je časová složitost cca $100000\times$ vyšší).

\begin{figure}
\centering
\includegraphics[width=\textwidth]{kvan}
\caption{Kvantitativní porovnání regulátorů pro různé hodnoty $(\hat{\theta}_0,P_0)$ - Na ose $y$ jsou relativní hodnoty průměrné ztráty, které dosáhly uvedené algoritmy (osa $x$). Referenční hodnotou byly průměrné ztráty, kterých dosáhl DP regulátor pro tytéž parametry.}
\label{kvan}
\end{figure}

Pro ilustraci závislosti výsledků na rozptylu šumu $\sigma$ a počátečním stavu $y_0$ byl zopakován výše uvedený pokus pro $\sigma=10$ a poté ještě pro $y_0=100$. Bylo použito stejné řízení, jako prve. Výsledky byly velmi podobné těm na obrázku \ref{kvan}. Celkově vykazovalo řízení pomocí popsaného algoritmu v průměru srovnatelně kvalitní řízení jako DP regulátor a to pro všechny testované $(\hat{\theta}_0,P_0,y_0,\sigma)$ se zanedbatelnými výpočetními nároky.

\subsection{Kvalitativní srovnání}
Při kvalitativním srovnání byly zkoumány četnosti konkrétních realizací ztráty. Obrázky \ref{kval0}, \ref{kval1}, \ref{kval10} udávají četnosti realizací relativních ztrát, kterých bylo nabyto při použití konkrétního regulátoru pro různé hodnoty $\hat{\theta}_0$ ($P_0=10$ a 1000 opakování). Referenční hodnotou byla průměrné ztráta DP regulátoru. Ostatní nastavení bylo stejné jako v případě kvantitativního srovnání. Uvedené histrgramy dokumentují, že kvalitativně je řízení získané novým algoritmem velmi odobné jako řízení DP regulátorem.
 
\begin{figure}
\centering
\includegraphics[width=\textwidth]{kval0}
\caption{Kvalitativní porovnání regulátorů pro $\hat{\theta}_0=0$ - Na ose $x$ jsou hodnoty relativní ztráty (referenční hodnotou je průměrná ztráta DP regulátoru), na osa $y$ pak četnosti relativní ztráty z 1000 opakování. Zvýrazněna je hodnota průměrné ztráty (červená) a mediánu (zelená).}
\label{kval0}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{kval1}
\caption{Kvalitativní porovnání regulátorů pro $\hat{\theta}_0=1$ - Na ose $x$ jsou hodnoty relativní ztráty (referenční hodnotou je průměrná ztráta DP regulátoru), na osa $y$ pak četnosti relativní ztráty z 1000 opakování.Zvýrazněna je hodnota průměrné ztráty (červená) a mediánu (zelená).}
\label{kval1}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{kval10}
\caption{Kvalitativní porovnání regulátorů pro $\hat{\theta}_0=0$ - Na ose $x$ jsou hodnoty relativní ztráty (referenční hodnotou je průměrná ztráta DP regulátoru), na osa $y$ pak četnosti relativní ztráty z 1000 opakování.Zvýrazněny jsou pak hodnoty průměrné ztráty (červená) a mediánu (zelená).}
\label{kval10}
\end{figure}

\section{Charakter konvergence SIDPS}
Obrázek \ref{konv} dává představu o charakteru konvergence $\pi_i$. První část obrázku zachycuje jednotlivé iterace $\pi_i$ (hodnoty $u^{(2)}_0$ v závislosti na $\beta$). Druhá pak relativní průměrnou ztrátu při řízení za použití $\pi_i$ (průměr přes 1000 opakování). Při simulacích bylo $(\hat{\theta}_0,P_0)=(1,10)$, ostatní parametry zůstaly stejné jako v předešlé sekci. Referenční hodnotou byla opět průměrná ztráta DP regulátoru.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{konv}
\caption{Charakter konvergence SIDPS - Průběh $\pi_i$ (první část) a relativní průměrná ztráta, které dosahuje (druhá část). Pro výpočet průměrné ztráty při řízení bylo $(\hat{\theta}_0,P_0)=(1,10)$ a jako referenční hodnota byla zvolena průměrná ztráta DP regulátoru. Ostatní parametry simulace viz. sekce \ref{kvansrov}}
\label{konv}
\end{figure}

\section{Možné volby $\tilde{H}_0^{(c)}$}
V předchozí srovnávací části jsme volili konstantně $\tilde{H}_0^{(c)}=\eta_0=0$. Pokud se systém nachází poblíž této hodnoty ($\eta_t \approx 0$), řízení by mělo být (a dle pokusů také je) dobré. Pokud ale systém bude v jiném stavu (to je možné zejména v počáteční fázi při volbě $y_0 \gg 0$), je zde nebezpečí, že se regulátor dopustí zbytečně velké chyby. Ilustrací tohoto jevu se zabývá následující sekce.

Při nastavení ze sekce \ref{implSIDPS} bylo algoritmem SIDPS napočítáno optimální řízení pro následující různé volby $\tilde{H}_0^{(c)}=\eta_0$
\begin{enumerate}
\item[0)]  $\eta_0=0$,
\item[1)]  $\eta_0=\eta_{max}*\text{\emph{rand}}$,
\label{2}
\item[2)]  $\eta_0=\frac{\eta_{max}}{10^{\beta_0}}$,
\label{3}
\end{enumerate}
kde $\eta_{max}$ je maximimální uvažovaná hodnota $\eta_0$, pro kterou budeme chtít řízení použít a \emph{rand} je náhodné číslo z intervalu $<0,1>$. Algoritmem SIDPS jsme tak získali tři (různé) návrhy na řízení. Volba 0) odpovídá původnímu nastavení (tj. předpokladu, že systém se bude nacházet poblíž optimální hodnoty). Předpis 1) počítá s libovolnou hodnotu $\eta_0$ z předpokládaného rozsahu a 2) vyjadřuje předpoklad souvislosti mezi počátešním stavem $\eta_0$ a neurčitostí $\beta_0$. Předpoklad se dá interpretovat tak, že čím nižší neurčitost je, tím více o systému víme a tedy jsme (pravděpodobně) blíže optimální hodnotě. 

Kvantitativní rozdíl řízení je diskutován. Parametry pokusu jsou převzaty ze sekce \ref{kvansrov} ($\sigma=1$, čili $\eta_t=y_t$). Výsledné řízení bylo aplikováno pro různé hodnoty $y_0 \in <0,100>$ (tedy $\eta_{max}=100$), pro $P_0=10$ a $\hat{\theta}_0 \in \{0,1,10\}$. Pro $P_0=1$ vyšly výsledky obdobně. Zajimavé bylo srovnání zejména pro $\hat{\theta}_0 =0$, v ostatních variantách bylo řízení vždy velmi dobré (průměrná relativní ztráta $\doteq 1$).

Obdržené výsledky ilustruje obrázek \ref{Hc}. 

\begin{figure}
\centering
\includegraphics[width=\textwidth]{Hc}
\caption{Závislost průměrné relativní ztráty při různých volbách $\tilde{H}_0^{(c)}$ (viz. text) -  Na ose $x$ jsou hodnoty skutečné počáteční výchylky systému $y_0$, na ose $y$ pak hodnoty průměrné relativní ztráty z 1000 opakování. Referenční hodnotou je průměrná ztráta DP regulátoru}
\label{Hc}
\end{figure}

Dle očekávání je pro $y_0\approx0$ nejvýhodnější použít volbu 0) (předpoklad $\eta_t=0$, při kterém  bylo řízení navrženo, platí s dobrou přesností vždy). Pro vyšší hodnoty $y_0$ ($>20$) je předpokladané $\eta_t=0$ pro $t=0$ velmi špatný odhad a je dosaženo vyšší ztráty.

V případě volby 1) byla situace opačná. Řízení bylo úspěšné zejména pro vyšší hodnoty $y_0$. Pro $y_0\approx0$ je budící složka zbytečně veliká (zhruba 2x oproti 0)) a je dosaženo vyšší ztráty.

Kompromisnisním řešením volby $\tilde{H}_0^{(c)}=\eta_0$ je možnost 2). Ta poskytuje v průměru dobré řízení v celém intervalu uvažovaných $y_0$.

V případě $y_0 \in <0,1>$ (tedy $\eta_{max}=1$), poskytly všechny volby přibližně stejné řízení (rozdíl v průměrné relativní ztrátě, $<1\%$, je způsobený náhodnými realizacemi poruch). Hodnota průměrné relativní ztráty pak byla pro všechny hodnoty $y_0 \in <0,1>$ $\doteq1$.

\section{Závěr}
Na základě provedených pokusů se řízení získané pomocí SIDPS jeví jako velmi dobré. Dává srovnatelné výsledky s DP regulátorem, navíc v zanedbatelném čase. Nevýhodou algoritmu SIDPS je nutnost vybrat veličiny, které budou tvořit $\tilde{H}_0$, předepsat $u^{(1)}_t$ a zvolit $\tilde{H}^{(c)}_0$. Dimenze $\tilde{H}_0$ má vliv zejména na náročnost výpočtu, $u^{(1)}_t$ a $\tilde{H}^{(c)}_0$ pak ovlivňují kvalitu výsledného řízení. Volba výše uvedených komponent algoritmu by měla vycházet z rozboru konkrétního systému, pro který má být algoritmus SIDPS implementován.
\end{document}
