數學建模方法-蟻群算法
一、引言
哈嘍大家好，今天要給大家介紹的是“蟻群算法”。跟粒子群算法一樣，蟻群算法也是基於對生物行為的研究所受到啟發而產生的。它的誕生比粒子群算法還要早3年，在1992年的某一天，一位叫Marco Dorigo的在他的博士論文中提出了蟻群算法，並稱其靈感來源於螞蟻在尋找食物過程中發現路徑的行為。。。好了歷史背景介紹到止，接下來就認真講一下何謂蟻群算法吧。
二、螞蟻尋食
2.1 科普知識
很久以前，生物學家就對螞蟻捕食行為很好奇。因為，可能大家不知道，螞蟻的視力特別差，你可以把它當成瞎子，可是每次卻都能找到食物，並且從它所找到的路徑，還是從洞穴到食物的最短路徑。大家可以想一想，每次你的桌子上有掉一顆糖，過不了多久，就立馬有螞蟻出現了，而且還是一成群的螞蟻沿著同一條線路過來的。這是為什麽呢？難道是螞蟻的鼻子很靈敏嗎？
後來生物學家們經過長時間的研究發現，並不是螞蟻的鼻子很靈敏，而是他們之間有一種很特別的信息傳遞方式——信息素，正是這種特別的信息傳遞方式，使得螞蟻這種生物，不會迷路哈哈。
好言歸正傳，這個信息素是怎麽用的呢？其實是這樣的，螞蟻在行走的過程會釋放一種化學物質，這種化學物質就叫“信息素”，這是用來標識自己的行走路徑。而在尋找食物的過程中，螞蟻根據這個信息素的濃度來選擇行走的方向，最終就到達了食物的所在地。
2.2 講個故事
有一天，有2個螞蟻（分別叫小黃和大黃吧）出來找食物了。從洞穴出來，有兩條路可以走，選擇哪條路的概率是一致的，因為兩條路目前都沒有信息素。我們假設結果是大黃走了直線路，小黃走了曲線路。如下圖
兩個螞蟻在行走過程中分泌信息素，用於標識路線。接著，過了一段時間，大黃找到了奶酪，此時小黃還在行走中。如下圖
接著大黃要回去的時候，由於直線路含有信息素而曲線路沒有信息素（因為當前小黃還沒走到奶酪附近），於是大黃選擇走有信息素的道路也就是直線路。又過了一段時間，大黃回到了路的起點，而小黃剛走到奶酪附近。可以看出現在兩條路的信息素濃度是2: 1。如圖
接著小黃要回去的時候，由於直線路含有的信息素濃度比曲線路高，於是小黃選擇走有信息素濃度高的道路也就是直線路。最終當小黃也回到路的起點的時候，兩條路的信息素濃度比例是3: 1。
以後，每當有螞蟻遇到起點岔路的時候，都會優先選擇信息素濃度高的道路。這樣，一大批的螞蟻就能夠準確的找到最短路徑了。
2.3 啟發
我們從這個螞蟻找食物的過程中，發現其實，這個跟我們之前提到的旅行商問題很像。大家還記得旅行商問題是什麽麽？不記得沒關系，我再講一下就好了。有這麽一個旅行商人，他要拜訪n個城市，但是每次走的路徑是有限制的，即每個城市只能拜訪一次，並且最後要回到原來出發的城市。如何選擇路徑使得總路程為最小值。螞蟻找食物和旅行商問題都是為了求得最短路徑。那接下來，我們就根據利用蟻群算法來解決我們的旅行商問題。誒博主等等你還沒講蟻群算法。別急，我兩者一起講，蟻群算法這東西需要實例比較好解釋。
三、蟻群算法的基本原理
本節以TSP問題為例介紹蟻群算法的原理。
1. 每只螞蟻從一個城市走到另一個城市的過程中都會在路徑上釋放信息素，並且螞蟻選擇下一個城市的依據是一個概率公式，如下：
$ P_{ij}^{k}(t) = \left\{\begin{matrix} \frac{\tau_{ij}^{\alpha}(t)\cdot \eta_{ij}^{\beta}(t)}{\sum_{s \notin tabu_{k}}\tau_{is}^{\alpha}(t)\cdot \eta_{is}^{\beta}(t)} & j \notin tabu_{k}\\ 0 & other \end{matrix}\right.$
上式中，各個符號的意義如下：
- $\alpha$ -- 信息素啟發式因子，它反映了信息素對螞蟻路徑選擇的作用；
- $\beta$ -- 期望啟發式因子，它反映了信息素在螞蟻選擇路徑時被重視的程度；
- $d_{ij}$ -- 城市i和j之間的距離;
- $\eta_{ij}(t)$ -- 啟發函數，表達式為$\eta_{ij}(t) = \frac{1} {d_{ij}}$；
- $tabu_{k}$ -- 禁忌表，記錄螞蟻k當前所走過的城市；
- $\tau_{ij}$ -- 城市i到城市j的路徑上的信息素的量；
2. 螞蟻留下的信息素，因為是化學物質，因此隨著時間的過去信息素也應該要以一定的速率揮發。根據不同的規則我們可以將蟻群算法分為三種模型——蟻周模型(Ant-Cycle)、蟻量模型(Ant-Quantity)和蟻密模型(Ant-Density)。通常使用的是蟻周模型，故本文只介紹蟻周模型，其他兩種模型大家自行查閱下哈。規則是：完成一次路徑循環後，螞蟻才釋放信息素。有了這麽一個規則後，我們可以來想想，經過一次路徑循環後，路徑上的信息素為多少？
根據上面所提供的信息，我們可以知道，當所有的螞蟻完成一次路徑循環後，才更新信息素。因此路徑上的信息素應該分為兩部分：之前未揮發所殘留的信息素 + 經過當前循環所有螞蟻在經過該路徑後所留下的信息素。用公式表述如下：
$\tau_{ij}(t + n) = (1 - \rho ) \cdot \tau_{ij}(t) + \Delta \tau_{ij}(t, t + n )$
$\Delta \tau_{ij}(t, t + n) = \sum_{k = 1}^{m}\Delta \tau_{ij}^{k}(t, t + n)$
其中，
- $\rho$ -- 信息素揮發因子，$\rho \in [0, 1)$；
- $\Delta \tau_{ij}(t, t + n)$ -- 經過一次循環後城市i到城市j的路徑上的信息素的增量；
- $(t, t+n)$ -- 走過n步以後螞蟻即完成一次循環；
- $\Delta \tau_{ij}^{k}(t, t + n)$ -- 表示經過一次循環後螞蟻k在它走過的路上的信息素增量。
好了，現在我們的未揮發所殘留的信息素很好解，定義一個信息素揮發因子$\rho$便能解決。但是，經過一次循環所有螞蟻留下的信息素該怎麽定義？在蟻周模型中，它被定義如下：
$\Delta \tau_{ij}^{k}(t, t + n) = \left\{\begin{matrix} \frac{Q}{L_{k}} & Ant \ k \ walk \ though \ path \ (i,j)\\ 0 & other \end{matrix}\right.$
它表明，螞蟻留下的信息素跟它走過的完整路徑的總長度有關。越長則留下的信息素越少，這個應該很好理解。為了找到更短的路徑，就應該讓短的路徑信息素濃度高些。其中這個$Q$是控制比例常數，一般定為1。
3. 有了數學公式後，我們可以寫出蟻群算法的基本步驟：
(1) 參數初始化：
基本上，可以按照以下策略來進行參數組合設定：
確定螞蟻數目：螞蟻數目與城市規模之比約為1.5；
參數粗調：即調整取值範圍較大的$\alpha$，$\beta$及$Q$；
參數微調：即調整取值範圍較小的$\rho$；
% 參數初始化 NC_max = 200; % 最大叠代次數，取100~500之間 m = 50; % 螞蟻的個數，一般設為城市數量的1.5倍 alpha = 1; % α 選擇[1, 4]比較合適 beta = 5; % β 選擇[3 4 5]比較合適 rho = 0.1; % ρ 選擇[0.1, 0.2, 0.5]比較合適 Q = 1; NC = 1; % 叠代次數，一開始為1 Eta = 1 ./ D; % η = 1 / D(i, j) ,這裏是矩陣 Tau = ones(n, n); % Tau(i, j)表示邊(i, j)的信息素量，一開始都為1 Table = zeros(m, n); % 路徑記錄表 rBest = zeros(NC_max, n); % 記錄各代的最佳路線 lBest = inf .* ones(NC_max, 1); % 記錄各代的最佳路線的總長度 lAverage = zeros(NC_max, 1); % 記錄各代路線的平均長度
(2) 隨機將螞蟻放於不同出發點，對每只螞蟻計算其下個訪問城市，直到有螞蟻訪問完所有城市。
% 第1步：隨機產生各個螞蟻的起點城市
start = zeros(m, 1);
for i = 1: m
temp = randperm(n);
start(i) = temp(1);
end
Table(:, 1) = start; % Tabu表的第一列即是所有螞蟻的起點城市
citys_index = 1: n; % 所有城市索引的一個集合
% 第2步：逐個螞蟻路徑選擇
for i = 1: m
% 逐個城市路徑選擇
for j = 2: n
tabu = Table(i, 1: (j - 1)); % 螞蟻i已經訪問的城市集合（稱禁忌表）
allow_index = ~ismember(citys_index, tabu);
Allow = citys_index(allow_index); % Allow表：存放待訪問的城市
P = Allow;
% 計算從城市j到剩下未訪問的城市的轉移概率
for k = 1: size(Allow, 2) % 待訪問的城市數量
P(k) = Tau(tabu(end), Allow(k))^alpha * Eta(tabu(end), Allow(k))^beta;
end
P = P / sum(P); % 歸一化
% 輪盤賭法選擇下一個訪問城市（為了增加隨機性）
Pc = cumsum(P);
target_index = find(Pc >= rand);
target = Allow(target_index(1));
Table(i, j) = target;
end
end
(3) 計算各螞蟻經過的路徑長度$L_{k}$，記錄當前叠代次數最優解（即最短路徑和最短距離），同時對路徑上的信息素濃度進行更新
% 第3步：計算各個螞蟻的路徑距離
length = zeros(m, 1);
for i = 1: m
Route = Table(i, :);
for j = 1: (n - 1)
length(i) = length(i) + D(Route(j), Route(j + 1));
end
length(i) = length(i) + D(Route(n), Route(1));
end
% 第4步：計算最短路徑距離及平均距離
if NC == 1
[min_Length, min_index] = min(length);
lBest(NC) = min_Length;
lAverage(NC) = mean(length);
rBest(NC, :) = Table(min_index, :);
else
[min_Length, min_index] = min(length);
lBest(NC) = min(lBest(NC - 1), min_Length);
lAverage(NC) = mean(length);
if lBest(NC) == min_Length
rBest(NC, :) = Table(min_index, :);
else
rBest(NC, :) = rBest((NC - 1), :);
end
end
% 第5步：更新信息素
Delta_tau = zeros(n, n);
for i = 1: m
for j = 1: (n - 1)
Delta_tau(Table(i, j), Table(i, j + 1)) = Delta_tau(Table(i, j), Table(i, j + 1)) + Q / length(i);
end
Delta_tau(Table(i, n), Table(i, 1)) = Delta_tau(Table(i, n), Table(i, 1)) + Q / length(i);
end
Tau = (1 - rho) .* Tau + Delta_tau;
% 第6步：叠代次數加1，並且清空路徑記錄表
NC = NC + 1;
Table = zeros(m, n);
(4) 判斷是否達到最大叠代次數，若否，返回步驟2；是，結束程序。
(5) 輸出結果，並根據需要輸出尋優過程中的相關指標，如運行時間、收斂叠代次數等。
四、完整Matlab代碼
以下是利用蟻群算法計算TSP問題的代碼：
%% I. 清空環境
clc
clear all
%% II. 符號說明
% C -- n個城市的坐標
% NC_max -- 最大叠代次數
% m -- 蟻群中螞蟻的數量，一般設置為城市的1.5倍
% D(i, j) -- 兩城市i和之間的距離
% Eta(i, j) = 1 ./ D(i, j) -- 啟發函數
% alpha -- 表征信息素重要程度的參數
% beta -- 表征啟發函數重要程度的參數
% rho -- 信息素揮發因子
% Q --
% rBest -- 各代最佳的路線
% lBest -- 各代最佳路線的長度
% lAverage --各代的平均長度
%% III. 導入城市位置數據
citys = [16.4700 96.1000
16.4700 94.4400
20.0900 92.5400
22.3900 93.3700
25.2300 97.2400
22.0000 96.0500
20.4700 97.0200
17.2000 96.2900
16.3000 97.3800
14.0500 98.1200
16.5300 97.3800
21.5200 95.5900
19.4100 97.1300
20.0900 92.5500];
%% IV. 計算距離矩陣
D = Distance(citys); % 計算距離矩陣
n = size(D, 1); % 城市的個數
%% V. 初始化參數
NC_max = 200; % 最大叠代次數，取100~500之間
m = 22; % 螞蟻的個數，一般設為城市數量的1.5倍
alpha = 1; % α 選擇[1, 4]比較合適
beta = 5; % β 選擇[3 4 5]比較合適
rho = 0.1; % ρ 選擇[0.1, 0.2, 0.5]比較合適
Q = 1;
NC = 1; % 叠代次數，一開始為1
Eta = 1 ./ D; % η = 1 / D(i, j) ,這裏是矩陣
Tau = ones(n, n); % Tau(i, j)表示邊(i, j)的信息素量，一開始都為1
Table = zeros(m, n); % 路徑記錄表
rBest = zeros(NC_max, n); % 記錄各代的最佳路線
lBest = inf .* ones(NC_max, 1); % 記錄各代的最佳路線的總長度
lAverage = zeros(NC_max, 1); % 記錄各代路線的平均長度
%% VI. 叠代尋找最佳路徑
while NC <= NC_max
% 第1步：隨機產生各個螞蟻的起點城市
start = zeros(m, 1);
for i = 1: m
temp = randperm(n);
start(i) = temp(1);
end
Table(:, 1) = start; % Tabu表的第一列即是所有螞蟻的起點城市
citys_index = 1: n; % 所有城市索引的一個集合
% 第2步：逐個螞蟻路徑選擇
for i = 1: m
% 逐個城市路徑選擇
for j = 2: n
tabu = Table(i, 1: (j - 1)); % 螞蟻i已經訪問的城市集合（稱禁忌表）
allow_index = ~ismember(citys_index, tabu);
Allow = citys_index(allow_index); % Allow表：存放待訪問的城市
P = Allow;
% 計算從城市j到剩下未訪問的城市的轉移概率
for k = 1: size(Allow, 2) % 待訪問的城市數量
P(k) = Tau(tabu(end), Allow(k))^alpha * Eta(tabu(end), Allow(k))^beta;
end
P = P / sum(P); % 歸一化
% 輪盤賭法選擇下一個訪問城市（為了增加隨機性）
Pc = cumsum(P);
target_index = find(Pc >= rand);
target = Allow(target_index(1));
Table(i, j) = target;
end
end
% 第3步：計算各個螞蟻的路徑距離
length = zeros(m, 1);
for i = 1: m
Route = Table(i, :);
for j = 1: (n - 1)
length(i) = length(i) + D(Route(j), Route(j + 1));
end
length(i) = length(i) + D(Route(n), Route(1));
end
% 第4步：計算最短路徑距離及平均距離
if NC == 1
[min_Length, min_index] = min(length);
lBest(NC) = min_Length;
lAverage(NC) = mean(length);
rBest(NC, :) = Table(min_index, :);
else
[min_Length, min_index] = min(length);
lBest(NC) = min(lBest(NC - 1), min_Length);
lAverage(NC) = mean(length);
if lBest(NC) == min_Length
rBest(NC, :) = Table(min_index, :);
else
rBest(NC, :) = rBest((NC - 1), :);
end
end
% 第5步：更新信息素
Delta_tau = zeros(n, n);
for i = 1: m
for j = 1: (n - 1)
Delta_tau(Table(i, j), Table(i, j + 1)) = Delta_tau(Table(i, j), Table(i, j + 1)) + Q / length(i);
end
Delta_tau(Table(i, n), Table(i, 1)) = Delta_tau(Table(i, n), Table(i, 1)) + Q / length(i);
end
Tau = (1 - rho) .* Tau + Delta_tau;
% 第6步：叠代次數加1，並且清空路徑記錄表
NC = NC + 1;
Table = zeros(m, n);
end
%% VII. 結果顯示
[shortest_Length, shortest_index] = min(lBest);
shortest_Route = rBest(shortest_index, :);
disp([‘最短距離:‘ num2str(shortest_Length)]);
disp([‘最短路徑:‘ num2str([shortest_Route shortest_Route(1)])]);
%% VIII. 繪圖
figure(1)
plot([citys(shortest_Route,1); citys(shortest_Route(1),1)],...
[citys(shortest_Route,2); citys(shortest_Route(1),2)],‘o-‘);
grid on
for i = 1: size(citys, 1)
text(citys(i, 1), citys(i, 2), [‘ ‘ num2str(i)]);
end
text(citys(shortest_Route(1), 1), citys(shortest_Route(1), 2), ‘ 起點‘);
text(citys(shortest_Route(end), 1), citys(shortest_Route(end), 2), ‘ 終點‘);
xlabel(‘城市位置橫坐標‘)
ylabel(‘城市位置縱坐標‘)
title([‘蟻群算法優化路徑(最短距離: ‘ num2str(shortest_Length) ‘)‘])
figure(2)
plot(1: NC_max, lBest, ‘b‘, 1: NC_max, lAverage, ‘r:‘)
legend(‘最短距離‘, ‘平均距離‘)
xlabel(‘叠代次數‘)
ylabel(‘距離‘)
title(‘各代最短距離與平均距離對比‘)
其中，求兩城市之間的路徑距離被封裝到Distance.m中，如下：
function D = Distance(citys)
%% 計算兩兩城市之間的距離
% 輸入:各城市的位置坐標(citys)
% 輸出:兩兩城市之間的距離(D)
n = size(citys, 1);
D = zeros(n, n);
for i = 1: n
for j = i + 1: n
D(i, j) = sqrt((citys(i, 1) - citys(j, 1))^2 + (citys(i, 2) - citys(j, 2))^2);
D(j, i) = D(i, j);
end
D(i, i) = 1e-4; %對角線的值為0，但由於後面的啟發因子要取倒數，因此用一個很小數代替0
end
運行上述代碼，獲得的結果如下：
最短距離:29.6889 最短路徑:8 1 11 9 10 2 14 3 4 5 6 12 7 13 8
獲得的優化路徑圖如下：
數學建模方法-蟻群算法