模拟退火算法（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: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库)。
