模拟退⽕算法（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
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 // 城市数量
19int city_list[N]; // ⽤于存放⼀个解
20double 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//函数声明
24double distance(double *,double *); // 计算两个城市距离
25double path_len(int *); // 计算路径长度
26void init(); //初始化函数
27void create_new(); // 产⽣新解
28// 距离函数
29double distance(double * city1,double * city2)
30 {
31double x1 = *city1;
32double y1 = *(city1+1);
33double x2 = *(city2);
34double y2 = *(city2+1);
35double dis = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
36return dis;
37 }
38
39// 计算路径长度
40double path_len(int * arr)
41 {
42double path = 0; // 初始化路径长度
43int index = *arr; // 定位到第⼀个数字(城市序号)
44for(int i=0;i<N-1;i++)
45 {
46int index1 = *(arr+i);
47int index2 = *(arr+i+1);
48double dis = distance(city_pos[index1-1],city_pos[index2-1]);
49 path += dis;
50 }
51int last_index = *(arr+N-1); // 最后⼀个城市序号
52int first_index = *arr; // 第⼀个城市序号
53double last_dis = distance(city_pos[last_index-1],city_pos[first_index-1]);
54 path = path + last_dis;
55return path; // 返回总的路径长度
56 }
57
58// 初始化函数
59void init()
60 {
61for(int i=0;i<N;i++)
62 city_list[i] = i+1; // 初始化⼀个解
63 }
64
65// 产⽣⼀个新解
66// 此处采⽤随机交叉两个位置的⽅式产⽣新的解
67void create_new()
68 {
69double r1 = ((double)rand())/(RAND_MAX+1.0);
70double r2 = ((double)rand())/(RAND_MAX+1.0);
71int pos1 = (int)(N*r1); //第⼀个交叉点的位置
72int pos2 = (int)(N*r2);
73int temp = city_list[pos1];
74 city_list[pos1] = city_list[pos2];
75 city_list[pos2] = temp; // 交换两个点
76 }
77
78// 主函数
79int main(void)
80 {
81 srand((unsigned)time(NULL)); //初始化随机数种⼦
82 time_t start,finish;
83 start = clock(); // 程序运⾏开始计时
84double T;
85int count = 0; // 记录降温次数
86 T = T0; //初始温度
87 init(); //初始化⼀个解
88int city_list_copy[N]; // ⽤于保存原始解
89double f1,f2,df; //f1为初始解⽬标函数值，f2为新解⽬标函数值，df为⼆者差值 90double r; // 0-1之间的随机数，⽤来决定是否接受新解
91while(T > T_end) // 当温度低于结束温度时，退⽕结束
92 {
93for(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准则
101if(df >= 0)
102 {
103 r = ((double)rand())/(RAND_MAX);
104if(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(); // 退⽕过程结束
114double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 计算时间
115 printf("采⽤模拟退⽕算法，初始温度T0=%.2f，降温系数q=%.2f，每个温度迭代%d次，共降温%d次，得到的TSP最优路径为:\n",T0,q,L,count); 116for(int i=0;i<N-1;i++) // 输出最优路径
117 {
118 printf("%d--->",city_list[i]);
119 }
120 printf("%d\n",city_list[N-1]);
121double len = path_len(city_list); // 最优路径长度
122 printf("最优路径长度为:%lf\n",len);
明解c语言123 printf("程序运⾏耗时:%lf秒.\n",duration);
124return0;
125 }
运⾏结果截图如下：
注：如果使⽤gcc编译报错，可能需要加上 -std=c99选项以及在末尾加上 -lm(链接math库)。
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论