[系統理論] Picard Iteration

考慮狀態空間
\[{\bf{\dot x}}\left( t \right) = f\left( {t,{\bf{x}}\left( t \right)} \right),{\bf{x}}\left( 0 \right) = {{\bf{x}}^0}\]則 Picard Iteration 給定
1. Initial guess
\[
x(t) = x^0
\]2.. Update Step
\[{{\bf{x}}^{k + 1}}\left( t \right) = {{\bf{x}}^0} + \int_0^t {f\left( {\tau ,{{\bf{x}}^k}\left( \tau  \right)} \right)d\tau } \]
注意到上述迭代式為 sequence of $\{{\bf x}^k\}$


Example 1:
考慮下列非線性系統
\[{\bf{\dot x}}\left( t \right) = \left[ \begin{array}{l}
{{\dot x}_1}\\
{{\dot x}_2}
\end{array} \right] = \left[ \begin{array}{l}
\cos {x_1}\\
t{x_1} + {e^{ - t}}{x_2}
\end{array} \right]
\]且 $x_1(0) = 2$ 與 $x_2(0) = -1$。
試透過 Picard iteration 求取 ${\bf x }^1(t)$

Solution:
\[\begin{array}{l}
{{\bf{x}}^{k + 1}}\left( t \right) = {{\bf{x}}^0} + \int_0^t {f\left( {\tau ,{{\bf{x}}^k}\left( \tau  \right)} \right)d\tau } \\
{{\bf{x}}^1}(t) = {{\bf{x}}^0}(t) + \int_0^t {\left[ \begin{array}{l}
\cos {x_1}\left( \tau  \right)\\
\tau {x_1} + {e^{ - \tau }}{x_2}\left( \tau  \right)
\end{array} \right]} d\tau \\
 \Rightarrow {{\bf{x}}^1}(t) = \left[ \begin{array}{l}
2\\
 - 1
\end{array} \right] + \int_0^t {\left[ \begin{array}{l}
\cos 2\\
\tau 2 + {e^{ - \tau }}\left( { - 1} \right)
\end{array} \right]} d\tau \\
 \Rightarrow {{\bf{x}}^1}(t) = \left[ \begin{array}{l}
2\\
 - 1
\end{array} \right] + \left[ \begin{array}{l}
t\cos 2\\
{t^2} + \left( {{e^{ - t}} - 1} \right)
\end{array} \right] = \left[ \begin{array}{l}
2 + t\cos 2\\
{t^2} + {e^{ - t}} - 2
\end{array} \right]

\end{array}\]

Exercise: 讀者可自行嘗試透過 Picard Iteration 求解 $x^2(t)$ 與 $x^3(t)$


Example 2: Saturation 
考慮非線性系統 $\dot{x} = f(x)$ 且 $x(0) =1$ 其中
\[f\left( x \right) = \left\{ \begin{array}{l}
1,\begin{array}{*{20}{c}}
{}&{}
\end{array}x > 2\\
\frac{1}{2}x,\begin{array}{*{20}{c}}
{}&{}
\end{array}\left| x \right| \le 2\\
 - 1,\begin{array}{*{20}{c}}
{}&{}
\end{array}x <  - 2
\end{array} \right.\]試透過 Picard iteration 求取 $ x^1(t)$ 與 $x^2(t)$
Solution:
首先計算 $x^1(t)$:
\[\begin{array}{l}
{x^1}\left( t \right) = {x^0}\left( t \right) + \int_0^t {f\left( {{x^0}\left( \tau  \right)} \right)d\tau } \\
 \Rightarrow {x^1}\left( t \right) = 1 + \int_0^t {f\left( 1 \right)d\tau } \\
 \Rightarrow {x^1}\left( t \right) = 1 + \frac{1}{2}t
\end{array}\]
接著我們計算 $x^2(t)$:
\[\begin{array}{l}
{x^2}\left( t \right) = {x^0}\left( t \right) + \int_0^t {f\left( {{x^1}\left( \tau  \right)} \right)d\tau } \\
 \Rightarrow {x^2}\left( t \right) = 1 + \left\{ \begin{array}{l}
\int_0^t {\frac{1}{2}\left( {1 + \frac{1}{2}\tau } \right)d\tau ,\begin{array}{*{20}{c}}
{}&{}
\end{array}0 \le t \le 2} \\
\frac{1}{2}\int_0^t {1d\tau } ,\begin{array}{*{20}{c}}
{}&{}
\end{array}t > 2
\end{array} \right.\\
 \Rightarrow {x^2}\left( t \right) = 1 + \left\{ \begin{array}{l}
\frac{1}{2}\left( {t + \frac{1}{4}{t^2}} \right),\begin{array}{*{20}{c}}
{}&{}
\end{array}0 \le t \le 2\\
\frac{1}{2}t,\begin{array}{*{20}{c}}
{}&{}
\end{array}t > 2
\end{array} \right.\\
 \Rightarrow {x^2}\left( t \right) = \left\{ \begin{array}{l}
1 + \frac{t}{2} + \frac{1}{8}{t^2},\begin{array}{*{20}{c}}
{}&{}
\end{array}0 \le t \le 2\\
1 + \frac{t}{2},\begin{array}{*{20}{c}}
{}&{}
\end{array}t > 2
\end{array} \right.
\end{array}\]

[控制理論] 線性化(Linearization)

這次要跟大家介紹的是 線性化 (Linearization) 的概念，讀者建議須先具備基本 Taylor Series 概念，如果不熟悉的讀者可先參閱 [微積分] 泰勒展開式 與 泰勒級數 。

為何要做線性化?
其實線性化的動機很簡單，主要是因為一般在分析動態系統的時候，大部分系統行為都是呈現非線性(EX: 電路系統(二極體 I/V curve)，倒單擺、撓性機構、機器人、生物細胞、金融模型...)，但這些非線性行為會有一個大的困難，就是難以直接求解其動態行為。且發展成熟的線性系統理論沒有辦法(有效的)應用在上面，但如果能夠透過一些假設/機制，我們可以把原本非線性的系統轉成線性系統，如此一來原本沒辦法使用的線性系統理論便可以派上用場!!

如何做線性化?
至於實際如何做到對任意 非線性函數 (e.g., $\sin, \cos, \exp, x^n$, ...)線性化呢? 簡單來說，就是採用切線 (微分) 的概念，如果我們對關心的某一點對該點取導數，則我們可以得到一條對該點的切線，此切線可以在某種程度上用來近似 該點附近的函數行為。

https://controls.engin.umich.edu/wiki/index.php/LinearizingODEs

----- 以下進入正題 ----

若用數學來描述非線性的系統可以寫成
\[
\dot x(t) = f(x)
\]其中 $x(t) \in \mathbb{R}^n$ 稱作系統狀態(state variable) (這邊考慮 $n$ 維空間，故有 $n$ 個系統狀態變數)； $\dot x(t)$ 為系統狀態的一階導數； $f$ 為用以描述動態系統的任意函數

在此我們考慮系統狀態為 $n$ 階。意思就是有 $n$ 個不同的系統狀態，記做  $x \in \mathbb{R}^n$


在介紹線性化之前，我們得先介紹 "平衡點(equilibrium point)"：

=====================
Definition: Equilibrium point
若 $f(\bar{x})=0$ ，則系統狀態 $\bar{x} \in \mathbb{R}^n$ 被稱作 平衡點(equilibrium point)。
=====================

Comments
由上述定義可以推知，如果 $x(0)=\bar{x}$ ，則 $x(t)=\bar{x}, \forall t \geq 0$ ；
也就是說一旦 在最一開始( $t=0$ )的時候，系統就處在平衡點的狀態，則對任意時刻 $t \geq 0$，系統狀態會持續處在平衡點的狀態。

====================
Definition: 穩定平衡點 (Stable equilibrium point)
平衡點若被稱為穩定的，或稱 穩定平衡點(stable equilibrium point)，則其必須滿足 在任意時刻 之狀態 $x(t)$ 都需收斂到平衡點 $\bar{x}$，亦即
\[
x(t) \rightarrow \bar{x}
\](  $|| x(0)-\bar{x} ||$ 為足夠小 )
反之，若不收斂則稱為 不穩定的平衡點(unstable equilibrium) 其中
 \[
|| x(0)-\bar{x} || := \left ( \displaystyle \sum_{i} (x_i(0)-\bar{x}_i)^2 \right )^{\frac{1}{2}}\] 為 2-norm
==================

Comments：
上述定義指明 所謂的 穩定平衡點是指 考慮 任意時刻的狀態 $x(t), \forall t$，若此狀態都會回到 某個平衡點 $\bar{x}$ 則我們說他是一個 穩定的平衡點。

在介紹完平衡點之後，我們便可介紹所謂的 線性化，誠如先前所說，線性化的基本概念是微分，所以在此我們會假設動態系統充分可微(smooth)，故我們可以進行 泰勒展開 (微分近似)。

=================
線性化(Linearization)：
注意：我們僅對 穩定平衡點 做線性化。(不穩定的平衡點亦可線性化只是實際用處不大)

現在我們回頭考慮 $n$ 階非線性系統
\[
\dot x(t) = f(x)
\]其系統狀態 $x(t)$ 可表為 平衡點狀態 $\bar{x}$ 加上 狀態(小擾動)增量( $\Delta x(t)$ )；注意。在此我們假設擾動 $\Delta x$ 不能太大。
\[
x(t) = \bar{x} + \Delta x(t)
\]則我們可寫下
\[
f_1 (x) = f_1 (x+\Delta x)\]
若此 $f$ 為 平滑函數(smooth) (也就是說可以對其做泰勒展開)，則我們可改寫上式如下:

對 $f_1$ 可寫出其泰勒展開式 (對 $0$ 點展開)

$\Rightarrow f_1(x) = f_1(\bar{x}) + \frac{\partial f_1}{\partial x_1}|_{x=\bar{x}} \Delta x_1  + ... + \frac{\partial f_1}{\partial x_n}|_{x=\bar{x}} \Delta x_n + H.O.T$ ....(1)

同樣的，我們也可以對 $f_2...f_n$ 展開。

$\Rightarrow f_2(x) = f_2(\bar{x}) + \frac{\partial f_2}{\partial x_1}|_{x=\bar{x}} \Delta x_1  + ... + \frac{\partial f_2}{\partial x_n}|_{x=\bar{x}} \Delta x_n + H.O.T.$

$\vdots$

$\Rightarrow f_n(x) = f_n(\bar{x}) + \frac{\partial f_n}{\partial x_1}|_{x=\bar{x}} \Delta x_1  + ... + \frac{\partial f_n}{\partial x_n}|_{x=\bar{x}} \Delta x_n + H.O.T$

其中 $H.O.T$ 表示 高階項(Higher Order Terms)

然後因為增量 $\Delta x_1, \Delta x_2...$ 假設為很小的擾動，在高階項的影響可被忽略

現在，回憶我們手邊有的狀態
\[
x(t) = \bar{x} + \Delta x(t)\]
對上式兩邊對時間微分，可得
\[
\dot x(t) = \Delta \dot x(t)\]

再者，因為我們知道 系統為 $\dot x(t) = f(x)$， 由式 (1) 我們可以帶入泰勒展開到 $f(x)$ 之中，最後整理可得線性化之後的 增量(擾動)系統

\[
 \Delta \dot x(t) = A \cdot \Delta x(t)  \ \ \ \  (2) \]

其中 $A$ 為矩陣其第 (i,j) 元素由下式表示

\[
a_{ij} = \frac{\partial f_i}{\partial x_j}|_{x=\bar{x}}\]

上式 $(2)$ 即稱為 線性化後的動態系統。由於此為線性，故所有的線性系統理論 (eigenvalue, controllability, observability) 便可以在其上進行討論

[線性系統] 動態方程式的求解(3) - LTV state equation- Total Solution

延續前篇文章 [線性系統] 動態方程式的求解(2) - LTV state equation- Homogeneous solution，這次要介紹線性時變 (Linear Time Varying, LTV ) 系統的狀態方程的全解。


考慮下列 LTV 動態系統
\[\left\{ {\begin{array}{*{20}{l}}
{{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)}\\
{{\bf{y}}\left( t \right) = {\bf{C}}\left( t \right){\bf{x}}\left( t \right) + {\bf{D}}\left( t \right){\bf{u}}\left( t \right)}
\end{array}} \right.
\] 且假設  ${\bf{A}}\left( t \right)$ 為 $n \times n$ 且矩陣中每一項元素 都為對時間 $t$ 連續函數。

Comment:
1. 上式中 ${{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)}$ 稱為狀態方程 (State equation)
2. ${{\bf{y}}\left( t \right) = {\bf{C}}\left( t \right){\bf{x}}\left( t \right) + {\bf{D}}\left( t \right){\bf{u}}\left( t \right)}$ 稱為 輸出方程 (Output equation)

===================
Claim:
給定初始狀態 ${\bf{x}}\left( {{t_0}} \right)$ 與 輸入 ${{\bf{u}}\left( t \right)}$，則狀態方程 ${{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)}$ 的解為
\[
{\bf{x}}\left( t \right) = {\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } \ \ \ \ (*)
\]其中  ${\bf{\Phi }}\left( {t,\tau } \right): = {\bf{X}}\left( t \right){{\bf{X}}^{ - 1}}\left( \tau  \right)$ 為 ${\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right)$ 的 State Transition matrix 滿足\[\frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,{t_0}} \right) = {\bf{A}}\left( t \right){\bf{\Phi }}\left( {t,{t_0}} \right)\]且 初始條件為 ${\bf{\Phi }}\left( {{t_0},{t_0}} \right) = {\bf{I}}$。
===================

Proof:
首先證明 $(*)$ 滿足初始條件：
\[\begin{array}{l}
{\bf{x}}\left( t \right) = {\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } \\
 \Rightarrow {\bf{x}}\left( {{t_0}} \right) = {\bf{\Phi }}\left( {{t_0},{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \underbrace {\int_{{t_0}}^{{t_0}} {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } }_{ = 0}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = {\bf{X}}\left( {{t_0}} \right){{\bf{X}}^{ - 1}}\left( {{t_0}} \right){\bf{x}}\left( {{t_0}} \right) = {\bf{Ix}}\left( {{t_0}} \right) = {\bf{x}}\left( {{t_0}} \right)
\end{array}\]接著我們證明 $(*)$ 確實滿足狀態方程。
\[\begin{array}{l}
{\bf{x}}\left( t \right) = {\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } \\
\frac{d}{{dt}}{\bf{x}}\left( t \right) = \frac{d}{{dt}}\left[ {{\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } } \right]\\
 \Rightarrow {\bf{\dot x}}\left( t \right) = \frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \frac{\partial }{{\partial t}}\left[ {\int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } } \right]
\end{array}
\] 利用 Fundamental Theorem of Calculus:
\[\frac{\partial }{{\partial t}}\int_{{t_0}}^t {f\left( {t,\tau } \right)d\tau }  = \left. {f\left( {t,\tau } \right)} \right|_{\tau  = t}^{} + \int_{{t_0}}^t {\left( {\frac{\partial }{{\partial t}}f\left( {t,\tau } \right)} \right)d\tau }
\] 我們得知
\[\begin{array}{l}
{\bf{\dot x}}\left( t \right) = \frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) \\
\ \ \ \ \ \ \ \ \ \ \ + \left[ {{\bf{\Phi }}\left( {t,t} \right){\bf{B}}\left( t \right){\bf{u}}\left( t \right) + \int_{{t_0}}^t {\left( {\frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)} \right)d\tau } } \right]\\
 \Rightarrow {\bf{\dot x}}\left( t \right) = \frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + {\bf{\Phi }}\left( {t,t} \right){\bf{B}}\left( t \right){\bf{u}}\left( t \right) \\
\ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \  \ \ \  + \int_{{t_0}}^t {\left( {\frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)} \right)d\tau }
\end{array}
\]再由 State Transition Matrix 定義 $\frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,{t_0}} \right) = {\bf{A}}\left( t \right){\bf{\Phi }}\left( {t,{t_0}} \right)$我們知道
\[\begin{array}{l}
{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + {\bf{\Phi }}\left( {t,t} \right){\bf{B}}\left( t \right){\bf{u}}\left( t \right) \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ + \int_{{t_0}}^t {\left( {\frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,\tau } \right)} \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } \\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + {\bf{\Phi }}\left( {t,t} \right){\bf{B}}\left( t \right){\bf{u}}\left( t \right)\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  + \int_{{t_0}}^t {{\bf{A}}\left( t \right){\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } \\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right)\underbrace {\left[ {{\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right)  + \int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } } \right]}_{ = {\bf{x}}\left( t \right)} \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ + \underbrace {{\bf{X}}\left( t \right){{\bf{X}}^{ - 1}}\left( t \right)}_{ = {\bf{I}}}{\bf{B}}\left( t \right){\bf{u}}\left( t \right)\\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)
\end{array}
\]

有了上述結果之後，我們便可以進一步求得 輸入輸出之間關係，將
\[
{\bf{x}}\left( t \right) = {\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau } \] 帶回輸出方程 ${{\bf{y}}\left( t \right) = {\bf{C}}\left( t \right){\bf{x}}\left( t \right) + {\bf{D}}\left( t \right){\bf{u}}\left( t \right)}$，可得
\[\begin{array}{l}
 \Rightarrow {\bf{y}}\left( t \right) = {\bf{C}}\left( t \right){\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} + {\bf{C}}\left( t \right)\int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)d\tau }  + {\bf{D}}\left( t \right){\bf{u}}\left( t \right)
\end{array}\]

[線性系統] 動態方程式的求解(2) - LTV state equation- Homogeneous solution

這次要介紹線性時變 (Linear Time Varying, LTV ) 系統的狀態方程求解。

考慮下列 LTV 動態系統
\[\left\{ {\begin{array}{*{20}{l}}
{{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)}\\
{{\bf{y}}\left( t \right) = {\bf{C}}\left( t \right){\bf{x}}\left( t \right) + {\bf{D}}\left( t \right){\bf{u}}\left( t \right)}
\end{array}} \right.
\] 且假設  ${\bf{A}}\left( t \right)$ 為 $n \times n$ 且矩陣中每一項元素 都為對時間 $t$ 連續函數。

NOTE: 若上述對 ${\bf{A}}\left( t \right)$ 時變矩陣的連續性假設成立，則對任意初始狀態 ${\bf{x}}\left( {{t_0}} \right)$ 與任意輸入 ${{\bf{u}}\left( t \right)}$， 狀態方程有唯一解。
(Proof ommitted)

在我們進行求解之前，我們首先求解
\[
{{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right)}
\] 其中 ${\bf{A}}\left( t \right)$ 為 $n \times n$ 且每一個 entry 都為 對時間 $t$ 連續的函數。故對任意初始狀態 ${\bf{x}}_i\left( {{t_0}} \right)$  狀態方程存在唯一解 $ {{\bf{x}}_i}\left( t \right),\forall i = 1,2,...,n$ 。

我們可以將這些 $n$ 個解蒐集起來寫作矩陣形式如下：
\[{\bf{X}}\left( t \right): = \left[ {\begin{array}{*{20}{c}}
{{{\bf{x}}_1}\left( t \right)}&{{{\bf{x}}_2}\left( t \right)}& \cdots &{{{\bf{x}}_n}\left( t \right)}
\end{array}} \right]
\]由於 每一個 $ {{\bf{x}}_i}\left( t \right)$ 都滿足 ${{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right)}$ 故我們有
\[{\bf{\dot X}}\left( t \right) = {\bf{A}}\left( t \right){\bf{X}}\left( t \right)
\]

現在我們給出下面的定義：
====================
Definition: (Fundamental Matrix)
若 ${\bf{X}}\left( {{t_0}} \right)$ 為 nonsingular 或者 $n$ 個初始狀態彼此之間為線性獨立，則時變矩陣 ${\bf{X}}\left( t \right)$ 稱作 ${{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right)}$ 的 Fundamental matrix 。
====================

Comment:
Fundamental matrix 並無唯一表示式 (因為初始狀態可以任選).

接著我們定義 狀態轉移矩陣 (State Transition Matrix)

====================
Definition: (State Transition Matrix)
令 ${\bf{X}}\left( {{t}} \right)$ 為 ${{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right)}$ 的 Fundamental matrix，則我們定義其對應的 狀態轉移矩陣 (State Transition Matrix) ${\bf{\Phi }}\left( {t,{t_0}} \right)$ 如下：
\[
{\bf{\Phi }}\left( {t,{t_0}} \right): = {\bf{X}}\left( t \right){{\bf{X}}^{ - 1}}\left( {{t_0}} \right)
\] 且 此狀態轉移矩陣 ${\bf{\Phi }}\left( {t,{t_0}} \right)$ 為 下列狀態方程的唯一解
\[\frac{\partial }{{\partial t}}{\bf{\Phi }}\left( {t,{t_0}} \right) = {\bf{A}}\left( t \right){\bf{\Phi }}\left( {t,{t_0}} \right)\]且 初始條件為 ${\bf{\Phi }}\left( {{t_0},{t_0}} \right) = {\bf{I}}$。
====================

====================
Theorem:
給定任意初始狀態 $t_0$，狀態方程 ${{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right)}$  的解為
\[{\bf{x}}\left( t \right) = {\bf{\Phi }}\left( {{t},{t_0}} \right){\bf{x}}\left( {{t_0}} \right)
\]====================

Proof: Omitted.

下面我們看個例子看看給定狀態方程  ${{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right)}$  如何求出對應的 Fundamental matrix。以及 State Transition Matrix。

====================
Example
考慮下列狀態方程
\[{\bf{\dot x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
0&0\\
t&0
\end{array}} \right]{\bf{x}}\left( t \right)
\]試求對應的 Fundamental matrix，State Transition Matrix，與 ${\bf{x}}\left( t \right)$。
====================

Solution:
由於 時變矩陣 ${\bf{A}}\left( t \right)$ 符合連續性假設，故我們有對任意初始狀態 ${\bf{x}}\left( {{t_0}} \right)$，存在唯一解。

現在我們觀察
\[{\bf{\dot x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
0&0\\
t&0
\end{array}} \right]{\bf{x}}\left( t \right) \Rightarrow \left[ {\begin{array}{*{20}{c}}
{{{\dot x}_1}\left( t \right)}\\
{{{\dot x}_2}\left( t \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
0&0\\
t&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_1}\left( t \right)}\\
{{x_2}\left( t \right)}
\end{array}} \right]
\]亦即
\[\left\{ \begin{array}{l}
{{\dot x}_1}\left( t \right) = 0\\
{{\dot x}_2}\left( t \right) = t{x_1}\left( t \right)
\end{array} \right.
\]故給定初始時間 $t_0 =0$ 我們可求解 $x_1(t)$ 與 $x_2(t)$ 如下
\[\begin{array}{l}
{{\dot x}_1}\left( t \right) = 0\\
 \Rightarrow \int_0^t {d{x_1}\left( \tau  \right)}  = 0\\
 \Rightarrow {x_1}\left( t \right) = {x_1}\left( 0 \right)
\end{array}\]與
\[\begin{array}{l}
{{\dot x}_2}\left( t \right) = t{x_1}\left( t \right)\\
 \Rightarrow \int_0^t {d{x_2}\left( \tau  \right)}  = \int_0^t {\tau {x_1}\left( 0 \right)d\tau } \\
 \Rightarrow {x_2}\left( t \right) = {x_1}\left( 0 \right)\frac{{{t^2}}}{2} + {x_2}\left( 0 \right)
\end{array}
\]為了建構 Fundamental matrix，我們可任意選定 等同時變矩陣階數數目的初始狀態，在此例中由於 $\bf{A}$ 為 $2 \times 2$ 時變矩陣，故我們可任選兩個 線性獨立的 初始狀態 來建構 Fundamental matrix ${\bf{X}}\left( t \right)$，比如說選 ${\bf{x}}\left( 0 \right)=[1 \; 0]^T$ 與 ${\bf{x}}\left( 0 \right)=[0 \; 1]^T$，則我們有
\[\begin{array}{l}
{\bf{x}}\left( 0 \right) = \left[ {\begin{array}{*{20}{c}}
{{x_1}\left( 0 \right)}\\
{{x_2}\left( 0 \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
1\\
0
\end{array}} \right]\\
 \Rightarrow {\bf{x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
{{x_1}\left( t \right)}\\
{{x_2}\left( t \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{x_1}\left( 0 \right)}\\
{{x_1}\left( 0 \right)\frac{{{t^2}}}{2} + {x_2}\left( 0 \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
1\\
{\frac{{{t^2}}}{2}}
\end{array}} \right]
\end{array}
\]與
\[\begin{array}{l}
{\bf{x}}\left( 0 \right) = \left[ {\begin{array}{*{20}{c}}
{{x_1}\left( 0 \right)}\\
{{x_2}\left( 0 \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
0\\
1
\end{array}} \right]\\
 \Rightarrow {\bf{x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
{{x_1}\left( t \right)}\\
{{x_2}\left( t \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{x_1}\left( 0 \right)}\\
{{x_1}\left( 0 \right)\frac{{{t^2}}}{2} + {x_2}\left( 0 \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
0\\
1
\end{array}} \right]
\end{array}
\] 由於 $[1 \; 0]^T$ 與 $[0 \; 1]^T$ 彼此線性獨立，故由 Fundamental matrix 的定義，我們確實得到了 一組 (不唯一) Fundamental Matrix 如下：
\[{\bf{X}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{{t^2}}}{2}}&1
\end{array}} \right]
\] 故由 State Transition Matrix 定義，我們可知
\[
{\bf{\Phi }}\left( {t,{t_0}} \right): = {\bf{X}}\left( t \right){{\bf{X}}^{ - 1}}\left( {{t_0}} \right)
\]其中
\[{{\bf{X}}^{ - 1}}\left( {{t_0}} \right) = {\left. {\frac{1}{1} \cdot \left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{ - {t^2}}}{2}}&1
\end{array}} \right]} \right|_{t = {t_0}}} = \left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{ - {t_0}^2}}{2}}&1
\end{array}} \right]
\]故 State Transition Matrix 為
\[\begin{array}{l}
{\bf{\Phi }}\left( {{t},{t_0}} \right) = {\bf{X}}\left( t \right){{\bf{X}}^{ - 1}}\left( {{t_0}} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{{t^2}}}{2}}&1
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{ - {t_0}^2}}{2}}&1
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{{t^2} - {t_0}^2}}{2}}&1
\end{array}} \right]
\end{array}
\] 由前述 Theorem 可知，狀態方程 ${\bf{\dot x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
0&0\\
t&0
\end{array}} \right]{\bf{x}}\left( t \right)$ 的解為
\[
{\bf{x}}\left( t \right) = {\bf{\Phi }}\left( {t,{t_0}} \right){\bf{x}}\left( {{t_0}} \right) = \left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{{t^2} - {t_0}^2}}{2}}&1
\end{array}} \right]{\bf{x}}\left( {{t_0}} \right) = \left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{{t^2}}}{2}}&1
\end{array}} \right]{\bf{x}}\left( 0 \right)
\]現在我們帶回驗證 上式 ${\bf{\Phi }}\left( {t,{t_0}} \right)$ 確實為 ，亦即對其微分
\[\begin{array}{l}
\frac{d}{{dt}}{\bf{x}}\left( t \right) = \frac{d}{{dt}}\left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{{t^2}}}{2}}&1
\end{array}} \right]{\bf{x}}\left( 0 \right) = \left[ {\begin{array}{*{20}{c}}
0&0\\
t&0
\end{array}} \right]\underbrace {{{\left[ {\begin{array}{*{20}{c}}
1&0\\
{\frac{{{t^2}}}{2}}&1
\end{array}} \right]}^{ - 1}}{\bf{x}}\left( t \right)}_{{\bf{x}}\left( 0 \right)}\\
 \Rightarrow {\bf{\dot x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
0&0\\
t&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
1&0\\
{ - \frac{{{t^2}}}{2}}&1
\end{array}} \right]{\bf{x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
0&0\\
t&0
\end{array}} \right]{\bf{x}}\left( t \right)
\end{array}
\]故得證。

[線性規劃] Discrete Time Dynamic System Linear Programming

這次要介紹的是如何把一個 有拘束條件 離散時間的動態系統
透過 線性規劃問題來求解。

考慮 離散時間系統的狀態空間表示(state space representation)：
\[
x(k+1) = A x(k) + B u(k), \ \ \ \ \text{given $x(0)$}
\]其中 $x(k) \in \mathbb{R}^n$ 表系統狀態， $u(k)$ 表控制力。

接著我們考慮一些實際的情況：

1. 我們希望系統狀態 隨著時間流逝 與 控制力的幫助，可以在某個時間點達到某個 指定的目標狀態 $x^1$，亦即考慮 終止狀態 $x(N) := x^1$
2. 由於我們是透過控制力來幫助我們讓系統狀態逐步移動到 指定的目標狀態，考量到一班情況，控制力的出力大小須受到一定程度拘束(不可以有無限大的控制力)，亦即
 \[
|u(k)| \leq M \Leftrightarrow -M \leq u(k) \leq M
\]

現在我們定義目標成本函數 (cost function)
\[
J(u) = \displaystyle \sum_{k=0}^{N-1} | u(k)|
\]

Comments:
1. 上述目標為 讓 控制力最小。簡單來說就是要 可以達到目標狀態的情況下 盡可能省力 (可以想像成要要開車 從 A 到 B地點，且盡可能省油。)

2. 上述cost function 最多只能是 取絕對值 (儘管是非線性但可透過一些技巧將其轉換成LP問題)，超過  一般 2階 (or 以上)的 cost function，比如
\[
J(u) = \displaystyle \sum_{k=0}^{N-1} | u(k)|^2
\]不可使用 Linear Programming。因為上述cost function 不再是線性! 。Forget LP in this case..

故我們手上有了一個標準的最佳化問題：
\[\begin{array}{l}
\min J\left( u \right) = \sum\limits_{k = 0}^{N - 1} {\left| {u\left( k \right)} \right|} \\
s.t.\\
\begin{array}{*{20}{c}}
{}
\end{array}x(k + 1) = Ax(k) + Bu(k)\\
\begin{array}{*{20}{c}}
{}
\end{array}x(N): = {x^1}\\
\begin{array}{*{20}{c}}
{}
\end{array}|u(k)| \le M
\end{array}
\]

首先注意到因為我們的成本函數是對 $u(k)$ 最小化，故上式中的拘束應改寫成由 $u(k)$表示：
故我們先改寫 目標狀態拘束 $x(N) = x^1$ 用 $u(k)$ 表示：

怎麼做呢!? 由 離散時間系統的狀態方程
\[
x(k+1) = A x(k) + B u(k), \ \ \ \ \text{given $x(0)$}
\]我們可以進行 跌代求解如下：
\[\begin{array}{l}
x(1) = Ax(0) + Bu(0)\\
x(2) = Ax(1) + Bu(1) = A\left[ {Ax(0) + Bu(0)} \right] + Bu(1)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = {A^2}x(0) + ABu(0) + Bu(1)\\
x(3) = Ax(2) + Bu(2) = A\left[ {{A^2}x(0) + ABu(0) + Bu(1)} \right] + Bu(2)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = {A^3}x(0) + {A^2}Bu(0) + ABu(1) + Bu(2)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} \vdots \\
x(N) = {A^N}x(0) + {A^{N - 1}}Bu(0) +  \cdots  + ABu(N - 2) + Bu(N - 1)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = {A^N}x(0) + \sum\limits_{k = 0}^{N - 1} {{A^{N - 1 - k}}Bu(k)}
\end{array}
\]又因為 $x(N) := x^1$，故我們得到
\[\begin{array}{l}
{x^1} = {A^N}x(0) + \sum\limits_{k = 0}^{N - 1} {{A^{N - 1 - k}}Bu(k)} \\
 \Rightarrow {x^1} - {A^N}x(0) = \left[ {\begin{array}{*{20}{c}}
{{A^{N - 1}}B}&{{A^{N - 2}}B}& \cdots &{{A^1}B}&{{A^0}B}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{u(0)}\\
{u(1)}\\
 \vdots \\
{u(N - 1)}
\end{array}} \right]
\end{array}
\] 由於上式左方 $x^1$ 與 $x(0)$ 已給定，故我們可以聯合上式與 控制力拘束條件 $ |u(k)| \leq M$，將其改寫成LP問題：

對於每一個固定的 $N$，求解 進行求解 LP問題:
\[\begin{array}{l}
\min J\left( u \right) = \sum\limits_{k = 0}^{N - 1} {\left| {u\left( k \right)} \right|} \\
s.t.\\
\left\{ \begin{array}{l}
\begin{array}{*{20}{c}}
{}
\end{array}{x^1} - {A^N}x(0) = \left[ {\begin{array}{*{20}{c}}
{{A^{N - 1}}B}&{{A^{N - 2}}B}& \cdots &{{A^1}B}&{{A^0}B}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{u(0)}\\
{u(1)}\\
 \vdots \\
{u(N - 1)}
\end{array}} \right]\\
\begin{array}{*{20}{c}}
{}
\end{array}|u(k)| \le M
\end{array} \right.
\end{array}\]

[線性系統] 離散時間狀態空間模型 與其解

一般若需要將控制系統透過 電腦 實現控制力 或者 對連續時間系統進行取樣，則我們稱此類系統為 電腦控制系統 或稱 數位控制系統。這次我們要介紹如何從連續時間模型 將其 透過適當的數學操作，從而獲得其對應的離散化的模型。

現在考慮有限維 線性非時變 連續時間狀態空間模型如下：
\[\left\{ \begin{array}{l}
\dot x\left( t \right) = A_cx\left( t \right) + B_cu\left( t \right);\begin{array}{*{20}{c}}
{}&{}
\end{array}x\left( 0 \right) = {x_0}\\
y\left( t \right) = C_c x\left( t \right) + D_cu\left( t \right)
\end{array} \right. \ \ \ \ \ \ \ \ \ (\star)
\]注意到上述狀態空間模型含有微分項 $\dot x$，若我們想要透過電腦實現微分方程，則我們需將其進行 離散化(Discretization)，
\[\dot x\left( t \right): = \mathop {\lim }\limits_{\Delta  \to 0} \frac{{x\left( {t + \Delta } \right) - x\left( t \right)}}{\Delta }\]則前述連續時間的狀態方程 $(\star)$ 可表為
\[\begin{array}{l}
\dot x\left( t \right) = {A_c}x\left( t \right) + {B_c}u\left( t \right)\\
 \Rightarrow x\left( {t + \Delta } \right) - x\left( t \right) = {A_c}x\left( t \right)\Delta  + {B_c}u\left( t \right)\Delta \\
 \Rightarrow x\left( {t + \Delta } \right) = x\left( t \right) + {A_c}x\left( t \right)\Delta  + {B_c}u\left( t \right)\Delta
\end{array}\]前述的 $\Delta$ 一般在實現上稱為取樣時間 (sampling time)，且對 $k \in \mathbb{Z}$， $t := k \Delta$，我們可改寫
\[\begin{array}{l}
\left\{ {\begin{array}{*{20}{l}}
{x\left( {k\Delta  + \Delta } \right) = \left( {I + {A_c}\Delta } \right)x\left( {k\Delta } \right) + {B_c}u\left( {k\Delta } \right)\Delta }\\
{y\left( {k\Delta } \right) = {C_c}x\left( {k\Delta } \right) + {D_c}u\left( {k\Delta } \right)}
\end{array}} \right.\\
 \Rightarrow \left\{ {\begin{array}{*{20}{l}}
{x\left( {\left( {k + 1} \right)\Delta } \right) = \left( {I + {A_c}\Delta } \right)x\left( {k\Delta } \right) + {B_c}u\left( {k\Delta } \right)\Delta }\\
{y\left( {k\Delta } \right) = {C_c}x\left( {k\Delta } \right) + {D_c}u\left( {k\Delta } \right)}
\end{array}} \right.
\end{array}\]上式現在可以很容易地透過 MATLAB 實現。

更進一步，若考慮 $u(t)$ 也透過電腦計算而得，再透過數位類比轉換 (Digital-Analog converter ) 或 零階保持(Zero-Order Hold)，則 $u(t)$ 將會是 piecewise constant 。也就是說對任意 $t$ 滿足 $k\Delta \le t < (k+1) \Delta$，我們的控制力可定義為
$$u(t) := u(k \Delta) := u(k) \;\;\; \forall k=0,1,2,...$$ 此時控制力 $u(k)$ 將只會在 取樣點 $k\Delta$ 上有變動，且在取樣點上，我們連續時間系統的解仍然成立，故對於 $t = k \Delta$ ，我們可定義狀態方程的解
\[\begin{array}{*{20}{l}}
{x\left( k \right): = x\left( {k\Delta } \right) = {{\left. {x\left( t \right)} \right|}_{t = k\Delta }}}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {{\left. {\left( {{e^{{A_c}t}}{x_0} + \int_0^t {{e^{{A_c}\left( {t - \tau } \right)}}} Bu\left( \tau  \right)d\tau } \right)} \right|}_{t = k\Delta }}}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {e^{{A_c}k\Delta }}{x_0} + \int_0^{k\Delta } {{e^{{A_c}\left( {k\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau }
\end{array}
\]且 對於 $t = (k+1) \Delta$ 而言，我們亦可定義
\[\small \begin{array}{*{20}{l}}
{x\left( {k + 1} \right): = x\left( {\left( {k + 1} \right)\Delta } \right) = {{\left. {x\left( t \right)} \right|}_{t = \left( {k + 1} \right)\Delta }}}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {{\left. {\left( {{e^{{A_c}t}}{x_0} + \int_0^t {{e^{{A_c}\left( {t - \tau } \right)}}} Bu\left( \tau  \right)d\tau } \right)} \right|}_{t = \left( {k + 1} \right)\Delta }}}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {e^{{A_c}\left( {k + 1} \right)\Delta }}{x_0} + \int_0^{\left( {k + 1} \right)\Delta } {{e^{{A_c}\left( {\left( {k + 1} \right)\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau }\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {e^{{A_c}\left( {k + 1} \right)\Delta }}{x_0} + \int_0^{k\Delta } {{e^{{A_c}\left( {\left( {k + 1} \right)\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau  + \int_{k\Delta }^{\left( {k + 1} \right)\Delta } {{e^{{A_c}\left( {\left( {k + 1} \right)\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau }\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {e^{{A_c}\Delta }}\underbrace {\left[ {{e^{{A_c}k\Delta }}{x_0} + \int_0^{k\Delta } {{e^{{A_c}\left( {k\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau } \right]}_{ = x\left( k \right)} + \int_{k\Delta }^{\left( {k + 1} \right)\Delta } {{e^{{A_c}\left( {\left( {k + 1} \right)\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau }\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {e^{{A_c}\Delta }}x\left( k \right) + \int_{k\Delta }^{\left( {k + 1} \right)\Delta } {{e^{{A_c}\left( {\left( {k + 1} \right)\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau \;\;\;\;\;(**)}
\end{array}\]現在若令新的變數 $\alpha : = \left( {k + 1} \right)\Delta  - \tau $ 則 $(**)$ 可被改寫為
\[ \begin{array}{l}
x\left( {k + 1} \right) = {e^{A\Delta }}x\left( k \right) + \int_{\left( k \right)\Delta }^{\left( {k + 1} \right)\Delta } {{e^{A\left( {\left( {k + 1} \right)\Delta  - \tau } \right)}}} Bu\left( \tau  \right)d\tau \\
 \Rightarrow x\left( {k + 1} \right) = {e^{A\Delta }}x\left( k \right) + \left( {\int_0^\Delta  {{e^{A\alpha }}} d\alpha } \right)Bu\left( k \right)
\end{array}\]注意到上式中第二行 $u(k)$ 被提到積分外面是因為 $u(t)$ 在 $k\Delta \le t < (k+1)\Delta$ 之間為 constant。

故總結上式，我們有
\[\left\{ \begin{array}{l}
x\left( {k + 1} \right) = \underbrace {{e^{{A_c}\Delta }}}_Ax\left( k \right) + \underbrace {\left( {\int_0^\Delta  {{e^{{A_c}\alpha }}} d\alpha } \right)B}_Bu\left( k \right)\\
y\left( k \right) = Cx\left( k \right) + Du\left( k \right)
\end{array} \right.\] 上式即稱為對應的 零階自保持的 有限維 線性非時變 離散時間狀態空間模型 。

若 $A_c$ 反矩陣存在，則上述的 $B$ 矩陣 可更進一步計算如下，
\[\begin{array}{l}
B = \left( {\int_0^\Delta  {{e^{{A_c}\alpha }}d\alpha } } \right)B\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left( {\int_0^\Delta  {\left( {I + {A_c}\alpha  + \frac{{{{\left( {{A_c}\alpha } \right)}^2}}}{{2!}} + ...} \right)d\alpha } } \right)B\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left( {\int_0^\Delta  {\left( I \right)d\alpha }  + \int_0^\Delta  {\left( {{A_c}\alpha } \right)d\alpha }  + \int_0^\Delta  {\left( {\frac{{{{\left( {{A_c}\alpha } \right)}^2}}}{{2!}}} \right)d\alpha }  + ...} \right)B\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left( {I\Delta  + {A_c}\frac{{{\Delta ^2}}}{2} + {A_c}^2\frac{{{\Delta ^3}}}{{3!}} + ...} \right)B\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left( {\sum\limits_{k = 1}^\infty  {{A_c}^{k - 1}\frac{{{\Delta ^k}}}{{k!}}} } \right)B\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \left( {{A_c}^{ - 1}\sum\limits_{k = 1}^\infty  {\frac{{{{\left( {{A_c}\Delta } \right)}^k}}}{{k!}}} } \right)B = \left( {{A_c}^{ - 1}\left( {{e^{{A_c}\Delta }} - I} \right)} \right)B
\end{array}\]


現在我們可以開始求解離散化的狀態空間方程
\[\left\{ {\begin{array}{*{20}{l}}
{x\left( {k + 1} \right) = Ax\left( k \right) + Bu\left( k \right)}\\
{y\left( k \right) = Cx\left( k \right) + Du\left( k \right)}
\end{array}} \right.\]首先觀察
\[\left\{ \begin{array}{l}
x\left( 1 \right) = Ax\left( 0 \right) + Bu\left( 0 \right)\\
x\left( 2 \right) = Ax\left( 1 \right) + Bu\left( 1 \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {A^2}x\left( 0 \right) + ABu\left( 0 \right) + Bu\left( 1 \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {A^2}x\left( 0 \right) + \sum\limits_{j = 0}^1 {{A^{1 - j}}Bu\left( j \right)} \\
x\left( 3 \right) = Ax\left( 2 \right) + Bu\left( 2 \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {A^3}x\left( 0 \right) + {A^2}Bu\left( 0 \right) + ABu\left( 1 \right) + Bu\left( 2 \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = {A^3}x\left( 0 \right) + \sum\limits_{j = 0}^2 {{A^{2 - j}}Bu\left( j \right)} \\
 \vdots
\end{array} \right.\]故
\[x\left( k \right) = {A^k}x\left( 0 \right) + \sum\limits_{j = 0}^{k - 1} {{A^{k - 1 - j}}Bu\left( j \right)} \]

[線性系統] 離散時間 LTI 系統的漸進穩定度

考慮離散時間系統
\[
x(k+1) = Ax(k) + Bu(k)
\]若 $A$ 為 穩定矩陣，且 $u(k) \to 0$ 則 $x(k) \to 0$

Proof:
我們要證明  $x(k) \to 0$，故給定任意 $\varepsilon>0$，要證明 存在 $M>0$ 使得 對任意 $k \ge M$，我們有
\[
|x(k)| \le \varepsilon
\]
注意到該系統 $x(k+1) = Ax(k) + Bu(k)$ 的解為
\[{x(k) = {A^k}x(0) + \sum\limits_{i = 0}^{k - 1} {{A^{k - 1 - j}}Bu(j)} }
\]兩邊同取 norm 並利用三角不等式 可得
\[\begin{array}{*{20}{l}}
{\left| {x(k)} \right| = \left| {{A^k}x(0) + \sum\limits_{i = 0}^{k - 1} {{A^{k - 1 - j}}Bu(j)} } \right|}\\
{\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} \le \left| {{A^k}} \right|\left| {x(0)} \right| + \sum\limits_{i = 0}^{k - 1} {\left| {{A^{k - 1 - j}}} \right|\left| B \right|} \left| {u(j)} \right|}
\end{array}\ \ \ \ \ (*)
\]回憶 Horn 與 Johnson (1985) 的結果：
==================
FACT:
\[
|A^k| \le c \lambda^k, \; c>0 \;\; \max_i |eig_i(A)| < \lambda <1
\]==================
故
\[\left| {x(k)} \right| \le c{\lambda ^k}\left| {x(0)} \right| + c\left| B \right|\sum\limits_{j = 0}^{\infty} {{\lambda ^{k - 1 - j}}\left| {u(j)} \right|} \]現在利用已知假設 $u(k) \to 0$ 我們可推知必存在 $ N > 0$ 使得 對任意 $k \ge N$，我們有
\[\left| {u(k)} \right| \le \frac{{\varepsilon \left( {1 - \lambda } \right)}}{{2c{\lambda ^k}\left| B \right|}}\]將上述結果代入 $(*)$ 可得
\[\begin{array}{l}
\left| {x(k)} \right| \le c{\lambda ^k}\left| {x(0)} \right| + c\left| B \right|\sum\limits_{j = 0}^\infty  {{\lambda ^{k - 1 - j}}\frac{{\varepsilon \left( {1 - \lambda } \right)}}{{2c{\lambda ^k}\left| B \right|}}} \\
 \Rightarrow \left| {x(k)} \right| \le c{\lambda ^k}\left| {x(0)} \right| + \frac{\varepsilon }{2}
\end{array}\]由於 $\lambda <1$ 我們可選 $M'>0$ 使得 對任意 $k \ge \max (M', N)$
\[
c\lambda^k |x(N)| \le \varepsilon/2, \;\; \forall k \ge m'
\]現在結合前述結果我們可推得 對任意 $k \ge M = \max(M', N)$，
\[\left| {x(k)} \right| \le \varepsilon
\]即為所求。 $\square$

[線性系統] 淺談動態系統的可控制性(1)

一般而言基礎控制理論中有三大重要性質

1. 可控制性 (controllability)
2. 可觀測性 (observability)
3. 穩定性 (stability)

這次主要是介紹 動態系統 的 可控制性(Controllability)。

考慮 $n$ 個狀態 且 $p$ 組輸入的狀態方程
\[{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)
\]其中 ${\bf A}(t)$ 為 $n \times n$ 時變矩陣，${\bf B}(t)$ 為 $n \times p$ 時變矩陣。

控制性的基本想法如下：
若系統某狀態 $\bf x$ 可以透過某對應的控制力 $\bf u$ 來影響 (在有限時間 $t$ 中從任意狀態 $x_0$ 被移動到指定狀態 $x(t)$)，則我們稱此狀態為可控制。

由於輸出方程與系統控制性無關，我們這邊只考慮狀態方程。

以下先給出對於線性非時變系統 其 系統可控制性的定義

=====================
Definition: (Controllability for LTI system)
我們稱狀態方程
\[
{\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)
\] 或者一組 $\left( {{\bf{A,B}}} \right)$ 為在時刻 $t_0$ 可控制(controllable)，若下列條件成立:
對任意初始狀態 ${\bf{x}}\left( t_0 \right) = {\bf x}_0$ 與 終止狀態 ${\bf x}_1$，存在時刻 $t_1 \ge t_0$ 與 一組輸入訊號 ${\bf u}(t)$ 使得可以在有限時間 $t_1$ 內將 ${{\bf{x}}_0}$ 送至 ${{\bf{x}}_1}$。反之我們稱此 $\left( {{\bf{A,B}}} \right)$ 不可控制 (uncontrollable)。
======================

那麼該如何檢驗控制性? 我們給初以下方法：
======================
Theorem: Controllability and Controllability Matrix Test
考慮 $n$ 階 LTI System 表為
\[
{\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)
\] 我們說 系統為 完全狀態可控制 若且為若 下列控制性矩陣(controllability matrix)
\[
\mathcal{C} :=[{\bf B}\;\; {\bf AB}\;\; {\bf A}^2{\bf B}\;\;...\;\; {\bf A}^{n-1}{\bf B}]
\]為 full rank 。(亦即 $\text{rank}({\cal C}) = n$)
======================

Comment: 若 $\text{rank}({\cal C}) = n -1$ 則表示系統有1個狀態不可控；若 $\text{rank}({\cal C}) = n -2$ 表系統有2個狀態不可控；以此類推；若 $\text{rank}({\cal C}) = 0$ 表系統所有狀態不可控。

另外我們還有另一個等價的方法可以檢驗控制性：
不過在介紹之前我們需要先介紹 線性獨立 (linearly independent) 的時間函數 $f_i(t)$ 以及如何判別一組時間函數是否彼此為線性獨立：

======================
Definition: (Linear independence of time functions)
考慮 $f_1(t), f_2(t), ..., f_n(t)$ 為連續函數。我們說此組連續函數為彼此 線性相依 (linear dependent) 若下列條件成立：
存在一組不全為零的純量 $\alpha_1, \alpha_2, ..., \alpha_n$ ，使得
\[
\alpha_1 f_1(t) + \alpha_2 f_2(t) + ... + \alpha_n f_n(t) = 0, \; \forall t \in [t_0,t_1]
\]反之，此組 $f_1(t), f_2(t), ..., f_n(t)$ 稱為彼此線性獨立(linear independent) (此時 所有的係數 $\alpha_i =0$)。
=====================

=====================
Theorem: 
給定 $t_0$ 為初始時刻，$t_1$ 為終止時刻，現在定義
\[
{\bf F}\left( t \right): = \left[ {\begin{array}{*{20}{c}}
{{f_1}\left( t \right)}\\
{{f_2}\left( t \right)}\\
 \vdots \\
{{f_n}\left( t \right)}
\end{array}} \right]
\] 與
\[
{\bf W}(t_1) := \int_{t_0}^{t_1} {\bf F}(\tau) {\bf F}^T(\tau) d \tau,
\]則 $f_i$ 在區間 $[t_0,t_1]$ 之間彼此互為線性獨立 若且為若 $W(t_1)$ 為 nonsingular 矩陣。
=====================

Proof:
首先證明 $(\Rightarrow)$
利用歸謬法(Suppose toward to contradiction)，假設 $W(t_1)$ 為 singular。則 $W(t_1)$ 有線性相依的行或者列向量，故存在一組不全為零的向量 $\alpha$ 使得
\[\alpha^T {\bf W}(t_1) \alpha=0,
\]亦即
\[\begin{array}{l}
{\alpha ^T}{\bf{W}}({t_1})\alpha  = 0\\
 \Rightarrow {\alpha ^T}{\bf{W}}({t_1})\alpha  = {\alpha ^T}\left( {\int_{{t_0}}^{{t_1}} {{\bf{F}}(\tau ){{\bf{F}}^T}(\tau )} d\tau } \right)\alpha  = 0\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \int_{{t_0}}^{{t_1}} {\underbrace {{\alpha ^T}{\bf{F}}(\tau )}_{: = {\bf{\xi }}\left( \tau  \right)}\underbrace {{{\bf{F}}^T}(\tau )\alpha }_{: = {{\bf{\xi }}^T}\left( \tau  \right)}} d\tau  = 0\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \int_{{t_0}}^{{t_1}} {{\bf{\xi }}\left( \tau  \right){{\bf{\xi }}^T}\left( \tau  \right)} d\tau  = \int_{{t_0}}^{{t_1}} {{{\left\| {{\bf{\xi }}\left( \tau  \right)} \right\|}^2}} d\tau  = 0\\
 \Rightarrow {\left\| {{\bf{\xi }}\left( \tau  \right)} \right\|^2} = 0 \Rightarrow {\bf{\xi }}\left( \tau  \right) = 0 \Rightarrow {\alpha ^T}{\bf{F}}(\tau ) = 0
\end{array}
\]
亦即我們有 $f_i$ 之間彼此線性相依。此結果與原本 $f_i$ 線性獨立的假設矛盾。

接著證明 $(\Leftarrow)$
同樣採用歸謬法，假設 $f_i$ 彼此線性相依，故我們知道存在一組不全為零的純量 $\alpha_i, \; i=1,2,...,n$ 使得 $\sum_i \alpha_i f_i =0$。我們可將此結果可寫成向量形式
\[
\alpha^T {\bf F}(t) =0
\]其中 $\alpha^T :=[\alpha_1\;\alpha_2\;...\;\alpha_n]$。現在我們觀察矩陣 ${\bf{W}}({t_1})$，並且對其左邊乘上向量 $\alpha^T$可得
\[\begin{array}{l}
{\alpha ^T}{\bf{W}}({t_1}) = {\alpha ^T}\left( {\int_{{t_0}}^{{t_1}} {{\bf{F}}(\tau ){{\bf{F}}^T}(\tau )} d\tau } \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \int_{{t_0}}^{{t_1}} {\underbrace {{\alpha ^T}{\bf{F}}(\tau )}_{ = 0}{{\bf{F}}^T}(\tau )} d\tau  = 0\\
 \Rightarrow {\alpha ^T}{\bf{W}}({t_1}) = 0
\end{array}
\]上式暗示了 ${\bf{W}}({t_1})$ 矩陣為 Singular，此結果與假設 ${\bf{W}}({t_1})$ 為 nonsingular 矩陣矛盾。

至此證畢。$\square$


那麼現在讓我們回到 線性時變 狀態方程
\[{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)
\]我們知道此狀態方程的解為
\[{\bf{x}}\left( t \right) = {\bf{\Phi }}\left( {{t_1},{t_0}} \right){\bf{x}}\left( {{t_0}} \right) + \int_{{t_0}}^t {{\bf{\Phi }}\left( {t,\tau } \right){\bf{B}}\left( \tau  \right){\bf{u}}\left( \tau  \right)} d\tau
\]其中 ${\bf \Phi}(t,\tau)$ 稱作狀態轉移矩陣 (State transition matrix)。

那麼我們現在便可以把控制性與之前所討論的線性獨立性質聯合再一起，我們將此結果記做下面的重要定理。

======================
Theorem: (Controllability criterion for LTV system)
我們稱 線性時變(LTV) 狀態方程
\[{\bf{\dot x}}\left( t \right) = {\bf{A}}\left( t \right){\bf{x}}\left( t \right) + {\bf{B}}\left( t \right){\bf{u}}\left( t \right)
\]為 在時刻 $t_0$可控制 若且為若 存在一時刻 $t_1 \ge t_0$ 使得 矩陣 ${\bf \Phi(t,\tau)}B(\tau)$ 在區間 $[t_0, t_1]$ 之間 具有線性獨立的 行向量 (row vectors)；或者等價地說，下列 $n \times n$矩陣
\[{\bf{W}}({t_1}): = \int_{{t_0}}^{{t_1}} {{\bf{\Phi }}\left( {{t_1},\tau } \right){\bf{B}}\left( \tau  \right){{\bf{B}}^T}\left( \tau  \right){{\bf{\Phi }}^T}\left( {{t_1},\tau } \right)} d\tau \]為 nonsingular。
=====================

Comment
1. 上述矩陣\[{\bf{W}}({t_1}): = \int_{{t_0}}^{{t_1}} {{\bf{\Phi }}\left( {{t_1},\tau } \right){\bf{B}}\left( \tau  \right){{\bf{B}}^T}\left( \tau  \right){{\bf{\Phi }}^T}\left( {{t_1},\tau } \right)} d\tau \]稱作 Controllability Grammian。

2. 上述定理亦可適用於 LTI system。我們將定理記做下面

======================
Theorem: (Controllability criterion for LTI system)
考慮系統為 LTI，其狀態方程
\[{\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)\]為可控制 若且為若 對任意 $t>0$，下列矩陣
\[{\bf{W}}(t): = {\int_{\rm{0}}^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{B}}{{\bf{B}}^T}\left( {{e^{{\bf{A}}\left( {t - \tau } \right)}}} \right)} ^T}d\tau \]為 nonsingular。
=====================

3. 關於狀態轉移矩陣相關議題，請參考BLOG相關文章
[線性系統] 動態方程式的求解(1) - LTI state equation
[線性系統] 動態方程式的求解(2) - LTV state equation- Homogeneous solution
[線性系統] 動態方程式的求解(3) - LTV state equation- Total Solution

[線性系統] 漸進穩定度 與 Lyapunov Theorem

這次要介紹如何 透過 Lyapunov Theorem 來檢驗線性系統 ${\bf{\dot x}} = {\bf{Ax}}$ 的漸進穩定度 (Asymptotic Stability)。關於非線性系統的漸進穩定度讀者可參考下列兩篇文章：

•  [系統理論] 離散時間系統的穩定度理論 (0) - 先備概念
• [系統理論] 離散時間系統的穩定度理論 (1) - Lyapunov Stability Theory

現在回憶我們先前提過控制系統的兩種絕對穩定度：BIBO穩定 與 漸進穩定度。
概念上 BIBO穩定為插上電源看看系統會不會壞掉，漸進穩定則是測試拔掉電源看看系統會不會停止。

Lyapunov Energy Ideas
一般而言，Lyapunov 觀點是透過能量的角度看系統穩定度。也就是說考慮系統狀態 ${\bf{x}}\left( t \right) $，那麼
 \[
{\bf{x}}\left( t \right) \to 0 \Leftrightarrow {{\bf{x}}^T}\left( t \right){\bf{x}}\left( t \right) \to 0
\] 注意到上述 ${{\bf{x}}^T}\left( t \right){\bf{x}}\left( t \right)$ 可看成能量。那麼為了達成上式，我們可以透過 能量對時間的變化率 (系統狀態能量對時間微分) 若為負值，則表示能量在逐漸溢散(decaying energy)，亦即可透過
\[
\frac{d}{dt} {{\bf{x}}^T}\left( t \right){\bf{x}}\left( t \right) <0
\] 達成 ${\bf{x}}\left( t \right) \to 0 \Leftrightarrow {{\bf{x}}^T}\left( t \right){\bf{x}}\left( t \right) \to 0$

注意：這邊我們說 ${\bf{A}}$ 矩陣為穩定若下面條件成立：
對 ${\bf{A}}$ 的所有 eigenvalue 有負實部。

現在我們看一個例子來展示 Lyapunov Energy Idea，

Example
考慮
\[{\bf{\dot x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
{ - 1}&0\\
2&{ - 3}
\end{array}} \right]{\bf{x}}\left( t \right)
\]現試著找出系統能量是否 decaying?

Solution
注意到此系統 $\bf A$ 矩陣 為常數下三角矩陣，其 eigenvalue 為 $-1, -3$ 由穩定度定理可知系統為穩定系統。現在我們看看 Lyapunov Energy Idea 是否也可以幫助我們判別系統穩定度。

首先觀察系統狀態能量的微分
\[\begin{array}{l}
\frac{d}{{dt}}{{\bf{x}}^T}\left( t \right){\bf{x}}\left( t \right) = \frac{d}{{dt}}\left( {\left[ {\begin{array}{*{20}{c}}
{{x_1}\left( t \right)}&{{x_2}\left( t \right)}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_1}\left( t \right)}\\
{{x_2}\left( t \right)}
\end{array}} \right]} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \frac{d}{{dt}}\left( {{x_1}^2\left( t \right) + {x_2}^2\left( t \right)} \right) = 2{x_1}\left( t \right){{\dot x}_1}\left( t \right) + 2{x_2}\left( t \right){{\dot x}_2}\left( t \right)
\end{array}
\]又因為
\[{\bf{\dot x}}\left( t \right) = \left[ {\begin{array}{*{20}{c}}
{ - 1}&0\\
2&{ - 3}
\end{array}} \right]{\bf{x}}\left( t \right) \Leftrightarrow \left[ {\begin{array}{*{20}{c}}
{{{\dot x}_1}\left( t \right)}\\
{{{\dot x}_2}\left( t \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{ - {x_1}\left( t \right)}\\
{2{x_1}\left( t \right) - 3{x_2}\left( t \right)}
\end{array}} \right]
\]故我們可得
\[\begin{array}{l}
\frac{d}{{dt}}{{\bf{x}}^T}\left( t \right){\bf{x}}\left( t \right) = 2{x_1}\left( t \right){{\dot x}_1}\left( t \right) + 2{x_2}\left( t \right){{\dot x}_2}\left( t \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = 2{x_1}\left( t \right)\left( { - {x_1}\left( t \right)} \right) + 2{x_2}\left( t \right)\left( {2{x_1}\left( t \right) - 3{x_2}\left( t \right)} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} =  - 2{x_1}^2\left( t \right) + 4{x_2}\left( t \right){x_1}\left( t \right) - 6{x_2}^2\left( t \right) \ \ \ \ (*)
\end{array}
\]現在若上式 $<0$ 則我們由 Lyapunov Energy Ideas 即可斷定系統狀態 ${\bf{x}}\left( t \right) \rightarrow 0$。故我們進一步改寫 $(*)$ 成矩陣形式：
\[\frac{d}{{dt}}{{\bf{x}}^T}\left( t \right){\bf{x}}\left( t \right) = {{\bf{x}}^T}\left( t \right)\underbrace {\left[ {\begin{array}{*{20}{c}}
{ - 2}&2\\
2&{ - 6}
\end{array}} \right]}_Q{\bf{x}}\left( t \right)
\]注意到上述矩陣 $Q$ 為 對稱 負定矩陣(negative definite ) 因為 由對稱負定矩陣定義 ： $-Q$ 必須為正定矩陣。由於
\[-Q = \left[ {\begin{array}{*{20}{c}}
2&{ - 2}\\
{ - 2}&6
\end{array}} \right]
\]其對應的
1st Leading principal minor: $= 2 >0$,
2nd Leading principal minor: $= 2 \times 6 - (-2) \times (-2) = 8>0$，故 $Q$ 為負定矩陣。且我們知道此系統能量會溢散。亦即 Lyapunov Energy Ideas 確實可以幫助我們判斷系統穩定度。

現在看下面這個定理:
=======================
Theorem: Lyapunov Theorem
 ${\bf{A}}$ 矩陣的 所有 eigenvalue 有 負實部 ( ${\bf{A}}$ 矩陣 為穩定) 若且為若
對任意給定 正定對稱 (Positive definite symmetric) 矩陣 ${\bf {Q}}$，其 Lyapunov equation
\[
{{\bf{A}}^T}{\bf{P + PA = }} - {\bf{Q}}
\]有 唯一解 ${\bf {P}}$， 且此唯一解 ${\bf {P}}$ 為 正定對稱矩陣。
=======================

上述定理與 Lyapunov Energy 能量觀點可以整合

對 ${\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right)$，現在定義 Energy-like function
\[V\left( {\bf{x}} \right){\rm{ }}: = {{\bf{x}}^T}\left( t \right){\bf{Px}}\left( t \right)\]
其中 $\bf P$ 為 Lyapunov equation ${{\bf{A}}^T}{\bf{P}} + {\bf{PA}} =  - {\bf{Q}} $ 的解。則
\[
\frac{d}{dt} {\bf V(x(t))} <0
\]為系統漸進穩定度的判別條件。

[線性系統] LTI 系統的輸入輸出 BIBO 穩定度

這次要介紹線性系統的穩定度。一般而言在設計控制系統的時候第一步就是要檢驗系統是否穩定，如果不穩定則往往導致系統損毀。不可不慎。

一般而言控制系統穩定度可區分兩類

1. 絕對穩定度：
   指系統是否穩定的指標：一般而言有 BIBO 穩定 與 漸進穩定。
2. 相對穩定度：指系統穩定程度的指標: 一般而言由 pole location，Phase Margin, Gain Margin 決定

而一般穩定度的判別方法也有兩種

1. Routh-Hurwitz criterion
   只適用於線性系統，有興趣讀者可自行參閱任何一本自動控制教科書都會有詳細介紹。
2. Lyapunov energy approach
   對線性/非線性系統皆適用。

Comment:

給不關心理論的讀者：事實上，在實用面上，大多時候我們可以直接使用 MATLAB 等套裝軟體直接求解 eigenvalue 並且判斷是否落在 s-plane 的左半面即可 (如果落在左半面不含虛軸，我們稱此系統 "穩定" )。

考慮一個 SISO LTI 系統描述如下：
\[
y(t) = \int_0^t g(t-\tau)u(\tau)d \tau = \int_0^t g(\tau) u(t- \tau)d\tau
\]其中 $g(t)$ 為系統脈衝響應(impulse response)

現在我們給出下面定義
====================
Definition: Bounded function
一個輸入函數 $u(t)$ 稱作有界 (bounded) 如果下列條件成立：
若存在一個夠大的常數 $u_M$ 使得
\[
|u(t)| \le u_M < \infty, \; \forall t \ge 0
\]====================

有了有界函數的定義，我們可以定義何謂 BIBO 穩定
====================
Definition: BIBO stability
一個系統被稱為 BIBO stable (Bounded input bounded output stable) 若下列條件成立：
對任意有界輸入，系統都產生有界輸出。則此系統為 BIBO 穩定。
====================

Comment:
1. BIBO 穩定度 只定義在 zero-state response 且系統假設為鬆弛系統亦即初始狀態為零。

2. 上述 BIBO 定義 陳述等價如下：
考慮系統 輸入為 $u(t)$，輸出為 $y(t)$，現若存在 $M, N >0$ 使得
\[
|u(t)| \le M < \infty \Rightarrow |y(t)| \le N < \infty
\]則系統稱為 BIBO穩定。


現在我們看個重要的定理:

=====================
Theorem: BIBO stability criterion of LTI system 
考慮 SISO LTI 系統
\[
y(t) = \int_0^t g(t-\tau)u(\tau)d \tau = \int_0^t g(\tau) u(t- \tau)d\tau
\]為 BIBO stable 若且為若 存在一個夠大的常數 $M$ 使得
\[
\int_0^\infty |g(t)|dt \le M < \infty
\]=====================
Proof
我們首先證明 $\Leftarrow$，亦即
假設
\[
\int_0^\infty |g(t)|dt \le M < \infty
\] 欲證 LTI 系統 為 BIBO穩定。

故由BIBO定義，給定 $u(t)$ 滿足 $|u(t)| \le u_M < \infty$，要證明 $|y(t)| < \infty$。

由於
\[\begin{array}{l}
\left| {y(t)} \right| = \left| {\int_0^t g (\tau )u(t - \tau )d\tau } \right|\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} \le \int_0^t {\left| {g(\tau )u(t - \tau )} \right|} d\tau  \le {u_M}\int_0^t {\left| {g(\tau )} \right|} d\tau
\end{array}
\]現在由假設 $\int_0^\infty |g(t)|dt \le M < \infty$，我們可得
\[\left| {y(t)} \right| \le {u_M}\int_0^t {\left| {g(\tau )} \right|} d\tau  \le {u_M}\int_0^\infty {\left| {g(\tau )} \right|} d\tau  \le {u_M}M < \infty
\] 故得證。

接著我們證明 $\Rightarrow$。亦即需要證明
若 SISO LTI 系統
\[
y(t) = \int_0^t g(t-\tau)u(\tau)d \tau = \int_0^t g(\tau) u(t- \tau)d\tau
\]為 BIBO stable，則 存在一個夠大的常數 $M$ 使得
\[
\int_0^\infty |g(t)|dt \le M < \infty
\]
現在我們取非，利用反證法 改證
對所有常數 $M$，\[
\int_0^\infty |g(t)|dt >M
\]，則 SISO LTI 系統
\[
y(t) = \int_0^t g(t-\tau)u(\tau)d \tau = \int_0^t g(\tau) u(t- \tau)d\tau
\]不為 BIBO stable (存在一組輸入 $u(t)$ 使得 輸出 $y(t) = \infty$)

首先注意到由於對所有 $M$，
\[
\int_0^\infty |g(t)|dt >M
\]，此陳述等價於 存在一個 $t_1$ 使得
\[
\int_0^{t_1} |g(t)|dt = \infty
\]接著我們要證明 LTI 系統不為 BIBO stable，故現在建構一組輸入 $u(t)$ 滿足下式
\[u\left( t_1 - \tau \right) := {\rm sgn}(g(\tau) )= \left\{ \begin{array}{l}
1,\begin{array}{*{20}{c}}
{}
\end{array}{\rm if}\begin{array}{*{20}{c}}
{}
\end{array}g\left( \tau  \right) \ge 0\\
 - 1,\begin{array}{*{20}{c}}
{}
\end{array}{\rm if}\begin{array}{*{20}{c}}
{}
\end{array}g\left( \tau  \right) < 0
\end{array} \right.
\]則其對應的輸出
\[y({t_1}) = \int_0^{{t_1}} g (\tau )u({t_1} - \tau )d\tau  = \int_0^{{t_1}} {\left| {g(\tau )} \right|} d\tau  = \infty \]故此系統不為 BIBO stable。至此得證。$\square$

Comment:
1. 上述結果等價如下：
若系統為 LTI 系統，則 系統 BIBO 穩定的充分必要條件為
其系統轉移函數對應的極點 $p_i$ 具有負實部；亦即
\[
Re[p_i] <0, i=1,2,...,n
\]

2. 若考慮系統為 SISO LTI 系統，且轉移函數 沒有發生極零點對消 (pole-zero cancelation)，則系統特性根 poles = Eigenvalues ；亦即 BIBO穩定 = 漸進穩定度。

3. 若轉移函數 發生極零點對消 (此時 eigenvalue 的數目會多於 pole 的數目 (因為 pole 被 zero消了!!) ) 且 極零點對消在左半面，BIBO = 漸進穩定；若對消在右半面 或者 虛軸上，則 BIBO 不等於 漸進穩定。亦即 BIBO 穩定不一定為 漸進穩定。

4. 若 $\dot x = Ax + Bu;\;\; y=Cx$ 則取拉式轉換我們有
\[\begin{array}{l}
\left\{ \begin{array}{l}
X\left( s \right) = {\left( {sI - A} \right)^{ - 1}}x\left( 0 \right) + {\left( {sI - A} \right)^{ - 1}}BU\left( s \right)\\
Y\left( s \right) = CX\left( x \right)
\end{array} \right.\\
 \Rightarrow Y\left( s \right) = \underbrace {\left[ {C{{\left( {sI - A} \right)}^{ - 1}}} \right]}_{{\rm{Free}}\begin{array}{*{20}{c}}
{}
\end{array}{\rm{Response}}}x\left( 0 \right) + \underbrace {\left[ {C{{\left( {sI - A} \right)}^{ - 1}}B} \right]}_{{\rm{Force}}\begin{array}{*{20}{c}}
{}
\end{array}{\rm{Response}}}U\left( s \right)
\end{array}\]

我們以下用個例子來看看剛剛討論的 pole-zero cancellation 現象 與 BIBO/Asymptotic stability。

Example: (BIBO stability does NOT imply Asymptotic stability)
考慮系統
\[\left\{ \begin{array}{l}
\dot x = Ax + Bu = \left[ {\begin{array}{*{20}{c}}
0&1\\
2&{ - 1}
\end{array}} \right]x + \left[ {\begin{array}{*{20}{c}}
0\\
1
\end{array}} \right]u\\
y = Cx = \left[ {\begin{array}{*{20}{c}}
{ - 2}&2
\end{array}} \right]x
\end{array} \right.
\](a) 試求 $A$ 矩陣 eigenvalue。此系統是否為 asymptotic stable?
(b) 試求此系統轉移函數?
(c) 試求此系統 unit step response，試問其對應的系統輸出 是否為 BIBO?

Solution: (a)
首先求 $A$ 矩陣的 eigenvalue：計算特性方程 $\det(sI-A) =0$ 可得
\[\begin{array}{l}
\det \left( {sI - A} \right) = 0\\
 \Rightarrow \det \left( {\left[ {\begin{array}{*{20}{c}}
s&0\\
0&s
\end{array}} \right] - \left[ {\begin{array}{*{20}{c}}
0&1\\
2&{ - 1}
\end{array}} \right]} \right) = \det \left[ {\begin{array}{*{20}{c}}
s&{ - 1}\\
{ - 2}&{s + 1}
\end{array}} \right] = s\left( {s + 1} \right) - 2 = 0\\
 \Rightarrow {s^2} + s - 2 = 0 \Rightarrow \left\{ \begin{array}{l}
s =  - 2\\
s = 1
\end{array} \right.
\end{array}
\]注意到此系統有一個 eigenvalue 為 $s=1$ (落在 s-plane 右半面)，故系統不為漸進穩定。

(b):
接著我們計算轉移函數
\[\begin{array}{l}
G\left( s \right) = C{\left( {sI - A} \right)^{ - 1}}B\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \left[ {\begin{array}{*{20}{c}}
{ - 2}&2
\end{array}} \right]{\left[ {\begin{array}{*{20}{c}}
s&{ - 1}\\
{ - 2}&{s + 1}
\end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}}
0\\
1
\end{array}} \right]\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \frac{{2s - 2}}{{{s^2} + s - 2}} = \frac{{2\left( {s - 1} \right)}}{{\left( {s + 2} \right)\left( {s - 1} \right)}} = \frac{2}{{s + 2}}
\end{array}\]注意到此時發生 pole-zero cancellation ($s=1$ 對消在右半面) !!

(c):
接著我們計算 unit-step response (注意到 unit-step 為 Bounded-input)
\[\begin{array}{l}
Y\left( s \right) = G\left( s \right)U\left( s \right) = \left( {\frac{2}{{s + 2}}} \right)\left( {\frac{1}{s}} \right) = \frac{2}{{s\left( {s + 2} \right)}} = \frac{1}{s} + \frac{{ - 1}}{{s + 2}}\\
 \Rightarrow y\left( t \right) = 1 - {e^{ - 2t}}
\end{array}\]可發現系統輸出 $|y(t)| \le 1$ 亦即為 Bounded output。故系統為 BIBO 系統 under unit-step input。

[線性系統] 動態方程式的求解(1) - LTI state equation

延續上篇，這次我們要介紹 線性非時變系統 (Linear Time-Invariant (LTI) System) 的求解。

考慮 LTI 動態系統的 狀態空間表示：
\[\left\{ {\begin{array}{*{20}{l}}
{{\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)}\\
{{\bf{y}}\left( t \right) = {\bf{Cx}}\left( t \right) + {\bf{Du}}\left( t \right)}
\end{array}} \right.
\]其中 $\bf{A}(\cdot), \bf{B}(\cdot), \bf{C}(\cdot),$ 與 $\bf{E}(\cdot)$ 為 $n \times n, n \times p, q \times n,$ 與 $q \times p$  常數矩陣。


我們的目標： 求解 $\bf{x}(t)$。

在求解之前我們需要一些 exponential function ${e^{  {\bf{A}}t}}$的 FACTs:
首先回憶若令 $\bf A = a$ 亦即不再是矩陣而是一個常數 $a$，則 $e^{at}$ 具有 Taylor series 如下
\[{e^{at}}: = 1 + at + \frac{{{a^2}{t^2}}}{{2!}} + ... + \frac{{{a^n}{t^n}}}{{n!}} + ...\]故若現在讓 $a$ 變回矩陣 $\bf A$ 則我們有
\[{e^{{\bf{A}}t}}: = {\bf{I}} + {\bf{A}}t + \frac{1}{{2!}}{{\bf{A}}^2}{t^2} + ... + \frac{1}{{n!}}{{\bf{A}}^n}{t^n} + ...
\]那麼下面幾個性質，讀者可以使用上述定義直接驗證。
=====================
FACT 1: Inverse property
\[
{e^{  {\bf{A}}t}}{e^{ - {\bf{A}}t}} = \bf{I}.
\]FACT 2: Identity matrix
\[
e^{\bf{0}} = \bf{I}.
\]FACT 3: Derivative property
\[\frac{d}{{dt}}{e^{{\bf{A}}t}} = {\bf{A}}{e^{{\bf{A}}t}} = {e^{{\bf{A}}t}}{\bf{A}}.
\]=====================

有了上述的 FACT，我們現在可以開始求解 LTI 動態系統，不過求解之前我們先試著計算看看 $e^{{\bf A}{t}}$ 該怎麼計算  (這邊我們採用定義求解法，之後會用更直接的方法求解)。

-----------
Example 1
考慮矩陣
\[{\bf{A}}: = \left[ {\begin{array}{*{20}{c}}
0&1\\
0&0
\end{array}} \right]\]試求 ${e^{{\bf{A}}t}} = ?$
-----------

Solution
由 ${e^{{\bf{A}}t}} $ 定義可知
\[\begin{array}{l}
{e^{{\bf{A}}t}}: = {\bf{I}} + {\bf{A}}t + \frac{1}{{2!}}{{\bf{A}}^2}{t^2} + ... + \frac{1}{{n!}}{{\bf{A}}^n}{t^n} + ...\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \left[ {\begin{array}{*{20}{c}}
1&0\\
0&1
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
0&1\\
0&0
\end{array}} \right]t + \frac{1}{{2!}}\underbrace {{{\left[ {\begin{array}{*{20}{c}}
0&1\\
0&0
\end{array}} \right]}^2}}_{ = {\bf{0}}}{t^2} + ...\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \left[ {\begin{array}{*{20}{c}}
1&t\\
0&1
\end{array}} \right]. \ \ \ \ \ \ \square
\end{array}\]

現在我們可以開始推導完整的 LTI 系統的狀態方程解，首先 考慮狀態方程
\[
{{\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)}
\] 對上式左右同乘 ${e^{{-\bf{A}}t}}$ 可得
\[\begin{array}{l}
{e^{ - {\bf{A}}t}}{\bf{\dot x}}\left( t \right) = {e^{ - {\bf{A}}t}}{\bf{Ax}}\left( t \right) + {e^{ - {\bf{A}}t}}{\bf{Bu}}\left( t \right)\\
 \Rightarrow {e^{ - {\bf{A}}t}}{\bf{\dot x}}\left( t \right) - {e^{ - {\bf{A}}t}}{\bf{Ax}}\left( t \right) = {e^{ - {\bf{A}}t}}{\bf{Bu}}\left( t \right)
\end{array}
\]觀察上式，我們可得
\[\frac{d}{{dt}}\left( {{e^{ - {\bf{A}}t}}{\bf{x}}\left( t \right)} \right) = {e^{ - {\bf{A}}t}}{\bf{Bu}}\left( t \right)
\]現在對兩邊同取積分從 $0$ 到 $t$ ，可得
\[\begin{array}{l}
\int_0^t {d\left( {{e^{ - {\bf{A}}t}}{\bf{x}}\left( t \right)} \right)}  = \int_0^t {{e^{ - {\bf{A}}\tau }}{\bf{Bu}}\left( \tau  \right)d\tau } \\
 \Rightarrow {e^{ - {\bf{A}}t}}{\bf{x}}\left( t \right) - {\bf{x}}\left( 0 \right) = \int_0^t {{e^{ - {\bf{A}}\tau }}{\bf{Bu}}\left( \tau  \right)d\tau } \\
 \Rightarrow {\bf{x}}\left( t \right) = {e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + {e^{{\bf{A}}t}}\int_0^t {{e^{ - {\bf{A}}\tau }}{\bf{Bu}}\left( \tau  \right)d\tau } \\
 \Rightarrow {\bf{x}}\left( t \right) = {e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + \int_0^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } \ \ \ \ (*)
\end{array}
\] 現在我們回頭檢驗上式 $(*)$ 確實為 我們狀態方程的解，我們首先檢驗其確實符合初始條件，亦即當 $t=0 $時，${\bf{x}}\left( t \right) = {\bf{x}}\left( 0 \right)$

現令 $t = 0$，我們可得
\[\begin{array}{l}
{\bf{x}}\left( t \right) = {e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + \int_0^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } \\
 \Rightarrow {\bf{x}}\left( 0 \right) = {e^{{\bf{A}}0}}{\bf{x}}\left( 0 \right) = {\bf{Ix}}\left( 0 \right) = {\bf{x}}\left( 0 \right)
\end{array}
\]故滿足初始條件，現在我們來檢驗 $(*)$ 滿足 狀態方程 ${{\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)}$，
故現在對 $(*)$ 微分可得
\[\frac{d}{{dt}}{\bf{x}}\left( t \right) = \frac{d}{{dt}}\left[ {{e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + \int_0^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } } \right]
\]由 Fundamental Theorem of Calculus，我們可知
\[\frac{\partial }{{\partial t}}\int_{{t_0}}^t {f\left( {t,\tau } \right)d\tau }  = \left. {f\left( {t,\tau } \right)} \right|_{\tau  = t}^{} + \int_{{t_0}}^t {\left( {\frac{\partial }{{\partial t}}f\left( {t,\tau } \right)} \right)d\tau }
\]故
\[\begin{array}{l}
\frac{d}{{dt}}{\bf{x}}\left( t \right) = \frac{d}{{dt}}\left[ {{e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + \int_0^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } } \right]\\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{A}}{e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + \frac{d}{{dt}}\left[ {\int_0^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } } \right]\\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{A}}{e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + {e^{{\bf{A}}\left( {t - t} \right)}}{\bf{Bu}}\left( t \right) + \int_0^t {\left( {\frac{d}{{dt}}{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)} \right)d\tau } \\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{A}}{e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + {\bf{Bu}}\left( t \right) + \int_0^t {{\bf{A}}{e^{{\bf{A}}\left( t \right)}}{e^{{\bf{A}}\left( { - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } \\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{A}}\underbrace {\left[ {{e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + \int_0^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } } \right]}_{{\rm{ = }}{\bf{x}}\left( t \right)} + {\bf{Bu}}\left( t \right)\\
 \Rightarrow {\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)
\end{array}
\] 故得證 $(*)$ 即稱為 LTI 狀態方程的解。

我們將上述結果紀錄在下方：
考慮 LTI 動態系統的 狀態方程：
\[{\bf{\dot x}}\left( t \right) = {\bf{Ax}}\left( t \right) + {\bf{Bu}}\left( t \right)
\]其解為
\[{\bf{x}}\left( t \right) = {e^{{\bf{A}}t}}{\bf{x}}\left( 0 \right) + \int_0^t {{e^{{\bf{A}}\left( {t - \tau } \right)}}{\bf{Bu}}\left( \tau  \right)d\tau } \]


那麼現在問題變成，如果我們拓展LTI系統到線性時變 (Linear Time Varying, LTV) 系統，則上述方法是否仍然可行?
答案是 "否定" 的，故我們需要令求其他方法來幫助我們，在下一篇文章我們將會介紹線性時變 (LTV) 系統的狀態方程該如何求解:
[線性系統] 動態方程式的求解(2) - LTV state equation, Fundamental Matrix, and State Transition Matrix

[線性系統] 動態方程式的求解(0) - Review 1st ODE, DE, & $e^{At}$

在求解動態方程是之前，我們先回顧一下基本 一階常微分方程的求解：

Example 1
試求解
\[
\dot {x} = a x, \ x(t_0) = x_0
\]其中 $a$ 為常數。
Solution
利用變數分離法：
\[\begin{array}{*{20}{l}}
{\frac{d}{{dt}}x\left( t \right) = ax\left( t \right)}\\
{ \Rightarrow \int_{{t_0}}^t {\frac{1}{{x\left( \tau  \right)}}dx\left( \tau  \right)}  = \int_{{t_0}}^t {ad\tau } }\\
{ \Rightarrow \left. {\ln \left( {x\left( \tau  \right)} \right)} \right|_{{t_0}}^t = a\left( {t - {t_0}} \right)}\\
{ \Rightarrow \ln \left( {x\left( t \right)} \right) - \ln \left( {x\left( {{t_0}} \right)} \right) = a\left( {t - {t_0}} \right)}\\
{ \Rightarrow \ln \left( {\frac{{x\left( t \right)}}{{{x_0}}}} \right) = a\left( {t - {t_0}} \right)}\\
{ \Rightarrow x\left( t \right) = {x_0}{e^{a\left( {t - {t_0}} \right)}}}. \ \ \ \ \square
\end{array}
\]

Example 2
試求解
\[
\dot {x} = a x + bu, \ x(t_0) = x_0
\]其中 $a, b \neq 0$ 為常數。
Solution
上式為標準一階常微分方程，一般而言，微分方程的解可分為兩個部分

1. 自由響應 (free response)又稱 零輸入響應(zero input response) ；亦即令 $u(t) = 0$ 所求得的解
2. 外力響應 (forced response) 又稱 零狀態響應 (zero state response)；亦即令 $x(t_0) = 0$所求得的解

整個解 $x(t) $  = 自由響應 + 外力響應；

故首先求解 自由響應  ，令 $u=0$，我們得到
\[
\dot {x} = a x , \ x(t_0) = x_0
\]
由 Example 1，我們知道對於 $\dot{x} = ax$ 的微分方程，其對應的解為 $x(t) = x_0e^{a(t-t_0)}$。

現在我們求解外力響應：令 $x(t_0) = 0$
故我們可以猜對於 $\dot {x} = a x + bu$ ，其解具有下列型態：
\[
 x(t) = e^{at} F(t) \ \ \ \ (*)
\] 其中 函數 $F(t)$ 待定。
現在將上式 待入微分方程 $\dot{x} = ax + bu$ 之中我們得到
\[\begin{array}{l}
 a{e^{at}}F(t) + {e^{at}}\dot F(t) = a{e^{at}}F(t) + bu\\
 \Rightarrow {e^{at}}\dot F(t) = bu
\end{array}\] 現在我們可求解 $F(t)$ 如下
\[\begin{array}{l}
\frac{{dF(t)}}{{dt}} = {e^{ - at}}bu\left( t \right)\\
 \Rightarrow \int_{{t_0}}^t {dF(\tau )}  = \int_{{t_0}}^t {{e^{ - a\tau }}bu\left( \tau  \right)d\tau } \\
 \Rightarrow F(t) - F\left( {{t_0}} \right) = \int_{{t_0}}^t {{e^{ - a\tau }}bu\left( \tau  \right)d\tau } \\
 \Rightarrow F(t) = \int_{{t_0}}^t {{e^{ - a\tau }}bu\left( \tau  \right)d\tau }  + F\left( {{t_0}} \right) \ \ \ \ (\star)
\end{array}
\] 注意到 $x(t_0)=0$，故由 $(*)$ 我們亦可知
\[
 x(t) = e^{at} F(t) \Rightarrow  x(t_0) = e^{at_0} F(t_0) \Rightarrow 0 = e^{at_0} F(t_0)
\]由於 $e^{a t_0} \neq 0$，故 $ F(t_0) = 0$，所以 $(\star)$ 變成：
\[
F(t) = \int_{{t_0}}^t {{e^{ - a\tau }}bu\left( \tau  \right)d\tau }
\] 又因為我們原本由 $(*)$ 定義為 $ x(t) = e^{at} F(t)$，故外力響應 為
\[
x(t) =\int_{{t_0}}^t {{e^{a\left( {t - \tau } \right)}}bu\left( \tau  \right)d\tau }
\]

現在將 自由響應 與 外力響應合併起來，即可求得整個 一階常微分方程的總解

\[ x(t) =  x_0e^{a(t-t_0)} + \int_{{t_0}}^t {{e^{a\left( {t - \tau } \right)}}bu\left( \tau  \right)d\tau }.  \ \ \ \ \square \]

Example 3
試求解 下列一階 時變(time-varying) 微分方程
\[
\dot{x} = a(t) x(t), \ x(t_0) =x_0
\]
Solution
同 Example 1，我們可將其改寫並利用變數分離法求解
\[\begin{array}{l} \frac{d}{{dt}}x\left( t \right) = a(t)x(t)\\ \Rightarrow \frac{1}{{x(t)}}dx\left( t \right) = a(t)dt\\ \Rightarrow \int_{{t_0}}^t {\frac{1}{{x(\tau )}}dx\left( \tau \right)} = \int_{{t_0}}^t {a(\tau )d\tau } \\ \Rightarrow \ln \left. {x(\tau )} \right|_{{t_0}}^t = \int_{{t_0}}^t {a(\tau )d\tau } \\ \Rightarrow \ln x(t) - \ln x({t_0}) = \int_{{t_0}}^t {a(\tau )d\tau } \\ \Rightarrow \ln \frac{{x(t)}}{{x({t_0})}} = \int_{{t_0}}^t {a(\tau )d\tau } \\ \Rightarrow x(t) = {x_0}{e^{\int_{{t_0}}^t {a(\tau )d\tau } }} \end{array}
\]

有了上述的三個例子，我們現在可以把目標轉向 矩陣形式：亦即 $n$ 維 動態方程的求解：

考慮 $n$ 維 線性 時變 動態方程 (Linear Time-varying dynamical equation)
\[
\left\{ \begin{array}{l}
{\bf{\dot x}}\left( t \right){\bf{ = A}}\left( t \right){\bf{x}}\left( t \right){\bf{ + B}}\left( t \right){\bf{u}}\left( t \right)\\
{\bf{y}}\left( t \right){\bf{ = C}}\left( t \right){\bf{x}}\left( t \right){\bf{ + D}}\left( t \right){\bf{u}}\left( t \right)
\end{array} \right.
\] 其中 $\bf{A}(\cdot), \bf{B}(\cdot), \bf{C}(\cdot),$ 與 $\bf{E}(\cdot)$ 為 $n \times n, n \times p, q \times n,$ 與 $q \times p$  在時刻 $t$ 定義在 $(-\infty, \infty)$ 的 實數 連續函數矩陣

上式的求解需要一些數學工具的幫忙，我們首先介紹特殊矩陣函數 $e^{\bf{A}(t)}$：

首先回憶指數函數 $e^{\lambda t}$ 可由 Taylor Series 無窮級數展開：
\[
e^{\lambda t} = 1 + \lambda t + \frac{\lambda^2 t}{2!} + ... + \frac{\lambda^n t^n}{n!}+...
\]同樣的若現在考慮常數矩陣 ${\bf{A}}$，由 Caley-Hamilton Theorem ，我們亦可將 ${e^{{\bf{A}}(t)}}$ 透過類似 Taylor Series 無窮級數展開將其表為：
\[
{e^{{\bf{A}}{t}}} = {\bf{I}} + t{\bf{A}} + \frac{t}{{2!}}{{\bf{A}}^2} + ... + \frac{{{t^n}}}{{n!}}{{\bf{A}}^n} + ... = \sum\limits_{k = 0}^\infty  {\frac{{{t^k}}}{{k!}}{{\bf{A}}^k}}
\]

現在先看幾個性質：
1.  $e^{\bf{0}} = \bf {I} $。
(Proof: 由上述無窮級數展開令 $\bf{A} = 0$)

2. ${e^{{\bf{A}}({t_1} + {t_2})}} = {e^{{\bf{A}}({t_1})}}{e^{{\bf{A}}({t_2})}}$
 (Proof: 由上述無窮級數展開令 $t = t_1 + t_2$ 即可得證)

3.  ${\left[ {{e^{{\bf{A}}t}}} \right]^{ - 1}} = {e^{ - {\bf{A}}t}}$
Proof:
由上述性質2，令 $t_2 = -t_1$，我們可得
\[{e^{{\bf{A}}({t_1} - {t_1})}} = {e^{{\bf{A}}({t_1})}}{e^{ - {\bf{A}}({t_1})}} = {e^{{\bf{A}} \cdot 0}} = {e^{\bf{0}}} = {\bf{I}}
\] 由反矩陣定義，可知 ${\left[ {{e^{{\bf{A}}t}}} \right]^{ - 1}} = {e^{ - {\bf{A}}t}}$。

4. $\frac{d}{{dt}}{e^{{\bf{A}}t}} = {\bf{A}}{e^{{\bf{A}}t}} = {e^{{\bf{A}}t}}{\bf{A}} $

有了上述性質之後我們便可以開始討論狀態方程求解
[線性系統] 動態方程式的求解(1) - LTI state equation

[最佳控制] 離散時間 穩態 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)

[數學分析] 淺談各種基本範數 (Norm)

這次要介紹的是數學上一個重要的概念：

Norm：一般翻譯成範數
(在英語中 norm 有規範的意思，比如我們說normalization就是把某種東西/物品/事件 做 正規化，也就是加上規範使其正常化)，不過個人認為其實翻譯成 範數 也是看不懂的...這邊建議把 Norm 想成長度就好 (事實上norm是長度的抽象推廣)，

也許讀者會認為好端端的長度不用，為何又要發明一個 norm 來自討苦吃?? 既抽象又艱澀。

事實上想法是這樣的：
比如說現在想要比較兩個數字 $3$ , $5$ 之間的大小，則我們可以馬上知道 $ 3 < 5 $；同樣的，如果再考慮小數與無理數如 $1.8753$ 與 $\pi$，我們仍然可以比較大小 $1.8753 < \pi = 3.1415...$ 故可以發現我們有辦法對 "純量" 做明確的比大小，WHY? 因為前述例子中 $3$, $5$, $1.8753$ or $\pi$ 其各自的大小有辦法被 "measure "!

但是如果是現在考慮的是一組數字 我們如何去measure 其大小呢?? 比如說
\[x:=[1, -2, 0.1, 0 ]^T
\]上式的大小該是多少? 是 $1$? $-2$? $0.1$???
再者如果更過分一點，我們考慮一個矩陣
\[A = \left[ {\begin{array}{*{20}{c}}
1&2\\
3&4
\end{array}} \right]
\]，想要知道這個矩陣的大小又該怎麼辦?? 是 $1$ ? $2$ 還是 $4$ ?..其實現階段我們說不清楚。

也正是如此，可以發現我們確實需要新的 "長度" 的定義來幫助我們如何去 measure 矩陣/向量/甚至是函數的大小。

故此，我們首先定義甚麼是Norm，(也就是把 "長度" or "大小" 的本質抽離出來)

==================
Definition: Norm
考慮 $V$ 為一個向量空間(Vector space)，則我們說  Norm 為一個函數 $||\cdot|| : V \rightarrow \mathbb{R}$ 且滿足下列性質：

(a) $||v|| \geq 0$, $||v||=0 \Leftrightarrow v=0$

(b) 對所有的 $\alpha \in \mathbb{R}, v \in V$, $||\alpha v|| = |\alpha| \cdot ||v||$

(c) 對任意 $v_1, v_2 \in V$ ， $||v_1 + v_2 || \le ||v_1|| + ||v_2||$ (三角不等式)
==================
Comments:
1. 觀察上述定義確實可以發現norm是長度的推廣。比如說性質(a) (positivity) 要求norm非負，可以發現直觀上我們的長度亦為非負；接著性質 (b) , (c) 亦可由原本長度定義中觀察而得。

2. 注意到上述 norm是定義在 向量空間 (vector space)，這是一個相當抽象的概念，一般視為是具有加法與乘法 ( $+, \cdot$) 的無窮維空間。簡單說可以把它想成是 $\mathbb{R}^n$的推廣，下面我們介紹的norm所定義的空間都為一個Vector space ( $\mathbb{R}^n$ space, $L^p$ space, 或者所有連續函數所形成的空間...)；如果對空間概念不熟的讀者可以暫時先略過此comment

==================
再來我們要考慮三大類常見的 norm  (也就是可以測量其大小但又滿足norm的3個基本性質)：
一類是有限維度 $\mathbb{R}^n$ 空間中的向量怎麼定義norm，
一類是如果拓展到無窮維度該怎麼定義norm?
最後則是如果考慮的是函數空間(function space) 該怎麼定義norm


首先是 $\mathbb{R}^n$ 空間 (有限維度空間) 的norm

考慮一個向量 $v \in \mathbb{R}^n$，亦即
\[v = \left[ {\begin{array}{*{20}{c}}
{{v_1}}\\
{{v_2}}\\
 \vdots \\
{{v_n}}
\end{array}} \right]
\]則我們可以定義下列各種不同的norm
\[\begin{array}{l}
{\left\| v \right\|_1}: = \sum\limits_{i = 1}^n {\left| {{v_i}} \right|}, \ \ \text{1-norm}\\
{\left\| v \right\|_2}: = \sqrt {\sum\limits_{i = 1}^n {{{\left| {{v_i}} \right|}^2}} }  = \sqrt {{v^T}v}, \ \ \text{2-norm}\\
{\left\| v \right\|_p}: = {\left( {\sum\limits_{i = 1}^n {{{\left| {{v_i}} \right|}^p}} } \right)^{\frac{1}{p}}}, \ \ \text{p-norm}\\
{\left\| v \right\|_\infty }: = \mathop {\max }\limits_i \left| {{v_i}} \right| = \max \{ |v_1|, |v_2|, ..., |v_n| \}, \ \ \text{Infinity norm}
\end{array}
\]

Comments:
1. 在 $\mathbb{R}^n$ 空間中，任意兩個norm的定義等價 (in convergence sense)。此稱作Two-norm theorem：我們將其紀錄如下： (注意如果已經離開有限維度空間，Two norm Theorem失效!!!)

Theorem: (Two-norm Theorem)
若 $||\cdot||_a$ 與 $||\cdot||_b$ 為 $\mathbb{R}^n$ 空間中的norm，則存在一個常數 $C < \infty$使得 對任意 $x \in \mathbb{R}^n$
\[
 ||x||_a \leq C ||x||_b
\]
Proof: omitted.

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

上述在有限維度空間的norm定義在一般情況仍然不足，因為如果我們考慮的是 一個sequence $\{x_k\}_{k=0}^{\infty}$，由於 sequence 有無窮多項 比如說 $\{x_k\}_{k=0}^{\infty}= 1/k $，此時上述的有限維度空間的norm定義無法使用。故我們進入了無窮維度的世界，此時無窮維度空間 的norm 可以仿照有限維度空間的定義如下
\[\begin{array}{*{20}{l}}
{{{\left\| x \right\|}_{{\ell _1}}}: = \sum\limits_{i = 1}^\infty  {\left| {{x_i}} \right|} }\\
{{{\left\| x \right\|}_{{\ell _2}}}: = \sqrt {\sum\limits_{i = 1}^\infty  {{{\left| {{x_i}} \right|}^2}} }  = \sqrt {{x^T}x} }\\
{{{\left\| x \right\|}_{{\ell _p}}}: = {{\left( {\sum\limits_{i = 1}^\infty  {{{\left| {{x_i}} \right|}^p}} } \right)}^{\frac{1}{p}}}}\\
{{{\left\| x \right\|}_{{\ell _\infty }}}: = \mathop {\sup }\limits_{i \ge 1} \left| {{x_i}} \right|}
\end{array}\]
Comment
注意到關於maximum norm這邊我們用 $\sup$ 取代 有限維度空間所採用的 $\max$，一般表示我們並不確定最大元素是否存在，故此時採用 $\sup$。

=======
現在如果我們考慮的是更複雜的情況，比如說我們考慮的是函數，則該怎麼對函數定義norm呢? (此時norm定義在函數空間 $L^p$-space)。考慮函數 $f : [a,b] \rightarrow \mathbb{R}$則我們有如下的norm
\[\begin{array}{l}
{\left\| {f\left( t \right)} \right\|_{{L_1}}}: = \int_a^b {\left| {f\left( t \right)} \right|dt} \\
{\left\| {f\left( t \right)} \right\|_{{L_2}}}: = {\left( {\int_a^b {{{\left| {f\left( t \right)} \right|}^2}dt} } \right)^{\frac{1}{2}}}\\
{\left\| {f\left( t \right)} \right\|_{{L_p}}}: = {\left( {\int_a^b {{{\left| {f\left( t \right)} \right|}^p}dt} } \right)^{\frac{1}{p}}}\\
{\left\| {f\left( t \right)} \right\|_{{L_\infty }}}: = \mathop {\sup }\limits_{t \in \left[ {a,b} \right]} \left| {f\left( t \right)} \right|
\end{array}
\]

Comment:
1. 上述三大類的norm皆滿足 norm本身的定義 (4個性質皆滿足，這需要證明但我們這邊不贅述)。
2. 對於連續函數 $f \in \mathcal{C} [[a,b], \mathbb{R}]$而言，其norm都選maximum norm：亦即 $||f|| := \max_{t \in [a,b]} |f(t)|$
3. 讀者也許會想問 $L_p$ 這個 $L$ 的意義為何，一般而言我們指 $f \in L_p$ 表示此函數 $f$ 為 Lebesgue integrable 且滿足上述 norm的條件。當然，如果有需要我們可以限制 $f \in R_p$ 其中 $R$ 表示的是函數 $f$ 為 Riemann integrable。在此不再贅述。


另外如果考慮的是矩陣的norm該怎麼定義??
因為我們手上只有對向量 與 對函數的norm，故方法就是先把矩陣變成向量，再算他的norm。我們稱此種算法得到的norm為 induced matrix norm。就是用一組向量把矩陣的norm 引(induce)出來：

Example:
考慮
\[\left[ {\begin{array}{*{20}{c}}
1&2\\
3&4
\end{array}} \right]\]那麼此矩陣的 norm 該怎麼計算?
我們可透過取向量 $[1 \; 0]^T$ 當作輸入，則
\[\left[ {\begin{array}{*{20}{c}}
1&2\\
3&4
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
1\\
0
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
1\\
3
\end{array}} \right]\]接著我們便可用 $\mathbb{R}^n$ 空間定義的向量來定義其 norm。


以下為正式的矩陣norm的定義

Definition: Induced matrix p-norm
令 $A = [a_{ij}] \in \mathbb{R}^{m \times n}$則此矩陣由 向量 $x$ 引出的 induced p-norm
\[
||A||_p := \sup_{x \neq 0}\frac{||Ax||_p}{||x||_p}
\]

注意到如果 $p=2$，則 induced matrix 2-nrom ，則我們有 $||A||_2$ 可由其最大的eigenvalue求得：
\[
||A||_2 := \sqrt{\lambda_{max} (A^T A)}
\]

Norm idea in System Theory (The Energy idea)
在系統理論中，我們另外亦常用 能量(Energy) 與 功率(power)來表示我們關心的物理量的大小：
考慮 $x(t)$ 為一維連續時間"訊號"，則定義在 $t\in [t_1,t_2]$ 之間 的 total energy E為
\[
E:=\int_{t_1}^{t_2} |x(t)|^2 dt
\]同樣，我們亦可擴展訊號的 total energy 到時間為無窮大 $t \in (-\infty, \infty)$，此定義為上式的極限，記做 $E_{\infty}$：
\[
E_{\infty} := \lim_{T \rightarrow \infty}\int_{-T}^{T} |x(t)|^2dt = \int_{-\infty}^{\infty}|x(t)|^2dt
\]若上式積分存在。

[動態規劃] 淺談 離散時間動態規劃 (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