模拟退火算法

模拟退火算法

​ 模拟退火算法最早的思想由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];
            else
                candidate[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::vector

    template<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