%帶精英策略的蟻群算法求解TSP問題的matlab程序
clear all
close all
clc

%初始化蟻群
m=31;%蟻群中螞蟻的數量，當m接近或等于城市個數n時，本算法可以在最少的迭代次數內找到最優解
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];%城市的坐標矩陣
Nc_max=200;%最大循環次數，即算法迭代的次數，亦即螞蟻出動的撥數（每撥螞蟻的數量當然都是m）
alpha=1;%螞蟻在運動過程中所積累信息（即信息素）在螞蟻選擇路徑時的相對重要程度，alpha過大時，算法迭代到一定代數后將出現停滯現象
beta=5;%啟發式因子在螞蟻選擇路徑時的相對重要程度
rho=0.5;%0<rho<1,表示路徑上信息素的衰減系數（亦稱揮發系數、蒸發系數），1-rho表示信息素的持久性系數
Q=100;%螞蟻釋放的信息素量，對本算法的性能影響不大

%變量初始化
n=size(C,1);%表示TSP問題的規模，亦即城市的數量
D=ones(n,n);%表示城市完全地圖的賦權鄰接矩陣，記錄城市之間的距離
for i=1:n
    for j=1:n
        if i<j
            D(i,j)=sqrt((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2);
        end
    D(j,i)=D(i,j);
    end
end
eta=1./D;%啟發式因子，這里設為城市之間距離的倒數
pheromone=ones(n,n);%信息素矩陣,這里假設任何兩個城市之間路徑上的初始信息素都為1
tabu_list=zeros(m,n);%禁忌表，記錄螞蟻已經走過的城市，螞蟻在本次循環中不能再經過這些城市。當本次循環結束后，禁忌表被用來計算螞蟻當前所建立的解決方案，即經過的路徑和路徑的長度
Nc=0;%循環次數計數器
routh_best=zeros(Nc_max,n);%各次循環的最短路徑
length_best=ones(Nc_max,1);%各次循環最短路徑的長度
length_average=ones(Nc_max,1);%各次循環所有路徑的平均長度

while Nc<Nc_max
    %將m只螞蟻放在n個城市上
    rand_position=[];
    for i=1:ceil(m/n)
        rand_position=[rand_position,randperm(n)];
    end
    tabu_list(:,1)=(rand_position(1:m))';%將螞蟻放在城市上之后的禁忌表，第i行表示第i只螞蟻，第i行第一列元素表示第i只螞蟻所在的初始城市
    %m只螞蟻按概率函數選擇下一座城市，在本次循環中完成各自的周游
    for j=2:n
        for i=1:m
            city_visited=tabu_list(i,1:(j-1));%已訪問的城市
            city_remained=zeros(1,(n-j+1));%待訪問的城市
            probability=city_remained;%待訪問城市的訪問概率
            cr=1;
            for k=1:n%for循環用于求待訪問的城市。比如如果城市個數是5，而已訪問的城市city_visited=[2 4],則經過此for循環后city_remanied=[1 3 5]
                if length(find(city_visited==k))==0
                    city_remained(cr)=k;
                    cr=cr+1;
                end
            end
            %for循環計算待訪問城市的訪問概率分布，此概率和兩個參數有關，一是螞蟻當前所在城市（即city_visited(end))和待訪問城市（即city_remained（k))路徑上的信息素，一是這兩者之間的啟發因子即距離的倒數            
            for k=1:length(city_remained)
                probability(k)=(pheromone(city_visited(end),city_remained(k)))^alpha*(eta(city_visited(end),city_remained(k)))^beta;
            end
            probability=probability/sum(probability);
            %按概率選取下一個要訪問的城市????????????????????????????????
            pcum=cumsum(probability);
            select=find(pcum>=rand);
            to_visit=city_remained(select(1));
            tabu_list(i,j)=to_visit;
        end
    end
    if Nc>0
        tabu_list(1,:)=routh_best(Nc,:);%將上一代的最優路徑（最優解）保留下來，保證上一代中的最適應個體的信息不會丟失
    end
    %記錄本次循環的最佳路線
    total_length=zeros(m,1);%m只螞蟻在本次循環中分別所走過的路徑長度
    for i=1:m
        r=tabu_list(i,:);%取出第i只螞蟻在本次循環中所走的路徑
        for j=1:(n-1)
            total_length(i)=total_length(i)+D(r(j),r(j+1));%第i只螞蟻本次循環中從起點城市到終點城市所走過的路徑長度
        end
        total_length(i)=total_length(i)+D(r(1),r(n));%最終得到第i只螞蟻在本次循環中所走過的路徑長度
    end
    length_best(Nc+1)=min(total_length);%把m只螞蟻在本次循環中所走路徑長度的最小值作為本次循環中最短路徑的長度
    position=find(total_length==length_best(Nc+1));%找到最短路徑的位置，即最短路徑是第幾只螞蟻或哪幾只螞蟻走出來的
    routh_best(Nc+1,:)=tabu_list(position(1),:);%把第一個走出最短路徑的螞蟻在本次循環中所走的路徑作為本次循環中的最優路徑
    length_average(Nc+1)=mean(total_length);%計算本次循環中m只螞蟻所走路徑的平均長度
    %更新信息素
    delta_pheromone=zeros(n,n);
    for i=1:m
        for j=1:(n-1)
            delta_pheromone(tabu_list(i,j),tabu_list(i,j+1))=delta_pheromone(tabu_list(i,j),tabu_list(i,j+1))+Q/total_length(i);%total_length(i)為第i只螞蟻在本次循環中所走過的路徑長度（蟻周系統區別于蟻密系統和蟻量系統的地方）
        end
        delta_pheromone(tabu_list(i,n),tabu_list(i,1))=delta_pheromone(tabu_list(i,n),tabu_list(i,1))+Q/total_length(i);%至此把第i只螞蟻在本次循環中在所有路徑上釋放的普通信息素（區別于下面的額外信息素）已經累加上去
    end%至此把m只螞蟻在本次循環中在所有路徑上釋放的普通信息素已經累加上去
    
    %%給予最優路徑以額外的信息素
    i=position(1);%找出在本次循環中走出最短路徑的第一只螞蟻
    m_elitist=length(position);%找出精英螞蟻的個數，即本次循環中走出最短路徑的螞蟻個數
    delta_pheromone_elitists=zeros(n,n);
    for j=1:(n-1)
        delta_pheromone_elitists(tabu_list(i,j),tabu_list(i,j+1))=delta_pheromone_elitists(tabu_list(i,j),tabu_list(i,j+1))+m_elitist*Q/length_best(Nc+1);%Q/length_best(Nc+1)為本次循環中精英螞蟻找出的最優路徑的長度
    end
    delta_pheromone_elitists(tabu_list(i,n),tabu_list(i,1))=delta_pheromone_elitists(tabu_list(i,n),tabu_list(i,1))+m_elitist*Q/length_best(Nc+1);%至此在最優路徑的每條邊上釋放了額外的信息素
  
    pheromone=(1-rho).*pheromone+delta_pheromone+delta_pheromone_elitists;%本次循環后所有路徑上的信息素
    
    %循環次數加1，禁忌表清零，準備下一次循環，螞蟻在下一次循環中又可以自由地進行選擇
    Nc=Nc+1
    tabu_list=zeros(m,n);
end

%輸出結果，繪制圖形
position=find(length_best==min(length_best));
shortest_path=routh_best(position(1),:)
shortest_length=length_best(position(1))
%繪制最短路徑
figure(1)
set(gcf,'Name','Ant Colony Optimization——Figure of shortest_path','Color','r')
N=length(shortest_path);
scatter(C(:,1),C(:,2),50,'filled');
hold on
plot([C(shortest_path(1),1),C(shortest_path(N),1)],[C(shortest_path(1),2),C(shortest_path(N),2)])
set(gca,'Color','g')
hold on
for i=2:N
    plot([C(shortest_path(i-1),1),C(shortest_path(i),1)],[C(shortest_path(i-1),2),C(shortest_path(i),2)])
    hold on
end
%繪制每次循環最短路徑長度和平均路徑長度
figure(2)
set(gcf,'Name','Ant Colony Optimization——Figure of length_best and length_average','Color','r')
plot(length_best,'r')
set(gca,'Color','g')
hold on
plot(length_average,'k')