企业官网建设流程全解析

简介：MOPSO（多目标粒子群优化）算法的MATLAB实现源码，面向需要学习多目标进化计算的研究生、工程师及竞赛爱好者，可用于求解帕累托前沿等典型优化问题。该程序结合粒子群全局搜索与帕累托支配机制，能够同时优化多个冲突目标并获得一组分布均匀的非支配解。压缩包大小仅5KB，由10个.m脚本文件组成，除主程序外还包含网格创建、个体支配判断、外部档案去重等辅助模块；该资源已有517人浏览学习。源码覆盖初始化、适应度评价、速度与位置更新、外部档案维护等核心环节，并附有fun.m测试函数，便于读者替换为自定义目标函数，直观观察非支配排序和网格法对多样性的影响。整体代码结构清晰、模块化程度高，适合初学者逐段调试，也适合在此基础上扩展改进，进行多目标优化实验与研究。

1. 为什么MOPSO不是多跑几次PSO那么简单

多目标粒子群优化这些年常被当成“PSO加个权重”来用，实际效果往往只在凸前沿上过得去，凹前沿直接漏掉一片解。MOPSO的核心变化在于把“唯一最优”换成“一组互不支配的解”，用支配关系决定粒子保留方向，再用外部档案维护整个帕累托前沿。这套 MATLAB 程序不是玩具代码，它把非支配排序、网格划分、去重和个体选择拆成了独立函数，适合做结构优化、资源配置、参数标定这类需要多个冲突目标平衡的人。对刚接触多目标优化的工程师，看完主循环就能把 PSO 的直觉平移过来；想深入研究的人，也可以拿 ZDT 系列测试函数逐行验证。

2. 从粒子群到帕累托：MOPSO的数学基础与支配判断

2.1 PSO的更新公式在多目标里卡在哪

经典粒子群在单目标问题中的迭代核心是： v_i^(t+1) = w v_i^t + c1 r1 (pBest_i - x_i^t) + c2 r2 (gBest - x_i^t) x_i^(t+1) = x_i^t + v_i^(t+1)

w是惯性权重，c1/c2是学习因子，r1/r2是[0,1]随机数。这个式子里只有一个gBest，所以每次迭代所有粒子都往同一个最优位置拉。把多目标问题加权成单目标，再跑多次不同权重，看起来像是“多跑几次PSO”，但它只能得到凸Pareto前沿上的解，凹前沿里的点无论取什么权重都不会成为单目标加权和的最优解。所以MOPSO必须在算法内部保留一组解，而不是外部枚举权重。

MOPSO把gBest从单一变量换成了外部档案中的一个引导者。每个粒子在每一代从档案中按概率选出一个gBest，于是不同粒子可以飞向不同区域，前沿才能铺开。pBest的定义也在变：当新目标向量支配旧pBest时，直接替换；当两者互不支配时，常见做法是保留离当前粒子更近的那个，或者随机保留一个，目的是维持粒子个体方向的稳定性。

2.2 支配关系与非支配排序

最小化问题里，解A支配解B的条件是：对每个目标都有A的目标值小于等于B，且至少有一个目标严格小于。代码里最基础的就是Domination.m，它一次只判断两个解：

function d = Domination(x, y) % 输入 x, y：两个解的目标函数值，例如 [f1, f2] % 输出 d = 1 表示 x 支配 y，d = 0 表示不支配 better = x < y; % 严格更优的目标 worse = x > y; % 严格更差的目标 if any(better) && ~any(worse) d = 1; else d = 0; end end

这里默认所有目标都是最小化；如果你的问题里有最大化目标，可以直接取负号再传入，避免在函数里反复判断。any(better) && ~any(worse)的意思是：至少在一个目标上严格更优，同时没有任何一个目标更差，这才是支配。两个条件缺一不可。

JudgePopDomination.m 的作用是把这层判断放到整个种群上。它的常见实现是：遍历每个粒子A，与其余所有粒子B逐个调用Domination，如果B支配A，就把A标记为“被支配”；如果A不支配所属的集合，那就不写入非支配层。更规范的做法是会做一次非支配排序，分出第一层、第二层前沿，但这套代码里只需要“是/不是非支配解”两步判断，所以JudgePopDomination这个名字足够说明问题：它只返回每个粒子的支配状态，而不是完整排序。

为了保证档案里的解彼此不重复，uniqueRep.m 会把目标函数向量完全一样的粒子合并。这里有个工程细节：用浮点数比较相等会出问题，因为0.99999999和1.00000001明明应该算同一个点。我一般会引入一个容忍度，比如round(cost*1e6)/1e6之后再比对，或者把目标值量化到网格精度。DeleteOneRepMember.m 则在档案容量超过设定值时，从拥挤网格中删除冗余成员，目标也很直白：别让相同的解占满档案。

2.3 网格法与外部档案的均匀性控制

只有一个放满非支配解的数组还不够，gBest到底从数组里哪一项选，直接决定前沿是否均匀。这套程序使用了自适应网格：CreateGrids.m 先根据当前档案的目标值范围，把目标空间切成若干个格子；FindGridIndex.m 再计算每个粒子落在第几个格子。Selectzbest.m 在选全局引导时，优先从粒子数量少的格子里随机挑一个解作为gBest。这样做的直觉是：稀疏区域的解更容易得到引导，粒子会陆续飞过去，把前沿的空洞补起来。

网格选择不是唯一方案，几种常见策略可以放在一起对比：

网格法最大的优点是不用算距离矩阵，计算量集中在网格索引上，在大档案场景下更划算。缺点集中在两个地方：一是网格数量nGrid设置过小，格子太粗，选出来的gBest差异性差；二是nGrid设置过大，很多格子只有一两个解，稀疏判断会失效。实际操作里，两个目标取5到15个格子，三个目标取4到8个格子，后续跑ZDT1时会更快看到差别。

另外，CreateGrids 在取目标上下界时，常见做法是在min和max基础上扩大一个alpha*span的边界。原因是一旦某个解恰好落在边界上，找邻居格子时可能会越界。有的实现会把alpha设成0.1或者0.5，具体要看目标范围是否归一化；不扩边界也经常会跑错，只是偶尔报错不够明显，结果前沿少了一块。

3. 逐文件拆解这套MATLAB版MOPSO程序：函数职责与调用链

当主程序 PSO_Multipleob2j.m 一打开，扑面而来的是一堆.m文件。这套代码和教学版PSO最大的差别在于，它把每个可复用的判断都拆成了独立函数。对学习来说这很友好，调通之后你完全可以只改fun.m和一个参数段，就能解决另一个领域的问题。

3.1 文件清单与职责

先把整个项目的文件对应关系列清：

主程序并不直接调用每个函数，调用关系是：主循环 → Selectzbest/CreateGrids/FindGridIndex/JudgePopDomination/uniqueRep/DeleteOneRepMember。fun.m 和 Domination.m 在最底层，几乎被其他函数共用。

3.2 主循环结构与速度更新

主程序里最值得抄的部分是迭代主循环。我把它简化成可读版本：

% PSO_Multipleob2j.m 主流程（简化后可照抄骨架） nVar = 30; % 决策变量维度 VarMin = zeros(1, nVar); VarMax = ones(1, nVar); N = 100; % 粒子数量 MaxIt = 200; % 最大迭代次数 w = 0.5; c1 = 1.5; c2 = 1.5; nGrid = 10; nRep = 100; % 初始化 particle(N).Position = []; particle(N).Velocity = []; for i = 1:N particle(i).Position = unifrnd(VarMin, VarMax); particle(i).Velocity = zeros(1, nVar); particle(i).Cost = fun(particle(i).Position); particle(i).Best.Position = particle(i).Position; particle(i).Best.Cost = particle(i).Cost; end % 外部档案初始化为当前所有非支配解 rep = particle; repCosts = vertcat(rep.Cost); keep = JudgePopDomination(repCosts); rep = rep(logical(keep)); rep = uniqueRep(rep); rep = DeleteOneRepMember(rep, nRep); Grid = CreateGrids(rep, nGrid); for it = 1:MaxIt for i = 1:N gBest = Selectzbest(rep, Grid); r1 = rand(1, nVar); r2 = rand(1, nVar); particle(i).Velocity = w * particle(i).Velocity ... + c1 * r1 .* (particle(i).Best.Position - particle(i).Position) ... + c2 * r2 .* (gBest.Position - particle(i).Position); particle(i).Position = particle(i).Position + particle(i).Velocity; particle(i).Position = min(max(particle(i).Position, VarMin), VarMax); particle(i).Cost = fun(particle(i).Position); if Domination(particle(i).Cost, particle(i).Best.Cost) particle(i).Best.Position = particle(i).Position; particle(i).Best.Cost = particle(i).Cost; end end % 把个体最优并入档案，再统一做支配筛选、去重、容量控制 rep = [rep, particle]; repCosts = vertcat(rep.Cost); keep = JudgePopDomination(repCosts); rep = rep(logical(keep)); rep = uniqueRep(rep); rep = DeleteOneRepMember(rep, nRep); Grid = CreateGrids(rep, nGrid); end

这份代码里有两个容易出错的地方：一是unifrnd(VarMin, VarMax)在两个输入都是向量时，直接生成的是对应维度的随机矩阵，要保证VarMin/VarMax的维度对齐；二是更新位置后的min(max(...))只是硬边界，粒子在边界处速度不归零，容易反复撞边界，实际中我会在边界处让速度乘以 -0.5，避免粒子一直贴在边界上。

3.3 Selectzbest、CreateGrids与FindGridIndex的配合

Selectzbest 每次被调用时，传入的是外部档案rep和当前网格Grid。它先调用 FindGridIndex 算出每个档案粒子的网格编号，再统计每个网格里的粒子数。比如三维目标空间里nGrid=10时，理论上有1000个格子，但实际只有少数格子有解。Selectzbest 会找出包含非支配解且粒子数最少的格子，在其中一个解里随机选一个返回。这段逻辑保证了稀疏区域有更高的概率被选为全局最优，粒子也就更容易往那里飞。

CreateGrids 则承担一个容易被忽略的任务：每个维度上设置nGrid等分。两个目标时它返回两个向量下边界和上边界，或者返回一个包含网格边界与编号规则的结构体。FindGridIndex 内部常见写法是：

function idx = FindGridIndex(cost, Grid) % 假设 Grid.lb 和 Grid.ub 是每维边界向量的 cell 数组 nObj = numel(cost); idx = zeros(1, nObj); for j = 1:nObj id = find(cost(j) >= Grid.lb{j} & cost(j) < Grid.ub{j}, 1); if isempty(id) idx(j) = -1; % 超出边界，说明网格没有外扩 else idx(j) = id; end end end

这里麻烦的地方在于，当某个目标值恰好等于上边界时，<会把最后一个点排除在网格外。所以 CreateGrids 里一定要把最后一格的上边界比实际最大值放大一点，常见做法是ub(end) = ub(end) + eps，或者用alpha*range外扩。遇到前沿呈长尾分布时，这个细节直接影响FindGridIndex能否正确索引。

可以顺带提一句：如果自己改问题，目标数量从2变成3时，CreateGrids和FindGridIndex的维度逻辑要同步调整，否则索引维度对不上。这套代码的目标函数fun.m如果返回两列，那nGrid可以设成10；如果返回三列，nGrid要降下来，否则网格数量指数增长，大部分格子都是空的。

4. 跑通并检验：ZDT系列测试函数上的参数设置与结果可视化

MOPSO写得好不好，不能光看能不能跑，要看它能不能在标准测试函数上逼近真实帕累托前沿。ZDT1是最常用的入门测试函数，有两个优化目标、变量维度可以自由设置。拿它来验证这套MATLAB程序，半小时内就能得到一张可以写进报告的前沿图。

4.1 如何把ZDT1写进fun.m

ZDT1的数学定义是，决策变量x1到xn都取值[0,1]： f1(x)=x1 g(x)=1+9/(n-1)sum(x2..xn) f2(x)=g(1-sqrt(f1/g))

它的真实Pareto前沿是f2=1-sqrt(f1)，对应的条件是x2到xn全等于0。代码实现可以这样写：

function z = fun(x) % ZDT1 测试函数，两个目标均为最小化 n = numel(x); f1 = x(1); g = 1 + 9 / (n - 1) * sum(x(2:end)); f2 = g * (1 - sqrt(f1 / g)); z = [f1, f2]; end

这段代码返回的 z 是一个1×2的行向量。主程序里判断支配关系时，用的就是这种目标向量。变量n取多少会影响求解难度：n=30是文献常用设置，但收敛慢，适合正式对比；n=10能快速看到趋势，适合调参数。第一次跑建议先用n=10。

4.2 主参数设置与影响

MOPSO不像单目标PSO那样只有“收敛”一个指标，它要同时看收敛性和分布性。所以参数表格比普通PSO多出nGrid和nRep两行。

这里我最常遇到的问题是 nGrid 和 nRep 不匹配：nGrid很小（比如5），nRep却设到200，最后会有大量粒子被塞进同一个格子，Selectzbest 的选择压力失效；反过来 nGrid很大、nRep只有30，档案里每个格子只有一两个解，稀疏判断会变成随机选择。先按 nGrid = 10、nRep = 100 起步，再看前沿图调整。

4.3 运行步骤与可视化对比

一套可复现的运行步骤是：

1. 解压整个MOPSO程序包，确认所有.m文件在同一个目录下，路径上不要出现中文文件夹名，否则MATLAB读取容易出问题。
2. 用上面的ZDT1代码覆盖fun.m。
3. 打开PSO_Multipleob2j.m，把nVar设为10，MaxIt设为200，其余按默认参数。
4. 运行脚本，主循环完成后在工作区里应看到rep变量。
5. 执行Plotfitness(rep)，观察Pareto前沿是否从左下角延伸到右上角。

如果主脚本运行完直接结束，没留下rep，就在末尾追加一行保存命令：

save('MOPSO_result.mat', 'rep');

然后用单独的画图脚本对比真实前沿：

load('MOPSO_result.mat'); costs = vertcat(rep.Cost); % rep中的Cost字段组装成 N×2 矩阵 real_x = linspace(0, 1, 200); real_y = 1 - sqrt(real_x); plot(real_y, real_x, 'k--', 'LineWidth', 1.5); hold on; scatter(costs(:,2), costs(:,1), 15, 'filled', 'MarkerFaceAlpha', 0.6); xlabel('f2'); ylabel('f1'); legend('真实 Pareto 前沿', 'MOPSO 结果', 'Location', 'best');

注意我画图时把f1放在纵轴、f2放在横轴，是为了和ZDT1的常见图形保持一致。如果你习惯f1横轴，就把real_x、real_y和scatter里的列对调。重点是：MOPSO求出的点应该落在真实前沿附近，而不是聚在某一小段。如果点全部挤在左侧，说明w太大或者MaxIt不够；如果点分散但离前沿很远，优先检查fun.m的g表达式里的除法是否写成了矩阵除法。

热搜词里常出现“mopso优化zdt1的pareto前沿”，其实就是在做这一件事：把标准测试函数跑出接近真实前沿的一组解，再用图形把收敛性和分布性展示出来。运行这一步时，会发现档案里有些点的f1分布并不均匀。这不一定代表算法错了，如果nGrid太小，网格选择无法区分稀疏区域，前沿中央会缺一段。把nGrid调到15再看一次，通常能看到明显变化。

5. 把MOPSO用到自己的问题：目标函数改写、网格扰动与性能验证

5.1 替换fun.m时的三个检查点

换掉ZDT1，改成自己的工程问题，只需要动fun.m和变量范围。三个地方最容易出错：第一，目标函数必须返回一行向量，内部有多少个目标就返回多少列，MOPSO不自动识别向量长度；第二，VarMin和VarMax的长度必须和决策变量维度一致，否则unifrnd会报错；第三，如果有约束条件，不要试图把约束写进支配判断里，常见做法是在fun.m返回值上叠加惩罚项，让违反约束的解在目标空间中处于支配劣势。

5.2 网格扰动实验

我看一个MOPSO实现是否稳定，会先做网格扰动实验：固定其他参数，把nGrid从7改成15，再改成25，观察同一问题下的Pareto前沿变化。nGrid=25时常常出现“前沿看起来很散，但IGD反而变好”的情况，这是因为细网格让更多粒子去填补稀疏区域，虽然输出点数量不变，位置更贴前沿。这个扰动实验也能用来判断nGrid是否过大：如果nGrid从15到25，结果几乎没变化，说明再细分已经没有信息量。

5.3 用IGD指标代替肉眼判断

视觉对比只能定性，定量验证可以算IGD（Inverted Generational Distance）。IGD衡量MOPSO解集到真实前沿的逼近程度，值越小越好。真实前沿通常用密集采样点表示，ZDT1可以直接用linspace生成：

function igd = IGD(PF, PFtrue) % PF：MOPSO求得的非支配解目标值矩阵（N行nObj列） % PFtrue：真实Pareto前沿密集采样点（M行nObj列） d = 0; for i = 1:size(PFtrue, 1) dist = sqrt(sum((PF - PFtrue(i, :)).^2, 2)); d = d + min(dist); end igd = d / size(PFtrue, 1); end

调用时传入真实前沿点和MOPSO结果，得到的是一个标量。我通常跑10次取中位数，避免单次随机性误导结果。调参时有个简单技巧：固定其他参数，只改变nRep，观察IGD曲线；多数情况下nRep越大IGD先变好后变差，因为档案过大时网格选择压力被分散，所以可以把nRep设为三倍目标解数量作为起点，再用这个小实验找到拐点。如果你发现前沿某一段始终缺失，优先检查CreateGrids的边界放大系数，我一般把alpha从0.1加到0.3，缺失段往往会自己修复。

本文还有配套的精品资源，点击获取