簡介上海交大儀器系研究生課程“線性系統(tǒng)分析與設(shè)計”的期末大作業(yè)以MATLAB為主要實現(xiàn)工具圍繞連續(xù)/離散線性系統(tǒng)的建模、分析與校正展開。整個項目打包為rar格式共3個文件包括2個m腳本與1份docx報告壓縮后僅193KB體量輕但內(nèi)容完整便于快速檢索。m代碼側(cè)重系統(tǒng)建模與仿真覆蓋傳遞函數(shù)構(gòu)建、Routh-Hurwitz與根軌跡穩(wěn)定性分析、控制器參數(shù)整定以及階躍/脈沖響應(yīng)對比等典型環(huán)節(jié)docx報告則詳細(xì)闡述了問題背景、設(shè)計流程和結(jié)果分析代碼與文字相互印證便于復(fù)現(xiàn)和舉一反三。已有680人學(xué)習(xí)下載對于正在修讀線性系統(tǒng)理論、需要參考完整MATLAB課程設(shè)計或希望熟悉Control System Toolbox的讀者是一份很實用的學(xué)習(xí)素材。1. 從公式到能跑的Matlab仿真這份線性系統(tǒng)大作業(yè)到底解決什么問題拿到某高校線性系統(tǒng)課程的大作業(yè)題目時我手里只有幾個狀態(tài)矩陣、一個“設(shè)計狀態(tài)反饋使系統(tǒng)穩(wěn)定”的需求和一張A4紙的評分標(biāo)準(zhǔn)。真正動手才發(fā)現(xiàn)公式推導(dǎo)是一回事把這些公式變成一段能在Matlab里跑通、能畫圖、能應(yīng)付答辯演示的完整程序完全是另一回事。這份課程大作業(yè)資源就是這么一套東西把狀態(tài)空間建模、能控能觀性分析、狀態(tài)反饋、觀測器設(shè)計和離散仿真串成了一條完整鏈路適合正被作業(yè)卡住的學(xué)生和想快速上手狀態(tài)空間方法的工程師。下文按復(fù)現(xiàn)順序講每段代碼可獨立執(zhí)行。2. 建模與狀態(tài)空間化寫Matlab之前先想清楚這三件事狀態(tài)反饋聽起來高大上但第一步是把被控對象變成A、B、C、D四個矩陣。這一章先把對象模型選型講清楚再落到Matlab代碼最后說采樣周期的坑——這三件事沒理順后面算出來的K再漂亮也是空中樓閣。2.1 為什么是帶執(zhí)行機構(gòu)動態(tài)的三階模型而不是教材雙積分器很多課本例題喜歡用二階雙積分器因為數(shù)學(xué)上清爽能控能觀矩陣看一眼就能求畫個根軌跡就收工。但真實的大作業(yè)場景里這種模型太“嬌貴”它只有兩個積分環(huán)節(jié)狀態(tài)反饋一加上去閉環(huán)特征方程簡單到讓人看不出設(shè)計空間。實際上課程大作業(yè)為了考察完整的現(xiàn)代控制流程通常會引入執(zhí)行機構(gòu)動態(tài)把對象變成三階也就是雙積分器串聯(lián)一個一階慣性環(huán)節(jié)。我用的對象模型是G(s) 1 / [s2(τs 1)]其中時間常數(shù)τ1s代表一個響應(yīng)偏慢的執(zhí)行機構(gòu)。選這個結(jié)構(gòu)出于三點考慮第一它包含積分環(huán)節(jié)能模擬位置伺服系統(tǒng)這類真實對象后邊算狀態(tài)反饋時能直接對應(yīng)到“控制器輸出到執(zhí)行器再到機械位置”的物理鏈條第二三階系統(tǒng)開環(huán)不穩(wěn)定必須靠反饋拉穩(wěn)能明顯體現(xiàn)狀態(tài)反饋的作用第三三個狀態(tài)對應(yīng)位置、速度、執(zhí)行機構(gòu)輸出答辯時每個狀態(tài)都能講出物理意義不至于被追問時露怯。對應(yīng)的狀態(tài)空間表示如下狀態(tài)變量定義為 x?位置x?速度x?執(zhí)行機構(gòu)輸出的加速度指令? Ax Bu y Cx Du矩陣取值矩陣取值含義A[0 1 0; 0 0 1; 0 0 -1]位置→速度、速度→加速度、執(zhí)行機構(gòu)一階滯后B[0; 0; 1]控制電壓直接作用在執(zhí)行機構(gòu)輸入端C[1 0 0]測量輸出為位置D0無直接饋通這里x?是執(zhí)行機構(gòu)輸出不是加速度本身也能講但把它當(dāng)成“經(jīng)過一階滯后后的控制作用”最直觀u經(jīng)過一個時間常數(shù)為1s的慣性環(huán)節(jié)后變成x?再積分兩次得到位置。這種狀態(tài)選取不是唯一答案但它是與物理過程對應(yīng)最直接的一組后邊設(shè)計觀測器時也容易設(shè)置初值。2.2 把題干寫進(jìn)Matlab狀態(tài)空間對象與交叉驗證建模這一步我習(xí)慣直接寫矩陣而不是用傳遞函數(shù)轉(zhuǎn)換因為矩陣寫出來能自己核對每一行物理意義。最小工作代碼% 定義三階對象模型 A [0 1 0; 0 0 1; 0 0 -1]; % 系統(tǒng)矩陣最后一行為執(zhí)行機構(gòu)動態(tài) B [0; 0; 1]; % 輸入矩陣控制量作用在第3個狀態(tài) C [1 0 0]; % 輸出矩陣測量位置x1 D 0; % 直接饋通為零 % 生成連續(xù)時間狀態(tài)空間對象 sys ss(A, B, C, D); % 檢查維度一致性并顯示對象信息 n size(A, 1); % 狀態(tài)數(shù)應(yīng)為3 m size(B, 2); % 輸入數(shù)應(yīng)為1 p size(C, 1); % 輸出數(shù)應(yīng)為1 disp([系統(tǒng)階數(shù) n , num2str(n)]);這里的size檢查很多人會跳過但我建議保留后續(xù)代碼里如果出現(xiàn)B矩陣維度寫錯比如忘了給第三行寫1ss會因為維度不匹配直接報錯到時候排錯比這多花十分鐘。sys對象建立后后續(xù)階躍響應(yīng)、bode圖、c2d離散化都可以直接復(fù)用它不需要重復(fù)寫矩陣。除了直接寫矩陣我還會做一個交叉驗證把同一對象寫成傳遞函數(shù)再用tf2ss轉(zhuǎn)換對比特征根是否一致。這一步不能證明狀態(tài)選取唯一但能確認(rèn)建模沒有筆誤。% 用傳遞函數(shù)交叉驗證 G tf(1, [1 1 0 0]); % 分母: s^3 s^2對應(yīng) 1/(s^2(s1)) sys_tf ss(G); % 轉(zhuǎn)換為狀態(tài)空間 disp(eig(sys.A)); % 應(yīng)該為 0, 0, -1 disp(eig(sys_tf.A)); % 與上行一致即可說明一下tf的分子為1分母寫成 [1 1 0 0]對應(yīng) s3s2展開就是 s2(s1)。eig比較的是狀態(tài)矩陣A的特征值正常情況下都是0、0、-1。如果兩次結(jié)果不一致說明矩陣符號或行序有誤。這種“從傳遞函數(shù)換算一遍”的交叉驗證是建模階段最便宜的后悔藥。2.3 連續(xù)域還是離散域采樣周期的約定狀態(tài)反饋算出來之后很多人習(xí)慣直接在連續(xù)域里跑Simulink然后把控制器原樣搬進(jìn)數(shù)字系統(tǒng)結(jié)果走了個大彎路。我的做法是一開始就定好兩個仿真層級連續(xù)域用于驗證控制算法的數(shù)學(xué)正確性離散域用于驗證工程可實現(xiàn)性。兩者差一個環(huán)節(jié)零階保持器。離散化這一步的核心參數(shù)是采樣周期Ts。對三階對象我一般先按系統(tǒng)帶寬的10到20倍來選。以本對象為例開環(huán)有一個零特征根閉環(huán)期望帶寬大約在25rad/s那么Ts取0.01s到0.05s都是合理的。太大會引入明顯相位滯后導(dǎo)致離散控制器與連續(xù)設(shè)計脫節(jié)太小則仿真步數(shù)多運行慢且沒有實際意義。% 采樣周期選擇與離散化 Ts 0.01; % 采樣周期10ms工程中常見 sys_d c2d(sys, Ts, zoh); % 零階保持離散化 % 觀察離散能控性矩陣是否保持滿秩 Co_d ctrb(sys_d.A, sys_d.B); disp([離散系統(tǒng)能控性秩: , num2str(rank(Co_d))]);關(guān)于zoh和foh的選擇控制系統(tǒng)仿真里默認(rèn)用zoh因為它模擬的是真實數(shù)字控制器的信號保持行為。foh在某些信號重建場景更平滑但大作業(yè)沒必要引入這種復(fù)雜度。離散化后必須重新檢查能控性這一點與連續(xù)系統(tǒng)不是一回事具體原因在第3章末尾細(xì)說。3. 能控能觀性檢驗rank命令算出1之前先理解矩陣在測什么狀態(tài)反饋設(shè)計的前提是系統(tǒng)能控觀測器設(shè)計的前提是系統(tǒng)能觀。這一步教材上有現(xiàn)成判據(jù)但Matlab里rank函數(shù)的結(jié)果經(jīng)常讓人摸不著頭腦明明數(shù)學(xué)上滿秩顯示出來卻是2。這章先把矩陣判據(jù)講透再給出帶容差處理的實操代碼。3.1 能控矩陣的列向量到底在測什么能控性矩陣定義是 Co [B AB A2B ...]對三階系統(tǒng)就是 [B AB A2B]。每一列代表的含義B列是控制量直接作用的狀態(tài)方向AB列是控制量經(jīng)過一步系統(tǒng)動態(tài)后影響的狀態(tài)方向A2B列是經(jīng)過兩步后影響的方向。如果這個矩陣秩等于狀態(tài)數(shù)n說明控制量通過最多n步的動態(tài)傳播能覆蓋整個狀態(tài)空間。反過來秩不足就意味著存在某個狀態(tài)方向永遠(yuǎn)不在控制量的可達(dá)范圍內(nèi)。物理上常表現(xiàn)為某個狀態(tài)變量與控制通道之間存在結(jié)構(gòu)性的解耦是零極點對消造成的。比如對象傳遞函數(shù)分子上有一個與系統(tǒng)極點相同的零點那個模態(tài)就被輸入“隔離”了。能觀性矩陣 Ob [C; CA; CA2]測的是另一個方向從輸出y的當(dāng)前值和各階導(dǎo)數(shù)值中能否反推出全部狀態(tài)初始值。對三階系統(tǒng)如果Ob滿秩就能從位置測量值及其一階、二階導(dǎo)數(shù)信息中恢復(fù)速度與執(zhí)行機構(gòu)狀態(tài)。這在后邊設(shè)計觀測器時是必要前提。3.2 Matlab實操判定代碼與閾值容差直接寫rank(ctrb(sys.A, sys.B))是最常見的做法但我不推薦只用這一句交差。原因是Matlab的rank默認(rèn)帶了一個與機器精度相關(guān)的容差默認(rèn)約等于 max(size(A)) * eps * norm(A)。若系統(tǒng)矩陣數(shù)值剛好讓某個奇異值落在容差附近rank會返回一個偏小的值造成誤判。更穩(wěn)妥的寫法是同時計算秩和最小奇異值并做一個顯式的閾值判斷% 能控能觀性完整檢驗 Co ctrb(A, B); % 能控性矩陣 3x3 Ob obsv(A, C); % 能觀性矩陣 3x3 % 分別計算奇異值 s_Co svd(Co); s_Ob svd(Ob); % 相對閾值小于最大奇異值的1e-10則視為數(shù)值零 tol 1e-10; r_Co sum(s_Co max(s_Co) * tol); r_Ob sum(s_Ob max(s_Ob) * tol); disp([能控性秩: , num2str(r_Co), / , num2str(size(A,1))]); disp([能觀性秩: , num2str(r_Ob), / , num2str(size(A,1))]);這里tol用的是相對閾值比rank默認(rèn)容差更直觀理解如果一個奇異值比最大的小10個數(shù)量級以上說明對應(yīng)方向上的能控/能觀性弱到工程上可以忽略。對本文三階模型Co和Ob的理論秩均為3上面代碼輸出應(yīng)為3/3和3/3。如果你計算出的最小奇異值在1e-3量級而最大奇異值是1e2比值1e-5系統(tǒng)雖然數(shù)學(xué)滿秩但工程意義上已經(jīng)接近不可控控制器增益會很大——這是后邊第5章要展開的坑。3.3 Gram矩陣與PBH判據(jù)答辯備用驗證當(dāng)rank結(jié)果模棱兩可時另一個更平滑的判斷工具是能控性格拉姆矩陣。連續(xù)系統(tǒng)能控性等價于格拉姆矩陣正定它能給出“控制能量”的量化視角最小特征值小說明讓某個狀態(tài)方向收斂需要很大的控制能量。% 能控性格拉姆矩陣 Wc gram(sys, c); disp(eig(Wc)); % 三個特征值均應(yīng)顯著大于0 % PBH判據(jù): 對每個特征值 lambda檢查 [lambda*I-A B] 行滿秩 lambda eig(A); for i 1:length(lambda) M [lambda(i)*eye(n) - A, B]; if rank(M) n disp([特征值 , num2str(lambda(i)), 不可控]); end endPBH判據(jù)的實用價值在于定位不可控模態(tài)哪條特征值對應(yīng)的矩陣缺秩就是哪個模態(tài)失控。應(yīng)付答辯時老師問“你怎么證明極點都能任意配置”用PBH判據(jù)作答比只貼rank結(jié)果更有說服力因為結(jié)論直接與特征值綁定。末尾提醒一個容易忽略的點系統(tǒng)在連續(xù)域能控離散化后不一定能控。當(dāng)采樣周期Ts恰好使某個連續(xù)模態(tài)變成離散域中的采樣零點離散能控矩陣秩就會下降。這就是上一章結(jié)尾我堅持要重算ctrb(sys_d.A, sys_d.B)的原因。對本文Ts0.01s不存在這個問題但若你把Ts換成1s會發(fā)現(xiàn)原本可控的系統(tǒng)離散后接近失控。4. 狀態(tài)反饋與觀測器極點怎么選、增益怎么算、分離定理怎么用能控能觀滿秩后真正的大作業(yè)重頭戲才開始極點配置。大多數(shù)人卡在不是不會調(diào)acker而是不知道怎么選極點以及為什么選了觀測器極點后整個系統(tǒng)就穩(wěn)定了。這章把選擇依據(jù)、代碼實現(xiàn)和分離定理串在一起講。4.1 極點位置與系統(tǒng)響應(yīng)的權(quán)衡狀態(tài)反饋的本質(zhì)是把閉環(huán)極點移動到設(shè)計者指定的位置。極點實部絕對值越大響應(yīng)越快但控制量峰值也越大虛部不為零時產(chǎn)生振鈴阻尼比越小超調(diào)越大。我整理三組參考極點極點組合阻尼比/形式預(yù)期效果控制量特點-2±2j, -50.707超調(diào)約5%響應(yīng)中等K元素適中執(zhí)行器壓力小-3±1.5j, -60.894超調(diào)3%以內(nèi)調(diào)節(jié)時間更短K元素較大控制量峰值翻倍-1±3j, -40.316明顯振蕩第一個超調(diào)可能上20%反饋增益偏小但系統(tǒng)動態(tài)差初學(xué)者最容易犯的錯是把極點全部設(shè)在實軸上并且離原點很遠(yuǎn)比如[-10, -11, -12]。這樣閉環(huán)響應(yīng)確實“快”但K矩陣各分量常常達(dá)到幾十上百控制量遠(yuǎn)超執(zhí)行器限幅仿真圖上響應(yīng)快得離譜實際系統(tǒng)一臺機都動不了。我一般先按阻尼比0.707、自然頻率23rad/s來定主導(dǎo)極點第三個非主導(dǎo)極點取-5保證距離主導(dǎo)極點約2倍以上避免第三極點影響主導(dǎo)模式。4.2 用acker還是place計算狀態(tài)反饋增益的姿勢Matlab提供兩個核心函數(shù)acker單輸入系統(tǒng)和place通用多輸入。對單輸入大作業(yè)acker夠用且代碼短但它在極點非常接近時數(shù)值上容易崩。place在MIMO和重極點場景更穩(wěn)但要求期望極點不能重復(fù)得太嚴(yán)苛。我的習(xí)慣是先acker出現(xiàn)數(shù)值告警立刻切place。% 期望極點主導(dǎo)極點阻尼0.707自然頻率3rad/s p_dom -3*0.707 1i*3*sqrt(1-0.707^2); % -2.12 2.12j p3 -5; % 非主導(dǎo)極點 poles [p_dom; conj(p_dom); p3]; % 狀態(tài)反饋增益計算 K acker(A, B, poles); % 單輸入系統(tǒng)直接算行向量 % K place(A, B, poles); % 若出現(xiàn)重極點警告時換成這個 % 閉環(huán)穩(wěn)定性驗證 Acl A - B*K; disp(閉環(huán)特征值:); disp(eig(Acl)); % 檢查K元素幅度是否合理 disp([反饋增益 max(abs(K)) , num2str(max(abs(K)))]);代碼中conj(p_dom)是補齊共軛極點這點不能漏否則特征多項式系數(shù)是復(fù)數(shù)物理不可實現(xiàn)。K的計算結(jié)果是一個1×3的行向量對應(yīng)三個狀態(tài)各自的反饋通道。注意eig(Acl)必須與poles一致差一位說明符號寫錯A-BK和ABK只差一個負(fù)號輸出完全不同。K元素幅度是常被忽略的量級檢查。對本文模型合理期望是max(abs(K))不超過20。如果你算出來K里有三位數(shù)多半是極點放得過遠(yuǎn)或系統(tǒng)病態(tài)回頭重新權(quán)衡響應(yīng)速度。4.3 觀測器設(shè)計與分離定理的工程含義狀態(tài)反饋設(shè)計假定全部狀態(tài)可測但大作業(yè)通常只會告訴你輸出是位置速度和執(zhí)行機構(gòu)狀態(tài)要靠狀態(tài)觀測器重建。觀測器增益L的求法本質(zhì)上與K對稱只不過把極點配置用在誤差動態(tài)上% 觀測器極點比反饋極點快3倍左右 obs_poles 3 * real(poles); % 實部放大3倍 obs_poles obs_poles - 1i * imag(poles)*0.5; % 適當(dāng)減小虛部防噪聲放大 % 求觀測器增益列向量 L place(A, C, obs_poles); % 轉(zhuǎn)置是經(jīng)典坑少一個轉(zhuǎn)置直接算錯 % 驗證觀測器誤差動態(tài) A_obs A - L*C; disp(觀測器誤差特征值:); disp(eig(A_obs));place用于觀測器時第一個參數(shù)是A第二個是C算完結(jié)果再轉(zhuǎn)置回來得到L。原因是誤差動態(tài)方程的特征值與A-LC的轉(zhuǎn)置A-CL相同寫成轉(zhuǎn)置形式才能復(fù)用place接口。obs_poles取反饋極點的3倍經(jīng)驗上能保證觀察誤差收斂速度明顯快于控制器響應(yīng)但又不會把測量噪聲放大到不可接受。分離定理說組合系統(tǒng)的閉環(huán)特征值等于反饋極點與觀測誤差極點的并集設(shè)計上可分別獨立進(jìn)行。為了向自己證明這一點可以把兩個子系統(tǒng)拼成四階組合矩陣驗證% 組合系統(tǒng)矩陣驗證分離定理 A_full [A - B*K, B*K; zeros(n,n), A - L*C]; eig_full eig(A_full); % 分離后對比 eig_sep [eig(A - B*K); eig(A - L*C)]; disp(組合特征值誤差:); disp(max(abs(sort(eig_full) - sort(eig_sep))));數(shù)值上兩者必須逐位相同誤差通常在1e-12左右。這個大括號里的驗證代碼答辯時直接展示比空口解釋分離定理有說服力得多。4.4 備選方案時間緊直接上LQR如果不想花時間挑極點LQR是一條兜底路線。LQR本質(zhì)上是通過最小化二次型代價函數(shù)來確定K不需要指定極點位置但需要設(shè)定Q和R矩陣。我建議Q取狀態(tài)權(quán)重R取控制權(quán)重% LQR 設(shè)計備選 Q_lqr C*C 0.1*eye(n); % 輸出權(quán)重輕微狀態(tài)懲罰 R_lqr 1; % 控制能量權(quán)重 [K_lqr, ~, ~] lqr(A, B, Q_lqr, R_lqr); disp(LQR反饋增益:); disp(K_lqr);Q矩陣中C*C是輸出項yy的擴展意思是重點懲罰輸出偏差加0.1eye(n)是讓三個狀態(tài)均有一定懲罰防止某個狀態(tài)被“無視”。R越小控制越激進(jìn)R越大增益越小響應(yīng)越慢。LQR的好處是不用擔(dān)心極點選錯缺點是答辯時難以講清“為什么Q取0.1而不是0.2”所以做演示我還是推薦用ackerplace的主線LQR留作補充對比。5. 大作業(yè)最常翻車的五個坑仿真發(fā)散不等于你算錯了把狀態(tài)反饋和觀測器放進(jìn)仿真之前先把自己從“下一步做什么”切換到“哪里可能出錯”。這五個坑是我反復(fù)在類似任務(wù)里遇到的每一條都真實發(fā)生過且每條都有明確的排查路線。5.1 現(xiàn)象rank(ctrb(A,B)) 顯示2但理論上明明滿秩查代碼半天發(fā)現(xiàn)矩陣沒寫錯能控性矩陣算出來也是3×3可是rank返回2。原因在于Matlab的rank默認(rèn)容差是相對機器精度的當(dāng)A矩陣較大或狀態(tài)變量間數(shù)值相差懸殊時最小奇異值落在容差附近會被當(dāng)成零。解決方法是顯式指定容差或像第3章那樣用svd看完整奇異值分布rank(ctrb(A,B), 1e-8) % 顯式指定容差 % 或 s svd(ctrb(A,B)); % 看全部奇異值 disp(min(s) / max(s)); % 比值小于1e-12再考慮判定不可控如果比值在1e-4到1e-8之間雖然rank滿秩也要小心對應(yīng)的狀態(tài)方向能控性很弱反饋增益會異常大。這個坑的教訓(xùn)是矩陣判據(jù)在紙上是個0/1問題在Matlab里是個數(shù)值精度問題不能只看rank返回值。5.2 現(xiàn)象acker 報錯“極點太近”place 返回的 K 大得離譜直接套用期望極點集合[-5, -5.0001, -6]時acker會因極點接近無法求解place能用但返回的K分分鐘上百。原因是極點接近會使得特征多項式在所求區(qū)域ill-conditioned閉環(huán)系統(tǒng)對參數(shù)擾動極度敏感。解決思路是把極點距離拉開要么均勻分布在一條圓弧上要么保留至少0.51的實部間隔。% 拉開間距的極點設(shè)置 r 3; % 半徑 p r * exp(1i*(pi (0:2)*2*pi/3)); % 在半徑3的圓上均布三個極點 K place(A, B, p); disp(max(abs(K)));均布極點的好處是任意兩個極點之間角度差120°最小間距大數(shù)值求解穩(wěn)定。算完務(wù)必看一眼max(abs(K))超過30就該減小半徑重新來。別為了追求“理想超調(diào)”犧牲控制量。5.3 現(xiàn)象Simulink 一運行就刷 NaN頻率高到屏幕閃仿真發(fā)散最常見原因不是公式錯而是初始條件不一致觀測器初值取0但真實系統(tǒng)初值在[1;0;0]誤差動態(tài)雖穩(wěn)定但初值誤差大觸發(fā)非線性環(huán)節(jié)數(shù)值溢出。另一個高頻原因是A-B*K本身不穩(wěn)定極點沒配上或狀態(tài)反饋加進(jìn)了Simulink時把K矩陣的乘號順序?qū)懛础Ｅ挪槁窂桨错樞蜃鱿仍诿钚信躶tep(ss(A-B*K, B, C, D))如果這里也發(fā)散是極點配置的問題先回頭查K如果這里穩(wěn)定而Simulink發(fā)散再去Simulink里核對State Space模塊的A、B、C、D參數(shù)和初值。我的習(xí)慣是Simulink里的狀態(tài)空間初值與真實系統(tǒng)設(shè)成一致觀測器初值單獨設(shè)到另一個Constant模塊避免手滑。5.4 現(xiàn)象離散化后控制電壓高頻抖動連續(xù)域仿真卻一切正常連續(xù)域仿真里K算好、響應(yīng)平滑換上c2d離散對象后控制量u-K*x在每個采樣點高頻振蕩。原因是采樣周期Ts太大離散控制器的相位滯后讓閉環(huán)在高頻段失去穩(wěn)定裕度另一個常見原因是觀測器極點選得太快比如取反饋極點的10倍以上導(dǎo)致離散化后觀測器放大高頻噪聲經(jīng)K反饋后出現(xiàn)抖振。解決方向把Ts從0.05s降到0.01s或0.005s觀測器極點與反饋極點的倍數(shù)控制在35倍以內(nèi)并給測量輸出加一個一階低通濾波器。離散化后一定要檢查離散閉環(huán)極點是否都在單位圓內(nèi)并且離單位圓周有足夠余量eig_d eig(sys_d.A - sys_d.B*K); disp(max(abs(eig_d))); % 必須小于1且小于0.95更穩(wěn)妥如果最大模在0.995說明離散閉環(huán)已經(jīng)很接近失穩(wěn)邊界無論連續(xù)域設(shè)計多完美上數(shù)字系統(tǒng)都會抖。這一步是離散域仿真的核心驗收線。5.5 現(xiàn)象仿真曲線很漂亮答辯被問“控制量峰值是多少”直接卡住這條不屬于軟件報錯而是設(shè)計與工程脫節(jié)。大作業(yè)評分常見追問就是“執(zhí)行器限幅/-10V你的控制量達(dá)到多少”。很多人的K設(shè)計只盯響應(yīng)曲線u-Kx的結(jié)果一眼沒看結(jié)果u峰值超過100一看就是不可實現(xiàn)的。解決方法是固定一個飽和模塊把輸出限幅接到被控對象前再做仿真然后觀察閉環(huán)是否還能穩(wěn)住。% 控制峰值快速計算 x0 [1; 0; 0]; t 0:0.001:5; u_max 0; for i 1:length(t) ut -K * (expm(Acl*t(i)) * x0); % 假設(shè)無觀測器動態(tài) u_max max(u_max, abs(ut)); end disp([控制峰值約: , num2str(u_max)]);這個循環(huán)本質(zhì)是在無噪聲、初值1的情況下掃一遍控制量軌跡。若超限第一步想的不應(yīng)該是增大限幅而是收回極點半徑或調(diào)小Q中狀態(tài)權(quán)重。反饋設(shè)計除了數(shù)學(xué)穩(wěn)定還要在物理約束內(nèi)穩(wěn)定否則答辯一問就翻車。下面這個速查表可以打印出來貼屏幕邊上現(xiàn)象首要懷疑優(yōu)先排查方式rank顯示缺秩數(shù)值容差svd看奇異值比值acker崩潰極點過近換place或拉開極點間距Simulink發(fā)散初值或K符號先跑step(A-B*K)隔離問題離散后抖動Ts或觀測器極點查看離散極點最大?？刂屏砍迾O點過遠(yuǎn)計算u峰值并減極點半徑6. 驗證與擴展從跑通到能答辯的一次參數(shù)掃描跑通仿真只是拿到了及格線。大作業(yè)答辯時老師真正想聽的是你的設(shè)計在參數(shù)變化下是否依然成立為什么這組極點合適。這就需要一個可復(fù)現(xiàn)的驗證腳本我的習(xí)慣是做一次帶增益掃描和極點掃描的批處理最后生成一張包絡(luò)圖。增益裕度檢驗是第一關(guān)把反饋增益乘以0.5倍和1.5倍觀察閉環(huán)特征值和響應(yīng)曲線是否仍穩(wěn)定。這個操作等價于模擬系統(tǒng)模型誤差15%甚至50%時控制器的強壯性。% 增益裕度掃描: 驗證魯棒性 alphas 0.5:0.1:1.5; % 增益縮放系數(shù) max_real zeros(size(alphas)); for i 1:length(alphas) K_test alphas(i) * K; e eig(A - B*K_test); max_real(i) max(real(e)); end plot(alphas, max_real); grid on; xlabel(增益縮放系數(shù)); ylabel(閉環(huán)最大實部); title(閉環(huán)極點隨增益變化);當(dāng)縮放系數(shù)從0.5升到1.5時max_real始終小于0說明系統(tǒng)至少具備±50%的增益裕度。若某一段越過0線說明設(shè)計剛好壓著穩(wěn)定邊界答辯時可能會被追問“模型不準(zhǔn)確怎么辦”。這個腳本是整個設(shè)計收口前必跑一遍的標(biāo)準(zhǔn)動作。參數(shù)掃描進(jìn)階玩法是對阻尼比和自然頻率做二維網(wǎng)格搜索。對每組參數(shù)算K、跑一次sim將超調(diào)量、調(diào)節(jié)時間、控制峰值記錄成表格最終選出滿足所有約束的最優(yōu)折中。% 極點參數(shù)掃描阻尼比 zeta 與自然頻率 omega zeta_range [0.5 0.707 0.9]; omega_range [2 4 8]; result []; for zeta zeta_range for omega omega_range p1 -zeta*omega 1i*omega*sqrt(1-zeta^2); p [p1; conj(p1); -5]; K_t acker(A, B, p); % 直接解析計算第2個狀態(tài)的零輸入響應(yīng)過沖 Acl_t A - B*K_t; % 用step獲取單位階躍響應(yīng)信息 y step(ss(Acl_t, B, C, D), 0:0.01:5); overshoot (max(y) - y(end)) / y(end) * 100; settle_idx find(abs(y - y(end)) 0.02*y(end), 1, last); settle_time 0.01 * (length(y) - settle_idx 1); result [result; zeta, omega, overshoot, settle_time, max(abs(K_t))]; end end % 按超調(diào)量排序打印候選設(shè)計 disp(array2table(result, ... VariableNames, {zeta, omega, overshoot, settle, Kmax}));說明一下循環(huán)邏輯外循環(huán)換阻尼比內(nèi)循環(huán)換自然頻率每組重新用acker解一次增益再基于閉環(huán)狀態(tài)空間算階躍響應(yīng)指標(biāo)。結(jié)果矩陣每一行是一組候選設(shè)計最后一列是K的最大絕對值用于評估控制量負(fù)擔(dān)。我通常把超調(diào)小于10%、調(diào)節(jié)時間小于3s、Kmax小于30的行勾出來那些行就是答辯PPT上的主推方案。有一次我趕進(jìn)度把增益掃描跑了但沒仔細(xì)看結(jié)果閉環(huán)保有極點剛好壓在-0.05虛線一畫看似穩(wěn)定結(jié)果模型參數(shù)改動10%后直接翻到右半平面。從那以后每次交這類作業(yè)我強制自己把“增益掃描極點掃描控制量峰值”三張圖全部生成并檢查一遍確保所有曲線在最惡劣參數(shù)下仍然穩(wěn)定。這套操作也能順帶回答答辯老師最常問的一個問題“你沒有沒做過魯棒性分析”希望幫到你。本文還有配套的精品資源點擊獲取