## 0 賽題思路

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

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

1 退火算法原理

1.1 物理背景

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

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

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

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

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

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

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

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

2.1 算法流程

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

2.2算法實(shí)現(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                                      # 定溫內(nèi)部迭代次數(shù)self.rate = rate                                            # 退火降溫速度#############################################################self.x_seed = random.uniform(interval[0], interval[1])      # 解空間內(nèi)的種子self.tab = tab.strip()                                      # 求解最大值還是最小值的標(biāo)簽: 'min' - 最小值；'max' - 最大值#############################################################self.solve()                                                # 完成主體的求解過程self.display()                                              # 數(shù)據(jù)可視化展示def solve(self):temp = 'deal_' + self.tab                                   # 采用反射方法提取對(duì)應(yīng)的函數(shù)if hasattr(self, temp):deal = getattr(self, temp)else:exit('>>>tab標(biāo)簽傳參有誤："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]:   # 將隨機(jī)解束縛在給定狀態(tài)空間內(nè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):                                      # 計(jì)算最小值時(shí)，容忍解的狀態(tài)遷移概率probability = np.exp(-delta/T)return probabilitydef p_max(self, delta, T):probability = np.exp(delta/T)                               # 計(jì)算最大值時(shí)，容忍解的狀態(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')

實(shí)現(xiàn)結(jié)果