模拟退火算法

模拟退火算法

​ 模拟退火算法最早的思想由Metropolis 等（ 1953 ）提出， 1983 年 Kirkpatrick 等将其应用于组合优化。

算法的目的：

克服优化过程陷入局部极小；

克服初值依赖性。

物理退火过程

在物理学中，退火是将金属加热到极高温度，然后让其极其缓慢地冷却的过程。

• 高温状态：原子运动剧烈，处于无序状态（高能量）。
• 等温过程 对于与环境换热而温度不变的封闭系统，系统状态的自发变化总是朝自由能减少的方向进行，当自由能达到最小时，系统达到平衡态；
• 缓慢冷却：原子逐渐找到最稳定的位置，形成整齐的晶体结构。
• 最终状态：系统的内能最低（全局最优）。

如果冷却得太快（淬火，Quenching），原子来不及调整位置就被“冻结”在杂乱的状态，系统处于亚稳态（局部最优，内能较高，材料脆）。

温度越低，物体的能量状态越低，到达足够的低点时，液体开始冷凝与结晶，在结晶状态时，系统的能量状态最低。缓慢降温（退火时，可达到最低能量状态；但如果快速降温（淬火，会导致不是最低能态的非晶形。

Boltzmann概率分布

在温度T，分子停留在状态r满足Boltzmann概率分布

在同一个温度T，选定两个能量E1<E2，有

（1）在同一个温度，分子停留在能量小状态的概率大于停留在能量大状态的概率

（2）温度越高，不同能量状态对应的概率相差越小；温度足够高时，各状态对应概率基本相同。

（3）随着温度的下降，能量最低状态对应概率越来越大；温度趋于0时，其状态趋于1

Metropolis准则

以概率接受新状态

假设当前状态为 x，能量为 E(x)。

我们随机生成了一个邻域新解 x'，能量为 E(x')。

定义能量差为：

简单来说，它的核心数学思想是：以一定的概率接受一个“更差”的解，从而跳出局部最优陷阱，最终趋向全局最优。

情况 A：新解更好 (ΔE<0)

如果新解的能量更低（比如在下山），我们100% 接受这个新解。

P(accept)=1

这对应了贪心算法（Gradient Descent）的部分。

情况 B：新解更差 (ΔE>0)

如果新解的能量更高（比如要爬坡，反方向），我们不是直接拒绝，而是以一定的概率接受它。这个概率 P 服从 玻尔兹曼分布 (Boltzmann Distribution)：

• ΔE：能量差（肯定为正）。

• T：当前的温度。

• k：物理中的玻尔兹曼常数（在算法中通常设为 1）。

(1) 温度 T 极高时 (探索阶段)

接受概率接近 100%。哪怕新解比旧解差很多，算法也会接受。

算法在搜索空间中随机游走 (Random Walk)，像个醉汉。这保证了它能翻越极高的山峰，从深坑里跳出来，遍历整个空间。

(2) 温度 T 降低时 (过渡阶段)

随着 T 变小，分母变小，指数部分变成较大的负数。

如果Delta E 很大（解变差很多），概率 P\(就会很小；如果 \Delta E\)很小（只差一点点），概率 P$还比较大。

算法开始变得挑剔。它仍然允许跳出浅坑（局部最优），但不再接受那些太离谱的差解。

(3) 温度 T 极低时 (收敛阶段)

接受更差解的概率几乎为 0。

算法退化为贪心算法 (Hill Climbing)。它只接受好解，不再爬坡。这时候它应该已经落入了全局最优

降温系数

在每一个固定的温度 T下，算法进行多次迭代。这实际上是在生成一个马尔可夫链。 如果迭代次数足够多，系统会达到服从玻尔兹曼分布热平衡分布 (Stationary Distribution)。

当 T缓慢下降时，概率分布图会变得越来越尖（Peaked），大部分概率密度会集中在全局最小值的附近。

引入冷却系数：

数学上的最优降温 (对数冷却)： Geman 在 1984 年证明，如果降温速度足够慢，满足：

那么模拟退火以概率 1 收敛到全局最优解。 缺点：这个速度太慢了，慢到实际上无法使用（可能需要几百年）。

工程上的降温 (指数冷却)：

这是对收敛速度和求解质量的折衷。虽然理论上不保证 100% 找到全局最优，但在有限时间内能找到“足够好”的解。

Rosenbrock 函数验证

N 维 Rosenbrock 函数的通常定义如下：

\[f(\mathbf{x}) = \sum_{i=1}^{N-1} [100 (x_{i+1} - x_i^2)^2 + (1 - x_i)^2] \]

c++库代码如下

sa.hpp

#ifndef SA_SA_HPP
#define SA_SA_HPP#include "params.hpp"
#include "policies.hpp"
#include "detail/solver.hpp"namespace sa {/*** @brief 模拟退火通用求解函数* * @tparam State 状态类型 (自动推导)* @tparam EnergyFunc 能量函数类型 (自动推导)* @tparam NeighborFunc 邻域函数类型 (可选)* @tparam CoolingPolicy 降温策略 (可选)* @tparam ConstraintPolicy 约束策略 (可选)* * @param initial_state 初始状态值* @param energy_func 能量函数句柄* @param params 算法参数配置* @param neighbor 邻域生成器实例* @param cooling 降温器实例* @param constraint 约束器实例* @return std::pair<State, double> {最优状态, 最优能量值}*/template<typename State,typename EnergyFunc,typename NeighborFunc = DefaultNeighbor<State>,typename CoolingPolicy = ExponentialCooling,typename ConstraintPolicy = std::nullptr_t>auto solve(const State& initial_state,EnergyFunc energy_func,Params params = Params{},NeighborFunc neighbor = NeighborFunc{},CoolingPolicy cooling = CoolingPolicy{},ConstraintPolicy constraint = ConstraintPolicy{}) {using AcceptancePolicy = MetropolisAcceptance;detail::Solver<State, EnergyFunc, NeighborFunc, CoolingPolicy, AcceptancePolicy, ConstraintPolicy>solver(params, energy_func, neighbor, cooling, AcceptancePolicy{}, constraint);return solver.solve(initial_state);}} // namespace sa#endif // SA_SA_HPP

policies.hpp

#ifndef SA_POLICIES_HPP
#define SA_POLICIES_HPP#include "params.hpp"
#include "detail/traits.hpp"#include <cmath>
#include <random>
#include <algorithm>
#include <stdexcept>
#include <vector>namespace sa {// 默认降温策略struct ExponentialCooling {inline double operator()(double T, const Params& p) const noexcept {return T * p.alpha;}};// 默认接受策略 (Metropolis 准则)struct MetropolisAcceptance {template<typename RNG>bool operator()(double delta_E,double T,RNG& rng,std::uniform_real_distribution<double>& dist) const{if (delta_E < 0.0) return true;return std::exp(-delta_E / T) > dist(rng);}};// 默认连续邻域生成策略template<typename State>struct DefaultNeighbor {double sigma = 1.0;State operator()(const State& current,double T,std::mt19937& rng) const{if constexpr (std::is_arithmetic_v<State>) {std::normal_distribution<double> dist(0.0, sigma * T);return static_cast<State>(current + dist(rng));}else if constexpr (detail::is_std_vector_v<State>) {using ValueType = typename State::value_type;static_assert(std::is_arithmetic_v<ValueType>,"vector value type must be arithmetic");std::normal_distribution<double> dist(0.0, sigma * T);State candidate = current;for (auto& v : candidate)v = static_cast<ValueType>(v + dist(rng));return candidate;}else {static_assert(sizeof(State) == 0, "No default neighbor for this State type");return current;}}};// 离散翻转邻域策略 (针对 vector<bool> 或 vector<int>)template<typename State>struct DiscreteFlipNeighbor {State operator()(const State& current,double T,std::mt19937& rng) const{static_assert(detail::is_std_vector_v<State>, "DiscreteFlipNeighbor requires std::vector");using ValueType = typename State::value_type;static_assert(std::is_same_v<ValueType, int> || std::is_same_v<ValueType, bool>,"DiscreteFlipNeighbor requires vector<int> or vector<bool>");State candidate = current;std::uniform_int_distribution<std::size_t> dist(0, candidate.size() - 1);std::size_t idx = dist(rng);if constexpr (std::is_same_v<ValueType, bool>)candidate[idx] = !candidate[idx];elsecandidate[idx] = 1 - candidate[idx];return candidate;}};// 边界约束策略 (Box Constraint)template<typename State>class BoxConstraintPolicy {public:using ValueType = std::conditional_t<std::is_arithmetic_v<State>,State,typename State::value_type>;BoxConstraintPolicy(ValueType lower, ValueType upper): lower_(lower), upper_(upper) {}void apply(State& state) const noexcept {if constexpr (std::is_arithmetic_v<State>) {state = std::clamp(state, lower_, upper_);}else {for (auto& v : state)v = std::clamp(v, lower_, upper_);}}private:ValueType lower_;ValueType upper_;};} // namespace sa#endif // SA_POLICIES_HPP

params.hpp

//
// Created by 31007 on 2026/2/12.
//#ifndef MATH_TYPES_HPP
#define MATH_TYPES_HPP
#include <cstddef>
#include <cstdint>
namespace sa {struct Params {double      initial_temp     = 100.0;       // 初始温度double      final_temp       = 1e-6;        // 终止温度double      alpha            = 0.98;        // 降温系数std::size_t iter_per_temp    = 100;         // 每个温度下的迭代次数std::size_t max_total_iters  = 1'000'000;   // 最大总迭代次数 (防止死循环)std::uint32_t seed           = 0;           // 随机种子 (0表示随机)};} // namespace sa
#endif //MATH_TYPES_HPP

solver.hpp

#ifndef SA_DETAIL_SOLVER_HPP
#define SA_DETAIL_SOLVER_HPP#include "../params.hpp"
#include <random>
#include <utility>
#include <stdexcept>namespace sa::detail {template<typename State,typename EnergyFunc,typename NeighborFunc,typename CoolingPolicy,typename AcceptancePolicy,typename ConstraintPolicy>class Solver {public:Solver(Params params,EnergyFunc energy,NeighborFunc neighbor,CoolingPolicy cooling,AcceptancePolicy acceptance,ConstraintPolicy constraint): params_(params),energy_(energy),neighbor_(neighbor),cooling_(cooling),acceptance_(acceptance),constraint_(constraint),dist_(0.0, 1.0){validate_params();if (params_.seed == 0) {std::random_device rd;rng_ = std::mt19937(rd());} else {rng_ = std::mt19937(params_.seed);}}std::pair<State, double> solve(const State& initial_state) {State current = initial_state;double current_energy = energy_(current);State best = current;double best_energy = current_energy;double T = params_.initial_temp;std::size_t total_iters = 0;while (T > params_.final_temp && total_iters < params_.max_total_iters) {for (std::size_t i = 0;i < params_.iter_per_temp && total_iters < params_.max_total_iters;++i, ++total_iters){State candidate = neighbor_(current, T, rng_);// 编译期判断是否存在约束策略if constexpr (!std::is_same_v<ConstraintPolicy, std::nullptr_t>) {constraint_.apply(candidate);}double candidate_energy = energy_(candidate);double delta = candidate_energy - current_energy;if (acceptance_(delta, T, rng_, dist_)) {current = std::move(candidate);current_energy = candidate_energy;if (current_energy < best_energy) {best = current;best_energy = current_energy;}}}T = cooling_(T, params_);}return {std::move(best), best_energy};}private:void validate_params() {if (params_.initial_temp <= 0.0) throw std::invalid_argument("initial_temp must be > 0");if (params_.final_temp <= 0.0) throw std::invalid_argument("final_temp must be > 0");if (params_.alpha <= 0.0 || params_.alpha >= 1.0) throw std::invalid_argument("alpha must be in (0,1)");if (params_.iter_per_temp == 0) throw std::invalid_argument("iter_per_temp must be > 0");}Params params_;EnergyFunc energy_;NeighborFunc neighbor_;CoolingPolicy cooling_;AcceptancePolicy acceptance_;ConstraintPolicy constraint_;std::mt19937 rng_;std::uniform_real_distribution<double> dist_;};} // namespace sa::detail#endif // SA_DETAIL_SOLVER_HPP

traits.hpp

#ifndef SA_DETAIL_TRAITS_HPP
#define SA_DETAIL_TRAITS_HPP#include <vector>
#include <type_traits>namespace sa::detail {// 类型萃取：判断是否为 std::vectortemplate<typename T>struct is_std_vector : std::false_type {};template<typename T, typename Alloc>struct is_std_vector<std::vector<T, Alloc>> : std::true_type {};template<typename T>inline constexpr bool is_std_vector_v = is_std_vector<T>::value;} // namespace sa::detail#endif // SA_DETAIL_TRAITS_HPP