[數學分析] 連續函數族的逐點上包絡函數不一定連續

連續函數有諸多用途，一般在參數最佳化領域中常見的情況是考慮所謂的上包絡函數(upper envelope function)。

Definition: 定義函數族 \(\{f_t : t \in T\} \) 其中 \(T\) 為 index set 並考慮對任意 \(x \in X\)，現在定義上包絡函數(upper envelope function) 或者 逐點上確界函數(pointwise supremum function)
$$ F(x) := \sup_t f_t(x)$$

Question: 一個有趣的問題是如果這些函數族成員都是連續函數，那麼取 supremum 之後所得到的新函數 \(F\) 是否仍為連續呢？

答案是否定的。以下例子說明甚至是定義在緊緻集合上的連續函數族也沒有保證上包絡函數連續。

Example: 考慮一連續函數族 \( \{f_t: t \in [0, T]\} \) 其中 \(f_t(x) := x^t\) 對 \(x \in [0,1]\) 且 \(t \in [0,1]\) 並定義 \(f_0 = 0\)。 則函數族的上包絡函數為 $$ \sup_t f_t(x) = \begin{cases} 1, & x \in (0, 1] \\ 0, & x = 0\end{cases}$$ 讀者不難發現此函數在 \(x=0\) 處有不連續跳點。

Comment: 在最佳控制與數理經濟中有個非常有用的定理可以刻畫上包絡函數的連續性稱作 Berge's Maximum Theorem 有興趣的讀者可以自行查閱。

[控制理論] 非線性模型預測控制 (0) - 引論

此文我們針對 非線性模型預測控制 (Nonlinear Model Predictive Control, NMPC) 做簡單的介紹： 是一套針對受控廠為 非線性系統 所發展的 回授控制理論，其控制律透過不斷求解最佳化問題而得。以下我們將其簡稱 NMPC。


基本想法：
假設你是一個圍棋高手正在與旗鼓相當的對手對弈，那麼你可能在心中盤算好幾步可能的走法 並同時 試圖來 "預測" 對手的路數，但輪到自己下子的時候，我們只能選取剛才在心中盤算出的所有走法中 " 最佳" 的那一步，並且只動那一步旗，動了之後都必須重新盤算上述過程。這個精神大概與 模型預測控制 本質幾乎相同。就是逐步最佳化，慢慢朝向目標前進。以下我們會逐步將此概念規範化，讀者可以讀完之後再回頭瞧瞧這個基本想法，也許會發現異曲同工之妙。


對 $n=0,1,2,...$，令 $x(n)$ 為當前系統狀態 且 $x^{ref}(n)$ 為 給定的參考目標軌跡。

模型預測控制 的 主要目的：
與一般控制問題相仿，模型預測控制的主要目的不外乎以下兩大類問題：
1. 鎮定(stabilizing)問題。使預測的狀態軌跡 盡可能趨近 零點。
2. 追蹤(tracking)問題。使 預測的狀態軌跡 盡可能跟隨給定的 參考目標軌跡 。

Comments:
1. 若稍有控制背景的讀者不難發現鎮定問題其實是追蹤問題的特例，更進一步地說，所謂追蹤控制問題是指：決定一組控制輸入 $u(n)$ 使得 $x(n)$ 能盡可能緊跟給定的參考狀態軌跡 $x^{ref}(n)$。換而言之，若當前狀態 $x(n)$ 與 $x^{ref}(n)$ 所差甚遠，則我們將盡可能控制此系統軌跡趨近 $x^{ref}(n)$：若當前狀態 $x(n) = x^{ref}(n)$ 則我們將盡可能控制此系統"維 持"在該狀態。

2. 關於鎮定問題提及的對零點追蹤 或者 追蹤問題的參考目標軌跡，都必須滿足隱藏的前提： 零點 或者 參考目標軌跡 必須是系統的 平衡點 (equilibrium point)。否則該系統無法進行鎮定 or 追蹤。


以下我們將介紹一類簡化的 NMPC 問題，並藉此展示此控制方法的一些基本精神。

一類簡化的非線性模型預測控制問題：
令 $x(n) \in X := \mathbb{R}^d$ 且 $u(n) \in U := \mathbb{R}^m$ 考慮 參考軌跡狀態(reference state trajectory)，記作 $x^{ref}(\cdot)$，寫作
\[
x^{ref}(n) :=x_* :=0,\;\;\;\;\; n\geq 0
\]上述 對零點狀態追蹤問題 退化為 鎮定問題。

以控制理論的基本想法，我們很自然想要對 當前狀態 $x(n)$ 與 目標 $x^{ref} = 0$ 之誤差 做回授控制，為了達成此目標，我們期望控制力 $u$ 必須具備以下形式：
$$
u(n) := \mu (x(n))
$$其中 $\mu $ 為 將狀態 $x \in X$ 映到控制輸出值之集合 $U$ 之映射。

系統模型動態(System Model Dynamics)：
模型預測控制的想法在於利用 系統的模型來進行未來狀態軌跡之預測，並且給出適當的目標函數後對其做最佳化求解控制力。更精確的說，系統動態 一般我們寫成下列方程
\[
x^+ = f(x,u) \;\;\;\;\; (*)
\]其中 $f: X \times U \to X$ 為已知且非線性 滿足 將狀態 $x$ 與控制輸出值 $u$ 映到下一個時刻的狀態 $x^+$ 。

預測模型動態(Predictive Model Dynamics)：
現在給定當前狀態 $x(n)$，對任意給定控制序列
\[
u(0),u(1),...,u(N-1)
\]滿足控制區間 $N \geq 2$，我們可以用此 控制序列 以及 給定的當前狀態 作為初始狀態來迭代 $(*)$，如此可得 預測狀態軌跡 $x_u$ 滿足下式：
\[
x_u(0) = x(n); \;\;\;\; x_u(k+1) = f(x_u(k), u(k)),\;\;\;\; k=0,...,N-1 \;\;\;\; (**)
\]由於 $x(n)$ 被作為初始狀態來進行預測，故我們 $x_u(k)$ 實際上為狀態 $x(n+k)$ (對時刻 $t_{n+k}$ ) 的 預測狀態。因此我們得到在時刻 $t_n, t_{n+1},...,t_{n+N}$ 對系統狀態的預測 (與狀態輸入序列 $u(0),...,u(N-1)$有關)。

每時刻的最佳控制力求解：
現在我們使用最佳控制的想法來決定 $u(0),...,u(N-1)$ 使得預測狀態 $x_u$ 能盡可能靠近 $x^{ref} = 0$。為了達成此目標，我們可選定 一距離函數
\[
l(x_u(k),u(k)) : = ||x_u(k)||^2 + \lambda ||u(k)||^2
\]其中 $||\cdot||$ 為 Euclidean norm 且 $\lambda \geq 0$為控制力的權重。則現在我們可以給出最佳控制問題如下
\[\min J\left( {x\left( n \right),u\left( . \right)} \right): = \sum\limits_{k = 0}^{N - 1} {l\left( {{x_u}\left( k \right),u\left( k \right)} \right)} \]對 $u(0),...,u(N-1)$ 為admissible 且 $x_u$ 滿足 (**)。

讓我們假設上述最佳控制問題具有解，且此解由一組最佳控制數列 記作
\[
u^*(0),...,u^*(N-1)
\]亦即我們有
\[\mathop {\min }\limits_{u\left( 0 \right),...,u\left( {N - 1} \right)} J\left( {x\left( n \right),u\left( . \right)} \right) = \sum\limits_{k = 0}^{N - 1} {l\left( {{x_{{u^*}}}\left( k \right),{u^*}\left( k \right)} \right)} \]現在為了讓控制力 $u$ 具有我們需要的回授形式 $\mu(x(n))$ 我們令
\[
\mu(x(n) ) := u^*(0)
\]亦即我們只使用最佳控制序列的第一元素。接著我們重複上述過程：在時刻 $t_{n+1}$ 可量得 $x(n+1)$ 並依此再度執行最佳化求得 $\mu(x(n+1))$。同理，對於時刻 $t_{n+2}$ 可量得 $x_{n+2}$ 並依此執行最佳化得到 $\mu(x(n+2))$...


Comments:
1. 因為受控系統為非線性，一般而言為我們所得到的 回授控制律 $\mu(\cdot)$ 必須透過迭代最佳化演算 而得，沒有解析解，故計算複雜度會是模型預測控制的一大挑戰。
2. 就預測模型觀點而言，我們可看出 預測狀態軌跡 $x_u(k), \;\; k=0,1,2,...,N$ 提供了對原本系統在已知時刻 $t_n$ 對於 $t_n,...,t_{n+N}$ 的預測，以及在時刻 $t_{n+1}$ 對於 $t_{n+1},...,t_{n+N+1}$ 的預測，以及 在時刻 $t_{n+2}$ 對於 $t_{n+2},...,t_{n+N+2}$ 的預測。因此我們可看出預測區間是 "移動" 的，故模型預測控制一般又稱為 移動區間控制 (Moving Horizon Control, MHC) 或稱 (Receding Horizon Control, RHC)
3. 若上述 $f$ 為線性，則上述討論稱之為  線性模型預測控制(Linear Model Predictive Control, LMPC)或者簡稱 模型預測控制 (Model Predictive Control, MPC)

[變分法] 離散泛函極值的必要條件

此文主要討論離散泛函的極值與其必要條件，也就是所謂的離散版本的 Euler-Lagrange Equation，推薦讀者可先複習先前介紹過的 連續泛函 的極值與必要條件的相關知識，整體推導而言可謂非常類似。

考慮離散泛函
\[\left\{ \begin{align*}
  &J\left( x \right): = \sum\limits_{k = 0}^{N - 1} {F\left( {x\left( k \right),x\left( {k + 1} \right),k} \right)} ; \hfill \\
  &x\left( {{0}} \right) = {x_0};x\left( {{N}} \right) = {x_1} \hfill \\
\end{align*}  \right.
\]其中 $F(x,y,t), \frac{{\partial F}}{{\partial x}}, \frac{{\partial F}}{{\partial y}}$ 在其定義域上連續函數。我們的目標是求序列 $x(0), x(1),...,x(N)$ 使得上述泛函 $J(x)$ 達到極值。

Comments:
前述設定中的離散狀態 $x(k) := x(t_k)$ 其中 $t_k = kT$ 且 $T$ 為取樣週期 (sampling period)

=====================================
Theorem: 離散版本的 Euler-Lagrange Equation
若 $x(1),...,x(N - 1) $ 使得上述泛函 $J(x)$ 達到極值，則對任意 $k=1,2,...,N-1$
\[\frac{{\partial F\left( {x(k),x(k + 1),k} \right)}}{{\partial x\left( k \right)}} + \frac{{\partial F\left( {x(k - 1),x(k),k - 1} \right)}}{{\partial x(k)}} =0\]=====================================

Proof: 令 $\delta x(k)$ 為 $x(k)$ 的變分，由於 $x(0) = x_0$ 與 $x(N)=x_1$，故 $\delta x(0) = \delta x(N) =0$ ，現在我們觀察 $J(x)$ ，由假設可知  $x(1),...,x(N - 1) $ 使得泛函 $J$ 達到極值，故
\[
J(x + \alpha \delta x) \geq J(x)
\]且此表明 $ J\left( {x + \alpha \delta x} \right)$ 在 $\alpha =0$ 處達到極值，由變分與泛函極值關係可知
\[
\delta J\left( {x\left( k \right)} \right) = {\left. {\frac{\partial }{{\partial \alpha }}J\left( {x\left( k \right) + \alpha \delta x\left( k \right)} \right)} \right|_{\alpha  = 0}} = 0
\]其中
\[J\left( {x(k) + \alpha \delta x(k)} \right) = \sum\limits_{k = 0}^{N - 1} {F\left( {x(k) + \alpha \delta x(k),x(k + 1) + \alpha \delta x(k + 1),k} \right)}
\]因此
\[
\delta J\left( {x\left( k \right)} \right) = {\left. {\frac{\partial }{{\partial \alpha }}\sum\limits_{k = 0}^{N - 1} {F\left( {x(k) + \alpha \delta x(k),x(k + 1) + \alpha \delta x(k + 1),k} \right)} } \right|_{\alpha  = 0}} = 0
\]故我們可推得
\[\begin{align*}
  &{\left. {\frac{\partial }{{\partial \alpha }}\sum\limits_{k = 0}^{N - 1} {F\left( {x(k) + \alpha \delta x(k),x(k + 1) + \alpha \delta x(k + 1),k} \right)} } \right|_{\alpha  = 0}} = 0 \hfill \\
 &  \Rightarrow {\left. {\sum\limits_{k = 0}^{N - 1} {\frac{\partial }{{\partial \alpha }}F\left( {x(k) + \alpha \delta x(k),x(k + 1) + \alpha \delta x(k + 1),k} \right)} } \right|_{\alpha  = 0}} \hfill \\
 &  \Rightarrow {\left. {\sum\limits_{k = 0}^{N - 1} {\frac{{\partial F}}{{\partial x\left( k \right)}}\delta x(k) + \frac{{\partial F}}{{\partial x(k + 1)}}\delta x(k + 1)} } \right|_{\alpha  = 0}} = 0  \;\;\;\; (*)
\end{align*}
\] 現在觀察上式的第二項，利用變數變換 定義 $k:=m-1$ 則我們可改寫為
\[\begin{gathered}
  \sum\limits_{k = 0}^{N - 1} {\frac{{\partial F\left( {x(k),x(k + 1),k} \right)}}{{\partial x(k + 1)}}\delta x(k + 1)}  = \sum\limits_{m = 1}^N {\frac{{\partial F\left( {x(m - 1),x(m),m - 1} \right)}}{{\partial x(m)}}\delta x(m)}  \hfill \\
   = \frac{{\partial F\left( {x(N - 1),x(N),N - 1} \right)}}{{\partial x(N)}}\delta x(N) + \sum\limits_{m = 1}^{N - 1} {\frac{{\partial F\left( {x(m - 1),x(m),m - 1} \right)}}{{\partial x(m)}}\delta x(m)}  \hfill \\
   = \frac{{\partial F\left( {x(N - 1),x(N),N - 1} \right)}}{{\partial x(N)}}\delta x(N) + \sum\limits_{k = 1}^{N - 1} {\frac{{\partial F\left( {x(k - 1),x(k),k - 1} \right)}}{{\partial x(k)}}\delta x(k)}  \hfill \\
\end{gathered} \]現在將上述結果代回 $(*)$，故可得
\[\small \begin{align*}
  \delta J\left( {x\left( k \right)} \right) &= \sum\limits_{k = 0}^{N - 1} {\frac{{\partial F}}{{\partial x\left( k \right)}}\delta x(k) + \frac{{\partial F}}{{\partial x(k + 1)}}\delta x(k + 1)}  \hfill \\
  & = \sum\limits_{k = 0}^{N - 1} {\frac{{\partial F}}{{\partial x\left( k \right)}}\delta x(k) + \frac{{\partial F\left( {x(N - 1),x(N),N - 1} \right)}}{{\partial x(N)}}\delta x(N) + \frac{{\partial F\left( {x(k - 1),x(k),k - 1} \right)}}{{\partial x(k)}}\delta x(k)}  \hfill \\
 &  = \frac{{\partial F\left( {x(N - 1),x(N),N - 1} \right)}}{{\partial x(N)}}\delta x(N) + \sum\limits_{k = 0}^{N - 1} {\left( {\frac{{\partial F}}{{\partial x\left( k \right)}} + \frac{{\partial F\left( {x(k - 1),x(k),k - 1} \right)}}{{\partial x(k)}}} \right)\delta x(k)}  \hfill \\
\end{align*}
\]注意到上式中 $\delta x(N) =0$ 且由於 $\delta x(k)$ 為任意變分，故由 $\delta J = 0$ 我們可知
\[\begin{align*}
 & \frac{{\partial F\left( {x(N - 1),x(N),N - 1} \right)}}{{\partial x(N)}}\delta x(N) +  \hfill \\
  \begin{array}{*{20}{c}}
  {}&{}
\end{array}&\;\;\;\; \sum\limits_{k = 0}^{N - 1} {\left( {\frac{{\partial F\left( {x(k),x(k + 1),k} \right)}}{{\partial x\left( k \right)}} + \frac{{\partial F\left( {x(k - 1),x(k),k - 1} \right)}}{{\partial x(k)}}} \right)\delta x(k)}  = 0
\end{align*} \]亦即對任意 $k=0,1,...,N-1$，
\[\frac{{\partial F\left( {x(k),x(k + 1),k} \right)}}{{\partial x\left( k \right)}} + \frac{{\partial F\left( {x(k - 1),x(k),k - 1} \right)}}{{\partial x(k)}} = 0\;\;\;\; \square
\]

[變分法] 連續泛函極值的必要條件

這次要介紹最簡單形式的 泛函極值問題的 必要條件，此條件一般又稱之為 Euler-Largrange Eqution。此方程可謂泛函極值的房角石，亦為之後在最佳控制理論中的最大值原理扮演開路先鋒，是極為重要的角色。在介紹之前，我們先做一般性的用語與基本性質介紹。


======================
Definition: 泛函
令 $\Omega$ 為 賦範函數空間 (normed function space)，若 對任意函數 $x(t) \in \Omega$ 都存在一個實數與之對應，則我們稱 $J$ 是定義在 $\Omega$ 上的 泛函 (functional)，記作 $J(x(t))$
======================

Comment:
1. 簡而言之，泛函 一詞即表示為由 函數空間 映射到 實數軸 上的函數 $J: \Omega \to \mathbb{R}$ 。
2. 再以下的討論中，集合 $\Omega$ 又稱為 泛函 $J$ 的 容許集 (admissible set)。


現取 $x_1, x \in \Omega$ 且 $\delta x := x_1 - x$ ，我們定義 關於 $\delta x$ 的 泛函增量 (increment) 如下
\[
\Delta J(\delta x)  := J(x_1) - J(x) =  J( x + \delta x) - J(x)
\]則由此 泛函增量，我們可以定義何謂泛函的變分。

======================
Definition: 泛函的變分
給定泛函 $J : \Omega \to \mathbb{R}$，若存在 一線性泛函 $L(x, \delta x)$ 使得泛函增量可被表為
\[
\Delta J(\delta x) = L(x, \delta x) + r(x, \delta x) \cdot | |\delta x||
\]其中 $r(x, \delta x)$ 為 其他高階剩餘項(remainder) 滿足 當 $| |\delta x|| \to 0 \Rightarrow r(x, \delta x) \to 0$，則我們稱上式中的 $L(x, \delta x)$ 為 $J(x)$ 的 變分 (variation)，記作 $\delta J := L(x, \delta x)$
======================

Comment:

1. 上述定義中的 線性泛函項 $L$ 與 其他高階剩餘項 $r$，可視為透過 Taylor 級數展開而得。

2. 變分 (variation) 一詞在文獻中又稱 differential

3. 若泛函變分存在，則該 變分 為唯一，在此不證明，有興趣讀者可參閱 [1]。
4. 關於線性泛函及其相關定義請讀者可參閱 [變分法] 淺論 線性泛函 
5. 有些文獻定義的泛函是透過所謂 Gateaux differentials 與 Freshet differential，但為求論述簡潔，在此不多作介紹，有興趣的讀者可以參閱 [2]

======================
Theorem: 泛函極值與變分關係
給定泛函 $J : \Omega \to \mathbb{R}$，若其變分存在，則 其變分可表為參數 $\alpha$ 的方向導數，亦即 變分滿足下式
\[\delta J(x(t)) = {\left. {\frac{\partial }{{\partial \alpha }}J(x(t) + \alpha \delta x(t))} \right|_{\alpha  = 0}}\]======================

Proof: 給定泛函 $J : \Omega \to \mathbb{R}$ 且假設其變分存在，我們要證明
\[\delta J(x(t)) = {\left. {\frac{\partial }{{\partial \alpha }}J(x(t) + \alpha \delta x(t))} \right|_{\alpha  = 0}}
\]首先由 $\delta J$ 存在可知：存在一線性泛函 $L$ 始得 泛函增量 $\Delta J$滿足
\[\begin{align*}
  \Delta J &= J\left( {x + \alpha \delta x} \right) - J\left( x \right) \hfill \\
   &= L(x,\alpha \delta x) + r(x,\alpha \delta x) \cdot || \alpha \delta x || \hfill \\
\end{align*}
\]由於 $L$ 為線性泛函，故 $L(x,\alpha \delta x) = \alpha L(x,\delta x)$，現在觀察
\[\begin{align*}
  {\left. {\frac{\partial }{{\partial \alpha }}J(x(t) + \alpha \delta x(t))} \right|_{\alpha  = 0}} &= \mathop {\lim }\limits_{\alpha  \to 0} \frac{{J(x + \alpha \delta x) - J\left( x \right)}}{\alpha }\\
   &= \mathop {\lim }\limits_{\alpha  \to 0} \frac{{L(x,\alpha \delta x) + r(x,\alpha \delta x)||\alpha \delta x||}}{\alpha } \hfill \\
   &= \mathop {\lim }\limits_{\alpha  \to 0} \frac{{L(x,\alpha \delta x)}}{\alpha }  + \underbrace {\mathop {\lim }\limits_{\alpha  \to 0} \frac{{r(x,\alpha \delta x) ||\alpha \delta x||}}{\alpha }}_{ = 0}  \hfill \\
   &= \mathop {\lim }\limits_{\alpha  \to 0} \frac{{\alpha L(x,\delta x)}}{\alpha }  \hfill \\
   &= L(x,\delta x) = \delta J(x) \;\;\;\;\; \square
\end{align*}
\]


======================
Theorem:
令 $J$ 為泛函且其變分存在，若 $J(x)$ 在 $x_0 \in \Omega$ 有(局部)極值，則其在 $x_0$ 之變分
\[
\delta J(x_0) =0
\] ======================

Comment: 上述定理中的 $x_0$ 又稱為 泛函 $J$ 的臨界點(critical point) 或者稱 不動點 (stationary point)。

Proof: 由於變分存在，我們可將變分用 以單變數參數 $\alpha$ 的方向導數表示
\[{\left. {\delta J\left( x \right) = \frac{\partial }{{\partial \alpha }}J(x + \alpha \delta x)} \right|_{\alpha  = 0}}\]由於 $J(x)$ 在 $x_0 \in \Omega$ 有局部極值，故我們可知 $\alpha =0$ 為$J(x_0 + \alpha \delta x)$ 的局部極值 (以極小值為例，可知對任意 $\alpha \in \mathbb{R}$， $J(x_0) \leq J(x_0 + \alpha \delta x)$，且極小值發生在 $\alpha = 0$)，故
\[{\left. {\delta J\left( {{x_0}} \right) = \frac{\partial }{{\partial \alpha }}J({x_0} + \alpha \delta x)} \right|_{\alpha  = 0}} = 0\;\;\;\; \square
\]


在討論一般設定之後，以下我們開始針對特殊形式的泛函來建構必要條件：考慮泛函
\[
J(x(t)) := \int_{t_0}^{t_1} F(t,x,\dot{x}) dt; \;\;\; x(t_0) :=x_0; \;\;\; x(t_1) \doteq x_1
\]且令其  admissible set 為
\[
\Omega := \{x(t) : x(t) \in C^2[t_0, t_1], \; x(t_0) = x_0, x(t_1) = x_1\}
\]且 $F(t, x, \dot{x})$ 為 $C^2$  (二階可導且連續)，我們欲求上述泛函極值的必要條件，此結果極為鼎鼎大名的 Euler-Largrange 方程，但在我們證明主要定理之前，底下我們先給個前置定理，此定理又稱為變分基本定理。


======================
Lemma: 變分基本引理
設函數 $F(t)$ 在區間 $[t_0, t_1]$ 上連續，若對於任意滿足 $\eta(t_0) = \eta(t_1) =0$ 的充分光滑函數 $\eta(t)$ 我們都有
\[
\int_{t_0}^{t_1} F(t) \eta(t) dt =0
\]則 $F(t) = 0$ 對 $t \in [t_0,t_1]$
======================

Proof: 利用反證法，假設 存在 $\xi \in (t_0,t_1)$ 使得 $F(\xi) \neq 0$，欲證明矛盾。不失一般性情況下我們假設 $F(\xi) >0$ 則由於 $F$的連續性，可知必存在 以 $\xi$ 為中心的鄰域 $N_\xi :=(\xi_1,\xi_2) \subset (t_0, t_1)$ 使得 對任意 $t \in N_\xi$，我們有 $F(t) > 0$。現在我們構造 $\eta(t)$ 函數如下
\[\eta \left( t \right): = \left\{ \begin{gathered}
  0,\begin{array}{*{20}{c}}
  {}&{}&{}&{}
\end{array}t \in \left[ {{t_0},{\xi _1}} \right) \hfill \\
  {\left[ {\left( {t - {\xi _1}} \right)\left( {t - {\xi _2}} \right)} \right]^2},\begin{array}{*{20}{c}}
  {}&{}
\end{array}t \in \left[ {{\xi _1},{\xi _2}} \right] \hfill \\
  0,\begin{array}{*{20}{c}}
  {}&{}&{}&{}
\end{array}t \in \left( {{\xi _2},{t_1}} \right] \hfill \\
\end{gathered}  \right.\]且注意到上述 $\eta(t)$ 函數滿足 $\eta(t_0) = \eta(t_1) = 0$ 且為連續函數，然而若我們觀察
\[
\int_{t_0}^{t_1} F(t) \eta(t) dt = \int_{\xi_1}^{\xi_2} F(t) \eta(t) dt > 0
\]此結果與我們的假設矛盾。$\square$



======================
Theorem: 泛函極值的必要條件 Euler-Lagrange Equation
設函數 $F(t, x, \dot{x})$ 具有連續二階偏導數，且設泛函\[
J(x(t)) := \int_{t_0}^{t_1} F(t,x,\dot{x}) dt; \;\;\; x(t_0) :=x_0; \;\;\; x(t_1) \doteq x_1
\]在 $x(t) \in \Omega$ 達到極值，則 $x(t)$ 滿足下列方程
\[\frac{\partial }{{\partial x}}F\left( {t,x,\dot x} \right) - \frac{d}{{dt}}\left( {\frac{\partial }{{\partial \dot x}}F\left( {t,x,\dot x} \right)} \right) = 0\]
======================

Proof: 首先令 $\phi(t) := \delta x(t)$，則由於 $x(t_0)=x_0$與 $x(t_1) = x_1$ 可知，$\phi(t)$ 滿足 $\phi(t_0) = \phi(t_1)=0$，現在由泛函極值與變分的關係可知下式必定成立：
\[\delta J\left( {x\left( t \right)} \right) = \left. {\frac{\partial }{{\partial \alpha }}J\left( {x\left( t \right) + \alpha \phi \left( t \right)} \right)} \right|_{\alpha = 0} = 0 \;\;\;(\star)
\]現在觀察
\[
J\left( {x\left( t \right) + \alpha \phi \left( t \right)} \right) = \int_{{t_0}}^{{t_1}} F (t,x + \alpha \phi ,\dot x + \alpha \dot \phi )dt
\]故我們可先行計算
\[{\left. {\frac{\partial }{{\partial \alpha }}J\left( {x\left( t \right) + \alpha \phi \left( t \right)} \right)} \right|_{\alpha  = 0}} = {\left. {\frac{\partial }{{\partial \alpha }}\int_{{t_0}}^{{t_1}} F (t,x + \alpha \phi ,\dot x + \alpha \dot \phi )dt} \right|_{\alpha  = 0}}
\]由 Libneiz Rule 可得
\[\begin{align*}
  {\left. {\frac{\partial }{{\partial \alpha }}J\left( {x\left( t \right) + \alpha \phi \left( t \right)} \right)} \right|_{\alpha  = 0}} &= {\left. {\frac{\partial }{{\partial \alpha }}\int_{{t_0}}^{{t_1}} F (t,x + \alpha \phi ,\dot x + \alpha \dot \phi )dt} \right|_{\alpha  = 0}} \hfill \\
   &= {\left. {\int_{{t_0}}^{{t_1}} {\frac{\partial }{{\partial \alpha }}F} (t,x + \alpha \phi ,\dot x + \alpha \dot \phi )dt} \right|_{\alpha  = 0}} \hfill \\
   &=  {\left. {\int_{{t_0}}^{{t_1}} {\left[ {\frac{{\partial F}}{{\partial x}}\phi  + \frac{{\partial F}}{{\partial \dot x}}\dot \phi } \right]} dt} \right|_{\alpha  = 0}}\;\;\;\; (*)
\end{align*}
\]注意到上述積分第二項可透過 integration by part 求得
\[\int_{{t_0}}^{{t_1}} {\frac{{\partial F}}{{\partial \dot x}}\dot \phi dt}  = \left. {\frac{{\partial F}}{{\partial \dot x}}\phi } \right|_{{t_0}}^{{t_1}} - \int_{{t_0}}^{{t_1}} {\phi \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} dt\]由於 $\phi(t_0) = \phi(t_1) = 0$，故我們得
\[\begin{gathered}
  \int_{{t_0}}^{{t_1}} {\frac{{\partial F}}{{\partial \dot x}}\dot \phi dt}  = \underbrace {\left. {\frac{{\partial F}}{{\partial \dot x}}\phi } \right|_{{t_0}}^{{t_1}}}_{ = 0} - \int_{{t_0}}^{{t_1}} {\phi \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} dt \hfill \\
   \Rightarrow \int_{{t_0}}^{{t_1}} {\frac{{\partial F}}{{\partial \dot x}}\dot \phi dt}  =  - \int_{{t_0}}^{{t_1}} {\phi \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} dt \hfill \\
\end{gathered}
\]現在將其帶回 $(*)$ 我們得到
\[\begin{align*}
  {\left. {\frac{\partial }{{\partial \alpha }}J\left( {x\left( t \right) + \alpha \phi \left( t \right)} \right)} \right|_{\alpha  = 0}}
   &= {\left. {\int_{{t_0}}^{{t_1}} {\left[ {\frac{{\partial F}}{{\partial x}}\phi  - \phi \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} \right]} dt} \right|_{\alpha  = 0}} \hfill \\
   &= {\left. {\int_{{t_0}}^{{t_1}} {\left[ {\frac{{\partial F}}{{\partial x}} - \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} \right]} \phi dt} \right|_{\alpha  = 0}} \hfill \\
\end{align*}
\]由於 $(\star)$ 可知，
\[\begin{align*}
  {\left. {\frac{\partial }{{\partial \alpha }}J\left( {x\left( t \right) + \alpha \phi \left( t \right)} \right)} \right|_{\alpha  = 0}}
&= {\left. {\int_{{t_0}}^{{t_1}} {\left[ {\frac{{\partial F}}{{\partial x}}\phi  + \frac{{\partial F}}{{\partial \dot x}}\dot \phi } \right]} dt} \right|_{\alpha  = 0}} \hfill \\
   &= {\left. {\int_{{t_0}}^{{t_1}} {\left[ {\frac{{\partial F}}{{\partial x}}\phi  - \phi \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} \right]} dt} \right|_{\alpha  = 0}} \hfill \\
   &= {\left. {\int_{{t_0}}^{{t_1}} {\left[ {\frac{{\partial F}}{{\partial x}} - \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} \right]} \phi dt} \right|_{\alpha  = 0}} = 0 \hfill \\
\end{align*}
\]由於 ${\frac{{\partial F}}{{\partial x}} - \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}}$ 在區間 $[t_0,t_1]$ 連續，且 $\phi$ 滿足 $\phi(t_0) = \phi(t_1) =0$ 且 $\phi \in C^2$，利用前述引理可知在 $[t_0,t_1]$ 上，
\[{\frac{{\partial F}}{{\partial x}} - \frac{d}{{dt}}\frac{{\partial F}}{{\partial \dot x}}} = 0\;\;\;\; \square\]


[1] I. M. Gelfand and S. V. Fomin, Calculus of Variations, 2000
[2] David G. Luenberger, Optimization By Vector Space Methods, 1997

[最佳控制] 離散時間 穩態 LQR 控制問題 (2) - LQR convergence

回憶離散時間穩態LQR問題： 定義 cost function
\[\begin{array}{l}
V\left( {x,{\bf{u}}} \right): = \frac{1}{2}\sum\limits_{k = 0}^\infty  {{x^T}\left( k \right)Qx\left( k \right) + {u^T}\left( k \right)Ru\left( k \right)} \\
s.t.\\
x\left( {k + 1} \right) = Ax\left( k \right) + Bu\left( k \right)\\
x\left( 0 \right) = x
\end{array}\]且 $Q, R >0$ 為正定矩陣。若 $(A,B)$ 為可控制 則 最佳化問題
\[
\min_{\bf u} V(x,{\bf u})
\]解存在 且唯一 (why?)。

現在我們將此解 (控制力序列) 記做 ${\bf u}^*(x)$；其中第一項為 $u^*(x)$。且對應的最佳控制律 $K_{\infty}(\cdot)$  故我們有
\[
K_\infty(x) = u^*(x) = {\bf u}^* (0;x)
\]

以下定理說明若可控制系統 採用 穩態 LQR 控制，則保證閉迴路穩定度：

=====================
Theorem: Steady State LQR Guarantees the Closed-Loop Stability
對 $(A,B)$ 可控制，穩態LQR問題 且 $Q,R >0$ 保證閉迴路系統
\[
x(k+1) = A x(k) + B K_\infty (x)
\]穩定。
=====================

Comment:
線性系統中，漸進穩定度 等價 漸進收斂 $x(k) \to 0$

Proof:
首先證明 穩態 cost function 有界。注意到由於  $(A,B)$ 可控制，由定義可知存在 一組控制力
\[
{\bf u} = \{u(0), u(1),...,u(n-1)\}
\] 使得 控制後的系統能在有限時間 $n$ 步之內，可由任意給定初始狀態 $x(0)$ 移動到 給定的終止狀態 $x(n) = 0$；且 在時間 $ n$ 之後的控制力 $ \{u(n+1), u(n+2),...\} = \{0,0,0,...\}$ 故原本無窮和的 cost 函數變為有限和
\[\begin{array}{l}
V\left( {x,{\bf{u}}} \right): = \frac{1}{2}\sum\limits_{k = 0}^\infty  {{x^T}\left( k \right)Qx\left( k \right) + {u^T}\left( k \right)Ru\left( k \right)} \\
 \Rightarrow V\left( {x,{\bf{u}}} \right) = \frac{1}{2}\sum\limits_{k = 0}^n {{x^T}\left( k \right)Qx\left( k \right) + {u^T}\left( k \right)Ru\left( k \right)}
\end{array}\]且由於 成本函數為 strictly convex in $\bf u$ 且 $R>0$ ( $u$ 沒有 vanish) 故 此穩態LQR最佳問題的解為唯一。

現在我們觀察 costs to go 的數列滿足
\[
V_{k+1} = V_k - 1/2 (x(k)'Qx(k) + u(k)' R u(k))
\]其中 $V_k := V^0(x(k))$ 為 在時刻 $k$ 對應 狀態為 $x(k)$ 與 對應最佳控制力 $u(k) = u^0(x(k))$ 的 cost。故此數列 $\{V_k\}$ 為非遞增 且 有下界 (為 $0$)。可推知此數列必定收斂；故
\[
|V_{k+1} - V_k| \to 0\;\;\;\; \text{as $k \to \infty$}
\]因此
\[\left\{ \begin{array}{l}
x\left( k \right)'Qx\left( k \right)\; \to 0\\
u\left( k \right)'{\rm{ }}R{\rm{ }}u\left( k \right) \to 0
\end{array} \right.\]又因為 $Q,R >0$ 故可推知
\[\left\{ \begin{array}{l}
x\left( k \right)\; \to 0\\
u\left( k \right) \to 0
\end{array} \right.\]亦即閉迴路系統狀態收斂 (閉迴路穩定!)。 $\square$

Comments 
1. $R>0$ 條件用以保證控制律唯一性。實際上若使用者不關心控制律，可將其設為參數非常小的正定矩陣
2. 上述定裡假設可從 controllability 推廣到 stabilizability。
3. $Q>0$ 的條件可被推廣到 $Q \ge 0$ 與 $(A,Q)$ detectable。

[最佳控制] 線性系統的最佳參數估計 Kalman Filter (0)- 預備知識

這次要介紹 Kalman filter 或稱 Optimal Linear State Estmator，以下我們將簡單介紹一些在下一篇文章需要使用的一些結果。

Preliminary
回憶若 $x$ 為 mean $m$ 且 variance $\sigma^2$ 的 Normal 隨機變數 則 我們表示為 $x \sim N(m, \sigma^2)$ 。其 probability density function 可記做
\[
\frac{1}{\sqrt{2 \pi} \sigma} \exp(\frac{-1}{2 \sigma^2}(x-m)^2)
\]現在若我們拓展上述結果到 多個隨機變數 (又稱 random vector 隨機向量) 的情況，令 $x$ 為Normal 隨機向量 記做
\[
x \sim N(m,P); \;\;
p_x(x) := n(x,m,P)
\]上述符號 表示 $x$ 為 normal distributed 且 mean vector $m$ 與 convariance matrix $P$。另外 $n(x,m,P)$ 表示 normal probability density function
\[\;n\left( {x,m,P} \right): = \frac{1}{{{{\left( {2\pi } \right)}^{n/2}}{{\left( {\det P} \right)}^{1/2}}}}\exp \left( { - \frac{1}{2}{{\left( {x - m} \right)}^T}{P^{ - 1}}\left( {x - m} \right)} \right)\]
Comment:
若 $x \in \mathbb{R}^n$ 則 mean vector $m \in \mathbb{R}^n$ 且 convariance matrix $P \in \mathbb{R}^{n \times n}$ 且為 實數 對稱 正定 矩陣 (正定條件用以確保 $P^{-1}$ 存在，使得上述的 probability density function 可以被定義)。若 $P$ 不為對 正定 我們稱為 singular normal distribution 或稱 degenerate normal。

Example 
若 $n=2$ 則我們可以繪製 normal density function； 比如說
\[m = \left[ {\begin{array}{*{20}{c}}
0\\
0
\end{array}} \right];\begin{array}{*{20}{c}}
{}&{}
\end{array}{P^{ - 1}} = \left[ {\begin{array}{*{20}{c}}
{3.5}&{2.5}\\
{2.5}&{4.0}
\end{array}} \right]\]則我們可繪製

現在我們看幾個 之後會使用到的基本結果：

================
FACT 1: Joint independent normals
若隨機向量 $x  \sim N(m_x,P_x)$ 與 $y  \sim N(m_y,P_y)$ 為 normally distributed 且 彼此互為獨立，則其 joint density $p_{x,y}(x,y)$ 如下
\[
p_{x,y}(x,y) = p_{x}(x)p_{y}(y)=n(x,m_x,P_x) \cdot n(y,m_y,P_y)
\]且
\[\left[ {\begin{array}{*{20}{c}}
x\\
y
\end{array}} \right] \sim N\left( {\left[ {\begin{array}{*{20}{c}}
{{m_x}}\\
{{m_y}}
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{P_x}}&0\\
0&{{P_y}}
\end{array}} \right]} \right)\]================


================
FACT 2: Linear Transformation of a normal
若 $x \sim N(m, P)$ 且 $y$ 為 $x$ 的線性轉換；亦即對任意矩陣 $A$， $y = Ax$ 則
\[
y \sim N(Am , APA^T)
\]================

================
FACT 3: Conditional of a joint normal
若 $x,y$ 為 jointly normal distributed (no independent assumption)
\[\left[ {\begin{array}{*{20}{c}}
x\\
y
\end{array}} \right]\sim N\left( {\left[ {\begin{array}{*{20}{c}}
{{m_x}}\\
{{m_y}}
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)\]則 conditional density of $x$ given $y$ 仍為 Normal 亦即
\[
(x|y) \sim N(m,P)
\]其 probability density function 為 $ p_{x|y} (x|y) = n(x,m,P)$ 其中 conditional mean vector $m$ 與 conditional convariance matrix $P$ 分別為
\[\begin{array}{l}
m = {m_x} + {P_{xy}}P_y^{ - 1}(y - {m_y})\\
P = {P_x} - {P_{xy}}P_y^{ - 1}{P_{yx}}
\end{array}\]================
注意到上式中 conditional mean $m$ 為 random vector (depends on $y$) 。

Proof:
回憶 conditional density of $x$ given y 定義
\[
p_{x|y}(x,y) := \frac{p_{x,y}(x,y)}{p_y(y)}
\]注意到 $(x,y)$ 為 joint normal 故我們有
\[\small \begin{array}{l}
{p_y}\left( y \right): = \frac{1}{{{{\left( {2\pi } \right)}^{{n_y}/2}}{{\left( {\det {P_y}} \right)}^{1/2}}}}\exp \left( { - \frac{1}{2}{{\left( {y - {m_y}} \right)}^T}{P_y}^{ - 1}\left( {y - {m_y}} \right)} \right)\\
{p_{x,y}}\left( {x,y} \right): = \frac{1}{{{{\left( {2\pi } \right)}^{\left( {{n_x} + {n_y}} \right)/2}}{{\left( {\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)}^{1/2}}}}\exp \left( { - \frac{1}{2}{{\left[ {\begin{array}{*{20}{c}}
{x - {m_x}}\\
{y - {m_y}}
\end{array}} \right]}^T}{{\left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]}^{ - 1}}\left[ {\begin{array}{*{20}{c}}
{x - {m_x}}\\
{y - {m_y}}
\end{array}} \right]} \right)
\end{array}\]亦即
\[\begin{array}{l}
{p_{x|y}}\left( {x|y} \right) = \frac{{{p_{x,y}}\left( {x,y} \right)}}{{{p_y}\left( y \right)}}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \frac{{\frac{1}{{{{\left( {2\pi } \right)}^{\left( {{n_x} + {n_y}} \right)/2}}{{\left( {\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)}^{1/2}}}}\exp \left( { - \frac{1}{2}{{\left[ {\begin{array}{*{20}{c}}
{x - {m_x}}\\
{y - {m_y}}
\end{array}} \right]}^T}{{\left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]}^{ - 1}}\left[ {\begin{array}{*{20}{c}}
{x - {m_x}}\\
{y - {m_y}}
\end{array}} \right]} \right)}}{{\frac{1}{{{{\left( {2\pi } \right)}^{{n_y}/2}}{{\left( {\det {P_y}} \right)}^{1/2}}}}\exp \left( { - \frac{1}{2}{{\left( {y - {m_y}} \right)}^T}{P_y}^{ - 1}\left( {y - {m_y}} \right)} \right)}}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \frac{{{{\left( {\det {P_y}} \right)}^{1/2}}\exp \left( { - \frac{1}{2}\left( {{{\left[ {\begin{array}{*{20}{c}}
{x - {m_x}}\\
{y - {m_y}}
\end{array}} \right]}^T}{{\left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]}^{ - 1}}\left[ {\begin{array}{*{20}{c}}
{x - {m_x}}\\
{y - {m_y}}
\end{array}} \right] - {{\left( {y - {m_y}} \right)}^T}{P_y}^{ - 1}\left( {y - {m_y}} \right)} \right)} \right)}}{{{{\left( {2\pi } \right)}^{{n_x}/2}}{{\left( {\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)}^{1/2}}}}
\end{array}\]注意到若我們取 $P:= P_x - P_{xy}P_y^{-1}P_{yx}$ 則 利用 Matrix inversion Lemma 可知
\[{\left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]^{ - 1}} = \left[ {\begin{array}{*{20}{c}}
{{P^{ - 1}}}&{ - {P^{ - 1}}{P_{xy}}P_y^{ - 1}}\\
{ - P_y^{ - 1}{P_{yx}}{P^{ - 1}}}&{P_y^{ - 1} + P_y^{ - 1}{P_{yx}}{P^{ - 1}}{P_{xy}}P_y^{ - 1}}
\end{array}} \right]\]將此結果帶回我們可得
\[\begin{array}{l}
{p_{x|y}}\left( {x|y} \right) = \frac{{{p_{x,y}}\left( {x,y} \right)}}{{{p_y}\left( y \right)}}\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \frac{{{{\left( {\det {P_y}} \right)}^{1/2}}\exp \left( { - \frac{1}{2}\left( \begin{array}{l}
{\left( {x - {m_x}} \right)^T}{P^{ - 1}}\left( {x - {m_x}} \right) - 2{\left( {y - {m_y}} \right)^T}P_y^{ - 1}{P_{yx}}{P^{ - 1}}\left( {x - {m_x}} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} + {\left( {y - {m_y}} \right)^T}\left( {P_y^{ - 1}{P_{yx}}{P^{ - 1}}{P_{xy}}P_y^{ - 1}} \right)\left( {y - {m_y}} \right)
\end{array} \right)} \right)}}{{{{\left( {2\pi } \right)}^{{n_x}/2}}{{\left( {\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)}^{1/2}}}}\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \frac{{{{\left( {\det {P_y}} \right)}^{1/2}}\exp \left( { - \frac{1}{2}\left( {{{\left[ {{{\left( {x - {m_x}} \right)}^T} - {P_{xy}}P_y^{ - 1}\left( {y - {m_y}} \right)} \right]}^T}{P^{ - 1}}\left[ {{{\left( {x - {m_x}} \right)}^T} - {P_{xy}}P_y^{ - 1}\left( {y - {m_y}} \right)} \right]} \right)} \right)}}{{{{\left( {2\pi } \right)}^{{n_x}/2}}{{\left( {\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)}^{1/2}}}}\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \frac{{{{\left( {\det {P_y}} \right)}^{1/2}}\exp \left( { - \frac{1}{2}\left( {{{\left[ {\left( {x - {m_x}} \right) - {P_{xy}}P_y^{ - 1}\left( {y - {m_y}} \right)} \right]}^T}{P^{ - 1}}\left[ {\left( {x - {m_x}} \right) - {P_{xy}}P_y^{ - 1}\left( {y - {m_y}} \right)} \right]} \right)} \right)}}{{{{\left( {2\pi } \right)}^{{n_x}/2}}{{\left( {\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)}^{1/2}}}}
\end{array}\]接著令 $m := m_x + P_{xy} P_y^{-1}(y-m_y)$ 可得
\[{p_{x|y}}\left( {x|y} \right) = \frac{{{p_{x,y}}\left( {x,y} \right)}}{{{p_y}\left( y \right)}} = \frac{{{{\left( {\det {P_y}} \right)}^{1/2}}\exp \left( { - \frac{1}{2}{{\left( {x - m} \right)}^T}{P^{ - 1}}\left( {x - m} \right)} \right)}}{{{{\left( {2\pi } \right)}^{{n_x}/2}}{{\left( {\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)}^{1/2}}}}\]注意到
\[\det \left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right] = \det {P_y}\det P\]故可得
\[{p_{x|y}}\left( {x|y} \right) = \frac{1}{{{{\left( {2\pi } \right)}^{{n_x}/2}}{{\left( {\det P} \right)}^{1/2}}}}\exp \left( { - \frac{1}{2}{{\left( {x - m} \right)}^T}{P^{ - 1}}\left( {x - m} \right)} \right) = n(x,m,P) \ \ \ \ \ \square
\]


如果要推導 最佳估測器 我們需要上述結果衍生：

================
FACT 1': Joint independent normals
若 $p_{x|z} (x|z) = n(x, m_x, P_x)$ 為 normal ，令 $y \sim N(m_y, P_y)$ 且 與 $x,z$ 彼此獨立 則 conditional joint density of $(x,y)$ given $z$ 為
\[{p_{x,y|z}}\left( {\left[ {\begin{array}{*{20}{c}}
x\\
y
\end{array}} \right]|z} \right) = n\left( {\left[ {\begin{array}{*{20}{c}}
x\\
y
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{m_x}}\\
{{m_y}}
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{P_x}}&0\\
0&{{P_y}}
\end{array}} \right]} \right)\]================

================
FACT 2': Linear Transformation of a Normal
若 $p_{x|z}(x|z)=  n(x,m, P)$ 且 $y$ 為 $x$ 的線性轉換；亦即 $y = Ax$ 則
\[{p_{y|z}}(y|z) = n(y,Am,AP{A^T})\]================

================

FACT 3': Conditional of a joint normal
若 $x,y$ 為 jointly normal distributed (no independent assumption)
\[{p_{x,y|z}}\left( {\left[ {\begin{array}{*{20}{c}}
x\\
y
\end{array}} \right]|z} \right) = n\left( {\left[ {\begin{array}{*{20}{c}}
x\\
y
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{m_x}}\\
{{m_y}}
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{P_x}}&{{P_{xy}}}\\
{{P_{yx}}}&{{P_y}}
\end{array}} \right]} \right)\]則 conditional density of $x$ given $y,z$ 仍為 Normal ，記為
\[
p_{x|y,z} (x|y,z) = n(x,m,P)
\]其中 conditional mean vector $m$ 與 conditional convariance matrix $P$ 分別為
\[\begin{array}{l}
m = {m_x} + {P_{xy}}P_y^{ - 1}(y - {m_y})\\
P = {P_x} - {P_{xy}}P_y^{ - 1}{P_{yx}}
\end{array}\]
================

Proof:
由 conditional density of $x$ given $y,z$ 定義可知
\[{p_{x|y,z}}\left( {x|y,z} \right) = \frac{{{p_{x,y,z}}\left( {x,y,z} \right)}}{{{p_{y,z}}\left( {y,z} \right)}}\]現在對等號右方 分子分母同乘 $p(z)$ 可得
\[\begin{array}{l} {p_{x|y,z}}\left( {x|y,z} \right) = \frac{{{p_{x,y,z}}\left( {x,y,z} \right)}}{{{p_{y,z}}\left( {y,z} \right)}}\\ \begin{array}{*{20}{c}} {}&{}&{} \end{array} = \frac{{{p_{x,y,z}}\left( {x,y,z} \right)}}{{{p_z}\left( z \right)}}\frac{{{p_z}\left( z \right)}}{{{p_{y,z}}\left( {y,z} \right)}}\\ \begin{array}{*{20}{c}} {}&{}&{} \end{array} = {p_{x,y|z}}\left( {x,y|z} \right) \cdot \frac{1}{{{p_{y|z}}\left( {y|z} \right)}} \end{array}\]則由先前 FACT 3 可計算 $p(y,z)$ 並且帶入 $p(x,y|z)$ 即可求得所求。$\square$


 ref: J. B. Rawlings and D. Q. Mayne, "Model Predictive Control: Theory and Design".

[變分法] 淺論 線性泛函

這次要介紹一些基本的 Functional

我們首先定義 Continuous Functional

給定一個 normed linear space $X$
============================
Definition:  Continuous Functional

考慮 一個 functional $J: X \rightarrow \mathbb{R}$, 令 $y \in X$，我們說 $J[y]$ 被稱作 在點 $\hat y \in X$ continuous 若下列條件成立

對任意 $\varepsilon >0$, 存在 $\delta >0$ 使得

\[||y- \hat y|| < \delta  \Rightarrow |J[y] - J[\hat y]| < \varepsilon
\]============================

Comment
上述連續性等價為
\[
||y(x) - \hat y(x)|| \rightarrow 0 \Rightarrow |f(y(x)) - f(\hat{y}(x))| \rightarrow 0
\]

============================
Definition:  Linear Functional
給定一個 normed linear space $X$，且 $y \in X$ 的元素，現在定義 $J[y] : X \rightarrow \mathbb{R}$ 為在 $X$ 上的 functional ，則我們說 $J[y]$ 為 linear functional 若下列條件成立
1. 對任意 $y\in X$ 與 $ \alpha \in \mathbb{R}$，$J[\alpha y] = \alpha J[y]$
2. 對任意 $y_1, y_2 \in X$，$J[y_1 + y_2] = J[y_1] + J[y_2]$
============================

Comment:
上述定義的兩個條件可簡化為
對任意  $y_1, y_2 \in X$ 與 $ \alpha, \beta \in \mathbb{R}$
\[
J[\alpha y_1 + \beta y_2] =\alpha  J[y_1] + \beta J[y_2]
\]

============================

Definition: Continuous Linear Functional
我們稱 $J[y]$ 為 continuous linear functional 若 $J[y]$ 為 linear functional，且 對任意 $h \in X$， $J[y]$ 為 連續。
============================

現在我們看一些例子：
Example 1
令 $X := \cal{C}^1[0,1]$，$y:[0,1] \rightarrow \mathbb{R}$，且考慮 $||x||:=||x||_{\infty}$現考慮
\[
f(y) := \frac{d}{dx} y(0)
\]則 此 $f(y)$ 為 Linear Functional 但不為連續。

Proof:
線性：
\[\begin{array}{l}
f(\alpha {y_1} + \beta {y_2}) = \frac{d}{{dx}}\left[ {\alpha {y_1}\left( 0 \right) + \beta {y_2}\left( 0 \right)} \right]\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \alpha \frac{d}{{dx}}{y_1}\left( 0 \right) + \beta \frac{d}{{dx}}{y_2}\left( 0 \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \alpha f({y_1}) + \beta f({y_2})
\end{array}
\]

接著考慮連續性
我們只要舉出一個反例即可說明其不連續
令 $\hat y =0$ 且 $y: = A  \cdot \sin \frac{x}{A}$, $(A \in \mathbb{R})$，則由 $f$ 定義
\[\begin{array}{l}
f(y) = \frac{d}{{dx}}y(0)\\
 \Rightarrow f(y) = \frac{d}{{dx}}{\left. {\left( {A \cdot \sin \frac{x}{A}} \right)} \right|_{x = 0}} = {\left. {\cos \frac{x}{A}} \right|_{x = 0}} = 1
\end{array}
\] 現在觀察連續性，我們需要
\[
||y(x) - \hat y(x)|| \rightarrow 0 \Rightarrow |f(y(x)) - f(\hat{y}(x))| \rightarrow 0
\]
故現在計算
\[{\left\| y \right\|_\infty } = \left| A \right|{\left\| {\sin \frac{x}{A}} \right\|_\infty } = \left| A \right|
\]現若讓 $A \rightarrow 0 $ 則 $||y|| \rightarrow 0$
但是此時對應的
 \[|f(y(x)) - f(\hat y(x))| = |f(y(x)) - f(0)| = |f(y(x))| = 1 \neq 0
\] 故此說明了 functional $f$ 在 $0$ 處不連續。 $\square$


Example 2
令 $X := \cal{C}[a,b]$，$ y : [a,b] \rightarrow \mathbb{R}$，現考慮下列積分
\[
J[y] := \int_a^b y(x)dx
\] Claim: 上述積分為 Linear Functional on $\cal{C}[a,b]$
Proof
1. 上式積分 為 Functional 因為 $J: \cal{C}[a,b] \rightarrow \mathbb{R}$，
2. 檢驗線性：
觀察
\[\begin{array}{l}
J[\alpha {y_1} + \beta {y_2}] = \int_a^b {\left( {\alpha {y_1}\left( x \right) + \beta {y_2}\left( x \right)} \right)} dx\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \alpha \int_a^b {{y_1}\left( x \right)} dx + \beta \int_a^b {{y_2}\left( x \right)} dx\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \alpha J[{y_1}] + \beta J[{y_2}]
\end{array}
\]故 上式積分確實為 Linear Functional。 $\square$

Example 2
令 $X := \cal{C}^n[a,b]$考慮下列積分
\[J[y]: = \int_a^b {\left[ {{\alpha _0}\left( x \right)y(x) + {\alpha _1}\left( x \right)y'(x) + ... + {\alpha _n}\left( x \right){y^{\left( n \right)}}(x)} \right]} dx
\]上式亦為 Linear functional 其中 $\alpha_i(x)$ 為 $\cal{C}[a,b]$ 上固定函數。

注意到 對任意 $y(x) $ 在特定 function space，現若讓上述積分 $=0$，則我們想問 對於 $\alpha_i(x)$ 會發生甚麼事情？

下面的 Lemma 可以回答此問題：
=======================
Lemma: 
若 $\alpha(x) \in \cal{C}[a,b]$ 且 對任意 $y(x) \in \cal{C}[a,b]$ 滿足 $y(a) = y(b) =0$ ，
\[
J[y]: = \int_a^b {\alpha \left( x \right)y(x)} dx = 0
\] 則對任意 $x \in [a,b]$
\[
\alpha(x) \equiv 0
\]=======================

Proof
利用歸謬法(Suppose toward to contradiction)，假設 存在 $x \in [a,b]$ 使得 $\alpha(x) \neq 0$。在不失一般性的情況下我們可設 $\alpha(x) >0$

現由於 $\alpha(x) \in \cal{C}[a,b]$ ，由連續性可知必存在一區間 $[x_1,x_2] \subset [a,b]$ 使得 對 $x \in [x_1,x_2]$
\[
\alpha(x) >0
\]
現在我們讓
\[y(x): = \left\{ \begin{array}{l}
(x - {x_1})({x_2} - x)\begin{array}{*{20}{c}}
{}
\end{array},x \in \left[ {{x_1},{x_2}} \right]\\
0\begin{array}{*{20}{c}}
{}&{}
\end{array},o.w.
\end{array} \right.
\] 則此 $y(x)$亦滿足我們假設條件 $y(x) \in \cal{C}[a,b]$ 滿足 $y(a) = y(b) =0$，但
\[J[y]: = \int_a^b {\alpha \left( x \right)y(x)} dx = \int_{{x_1}}^{{x_2}} {\alpha \left( x \right)(x - {x_1})({x_2} - x)} dx > 0\] 與假設 $J[y] =0$ 矛盾。 $\square$

[變分法] 基本變分問題

這次要與大家介紹 變分法(Variational Calculus) ，在一般傳統數學分析求極值時，微積分 (Calculus)處理的對象為函數 (Function)的極值問題，而變分法 (Variational Calculus) 所處理的對象為泛函 (Functional) 的極值問題。

所謂泛函通常是指一種 domain 為函數空間(亦即無窮維的向量空間)，而 codmain 為實數 或者 Euclidean 空間的 函數。簡而言之，泛函可視為 函數的函數。

我們給出泛函定義如下：
==================
Definition: Functional
令 $X$ 為任意 Vector Space，我們稱函數 $J: X \to \mathbb{R}$ 為一個 泛函 (Functional)
==================

Comments:
考慮 $X, Y$ 為 Vector Space，則
1. $ g: X \rightarrow \mathbb{R}$ 為泛函
2. $ g: X \rightarrow \mathbb{R}^n$ 為泛函
3. $ g: X \rightarrow Y$ :此稱為 operator 不稱為泛函


以下我們給出幾個 functional 例子：


Examples of Functional
0. 給定 $X$ 任意 賦範空間，則對任意 $x \in X$，定義 $f(x):=||x||$ 為一個 functional。

1.考慮 $y(x)$ 為定義在 $[a,b]$上 一階連續可微函數，則下式
\[
J[y] = \int_a^b y'^2(x)dx
\] 為一個 Functional

2. 考慮一平面上任兩點 $A,B$ ，現在假設有一 質點(particle) 具有固定速度 $v$ ，且此 質點 可以沿著任意平面上任意路徑從 $A$ 移動到 $B$，則如果想描述此質點 花多少時間來通過上述的(任意)路徑，則描述結果會是用一個積分 以 Functional  來表示。

3. 令 $F(\alpha, \beta, \gamma)$ 為三變數的連續函數，則下式
\[
J[y] = \int_a^b F[x,y(x),y'(x)]dx
\] 為一個 Functional，其中 $y(x)$ 定義在 $[a,b]$ 上一階連續可微函數。


基本變分問題
整個變分法主要處理的問題為試圖 "找出 某 Functional 的極值"，下面是一些經典的變分基本問題：

1. [最短曲線問題] 
找出一曲線 $y = y(x)$ 使得
\[
\int_a^b \sqrt{1 + y'^2}dx
\]為最小。

2. [最速下降問題 (Brachistochrone problem)] 
令 $A, B$ 為固定兩點，現在考慮一質點透過重力從 $A$ 滑向 $B$ 點的時間 是與其滑動的路徑有關，故我們的目標是找出一個曲線使得 此質點由 $A$ 滑向 $B$ 的時間最短。

3. [最大面積問題] 
給定固定長度曲線段，試找出此曲線可圍成的最大面積。

事實上，上述所有問題皆可由下列 Functional 表示，再透過變分法求其極值。
\[\int_a^b F (x,y,y')dx\]


那麼我們該如何才能求解 Functional 的極值問題呢?? 事實上我們可以借鏡 數學分析中對於函數的極值問題求解方法。也就是說如果有辦法將 Functional 轉化為 Function 則我們便可以利用傳統數學分析的極值問題來對付它們：

透過 多變數函數分析 近似 Functional

首先回憶我們關心的 Functional 如下
\[
J[y] = \int_a^b F(x,y,y')dx, \;\; y(a) =A, \;\; y(b) = B
\] 我們現在利用 Rieman Integral 的想法來對付 上述 Functional，現在我們將區間 $[a,b]$ 分割成 $n+1$ 等分；亦即令
\[
a=x_0, \; x_1, \; ..., \; x_n, \;x_{n+1} = b
\]那麼我們可以將曲線 $y = y(x)$ 用 多邊線段連線，且對應的多邊線段各端點可寫為
\[
({x_0},\underbrace A_{y({x_0})}),\;({x_1},y({x_1})),\;...,\;({x_n},y({x_n})),\;({x_{n + 1}},\underbrace B_{y({x_{n + 1}})})
\]透過上述分割，我們可以將上述 Functional $J[y]$ 透過下面累加近似
\[J\left( {{y_1},{y_2},...,{y_n}} \right) = \sum\limits_{i = 1}^{n + 1} {F\left( {{x_i},y\left( {{x_i}} \right),\frac{{y\left( {{x_i}} \right) - y\left( {{x_{i - 1}}} \right)}}{h}} \right)} h\]其中  $h = x_i - x_{i-1}$ 且 $y_i := y(x_i)$

最後，我們讓 $n \rightarrow \infty$ 使上述近似還原回 $J[y]$ 。

也就是說，基本變分問題  或者 泛函極值問題 (Functional Extrema Problem) 可以被轉換成 對 $J(y_1,y_2,...,y_n)$ 的 $n$ 變數函數極值問題 再取極限。

Comments
1. 注意到若 $n \rightarrow \infty$ 可看出 $J[y] = J(y_1, y_2, ....)$，亦即泛函可以視為是具有無窮多變數的函數。而我們採用的變分法則可視為是對應此無窮多變數函數 類比於微積分的數學工具。

2. 由於變分問題處理的對象為 Functional，又由前述comment可知 Functional 可視為無窮多變數的函數，故我們處理此問題需要在無窮維的函數空間 (function space)。


延伸閱讀：
[變分法] 泛函極值的必要條件

[最佳控制] Optimizing Multistage Functions - Forward/Backward Dynamic Programming

Life can only be understood going backwards, but it must be lived going forwards. --- Kierkegaard.

考慮一組變數 $w,x,y,z$ 且我們希望最佳化下列的成本函數
\[
f(w,x) + g(x,y) + h(y,z)
\]讀者可已注意到上述的成本函數中有特殊的結構，亦即每一項只有兩個變數。

Backward Dynamic Programming
若我們考慮 $w$ 固定，則上述最佳化問題
\[
\min_{x,y,z} f(w,x) + g(x,y) + h(y,z)
\]可以改寫為 分別最佳化的子問題
\[\mathop {\min }\limits_x \left[ {f(w,x) + \mathop {\min }\limits_y \left[ {g(x,y) + \mathop {\min }\limits_z h(y,z)} \right]} \right]\]故我們可以先解最內部的最佳化問題並得到對應的最佳解 $z^*$ 與 最佳解對應的成本值 $h^*$
\[\left\{ \begin{array}{l}
\mathop {\min }\limits_z h(y,z): = {h^*}\left( y \right)\\
\arg \mathop {\min }\limits_z h(y,z): = {z^*}\left( y \right)
\end{array} \right.\]接著我們求解第二部分的 最佳化的子問題
\[\mathop {\min }\limits_y \left[ {g(x,y) + {h^*}\left( y \right)} \right]\]其對應的解
\[\left\{ \begin{array}{l}
\mathop {\min }\limits_y \left[ {g(x,y) + {h^*}\left( y \right)} \right]: = {g^*}\left( x \right)\\
\arg \mathop {\min }\limits_y \left[ {g(x,y) + {h^*}\left( y \right)} \right]: = {y^*}\left( x \right)
\end{array} \right.\]最後我們求解 最後一部分的 最佳化子問題
\[\mathop {\min }\limits_x \left[ {f(w,x) + {g^*}\left( x \right)} \right]\]其對應的最佳解
\[\left\{ \begin{array}{l}
\mathop {\min }\limits_x \left[ {f(w,x) + {g^*}\left( x \right)} \right]: = {f^*}\left( w \right)\\
\mathop {\arg \min }\limits_x \left[ {f(w,x) + {g^*}\left( x \right)} \right]: = {x^*}\left( w \right)
\end{array} \right.\]

Forward Dynamic Programming
若我們改考慮 $z$ 固定，則上述最佳化問題
\[
\min_{x,y,z} f(w,x) + g(x,y) + h(y,z)
\]可以改寫為 分別最佳化的子問題，下列形式稱為 Forward Dynamic Programming
\[\mathop {\min }\limits_y \left[ {h(y,z) + \mathop {\min }\limits_x \left[ {g\left( {x,y} \right) + \mathop {\min }\limits_w f\left( {w,x} \right)} \right]} \right]\]故可得
\[\underbrace {\mathop {\min }\limits_y \left[ {h(y,z) + \underbrace {\mathop {\min }\limits_x \left[ {g\left( {x,y} \right) + \underbrace {\mathop {\min }\limits_w f\left( {w,x} \right)}_{{f^*}\left( x \right),\begin{array}{*{20}{c}}
{}
\end{array}{w^*}\left( x \right)}} \right]}_{{g^*}\left( {x,y} \right),\begin{array}{*{20}{c}}
{}
\end{array}{x^*}\left( y \right)}} \right]}_{{h^*}(y,z),\begin{array}{*{20}{c}}
{}
\end{array}{y^*}\left( z \right)}\]

[最佳控制] 離散時間 穩態 LQR 控制問題 (1)

延續前篇，這次要介紹的是 Discrete Time Linear Quadratic Regulator in Infinite Horizon 或稱 Steady State LQR。

================
LQR Problem (Infinite Horizon LQR)：
考慮離散狀態方程：
\[
x(k+1) = A x(k) + B u(k)
\]其中 $x(k) \in \mathbb{R}^n, A\in \mathbb{R}^{n \times n}, B \in \mathbb{R}^{n \times m}, u(k) \in \mathbb{R}^{m \times 1}$且 $(A,B)$ controllable。
定義 Performance index：
\[
J(u) = \displaystyle \sum_{k=0}^{\infty} x^T(k+1) Q x(k+1) + u^T(k) R u(k)
\] 其中 $Q, R$ 必須滿足 $Q^T = Q, Q \succ 0$， $R^T = R, R \succ 0$。 (亦即 $Q, R$ 必須為 對稱 + 正定 矩陣)

試求出一組最佳控制力序列 $u^*$ 使得成本函數 $J(u)$ 最小。
================

Comment:
讀者須注意到 Infinite Horizon 的 LQR問題要求計算 Performance index 為無窮級數和，此解必須保證收斂。以下定理告訴我們何時 此 Performance index 收斂

Lemma
考慮離散系統 $x(k+1) = A x(k) + B u(k)$，若 $(A,B)$ 可控制，且選 $Q, R >0$ 為正定矩陣，則上述 infinite horizon LQR 問題保證 閉迴路系統 狀態收斂到 $0$ 且 cost 為有界。

Proof: omitted. (see J. B. Rawlings and D. Q. Mayne, "Model Predictive Control: Theory and Design, p. 24", 2009)


現在我們可以開始求解 Infinite Horizon LQR問題：
Solution
回憶 Steady State Bellman Equation，為了符號簡便起見，我們寫成 functional equation 形式，
\[
I(x) = \displaystyle \min_{u \in \Omega} \{J(x,u) + I(f(x,u)) \}
\] 上式中 $J(x,u)$ 為 Branch cost，亦即 $J(x,u) = x^T Q x + u^T R u$ (並非 $\sum_{k=0}^{\infty} (\cdot)...$)

首先我們猜一組解 $I(x) = x^T P x$ 且矩陣 $P$ 為對稱正定矩陣，亦即滿足 $P^T = P, P \succ 0$。我們之後會找到此 $P$ 應該長甚麼樣子。

將猜測的解代入上述的 Steady State Bellman Equation，故現在我們得到
\[
I(x) = \min_{u \in \Omega} \{J(x,u) + I(f(x,u)) \}
\]注意到 $I(f(x,u) = f(x,u)^T P f(x,u) = (Ax+Bu)^TP(Ax+Bu)$，故我們可得
 \[
\begin{array}{l} I(x) = \mathop {\min }\limits_{u \in \Omega } \{ J(x,u) + I(f(x,u))\} \\ \Rightarrow {x^T}Px = \mathop {\min }\limits_u \left\{ {{x^T}Qx + {u^T}Ru + {{\left( {Ax + Bu} \right)}^T}P\left( {Ax + Bu} \right)} \right\}\\ \Rightarrow {x^T}Px = \mathop {\min }\limits_u \left\{ {{x^T}\left( {Q + {A^T}PA} \right)x + 2{x^T}{A^T}PBu + {u^T}Ru + {u^T}{B^T}PBu} \right\} \end{array}
\]透過一階必要條件 FONC: $ \frac{\partial }{{\partial u}} = 0$ 對上式右邊求解
 \[\begin{array}{l} 2{\left( {{x^T}{A^T}PB} \right)^T} + 2Ru + 2{B^T}PBu = 0\\ \Rightarrow {u^*} = - {\left( {R + {B^T}PB} \right)^{ - 1}}{B^T}PAx \end{array}
\]現在將 $u^*$ 代回 $(*)$  可得 \[\begin{array}{l} {x^T}Px = \mathop {\min }\limits_u \left\{ {{x^T}\left( {Q + {A^T}PA} \right)x + 2{x^T}{A^T}PBu + {u^T}Ru + {u^T}{B^T}PBu} \right\}\\ \Rightarrow {x^T}Px = \left\{ {{x^T}\left\{ {Q + {A^T}PA - {A^T}PB{{\left( {R + {B^T}PB} \right)}^{ - 1}}{B^T}PA} \right\}x} \right\} \end{array}
\]比較左右兩邊可得到 $P$ 必須滿足下式： \[P = Q + {A^T}PA - {A^T}PB{\left( {R + {B^T}PB} \right)^{ - 1}}{B^T}PA\] 此式稱為 Discrete Time Algebraic Ricatti Equation (ARE)，一般而言，可利用 MATLA 指令 dare(A,B,Q, R) 求解 P。

由於 $u^* = - {\left( {R + {B^T}PB} \right)^{ - 1}}{B^T}PAx$，其中除了 $P$ 未定之外，其餘所需要的參數都已知且皆與跌代時間無關，故此無窮時間LQR問題得到的 最佳控制力為 Time invariant。

現在我們總結如下：求解無窮時間的LQR問題只要做兩個步驟即可

STEP 1: 求解一次 Algebraic Ricatti Equation 得到 $P$ (利用 MATLAB: dare.m 或者徒手計算)

\[
P = Q + {A^T}PA - {A^T}PB{\left( {R + {B^T}PB} \right)^{ - 1}}{B^T}PA
\]STEP2 : 將 $P$ 代入 ${u^*} =  - {\left( {R + {B^T}PB} \right)^{ - 1}}{B^T}PAx$

下面我們看個例子：

Example:
考慮一個離散時間線性系統狀態方程：
\[
x_1(k+1) = x_2(k) \\
x_2(k+1) = x_1(k) + u(k)
\]且考慮 Cost function：
\[
J = \sum_{k=0}^{\infty}2x_1^2(k) + 2x_1(k)x_2(k) + x_2^2(k) + 3u^2(k)
\] 且控制力具有如下形式：
\[
u(k) = K_1 x_1(k) + K_2 x_2(k)

\] 試求 $K_1, K_2$ 使 上述 Cost function 最小：

Solution
首先定義  $x(k) := [x_1(k), x_2(k)]^T$ ，則我們有
\[x\left( {k + 1} \right) = \underbrace {\left[ {\begin{array}{*{20}{c}}
0&1\\
1&0
\end{array}} \right]}_A\left[ {\begin{array}{*{20}{c}}
{{x_1}(k)}\\
{{x_2}(k)}
\end{array}} \right] + \underbrace {\left[ {\begin{array}{*{20}{c}}
0\\
1
\end{array}} \right]}_Bu(k)
\] 與 Cost function
\[\begin{array}{l}
J = \sum\limits_{k = 0}^\infty {(2x_1^2(} k) + 2{x_1}(k){x_2}(k) + x_2^2(k) + 3{u^2}(k))\\
\begin{array}{*{20}{c}}
{}
\end{array} = x{\left( k \right)^T}\underbrace {\left[ {\begin{array}{*{20}{c}}
2&1\\
1&1
\end{array}} \right]}_Qx\left( k \right) + \underbrace 3_R{u^2}(k)
\end{array}
\] 那麼現在此問題變成 Steady-state LQR problem，故由前述討論可知我們有 Optimal feedback control 為
\[
u^*(k) = -(R+B^T P B)^{-1} B^T PA \cdot x(k)

\] 其中 $P$ 滿足 $P=P^T, P \succ 0$ 可由 ARE
\[
P= A^TPA - A^T PB (R+ B^TPB)^{-1}B^TPA+Q

\]利用 MATLAB 指令 dare(A,B,Q,R) 解得 $P = \left[ {\begin{array}{*{20}{c}}
{3.7841}&{1.6815}\\
{1.6815}&{4.4022}
\end{array}} \right]$ 現在將 $P$ 帶回 $u^*$中
\[u\left( k \right) = - \left[ {\begin{array}{*{20}{c}}
{0.5947}&{0.2272}
\end{array}} \right]x\left( k \right)

\]i.e., $K_1 = -05947, K_2 =-0.2272$. $\square$

[最佳控制] 離散時間 LQR- Finite Time Horizon

這次要介紹 控制理論中一個重要的結果：
 Discrete Time Linear Quadratic Regulator (LQR) in Finite Time Horizon，
中文翻譯為 離散時間線性二次調節器，我們這邊會針對此問題利用 Dynamic Programming 的方法來逐步求解

================
LQR Problem (Finite Horizon)：
考慮狀態方程：
\[
x(k+1) = A x(k) + B u(k)
\]其中 $x(k) \in \mathbb{R}^n, A\in \mathbb{R}^{n \times n}, B \in \mathbb{R}^{n \times m}, u(k) \in \mathbb{R}^{m \times 1}$
且 考慮 Performance index：
\[
J(u) = \displaystyle \sum_{k=0}^{N-1} x^T(k+1) Q x(k+1) + u^T(k) R u(k)
\] 其中 $Q, R$ 必須滿足 $Q^T = Q, Q \succ 0$， $R^T = R, R \succ 0$。 (亦即 $Q, R$ 必須為 對稱 + 正定 矩陣)

試求出一組最佳控制力序列 $u(N-1), u(N-2),... u(0)$ 使得成本函數 $J(u)$ 最小。
================

Comment:
1. LQR 顧名思義是其具有系統狀態方程為線性 $x(k+1) = Ax(k)+Bu(k)$ 與  Performance Index 中的項都為二次式。
\[
J(u) = \displaystyle \sum_{k=0}^{N-1} x^T(k+1) Q x(k+1) + u^T(k) R u(k)
\]
2. 上述對於矩陣 $Q, R$ 的對稱與正定假設是必須的 (之後需求在求解最佳控制力序列的時候需要求解反矩陣，故需要這些性質。)

3. 上述 LQR in Finite Horizon的問題所求得的最佳控制力序列為 Time Varying 。(此性質會在下面求解的時候再度強調。)

4. 若 Performance index 考慮 $N \rightarrow \infty$，亦即
\[
J(u) = \displaystyle \sum_{k=0}^{\infty} x^T(k+1) Q x(k+1) + u^T(k) R u(k)
\]則我們說此問題為 LQR in Infinite Horizon 或稱 Steady State LQR。這類問題會在之後再做介紹。(此類問題所得到的最佳控制序列將不再是 Time Varying，且只需求解一次 Algebraic Ricatti Equation 即可獲得最佳控制力序列)



Solution of Discrete Time LQR problem in Finite Horizon 
現在我們開始進行求解。

這邊我們使用 Dynamic Programming 方法來求解上述 LQR問題。回憶 Dynamic Programming Equation 的定義：
\[I\left( {x\left( l \right),N - l} \right) = \mathop {\min }\limits_{u\left( l \right) \in {\Omega _l}} \left\{ {J\left( {x\left( l \right),u\left( l \right)} \right) + I\left( {x\left( {l + 1} \right),N - \left( {l + 1} \right)} \right)} \right\}
\] 考慮 Optimal Cost of one-step-to-go ($l=N-1$):
\[\begin{array}{l}
 \Rightarrow I\left( {x\left( {N - 1} \right),1} \right) = \min \left\{ {J\left( {x\left( {N - 1} \right),u\left( {N - 1} \right)} \right) + I\left( {x\left( N \right),0} \right)} \right\}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}
\end{array} = \min \left\{ {{x^T}(N)Qx(N) + {u^T}(N - 1)Ru(N - 1)} \right\}
\end{array}
\]將系統狀態方程 $x(k+1) = A x(k) + B u(k) $ 代入：
\[\begin{array}{l}
 \Rightarrow I\left( {x\left( {N - 1} \right),1} \right) = \\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array}\min \left\{ \begin{array}{l}
{\left[ {Ax\left( {N - 1} \right) + Bu\left( {N - 1} \right)} \right]^T}Q\left[ {Ax\left( {N - 1} \right) + Bu\left( {N - 1} \right)} \right]\\
 + {u^T}(N - 1)Ru(N - 1)
\end{array} \right\}\\
 \Rightarrow I\left( {x\left( {N - 1} \right),1} \right) = \min \left\{ \begin{array}{l}
{x^T}\left( {N - 1} \right){A^T}QAx\left( {N - 1} \right)\\
\begin{array}{*{20}{c}}
{}
\end{array} + 2{x^T}\left( {N - 1} \right){A^T}QBu\left( {N - 1} \right)\\
\begin{array}{*{20}{c}}
{}
\end{array} + {u^T}\left( {N - 1} \right){B^T}QBu\left( {N - 1} \right)\\
\begin{array}{*{20}{c}}
{}
\end{array} + {u^T}(N - 1)Ru(N - 1)
\end{array} \right\}
\end{array}
\] 由一階必要條件 FONC： ${\nabla _{u\left( {N - 1} \right)}}I\left( {x\left( {N - 1} \right),1} \right) = 0$ 可知
\[\begin{array}{l}
2{B^T}{Q^T}Ax\left( {N - 1} \right) + 2{B^T}QBu\left( {N - 1} \right) + 2Ru(N - 1) = 0\\
 \Rightarrow \left[ {R + {B^T}QB} \right]u(N - 1) =  - {B^T}{Q^T}Ax\left( {N - 1} \right)\\
 \Rightarrow u^*(N - 1) =  - \underbrace {{{\left[ {R + {B^T}QB} \right]}^{ - 1}}{B^T}QA}_{K\left( {N - 1} \right)}x\left( {N - 1} \right)
\end{array}
\]接著我們把上述求得的 1-step-to-go 的最佳控制力 $u^*(k)$ 代回 Optimal cost of 1-step-to-go，我們得到 (透過一些代數運算)
\[I\left( {x\left( {N - 1} \right),1} \right) = {x^T}\left( {N - 1} \right)\underbrace {{A^T}\left\{ {Q - QB{{\left[ {R + {B^T}QB} \right]}^{ - 1}}^T{B^T}Q} \right\}A}_{: = M\left( {N - 1} \right)}x\left( {N - 1} \right)
\]在做下一步跌代之前我們先給個 comments

Comments :
1.注意到上式中我們求 ${u^*}(N - 1) =  - {\left[ {R + {B^T}QB} \right]^{ - 1}}{B^T}QAx\left( {N - 1} \right)$ 需要計算反矩陣 ${\left[ {R + {B^T}QB} \right]^{ - 1}}$ ，故需檢驗反矩陣是否存在：不過這個問題可以被證明反矩陣確實存在： (因為 $R$ 為對稱正定矩陣，$B^TQB$ 為對稱半正定矩陣，由FACT: 正定矩陣+ 半正定矩陣 = 正定矩陣，且正定矩陣具有 eigenvalue 全為正，故可推知反矩陣存在。)

2. 上式所求得的最佳控制力 ${u^*}(N - 1) =  - {\left[ {R + {B^T}QB} \right]^{ - 1}}{B^T}QAx\left( {N - 1} \right) =  - K\left( {N - 1} \right)x\left( {N - 1} \right)$ 其中的 $K_{N-1}$ 稱作 Optimal Gain Matrix。此控制力具有回授控制形式 Feedback form。 (類似 $u=-Kx$ 的形式)

3. 最佳控制力儘管具有回授控制 (Feedback form)的形式，但在此問題中為時變得 Gain。此性質在下一步跌代會顯現出來。

4. 當我們有了上述第一步跌代的結果，整個 LQR問題就簡單許多，因為之後的跌代只有 $Q$ 矩陣會改變其餘結果均不變。我們直接看下一步跌代便會發現此性質：

Back to computation :
考慮 Optimal Cost of 2-steps-to-go： ($l=N-2$):
\[\begin{array}{l}
I\left( {x\left( {N - 2} \right),2} \right) = \min \left\{ {J\left( {x\left( {N - 2} \right),u\left( {N - 2} \right)} \right) + I\left( {x\left( {N - 1} \right),1} \right)} \right\}\\
 \Rightarrow I\left( {x\left( {N - 2} \right),2} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \min \left\{ \begin{array}{l}
{x^T}(N - 1)Qx(N - 1) + {u^T}(N - 2)Ru(N - 2)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} + {x^T}\left( {N - 1} \right)M\left( {N - 1} \right)x\left( {N - 1} \right)
\end{array} \right\}\\
 \Rightarrow I\left( {x\left( {N - 2} \right),2} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \min \left\{ {{x^T}(N - 1)\underbrace {\left[ {Q + M\left( {N - 1} \right)} \right]}_{: = \tilde Q}x\left( {N - 1} \right) + {u^T}(N - 2)Ru(N - 2)} \right\}
\end{array}
\]可以發現上式中只有 $Q$ 變成了 $\tilde{Q}$ 其餘參數均固定。故我們可以馬上寫下對應的最佳控制力 $u^*(N-2)$:
\[{u^*}(N - 2) =  - {\left[ {R + {B^T}\tilde QB} \right]^{ - 1}}{B^T}\tilde QAx\left( {N - 2} \right) = K (N-2) x(N-2)
\]其對應的 Optimal cost of 2 steps to go：
\[I\left( {x\left( {N - 2} \right),2} \right) = {x^T}\left( {N - 2} \right)\underbrace {{A^T}\left\{ {\tilde Q - \tilde QB{{\left[ {R + {B^T}\tilde QB} \right]}^{ - 1}}^T{B^T}\tilde Q} \right\}A}_{: = M\left( {N - 2} \right)}x\left( {N - 2} \right)
\]故我們只需要持續重複上述步驟，即可依序解得 $K(N-3), K(N-4),..., K(0)$ 與 $M(N-3), M(N-4),...,M(0)$ 與 $J^*$


Comment:
注意到 $u^*(N-2)$ 的 Gain $K(N-2)$ 與 $u^*(N-1)$ 的 Gain $K(N-1)$ 並不相同，此說明了之前所敘的時變特性(Time Varying Property)

[最佳控制] Finite Horizon LQR with Penalty on Final State

考慮 LQR 問題如下

成本函數
\[{V_N}({x_0},u): = \sum\limits_{i = 1}^{N - 1} {\underbrace {l(x\left( i \right),u\left( i \right))}_{{\rm{stage}}\begin{array}{*{20}{c}}
{}
\end{array}{\rm{cost}}}}  + \underbrace {{l_N}({x_N})}_{{\rm{terminal}}\begin{array}{*{20}{c}}
{}
\end{array}{\rm{cost}}}\]其中 \[\left\{ \begin{array}{l}
u: = \{ u(0),u(1),...,u(N - 1)\} \\
l(x,u): = \frac{1}{2}({x^T}Qx + {u^T}Ru)\\
{l_N}(x): = \frac{1}{2}{x^T}{P_f}x
\end{array} \right.\]另外假設 $Q \ge 0$  positive semi-definite 且 $R >0$ positive definite。

我們的目標：找到 $u^*=\{u^*(0),u^*(1),...,u^*(N-1)\}$ 使得
\[\begin{array}{l}
\mathop {\min }\limits_u {V_N}\left( {x\left( 0 \right),u} \right)\\
s.t.\;\;{x^ + } = Ax + Bu
\end{array}\]

現在利用 Backward Dynamic Programming，我們從最後的狀態 $x(N)$ 逐步回推最佳解，亦即觀察
\[{V_N}({x_0},u): = l({x_0},{u_0}) + l({x_1},{u_1}) + ... + \underbrace {l({x_{N - 1}},{u_{N - 1}}) + {l_N}({x_N})}_{{u_{N - 1}}\begin{array}{*{20}{c}}
{}
\end{array}{\rm{affects}}\begin{array}{*{20}{c}}
{}
\end{array}{\rm{only}}\begin{array}{*{20}{c}}
{}
\end{array}{\rm{there!}}}\]由 Optimality Principle 逐步求解的最佳解 必為整體最佳的一環，故此我們先最佳化 下列子問題
\[ \begin{array}{l}
\mathop {\min }\limits_{{u_{N - 1}}} l({x_{N - 1}},{u_{N - 1}}) + {l_N}({x_N})\\
s.t.\begin{array}{*{20}{c}}
{}
\end{array}{x_N} = A{x_{N - 1}} + B{u_{N - 1}}
\end{array}\]現在將拘束條件帶入，我們可得到以下無拘束最佳化問題：
\[\small \begin{array}{l}
\left\{ \begin{array}{l}
\mathop {\min }\limits_{{u_{N - 1}}} l({x_{N - 1}},{u_{N - 1}}) + {l_N}({x_N})\\
s.t.\begin{array}{*{20}{c}}
{}
\end{array}{x_N} = A{x_{N - 1}} + B{u_{N - 1}}
\end{array} \right.\\
 \Rightarrow \mathop {\min }\limits_{{u_{N - 1}}} \frac{1}{2}\left( {x_{N - 1}^TQx_{N - 1}^{} + u_{N - 1}^TRu_{N - 1}^{}} \right) + \frac{1}{2}{\left( {A{x_{N - 1}} + B{u_{N - 1}}} \right)^T}{P_f}\left( {A{x_{N - 1}} + B{u_{N - 1}}} \right)\\
 \Rightarrow \mathop {\min }\limits_{{u_{N - 1}}} \frac{1}{2}\left[ {x_{N - 1}^TQx_{N - 1}^{} + \underbrace {u_{N - 1}^TRu_{N - 1}^{}}_{{V_1}} + \underbrace {{{\left( {A{x_{N - 1}} + B{u_{N - 1}}} \right)}^T}{P_f}\left( {A{x_{N - 1}} + B{u_{N - 1}}} \right)}_{{V_2}}} \right] \ \ \ \ \ (*)
\end{array}\]現在我們可利用下列結果
------------
FACT: 若
\[\left\{ \begin{array}{l}
{V_1} = \frac{1}{2}{\left( {x - a} \right)^T}\Phi \left( {x - a} \right)\\
{V_2} = \frac{1}{2}{\left( {\Theta x - b} \right)^T}\Gamma \left( {\Theta x - b} \right)
\end{array} \right.\]且 $\Phi $ positive definite 與 $\Theta $ 為 positive semi-definite。則其和可表為\[V = {V_1} + {V_2} = \frac{1}{2}{\left( {x - v} \right)^T}\Omega \left( {x - v} \right) + d\]其中
\[\left\{ \begin{array}{l}
\Omega  = \Phi  + {\Theta ^T}\Gamma \Theta \\
v = {\left( {\Phi  + {\Theta ^T}\Gamma \Theta } \right)^{ - 1}}\left( {\Phi a + {\Theta ^T}\Gamma b} \right)\\
d = V\left( v \right)
\end{array} \right.\]
------------

觀察式 $(*)$，我們可利用上述 FACT 比較係數
\[\Phi  = R;\Theta  = B;\Gamma  = {P_f};v = u_{N - 1}^*;b =  - A{x_{N - 1}};a = 0
\]故在此階段的最佳解 $u_{N-1}^*$ 為

\[u_{N - 1}^* = \underbrace { - {{\left( {R + {B^T}{P_f}B} \right)}^{ - 1}}{B^T}{P_f}A}_{{K_{N - 1}}}{x_{N - 1}} = {K_{N - 1}}{x_{N - 1}}
\]且在此階段的 optimal cost (又稱 optimal cost to go) 為
\[V_{N - 1}^*\left( {{x_{N - 1}}} \right) = d + \frac{1}{2}x_{N - 1}^TQx_{N - 1}^{}
\]其中 $d$ 為前述 FACT 中的常數項，我們計算如下
\[\small \begin{array}{*{20}{l}}
{d = V\left( v \right) = {V_1}\left( v \right) + {V_2}\left( v \right)}\\
{\begin{array}{*{20}{c}}
{}&{}
\end{array} = \frac{1}{2}{{\left( {x - a} \right)}^T}\Phi \left( {x - a} \right) + \frac{1}{2}{{\left( {\Theta x - b} \right)}^T}\Gamma \left( {\Theta x - b} \right)}\\
{\begin{array}{*{20}{c}}
{}&{}
\end{array} = \frac{1}{2}{{\left( {{K_{N - 1}}{x_{N - 1}}} \right)}^T}R\left( {{K_{N - 1}}{x_{N - 1}}} \right) + \frac{1}{2}{{\left( {\left( {A + B{K_{N - 1}}} \right){x_{N - 1}}} \right)}^T}{P_f}\left( {\left( {A + B{K_{N - 1}}} \right){x_{N - 1}}} \right)}\\
{\begin{array}{*{20}{c}}
{}&{}
\end{array} = \frac{1}{2}x_{N - 1}^T\left[ {K_{N - 1}^TR{K_{N - 1}} + {{\left( {A + B{K_{N - 1}}} \right)}^T}{P_f}\left( {A + B{K_{N - 1}}} \right)} \right]{x_{N - 1}}}
\end{array}
\]故 optimal cost to go
\[\small \begin{array}{l}
V_{N - 1}^*\left( {{x_{N - 1}}} \right) = d + \frac{1}{2}x_{N - 1}^TQx_{N - 1}^{}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \frac{1}{2}x_{N - 1}^T\left[ {K_{N - 1}^TR{K_{N - 1}} + {{\left( {A + B{K_{N - 1}}} \right)}^T}{P_f}\left( {A + B{K_{N - 1}}} \right)} \right]{x_{N - 1}} + \frac{1}{2}x_{N - 1}^TQx_{N - 1}^{}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \frac{1}{2}x_{N - 1}^T\left[ {Q + K_{N - 1}^TR{K_{N - 1}} + {{\left( {A + B{K_{N - 1}}} \right)}^T}{P_f}\left( {A + B{K_{N - 1}}} \right)} \right]{x_{N - 1}}
\end{array}\]如果我們把最佳解 $u_{N - 1}^* = {K_{N - 1}}{x_{N - 1}} =  - {\left( {R + {B^T}{P_f}B} \right)^{ - 1}}{B^T}{P_f}A{x_{N - 1}}$ 帶入 optimal cost to go，可得
\[\small \begin{array}{l}
V_{N - 1}^*\left( {{x_{N - 1}}} \right) = \frac{1}{2}x_{N - 1}^T\left[ {Q + K_{N - 1}^TR{K_{N - 1}} + {{\left( {A + B{K_{N - 1}}} \right)}^T}{P_f}\left( {A + B{K_{N - 1}}} \right)} \right]{x_{N - 1}}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \frac{1}{2}x_{N - 1}^T\left[ {Q + {A^T}{P_f}A + K_{N - 1}^T\left( {R + {B^T}{P_f}B} \right){K_{N - 1}} + 2K_{N - 1}^T{B^T}{P_f}A} \right]{x_{N - 1}}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \frac{1}{2}x_{N - 1}^T\underbrace {\left[ {Q + {A^T}{P_f}A - {A^T}{P_f}B{{\left( {R + {B^T}{P_f}B} \right)}^{ - 1}}{B^T}{P_f}A} \right]}_{{M_{N - 1}}}{x_{N - 1}}
\end{array}\]但注意到我們並不曉得 $x_{N-1}$ 故我們需要繼續回頭解 $x_{N-2}$ 亦即 接著求解
\[\begin{array}{l}
\left\{ \begin{array}{l}
\mathop {\min }\limits_{{u_{N - 2}}} l\left( {{x_{N - 2}},{u_{N - 2}}} \right) + V_{N - 1}^*\left( {{x_{N - 1}}} \right)\\
s.t.\begin{array}{*{20}{c}}
{}
\end{array}{x_{N - 1}} = A{x_{N - 2}} + B{u_{N - 2}}
\end{array} \right.\\
 \Rightarrow \left\{ \begin{array}{l}
\mathop {\min }\limits_{{u_{N - 2}}} l\left( {{x_{N - 2}},{u_{N - 2}}} \right) + \frac{1}{2}x_{N - 1}^T{M_{N - 1}}{x_{N - 1}}\\
s.t.\begin{array}{*{20}{c}}
{}
\end{array}{x_{N - 1}} = A{x_{N - 2}} + B{u_{N - 2}}
\end{array} \right.
\end{array}\]但注意到先前我們已經解了
\[\left\{ \begin{array}{l}
\mathop {\min }\limits_{{u_{N - 1}}} l\left( {{x_{N - 1}},{u_{N - 1}}} \right) + \frac{1}{2}x_N^T{P_f}{x_N}\\
s.t.\begin{array}{*{20}{c}}
{}
\end{array}{x_N} = A{x_{N - 1}} + B{u_{N - 1}}
\end{array} \right.\]故讀者可直接比對上述兩者，即可得知只要將 前述我們所推出的解其中的 $P_f$ 矩陣換成 $M_{n-1}$ 即可馬上得到 $N-2$ stage的最佳解 與 optimal cost to go。亦即
\[\left\{ \begin{array}{l}
u_{N - 2}^* = {K_{N - 2}}{x_{N - 2}} =  - {\left( {R + {B^T}{M_{N - 1}}B} \right)^{ - 1}}{B^T}{M_{N - 1}}A{x_{N - 2}}\\
V_{N - 2}^*\left( {{x_{N - 2}}} \right) = \frac{1}{2}x_{N - 2}^T{M_{N - 2}}{x_{N - 2}}\\
{M_{N - 2}} = Q + {A^T}{M_{N - 1}}A - {A^T}{M_{N - 1}}B{\left( {R + {B^T}{M_{N - 1}}B} \right)^{ - 1}}{B^T}{M_{N - 1}}A
\end{array} \right.\]

Summary
透過跌代求解 Riccati equation 進而得到一組控制力序列$u(0),u(1),...u(N-1)$；亦即透過下列跌代式：
\[\begin{array}{l}
for\begin{array}{*{20}{c}}
{}
\end{array}k = N - 1,N - 2,...,0\\
\left\{ \begin{array}{l}
u_{}^*\left( k \right) = K\left( k \right)x\left( k \right),\begin{array}{*{20}{c}}
{}&{}
\end{array}\\
K\left( k \right) =  - {\left( {{B^T}M\left( {k + 1} \right)B + R} \right)^{ - 1}}{B^T}M\left( {k + 1} \right)A,
\end{array} \right.\\
\\
for\begin{array}{*{20}{c}}
{}
\end{array}k = N,N - 1,...,0\\
\left\{ \begin{array}{l}
M\left( {k - 1} \right): = Q + {A^T}M\left( k \right)A - {A^T}M\left( k \right)B{\left( {{B^T}M\left( k \right)B + R} \right)^{ - 1}}{B^T}M\left( k \right)A\\
M\left( N \right): = {P_f}\\
V_{}^*\left( k \right) = \frac{1}{2}{x^T}\left( k \right)M\left( {k + 1} \right)x\left( k \right)
\end{array} \right.
\end{array}\]

但注意到最佳解並不一定保證系統穩定 (Kalman (1960b, p.113) 指出系統的 optimality 並不保證 stability，亦即 系統採用 optimal control law 並不保證系統閉迴路穩定。)，現在我們將透過以下例子顯示 儘管 $Q >0, R>0$ 且 $N \ge 1$ 所求得的最佳控制力並不保證系統閉迴路穩定。

Example: Finite Horizon LQ control
考慮離散系統 $x(k+1) = Ax(k) + B u(k)$ 且 $y(k) = Cx(k)$；其中
\[A = \left[ {\begin{array}{*{20}{c}}
{4/3}&{ - 2/3}\\
1&0
\end{array}} \right];B = \left[ {\begin{array}{*{20}{c}}
1\\
0
\end{array}} \right];C = \left[ {\begin{array}{*{20}{c}}
{ - 2/3}&1
\end{array}} \right]\]現在考慮  $N=5$ 與 $N=7$ ，試設計 有限維度 LQ 控制器時的最佳解。
Solution
上述系統對應的 轉移函數如下
\[G\left( z \right) = \frac{{ - 2/3z + 1}}{{{z^2} - 4/3z + 2/3}}\]注意到此系統的 開迴路 零點 (zero) 落在 $z = 3/2 = 1.5$；亦即具有 不穩定 零點 ( 因為離散系統穩定範圍落在單位圓 $|z| \le 1$ 之內。) 現在我們利用 LQ 控制器，令參數矩陣 $R:=0.001$ 且 $Q := C^TC = P_f$ 且 $Q \ge 0$ 為半正定矩陣；此時若我們額外再加入一點小擾動
\[\begin{array}{l}
Q: = {C^T}C + 0.001I\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left[ {\begin{array}{*{20}{c}}
{ - 2/3}\\
1
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{ - 2/3}&1
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
{0.001}&0\\
0&{0.001}
\end{array}} \right]\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left[ {\begin{array}{*{20}{c}}
{4/9}&{ - 2/3}\\
{ - 2/3}&1
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
{0.001}&0\\
0&{0.001}
\end{array}} \right]\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left[ {\begin{array}{*{20}{c}}
{4/9 + 0.001}&{ - 2/3}\\
{ - 2/3}&{1.001}
\end{array}} \right]
\end{array}\]接著我們逐步跌代求解最佳控制力序列 $u(4),u(3),u(2),u(1),u(0)$ 如下：
$k = N = 5$，我們有
\[\begin{array}{l}
 \Rightarrow \left\{ \begin{array}{l}
u_{}^*\left( 4 \right) = K\left( 4 \right)x\left( 4 \right),\begin{array}{*{20}{c}}
{}&{}
\end{array}\\
K\left( 4 \right) = \left[ {\begin{array}{*{20}{c}}
{0.1629}&{0.6652}
\end{array}} \right]\\
M\left( 5 \right) = Q\\
V_{}^*\left( 5 \right) = \frac{1}{2}{x^T}\left( 5 \right)Qx\left( 5 \right)
\end{array} \right. \Rightarrow \left\{ \begin{array}{l}
u_{}^*\left( 3 \right) = K\left( 3 \right)x\left( 3 \right),\begin{array}{*{20}{c}}
{}&{}
\end{array}\\
K\left( 3 \right) = \left[ {\begin{array}{*{20}{c}}
{0.1518}&{0.6652}
\end{array}} \right],\\
M\left( 4 \right): = \left[ {\begin{array}{*{20}{c}}
{0.4489}&{ - 0.666}\\
{ - 0.666}&{1.0014}
\end{array}} \right]\\
V_{}^*\left( 4 \right) = \frac{1}{2}{x^T}\left( 4 \right)M\left( 4 \right)x\left( 4 \right)
\end{array} \right.\\
 \Rightarrow \left\{ \begin{array}{l}
u_{}^*\left( 2 \right) = K\left( 2 \right)x\left( 2 \right),\begin{array}{*{20}{c}}
{}&{}
\end{array}\\
K\left( 2 \right) = \left[ {\begin{array}{*{20}{c}}
{0.1257}&{0.6652}
\end{array}} \right]\\
M\left( 3 \right): = \left[ {\begin{array}{*{20}{c}}
{0.4568}&{ - 0.666}\\
{ - 0.666}&{1.0014}
\end{array}} \right]\\
V_{}^*\left( 3 \right) = \frac{1}{2}{x^T}\left( 3 \right)M\left( 3 \right)x\left( 3 \right),
\end{array} \right. \Rightarrow \left\{ \begin{array}{l}
u_{}^*\left( 1 \right) = K\left( 1 \right)x\left( 1 \right),\begin{array}{*{20}{c}}
{}&{}
\end{array}\\
K\left( 1 \right) = \left[ {\begin{array}{*{20}{c}}
{0.0724}&{0.6653}
\end{array}} \right],\\
M\left( 2 \right): = \left[ {\begin{array}{*{20}{c}}
{0.4741}&{ - 0.666}\\
{ - 0.666}&{1.0014}
\end{array}} \right]\\
V_{}^*\left( 2 \right) = \frac{1}{2}{x^T}\left( 2 \right)M\left( 2 \right)x\left( 2 \right)
\end{array} \right.\\
 \Rightarrow \left\{ \begin{array}{l}
u_{}^*\left( 0 \right) = K\left( 0 \right)x\left( 0 \right),\begin{array}{*{20}{c}}
{}&{}
\end{array}\\
K\left( 0 \right) = \left[ {\begin{array}{*{20}{c}}
{ - 0.0256}&{0.6653}
\end{array}} \right],\\
M\left( 1 \right): = \left[ {\begin{array}{*{20}{c}}
{0.5098}&{ - 0.666}\\
{ - 0.666}&{1.0014}
\end{array}} \right]\\
V_{}^*\left( 1 \right) = \frac{1}{2}{x^T}\left( 1 \right)M\left( 1 \right)x\left( 1 \right)
\end{array} \right.
\end{array}\]注意到若我們計算閉迴路系統的極點可得
$$eig(A+ BK_{N=5}(0)) = \{1.307, 0.001\}$$ 亦即 使用此最佳控制力會導致閉迴路系統不穩定 $(z = 1.307 \ge 1)$ 但讀者可自行試驗若改考慮 $N=7$時候，閉迴路極點會變成
\[
eig(A+ BK_{N=7}(0)) = \{0.989, 0.001\}
\]亦即系統被穩定化。若我們考慮 $N=\infty$ 時可以得到
\[
eig(A+ BK_{N=\infty}(0)) = \{0.664, 0.001\}
\]此結果稱作 infinite horizon control law。會在之後再行介紹。

[線性代數] 矩陣二次式的等價運算

假設 $A$ 為 對稱 正定矩陣  (亦即 $A^T = A$ 且 $A$ 的 eigenvalue 全為正值)，現在考慮一個矩陣二次函數：
\[
V(x) = x^T A x + c^T x + d
\]上述矩陣二次項為 $x^TAx$ 且 線性項為 $c^T x$ 常數項為 $d$。

注意到上式可改寫為 $ V(x) = (x-v)^T H (x-v) +d$。WHY? 因為改寫成此形式之後，最小值一目了然，亦即 $x=v$ 可得最小值。


現在若考慮兩組矩陣二次式
\[\left\{ \begin{array}{l}
{V_1}(x) = \frac{1}{2}{(x - a)^T}A(x - a)\\
{V_2}(x) = \frac{1}{2}{(x - b)^T}B(x - b)
\end{array} \right.\]且假設 $A >0$ 為 正定矩陣，$B$為半正定矩陣。

===============
FACT:  $V_1, V_2$ 皆為矩陣二次式，其和亦為矩陣二次式；亦即
$$V(x) = \frac{1}{2} (x-v)^T H (x-v) + d = V_1(x) + V_2(x)
$$===============
故現在問題變成如何找出 $d, H, v $ 用 $A,B,a,b$表示?

===============
FACT: 考慮兩組矩陣二次式
\[\left\{ \begin{array}{l}
{V_1}(x) = \frac{1}{2}{(x - a)^T}A(x - a)\\
{V_2}(x) = \frac{1}{2}{(x - b)^T}B(x - b)
\end{array} \right.\]若 $A^T = A$ 且 $B^T = B$，則 $$V(x) = \frac{1}{2} (x-v)^T H (x-v) + d = V_1(x) + V_2(x)$$ 且
\[\left\{ \begin{array}{l}
H = A + B\\
v = {\left( {A + B} \right)^{ - 1}}\left( {Aa + Bb} \right)\\
d = {V_1}\left( v \right) + {V_2}\left( v \right)
\end{array} \right.\]===============

Proof:
注意到
\[V(x) = \frac{1}{2}\left( {{x^T}Hx - 2{x^T}Hv + {v^T}Hv} \right) + d\]對 $V(x)$ 取 一階導數 與 二階導數，可得
\[\left\{ \begin{array}{l}
\frac{d}{{dx}}V(x) = \frac{1}{2}\left( {Hx + {H^T}x - 2Hv} \right) = H\left( {x - v} \right)\\
\frac{{{d^2}}}{{d{x^2}}}V(x) = \frac{d}{{dx}}\left( {\frac{d}{{dx}}V(x)} \right) = \frac{d}{{dx}}\left( {H\left( {x - v} \right)} \right) = {H^T} = H
\end{array} \right.
\]又因為
\[\left\{ \begin{array}{l}
\frac{d}{{dx}}V(x) = \frac{d}{{dx}}{V_1}(x) + \frac{d}{{dx}}{V_2}(x)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \frac{1}{2}(Ax + {A^T}x - 2Aa) + \frac{1}{2}(Bx + {B^T}x - 2Bb)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = A(x - a) + B(x - b)\\
\frac{{{d^2}}}{{d{x^2}}}V(x) = \frac{d}{{dx}}\left( {A(x - a) + B(x - b)} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = A + B
\end{array} \right.\]現在比較手邊結果可得
\[\begin{array}{l}
\left\{ \begin{array}{l}
\frac{d}{{dx}}V(x) = H\left( {x - v} \right) = A(x - a) + B(x - b)\\
\frac{{{d^2}}}{{d{x^2}}}V(x) = H = A + B
\end{array} \right.\\
 \Rightarrow \left\{ \begin{array}{l}
v = {\left( {A + B} \right)^{ - 1}}\left( {Aa + Bb} \right)\\
H = A + B
\end{array} \right.
\end{array}\]上述 $v$ 的求解 需要 $(A+B)^{-1}$ 但因為我們假設 $A$ 為正定 且 $B$ 為半正定，故 $(A+B)$ 為正定矩陣，反矩陣存在。

接著我們計算常數項 $d$：注意到
\[V\left( v \right) = d = {V_1}\left( v \right) + {V_2}\left( v \right)\]


現在我們進一步推廣上述結果：
================
FACT: 
考慮兩組矩陣二次式
\[\left\{ {\begin{array}{*{20}{l}}
{{V_1}(x) = \frac{1}{2}{{(x - a)}^T}A(x - a)}\\
{{V_2}(x) = \frac{1}{2}{{(Cx - b)}^T}B(Cx - b)}
\end{array}} \right.\]
若 $A^T = A$ 且 $B^T = B$，則 $$V(x) = \frac{1}{2} (x-v)^T H (x-v) + d = V_1(x) + V_2(x)$$  且
\[{\left\{ \begin{array}{l}
\begin{array}{*{20}{l}}
{v = {{\left( {A + {C^T}BC} \right)}^{ - 1}}\left( {Aa + {C^T}Bb} \right)}\\
{H = A + {C^T}BC}
\end{array}\\
d = {V_1}\left( v \right) + {V_2}\left( v \right)
\end{array} \right.}\]================

Proof:
證明同前述 FACT，注意到
\[V(x) = \frac{1}{2}\left( {{x^T}Hx - 2{x^T}Hv + {v^T}Hv} \right) + d\]現在分別對 $V(x)$ 取 一階導數 與 二階導數，可得
\[\left\{ \begin{array}{l}
\frac{d}{{dx}}V(x) = \frac{1}{2}\left( {Hx + {H^T}x - 2Hv} \right) = H\left( {x - v} \right)\\
\frac{{{d^2}}}{{d{x^2}}}V(x) = \frac{d}{{dx}}\left( {\frac{d}{{dx}}V(x)} \right) = \frac{d}{{dx}}\left( {H\left( {x - v} \right)} \right) = {H^T} = H
\end{array} \right.
\]又因為
\[\left\{ {\begin{array}{*{20}{l}}
{\frac{d}{{dx}}V(x) = \frac{d}{{dx}}{V_1}(x) + \frac{d}{{dx}}{V_2}(x)}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \frac{1}{2}(Ax + {A^T}x - 2Aa) + \frac{1}{2}({C^T}BCx + {C^T}BCx - 2{C^T}Bb)}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = A(x - a) + {C^T}B(Cx - b)}\\
{\frac{{{d^2}}}{{d{x^2}}}V(x) = \frac{d}{{dx}}\left( {A(x - a) + {C^T}B(Cx - b)} \right)}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = A + {C^T}BC}
\end{array}} \right.\]現在比較手邊結果可得
\[\begin{array}{*{20}{l}}
{\left\{ {\begin{array}{*{20}{l}}
{\frac{d}{{dx}}V(x) = H\left( {x - v} \right) = A(x - a) + {C^T}B(Cx - b)}\\
{\frac{{{d^2}}}{{d{x^2}}}V(x) = H = A + {C^T}BC}
\end{array}} \right.}\\
{ \Rightarrow \left\{ {\begin{array}{*{20}{l}}
{v = {{\left( {A + {C^T}BC} \right)}^{ - 1}}\left( {Aa + {C^T}Bb} \right)}\\
{H = A + {C^T}BC}
\end{array}} \right.}
\end{array}\]上述 $v$ 的求解 需要 $(A+C^TBC)^{-1}$ 但因為我們假設 $A$ 為正定 且 $C^TBC$ 永遠為半正定，故 $(A+C^TBC)$ 為正定矩陣，反矩陣存在。

接著我們計算常數項 $d$：注意到
\[V\left( v \right) = d = {V_1}\left( v \right) + {V_2}\left( v \right)\]

[動態規劃] 淺談 離散時間動態規劃 (1) - Bellman equation in Infinite Horizon

這次要介紹 無窮時間的 Bellman Equation，亦即我們的 Cost function 為 branch cost 加到無窮大的情況：
\[
J(u) := \displaystyle \sum_{k=0}^{\infty} J(x(k), u(k)) + \Phi(x(N))
\]
那麼現在問題變成 儘管我們手上有 有限時間的 Bellman Equation (請參考前篇文章)，但對於此類無窮時間的問題該如何處理??

亦即如果今天我們 cost function 的 $N \rightarrow \infty$ 該怎麼處理? 我們稱這一類問題叫做 Steady State Dynamic Programming 或稱 Dynamic Programming in Infinite Horizon。

無窮時間的動態規劃問題 (Dynamic Programming Problem in Infinite Horizon)：

考慮 Performance index (cost function)
\[
J(u) := \displaystyle \sum_{k=0}^{\infty} J(x(k), u(k)) + \Phi(x(N))
\]狀態方程( state equation)
\[
x(k+1) = f(x(k), u(k)), \ x(0) \ \text{is given} \\
\]其中 $x(k)$ 為系統在第 $k$ 時刻的 狀態，$f:\mathbb{R}^n \times \mathbb{R}^m \rightarrow \mathbb{R}^n$，
與控制力拘束條件
\[
u(k) \in \Omega
\]其中$\Omega$ 為拘束條件

此時對應的 Bellman Equation (or Dynamic Programming Equation) 寫為如下的 functional form
\[
I(x) = \min_{u \in \Omega} \{ J(x,u) + I(f(x,u)) \}
\]亦即與 $N$ 無關
上式稱為 Steady State Bellman Equation 或者 Bellman Equation in Infinite Horizon。


現在我們看個例子：

Example
考慮系統狀態方程表示如下：
\[x(k+1) = x(k) - u(k)
\]且 cost function 為
\[
J(u) = \displaystyle \sum_{k=0}^{\infty} (x(k+1)-u(k))^2 + x^2(k+1)
\]試求 Optimal $u^*$ 與其對應的 optimal cost to go

Solution
首先我們寫下 Steady State Bellman Equation:
\[
I(x) = \min_{u \in \Omega} \{J(x,u) + I(f(x,u)) \}
\]現在我們要求解上式，故我們首先猜一個解 (事實上此問題在之前文章中我們解過有限時間的 Bellman equation，當時我們解出 Optimal cost to go 具有 $I(x) = \alpha x^2$ 的形式，有興趣的讀者可前往前篇文章檢驗)，故我們現在很合理的可以猜這 $I(x) = \alpha x^2$ 其中 $\alpha$ 為代定係數 $\alpha \geq 0$。
故我們將此解代入 Steady State Bellman Equation，可得
\[\begin{array}{l}
I(x) = {\min _{u \in \Omega }}\{ J(x,u) + I(f(x,u))\} \\
 \Rightarrow \alpha {x^2} = {\min _{u \in \Omega }}\{ {(x - u)^2} + {x^2} + \alpha {\left( {x - u} \right)^2}\} \\
 \Rightarrow \alpha {x^2} = {\min _{u \in \Omega }}\{ {x^2} + \left( {1 + \alpha } \right){\left( {x - u} \right)^2}\}
\end{array}
\]由 FONC: $\frac{\partial }{{\partial u}} = 0$，我們可求解上式右方最佳控制力 $u^*$ ：
\[
u^* = x
\]故代回上式我們得到
\[\begin{array}{l}
\alpha {x^2} = {x^2}\\
 \Rightarrow \alpha  = 1
\end{array}
\] 故如果我們選 $\alpha =1$即可解得 Steady-State Bellman Equation。 $\square$

上述問題可以進一步拓展到 $\mathbb{R}^n$ 空間，這一類問題在控制理論中稱為無窮時間的線性二次調節器 (Linear Quadratic Regulator, LQR) 問題，有興趣的讀者請參考：
[最佳控制] 離散時間 穩態線性二次調節器 Discrete Time Linear Quadratic Regulator in Infinite Horizon via Dynamic Programming