蚁群算法及matlab实现TSP问题源代码

1 蚁群算法原理

自1991年由意大利学者 M. Dorigo，V. Maniezzo 和 A. Colorni 通过模拟蚁群觅食行为提出了一种基于种群的模拟进化算法——蚁群优化。该算法的出现引起了学者们的极大关注，蚁群算法的特点：
① 其原理是一种正反馈机制或称增强型学习系统; 它通过【最优路径上蚂蚁数量的增加→信息素强度增加→后来蚂蚁选择概率增大→最优路径上蚂蚁数量更大增加】达到最终收敛于最优路径上L
② 它是一种通用型随机优化方法, 它吸收了蚂蚁的行为特(内在搜索机制) , 它是使用人工蚂蚁仿真(也称蚂蚁系统) 来求解问题L但人工蚂蚁决不是对实际蚂蚁的一种简单模拟, 它融进了人类的智能L人工蚂蚁有一定的记忆; 人工蚂蚁不完全是瞎的; 人工蚂蚁生活的时空是离散的L
③ 它是一种分布式的优化方法, 不仅适合目前的串行计算机, 而且适合未来的并行计算机L
④ 它是一种全局优化的方法, 不仅可用于求解单目标优化问题, 而且可用于求解多目标优化问题L
⑤ 它是一种启发式算法, 计算复杂性为o (Nc*n2*m) , 其中Nc 是迭代次数, m 是蚂蚁数目, n 是目的节点数目L
蚁群发现最短路径的原理和机制[1]
下面用图 1解释蚁群发现最短路径的原理和机制。
如图 1(a)所示，在蚁巢和食物源之间有两条道路 Nest-A-B-D-Food 和Nest-A-C-D-Food，其长度分别为 4 和 6。单位时间内蚂蚁可移动一个单位长度的距离。开始时所有路径上都没有外激素。
如图 1(b)，在 t=0 时刻，20 只蚂蚁从蚁巢出发移动到 A。由于路径上没有外激素，它们以相同概率选择左侧或右侧道路，因此平均有 10 只蚂蚁走左侧，另外 10 只走右侧。
如图 1(c)，在 t=4 时刻，第一组先到达食物源的蚂蚁将折回。
如图 1(d)，在 t=5 时刻，两组蚂蚁将在 D 点相遇。此时 BD 上的外激素数量与 CD 上的相同，因此返回的 10 只蚂蚁中有 5 只选择 BD 而另 5 只选择 CD。
如图 1(e)，在 t=8 时刻，前 5 个蚂蚁将返回巢穴，而在 AC、CD 和 AB 上各有 5 个蚂蚁。
如图 1(f)，在 t=9 时刻，前 5 个蚂蚁又回到 A 并且再次面对往左还是往右的选择。这时，AB 上的轨迹数是 20 而 AC 上是 15，因此将有较为多数的蚂蚁选择往右，从而增强了 AB 上外激素的量。随着该过程的继续，两条道路上外激素数量的差距将越来越大，直至绝大多数蚂蚁都选择了最短的路径。正是由于一条道路要比另一条道路短，因此，在相同的时间间隔内，短的路线会有更多的机会被选择。

根据仿生学家的研究结果,蚂蚁凭借路径寻优的能力能够找到蚁巢与食物之间的最短路径,其原理在于:蚂蚁在所经过的路径上留下一种挥发性分泌物(pheromone,以下称为信息素),信息素随着时间的推移会逐渐挥发消失.蚂蚁在觅食过程中能够感知这种物质的存在及其强度,并以此来指导自己的运动方向,倾向于朝着这种物质强度高的方向移动,即选择该路径的概率与当时这条路径上该物质的强度成正比.信息素强度越高的路径,选择它的蚂蚁就越多,则在该路径上留下的信息素的强度就更大,而强度大的信息素又吸引更多的蚂蚁,从而形成一种正反馈.通过这种正反馈,蚂蚁最终可以发现最佳路径,导致大部分的蚂蚁都会走此路径.
以求解ｎ个城市的ＴＳＰ旅行商问题为例说明ＡＣＡ模型.
设蚁群中蚂蚁的数量为ｍ,ｄij (ｉ,ｊ=1,2,…,ｎ)表示城市ｉ和城市ｊ之间的距离,ｂi(ｔ)表示ｔ时刻位于城市ｉ的蚂蚁的个数,则有 表示ｔ时刻在城市ｉ,ｊ连线上残留的信息量.初始时刻,各条路径上信息量相等,设τｉｊ(0)=Ｃ(Ｃ为常数).蚂蚁ｋ(ｋ=1,2,…,ｍ)在运动过程中,根据各条路径上的信息量决定转移方向.
表示在ｔ时刻蚂蚁ｋ由城市ｉ转移到城市ｊ的概率.
（1）
残留信息的重要程度;β——启发信息的重要程度;tabuk——记录蚂蚁ｋ当前所走过的城市,称为记忆列表,ｋ=1,2,…,ｍ,集合tabuk随着进化过程作动态调整.经过ｎ个时刻,所有蚂蚁都完成了一次遍历.此时,计算每一只蚂蚁所走过的路径Ｌｋ,并保存最短路径Ｌｍｉｎ=ｍｉｎ{Ｌｋ︱ｋ=1,2,…,ｍ}.在蚂蚁完成一次循环以后,各路径上的信息量进行如下调整
τｉｊ(ｔ+1)=(1-ρ)τｉｊ(ｔ)+Δτｉｊ (2)
式中ρ∈(0,1),表示信息素τｉｊ(ｔ)随时间的推移而衰减的程度.所以1-ρ为信息素残留因子，开始时Δτｉｊ（0）=0，
信息素增量Δτｉｊ可表示 (3)
式中Δτｋｉｊ为蚂蚁ｋ在本次循环中在城市ｉ和ｊ之间留下的信息量,它的计算公式根据具体问题而定.Dorigo曾给出Δτｋｉｊ3种不同的模型,分别称为Ant-Cycle模型、Ant-Quantity模型、Ant-Density模型，它们的区别就在于信息素的更新机制,即其差别在于Δτｋｉｊ
在Ant-Cycle模型中：
（4） 式中，Q表示信息素强度，它在一定程度上影响算法的收敛速度；Lk表示第K只蚂蚁在本次循环中所奏路径的总长度。
在Ant-Quantity模型中：
（5） 式中，Q表示信息素强度，它在一定程度上影响算法的收敛速度；dij表示第K只蚂蚁在t和t+1之间经过的（ i, j ）
在Ant-Density模型中：
（6） 区别:式(5)式(6)中利用的是局部信息，即蚂蚁完成一步后更新路径上的信息素；而式(4)中利用的是整体信息，即蚂蚁完成一个循环后所有路径上的信息素。经过大量试验总结研究，采用式(4)性能较好，所以 Ant-Cycle模型是最优的。
以上说明了信息素残留因子1-ρ、信息启发式因子α、期望启发式因子β、信息素强度Q、蚂蚁数目M等都是非常重要的参数，其选区方式和选区原则直接影响到蚁群算法的全局收敛性和求解效率。我们学习到这种“三步走”[2]选择蚁群算法最优组合参数的有效方法：
（1） 确定蚂蚁数目M，根据 城市规模 / 蚂蚁数目 ≈1.5的选择策略来确定蚂蚁的总数目。
（2） 参数粗调，即调整数值范围较大的信息启发式因子α、期望启发式因子β、信息素强度Q等参数，已得到较理想的解。
（3） 参数微调，即调整数值范围较小的信息素残留因子1-ρ。
2 目前蚁群算法的应用
虽然对蚁群算法的研究时间不长, 但是初步研究已显示出它在求解复杂优化问题方***有很大的优势, 特别是1998 年在比利时布鲁塞尔专门召开了第一届蚂蚁优化国际研讨会后, 现在每两年召开一次这样的蚂蚁优化国际研讨会。这标志着蚁群算法的研究已经得到了国际上的广泛支持，使得这种新兴的智能进化仿生算法展现出了勃勃生机[3]。
以蚁群算法为代表的群体智能已成为当今分布式人工智能研究的一个热点,许多源于蜂群和蚁群模型设计的算法已越来越多地被用于企业的运转模式的研究。美国五角大楼正在资助关于群体智能系统的研究工作--群体战略(SWARM STRATEGY),它的一个实战用途是通过运用成群的空中无人驾驶飞行器和地面车辆来转移敌人的注意力,让自己的军队在敌人后方不被察觉地安全行进。英国电信公司和美国世界通信公司以电子蚂蚁为基础,对新的电信网络管理方法进行了试验。群体智能还被应用于工厂生产计划的制定和运输部门的后勤管理。美国太平洋西南航空公司采用了一种直接源于蚂蚁行为研究成果的运输管理软件,结果每年至少节约了1000万美元费用开支。英国联合利华公司已率先利用群体智能技术改善其一家牙膏厂的运转状况。美国通用汽车公司,法国液气公司,荷兰公路交通部和美国一些移民事务机构也都采用这种技术来改善其运转的机能。又如美国MCIWorld.com公司一直研究人工蚂蚁，并用于管理公司的电话网，对用户记账收费等工作。另外，还设计“人工蚂蚁”打算用于因特网的路由管理。鉴于群体智能广阔的应用前景,美国和欧洲联盟均于近几年开始出资资助基于群体智能模拟的相关研究项目,
关在一些院校开设群体智能的相关课程.牛津大学出版社1999年版的E.Bonabeau和M.Dorigo等人编写的专著《群体智能:从自然到人工系统》(Swarm Intelligence:From Natural to Artificial System),以及2001年出版的J.Kennedy和R.Eberhart编著的《群体智能》(Swarm Intelligence)进一步扩大了群体智能的影响.IEEE进化计算会刊也于2002年8月出版了蚁群优化算特刊。国内也有研究者用蚂蚁算法求解全国144个城市的最短回路问题，求得的解同其它方法求到得解一样精确，这说明蚂蚁算法不但是求解组合优化问题的可行方法，而且是一种很有竞争力的算法。国家自然科学基金"十五"期间学科交叉类优先资助领域中的认知科学及其信息处理的研究内容中也明确列出了群体智能领域的进化,自适应与现场认知主题[4]。而且从1999年开始，几乎每年都会有几项相关项目获得资助。蚁群算法是一种新型的模拟进化算法，其在数据挖掘中的应用正逐步引起人们的关注。目前，人工蚁群在知识发现的过程中主要用于发掘聚类模型和分类模型。
2．1蚁群算法在数据挖掘中的应用
聚类是将一组对象分成若干个群体，每个群体构成一个簇，使得簇内的对象尽可能具有最大的相似性，不同簇之间的对象尽可能有最大的相异性。目前，聚类方法主要有K均值法，模糊聚类、神经网络聚类、基于遗传算法的聚类、小波变换聚类以及将这些算法有效结合而形成的改进方法。随着蚁群算法研究的兴起，人们发现在某些方面采用蚁群模型进行聚类更加接近实际的聚类问题。将蚁群算法用于聚类分析,灵感源于蚂蚁堆积他们的尸体和分类他们的幼体。基于蚁群算法的聚类方法从原理上可分为两种：一种是基于蚁堆形成原理来实现数据聚类，另一种是运用蚂蚁觅食的原理，利用信息来实现聚类分析。
而数据是数据挖掘的另一个重要主题，它是在数据库对象集合中寻找属性，并根据分类模式将其划分为不同类别的过程。分类过程利用历史数据记录自动推导出对给定数据的分类树。分类器构造方法有统计学方法、机器学习法、神经网络、决策树等。从知识发现的观点来看，分类规则的表达方式形如if<条件>then<类>规则前件（if
部分）包含一组条件集合，一般由逻辑连接符连接；规则结论（then部分）定义了样本的预测类，这些样本的预测属性满足规则前件所定义的所有条件[5]。将蚁群算法引入分类规则的发现，是利用蚁群觅食原理在数据库中进行搜索，对随机产生的一组规则进行选择优化，直到数据库能被该组规则覆盖，从而挖掘出隐含在数据库中的规则，建立最优的分类模型。蚁群算法搜索的初始条件为发现规则的集合为空，且训练集包含所有的训练样本。蚂蚁搜索一次要完成规则生成、规则剪枝、信息素更新三个任务。一次搜索生成一条规则，并且将这条规则加入发现规则集合，同时将该条规则所覆盖的训练样本从训练集中删除。如果未覆盖训练样本的数目大于用户定义的阈值，即最大未覆盖样本数，就反复执行上述过程，最终算法将得到一组最优分类规则集合[5]。
最早在这一领域开展工作的是Deneubourg 等[6],他们根据数据对象与其周围对象的相似性,让蚂蚁随机地移动、拾起或放下数据对象,以达到聚类数据的目的,这个基本模型已成功地应用于机器人领域。Lumer 等首先改进此算法,提出了LF算法。Wu 等、Ramos等、Yang等[7]从不同角度对LF算法进行了改进,在用蚁群算法进行聚类分析方面取得了一定成效。近几年，学者在这方面的研究从来没有间断过，也取得了一定的研究成果。
2．2 结论
不过，将蚁群算法运用于数据发掘还存在一些问题，需要进一步研究：
（1）如何将现实的挖掘任务转换成蚁群求解的问题空间，并用适当的方式表达。如何定义“人工蚂蚁”以及蚂蚁间的非直接通信方式（如路径上的信息素、对象的分布状态等）的选择。
（2）如何建立正反馈机制，定义启发函数，递增地进行问题求解，并且使得到的解与问题定义中现实世界的情况相对应。
（3）基于蚁群的算法要初始化大量的参数，这些参数的选择会对算法的性能产生较大的影响，但其选取的方法和原则目前尚无理论上的依据，只能通过多次实验调优，因此参数的最佳设置原则还有待进一步研究。
（4）蚁群算法的搜索时间较长，如何将蚁群算法与遗传算法、免疫算法等优化算法相结合，改善和提高算法性能，以适应海量数据库的知识发现。
所以如何在数据挖掘中运用蚁群算法快速、高效地获得高质量的知识越来越受到人们的关注，逐渐成为近期的研究热点[5]。

以下是解放军信息工程大学一个老师编的matlab程序，请尊重原作者劳动。
我经过修改增加了注释，已经运行过，无误，

function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)
%%-------------------------------------------------------------------------
%% 主要符号说明
%% C n个城市的坐标，n×2的矩阵
%% NC_max 最大迭代次数
%% m 蚂蚁个数
%% Alpha 表征信息素重要程度的参数
%% Beta 表征启发式因子重要程度的参数
%% Rho 信息素蒸发系数
%% Q 信息素增加强度系数
%% R_best 各代最佳路线
%% L_best 各代最佳路线的长度
%%=========================================================================
%%第一步：变量初始化
n=size(C,1);%n表示问题的规模（城市个数）
D=zeros(n,n);%D表示完全图的赋权邻接矩阵
for i=1:n
for j=1:n
if i~=j
D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
else
D(i,j)=eps; %i=j时不计算，应该为0，但后面的启发因子要取倒数，用eps（浮点相对精度）表示
end
D(j,i)=D(i,j); %对称矩阵
end
end
Eta=1./D; %Eta为启发因子，这里设为距离的倒数
Tau=ones(n,n); %Tau为信息素矩阵
Tabu=zeros(m,n); %存储并记录路径的生成
NC=1; %迭代计数器，记录迭代次数
R_best=zeros(NC_max,n); %各代最佳路线
L_best=inf.*ones(NC_max,1); %各代最佳路线的长度
L_ave=zeros(NC_max,1); %各代路线的平均长度
while NC<=NC_max %停止条件之一：达到最大迭代次数，停止
%%第二步：将m只蚂蚁放到n个城市上
Randpos=[]; %随即存取
for i=1:(ceil(m/n))
Randpos=[Randpos,randperm(n)];
end
Tabu(:,1)=(Randpos(1,1:m))'; %此句不太理解？
%%第三步：m只蚂蚁按概率函数选择下一座城市，完成各自的周游
for j=2:n %所在城市不计算
for i=1:m
visited=Tabu(i,1:(j-1)); %记录已访问的城市，避免重复访问
J=zeros(1,(n-j+1)); %待访问的城市
P=J; %待访问城市的选择概率分布
Jc=1;
for k=1:n
if length(find(visited==k))==0 %开始时置0
J(Jc)=k;
Jc=Jc+1; %访问的城市个数自加1
end
end
%下面计算待选城市的概率分布
for k=1:length(J)
P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta);
end
P=P/(sum(P));
%按概率原则选取下一个城市
Pcum=cumsum(P); %cumsum，元素累加即求和
Select=find(Pcum>=rand); %若计算的概率大于原来的就选择这条路线
to_visit=J(Select(1));
Tabu(i,j)=to_visit;
end
end
if NC>=2
Tabu(1,:)=R_best(NC-1,:);
end
%%第四步：记录本次迭代最佳路线
L=zeros(m,1); %开始距离为0，m*1的列向量
for i=1:m
R=Tabu(i,:);
for j=1:(n-1)
L(i)=L(i)+D(R(j),R(j+1)); %原距离加上第j个城市到第j+1个城市的距离
end
L(i)=L(i)+D(R(1),R(n)); %一轮下来后走过的距离
end
L_best(NC)=min(L); %最佳距离取最小
pos=find(L==L_best(NC));
R_best(NC,:)=Tabu(pos(1),:); %此轮迭代后的最佳路线
L_ave(NC)=mean(L); %此轮迭代后的平均距离
NC=NC+1 %迭代继续
%%第五步：更新信息素
Delta_Tau=zeros(n,n); %开始时信息素为n*n的0矩阵
for i=1:m
for j=1:(n-1)
Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);
%此次循环在路径（i，j）上的信息素增量
end
Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);
%此次循环在整个路径上的信息素增量
end
Tau=(1-Rho).*Tau+Delta_Tau; %考虑信息素挥发，更新后的信息素
%%第六步：禁忌表清零
Tabu=zeros(m,n); %%直到最大迭代次数
end
%%第七步：输出结果
Pos=find(L_best==min(L_best)); %找到最佳路径（非0为真）
Shortest_Route=R_best(Pos(1),:) %最大迭代次数后最佳路径
Shortest_Length=L_best(Pos(1)) %最大迭代次数后最短距离
subplot(1,2,1) %绘制第一个子图形
DrawRoute(C,Shortest_Route) %画路线图的子函数
subplot(1,2,2) %绘制第二个子图形
plot(L_best)
hold on %保持图形
plot(L_ave,'r')
title('平均距离和最短距离') %标题
function DrawRoute(C,R)
%%=========================================================================
%% DrawRoute.m
%% 画路线图的子函数
%%-------------------------------------------------------------------------
%% C Coordinate 节点坐标，由一个N×2的矩阵存储
%% R Route 路线
%%=========================================================================
N=length(R);
scatter(C(:,1),C(:,2));
hold on
plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)],'g')
hold on
for ii=2:N
plot([C(R(ii-1),1),C(R(ii),1)],[C(R(ii-1),2),C(R(ii),2)],'g')
hold on
end
title('旅行商问题优化结果 ')