模拟退火(SA)算法介紹和應用細節-附SA結合登山算法求解VRPTW問題C++代碼

文章目錄

• 1、模拟退火算法介紹
• • 1.1 精确SA算法
  • • 1.1.1 Metropolis Acceptance Criterion
    • 1.1.2 Boltzmann distribution
  • 1.2 啟發式SA算法
  • • 1.2.1 溫度處理
• 2、結合SA的登山算法
• • 2.1 Hill Climbing Pseudo Code
  • 2.2 應用于VRPTW問題
  • 2.3 VRPTW求解c++代碼已上傳至GitHub，下面給對外連結接：

1、模拟退火算法介紹

\qquad 模拟退火算法(Simulated Annealing)是一種用于求解離散優化和部分連續優化問題的有效的元啟發式算法。SA算法的核心特點在于它通過接收更差的解來輔助跳出局部最優解，以求得到全局最優解。

\qquad SA算法起源于晶體降溫的實體實驗，将傳熱動力學的行為和離散優化問題最優解的尋優進行互相關聯，最終将這種關聯成功應用到了優化問題上。在SA求解的過程中，較優的解總是被接受，但是一部分次優的解也會被接受，接收次優解的機率取決于溫度參數，而溫度參數在疊代的過程中通常是一個非增(non-increasing)的參數。

\qquad 模拟退火算法可以使用Makov鍊來證明最優收斂性，下面介紹精确SA算法的理論内容。

1.1 精确SA算法

1.1.1 Metropolis Acceptance Criterion

\qquad 首先定義 Ω \Omega Ω表示求解空間；定義 f f f表示目标函數；定義 ω ∗ \omega^* ω∗表示最優解，則 f ( ω ) ≥ f ( ω ∗ ) ， ω ∈ Ω f(\omega)\geq f(\omega^*)，\omega \in \Omega f(ω)≥f(ω∗)，ω∈Ω；定義鄰域空間 N ( ω ) ， ω ∈ Ω N(\omega)，\omega \in \Omega N(ω)，ω∈Ω表示LS算法通過一步疊代可以找到的鄰域解集合。SA算法基于Metropolis Acceptance Criterion來判斷是否接受新的較差的鄰域解 ω \omega ω來替代目前解。Metropolis Acceptance Criterion的機率計算如下所示： P ( A c c e p   ω ′   a s   n e x t   S o l u t i o n ) = { e x p [ − ( f ( ω ′ ) − f ( ω ) ) / t k ] , if  f ( ω ′ ) − f ( ω ) > 0 1 , if  f ( ω ′ ) − f ( ω ) ≤ 0 P(Accep \ \omega'\ as\ next\ Solution)=\begin{cases} exp[-(f(\omega')-f(\omega))/t_k], & \text{if $f(\omega')-f(\omega)>0$} \\ 1, & \text{if $f(\omega')-f(\omega)\leq0$} \\ \end{cases} P(Accep ω′ as next Solution)={exp[−(f(ω′)−f(ω))/tk​],1,​if f(ω′)−f(ω)>0if f(ω′)−f(ω)≤0​ \qquad 其中 t k t_k tk​表示和外循環次數 k k k相關的溫度參數，滿足： t k > 0       l i m k → + ∞ t k = 0 t_k>0 \ \ \ \ \ \underset{k→+\infty}{lim}t_k=0 tk​>0     k→+∞lim​tk​=0 \qquad 上式是基本SA算法應用的接收準則。

1.1.2 Boltzmann distribution

\qquad 進一步的搜尋方向選擇需要用到Boltzmann distribution等式，其描述了在溫度 T T T情況下，鄰域解 ω \omega ω出現的機率： P ( S y s t e m   i s   i n   s t a t e   ω   a t   t e m p e r a t u r e   T ) = e x p ( − f ( ω ) / t k ) ∑ ω ′ ′ ∈ Ω e x p ( − f ( ω ′ ′ ) / t k ) P(System\ is\ in \ state\ \omega\ at\ temperature\ T)=\frac{exp(-f(\omega)/t_k)}{\sum_{\omega''\in\Omega}exp(-f(\omega'')/t_k)} P(System is in state ω at temperature T)=∑ω′′∈Ω​exp(−f(ω′′)/tk​)exp(−f(ω)/tk​)​ \qquad 還需要用到選擇鄰域解 ω ′ \omega' ω′的機率 g k ( ω ,   ω ′ ) g_k(\omega,\ \omega') gk​(ω, ω′)，其滿足： ∑ ω ′ ∈ N ( ω ) g k ( ω ,   ω ′ ) = 1 ,   f o r a l l   ω ∈ Ω ,   k = 1 , 2 , . . . \underset{\omega'\in N(\omega)}{\sum} g_k(\omega,\ \omega')=1,\ forall\ \omega \in \Omega,\ k=1,2,... ω′∈N(ω)∑​gk​(ω, ω′)=1, forall ω∈Ω, k=1,2,... \qquad 則在某一次疊代時，選擇某個鄰域解 ω \omega ω進行更新目前解的機率為： P k ( ω ,   ω ′ ) = { g k ( ω ,   ω ′ ) e x p ( − Δ ω , ω ′ / t k ) ,    i f   ω ∈ N ( ω ) ,   ω ≠ ω 0 ,    i f   ω ∉ N ( ω ) ,   ω ≠ ω 1 − ∑ ω ′ ′ ∈ N ( ω ) ω ′ ′ ≠ ω p k ( ω ,   ω ′ ′ ) ,       i f   ω = ω P_k(\omega,\ \omega')=\begin{cases} g_k(\omega,\ \omega')exp(-\Delta_{\omega,\omega'}/t_k),\ \ if\ \omega \in N(\omega), \ \omega \neq \omega \\ 0,\qquad \qquad\qquad\qquad\qquad\ \ if\ \omega \notin N(\omega), \ \omega \neq \omega \\ 1-\underset{\omega''\neq\omega}{\underset{\omega''\in N(\omega)}{\sum}}p_k(\omega,\ \omega''),\quad \ \ \ \ \ if\ \omega=\omega \end{cases} Pk​(ω, ω′)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​gk​(ω, ω′)exp(−Δω,ω′​/tk​),  if ω∈N(ω), ω​=ω0,  if ω∈/​N(ω), ω​=ω1−ω′′​=ωω′′∈N(ω)∑​​pk​(ω, ω′′),     if ω=ω​ \qquad 從上式獲得每個鄰域解的選擇機率之後，選擇機率最大的鄰域解作為此次搜尋的方向，繼續執行搜尋，這種精确SA算法可以通過Makov Chain證明其可以收斂到最優解的最優性。

1.2 啟發式SA算法

\qquad 上述精确SA算法最大的缺陷在于每一步計算所有鄰域解會耗費大量的時間，導緻其求解時間很長。是以在實際應用之中，通常将SA算法和其他啟發式算法進行結合，借助SA的思想來跳出局部最優解，通過快速疊代尋優，最終得到一個較好的近似最優解。

\qquad 是以啟發式SA算法隻需要用到上述Metropolis Acceptance Criterion的準則便已經足夠，在每一次建構出一個新的鄰域解之後，判斷其是否優于目前解，若是則直接用其替換目前解；若不是，則計算一個接收次優解的機率 P P P，之後再随機生成一個機率 P r P^r Pr，若 P r ≥ P P^r\geq P Pr≥P，則接受目前次優解，否則放棄目前次優解。

1.2.1 溫度處理

\qquad 在應用SA算法時，一個重要的考察項是溫度相關的參數設定，包括初始溫度，降溫速率的選擇和停止降溫的标準。本人在應用SA時設定初始溫度 t 0 t_0 t0​的方式如下： t 0 = p v ∗ f ( i n i t ) / l n ( 0.5 ) t_0=p^v*f(init)/ln(0.5) t0​=pv∗f(init)/ln(0.5) \qquad 其中 p v p^v pv表示設定的降溫速率參數，一般取值為0.005，這樣降溫速率即為0.99975；本人在實際應用中設定停止降溫的準則為溫度小于初始溫度的 1 / 3 1/3 1/3，之後會使用目前最優解 f ( p r e ) f(pre) f(pre)重新計算一個初始溫度： t 0 = p v ∗ f ( p r e ) / l n ( 0.5 ) t_0=p^v*f(pre)/ln(0.5) t0​=pv∗f(pre)/ln(0.5) \qquad 上述退火速率是一個靜态的調整方式，但更加有效的方式可以使用自适應的方式，根據目标函數值等資訊來動态調整退火速率。

2、結合SA的登山算法

\qquad 登山算法為LS提供了一個通用的架構來求解棘手的離散優化問題。所有的登山算法都有相似的結構，隻在下述兩個部分稍有改動。第一部分為登山随機變量，它用來判斷是接收一個解還是拒絕一個解；第二部分為鄰域結構，這個跟問題的特征息息相關。下面使用SA作為登山随機變量，給出嵌入SA的登山算法僞代碼。

2.1 Hill Climbing Pseudo Code

1：選擇一個初始解 ω ∈ Ω \omega \in \Omega ω∈Ω，初始化初始溫度 t 0 = − 0.005 ∗ f ( ω ) / l n ( 0.5 ) t_0=-0.005*f(\omega)/ln(0.5) t0​=−0.005∗f(ω)/ln(0.5)

2：設定外循環最大疊代次數 K K K 和内循環最大疊代次數數組 M k ,   k = 1 , 2 , . . . , K M_k,\ k=1,2,...,K Mk​, k=1,2,...,K

3：設定外/内循環疊代計數變量 k = m = 1 k=m=1 k=m=1

4：Repeat while k ≤ K k \leq K k≤K

5： \qquad Repeat while m ≤ M k m \leq M_k m≤Mk​

6： \qquad 生成一個鄰域解 ω ′ ∈ N ( ω ) \omega' \in N(\omega) ω′∈N(ω)

7： \qquad 計算目前解和鄰域解之間的內插補點 Δ ω , ω ′ = f ( ω ′ ) − f ( ω ) \Delta_{\omega, \omega'}=f(\omega')-f(\omega) Δω,ω′​=f(ω′)−f(ω)

8： \qquad If Δ ω , ω ′ ≤ 0 \Delta_{\omega, \omega'} \leq0 Δω,ω′​≤0 then ω ← ω ′ \omega ←\omega' ω←ω′

9： \qquad If Δ ω , ω ′ > 0 \Delta_{\omega, \omega'} >0 Δω,ω′​>0 then ω ← ω ′ \omega ←\omega' ω←ω′ 依機率 e x p ( − Δ ω , ω ′ / t k ) exp(-\Delta_{\omega, \omega'}/t_k) exp(−Δω,ω′​/tk​)

10: m ← m + 1 \qquad m←m+1 m←m+1

11: \qquad Until m = M k m=M_k m=Mk​

12: 更新溫度參數 t k ← t k ∗ 0.99975 t_k←t_k*0.99975 tk​←tk​∗0.99975

13: If t k < 1 3 ∗ t 0 t_k<\frac{1}{3}*t_0 tk​<31​∗t0​ Then 重新設定 t 0 = − 0.005 ∗ f ( ω ) / l n ( 0.5 ) t_0=-0.005*f(\omega)/ln(0.5) t0​=−0.005∗f(ω)/ln(0.5)

14: k ← k + 1 k←k+1 k←k+1

15: Until k = K k=K k=K 或者其他終止準則滿足結束算法

2.2 應用于VRPTW問題

\qquad 在SA+HC算法應用到VRPTW問題時，鄰域結構選擇移除和插入一個客戶點，LS算法細節可見：

https://blog.csdn.net/weixin_43160744/article/details/119107248

\qquad 之後套用上述算法架構即可。試驗證明，加入了SA的HC比單純的HC效果要好許多，在收斂速度和求解品質上都有明顯的改善。

2.3 VRPTW求解c++代碼已上傳至GitHub，下面給對外連結接：

https://github.com/Dragon-wang-fei-long/Huristic-code/tree/Dragon-wang-fei-long-HC+SA