0.写在前面

  本文要使用的知识包括C语言、一点STL和一点点数据结构知识就可以,实现起来并不困难。

1.旅行商问题

  旅行商问题,即TSP问题(Travelling Salesman Problem)又译为旅行推销员问题、货郎担问题,是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市,他必须选择所要走的路径,路径的限制是每个城市只能拜访一次,而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。
  旅行商问题是一个 NP 完全问题,目前求解 TSP 问题的主要方法有模拟退火算法、遗传算法、启发式搜索法、蚁群算法等。

2.模拟退火算法

  模拟退火算法(Simulated Annealing)来源于固体退火原理,是一种基于概率的算法。将固体加温至充分高的温度,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,分子和原子越不稳定。而徐徐冷却时粒子渐趋有序,能量减少,原子越稳定。在冷却(降温)过程中,固体在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。
  模拟退火算法从某一较高初温出发,伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解,即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性,从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。
  模拟退火算法包含两个部分即Metropolis算法和退火过程,,分别对应内循环和外循环。外循环就是退火过程,将固体达到较高的温度(初始温度T0),然后按照降温系数alpha使温度按照一定的比例下降,当达到终止温度Tf时,冷却结束,即退火过程结束。
  Metropolis算法是内循环,即在每次温度下,迭代L次,寻找在该温度下能量的最小值(即最优解)。下图中所示即为在一次温度下,迭代L次,能量发生的变化。在该温度下,整个迭代过程中温度不发生变化,能量发生变化,当前一个状态x(n)的能量大于后一个状态x(n+1)的能量时,状态x(n)的解没有状态x(n+1)的解好,所以接受状态x(n+1)。但是如果下一状态的能量比前一个状态的能量高时,那就按照Metropolis准则对新解按照一定概率接受,这就避免了爬山算法中陷入局部最优解。

全局最优与局部最优

Metropolis准则:

Metropolis准则

3.相关数据

软件:DevC++ 5.10
硬件:任意计算机都可以
数据集:TSPLIP上的一个小数据集

4.代码实现过程

需要使用的头文件:
1
2
3
4
5
6
7
8
9
10
11
#include<stdio.h>
#include<iostream>
#include<map>
#include<string>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<math.h>
#include<algorithm>
#include<vector>
using namespace std;

变量声明:

1
2
3
4
5
6
7
8
9
10
11
12
map<int,string> mp;  
const int N=48; //城市数量
const int MIN=0x3fffffff; //初始化最小值
const double T0=50000,T_end=1e-8; //初始温度和最终温度
double q=0.997; //降温系数
const int L=1000; //链长
int List[N]; //经过的城市序号顺序
int temp[N]; //临时变量,中途经过的城市序号顺序
double minp=MIN;
struct node{
double x,y;
}; //城市坐标结构体

本文使用ATT48@TSPLIB数据集,数据集如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
node City[48]= {
{6734,1453},
{2233,10},
{5530,1424},
{401,841},
{3082,1644},
{7608,4458},
{7573,3716},
{7265,1268},
{6898,1885},
{1112,2049},
{5468,2606},
{5989,2873},
{4706,2674},
{4612,2035},
{6347,2683},
{6107,669},
{7611,5184},
{7462,3590},
{7732,4723},
{5900,3561},
{4483,3369},
{6101,1110},
{5199,2182},
{1633,2809},
{4307,2322},
{675,1006},
{7555,4819},
{7541,3981},
{3177,756},
{7352,4506},
{7545,2801},
{3245,3305},
{6426,3173},
{4608,1198},
{23,2216},
{7248,3779},
{7762,4595},
{7392,2244},
{3484,2829},
{6271,2135},
{4985,140},
{1916,1569},
{7280,4899},
{7509,3239},
{10,2676},
{6807,2993},
{5185,3258},
{3023,1942}};

绘图如下:

ATT48

距离计算公式采用伪欧几里得公式:

ATT48

1
2
3
4
5
6
7
8
double Distance(int num1,int num2)
{
double x1=City[num1].x,x2=City[num2].x,y1=City[num1].y,y2=City[num2].y;
double dis=sqrt(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))/10);
int dis2=dis;
if(dis>dis2) dis2+=1;
return dis2;
}

计算路径总长度的函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
double Path(int t[])
{
double path = 0; // 初始化路径长度
for(int i=0;i<N-1;i++)
{
int index1 = t[i];
int index2 = t[i+1];
double dis = Distance(index1,index2);
path += dis;
}
int last_index =t[N-1]; // 最后一个城市序号
int first_index =t[0]; // 第一个城市序号
double last_dis = Distance(last_index,first_index);
path = path + last_dis;
return path; // 返回总的路径长度
}

初始化函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
void init()
{
srand((unsigned)time(NULL));
bool visit[N];
fill(visit,visit+N,false);
double rr=((double)rand())/(RAND_MAX+1.0);
temp[0]=rr*N;
visit[temp[0]]=true;
for(int i=1;i<N;i++)
{
int mind=MIN,minn=-1;
for(int j=0;j<N;j++)
{
if(visit[j]==false&&Distance(temp[i-1],j)<mind)
{
mind=Distance(temp[i-1],j);minn=j;
}
}
temp[i]=minn; // 初始化一个解
visit[minn]=true;
}
for(int i=0;i<N;i++) List[i]=temp[i];
minp=Path(temp);
}

产生新解的方案(一共有6种):

产生新解的方案

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
void create_new1()//方案一,随机找到两个点,逆转两点间的其他点 
{
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 first=min(pos1,pos2);
int end=max(pos1,pos2);
while(end>first)
{
int tempnum=temp[first];
temp[first]=temp[end];
temp[end]=tempnum;
first++;end--;
}
}
void create_new2()//方案二,随机交换两个点
{
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 tempnum=temp[pos1];
temp[pos1]=temp[pos2];
temp[pos2]=tempnum;
return ;
}
void create_new3()//方案三,随机生成一点,并与其右侧点交换,若为末尾,则与右侧点交换
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
if(pos1==N-1)
{
int tempnum=temp[pos1];
temp[pos1]=temp[pos1-1];
temp[pos1-1]=tempnum;
}
else
{
int tempnum=temp[pos1];
temp[pos1]=temp[pos1+1];
temp[pos1+1]=tempnum;
}
return;
}
void create_new4()//方案四,随机生成一点,并与其对称位置点交换
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
int pos2;
if(pos1<N/2) pos2=pos1+N/2;
else pos2=pos1-N/2;
// cout<<endl<<pos1<<" "<<pos2<<endl;
int tempnum=temp[pos1];
temp[pos1]=temp[pos2];
temp[pos2]=tempnum;
return ;
}
void create_new5()//方案五,反序前插排序
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
int templist[N];
memcpy(templist,temp,N*sizeof(int));
int pos2=N-1,temp1=0;
while(pos2>pos1)
{
temp[temp1++]=templist[pos2--];
}
pos2=0;
while(pos2<=pos1) temp[temp1++]=templist[pos2++];

return ;
}
void create_new6()//方案六,后向前插排序
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
double r2= ((double)rand())/(RAND_MAX+1.0);
int pos2 = (int)((N)*r2);
int first1=0,first2=min(pos1,pos2),first3=max(pos1,pos2),end1=first2-1,end2=first3-1,end3=N-1;
while(end1>first1)
{
int tempnum=temp[first1];
temp[first1]=temp[end1];
temp[end1]=tempnum;
first1++;end1--;
}
while(end2>first2)
{
int tempnum=temp[first2];
temp[first2]=temp[end2];
temp[end2]=tempnum;
first2++;end2--;
}
while(end3>first3)
{
int tempnum=temp[first3];
temp[first3]=temp[end3];
temp[end3]=tempnum;
first3++;end3--;
}
return ;
}

main函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
int main()
{

srand((unsigned)time(NULL)); //初始化随机数种子
time_t start,finish;
int chflag=0;
int tempcopy[N];
start = clock(); // 程序运行开始计时
double T;
int count = 0; // 记录降温次数
T = T0; //初始温度
init(); //初始化一个解
double f1,f2,df; //f1为初始解目标函数值,f2为新解目标函数值,df为二者差值
double r; // 0-1之间的随机数,用来决定是否接受新解
while(T>T_end) // 当温度低于结束温度时,退火结束
{ if(chflag>0)
{
chflag--;
if(chflag==0) q=0.997;
}
int cf=0;
for(int i=0;i<L;i++)
{
memcpy(tempcopy,temp,N*sizeof(int)); // 复制数组
int choosenum=rand()%2; // 产生新解,此处为1,5方案随机,可以通过修改取余的数字和if-else语句修改所需要的方案类型
// int choosenum=5;
if(choosenum==0) create_new1();
// else if(choosenum==1) create_new2();
// else if(choosenum==2) create_new3();
// else if(choosenum==3) create_new4();
else if(choosenum==1) create_new5();
// else if(choosenum==2) create_new6();

f1 = Path(tempcopy);
f2 = Path(temp);
df = f2 - f1;
if(f2<minp)//若新解代价小于最小值,直接接受
{
minp=f2;
memcpy(List,temp,N*sizeof(int));
cf=1;
}
// 根据Metropolis准则是否接受新解
if(df >= 0)
{
r = ((double)rand())/(RAND_MAX);
if(exp(-df/T) <= r) // 保留原来的解
{
memcpy(temp,tempcopy,N*sizeof(int));
}
}
}
if(cf==1)
{
chflag=700;q=0.99999;
}
//cout<<q<<endl;
T*=q; // 降温
printf("当前温度为:%.18f 当前最短路径长度为:%.0f\n",T,minp);
count++;

/* for(int i=0;i<N;i++) // 输出最优路径
{
printf("%d--->",temp[i]);
}
cout<<endl;*/
}
finish = clock(); // 退火过程结束
double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 计算时间
printf("采用模拟退火算法,初始温度T0=%.2f,降温系数q=%.3f,每个温度迭代%d次,共降温%d次,得到的TSP最优路径为:\n",T0,q,L,count);
for(int i=0;i<N-1;i++) // 输出最优路径
{
printf("%d--->",List[i]);
}
printf("%d\n",List[N-1]);
// double len = path_len(city_list); // 最优路径长度
printf("最优路径长度为:%lf\n",minp);
printf("程序运行耗时:%lf秒.\n",duration);
return 0;
}

5.结果

某次结果如图:

结果

使用Origin 2019b绘图,可以更直观的看到经过的路线:

结果可视化

在TSPLIB上显示该数据集的最优解为10628,该代码可以达到这个效果,不过由于退火的随机性,每次的结果都会有差异,如果没有获得最优解可以多尝试几次。

附录

完整代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
#include<stdio.h>
#include<iostream>
#include<map>
#include<string>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<math.h>
#include<algorithm>
#include<vector>
using namespace std;
map<int,string> mp;
const int N=48; //城市数量
const int MIN=0x3fffffff; //初始化最小值
const double T0=50000,T_end=1e-8; //初始温度和最终温度
double q=0.997; //降温系数
const int L=1000; //链长
int List[N]; //经过的城市序号顺序
int temp[N]; //临时变量,中途经过的城市序号顺序
double minp=MIN;
struct node{
double x,y;
}; //城市坐标结构体
node City[N]= {
{6734,1453},
{2233,10},
{5530,1424},
{401,841},
{3082,1644},
{7608,4458},
{7573,3716},
{7265,1268},
{6898,1885},
{1112,2049},
{5468,2606},
{5989,2873},
{4706,2674},
{4612,2035},
{6347,2683},
{6107,669},
{7611,5184},
{7462,3590},
{7732,4723},
{5900,3561},
{4483,3369},
{6101,1110},
{5199,2182},
{1633,2809},
{4307,2322},
{675,1006},
{7555,4819},
{7541,3981},
{3177,756},
{7352,4506},
{7545,2801},
{3245,3305},
{6426,3173},
{4608,1198},
{23,2216},
{7248,3779},
{7762,4595},
{7392,2244},
{3484,2829},
{6271,2135},
{4985,140},
{1916,1569},
{7280,4899},
{7509,3239},
{10,2676},
{6807,2993},
{5185,3258},
{3023,1942}};
//函数申明
double Distance(int num1,int num2);//计算两城市直接距离




double Distance(int num1,int num2)
{
double x1=City[num1].x,x2=City[num2].x,y1=City[num1].y,y2=City[num2].y;
double dis=sqrt(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))/10);
int dis2=dis;
if(dis>dis2) dis2+=1;
return dis2;
}

double Path(int t[])
{
double path = 0; // 初始化路径长度
for(int i=0;i<N-1;i++)
{
int index1 = t[i];
int index2 = t[i+1];
double dis = Distance(index1,index2);
path += dis;
}
int last_index =t[N-1]; // 最后一个城市序号
int first_index =t[0]; // 第一个城市序号
double last_dis = Distance(last_index,first_index);
path = path + last_dis;
return path; // 返回总的路径长度
}

// 初始化函数
void init()
{
srand((unsigned)time(NULL));
bool visit[N];
fill(visit,visit+N,false);
double rr=((double)rand())/(RAND_MAX+1.0);
temp[0]=rr*N;
visit[temp[0]]=true;
for(int i=1;i<N;i++)
{
int mind=MIN,minn=-1;
for(int j=0;j<N;j++)
{
if(visit[j]==false&&Distance(temp[i-1],j)<mind)
{
mind=Distance(temp[i-1],j);minn=j;
}
}
temp[i]=minn; // 初始化一个解
visit[minn]=true;
}
for(int i=0;i<N;i++) List[i]=temp[i];
minp=Path(temp);
}

// 产生一个新解

void create_new1()//方案一,随机找到两个点,逆转两点间的其他点
{
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 first=min(pos1,pos2);
int end=max(pos1,pos2);
while(end>first)
{
int tempnum=temp[first];
temp[first]=temp[end];
temp[end]=tempnum;
first++;end--;
}
}
void create_new2()//方案二,随机交换两个点
{
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 tempnum=temp[pos1];
temp[pos1]=temp[pos2];
temp[pos2]=tempnum;
return ;
}
void create_new3()//方案三,随机生成一点,并与其右侧点交换,若为末尾,则与右侧点交换
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
if(pos1==N-1)
{
int tempnum=temp[pos1];
temp[pos1]=temp[pos1-1];
temp[pos1-1]=tempnum;
}
else
{
int tempnum=temp[pos1];
temp[pos1]=temp[pos1+1];
temp[pos1+1]=tempnum;
}
return;
}
void create_new4()//方案四,随机生成一点,并与其对称位置点交换
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
int pos2;
if(pos1<N/2) pos2=pos1+N/2;
else pos2=pos1-N/2;
// cout<<endl<<pos1<<" "<<pos2<<endl;
int tempnum=temp[pos1];
temp[pos1]=temp[pos2];
temp[pos2]=tempnum;
return ;
}
void create_new5()//方案五,反序前插排序
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
int templist[N];
memcpy(templist,temp,N*sizeof(int));
int pos2=N-1,temp1=0;
while(pos2>pos1)
{
temp[temp1++]=templist[pos2--];
}
pos2=0;
while(pos2<=pos1) temp[temp1++]=templist[pos2++];

return ;
}
void create_new6()//方案六,后向前插排序
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)((N)*r1);
double r2= ((double)rand())/(RAND_MAX+1.0);
int pos2 = (int)((N)*r2);
int first1=0,first2=min(pos1,pos2),first3=max(pos1,pos2),end1=first2-1,end2=first3-1,end3=N-1;
while(end1>first1)
{
int tempnum=temp[first1];
temp[first1]=temp[end1];
temp[end1]=tempnum;
first1++;end1--;
}
while(end2>first2)
{
int tempnum=temp[first2];
temp[first2]=temp[end2];
temp[end2]=tempnum;
first2++;end2--;
}
while(end3>first3)
{
int tempnum=temp[first3];
temp[first3]=temp[end3];
temp[end3]=tempnum;
first3++;end3--;
}


return ;
}

int main()
{

srand((unsigned)time(NULL)); //初始化随机数种子
time_t start,finish;
int chflag=0;
int tempcopy[N];
start = clock(); // 程序运行开始计时
double T;
int count = 0; // 记录降温次数
T = T0; //初始温度
init(); //初始化一个解
double f1,f2,df; //f1为初始解目标函数值,f2为新解目标函数值,df为二者差值
double r; // 0-1之间的随机数,用来决定是否接受新解
while(T>T_end) // 当温度低于结束温度时,退火结束
{ if(chflag>0)
{
chflag--;
if(chflag==0) q=0.997;
}
int cf=0;
for(int i=0;i<L;i++)
{
memcpy(tempcopy,temp,N*sizeof(int)); // 复制数组
int choosenum=rand()%2; // 产生新解,此处为1,5方案随机,可以通过修改取余的数字和if-else语句修改所需要的方案类型
// int choosenum=5;
if(choosenum==0) create_new1();
// else if(choosenum==1) create_new2();
// else if(choosenum==2) create_new3();
// else if(choosenum==3) create_new4();
else if(choosenum==1) create_new5();
// else if(choosenum==2) create_new6();

f1 = Path(tempcopy);
f2 = Path(temp);
df = f2 - f1;
if(f2<minp)//若新解代价小于最小值,直接接受
{
minp=f2;
memcpy(List,temp,N*sizeof(int));
cf=1;
}
// 根据Metropolis准则是否接受新解
if(df >= 0)
{
r = ((double)rand())/(RAND_MAX);
if(exp(-df/T) <= r) // 保留原来的解
{
memcpy(temp,tempcopy,N*sizeof(int));
}
}
}
if(cf==1)
{
chflag=700;q=0.99999;
}
//cout<<q<<endl;
T*=q; // 降温
printf("当前温度为:%.18f 当前最短路径长度为:%.0f\n",T,minp);
count++;

/* for(int i=0;i<N;i++) // 输出最优路径
{
printf("%d--->",temp[i]);
}
cout<<endl;*/
}
finish = clock(); // 退火过程结束
double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 计算时间
printf("采用模拟退火算法,初始温度T0=%.2f,降温系数q=%.3f,每个温度迭代%d次,共降温%d次,得到的TSP最优路径为:\n",T0,q,L,count);
for(int i=0;i<N-1;i++) // 输出最优路径
{
printf("%d--->",List[i]);
}
printf("%d\n",List[N-1]);
// double len = path_len(city_list); // 最优路径长度
printf("最优路径长度为:%lf\n",minp);
printf("程序运行耗时:%lf秒.\n",duration);
return 0;
}