模拟退火算法

本文最后更新于：2024年8月27日 上午

模拟退火算法（Simulated Annealing）

1.是什么？怎么来的？

模拟退火是模拟物理上退火方法，通过N次迭代（退火），逼近函数的上的一个最值（最大或者最小值）。尽管可能不是最优解，但得到的是一个可行解。

很多优秀算法的设计灵感都源于自然界和社会生活当中，比如蚁群算法，遗传算法等等，模拟退火算法也不例外，在物理退火当中，先将物体加温，最后让物体自然冷却，缓缓降温，冷却后可达到最低能量状态，较为柔韧，变得更容易弯折，适用于制作各种饰品。哈哈哈大自然都知道慢工出细活😁！

物理退火分为三个过程

  • 加温过程——增强粒子的热运动，消除系统原先可能存在的非均匀态；
  • 等温过程——对于与环境换热而温度不变的封闭系统，系统状态的自发变化总是朝自由能减少的方向进行，当自由能达到最小时，系统达到平衡态；
  • 冷却过程——使粒子热运动减弱并渐趋有序，系统能量逐渐下降，从而得到低能的晶体结构。

最早的想法是由N.Metropolis 等人于1953 年所提出，在当时并沒有受到重视。直到1983 年由Kirkpatrick et al. 提出蒙特卡罗模拟（MonteCarlo Simulated）概念的随机搜寻技巧，利用此方法来求解的组合优化问题时，才使此算法受到重视。

2.怎么做？

大方向

首先，理解一下大方向：

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

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

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

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

退火详解

  1. 第一步，我们需要在找到最值的那个函数里随机找到一点X0，作为我们的初始值。

  2. 第二步，其实我们可以把寻找最值的过程比喻成在山峦如聚中寻找最高的山峰，我们在上课时，老师也是这么说的，那么问题来了，下一步我们该往哪边走呢？假设只有一个二维的平面，往左还是往右呢？答案应该是随机的吧，因为你也不知道哪边的山会更高。那到底是怎么一个随机法呢？上课的时候这个问题难倒了不少人，同学们各抒己见。结论是：移动的幅度和当前的温度T有关。温度T越大，移动的幅度越大。温度T越小，移动的幅度就越小。这是在模拟粒子无序运动的状态。关于这一点的公式解释，我们会在后面说到。

  3. 第三步，我们需要接受更好的状态，假设我们移动到了X1处，此时的山峰比刚才的山峰要高，那我们就接受这个值，并将此时的值设为目前的最优解。

  4. 第四步，如果我们不小心到达了一个更低的山峰，那我们就一定抛弃这个结果吗？不一定，因为可能低谷的再往同一个方向走，就会出现更高的山峰。如果我们鼠目寸光，只盯着当前的局部最优解，很有可能随着温度的下降，移动幅度减小而陷入小山丘中。这也是人生的大智慧啊！而我们如果到了那个更低的山峰，以一定的概率接受了它（概率大小和温度以及当前的值的关键程度有关），会在跳转幅度减少之前，尽可能找到最优点，也就是那更高的山峰。

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

  5. 第五步，$\Delta f$必须是一个负数，因为概率是属于0到1的，e的0次方等于1，只要是次数大于0，这个值都会大于1，而$k$和$T$都是大于0的。所以只有当次数小于0时，这个值域才和概率相符。

所以总结一下就是：

  • 随机后的函数值如果结果更好，我们一定选择它
  • 随机后的函数值如果结果更差，我们以一定的概率接受它
  • 在前期，由于T的基础数据十分大，分子十分活跃，所以即使我们取到了一个不是那么好的数据，我们接受它的概率也比较大，因为T是分母，k是常数，T变小时，$ \huge e^{\frac{\Delta f }{kT}} $就趋近于0，而且是始终小于1的

3.例题：求$ \sqrt {n} $

  1. rand()函数可以默认取到[0, 32767]内的随机整数
  2. RAND_MAX = 32767,是常量
  3. rand()*2-RAND_MAX的范围是[-32767, 32767]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#include <bits/stdc++.h>
using namespace std;
//n代表我们最后函数要逼近的值
double n;
//x表示我们随机产生的那个数的平方和n的靠近程度
double func(double x) {
    return fabs(x * x - n);
}
double T = 20000; //初始温度，初始温度主要是看题目开始需要跳转的幅度。
double dT = 0.993; //变化率，这里需要速度稍微慢一点，写0.995 或者 0.997都可以，但是越靠近1速度就越慢 
const double eps = 1e-14; //10的-14次方已经非常小了，写这个大概率没问题
void SA() {
    //首先随机生成一个点x0，这里我用0代替。
    double x = 0;
    //算出x平方和n的差距f(x0)
    double f = func(x);
    while(T > eps) {
        //这里x0既可以变小，也可以变大，所以我们正负都要进行一个跳转,算出变换后的点dx
        double dx = x + (2 * rand() - RAND_MAX) * T;
        //但是请注意，dx很明显要保证 >= 0才行，因为算术平方根的定义域是>=0，因此小于0就重新随机
        while(dx < 0) dx = x + (2 * rand() - RAND_MAX) * T;
        //算出变换后的点dx的平方和n的差距，f(dx)
        double df = func(dx);
        //这里就是关键的地方了，很明显我们需要算出来的function值越小，自变量x更加接近那个根号值。
        //所以如果新来的值df 比 f更小，我们百分百接受
        if(df < f) {
            //注意更新所有变量
            f = df; x = dx;
        }
        //否则我们概率接受，这里的需要写 f - df了，因为这样才是负值。负值说明我们并不是贪心接受的，他是不太好的值。
        else if(exp((f - df)/T) * RAND_MAX > rand()) {
            //把概率扩大32767倍，使得两边函数值域相同
            //注意更新所有变量
            f = df; x = dx;
        }
        //温度下降一下
        T *= dT;
    }
    printf("%.8lf",x);
}
int main() 
{
    cin >> n;
    SA();
}

参考文章：速通模拟退火 - Don’t Cry！ (tocry.cn)

笔记
算法

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！