scikit-opt 模拟退火三变体实战：Fast、Cauchy、Boltzmann 的数学原理与 Python 实现

本文以 docs/en/more_sa.md 为主线，深入讲解 scikit-opt 中模拟退火（Simulated Annealing, SA）算法的三种温度退火变体——Fast、Cauchy、Boltzmann。文中不仅给出每个变体的完整可运行示例（含边界约束版本），还结合 sko/SA.py 源码逐行拆解其状态更新公式、降温策略与 Metropolis 接受准则，帮助你理解不同退火计划的差异、参数取值的影响，以及如何在自己的最优化问题中直接落地使用。读完本文，你将能够：1）准确区分三种变体的采样与降温机制；2) 熟练使用 SAFast、SABoltzmann、SACauchy 三类完成带/不带边界的连续函数优化；3) 通过 best_x、best_y、best_y_history 等属性分析与可视化收敛过程。

一、模拟退火算法速览：三类变体的共性与差异

模拟退火的核心思想源自金属退火工艺：在高温下系统剧烈扰动，允许接受"变差"的解以跳出局部最优；随着温度下降，扰动幅度逐渐收敛，最终稳定在全局最优附近。scikit-opt 在 sko/SA.py 中实现了同一套算法骨架，只是**新解生成方式（get_new_x）与温度下降方式（cool_down）**不同，由此派生出三种变体：

变体	类名	新解扰动分布	降温公式
Fast（快速）	SAFast	基于符号函数的指数式跳跃	T_new = T0 * exp(-c * k^quench)
Cauchy（柯西）	SACauchy	柯西分布（tan(u) 重尾采样）	T_new = T0 / (1 + k)
Boltzmann（玻尔兹曼）	SABoltzmann	正态分布扰动	T_new = T0 / log(1 + k)

三个类都继承自 SimulatedAnnealingValue，而后者又继承自 SimulatedAnnealingBase（实现公共的 run() 主循环、Metropolis 接受准则与停止条件）。因此，只需要替换类名，即可在三类退火策略之间自由切换，其余 API 完全一致。

另外需要说明：在 sko/SA.py 末尾，SA = SAFast，即不带具体后缀的 SA 类默认就是 Fast 变体，这是示例中最常用的入门写法。

二、三种变体的数学原理与源码对应

原文档 docs/en/more_sa.md 用伪代码给出了三种退火计划的更新规则。下面逐一对照源码实现进行解读（相关实现均在 sko/SA.py 中）。

2.1 Fast 快速退火

原文档公式（Fast schedule）：

u ~ Uniform(0, 1, size = d)
y = sgn(u - 0.5) * T * ((1 + 1/T)**abs(2*u - 1) - 1.0)

xc = y * (upper - lower)
x_new = x_old + xc

c = n * exp(-n * quench)
T_new = T0 * exp(-c * k**quench)

源码对应 SAFast.get_new_x 与 SAFast.cool_down：

def get_new_x(self, x):
    r = np.random.uniform(-1, 1, size=self.n_dim)
    xc = np.sign(r) * self.T * ((1 + 1.0 / self.T) ** np.abs(r) - 1.0)
    x_new = x + xc * self.hop          # hop 相当于公式中的 (upper - lower)
    if self.has_bounds:
        return np.clip(x_new, self.lb, self.ub)
    return x_new

def cool_down(self):
    self.T = self.T_max * np.exp(-self.c * self.iter_cycle ** self.quench)

两者本质上是同一套公式：文档中的 u ~ Uniform(0,1) 与源码的 r ~ Uniform(-1,1) 可通过 r = 2u - 1 互相换算，sgn(u - 0.5) 对应 np.sign(r)；文档里的 y * (upper - lower) 在源码中体现为 xc * hop（当给定 lb/ub 时，hop 默认等于 ub - lb，见下文边界约束一节）。Fast 的扰动在高温时跨度很大（能大范围跳跃探索），随温度下降迅速收敛，因此得名"快速"。

关键点是 Fast 变体的降温和扰动是"指数-幂"耦合：源码构造函数中：

self.m, self.n, self.quench = kwargs.get('m', 1), kwargs.get('n', 1), kwargs.get('quench', 1)
self.c = self.m * np.exp(-self.n * self.quench)

即退火常数 c = m * exp(-n * quench)，默认 m = n = quench = 1。从源码结构看，Fast 变体实际读取的命名参数是 m、n、quench，而原文档示例中传入的 q=0.99 并不会参与 Fast 变体的计算（会被静默忽略），这也是使用文档示例时需要注意的实现细节。

2.2 Cauchy 柯西退火

原文档公式（Cauchy schedule）：

u ~ Uniform(-pi/2, pi/2, size=d)
xc = learn_rate * T * tan(u)
x_new = x_old + xc

T_new = T0 / (1 + k)

源码对应 SACauchy.get_new_x 与 SACauchy.cool_down：

def get_new_x(self, x):
    u = np.random.uniform(-np.pi / 2, np.pi / 2, size=self.n_dim)
    xc = self.learn_rate * self.T * np.tan(u)
    x_new = x + xc
    if self.has_bounds:
        return np.clip(x_new, self.lb, self.ub)
    return x_new

def cool_down(self):
    self.T = self.T_max / (1 + self.iter_cycle)

Cauchy 变体利用 tan(u)（u 在 (-π/2, π/2) 均匀采样）产生重尾分布的扰动：绝大多数扰动步长很小，但偶尔会产生非常大的跳跃，这种长尾特性使其在高维、多峰问题上更容易翻越"山脊"。其降温方式是最简单的倒数式降温 T0 / (1 + k)，速度较慢但稳健。learn_rate 是学习率系数，源码默认值为 0.5，它直接放缩扰动步长。

2.3 Boltzmann 玻尔兹曼退火

原文档公式（Boltzmann schedule）：

std = minimum(sqrt(T) * ones(d), (upper - lower) / (3*learn_rate))
y ~ Normal(0, std, size = d)
x_new = x_old + learn_rate * y

T_new = T0 / log(1 + k)

源码对应 SABoltzmann.get_new_x 与 SABoltzmann.cool_down：

def get_new_x(self, x):
    a, b = np.sqrt(self.T), self.hop / 3.0 / self.learn_rate
    std = np.where(a < b, a, b)          # 即 minimum(sqrt(T), hop/(3*learn_rate))
    xc = np.random.normal(0, 1.0, size=self.n_dim)
    x_new = x + xc * std * self.learn_rate
    if self.has_bounds:
        return np.clip(x_new, self.lb, self.ub)
    return x_new

def cool_down(self):
    self.T = self.T_max / np.log(self.iter_cycle + 1.0)

Boltzmann 变体使用正态分布扰动：扰动步长的标准差 std 取 sqrt(T) 与 hop/(3·learn_rate) 中的较小值（前者在高温时起主导、保证足够的探索幅度；后者在低温时限制扰动不超过可行域宽度的约三分之一）。其降温采用对数式 T0 / log(1 + k)，是所有变体中降温最慢的一种，因此在同等迭代次数下会保留更长时间的"高温探索"，适合对精度要求较高或对初值不敏感的场景。

2.4 三种变体小结

• 扰动来源：Fast 用符号指数式跳跃（有界大步长）、Cauchy 用重尾柯西采样、Boltzmann 用正态采样，三者分别对应"快速收敛""长尾穿越"与"精细探索"三种搜索风格；
• 降温速度：Boltzmann（对数）< Cauchy（倒数）< Fast（指数-幂），速度越快，算法越早进入局部细化阶段；
• 共性框架：三者共享 SimulatedAnnealingBase.run() 主循环，即每个温度下执行 L 次邻域搜索（链长），经 Metropolis 准则筛选后降温一轮，如此循环直到温度低于 T_min 或最优值连续 max_stay_counter 轮不变。

三、公共参数与运行机制（源码级解读）

无论选择哪种变体，构造与运行方式都一致。以 sko/SA.py 中 SimulatedAnnealingBase.__init__ 为准，核心参数如下：

参数	默认值	含义
func	必填	目标函数，输入为长度 n_dim 的数组，返回标量（求最小值）
x0	必填	初始解，数组长度为 n_dim，维度由 len(x0) 自动推断
T_max	100	初始温度，构造时强制满足 T_max > T_min > 0
T_min	1e-7	终止温度，当 self.T < self.T_min 时停止
L	300	每个温度下的迭代次数（链长，即内循环长度）
max_stay_counter	150	最优值连续不变的最大轮数，超过则提前停止（"冷却"机制）
lb / ub	None	变量下界/上界，两者要么同时给出、要么同时省略，否则抛出 ValueError
hop	见下文	扰动步长基准：有边界时默认 ub - lb，无边界时默认 10
learn_rate	0.5	仅 Cauchy 与 Boltzmann 变体使用，放缩扰动步长
m / n / quench	1	仅 Fast 变体使用，决定退火常数 c = m * exp(-n * quench)

运行主循环的核心逻辑（run() 方法）：

while True:
    for i in range(self.L):
        x_new = self.get_new_x(x_current)
        y_new = self.func(x_new)
        # Metropolis 准则
        df = y_new - y_current
        if df < 0 or np.exp(-df / self.T) > np.random.rand():
            x_current, y_current = x_new, y_new
            if y_new < self.best_y:
                self.best_x, self.best_y = x_new, y_new
    self.iter_cycle += 1
    self.cool_down()
    ...
    if self.T < self.T_min:
        break
    if stay_counter > self.max_stay_counter:
        break

要点解读：

• Metropolis 准则：df < 0（新解更优）必然接受；否则以概率 exp(-df / T) 接受更差的解。温度越高、恶化幅度越小，越容易接受差解——这正是跳出局部最优的机制；
• 双循环结构：外层循环负责降温，内层 L 次迭代在同一温度下反复采样；每轮结束记录 generation_best_Y、generation_best_X，并兼容历史属性 best_y_history、best_x_history；
• 双停止条件：温度低于 T_min，或最优值连续超过 max_stay_counter 轮未改善；run() 返回 (best_x, best_y) 元组；
• 结果属性：best_x / best_y 为历史最优解，best_y_history 为每轮最优值的轨迹，可用于绘制收敛曲线；
• 类还继承自 SkoBase（sko/base.py），其中 fit() 等价于 run()（但标记为即将弃用，建议统一使用 run()）。

四、Fast 模拟退火实战示例

原文档的示例代码取自 examples/demo_sa.py，目标函数为三维二次函数：

demo_func = lambda x: x[0] ** 2 + (x[1] - 0.05) ** 2 + x[2] ** 2

该函数理论最小值位于 (0, 0.05, 0) 附近，值约为 0，适合作为三类算法的对照基准。

4.1 基础版 Fast Simulated Annealing

from sko.SA import SAFast

sa_fast = SAFast(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_fast.run()
print('Fast Simulated Annealing: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)

说明：

• x0=[1, 1, 1] 同时确定了维度 n_dim=3；
• T_max=1, T_min=1e-9 指定了温度区间，由于 Fast 是指数降温，这样设置可以在较少轮内完成"高温探索 → 低温细化"；
• L=300 保证每个温度下充分采样；max_stay_counter=150 允许最优值最多 150 轮不变，避免无限循环；
• 如前面源码分析所述，示例中的 q=0.99 在 SAFast 中不会被读取（该变体读取的是 m/n/quench），因此该参数对结果无影响；实际需要调节 Fast 的退火曲线时，请使用 m、n、quench。

4.2 带边界约束的 Fast Simulated Annealing

from sko.SA import SAFast

sa_fast = SAFast(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
                 lb=[-1, 1, -1], ub=[2, 3, 4])
sa_fast.run()
print('Fast Simulated Annealing with bounds: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)

边界约束的处理逻辑位于 SimulatedAnnealingValue.__init__：

if lb is not None and ub is not None:
    self.has_bounds = True
    self.lb, self.ub = np.array(lb) * np.ones(self.n_dim), np.array(ub) * np.ones(self.n_dim)
    assert self.n_dim == len(self.lb) == len(self.ub), 'dim == len(lb) == len(ub) is not True'
    assert np.all(self.ub > self.lb), 'upper-bound must be greater than lower-bound'
    self.hop = kwargs.get('hop', self.ub - self.lb)

• 传入的 lb/ub 会被广播为长度 n_dim 的数组，维度必须与 x0 一致，且要求 ub > lb；
• 当给定边界时，hop 默认取 ub - lb——这正是原文档公式中 xc = y * (upper - lower) 的出处：扰动幅度按可行域宽度放缩；
• 每次生成新解后通过 np.clip(x_new, self.lb, self.ub) 将越界分量直接截断回边界内。

五、Boltzmann 模拟退火实战示例

5.1 基础版 Boltzmann Simulated Annealing

from sko.SA import SABoltzmann

sa_boltzmann = SABoltzmann(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_boltzmann.run()
print('Boltzmann Simulated Annealing: best_x is ', sa_boltzmann.best_x, 'best_y is ', sa_boltzmann.best_y)

Boltzmann 变体可额外传入 learn_rate（默认 0.5）来控制扰动步长的整体放缩：x_new = x_old + learn_rate * y。对数降温 T0 / log(1 + k) 意味着温度下降非常平缓，在同样的 T_min 条件下需要更多轮次，适合对解质量要求高、可接受更长运行时间的场景。同理，q=0.99 在 SABoltzmann 中同样不会被读取（该变体读取的是 learn_rate）。

5.2 带边界约束的 Boltzmann Simulated Annealing

from sko.SA import SABoltzmann

sa_boltzmann = SABoltzmann(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
                           lb=-1, ub=[2, 3, 4])
sa_boltzmann.run()
print('Boltzmann Simulated Annealing with bounds: best_x is ', sa_boltzmann.best_x, 'best_y is ', sa_boltzmann.best_y)

注意此处的边界写法：lb=-1 是标量，ub=[2, 3, 4] 是数组。得益于构造函数中的 np.array(lb) * np.ones(self.n_dim) 广播机制，标量与数组可以混用，标量会被自动扩展为所有维度的统一下界。Boltzmann 的扰动标准差会被 hop / (3 * learn_rate)（即 (upper - lower) / (3*learn_rate)）封顶，保证低温阶段不会越过可行域，因此配合边界使用时收敛行为更稳定。

六、Cauchy 模拟退火实战示例

6.1 基础版 Cauchy Simulated Annealing

from sko.SA import SACauchy

sa_cauchy = SACauchy(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_cauchy.run()
print('Cauchy Simulated Annealing: best_x is ', sa_cauchy.best_x, 'best_y is ', sa_cauchy.best_y)

Cauchy 变体同样支持 learn_rate（默认 0.5）参数。由于 tan(u) 的重尾特性，即使温度已经较低，仍有小概率产生大跨度跳跃，这对多峰、强非线性目标函数特别有价值——它比正态扰动更不容易困在局部极值中。其倒数降温 T0 / (1 + k) 介于 Fast 与 Boltzmann 之间。

6.2 带边界约束的 Cauchy Simulated Annealing

from sko.SA import SACauchy

sa_cauchy = SACauchy(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
                     lb=[-1, 1, -1], ub=[2, 3, 4])
sa_cauchy.run()
print('Cauchy Simulated Annealing with bounds: best_x is ', sa_cauchy.best_x, 'best_y is ', sa_cauchy.best_y)

与 Fast、Boltzmann 相同，Cauchy 在给定边界时同样通过 np.clip 将新解限制在 [lb, ub] 内，因此即使 tan(u) 偶尔产生极端值，也不会让搜索跑出可行域。

七、结果分析：如何评估与可视化收敛

run() 返回 (best_x, best_y)，同时对象保存了完整的历史轨迹。利用这些属性可以直观对比三种变体的收敛行为：

import matplotlib.pyplot as plt
import pandas as pd

# 以任意一个 sa_xxx 对象为例，绘制每轮最优值的最小值累计曲线
plt.plot(pd.DataFrame(sa_fast.best_y_history).cummin(axis=0))
plt.show()

• best_y_history：每轮外层循环结束后历史最优值构成的序列，cummin 后即"最优值随轮次下降"的单调收敛曲线（原文档与 examples/demo_sa.py 均采用这一画法）；
• 横向对比 SAFast、SABoltzmann、SACauchy 三条曲线，可以看到 Boltzmann 下降平缓、Fast 前期下降快、Cauchy 介于两者之间，由此可在实践中根据"时间预算 vs 解质量"选择变体；
• 迭代轮次由停止条件决定（达到 T_min 或连续 max_stay_counter 轮无改善），因此无需手动指定迭代数，适合黑盒调用。

八、补充：默认别名与更多 SA 能力

• SA = SAFast：在 sko/SA.py 中，无后缀的 SA 即 Fast 变体别名，from sko.SA import SA 的入门示例与 SAFast 完全等价；
• SA_TSP：同一模块还提供了面向旅行商问题（TSP）的 SA_TSP 类（使用 swap/reverse/transpose 三类置换变异生成新路线，见 sko/operators/mutation.py），完整示例在 examples/demo_sa_tsp.py，本文重点的三种连续优化变体与其共用 run() 主循环与停止逻辑；
• 安装方式：pip install scikit-opt；本地运行示例前请确保已安装 numpy、scipy，绘图示例还需 matplotlib 与 pandas。

九、小结：如何为你的问题选择合适的变体

场景	推荐变体	理由
快速出结果、目标函数平滑单峰	SAFast	指数降温收敛最快，所需轮次最少
多峰、强非线性、易陷入局部最优	SACauchy	重尾扰动持续提供越障能力，倒数降温适中
对解精度要求高、预算充足	SABoltzmann	对数降温最平缓，后期细化最充分
需要限定变量范围	任一变体 + lb/ub	统一的 np.clip 截断 + hop = ub - lb 自适应步长

实践建议：先用默认参数跑通 SAFast 建立基线，再按需切换 SABoltzmann/SACauchy 并微调 L（链长）、T_max/T_min（温度区间）与 max_stay_counter（冷却阈值）。如需进一步研究参数细节，可查阅 docs/en/args.md 与 docs/en/README.md，并在 examples/demo_sa.py 中直接修改参数进行对照实验。