---

1、最優化與線性規劃

最優化問題的三要素是決策變量、目標函數和約束條件。

線性規劃（Linear programming），是研究線性約束條件下線性目標函數的極值問題的優化方法，常用于解決利用現有的資源得到最優決策的問題。

簡單的線性規劃問題可以用 Lingo軟件求解，Matlab、Python 中也有求解線性規劃問題的庫函數或求解器，很容易學習和使用，并不需要用模擬退火算法。但是，由一般線性規劃問題所衍生的整數規劃、混合規劃、0/1規劃、二次規劃、非線性規劃、組合優化問題，則并不是調用某個庫函數都能處理的。而模擬退火算法在很多復雜問題中具有較好的適應性，可以作為一種入門的通用智能算法來學習。

也就是說，如果只是處理線性規劃問題，就不要用模擬退火算法了。但如果是現有方法無法處理的復雜優化問題，或者對某類、某個優化問題你不知道用什么方法處理了，這時用模擬退火算法還是能解決的。

本文使用懲罰函數法，分析模擬退火算法處理線性規劃問題，相關內容也適用于非線性規劃問題。

歡迎關注 Youcans 原創系列，每周更新數模筆記

Python數模筆記-PuLP庫
Python數模筆記-StatsModels統計回歸
Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火算法

2、模擬退火算法處理約束條件

線性規劃問題是約束優化問題，而模擬退火算法則更適合處理無約束優化問題。對于優化問題中的約束條件，模擬退火算法有幾種常用的處理方法：

• 1.決策變量取值的上下限約束。
  此類約束條件比較容易處理，只要設定初始解、新解在決策變量取值的上下限之間就可以解決。例如，（1）設置產生新解的隨機數的上下限為決策變量的上下限，即 [Xmin, Xmax]；（2）設置產生新解的隨機數的上下限為當前解與決策變量的上下限，即 [Xnow, Xmax]；（3）通過條件判斷，當新解超出決策變量上下限，則令其取上下限，即 xNew = max(min(xNew, xMax), xMin)。當然，這些處理方式，都會影響隨機數的概率分布，因而也影響模擬退火算法的優化性能，在此不做深入討論。

• 2.檢驗法處理不等式約束問題。
  在模擬退火算法的迭代過程中，將每次產生的新解代入每個不等式約束函數，判斷是否滿足約束條件；如果新解不滿足約束條件，則舍棄這個新解，返回重新產生一個新解進行檢驗，直到產生的新解滿足全部約束條件為止。這個方法的思路簡單，每次迭代都在可行域內進行，但是對于約束條件眾多、苛刻的復雜問題，多次產生的新解都不能滿足約束條件，會使計算時間很長，甚至停滯不前。

• 3.消元法處理等式約束問題。
  對于等式約束，很難通過隨機產生的新解滿足約束條件，通常不能使用檢驗法處理。消元法是通過解方程將等式約束中的某個決策變量表示為其它決策變量的線性關系后，代入目標函數和不等式約束條件中，從而消去該約束條件。消元法不僅解決了等式約束問題，而且減少了決策變量的數量，從而有效簡化了優化問題的復雜度，是一舉兩得的處理方法。但是，對于非線性等式約束，求解非線性方程組也是非常困難的，消元法并不是普遍都能適用的。

• 4.更為通用的處理約束條件的方法是懲罰函數法，以下進行介紹。YouCans, XUPT

3、懲罰函數法

懲罰函數法是一類常用的處理約束條件的技術，在模擬退火算法中處理約束條件非常有效。方法的思想是將約束條件轉化為懲罰函數，附加在原有的目標函數上構造新的目標函數；當不滿足約束條件時，通過懲罰函數使新的目標函數變差而被舍棄。

懲罰函數法有外點法和內點法。外點法對可行域外的點（即不滿足約束的點）施加懲罰，對可行域內部的點不懲罰，從而使迭代點向可行域D逼近。 內點法是在可行域內部進行搜索，約束邊界起到類似圍墻的作用，使目標函數無法穿過，就把搜索點限制在可行域內了，因此只適用于不等式約束。

考慮約束優化問題：　　
$$min:f(X)$$
滿足限制
$$c_i(X)\le 0$$
懲罰函數法將問題轉化成如下無約束問題
$$min:L_k(X)=f(X)+{\sigma }_k\sum _{i}g(c_i(X))$$
其中
$$g(c_i(X))=max(0,c_i(X){)}^2$$

在上述方程，$g(ci(X))$ 稱為外部罰函數，$sigma(k)$ 稱為懲罰因子。在每一次迭代中，$sigma(k)$ 都增大，然后求解該無約束問題。將每一次迭代的結果將組成一個序列，此序列的極限即為原約束問題的解。

4、數模案例

雖然對于線性規劃問題并不推薦使用模擬退火算法求解。但為了便于理解，本文仍使用之前的線性規劃問題作為處理約束條件的案例。對于非線性規劃問題，以及非線性約束問題，處理方法都是類似的，將在后續進行介紹。

4.1 問題描述：

某廠生產甲乙兩種飲料，每百箱甲飲料需用原料6千克、工人10名，獲利10萬元；每百箱乙飲料需用原料5千克、工人20名，獲利9萬元。
今工廠共有原料60千克、工人150名，又由于其他條件所限甲飲料產量不超過8百箱。
　　（1）問如何安排生產計劃，即兩種飲料各生產多少使獲利最大？

4.2 問題建模：

決策變量：
　　　　x1：甲飲料產量（單位：百箱）
　　　　x2：乙飲料產量（單位：百箱）
　　目標函數：
　　　　max fx = 10x1 + 9x2
　　約束條件：
　　　　6x1 + 5x2 <= 60
　　　　10x1 + 20x2 <= 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

4.3 懲罰函數法求解約束優化問題：

構造懲罰函數：
　　　　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) 隨外循環迭代次數逐漸增大，但在內循環中應保持不變。

5、模擬退火算法 Python 程序：懲罰函數法求解約束優化線性規劃問題

# 模擬退火算法 程序：懲罰函數法求解線性規劃問題 # Program: SimulatedAnnealing_v2.py # Purpose: Simulated annealing algorithm for function optimization # v2.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))**2p2 = (max(0, 10*X[0]+20*X[1]-150))**2fx = -(10*X[0]+9*X[1])return fx+mk*(p1+p2)# 子程序：模擬退火算法的參數設置 def ParameterSetting():cName = "funcOpt" # 定義問題名稱 YouCans, XUPTnVar = 2 # 給定自變量數量，y=f(x1,..xn)xMin = [0, 0] # 給定搜索空間的下限，x1_min,..xn_minxMax = [8, 7.5] # 給定搜索空間的上限，x1_max,..xn_maxtInitial = 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):# random.uniform(min,max) 在 [min,max] 范圍內隨機生成一個實數xInitial[v] = random.uniform(xMin[v], xMax[v])# 調用子函數 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] = xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1)# random.normalvariate(0, 1)：產生服從均值為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 = TruekBetter += 1else: # 容忍解：如果新解的目標函數比當前解差，則以一定概率接受新解pAccept = math.exp(-deltaE / tNow) # 計算容忍解的狀態遷移概率if pAccept > random.random():accept = True # 接受劣質解kBadAccept += 1else:accept = False # 拒絕劣質解kBadRefuse += 1# 保存新解if accept == True: # 如果接受新解，則將新解保存為當前解xNow[:] = xNew[:]fxNow = fxNewif fxNew < fxBest: # 如果新解的目標函數好于最優解，則將新解保存為最優解fxBest = fxNewxBest[:] = xNew[:]totalImprove += 1scale = 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 * alfakIter = kIter + 1fxBest = 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!")returnelse: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# 主程序= 關注 Youcans，分享原創系列 https://blog.csdn.net/youcans = 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)# = 關注 Youcans，分享原創系列 https://blog.csdn.net/youcans = if __name__ == '__main__':main()

6、運行結果

Optimization by simulated annealing algorithm:x[0] = 6.577964x[1] = 4.111469f(x):-102.782857

參考文獻：

（1）胡山鷹，陳丙珍，非線性規劃問題全局優化的模擬退火法，清華大學學報，1997,37(6):5-9
（2）李歧強，具有約束指導的模擬退火算法，系統工程，2001,19(3):49-55

版權說明：

原創作品 = 關注 Youcans，分享原創系列 https://blog.csdn.net/youcans =
Copyright 2021 YouCans, XUPT
Crated：2021-05-01

關注 Youcans，分享原創系列 https://blog.csdn.net/youcans

Python數模筆記-PuLP庫（1）線性規劃入門
Python數模筆記-PuLP庫（2）線性規劃進階
Python數模筆記-PuLP庫（3）線性規劃實例
Python數模筆記-StatsModels 統計回歸（1）簡介
Python數模筆記-StatsModels 統計回歸（2）線性回歸
Python數模筆記-StatsModels 統計回歸（3）模型數據的準備
Python數模筆記-StatsModels 統計回歸（4）可視化
Python數模筆記-Sklearn （1）介紹
Python數模筆記-Sklearn （2）聚類分析
Python數模筆記-Sklearn （3）主成分分析
Python數模筆記-Sklearn （4）線性回歸
Python數模筆記-Sklearn （5）支持向量機
Python數模筆記-模擬退火算法（1）多變量函數優化
Python數模筆記-模擬退火算法（2）約束條件的處理
Python數模筆記-模擬退火算法（3）整數規劃問題
Python數模筆記-模擬退火算法（4）旅行商問題

總結

以上是生活随笔為你收集整理的Python数模笔记-模拟退火算法（2）约束条件的处理的全部內容，希望文章能夠幫你解決所遇到的問題。

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