模擬退火算法原理

????退火這個詞，其實是鐵匠發明的。它的意思很簡單，就是將鐵匠爐燒熱后，再把下邊的火撤掉，讓金屬在爐子里邊慢慢冷卻。人們發現，這個緩慢的降溫過程能消除金屬內部的各種缺陷，使得其恢復能量最低的狀態。后來，受到退火工藝的啟發，研究人員把將退火的緩慢的降溫思路推廣開來，用于尋找復雜函數的全局最優值。舉個簡單的例子，氦原子在金屬中空位附近運動時，其勢能函數大概長這樣：
????舉個簡單的例子，氦原子在金屬中空位附近運動時，其勢能函數大概長這樣：

????把勢能函數看成一個軌道，再把氦看成一個小球放到這個軌道上。那么，小球的高度就是它的重力勢能。求這個函數的全局最小值，其實就是讓小球滾到最深那個坑過程。

????由于函數中存在非常多局域極小值（淺坑），常見的貪心算法（如梯度下降法，就是閉著眼睛沿坡往下滾）會陷入最近的局域極小值，很錯過到全局最小值（深坑）。
????那么，怎么才能讓滾進深坑呢？
????我們先晃一晃整個體系，讓小球隨機動起來。我們把晃動的劇烈程度稱為溫度T。在足夠長的時間內。小球出現在i點的概率呈玻爾茲曼分布，與i點的高度(勢能Ei)和溫度T密切相關：

????溫度為0時，小球出現在全局最小值的概率為100%。
????那么，把溫度設為0不就行了？
????當然不行。上面的概率成立的前提是提供無窮長的時間，使得小球能夠遍歷整個函數。
????實際中，無窮長的時間當然是不現實的。小球從低概率區向高概率區轉移，這個過程是需要消耗時間的。
????高溫下，轉移速度快。但另一方面，高溫下，小球落在各個坑里的概率差別比較小。
????為了讓小球能夠盡快的轉移到最低點，我們可以采取一個降溫策略：

• 先設定一個較高的溫度T，使得小球有足夠的能量越過能壘，小球到處亂跑的概率比較大。
• 隨機移動小球的位置，檢查其移動前后的能量差 ΔE。
• 若 ΔE<0 ，小球能量降低，則接受此次移動。
• 若 ΔE>0，小球能量升高，有一定概率接受此次移動，概率為：exp(-ΔE/(K_BT)) 。這個判定方法稱為Metropolis判據。
• 逐漸降低溫度T，使得小球滑向低處的概率逐漸增大，直到小球緩緩收斂到最低點，不再晃動。
  
  

????通過這樣一個逐漸降溫的過程，小球最終有很大概率（并不是一定）能找到全局最低值。這個算法就是模擬退火算法。
????模擬退火算法最大的優勢在于其通用性：函數形式如何復雜我不管，我就一步步瞎跑，然后按Metropolis判據概率性的接受這一步就行。

流程圖見下：

---

算例

算例1

????計算函數f(x)=∑x^2(-20≤x≤20)的最小值，其中個體x的維數n= 10。 這是一個簡單的平方和函數，只有一個極小點x=(0, 0, …0)， 理論最小值f(0，0, .，. 0)=0

MATLAB版求解

%%%%%%%%%%%%%%%%%%%%%%模擬退火算法解決函數極值%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; %清除所有變量 close all; %清圖 clc; %清屏 D=10; %變量維數 Xs=20; %上限 Xx=-20; %下限 %%%%%%%%%%%%%%%%%%%%%%%%%%%冷卻表參數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% L = 200; %馬可夫鏈長度 K = 0.998; %衰減參數 S = 0.01; %步長因子 T=100; %初始溫度 YZ = 1e-8; %容差 P = 0; %Metropolis過程中總接受點 %%%%%%%%%%%%%%%%%%%%%%%%%%隨機選點 初值設定%%%%%%%%%%%%%%%%%%%%%%%%% PreX = rand(D,1)*(Xs-Xx)+Xx; PreBestX = PreX; PreX = rand(D,1)*(Xs-Xx)+Xx; BestX = PreX; %%%%%%%%%%%每迭代一次退火一次(降溫), 直到滿足迭代條件為止%%%%%%%%%%%% deta=abs( func1( BestX)-func1(PreBestX)); while (deta > YZ) && (T>0.001)T=K*T; %%%%%%%%%%%%%%%%%%%%%在當前溫度T下迭代次數%%%%%%%%%%%%%%%%%%%%%%for i=1:L %%%%%%%%%%%%%%%%%在此點附近隨機選下一點%%%%%%%%%%%%%%%%%%%%%NextX = PreX + S* (rand(D,1) *(Xs-Xx)+Xx);%%%%%%%%%%%%%%%%%邊界條件處理%%%%%%%%%%%%%%%%%%%%%%%%%%for ii=1:Dif NextX(ii)>Xs | NextX(ii)<XxNextX(ii)=PreX(ii) + S* (rand *(Xs-Xx)+Xx);endend %%%%%%%%%%%%%%%%%%%%%%%是否全局最優解%%%%%%%%%%%%%%%%%%%%%%if (func1(BestX) > func1(NextX))%%%%%%%%%%%%%%%%%%保留上一個最優解%%%%%%%%%%%%%%%%%%%%%PreBestX = BestX;%%%%%%%%%%%%%%%%%%%此為新的最優解%%%%%%%%%%%%%%%%%%%%%%BestX=NextX;end%%%%%%%%%%%%%%%%%%%%%%%% Metropolis過程%%%%%%%%%%%%%%%%%%%if( func1(PreX) - func1(NextX) > 0 )%%%%%%%%%%%%%%%%%%%%%%%接受新解%%%%%%%%%%%%%%%%%%%%%%%%PreX=NextX;P=P+1;elsechanger = -1*(func1(NextX)-func1(PreX))/ T ;p1=exp(changer);%%%%%%%%%%%%%%%%%%%%%%%%接受較差的解%%%%%%%%%%%%%%%%%%%%if p1 > rand PreX=NextX;P=P+1; endendtrace(P+1)=func1( BestX); enddeta=abs( func1( BestX)-func1 (PreBestX)); end disp('最小值在點:'); BestX disp( '最小值為:'); func1(BestX) figure plot(trace(2:end)) xlabel('迭代次數') ylabel('目標函數值') title('適應度進化曲線')%%%%%%%%%%%%%%%%%%%%%%%%%適應度函數%%%%%%%%%%%%%%%%%%%%%%% function result=func1(x) summ=sum(x.^2); result=summ; end

從圖中可以看到迭代次數太多，是個體尋優，需要花費很多時間。相比其他群體尋優方法，模擬退火算法很差。

python版求解

#!/usr/bin/env python3 # -*- coding: utf-8 -*- # @Author: yudengwu(余登武） # @Date : 2021/6/21 #@email:1344732766@qq.comimport numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import matplotlib; matplotlib.use('TkAgg') mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負號'-'顯示為方塊的問題#=================初始化參數============ D=10 #變量維數 Xs=20 #上限 Xx=-20 #下限 #==冷卻表參數== L = 200 #馬可夫鏈長度 #在溫度為t情況下的迭代次數 K = 0.998 #衰減參數 S = 0.01 #步長因子 T=100 #初始溫度 YZ = 1e-7 #容差 P = 0 #Metropolis過程中總接受點 #====隨機選點 初值設定==== PreX = np.random.uniform(size=(D,1))*(Xs-Xx)+Xx PreBestX = PreX#t-1代的全局最優X PreX = np.random.uniform(size=(D,1))*(Xs-Xx)+Xx BestX = PreX#t時刻的全局最優X#==============目標函數============= def func1(x):return np.sum([i**2 for i in x])#====每迭代一次退火一次(降溫), 直到滿足迭代條件為止=== deta=np.abs(func1(BestX)-func1(PreBestX))#前后能量差trace=[]#記錄 while (deta > YZ) and (T>0.1):#如果能量差大于允許能量差 或者溫度大于閾值T = K * T#降溫print(T)#=在當前溫度T下迭代次數=for i in range(L):##=在此點附近隨機選下一點==NextX = PreX + S * (np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx)#===邊界條件處理for ii in range(D):#遍歷每一個維度while NextX[ii]>Xs or NextX[ii]<Xx:NextX[ii]=PreX[ii] + S* (np.random.random() *(Xs-Xx)+Xx)#===是否全局最優解 ===if (func1(BestX) > func1(NextX)):#保留上一個最優解PreBestX = BestX#此為新的最優解BestX = NextX#====Metropolis過程===if (func1(PreX) - func1(NextX) > 0):#后一個比前一個好#接受新解PreX = NextXP = P + 1else:changer = -1 * (func1(NextX) - func1(PreX)) / Tp1 = np.exp(changer)#接受較差的解if p1 > np.random.random():PreX = NextXP = P + 1trace.append(func1(BestX))deta = np.abs(func1(BestX)- func1(PreBestX))#修改前后能量差print('最小值點\n',BestX) print('最小值\n',func1(BestX)) plt.plot(trace,label='迭代曲線') plt.xlabel('迭代次數') plt.ylabel('目標函數值') plt.show()

---

復雜約束求解：算例2

復雜約束求解：python版求解

在約束求解中，添加了懲罰項。

import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import matplotlib; matplotlib.use('TkAgg') mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負號'-'顯示為方塊的問題#=================初始化參數============ D=2 #變量維數 ,對應x,y Xs=1 #搜索區間上限 Xx=-1 #搜索區間下限 #==冷卻表參數== L = 300 #馬可夫鏈長度 #在溫度為t情況下的迭代次數 K = 0.98 #衰減參數 S = 0.01 #步長因子 T=100 #初始溫度 YZ = 1e-8 #容差(運行t-1,，t時刻的誤差） P = 0 #Metropolis過程中總接受點 eloss=0.1 #允許的懲罰項誤差 #====隨機選點 初值設定==== PreX = np.zeros(shape=(D,1)) PreX[0]=np.random.uniform(1,2,1)#X范圍【1,2】 PreX[1]=np.random.uniform(-1,0,1)#y范圍【-1,0】 PreBestX = PreX#t-1代的全局最優XPreX = np.zeros(shape=(D,1)) PreX[0]=np.random.uniform(1,2,1)#X范圍【1,2】 PreX[1]=np.random.uniform(-1,0,1)#y范圍【-1,0】 BestX = PreX#t時刻的全局最優X#==============目標函數============= def func1(X):A = 10pi = np.pix = X[0]y = X[1]return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y) #==========懲罰項=========== def calc_e(X):"""計算個體的目懲罰項，X 的維度是 size * 2 """ee = 0"""計算第一個約束的懲罰項"""e1 = X[0] + X[1] - 6ee += max(0, e1)"""計算第二個約束的懲罰項"""e2 = 3 * X[0] - 2 * X[1] - 5ee += max(0, e2)return ee#====每迭代一次退火一次(降溫), 直到滿足迭代條件為止=== deta=np.abs(func1(BestX)-func1(PreBestX))#前后能量差trace=[]#記錄 while (deta > YZ) and (T>0.1):#如果能量差大于允許能量差 或者溫度大于閾值T = K * T#降溫print(T)# =在當前溫度T下迭代次數=for i in range(L): ## =在此點附近隨機選下一點==NextX = PreX + S * (np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx)# ===邊界條件處理while NextX[0] > 2 or NextX[0] < 1:#x屬于【1,2】NextX[0] =PreX[0] + S* (np.random.uniform() * (Xs - Xx) + Xx)while NextX[1] > 0 or NextX[0] < -1: # y屬于【-1,0】NextX[1] = PreX[1] + S* (np.random.uniform() * (Xs - Xx) + Xx)# ===是否全局最優解 ===if (func1(BestX) > func1(NextX)) and (calc_e(BestX)<eloss) and (calc_e(NextX)<eloss):#兩者都沒有違反約束，下一代適應度更好時選下一代# 保留上一個最優解PreBestX = BestX# 此為新的最優解BestX = NextXelif (calc_e(BestX)>eloss) and (calc_e(NextX)>eloss) and (func1(BestX) > func1(NextX)):#兩者都違反約束，下一代適應度更好時（選下一代）# 保留上一個最優解PreBestX = BestX# 此為新的最優解BestX = NextXelif (calc_e(BestX)>eloss) and (calc_e(NextX)<eloss) :#上一代懲罰項不滿足約束，下一代懲罰項滿足約束（選下一代）# 保留上一個最優解PreBestX = BestX# 此為新的最優解BestX = NextX# ====Metropolis過程 當代溫度T 前后個體判斷===if (func1(PreX) - func1(NextX) > 0) and (calc_e(PreX)<eloss) and(calc_e(NextX)<eloss) : # 兩個懲罰項都滿足約束，后一個的適應度更小（選后一個）# 接受新解PreX = NextXP = P + 1elif (calc_e(PreX)>eloss) and(calc_e(NextX)>eloss) and (func1(PreX) - func1(NextX) > 0) :#兩個懲罰項都不滿足時，后一個適應度比前一個小（選后一個）# 接受新解PreX = NextXP = P + 1elif (calc_e(PreX)>eloss) and(calc_e(NextX)<eloss) :#前一個懲罰項不滿足，后一個滿足（選后一個）# 接受新解PreX = NextXP = P + 1else:changer = -1 * (func1(NextX) - func1(PreX)) / Tp1 = np.exp(changer)# 接受較差的解if p1 > np.random.random():PreX = NextXP = P + 1trace.append(func1(BestX))deta = np.abs(func1(BestX) - func1(PreBestX)) # 修改前后能量差print('最優值X\n',BestX) print('最優目標函數值\n',func1(BestX)) print('最優值對應的懲罰項',calc_e(BestX)) plt.plot(trace,label='迭代曲線') plt.xlabel('迭代次數') plt.ylabel('目標函數值') plt.show()

作者：電氣-余登武

總結

以上是生活随笔為你收集整理的万字长文了解模拟退火算法原理及求解复杂约束问题（源码实现）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。