尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

IPMSM标定核心:Ld/Lq曲面拟合与MTPA/MTPV轨迹求解

IPMSM标定核心:Ld/Lq曲面拟合与MTPA/MTPV轨迹求解 简介在新能源汽车驱动系统中内置式永磁同步电机IPMSM凭借高功率密度与宽调速范围被广泛采用而MTPA与MTPV控制是发挥其最高效运行能力的关键。这份资料聚焦IPMSM参数标定场景面向具备永磁电机数学方程基础的电机控制工程师提供一套完整的MATLAB标定算法用于解决不同电流下电机参数辨识及MAP表生成问题。算法基于已知电流下的Ld、Lq值通过拟合得到未知参数或利用台架标定数据快速生成MTPA和MTPV对应的id、iq MAP表从而降低手动标定工作量已在相关电机产品上得到实际应用。压缩包共15个文件其中11个m脚本覆盖LdLq拟合、MTPA/MTPV计算、MAP表获取、转矩磁链拟合等核心流程另含2个txt操作说明和2个xlsx结果数据表便于对照执行。整体仅144KB轻量易部署。目前已有64人参与学习适合需要理解IPMSM高效控制标定流程的工程师参考有助于系统掌握从参数辨识到控制策略落地的完整链路。1. IPMSM标定绕不开的核心矛盾Ld/Lq不是一个常数台架标定时常会遇到这种情况同一台电机空载点电流角算得挺准加载到额定转矩附近MTPA预测的电流角比实测小三四度转矩误差超过5%。问题不在公式(Ld-Lq)项就是磁阻转矩来源但Ld和Lq本身会随电流变磁路饱和让Ld随去磁电流增大明显下降交叉耦合又让iq反过来影响Ld。固定Ld/Lq的模型只在空载附近成立。真正能上车的标定得先把不同(id, iq)下的Ld/Lq拟合成曲面再在这个变参数模型上解MTPA和MTPV轨迹。这套MATLAB算法做的就是这件事入口顺序在matlab指令说明.txt里写清楚了。适合正在做eMotor台架标定或需要自己搭IPMSM控制模型的控制工程师和研究生。2. FitLdLq.m把散点Ld/Lq数据拟合成可用的连续模型2.1 先理解为什么不能拿额定点的Ld/Lq直接算IPMSM的转矩方程是 T 1.5p[ψm·iq (Ld-Lq)·id·iq]磁阻转矩完全由(Ld-Lq)这一项贡献。很多人在仿真阶段用供应商给的额定点Ld/Lq跑出来效率很高一上台架就露馅。原因在于主磁路饱和后磁导率下降Ld随去磁电流(id为负)显著减小q轴电流增大时Lq同样下降而且下降斜率往往比Ld还大。结果是(Ld-Lq)这个差值在全电流范围内可能是额定点的0.8倍到1.5倍拿固定值算出来的MTPA电流角在深度饱和区域会整体偏小这才出现开头说的三四度偏差。所以标定第一步不是上来就解MTPA而是先把Ld/Lq表达成关于id和iq的连续函数。算法包里的FitLdLq.m做的是这件事它读入你已经准备好的电流-电感离散表输出拟合函数对象供后面MTPA1.m、GetMtpvIdIq.m反复调用。输入数据来源可以是台架也可以是FEA仿真算法本身不区分区分的是覆盖范围够不够。2.2 输入数据的组织和坐标系约定要跑通这个环节先准备一份四列数据id、iq、Ld、Lq。单位统一用A和mH脚本内部不做单位换算。包里的results.xlsx就是模板典型sheet结构如下列符号单位典型范围说明AidA-Imax ~ 0d轴去磁电流负值BiqA0 ~ Imaxq轴转矩电流正值CLdmH0.03 ~ 0.25随-id增大而减小DLqmH0.08 ~ 0.45随iq增大而减小扫描步长我一般取Imax/20乘用车主驱Imax在300A到600A之间步长15A到30A这个量级就够了。步长太细台架数据相邻点的电感差异会被测量噪声淹没步长太粗MTPV轨迹到了高转速区会抖动。数据准备好后Loaddata.m负责把xlsx读成matlab表格变量再交给FitLdLq.m。2.3 拟合核心代码与阶数选择FitLdLq.m内部做的事情可以用下面这段代码概括实际脚本里多加了边界判定但核心逻辑一致% 读入四列原始数据Ld和Lq分开拟合 data readtable(results.xlsx, Sheet, LdLq_raw); id data.id; iq data.iq; Ld data.Ld; Lq data.Lq; % poly22 形式: p00 p10*x p01*y p20*x^2 p11*x*y p02*y^2 % 对车用 IPMSM 来说二阶已经能抓住 Ld/Lq 随电流变化的主要趋势 fit_ld fit([id, iq], Ld, poly22); fit_lq fit([id, iq], Lq, poly22);说明fit的第一个输入是[id, iq]两列矩阵第二个输入是被拟合的Ld列向量。poly22共6个系数对车用IPMSM的Ld/Lq曲面足够。拟合完成后任意(id, iq)点都能即时算出Ld/Lq这就是后续求导和优化能够成立的前提。poly33、poly44不是不能用但台架数据本身带噪声时高阶多项式会把测量噪声当成真实趋势拟合进去残差曲线出现波浪形反而不如poly22稳。如果你手里的数据是FEA仿真的理想平滑曲面再考虑poly33。2.4 拟合质量检查别只看R²脚本末尾会输出拟合统计量。我一般看两个指标R²和最大绝对残差。R²做到0.99以上不稀奇因为Ld/Lq曲面整体趋势很平滑真正要盯的是最大残差出现的区域如果出现在大电流边界说明边界点被当成了离群点剔除。常见做法是保留边界值并在FitLdLq.m末尾做一步饱和外插保护当查表电流超出拟合范围时电感取边界值而不是继续用多项式外推。多项式外推在超出范围100A以上时会出现非物理的负电感这个坑在MTPV高转速区特别容易踩到。提示不需要Simulink和额外工具箱脚本只用到fit和readtable这些基础函数R2019b之后任何matlab安装版本都能直接跑包里的matlab指令说明.txt里有逐行解释。新版matlab对xlsread有弃用警告统一用readtable就不会有兼容问题。3. MTPA1.m与MTPA2.m两种思路解最大转矩电流比轨迹3.1 MTPA的数学条件MTPA的本质是在电流幅值约束 id²iq²Is² 下最大化 T(id,iq) 1.5p[ψm·iq(Ld-Lq)·id·iq]。用拉格朗日乘子法可以得到最优性条件在最优轨迹点上转矩梯度与约束圆的法向平行写成 ∂T/∂id·iq - ∂T/∂iq·id 0。注意这里的Ld、Lq已经是2.3节拟合出来的关于(id,iq)的连续函数求偏导时Ld/Lq本身的导数项不能丢这是固定参数MTPA公式经常漏掉的部分。有了这个条件问题变成一个求根问题对给定的请求转矩T*找到一组(id,iq)使转矩方程和MTPA条件同时成立。算法包提供了两种风格不同的解法MTPA1.m走解析数值求根MTPA2.m走网格扫描。两种结果理论上一致实际用哪个看你对数据噪声的容忍度。3.2 MTPA1.m解析求根路线MTPA1.m的思路是把电流写成幅值和角度的极坐标形式然后用fzero找满足转矩条件的角度。示意代码如下% 给定请求转矩 Treq求解对应的电流角 beta电角度 % 约定: id Is*cos(beta), iq Is*sin(beta), beta 范围 [0, pi/2] for k 1:length(Treq) beta_sol(k) fzero((b) ... TorqueByBeta(b, Is_guess(k)) - Treq(k), [0.01, pi/2-0.01]); id_mtpa(k) Is_guess(k) * cos(beta_sol(k)); iq_mtpa(k) Is_guess(k) * sin(beta_sol(k)); end参数说明TorqueByBeta是脚本内置函数内部调用FitLdLq得到的拟合对象把极坐标电流代回转矩公式。fzero的搜索区间给[0.01, pi/2-0.01]是为了避开pi/2端点那里iq接近零、转矩对角度不敏感数值上容易抖动。Is_guess是电流幅值的初始猜测从空载到峰值电流分段给定每段用上一次的解做初值收敛快很多。实际运行时你会发现beta初值给得太远fzero会走到非物理解上去。常见做法是先用固定Ld/Lq的解析公式算一个初值beta0再在这个初值附近展开搜索。这个细节在MTPA1.m里做了属于典型的工程性保护。3.3 MTPA2.m网格扫描路线如果交叉耦合严重解析求偏导表达式很长漏一项就难查。MTPA2.m换了个思路在d-q电流平面上直接铺网格用拟合好的Ld/Lq计算每个网格点的转矩然后对每个请求转矩T*找电流幅值最小的网格点作为MTPA工作点。% 在 [ -Imax, 0 ] x [ 0, Imax ] 上生成电流网格 id_grid linspace(-Imax, 0, 51); iq_grid linspace(0, Imax, 51); [ID, IQ] meshgrid(id_grid, iq_grid); % 用拟合对象批量算转矩矩阵 T_grid calcTorqueGrid(ID, IQ, fit_ld, fit_lq); IS_grid sqrt(ID.^2 IQ.^2); % 对每个转矩点取 IS 最小的网格点作为 MTPA 工作点 for k 1:length(Tstar) mask abs(T_grid - Tstar(k)) T_tol; [~, idx] min(IS_grid(mask)); id_map(k) ID(mask); iq_map(k) IQ(mask); end逻辑说明网格遍历覆盖了全部允许电流区域不需要担心fzero初值问题。代码里ID(mask)这种取法可能取到多个满足容差的点取幅值最小的那一个。实际脚本要控制转矩容差T_tol取请求转矩的0.5%比较合适太大会把相邻转矩等级的轨迹点混在一起。网格法的缺点是精度受步长限制。我在项目里一般用51×51网格对应Imax500A时步长10A得到的MTPA轨迹在中转矩区误差在1%以内如果map表用于最终量产刷写建议把网格加密到101×101计算量虽大一点但matlab矩阵运算几秒就能跑完。3.4 合成MTPA表并与转矩公式自检MTPA1和MTPA2跑完之后各自生成一组(id_mtpa, iq_mtpa, T)表。两张表直接画在同一张图里理想情况下曲线重叠如果出现局部偏差差异点大概率落在电流幅值接近Imax的区域此时以MTPA2的扫描结果为准因为解析法在边界处的fzero收敛半径会缩小。最后用TorqueFluxSpeedCmd.m把id/iq变换成磁链和转矩命令输出格式如下列名含义单位T_cmd请求转矩N·mid_cmdd轴电流命令Aiq_cmdq轴电流命令Aflux_cmd定子磁链命令WbU_ref参考电压幅值V这里有个自检技巧把生成的id_cmd、iq_cmd代回转矩公式算一遍和下发的T_cmd对比偏差应小于0.5%。如果偏差大优先检查拟合残差大的区域是不是刚好落在MTPA轨迹上。4. GetMtpvIdIq.m高转速区MTPV轨迹与电压极限的求解4.1 为什么MTPA标完之后还要标MTPVMTPA只保证电流幅值受限下的最高效率车速升高后逆变器输出电压受限MTPA点对应的电压幅值会超出调制极限必须把工作点往去磁方向压进入弱磁区。弱磁I区沿电流极限圆走弱磁II区沿MTPV轨迹走。很多标定脚本只做MTPA和弱磁I区不碰MTPV结果就是最高转速区扭矩输出不足或者电流角乱跳。GetMtpvIdIq.m解决的就是弱磁II区的轨迹计算。电压约束忽略电阻后写成(Lq·iq)² (Ld·idψm)² ≤ (Umax/ω)²。在id-iq平面上这是一个中心在(-ψm/Ld, 0)的椭圆转速升高时椭圆按1/ω比例缩小。MTPV的点就是恒转矩曲线族与该电压椭圆相切的切点集合。4.2 用fmincon把带约束的转矩最大化写成标准优化GetMtpvIdIq.m里用的是约束优化思路而不是像MTPA那样直接求切线方程。这在工程上更容易维护改电压约束或电流约束时只改约束函数就行。核心代码段% 输入: 转速 w, 母线电压 Udc, 峰值电流 Imax for k 1:length(w) Umax Udc / sqrt(3) - Udead; % Udead 为死区压降折算 lb [-Imax, 0]; ub [0, Imax]; % id 非正, iq 非负 nonlcon (x) voltageEllipse(x, Umax/w(k), fit_ld, fit_lq); x0 [id_mtpa(end), iq_mtpa(end)]; % 用上一转速的解做初值 opt optimoptions(fmincon, Algorithm, sqp, Display, off); x fmincon((x) -calcTorque(x(1), x(2), fit_ld, fit_lq), ... x0, [], [], [], [], lb, ub, nonlcon, opt); id_mtpv(k) x(1); iq_mtpv(k) x(2); end逻辑说明目标函数取负转矩fmincon最小化负转矩等价于最大化转矩。nonlcon返回电压不等式约束c(x) (Lq·x(2))² (Ld·x(1)ψm)² - (Umax/ω)² ≤ 0fmincon会通过增广拉格朗日法处理这个不等式。初值x0用上一转速的MTPV解因为相邻转速下最优工作点离得很近收敛快且稳定第一圈转速初值用的是MTPA末点正好从弱磁I区切进II区。这里有一个必须注意的点死区压降Udead的处理。我一般取Udc的3%到5%全桥逆变器的死区时间、管压降都会实时吃掉一部分电压Umax取理论值Udc/√3的话高速区计算出来的可用转矩会高估3%以上台架上一对比就露馅。如果整车控制用了过调制Umax上限可以适当放宽但电压余量要给足。4.3 如何从优化结果里区分MTPV段和电流极限段过弱磁拐点后并不是所有点都落在MTPV轨迹上。判断方法很简单看fmincon求得的解是否碰到电流幅值边界Imax。如果解的电流幅值小于Imax说明这是MTPV切点如果刚好等于Imax说明这是电流极限圆与电压椭圆的交点属于弱磁I区。GetMtpvIdIq.m输出的表里同时包含这两段作图标出来分界线一目了然。工程上拐点速度附近的切换区要做一次低通滤波防止查表时id在切换点附近跳变。4.4 MaxTMtpvWcontourFit.m把离散点拟合成连续命令表MTPV按转速逐点计算后得到若干组(n, id, iq)离散点。MaxTMtpvWcontourFit.m的作用是把这些离散点拟合成n-T平面上连续可查的id(n,T)、iq(n,T)命令面。脚本内部用的还是fit家族函数但x、y换成了转速n和转矩T。拟合前先把MTPV段和弱磁I段的接缝切开分别拟合接缝处做重叠平滑否则中间会出现马鞍形畸变。转速网格的取法有一点讲究低转速区电压椭圆很大相邻转速下id/iq变化小网格可以粗高转速区电压椭圆急剧收缩网格要加密。我通常用非线性间距转速区间网格步长说明0 ~ 3000 rpm100 rpm电压椭圆大工作点变化慢3000 rpm 以上25 rpm椭圆快速收缩需要加密配合起来matlab跑完生成results.xlsx的最终sheet再回填到整车控制器查表。5. DebugMFPT.m一张等高线图验证整套标定结果5.1 把转矩等高线、电压椭圆和MTPA/MTPV轨迹叠加画出来DebugMFPT.m是整个包里最值得先跑的脚本。它把前面几步的结果统一画在id-iq平面上转矩等高线、恒定转速下的电压椭圆、MTPA轨迹、MTPV轨迹以及电流极限圆。画法如下figure; contourf(ID, IQ, T_grid, 30); hold on; plot(id_mtpa, iq_mtpa, r-, LineWidth, 2); plot(id_mtpv, iq_mtpv, b--, LineWidth, 2); % 画电压椭圆: 中心(-psi/ld, 0)半轴 Umax/(w*lq) 和 Umax/(w*ld) phi linspace(0, 2*pi, 200); plot(-psi/ld (Umax/w)*cos(phi)/lq, (Umax/w)*sin(phi)/ld, k--); axis equal; xlabel(id (A)); ylabel(iq (A));验证逻辑很简单正确标定出的MTPA/MTPV轨迹应当和转矩等高线族处处相切。拿固定转速的电压椭圆去比红色线段伸出椭圆之外的那一段就是这个转速下的不可用点输出前需要截断。5.2 三个常见的坑坑一是Ld/Lq拟合范围不够。有的工程师在台架上只扫到id-200A实际最大去磁电流是-300A拟合曲面在大电流区纯靠外推。表现是Debug图里高转矩区域的MTPA轨迹明显偏离等高线切点。解决方式是回台架补点或者用FEA结果补充大电流边界。坑二是MTPV是否存在被误判。当电压椭圆中心到原点的距离ψm/Ld大于椭圆半短轴Umax/(ω·Ld)时电压椭圆不包含原点MTPV才有实际意义否则弱磁II区就是电流极限圆直接和电压椭圆相交整个弱磁段都该走电流极限。没做这个判断的脚本会在部分转速段错误地输出MTPV点电流角乱跳就是这么来的。坑三是map表分辨率不够导致插值振荡。整车控制器查表用线性插值标定表在弱磁区步长太粗时插值出的id会出现±5A以上的抖动Debug图上蓝线出现锯齿。处理办法是回到MTPA2加密网格或对生成的map做一次savgol滤波滤波窗口取5个点、阶数2就行窗口太大会把真实的弱磁拐点抹平。把这三个问题排掉DebugMFPT.m画出来的四层曲线就能和台架实测逐点对上也。产线标定数据刷写前拿这组图做验收比只看数值表直观得多。本文还有配套的精品资源点击获取
返回列表