(資料圖片)

超級(jí)簡(jiǎn)單的模擬退火算法實(shí)現(xiàn)ε?(?> ? <)?з搭配最簡(jiǎn)單的線性規(guī)劃模型進(jìn)行講解！但是如果需要的話可以直接修改編程非線性問(wèn)題哦(′つヮ??)

\(s.t.\)

對(duì)約束條件引入懲罰函數(shù)：

對(duì)約束條件(1)，懲罰函數(shù)為：\(p_1=max(0,6x_1+5x_2-60)^2\)

對(duì)約束條件(2)，懲罰函數(shù)為：\(p_2=max(0,10x_1+20x_2-150)^2\)

那么，該問(wèn)題的懲罰函數(shù)可以表示為：

由此，可將該問(wèn)題的約束條件放入目標(biāo)函數(shù)中，此時(shí)模型變?yōu)椋?/p>\[min\,g(x)=-(10x_1+9x_2)+P(x)\quad\forall{x_1,x_2}\in{[0,8]}\]2. 程序?qū)崿F(xiàn)

# 模擬退火算法 程序：求解線性規(guī)劃問(wèn)題（整數(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, XUPTimport numpy as np                  # 導(dǎo)入模塊 numpy，并簡(jiǎn)寫(xiě)成 npimport matplotlib.pyplot as pltfrom datetime import datetime # 子程序：定義優(yōu)化問(wèn)題的目標(biāo)函數(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"           # 定義問(wèn)題名稱(chēng) YouCans, XUPT    nVar = 2                    # 給定自變量數(shù)量，y=f(x1,..xn)    xMin = [0, 0]               # 給定搜索空間的下限，x1_min,..xn_min    xMax = [8, 8]               # 給定搜索空間的上限，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鏈長(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)化問(wèn)題的初始解 ======    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鏈    # ====== 開(kāi)始模擬退火優(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)生新解：通過(guò)在當(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 = True                kBetter += 1            else:  # 容忍解：如果新解的目標(biāo)函數(shù)比當(dāng)前解差，則以一定概率接受新解                pAccept = math.exp(-deltaE / tNow)  # 計(jì)算容忍解的狀態(tài)遷移概率                if pAccept > random.random():                    accept = True  # 接受劣質(zhì)解                    kBadAccept += 1                else:                    accept = False  # 拒絕劣質(zhì)解                    kBadRefuse += 1            # 保存新解            if accept == True:  # 如果接受新解，則將新解保存為當(dāng)前解                xNow[:] = xNew[:]                fxNow = fxNew                if fxNew < fxBest:  # 如果新解的目標(biāo)函數(shù)好于最優(yōu)解，則將新解保存為最優(yōu)解                    fxBest = fxNew                    xBest[:] = xNew[:]                    totalImprove += 1                    scale = 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 * alfa        kIter = kIter + 1        fxBest = cal_Energy(xBest, nVar, kIter)  # 由于迭代后懲罰因子增大，需隨之重構(gòu)增廣目標(biāo)函數(shù)        # ====== 結(jié)束模擬退火過(guò)程 ======    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!")        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    # 參數(shù)設(shè)置，優(yōu)化問(wèn)題參數(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()

輸出結(jié)果：

x_Initial:0.000000,4.000000,f(x_Initial):-36.000000i:0,t(i):100.00, badAccept:0.925373, f(x)_best:-152.000000i:10,t(i):81.71, badAccept:0.671053, f(x)_best:-98.000000i:20,t(i):66.76, badAccept:0.722892, f(x)_best:-98.000000i:30,t(i):54.55, badAccept:0.704225, f(x)_best:-98.000000i:40,t(i):44.57, badAccept:0.542169, f(x)_best:-98.000000i:50,t(i):36.42, badAccept:0.435294, f(x)_best:-98.000000i:60,t(i):29.76, badAccept:0.359551, f(x)_best:-98.000000i:70,t(i):24.31, badAccept:0.717647, f(x)_best:-98.000000i:80,t(i):19.86, badAccept:0.388235, f(x)_best:-98.000000i:90,t(i):16.23, badAccept:0.555556, f(x)_best:-98.000000i:100,t(i):13.26, badAccept:0.482353, f(x)_best:-98.000000i:110,t(i):10.84, badAccept:0.527473, f(x)_best:-98.000000i:120,t(i):8.85, badAccept:0.164948, f(x)_best:-98.000000i:130,t(i):7.23, badAccept:0.305263, f(x)_best:-98.000000i:140,t(i):5.91, badAccept:0.120000, f(x)_best:-98.000000i:150,t(i):4.83, badAccept:0.422680, f(x)_best:-98.000000i:160,t(i):3.95, badAccept:0.111111, f(x)_best:-98.000000i:170,t(i):3.22, badAccept:0.350000, f(x)_best:-98.000000i:180,t(i):2.63, badAccept:0.280000, f(x)_best:-98.000000i:190,t(i):2.15, badAccept:0.310000, f(x)_best:-98.000000i:200,t(i):1.76, badAccept:0.390000, f(x)_best:-98.000000i:210,t(i):1.44, badAccept:0.390000, f(x)_best:-98.000000i:220,t(i):1.17, badAccept:0.380000, f(x)_best:-98.000000improve:10Optimization by simulated annealing algorithm:x[0] = 8.0x[1] = 2.0f(x) = -98.0