hello,大家好,今天为各位讲解使用头脑风暴优化(BSO)算法求解旅行商问题(TSP)。前期在头脑风暴优化(BSO)算法(附MATLAB代码)这篇推文中讲解了BSO的基本思想,在一种新颖的交叉算子——头脑风暴优化(BSO)算法求解旅行商问题(TSP)预备知识这篇推文中讲解了使用BSO求解TSP的预备知识。
今天的推文在上述两篇推文基础上进一步深入讲解如何使用BSO求解TSP?
一 | BSO基本思想回顾
我们已经学习过BSO是和遗传算法类似的群智能优化算法,但是BSO有其独特的更新解的方式。假设现在种群数目为NIND,那么BSO更新解的方法如下所述:
STEP1:将种群中的NIND个个体分为m类,标准的BSO是采用k-means聚类的方法进行分类的。在matlab可以使用kmeans函数来实现k-means聚类,使用方法如下:
[idx,C,sumD,D] = kmeans(X,K,'Distance','cityblock','MaxIter',1000,'Replicates',5);
其中输入数据:
X——一个N行D列的矩阵,N是数据个数,D是数据维数
K——将X划分为几类
Distance——选择计算一个个体与聚类中心之间距离的方式(默认为sqeuclidean)
cityblock——计算一个个体与聚类中心之间距离的一种方式
MaxIter——最大迭代次数(默认为100)
Replicates——使用新的初始群集质心位置重复群集的次数(默认为1)
输出数据:
idx——N行1列的矩阵,存储的是每个点的聚类标号
C——K行D列的矩阵,存储K个聚类中心
sumD——K行1列的矩阵,存储当前聚类中,所有点与这个聚类中心的距离之和
D——N行D列的矩阵,存储每个点与所有聚类中心的距离
既然标准的BSO采用的是k-means方法,那么这里的分类方法其实是可以进行改进的。本次推文使用的是k-means方法,如果有小伙伴想学习改进的分类方法可以在本次推文的留言小程序给我们留言。
STEP2:将这m类中的每一类中的个体按照目标函数值进行排序,并将每一类中目标函数值最优的那个个体作为聚类中心(PS:虽然在k-means聚类的时候已经选好聚类中心,但是我们并没有使用,而是重新选择聚类中心),同时存储排序后的每一类中的个体及其对应的目标函数值。
STEP3:以一定的概率随机从m个聚类中心中选择出一个聚类中心,并用一个新产生的随机解更新这个被选中的聚类中心。
STEP4:遍历这NIND个个体,想办法使这个NIND个个体向目标函数值更好的方向更新,具体更新方法如下。
在介绍更新方法,先进行定义。
rand——0~1之间的随机数
p_one——选择1个聚类的概率,0~1之间的数
p_two——选择2个聚类的概率,p_two=1-p_one
p_one_center——选择1个聚类中聚类中心的概率,0~1之间的数
p_two_center——选择2个聚类中聚类中心的概率,0~1之间的数
Xselect——选出的个体
omega1——0~1之间的数
omega2——0~1之间的数
因此,更新这NIND个个体的方法有两种:
如果rand<p_one,
1)随机选择1个聚类:
如果rand<p_one_center
a)Xselect=这个聚类中的聚类中心
否则
b) Xselect=从这个聚类中随机选出一个个体
否则
2)随机选择2个聚类:
如果rand<p_two_center
a)Xselect=omega1*聚类1中的聚类中心+omega2*聚类2中的聚类中心
否则
b) Xselect=omega1*聚类1中随机选出的个体1+omega2*聚类2中随机选出的个体2
经过上述选择过程后,得到Xselect,然后使用如下公式对进行Xselect更新。
其中D表示的是数据的维数,normrnd在matlab中是产生正太分布的随机数的函数,logsig在matlab中是逻辑回归中的sigmoid函数。
STEP5:将Xnew与个体i进行比较,如果Xnew的目标函数值优于个体i的目标函数值,则更新个体i=Xnew。否则,不更新个体i。
STEP6:当迭代到最大迭代次数后,迭代过程终止。
二 | 如何设计求解TSP的BSO算法
01 | 种群初始化
BSO求解TSP的种群初始化方法和遗传算法求解TSP的种群初始化方法是相同的。如果城市数目为5,那么BSO求解TSP中的一个个体就可以表示为1~5的一个随机排列,在matlab中用randperm(5)这个函数。
现假设种群数目为NIND,城市数目为N。那么种群初始化函数如下:
%
% @作者:随心390
% @微信公众号:优化算法交流地
%
%% 种群初始化
%输入NIND:种群大小
%输入N:城市数目
%输出Population:随机生成的初始种群
function Population=InitPop(NIND,N)
Population=zeros(NIND,N); %种群初始化为NIND行N列的零矩阵
for i=1:NIND
Population(i,:)=randperm(N); %每个个体为1~N的随机排列
end
end
02 | 计算个体的目标函数值
个体的目标函数值就是这个个体的总距离,具体代码如下:
%
% @作者:随心390
% @微信公众号:优化算法交流地
%
%输入route:路线
%输入dist:距离矩阵
%输出L:该路线总距离
function L=RouteLength(route,dist)
n=length(route);
route=[route route(1)];
L=0;
for k=1:n
i=route(k);
j=route(k+1);
L=L+dist(i,j);
end
end
03 | 计算种群的目标函数值
上述代码是求解一个个体的目标函数值,求解种群中所有个体的目标函数值使用的代码如下:
%
% @作者:随心390
% @微信公众号:优化算法交流地
%
%% 计算种群目标函数值,即每个个体的总距离
%输入Population:种群
%输入dist:距离矩阵
%输出Obj:每个个体的目标函数值,即每个个体的总距离
function Obj=ObjFunction(Population,dist)
NIND=size(Population,1); %种群大小
Obj=zeros(NIND,1); %目标函数初始化为0
for i=1:NIND
route=Population(i,:); %当前个体
Obj(i,1)=RouteLength(route,dist); %计算当前个体的总距离
end
end
04 | 使用kmeans函数将种群中的个体进行分类
kmeans函数的第一个输入参数是要聚类的对象。就TSP而言,如果将一个种群所包含的NIND个个体全作为kmeans函数的第一个输入参数,那么会产生一个问题,最终的分类结果究竟代表什么含义呢?举个例子,如果城市数目为3,假设此时种群有3个个体,分别为:
个体1:1 2 3
个体2:1 3 2
个体3:2 1 3
只看这3个个体,很难说那几个个体是一类的。
为了解决这一问题,我们将NIND个个体的目标函数值作为kmeans函数的第一个输入参数,因为目标函数值是一维的,很容易看出哪一段聚集。
05 | 更新个体的方法
STEP2:将这m类中的每一类中的个体按照目标函数值进行排序,并将每一类中目标函数值最优的那个个体作为聚类中心(PS:虽然在k-means聚类的时候已经选好聚类中心,但是我们并没有使用,而是重新选择聚类中心),同时存储排序后的每一类中的个体及其对应的目标函数值。
这一步保持不变。
STEP3:以一定的概率随机从m个聚类中心中选择出一个聚类中心,并用一个新产生的随机解更新这个被选中的聚类中心。
这里用randperm(N)函数来产生新的随机解。
对STEP4进行较大改动,具体如下:
STEP4:在介绍更新方法,先进行定义。
rand——0~1之间的随机数
p_one——选择1个聚类的概率,0~1之间的数
p_two——选择2个聚类的概率,p_two=1-p_one
p_one_center——选择1个聚类中聚类中心的概率,0~1之间的数
p_two_center——选择2个聚类中聚类中心的概率,0~1之间的数
Xselect——选出的个体
Xnew——将选出的个体Xselect进行交换或交叉操作后得到的新个体
更新这NIND个个体的方法有两种:
如果rand<p_one,
1)随机选择1个聚类:
如果rand<p_one_center
a)Xselect=这个聚类中的聚类中心
否则
b) Xselect=从这个聚类中随机选出一个个体
对Xselect进行交换操作得到Xnew,交换操作即交换两个位置上的城市。
否则
2)随机选择2个聚类:
如果rand<p_two_center
a)令select_ind1=聚类1中的聚类中心;令select_ind2=聚类2中的聚类中心;
Xselect=omega1*聚类1中的聚类中心+omega2*聚类2中的聚类中心
否则
b) 令select_ind1=聚类1中随机选出的个体1;令select_ind2=聚类2中随机选出的个体2;
对select_ind1和select_ind2进行交叉操作,得到两个个体,将目标函数值更小的那个个体赋值给Xnew。关于交叉操作的详细内容可以看一种新颖的交叉算子——头脑风暴优化(BSO)算法求解旅行商问题(TSP)预备知识这篇推文。
STEP5:将Xnew与个体i进行比较,如果Xnew的目标函数值优于个体i的目标函数值,则更新个体i=Xnew。否则,不更新个体i。
STEP6:当迭代到最大迭代次数后,迭代过程终止。
三 | BSO求解TSP效果展示
代码中求解的是TSP标准测试算例中pr226算例,TSP标准测试算例网址为:
http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/tsp/
各个算例所对应的最优距离网址为:
http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/STSP.html
pr226算例目前最优距离为80369。
我们将算法的种群数目设为50,迭代次数设为1000。
最终得到的结果如下:
pr226最终路线图
pr226优化过程图
知乎 | bilibili:随心390