手把手教你学Simulink——基于Simulink的扩展卡尔曼滤波（EKF）PMSM状态估计

手把手教你学Simulink

——基于Simulink的扩展卡尔曼滤波（EKF）PMSM状态估计

一、问题背景

在永磁同步电机（PMSM）无位置传感器控制中，需实时估计转子位置 (\theta_r) 和转速 (\omega_r)。传统龙伯格观测器依赖线性模型，在非线性、噪声环境下性能受限。

扩展卡尔曼滤波（Extended Kalman Filter, EKF）通过一阶泰勒展开处理非线性系统，并最优融合模型预测与测量值，具有：

• 强非线性处理能力
• 显式考虑过程/测量噪声
• 自适应增益调节

本教程将在 Simulink 中实现 基于 EKF 的 PMSM 全状态估计器，涵盖：

1. 非线性状态空间建模
2. EKF 算法推导与离散化
3. 噪声协方差整定
4. 与真实状态对比验证

先决条件：需理解卡尔曼滤波基础；推荐使用 MATLAB R2020a 及以上版本

二、PMSM 非线性状态空间模型（αβ 坐标系）

为避免角度依赖，采用静止 αβ 坐标系，并选择以下状态向量：

[
\mathbf{x} = \begin{bmatrix}
i_\alpha \
i_\beta \
\omega_r \
\theta_r
\end{bmatrix}
]

连续时间状态方程：

[
\begin{aligned}
\frac{di_\alpha}{dt} &= -\frac{R_s}{L} i_\alpha + \frac{1}{L} v_\alpha + \frac{\psi_f \omega_r}{L} \sin(p\theta_r) \
\frac{di_\beta}{dt} &= -\frac{R_s}{L} i_\beta + \frac{1}{L} v_\beta - \frac{\psi_f \omega_r}{L} \cos(p\theta_r) \
\frac{d\omega_r}{dt} &= \frac{1}{J} \left( \frac{3}{2} p \psi_f [i_q] - T_L \right) \
\frac{d\theta_r}{dt} &= \omega_r
\end{aligned}
]

其中：

• (i_q = -i_\alpha \sin(p\theta_r) + i_\beta \cos(p\theta_r))
• (T_L)：负载转矩（视为常数或慢变扰动）
• (p)：极对数，(J)：转动惯量

关键非线性项：(\sin(p\theta_r)), (\cos(p\theta_r)), (i_q)

输出方程（可测）：

[
\mathbf{y} = \begin{bmatrix} i_\alpha \ i_\beta \end{bmatrix} = h(\mathbf{x}) = \begin{bmatrix} x_1 \ x_2 \end{bmatrix}
]

三、EKF 算法设计（离散时间）

设采样周期 (T_s = 100,\mu\text{s})，离散化步骤如下：

1. 状态预测（Time Update）

[
\hat{\mathbf{x}}{k|k-1} = f(\hat{\mathbf{x}}{k-1|k-1}, \mathbf{u}_{k-1})
]

其中 (f(\cdot)) 为离散化状态转移函数（可用欧拉法）：

% 欧拉离散化示例
i_alpha_next = i_alpha + Ts*(-Rs/L*i_alpha + 1/L*v_alpha + psi_f*omega_r/L*sin(p*theta_r));
...

预测误差协方差：
[
\mathbf{P}_{k|k-1} = \mathbf{F}k \mathbf{P}{k-1|k-1} \mathbf{F}_k^\top + \mathbf{Q}
]

其中 (\mathbf{F}k = \left. \frac{\partial f}{\partial \mathbf{x}} \right|{\hat{\mathbf{x}}_{k-1|k-1}}) 为雅可比矩阵。

2. 量测更新（Measurement Update）

计算卡尔曼增益：
[
\mathbf{K}k = \mathbf{P}{k|k-1} \mathbf{H}_k^\top (\mathbf{H}k \mathbf{P}{k|k-1} \mathbf{H}_k^\top + \mathbf{R})^{-1}
]

状态修正：
[
\hat{\mathbf{x}}{k|k} = \hat{\mathbf{x}}{k|k-1} + \mathbf{K}_k (\mathbf{y}k - h(\hat{\mathbf{x}}{k|k-1}))
]

协方差更新：
[
\mathbf{P}_{k|k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}k) \mathbf{P}{k|k-1}
]

其中 (\mathbf{H}k = \left. \frac{\partial h}{\partial \mathbf{x}} \right|{\hat{\mathbf{x}}_{k|k-1}} = \begin{bmatrix} 1 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 \end{bmatrix})

四、关键矩阵计算（雅可比）

状态转移雅可比 (\mathbf{F} = \partial f / \partial \mathbf{x})

经推导（略去中间步骤），得：

[
\mathbf{F} = \begin{bmatrix}
-\frac{R_s}{L} & 0 & a & b \
0 & -\frac{R_s}{L} & c & d \
0 & 0 & 0 & e \
0 & 0 & 1 & 0
\end{bmatrix}
]

其中：

• (a = \frac{\psi_f}{L} \sin(p\theta_r))
• (b = \frac{\psi_f \omega_r p}{L} \cos(p\theta_r))
• (c = -\frac{\psi_f}{L} \cos(p\theta_r))
• (d = \frac{\psi_f \omega_r p}{L} \sin(p\theta_r))
• (e = \frac{3}{2} \frac{p^2 \psi_f}{J} [-i_\alpha \cos(p\theta_r) - i_\beta \sin(p\theta_r)])

注意：(\mathbf{F}) 在每个采样周期需重新计算！

五、Simulink 建模步骤

第一步：搭建 PMSM 主系统

• 使用 Simscape Electrical 的 Permanent Magnet Synchronous Motor
• 参数示例：
  • (R_s = 1.2,\Omega), (L_d = L_q = 8,\text{mH})
  • (\psi_f = 0.175,\text{Wb}), (p = 4), (J = 0.0008,\text{kg·m}^2)
• 逆变器 + 电流采样（abc → αβ）

第二步：实现 EKF 核心模块

创建子系统 EKF_Estimator，内部使用 MATLAB Function 模块实现完整 EKF 循环：

function [theta_hat, omega_hat] = fcn(v_alpha, v_beta, i_alpha, i_beta, Ts)
persistent x_hat P Q R F H I

if isempty(x_hat)
    % 初始化
    x_hat = [0; 0; 0; 0];  % [i_a; i_b; w_r; theta_r]
    P = diag([1e-3, 1e-3, 100, 1]);  % 初始协方差
    Q = diag([1e-5, 1e-5, 0.1, 0.01]); % 过程噪声
    R = diag([1e-4, 1e-4]);            % 测量噪声
    I = eye(4);
end

% --- Time Update ---
% 状态预测 (欧拉法)
w_r = x_hat(3); theta_r = x_hat(4);
i_a_pred = x_hat(1) + Ts*(-1.2/0.008*x_hat(1) + 1/0.008*v_alpha + 0.175*w_r/0.008*sin(4*theta_r));
i_b_pred = x_hat(2) + Ts*(-1.2/0.008*x_hat(2) + 1/0.008*v_beta - 0.175*w_r/0.008*cos(4*theta_r));
w_r_pred = x_hat(3);  % 忽略负载转矩变化
theta_r_pred = x_hat(4) + Ts * x_hat(3);
x_pred = [i_a_pred; i_b_pred; w_r_pred; theta_r_pred];

% 计算 F 矩阵
p = 4; L = 0.008; psi = 0.175; J = 0.0008;
a = psi/L * sin(p*theta_r_pred);
b = psi*w_r_pred*p/L * cos(p*theta_r_pred);
c = -psi/L * cos(p*theta_r_pred);
d = psi*w_r_pred*p/L * sin(p*theta_r_pred);
e = 1.5 * p^2 * psi / J * (-x_pred(1)*cos(p*theta_r_pred) - x_pred(2)*sin(p*theta_r_pred));
F = [-1.2/L, 0,      a, b;
      0,     -1.2/L, c, d;
      0,      0,     0, e;
      0,      0,     1, 0];

% 预测协方差
P_pred = F * P * F' + Q;

% --- Measurement Update ---
y = [i_alpha; i_beta];
y_pred = [x_pred(1); x_pred(2)];
H = [1, 0, 0, 0; 0, 1, 0, 0];
S = H * P_pred * H' + R;
K = P_pred * H' / S;
x_hat = x_pred + K * (y - y_pred);
P = (I - K*H) * P_pred;

% 输出
theta_hat = x_hat(4);
omega_hat = x_hat(3);
end

注意：实际应用中应将电机参数设为输入，而非硬编码。

第三步：闭环控制集成

• 将 theta_hat 用于 Park/Clarke 变换
• 将 omega_hat 用于速度环反馈
• 构建完整无感矢量控制系统

六、仿真设置与结果

测试场景：加速 → 恒速 → 负载突变

• 0–0.2 s：0 → 1500 rpm 加速
• 0.2–0.4 s：恒速 1500 rpm
• t=0.3 s：负载从 2 N·m → 5 N·m

关键波形：

性能指标：

优势：相比龙伯格观测器，EKF 在低速和噪声下表现更优。

七、EKF vs 其他状态估计算法

EKF 优势：在“精度”与“鲁棒性”之间取得最佳平衡，适合高性能驱动。

八、工程优化建议

1. 噪声协方差整定：

   • (Q) 大 → 信任测量，响应快但噪
   • (R) 大 → 信任模型，平滑但滞后
   • 建议：通过实验试凑，或使用自适应 EKF

2. 离散化方法升级：

   • 用 四阶龙格-库塔（RK4）替代欧拉法，提升精度

3. 负载转矩估计：

   • 将 (T_L) 加入状态向量，实现扰动抑制

4. 定点实现：

   • 使用 Fixed-Point Designer 避免浮点溢出
   • 对矩阵求逆用 Cholesky 分解 提升数值稳定性

九、总结

本教程完成了：

1. 建立了 PMSM 的非线性状态空间模型
2. 推导了离散 EKF 算法及关键雅可比矩阵
3. 在 Simulink 中通过 MATLAB Function 实现了完整 EKF
4. 验证了其在宽速域、负载扰动下的高精度估计能力

该技术适用于：

• 电动汽车主驱电机
• 工业伺服系统
• 航空航天作动器

核心思想：
“以概率之眼，观非线性之变” —— 在不确定中寻找最优估计，于噪声里提取真实状态。

十、动手建议

1. 尝试 无迹卡尔曼滤波（UKF）进一步提升非线性精度
2. 对比 EKF vs SMO 在高频噪声下的表现
3. 加入 参数在线辨识（如 (\psi_f) 漂移）
4. 使用 Simulink Coder 生成嵌入式 C 代码部署到 DSP

通过本模型，你已掌握现代最优估计理论在电机控制中的高级应用——EKF 无感智能驱动。