[微分方程] 人口增長模型

考慮以下初始值問題
\[
\frac{dy}{dt} = ky - cy^2, \;\; y(0) = A \in \mathbb{R}
\]上述微分方程用以描述 具備外在競爭關係 $(-cy^2)$ 的 (人口)成長模型 $(ky)$。一般稱之為 Logistic Equation of Population Growth

Comments:
1. 上述 Logistic Equation 具備 $y' = f(y)$ 形式，其中 $f(y) := ky - cy^2$ ，亦即函數 $f$ 僅與 $y$ 有關而無直接與時間 $t$相關，我們稱此類微分方程為 autonomous equation。


現在我們開始求解 Logistic Equation

Solution:
令 $\phi(t) \neq 0$ 為上述 IVP 在某區間 $I$ 包含 $t=0$ 的解，則在此區間 $I$ 上，我們有
\[
\phi '\left( t \right) = k\phi \left( t \right) - c{\phi ^2}\left( t \right) \;\;\;\;\; (*)
\]注意到儘管上式為非線性，但其具備分離變數形式，故若
\[\begin{array}{l}
k\phi \left( t \right) - c{\phi ^2}\left( t \right) \ne 0\\
 \Rightarrow \left( {k - c\phi \left( t \right)} \right)\phi \left( t \right) \ne 0
\end{array}
\]亦即 $\phi \left( t \right) \ne 0 $ 且 ${k - c\phi \left( t \right)} \neq 0$(或者等價 $\phi(t) \neq k/c$) 則我們可以將 $(*)$ 改寫為
\[\begin{array}{l}
\phi '\left( t \right) = k\phi \left( t \right) - c{\phi ^2}\left( t \right)\\
 \Rightarrow \frac{{\phi '\left( t \right)}}{{k\phi \left( t \right) - c{\phi ^2}\left( t \right)}} = 1
\end{array}
\]對上式兩邊同積分從 $0$ 到 $t$ 可得
\[\int_0^t {\frac{{\phi '\left( s \right)}}{{k\phi \left( s \right) - c{\phi ^2}\left( s \right)}}ds}  = \underbrace {\int_0^t 1 dt}_{ = t}
\]現在令變數變換 $u:= \phi(s)$ 且使用初始條件 $\phi(0) = A$ 我們可將上述積分改寫如下
\[\begin{array}{l}
\int_A^{\phi \left( t \right)} {\frac{1}{{ku - c{u^2}}}} du = t\\
 \Rightarrow \int_A^{\phi \left( t \right)} {\frac{1}{{u\left( {k - cu} \right)}}} du = t \;\;\; (**)
\end{array}
\]注意到左式被積函數可透過部分分式展開如下
\[\frac{1}{{u\left( {k - cu} \right)}} = \frac{C}{u} + \frac{D}{{k - cu}}
\]其中 $C = 1/k$ 且 $D =c/k$ 故 $(**)$ 式子改寫為
\[\begin{array}{l}
\frac{1}{k}\int_A^{\phi \left( t \right)} {\frac{1}{u}} du + \frac{c}{k}\int_A^{\phi \left( t \right)} {\frac{1}{{k - cu}}} du = t\\
 \Rightarrow \frac{1}{k}\left( {\ln \left| {\phi \left( t \right)} \right| - \ln |A|} \right) + \left( {\frac{1}{{ - k}}} \right)\left( {\ln \left| {k - c\phi \left( t \right)} \right| - \ln \left| {k - cA} \right|} \right) = t\\
 \Rightarrow \frac{1}{k}\left( {\ln \left| {\frac{{\phi \left( t \right)}}{A}} \right|} \right) - \frac{1}{k}\left( {\ln \left| {\frac{{k - c\phi \left( t \right)}}{{k - cA}}} \right|} \right) = t\\
 \Rightarrow \frac{1}{k}\ln \left| {\frac{{\phi \left( t \right)\left( {k - cA} \right)}}{{A\left( {k - c\phi \left( t \right)} \right)}}} \right| = t
\end{array}
\]由上述結果不難求解 $\phi(t)$ 如下
\[\begin{array}{l}
\phi \left( t \right) = \frac{{kA{e^{kt}}}}{{\left( {k - cA} \right) + cA{e^{kt}}}}\\
 \Rightarrow \phi \left( t \right) = \frac{{kA}}{{\left( {k - cA} \right){e^{ - kt}} + cA}}
\end{array}\]自此我們求解完畢，讀者可自行驗證上述解 $\phi(t)$ 確實滿足 Logistic Equaiton 與初始條件。

Comments:
1. 若 $c = 0$，我們可得
\[{\left. {\phi \left( t \right)} \right|_{c = 0}} = {\left. {\frac{{kA}}{{\left( {k - cA} \right){e^{ - kt}} + cA}}} \right|_{c = 0}} = A{e^{kt}}\]亦即滿足 $y' = ky$ 之解。

2. 若 $c \neq 0$ 且 $A \neq 0$ 則我們可推得系統其中的一個穩態解
\[\mathop {\lim }\limits_{t \to \infty } \phi \left( t \right) = \mathop {\lim }\limits_{t \to \infty } \frac{{kA}}{{\left( {k - cA} \right){e^{ - kt}} + cA}} = \frac{k}{c}\]
3. 若 初始條件 $y(0) = A = 0$ 則 $\phi(t) = 0$ 為第二個穩態解

4. 一般我們稱邏輯函數 (Logistic Function) 為
\[L\left( t \right) = \frac{1}{{1 + {e^{ - t}}}}\]讀者可驗證其函數圖形具有 $S$ 形。現在注意到我們上述所求得的解 $\phi(t)$ 形式具備此邏輯函數的形式，故其 ODE
\[
\frac{dy}{dt} = ky - cy^2
\]我們稱其為 Logistic Equation (of Population Growth)

[系統理論] 離散時間系統的穩定度理論 (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

[動態系統] 非線性一階動態系統

考慮 1階  ODE 動態系統 可表示 如下
\[
\dot x = f(x)
\]其中 $x \in \mathbb{R}^d$ 為狀態變數, $ f:\mathbb{R}^d \to \mathbb{R}^d$，亦即我們可將其改寫為
\[\begin{array}{l}
\dot x = f\left( x \right)\\
 \Rightarrow \left\{ \begin{array}{l}
{{\dot x}_1} = {f_1}\left( {{x_1},...,{x_d}} \right)\\
{{\dot x}_2} = {f_2}\left( {{x_1},...,{x_d}} \right)\\
 \vdots \\
{{\dot x}_d} = {f_d}\left( {{x_1},...,{x_d}} \right)
\end{array} \right.
\end{array}\] 上述系統我們稱做 autonomous system。

Comments:
1. autonomous system 的 autonomous 源自希臘，意指與時間無關 (independent of time $t$ )
2. 如果 $\dot x = f(x, t)$ 此時系統稱作 non-autonomous system。但我們可以透過引入新的 狀態變數 將系統改寫回 autonomous system。


Example
考慮 non-autonomous system 如下
\[
\dot x = f(x,t)
\]其中 $x \in \mathbb{R}^d$ 且 $f: \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}^d$。則我們可令 $y:= t$ 且 $\dot y =1$ 改寫原式
\[\dot x = f(x,t) \Rightarrow \left\{ \begin{array}{l}
\dot x = f(x,y)\\
\dot y = 1
\end{array} \right.\]此時我們可發現若將 提高系統 1維度 (引入額外狀態變數 $y := t$ )作為代價，則可以將 non-autonomous system 改寫為 autonomous system。

接著我們想探討儘管低維度情況，動態系統仍然非常複雜。

Example: Nonlinear Lorenz System
考慮 $(x_1,x_2,x_3) \in \mathbb{R}^3$，
\[\left\{ \begin{array}{l}
{{\dot x}_1} = \sigma \left( {{x_2} - {x_1}} \right)\\
{{\dot x}_2} = r{x_1} - {x_2} - {x_3}{x_1}\\
{{\dot x}_3} = {x_1}{x_2} - \beta {x_3}
\end{array} \right.\]上述系統稱為 Lorenz System ，其中 $\sigma, r,\beta$ 為系統參數。一般而言常見的情況為 $\sigma = 10; r=28\; \beta = 4/3$。Lorenz 發現上述非線性系統方程的解存在 混沌(Chaotic)現象 (亦即對初始條件十分敏感)，並用此系統加以闡述天氣的不可預測性 (蝴蝶效應)：

"Does the Flap of a Butterfly's Wings in Brazil Set off a Tornado in Texas?" E. Lorenz, 1972

這說明了儘管在低維度情況，非線性系統的行為仍然可能非常複雜。那麼在複雜之中是否可以看出一些端倪? 我們需要新的想法來幫助我們：比如說系統初始行為可能很複雜，但是否可以僅僅觀察其穩態的情況?


平衡點 (equilibrium point)
現在我們回頭考慮一般的動態方程： $\dot x = f(x)$，若 $\bar x \in \mathbb{R}^d$ 為 $f$ 的零點 (zero)，亦即
\[
f(\bar x) =0
\]則 動態方程有常數解 $x(t) = \bar x$，我們稱此點 $\bar x$ 為 平衡解 或稱 平衡點 (equilibrium point) 或稱 穩態解 (steady-state solution) 或者 固定點 (fixed point)。

注意到 平衡點可以是 穩定的平衡點  (受小擾動後處在平衡點的系統平衡狀態被打破但可再度回到同一個平衡點) 亦可為 不穩定的平衡點 (受小擾動後處在平衡點的系統平衡狀態被打破且不再回到平衡點 或者 產生來回震盪)。


穩定平衡點判別
考慮一 autonomous 系統為 $\frac{dx}{dt} = f(x)$ 且 $x = \bar x$ 為此系統之一組平衡點。若 \[{\left. {\frac{{df(x)}}{{dx}}} \right|_{x = \bar x}} < 0\]則我們說 $\bar x$ 為穩定平衡點