模拟退火学习笔记

· · 个人记录

退火定义

退火是一种金属热处理工艺，指的是将金属缓慢加热到一定温度，保持足够时间，然后以适宜速度冷却。目的是降低硬度，改善切削加工性；消除残余应力，稳定尺寸，减少变形与裂纹倾向；细化晶粒，调整组织，消除组织缺陷。准确的说，退火是一种对材料的热处理工艺，包括金属材料、非金属材料。而且新材料的退火目的也与传统金属退火存在异同。

模拟退火定义

模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法， 用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。

模拟退火原理

模拟退火的原理也和金属退火的原理近似：将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。

——百度百科

算法描述

如果新状态的解更优则修改答案，否则以一定概率接受新状态

亦即：

模拟退火时有三个参数，分别是初始温度 T_{0} 、降温系数 \Delta 、终止温度 T_{k}

我们先让温度 $T=T_{0}$ ，然后每次降温时 $T=T\cdot\Delta$ ，直到 $T\leq T_{k}$ 为止。 设这个新的解与最优解的差为 $\Delta E$，温度为 $T$ ，那么这个概率为 $e^{\frac{-\Delta E}{T}}

大致过程如下：

—— Wiki - Simulated annealing

可以看出，随着温度的降低，解逐渐稳定下来，并逐渐集中在最优解附近。

步骤

  1. 选定一个初始状态（比如选定所有点坐标的平均数）& 初始温度 T_{0}

  2. 当温度大于一个边界值时：

    {
　　随机变化坐标，变化幅度为 T 。
　　计算新解与当前解的差 ΔE。
　　如果新解比当前解优(ΔE > 0)，就用新解替换当前解。
　　否则以 e ^ (ΔE / T) 的概率用新解替换当前解。
　　温度乘上一个小于 1 的系数，即"降温"。
}

如何生成新解？

  1. 坐标系内：随机生成一个点，或者生成一个向量。
  2. 序列问题： random_shuffle 或者随机交换两个数。
  3. 网格问题：可以看做二维序列，每次交换两个格子即可。

注意事项

1.温度T的初始值设置问题
初始温度高，搜索到全局最优解的可能性大，但因此要花费大量的计算时间；
反之，则可节约计算时间，但全局搜索性能可能受到影响。

2.退火速度问题
模拟退火算法的全局搜索性能也与退火速度密切相关，
同一温度下的“充分”搜索(退火)是相当必要的，但这需要计算时间。

3.温度管理问题 降温系数应为正的略小于1.00的常数。

例题

P1337 [JSOI2004]平衡点 / 吊打XXX

此题自然变化进行的方向都是使能量降低，因为能量较低的状态比较稳定。

因为物重一定，绳子越短，重物越低，势能越小，势能又与物重成正比，所以，只要使得 \sum_{i=1}^ndist_{i} \times weight_{i} 也就是总的重力势能最小，就可以使系统平衡

Code

#include <bits/stdc++.h>
using namespace std;

const int N = 10005;
int n, x[N], y[N], w[N];
double ansx, ansy, dis;//全局最优解坐标 & 全局最优解

double Rand() { return (double)rand() / RAND_MAX; } //RAND_MAX 是 <stdlib.h> 中伪随机数生成函数 rand 所能返回的最大数值

int read()
{
  int s=0,w=1;char g=getchar();
  while(g>'9'||g<'0'){if(g=='-')w=0;g=getchar();}
  while(g>='0'&&g<='9'){s=(s<<1)+(s<<3)+(g^48);g=getchar();}
  return w?s:-s;
}

double calc(double xx, double yy) { //计算整个系统的能量
  double res = 0;
  for (int i = 1; i <= n; ++i) {
    double dx = x[i] - xx, dy = y[i] - yy;
    res += sqrt(dx * dx + dy * dy) * w[i];
  }
  if (res < dis) dis = res, ansx = xx, ansy = yy;
  return res;
}

void SA() { //SA核心(代码与题面无关)
  double t = 2000; //参数1：初始温度
  double nowx = ansx, nowy = ansy;
  while (t > 1e-16) {//参数2：目标温度
    double X = nowx + t * (Rand() * 2 - 1);
    double Y = nowy + t * (Rand() * 2 - 1); //得出一个新的坐标
    double delta = calc(X, Y) - calc(nowx, nowy); //降温系数delta
    if (exp(-delta / t) > Rand()) nowx = X, nowy = Y; //接受新的答案
    t *= 0.998;//最重要参数 参数3：徐徐退火
  }
}

int main() {
  srand(20030624);//机房dalao生日
  n=read();
  for (int i = 1; i <= n; ++i) {
    x[i]=read(), y[i]=read(), w[i]=read();
    ansx += x[i], ansy += y[i];
  }
  ansx /= n, ansy /= n, dis = calc(ansx, ansy);
  while ((double)clock()/CLOCKS_PER_SEC < 0.6) SA();//今天你洗脸了没有
  printf("%.3lf %.3lf\n", ansx, ansy);
  return 0;
}

时间复杂度

O($玄学$)

卡时

有一个 clock() 函数，返回程序运行时间。

可以把主程序中的 simulateAnneal(); 换成 while ((double)clock()/CLOCKS_PER_SEC < MAX_TIME) simulateAnneal(); 。这样子就会一直跑模拟退火，直到用时即将超过时间限制。

这里的 MAX_TIME 是一个自定义的略小于时限的数。

*分块模拟退火

有时函数的峰很多，模拟退火难以跑出最优解。

此时可以把整个值域分成几段，每段跑一遍模拟退火，然后再取最优解。