python模拟退火算法求解规划问题

import math
import random
import numpy as np
import pandas as pd
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"  # 定义问题名称
    nVar = 2  # 自变量数量 y=f（x1,x2,...xn)
    xMin = [0, 0]  #给定搜索空间的下限， x1_min, ... xn_min
    xMax = [8, 7.5 ]  # 给定搜索空间的上限 x1_max,..xn_max

    tInitial = 100.0  # 设置初始退火温度
    tFinal = 1  # 设置种子退火温度
    alfa = 0.98
    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)) #初始化创建数组 【0， 0】
    for v in range(nVar):
        # random.uniform(min, max) 在[min, max] 范围内随机生成一个实数
        xInitial[v] = random.uniform(xMin[v], xMax[v])# x1, x2

    #调用子函数， cal——Energy(xInitial, nVar, 1) # m(k) ， 惩罚因子，初值为1
    fxInitial = cal_Energy(xInitial, nVar, 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 # 总计Mrakov链长度
    totalImprove = 0 #fxBest改善次数
    nMarkov = meanMarkov  # 固定长度Markov链


    #=====开始模拟退火优化====
    #外循环，直到当前温度达到终止温度时结束
    tNow = tInitial  #初始化当前温度
    while tNow >=tFinal: #外循环，直到当前温度达到终止温度时结束  100 98  , 1
        # 在当前温度下，进行充分次数（nMarkov）的状态转移以达到热平衡
        kBetter = 0 #获得优值解的次数
        kBadAccept = 0 #接受劣质解的次数
        kBadRefuse = 0# 拒绝劣质解的次数

        #----内循环，循环次数为Markov链长度 100
        for k in range(nMarkov):  #内循环， 循环次数为Markov链长度
            totalMar +=1  #总Markov链长度计数器
        # x1,x2
            #----产生新解
            # 产生新解，通过在当前解附近随机扰动而产生新解 ，新解必须在[min,max]范围内
            #方案1，只对n元变量中的一个进行扰动，其他n-1个变量保持不变
            xNew[:] = xNow[:]
            v = random.randint(0, nVar-1) #产生[ 0, nVar]之间的随机数 8+0.5*
            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准则接受新解
            #接受判别： 按照Metropolics 准则决定是否接受新解

            if fxNew  random.random():
                    accept = True#接受劣质解
                    kBadAccept +=1

            #保存新解
            if accept == True: #如果接受新解， 则将新解保存为当前解
                xNow[:] = xNew[:]
                fxNow = fxNew
                if fxNew  1e-3:  #检验目标函数
        print("Error 2: Wrong total millage")
        return
    else:
        print("n Optomozation by simulated annealing algorithm")
        for i in range(nVar):
            print("t x[{}] = {:.6f}".format(i, xBest[i]))
        print("nt tf(x):{:.6f}".format(cal_Energy(xBest, nVar, 0)))

        #优化结果写入数据文件=====

    nowTime = datetime.now().strftime('%m%d%H%M')
    fileName = '{}_{}.dat'.format(cName, nowTime )
    optRecord = {
        'iter': recordIter,
        "FxNow": recordFxNow,
        'FxBest': recordxBest,
        "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('n Optimization by simulated annealing algorithm:')
        for i in range(nVar):
            fid.write('ntx[{}] = {:.6f}'.format(i,xBest[i]))

        fid.write('ntf(x):{:.6f}'.format(cal_Energy(xBest, nVar, 0)))
    print("写入数据文件：%s完成。" %fileName)

    #0-----优化结果图形输出
    plt.figure(figsize=(6,4), facecolor='#FFFFFF')
    plt.title("Optimization result: {}".format(cName))
    plt.xlim((0, kInter))
    plt.xlabel("iter")
    plt.ylabel("f(x)")
    plt.plot(recordIter, recordFxNow, 'b-', label='FxNow')
    plt.plot(recordIter, recordxBest, 'r-', label='FxBest')
    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)


if __name__ == '__main__':
    main()