Python數模筆記-模擬退火算法（1）多變量函數優化

時間 2021-07-23

標籤 html python 算法數組 app dom 函數性能測試優化欄目 Python 简体版

原文原文鏈接

一、模擬退火算法

　　模擬退火算法借鑑了統計物理學的思想，是一種簡單、通用的啓發式優化算法，並在理論上具備機率性全局優化性能，於是在科研和工程中獲得了普遍的應用。
　　退火是金屬從熔融狀態緩慢冷卻、最終達到能量最低的平衡態的過程。模擬退火算法基於優化問題求解過程與金屬退火過程的類似性，以優化目標爲能量函數，以解空間爲狀態空間，以隨機擾動模擬粒子的熱運動來求解優化問題（[1] KIRKPATRICK,1988）。
　　模擬退火算法結構簡單，由溫度更新函數、狀態產生函數、狀態接受函數和內循環、外循環終止準則構成。html

　　溫度更新函數是指退火溫度緩慢下降的實現方案，也稱冷卻進度表；
　　狀態產生函數是指由當前解隨機產生新的候選解的方法；
　　狀態接受函數是指接受候選解的機制，一般採用Metropolis準則；
　　外循環是由冷卻進度表控制的溫度循環；
　　內循環是在每一溫度下循環迭代產生新解的次數，也稱Markov鏈長度。python

　　模擬退火算法的基本流程以下：算法

　　（1）初始化：初始溫度T，初始解狀態s，迭代次數L；
　　（2）對每一個溫度狀態，重複 L次循環產生和機率性接受新解：
　　（3）經過變換操做由當前解s 產生新解s′；
　　（4）計算能量差 ∆E，即新解的目標函數與原有解的目標函數的差；
　　（5）若∆E <0則接受s′做爲新的當前解，不然以機率exp(-∆E/T) 接受s′ 做爲新的當前解；
　　（6）在每一個溫度狀態完成 L次內循環後，下降溫度 T，直到達到終止溫度。數組

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

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

二、多變量函數優化問題

　　選取經典的函數優化問題和組合優化問題做爲測試案例。函數

　　問題 1：Schwefel 測試函數，是複雜的多峯函數，具備大量局部極值區域。
　　F(X)=418.9829×n-∑(i=1,n)〖xi* sin⁡(√(|xi|)) 〗性能

　　本文取 d=10, x=[-500,500]，函數在 X=(420.9687,...420.9687)處爲全局最小值 f(X)=0.0。測試

　　使用模擬退火算法的基本方案：控制溫度按照 T(k) = a * T(k-1) 指數衰減，衰減係數取 a；如式（1）按照 Metropolis 準則接受新解。對於問題 1（Schwefel函數），經過對當前解的一個自變量施加正態分佈的隨機擾動產生新解。優化

三、模擬退火算法 Python 程序

# 模擬退火算法 程序：多變量連續函數優化
# Program: SimulatedAnnealing_v1.py
# Purpose: Simulated annealing algorithm for function optimization
# v1.0:
#   (1) 基本算法：單變量連續函數優化問題
#   (2) 文件輸出優化結果和中間過程數據
#   (3) 設置指標參數計數器
#   (4) 圖形輸出壞解接受機率
# Copyright 2021 YouCans, XUPT
# Crated：2021-04-30

#  -*- coding: utf-8 -*-
import math                         # 導入模塊
import random                       # 導入模塊
import pandas as pd                 # 導入模塊 YouCans, XUPT
import numpy as np                  # 導入模塊 numpy, 並簡寫成 np
import matplotlib.pyplot as plt     # 導入模塊 matplotlib.pyplot, 並簡寫成 plt
from datetime import datetime


# 子程序：定義優化問題的目標函數
def cal_Energy(X, nVar):
    # 測試函數 1： Schwefel 測試函數
    # -500 <= Xi <= 500
    # 全局極值：(420.9687,420.9687,...),f(x)=0.0
    sum = 0.0
    for i in range(nVar):
        sum += X[i] * np.sin(np.sqrt(abs(X[i])))
    fx = 418.9829 * nVar - sum
    return fx


# 子程序：模擬退火算法的參數設置
def ParameterSetting():
    cName = "funcOpt"           # 定義問題名稱
    nVar = 2                    # 給定自變量數量，y=f(x1,..xn)
    xMin = [-500, -500]         # 給定搜索空間的下限，x1_min,..xn_min
    xMax = [500, 500]           # 給定搜索空間的上限，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):
        # random.uniform(min,max) 在 [min,max] 範圍內隨機生成一個實數
        xInitial[v] = random.uniform(xMin[v], xMax[v])
    # 調用子函數 cal_Energy 計算當前解的目標函數值
    fxInitial = cal_Energy(xInitial, nVar)

    # ====== 模擬退火算法初始化 ======
    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)
            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
        # ====== 結束模擬退火過程 ======

    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)
    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[{}] = {:.6f}'.format(i,xBest[i]))
        print('\n\tf(x):{:.6f}'.format(fxBest))

    # ====== 優化結果寫入數據文件 ======
    nowTime = datetime.now().strftime('%m%d%H%M')       # '02151456'
    fileName = "..\data\{}_{}.dat".format(cName,nowTime)# 數據文件的地址和文件名
    optRecord = {
        "iter":recordIter,
        "FxNow":recordFxNow,
        "FxBest":recordFxBest,
        "PBad":recordPBad}
    df_Record = pd.DataFrame(optRecord)
    df_Record.to_csv(fileName, index=False, encoding="utf_8_sig")
    with open(fileName, 'a+', encoding="utf_8_sig") as fid:
        fid.write("\nOptimization by simulated annealing algorithm:")
        for i in range(nVar):
            fid.write('\n\tx[{}] = {:.6f}'.format(i,xBest[i]))
        fid.write('\n\tf(x):{:.6f}'.format(fxBest))
    print("寫入數據文件: %s 完成。" % fileName)

    # ====== 優化結果圖形化輸出 ======
    plt.figure(figsize=(6, 4), facecolor='#FFFFFF')     # 建立一個圖形窗口
    plt.title('Optimization result: {}'.format(cName))  # 設置圖形標題
    plt.xlim((0, kIter))                                # 設置 x軸範圍
    plt.xlabel('iter')                                  # 設置 x軸標籤
    plt.ylabel('f(x)')                                  # 設置 y軸標籤
    plt.plot(recordIter, recordFxNow,'b-', label='FxNow')     # 繪製 FxNow 曲線
    plt.plot(recordIter, recordFxBest, 'r-', label='FxBest')  # 繪製 FxBest 曲線
    # plt.plot(recordIter,recordPBad,'r-',label='pBadAccept')  # 繪製 pBadAccept 曲線
    plt.legend()  # 顯示圖例
    plt.show()

    return


# 主程序
def main():

    # 參數設置，優化問題參數定義，模擬退火算法參數設置
    [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://www.cnblogs.com/youcans/ ===
if __name__ == '__main__':
    main()

版權說明：
原創做品
Copyright 2021 YouCans, XUPT
Crated：2021-05-01

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