1、模擬退火演算法
模擬退火演算法借鑑了統計物理學的思想，是一種簡單、通用的啟發式優化演算法，並在理論上具有概率性全域性優化效能，因而在科研和工程中得到了廣泛的應用。
退火是金屬從熔融狀態緩慢冷卻、最終達到能量最低的平衡態的過程。模擬退火演算法基於優化問題求解過程與金屬退火過程的相似性，以優化目標為能量函式，以解空間為狀態空間，以隨機擾動模擬粒子的熱運動來求解優化問題（[1] KIRKPATRICK,1988）。
模擬退火演算法結構簡單，由溫度更新函式、狀態產生函式、狀態接受函式和內迴圈、外迴圈終止準則構成。
溫度更新函式是指退火溫度緩慢降低的實現方案，也稱冷卻進度表；
狀態產生函式是指由當前解隨機產生新的候選解的方法；
狀態接受函式是指接受候選解的機制，通常採用Metropolis準則；
外迴圈是由冷卻進度表控制的溫度迴圈；
內迴圈是在每一溫度下迴圈迭代產生新解的次數，也稱Markov鏈長度。
模擬退火演算法的基本流程如下：
（1）初始化：初始溫度T，初始解狀態s，迭代次數L；
（2）對每個溫度狀態，重複 L次迴圈產生和概率性接受新解：
（3）通過變換操作由當前解s 產生新解s′；
（4）計算能量差 ∆E，即新解的目標函式與原有解的目標函式的差；
（5）若∆E <0則接受s′作為新的當前解，否則以概率exp(-∆E/T) 接受s′ 作為新的當前解；
（6）在每個溫度狀態完成 L次內迴圈後，降低溫度 T，直到達到終止溫度。
2、多變數函式優化問題
選取經典的函式優化問題和組合優化問題作為測試案例。
問題 1：Schwefel 測試函式，是複雜的多峰函式，具有大量區域性極值區域。
F(X)=418.9829×n-∑(i=1,n)〖xi* sin(√(|xi|)) 〗
本文取 d=10, x=[-500,500]，函式在 X=(420.9687,...420.9687)處為全域性最小值 f(X)=0.0。
使用模擬退火演算法的基本方案：控制溫度按照 T(k) = a * T(k-1) 指數衰減，衰減係數取 a；如式（1）按照 Metropolis 準則接受新解。對於問題 1（Schwefel函式），通過對當前解的一個自變數施加正態分佈的隨機擾動產生新解。
3、模擬退火演算法 Python 程式
# 模擬退火演算法 程式：多變數連續函式優化
# Program: SimulatedAnnealing_v1.py
# Purpose: Simulated annealing algorithm for function optimization
# v1.0:
# (1) 基本演算法：單變數連續函式優化問題
# (2) 檔案輸出優化結果和中間過程資料
# (3) 設定指標引數計數器
# (4) 圖形輸出壞解接受概率
# Copyright 2021 YouCans, XUPT
# Crated：2021-04-30
# -*- coding: utf-8 -*-
import math # 匯入模組
import random # 匯入模組
import pandas as pd # 匯入模組 YouCans, XUPT
import numpy as np # 匯入模組 numpy, 並簡寫成 np
import matplotlib.pyplot as plt # 匯入模組 matplotlib.pyplot, 並簡寫成 plt
from datetime import datetime
# 子程式：定義優化問題的目標函式
def cal_Energy(X, nVar):
# 測試函式 1： Schwefel 測試函式
# -500 <= Xi <= 500
# 全域性極值：(420.9687,420.9687,...),f(x)=0.0
sum = 0.0
for i in range(nVar):
sum += X[i] * np.sin(np.sqrt(abs(X[i])))
fx = 418.9829 * nVar - sum
return fx
# 子程式：模擬退火演算法的引數設定
def ParameterSetting():
cName = "funcOpt" # 定義問題名稱
nVar = 2 # 給定自變數數量，y=f(x1,..xn)
xMin = [-500, -500] # 給定搜尋空間的下限，x1_min,..xn_min
xMax = [500, 500] # 給定搜尋空間的上限，x1_max,..xn_max
tInitial = 100.0 # 設定初始退火溫度(initial temperature)
tFinal = 1 # 設定終止退火溫度(stop temperature)
alfa = 0.98 # 設定降溫引數，T(k)=alfa*T(k-1)
meanMarkov = 100 # Markov鏈長度，也即內迴圈執行次數
scale = 0.5 # 定義搜尋步長，可以設為固定值或逐漸縮小
return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale
# 模擬退火演算法
def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale):
# ====== 初始化隨機數發生器 ======
randseed = random.randint(1, 100)
random.seed(randseed) # 隨機數發生器設定種子，也可以設為指定整數
# ====== 隨機產生優化問題的初始解 ======
xInitial = np.zeros((nVar)) # 初始化，建立陣列
for v in range(nVar):
# random.uniform(min,max) 在 [min,max] 範圍內隨機生成一個實數
xInitial[v] = random.uniform(xMin[v], xMax[v])
# 呼叫子函式 cal_Energy 計算當前解的目標函式值
fxInitial = cal_Energy(xInitial, nVar)
# ====== 模擬退火演算法初始化 ======
xNew = np.zeros((nVar)) # 初始化，建立陣列
xNow = np.zeros((nVar)) # 初始化，建立陣列
xBest = np.zeros((nVar)) # 初始化，建立陣列
xNow[:] = xInitial[:] # 初始化當前解，將初始解置為當前解
xBest[:] = xInitial[:] # 初始化最優解，將當前解置為最優解
fxNow = fxInitial # 將初始解的目標函式置為當前值
fxBest = fxInitial # 將當前解的目標函式置為最優值
print('x_Initial:{:.6f},{:.6f},\tf(x_Initial):{:.6f}'.format(xInitial[0], xInitial[1], fxInitial))
recordIter = [] # 初始化，外迴圈次數
recordFxNow = [] # 初始化，當前解的目標函式值
recordFxBest = [] # 初始化，最佳解的目標函式值
recordPBad = [] # 初始化，劣質解的接受概率
kIter = 0 # 外迴圈迭代次數，溫度狀態數
totalMar = 0 # 總計 Markov 鏈長度
totalImprove = 0 # fxBest 改善次數
nMarkov = meanMarkov # 固定長度 Markov鏈
# ====== 開始模擬退火優化 ======
# 外迴圈，直到當前溫度達到終止溫度時結束
tNow = tInitial # 初始化當前溫度(current temperature)
while tNow >= tFinal: # 外迴圈，直到當前溫度達到終止溫度時結束
# 在當前溫度下，進行充分次數(nMarkov)的狀態轉移以達到熱平衡
kBetter = 0 # 獲得優質解的次數
kBadAccept = 0 # 接受劣質解的次數
kBadRefuse = 0 # 拒絕劣質解的次數
# ---內迴圈，迴圈次數為Markov鏈長度
for k in range(nMarkov): # 內迴圈，迴圈次數為Markov鏈長度
totalMar += 1 # 總 Markov鏈長度計數器
# ---產生新解
# 產生新解：通過在當前解附近隨機擾動而產生新解，新解必須在 [min,max] 範圍內
# 方案 1：只對 n元變數中的一個進行擾動，其它 n-1個變數保持不變
xNew[:] = xNow[:]
v = random.randint(0, nVar-1) # 產生 [0,nVar-1]之間的隨機數
xNew[v] = xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1)
# random.normalvariate(0, 1)：產生服從均值為0、標準差為 1 的正態分佈隨機實數
xNew[v] = max(min(xNew[v], xMax[v]), xMin[v]) # 保證新解在 [min,max] 範圍內
# ---計算目標函式和能量差
# 呼叫子函式 cal_Energy 計算新解的目標函式值
fxNew = cal_Energy(xNew, nVar)
deltaE = fxNew - fxNow
# ---按 Metropolis 準則接受新解
# 接受判別：按照 Metropolis 準則決定是否接受新解
if fxNew < fxNow: # 更優解：如果新解的目標函式好於當前解，則接受新解
accept = True
kBetter += 1
else: # 容忍解：如果新解的目標函式比當前解差，則以一定概率接受新解
pAccept = math.exp(-deltaE / tNow) # 計算容忍解的狀態遷移概率
if pAccept > random.random():
accept = True # 接受劣質解
kBadAccept += 1
else:
accept = False # 拒絕劣質解
kBadRefuse += 1
# 儲存新解
if accept == True: # 如果接受新解，則將新解儲存為當前解
xNow[:] = xNew[:]
fxNow = fxNew
if fxNew < fxBest: # 如果新解的目標函式好於最優解，則將新解儲存為最優解
fxBest = fxNew
xBest[:] = xNew[:]
totalImprove += 1
scale = scale*0.99 # 可變搜尋步長，逐步減小搜尋範圍，提高搜尋精度
# ---內迴圈結束後的資料整理
# 完成當前溫度的搜尋，儲存資料和輸出
pBadAccept = kBadAccept / (kBadAccept + kBadRefuse) # 劣質解的接受概率
recordIter.append(kIter) # 當前外迴圈次數
recordFxNow.append(round(fxNow, 4)) # 當前解的目標函式值
recordFxBest.append(round(fxBest, 4)) # 最佳解的目標函式值
recordPBad.append(round(pBadAccept, 4)) # 最佳解的目標函式值
if kIter%10 == 0: # 模運算，商的餘數
print('i:{},t(i):{:.2f}, badAccept:{:.6f}, f(x)_best:{:.6f}'.\
format(kIter, tNow, pBadAccept, fxBest))
# 緩慢降溫至新的溫度，降溫曲線：T(k)=alfa*T(k-1)
tNow = tNow * alfa
kIter = kIter + 1
# ====== 結束模擬退火過程 ======
print('improve:{:d}'.format(totalImprove))
return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad
# 結果校驗與輸出
def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter):
# ====== 優化結果校驗與輸出 ======
fxCheck = cal_Energy(xBest,nVar)
if abs(fxBest - fxCheck)>1e-3: # 檢驗目標函式
print("Error 2: Wrong total millage!")
return
else:
print("\nOptimization by simulated annealing algorithm:")
for i in range(nVar):
print('\tx[{}] = {:.6f}'.format(i,xBest[i]))
print('\n\tf(x):{:.6f}'.format(fxBest))
# ====== 優化結果寫入資料檔案 ======
nowTime = datetime.now().strftime('%m%d%H%M') # '02151456'
fileName = "..\data\{}_{}.dat".format(cName,nowTime)# 資料檔案的地址和檔名
optRecord = {
"iter":recordIter,
"FxNow":recordFxNow,
"FxBest":recordFxBest,
"PBad":recordPBad}
df_Record = pd.DataFrame(optRecord)
df_Record.to_csv(fileName, index=False, encoding="utf_8_sig")
with open(fileName, 'a+', encoding="utf_8_sig") as fid:
fid.write("\nOptimization by simulated annealing algorithm:")
for i in range(nVar):
fid.write('\n\tx[{}] = {:.6f}'.format(i,xBest[i]))
fid.write('\n\tf(x):{:.6f}'.format(fxBest))
print("寫入資料檔案: %s 完成。" % fileName)
# ====== 優化結果圖形化輸出 ======
plt.figure(figsize=(6, 4), facecolor='#FFFFFF') # 建立一個圖形視窗
plt.title('Optimization result: {}'.format(cName)) # 設定圖形標題
plt.xlim((0, kIter)) # 設定 x軸範圍
plt.xlabel('iter') # 設定 x軸標籤
plt.ylabel('f(x)') # 設定 y軸標籤
plt.plot(recordIter, recordFxNow,'b-', label='FxNow') # 繪製 FxNow 曲線
plt.plot(recordIter, recordFxBest, 'r-', label='FxBest') # 繪製 FxBest 曲線
# plt.plot(recordIter,recordPBad,'r-',label='pBadAccept') # 繪製 pBadAccept 曲線
plt.legend() # 顯示圖例
plt.show()
return
# 主程式
def main():
# 引數設定，優化問題引數定義，模擬退火演算法引數設定
[cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale] = ParameterSetting()
# print([nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale])
# 模擬退火演算法
[kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad] \
= OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale)
# print(kIter, fxNow, fxBest, pBadAccept)
# 結果校驗與輸出
ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter)
if __name__ == '__main__':
main()
版權說明：
原創作品
Copyright 2021 YouCans, XUPT
Crated：2021-05-01