0 賽題思路

（賽題出來以后第一時間在群內(nèi)分享）

2023高教社杯全國大學生數(shù)學建模競賽

資料思路分享Q群：575833716

1 退火算法原理

1.1 物理背景

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

如下圖所示，首先（左圖）物體處于非晶體狀態(tài)。我們將固體加溫至充分高（中圖），再讓其徐徐冷卻，也就退火（右圖）。加溫時，固體內(nèi)部粒子隨溫升變?yōu)闊o序狀，內(nèi)能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態(tài)，最后在常溫時達到基態(tài)，內(nèi)能減為最?。ù藭r物體以晶體形態(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時的內(nèi)能，Δ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)化問題，將內(nèi)能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，然后轉(zhuǎn)第2

2.2算法實現(xiàn)

import numpy as np
import matplotlib.pyplot as plt
import random

class 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() ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 求解最大值還是最小值的標簽: '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_seed
? ? ? ?T = self.T_max
? ? ? ?while T >= self.T_min:
? ? ? ? ? ?for i in range(self.iterMax):
? ? ? ? ? ? ? ?f1 = self.func(x1)
? ? ? ? ? ? ? ?delta_x = random.random() * 2 - 1
? ? ? ? ? ? ? ?if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]: ? # 將隨機解束縛在給定狀態(tài)空間內(nèi)
? ? ? ? ? ? ? ? ? ?x2 = x1 + delta_x
? ? ? ? ? ? ? ?else:
? ? ? ? ? ? ? ? ? ?x2 = x1 - delta_x
? ? ? ? ? ? ? ?f2 = self.func(x2)
? ? ? ? ? ? ? ?delta_f = f2 - f1
? ? ? ? ? ? ? ?x1 = deal(x1, x2, delta_f, T)
? ? ? ? ? ?T *= self.rate
? ? ? ?self.x_solu = x1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 提取最終退火解

? ?def func(self, x): ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 狀態(tài)產(chǎn)生函數(shù) - 即待求解函數(shù)
? ? ? ?value = np.sin(x**2) * (x**2 - 5*x)
? ? ? ?return value

? ?def p_min(self, delta, T): ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 計算最小值時，容忍解的狀態(tài)遷移概率
? ? ? ?probability = np.exp(-delta/T)
? ? ? ?return probability

? ?def p_max(self, delta, T):
? ? ? ?probability = np.exp(delta/T) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 計算最大值時，容忍解的狀態(tài)遷移概率
? ? ? ?return probability

? ?def deal_min(self, x1, x2, delta, T):
? ? ? ?if delta < 0: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 更優(yōu)解
? ? ? ? ? ?return x2
? ? ? ?else: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 容忍解
? ? ? ? ? ?P = self.p_min(delta, T)
? ? ? ? ? ?if P > random.random(): return x2
? ? ? ? ? ?else: return x1

? ?def deal_max(self, x1, x2, delta, T):
? ? ? ?if delta > 0: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 更優(yōu)解
? ? ? ? ? ?return x2
? ? ? ?else: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 容忍解
? ? ? ? ? ?P = self.p_max(delta, T)
? ? ? ? ? ?if P > random.random(): return x2
? ? ? ? ? ?else: return x1

? ?def 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)結果

建模資料

資料分享: 最強建模資料

2023高教社杯全國大學生數(shù)學建模競賽

資料思路分享Q群：575833716