焼きなまし法のコツ Ver. 1.3

TopCoder マラソンマッチ

焼きなまし法そのものの解説は少しだけにして、焼きなまし法の使い方のコツをメインに書こうと思います。

焼きなまし法の概要

最適化（=最大値か最小値を求める）の手法の１つです。以下は擬似コードです。

あらかじめ、計算する時間を決める。
以下のスコープ内のコードを時間経過するまで繰り返す。
{
    適当に、次の状態（近傍）を選ぶ。
    (1) 次の状態のスコアが、今の状態のスコアより良くなった場合は、常に、状態を次の状態にする（遷移）。
    (2) 次の状態のスコアが、今の状態のスコアより悪くなった場合でも、ある確率で、状態を次の状態にする。
        確率は以下のように決める。
       (2.A) 序盤(経過時間が0に近いとき)ほど高確率で、終盤ほど低確率。
       (2.B) スコアが悪くなった量が大きければ大きいほど、低確率。
}
最後に、今までで一番スコアのよかった状態を選ぶ。

次の状態が悪くてもそれを選ぶ可能性がある(2)の部分が、焼きなまし法の最大の特徴です。

「状態」「スコア」というのがピンとこないかもしれませんが、

  • もし、巡回セールスマン問題であれば、
    • 状態 →ルート
    • 次の状態（近傍）→元のルートをちょっとだけかえてみたルート
    • スコア →総距離
  • もし、y=f(x)の最大値を求める問題であれば、
    • 状態 →x
    • 次の状態（近傍）→ x' = x + a (aは小さめの値)
    • スコア →y

といったかんじです。

サンプルソースコード

AtCoderのコンテストIntroduction to Heuristics Contest - AtCoderで、スケジューリング問題が出題されました。 この問題は焼きなまし法が効果的です。C++でのサンプルを2つ用意しました。変数名は問題に合わせています。

サンプル1

単純な焼きなまし法。1165人中400位相当の解法です。焼きなまし法の本質は、// 焼きなまし法（ここから）の部分です。
AtCoderの提出コード


ソースコードを開く

#include <iostream>
#include <vector>
#include <algorithm>
#include <chrono>
#include <cassert>
#include <climits>
#include <cmath>

using namespace std;
using namespace chrono;

// 0以上UINT_MAX(0xffffffff)以下の整数をとる乱数 xorshift https://ja.wikipedia.org/wiki/Xorshift
static uint32_t randXor()
{
    static uint32_t x = 123456789;
    static uint32_t y = 362436069;
    static uint32_t z = 521288629;
    static uint32_t w = 88675123;
    uint32_t t;

    t = x ^ (x << 11);
    x = y; y = z; z = w;
    return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
}

// 0以上1未満の小数をとる乱数
static double rand01()
{
    return (randXor() + 0.5) * (1.0 / UINT_MAX);
}

const static int NUM_TYPES = 26; // コンテストの種類数

// 全体のスコア（今回の問題のスコア）
static int getFullScore(const int D, const vector <int>& C, const vector <int>& T, const vector<vector<int>>& S)
{
    int score = 0;

    vector <int> last(NUM_TYPES, 0);

    for (int d = 0; d < D; d++)
    {
        const int j = T[d];
        last[j] = d + 1;
        for (int t = 0; t < NUM_TYPES; t++)
        {
            score -= (d + 1 - last[t]) * C[t];
        }
        score += S[d][j];
    }

    return score;
}

int main()
{
    //  入力
    int D;
    cin >> D;
    vector<int> C(NUM_TYPES);
    for (int &c : C)
    {
        cin >> c;
    }

    vector<vector<int>> S(D, vector<int>(NUM_TYPES));
    for (auto &s : S)
        for (auto &x : s)
        {
            cin >> x;
        }

    // 初期解
    vector<int> T(D);   // T[d] d日目に開くコンテスト
    for (int d = 0; d < D; d++)
    {
        T[d] = d % NUM_TYPES;
    }

    // 焼きなまし法（ここから）
    auto startClock = system_clock::now();
    int score = getFullScore(D, C, T, S);   // 現在のスコア
    int bestScore = INT_MIN;                // 最高スコア
    vector <int> bestT;                     // 最高スコアを出したときのT

    const static double START_TEMP = 1500; // 開始時の温度
    const static double END_TEMP = 100; // 終了時の温度
    const static double END_TIME = 1.8; // 終了時間（秒）

    double time = 0.0;   // 経過時間（秒）
    do
    {
        const double progressRatio = time / END_TIME;   // 進捗。開始時が0.0、終了時が1.0
        const double temp = START_TEMP + (END_TEMP - START_TEMP) * progressRatio;

        // 近傍 : d日目のコンテストをtに変えてみる。d,tはランダムに選ぶ。
        const int d = randXor() % D;
        const int t = randXor() % NUM_TYPES;
        const int backupT = T[d];
        int deltaScore = 0; // スコアの差分 = 変更後のスコア - 変更前のスコア
        deltaScore -= getFullScore(D, C, T, S);
        T[d] = t;   // 変更
        deltaScore += getFullScore(D, C, T, S);
        const double probability = exp(deltaScore / temp); // 焼きなまし法の遷移確率

        if (probability > rand01())
        {
            // 変更を受け入れる。スコアを更新
            score += deltaScore;
            if (score > bestScore)
            {
                bestScore = score;
                bestT = T;
            }
        }
        else
        {
            // 変更を受け入れないので、元に戻す
            T[d] = backupT;
        }

        time = (duration_cast<microseconds>(system_clock::now() - startClock).count() * 1e-6);
    } while (time < END_TIME);
    // 焼きなまし法（ここまで）


    // 出力
    for (int t : bestT)
    {
        cout << t + 1 << endl;
    }
    cerr << "bestScore: " << max(0, bestScore + 1000000) << endl;

    return 0;
}


サンプル2

やや複雑ですが、かなり改善した焼きなまし法。1位相当の解法です。
AtCoderの提出コード


ソースコードを開く

#include <iostream>
#include <vector>
#include <algorithm>
#include <chrono>
#include <cassert>
#include <climits>
#include <cmath>

using namespace std;
using namespace chrono;

// 0以上UINT_MAX(0xffffffff)以下の整数をとる乱数 xorshift https://ja.wikipedia.org/wiki/Xorshift
static uint32_t randXor()
{
    static uint32_t x = 123456789;
    static uint32_t y = 362436069;
    static uint32_t z = 521288629;
    static uint32_t w = 88675123;
    uint32_t t;

    t = x ^ (x << 11);
    x = y; y = z; z = w;
    return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
}

// 0以上1未満の小数をとる乱数
static double rand01()
{
    return (randXor() + 0.5) * (1.0 / UINT_MAX);
}

const static int NUM_TYPES = 26; // コンテストの種類数

// 全体のスコア（今回の問題のスコア）
static int getFullScore(const int D, const vector <int>& C, const vector <int>& T, const vector<vector<int>>& S)
{
    int score = 0;

    vector <int> last(NUM_TYPES, 0);

    for (int d = 0; d < D; d++)
    {
        const int j = T[d];
        last[j] = d + 1;
        for (int t = 0; t < NUM_TYPES; t++)
        {
            score -= (d + 1 - last[t]) * C[t];
        }
        score += S[d][j];
    }

    return score;
}

// d日目のコンテストT[d]を開いたときのスコア増分
static int getContestScore(const int d, const int D, const vector <int>& C, const vector <int>& T, const vector<vector<int>>& S)
{
    const int t = T[d];
    int score = S[d][t];

    int prev = d - 1;   // d日目以前に開かれる同じコンテストの日
    while (prev >= 0 && T[prev] != t)
    {
        prev--;
    }

    int next = d + 1;   // d日目以降に開かれる同じコンテストの日
    while (next < D && T[next] != t)
    {
        next++;
    }

    // sum(x) = x*(x+1)/2 つまり1,2...,xまでの等差数列の和とすると
    //	score -= sum(next - prev -1) * (-C[t]);                       // ... (1)
    //	score += sum(day  - prev -1) + sum(next - day  -1) * (-C[t]); // ... (2)

    // p---------n コンテストをおく前
    //  123456789 // ... (1)

    // p---d-----n コンテストをおいた後
    //  123 12345 // ... (2) 

    score += (next - d) * (d - prev) * C[t];    // (1)と(2)を計算すると、簡単になる。

    return score;
}

int main()
{
    auto startClock = system_clock::now();

    //  入力
    int D;
    cin >> D;
    vector<int> C(NUM_TYPES);
    for (int &c : C)
    {
        cin >> c;
    }

    vector<vector<int>> S(D, vector<int>(NUM_TYPES));
    for (auto &s : S)
        for (auto &x : s)
        {
            cin >> x;
        }

    vector<int> T(D);   // T[d] d日目に開くコンテスト
    for (int d = 0; d < D; d++)
    {
        T[d] = d % NUM_TYPES;
    }

    int score = getFullScore(D, C, T, S);   // 現在のスコア
    int bestScore = INT_MIN;                // 最高スコア
    vector <int> bestT;                     // 最高スコアを出したときのT

    const static double START_TEMP = 1500; // 開始時の温度
    const static double END_TEMP   =  100; // 終了時の温度
    const static double END_TIME   =  1.8; // 終了時間（秒）

    double temp = START_TEMP;   // 現在の温度

    long long steps;    // 試行回数
    for (steps = 0; ; steps++)
    {
        if (steps % 10000 == 0)
        {
            const double time = duration_cast<microseconds>(system_clock::now() - startClock).count() * 1e-6;   // 経過時間（秒）
            if (time >= END_TIME)
            {
                break;
            }

            const double progressRatio = time / END_TIME;   // 進捗。開始時が0.0、終了時が1.0
            temp = START_TEMP + (END_TEMP - START_TEMP) * progressRatio;
        }

        if (rand01() > 0.5)
        {
            // 近傍1 : d日目のコンテストをtに変えてみる
            const int d = randXor() % D;
            const int t = randXor() % NUM_TYPES;
            const int backupT = T[d];
            int deltaScore = 0;

            // 変更前のd日目のコンテストのスコアを引き、変更して、変更後のコンテストのスコアtを足す。
            deltaScore -= getContestScore(d, D, C, T, S);
            T[d] = t;
            deltaScore += getContestScore(d, D, C, T, S);
            const double probability = exp(deltaScore / temp); // 焼きなまし法の遷移確率

            if (probability > rand01())
            {
                // 変更を受け入れる。スコアを更新
                score += deltaScore;
                if (score > bestScore)
                {
                    bestScore = score;
                    bestT = T;
                }
            }
            else
            {
                // 変更を受け入れないので、元に戻す
                T[d] = backupT;
            }
        }
        else
        {
            // 近傍2 : d0日目のコンテストとd1日目コンテストを交換してみる
            const int diff = randXor() % 10 + 1;    // d0とd1は10日差まで
            const int d0 = randXor() % (D - diff);
            const int d1 = d0 + diff;

            int deltaScore = 0;

            // 変更前のd0日目とd1日目のコンテストのスコアを引き、変更後のd0日目とd1日目のコンテストのスコアを足す
            deltaScore -= getContestScore(d0, D, C, T, S);
            deltaScore -= getContestScore(d1, D, C, T, S);
            swap(T[d0], T[d1]);
            deltaScore += getContestScore(d0, D, C, T, S);
            deltaScore += getContestScore(d1, D, C, T, S);

            const double probability = exp(deltaScore / temp); // 焼きなまし法の遷移確率
            if (probability > rand01())
            {
                // 変更を受け入れる。スコアを更新
                score += deltaScore;
                if (score > bestScore)
                {
                    bestScore = score;
                    bestT = T;
                }
            }
            else
            {
                // 変更を受け入れないので、元に戻す
                swap(T[d0], T[d1]);
            }
        }
    }

    // 出力
    for (int t : bestT)
    {
        cout << t + 1 << endl;
    }
    cerr << "steps: " << steps << endl;
    cerr << "bestScore: " << max(0, bestScore + 1000000) << endl;

    return 0;
}

焼きなまし法のコツ

★が多いほど、重要だと思っています。

1. 高速化

試行回数が多いほど、スコアが伸びる可能性があります。基本的にはボトルネックを調べて、そこを高速化しましょう。だいたい以下のところがボトルネックになるでしょう。

★★★★ プロファイルして、計算時間がかかっているボトルネックを探す

2021/03/06 追加しました。
焼きなまし法に限らず、マラソンマッチに限らず、プログラミング一般であらゆる高速化の前にやるべきことです。プロファイルの仕方はOSやツールにより違いますが、例えばVisual Studioなら「デバッグ」→「パフォーマンス プロファイラ」で、一発でできます。思わぬところがボトルネックになってることも多いので、絶対しましょう。
f:id:shindannin:20210306115133p:plain

★★★★ 高速化がどれだけ効果的か検討する。

重要です。ためしに10秒だけやってる計算を、あえて100秒とかに伸ばしてみて、結果がどうなるか見てみましょう。

  • 時間をかければかけるほどスコアが伸びる
    • 高速化しよう＆スコアが素早く伸びるように次の状態を決めよう。山登り法に変えるのも良い。
  • 時間をかけてもスコアが伸びない & スコアが高い
    • 最適解がでてる可能性が高い。計算を打ち切って、余った時間を他に使う
  • 時間をかけてもスコアが伸びない & スコアが低い
    • 次の状態の決め方、スコア関数の見直しを再検討する。焼きなまし法を使わない法がよいときもあるのも忘れずに。

こういった実験は自分のPCで行うわけなので、並列化などを行い計算リソースも存分に使ってよいです。

★★★★ 変更がある変数だけを保存して、戻す

先ほどのコードだと、

        const int backupT = T[d];

        // （略）    

        {
            // 変更を受け入れないので、元に戻す
            T[d] = backupT;
        }

の部分です。変更した、T[d]の部分だけをバックアップしておけば、元に戻せます。T全体を保存する必要はありません。

★★ 逆の手で、戻す

逆の手をやれば、前の状態に高速に戻せるのであれば、そうしましょう。例えば、上下左右に駒をうごかせるというパズル問題であれば、左に動かした後であれば、右に動かせば、戻せます。また、2つをswapするといった場合は、もう1回swapすれば戻せます。

★★★★ スコア関数の高速化

スコア関数の求め方は問題によって大きく異なるので、計算時間も異なってきます。ボトルネックになってたら高速化しましょう。
また、スコア関数の正確さを犠牲にしてでも、高速化したほうが良い場合もあります。１つのコツとして覚えて置きましょう。

★★★★ スコア差分計算時の打ち切り

2021/03/06 追加しました。
スコア差分計算の途中で、差分がすでに大きくマイナスになっていて、もうその後計算し続けていても絶対に遷移の受理確率がほぼ0っぽいのであれば、差分計算を打ち切って高速化することができます。焼きなまし法は遷移が成功するケースより遷移が失敗するケースが多く、失敗するケースをさっさと打ち切ることができれば高速化につながることが多いです。
また問題によっては、スコア差分の上界・下界が分かるケースもあり、これらが短時間で求められるのであれば、それもスコア差分計算の打ち切りに使用可能です。例えば、Google Hash Code 2021の練習問題では、10000ビットのORを取り、(ビットが1の数)の2乗がスコアになる問題がありました。10000ビットのORはそれなりに時間がかかりますが、

（ほんとは10000ビットある）

1110010101001   （1が7個）
0011000101000   （1が4個）  
↓ OR
1111010101001   （1が9個。つまり81点）

と実際にORをとらずとも、

1110010101001   （1が7個）
0011000101000   （1が4個）  
↓ OR
???????????????    （1が最小で7個、最大で11個、つまり49点～121点のあいだ）

というのが分かれば、打ち切り時の計算の参考になります。

★★★ 状態を変えずにスコア差分を計算する

2021/03/06 追加しました。
状態の更新をする部分がボトルネックとなっている場合は、
サンプル2の以下のような書き方よりも、

状態を変える
差分スコアを求める
遷移確率probabilityを求める
if (probability > rand01()) // 遷移条件を満たすか
{
  // 変更を受け入れる。スコアを更新
}
else
{
  // 変更を受け入れないので、状態を元に戻す
}

以下のような形のほうが良いです。なぜなら状態遷移しないときに状態更新の計算を省けるからです。また、スコア差分計算時の打ち切りを行う際も、打ち切った際に状態を戻す必要がないため、この形のほうが恩恵が大きいです。

状態を変えずに差分スコアを求める
遷移確率probabilityを求める
if (probability > rand01()) // 遷移条件を満たすか
{
  // 変更を受け入れる。スコアを更新。状態も更新
}

ただ、短時間で書けるほうは大抵サンプル2の形なので、まずはサンプル2の形で書いて、その後この形に持っていくのが良いと思います。

★★ C言語のrand()を使わない。

C言語のrand()は遅いので、高速なxor128を使いましょう。

unsigned int randxor()
{
    static unsigned int x=123456789,y=362436069,z=521288629,w=88675123;
    unsigned int t;
    t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}
★ 時間計測の関数を呼びすぎない。

サンプル2のように、10000回に1回だけ時間計測をするという方法があります（時間計測の関数はsystem_clockを使えば近年は十分速いですが）。時間計測以外の値でも、毎回計算する必要がないものがあればここで計算してしまいましょう。

    long long steps;    // 試行回数
    for (steps = 0; ; steps++)
    {
        if (steps % 10000 == 0)
        {
            const double time = duration_cast<microseconds>(system_clock::now() - startClock).count() * 1e-6;   // 経過時間（秒）
            if (time >= END_TIME)
            {
                break;
            }

ただ、TopCoderだとこの回数を増やしすぎて、制限時間オーバーしてしまうという可能性もあります。そういった点も考えて、適切な回数おきに時間をはかるようにしましょう。残り時間が少なくなったら、ループ回数を減らすといった安全策をとる方法もあります。

2. 次の状態（状態近傍）をどう決めるか？

ここは鍵になる部分ですが、コツとしてまとめるのは難しいところですね…。

★★★★ スコアが改善しやすくなるように、次の状態を選ぶ。


これは、問題ごとによく考えて対処しなければいけないので、あまり「コツ」というかんじではないのですが…。

[問題]
タスク数aの仕事Aと、タスク数bの仕事Bがある。
作業員はx人いて、それぞれ、1時間あたりにこなせる仕事Aのタスク数と、仕事Bのタスク数は決まっている。
ただし、仕事Aの仕事Bのどちらかしか行うことができない。
仕事Aと仕事Bの両方が終わるのに時間が最小になるように、作業員を仕事Aか仕事Bのどちらかに割り当てよ。

焼きなまし法で解くなら、まず最初にてきとーに作業員をA,Bに割り振り、ここからスコア＝作業時間として、改善させていきます。こういった問題を解くときに、次の状態として以下の２つを考えてみましょう。

  1. 移動：ある作業員をA→B または B→A に移動させる。
  2. 交換：Aの作業員の誰かと、Bの作業員の誰かを交換する。

移動２回は交換になることもあるので、とりあえず移動をさせておけば、交換もたまにしてくれるし良さそうに見えますが、そうとは限りません。この場合は、交換を主に行い、たまに移動させるぐらいがちょうど良いと思います。移動では、人数が偏り、片方の仕事の作業時間が増えやすくなります。２回移動すれば、交換になって、作業時間が減る可能性がありますが、

  • 移動：良スコア → 悪いスコア → とても良いスコア
  • 交換：良スコア 　　 → とても良いスコア

と考えると、移動だけだと、スコアが悪い状態を経由しないとダメなので、運がよくないと、とても良いスコアまで辿りつけません。その点、交換ならば、人員が偏った悪いスコアになりやすい状態をスキップできるので、スムーズにとても良いスコアまで改善できます。

もちろん、正しい人数割り当てもしなければいけないので、移動も少しは試したほうがいいですが、こういった点も考えて次の状態を決めると、スコアが改善しやすくなります。

★★★★ 近傍の種類を増やす

2021/03/06 追加しました。

近傍の種類は多くあればあるほど良いです。上記のサンプル2のコードには、近傍が２種類用意されています。

            // 近傍1 : d日目のコンテストをtに変えてみる
            const int d = randXor() % D;
            const int t = randXor() % NUM_TYPES;

            // 近傍2 : d0日目のコンテストとd1日目コンテストを交換してみる
            const int diff = randXor() % 10 + 1;    // d0とd1は10日差まで
            const int d0 = randXor() % (D - diff);
            const int d1 = d0 + diff;

サンプル2では、近傍1も近傍2も50%ずつの確率で選んでいますが、ここをどの近傍をどれぐらいの確率で試すかを調整すれば、性能を改善できるチャンスがあります。また、統計をとるのも大事です。遷移が成功した回数・失敗した回数のデータを取っておけば、どの近傍が有効なのかが分かり、問題を理解するヒントにもなります。さらに、「近傍1で遷移が成功すると、その直後で近傍2で遷移が成功しやすい」などの事実が分かれば、あらたに近傍3「近傍1のあとに近傍2」という新たな近傍を用意することもできるでしょう。

★★ Greedyに、次の状態を選ぶ。

先ほどの問題であれば、「仕事Aと仕事Bそれぞれから、常に作業時間が一番多くかかっている人を1人ずつ選び、交換」とGreedyに、次の状態を選ぶことができます。こういう選び方をすると、ベストスコアは序盤はどんどん伸びていきますが、局所最適に陥りやすいです。また、ありがちなのが、Greedyに選ぶアルゴリズムが、計算時間のボトルネックになってしまうパターンです。複雑なアルゴリズムより、ランダムでとにかくいろいろ試すのが焼きなまし法の強みです。Greedyにやるにせよ、多少ランダムを入れたり、計算時間はかからない単純な選び方をしたりする工夫が必要です。

★★ 序盤は大きく変化、終盤は小さく変化させたほうがいい。

単純にy=f(x)のyを最大化する問題なら、

  • 序盤(t=0に近いとき) 次のx'は x'=x+rand()%1001 - 500
  • 終盤(t=Tに近いとき) 次のx'は x'=x+rand()% 11 - 5

とか、
作業員割り当ての問題なら、

  • 序盤(t=0に近いとき) 交換を100回やる
  • 終盤(t=Tに近いとき) 交換を 1回やる

といったかんじです。スコアが局所最適でハマっている場合は有効です。ただ、焼きなまし法自体が、局所最適にはまりにくいアルゴリズムなので、あまり効果がない場合も多いです。

★★★★ 状態が大きく変化するわりには、スコアが変わらない近傍は良い。

2021/03/06 追加しました。

  • 状態変化が大きいけど、スコアが大きい変化しかしない（悪化する確率が圧倒的に高い）
  • 状態変化が小さいけど、スコアが小さい変化しかしない近傍

となりがちなのですが、問題によっては状態が大きく変化するわりには、スコアがあまり変化がしないケースがあります。こういう近傍があると、遷移が小さい近傍だけは辿りつけない状態にも短時間で行くことができ、焼きなましの性能がアップします。見逃さないようにしましょう。

★★★ 初期状態からすべての状態に行けるようにする。

2021/03/06 追加しました。
選んだ初期状態によっては、どんなに運がよくても用意した近傍ではすべての状態に辿りつけない場合があります。特に解に制約条件があるような問題（例：パズル問題など）で、このようなことが起こりやすいです。すべての状態にいけるような近傍を用意しましょう。

  • 無理にランダムに大きく状態が変化するような近傍を入れてしまうと、それ自体が遷移の成功確率を大きく減らすことにもなり、全体の性能を落とす可能性があるので、簡単ではありません。
  • 多点スタート（後述）にすれば、ある程度緩和されますが、多点スタートでは１つの解を求めるのにかけられる計算時間が減ってしまい、あまりよくありません。
  • スコア関数にペナルティ（後述）を入れる方法もあります。
★★ 良い初期状態からスタートする

これは問題に大きく依存する項目です。初期状態はてきとーでも良い場合と、すこし時間をかけて計算して良い状態から始めたほうが良い場合があります。特に、初期状態が悪いと状態を改善する近傍の選択肢が少なすぎる、逆に言えば、状態を改善すればするほどさらに改善できる可能性が増えるタイプの問題は、初期状態に少し時間を書けて求めて、それから焼きなましたほうが良い結果がでる場合があります。

3. スコア関数

スコア関数（scoring function）とも、評価関数（evaluation function）とも呼びます。

★★★★★ 良い状態を正しく評価できるような、スコア関数を新たに用意する

問題のスコア関数をそのまま使うと、焼きなまし法がうまく行かないことがあります。良い状態を正しく評価できるような、新たなスコア関数を用意しましょう。

[問題]
 （大きく略）…  i番目の人の満足値をs(i)とする。「満足度の最小値 min s(i) 」を、最大化しなさい。

この問題は「満足度の最小値 min s(i) 」を最大化する問題なので、素直に考えると、スコア＝「満足度の最小値 min s(i) 」としたくなりそうですが、これは良くありません。仮にそうした場合、

  状態A s(0)=  10, s(1)=10, s(2)=  10 のとき、スコアはmin{s(0),s(1),s(2)} = 10
  状態B s(0)=1000, s(1)=10, s(2)=1000 のとき、スコアはmin{s(0),s(1),s(2)} = 10

と状態Aと状態Bのスコアは同じになってしまいます。これはスコアは同じといえ、状態Bのほうが明らかに良い状態です。なぜなら、s(1)の点が仮に300点に増えたとしたら

  状態A' s(0)=  10, s(1)=300, s(2)=  10 のとき、スコアはmin{s(0),s(1),s(2)} =  10
  状態B' s(0)=1000, s(1)=300, s(2)=1000 のとき、スコアはmin{s(0),s(1),s(2)} = 300

と状態Bのほうは一気にスコアが増えます。しかし元のスコア関数のままだと状態Aと状態Bは同じスコアで、状態A→状態Bはスコア改善になりません。つまり焼きなまし法を行っても、状態A→状態Bには運が良くないと移動しません。
この問題であれば、例えば、
スコア＝1番小さいs(i) + 0.0001 * 2番目に小さいs(i)
とすれば

  状態A s(0)=  10, s(1)=10, s(2)=  10 のとき、スコアは 10 + 0.0001 *   10  = 10.001
  状態B s(0)=1000, s(1)=10, s(2)=1000 のとき、スコアは 10 + 0.0001 * 1000  = 10.1

と状態Bの方が好スコアになり、状態A→状態Bに移動しやすくなります。
ここは問題によって大きく有効かどうか大きく異なってきますが、必ず考慮に入れたほうがいいでしょう。結果に大差がつきやすいポイントです。

★★★スコアが改善しやすくなるようなスコア関数を決める。

上記の言い換えともいえるのですが、わりと焼きなまし法のスコア関数をつくるときに、問題のスコア関数をそのまま使ったり、違う状態を正しく評価しなかったりすると、同じスコアになりがちで、青い線のような飛び飛びの平坦な関数になることが多いです。赤いように、山登りしやすいスコア関数を作れないか考えましょう。必ずしも、問題のスコア関数の最大値と、焼きなまし法で使うスコア関数の最大値が一致する必要はありません。最大値をとりうるような状態（argmax）を一致させることが重要です。
f:id:shindannin:20150526014955p:plain

★★★スコア増分が0のときを、無くす。

2021/03/06 追加しました。
スコア増分が0のときは、焼きなまし法をやっても状態が行ったりきたりブラウン運動するだけで、解がなかなか改善しません。また、遷移するときの処理のほうが遷移しないときより重くなる場合、さらにひどいことになります。必ずスコア増分の履歴をチェックし、スコア増分が0になっているときはスコアの差をつけましょう。

★★★ 制約条件がある問題では、スコア関数にペナルティをつける。

2021/03/06 追加しました。
制約条件がある問題で、「制約条件を満たさない状態への遷移を許さない」とすると、ほとんど状態遷移されないままになってしまうことがあります。こういうときの対策としては、制約を満たしていない解も許して、ペナルティをつけて、スコア関数から減点します。

  • 制約条件を満たしている度合いがひどいほど、ペナルティは大きく減点します。
  • 計算の序盤ではペナルティは小さく、徐々にペナルティを大きくして、最後には絶対に制約条件を満たすようにします。
  • 最後にでてきた最適解が制約条件を満たしているか確認しましょう（ただ、そうならないようなペナルティのつけ方を見つけるのが理想です。）

4. 遷移確率と温度

たいていの焼きなまし法の解説書には、expを使った方法が載っています。
焼きなまし法 - Wikipedia

遷移確率 P(e, e' , T) は e' < e の時に 1 を返すよう定義されている（すなわち下りの移動は必ず実行される）。
そうでない場合の確率はexp((e−e' )/T) と定義 

注：wikipediaの引用なので、本文の設定と違います。
    ここでのTは温度（=時間とともに減っていく値）です。
    最小化問題なのでe−e'となってます（最大化問題ならe'-eです）。
★時刻tに置ける温度Tを調整する

最初の例の以下のコードの部分です。

     const bool FORCE_NEXT = R*(T-t)>T*(rand()%R);

上の例だと、

  • t=0のとき、常にtrue
  • t=Tのとき、常にfalse

になるような線形関数にしています。だいたいの場合は、これで問題ありませんが、バリエーションはいろいろ考えられます。実際には、ベストスコアの更新状況をみて、更新があまりに序盤と終盤に固まっているときは、２乗してみるとか、ロジスティック曲線をつかうとか、その程度で良いと思います。

★★★★exp(e-e')/Tを使う。（e-e':スコアの悪化量）

wikipediaでは以下のように書いていますが、exp(e-e')/Tの形にしたほうが良いでしょう。（理由は、次回の改定時に書きます）

焼きなまし法 - Wikipedia

関数 P は e' −e の値が増大する際には確率を減らす値を返すように設定される。
つまり、ちょっとしたエネルギー上昇の向こうに極小がある可能性の方が、
どんどん上昇している場合よりも高いという考え方である。しかし、この条件は
必ずしも必須ではなく、上記の必須条件が満たされていればよい。

ただ、このexpの式、分子はスコア差なので分かりますが、分母の温度をどう決めれば良いかは分かりづらいです。目安としては、たまには許容しても良いスコアの悪化量を温度とすると良いです。 最大の山を越えられるようにと考えると、考えられるスコアの最大値-最小値を初期温度を設定したほうが良いのですが、実際には温度調整して最終スコアがどうなるか実験してみてください。

double startTemp   =  100;
double endTemp     =   10;
double temp        = startTemp + (endTemp - startTemp) * t / T;
double probability = exp((nextScore - currentScore) / temp);  // 最大化なので、分子はnextScore - currentScoreで良い。
bool FORCE_NEXT    = probability > (double)(mXor.random() % R) / R;

時刻t=0（最初）でtemp=startTemp、時刻t=T（最後）でtemp=endTempとなります。

例えばスコアが100悪化したとき、nextScore - currentScore = -100なので、startTempを100にしてみます。endTempはとりあえず10にしてみましょう。
時間t=0であれば、exp(-100/100) = exp(- 1) = 0.36787944117と、3回に1回ぐらいは遷移します。
時間t=Tであれば、exp(-100/ 10) = exp(-10) = 0.00004539992と、まず遷移はしません。
このような感じでやると、良いと思います。

5. アレンジ

★★★★ 山登り

最初の例をFORCE_NEXT=0にすれば、山登り法になります。そもそも関数が凸であれば、山登り法も焼きなまし法も結果は一緒で、しかもrand()を使わない分、山登り法のほうが高速で、スコア改善のチャンスも多くなります。簡単に試せるので、試しましょう。

★ 多点スタート

焼きなまし法をつかっても、局所解にはまりやすい場合は、初期値を毎回かえた、多点スタートにする方法もあります。山登りと同時に使う方法もあります。ただ、多点スタートで結果がよくなるケースは、そもそもスコア関数や近傍が良くないので、そこを見直したほうがよいかもしれません。

★★★ 1変数ごとに焼きなまし

焼きなまし法で、たとえば、y = f(x1,x2,x3,x4,x5)と、状態変数が多いとき、次の状態のパターンが多くなり、結果的にあまり多くの状態を試せず、良い解が求められないときがあります。このような場合は、x2=x3=x4=x5を固定して、x1を焼きなまし、x1=x3=x4=x5を固定して、x2を焼きなまし、x1=x2=x4=x5を固定して、x3を焼きなましというふうに１変数ずつ焼きなまししていく方法があります。xの値がそれぞれ独立性が高い(?)場合は、この方法は有効です。

6. デバッグ

★ スコア差分の検算

2021/03/06 追加しました。
例えばサンプル2ではスコアの差分を計算していますが、いきなり書くとバグらせる可能性もあります。以下のような手順をとれば、バグがないか確認することができます。

  1. サンプル1のgetFullScoreのように、全体のスコアの関数を用意して、それが正しいか確認します。通常のコンテストなら提出することで確認できます。
  2. 差分スコアrealDeltaScore = getFullScore(遷移後) - getFullScore(遷移前)として計算します。
  3. サンプル2のようにdeltaScoreを直接求める方法ができたら、realDeltaScoreと一致するか確認します。
★ デバッグの際には、決定的にする

2021/03/06 追加しました。
状態やスコアが、まれにおかしくなるといったバグがあった場合、通常の焼きなまし法は非決定的なのでバグが再現ができない場合があります。そういったときは、デバッグ時には決定的にしましょう。

  • 温度を求めるときにsystem clockを使わずに、最大ループ回数を固定して進行割合で温度を求めるようにしましょう。
    const static long long MAX_STEPS   =  10000000; // 終了までの回数
    for (long long steps = 0; steps < MAX_STEPS ; steps++) // 試行回数
    {
            const double progressRatio = static_cast<double>(steps) / MAX_STEPS;   // 進捗。開始時が0.0、終了時が1.0
            temp = START_TEMP + (END_TEMP - START_TEMP) * progressRatio;
  • 乱数はシードさえ固定すれば決定的になるので問題ありません。
  • もしマルチスレッドを使っている場合は、グローバル変数とstatic変数に気をつけましょう。とくに乱数でstaticを使っている場合は乱数をクラス化するか、シングルスレッドで実行しましょう。

7. まだ書いていない内容

  • 並列化
    • 排他制御つきの焼きなまし
    • ロックフリーの焼きなまし
  • タブーサーチ + 焼きなまし
  • レプリカ交換法

リンク集

基礎編

AtCoder wataさんによる公式解説
https://img.atcoder.jp/intro-heuristics/editorial.pdf

Tsukammoさんの記事
qiita.com

gacinさんの記事
gasin.hatenadiary.jp

応用編

ats5515さんは、巨大近傍法（状態遷移の回数は少なめになるけど、計算時間をかけてより良い次の状態を探す）が特に強いです。近傍の計算時間に。焼きなまし法の近傍を探すのに2段目の焼きなまし、制約条件をほぼ満たす準validな状態、温度遷移の工夫など、他の人があまり工夫しないもところでもいろいろバリエーションがあります。
ats5515.hatenablog.com

hakomoさんの記事は、まさに詳細という内容です。探索空間の話は、この記事が一番詳しいと思います。事例へのリンクも多いです。
github.com