最近在搞外转子开关磁阻电机ER-SRM的优化设计拿NSGA-II跑多目标过程中踩了不少坑也整理了不少好用的Matlab画图代码。这篇文章就把整个从问题建模到算法实现再到结果可视化的过程完整拆一遍重点解决两个事一是怎么把ER-SRM的结构参数和目标函数跟NSGA-II串起来二是你们最常问的Matlab画图代码到底怎么写才好看、能直接放进论文和报告里。不管你是刚接触电机优化的研究生还是已经在用有限元做设计的老工程师这份实践笔记应该都能帮你省下不少试错时间。1. 先搞清楚ER-SRM为什么要用NSGA-II优化1.1 外转子开关磁阻电机的特殊之处外转子开关磁阻电机External Rotor Switched Reluctance Motor简称ER-SRM最大的特点就是转子在定子外面定子在里面两者都是凸极结构依靠磁阻最小原理产生转矩。这个结构在直驱应用里特别吃香比如电动自行车轮毂电机、洗衣机直驱电机、风机类负载因为外转子可以直接和负载连接省掉了减速机构整个传动链刚度高、效率损失小。但外转子结构也带来一个麻烦转子转动惯量比内转子大转矩脉动更容易传导到负载端。我在实际测试中遇到过空载工况下电机振动明显噪声频谱里低频成分占主导。这种问题不是靠加几个PI参数就能解决的根子上还是电机电磁结构参数设计得不够匹配。比如定子极弧系数、转子极弧系数、气隙长度、绕组匝数这些参数稍微动一点转矩波形就会变化。传统做法是用有限元软件一个个参数扫描一个方案算一次瞬态场跑一组参数要一两个小时。想同时把转矩脉动、效率、材料成本都照顾到基本要几百组方案手工试错根本做不完。所以ER-SRM这种多参数多目标设计天然就是多目标优化算法的用武之地。1.2 多目标优化的矛盾点做电机优化的人都知道性能指标之间很少是共赢关系。拿ER-SRM来说你可能想提高平均转矩但为了达到高转矩往往要加大电流和极弧占比结果转矩脉动随之上升想降低铜耗提高效率就得把绕组匝数和电流密度往合理方向调整但气隙磁密又会变化转矩特性跟着变。我做过一个测试保持电流密度不变只把转子极弧系数从0.35调到0.45平均转矩升高了8%但转矩脉动系数从28%飙到了41%。这种此消彼长的关系在电机设计里几乎是无处不在的。所以问题就变成了在满足几何尺寸和热约束的前提下找到一组结构参数使多个目标同时尽可能好。这时候就需要帕累托最优的概念。多目标优化算法给出的不是一个解而是一组解每个解代表一个权衡方案。比如方案A效率高但转矩脉动大方案B转矩脉动小但效率略低两个方案在帕累托前沿上是非支配关系——没有哪个绝对优于另一个。工程师最后是挑一个符合实际工况的折中解。1.3 为什么选NSGA-II而不是其他算法NSGA-II全称是带精英策略的非支配排序遗传算法2002年由Deb提出。在多目标进化算法里它一直是应用最广泛、工程验证最充分的老大哥。我为什么没有用多目标粒子群或者SPEA2主要原因有三点。第一NSGA-II的处理逻辑非常透明非支配排序加拥挤度距离没有什么黑盒机制。一旦优化结果不理想很容易从种群分布和排序结果反推问题出在目标函数还是参数范围上。第二NSGA-II对决策变量的类型宽容度大连续变量和整数变量混合都没问题。ER-SRM里绕组匝数就是典型整数变量粒子群处理起来还得专门做藕合转换NSGA-II在遗传算子里处理起来更自然。第三Matlab里无论是内置的gamultiobj还是各种开源的NSGA-II工具箱实现都比较成熟代码可读性好。对工程师来说能用最短时间调通并做出可复现的结果这才是最关键的。2. 优化问题的定义与决策变量2.1 目标函数怎么定很多人拿到优化项目第一件事就是找算法但我建议先花时间把目标函数想清楚。ER-SRM优化最常用的三个目标是平均转矩最大化、转矩脉动最小化和效率最大化。注意NSGA-II内部默认是统一最小化所以最大化目标要取负号或者用倒数形式。我自己的Matlab目标函数接口长这样function f ER_SRM_eval(x) % x为结构参数向量 % f [f1, f2, f3] f1 -calc_average_torque(x); % 平均转矩取负因为要最小化 f2 calc_torque_ripple(x); % 转矩脉动越小越好 f3 -calc_efficiency(x); % 效率取负 end这样三个目标都是求最小NSGA-II直接处理不需要改算法底层代码。另外我建议加一个材料成本或者电机体积目标不过要注意目标数量太多时帕累托前沿的可视化会变难解的分布也不容易均匀。三个目标是我个人觉得工程应用里比较合适的折中。2.2 决策变量范围的选择决策变量的选择直接决定优化空间的形状。ER-SRM里影响电磁性能的结构参数不少但每个都加进去会导致维度灾难。我一般选下面这6个作为核心变量变量物理含义范围备注βs定子极弧系数0.35~0.55影响气隙磁导和转矩βr转子极弧系数0.30~0.50对转矩脉动影响大g气隙长度0.2~0.6 mm太小难加工太大转矩下降N每极绕组匝数20~60整数整数变量影响磁动势hys定子轭厚8~20 mm影响磁路饱和和重量hyr转子轭厚8~20 mm影响转子刚度和磁路变量范围不能拍脑袋。我建议先做一次多个参数的单因素扫描把每个变量的可行区间大致摸出来再放到优化里。否则NSGA-II的初始种群会散落在不可行区域算法前几十代都在找可行边界效率很差。2.3 约束条件处理约束条件是电机设计里绕不开的。常见的约束包括最大电流密度不超过5A/mm²、定子齿磁密不超过1.8T、外径限定在200mm以内、绕组温升不超过120℃。在NSGA-II里处理约束我优先用罚函数法。具体做法是如果某个体不满足约束在目标函数后面加一个和违反程度成正比的惩罚项让它很难在非支配排序中胜出。if violation eps penalty 1e6 * violation; f f penalty; end但这里有个细节罚函数系数不能太大否则会导致所有边界个体都变成支配解种群快速收敛到某个角落。我建议把惩罚项的量级设置成和正常目标值量级一致比如目标值都是两位数惩罚项基值取1000就足够了。还有一种更稳妥的做法是把约束违反程度单独作为一个目标函数但这样目标数量会多一个不太推荐除非你实在调不好罚权重。3. Matlab实现NSGA-II的完整流程3.1 NSGA-II算法核心流程回顾先快速过一遍NSGA-II的主循环因为后面代码框架都是围绕这个流程写的。算法流程是这样的初始化一个种群种群中每个个体对应一组ER-SRM结构参数。计算每个个体的目标函数值。对当前种群做快速非支配排序得到每个个体所处的帕累托层rank。计算同一层内个体的拥挤度距离用于区分优先保留哪些个体。通过锦标赛选择选出父代进行交叉和变异生成子代种群。将父代和子代合并重新做非支配排序和时间距离排序选出下一代种群。重复2~6直到达到最大代数。这个流程里最关键的是第3步和第6步因为非支配排序的质量决定了解集的收敛性拥挤度距离决定了解集的均匀分布性。3.2 代码框架与关键函数设计在Matlab里我习惯把算法拆成几个独立的函数主脚本、目标函数、非支配排序函数、锦标赛选择函数、交叉变异函数。这里给出一个常用的主循环骨架你可以直接在这个基础上去改% main_ER_SRM_NSGAII.m clc; clear; close all; n_var 6; n_obj 3; pop_size 80; max_gen 200; bounds [0.35 0.55; 0.30 0.50; 0.2e-3 0.6e-3; 20 60; 8e-3 20e-3; 8e-3 20e-3]; % 初始化种群 pop zeros(pop_size, n_var); for i 1:pop_size for j 1:n_var if j 4 % 匝数取整数 pop(i,j) round(bounds(j,1) rand*(bounds(j,2)-bounds(j,1))); else pop(i,j) bounds(j,1) rand*(bounds(j,2)-bounds(j,1)); end end end % 评估初始种群 obj zeros(pop_size, n_obj); for i 1:pop_size obj(i,:) ER_SRM_eval(pop(i,:)); end % 代数循环 for gen 1:max_gen % 非支配排序 [rank, dist] non_dominated_sort(obj); % 锦标赛选择 parent_idx tournament_selection(rank, dist, pop_size); % 模拟二进制交叉和多项式变异 offspring crossover_mutation(pop(parent_idx,:), bounds); % 评估子代 offspring_obj zeros(pop_size, n_obj); for i 1:pop_size offspring_obj(i,:) ER_SRM_eval(offspring(i,:)); end % 父子合并 combined_pop [pop; offspring]; combined_obj [obj; offspring_obj]; % 重插入选择 [pop, obj] select_next_generation(combined_pop, combined_obj, pop_size); % 记录每代最优值以目标1最小为例 best_hist(gen) min(obj(:,1)); end这个框架很通用你只需要把ER_SRM_eval替换成你自己的评估函数。非支配排序和时间距离函数是实现难点但我建议一开始先用经典的O(N^2)实现别急着写太复杂的高效版本。种群规模80一次排序也就是6400次比较Matlab完全顶得住。3.3 目标函数与有限元/解析模型结合的方式目标函数是ER-SRM优化里最耗时的部分。如果你打算在NSGA-II里直接调用Maxwell或者Flux的瞬态仿真脚本那你得做好等一整天的心理准备。我最初就是直接用Maxwell脚本算每一代一个个体算一次3秒转速下的瞬态场大概40秒种群80个个体跑200代总时间就是8020040/3600约178小时这谁顶得住。所以工程上必须用代理模型或者简化解析模型。最简单的是先离线采样几十个参数组合用Maxwell算出对应的平均转矩、转矩脉动和效率然后用Matlab的fitrgp做高斯过程回归或者用Kriging工具箱把采样点拟合成一个响应面。之后ER_SRM_eval内部直接调用这个代理模型几毫秒就能出结果。function f ER_SRM_eval(x) % 加载训练好的代理模型 gpr_torque, gpr_ripple, gpr_eff persistent gpr_torque gpr_ripple gpr_eff if isempty(gpr_torque) S load(srm_surrogate.mat); gpr_torque S.gpr_torque; gpr_ripple S.gpr_ripple; gpr_eff S.gpr_eff; end f(1) -predict(gpr_torque, x); f(2) predict(gpr_ripple, x); f(3) -predict(gpr_eff, x); end用代理模型优化完以后一定要记得对最终选出来的几个帕累托解再跑一次真实有限元验证。代理模型毕竟有误差误差小于3%一般可以接受如果发现某个候选解实际转矩脉动比预测高很多就要补充采样点重新训练模型然后重跑优化。4. 画图代码与结果可视化4.1 Pareto前沿绘制先说说二维Pareto前沿怎么画。这是最基础的图适合展示两个目标之间的关系比如转矩脉动和平均转矩。最后一代的非支配解集可以从obj变量里筛选出来然后画散点图。figure(Color, w, Position, [100 100 700 500]); plot(-obj(:,1), obj(:,2), o, MarkerSize, 6, ... MarkerEdgeColor, [0.2 0.4 0.8], MarkerFaceColor, [0.6 0.8 1]); xlabel(平均转矩 (N·m), FontSize, 12, FontWeight, bold); ylabel(转矩脉动 (N·m), FontSize, 12, FontWeight, bold); title(ER-SRM帕累托前沿转矩脉动 vs 平均转矩, FontSize, 14); grid on; box on;这里有个坑NNSGA-II里我初始是把平均转矩取负了所以画图时要记得再次取负不然横轴是负数看着就奇怪。另外如果你用内置gamultiobj函数返回的是求最小化的目标值也需要做转换。如果优化了三个目标建议用三维散点图figure(Color, w); scatter3(-obj(:,1), obj(:,2), -obj(:,3), 30, -obj(:,3), filled); xlabel(平均转矩); ylabel(转矩脉动); zlabel(效率); colormap(jet); colorbar;三维图最重要的是观察解的分布是否延展开。如果所有点挤成一团说明算法早熟或者目标函数区分度不够需要回去调种群大小或目标函数权重。4.2 收敛曲线绘制多目标优化的收敛曲线不能只看某一个目标的最小值因为非支配解集整体收敛才是关键。但工程上简单起见我常画种群中最小的一个目标随时间的变化曲线比如最小转矩脉动随代数下降。figure; plot(1:max_gen, best_hist(:,1), LineWidth, 1.8, Color, [0.85 0.32 0.1]); xlabel(代数, FontSize, 12); ylabel(最小转矩脉动 (N·m), FontSize, 12); title(NSGA-II收敛曲线, FontSize, 14); grid on; set(gca, FontSize, 11);如果想更专业一点可以画超体积指标Hypervolume随代数增长。但超体积计算需要写个工具函数而且对三维目标计算较慢。我一般只用目标值下降曲线来确认有没有收敛够用了。4.3 结构参数分布图与平行坐标帕累托前沿上的解对应的结构参数分布往往是工程师最关心的。比如哪些参数在解集里变化大哪些参数几乎恒定。平行坐标图是展示高维参数分布的好工具。% pareto_param是Pareto解对应的决策变量矩阵已经按列归一化到[0,1] parallelcoords(pareto_param, Labels, {定子极弧系数,转子极弧系数,气隙长度,绕组匝数,定子轭厚,转子轭厚}, ... Color, k, LineWidth, 0.8); set(gca, FontSize, 10);画出来以后你会看到有些线在某一列特别聚拢说明这些参数对优化结果不那么敏感可以在后续设计里定为固定值把精力放到敏感参数上。比如我那个项目里转子轭厚hyr在Pareto解集里基本都在9~10mm之间后来我就直接定10mm把变量维度减到5个优化速度明显加快。4.4 优化前后转矩波形对比图这是一个很有效的展示图把初始设计和优化后的最优折中解的转矩波形画在同一张图里能直观看到脉动降低效果。转矩波形数据可以从有限元或解析模型中得到。% angle_deg是电角度torque_initial和torque_opt是两个方案的转矩波形 figure; plot(angle_deg, torque_initial, --, LineWidth, 1.5, Color, [0.2 0.5 0.2]); hold on; plot(angle_deg, torque_opt, -, LineWidth, 1.5, Color, [0.9 0.2 0.2]); xlabel(转子位置 (deg), FontSize, 12); ylabel(转矩 (N·m), FontSize, 12); legend(初始设计, NSGA-II优化解, Location, best); grid on; box on; title(ER-SRM转矩波形优化对比, FontSize, 13);对比图的重点是要把坐标轴范围对齐否则视觉上脉动变化被放大或缩小。建议统一用ylim固定为优化前波形上下浮动范围的1.5倍这样两张图才有可比性。5. 实操中的坑与优化技巧5.1 种群规模和代数怎么设置很多人一上来就设1000个个体跑500代以为这样能全局收敛。我试过在ER-SRM优化里这样玩不仅慢而且容易出现数值不稳定的个体。经验上决策变量6个目标3个种群规模80~120就足够了。代数也未必需要很大我用200代跑出来的解集跟400代的结果对比帕累托前沿几乎重叠。过度加大种群只是精细化边缘解对工程选型没啥意义。我建议先用小规模快速验证种群40代数80看看目标函数有没有明显问题。确认没问题后再把规模加到120代数200。这样总调试时间能省一大半。5.2 计算代价太大怎么办前文提到了代理模型这里再补充一个细化方案多阶段采样。先用30个样本建一个粗糙的响应面跑完NSGA-II后从Pareto解里挑3个代表点做真实有限元仿真把仿真结果加入训练集重新训练再运行一次NSGA-II。这样循环2~3次代理模型精度会逐步逼近真实模型。我实测第二轮之后的预测误差能从8%降到2%以内。如果完全不想用代理模型可以退而求其次用解析法快速估算转矩。比如简化磁路法和线性电感模型虽然精度差点但作为优化过程中的排序依据是够的。等优化完再对选出的几个解做有限元核实这是一种比较经济的策略。5.3 常见错误和排查我整理几个我踩过而且周围不少人也在踩的坑。第一个坑是目标函数方向写反。NSGA-II里个体A支配个体B的定义是A的所有目标都不比B差且至少一个目标更优。如果你把最大化目标直接当最小化用算法就会往反方向跑。所以强烈建议在评估函数里写一行注释标明每个目标的物理意义和正向趋势。第二个坑是边界处理。交叉和变异之后子代可能超出变量范围。有些NSGA-II实现里没有自动修复导致目标函数里出现负的气隙长度仿真直接报错。解决办法很简单在交叉变异函数末尾加一个截断% 越界修复 offspring(offspring lb) lb(offspring lb); offspring(offspring ub) ub(offspring ub); offspring(:,4) round(offspring(:,4)); % 匝数取整第三个坑是种群中个别个体目标函数返回NaN。通常是因为代理模型预测时遇到输入量落在外推区域或者有限元网格剖分失败。一定要在评估函数里加防护if ~all(isfinite(f)) f 1e10 * ones(1, n_obj); % 惩罚 end否则NaN会在非支配排序里传播导致排序结果崩掉。5.4 算法参数调整心得NSGA-II里有几个参数很关键模拟二进制交叉的分布指数θ有的代码里叫eta_c和变异分布指数θ_m。我常用值θ_c15θ_m20。交叉概率通常0.9变异概率设为1/n_var大概0.15~0.2。如果发现Pareto前沿分布不均匀比如都聚在一端可以适当调大变异概率让搜索更发散如果发现收敛太慢可以增大交叉概率让种群更快混合。还有一个小技巧把随机种子固定住比如rng(42)这样每次跑出来的结果完全一致方便调试和对比不同参数的效果。等最终参数确定后再放开随机种子跑多次挑选最好的一次结果作为论文数据。6. 个人经验总结这块说说我自己的体会。NSGA-II虽然经典但真正的难点从来不是算法本身而是你前面把问题定义得清不清楚。ER-SRM优化里目标函数和决策变量的选择比算法的精妙度重要得多。我在最初一个版本里把转矩脉动定义成了“最大值减最小值”后来发现这个定义受异常点影响很大优化出来的解在平均转矩上表现很差。后来我改成了“转矩方差的开方除以平均转矩”结果帕累托前沿质量明显提升。所以建议大家花时间验证自己的目标指标是否与物理现象一致。另一个体会是画图代码要尽早写不要等优化跑完了再写。我就是跑了三天才想起来忘了保存每代的种群和目标的快照结果没法画收敛曲线只能重新跑一次。后来我在主循环里每代结束都存一份mat文件哪怕存不下全部群体也至少保存每代的最优目标值历史。这笔时间花得很值。最后一个小技巧优化完不要只盯着帕累托前沿上最靠角落的解。那个解通常很极端比如转矩脉动最低但平均转矩很小根本不实用。建议在解集里用“折中解”选取方法比如求每个目标的最大值最小值归一化后找与理想点各目标单独最优欧氏距离最近的那个解。我一般会在图上用不同颜色标记这个折中解作为后续详细设计的基准方案。这个方法简单有效比人眼在散点图上挑点靠谱得多。