做伺服驱动调试这些年转动惯量辨识这个问题算是绕不开的硬骨头。速度环PI参数整定得准不准很大程度取决于你对负载惯量掌握得怎么样尤其是永磁同步电机PMSM拖动的负载经常变化时固定一组参数很难兼顾动态响应和稳态精度。扩展卡尔曼滤波EKF作为非线性状态估计的经典工具被不少人拿来做惯量在线辨识而在MATLAB/Simulink环境中用S函数实现EKF算法原型是一条被验证过很多次的成熟路径。这篇文章结合我实际调试中的理解和踩过的坑把EKF做PMSM转动惯量辨识的思路、状态方程构建、S函数实现细节、以及阅读文献时需要留意的坑位一次性讲清楚。适合刚接触参数辨识的工程师、做算法验证的研究生以及想把仿真算法搬到嵌入式代码里的开发者参考。1. 转动惯量辨识为什么它值得专门做1.1 速度环整定绕不开的物理量先回到最基本的运动方程。永磁同步电机拖动负载时机械运动方程可以写成J * dω_m/dt T_e - T_L - B * ω_m这里的J是折算到电机轴上的总转动惯量包含转子自身惯量加上负载折算值T_e是电磁转矩T_L是负载转矩B是粘滞摩擦系数。从控制角度看转动惯量决定了一个最基础的输入输出关系在同样的电磁转矩作用下J越大速度变化越慢J越小速度变化越陡。这个特性直接决定了速度环PI参数整定的方向和取值。速度环在设计PI控制器的时候常见做法是把电流环近似成一阶惯性环节甚至直接视为理想环节那么速度环的开环传递函数可以化简成类似 Kp*(sKi/Kp 1) / (Js) 的形式。为了让系统获得较好的阶跃响应工程上通常会让速度环积分时间常数与机电时间常数匹配最终Kp和Ki的取值基本与J成正比。换句话说如果你把J估小了一半速度环实际增益会比设计值大一倍系统的相位裕度下降带轻载时只是反应快一点带重载时可能就直接振荡反过来把J估大了一倍响应会明显变钝电机给人的感觉就是“软绵绵”的跟踪误差变大。这个物理过程一点也不抽象。我曾经遇到一个数控转台的项目空载惯量和满载惯量差了将近三倍用同一组速度环参数分别跑空载和满载空载时电流谐波明显、机械噪音刺耳满载时跟踪误差增加了一倍多。后来就是靠在线惯量辨识配合增益调度把这个问题解决的。所以惯量辨识不是实验室里的理论游戏它是工程整定的前置条件。1.2 传统辨识方法为什么不够用转动惯量辨识不是新问题早在线辨识概念流行之前就有多种离线方案。最朴素的是加减速法给电机一段匀加速指令测量转速变化量和电磁转矩利用运动方程反推J。具体做法是让电机从静止加速到某个目标速度记录q轴电流和转速斜率。这个方法在实验室环境没问题但在线场合就很尴尬——系统运行中很难刻意去创造匀加速区间而且负载转矩本身在波动时单点测得的转矩差里混着负载扰动辨识精度得不到保障。实际实现中还需要对转速做差分来获取加速度数字微分会把高频噪声放大结果往往很毛糙。模型参考自适应MRAS也是常见方案。它构造一个参考模型和一个可调模型通过自适应律推动可调模型参数逼近参考模型。MRAS做PMSM参数辨识的资料非常多但它的收敛速度和稳定性受自适应增益影响很大增益取小了收敛慢取大了又容易在噪声激励下发散。需要快速响应变负载的场合MRAS在调参上花的心思一点不比EKF少。EKF的思路则不同。它本质上是把标准卡尔曼滤波推广到非线性系统状态方程和观测方程都允许是非线性的而PMSM的电气方程和机械方程恰恰是典型的非线性结构。EKF能同时处理过程噪声和传感器测量噪声把它们各自的统计特性融进增益计算中这意味着噪声不再只是被“滤掉”的东西而是成为估计器做决策时的信息依据。把它用在惯量辨识里把J作为扩展状态放进状态向量在线估计就不需要额外注入激励信号也不打断正常控制运行这是它相比前两类方法最舒服的一点。需要强调的是EKF对非线性系统做的是局部线性化——在每个采样周期用一阶泰勒展开的雅可比矩阵来近似状态转移所以它是“次优”估计器。好在PMSM的机械动态相对电气动态慢得多采样周期内线性化误差通常可以接受工程应用效果足够。2. EKF辨识方案的整体设计与状态方程构建2.1 增广状态向量的两种构建思路做EKF的第一步是确定状态向量。在PMSM惯量辨识这个题目下大体上有两种层次的选择。第一种是只做机械子系统辨识状态向量取为 x [ω_m, T_L, J]即机械角速度、负载转矩、转动惯量。这种做法的前提是电流环闭环带宽够高可以把电磁转矩直接用q轴电流指令换算得到相当于把电流环当成一个理想环节。好处是状态维数低只有3维雅可比矩阵是3x3Q矩阵和R矩阵整定起来直观很多代码也短。坏处是如果电流环带宽不高或者运行转速很高导致电流环动态明显滞后用电流指令代替实际转矩就会引入相位误差惯量估计值出现系统性偏差。所以这个方案一般适合中低速、电流环带宽充裕的场合。第二种是电气量和机械量耦合在一起的全阶模型状态向量取 x [i_d, i_q, ω_m, J]有的还会把定子电阻R_s、永磁体磁链ψ_f甚至负载转矩T_L一起增广上去。全阶模型把电气方程和机械方程放在同一个系统里理论上最完整能直接利用电流残差和速度残差共同修正惯量高频信息更丰富。但状态维数一提高雅可比矩阵的推导和Q/R整定的难度都会显著上升而且各状态之间的量纲差异很大——电流是安培级别惯量可能小于10的负四次方数量级数值上容易病态。对刚接触EKF的读者我的建议很直接先做机械子系统的三阶模型把整个流程跑通理解清楚EKF每一行代码在干什么再考虑扩展成全阶模型。一上来就做六维状态估计很容易被调参和矩阵推导淹没反而不利于掌握算法的本质。2.2 状态方程与雅可比矩阵的推导拿机械子系统模型来列写。在dq同步旋转坐标系下表贴式PMSML_d L_q L_s的电磁转矩为T_e 1.5 * p * ψ_f * i_q其中p是极对数ψ_f是永磁体磁链。如果把i_q当作已知输入来自电流环指令状态方程写为dω/dt (1.5*p*ψ_f*i_q - T_L - B*ω) / J dT_L/dt 0 dJ/dt 0负载转矩和转动惯量都假设为随机游走模型慢变化量导数为零但伴随过程噪声。测量方程取 y ω也就是状态向量的第一个分量这个三阶系统是EKF最标准的结构形式。看起来简单但雅可比矩阵推导时最容易栽跟头。对f(x)求偏导得到F [ -B/J, -1/J, -(1.5*p*ψ_f*i_q - T_L - B*ω) / J^2 ] [ 0, 0, 0 ] [ 0, 0, 0 ]注意F的第一行第三列这是对J求偏导得到的分式数值上等于“净转矩除以J的平方”再取负号。这一项是EKF里J的修正能够与转速残差耦合起来的关键通道漏掉它J的估计虽然不会立刻发散但收敛特性和稳态方差都会明显不对属于最难排查的隐性错误。如果做全阶模型状态方程还要在i_d和i_q上补两行di_d/dt (-R_s*i_d ω_e*L_s*i_q u_d) / L_s di_q/dt (-R_s*i_q - ω_e*L_s*i_d - ω_e*ψ_f u_q) / L_s dω/dt (1.5*p*ψ_f*i_q - T_L - B*ω) / J这里的ω_e pω是电角速度。注意di_q/dt的反电动势项里ψ_f乘上ω_e线性化时对ω求偏导会产生 pψ_f 项这很好理解就是反电动势与转速耦合的直接体现。我建议推导雅可比矩阵时先用符号计算工具验证一遍再手推核对两种途径结果一致后再写进代码能省下大量排查时间。2.3 EKF离散递推流程与数值积分EKF在计算机里实现时是离散递推的。预测步根据系统模型向前推一步x̂_k|k-1 x̂_{k-1} T_s * f(x̂_{k-1}) P_k|k-1 F_k * P_{k-1} * F_k^T Q这里的积分本质上是数值积分。采样周期T_s相对系统动态来说已经很小欧拉法凑合能用但实际调试中我发现在高速段比如额定转速附近甚至进入弱磁区时反电动势项变化剧烈欧拉法一步积分产生的截断误差会污染残差惯量估计值的抖动会明显增大。改用RK2中点法甚至RK4能让预测步更贴近真实非线性动态计算开销对现代控制芯片来说完全在可接受范围内我自己的工程代码里至少会用RK2。更新步就是标准的卡尔曼滤波更新公式S_k H * P_k|k-1 * H^T R K_k P_k|k-1 * H^T * S_k^(-1) x̂_k x̂_k|k-1 K_k * (y_k - h(x̂_k|k-1)) P_k (I - K_k * H) * P_k|k-1整个EKF代码确实不长核心难点不在公式本身而在三个矩阵的取值过程噪声协方差Q怎么给、测量噪声协方差R怎么给、初始协方差P0怎么给。这三个矩阵直接决定辨识的动态性能和稳态精度后面我会专门用一节来展开。3. S函数实现EKF的核心细节3.1 回调结构与EKF的对应关系在MATLAB/Simulink里实现EKF最直接的方式是写Level-2 MATLAB S-Function。很多初学者一打开S函数的模板就被一堆回调函数吓住其实只需要理清楚四个关键回调。mdlInitializeSizes负责声明输入输出端口的数量、尺寸和直接馈通属性mdlInitializeSampleTimes负责设置采样时间mdlUpdate是算法主体所在在每个固定步长触发EKF的预测和更新都放这里mdlOutputs负责把估计状态写到输出端口。有一个非常容易踩的坑需要提前说明EKF的预测步本质上是在做数值积分而Simulink里的S函数如果启用了连续状态在mdlDerivatives里返回状态导数那积分步长会由Simulink求解器控制不一定等于你的控制周期。如果EKF内部还有自己的一套离散递推两者就会互相打架。我见过不少初学者在同一个S函数里既写了mdlDerivatives又写了mdlUpdate结果仿真时收敛时而不收敛根源就是求解器步长和算法离散周期混在了一起。最干净的做法是把EKF完全写成离散形式在mdlUpdate里手动完成预测和更新采样时间设为固定步长且与仿真模型的步长保持一致。3.2 Level-2 S-Function框架代码下面给出一个基于三阶机械子系统模型的EKF辨识S函数骨架。重点展示回调函数如何组织EKF的预测与更新代码里我做了一些简化但整体结构是按实际调试经验拆的。function EKF_Inertia(block) % Level-2 MATLAB S-Function: EKF based PMSM inertia identification % 输入: u(1)iq电流指令, u(2)实际机械转速 % 输出: y(1)估计转速, y(2)估计负载转矩, y(3)估计转动惯量 % 对话框参数: [1] Rs, [2] Ls, [3] 磁链, [4] 极对数, [5] 粘滞系数B, [6] 采样周期Ts setup(block); function setup(block) block.NumInputPorts 2; block.NumOutputPorts 3; block.InputPort(1).Dimensions 1; block.InputPort(2).Dimensions 1; block.OutputPort(1).Dimensions 1; block.OutputPort(2).Dimensions 1; block.OutputPort(3).Dimensions 1; block.NumContStates 0; block.NumDworks 2; % Dwork1: 状态向量x, Dwork2: 协方差矩阵P按列展开 block.SampleTimes [0.0001 0]; % 与仿真步长保持一致 block.NumDialogPrms 6; block.RegDwork(DworkState, 3); % x [wm; TL; J] block.RegDwork(DworkCov, 9); % P 按列存放 3x3 block.SetPreCompOutPortParam(1, DworkState, 1); block.SetPreCompOutPortParam(2, DworkState, 2); block.SetPreCompOutPortParam(3, DworkState, 3); end function mdlInitializeSampleTimes(block) block.SampleTimes [block.DialogPrms(6) 0]; end function mdlUpdate(block) prm.Rs block.DialogPrms(1); prm.Ls block.DialogPrms(2); prm.flux block.DialogPrms(3); prm.p block.DialogPrms(4); prm.B block.DialogPrms(5); Ts block.DialogPrms(6); x block.Dwork(1).Data; % 3x1 P reshape(block.Dwork(2).Data, 3, 3); % 3x3 u block.InputPort(1).Data(1); % iq wm_meas block.InputPort(2).Data(1); % --- 预测步RK2离散化 --- [f1, F1] sysDyn(x, u, prm); x_mid x 0.5 * Ts * f1; [f2, ~] sysDyn(x_mid, u, prm); x_pred x Ts * f2; F eye(3) Ts * F1; % 一阶离散化雅可比工程上可近似 P_pred F * P * F Qmatrix(x, prm); % --- h(x) wm --- H [1 0 0]; y_pred x_pred(1); S H * P_pred * H Rmatrix(wm_meas, prm); K P_pred * H / S; innovation wm_meas - y_pred; x_upd x_pred K * innovation; P_upd (eye(3) - K * H) * P_pred; % 保证数值合理性 if x_upd(3) 1e-6 x_upd(3) 1e-6; end block.Dwork(1).Data x_upd; block.Dwork(2).Data P_upd(:); end function mdlOutputs(block) x block.Dwork(1).Data; block.OutputPort(1).Data x(1); block.OutputPort(2).Data x(2); block.OutputPort(3).Data x(3); end function [f, F] sysDyn(x, iq, prm) wm x(1); TL x(2); J x(3); Te 1.5 * prm.p * prm.flux * iq; f zeros(3,1); f(1) (Te - TL - prm.B * wm) / J; f(2) 0; f(3) 0; F zeros(3,3); F(1,1) -prm.B / J; F(1,2) -1 / J; F(1,3) -(Te - TL - prm.B * wm) / (J * J); end这段代码里Qmatrix和Rmatrix需要根据被控对象和采样周期具体填写后面我会讲怎么调。Dwork的用法很重要它相当于S函数内部自己的持久化存储不像全局变量那样容易被其他模块污染也比连续状态更可控推荐用Dwork来保存EKF的状态向量和协方差矩阵。3.3 Simulink仿真搭建与信号处理要点搭建仿真模型时我习惯分几个层次来组织底层放PMSM的电气机械模型可以用Simulink自带的电机模块也可以自己搭上层放电流环PI控制器和SVPWM逆变器再上层放速度环EKF辨识模块从电流环和速度环取信号输出估计惯量再接一个增益调度模块在线更新速度环PI参数。一个容易犯的错误是信号类型不匹配。EKF模块内部做矩阵运算默认是double类型但Simulink模型里如果混入了布尔信号或者整数信号S函数端口会报类型错误或者悄悄做类型转换影响结果。建议在S函数的每个输入端口前加Data Type Conversion模块强制转成double。另外各种PI模块的输出限幅一定要处理好。我踩过这样一个坑电流环PI在启动阶段饱和输出的iq指令被限幅但反馈给EKF的“电流”信号取的是PI输出未限幅前的值结果算法以为电机输出了很大的转矩实际上电机并没有输出那么多J的估计值在启动瞬间就冲偏了。等限幅结束残差才慢慢把估计拉回来。这种问题在仿真曲线里可能只是一点点毛刺但放到真实平台就是启动阶段一次不小的扰动。统一做法是EKF的输入信号必须取饱和限幅之后的实际指令值。4. 参考资料解析如何把论文方法搬进工程4.1 建模粒度差异决定了调参难度搜索“EKF永磁同步电机惯量辨识”这类关键词能看到不少期刊论文和学位论文。大部分论文都会给出状态方程、雅可比矩阵和仿真波形但论文之间的状态选择差异很大。常见的取法就是前面说的三类只含机械量的三阶模型、含电气电流的四阶模型、以及把电阻磁链一起增广的更高阶模型。我读文献的习惯是第一眼先看状态向量维数因为维数直接决定了Q/R矩阵整定的复杂度。不要被“状态越多越全面”这种说法带着走。对惯量辨识这个具体目标来说如果电流环动态确实够快三阶模型往往比六阶模型稳得多调试成本低一个量级。我实际接触过一个团队用六阶EKF调了一个多月也没稳定收敛换回三阶模型一两天就出了可用的结果。这并不是说高阶模型不好而是说在工程节奏里你要先分清自己的需求边界你是在做控制器参数整定需要惯量还是在做电机参数全面自学习需求不同模型选择的逻辑完全不同。4.2 仿真结果与实验平台的差距这是资料解析里最重要的一条经验。论文里的仿真结果绝大多数是在标称参数、理想传感器、高斯白噪声的假设下跑出来的。仿真环境里PMSM的电阻、电感、磁链都是固定值传感器没有零漂逆变器没有死区那么EKF的收敛曲线当然是又光滑又漂亮的。但真实平台上有四件事几乎一定会破坏这种理想状态电流采样零漂、码盘量化误差、逆变器死区导致的电压畸变、以及温度变化引起的电机参数漂移。这些干扰叠加在EKF的“测量噪声”里但它们不是高斯白噪声而是有偏噪声。尤其电流采样零漂会让dq变换后的i_d/i_q带上常值偏置EKF基于偏置后的值做递推惯量估计就会带着偏置逆变器死区效应则会让实际加载电机上的电压与控制器输出的参考电压不一致这种电压误差在低速时尤其明显。所以看论文的时候不能只抄状态方程和增益矩阵还要注意它的噪声模型是怎么建立的是否考虑了死区补偿是否对采样偏置做了校正。工程上我的经验是先做死区补偿或者在辨识算法里对电压误差建模否则直接照抄论文方案上机结果往往比仿真差一大截。4.3 复现论文时最容易被忽略的初始化与数据链路问题论文对初始状态x0和协方差矩阵P0的处理通常就一句话“初始协方差取适当值”但这恰恰是复现时最坑人的地方。惯量J的物理量级可能只有10的负四次方转速是几百rad/s电流是几安培这些状态量之间差了好几个数量级。如果P0对角线统一设成1那么预测步里P矩阵与F矩阵、Q矩阵耦合后增益K会先被初始协方差主导J的初值修正会产生一个很大的瞬时跳变严重时会激起速度环振荡。规范做法是按物理量级分别初始化P0对角线。比如转速对应的P0(1,1)取1到10负载转矩对应的P0(2,2)取10到100惯量对应的P0(3,3)取0.0001量级。原则就是不要让P0各分量拉开巨大数量级差也不要让不同物理量共用同一个数值基准。论文里很少写这种细节但调试时几乎都会碰到。另一个被忽视的问题是数据链路一致性。EKF的观测方程是 y ω那实际送到算法端口的信号就必须是纯粹的机械角速度如果这个信号经过滤波、经过标幺变换、或者经过了单位换算观测方程也要同步修改。前面提到转速低通滤波带来的相位延迟问题就是典型例子模型和信号链不一致EKF的残差就永远带着系统误差惯量估计值自然不准。5. 参数调试与常见问题排查5.1 Q/R矩阵整定的实际操作顺序Q和R是EKF里最核心的自由参数也是最让人头疼的参数。先把物理含义说透Q代表你对系统模型不确定性的估计R代表你对传感器测量不确定性的估计。在惯量辨识这个场景里Q矩阵中对应J的那个元素随机游走方差直接决定J的跟踪速度R中对应速度测量的元素决定滤波对速度残差的信任程度两者之间的比值决定卡尔曼增益的整体水平。我调试的顺序一直很固定。第一步先定R_ω让电机在匀速段运行用示波器或数据记录仪统计转速反馈的标准差这个数值就是R_ω的最小合理下界。第二步扫Q_J从小到大逐步增加Q_J观察J估计值的阶跃响应时间和稳态抖动幅度找到一个“响应够快、抖动可接受”的交点。第三步定Q_TL负载转矩对应的过程噪声通常取0.01到1之间的值太小则负载突变时T_L跟不上J会代偿性地跳变太大会让T_L估计里混入大量噪声J反而受影响。整体逻辑是先让T_L估计稳定再观察J的响应。一个很实用的离线调试手段是用实验数据回放来扫参。在Simulink里搭一个开环回放模型把实际平台录制的iq指令、转速、电流采样作为信号源喂给EKF模块然后批量扫Q_J和R_ω看J估计值的稳态方差和响应时间选出折中值后再上在线平台。这个方法能省掉大量在线调试时间而且排查问题也更安全。5.2 典型故障现象速查表我把调试EKF惯量辨识时最常遇到的故障现象、可能原因和排查方向整理成一张表方便对照。现象可能原因排查方向J收敛到负值或极小值运动方程符号错误、雅可比矩阵错项、坐标变换方向不一致核对运动方程正负号检查坐标变换方向J收敛正常但响应慢Q_J过小或T_L估计没有跟上增大Q_J检查T_L估计跟踪情况J估计值抖动剧烈Q_J过大或R_ω过小适当降低Q_J或增大R_ω阶跃负载后J先跳变后恢复T_L的过程噪声设置过小T_L跟不上负载突变增大Q_TL或者把T_L增广为状态单独配置方差高速段J漂移明显反电动势线性化误差累积、电流环带宽不足改用RK2/RK4积分检查电流环带宽稳态时J与标称值偏差较大电流采样零偏、逆变器死区、电阻磁链温漂增加死区补偿、校准电流零偏、考虑参数在线修正这张表里的现象我基本都亲手遇到过。如果读者遇到的情况不在表里建议按照“先核对模型与代码的一致性、再查单位与坐标变换、最后调Q/R”的顺序排查。模型错了去调Q/R等于在错误的地基上折腾装修白费力气。5.3 离线回放验证与交叉校验流程最后说验证流程。做完仿真后不要直接上实物盲调我惯用的流程是先在运行中的平台录制一段输入信号回放数据把q轴电流指令、实测转速、电流采样原始数据都保存下来然后在Simulink里搭开环回放模型用真实数据喂给EKF辨识模块对比J的估计曲线和理论预期。这样做的好处是EKF模块的行为可以完全在离线环境里调试不占用平台资源而且每次改动之后的结果都可复现不会因为现场环境变化而干扰判断。交叉验证方面最省事的校对手段是用加减速法离线测一次惯量作为参考基准。如果在线的辨识结果与离线基准偏差在5%以内基本可以认定EKF的模型和参数是靠谱的偏差偏大就回头检查状态方程、坐标轴定义和传感器数据链路而不是继续盲目调Q/R。这里有个细节值得说加减速法离线测量时负载状态必须和在线测试时保持一致。有一次我做交叉校验时偏差超过10%排查了很久才发现离线测试时转台上多装了一个夹具惯量本身就变了。测量条件不一致导致的偏差怪不到算法头上。这种低级错误最容易在赶进度的时候出现先统一测试条件再下结论。我个人在实际操作中的体会是EKF惯量辨识这个任务本身并不玄学绝大多数问题出在模型一致性、单位换算和数值处理上。调试中我曾经花很长时间排查J收敛后数值偏小的问题最后发现是转速反馈信号在前面加了一个低通滤波器滤波造成的相位延迟让EKF的观测残差混入了动态误差J的估计便被系统性拉低。这也再次印证了那条原则EKF的观测模型必须和你实际送到算法里的测量值严格匹配任何预处理环节都不能随意加。后来我把这个S函数模块的参数全部接口化换不同功率等级的PMSM模型时只需要改电机参数和Q/R初值代码结构完全不用动项目复用性一下子提升了很多。如果你也正卡在雅可比矩阵推导或者Q/R整定上我建议先从机械子系统三阶模型入手把离线回放验证跑通再考虑全阶扩展。这是我自己验证过最稳妥的一条路。
