简介：这份资源面向需要在 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 还是 paretosearch
Matlab 里做多目标优化主要有两个入口：gamultiobj基于 NSGA-II 的变体，适合目标函数不连续、非凸、甚至带噪声的情况；paretosearch基于模式搜索，适合目标函数光滑但梯度不好求、或者设计变量维度较低的场景。选型不是拍脑袋，看三个指标：变量维度、函数评估代价、是否需要整数变量。
| 对比项 | gamultiobj | paretosearch |
|---|---|---|
| 算法基础 | 遗传算法（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 聚类，k=5 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 当成了大目标值，结果前沿上混进了几个被支配的假解。希望帮到你。
本文还有配套的精品资源，点击获取
