如何用 python 代碼復現 MATLAB simulink 的 PID

MATLAB
- 在 Simulink 里做以下設置
- MATLAB 腳本調用示例
python 實現
- 離散 PID 實現（并行形式）
Simulink 中兩種 PID 結構（并聯形式, I-形式）下連續/離散時域里積分增益 I 的表示
- 并聯（Parallel） vs 理想（I-Form），use I*Ts 的勾選
- 什么時候用連續域 vs 離散域？

MATLAB

首先，MATLAB simulink 的 PID 長這樣：
在這里插入圖片描述

其中，PID 控制器以 1e-4 的步長運行，而物理模型用固定步長 1e-5 的求解器來集成。

這樣，每仿真一次 StopTime = 1e-4，就得到了“10 步模型響應 + 1 次 PID 更新”之后的 $Z$ 和 $V\_IC1$ ，非常符合“每個控制電壓需要響應 10 次”的需求。

在 Simulink 里做以下設置

Solver
- 打開 Model Settings (Ctrl+E) → Solver
- Type 選 Fixed-step
- Solver 選 ode4 (Runge–Kutta)（或者其他固定步長求解器）
- Fixed-step size 填 1e-5
PID Controller
- 雙擊 PID 塊（或 PID Controller），在 Discrete 模式下
- Sample time (Ts) 填 1e-4
Rate Transition（可選，但推薦）
- 在 PID 輸出和物理模型輸入之間加一個 Rate Transition 塊，保證兩條不同速率信號安全切換。
- 默 1e-4 → 1e-5，不用改參數，Simulink 會自動用最近鄰保持。
Initial Condition
- 把物理模型里初始化狀態的那一支“Initial Condition”塊的初值寫成變量 x0（工作區定義）。
To Workspace
- 用 To Workspace 塊把 outZ（物理模型的 Z）和 outV_IC1（PID 輸出）導出。

MATLAB 腳本調用示例

%%—— 預先在 base workspace 定義好初始狀態和參考值 ——%%
x0 = [...];             % 初始狀態向量
Z_reference = 0;        % 參考值恒為 0assignin('base','x0',x0);
assignin('base','Z_reference',Z_reference);%%—— 設置仿真——%%
Ts_ctrl   = 1e-4;       % 控制器采樣周期
Ts_solver = 1e-5;       % 求解器步長mdl = 'PID';
open_system(mdl);% 1) 全局求解器
set_param(mdl, ...'Solver',        'ode4', ...         % 固定步長 RK4'FixedStep',     num2str(Ts_solver), ...'StopTime',      num2str(Ts_ctrl)    % 仿真到一個控制步
);% 2) PID 塊采樣時間
set_param([mdl '/PID Controller'], ...'SampleTime', num2str(Ts_ctrl) ...
);%%—— 運行仿真 ——%%
simOut = sim(mdl, ...'SrcWorkspace','base', ...'SaveOutput','on', ...'ReturnWorkspaceOutputs','on');%%—— 提取最后一步 ——%%
Z_all     = simOut.get('outZ');        % time series
V_IC1_all = simOut.get('outV_IC1');Z     = Z_all(end);
V_IC1 = V_IC1_all(end);fprintf('t=%.0e 時刻： Z = %.6f,  V_IC1 = %.6f\n', Ts_ctrl, Z, V_IC1);

python 實現

對于相同的求解器初始值，如果 MATLAB PID 和 python PID 的 第一步輸出相差一個數量級，最常見的坑就是：

在這里插入圖片描述

MATLAB PID 塊的參數含義可能和你認為的不一樣（詳見下文 Simulink 中兩種 PID 結構中 I 增益的表示）

Simulink 默認的 PID 塊有好幾種形式（Parallel / I-form / P-only …），比如，連續時域中 它的 “I 增益” 在 Parallel 形式下 表示的是 $\frac{K_i}{s}$ 。如果直接把塊里的 “I 增益” 當成 Python 里的 $K_i$ ，就會差好幾倍。
檢查：雙擊 Simulink PID 塊，看它 是連續時間還是離散時間，是 Parallel 形式還是 I-form，然后 把參數轉換成 Python 里的 $K_p,\,K_i,\,K_d,\,N)$ 平行形式。

對號入座，才能在 Python 里準確還原 MATLAB/Simulink 的 PID 積分行為。

離散 PID 實現（并行形式）

Simulink 的離散 PID Controller 塊并不是把連續‐時間的

$u(t)=-\bigl(K_p e + K_i\!\int e\,dt + K_d\dot e\bigr)$

簡單地用歐拉積分＋一階濾波來近似。它用的是一整套 “ $z$ 域” 差分方程：

$C(z)\;=\;P \;+\; I\,T_s\frac{1}{z-1} \;+\; D\,\frac{N}{1+N\,T_s}\,\frac{z-1}{z}$

積分項

$u_I[k]=u_I[k-1]+I\,T_s\;e[k]$
微分項（帶濾波）
令 $N_d=N\,T_s$ ，有

$y_D[k] =\frac{1}{1+N_d}\,y_D[k-1] \;+\;\frac{K_d\,N_d}{1+N_d}\bigl(e[k]-e[k-1]\bigr)$
總輸出

$u[k] = K_p\,e[k] + u_I[k] + y_D[k]$

而如果在 Python 里用的是連續時間那套，

new_integral = integral + e*dt
derivative  = (N*dt*raw_derivative + prev_filtered)/(1+N*dt)
u =  Kp*e + Ki*integral + Kd*derivative

它們并不等價。

下面這段代碼給出了一個跟 Simulink 離散 PID Controller 塊 一模一樣的差分實現。它對應于塊上方顯示的 并行形式，

$\;+\; I\,T_s\frac{1}{z-1} \;+\; D\frac{N}{1+N\,T_s}\frac{z-1}{z},$

其中 “積分(I)”框里顯示的就是 $I\,T_s$ （注意勾選），所以在代碼里直接用它來做積分增量。積分輸出 和 最終的控制量 都被鉗位在 $[-V_{\rm lim},\,V_{\rm lim}]$ 以內。

$P$ 項沒有啥好說的，首先根據上面說的實現 $I$ 項，

# 內部狀態
self.uI     = 0.0   # 上一步積分項 u_I[k-1]
self.yD     = 0.0   # 上一步濾波微分 y_D[k-1]
self.e_prev = 0.0   # 上一次誤差 e[k-1]

# 積分項（差分形式）
uI_candidate  = self.uI_prev + self.Ki * error   # use I*Ts in simulink
# Anti-windup: limit integral term 將 integral 限在 [-V_lim, +V_lim]
uI_candidate  = torch.clamp(uI_candidate, -self.V_lim, self.V_lim)

注意，在 Simulink 離散 PID Controller 塊中，積分項的更新總是滯后一拍，它遵循“先用上一步（或初始）的狀態算輸出，再更新狀態” 的時序。

可以在 MATLAB 中只打開 I 控制器實驗，會發現第一步的控制器輸出為 0。

初始時刻 $t = 0$ 的輸出

塊會先把“初始積累值” $u_I[0]$ （默認由 Integrator Initial Condition 參數決定，缺省為 0）和“初始濾波值” $y_D[0]$ （同樣缺省為 0）帶入控制律中。
此時即便誤差 $e[0]\neq0$ ，積分項 還沒有 加上 $\,K_iT_s\cdot e[0]$ 的那部分增量，因此積分輸出項是 0；

在第一個采樣周期結束（ $t=T_s$ ）時

Simulink 才會用 $e [0]$ 來更新積分器：
$u_I[1] = u_I[0] + (K_iT_s)\,e[0],$
并在此刻把新的 $u_I[1]$ 參與下一次輸出計算。

把 forward() 實現改成這樣，就能和 Simulink 保持一致：

# 先用 P項、舊濾波、微分項 計算輸出
u_unsat = self.Kp*e + self.uI_prev + self.yD# 再更新積分狀態，為下一步準備
uI_candidate  = self.uI_prev + self.Ki * error   # use I*Ts in simulink

接下來是微分項 $D$ 的實現，

在這里插入圖片描述

 # 更新 D 項 —— 前向歐拉離散濾波, 差分方程來源于： yD[k] = yD[k-1] + Ts*( Kd*N*(e[k]-e[k-1])/Ts - N*yD[k-1] )
raw_deriv = (error - self.e_prev) / self.dt
self.yD = self.yD_prev + self.dt * ( self.Kd * self.N * raw_deriv - self.N * self.yD_prev )

Simulink 中兩種 PID 結構（并聯形式, I-形式）下連續/離散時域里積分增益 I 的表示

先把幾個概念理清：

連續域的 $s$ （拉普拉斯算子）
- 把時間域信號 $f (t)$ 變換到復頻域后，微分 $\tfrac{d}{dt}$ 對應乘以 $s$ ：
  
  $\mathcal{L}\{\,\dot x(t)\} = s X(s).$
- 因此， $\dfrac{K_i}{s}$ 就是“對誤差做積分”—— $K_i$ 是積分增益，單位大致是“控制量／（誤差×秒）”。
離散域的 $\tfrac{1}{z-1}$ （差分算子）
- 把時間序列 $x [k]$ 做 $z$ 變換后，
  
  $\mathcal{Z}\{\,x[k]-x[k-1]\} = (1 - z^{-1})X(z) \quad\Longrightarrow\quad \mathcal{Z}\Bigl\{\sum_{i=0}^k x[i]\Bigr\} = \frac{1}{1-z^{-1}}X(z) = \frac{z^{-1}}{1-z^{-1}}\,X(z)\,.$
- 換一種常見寫法 $\tfrac{1}{z-1} = -\tfrac{z^{-1}}{1-z^{-1}}$ ，它就對應“累加求和”那一塊。離散積分項一般寫成
  
  $K_i\,T_s\;\frac{1}{z-1} \quad\Longleftrightarrow\quad u_I[k] = u_I[k-1] + K_i\,T_s\,e[k].$
- 這里 $T_s$ 是采樣周期，把每步的積分增益拉成了和連續域里 $K_i/s$ 數值等價的 “每步累加” 形式。

并聯（Parallel） vs 理想（I-Form），use I*Ts 的勾選

特性	并聯（Parallel）形式	I-Form（Ideal）形式
連續時域公式	$C(s)=K_p + \dfrac{K_i}{s} + K_d\,\frac{N}{1+N \frac{1}{s}}$	$C(s)=K_p + \dfrac{K_p}{T_i}\,\dfrac{1}{s} + K_p K_d \frac{N}{1+N \frac{1}{s}}$
離散時域公式	$C(z)=K_p + K_i T_s \dfrac{1}{z-1} + K_d \frac{N}{1 + N T_s\frac{1}{z-1}}$	$C(z)=K_p + \dfrac{K_p T_s}{T_i}\,\dfrac{1}{z-1} + K_p K_d \frac{N}{1 + N T_s\frac{1}{z-1}}$
離散時域差分實現	$u_I[k]=u_I[k-1]+K_i\,T_s\,e[k]$	$\;u_I[k]=u_I[k-1]+\tfrac{K_p\,T_s}{T_i}\,e[k]$

在 離散模式 下，
- 勾選 “Use I*Ts” 時，PID Controller 塊中的“積分 (I)”文本框顯示的就是 $K_i\,T_s$ ，這個數值直接用于差分方程 $u_I[k]=u_I[k-1]+(K_iT_s)\,e[k]$ 。
- 若 不勾選 “Use I*Ts”，文本框則顯示純粹的離散積分增益 $K_i$ ，塊會在后臺自動再乘以 $T_s$ 來進行積分運算，即等效地使用 $K_iT_s$ 。
在 連續模式 下，無論是否有“Use I*Ts”，文本框顯示的都是連續域的積分增益
$\dfrac{K_i}{s}$ 中的 $K_i$ ，即單位為“控制量/（誤差·秒）”的參數。

結構類型	文本框顯示（離散）	差分方程增量	Simulink 幫你做的事
并聯（Parallel）	$K_iT_s$ （勾選 Use I*Ts） $K_i$ （不勾選）	$u_I[k]=u_I[k-1]+(K_iT_s)\,e[k]$	若不勾選，塊會自動乘以 $T_s$ ；勾選則用文本框數值直接累加。
理想（I-Form）	$T_i$ （始終）	$\;u_I[k]=u_I[k-1]+\tfrac{K_p\,T_s}{T_i}\,e[k]$	勾選后內部把 $T_i$ 換算為 $\tfrac{K_pT_s}{T_i}$ 。

根據 MATLAB 官方文檔，

在這里插入圖片描述

【Tip】

如果你已經在并聯下直接獲得了離散 $K_iT_s$ （比如 Simulink 給你的那一欄），就直接拿來用，不要再二次換算；
在并聯（Parallel）結構下，Simulink 參數欄 “積分(I)”框里顯示的正是 $K_i\,T_s$ （離散）或 $K_i$ （連續），直接把 Simulink 中讀到的 “積分(I)” 值當做 $K_i\,T_s$ （離散）或 $K_i$ （連續）使用即可。
如果是在 理想（I-Form） 下得到一個積分時間 $T_i$ ，“積分( I )” 填的是 $T_i$ ，它會自動在內部換算成 $K_pT_s/T_i$ ；python 實現需要再用
$K_{i,\mathrm{discrete}} = \frac{K_p\,T_s}{T_i}$
換算成“每步累加增量”。

什么時候用連續域 vs 離散域？

如果把整個閉環模型當成一個離散‐時間系統 —— 比如用

$x_{k+1} = e^{A\,T_s}\,x_k \;+\; A^{-1}(e^{A\,T_s}-I)\,B\,u_k$

這樣的精確離散更新公式（基于矩陣指數）來推進狀態，那么整個控制器的設計、仿真與實現，就更自然地放在離散時域里：
- 在每個采樣點 $k$ 計算誤差 $e [k]$ ，更新離散 PID 差分方程，輸出 $u_k$ ，然后用上面的公式算 $x_{k+1}$ 。
- 這時把 Simulink 里“并聯形式”的離散 PID 差分方程直接拿過來用，就保證“數值上一致”。
- 連續域的 $K_i/s$ 并不直接作用于這個仿真步長里，它只是一個設計思想——真正跑的時候都要離散化。
如果你打算在連續‐時間下設計 PID（比如先用頻域調參工具、根軌跡、Bode 曲線，得到 $K_i,K_p,K_d$ ），再把它離散化（如向后差分、雙線性變換 Tustin 法、零階保持 ZOH），那么就先在連續時域下選好參數，然后再用某種離散化方法轉換成離散差分形式。