多目标粒子群优化（MOPSO）MATLAB

多目标粒子群优化（MOPSO）MATLAB

1. 算法原理与伪代码
2. 核心函数（外部存档、非支配排序、拥挤距离、速度-位置更新）
3. 多目标测试函数（ZDT1）
4. 结果可视化（Pareto 前沿、IGD 指标）
5. 如何套用到自己的多目标问题

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)>pop
            cd = crowdingDist(REPf);
            [~,rankCd] = sort(cd,'descend');
            REP  = REP(rankCd(1:pop),:);
            REPf = REPf(rankCd(1:pop),:);
        end
    end
    % 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:n
    for j = 1:n
        if i~=j
            if all(F(i,:)<=F(j,:)) && any(F(i,:)<F(j,:))
                S{i} = [S{i},j];      % i 支配 j
            elseif all(F(j,:)<=F(i,:)) && any(F(j,:)<F(i,:))
                nDom(i) = nDom(i)+1;  % j 支配 i
            end
        end
    end
end
% 分层
current = 0;
F1 = find(nDom==0);
front(F1) = 1;
rank(F1)  = 1;
while ~isempty(F1)
   current = current+1;
   Q = [];
   for i = F1
       for j = S{i}
           nDom(j) = nDom(j)-1;
           if nDom(j)==0
               rank(j) = current+1;
               Q = [Q,j];
           end
       end
   end
   F1 = 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>fmin
        for i = 2:length(ord)-1
            cd(ord(i)) = cd(ord(i)) + (F(ord(i+1),m)-F(ord(i-1),m))/(fmax-fmin);
        end
    end
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 如何换成自己的问题

1. 改 nVar, lb, ub
2. 新建 myProblem.m：

function [f1,f2] = myProblem(x)
% x: N×nVar
f1 = ...;   % 目标 1
f2 = ...;   % 目标 2

3. 在 main.m 把 @ZDT1 换成 @myProblem 即可。

参考代码 用于多目标优化的粒子群优化算法MOPSO www.youwenfan.com/contentcsf/46195.html

6 改进方向