简介基于模型预测控制MPC的永磁同步电机Matlab代码面向电机控制方向的研究生、工程师及算法学习者解决MPC策略在PMSG永磁同步电机调速系统中的快速仿真验证问题。压缩包共2个m文件、容量仅2KBPMSG_MPC1.m实现单电流环MPC速度环采用PI调节器PMSG_MPC2.m实现速度环与电流环均用MPC的非级联全MPC策略。电机本体以数学模型代替Simulink物理模型便于查看算法细节和修改参数轻量紧凑适合MPC入门与对比实验。目前已有826人学习下载。读者可直观对比两种MPC架构的控制差异学习预测控制器与电机数学模型结合的程序实现思路为后续扩展自定义优化目标提供基础。1. 模型预测控制MPC在永磁同步电机里到底控制什么先给一个反直觉的结论MPC 在电机控制里并不是一个“全能最优”的黑盒它的强项是把“约束”和“目标”直接塞进控制器里代价是要在每个采样周期内反复做预测和寻优。和 PI 相比MPC 的调参逻辑完全不同——PI 整定的是增益MPC 整定的是代价函数里的权重。所以你会看到很多人说“MPC 效果不错但到底怎么调权重、怎么在 Matlab 里落地代码和理论对不上”。这篇文章要讲的就是一套能直接放到 Matlab/Simulink 里跑的永磁同步电机电流环 MPC 方案包含离散预测模型推导、代价函数设计、8 个电压矢量的扫描式求最小值以及采样周期、权重系数和延时补偿这些真正影响仿真能不能收敛的关键点。无论你是刚接触 MPC 的从业者还是想把仿真结果搬上 DSP 做硬件迭代的工程师下面的内容都按“理论 → 代码 → 排错”的顺序展开每段代码都有对应的参数说明。2. 永磁同步电机离散预测模型与代价函数设计2.1 永磁同步电机在 dq 旋转坐标系下的连续模型先把被控对象写清楚。表贴式永磁同步电机在 dq 旋转坐标系下的电压方程是u_d Rs * i_d Ld * di_d/dt - ωe * Lq * i_q u_q Rs * i_q Lq * di_q/dt ωe * (Ld * i_d ψf)其中 Rs 是定子电阻Ld、Lq 是 d、q 轴电感ψf 是永磁磁链ωe 是电角速度。电磁转矩方程为Te 1.5 * p * ψf * i_q表贴式忽略磁阻转矩机械方程为J * dωm/dt Te - B * ωm - TL。大部分 MPC 仿真文章用的是表贴式电机也就是 Ld ≈ Lq。这有个好处预测模型里不需要处理磁阻转矩项代价函数里 d、q 轴电流可以独立评估初学阶段不至于被耦合项绕晕。如果用的是内置式电机IPMSM模型里的电感矩阵和转矩方程都要加上(Ld - Lq)项代价函数也会更复杂后面代码里的预测方程需要同步换掉。2.2 前向欧拉离散化把连续方程变成可迭代的预测公式MPC 的核心是“预测”也就是用当前时刻的电流和电压估算下一个采样周期的电流。微控制器和 Simulink 仿真里最常见的做法是前向欧拉离散化把di/dt近似成(i(k1) - i(k)) / Ts代入上面的电压方程就可以解出预测电流i_d(k1) i_d(k) Ts/Ld * (u_d(k) - Rs*i_d(k) ωe(k)*Lq*i_q(k)) i_q(k1) i_q(k) Ts/Lq * (u_q(k) - Rs*i_q(k) - ωe(k)*(Ld*i_d(k) ψf))注意两个容易写错的地方。第一个是交叉耦合项d 轴方程里是 ωe * Lq * i_qq 轴方程里是- ωe * Ld * i_d符号写反了会导致预测的电流方向错误仿真里表现为电流发散。第二个是反电动势项ωe * ψf只在 q 轴方程里出现它代表转子旋转带来的反电动势如果漏掉这一项高速工况下 MPC 选择的电压矢量会不准确。另一个容易被忽略的是角度同步问题。上面公式里的 u_d、u_q 是同一个时刻的电压在 dq 轴上的投影而 dq 轴是随转子旋转的。实际实现时你拿到的是直流母线电压和逆变器开关状态要先通过 Clark 变换得到 αβ 轴电压再用当前转子角度做 Park 变换得到 u_d、u_q。也就是说MPC 每一拍都要读一次角度 θe这也是它比 FOC 对编码器更敏感的原因。2.3 有限控制集 MPC 与代价函数设计永磁同步电机驱动用的是电压源逆变器三相桥臂只有 8 种开关状态不能任意输出连续的电压幅值。所以这里用的 MPC 叫有限控制集 MPCFCS-MPC它不依赖传统意义上的 MPC 求解器去解二次规划而是把 8 个候选电压矢量全部预测一遍选代价函数最小的那个直接输出。8 个开关状态对应的 αβ 轴电压可以用下面公式算出u_alpha Vdc / 3 * (2*Sa - Sb - Sc); u_beta Vdc / sqrt(3) * (Sb - Sc);其中 Sa、Sb、Sc 是三个桥臂的开关状态取值 0 或 1。把 8 种组合代入会得到两个零矢量000 和 111和 6 个有效矢量幅值为2/3 * Vdc相邻矢量夹角 60 度。下面这张表是 8 个开关状态在 αβ 轴下的电压分量表达式仿真时可以直接用这个表做信号转换开关状态 (Sa,Sb,Sc)uαuβ对应空间矢量00000V0零矢量1002Vdc/30V1110Vdc/3Vdc/√3V2010-Vdc/3Vdc/√3V3011-2Vdc/30V4001-Vdc/3-Vdc/√3V5101Vdc/3-Vdc/√3V611100V7零矢量代价函数是所有 FCS-MPC 的“决策中心”控制目标越复杂代价函数项就越多。一个最朴素、也最能反映问题本质的代价函数是 d、q 轴电流跟踪误差的平方和J(i) (i_d_ref - i_d_pred(i))^2 (i_q_ref - i_q_pred(i))^2上面的 i 表示第 i 个候选电压矢量i_d_pred(i)、i_q_pred(i) 是套用 2.2 节离散公式预测出的下一个采样周期电流值。代价函数只选最小值对应的电压矢量这就是“滚动优化”在最简单形式下的落地方式。后面如果要限制开关频率可以在代价函数后面加一项λ * n_switchn_switch 是本周期开关状态与上一周期相比发生动作的次数权值 λ 越大开关频率越低但电流跟踪性能会变差这个权衡在第 4 章详细说。3. 从零写一套可运行的 Matlab 电流环 MPC 代码3.1 参数初始化脚本在 Matlab 里做永磁同步电机 MPC 仿真我一般先建立参数脚本再写预测函数最后搭主循环或 Simulink 模型。参数脚本是一切仿真的基础建议把电机参数和控制器参数分开方便后续做批量仿真。%% PMSM 参数 Rs 0.958; % 定子电阻Ω Ld 0.00525; % d 轴电感H Lq 0.00525; % q 轴电感H表贴式电机Ld Lq psif 0.1827; % 永磁磁链Wb p 4; % 极对数 J 0.000632; % 转动惯量kg·m^2 B 0.0003035; % 阻尼系数N·m·s Vdc 310; % 直流母线电压V %% MPC 采样参数 Ts 1e-5; % 电流环采样周期s对应 10 kHz 控制频率这里 Ts 是 MPC 控制周期不是 Simulink 里整个仿真模型的定步长。如果 Ts 取 1e-410 kHz那么每个周期有 100 微秒用来做预测、求代价和选矢量对 Matlab 仿真来说没问题但对后续移植到 DSP 会紧一些。通常我会先在仿真里用 1e-5 验证控制性能再换成 5e-5 做代码生成前的测试看看性能损失多少这个步骤能提前暴露延时问题。3.2 预测函数把离散方程写成可直接调用的函数预测函数是 MPC 最核心的计算单元不管你后续用脚本写主循环还是在 Simulink 的 MATLAB Function 块里实现都建议先把这一步封装成独立的函数方便测试和换参数。function [id_next, iq_next] pmsm_predict(id, iq, ud, uq, we, Ts, Rs, Ld, Lq, psif) % 前向欧拉离散 PMSM 电流模型 % 输入: 当前采样时刻电流、电压、电角速度 % 输出: 下一个采样时刻的预测电流 did (ud - Rs*id we*Lq*iq) / Ld; diq (uq - Rs*iq - we*(Ld*id psif)) / Lq; id_next id Ts*did; iq_next iq Ts*diq; end函数体只有四行但参数顺序和符号必须和 2.2 节公式完全一致。我在硬件移植时踩过的一个坑是仿真里电角速度 we 用的是弧度每秒而部分调试界面显示的是转每分rpm换算关系是we p * 2 * pi * rpm / 60。如果单位不统一预测模型算出的电流值会偏大代价函数的选择结果直接偏向某一个固定电压矢量电流波形看起来像“单极性调制”而不是正常的 MPC 对称输出。3.3 主循环里的 8 矢量扫描与代价函数评估脚本式仿真适合理解 MPC 的完整流程每拍读电流、角度、转速生成 8 个候选电压矢量逐个计算预测电流和代价函数找到最小值并输出对应的开关状态。% 8 个开关状态 states [0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 1 1; 0 0 1; 1 0 1; 1 1 1]; num_steps 50000; % 仿真步数 sa zeros(num_steps, 1); sb zeros(num_steps, 1); sc zeros(num_steps, 1); for k 1:num_steps % 读取当前状态这里用变量占位实际从被控对象模型取 id x_id(k); iq x_iq(k); theta x_theta(k); we x_we(k); % 转速外环给定 iq_ref第 5 章展开 id_ref 0; iq_ref speed_pi_output(k); J zeros(8, 1); for i 1:8 % 开关状态转 αβ 电压 ua Vdc/3 * (2*states(i,1) - states(i,2) - states(i,3)); ub Vdc/sqrt(3) * (states(i,2) - states(i,3)); % Park 变换到 dq 轴 ud ua * cos(theta) ub * sin(theta); uq -ua * sin(theta) ub * cos(theta); % 预测下一拍电流 [id_p, iq_p] pmsm_predict(id, iq, ud, uq, we, Ts, Rs, Ld, Lq, psif); % 代价函数 J(i) (id_ref - id_p)^2 (iq_ref - iq_p)^2; end [~, idx] min(J); sa(k) states(idx, 1); sb(k) states(idx, 2); sc(k) states(idx, 3); end这段代码里的 Park 变换方向要特别注意。这里的变换方向是ud ua*cos(theta) ub*sin(theta)对应的是 d 轴与 α 轴夹角为 theta 的常用坐标定义。如果你的电机模型中转子角定义不一样整个 MPC 的电压矢量选择都会错位表现出来的现象是低速能转但噪声大高速直接电流失控。建议先用一个固定角度的测试用例验证变换方向对不对比如 theta 0 时ud 应该等于 uauq 等于 -ub。3.4 Simulink 里用 MATLAB Function 块封装 MPC 控制器脚本仿真验证算法逻辑后最自然的下一步是搭 Simulink 模型把电机本体模型、逆变器和负载放到一起做闭环验证。Simulink 里的 PMSM 模块通常在 Simscape 或 Specialized Power Systems 库中逆变器可以用三个受控开关加母线电源实现也可以用简化的平均模型。MPC 控制部分我习惯放在 MATLAB Function 块里输入当前电流、角度、速度输出三路开关信号和脚本式实现的逻辑完全一致。function [sa, sb, sc] mpc_current_ctrl(id, iq, id_ref, iq_ref, theta, we, Ts, Rs, Ld, Lq, psif, Vdc) states [0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 1 1; 0 0 1; 1 0 1; 1 1 1]; J zeros(8, 1); for i 1:8 ua Vdc/3 * (2*states(i,1) - states(i,2) - states(i,3)); ub Vdc/sqrt(3) * (states(i,2) - states(i,3)); ud ua * cos(theta) ub * sin(theta); uq -ua * sin(theta) ub * cos(theta); [id_p, iq_p] pmsm_predict(id, iq, ud, uq, we, Ts, Rs, Ld, Lq, psif); J(i) (id_ref - id_p)^2 (iq_ref - iq_p)^2; end [~, idx] min(J); sa states(idx, 1); sb states(idx, 2); sc states(idx, 3); end在 MATLAB Function 块里调用 pmsm_predict 这个外部函数有个限制Simulink 的代码生成器要求在同一个文件夹下能找到函数文件或者把预测函数直接写成 MATLAB Function 块内的嵌套函数。更省事的方式是把预测公式直接复制进这个函数体避免文件路径问题。我习惯的做法是保留 pmsm_predict 函数做纯脚本调试Simulink 块里用内嵌版本。这里有几个 Simulink 模型层面的设置容易踩坑。一是步长仿真器的解算器要设成固定步长Ts否则变步长求解器会为了过零检测动态改变步长导致 MPC 的采样周期不稳定仿真结果看起来误码率很高。二是代数环如果把 PWM 开关信号直接拉回重复计算的模块Simulink 会产生代数环报错解决办法是在反馈通道里加一个 Unit Delay 块代表一拍计算延时这和真实数字控制器的行为一致。4. 采样周期、权重系数、延时补偿与排错4.1 采样周期 Ts 和时间尺度约束FCS-MPC 的采样周期直接决定了控制性能和计算负载Ts 越小电流纹波和总谐波失真越低但每个控制周期允许的运算时间越短。工程上有个经验范围电流环 MPC 的 Ts 通常取 50 微秒到 20 微秒对应 20 kHz 到 50 kHz 的控制频率这个量级既能保证电流环带宽又不至于让数字信号处理器算不过来。在 Matlab 仿真中Ts 还必须和电机模型的离散方式一致。如果使用 Simscape 模块求解器需要支持连续状态但 MPC 控制器只在每个 Ts 时刻更新一次输出。纯脚本仿真里被控对象如果也采用前向欧拉那么 Ts 过大时会导致预测发散一个直观的判据是Ts 2 * Ld / Rs。用上表中的参数算一下2 * 0.00525 / 0.958 ≈ 0.011秒也就是 11 毫秒1e-5 秒的 Ts 离发散边界很远但如果把电感改小到 5e-4 H同样的 Ts 计算就会出现数值不稳定此时必须减小 Ts 或改用隐式欧拉离散。4.2 代价函数权重 λ 的整定边界第 2 章提到过抑制开关频率需要在代价函数里加动作次数惩罚项。加入开关频率惩罚后的代价函数如下% J_cost 里追加开关次数项 J(i) (id_ref - id_p)^2 (iq_ref - iq_p)^2 lambda_sw * n_switch(i);其中n_switch(i)表示从上一周期开关状态切换到第 i 个状态时发生动作的桥臂数量取值为 0 到 3。lambda_sw 的取值范围需要根据电流纹波的实际幅度来定而不是拍脑袋取 0.1 或者 1。这里给一组参考做法先让lambda_sw 0跑一次记录稳态时 id、iq 纹波的均方根值再把lambda_sw从 0 逐渐增大直到开关频率下降约 20%观察电流跟踪误差的均方根上升了多少。这个斜率就是你可接受的损失边界。超过这个边界代价函数里电流项权重相对变小iq 实际值会明显偏离给定值出现“静差”般的效果一般调参时把这个当成权重过大的信号。还有一种容易误用的做法是把 d、q 轴的电流误差项单独设权重。比如J λd * delta_id^2 λq * delta_iq^2通过让 λd 大于 λq 来优先保证 d 轴电流为 0。实际在表贴式电机里d 轴电流对转矩贡献很弱加大 λd 只会让最优矢量选择偏向 d 轴电压效果是 q 轴电流跟踪变慢转矩响应变迟钝所以建议表贴式电机里不要单独调两个电流权重统一用相等的权重约束效果通过限幅和过调制逻辑实现。4.3 一拍延时补偿预测从 k1 时刻开始数字控制器从采样到 PWM 更新有计算延时这个延时等于一个控制周期。如果 MPC 仍然只预测一拍实际作用到电机的开关状态晚了一拍等效于引入一个 Ts 的纯滞后会显著降低电流环相位裕度。标准补偿做法是用当前拍周期实际作用的电压矢量先预测出一个虚拟的 k1 时刻电流再从这个虚拟电流出发评估 8 个候选矢量。function [sa, sb, sc] mpc_pmsm_delay_comp(id, iq, id_ref, iq_ref, theta, we, ...) ud_prev ...; % 当前正在作用的电压矢量 dq 分量上一拍选择的矢量 [id1, iq1] pmsm_predict(id, iq, ud_prev, uq_prev, we, Ts, Rs, Ld, Lq, psif); for i 1:8 [id2, iq2] pmsm_predict(id1, iq1, ud(i), uq(i), we, Ts, Rs, Ld, Lq, psif); J(i) (id_ref - id2)^2 (iq_ref - iq2)^2; end % 后续选优逻辑和普通 MPC 一致 end注意两次预测都用同一个电角速度 we这是近似处理因为在 Ts 量级下转速变化很小可以忽略。如果电机处于极低速或者零速启动阶段角度变化率低这个近似没有问题高速时 we 本身较大但 Ts 内的增量依然很小。The key is the delay compensation has to be symmetric偶数拍每条路径都补偿一次否则补偿一半会导致预测起点错位。4.4 仿真里最常见的四类异常现象可能原因排查与处理NaN 或电流发散Ts 超过数值稳定边界降低 Ts或检查电感和电阻参数数量级代数环错误MATLAB Function 输出直连到了自身输入插入 Unit Delay模拟一拍延时电流波形呈现等宽脉冲Park 变换角度方向不匹配用 theta0 的固定角度测试电压变换方向每次开关动作都切到同一个矢量代价函数权重失衡先令 lambda_sw0 验证跟踪性能这些异常里NaN 数据可以打开 Simulink 的数据记录器定位第一次出现 NaN 的时刻往回推几个周期检查输入信号是否越界。特别要注意的是 MATLAB Function 块内部如果有全局变量或工作区共享变量代码生成后会静默失效信号表现为上电后第一次调用时输出为 0之后全部正常。这种情况经常在使用persistent变量保存上一拍开关状态时出现解决思路是把所有状态都建模为输入端口用 Unit Delay 模块保存状态。5. 从电流环扩展到转速环闭环和双矢量 MPC5.1 转速外环 PI 与内环 MPC 的级联结构电流环 MPC 只解决“给定 id_ref 和 iq_ref如何选电压矢量”的问题转速指令需要由上层的转速环生成。工程上最可靠的结构是外环用整数阶 PI内环用 MPC这样既发挥 MPC 在电流快速跟踪上的优势又保留 PI 转速环的工程成熟度。转速环的采样时间一般是电流环的 5 到 10 倍比如 Ts_speed 1e-4 秒而电流环仍为 1e-5 秒。转速环程序每 10 个电流周期更新一次 iq_ref离散 PI 实现如下%% 转速外环 PI每 Ts_speed 执行一次 err_w w_ref - w_mech; iq_ref Kp_w * err_w i_w_integral; i_w_integral i_w_integral Ki_w * err_w * Ts_speed; % 转矩电流限幅 iq_ref min(max(iq_ref, -iq_max), iq_max);注意积分项必须在限幅之后处理一种常用做法是用当前实际输出 iq_ref 反推积分项防止积分饱和导致转速超调后难以回落。MPC 的内环响应速度很快外环 PI 的带宽反而不能设太高否则两级控制器的相位叠加会造成转速抖动。经验值是转速环带宽设为电流环带宽的 1/10 到 1/20也就是电流环按 1 kHz 设计时转速环只到 50100 Hz。如果想把转速环也换成 MPC代价函数会加入转矩跟踪项并需要考虑转速预测模型复杂度明显上升但仿真的整体结构和电流环 MPC 完全一致预测、选优、输出。这种模型预测速度控制的优势是能直接处理转速限制和电流限制的耦合约束代价是每个控制周期要多算一组转速预测方程。5.2 双矢量 MPC用占空比换取更低的电流纹波单矢量 MPC 每个采样周期只输出一个有效电压矢量或零矢量电流纹波相对较大。双矢量 MPC 在一个周期内分配两个矢量的作用时间最常用的是一个最优有效矢量加一个零矢量通过调整两者占空比来逼近目标电压。占空比计算的近似公式是% J1 为最优矢量作用整个周期时的代价 % J0 为零矢量作用整个周期时的代价 d J0 / (J1 J0); % 最优矢量的作用时间比例 % 实际输出先施加 d*Ts 的 v_opt再施加 (1-d)*Ts 的零矢量这个公式的物理解释很直观如果零矢量的代价非常小说明当前电流误差靠近目标应该多用零矢量让电流自然续流反之如果最优矢量的代价小到接近 0说明必须全力加速跟踪此时 d 接近 1。实际仿真里发现d 的表达式分母加一个极小数 ε 可以避免 J1 和 J0 同时为 0 时的除零错误这个情况通常发生在电流已经精确到达参考值时。5.3 验证比对的三个具体指标把单矢量 MPC 和双矢量 MPC 在同一模型上跑完 0.1 秒稳态仿真后可以用三个指标量化差异iq 电流纹波的峰峰值、相电流的总谐波失真、以及平均开关频率。Matlab 里可以用mean(abs(diff(sa))) / Ts计算平均开关频率用thd函数直接分析相电流的谐波成分。双矢量 MPC 通常能把 iq 纹波降低 30% 左右但代价函数计算量翻倍。最后给一个调试建议跑通电流环 MPC 之后先把id_ref设成 0iq_ref给一个 2 A 左右的阶跃观察响应曲线。理想的 MPC 电流响应应该在一到两个采样周期内跟上给定不会有超调这比调 PI 的时间常数直观得多。确认这一步静态跟踪没有静差后再叠加转速外环。如果这一步出现电流静差优先检查预测模型里电感参数是否和电机对象一致不匹配时预测误差会表现为固定的稳态误差这和 PI 控制器的积分作用逻辑完全不同。本文还有配套的精品资源点击获取