蚁群算法原理及在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()n

i 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 进化过程作动态调整。在搜索过程中，蚂蚁根据各条路径上的信息量及路

径的启发信息来计算状态转移概率。()k

ij p t 表示在t 时刻蚂蚁k 由元素(城市)i 转移

到元素(城市)j 的状态转移概率。

()*()()*()()0k ij ij k k

ij ij ij s allowed t t j allowed t t p t αβ

αβτητη???????????

∈?????=?????

???

∑若否则

(1)

式中，{}k k allowed C tabuk =-表示蚂蚁k 下一步允许选择的城市；α为信息启发式因子，表示轨迹的相对重要性，反映了蚂蚁在运动过程中所积累的信息在蚂蚁运动时所起作用，其值越大，则该蚂蚁越倾向于选择其他蚂蚁经过的路径，蚂蚁之间协作性越强；β为期望启发式因子，表示能见度的相对重要性，反映了蚂蚁在运动过程中启发信息在蚂蚁选择路径中的重视程度，其值越大，则该状态转移概率越接近于贪心规则；()ij t η为启发函数，其表达式如下：

1

()ij ij

t d η=

(2)

式中，ij d 表示相邻两个城市之间的距离。对蚂蚁k 而言，ij d 越小，则()ij t η越

大，()k

ij p t 也就越大。显然，该启发函数表示蚂蚁从元素(城市) i 转移到元素(城

市) j 的期望程度。

为了避免残留信息素过多引起残留信息淹没启发信息，在每只蚂蚁走完一步或者完成对所有n 个城市的遍历(也是一个循环结束)后，要对残留信息进行更新处理。这种更新策略模仿了人类大脑记忆的特点，在新信息不断存入大脑的同时，存储在大脑中的旧信息随着时间的推移逐渐淡化，甚至忘记。由此，t n +时刻在路径(,)i j 上的信息量可按如下规则进行调整：

()(1)*()()ij ij ij t n t t τρττ+=-+? (3)

1()()m

k

ij ij k t t ττ=?=?∑ (4)

式中ρ表示信息挥发系数，则1ρ-表示信息素残留因子，为了防止信息的无限积累，ρ的取值范围为：[0,1)ρ?；()ij t τ?表示本次循环中路径(,)i j 上的信

息素增量，初始时刻()0ij t τ?=，()k

ij t τ?表示第k 只蚂蚁在本次循环中留在路径

(,)i j 上的信息量。

根据信息素更新策略的不同，Dorigo M 提出了三种不同的基本蚁群算法模型，分别称之为Ant-Cycle 模型、Ant-Quantity 模型及Ant-Density 模型，其

差别在于()k ij t τ?求法的不同。

在Ant-Cycle 模型中

,k i j ()0k

k

ij Q

L t τ???=???

若第只蚂蚁在本次循环中经过（，）

，否则

(5)

式中，Q 表示信息素强度，它在一定程度上影响算法的收敛速度；k L 表示k 只蚂蚁在本次循环中所走路径的总长度。

在Ant-Quantity 模型中

,k i j ()0k

ij

ij Q d t τ???=???

若第只蚂蚁在t 和t+1之间经过（，）

，否则

(6)

在Ant-Density 模型中

,k i j ()0

k

ij 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

(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.m

clc; 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.33

37.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; toc

B: 函数文件1

%% 利用蚁群算法解决TSP问题的函数文件ACATSP.m

function

[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: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; % eps = 2.2204e-016，i=j,则距离为0

end

D(j,i)=D(i,j); % 距离矩阵为对称均值（n*n的矩阵）

end

end

Eta=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)：随机产生一个整数排列end

Tabu(:,1)=(Randpos(1,1:m))'; %每只蚂蚁（m只）都对应有一个位置，Tabu(:,1)为每只蚂蚁的走过的第一个城市

%% 第三步：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 %利用循环求解待访问的城市，如果第k个城市不属于已访问的城市，则其为待访问的城市

if isempty(find(visited==k, 1))

% if length(find(visited==k))==0

J(Jc)=k;

Jc=Jc+1; %下标加1，便于下一步存储待访问的城市

end

end

%下面计算待访问城市的概率分布

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; %将访问过的城市加入禁忌表中

end

end

if NC>=2 %如果迭代次数大于等于2，则将上一次迭代的最佳路线存入Tabu的第一行中

Tabu(1,:)=R_best(NC-1,:);

end

%% 第四步：记录本次迭代最佳路线

L=zeros(m,1);

for i=1:m

R=Tabu(i,:);

for j=1:(n-1)

L(i)=L(i)+D(R(j),R(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); %求第k次迭代的平均距离

NC=NC+1;

%% 第五步：更新信息素

Delta_T au=zeros(n,n); %Delta_Tau（i,j）表示所有的蚂蚁留在第i个城市到第j个城市路径上的信息素增量

for i=1:m

for 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);

end

Delta_T au(Tabu(i,n),T abu(i,1))=Delta_T au(Tabu(i,n),T abu(i,1))+Q/L(i);

end

Tau=(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.m

function 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; end

for 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 %绘制网络线