基于MATLAB的PMSM abc三相系统仿真:从状态方程到S-Function实现
2026/9/13 1:55:14 网站建设 项目流程

简介:一套基于MATLAB的电机ABC系统仿真程序,面向电气工程、自动化及电机控制领域的学习者与工程师,用于模拟交流电机在三相坐标系下的运行特性,帮助理解速度控制、扭矩输出与效率分析。压缩包共4个M文件,包含主仿真脚本、ODE求解与同步发电机模型等,整体仅4KB,代码紧凑,便于直接修改与二次开发。目前已有186人浏览学习,适合初次接触电机建模仿真的研究生、工程师及课程设计使用者。资源内的仿真程序覆盖三相电压源、定子绕组、转子绕组与负载等典型环节,支持设定三相参数及控制器策略,运行后可得到电流波形、转速、转矩与效率等关键结果,方便进行启动性能评估、稳态分析、过载能力验证和故障诊断。可在Simulink中进一步结合电压矢量控制或直接转矩控制算法优化设计,并借助MATLAB代码生成部署到嵌入式实时环境,缩短电机驱动方案开发周期。

1. 为什么电机仿真要单开一套 abc 系统而不是只套 dq 模型

在电机控制从业者的工具箱里,dq 坐标是分析的主力,但真正做逆变器接入、死区补偿、不平衡电网电压这类问题的工程师,往往更需要一份能直接吃三相电压的 abc 仿真程序。所谓“matlab_电机abc系统仿真程序”,就是在 MATLAB/Simulink 里把永磁同步电机(PMSM)或异步电机的三相自然坐标系模型写成独立程序,输入是 ua/ub/uc,输出是 ia/ib/ic、转矩和转子位置。它不替你去掉坐标变换,而是让你在坐标变换之外,多一个能验证相电压与相电流一致性的参照物。适合两类人:一类是刚学 FOC 想搞清楚电流波形从哪来的读者;另一类是已经在用 dq 控制,却想排查零序、死区、负载突变影响的工程师。仿真程序的骨架并不复杂:四个状态、三路电压输入加一个机械方程,真正决定结果质量的是参数换算和中性点处理方式,下文按这个顺序展开。

2. 在 MATLAB 里把 PMSM 的 abc 电压磁链方程整理成可积分的状态方程

2.1 表贴式永磁电机的相电压方程为什么能写成每相独立形式

做“电机abc系统仿真程序”之前,第一步是把物理方程落到状态方程。这里以表贴式永磁同步电机为例,忽略磁路饱和、涡流损耗和温度变化,定子三相绕组的电压方程写为:

u_a = R_s * i_a + d(psi_a)/dt psi_a = L_s * i_a + M * i_b + M * i_c + psi_f * cos(theta_e)

b、c 两相把角度依次平移 2π/3。由于电机中性点悬空,三相电流天然满足 ia+ib+ic=0,互感项 M(ib+ic) = -M·ia 被吸收进“一相等效电感” Lph = L_s - M。于是每相电流独立,右端函数里不再需要同时求解另外两相的耦合项。反电动势由永磁磁链对时间求导得到,θe 随时间变化时,ea = -ψf·ωe·sin(θe),其中 ωe = dθe/dt = p·ωm。

这里要注意两点。第一,这个“独立”结论依赖星形连接且无零序通路,如果电机改成四线制或直接把中性点接地,仿真程序就必须回到完整的 3×3 电感矩阵。第二,方程里没有出现凸极项,是因为表贴式转子的自感不随转子位置变化;内置式 PMSM 的 Ld 与 Lq 不等,abc 模型的自感和互感都会随 θe 变化,不能继续用常数 Lph,这点在第 5 章会单独讲。表贴式模型用到的主要符号如下表。

符号物理含义常用单位
R_s相电阻Ω
L_s相自感H
M相间互感,代数值为负H
Lph一相等效电感 L_s - MH
ψf永磁磁链幅值Wb
θe转子电角度rad
ωe转子电角速度rad/s
p极对数1

2.2 机械方程、电角度换算与 4 状态变量选择

机械侧的方程是转动惯量乘角加速度等于电磁转矩减负载转矩和阻尼转矩,转子位置角再对机械角速度积分。电磁转矩在 abc 坐标下写成三相电流与反电动势系数直接相乘的形式:

Te = -p * psi_f * [ia*sin(θe) + ib*sin(θe - 2π/3) + ic*sin(θe + 2π/3)]

这个式子看起来和 dq 下的 Te = 1.5·p·ψf·iq 差别很大,做一次 Park 变换就能互相推出。把 θe = p·θm 代入后,状态方程里的非线性项只剩 sin/cos,适合用 ode45 直接处理。状态变量一般选 x = [ia; ib; ωm; θm] 四个,而不是把三个电流都放进去。原因很直接:ic 被电流约束锁死,它在代数上是 ia、ib 的线性组合;如果强行积分三个电流,初始值和每一步的导数都必须满足 ia+ib+ic=0,ode45 在误差估计里会把这个约束当成额外自由度,既拖慢步长又容易在事件触发时出现毫安级的零序残差。实际工程里把 ic 作为输出计算即可。

至于角度状态,积分 θm 而不是 θe,是因为机械方程本身需要机械角速度,而 θe 在 θm 乘 p 之后可以直接得到。同时把 θe 的数值范围控制在 2π 以内,后续如果做 Park/Clarke 变换,角度查表时不会出现大角度精度问题。另外还要在每次采样时对 θe 做 mod(θe, 2π) 归一化,避免长时间仿真后角度值过大影响正弦计算的数值一致性。

2.3 把电机手册参数换算成仿真程序输入的三张表

仿真程序吃进去的是 R_s、Lph、ψf、J、B、p 这六个量,电机手册通常不会直接给出全部。我一般按下面这张表做映射:

程序参数手册来源说明
R_s绕组相电阻测试时注意环境温度
Lph表贴式取 Ld换算规则见 5.1 节
ψf反电动势系数换算换算表见 5.2 节
J转子惯量或转动惯量带负载时要把负载折算进去
B阻尼系数小功率电机可先置 0
p极对数电机极数除以 2

换算中最容易出问题的是 ψf。一个常见误区是拿线电压有效值直接除以电角速度,结果偏大 √3 倍。正确做法是先明确手册给的是相反电动势峰值、线电压有效值还是线电压峰值,三者之间系数分别是 √2、√3 的组合。这一节的数值规则放到第 5 章集中展开,因为最后一个换算式里还要用到极对数 p 和额定转速。如果手头只有 dq 坐标系下的转矩常数 Kt,也可以直接用 ψf = 2·Kt/(3p) 反算,前提是 Kt 的定义是峰值转矩与峰值 q 轴电流的比值。

这段代码是参数初始化的最小模板,后面所有脚本都从这里改:

prm.Rs = 0.5; % 相电阻,欧姆 prm.Lph = 5e-3; % 一相等效电感,亨利(表贴式电机取 Ld) prm.psi_f = 0.11; % 永磁磁链,韦伯 prm.J = 0.002; % 转动惯量,kg*m^2 prm.B = 0; % 阻尼系数 prm.p = 4; % 极对数

prm 是 MATLAB 结构体,用 prm 前缀而不是一堆全局变量的好处是:ode45 的右端函数只需要多接一个参数,后续做参数扫描时可以循环修改结构体字段而不影响函数签名。注意 Lph 的单位是亨利而不是毫亨,读手册时 5 mH 要写成 5e-3,这类单位错误在仿真里表现不是报错,而是电流整体变大 1000 倍,极难排查。

3. 用 MATLAB 脚本把 abc 仿真程序跑起来:ode45 右端函数与启动脚本

3.1 最小可运行的右端函数 pmsm_abc_rhs.m

基于第 2 章的状态方程,先写一个能被 ode45 直接调用的右端函数。这个函数只做四件事:取状态、算反电动势、削掉中性点电位、返回导数。

function dx = pmsm_abc_rhs(t, x, prm, ufun) % 状态: x = [ia; ib; wm; theta_m] % 输出: dx 为对应导数 ia = x(1); ib = x(2); wm = x(3); theta_m = x(4); ic = -ia - ib; % 由约束得到第三相电流 theta_e = prm.p * theta_m; we = prm.p * wm; u = ufun(t, prm); % 三相端电压,列向量 [ua0;ub0;uc0] u_n = (u(1)+u(2)+u(3)) / 3; % 悬空中性点漂移电位 ua = u(1) - u_n; ub = u(2) - u_n; uc = u(3) - u_n; ea = -prm.psi_f * we * sin(theta_e); eb = -prm.psi_f * we * sin(theta_e - 2*pi/3); ec = -prm.psi_f * we * sin(theta_e + 2*pi/3); dia = (ua - prm.Rs*ia - ea) / prm.Lph; dib = (ub - prm.Rs*ib - eb) / prm.Lph; Te = -prm.p * prm.psi_f * ... (ia*sin(theta_e) + ib*sin(theta_e - 2*pi/3) + ic*sin(theta_e + 2*pi/3)); dwm = (Te - prm.TL - prm.B*wm) / prm.J; dtheta = wm; dx = [dia; dib; dwm; dtheta]; end

代码里出现了一个前面推导没重点提的分支:u_nufun返回的是三个端子相对“电源参考地”的电压,不是相对电机中性点的相电压。当三相端电压之和不为 0,比如死区时间、母线波动、逆变器上下管不一致都会导致这种情况,悬空的中性点会整体漂移到三者平均值。把参考电压减掉u_n之后,等效相电压之和才严格为 0,电流约束才不会被破坏。

3.1.1 为什么不让 ode45 看到 ic 和 ec

从微分角度看,只要 didt_a + didt_b + didt_c = 0,约束就能保持。把 ic 作为代数式输出,ec 只用于转矩计算,不进入电流导数,就自动保证了这一关系。如果反过来在右端函数里写出 dic/dt = (uc - R_s·ic - ec)/Lph,那么零序电压会毫无阻碍地进入积分,最终结果在毫秒级内出现明显的直流量。这个现象在一部分公开代码里经常能看到,写成 abc 模型却不处理中性点,仿真曲线看着像模像样,换成线电压驱动后在零点附近全是毛刺。下面这个 ufun 用于空载启动测试:

Vmag = 30; f0 = 20; ufun = @(t, prm) [Vmag*cos(2*pi*f0*t); Vmag*cos(2*pi*f0*t - 2*pi/3); Vmag*cos(2*pi*f0*t + 2*pi/3)];

Vmag 是相电压幅值,f0 是电气频率。这里没加斜坡,所以启动瞬间电流会有一到两个周期的冲击。正式跑之前建议把 Vmag 用 ramp 从 0 抬升,或者直接给个限幅器。

3.2 空载启动的驱动脚本 run_abc_sim.m

有了右端函数和驱动电压,主脚本只需要设置参数、调 ode45、画三条电流曲线:

x0 = [0; 0; 0; 0]; tspan = [0 0.5]; opts = odeset('RelTol',1e-3,'AbsTol',1e-5,'MaxStep',1e-4); [t, x] = ode45(@(t,x) pmsm_abc_rhs(t, x, prm, ufun), tspan, x0, opts); ia = x(:,1); ib = x(:,2); ic = -ia-ib; wm = x(:,3); theta_m = x(:,4); figure; plot(t*1000, ia, t*1000, ib, t*1000, ic); legend('ia','ib','ic'); xlabel('t/ms'); ylabel('相电流/A');

MaxStep设置了 1e-4 秒,也就是 100 微秒。对于 20 Hz 的基波和 5 mH 电感,这个步长足够保证 ode45 不会在电流峰值附近来回试探。如果把 MaxStep 改成 1e-2,代码仍能跑完,但波形上会出现明显的锯齿,那是数值误差而不是真实电流纹波。几个关键量的取值和含义如下表。

名称取值作用
Vmag30相电压幅值 V
f020电气频率 Hz
MaxStep1e-4最大积分步长 s
RelTol1e-3相对误差容差
AbsTol1e-5绝对误差容差

3.3 快速健康检查:把反电动势和电流画在同一张图上查相位

很多人在 abc 仿真里第一步就是看三条电流长什么样,这其实不够。电流是电压、电阻、电感、反电动势共同作用的结果,单独看电流只能判断“跑没跑”,不能判断“对不对”。一个更快的体检项目是把反电动势画出来对比相位:

win = t > 0.3 & t < 0.32; we_end = prm.p * wm(end); ea_plot = -prm.psi_f * we_end * sin(prm.p * theta_m(win)); plot(t(win)*1000, ia(win), t(win)*1000, ea_plot); legend('ia(A)', 'ea(V)');

正常空载且电压频率刚好等于转子转速对应的电气频率时,电流会衰减到接近 0 附近的小值,并与反电动势保持 90° 左右的相位差。如果电流幅值远大于 Vmag/(Rs + jωLph) 的估算值,先查参数单位;如果电流完全不衰减,先查 ufun 里的频率填的是电气频率还是机械频率,给定频率应等于 p 倍机械频率,只填机械转速对应的 Hz 会让反电动势跑在电压前面,电机始终无法进入同步。

4. 把 abc 仿真程序接入 Simulink:Level-2 S-Function 的通用写法

4.1 为什么不自带模块而要改用 S-Function

Simulink 自带的 Permanent Magnet Synchronous Motor 模块效率高,但它的三相接口是经过坐标变换后的 dq 电压,死区、反电动势谐波、中性点漂移这些信息都被内部封装吃掉了。想验证 abc 仿真程序,最可靠的做法是把它封装成 Level-2 MATLAB S-Function,让 Simulink 里的 Scope、To Workspace 直接接在这个块外面。这样做的好处是第 3 章调通的右端函数原样复用,不会出现“脚本里能跑、模型里跑不出”的两套代码。

4.2 setup 里的端口定义和连续状态声明

S-Function 模板长这样:

function pmsm_abc_sfun(block) setup(block); function setup(block) block.NumInputPorts = 1; block.NumOutputPorts = 1; block.NumContStates = 4; block.SetPreCompInpPortInfoToDynamic; block.InputPort(1).Dimensions = 3; block.InputPort(1).DirectFeedthrough = true; % Derivatives 用到输入 block.InputPort(1).SamplingMode = 'sample'; block.OutputPort(1).Dimensions = 5; block.OutputPort(1).SamplingMode = 'sample'; block.NumDialogPrms = 1; block.SampleTimes = [0 0]; % 连续系统 block.RegBlockMethod('InitializeConditions', @Init); block.RegBlockMethod('Outputs', @Output); block.RegBlockMethod('Derivatives', @Deriv); function Init(block) block.ContStates.Data(1:4) = zeros(4,1); function Output(block) ia = block.ContStates.Data(1); ib = block.ContStates.Data(2); wm = block.ContStates.Data(3); th = block.ContStates.Data(4); block.OutputPort(1).Data = [ia; ib; -ia-ib; wm; th]; function Deriv(block) prm = block.DialogPrm(1).Data; u = block.InputPort(1).Data; % 3x1 x = block.ContStates.Data; dx = pmsm_abc_rhs_sfun(x, u, prm); block.Derivatives.Data = dx;

pmsm_abc_rhs_sfun(x,u,prm)和上一节的pmsm_abc_rhs(t,x,prm,ufun)内容一样,区别只在于把计算u的那一行换成直接使用入参。文件名要和函数名一致,放在当前路径或 MATLAB 搜索路径下,模块对话框里填结构体变量名prm而不是大括号内容。DirectFeedthrough最好置 true,因为Derivatives回调里读取了输入口数据,置 false 在少数求解器配置下会警告。

提示:把第 3 章的 pmsm_abc_rhs 函数签名从 (t, x, prm, ufun) 改成 (x, u, prm),删掉u = ufun(t,prm)一行,其余不动。t 在右端函数里本来就没用到。

4.3 端口信号约定和三种电压源接法

为了不让后来接手的人对着端口发呆,信号顺序固定成下面这样:

S函数端口信号含义典型上游模块
输入1(1)ua0 端电压参考地三相平均电压/逆变器输出电压
输入1(2)ub0同上
输入1(3)uc0同上
输出1(1:3)ia/ib/icScope/To Workspace
输出1(4)ωm速度环反馈
输出1(5)θm位置环/编码器解算

三相电压来源有三种接法:查表方式、正弦平均值方式和 PWM 脉冲方式。前两种直接给连续电压值,仿真步长可以放到 1e-5 秒以上;第三种接 PWM 发生器时要小心,为了分辨 10 kHz 载波,连续积分步长会压到微秒级,0.1 秒的仿真可能要跑十几分钟。我一般先用平均电压调通逻辑,再单独起一个开关级模型验证死区波形。如果把 S-Function 和自带 PMSM 模块并联做对比,两个块的转子初始角度必须对齐。自带模块内部通常从 d 轴对准 A 相绕组开始计数,而 abc 模型里 θm=0 对应的是永磁磁链最大值的位置,两者如果直接比较,电流相位会差一个电角度偏移。处理办法是在对比前做一次角度偏置校正,或者先把两路电流都变换到 dq 坐标再看。

5. 参数标定里最容易错的四个点:电感、Ke、零序和数值刚性

5.1 Lph 从 Ld/Lq 换算的“三分之二”之争

第 2 章从物理侧得到一相等效电感 Lph = L_s - M。到了工程手册上,厂商通常只给 Ld、Lq。对表贴式电机,Ld = Lq,很多教材又会把 dq 方程中的电感写成 Ld,读者在换算时经常纠结要不要乘 2/3。先说结论:如果仿真程序里的坐标系和 FOC 控制里用的坐标变换一致,且 Ld 是按“同步电感”标定的,注意 MATLAB 自带的 PMSM 模块输入对话框也是这样,那么在 abc 方程里直接把 Lph 填成 Ld 即可,不需要额外乘系数。原因在于 Clarke 变换有两种定义,功率不变形式会把 αβ 轴电感变成物理相电感的 1.5 倍,而大部分工程模型内部做的是去除零序分量的等幅值变换,两种定义在 dq 轴看到的电压方程最终会回到同一个数值,但 abc 方程里必须与你的变换约定保持一致。

怎么验证自己填对了:在 dq 模型里做一个电流环阶跃,看 d 轴电流时间常数 τd = Ld/Rs;再回到 abc 模型,给相同的相电压阶跃,测 ia 的上升时间。两个时间常数一致,Lph 就对了。

5.2 反电动势系数 Ke 的四种写法与 ψf 换算

永磁磁链 ψf 在 abc 方程里直接决定反电动势幅值,而手册里的“反电动势系数”单位极其混乱。我整理了一张换算表:

手册给的形式含义ψf 换算
E_ph_pk / ωe相反电动势峰值除以电角速度ψf = E_ph_pk / ωe
E_line_rms / ωe线电压有效值除以电角速度ψf = sqrt(2/3) * E_line_rms / ωe
Ke_Vpk_krpm线电压峰值 @ 1000 rpmψf = Ke / (sqrt(3) * (100π/3) * p)
Kt 峰值(N·m/A)峰值转矩常数ψf = 2Kt / (3p)

第二行那个 sqrt(2/3) 很反直觉:线电压 RMS 要先变相电压 RMS,除以 √3,再变相电压峰值,乘 √2,合并起来就是 √(2/3)。第四行来源于 dq 转矩公式 Te = 1.5·p·ψf·iq,这里的 Kt 是“峰值转矩/峰值 q 轴电流”,如果手册给的是 RMS 值,还得再乘 √2。这些都是老生常谈,但真到借别人的参数表时,依然有相当一部分写错。

5.3 悬空中性点电压 u_n:取平均不是玄学而是代数约束

这一节对应第 3 章代码里的u_n = (u(1)+u(2)+u(3))/3。为什么不能省?以三相半桥逆变器为例,当一个开关周期内 A 相上管占空比 50%、B 相 51%、C 相 49% 时,三个端电压的平均值不为 0,电机中性点会整体偏移。仿真程序里如果不削掉这个偏移,积分器会把零序电压直接加到相电流上,表现为三相电流在直流偏置附近同时抖动。这个偏置在有死区补偿的模型里尤其致命。如果电机中性点实际接地,则 u_n 强制为 0,代码相应改成:

if prm.isGrounded u_n = 0; else u_n = (u(1)+u(2)+u(3)) / 3; end

多一个isGrounded开关不会让计算量增加,但能让程序同时覆盖星形悬空和四线制两种拓扑。默认值设 false,因为绝大多数工业驱动的电机都是星形且中性点不引出。

提示:运行结果里 ia+ib+ic 的残差如果长时间大于 1e-6 A,先不要查积分器,第一嫌疑就是 u_n 被注释掉了。

5.4 步长、容差和数值刚性的边界

当输入是 PWM 脉冲电压时,右端函数里的 u 每微秒都在跳变,ode45 的误差估计器会不断缩步长,最终陷入“每步都失败重算”的循环。这时把所有责任推给求解器没用,要从模型层面压缩高频分量。三个可用方法,从最常用到最狠:

第一,用odesetMaxStep设置为开关周期的 1/20 到 1/50,并同时放大RelTolRelTol从默认 1e-3 放宽到 3e-3 通常能快一倍:

opts = odeset('RelTol',3e-3,'AbsTol',1e-5,'MaxStep',1e-4);

第二,改用求解器ode15sode23t,它们对“快电流慢机械”的多时间尺度问题更稳定,但代价是每次迭代都要算 Jacobian,PMSM 这个 4 阶系统完全承担得起。第三,如果只是验证控制逻辑,用平均电压模型。PWM 的开关纹波不属于控制设计关心的频段,把占空比乘以母线电压得到平均相电压,再把平均电压接到 S-Function 的输入口,仿真时间直接从“分钟级”降到“秒级”。需要看纹波时,再单独开一个开关级模型,两条支路共用同一个参数结构体。

6. 一个验证技巧:用复矢量轨迹和零序分量检查 abc 仿真程序对不对

6.1 把 ia/ib/ic 合成为复矢量看轨迹圆不圆

abc 仿真的输出最容易骗人:三条曲线看起来都是正弦,但幅值、相位、直流偏置只要错一个,整体依然是正弦。判断参数对不对,我一般不看三条线,而是算空间矢量:iα = ia,iβ = (ia + 2·ib)/√3,然后把 (iα, iβ) 画在复平面上,MATLAB 里做:

figure; plot(i_alpha, i_beta); axis equal; grid on;

稳态且参数正确时看到的是一个以原点为中心的圆;椭圆代表三相幅值或相位不对称;圆上有明显尖刺代表反电动势谐波或数值步长过大;圆心不在原点说明有直流偏置,优先查中性点处理。

6.2 零序分量 i0 的大小是程序健康的体温计

零序电流的定义是 (ia+ib+ic)/3,理论上无论星形还是四线制,这个量在相电流里都应该是 0。直接在脚本里加一行:

i0 = (ia+ib+ic)/3; rms_i0 = rms(i0); fprintf('零序电流 RMS = %.3e A\n', rms_i0);

这个数值如果在 1e-6 量级,说明状态选取和中性点处理正确;如果到了 1e-2 量级,说明右端函数里用了多余的第三相导数,或者 u_n 被去掉。这个方法同样适合检查从外部录波器导入的数据,先用 readmatrix 读入 CSV 再算同一表达式,能快速分辨现场数据和仿真程序的差异。

6.3 用 FFT 核对反电动势相位,揪出坐标变换的符号错误

最后补充一个和坐标变换联调的小技巧:在 abc 仿真程序里给永磁磁链的初始角度偏置 30°,接着把反电动势 ea 和编码器位置同时导出,做一次简短的 FFT,检查 ea 的相位是否与 sin(θe) 对齐。之前遇到过把 θe 初始位置设成 90° 导致 dq 变换后 id、iq 互换的案例,靠这条检查一眼就能看出来。数字上 ea = -ψf·ωe·sin(θe) 中的负号不要丢,丢了反电动势会变成励磁方向,电流被拉偏,控制器的 q 轴电流指令会变成奇怪的直流偏置。这个技巧不复杂,但它能在一分钟内回答“程序里的 abc 模型和控制器的坐标变换是不是同一个世界”。如果再把 iα、iβ 送进锁相环或 Park 变换,得到的 dq 分量与理想值之差,也可以反过来作为模型标定的误差指标。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询