跳转至

模擬退火

引入

模擬退火是一種隨機化算法。當一個問題的方案數量極大（甚至是無窮的）而且不是一個單峯函數時，我們常使用模擬退火求解。

解釋

根據 爬山算法 的過程，我們發現：對於一個當前最優解附近的非最優解，爬山算法直接捨去了這個解。而很多情況下，我們需要去接受這個非最優解從而跳出這個局部最優解，即為模擬退火算法。

什麼是退火？（選自 百度百科）

退火是一種金屬熱處理工藝，指的是將金屬緩慢加熱到一定温度，保持足夠時間，然後以適宜速度冷卻。目的是降低硬度，改善切削加工性；消除殘餘應力，穩定尺寸，減少變形與裂紋傾向；細化晶粒，調整組織，消除組織缺陷。準確的説，退火是一種對材料的熱處理工藝，包括金屬材料、非金屬材料。而且新材料的退火目的也與傳統金屬退火存在異同。

由於退火的規律引入了更多隨機因素，那麼我們得到最優解的概率會大大增加。於是我們可以去模擬這個過程，將目標函數作為能量函數。

過程

先用一句話概括：如果新狀態的解更優則修改答案，否則以一定概率接受新狀態。

我們定義當前温度為 \(T\)，新狀態 \(S'\) 與已知狀態 \(S\)（新狀態由已知狀態通過隨機的方式得到）之間的能量（值）差為 \(\Delta E\)（\(\Delta E\geqslant 0\)），則發生狀態轉移（修改最優解）的概率為

\[ P(\Delta E)= \begin{cases} 1, & S' \text{ is better than } S,\\ \mathrm{e}^\frac{-\Delta E}{T}, & \text{otherwise}. \end{cases} \]

注意：我們有時為了使得到的解更有質量，會在模擬退火結束後，以當前温度在得到的解附近多次隨機狀態，嘗試得到更優的解（其過程與模擬退火相似）。

如何退火（降温）？

模擬退火時我們有三個參數：初始温度 \(T_0\)，降温係數 \(d\)，終止温度 \(T_k\)。其中 \(T_0\) 是一個比較大的數，\(d\) 是一個非常接近 \(1\) 但是小於 \(1\) 的數，\(T_k\) 是一個接近 \(0\) 的正數。

首先讓温度 \(T=T_0\)，然後按照上述步驟進行一次轉移嘗試，再讓 \(T=d\cdot T\)。當 \(T<T_k\) 時模擬退火過程結束，當前最優解即為最終的最優解。

注意為了使得解更為精確，我們通常不直接取當前解作為答案，而是在退火過程中維護遇到的所有解的最優值。

引用一張 Simulated annealing - Wikipedia 的圖片（隨着温度的降低，跳躍越來越不隨機，最優解也越來越穩定）。

實現

此處代碼以 「BZOJ 3680」吊打 XXX（求 \(n\) 個點的帶權類費馬點）為例。

 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
46
47
48
49
50
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>

const int N = 10005;
int n, x[N], y[N], w[N];
double ansx, ansy, dis;

double Rand() { return (double)rand() / RAND_MAX; }

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 simulateAnneal() {
  double t = 100000;
  double nowx = ansx, nowy = ansy;
  while (t > 0.001) {
    double nxtx = nowx + t * (Rand() * 2 - 1);
    double nxty = nowy + t * (Rand() * 2 - 1);
    double delta = calc(nxtx, nxty) - calc(nowx, nowy);
    if (exp(-delta / t) > Rand()) nowx = nxtx, nowy = nxty;
    t *= 0.97;
  }
  for (int i = 1; i <= 1000; ++i) {
    double nxtx = ansx + t * (Rand() * 2 - 1);
    double nxty = ansy + t * (Rand() * 2 - 1);
    calc(nxtx, nxty);
  }
}

int main() {
  srand(0);  // 注意，在实际使用中，不应使用固定的随机种子。
  scanf("%d", &n);
  for (int i = 1; i <= n; ++i) {
    scanf("%d%d%d", &x[i], &y[i], &w[i]);
    ansx += x[i], ansy += y[i];
  }
  ansx /= n, ansy /= n, dis = calc(ansx, ansy);
  simulateAnneal();
  printf("%.3lf %.3lf\n", ansx, ansy);
  return 0;
}

一些技巧

分塊模擬退火

有時函數的峯很多，模擬退火難以跑出最優解。

此時可以把整個值域分成幾段，每段跑一遍模擬退火，然後再取最優解。

卡時

有一個 clock() 函數，返回程序運行時間。

可以把主程序中的 simulateAnneal(); 換成 while ((double)clock()/CLOCKS_PER_SEC < MAX_TIME) simulateAnneal();。這樣子就會一直跑模擬退火，直到用時即將超過時間限制。

這裏的 MAX_TIME 是一個自定義的略小於時限的數（單位：秒）。

習題

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用