模拟退火算法+Python实现与可视化

本文借鉴1：模拟退火算法详细讲解（含实例python代码）
本文借鉴2：清风建模模拟退火讲解

一、退火的物理过程：

1. 加温过程。其目的是增强粒子的热运动，使其偏离平衡位置。当温度足够高时，固体将熔为液体，从而消除系统原先存在的非均匀状态。
2. 等温过程。对于与周围环境交换热量而温度不变的封闭系统，系统状态的自发变化总是朝自由能减少的方向进行的，当自由能达到最小时，系统达到平衡状态。
3. 冷却过程。使粒子热运动减弱，系统能量下降，得到晶体结构。

总结：
先在一个高温状态下（设定初值），然后逐渐退火（控制温度T下降），在每个温度下慢慢冷却（迭代+Metropolis算法），最终达到物理基态（相当于算法找到最优解）。

其中Metropolis 准则是SA算法收敛于全局最优解的关键所在——Metropolis准则以一定的概率接受恶化解，这样就使算法跳离局部最优的陷阱。

二、退火算法的步骤：

求解函数最小值问题：

上面包含两个部分：Metropolis算法和退火过程,，分别对应内循环和外循环。

三.退火算法的过程：

在一切开始之前先在区间里随机取一个初值x(0)。
下面进入内循环。

1.内循环

内循环就是Metropolis算法：上面步骤的(3)至(5)步。即在每次温度下，迭代L次，寻找在该温度下能量的最小值（即最优解）。

1. 扰动产生新解
   迭代过程的新解x（n+1）是由x(n)随机扰动产生（形象理解就是在x(n)基础上向左或向右随机移动），并且在温度越高的时候，移动的幅度越大，这模拟了物理中粒子的无序运动。

2. Metropolis准则
   （1）在固定温度下，整个迭代过程中温度不发生变化，能量发生变化，当前一个状态x(n)的能量大于后一个状态x(n+1)的能量时，说明状态x(n)的解没有状态x(n+1)的解好，所以接受状态x(n+1)。（这里说的能量即为实际问题中的目标函数值）

   （2）但是如果下一状态的能量比前一个状态的能量高时，该不该接受下一状态呢？在这里设置一个接受概率P,即如果下一状态的能量比前一个状态的能量高（即比前一个解差），则接受下一状态的概率为P,这就是Metropolis准则。

具体实现的时候，会随机生成一个在区间[0,1]上服从均匀分布的随机数r, 如果r<p则接受新解

3. 关于概率P

（1）温度一定时， ∆E越小，概率越大，即目标函数相差越小接受的可能性越大。
（2）∆E一定时，温度越高，概率越大，即搜索前期温度较高时更有可能接受新解。

当迭代次数足够时，进入外循环

2.外循环

外循环就是退火过程，将固体达到较高的温度（初始温度T（0）），控制温度不断下降，当达到终止温度Tf时，冷却结束，即退火过程结束

1. 冷却进度表
   
   （1）初始温度最常取100。
   （2）衰减函数最常用的是：

（3）这里说的停止准则是指当温度下降到一个特别小的数时就退出外循环。（T<1e-2） （4）这里提到了Markov链的长度, 这个是随机过程中马尔科夫链里面的概念，它和我们前面步骤(5): “在温度T下，将步骤(3)和(4)重复L次; ” 是一个意思,即在当前温度下需要内循环的次数。从经验来说，对于简单的情况可以令Lk=100n，其中n为自变量的个数

PYTHON代码实现

整个代码是我在浏览其他博友时看到的，并非我本人原创，只是我比较困惑内循环迭代为何是先随机生成L个点，每个迭代一次最终取最好的，而不是选取一个点迭代L次。

所以我按照后者改了一下（如果有明白为何按前者来做的朋友麻烦给我解答一下疑惑）

最后还画了能量随温度下降曲线和整个过程的一个动画图。

import math
from random import random
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation


# 目标函数
def func(x):
    """
    函数优化问题: min f(x)=x^4-4*x^2-2，范围[-5,5]
    """
    res = x ** 4 - 14 * x ** 2 + +10 * x + 4
    return res

#  SA类
class SA:
    """模拟退火算法的类"""

    def __init__(self, func, iter=100, T0=100, Tf=0.01, alpha=0.95):
        self.func = func
        self.iter = iter  # 内循环迭代次数
        self.alpha = alpha  # 降温系数
        self.T0 = T0  # 初始温度
        self.Tf = Tf  # 温度终值
        self.T = T0  # 当前温度
        self.x = random() * 4 - 2  # 随机生成(-2,2)之间的一个值
        self.most_best = []
        self.history = {'x': [], 'f': [], 'T': [],'f1':[],'T1':[]}  # 画图

    def generate_new(self, x):  # 扰动产生新解的过程
        while True:
            x_new = x + self.T * (random() - random())
            if -5 <= x_new <= 5:
                break  # 重复得到新解，直到产生的新解满足约束条件
        return x_new

    # Metropolis准则
    def Metrospolis(self, f, f_new):
        if f_new <= f:
            return 1  # 1表示接受，0表示拒绝
        else:
            p = math.exp((f - f_new) / self.T)
            if random() < p:  # random()生成一个[0,1)之间的随机数
                return 1
            else:
                return 0

    # 外循环迭代
    def run(self):
        while self.T > self.Tf:  # 当温度大于终止温度的阈值
            # 内循环迭代L次
            for i in range(self.iter):
                f = self.func(self.x)  # 旧值
                x_new = self.generate_new(self.x)  # 产生新解
                f_new = self.func(x_new)  # 新值
                if self.Metrospolis(f, f_new):  # 判断是否接受新值
                    self.x = x_new  # 如果接受新解，则把新解存入x
                ft = self.func(self.x)  # 迭代L次后最优值
                xt = self.x
                Tt = self.T
                self.history['x'].append(xt)
                self.history['f'].append(ft)
                self.history['T'].append(Tt)
            self.history['f1'].append(ft)
            self.history['T1'].append(Tt)
            # 温度按照一定的比例下降（冷却）
            self.T = self.T * self.alpha
        # 得到最优解
        print(f"F={ft}, x={xt}")


# 运行程序
sa = SA(func)
sa.run()

# 能量变化
plt.plot(sa.history['T1'], sa.history['f1'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()

# # 绘制随机点移动图像
def update_points(num):
    """
    更新数据点
    """
    point_ani.set_data(X[num], Y[num])
    text_pt.set_text("T=%.3f" % T[num])
    return point_ani, text_pt,


x = np.linspace(-5, 5, 100)  # 生成[0,10]范围内，100个等距离的数组
y = func(x)
X = sa.history['x']
Y = sa.history['f']
T = sa.history['T']

fig = plt.figure(tight_layout=True)
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('f(x)')

point_ani, = plt.plot(X[0], Y[0], "ro")
plt.grid(ls="--")
text_pt = plt.text(4, 0.8, '', fontsize=10)

ani = animation.FuncAnimation(fig, update_points, np.arange(0, len(X)), interval=0.01, blit=True, repeat=False)

plt.show()