Python數(shù)學(xué)建模學(xué)習(xí)模擬退火算法多變量函數(shù)優(yōu)化示例解析
更新時間：2021年10月18日 11:42:46 作者：youcans
模擬退火算法借鑒了統(tǒng)計物理學(xué)的思想，是一種簡單、通用的啟發(fā)式優(yōu)化算法，并在理論上具有概率性全局優(yōu)化性能，因而在科研和工程中得到了廣泛的應(yīng)用
1、模擬退火算法
退火是金屬從熔融狀態(tài)緩慢冷卻、最終達到能量最低的平衡態(tài)的過程。模擬退火算法基于優(yōu)化問題求解過程與金屬退火過程的相似性，以優(yōu)化目標為能量函數(shù)，以解空間為狀態(tài)空間，以隨機擾動模擬粒子的熱運動來求解優(yōu)化問題（[1] KIRKPATRICK,1988）。
模擬退火算法結(jié)構(gòu)簡單，由溫度更新函數(shù)、狀態(tài)產(chǎn)生函數(shù)、狀態(tài)接受函數(shù)和內(nèi)循環(huán)、外循環(huán)終止準則構(gòu)成。
溫度更新函數(shù)是指退火溫度緩慢降低的實現(xiàn)方案，也稱冷卻進度表；
狀態(tài)產(chǎn)生函數(shù)是指由當(dāng)前解隨機產(chǎn)生新的候選解的方法；
狀態(tài)接受函數(shù)是指接受候選解的機制，通常采用Metropolis準則；
外循環(huán)是由冷卻進度表控制的溫度循環(huán)；
內(nèi)循環(huán)是在每一溫度下循環(huán)迭代產(chǎn)生新解的次數(shù)，也稱Markov鏈長度。
模擬退火算法的基本流程如下：
（1）初始化：初始溫度T，初始解狀態(tài)s，迭代次數(shù)L；
（2）對每個溫度狀態(tài)，重復(fù) L次循環(huán)產(chǎn)生和概率性接受新解：
（3）通過變換操作由當(dāng)前解s 產(chǎn)生新解s′；
（4）計算能量差 ∆E，即新解的目標函數(shù)與原有解的目標函數(shù)的差；
（5）若∆E <0則接受s′作為新的當(dāng)前解，否則以概率exp(-∆E/T) 接受s′ 作為新的當(dāng)前解；
（6）在每個溫度狀態(tài)完成 L次內(nèi)循環(huán)后，降低溫度 T，直到達到終止溫度。
2、多變量函數(shù)優(yōu)化問題
選取經(jīng)典的函數(shù)優(yōu)化問題和組合優(yōu)化問題作為測試案例。
問題 1：Schwefel 測試函數(shù)，是復(fù)雜的多峰函數(shù)，具有大量局部極值區(qū)域。
本文取 d=10, x=[-500,500]，函數(shù)在 X=(420.9687,…420.9687)處為全局最小值 f(X)=0.0。
F(X)=418.9829×n-∑(i=1,n)〖xi* sin(√(|xi|)) 〗
使用模擬退火算法的基本方案：控制溫度按照 T(k) = a * T(k-1) 指數(shù)衰減，衰減系數(shù)取 a；如式（1）按照 Metropolis 準則接受新解。對于問題 1（Schwefel函數(shù)），通過對當(dāng)前解的一個自變量施加正態(tài)分布的隨機擾動產(chǎn)生新解。
3、模擬退火算法 Python 程序
# 模擬退火算法 程序：多變量連續(xù)函數(shù)優(yōu)化
# Program: SimulatedAnnealing_v1.py
# Purpose: Simulated annealing algorithm for function optimization
# Copyright 2021 YouCans, XUPT
# Crated：2021-04-30
# = 關(guān)注 Youcans，分享原創(chuàng)系列 https://blog.csdn.net/youcans =
# -*- coding: utf-8 -*-
import math # 導(dǎo)入模塊
import random # 導(dǎo)入模塊
import pandas as pd # 導(dǎo)入模塊
import numpy as np # 導(dǎo)入模塊 numpy，并簡寫成 np
import matplotlib.pyplot as plt # 導(dǎo)入模塊 matplotlib.pyplot，并簡寫成 plt
from datetime import datetime
# 子程序：定義優(yōu)化問題的目標函數(shù)
def cal_Energy(X, nVar):
# 測試函數(shù) 1： Schwefel 測試函數(shù)
# -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
# 子程序：模擬退火算法的參數(shù)設(shè)置
def ParameterSetting():
cName = "funcOpt" # 定義問題名稱
nVar = 2 # 給定自變量數(shù)量，y=f(x1,..xn)
xMin = [-500, -500] # 給定搜索空間的下限，x1_min,..xn_min
xMax = [500, 500] # 給定搜索空間的上限，x1_max,..xn_max
tInitial = 100.0 # 設(shè)定初始退火溫度(initial temperature)
tFinal = 1 # 設(shè)定終止退火溫度(stop temperature)
alfa = 0.98 # 設(shè)定降溫參數(shù)，T(k)=alfa*T(k-1)
meanMarkov = 100 # Markov鏈長度，也即內(nèi)循環(huán)運行次數(shù)
scale = 0.5 # 定義搜索步長，可以設(shè)為固定值或逐漸縮小
return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale
# 模擬退火算法
def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale):
# ====== 初始化隨機數(shù)發(fā)生器 ======
randseed = random.randint(1, 100)
random.seed(randseed) # 隨機數(shù)發(fā)生器設(shè)置種子，也可以設(shè)為指定整數(shù)
# ====== 隨機產(chǎn)生優(yōu)化問題的初始解 ======
xInitial = np.zeros((nVar)) # 初始化，創(chuàng)建數(shù)組
for v in range(nVar):
# random.uniform(min,max) 在 [min,max] 范圍內(nèi)隨機生成一個實數(shù)
xInitial[v] = random.uniform(xMin[v], xMax[v])
# 調(diào)用子函數(shù) cal_Energy 計算當(dāng)前解的目標函數(shù)值
fxInitial = cal_Energy(xInitial, nVar)
# ====== 模擬退火算法初始化 ======
xNew = np.zeros((nVar)) # 初始化，創(chuàng)建數(shù)組
xNow = np.zeros((nVar)) # 初始化，創(chuàng)建數(shù)組
xBest = np.zeros((nVar)) # 初始化，創(chuàng)建數(shù)組
xNow[:] = xInitial[:] # 初始化當(dāng)前解，將初始解置為當(dāng)前解
xBest[:] = xInitial[:] # 初始化最優(yōu)解，將當(dāng)前解置為最優(yōu)解
fxNow = fxInitial # 將初始解的目標函數(shù)置為當(dāng)前值
fxBest = fxInitial # 將當(dāng)前解的目標函數(shù)置為最優(yōu)值
print('x_Initial:{:.6f},{:.6f},\tf(x_Initial):{:.6f}'.format(xInitial[0], xInitial[1], fxInitial))
recordIter = [] # 初始化，外循環(huán)次數(shù)
recordFxNow = [] # 初始化，當(dāng)前解的目標函數(shù)值
recordFxBest = [] # 初始化，最佳解的目標函數(shù)值
recordPBad = [] # 初始化，劣質(zhì)解的接受概率
kIter = 0 # 外循環(huán)迭代次數(shù)，溫度狀態(tài)數(shù)
totalMar = 0 # 總計 Markov 鏈長度
totalImprove = 0 # fxBest 改善次數(shù)
nMarkov = meanMarkov # 固定長度 Markov鏈
# ====== 開始模擬退火優(yōu)化 ======
# 外循環(huán)，直到當(dāng)前溫度達到終止溫度時結(jié)束
tNow = tInitial # 初始化當(dāng)前溫度(current temperature)
while tNow >= tFinal: # 外循環(huán)，直到當(dāng)前溫度達到終止溫度時結(jié)束
# 在當(dāng)前溫度下，進行充分次數(shù)(nMarkov)的狀態(tài)轉(zhuǎn)移以達到熱平衡
kBetter = 0 # 獲得優(yōu)質(zhì)解的次數(shù)
kBadAccept = 0 # 接受劣質(zhì)解的次數(shù)
kBadRefuse = 0 # 拒絕劣質(zhì)解的次數(shù)
# ---內(nèi)循環(huán)，循環(huán)次數(shù)為Markov鏈長度
for k in range(nMarkov): # 內(nèi)循環(huán)，循環(huán)次數(shù)為Markov鏈長度
totalMar += 1 # 總 Markov鏈長度計數(shù)器
# ---產(chǎn)生新解
# 產(chǎn)生新解：通過在當(dāng)前解附近隨機擾動而產(chǎn)生新解，新解必須在 [min,max] 范圍內(nèi)
# 方案 1：只對 n元變量中的一個進行擾動，其它 n-1個變量保持不變
xNew[:] = xNow[:]
v = random.randint(0, nVar-1) # 產(chǎn)生 [0,nVar-1]之間的隨機數(shù)
xNew[v] = xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1)
# random.normalvariate(0, 1)：產(chǎn)生服從均值為0、標準差為 1 的正態(tài)分布隨機實數(shù)
xNew[v] = max(min(xNew[v], xMax[v]), xMin[v]) # 保證新解在 [min,max] 范圍內(nèi)
# ---計算目標函數(shù)和能量差
# 調(diào)用子函數(shù) cal_Energy 計算新解的目標函數(shù)值
fxNew = cal_Energy(xNew, nVar)
deltaE = fxNew - fxNow
# ---按 Metropolis 準則接受新解
# 接受判別：按照 Metropolis 準則決定是否接受新解
if fxNew < fxNow: # 更優(yōu)解：如果新解的目標函數(shù)好于當(dāng)前解，則接受新解
accept = True
kBetter += 1
else: # 容忍解：如果新解的目標函數(shù)比當(dāng)前解差，則以一定概率接受新解
pAccept = math.exp(-deltaE / tNow) # 計算容忍解的狀態(tài)遷移概率
if pAccept > random.random():
accept = True # 接受劣質(zhì)解
kBadAccept += 1
else:
accept = False # 拒絕劣質(zhì)解
kBadRefuse += 1
# 保存新解
if accept == True: # 如果接受新解，則將新解保存為當(dāng)前解
xNow[:] = xNew[:]
fxNow = fxNew
if fxNew < fxBest: # 如果新解的目標函數(shù)好于最優(yōu)解，則將新解保存為最優(yōu)解
fxBest = fxNew
xBest[:] = xNew[:]
totalImprove += 1
scale = scale*0.99 # 可變搜索步長，逐步減小搜索范圍，提高搜索精度
# ---內(nèi)循環(huán)結(jié)束后的數(shù)據(jù)整理
# 完成當(dāng)前溫度的搜索，保存數(shù)據(jù)和輸出
pBadAccept = kBadAccept / (kBadAccept + kBadRefuse) # 劣質(zhì)解的接受概率
recordIter.append(kIter) # 當(dāng)前外循環(huán)次數(shù)
recordFxNow.append(round(fxNow, 4)) # 當(dāng)前解的目標函數(shù)值
recordFxBest.append(round(fxBest, 4)) # 最佳解的目標函數(shù)值
recordPBad.append(round(pBadAccept, 4)) # 最佳解的目標函數(shù)值
if kIter%10 == 0: # 模運算，商的余數(shù)
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
# ====== 結(jié)束模擬退火過程 ======
print('improve:{:d}'.format(totalImprove))
return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad
# 結(jié)果校驗與輸出
def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter):
# ====== 優(yōu)化結(jié)果校驗與輸出 ======
fxCheck = cal_Energy(xBest,nVar)
if abs(fxBest - fxCheck)>1e-3: # 檢驗?zāi)繕撕瘮?shù)
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))
return
# 加粗樣式
def main():
# 參數(shù)設(shè)置，優(yōu)化問題參數(shù)定義，模擬退火算法參數(shù)設(shè)置
[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)
# 結(jié)果校驗與輸出
ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter)
if __name__ == '__main__':
main()
4、程序運行結(jié)果
x_Initial:-143.601793,331.160277, f(x_Initial):959.785447 i:0,t(i):100.00, badAccept:0.469136, f(x)_best:300.099320 i:10,t(i):81.71, badAccept:0.333333, f(x)_best:12.935760 i:20,t(i):66.76, badAccept:0.086022, f(x)_best:2.752498 ... i:200,t(i):1.76, badAccept:0.000000, f(x)_best:0.052055 i:210,t(i):1.44, badAccept:0.000000, f(x)_best:0.009448 i:220,t(i):1.17, badAccept:0.000000, f(x)_best:0.009448 improve:18 Optimization by simulated annealing algorithm: x[0] = 420.807471 x[1] = 420.950005 f(x):0.003352
以上就是Python數(shù)學(xué)建模學(xué)習(xí)模擬退火算法多變量函數(shù)優(yōu)化示例解析的詳細內(nèi)容，更多關(guān)于數(shù)學(xué)建模模擬退火算法的資料請關(guān)注腳本之家其它相關(guān)文章！
您可能感興趣的文章:
相關(guān)文章
python實現(xiàn)飛機大戰(zhàn)游戲(pygame版)
這篇文章主要為大家詳細介紹了python實現(xiàn)pygame版的飛機大戰(zhàn)游戲，文中示例代碼介紹的非常詳細，具有一定的參考價值，感興趣的小伙伴們可以參考一下2020-02-02
Python 的 httpx 包是一個用于 HTTP 交互的一個優(yōu)秀且靈活的模塊。本文進行詳細的講解，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步，早日升職加薪2022-03-03
Python工程實踐之np.loadtxt()讀取數(shù)據(jù)
為了方便使用和記憶,有時候我們會把numpy.loadtxt()縮寫成np.loadtxt(),這篇文章主要給大家介紹了關(guān)于Python工程實踐之np.loadtxt()讀取數(shù)據(jù)的相關(guān)資料,需要的朋友可以參考下2023-07-07
Windows系統(tǒng)下pycharm中的pip換源
這篇文章主要介紹了Windows系統(tǒng)下pycharm中的pip換源,文中通過示例代碼介紹的非常詳細，對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2020-02-02
在本篇文章里小編給大家整理的是一篇關(guān)于python文件路徑操作方法總結(jié)內(nèi)容，有需要的朋友們可以學(xué)習(xí)下。2020-12-12
本文主要介紹了Python?json模塊常用方法小結(jié)，文中通過示例代碼介紹的非常詳細，對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2023-03-03
這篇文章主要為大家介紹了Python dict的使用誤區(qū)，具有一定的參考價值，感興趣的小伙伴們可以參考一下，希望能夠給你帶來幫助2022-01-01
這篇文章主要介紹了Python面向?qū)ο缶幊讨惖睦^承，繼承Inheritance是指代碼復(fù)用的高級抽象，繼承是面向?qū)ο笤O(shè)計的精髓之一，實現(xiàn)了以類為單位的高級抽象級別代碼復(fù)用，下面進入文章看該內(nèi)容的下詳情2021-11-11
這篇文章主要為大家介紹了Python進行有限元仿真的創(chuàng)建及使用，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步，早日升職加薪2023-10-10
下面小編就為大家分享一篇pandas數(shù)據(jù)分組和聚合操作方法，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧2018-04-04