模拟退火算法Python

1.是什么？为什么？怎么做？

Q1:模拟退火是什么算法？

模拟退火是模拟物理上退火方法，通过N次迭代（退火），逼近函数的上的一个最值（最大或者最小值）。

比如逼近这个函数的最大值C点：

Q2:模拟退火为什么可行？

讨论这个问题需要理解一下物理原型是怎么样的，也就是原来是怎么“退火”的：

模拟退火算法的思想借鉴于固体的退火原理，当固体的温度很高的时候，内能比较大，固体的内部粒子处于快速无序运动，当温度慢慢降低的过程中，固体的内能减小，粒子的慢慢趋于有序，最终，当固体处于常温时，内能达到最小，此时，粒子最为稳定。

注意标粗字体：

1. 温度高->运动速度快（温度低->运动速度慢）
2. 温度是缓慢（想象成特别慢的那种）降低的
3. 温度基本不再变化后，趋于有序(最后内能达到最小，也就是接近最优)

我们通过模拟这个操作，使得我们需要的答案“趋于有序”，也就是最靠近需要的值（最值）。

Q3:怎么做？

大方向

首先，理解一下大方向：

模拟退火就是一种循环算法。

1. 我们先设定一个初始的温度$T$（这个温度会比较高，比如2000）
2. 每次循环都退火一次。（具体怎么操作后面详解）
3. 然后降低$T$的温度，我们通过让$T$和一个“降温系数”$DeltaT$（一个接近1的小数，比如$0.99$)相乘，达到慢慢降低温度的效果，直到接近于0（我们用$eps$来代表一个接近0的数(比如0.00001)，只要$T>eps$就可以退出循环了）

所以总的来说，用伪代码表示退火的流程是这样的：

T = 2000 # 代表开始的温度
dT = 0.99 # 代表系数delta T
eps = 1e-14 # 相当于0.0000000000000001
while T > eps：
    # --------------
    # 这里是每一次退火的操作
    # --------------
    T = T * dT # 温度每次下降一点点， T * 0.99

退火详解

我们要求解的答案无非是两个：自变量$x$和对应函数的最大值$f(x)$

那么模拟退火的是怎么做到的呢？【下面车速很快，请系好安全带！】

1. 我们先随机找一点$x\_0$ ，不论是哪个点都可以，随机！（不超过定义域就行）。

   这个点作为我们的初始值（相当于物体里的一个粒子）

2. 再找到一点$f(x\_0)$,来代表$x\_0$所对应的函数值

3. 现在正式开始退火！

   刚才我们说了$x\_0$相当于是一个粒子，所以我们会进行一个无序运动，也就是向左或者向右随机移动

   是的，是随机移动，可能向左，也可能向右，但是请记住一个关键点：移动的幅度和当前的温度$T$有关。

   温度$T$越大，移动的幅度越大。温度$T$越小，移动的幅度就越小。这是在模拟粒子无序运动的状态。

4. 接受(Accept)更"好"的状态

假设我们移动到了$x\_1$处，那么这个点对应的$f(x\_1)$很明显答案是优于（大于）当前的$f(x\_0)$的

因此我们将答案进行更新。也就是将初始值进行替换：$x\_0=x\_1,f(x\_0)=f(x\_1)$。这是一种贪心的思想。

5. 以一定概率接受(Accept)更差的状态

   这是退火最精彩的部分。

   为什么我们要接受一个更加差的状态呢？因为可能在一个较差的状态旁边会出现一个更加高的山峰

如果我们鼠目寸光，只盯着右半区，很容易随着温度的下降、左右跳转幅度的减小而迷失自己，最后困死在小山丘中。
而我们如果找到了左边山峰的低点，以一定的概率接受了它（概率大小和温度以及当前的值的关键程度有关），会在跳转幅度减少之前，尽可能找到最优点。
那么我们以多少的概率去接受它呢？我们用一个公式表示（这个公式我们只需记住，这是科学家推导出来的结论）：
$$\begin{gathered} Huge{e}^{ \\ frac \\ DeltafkT} \end{gathered}$$
别慌！很简单！我们来理解一下这里面的变量：
1. $e$是自然对数，约等于2.71。我们可以把右上角这一坨值$$\begin{gathered} frac \\ DeltafkT \end{gathered}$$看成一个整体$x$:
$e^x$的图形画出来是这样的：

因为我们想要函数$e^x$来代表一个概率值，所以我们只需要关注x为负数的部分即可：

负数部分的值域是在$(0,1)$开区间内，x越小，越接近0，越大越靠近1。

因为在0到1之间，所以这个值相当于是概率了。比如$e^x=0.97$，那么我们接受的概率就是$97$

而正数部分的值域会大于1，也就是说概率会超过$100$，所以会一定选（其实是上一种找到更优的情况）

2. $kT$

   $k$其实是个物理学常数，我们在代码中不会用到。

   $T$很简单，就是当前的温度。所以实际上这个分母就是$T$，$k$当做$1$使用。

3. $Deltaf$

我们着重讲一下什么是$Deltaf$。

其实从前面的函数$e^x$中可以发现，$Deltaf$必须是个负数！

我们想要函数$e^x$来代表一个概率值，一定要让它的值域属于$(0,1)$，所以$$\begin{gathered} frac \\ DeltafkT \end{gathered}$$必须是个负数。但是$kT$在我们的模拟中一定是正数，那么$Deltaf$必须是个负数！

其实$Deltaf$就是当前解的函数值与目标解函数值之差，$$\begin{gathered} Deltaf=- \\ left|f \\ left(x\_0 \\ right)-f \\ left(x\_1 \\ right) \\ right| \end{gathered}$$，并且一定是个负数。这个需要具体问题具体分析。

比如现在我们求一个函数的最大值，那么如果$f(x\_0)<f(x\_1)$了，那就说明结果变好了，我们肯定选择它（见第4点）

如果$f(x\_0)>f(x\_1)$，那就说明结果变差了，我们需要概率选择它，因此$Deltaf=-(f(x\_0)–f(x\_1))$

所以总结一下就是：

随机后的函数值如果结果更好，我们一定选择它(即$x\_0=x\_1,f(x\_0)=f(x\_1)$)
随机后的函数值如果结果更差，我们以$$\begin{gathered} LARGE{e}^{ \\ frac \\ DeltafkT} \end{gathered}$$的概率接受它

伪代码流程

注：对代码中的函数作出解释：
import random
①对random函数

1. random.randint(0,32767)拿到[0,32767]内的随机整数

2. random.randint(0,32767)2的范围是[0,32767 2]

3. random.randint(0,32767)*2-32767 的范围是[-32767, 32767]

import math
②对math.exp()函数

1. math.exp(x)代表$e^x$

③关于math.exp((f - df) / T) * 32767 > random.randint(0, 32767)

1. 目的是要概率接受，但是$e^x$是个准确值，所以从理论上我们可以生成一个（0,1）的随机数，如果$e^x$比(0,1)这个随机数要大，那么我们就接受。
2. 但是由于rand()只能产生[0,32767]内的随机整数，化成小数太过麻烦。所以我们可以把左边乘以32767（也就是把概率同时扩大32767倍），效果就等同于$e^x$比(0,1)了。

import random
import math

T = 2000 # 代表开始的温度
dT = 0.99 # 代表系数delta T
eps = 1e-14 # 相当于0.0000000000000001

# 用自变量计算函数值,这里可能存在多个自变量对应一个函数值的情况，比如f(x,y)
def func(x, ... ) :
    # 这里是对函数值进行计算
    ans = .......
    return ans
# 原始值
x = random.randint(0,32767) # x0取随机值
f = func(x,...) # 通过自变量算出f(x0)的值
while T > eps: 
    # --------------
    # 这里是每一次退火的操作

    # x1可以左右随机移动，幅度和温度T正相关，所以*T
    # 注意这里移动可以左右移动，但是也可以单向移动
    dx = x + (random.randint(0, 32767) * 2 - 32767) * T 

    # 让x落在定义域内，如果没在里面，就重新随机。题目有要求需要写，否则不用写
    # ================
    while x > ? || x < ? ... :
        dx = (2*rand() - RAND_MAX) * T 
    # ================

    # 求出f(x1)的值
    df = func(dx)
    # 这里需要具体问题具体分析，我们要接受更加优秀的情况。可能是df < f(比如求最小值）
    if f < df :
        f = df
        x = dx  [...,y = dy] # 接受，替换值，如果多个自变量，那么都替换
    # 否则概率接受，注意这里df-f也要具体问题具体分析。
    elif math.exp((f - df) / T) * 32767 > random.randint(0, 32767):
        f = df
        x = dx  [...y = dy] # 接受，替换值，如果多个自变量，那么都替换
    # --------------
    T = T * dT # 温度每次下降一点点， T * 0.99

2.通过模拟退火算出 根号n的值

思路是这样的：我们试图通过退火找出一个值$x\_0$，使得${x\_0}^2$的值更加接近于$$\begin{gathered} { \\ sqrtn}^2 \end{gathered}$$。（记住退火是让一个随机值去逼近最后的答案）

因为${x\_0}^2$的值更加接近于$$\begin{gathered} { \\ sqrtn}^2 \end{gathered}$$。因此$x\_0$值就更加接近于$sqrtn$

#完整代码和注释供大家调试:

# 通过模拟退火算出根号n的值
import random
import math

n = int(input())  # n代表我们最后函数要逼近的值
dT = 0.993  # 代表系数delta T，变化率，这里需要速度稍微慢一点，写0.995 或者 0.997都可以，但是越靠近1速度就越慢
eps = 1e-14  # 10的-14次方已经非常小了，写这个大概率没问题


def func(x):
    return math.fabs(x * x - n)


def SA():
    x = 0  # 首先随机生成一个点x0，这里我用0代替。
    T = 20000  # 代表开始的温度，初始温度主要是看题目开始需要跳转的幅度。
    f = func(x)  # 算出x平方和n的差距f(x0)
    while T > eps:
        # 这里x0既可以变小，也可以变大，所以我们正负都要进行一个跳转,算出变换后的点dx
        dx = x + (random.randint(0, 32767) * 2 - 32767) * T
        # 但是请注意，dx很明显要保证 >= 0才行，因为算术平方根的定义域是>=0，因此小于0就重新随机
        while dx < 0:
            dx = x + (random.randint(0, 32767) * 2 - 32767) * T
        df = func(dx)  # 算出变换后的点dx的平方和n的差距，f(dx)
        # 这里就是关键的地方了，很明显我们需要算出来的function值越小，自变量x更加接近那个根号值。所以如果新来的值df 比 f更小，我们百分百接受
        if df < f:
            x = dx
            f = df
        # 否则我们概率接受，这里的需要写 f - df了，因为这样才是负值。负值说明我们并不是贪心接受的，他是不太好的值。
        elif math.exp((f - df) / T) * 32767 > random.randint(0, 32767):
            x = dx
            f = df
        # 温度下降一下
        T *= dT
    print("%.8f"% x)


SA()

假设输入为2，则输出为：

参考：
B站视频：大学生速通模拟退火算法_哔哩哔哩_bilibili
博客：速通模拟退火 - Don’t Cry！