模拟退火算法（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）步。

　　　　　　不廢話了，直接上代碼吧。TSP的資料和之前的資料一樣，使用C語言實作。代碼如下：

　　　　　　

1 /*
  2  * 使用模拟退火算法(SA)求解TSP問題(以中國TSP問題為例)
  3  * 參考自《Matlab 智能算法30個案例分析》
  4  * 模拟退火的原理這裡略去，可以參考上書或者相關論文
  5  * update: 16/12/11
  6  * author:lyrichu
  7  * email:[email protected]
  8  */
  9 #include<stdio.h>
 10 #include<stdlib.h>
 11 #include<string.h>
 12 #include<time.h>
 13 #include<math.h>
 14 #define T0 50000.0  // 初始溫度
 15 #define T_end (1e-8)
 16 #define q  0.98   // 退火系數
 17 #define L 1000  // 每個溫度時的疊代次數，即鍊長
 18 #define N 31  // 城市數量
 19 int city_list[N]; // 用于存放一個解
 20 double city_pos[N][2] = {{1304,2312},{3639,1315},{4177,2244},{3712,1399},{3488,1535},{3326,1556},{3238,1229},{4196,1004},{4312,790},
 21     {4386,570},{3007,1970},{2562,1756},{2788,1491},{2381,1676},{1332,695},{3715,1678},{3918,2179},{4061,2370},{3780,2212},{3676,2578},{4029,2838},
 22     {4263,2931},{3429,1908},{3507,2367},{3394,2643},{3439,3201},{2935,3240},{3140,3550},{2545,2357},{2778,2826},{2370,2975}}; // 中國31個城市坐标
 23 //函數聲明
 24 double distance(double *,double *); // 計算兩個城市距離
 25 double path_len(int *);  // 計算路徑長度
 26 void  init();  //初始化函數
 27 void create_new(); // 産生新解
 28 // 距離函數
 29 double distance(double * city1,double * city2)
 30 {
 31     double x1 = *city1;
 32     double y1 = *(city1+1);
 33     double x2 = *(city2);
 34     double y2 = *(city2+1);
 35     double dis = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 36     return dis;
 37 }
 38 
 39 // 計算路徑長度
 40 double path_len(int * arr)
 41 {
 42     double path = 0; // 初始化路徑長度
 43     int index = *arr; // 定位到第一個數字(城市序号)
 44     for(int i=0;i<N-1;i++)
 45     {
 46         int index1 = *(arr+i);
 47         int index2 = *(arr+i+1);
 48         double dis = distance(city_pos[index1-1],city_pos[index2-1]);
 49         path += dis;
 50     }
 51     int last_index = *(arr+N-1); // 最後一個城市序号
 52     int first_index = *arr; // 第一個城市序号
 53     double last_dis = distance(city_pos[last_index-1],city_pos[first_index-1]);
 54     path = path + last_dis;
 55     return path; // 傳回總的路徑長度
 56 }
 57 
 58 // 初始化函數
 59 void init()
 60 {
 61     for(int i=0;i<N;i++)
 62         city_list[i] = i+1;  // 初始化一個解
 63 }
 64 
 65 // 産生一個新解
 66 // 此處采用随機交叉兩個位置的方式産生新的解
 67 void create_new()
 68 {
 69     double r1 = ((double)rand())/(RAND_MAX+1.0);
 70     double r2 = ((double)rand())/(RAND_MAX+1.0);
 71     int pos1 = (int)(N*r1); //第一個交叉點的位置
 72     int pos2 = (int)(N*r2);
 73     int temp = city_list[pos1];
 74     city_list[pos1] = city_list[pos2];
 75     city_list[pos2] = temp;   // 交換兩個點
 76 }
 77 
 78 // 主函數
 79 int main(void)
 80 {
 81     srand((unsigned)time(NULL)); //初始化随機數種子
 82     time_t start,finish;
 83     start = clock(); // 程式運作開始計時
 84     double T;
 85     int count = 0; // 記錄降溫次數
 86     T = T0; //初始溫度
 87     init(); //初始化一個解
 88     int city_list_copy[N]; // 用于儲存原始解
 89     double f1,f2,df; //f1為初始解目标函數值，f2為新解目标函數值，df為二者內插補點
 90     double r; // 0-1之間的随機數，用來決定是否接受新解
 91     while(T > T_end) // 當溫度低于結束溫度時，退火結束
 92     {
 93         for(int i=0;i<L;i++)
 94         {
 95             memcpy(city_list_copy,city_list,N*sizeof(int)); // 複制數組
 96             create_new(); // 産生新解
 97             f1 = path_len(city_list_copy);
 98             f2 = path_len(city_list);
 99             df = f2 - f1;
100             // 以下是Metropolis準則
101             if(df >= 0)
102             {
103                 r = ((double)rand())/(RAND_MAX);
104                 if(exp(-df/T) <= r) // 保留原來的解
105                 {
106                     memcpy(city_list,city_list_copy,N*sizeof(int));
107                 }
108             }
109         }
110         T *= q; // 降溫
111         count++;
112     }
113     finish = clock(); // 退火過程結束
114     double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 計算時間
115     printf("采用模拟退火算法，初始溫度T0=%.2f,降溫系數q=%.2f,每個溫度疊代%d次,共降溫%d次，得到的TSP最優路徑為:\n",T0,q,L,count);
116     for(int i=0;i<N-1;i++)  // 輸出最優路徑
117     {
118         printf("%d--->",city_list[i]);
119     }
120     printf("%d\n",city_list[N-1]);
121     double len = path_len(city_list); // 最優路徑長度
122     printf("最優路徑長度為:%lf\n",len);
123     printf("程式運作耗時:%lf秒.\n",duration);
124     return 0;
125 }      

運作結果截圖如下：

注：如果使用gcc編譯報錯，可能需要加上 -std=c99選項以及在末尾加上 -lm(連結math庫)。

熱愛程式設計，熱愛機器學習！

github:http://www.github.com/Lyrichu

github blog:http://Lyrichu.github.io

個人部落格站點:http://www.movieb2b.com(不再維護)