控制系統/線性系統解決方案

狀態方程解

本章的解決方案深深紮根於對 常微分方程 的先驗知識。讀者在閱讀本章之前應該先了解該主題。

狀態方程是一個一階線性微分方程，或更準確地說，是一個線性微分方程組。由於這是一個一階方程，我們可以使用 常微分方程 的結果來找到關於狀態變數x 的方程的通解。一旦狀態方程對x 求解，該解就可以代入輸出方程。得到的方程將顯示系統輸入和系統輸出之間的直接關係，無需明確考慮系統的內部狀態。本章中的各節將討論狀態空間方程的解，從最簡單的案例（時不變、無輸入）開始，到最複雜的案例（時變系統）結束。

求解無輸入時的 x(t)

再次看看狀態方程

$x'=Ax(t)+Bu(t)$

我們可以看到，這個方程是一個一階微分方程，只是變數是向量，係數是矩陣。然而，由於矩陣微積分的規則，這些區別並不重要。我們可以忽略輸入項（現在），並將該方程改寫為以下形式

${\frac {dx(t)}{dt}}=Ax(t)$

我們可以像這樣分離變數

${\frac {dx(t)}{x(t)}}=Adt$

對兩邊積分，並將兩邊都乘以e 的冪，我們得到結果

$x(t)=e^{At+C}$

其中C 是一個常數。我們可以將D = e^C 賦值以使方程更容易，但我們也知道，D 然後將成為系統的初始條件。如果我們將變數t 中的值零代入，這將變得很明顯。該方程的最終解如下所示

$x(t)=e^{A(t-t_{0})}x(t_{0})$

我們將矩陣指數e^At 稱為狀態轉移矩陣，計算它雖然有時很困難，但對於分析和操縱系統至關重要。我們將在下面詳細討論計算矩陣指數。

求解有非零輸入時的 x(t)

然而，如果我們的輸入是非零的（通常情況下，任何有趣的系統都是這樣），我們的解決方案會稍微複雜一些。請注意，現在我們有了輸入項，我們不能再輕鬆地分離變數並對兩邊積分。

$x'(t)=Ax(t)+Bu(t)$

我們減去以獲得 $Ax(t)$ 在左側，然後我們做一些奇怪的事情；我們用逆狀態轉移矩陣在兩邊預乘

$e^{-At}x'(t)-e^{-At}Ax(t)=e^{-At}Bu(t)$

這最後一步的理由可能看起來很模糊，所以我們將用一個例子來說明這一點

示例

對以下表達式關於時間求導

$e^{-At}x(t)$

微分中的乘積法則提醒我們，如果我們有兩個函式相乘

$f(t)g(t)$

然後我們對 t 求導，結果為

$f(t)g'(t)+f'(t)g(t)$

如果我們相應地設定我們的函式

$f(t)=e^{-At}\qquad f'(t)=-Ae^{-At}$

$g(t)=x(t)\qquad g'(t)=x'(t)$

則輸出結果為

$e^{-At}x'(t)-e^{-At}Ax(t)$

如果我們看一下這個結果，它與我們上面方程的結果相同。

使用我們例子中的結果，我們可以將我們方程左側濃縮成一個導數

${\frac {d(e^{-At}x(t))}{dt}}=e^{-At}Bu(t)$

現在我們可以對兩邊積分，從初始時間 (t₀) 到當前時間 (t)，使用一個啞變數 τ，我們將更接近我們的結果。最後，如果我們預乘以 e^At，我們將得到最終結果

[一般狀態方程解]

$x(t)=e^{A(t-t_{0})}x(t_{0})+\int _{t_{0}}^{t}e^{A(t-\tau )}Bu(\tau )d\tau$

如果我們將這個解代入輸出方程，我們將得到

[一般輸出方程解]

$y(t)=Ce^{A(t-t_{0})}x(t_{0})+C\int _{t_{0}}^{t}e^{A(t-\tau )}Bu(\tau )d\tau +Du(t)$

這是狀態空間方程的非零輸入的廣義時不變解。這些方程是重要的結果，對控制系統有進一步研究興趣的學生應該記住這些方程。

狀態轉移矩陣

有關矩陣指數的更多資訊，請參閱
工程分析

狀態轉移矩陣，e^At，是上述時不變情況下廣義狀態空間解的重要組成部分。計算這個矩陣指數函式是分析新系統時首先要做的幾件事之一，計算結果將揭示有關所討論系統的關鍵資訊。

可以使用泰勒級數展開直接計算矩陣指數

$e^{At}=\sum _{n=0}^{\infty }{\frac {(At)^{n}}{n!}}$

有關對角矩陣和約旦型矩陣的更多資訊，請參閱
工程分析

此外，我們可以嘗試將矩陣 A 對角化為對角矩陣或約旦標準型矩陣。對角矩陣的指數只是對角元素分別乘以該指數。約旦標準型矩陣的指數稍微複雜一些，但存在一個有用的模式可以利用來快速找到解。有興趣的讀者應該閱讀工程分析中的相關段落。

狀態轉移矩陣，以及一般意義上的矩陣指數，是控制工程中非常重要的工具。

對角矩陣

如果矩陣是對角矩陣，則狀態轉移矩陣可以透過將矩陣的每個對角元素作為e的冪次方來計算。

約旦標準型

如果A矩陣是約旦標準型，則可以使用以下公式快速生成矩陣指數

$e^{Jt}=e^{\lambda t}{\begin{bmatrix}1&t&{\frac {1}{2!}}t^{2}&\cdots &{\frac {1}{n!}}t^{n}\\0&1&t&\cdots &{\frac {1}{(n-1)!}}t^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\cdots &1\end{bmatrix}}$

其中λ是約旦標準型矩陣的特徵值（對角線上的值）。

拉普拉斯反變換方法

我們可以透過對以下拉普拉斯反變換求解來計算狀態轉移矩陣（或任何矩陣指數函式）

$e^{At}={\mathcal {L}}^{-1}[(sI-A)^{-1}]$

如果A是高階矩陣，則求逆運算可能難以解決。

如果A矩陣是約旦標準型，則可以使用以下公式快速生成矩陣指數

   $e^{Jt}=e^{\lambda t}{\begin{bmatrix}1&t&{\frac {1}{2!}}t^{2}&\cdots &{\frac {1}{n!}}t^{n}\\0&1&t&\cdots &{\frac {1}{(n-1)!}}t^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\cdots &1\end{bmatrix}}$

其中λ是約旦標準型矩陣的特徵值（對角線上的值）。

譜分解

如果我們知道A的所有特徵值，我們可以建立我們的轉移矩陣T和逆轉移矩陣T^-1。這些矩陣將分別是右特徵向量和左特徵向量的矩陣。如果我們同時擁有左特徵向量和右特徵向量，我們可以計算出狀態轉移矩陣為

[譜分解]

$e^{At}=\sum _{i=1}^{n}e^{\lambda _{i}t}v_{i}w_{i}'$

請注意，w_i' 是第i個左特徵向量的轉置，而不是它的導數。我們將在後面的章節中詳細討論“特徵值”、“特徵向量”以及譜分解技術。

凱萊-哈密頓定理

有關凱萊-哈密頓定理的更多資訊，請參閱
工程分析

凱萊-哈密頓定理也可以用來找到矩陣指數的解。對於系統矩陣A的任何特徵值λ，我們可以證明這兩個方程是等價的

$e^{\lambda t}=a_{0}+a_{1}\lambda t+a_{2}\lambda ^{2}t^{2}+\cdots +a_{n-1}\lambda ^{n-1}t^{n-1}$

一旦我們解出方程的係數a，就可以將這些係數代入以下方程

$e^{At}=a_{0}I+a_{1}At+a_{2}A^{2}t^{2}+\cdots +a_{n-1}A^{n-1}t^{n-1}$

示例：非對角矩陣

給定以下矩陣A，求狀態轉移矩陣

$A={\begin{bmatrix}0&1\\-1&0\end{bmatrix}}$

我們可以找到該矩陣的特徵值為λ = i，-i。如果我們將這些值代入我們的特徵向量方程，得到

${\begin{vmatrix}i&-1\\1&i\end{vmatrix}}v_{1}=0$

${\begin{vmatrix}-i&-1\\1&-i\end{vmatrix}}v_{2}=0$

我們可以解出特徵向量

$v_{1}={\begin{bmatrix}1\\i\end{bmatrix}}$

$v_{2}={\begin{bmatrix}1\\-i\end{bmatrix}}$

有了特徵向量，我們就可以解出左特徵向量

$w_{1}={\begin{bmatrix}1\\-i\end{bmatrix}}$

$w_{2}={\begin{bmatrix}1\\i\end{bmatrix}}$

現在，利用譜分解，我們可以構建狀態轉移矩陣

$e^{At}=e^{it}{\begin{bmatrix}1\\i\end{bmatrix}}{\begin{bmatrix}1&-i\end{bmatrix}}+e^{-it}{\begin{bmatrix}1\\-i\end{bmatrix}}{\begin{bmatrix}1&i\end{bmatrix}}$

如果我們記得尤拉恆等式，我們可以將復指數分解成正弦函式。執行向量相乘後，所有的虛數項都會抵消，剩下的就是我們的結果

$e^{At}={\begin{bmatrix}\cos t&\sin t\\-\sin t&\cos t\end{bmatrix}}$

鼓勵讀者進行乘法運算，並嘗試推匯出該結果。

示例：Sympy 計算

使用免費的 Python 庫 'sympy'，我們可以非常容易地自動計算狀態轉移矩陣

>>> from sympy import *
>>> t = symbols('t', positive = true)
>>> A = Matrix([[0,1],[-1,0]])
>>> exp(A*t).expand(complex=True)

⎡cos(t)   sin(t)⎤
⎢               ⎥
⎣-sin(t)  cos(t)⎦

你也可以在這個網站上嘗試一下

sympy live

示例：MATLAB 計算

使用 MATLAB 中的符號工具箱，我們可以編寫 MATLAB 程式碼來自動生成給定輸入矩陣A的狀態轉移矩陣。以下是一個可以執行此任務的 MATLAB 程式碼示例

function [phi] = statetrans(A)
   t = sym('t');
   phi = expm(A * t);
end

使用此 MATLAB 函式來查詢以下矩陣的狀態轉移矩陣（警告，計算可能需要一些時間）

1. $A_{1}={\begin{bmatrix}2&0\\0&2\end{bmatrix}}$
2. $A_{2}={\begin{bmatrix}0&1\\-1&0\end{bmatrix}}$
3. $A_{3}={\begin{bmatrix}2&1\\0&2\end{bmatrix}}$

矩陣 1 是一個對角矩陣，矩陣 2 具有復特徵值，矩陣 3 是約旦標準型。這三個矩陣應該代表系統矩陣的一些常見形式。以下程式碼片段是輸入到 MATLAB 中以生成這些矩陣的輸入命令以及輸出結果

矩陣 A1

>> A1 = [2 0 ; 0 2];
>> statetrans(A1)
 
ans =
 
[ exp(2*t),        0]
[        0, exp(2*t)]

矩陣 A2

>> A2 = [0 1 ; -1 0];
>> statetrans(A1)
 
ans =
 
[  cos(t),  sin(t)]
[ -sin(t),  cos(t)]

矩陣 A3

>> A1 = [2 1 ; 0 2];
>> statetrans(A1)
 
ans =
 
[   exp(2*t), t*exp(2*t)]
[          0,   exp(2*t)]

示例：MATLAB 中的多種方法

在 MATLAB 中，有多種方法可以從標量（時不變）矩陣 A 計算狀態轉移矩陣。以下方法都將依賴於 *符號工具箱* 來執行方程操作。在每個程式碼片段的末尾，變數 **eAt** 包含矩陣 A 的狀態轉移矩陣。

直接法

t = sym('t');
eAt = expm(A * t);

拉普拉斯變換法

s = sym('s');
[n,n] = size(A);
in = inv(s*eye(n) - A);
eAt = ilaplace(in);

譜分解

t = sym('t');
[n,n] = size(A);
[V, e] = eig(A);
W = inv(V);
sum = [0 0;0 0];
for I = 1:n
   sum = sum + expm(e(I,I)*t)*V(:,I)*W(I,:);
end;
eAt = sum;

這三種方法都應該產生相同的答案。鼓勵學生驗證這一點。

← 狀態空間方程

控制系統

時變系統解 →