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
%%-------------------------------------------------------------------------
%% 主要符號(hào)說(shuō)明
%% C n個(gè)城市的坐標(biāo)，n×2的矩陣
%% NC_max 最大迭代次數(shù)
%% m 螞蟻個(gè)數(shù)
%% Alpha 表征信息素重要程度的參數(shù)
%% Beta 表征啟發(fā)式因子重要程度的參數(shù)
%% Rho 信息素蒸發(fā)系數(shù)
%% Q 信息素增加強(qiáng)度系數(shù)
%% R_best 各代最佳路線
%% L_best 各代最佳路線的長(zhǎng)度
%%=========================================================================

%%第一步：變量初始化
n=size(C,1);%n表示問(wèn)題的規(guī)模（城市個(gè)數(shù)）
D=zeros(n,n);%D表示完全圖的賦權(quán)鄰接矩陣
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為啟發(fā)因子，這里設(shè)為距離的倒數(shù)
Tau=ones(n,n);%Tau為信息素矩陣
Tabu=zeros(m,n);%存儲(chǔ)并記錄路徑的生成
NC=1;%迭代計(jì)數(shù)器
R_best=zeros(NC_max,n);%各代最佳路線
L_best=inf.*ones(NC_max,1);%各代最佳路線的長(zhǎng)度
L_ave=zeros(NC_max,1);%各代路線的平均長(zhǎng)度

while NC<=NC_max%停止條件之一：達(dá)到最大迭代次數(shù)
%%第二步：將m只螞蟻放到n個(gè)城市上
Randpos=[];
for i=1:(ceil(m/n))
Randpos=[Randpos,randperm(n)];
end
Tabu(:,1)=(Randpos(1,1:m))';

%%第三步：m只螞蟻按概率函數(shù)選擇下一座城市，完成各自的周游
for j=2:n
for i=1:m
visited=Tabu(i,1:(j-1));%已訪問(wèn)的城市
J=zeros(1,(n-j+1));%待訪問(wèn)的城市
P=J;%待訪問(wèn)城市的選擇概率分布
Jc=1;
for k=1:n
if length(find(visited==k))==0
J(Jc)=k;
Jc=Jc+1;
end
end
%下面計(jì)算待選城市的概率分布
for k=1:length(J)
P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta);
end
P=P/(sum(P));
%按概率原則選取下一個(gè)城市
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

%%第七步：輸出結(jié)果
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
%% 畫路線圖的子函數(shù)
%%-------------------------------------------------------------------------
%% C Coordinate 節(jié)點(diǎn)坐標(biāo)，由一個(gè)N×2的矩陣存儲(chǔ)
%% 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

設(shè)置初始參數(shù)如下：
m=31;Alpha=1;Beta=5;Rho=0.1;NC_max=200;Q=100;
31城市坐標(biāo)為：
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