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 
%% by wangxuefeng
%%------------------------------------------------------------------------- 
%% 主要符號說明 
%% C n個城市的坐標，n×2的矩陣 
%% NC_max 最大迭代次數 
%% m 螞蟻個數 
%% Alpha 表征信息素重要程度的參數 
%% Beta 表征啟發式因子重要程度的參數 
%% Rho 信息素蒸發系數 
%% Q 信息素增加強度系數 
%% R_best 各代最佳路線 
%% L_best 各代最佳路線的長度 
%%========================================================================= 
m=31;Alpha=1;Beta=5;Rho=0.1;NC_max=200;Q=100; 
C=[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]  

%%第一步：變量初始化 
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)