%初始化
clear;

Alpha=1; %信息素重要程度的參數(shù)
Beta=5; %啟發(fā)式因子重要程度的參數(shù) 
Rho=0.95; %信息素蒸發(fā)系數(shù)
NC_max=200; %最大迭代次數(shù)
Q=100; %信息素增加強度系數(shù)
CityNum=30;  %問題的規(guī)模（城市個數(shù)）
[dislist,Clist]=tsp(CityNum);
m=CityNum; %螞蟻個數(shù)
Eta=1./dislist;%Eta為啟發(fā)因子，這里設(shè)為距離的倒數(shù)
Tau=ones(CityNum,CityNum);%Tau為信息素矩陣
Tabu=zeros(m,CityNum);%存儲并記錄路徑的生成
NC=1;%迭代計數(shù)器
R_best=zeros(NC_max,CityNum); %各代最佳路線
L_best=inf.*ones(NC_max,1);%各代最佳路線的長度
L_ave=zeros(NC_max,1);%各代路線的平均長度

figure(1);
while NC<=NC_max %停止條件之一：達到最大迭代次數(shù)
    %將m只螞蟻放到CityNum個城市上
    Randpos=[];
    for i=1:(ceil(m/CityNum))
        Randpos=[Randpos,randperm(CityNum)];
    end
    Tabu(:,1)=(Randpos(1,1:m))';
    
    %m只螞蟻按概率函數(shù)選擇下一座城市，完成各自的周游
    for j=2:CityNum
        for i=1:m
            visited=Tabu(i,1:(j-1)); %已訪問的城市
            J=zeros(1,(CityNum-j+1));%待訪問的城市
            P=J;%待訪問城市的選擇概率分布
            Jc=1;
            for k=1:CityNum
                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,:);
        L(i)=CalDist(dislist,R);
    end
    L_best(NC)=min(L);
    pos=find(L==L_best(NC));
    R_best(NC,:)=Tabu(pos(1),:);
    L_ave(NC)=mean(L);
    drawTSP(Clist,R_best(NC,:),L_best(NC),NC,0);
    NC=NC+1;
    %更新信息素
    Delta_Tau=zeros(CityNum,CityNum);
    for i=1:m
        for j=1:(CityNum-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,CityNum),Tabu(i,1))=Delta_Tau(Tabu(i,CityNum),Tabu(i,1))+Q/L(i);
    end
    Tau=(1-Rho).*Tau+Delta_Tau;
    Tabu=zeros(m,CityNum); %禁忌表清零
    %pause;
    tauji(NC)=Tau(1,2);
end

%輸出結(jié)果
Pos=find(L_best==min(L_best));
Shortest_Route=R_best(Pos(1),:);
Shortest_Length=L_best(Pos(1));
figure(2);
plot([L_best L_ave]);
legend('最短距離','平均距離');