焼きなまし法（Simulated Annealing）の仕組みとPython実装

import numpy as np

np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100


def tour_length(tour, cities):
    n = len(tour)
    d = 0.0
    for i in range(n):
        d += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return d


def two_opt_delta(tour, cities, i, j):
    """2-optスワップ（tour[i:j+1]を反転）による距離変化量を差分計算"""
    n = len(tour)
    a, b = tour[i - 1], tour[i]
    c, d = tour[j], tour[(j + 1) % n]
    old = np.linalg.norm(cities[a] - cities[b]) + np.linalg.norm(cities[c] - cities[d])
    new = np.linalg.norm(cities[a] - cities[c]) + np.linalg.norm(cities[b] - cities[d])
    return new - old


def two_opt_swap(tour, i, j):
    new_tour = tour.copy()
    new_tour[i : j + 1] = tour[i : j + 1][::-1]
    return new_tour


def run_with_rule(cities, init_tour, proposals, accept_rule):
    """proposals: (i, j, u, T) の共通提案列。accept_rule だけを差し替えて比較する"""
    tour = init_tour.copy()
    current = tour_length(tour, cities)
    best_dist = current
    decisions = []
    for i, j, u, T in proposals:
        delta = two_opt_delta(tour, cities, i, j)
        accept = accept_rule(delta, T, u)
        decisions.append(accept)
        if accept:
            tour = two_opt_swap(tour, i, j)
            current += delta
            best_dist = min(best_dist, current)
    return best_dist, decisions


# --- 共通の提案列を事前生成 ---
max_iter = 20000
rng = np.random.RandomState(7)
init_tour = list(rng.permutation(n_cities))
proposals_template = []
for _ in range(max_iter):
    i, j = sorted(rng.choice(n_cities, 2, replace=False))
    if i == 0:
        i = 1
    if j <= i:
        j = min(i + 1, n_cities - 1)
    u = rng.random()
    proposals_template.append((i, j, u))

# --- 温度 T を固定したSA vs 貪欲法（同一提案列） ---
for T_const in [1.0, 1e-3, 1e-6, 1e-9]:
    proposals = [(i, j, u, T_const) for (i, j, u) in proposals_template]

    def sa_rule(delta, T, u):
        return delta <= 0 or u < np.exp(-delta / T)

    def greedy_rule(delta, T, u):
        return delta <= 0

    sa_dist, sa_decisions = run_with_rule(cities, init_tour, proposals, sa_rule)
    greedy_dist, greedy_decisions = run_with_rule(cities, init_tour, proposals, greedy_rule)
    mismatches = sum(a != b for a, b in zip(sa_decisions, greedy_decisions))

    print(
        f"T={T_const:.0e}  SA={sa_dist:.4f}  greedy={greedy_dist:.4f}  mismatches={mismatches}/{max_iter}"
    )

import numpy as np

np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100


def tour_length(tour, cities):
    n = len(tour)
    d = 0.0
    for i in range(n):
        d += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return d


def two_opt_delta(tour, cities, i, j):
    n = len(tour)
    a, b = tour[i - 1], tour[i]
    c, d = tour[j], tour[(j + 1) % n]
    old = np.linalg.norm(cities[a] - cities[b]) + np.linalg.norm(cities[c] - cities[d])
    new = np.linalg.norm(cities[a] - cities[c]) + np.linalg.norm(cities[b] - cities[d])
    return new - old


def two_opt_swap(tour, i, j):
    new_tour = tour.copy()
    new_tour[i : j + 1] = tour[i : j + 1][::-1]
    return new_tour


def sa_tsp_track_acceptance(cities, T0, alpha, max_iter, seed, n_bins=10):
    rng = np.random.RandomState(seed)
    n = len(cities)
    tour = list(rng.permutation(n))
    current = tour_length(tour, cities)
    best_dist = current
    T = T0
    bin_size = max_iter // n_bins
    accept_counts = np.zeros(n_bins)
    propose_counts = np.zeros(n_bins)

    for k in range(max_iter):
        i, j = sorted(rng.choice(n, 2, replace=False))
        if i == 0:
            i = 1
        if j <= i:
            continue
        delta = two_opt_delta(tour, cities, i, j)
        b = min(k // bin_size, n_bins - 1)
        propose_counts[b] += 1
        if delta <= 0 or rng.random() < np.exp(-delta / T):
            tour = two_opt_swap(tour, i, j)
            current += delta
            accept_counts[b] += 1
            best_dist = min(best_dist, current)
        T *= alpha

    acc_rate_per_bin = accept_counts / np.maximum(propose_counts, 1)
    overall_acc = accept_counts.sum() / propose_counts.sum()
    return best_dist, overall_acc, acc_rate_per_bin


max_iter = 20000
alpha = 0.9995
for T0 in [0.01, 1, 10, 100, 1000, 10000, 100000]:
    dists, overalls, bins_all = [], [], []
    for seed in range(5):
        best_dist, overall_acc, acc_bins = sa_tsp_track_acceptance(cities, T0, alpha, max_iter, seed)
        dists.append(best_dist)
        overalls.append(overall_acc)
        bins_all.append(acc_bins)
    bins_mean = np.mean(bins_all, axis=0)
    print(
        f"T0={T0:g}  overall_acc={np.mean(overalls)*100:.2f}%  "
        f"final_dist={np.mean(dists):.2f}  "
        f"bins(0-10%,40-50%,90-100%)={bins_mean[0]*100:.1f}%,{bins_mean[4]*100:.1f}%,{bins_mean[9]*100:.1f}%"
    )

import numpy as np
import matplotlib.pyplot as plt

def rastrigin(x):
    """Rastrigin関数"""
    return 10 * len(x) + np.sum(x**2 - 10 * np.cos(2 * np.pi * x))

def simulated_annealing(objective, dim, bounds=(-5.12, 5.12),
                        T_init=100.0, alpha=0.995, max_iter=10000):
    """焼きなまし法による連続最適化"""
    # 初期解
    x = np.random.uniform(*bounds, size=dim)
    best_x = x.copy()
    best_f = objective(x)
    current_f = best_f

    T = T_init
    history = [best_f]

    for i in range(max_iter):
        # 近傍解の生成（ガウスノイズ）
        sigma = 0.5 * (T / T_init)  # 温度に比例してステップ幅を調整
        x_new = x + np.random.normal(0, sigma, size=dim)
        x_new = np.clip(x_new, *bounds)

        new_f = objective(x_new)
        delta = new_f - current_f

        # メトロポリス基準
        if delta <= 0 or np.random.random() < np.exp(-delta / T):
            x = x_new
            current_f = new_f

            if current_f < best_f:
                best_x = x.copy()
                best_f = current_f

        T *= alpha
        history.append(best_f)

    return best_x, best_f, history

# --- 実行 ---
np.random.seed(42)
best_x, best_f, history = simulated_annealing(rastrigin, dim=10)

print(f"最良解の適応度: {best_f:.6f}")

# --- 収束曲線のプロット ---
plt.figure(figsize=(10, 5))
plt.plot(history)
plt.xlabel('Iteration')
plt.ylabel('Best Objective Value')
plt.title('Simulated Annealing on Rastrigin Function (10D)')
plt.yscale('log')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def total_distance(tour, cities):
    """巡回路の総距離を計算"""
    n = len(tour)
    dist = 0
    for i in range(n):
        dist += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return dist

def two_opt_swap(tour, i, j):
    """2-optスワップ: tour[i:j+1]を反転"""
    new_tour = tour.copy()
    new_tour[i:j+1] = tour[i:j+1][::-1]
    return new_tour

def sa_tsp(cities, T_init=1000.0, alpha=0.9995, max_iter=100000):
    """焼きなまし法によるTSP"""
    n = len(cities)
    tour = list(range(n))
    np.random.shuffle(tour)

    current_dist = total_distance(tour, cities)
    best_tour = tour.copy()
    best_dist = current_dist
    T = T_init
    history = [best_dist]

    for _ in range(max_iter):
        # ランダムな2-optスワップ
        i, j = sorted(np.random.choice(n, 2, replace=False))
        new_tour = two_opt_swap(tour, i, j)
        new_dist = total_distance(new_tour, cities)

        delta = new_dist - current_dist
        if delta <= 0 or np.random.random() < np.exp(-delta / T):
            tour = new_tour
            current_dist = new_dist
            if current_dist < best_dist:
                best_tour = tour.copy()
                best_dist = current_dist

        T *= alpha
        history.append(best_dist)

    return best_tour, best_dist, history

# --- 実行 ---
np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100

best_tour, best_dist, history = sa_tsp(cities)
print(f"最短距離: {best_dist:.2f}")

# --- 結果の可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 巡回路
tour_cities = cities[best_tour + [best_tour[0]]]
axes[0].plot(tour_cities[:, 0], tour_cities[:, 1], 'b-o', markersize=5)
axes[0].set_title(f'Best Tour (distance={best_dist:.2f})')
axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)

# 収束曲線
axes[1].plot(history)
axes[1].set_xlabel('Iteration')
axes[1].set_ylabel('Best Distance')
axes[1].set_title('SA Convergence on TSP')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np

np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100


def tour_length(tour, cities):
    n = len(tour)
    d = 0.0
    for i in range(n):
        d += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return d


def two_opt_delta(tour, cities, i, j):
    n = len(tour)
    a, b = tour[i - 1], tour[i]
    c, d = tour[j], tour[(j + 1) % n]
    old = np.linalg.norm(cities[a] - cities[b]) + np.linalg.norm(cities[c] - cities[d])
    new = np.linalg.norm(cities[a] - cities[c]) + np.linalg.norm(cities[b] - cities[d])
    return new - old


def two_opt_swap(tour, i, j):
    new_tour = tour.copy()
    new_tour[i : j + 1] = tour[i : j + 1][::-1]
    return new_tour


def sa_tsp_schedule(cities, T0, schedule, max_iter, seed, T_min=1e-6):
    rng = np.random.RandomState(seed)
    n = len(cities)
    tour = list(rng.permutation(n))
    current = tour_length(tour, cities)
    best_dist = current
    checkpoints = {}
    cp_iters = sorted(set(int(max_iter * f) for f in (0.1, 0.5, 1.0)))

    for k in range(max_iter):
        if schedule == "log":
            c = T0 * np.log(2)
            T = max(c / np.log(k + 2), T_min)
        elif schedule == "exp":
            T = max(T0 * 0.9995**k, T_min)
        elif schedule == "linear":
            T = max(T0 - (T0 / max_iter) * k, T_min)

        i, j = sorted(rng.choice(n, 2, replace=False))
        if i == 0:
            i = 1
        if j <= i:
            continue
        delta = two_opt_delta(tour, cities, i, j)
        if delta <= 0 or rng.random() < np.exp(-delta / T):
            tour = two_opt_swap(tour, i, j)
            current += delta
            best_dist = min(best_dist, current)

        if (k + 1) in cp_iters:
            checkpoints[k + 1] = best_dist

    return best_dist, checkpoints


T0 = 1000.0
max_iter = 100000
for schedule in ["log", "exp", "linear"]:
    finals, cps_list = [], []
    for seed in range(5):
        best_dist, cps = sa_tsp_schedule(cities, T0, schedule, max_iter, seed)
        finals.append(best_dist)
        cps_list.append(cps)
    finals = np.array(finals)
    iters = sorted(cps_list[0].keys())
    cp_means = {it: np.mean([cps_list[s][it] for s in range(5)]) for it in iters}
    print(f"{schedule}: final={finals.mean():.3f}±{finals.std():.3f}  checkpoints={cp_means}")