各位小伙伴可在闲鱼搜索 优化算法交流地，即可搜索到官方闲鱼账号，谨防上当受骗。
hello，大家立秋快乐！今天为各位免费分享头脑风暴优化（BSO）算法求解旅行商问题（TSP）（附MATLAB代码）这篇推文的源代码。
这份源代码包含7个函数，分别如下：
01 | 主函数
主函数的输入是文本文件pr226.txt（第1列是序号，第2列是x坐标，第3列是y坐标），输出是最优路线。文本文件可根据自己需要进行替换，只要保持3列的这种形式即可，文本文件格式如下。
序号 | x坐标 | y坐标 |
%
% @作者：随心390
% @微信公众号：优化算法交流地
%
tic
clear
clc
%% 导入数据
pr226=importdata('pr226.txt');
N=size(pr226,1); %城市数目
vertexs=pr226(:,2:3); %各点xy坐标
x=vertexs(:,1); %x坐标
y=vertexs(:,2); %y坐标
h=pdist(vertexs);
dist=squareform(h); %距离矩阵
%% 参数初始化
MAXGEN=1000; %最大迭代次数
NIND=50; %种群数目
cluster_num=5; %聚类数目
p_replace=0.1; %用随机解替换一个聚类中心的概率
p_one=0.5; %选择1个聚类的概率
p_two=1-p_one; %选择2个聚类的概率，p_two=1-p_one
p_one_center=0.3; %选择1个聚类中聚类中心的概率
p_two_center=0.2; %选择2个聚类中聚类中心的概率
%% 种群初始化
Population=InitPop(NIND,N);
%% 主循环
gen=1; %计数器初始化
BestPop=zeros(MAXGEN,N); %记录每次迭代过程中全局最优个体
BestL=zeros(MAXGEN,1); %记录每次迭代过程中全局最优个体的总距离
while gen<=MAXGEN
%% 计算目标函数值
Obj=ObjFunction(Population,dist);
%% K-means聚类
[Idx,C,sumD,D]=kmeans(Obj,cluster_num,'Distance','cityblock','Replicates',2);
cluster=cell(cluster_num,2); %将解储存在每一个聚类中
order_cluster=cell(cluster_num,2); %将储存在每一个聚类中的解按照目标函数值排序
for i=1:cluster_num
cluster{i,1}=Population(Idx==i,:); %将个体按照所处的聚类编号储存到对应的聚类中
cluster_row(i)=size(cluster{i,1},1); %计算当前聚类中个体数目
for j=1:cluster_row(i)
Individual=cluster{i,1}(j,:); %当前聚类中第j个个体
cluster{i,2}(j,1)=ObjFunction(Individual,dist); %计算当前聚类中第j个个体的目标函数值
end
[order_cluster{i,2},order_index]=sort(cluster{i,2}) ; %将当前聚类中的所有个体按照目标函数值从小到大的顺序进行排序
order_cluster{i,1}=cluster{i,1}(order_index,:); %将当前聚类中的所有个体按照排序结果重新排列
order_index=0; %重置排序序号
end
cluster_fit=cell2mat(order_cluster); %将聚类的元胞数组转换为矩阵，最后一列为个体的目标函数值
%% 以一定的概率随机从m个聚类中心中选择出一个聚类中心，并用一个新产生的随机解更新这个被选中的聚类中心
R1=rand(1,1);
if R1<=p_replace
%随机选择一个聚类中心
repalce_cluster_num=randi([1,cluster_num],1,1);
%随机产生一个解
replace_solution=randperm(N);
%并用这个新产生的随机解更新这个被选中的聚类中心
order_cluster{repalce_cluster_num,1}(1,:)=replace_solution;
%计算新解的目标函数值
replace_solution_fitness=ObjFunction(replace_solution,dist);
%将新解的目标函数值储存到order_cluster中
order_cluster{repalce_cluster_num,2}(1,:)=replace_solution_fitness;
end
%% 更新这n个个体
for i=1:NIND
%如果随机数小于选择1个聚类的概率，则随机选择1个聚类
if rand()<p_one
select_one_cluster=randi([1,cluster_num],1,1); %选择选择一个聚类
%如果随机数小于选择1个聚类中聚类中心的概率，或当前聚类中只有一个个体
if rand()<p_one_center||cluster_row(select_one_cluster)==1
select_ind=order_cluster{select_one_cluster,1}(1,:); %选择当前解的聚类中心（只有一个个体时，就选择该个体）
else
r_1= randi([2,cluster_row(select_one_cluster)],1,1); %随机选择当前聚类中除聚类中心外的其它个体序号
select_ind=order_cluster{select_one_cluster,1}(r_1,:); %随机选择当前聚类中除聚类中心外的其它个体
end
indi_temp=Swap(select_ind); %将该个体进行交换操作
else %如果随机数不小于选择1个聚类的概率，则随机选择2个聚类
cluster_two=[0,0]; %随机产生两个聚类的序号
while cluster_two(1,1)==cluster_two(1,2)
cluster_two=randi([1,cluster_num],1,2); %这两个聚类序号不能相同
end
%如果随机数小于选择2个聚类中聚类中心的概率，或当前两个聚类中都只有一个个体
if (rand()<p_two_center)||(cluster_row(cluster_two(1,1))==1&&cluster_row(cluster_two(1,2))==1)
%选择这2个聚类中聚类中心
select_ind1=order_cluster{cluster_two(1,1),1}(1,:); %第一个被选择的聚类中心
select_ind2=order_cluster{cluster_two(1,2),1}(1,:); %第二个被选择的聚类中心
else
%如果第1个选择的聚类中只有一个个体
if cluster_row(cluster_two(1,1))==1
r_2=randi([2,cluster_row(cluster_two(1,2))],1,1); %选择第2个聚类中除聚类中心外的其它个体
select_ind1=order_cluster{cluster_two(1,1),1}(1,:); %第一个聚类中被选择的个体
select_ind2=order_cluster{cluster_two(1,2),1}(r_2,:); %第二个聚类中被选择的个体
elseif cluster_row(cluster_two(1,2))==1 %如果第2个选择的聚类中只有一个个体
r_3=randi([2,cluster_row(cluster_two(1,1))],1, 1); %选择第1个聚类中除聚类中心外的其它个体
select_ind1=order_cluster{cluster_two(1,1),1}(r_3,:); %第一个聚类中被选择的个体
select_ind2=order_cluster{cluster_two(1,2),1}(1,:); %第二个聚类中被选择的个体
elseif cluster_row(cluster_two(1,1))>1&&cluster_row(cluster_two(1,2))>1
%%选择这2个聚类中除聚类中心外的其它个体
r_4=randi([2,cluster_row(cluster_two(1,1))],1,1);
r_5=randi([2,cluster_row(cluster_two(1,2))],1,1);
select_ind1=order_cluster{cluster_two(1,1),1}(r_4,:); %第一个聚类中被选择的个体
select_ind2=order_cluster{cluster_two(1,2),1}(r_5,:); %第二个聚类中被选择的个体
end
end
[child1,child2,min_index,start_c]=heuristic_crossover(select_ind1,select_ind2,dist); %启发式交叉算子
%如果子代个体1目标函数值更小，则将indi_temp赋值为child1，否则赋值为child2
if min_index==1
indi_temp=child1;
else
indi_temp=child2;
end
end
fit_indi_temp=ObjFunction(indi_temp,dist);
%如果fit_indi_temp比原来位置上的个体目标函数值更小，则更新该位置上的个体
if fit_indi_temp<order_cluster{Idx(i),2}(1,:)
Population(i,:) =indi_temp(1,:);
end
end
%% 计算目标函数值
Obj=ObjFunction(Population,dist);
[min_len,min_index]=min(Obj); %当前种群中最优个体以及所对应的序号
BestL(gen,1)=min_len; %记录当前迭代中最优个体
BestPop(gen,:)=Population(min_index,:); %记录当前迭代中最优个体的总距离
%% 打印各代最优解
disp(['第',num2str(gen),'代最优个体的总距离为',num2str(min_len)]);
%% 计数器加1
gen=gen+1;
end
%% 计算最终种群的目标函数值
Obj=ObjFunction(Population,dist);
[best_len,best_index]=min(Obj);
bestR=Population(best_index,:); %全局最优个体
PlotRoute(bestR,x,y);
%% 打印每次迭代的全局最优个体的总距离变化趋势图
figure;
plot(BestL,'LineWidth',1);
title('优化过程')
xlabel('迭代次数');
ylabel('总距离');
toc
02 | 种群初始化函数
种群初始化函数InitPop的输入是种群大小NIND、城市数目N，输出是随机生成的初始种群Population。
%
% @作者：随心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
03 | 目标函数
目标函数ObjFunction的输入是种群Population、距离矩阵dist，输出是每个个体的目标函数值Obj，即每个个体的总距离。
%
% @作者：随心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 | 路线长度计算函数
路线长度计算函数RouteLength的输入是路线route、距离矩阵dist，输出是该路线总距离L。
%
% @作者：随心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
05 | 交换函数
交换函数Swap的输入是路线1 route1，输出是经过交换结构变换后的路线2 route2。
%
% @作者：随心390
% @微信公众号：优化算法交流地
%
%比如说有6个城市，当前解为123456，我们随机选择两个位置，然后将这两个位置上的元素进行交换。
%比如说，交换2和5两个位置上的元素，则交换后的解为153426。
%输入route1： 路线1
%输出route2： 经过交换结构变换后的路线2
function route2=Swap(route1)
n=length(route1);
seq=randperm(n);
I=seq(1:2);
i1=I(1);
i2=I(2);
route2=route1;
route2([i1 i2])=route1([i2 i1]);
end
06 | 启发式交叉函数
启发式交叉函数heuristic_crossover的输入是父代个体1 parent1、父代个体2 parent2、距离矩阵 dist，输出是子代个体1 child1、子代个体2 child2、目标函数值更小的个体序号min_index、初始选择一个城市作为子代个体的起点start_c。
%
% @作者：随心390
% @微信公众号：优化算法交流地
%
%% 启发式交叉算子
%输入parent1： 父代个体1
%输入parent2： 父代个体2
%输入dist： 距离矩阵
%输出child1： 子代个体1
%输出child2： 子代个体2
%输出min_index： 目标函数值更小的个体序号
%输出start_c： 初始选择一个城市作为子代个体的起点
function [child1,child2,min_index,start_c]=heuristic_crossover(parent1,parent2,dist)
n=numel(parent1); %城市数目
child1=zeros(1,n); %初始化子代个体1
child2=zeros(1,n); %初始化子代个体2
start=randperm(n,1); %随机选择一个城市作为子代个体的起点
start_c=start; %复制start
child1(1)=start; %更新子代个体1的起点为start
child2(1)=start; %更新子代个体2的起点为start
%% 终止条件是parent1_c1只剩一个点时
parent1_c1=parent1; %复制parent1
parent2_c1=parent2; %复制parent2
pos=2; %记录child中的位置
while numel(parent1_c1)~=1
start1=find(parent1_c1==start,1,'first'); %找出start在parent1的位置
start2=find(parent2_c1==start,1,'first'); %找出start在parent2的位置
%% 找出start在parent1_c1的右紧邻点
if start1==numel(parent1_c1)
right1=parent1_c1(1); %如果start1是parent1最后一个点，则start1右边紧邻的点是parent1第一个点
else
right1=parent1_c1(start1+1); %start1右边紧邻的点
end
%% 找出start在parent2_c1的右紧邻点
if start2==numel(parent2_c1)
right2=parent2_c1(1); %如果start2是parent2最后一个点，则start2右边紧邻的点是parent2第一个点
else
right2=parent2_c1(start2+1); %start2右边紧邻的点
end
%% 比较dist(start,right1)和dist(start,right2)，目的是更新chlid1
if dist(start,right1)<=dist(start,right2)
%如果dist(start,right1)<=dist(start,right2)，则说明start离right1更近
%则将child1的下一个位置更新为right1
child1(pos)=right1;
else
%如果dist(start,right1)>dist(start,right2)，则说明start离right2更近
%则将child1的下一个位置更新为right2
child1(pos)=right2;
end
%% 更新parent1_c1和parent2_c1
parent1_c1(parent1_c1==start)=[];
parent2_c1(parent2_c1==start)=[];
%% 更新start
start=child1(pos);
%% 更新pos
pos=pos+1;
end
%% 终止条件是parent1_c2只剩一个点时
parent1_c2=parent1; %复制parent1
parent2_c2=parent2; %复制parent2
start=start_c; %更新start为初始start
pos=2; %记录child中的位置
while numel(parent1_c2)~=1
start1=find(parent1_c2==start,1,'first'); %找出start在parent1的位置
start2=find(parent2_c2==start,1,'first'); %找出start在parent2的位置
%% 找出start在parent1_c1的左紧邻点
if start1==1
left1=parent1_c2(end); %如果start1是parent1第一个点，则start1左边紧邻的点是parent1最后一个点
else
left1=parent1_c2(start1-1); %start1左边紧邻的点
end
%% 找出start在parent2_c1的左紧邻点
if start2==1
left2=parent2_c2(end); %如果start2是parent2第一个点，则start2左边紧邻的点是parent2最后一个点
else
left2=parent2_c2(start2-1); %start2左边紧邻的点
end
%% 比较dist(left1,start)和dist(left2.start)，目的是更新chlid2
if dist(left1,start)<=dist(left2,start)
%如果dist(left1,start)<=dist(left2.start)，则说明start离left1更近
%则将child2的下一个位置更新为left1
child2(pos)=left1;
else
%如果dist(left1,start)>dist(left2.start)，则说明start离left2更近
%则将child2的下一个位置更新为left2
child2(pos)=left2;
end
%% 更新parent1_c2和parent2_c2
parent1_c2(parent1_c2==start)=[];
parent2_c2(parent2_c2==start)=[];
%% 更新start
start=child2(pos);
%% 更新pos
pos=pos+1;
end
len1=ObjFunction(child1,dist);
len2=ObjFunction(child2,dist);
[~,min_index]=min([len1,len2]);
end
07 | 画图函数
画图函数PlotRoute的输入是路线 route、横纵坐标x和y，输出是最优路线图。
%
% @作者：随心390
% @微信公众号：优化算法交流地
%
%输入route： 路线
%输入x,y： x,y坐标
function PlotRoute(route,x,y)
route=[route route(1)];
plot(x(route),y(route),'k-o','MarkerSize',10,'MarkerFaceColor','w','LineWidth',1.5);
xlabel('x');
ylabel('y');
end
最后，各位小伙伴有什么疑问可点击 写留言 小程序为我们留言，我们会尽可能为各位解答。
