Python數(shù)學(xué)建模學(xué)習(xí)之模擬退火算法約束條件處理的實現(xiàn)方法

這篇文章主要講解了“Python數(shù)學(xué)建模學(xué)習(xí)之模擬退火算法約束條件處理的實現(xiàn)方法”，文中的講解內(nèi)容簡單清晰，易于學(xué)習(xí)與理解，下面請大家跟著小編的思路慢慢深入，一起來研究和學(xué)習(xí)“Python數(shù)學(xué)建模學(xué)習(xí)之模擬退火算法約束條件處理的實現(xiàn)方法”吧！

1、最優(yōu)化與線性規(guī)劃

最優(yōu)化問題的三要素是決策變量、目標函數(shù)和約束條件。

線性規(guī)劃（Linear programming），是研究線性約束條件下線性目標函數(shù)的極值問題的優(yōu)化方法，常用于解決利用現(xiàn)有的資源得到最優(yōu)決策的問題。

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

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

本文使用懲罰函數(shù)法，分析模擬退火算法處理線性規(guī)劃問題，相關(guān)內(nèi)容也適用于非線性規(guī)劃問題。

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

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

1.決策變量取值的上下限約束。
此類約束條件比較容易處理，只要設(shè)定初始解、新解在決策變量取值的上下限之間就可以解決。例如：

（1）設(shè)置產(chǎn)生新解的隨機數(shù)的上下限為決策變量的上下限，即 [Xmin, Xmax]；

（2）設(shè)置產(chǎn)生新解的隨機數(shù)的上下限為當前解與決策變量的上下限，即 [Xnow, Xmax]；

（3）通過條件判斷，當新解超出決策變量上下限，則令其取上下限，即 xNew = max(min(xNew, xMax), xMin)。當然，這些處理方式，都會影響隨機數(shù)的概率分布，因而也影響模擬退火算法的優(yōu)化性能，在此不做深入討論。

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

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

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

3、懲罰函數(shù)法

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

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

4、數(shù)模案例

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

4.1 問題描述：

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

4.2 問題建模：

決策變量：
　　　　x1：甲飲料產(chǎn)量（單位：百箱）
　　　　x2：乙飲料產(chǎn)量（單位：百箱）
　　目標函數(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)化問題：

構(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)造增廣目標函數(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ī)劃問題

# 模擬退火算法 程序：懲罰函數(shù)法求解線性規(guī)劃問題
# Program: SimulatedAnnealing_v2.py
# Purpose: Simulated annealing algorithm for function optimization
# v2.0: 使用懲罰函數(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，并簡寫成 np
import matplotlib.pyplot as plt     
from datetime import datetime
# 子程序：定義優(yōu)化問題的目標函數(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"           # 定義問題名稱 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鏈長度，也即內(nèi)循環(huán)運行次數(shù)
    scale   = 0.5               # 定義搜索步長，可以設(shè)為固定值或逐漸縮小
    return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale
# 模擬退火算法
def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale):
    # ====== 初始化隨機數(shù)發(fā)生器 ======
    randseed = random.randint(1, 100)
    random.seed(randseed)  # 隨機數(shù)發(fā)生器設(shè)置種子，也可以設(shè)為指定整數(shù)
    # ====== 隨機產(chǎn)生優(yōu)化問題的初始解 ======
    xInitial = np.zeros((nVar))   # 初始化，創(chuàng)建數(shù)組
    for v in range(nVar):
        # random.uniform(min,max) 在 [min,max] 范圍內(nèi)隨機生成一個實數(shù)
        xInitial[v] = random.uniform(xMin[v], xMax[v])
    # 調(diào)用子函數(shù) cal_Energy 計算當前解的目標函數(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[:]          # 初始化當前解，將初始解置為當前解
    xBest[:] = xInitial[:]          # 初始化最優(yōu)解，將當前解置為最優(yōu)解
    fxNow  = fxInitial              # 將初始解的目標函數(shù)置為當前值
    fxBest = fxInitial              # 將當前解的目標函數(shù)置為最優(yōu)值
    print('x_Initial:{:.6f},{:.6f},\tf(x_Initial):{:.6f}'.format(xInitial[0], xInitial[1], fxInitial))
    recordIter = []                 # 初始化，外循環(huán)次數(shù)
    recordFxNow = []                # 初始化，當前解的目標函數(shù)值
    recordFxBest = []               # 初始化，最佳解的目標函數(shù)值
    recordPBad = []                 # 初始化，劣質(zhì)解的接受概率
    kIter = 0                       # 外循環(huán)迭代次數(shù)，溫度狀態(tài)數(shù)
    totalMar = 0                    # 總計 Markov 鏈長度
    totalImprove = 0                # fxBest 改善次數(shù)
    nMarkov = meanMarkov            # 固定長度 Markov鏈
    # ====== 開始模擬退火優(yōu)化 ======
    # 外循環(huán)，直到當前溫度達到終止溫度時結(jié)束
    tNow = tInitial                 # 初始化當前溫度(current temperature)
    while tNow >= tFinal:           # 外循環(huán)，直到當前溫度達到終止溫度時結(jié)束
        # 在當前溫度下，進行充分次數(shù)(nMarkov)的狀態(tài)轉(zhuǎn)移以達到熱平衡
        kBetter = 0                 # 獲得優(yōu)質(zhì)解的次數(shù)
        kBadAccept = 0              # 接受劣質(zhì)解的次數(shù)
        kBadRefuse = 0              # 拒絕劣質(zhì)解的次數(shù)
        # ---內(nèi)循環(huán)，循環(huán)次數(shù)為Markov鏈長度
        for k in range(nMarkov):    # 內(nèi)循環(huán)，循環(huán)次數(shù)為Markov鏈長度
            totalMar += 1           # 總 Markov鏈長度計數(shù)器
            # ---產(chǎn)生新解
            # 產(chǎn)生新解：通過在當前解附近隨機擾動而產(chǎn)生新解，新解必須在 [min,max] 范圍內(nèi)
            # 方案 1：只對 n元變量中的一個進行擾動，其它 n-1個變量保持不變
            xNew[:] = xNow[:]
            v = random.randint(0, nVar-1)   # 產(chǎn)生 [0,nVar-1]之間的隨機數(shù)
            xNew[v] = xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1)
            # random.normalvariate(0, 1)：產(chǎn)生服從均值為0、標準差為 1 的正態(tài)分布隨機實數(shù)
            xNew[v] = max(min(xNew[v], xMax[v]), xMin[v])  # 保證新解在 [min,max] 范圍內(nèi)
            # ---計算目標函數(shù)和能量差
            # 調(diào)用子函數(shù) cal_Energy 計算新解的目標函數(shù)值
            fxNew = cal_Energy(xNew, nVar, kIter)
            deltaE = fxNew - fxNow
            # ---按 Metropolis 準則接受新解
            # 接受判別：按照 Metropolis 準則決定是否接受新解
            if fxNew < fxNow:  # 更優(yōu)解：如果新解的目標函數(shù)好于當前解，則接受新解
                accept = True
                kBetter += 1
            else:  # 容忍解：如果新解的目標函數(shù)比當前解差，則以一定概率接受新解
                pAccept = math.exp(-deltaE / tNow)  # 計算容忍解的狀態(tài)遷移概率
                if pAccept > random.random():
                    accept = True  # 接受劣質(zhì)解
                    kBadAccept += 1
                else:
                    accept = False  # 拒絕劣質(zhì)解
                    kBadRefuse += 1
            # 保存新解
            if accept == True:  # 如果接受新解，則將新解保存為當前解
                xNow[:] = xNew[:]
                fxNow = fxNew
                if fxNew < fxBest:  # 如果新解的目標函數(shù)好于最優(yōu)解，則將新解保存為最優(yōu)解
                    fxBest = fxNew
                    xBest[:] = xNew[:]
                    totalImprove += 1
                    scale = scale*0.99  # 可變搜索步長，逐步減小搜索范圍，提高搜索精度                    
        # ---內(nèi)循環(huán)結(jié)束后的數(shù)據(jù)整理
        # 完成當前溫度的搜索，保存數(shù)據(jù)和輸出
        pBadAccept = kBadAccept / (kBadAccept + kBadRefuse)  # 劣質(zhì)解的接受概率
        recordIter.append(kIter)  # 當前外循環(huán)次數(shù)
        recordFxNow.append(round(fxNow, 4))  # 當前解的目標函數(shù)值
        recordFxBest.append(round(fxBest, 4))  # 最佳解的目標函數(shù)值
        recordPBad.append(round(pBadAccept, 4))  # 最佳解的目標函數(shù)值
        if kIter%10 == 0:                           # 模運算，商的余數(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)增廣目標函數(shù)
        # ====== 結(jié)束模擬退火過程 ======
    print('improve:{:d}'.format(totalImprove))
    return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad
# 結(jié)果校驗與輸出
def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter):
    # ====== 優(yōu)化結(jié)果校驗與輸出 ======
    fxCheck = cal_Energy(xBest, nVar, kIter)
    if abs(fxBest - fxCheck)>1e-3:   # 檢驗?zāi)繕撕瘮?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)化問題參數(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é)果校驗與輸出
    ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter)
if __name__ == '__main__':
    main()

6、運行結(jié)果

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

感謝各位的閱讀，以上就是“Python數(shù)學(xué)建模學(xué)習(xí)之模擬退火算法約束條件處理的實現(xiàn)方法”的內(nèi)容了，經(jīng)過本文的學(xué)習(xí)后，相信大家對Python數(shù)學(xué)建模學(xué)習(xí)之模擬退火算法約束條件處理的實現(xiàn)方法這一問題有了更深刻的體會，具體使用情況還需要大家實踐驗證。這里是億速云，小編將為大家推送更多相關(guān)知識點的文章，歡迎關(guān)注！