微分方程

（Differential Equation）是描述一個動態系統的時域數學模型，我們藉此來間接研究系統的行為。微分方程的性質和解結構，與系統的響應特點有本質上的聯繫，這點需要在進行以傳遞函數為主的頻域模型講解前說明。

本文探究三個方面，ODE與LTI系統響應的關係，以及線性化的方法與合理性。後者是本文的重點內容。文中多次出現系統響應，系統輸出，解軌跡都是指對應微分方程的解。

微分方程有很多種，我們在控制理論中多見的是常微分方程（Ordinary Differential Equation）。我們可以通過機理分析或者系統辨識來得到ODEs，從而進行基於模型的控制設計。

ODE分為線性(Linear)和非線性(Nonlinear)的。在經典控制的範疇中，我們關注點主要在線性常係數ODE，也會設計到非常係數的ODE。前者描述的系統我們稱之為線性時不變系統，或者常稱為線性定常系統（Linear Time-Invariant System，LTI System）。下面所述的ODE均為線性常係數ODE（Constant-Coefficient Linear ODE）

既然不是所有系統都是線性的，或者說，現實中的系統幾乎都存在非線性。利用線性常係數ODEs來描述系統總是有一定誤差的。因此我們自然想知道使用線性常係數ODEs的合理性。

ODE與系統響應

ODEs描述了系統的運動規律，ODE的解對應的就是某個系統的響應。我們所熟知的牛頓第二運動定律（Newtons second law of motion）就是由一個微分方程來刻畫的

F=mddot{x}

這是一個十分簡單的二階ODE，描述力物體受力後的位移變化規律。我們可以直接通過積分來得到其解。不過從一般解微分方程的角度來看，此方程的解依舊由兩部分組成： 通解（general solution）以及特解（particular solution）。

首先解得 mddot{x}=0的通解：特徵方程 {{r}^{2}} ext{=}0 。0為重根，因此通解為： x=A+Bt

特解： 目測可得其中一個特解 x=frac{F}{2m}{{t}^{2}}

由此可知全解為 x=A+Bt ext{+}frac{F}{2m}{{t}^{2}}

給定初始條件 x(0)=A_{t},dot{x}(0)=B_{t} ，那麼

x(0)=A=A_{t}

dot{x}(0)=B ext{+}frac{F}{m}cdot 0=B_{t}

此微分方程的解為： x=A_{t}+B_{t}t ext{+}frac{F}{2m}{{t}^{2}} (1)

特別的，當零初始條件時，解也就是高中物理裏常用的公式 x=frac{F}{2m}{{t}^{2}}

那麼當初始狀態為0時，我們稱這樣的解是被該ODE描述的系統的零狀態響應（Zero State Response ），此時系統的響應變化完全是由輸入F所驅動的，所以也可以理解為純輸入響應。 零狀態響應只需要將零初始條件賦予系統全解形式即可得到，本質上還是通解與特解之和，這裡的例子比較特殊，通解部分正好為0。

現在在某一時刻，我們突然撤走F，系統該如何變化呢？ 此時描述系統的ODE的右端F變成了0，由非齊次ODE變成了齊次ODE

mddot{x}=0

我們知道此時齊次系統的全解，也就是原非齊次的通解。假設初始條件為

x(0)=A_{t},dot{x}(0)=B_{t}

那麼對應的解為

egin{align} & x(t)={{A}_{t}}+{{B}_{t}}t \ end{align}

我們稱這樣的解為該ODE描述的系統的零輸入響應（Zero Input Response）。從撤走外力的那一刻起，系統的響應完全由系統自身的狀態決定。試想如果此時系統的初始狀態也為0，那麼x就會從一開始就保持為0。

ODE解的疊加性

我們把 mddot{x}=0 在非零初始條件 x(0)=A_{t},dot{x}(0)=B_{t} 下的解 egin{align} & x(t)={{A}_{t}}+{{B}_{t}}t \ end{align}

和 mddot{x}=F 的零初始條件解 x=frac{F}{2m}{{t}^{2}} 相加

得到的解：

x=A_{t}+B_{t}t ext{+}frac{F}{2m}{{t}^{2}} (2)

對比(1)(2) 完全一致。從這個例子不足以推廣到所有線性ODE，但是還是不妨礙我們搬出這個結論：

非齊次線性ODE的解，由其零初始條件下的非齊次方程的解 和 非零初始條件下齊次方程的解組成。其對應的系統的全響應，則由這兩部分對應的零狀態響應與零輸入響應相加而得。

線性常係數ODE 全解 → 齊次方程非零初始條件解 + 非齊次方程零初始條件解

LTI系統Complete response → Zero State Response+ Zero Input Response

對於其他線性常係數ODE來說，這一疊加性也都是滿足的。系統的解與ODE之間的關係，我們在下一章的傳遞函數和穩定性中會用到，對於理解傳遞函數的一些概念尤其重要。

ODE線性化方法

我們會遇到一些非線性的ODE，比如加入了一些三角函數或高階多項式

ddot{x} ext{+cos}(x)dot{x}+x+{{x}^{3}}=F

雖然這仍然是一個二階ODE，但由於是非線性的，並不滿足疊加原理。對於這個方程描述的非線性系統，我們仍然可以通過線性化的方法來找到一個它的線性版本。

線性化的理論基礎是泰勒展開（Taylor expansion），任意函數在定義域內某個點可以寫成多項式的求和形式。於是我們只取泰勒展開的一階導數項和常數項，作為原函數的近似。比如取展開點為x=0，上面式子中的cos(x)和x3根據泰勒公式有：

cos (x)=cos (0)-sin (0)cdot x-frac{1}{2}cos (0)cdot {{x}^{2}}...

{{x}^{3}}=0+0cdot x+0cdot {{x}^{2}}+{{x}^{3}}

略去高次多項式，只取其線性部分，我們有如下近似

[cos (x)approx cos (0)-sin (0)cdot x=1]

{{x}^{3}}approx 0

這樣上述非線性系統就近似成了x=0附近的線性系統

ddot{x}+dot{x}+x=F

注意，我們這裡展開時沒有指定 dot{x} 的值，那是因為 dot{x} 自身沒有非線性項存在，我們只要線性化其係數。但是如果出現 {{dot{x}}^{2}} 類似關於 dot{x} 的非線性函數，則 dot{x} 展開值必須考慮。即展開點一般情況下必須是 x={{C}_{1}},dot{x}={{C}_{2}} ，它是ODE解空間中的一個點。我們這個例子中隱含了 dot{x} 的初值與近似結果並沒有關係。更高階的n階系統，展開點將是n維解空間中的一個向量，也可以認為是一個n維空間中的點。

ODE線性化的合理性

線性系統的這種線性化（Linearization）近似合不合理呢？ 從結果上來看是合理的，線性理論造就了一大批控制系統，連飛機都用到了線性化模型，實踐看來是ok的。那麼我們自然要問，這種近似是任何時候都成立的嗎？ 答案是，No。

數學上我們知道只要x取值在展開點附近，泰勒展開的線性部分對原函數的近似還是令人滿意的。從上節線性化的方法中可以看，展開點似乎是可以任意選擇的，理論上可以把一個非線性ODE在不同的展開點處近似成無數個線性ODE。但線性方程與原非線性方程之間的近似關係只有在系統的解軌跡（trajectory）不遠離展開點附近的情況下才成立。根據泰勒公式我們知道，當 [left| x(t)-{{x}_{0}} ight|] ，即系統某時刻的解與展開點差的模，為一個比較小的數時，它的高階多項式才能忽略，線性化的近似關係才能成立。否則，高次項的值遠大於線性部分，線性方程描述的系統與原非線性方程描述的系統已經相去甚遠。

這樣產生了兩個問題：

  1. 原非線性方程解空間中這麼多點，究竟要在哪個點進行線性化？
  2. 線性化後得到的線性方程，怎麼保證解就能在原來展開點的附近呢？

我們嘗試回答這兩個問題

一，哪個點進行線性化，取決於我們希望系統穩態時所處的工作點（set point）。

在regulation問題中我們希望系統的輸出能夠穩定到set point。一旦有了外來擾動，那麼控制器就會給出控制信號來繼續穩定輸出。控制器設計指標之一就是不讓輸出有大幅度變化，要又快，又穩，又準地把輸出再維持到設定值。既然regulation中我們輸出一旦穩定便不再大幅度變化，我們就可以把線性化的展開點放到這個set point上，這個穩態時的點，也叫做系統的穩定平衡點（stable equilibrium）。具體可以這麼做：

  1. 選擇一個期望的系統穩定時的set point
  2. 對系統的非線性ODE在該工作點處進行線性化。
  3. 按照該線性系統的模型設計控制器，控制非線性系統的輸出（output），或者也說響應（response）穩定到這個工作點。

在tracking問題中，我們希望的是系統output與reference input的誤差能夠穩定到減小，正常情況下，我們當然是要求誤差越小越好，那樣就能很好地跟蹤輸入信號。於是我們將線性化誤差的微分方程，如之前的例子，我們想讓輸出追蹤一個光滑的reference x^{*} ，令 e=x-{{x}^{*}}Rightarrow x=e+{{x}^{*}} ,帶入原方程得到

[egin{align} & ddot{x} ext{+cos}(x)dot{x}+x+{{x}^{3}}=F \ & Rightarrow \ & ddot{e}+{{{ddot{x}}}^{*}}+cos (e+{{x}^{*}})+e+{{x}^{*}}+{{(e+{{x}^{*}})}^{3}}=F \ end{align}]

同樣選擇讓誤差e穩定到你希望的set point，一般都是0，來線性化這個非線性ODE。接下來的工作就是設計控制器。

二，控制器的設計以及擾動對線性化模型的有效性起了很大影響。

我們注意到上小節中的第3步，按照線性化模型設計控制器時，我們就假定了線性模型和原系統是可以近似等價的，但等價是有條件的。反過來設計出來的控制器，如果無法滿足近似等價的條件，那這時候控制器就可能要失效了，或者說是因為線性化模型已經不能很好地描述原來的非線性系統了。

我們說，選擇了某個set point進行線性化後，線性化模型與原系統之間擁有良好近似關係的關鍵就在於控制器能否把系統的輸出都控制在set point附近。另一方面，控制器controller的一大任務就是抗幹擾。如果幹擾特別大，可以想像輸出很容易大幅度偏離set point，也會威脅線性化模型效果。

控制器如何保證受擾後的系統依舊能夠把解穩定在set point附近呢？ 我們先說一個不加控制作用的系統，它天然地也可以具有把解穩定到自身穩定平衡點的作用。

來看非線性系統的例子，不加任何外力F的系統，則變成了齊次方程：

ddot{x} ext{+cos}(x)dot{x}+x+{{x}^{3}}=0

這個非線性系統原本就存在一個穩定平衡點 x=0,dot{x}=0 ，我們可以在MATLAB中畫一個解空間的軌跡圖。我取了不同的初值，讓 dot{x}(0)=0 都為0， x(0) 從0變化到4，步長為0.1。就有了下圖。

Solution of Nolinear System with multiple inital conditons

「這蜜汁漩渦…」 好像很密很亂，但是我們可以看到有一些解的軌跡最後「轉轉轉」，落入了原點（origin）。而明顯，我只取了x1 為0到4，那外圍那些大圈圈是個什麼鬼？

於是我讓 x(0) 從3變化到4，步長取0.5，只有三條軌跡，這下總清晰了吧。

我們可以看到初值為(3,0)的解隨著時間變化有落入原點（0,0）的趨勢，好像被它吸引了，嗯，那就叫它吸引子（attractor）好了 ^{[Note1]} 。初值為(3.5,0)和(4,0)的值都沒有被origin吸引，反而都往「外面」發散了。如果把坐標軸延長，這種發散將延伸到解空間的無窮遠處。回過頭來看第一張圖，直覺上告訴我們，似乎初值 x(0) 值小於3， dot{x}(0) 為0時的方程解最終都能被吸引到原點，最後就停留在了原點。而大於3的那些解都跑得遠遠的，再也回不來了。

我們的判斷基本是正確的， ddot{x} ext{+cos}(x)dot{x}+x+{{x}^{3}}=0 的穩定平衡點就在原點origin，因為在穩定平衡點附近的解都被吸引到平衡點上去了。我們將來可以證明，origin在這裡就是stable equilibrium。

那麼這個「附近」的概念是多大呢？初值選取不當會導致解解的發散。我們從圖像中也能有個大概估計。這次我讓 dot{x} 也變化，選取了一些初值的區域，粗略地估計了這個非線性系統的吸引域（region of attraction）。

Estimation of region of attraction （Red part）

紅色部分就是吸引域的估計，只要初值落在紅色區域內 ^{[Note 2]} ，那麼方程的解最後就會落入原點origin。在外面的初值會導致解發散到無窮遠。

那麼我們能否改變系統的穩定平衡點，來滿足我們的需求呢？

這個例子讓我們意識到，設計控制器的目的其實就是改變系統結構，控制器的任務之一就是把系統的穩定平衡點給「挪」到我們自己設定的set point上。如果系統本來就沒有這樣的穩定平衡點，控制器的任務之一就是「創造」這樣的穩定平衡點 ^{[Note3]} 。加了控制器的系統能夠擁有上述非線性系統的特性，即把set point作為系統的穩定平衡點，把set point附近的解通通吸引過去。這樣就能做到解總是在set point附近，從而保證了按照線性系統設計的控制器對非線性系統能夠有效地控制。

線性控制系統的響應一般都會控制在小範圍的波動。一旦由於外界幹擾，或者參考輸入信號突然的大幅度改變，導致誤差驟增，就容易導致系統的不穩定，這十分考驗控制器的能力。這是因為一旦誤差驟增，與控制器自身的增益相乘後，將會得到一個數值很大的控制信號。這樣的信號是為了讓系統能夠迅速修正誤差，重新進入穩態。這點很容易理解，就好像你猛踩油門，為了達到120碼一樣。但這可能造成系統輸出在短時間內的大幅度的波動，從而讓解在某一時刻跑出被控系統的吸引域。一旦跑出吸引域，麻煩就大了，這時候解可能就再也跑不回set point了，而解持續增大，從而線性化模型就失去了近似的意義，由線性化模型設計的控制器就很有可能會失效。

實際系統的輸入突然變得很大，不管何種類型的輸入：控制器增益K太大導致控制器輸出的信號過大，還是外界突然來了一個強幹擾繞使得幹擾輸入過大，又或者一個參考信號突然丟失從而誤差變大，這都會導致plant或者process的輸入驟然變大，從而使得輸出跑出其穩定平衡點的吸引範圍。即便這些問題修復了，系統從那一時刻起的初始條件已然不會再讓系統的輸出返回穩定平衡點，也就是set point了。系統的非線性特性可能隨著這種輸入的突然增大而被激發，很有可能會讓線性控制器失效。對於更複雜的非線性系統，情況會更加複雜，它們可能會擁有多個平衡點，我們不希望解最後落入我們不想讓它去的地方。這些日後有機會再提。

間接地我們得到一個信息，控制器的一大設計目標，自然就是擴大吸引域的範圍了。換句話說，這也就是控制系統魯棒性（robustness）的表現，因為外界激勵的增大導致線性模型失效其實也等同於模型參數的改變。如果不能實時地更新線性模型的參數和控制器的參數，那麼這樣的控制器就不能適應較大模型變動。好點的情況是控制器依舊能夠穩定系統，但是控制性能指標得不到滿足，差點的情況就是直接無法穩定而導致系統fail。

不光光是線性控制器，即便是非線性控制器，也需要盡量擴大吸引域。很多時候很難通過設計來賦予系統一個全局穩定平衡點（globally stable equilibrium） ^{[Note4]} ，即無論從哪個初值出發都最後落入穩定平衡點。這些就涉及到非線性控制（nonlinear control）和非線性系統（nonlinear system）上去了，這裡不再多展開了。

總結一下

線性常係數ODE 的解描述了LTI系統的響應，他們之間有對應關係。

系統的線性化理論基礎來自泰勒展開和小誤差近似。

非線性系統的近似線性化，其初值條件一般不能偏離展開點太遠，否則容易產生大偏差，而導致線性化模型不再準確，導致線性控制器的控制效果變差。甚至會激發原系統的非線性特性，控制器失效，導致系統變得不穩定。較大的幹擾往往容易引發線性控制器的失效，此時引入非線性控制器可以較好地控制幹擾，幫助系統的工作區域能擴展到更大的範圍而保持穩定。

Note：

[1] 這個名字取得並不草率，而是相當專業的非線性動力學名詞。

[2] 嚴格地說不能隨便下這種結論，但是對於這裡例子平衡點確實只有一個，所以這些紅點大致覆蓋了吸引域。對於有多個平衡點的非線性系統，粗略地畫圖是不能判斷這片紅色區域裏還有沒有其他不穩定平衡點的。

[3] 這裡涉及到了一個是不是任何系統都能這樣挪的問題，要考慮能控性，這裡先不管。

[4] 線性系統本身只有一個平衡點，確實是全局穩定的平衡點，即吸引域無窮大。我們之後會再提。

Reference

[1] Hassan K.Khali，Nonlinear Systems , 3rd Edition, 2014, Pearson

[2] JJ E.Slotine, Weiping Li, Applied Nonlinear Control, 1991, Prentice Hall

[3] 胡壽松，自動控制原理（第六版），2013，科學出版社

Appendix

MATLAB codes for estimation of region of attraction(ROA)

% SJM Frank
% estimation of ROA
% 06-22-2019
dt = 0.01; %step size
t_end = 50;
t = 0:dt:t_end; %time sequence
N = size(t,2);

y =zeros(2,N);
for j = -10:0.2:10 % initial value for x1
for k = -10:0.2:10 % initial value for x2
x =zeros(2,N);
x(:,1)=[j;k];

for i = 1 : N-1 %Euler method
x(1,i+1) = x(1,i) + x(2,i) * dt;
x(2,i+1) = x(2,i)+ (-cos(x(1,i)).*x(2,i)-x(1,i)-x(1,i).^3) * dt;
end
hold on
if abs (x(1,end))<0.1 ||abs( x(2,end))<0.1 % if converges. 0.1 is changable.
plot(j,k,xr)
end
end
end
ylim([-10 10])
xlim([-10 10])
xlabel(x1)
ylabel(x2)
% legend(x(0)=3,x(0)=3.5,x(0)=4)
% figure;
% plot(y(1,:),y(2,:))

若有紕漏，煩請指出。轉載還請私信聯繫，請勿私自轉載。O(∩_∩)O


推薦閱讀：
相關文章