模拟退火算法

模拟退火算法​模拟退火算法最早的思想由Metropolis等1953提出1983年Kirkpatrick等将其应用于组合优化。算法的目的克服优化过程陷入局部极小克服初值依赖性。物理退火过程在物理学中退火是将金属加热到极高温度然后让其极其缓慢地冷却的过程。高温状态原子运动剧烈处于无序状态高能量。等温过程对于与环境换热而温度不变的封闭系统系统状态的自发变化总是朝自由能减少的方向进行当自由能达到最小时系统达到平衡态缓慢冷却原子逐渐找到最稳定的位置形成整齐的晶体结构。最终状态系统的内能最低全局最优。如果冷却得太快淬火Quenching原子来不及调整位置就被“冻结”在杂乱的状态系统处于亚稳态局部最优内能较高材料脆。温度越低物体的能量状态越低到达足够的低点时液体开始冷凝与结晶在结晶状态时系统的能量状态最低。缓慢降温退火时可达到最低能量状态但如果快速降温淬火会导致不是最低能态的非晶形。Boltzmann概率分布在温度T分子停留在状态r满足Boltzmann概率分布在同一个温度T选定两个能量E1E2有1在同一个温度分子停留在能量小状态的概率大于停留在能量大状态的概率2温度越高不同能量状态对应的概率相差越小温度足够高时各状态对应概率基本相同。3随着温度的下降能量最低状态对应概率越来越大温度趋于0时其状态趋于1Metropolis准则以概率接受新状态假设当前状态为 x能量为 E(x)。我们随机生成了一个邻域新解 x能量为 E(x)。定义能量差为简单来说它的核心数学思想是以一定的概率接受一个“更差”的解从而跳出局部最优陷阱最终趋向全局最优。情况 A新解更好 (ΔE0)如果新解的能量更低比如在下山我们100% 接受这个新解。P(accept)1这对应了贪心算法Gradient Descent的部分。情况 B新解更差 (ΔE0)如果新解的能量更高比如要爬坡反方向我们不是直接拒绝而是以一定的概率接受它。这个概率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_{i1}^{N-1} [100 (x_{i1} - x_i^2)^2 (1 - x_i)^2] \]c库代码如下sa.hpp/* by 01130.hk - online tools website : 01130.hk/zh/guid.html */ #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::pairState, double {最优状态, 最优能量值} */ template typename State, typename EnergyFunc, typename NeighborFunc DefaultNeighborState, 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::SolverState, EnergyFunc, NeighborFunc, CoolingPolicy, AcceptancePolicy, ConstraintPolicy solver(params, energy_func, neighbor, cooling, AcceptancePolicy{}, constraint); return solver.solve(initial_state); } } // namespace sa #endif // SA_SA_HPPpolicies.hpp/* by 01130.hk - online tools website : 01130.hk/zh/guid.html */ #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 { templatetypename RNG bool operator()(double delta_E, double T, RNG rng, std::uniform_real_distributiondouble dist) const { if (delta_E 0.0) return true; return std::exp(-delta_E / T) dist(rng); } }; // 默认连续邻域生成策略 templatetypename State struct DefaultNeighbor { double sigma 1.0; State operator()(const State current, double T, std::mt19937 rng) const { if constexpr (std::is_arithmetic_vState) { std::normal_distributiondouble dist(0.0, sigma * T); return static_castState(current dist(rng)); } else if constexpr (detail::is_std_vector_vState) { using ValueType typename State::value_type; static_assert(std::is_arithmetic_vValueType, vector value type must be arithmetic); std::normal_distributiondouble dist(0.0, sigma * T); State candidate current; for (auto v : candidate) v static_castValueType(v dist(rng)); return candidate; } else { static_assert(sizeof(State) 0, No default neighbor for this State type); return current; } } }; // 离散翻转邻域策略 (针对 vectorbool 或 vectorint) templatetypename State struct DiscreteFlipNeighbor { State operator()(const State current, double T, std::mt19937 rng) const { static_assert(detail::is_std_vector_vState, DiscreteFlipNeighbor requires std::vector); using ValueType typename State::value_type; static_assert( std::is_same_vValueType, int || std::is_same_vValueType, bool, DiscreteFlipNeighbor requires vectorint or vectorbool ); State candidate current; std::uniform_int_distributionstd::size_t dist(0, candidate.size() - 1); std::size_t idx dist(rng); if constexpr (std::is_same_vValueType, bool) candidate[idx] !candidate[idx]; else candidate[idx] 1 - candidate[idx]; return candidate; } }; // 边界约束策略 (Box Constraint) templatetypename State class BoxConstraintPolicy { public: using ValueType std::conditional_t std::is_arithmetic_vState, 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_vState) { 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_HPPparams.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 1000000; // 最大总迭代次数 (防止死循环) std::uint32_t seed 0; // 随机种子 (0表示随机) }; } // namespace sa #endif //MATH_TYPES_HPPsolver.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::pairState, 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_vConstraintPolicy, 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_distributiondouble dist_; }; } // namespace sa::detail #endif // SA_DETAIL_SOLVER_HPPtraits.hpp#ifndef SA_DETAIL_TRAITS_HPP #define SA_DETAIL_TRAITS_HPP #include vector #include type_traits namespace sa::detail { // 类型萃取判断是否为 std::vector templatetypename T struct is_std_vector : std::false_type {}; templatetypename T, typename Alloc struct is_std_vectorstd::vectorT, Alloc : std::true_type {}; templatetypename T inline constexpr bool is_std_vector_v is_std_vectorT::value; } // namespace sa::detail #endif // SA_DETAIL_TRAITS_HPP