Python数模笔记-模拟退火算法（3）整数规划问题
1、整數(shù)規(guī)劃問題
整數(shù)規(guī)劃問題在工業(yè)、經(jīng)濟(jì)、國(guó)防、醫(yī)療等各行各業(yè)應(yīng)用十分廣泛，是指規(guī)劃中的變量（全部或部分）限制為整數(shù)，屬于離散優(yōu)化問題（Discrete Optimization）。
線性規(guī)劃問題的最優(yōu)解可能是分?jǐn)?shù)或小數(shù)。但很多實(shí)際問題常常要求某些變量必須是整數(shù)解，例如：機(jī)器的臺(tái)數(shù)、工作的人數(shù)或裝貨的車數(shù)。根據(jù)對(duì)決策變量的不同要求，整數(shù)規(guī)劃又可以分為：純整數(shù)規(guī)劃、混合整數(shù)規(guī)劃、0-1整數(shù)規(guī)劃、混合0－1規(guī)劃。
整數(shù)規(guī)劃與線性規(guī)劃的差別只在于增加了整數(shù)約束。初看起來似乎只要把線性規(guī)劃得到的非整數(shù)解舍入化整就可以得到整數(shù)解，但是這樣化整后的整數(shù)解不一定是最優(yōu)解，甚至可能不是可行解。因此，通常需要采用特殊的方法來求解整數(shù)規(guī)劃，這比求解線性規(guī)劃問題復(fù)雜的多，以至于至今還沒有一般的多項(xiàng)式解法。因此，整數(shù)規(guī)劃問題被看作數(shù)學(xué)規(guī)劃中、甚至是數(shù)學(xué)中最困難的問題之一。
求解整數(shù)規(guī)劃比較成功又流行的方法是分支定界法和割平面法。核心思想是把整數(shù)規(guī)劃問題分解為一系列線性規(guī)劃問題，并追蹤整數(shù)規(guī)劃問題的上界（最優(yōu)可行解）和下界（最優(yōu)線性松弛解），逐步迭代收斂到最優(yōu)解。由于精確算法為指數(shù)復(fù)雜度，因此在有限時(shí)間內(nèi)也不能獲得全局最優(yōu)解，只能獲得近似最優(yōu)解。YouCans
目前整數(shù)規(guī)劃問題的優(yōu)化求解器主要有：IBM Cplex，Gurobi，FICO Xpress，SCIP，2018年中科院發(fā)布了CMIP混合整數(shù)規(guī)劃求解器。使用 Lingo 可以求解整數(shù)規(guī)劃問題，使用 Matlab 也可以用intlinprog 函數(shù)求解整數(shù)規(guī)劃問題，實(shí)際上都是使用軟件中內(nèi)建的求解器。Python 也可以使用第三方庫(kù)求解整數(shù)規(guī)劃問題，例如 Cvxpy、PuLp 都可以求解整數(shù)規(guī)劃問題，Cplex、Gurobi也有自己的python API。
歡迎關(guān)注 Youcans 原創(chuàng)系列，每周更新數(shù)模筆記
Python數(shù)模筆記-PuLP庫(kù)
Python數(shù)模筆記-StatsModels統(tǒng)計(jì)回歸
Python數(shù)模筆記-Sklearn
Python數(shù)模筆記-NetworkX
Python數(shù)模筆記-模擬退火算法
2、模擬退火算法處理整數(shù)約束
由于整數(shù)規(guī)劃問題在有限時(shí)間內(nèi)不能獲得全局最優(yōu)解，啟發(fā)式算法就有了用武之地。下面我們討論模擬退火算法處理整數(shù)約束，求解整數(shù)規(guī)劃問題。
上一篇文章中我們討論模擬退火算法處理線性規(guī)劃的約束條件時(shí)，方法比其它常用算法復(fù)雜的多。但是，模擬退火算法在處理整數(shù)約束時(shí)，方法卻極其簡(jiǎn)單：
對(duì)于決策變量為連續(xù)變量的一般優(yōu)化問題，基本的模擬退火算法在決策變量的取值范圍隨機(jī)產(chǎn)生初始解，新解則是在現(xiàn)有解的鄰域施加擾動(dòng)產(chǎn)生，算法上通過均勻分布或正態(tài)分布的隨機(jī)數(shù)來實(shí)現(xiàn)：
xInitial = random.uniform(xMin, xMax)
# random.uniform(min,max) 在 [min,max] 范圍內(nèi)隨機(jī)生成一個(gè)實(shí)數(shù)
xNew = xNow + scale * (xMax-xMin) * random.normalvariate(0, 1)
# random.normalvariate(0, 1)：產(chǎn)生服從均值為0、標(biāo)準(zhǔn)差為 1 的正態(tài)分布隨機(jī)實(shí)數(shù)
xNew = max(min(xNew, xMax), xMin) # 保證新解在 [min,max] 范圍內(nèi)
對(duì)于整數(shù)規(guī)劃問題，只要將產(chǎn)生初值/新解的隨機(jī)實(shí)數(shù)發(fā)生器 random.uniform、random.normalvariate 改為隨機(jī)整數(shù)發(fā)生器 random.randint即可：
xInitial = random.randint(xMin, xMax)
# random.randint(xMin, xMax) 產(chǎn)生 [min,max]之間的隨機(jī)整數(shù)
由于模擬退火算法與問題無關(guān)（Problem-independent），所以通常來說這樣處理并不會(huì)影響算法的性能：既不會(huì)引起不可行解，也不用擔(dān)心得不到最優(yōu)解——近似算法只能得到近似最優(yōu)解的，而且可以得到近似最優(yōu)解。
既然如此，更簡(jiǎn)單的處理方法，連隨機(jī)整數(shù)發(fā)生器都不需要，直接把線性規(guī)劃得到的非整數(shù)解舍入化整就可以了：
xNew = round(xNow + scale * (xMax-xMin) * random.normalvariate(0, 1))
# random.normalvariate(0, 1)：產(chǎn)生服從均值為0、標(biāo)準(zhǔn)差為 1 的正態(tài)分布隨機(jī)實(shí)數(shù)
xNew = max(min(xNew, xMax), xMin) # 保證新解在 [min,max] 范圍內(nèi)
這樣處理的好處是：（1）簡(jiǎn)單、直接，（2）便于實(shí)現(xiàn)所需的概率分布。
3、數(shù)模案例
為了便于理解，本文仍使用之前的案例。
3.1 問題描述：
某廠生產(chǎn)甲乙兩種飲料，每百箱甲飲料需用原料6千克、工人10名，獲利10萬(wàn)元；每百箱乙飲料需用原料5千克、工人20名，獲利9萬(wàn)元。
今工廠共有原料60千克、工人150名，又由于其他條件所限甲飲料產(chǎn)量不超過8百箱。
（5）若不允許散箱（按整百箱生產(chǎn)），如何安排生產(chǎn)計(jì)劃，即兩種飲料各生產(chǎn)多少使獲利最大？
3.2 問題分析：
問題（5）要求按整百箱生產(chǎn)，即要求決策變量為整數(shù)，是整數(shù)規(guī)劃問題。
對(duì)于模擬退火算法，基本算法中的初值/新解都是隨機(jī)生成的浮點(diǎn)實(shí)數(shù)（均勻分布或正態(tài)分布）。對(duì)于整數(shù)規(guī)劃問題，只要將產(chǎn)生初值/新解的隨機(jī)實(shí)數(shù)發(fā)生器改為隨機(jī)整數(shù)發(fā)生器即可，或者把線性規(guī)劃得到的非整數(shù)解舍入化整。
3.3 問題建模：
決策變量：
x1：甲飲料產(chǎn)量，正整數(shù)（單位：百箱）
x2：乙飲料產(chǎn)量，正整數(shù)（單位：百箱）
目標(biāo)函數(shù)：
max fx = 10*x1 + 9*x2
約束條件：
6*x1 + 5*x2 <= 60
10*x1 + 20*x2 <= 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
3.4 懲罰函數(shù)法求解約束優(yōu)化問題：
構(gòu)造懲罰函數(shù)：
p1 = (max(0, 6*x1+5*x2-60))**2
p2 = (max(0, 10*x1+20*x2-150))**2
說明：如存在等式約束，例如：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)保持不變。
4、模擬退火算法 Python 程序：求解整數(shù)規(guī)劃問題
# 模擬退火算法 程序：求解線性規(guī)劃問題（整數(shù)規(guī)劃） # Program: SimulatedAnnealing_v4.py # Purpose: Simulated annealing algorithm for function optimization # v4.0: 整數(shù)規(guī)劃：滿足決策變量的取值為整數(shù)（初值和新解都是隨機(jī)生成的整數(shù)） # 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)寫成 np import matplotlib.pyplot as plt from datetime import datetime# 子程序：定義優(yōu)化問題的目標(biāo)函數(shù) def cal_Energy(X, nVar, mk): # m(k)：懲罰因子，隨迭代次數(shù) k 逐漸增大p1 = (max(0, 6*X[0]+5*X[1]-60))**2p2 = (max(0, 10*X[0]+20*X[1]-150))**2fx = -(10*X[0]+9*X[1])return fx+mk*(p1+p2)# 子程序：模擬退火算法的參數(shù)設(shè)置 def ParameterSetting():cName = "funcOpt" # 定義問題名稱 YouCans, XUPTnVar = 2 # 給定自變量數(shù)量，y=f(x1,..xn)xMin = [0, 0] # 給定搜索空間的下限，x1_min,..xn_minxMax = [8, 8] # 給定搜索空間的上限，x1_max,..xn_maxtInitial = 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)化問題的初始解 ======xInitial = np.zeros((nVar)) # 初始化，創(chuàng)建數(shù)組for v in range(nVar):# xInitial[v] = random.uniform(xMin[v], xMax[v]) # 產(chǎn)生 [xMin, xMax] 范圍的隨機(jī)實(shí)數(shù)xInitial[v] = random.randint(xMin[v], xMax[v]) # 產(chǎn)生 [xMin, xMax] 范圍的隨機(jī)整數(shù)# 調(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鏈# ====== 開始模擬退火優(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)生新解：通過在當(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] = round(xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1))# 滿足決策變量為整數(shù)，采用最簡(jiǎn)單的方案：產(chǎn)生的新解按照四舍五入取整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 = TruekBetter += 1else: # 容忍解：如果新解的目標(biāo)函數(shù)比當(dāng)前解差，則以一定概率接受新解pAccept = math.exp(-deltaE / tNow) # 計(jì)算容忍解的狀態(tài)遷移概率if pAccept > random.random():accept = True # 接受劣質(zhì)解kBadAccept += 1else:accept = False # 拒絕劣質(zhì)解kBadRefuse += 1# 保存新解if accept == True: # 如果接受新解，則將新解保存為當(dāng)前解xNow[:] = xNew[:]fxNow = fxNewif fxNew < fxBest: # 如果新解的目標(biāo)函數(shù)好于最優(yōu)解，則將新解保存為最優(yōu)解fxBest = fxNewxBest[:] = xNew[:]totalImprove += 1scale = 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 * alfakIter = kIter + 1fxBest = cal_Energy(xBest, nVar, kIter) # 由于迭代后懲罰因子增大，需隨之重構(gòu)增廣目標(biāo)函數(shù)# ====== 結(jié)束模擬退火過程 ======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!")returnelse: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# = 關(guān)注 Youcans，分享原創(chuàng)系列 https://blog.csdn.net/youcans = # 主程序 def main(): # YouCans, XUPT# 參數(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é)果校驗(yàn)與輸出ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter)if __name__ == '__main__':main()5、運(yùn)行結(jié)果
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:7Optimization by simulated annealing algorithm:x[0] = 8.0x[1] = 2.0f(x) = -98.0參考文獻(xiàn)：
（1）田澎，楊自厚，張嗣瀛，一類非線性整數(shù)規(guī)劃的模擬退火求解，1993年控制理論及其應(yīng)用年會(huì)論文集，海洋出版社，1993，533-537.
版權(quán)說明：
= 關(guān)注 Youcans，分享原創(chuàng)系列 https://blog.csdn.net/youcans =
原創(chuàng)作品
Copyright 2021 YouCans, XUPT
Crated：2021-05-01
關(guān)注 Youcans，分享原創(chuàng)系列 https://blog.csdn.net/youcans
Python數(shù)模筆記-PuLP庫(kù)（1）線性規(guī)劃入門
Python數(shù)模筆記-PuLP庫(kù)（2）線性規(guī)劃進(jìn)階
Python數(shù)模筆記-PuLP庫(kù)（3）線性規(guī)劃實(shí)例
Python數(shù)模筆記-StatsModels 統(tǒng)計(jì)回歸（1）簡(jiǎn)介
Python數(shù)模筆記-StatsModels 統(tǒng)計(jì)回歸（2）線性回歸
Python數(shù)模筆記-StatsModels 統(tǒng)計(jì)回歸（3）模型數(shù)據(jù)的準(zhǔn)備
Python數(shù)模筆記-StatsModels 統(tǒng)計(jì)回歸（4）可視化
Python數(shù)模筆記-Sklearn （1）介紹
Python數(shù)模筆記-Sklearn （2）聚類分析
Python數(shù)模筆記-Sklearn （3）主成分分析
Python數(shù)模筆記-Sklearn （4）線性回歸
Python數(shù)模筆記-Sklearn （5）支持向量機(jī)
Python數(shù)模筆記-模擬退火算法（1）多變量函數(shù)優(yōu)化
Python數(shù)模筆記-模擬退火算法（2）約束條件的處理
Python數(shù)模筆記-模擬退火算法（3）整數(shù)規(guī)劃問題
Python數(shù)模筆記-模擬退火算法（4）旅行商問題
總結(jié)
以上是生活随笔為你收集整理的Python数模笔记-模拟退火算法（3）整数规划问题的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 【OpenCV 例程200篇】10. 图
- 下一篇： Python数模笔记-模拟退火算法（4）