模擬退火理論（SA——Simulated Annealing）
和局部束搜索相反，模擬退火將最優(yōu)化策略改變，引入隨機噪聲，不一定每次都是最優(yōu)，但是內部機制保證了最終的走向是最優(yōu)，總的過程可以理解為初期廣泛探索（Exploration），逐步過渡到深挖（Exploitation）。其中機理比較復雜，我們逐步去理解。
基本思路
首先聲明，我們這里還是要將函數(shù)優(yōu)化到最大，如果要優(yōu)化最小，符號都要變。
所謂退火就是從高溫逐步降溫的過程。在這個過程中，初期高溫，狀態(tài)容易改變，不穩(wěn)定，后期低溫，狀態(tài)穩(wěn)定。我們模擬退火中進行隨機變化的概率就和溫度有關。這里引入溫度變量T，T怎么變呢？一般來說是指數(shù)變化 T = 0.9 9 t T 0 T=0.99^tT_0 T=0.99tT0?，這個0.99是退火率，可以進行修改，t是時間。之所以采用指數(shù)變化而不是什么線性遞減，也有門道，后面解釋。
看代碼，隨著循環(huán)，t增大，T以指數(shù)速率先快后慢變小。之后進行鄰域中的隨機選擇，而不是在鄰域中找最優(yōu)。這時結果就只能有兩種，我們通過新結果和原結果的差距 △ E \triangle E △E判斷：
- △ E > 0 \triangle E>0 △E>0，更優(yōu)。那我們就直接笑納。這個效果上相當于登山算法。
- △ E ≤ 0 \triangle E\leq 0 △E≤0, 相等或者更差，總之要么確實是不好的結果，要么再走走可能會碰到更好的結果，先苦后甜。那到底怎么抉擇呢，這可能是個機遇，但是也有風險，所以我們要計算概率。
概率公式初步理解
接受壞結果的概率 P = e △ E T P=e^{\frac{\triangle E}{T}} P=eT△E?，解析一下這個式子。
- 首先指數(shù)部分一定是負的，所以保證P是介于0-1之間的。
- 初期T大，不管E怎么變，比值總是接近于0，概率較大，所以很容易接受所有改動
- 中期T適中，由 ∣ △ E ∣ |\triangle E| ∣△E∣的大小決定是否接受，delta絕對值越大，接受概率越小，因為巨大的退步代表巨大的風險，合理。
- 后期T小，一般情況比值都更接近于負無窮，只有在 ∣ △ E ∣ |\triangle E| ∣△E∣特別小時能有接近于0的比值，才能接受，也就是說后期要么就是更優(yōu)，直接接受，要么就是產生一個更差結果，我大概率不接受而保持原狀，所以后期就相當于一個局部搜索。這里可見指數(shù)優(yōu)點，在低溫區(qū)變化夠慢，可以保證足夠的時間達到局部最優(yōu)，保證收斂。其實還可以采取別的，比如現(xiàn)實世界的溫度變化函數(shù)，但是規(guī)律一定是先快后慢。
如何理解步長逐漸縮小
我們進一步從步長的角度理解，這個步長是指 ∣ △ E ∣ |\triangle E| ∣△E∣，而這個常常和解空間中點的距離有較大的關系，所以我們就叫他步長：
- 剛開始T大，所以 ∣ △ E ∣ |\triangle E| ∣△E∣即使很大，也可以走出去，所以初期步長可大可小，可以接受各種結果。
- 但是到了后期，T變小， ∣ △ E ∣ |\triangle E| ∣△E∣一旦變大，給的概率就會變小，步子越大，越不容易接受，所以T越小，小步長的概率就越大，從這個角度理解，隨著T減小，平均步長是會越來越小的。
宏觀理解
再進行宏觀理解，理解下這個算法的合理性與有效性：
- 從宏觀角度理解，如果你把每次計算出來的概率輸出，你會發(fā)現(xiàn)接受較差結果的概率是越來越小。
- 合理性。不需要擔心算法質量越來越差，因為概率保證了質量。因為是隨機搜索的，所以假定產生好結果的概率是0.5，那么產生壞結果的前提下再接受壞結果的概率肯定是小于0.5的，這就保證了在任何時候，總體上接受好結果的概率要大于接受壞結果的概率，保證方向是正確的。
- 收斂性。初期對壞結果的接受度很高，所以屬于橫沖直撞的廣泛探索，不會收斂。后期溫度降低，且保持足夠長時間的低溫。低溫狀態(tài)下，基本只會朝著優(yōu)化方向走，所以在低溫狀態(tài)下相當于一個登山算法，而登山算法只要迭代次數(shù)夠多，必然能達到局部最優(yōu)。所以一定會收斂。
- 收斂較慢的特性保證了解的優(yōu)秀性。
局部搜索實踐：北理工最優(yōu)化大作業(yè)——流水車間調度
流水車間問題入門
基于改進模擬退火算法的大規(guī)模置換流水車間調度_黎陽
n個器件，在m個機器上加工。
每個器件在一個機器上加工的時間都不同，會給出一個矩陣，告訴你哪個器件在哪個機器上的耗時。
求最好的順序。
為什么順序會影響時間呢？那是因為一個器件的加工，要等自己在前一個機器的加工完成，又要等上一個器件在當前機器的加工完成，這兩個等待時間和順序有關。
對于一種順序，時間計算是很容易的，我后面也有現(xiàn)成的代碼。
初探：鄰域如何確定
問題來了，狀態(tài)轉移比較難。人家八皇后可以通過沖突數(shù)來確定狀態(tài)變化，我們通過什么確定呢？或者準確說，我們在什么鄰域里搜索呢？
既然是順序，那么我每次變化就只能變化順序，必須保證不沖突，就沒辦法對一個器件進行單獨變化。但是可以明確的是，每次變化一定是交換兩個工件的順序，所以如何選定兩個工件就是問題。是否存在一種啟發(fā)性算法能找出合理的工件呢？我們還得從目標上下手。
經(jīng)過我的一通查找，也看了一些論文，最后發(fā)現(xiàn)是有啟發(fā)性算法的，但是比較麻煩，而且效果也就好那么一丟丟，性價比太低。所以最后的結論就是：不需要使用啟發(fā)性算法，直接就進行鄰域隨機搜索。
你可能會懷疑，這和領域有什么關系呢？其實你只交換兩個，或者只交換一部分，這已經(jīng)是就近原則了。真正的全局純隨機，是利用np.random.shuffle函數(shù)將順序徹底打亂，我試過這個方法，十分影響表現(xiàn)，所以不能用。這里也看出了，我們的三種交換策略確實和全局純隨機不一樣。
基本函數(shù)結構
這里給出導入的包，以及要用的函數(shù)。
其中有三種交換策略，基礎的使用點交換，論文中還提到了區(qū)間顛倒和區(qū)間交換，我也都寫了，經(jīng)過測試保證有正確的結果。
# 解流水車間調度問題
# 函數(shù)編寫
# 讀取一個用例
def read():
with open('flowshop-test-10-student.txt') as file_obj:
#with open('one_instance.txt') as file_obj:
line=file_obj.readline()#去除頭兩行
line=file_obj.readline()
line=file_obj.readline().split()# 獲取n，m
n=int(line[0])
m=int(line[1])
data=[]# 逐行讀取內容
for i in range(n):
part=file_obj.readline().split()# 將工件數(shù)據(jù)提取出來
part=[int(x) for x in part]# 類型裝換
part=part[1::2]# 切片
data.append(part)
return data
# 時間到溫度
def time_to_temp(time):
return 0.99**time*100
# 動態(tài)規(guī)劃計算調度時間,schedule儲存順序
def all_time(data,schedule):
C=np.zeros((len(data),len(data[0])))
C[0][0]=data[schedule[0]][0]
for j in range(1,len(data[0])):# 計算1，m
C[0][j]=C[0][j-1]+data[schedule[0]][j]
for i in range(1,len(data)): # 計算n，1
C[i][0]=C[i-1][0]+data[schedule[i]][0]
for j in range(1,len(data[0])): # 填充
for i in range(1,len(data)):
C[i][j]=max(C[i-1][j],C[i][j-1])+data[schedule[i]][j]
#return C
return C[-1][-1]
# 鄰域搜索
def two_exchange(pre_schedule): # 將選定區(qū)域中工件順序顛倒
schedule=pre_schedule[:]
n=len(pre_schedule)
i=np.random.randint(0,n) # [low,high)，恰好覆蓋所有可能
j=np.random.randint(0,n)
if j<i: # 保持ij大小關系
temp=j
j=i
i=temp
for k in range(j-i+1):# 顛倒區(qū)間
schedule[j-k]=pre_schedule[i+k]
return schedule
def three_exchange(pre_schedule): # 交換兩段子序列位置
n=len(pre_schedule)
nodes=np.random.randint(0,n,3)
nodes=sorted(nodes)
print(nodes)
schedule=pre_schedule[:nodes[0]] # 兩邊要邊界，中間分給右邊
schedule.extend(pre_schedule[nodes[1]:nodes[2]+1])
schedule.extend(pre_schedule[nodes[0]:nodes[1]])
schedule.extend(pre_schedule[nodes[2]+1:])
return schedule
def point_exchange(pre_schedule):# 隨機兩點交換
schedule=pre_schedule[:]
n=len(schedule)
i=np.random.randint(0,n) # [low,high)，恰好覆蓋所有可能
j=np.random.randint(0,n)
#print(f"exchange {i},{j}")
t=schedule[i]
schedule[i]=schedule[j]
schedule[j]=t
return schedule
登山搜索框架
用例：
這是一個簡單的用例，用的時候保存成txt放在目錄下，自己調整一下代碼里的文件路徑就好，
1
1
11
11^{11}
1111的解空間也還不算特別大，勉強可以遍歷，實際后面我們會在50個工件級別的問題中計算，所以趁早死了遍歷的心。
instance 0
optimal-result:7038
11 5
0 375 1 12 2 142 3 245 4 412
0 632 1 452 2 758 3 278 4 398
0 12 1 876 2 124 3 534 4 765
0 460 1 542 2 523 3 120 4 499
0 528 1 101 2 789 3 124 4 999
0 796 1 245 2 632 3 375 4 123
0 532 1 230 2 543 3 896 4 452
0 14 1 124 2 214 3 543 4 785
0 257 1 527 2 753 3 210 4 463
0 896 1 896 2 214 3 258 4 259
0 532 1 302 2 501 3 765 4 988
代碼：
# 登山搜索
# 編碼和初始化
data=read()
print(data)
schedule=list(range(len(data)))# 初始順序調度
# 循環(huán)
t=0
T=time_to_temp(t)
epsilon=0.1 # 閾值溫度
min_time=[] # 記錄每個變化節(jié)點
min_t=[]
min_time.append(all_time(data,schedule))
min_t.append(0)
while T>epsilon:
# 多規(guī)則鄰域搜索生成新解，協(xié)同比較
new_schedule=point_exchange(schedule)
# 判斷新解情況，決定新狀態(tài)
now_time=all_time(data,new_schedule)
if(now_time<min_time[-1]): # 耗時更少
min_time.append(now_time)
min_t.append(t)
schedule=new_schedule
# 更新時間和溫度
t=t+1
T=time_to_temp(t)
min_time.append(min_time[-1])
min_t.append(t)
print(f"after {t} epoch",min_time[-1],schedule)
plt.xlabel('迭代次數(shù)')
plt.ylabel('調度總時間')
plt.plot(min_t,min_time)
一次結果，優(yōu)化到7038，有時候會有一點差距，這個和運氣有關系：
模擬退火框架
數(shù)據(jù)還是用的上面給出的數(shù)據(jù)，程序其實也僅僅是在對于新解判斷的一步進行改動。
# 模擬退火
# 編碼和初始化
data=read()
print(data)
schedule=list(range(len(data)))# 初始順序調度
# 循環(huán)
t=0
T=time_to_temp(t)
epsilon=0.1 # 閾值溫度
min_time=[] # 記錄每個變化節(jié)點
min_t=[]
min_time.append(all_time(data,schedule))
min_t.append(0)
while T>epsilon:
# 多規(guī)則鄰域搜索生成新解，協(xié)同比較
new_schedule=point_exchange(schedule)
# 判斷新解情況，決定新狀態(tài)
now_time=all_time(data,new_schedule)
if(now_time<min_time[-1]): # 耗時更少
min_time.append(now_time)
min_t.append(t)
schedule=new_schedule
else: # 概率接受更差結果
P=np.exp((min_time[-1]-now_time)/T)
if np.random.rand()<P: # 接受差結果
min_time.append(now_time)
min_t.append(t)
schedule=new_schedule
# 更新時間和溫度
t=t+1
T=time_to_temp(t)
min_time.append(min_time[-1])
min_t.append(t)
print(f"after {t} epoch",min_time[-1],schedule)
plt.xlabel('迭代次數(shù)')
plt.ylabel('調度總時間')
plt.plot(min_t,min_time)
跑出來的圖大概是這樣，可以看到波動還是比較大的，但是最后可以達到更好的效果。
完全版本項目文件
項目包含加強版源代碼，數(shù)據(jù)集，實驗結果圖，以及我寫的論文 。
百度云
提取碼：cyyy
上面給出的倆程序都很簡單，實際上我們還可以通過一些技巧，提高5%左右的表現(xiàn)。我自己就是從3500優(yōu)化到了3280的，基本上優(yōu)化了 220 3250 = 6.8 % \dfrac{220}{3250}=6.8\% 3250220?=6.8%，這么一看還是挺可觀的。
具體的改進：
- 程序架構。其實有11組數(shù)據(jù)，所以要用新的read函數(shù)，一次性讀取出11個矩陣，然后遍歷逐個處理。
- 10組隨機重啟搜索。每個測試用例都要進行10次隨機重啟搜索，互不干擾
- 增加記憶機制，選擇最好的組作為結果。這樣可以降低結果對出生點的依賴。
- 實現(xiàn)三種鄰域搜索，每次搜索都要從當前狀態(tài)進行三種搜索，選擇其中最好的結果，再進行判斷是否保留。
- 圖像繪制。采用1000dpi的清晰度，記錄每個優(yōu)化點的時間與結果，將10組都畫在一張圖里。最后保存到目錄里，用于論文撰寫。
看看我最后的成果，從圖中可以看出：
- 小規(guī)模問題上，登山搜索和模擬退火不會有太大差異
- 但是稍微上升到50個工件的規(guī)模，模擬退火就逐漸體現(xiàn)出優(yōu)勢，登山搜索10組才出來一個3280級別的，而模擬退火一眼看過去全是3280級別的，有理由認為登山搜索那個3280是偶然，也就是隨機重啟帶來的好處，而模擬退火即使是只跑一組實驗，結果也很穩(wěn)定
- 如果上升到100工件，在結果上模擬退火應該是會更優(yōu)的。
文章來源：http://www.zghlxwxcb.cn/news/detail-432375.html
文章來源地址http://www.zghlxwxcb.cn/news/detail-432375.html
到了這里，關于【數(shù)學建模】模擬退火全解析的文章就介紹完了。如果您還想了解更多內容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關文章，希望大家以后多多支持TOY模板網(wǎng)！
