數學建模方法-模擬退火算法
一、引言
哈嘍大家好，看到標題大家應該知道我今天要講什麽了吧。“模擬退火算法”，怎麽聽起來很燃的感覺，哈哈並沒有啦，一點都不燃，但是很有用！！看完這篇文章你就懂我什麽意思了。
二、退火現象
首先，我們了解一下什麽是“退火”。是指妖魔鬼怪快離開，火也快離開的意思嗎？差不多哈哈，因為火離開了溫度就低了嘛。好啦不開玩笑，講正經點，退火就是就是將金屬緩慢加熱到一定溫度，然後讓其緩慢的冷卻到室溫。比如說，博主有一碗面，可是溫度太高好燙吃不了，那就放在桌子上，等它慢慢降溫到合適的溫度就可以吃啦，這個過程就叫退火啦。
誒那我們這個模擬退火算法跟這個退火現象有什麽關系呢？其實還是有點關系的，模擬退火算法的思想，其實跟退火的現象有關。具體怎麽有關聽我慢慢描述哈。
三、模擬退火算法
首先看下面這個圖，假如我有一個函數$y = f(x)$，畫出來的圖就跟下圖那個一樣。現在，我想找到這個函數的最大值。那麽該怎麽做呢？模擬退火算法是這麽認為的：
我們先在$x$的定義域內，取一個起始點$x = i$，如圖紅色虛線，得到$y = f(i)$。接著，我們可以采取如下的策略：
讓紅色線往右移動一個單位，然後作如下判定：
(1) 如果 f(i + 1) > f(i) ，則接收該移動
這很好理解，如果你發現移動到的新點求出來的值更大，肯定是接收這個移動啦。
現在，我們的紅色線在不斷重復(1)操作的過程中，來到了函數的第一個極大值。如下圖：
現在，我再進行(1)步驟，就會發現：
f(i + 1) < f(i) ，那到了這裏我們是不是就不接受移動了。可是在這裏因為我們是知道這個圖是長什麽樣的，我們知道只要這個紅色線再經過“低谷”後就能找到一個真正的極大值了。可是電腦是不會知道圖的樣子的，那麽該怎麽做呢？
模擬退火算法的核心就在這裏，它不會放棄這個新的f(i + 1)，它說，f(i + 1)你可以比f(i)小，我給你機會，給你一段時間，看這個f(i + 1)能不能蛻變成“最大值”。因此對於f(i + 1) < f(i) ，執行如下操作：
(2) if f(i + 1) < f(i)， 以一定概率接收該移動，並隨時間推移概率降低。
而這個概率，就是參考了金屬冶煉的退火過程，這也正是這個算法被稱為模擬退火算法的原因。
四、Metropolis準則
在金屬冶煉的退火過程中，在溫度為T時，出現能量差為dE的粒子的降溫的概率為P(dE)，表示為：
$P \left( dE \right ) = exp\left({\frac{dE}{kT}}\right )$
其中k是一個常數，exp表示自然指數，且dE < 0。這條公式的含義是：溫度越高，出現一次能量差為dE的降溫的概率就越大；溫度越低，則出現降溫的概率就越小。又由於dE總是小於0（否則就不叫退火了），因此dE/kT < 0 ，所以P(dE)的函數取值範圍是(0,1) 。
在模擬退火的算法中，我們的一定概率是這樣定義的：
dC = f(i + 1) - f(i) < 0;
if exp(dC / T) >= rand，則接收該移動，否則不接收
上面這個定義稱為Metropolis準則。
因此。整個算法的流程是這樣的（偽代碼）
//目的：求系統的最小值
S(i): 系統在狀態y時的評價函數值
i： 系統當前狀態
i + 1: 系統的新狀態
rate: 控制降溫的速率
T0: 系統的初始溫度（高溫狀態）
T_min: 溫度的下限
while (T0 > T_min)
{
dE = S(i + 1) - S(i);
if dE < 0
//接收從S(i)到S(i+1)的移動
else if (exp(-dE / T) > random(0,1))
//接收從S(i)到S(i+1)的移動
else
//不接收從S(i)到S(i+1)的移動
T0 = r * T0; //降溫退火(0 < r < 1 ; r越大, 降溫越慢; r越小, 降溫越快)
i++;
}
五、模擬退火算法的應用
有這麽一個古老問題，被稱為旅行商問題 ( TSP , Traveling Salesman Problem ) ，是數學領域中著名的問題，問題內容是這樣的：有這麽一個旅行商人，他要拜訪n個城市，但是每次走的路徑是有限制的，即每個城市只能拜訪一次，並且最後要回到原來出發的城市。如何選擇路徑使得總路程為最小值。
之所以著名是因為這個問題至今還沒有找到一個有效的算法。如果想精確的求解只能只能通過窮舉所有的路徑組合但是隨著城市數量的增加其復雜度會很高。
不過還好，我們有模擬退火算法，當然它不能精確的求解這個問題，不過他能近似的求解這個問題。具體做法如下：
1. 設定初始溫度T0，並且隨機選擇一條遍歷路徑P(i)作為初始路徑，算出其長度L(P(i))；
2. 隨機產生一條新的遍歷路徑P(i+1)，算出其長度L(P(i + 1));
3. 若L(P(i + 1)) < L(P(i))，則接收P(i + 1)為新路徑，否則以模擬退火的概率接收P(i + 1)，然後降溫
4. 重復步驟1和2直至溫度達到最低值Tmin。
產生新的遍歷路徑的方法有很多，下面列舉其中3種：
1. 隨機選擇2個節點，交換路徑中的這2個節點的順序。
2. 隨機選擇2個節點，將路徑中這2個節點間的節點順序逆轉。
3. 隨機選擇3個節點m，n，k，然後將節點m與n間的節點移位到節點k後面。
六、模擬退火算法Matlab代碼
以下代碼是模擬退火算法用於TSP問題的。主函數如下：
%% I. 清空環境變量
clear all
clc
%% II. 導入城市位置數據
X = [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];
%% III. 計算距離矩陣
D = Distance(X); %計算距離矩陣
n = size(D,1); %城市的個數
%% IV. 初始化參數
T0 = 1e10; % 初始溫度
Tf = 1e-30; % 終止溫度
%L = 2; % 各溫度下的叠代次數
q = 0.9; % 降溫速率
Time = ceil(double(solve([num2str(T0) ‘*(0.9)^x = ‘, num2str(Tf)]))); % 計算叠代的次數 T0 * (0.9)^x = Tf
count = 0; % 初始化叠代計數
Obj = zeros(Time, 1); % 目標值矩陣初始化
path = zeros(Time, n); % 每代的最優路線矩陣初始化
%% V. 隨機產生一個初始路線
S1 = randperm(n);
DrawPath(S1, X) % 畫出初始路線
disp(‘初始種群中的一個隨機值:‘)
OutputPath(S1); % 用箭頭的形式表述初始路線
rlength = PathLength(D, S1); % 計算初始路線的總長度
disp([‘總距離：‘, num2str(rlength)]);
%% VI. 叠代優化
while T0 > Tf
count = count + 1; % 更新叠代次數
%temp = zeros(L,n+1);
% 1. 產生新解
S2 = NewAnswer(S1);
% 2. Metropolis法則判斷是否接受新解
[S1, R] = Metropolis(S1, S2, D, T0); % Metropolis 抽樣算法
% 3. 記錄每次叠代過程的最優路線
if count == 1 || R < Obj(count-1) % 如果當前溫度下的距離小於上個溫度的距離，記錄當前距離
Obj(count) = R;
else
Obj(count) = Obj(count-1);
end
path(count,:) = S1; % 記錄每次叠代的路線
T0 = q * T0; % 以q的速率降溫
end
%% VII. 優化過程叠代圖
figure
plot(1: count, Obj)
xlabel(‘叠代次數‘)
ylabel(‘距離‘)
title(‘優化過程‘)
grid on
%% VIII. 繪制最優路徑圖
DrawPath(path(end, :), X) % path矩陣的最後一行一定是最優路線
%% IX. 輸出最優解的路線和總距離
disp(‘最優解:‘)
S = path(end, :);
p = OutputPath(S);
disp([‘總距離：‘, num2str(PathLength(D, S))]);
其中為了更好的維護代碼和使代碼更加易讀，對一些功能進行了函數包裝，如下：
(1) 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
end
(2) DrawPath.m：將路徑用圖的形式表示出來
function DrawPath(Route, citys)
%% 畫路徑函數
% 輸入
% Route 待畫路徑
% citys 各城市坐標位置
%畫出初始路線
figure
plot([citys(Route, 1); citys(Route(1), 1)], [citys(Route, 2); citys(Route(1), 2)], ‘o-‘);
grid on
%給每個地點標上序號
for i = 1: size(citys, 1)
text(citys(i, 1), citys(i, 2), [‘ ‘ num2str(i)]);
end
text(citys(Route(1), 1), citys(Route(1), 2), ‘ 起點‘);
text(citys(Route(end), 1), citys(Route(end), 2), ‘ 終點‘);
(3) Metropolis.m：Metropolis準則判定
function [S,R] = Metropolis(S1, S2, D, T)
%% 輸入
% S1： 當前解
% S2: 新解
% D: 距離矩陣（兩兩城市的之間的距離）
% T: 當前溫度
%% 輸出
% S： 下一個當前解
% R： 下一個當前解的路線距離
R1 = PathLength(D, S1); % 計算S1路線長度
n = length(S1); % 城市的個數
R2 = PathLength(D,S2); % 計算S2路線長度
dC = R2 - R1; % 計算能力之差
if dC < 0 % 如果能力降低 接受新路線
S = S2; % 接受新解
R = R2;
elseif exp(-dC/T) >= rand % 以exp(-dC/T)概率接受新路線
S = S2;
R = R2;
else % 不接受新路線
S = S1;
R = R1;
end
end
(4) NewAnswer.m: 隨機產生新解的方式
function S2 = NewAnswer(S1) %% 輸入 % S1: 當前解 %% 輸出 % S2： 新解 n = length(S1); S2 = S1; a = round(rand(1, 2) * (n - 1) + 1); %產生兩個隨機位置，用於交換 W = S2(a(1)); S2(a(1)) = S1(a(2)); S2(a(2)) = W;
(5) OutputPath.m: 用箭頭的形式描繪出路徑方向
function p = OutputPath(Route)
%% 輸出路徑函數
% 輸入: 路徑(R)
%給R末端添加起點即R(1)
Route = [Route, Route(1)];
n = length(Route);
p = num2str(Route(1));
for i = 2: n
p = [p, ‘ —> ‘, num2str(Route(i))];
end
disp(p)
(6) PathLength.m: 算出總路徑的長度
function Length = PathLength(D, Route)
%% 計算起點與終點之間的路徑長度
% 輸入：
% D 兩兩城市之間的距離
% Route 路線
Length = 0;
n = length(Route);
for i = 1: (n - 1)
Length = Length + D(Route(i), Route(i + 1));
end
Length = Length + D(Route(n), Route(1));
運行了上述代碼後，得到的結果如下：
初始種群中的一個隨機值: 11 —> 10 —> 14 —> 9 —> 5 —> 3 —> 7 —> 4 —> 8 —> 1 —> 2 —> 13 —> 12 —> 6 —> 11 總距離：62.7207 最優解: 9 —> 11 —> 8 —> 13 —> 7 —> 12 —> 6 —> 5 —> 4 —> 3 —> 14 —> 2 —> 1 —> 10 —> 9 總距離：29.3405
也就是通過這個算法，它給我們找到的一個最短的距離是29.3405。我們不能保證這個答案一定是最短的，這就是這個算法的缺點了。如果你叠代的次數夠多，得到的結果就越接近正確的結果。但是叠代次數過多就會需要很長的計算時間，所以要適當選擇這個叠代次數。
最開始隨機選的路徑和算出來的最優解的一個路徑分別如下：
另外，隨著叠代次數的增加，總長度的變化圖如下：
可以看出，隨著我們不斷的叠代，距離是一直在不斷下降的。因此模擬退火算法，從某種意義上而言，在近似解決一些問題方面，還是起了很大的作用。。
數學建模方法-模擬退火算法