Python数模笔记-模拟退火算法（1）多变量函数优化
1、模拟退火算法
模拟退火算法借鉴了统计物理学的思想，是一种简单、通用的启发式优化算法，并在理论上具有概率性全局优化性能，因而在科研和工程中得到了广泛的应用。
退火是金属从熔融状态缓慢冷却、最终达到能量最低的平衡态的过程。模拟退火算法基于优化问题求解过程与金属退火过程的相似性，以优化目标为能量函数，以解空间为状态空间，以随机扰动模拟粒子的热运动来求解优化问题（[1] KIRKPATRICK,1988）。
模拟退火算法结构简单，由温度更新函数、状态产生函数、状态接受函数和内循环、外循环终止准则构成。
温度更新函数是指退火温度缓慢降低的实现方案，也称冷却进度表；
状态产生函数是指由当前解随机产生新的候选解的方法；
状态接受函数是指接受候选解的机制，通常采用Metropolis准则；
外循环是由冷却进度表控制的温度循环；
内循环是在每一温度下循环迭代产生新解的次数，也称Markov链长度。
模拟退火算法的基本流程如下：
（1）初始化：初始温度T，初始解状态s，迭代次数L；
（2）对每个温度状态，重复 L次循环产生和概率性接受新解：
（3）通过变换操作由当前解s 产生新解s′；
（4）计算能量差 ∆E，即新解的目标函数与原有解的目标函数的差；
（5）若∆E <0则接受s′作为新的当前解，否则以概率exp(-∆E/T) 接受s′ 作为新的当前解；
（6）在每个温度状态完成 L次内循环后，降低温度 T，直到达到终止温度。
欢迎关注 Youcans 原创系列，每周更新数模笔记
Python数模笔记-PuLP库
Python数模笔记-StatsModels统计回归
Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法
2、多变量函数优化问题
选取经典的函数优化问题和组合优化问题作为测试案例。
问题 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函数），通过对当前解的一个自变量施加正态分布的随机扰动产生新解。
3、模拟退火算法 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 原创系列，每周更新数模笔记
Python数模笔记-PuLP库（1）线性规划入门
Python数模笔记-PuLP库（2）线性规划进阶
Python数模笔记-PuLP库（3）线性规划实例
Python数模笔记-Scipy库（1）线性规划问题
Python数模笔记-StatsModels 统计回归（1）简介
Python数模笔记-StatsModels 统计回归（2）线性回归
Python数模笔记-StatsModels 统计回归（3）模型数据的准备
Python数模笔记-StatsModels 统计回归（4）可视化
Python数模笔记-Sklearn （1）介绍
Python数模笔记-Sklearn （2）聚类分析
Python数模笔记-Sklearn （3）主成分分析
Python数模笔记-Sklearn （4）线性回归
Python数模笔记-Sklearn （5）支持向量机
Python数模笔记-NetworkX（1）图的操作
Python数模笔记-NetworkX（2）最短路径
Python数模笔记-NetworkX（3）条件最短路径
Python数模笔记-模拟退火算法（1）多变量函数优化
Python数模笔记-模拟退火算法（2）约束条件的处理
Python数模笔记-模拟退火算法（3）整数规划问题
Python数模笔记-模拟退火算法（4）旅行商问题
