模拟退火算法（c++实现）

模拟退火算法

---

算法简介

  模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温(这个过程被称为退火)，粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

  假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。

---

算法过程

  考虑这样一个组合优化问题：我们现在要求函数f(x)的最小值，其中x∈S，S表示函数的定义域，N(x)⊆S表示x的一个邻域集合。

  首先给定一个初始温度T0和该优化问题的一个初始解x(0)，并由x(0)生成下一个解x′∈N(x(0)),是否接受x′作为新解x(1)依赖于下面概率：

P(x(0)→x′)=⎧⎩⎨⎪⎪1若f(x′)<f(x(0))e−f(x′)−f(x(0))T0其它

换句话说，如果生成的解x′的函数值比前一个解的函数值更小，则接受x(1)=x′作为一个新解。否则以概率e−f(x′)−f(x(0))T0接受x′作为一个新解。

  同样的，对于模拟退火过程中任意一个解x(k)，接受x′作为下一个新解x(k+1)的概率和上面公式类似：

P(x(k)→x′)=⎧⎩⎨⎪⎪1若f(x′)<f(x(k))e−f(x′)−f(x(k))Ti其它

  整个过程中，温度Ti不断下降。模拟退火有概率接受一个更差的解，目的是为了能够跳出局部最优，从而寻找全局最优解。在每个Ti下，所得到的一个新状态x(k+1)完全依赖于前一个状态x(k)，与前面的状态x(0),x(1)…,x(k-1)无关，因此是一个马尔科夫过程。

  这里还有一个结论：如果温度下降十分缓慢，而且在每个温度都有足够多次的状态转移，使其在每一个温度下达到热平衡，则全局最优解将以概率1被找到。

---

模拟退火解决TSP问题

  问题：给定平面上20个点的名称与坐标，两个点之间的距离为它们的欧几里得距离。求一条路径，刚好经过每个点1次，使其路径长度最短。

  输入文件：in.txt

20A 2.5 4.0B 1.2 -2.4C 8.7 1.2D 3.6 12.1E -5.5 0.94F -6.6 -12.6G 0.18 5.219H 12.5 14.3609I 22.5 -5.26J 1.61 4.5K 2.1 -5.6L 0 25M 9.2 -32N -1 7O -5 -8P 21 35Q 16 7.5R -21 5S -7 -25.5T 12 17.5

  以下是c++的实现代码：

#include <algorithm>#include <cmath>#include <cstdio>#include <cstring>#include <deque>#include <iostream>#include <map>#include <queue>#include <set>#include <stack>#include <string>#include <vector>using namespace std;typedef long long LL;const int maxn = 1e2 + 7;const int INF = 0x7fffffff;const double PI = acos(-1);struct Point { //点类    string name;    double x, y;    int i;  //编号};vector<Point> p;double d[maxn][maxn]; //距离矩阵int n;double sum=0; //当前最短路径长度double dist(Point a, Point b) { //计算两点距离    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));}double get_sum(vector<Point> a){  //返回路径长度    double sum=0;    for(int i=1;i<a.size();i++){        sum+=d[a[i].i][a[i-1].i];    }    sum += d[a[0].i][a[a.size()-1].i];    return sum;}void init() { //初始化    srand((unsigned)time(NULL)); //设置随机数种子    cin >> n;    p.clear();    for (int i = 0; i < n; i++) {        Point t;        cin >> t.name >>t.x >>t.y;        t.i=i;        p.push_back(t);    }    for (int i = 0; i < n; i++) {        for (int j = i + 1; j < n; j++) {            d[i][j] = d[j][i] = dist(p[i], p[j]);        }    }    sum=get_sum(p);}void Monte_Carlo(){  //蒙特卡洛得到一个较好的初始解    vector<Point> cur=p;    for(int t=0;t<8000;t++){        for(int i=0;i<n;i++){            int j = rand()%n;            swap(cur[i],cur[j]);        }        double new_sum = get_sum(cur);        if(new_sum<sum){            sum=new_sum;            p=cur;        }    }}double e=1e-16,at=0.99999999,T=1.0;//e表示终止温度  at表示温度的变化率  T是初始温度int L = 200000; //L为最大迭代次数void Simulated_Annealing(){  //模拟退火     while(L--){  //最多迭代L次        vector<Point> cur=p;        int c1=rand()%n,c2=rand()%n;        if(c1==c2){            L++;continue;        }        swap(cur[c1],cur[c2]);        double df = get_sum(cur) - sum;  //计算代价        double sj=rand()%10000;     //sj为0~1之间的随机数        sj/=10000;        if(df < 0){  //如果结果更优，直接接受            p = cur;            sum+=df;        }else if(exp(-df/T)>sj){  //否则以一定概率接受            p = cur;            sum+=df;        }        T*=at;  //降温        if(T<e) break;  //T降到终止温度后退出    }}void show(){  //显示当前结果    cout<<"路径长度: "<<sum<<endl;    cout<<"路径:";    for(int i=0;i<n;i++)        cout<<' '<<p[i].name;    puts("");}int main() {#ifndef ONLINE_JUDGE    freopen("in.txt", "r", stdin);#endif    init();    Monte_Carlo();    Simulated_Annealing();    show();    return 0;}

  运行后得到结果（每次运行结果可能不同）：

路径长度: 221.983路径: I M S F O R E B K C A J G N D L P T H Q