模擬退火的實作思路、搜索模型及其可視化

原理與特點

模擬退火(Simulated Annealing) 是一種啟發式的優化演算法，可以針對較復雜的決議式實作快速搜索其最值，優化目標函式，演算法的實作程序與金屬的退火程序類似，

退火 是指將金屬加熱到一定溫度直到液化再冷卻下來的程序，在退火程序中原子按照區域能量最小化進行排列，在冶金學中常利用金屬的退火進行原子重排，以得到韌性和硬度都更好的金屬，¹

模擬退火實質上是一種隨機搜索演算法，其核心在于溫度這一變數，在退火初期溫度較高時，目標函式的狀態變數通過隨機無序運動盡可能探索狀態空間，尋找目標函式更優值，并一定概率接受更差狀態以便于跳出區域最優解，而隨著溫度降低，探索幅度下降，狀態變數趨于穩定，

模擬退火常常用于解當演算法無法遍歷整個狀態空間去尋找目標函式的最優解時的問題，也就是組合爆炸(combination explosion) 的情況，這個時候可以采用貪心演算法去快速尋找一個近似解，但同樣也會有概率陷入區域最優解，模擬退火演算法一方面能夠通過隨機搜索的方式節省算力，為NP復雜性問題提供有效的近似解，另一方面，也能實作通過溫度這一變數實作盡可能廣的探索，使得收斂結果能夠跳出區域最優解，‘

演算法實作思路

1.冷卻程序

在模擬退火中，溫度作為一個全域變數使用，首先需要明確溫度下降函式，（溫度T減小的模型也稱為冷卻進度表）

溫度的下降一般有幾何冷卻和線性冷卻兩種方式：

線性冷卻

T ′ = T ? α (1) T'=T-\alpha \tag{1} T′=T?α(1)

幾何冷卻

T ′ = T α (2) T'=T\alpha \tag{2} T′=Tα(2)

其中 α \alpha α為0到1之間的數
一般在模擬退火演算法中采用幾何冷卻的方式模擬溫度的下降

退火代碼總體框架

double T0=5000; //初始溫度
double alpha=0.01; //每次退火1%  
while(T>1e-4)
{
            //退火程序 Process of annealing

T=T*（1-alpha）; //冷卻程序 Process of cooling
}

2.退火程序

（1）選擇起點

也就是初始化自變數 x 0 x_0 x0?

（2）與溫度正相關的隨機搜索

隨機搜索的幅度一般選擇與溫度線性相關，用公式表示如下

d x = R a n d ? T (3) dx=Rand*T \tag{3} dx=Rand?T(3)
x n e w = x o l d + d x (4) x_{new}=x_{old}+dx \tag{4} xnew?=xold?+dx(4)

其中 d x dx dx表示隨機搜索的幅度， R a n d Rand Rand為使得 x n e w x_{new} xnew?生成在定義域范圍內滿足平均分布的亂數,可以設定為未生成在定義域內則重新生成，但演算法往往會在生成規范的亂數上花一定的時間，降低了退火程序的效率，也可以通過設計搜索模型實作一次性生成滿足要求的亂數，

隨機搜索是模擬退火的關鍵，決定了最終解的精確度，其核心是: 滿足在退火的前期隨機搜索的程度較大使得盡可能搜索狀態空間，優化目標函式并跳出區域最優解，而在退火的后期應該減小隨機搜索程度使得自變數在最優解附近的空間搜索，使得解更加精確，

（3）接受更好的狀態或以一定概率接受更壞的狀態

退火演算法中，跳出區域最優解的核心在于利用Metropolis抽樣法確定對新的狀態的選擇準則，當系統能量狀態將要變化時，系統能量由E(n)變化到E(n+1)的概率為
r = 1 , E ( n + 1 ) < E ( n ) (5) r=1, E(n+1)<E(n) \tag{5} r=1,E(n+1)<E(n)(5)

r = e ? E ( n ) ? E ( n + 1 ) T , E ( n + 1 ) > E ( n ) (6) r=e^{-\frac{E(n)-E(n+1)}{T}} ,E(n+1)>E(n) \tag{6} r=e?TE(n)?E(n+1)?,E(n+1)>E(n)(6)

也就是系統狀態向更有利的方向（能量更小）發展時，一定接受新的狀態；如果向不利方向發展，則以一定概率( e ? E ( n ) ? E ( n + 1 ) T e^{-\frac{E(n)-E(n+1)}{T}} e?TE(n)?E(n+1)?)接受，
（E(n)可以看作就是目標函式f(x)，E(n)越小，代表f(x)更優，對應的x更優)

可以看出，溫度影響了對新的狀態接受概率，當溫度較高時，可以接受較差情況的解，而當溫度降低時，接收更差情況的解概率更小，以便于穩定當前解，

實體分析

例如利用MATLAB求解函式 f ( x ) = s i n ( x ) + 6 s i n ( 2 x ) f(x)=sin(x)+6sin(2x) f(x)=sin(x)+6sin(2x) 在（-3，3）之間的最大值，

當然在數學建模應用中，目標函式的實際決議式會比較復雜，難以通過遍歷的方式求得最值時，可以用模擬退火得到解的近似值，

以下的代碼均為MATLAB代碼格式

（1）初始化設定

%目標函式： 在（-3，3）區間找到f(x)最大值對應的x
f=@(x)(sin(x)+3*sin(2*x)); 

%退火初始化設定
T0=12000;%初始溫度
T=T0;
dT=0.01;%每次退火當前溫度的1%
x=0;%退火起點
x_final=x;%最優解   精確解為x=0.8411

%亂數初始化設定
rng shuffle; %為每次計算設定不同的亂數種子

(2)與溫度正相關的隨機搜索

%退火程序
while(T>0.0001）
    %以溫度的正相關概率隨機移動（難點）
    dx=(((-3-x)/T0)+(6/T0)*rand())*(T0+0.8*(T-T0)*exp(-T));  
    %在((-3-x)/T0，(3-x)/T0)范圍內的隨機移動

這里有幾種方法設計dx=f(T)的函式，此代碼中所設計的函式為：
d x = [ ? 3 ? x T 0 + 6 T 0 r a n d ( ) ] ? [ T 0 + K ( T ? T 0 ) e ? T ] (7) dx= [\frac{-3-x}{T_0}+\frac{6}{T_0}rand()]*[T_0+K(T-T_0)e^{-T}] \tag{7} dx=[T0??3?x?+T0?6?rand()]?[T0?+K(T?T0?)e?T](7)

T₀ 為初始溫度，x 為當前的解，rand() 為(0,1)之間的亂數，T 為當前溫度，K 為常數，

注：生成范圍在（a,b）之間的亂數可以表示為 a + ( b ? a ) ? r a n d ( ) a+(b-a)*rand() a+(b?a)?rand()表示，例如 ? 3 ? x T 0 + 6 T 0 r a n d ( ) \frac{-3-x}{T_0}+\frac{6}{T_0}rand() T0??3?x?+T0?6?rand()表示生成在（ ? 3 ? x T 0 , 3 ? x T 0 \frac{-3-x}{T_0},\frac{3-x}{T_0} T0??3?x?,T0?3?x?）之間的亂數，通過此方法生成的 d x dx dx有 x + d x x+dx x+dx滿足定義域（-3，3）范圍

此搜索模型可以一次性生成在定義域范圍內的亂數，提高退火效率，也能控制搜索的衰減程序，
引數K可以控制無序運動的衰減程度，設定合適的引數K可以更有效靈活地控制退火的搜索程序，使得在退火前期充分搜索退火后期處于穩定狀態，注意此處K的取值為（0，1）

當然也可以利用重復生成與溫度正相關的亂數，直到生成定義域范圍的方法，代碼如下：

  dx=inf;
  while(dx+x>b||dx+x<a) %當新生成的點在定義域外時
        dx=(2*rand()-1)*T; %重新生成亂數直到滿足要求
  end

（3）接受更好的狀態或以一定概率接受更壞的狀態

如果隨機搜索得到的新的自變數對應的目標函式 f ( x n e w ) f(x_{new}) f(xnew?)狀態更優，則接受這個自變數作為當前解( x f i n a l x_{final} xfinal?)，如果新的狀態 f ( x n e w ) f(x_{new}) f(xnew?)更差則以一定概率接受，

%接受更好的狀態
    if(f(x)>=f(x_final))
        x_final=x;  
    %以一定概率接受更差的狀態
    else
        df=f(x)-f(x_final); %計算df 即兩者的狀態差距大小
        if(rand()<exp(df/T))%（0，1）的亂數小于exp(df/T)的概率為exp(df/T)
            x_final=x;
        end
    end

（4）降溫

    T=T*(1-dT);
    end

完整代碼

%模擬退火
clear all

%目標函式： 在（-3，3）區間找到f(x)最大值對應的x
f=@(x)(sin(x)+3*sin(2*x)); 

%退火初始化設定
T0=12000;%初始溫度
T=T0;
dT=0.01;%每次退火1%
x=0;%退火起點
x_final=x;%最優解   精確解為x=0.8411

%亂數初始化設定
rng shuffle; %為每次計算設定不同的亂數種子

%退火程序
while(T>0.0001)

    %以溫度的正相關概率隨機移動（難點）
    dx=(((-3-x)/T0)+(6/T0)*rand())*(T0+0.8*(T-T0)*exp(-T));  %在((-3-x)/T0，(3-x)/T0)范圍內的隨機移動
    x=x+dx;
    
    %接受更好的狀態
    if(f(x)>=f(x_final))
        x_final=x;  
    %以一定概率接受更差的狀態
    else
        df=f(x)-f(x_final); %計算df 即兩者的狀態差距大小
        if(rand()<exp(df/T))%（0，1）的亂數小于exp(df/T)的概率為exp(df/T)
            x_final=x;
        end
    end   
    
    %退火
    T=T*(1-dT);
end

求解結果

進行十次求解，得到的結果如下表所示，
注：精確解為x=0.8411

實作較準確的結果

（1）初始溫度與終止溫度

通過設定較高的初始溫度 T 0 T_0 T0?和盡可能接近0的終止溫度提高搜索的精度，

（2）溫度下降模型

一般采用幾何冷卻的模型，也就是溫度T按照指數函式規律下降，使得在搜索后期能夠在精確解附近進行搜索，

（3）搜索模型

關于隨機搜索 d x = f ( T ) dx=f(T) dx=f(T)有幾種模型可供參考：

3.1衰減可控模型

d x = [ a ? x T 0 + b ? a T 0 r a n d ( ) ] ? [ T 0 + K ( T ? T 0 ) e ? T ] (8) dx= [\frac{a-x}{T_0}+\frac{b-a}{T_0}rand()]*[T_0+K(T-T_0)e^{-T}] \tag{8} dx=[T0?a?x?+T0?b?a?rand()]?[T0?+K(T?T0?)e?T](8)
T 0 T_0 T0? 為初始溫度， x x x為當前的解， r a n d ( ) rand() rand()為(0,1)之間的亂數，模型滿足 d x + x dx+x dx+x（新的搜索點）的范圍為[a,b]， T T T為當前溫度， K K K 為常數，

可以通過控制引數 K K K實作對搜索范圍衰減的控制，以下是在初始溫度 T 0 = 12000 T_0=12000 T0?=12000、衰減率 d T = 0.01 dT=0.01 dT=0.01，終止溫度為0.001時的測驗結果，

1）當K=0.98時，隨機運動的衰減程度較高，容易過早地收斂從而無法找到更精確的全域解，

2）當K=0.8時，能夠實作退火前期充分搜索，后期在精確解附近搜索，從而準確收斂，

3）當K=0.5時，搜索程序衰減較小，能夠實作對自變數狀態空間較充分的搜索，不會陷入區域最優解，到了退火后期，當前解在精確解附近，但由于搜索幅度仍然較大，導致有一定幾率無法搜索到精確解附近的狀態空間，

實體求解中采用此模型

2.線性相關模型

d x = T ? R a n d (9) dx=T*Rand \tag{9} dx=T?Rand(9)

采用重復生成的方法，直到生成滿足定義域的亂數，大小與溫度線性相關，由于溫度是指數衰減，也就是實際上無序運動按照指數規律衰減，

  dx=inf;
  while(dx+x>b||dx+x<a) %當新生成的點在定義域外時
        dx=(2*rand()-1)*T; %重新生成亂數直到滿足要求
  end

同樣此模型會一定程度降低演算法效率
但如果采用初始溫度標幺模型：

d x = [ a ? x T 0 + b ? a T 0 r a n d ( ) ] ? T (10) dx= [\frac{a-x}{T_0}+\frac{b-a}{T_0}rand()]*T \tag{10} dx=[T0?a?x?+T0?b?a?rand()]?T(10)

雖然能夠使得生成的亂數一次性滿足定義域 （ a , b ） （a,b） （a,b）,但通過收斂表可以看到搜索程序隨著溫度衰減而過早地衰減，以至于還沒有搜索到精確解附近的狀態空間時就已經陷入了區域最優解，因此不推薦此類模型，

模擬退火可視化

實體：求解 f ( x ) = s i n ( x ) + 6 s i n ( 2 x ) f(x)=sin(x)+6sin(2x) f(x)=sin(x)+6sin(2x)在（-3，3）之間的最大值

求解程序如圖所示：

可視化之前，需要對每個搜索程序的變數保存到作業區，最后統一畫圖，在畫圖腳本里，可以使用matlab命令hold on和hold off，使得每一幀影像只包含函式影像和當前解，之前的解不再出現，然后通過儲存每一幀影像合成動圖gif格式，具體可參考命令：help(imwrite)，

可視化代碼：
注：x_final_temp變數儲存了每一次的當前解

%畫動圖
    figure(4)
    i=-3:0.01:3;
    
    %畫動態影像
    for k=1:count-1
    
    %畫整個函式影像
    plot(i,f(i));
    hold on        
    grid on 
    title(['溫度Temperature=',num2str(T_temp(k)),"℃"])%題目
    
    %用.標記出當前解，大小為15
    plot(x_final_temp(k),f(x_final_temp(k)),'.','MarkerSize',15);
    legend('f(x)','當前解') %圖注
    hold off
    
    %儲存每一幀影像到檔案
    frame = getframe(figure(4));
    im{k} = frame2im(frame);
    filename = 'SA.gif'; 
    [A,map] = rgb2ind(im{k},256);
    if k == 1
        imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
    else
        imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
    end
    end 

當然，目標函式的自變數也可以拓展到二維，例如將 f ( x ) = s i n ( x ) + 6 s i n ( 2 x ) f(x)=sin(x)+6sin(2x) f(x)=sin(x)+6sin(2x)繞Z軸旋轉得到曲面，模擬退火求解程序類似，不再敘述，由于gif檔案過大，這里不給出可視化，