[測度論] 期望值下確界與函數值下確界之恆等式

 Claim: 令 $(X, \mathcal{F})$ 為可測空間。令 $g: X \to \mathbb{R}$ 為可測函數，則 $$\inf_{\mathbb{P} \in \mathcal{P}(X)} \int_X g(x) d\mathbb{P}(x) = \inf_{x \in X} g(x)$$ 其中 $\mathcal{P}(X)$ 為 $(X, \mathcal{F})$ 上所有機率測度所成之集合。

Proof: 先證明 $\geq$: 對任意機率測度 $\mathbb{P}$，我們有 $$ g(x) \geq \inf_{x \in X}g(x) $$ 故取期望值不等式仍成立，亦即 $$ \mathbb{E}^\mathbb{P}[g(X)] = \int_X g(x) d\mathbb{P}(x) \geq \inf_{x \in X} g(x) $$  

以下接著證明 $\leq$: 固定 $\varepsilon > 0$，則由 infimum 定義，存在 $x_\varepsilon \in X$ 滿足 $$ g(x_\varepsilon) \leq \inf_x g(x) + \varepsilon \qquad (*) $$ 令 $\mathbb{P}:=\delta_{x_\varepsilon}$ (Dirac at $x_\varepsilon$ 滿足 $\delta_x(A):=1_{x \in A}$, $A \in \mathcal{F}$ ) 則 $$ \mathbb{E}^\mathbb{P}[g(X)] = \int_X g(x) d\delta_{x_\varepsilon} = g(x_\varepsilon) $$ 由$(*)$我們進一步得到 $$ \int_X g(x) d\delta_{x_\varepsilon} = g(x_\varepsilon) \leq \inf_x g(x) + \varepsilon $$ 對兩邊同取 $\inf_\mathbb{P}$ 可得 $$ \inf_\mathbb{P} \int_X g(x) d\mathbb{P}(x) \leq \inf_x g(x) + \varepsilon $$ 令 $\varepsilon \downarrow 0$ 得到 $\inf_{\mathbb{P} \in \mathcal{P}(X)} \int_X g(x) d\mathbb{P}(x) \leq \inf_{x \in X} g(x)$

Remark: (Dirac 測度不需單點可測）：在任意可測空間 $(X,\mathcal F) $上，對每個 $x\in X$ 定義 $\delta_x(A)=\mathbf 1_{\{x\in A\}}$ 其中 $A \in \mathcal F$，則 $\delta_x$ 是機率測度，且對一切 $\mathcal F$-可測 $g$ 有 $\int g\,d\delta_x=g(x)$。因此上述證明中以 $\delta_{x_\varepsilon}$ 作為選擇的測度不需要額外假設 $\{x\}\in\mathcal F$。

[最佳化] C^2 函數一階逼近的餘項積分表示

令 $f: \mathbb{R}^m \to \mathbb{R}$ 為 $C^2$-函數。對 $f$ 在 $y$ 附近使用一階泰勒展開：

\[ T_y(x) := f(y) + \nabla f(y)^\top (x - y) \]

則其餘項 $R(x,y)$ 訂為 $$R(x,y ):= T_y(x) - f(x)$$

現在定義單變數輔助函數 $g: [0,1] \to \mathbb{R}$ 滿足 $$g(t) : = f(y + t(x - y))$$現在觀察 $g(0) = f(y)$ 且 $g(1) = f(x)$。我們可以計算 $g(t)$ 導數透過多變數鏈鎖律：

$$g'(t) = \nabla f(y + t(x - y))^\top (x-y)$$且

$$g''(t) = (x-y)^\top \nabla^2 f(y + t(x-y)) (x-y)$$ 其中 $\nabla^2 f$ 為 $f$ 的 Hessian matrix。那麼由 Lemma 1可知單變數Taylor Theorem 對 $g(t)$ 在 $t=0$處展開有

$$g(1) = g(0) + g'(0)(1-0) + \int_0^1 g''(t) (1-t) dt \qquad (*)$$

現在代入 $g(1) = f(x), g(0)=f(y)$ 與 $g'(0) = \nabla f(y)^\top (x-y)$，上述 式$(*)$ 可改寫為

$$f(x) = \underbrace{ f(y) + \nabla f(y)^\top (x-y) }_{T_y(x)}+ \int_0^1  (1-t) (x-y)^\top \nabla^2 f(y + t(x-y)) (x-y) dt $$

因此，我們得到

$$f(x) - T_y(x) =  \int_0^1  (1-t) (x-y)^\top \nabla^2 f(y + t(x-y)) (x-y) dt $$

回憶餘項定法為 $R(x,y ):= T_y(x) - f(x)$，故我們有 

$$R(x,y) = -\int_0^1  (1-t) (x-y)^\top \nabla^2 f(y + t(x-y)) (x-y) dt$$

Lemma 1:

令 $g \in C^2([0,1])$，則單變數Taylor Theorem 對 $g(t)$ 在 $t=0$處展開有

$$g(1) = g(0) + g'(0) + \int_0^1 g''(t) (1-t) dt$$

Proof:

給定$g \in C^2$，回憶微積分基本定理對 $g$ 函數而言，

$$g(1) - g(0) = \int_0^1 g'(t) dt \qquad (**)$$

對於 $g'(t)$ 而言，我們亦可在次使用微積分基本定理：

$$g'(1) - g'(0) = \int_0^1 g''(s) ds$$

故對任意 $t \in [0,1]$ 我們有

$$g'(t) = g'(0) + \int_0^t g''(s) ds \qquad (@)$$

將 $(@)$ 代入 $(**)$ 得到

\begin{align*}g(1) - g(0) &= \int_0^1 [g'(0) + \int_0^t g''(s) ds] dt \\ &= \int_0^1 g'(0) dt + \int_0^1 \left( \int_0^t g''(s) \right) ds dt \qquad (@@)\end{align*}

其中

$$ \int_0^1 g'(0) dt = g'(0) \cdot t|_0^1 = g'(0)$$

且

\begin{align*} \int_0^1 \left( \int_0^t g''(s) \right) ds dt &= \int_0^1 \left( \int_0^t g''(s) \right) ds dt \\ &= \int_0^1 \int_0^1 1_{\{s < t\}} g''(s) ds dt \\ &= \int_0^1 \int_0^1 1_{\{s < t\}} g''(s) dt ds \\ &= \int_0^1 g''(s) \int_s^1 dt ds \\ &= \int_0^1 g''(s) (1-s) ds \\ \end{align*}

也就是說，$(@@)$ 可改寫為

$$g(1) - g(0)  =  g'(0) + \int_0^1 g''(s) (1-s) ds $$至此得證

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

連續函數有諸多用途，一般在參數最佳化領域中常見的情況是考慮所謂的上包絡函數(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 有興趣的讀者可以自行查閱。

[凸分析] 凸優化最佳解所成之集合為凸集

在凸優化問題中，僅管凸性保證了任意局部最優解 (local minimizer) 就是 全局最優解 (global minimizer)，但凸性並沒有保證所考慮的凸優化問題 一定 存在 最優極小解。下面的結果刻劃了凸優化最佳解的性質，常被用來檢驗最佳解的存在性，是個十分有用的結果。

===========
Theorem: (凸優化最佳解的集合為凸集)
令 $S \subseteq \mathbb{R}^n$ 為 凸集合 且 $f: S\to \mathbb{R}$ 為 凸函數。令 $S^*$ 為所有極小點所成之集合亦即
\[
S^* := \{x\in S: f(x) \leq f(y), \forall \; y \in S\; \}
\] 則 $S^*$ 為凸集。
===========

Proof: 若 $S^* = \emptyset$ 則上述定理陳述自動成立。若 $S^* \neq \emptyset$，則存在 $x_0 \in S^*$。考慮 level set
\[
S_{\leq f(x_0)} := \{x\in S: f(x) \leq f(x_0)\}
\] 則不難驗證 $S_{\leq f(x_0)} = S^*$。接著由下述 Lemma 可知 $S^*$ 為 convex。至此證明完畢。$\square$


===========
Lemma: 令 $S \subseteq \mathbb{R}^n$為凸集，且 $f:S \to \mathbb{R}$  為凸函數。則對任意 $\alpha \in \mathbb{R}$， (lower) level set
$$
S_{\leq \alpha}:= \{x \in \mathbb{R}^n: f(x) \leq \alpha\}
$$ 為 凸集。
===========

Proof: 若 level set $S_{\leq \alpha}$ 為空集合或者單點集，則陳述自動成立。若不然，取 $x_1,x_2 \in S_{\leq \alpha}$ ，則 $f(x_1) \leq \alpha$ 且 $f(x_2) \leq \alpha$。我們要證明 convex combination of $x_1$ 與 $x_2$ 仍落在 $S_{\leq \alpha}$ 之中，亦即我們要證明 $\lambda x_1 + (1-\lambda)x_2 \in S_{\leq \alpha}$。為此，我們利用 $f$ 的凸性，對任意 $\lambda \in (0,1)$，
\[
f(\lambda x_1 + (1-\lambda) x_2) \leq \lambda f(x_1) + (1-\lambda) f(x_2) \leq \alpha
\]故 $\lambda x_1 + (1-\lambda)x_2 \in S_{\leq \alpha}$，至此得證。$\square$


Remark:
對於 concave 函數我們仍有類似的結果記錄如下：

令 $S \subseteq \mathbb{R}^n$ 為凸集，且 $f: S \to \mathbb{R}$ 為 concave 函數，則所有極大點所成的集合 $S^*$ 為 convex set。

[最佳化] 對原最佳化問題的解是否能"回收"使用到新最佳化問題

令 $J: \mathbb{R}^n \to \mathbb{R}$ ，考慮以下最佳化問題
\[
\min_{x_1,x_2,...x_n} J(x_1,x_2,...x_n)  := J(x_1^*,x_2^*,...,x_n^*)
\]上述 $x_i^*$ 表示最佳解。現在考慮新的目標函數 $G: \mathbb{R}^n \to \mathbb{R}$ 為 上述的 $J:\mathbb{R}^n \to \mathbb{R}$ 額外加上新的函數 $F: \mathbb{R}^n \to \mathbb{R}$，亦即
\[
 G(x_1,x_2,...,x_n) :=J(x_1,x_2,...,x_n) + F(x_1,x_2,...,x_n)
\]我們想問前述獲得的最佳解 $x_1^*,x_2^*,.., x_n^*$ 是否仍然對新的目標函數成立？換句話說，是否能夠 "回收" 之前已經算好的最佳解  $ x_i^*$ 用在新的目標函數 $G$ 上呢。答案是否定的。考慮以下一個簡單的反例

Example:
對 $i=1,2,$，令 $x_i \in [-1,1]$並且將所有符合此條件的 $x_i$ 所成之集合記作 $\mathcal{X}$。現在考慮目標函數 $J(x_1,x_2) := x_1^2+x_2^2$ 並且 我們要求
$$
\min_{x_1,x_2 \in \mathcal{X}} J(x_1,x_2) =  \min_{x_1,x_2 \in \mathcal{X}}x_1^2+x_2^2
$$則最佳解不難發現為 $x_1^*=x_2^*=0$。現在我們考慮新的目標函數，將其記作
$$
G(x_1,x_2) :=J(x_1,x_2) + x_2
$$亦即 $G$ 為舊的目標函數 $J$ 額外加上 線性函數 $x_2 $。我們要求
$$
\min_{x_1,x_2 \in \mathcal{X}} G(x_1,x_2) = \min_{x_1,x_2 \in \mathcal{X}} x_1^2+x_2^2 + x_2
$$其最佳解變成 $x_1^* = 0$ 但 $x_2^* = -1/2 \;\;\; ( \neq 0)$。亦即舊的最佳解不能被"回收"使用。

Comments:
1. 上述謬誤偶爾能在文獻中發現。讀者應小心並盡量避免犯此錯誤。
2. 上述例子中若要使原最佳解可以被回收使用到新最佳解有很多方法，比如限制 可行集 $\mathcal{X}$ 將其改為 $0 \leq x_i \leq 1$ 便是一種。但是否符合需求又是另外一層考量。
3. 上述例子中若把 $\mathcal{X} := \mathbb{R}^2$，則有拘束最佳化問題變成無拘束最佳化問題，但反例仍然成立。

[凸分析] 定義在凸集上的凸函數其相對極小即為全域極小

這次要介紹凸分析 或者 凸優化 中可以說是最重要的結果：

=====================
Theorem:
令 $f$ 為 convex on convex set $\Omega$。則
1. $f$ 的任意 相對極小點 $x^*$ (local minimum of $f$ on $\Omega$) 必為 全域極小點 (global minimum of $f$ on $\Omega$)
2. 上述 凸函數的極小點所成之集合為凸集，亦即
\[
S:= \{x \in \Omega: f(x) = f(x^*)\}
\]為凸集。
=====================

Proof (1):
利用反證法：令  $y \in  \Omega$ 為 $f$在 $\Omega$ 上的相對極小點，亦即存在 $\delta>0$ 使得鄰域 $$N_\delta(y):=\{z: ||z-y|| < \delta\} \subset \Omega$$ 且
\[
f(y) \leq f(x), \;\;\; \forall x \in N_\delta(y)
\] 但 $y$ 不為 global minimum ：亦即 存在 $x^* \in \Omega$ 使得 \[
f(x^*) < f(y)
\]我們要證明此假設矛盾。

首先注意到 $x^* \notin N_\delta(y)$ 因為若不然，則 $f(y) \leq f(x^*)$ 此與假設不符。

現在，我們利用 $x^*$ 與 $y$ 來定義一個 新的點
\[
z(\alpha):= (1-\alpha)y +   \alpha x^*, \;\;\; \alpha := \frac{\delta}{2||y-x^*||} \in (0,1)
\]注意到 $z(\alpha) \in \Omega$ (因為 $\Omega$ 為凸集) 且我們觀察
\begin{align*}
  |z\left( \alpha  \right) - y|| &= \left\| {  (1 - \alpha ) y + \alpha x^* - y} \right\| \hfill \\
   &= \left\| {\alpha \left( {y - {x^*}} \right)} \right\| \hfill \\
   &\leqslant \frac{\delta }{{2\left\| {y - {x^*}} \right\|}}\left\| {y - {x^*}} \right\| < \delta  \hfill \\
\end{align*} 此結果表明
\[
z(\alpha) \in N_\delta(y)
\]
現在利用 $f$ 為凸函數性質，我們可知
\begin{align*}
  f(z(\alpha )) &= f\left( {\left( {1 - \alpha } \right) y + \alpha  x^*} \right) \hfill \\
   &\leqslant \alpha f\left( y \right) + \left( {1 - \alpha } \right)\underbrace {f\left( {{x^*}} \right)}_{f\left( y \right)} < f\left( y \right) \hfill \\
\end{align*} 上述不等式表明我們找到一個新的點 $z(\alpha) \neq y$  且 $z(\alpha) \in N_\delta(y)$ 使得
\[
f(z(\alpha)) < f(y)
\]此違反了 $y$ 在 $\Omega$ 上為相對極小的假設，得到矛盾。 故 $y$ 必為 global minimum。

Proof:(2)
接著我們證明上述 全域極小點 所成的集合為凸集，亦即
\[
S:= \{x \in \Omega: f(x) = f(x^*)\}
\]
首先觀察凸函數 $f$ 的 $c$-level set
\[
S_c:= \{x \in \Omega: f(x) \leq c\}
\]由下方的 FACT 可知 凸函數的 level set 為凸集。但因為 $x^*$ 為全域極小，故若取 $c:=f(x^*)$ 則
\[
S_c|_{c=f(x^*)}  =  \{x \in \Omega: f(x) \leq f(x^*)\} = \{x \in \Omega: f(x) = f(x^*)\} =S
\]由於為等式左方 $S_c$ 為凸集，故 $S$ 亦為凸集。 $\square$


=====================
FACT: 凸函數的 Sublevel Set 為凸集
令 $f: \Omega \subset \mathbb{R}^n \to \mathbb{R}$ 為凸函數，則對任意 $c \in \mathbb{R}$ ，其 $c$-level set
\[
S_c:=\{x: f(x) \leq c\}
\]為凸集
=====================

Proof: 利用凸函數定義，取 $x,y \in S_c$ 且 $\alpha \in [0,1]$ ，觀察
\[
f(\alpha x + (1-\alpha) y) \leq \alpha f(x) + (1-\alpha) f(y) \leq c 
\]故 $\alpha x + (1-\alpha) y\in S_c$，此表明 $S_c$ 為凸集。$\square$

[凸分析] 一階可導凸函數利用單點近似必定低估

Theorem: 
令 $f \in C^1$ 且 $f$ 為 convex on convex set $\Omega \subset \mathbb{R}^n$ 若且唯若 對任意 $x,y \in \Omega$ 而言，
\[
f(y) \geq f(x) + \nabla f(x) \cdot (y-x)
\]其中 $\nabla f(x) \cdot (y-x) := \nabla f(x)^T (y-x)$

給出證明之前我們先給一些直觀上的看法：

Comments:
1. 上述定理算是相當直覺，簡而言之就是說 affine (in $y$) function：
$ f(x) + \nabla f(x) (y-x)$ 可以做為 凸函數 $f$ 在 $x$ 點附近的 1 階 Taylor 近似，如下圖所示：

2. 注意到上述定理闡述的不等式對於所有 $x,y \in \Omega$ 都成立，也就是說透過 對$x$ 一階 Taylor 近似必定低估，一般 $f(x) + \nabla f(x) (y-x)$ 又稱作 global underestimaotr  of $f$。
3. 上述結果指出利用局部資訊 (一階導數) 可以得到 全域資訊 (global understametor )。
4. 若 $\nabla f(x) = 0$ 則對任意 $y \in \Omega$，我們有
\[
f(y) \geq f(x)
\]亦即 $x$ 為 全域及小點 (global minimizer) of $f$


以下我們給出證明

Proof: 先證明 $(\Rightarrow)$
令 $f \in C^1$ 且 $f$ 為 convex on convex set $\Omega \subset \mathbb{R}^n$，給定任意 $x,y \in \Omega$ ，我們要證
\[
f(y) \geq f(x) + \nabla f(x) (y-x)
\]
由於  $f$ 為 convex，令 $\alpha \in (0,1)$ 且定義
$$
z(\alpha) := \alpha x + (1-\alpha) y
$$則 $z(\alpha) \in \Omega$ 且由 $f$的凸性，我們有
\begin{align*}
  f\left( {z(\alpha )} \right) &= f\left( {\alpha x + \left( {1 - \alpha } \right)y} \right) \hfill \\
   &\leqslant \alpha f\left( x \right) + \left( {1 - \alpha } \right)f\left( y \right) \hfill \\
\end{align*} 由於 $\alpha \neq 0$ 我們可整理上式得到
\[\frac{{f\left( {\alpha x + \left( {1 - \alpha } \right)y} \right) - f\left( y \right)}}{\alpha } \leqslant f\left( x \right) - f\left( y \right)\]或者
\[\frac{{f\left( {y - \alpha \left( {y - x} \right)} \right) - f\left( y \right)}}{\alpha } \leqslant f\left( x \right) - f\left( y \right)\]取 $\alpha \to 0$，由於 $f\in C^1$ 我們不難看出上述不等式左方 為沿著 $y-x$ 的方向導數，故我們有
\[
\nabla f(y) \cdot (y-x) \leq f(x) -f(y)
\]或者
\[
f(x) \geq f(y) + \nabla f(y) \cdot (y-x)
\]上述結果對 任意 $x,y \in \Omega$ 成立，故我們將 $x,y$ 角色對換即得到定理要求的陳述。

接著我們證明$(\Leftarrow)$：
假設  對任意 $x,y \in \Omega$ 而言，
\[
f(y) \geq f(x) + \nabla f(x) (y-x) \;\;\;\;\; (**)
\]我們要證明 $f$ 為 convex。故令 $x_1, x_2 \in \Omega$ 與 $\alpha \in [0,1]$ ，並且我們額外定義
\[
\bar{x} := \alpha x_1 + (1- \alpha) x_2
'\]
則由假設可知 $x_1, x_2, \bar{x}$ 必定滿足 $(**)$，我們可寫下
\[\begin{gathered}
  f({x_1}) \geqslant f(\bar x) + \nabla f(\bar x)({x_1} - \bar x) \hfill \\
  f({x_2}) \geqslant f(\bar x) + \nabla f(\bar x)({x_2} - \bar x) \hfill \\
\end{gathered} \]現在對上述 第一條不等式 兩邊同乘上 $\alpha$ ，對 第二條不等式 兩邊乘上 $1- \alpha$ ，亦即
\begin{align*}
 & \alpha f({x_1}) \geqslant \alpha f(\bar x) + \alpha \nabla f(\bar x)({x_1} - \bar x) \hfill \\
  &\left( {1 - \alpha } \right)f({x_2}) \geqslant \left( {1 - \alpha } \right)f(\bar x) + \left( {1 - \alpha } \right)\nabla f(\bar x)({x_2} - \bar x) \hfill \\
\end{align*}
現在觀察
\begin{align*}
  \alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) &\geqslant \alpha f(\bar x) + \alpha \nabla f(\bar x)({x_1} - \bar x) \hfill \\
   & \hspace{10mm}+ \left( {1 - \alpha } \right)f(\bar x) + \left( {1 - \alpha } \right)\nabla f(\bar x)({x_2} - \bar x)
\end{align*}
將上式稍微做一下整理可得
\begin{align*}
  &\alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) \geqslant f(\bar x) + \nabla f(\bar x)\left( {\alpha ({x_1} - \bar x) + \left( {1 - \alpha } \right)({x_2} - \bar x)} \right) \hfill \\
   &\Rightarrow \alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) \geqslant f(\bar x) + \nabla f(\bar x)\underbrace {\left( {\alpha {x_1} + \left( {1 - \alpha } \right){x_2} - \bar x} \right)}_{ = 0} \hfill \\
   &\Rightarrow \alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) \geqslant f(\bar x) \hfill \\
\end{align*} 上述不等式表明 $f$ 為凸函數。$\square$


Comments:
1. 若 $f$ 為 $C^1$ strict convex 函數 on $\Omega$，則對任意 $x,y \in \Omega$ 而言，
\[
f(y) >f(x) + \nabla f(x) (y-x)
\]
2. 若 $f$ 為 concave 則利用 $-f$ 為 convex 特性可知 對於 concave 函數而言，定理的不等式變成: 對任意 $x,y \in \Omega$ 而言，
\[
f(y) \leq f(x) + \nabla f(x) (y-x)
\]

[最佳化理論] 有限維空間 無拘束最佳化問題 的二階充分必要條件

回顧先前我們討論過的一階必要條件，以下我們介紹有限維空間 無拘束最佳化問題的 二階必要與充分條件，首先是二階必要條件：

====================
Theorem: Second-Order Necessary Condition, SONC:
令 $S \subset \mathbb{R}^n$ 且令 $f \in C^2(S)$。若 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 則 對 在點 ${\bf x}^*$的任意可行方向 ${\bf d} \in \mathbb{R}^n$，我們有
1. $\nabla f({\bf x}^*) \cdot {\bf d} \geq 0$
2. 若 $\nabla f({\bf x}^*) \cdot {\bf d} =0$ 則 ${\bf d}^T \nabla^2 f({\bf x}^*) {\bf d} \geq 0$
====================

Proof: 由於 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 對條件 1 自動成立。我們僅需證明條件2。令  ${\bf d} \in \mathbb{R}^n$ 為 在點 ${\bf x}^*$的任意可行方向，故存在 $\bar \alpha >0$ 使得
\[
{\bf x}(\alpha) := {\bf x}^* + \alpha {\bf d} \in S, \;\;\; \alpha \in [0, \bar \alpha]
\] 利用 Taylor Theorem 對 ${\bf x}^*$ 展開到二階項，我們可得
\begin{align*}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right)
&= f\left( {{{\mathbf{x}}^*}} \right) + \nabla f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) \\
& \hspace{10mm}+ \frac{1}{2}{\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right)^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
 &  = f\left( {{{\mathbf{x}}^*}} \right) + \alpha \nabla f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + \frac{1}{2}{\alpha ^2}{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
\end{align*} 若 $\nabla f({\bf x}^*) \cdot {\bf d} =0$，則上式變成
\[\begin{gathered}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right) = f\left( {{{\mathbf{x}}^*}} \right) + \alpha \underbrace {\nabla f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}}_{ = 0} + \frac{1}{2}{\alpha ^2}{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
   = f\left( {{{\mathbf{x}}^*}} \right) + \frac{1}{2}{\alpha ^2}{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
\end{gathered} \]注意到若 ${{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}<0$ 對在 ${\bf x}^*$任意足夠小的鄰域成立，則我們得到
\[f\left( {{\mathbf{x}}\left( \alpha  \right)} \right) \leqslant f\left( {{{\mathbf{x}}^*}} \right) \]此與 ${\bf x}^*$是 local minimum 矛盾，故
\[{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}\geq 0\;\;\;\;\; \square\]

Comments:
讀者大概不難發現不論是一階 或者 二階 條件都是依賴 Taylor 定理展開式，依此來進行估計。


====================
Theorem: Second-Order Sufficient Condition, SOSC: 
令 $S \subset \mathbb{R}^n$ 且令 $f \in C^2(S)$。若 ${\bf x}^* \in S$ 為interior point 且$f$ 在該點有定義。若
1. $\nabla f({\bf x}^*) = {\bf 0}$
2. 對任意 ${\bf d} \in \mathbb{R}^n$，${\bf d}^T \nabla^2 f({\bf x}^*) {\bf d} >0$ (亦即 $\nabla^2 f({\bf x}^*)$ 為 positive definite)，則 ${\bf x}^*$ 為 strict local minimum of $f$。
====================

Proof: 要證明 ${\bf x}^*$ 為 strict local minimum of $f$，我們僅需證明
\[
f({\bf x}) > f({\bf x^*}),\;\;\; {\bf x} \in \mathcal{N}({\bf x}^*)
\]其中 $\mathcal{N}({\bf x}^*)$為以 ${\bf x}^*$為中心所成之鄰域。故令  ${\bf d} \in \mathbb{R}^n$ 為 在點 ${\bf x}^*$的任意可行方向，故存在 $\bar \alpha >0$ 使得
\[
{\bf x}(\alpha) := {\bf x}^* + \alpha {\bf d} \in \mathcal{N}({\bf x}^*), \;\;\; \alpha \in [0, \bar \alpha]
\] 利用 $\nabla f({\bf x}^*) = {\bf 0}$ 與 Taylor Theorem 對 ${\bf x}^*$ 展開到二階項，我們可得
\begin{align*}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right)
&= f\left( {{{\mathbf{x}}^*}} \right) + \nabla f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) \\
& \hspace{10mm}+ \frac{1}{2}{\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right)^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
 &  =f\left( {{{\mathbf{x}}^*}} \right) + \frac{1}{2}{\alpha ^2}\underbrace {{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}}_{ > 0} + o\left( {{{\left\| \alpha  \right\|}^2}} \right)\\
\end{align*}
在 ${\bf x}^*$ 附近鄰域而言，我們可推知
\[\begin{gathered}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right) - f\left( {{{\mathbf{x}}^*}} \right) = \frac{1}{2}{\alpha ^2}\underbrace {{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}}_{ > 0} > 0 \hfill \\
\end{gathered} \]則 ${\bf x}^*$ 為 strict local minimum of $f$ 。$\square$

[最佳化理論] 有限維空間 無拘束最佳化問題 的一階必要條件

令 $f: \mathbb{R}^n \to \mathbb{R}$ 且 $S \subset \mathbb{R}^n$ 為 feasible set，現在我們考慮以下的 有限維度 (無拘束)最佳化問題
\[\begin{gathered}
  \min f\left( {\mathbf{x}} \right) \hfill \\
  s.t.{\mathbf{x}} \in S \subseteq {R^n} \hfill \\
\end{gathered} \]
我們想知道上述最佳化問題是否有解？ 若有解則是哪一種解 e.g., 局部最佳解(local optimum)或者 全域最佳解(global optimum)？)，以及上述的解是否能夠透過某種方法來將其描述。

對於上述有限維度最佳解的存在性問題一般可由 Weierstrass Extremum Theorem處理，在此不做贅述。以下我們討論 最佳解 存在的必要條件：更近一步地說是 最佳化理論中的求取 "局部最佳解" 的一階必要條件，在給出結果之前我們首先定義 可行方向 (Feasible Direction)如下：

======================
Definition: Feasible Direction
給定 ${\bf x} \in S \subset \mathbb{R}^n$，我們說向量 ${\bf d} \in \mathbb{R}^n$ 為 在 ${\bf x}$ 處的可行方向 (feasible direction) 若下列條件成立：存在常數 $ \bar \theta >0$ 使得對任意 $\theta \in [0, \bar\theta]$而言，
\[
{\bf x} + \theta {\bf d} \in S
\]======================

有了以上的可行方向的想法，其局部最佳解的一階必要條件有如下陳述：

======================
Theorem: (First Order Necessary Condition, FONC): 令 $S \subset \mathbb{R}^n$ 且 $f \in C^1$ on $S$ ( $f$為一階可導且導數連續)。若 ${\bf x}^*$ 為 局部極小解 (local minimum point) of $f$ over $S$ 則 對任意在 ${\bf x}^*$點上的可行方向 ${\bf d} \in \mathbb{R}^n$，我們有
\[
\nabla f({\bf x}^*) \cdot {\bf d} \geq 0
\]======================

Proof: 用反證法：假設存在 對 ${\bf x}^*$點上的可行方向 ${\bf d} \in \mathbb{R}^n$ 使得
\[
\nabla f({\bf x}^*) \cdot {\bf d} <0
\]我們要證明矛盾。由 ${\bf d}$為可行方向之定義可知： 存在 $\bar \theta>0$ 使得 對任意 $\theta \in [0, \bar\theta]$而言，
\[
{\bf x}(\theta) = {\bf x}^* + \theta {\bf d} \in S
\]現在利用 Taylor 定理對 ${\bf x}^*$ 展開 且利用已知假設 $\nabla f({\bf x}^*) \cdot {\bf d} <0$ 可得
\begin{align*}
  f({\bf x}(\theta )) &= f({{\bf x}^*}) + \nabla f({{\bf x}^*})({\bf x}(\theta ) - {{\bf x}^*}) + o(||{\theta}||) \hfill \\
   &= f({{\bf x}^*}) + \nabla f({ {\bf x}^*})\theta {\bf d} + o(||{\theta}||) \hfill \\
  &= f({{\mathbf{x}}^*}) + \theta \underbrace {\nabla f({{\mathbf{x}}^*}) \cdot {\mathbf{d}}}_{ < 0} + o(||\theta ||)\\
   &< f({{\bf x}^*}) \hfill \\
\end{align*} 上述最後一條不等式當 $\theta$ 足夠小的時候成立。此與我們假設 $ f({ {\bf x}^*})$最小 矛盾。$\square$


Comments:
1. 上述一階必要條件說明了若我們已經處在局部最佳解的位置 ${\bf x}^*$ 則 沿著任何其他可行方向移動都會增加 目標函數值，亦即 \[
\nabla f({\bf x}^*) \cdot {\bf d} \geq  0
\]
2. 關於上述使用的 Taylor Theorem 與 little-oh 符號可以參考下方補充說明或者BLOG其餘相關文章。
3. 上述 一階必要條件 使用到了 Taylor Theorem 的一階項，一般而言可視為對 ${\bf x}^*$ 的 一階近似估計。


若為無拘束情況，則我們有以下的衍生結果：

=======================
Corollary: 令 $S \subset \mathbb{R}^n$ 且 $f \in C^1$ on $S$。若 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 且若 ${\bf x}^*$ 為 $S$ 的 interior point 則
\[
\nabla f({\bf x}^*) = {\bf 0}
\]=======================

Proof Sketch: 設 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 則由 FONC可知 對任意 ${\bf d} \in \mathbb{R}^n$ 為在 ${\bf x}^*$點上的可行方向，我們有
\[
\nabla f({\bf x}^*) \cdot {\bf d} \geq 0
\]但因為 若 ${\bf x}^*$ 為 $S$ 的 interior point 故在該點${\bf x}^*$ 之可行方向 ${\bf d}$ 可在足夠小的區域選為任意方向來移動， 若 ${\bf d} \neq {\bf 0}$ 則 為了使 $\nabla f({\bf x}^*) \cdot {\bf d} = 0$ 成立，我們必定要求
\[
\nabla f({\bf x}^*)  = {\bf 0} \;\;\;\; \square
\]

Comments:
1. 上述 FONC 無拘束情況的 FONC 將原本最佳問題轉成求解有 $n$ 未知數的 $n$ 個系統方程問題。

2. 上述結果只保證 局部最佳解 的必要條件，如果想要得到充分條件，通常需要 $f$ 的二階導數的資訊也就是 需要 Hessian matrix，一般稱作 Second-Order Sufficient Condition我們在以後的文章會在提及。

3. 局部最佳解落在集合 $S$ 上，若此集合為 convex 集合且 $f$ 為凸函數，則局部最佳解 為 全域最佳解，相關結果可翻閱本 BLOG關於最佳化與凸分析的文章。

4. 上述最佳化問題可以與 變分不等式問題 (Variational Inequality Problem)等價，以下給出此結果：

=======================
Corollary 2: Relationship Between Optimization and Variational Inequality
令 $S \subset \mathbb{R}^n$ 為 convex 且 $f : \mathbb{R}^n \to \mathbb{R}$ 且 $f \in C^1(S)$。若 ${\bf x}^*$ 為 local minimum point of $f$ over $S$，則 ${\bf x}^*$ 為下列變分不等式問題的解：求 ${\bf x} \in S$ 使得 對任意 ${\bf x}' \in S$ 而言，
\[
\langle {\bf x}' - {\bf x}, \nabla f({\bf x}) \rangle \geq 0
\]=======================
Proof:
取 ${\bf d} := {\bf x}' - {\bf x}$ 即可。



補充定理
==================
Taylor Theorem with Little-oh Notation:
假設 $X \subset \mathbb{R}^n$ 為 open，且 ${\bf x} \in X$ 與 $f: X \to \mathbb{R}$ 為 $C^1$。則
\[
f({\bf x}+{\bf h}) = f({\bf x}) + Df({\bf x}) {\bf h} + o(||{\bf h}||),\;\;\; ||h|| \to 0
\]==================

[最佳化] 無窮維 Weierstrass 極值定理

===================
Infinite Dimensional Weierstrass Extreme Theorem: 令 $X$ 為 normed vector space 且 $K \subset X$ 為 compact set 。若 $f$ 為 upper semicontinuous on $K$ 則 $f$ 在 $K$上可達到最大值，亦即 存在 $x \in K$ 使得
\[
f(x) = \sup_{z \in K} f(z)
\]=================

Proof: 要證明 $f$ 在 $K$上可達到最大值，亦即 存在 $x \in K$ 使得
\[
f(x) = \sup_{z \in K} f(z)
\]上式等價為
\[
f(x) \geq \sup_{z \in K} f(z) \;\;\; \text{and } \; f(x) \leq \sup_{z \in K} f(z)
\] 注意到 $f(x) \leq \sup_{z\in K} f(z)$ 為顯然，故以下僅需證明\[
f(x) \geq \sup_{z \in K} f(z)
\]

令 $M:= \sup_{z \in K} f(z)$ 則由 supremum 性質可知存在數列 $\{x_n\} \subset K$ 使得
\[
f(x_n) \to M \;\;\;\; (*)
\] 由於 $K$ 為 compact 故必定存在 $\{x_n\} $ 的子數列 記作 $\{x_{n_k}\} \subset K$  使得其收斂在  $x \in K$ ，另外由於 式子 $(*)$ ，我們可推知子數列亦滿足
\[
f( x_{n_k} ) \to M \;\;\;\; \text{ as $k \to \infty$}
 \]由於 $f$ 為 upper semicontinuous on $K$，可知
\[
\limsup_{z \to x} f(z) \leq f(x)
\]由 $\limsup$ 的 數列與函數數列性質 ，我們有
\[
\limsup_{k \to \infty} f(x_{n_k}) \leq \limsup_{z \to x} f(z)  \;\;\;\; (\star)
\]又因為 $f(x_{n_k}) \to M$ (as $k \to \infty$) 故
\[
\limsup_{k \to \infty} f(x_{n_k}) = \lim_{k \to \infty} f(x_{n_k}) = M \;\;\;\; (**)
\]由 $(\star)$ 與 $(**)$，我們有
\[
M=\lim_{k \to \infty} f(x_{n_k}) \leq f(x)
\]即為所求 $\square$

[投資組合理論] Markowitz 最小變異 投資組合 之解

此文將討論 Markowitz 在1952年針對單期 報酬-風險 投資策略所建構的最小變異投資組合理論，簡而言之就是以期望報酬為報酬，風險變異（或者標準差）為風險，試問如何建構一組投資組合並對各資產給定適當的權重使得風險變異被最小化。推導過程會用到一些必要的最佳化與線性代數的知識。


Markowitz 的單期投資組合的描述：
假設在期初手上有 $V(0)>0$ 資產，現在我們打算在期初時購入 $n$ 種 互為相關 的風險資產 (correlated risky asset) 用以建構投資組合，其個別資產之隨機報酬 表為 $r_1,r_2,....,r_n$ 且對應的 期望報酬 為 $E[r_1], E[r_2],...,E[r_n]$ 與 風險變異 $Var(r_1), Var(r_2),...,Var(r_n)$ 且 資產之間的共變異 為 $ Cov(r_i,r_j),\;\; \forall \;  i,j=1,2,...,n$。

Comments:
為求分析簡便在以下分析中，我們建構的投資組合不考慮無風險資產 (risk-free asset)，亦即 $Var(r) = 0$ 的資產我們不考慮。


投資策略： 
對於第 $i$ 資產之投資策略為對 $i=1,2,...,n$，令在期初之投資策略為 $I(0)$ 滿足
$$I_i(0) = K_i V(0)$$ 其中我們要求 $\sum_{i=1}^n K_i = 1$ （但允許 $K_i$ 為負值，亦即我們允許賣空）。

Comments:
一般投資書籍在討論上述投資策略或者廣義的資產配置問題時，多半僅稱呼 $K_i$ 為權重，且對於整體投資策略不多著墨。不過事實上，此類問題可以透過引入 控制理論 觀點，將投資策略視為標準回授控制 $I=KV$。


單期資產動態模型：
則我們打算持有單期 (比如說 一年) 則期末資產為
\[
V(1) = V(0) + \sum_{i=1}^n K_i r_i V(0)
\]則不難得知我們投資組合的期末報酬，記作 $r_p$ 可由上式推得為
\[{r_p}: = \frac{{V(1) - V(0)}}{{V(0)}} = \sum\limits_{i = 1}^n {{K_i}} {r_i}\]
並且回憶投資組合的期望收益率為
\[
E[{r_p}]  =  \sum\limits_{i = 1}^n {{K_i}} E[{r_i}]\]且對應的變異可表為
\begin{align*}
  Var[{r_p}] = \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{K_i}{K_j}Cov\left( {{r_i},{r_j}} \right)} }
\end{align*}


Markowitz 最小變異投資組合 的等價 最佳化問題：
為求取 Markowitz 最小變異投資組合，我們首先給定任意投資組合之期望報酬 $\widehat{r}$，並接著建構以下最佳化問題：
\begin{align*}
  &\min \frac{1}{2}Var\left( {{r_p}} \right) \hfill \\
  &s.t. \hfill \\
  &E\left[ {{r_p}} \right] = \widehat r\;\;; \hfill \\
  &\sum\limits_{i = 1}^n {{K_i}}  = 1 \hfill \\
\end{align*}
利用前面推得的 $E[r_p], Var(r_p)$ ，我們知道上述最佳化問題等價為
\begin{align*}
 & \min \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{K_i}{K_j}Cov\left( {{r_i},{r_j}} \right)} }  \hfill \\
  &s.t. \hfill \\
 & \sum\limits_{i = 1}^n {{K_i}} E[{r_i}] = \widehat{r}\;\;\;; \hfill \\
  &\sum\limits_{i = 1}^n {{K_i}}  = 1 \hfill \\
\end{align*}
注意到上述問題為具有等式拘束的最佳化問題。常用的工具為透過 Lagrange Multiplier 來求解必要條件。

Comments:
熟習線性代數與最佳化的讀者，應不難看出上述最佳化問題可被化約為二次規劃 (Quadratic Programming)問題，亦即目標函數為二次函數，且具有線性拘束的最佳化問題：
\[\begin{gathered}
  \min {K^T}CK \hfill \\
  s.t. \;\; AK = b \hfill \\
\end{gathered} \]其中 $K:= [K_1,K_2,...,K_n]$ 且 $C$ 為共變異矩陣其中第 $ij$ 個元素為 $Cov(r_i,r_j)$ 且拘束條件為
\[\underbrace {\left[ {\begin{array}{*{20}{c}}
  {E[{r_1}]}&{E[{r_2}]}& \cdots &{E[{r_2}]}&{E[{r_2}]} \\
  1&1& \cdots &1&1
\end{array}} \right]}_A\underbrace {\left[ \begin{gathered}
  {K_1} \hfill \\
  {K_2} \hfill \\
   \vdots  \hfill \\
  {K_n} \hfill \\
\end{gathered}  \right]}_K = \underbrace {\left[ \begin{gathered}
  {\hat r} \hfill \\
  1 \hfill \\
\end{gathered}  \right]}_b\]注意到上述討論中， $A$ 矩陣的變數比方程多，這一類特殊問題又稱 minimum norm problem，(亦即我們將 $min K^TCK := min ||K||_C^2 $) 。若 $rank(A) = 2$ 且共變異矩陣為正定矩陣 則有立刻的唯一解滿足拘束 $AK=b$ 且最小化 $K^TCK$ ，記作 $K^*$，如下
\[{K^*} = {C^{ - 1}}{A^T}{(A{C^{ - 1}}{A^T})^{ - 1}}b\]有興趣讀者可以自行驗證。以下我們將採用 Lagrange Muliplier 方式來求解此問題。


求解最小變異投資組合 (利用 Lagrange Multiplier)：
因為我們有兩條等式拘束，故可取 $\lambda, \mu$ 為 Lagrange Multiplier 並建構 Largangian 函數 $L$ 如下：
\[L: = \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{K_i}{K_j}Cov\left( {{r_i},{r_j}} \right)} }  - \lambda \left( {\sum\limits_{i = 1}^n {{K_i}} E[{r_i}] - \widehat r} \right) - \mu \left( {\sum\limits_{i = 1}^n {{K_i}}  - 1} \right)\]
欲求最佳解的必要條件，我們對求 $L$ 每一個 $K_k$ ($k=1,2,...,n$) 之偏導數並令其為零，注意到在此我們須對雙重加總求導，一般常用的結果為以下 FACT:

==========
FACT: 有限雙重加總的求導
$$
\frac{\partial}{\partial x_k} \sum_{i, j} a_{i j} x_i x_j
   = \sum_{i, j} a_{i j}
             \left( \frac{\partial x_i}{\partial x_k} x_j
                      + x_i \frac{\partial x_j}{\partial x_k} \right)
   = \sum_j a_{k j} x_j + \sum_i a_{i k} x_i
$$==========

故利用上述 FACT ，我們首先對 Lagragian $L$ 的第一項雙重加總取導可得
\begin{align*} \frac{\partial }{{\partial {K_k}}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{K_i}{K_j}Cov\left( {{r_i},{r_j}} \right)} } &= \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {Cov\left( {{r_i},{r_j}} \right)\left( {\frac{\partial }{{\partial {K_k}}}{K_i}{K_j}} \right)} } \hfill \\ &= \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {Cov\left( {{r_i},{r_j}} \right)\left( {\frac{{\partial {K_i}}}{{\partial {K_k}}}{K_j} + {K_i}\frac{{\partial {K_j}}}{{\partial {K_k}}}} \right)} } \hfill \\ &= \sum\limits_{j = 1}^n {Cov\left( {{r_k},{r_j}} \right){K_j}} + \sum\limits_{i = 1}^n {Cov\left( {{r_i},{r_k}} \right){K_i}} \hfill \\ \end{align*}
利用 $Cov(r_i, r_j) = Cov(r_j, r_i)$ 的對稱性，我們可以計算出對 Lagragian  $L$ 的求導：
\[\begin{gathered}
  \frac{\partial }{{\partial {K_k}}}L = \frac{1}{2}\left( {\sum\limits_{j = 1}^n {Cov\left( {{r_k},{r_j}} \right){K_j}}  + \sum\limits_{i = 1}^n {Cov\left( {{r_i},{r_k}} \right){K_i}} } \right) - \lambda E[{r_k}] - \mu  \hfill \\
   = \sum\limits_{i = 1}^n {Cov\left( {{r_k},{r_i}} \right){K_i}}  - \lambda E[{r_k}] - \mu  \hfill \\
\end{gathered} \]令此式為零可得 $n$ 條等式：
\[\frac{\partial }{{\partial {K_k}}}L = \sum\limits_{i = 1}^n {Cov\left( {{r_k},{r_i}} \right){K_i}}  - \lambda E[{r_k}] - \mu : = 0, k=1,2,...,n\]
現在總結以上討論，我們得到以下結果：


給定期望報酬 $\widehat{r}$ 對各項資產之權重 $K_i$ $(i=1,2,..,n)$ 組成之投資組合可達成最小化變異的必要條件為 $K_i, \lambda, \mu$ 同時滿足下列 $n+2$ 條等式：
\[\left\{ \begin{align*}
  & \sum\limits_{i = 1}^n {Cov\left( {{r_k},{r_i}} \right){K_i}}  - \lambda E[{r_k}] - \mu  = 0,\;\;\; k=1,2,...,n \hfill \\
 & \sum\limits_{i = 1}^n {{K_i}} E[{r_i}] = \hat r \hfill \\
 & \sum\limits_{i = 1}^n {{K_i}}  = 1 \hfill \\
\end{align*}  \right.\]
一般而言，我們有 $n+2$ 條等式與 $n+2$ 個變數，依照線性代數基本定理可知若這些式子滿足 full row & colum rank條件，則上述方程組有解 $\{(K_1,K_2,...,K_n), \lambda, \mu\}$ 且此組解為唯一解。


Example: 不相關資產的例子：
假設有三種互不相關的資產，且假設
$$E[r_1] = 1, E[r_2] = 2, E[r_3] = 3
$$與
$$Var(r_1) = Var(r_2) = Var(r_3) = 1$$
且因為互不相關，各資產之間共變異為 $Cov(r_i,r_j) =0, \; \forall i,j=1,2,3, i \neq j$。

(a) 試求投資組合的風險變異 $Var[r_p]$ 與期望報酬 $E[r_p]$
(b) 現在給定任意 $\widehat{r}$ 試求出 $K_1^*,K_2^*,K_3^*$ 使得我們達成最小變異投資組合且滿足 $\sum_{i=1}^3 K_i = 1$ 且 $E[r_p] = \widehat{r}$。
(c) 利用part(b) 之解，試問最小變異之值為何？

Solution (a)
計算投資組合的期望報酬與變異如下：\[\begin{array}{l}
E[{r_p}] = {K_1}E[{r_1}] + {K_1}E[{r_1}] + {K_1}E[{r_1}]\\
 = {K_1} + 2{K_1} + 3{K_1}
\end{array}\]
同理
\[\begin{array}{l}
Var[{r_p}] = \sum\limits_{i = 1}^3 {\sum\limits_{j = 1}^3 {{K_i}{K_j}Cov\left( {{r_i},{r_j}} \right)} } \\
 = {K_1}{K_1}Cov\left( {{r_1},{r_1}} \right) + {K_2}{K_2}Cov\left( {{r_2},{r_2}} \right) + {K_3}{K_3}Cov\left( {{r_3},{r_3}} \right)\\
 = {K_1}^2 + {K_2}^2 + {K_3}^2
\end{array}\]

Solution (b)
定義 Lagrange Multiplier $\mu,\lambda$ ，並使用前述討論的結果
\[\begin{array}{l}
\left\{ \begin{array}{l}
\sum\limits_{i = 1}^n {Cov\left( {{r_k},{r_i}} \right){K_i}}  - \lambda E[{r_k}] - \mu  = 0,\;\;\;k = 1,2,...,n\\
\sum\limits_{i = 1}^n {{K_i}} E[{r_i}] = \hat r\\
\sum\limits_{i = 1}^n {{K_i}}  = 1
\end{array} \right.\\
 \Rightarrow \left\{ \begin{array}{l}
Cov\left( {{r_1},{r_1}} \right){K_1} - \lambda E[{r_1}] - \mu  = 0\\
Cov\left( {{r_2},{r_2}} \right){K_2} - \lambda E[{r_2}] - \mu  = 0\\
Cov\left( {{r_3},{r_3}} \right){K_3} - \lambda E[{r_3}] - \mu  = 0\\
{K_1}E[{r_1}] + {K_2}E[{r_2}] + {K_3}E[{r_3}] = \hat r\\
{K_1} + {K_2} + {K_3} = 1
\end{array} \right.\\
 \Rightarrow \left\{ \begin{array}{l}
{K_1} - \lambda  - \mu  = 0\\
{K_2} - 2\lambda  - \mu  = 0\\
{K_3} - 3\lambda  - \mu  = 0\\
{K_1} + 2{K_2} + 3{K_3} = \hat r\\
{K_1} + {K_2} + {K_3} = 1
\end{array} \right.
\end{array}\]上式可改寫為矩陣形式如下
\[\left\{ \begin{array}{l}
{K_1} - \lambda  - \mu  = 0\\
{K_2} - 2\lambda  - \mu  = 0\\
{K_3} - 3\lambda  - \mu  = 0\\
{K_1} + 2{K_2} + 3{K_3} = \hat r\\
{K_1} + {K_2} + {K_3} = 1
\end{array} \right. \Rightarrow \left[ {\begin{array}{*{20}{c}}
1&0&0&{ - 1}&{ - 1}\\
0&1&0&{ - 2}&{ - 1}\\
0&0&1&{ - 3}&{ - 1}\\
1&2&3&0&0\\
1&1&1&0&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{K_1}}\\
{{K_2}}\\
{{K_3}}\\
\lambda \\
\mu
\end{array}} \right] = \left[ \begin{array}{l}
0\\
0\\
0\\
{\hat r}\\
1
\end{array} \right]\]上述等式左方之矩陣為 full rank 反矩陣存在，故
\[\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{{K_1}}\\
{{K_2}}\\
{{K_3}}\\
\lambda \\
\mu
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{1/6}&{ - 1/3}&{1/6}&{ - 1/2}&{4/3}\\
{ - 1/3}&{2/3}&{ - 1/3}&0&{1/3}\\
{1/6}&{ - 1/3}&{1/6}&{1/2}&{ - 2/3}\\
{1/2}&0&{ - 1/2}&{1/2}&{ - 1}\\
{ - 4/3}&{ - 1/3}&{2/3}&{ - 1}&{7/3}
\end{array}} \right]\left[ \begin{array}{l}
0\\
0\\
0\\
{\hat r}\\
1
\end{array} \right]\\
 \Rightarrow \left\{ \begin{array}{l}
{K_1}^* =  - \hat r/2 + 4/3\\
{K_2}^* = 1/3\\
{K_3}^* = \hat r/2 - 2/3
\end{array} \right.
\end{array}\]上述 $K_1,K_2,K_3$ 即為最佳解 (最小變異解)

Solution (c)
將 par(b) 的最佳解結果代回組合變異：
\[\begin{array}{l}
Var[{r_p}] = {K_1}^2 + {K_2}^2 + {K_3}^2\\
 \Rightarrow Var{[{r_p}]^*} = {\left( {\frac{{ - \hat r}}{2} + \frac{4}{3}} \right)^2} + {\left( {\frac{1}{3}} \right)^2} + {\left( {\frac{{\hat r}}{2} - \frac{2}{3}} \right)^2} = \frac{1}{2}(\hat r - 4)\hat r + \frac{7}{3}
\end{array}\]我們可進一步繪製 $\bar{r}$ vs $\sigma$ 如下圖所示 (這裡我們取標準差 $\sigma := \sqrt{Var(\cdot)}$ )

觀察上圖可以發現當 $\widehat{r} = 2$，我們可得到最小變異。

Comments:
一般金融文獻或者投資學文獻通常不繪製上圖，他們會將 x 軸定為風險 (亦即 $\sigma$) 並將 $y$ 軸定為 報酬。亦即一般繪製成類似 "子彈" 的形狀如下：

注意到此曲線的上半部稱為 efficient frontier，亦即給定任意報酬 $\widehat{r}$，上半部可得到較小風險。

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

此文主要討論離散泛函的極值與其必要條件，也就是所謂的離散版本的 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

[系統理論] 離散時間系統的穩定度理論 (1) - Lyapunov Stability Theory

延續前篇 [系統理論] 離散時間系統的穩定度理論 (0) - 先備概念，我們現在可以開始介紹 Lyapunov Stability Theory。

Definition: Lyapunov Function
一個函數 $V: \mathbb{R}^n \to \mathbb{R}_{\ge 0}$ 被稱作為 Lyapunov function for system $x^+ = f(x)$ 與 集合 $\mathcal{A}$ 若下列條件成立：
存在函數 $\alpha_1(\cdot), \alpha_2(\cdot),  \alpha_3(\cdot) \in \mathcal{K}_\infty$ 使得對任意 $x \in \mathbb{R}^n$，

1. $V(x) \ge \alpha_1(|x|_\mathcal{A})$
2. $V(x) \le \alpha_2(|x|_\mathcal{A})$
3. $V(f(x)) - V(x) \le -\alpha_3(|x|_\mathcal{A})$

給定 $\mathcal{A}$ 為 closed positive invariant for $x^+ = f(x)$ 且 $\mathcal{A} \subset X$ 我們說函數  $V(\cdot)$ 為 Lyapunov function in $X$ for system $x^+ = f(x)$ 與 集合 $\mathcal{A}$ 若下列條件成立：

對任意 $x \in X$，$V(\cdot)$ 滿足上述三條不等式。

上述 Lyapunov function 與 globally asymptotically stable 息息相關，事實上此Lyapunov function 的存在性為 globally asymptotically stable 的充分條件。我們將此記做下方結果

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

Theorem 1: Existence of Lyapunov Function Implies Globally Asymptotically Stability

假設 $V(\cdot)$ 為 Lyapuonv function for $x^+ = f(x)$ 與 $\mathcal{A}$，則 $\mathcal{A}$ 為 globally asymptotically stable。

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

Proof:

我們要證 $\mathcal{A}$ 為 globally asymptotically stable。故須證明

1. $\mathcal{A}$ 為 locally stable。

2. $\mathcal{A}$ 為 globally attractive。

先證 $\mathcal{A}$ 為 locally stable：給定 $\varepsilon >0$ 我們要找出 $\delta >0$ 使得對任意 $i \in \mathbb{Z}_{\ge 0}$ $|x|_\mathcal{A}< \delta \Rightarrow |\phi(i;x)|_\mathcal{A} <\varepsilon $

由於 $V$ 為 Lyapunov function 故由其定義可繪製下圖幫助我們選擇 $\delta$

令 $\delta : = \alpha _2^{ - 1}\left( {{\alpha _1}\left( \varepsilon  \right)} \right)$；現在給定 $i \in \mathbb{Z}_{\ge 0}$， 我們要證明 $ |x|_\mathcal{A}< \delta \Rightarrow |\phi(i;x)|_\mathcal{A} <\varepsilon $ ；故假設 $|x|_\mathcal{A}< \delta$ 則結合前述 $\delta$ 定義 我們有
\[{\left| x \right|_{{\cal A}}} < \delta  = \alpha _2^{ - 1}\left( {{\alpha _1}\left( \varepsilon  \right)} \right)\]故
\[\Rightarrow \alpha _2^{}\left( {{{\left| x \right|}_{{\cal A}}}} \right) < {\alpha _1}\left( \varepsilon  \right)\]由Lyapunov function 定義第2條不等式可知
\[V(x) \le {\alpha _2}(|x{|_{{\cal A}}}) \Rightarrow V(x) \le {\alpha _1}\left( \varepsilon  \right) \ \ \ \ (*)
\]現在觀察 Lyapunov function 定義第3條不等式，並且令 $\phi(i,x)$ 為 $x^+ = f(x)$ 之解，且注意到 $\alpha_i(\cdot)$ 其中 $i=1,2,3$ 皆為 $\mathcal{K}_\infty$ 函數，故我們可寫下
\[\begin{array}{*{20}{l}}
{V(f(x)) - V(x) \le  - {\alpha _3}\left( {|x|} \right)}\\
{ \Rightarrow V(f(x)) \le V(x)}\\
{ \Rightarrow V(f(x)) \le V(x) \le {\alpha _1}\left( \varepsilon  \right)\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array}by\begin{array}{*{20}{c}}
{}
\end{array}\left( * \right).}\\
{ \Rightarrow V(\phi \left( {i;x} \right)) \le {\alpha _1}\left( \varepsilon  \right)}\\
{ \Rightarrow {\alpha _1}(|\phi \left( {i;x} \right){|_A}) \le V(\phi \left( {i;x} \right)) \le {\alpha _1}\left( \varepsilon  \right)\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array}by\begin{array}{*{20}{c}}
{}
\end{array}1.}\\
{ \Rightarrow |\phi \left( {i;x} \right){|_A} \le \varepsilon }
\end{array}\]至此我們得到 $|x|_\mathcal{A}< \delta \Rightarrow |\phi(i;x)|_\mathcal{A} <\varepsilon $ 即為所求。

接著我們證明 global attractivity：故給定任意 $x \in X$， 要證明 $|\phi(i;x)|_\mathcal{A} \to 0 \text{ as $i \to \infty$}$。基本想法為比較不同時間點 $\phi$。 現在給定 $\phi(i;x)$ 為 $x^+ = f(x)$ 之解，由  Lyapunov function 定義第3條不等式可知
\[\begin{array}{l}
V(f(x)) - V(x) \le  - {\alpha _3}\left( {\left| x \right|} \right)\\
 \Rightarrow V(\phi \left( {i + 1;x} \right)) - V(\phi \left( {i;x} \right)) \le  - {\alpha _3}\left( {\left| {\phi \left( {i;x} \right)} \right|} \right)
\end{array}\]令 $V_{i+1}:= V(\phi \left( {i + 1;x} \right))$ 且 $V_i := V(\phi \left( {i;x} \right))$ 則由上式可推知 對任意 $x$ 而言， 數列 $\{V\}_i$ 為非遞增數列 且 有下界為 $0$，故可知 $V_i$ 收斂亦即
\[{V_{i + 1}} - {V_i} \le  - {\alpha _3}\left( {\left| {\phi \left( {i;x} \right)} \right|} \right) \to 0\begin{array}{*{20}{c}}
{}
\end{array}as\begin{array}{*{20}{c}}
{}
\end{array}i \to \infty \]故${\alpha _3}\left( {\left| {\phi \left( {i;x} \right)} \right|} \right) \to 0$ 又由於 $\alpha_3 \in \mathcal{K}_\infty$，故
\[\begin{array}{l}
\left| {\phi \left( {i;x} \right)} \right| = \alpha _3^{ - 1}\underbrace {\left( {{\alpha _3}\left( {\left| {\phi \left( {i;x} \right)} \right|} \right)} \right)}_{ \to 0}\\
 \Rightarrow \left| {\phi \left( {i;x} \right)} \right| = \alpha _3^{ - 1}\left( {{\alpha _3}\left( {\left| {\phi \left( {i;x} \right)} \right|} \right)} \right) \to 0\begin{array}{*{20}{c}}
{}&{}&{}
\end{array}since\begin{array}{*{20}{c}}
{}
\end{array}\alpha _3^{ - 1} \in {\mathcal{K}_\infty }
\end{array}\]至此證畢。 $\square$

上述定理告訴我們 asymptotically stable 的充分條件，但對於 必要條件 並無著墨。所幸透過適度的增強假設，我們仍可得到 asymptotically stable 必要條件，在此紀錄如下：

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

Theorem 2: Converse Theorem for Asymptotic Stability

令 $\mathcal{A}$ 為 compact 且 $f(\cdot)$ 為連續函數。假設 $\mathcal{A}$ 為 globally asymptotic stable for the system $x^+ = f(x)$ 則 存在 平滑 (smooth) Lyapunov function for system $x^+ = f(x)$ 與 集合 $\mathcal{A}$。
=============

Proof: omitted

故我們將 Theorem 1 與 Theorem 2 整合可得如下充分必要條件：

=============
Theorem 3 
若 $\mathcal{A}$ 為 compact 且 $f(\cdot)$ 為連續函數。則集合  $\mathcal{A}$ 為 globally asymptotic stable for the system $x^+ = f(x)$ 若且唯若  存在 平滑 (smooth) Lyapunov function for system $x^+ = f(x)$ 與 集合 $\mathcal{A}$。
=============

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

穩定度理論可追朔至 Aleksandr Lyapunov在 1892 出版的 The General Problem of Stability of Motion 提出，主要是透過建構 Lyapunov 函數 來判別動態系統是否穩定。以下討論我們將以 離散時間 非線性動態系統為主。

現在考慮以下 離散時間 非線性動態系統
\[
x^+ = f(x,u)
\]其中 $x \in \mathbb{R}^n$ 為當前系統狀態 且  $u \in \mathbb{R}^m$ 為當前的控制力；$x^+$ 為下個時刻的系統狀態。且假設 $f: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$ 為連續函數。

定義 $\phi(k; x, {\bf u}) $為在時刻 $k$，對於動態系統 $x^+ = f(x,u)$ 的解 (初始值為 $x(0)=x$ ； 控制力序列 $\bf u$ $:=\{u(0), u(1), ...\}$)

若 控制律 $u := \kappa (x)$ 決定，則系統閉迴路可表為
\[
x^+ = f(x,\kappa(x)):=f_c(x)
\]注意到 $\kappa(\cdot)$ 不一定為連續函數，此時對應的 $f(x, \kappa(\cdot))$ 亦不一定為連續。對此不連續的情況我們額外假設 $f_c(\cdot)$ 為 局部有界(locally bounded)。

目標：我們希望 控制系統 要"穩定"。

在此所謂的穩定 意指 控制系統對於 初始狀態 的小擾動 不會 導致 閉迴路系統響應 大幅度擾動 且 系統狀態能夠收斂到指定的狀態 或者 收斂到指定的 狀態集合 (此情況多半發生在有外部干擾的時候)。


以下我們會針對定義 系統的 穩定度 與 漸進穩定度；在介紹之前我們需要先定義一些名詞：首先是 如何指出系統狀態的收斂

============
Definition: Equilibrium Point or Steady-State
狀態 $x^*$ 被稱作 $x^+ = f(x)$ 的 平衡點(equilibrium point) 若 $$x(0) = x^* \Rightarrow x(k) = \phi(k;x^*) = x^*, \;\; \forall k \ge 0$$
============

Comment:
1. 上述定義表示 $x^*$ 為 平衡點 若其滿足 $x^* = f(x^*)$
2. equilibrium point 為 被隔離的(isolated) 若在 $x^*$ 附近沒有其他的平衡點。
3. 非線性系統可能有多個 被隔離的平衡點


Example: Equilibrium Point of Linear System
考慮離散時間線性系統
\[
x^+  = Ax +b
\]具有平衡點 $x^*$ 則
\[\begin{array}{l}
{x^ + } = Ax + b\
 \Rightarrow {x^*} = A{x^*} + b\\
 \Rightarrow \left( {I - A} \right){x^*} = b
\end{array}\]若 $I-A$ 反矩陣存在，則我們說 此線性系統有 unique (isolated) 平衡點
\[{x^*} = {\left( {I - A} \right)^{ - 1}}b\]若 $I-A$ 反矩陣不存在，則我們說此線性系統有 連續統 (continuum) $\{x: (I-A) x=b\}$ 的平衡點。


若我們考慮 震盪系統 的穩定度，則此時不再是討論 是否收斂到某個狀態 (平衡點)；而是討論收斂到某個集合。以下我們給出此類集合所需的定義：

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

Definition: Positive Invariant Set
一個集合 $\mathcal{A}$ 稱作 positive invariant for system $x^+ = f(x)$ 若下列條件成立：
\[
x \in A \Rightarrow f(x) \in \mathcal{A}
\]=============

Comment:

1. Positive 來自於 $x^+ = f(x)$ 為動態系統隨時間 $k$ "增加" 而持續變動。
2. 考慮 closed set $\mathcal{A}:=\{x^*\}$ 且 $x^*$ 為系統 $x^+ = f(x)$ 的平衡點，則

\[
x \in \mathcal{A}\;\; (\text{since} \;x^* \in \mathcal{A}) \Rightarrow f(x) \in \mathcal{A} \;\;(\text{since} \;f(x) = x^*)
\]

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

Definition: K, K infinity, KL function
一個函數 $g: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$ 為 $\mathcal{K}$ 類函數若下列條件滿足：

1. $g$ 為連續
2. $g(0) = 0$
3. 嚴格遞增(strictly increasing)；亦即 $\forall x,y$，$y > x \Rightarrow g(y) > g(x)$

我們說 一個函數 $g: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$ 為 $\mathcal{K}_\infty$ 類函數若下列條件滿足：

1. $g$ 為 $\mathcal{K}$ 類函數
2. 當 $t \to \infty$，$g(t) \to \infty$

我們說一個函數 $h: \mathbb{R}_{\ge 0} \times \mathbb{Z}_{\ge 0} \to \mathbb{R}_{\ge 0}$ 為 $\mathcal{KL}$ 類函數若下列條件滿足：

1. 對任意 $t \ge 0$，$h(\cdot, t)$ 為 $\mathcal{K}$ 類函數
2. 對任意 $s \ge 0$，$h(s, \cdot)$ 為非遞增(nonincreasing) 且 滿足 $\lim_{t\to \infty} h(s,t) =0$

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

Example

1. $g(x) := x$  為 $\mathcal{K}$類函數 (亦為 $\mathcal{K}_\infty $ 函數)

2. $erf(x)$ 為  $\mathcal{K}$類函數

以下我們將前述 $\mathcal{K}$類函數的重要性質：

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

FACT 1: Inverse K function is a K function
若 $\alpha_1(\cdot), \alpha_2(\cdot)$ 為 $\mathcal{K}$ 類函數 (或者 $\mathcal{K}_\infty$ 函數)，則其反函數 $\alpha_1^{-1}(\cdot), \alpha_1^{-1}(\cdot)$ 亦仍為 $\mathcal{K}$ 類函數 (或者 $\mathcal{K}_\infty$ 函數)

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

Proof: omitted

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

FACT 2: 
若 $\alpha_1(\cdot)$ 與 $\alpha_2(\cdot)$ 為  $\mathcal{K}$ 類函數 且 $\beta(\cdot)$ 為 $\mathcal{KL}$ 函數，則 $\sigma (r,s): = {\alpha _1}(\beta \left( {{\alpha _2}\left( r \right)} \right),s)$ 為 $\mathcal{KL}$函數。

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

Proof: omitted


有了以上定義我們可以開始引入 穩定度 的嚴格定義。以下我們考慮 $x^+ = f(x)$ 且假設 $f(\cdot)$ 為 局部有界(locally bounded) 且集合 $A$ 為 closed 與 positive invariant 。

==================
Definition: Local Stability (Stability in Lyapunov Sense)
給定 closed positive invariant 集合 $\mathcal{A}$ 。我們稱 此集合 $\mathcal{A}$ 為 locally stable for $x^+ = f(x)$ 若下列條件成立：
對任意 $\varepsilon>0$ 存在 $\delta >0$ 使得對任意 $i \in Z_{\ge 0}$， $|x|_\mathcal{A} < \delta \Rightarrow |\phi(i; x)|_\mathcal{A} < \varepsilon$
其中 $|x|_\mathcal{A} := \inf_{z \in \mathcal{A}} |x - z|$ 

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

Example:
考慮 $A:= \{0\}$ 則 Local Stability 可由下圖得知

上圖顯示了若給定任意初始位置 $x$ 且此 $x$ 與原點 $A:=\{0\}$ 距離落在 開球 $B_{\delta}$ 之中，且若系統 $x^+ = f(x)$ 隨時間變化演進，其解 $\phi(i,x)$ 到原點距離 $A=\{0\}$ 持續落在另一開球 $B_\varepsilon$之中，故此系統稱為 Local stable。

======================
Definition: Global Attraction
給定 closed positive invariant 集合 $\mathcal{A}$ 。我們說此集合 $A$ 為 globally attractive for system $x^+ = f(x)$ 若下列條件成立：
\[
|\phi(i;x)|_\mathcal{A} \to 0 \text{ as $i \to \infty$} \;\; \forall x\in \mathbb{R}^n
\]======================

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

Definition: (Global Asymptotic Stability (GAS))

給定  closed positive invariant 集合 $\mathcal{A}$ 為 globally asymptotically stable for system $x^+ = f(x)$ 若下列條件成立：

$\mathcal{A}$ 為 locally stable 且 globally attractive
======================

Comment:
考慮 $\mathcal{A}:=\{0\}$，有可能 globally attractive 但並非 locally stable。比如說考慮\[
x^+ = Ax + \phi(x)
\]其中 $A$ 有 eigenvalue $\lambda_1 = 0.5$ 與 $\lambda_2 = 2$ 且對應的 eigenvector 為 $w_1, w_2$且 $\phi(\cdot)$為 平滑函數 滿足 $\phi(0) = 0$ 與 ${\left. {\frac{\partial }{{\partial x}}\phi (x)} \right|_{x = 0}} = 0$

故在 $0$ 附近，$x^+ = Ax + \phi(x)$ 行為將會非常接近 $x^+ = Ax$ ；故若 $\phi(x) =0$ 則 特徵向量$w_1$ 因為具有特徵值 $\lambda_1 = 0.5$ (落在 unit circle 之中)故此特徵向量會迫使狀態收斂到 $0$點，但 特徵向量 $w_2$ 具有不穩定的特徵值 $\lambda_2 = 2$ 故此特徵向量會迫使狀態發散。故總和此兩者，可知儘管有 globally attractive 但卻沒有 stable origin ($\mathcal{A}:=\{0\}$無法滿足 local stability 定義)


以下我們將相關的穩定度定義總結如下：
=============================
Definition: Stability without constraint
給定 closed positive invariant 集合 $\mathcal{A}$ 為

1. locally stable 若 對任意 $\varepsilon >0$ 存在 $\delta >0$ 使得對任意 $i \in \mathbb{Z}_{\ge 0}$ $|x|_\mathcal{A} < \delta \Rightarrow |\phi(i;x)|_\mathcal{A} <\varepsilon $
2. unstable 若 其 不為 locally stable
3. locally attractive 若 存在 $\eta >0$ 使得 $|x|_\mathcal{A} < \eta \Rightarrow |\phi(i;x)|_\mathcal{A} \to 0\;\; \text{as $i \to \infty$} $
4. globally attractive 若對任意 $x \in \mathbb{R}^n$， $|\phi(i;x)|_\mathcal{A} \to 0\;\; \text{as $i \to \infty$} $ 
5. locally asymptotically stable 若其為 locally stable 與 locally attractive
6. globally asymptotically stable 若其為 locally stable 與 globally attractive
7. locally exponentially stable 若 存在 $\eta >0, c>0$ 與 $\gamma \in (0,1)$ 使得 對任意 $i \in \mathbb{Z}_{\ge 0}$ 而言， $|x|_\mathcal{A} < \eta \Rightarrow |\phi(i;x)|_{\mathcal{A}} \le c |x|_\mathcal{A} \gamma^i$
8. globally exponentially stable 若 存在 $c>0$ 與 $\gamma \in (0,1)$ 使得 對任意 $i \in \mathbb{Z}_{\ge 0}$ 而言， $|\phi(i;x)|_{\mathcal{A}} \le c |x|_\mathcal{A} \gamma^i$

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

以下結果將前述穩定度定義 與 KL 函數做連結：
=============================
FACT: (Globally Asymptotic Stable and KL function)
令集合 $\mathcal{A}$ 為 compact 且 positive invariant； $f(\cdot)$ 為連續函數。則 $\mathcal{A}$ 為 globally asymptotic stable for $x^+ = f(x)$ 若且唯若 存在 $\mathcal{KL}$ 函數 $\beta(\cdot)$ 使得 對任意 $x \in \mathbb{R}^n$，
\[
|\phi(i;x)|_\mathcal{A} \le \beta(|x|_\mathcal{A},i)\;\; \forall i \in \mathbb{Z}_{\ge 0}
\]=============================

另外，實際上若考慮系統狀態有拘束的情形，則 globally asymptotic stability 並不保證能夠達成，此時我們需要再次拓展前述定義來滿足有拘束的情況：

=============================
Definition: Stability with Constraint Set X
假設狀態拘束集合  $X \subset \mathbb{R}^n$ 為 positive invariant for $x^+ = f(x)$ 且集合 $\mathcal{A}$  closed positive invariant for $x^+ = f(x)$  且 $\mathcal{A} \subset int(X)$  ( $int(X) :=$ interior of $X$) 則我們說集合 $\mathcal{A}$ 為

1. locally stable in $X$ 若 對任意 $\varepsilon >0$ 存在 $\delta >0$ 使得對任意 $i \in \mathbb{Z}_{\ge 0}$ $x \in X \cap (\mathcal{A} \oplus  B_\delta) \Rightarrow |\phi(i;x)|_\mathcal{A} <\varepsilon $
2. locally attractive in $X$ 若 存在 $\eta >0$ 使得 $x \in X \cap (\mathcal{A} \oplus  B_\delta) \Rightarrow |\phi(i;x)|_\mathcal{A} \to 0\;\; \text{as $i \to \infty$} $
3. attractive in $X$ 若 $ |\phi(i;x)|_\mathcal{A} \to 0\;\; \text{as $i \to \infty$} \; \forall x \in X$
4. locally asymptotically stable in $X$ 若其為 locally stable in $X$ 與 locally attractive in $X$
5. asymptotically stable  with region of attraction $X$ 若其為 locally stable in $X$ 與 attractive in $X$。
6. locally exponentially stable with region of attraction $X$ 若 存在 $\eta >0, c>0$ 與 $\gamma \in (0,1)$ 使得 對任意 $i \in \mathbb{Z}_{\ge 0}$ 而言， $x \in X \cap (\mathcal{A} \oplus  B_\eta) \Rightarrow |\phi(i;x)|_{\mathcal{A}} \le c |x|_\mathcal{A} \gamma^i$
7. globally exponentially stable with region of attraction $X$ 若 存在 $c>0$ 與 $\gamma \in (0,1)$ 使得 對任意 $i \in \mathbb{Z}_{\ge 0}$ 而言， $|\phi(i;x)|_{\mathcal{A}} \le c |x|_\mathcal{A} \gamma^i$

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

延伸閱讀
[系統理論] 離散時間系統的穩定度理論 (1) - Lyapunov Stability Theory


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

[控制理論] 離散線性系統的 追蹤 與 調節 問題

Setpoint Tracking
考慮離散線性系統
\[\left\{ \begin{array}{l}
x\left( {k + 1} \right) = Ax\left( k \right) + Bu\left( k \right)\\
y\left( k \right) = Cx\left( k \right)
\end{array} \right.\] $x \in \mathbb{R}^n, y \in \mathbb{R}^p, u \in \mathbb{R}^m$；一般而言，控制系統中常見的 追蹤(setpoint tracking)問題 (或稱 servo problem) 如下：
給定 setpoint  $y_{sp}$ (e.g., 步階訊號 or 常數)，我們希望在系統達到穩態(steady-state) 之後，系統的輸出 $ y  = y_{sp} $。

若我們考慮給定的 setpoint $y_{sp} = 0$ 則稱此類問題為 regulation problem。

回憶在最佳控制理論中的 LQR 方法給予我們對於線性系統的 regulation 問題提供一組最佳解，故我們的想問是否能將此法應用在 一般的追蹤問題?

答案是肯定的，僅需引入 deviation variable 做基本 座標轉換 即可。

現在我們令 $y_{sp}$ 為 output setpoint，定義系統穩態時候的狀態 與 控制力為 $(x_s, u_s)$ 則線性系統方程可改寫為
\[\left\{ \begin{array}{l}
{x_s} = A{x_s} + B{u_s}\\
{y_s} = C{x_s}
\end{array} \right.\]注意到對於穩態時 我們希望 $y_s = y_{sp}$ 故 我們有 $y_{sp} = C x_s$ 現在我們將上式改寫為矩陣形式
\[\left\{ \begin{array}{l}
{x_s} = A{x_s} + B{u_s}\\
{y_{sp}} = C{x_s}
\end{array} \right. \Rightarrow \left[ {\begin{array}{*{20}{c}}
{I - A}&{ - B}\\
C&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_s}}\\
{{u_s}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
0\\
{{y_{sp}}}
\end{array}} \right] \ \ \ \ (*)
\]若上述矩陣方程有解 (亦即可以解出 $x_s, u_s$)，且系統控制力無拘束，則我們可以定義 deviation variables 如下
\[\left\{ \begin{array}{l}
\tilde x\left( k \right): = x\left( k \right) - {x_s}\\
\tilde u\left( k \right): = u\left( k \right) - {u_s}
\end{array} \right.\]現在觀察
\[\begin{array}{l}
\tilde x\left( {k + 1} \right) = x\left( {k + 1} \right) - {x_s}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = Ax\left( k \right) + Bu\left( k \right) - \left( {A{x_s} + B{u_s}} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = A\left( {x\left( k \right) - {x_s}} \right) + B\left( {u\left( k \right) - {u_s}} \right)\\
 \Rightarrow \tilde x\left( {k + 1} \right) = A\tilde x\left( k \right) + B\tilde u\left( k \right)
\end{array}\]此表示我們的 deviation variable 仍然滿足原本給定的線性系統。且此時若考慮使用 deviation variable $\tilde{x}, \tilde{u}$ 則 我們的 setpoint tracking problem 被改寫為 regulation problem ；也就是說我們要找到一組 $\tilde{u}(k)$ 使得 $\tilde x(k) \to 0$ (此等價為 $x(k) \to x_s$，故 $C x(k) \to C x_s = y_{sp}$)。解完此 regulation problem 之後，真實的控制力 $u(k) = \tilde u(k) + u_s$。

Comment:
注意到上述論述建立在 $(n + p) \times (n + m)$ 的矩陣方程
\[\left\{ \begin{array}{l}
{x_s} = A{x_s} + B{u_s}\\
{y_{sp}} = C{x_s}
\end{array} \right. \Rightarrow \left[ {\begin{array}{*{20}{c}}
{I - A}&{ - B}\\
C&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_s}}\\
{{u_s}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
0\\
{{y_{sp}}}
\end{array}} \right] \ \ \ \ (*)
\]有解。但何時才有解/與此解是否唯一的問題 我們並未解決 !! 以下我們將對此點進行討論。


對任意 setpoint $y_{sp}$， 上述矩陣方程 $(*)$ 有解 的充分條件為：$(*)$ 要有 linearly independent rows；亦即 $p \le m$ 也就是說我們需要 控制力的數目 $m$ 與 量測輸出的數目 $p$ 至少要相等)。但在實際情況上卻是常常相反，我們可能會獲得非常多的量測輸出數目 (因為裝了很多 sensor)，但實際可以調控的變數卻少於 sensor 數目。為了要解決此問題，我們可選定某矩陣 $H$ 並引入新的 控制變數 $r \in \mathbb{R}^{n_c}$ 為 量測輸出的線性組合；亦即
\[
r := Hy
\]如此一來，若 $p > m$ 情況發生時，我們可選一部份的 輸出 $n_c \le m$ 作為控制變數 並且 對此引入的控制變數 $r$ 也給定所需的 setpoint $r_{sp}$。

另外若 $m > p$ (控制力數目 大於 量測輸出的數目) 時，則對某些 $H$ 與 $r_{sp}$，矩陣方程 $(*)$ 有解 但此時唯一性並不被保證，故我們會希望有唯一解，此時需要對 穩態控制力  $u_s$ 也給定 setpoint $u_{sp}$。


總和以上所述，我們可以建構以下 Steady-State Target Problem
========================
Steady-State Target Problem
考慮最佳化問題
\[
\min_{x_s, u_s} \frac{1}{2} (|u_s - u_{sp}|_{R_s}^2 + |y_s - y_{sp}|_{Q_s}^2)
\]subject to
\[\left\{ \begin{array}{l}
{x_s} = A{x_s} + B{u_s}\\
{y_{sp}} = C{x_s}
\end{array} \right. \Rightarrow \left[ {\begin{array}{*{20}{c}}
{I - A}&{ - B}\\
C&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_s}}\\
{{u_s}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
0\\
{{y_{sp}}}
\end{array}} \right]
\]與 $E u_s \le e$ 與 $F C x_s \le f$。其中 $R_s$ 為 正定矩陣。且對 控制變數的 setpoints $r_{sp}$，上述 target problem 為 feasible 。
========================

有了 steady-state 的解之後我們可以便可以求解 當初我們引入 deviation variable 的 regulation problem:
========================
Dynamic Regulation Problem
考慮以下 cost function
\[\begin{array}{l}
V\left( {\tilde x\left( 0 \right),{\bf{\tilde u}}} \right) = \frac{1}{2}\sum\limits_{k = 0}^{N - 1} {\left| {\tilde x\left( k \right)} \right|_Q^2 + \left| {\tilde u\left( k \right)} \right|_R^2} \\
s.t.\\
\tilde x\left( {k + 1} \right) = A\tilde x\left( k \right) + B\tilde u\left( k \right)
\end{array}\]其中 $\tilde x(0) = \hat x(k) - x_s$ (此初始條件表示  透過 steady-state $x_s$ 平移估計狀態 $\hat x$ 而得。)；且所求得的 regulator 將會求解以下的 regulation problem
 \[\begin{array}{l}
\mathop {\min }\limits_{{\bf{\tilde u}}} V\left( {\tilde x\left( 0 \right),{\bf{\tilde u}}} \right)\\
s.t.\\
E\tilde u \le e - E{u_s}\\
FC\tilde x \le f - FC{x_s}
\end{array}\]===================
上述 regulation problem 的 optimal cost 為 $V^*(\tilde x(0))$ 且對應的控制力為 $\tilde {u}^0 (\tilde{x}(0))$

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

[最佳控制] 線性系統的最佳參數估計 - Kalman Filter (1)

此文章基於前篇 [最佳控制] 線性系統的最佳參數估計 - Kalman Filter (0)；強烈建議讀者先參閱前篇再行閱讀此文。此文將利用 計算 probability density 方式來求解 Kalman filter。


考慮離散時間動態系統
\[
x(k+1) = Ax(k) + w(k);\;\; y(k) = Cx(k) + v(k)
\]且 $w \sim N(0,Q)$ 為製程雜訊 process noise 或稱 干擾 process disturbance；$v \sim N(0,R)$ 為量測雜訊 (measurement noise) ； $x(0) \sim N(\bar x(0), Q(0))$；$x \in \mathbb{R}^n, A \in \mathbb{R}^{n \times n}, C \in \mathbb{R}^{p \times n}, y \in \mathbb{R}^p$，Kalman filter 便是在試圖回答：假設狀態未知我們只能拿到量測輸出，則最佳的狀態估計該是如何?


First Step : $k=0$
假設 初始狀態 $x(0)$ 為 mean $\bar x(0)$ 且 convariance matrix $Q(0)$ 的 normal distrbuted 隨機向量，亦即
\[
x(0) \sim N(\bar x(0), Q(0))
\]接著獲得 初始 (受雜訊污染的) 量測輸出 $y(0)$ 滿足下式
\[
y(0) = C x(0) + v(0)
\]其中 $v(0) \sim N(0, R)$ 為 雜訊 (measurement noise)。

我們的目標：獲得 conditional density $p_{x(0)|y(0)}(x(0) | y(0))$ ，則我們的狀態估計 $\hat x$ 即可透過此 conditional density 求得
$$\hat x := \arg \max_x p_{x(0)|y(0)}(x(0) | y(0))$$

Comments:
1. 若未知 $\bar x(0)$ 或者 $Q(0)$ 則我們通常選 $\bar x(0) :=0$ 且 $Q(0)$ 很大 來表示我們對初始狀態所知甚少 (noninformative prior)。
2. 若量測過程中受到大的雜訊，則我們會給予較大的 $R$。若量測過程十分精確沒有太多雜訊汙染我們的輸出 $y$ 則 $R$ 較小。
3. Conditional density 描述了在我們獲得 初始量測輸出 $y(0)$ 之後我們對於 $x(0)$ 的了解。


現在回歸我們的目標，究竟該如何推得 $p_{x(0)|y(0)}(x(0) | y(0))$ ?
首先考慮 $(x(0), y(0))$ 如下
\[\left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{y\left( 0 \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{v\left( 0 \right)}
\end{array}} \right]
\]假設 $v(0)$ 與 $x(0)$ 彼此互為獨立。注意到  $x(0), v(0)$ 為 joint normal 亦即
\[\left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{v\left( 0 \right)}
\end{array}} \right] \sim N\left( {\left[ {\begin{array}{*{20}{c}}
{\bar x\left( 0 \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{Q\left( 0 \right)}&0\\
0&{{R}}
\end{array}} \right]} \right)\]
上述 $[x(0) \;\; y(0)]$ 為  $[x(0) \;\; v(0)]$ 線性轉換 且 故我們可以馬上知道
\[\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{y\left( 0 \right)}
\end{array}} \right]\sim N\left( {\left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{\bar x\left( 0 \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{Q\left( 0 \right)}&0\\
0&R
\end{array}} \right]{{\left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]}^T}} \right)\\
 \Rightarrow \left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{y\left( 0 \right)}
\end{array}} \right] \sim N\left( {\left[ {\begin{array}{*{20}{c}}
{\bar x\left( 0 \right)}\\
{C\bar x\left( 0 \right)}
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{Q\left( 0 \right)}&0\\
0&R
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
I&{{C^T}}\\
0&I
\end{array}} \right]} \right)\\
 \Rightarrow \left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{y\left( 0 \right)}
\end{array}} \right]\sim N\left( {\left[ {\begin{array}{*{20}{c}}
{\bar x\left( 0 \right)}\\
{C\bar x\left( 0 \right)}
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{Q\left( 0 \right)}&{Q\left( 0 \right){C^T}}\\
{CQ\left( 0 \right)}&{CQ\left( 0 \right){C^T} + R}
\end{array}} \right]} \right)
\end{array}\]現在有了 $x(0), y(0)$ 的 joint density， 注意到此 joint density 為 normal，故若要計算 conditonal density of $x(0)$ given $y(0)$，亦即 $p_{x(0)|y(0)}(x(0)|y(0))$ ；則我們可以使用前述文章討論的 FACT 3 來求得，亦即
\[
p_{x(0)|y(0)}(x(0)|y(0)) = n(x(0), m, P)
\]其中
\[\left\{ {\begin{array}{*{20}{l}}
\begin{array}{l}
m = \bar x\left( 0 \right) + L\left( 0 \right)(y\left( 0 \right) - C\bar x\left( 0 \right))\\
L\left( 0 \right): = Q\left( 0 \right){C^T}\left( {CQ\left( 0 \right){C^T} + R} \right)_{}^{ - 1}
\end{array}\\
{P = Q\left( 0 \right) - Q\left( 0 \right){C^T}\left( {CQ\left( 0 \right){C^T} + R} \right)_{}^{ - 1}CQ\left( 0 \right)}
\end{array}} \right.\]則 optimal state estimation  $\hat x$ 及為 具有最大 conditional density 的 $x(0)$；對於 normal distribution 而言，此 $\hat x$ 剛好為 mean；故我們選 $\hat x(0) := m$；且上式中 $P := P(0)$ 表示獲得 量測輸出 $y(0)$ 之後的 variance

Next Step: State Prediction
考慮現在狀態從 $k=0$ 移動到 $k=1$，我們有
\[
x(1) = Ax(0) + w(0)
\]亦可將其寫成線性轉換型式：
\[x\left( 1 \right) = \left[ {\begin{array}{*{20}{c}}
A&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{w\left( 0 \right)}
\end{array}} \right]\]其中 $w(0) \sim N(0, Q)$ 為 系統干擾 (disturbance) 或稱 製程雜訊 (process noise)。若 狀態受到大的外部擾動，則可以想見 會有較大的 $Q$ convaraince matrix。同理若外部干擾較小則有較小的雜訊。

故我們目標：要計算 conditional density $p_{x(1)|y(0)}(x(1),y(0))$

現在我們需要 conditional on joint density $(x(0),w(0))$ given $y(0)$，假設 $w(0)$ 與 $x(0),v(0)$ 彼此互為獨立，則我們可寫
\[\left( {\left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{w\left( 0 \right)}
\end{array}} \right]|y\left( 0 \right)} \right)\sim N\left( {\left[ {\begin{array}{*{20}{c}}
{\hat x\left( 0 \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{P\left( 0 \right)}&0\\
0&Q
\end{array}} \right]} \right)\]故現在利用 前述文章討論的 FACT 3' 來求得  conditional density，亦即
 \[\begin{array}{l}
\left( {\left[ {\begin{array}{*{20}{c}}
{x\left( 0 \right)}\\
{w\left( 0 \right)}
\end{array}} \right]|y\left( 0 \right)} \right)\sim  N\left( {\left[ {\begin{array}{*{20}{c}}
{\hat x\left( 0 \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{P\left( 0 \right)}&0\\
0&Q
\end{array}} \right]} \right)\\
 \Rightarrow \left( {x\left( 1 \right)|y\left( 0 \right)} \right)\sim N\left( {\left[ {\begin{array}{*{20}{c}}
A&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{\hat x\left( 0 \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
A&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{P\left( 0 \right)}&0\\
0&Q
\end{array}} \right]{{\left[ {\begin{array}{*{20}{c}}
A&I
\end{array}} \right]}^T}} \right)\\
 \Rightarrow \left( {x\left( 1 \right)|y\left( 0 \right)} \right)\sim N\left( {A\hat x\left( 0 \right),AP\left( 0 \right){A^T} + Q} \right)
\end{array}\]故 conditional density 仍為 normal
\[
p_{x(1)|y(0)}(x(1)|y(0)) = n(x(1), \hat x^-(1), P^-(1))
\]其中
\[\left\{ \begin{array}{l}
{{\hat x}^ - }\left( 1 \right) = A\hat x\left( 0 \right)\\
{P^ - }\left( 1 \right) = AP\left( 0 \right){A^T} + Q
\end{array} \right.\]接著我們僅需 遞迴重複上述步驟 $k=2,3,4,...$ 即可。以下我們給出總結：

Summary
定義 量測輸出從初始 直到 時間 $k$ 則
\[
{\bf y}(k) := \{y(0),y(1),...,y(k)\}
\]在 時間 $k$ 時，conditional density with data ${\bf y}(k-1)$ 為 normal
\[
p_{x(k)| {\bf y}(k-1)}(x(k) | {\bf y}(k-1)) = n(x(k), \hat x^-(k),  P^-(k))
\]上述 mean 與 covariance matrix 有上標 $^-$ 號表示此估計為僅僅透過 量測 ${\bf y}(k-1)$ 的輸出，並未獲得 $k$ 時刻的量測輸出。 (表示用過去 $k-1$ 資料 預測 $k$ 的狀態! ) 注意到在 $k=0$，遞迴起始於 $\hat x^- (0) = \bar x(0)$ 與 $P^-(0) = Q(0)$。

接著我們會獲得 $y(k)$ 滿足
\[\left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{y\left( k \right)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{v\left( k \right)}
\end{array}} \right]
\]上式表示線性轉換。由於 量測雜訊 $v(k)$ 與 $x(k)$ 以及 ${\bf y}(k-1)$ 彼此獨立，故我們有  density of $(x(k), v(k))$
\[\underbrace {\left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{v\left( k \right)}
\end{array}} \right]}_{ = \left( {\left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{v\left( k \right)}
\end{array}} \right]|{\bf{y}}\left( {k - 1} \right)} \right)}\sim N\left( {\left[ {\begin{array}{*{20}{c}}
{{{\hat x}^ - }\left( k \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{P^ - }\left( k \right)}&0\\
0&R
\end{array}} \right]} \right)\]透過線性轉換結果，可得
\[\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{y\left( k \right)}
\end{array}} \right]\sim N\left( {\left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{{\hat x}^ - }\left( k \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{P^ - }\left( k \right)}&0\\
0&R
\end{array}} \right]{{\left[ {\begin{array}{*{20}{c}}
I&0\\
C&I
\end{array}} \right]}^T}} \right)\\
 \Rightarrow \left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{y\left( k \right)}
\end{array}} \right]\sim N\left( {\left[ {\begin{array}{*{20}{c}}
{{{\hat x}^ - }\left( k \right)}\\
{C{{\hat x}^ - }\left( k \right)}
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{{P^ - }\left( k \right)}&{{P^ - }\left( k \right){C^T}}\\
{C{P^ - }\left( k \right)}&{C{P^ - }\left( k \right){C^T} + R}
\end{array}} \right]} \right)
\end{array}\]注意到 $\{{\bf y}(k-1), y(k)\} \equiv {\bf y}(k)$ 故利用 conditional density 結果可得
\[
p_{x(k)|{\bf y} (k)}(x(k)| {\bf y}(k)) = n(x(k), \hat x(k), P(k))
\]其中
\[\left\{ {\begin{array}{*{20}{l}}
\begin{array}{l}
\hat x\left( k \right) = {{\hat x}^ - }\left( k \right) + L\left( k \right)(y\left( k \right) - C{{\hat x}^ - }\left( k \right))\\
L\left( k \right) = {P^ - }\left( k \right){C^T}\left( {C{P^ - }\left( k \right){C^T} + R} \right)_{}^{ - 1}
\end{array}\\
{P\left( k \right) = {P^ - }\left( k \right) - {P^ - }\left( k \right){C^T}\left( {C{P^ - }\left( k \right){C^T} + R} \right)_{}^{ - 1}C{P^ - }\left( k \right)}
\end{array}} \right.\]現在我們用下列模型來預估 基於 $k$ 時刻量測值 預估 $k+1$ 時刻
\[x\left( {k + 1} \right) = \left[ {\begin{array}{*{20}{c}}
A&I
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{w\left( k \right)}
\end{array}} \right]\]由於 $w(k)$ 與 $x(k)$ 以及 $\bf y$$(k)$ 彼此獨立，故 joint density of $(x(k),w(k))$ 可寫為
\[\left[ {\begin{array}{*{20}{c}}
{x\left( k \right)}\\
{w\left( k \right)}
\end{array}} \right]\sim N\left( {\left[ {\begin{array}{*{20}{c}}
{\hat x\left( k \right)}\\
0
\end{array}} \right],\left[ {\begin{array}{*{20}{c}}
{P\left( k \right)}&0\\
0&Q
\end{array}} \right]} \right)\]故 conditional density 為
\[{p_{x\left( {k + 1} \right)|{\bf{y}}\left( k \right)}}\left( {x\left( {k + 1} \right)|{\bf{y}}\left( k \right)} \right) = n\left( {x\left( {k + 1} \right),{{\hat x}^ - }\left( {k + 1} \right),{P^ - }\left( {k + 1} \right)} \right)\]其中
\[\left\{ \begin{array}{l}
{{\hat x}^ - }\left( {k + 1} \right) = A\hat x\left( k \right)\\
{P^ - }\left( {k + 1} \right) = AP\left( k \right){A^T} + Q
\end{array} \right.\]

[最佳控制] 線性系統的最佳參數估計 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".