简介这份资源面向需要在 MATLAB 中处理多目标优化问题的学生、科研人员与工程技术人员核心是借助遗传算法求解相互冲突目标之间的帕累托最优解集。压缩包共 15 个文件以 8 个 m 脚本和 7 个 asv 自动备份文件为主整体约 8KB体积轻量便于快速查阅与二次修改。脚本覆盖双目标问题建模、种群初始化、选择、交叉、变异、解码、淘汰以及帕累托前沿计算等完整环节可帮助读者理解 gamultiobj 之外的手写实现思路掌握适应度评估与终止条件设置等关键细节。目前已有 4795 人学习下载适合作为课程作业、算法入门或工程原型开发的参考代码通过对照运行与调参读者能直观看到帕累托前沿的生成过程并在此基础上扩展至 NSGA-II、MOEA/D 等更复杂的多目标进化算法。1. 从一组 Pareto 解说起Matlab 多目标优化到底能干什么很多人第一次接触多目标优化是在一个很具体的场景里电机控制器调参既要响应快又要超调小或者结构件轻量化既要减重又要保证刚度。单目标优化只有一个最优解多目标优化给出的是一整条 Pareto 前沿——你没法说哪个点绝对最好但你能看到“再快一点要牺牲多少超调”。Matlab 在这件事上的价值是把遗传算法、粒子群、模式搜索这些求解器统一到gamultiobj、paretosearch这类函数接口下你不用自己从零写非支配排序和拥挤度计算。这份资源围绕的就是这套东西目标函数怎么写、约束怎么加、参数怎么调、结果怎么读。适合已经会 Matlab 基础语法、但被“多目标”三个字卡住的工程师也适合想从 fmincon 单目标迁移过来的熟手。2. 目标函数与约束的写法从数学描述到可运行代码2.1 先把“多目标”翻译成 Matlab 能吃的函数句柄多目标优化的第一步不是打开工具箱而是把工程语言翻译成数学语言。假设你要优化一个悬臂梁结构两个目标质量最小、端部挠度最小。设计变量是截面宽 b 和高 h约束是应力不超过许用值。Matlab 的gamultiobj要求目标函数返回一个行向量每个元素对应一个目标值且默认求最小值。如果你原本的指标是“越大越好”比如效率、增益必须取负号。function f beamObjectives(x) % x(1) b 截面宽, x(2) h 截面高 b x(1); h x(2); L 1.0; % 梁长 m rho 7850; % 密度 kg/m^3 E 2.1e11; % 弹性模量 Pa P 5000; % 端部载荷 N mass rho * b * h * L; % 目标1质量 deflection 4 * P * L^3 / (E * b * h^3); % 目标2端部挠度 f(1) mass; f(2) deflection; end这段代码里有两个容易翻车的地方。第一f必须预先分配成行向量如果写成f [mass; deflection]返回列向量gamultiobj会直接报维度错误。第二目标值的量级差异不能太大。质量可能是几十千克挠度可能是零点几毫米如果直接放一起遗传算法会偏向数量级大的那个目标Pareto 前沿会挤在一端。常见做法是对每个目标做归一化或者用optimoptions里的UseParallel配合自定义的适应度函数做尺度调整。约束的写法分两种。线性约束用矩阵形式传给A*x b、Aeq*x beq非线性约束写成一个单独函数返回[c, ceq]其中c 0、ceq 0。下面这个应力约束就是非线性约束的典型写法。function [c, ceq] beamConstraints(x) b x(1); h x(2); L 1.0; P 5000; sigma_allow 235e6; % 许用应力 Pa sigma 6 * P * L / (b * h^2); % 最大弯曲应力 c sigma - sigma_allow; % 不等式约束要求 0 ceq []; % 无等式约束 end参数说明c返回的是违反约束的程度负值表示满足正值表示违反。ceq如果为空也要写[]不能省略。很多人在这里踩坑把c写成sigma_allow - sigma结果符号反了算法会往违反约束的方向跑最后得到的解全部不可行。2.2 求解器选型gamultiobj 还是 paretosearchMatlab 里做多目标优化主要有两个入口gamultiobj基于 NSGA-II 的变体适合目标函数不连续、非凸、甚至带噪声的情况paretosearch基于模式搜索适合目标函数光滑但梯度不好求、或者设计变量维度较低的场景。选型不是拍脑袋看三个指标变量维度、函数评估代价、是否需要整数变量。对比项gamultiobjparetosearch算法基础遗传算法NSGA-II 变体模式搜索适合变量维度中高维5~50低维2~15函数评估次数多但可并行相对少整数变量支持支持IntCon不支持典型场景多目标参数整定、结构优化小规模光滑问题我一般会先用paretosearch快速探一下前沿形状如果发现前沿不连续或者有多个断档再换gamultiobj加种群规模重跑。反过来先跑gamultiobj再换paretosearch的情况也有但少因为遗传算法的评估次数通常更大先跑贵的没必要。调用时的核心参数在optimoptions里设置。下面是一段可复现的调用代码。nvars 2; lb [0.02, 0.05]; % 宽和高下限 m ub [0.10, 0.30]; % 宽和上限 m options optimoptions(gamultiobj, ... PopulationSize, 100, ... % 种群规模默认 50 MaxGenerations, 200, ... % 最大代数 ParetoFraction, 0.35, ... % Pareto 解占比 CrossoverFraction, 0.8, ... % 交叉概率 UseParallel, true, ... % 开启并行 Display, iter); % 迭代显示 [x, fval, exitflag, output] gamultiobj(... beamObjectives, nvars, [], [], [], [], lb, ub, ... beamConstraints, options);参数说明PopulationSize太小会导致前沿覆盖不全太大则评估次数暴涨。我一般从 100 起步如果前沿点少于 20 个就加到 200。ParetoFraction控制种群中保留多少非支配解默认 0.35如果前沿形状复杂可以降到 0.3 增加多样性。UseParallel需要 Parallel Computing Toolbox开启后每个种群的适应度评估并行执行评估代价高的目标函数能省不少时间。3. 结果解读与前沿可视化别只看散点图3.1 从 fval 矩阵里提取有效信息gamultiobj返回的fval是一个 N×M 的矩阵N 是 Pareto 解个数M 是目标个数。很多人拿到这个矩阵直接plot(fval(:,1), fval(:,2))就完事了但这样会漏掉两个关键信息解的分布均匀度和端点解。均匀度决定你能不能从前沿上挑到中间权衡点端点解往往对应某个目标的极端最优值。% 绘制 Pareto 前沿 figure; scatter(fval(:,1), fval(:,2), 40, filled); xlabel(质量 (kg)); ylabel(端部挠度 (m)); grid on; % 计算相邻解的距离评估分布均匀度 sorted_f sortrows(fval, 1); distances sqrt(sum(diff(sorted_f).^2, 2)); fprintf(前沿点平均间距: %.4f\n, mean(distances)); fprintf(前沿点间距标准差: %.4f\n, std(distances));如果标准差远大于均值说明前沿点分布不均匀有些区域点密集、有些区域空白。这时候可以调大PopulationSize或者改用paretosearch的MeshTolerance参数。另一个常被忽略的是output结构体里的output.message它会告诉你算法是因为达到最大代数停止还是因为前沿变化小于阈值提前停止。如果是前者说明代数不够需要加MaxGenerations。3.2 用决策变量反查设计参数Pareto 前沿上的每个点都对应一组设计变量x。如果你只看了fval没看x相当于拿到了菜单但没点菜。实际工程里你需要从前沿上挑一个折中点然后反查对应的 b 和 h 是多少。% 假设要挑质量不超过 15kg 且挠度最小的解 idx find(fval(:,1) 15); [~, best_idx] min(fval(idx, 2)); selected_x x(idx(best_idx), :); fprintf(选中解: 宽 %.4f m, 高 %.4f m\n, selected_x(1), selected_x(2)); fprintf(对应质量 %.2f kg, 挠度 %.6f m\n, ... fval(idx(best_idx),1), fval(idx(best_idx),2));这段代码的逻辑是先按质量筛选再在满足条件的解里找挠度最小的。注意idx(best_idx)这个写法因为idx是原始索引的子集best_idx是子集内的位置直接混用会取错行。血泪经验我见过有人在这里取错解把最优设计参数用到了错误的目标值上结果实物测试和仿真对不上。4. 避坑与排查多目标优化里最容易翻车的五件事4.1 现象算法跑完得到的解全部不可行原因通常是约束函数的符号写反了或者约束函数的输入维度不对。gamultiobj调用非线性约束时传入的x是行向量如果你在约束函数里按列向量索引比如x(:,1)会报错或者取到错误的值。解决方法是统一用x(1)、x(2)这种线性索引并且在函数开头加一句x x(:);强制转成行向量。4.2 现象Pareto 前沿只有两三个点原因可能是PopulationSize太小或者ParetoFraction设得太低。另一个容易被忽略的原因是目标函数返回了复数。比如挠度计算里b或h被遗传算法扰动成了负数开方或幂运算后产生复数gamultiobj会直接丢弃这些个体。解决方法是在目标函数里加边界检查或者用lb严格限制下限为正数。4.3 现象每次运行结果都不一样无法复现遗传算法本身是随机算法默认使用当前时间作为随机种子。如果你需要复现某次结果必须在调用前设置随机种子。注意gamultiobj内部可能使用多个随机流只设rng(42)有时不够还需要在optimoptions里指定UseVectorized, false并配合rng的twister生成器。rng(42, twister); % 然后再调用 gamultiobj4.4 现象并行加速后结果反而变差UseParallel开启后种群个体的评估顺序被打乱虽然不改变算法逻辑但如果你在目标函数里用了persistent变量或者依赖全局状态并行执行会导致状态不一致。解决方法是把目标函数写成纯函数所有输入通过参数传入不依赖任何外部可变状态。4.5 现象前沿上出现明显不合理的解比如质量很小但挠度也极小理论上不可能同时满足。这通常是因为约束函数没有正确传递或者nonlcon的位置参数写错了。gamultiobj的调用签名里nonlcon是第九个参数如果中间某个[]漏写了约束函数会被当成其他参数解析算法实际上在无约束运行。检查方法是看output里的constrviolation字段如果大于零说明最终解违反了约束。5. 进阶技巧把 Pareto 前沿变成可决策的工程方案5.1 用聚类从前沿上挑代表解前沿点太多的时候逐个看是不现实的。我一般用 k-means 对fval做聚类每个簇取一个代表点这样能把几十个解压缩到五六个方案。聚类前要对目标值做归一化否则量级大的目标会主导距离计算。% 归一化 fval_norm (fval - min(fval)) ./ (max(fval) - min(fval)); % k-means 聚类k5 k 5; [idx_cluster, C] kmeans(fval_norm, k); % 每个簇取距离簇中心最近的解 representatives zeros(k, size(fval,2)); for i 1:k cluster_points fval_norm(idx_cluster i, :); [~, min_idx] min(sum((cluster_points - C(i,:)).^2, 2)); original_idx find(idx_cluster i); representatives(i,:) fval(original_idx(min_idx), :); end disp(代表解目标值); disp(representatives);参数说明k的选择没有绝对标准一般取 3 到 7 之间对应你能向团队展示的方案数量。归一化用 min-max 是常见做法但如果某个目标存在极端离群值改用 z-score 更稳。5.2 用代理模型加速高代价优化如果单次目标函数评估需要跑几分钟的仿真gamultiobj动辄几千次评估是扛不住的。常见做法是先跑少量样本构建 Kriging 或 RBF 代理模型用代理模型代替真实仿真做优化再对 Pareto 解做真实评估验证。Matlab 的fitrgp可以快速构建高斯过程回归模型配合predict使用。% 假设已有样本 X_samples 和对应目标值 Y_samples gprMdl fitrgp(X_samples, Y_samples(:,1), KernelFunction, squaredexponential); % 用代理模型预测新点的目标值 predicted_mass predict(gprMdl, x_new);注意代理模型的精度高度依赖样本分布如果样本集中在设计空间一角预测外推区域会严重失真。我一般会先用 Latin Hypercube 采样铺满整个设计空间样本量至少是变量维度的 10 倍。5.3 验证 Pareto 解是否真的非支配算法返回的 Pareto 解理论上应该是非支配的但数值误差可能导致某个解被另一个解支配。写一个简单的非支配检查函数对fval逐对比较如果发现被支配的解说明算法收敛有问题需要检查目标函数是否有噪声或者约束是否一致。function is_dominated checkDominance(fval) n size(fval, 1); is_dominated false(n, 1); for i 1:n for j 1:n if i ~ j if all(fval(j,:) fval(i,:)) any(fval(j,:) fval(i,:)) is_dominated(i) true; break; end end end end end从那以后我每次跑完gamultiobj都会先跑一遍这个检查确认没有解被支配再看前沿形状。这个习惯帮我省过好几次返工——有一次目标函数里有个除法在边界处产生了 Inf算法把 Inf 当成了大目标值结果前沿上混进了几个被支配的假解。希望帮到你。本文还有配套的精品资源点击获取
