簡介本資源是一份面向自動化、控制工程專業(yè)學(xué)生及初階工程師的Python控制系統(tǒng)實踐資料包聚焦利用Python科學(xué)計算生態(tài)開展控制系統(tǒng)建模、仿真與分析。資源涵蓋SISO/MIMO系統(tǒng)建模、根軌跡與頻率響應(yīng)分析、PID控制器設(shè)計等核心內(nèi)容有效 bridging 控制理論與工程實現(xiàn)。壓縮包共135個文件以64個Python源碼含實驗?zāi)_本與算法實現(xiàn)和65張結(jié)果圖如系統(tǒng)響應(yīng)曲線、根軌跡圖、Bode圖等為主輔以XML配置、項目元數(shù)據(jù)文件整體僅1.64MB輕量易用且結(jié)構(gòu)清晰便于按實驗?zāi)K快速定位代碼與可視化結(jié)果。已有1164人學(xué)習(xí)下載適合課堂實驗延伸、課程設(shè)計參考或自學(xué)鞏固。讀者可直接運行代碼復(fù)現(xiàn)經(jīng)典控制分析流程獲取完整可執(zhí)行的仿真實驗鏈路、帶注釋的算法實現(xiàn)及對應(yīng)圖形化驗證結(jié)果顯著降低從理論到代碼落地的學(xué)習(xí)門檻。1. 這不是 Python 教程是控制系統(tǒng)工程師的「可執(zhí)行論文」用python-control庫從傳遞函數(shù)建模、根軌跡繪制、PID 整定到閉環(huán)仿真全鏈路跑通你手頭有一份控制系統(tǒng)課程設(shè)計任務(wù)書要求畫出某電機驅(qū)動系統(tǒng)的根軌跡、設(shè)計超調(diào)小于15%的 PID 控制器、給出階躍響應(yīng)曲線并標(biāo)注調(diào)節(jié)時間——但 MATLAB 許可證沒續(xù)費Simulink 模型打不開而你剛學(xué)完《自動控制原理》卻卡在“怎么把課本公式變成能跑的代碼”這一步。別急這不是理論推導(dǎo)題而是一份帶源碼、帶注釋、帶驗證邏輯的 Python 控制系統(tǒng)實戰(zhàn)論文它不講 Laplace 變換定義只做一件事——用python-control官方庫名control把傳遞函數(shù)G(s) 100/(s^2 5s 100)從建模、分析、設(shè)計到閉環(huán)驗證走完完整閉環(huán)。它面向的是正在趕課設(shè) deadline 的本科生、需要快速驗證控制器參數(shù)的嵌入式開發(fā)者以及想甩掉 MATLAB 依賴轉(zhuǎn)戰(zhàn) Python 工程棧的控制工程師。所有代碼塊均可直接復(fù)制粘貼運行所有圖表均附坐標(biāo)軸標(biāo)注與性能指標(biāo)計算邏輯所有參數(shù)修改點都標(biāo)出物理含義。這不是“Python 能不能做控制”的科普而是“今天下午三點前你就能交出一份帶仿真圖和數(shù)據(jù)表的 PDF 報告”的實操路徑。2. 從零構(gòu)建控制系統(tǒng)模型control.tf()與control.ss()的選型邏輯與建模實操控制系統(tǒng)建模是后續(xù)所有分析與設(shè)計的地基。python-control提供兩種核心建模方式傳遞函數(shù)TF與狀態(tài)空間SS。選哪個不是看誰更“高級”而是看你的原始問題來源——若題目給的是微分方程或框圖化簡結(jié)果如Y(s)/U(s) (2s1)/(s^34s^25s2)用control.tf()最直覺若你已有 A/B/C/D 矩陣比如從 Simulink 線性化導(dǎo)出、或從機械臂動力學(xué)推導(dǎo)而來則control.ss()是唯一選擇。二者本質(zhì)等價但轉(zhuǎn)換過程會暴露數(shù)值精度陷阱這點我們放在第 4 章深挖。現(xiàn)在先聚焦建模本身。2.1 用control.tf()構(gòu)建經(jīng)典二階系統(tǒng)并驗證零極點假設(shè)課程設(shè)計中某溫度控制系統(tǒng)被簡化為標(biāo)準二階模型G(s) ω_n2 / (s2 2ζω_n s ω_n2)其中自然頻率ω_n 10 rad/s阻尼比ζ 0.4。我們需要生成該傳遞函數(shù)對象并立刻驗證其零極點是否符合預(yù)期這是防止手敲系數(shù)出錯的第一道防線import control as ct import numpy as np # 定義參數(shù) wn 10.0 zeta 0.4 # 構(gòu)建傳遞函數(shù)分子為 [wn2]分母為 [1, 2*zeta*wn, wn2] num [wn**2] den [1, 2*zeta*wn, wn**2] G ct.tf(num, den) # 打印零極點驗證建模正確性 print(零點, G.zeros()) print(極點, G.poles())邏輯說明ct.tf()接收兩個列表num是分子多項式系數(shù)按 s 的降冪排列den是分母系數(shù)。此處num[100]對應(yīng)100den[1, 8, 100]對應(yīng)s2 8s 100。G.zeros()返回空數(shù)組因分子為常數(shù)無零點G.poles()應(yīng)返回共軛復(fù)數(shù)對-4 ± j9.165其模長即ω_n10實部絕對值除以ω_n即ζ0.4。參數(shù)說明wn和zeta是物理可調(diào)參數(shù)后續(xù)整定 PID 時將直接修改它們觀察響應(yīng)變化而非硬編碼系數(shù)。2.2 從狀態(tài)空間矩陣構(gòu)建同一系統(tǒng)并對比等效性現(xiàn)在我們用狀態(tài)空間形式重寫同一系統(tǒng)。對標(biāo)準二階系統(tǒng)常用可控規(guī)范型Controllable Canonical FormA [[0, 1], [-ω_n2, -2ζω_n]], B [[0], [1]], C [[ω_n2, 0]], D 0構(gòu)建后必須驗證其與 TF 模型完全等價——這是避免“以為建對了實則矩陣符號搞反”的關(guān)鍵步驟# 構(gòu)建狀態(tài)空間矩陣 A np.array([[0, 1], [-wn**2, -2*zeta*wn]]) B np.array([[0], [1]]) C np.array([[wn**2, 0]]) D 0 G_ss ct.ss(A, B, C, D) # 驗證 TF 與 SS 模型是否等價比較脈沖響應(yīng) t np.linspace(0, 2, 1000) t_tf, y_tf ct.impulse_response(G, Tt) t_ss, y_ss ct.impulse_response(G_ss, Tt) # 檢查響應(yīng)最大誤差應(yīng)遠小于 1e-10 max_error np.max(np.abs(y_tf - y_ss)) print(fTF 與 SS 脈沖響應(yīng)最大誤差{max_error:.2e}) assert max_error 1e-10, TF 與 SS 模型不等價檢查矩陣構(gòu)造邏輯說明ct.ss()接收四元組(A, B, C, D)順序不可顛倒。impulse_response()是最嚴苛的等價性檢驗——因為脈沖響應(yīng)完全由系統(tǒng)極點決定若兩者響應(yīng)一致則零極點必然相同。參數(shù)說明Tt指定仿真時間向量t必須是等間距一維數(shù)組np.linspace保證否則impulse_response會報錯或返回錯誤結(jié)果。此處t覆蓋 0~2 秒步長約 2ms足夠捕捉二階系統(tǒng)動態(tài)。2.3 處理實際工程中的非標(biāo)準模型含純延遲與多輸入多輸出MIMO真實系統(tǒng)常含純延遲e^(-sT)或多個輸入輸出通道。python-control原生不支持符號延遲但提供pade()函數(shù)進行帕德近似。例如某液位控制系統(tǒng)有 0.5 秒傳輸延遲# 添加 0.5 秒延遲3階帕德近似 delay_num, delay_den ct.pade(0.5, 3) # 返回延遲近似的分子分母系數(shù) G_delayed ct.series(G, ct.tf(delay_num, delay_den)) # 串聯(lián)延遲 # 驗證延遲效果比較階躍響應(yīng)起始時間 t_step, y_step ct.step_response(G, Tnp.linspace(0, 3, 1000)) t_delayed, y_delayed ct.step_response(G_delayed, Tnp.linspace(0, 3, 1000)) # 檢查延遲后響應(yīng)起始時間是否≈0.5秒取y0.01的第一個時刻 start_idx np.argmax(y_delayed 0.01) delay_estimated t_delayed[start_idx] print(f估計延遲時間{delay_estimated:.3f} 秒目標(biāo) 0.5 秒)邏輯說明ct.pade(T, n)生成n階帕德近似傳遞函數(shù)T為延遲時間。ct.series()實現(xiàn)模塊串聯(lián)。參數(shù)說明n3是工程常用折中——階數(shù)太低如 n1近似誤差大太高如 n10引入虛假高頻極點導(dǎo)致仿真不穩(wěn)定。此處n3在 0.5 秒延遲下誤差可控。對于 MIMO 系統(tǒng)ct.tf()支持二維分子分母列表如G_mimo ct.tf([[[1], [2]], [[3], [4]]], [[[1,1], [1,2]], [[1,3], [1,4]]])表示 2x2 傳遞函數(shù)矩陣但需確保每個子項維度匹配。3. 根軌跡與頻域分析root_locus()與bode_plot()的參數(shù)精調(diào)與物理意義解讀建模完成后必須回答兩個核心問題系統(tǒng)穩(wěn)定嗎如何通過調(diào)整參數(shù)如 PID 增益使其滿足性能指標(biāo)根軌跡Root Locus和伯德圖Bode Plot是控制工程師的“X 光機”——它們不告訴你具體怎么調(diào)但清晰顯示“往哪調(diào)、調(diào)多少、會變什么樣”。python-control的繪圖函數(shù)默認設(shè)置常導(dǎo)致關(guān)鍵信息丟失必須手動干預(yù)參數(shù)。3.1 繪制高保真根軌跡規(guī)避默認范圍失真與分支斷裂對前述二階系統(tǒng)G(s) 100/(s2 5s 100)添加比例控制器K后開環(huán)傳函為K*G(s)。根軌跡顯示所有可能的閉環(huán)極點位置。但默認root_locus()會自動縮放坐標(biāo)軸導(dǎo)致關(guān)鍵區(qū)域如虛軸附近被壓縮甚至遺漏分支import matplotlib.pyplot as plt # 創(chuàng)建開環(huán)系統(tǒng) K*G(s) G_open ct.tf([100], [1, 5, 100]) # 關(guān)鍵手動指定 K 的范圍和點數(shù)避免默認采樣過粗 K_range np.logspace(-2, 3, 1000) # K 從 0.01 到 1000對數(shù)均勻采樣 rlist, klist ct.root_locus(G_open, kvectK_range, plotFalse) # 手動繪制強制坐標(biāo)軸范圍覆蓋關(guān)鍵區(qū)域 plt.figure(figsize(8, 6)) plt.plot(rlist.real, rlist.imag, b-, linewidth1.2, label根軌跡) plt.axhline(y0, colork, linewidth0.8, linestyle--) # 實軸 plt.axvline(x0, colork, linewidth0.8, linestyle--) # 虛軸 plt.xlim(-15, 5) # 強制 x 范圍覆蓋左半平面穩(wěn)定區(qū)及右半平面不穩(wěn)定區(qū) plt.ylim(-15, 15) # 強制 y 范圍確保共軛極點對可見 plt.xlabel(實部) plt.ylabel(虛部) plt.title(根軌跡圖K 從 0.01 到 1000) plt.grid(True, alpha0.3) plt.legend() plt.show()邏輯說明kvectK_range顯式指定增益向量plotFalse禁用默認繪圖以便手動控制。plt.xlim/ylim是核心——默認范圍常以極點為中心縮放而我們關(guān)心的是極點是否穿越虛軸穩(wěn)定性邊界。參數(shù)說明K_range用np.logspace而非np.linspace因為根軌跡在小 K 和大 K 區(qū)域變化劇烈程度不同對數(shù)采樣能更好捕捉臨界點如分離點、與虛軸交點。此處(-2,3)覆蓋典型工程增益范圍。3.2 伯德圖增益裕度與相位裕度的自動化提取伯德圖用于判斷系統(tǒng)魯棒性。python-control的margin()函數(shù)可直接計算增益裕度GM和相位裕度PM但需注意其輸入必須是開環(huán)系統(tǒng)且margin()返回的頻率單位是 rad/s非 Hz易被忽略# 計算開環(huán)系統(tǒng) G_open 的頻域裕度 gm, pm, wg, wp ct.margin(G_open) print(f增益裕度 GM {gm:.2f} (線性值)即 {20*np.log10(gm):.1f} dB) print(f相位裕度 PM {pm:.1f}°) print(f增益穿越頻率 wg {wg:.2f} rad/s ({wg/(2*np.pi):.2f} Hz)) print(f相位穿越頻率 wp {wp:.2f} rad/s ({wp/(2*np.pi):.2f} Hz)) # 繪制伯德圖并標(biāo)注裕度點 plt.figure(figsize(10, 6)) ct.bode_plot(G_open, omega_limits(0.1, 100), dBTrue, degTrue) plt.show()邏輯說明ct.margin()返回四個值gm增益裕度線性值、pm相位裕度度數(shù)、wg增益穿越頻率 rad/s、wp相位穿越頻率 rad/s。20*np.log10(gm)將線性值轉(zhuǎn)為工程常用的 dB 值。參數(shù)說明omega_limits(0.1, 100)強制伯德圖橫軸范圍避免默認范圍遺漏低頻或高頻段dBTrue和degTrue確?？v軸為 dB 和度符合行業(yè)習(xí)慣。注意wg和wp是 rad/s若需 Hz必須除以2π。3.3 尼科爾斯圖Nichols Chart與閉環(huán)帶寬的關(guān)聯(lián)分析當(dāng)需要評估系統(tǒng)對參考信號的跟蹤能力如伺服系統(tǒng)時閉環(huán)帶寬Bandwidth比開環(huán)穿越頻率更具物理意義。尼科爾斯圖可直觀定位帶寬點閉環(huán)幅值下降 3dB 處# 計算閉環(huán)系統(tǒng) T(s) G(s)/(1G(s)) 的頻率響應(yīng) T_closed ct.feedback(G_open, 1) # 單位負反饋 # 獲取閉環(huán)頻率響應(yīng) omega np.logspace(-1, 2, 1000) # 0.1 到 100 rad/s mag, phase, _ ct.freqresp(T_closed, omega) # 找到閉環(huán)帶寬|T(jω)| 1/√2 ≈ 0.707 (即 -3dB) mag_db 20 * np.log10(mag) bandwidth_idx np.argmax(mag_db -3.0) # 第一個低于 -3dB 的點 bandwidth_omega omega[bandwidth_idx] print(f閉環(huán)帶寬 {bandwidth_omega:.2f} rad/s ({bandwidth_omega/(2*np.pi):.2f} Hz)) # 繪制尼科爾斯圖需安裝 control 庫 0.9.0 try: ct.nichols_plot(G_open, gridTrue) plt.title(尼科爾斯圖開環(huán)頻率響應(yīng)) plt.show() except AttributeError: print(當(dāng)前 control 版本不支持 nichols_plot使用 bode 替代)邏輯說明ct.feedback(G, 1)默認實現(xiàn)單位負反饋T G/(1G)。freqresp()返回閉環(huán)系統(tǒng)在指定omega下的幅值mag和相位phase。帶寬定義為閉環(huán)幅值首次降至 -3dB 的頻率。參數(shù)說明omega必須用np.logspace覆蓋足夠?qū)挼念l段否則可能找不到帶寬點np.argmax(mag_db -3.0)利用布爾索引找到第一個滿足條件的索引比np.interp插值更魯棒。4. PID 控制器設(shè)計與整定control.pid()的底層實現(xiàn)與手動整定避坑指南python-control不提供黑箱 PID 自動整定函數(shù)如 MATLAB 的pidtune但這恰是優(yōu)勢——它迫使你理解 PID 參數(shù)的物理作用。我們用control.pid()創(chuàng)建控制器對象再通過series()與被控對象串聯(lián)最后用feedback()構(gòu)成閉環(huán)。但直接套用教科書公式如 Ziegler-Nichols常翻車原因在于模型失配與離散化效應(yīng)。4.1 用control.pid()構(gòu)建標(biāo)準 PID 并驗證離散化一致性control.pid()默認生成連續(xù)時間 PID但實際部署常需離散化如單片機。python-control的c2d()函數(shù)支持多種離散化方法其中 Tustin雙線性變換在中頻段精度最高# 創(chuàng)建連續(xù) PID 控制器Kp10, Ki2, Kd1 C_cont ct.pid(Kp10, Ki2, Kd1, Tf0.1) # Tf 為微分濾波時間常數(shù) # 離散化采樣周期 Ts0.01sTustin 方法 Ts 0.01 C_disc ct.c2d(C_cont, Ts, methodtustin) # 驗證離散化后 DC 增益是否一致關(guān)鍵 s_cont ct.tf([1, 0], [1]) # s operator dc_gain_cont ct.dcgain(C_cont) # 連續(xù) PID 的 DC 增益Ki/0 → ∞故為 inf # 離散 PID 的 DC 增益應(yīng)為 sum of coefficients dc_gain_disc C_disc.dcgain() print(f離散 PID DC 增益 {dc_gain_disc:.2f}) # 更實用的驗證比較階躍響應(yīng)連續(xù) vs 離散 t_sim np.linspace(0, 2, 2000) t_cont, y_cont ct.step_response(C_cont, Tt_sim) t_disc, y_disc ct.step_response(C_disc, Tt_sim) # 計算響應(yīng)均方誤差MSE mse np.mean((y_cont - y_disc)**2) print(f連續(xù)與離散 PID 階躍響應(yīng) MSE {mse:.2e})邏輯說明ct.pid()的Tf參數(shù)至關(guān)重要——它限制微分作用的高頻增益防止噪聲放大。ct.c2d()的methodtustin是首選zoh零階保持在低頻好但中高頻相位滯后大。參數(shù)說明Ts0.01對應(yīng) 100Hz 采樣率是常見嵌入式系統(tǒng)上限Tf0.1意味著微分作用在 10 rad/s 以上被衰減平衡響應(yīng)速度與抗噪性。4.2 手動整定 PID基于根軌跡的“三步法”與性能指標(biāo)映射教科書 Ziegler-Nichols 法要求先找到臨界比例度但實際系統(tǒng)可能不穩(wěn)定。更穩(wěn)健的方法是1) 用根軌跡確定使閉環(huán)極點位于期望阻尼比 ζ0.70745°線和自然頻率 ω_n15 rad/s 的 Kp2) 加入積分消除穩(wěn)態(tài)誤差3) 加入微分抑制超調(diào)。以下為完整流程# 步驟1僅用比例控制根軌跡上找 Kp 使主導(dǎo)極點滿足 ζ0.707, ω_n15 G_open_p ct.tf([100], [1, 5, 100]) rlist_p, klist_p ct.root_locus(G_open_p, kvectnp.logspace(-1, 2, 500), plotFalse) # 尋找最接近期望極點 (-10.6, ±10.6) 的實際極點ζ0.707 σω_n/√210.6 target_real -10.6 target_imag 10.6 distances np.sqrt((rlist_p.real - target_real)**2 (rlist_p.imag - target_imag)**2) best_idx np.argmin(distances) Kp_desired klist_p[best_idx] print(f步驟1Kp {Kp_desired:.2f} 使主導(dǎo)極點接近 ζ0.707, ω_n15) # 步驟2加入積分控制PI設(shè)計 Ki 使穩(wěn)態(tài)誤差為0且不顯著惡化動態(tài) C_pi ct.tf([Kp_desired, Kp_desired*0.5], [1, 0]) # Ki Kp*0.5經(jīng)驗初值 G_open_pi ct.series(C_pi, G_open_p) t_pi, y_pi ct.step_response(ct.feedback(G_open_pi, 1), Tnp.linspace(0, 2, 1000)) ess_pi 1 - y_pi[-1] # 穩(wěn)態(tài)誤差 print(f步驟2PI 控制后穩(wěn)態(tài)誤差 {ess_pi:.4f}) # 步驟3加入微分PIDKd0.1*Kp 抑制超調(diào) C_pid ct.tf([Kp_desired*0.1, Kp_desired, Kp_desired*0.5], [1, 0.1]) # Kd0.1*Kp, Tf0.1 G_open_pid ct.series(C_pid, G_open_p) t_pid, y_pid ct.step_response(ct.feedback(G_open_pid, 1), Tnp.linspace(0, 2, 1000)) # 計算性能指標(biāo) overshoot (np.max(y_pid) - 1) * 100 ts_2percent t_pid[np.argmax(y_pid 0.98)] # 2% 調(diào)節(jié)時間 print(f步驟3PID 后超調(diào) {overshoot:.1f}%, 調(diào)節(jié)時間 {ts_2percent:.3f}s)邏輯說明此“三步法”將抽象指標(biāo)ζ, ω_n映射到具體參數(shù)。Kp由根軌跡確定Ki由穩(wěn)態(tài)誤差需求反推Kd由超調(diào)抑制需求設(shè)定。參數(shù)說明Ki Kp*0.5是經(jīng)驗起點若穩(wěn)態(tài)誤差仍大可增大 KiKd 0.1*Kp是保守初值若超調(diào)仍高可增至0.2*Kp但需警惕噪聲放大。4.3 避坑PID 整定中五個血淚教訓(xùn)與排查方案現(xiàn)象 → 原因 → 解決現(xiàn)象加入 PID 后系統(tǒng)振蕩加劇甚至發(fā)散。原因Kd過大或Tf過小導(dǎo)致微分作用在高頻段增益過高放大測量噪聲并激發(fā)未建模模態(tài)。解決立即減小Kd至原值 1/3同時增大Tf至1/(10*Kd)量級用bode_plot(C_pid)檢查微分環(huán)節(jié)在 100 rad/s 以上是否仍具高增益?，F(xiàn)象階躍響應(yīng)中存在緩慢爬升的“尾巴”非單調(diào)收斂。原因積分時間常數(shù)Ti Kp/Ki過大導(dǎo)致積分作用太弱無法及時消除穩(wěn)態(tài)誤差。解決將Ki增大 2~3 倍重新仿真若超調(diào)隨之增大則同步微調(diào)Kp降低至原值 0.8 倍?，F(xiàn)象離散 PID 在采樣率Ts0.02s下響應(yīng)正常但Ts0.005s時出現(xiàn)高頻抖動。原因c2d()的 Tustin 方法在ω_s 10*ω_cω_c為系統(tǒng)帶寬時數(shù)值不穩(wěn)定Ts過小導(dǎo)致離散化后極點靠近單位圓外側(cè)。解決改用zoh方法重離散化或?qū)_cont先加一階低通濾波ct.tf([1], [0.01, 1])再離散化?，F(xiàn)象ct.feedback(G, C)與ct.feedback(C, G)結(jié)果不同且后者導(dǎo)致仿真崩潰。原因ct.feedback()默認sign-1負反饋但若控制器C本身含負號如某些狀態(tài)反饋結(jié)構(gòu)誤用ct.feedback(G, C)會形成正反饋。解決始終顯式指定ct.feedback(G, C, sign-1)用ct.minreal()對閉環(huán)系統(tǒng)降階檢查是否有極點在右半平面?，F(xiàn)象root_locus()繪制的軌跡在 K 很大時出現(xiàn)“斷點”或虛線。原因kvect采樣點過少高增益區(qū)極點移動劇烈線性插值失效。解決將kvect改為np.logspace(-1, 4, 2000)并用plt.plot(rlist.real, rlist.imag, b., markersize0.5)以散點圖替代連線避免誤導(dǎo)。5. 閉環(huán)系統(tǒng)仿真與性能驗證forced_response()與自定義性能指標(biāo)計算設(shè)計完控制器必須用真實激勵非理想階躍驗證其魯棒性。python-control的forced_response()支持任意輸入信號而不僅是step_response()或impulse_response()。更重要的是性能指標(biāo)超調(diào)、調(diào)節(jié)時間、穩(wěn)態(tài)誤差必須用代碼精確計算而非目測圖表——這是工程報告可信度的底線。5.1 用forced_response()仿真真實工況含噪聲與擾動的輸入實際系統(tǒng)面臨傳感器噪聲與負載擾動。我們用forced_response()注入白噪聲與階躍擾動觀察閉環(huán)系統(tǒng)抑制能力# 構(gòu)建閉環(huán)系統(tǒng)PID 控制下的被控對象 C_final ct.tf([1.0, 10.0, 5.0], [1, 0.1]) # Kp10, Ki5, Kd1, Tf0.1 G_cl ct.feedback(ct.series(C_final, G_open_p), 1) # 生成時間向量10000 點覆蓋 0~10 秒 t_sim np.linspace(0, 10, 10000) dt t_sim[1] - t_sim[0] # 構(gòu)建復(fù)合輸入?yún)⒖夹盘?~2s 階躍2~4s 斜坡4~10s 正弦 擾動4s 時加入 0.2 幅值階躍 u_ref np.zeros_like(t_sim) u_ref[(t_sim 0) (t_sim 2)] 1.0 u_ref[(t_sim 2) (t_sim 4)] (t_sim[(t_sim 2) (t_sim 4)] - 2) * 0.5 1.0 u_ref[t_sim 4] 1.0 0.5 * np.sin(2*np.pi*0.5*(t_sim[t_sim 4] - 4)) u_dist np.zeros_like(t_sim) u_dist[t_sim 4] 0.2 # 4s 時加入 0.2 幅值擾動 # 總輸入 參考 擾動注意擾動通常加在被控對象輸入端故需重構(gòu)系統(tǒng) # 此處簡化將擾動作為額外輸入用 superposition 原理疊加響應(yīng) _, y_ref ct.forced_response(G_cl, Tt_sim, Uu_ref) _, y_dist ct.forced_response(G_cl, Tt_sim, Uu_dist) # 注意此用法需 G_cl 為 SISO且擾動通道已建模 y_total y_ref y_dist # 線性系統(tǒng)疊加原理 # 添加傳感器白噪聲SNR40dB noise_power np.var(y_total) / (10**(40/10)) y_noisy y_total np.random.normal(0, np.sqrt(noise_power), len(y_total)) # 繪制結(jié)果 plt.figure(figsize(12, 8)) plt.subplot(2,1,1) plt.plot(t_sim, u_ref, g-, label參考輸入, linewidth1.2) plt.plot(t_sim, y_noisy, b-, label含噪輸出, linewidth1.2) plt.grid(True, alpha0.3) plt.legend() plt.ylabel(幅值) plt.title(復(fù)合工況仿真參考信號 擾動 噪聲) plt.subplot(2,1,2) plt.plot(t_sim, y_noisy - u_ref, r-, label跟蹤誤差, linewidth1.2) plt.axhline(y0, colork, linestyle--, alpha0.5) plt.grid(True, alpha0.3) plt.xlabel(時間 (s)) plt.ylabel(誤差) plt.legend() plt.show()邏輯說明forced_response()的U參數(shù)接受一維數(shù)組長度必須與T相同。superposition原理在此成立因系統(tǒng)線性。參數(shù)說明dt由t_sim間隔隱含影響數(shù)值積分精度SNR40dB是典型工業(yè)傳感器信噪比noise_power計算確保噪聲功率與信號功率比為10^(-40/10)0.01。5.2 自定義性能指標(biāo)函數(shù)超調(diào)、調(diào)節(jié)時間、穩(wěn)態(tài)誤差的魯棒計算目測圖表誤差大必須用代碼量化。以下函數(shù)處理常見異常響應(yīng)未達穩(wěn)態(tài)、存在振蕩、噪聲干擾def compute_performance_metrics(t, y, ref_funcNone, threshold0.02): 計算閉環(huán)系統(tǒng)性能指標(biāo) :param t: 時間向量 :param y: 輸出響應(yīng)向量 :param ref_func: 參考信號函數(shù)若為 None 則假設(shè)為單位階躍 :param threshold: 調(diào)節(jié)時間閾值如 0.02 表示 2% :return: dict 包含 overshoot, ts, ess, rise_time if ref_func is None: y_ref np.ones_like(y) else: y_ref ref_func(t) error y - y_ref # 超調(diào)僅對階躍響應(yīng)定義 if ref_func is None or (callable(ref_func) and ref_func(0)ref_func(1)): overshoot max(0, (np.max(y) - 1) * 100) else: overshoot 0.0 # 上升時間10% 到 90% y_min, y_max np.min(y_ref), np.max(y_ref) y10 y_min 0.1*(y_max - y_min) y90 y_min 0.9*(y_max - y_min) idx_10 np.argmax(y y10) idx_90 np.argmax(y y90) rise_time t[idx_90] - t[idx_10] if idx_90 idx_10 else np.nan # 調(diào)節(jié)時間進入 ±threshold 帶內(nèi)不再出來 band_upper y_ref threshold * np.abs(y_ref) band_lower y_ref - threshold * np.abs(y_ref) in_band (y band_upper) (y band_lower) # 找到最后一個連續(xù) in_band 的起始點 last_out_idx np.where(~in_band)[0] if len(last_out_idx) 0: ts t[-1] else: last_out last_out_idx[-1] ts t[np.argmax(in_band[last_out:]) last_out] if np.any(in_band[last_out:]) else np.nan # 穩(wěn)態(tài)誤差最后 10% 時間段的平均誤差 tail_len int(0.1 * len(y)) ess np.mean(np.abs(error[-tail_len:])) return { overshoot_percent: round(overshoot, 2), rise_time_sec: round(rise_time, 3), settling_time_sec: round(ts, 3), steady_state_error: round(ess, 4) } # 使用示例對前述復(fù)合仿真結(jié)果計算指標(biāo) metrics compute_performance_metrics(t_sim, y_noisy, ref_funclambda t: np.piecewise(t, [t2, (t2)(t4), t4], [1, lambda t: (t-2)*0.51, lambda t: 10.5*np.sin(2*np.pi*0.5*(t-4))]), threshold0.02) print(性能指標(biāo), metrics)邏輯說明該函數(shù)處理三種典型場景純階躍、斜坡、正弦參考。threshold0.02對應(yīng) 2% 調(diào)節(jié)時間標(biāo)準。參數(shù)說明ref_func用np.piecewise構(gòu)建分段函數(shù)compute_performance_metrics內(nèi)部用np.argmax和布爾索引確保計算魯棒即使響應(yīng)未完全收斂也返回合理值如tsnp.nan。5.3 生成符合課程設(shè)計要求的 PDF 報告Matplotlib 導(dǎo)出與 LaTeX 公式嵌入最終交付物常為本文還有配套的精品資源點擊獲取