溫馨提示×
您好，登錄后才能下訂單哦！
點擊 登錄注冊 即表示同意《億速云用戶服務條款》
您好，登錄后才能下訂單哦！
本篇內容主要講解“Python數學建模學習之如何模擬退火算法整數規劃”，感興趣的朋友不妨來看看。本文介紹的方法操作簡單快捷，實用性強。下面就讓小編來帶大家學習“Python數學建模學習之如何模擬退火算法整數規劃”吧!
線性規劃問題的最優解可能是分數或小數。但很多實際問題常常要求某些變量必須是整數解，例如：機器的臺數、工作的人數或裝貨的車數。根據對決策變量的不同要求，整數規劃又可以分為：純整數規劃、混合整數規劃、0-1整數規劃、混合0－1規劃。
整數規劃與線性規劃的差別只在于增加了整數約束。初看起來似乎只要把線性規劃得到的非整數解舍入化整就可以得到整數解，但是這樣化整后的整數解不一定是最優解，甚至可能不是可行解。因此，通常需要采用特殊的方法來求解整數規劃，這比求解線性規劃問題復雜的多，以至于至今還沒有一般的多項式解法。因此，整數規劃問題被看作數學規劃中、甚至是數學中最困難的問題之一。
求解整數規劃比較成功又流行的方法是分支定界法和割平面法。核心思想是把整數規劃問題分解為一系列線性規劃問題，并追蹤整數規劃問題的上界（最優可行解）和下界（最優線性松弛解），逐步迭代收斂到最優解。由于精確算法為指數復雜度，因此在有限時間內也不能獲得全局最優解，只能獲得近似最優解。YouCans
目前整數規劃問題的優化求解器主要有：IBM Cplex，Gurobi，FICO Xpress，SCIP，2018年中科院發布了CMIP混合整數規劃求解器。使用 Lingo 可以求解整數規劃問題，使用 Matlab 也可以用intlinprog 函數求解整數規劃問題，實際上都是使用軟件中內建的求解器。Python 也可以使用第三方庫求解整數規劃問題，例如 Cvxpy、PuLp 都可以求解整數規劃問題，Cplex、Gurobi也有自己的python API。
由于整數規劃問題在有限時間內不能獲得全局最優解，啟發式算法就有了用武之地。下面我們討論模擬退火算法處理整數約束，求解整數規劃問題。
上一篇文章中我們討論模擬退火算法處理線性規劃的約束條件時，方法比其它常用算法復雜的多。但是，模擬退火算法在處理整數約束時，方法卻極其簡單：
對于決策變量為連續變量的一般優化問題，基本的模擬退火算法在決策變量的取值范圍隨機產生初始解，新解則是在現有解的鄰域施加擾動產生，算法上通過均勻分布或正態分布的隨機數來實現：
xInitial = random.uniform(xMin, xMax)
# random.uniform(min,max) 在 [min,max] 范圍內隨機生成一個實數xNew = xNow + scale * (xMax-xMin) * random.normalvariate(0, 1)
# random.normalvariate(0, 1)：產生服從均值為0、標準差為 1 的正態分布隨機實數
xNew = max(min(xNew, xMax), xMin) # 保證新解在 [min,max] 范圍內
對于整數規劃問題，只要將產生初值/新解的隨機實數發生器 random.uniform、random.normalvariate 改為隨機整數發生器 random.randint即可：
xInitial = random.randint(xMin, xMax)
# random.randint(xMin, xMax) 產生 [min,max]之間的隨機整數
由于模擬退火算法與問題無關（Problem-independent），所以通常來說這樣處理并不會影響算法的性能：既不會引起不可行解，也不用擔心得不到最優解——近似算法只能得到近似最優解的，而且可以得到近似最優解。
既然如此，更簡單的處理方法，連隨機整數發生器都不需要，直接把線性規劃得到的非整數解舍入化整就可以了：
xNew = round(xNow + scale * (xMax-xMin) * random.normalvariate(0, 1))
# random.normalvariate(0, 1)：產生服從均值為0、標準差為 1 的正態分布隨機實數
xNew = max(min(xNew, xMax), xMin) # 保證新解在 [min,max] 范圍內
這樣處理的好處是：（1）簡單、直接，（2）便于實現所需的概率分布。
為了便于理解，本文仍使用之前的案例。
某廠生產甲乙兩種飲料，每百箱甲飲料需用原料6千克、工人10名，獲利10萬元；每百箱乙飲料需用原料5千克、工人20名，獲利9萬元。
今工廠共有原料60千克、工人150名，又由于其他條件所限甲飲料產量不超過8百箱。
（5）若不允許散箱（按整百箱生產），如何安排生產計劃，即兩種飲料各生產多少使獲利最大？
問題（5）要求按整百箱生產，即要求決策變量為整數，是整數規劃問題。
對于模擬退火算法，基本算法中的初值/新解都是隨機生成的浮點實數（均勻分布或正態分布）。對于整數規劃問題，只要將產生初值/新解的隨機實數發生器改為隨機整數發生器即可，或者把線性規劃得到的非整數解舍入化整。
決策變量：
x1：甲飲料產量，正整數（單位：百箱）
x2：乙飲料產量，正整數（單位：百箱）
目標函數：
max fx = 10*x1 + 9*x2
約束條件：
6*x1 + 5*x2 <= 60
10*x1 + 20*x2 <= 150
取值范圍：
給定條件：x1, x2 >= 0，x1 <= 8
推導條件：由 x1,x2>=0 和 10*x1+20*x2<=150 可知：0<=x1<=15；0<=x2<=7.5
因此，0 <= x1<=8，0 <= x2<=7.5
構造懲罰函數：
p1 = (max(0, 6*x1+5*x2-60))**2
p2 = (max(0, 10*x1+20*x2-150))**2
說明：如存在等式約束，例如：x1 + 2*x2 = m，也可以轉化為懲罰函數：
p3 = (x1+2*x2-m)**2
P(x) = p1 + p2 + …
構造增廣目標函數：
L(x,m(k)) = min(fx) + m(k)*P(x)
m(k)：懲罰因子，隨迭代次數 k 逐漸增大
在模擬退火算法中，m(k) 隨外循環迭代次數逐漸增大，但在內循環中應保持不變。
# 模擬退火算法 程序：求解線性規劃問題（整數規劃）
# Program: SimulatedAnnealing_v4.py
# Purpose: Simulated annealing algorithm for function optimization
# v4.0: 整數規劃：滿足決策變量的取值為整數（初值和新解都是隨機生成的整數）
# Copyright 2021 YouCans, XUPT
# Crated：2021-05-01
# = 關注 Youcans，分享原創系列 https://blog.csdn.net/youcans =
# -*- coding: utf-8 -*-
import math # 導入模塊
import random # 導入模塊
import pandas as pd # 導入模塊 YouCans, XUPT
import numpy as np # 導入模塊 numpy，并簡寫成 np
import matplotlib.pyplot as plt
from datetime import datetime
# 子程序：定義優化問題的目標函數
def cal_Energy(X, nVar, mk): # m(k)：懲罰因子，隨迭代次數 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)
# 子程序：模擬退火算法的參數設置
def ParameterSetting():
cName = "funcOpt" # 定義問題名稱 YouCans, XUPT
nVar = 2 # 給定自變量數量，y=f(x1,..xn)
xMin = [0, 0] # 給定搜索空間的下限，x1_min,..xn_min
xMax = [8, 8] # 給定搜索空間的上限，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):
# xInitial[v] = random.uniform(xMin[v], xMax[v]) # 產生 [xMin, xMax] 范圍的隨機實數
xInitial[v] = random.randint(xMin[v], xMax[v]) # 產生 [xMin, xMax] 范圍的隨機整數
# 調用子函數 cal_Energy 計算當前解的目標函數值
fxInitial = cal_Energy(xInitial, nVar, 1) # m(k)：懲罰因子，初值為 1
# ====== 模擬退火算法初始化 ======
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] = round(xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1))
# 滿足決策變量為整數，采用最簡單的方案：產生的新解按照四舍五入取整
xNew[v] = max(min(xNew[v], xMax[v]), xMin[v]) # 保證新解在 [min,max] 范圍內
# ---計算目標函數和能量差
# 調用子函數 cal_Energy 計算新解的目標函數值
fxNew = cal_Energy(xNew, nVar, kIter)
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
fxBest = cal_Energy(xBest, nVar, kIter) # 由于迭代后懲罰因子增大，需隨之重構增廣目標函數
# ====== 結束模擬退火過程 ======
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, kIter)
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[{}] = {:.1f}'.format(i,xBest[i]))
print('\n\tf(x) = {:.1f}'.format(cal_Energy(xBest,nVar,0)))
return
# 主程序
def main(): # YouCans, XUPT
# 參數設置，優化問題參數定義，模擬退火算法參數設置
[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()x_Initial:2.000000,7.000000, f(x_Initial):17.000000 i:0,t(i):100.00, badAccept:0.814286, f(x)_best:-152.000000 i:10,t(i):81.71, badAccept:0.635135, f(x)_best:-98.000000 i:20,t(i):66.76, badAccept:0.782051, f(x)_best:-98.000000 ... i:200,t(i):1.76, badAccept:0.090000, f(x)_best:-98.000000 i:210,t(i):1.44, badAccept:0.120000, f(x)_best:-98.000000 i:220,t(i):1.17, badAccept:0.130000, f(x)_best:-98.000000 improve:7 Optimization by simulated annealing algorithm: x[0] = 8.0 x[1] = 2.0 f(x) = -98.0
到此，相信大家對“Python數學建模學習之如何模擬退火算法整數規劃”有了更深的了解，不妨來實際操作一番吧！這里是億速云網站，更多相關內容可以進入相關頻道進行查詢，關注我們，繼續學習！
免責聲明：本站發布的內容（圖片、視頻和文字）以原創、轉載和分享為主，文章觀點不代表本網站立場，如果涉及侵權請聯系站長郵箱：is@yisu.com進行舉報，并提供相關證據，一經查實，將立刻刪除涉嫌侵權內容。
