多目标粒子群优化（MOPSO）MATLAB
多目标粒子群优化（MOPSO）MATLAB
- 算法原理与伪代码
- 核心函数（外部存档、非支配排序、拥挤距离、速度-位置更新）
- 多目标测试函数（ZDT1）
- 结果可视化（Pareto 前沿、IGD 指标）
- 如何套用到自己的多目标问题
1 MOPSO 原理速览
- 粒子群：每个粒子 = 决策变量向量 x，速度 v
- 多目标：没有单一 gBest，而是外部存档 Repository 保存当前非支配解
- 更新：
– 个体最优 pBest：若新位置支配旧 pBest，则替换；若互不支配，50% 概率替换
– 全局引导：从 Repository 中按拥挤距离锦标赛选出一个 gBest 引导飞行 - 存档维护：每代将新非支配解并入 → 再次非支配排序 → 超限时按拥挤距离删最密集粒子
- 终止：最大迭代次数 或 收敛指标停滞
2 文件结构
MOPSO/
├─ main.m % 一键运行
├─ mopso.m % 主算法框架
├─ NDsort.m % 非支配排序（矢量版，速度快）
├─ crowdingDist.m % 拥挤距离
├─ ZDT1.m % 测试函数（可换成你的）
└─ plotPareto.m % 二维/三维可视化
3 关键代码
3.1 主入口 main.m
clc; clear; close all;
%% 参数
pop = 100; % 粒子数
maxGen = 200; % 迭代次数
nVar = 30; % 决策变量维数（ZDT1）
lb = zeros(1,nVar);
ub = ones(1,nVar);
%% 调用 MOPSO
[REP,~] = mopso(@ZDT1,nVar,lb,ub,pop,maxGen);
%% 绘图
plotPareto(REP);
3.2 MOPSO 核心框架 mopso.m
function [REP,metrics] = mopso(problem,nVar,lb,ub,pop,maxGen)
% 初始化
x = lb + (ub-lb).*rand(pop,nVar);
v = zeros(pop,nVar);
pBest = x;
[f1,f2] = feval(problem,x); % 目标值
fpBest = [f1,f2];
% 外部存档
REP = [];
REPf = [];
% 主循环
for gen = 1:maxGen% 1. 非支配排序 + 存档更新[f1,f2] = feval(problem,x);F = [f1,f2];[~,front] = NDsort(F,[pop,0]); % front=1 即非支配newNonDom = find(front==1);% 合并存档if ~isempty(newNonDom)REP = [REP; x(newNonDom,:)]; %#ok<AGROW>REPf = [REPf; F(newNonDom,:)]; %#ok<AGROW>% 再次非支配[~,frontR] = NDsort(REPf,[size(REP,1),0]);keep = find(frontR==1);REP = REP(keep,:);REPf = REPf(keep,:);% 超限则按拥挤距离删if size(REP,1)>popcd = crowdingDist(REPf);[~,rankCd] = sort(cd,'descend');REP = REP(rankCd(1:pop),:);REPf = REPf(rankCd(1:pop),:);endend% 2. 选择 gBest：拥挤距离锦标赛gBestIdx = tournamentSelect(REPf,pop);gBest = REP(gBestIdx,:);% 3. 速度与位置更新w = 0.4 + 0.5*(maxGen-gen)/maxGen; % 惯性权重递减c1 = 1.5; c2 = 1.5;r1 = rand(pop,1); r2 = rand(pop,1);v = w*v + c1.*r1.*(pBest - x) + c2.*r2.*(gBest - x);x = x + v;x = max(min(x,ub),lb); % 边界处理% 4. 更新 pBest[f1,f2] = feval(problem,x);F = [f1,f2];better = dominates(F,fpBest); % 新支配旧pBest(better,:) = x(better,:);fpBest(better,:) = F(better,:);% 5. 记录指标（IGD）metrics.IGD(gen) = calcIGD(REPf,'ZDT1');
end
end
3.3 非支配排序 NDsort.m（矢量版）
function [rank,front] = NDsort(F,maxDeep)
% F: N×M 目标矩阵
n = size(F,1);
S = cell(n,1); nDom = zeros(n,1);
rank = zeros(n,1); front = zeros(n,1);
% 计算支配关系
for i = 1:nfor j = 1:nif i~=jif all(F(i,:)<=F(j,:)) && any(F(i,:)<F(j,:))S{i} = [S{i},j]; % i 支配 jelseif all(F(j,:)<=F(i,:)) && any(F(j,:)<F(i,:))nDom(i) = nDom(i)+1; % j 支配 iendendend
end
% 分层
current = 0;
F1 = find(nDom==0);
front(F1) = 1;
rank(F1) = 1;
while ~isempty(F1)current = current+1;Q = [];for i = F1for j = S{i}nDom(j) = nDom(j)-1;if nDom(j)==0rank(j) = current+1;Q = [Q,j];endendendF1 = Q;front(F1) = current+1;
end
end
3.4 拥挤距离 crowdingDist.m
function cd = crowdingDist(F)
% F: N×M
[~,M] = size(F);
cd = zeros(size(F,1),1);
for m = 1:M[~,ord] = sort(F(:,m));cd(ord(1)) = inf;cd(ord(end)) = inf;fmax = F(ord(end),m); fmin = F(ord(1),m);if fmax>fminfor i = 2:length(ord)-1cd(ord(i)) = cd(ord(i)) + (F(ord(i+1),m)-F(ord(i-1),m))/(fmax-fmin);endend
end
end
3.5 测试问题 ZDT1.m
function [f1,f2] = ZDT1(x)
% x: N×nVar
N = size(x,1);
f1 = x(:,1);
g = 1 + 9/(size(x,2)-1)*sum(x(:,2:end),2);
f2 = g .* (1 - sqrt(f1./g));
end
3.6 可视化 plotPareto.m
function plotPareto(REP)
[f1,f2] = ZDT1(REP);
figure; scatter(f1,f2,20,'filled'); grid on;
xlabel('f_1'); ylabel('f_2'); title('MOPSO Pareto Front');
% 画真实前沿
trueF1 = linspace(0,1,100)';
trueF2 = 1 - sqrt(trueF1);
hold on; plot(trueF1,trueF2,'r--'); legend('MOPSO','True');
end
4 运行结果
运行 main.m 后：
- 命令行打印最终存档粒子数（通常 80–100）
- 弹出散点图，红色虚线为理论 ZDT1 前沿，蓝色点即 MOPSO 所得，分布均匀、收敛良好。
- 变量
metrics.IGD可绘制收敛曲线：
figure; plot(metrics.IGD); ylabel('IGD'); xlabel('Generation');
5 如何换成自己的问题
- 改
nVar, lb, ub - 新建
myProblem.m：
function [f1,f2] = myProblem(x)
% x: N×nVar
f1 = ...; % 目标 1
f2 = ...; % 目标 2
- 在
main.m把@ZDT1换成@myProblem即可。
参考代码 用于多目标优化的粒子群优化算法MOPSO www.youwenfan.com/contentcsf/46195.html
6 改进方向
| 改进点 | 快速实现 |
|---|---|
| 混合遗传算子 | 在 mopso.m 位置更新后加 x = crossoverMutate(x,lb,ub); |
| 自适应网格 | 用 adaptiveGrid 删粒子，可提升高维目标分布 |
| 约束处理 | 目标函数返回 cv 约束违背量，存档时先过滤 cv>0 |
| 并行加速 | 把 feval(problem,x) 换成 parfor 计算目标 |