這篇文章是之前寫的智能算法（遺傳算法（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:919987476@qq.com 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庫)。