【TSP】基于matlab蜜蜂算法求解旅行商问题【含matlab源码 1248期】

一、获取代码方式

获取代码方式1：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

获取代码方式2：
完整代码已上传我的资源：【TSP】基于matlab蜜蜂算法求解旅行商问题【含matlab源码 1248期】

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、TSP简介

旅行商问题，即TSP问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。
TSP的数学模型

三、蜜蜂算法简介

1 什么是人工蜂群算法？
人工蜂群算法是模仿蜜蜂行为提出的一种优化方法，是集群智能思想的一个具体应用，它的主要特点是不需要了解问题的特殊信息，只需要对问题进行优劣的比较，通过各人工蜂个体的局部寻优行为，最终在群体中使全局最优值突现出来，有着较快的收敛速度。为了解决多变量函数优化问题，Karaboga提出了人工蜂群算法ABC模型(artificial bee colony algorithm)。
人工蜂群算法属于群智能算法的一种，其受启发于蜜蜂的寻蜜和采蜜过程，相比于常见的启发式算法，它的优点在于其使用了较少的控制参数，并且鲁棒性强，在每次迭代过程中都会进行全局和局部的最优解搜索，因此能够找到最优解的概率大大增加。

相比于遗传算法来说，人工蜂群算法在局部的收敛和寻优能力上要更为出色，不会出现遗传算法的“早熟”现象，并且算法的复杂度也较低。但由于遗传算法有交叉以及变异的操作，因此遗传算法在全局最优值的搜索上要优于人工蜂群算法。此外，人工蜂群算法适用于进行连续函数的全局优化问题，而不适用于一些离散函数。

2 蜜蜂的采蜜机制
蜜蜂具有群智能应必备的两个条件：自组织性和分工合作性。虽然单个蜜蜂的行为很简单，但是由单个蜜蜂所组成的群体却能够表现出极其复杂的行为，它们可以在任何复杂的环境下以很高的效率从花朵中采集花蜜，同时还能够很快的适应环境的改变。

任务分工
人工蜂群算法就是模拟蜜蜂的采蜜过程而提出的一种新型智能优化算法，它也是由食物源、雇佣蜂和非雇佣蜂三部分组成。

食物源：食物源即为蜜源。在任何一个优化问题中，问题的可行解都是以一定形式给出的。在人工蜂群算法中，食物源就是待求优化问题的可行解，是人工蜂群算法中所要处理的基本对象。食物源的优劣即可行解的好坏是用蜜源花蜜量的大小即适应度来评价的。

雇佣蜂：雇佣蜂即为引领蜂(采蜜蜂)与食物源的位置相对应，一个食物源对应一个引领蜂。在人工蜂群算法中，食物源的个数与引领蜂的个数相等；引领蜂的任务是发现食物源信息并以一定的概率与跟随蜂分享；概率的计算即为人工蜂群算法中的选择策略，一般是根据适应度值以轮盘赌的方法计算。

非雇佣蜂：非雇佣蜂包括跟随蜂(观察蜂)和侦査蜂跟随蜂在蜂巢的招募区内根据引领蜂提供的蜜源信息来选择食物源，而侦查蜂是在蜂巢附近寻找新的食物源。在人工蜂群算法中，跟随蜂依据引领蜂传递的信息，在食物源附近搜索新食物源，并进行贪婪选择。若一个食物源在经过次后仍未被更新，则此引领蜂变成侦査蜂，侦查蜂寻找新的食物源代替原来的食物源。

采蜜机制
在蜜蜂群体智能形成的过程中，蜜蜂之间的信息交流是最重要的环节，而舞蹈区是蜂巢中最重要的信息交换地。采蜜蜂在舞蹈区通过跳摇摆舞与其他蜜蜂共同分享食物源的信息，观察蜂则是通过采蜜蜂所跳的摇摆舞来获得当前食物源的信息的，所以，观察蜂要以最小的资源耗费来选择到哪个食物源采蜜。因此，蜜蜂被招募到某个食物源的概率与食物源的收益率成正比。

初始时刻，蜜蜂的搜索不受任何先验知识的决定，是完全随机的。此时的蜜蜂有以下两种选择：
它转变成为侦察蜂，并且由于一些内部动机或可能的外部环境自发地在蜂巢附近搜索食物源；
在观看了摇摆舞之后，它可能被招募到某个食物源，并且开始开采食物源。
在蜜蜂确定食物源后，它们利用自己本身的存储能力来记忆位置信息并开始采集花蜜。此时，蜜蜂将转变为“雇佣蜂”。蜜蜂在食物源处采集完花蜜，回到蜂巢并卸下花蜜后有如下选择：
放弃食物源成为非雇佣蜂；
跳摇摆舞为所对应的食物源招募更多的蜜蜂，然后回到食物源采蜜；
继续在同一食物源采蜜而不进行招募。
蜜蜂在采蜜时所表现出来的这种自组织性和合理分配性主要由其自身的基本性质所决定的，它们所特有的基本性质如下：
正反馈性：食物源的花蜜量与食物源被选择的可能性成正比；
负反馈性：蜜蜂停止对较差食物源的开采过程；
波动性：在某个食物源被放弃时，随机搜索一个食物源替代原食物源；
互动性：蜜蜂在舞蹈区与其他蜜蜂共同分享食物源的相关信息。

3 蜂群算法的原理及流程
标准的ABC算法通过模拟实际蜜蜂的采蜜机制将人工蜂群分为3类: 采蜜蜂、观察蜂和侦察蜂。整个蜂群的目标是寻找花蜜量最大的蜜源。在标准的ABC算法中，采蜜蜂利用先前的蜜源信息寻找新的蜜源并与观察蜂分享蜜源信息；观察蜂在蜂房中等待并依据采蜜蜂分享的信息寻找新的蜜源；侦查蜂的任务是寻找一个新的有价值的蜜源，它们在蜂房附近随机地寻找蜜源。

算法原理
假设问题的解空间是D维的，采蜜蜂与观察蜂的个数都是SN，采蜜蜂的个数或观察蜂的个数与蜜源的数量相等。则标准的ABC算法将优化问题的求解过程看成是在D维搜索空间中进行搜索。每个蜜源的位置代表问题的一个可能解，蜜源的花蜜量对应于相应的解的适应度。一个采蜜蜂与一个蜜源是相对应的。与第i个蜜源相对应的采蜜蜂依据如下公式寻找新的蜜源：

算法流程
人工蜂群算法具体实现步骤：
初始化种群解；
计算种群中各个蜜蜂的适应值；
重复一下步骤，知道满足终止条件：
采蜜蜂根据(1)产生新的解并计算适应值；
采蜜蜂根据贪心策略选择蜜源；
根据(2)式计算选择蜜源的概率；
观察蜂根据概率选择蜜源，根据(3)式在该蜜源附近产生新的蜜源，并计算新蜜源的适应值；
观察蜂根据贪心策略选择蜜源；
决定是否存在需要放弃的蜜源，如果存在，则随机产生一个蜜源替代它；
记录最优解直至满足终止条件输出最优解；

四、部分源代码

%蜜蜂求解TSP问题
clear
clc
close all
%% 加载数据
data=load('berlin52.txt');
X = data(:,2:3);
D=Distanse(X);  %生成距离矩阵
N=size(D,1);    %城市个数
%% 算法参数
MAXGEN=600;     %最大迭代次数
NIND=400;       %种群大小
PScouts=0.2;    %侦察蜂比例
NScouts=NIND*PScouts;    %侦察蜂数量
NWorkers=NIND-NScouts;   %工蜂数量
NElites=floor(N/10);     %精英解数量
NbRange=5;       %局部搜索的邻域范围，即可改变的最大位数
NWorkersElites = round(NWorkers/NElites);%各精英解邻域内局部搜索解的数量%% 初始化种群,即首先执行一次全局搜索
Bees=GlobalSearch(NIND,N);
%% 画出随机解的路径图
DrawPath(Bees(1,:),X)
pause(0.0001)
%% 输出随机解的路径和总距离
disp('初始种群中的一个随机值:')
OutputPath(Bees(1,:));
Rlength=PathLength(D,Bees(1,:));
disp(['总距离：',num2str(Rlength)]);
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')gen=0;
figure;
hold on;box on
xlim([0,MAXGEN])
title('优化过程')
xlabel('代数')
ylabel('最优值')
ObjV=PathLength(D,Bees);  %计算本代个体路径长度
[preObjV,preObjVPos]=min(ObjV);
FitnV=Fitness(ObjV);
%% 选择精英解
Elites=Select(Bees,FitnV,NElites); %赌轮盘选择
Elites(1,:)=Bees(preObjVPos,:);  %确保保留每代最优解while gen<MAXGEN%% 更新迭代次数gen=gen+1 ;%% 局部搜索BeesLocalSearch=LocalSearch(Elites,NbRange,NWorkersElites);%% 全局搜索BeesGlobalSearch=GlobalSearch(NScouts,N); %% 生成新种群Bees=[Elites;BeesLocalSearch;BeesGlobalSearch]; %% 计算适应度ObjV=PathLength(D,Bees);  %计算路径长度fprintf('%d   %1.10f\n',gen,min(ObjV))line([gen-1,gen],[preObjV,min(ObjV)]);
%     pause(0.0001)[preObjV,preObjVPos]=min(ObjV);FitnV=Fitness(ObjV);%% 选择精英解Elites=Select(Bees,FitnV,NElites); %赌轮盘选择Elites(1,:)=Bees(preObjVPos,:);  %确保保留每代最优解
end
%% 画出最优解的路径图
ObjV=PathLength(D,Bees);  %计算路径长度
[minObjV,minInd]=min(ObjV);
DrawPath(Bees(minInd(1),:),X)
%% 输出最优解的路径和总距离
disp('最优解:')
p=OutputPath(Bees(minInd(1),:));
disp(['总距离：',num2str(ObjV(minInd(1)))]);
disp('-------------------------------------------------------------')
%% 交叉操作
% 输入
%SelCh  被选择的个体
%Pc     交叉概率
%输出：
% SelCh 交叉后的个体
function SelCh=Recombin(SelCh,Pc)
NSel=size(SelCh,1);
for i=1:2:NSel-mod(NSel,2)if Pc>=rand %交叉概率Pc[SelCh(i,:),SelCh(i+1,:)]=intercross(SelCh(i,:),SelCh(i+1,:));end
end%输入：
%a和b为两个待交叉的个体
%输出：
%a和b为交叉后得到的两个个体function [a,b]=intercross(a,b)
L=length(a);
r1=randsrc(1,1,[1:L]);
r2=randsrc(1,1,[1:L]);
if r1~=r2a0=a;b0=b;s=min([r1,r2]);e=max([r1,r2]);for i=s:ea1=a;b1=b;a(i)=b0(i);b(i)=a0(i);%解决因交叉引起的城市重复出现问题，逐位修改x=find(a==a(i));y=find(b==b(i));i1=x(x~=i);i2=y(y~=i);if ~isempty(i1)a(i1)=a1(i);endif ~isempty(i2)b(i2)=b1(i);endend
end%
% %交叉算法采用部分匹配交叉%交叉算法采用部分匹配交叉
% function [a,b]=intercross(a,b)
% L=length(a);
% r1=ceil(rand*L);
% r2=ceil(rand*L);
% r1=4;r2=7;
% if r1~=r2
%     s=min([r1,r2]);
%     e=max([r1,r2]);
%     a1=a;b1=b;
%     a(s:e)=b1(s:e);
%     b(s:e)=a1(s:e);
%     for i=[setdiff(1:L,s:e)]
%         [tf, loc] = ismember(a(i),a(s:e));
%         if tf
%             a(i)=a1(loc+s-1);
%         end
%         [tf, loc]=ismember(b(i),b(s:e));
%         if tf
%             b(i)=b1(loc+s-1);
%         end
%     end
% end
%% 初始化种群
%输入：
% NIND：种群大小
% N：   个体染色体长度（这里为城市的个数）
%输出：
%初始种群
function Chrom=GlobalSearch(NIND,N)
Chrom=zeros(NIND,N);%用于存储种群
for i=1:NINDChrom(i,:)=randperm(N);%随机生成初始种群
end%% 重插入子代的新种群%输入：%Chrom  父代的种群%SelCh  子代种群%ObjV   父代适应度%输出% Chrom  组合父代与子代后得到的新种群
function Chrom=Reins(Chrom,SelCh,ObjV)
NIND=size(Chrom,1);
NSel=size(SelCh,1);
[TobjV,index]=sort(ObjV);
Chrom=[Chrom(index(1:NIND-NSel),:);SelCh];
%子代新解全部留下，父代仅留下较好的解
%可以推断，当选择率、交叉率和变异率变小时，
%出现新解的机会变小，则在迭代初期进化变慢。但在
%后期，由于稳定性提高（每一代解中保留上一代优势解的比例较大，
%随机产生的新解的比例变小)，则最终解的收敛速度和质量并不差。
%% 进化逆转函数
%输入
%SelCh 被选择的个体
%D     各城市的距离矩阵
%输出
%SelCh  进化逆转后的个体
function SelCh=Reverse(SelCh,D)
[row,col]=size(SelCh);
ObjV=PathLength(D,SelCh);  %计算路径长度
SelCh1=SelCh;
for i=1:rowr1=randsrc(1,1,[1:col]); %1到col之间随机取一个整数r2=randsrc(1,1,[1:col]);mininverse=min([r1 r2]);maxinverse=max([r1 r2]);SelCh1(i,mininverse:maxinverse)=SelCh1(i,maxinverse:-1:mininverse);
end
ObjV1=PathLength(D,SelCh1);  %计算路径长度
index=ObjV1<ObjV; %找出逆转之后得到改进的方案
SelCh(index,:)=SelCh1(index,:);

五、运行结果

六、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.