基于蚁群算法的TSP问题
基于蚁群算法的TSP问题求解
1引言
1.1 问题描述
设计求解以下两个TSP问题的蚁群优化（ACO）算法。其中城市的坐标见附件（kroA100.tsp和kroB100.tsp）。
1.2 理论基础
1.2.1 蚁群算法简介
蚁群算法是由意大利学者M.Dorigo等人于20世纪90年代初提出的一种新的模拟进化算法，其真实地模拟了自然界蚂蚁群体的觅食行为。M.Dorigo等人将其用于解决旅行商问题（traveling salesman problem, TSP），并取得了较好的实验结果。
近年来，许多专家学者致力于蚁群算法的研究，并将其应用于交通、通信、化工、电力等领域，成功解决了许多组合优化问题，如调度问题（job–shop scheduling problem）、指派问题（quadratic assignment problem）、旅行商问题（traveling salesman problem）等。
1.2.2 蚁群算法基本思想
生物学家研究发现，自然界中的蚂蚁觅食是一种群体性行为，并非单只蚂蚁自行寻找食物源。蚂蚁在寻找食物源时，会在其经过的路径上释放一种信息素，并能够感知其它蚂蚁释放的信息素。信息素浓度的大小表征路径的远近，信息素浓度越高，表示对应的路径距离越短。通常，蚂蚁会以较大的概率优先选择信息素浓度较高的路径，并释放一定量的信息素，以增强该条路径上的信息素浓度，这样会形成一个正反馈。最终，蚂蚁能够找到一条从巢穴到食物源的最佳路径，即是最短距离。值得一提的是，生物学家同时发现，路径上的信息素浓度会随着时间的推进而逐渐衰减。
将蚁群算法应用于解决优化问题的基本思路为：用蚂蚁的行走路径表示待优化问题的可行解，整个蚂蚁群体的所有路径构成待优化问题的解空间。较短的路径上蚂蚁释放的信息素量较多，随着时间的推进，较短的路径上积累的信息素浓度逐渐增高，选择该路径的蚂蚁个数也愈来愈多。最终，整个蚂蚁会在正反馈的作用下集中到最佳的路径上，此时对应的便是待优化问题的最优解。
1.2.3 蚁群算法解决TSP 问题的基本原理
为不失一般性，设整个蚂蚁群体中蚂蚁的数量为m ，城市的数量为n ，城市
i 与城市j 之间的距离为ij d (,1,2,...,)i j m =，t 时刻城市i 与城市j 连接路径上的信息素浓度为()ij t τ。初始时刻，各个城市连接路径上的信息素浓度相同，不妨设为0(0)ij ττ=。
蚂蚁(1,2...,)k k m =根据各个城市间连接路径上的信息素浓度决定其下一个
访问城市，设()k
ij
p t 表示t 时刻蚂蚁k 从城市j 的概率，其计算公式为：
[][]()(),()()
0,k ij ij k k
ij is is s allow k t t s allow p t t s allow αβ
αβτητη∈???????????∈?=??????
∑ （1-1）
其中()ij t η为启发函数，()1/ij ij t d η=，表示蚂蚁从城市i 转移到城市j 的期望
程度；(1,2,...,)k allow k m =为蚂蚁k 待访问城市的集合，开始时，
k allow 中有(1)n -个元素，即包括除了蚂蚁k 出发城市的其它所有城市，随着时间的推进，k allow 中的元素不断减少，直至为空，即表示所有的城市均访问完毕；α为信息素重要程度因子，其值越大，表示信息素的浓度在转移中起的作用越大；β为启发函数重要程度因子，其值越大，表示启发函数在转移中的作用越大，即蚂蚁会以较大的概率转移到距离短的城市。
如前文所述，在蚂蚁释放信息素的同时，各个城市间连接路径上的信息素逐渐消失，设参数(01)ρρ<<表示信息素的挥发程度。因此，当所有蚂蚁完成一次循环后，各个城市间连接路径上的信息素浓度需进行实时更新，即：
1(1)(1)(),01ij ij ij n k
ij ij k t t τρττρττ=+=-+???<
?
∑ （1-2）
其中，k ij τ?表示第k 只蚂蚁在城市i 与城市j 连接路径上释放的信息素浓度；
ij τ?表示所有蚂蚁在城市i 与城市j 连接路径上释放的信息素浓度之和。
针对蚂蚁释放信息素问题，M.Dorigo 等人曾给出了三种不同的模型，分别称之为Ant cycle system ，Ant quantity system 和Ant density system ，其计算公式分别如下。
（1） Ant cycle system 模型
Ant cycle system 模型中，k ij τ?的计算公式为： /,0,k k
ij
Q L k i j τ??=??
第只蚂蚁从城市访问城市其他 （1-3）
其中，Q 为常数，表示蚂蚁循环一次所释放的信息素总量；k L 为第k 只蚂蚁经过路径的长度。
（2） Ant quantity system 模型
Ant quantity system 模型中，k ij τ?的计算公式为：
/,=0,ij k ij Q d k i j
τ????
第只蚂蚁从城市访问城市其他 （1-4）
（3） Ant density system 模型
Ant density system 模型中，k ij τ?的计算公式为：
,=0,k ij Q k i j
τ????
第只蚂蚁从城市访问城市其他 （1-5）
上述3种模型中，Ant cycle system 模型利用蚂蚁经过路径的整体信息（经过路径的总长）计算释放的信息素浓度；Ant quantity system 模型则利用蚂蚁经过路径的局部信息（经过各个城市间的距离）计算释放的信息素浓度；而Ant density system 模型则更为简单地将信息素释放的浓度取为恒值，并没有考虑不同蚂蚁经过路径长短的影响。因此，一般选用Ant cycle system 模型计算释放的信息素浓度，即蚂蚁经过的路径越短，释放的信息素浓度越高。
1.2.4 蚁群算法的特点
基于蚁群算法的基本思想及解决TSP 问题的基本原理，不难发现，与其它
优化算法相比，蚁群算法具有以下几个特点：
（1）采用正反馈机制，使得搜索过程不断收敛，最终逼近最优解。
（2）每个个体可以通过释放信息素来改变周围的环境，且每个个体能够感知周围环境的实时变化，个体间通过环境进行间接地通讯。
（3）搜索过程采用分布式计算方式，多个个体同时进行并计算，大大提高了算法的计算能力和运行效率。
（4）启发式的概率搜索方式不容易陷入局部最优，易于寻找到全局最优解。
1 算法设计
蚁群算法应用于解决TSP 问题一般需要以下几个步骤，如图1所示。
图1 蚁群算法解决TSP 问题的基本步骤
（1） 初始化参数
在计算之初，需要对相关的参数进行初始化，如蚁群规模（蚂蚁数量）m 、信息素重要程度因子α、启发函数重要程度因子β、信息素挥发因子ρ、信息素释放总量Q 、最大迭代次数iter_max 、迭代次数初值iter=1。
（2） 构建解空间
将各个蚂蚁随机地置于不同出发点，对每个蚂蚁(1,2...,)k k m =，按式（1-1）计算其下一个待访问的城市，直到所有蚂蚁访问完所有的城市。 （3） 更新信息素
计算各个蚂蚁经过的路径长度(1,2...,)k L k m =，记录当前迭代次数中的最优
解（最短路径）。同时，根据式（1-2）和式（1-3）对各个城市连接路径上的信息素浓度进行更新。
（4） 判断是否终止
若iter 2 算例设计 依据蚁群算法解决TSP 问题的基本原理及步骤，实现该TSP 问题求解大体上可以分为以下几个步骤，如图2所示。 图2 蚁群算法求解TSP 问题一般步骤 （1） 计算城市间相互距离 根据城市的坐标位置，计算两两城市间的相互距离，从而得到对称的距离矩阵（维数为100的方阵）。需要说明的是，计算出的矩阵对角线上的元素为0，然而如前文所述，由于启发函数()1/ij ij t d η=，因此，为了保证分母不为零，将对角线上的元素修正为一个非常小的正数。 （2） 初始化参数 如前文所述，在计算之前，需要对相关的参数进行初始化。 （3） 迭代寻找最佳路径 首先构建解空间，即各个蚂蚁根据转移概率公式（1-1）访问所有的城市。然后计算各个蚂蚁经过路径的长度，并在每次迭代后根据式（1-2）和式（1-3）实时更新各个城市连接路径上的信息素浓度。经过循环迭代，记录下最优的路径 及长度。 （4）结果分析 找到最优路径后，可以将之与其它方法得出的结果进行比较，从而对蚁群算法的性能进行评价。同时，也可以研究不同取值的参数对优化结果的影响，从而找到一组最佳或者较佳的参数组合。 3仿真实验设计与结果分析 4.1 仿真实验结果 编写MATLAB程序，分别载入两个附件中100个城市的坐标，得到基于蚁群算法的TSP问题的最优解如下： （1）KroA100 最短距离：22756.0988 最短路径：91 98 23 45 32 11 17 15 59 74 21 72 10 84 36 99 38 24 18 79 53 88 16 22 94 70 66 65 4 26 19 75 97 56 80 31 89 42 8 92 1 63 6 49 90 47 93 28 67 58 61 25 81 69 64 40 54 2 44 50 73 68 85 39 30 96 78 52 5 37 13 76 33 95 82 48 100 71 41 14 3 43 46 29 34 83 55 12 27 86 35 20 57 7 9 87 51 77 62 60 91 （2）KorB100 最短距离：23537.4394 最短路径：96 92 19 44 41 17 45 36 18 24 77 16 2 13 78 63 31 48 51 64 14 42 82 33 15 6 4 83 87 60 74 66 89 43 52 54 58 84 7 57 94 61 35 71 12 98 32 59 76 99 8 29 97 91 28 3 11 93 85 73 53 70 39 40 67 5 62 26 100 56 69 75 30 80 20 38 86 49 68 10 21 90 1 95 46 25 9 27 34 72 37 65 79 81 47 88 22 23 55 50 96 运行结果对应的路径如图3和图4所示，各代的最短距离与平均距离如图5和图6所示。从图中不难发现，随着迭代次数的增加，最短距离与平均距离均呈现下降趋势。当迭代次数大于某一数值后最短距离已不再变化，说明已经找到了最优解。 图3 KroA优化路径 图4 KroB优化路径 图5 KroA各代最短距离与平均距离 图6 KroB各代最短距离与平均距离 4.2 参数对优化结果的影响 （1）蚂蚁数量m的影响 为研究蚂蚁数m对计算结果的影响，比较10组不同蚂蚁数的计算结果，每组计算10次取平均值。同时，由于计算次数多达100次，为提高程序运行速度，暂将最大迭代次数减少为20次，仅对m变化引起结果变化的趋势作以研究（下同）。取α=1，β=5，ρ=0.1，Q=1，研究结果如图7和图8所示。 图7 蚂蚁数对计算结果的影响 图8 蚂蚁数对平均最短距离的影响 可见，随着蚂蚁数的增加，平均最短距离呈减小的趋势，图中蚂蚁数为45时，得到最小的平均最短距离为27056。 为研究信息素重要程度因子α对计算结果的影响，比较10组不同α值的计算结果，每组计算10次取平均值。取m=45，β=5，ρ=0.1，Q=1，研究结果如图9和图10所示。 图9 信息素重要程度因子对计算结果的影响 图10 信息素重要程度因子对平均最短距离的影响 可见，随α增大，平均最短距离大体呈下降趋势，当α=10时，得到最小的平均最短距离为27004。 为研究启发函数重要程度因子β对计算结果的影响，比较10组不同β值的计算结果，每组计算10次取平均值。取m=45，α=10，ρ=0.1，Q=1，研究结果如图11和图12所示。 图11 启发函数重要程度因子对计算结果的影响 图12 启发函数重要程度因子对平均最短距离的影响 可见，随β增大，平均最短距离呈下降趋势，而当β=10时，得到最小的平均最短距离为24438。 （4）信息素挥发因子ρ的影响 为研究信息素挥发因子ρ对计算结果的影响，比较10组不同ρ值的计算结果，每组计算10次取平均值。取m=45，α=10，β=10，Q=1，研究结果如图13和图14所示。 图13 信息素挥发因子对计算结果的影响 图14 信息素挥发因子对平均最短距离的影响 可见，随ρ增大，平均最短距离呈先下降后上升趋势，而当ρ=0.2时，得到最小的平均最短距离为24068。 （5）信息素释放量Q的影响 为研究信息素释放总量Q对计算结果的影响，比较10组不同Q值的计算结果，每组计算10次取平均值。取m=45，α=10，β=10，ρ=0.2，研究结果如图15 和图16所示。 图15 信息素释放量对计算结果的影响 图16 信息素释放量对平均最短距离的影响 可见，随Q增大，平均最短距离大致呈上升趋势，而当Q=10时，得到最小的平均最短距离为24179。 综合分析以上参数的影响，可以看出，随着研究的进行，通过控制变量法，得到的平均最短距离有下降趋势，但最后一次的24179却比倒数第二次的24068要大。我们分析这里可能有两个原因，一是因为我们为尽快运行程序，降低了迭代次数，算法寻优过程的偶然性增大；二是由于蚁群算法本身的偶然性，即使将 迭代次数调整回200次，运行结果有时也会偏大，但增大迭代次数，得到最优解的概率会越大。 4总结与展望 本文详细叙述了应用蚁群优化算法解决TSP问题的一般过程，并将其应用于题中所给出的具体问题。算法使用MATLAB编写程序实现，给出了程序运行的仿真结果并对其进行了分析。考虑程序运行的结果与效率，可以看出蚁群优化算法可以很好的解决TSP等组合优化问题。但是算法中的参数一般由人为确定，导致方法的优化性能往往与人的经验密切相关。因此本文接下来分析了蚁群算法中参数的选取对结果的影响，通过仿真实验得到了求解该TSP问题的最佳参数。 蚁群算法以其分布并行式计算、启发式搜索方式等特点，在各个领域中取得了广泛的应用，解决了诸多组合优化方面的难题。针对其可能陷入局部最优的缺点，不少专家和学者提出了许多改进的方法，并将其它算法（如遗传算法、粒子群算法等）与蚁群算法相结合，对蚂蚁个数、信息素重要程度因子、启发函数重要程度因子、信息素挥发因子、信息素释放总量等参数进行优化选择，取得了不错的效果。这些都可以作为我们以后的研究课题。 附录 MATLAB源程序代码clear all close all clc load Data2.txt; City=Data2(:,2:end); n=size(City,1); D=zeros(n,n); %%Distance for i=1:n for j=1:n if i~=j D(i,j)=sqrt(sum((City(i,:)-City(j,:)).^2)); else D(i,j)=1e-4; end end end %%Initialize m=100;%蚂蚁数量 alpha=1;%信息素重要程度因子 beta=5;%启发函数重要程度因子 rho=0.1;信息素挥发因子 Q=1;%常系数 Eta=1./D;%启发函数 Tau=ones(n,n);%信息素矩阵 Table=zeros(m,n);%路径记录 iter=1;%迭代次数 iter_max=200; Route_best=zeros(iter_max,n); Length_best=zeros(iter_max,1); Length_ave=zeros(iter_max,1); %%iteration while iter<=iter_max %随机产生各蚂蚁起点 start=zeros(m,1); for i=1:m temp=randperm(n); start(i)=temp(i); end Table(:,1)=start; %构建解空间 citys_index=1:n; for j=2:n tabu=Table(i,1:(j-1)); allow_index=~ismember(citys_index,tabu); allow=citys_index(allow_index); P=allow; %城市间转移概率 for k=1:length(allow) P(k)=Tau(tabu(end),allow(k))^alpha*Eta(tabu(end),allow(k))^beta; end P=P/sum(P); Pc=cumsum(P); target_index=find(Pc>=rand); target=allow(target_index(1)); Table(i,j)=target; end end %各蚂蚁路径 Length=zeros(m,1); for i=1:m Route=Table(i,:); for j=1:(n-1) Length(i)=Length(i)+D(Route(j),Route(j+1)); end Length(i)=Length(i)+D(Route(n),Route(1)); end %最短路径距离和平均距离 if iter==1 [min_Length,min_index]=min(Length); Length_best(iter)=min_Length; Length_ave(iter)=mean(Length); Route_best(iter,:)=Table(min_index,:); else [min_Length,min_index]=min(Length); Length_best(iter)=min(Length_best(iter-1),min_Length); Length_ave(iter)=mean(Length); if Length_best(iter)==min_Length Route_best(iter,:)=Table(min_index,:); else Route_best(iter,:)=Route_best(iter-1,:); end end %更新信息素 Delta_Tau=zeros(n,n); for j=1:n-1 Delta_Tau(Table(i,j),Table(i,j+1))=Delta_Tau(Table(i,j),Table(i,j+1))+Q/Length(i); end Delta_Tau(Table(i,n),Table(i,1))=Delta_Tau(Table(i,n),Table(i,1))+Q/Length(i); end Tau=(1-rho)*Tau+Delta_Tau; iter=iter+1; Table=zeros(m,n); end %result [Shortest_Length,index]=min(Length_best); Shortest_Route=Route_best(index,:); disp(['The Shortest Length￡o',num2str(Shortest_Length)]); disp(['The Shortest Route￡o',num2str([Shortest_Route Shortest_Route(1)])]); figure(1) plot([City(Shortest_Route,1);City(Shortest_Route(1),1)],[City(Shortest_Route,2);City(Shortest_Route( 1),2)],'o-'); grid on; for i=1:size(City,1) text(City(i,1),City(i,2),[' ',num2str(i)]); end text(City(Shortest_Route(1),1),City(Shortest_Route(1),2),' starting'); text(City(Shortest_Route(end),1),City(Shortest_Route(end),2),' end'); xlabel('The X-coordinate of City'); ylabel('The Y-coordinate of City'); title(['ACA optimal result:',num2str(Shortest_Length)]); figure(2) plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r') legend('The Shortest Length','The Average Length') xlabel('Iterations') ylabel('Length') title('The Shortest Length and Average Length of each iteration')