維基百科

模擬退火

模擬退火（英語：Simulated annealing，縮寫作SA）是一種逼近給定函數全局最佳的通用機率演算法，具體來說，它是一種元啟發算法，常用來在一定時間內，尋找在一個很大搜尋空間中的近似全局最優解。在有大量局部最佳解時，模擬退火算法可以找到全局最佳解。[1] 模擬退火常用於搜索空間離散的情形（如旅行推銷員問題、布林可滿足性問題、蛋白質結構預測、作業車間調度問題等）。對於在固定時間內找到近似全局最佳優先於找到精確局部最佳的問題，模擬退火算法可能優於梯度下降法或分支定界等精確方法。

模擬退火算法可用於求解組合問題。此處將模擬退火算法應用於求解旅行商問題，求出連接125個點的最小路線長度
使用模擬退火算法求解120個點的三維旅行商問題

模擬退火算法解決的問題包含多元目標函數與若干約束。實踐中，約束可作為目標函數的一部分進行懲罰。

Pincus (1970)、[2]Khachaturyan et al (1979,[3] 1981[4])、Kirkpatrick、Gelatt及Vecchi (1983)、Cerny (1985)等人先後提出過類似技術。[5]現在的「模擬退火」算法在1983年為S. Kirkpatrick, C. D. Gelatt和M. P. Vecchi解決旅行推銷員問題所發明[6]，V. Černý也在1985年獨立發明此演算法。

模擬退火算法中的慢冷卻是指在探索解空間的過程中，接受較差解的機率會緩慢下降。接受較差解可以更廣泛地搜索全局最佳解。總的來說，模擬退火算法的工作原理如下：溫度從初始值逐漸降低到0，每個時間步長內，算法隨機選擇一個與當前解接近的解，並根據與溫度相關的機率選擇更優的。搜索過程中，這機率會趨近於1。

可通過求解機率密度函數的動力方程式[7][8]或隨機採樣法進行模擬。[6][9]這種方法是N. Metropolis et al. (1953)發表的梅特羅波利斯-黑斯廷斯算法的改進版，是一種生成熱力學系統樣本狀態的蒙特卡羅方法。[10]

簡介

「模擬退火」來自冶金學術語退火，是將材料加熱後再經特定速率冷卻的技術，目的是增大晶粒的體積，並且減少晶格中的缺陷，以改變材料的物理性質。材料中的原子原來會停留在使內能有局部最小值的位置，加熱使能量變大，原子會離開原來位置，而隨機在其他位置中移動。退火冷卻時速度較慢，使得原子有較多可能可以找到內能比原先更低的位置。

模擬退火的原理也和金屬退火的原理近似：我們將熱力學的理論套用到統計學上，將搜尋空間內每一點想像成空氣內的分子；分子的能量，就是它本身的動能；而搜尋空間內的每一點，也像空氣分子一樣帶有「能量」，以表示該點對命題的合適程度。演算法先以搜尋空間內一個任意點作起始：每一步先選擇一個「鄰居」，然後再計算從現有位置到達「鄰居」的機率。

可以證明，模擬退火算法所得解依機率收斂到全局最佳解。

模擬退火法可用於精確算法失效的高難度計算最佳化問題，雖然通常只能獲得全局最佳的近似，但對很多實際問題已經足夠。

描述

物理系統的熱力學狀態s，以及要最小化的函數 （類似於系統當前狀態下的內能）。目標是將系統從任意初始狀態帶到內能儘可能小的狀態。

 
運用模擬退火算法搜索最大值，目標是達到最大值點。本例中，簡單的爬山算法是不夠的，因為有許多局部最值。緩慢降溫可以找到全局最值。

基本疊代

模擬退火啟發式在每一步都會考慮當前狀態s的某鄰態 ，並以機率決定是否移動到狀態 。機率最終引導系統進入低能狀態。這步驟一般會重複進行，直到系統達到滿足需求的狀態，或達到預設的計算量。

狀態的鄰域

解的最佳化包括計算問題狀態的鄰域，即由保守改變現狀態產生的新狀態。例如，在旅行推銷員問題中，狀態是待訪問城市的排列，狀態的鄰域是交換任意兩城市產生的排列集合。良定義的到鄰態的方法稱為移動，不同移動會產生不同的鄰態集。

爬山算法之類啟發法逐個尋找更好的鄰態來移動，並在無更好鄰態時停止，顯然這很容易陷入局部最佳。元啟發算法利用解的鄰域作為探索解空間的一種方式，雖然更喜歡較好的鄰態，但也接受較差的鄰態，以免陷入局部最佳。若運行時間夠長，則可以找到全局最佳。

接受機率

從現狀態s轉移到候選的新狀態 的機率由接受機率函數（acceptance probability function） 確定，其取決於兩狀態的能量 ，以及稱作溫度的全局時變參數T。轉移到能量更小的狀態的機率更大。機率函數P在 時也是正的，這可以防止算法陷入局部最佳。

T趨近於0時，若 ，機率 也要趨近於0或某正值。對足夠小的T，系統會越來越傾向於「下山」（向低能值移動）而避免「上山」。 時，程序將簡化為貪心算法，只進行下山轉移。

在模擬退火的原始描述中， 時機率 ，即無論溫度如何，程序找到下山的方法時總會下山。許多模擬退火算法的描述和實現仍將此條件作為方法定義的一部分，但這條件實際上並非必須。

P函數通常是這樣選擇的：當差值 增加時，接受較小上山轉移的機率就會比較大的更大。不過，在滿足上述要求的前提下，這要求並非絕對必要。

鑑於這些特性，溫度T在控制系統狀態s演化方面起到關鍵作用，因為它對系統能量變化非常敏感。確切地說，T的大小決定著s演化敏感的「粒度」。

退火歷程

 
快
 
慢
冷卻時間對模擬退火性能的影響。待解決問題是重新排列一幅圖像的像素，使某勢能函數最小，它會使較近的相似顏色相互吸引，較遠的則相斥，基本移動是像素交換位置。快速冷卻（左）與較慢冷卻（右）時，結果分別類似於無定形體與晶體。

算法名稱與靈感要求嵌入一個與溫度有關的有趣特徵。開始時，溫度T被設為一個較大值（或無窮大），按用戶指定的退火歷程，每次疊代降溫，但必須在一定時間預算內結束為 。這樣，預計系統最初會在搜索空間中包含良好解的廣闊區域內遊蕩，忽略能量函數的微小特徵；然後，向能量更低的區域漂移，最後根據梯度下降法啟發式向下移動。

對任意給定的有限問題，隨著退火歷程延長，模擬退火算法以全局最佳解終止的機率趨近於1。[11]但這理論結果並不很有用，因為確保顯著成功機率的耗時通常大於暴力搜索整個解空間的耗時。[12]

偽代碼

尋找能量   最低的狀態   

s := s0; e := E(s)                  // 设定目前状态为 s0，其能量 E(s0)
k := 0                              // 评估次数 k
while k < kmax and e > emin         // 若还有时间（评估次数 k 还不到 kmax）且结果还不够好（能量 e 不够低）则：
    sn := neighbour(s)              // 随机选取一邻近状态 sn
    en := E(sn)                     // sn 的能量为 E(sn)
    if random() < P(e, en, temp(k/kmax)) // 决定是否移至邻近状态 sn
        s := sn; e := en            // 移至邻近状态 sn
    k := k + 1                      // 评估完成，次数 k 加一
return s                            // 返回状态 s

下面以自然語言解說模擬退火算法的演算步驟。

初始化

由一個產生函數從當前解產生一個位於解空間的新解，並定義一個足夠大的數值作為初始溫度。

疊代過程

疊代過程是模擬退火算法的核心步驟，分為新解的產生和接受新解兩部分：

  1. 由一個產生函數從當前解產生一個位於解空間的新解；為便於後續的計算和接受，減少算法耗時，通常選擇由當前新解經過簡單地轉換即可產生新解的方法，如對構成新解的全部或部分元素進行置換、互換等，注意到產生新解的轉換方法決定了當前新解的鄰域結構，因而對冷卻進度表的選取有一定的影響。
  2. 計算與新解所對應的目標函數差。因為目標函數差僅由轉換部分產生，所以目標函數差的計算最好按增量計算。事實表明，對大多數應用而言，這是計算目標函數差的最快方法。
  3. 判斷新解是否被接受，判斷的依據是一個接受準則，最常用的接受準則是Metropolis準則：若Δt′<0則接受S′作為新的當前解S，否則以機率exp（-Δt′/T）接受S′作為新的當前解S。
  4. 當新解被確定接受時，用新解代替當前解，這只需將當前解中對應於產生新解時的轉換部分予以實現，同時修正目標函數值即可。此時，當前解實現了一次疊代。可在此基礎上開始下一輪試驗。而當新解被判定為捨棄時，則在原當前解的基礎上繼續下一輪試驗。

模擬退火算法與初始值無關，算法求得的解與初始解狀態S（是算法疊代的起點）無關；模擬退火算法具有漸近收斂性，已在理論上被證明是一種以機率1收斂於全局最佳解的全局最佳化算法；模擬退火算法具有並行性。

停止準則

疊代過程的一般停止準則：溫度T降低至某閾值時，或連續若干次疊代均未接受新解時，停止疊代，接受當前尋找的最佳解為最終解。

退火方案

在某個溫度狀態T下，當一定數量的疊代操作完成後，降低溫度T，在新的溫度狀態下執行下一個批次的疊代操作。

選擇參數

將模擬退火算法用於實際問題，需要指定狀態空間、能量（目標）函數E()、候選解生成器neighbour()、接受機率函數P()、退火歷程temperature()以及初始溫度init_temp。這些選擇會對算法效率產生重大影響，遺憾的是它們的選擇不適合所有問題，也沒有針對特定問題找到最佳選擇的通用方法。下面幾節將給出一些通用指導。

充分近的鄰態

模擬退火可建模為搜索圖上的隨機漫步，其頂點是所有狀態，邊為候選的移動。neighbour()的基本要求是提供從初態到任何可能是全局最佳態的足夠短路徑，這要求搜索圖的直徑足夠小。以上面的旅行推銷員問題為例，n = 20個城市的搜索空間有n! = 2,432,902,008,176,640,000（243億億）個狀態，而每個頂點有 條邊（來自n選2），圖直徑為 。

轉移機率

要研究模擬退火算法在特定問題上的行為，可考慮算法實施過程中各種設計選擇產生的轉移機率。對搜索圖中的每條邊 ，轉移機率定義為現態為s時轉移到 的機率，取決於temperature()指定的當前溫度、neighbour()生成的候選移動的順序及接受機率函數P()（注意轉移機率不是簡單的 ，因為候選轉移是以序列測試的）。

接受機率

neighbour()、P()、temperature()的指定是多餘的，因為實踐中通常會使用相同的接受函數P()，並根據具體問題調整其他函數。

在Kirkpatrick et al.的方法中，若 則接受機率函數 ，否則為 。表面上看，這個公式是通過類比物理系統的轉移來證明的；在 且梅特羅波利斯-黑斯廷斯的提議分布對稱時，對應梅特羅波利斯-黑斯廷斯算法。即使與梅特羅波利斯-黑斯廷斯算法中的提議分布類似的neighbour()函數不對稱或根本不是機率分布，也常用於模擬退火。因此，轉移機率同類似物理系統的轉移機率並不對應，恆定溫度T下的長期狀態分布也不必同物理系統在任何溫度下的熱力平衡態分布有任何相似處。儘管如此，大多數模擬退火的描述都會假定原接受函數，而這種函數可能在許多模擬退火的實現中都是硬編碼（hard-coded）的。

Moscato and Fontanari (1990)[13]與Dueck and Scheuer[14]分別提出，確定性更新（即不基於機率接受規則的）可加速最佳化，且不影響最終質量。Moscato and Fontanari通過觀察「閾值更新」退火的類「比熱」曲線得出結論：「在模擬退火算法中，Metropolis更新的隨機性在尋找近優最小值的過程中沒有發揮主要作用。」相反，他們認為「高溫下成本函數景觀的平滑化與冷卻過程中最小值的逐步確定是模擬退火算法成功的基本要素。」後來，由於Dueck and Scheuer的命名，這種方法也稱作「閾值接受法」。Franz, Hoffmann and Salamon (2001)表明，在一大類模擬成本/能量景觀隨機漫步的算法中，確定性更新策略確實是最佳策略。[15]

高效生成鄰態

選擇候選生成器neighbour()時，必須考慮到算法疊代幾次後，現態的能量將低於隨機狀態。因此一般來說，生成器應偏向更接近目標狀態 能量的鄰域。這種啟發式（也是梅特羅波利斯-黑斯廷斯算法的主要原理）往往會排除「非常好」與「非常壞」的轉移，不過前者通常更少見，所以啟發式往往相當有效。

例如，在上述旅行商問題中，低能旅行中交換兩連續城市對能量（總距離）影響不大，而交換兩任意城市更可能增加能量。因此，連續交換優於隨機交換，儘管後者能提供更優路徑（ 次交換，而非 次）。

啟發式一個更準確的詮釋是，應先嘗試 較大的鄰態 。對上述「標準」接受函數P來說，這意味著 在T或更小的數量級上因此在上述旅行商問題中，可用交換兩隨機城市的neighbour()，並設置兩城市距離超過T時不予選擇。

避免障礙

選擇候選生成器neighbour()時，還必須儘量減少能量遠低於所有鄰態的「深」局部最佳（或相連狀態集）數，這種「能量函數盆地」可能會以較高機率（與盆地中狀態數大致呈正比）與較長時間（與周圍狀態同盆地底部的能量差大致呈指數關係）困住模擬退火算法。

一般來說設計不出既能滿足這一目標，又能優先處理能量相近的候選粒子的候選發生器。另一方面，往往可對生成器進行相對簡單的修改，大大提高模擬退火的效率。例如，在旅行商問題中，不難找到兩條長度幾乎相等的路線 ，使得(1)  最佳、(2) 將 轉換為 的城市交換序列要經歷比二者長得多的距離、(3)   可通過翻轉（顛倒順序）一組連續的城市序列變為 。此例中， 就處於不同的「盆地」中，但若生成器執行隨機分段翻轉，則將位於同一個盆地。

冷卻歷程

模擬退火的物理類比假設冷卻速率足夠低，當前狀態的機率分布比任意時刻都接近熱力平衡。 然而，弛豫時間（relaxation time，溫度變化後恢復平衡耗時）很大程度上取決於能量函數的「地形」和當前溫度。模擬退火算法中，弛豫時間還以非常複雜的方式取決於候選生成器。注意所有參數通常作為黑盒函數提供給模擬退火算法。因此，理想的冷卻速度無法事先確定，只能根據經驗分析具體問題。自適應模擬退火算法將冷卻歷程同搜索進度聯繫起來，解決了這一問題。其他自適應法有熱力學模擬退火，[16]會根據熱力學定律，依兩狀態能量差自動調整每步的溫度。

重啟

有時，從現態出發不如回到明顯更優的解，這一過程稱作模擬退火的重啟（restarting）。為此，置s、e為sbest、ebest，然後重啟退火歷程。重啟的決定可基於多個標準，其中值得注意的是包括基於固定步數、基於當前能量與迄今最佳的能量、隨機重啟等等。

相關算法

  • 交互式梅特羅波利斯-黑斯廷斯算法（或順序蒙地卡羅[17]）將模擬退火轉移同配備交互式回收機制的最佳擬合個體的接受-剔除相結合。
  • 量子退火用「量子擾動」而非熱波動來穿越目標函數中較薄的障礙。
  • 隨機隧道法試圖以隧道穿越障礙，克服了隨溫度降低，模擬退火越來越難以擺脫局部最佳的問題。
  • 禁忌搜索通常會移動到能量較低的鄰態，發現陷入局部最佳時，就會採取上坡轉移，並會保留已知解列表（即「禁忌表」）以避免循環。
  • 雙相演化是一組算法與過程（SA屬於其中），利用搜索空間的相變，在局部搜索與全局搜索之間進行調節。
  • 反應式搜索最佳化將機器學習與最佳化相結合，增加內部反饋迴路，根據問題、實例與當前解周圍的情形，對算法的自由參數進行自我調整。
  • 遺傳算法保留一組解（「池」，pool），而非一個解。新的候選解可由突變（如SA）和「重組」從池中產生。機率標準和SA中所用的類似，用於選擇突變或重組的解，並剔除多餘解。
  • 文化基因算法使用一組既合作又競爭的搜索器（agent）搜索解。有時搜索器的策略設計模擬退火，以便在重組之前獲得高質量解。[18]退火也被認為是增加搜索多樣性的一種機制。[19]
  • 漸進最佳化在最佳化過程中對目標函數進行「平滑化」處理。
  • 蟻群算法（ACO）用許多搜索器遍歷解空間，並找到局部較優的區域。
  • 交叉熵法（CE）通過參數化機率分布生成候選解。參數通過交叉熵最小化進行更新，以便在下次疊代中生成更好的樣本。
  • 和聲搜索模擬音樂家即興演奏的過程，每個音樂家演奏一個音符，共同尋找最佳的和聲。
  • 隨機最佳化是一系列方法的總稱，包括模擬退火等很多方法。
  • 粒子群最佳化是一種以群體智能為模型的算法，能在搜索空間中找到最佳化問題的解，或在有目標的情形下模擬、預測群體行為。
  • 莖根算法（runner-root algorithm，RRA）是一種元啟發最佳化算法，用於解決單模態和多模態問題，靈感來自植物的莖與葉的關係。
  • 智能水滴算法（Intelligent water drops algorithm，IWD）模仿自然水滴的行為解決最佳化問題。
  • 並行退火在不同溫度（或哈密頓量）下模擬模型副本，以克服潛在障礙。
  • 多目標模擬退火用於多目標最佳化。[20]

另見

參考資料

  1. ^ What is Simulated Annealing?. www.cs.cmu.edu. [2023-05-13]. （原始內容存檔於2024-08-29）. 
  2. ^ Pincus, Martin. A Monte-Carlo Method for the Approximate Solution of Certain Types of Constrained Optimization Problems. Journal of the Operations Research Society of America. Nov–Dec 1970, 18 (6): 967–1235. doi:10.1287/opre.18.6.1225. 
  3. ^ Khachaturyan, A.: Semenovskaya, S.: Vainshtein B., Armen. Statistical-Thermodynamic Approach to Determination of Structure Amplitude Phases. Soviet Physics, Crystallography. 1979, 24 (5): 519–524. 
  4. ^ Khachaturyan, A.; Semenovskaya, S.; Vainshtein, B. The Thermodynamic Approach to the Structure Analysis of Crystals. Acta Crystallographica. 1981, A37 (5): 742–754. Bibcode:1981AcCrA..37..742K. doi:10.1107/S0567739481001630. 
  5. ^ Laarhoven, P. J. M. van (Peter J. M.). Simulated annealing : theory and applications. Aarts, E. H. L. (Emile H. L.). Dordrecht: D. Reidel. 1987. ISBN 90-277-2513-6. OCLC 15548651. 
  6. ^ 6.0 6.1 Kirkpatrick, S.; Gelatt Jr, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science. 1983, 220 (4598): 671–680. Bibcode:1983Sci...220..671K. CiteSeerX 10.1.1.123.7607  . JSTOR 1690046. PMID 17813860. S2CID 205939. doi:10.1126/science.220.4598.671. 
  7. ^ Khachaturyan, A.; Semenovskaya, S.; Vainshtein, B. Statistical-Thermodynamic Approach to Determination of Structure Amplitude Phases. Sov.Phys. Crystallography. 1979, 24 (5): 519–524. 
  8. ^ Khachaturyan, A.; Semenovskaya, S.; Vainshtein, B. The Thermodynamic Approach to the Structure Analysis of Crystals. Acta Crystallographica. 1981, 37 (A37): 742–754. Bibcode:1981AcCrA..37..742K. doi:10.1107/S0567739481001630. 
  9. ^ Černý, V. Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm. Journal of Optimization Theory and Applications. 1985, 45: 41–51. S2CID 122729427. doi:10.1007/BF00940812. 
  10. ^ Metropolis, Nicholas; Rosenbluth, Arianna W.; Rosenbluth, Marshall N.; Teller, Augusta H.; Teller, Edward. Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics. 1953, 21 (6): 1087. Bibcode:1953JChPh..21.1087M. OSTI 4390578. S2CID 1046577. doi:10.1063/1.1699114. 
  11. ^ Granville, V.; Krivanek, M.; Rasson, J.-P. Simulated annealing: A proof of convergence. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1994, 16 (6): 652–656. doi:10.1109/34.295910. 
  12. ^ Nolte, Andreas; Schrader, Rainer, A Note on the Finite Time Behaviour of Simulated Annealing, Operations Research Proceedings 1996 1996 (Berlin, Heidelberg: Springer Berlin Heidelberg), 1997, 1996: 175–180 [2023-02-06], ISBN 978-3-540-62630-5, doi:10.1007/978-3-642-60744-8_32 
  13. ^ Moscato, P.; Fontanari, J.F., Stochastic versus deterministic update in simulated annealing, Physics Letters A, 1990, 146 (4): 204–208, Bibcode:1990PhLA..146..204M, doi:10.1016/0375-9601(90)90166-L 
  14. ^ Dueck, G.; Scheuer, T., Threshold accepting: A general purpose optimization algorithm appearing superior to simulated annealing, Journal of Computational Physics, 1990, 90 (1): 161–175, Bibcode:1990JCoPh..90..161D, ISSN 0021-9991, doi:10.1016/0021-9991(90)90201-B 
  15. ^ Franz, A.; Hoffmann, K.H.; Salamon, P, Best optimal strategy for finding ground states, Physical Review Letters, 2001, 86 (3): 5219–5222, PMID 11384462, doi:10.1103/PhysRevLett.86.5219 
  16. ^ De Vicente, Juan; Lanchares, Juan; Hermida, Román. Placement by thermodynamic simulated annealing. Physics Letters A. 2003, 317 (5–6): 415–423. Bibcode:2003PhLA..317..415D. doi:10.1016/j.physleta.2003.08.070. 
  17. ^ Del Moral, Pierre; Doucet, Arnaud; Jasra, Ajay. Sequential Monte Carlo samplers. Journal of the Royal Statistical Society, Series B. 2006, 68 (3): 411–436. S2CID 12074789. arXiv:cond-mat/0212648  . doi:10.1111/j.1467-9868.2006.00553.x. 
  18. ^ Moscato, Pablo. An introduction to population approaches for optimization and hierarchical objective functions: A discussion on the role of tabu search. Annals of Operations Research. June 1993, 41 (2): 85–121. S2CID 35382644. doi:10.1007/BF02022564. 
  19. ^ Moscato, P. On Evolution, Search, Optimization, Genetic Algorithms and Martial Arts: Towards Memetic Algorithms. Caltech Concurrent Computation Program. 1989, (report 826). 
  20. ^ Deb, Bandyopadhyay. A Simulated Annealing-Based Multiobjective Optimization Algorithm: AMOSA. IEEE Transactions on Evolutionary Computation. June 2008, 12 (3): 269–283. S2CID 12107321. doi:10.1109/TEVC.2007.900837. 

延伸閲讀

參閲

外部連結

取自「https://zh.wikipedia.org/w/index.php?title=模拟退火&oldid=84243853」