模擬退火演算法（SA）求解TSP 問題（C語言實現）

　　　　這篇文章是之前寫的智慧演算法（遺傳演算法（GA）、粒子群演算法（PSO））的補充。其實程式碼我老早之前就寫完了，今天恰好重新翻到了，就拿出來給大家分享一下，也當是回顧與總結了。

　　　　首先介紹一下模擬退火演算法（SA）。模擬退火演算法（simulated annealing，SA）演算法最早是由Metropolis等人提出的。其出發點是基於物理中固體物質的退火過程與一般組合優化問題之間的相似性。模擬退火演算法是一種通用的優化演算法，其物理退火過程由以下三部分組成：

　　　　（1）加溫過程

　　　　（2）等溫過程

　　　　（3）冷卻過程

　　　　　其中加溫過程對應演算法設定的初始溫度，等溫過程對應演算法的Metropolis抽樣過程，冷卻過程對應控制引數的下降。這裡能量的變化就是目標函式，要得到的最優解就是能量最低狀態。Metropolis準則是SA演算法收斂於全域性最優解的關鍵所在，Metropolis準則以一定的概率接受惡化解，這樣就使得演算法可以跳離區域性最優解的陷阱。

　　　　　模擬退火演算法為求解傳統方法難以處理的TSP問題提供了一個有效的途徑和通用的處理框架，並逐漸發展成為一種迭代自適應啟發式概率搜尋演算法。模擬退火演算法可以用於求解不同的非線性問題，對於不可微甚至不連續函式的優化，能以較大概率求得全域性最優解，該演算法還具有較強的魯棒性、全域性收斂性、隱含並行性以及廣泛的適應性，對目標函式以及約束函式沒有任何要求。

　　　　　SA 演算法實現的步驟如下：（下面以最小化問題為例）

　　　　　（1）初始化：取溫度T0足夠大，令T = T0，取任意解S1，確定每個T時 的迭代次數，即 Metropolis鏈長L。

　　　　　（2）對當前溫度T和k=1,2,3，...,L,重複步驟（3）~（6）

　　　　　（3）對當前解S1隨機產生一個擾動得到一個新解 S2.

　　　　　（4） 計算S2的增量df = f（S2） - f（S1），其中f（S1）為S1的代價函式。

　　　　　 （5）若df < 0,接受S2作為新的當前解，即S1 = S2；否則S2的接受概率為 exp（-df/T）,即隨機產生（0,1）上的均勻分佈的隨機數rand，若 exp（-df/T）>rand

,則接受S2作為新的當前解，S1 = S2；否則保留當前解。

                （6）如果滿足最終的終止條件，Stop，則輸出當前解S1作為最優解，結束程式。終止條件Stop通常為：在連續若干個Metropolis鏈中新解S2都沒有被接受時終止演算法，或者是設定結束溫度。否則按衰減函式衰減T後返回（2）

　　　　　　以上的步驟稱之為Metropolis過程。逐步降低控制溫度，重複Metropolis過程，直至滿足結束準則Stop求出最優解。可以看出SA整體的步驟相比GA以及PSO還是簡單了很多了，而且親測效果還不錯，所以屬於價效比較高的演算法。關鍵的步驟在第（6）步。

原始碼如下：

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<string.h>
#include<math.h>

#define T0 50000.0    //初始溫度
#define T_end (1e-8)
#define q 0.98   //退火係數
#define L 1000    //每個溫度時的迭代次數，即鏈長
#define N 31     //城市數量
int city_list[N];  //用於存放一個解
double city_pos[N][2] = {{1304,2312},{3639,1315},{4177,2244},{3712,1399},{3488,1535},{3326,1556},{3238,1229},{4196,1004},{4312,790}, 
{4386,570},{3007,1970},{2562,1756},{2788,1491},{2381,1676},{1332,695},{3715,1678},{3918,2179},{4061,2370},{3780,2212},{3676,2578},{4029,2838},
{4263,2931},{3429,1908},{3507,2367},{3394,2643},{3439,3201},{2935,3240},{3140,3550},{2545,2357},{2778,2826},{2370,2975}}; 
 // 中國31個城市座標

//函式宣告
double distance(double *,double*);   //計算兩個城市的距離
double path_len(int *);       //計算路徑長度
void  init();  //初始化函式
void create_new(); // 產生新解

//距離函式
double distance(double *city1,double *city2)
{

     double x1=*city1;
     double y1=*(city1+1);
     double x2=*city2;
     double y2=*(city2+1);
     double dis=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
     return dis;

}

//計算路徑長度
double path_len(int *arr)
{
  
    double path=0;  //初始化路徑長度
    int index=*arr;  //定位到第一個數字（城市序號）
    for(int i=0;i<N;i++)
    {
        int index1=*(arr+i);
        int index2=*(arr+i+1);
        double dis=distance(city_pos[index1-1],city_pos[index2-1]);
        path+=path+dis;
    }
    int last_index=*(arr+N-1);   //最後一個城市序號
    int first_index=*arr;         //第一個城市序號
    double last_dis=distance(city_pos[last_index-1],city_pos[first_index-1]);
    path=path+last_dis;
    return path;  //返回路徑總長度
}

//初始化函式
void init()
{
   for(int i=0;i<N;i++)
       city_list[i]=i+1;   //初始化一個解
}

//產生一個新解， 此處採用隨機交叉兩個位置的方式產生新的解
void create_new()
{
     double r1=((double)rand())/(RAND_MAX+1.0);
     double r2=((double)rand())/(RAND_MAX+1.0);
     int pos1=(int)(N*r1);   //第一個交叉點的位置
     int pos2=(int)(N*r2);
     int temp=city_list[pos1];
     city_list[pos1]=city_list[pos2];
     city_list[pos2]=temp;      //交換兩個點
}

//主函式
int main()
{
    srand((unsigned)time(NULL));   //初始化隨機數種子
    time_t start,finish;
    start=clock();  //程式執行開始計時
    double T;
    int count=0;
    T=T0;  //初始溫度
    init(); //初始化一個解
    int city_list_copy[N]; //用於儲存原始解
    double f1,f2,df; //f1為初始解目標函式值，f2為新解目標函式值，df為二者差值
    double r;   //0-1之間的隨機數，用來決定是否接受新解
    while(T>T_end)   //當溫度低於結束溫度時，退火結束
    {
        for(int i=0;i<L;i++)
        {
            memcpy(city_list_copy,city_list,N*sizeof(int)); //複製陣列
            create_new();
            f1=path_len(city_list_copy);
            f2=path_len(city_list);
            df=f2-f1;
            //以下是metropolis準則
            if(df>=0)
            {
                r=((double)rand())/(RAND_MAX);
                if(exp(-df/T)<=r)  //保留原來的解
                {
                     memcpy(city_list,city_list_copy,N*sizeof(int));
                }
            }
        }
        T*=q;  //降溫
        count++;
    }
    finish=clock();  //退火過程結束
    double duration=((double)(finish-start))/CLOCKS_PER_SEC;  //計算時間
    printf("採用模擬退火演算法，初始溫度T0=%.2f,降溫係數q=%.2f,每個溫度迭代%d次,共降溫%d次，得到的TSP最優路徑為:\n",T0,q,L,count);
    for(int i=0;i<N-1;i++)  //輸出最優路徑
    {
        printf("%d--->",city_list[i]);
    }
    printf("%d\n",city_list[N-1]);
    double len = path_len(city_list); // 最優路徑長度
    printf("最優路徑長度為:%lf\n",len);
    printf("程式執行耗時:%lf秒.\n",duration);
    return 0;
}
 

結果如下截圖：