模拟退火算法小谈
模拟退火算法：
为了解决局部最优解问题， 1983年，Kirkpatrick等提出了模拟退火算法（SA）能有效的解决局部最优解问题。我们知道在分子和原子的世界中，能量越大，意味着分子和原子越不稳定，当能量越低时，原子越稳定。
‘退火’是物理学术语，指对物体加温在冷却的过程。
模拟退火算法来源于晶体冷却的过程，如果固体不处于最低能量状态，给固体加热再冷却，随着温度缓慢下降，固体中的原子按照一定形状排列，形成高密度、低能量的有规则晶体，对应于算法中的全局最优解。
而如果温度下降过快，可能导致原子缺少足够的时间排列成晶体的结构，结果产生了具有较高能量的非晶体，这就是局部最优解。因此就可以根据退火的过程，给其在增加一点能量，然后在冷却，如果增加能量，跳出了局部最优解，这本次退火就是成功的.模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成.
这个算法来源于金属热处理中的退火（annealing）:
In condensed matter physics, annealing denotes a physical process in which a solid in a heat bathis heated up by increasing the temperature of the heat bath to a maximum value at which all particles of the solid randomly arrange themselves in the liquid phase, followed by cooling through slowly lowering the temperature of the heat bath. In this way, all particles arrange themselves in the low energy ground state of a corresponding lattice, provided the maximum temperature is sufficiently high and the cooling is carried out sufficiently slowly. 在凝聚态物理中，退火是指这样一个物理过程：将热浴的温度升高到最大值来加热热浴中的固体的，在该最大值处，所有固体颗粒随机排列在液相中，然后通过缓慢降低热浴的温度进行冷却。通过这种方式，只要最高温度足够高并且冷却进行得足够慢，所有颗粒会自身排列在相应晶格的低能量基态。
同时我们认为在每一个温度
(temperature)下，固体能够达到某个热平衡；为了模拟固体在某温度下达到热平衡的过程，我们引入Metropolis criterion：
假设固体此时的状态能量为
，我们细微地扰动固体的内在结构从而得到新的状态，其能量为
，为了后面表示方便，我们用
来表示扰动前后的能量差，然后我们接受这个扰动（即接受扰动后的状态）的概率（acceptance probability function）为
。不断大量地重复这个过程，我们认为该固体将会达到热平衡。在达到热平衡后，降低温度的值，重复上述过程，直到达到低能量基态。通过观察
，我们可以发现如果新的状态比旧的状态更“优”，则我们100%接受他，否则我们会以一定的概率接受新状态（往往比较小)。
那么从计算机科学的角度来说，我们往往把某一个问题的潜在解认做不同的状态；在搜索解空间的过程中，我们不断的比较当前解和新解的“能量”（需要引入目标函数energy/goal/objective function来计算该状态下的“能量”值）,根据acceptance probability function来决定接受哪一个解。降温的过程（cooling scheme）也就对应着调整
的过程。
基本的流程如下：
- 初始解，目标函数，cooling scheme，最大迭代次数；
- 由目标函数求得当前解的能量值；
- 随机扰动获得新解，并计算得到其能量值；
- 根据acceptance probability function来决定接受哪一个解；
- 按照cooling scheme调整参数；
- 将接受解作为1中的当前解，重复2-3，并记录迭代次数；
- 达到最大迭代次数（或满足其他退出条件）。
在最基本的模拟退火算法中，cooling scheme会根据迭代的次数不断降低
的值，以使得worse解的接受概率变小。显而易见的，cooling scheme也好，最大迭代次数也好，都会影响最终的结果，比方说温度下降太快，那么搜索过程会很快结束，同时得到局部最优解。同时，如何设定某一状态的neighborhood（即第1.2步，扰动即意味着选择另外一个“邻近”的解），如何设置目标函数等等也需要考虑。
总体来看，模拟退火通过引入概率来避免局部最优解，并尝试搜寻全局最优解，其基本想法和实现都比较简单。[1]
下面我们就详细讲讲他是如何在局部最优解跳出来到全局最优解的：
模拟退火算法包含两个部分即Metropolis算法和退火过程。Metropolis算法就是如何在局部最优解的情况下让其跳出来，是退火的基础。1953年Metropolis提出重要性采样方法，即以概率来接受新状态，而不是使用完全确定的规则，称为Metropolis准则，计算量较低。下面先形象的说一下，然后在用数学公式：
假设开始状态在A，随着迭代次数更新到B的局部最优解，这时发现更新到B时，能力比A要低，则说明接近最优解了，因此百分百转移，状态到达B后，发现下一步能量上升了，如果是梯度下降则是不允许继续向前的，而这里会以一定的概率跳出这个坑，这各概率和当前的状态、能量等都有关系，下面会详细说，如果B最终跳出来了到达C，又会继续以一定的概率跳出来，可能有人会迷惑会不会跳回之前的B呢？下面会解释，直到到达D后，就会稳定下来。所以说这个概率的设计是很重要的，下面从数学方面进行解释。
讲真,微信没有好的数学公式编辑器,好多公式就不放了,等我后面有自己的个人博客再放吧,微信太浪费时间了.
模拟退火算法可以分解为解空间、目标函数和初始解三部分。
模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率exp(-ΔT/T)接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。
退火算法程序流程图;
模拟退火算法的优缺点
- 迭代搜索效率高，并且可以并行化；
- 算法中有一定概率接受比当前解较差的解，因此一定程度上可以跳出局部最优；
- 算法求得的解与初始解状态S无关，因此有一定的鲁棒性；
- 具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法
说到了这里,有一个小重点:
只是接近最优解？
是的！对于一些问题，如TSP问题，如果把所有可能的解都遍历一遍，需要的时间是随着城市的数量增加而呈爆炸性增长的，因此对于多个城市的TSP问题，老老实实去解，会花费特别多的时间。若是对于解的精度没有太高的要求，则可以损失一些精度来换取速度。这个其实是很划算的,总比你没有的强,对吧.
再谈核心思想
模拟退火算法的核心思想是：首先随机选择一个解作为开始，接下来产生一个随机扰动，如果找到比上一个解更接近最优解的解，那么就直接接受这个解。而如果找到的解离得更远了，没关系，以一定的概率接受。
关于疑问
这个算法怎么感觉概率的地方基本不受控制，温度降到最后还是大量接受随机解,等我再学两年可能会解释吧.
上代码
import numpy as np
import matplotlib.pyplot as plt
import random
class SA(object):
def __init__(self, interval, tab='min', T_max=10000, T_min=1, iterMax=1000, rate=0.95):
self.interval = interval # 给定状态空间 - 即待求解空间
self.T_max = T_max # 初始退火温度 - 温度上限
self.T_min = T_min # 截止退火温度 - 温度下限
self.iterMax = iterMax # 定温内部迭代次数
self.rate = rate # 退火降温速度
#############################################################
self.x_seed = random.uniform(interval[0], interval[1]) # 解空间内的种子
self.tab = tab.strip() # 求解最大值还是最小值的标签: 'min' - 最小值；'max' - 最大值
#############################################################
self.solve() # 完成主体的求解过程
self.display() # 数据可视化展示
def solve(self):
temp = 'deal_' + self.tab # 采用反射方法提取对应的函数
if hasattr(self, temp):
deal = getattr(self, temp)
else:
exit('>>>tab标签传参有误："min"|"max"<<<')
x1 = self.x_seed
T = self.T_max
while T >= self.T_min:
for i in range(self.iterMax):
f1 = self.func(x1)
delta_x = random.random() * 2 - 1
if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]: # 将随机解束缚在给定状态空间内
x2 = x1 + delta_x
else:
x2 = x1 - delta_x
f2 = self.func(x2)
delta_f = f2 - f1
x1 = deal(x1, x2, delta_f, T)
T *= self.rate
self.x_solu = x1 # 提取最终退火解
def func(self, x): # 状态产生函数 - 即待求解函数
value = np.sin(x**2) * (x**2 - 5*x)
return value
def p_min(self, delta, T): # 计算最小值时，容忍解的状态迁移概率
probability = np.exp(-delta/T)
return probability
def p_max(self, delta, T):
probability = np.exp(delta/T) # 计算最大值时，容忍解的状态迁移概率
return probability
def deal_min(self, x1, x2, delta, T):
if delta < 0: # 更优解
return x2
else: # 容忍解
P = self.p_min(delta, T)
if P > random.random(): return x2
else: return x1
def deal_max(self, x1, x2, delta, T):
if delta > 0: # 更优解
return x2
else: # 容忍解
P = self.p_max(delta, T)
if P > random.random(): return x2
else: return x1
def display(self):
print('seed: {}\nsolution: {}'.format(self.x_seed, self.x_solu))
plt.figure(figsize=(6, 4))
x = np.linspace(self.interval[0], self.interval[1], 300)
y = self.func(x)
plt.plot(x, y, 'g-', label='function')
plt.plot(self.x_seed, self.func(self.x_seed), 'bo', label='seed')
plt.plot(self.x_solu, self.func(self.x_solu), 'r*', label='solution')
plt.title('solution = {}'.format(self.x_solu))
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.savefig('SA.png', dpi=500)
plt.show()
plt.close()
if __name__ == '__main__':
SA([-5, 5], 'max')
f(x)=(x2−5x)sin(x2
MATLAB版
%% SA 模拟退火: 求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值(动画演示）
clear; clc
%% 绘制函数的图形
x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);
figure(1)
plot(x,y,'b-')
hold on % 不关闭图形，继续在上面画图
%% 参数初始化
narvs = 1; % 变量个数
T0 = 100; % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 200; % 最大迭代次数
Lk = 100; % 每个温度下的迭代次数
alfa = 0.95; % 温度衰减系数
x_lb = -3; % x的下界
x_ub = 3; % x的上界
%% 随机生成一个初始解
x0 = x_lb + (x_ub - x_lb)*rand(1,narvs);
h = scatter(x0,Obj_fun1(x0),'*r'); % scatter是绘制二维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）
%% 初始化用来保存中间结果的x和y的取值
iter_x = x0;
iter_y = Obj_fun1(x0);
%% 模拟退火过程
for iter = 1 : maxgen % 外循环, 我这里采用的是指定最大迭代次数
for i = 1 : Lk % 内循环，在每个温度下开始迭代
% title(['当前温度：',num2str(T)])
y0 = Obj_fun1(x0); % 计算当前解的函数值
% disp('**************************分**割**线**************************')
% disp(['当前解的函数值：',num2str(y0)])
y = randn(1,narvs); % 生成1行narvs列的N(0,1)随机数
z = y / sqrt(sum(y.^2)); % 根据新解的产生规则计算z
x_new = x0 + z*T; % 根据新解的产生规则计算x_new的值
% 如果这个新解的位置超出了定义域，就对其进行调整
for j = 1: narvs
if x_new(j) < x_lb(j)
r = rand(1);
x_new(j) = r*x_lb(j)+(1-r)*x0(j);
elseif x_new(j) > x_ub(j)
r = rand(1);
x_new(j) = r*x_ub(j)+(1-r)*x0(j);
end
end
x1 = x_new; % 将调整后的x_new赋值给新解x1
y1 = Obj_fun1(x1); % 计算新解的函数值
% disp(['新解的函数值：',num2str(y1)])
if y1 > y0 % 如果新解函数值大于当前解的函数值
% disp('更新当前解为新解')
x0 = x1; % 更新当前解为新解
iter_x = [iter_x; x1]; % 将新找到的x1添加到中间结果中
iter_y = [iter_y; y1]; % 将新找到的y1添加到中间结果中
else
p = exp(-(y0 - y1)/T); % 根据Metropolis准则计算一个概率
% disp(['新解的函数值小于当前解的函数值，接受新解的概率为：',num2str(p)])
if rand(1) < p % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
% disp('该随机数小于这个概率，那么我们接受新解')
x0 = x1; % 更新当前解为新解
end
end
end
h.XData = x0; % 更新散点图句柄的x轴的数据（此时解的位置在图上发生了变化）
h.YData = Obj_fun1(x0); % 更新散点图句柄的y轴的数据（此时解的位置在图上发生了变化）
T = alfa*T; % 温度下降
pause(0.01) % 暂停一段时间(单位：秒)后再接着画图
end
[best_y, ind] = max(iter_y); % 找到最大的y的值，以及其在向量中的下标
best_x = iter_x(ind,:); % 根据下标找到此时x的值
disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(best_y)
pause(1)
h.XData = []; h.YData = []; % 将原来的散点删除
scatter(best_x,best_y,'*r'); % 在最大值处重新标上散点
title(['模拟退火找到的最大值为', num2str(best_y)]) % 加上图的标题
结果
写到这里还是底虚,因为这个东西属于概率论的东西,我正统概率论没有系统的学,所以只能找各种资料自学.就这样吧.日后有能力重构本文.