多目标粒子群优化（MOPSO）解决ZDT1问题

前言

提醒：
文章内容为方便作者自己后日复习与查阅而进行的书写与发布，其中引用内容都会使用链接表明出处（如有侵权问题，请及时联系）。
其中内容多为一次书写，缺少检查与订正，如有问题或其他拓展及意见建议，欢迎评论区讨论交流。

内容由AI辅助生成，仅经笔者审核整理，请甄别食用。

---

matlab代码

clc; clear; close all;% ---------------- 参数设置 ----------------
nVar = 30;         % 决策变量个数
VarMin = 0;        % 最小值--变量范围（ZDT1 要求 x_i ∈ [0,1]）
VarMax = 1;        % 最大值
MaxIt = 100;       % 最大迭代次数
nPop = 50;         % 粒子群规模
nRep = 50;         % 外部存档大小w = 0.5;           % 惯性权重
c1 = 1.5;          % 个人引导因子
c2 = 2.0;          % 群体引导因子
mu = 0.1;          % 变异概率% ---------------- 粒子结构体 ----------------
particle(nPop).Position = [];
particle(nPop).Velocity = [];
particle(nPop).Cost = [];
particle(nPop).Best.Position = [];
particle(nPop).Best.Cost = [];for i = 1:nPopparticle(i).Position = rand(1, nVar);particle(i).Velocity = zeros(1, nVar);particle(i).Cost = ZDT1(particle(i).Position);particle(i).Best.Position = particle(i).Position;particle(i).Best.Cost = particle(i).Cost;
end% ---------------- 初始化外部存档 ----------------
Rep = GetNonDominated(particle);
Rep = ReduceArchive(Rep, nRep);% ---------------- 主循环 ----------------
for it = 1:MaxItfor i = 1:nPopleader = Rep(randi(length(Rep)));% PSO更新particle(i).Velocity = w*particle(i).Velocity ...+ c1*rand(1,nVar).*(particle(i).Best.Position - particle(i).Position) ...+ c2*rand(1,nVar).*(leader.Position - particle(i).Position);particle(i).Position = particle(i).Position + particle(i).Velocity;particle(i).Position = min(max(particle(i).Position, VarMin), VarMax);% 变异if rand < muk = randi([1, nVar]);particle(i).Position(k) = rand();end% 评价目标函数particle(i).Cost = ZDT1(particle(i).Position);% 更新个体最优if Dominates(particle(i).Cost, particle(i).Best.Cost)particle(i).Best.Position = particle(i).Position;particle(i).Best.Cost = particle(i).Cost;endend% 更新外部存档Rep = [Rep, particle];Rep = GetNonDominated(Rep);Rep = ReduceArchive(Rep, nRep);% 可视化Costs = reshape([Rep.Cost], 2, [])';% 绘制理论 Pareto 前沿f1_theory = linspace(0, 1, 100);f2_theory = 1 - sqrt(f1_theory);plot(f1_theory, f2_theory, 'b-', 'LineWidth', 2); hold on;% 绘制当前粒子 Pareto 解plot(Costs(:,1), Costs(:,2), 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r');% 图形设置title(['MOPSO 第 ' num2str(it) ' 代']);xlabel('f_1'); ylabel('f_2'); grid on;legend('理论 Pareto 前沿','MOPSO 非支配解','Location','SouthWest');axis([0 1 0 1]); drawnow; hold off;end% ---------------- 多目标函数：ZDT1 ----------------
function y = ZDT1(x)f1 = x(1);g = 1 + 9 * sum(x(2:end)) / (length(x)-1);f2 = g * (1 - sqrt(f1/g));y = [f1 f2];
end% ---------------- 支配关系判断 ----------------
function b = Dominates(x, y)b = all(x <= y) && any(x < y);
end% ---------------- 非支配筛选 ----------------
function nd = GetNonDominated(pop)n = length(pop);keep = true(1, n);for i = 1:nfor j = 1:nif i ~= j && Dominates(pop(j).Cost, pop(i).Cost)keep(i) = false;break;endendendnd = pop(keep);
end% ---------------- 拥挤度计算并裁剪 ----------------
function pop = ReduceArchive(pop, N)if length(pop) <= Nreturn;endCosts = reshape([pop.Cost], 2, [])';D = zeros(length(pop));for i = 1:length(pop)for j = i+1:length(pop)D(i,j) = norm(Costs(i,:) - Costs(j,:));D(j,i) = D(i,j);endendD(1:size(D,1)+1:end) = inf;crowd = mean(D, 2);[~, idx] = sort(crowd, 'descend');pop = pop(idx(1:N));
end

运行结果

以下从多目标粒子群优化（MOPSO）的数学原理出发，结合代码逐部分拆解，帮你理解 MOPSO 如何优化 ZDT1 多目标函数：

---

代码分析

相关引用：粒子群优化算法（Particle Swarm Optimization, PSO） 求解二维 Rastrigin 函数最小值问题

一、多目标优化的核心概念

多目标优化（MOP）需同时优化多个冲突目标，通常不存在单一最优解，而是寻找Pareto 前沿（非支配解的集合）。

1. Pareto 支配（Dominance）

• 若解 $\mathbf{x}$ 的所有目标值不差于解 $\mathbf{y}$，且至少一个目标值严格更优，则称 $\mathbf{x}$ 支配 $\mathbf{y}$，记为 $\mathbf{x}\prec \mathbf{y}$。
• 数学定义（代码中 Dominates 函数）：
  $$\text{Dominates}(\mathbf{x},\mathbf{y})⟺\left(\forall i,x_i\le y_i\right)\wedge \left(\exists i,x_i<y_i\right)$$

2. 非支配解（Non-dominated Solution）

不被任何其他解支配的解，构成Pareto 前沿（理论最优边界）。

---

二、ZDT1 多目标函数（测试基准）

ZDT1 是经典多目标测试函数，数学形式：
$$\{\begin{array}{l}{f}_1(\mathbf{x})=x_1\\ f_2(\mathbf{x})=g(\mathbf{x})\cdot \left(1-\sqrt{\frac{f_1(\mathbf{x})}{g(\mathbf{x})}}\right)\\ g(\mathbf{x})=1+9\cdot \frac{{\sum }_{i=2}^n{x}_i}{n-1}\end{array}$$

• $\mathbf{x}\in [0,1{]}^n$（$n$ 为决策变量个数，代码中 $n=30$）。
• 理论 Pareto 前沿：当 $g(\mathbf{x})=1$（即 $x_2=\cdots =x_n=0$），此时 $f_2=1-\sqrt{f_1}$，对应代码中 f2_theory = 1 - sqrt(f1_theory)。

---

三、MOPSO 算法的核心数学模型

MOPSO 扩展了单目标 PSO，引入外部存档（Repository保存非支配解，并通过拥挤度（Crowding Distance维护存档多样性。

1. 粒子状态更新（类似单目标 PSO，但目标是向量）

粒子速度更新公式：
$${\mathbf{v}}_i^{(t+1)}=w\cdot {\mathbf{v}}_i^{(t)}+c_1\cdot r_1\cdot \left({\text{pbest}}_i-{\mathbf{x}}_i^{(t)}\right)+c_2\cdot r_2\cdot \left(\text{lbest}-{\mathbf{x}}_i^{(t)}\right)$$

• $\text{lbest}$：从外部存档中随机选的“引导粒子”（替代单目标的 gbest，保持多样性）。

2. 外部存档（Repository）

• 保存迭代中发现的所有非支配解，是 Pareto 前沿的近似。
• 通过拥挤度裁剪（ReduceArchive）控制存档大小，避免爆炸。

---

四、代码解析

（1）参数设置

nVar = 30;         % 决策变量个数（ZDT1 的维度，x ∈ ℝ^30）
VarMin = 0; VarMax = 1; % 变量范围（ZDT1 要求 x_i ∈ [0,1]）
MaxIt = 100;       % 最大迭代次数
nPop = 50;         % 粒子群规模
nRep = 50;         % 外部存档最大容量w = 0.5;           % 惯性权重（平衡探索与利用）
c1 = 1.5; c2 = 2.0;% 学习因子（个体/群体引导）
mu = 0.1;          % 变异概率（增加多样性）

（2）粒子初始化

particle(nPop).Position = []; % 粒子位置（决策变量）
particle(nPop).Velocity = []; % 粒子速度
particle(nPop).Cost = [];     % 目标函数值（2 维向量 [f1,f2]）
particle(nPop).Best.Position = []; % 个体最优位置
particle(nPop).Best.Cost = [];     % 个体最优目标值for i = 1:nPop% 随机初始化位置：x_i ∈ [0,1]particle(i).Position = rand(1, nVar); particle(i).Velocity = zeros(1, nVar); % 速度初始化为 0% 计算目标函数（ZDT1）particle(i).Cost = ZDT1(particle(i).Position); % 个体最优初始化为当前位置particle(i).Best.Position = particle(i).Position; particle(i).Best.Cost = particle(i).Cost;
end

• 数学对应：
  位置初始化：${\mathbf{x}}_i\sim \text{Uniform}(0,1)$
  目标函数：$\text{Cost}({\mathbf{x}}_i)=\text{ZDT1}({\mathbf{x}}_i)=[f_1({\mathbf{x}}_i),f_2({\mathbf{x}}_i)]$

（3）外部存档初始化

% 筛选粒子群中的非支配解
Rep = GetNonDominated(particle); 
% 裁剪存档到 nRep 大小（保持多样性）
Rep = ReduceArchive(Rep, nRep); 

• GetNonDominated：遍历所有粒子，保留不被任何其他粒子支配的解（非支配解）。
• ReduceArchive：若存档超过 nRep，用拥挤度裁剪（保留分布均匀的解）。

（4）主循环（MOPSO 核心迭代）

for it = 1:MaxItfor i = 1:nPop% 从外部存档随机选引导粒子（lbest）leader = Rep(randi(length(Rep))); % 速度更新（PSO 公式）particle(i).Velocity = w*particle(i).Velocity ...+ c1*rand(1,nVar).*(particle(i).Best.Position - particle(i).Position) ...+ c2*rand(1,nVar).*(leader.Position - particle(i).Position);% 位置更新（限制在 [0,1] 内）particle(i).Position = particle(i).Position + particle(i).Velocity; particle(i).Position = min(max(particle(i).Position, VarMin), VarMax);% 变异（增加多样性，避免早熟收敛）if rand < muk = randi([1, nVar]); % 随机选一个维度particle(i).Position(k) = rand(); % 随机重置该维度end% 计算新位置的目标函数particle(i).Cost = ZDT1(particle(i).Position);% 更新个体最优（若新解支配历史最优）if Dominates(particle(i).Cost, particle(i).Best.Cost)particle(i).Best.Position = particle(i).Position;particle(i).Best.Cost = particle(i).Cost;endend% 更新外部存档：合并粒子群和原存档，筛选非支配解，再裁剪Rep = [Rep, particle];        % 合并Rep = GetNonDominated(Rep);   % 保留非支配解Rep = ReduceArchive(Rep, nRep); % 裁剪到 nRep 大小% 可视化（绘制理论前沿和当前非支配解）Costs = reshape([Rep.Cost], 2, [])';f1_theory = linspace(0, 1, 100);f2_theory = 1 - sqrt(f1_theory);plot(f1_theory, f2_theory, 'b-', 'LineWidth', 2); hold on;plot(Costs(:,1), Costs(:,2), 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r');title(['MOPSO 第 ' num2str(it) ' 代']);xlabel('f_1'); ylabel('f_2'); grid on;legend('理论 Pareto 前沿','MOPSO 非支配解','Location','SouthWest');axis([0 1 0 1]); drawnow; hold off;
end

• 核心逻辑：
  1. 速度/位置更新：与单目标 PSO 类似，但用外部存档的随机粒子替代 gbest，避免过早收敛到局部前沿。
  2. 变异：随机重置粒子的一个维度，增加种群多样性。
  3. 存档维护：不断纳入新解，筛选非支配解并裁剪，保证存档是 Pareto 前沿的均匀近似。

（5）关键函数解析

① ZDT1 函数

function y = ZDT1(x)f1 = x(1); % 第一个目标：f1 = x1% 计算 g(x)：1 + 9*(x2+...+xn)/(n-1)g = 1 + 9 * sum(x(2:end)) / (length(x)-1); % 第二个目标：f2 = g*(1 - sqrt(f1/g))f2 = g * (1 - sqrt(f1/g)); y = [f1 f2]; % 返回 2 维目标向量
end

• 数学对应：严格实现 ZDT1 的公式，注意 $g(\mathbf{x})$ 是除 $x_1$ 外其他变量的线性组合。

② 支配关系判断

function b = Dominates(x, y)% x 支配 y 的条件：x 的所有目标 ≤ y，且至少一个 < yb = all(x <= y) && any(x < y); 
end

• 数学对应：直接翻译 Pareto 支配的定义。

③ 非支配筛选

function nd = GetNonDominated(pop)n = length(pop);keep = true(1, n); % 标记是否保留for i = 1:nfor j = 1:n% 若存在 j 支配 i，则 i 不保留if i ~= j && Dominates(pop(j).Cost, pop(i).Cost)keep(i) = false;break;endendendnd = pop(keep); % 返回所有非支配解
end

• 逻辑：双重循环遍历所有粒子，标记不被任何其他粒子支配的解。

④ 拥挤度裁剪

function pop = ReduceArchive(pop, N)if length(pop) <= Nreturn; % 数量不足，直接返回end% 提取所有目标值（2 维）Costs = reshape([pop.Cost], 2, [])'; % 初始化拥挤度矩阵D = zeros(length(pop)); %初始化一个用于存储目标空间中解之间距离的方阵(N*N)for i = 1:length(pop)for j = i+1:length(pop)% 计算解 i 和 j 的欧氏距离（目标空间）D(i,j) = norm(Costs(i,:) - Costs(j,:)); D(j,i) = D(i,j); % 对称矩阵endendD(1:size(D,1)+1:end) = inf; % 对角线设为无穷大（自身距离不计）crowd = mean(D, 2); % 每个解的平均距离（拥挤度：值越大，周围越空）% 按拥挤度降序排序，保留前 N 个（优先保留稀疏区域的解）[~, idx] = sort(crowd, 'descend'); %拥挤度越大越稀疏pop = pop(idx(1:N));
end

reshape(A, m, n)，表示将数组 A 重塑为 m 行 n 列的矩阵（元素总数不变）
reshape([pop.Cost], 2, [])，2 表示目标函数的数量（双目标），[] 表示让 MATLAB 自动计算列数（列数 = 总元素数 ÷ 行数 = 2×nPop ÷ 2 = nPop）
mean(D, 2) 的作用是计算距离矩阵 $\mathbf{D}$中每一行元素的平均值
mean(D, 2) 中的参数 2 表示按行计算均值（MATLAB 中 mean(A, dim) 表示沿第 dim 维求均值）。由于对角线元素为 $\infty$，计算时自动排除自身（$j=i$），因此分母为 $N-1$。

• 数学逻辑：
  • 拥挤度：衡量解在目标空间的“稀疏程度”，距离越远（crowd 越大），说明该解周围的解越少，需要保留以维持多样性。
  • 裁剪：保留拥挤度高的解，保证 Pareto 前沿的均匀分布。

---

五、MOPSO 优化 ZDT1 的关键逻辑

1. 平衡多样性与收敛性：

   • 通过外部存档保存非支配解，跟踪 Pareto 前沿。
   • 通过拥挤度裁剪保证存档解分布均匀，避免陷入局部前沿。

2. 处理多目标冲突：

   • 用“支配关系”替代单目标的“适应度值”，不要求解严格优于所有其他解，而是保留“互不支配”的解集合。

3. 应对 ZDT1 的特性：

   • ZDT1 的 Pareto 前沿是连续曲线 $f_2=1-\sqrt{f_1}$，MOPSO 通过粒子迭代和存档维护，逐步逼近该曲线。

六、结果分析

• 可视化：每代绘制“理论 Pareto 前沿”和“MOPSO 非支配解”，理想情况下，红色点会逐渐覆盖蓝色曲线，形成均匀分布的近似前沿。
• 存档意义：外部存档 Rep 最终保存的非支配解，是多目标优化的“折中解集合”，可根据实际需求（如侧重 $f_1$ 或 $f_2$）选择。

代码展示了如何用 MOPSO 的数学模型优化 ZDT1 函数，以及每个步骤如何服务于“寻找 Pareto 前沿”的目标。MOPSO 的核心是群体协作 + 存档维护，适合解决多目标优化问题。