简介：面向 MATLAB 多目标优化应用场景，这份资源围绕 gamultiobj 函数与 Pareto 最优前沿展开，适合需要求解多目标问题、理解帕累托解的工程师和研究人员。压缩包体积小巧，共含三个文件：两个脚本文件分别用于多目标问题的定义与优化过程执行，一个文本文件记录求解结果与目标函数值，便于对照分析，整体大小仅 2KB。目前已有 851 人学习下载，说明其具有一定的实践参考价值。通过这份代码示例，读者可以掌握 gamultiobj 的基本调用方式，学习设置变量边界与遗传算法参数，理解 Pareto 最优解集的生成与筛选方法，并直接借鉴代码模板快速部署到自己的多目标优化任务中。对于需要深入了解遗传算法原理与多目标权衡决策的 MATLAB 用户，该示例具有很强的启发性，可作为教学或科研的基础工具。
1. 多目标优化不是“选一个解”，而是“找一整条前沿”
如果你还停留在“优化 = 求最大值 / 最小值”的阶段，gamultiobj 会推翻这个直觉：它一次给出成一个解集，而不是一个点。这个解集叫 Pareto 最优前沿，也就是在相互冲突的多个目标之间，所有“再改善一个目标就必然恶化另一个目标”的点的集合。当目标函数是连续可求值的黑箱，决策变量维度又不算太高（几十个以内）时，用 MATLAB 优化工具箱里的 gamultiobj 函数配合遗传算法，是在不自己实现 NSGA-II 的前提下最靠谱的落地路线。本文面向两类人：一类是论文或课题里需要真正产出 Pareto 前沿图的工程师，另一类是在生产排产、参数标定、结构轻量化中反复被“加权求和”坑过、想把多目标问题一次算干净的研究者。接下来直接讲透这个函数的原理、参数与坑。
2. 从 Pareto 最优前沿到 gamultiobj 的求解框架
2.1 非支配排序与拥挤度：Pareto 前沿是怎么被“筛”出来的
gamultiobj 的底层是 NSGA-II 的变体，核心机制是两个：
第一个是非支配排序。一个解 a 支配另一个解 b，当且仅当 a 在所有目标上都不比 b 差，且至少在一个目标上严格优于 b。把所有互不支配的解放在一起，就是第一层非支配前沿，也就是当前种群中的 Pareto 最优前沿。去掉第一层后重复这个过程，就得到第二层、第三层。种群选择时，优先保留层级靠前的个体。
第二个是拥挤度距离。同一层前沿上的解无法用支配关系区分优劣，这时需要衡量每个解周围解的密集程度。拥挤度越大，说明周围越空旷，保留它有利于维持前沿的延展性。gamultiobj 在每一代的选择、交叉、变异后都执行这两个步骤，迭代若干代后收敛到近似 Pareto 前沿。
这里有个关键认知：gamultiobj 输出的并不是数学上的精确 Pareto 最优集，而是一个分布良好、收敛到接近真实前沿的近似集。因此评估结果时，不仅要看解的质量，还要看解的分布是否覆盖了完整的目标区间。如果你的前沿明显缺失某一端目标函数的取值段，说明种群多样性不够，通常要加大 PopulationSize 或 MaxGenerations。
2.2 最小可运行代码：从“Himmelblau 双目标测试问题”看完整流程
不必先啃进化算法的数学细节，先跑通一个最小例子。考虑经典双目标问题：
- min f1(x) = (x1^4 + x2^4 - 4x1^2 - 4x2^2 + x1*x2)
- min f2(x) = x1^2 + x2^2 - 4
变量范围均为 [-3, 3]。用 gamultiobj 求解：
fun = @(x) [x(1)^4 + x(2)^4 - 4*x(1)^2 - 4*x(2)^2 + x(1)*x(2); x(1)^2 + x(2)^2 - 4]; nvars = 2; lb = [-3, -3]; ub = [3, 3]; options = optimoptions('gamultiobj', ... 'PopulationSize', 100, ... 'MaxGenerations', 200, ... 'Display', 'iter', ... 'UseParallel', true); [x, fval, exitflag, output] = gamultiobj(fun, nvars, [], [], [], [], lb, ub, [], options); plot(fval(:,1), fval(:,2), 'o'); xlabel('f1'); ylabel('f2'); grid on;说明几点：
- 目标函数句柄 fun 返回一个列向量，每个元素对应一个目标值。这是 gamultiobj 的硬性输入格式
- 线性约束 A、b 和 Aeq、beq 都可以是空矩阵，但位置不能省
- 返回的 fval 是 Pareto 前沿上的目标函数值矩阵，行数与 x 相同，列数等于目标数。
plot(fval(:,1), fval(:,2))画的即是双目标下的 Pareto 最优前沿分布
2.3 从输出结构体里把 Pareto 解集捞出来
很多人只拿 fval 画了图就交差，但后续做决策时还需要决策变量本身。gamultiobj 的输出里，x 就是与 fval 一一对应的 Pareto 解集，每一行是一组决策变量。output 结构体里的信息也值得拆开看：
output.message % 终止信息，可判断是否达到 TolFun 或 MaxGenerations output.funccount % 目标函数总调用次数，评估计算量用的核心指标 output.rngstate % 当前随机数状态，复现实验全靠它在我的实践中，一个常见的失误是只看 exitflag。gamultiobj 的 exitflag 为 1 只代表达到 MaxGenerations 被正常终止，并不代表收敛质量好。要判断是否收敛，看output.funccount是否有大量浪费，以及最后前沿随迭代的变化是否趋于稳定。更保险的做法是重跑一次，对比两次前沿的分布差异。
3. 把工程问题翻译成 gamultiobj 可解的约束与目标形式
3.1 目标函数怎么写：这是一个纯“黑箱优化”问题吗
gamultiobj 对目标函数的要求比 fmincon 宽松得多——不需要梯度，甚至不需要目标函数连续可微。这意味着你可以把 Simulink 仿真、有限元计算脚本、外部可执行程序封装在自定义函数里作为目标。
常见的做法是写一个独立 .m 文件：
function y = cost_fun(x) % x 为决策变量行向量 % 调用外部仿真，返回两个目标值 sim_result = run_simulation(x); % 自定义的仿真封装 y(1) = sim_result.energy_consumption; y(2) = sim_result.processing_time; end注意，gamultiobj 默认认为所有目标都是最小化问题。如果你的某个目标实际要最大化，最简单的方式是取其相反数，并在画图时把坐标轴翻转。工具箱里的optimoptions并不支持对单个目标设置方向，所以取负是目前最干净的做法。
3.2 约束怎么给：线性约束优先，非线性约束要注意计算成本
gamultiobj 支持三类约束：
- 边界约束 lb、ub：必须显式给出，能大幅加速收敛，也能避免遗传操作产生明显不可行解
- 线性不等式与等式约束 A、b 和 Aeq、beq：写在调用参数第 4 到第 7 个位置
- 非线性约束 nonlcon：返回 [c, ceq]，c ≤ 0 表示可行
非线性约束的坑在于代价高。遗传算法中每一代都要对大量个体调用约束函数，如果约束本身依赖另一个仿真，计算成本会爆炸。我一般先用简单的数学表达式粗略过滤，再用真实仿真只验证 Pareto 前沿上的候选解。
3.3 五个必调参数：不要只用默认值
gamultiobj 的默认参数对玩具问题够用，但工程问题上常常需要手工调整。我整理出优先级最高的几个：
| 参数 | 默认值 | 建议设置 | 理由 |
|---|---|---|---|
| PopulationSize | 50 | 100–200 | 种群太小前沿稀疏；太大则每代计算量线性上升 |
| MaxGenerations | 100 | 200–1000 | 复杂前沿需要更多代收敛，尤其在存在约束时 |
| ParetoFraction | 0.35 | 0.3–0.7 | 控制前沿保留比例；比例越小越偏向收敛，越大越偏向多样性 |
| FunctionTolerance | 1e-4 | 1e-5–1e-6 | 用于判断前沿平均变化，设太大会提前终止 |
| UseParallel | false | true | 多核并行评估种群，大幅缩短计算时间 |
以下是一个启动并行池并设置参数的典型写法：
parpool('local', 4); % 按 CPU 核心数调整 options = optimoptions('gamultiobj', ... 'PopulationSize', 150, ... 'MaxGenerations', 500, ... 'ParetoFraction', 0.5, ... 'FunctionTolerance', 1e-5, ... 'UseParallel', true, ... 'Display', 'iter');需要特别说明 ParetoFraction：它控制每一代保留的非支配个体比例，比例越大保留的多样性越好。但如果你发现最终前沿出现大量相距极近的点，说明这个值偏高，前沿拥挤且计算浪费在重复区域，这时适当调小。
3.4 混合整数多目标：gamultiobj 的限制与变通
gamultiobj 的官方实现要求决策变量是连续的。如果问题中混有整数变量（例如设备启停的 0/1 变量），工具箱没有直接的整数支持。常见变通方案有三种：
- 把整数变量传入目标函数后自行四舍五入，但要注意这会造成目标函数在整数边界附近不连续，遗传算法仍然可以处理
- 使用
ga单目标的混合整数支持，配合权重法转成多个单目标分别求解 - 使用全局优化工具箱 3.x 版本中利用模式搜索的多目标算法 paretosearch，可以处理整数
第三种是 R2019a 之后比较省事的方案。我把具体对比放在后面第 5 章。
4. 从 Pareto 前沿到决策：画图、选解与评价指标的完整闭环
4.1 二维与三维前沿的可视化方法
双目标问题直接散点图，三目标问题可以散点 + 颜色映射。
fval_sorted = sortrows(fval, 1); % 按第一个目标排序，有利于连线 plot(fval_sorted(:,1), fval_sorted(:,2), 'o-');注意排序的目的：如果直接 plot 不排序，连线会乱序交叉。工业界还常用平行坐标图看解集特征：
parallelcoords([x, fval]);这个图可以直观看到不同 Pareto 解对应的决策变量分布区间，是判断“哪个变量对目标影响最敏感”的有力工具：若某列平行线束很窄，说明该决策变量可以被锁定到一个小区间。
4.2 从 Pareto 集中挑选“最终解”的三种实用方法
得到 Pareto 前沿只是第一步，实际工程必须从中选一个点去落地。三种方法我都在用：
方法一：拐点法。画双目标前沿时，曲率变化最剧烈的点往往比两端极端点更均衡。
% 把 fval 归一化到 [0,1] fval_norm = (fval - min(fval)) ./ (max(fval) - min(fval)); % 计算每个点到原点的欧氏距离 dist = sqrt(sum(fval_norm.^2, 2)); [min_dist, idx] = min(dist); best_solution = x(idx, :);它选的是归一化后距离原点最近的点，也就是两个目标折衷最均衡的方案。
方法二：加权 TOPSIS。当决策者对各目标有偏好权重时，可以通过加权贴近度排序。
weights = [0.6, 0.4]; score = sum(fval_norm .* weights, 2); [~, idx] = min(score);本质是线性加权法作用于离散 Pareto 集上，而不是对目标函数连续加权。它不破坏 Pareto 前沿，只是提供排序。
方法三：交互式决策。把前沿展示给决策者，从中选点。这看起来不“自动化”，却往往是工业项目里最容易被接受的方式，因为决策者能看到取舍的代价。工程师要做的只是把选中的点的决策变量输出为参数文件，并标明该点的目标值。
4.3 检查约束满足情况：一个容易忽视的步骤
gamultiobj 返回的所有解理论上都应满足约束，但非线性约束在遗传操作中可能生成轻微违约个体，由于约束违反量很小而未被淘汰。稳妥做法是对选中的最终解单独验证：
[c, ceq] = nonlinear_constraints(best_solution); if any(c > 1e-6) || any(abs(ceq) > 1e-6) warning('候选解约束违约，请重新选点'); else accept_solution = best_solution; end这里的 1e-6 是容差，根据实际问题量级可以放大或缩小。这个步骤在写论文时尤其重要，审稿人如果抽查到违约解，会直接质疑整个算法实现。
4.4 评价前沿质量：不只是看一眼图
客观评价 Pareto 前沿质量通常用两个指标，在 MATLAB 里可以自己写几行代码计算：
- IGD（反世代距离）：度量近似前沿与真实前沿的距离，前提是已知真实前沿（测试函数可解析得到或通过网格采样获得）
- Spread（分布度）：度量解在前沿上是否均匀分布，gamultiobj 输出中没有直接给出，但可以自己计算相邻点间的距离方差
% 计算相邻前沿点间距的标准差，越小说明分布越均匀 fval_norm = (fval - min(fval)) ./ (max(fval) - min(fval)); d = sqrt(sum(diff(fval_norm).^2, 2)); spread = std(d);这个指标可以用于调参前后对比：如果调参后 spread 明显变小，说明解集分布更均匀。
5. 用 paretoset 提取最终前沿：一个值得固化的收尾技巧
工作收尾时，有一个小技巧我每次都会用：gamultiobj 的输出里其实已经包含了当前种群所有非支配个体，但默认的返回结果里，x 和 fval 已经过滤到了 Pareto 前沿。不过，在较老版本中，gamultiobj还遗留了一个不对外开放的成员output.paretoset。它可以返回前端解对应的索引，再配合fval和x，能够反查哪些解属于“同一层前沿”，在多目标数目大于 3 时，可以把原始输出进一步精简：
if isfield(output, 'paretoset') idx = output.paretoset; pareto_x = x(idx, :); pareto_fval = fval(idx, :); else % 老版本手动筛选非支配解 n = size(fval, 1); dominated = false(n, 1); for i = 1:n for j = 1:n if i ~= j && all(fval(j,:) <= fval(i,:)) && any(fval(j,:) < fval(i,:)) dominated(i) = true; break; end end end pareto_x = x(~dominated, :); pareto_fval = fval(~dominated, :); end这段代码的逻辑说明：新版本（R2015a 之后）output.paretoset已经存在，它记录的是 Pareto 前沿上的个体在原种群中的索引；如果手动实现非支配筛选，关键在于嵌套循环里用all(... <= ...)和any(... < ...)两者同时成立来判断被支配。数据量大时，这段双重循环可以向量化优化，但在 Pareto 解数量几千以内时，直接运行完全可接受。
另外，如果你在 R2019a 及之后版本中使用了paretosearch函数，它可以直接与 gamultiobj 做交叉验证：同一问题上跑两个算法，对比前沿中点分布密度的差异。对新问题建模时，我会先用gamultiobj暴力探索解空间，再用paretosearch做局部精修，两者配合往往能省下大量调参时间。这个组合并不复杂，但能显著提升结果的说服力与稳定性。
本文还有配套的精品资源，点击获取
