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

蚁群算法原理及在TSP中的应用

1 蚁群算法（ACA）原理

1.1 基本蚁群算法的数学模型

以求解平面上一个n阶旅行商问题(Traveling Salesman Problem，TSP)为例来说明蚁群算法ACA（Ant Colony Algorithm）的基本原理。对于其他问题，可以对此模型稍作修改便可应用。TSP问题就是给定一组城市，求一条遍历所有城市的最短回路问题。

设bi(t)表示t时刻位于元素i的蚂蚁数目，ij(t)为t时刻路径(i,j)上的信息量，n表示TSP规模，m为蚁群的总数目，则mbi(t)；{ij(t)ci,ciC}是

n

t时刻集合C中元素(城市)两两连接tij上残留信息量的集合。在初始时刻各条路径上信息量相等，并设 ij(0)const，基本蚁群算法的寻优是通过有向图

g(C,L,)实现的。

i1

蚂蚁k(k1,2,...,m)在运动过程中，根据各条路径上的信息量决定其转移方向。这里用禁忌表tabuk(k1,2,...,m)来记录蚂蚁k当前所走过的城市，集合随着

tabuk进化过程作动态调整。在搜索过程中，蚂蚁根据各条路径上的信息量及路

k

径的启发信息来计算状态转移概率。pij(t)表示在t时刻蚂蚁k由元素(城市)i转移

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



(t)*(t)ijij

k

(t)*(t)pij(t)ijij

sallowedk0

若jallowedk否则

(1)

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

ij(t)

1

(2) dij

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

k

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

市) j的期望程度。

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

ij(tn)(1)*ij(t)ij(t) (3)

k

ij(t)ij(t) (4)

k1m

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

k

息素增量，初始时刻ij(t)0，ij(t)表示第k只蚂蚁在本次循环中留在路径

(i,j)上的信息量。

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

k在于ij(t)求法的不同。

在Ant-Cycle模型中

Q

,k

ij(t)Lk

0

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

(5)

，否则

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

在Ant-Quantity模型中

Qd,k

ij(t)ij

0

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

(6)

，否则

在Ant-Density模型中

Q,k

ij(t)

0

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

(7)

，否则

区别：式(6)和式(7)中利用的是局部信息，即蚂蚁完成一步后更新路径上的

信息素；而式(5)中利用的是整体信息，即蚂蚁完成一个循环后更新所有路径上的信息素，在求解TSP时性能较好，因此通常采用式(5)作为蚁群算法的基本模型。

1.2 基本蚁群算法的实现

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

(1) 参数初始化。令时间t=0和循环次数Nc0，设置最大循环次数Ncmax，将m蚂蚁置于n个元素(城市)上，另有向图上每条边(i,j)的初始化信息量

ij(t)const，其中const表示常数，且初始时刻ij(0)0。

(2) 循环次数NcNc1。 (3) 蚂蚁等禁忌表索引号k=1。 (4) 蚂蚁数目kk1。

(5) 蚂蚁个体根据状态转移概率公式(1)计算的概率选择元素(城市)j并前进，jCtabuk

(6) 修改禁忌表指针，即选择好之后将蚂蚁移动到新的元素(城市)，并把该元素(城市)移动到该蚂蚁个体的禁忌表中。

(7) 若集合C中元素(城市)未遍历完，即k

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

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

2 程序实现（城市为中国各省省会城市）

2.1数据列表

TSP问题中的城市选为中国各省省会城市，其实际地理位置的经纬度如表1所示。

表1 中国各省省会城市所在地理位置的经纬度

2.2程序实现

程序说明：程序是在GreenSim团队编写的程序的基础上进行修改得到，区别在于在程序中加入了更多的说明，便于MATLAB初学者以及对蚁群算法不是很熟悉的读者理解；对程序进行了略微的修改，加入了主程序和真实城市坐标数据。其中主程序名称为ACAmain_city.m，其是可直接运行的程序。运行结果如图1-图3所示，结果经过略微调整。图2为将省会城市改成相应的省（或直辖市、特别行政区）后的图形。结果均为一次运行所得，但由于算法的随机性，每次运行结果可能会不同。另外，算法的初始条件对运行结果有很大影响，各种初始条

件取值在附录中的程序已给出。

图1 ACA找到的最佳的中国环游路径（省会城市）

图2 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

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_Tau=zeros(n,n); %Delta_Tau（i,j）表示所有的蚂蚁留在第i个城市到第j个城市路径上的信息素增量 for i=1:m

for j=1:(n-1) %建立了完整路径后路径后在释放信息素：蚁周系统Q/L Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i); 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)); %找到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 %绘制网络线

蚁群算法原理及在TSP中的应用

1 蚁群算法（ACA）原理

1.1 基本蚁群算法的数学模型

以求解平面上一个n阶旅行商问题(Traveling Salesman Problem，TSP)为例来说明蚁群算法ACA（Ant Colony Algorithm）的基本原理。对于其他问题，可以对此模型稍作修改便可应用。TSP问题就是给定一组城市，求一条遍历所有城市的最短回路问题。

设bi(t)表示t时刻位于元素i的蚂蚁数目，ij(t)为t时刻路径(i,j)上的信息量，n表示TSP规模，m为蚁群的总数目，则mbi(t)；{ij(t)ci,ciC}是

n

t时刻集合C中元素(城市)两两连接tij上残留信息量的集合。在初始时刻各条路径上信息量相等，并设 ij(0)const，基本蚁群算法的寻优是通过有向图

g(C,L,)实现的。

i1

蚂蚁k(k1,2,...,m)在运动过程中，根据各条路径上的信息量决定其转移方向。这里用禁忌表tabuk(k1,2,...,m)来记录蚂蚁k当前所走过的城市，集合随着

tabuk进化过程作动态调整。在搜索过程中，蚂蚁根据各条路径上的信息量及路

k

径的启发信息来计算状态转移概率。pij(t)表示在t时刻蚂蚁k由元素(城市)i转移

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



(t)*(t)ijij

k

(t)*(t)pij(t)ijij

sallowedk0

若jallowedk否则

(1)

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

ij(t)

1

(2) dij

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

k

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

市) j的期望程度。

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

ij(tn)(1)*ij(t)ij(t) (3)

k

ij(t)ij(t) (4)

k1m

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

k

息素增量，初始时刻ij(t)0，ij(t)表示第k只蚂蚁在本次循环中留在路径

(i,j)上的信息量。

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

k在于ij(t)求法的不同。

在Ant-Cycle模型中

Q

,k

ij(t)Lk

0

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

(5)

，否则

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

在Ant-Quantity模型中

Qd,k

ij(t)ij

0

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

(6)

，否则

在Ant-Density模型中

Q,k

ij(t)

0

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

(7)

，否则

区别：式(6)和式(7)中利用的是局部信息，即蚂蚁完成一步后更新路径上的

信息素；而式(5)中利用的是整体信息，即蚂蚁完成一个循环后更新所有路径上的信息素，在求解TSP时性能较好，因此通常采用式(5)作为蚁群算法的基本模型。

1.2 基本蚁群算法的实现

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

(1) 参数初始化。令时间t=0和循环次数Nc0，设置最大循环次数Ncmax，将m蚂蚁置于n个元素(城市)上，另有向图上每条边(i,j)的初始化信息量

ij(t)const，其中const表示常数，且初始时刻ij(0)0。

(2) 循环次数NcNc1。 (3) 蚂蚁等禁忌表索引号k=1。 (4) 蚂蚁数目kk1。

(5) 蚂蚁个体根据状态转移概率公式(1)计算的概率选择元素(城市)j并前进，jCtabuk

(6) 修改禁忌表指针，即选择好之后将蚂蚁移动到新的元素(城市)，并把该元素(城市)移动到该蚂蚁个体的禁忌表中。

(7) 若集合C中元素(城市)未遍历完，即k

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

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

2 程序实现（城市为中国各省省会城市）

2.1数据列表

TSP问题中的城市选为中国各省省会城市，其实际地理位置的经纬度如表1所示。

表1 中国各省省会城市所在地理位置的经纬度

2.2程序实现

程序说明：程序是在GreenSim团队编写的程序的基础上进行修改得到，区别在于在程序中加入了更多的说明，便于MATLAB初学者以及对蚁群算法不是很熟悉的读者理解；对程序进行了略微的修改，加入了主程序和真实城市坐标数据。其中主程序名称为ACAmain_city.m，其是可直接运行的程序。运行结果如图1-图3所示，结果经过略微调整。图2为将省会城市改成相应的省（或直辖市、特别行政区）后的图形。结果均为一次运行所得，但由于算法的随机性，每次运行结果可能会不同。另外，算法的初始条件对运行结果有很大影响，各种初始条

件取值在附录中的程序已给出。

图1 ACA找到的最佳的中国环游路径（省会城市）

图2 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

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_Tau=zeros(n,n); %Delta_Tau（i,j）表示所有的蚂蚁留在第i个城市到第j个城市路径上的信息素增量 for i=1:m

for j=1:(n-1) %建立了完整路径后路径后在释放信息素：蚁周系统Q/L Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i); 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)); %找到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 %绘制网络线