function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)
%%=========================================================================
%% ACATSP.m
%% Ant Colony Algorithm for Traveling Salesman Problem
%% ChengAihua,PLA Information Engineering University,ZhengZhou,China
%% Email:aihuacheng@gmail.com
%% All rights reserved
%%-------------------------------------------------------------------------
%% 主要符號說明
%% 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;
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
J(Jc)=k;
Jc=Jc+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);
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);
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);
NC=NC+1

%%第五步：更新信息素
Delta_Tau=zeros(n,n);
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);
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));
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)

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)])
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)])
hold on
end

設置初始參數如下：
m=31;Alpha=1;Beta=5;Rho=0.1;NC_max=200;Q=100;
31城市坐標為：
1304  2312
3639  1315
4177  2244
3712  1399
3488  1535
3326  1556
3238  1229
4196  1004
4312  790
4386  570
3007  1970
2562  1756
2788  1491
2381  1676
1332  695
3715  1678
3918  2179
4061  2370
3780  2212
3676  2578
4029  2838
4263  2931
3429  1908
3507  2367
3394  2643
3439  3201
2935  3240
3140  3550
2545  2357
2778  2826
2370  2975