【Matlab】智能优化算法_鲸鱼优化算法WOA
【Matlab】智能优化算法_鲸鱼优化算法WOA
- 1.背景介绍
- 2.基本思想
- 3.公式推导
- 3.1 包围猎物
- 3.2 发泡网攻击
- 3.2.1 收缩包围机制
- 3.2.2 螺旋更新位置
- 3.3 搜索捕食
- 4.算法流程图
- 5.文件结构
- 6.伪代码
- 7.详细代码及注释
- 7.1 func_plot.m
- 7.2 Get_Functions_details.m
- 7.3 initialization.m
- 7.4 main.m
- 7.5 WOA.m
- 8.运行结果
- 9.参考文献
1.背景介绍
鲸鱼优化算法提出的背景是为了模拟座头鲸的捕猎行为和螺旋线运动，从而设计一种新的群体智能优化算法。Mirjalili 等人通过观察座头鲸的社会行为和迁徙模式，发现它们具有以下三种特征：
- 座头鲸可以通过声波或其他方式感知猎物的位置，并向其靠近。
- 座头鲸可以通过环绕猎物产生气泡网来围捕猎物，并逐渐缩小气泡网的半径。
- 座头鲸可以通过螺旋线运动来追逐猎物或探索新的领域。
基于这些特征，Mirjalili 等人提出了一种数学模型来描述座头鲸的运动，并将其应用于优化问题。
2.基本思想
- 鲸鱼优化算法是模仿自然界中鲸鱼捕食行为的新型群体智能优化算法，而鲸鱼的捕食行为主要分为3类：1）包围猎物；2）发泡网攻击；3）搜索捕食。因此，在使用WOA求解问题之前，需要对上述3类捕食行为进行建模，即用数学公式表达上述3类捕食行为。
- 鲸鱼捕食行为的目的是捕获猎物，一群鲸鱼在共同寻找猎物时，一定会存在某条鲸鱼先发现猎物的情况，这时候其它鲸鱼一定会向这条发现猎物的鲸鱼游来争抢猎物。
- 将上述捕食过程可以应用到WOA求解问题的过程中，即一个解就可以用一个鲸鱼个体表示，若干个解就可以用若干个鲸鱼个体表示。在使用WOA搜索问题解的过程就可以看作是若干个鲸鱼个体不断更新个体位置，直至搜索到满意的解为止。
3.公式推导
在本节中，首先提供了包围猎物、螺旋泡网喂食策略和搜索猎物的数学模型。然后提出WOA算法
3.1 包围猎物
座头鲸能够识别猎物的位置并环绕它们。由于最优设计在搜索空间中的位置先验未知，WOA算法假设当前最佳候选解是目标猎物或接近最优。在定义了最佳搜索代理之后，其他搜索代理将因此尝试向最佳搜索代理更新其位置。该行为由以下等式表示：
其中t表示当前迭代，A和C是系数向量，X * 是迄今为止获得的最佳解的位置向量，X是位置向量，||是绝对值，·是逐元素乘法。这里值得一提的是，如果有更好的解决方案，应在每次迭代中更新X * 。矢量A和C的计算如下：
其中a在其操作过程中（在勘探和开采阶段）从2线性减少到0，r是[0,1]中的随机向量。
例如。3（a）说明了方程式（2.2）对于2D问题的原理。搜索代理的位置（X，Y）可以根据当前最佳记录（X*，Y*）的位置进行更新。通过调整A和C向量的值，可以相对于当前位置实现最佳代理周围的不同位置。搜索代理在3D空间中的可能更新位置也如图3（b）所示。应该注意，通过定义随机向量（r），可以到达位于图3所示关键点之间的搜索空间中的任何位置。因此，等式（2.2）允许任何搜索代理更新其在当前最佳解附近的位置，并模拟包围猎物。
3.2 发泡网攻击
为了对座头鲸的气泡网行为进行数学建模，设计了以下两种方法：
3.2.1 收缩包围机制
通过减小方程（2.3）中a的值，可以实现这种行为。注意，a的波动范围也会减小。A在区间[−a，a]中的随机值，其中a在迭代过程中从2减少到0。在[−1,1]中为A设置随机值，可以在代理的原始位置和当前最佳代理的位置之间的任何位置定义搜索代理的新位置。图4（a）显示了二维空间中0≤a≤1可实现的从（X，Y）到（X*，Y*）的可能位置。
3.2.2 螺旋更新位置
如图4（b）所示，该方法首先计算位于（X，Y）的鲸鱼和位于（X*，Y*）的猎物之间的距离。然后在鲸鱼和猎物的位置之间建立螺旋方程，以模拟座头鲸的螺旋形运动，如下所示：
其中D=|X*（t）−X（t）|，表示第i头鲸鱼与猎物的距离（目前获得的最佳解），b是定义对数螺旋形状的常数，l是[−1,1]中的随机数和·是逐元素乘法。
请注意，座头鲸在一个缩小的圆圈内和沿着螺旋形的路径同时围绕猎物游动。为了模拟这种同时发生的行为，我们假设有50%的可能性在收缩环绕机制或螺旋模型之间进行选择，以在优化过程中更新鲸鱼的位置。数学模型如下：
其中p是[0,1]中的随机数。
3.3 搜索捕食
其中Xrand是从当前种群中选择的随机位置向量（随机鲸鱼）。图5描述了a>1的特定溶液周围的一些可能位置。
WOA算法从一组随机解开始。在每一次迭代中，搜索代理更新其相对于随机选择的搜索代理或迄今为止获得的最佳解决方案的位置。a参数从2减小到0，以便分别提供勘探和开采。当|A|＞1时选择随机搜索代理，而当|A|<1时选择最佳解决方案以更新搜索代理的位置。根据p的值，WOA可以在螺旋或圆形运动之间切换。最后，WOA算法通过满足终止标准而终止.
4.算法流程图
5.文件结构
func_plot.m % 绘图函数
Get_Functions_details.m % 获取函数信息
initialization.m % 初始化
main.m % 主函数
WOA.m % 鲸鱼优化算法
6.伪代码
7.详细代码及注释
7.1 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:Lfor j=1:Lif 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')==0f(i,j)=fobj([x(i),y(j)]);endif strcmp(func_name,'F15')==1f(i,j)=fobj([x(i),y(j),0,0]);endif strcmp(func_name,'F19')==1f(i,j)=fobj([x(i),y(j),0]);endif strcmp(func_name,'F20')==1f(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')==1f(i,j)=fobj([x(i),y(j),0,0]);end end
endsurfc(x,y,f,'LineStyle','none');end7.2 Get_Functions_details.m
function [lb,ub,dim,fobj] = Get_Functions_details(F)
switch Fcase '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;
endend% F1function o = F1(x)
o=sum(x.^2);
end% F2function o = F2(x)
o=sum(abs(x))+prod(abs(x));
end% F3function o = F3(x)
dim=size(x,2);
o=0;
for i=1:dimo=o+sum(x(1:i))^2;
end
end% F4function o = F4(x)
o=max(abs(x));
end% F5function 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% F6function o = F6(x)
o=sum(abs((x+.5)).^2);
end% F7function o = F7(x)
dim=size(x,2);
o=sum([1:dim].*(x.^4))+rand;
end% F8function o = F8(x)
o=sum(-x.*sin(sqrt(abs(x))));
end% F9function o = F9(x)
dim=size(x,2);
o=sum(x.^2-10*cos(2*pi.*x))+10*dim;
end% F10function 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% F11function o = F11(x)
dim=size(x,2);
o=sum(x.^2)/4000-prod(cos(x./sqrt([1:dim])))+1;
end% F12function 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% F13function 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% F14function 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:25bS(j)=sum((x'-aS(:,j)).^6);
end
o=(1/500+sum(1./([1:25]+bS))).^(-1);
end% F15function 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% F16function 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% F17function 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% F18function 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% F19function 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:4o=o-cH(i)*exp(-(sum(aH(i,:).*((x-pH(i,:)).^2))));
end
end% F20function 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:4o=o-cH(i)*exp(-(sum(aH(i,:).*((x-pH(i,:)).^2))));
end
end% F21function 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:5o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end% F22function 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:7o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end% F23function 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:10o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
endfunction o=Ufun(x,a,k,m)
o=k.*((x-a).^m).*(x>a)+k.*((-x-a).^m).*(x<(-a));
end
7.3 initialization.m
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==1Positions=rand(SearchAgents_no,dim).*(ub-lb)+lb;
end% If each variable has a different lb and ub
if Boundary_no>1for i=1:dimub_i=ub(i);lb_i=lb(i);Positions(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;end
end
7.4 main.m
clc
clear
close all
SearchAgents_no=30; % 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);
[Best_score,Best_pos,WOA_cg_curve]=WOA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj);
figure('Position',[269 240 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(WOA_cg_curve,'Color','r')
title('Objective space')
xlabel('Iteration');
ylabel('Best score obtained so far');axis tight
grid on
box on
legend('WOA')display(['The best solution obtained by WOA is : ', num2str(Best_pos)]);
display(['The best optimal value of the objective funciton found by WOA is : ', num2str(Best_score)]);
7.5 WOA.m
function [Leader_score,Leader_pos,Convergence_curve]=WOA(SearchAgents_no,Max_iter,lb,ub,dim,fobj)% initialize position vector and score for the leader
Leader_pos=zeros(1,dim);
Leader_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);t=0;% Loop counter% Main loop
while t<Max_iterfor i=1:size(Positions,1)% Return back the search agents that go beyond the boundaries of the search spaceFlag4ub=Positions(i,:)>ub;Flag4lb=Positions(i,:)<lb;Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;% Calculate objective function for each search agentfitness=fobj(Positions(i,:));% Update the leaderif fitness<Leader_score % Change this to > for maximization problemLeader_score=fitness; % Update alphaLeader_pos=Positions(i,:);endenda=2-t*((2)/Max_iter); % a decreases linearly fron 2 to 0 in Eq. (2.3)% a2 linearly dicreases from -1 to -2 to calculate t in Eq. (3.12)a2=-1+t*((-1)/Max_iter);% Update the Position of search agents for i=1:size(Positions,1)r1=rand(); % r1 is a random number in [0,1]r2=rand(); % r2 is a random number in [0,1]A=2*a*r1-a; % Eq. (2.3) in the paperC=2*r2; % Eq. (2.4) in the paperb=1; % parameters in Eq. (2.5)l=(a2-1)*rand+1; % parameters in Eq. (2.5)p = rand(); % p in Eq. (2.6)for j=1:size(Positions,2)if p<0.5 if abs(A)>=1rand_leader_index = floor(SearchAgents_no*rand()+1);X_rand = Positions(rand_leader_index, :);D_X_rand=abs(C*X_rand(j)-Positions(i,j)); % Eq. (2.7)Positions(i,j)=X_rand(j)-A*D_X_rand; % Eq. (2.8)elseif abs(A)<1D_Leader=abs(C*Leader_pos(j)-Positions(i,j)); % Eq. (2.1)Positions(i,j)=Leader_pos(j)-A*D_Leader; % Eq. (2.2)endelseif p>=0.5distance2Leader=abs(Leader_pos(j)-Positions(i,j));% Eq. (2.5)Positions(i,j)=distance2Leader*exp(b.*l).*cos(l.*2*pi)+Leader_pos(j);endendendt=t+1;Convergence_curve(t)=Leader_score;[t Leader_score]
end8.运行结果
9.参考文献
[1]Seyedali Mirjalili,Andrew Lewis. The Whale Optimization Algorithm[J]. Advances in Engineering Software,2016,95.
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！