本文轉自：http://www.kancloud.cn/digest/prandmethod/102851

        本節的目的是記錄以下學習和掌握模擬退火(Simulated Annealing，簡稱SA演算法)過程。模擬退火演算法是一種通用概率演算法，用來在一個大的搜尋空間內尋找命題的最優解。這裡分別使用隨機模擬退火演算法和確定性模擬退火演算法，在MATLAB平臺上進行程式設計，以尋找一個6-單元全連線網路的能量最小化模型。

參考書籍：Richard O. Duda, Peter E. Hart, David G. Stork 著《模式分類》

一、技術論述

1.隨機方法

學習在構造模式分類器中起著中心的作用。一個通常的做法是先假設一個單引數或多引數的模型，然後根據訓練樣本來估計各引數的取值。當模型相當簡單並且低維時，可以採用解析的方法，如求函式導數，來顯式求解方程以獲得最優引數。如果模型相對複雜一些，則可以通過計算區域性的導數而採用梯度下降演算法來解，如人工神經網路或其他一些最大似然方法。而對於更復雜的模型，經常會出現許多區域性極值，上述方法的效果往往不盡人意。

如果一個問題越複雜，或者先驗知識和訓練樣本越少，我們對能夠自動搜尋可行解的複雜搜尋演算法的依賴性就越強，如基於引數搜尋的隨即方法。一個通常的做法是使搜尋朝著預期最優解的區域前進，同時允許一定程度的隨機擾動，以發現更好的解。

2.隨機搜尋

這裡主要研究在多個候選解中搜索最優解的方法。假設給定多個變數si,i=1,2,…,N，其中每個變數的數值都取兩個離散值之一（如-1和1）。優化問題是這樣描述的：確定N個si的合適取值，時下述能量函式（又稱為代價函式）最小：

其中w_ij是一個實對稱矩陣，取值可正可負，其中令到自身的反饋權為零（即w_ii=0），這是因為非零w_ii只是在能量函式E上增加一個與si無關的常數，不影響問題的本質。這個優化問題可以用網路和節點的方式表示，如下圖所示，其中節點之間的連結對應每個權值w_ij。

3.貪心演算法的侷限性

如上所述，對於求解有N個變數si的能量E最小化問題，除非N的取值很小，否則往往無法直接求解，因為其構型數目高達N^2。如果使用貪心演算法搜尋最優的構型，需要先隨機選取每個節點的起始狀態，然後順序考查每個節點從而計算與之相聯絡的si=1狀態和si=-1狀態的能量，最後選取能夠降低能量的狀態遷移。這種判斷只用到了直接與之相連的具有非零權值的相鄰節點。

這種貪心搜尋演算法成功找到最優解的可能性很小，因為系統常常會陷入區域性能量最小值，或根本就不收斂，因此需要選擇其他搜尋方法。

4.模擬退火

在熱力學中，固體的退火過程主要由以下三部分組成：

1. 加溫過程。其目的是增強粒子的熱運動，使其偏離平衡位置。當溫度足夠高時，固體將熔解為液體，從而消除系統原先可能存在的非均勻態，使隨後進行的冷卻過程以某一平衡態為起點。熔解過程實際是系統的熵增過程，系統能量也隨溫度的升高而增大。
2. 等溫過程。物理學的知識告訴我們，對於與周圍環境交換熱量而溫度不變的封閉系統，系統狀態的自發變化總是朝自由能減少的方向進行，當自由能達到最小時，系統達到平衡態。
3. 冷卻過程。其目的是使粒子的熱運動減弱並漸趨有序，系統能量逐漸下降，從而得到低能的晶體結構。

Metropolis等在1953年提出了重要性取樣法，即以概率大小接受新狀態。具體而言，在溫度T時, 由當前狀態i產生新狀態j，兩者的能量分別為Ei和Ej，若Ej小於Ei，則接受新狀態j為當前狀態；否則，計算概率p(∆E)：

若p(∆E)大於[0,1]區間內的隨機數，則仍舊接受新狀態j為當前狀態；若不成立則保留i為當前狀態，其中k為玻爾茲曼常數，T為系統溫度。上述重要性取樣過程通常稱為Metropolis準則：

1. 在高溫下可接受與當前狀態能量差較大的新狀態；
2. 而在低溫下基本只接受與當前能量差較小的新狀態；
3. 且當溫度趨於零時，就不能接受比當前狀態能量高的新狀態。

1983年Kirkpatrick 等意識到組合優化與物理退火的相似性，並受到Metropolis 準則的啟迪，提出了模擬退火演算法。模擬退火演算法是基於Monte Carlo 迭代求解策略的一種隨機尋優演算法，其出發點是基於物理退火過程與組合優化之間的相似性，模擬退火方法由某一較高初溫開始，利用具有概率突跳特性的Metropolis抽樣策略在解空間中進行隨機搜尋，伴隨溫度的不斷下降，重複抽樣過程，最終得到問題的全域性最優。對比貪心演算法，模擬退火演算法主要的優勢在於它使系統有可能從區域性最小處跳出。

對於一個優化問題：

把優化問題的求解過程與統計熱力學的熱平衡問題進行類比，通過模擬高溫物體退火過程的方法，試圖找到優化問題的全域性最優或近似全域性最優解；

允許隨著引數的調整，目標函式可以偶爾向能量增加的方向發展（對應於能量有時上升），以利於跳出區域性極小的區域，隨著假想溫度的下降（對應於物體的退火），系統活動性降低，最終穩定在全域性最小所在的區域。

5.兩種模擬退火演算法

兩種模擬退火演算法，即隨機模擬退火和和確定性模擬退火演算法的實現步驟如下所示：

隨機模擬退火演算法收斂很慢，部分原因在於其中搜索的全部的構型空間的離散本質，即構型空間是一個N維超立方體。每一次搜尋軌跡都只能沿著超立方體的一條邊，狀態只能落在超立方體的頂點上，因此失去了完整的梯度資訊。而梯度資訊是可以用超立方體內部的連續狀態值提供的。一種更快的方法就是以下的確定性模擬退火演算法：

二、實驗結果討論

構造一個6-單元全連線網路，能量函式使用公式：

其中網路的連線權值矩陣如下：

設計步驟主要包括以下幾個部分： 
編寫程式[E, s_out] = RandomSimulatedAnnealing(T_max, time_max, c, s, w)，實現以上演算法1所述的隨機模擬退火演算法。這裡需要設定以下引數： T_max=10，T(m+1) =c*T(m)，c=0.9，進行實驗，能量隨溫度下降次數的變化曲線如圖2所示（由於模擬退火演算法所得到的結果有一定的隨機性，因此以下步驟均執行四次演算法進行觀察），四次所得到的最終構型s如圖3所示。

改變引數：初始溫度：T_max=5，T(m+1) =c*T(m)，c=0.5，進行實驗，能量隨溫度下降次數的變化曲線如圖4所示，四次所得到的最終構型s如圖5所示。

編寫程式[E, s_out] = DeterministicAnnealing(T_max, time_max, c, s, w)，實現以上演算法2所述的確定性模擬退火演算法。這裡需要設定以下引數： T_max=10，T(m+1) =c*T(m)，c=0.9，進行實驗，能量隨溫度下降次數的變化曲線如圖6所示，四次所得到的最終構型s如圖7所示。

改變引數：初始溫度：T_max=5，T(m+1) =c*T(m)，c=0.5，進行實驗，能量隨溫度下降次數的變化曲線如圖8所示，四次所得到的最終構型s如圖9所示。

結論：圖2、3給出了多次隨機模擬退火演算法的執行結果，可以看到構型s不一定完全一樣；能量函式E的波形在經過若干次逐漸遞減的震盪後基本收斂到全域性最小值-19。當改變T(1)=5,c=0.5時，從圖4中可觀察到能量函式E的波形極速下降，並達到較小的值，中間少了一個溫度漸變和震盪調整的過程。

圖6、7給出多次確定性模擬退火演算法的執行結果，且每次得到的最終構型s均一致；能量函式E的波形在經過平緩遞減後收斂到全域性最小值-19，不出現隨機模擬退火中劇烈震盪的情況。當改變T(1)=5,c=0.5時，從圖8中可觀察到能量函式E的波形同樣呈現出極速下降的態勢，並達到較小的值，中間少了一個溫度漸變和調整的過程。

綜合以上的實驗結果，我們發現隨機退火和確定性退火均能給出相似的最終解，但對於一些大規模的現實問題，隨機模擬退火的執行速度很慢，而相比之下確定性退火演算法要快很多，有時可以快2~3個數量級。

另外，初溫T(1)和溫度下降係數c的選擇對演算法效能也有很大影響。初溫的確定應折衷考慮優化質量和優化效率，常用方法包括以下幾種：

• 均勻抽樣一組狀態，以各狀態目標值的方差為初溫；
• 隨機產生一組狀態，確定兩兩狀態間的最大目標值差|∆max|，然後依據差值，利用一定的函式確定初溫。譬如T(1)=-∆/ln 
  p_r，其中p_r為初始接受概率；
• 利用經驗公式給出。

模擬退火演算法設計中包括三個重要的函式：狀態產生函式、狀態接受函式、溫度更新函式；同時在程式設計時，需遵循內迴圈終止準則、外迴圈終止準則。這些環節的設計將決定模擬退火演算法的優化效能。

三、實驗結果

四、簡單程式碼實現

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 隨機模擬退火函式
% 輸入引數：
%   T_max：初始溫度
%   time_max：最大迭代次數
%   c：溫度下降比率
%   s：初始構型
%   w：權值矩陣
% 輸出引數：
%   E：能量變化矩陣
%   s_out：經過演算法計算後的構型
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E, s_out] = RandomSimulatedAnnealing(T_max, time_max, c, s, w)

[x, y] = size(s);

time = 1;        % 迭代次數
T(time) = T_max; % 初始溫度設定

while (time < (time_max + 1)) % (T(time) > T_min) 
    for i = 1:1000
        num = ceil(rand(1) * y); % 選擇產生一個1到y之間的隨機數
        for j = 1:y
            e(j) = w(num, j) * s(num) * s(j);
        end
        Ea(time) = -1 / 2 * sum(e);
        Eb(time) = -Ea(time);
        if Eb(time) < Ea(time)
            s(num) = -s(num);
        elseif (exp(-(Eb(time) - Ea(time)) / T(time)) > rand())
            s(num) = -s(num);
        else
            s(num) = s(num);
        end
    end
    % 計算能量E
    E(time) = 0;
    for it = 1:6
        for jt = 1:6
            E(time) = E(time) + w(it, jt) * s(it) * s(jt);
        end
    end
    E(time) = E(time) * (-0.5);
    s_out(time,:) = s;  % 每次形成的構型
    time = time + 1;
    T(time) = T(time - 1) * c;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 確定性模擬退火函式
% 輸入引數：
%   T_max：初始溫度
%   time_max：最大迭代次數
%   c：溫度下降比率
%   s：初始構型
%   w：權值矩陣
% 中間函式：
%   tanh(l / T)：響應函式，該函式有一個隱含的重新規格化的作用
% 輸出引數：
%   E：能量變化矩陣
%   s_out：經過演算法計算後的構型
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E, s_out] = DeterministicAnnealing(T_max, time_max, c, s, w)

[x, y] = size(s);

time = 1;        % 迭代次數
T(time) = T_max; % 初始溫度設定

while (time < (time_max + 1))
    num = ceil(rand(1) * y); % 選擇產生一個1到y之間的隨機數
        for j = 1:y
            e(j) = w(num, j) * s(j);
        end
    l(time) = sum(e);
    s(num) = tanh(l(time) / T(time));
    % 計算能量E
    E(time) = 0;
    for it = 1:6
        for jt = 1:6
            E(time) = E(time) + w(it, jt) * s(it) * s(jt);
        end
    end
    E(time) = E(time) * (-0.5);
    s_out(time,:) = s; % 每次形成的構型
    time = time + 1;
    T(time) = T(time - 1) * c;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%模擬退火演算法實驗
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all;

% 網路的連線權值矩陣
w = [ 0 5 -3 4 4 1;...
      5 0 -1 2 -3 1;...
     -3 -1 0 2 2 0;...
      4 2 2 0 3 -3;...
      4 -3 2 3 0 5;...
      1 1 0 -3 5 0];

num = 6; % 總共產生6個數 
s_in = rand(1,num); % 生成1和-1的隨機序列 
s_in(s_in > 0.5) = 1; 
s_in(s_in < 0.5) = -1;
disp(['初始構型S為：',num2str(s_in)]);

% 以下是隨機模擬退火演算法
T_max = 10;        % 初始溫度設定
time_max = 100;    % 最大迭代次數
c = 0.9;           % 溫度變化比率
[E1, s_out1] = RandomSimulatedAnnealing(T_max, time_max, c, s_in, w);
subplot(221),plot(E1);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',隨機模擬退火演算法能量變化曲線']);
disp(['T(1) = 10,c = 0.9,隨機模擬退火演算法最終構型S為：',num2str(s_out1(time_max,:))]);

T_max = 5;         % 初始溫度設定
time_max = 100;    % 最大迭代次數
c = 0.5;           % 溫度變化比率
[E2, s_out2] = RandomSimulatedAnnealing(T_max, time_max, c, s_in, w);
subplot(222),plot(E2);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',隨機模擬退火演算法能量變化曲線']);
disp(['T(1) = 5,c = 0.5,隨機模擬退火演算法最終構型S為：',num2str(s_out2(time_max,:))]);

% 以下是確定性模擬退火演算法
T_max = 10;        % 初始溫度設定
time_max = 100;    % 最大迭代次數
c = 0.9;           % 溫度變化比率
[E3, s_out3] = DeterministicAnnealing(T_max, time_max, c, s_in, w);
subplot(223),plot(E3);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',確定性模擬退火演算法能量變化曲線']);
disp(['T(1) = 10,c = 0.9,確定性模擬退火演算法最終構型S為：',num2str(s_out3(time_max,:))]);

T_max = 5;         % 初始溫度設定
time_max = 100;    % 最大迭代次數
c = 0.5;           % 溫度變化比率
[E4, s_out4] = DeterministicAnnealing(T_max, time_max, c, s_in, w);
subplot(224),plot(E4);grid on;
title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',確定性模擬退火演算法能量變化曲線']);
disp(['T(1) = 5,c = 0.5,確定性模擬退火演算法最終構型S為：',num2str(s_out4(time_max,:))]);