%% 主要符號說明%% C n個城市的坐標，n×2的矩陣%% NC_max 最大迭代次數%% m 螞蟻個數%% Alpha 表征信息素重要程度的參數%% Beta 表征啟發式因子重要程度的參數%% Rho 信息素蒸發系數%% Q 信息素增加強度系數%% R_best 各代最佳路線%% L_best 各代最佳路線的長度%%=========================================================================%%第一步：變量初始化C=[1304  23123639  13154177  22443712  13993488  15353326  15563238  12294196  10044312  7904386  5703007  19702562  17562788  14912381  16761332  6953715  16783918  21794061  23703780  22123676  25784029  28384263  29313429  19083507  23673394  26433439  32012935  32403140  35502545  23572778  28262370  2975];%城市坐標m=31;Alpha=1;Beta=5;Rho=0.1;NC_max=200;Q=100;n=size(C,1);%n表示問題的規模（城市個數）D=zeros(n,n);%D表示完全圖的賦權鄰接矩陣for i=1:nfor j=1:nif i~=jD(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;elseD(i,j)=eps;%eps為計算機中的最小數endD(j,i)=D(i,j);endendEta=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)];endTabu(:,1)=(Randpos(1,1:m))';%%第三步：m只螞蟻按概率函數選擇下一座城市，完成各自的周游for j=2:nfor i=1:mvisited=Tabu(i,1:(j-1));%已訪問的城市J=zeros(1,(n-j+1));%待訪問的城市P=J;%待訪問城市的選擇概率分布Jc=1;for k=1:nif length(find(visited==k))==0J(Jc)=k;Jc=Jc+1;endend%下面計算待選城市的概率分布for k=1:length(J)P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta);endP=P/(sum(P));%按概率原則選取下一個城市Pcum=cumsum(P);Select=find(Pcum>=rand);to_visit=J(Select(1));Tabu(i,j)=to_visit;endendif NC>=2Tabu(1,:)=R_best(NC-1,:);end%%第四步：記錄本次迭代最佳路線L=zeros(m,1);for i=1:mR=Tabu(i,:);for j=1:(n-1)L(i)=L(i)+D(R(j),R(j+1));endL(i)=L(i)+D(R(1),R(n));endL_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:mfor 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);endDelta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);endTau=(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 onplot(L_ave)