blog1: Click here, blog2: Click here, video: Click here

模拟退火算法(SA)来源于固体退火原理，是一种基于概率的算法。将固体加温至充分高的温度，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，分子和原子越不稳定。而徐徐冷却时粒子渐趋有序，能量减少，原子越稳定。在冷却（降温）过程中，固体在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

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

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

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

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

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

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

注意标粗字体：

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

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

Q3:怎么做？

大方向

首先，理解一下大方向：

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

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

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

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

退火详解

我们要求解的答案无非是两个：自变量𝑥

和对应函数的最大值 𝑓(𝑥)

那么模拟退火的是怎么做到的呢？

我们先随机找一点𝑥0 ，不论是哪个点都可以，随机！（不超过定义域就行）。这个点作为我们的初始值（相当于物体里的一个粒子），再计算𝑓(𝑥0), 来代表𝑥0所对应的函数值。

现在正式开始退火！刚才我们说了𝑥0，相当于是一个粒子，所以我们会进行一个无序运动，也就是向左或者向右随机移动。是的，是随机移动，可能向左，也可能向右，但是请记住一个关键点：移动的幅度和当前的温度 𝑇 有关。温度T越大，移动的幅度越大。温度T越小，移动的幅度就小。这是在模拟粒子无序运动的状态。接受(Accept)更”好”的状态。

假设我们移动到了𝑥1处，那么这个点对应的𝑓(𝑥1)很明显答案是优于（大于）当前的𝑓(𝑥0)的。

因此我们将答案进行更新。也就是将初始值进行替换：𝑥0=𝑥1,𝑓(𝑥0)=𝑓(𝑥1)。这是一种贪心的思想。以一定概率接受(Accept)更差的状态这是退火最精彩的部分。为什么我们要接受一个更加差的状态呢？因为可能在一个较差的状态旁边会出现一个更加高的山峰。

如果我们鼠目寸光，只盯着右半区，很容易随着温度的下降、左右跳转幅度的减小而迷失自己，最后困死在小山丘中。

而我们如果找到了左边山峰的低点，以一定的概率接受了它（概率大小和温度以及当前的值的关键程度有关），会在跳转幅度减少之前，尽可能找到最优点。

那么我们以多少的概率去接受它呢？我们用一个公式表示（这个公式我们只需记住，这是科学家推导出来的结论）：𝑒Δ𝑓𝑘𝑇。

k为常数，这里可以看作1
T为温度

Δ𝑓就是当前解的函数值与目标解函数值之差，Δ𝑓=−

𝑓(𝑥0)−𝑓(𝑥1)

，并且一定是个负数。这个需要具体问题具体分析。

比如现在我们求一个函数的最大值，那么如果𝑓(𝑥0)<𝑓(𝑥1)了，那就说明结果变好了，我们肯定选择它（见第4点）。

如果𝑓(𝑥0)>𝑓(𝑥1)，那就说明结果变差了，我们需要概率选择它，因此Δ𝑓=−(𝑓(𝑥0)−𝑓(𝑥1))。

所以总结一下就是：

随机后的函数值如果结果更好，我们一定选择它(即𝑥0=𝑥1,𝑓(𝑥0)=𝑓(𝑥1))
随机后的函数值如果结果更差，我们以𝑒Δ𝑓𝑘𝑇的概率接受它

当𝑘𝑇越大时(温度越高)，由于Δ𝑓是一个负数，所以算出来的值会越大。比如-1/1000 等于-0.001，很明显-0.001>-1，由于𝑒^𝑥是个单调递增函数，所以温度越高时，接受的概率就越大。Δ𝑓越小，说明答案越糟糕，那么接受的概率就越小。

伪代码流程

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

①对rand()函数

rand()函数可以默认拿到[0,32767]内的随机整数

RAND_MAX = 32767，可以看作常量。本质是宏定义: #define RAND_MAX 32767

rand() * 2 的范围是[0,32767 * 2]

rand() * 2 - RAND_MAX 的范围是[-32767, 32767]

②对exp()函数

exp(x)代表𝑒^𝑥

关于exp((df - f) / T) * RAND_MAX > rand()，目的是要概率接受，但是𝑒^𝑥是个准确值，所以从理论上我们可以生成一个（0,1）的随机数，如果𝑒^𝑥比(0,1)这个随机数要大，那么我们就接受。

但是由于rand()只能产生[0,32767]内的随机整数，化成小数太过麻烦。所以我们可以把左边乘以RAND_MAX（也就是把概率同时扩大32767倍），效果就等同于𝑒^𝑥

C++ Code

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

//用自变量计算函数值,这里可能存在多个自变量对应一个函数值的情况，比如f(x,y)
double func(int x, ... ) {
    //这里是对函数值进行计算
    double ans = .......
    return ans;
}
//原始值
double x = rand(); //x0取随机值
double f = func(x,...); //通过自变量算出f(x0)的值
while(T > eps) {
    //--------------
    //这里是每一次退火的操作
    
    //x1可以左右随机移动，幅度和温度T正相关，所以*T
    //注意这里移动可以左右移动，但是也可以单向移动
    //关于rand()详细见开头注的①
    double dx = x + (2*rand() - RAND_MAX) * T; 
    
    //让x落在定义域内，如果没在里面，就重新随机。题目有要求需要写，否则不用写
    // ================
    while(x > ? || x < ? ...) {
        double dx = x + (2*rand() - RAND_MAX) * T; 
    }
    // ================
    
    //求出f(x1)的值
    double df = func(dx);
    //这里需要具体问题具体分析，我们要接受更加优秀的情况。可能是df < f(比如求最小值）
    if(f < df) {
        f = df; x = dx;  [...,y = dy;] // 接受，替换值，如果多个自变量，那么都替换
    }
    //否则概率接受，注意这里df-f也要具体问题具体分析。
    //详细见开头注的②③
    else if(exp((df - f) / T) * RAND_MAX > rand()) {
        f = df; x = dx;  [...y = dy;] // 接受，替换值，如果多个自变量，那么都替换
    }
	//--------------
    T = T * dT; //温度每次下降一点点， T * 0.99
}
//最后输出靠近最优的自变量x值，和函数值f(x)
cout << x << " " << f << endl;

Python Code

        
      
import math
from random import random
import matplotlib.pyplot as plt

def func(x, y):                  #函数优化问题
    res= 4*x**2-2.1*x**4+x**6/3+x*y-4*y**2+4*y**4
    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)
            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:       
           
            #内循环迭代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()
        print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")

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()

模拟退火 —— Simulated Annealing