帶約束條件的運籌規劃問題求解(模擬退火演算法實現)

碼頭牛牛發表於2023-04-18

0. 寫在前面

超級簡單的模擬退火演演算法實現ε٩(๑> ₃ <)۶з搭配最簡單的線性規劃模型進行講解!但是如果需要的話可以直接修改程式設計非線性問題哦(´つヮ⊂︎)

1. 模型描述及處理

1.1 線性規劃模型

\[max\,f(x)=10x_1+9x_2 \]

\(s.t.\)

\[6x_1+5x_2\leq{60}\tag{1} \]

\[10x_1+20x_2\leq{150}\tag{2} \]

\[0\leq{x_1}\leq{8}\tag{3} \]

\[0\leq{x_2}\leq{8}\tag{4} \]

1.2 引入懲罰函式處理模型

對約束條件引入懲罰函式:

  • 對約束條件(1),懲罰函式為:\(p_1=max(0,6x_1+5x_2-60)^2\)

  • 對約束條件(2),懲罰函式為:\(p_2=max(0,10x_1+20x_2-150)^2\)

那麼,該問題的懲罰函式可以表示為:

\[P(x)=p_1+p_2 \]

由此,可將該問題的約束條件放入目標函式中,此時模型變為:

\[min\,g(x)=-(10x_1+9x_2)+P(x)\quad\forall{x_1,x_2}\in{[0,8]} \]

2. 程式實現

# 模擬退火演演算法 程式:求解線性規劃問題(整數規劃)
# 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:0.000000,4.000000,	f(x_Initial):-36.000000
i:0,t(i):100.00, badAccept:0.925373, f(x)_best:-152.000000
i:10,t(i):81.71, badAccept:0.671053, f(x)_best:-98.000000
i:20,t(i):66.76, badAccept:0.722892, f(x)_best:-98.000000
i:30,t(i):54.55, badAccept:0.704225, f(x)_best:-98.000000
i:40,t(i):44.57, badAccept:0.542169, f(x)_best:-98.000000
i:50,t(i):36.42, badAccept:0.435294, f(x)_best:-98.000000
i:60,t(i):29.76, badAccept:0.359551, f(x)_best:-98.000000
i:70,t(i):24.31, badAccept:0.717647, f(x)_best:-98.000000
i:80,t(i):19.86, badAccept:0.388235, f(x)_best:-98.000000
i:90,t(i):16.23, badAccept:0.555556, f(x)_best:-98.000000
i:100,t(i):13.26, badAccept:0.482353, f(x)_best:-98.000000
i:110,t(i):10.84, badAccept:0.527473, f(x)_best:-98.000000
i:120,t(i):8.85, badAccept:0.164948, f(x)_best:-98.000000
i:130,t(i):7.23, badAccept:0.305263, f(x)_best:-98.000000
i:140,t(i):5.91, badAccept:0.120000, f(x)_best:-98.000000
i:150,t(i):4.83, badAccept:0.422680, f(x)_best:-98.000000
i:160,t(i):3.95, badAccept:0.111111, f(x)_best:-98.000000
i:170,t(i):3.22, badAccept:0.350000, f(x)_best:-98.000000
i:180,t(i):2.63, badAccept:0.280000, f(x)_best:-98.000000
i:190,t(i):2.15, badAccept:0.310000, f(x)_best:-98.000000
i:200,t(i):1.76, badAccept:0.390000, f(x)_best:-98.000000
i:210,t(i):1.44, badAccept:0.390000, f(x)_best:-98.000000
i:220,t(i):1.17, badAccept:0.380000, f(x)_best:-98.000000
improve:10

Optimization by simulated annealing algorithm:
	x[0] = 8.0
	x[1] = 2.0

	f(x) = -98.0

相關文章