模拟退火算法的 C++ 实现

最近的一个项目中需要实现个路径规划的算法，需要求得的路径的总长度尽可能的短。这就是典型的旅行商(TSP)问题了。解决这个问题的一个比较好用的方法就是模拟退火算法。网上关于用模拟退火算法解决 TSP 问题的文章挺多的，其中也有不少号称给出了 C++ 代码。但是说句实话，这些代码中没有一个是按照面向对象的思想来实现了，并没有把模拟退火算法的框架封装好。因此，遇到一个新问题时，改写算法就很麻烦。

看了几个网上的实现后，发现 gsl 库给出的框架还是不错的。但是 gsl 是 C 语言的库，虽然是按照面向对象的思想来设计的，但是实现的接口真是挺繁琐的，用起来并不好用。因此，我就花了半天时间，把 gsl 里的相关代码扒了出来，重新封装到了一个类中。

模拟退火算法最早的思想是由 Metropolis 等人提出的。1983 年, Kirkpatrick 成功地将退火思想引入到组合优化领域。简单的说，模拟退火算法就是一种基于 Monte-Carlo 迭代求解策略的随机寻优算法。关于模拟退火算法的理论网上有很多文章介绍。这里就不多写了，今天主要说说如何C++ 实现。

模拟退火算法的核心功能封装到了 SimulatedAnnualingSolver 类中。

#ifndef SIMULATEDANNUALINGSOLVER_H#define SIMULATEDANNUALINGSOLVER_H#include<random>class SimulatedAnnualingSolver{public:    SimulatedAnnualingSolver(int rand_seed);    void setParameters(double t_initial, double mu_t, double k = 1, double t_min = 1e-5, int iters_fixed_T = 100, double step_size = 1);    void enablePrint(bool on) {m_print_position = on;}    void solve();protected:    virtual double take_step(double step_size) = 0;    virtual void undo_step() = 0;    virtual void save_best() = 0;    virtual double energy() = 0;    virtual void print() = 0;private:    bool m_print_position;    int m_iters_fixed_T;    double m_step_size;    double m_k;    double m_t_initial;    double m_mu_t;    double m_t_min;    std::mt19937 m_randGenerator;    double boltzmann(double E, double new_E, double T, double k);};#endif // SIMULATEDANNUALINGSOLVER_H

#include "SimulatedAnnualingSolver.h"#include <cmath>#include <iostream>#define GSL_LOG_DBL_MIN   (-7.0839641853226408e+02)SimulatedAnnualingSolver::SimulatedAnnualingSolver(int rand_seed)    :m_t_initial(1),      m_mu_t(1.01),      m_k(1),      m_t_min(0.01),      m_iters_fixed_T(100),      m_print_position(false),      m_step_size(1),      m_randGenerator(rand_seed){}void SimulatedAnnualingSolver::setParameters(double t_initial, double mu_t, double k, double t_min, int iters_fixed_T, double step_size){    m_t_initial = t_initial;    m_mu_t = mu_t;    m_k = k;    m_t_min = t_min;    m_iters_fixed_T = iters_fixed_T;    m_step_size = step_size;}void SimulatedAnnualingSolver::solve(){    int n_evals = 1, n_iter = 0;    double E = energy();    double best_E = E;    save_best(); // 将当前状态存为最佳解    double T = m_t_initial;    double T_factor = 1.0 / m_mu_t;    std::uniform_real_distribution<> dis(0, 1);    while (1)    {        int n_accepts = 0;        int n_rejects = 0;        int n_eless = 0;        for (int i = 0; i < m_iters_fixed_T; ++i)        {            double new_E = take_step (m_step_size);            if(new_E <= best_E)            {                best_E = new_E;                save_best();            }            ++n_evals;            if (new_E < E)            {                if (new_E < best_E)                {                    best_E = new_E;                    save_best();                }                /* yay! take a step */                E = new_E;                ++n_eless;            }            else if (dis(m_randGenerator) < boltzmann(E, new_E, T, m_k))            {                /* yay! take a step */                E = new_E;                ++n_accepts;            }            else            {                undo_step(); // 回退到上一个状态                ++n_rejects;            }        }        if (m_print_position)        {            printf ("%5d   %7d  %12g", n_iter, n_evals, T);            print();            printf ("  %12g  %12g\n", E, best_E);        }        /* apply the cooling schedule to the temperature */        T *= T_factor;        ++n_iter;        if (T < m_t_min)        {            break;        }    }}inline double SimulatedAnnualingSolver::boltzmann(double E, double new_E, double T, double k){    double x = -(new_E - E) / (k * T);    /* avoid underflow errors for large uphill steps */    return (x < GSL_LOG_DBL_MIN) ? 0.0 : exp(x);}

这个类有几个纯虚函数。

• take_step 用来在解空间中随机游走到一个新的状态。
• undo_step 撤销最近的一次游走。
• save_best 将当前的状态保存下来，最为当前最佳的解。
• energy 能量函数，返回当前状态的目标函数值。（模拟退火算法就是求这个能量函数最小值所对应的状态）
• print 输出当前状态。这个函数只是起辅助作用，如果不需要，可以实现为空函数。

对于我们一个具体的问题，我们只需要定义一个派生自这个类的新类，在新类中实现这些虚函数。

比如对于我的问题，一个简化了的旅行商问题，我建立了一个 TspSolver 的类。

#ifndef TSPSOLVER_H#define TSPSOLVER_H#include "SimulatedAnnualingSolver.h"#include <QDebug>#include <QVector>#include <QPoint>#include <random>class tspSolver: public SimulatedAnnualingSolver{public:    tspSolver(int rand_seed, QVector<QPoint> points);    QVector<QPoint> result();protected:    double take_step(double step_size) override;    void undo_step() override;    void save_best() override;    double energy() override;    void print() override;private:    std::mt19937 m_randGen;    QVector<QPoint> m_points;    QVector<QPoint> m_best;    int m_N;    int m_n1;    int m_n2;};#endif // TSPSOLVER_H

#include "tspsolver.h"tspSolver::tspSolver(int rand_seed, QVector<QPoint> points)    :SimulatedAnnualingSolver(rand_seed){     std::random_device r;     std::seed_seq seed2{r(), r(), r(), r(), r(), r(), r(), r()};    m_randGen.seed(seed2);    m_points = points;    m_N = m_points.size();}void tspSolver::save_best(){    m_best = m_points;}QVector<QPoint> tspSolver::result(){    return m_best;}void tspSolver::print(){    qDebug() << m_points;}double tspSolver::take_step(double step_size){    std::uniform_int_distribution<> dis(1, m_N - 1);    m_n1 = dis(m_randGen);    do    {        m_n2 = dis(m_randGen);    }    while(m_n2 == m_n1);//    qDebug() << "n1 = " << m_n1 << "n2 = " << m_n2;    std::swap(m_points[m_n1], m_points[m_n2]);    return energy();}void tspSolver::undo_step(){    std::swap(m_points[m_n1], m_points[m_n2]);}double tspSolver::energy(){    double e = 0;    for(int i = 1; i < m_N; i ++)    {        QPoint a1 = m_points.at(i - 1);        QPoint a2 = m_points.at(i);        double x = abs(a1.x() - a2.x());        double y = abs(a1.y() - a2.y());        e += qMax(x, y);    }    return e;}

大家可以看到，我这个 TSP 问题中的两点间的距离不是普通的直线距离，而是 X Y 方向投影距离的最大值。
take_step 函数实现的也很简答，就是随机选两个点，交换这两个点的坐标。undo_step 就是把这两个点的坐标再交换回去。

下面是一个算例，随便选了 9 个点。规划出的路径看起来还不错。

这个算例选了 32 个点。规划出的路径挺复杂，但是每一步的长度都是 1，所以是最优解。当然这个问题的最优解不止一个。这里得到的解是最优解中较复杂的一个情形。