模擬退火法の概要

模擬退火（Simulated Annealing, SA）は、大域的最適解近似アルゴリズムの一種であり、組み合わせ最適化問題などに広く適用される確率的探索手法です。解空間が非常に広大であり、かつ目的関数が単峰性を持たない場合などに有効です。

この手法は、金属加工における「焼きなまし（Annealing）」プロセスに由来します。金属を高温で加熱し、徐々に冷却することで結晶構造を整え、エネルギー状態を安定させる物理現象を模倣しています。アルゴリズムにおいては、目的関数の値をエネルギーとみなし、温度パラメータを制御しながら解の探索を行います。

アルゴリズムの仕組み

模擬退火の核心は、現在の状態よりも悪い状態であっても、一定の確率で受け入れる点にあります。これにより、局所最適解に陥ることを防ぎ、大域的最適解への到達確率を高めます。

現在の状態を $S$、新たに生成された状態を $S'$、それぞれのエネルギー（目的関数値）の差を $\Delta E = E(S') - E(S)$ とします。現在の温度を $T$ としたとき、状態遷移の受諾確率 $P$ はメトロポリス基準により以下のように定義されます。

温度 $T$ が高い初期段階では悪い解も受け入れやすく、探索範囲が広がります。温度が下がるにつれて受諾確率が低下し、最終的には勾配降下法のように動作して解が収束します。

冷却スケジュールとパラメータ

実装においては、以下の 3 つのパラメータを設定する必要があります。

• 初期温度 ($T_{\text{start}}$): 十分に大きな値を設定し、初期段階で広範囲を探索できるようにします。
• 冷却係数 ($D$): 1 よりわずかに小さい実数（例：0.98 〜 0.99）を設定し、温度の減衰率を制御します。
• 終了温度 ($T_{\text{end}}$): 0 に近い非常に小さな値（例：1e-3）を設定し、これ以下になったら探索を終了します。

アルゴリズムは、温度が $T_{\text{end}}$ を下回るまで、状態の生成・評価・遷移を繰り返します。また、精度向上のため、探索过程中に見つけた最良の解を別途保持しておくことが推奨されます。

実装事例

以下に、模擬退火法を用いた代表的な問題の解決例を示します。

事例 1: 重心バランス問題 (P1337)

複数の重りが座標平面上に配置されており、糸でつながれているときの平衡点（重心に近い点）を求める問題です。全系のポテンシャルエネルギーが最小となる点を探します。

#include <iostream>
#include <cmath>
#include <algorithm>
#include <ctime>
#include <vector>
#include <iomanip>

using namespace std;

struct WeightPoint {
    double x, y, w;
};

int n;
vector<WeightPoint> points;
double best_x, best_y, min_energy;

// エネルギー計算関数
double compute_energy(double x, double y) {
    double energy = 0.0;
    for (const auto& p : points) {
        double dx = p.x - x;
        double dy = p.y - y;
        energy += sqrt(dx * dx + dy * dy) * p.w;
    }
    if (energy < min_energy) {
        min_energy = energy;
        best_x = x;
        best_y = y;
    }
    return energy;
}

void simulated_annealing() {
    double temperature = 100000.0;
    double decay = 0.998;
    double current_x = best_x;
    double current_y = best_y;

    while (temperature > 1e-3) {
        double next_x = current_x + temperature * (2.0 * rand() / RAND_MAX - 1.0);
        double next_y = current_y + temperature * (2.0 * rand() / RAND_MAX - 1.0);
        
        double delta = compute_energy(next_x, next_y) - compute_energy(current_x, current_y);
        
        if (delta < 0 || exp(-delta / temperature) > (double)rand() / RAND_MAX) {
            current_x = next_x;
            current_y = next_y;
        }
        temperature *= decay;
    }
    
    // 終端処理：低温付近で微調整
    for(int i=0; i<1000; ++i){
        double next_x = best_x + temperature * (2.0 * rand() / RAND_MAX - 1.0);
        double next_y = best_y + temperature * (2.0 * rand() / RAND_MAX - 1.0);
        compute_energy(next_x, next_y);
    }
}

int main() {
    srand(time(0));
    cin >> n;
    points.resize(n);
    double sum_x = 0, sum_y = 0;
    for (int i = 0; i < n; ++i) {
        cin >> points[i].x >> points[i].y >> points[i].w;
        sum_x += points[i].x;
        sum_y += points[i].y;
    }
    // 初期値を重心に設定
    best_x = sum_x / n;
    best_y = sum_y / n;
    min_energy = compute_energy(best_x, best_y);
    
    simulated_annealing();
    
    cout << fixed << setprecision(3) << best_x << " " << best_y << endl;
    return 0;
}

事例 2: 爆撃範囲の最適化 (P5544)

建物を避けて爆弾を投げるとき、最も多くの敵を巻き込める位置を求める問題です。建物との距離制限内で敵の数を最大化します。

#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <ctime>

using namespace std;

struct Building { double x, y, r; };
struct Enemy { double x, y; };

int n, m;
double radius_limit;
vector<Building> buildings;
vector<Enemy> enemies;
int max_killed = 0;

double get_dist(double x1, double y1, double x2, double y2) {
    return sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2));
}

int count_enemies(double x, double y) {
    double safe_dist = radius_limit;
    for (const auto& b : buildings) {
        safe_dist = min(safe_dist, get_dist(x, y, b.x, b.y) - b.r);
    }
    if (safe_dist < 0) return 0;
    
    int count = 0;
    for (const auto& e : enemies) {
        if (get_dist(x, y, e.x, e.y) <= safe_dist) count++;
    }
    return count;
}

void solve_sa(double &x, double &y) {
    double temp = radius_limit;
    int current_score = count_enemies(x, y);
    
    while (temp > 1e-3) {
        double nx = x + 2 * temp * (2.0 * rand() / RAND_MAX - 1.0);
        double ny = y + 2 * temp * (2.0 * rand() / RAND_MAX - 1.0);
        int next_score = count_enemies(nx, ny);
        
        // 受諾判定
        if (next_score > current_score || exp(10000.0 * (next_score - current_score) / temp) > (double)rand() / RAND_MAX) {
            x = nx;
            y = ny;
            current_score = next_score;
        }
        max_killed = max(max_killed, current_score);
        temp *= 0.989;
    }
}

int main() {
    srand(time(0));
    cin >> n >> m >> radius_limit;
    buildings.resize(n);
    enemies.resize(m);
    for (int i = 0; i < n; ++i) cin >> buildings[i].x >> buildings[i].y >> buildings[i].r;
    for (int i = 0; i < m; ++i) cin >> enemies[i].x >> enemies[i].y;
    
    // 複数回試行して精度を上げる
    for (int i = 0; i < 25; ++i) {
        double sx = (2.0 * rand() / RAND_MAX - 1.0) * radius_limit;
        double sy = (2.0 * rand() / RAND_MAX - 1.0) * radius_limit;
        solve_sa(sx, sy);
    }
    cout << max_killed << endl;
    return 0;
}

事例 3: 球形空間生成器 (P4035)

n 次元空間上の点群から、それらへの距離が等しくなる中心点を探す問題です。ここでは勾配のような力を模擬退火で適用します。

#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>

using namespace std;

int n;
vector<vector<double>> coords;
vector<double> center;

void run_sa() {
    double temp = 10000.0;
    while (temp > 1e-4) {
        vector<double> forces(n, 0.0);
        double total_dist = 0.0;
        vector<double> dists;
        
        // 各点までの距離を計算
        for (int i = 0; i < n + 1; ++i) {
            double d = 0.0;
            for (int j = 0; j < n; ++j) {
                d += pow(coords[i][j] - center[j], 2);
            }
            d = sqrt(d);
            dists.push_back(d);
            total_dist += d;
        }
        double avg_dist = total_dist / (n + 1);
        
        // 合力を計算
        for (int i = 0; i < n + 1; ++i) {
            for (int j = 0; j < n; ++j) {
                forces[j] += (dists[i] - avg_dist) * (coords[i][j] - center[j]) / avg_dist;
            }
        }
        
        // 温度に応じて移動
        for (int j = 0; j < n; ++j) {
            center[j] += forces[j] * temp;
        }
        temp *= 0.99995;
    }
}

int main() {
    cin >> n;
    coords.resize(n + 1, vector<double>(n));
    center.assign(n, 0.0);
    
    for (int i = 0; i < n + 1; ++i) {
        for (int j = 0; j < n; ++j) {
            cin >> coords[i][j];
            center[j] += coords[i][j];
        }
    }
    // 初期値を平均に
    for (int j = 0; j < n; ++j) center[j] /= (n + 1);
    
    run_sa();
    
    cout << fixed << setprecision(3);
    for (int j = 0; j < n; ++j) cout << center[j] << " ";
    cout << endl;
    return 0;
}

事例 4: 金貨の分配 (P3878)

金貨の山を 2 つのグループに分け、その合計価値の差を最小化する問題です。配列のシャッフルと要素の交換を繰り返します。

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <ctime>
#include <climits>

using namespace std;

int n;
vector<long long> coins;
long long min_diff = LLONG_MAX;

long long calculate_diff() {
    long long sum1 = 0, sum2 = 0;
    for (int i = 0; i < n / 2; ++i) sum1 += coins[i];
    for (int i = n / 2; i < n; ++i) sum2 += coins[i];
    return abs(sum1 - sum2);
}

void execute_sa() {
    double temp = 100.0;
    long long current_diff = calculate_diff();
    min_diff = min(min_diff, current_diff);
    
    while (temp > 1e-8) {
        int i = rand() % n;
        int j = rand() % n;
        if (i == j) continue;
        
        swap(coins[i], coins[j]);
        long long next_diff = calculate_diff();
        min_diff = min(min_diff, next_diff);
        
        if (exp((double)(current_diff - next_diff) / temp) < (double)rand() / RAND_MAX) {
            swap(coins[i], coins[j]); // 戻す
        } else {
            current_diff = next_diff;
        }
        temp *= 0.989;
    }
}

int main() {
    srand(time(0));
    int t;
    cin >> t;
    while (t--) {
        cin >> n;
        coins.resize(n);
        min_diff = LLONG_MAX;
        for (int i = 0; i < n; ++i) cin >> coins[i];
        
        if (n == 1) {
            cout << coins[0] << endl;
            continue;
        }
        
        // 複数回試行
        for (int i = 0; i < 50; ++i) {
            random_shuffle(coins.begin(), coins.end());
            execute_sa();
        }
        cout << min_diff << endl;
    }
    return 0;
}

事例 5: マトリックスの彩色 (P3936)

グリッド上のマスを指定された色の数で塗り分け、隣接するマスで色が異なる境界線の数を最小化する問題です。交換によるコスト計算の高速化が重要です。

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <ctime>

using namespace std;

int H, W, C;
vector<int> color_counts;
vector<vector<int>> grid;
vector<vector<int>> best_grid;
int min_cost = 1e9;
int dx[4] = {0, 1, 0, -1};
int dy[4] = {1, 0, -1, 0};

void init_grid() {
    int current_color = 0;
    int count = 0;
    for (int i = 0; i < H; ++i) {
        for (int j = 0; j < W; ++j) {
            if (count < color_counts[current_color]) {
                grid[i][j] = current_color;
                count++;
            } else {
                current_color++;
                grid[i][j] = current_color;
                count = 1;
            }
        }
    }
}

int evaluate_cost() {
    int cost = 0;
    for (int i = 0; i < H; ++i) {
        for (int j = 0; j < W; ++j) {
            for (int k = 0; k < 4; ++k) {
                int ni = i + dx[k];
                int nj = j + dy[k];
                if (ni >= 0 && ni < H && nj >= 0 && nj < W) {
                    if (grid[i][j] != grid[ni][nj]) cost++;
                }
            }
        }
    }
    return cost / 2;
}

void update_best() {
    best_grid = grid;
    min_cost = evaluate_cost();
}

void run_sa() {
    double temp = 1.0;
    int current_cost = min_cost;
    
    while (temp > 1e-15) {
        int x1 = rand() % H;
        int y1 = rand() % W;
        int x2 = rand() % H;
        int y2 = rand() % W;
        
        if (grid[x1][y1] == grid[x2][y2]) continue;
        
        // 交換前の境界コスト計算（局所的）
        int delta_before = 0;
        for(int k=0; k<4; ++k){
            int nx = x1 + dx[k], ny = y1 + dy[k];
            if(nx>=0 && nx<H && ny>=0 && ny<W && grid[nx][ny] != grid[x1][y1]) delta_before++;
            nx = x2 + dx[k], ny = y2 + dy[k];
            if(nx>=0 && nx<H && ny>=0 && ny<W && grid[nx][ny] != grid[x2][y2]) delta_before++;
        }
        
        swap(grid[x1][y1], grid[x2][y2]);
        
        // 交換後の境界コスト計算
        int delta_after = 0;
        for(int k=0; k<4; ++k){
            int nx = x1 + dx[k], ny = y1 + dy[k];
            if(nx>=0 && nx<H && ny>=0 && ny<W && grid[nx][ny] != grid[x1][y1]) delta_after++;
            nx = x2 + dx[k], ny = y2 + dy[k];
            if(nx>=0 && nx<H && ny>=0 && ny<W && grid[nx][ny] != grid[x2][y2]) delta_after++;
        }
        
        int new_cost = current_cost - delta_before + delta_after;
        
        if (new_cost <= current_cost || exp((double)(current_cost - new_cost) / temp) > (double)rand() / RAND_MAX) {
            current_cost = new_cost;
            if (current_cost < min_cost) update_best();
        } else {
            swap(grid[x1][y1], grid[x2][y2]); // 戻す
        }
        temp *= 0.99999;
    }
}

int main() {
    srand(time(0));
    cin >> H >> W >> C;
    color_counts.resize(C);
    grid.assign(H, vector<int>(W));
    best_grid.assign(H, vector<int>(W));
    
    for (int i = 0; i < C; ++i) cin >> color_counts[i];
    
    init_grid();
    update_best();
    
    for (int i = 0; i < 3; ++i) run_sa();
    
    for (int i = 0; i < H; ++i) {
        for (int j = 0; j < W; ++j) {
            cout << best_grid[i][j] << " ";
        }
        cout << endl;
    }
    return 0;
}