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