多目标灰狼算法(MOGWO)的Matlab代码详细注释及难点解释(佳点集改进初始种群的MOGWO)

一、外部种群Archive机制

二、领导者选择机制

三、多目标灰狼算法运行步骤

四、MOGWO的Matlab部分代码详细注释

五、MOGWO算法难点解释

5.1 网格与膨胀因子

5.2 轮盘赌方法选择每个超立方体概率

为了将灰狼算法应用于多目标优化问题,在灰狼算法中引入外部种群Archive，用于存储非支配最优解。采用领导者选择策略，从外部种群Archive选择捕食过程中的领导者 $eq?%5Calpha$ 狼、 $eq?%5Cbeta$ 狼及 $eq?%5Cdelta$ 狼。

一、外部种群Archive机制

算法每次迭代会产生新的个体，当个体要加入外部种群Archive时，将这些个体逐一与Archive中的个体进行比较，将会出现4种情况：

(1)如果新个体被至少一个存档中的个体支配，则该个体不加入存档。 (2)如果新个体支配存档中的一个或多个个体，则新个体加入存档，同时将被其支配的个体从存档中删除。 (3)如果新个体与存档中的任一个体互不支配，则将该个体加入存档。 (4)如果存档中已满,应首先运行网格机制，重新安排目标空间的划分，并在最拥挤组中随机删除某些个体，新的个体被插入到不拥挤组中，提高接近最优前沿的多样性。

二、领导者选择机制

在存档中存放所有的非支配最优解，采用轮盘赌的方式从存档中选取头狼(包括 $eq?%5Calpha$ 狼、 $eq?%5Cbeta$ 狼及 $eq?%5Cdelta$ 狼)，选择搜索空间中最不拥挤的区域，并给出其非支配解中的一个作为 $eq?%5Calpha$ 狼、 $eq?%5Cbeta$ 狼或 $eq?%5Cdelta$ 狼。为了提高算法的探索能力，每一个个体被选择的概率与其所在组的个体数成反比。

如果在最不拥挤的超立方体中少于三个解，也会找到第二个不拥挤的超立方体来选择其它领导者，暂时将Delta排除在存档之外，以避免选择相同领导者；如果第二个最不拥挤的超立方体有一个解，这种情况也是一样的，因此应该从第三个最不拥挤的超立方体，选择其它领导者，暂时将Beta排除在存档之外，以避免选择相同领导者。

三、多目标灰狼算法运行步骤

四、MOGWO的Matlab部分代码详细注释

完整注释代码私信获取

clear all
clc
drawing_flag = 1;
nVar=5;        %变量维数
fobj=@(x) ZDT1(x);
lb=zeros(1,5); %下限
ub=ones(1,5);  %上限
VarSize=[1 nVar];
GreyWolves_num=100; %种群规模
for i=1:GreyWolves_num
    for j=1:nVar
        p1=2*nVar+3;
        r=2*cos(2*pi*j/p1);
        G(i,j)=mod(r*i,1);
    end
end
Boundary_no= size(ub,2); % numnber of boundaries

% If the boundaries of all variables are equal and user enter a signle
% number for both ub and lb
if Boundary_no==1
    X=G.*(ub-lb)+lb;
end

% If each variable has a different lb and ub
if Boundary_no>1
    for i=1:nVar
        ub_i=ub(i);
        lb_i=lb(i);
        X(:,i)=G(:,i).*(ub_i-lb_i)+lb_i;
    end
end

MaxIt=100;  % 最大迭代次数
Archive_size=100;   % 存档规模
alpha=0.1;  % 网格膨胀因子
nGrid=10;   % 每一维上的网格数
beta=4;     % 领导者选择压力因子
gamma=2;    % 额外的（待删除）存档个体选择压力
% 初始化
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)     %对于第i个个体，第j个目标
    [Archive(i).GridIndex, Archive(i).GridSubIndex]=GetGridIndex(Archive(i),G);
end


%% MOGWO main loop
for it=1:MaxIt
    a=2-it*((2)/MaxIt);     %a从2线性减小到0
    for i=1:GreyWolves_num  %对于每个搜索个体     
        clear rep2
        clear rep3       
        % Choose the alpha, beta, and delta grey wolves
        %领导者选择组件选择搜索空间中最不拥挤的区域，并给出其非支配解中一个作为alpha,、beta或delta狼

        Delta=SelectLeader(Archive,beta);
        Beta=SelectLeader(Archive,beta);
        Alpha=SelectLeader(Archive,beta);  

        % 如果在最不拥挤的超立方体中少于三个解，也会找到第二个不拥挤的超立方体来选择其它领导者
       
        if size(Archive,1)>1           %如果存档中不止一个个体
            counter=0;
            for newi=1:size(Archive,1) %对于存档中每个个体
                if sum(Delta.Position~=Archive(newi).Position)~=0 %如果Delta是存档中的个体
                    counter=counter+1;
                    rep2(counter,1)=Archive(newi);                %rep2是排除Delta后的存档
                end
            end
            Beta=SelectLeader(rep2,beta); %在排除Delta后的存档中用轮盘赌法选择Beta
        end  

.
.
.
.
.
.
.

function [occ_cell_index ,occ_cell_member_count]=GetOccupiedCells(pop)
    GridIndices=[pop.GridIndex];         %存档中每个个体在网格中的位置
    occ_cell_index=unique(GridIndices);  %GridIndices剔除一个小网格中的重复值，从小到大排序 
    occ_cell_member_count=zeros(size(occ_cell_index));
    m=numel(occ_cell_index);             %存档中个体占据的小网格数
    for k=1:m
        occ_cell_member_count(k)=sum(GridIndices==occ_cell_index(k)); %计算小网格上解数量
    end    
end


function i=RouletteWheelSelection(p)
    r=rand;
    c=cumsum(p);            %累计概率
    i=find(r<=c,1,'first'); %r首次<=c的第几个值
end


function z=ZDT1(x)
n=numel(x);
f1=x(1);
g=1+9/(n-1)*sum(x(2:end));
h=1-sqrt(f1/g);
f2=g*h;
z=[f1
    f2];
end

function z=ZDT2(x)
n=numel(x);
f1=x(1);
g=1+9/(n-1)*sum(x(2:end));
h=1-(f1/g)^2;
f2=g*h;
z=[f1
    f2];
end

function z=ZDT3(x)
n=numel(x);
f1=x(1);
g=1+9/(n-1)*sum(x(2:end));
h=1-f1/g-(f1/g)*sin(10*pi*x(1));
f2=g*h;
z=[f1
    f2];
end

function z=ZDT4(x)
n=numel(x);
f1=x(1);
sum=0;
for i=2:n
    sum = sum+(x(i)^2-10*cos(4*pi*x(i)));
end
g=1+(n-1)*10+sum;
f2=g*(1-(f1/g)^0.5);
z=[f1
    f2];
end

五、MOGWO算法难点解释

5.1 网格与膨胀因子

算法中的网格其实就是将目标空间划分为一个个网格，划分网格后就可以知道每个网格中有几个粒子，同理也就能知道某个粒子在哪个网格中。网格膨胀因子就是将网格膨胀的参数。图一中第一坐标外的就是膨胀出的网格。 $eq?mincj$ 就是膨胀后的网格下边界， $eq?maxcj$ 就是膨胀后的网格上边界。

图1 网格膨胀因子解释图

 mincj=min(costs(j,:));  %每个目标上的存档个体目标函数最小值
 maxcj=max(costs(j,:));  %每个目标上的存档个体目标函数的最大值      
 dcj=alpha*(maxcj-mincj);      
 mincj=mincj-dcj;          
 maxcj=maxcj+dcj;

5.2 轮盘赌方法选择每个超立方体概率

 beta=4;     % 领导者选择压力因子
.
.
.
 p=occ_cell_member_count.^(-beta);
 p=p/sum(p); %通过轮盘赌方法进行选择，每个超立方体概率 
 selected_cell_index=occ_cell_index(RouletteWheelSelection(p)); %用轮盘赌选择的所在区域 
.
.
.
 function i=RouletteWheelSelection(p)
    r=rand;
    c=cumsum(p);            %累计概率
    i=find(r<=c,1,'first'); %r首次<=c的第几个值
 end

例: occ_cell_member_count=[2 1 1 1 1 1 1]，则 p = [0.0103 0.1649 0.1649 0.1649 0.1649 0.1649 0.1639]，因为这里 p = occ_cell_member_count.^(-beta))。则累计概率c = [0.0103 0.1753 0.3402 0.5052 0.6701 0.8351 1]。因为这里 $eq?c%3Dcumsum%5Cleft%20%28%20p%20%5Cright%20%29$ 。

如果运行后随机数 $eq?r%3Drand%3D0.3828$ ，则 $eq?r$ 首次 $eq?%5Cleq%20c$ 时位置为4，所以 $eq?i%3D4$ 。选择的区域selected_cell_index为occ_cell_index数组中第4个位置的区域。这里假设occ_cell_index=[17 19 26 33 43 62 82]，就有selected_cell_index=33。