速通模拟退火(Simulated Annealing)
作者:郭程(judgemeido@163.com)
本教程为原创教程,转载请标明出处,禁止商业用途。
本教程属于快速入门(Quick-Start)系列,不对具体退火细节深究,只论算法竞赛题目研究。
1.是什么?为什么?怎么做?
Q1:模拟退火是什么算法?
模拟退火是模拟物理上退火方法,通过N次迭代(退火),逼近函数的上的一个最值(最大或者最小值)。
比如逼近这个函数的最大值C点:
Q2:模拟退火为什么可行?
讨论这个问题需要理解一下物理原型是怎么样的,也就是原来是怎么“退火”的:
模拟退火算法的思想借鉴于固体的退火原理,当固体的温度很高的时候,内能比较大,固体的内部粒子处于快速无序运动,当温度慢慢降低的过程中,固体的内能减小,粒子的慢慢趋于有序,最终,当固体处于常温时,内能达到最小,此时,粒子最为稳定。
注意标粗字体:
- 温度高->运动速度快(温度低->运动速度慢)
- 温度是缓慢(想象成特别慢的那种)降低的
- 温度基本不再变化后,趋于有序(最后内能达到最小,也就是接近最优)
我们通过模拟这个操作,使得我们需要的答案“趋于有序”,也就是最靠近需要的值(最值)。
Q3:怎么做?
大方向
首先,理解一下大方向:
模拟退火就是一种循环算法。
- 我们先设定一个初始的温度$T$(这个温度会比较高,比如2000)
- 每次循环都退火一次。(具体怎么操作后面详解)
- 然后降低$T$的温度,我们通过让$T$和一个“降温系数”$\Delta T$(一个接近1的小数,比如$0.99$)相乘,达到慢慢降低温度的效果,直到接近于0(我们用$eps$来代表一个接近0的数(比如0.00001),只要$T>eps$就可以退出循环了)
所以总的来说,用伪代码表示退火的流程是这样的:
double T = 2000; //代表开始的温度
double dT = 0.99; //代表系数delta T
double eps = 1e-14; //相当于0.0000000000000001
while(T > eps) {
//--------------
//这里是每一次退火的操作
//--------------
T = T * dT; //温度每次下降一点点, T * 0.99
}
退火详解
我们要求解的答案无非是两个:自变量$x$和对应函数的最大值$f(x)$
那么模拟退火的是怎么做到的呢?【下面车速很快,请系好安全带!】
我们先随机找一点$x_0$ ,不论是哪个点都可以,随机!(不超过定义域就行)。
这个点作为我们的初始值(相当于物体里的一个粒子)
- 再找到一点$f(x_0)$,来代表$x_0$所对应的函数值
现在正式开始退火!
刚才我们说了$x_0$相当于是一个粒子,所以我们会进行一个无序运动,也就是向左或者向右随机移动
是的,是随机移动,可能向左,也可能向右,但是请记住一个关键点:移动的幅度和当前的温度$T$有关。
温度$T$越大,移动的幅度越大。温度$T$越小,移动的幅度就越小。这是在模拟粒子无序运动的状态。
- 接受(Accept)更"好"的状态
假设我们移动到了$x_1$处,那么这个点对应的$f(x_1)$很明显答案是优于(大于)当前的$f(x_0)$的
因此我们将答案进行更新。也就是将初始值进行替换:$x_0=x_1,f(x_0)=f(x_1)$。这是一种贪心的思想。
以一定概率接受(Accept)更差的状态
这是退火最精彩的部分。
为什么我们要接受一个更加差的状态呢?因为可能在一个较差的状态旁边会出现一个更加高的山峰
如果我们鼠目寸光,只盯着右半区,很容易随着温度的下降、左右跳转幅度的减小而迷失自己,最后困死在小山丘中。
而我们如果找到了左边山峰的低点,以一定的概率接受了它(概率大小和温度以及当前的值的关键程度有关),会在跳转幅度减少之前,尽可能找到最优点。
那么我们以多少的概率去接受它呢?我们用一个公式表示(这个公式我们只需记住,这是科学家推导出来的结论):
$$ \Huge e^{\frac{\Delta f}{kT}} $$
别慌!很简单!我们来理解一下这里面的变量:
$e$是自然对数,约等于2.71。我们可以把右上角这一坨值$\frac{\Delta f}{kT}$看成一个整体$x$:
$e^x$的图形画出来是这样的:
因为我们想要函数$e^x$来代表一个概率值,所以我们只需要关注x为负数的部分即可:
负数部分的值域是在$(0,1)$开区间内,x越小,越接近0,越大越靠近1。
因为在0到1之间,所以这个值相当于是概率了。比如$e^x=0.97$,那么我们接受的概率就是$97\%$
而正数部分的值域会大于1,也就是说概率会超过$100\%$,所以会一定选(其实是上一种找到更优的情况)
- $kT$
$k$其实是个物理学常数,我们在代码中不会用到。
$T$很简单,就是当前的温度。所以实际上这个分母就是$T$,$k$当做$1$使用。
- $\Delta f$
我们着重讲一下什么是$\Delta f$。
其实从前面的函数$e^x$中可以发现,$\Delta f$必须是个负数!
我们想要函数$e^x$来代表一个概率值,一定要让它的值域属于$(0,1)$,所以$\frac{\Delta f}{kT}$必须是个负数。但是$kT$在我们的模拟中一定是正数,那么$\Delta f$必须是个负数!
其实$\Delta f$就是当前解的函数值与目标解函数值之差,$\Delta f=-\left|f\left(x_0\right)-f\left(x_1\right)\right|$,并且一定是个负数。这个需要具体问题具体分析。
比如现在我们求一个函数的最大值,那么如果$f(x_0) < f(x_1)$了,那就说明结果变好了,我们肯定选择它(见第4点)
如果$f(x_0) > f(x_1)$,那就说明结果变差了,我们需要概率选择它,因此$\Delta f=-(f(x_0) - f(x_1))$
所以总结一下就是:
随机后的函数值如果结果更好,我们一定选择它(即$x_0=x_1,f(x_0)=f(x_1)$)
随机后的函数值如果结果更差,我们以$\LARGE e^{\frac{\Delta f}{kT}}$的概率接受它
伪代码流程
注:对代码中的函数作出解释:
①对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)代表$e^x$
③关于exp((df - f) / T) * RAND_MAX > rand()
- 目的是要概率接受,但是$e^x$是个准确值,所以从理论上我们可以生成一个(0,1)的随机数,如果$e^x$比(0,1)这个随机数要大,那么我们就接受。
- 但是由于rand()只能产生[0,32767]内的随机整数,化成小数太过麻烦。所以我们可以把左边乘以RAND_MAX(也就是把概率同时扩大32767倍),效果就等同于$e^x$比(0,1)了。
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 = (2*rand() - RAND_MAX) * T;
//让x落在定义域内,如果没在里面,就重新随机。题目有要求需要写,否则不用写
// ================
while(x > ? || x < ? ...) {
double dx = (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;
2.通过模拟退火算出 根号n的值
思路是这样的:我们试图通过退火找出一个值$x_0$,使得${x_0}^2$的值更加接近于${\sqrt{n}}^2$。(记住退火是让一个随机值去逼近最后的答案)
因为${x_0}^2$的值更加接近于${\sqrt{n}}^2$。因此$x_0$值就更加接近于$\sqrt{n}$
- 所以我们需要一个函数$f(x)$,算出${x_0}^2$和${\sqrt{n}}^2$的接近程度,那么毋庸置疑,我们需要算出他们绝对值的差。
$$ f(x)=\left |x^2 - n \right | $$
也就是说我们的函数表达式就有了
//n代表我们最后函数要逼近的值
double n;
//x表示我们随机产生的那个数的平方和n的靠近程度
double func(double x) {
return fabs(x * x - n);
}
写出退火函数SA()
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()) { //注意更新所有变量 f = df; x = dx; } //温度下降一下 T *= dT; } printf("%.8lf",x); }
最后贴上完整代码和注释供大家调试。
#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()) {
//注意更新所有变量
f = df; x = dx;
}
//温度下降一下
T *= dT;
}
printf("%.8lf",x);
}
int main()
{
cin >> n;
SA();
}
3.例题POJ - 2420
题目描述
给出平面上N(N<=100)个点,你需要找到一个这样的点,使得这个点到N个点的距离之和尽可能小。输出这个最小的距离和(四舍五入到最近的整数)。
Input输入
第一行N,接下来N行每行两个整数,表示N个点
Output输出
一行一个正整数,最小的距离和。
Sample Input样例输入
4
0 0
0 10000
10000 10000
10000 0
Sample Output样例输出
28284
题目解析
学过平面几何的一定知道,假定两点的坐标为$a(x_1,y_1),b(x_2,y_2)$,那么他们之间的距离$d$可以用勾股定理计算:
$$ d = \sqrt{(x_1-x_2)^2-(y_1-y_2)^2} $$
所以我们的函数$func()$就是:求出一个随机点 $A=(x_0,y_0)$到$N$个点的距离之和$D$,就是把这个点A和所有N个点的距离相加,即:
$$ \large D=\sum_{i=1}^{N} {\sqrt{(x_0-x_i)^2-(y_0-y_i)^2}} $$
我们将N个点用结构体存起来
const int N = 1e5 + 5;
struct Point {
double x, y;
}e[N];
//main函数中输入N个点
int main()
{
cin >> n;
for(int i = 1; i <= n; i++)
scanf("%lf%lf", &e[i].x, &e[i].y);
}
那么我们最关键的函数(求出一个随机点$A=(x_0,y_0)$到$N$个点的距离之和$D$)就是:
//输入我们随机生成的点A 坐标 x0, y0
double func(double x, double y) {
double ans = 0;
for(int i = 1; i <= n; i++) {
double xx = e[i].x, yy = e[i].y; // xx 就是 xi ,yy 就是 yi
double p = fabs(e[i].x - x); // x0 - xi
double q = fabs(e[i].y - y); // y0 - yi
ans += (sqrt((p * p + q * q))); // ans加上到随机点A到点i的距离
}
return ans;
}
然后套用上面的公式,我们写出SA()模拟退火
double T = 2000; //初始温度
double dT = 0.993; //变化率,这里需要速度稍微慢一点,写0.995 或者 0.997都可以,但是越靠近1速度就越慢
const double eps = 1e-14; //10的-14次方已经非常小了,写这个大概率没问题
void SA() {
//随机生成点A(x0,y0),这里我直接赋值为了0。当然也可以写rand(),差别都不大。
double x = 0, y = 0;
//生成的点A(x0,y0)对应的函数值f(x0, y0)
double f = func(x, y);
while(T > eps) {
//这里对x和y同时进行一个偏移,很明显在一个平面中,上下左右都可以移动,所以我们选择了rand() * 2 - RAND_MAX
//对于每个点,我们的偏移幅度都要 * T,温度越高,偏移量越大
double dx = x + (rand() * 2 - RAND_MAX) * T;
double dy = y + (rand() * 2 - RAND_MAX) * T;
//算出偏移后的点B(dx, dy)对应的函数值f(dx, dy)
double df = func(dx, dy);
//这里就是关键的地方了,很明显我们需要算出来的function值越小,就更加优秀
//所以如果新来的值df 比 f更小,我们百分百接受
if(df < f) {
//注意更新所有变量
f = df; y = dy; x = dx;
}
//否则我们概率接受,这里的需要写 f - df了,因为这样才是负值
else if(exp((f - df) / T) * RAND_MAX > rand()) {
//注意更新所有变量
f = df; y = dy; x = dx;
}
//降低温度
T = T * dT;
}
//输出答案
printf("%.f\n", f);
}
最后贴上完全版代码
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <ctime>
using namespace std;
const int N = 1e5 + 5;
double dT = 0.993;
double T = 2000;
const double eps = 1e-14;
struct vega {
double x, y;
}e[N];
int n;
//输入我们随机生成的点A 坐标 x0, y0
double func(double x, double y) {
double ans = 0;
for(int i = 1; i <= n; i++) {
double xx = e[i].x, yy = e[i].y; // xx 就是 xi ,yy 就是 yi
double p = fabs(e[i].x - x); // x0 - xi
double q = fabs(e[i].y - y); // y0 - yi
ans += (sqrt((p * p + q * q))); // ans加上到随机点A到点i的距离
}
return ans;
}
void SA() {
//随机生成点A(x0,y0),这里我直接赋值为了0。当然也可以写rand(),差别都不大。
double x = 0, y = 0;
//生成的点A(x0,y0)对应的函数值f(x0, y0)
double f = func(x, y);
while(T > eps) {
//这里对x和y同时进行一个偏移,很明显在一个平面中,上下左右都可以移动,所以我们选择了rand() * 2 - RAND_MAX
//对于每个点,我们的偏移幅度都要 * T,温度越高,偏移量越大
double dx = x + (rand() * 2 - RAND_MAX) * T;
double dy = y + (rand() * 2 - RAND_MAX) * T;
//算出偏移后的点B(dx, dy)对应的函数值f(dx, dy)
double df = func(dx, dy);
//这里就是关键的地方了,很明显我们需要算出来的function值越小,就更加优秀
//所以如果新来的值df 比 f更小,我们百分百接受
if(df < f) {
//注意更新所有变量
f = df; y = dy; x = dx;
}
//否则我们概率接受,这里的需要写 f - df了,因为这样才是负值
else if(exp((f - df) / T) * RAND_MAX > rand()) {
//注意更新所有变量
f = df; y = dy; x = dx;
}
//降低温度
T = T * dT;
}
//输出答案
printf("%.f\n", f);
}
int main()
{
cin >> n;
for(int i = 1; i <= n; i++) {
scanf("%lf%lf", &e[i].x, &e[i].y);
}
SA();
}
求两个坐标之间的距离好像写错了,应该是sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))