1.算法简介
模拟退火算法（Simulated Annealing，SA）是一种模拟物理退火的过程而设计的随机优化算法，结合爬山法和随机行走算法，同时避免算法进入局部最优，早期用于组合优化，后来发展成一种通用的优化算法。它的基本思想最早在1953年就被Metropolis提出，但直到1983年Kirkpatrick等人才设计出真正意义上的模拟退火算法并进行应用。
该算法采用类似于物理退火的过程，先在一个高温状态下（相当于算法随机搜索），然后逐渐退火，在每个温度下（相当于算法的每一次状态转移），徐徐冷却（相当于算法局部搜索），最终达到物理基态（相当于算法找到最优解）。
高温过程——增强粒子的热运动，使其偏离平衡位置，目的是消除系统原先可能存在的非均匀态；
等温过程——退火过程中要让温度慢慢降低，在每一个温度下要达到热平衡状态，对于与环境换热而温度不变的封闭系统满足自由能较少定律，系统状态的自发变化总是朝自由能减少的方向进行，当自由能达到最小时，系统达到平衡态；
冷却过程——使粒子热运动减弱并渐趋有序，系统能量逐渐下降，从而得到低能的晶体结构。当液体凝固为固体的晶态时退火过程完成。
因此模拟退火算法从某一高温出发，在高温状态下计算初始解，然后以预设的邻域函数产生一个扰动量，从而得到新的状态，即模拟粒子的无序运动，比较新旧状态下的能量，即目标函数的解。如果新状态的能量小于旧状态，则状态发生转化；如果新状态的能量大于旧状态，则以Metropolis接受准则发生转化。当状态稳定后，便可以看作达到了当前状态的最优解，便可以开始降温，在下一个温度继续迭代，最终达到低温的稳定状态，便得到了模拟退火算法产生的结果
该算法的关键点如下：
1、对固体退火过程的模拟；
2、采用Metropolis接受准则；
3、用冷却进度表控制算法进程，使算法在多项式时间里给出一个近似解。
固体退火过程是SAA的物理背景；Metropolis接受准则使算法跳离局部最优 “险井”；而冷却进度表的合理选择是算法应用的前提。
下面依次介绍上述三个关键点
1.1 固体退火过程：
对于该部分如果不理解，可以先简单跳过。
固体退火过程的数学表述：
在温度T，分子停留在状态r满足Boltzmann概率分布
1.2 Metropolis准则
Metropolis准则（1953）——以概率接受新状态，固体在恒定温度下达到热平衡的过程可以用Monte Carlo方法（计算机随机模拟方法）加以模拟。
Monte Carlo模拟退火过程：蒙特卡罗(Monte Carlo)方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。
若在温度T，当前状态i → 新状态j
若Ej<Ei，则接受 j 为当前状态；
否则，若概率 p=e x p [ − ( E j − E i ) / k B T ] exp[-(Ej-Ei)/k_BT]exp[−(Ej−Ei)/k B
T](也就是状态j与状态i概率的比值) 大于[0,1)区间的随机数，则仍接受状态 j 为当前状态；若不成立则保留状态 i 为当前状态。
可以看出，高温下可接受与当前状态能差较大的状态为重要状态，而在低温下只能接受与当前状态能差较小的新状态为重要状态。在温度趋于零时，就不再接受Ej>Ei的新状态j了.如此反复，达到系统在此温度下的热平衡。这个过程称作Metropolis抽样过程。上面这种方法被称为重要性采样法，该方法以概率形式接受新状态，避免局部最优。
1.3 冷却进度表
在Metropolis抽样过程中温度T缓慢的降低，模拟退火过程就是通过T参数的变化使状态收敛于最小能量处。因而，T参数的选择对于算法最后的结果有很大影响。初始温度和终止温度设置的过低或过高都会延长搜索时间。降温步骤太快，往往会漏掉全局最优点，使算法收敛至局部最优点。降温步骤太慢，则会大大延长搜索全局最优点的计算时间，从而难以实际应用。因此选择一个合适的降温函数非常重要，冷却进度表就是指的从某一高温状态T向低温状态冷却时的降温函数。
import math
from random import random
import matplotlib.pyplot as plt
"""
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率P接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
"""
def func(x, y): # 函数优化问题
res = 4 * x ** 2 +2*y**2
return res
# x为公式里的x1,y为公式里面的x2
class SA:
def __init__(self, func, iter=100, T0=100, Tf=0.01, alpha=0.99):
self.func = func
self.iter = iter # 内循环迭代次数,即为L =100
self.alpha = alpha # 降温系数，alpha=0.99
self.T0 = T0 # 初始温度T0为100
self.Tf = Tf # 温度终值Tf为0.01
self.T = T0 # 当前温度
self.x = [random() * 11 - 5 for i in range(iter)] # 随机生成100个x的值列表形式
self.y = [random() * 11 - 5 for i in range(iter)] # 随机生成100个y的值
self.most_best = []
"""
random()这个函数取0到1之间的小数
如果你要取0-10之间的整数（包括0和10）就写成 (int)random()*11就可以了，11乘以零点多的数最大是10点多，最小是0点多
该实例中x1和x2的绝对值不超过5（包含整数5和-5），（random() * 11 -5）的结果是-6到6之间的任意值（不包括-6和6）
（random() * 10 -5）的结果是-5到5之间的任意值（不包括-5和5），所有先乘以11，取-6到6之间的值，产生新解过程中，用一个if条件语句把-5到5之间（包括整数5和-5）的筛选出来。
"""
self.history = {'f': [], 'T': []}
def generate_new(self, x, y): # 扰动产生新解的过程
while True:
x_new = x + self.T * (random() - random())
y_new = y + self.T * (random() - random())
if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):
break # 重复得到新解，直到产生的新解满足约束条件
return x_new, y_new
def Metrospolis(self, f, f_new): # Metropolis准则，新的小于旧的
if f_new <= f:
return 1
else:
p = math.exp((f - f_new) / self.T) #返回e^x,x表示参数,e是欧拉常数
if random() < p:
return 1
else:
return 0
def best(self): # 获取最优目标函数值
f_list = [] # f_list数组保存每次迭代之后的值
for i in range(self.iter):
f = self.func(self.x[i], self.y[i])
f_list.append(f)
f_best = min(f_list)
idx = f_list.index(f_best)
return f_best, idx # f_best,idx分别为在该温度下，迭代L次之后目标函数的最优解和最优解的下标
def run(self):
count = 0
# 外循环迭代，当前温度小于终止温度的阈值
while self.T > self.Tf: #T根据降温系数不断减小，但为防止直接达到最优,内部循环100次
# 内循环迭代100次
for i in range(self.iter):
f = self.func(self.x[i], self.y[i]) # f为迭代一次后的值
x_new, y_new = self.generate_new(self.x[i], self.y[i]) # 产生新解
f_new = self.func(x_new, y_new) # 产生新值
if self.Metrospolis(f, f_new): # 判断是否接受新值
self.x[i] = x_new # 如果接受新值，则把新值的x,y存入x数组和y数组
self.y[i] = y_new
# 迭代L次记录在该温度下最优解
ft, _ = self.best()
self.history['f'].append(ft)
self.history['T'].append(self.T)
# 温度按照一定的比例下降（冷却）
self.T = self.T * self.alpha
count += 1
# 得到最优解
f_best, idx = self.best() #分别为在该温度下，迭代L次之后目标函数的最优解和最优解的下标
print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")
print(count)
sa = SA(func)
sa.run()
plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()结果：F=3.591534016505581e-05, x=0.0010217709965204412, y=-0.003983671440250457
def demo2_func(x):
x1, x2, x3 = x
return x1**2+x2**2+x3**2
#模拟退火算法
from sko.SA import SAFast
sa_fast = SAFast(func=demo2_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_fast.run()
print('Fast Simulated Annealing: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)
