【Matlab】智能优化算法_蝗虫算法GOA
- 1.背景介绍
- 2.基本思想
- 3.公式推导
- 4.算法流程图
- 5.文件结构
- 6.伪代码
- 7.详细代码及注释
- 7.1 distance.m
- 7.2 func_plot.m
- 7.3 Get_Functions_details.m
- 7.4 GOA.m
- 7.5 initialization.m
- 7.6 main.m
- 7.7 S_func.m
- 8.运行结果
- 9.参考文献
1.背景介绍
蝗虫是昆虫。由于它们对作物生产和农业造成破坏，因此被认为是一种害虫。蝗虫的生命周期如图1所示。尽管蝗虫通常在自然界中单独出现，但它们是所有生物中最大的群体之一。蜂群的规模可能是持续的，对农民来说是一场噩梦。蝗虫群的独特之处在于，其群集行为在若虫和成年期都有。数以百万计的若虫蝗虫像滚动的圆柱体一样跳跃和移动。在它们的路径上，它们几乎吃掉了所有的植被。在这种行为之后，当它们成年后，就会在空中形成一个群体。这就是蝗虫长途迁徙的方式。
2.基本思想
蝗虫群在幼虫期的主要特征是运动缓慢，步伐很小。相比之下，长距离和突然的运动是成年蜂群的基本特征。寻找食物来源是蝗虫群集的另一个重要特征。正如引言中所讨论的，受自然启发的算法在逻辑上将搜索过程分为两种趋势：探索和开发。在解释中，搜索代理被鼓励突然移动，而他们在利用过程中往往会在本地移动。这两项功能以及目标搜索都是由蝗虫natu集会执行的。因此，如果我们找到一种对这种行为进行数学建模的方法，我们就可以设计一种新的受自然启发的算法。
3.公式推导
用于模拟蝗虫群集行为的数学模型如下所示：
其中，Xi定义了第i只蝗虫的位置，Si是社会相互作用，Gi是第i只蚱蜢的重力，Ai表示风平流。注意，为了提供随机行为，方程式可以写成Xi=r1Si+r2Gi+r3Ai，其中r1、r2和r3是[0,1]中的随机数。
其中dij是第i和第j蝗虫之间的距离，计算为dij=| xj−xi |，s是定义社会力量强度的函数，例如公式（2.3），而dij=(x j−xi )/dij是从第i蝗虫到第j蝗的单位向量。
定义社会力量的s函数计算如下：
其中f表示吸引力的强度，l是吸引力的长度尺度。
函数s如图2，显示了它如何影响蝗虫的社会互动（吸引和排斥）。
从图中可以看出，考虑了从0到15的距离。排斥发生在区间[0 2.079]。当一只蝗虫与另一只蝗虫相距2.079个单位时，既没有吸引力，也没有排斥力。这被称为舒适区或舒适距离。图2还显示，吸引力从2.079距离单位增加到接近4，然后逐渐减小。改变方程（2.3）中的参数l和f导致人工蝗虫的不同社会行为。为了了解这两个参数的影响，在图中3重新绘制了函数s，独立地改变l和f。
该图显示，参数l和f显著改变了舒适区、吸引区和排斥区。应该注意的是，对于某些值（例如l＝1.0或f＝1.0），吸引或再脉动区域非常小。从所有这些值中，我们选择了l=1.5和f=0.5。
使用函数s的蝗虫和舒适区之间相互作用的概念模型如图4所示。可以注意到，在简化的形式中，这种社会互动是一些早期蝗虫群集模型的动力。
尽管函数s能够将两只蝗虫之间的空间划分为排斥区、舒适区和吸引区，但如图2和图3所示，当距离大于10时，该函数会返回接近零的值。因此，该功能无法在蝗虫之间施加较大距离的强大力。为了解决这个问题，我们绘制了蝗虫在[1,4]间隔内的距离图。该区间中函数s的形状如图2（右）所示。
公式（2.1）中的G分量计算如下：
其中g是引力常数，eg表示朝向地球中心的单位向量。
公式（2.1）中的A分量计算如下：
其中u是一个常数漂移，ew是风向上的单位向量。
Nymph蝗虫没有翅膀，所以它们的运动与风向高度相关。
将方程（2.1）中的S、G和A代入，该方程可以导出如下：
由于若虫蝗虫降落在地面上，它们的位置不应该低于阈值。然而，我们不会在群模拟和优化算法中使用这个方程，因为它会阻止算法探索和利用解决方案周围的搜索空间。事实上，用于蜂群的模型是在自由空间中。因此，使用了方程（2.6），可以模拟蝗虫群中蝗虫之间的相互作用。图5和图6显示了使用该方程的两个蜂群在2D和3D空间中的行为。在这两张图中，20只人工蝗虫需要移动10个单位的时间。
图5显示了方程（2.6）如何使初始随机种群更接近，直到它们形成一个统一的、受调节的种群。经过10个单位的时间后，所有的蝗虫都到达了舒适区，不再移动。在图6中的3D空间中观察到了相同的行为。这表明该数学模型能够在二维、三维和超维空间中模拟蝗虫群。
然而，这个数学模型不能直接用于解决优化问题，主要是因为蝗虫很快到达舒适区，而且蝗虫群不会收敛到指定的点。该方程的修改版本如下所示，以解决优化问题。
等式（2.7）表明，蝗虫的下一个位置是基于其当前位置、目标的位置和所有其他蝗虫的位置来定义的。注意，这个方程的第一个分量考虑了当前蝗虫相对于其他蝗虫的位置。事实上，我们已经考虑了所有蝗虫的状态，以确定搜索代理在目标周围的位置。
这里还值得一提的是，自适应参数c已经在等式（2.7）中使用了两次，原因如下：
应该注意的是，内部c有助于减少蝗虫之间的排斥力/吸引力，与迭代次数成比例，而外部c随着迭代次数的增加减少了目标周围的搜索范围。
总之，方程的第一项（2.7），即总和，考虑了其他蝗虫的位置，并实现了蝗虫在自然界中的相互作用。第二个术语T d模拟了它们向食物来源移动的十种倾向。此外，参数c模拟蝗虫接近食物来源并最终消耗食物的减速。为了提供更随机的行为，作为一种替代方案，等式（2.7）中的两个项都可以用随机值多重叠加。此外，可以将单个项与随机值相乘，以提供蝗虫相互作用或食物来源倾向的随机行为。
4.算法流程图
5.文件结构
distance.m % 距离
func_plot.m % 绘制的基准函数
Get_Functions_details.m % 基准的全部信息和实现
GOA.m % 蝗虫优化算法
initialization.m % 初始化
main.m % 主函数
S_func.m % S函数6.伪代码
7.详细代码及注释
7.1 distance.m
function d = distance(a,b)
d=sqrt((a(1)-b(1))^2+(a(2)-b(2))^2);7.2 func_plot.m
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');
end7.3 Get_Functions_details.m
function [lb,ub,dim,fobj] = Get_Functions_details(F)
switch F
case 'F1'
fobj = @F1;
lb=-100;
ub=100;
dim=5;
case 'F2'
fobj = @F2;
lb=-10;
ub=10;
dim=5;
case 'F3'
fobj = @F3;
lb=-100;
ub=100;
dim=5;
case 'F4'
fobj = @F4;
lb=-100;
ub=100;
dim=5;
case 'F5'
fobj = @F5;
lb=-30;
ub=30;
dim=5;
case 'F6'
fobj = @F6;
lb=-100;
ub=100;
dim=5;
case 'F7'
fobj = @F7;
lb=-1.28;
ub=1.28;
dim=5;
case 'F8'
fobj = @F8;
lb=-500;
ub=500;
dim=5;
case 'F9'
fobj = @F9;
lb=-5.12;
ub=5.12;
dim=5;
case 'F10'
fobj = @F10;
lb=-32;
ub=32;
dim=5;
case 'F11'
fobj = @F11;
lb=-600;
ub=600;
dim=5;
case 'F12'
fobj = @F12;
lb=-50;
ub=50;
dim=5;
case 'F13'
fobj = @F13;
lb=-50;
ub=50;
dim=5;
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));
end7.4 GOA.m
function [TargetFitness,TargetPosition,Convergence_curve,Trajectories,fitness_history, position_history]=GOA(N, Max_iter, lb,ub, dim, fobj)
disp('GOA is now estimating the global optimum for your problem....')
flag=0;
if size(ub,1)==1
ub=ones(dim,1)*ub;
lb=ones(dim,1)*lb;
end
if (rem(dim,2)~=0) % this algorithm should be run with a even number of variables. This line is to handle odd number of variables
dim = dim+1;
ub = [ub; 100];
lb = [lb; -100];
flag=1;
end
%Initialize the population of grasshoppers
GrassHopperPositions=initialization(N,dim,ub,lb);
GrassHopperFitness = zeros(1,N);
fitness_history=zeros(N,Max_iter);
position_history=zeros(N,Max_iter,dim);
Convergence_curve=zeros(1,Max_iter);
Trajectories=zeros(N,Max_iter);
cMax=1;
cMin=0.00004;
%Calculate the fitness of initial grasshoppers
for i=1:size(GrassHopperPositions,1)
if flag == 1
GrassHopperFitness(1,i)=fobj(GrassHopperPositions(i,1:end-1));
else
GrassHopperFitness(1,i)=fobj(GrassHopperPositions(i,:));
end
fitness_history(i,1)=GrassHopperFitness(1,i);
position_history(i,1,:)=GrassHopperPositions(i,:);
Trajectories(:,1)=GrassHopperPositions(:,1);
end
[sorted_fitness,sorted_indexes]=sort(GrassHopperFitness);
% Find the best grasshopper (target) in the first population
for newindex=1:N
Sorted_grasshopper(newindex,:)=GrassHopperPositions(sorted_indexes(newindex),:);
end
TargetPosition=Sorted_grasshopper(1,:);
TargetFitness=sorted_fitness(1);
% Main loop
l=2; % Start from the second iteration since the first iteration was dedicated to calculating the fitness of antlions
while l<Max_iter+1
c=cMax-l*((cMax-cMin)/Max_iter); % Eq. (2.8) in the paper
for i=1:size(GrassHopperPositions,1)
temp= GrassHopperPositions';
for k=1:2:dim
S_i=zeros(2,1);
for j=1:N
if i~=j
Dist=distance(temp(k:k+1,j), temp(k:k+1,i)); % Calculate the distance between two grasshoppers
r_ij_vec=(temp(k:k+1,j)-temp(k:k+1,i))/(Dist+eps); % xj-xi/dij in Eq. (2.7)
xj_xi=2+rem(Dist,2); % |xjd - xid| in Eq. (2.7)
s_ij=((ub(k:k+1) - lb(k:k+1))*c/2)*S_func(xj_xi).*r_ij_vec; % The first part inside the big bracket in Eq. (2.7)
S_i=S_i+s_ij;
end
end
S_i_total(k:k+1, :) = S_i;
end
X_new = c * S_i_total'+ (TargetPosition); % Eq. (2.7) in the paper
GrassHopperPositions_temp(i,:)=X_new';
end
% GrassHopperPositions
GrassHopperPositions=GrassHopperPositions_temp;
for i=1:size(GrassHopperPositions,1)
% Relocate grasshoppers that go outside the search space
Tp=GrassHopperPositions(i,:)>ub';Tm=GrassHopperPositions(i,:)<lb';GrassHopperPositions(i,:)=(GrassHopperPositions(i,:).*(~(Tp+Tm)))+ub'.*Tp+lb'.*Tm;
% Calculating the objective values for all grasshoppers
if flag == 1
GrassHopperFitness(1,i)=fobj(GrassHopperPositions(i,1:end-1));
else
GrassHopperFitness(1,i)=fobj(GrassHopperPositions(i,:));
end
fitness_history(i,l)=GrassHopperFitness(1,i);
position_history(i,l,:)=GrassHopperPositions(i,:);
Trajectories(:,l)=GrassHopperPositions(:,1);
% Update the target
if GrassHopperFitness(1,i)<TargetFitness
TargetPosition=GrassHopperPositions(i,:);
TargetFitness=GrassHopperFitness(1,i);
end
end
Convergence_curve(l)=TargetFitness;
disp(['In iteration #', num2str(l), ' , target''s objective = ', num2str(TargetFitness)])
l = l + 1;
end
if (flag==1)
TargetPosition = TargetPosition(1:dim-1);
end7.5 initialization.m
function [X]=initialization(N,dim,up,down)
if size(up,1)==1
X=rand(N,dim).*(up-down)+down;
end
if size(up,1)>1
for i=1:dim
high=up(i);low=down(i);
X(:,i)=rand(1,N).*(high-low)+low;
end
end7.6 main.m
clear all
clc
SearchAgents_no=100; % Number of search agents
Function_name='F1'; % 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);
[Target_score,Target_pos,GOA_cg_curve, Trajectories,fitness_history, position_history]=GOA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj);
figure('Position',[454 445 894 297])
%Draw search space
subplot(1,5,1);
func_plot(Function_name);
title('Parameter space')
xlabel('x_1');
ylabel('x_2');
zlabel([Function_name,'( x_1 , x_2 )'])
box on
axis tight
subplot(1,5,2);
hold on
for k1 = 1: size(position_history,1)
for k2 = 1: size(position_history,2)
plot(position_history(k1,k2,1),position_history(k1,k2,2),'.','markersize',1,'MarkerEdgeColor','k','markerfacecolor','k');
end
end
plot(Target_pos(1),Target_pos(2),'.','markersize',10,'MarkerEdgeColor','r','markerfacecolor','r');
title('Search history (x1 and x2 only)')
xlabel('x1')
ylabel('x2')
box on
axis tight
subplot(1,5,3);
hold on
plot(Trajectories(1,:));
title('Trajectory of 1st grasshopper')
xlabel('Iteration#')
box on
axis tight
subplot(1,5,4);
hold on
plot(mean(fitness_history));
title('Average fitness of all grasshoppers')
xlabel('Iteration#')
box on
axis tight
%Draw objective space
subplot(1,5,5);
semilogy(GOA_cg_curve,'Color','r')
title('Convergence curve')
xlabel('Iteration#');
ylabel('Best score obtained so far');
box on
axis tight
set(gcf, 'position' , [39 479 1727 267]);
display(['The best solution obtained by GOA is : ', num2str(Target_pos)]);
display(['The best optimal value of the objective funciton found by GOA is : ', num2str(Target_score)]);7.7 S_func.m
function o=S_func(r)
f=0.5;
l=1.5;
o=f*exp(-r/l)-exp(-r); % Eq. (2.3) in the paper
end8.运行结果
