轉載自http://cighao.com/2015/12/03/introduction-of-SA/

1. 簡介

模擬退火演算法 (Simulated Annealing，SA) 最早的思想是由 N. Metropolis 等人於1953年提出。1983年, S. Kirkpatrick 等成功地將退火思想引入到組合優化領域。它是基於 Monte-Carlo 迭代求解策略的一種隨機尋優演算法，其出發點是基於物理中固體物質的退火過程與一般組合優化問題之間的相似性。模擬退火演算法從某一較高初溫出發，伴隨溫度引數的不斷下降,結合概率突跳特性在解空間中隨機尋找目標函式的全域性最優解，即在區域性最優解能概率性地跳出並最終趨於全域性最優。

2. 相關概念

退火

退火：將固體加熱到足夠高的溫度，使分子呈隨即排列狀態，然後逐步降溫使之冷卻，最後分子以低能狀態排列，固體達到某種穩定狀態。

退火過程：

加溫過程：增強粒子運動，消除系統原先可能存在的非均勻態。
等溫過程：對於與環境換熱而溫度不變的封閉系統，系統狀態的自發變化總是朝自由能減少的方向進行，當自由能達到最小時，系統達到平衡。
冷卻過程：使粒子熱運動減弱並漸趨有序，系統能量逐漸下降，從而得到低能的晶體結構。
Boltzman 概率分佈
Boltzman 概率分佈告訴我們：

同一個溫度，分子停留在能量小狀態的概率大於停留在能量大狀態的概率
溫度越高，不同能量狀態對應的概率相差越小，溫度足夠高時，各狀態對應概率基本相同。
隨著溫度的下降，能量最低狀態對應概率越來越大，溫度趨於0時，其狀態趨於1

Metropolis 準則
固體在恆定溫度下達到熱平衡的過程可以用 Metropolis 方法加以模擬。

溫度恆定時：若在溫度T，當前狀態i –> 新狀態j，若Ej小於Ei,則接受j為當前狀態；否則，若概率p=e^{-(Ej-Ei)/(k*T)} 大於[0,1)區間的隨機數，則仍接受狀態j為當前狀態；若不成立則保留狀態i為當前狀態。

溫度變化時：p=e^{-(Ej-Ei)/(k*T)} ,在高溫下，可接受當前狀態能量差較大的新狀態；在低溫下，只接受與當前狀態能量差較小的新狀態。

大白話精髓函式：針對模擬退火的這個精髓函式，其實它的最基本的原則就只有一個，就是如果交換後距離變小了，那麼一定接受這個交換結果，（重點來了）如果交換後距離變大了，那麼它不是立馬就拒絕這個結果，而是以一定的概率接受這個結果，這個概率的計算方法是 p=exp[-(Ej-Ei)/T],其中Ej相當於交換後的距離值，Ei相當於交換前的距離值，而T就是我們所說的溫度了，那麼我們說當Ej>Ei時，是不是以這個概率p接受這個解呀，想想，此時的(Ej-Ei)是不是一個小的正值呀（相當於誤差），而T溫度，肯定也是一個正值，只不過在迭代的過程中我們在一直減小這個T值，所以叫降溫嘛。好了我們來看看這個p值怎麼變化的吧，再想想，我們說開始的時候高溫，也就是T很大是不是，那麼exp[-(Ej-Ei)/T]想想大概是多少，假設(Ej-Ei)也就是誤差波動不大的話，是不是應為T很大導致-(Ej-Ei)/T從負向趨近於0呀，那麼exp[-(Ej-Ei)/T]趨近於1是吧（exp函式性質，不懂的話就沒辦法了），也就是這個概率p趨近於1，什麼概念呀，是不是基本上溫度T高的時候，即使碰到了交換後不好的結果也接受呀。好了，再想想，隨著T迭代慢慢的下降，T作為分母是不是使得整個P值在慢慢下降呀，也就是接受壞的結果的概率慢慢下降，到最後當T下降到很小的時候，-(Ej-Ei)/T是不是就負向很大了呀，那麼p=-(Ej-Ei)/T基本上趨近於0了，也就是再碰到壞的解就基本上不接受了吧。理解了這個精髓部分，這個演算法基本上差不多了，剩下的就是發揮部分了，怎麼去選擇初始溫度呀，迭代次數是多少合適呀，每次交換幾個點的位置比較好呀，再有就是如何降溫，是以一定的比例降溫還是怎麼降溫呀，這一切一切都決定著你的演算法是否能很好的找到最優解

3. 原理描述

模擬退火演算法來源於固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時，固體內部粒子隨溫升變為無序狀，內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最後在常溫時達到基態，內能減為最小。根據 Metropolis 準則，粒子在溫度T時趨於平衡的概率為e(-ΔE/(kT))，其中E為溫度T時的內能，ΔE為其改變數，k為 Boltzmann 常數。用固體退火模擬組合優化問題，將內能E模擬為目標函式值f，溫度T演化成控制引數t，即得到解組合優化問題的模擬退火演算法：由初始解i和控制引數初值t開始，對當前解重複“產生新解→計算目標函式差→接受或捨棄”的迭代，並逐步衰減t值，演算法終止時的當前解即為所得近似最優解，這是基於蒙特卡羅迭代求解法的一種啟發式隨機搜尋過程。退火過程由冷卻進度表(Cooling Schedule)控制，包括控制引數的初值t及其衰減因子Δt、每個t值時的迭代次數L和停止條件S。

4.演算法步驟

step1
由一個產生函式從當前解產生一個位於解空間的新解；為便於後續的計算和接受，減少演算法耗時，通常選擇由當前新解經過簡單地變換即可產生新解的方法，如對構成新解的全部或部分元素進行置換、互換等，注意到產生新解的變換方法決定了當前新解的鄰域結構，因而對冷卻進度表的選取有一定的影響。

step2
計算與新解所對應的目標函式差。因為目標函式差僅由變換部分產生，所以目標函式差的計算最好按增量計算。事實表明，對大多數應用而言，這是計算目標函式差的最快方法。

step3
判斷新解是否被接受,判斷的依據是一個接受準則，最常用的接受準則是 Metropolis 準則: 若Δt′<0則接受S′作為新的當前解S，否則以概率exp(-Δt′/T)接受S′作為新的當前解S。

step4
當新解被確定接受時，用新解代替當前解，這隻需將當前解中對應於產生新解時的變換部分予以實現，同時修正目標函式值即可。此時，當前解實現了一次迭代。可在此基礎上開始下一輪試驗。而當新解被判定為捨棄時，則在原當前解的基礎上繼續下一輪試驗。

程式碼

function [xo,fo] = Opt_Simu(f,x0,u,l,kmax,q,TolFun)
% 模擬退火演算法求函式 f(x)的最小值點， 且 l <= x <= u
% f為待求函式，x0為初值點，l，u分別為搜尋區間的上下限，kmax為最大迭代次數
% q為退火因子，TolFun為函式容許誤差
%%%%演算法第一步根據輸入變數數，將某些量設為預設值
if nargin < 7
    TolFun = 1e-8;
end
if nargin < 6
    q = 1;
end
if nargin < 5
    kmax = 100;
end
%%%%演算法第二步，求解一些基本變數
N = length(x0); %自變數維數
x = x0;
fx = feval(f,x); %函式在初始點x0處的函式值
xo = x;
fo = fx;
%%%%%演算法第三步，進行迭代計算，找出近似全域性最小點
for k =0:kmax
    Ti = (k/kmax)^q;
    mu = 10^(Ti*100); % 計算mu
    dx = Mu_Inv(2*rand(size(x))-1,mu).*(u - l);%步長dx
    x1 = x + dx; %下一個估計點
    x1 = (x1 < l).*l +(l <= x1).*(x1 <= u).*x1 +(u < x1).*u; %將x1限定在區間[l,u]上
    fx1 = feval(f,x1);
    df = fx1- fx;
    if df < 0||rand < exp(-Ti*df/(abs(fx) + eps)/TolFun) %如果fx1<fx或者概率大於隨機數z
        x = x1;
        fx = fx1;
    end
    if fx < fo
        xo = x;
        fo = fx1;
    end
end
function x = Mu_Inv(y,mu)
x = (((1+mu).^abs(y)- 1)/mu).*sign(y);