模擬退火(Simulated Annealing, SA)演算法簡介與MATLAB實現
目錄
模擬退火演算法概述
-
模擬退火演算法(Simulated Annealing，簡稱SA)的思想最早是由Metropolis等提出的。其出發點是基於物理中固體物質的退火過程與一般的組合優化問題之間的相似性。模擬退火法是一種通用的優化演算法，其物理退火過程由以下三部分組成:
-
加溫過程。其目的是增強粒子的熱運動，使其偏離平衡位置。當溫度足夠高時，固體將熔為液體，從而消除系統原先存在的非均勻狀態。
-
等溫過程。對於與周圍環境交換熱量而溫度不變的封閉系統，系統狀態的自發變化總是朝自由能減少的方向進行的，當自由能達到最小時，系統達到平衡狀態
。 -
冷卻過程。使粒子熱運動減弱，系統能量下降，得到晶體結構。
-
-
加溫過程相當於對演算法設定初值，等溫過程對應演算法的Metropolis抽樣過程，冷卻過程對應控制引數的下降。這裡能量的變化就是目標函式，我們要得到的最優解就是能量最低態。其中Metropolis準則是SA演算法收斂於全域性最優解的關鍵所在，Metropolis準則以一定的概率接受惡化解，這樣就使演算法跳離區域性最優的陷阱。
-
SA演算法的Metropolis準則允許接受一定的惡化解，具體來講，是以一定概率來接受非最優解。舉個例子，相當於保留一些“潛力股”，使解空間裡有更多的可能性。對比輪盤賭法，從概率論來講，它是對非最優解給予概率0，即全部拋棄。
-
模擬退火本身是求一個最小值問題，但可以轉化為求最大值問題，只需要對目標函式加個負號或者取倒數。
演算法步驟
-
初始化:取初始溫度T0足夠大，令T = T0，任取初始解S1。
對比GA、ACA、PSO之類的群優化演算法，需要在解空間中產生多個個體，再繼續尋優。而SA演算法只需要一個點即可。 -
對當前溫度T，重複第(3)~(6)步。
-
對當前解S1隨機擾動產生一個新解S2。
此處隨機的擾動沒有定義。結合實際例子做選擇。 -
計算S2的增量df = f(S2) - f(S1)，其中f(S1)為S1的代價函式。
代價函式相當於之前群優化演算法中講的適應度函式。 -
若df < 0，則接受S2作為新的當前解，即S1 = S2；否則，計算S2的接受概率exp(-df/T)， T是溫度。隨機產生(0,1)區間上均勻分佈的隨機數rand，若exp(-df/T) > rand，也接受S2作為新的當前解S1 = S2，否則保留當前解S1。
這是SA演算法的核心部分，即有一定概率接受非最優解。 -
如果滿足終止條件Stop，則輸出當前解S1為最優解，結束程式，終止條件Stop通常取為在連續若干個Metropolis鏈中新解S2都沒有被接受時終止演算法或者是設定結束溫度。否則按衰減函式衰減T後返回第(2)步，即被接受的新的解一直在產生，則我們要對問題進行降溫，使得非最優解被接受的可能不斷降低，結果愈發收斂於最優解。
演算法特點
• 與遺傳演算法、粒子群優化演算法和蟻群演算法等不同，模擬退火演算法不屬於群優化演算法，不需要初始化種群操作。
• 收斂速度較慢。因為1）它初始溫度一般設定得很高，而終止溫度設定得低，這樣才符合物體規律，認為物質處於最低能量平衡點；2）它接受惡化解，並不是全程都在收斂的過程中。這一點可以類比GA中的變異，使得它不是持續在收斂的，所以耗時更多一些。
• 溫度管理（起始、終止溫度）、退火速度（衰減函式）等對尋優結果均有影響。比如T的衰減速度如果太快，就會導致可能尋找不到全域性最優解。
模擬退火演算法MATLAB實現
MATLAB自帶模擬退火演算法工具箱。在本文中，我們採用自帶工具箱來解決函式優化的問題，再使用自己編寫的程式來解決一個TSP問題。
要使用自帶的工具箱，我們先安裝：
【例1】一元/多元函式優化
我們要找：
一元函式：x = [1,2]範圍內 y = sin(10*pi*x) / x 的極值
二元函式：在x,y都是[-5,5]範圍內找z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20 的極值
上面是我們提前獲得的ground-truth（用main.m，程式碼在下文）。下面我們假裝不知道這個結果的，用模擬退火方法來搜尋：
首先定義我們在SA演算法中需要用到的代價函式fitness.m：
function fitnessVal = fitness( x )
%一元函式優化：
fitnessVal = sin(10*pi*x) / x; %求最小值
% fitnessVal = -1 * sin(10*pi*x) / x; 用模擬退火求最大值，可以加個負號或者弄個倒數！
%二元函式優化：
% fitnessVal = -1 * (x(1)^2 + x(2).^2 - 10*cos(2*pi*x(1)) - 10*cos(2*pi*x(2)) + 20);
end
主程式main.m，用於我們直觀看圖先獲得ground-truth：
%% I. 清空環境變數
clear all
clc
%% II. 一元函式優化
x = 1:0.01:2;
y = sin(10*pi*x) ./ x;
figure
plot(x,y,'linewidth',1.5)
ylim([-1.5, 1.5])
xlabel('x')
ylabel('y')
title('y = sin(10*pi*x) / x')
hold on
%%
% 1. 標記出最大值點
[maxVal,maxIndex] = max(y);
plot(x(maxIndex), maxVal, 'r*','linewidth',2)
text(x(maxIndex), maxVal, {[' X: ' num2str(x(maxIndex))];[' Y: ' num2str(maxVal)]})
hold on
%%
% 2. 標記出最小值點
[minVal,minIndex] = min(y);
plot(x(minIndex), minVal, 'ks','linewidth',2)
text(x(minIndex), minVal, {[' X: ' num2str(x(minIndex))];[' Y: ' num2str(minVal)]})
%% III. 二元函式優化
[x,y] = meshgrid(-5:0.1:5,-5:0.1:5);
z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20;
figure
mesh(x,y,z)
hold on
xlabel('x')
ylabel('y')
zlabel('z')
title('z = x^2 + y^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20')
%%
% 1. 標記出最大值點
maxVal = max(z(:));
[maxIndexX,maxIndexY] = find(z == maxVal);
for i = 1:length(maxIndexX)
plot3(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)), maxVal, 'r*','linewidth',2)
text(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)), maxVal, {[' X: ' num2str(x(maxIndexX(i),maxIndexY(i)))];[' Y: ' num2str(y(maxIndexX(i),maxIndexY(i)))];[' Z: ' num2str(maxVal)]})
hold on
end
現在正式開始模擬退火演算法環節。在MATLAB命令中輸入：optimtool，開啟工具箱。
boltzmann就是對應Metropolis準則的退火方式。
選擇視覺化的輸出的專案：
對於二元函式，在工具箱中的設定方法大同小異：記得先在fitness.m中修改我們的目標函式！
執行後即可得到結果。
【例2】TSP問題
這裡和上一篇講蟻群演算法的博文不同，我們的TSP中城市是自己虛擬的14座城市。
main.m:
%% 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; % 初始溫度，10的10次方！需要設定一個很大的溫度。
Tend = 1e-30; % 終止溫度
L = 2; % 各溫度下的迭代次數
q = 0.9; %降溫速率
Time = ceil(double(solve([num2str(T0) '*(0.9)^x = ',num2str(Tend)]))); % 計算迭代的次數
% Time = 132;
count = 0; %迭代計數
Obj = zeros(Time,1); %目標值矩陣初始化
track = zeros(Time,N); %每代的最優路線矩陣初始化
%% V. 隨機產生一個初始路線
S1 = randperm(N);
DrawPath(S1,X)
disp('初始種群中的一個隨機值:')
OutputPath(S1);
Rlength = PathLength(D,S1);
disp(['總距離：',num2str(Rlength)]);
%% VI. 迭代優化
while T0 > Tend
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
track(count,:) = S1;
T0 = q * T0; %降溫
end
%% VII. 優化過程迭代圖
figure
plot(1:count,Obj)
xlabel('迭代次數')
ylabel('距離')
title('優化過程')
%% VIII. 繪製最優路徑圖
DrawPath(track(end,:),X)
%% IX. 輸出最優解的路線和總距離
disp('最優解:')
S = track(end,:);
p = OutputPath(S);
disp(['總距離：',num2str(PathLength(D,S))]);
計算距離的函式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(sum((citys(i,:) - citys(j,:)).^2));
D(j,i) = D(i,j);
end
end
畫出路徑的函式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),' 終點');
輸出路徑函式OutputPath.m：
function p = OutputPath(R)
%% 輸出路徑函式
% 輸入：R 路徑
R = [R,R(1)];
N = length(R);
p = num2str(R(1));
for i = 2:N
p = [p,'―>',num2str(R(i))];
end
disp(p)增加隨機擾動產生新解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)) = S2(a(2));
S2(a(2)) = W; %得到一個新路線我們的做法是隨機產生兩個城市讓他們交換位置，從而得到一個新的路徑。當然，這只是這個問題的一個做法，也有其他“增加隨機擾動”的做法，而且對於多元函式問題更加簡單，只要在當前解的附近增加一些小的值即可。
Metropolis準則的實現：
function [S,R] = Metropolis(S1,S2,D,T)
%% 輸入
% S1： 當前解
% S2: 新解
% D: 距離矩陣（兩兩城市的之間的距離）
% T: 當前溫度
%% 輸出
% S： 下一個當前解
% R： 下一個當前解的路線距離
R1 = PathLength(D,S1); %計算路線長度
N = length(S1); %得到城市的個數
R2 = PathLength(D,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程式結果：
