【Matlab】智能优化算法_灰狼优化算法GWO

1.背景介绍

灰狼算法的提出背景是受到了灰狼群体捕猎行为的启发。灰狼是一种高度社会化的动物，它们有着严格的等级制度和协作机制3。灰狼算法模拟了自然界灰狼的领导层级和狩猎机制，利用四种类型的灰狼（α、β、δ、ω）来模拟领导阶层，并实现了寻找猎物、包围猎物和攻击猎物的三个主要步骤。
灰狼优化算法（Grey Wolf Optimizer，GWO）是由自然界中灰狼群体的社会等级机制和捕猎行为而衍生出来的一种群体优化智能算法，目前已成功运用到车间调度、参数优化、图像分类等领域中。

2.基本思想

2.1 等级制度

1. 灰狼优化算法的基本思想是模仿灰狼的狩猎和社会等级行为，利用一小部分拥有绝对话语权的灰狼带领一群灰狼向猎物前进。灰狼优化算法中，每只灰狼的位置代表一个可行解，而最优解就是最佳的猎物位置。灰狼优化算法模拟了自然界灰狼的领导层级和狩猎机制，分为四种类型的灰狼（α、β、δ、ω），并实现了寻找猎物、包围猎物和攻击猎物的三个主要步骤。
2. 灰狼（犬狼疮）属于犬科。灰狼被认为是顶级捕食者，这意味着它们处于食物链的顶端。灰狼大多喜欢群居。平均人数为5-12人。特别令人感兴趣的是，他们有一个非常严格的社会支配等级制度，如图1所示。

1. 灰狼群一般分为4个等级：领导者是一男一女，称为α。α主要负责决定狩猎、睡觉地点、醒来时间等。α的决定是由狼群决定的。然而，也观察到了某种民主行为，其中α跟随狼群中的其他狼。在聚会中，整个人都会低头认输。α狼也被称为优势狼，因为他/她的命令应该由狼群跟随。α狼只允许在狼群中交配。有趣的是，α不一定是团队中最强的成员，但在管理团队方面却是最好的。这表明，一个团队的组织和纪律比其力量更重要。
2. 灰狼等级体系中的第二级是β。β狼是从属的狼，在决策或其他狼群活动中帮助α狼。β狼可以是雄性或雌性，如果α狼中的一只去世或变得很老，他/她可能是成为α狼的最佳人选。β狼应该尊重阿尔法狼，但也指挥其他较低级别的狼。它扮演着α的顾问和团队的惩戒者的角色。
3. 处于第三阶段的灰狼用δ表示，侦察狼、守卫狼、老狼和捕食狼都是这一类。
4. 排名最低的灰狼是ω。ω扮演了替罪羊的角色。ω狼总是要屈从于所有其他占优势的狼。它们是最后被允许进食的狼。ω似乎不是一个重要的个体，但据观察，如果失去ω，整个群体都会面临内部斗争和问题。这是由于ω人发泄了所有狼群的暴力和沮丧。这有助于满足整个团队的需求，并保持优势结构。在某些情况下，ω也是团队中的保姆。

2.2 狩猎方式

除了狼的社会等级之外，群体狩猎是灰狼的另一种有趣的社会行为。灰狼狩猎的主要阶段如下：

1. 跟踪、追逐和接近猎物
2. 追逐、包围和骚扰猎物，直到其停止移动
3. 攻击猎物

3.公式推导

依据上述灰狼等级制度，可以对灰狼的社会等级进行数学建模，认为最合适的解是α，那么第二和第三最优解分别表示为β和δ，而剩余其他解都假定为ω。在GWO中，通过α、β和δ来导引捕食（优化），ω听从于这三种狼。在这项工作中，对灰狼的狩猎技术和社会等级进行了数学建模，以设计GWO并进行优化。提供了社会等级、追踪、包围和攻击猎物的数学模型。然后概述了GWO算法。

3.1 社会等级制度

为了在设计GWO时对狼的社会等级进行数学建模，我们将最适合的解决方案视为α。因此，第二和第三最佳解决方案分别命名为β和δ。其余的候选解决方案假定为ω。在GWO算法中，搜索（优化）由α、β和δ指导。

3.2 包围猎物

如上所述，灰狼在狩猎过程中包围猎物。为了对环绕行为进行数学建模，提出了以下方程：

其中t表示当前迭代，A和C是系数向量，Xp猎物的位置向量，X表示灰狼的位置向量。
矢量A和C计算如下

其中的分量a在迭代过程中从2线性减少到0，r1，r2是[0,1]中的随机向量。

为了观察方程（3.1）和（3.2）的效果，图3（a）中说明了二维位置向量和一些可能的邻居。如图所示，处于（X，Y）位置的灰狼可以根据猎物的位置（X*，Y*）更新其位置。通过调整A和C向量的值，可以相对于当前位置到达最佳代理周围的不同位置。例如，通过设置A=（1,0）和C=（1,1）可以达到（X*-X，Y*）。灰狼在3D空间中可能的更新位置如图3（b）所示。注意，随机向量r1，r2允许狼到达图3所示的点之间的任何位置。因此，灰狼可以通过使用等式（3.1）和（3.2）来更新其在任何随机位置的猎物周围空间内的位置。

3.3 包围猎物

灰狼有能力识别猎物的位置并包围它们。狩猎通常由α指导。β和δ也可能偶尔参与狩猎。然而，在抽象的搜索空间中，我们不知道最佳（猎物）的位置。为了数学模拟灰狼的狩猎行为，我们假设α（最佳候选解）β和δ对猎物的潜在位置有更好的了解。因此，我们保存到目前为止获得的前三个最佳解决方案，并强制其他搜索代理（包括omegas）根据最佳搜索代理的位置更新其位置。在这方面提出了以下公式。

图4显示了搜索代理如何根据2D搜索空间中的α，β和δ更新其位置。可以观察到，最终位置将位于由搜索空间中的α、β和δ的位置定义的圆内的随机位置。换句话说，α、β和δ估计猎物的位置，而其他狼则随机更新猎物周围的位置。

3.4 攻击猎物

如上所述，灰狼在猎物停止移动时攻击猎物，从而完成狩猎。为了对接近猎物进行数学建模，我们降低了a的值。注意，A的波动范围也减小了。换句话说，A是区间[-a，a]中的随机值，其中a在迭代过程中从2减少到0。当的随机值A为[-1,1]时，搜索代理的下一个位置可以位于其当前位置和猎物位置之间的任何位置。图5（a）显示|A|<1迫使狼向猎物攻击。

对于目前提出的算子，GWO算法允许其搜索代理基于α、β和δ的位置更新其位置；并攻击猎物。然而，GWO算法在这些算子的局部解中容易停滞。的确，所提出的包围机制在某种程度上显示了探索，但GWO需要更多的运营商来强调探索。

3.5 搜索猎物

灰狼主要根据α、β和δ的位置进行搜索。它们彼此分开寻找猎物，并会聚攻击猎物。为了对发散进行数学建模，我们利用A大于1或小于-1的随机值迫使搜索代理偏离猎物。这强调了探索，并允许GWO算法进行全局搜索。图5（b）还显示|A|>1迫使灰狼偏离猎物，希望找到更合适的猎物。GWO的另一个有利于勘探的组成部分是C。如方程（3.4）所示，向量C包含[0，2]中的随机值。该组件为猎物提供随机权重，以便在定义方程（3.1）中的距离时随机强调（C>1）或弱化（C<1）猎物的影响。这有助于GWO在整个优化过程中表现出更随机的行为，有利于探索和避免局部最优。这里值得一提的是，与A相比，C并不是线性减少的。我们故意要求C始终提供随机值，以便不仅在初始迭代期间而且在最终迭代期间强调探索。该组件在局部最优停滞的情况下非常有用，特别是在最终迭代中。
C矢量也可以被认为是自然界中接近猎物的障碍物的影响。一般来说，自然界中的障碍出现在狼的狩猎路径上，事实上阻止它们快速方便地接近猎物。这正是向量C所做的。根据狼的位置，它可以随机地给猎物一个重量，使其更难和更远地接近狼，反之亦然。
总之，搜索过程从在GWO算法中创建灰狼（候选解决方案）的随机种群开始。在迭代过程中，α、β和δ狼估计猎物的可能位置。每个候选解决方案都会更新其与猎物的距离。参数a从2减小到0，以分别强调勘探和开采。当|A|>1时，候选解倾向于偏离猎物；当|A|<1时，向猎物汇聚。最后，GWO算法通过满足结束标准而终止。

4.算法流程图

5.文件结构

func_plot.m					% 绘制的基准函数
Get_Functions_details.m		% 基准的全部信息和实现
GWO.m						% 灰狼优化算法
initialization.m			% 初始化
main.m						% 主函数

6.伪代码

7.详细代码及注释

7.1 func_plot.m

% This function draw the benchmark functions

function func_plot(func_name)

[lb,ub,dim,fobj]=Get_Functions_details(func_name);

switch func_name 
    case 'F1' 
        x=-100:2:100; y=x; %[-100,100]
        
    case 'F2' 
        x=-100:2:100; y=x; %[-10,10]
        
    case 'F3' 
        x=-100:2:100; y=x; %[-100,100]
        
    case 'F4' 
        x=-100:2:100; y=x; %[-100,100]
    case 'F5' 
        x=-200:2:200; y=x; %[-5,5]
    case 'F6' 
        x=-100:2:100; y=x; %[-100,100]
    case 'F7' 
        x=-1:0.03:1;  y=x  %[-1,1]
    case 'F8' 
        x=-500:10:500;y=x; %[-500,500]
    case 'F9' 
        x=-5:0.1:5;   y=x; %[-5,5]    
    case 'F10' 
        x=-20:0.5:20; y=x;%[-500,500]
    case 'F11' 
        x=-500:10:500; y=x;%[-0.5,0.5]
    case 'F12' 
        x=-10:0.1:10; y=x;%[-pi,pi]
    case 'F13' 
        x=-5:0.08:5; y=x;%[-3,1]
    case 'F14' 
        x=-100:2:100; y=x;%[-100,100]
    case 'F15' 
        x=-5:0.1:5; y=x;%[-5,5]
    case 'F16' 
        x=-1:0.01:1; y=x;%[-5,5]
    case 'F17' 
        x=-5:0.1:5; y=x;%[-5,5]
    case 'F18' 
        x=-5:0.06:5; y=x;%[-5,5]
    case 'F19' 
        x=-5:0.1:5; y=x;%[-5,5]
    case 'F20' 
        x=-5:0.1:5; y=x;%[-5,5]        
    case 'F21' 
        x=-5:0.1:5; y=x;%[-5,5]
    case 'F22' 
        x=-5:0.1:5; y=x;%[-5,5]     
    case 'F23' 
        x=-5:0.1:5; y=x;%[-5,5]  
end    

    

L=length(x);
f=[];

for i=1:L
    for j=1:L
        if strcmp(func_name,'F15')==0 && strcmp(func_name,'F19')==0 && strcmp(func_name,'F20')==0 && strcmp(func_name,'F21')==0 && strcmp(func_name,'F22')==0 && strcmp(func_name,'F23')==0
            f(i,j)=fobj([x(i),y(j)]);
        end
        if strcmp(func_name,'F15')==1
            f(i,j)=fobj([x(i),y(j),0,0]);
        end
        if strcmp(func_name,'F19')==1
            f(i,j)=fobj([x(i),y(j),0]);
        end
        if strcmp(func_name,'F20')==1
            f(i,j)=fobj([x(i),y(j),0,0,0,0]);
        end       
        if strcmp(func_name,'F21')==1 || strcmp(func_name,'F22')==1 ||strcmp(func_name,'F23')==1
            f(i,j)=fobj([x(i),y(j),0,0]);
        end          
    end
end

surfc(x,y,f,'LineStyle','none');

end

7.2 Get_Functions_details.m

% This function containts full information and implementations of the benchmark 
% functions in Table 1, Table 2, and Table 3 in the paper

% lb is the lower bound: lb=[lb_1,lb_2,...,lb_d]
% up is the uppper bound: ub=[ub_1,ub_2,...,ub_d]
% dim is the number of variables (dimension of the problem)

function [lb,ub,dim,fobj] = Get_Functions_details(F)


switch F
    case 'F1'
        fobj = @F1;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F2'
        fobj = @F2;
        lb=-10;
        ub=10;
        dim=30;
        
    case 'F3'
        fobj = @F3;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F4'
        fobj = @F4;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F5'
        fobj = @F5;
        lb=-30;
        ub=30;
        dim=30;
        
    case 'F6'
        fobj = @F6;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F7'
        fobj = @F7;
        lb=-1.28;
        ub=1.28;
        dim=30;
        
    case 'F8'
        fobj = @F8;
        lb=-500;
        ub=500;
        dim=30;
        
    case 'F9'
        fobj = @F9;
        lb=-5.12;
        ub=5.12;
        dim=30;
        
    case 'F10'
        fobj = @F10;
        lb=-32;
        ub=32;
        dim=30;
        
    case 'F11'
        fobj = @F11;
        lb=-600;
        ub=600;
        dim=30;
        
    case 'F12'
        fobj = @F12;
        lb=-50;
        ub=50;
        dim=30;
        
    case 'F13'
        fobj = @F13;
        lb=-50;
        ub=50;
        dim=30;
        
    case 'F14'
        fobj = @F14;
        lb=-65.536;
        ub=65.536;
        dim=2;
        
    case 'F15'
        fobj = @F15;
        lb=-5;
        ub=5;
        dim=4;
        
    case 'F16'
        fobj = @F16;
        lb=-5;
        ub=5;
        dim=2;
        
    case 'F17'
        fobj = @F17;
        lb=[-5,0];
        ub=[10,15];
        dim=2;
        
    case 'F18'
        fobj = @F18;
        lb=-2;
        ub=2;
        dim=2;
        
    case 'F19'
        fobj = @F19;
        lb=0;
        ub=1;
        dim=3;
        
    case 'F20'
        fobj = @F20;
        lb=0;
        ub=1;
        dim=6;     
        
    case 'F21'
        fobj = @F21;
        lb=0;
        ub=10;
        dim=4;    
        
    case 'F22'
        fobj = @F22;
        lb=0;
        ub=10;
        dim=4;    
        
    case 'F23'
        fobj = @F23;
        lb=0;
        ub=10;
        dim=4;            
end

end

% F1

function o = F1(x)
o=sum(x.^2);
end

% F2

function o = F2(x)
o=sum(abs(x))+prod(abs(x));
end

% F3

function o = F3(x)
dim=size(x,2);
o=0;
for i=1:dim
    o=o+sum(x(1:i))^2;
end
end

% F4

function o = F4(x)
o=max(abs(x));
end

% F5

function o = F5(x)
dim=size(x,2);
o=sum(100*(x(2:dim)-(x(1:dim-1).^2)).^2+(x(1:dim-1)-1).^2);
end

% F6

function o = F6(x)
o=sum(abs((x+.5)).^2);
end

% F7

function o = F7(x)
dim=size(x,2);
o=sum([1:dim].*(x.^4))+rand;
end

% F8

function o = F8(x)
o=sum(-x.*sin(sqrt(abs(x))));
end

% F9

function o = F9(x)
dim=size(x,2);
o=sum(x.^2-10*cos(2*pi.*x))+10*dim;
end

% F10

function o = F10(x)
dim=size(x,2);
o=-20*exp(-.2*sqrt(sum(x.^2)/dim))-exp(sum(cos(2*pi.*x))/dim)+20+exp(1);
end

% F11

function o = F11(x)
dim=size(x,2);
o=sum(x.^2)/4000-prod(cos(x./sqrt([1:dim])))+1;
end

% F12

function o = F12(x)
dim=size(x,2);
o=(pi/dim)*(10*((sin(pi*(1+(x(1)+1)/4)))^2)+sum((((x(1:dim-1)+1)./4).^2).*...
(1+10.*((sin(pi.*(1+(x(2:dim)+1)./4)))).^2))+((x(dim)+1)/4)^2)+sum(Ufun(x,10,100,4));
end

% F13

function o = F13(x)
dim=size(x,2);
o=.1*((sin(3*pi*x(1)))^2+sum((x(1:dim-1)-1).^2.*(1+(sin(3.*pi.*x(2:dim))).^2))+...
((x(dim)-1)^2)*(1+(sin(2*pi*x(dim)))^2))+sum(Ufun(x,5,100,4));
end

% F14

function o = F14(x)
aS=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;,...
-32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];

for j=1:25
    bS(j)=sum((x'-aS(:,j)).^6);
end
o=(1/500+sum(1./([1:25]+bS))).^(-1);
end

% F15

function o = F15(x)
aK=[.1957 .1947 .1735 .16 .0844 .0627 .0456 .0342 .0323 .0235 .0246];
bK=[.25 .5 1 2 4 6 8 10 12 14 16];bK=1./bK;
o=sum((aK-((x(1).*(bK.^2+x(2).*bK))./(bK.^2+x(3).*bK+x(4)))).^2);
end

% F16

function o = F16(x)
o=4*(x(1)^2)-2.1*(x(1)^4)+(x(1)^6)/3+x(1)*x(2)-4*(x(2)^2)+4*(x(2)^4);
end

% F17

function o = F17(x)
o=(x(2)-(x(1)^2)*5.1/(4*(pi^2))+5/pi*x(1)-6)^2+10*(1-1/(8*pi))*cos(x(1))+10;
end

% F18

function o = F18(x)
o=(1+(x(1)+x(2)+1)^2*(19-14*x(1)+3*(x(1)^2)-14*x(2)+6*x(1)*x(2)+3*x(2)^2))*...
    (30+(2*x(1)-3*x(2))^2*(18-32*x(1)+12*(x(1)^2)+48*x(2)-36*x(1)*x(2)+27*(x(2)^2)));
end

% F19

function o = F19(x)
aH=[3 10 30;.1 10 35;3 10 30;.1 10 35];cH=[1 1.2 3 3.2];
pH=[.3689 .117 .2673;.4699 .4387 .747;.1091 .8732 .5547;.03815 .5743 .8828];
o=0;
for i=1:4
    o=o-cH(i)*exp(-(sum(aH(i,:).*((x-pH(i,:)).^2))));
end
end

% F20

function o = F20(x)
aH=[10 3 17 3.5 1.7 8;.05 10 17 .1 8 14;3 3.5 1.7 10 17 8;17 8 .05 10 .1 14];
cH=[1 1.2 3 3.2];
pH=[.1312 .1696 .5569 .0124 .8283 .5886;.2329 .4135 .8307 .3736 .1004 .9991;...
.2348 .1415 .3522 .2883 .3047 .6650;.4047 .8828 .8732 .5743 .1091 .0381];
o=0;
for i=1:4
    o=o-cH(i)*exp(-(sum(aH(i,:).*((x-pH(i,:)).^2))));
end
end

% F21

function o = F21(x)
aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];

o=0;
for i=1:5
    o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end

% F22

function o = F22(x)
aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];

o=0;
for i=1:7
    o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end

% F23

function o = F23(x)
aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];

o=0;
for i=1:10
    o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end

function o=Ufun(x,a,k,m)
o=k.*((x-a).^m).*(x>a)+k.*((-x-a).^m).*(x<(-a));
end

7.3 GWO.m

% Grey Wolf Optimizer
function [Alpha_score,Alpha_pos,Convergence_curve]=GWO(SearchAgents_no,Max_iter,lb,ub,dim,fobj)

% initialize alpha, beta, and delta_pos
Alpha_pos=zeros(1,dim);
Alpha_score=inf; %change this to -inf for maximization problems

Beta_pos=zeros(1,dim);
Beta_score=inf; %change this to -inf for maximization problems

Delta_pos=zeros(1,dim);
Delta_score=inf; %change this to -inf for maximization problems

%Initialize the positions of search agents
Positions=initialization(SearchAgents_no,dim,ub,lb);

Convergence_curve=zeros(1,Max_iter);

l=0;% Loop counter

% Main loop
while l<Max_iter
    for i=1:size(Positions,1)  
        
       % Return back the search agents that go beyond the boundaries of the search space
        Flag4ub=Positions(i,:)>ub;
        Flag4lb=Positions(i,:)<lb;
        Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               
        
        % Calculate objective function for each search agent
        fitness=fobj(Positions(i,:));
        
        % Update Alpha, Beta, and Delta
        if fitness<Alpha_score 
            Alpha_score=fitness; % Update alpha
            Alpha_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness<Beta_score 
            Beta_score=fitness; % Update beta
            Beta_pos=Positions(i,:);
        end
        
        if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score 
            Delta_score=fitness; % Update delta
            Delta_pos=Positions(i,:);
        end
    end
    
    
    a=2-l*((2)/Max_iter); % a decreases linearly fron 2 to 0
    
    % Update the Position of search agents including omegas
    for i=1:size(Positions,1)
        for j=1:size(Positions,2)     
                       
            r1=rand(); % r1 is a random number in [0,1]
            r2=rand(); % r2 is a random number in [0,1]
            
            A1=2*a*r1-a; % Equation (3.3)
            C1=2*r2; % Equation (3.4)
            
            D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); % Equation (3.5)-part 1
            X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1
                       
            r1=rand();
            r2=rand();
            
            A2=2*a*r1-a; % Equation (3.3)
            C2=2*r2; % Equation (3.4)
            
            D_beta=abs(C2*Beta_pos(j)-Positions(i,j)); % Equation (3.5)-part 2
            X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2       
            
            r1=rand();
            r2=rand(); 
            
            A3=2*a*r1-a; % Equation (3.3)
            C3=2*r2; % Equation (3.4)
            
            D_delta=abs(C3*Delta_pos(j)-Positions(i,j)); % Equation (3.5)-part 3
            X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3             
            
            Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7)
            
        end
    end
    l=l+1;    
    Convergence_curve(l)=Alpha_score;
end

7.4 initialization.m

% This function initialize the first population of search agents
function Positions=initialization(SearchAgents_no,dim,ub,lb)

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
    Positions=rand(SearchAgents_no,dim).*(ub-lb)+lb;
end

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

7.5 main.m

clear all 
clc

SearchAgents_no=30; % Number of search agents

Function_name='F2'; % Name of the test function that can be from F1 to F23 (Table 1,2,3 in the paper)

Max_iteration=500; % Maximum numbef of iterations

% Load details of the selected benchmark function
[lb,ub,dim,fobj]=Get_Functions_details(Function_name);

[Best_score,Best_pos,GWO_cg_curve]=GWO(SearchAgents_no,Max_iteration,lb,ub,dim,fobj);
figure('Position',[500 500 660 290])
%Draw search space
subplot(1,2,1);
func_plot(Function_name);
title('Parameter space')
xlabel('x_1');
ylabel('x_2');
zlabel([Function_name,'( x_1 , x_2 )'])

%Draw objective space
subplot(1,2,2);
semilogy(GWO_cg_curve,'Color','r')
title('Objective space')
xlabel('Iteration');
ylabel('Best score obtained so far');
axis tight
grid on
box on
display(['The best solution obtained by GWO is : ', num2str(Best_pos)]);
display(['The best optimal value of the objective funciton found by GWO is : ', num2str(Best_score)]);

8.运行结果

9.参考文献

[1]Seyedali Mirjalili,Seyed Mohammad Mirjalili,Andrew Lewis. Grey Wolf Optimizer[J]. Advances in Engineering Software,2014,69.