矩陣的核心原理、計算與MATLAB實踐)
1. 從“單輸入單輸出”到“多輸入多輸出”為什么我們需要傳遞函數(shù)矩陣如果你是從經(jīng)典控制理論比如拉普拉斯變換、伯德圖、奈奎斯特判據(jù)一路學(xué)過來的那么“傳遞函數(shù)”這個概念對你來說應(yīng)該像呼吸一樣自然。一個輸入一個輸出一個函數(shù)關(guān)系清晰明了。我們分析系統(tǒng)的穩(wěn)定性、動態(tài)性能都圍繞著這個單一的G(s)展開。但現(xiàn)實世界的系統(tǒng)尤其是現(xiàn)代工程系統(tǒng)很少有這么“單純”的。想想看一臺無人機你需要同時控制它的俯仰、橫滾、偏航和高度。你的操縱指令四個電機的轉(zhuǎn)速是多個輸入它的姿態(tài)和位置是多個輸出。一個化工反應(yīng)釜你需要控制溫度、壓力、pH值和進料流量。加熱功率、冷卻水閥、酸堿泵和進料泵是你的輸入傳感器讀數(shù)是你的輸出。一臺機械臂每個關(guān)節(jié)的電機扭矩是輸入末端執(zhí)行器的位置和姿態(tài)是輸出。這些系統(tǒng)都有一個共同點多輸入多輸出英文簡稱MIMO。在經(jīng)典控制里我們可能會嘗試為每一個輸出單獨設(shè)計一個控制器去調(diào)節(jié)某一個輸入這叫做“單回路控制”。但問題很快就來了——這些回路之間是耦合的。你調(diào)節(jié)無人機第一個電機想讓它抬頭結(jié)果它可能不僅抬頭還開始往左偏航。你加大反應(yīng)釜的加熱功率想升溫壓力可能也跟著飆升。這種耦合意味著你不能把系統(tǒng)簡單地拆成幾個獨立的單變量系統(tǒng)來處理。它們是一個整體牽一發(fā)而動全身。這時沿用單個傳遞函數(shù)的思路就捉襟見肘了。我們需要一個數(shù)學(xué)工具能夠同時描述所有輸入對所有輸出的影響關(guān)系。這個工具就是傳遞函數(shù)矩陣。你可以把它想象成一張關(guān)系網(wǎng)或者一個表格。假設(shè)系統(tǒng)有p個輸入q個輸出那么這個傳遞函數(shù)矩陣G(s)就是一個q x p的矩陣。矩陣中的每一個元素G_ij(s)都是一個傳遞函數(shù)它單獨地描述了第j個輸入對第i個輸出的影響而此時其他所有輸入都為零這是線性系統(tǒng)疊加原理的前提。所以當我們從“現(xiàn)控理論”的視角重新審視系統(tǒng)時傳遞函數(shù)矩陣就是我們連接系統(tǒng)外部描述輸入輸出關(guān)系和內(nèi)部狀態(tài)描述狀態(tài)空間方程的一座關(guān)鍵橋梁也是處理MIMO系統(tǒng)分析與設(shè)計問題的起點。它保留了傳遞函數(shù)直觀的頻率域特性又將系統(tǒng)的多變量特性用矩陣這一強大的數(shù)學(xué)語言清晰地表達了出來。2. 傳遞函數(shù)矩陣的核心定義與數(shù)學(xué)表達理解了為什么需要它我們來看看它具體是什么。我們從最通用的線性時不變系統(tǒng)的狀態(tài)空間描述出發(fā)這是現(xiàn)代控制理論的基礎(chǔ)模型狀態(tài)方程 ?(t) A x(t) B u(t) 輸出方程 y(t) C x(t) D u(t)其中x(t)是n維狀態(tài)向量。u(t)是p維輸入向量。y(t)是q維輸出向量。A, B, C, D是維數(shù)匹配的常數(shù)矩陣。我們對上述方程兩邊進行拉普拉斯變換并假設(shè)初始狀態(tài)x(0) 0這是為了聚焦于輸入輸出的傳遞特性就像經(jīng)典控制里一樣。得到sX(s) A X(s) B U(s) Y(s) C X(s) D U(s)由第一個式子解出X(s)(sI - A) X(s) B U(s) X(s) (sI - A)^{-1} B U(s)這里I是n x n的單位矩陣(sI - A)^{-1}就是預(yù)解矩陣它在系統(tǒng)分析中至關(guān)重要。將X(s)代入輸出方程Y(s) C (sI - A)^{-1} B U(s) D U(s) [C (sI - A)^{-1} B D] U(s)于是我們得到了輸入U(s)和輸出Y(s)在復(fù)頻域的關(guān)系Y(s) G(s) U(s)其中傳遞函數(shù)矩陣G(s)被定義為G(s) C (sI - A)^{-1} B D這是一個q x p的矩陣。它的每一個元素G_ij(s)都是一個關(guān)于復(fù)變量s的有理分式函數(shù)分子分母都是s的多項式。G_ij(s)表示的是當只有第j個輸入U_j(s)作用時所產(chǎn)生第i個輸出Y_i(s)的傳遞函數(shù)。一個關(guān)鍵的理解點(sI - A)^{-1}包含了系統(tǒng)的所有動態(tài)模態(tài)信息極點而矩陣C和B則決定了這些模態(tài)如何被輸入激發(fā)以及如何被輸出觀測。矩陣D代表了輸入到輸出的直接饋通在很多物理系統(tǒng)中如電路、機械D常常是零矩陣因為輸入通常不會瞬間無延遲地影響輸出。注意這里假設(shè)了(sI - A)是可逆的這要求s不是系統(tǒng)矩陣A的特征值。系統(tǒng)矩陣A的特征值正是傳遞函數(shù)矩陣G(s)的極點分母為零的點它們決定了系統(tǒng)的穩(wěn)定性。3. 如何計算與解讀傳遞函數(shù)矩陣一個詳細的計算實例定義看起來有點抽象我們通過一個具體的、簡化的雙輸入雙輸出系統(tǒng)來親手算一遍感受一下這個過程?？紤]一個耦合的彈簧-質(zhì)量塊系統(tǒng)或者一個簡單的電路網(wǎng)絡(luò)我們可以抽象出如下狀態(tài)空間模型為了計算方便數(shù)字是假設(shè)的設(shè)系統(tǒng)方程為A [ -2 1 0; 1 -3 1; 0 1 -1 ]; B [ 1 0; 0 1; 1 0 ]; C [ 1 0 0; 0 1 0 ]; D [ 0 0; 0 0 ];這里狀態(tài)x是3維輸入u是2維輸出y是2維。D矩陣為零表示沒有直接饋通。我們的目標是求出G(s) C (sI - A)^{-1} B。第一步構(gòu)造 (sI - A)sI - A [ s2 -1 0; -1 s3 -1; 0 -1 s1 ];第二步求 (sI - A) 的逆矩陣求一個3x3矩陣的逆我們可以用伴隨矩陣法(sI - A)^{-1} adj(sI - A) / det(sI - A)。 先求行列式Δ(s) det(sI - A)Δ(s) (s2)*[(s3)(s1) - (-1)(-1)] - (-1)*[(-1)(s1) - (0)(-1)] 0*... (s2)[(s3)(s1) - 1] 1*[-(s1)] (s2)(s^2 4s 3 - 1) - (s1) (s2)(s^2 4s 2) - s - 1 s^3 4s^2 2s 2s^2 8s 4 - s - 1 s^3 6s^2 9s 3所以系統(tǒng)的特征多項式是s^3 6s^2 9s 3它的根就是系統(tǒng)的極點。再求伴隨矩陣adj(sI - A)。這需要計算所有代數(shù)余子式過程略這是基本功練習(xí)。假設(shè)我們算得adj(sI - A) [ (s3)(s1)-1 (s1) 1; (s1) (s2)(s1) (s2); 1 (s2) (s2)(s3)-1 ] [ s^24s2 s1 1; s1 s^23s2 s2; 1 s2 s^25s5 ]因此(sI - A)^{-1} (1 / Δ(s)) * [ s^24s2, s1, 1; s1, s^23s2, s2; 1, s2, s^25s5 ]第三步計算 C (sI - A)^{-1} B這是一個連續(xù)的矩陣乘法。C是 2x3逆矩陣是 3x3B是 3x2。所以結(jié)果G(s)是 2x2。 我們先計算中間結(jié)果T(s) (sI - A)^{-1} B。B矩陣有兩列我們分別計算。 令B1 [1; 0; 1](B的第一列)B2 [0; 1; 0](B的第二列)。計算T1(s) (sI - A)^{-1} * B1T1(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [1; (1/Δ(s)) * [ (s^24s2)*1 (s1)*0 1*1; s1, s^23s2, s2; 0; (s1)*1 (s^23s2)*0 (s2)*1; 1, s2, s^25s5 ] 1] 1*1 (s2)*0 (s^25s5)*1 ] (1/Δ(s)) * [ s^24s2 1; s1 s2; 1 s^25s5 ] (1/Δ(s)) * [ s^24s3; 2s3; s^25s6 ]同理計算T2(s) (sI - A)^{-1} * B2T2(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [0; (1/Δ(s)) * [ (s^24s2)*0 (s1)*1 1*0; s1, s^23s2, s2; 1; (s1)*0 (s^23s2)*1 (s2)*0; 1, s2, s^25s5 ] 0] 1*0 (s2)*1 (s^25s5)*0 ] (1/Δ(s)) * [ s1; s^23s2; s2 ]所以T(s) [T1(s), T2(s)] (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2; s^25s6, s2 ]現(xiàn)在左乘C矩陣。C矩陣只取前兩行因為輸出y只對應(yīng)前兩個狀態(tài)所以G(s) C * T(s)相當于取T(s)的前兩行。G(s) (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2 ]其中Δ(s) s^3 6s^2 9s 3。解讀G(s) 這個 2x2 的傳遞函數(shù)矩陣G(s)告訴我們G_11(s) (s^24s3) / (s^36s^29s3)輸入u1對輸出y1的傳遞函數(shù)。G_12(s) (s1) / (s^36s^29s3)輸入u2對輸出y1的傳遞函數(shù)。G_21(s) (2s3) / (s^36s^29s3)輸入u1對輸出y2的傳遞函數(shù)。G_22(s) (s^23s2) / (s^36s^29s3)輸入u2對輸出y2的傳遞函數(shù)。關(guān)鍵觀察共同分母所有四個傳遞函數(shù)都有相同的分母多項式Δ(s)。這個多項式正是矩陣A的特征多項式。這意味著整個MIMO系統(tǒng)的極點決定穩(wěn)定性和基本動態(tài)是由矩陣A決定的是所有輸入輸出通道所共享的。這是MIMO系統(tǒng)與多個獨立SISO系統(tǒng)本質(zhì)不同的地方。耦合體現(xiàn)非對角線元素G_12(s)和G_21(s)不為零。這說明輸入u1會影響輸出y2輸入u2也會影響輸出y1。系統(tǒng)是耦合的。相對階G_11和G_22的分子是二階分母是三階相對階為1。G_12和G_21的分子是一階相對階為2。不同通道的動態(tài)特性可能不同。實操心得對于維數(shù)高于3的系統(tǒng)手工計算逆矩陣會非常繁瑣且容易出錯。在實際工程或?qū)W習(xí)中我們強烈依賴數(shù)學(xué)工具。在MATLAB/Octave中計算傳遞函數(shù)矩陣極其簡單。假設(shè)你已經(jīng)定義了A, B, C, D矩陣只需一行命令G ss(A, B, C, D); G_tf tf(G);。ss創(chuàng)建狀態(tài)空間模型tf將其轉(zhuǎn)換為傳遞函數(shù)矩陣形式。對于我們的例子在MATLAB中驗證上述計算是很好的練習(xí)。手工推導(dǎo)的目的在于深刻理解其數(shù)學(xué)本質(zhì)和耦合關(guān)系的來源而不是用于解決大規(guī)模問題。4. 傳遞函數(shù)矩陣的性質(zhì)與在系統(tǒng)分析中的核心作用得到了傳遞函數(shù)矩陣我們能用它做什么它不僅僅是狀態(tài)空間模型的一種等價表示更是我們分析MIMO系統(tǒng)特性的強大工具。4.1 極點與零點MIMO系統(tǒng)的擴展定義對于SISO系統(tǒng)極點就是傳遞函數(shù)分母多項式的根零點就是分子多項式的根。對于MIMO系統(tǒng)由于G(s)是一個矩陣極點和零點的定義需要擴展。極點如前所述傳遞函數(shù)矩陣G(s)的所有元素的分母多項式在約去公因子后的根就是系統(tǒng)的極點。更本質(zhì)地說系統(tǒng)的極點就是系統(tǒng)矩陣A的特征值。這一定義從狀態(tài)空間模型出發(fā)更為根本。所有輸入輸出通道共享同一組極點這決定了系統(tǒng)的穩(wěn)定性所有極點均具有負實部則漸近穩(wěn)定和基本響應(yīng)模態(tài)如振蕩頻率、衰減速度。傳輸零點這是一個MIMO系統(tǒng)特有的、非常重要的概念。它不是單個元素分子為零的點。傳輸零點的定義是存在一個非零的復(fù)頻率z和一個非零的輸入向量U_0使得在零初始狀態(tài)下系統(tǒng)的輸出Y(s)恒為零。即G(z) U_0 0。這意味著在頻率z處存在某種特定的輸入信號組合其效果被系統(tǒng)內(nèi)部完全“抵消”了無法在輸出端被觀測到。計算對于D0的系統(tǒng)傳輸零點z是使得復(fù)合矩陣P(s) [sI-A, -B; C, 0]降秩的s值。在MATLAB中可以使用tzero或zero函數(shù)直接計算。物理意義傳輸零點反映了系統(tǒng)輸入輸出之間的阻塞特性。例如在飛機控制中某些特定的舵面偏轉(zhuǎn)組合可能無法引起飛機姿態(tài)的變化在某個頻率下這個頻率就是傳輸零點。傳輸零點會影響系統(tǒng)的可控制性、可觀測性以及控制性能的極限例如非最小相位系統(tǒng)的右半平面零點會限制控制帶寬。4.2 穩(wěn)定性分析基于極點的穩(wěn)定性判據(jù)在MIMO系統(tǒng)中依然直接適用線性時不變系統(tǒng)漸近穩(wěn)定的充要條件是其傳遞函數(shù)矩陣的所有極點即矩陣A的所有特征值都具有負實部。通過傳遞函數(shù)矩陣G(s)我們可以計算其特征多項式即行列式det(sI-A)或G(s)各元素分母的最小公倍式然后應(yīng)用勞斯判據(jù)、赫爾維茨判據(jù)或直接求根來判斷穩(wěn)定性。在MATLAB中使用pole(G)或eig(A)直接獲取極點。注意對于MIMO系統(tǒng)不能僅僅通過觀察G(s)某個對角元比如G_11(s)的穩(wěn)定性來判斷整個系統(tǒng)的穩(wěn)定性。因為非對角元的耦合可能引入不穩(wěn)定的隱藏模態(tài)與狀態(tài)空間中的能控性/能觀性相關(guān)。必須檢查系統(tǒng)矩陣A的全部特征值。4.3 頻域分析從伯德圖到奇異值圖這是傳遞函數(shù)矩陣威力巨大的地方。在SISO系統(tǒng)中我們繪制單個傳遞函數(shù)G(jω)的伯德圖幅頻和相頻特性。在MIMO系統(tǒng)中G(jω)是一個復(fù)數(shù)矩陣。我們?nèi)绾卫L制它的頻率響應(yīng)逐個元素法可以為G(s)的每一個元素G_ij(jω)繪制伯德圖。這能讓我們看清每一對輸入輸出通道在頻域的特性。這對于理解耦合的強度比如G_21相對于G_22的幅值大小很有幫助。在MATLAB中bode(G)命令會自動為所有通道生成伯德圖。奇異值分析法更強大的工具對于MIMO系統(tǒng)更本質(zhì)的頻域分析工具是奇異值。對于每個頻率點ω我們計算復(fù)數(shù)矩陣G(jω)的奇異值。奇異值總是非負的實數(shù)。假設(shè)G是q x p矩陣那么它在每個頻率點有min(p, q)個奇異值記為σ_1(ω) ≥ σ_2(ω) ≥ ... ≥ σ_min(p,q)(ω) ≥ 0。最大奇異值σ_max(ω)可以理解為系統(tǒng)在頻率ω處對所有可能輸入方向向量的最大增益。最小奇異值σ_min(ω)可以理解為系統(tǒng)在頻率ω處對所有可能輸入方向的最小增益。繪制σ_max(ω)和σ_min(ω)隨頻率變化的曲線就得到了MIMO系統(tǒng)的奇異值伯德圖。為什么重要在魯棒控制和系統(tǒng)性能分析中奇異值提供了關(guān)鍵的洞察。例如σ_min在低頻段的大小反映了系統(tǒng)抗干擾和解耦的能力σ_max在高頻段的衰減速率反映了系統(tǒng)對噪聲和未建模動態(tài)的魯棒性。閉環(huán)系統(tǒng)的穩(wěn)定裕度也可以用開環(huán)傳遞函數(shù)矩陣的奇異值來評估。在MATLAB中可以使用sigma(G)命令繪制奇異值圖。4.4 能控性與能觀性與傳遞函數(shù)矩陣的關(guān)聯(lián)能控性和能觀性是狀態(tài)空間模型的核心概念它們與傳遞函數(shù)矩陣有著深刻的聯(lián)系。能控性系統(tǒng)是否能在有限時間內(nèi)通過合適的輸入u(t)將狀態(tài)從任意初始點驅(qū)動到原點。狀態(tài)空間判據(jù)是能控性矩陣[B, AB, A^2B, ..., A^{n-1}B]滿秩。能觀性系統(tǒng)是否能在有限時間內(nèi)通過輸出的觀測y(t)唯一地確定初始狀態(tài)x(0)。狀態(tài)空間判據(jù)是能觀性矩陣[C; CA; CA^2; ...; CA^{n-1}]滿秩。從傳遞函數(shù)矩陣的角度看如果系統(tǒng)是狀態(tài)空間能控且能觀的那么其傳遞函數(shù)矩陣G(s)將是既約的即沒有零極點對消。如果發(fā)生了零極點對消則意味著系統(tǒng)不是完全能控和/或完全能觀的對消掉的模態(tài)在輸入輸出描述中“消失”了但它們可能隱藏在系統(tǒng)內(nèi)部影響實際動態(tài)甚至穩(wěn)定性。在計算G(s) C(sI-A)^{-1}B D時如果最終得到的各元素傳遞函數(shù)有公因子被約去這些被約去的因子就對應(yīng)著不能控或不能觀的模態(tài)。因此在由傳遞函數(shù)矩陣反推或?qū)崿F(xiàn)狀態(tài)空間模型時必須非常小心要確保得到的是最小實現(xiàn)即完全能控且完全能觀。5. 從理論到實踐在MATLAB/SciPy中操作傳遞函數(shù)矩陣理論學(xué)習(xí)之后我們必須能在工具中實現(xiàn)它。這里以MATLAB及其開源替代品GNU Octave和Python的SciPy庫為例展示基本操作。5.1 MATLAB/Octave 環(huán)境1. 定義系統(tǒng)與計算傳遞函數(shù)矩陣% 定義狀態(tài)空間矩陣 (沿用之前的例子) A [-2, 1, 0; 1, -3, 1; 0, 1, -1]; B [1, 0; 0, 1; 1, 0]; C [1, 0, 0; 0, 1, 0]; D zeros(2,2); % 創(chuàng)建狀態(tài)空間模型對象 sys_ss ss(A, B, C, D); % 轉(zhuǎn)換為傳遞函數(shù)矩陣形式 sys_tf tf(sys_ss); disp(傳遞函數(shù)矩陣 G(s):); disp(sys_tf); % 輸出會顯示一個2x2的tf數(shù)組每個元素都是tf對象。 % 例如從輸入 1 到輸出 1: s^2 4 s 3 / s^3 6 s^2 9 s 3 % 從輸入 2 到輸出 1: s 1 / s^3 6 s^2 9 s 3 ... 等等2. 提取特定通道的傳遞函數(shù)% 提取輸入1到輸出1的SISO傳遞函數(shù) G11 sys_tf(1, 1); % 提取輸入2到輸出2的SISO傳遞函數(shù) G22 sys_tf(2, 2);3. 計算極點和零點% 計算系統(tǒng)極點 (即A的特征值) poles pole(sys_ss); % 或 pole(sys_tf) disp(系統(tǒng)極點); disp(poles); % 計算系統(tǒng)傳輸零點 zeros tzero(sys_ss); % 或 zero(sys_tf) disp(系統(tǒng)傳輸零點); disp(zeros);4. 頻域分析繪圖% 繪制所有通道的伯德圖 (4個子圖) figure; bode(sys_ss); % 或 bode(sys_tf) grid on; title(各通道伯德圖); % 繪制奇異值圖 (MIMO系統(tǒng)頻域分析利器) figure; sigma(sys_ss); % 或 sigma(sys_tf) grid on; title(系統(tǒng)奇異值圖);5. 時域仿真% 定義時間向量 t 0:0.01:10; % 定義輸入信號u1是階躍u2是正弦波 u [ones(size(t)); sin(2*t)]; % 注意維度是 2 x length(t) % 進行仿真 [y, t_out, x] lsim(sys_ss, u, t); % 繪制輸出響應(yīng) figure; subplot(2,1,1); plot(t_out, y(:,1)); ylabel(y1); grid on; legend(輸出1); subplot(2,1,2); plot(t_out, y(:,2)); ylabel(y2); xlabel(時間 (s)); grid on; legend(輸出2);5.2 Python (SciPy Matplotlib) 環(huán)境Python在科學(xué)計算領(lǐng)域應(yīng)用廣泛控制庫雖不如MATLAB專業(yè)但基礎(chǔ)功能完備。import numpy as np import matplotlib.pyplot as plt from scipy import signal # 1. 定義狀態(tài)空間矩陣 A np.array([[-2, 1, 0], [1, -3, 1], [0, 1, -1]]) B np.array([[1, 0], [0, 1], [1, 0]]) C np.array([[1, 0, 0], [0, 1, 0]]) D np.array([[0, 0], [0, 0]]) # 2. 創(chuàng)建狀態(tài)空間系統(tǒng) sys_ss signal.StateSpace(A, B, C, D) # 3. 轉(zhuǎn)換為傳遞函數(shù)形式 (注意scipy.signal 的 tf 表示是針對SISO的) # 對于MIMO我們需要逐個元素計算或使用專門的庫如 control需安裝pip install control # 這里演示使用 control 庫 (如果已安裝) try: import control as ct # 創(chuàng)建控制系統(tǒng)庫的狀態(tài)空間對象 sys_ct ct.ss(A, B, C, D) # 計算傳遞函數(shù)矩陣 (返回一個 TransferFunction 對象列表的列表) sys_tf ct.tf(sys_ct) print(傳遞函數(shù)矩陣:) print(sys_tf) except ImportError: print(control 庫未安裝。對于MIMO傳遞函數(shù)建議安裝 control 庫。) # 手動計算其中一個通道作為示例 (例如 G11) # G(s) C(sI-A)^-1 B D, 我們計算第一個元素 # 這里省略手動計算逆矩陣的代碼較為復(fù)雜。 # 作為替代我們可以直接進行時域/頻域仿真。 # 4. 計算極點和零點 (使用 scipy.signal 或 control) # 極點就是 A 的特征值 poles np.linalg.eigvals(A) print(\n系統(tǒng)極點 (A的特征值):) print(poles) # 傳輸零點計算較復(fù)雜control庫有提供 try: import control as ct zeros ct.zero(sys_ct) print(\n系統(tǒng)傳輸零點 (來自 control 庫):) print(zeros) except: print(\n傳輸零點計算需要 control 庫。) # 5. 時域仿真 (使用 scipy.signal) t np.arange(0, 10, 0.01) # 輸入u1 階躍 u2 正弦波 u np.vstack([np.ones_like(t), np.sin(2*t)]) # 形狀 (2, len(t)) # lsim 需要輸入 u 的形狀為 (len(t), num_inputs) t_out, y, x signal.lsim(sys_ss, Uu.T, Tt) # 注意 u.T 轉(zhuǎn)置 # 繪制結(jié)果 fig, (ax1, ax2) plt.subplots(2, 1, figsize(10, 6)) ax1.plot(t_out, y[:, 0]) ax1.set_ylabel(y1) ax1.grid(True) ax1.legend([輸出1]) ax2.plot(t_out, y[:, 1]) ax2.set_ylabel(y2) ax2.set_xlabel(時間 [s]) ax2.grid(True) ax2.legend([輸出2]) plt.suptitle(系統(tǒng)時域響應(yīng)) plt.tight_layout() plt.show() # 6. 頻域分析 - 伯德圖 (scipy.signal 的 bode 是 SISO 的) # 我們可以循環(huán)繪制每個通道 try: import control as ct plt.figure() ct.bode_plot(sys_ct) # control庫可以繪制MIMO伯德圖 plt.suptitle(各通道伯德圖 (control庫)) plt.tight_layout() plt.show() except: print(繪制MIMO伯德圖需要 control 庫。)實操心得與避坑指南工具選擇對于嚴肅的控制系統(tǒng)分析與設(shè)計MATLAB及其控制系統(tǒng)工具箱仍然是行業(yè)標桿文檔齊全、函數(shù)豐富。Python的control庫功能日益完善對于學(xué)習(xí)、研究和輕量級應(yīng)用是完全足夠的且免費開源。SciPy.signal更側(cè)重于信號處理對MIMO支持較弱。維度匹配在定義矩陣A, B, C, D時務(wù)必反復(fù)檢查維度。(n x n), (n x p), (q x n), (q x p)。一個常見的錯誤是B或C矩陣的維度定義反了。仿真輸入使用lsim或signal.lsim進行仿真時輸入矩陣U的維度容易出錯。在MATLAB中U的列數(shù)等于輸入數(shù)行數(shù)等于時間點數(shù)量。在SciPy的signal.lsim中U的形狀是(len(t), num_inputs)。務(wù)必查閱文檔或打印矩陣形狀來確認。最小實現(xiàn)如果你是從物理方程推導(dǎo)或系統(tǒng)辨識得到的傳遞函數(shù)矩陣在轉(zhuǎn)換為狀態(tài)空間模型時使用minreal函數(shù)MATLAB或minreal方法control庫來獲取最小實現(xiàn)消除不能控或不能觀的模態(tài)避免后續(xù)分析與設(shè)計出現(xiàn)問題。6. 傳遞函數(shù)矩陣的局限性與狀態(tài)空間法的優(yōu)勢盡管傳遞函數(shù)矩陣是分析MIMO系統(tǒng)的有力工具但它并非萬能也存在固有的局限性。理解這些局限性正是我們深入學(xué)習(xí)現(xiàn)代控制理論狀態(tài)空間法的動力。1. 僅適用于線性時不變系統(tǒng)傳遞函數(shù)建立在拉普拉斯變換的基礎(chǔ)上其核心是疊加原理和常系數(shù)線性微分方程。對于非線性系統(tǒng)、時變系統(tǒng)傳遞函數(shù)矩陣的定義不再成立。而狀態(tài)空間方程? f(x, u, t), y g(x, u, t)在形式上可以描述更廣泛的系統(tǒng)。2. 丟失內(nèi)部狀態(tài)信息只反映輸入輸出關(guān)系這是傳遞函數(shù)方法最根本的局限。G(s)描述的是“黑箱”的外部特性。如果系統(tǒng)不是完全能控和完全能觀的那么傳遞函數(shù)矩陣就無法反映系統(tǒng)的全部動態(tài)。那些不能控或不能觀的模態(tài)在G(s)中會被對消掉。但從內(nèi)部看這些隱藏模態(tài)可能是不穩(wěn)定的會導(dǎo)致實際系統(tǒng)出現(xiàn)問題。狀態(tài)空間法通過直接研究狀態(tài)向量x(t)能完整揭示系統(tǒng)的內(nèi)部行為。3. 難以處理非零初始條件傳遞函數(shù)分析通常假設(shè)零初始條件。對于非零初始狀態(tài)的系統(tǒng)響應(yīng)傳遞函數(shù)方法處理起來比較麻煩。而狀態(tài)空間方程結(jié)合初始狀態(tài)x(0)可以非常自然地求解全響應(yīng)。4. 多變量系統(tǒng)分析與綜合的復(fù)雜性雖然傳遞函數(shù)矩陣將多變量關(guān)系封裝了起來但在進行控制器設(shè)計時如經(jīng)典的頻域設(shè)計法直接處理一個矩陣函數(shù)仍然非常復(fù)雜。例如如何為MIMO系統(tǒng)設(shè)計一個PID控制器每個通道單獨設(shè)計往往會因為耦合而失敗。而狀態(tài)空間法通過狀態(tài)反饋u -Kx和觀測器設(shè)計提供了系統(tǒng)化的多變量控制器設(shè)計框架如線性二次型調(diào)節(jié)器LQR、極點配置等概念上更統(tǒng)一。5. 對系統(tǒng)結(jié)構(gòu)的洞察較弱狀態(tài)空間表示法通過A, B, C, D矩陣清晰地分離了系統(tǒng)動態(tài)(A)、輸入影響(B)、輸出測量(C)和直接傳遞(D)。這種結(jié)構(gòu)化的表示更容易與物理模型對應(yīng)也便于進行能控性、能觀性等結(jié)構(gòu)性質(zhì)的分析。那么為什么我們還要學(xué)習(xí)傳遞函數(shù)矩陣因為它是連接經(jīng)典控制與現(xiàn)代控制的橋梁。它保留了頻域分析的直觀性伯德圖、奈奎斯特圖為理解系統(tǒng)的頻率響應(yīng)、帶寬、魯棒穩(wěn)定性提供了圖形化工具。許多先進的多變量頻域設(shè)計方法如H∞魯棒控制也是建立在傳遞函數(shù)矩陣的基礎(chǔ)之上。在實際工程中常常是頻域方法基于傳遞函數(shù)矩陣和時域方法基于狀態(tài)空間結(jié)合使用取長補短。因此傳遞函數(shù)矩陣是現(xiàn)代控制理論工具箱中不可或缺的一部分。它讓我們能夠用熟悉的頻域語言去理解和分析復(fù)雜的多變量系統(tǒng)同時又提醒我們其邊界所在引導(dǎo)我們走向更深刻、更強大的狀態(tài)空間分析領(lǐng)域。掌握了它你就拿到了理解和處理現(xiàn)實世界中耦合、多變量系統(tǒng)的第一把鑰匙。 從SISO到MIMO：傳遞函數(shù)矩陣的核心原理、計算與MATLAB實踐
1. 從“單輸入單輸出”到“多輸入多輸出”為什么我們需要傳遞函數(shù)矩陣如果你是從經(jīng)典控制理論比如拉普拉斯變換、伯德圖、奈奎斯特判據(jù)一路學(xué)過來的那么“傳遞函數(shù)”這個概念對你來說應(yīng)該像呼吸一樣自然。一個輸入一個輸出一個函數(shù)關(guān)系清晰明了。我們分析系統(tǒng)的穩(wěn)定性、動態(tài)性能都圍繞著這個單一的G(s)展開。但現(xiàn)實世界的系統(tǒng)尤其是現(xiàn)代工程系統(tǒng)很少有這么“單純”的。想想看一臺無人機你需要同時控制它的俯仰、橫滾、偏航和高度。你的操縱指令四個電機的轉(zhuǎn)速是多個輸入它的姿態(tài)和位置是多個輸出。一個化工反應(yīng)釜你需要控制溫度、壓力、pH值和進料流量。加熱功率、冷卻水閥、酸堿泵和進料泵是你的輸入傳感器讀數(shù)是你的輸出。一臺機械臂每個關(guān)節(jié)的電機扭矩是輸入末端執(zhí)行器的位置和姿態(tài)是輸出。這些系統(tǒng)都有一個共同點多輸入多輸出英文簡稱MIMO。在經(jīng)典控制里我們可能會嘗試為每一個輸出單獨設(shè)計一個控制器去調(diào)節(jié)某一個輸入這叫做“單回路控制”。但問題很快就來了——這些回路之間是耦合的。你調(diào)節(jié)無人機第一個電機想讓它抬頭結(jié)果它可能不僅抬頭還開始往左偏航。你加大反應(yīng)釜的加熱功率想升溫壓力可能也跟著飆升。這種耦合意味著你不能把系統(tǒng)簡單地拆成幾個獨立的單變量系統(tǒng)來處理。它們是一個整體牽一發(fā)而動全身。這時沿用單個傳遞函數(shù)的思路就捉襟見肘了。我們需要一個數(shù)學(xué)工具能夠同時描述所有輸入對所有輸出的影響關(guān)系。這個工具就是傳遞函數(shù)矩陣。你可以把它想象成一張關(guān)系網(wǎng)或者一個表格。假設(shè)系統(tǒng)有p個輸入q個輸出那么這個傳遞函數(shù)矩陣G(s)就是一個q x p的矩陣。矩陣中的每一個元素G_ij(s)都是一個傳遞函數(shù)它單獨地描述了第j個輸入對第i個輸出的影響而此時其他所有輸入都為零這是線性系統(tǒng)疊加原理的前提。所以當我們從“現(xiàn)控理論”的視角重新審視系統(tǒng)時傳遞函數(shù)矩陣就是我們連接系統(tǒng)外部描述輸入輸出關(guān)系和內(nèi)部狀態(tài)描述狀態(tài)空間方程的一座關(guān)鍵橋梁也是處理MIMO系統(tǒng)分析與設(shè)計問題的起點。它保留了傳遞函數(shù)直觀的頻率域特性又將系統(tǒng)的多變量特性用矩陣這一強大的數(shù)學(xué)語言清晰地表達了出來。2. 傳遞函數(shù)矩陣的核心定義與數(shù)學(xué)表達理解了為什么需要它我們來看看它具體是什么。我們從最通用的線性時不變系統(tǒng)的狀態(tài)空間描述出發(fā)這是現(xiàn)代控制理論的基礎(chǔ)模型狀態(tài)方程 ?(t) A x(t) B u(t) 輸出方程 y(t) C x(t) D u(t)其中x(t)是n維狀態(tài)向量。u(t)是p維輸入向量。y(t)是q維輸出向量。A, B, C, D是維數(shù)匹配的常數(shù)矩陣。我們對上述方程兩邊進行拉普拉斯變換并假設(shè)初始狀態(tài)x(0) 0這是為了聚焦于輸入輸出的傳遞特性就像經(jīng)典控制里一樣。得到sX(s) A X(s) B U(s) Y(s) C X(s) D U(s)由第一個式子解出X(s)(sI - A) X(s) B U(s) X(s) (sI - A)^{-1} B U(s)這里I是n x n的單位矩陣(sI - A)^{-1}就是預(yù)解矩陣它在系統(tǒng)分析中至關(guān)重要。將X(s)代入輸出方程Y(s) C (sI - A)^{-1} B U(s) D U(s) [C (sI - A)^{-1} B D] U(s)于是我們得到了輸入U(s)和輸出Y(s)在復(fù)頻域的關(guān)系Y(s) G(s) U(s)其中傳遞函數(shù)矩陣G(s)被定義為G(s) C (sI - A)^{-1} B D這是一個q x p的矩陣。它的每一個元素G_ij(s)都是一個關(guān)于復(fù)變量s的有理分式函數(shù)分子分母都是s的多項式。G_ij(s)表示的是當只有第j個輸入U_j(s)作用時所產(chǎn)生第i個輸出Y_i(s)的傳遞函數(shù)。一個關(guān)鍵的理解點(sI - A)^{-1}包含了系統(tǒng)的所有動態(tài)模態(tài)信息極點而矩陣C和B則決定了這些模態(tài)如何被輸入激發(fā)以及如何被輸出觀測。矩陣D代表了輸入到輸出的直接饋通在很多物理系統(tǒng)中如電路、機械D常常是零矩陣因為輸入通常不會瞬間無延遲地影響輸出。注意這里假設(shè)了(sI - A)是可逆的這要求s不是系統(tǒng)矩陣A的特征值。系統(tǒng)矩陣A的特征值正是傳遞函數(shù)矩陣G(s)的極點分母為零的點它們決定了系統(tǒng)的穩(wěn)定性。3. 如何計算與解讀傳遞函數(shù)矩陣一個詳細的計算實例定義看起來有點抽象我們通過一個具體的、簡化的雙輸入雙輸出系統(tǒng)來親手算一遍感受一下這個過程?？紤]一個耦合的彈簧-質(zhì)量塊系統(tǒng)或者一個簡單的電路網(wǎng)絡(luò)我們可以抽象出如下狀態(tài)空間模型為了計算方便數(shù)字是假設(shè)的設(shè)系統(tǒng)方程為A [ -2 1 0; 1 -3 1; 0 1 -1 ]; B [ 1 0; 0 1; 1 0 ]; C [ 1 0 0; 0 1 0 ]; D [ 0 0; 0 0 ];這里狀態(tài)x是3維輸入u是2維輸出y是2維。D矩陣為零表示沒有直接饋通。我們的目標是求出G(s) C (sI - A)^{-1} B。第一步構(gòu)造 (sI - A)sI - A [ s2 -1 0; -1 s3 -1; 0 -1 s1 ];第二步求 (sI - A) 的逆矩陣求一個3x3矩陣的逆我們可以用伴隨矩陣法(sI - A)^{-1} adj(sI - A) / det(sI - A)。 先求行列式Δ(s) det(sI - A)Δ(s) (s2)*[(s3)(s1) - (-1)(-1)] - (-1)*[(-1)(s1) - (0)(-1)] 0*... (s2)[(s3)(s1) - 1] 1*[-(s1)] (s2)(s^2 4s 3 - 1) - (s1) (s2)(s^2 4s 2) - s - 1 s^3 4s^2 2s 2s^2 8s 4 - s - 1 s^3 6s^2 9s 3所以系統(tǒng)的特征多項式是s^3 6s^2 9s 3它的根就是系統(tǒng)的極點。再求伴隨矩陣adj(sI - A)。這需要計算所有代數(shù)余子式過程略這是基本功練習(xí)。假設(shè)我們算得adj(sI - A) [ (s3)(s1)-1 (s1) 1; (s1) (s2)(s1) (s2); 1 (s2) (s2)(s3)-1 ] [ s^24s2 s1 1; s1 s^23s2 s2; 1 s2 s^25s5 ]因此(sI - A)^{-1} (1 / Δ(s)) * [ s^24s2, s1, 1; s1, s^23s2, s2; 1, s2, s^25s5 ]第三步計算 C (sI - A)^{-1} B這是一個連續(xù)的矩陣乘法。C是 2x3逆矩陣是 3x3B是 3x2。所以結(jié)果G(s)是 2x2。 我們先計算中間結(jié)果T(s) (sI - A)^{-1} B。B矩陣有兩列我們分別計算。 令B1 [1; 0; 1](B的第一列)B2 [0; 1; 0](B的第二列)。計算T1(s) (sI - A)^{-1} * B1T1(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [1; (1/Δ(s)) * [ (s^24s2)*1 (s1)*0 1*1; s1, s^23s2, s2; 0; (s1)*1 (s^23s2)*0 (s2)*1; 1, s2, s^25s5 ] 1] 1*1 (s2)*0 (s^25s5)*1 ] (1/Δ(s)) * [ s^24s2 1; s1 s2; 1 s^25s5 ] (1/Δ(s)) * [ s^24s3; 2s3; s^25s6 ]同理計算T2(s) (sI - A)^{-1} * B2T2(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [0; (1/Δ(s)) * [ (s^24s2)*0 (s1)*1 1*0; s1, s^23s2, s2; 1; (s1)*0 (s^23s2)*1 (s2)*0; 1, s2, s^25s5 ] 0] 1*0 (s2)*1 (s^25s5)*0 ] (1/Δ(s)) * [ s1; s^23s2; s2 ]所以T(s) [T1(s), T2(s)] (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2; s^25s6, s2 ]現(xiàn)在左乘C矩陣。C矩陣只取前兩行因為輸出y只對應(yīng)前兩個狀態(tài)所以G(s) C * T(s)相當于取T(s)的前兩行。G(s) (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2 ]其中Δ(s) s^3 6s^2 9s 3。解讀G(s) 這個 2x2 的傳遞函數(shù)矩陣G(s)告訴我們G_11(s) (s^24s3) / (s^36s^29s3)輸入u1對輸出y1的傳遞函數(shù)。G_12(s) (s1) / (s^36s^29s3)輸入u2對輸出y1的傳遞函數(shù)。G_21(s) (2s3) / (s^36s^29s3)輸入u1對輸出y2的傳遞函數(shù)。G_22(s) (s^23s2) / (s^36s^29s3)輸入u2對輸出y2的傳遞函數(shù)。關(guān)鍵觀察共同分母所有四個傳遞函數(shù)都有相同的分母多項式Δ(s)。這個多項式正是矩陣A的特征多項式。這意味著整個MIMO系統(tǒng)的極點決定穩(wěn)定性和基本動態(tài)是由矩陣A決定的是所有輸入輸出通道所共享的。這是MIMO系統(tǒng)與多個獨立SISO系統(tǒng)本質(zhì)不同的地方。耦合體現(xiàn)非對角線元素G_12(s)和G_21(s)不為零。這說明輸入u1會影響輸出y2輸入u2也會影響輸出y1。系統(tǒng)是耦合的。相對階G_11和G_22的分子是二階分母是三階相對階為1。G_12和G_21的分子是一階相對階為2。不同通道的動態(tài)特性可能不同。實操心得對于維數(shù)高于3的系統(tǒng)手工計算逆矩陣會非常繁瑣且容易出錯。在實際工程或?qū)W習(xí)中我們強烈依賴數(shù)學(xué)工具。在MATLAB/Octave中計算傳遞函數(shù)矩陣極其簡單。假設(shè)你已經(jīng)定義了A, B, C, D矩陣只需一行命令G ss(A, B, C, D); G_tf tf(G);。ss創(chuàng)建狀態(tài)空間模型tf將其轉(zhuǎn)換為傳遞函數(shù)矩陣形式。對于我們的例子在MATLAB中驗證上述計算是很好的練習(xí)。手工推導(dǎo)的目的在于深刻理解其數(shù)學(xué)本質(zhì)和耦合關(guān)系的來源而不是用于解決大規(guī)模問題。4. 傳遞函數(shù)矩陣的性質(zhì)與在系統(tǒng)分析中的核心作用得到了傳遞函數(shù)矩陣我們能用它做什么它不僅僅是狀態(tài)空間模型的一種等價表示更是我們分析MIMO系統(tǒng)特性的強大工具。4.1 極點與零點MIMO系統(tǒng)的擴展定義對于SISO系統(tǒng)極點就是傳遞函數(shù)分母多項式的根零點就是分子多項式的根。對于MIMO系統(tǒng)由于G(s)是一個矩陣極點和零點的定義需要擴展。極點如前所述傳遞函數(shù)矩陣G(s)的所有元素的分母多項式在約去公因子后的根就是系統(tǒng)的極點。更本質(zhì)地說系統(tǒng)的極點就是系統(tǒng)矩陣A的特征值。這一定義從狀態(tài)空間模型出發(fā)更為根本。所有輸入輸出通道共享同一組極點這決定了系統(tǒng)的穩(wěn)定性所有極點均具有負實部則漸近穩(wěn)定和基本響應(yīng)模態(tài)如振蕩頻率、衰減速度。傳輸零點這是一個MIMO系統(tǒng)特有的、非常重要的概念。它不是單個元素分子為零的點。傳輸零點的定義是存在一個非零的復(fù)頻率z和一個非零的輸入向量U_0使得在零初始狀態(tài)下系統(tǒng)的輸出Y(s)恒為零。即G(z) U_0 0。這意味著在頻率z處存在某種特定的輸入信號組合其效果被系統(tǒng)內(nèi)部完全“抵消”了無法在輸出端被觀測到。計算對于D0的系統(tǒng)傳輸零點z是使得復(fù)合矩陣P(s) [sI-A, -B; C, 0]降秩的s值。在MATLAB中可以使用tzero或zero函數(shù)直接計算。物理意義傳輸零點反映了系統(tǒng)輸入輸出之間的阻塞特性。例如在飛機控制中某些特定的舵面偏轉(zhuǎn)組合可能無法引起飛機姿態(tài)的變化在某個頻率下這個頻率就是傳輸零點。傳輸零點會影響系統(tǒng)的可控制性、可觀測性以及控制性能的極限例如非最小相位系統(tǒng)的右半平面零點會限制控制帶寬。4.2 穩(wěn)定性分析基于極點的穩(wěn)定性判據(jù)在MIMO系統(tǒng)中依然直接適用線性時不變系統(tǒng)漸近穩(wěn)定的充要條件是其傳遞函數(shù)矩陣的所有極點即矩陣A的所有特征值都具有負實部。通過傳遞函數(shù)矩陣G(s)我們可以計算其特征多項式即行列式det(sI-A)或G(s)各元素分母的最小公倍式然后應(yīng)用勞斯判據(jù)、赫爾維茨判據(jù)或直接求根來判斷穩(wěn)定性。在MATLAB中使用pole(G)或eig(A)直接獲取極點。注意對于MIMO系統(tǒng)不能僅僅通過觀察G(s)某個對角元比如G_11(s)的穩(wěn)定性來判斷整個系統(tǒng)的穩(wěn)定性。因為非對角元的耦合可能引入不穩(wěn)定的隱藏模態(tài)與狀態(tài)空間中的能控性/能觀性相關(guān)。必須檢查系統(tǒng)矩陣A的全部特征值。4.3 頻域分析從伯德圖到奇異值圖這是傳遞函數(shù)矩陣威力巨大的地方。在SISO系統(tǒng)中我們繪制單個傳遞函數(shù)G(jω)的伯德圖幅頻和相頻特性。在MIMO系統(tǒng)中G(jω)是一個復(fù)數(shù)矩陣。我們?nèi)绾卫L制它的頻率響應(yīng)逐個元素法可以為G(s)的每一個元素G_ij(jω)繪制伯德圖。這能讓我們看清每一對輸入輸出通道在頻域的特性。這對于理解耦合的強度比如G_21相對于G_22的幅值大小很有幫助。在MATLAB中bode(G)命令會自動為所有通道生成伯德圖。奇異值分析法更強大的工具對于MIMO系統(tǒng)更本質(zhì)的頻域分析工具是奇異值。對于每個頻率點ω我們計算復(fù)數(shù)矩陣G(jω)的奇異值。奇異值總是非負的實數(shù)。假設(shè)G是q x p矩陣那么它在每個頻率點有min(p, q)個奇異值記為σ_1(ω) ≥ σ_2(ω) ≥ ... ≥ σ_min(p,q)(ω) ≥ 0。最大奇異值σ_max(ω)可以理解為系統(tǒng)在頻率ω處對所有可能輸入方向向量的最大增益。最小奇異值σ_min(ω)可以理解為系統(tǒng)在頻率ω處對所有可能輸入方向的最小增益。繪制σ_max(ω)和σ_min(ω)隨頻率變化的曲線就得到了MIMO系統(tǒng)的奇異值伯德圖。為什么重要在魯棒控制和系統(tǒng)性能分析中奇異值提供了關(guān)鍵的洞察。例如σ_min在低頻段的大小反映了系統(tǒng)抗干擾和解耦的能力σ_max在高頻段的衰減速率反映了系統(tǒng)對噪聲和未建模動態(tài)的魯棒性。閉環(huán)系統(tǒng)的穩(wěn)定裕度也可以用開環(huán)傳遞函數(shù)矩陣的奇異值來評估。在MATLAB中可以使用sigma(G)命令繪制奇異值圖。4.4 能控性與能觀性與傳遞函數(shù)矩陣的關(guān)聯(lián)能控性和能觀性是狀態(tài)空間模型的核心概念它們與傳遞函數(shù)矩陣有著深刻的聯(lián)系。能控性系統(tǒng)是否能在有限時間內(nèi)通過合適的輸入u(t)將狀態(tài)從任意初始點驅(qū)動到原點。狀態(tài)空間判據(jù)是能控性矩陣[B, AB, A^2B, ..., A^{n-1}B]滿秩。能觀性系統(tǒng)是否能在有限時間內(nèi)通過輸出的觀測y(t)唯一地確定初始狀態(tài)x(0)。狀態(tài)空間判據(jù)是能觀性矩陣[C; CA; CA^2; ...; CA^{n-1}]滿秩。從傳遞函數(shù)矩陣的角度看如果系統(tǒng)是狀態(tài)空間能控且能觀的那么其傳遞函數(shù)矩陣G(s)將是既約的即沒有零極點對消。如果發(fā)生了零極點對消則意味著系統(tǒng)不是完全能控和/或完全能觀的對消掉的模態(tài)在輸入輸出描述中“消失”了但它們可能隱藏在系統(tǒng)內(nèi)部影響實際動態(tài)甚至穩(wěn)定性。在計算G(s) C(sI-A)^{-1}B D時如果最終得到的各元素傳遞函數(shù)有公因子被約去這些被約去的因子就對應(yīng)著不能控或不能觀的模態(tài)。因此在由傳遞函數(shù)矩陣反推或?qū)崿F(xiàn)狀態(tài)空間模型時必須非常小心要確保得到的是最小實現(xiàn)即完全能控且完全能觀。5. 從理論到實踐在MATLAB/SciPy中操作傳遞函數(shù)矩陣理論學(xué)習(xí)之后我們必須能在工具中實現(xiàn)它。這里以MATLAB及其開源替代品GNU Octave和Python的SciPy庫為例展示基本操作。5.1 MATLAB/Octave 環(huán)境1. 定義系統(tǒng)與計算傳遞函數(shù)矩陣% 定義狀態(tài)空間矩陣 (沿用之前的例子) A [-2, 1, 0; 1, -3, 1; 0, 1, -1]; B [1, 0; 0, 1; 1, 0]; C [1, 0, 0; 0, 1, 0]; D zeros(2,2); % 創(chuàng)建狀態(tài)空間模型對象 sys_ss ss(A, B, C, D); % 轉(zhuǎn)換為傳遞函數(shù)矩陣形式 sys_tf tf(sys_ss); disp(傳遞函數(shù)矩陣 G(s):); disp(sys_tf); % 輸出會顯示一個2x2的tf數(shù)組每個元素都是tf對象。 % 例如從輸入 1 到輸出 1: s^2 4 s 3 / s^3 6 s^2 9 s 3 % 從輸入 2 到輸出 1: s 1 / s^3 6 s^2 9 s 3 ... 等等2. 提取特定通道的傳遞函數(shù)% 提取輸入1到輸出1的SISO傳遞函數(shù) G11 sys_tf(1, 1); % 提取輸入2到輸出2的SISO傳遞函數(shù) G22 sys_tf(2, 2);3. 計算極點和零點% 計算系統(tǒng)極點 (即A的特征值) poles pole(sys_ss); % 或 pole(sys_tf) disp(系統(tǒng)極點); disp(poles); % 計算系統(tǒng)傳輸零點 zeros tzero(sys_ss); % 或 zero(sys_tf) disp(系統(tǒng)傳輸零點); disp(zeros);4. 頻域分析繪圖% 繪制所有通道的伯德圖 (4個子圖) figure; bode(sys_ss); % 或 bode(sys_tf) grid on; title(各通道伯德圖); % 繪制奇異值圖 (MIMO系統(tǒng)頻域分析利器) figure; sigma(sys_ss); % 或 sigma(sys_tf) grid on; title(系統(tǒng)奇異值圖);5. 時域仿真% 定義時間向量 t 0:0.01:10; % 定義輸入信號u1是階躍u2是正弦波 u [ones(size(t)); sin(2*t)]; % 注意維度是 2 x length(t) % 進行仿真 [y, t_out, x] lsim(sys_ss, u, t); % 繪制輸出響應(yīng) figure; subplot(2,1,1); plot(t_out, y(:,1)); ylabel(y1); grid on; legend(輸出1); subplot(2,1,2); plot(t_out, y(:,2)); ylabel(y2); xlabel(時間 (s)); grid on; legend(輸出2);5.2 Python (SciPy Matplotlib) 環(huán)境Python在科學(xué)計算領(lǐng)域應(yīng)用廣泛控制庫雖不如MATLAB專業(yè)但基礎(chǔ)功能完備。import numpy as np import matplotlib.pyplot as plt from scipy import signal # 1. 定義狀態(tài)空間矩陣 A np.array([[-2, 1, 0], [1, -3, 1], [0, 1, -1]]) B np.array([[1, 0], [0, 1], [1, 0]]) C np.array([[1, 0, 0], [0, 1, 0]]) D np.array([[0, 0], [0, 0]]) # 2. 創(chuàng)建狀態(tài)空間系統(tǒng) sys_ss signal.StateSpace(A, B, C, D) # 3. 轉(zhuǎn)換為傳遞函數(shù)形式 (注意scipy.signal 的 tf 表示是針對SISO的) # 對于MIMO我們需要逐個元素計算或使用專門的庫如 control需安裝pip install control # 這里演示使用 control 庫 (如果已安裝) try: import control as ct # 創(chuàng)建控制系統(tǒng)庫的狀態(tài)空間對象 sys_ct ct.ss(A, B, C, D) # 計算傳遞函數(shù)矩陣 (返回一個 TransferFunction 對象列表的列表) sys_tf ct.tf(sys_ct) print(傳遞函數(shù)矩陣:) print(sys_tf) except ImportError: print(control 庫未安裝。對于MIMO傳遞函數(shù)建議安裝 control 庫。) # 手動計算其中一個通道作為示例 (例如 G11) # G(s) C(sI-A)^-1 B D, 我們計算第一個元素 # 這里省略手動計算逆矩陣的代碼較為復(fù)雜。 # 作為替代我們可以直接進行時域/頻域仿真。 # 4. 計算極點和零點 (使用 scipy.signal 或 control) # 極點就是 A 的特征值 poles np.linalg.eigvals(A) print(\n系統(tǒng)極點 (A的特征值):) print(poles) # 傳輸零點計算較復(fù)雜control庫有提供 try: import control as ct zeros ct.zero(sys_ct) print(\n系統(tǒng)傳輸零點 (來自 control 庫):) print(zeros) except: print(\n傳輸零點計算需要 control 庫。) # 5. 時域仿真 (使用 scipy.signal) t np.arange(0, 10, 0.01) # 輸入u1 階躍 u2 正弦波 u np.vstack([np.ones_like(t), np.sin(2*t)]) # 形狀 (2, len(t)) # lsim 需要輸入 u 的形狀為 (len(t), num_inputs) t_out, y, x signal.lsim(sys_ss, Uu.T, Tt) # 注意 u.T 轉(zhuǎn)置 # 繪制結(jié)果 fig, (ax1, ax2) plt.subplots(2, 1, figsize(10, 6)) ax1.plot(t_out, y[:, 0]) ax1.set_ylabel(y1) ax1.grid(True) ax1.legend([輸出1]) ax2.plot(t_out, y[:, 1]) ax2.set_ylabel(y2) ax2.set_xlabel(時間 [s]) ax2.grid(True) ax2.legend([輸出2]) plt.suptitle(系統(tǒng)時域響應(yīng)) plt.tight_layout() plt.show() # 6. 頻域分析 - 伯德圖 (scipy.signal 的 bode 是 SISO 的) # 我們可以循環(huán)繪制每個通道 try: import control as ct plt.figure() ct.bode_plot(sys_ct) # control庫可以繪制MIMO伯德圖 plt.suptitle(各通道伯德圖 (control庫)) plt.tight_layout() plt.show() except: print(繪制MIMO伯德圖需要 control 庫。)實操心得與避坑指南工具選擇對于嚴肅的控制系統(tǒng)分析與設(shè)計MATLAB及其控制系統(tǒng)工具箱仍然是行業(yè)標桿文檔齊全、函數(shù)豐富。Python的control庫功能日益完善對于學(xué)習(xí)、研究和輕量級應(yīng)用是完全足夠的且免費開源。SciPy.signal更側(cè)重于信號處理對MIMO支持較弱。維度匹配在定義矩陣A, B, C, D時務(wù)必反復(fù)檢查維度。(n x n), (n x p), (q x n), (q x p)。一個常見的錯誤是B或C矩陣的維度定義反了。仿真輸入使用lsim或signal.lsim進行仿真時輸入矩陣U的維度容易出錯。在MATLAB中U的列數(shù)等于輸入數(shù)行數(shù)等于時間點數(shù)量。在SciPy的signal.lsim中U的形狀是(len(t), num_inputs)。務(wù)必查閱文檔或打印矩陣形狀來確認。最小實現(xiàn)如果你是從物理方程推導(dǎo)或系統(tǒng)辨識得到的傳遞函數(shù)矩陣在轉(zhuǎn)換為狀態(tài)空間模型時使用minreal函數(shù)MATLAB或minreal方法control庫來獲取最小實現(xiàn)消除不能控或不能觀的模態(tài)避免后續(xù)分析與設(shè)計出現(xiàn)問題。6. 傳遞函數(shù)矩陣的局限性與狀態(tài)空間法的優(yōu)勢盡管傳遞函數(shù)矩陣是分析MIMO系統(tǒng)的有力工具但它并非萬能也存在固有的局限性。理解這些局限性正是我們深入學(xué)習(xí)現(xiàn)代控制理論狀態(tài)空間法的動力。1. 僅適用于線性時不變系統(tǒng)傳遞函數(shù)建立在拉普拉斯變換的基礎(chǔ)上其核心是疊加原理和常系數(shù)線性微分方程。對于非線性系統(tǒng)、時變系統(tǒng)傳遞函數(shù)矩陣的定義不再成立。而狀態(tài)空間方程? f(x, u, t), y g(x, u, t)在形式上可以描述更廣泛的系統(tǒng)。2. 丟失內(nèi)部狀態(tài)信息只反映輸入輸出關(guān)系這是傳遞函數(shù)方法最根本的局限。G(s)描述的是“黑箱”的外部特性。如果系統(tǒng)不是完全能控和完全能觀的那么傳遞函數(shù)矩陣就無法反映系統(tǒng)的全部動態(tài)。那些不能控或不能觀的模態(tài)在G(s)中會被對消掉。但從內(nèi)部看這些隱藏模態(tài)可能是不穩(wěn)定的會導(dǎo)致實際系統(tǒng)出現(xiàn)問題。狀態(tài)空間法通過直接研究狀態(tài)向量x(t)能完整揭示系統(tǒng)的內(nèi)部行為。3. 難以處理非零初始條件傳遞函數(shù)分析通常假設(shè)零初始條件。對于非零初始狀態(tài)的系統(tǒng)響應(yīng)傳遞函數(shù)方法處理起來比較麻煩。而狀態(tài)空間方程結(jié)合初始狀態(tài)x(0)可以非常自然地求解全響應(yīng)。4. 多變量系統(tǒng)分析與綜合的復(fù)雜性雖然傳遞函數(shù)矩陣將多變量關(guān)系封裝了起來但在進行控制器設(shè)計時如經(jīng)典的頻域設(shè)計法直接處理一個矩陣函數(shù)仍然非常復(fù)雜。例如如何為MIMO系統(tǒng)設(shè)計一個PID控制器每個通道單獨設(shè)計往往會因為耦合而失敗。而狀態(tài)空間法通過狀態(tài)反饋u -Kx和觀測器設(shè)計提供了系統(tǒng)化的多變量控制器設(shè)計框架如線性二次型調(diào)節(jié)器LQR、極點配置等概念上更統(tǒng)一。5. 對系統(tǒng)結(jié)構(gòu)的洞察較弱狀態(tài)空間表示法通過A, B, C, D矩陣清晰地分離了系統(tǒng)動態(tài)(A)、輸入影響(B)、輸出測量(C)和直接傳遞(D)。這種結(jié)構(gòu)化的表示更容易與物理模型對應(yīng)也便于進行能控性、能觀性等結(jié)構(gòu)性質(zhì)的分析。那么為什么我們還要學(xué)習(xí)傳遞函數(shù)矩陣因為它是連接經(jīng)典控制與現(xiàn)代控制的橋梁。它保留了頻域分析的直觀性伯德圖、奈奎斯特圖為理解系統(tǒng)的頻率響應(yīng)、帶寬、魯棒穩(wěn)定性提供了圖形化工具。許多先進的多變量頻域設(shè)計方法如H∞魯棒控制也是建立在傳遞函數(shù)矩陣的基礎(chǔ)之上。在實際工程中常常是頻域方法基于傳遞函數(shù)矩陣和時域方法基于狀態(tài)空間結(jié)合使用取長補短。因此傳遞函數(shù)矩陣是現(xiàn)代控制理論工具箱中不可或缺的一部分。它讓我們能夠用熟悉的頻域語言去理解和分析復(fù)雜的多變量系統(tǒng)同時又提醒我們其邊界所在引導(dǎo)我們走向更深刻、更強大的狀態(tài)空間分析領(lǐng)域。掌握了它你就拿到了理解和處理現(xiàn)實世界中耦合、多變量系統(tǒng)的第一把鑰匙。