0. 写在前面
超级简单的模拟退火算法实现ε٩(๑> ₃
1. 模型描述及处理
1.1 线性规划模型
[max,f(x)=10x_1+9x_2
]
]
(s.t.)
[6x_1+5x_2leq{60}tag{1}
]
]
[10x_1+20x_2leq{150}tag{2}
]
]
[0leq{x_1}leq{8}tag{3}
]
]
[0leq{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)quadforall{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 random.random():
accept = True # 接受劣质解
kBadAccept += 1
else:
accept = False # 拒绝劣质解
kBadRefuse += 1
# 保存新解
if accept == True: # 如果接受新解，则将新解保存为当前解
xNow[:] = xNew[:]
fxNow = fxNew
if fxNew 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('ntf(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
服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net
机房租用，北京机房租用，IDC机房托管， http://www.e1idc.net
