蟻群演算法簡介

蟻群是自然界中常見的一種生物，人們對螞蟻的關注大都是因為“蟻群搬家，天要下雨”之類的民諺。然而隨著近代仿生學的發展，這種似乎微不足道的小東西越來越多地受到學者們地關注。 1991 年義大利學者 M. Dorigo 等人首先提出了蟻群演算法，人們開始了對蟻群的研究：相對弱小，功能並不強大的個體是如何完成複雜的工作的（如尋找到食物的最佳路徑並返回等）。在此基礎上一種很好的優化演算法逐步發展起來。
蟻群演算法的特點是模擬自然界中螞蟻的群體行為。科學家發現，蟻群總是能夠發現從蟻巢到食物源的最短路徑。經研究發現，螞蟻在行走過的路上留下一種揮發性的激素，螞蟻就是通過這種激素進行資訊交流。螞蟻趨向於走激素積累較多的路徑。找到最短路徑的螞蟻總是最早返回巢穴，從而在路上留下了較多的激素。由於最短路徑上積累了較多的激素，選擇這條路徑的螞蟻就會越來越多，到最後所有的螞蟻都會趨向於選擇這條最短路徑。基於螞蟻這種行為而提出的蟻群演算法具有群體合作，正反饋選擇，平行計算等三大特點，並且可以根據需要為人工蟻加入前瞻、回溯等自然蟻所沒有的特點。
在使用蟻群演算法求解現實問題時，先生成具有一定數量螞蟻的蟻群，讓每一隻螞蟻建立一個解或解的一部分，每隻人工蟻從問題的初始狀態出發，根據“激素”濃度來選擇下一個要轉移到的狀態，直到建立起一個解，每隻螞蟻根據所找到的解的好壞程度在所經過的狀態上釋放與解的質量成正比例的“激素”。之後，每隻螞蟻又開始新的求解過程，直到尋找到滿意解。為避免停滯現象，引入了激素更新機制。

摘來理論部分（我比較懶，圖片是摘自司守奎的《數學建模演算法與應用》）

程式碼說明

在網上無意間找到了這個程式碼，寫的非常不錯：依照原理論一板一眼的寫的程式碼，沒有偷工減料，很難得。但是看起來比較費勁，就自己改了下：

• 刪改了原始碼批註，並增加了每一行程式碼的註釋，相信再接觸這份程式碼理解起來會快一些
• 精簡了兩三處程式碼
• 改了下變數名稱，總覺得原來的看著亂

Matlab程式碼

主呼叫檔案：

City = rand(20,2).*randi([1,100],20,2);
Ite = 200;
Ant_num = 40;
Alpha = 0.7;
Beta = 0.9
 
;
Rho = 0.3;
Q = 100;
[R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(City,Ite,Ant_num,Alpha,Beta,Rho,Q)

主程式：

function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(City,Ite,Ant_num,Alpha,Beta,Rho,Q)
%% 主要符號說明
% City n個城市的座標，n×2的矩陣
% Ite 最大迭代次數
% Ant_num 螞蟻個數
% Alpha 表徵資訊素重要程度的引數
 

% Beta 表徵啟發式因子重要程度的引數
% Eta 啟發因子，這裡設為距離的倒數
% Tau 資訊素矩陣
% Rho 資訊素蒸發係數
% Q 資訊素增加強度係數
% R_best 各代最佳路線
% L_best 各代最佳路線的長度

%% 第一步：變數初始化
City_num = size(City,1);%n表示問題的規模（城市個數）
Distance = zeros(City_num,City_num);%D表示完全圖的賦權鄰接矩陣

for i = 1:City_num
    for j = 1:City_num
       if i ~= j
            Distance(i,j) = sqrt( (City(i,1)-City(j,1))^2 + (City(i,2)-City(j,2))^2 );% 計算任意兩點間距離
       else
            Distance(i,j) = eps;      % i=j時不計算，應該為0，但後面的啟發因子要取倒數，用eps（浮點相對精度）表示
       end
    Distance(j,i) = Distance(i,j);   % 對稱矩陣
    end
end

Eta = 1./Distance;          % Eta為啟發因子，這裡設為距離的倒數
Tau = ones(City_num,City_num);     % Tau為資訊素矩陣
Ite_num = 1;               % 迭代計數器，記錄迭代次數
R_rec = zeros(Ant_num,City_num);   % 儲存並記錄路徑的生成
R_best = zeros(Ite,City_num);       % 各代最佳路線
L_best = inf.*ones(Ite,1);   % 各代最佳路線的長度
L_ave = zeros(Ite,1);        % 各代路線的平均長度

while Ite_num<=Ite        % 停止條件之一：達到最大迭代次數，停止

%% 第二步：將Ant_num只螞蟻放到City_num個城市上
    Init_ant_position = [];   % 將初始狀態下的螞蟻隨機分配到各個城市的臨時矩陣。
    for i = 1:( ceil(Ant_num/City_num) )
        Init_ant_position = [Init_ant_position,randperm(City_num)]; 
        % 每次迭代將螞蟻隨機分配到所有城市。生成一個1行多列（由於ceiling向上取整，則列數大於等於螞蟻個數）的矩陣。
    end    
    R_rec(:,1) = (Init_ant_position(1,1:Ant_num))';   
    % 矩陣轉置後變成 Ant_num行-1列的一個一個初始矩陣，每一行代表一隻螞蟻，每個值代表當前螞蟻所在城市程式碼。

%% 第三步：Ant_num只螞蟻按概率函式選擇下一座城市，完成各自的周遊
    % 這裡說明下：這是個兩重的大迴圈，外層是城市，內層是螞蟻。
    % 也就是說每次迭代每隻螞蟻向前走一步，而不是一隻螞蟻走完所有城市再開始下一隻。
    for j = 2:City_num     % 所在城市不計算
        for i = 1:Ant_num    % 對每一隻螞蟻
            City_visited = R_rec(i,1:(j-1));   % 記錄已訪問的城市，避免重複訪問
            City_unvisited = zeros(1,(City_num-j+1));  % 待訪問的城市，初始為空
            P = City_unvisited;   % 待訪問城市的選擇概率分佈，我猜這裡作者弄了個簡便寫法，其實只是想弄一個同型矩陣。
            count = 1; % 待訪問城市 City_unvisited 的下標計數器

            % 統計未去過的城市
            for k = 1:City_num
                if isempty( find( City_visited == k, 1 ) ) 
                    % 如果去過k城市，則find不為空,find(x,1)的意思是找到第一個就結束，是一個提高運算效能的寫法。
                    % 這句話是為了避免重複去同一個城市。
                    City_unvisited(count) = k; % 如果if判斷為真，說明第k 個城市沒有去過。
                    count = count+1;      % 下標計數器加1
                end
            end

            % 下面計算待選城市的概率分佈
            for k = 1:length(City_unvisited)
                P(k) = ( Tau( City_visited(end), City_unvisited(k) )^Alpha )*...
                       ( Eta( City_visited(end), City_unvisited(k) )^Beta );
            end
            P=P/(sum(P));

            % 按概率原則選取下一個城市
            P_cum = cumsum(P);  % cumsum函式是一個比較特殊的求和函式，這裡是得到P 的累積概率矩陣。
            Select = find(P_cum>=rand); % 若計算的概率大於原來的就選擇這條路線
            To_visit = City_unvisited(Select(1)); % Select(1)的意思是選中第一個累積概率大於rand隨機數的城市
            R_rec(i,j) = To_visit;  % R_rec(i,j) 是指第i只螞蟻，第將要去j步將要去的那個城市，迴圈結束後得到每隻螞蟻的路徑
        end
    end

    % 如果不是第一次迴圈，則將最優路徑賦給路徑記錄矩陣的第一行
    if Ite_num >= 2
        R_rec(1,:) = R_best(Ite_num-1,:);
    end

%% 第四步：記錄本次迭代最佳路線
    Len = zeros(Ant_num,1);     % length 距離矩陣，初始為0。記錄每隻螞蟻當前路徑的總距離。

    for i=1:Ant_num
        R_temp = R_rec(i,:); % 取得第i 只螞蟻的路徑
        % 計算第i只螞蟻走過的總距離
        for j = 1:(City_num-1)
            Len(i) = Len(i) + Distance( R_temp(j),R_temp(j+1) );  % 原距離加上第j個城市到第j+1個城市的距離
        end
        Len(i)=Len(i)+Distance(R_temp(1),R_temp(City_num));      % 一輪下來後走過的距離
    end

    [L_best(Ite_num), index] = min(Len);     % 最佳距離取最小
    R_best(Ite_num,:) = R_rec(index(1), :); 
    % 此輪迭代後的最佳路線。為什麼是index(1)，這是嚴謹寫法：因為min求出後如果有多個最小值則index不唯一。
    L_ave(Ite_num) = mean(Len);           % 此輪迭代後的平均距離
    Ite_num=Ite_num+1;                      % 迭代繼續

%% 第五步：更新資訊素
    Delta_Tau = zeros(City_num, City_num);        % 開始時資訊素為n*n的0矩陣

    for i = 1:Ant_num
        for j = 1:(City_num-1)
            Delta_Tau(R_rec(i,j), R_rec(i,j+1)) = Delta_Tau(R_rec(i,j), R_rec(i,j+1))+Q/Len(i);   
            %此次迴圈在路徑（i，j）上的資訊素增量
        end
        Delta_Tau(R_rec(i,City_num), R_rec(i,1)) = Delta_Tau(R_rec(i,City_num), R_rec(i,1))+Q/Len(i);
        %此次迴圈在整個路徑上的資訊素增量
    end

    Tau = (1-Rho).*Tau + Delta_Tau; %考慮資訊素揮發，更新後的資訊素  Rho 資訊素蒸發係數

%% 第六步：禁忌表清零
    R_rec=zeros(Ant_num,City_num);             %%直到最大迭代次數

end

%% 第七步：輸出結果

Pos = find(L_best==min(L_best)); % 找到最佳路徑（非0為真）
Shortest_Route=R_best(Pos(1), :) % 最大迭代次數後最佳路徑
Shortest_Length=L_best(Pos(1) ) % 最大迭代次數後最短距離

figure                % 繪製第一個子圖形
DrawRoute(City,Shortest_Route)     % 畫路線圖的子函式

figure                % 繪製第二個子圖形
plot(L_best)
hold on                         % 保持圖形
plot(L_ave,'r')
title('平均距離和最短距離')     % 標題

繪圖檔案

function DrawRoute(C,R)
%% 畫路線圖的子函式
% C Coordinate 節點座標，由一個N×2的矩陣儲存
% 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)],'g')

for ii=2:N
    plot([C(R(ii-1),1),C(R(ii),1)],[C(R(ii-1),2),C(R(ii),2)],'g')
end
title('旅行商問題優化結果 ')
hold off