蚁群算法原理及在TSP中的应用(附程序)

蚁群算法原理及在TSP 中的应用1 蚁群算法（ACA ）原理1.1 基本蚁群算法的数学模型以求解平面上一个n 阶旅行商问题(Traveling Salesman Problem ，TSP)为例来说明蚁群算法ACA （Ant Colony Algorithm ）的基本原理。

对于其他问题，可以对此模型稍作修改便可应用。

TSP 问题就是给定一组城市，求一条遍历所有城市的最短回路问题。

设()i b t 表示t 时刻位于元素i 的蚂蚁数目，()ij t τ为t 时刻路径(,)i j 上的信息量，n 表示TSP 规模，m 为蚁群的总数目，则1()ni i m b t ==∑；{(),}ij i i t c c C τΓ=⊂是t 时刻集合C 中元素(城市)两两连接ij t 上残留信息量的集合。

在初始时刻各条路径上信息量相等，并设 (0)ij const τ=，基本蚁群算法的寻优是通过有向图(,,)g C L =Γ实现的。

蚂蚁(1,2,...,)k k m =在运动过程中，根据各条路径上的信息量决定其转移方向。

这里用禁忌表(1,2,...,)k tabu k m =来记录蚂蚁k 当前所走过的城市，集合随着k tabu 进化过程作动态调整。

在搜索过程中，蚂蚁根据各条路径上的信息量及路径的启发信息来计算状态转移概率。

()kij p t 表示在t 时刻蚂蚁k 由元素(城市)i 转移到元素(城市)j 的状态转移概率。

()*()()*()()0k ij ij k kij ij ij s allowed t t j allowed t t p t αβαβτητη⊂⎧⎡⎤⎡⎤⎣⎦⎣⎦⎪∈⎪⎡⎤⎡⎤=⎨⎣⎦⎣⎦⎪⎪⎩∑若否则(1)式中，{}k k allowed C tabuk =-表示蚂蚁k 下一步允许选择的城市；α为信息启发式因子，表示轨迹的相对重要性，反映了蚂蚁在运动过程中所积累的信息在蚂蚁运动时所起作用，其值越大，则该蚂蚁越倾向于选择其他蚂蚁经过的路径，蚂蚁之间协作性越强；β为期望启发式因子，表示能见度的相对重要性，反映了蚂蚁在运动过程中启发信息在蚂蚁选择路径中的重视程度，其值越大，则该状态转移概率越接近于贪心规则；()ij t η为启发函数，其表达式如下：1()ij ijt d η=(2)式中，ij d 表示相邻两个城市之间的距离。

对蚂蚁k 而言，ij d 越小，则()ij t η越大，()kij p t 也就越大。

显然，该启发函数表示蚂蚁从元素(城市) i 转移到元素(城市) j 的期望程度。

为了避免残留信息素过多引起残留信息淹没启发信息，在每只蚂蚁走完一步或者完成对所有n 个城市的遍历(也是一个循环结束)后，要对残留信息进行更新处理。

这种更新策略模仿了人类大脑记忆的特点，在新信息不断存入大脑的同时，存储在大脑中的旧信息随着时间的推移逐渐淡化，甚至忘记。

由此，t n +时刻在路径(,)i j 上的信息量可按如下规则进行调整：()(1)*()()ij ij ij t n t t τρττ+=-+∆ (3)1()()mkij ij k t t ττ=∆=∆∑ (4)式中ρ表示信息挥发系数，则1ρ-表示信息素残留因子，为了防止信息的无限积累，ρ的取值范围为：[0,1)ρ⊂；()ij t τ∆表示本次循环中路径(,)i j 上的信息素增量，初始时刻()0ij t τ∆=，()kij t τ∆表示第k 只蚂蚁在本次循环中留在路径(,)i j 上的信息量。

根据信息素更新策略的不同，Dorigo M 提出了三种不同的基本蚁群算法模型，分别称之为Ant-Cycle 模型、Ant-Quantity 模型及Ant-Density 模型，其差别在于()k ij t τ∆求法的不同。

在Ant-Cycle 模型中,k i j ()0kkij QL t τ⎧⎪∆=⎨⎪⎩若第只蚂蚁在本次循环中经过（，），否则(5)式中，Q 表示信息素强度，它在一定程度上影响算法的收敛速度；k L 表示k 只蚂蚁在本次循环中所走路径的总长度。

在Ant-Quantity 模型中,k i j ()0kijij Q d t τ⎧⎪∆=⎨⎪⎩若第只蚂蚁在t 和t+1之间经过（，），否则(6)在Ant-Density 模型中,k i j ()0kij Q t τ⎧∆=⎨⎩若第只蚂蚁在t 和t+1之间经过（，），否则(7)区别：式(6)和式(7)中利用的是局部信息，即蚂蚁完成一步后更新路径上的信息素；而式(5)中利用的是整体信息，即蚂蚁完成一个循环后更新所有路径上的信息素，在求解TSP 时性能较好，因此通常采用式(5)作为蚁群算法的基本模型。

1.2 基本蚁群算法的实现以下是解决TSP 问题的蚁群算法的基本流程描述，其中的参数设置来自于Dorigo 等人的试验。

基本蚁群算法的具体实现步骤如下：(1) 参数初始化。

令时间t=0和循环次数0c N =，设置最大循环次数max c N ，将m 蚂蚁置于n 个元素(城市)上，另有向图上每条边(,)i j 的初始化信息量()ij t const τ=，其中const 表示常数，且初始时刻(0)0ij τ∆=。

(2) 循环次数1c c N N ←+。

(3) 蚂蚁等禁忌表索引号k=1。

(4) 蚂蚁数目1k k ←+。

(5) 蚂蚁个体根据状态转移概率公式(1)计算的概率选择元素(城市)j 并前进，{}k j C tabu ∈-(6) 修改禁忌表指针，即选择好之后将蚂蚁移动到新的元素(城市)，并把该元素(城市)移动到该蚂蚁个体的禁忌表中。

(7) 若集合C 中元素(城市)未遍历完，即k<m ，则跳转到第(4)步，否则执行第(8)步。

(8) 根据公式(3)和(4)更新每条路径上的信息量。

(9) 若满足结束条件，即如果循环次数max c c N N ≥，则循环结束并输出程序计算结果，否则清空禁忌表并跳转到第(2)步。

2 程序实现（城市为中国各省省会城市）2.1数据列表TSP 问题中的城市选为中国各省省会城市，其实际地理位置的经纬度如表1所示。

表1 中国各省省会城市所在地理位置的经纬度2.2程序实现程序说明：程序是在GreenSim团队编写的程序的基础上进行修改得到，区别在于在程序中加入了更多的说明，便于MATLAB初学者以及对蚁群算法不是很熟悉的读者理解；对程序进行了略微的修改，加入了主程序和真实城市坐标数据。

其中主程序名称为ACAmain_city.m，其是可直接运行的程序。

运行结果如图1-图3所示，结果经过略微调整。

图2为将省会城市改成相应的省（或直辖市、特别行政区）后的图形。

结果均为一次运行所得，但由于算法的随机性，每次运行结果可能会不同。

另外，算法的初始条件对运行结果有很大影响，各种初始条件取值在附录中的程序已给出。

图1 ACA找到的最佳的中国环游路径（省会城市）图2 ACA找到的最佳的中国环游路径（省）图3 ACA迭代收敛曲线附录：（程序）注意：程序是直接从MATLAB中复制过来，拷贝到Word中可能会自动换行，再拷贝到MATLAB中运行时，应作相应调整。

函数文件应放到新的m文件中，并以相应的文件名命名，可运行的只有主程序。

A: 主程序%% 蚁群算法求解最短路径问题：ACAmain_city.mclc; close all; clear all; tic%% 中国各省省会城市23+5(直辖市4，特别行政区2)所在地理位置的经纬度（东经-北纬）C=[ 117.17 31.52; 119.18 26.05; 103.51 36.04; 113.14 23.08; 106.42 26.35; 110.20 20.02; 114.30 38.02; 113.40 34.46; 126.36 45.44; 114.17 30.35; 112.59 28.12; 125.19 43.54; 118.46 32.03; 115.55 28.40; 123.25 41.48; 101.48 36.38; 117.00 36.40; 112.3337.54; 108.57 34.17; 104.04 30.40; 102.42 25.04; 120.10 30.16; 121.30 25.03; 108.19 22.48; 111.41 40.48; 106.16 38.27; 90.08 29.39; 87.36 43.45; 121.29 31.14; 116.24 39.55; 117.12 39.02; 106.33 29.35; 115.07 21.33; 115.12 21.23];T=[ '合肥','福州','兰州','广州','贵阳','海口','石家','郑州','哈尔','武汉','长沙','长春','南京','南昌','沈阳','西宁','济南','太原','西安','成都','昆明','杭州','台北','南宁','呼和','银川','拉萨','乌鲁','上海','北京','天津','重庆','澳门','香港'];%% 参数初始化NC_max=150; m=34; Alpha=1; Beta=5; Rho=0.1; Q=100;%% 绘制找到的最优路径[R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha, Beta,Rho,Q); %函数调用figure(1); DrawCity(C,T,Shortest_Route); %绘制找到的最优路径toc %计算运行时间%% 绘制收敛曲线figure(2); iter=1:length(L_best);plot(iter,L_best,'-m*',iter,L_ave,':rp','LineWidth',2)xlabel('迭代次数'); legend('各代最佳路线的长度','各代路线的平均长度');grid on; tocB: 函数文件1%% 利用蚁群算法解决TSP问题的函数文件ACATSP.mfunction[R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)%% Ant Colony Algorithm for Traveling Salesman Problem%% 主要符号说明%% C n个城市的坐标，n×2的矩阵%% NC_max 最大迭代次数%% m 蚂蚁个数%% Alpha 表征信息素重要程度的参数%% Beta 表征启发式因子重要程度的参数%% Rho 信息素蒸发系数%% Q 信息素增加强度系数%% R_best 各代最佳路线%% L_best 各代最佳路线的长度%% NC_max=150; m=25; Alpha=1; Beta=5; Rho=0.1; Q=100;%%===================================================== ===========%% 第一步：变量初始化n=size(C,1); % n表示问题的规模（城市个数）D=zeros(n,n); % D表示完全图的赋权邻接矩阵for i=1:nfor j=1:nif i~=jD(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5; % 计算两城市之间的距离elseD(i,j)=eps; % eps = 2.2204e-016，i=j,则距离为0endD(j,i)=D(i,j); % 距离矩阵为对称均值（n*n的矩阵）endendEta=1./D; %Eta为启发因子，这里设为距离的倒数（n*n的矩阵）Tau=ones(n,n); %Tau为信息素矩阵，初始化全为1，Tabu=zeros(m,n); %存储并记录路径的生成tabu:踏步（停止，禁忌表）（m*n的矩阵）NC=1; %迭代计数器R_best=zeros(NC_max,n); %各代最佳路线(行数为最大迭代次数NC_max，列数为城市个数n)L_best=inf.*ones(NC_max,1); %各代最佳路线的长度（inf:无穷大）L_ave=zeros(NC_max,1); %各代路线的平均长度%%while NC<=NC_max %停止条件之一：达到最大迭代次数%% 第二步：将m只蚂蚁放到n个城市上Randpos=[];for i=1:(ceil(m/n)) %ceil表示向无穷方向取整Randpos=[Randpos,randperm(n)]; %randperm(n)：随机产生一个整数排列endTabu(:,1)=(Randpos(1,1:m))'; %每只蚂蚁（m只）都对应有一个位置，Tabu(:,1)为每只蚂蚁的走过的第一个城市%% 第三步：m只蚂蚁按概率函数选择下一座城市，完成各自的周游for j=2:n %城市从第二个开始for i=1:mvisited=Tabu(i,1:(j-1)); %已访问的城市J=zeros(1,(n-j+1)); %待访问的城市P=J; %待访问城市的选择概率分布（初始化）Jc=1; %循环下标for k=1:n %利用循环求解待访问的城市，如果第k个城市不属于已访问的城市，则其为待访问的城市if isempty(find(visited==k, 1))% if length(find(visited==k))==0J(Jc)=k;Jc=Jc+1; %下标加1，便于下一步存储待访问的城市endend%下面计算待访问城市的概率分布for k=1:length(J) %length(J)表示待访问的城市的个数P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta); %概率计算公式中的分子end %Tau为信息素矩阵,Eta为启发因子矩阵P=P/(sum(P)); %概率分布：长度为待访问城市个数%按概率原则选取下一个城市Pcum=cumsum(P); %cumsum求累加和: cumsum([1 1 1])= 1 2 3，求累加的目的在于使Pcum的值总有大于rand的数Select=find(Pcum>=rand); %当累积概率和大于给定的随机数，则选择个被加上的最后一个城市作为即将访问的城市to_visit=J(Select(1)); %to_visit表示即将访问的城市Tabu(i,j)=to_visit; %将访问过的城市加入禁忌表中endendif NC>=2 %如果迭代次数大于等于2，则将上一次迭代的最佳路线存入Tabu的第一行中Tabu(1,:)=R_best(NC-1,:);end%% 第四步：记录本次迭代最佳路线L=zeros(m,1);for i=1:mR=Tabu(i,:);for j=1:(n-1)L(i)=L(i)+D(R(j),R(j+1)); %求路径距离endL(i)=L(i)+D(R(1),R(n)); %加上最后一个城市与第一个城市之间的距离endL_best(NC)=min(L); %最优路径为距离最短的路径pos=find(L==L_best(NC)); %找出最优路径对应的位置：即是哪只个蚂蚁R_best(NC,:)=Tabu(pos(1),:); %确定最优路径对应的城市顺序L_ave(NC)=mean(L); %求第k次迭代的平均距离NC=NC+1;%% 第五步：更新信息素Delta_T au=zeros(n,n); %Delta_Tau（i,j）表示所有的蚂蚁留在第i个城市到第j个城市路径上的信息素增量for i=1:mfor j=1:(n-1) %建立了完整路径后路径后在释放信息素：蚁周系统Q/L Delta_T au(Tabu(i,j),T abu(i,j+1))=Delta_Tau(Tabu(i,j),T abu(i,j+1))+Q/L(i);endDelta_T au(Tabu(i,n),T abu(i,1))=Delta_T au(Tabu(i,n),T abu(i,1))+Q/L(i);endTau=(1-Rho).*T au+Delta_T au; %信息素更新公式%% 第六步：禁忌表清零Tabu=zeros(m,n); %每迭代一次都将禁忌表清零end%% 第七步：输出结果Pos=find(L_best==min(L_best)); %找到L_best中最小值所在的位置并赋给Pos Shortest_Route=R_best(Pos(1),:); %提取最短路径Shortest_Length=L_best(Pos(1)); %提取最短路径的长度C: 函数文件2%% 绘制各个位置坐标即利用ACA找到的最优路径的函数DrawCity.mfunction DrawCity(C,T,R)%% C Coordinate 节点坐标，由一个N×2的矩阵存储%% T text 各城市的说明%% R Route 路线%%================================================== ==================N=length(R);scatter(C(:,1),C(:,2),'*','LineWidth',3);hold on; %绘制散点图scatter(C(:,1),C(:,2),'r','LineWidth',2);hold on; %绘制散点图%[C(R(1),1),C(R(N),1)]最佳路径第一个城市和最后一个城市的横坐标plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)],'m','LineWidth',2);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)],'-','LineWidth',2); hold on; endfor k=1:length(C(:,1)) %城市标注text(C(k,1)+0.2,C(k,2)+0.3,T(2*k-1:2*k)); hold on;end% xlim([85 130]); ylim([15 50]); %限定横纵坐标显示范围xlabel('东经'); ylabel('北纬'); title('ACA找到的最佳路径'); grid on %绘制网络线。