鲸鱼优化算法WOA求解旅行商TSP优化问题(2022.6.2)
引言
????????作为世界上最大的哺乳动物,一条成年鲸鱼的身长可达到30m,体重能达到180t。鲸鱼大多被认为是食肉动物,它们从不睡觉,因为必须从海洋表面呼吸。事实证明,鲸鱼可以像人类一样思考、学习、判断、交流 ,它们也被认为是具有情感的高智商动物。另一件有趣的事情是鲸鱼的社交行为:它们独居或群居,但大多是处于群体当中。成年座头鲸(Megaptera novaeangliae) 几乎和校车一样大。它们最喜欢的猎物是磷虾和小鱼群。
1、鲸鱼优化算法WOA
WOA算法原始论文
Web of Science高被引论文
|
百度学术检索结果
|
????????2016年,来自澳大利亚的两位研究人员Mirjalili Seyedali 和Lewis Andrew 通过模拟自然界生物座头鲸的社会行为,提出了一种元启发式的优化算法——鲸鱼优化算法WOA,撰写的论文在Advances in Engineering Software 期刊上在线发表,被SCI 收录,引起了学者们的广泛关注,迄今为止论文的被引频次高达3496。
座头鲸的气泡网法捕食行为示意图
????????算法的核心在于座头鲸的气泡网式捕捉策略,算法利用三个操作算子来模拟座头鲸对猎物的搜索、包围猎物和气泡网觅食行为 。经过29个数学优化问题和6个结构设计问题的实际测试,结果证明,与其他传统的元启发式算法相比,WOA 具有十分明显的计算优势。(注:WOA算法源代码已公开,可点击下载)
1.1 WOA算法原理介绍
1.1.1 包围猎物
????????座头鲸能够识别猎物位置并且对其进行转圈包围,WOA 算法假设当前最优的候选位置就是目标猎物或者说更接近最优猎物的位置,在定义了最好的搜索智能体后,其他的搜索智能体就会向着最好的搜索智能体来尝试更新它们的位置,这个行为可以用下列方程来表示:
D
?
=
∣
C
?
?
X
?
?
(
t
)
?
X
?
(
t
)
∣
(公式1)
\vec{D}=|\vec{C}\cdot \vec{X^{*}}(t)-\vec{X}(t)| \tag{公式1}
D
=∣C
?X?
(t)?X
(t)∣(公式1)
X
?
(
t
+
1
)
=
X
?
?
(
t
)
?
A
?
?
D
?
(公式2)
\vec{X}(t+1)=\vec{X^{*}}(t)-\vec{A}\cdot \vec{D} \tag{公式2}
X
(t+1)=X?
(t)?A
?D
(公式2) ????????其中,
t
t
t 表示当前代数;
A
?
\vec{A}
A
和
C
?
\vec{C}
C
是系数向量;
X
?
X^{*}
X? 是目前为止所获取最优解的位置向量;
X
?
\vec{X}
X
是位置向量;
∣
∣
| |
∣∣ 是绝对值;
?
\cdot
? 是点积,对应元素逐点相乘。
A
?
=
2
a
?
?
r
?
?
a
?
(公式3)
\vec{A}=2\vec{a}\cdot \vec{r} -\vec{a} \tag{公式3}
A
=2a
?r
?a
(公式3)
C
?
=
2
?
r
?
(公式4)
\vec{C}=2\cdot{\vec{r}} \tag{公式4}
C
=2?r
(公式4) ????????
a
?
\vec{a}
a
随着代数的增加线性递减,取值范围为[2,1];
r
?
\vec{r}
r
是一个取值范围在[0,1]的随机向量。
????????如上图所示,无论是二维空间还是三维空间,通过对当前位置调整
A
?
\vec{A}
A
和
C
?
\vec{C}
C
的值可到达最好个体周围的不同位置,因此,公式2可用来根据当前最优解的邻域对一些搜索个体的位置进行更新,这样来实现对猎物的包围。
1.1.2 气泡网式攻击猎物(开发阶段)
????????为了对座头鲸的气泡网攻击猎物方法进行建模,故而设计了以下两种方法,如下图所示。
- 收缩环绕机制:这个行为可通过减小公式3中
a
?
\vec{a}
a
的值来实现,当a随着代数增加从2逐渐下降到0,
A
?
\vec{A}
A
的取值范围,为[-a,a]。如果设定
A
?
\vec{A}
A
的随机值在[-1,1],那么一个搜索个体的新位置就会位于其原来位置和当前最优个体位置之间。上图中(a)展示了2维空间
0
≤
A
≤
1
0≤A≤1
0≤A≤1 时从
(
X
,
Y
)
(X,Y)
(X,Y) 向
(
X
?
,
Y
?
)
(X^{*},Y^{*})
(X?,Y?) 的可能搜索位置。
- 螺旋线更新位置:如上图中(b)所示,首先计算鲸鱼所在位置
(
X
,
Y
)
(X,Y)
(X,Y)到猎物所在位置
(
X
?
,
Y
?
)
(X^{*},Y^{*})
(X?,Y?)的距离,然后在鲸鱼和猎物的位置之间建立一个螺旋方程来模拟座头鲸的螺旋形运动,如下所示:
X
?
(
t
+
1
)
=
D
′
?
?
e
b
l
?
c
o
s
(
2
π
l
)
+
X
?
?
(
t
)
(公式5)
\vec{X}(t+1)=\vec{D^{'}} \cdot e^{bl} \cdot cos(2\pi l) + \vec{X^{*}}(t) \tag{公式5}
X
(t+1)=D′
?ebl?cos(2πl)+X?
(t)(公式5)????????式中,
D
′
?
=
X
?
?
(
t
)
?
X
?
(
t
)
\vec{D^{'}}=\vec{X^{*}}(t) -\vec{X}(t)
D′
=X?
(t)?X
(t) 表示第
i
i
i 代鲸鱼到目前最优猎物的距离;
b
b
b 是一个定义对数螺旋线形状的常数;
l
l
l 是一个取值范围在[-1,1]区间的随机数。
????????注意到座头鲸在围绕猎物游动时会同时进行缩小包围圈和螺旋形位置更新这两种行为,为了建模这个同时进行的行为,假设在优化期间存在50%的概率 在缩小包围圈和螺旋形更新位置这两个行为中做一个随机选择,数学模型如下:
X
?
(
t
+
1
)
=
{
X
?
?
(
t
)
?
A
?
?
D
?
p
<
0.5
D
′
?
?
e
b
l
?
c
o
s
(
2
π
l
)
+
X
?
?
(
t
)
p
≥
0.5
(公式6)
\vec{X}(t+1)=\left\{\begin{matrix} \vec{X^{*}}(t)-\vec{A}\cdot \vec{D} & p < 0.5 \\ \vec{D^{'}} \cdot e^{bl} \cdot cos(2\pi l) + \vec{X^{*}}(t) &p \geq 0.5 \\ \end{matrix}\right. \tag{公式6}
X
(t+1)={X?
(t)?A
?D
D′
?ebl?cos(2πl)+X?
(t)?p<0.5p≥0.5?(公式6)????????其中,随机概率
p
∈
[
0
,
1
]
p∈[0,1]
p∈[0,1]。
1.1.3 寻找猎物(探索阶段)
????????座头鲸在寻找猎物时是按照随机的方式根据彼此的位置进行搜索,使用具有小于-1或者大于1的随机值
A
?
\vec{A}
A
来迫使搜索个体向远离参考鲸鱼个体的方向移动,根据随机选择的个体进行距离更新而不是选择最优个体,
A
?
>
1
\vec{A}>1
A
>1强调在探索阶段WOA算法可进行全局搜索,具体公式如下:
D
?
=
∣
C
?
?
X
r
a
n
d
?
(
t
)
?
X
?
(
t
)
∣
(公式7)
\vec{D}=|\vec{C}\cdot \vec{X_{rand}}(t)-\vec{X}(t)| \tag{公式7}
D
=∣C
?Xrand?
?(t)?X
(t)∣(公式7)
X
?
(
t
+
1
)
=
X
r
a
n
d
?
(
t
)
?
A
?
?
D
?
(公式8)
\vec{X}(t+1)=\vec{X_{rand}}(t)-\vec{A}\cdot \vec{D} \tag{公式8}
X
(t+1)=Xrand?
?(t)?A
?D
(公式8)
????????式中,
X
r
a
n
d
(
t
)
X_{rand}(t)
Xrand?(t)表示从当前代(第t代)随机选择的座鱼鲸个体的位置向量。
1.2 WOA算法流程
????????WOA算法的伪代码 和计算流程 如下图所示。
WOA算法的伪代码
|
WOA算法计算流程
|
1.3 WOA官方源码
WOA源代码主页
????????在WOA 官网中给出了源代码Source Codes of WOA和工具箱Source Codes of WOA toolbox的下载地址,下载至本地并解压,如下图所示。
源码文件夹内容
|
工具箱文件夹内容
|
2、WOA源码测试运行(Matlab)
????????测试运行环境为Matlab R2021a 64位
2.1 WOA代码简介
????????WOA 代码共包含5个.m 脚本文件:Get_Functions_details.m、func_plot.m、initialization.m、WOA.m和main.m,其中前两个包含23个函数的定义,initialization.m 是种群初始化文件,WOA.m 中执行WOA 算法整个流程,main.m 本质是调用WOA.m 执行,可进行传参,设定目标函数(F1~F23 )后即可运行。23个经典数学函数的表达式、自变量维度、变量范围、最小值的详细信息如下所示:
2.1.1 Get_Functions_details.m
%_________________________________________________________________________%
% Whale Optimization Algorithm (WOA) source codes demo 1.0 %
% %
% Developed in MATLAB R2011b(7.13) %
% %
% Author and programmer: Seyedali Mirjalili %
% %
% e-Mail: ali.mirjalili@gmail.com %
% seyedali.mirjalili@griffithuni.edu.au %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: S. Mirjalili, A. Lewis %
% The Whale Optimization Algorithm, %
% Advances in Engineering Software , in press, %
% DOI: http://dx.doi.org/10.1016/j.advengsoft.2016.01.008 %
% %
%_________________________________________________________________________%
% 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
2.1.2 func_plot.m
%_________________________________________________________________________%
% Whale Optimization Algorithm (WOA) source codes demo 1.0 %
% %
% Developed in MATLAB R2011b(7.13) %
% %
% Author and programmer: Seyedali Mirjalili %
% %
% e-Mail: ali.mirjalili@gmail.com %
% seyedali.mirjalili@griffithuni.edu.au %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: S. Mirjalili, A. Lewis %
% The Whale Optimization Algorithm, %
% Advances in Engineering Software , in press, %
% DOI: http://dx.doi.org/10.1016/j.advengsoft.2016.01.008 %
% %
%_________________________________________________________________________%
% 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
2.1.3 initialization.m
%_________________________________________________________________________%
% Whale Optimization Algorithm (WOA) source codes demo 1.0 %
% %
% Developed in MATLAB R2011b(7.13) %
% %
% Author and programmer: Seyedali Mirjalili %
% %
% e-Mail: ali.mirjalili@gmail.com %
% seyedali.mirjalili@griffithuni.edu.au %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: S. Mirjalili, A. Lewis %
% The Whale Optimization Algorithm, %
% Advances in Engineering Software , in press, %
% DOI: http://dx.doi.org/10.1016/j.advengsoft.2016.01.008 %
% %
%_________________________________________________________________________%
% 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
2.1.4 WOA.m
%_________________________________________________________________________%
% Whale Optimization Algorithm (WOA) source codes demo 1.0 %
% %
% Developed in MATLAB R2011b(7.13) %
% %
% Author and programmer: Seyedali Mirjalili %
% %
% e-Mail: ali.mirjalili@gmail.com %
% seyedali.mirjalili@griffithuni.edu.au %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: S. Mirjalili, A. Lewis %
% The Whale Optimization Algorithm, %
% Advances in Engineering Software , in press, %
% DOI: http://dx.doi.org/10.1016/j.advengsoft.2016.01.008 %
% %
%_________________________________________________________________________%
% The Whale Optimization Algorithm
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_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 the leader
if fitness<Leader_score % Change this to > for maximization problem
Leader_score=fitness; % Update alpha
Leader_pos=Positions(i,:);
end
end
a=2-t*((2)/Max_iter); % a decreases linearly fron 2 to 0 in Eq. (2.3)
a2=-1+t*((-1)/Max_iter);% a2 linearly dicreases from -1 to -2 to calculate t in Eq. (3.12)
% 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 paper
C=2*r2; % Eq. (2.4) in the paper
b=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)>=1
rand_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)<1
D_Leader=abs(C*Leader_pos(j)-Positions(i,j)); % Eq. (2.1)
Positions(i,j)=Leader_pos(j)-A*D_Leader; % Eq. (2.2)
end
elseif p>=0.5
distance2Leader=abs(Leader_pos(j)-Positions(i,j));
Positions(i,j)=distance2Leader*exp(b.*l).*cos(l.*2*pi)+Leader_pos(j); % Eq. (2.5)
end
end
end
t=t+1;
Convergence_curve(t)=Leader_score;
[t Leader_score]
end
2.1.5 main.m(运行文件)
%_________________________________________________________________________%
% Whale Optimization Algorithm (WOA) source codes demo 1.0 %
% %
% Developed in MATLAB R2011b(7.13) %
% %
% Author and programmer: Seyedali Mirjalili %
% %
% e-Mail: ali.mirjalili@gmail.com %
% seyedali.mirjalili@griffithuni.edu.au %
% %
% Homepage: http://www.alimirjalili.com %
% %
% Main paper: S. Mirjalili, A. Lewis %
% The Whale Optimization Algorithm, %
% Advances in Engineering Software , in press, %
% DOI: http://dx.doi.org/10.1016/j.advengsoft.2016.01.008 %
% %
%_________________________________________________________________________%
% You can simply define your cost in a seperate file and load its handle to fobj
% The initial parameters that you need are:
%__________________________________________
% fobj = @YourCostFunction
% dim = number of your variables
% Max_iteration = maximum number of generations
% SearchAgents_no = number of search agents
% lb=[lb1,lb2,...,lbn] where lbn is the lower bound of variable n
% ub=[ub1,ub2,...,ubn] where ubn is the upper bound of variable n
% If all the variables have equal lower bound you can just
% define lb and ub as two single number numbers
% To run WOA: [Best_score,Best_pos,WOA_cg_curve]=WOA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj)
%__________________________________________
clear all
clc
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)]);
2.2 运行结果(函数F1~F23)
函数F1
|
函数F2
|
函数F3
|
函数F4
|
函数F5
|
函数F6
|
函数F7
|
函数F8
|
函数F9
|
函数F10
|
函数F11
|
函数F12
|
函数F13
|
函数F14
|
函数F15
|
函数F16
|
函数F17
|
函数F18
|
函数F19
|
函数F20
|
函数F21
|
函数F22
|
函数F23
|
3、WOA求解TSP(C++ or Python)
3.1 WOA求解TSP(C++)
????????在利用优化算法求解TSP 问题时,关键在于构造问题的有效解 、确定目标函数 、求解最大值或最小值 ,在相关算子的作用下使得个体在解空间不断搜寻,进而找到相对较优的满意解。 ????????这里用到的C++ 环境为gcc 8.1.0 ,编写代码所使用的IDE为CodeBlocks 。
gcc 8.0.1
CodeBlocks 20.03
3.1.1 C++代码(main.cpp)
main.cpp
#include <iostream>
#include <fstream>
#include <string>
#include <random>
#include <ctime>
using namespace std;
const int citycount = 29;
#define PI 3.14159265358979323846
#define N 999
double round(double r)
{
return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
}
class City
{
public:
string ID;
double x, y;
void shuchu()
{
std::cout <<ID+":"<<"("<< x << "," << y <<")"<< endl;
}
};
class CityGraph
{
public:
City city[citycount];
double distance[citycount][citycount];
void Readcoordinatetxt(string txtfilename)
{
ifstream myfile(txtfilename, ios::in);
double x = 0, y = 0;
int z = 0;
if (!myfile.fail())
{
int i = 0;
while (!myfile.eof() && (myfile >> z >> x >> y))
{
city[i].ID = to_string(z);
city[i].x = x; city[i].y = y;
i++;
}
}
else
cout << "文件不存在";
myfile.close();
for (int i = 0; i < citycount; i++)
for (int j = 0; j < citycount; j++)
{
distance[i][j] = sqrt((pow((city[i].x - city[j].x), 2) + pow((city[i].y - city[j].y), 2))/10.0);
if (round(distance[i][j] < distance[i][j]))distance[i][j] = round(distance[i][j]) + 1;
else distance[i][j] = round(distance[i][j]);
}
}
void shuchu()
{
cout << "城市名称 " << "坐标x" << " " << "坐标y" << endl;
for (int i = 0; i < citycount; i++)
city[i].shuchu();
cout << "距离矩阵: " << endl;
for (int i = 0; i < citycount; i++)
{
for (int j = 0; j < citycount; j++)
{
if (j == citycount - 1)
cout << distance[i][j] << endl;
else
cout << distance[i][j] << " ";
}
}
}
};
CityGraph Map_City;
int * Random_N(int n)
{
int *geti;
geti = new int[n];
int j = 0;
while(j<n)
{
while (true)
{
int flag = -1;
int temp = rand() % n + 1;
if (j > 0)
{
int k = 0;
for(; k < j; k++)
{
if (temp == *(geti + k))break;
}
if (k == j)
{
*(geti + j) = temp;
flag = 1;
}
}
else
{
*(geti + j) = temp;
flag = 1;
}
if (flag == 1)break;
}
j++;
}
return geti;
}
class JingYu
{
public:
int *X;
double fitness;
void Init()
{
X = new int[citycount];
int *M = Random_N(citycount);
for (int j = 0; j < citycount; j++)
X[j] = *(M + j);
fitness=0;
}
void shuchu()
{
for (int j = 0; j < citycount; j++)
{
if (j == citycount -1) std::cout << X[j] << " "<<fitness<<endl;
else std::cout << X[j] << "->";
}
}
};
class WOA
{
public:
JingYu *Population;
int dimension;
double LB;
double UB;
int Pop_Size;
int Itetime;
JingYu Bestgeti;
double BestFitness;
void Init(int popsize,double dimen,double lb,double ub,int itetime,string filename)
{
Map_City.Readcoordinatetxt(filename);
Map_City.shuchu();
Pop_Size = popsize;
dimension = dimen;
LB = lb;
UB = ub;
Itetime = itetime;
Population = new JingYu[Pop_Size];
for (int i = 0; i < Pop_Size; i++)
{
Population[i].Init();
Population[i].fitness = Evaluate(Population[i]);
}
Bestgeti.Init();
for (int j = 0; j < citycount; j++)
Bestgeti.X[j] = Population[0].X[j];
Bestgeti.fitness = Evaluate(Bestgeti);
for (int i = 0; i < Pop_Size; i++)
{
if (Population[i].fitness < Bestgeti.fitness)
{
for (int j = 0; j < citycount; j++)
Bestgeti.X[j] = Population[i].X[j];
Bestgeti.fitness = Evaluate(Bestgeti);
}
}
cout<<"初始化鲸鱼种群如下:"<<endl;
for (int i = 0; i < Pop_Size; i++)
Population[i].shuchu();
cout<<"初始化最优的鲸鱼个体如下:"<<endl;
Bestgeti.shuchu();
}
void Adjuxt_validGeti(JingYu jy)
{
for(int j=0;j<dimension;j++)
{
if(jy.X[j] > LB && jy.X[j] < UB)
jy.X[j] = (int)jy.X[j];
else
jy.X[j] = int(rand() % dimension);
}
int route[citycount];
bool flag[citycount];
int biaoji[citycount];
for (int j = 0; j < citycount; j++)
{
route[j] = j + 1;
flag[j] = false;
biaoji[j] = 0;
}
for (int j = 0; j < citycount; j++)
{
int num = 0;
for (int k = 0; k < citycount; k++)
{
if (jy.X[k] == route[j])
{
biaoji[k] = 1;
num++; break;
}
}
if (num == 0) flag[j] = false;
else if (num == 1) flag[j] = true;
}
for (int k = 0; k < citycount; k++)
{
if (flag[k] == false)
{
int i = 0;
for (; i < citycount; i++)
{
if (biaoji[i] != 1)break;
}
jy.X[i] = route[k];
biaoji[i] = 1;
}
}
}
double Evaluate(JingYu jy)
{
double Sum_dist=0;
for (int j = 0; j < citycount-1; j++)
Sum_dist += Map_City.distance[jy.X[j]][jy.X[j + 1]];
Sum_dist += Map_City.distance[jy.X[citycount -1]][jy.X[0]];
jy.fitness = Sum_dist;
return Sum_dist;
}
void UpdateBestgeti()
{
for (int i = 0; i < Pop_Size; i++)
{
if (Population[i].fitness < Bestgeti.fitness)
{
for (int j = 0; j < citycount; j++)
Bestgeti.X[j] = Population[i].X[j];
}
}
Bestgeti.fitness = Evaluate(Bestgeti);
}
void ShuchuPopulation(int ite)
{
for (int i = 0; i < Pop_Size; i++)
{
cout << "第"<<ite<<"代种群中第"<<i + 1<<"个鲸鱼"<<"->";
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) cout << Population[i].X[j] <<")对应的适应度为: "<<Evaluate( Population[i])<< endl;
else if(j==0)
cout <<"("<< Population[i].X[j] << ", ";
else
cout << Population[i].X[j] << ", ";
}
}
}
void ShuchuBestgeti(int ite)
{
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) std::cout << Bestgeti.X[j] << ")" << "对应的适应度为:" << Evaluate(Bestgeti) << endl;
else if (j == 0) std::cout << "第"<<ite<<"代种群最优鲸鱼个体(" << Bestgeti.X[j] << ",";
else std::cout << Bestgeti.X[j] << ",";
}
}
void WOA_TSP(int popsize,double dimen,double lb,double ub,int Max_iter,string filename)
{
ofstream outfile;
outfile.open("result.txt",ios::trunc);
Init(popsize,dimen,lb,ub,Max_iter,filename);
outfile<<"初始化鲸鱼种群如下:"<<endl;
for (int i = 0; i < Pop_Size; i++)
{
outfile << "初始化种群中第"<<i + 1<<"个鲸鱼"<<"->";
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) outfile << Population[i].X[j] <<")对应的适应度为: "<<Evaluate( Population[i])<< endl;
else if(j==0)
outfile <<"("<< Population[i].X[j] << ", ";
else
outfile << Population[i].X[j] << ", ";
}
}
outfile<<"初始化最优的鲸鱼个体如下:"<<endl;
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) outfile << Bestgeti.X[j] << ")" << "对应的适应度为:" << Evaluate(Bestgeti) << endl;
else if (j == 0) outfile << "初始种群最优鲸鱼个体(" << Bestgeti.X[j] << ",";
else outfile << Bestgeti.X[j] << ",";
}
int t=0;
while (t < Itetime)
{
double a = 2 - t * (2 / Itetime);
double a2 = -1 + t * ((-1) / Itetime);
for (int i = 0; i < Pop_Size; i++)
{
double r1 = rand() % (N + 1) / (float)(N + 1);
double r2 = rand() % (N + 1) / (float)(N + 1);
double A = 2 * a * r1 - a;
double C = 2 * r2;
double b = 1;
double l = (a2 - 1) * (rand() % (N + 1) / (float)(N + 1)) +1;
double p = rand() % (N + 1) / (float)(N + 1);
for(int j=0;j<dimension;j++)
{
if(p < 0.5)
{
if(abs(A) >= 1)
{
int rand_leader_index = rand() % dimension;
JingYu X_rand;
X_rand.Init();
for(int k=0;k<dimension;k++)
X_rand.X[k] = Population[rand_leader_index].X[k];
double D_X_rand = abs(C * X_rand.X[j] - Population[i].X[j]);
Population[i].X[j] = int(X_rand.X[j] - A * D_X_rand);
}
else
{
double D_Leader = abs(C * Bestgeti.X[j] - Population[i].X[j]);
Population[i].X[j] = Bestgeti.X[j] - A * D_Leader;
}
}
else
{
double distance2Leader = abs(Bestgeti.X[j] - Population[i].X[j]);
Population[i].X[j] = int(distance2Leader * exp(b * l) * cos(l * 2 * PI) + Bestgeti.X[j]);
}
}
}
t = t + 1;
for (int i = 0; i < Pop_Size; i++)
{
Adjuxt_validGeti(Population[i]);
Population[i].fitness = Evaluate(Population[i]);
}
UpdateBestgeti();
ShuchuPopulation(t);
ShuchuBestgeti(t);
for (int i = 0; i < Pop_Size; i++)
{
outfile << "第"<<t<<"代种群中第"<<i + 1<<"个鲸鱼"<<"->";
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) outfile << Population[i].X[j] <<")对应的适应度为: "<<Evaluate( Population[i])<< endl;
else if(j==0)
outfile <<"("<< Population[i].X[j] << ", ";
else
outfile << Population[i].X[j] << ", ";
}
}
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) outfile << Bestgeti.X[j] << ")" << "对应的适应度为:" << Evaluate(Bestgeti) << endl;
else if (j == 0) outfile << "第"<<t<<"代种群最优鲸鱼个体(" << Bestgeti.X[j] << ",";
else outfile << Bestgeti.X[j] << ",";
}
}
cout<<"****************迭代结束!****************"<<endl;
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) std::cout << Bestgeti.X[j] << ")" << "对应的适应度为:" << Evaluate(Bestgeti) << endl;
else if (j == 0) std::cout << Max_iter<<"次迭代后最终得到的最优鲸鱼个体(" << Bestgeti.X[j] << ",";
else std::cout << Bestgeti.X[j] << ",";
}
outfile<<"****************迭代结束!****************"<<endl;
for (int j = 0; j < dimension; j++)
{
if (j == dimension - 1) outfile << Bestgeti.X[j] << ")" << "对应的适应度为:" << Evaluate(Bestgeti) << endl;
else if (j == 0) outfile << Max_iter<<"次迭代后最终得到的最优鲸鱼个体(" << Bestgeti.X[j] << ",";
else outfile << Bestgeti.X[j] << ",";
}
outfile.close();
}
};
int main()
{
srand((int)time(0));
std::cout << "****************鲸鱼优化算法求解旅行商问题!****************" << endl;
WOA woa;
woa.WOA_TSP(50,29,1,29,100,"D:\\program files\\CodeBlocks\\myworkspace\\WOA_TSP\\bayg29.tsp");
return 0;
}
3.1.2 运行结果
CodeBlocks 控制台运行结果
文件夹下生成的文件
3.2 WOA求解TSP(Python)
????????这里用到的Python 环境为Python 3.8 ,安装的依赖包有numpy 和matplotlib ,使用的IDE为PyCharm 。
Python 3.8
PyCharm 社区版
3.2.1 Python代码(WOA.py)
import random
import numpy as np
import math
from matplotlib import pyplot as plt
class City(object):
def __init__(self, ID, CoordinateX, CoordinateY):
self.ID = ID
self.CoordinateX = CoordinateX
self.CoordinateY = CoordinateY
def shuchu(self):
print(self.ID, self.CoordinateX, self.CoordinateY)
def getDistanceOfTwoCity(city1, city2):
return math.sqrt((pow((city1.CoordinateX - city2.CoordinateX), 2) + pow((city1.CoordinateY - city2.CoordinateY), 2)) / 10.0)
def ReplaceContinueousSpace(str):
n = len(str)
newstr = ''
for i in range(n):
if str[i] != ' ':
newstr += str[i]
if i+1 < n and str[i+1] == ' ':
newstr += ','
return newstr
Distances = []
class CityGraph(object):
def __init__(self, filename):
self.Citys = []
self.n = 0
self.readDataFile(filename)
self.computeDistances()
def readDataFile(self, filename):
f = open(filename)
line = f.readline()
while line:
line_str = line.strip('\n')
city_linestr = ReplaceContinueousSpace(line_str).split(',')
city_i = City(int(city_linestr[0]), float(city_linestr[1]), float(city_linestr[2]))
self.Citys.append(city_i)
line = f.readline()
f.close()
self.n = len(self.Citys)
def computeDistances(self):
for i in range(0, self.n):
dist_i = []
for j in range(0, self.n):
dist_ij = getDistanceOfTwoCity(self.Citys[i], self.Citys[j])
dist_i.append(dist_ij)
Distances.append(dist_i)
def shuchu(self):
print('城市ID-----城市横坐标X-----城市纵坐标Y')
for i in range(self.n):
self.Citys[i].shuchu()
print('----------各城市间的距离矩阵---------')
for i in range(0,self.n):
for j in range(0, self.n):
print(round(Distances[i][j],3), end=',')
print('\n')
class WOA(object):
def __init__(self, population_size=50, dimension=29, lb=1, ub=29, Max_iteration = 100):
self.dimension = dimension
self.lb = lb
self.ub = ub
self.PopSize = population_size
self.Max_iter = Max_iteration
self.Population = []
self.BestGeti = []
self.BestFitness = -100
self.Convergence_curve = []
def Init_Population(self):
self.Population = []
for i in range(0, self.PopSize):
geti_i = []
for j in range(0, self.dimension):
xj = j + 1
geti_i.append(xj)
random.shuffle(geti_i)
self.Population.append(geti_i)
self.BestGeti = self.Population[0]
self.BestFitness = self.Evaluate(self.Population[0])
for i in range(0,self.PopSize):
fitness_i = self.Evaluate(self.Population[i])
print('种群初始化第 {} 个服务链座鱼鲸个体 {} 的适应度为 {}'.format(i + 1, self.Population[i], fitness_i))
if fitness_i < self.BestFitness:
self.BestGeti = self.Population[i]
self.BestFitness = fitness_i
print('\n初始化种群最优服务链座鱼鲸个体 {} 的适应度为 {}'.format(self.BestGeti, self.BestFitness))
def Valid_Geti(self,geti):
for j in range(0, self.dimension):
if geti[j] > self.lb and geti[j] < self.ub:
geti[j] = int(math.floor(geti[j]))
else:
geti[j] = random.randint(int(self.lb), int(self.ub))
n = len(geti)
std_geti = []
std_flag = []
Count = []
for i in range(0, n):
std_geti.append(i + 1)
Count.append(0)
std_flag.append(True)
for i in range(0, n):
findNum = std_geti[i]
count = 0
for j in range(0, n):
if (findNum == geti[j]):
count = count + 1
if count > 1:
std_flag[j] = False
Count[i] = count
for i in range(0, n):
if Count[i] == 0:
for j in range(0, n):
if std_flag[j] == False:
geti[j] = std_geti[i]
std_flag[j] = True
break
def Evaluate(self, geti):
sumDist = 0
n = len(geti)
for i in range(0, n-1):
sumDist += Distances[geti[i]-1][geti[i+1]-1]
sumDist += Distances[geti[n-1]-1][geti[0]-1]
return int(sumDist)
def Start(self):
self.Init_Population()
t = 0
while t < self.Max_iter:
a = 2 - t * (2 / self.Max_iter)
a2 = -1 + t * ((-1) / self.Max_iter)
for i in range(0, self.PopSize):
r1 = random.random()
r2 = random.random()
A = 2 * a * r1 - a
C = 2 * r2
b = 1
l = (a2 - 1) * random.random() +1
p = random.random()
for j in range(0,self.dimension):
if p < 0.5:
if abs(A) >= 1:
rand_leader_index = math.floor(self.PopSize * random.random())
X_rand = self.Population[rand_leader_index]
D_X_rand = abs(C * X_rand[j] - self.Population[i][j])
self.Population[i][j] = X_rand[j] - A * D_X_rand
else:
D_Leader = abs(C * self.BestGeti[j] - self.Population[i][j])
self.Population[i][j] = self.BestGeti[j] - A * D_Leader
else:
distance2Leader = abs(self.BestGeti[j] - self.Population[i][j])
self.Population[i][j] = distance2Leader * np.exp(b * l) * np.cos(l * 2 * math.pi) + self.BestGeti[j]
t = t + 1
self.Convergence_curve.append(self.BestFitness)
for i in range(0,self.PopSize):
self.Valid_Geti(self.Population[i])
fitness_i = self.Evaluate(self.Population[i])
if fitness_i < self.BestFitness:
self.BestGeti = self.Population[i]
self.BestFitness = fitness_i
print('第 {} 代种群第 {} 个座鱼鲸个体 {} 的适应度为 {}'.format(t + 1, i+1, self.Population[i], fitness_i))
print('第 {} 代种群 最优座鱼鲸个体 {} 的适应度为 {}'.format(t + 1, self.BestGeti, self.BestFitness))
print('{} 次迭代后WOA种群最优座鱼鲸个体 {} 的适应度值为 {}'.format(self.Max_iter, self.BestGeti, self.BestFitness))
if __name__ == "__main__":
filename = 'bayg29.tsp'
citygraph = CityGraph(filename)
citygraph.shuchu()
woa = WOA(50, 29, 1, 29, 500)
woa.Start()
epochs = range(len(woa.Convergence_curve))
plt.figure()
plt.plot(epochs, woa.Convergence_curve, "b", label="woa")
plt.title('WOA_ServiceChain')
plt.xlabel("Epochs")
plt.ylabel("Fitness")
plt.legend()
plt.show()
3.2.2 运行结果
PyCharm控制台运行结果
4、总结
Web of Science官网与WOA主题相关的论文
WOA论文中对元启发式算法的分类
????????WOA 作为一种新型的启发式优化算法,通过对座头鲸的捕食行为进行数学建模,螺旋式上升 和收缩包围 机制易于理解,算法的执行时间短,求解速度快且效率高,在云计算、无线传感网、图像处理、交通、电力等领域都有了较为广泛的应用,对于求解旅行商、资源分配、组合优化等问题十分有效。
|