目录
原理简介
一、Pareto最优概念
二、单目标GWO
三、多目标GWO优化机制
四、整体算法流程
代码实现
今天为大家带来一期多目标灰狼算法（MOGWO）代码，该算法由 Seyedali Mirjalili 等人于 2016 年发表在SCI一区顶刊《Expert Systems With Applications》上！
目前，MOGWO已经被广泛应用于如能源系统优化、物流路径优化、参数调优等不同场景，相比于单目标算法，多目标算法考虑的内容更多，更容易受到审稿人的青睐。
本期代码免费赠送，需要代码的小伙伴可直接拉到最后！
原理简介
一、Pareto最优概念
多目标灰狼优化算法 (Multi-objective Grey Wolf Optimizer, MOGWO)是灰狼优化算法（GWO）的多目标版本，旨在解决多准则下无法比较多目标空间中解的优劣问题，因此引入了Pareto最优解集的概念。
以最小化为例，解A对解B在某个目标函数上存在f(A)<f(B)，则称解A支配解B。在解集内，找不到其他解在所有目标函数上都优于解A的解，则解A为Pareto最优解，这一类解组成的集合为Pareto最优解集，而Pareto前沿则由Pareto最优解的目标函数值组成。
二、单目标GWO
为了更好地了解MOGWO，首先介绍一下单目标GWO原理。单目标GWO通过模仿灰狼捕猎行为进行寻优，其数学模型如下：
式中：t为当前迭代次数；Xα、Xβ和Xδ为阿尔法狼、贝塔狼和德尔塔狼的位置向量；X为灰狼的位置向量；A和C是系数向量。其计算如下：
式中：a在迭代过程中线性地从2减少到0；r1和r2则是[0,1]的随机向量。
三、多目标GWO优化机制
不同于传统的单目标算法，多目标算法能够通过寻找帕累托解平衡多个相互竞争的目标。而MOGWO相比于其他多目标算法，则有两个较为明显的改进，一是引入存档机制，二是改进头狼选择方式。
第一，存档机制。外部存档Archive保存到目前为止获得Pareto最优解，在迭代中新得到的非支配解与存档中的常驻解采用以下处理方式：
（1）新个体被至少以一个存档中的常驻解支配时，新个体不被允许进入存档。
（2）新个体支配存档中的一个或多个解时，新个体进入存档，存档内被支配的解则被省略掉。
（3）如果新个体与存档内的解都不相互支配，则应将新个体加入存档。
（4）当存档已满时，运行网格机制重新安排目标空间的分割，去掉最拥挤的部分的一个解，将新解插入到最不拥挤的位置，以提高Pareto前沿的多样性。
第二，改进头狼选择方式。为选择出合适的三匹头狼（α狼、β狼、δ狼），通过轮盘赌法在Archive中最不拥挤的部分按照如下概率选择头狼：
式中：c为大于1的常数；Ni为该第i组中Pareto最优解个数。
四、整体算法流程
MOGWO的具体流程如下：
（1）设置算法的种群数量、最大迭代次数，设置外部存档Archive大小、轮盘赌法参数等。
（2）计算种群个体的目标参数值，确定支配关系，将非支配解存入Archive中。
（3）根据外部存档中的拥挤度，依据轮盘赌法确定头狼（α狼、β狼、δ狼）。
（4）利用得到的头狼更新种群个体位置并计算目标函数值。
（5）比较新的种群个体与存档中个体的支配关系，确定新的非支配解更新存档。
（6）对步骤（3）、步骤（4）和步骤（5）迭代运行，达到迭代上限停止，输出Archive解。
代码实现
MOGWO核心代码如下：
clear
clc
drawing_flag = 1;
nVar=5;
%% 测试函数
fobj=@(x) ZDT3(x);
%% MOGWO算法参数
lb=zeros(1,5);
ub=ones(1,5);
VarSize=[1 nVar];
GreyWolves_num=100; % 种群数量
MaxIt=50; % 迭代次数
Archive_size=100; % 存档数量
%% 网格机制的参数
alpha=0.1; % Grid Inflation Parameter
nGrid=10; % Number of Grids per each Dimension
beta=4; % Leader Selection Pressure Parameter
gamma=2; % Extra (to be deleted) Repository Member Selection Pressure
%% 种群初始化
GreyWolves=CreateEmptyParticle(GreyWolves_num);
for i=1:GreyWolves_num
GreyWolves(i).Velocity=0;
GreyWolves(i).Position=zeros(1,nVar);
for j=1:nVar
GreyWolves(i).Position(1,j)=unifrnd(lb(j),ub(j),1);
end
GreyWolves(i).Cost=fobj(GreyWolves(i).Position')';
GreyWolves(i).Best.Position=GreyWolves(i).Position;
GreyWolves(i).Best.Cost=GreyWolves(i).Cost;
end
%% 确定支配关系
GreyWolves=DetermineDomination(GreyWolves);
%% 非支配解存档
Archive=GetNonDominatedParticles(GreyWolves);
%% 网格机制
Archive_costs=GetCosts(Archive);
G=CreateHypercubes(Archive_costs,nGrid,alpha);
for i=1:numel(Archive)
[Archive(i).GridIndex Archive(i).GridSubIndex]=GetGridIndex(Archive(i),G);
end
%% 主程序迭代
for it=1:MaxIt
a=2-it*((2)/MaxIt);
for i=1:GreyWolves_num
clear rep2
clear rep3
% Choose the alpha, beta, and delta grey wolves
Delta=SelectLeader(Archive,beta);
Beta=SelectLeader(Archive,beta);
Alpha=SelectLeader(Archive,beta);
% If there are less than three solutions in the least crowded
% hypercube, the second least crowded hypercube is also found
% to choose other leaders from.
if size(Archive,1)>1
counter=0;
for newi=1:size(Archive,1)
if sum(Delta.Position~=Archive(newi).Position)~=0
counter=counter+1;
rep2(counter,1)=Archive(newi);
end
end
Beta=SelectLeader(rep2,beta);
end
% This scenario is the same if the second least crowded hypercube
% has one solution, so the delta leader should be chosen from the
% third least crowded hypercube.
if size(Archive,1)>2
counter=0;
for newi=1:size(rep2,1)
if sum(Beta.Position~=rep2(newi).Position)~=0
counter=counter+1;
rep3(counter,1)=rep2(newi);
end
end
Alpha=SelectLeader(rep3,beta);
end
% Eq.(3.4) in the paper
c=2.*rand(1, nVar);
% Eq.(3.1) in the paper
D=abs(c.*Delta.Position-GreyWolves(i).Position);
% Eq.(3.3) in the paper
A=2.*a.*rand(1, nVar)-a;
% Eq.(3.8) in the paper
X1=Delta.Position-A.*abs(D);
% Eq.(3.4) in the paper
c=2.*rand(1, nVar);
% Eq.(3.1) in the paper
D=abs(c.*Beta.Position-GreyWolves(i).Position);
% Eq.(3.3) in the paper
A=2.*a.*rand()-a;
% Eq.(3.9) in the paper
X2=Beta.Position-A.*abs(D);
% Eq.(3.4) in the paper
c=2.*rand(1, nVar);
% Eq.(3.1) in the paper
D=abs(c.*Alpha.Position-GreyWolves(i).Position);
% Eq.(3.3) in the paper
A=2.*a.*rand()-a;
% Eq.(3.10) in the paper
X3=Alpha.Position-A.*abs(D);
% Eq.(3.11) in the paper
GreyWolves(i).Position=(X1+X2+X3)./3;
% Boundary checking
GreyWolves(i).Position=min(max(GreyWolves(i).Position,lb),ub);
GreyWolves(i).Cost=fobj(GreyWolves(i).Position')';
end
GreyWolves=DetermineDomination(GreyWolves);
non_dominated_wolves=GetNonDominatedParticles(GreyWolves);
Archive=[Archive
non_dominated_wolves];
Archive=DetermineDomination(Archive);
Archive=GetNonDominatedParticles(Archive);
for i=1:numel(Archive)
[Archive(i).GridIndex Archive(i).GridSubIndex]=GetGridIndex(Archive(i),G);
end
if numel(Archive)>Archive_size
EXTRA=numel(Archive)-Archive_size;
Archive=DeleteFromRep(Archive,EXTRA,gamma);
Archive_costs=GetCosts(Archive);
G=CreateHypercubes(Archive_costs,nGrid,alpha);
end
disp(['In iteration ' num2str(it) ': Number of solutions in the archive = ' num2str(numel(Archive))]);
save results
% Results
costs=GetCosts(GreyWolves);
Archive_costs=GetCosts(Archive);
if drawing_flag==1
hold off
plot(costs(1,:),costs(2,:),'k.');
hold on
plot(Archive_costs(1,:),Archive_costs(2,:),'r*');
legend('灰狼种群','非支配解');
set(gcf,'color','w')
drawnow
end
end代码里提供了四种多目标函数，分别为ZDT1、ZDT2、ZDT3、ZDT4，大家可以自行切换，以ZDT3为例：
这是迭代过程图，图中可以很清晰的显示灰狼种群与各非支配解，在迭代完成后选择需要的非支配解即可。
