Python數(shù)學(xué)建模學(xué)習(xí)模擬退火算法約束條件處理示例解析
1、最優(yōu)化與線性規(guī)劃
最優(yōu)化問(wèn)題的三要素是決策變量、目標(biāo)函數(shù)和約束條件。
線性規(guī)劃（Linear programming），是研究線性約束條件下線性目標(biāo)函數(shù)的極值問(wèn)題的優(yōu)化方法，常用于解決利用現(xiàn)有的資源得到最優(yōu)決策的問(wèn)題。
簡(jiǎn)單的線性規(guī)劃問(wèn)題可以用 Lingo軟件求解，Matlab、Python 中也有求解線性規(guī)劃問(wèn)題的庫(kù)函數(shù)或求解器，很容易學(xué)習(xí)和使用，并不需要用模擬退火算法。但是，由一般線性規(guī)劃問(wèn)題所衍生的整數(shù)規(guī)劃、混合規(guī)劃、0/1規(guī)劃、二次規(guī)劃、非線性規(guī)劃、組合優(yōu)化問(wèn)題，則并不是調(diào)用某個(gè)庫(kù)函數(shù)都能處理的。而模擬退火算法在很多復(fù)雜問(wèn)題中具有較好的適應(yīng)性，可以作為一種入門(mén)的通用智能算法來(lái)學(xué)習(xí)。
也就是說(shuō)，如果只是處理線性規(guī)劃問(wèn)題，就不要用模擬退火算法了。但如果是現(xiàn)有方法無(wú)法處理的復(fù)雜優(yōu)化問(wèn)題，或者對(duì)某類(lèi)、某個(gè)優(yōu)化問(wèn)題你不知道用什么方法處理了，這時(shí)用模擬退火算法還是能解決的。
本文使用懲罰函數(shù)法，分析模擬退火算法處理線性規(guī)劃問(wèn)題，相關(guān)內(nèi)容也適用于非線性規(guī)劃問(wèn)題。
2、模擬退火算法處理約束條件
線性規(guī)劃問(wèn)題是約束優(yōu)化問(wèn)題，而模擬退火算法則更適合處理無(wú)約束優(yōu)化問(wèn)題。對(duì)于優(yōu)化問(wèn)題中的約束條件，模擬退火算法有幾種常用的處理方法：
1.決策變量取值的上下限約束。
此類(lèi)約束條件比較容易處理，只要設(shè)定初始解、新解在決策變量取值的上下限之間就可以解決。例如：
（1）設(shè)置產(chǎn)生新解的隨機(jī)數(shù)的上下限為決策變量的上下限，即 [Xmin, Xmax]；
（2）設(shè)置產(chǎn)生新解的隨機(jī)數(shù)的上下限為當(dāng)前解與決策變量的上下限，即 [Xnow, Xmax]；
（3）通過(guò)條件判斷，當(dāng)新解超出決策變量上下限，則令其取上下限，即 xNew = max(min(xNew, xMax), xMin)。當(dāng)然，這些處理方式，都會(huì)影響隨機(jī)數(shù)的概率分布，因而也影響模擬退火算法的優(yōu)化性能，在此不做深入討論。
2.檢驗(yàn)法處理不等式約束問(wèn)題。
在模擬退火算法的迭代過(guò)程中，將每次產(chǎn)生的新解代入每個(gè)不等式約束函數(shù)，判斷是否滿足約束條件；如果新解不滿足約束條件，則舍棄這個(gè)新解，返回重新產(chǎn)生一個(gè)新解進(jìn)行檢驗(yàn)，直到產(chǎn)生的新解滿足全部約束條件為止。這個(gè)方法的思路簡(jiǎn)單，每次迭代都在可行域內(nèi)進(jìn)行，但是對(duì)于約束條件眾多、苛刻的復(fù)雜問(wèn)題，多次產(chǎn)生的新解都不能滿足約束條件，會(huì)使計(jì)算時(shí)間很長(zhǎng)，甚至停滯不前。
3.消元法處理等式約束問(wèn)題。
對(duì)于等式約束，很難通過(guò)隨機(jī)產(chǎn)生的新解滿足約束條件，通常不能使用檢驗(yàn)法處理。消元法是通過(guò)解方程將等式約束中的某個(gè)決策變量表示為其它決策變量的線性關(guān)系后，代入目標(biāo)函數(shù)和不等式約束條件中，從而消去該約束條件。消元法不僅解決了等式約束問(wèn)題，而且減少了決策變量的數(shù)量，從而有效簡(jiǎn)化了優(yōu)化問(wèn)題的復(fù)雜度，是一舉兩得的處理方法。但是，對(duì)于非線性等式約束，求解非線性方程組也是非常困難的，消元法并不是普遍都能適用的。
4.更為通用的處理約束條件的方法是懲罰函數(shù)法，以下進(jìn)行介紹。YouCans, XUPT
3、懲罰函數(shù)法
懲罰函數(shù)法是一類(lèi)常用的處理約束條件的技術(shù)，在模擬退火算法中處理約束條件非常有效。方法的思想是將約束條件轉(zhuǎn)化為懲罰函數(shù)，附加在原有的目標(biāo)函數(shù)上構(gòu)造新的目標(biāo)函數(shù)；當(dāng)不滿足約束條件時(shí)，通過(guò)懲罰函數(shù)使新的目標(biāo)函數(shù)變差而被舍棄。
懲罰函數(shù)法有外點(diǎn)法和內(nèi)點(diǎn)法。外點(diǎn)法對(duì)可行域外的點(diǎn)（即不滿足約束的點(diǎn)）施加懲罰，對(duì)可行域內(nèi)部的點(diǎn)不懲罰，從而使迭代點(diǎn)向可行域D逼近。 內(nèi)點(diǎn)法是在可行域內(nèi)部進(jìn)行搜索，約束邊界起到類(lèi)似圍墻的作用，使目標(biāo)函數(shù)無(wú)法穿過(guò)，就把搜索點(diǎn)限制在可行域內(nèi)了，因此只適用于不等式約束。
4、數(shù)模案例
雖然對(duì)于線性規(guī)劃問(wèn)題并不推薦使用模擬退火算法求解。但為了便于理解，本文仍使用之前的線性規(guī)劃問(wèn)題作為處理約束條件的案例。對(duì)于非線性規(guī)劃問(wèn)題，以及非線性約束問(wèn)題，處理方法都是類(lèi)似的，將在后續(xù)進(jìn)行介紹。
4.1 問(wèn)題描述：
某廠生產(chǎn)甲乙兩種飲料，每百箱甲飲料需用原料6千克、工人10名，獲利10萬(wàn)元；每百箱乙飲料需用原料5千克、工人20名，獲利9萬(wàn)元。
今工廠共有原料60千克、工人150名，又由于其他條件所限甲飲料產(chǎn)量不超過(guò)8百箱。
（1）問(wèn)如何安排生產(chǎn)計(jì)劃，即兩種飲料各生產(chǎn)多少使獲利最大？
4.2 問(wèn)題建模：
決策變量：
x1：甲飲料產(chǎn)量（單位：百箱）
x2：乙飲料產(chǎn)量（單位：百箱）
目標(biāo)函數(shù)：
max fx = 10x1 + 9x2
約束條件：
6x1 + 5x2 <= 60
10x1 + 20x2 <= 150
取值范圍：
給定條件：x1, x2 >= 0，x1 <= 8
推導(dǎo)條件：由 x1,x2>=0 和 10*x1+20*x2<=150 可知：0<=x1<=15；0<=x2<=7.5
因此，0 <= x1<=8，0 <= x2<=7.5
4.3 懲罰函數(shù)法求解約束優(yōu)化問(wèn)題：
構(gòu)造懲罰函數(shù)：
p1 = (max(0, 6*x1+5*x2-60))**2
p2 = (max(0, 10*x1+20*x2-150))**2
說(shuō)明：如存在等式約束，例如：x1 + 2*x2 = m，也可以轉(zhuǎn)化為懲罰函數(shù)：
p3 = (x1+2*x2-m)**2
P(x) = p1 + p2 + …
構(gòu)造增廣目標(biāo)函數(shù)：
L(x,m(k)) = min(fx) + m(k)*P(x)
m(k)：懲罰因子，隨迭代次數(shù) k 逐漸增大
在模擬退火算法中，m(k) 隨外循環(huán)迭代次數(shù)逐漸增大，但在內(nèi)循環(huán)中應(yīng)保持不變。
5、模擬退火算法 Python 程序
懲罰函數(shù)法求解約束優(yōu)化線性規(guī)劃問(wèn)題
# 模擬退火算法 程序：懲罰函數(shù)法求解線性規(guī)劃問(wèn)題
# Program: SimulatedAnnealing_v2.py
# Purpose: Simulated annealing algorithm for function optimization
# v2.0: 使用懲罰函數(shù)法處理約束問(wèn)題
# Copyright 2021 YouCans, XUPT
# Crated：2021-05-01
# = 關(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)入模塊 YouCans, XUPT
import numpy as np # 導(dǎo)入模塊 numpy，并簡(jiǎn)寫(xiě)成 np
import matplotlib.pyplot as plt
from datetime import datetime
# 子程序：定義優(yōu)化問(wèn)題的目標(biāo)函數(shù)
def cal_Energy(X, nVar, mk): # m(k)：懲罰因子，隨迭代次數(shù) k 逐漸增大
p1 = (max(0, 6*X[0]+5*X[1]-60))**2
p2 = (max(0, 10*X[0]+20*X[1]-150))**2
fx = -(10*X[0]+9*X[1])
return fx+mk*(p1+p2)
# 子程序：模擬退火算法的參數(shù)設(shè)置
def ParameterSetting():
cName = "funcOpt" # 定義問(wèn)題名稱 YouCans, XUPT
nVar = 2 # 給定自變量數(shù)量，y=f(x1,..xn)
xMin = [0, 0] # 給定搜索空間的下限，x1_min,..xn_min
xMax = [8, 7.5] # 給定搜索空間的上限，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鏈長(zhǎng)度，也即內(nèi)循環(huán)運(yùn)行次數(shù)
scale = 0.5 # 定義搜索步長(zhǎng)，可以設(shè)為固定值或逐漸縮小
return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale
# 模擬退火算法
def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale):
# ====== 初始化隨機(jī)數(shù)發(fā)生器 ======
randseed = random.randint(1, 100)
random.seed(randseed) # 隨機(jī)數(shù)發(fā)生器設(shè)置種子，也可以設(shè)為指定整數(shù)
# ====== 隨機(jī)產(chǎn)生優(yōu)化問(wèn)題的初始解 ======
xInitial = np.zeros((nVar)) # 初始化，創(chuàng)建數(shù)組
for v in range(nVar):
# random.uniform(min,max) 在 [min,max] 范圍內(nèi)隨機(jī)生成一個(gè)實(shí)數(shù)
xInitial[v] = random.uniform(xMin[v], xMax[v])
# 調(diào)用子函數(shù) cal_Energy 計(jì)算當(dāng)前解的目標(biāo)函數(shù)值
fxInitial = cal_Energy(xInitial, nVar, 1) # m(k)：懲罰因子，初值為 1
# ====== 模擬退火算法初始化 ======
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 # 將初始解的目標(biāo)函數(shù)置為當(dāng)前值
fxBest = fxInitial # 將當(dāng)前解的目標(biāo)函數(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)前解的目標(biāo)函數(shù)值
recordFxBest = [] # 初始化，最佳解的目標(biāo)函數(shù)值
recordPBad = [] # 初始化，劣質(zhì)解的接受概率
kIter = 0 # 外循環(huán)迭代次數(shù)，溫度狀態(tài)數(shù)
totalMar = 0 # 總計(jì) Markov 鏈長(zhǎng)度
totalImprove = 0 # fxBest 改善次數(shù)
nMarkov = meanMarkov # 固定長(zhǎng)度 Markov鏈
# ====== 開(kāi)始模擬退火優(yōu)化 ======
# 外循環(huán)，直到當(dāng)前溫度達(dá)到終止溫度時(shí)結(jié)束
tNow = tInitial # 初始化當(dāng)前溫度(current temperature)
while tNow >= tFinal: # 外循環(huán)，直到當(dāng)前溫度達(dá)到終止溫度時(shí)結(jié)束
# 在當(dāng)前溫度下，進(jìn)行充分次數(shù)(nMarkov)的狀態(tài)轉(zhuǎn)移以達(dá)到熱平衡
kBetter = 0 # 獲得優(yōu)質(zhì)解的次數(shù)
kBadAccept = 0 # 接受劣質(zhì)解的次數(shù)
kBadRefuse = 0 # 拒絕劣質(zhì)解的次數(shù)
# ---內(nèi)循環(huán)，循環(huán)次數(shù)為Markov鏈長(zhǎng)度
for k in range(nMarkov): # 內(nèi)循環(huán)，循環(huán)次數(shù)為Markov鏈長(zhǎng)度
totalMar += 1 # 總 Markov鏈長(zhǎng)度計(jì)數(shù)器
# ---產(chǎn)生新解
# 產(chǎn)生新解：通過(guò)在當(dāng)前解附近隨機(jī)擾動(dòng)而產(chǎn)生新解，新解必須在 [min,max] 范圍內(nèi)
# 方案 1：只對(duì) n元變量中的一個(gè)進(jìn)行擾動(dòng)，其它 n-1個(gè)變量保持不變
xNew[:] = xNow[:]
v = random.randint(0, nVar-1) # 產(chǎn)生 [0,nVar-1]之間的隨機(jī)數(shù)
xNew[v] = xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1)
# random.normalvariate(0, 1)：產(chǎn)生服從均值為0、標(biāo)準(zhǔn)差為 1 的正態(tài)分布隨機(jī)實(shí)數(shù)
xNew[v] = max(min(xNew[v], xMax[v]), xMin[v]) # 保證新解在 [min,max] 范圍內(nèi)
# ---計(jì)算目標(biāo)函數(shù)和能量差
# 調(diào)用子函數(shù) cal_Energy 計(jì)算新解的目標(biāo)函數(shù)值
fxNew = cal_Energy(xNew, nVar, kIter)
deltaE = fxNew - fxNow
# ---按 Metropolis 準(zhǔn)則接受新解
# 接受判別：按照 Metropolis 準(zhǔn)則決定是否接受新解
if fxNew < fxNow: # 更優(yōu)解：如果新解的目標(biāo)函數(shù)好于當(dāng)前解，則接受新解
accept = True
kBetter += 1
else: # 容忍解：如果新解的目標(biāo)函數(shù)比當(dāng)前解差，則以一定概率接受新解
pAccept = math.exp(-deltaE / tNow) # 計(jì)算容忍解的狀態(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: # 如果新解的目標(biāo)函數(shù)好于最優(yōu)解，則將新解保存為最優(yōu)解
fxBest = fxNew
xBest[:] = xNew[:]
totalImprove += 1
scale = scale*0.99 # 可變搜索步長(zhǎng)，逐步減小搜索范圍，提高搜索精度
# ---內(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)前解的目標(biāo)函數(shù)值
recordFxBest.append(round(fxBest, 4)) # 最佳解的目標(biāo)函數(shù)值
recordPBad.append(round(pBadAccept, 4)) # 最佳解的目標(biāo)函數(shù)值
if kIter%10 == 0: # 模運(yùn)算，商的余數(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
fxBest = cal_Energy(xBest, nVar, kIter) # 由于迭代后懲罰因子增大，需隨之重構(gòu)增廣目標(biāo)函數(shù)
# ====== 結(jié)束模擬退火過(guò)程 ======
print('improve:{:d}'.format(totalImprove))
return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad
# 結(jié)果校驗(yàn)與輸出
def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter):
# ====== 優(yōu)化結(jié)果校驗(yàn)與輸出 ======
fxCheck = cal_Energy(xBest, nVar, kIter)
if abs(fxBest - fxCheck)>1e-3: # 檢驗(yàn)?zāi)繕?biāo)函數(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(cal_Energy(xBest,nVar,0)))
return
def main(): # YouCans, XUPT
# 參數(shù)設(shè)置，優(yōu)化問(wèn)題參數(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é)果校驗(yàn)與輸出
ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter)
if __name__ == '__main__':
main()
6、運(yùn)行結(jié)果
Optimization by simulated annealing algorithm:
x[0] = 6.577964
x[1] = 4.111469
f(x):-102.782857
參考文獻(xiàn)：
（1）胡山鷹，陳丙珍，非線性規(guī)劃問(wèn)題全局優(yōu)化的模擬退火法，清華大學(xué)學(xué)報(bào)，1997,37(6):5-9
（2）李歧強(qiáng)，具有約束指導(dǎo)的模擬退火算法，系統(tǒng)工程，2001,19(3):49-55
以上就是Python數(shù)學(xué)建模學(xué)習(xí)模擬退火算法約束條件處理示例解析的詳細(xì)內(nèi)容，更多關(guān)于數(shù)學(xué)建模學(xué)習(xí)模擬退火算法的資料請(qǐng)關(guān)注腳本之家其它相關(guān)文章！
相關(guān)文章
新一代爬蟲(chóng)利器Python?Playwright詳解
這篇文章主要為大家介紹了新一代爬蟲(chóng)利器Playwright，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下，希望能夠給你帶來(lái)幫助2021-12-12
Django實(shí)現(xiàn)web端tailf日志文件功能及實(shí)例詳解
這篇文章主要介紹了Django實(shí)現(xiàn)web端tailf日志文件功能，本文通過(guò)實(shí)例給大家介紹的非常詳細(xì)，具有一定的參考借鑒價(jià)值,需要的朋友可以參考下2019-07-07
Python讀取mat文件,并轉(zhuǎn)為csv文件的實(shí)例
今天小編就為大家分享一篇Python讀取mat文件,并轉(zhuǎn)為csv文件的實(shí)例，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2018-07-07
這篇文章主要介紹了Django框架下在URLconf中指定視圖緩存的方法,在Python豐富多彩的web框架中,Django是最具人氣的一個(gè),需要的朋友可以參考下2015-07-07
使用Python和百度語(yǔ)音識(shí)別生成視頻字幕的實(shí)現(xiàn)
這篇文章主要介紹了使用Python和百度語(yǔ)音識(shí)別生成視頻字幕，文中通過(guò)示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2020-04-04
在本篇文章里小編給大家整理的是一篇關(guān)于用Python繪制漫步圖實(shí)例講解內(nèi)容，需要的朋友們可以學(xué)習(xí)下。2020-02-02
python 讀文件,然后轉(zhuǎn)化為矩陣的實(shí)例
下面小編就為大家分享一篇python 讀文件,然后轉(zhuǎn)化為矩陣的實(shí)例，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2018-04-04
Python實(shí)現(xiàn)搜索算法的實(shí)例代碼
這篇文章主要介紹了Python實(shí)現(xiàn)搜索算法,文中通過(guò)實(shí)例代碼給大家分享高級(jí)搜索算法的想法，分為線性搜索和插值搜索，本文通過(guò)實(shí)例代碼給大家介紹的非常詳細(xì)，具有一定的參考借鑒價(jià)值，需要的朋友可以參考下2020-01-01
這篇文章主要介紹了python 為什么說(shuō)eval要慎用，文中通過(guò)示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2019-03-03
解決python寫(xiě)入帶有中文的字符到文件錯(cuò)誤的問(wèn)題
今天小編就為大家分享一篇解決python寫(xiě)入帶有中文的字符到文件錯(cuò)誤的問(wèn)題，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2019-01-01