模拟退火算法
模拟退火（Simulated Annealing, SA）是一种模拟金属退火物理过程的启发式优化算法，专门用来解决「容易陷入局部最优」的复杂优化问题。它最天才的地方在于一个简单的接受准则：变好的解一定接受，变差的解也有一定概率接受——正是这个「允许变差」的机制，让它能像高温金属里的原子一样翻越能量壁垒，最终冷却到全局能量最低的状态。本文从物理原理、Metropolis 准则、降温调度讲起，覆盖适用场景、评价指标、四张核心诊断图，并给出可运行的手写实现（连续函数优化 + TSP 组合优化，均与精确解/库实现对照验证）。
一、算法含义
1.1 物理背景：金属退火
把金属加热到高温再缓慢冷却，原子在高温下有足够动能随机移动（可以「跳」到能量更高的位置），随着温度下降，原子逐渐被束缚在能量最低的晶体结构上。如果冷却太快（淬火），原子来不及找到最优排列，就会冻结在亚稳态——这就是「局部最优」的物理原型。
1.2 核心：Metropolis 接受准则
把物理退火映射到优化问题：
| 物理概念 | 优化问题 |
|---|---|
| 原子的位置状态 | 解 |
| 系统的能量 | 目标函数值 （求最小） |
| 温度 | 控制参数（逐步下降） |
| 原子随机移动 | 邻域扰动 |
模拟退火的灵魂是 Metropolis 准则——产生新解 后，接受它的概率为：
解读这张核心公式：
- （新解更优）：，必接受——永远不拒绝进步；
- （新解更差）：以 的概率「容忍」变差。温度 越高，差解越容易被接受（高温=大范围探索）； 越低，差解越难被接受（低温=精细收敛）。
与贪心算法的本质区别：贪心「只进不退」，一旦掉进局部最优就出不来；SA 用「以一定概率变差」换取跳出局部最优的机会。
1.3 降温调度（Cooling Schedule）
温度按几何方式逐层下降：
- 初始温度 ：要足够高，使初期接受率接近 100%（高温充分「熔化」解空间）；
- 降温系数 ：越接近 1 降温越慢、搜索越充分、耗时越长（α 是速度与质量的权衡旋钮，见 3.4 节敏感性分析）；
- 终止温度 ：降到足够低后停止；
- 内循环次数：每个温度下做多次扰动采样（马尔可夫链），保证该温度下达到「热平衡」。
1.4 算法流程（伪代码）
1. 初始化：随机生成初始解 x0，温度 T = T0
2. while T > T_min: # 外层：逐层降温
3. repeat inner_iters 次: # 内层：同温采样
4. x' = neighbor(x) # 邻域扰动生成新解
5. ΔE = f(x') - f(x)
6. if ΔE ≤ 0 或 random() < exp(-ΔE/T):
7. x = x' # Metropolis 准则
8. if f(x) < f(x_best): x_best = x # 精英保留
9. T = α · T # 几何降温
10. return x_best
1.5 优缺点
- 优点：原理简单、实现容易；理论上能以概率 1 收敛到全局最优（降温足够慢时）；对目标函数零数学要求（黑盒即可）；连续与组合问题通吃。
- 缺点：收敛慢（需要大量迭代）；参数（、、内循环次数）敏感；单次运行有随机性，必须报告多次运行统计。
二、何时使用（适用场景与条件）
2.1 适用场景
| 场景 | 例子 |
|---|---|
| 组合优化 | 旅行商（TSP）、车辆路径（VRP）、排班调度、背包 |
| 连续多峰优化 | 神经网络参数寻优、复杂函数的参数拟合 |
| 黑盒优化 | 目标函数无法求导、甚至没有显式形式时 |
| 逃离局部最优 | 贪心/爬山法卡住后，用 SA 再优化 |
2.2 竞赛典型题目
- 物流配送路径规划（TSP/VRP 类）；
- 车间调度、人员排班（组合爆炸、难以精确求解）；
- 模型的非线性参数寻优（目标函数复杂、多峰）；
- 「先精确后启发」的混合题：小规模精确验证 + 大规模 SA 求解。
2.3 使用前提
- 能定义目标函数与邻域结构即可（几乎无数学要求）；
- 有足够的计算时间预算（SA 是「以时间换质量」的算法）。
2.4 不适用情形与对比选择
- 问题规模小、可精确求解时（如本文 TSP 的 10 城实例可直接穷举），优先精确算法；
- 实时性要求高（毫秒级决策）时，SA 太慢；
- 对参数调优无耐心时，可考虑「少参数」的 PSO（但 PSO 也易早熟，各有利弊）。
| 算法 | 适合问题 | 相对 SA 的特点 |
|---|---|---|
| 遗传算法 | 组合优化、多峰连续 | 种群并行、全局性强，但参数更多、实现更复杂 |
| 粒子群算法 | 连续优化 | 收敛快、实现简单，但易早熟 |
| 蚁群算法 | 路径类组合优化 | 正反馈强，仅适合路径/图类问题 |
| 模拟退火 | 连续+组合通吃、黑盒 | 参数最少、原理最直观，收敛较慢 |
三、算法指标
3.1 指标汇总表
| 指标 | 公式/含义 | 解读 |
|---|---|---|
| 历史最优目标值 | （精英保留） | 越小越好；是最终汇报的核心指标 |
| 与已知最优的差距 gap | （或相对量） | 小规模问题可与穷举/精确解对照，gap=0 证明实现正确 |
| 温度轨迹 | 高温段波动大（探索）、低温段平稳（收敛） | |
| 接受率 | 被接受提议数 / 总提议数 | 高温应接近 1；随降温逐渐下降，过低说明「退化成贪心」 |
| 差解接受率 | 被接受的变差提议 / 变差提议总数 | 体现「允许变差」机制是否有效 |
| 多次运行统计 | ， | 随机算法必须报告：μ 体现平均质量，σ 体现稳定性 |
| 收敛率 | 达标运行数 / 总运行数 | 如「 的运行比例」，比均值更直观 |
| 降温系数 α 敏感性 | 不同 α 下的 分布 | α 越大质量越好但层数越多（见 3.4） |
3.2 温度轨迹与目标值轨迹（本文图 1）
横轴为累计迭代次数，左轴温度（对数坐标）呈锯齿下降，右轴目标值：高温段目标值大幅波动（正在四处探索），低温段趋于平稳（已经收敛）。若低温段目标值还在剧烈波动，说明降温太快或内循环不足。
3.3 接受率曲线（本文图 2）
高温段接受率接近 1（几乎来者不拒），随着降温逐步下降；低温段差解接受率趋近 0（只接受改进解）。若全程接受率都很低 → 相当于贪心，失去全局搜索能力；若全程都很高 → 降温无效。
3.4 降温系数 α 的敏感性
本文实测（单轮降温、每个 α 独立运行 10 次）：
| α | 降温层数 | 平均 | 达标次数（） |
|---|---|---|---|
| 0.85 | 57 | 0.4550 | 1/10 |
| 0.90 | 88 | 0.1556 | 2/10 |
| 0.95 | 180 | 0.2005 | 2/10 |
| 0.99 | 917 | 0.0117 | 6/10 |
规律：α 越接近 1，降温越慢、搜索越充分、解质量越好——代价是降温层数与耗时成倍增长（917 层 vs 57 层）。实践中常取 0.9~0.99，配合「多轮重升温 + 精修」策略（见 6.2 节）平衡速度与质量。
四、可视化图表
| 图名 | 用途 | 关键解读点 |
|---|---|---|
① 温度-目标值双轴图（sa_temperature_curve.png） | 展示整个退火过程的探索→收敛节奏 | 温度锯齿下降、目标值「高温大波动→低温平稳」；精修阶段（橙色阴影）目标值进一步压实 |
② 接受率曲线（sa_acceptance_rate.png） | 检验 Metropolis 机制是否正常工作 | 高温段接受率≈1、低温段→0；差解接受率单调下降；异常时全程过低（退化成贪心）或过高（降温失效） |
③ Rastrigin 等高线+搜索路径（sa_rastrigin_path.png） | 直观展示「跳出局部最优」 | 白色轨迹点穿越多个局部极小盆地，最终停在原点（全局最优）；红点为最终解 |
④ TSP 路径对比图（sa_tsp_path.png） | 展示 SA 解的质量 | SA 路径（红色实线）与穷举最优一致、无交叉；最近邻路径（灰色虚线）有交叉且更长 |
好图特征：图 1 的锯齿结构清晰、图 2 的两条曲线「高开低走」、图 3 的轨迹横跨多个盆地、图 4 的 SA 路径无自交。 异常特征：图 1 目标值一路单调下降无波动 → 内循环过少/降温太快；图 3 轨迹只在一个盆地内打转 → 初始温度太低。
五、符号说明
| 符号 | 含义 | 示例/单位 |
|---|---|---|
| 当前解（连续问题为向量，TSP 为路径序列） | ||
| 目标函数值（求最小） | ||
| 能量差 | — | |
| / | 初始/终止温度 | ， |
| 降温系数 | ||
| 当前温度 | 逐层下降 | |
| inner_iters | 内循环次数（同温采样数） | 200 |
| Metropolis 接受概率 | ||
| 历史最优目标值（精英保留） | ||
| gap | 与已知最优解的差距 | TSP：0.0000% |
| / | 多次运行最优值的均值/标准差 |
六、可运行程序（完整代码）
环境要求：Python 3.12，依赖 numpy、scipy、matplotlib、pandas（本目录
.venv已装好）。以下代码块按顺序拼接保存为一个.py文件即可运行：控制台打印第三节全部指标与对照结果，并在figures/子目录生成第四节 4 张图（sa_前缀）。无图形界面时可用MPLBACKEND=Agg运行。
# -*- coding: utf-8 -*-
"""
============================================================
模拟退火算法完整示例：手写 SA（Rastrigin + TSP）+ 库实现对照
------------------------------------------------------------
实例 1：手写 SA 求解 2D Rastrigin 函数（[-5.12, 5.12]^2，全局最优 f(0,0)=0）
——含多轮重升温 + 最终低温精修、10 次独立运行统计、单轮降温随机性对照、
α 敏感性分析、与 scipy.optimize.dual_annealing 的对照
实例 2：手写 SA（2-opt 邻域）求解 10 城 TSP（np.random.seed(42) 生成坐标），
排列穷举精确最优对照 + 最近邻贪心对照
输出：控制台打印全部指标 + figures/ 目录下 4 张图（sa_ 前缀）
依赖：numpy、scipy、matplotlib、pandas（无外部文件依赖）
============================================================
"""
# ========== 0. 导入库与全局设置 ==========
import copy
import os
from itertools import permutations
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import dual_annealing
# ---- matplotlib 中文显示设置（防止图内中文乱码）----
plt.rcParams["font.sans-serif"] = ["PingFang SC", "Arial Unicode MS", "SimHei"]
plt.rcParams["axes.unicode_minus"] = False
# ---- 图片输出目录（相对当前工作目录的 figures/ 子目录）----
os.makedirs("figures", exist_ok=True)
# ========== 1. 目标函数与邻域生成 ==========
BOUND = 5.12 # Rastrigin 的常用搜索域边界
def rastrigin(x):
"""Rastrigin 函数：f(x) = 10n + Σ[x_i^2 - 10cos(2πx_i)]，n=2 时全局最优 f(0,0)=0。
搜索域内约有 100 个局部极小点（整数格点），是经典的"多峰"测试函数。
同时支持单点输入 (2,) 与网格输入 (H, W, 2)（画等高线用）。"""
x = np.asarray(x, dtype=float)
return 20.0 + np.sum(x ** 2 - 10.0 * np.cos(2.0 * np.pi * x), axis=-1)
def rastrigin_neighbor(x, T, rng, sigma_floor=0.004, sigma_scale=0.15):
"""连续问题的邻域生成：高斯扰动，步长 σ 随温度收缩。
σ = max(σ_floor, σ_scale·T)：高温大步长负责全局搜索、跨越壁垒；
低温小步长负责在盆地底部精细下降；越界点截断回搜索域内。
sigma_floor=0.004 用于主降温轮；最终精修轮改用更小的 0.002。"""
sigma = max(sigma_floor, sigma_scale * T)
return np.clip(x + rng.standard_normal(len(x)) * sigma, -BOUND, BOUND)
# ========== 2. 手写模拟退火求解器（完整类） ==========
class SimulatedAnnealing:
"""手写模拟退火求解器：邻域生成 + Metropolis 接受准则 + 几何降温。
objective : callable，f(x) -> float，目标函数（求最小值）
neighbor : callable，neighbor(x, T, rng) -> x_new，邻域生成（可依赖温度 T）
"""
def __init__(self, objective, neighbor):
self.objective = objective
self.neighbor = neighbor
def anneal(self, x0, T0, T_min, alpha, inner_iters, seed=None,
neighbor=None, record=False, collect_path=False):
"""执行一轮完整降温（外层：温度逐层几何下降；内层：同温度马尔可夫采样）。
返回字典：
x_final, f_final : 末态解与目标值
x_best, f_best : 本轮历史最优解（精英保留）
n_levels : 降温层数
n_proposals / n_accepted / n_worse / n_accepted_worse : 接受率统计
history（record=True）: 逐层记录 T、累计迭代数、f_cur、f_best、
每层接受率、每层差解接受率（画图 1、图 2 用）
path（collect_path=True）: 每个被接受解的坐标（画图 3 用）
"""
rng = np.random.default_rng(seed)
neighbor = neighbor if neighbor is not None else self.neighbor
x = copy.deepcopy(x0)
fx = float(self.objective(x))
x_best, f_best = copy.deepcopy(x), fx
T = T0
n_levels = 0
n_proposals = n_accepted = n_worse = n_accepted_worse = 0
path = [] if collect_path else None
history = ({"T": [], "evals": [], "f_cur": [], "f_best": [],
"acc_rate": [], "worse_acc_rate": []} if record else None)
while T > T_min: # 外层循环：温度每轮乘 α
acc_lv = worse_lv = acc_worse_lv = 0
for _ in range(inner_iters): # 内层循环：同一温度下的采样
x_new = neighbor(x, T, rng)
f_new = float(self.objective(x_new))
dE = f_new - fx # 能量差 ΔE = E(x_new) - E(x)
n_proposals += 1
if dE > 0: # 统计"差解"（变差的提议）
n_worse += 1
worse_lv += 1
# ---- Metropolis 准则：P(接受) = min(1, exp(-ΔE/T)) ----
# 变好（ΔE<=0）必接受；变差（ΔE>0）以 exp(-ΔE/T) 的概率接受
if dE <= 0 or rng.random() < np.exp(-dE / T):
n_accepted += 1
acc_lv += 1
if dE > 0:
n_accepted_worse += 1
acc_worse_lv += 1
x, fx = x_new, f_new
if collect_path:
path.append(copy.deepcopy(x))
if fx < f_best: # 精英保留：记住历史最优
x_best, f_best = copy.deepcopy(x), fx
if record:
history["T"].append(T)
history["evals"].append(n_proposals)
history["f_cur"].append(fx)
history["f_best"].append(f_best)
history["acc_rate"].append(acc_lv / inner_iters)
history["worse_acc_rate"].append(
acc_worse_lv / worse_lv if worse_lv > 0 else 0.0)
T *= alpha # 几何降温：T_{k+1} = α·T_k
n_levels += 1
return {"x_final": x, "f_final": fx, "x_best": x_best, "f_best": f_best,
"n_levels": n_levels, "n_proposals": n_proposals,
"n_accepted": n_accepted, "n_worse": n_worse,
"n_accepted_worse": n_accepted_worse,
"history": history, "path": path}
def solve(self, x0, n_cycles=12, T0=10.0, T_min=1e-3, alpha=0.95, inner_iters=200,
polish=True, polish_neighbor=None, polish_T0=0.5, polish_T_min=1e-5,
polish_alpha=0.9, polish_inner=100, seed=None, record=False,
collect_path_first_cycle=False):
"""多轮重升温求解 + 可选最终低温精修（实践中常用的组合策略）。
每轮降温都从"当前全局最优"重新升温开始：重升温给搜索多次跨越能量壁垒
的机会，显著降低陷入局部最优的概率；最后的低温精修（小步长、慢降温）
把解抛光到所在盆地的底部。返回值同 anneal，另附：
cycles : 每轮 anneal 的结果列表（含精修轮）
polished : 是否执行了精修轮
history（record=True）: 各轮历史按累计迭代数拼接后的总轨迹
"""
rng = np.random.default_rng(seed)
gbest_x, gbest_f = copy.deepcopy(x0), float(self.objective(x0))
cycles = []
for c in range(n_cycles):
r = self.anneal(gbest_x, T0, T_min, alpha, inner_iters,
seed=int(rng.integers(0, 2 ** 31)), record=record,
collect_path=(collect_path_first_cycle and c == 0))
cycles.append(r)
if r["f_best"] < gbest_f:
gbest_f = r["f_best"]
gbest_x = copy.deepcopy(r["x_best"])
polished = False
if polish:
r = self.anneal(gbest_x, polish_T0, polish_T_min, polish_alpha, polish_inner,
seed=int(rng.integers(0, 2 ** 31)),
neighbor=polish_neighbor, record=record)
cycles.append(r)
polished = True
if r["f_best"] < gbest_f:
gbest_f = r["f_best"]
gbest_x = copy.deepcopy(r["x_best"])
hist = None
if record: # 拼接各轮历史（迭代数加偏移）
hist = {"T": [], "evals": [], "f_cur": [], "f_best": [],
"acc_rate": [], "worse_acc_rate": []}
offset = 0
for r in cycles:
h = r["history"]
for k in ("T", "f_cur", "f_best", "acc_rate", "worse_acc_rate"):
hist[k].append(np.asarray(h[k], dtype=float))
hist["evals"].append(np.asarray(h["evals"], dtype=float) + offset)
offset += r["n_proposals"]
for k in hist:
hist[k] = np.concatenate(hist[k])
return {"x_best": gbest_x, "f_best": gbest_f,
"x_final": copy.deepcopy(cycles[-1]["x_final"]),
"f_final": cycles[-1]["f_final"],
"cycles": cycles, "polished": polished,
"n_proposals": sum(r["n_proposals"] for r in cycles),
"n_accepted": sum(r["n_accepted"] for r in cycles),
"n_worse": sum(r["n_worse"] for r in cycles),
"n_accepted_worse": sum(r["n_accepted_worse"] for r in cycles),
"history": hist, "path": cycles[0]["path"]}
def show(x):
"""小数值改用科学计数法显示，避免打印成误导性的 0.000000"""
return f"{x:.6g}" if (x == 0 or 1e-4 <= abs(x) < 1e6) else f"{x:.2e}"
# ========== 3. 实例 1：手写 SA 求解 2D Rastrigin 函数 ==========
print("=" * 68)
print("【实例 1】模拟退火求解 2D Rastrigin 函数")
print("搜索域 [-5.12, 5.12]^2，已知全局最优 f(0, 0) = 0")
print("主参数：T0=10，T_min=1e-3，α=0.95，内循环=200，重升温轮数=12，精修(T0=0.5，α=0.9)")
sa_global = SimulatedAnnealing(rastrigin, rastrigin_neighbor) # 主降温轮
polish_nb = lambda x, T, rng: rastrigin_neighbor(x, T, rng, sigma_floor=0.002) # 精修轮
rng = np.random.default_rng(42)
x0 = rng.uniform(-BOUND, BOUND, 2)
print(f"随机起点 x0 = [{x0[0]:.4f}, {x0[1]:.4f}]，f(x0) = {rastrigin(x0):.3f}")
res_main = sa_global.solve(x0, n_cycles=12, seed=42, record=True,
collect_path_first_cycle=True,
polish=True, polish_neighbor=polish_nb)
gb_f = np.inf
gb_x = None
for c, r in enumerate(res_main["cycles"]):
tag = "精修阶段" if (res_main["polished"] and c == len(res_main["cycles"]) - 1) \
else f"第 {c + 1} 轮降温"
if r["f_best"] < gb_f:
gb_f, gb_x = r["f_best"], r["x_best"]
print(f" {tag:8s}：本轮最优 f* = {show(r['f_best']):>10s}，"
f"全局最优 f* = {show(gb_f):>10s} @ x = [{gb_x[0]:+.4f}, {gb_x[1]:+.4f}]")
print("-" * 68)
print("【主运行指标汇总】")
print(f"历史最优目标值 f* = {show(res_main['f_best'])}，最优解 x* = "
f"[{res_main['x_best'][0]:.2e}, {res_main['x_best'][1]:.2e}]")
print(f"与已知全局最优 0 的差距 gap = {show(res_main['f_best'] - 0.0)}")
print(f"温度轨迹：T0 = {show(10.0)} 起，每层乘 α = 0.95，{res_main['cycles'][0]['n_levels']} 层/轮，"
f"前 6 层 T = " + ", ".join(f"{t:.2f}" for t in res_main["history"]["T"][:6]) + " ...")
print(f"总迭代（提议）次数 = {res_main['n_proposals']}"
f"（{sum(r['n_levels'] for r in res_main['cycles'])} 层降温）")
print(f"总接受率 = {res_main['n_accepted'] / res_main['n_proposals'] * 100:.1f}%，"
f"差解接受率 = {res_main['n_accepted_worse'] / max(res_main['n_worse'], 1) * 100:.1f}%")
# ---- 3.1 多次运行统计（完整协议，10 次独立运行） ----
print("=" * 68)
print("【多次运行统计】完整协议（12 轮重升温 + 精修），10 次独立运行（seed=0..9）")
f_vals = []
for s in range(10):
rng_s = np.random.default_rng(s)
x0_s = rng_s.uniform(-BOUND, BOUND, 2)
r = sa_global.solve(x0_s, n_cycles=12, seed=s, polish=True, polish_neighbor=polish_nb)
f_vals.append(r["f_best"])
f_vals = np.array(f_vals)
df_runs = pd.DataFrame({"seed": list(range(10)),
"f*": [show(v) for v in f_vals],
"是否达标(f*<0.01)": ["是" if v < 0.01 else "否" for v in f_vals]})
print(df_runs.to_string(index=False))
print(f"均值 = {show(f_vals.mean())}，标准差 = {show(f_vals.std(ddof=1))}，"
f"最小值 = {show(f_vals.min())}，最大值 = {show(f_vals.max())}")
print(f"收敛率（f* < 0.01） = {int((f_vals < 0.01).sum())}/10")
# ---- 3.2 对照：只做一轮降温（不重升温、不精修）----
print("=" * 68)
print("【随机性对照】只用一轮降温（无重升温、无精修），10 次独立运行（seed=100..109）")
single = []
for s in range(100, 110):
rng_s = np.random.default_rng(s)
x0_s = rng_s.uniform(-BOUND, BOUND, 2)
r = sa_global.anneal(x0_s, T0=10.0, T_min=1e-3, alpha=0.95, inner_iters=200, seed=s)
single.append(r["f_best"])
single = np.array(single)
print("各次 f*：" + ", ".join(show(v) for v in single))
print(f"收敛率（f* < 0.01） = {int((single < 0.01).sum())}/10，"
f"平均 f* = {show(single.mean())}（可见 SA 的随机性：部分运行陷入局部最优）")
# ---- 3.3 降温系数 α 敏感性 ----
print("=" * 68)
print("【降温系数 α 敏感性】单轮降温，每个 α 独立运行 10 次（同一组 seed=200..209）")
rows = []
for alpha in (0.85, 0.90, 0.95, 0.99):
fs = []
for s in range(200, 210):
rng_s = np.random.default_rng(s)
x0_s = rng_s.uniform(-BOUND, BOUND, 2)
r = sa_global.anneal(x0_s, T0=10.0, T_min=1e-3, alpha=alpha, inner_iters=200, seed=s)
fs.append(r["f_best"])
fs = np.array(fs)
rows.append({"降温系数 α": alpha,
"平均 f*": f"{fs.mean():.4f}",
"最好 f*": show(fs.min()),
"达标次数": f"{int((fs < 0.01).sum())}/10",
"降温层数": int(np.ceil(np.log(1e-3 / 10.0) / np.log(alpha)))})
print(pd.DataFrame(rows).to_string(index=False))
# ---- 3.4 库实现对照：scipy.optimize.dual_annealing ----
print("=" * 68)
print("【库实现对照】scipy.optimize.dual_annealing 求解同一问题（seed=42，maxiter=300）")
res_da = dual_annealing(rastrigin, bounds=[(-BOUND, BOUND)] * 2, seed=42, maxiter=300)
print(f"f* = {show(res_da.fun)}，x* = [{res_da.x[0]:.2e}, {res_da.x[1]:.2e}]，"
f"函数评估次数 nfev = {res_da.nfev}")
# ========== 4. 实例 2：手写 SA 求解 10 城 TSP（2-opt 邻域） ==========
print("=" * 68)
print("【实例 2】模拟退火（2-opt 邻域）求解 10 城 TSP")
np.random.seed(42)
cities = np.random.rand(10, 2) * 100 # 10 个城市坐标（seed 固定，可复现）
n_city = 10
d = np.linalg.norm(cities[:, None, :] - cities[None, :, :], axis=2) # 距离矩阵 d_ij
print("城市坐标（np.random.seed(42) 生成）：")
print(pd.DataFrame({"城市": list(range(n_city)),
"x": cities[:, 0].round(2),
"y": cities[:, 1].round(2)}).to_string(index=False))
# ---- 4.1 排列穷举求精确最优（固定城市 0 为起点，9! = 362880 条路径） ----
perms = np.array(list(permutations(range(1, n_city))), dtype=int)
full = np.hstack([np.zeros((len(perms), 1), dtype=int), perms])
legs = np.linalg.norm(cities[full[:, :-1]] - cities[full[:, 1:]], axis=2).sum(axis=1)
closing = np.linalg.norm(cities[full[:, -1]] - cities[full[:, 0]], axis=1)
total = legs + closing
idx = int(np.argmin(total))
exact_len = float(total[idx])
exact_route = [int(v) for v in full[idx]]
print(f"穷举精确最优：长度 = {exact_len:.4f}，"
f"路径 = {' → '.join(map(str, exact_route))}")
# ---- 4.2 最近邻贪心（对比基准） ----
def nn_route(start=0):
"""最近邻贪心：每次去往最近的未访问城市（局部最优，不一定全局最优）"""
unvisited = set(range(n_city))
unvisited.discard(start)
route, cur = [start], start
while unvisited:
nxt = min(unvisited, key=lambda j: d[cur, j])
route.append(nxt)
unvisited.discard(nxt)
cur = nxt
return route
nn_r = nn_route(0)
nn_len = sum(d[nn_r[i], nn_r[(i + 1) % n_city]] for i in range(n_city))
print(f"最近邻贪心：长度 = {nn_len:.4f}（比精确最优长 {(nn_len / exact_len - 1) * 100:.2f}%），"
f"路径 = {' → '.join(map(str, nn_r))}")
# ---- 4.3 手写 SA（2-opt 邻域） ----
def tour_length(route):
"""路径总长度（目标函数）"""
return sum(d[route[i], route[(i + 1) % n_city]] for i in range(n_city))
def two_opt_neighbor(route, T, rng):
"""2-opt 邻域：随机截取一段子路径并反转（等价于交换两条边、消除交叉）"""
r = list(route)
while True:
i, j = int(rng.integers(0, n_city)), int(rng.integers(0, n_city))
if i > j:
i, j = j, i
if j - i >= 2:
r[i:j] = r[i:j][::-1]
return r
sa_tsp = SimulatedAnnealing(tour_length, two_opt_neighbor)
print("SA 求解 TSP：T0=100，T_min=0.01，α=0.95，内循环=150，独立运行 5 次")
best_len, best_route = np.inf, None
for run in range(5):
r0 = list(np.random.default_rng(1000 + run).permutation(n_city)) # 随机初始路径
r = sa_tsp.anneal(r0, T0=100.0, T_min=0.01, alpha=0.95, inner_iters=150,
seed=2000 + run)
gap = (r["f_best"] - exact_len) / exact_len * 100
print(f" 第 {run + 1} 次：f* = {r['f_best']:.4f}，与穷举最优的差距 gap = {gap:.2f}%"
+ (f"，接受率 = {r['n_accepted'] / r['n_proposals'] * 100:.1f}%，"
f"差解接受率 = {r['n_accepted_worse'] / max(r['n_worse'], 1) * 100:.1f}%"
if run == 0 else ""))
if r["f_best"] < best_len:
best_len, best_route = r["f_best"], [int(v) for v in r["x_best"]]
print(f"SA 最优：长度 = {best_len:.4f}，gap = {(best_len / exact_len - 1) * 100:.4f}%，"
f"路径 = {' → '.join(map(str, best_route))}")
# ========== 5. 图 1：温度-迭代双轴图（左轴温度，右轴目标值） ==========
hist = res_main["history"]
evals = hist["evals"] / 1e4 # 横轴单位：万次迭代
gb_curve = np.minimum.accumulate(hist["f_best"]) # 全局历史最优曲线
stage_offsets = np.cumsum([r["n_proposals"] for r in res_main["cycles"]]) / 1e4
polish_start = stage_offsets[-2] if res_main["polished"] else None
fig, ax1 = plt.subplots(figsize=(10, 5))
ax1.plot(evals, hist["T"], color="#1f77b4", lw=1.2)
ax1.set_yscale("log")
ax1.set_xlabel("累计迭代次数（万次）")
ax1.set_ylabel("温度 T（对数坐标）", color="#1f77b4")
ax1.tick_params(axis="y", labelcolor="#1f77b4")
if polish_start is not None:
ax1.axvspan(polish_start, evals[-1], color="orange", alpha=0.15)
ax2 = ax1.twinx()
ax2.plot(evals, hist["f_cur"], color="gray", lw=0.6, alpha=0.55)
ax2.plot(evals, gb_curve, color="red", lw=1.6)
ax2.set_ylabel("目标值 f(x) / f*", color="red")
ax2.tick_params(axis="y", labelcolor="red")
lines = [plt.Line2D([0], [0], color="#1f77b4", lw=1.2, label="温度 T（左轴）"),
plt.Line2D([0], [0], color="gray", lw=0.8, alpha=0.6, label="当前目标值 f(x)（右轴）"),
plt.Line2D([0], [0], color="red", lw=1.6, label="历史最优 f*（右轴）"),
plt.Rectangle((0, 0), 1, 1, color="orange", alpha=0.15, label="精修阶段")]
ax1.legend(handles=lines, loc="upper left", fontsize=9)
ax1.set_title("模拟退火降温过程：温度（左轴）与目标值（右轴）轨迹，重升温形成锯齿")
plt.tight_layout()
plt.savefig("figures/sa_temperature_curve.png", dpi=150)
# ========== 6. 图 2：接受率随迭代变化曲线 ==========
acc = hist["acc_rate"]
wacc = hist["worse_acc_rate"]
w = 21 # 滑动平均窗口（每层 200 次迭代）
kernel = np.ones(w) / w
acc_s = np.convolve(acc, kernel, mode="valid")
wacc_s = np.convolve(wacc, kernel, mode="valid")
ev = evals[(w - 1) // 2: (w - 1) // 2 + len(acc_s)]
fig, ax = plt.subplots(figsize=(10, 4.5))
ax.plot(evals, acc, color="#1f77b4", lw=0.5, alpha=0.35)
ax.plot(ev, acc_s, color="#1f77b4", lw=1.8, label="整体接受率（滑动平均）")
ax.plot(evals, wacc, color="#d62728", lw=0.5, alpha=0.35)
ax.plot(ev, wacc_s, color="#d62728", lw=1.8, label="差解接受率（滑动平均）")
if polish_start is not None:
ax.axvspan(polish_start, evals[-1], color="orange", alpha=0.15, label="精修阶段")
ax.set_xlabel("累计迭代次数（万次）")
ax.set_ylabel("接受率")
ax.set_ylim(0, 1.02)
ax.set_title("接受率随迭代的变化：高温段大量接受差解，低温段几乎只接受改进解")
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig("figures/sa_acceptance_rate.png", dpi=150)
# ========== 7. 图 3：2D Rastrigin 等高线 + 搜索路径散点 ==========
fig, ax = plt.subplots(figsize=(8.5, 7.5))
g = np.linspace(-BOUND, BOUND, 400)
X, Y = np.meshgrid(g, g)
Z = rastrigin(np.stack([X, Y], axis=-1))
cf = ax.contourf(X, Y, Z, levels=40, cmap="viridis")
cb = fig.colorbar(cf, ax=ax, shrink=0.8)
cb.set_label("目标值 f(x)")
path = res_main["path"]
if path:
P = np.array(path)
ax.scatter(P[:, 0], P[:, 1], s=2, c="white", alpha=0.35,
label=f"搜索轨迹（第 1 轮降温，{len(path)} 个被接受点）")
ax.scatter(x0[0], x0[1], marker="o", s=90, color="lime", edgecolors="black",
zorder=5, label="起点")
ax.scatter(res_main["x_best"][0], res_main["x_best"][1], marker="*", s=260,
color="gold", edgecolors="black", zorder=6, label="历史最优解")
ax.scatter(res_main["x_final"][0], res_main["x_final"][1], marker="o", s=130,
color="red", edgecolors="black", zorder=6, label="最终解（红点）")
ax.set_xlim(-BOUND, BOUND)
ax.set_ylim(-BOUND, BOUND)
ax.set_xlabel("x1")
ax.set_ylabel("x2")
ax.set_title(f"2D Rastrigin 等高线与 SA 搜索路径（最终 f* = {show(res_main['f_best'])}）")
ax.legend(loc="upper right", fontsize=9)
plt.tight_layout()
plt.savefig("figures/sa_rastrigin_path.png", dpi=150)
# ========== 8. 图 4：小规模 TSP 最优路径图（与最近邻对比） ==========
fig, ax = plt.subplots(figsize=(8.5, 7.5))
def draw_route(route, color, lw, ls="-"):
"""按路径顺序画带箭头的折线"""
pts_r = cities[route]
for i in range(n_city):
j = (i + 1) % n_city
ax.annotate("", xy=pts_r[j], xytext=pts_r[i],
arrowprops=dict(arrowstyle="-|>", color=color, lw=lw,
linestyle=ls, shrinkA=7, shrinkB=7, alpha=0.9))
draw_route(nn_r, color="gray", lw=1.1, ls=(0, (5, 4)))
draw_route(best_route, color="#d62728", lw=2.0)
for i in range(n_city):
is_start = (i == best_route[0])
ax.scatter(cities[i, 0], cities[i, 1], s=170 if is_start else 120,
color="red" if is_start else "#1f77b4",
marker="s" if is_start else "o", zorder=5, edgecolors="white")
ax.annotate(str(i), (cities[i, 0], cities[i, 1]),
xytext=(cities[i, 0] + 1.5, cities[i, 1] + 1.5), fontsize=10)
handles = [plt.Line2D([0], [0], color="#d62728", lw=2.0,
label=f"SA 最优路径（长度 {best_len:.2f}，与穷举最优一致）"),
plt.Line2D([0], [0], color="gray", lw=1.1, ls=(0, (5, 4)),
label=f"最近邻路径（长度 {nn_len:.2f}，长 {(nn_len / exact_len - 1) * 100:.1f}%）"),
plt.Line2D([0], [0], marker="s", color="w", markerfacecolor="red",
markersize=9, label="起点城市（SA 路径首城）")]
ax.legend(handles=handles, loc="best", fontsize=9)
ax.set_xlabel("x 坐标")
ax.set_ylabel("y 坐标")
ax.set_title(f"10 城 TSP：SA（2-opt）最优路径 vs 最近邻贪心路径（穷举最优 = {exact_len:.4f}）")
ax.set_aspect("equal")
plt.tight_layout()
plt.savefig("figures/sa_tsp_path.png", dpi=150)
# ========== 9. 显示全部图形 ==========
plt.show()
print("=" * 68)
print("4 张图已保存到 figures/ 目录：")
for f in ("sa_temperature_curve.png", "sa_acceptance_rate.png",
"sa_rastrigin_path.png", "sa_tsp_path.png"):
print(" figures/" + f)
七、结果解读与注意事项
7.1 运行结果解读（对照真实输出）
运行第六节代码，关键输出如下（第七节所有数值均来自实际运行）：
实例 1（Rastrigin 主运行，seed=42）：随机起点 ，第 1 轮降温即从 31.887 降至 0.2354，第 4 轮降到 ，精修后最终 （与全局最优 0 的 gap 仅为 ）。总迭代 44.23 万次（2263 层降温），总接受率 41.9%，差解接受率 26.5%。
多次运行统计：完整协议（12 轮重升温 + 精修）10 次独立运行全部达标（收敛率 10/10）， 均值 、标准差 ——质量与稳定性俱佳。而单轮降温对照（无重升温、无精修）10 次运行只有 3/10 达标、平均 ——这正是「重升温 + 精修」策略的价值：重升温给搜索多次翻越壁垒的机会，低温精修把解抛光到盆地底部。
α 敏感性：单轮降温下 α 从 0.85 增大到 0.99，达标次数 1/10 → 6/10、平均 0.4550 → 0.0117，但降温层数从 57 层膨胀到 917 层——降温越慢质量越好，耗时同步增长。
库实现对照：scipy 的 dual_annealing 在同一问题上得到 （仅 1333 次评估），说明成熟库实现在效率上远优于朴素手写版——竞赛中若时间紧，直接用库即可；手写实现的价值在于理解原理与定制邻域。
实例 2（10 城 TSP）：穷举精确最优 290.3068；最近邻贪心 312.2388（长 7.55%，贪心局部最优的典型）；手写 SA（2-opt 邻域）5 次独立运行全部与穷举最优一致（gap = 0.0000%）——在小规模问题上验证了实现正确性，这是竞赛中「先小规模精确验证、再大规模启发求解」的标准打法。
7.2 常见坑
- 降温太快 → 退化成贪心：α 太小（如 0.7）或内循环太少，SA 失去全局搜索能力。判断信号：接受率曲线全程很低（见 3.3）。
- 初始温度过高 → 浪费时间： 太高时前期几十层都在随机游走。经验： 取使初始接受率在 0.8~1.0 附近即可。
- 内循环不足 → 未达热平衡：同一温度下采样次数太少，无法充分探索该温度的等能面。
- 邻域设计不当：扰动步长不随温度收缩（连续问题），低温阶段解永远「够不着」精确最优；TSP 若用「交换两城市」邻域而非 2-opt，收敛会慢很多。
- 只跑一次就下结论：SA 是随机算法，单次结果没有统计意义，必须报告多次运行的均值/标准差（本文 3.1 节）。
- 不报告 gap：小规模问题不做精确解对照，无法证明实现正确性。
7.3 竞赛论文写作建议（话术模板）
针对……问题建立模型后，目标函数为多峰非线性函数、传统算法易陷入局部最优，故采用模拟退火算法求解。算法参数设置：初始温度 ，降温系数 ，内循环 200 次，采用「多轮重升温 + 低温精修」策略。在 10 城小规模实例上，算法 5 次运行均与穷举最优解一致（gap=0），验证了实现的正确性；大规模实例独立运行 10 次，最优解均值 …、标准差 …，收敛率 10/10，算法稳定可靠。降温系数敏感性分析表明 α=0.95 兼顾了解质量与计算时间。
八、延伸阅读
- 爬山算法（Hill Climbing）：SA 的「0 温度」特例（永远不接受差解），可作对比基线；
- 禁忌搜索（Tabu Search）：用「禁忌表」记忆近期访问的解来避免循环，组合优化中与 SA 互补；
- 遗传模拟退火混合：种群进化 + Metropolis 接受准则结合；
- 量子退火（Quantum Annealing）：利用量子隧穿效应的物理退火（提及即可）；
- scipy 库：
scipy.optimize.dual_annealing（广义模拟退火，竞赛中直接可用的成熟实现）。