root/applications/doprava/texty/novotny_vyzk_LQ/Implementation/Minimalization.tex @ 1430

Revision 1430, 7.8 kB (checked in by jabu, 12 years ago)
Line 
1\subsection{Minimalizace kritéria}
2Podobně jako v předchozí kapitole budeme minimalizovat kvaratické kritérium $J$.
3Minimalizaci budeme provádět v proměnných $u(t)$ pro časový horizont $h$,
4tedy pro vektory $\Delta C(t_0), \Delta C(t_0 + 1), ...,u(t_0 + h)$. Pro zpřehlednění
5zápisu položíme bez újmy na obecnosti $t_0 = 0$.
6Matici $I_0(t)$, která vyjadřuje příjezd vozidel z okolí do sledované sítě,
7budeme v rámci minimalizačního horizontu považovat za konstantu značenou $I_0$ a
8vztah \ref{eq:my_trans_mat} se dá tedy přepsat do tvaru
9\begin{equation}\label{eq:prechod_subs_01}
10  \left( \begin{array}{c} q(t+1) \\ 1 \end{array} \right)  =
11    \left( \begin{array}{cc} A_0 & I_0 \\ 0 & 1 \end{array} \right)
12    \left( \begin{array}{c} q(t) \\ 1 \end{array} \right)
13    +
14    \left( \begin{array}{c} B_0 \\ 0  \end{array} \right)
15    u(t)
16\;.
17\end{equation}
18Po provedení substituce
19\begin{equation}
20  \begin{array}{cccc}
21    x(t) = \left( \begin{array}{c} q(t) \\ 1 \end{array} \right) \;, &
22    A = \left( \begin{array}{cc} A_0 & I_0 \\ 0 & 1 \end{array} \right) \;, &   
23    B = \left( \begin{array}{c} B_0 \\ 0  \end{array} \right) \;,
24  \end{array} 
25\end{equation}
26se rovnice \ref{eq:prechod_subs_01} zjednoduší na
27\begin{equation}\label{eq:prechod_mat_po_subs}
28 x(t+1) = Ax(t) + Bu(t) \;,
29\end{equation}
30kde chybí konstantní člen. Poté již můžeme použít metodu popsanou v kapitole \ref{sec:minim}.
31
32% Kvadratické kritérium ve tvaru
33% \begin{equation}\label{eq:J}
34%  J = \sum_{t=0}^h x(t)^T Q x(t) + u(t)^T R u(t) \;, 
35% \end{equation}
36% kde $Q$ a $R$ jsou diagonální pozitivně definitní matice vah,
37% rozepíšeme pomocí vztahu \ref{eq:prechod_mat_po_subs}
38% a budeme minimalizovat pro jednotlivá $t$ podle $u(t)$.
39% Proměnná $u(h)$ se v kritériu vyskytuje pouze ve členu
40% \begin{equation}
41% u(h)^T Ru(h)
42% \end{equation}
43% a tedy musí být $u(h) = 0$. $u(h - 1)$ se vyskytuje ve členech
44% \begin{equation}
45%  u(h - 1)^T Ru(h - 1) + (x(h - 1)^T A^T +u(h - 1)^T B^T ) Q ( A x(h - 1) + B u(h - 1) + I_0) \;,
46% \end{equation}
47% což lze pomocí složených vektorů a matic přepsat do tvaru
48% \begin{equation} \label{eq:J_sloz}
49%  (u(h-1)^T, x(h - 1)^T )
50%   \left( \begin{array}{cc}
51%           B^T \sqrt{Q}^T & \sqrt{R}^T \\
52%         A^T \sqrt{Q}^T & 0
53%         \end{array}
54%   \right) 
55%   \left( \begin{array}{cc}
56%            \sqrt{Q} B & \sqrt{Q} A \\
57%         \sqrt{R} & 0
58%         \end{array}
59%   \right)
60%   \left( \begin{array}{c}
61%          u(h-1) \\
62%           x(h - 1)
63%         \end{array}
64%   \right) \;,
65% \end{equation}
66% kde symbolem $\sqrt{Q}$ myslíme matici, pro níž platí $\sqrt{Q}^T \sqrt{Q} = Q$.
67% Matice $Q$ a $R$ jsou pozitivně definitní diagonální matice, takže $\sqrt{Q}$ bude také
68% pozitivně definitní diagonální a její diagonální prvky budou rovny odmocnině
69% příslušných prvků původní matice $Q$.
70% Složenou matici
71% \begin{equation} 
72%   M_0 =  \left( \begin{array}{cc}
73%            \sqrt{Q} B & \sqrt{Q} A \\
74%         \sqrt{R} & 0
75%         \end{array}
76%   \right)
77% \end{equation}
78% můžeme pomocí QR dekompozice rozložit na součin
79%
80% \begin{equation}
81%  M_0 = M_R M_Q ,
82% \end{equation}
83% kde $M_Q$ je horní trojůhelníková matice a $M_R$ je matice ortonormální.
84% Pro ortonormální matici $M_R$ platí
85% \begin{equation}
86%  M_R^T M_R = I \;,
87% \end{equation}
88% neboť její zloupce jsou vzájemně ortogonální a normované na jednotku.
89% Z toho plyne, že součin matic $M_0^T M_0$ vyskytující se v členu \ref{eq:J_sloz}
90% můžeme pomocí QR rozlkladu převést na tvar
91% \begin{equation}
92%  M_0^T M_0 = ( M_R M_Q )^T  M_R M_Q = M_Q^T M_R^T M_R M_Q = M_Q^T M_Q \;,
93% \end{equation}
94% kde
95% \begin{equation}
96%  M_Q =  \left( \begin{array}{cc}
97%           L_u & L \\
98%         0 & L_q
99%         \end{array}
100%   \right) \;,
101% \end{equation}
102% je horní trojůhelníková matice.
103% Člen \ref{eq:J_sloz} poté přejde na tvar
104% \begin{equation}
105%  (\Delta C(h-1)^T, x(h - 1)^T )
106%   \left( \begin{array}{cc}
107%           L_q^T & 0 \\
108%         L^T & L_u^T
109%         \end{array}
110%   \right) 
111%   \left( \begin{array}{cc}
112%           L_u & L \\
113%         0 & L_q
114%         \end{array}
115%   \right)
116%   \left( \begin{array}{c}
117%          u(h-1) \\
118%           x(h - 1)
119%         \end{array}
120%   \right) \;,
121% \end{equation}
122% který můžeme dále upravit na
123% \begin{equation}
124%  ( L_u u(h-1) +  L_q x(h-1) )^T ( L_uu(h-1) +  L_q x(h-1) ) + x(h-1)^T L_q^T L_q x(h-1)
125% \end{equation}. Pokud tento člen chceme minimalizovat v proměnné $x(h-1)$, musí nutně platit
126% \begin{equation}
127% u(h-1) = - L_u^{-1} L x(h-1)
128% \end{equation}. Zbývající nenulová část $x(h-1)^T L_q^T L_q x(h-1) $ lze přepsat pomocí rovnice \ref{eq:prechod} do tvaru
129% \begin{equation}
130%  x(h-1)^T L_q^T L_q x(h-1) = ( A x(h-2) + B u(h-2) )^T L_q^T L_q ( A x(h-2) + B u(h-2) ) \;,
131% \end{equation}
132% kde se vyskytuje $u(t)$ pouze pro $t = h-2$.
133% Při minimalizaci podle $u(h-2)$ musíme zahrnout i tuto část, člen ve složeném maticovém zápisu
134% tudíž nabyde tvaru
135% \begin{equation}
136%  (u(h-2)^T, x(h-2)^T )
137%   \left( \begin{array}{ccc}
138%           B^T \sqrt{Q}^T & \sqrt{R}^T & L_x A \\
139%         A^T \sqrt{Q}^T & 0 & L_x B
140%         \end{array}
141%   \right) 
142%   \left( \begin{array}{cc}
143%            \sqrt{Q} B & \sqrt{Q} A \\
144%         \sqrt{R} & 0 \\
145%         L_x A & L_x B
146%         \end{array}
147%   \right)
148%   \left( \begin{array}{c}
149%          u(h-2) \\
150%           x(h-2)
151%         \end{array}
152%   \right) \;.
153% \end{equation}
154% Matici
155% \begin{equation}
156%  M = \left( \begin{array}{cc}
157%            \sqrt{Q} B & \sqrt{Q} A \\
158%         \sqrt{R} & 0 \\
159%         L_x A & L_x B
160%         \end{array}
161%   \right) = \left( \begin{array}{c}
162%            M_0  \\
163%         L_x A \; L_x B
164%         \end{array}
165%   \right)
166% \end{equation}
167% opět převedeme QR dekompozicí do horního trojůhelníkového tvaru a dostaneme vztah pro $u(h-2)$.
168% Tento postup se analogicky aplikuje na každé $u(t)$. Matice $M_0$ je stále stjná, $L_x$ použijeme z předchozího kroku minimalizace.
169
170\subsubsection{Implementace minimalizace}sec:minim
171
172Minimaizace kvadratického kritéria je v simulaci prováděna metoduo \texttt{mat L( const int horizont )} ve třídě
173\texttt{ QuadraticMinimalizator }, které se jako parametry konstruktoru předávají matice \texttt{A}, \texttt{B}, \texttt{L} a \texttt{Q}.
174Pro numerické výpočty je použite knihovna IT++, která umožňuje jednoduché operace s maticemi a vektory syntakticky podobně jako v jazyku
175MATLAB. Jsou zde také implementovány jednoduché algoritmy, jako například pro naše účely potřebný QR rozklad.
176Níže je vložen zdrojový kód metody pro minimalizaci kritéria a získání hodnot řídících parametrů. Proměnné jsou v kódu značeny
177stejně jako v textu výše. Ty, které jsou v kódu použity navíc, jsou popsány v komentářích.
178
179\begin{lstlisting}
180 mat L( const int horizont ) {
181        int xdim = A.cols(); // Dimenze vektoru x
182        int udim = B.cols(); // Dimenze vektoru u
183        mat sqQ = sqrt(Q);   // Odmocninove matice
184        mat sqR = sqrt(R);
185        mat M0, M;
186        mat qrQ;
187        mat L, Ls, Lx, Lu;
188        // sestaveni matice M0
189        M0 = concat_vertical(                           
190                concat_horizontal( sqQ * B,     sqQ * A ),
191                concat_horizontal( sqR,         zeros( udim, xdim ) )
192            );
193        M = M0;
194        // hlavni cyklus
195        for (int h = horizont; h >= 0; h --) {
196            qr(M, L);  // QR rozklad matice M
197            // rozklad matice L na do slozek
198            Lu = L(0,    udim-1,       0,    udim-1);     
199            Lx = L(udim, udim+xdim-1,  udim, udim+xdim-1);
200            Ls = L(0,    udim-1,       udim, udim+xdim-1);
201            // kompozice matice M
202            M = concat_vertical(
203                    M0,
204                    concat_horizontal(Lx * B, Lx * A )
205                );
206        }
207        // navratova hodnota - vysledna matice L : u(t+1) = Lx(t)
208        return - inv(Lu) * Ls;
209    }
210\end{lstlisting}
Note: See TracBrowser for help on using the browser.