## 0 賽題思路

（賽題出來以后第一時間在CSDN分享）

https://blog.csdn.net/dc_sinor?type=blog

1 退火算法原理

1.1 物理背景

在熱力學上，退火（annealing）現(xiàn)象指物體逐漸降溫的物理現(xiàn)象，溫度愈低，物體的能量狀態(tài)會低；夠低后，液體開始冷凝與結晶，在結晶狀態(tài)時，系統(tǒng)的能量狀態(tài)最低。大自然在緩慢降溫（亦即，退火）時，可“找到”最低能量狀態(tài)：結晶。但是，如果過程過急過快，快速降溫（亦稱「淬煉」，quenching）時，會導致不是最低能態(tài)的非晶形。

如下圖所示，首先（左圖）物體處于非晶體狀態(tài)。我們將固體加溫至充分高（中圖），再讓其徐徐冷卻，也就退火（右圖）。加溫時，固體內部粒子隨溫升變?yōu)闊o序狀，內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態(tài)，最后在常溫時達到基態(tài)，內能減為最小（此時物體以晶體形態(tài)呈現(xiàn)）。

1.2 背后的數(shù)學模型

如果你對退火的物理意義還是暈暈的，沒關系我們還有更為簡單的理解方式。想象一下如果我們現(xiàn)在有下面這樣一個函數(shù)，現(xiàn)在想求函數(shù)的（全局）最優(yōu)解。如果采用Greedy策略，那么從A點開始試探，如果函數(shù)值繼續(xù)減少，那么試探過程就會繼續(xù)。而當?shù)竭_點B時，顯然我們的探求過程就結束了（因為無論朝哪個方向努力，結果只會越來越大）。最終我們只能找打一個局部最后解B。

根據(jù)Metropolis準則，粒子在溫度T時趨于平衡的概率為exp(-ΔE/(kT))，其中E為溫度T時的內能，ΔE為其改變數(shù),k為Boltzmann常數(shù)。Metropolis準則常表示為

Metropolis準則表明，在溫度為T時，出現(xiàn)能量差為dE的降溫的概率為P(dE)，表示為：P(dE) = exp( dE/(kT) )。其中k是一個常數(shù)，exp表示自然指數(shù)，且dE<0。所以P和T正相關。這條公式就表示：溫度越高，出現(xiàn)一次能量差為dE的降溫的概率就越大；溫度越低，則出現(xiàn)降溫的概率就越小。又由于dE總是小于0（因為退火的過程是溫度逐漸下降的過程），因此dE/kT < 0 ，所以P(dE)的函數(shù)取值范圍是(0,1) 。隨著溫度T的降低，P(dE)會逐漸降低。

我們將一次向較差解的移動看做一次溫度跳變過程，我們以概率P(dE)來接受這樣的移動。也就是說，在用固體退火模擬組合優(yōu)化問題，將內能E模擬為目標函數(shù)值 f，溫度T演化成控制參數(shù) t，即得到解組合優(yōu)化問題的模擬退火演算法：由初始解 i 和控制參數(shù)初值 t 開始，對當前解重復“產(chǎn)生新解→計算目標函數(shù)差→接受或丟棄”的迭代，并逐步衰減 t 值，算法終止時的當前解即為所得近似最優(yōu)解，這是基于蒙特卡羅迭代求解法的一種啟發(fā)式隨機搜索過程。退火過程由冷卻進度表(Cooling Schedule)控制，包括控制參數(shù)的初值 t 及其衰減因子Δt 、每個 t 值時的迭代次數(shù)L和停止條件S。

2 退火算法實現(xiàn)

2.1 算法流程

(1) 初始化：初始溫度T(充分大)，初始解狀態(tài)S(是算法迭代的起點)， 每個T值的迭代次數(shù)L
(2) 對k=1，……，L做第(3)至第6步：
(3) 產(chǎn)生新解S′
(4) 計算增量Δt′=C(S′)-C(S)，其中C(S)為評價函數(shù)
(5) 若Δt′<0則接受S′作為新的當前解，否則以概率exp(-Δt′/T)接受S′作為新的當前解.
(6) 如果滿足終止條件則輸出當前解作為最優(yōu)解，結束程序。
終止條件通常取為連續(xù)若干個新解都沒有被接受時終止算法。
(7) T逐漸減少，且T->0，然后轉第2

2.2算法實現(xiàn)

import numpy as np
import matplotlib.pyplot as plt
import randomclass SA(object):def __init__(self, interval, tab='min', T_max=10000, T_min=1, iterMax=1000, rate=0.95):self.interval = interval                                    # 給定狀態(tài)空間 - 即待求解空間self.T_max = T_max                                          # 初始退火溫度 - 溫度上限self.T_min = T_min                                          # 截止退火溫度 - 溫度下限self.iterMax = iterMax                                      # 定溫內部迭代次數(shù)self.rate = rate                                            # 退火降溫速度#############################################################self.x_seed = random.uniform(interval[0], interval[1])      # 解空間內的種子self.tab = tab.strip()                                      # 求解最大值還是最小值的標簽: 'min' - 最小值；'max' - 最大值#############################################################self.solve()                                                # 完成主體的求解過程self.display()                                              # 數(shù)據(jù)可視化展示def solve(self):temp = 'deal_' + self.tab                                   # 采用反射方法提取對應的函數(shù)if hasattr(self, temp):deal = getattr(self, temp)else:exit('>>>tab標簽傳參有誤："min"|"max"<<<')x1 = self.x_seedT = self.T_maxwhile T >= self.T_min:for i in range(self.iterMax):f1 = self.func(x1)delta_x = random.random() * 2 - 1if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]:   # 將隨機解束縛在給定狀態(tài)空間內x2 = x1 + delta_xelse:x2 = x1 - delta_xf2 = self.func(x2)delta_f = f2 - f1x1 = deal(x1, x2, delta_f, T)T *= self.rateself.x_solu = x1                                            # 提取最終退火解def func(self, x):                                              # 狀態(tài)產(chǎn)生函數(shù) - 即待求解函數(shù)value = np.sin(x**2) * (x**2 - 5*x)return valuedef p_min(self, delta, T):                                      # 計算最小值時，容忍解的狀態(tài)遷移概率probability = np.exp(-delta/T)return probabilitydef p_max(self, delta, T):probability = np.exp(delta/T)                               # 計算最大值時，容忍解的狀態(tài)遷移概率return probabilitydef deal_min(self, x1, x2, delta, T):if delta < 0:                                               # 更優(yōu)解return x2else:                                                       # 容忍解P = self.p_min(delta, T)if P > random.random(): return x2else: return x1def deal_max(self, x1, x2, delta, T):if delta > 0:                                               # 更優(yōu)解return x2else:                                                       # 容忍解P = self.p_max(delta, T)if P > random.random(): return x2else: return x1def display(self):print('seed: {}\nsolution: {}'.format(self.x_seed, self.x_solu))plt.figure(figsize=(6, 4))x = np.linspace(self.interval[0], self.interval[1], 300)y = self.func(x)plt.plot(x, y, 'g-', label='function')plt.plot(self.x_seed, self.func(self.x_seed), 'bo', label='seed')plt.plot(self.x_solu, self.func(self.x_solu), 'r*', label='solution')plt.title('solution = {}'.format(self.x_solu))plt.xlabel('x')plt.ylabel('y')plt.legend()plt.savefig('SA.png', dpi=500)plt.show()plt.close()if __name__ == '__main__':SA([-5, 5], 'max')

實現(xiàn)結果