[AI]模拟退火解决TSP问题（含源码）

模拟退火是一种通用概率算法，用来在固定时间内寻求在一个大的搜寻空间内找到的最优解。模拟退火是S.
Kirkpatrick, C. D. Gelatt和M. P. Vecchi在1983年所发明。而V. Černý在1985年也独立发明此算法。（摘自Wiki百科：https://zh.wikipedia.org/wiki/%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%AB）

本次以144个城市为例解决TSP问题。

一、算法描述：

Step1： 随机产生一组解s并计算其代价，并用其初始化Ans和BestAns，初始化起始温度temp=TEMP，连续DecMaintain个温度最优解没有得到更新，DecMaintain初始为0，

Step2： 如果temp低于TEMPTHRE或连续DecMaintain=DECTHRE个温度最优值没有的到更新，算法结束；

Step3： 记录下当前最佳解LastAns，初始化马尔科夫链长IterCnt=ITER*CITYNUM，最优解回带次数计数器RepeatCnt=0，温度不变下最优解连续未被更新计数器MaintainCnt=0，进入马尔科夫链循环；

Step4： 如果尝试过的次数达到马尔科夫链长或最优解回带次数RepeatCnt耗尽，跳至Step9：；

Step5： 随机对现有解s中城市进行交换，并计算其代价，如果新的代价小于Ans，更新Ans，跳至Step7：；否则，转Step6：；

Step6： 根据接受函数计算对于s.cost的接受概率ProbAccp，并随机一个0~1之间的小数MyChoose，如果MyChoose< ProbAccp，则接受这个新解s，并更新至Ans；否则将s更新回Ans；

Step7： 如果得到的新Ans.cost比BestAns.cost小，则更新BestAns,并将MaintainCnt清零，转Step4：；否则，转Step8：；

Step8： 将计数器MaintainCnt加1，如果更新后达到MAINTAINTHRE，将BestAns回带到s，并将回带次数计数器RepeatCnt加1，MaintainCnt清0；转Step4：；

Step9： 如果LastAns与当前BestAns一致，则将DecMaintain加1；否则将其清0；

Step10： 对当前温度进行降温，转Step2：；

二、优化：

a) 对于同一温度下，采取“变长”马尔科夫链:

调试过程中发现降温到一定程度后，对于给定温度，执行14400次无法得到一个更优的解，并且对于一个差解也从未接受过，这就产生了一些无效工作，通过多次调试得到两个阈值MAINTAINTHRE= (ITER*CITYNUM)*0.25，REPEATTHRE = 3。当有1/4的循环次数都没有得到更好的解后，将BestAns回带到s中，给它一定产生更优解的辅助操作，但这种回带操作有一定的次数限制，这里设为3（最起码要保证0.25*3<1，不然并没有对马尔科夫链有实际缩短效果）。

b) 对不同温度下，观察其最优解更新情况并选择性提前停止降温：

该优化操作需要提前对总降温次数有一个总体估计DECTIME = log(double(TEMPTHRE /TEMP)) / log(COEF)，这里提前停止降温的条件是连续DECTIME*10%次降温后BestAns都没有得到更新，则默认它将无法产生更优的解，算法提前结束。

c) 起始温度TEMP、截止温度TEMPTHRE及降温系数COEF的设定：

对于起始温度，设的低了没有足够大概率去接受差解，设的高了浪费时间，书上给出280℃作为起始温度，但即使我将其置为1000℃、2800℃，实质上相比280℃并没有增加很多操作，如b)中计算总降温次数，在截至温度相同（0.0001）下，不同起始温度降温次数如表 1所示

表 1 起始温度与降温次数

起始温度（℃）	280	1000	2800
降温次数	734	797	848

这是由于数值上看似差距很大，但进行降温次数计算时为对数运算，实质上数量级没有发生很大变化，所以适当提高起始温度并不会产生很大的效率上的影响。

除此之外，起始温度也可根据Metropolis准则去生成，基本思想是使得在起始温度t0下各个状态的接受概率趋近于1（等概率出现），P0为接受概率，可取0.9或0.8。



对于截至温度设置可以类比起始温度，由于我采用了选择性提前停止降温，所以我在这里尽可能的将截至温度设低，目的是——如果你有能力得到更优解，那我给你足够低温度让你进行。实际上大多数（可以说基本所有）操作在0.1（甚至1）℃左右时就不会产生更优解了，将提前停止降温。

对于降温系数COEF，它在一定程度上与差解的接受函数有关，调试过程中发现在15℃~3℃范围内解得到优化的效果最好，为了避免降温过快没有得到足够稳定以及降温过慢效率底下，0.98是在多次调试过程中得到的经验系数。

d) 随机交换函数Swap(Solu &s):

Swap函数中随机出两个点leftpos、rightpos，采用三种交换方式，每次随机选择其中一种：

i. 直接交换两点；

ii. 将两点间的城市序列反转；

iii. 将leftpos放在rightpos之后。

以上操作全部基于随机选择，因为在算法设计中对于大数据量操作时，随机一般是效率最高的，产生较优选择的概率也是最高的。

e) 距离表Cost：

该表在读入数据时建立，提前计算好两点间的距离，避免之后的重复计算。

三、源码及测试文件：

#include <iostream>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <algorithm>
#define CITYNUM 144//城市数
#define TEMP 1000//初始温度
#define TEMPTHRE 0.0001//截至温度
#define ITER 100//马尔科夫链长
#define COEF 0.98//降温系数
#define PRECISION 1000000//产生随机整数后转化为小数时的精度
#define TESTTIME 10//测试次数
#define INT_MAX 2147483647
using namespace std;
struct MyCity{
int ID;
int x, y;
MyCity(int id = 0, int xx = 0, int yy = 0){
ID = id;
x = xx;
y = yy;
}
};
struct Solu
{
double cost;
int path[CITYNUM];
};
struct my_unique
{
int val;
my_unique(){ val = -1; }
int operator()(){ return ++val; }
}UniqueNum;
MyCity City[CITYNUM];//读入的城市数据
Solu Ans = { INT_MAX, { 0 } };//SA过程中产生的解
Solu BestAns = { INT_MAX, { 0 } };//SA总过程中产生的最优解
double Cost[CITYNUM][CITYNUM];//每两个城市间的代价，读入数据后提前算好，提高效率
double MyAns[TESTTIME];//多次执行算法得到的结果数组
const int DECTIME = log(double(TEMPTHRE / TEMP)) / log(COEF);//得到大致的降温次数
const int DECTHRE = DECTIME*0.1;//当有连续10%次降温没有得到更优解时提前结束算法
const int MAINTAINTHRE = (ITER*CITYNUM)*0.25;//某个温度下最优解连续不更新的次数阈值
const int REPEATTHRE = 3;/*与MAINTAINTHRE的系数相乘小于1，
为某个温度下将最优解重新带回退火过程的机会，可以提高效率*/
void OutputPath(Solu &Ans)
{
for (int i = 0; i < CITYNUM; i++)
printf("%d->", Ans.path[i]+1);
printf("%d\n", Ans.path[0]+1);
}
void Eval(Solu &s)
{
s.cost = 0;
for (int i = 1; i < CITYNUM; i++)
s.cost += Cost[s.path[i - 1]][s.path[i]];
s.cost += Cost[s.path[CITYNUM - 1]][s.path[0]];
}
void Swap(Solu &s)
{
int leftpos = rand() % CITYNUM;
int rightpos = rand() % CITYNUM;
if (leftpos == rightpos)
return;
if (leftpos > rightpos)
{
leftpos ^= rightpos;
rightpos ^= leftpos;
leftpos ^= rightpos;
}
int type = rand() % 3;
switch (type)
{
case 0:
reverse(s.path + leftpos, s.path + rightpos);
break;
case 1:
s.path[leftpos] ^= s.path[rightpos];
s.path[rightpos] ^= s.path[leftpos];
s.path[leftpos] ^= s.path[rightpos];
break;
case 2://rightpos前移
{
if (rightpos == CITYNUM + 1)
return;
rotate(s.path + leftpos, s.path + rightpos, s.path + rightpos + 1);
}
default:
break;
}
Eval(s);
}
void SimulatedAnnealing(Solu &s)
{
memcpy(&BestAns, &s, sizeof(Solu));
memcpy(&Ans, &s, sizeof(Solu));
double temp = TEMP;
int DecMaintain = 0;
while (temp > TEMPTHRE && DecMaintain < DECTHRE)
{
Solu LastAns;
memcpy(&LastAns, &BestAns, sizeof(Solu));
int IterCnt = ITER*CITYNUM;
int MaintainCnt = 0, RepeatCnt = 0;
while (IterCnt--&&RepeatCnt<REPEATTHRE)
{
Swap(s);
double distin = Ans.cost - s.cost;
if (distin < 0)
{
double ProbAccp = exp(distin / temp);
double MyChoose = double(rand() % (PRECISION + 1)) / PRECISION;
if (MyChoose > ProbAccp)//Not Accept
memcpy(&s, &Ans, sizeof(Solu));
else
memcpy(&Ans, &s, sizeof(Solu));
}
else
memcpy(&Ans, &s, sizeof(Solu));
if (Ans.cost < BestAns.cost)
{
memcpy(&BestAns, &Ans, sizeof(Solu));
MaintainCnt = 0;
}
else
{
MaintainCnt++;
if (MaintainCnt >= MAINTAINTHRE)
{
RepeatCnt++;
MaintainCnt = 0;
memcpy(&s, &BestAns, sizeof(Solu));
}
}
}
if (LastAns.cost == BestAns.cost)
DecMaintain++;
else if (LastAns.cost > BestAns.cost)
DecMaintain = 0;
temp *= COEF;
}
}
int main()
{
freopen("Cities(144).txt", "r", stdin);
//freopen("SAResult(50).txt", "w", stdout);
srand((unsigned)time(NULL));
for (int i = 0; i < CITYNUM; i++)
cin >> City[i].ID >> City[i].x >> City[i].y;
for (int i = 0; i < CITYNUM; i++)
{
Cost[i][i] = INT_MAX;
for (int j = i + 1; j < CITYNUM; j++)
Cost[i][j] = Cost[j][i] = sqrt(pow(City[i].x - City[j].x, 2) + pow(City[i].y - City[j].y, 2));
}
int test_time = TESTTIME;
memset(MyAns, 0, sizeof(MyAns));
while (test_time--)
{
cout << "*************************************************************************" << endl;
Solu TmpAns;
generate(TmpAns.path, TmpAns.path + CITYNUM, UniqueNum);
random_shuffle(TmpAns.path,TmpAns.path+CITYNUM);
Eval(TmpAns);
cout << "the SA result is:" << endl;
SimulatedAnnealing(TmpAns);
OutputPath(Ans);
MyAns[test_time] = BestAns.cost;
cout << Ans.cost << endl;
cout << "the best result is:" << endl;
OutputPath(BestAns);
cout << BestAns.cost << endl;
}
sort(MyAns, MyAns + TESTTIME);
for (int i = 0; i < TESTTIME; i++)
cout << MyAns[i] << endl;
fclose(stdin);
//fclose(stdout);
return 0;
}

1 3639 1315
2 4177 2244
3 3712 1399
4 3569 1438
5 3757 1187
6 3493 1696
7 3904 1289
8 3488 1535
9 3791 1339
10 3506 1221
11 3374 1750
12 3376 1306
13 3237 1764
14 3326 1556
15 3188 1881
16 3089 1251
17 3258 911
18 3814 261
19 3238 1229
20 3646 234
21 3583 864
22 4172 1125
23 4089 1387
24 4297 1218
25 4020 1142
26 4196 1044
27 4116 1187
28 4095 626
29 4312 790
30 4252 882
31 4403 1022
32 4685 830
33 4386 570
34 4361 73
35 4720 557
36 4643 404
37 4634 654
38 4153 426
39 4784 279
40 2846 1951
41 2831 2099
42 3007 1970
43 3054 1710
44 3086 1516
45 1828 1210
46 2562 1756
47 2716 1924
48 2061 1277
49 2291 1403
50 2751 1559
51 2788 1491
52 2012 1552
53 1779 1626
54 2381 1676
55 682  825
56 1478 267
57 1777 892
58 518  1251
59 278  890
60 1064 284
61 1332 695
62 3715 1678
63 3688 1818
64 4016 1715
65 4181 1574
66 3896 1656
67 4087 1546
68 3929 1892
69 3918 2179
70 4062 2220
71 3751 1945
72 3972 2136
73 4061 2370
74 4207 2533
75 4029 2498
76 4201 2397
77 4139 2615
78 3766 2364
79 3777 2095
80 3780 2212
81 3896 2443
82 3888 2261
83 3594 2900
84 3796 2499
85 3678 2463
86 3676 2578
87 3478 2705
88 3789 2620
89 4029 2838
90 3810 2969
91 3862 2839
92 3928 3029
93 4167 3206
94 4263 2931
95 4186 3037
96 3486 1755
97 3492 1901
98 3322 1916
99 3334 2107
100 3479 2198
101 3429 1908
102 3587 2417
103 3318 2408
104 3176 2150
105 3507 2376
106 3296 2217
107 3229 2367
108 3264 2551
109 3394 2643
110 3402 2912
111 3360 2792
112 3101 2721
113 3402 2510
114 3439 3201
115 3792 3156
116 3468 3018
117 3526 3263
118 3142 3421
119 3356 3212
120 3012 3394
121 3130 2973
122 3044 3081
123 2935 3240
124 2765 3321
125 3140 3550
126 3053 3739
127 2545 2357
128 2769 2492
129 2284 2803
130 2611 2275
131 2348 2652
132 2577 2574
133 2860 2862
134 2778 2826
135 2592 2820
136 2801 2700
137 2126 2896
138 2401 3164
139 2370 2975
140 1890 3033
141 1304 2312
142 1084 2313
143 3538 3298
144 3470 3304