c++實現的模擬退火演算法

演算法簡介

  模擬退火演算法得益於材料的統計力學的研究成果。統計力學表明材料中粒子的不同結構對應於粒子的不同能量水平。在高溫條件下，粒子的能量較高，可以自由運動和重新排列。在低溫條件下，粒子能量較低。如果從高溫開始，非常緩慢地降溫(這個過程被稱為退火)，粒子就可以在每個溫度下達到熱平衡。當系統完全被冷卻時，最終形成處於低能狀態的晶體。

  假定我們要解決的問題是一個尋找最小值的優化問題。將物理學中模擬退火的思想應用於優化問題就可以得到模擬退火尋優方法。

---

演算法過程

  考慮這樣一個組合優化問題：我們現在要求函式f(x)的最小值，其中x∈Sx∈S，S表示函式的定義域，N(x)⊆SN(x)⊆S表示x的一個鄰域集合。

  首先給定一個初始溫度T0和該優化問題的一個初始解x(0)，並由x(0)生成下一個解x′∈N(x(0))x′∈N(x(0)),是否接受x′x′作為新解x(1)依賴於下面概率： 

P(x(0)→x′)={1若f(x′)<f(x(0))e−f(x′)−f(x(0))T0其它P(x(0)→x′)={1若f(x′)<f(x(0))e−f(x′)−f(x(0))T0其它

換句話說，如果生成的解x′x′的函式值比前一個解的函式值更小，則接受x(1)=x′x(1)=x′作為一個新解。否則以概率e−f(x′)−f(x(0))T0e−f(x′)−f(x(0))T0接受x′x′作為一個新解。

 

  同樣的，對於模擬退火過程中任意一個解x(k)，接受x′x′作為下一個新解x(k+1)的概率和上面公式類似： 

P(x(k)→x′)={1若f(x′)<f(x(k))e−f(x′)−f(x(k))Ti其它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

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