用模拟退火算法解旅行商问题

1.问题描述

旅行商问题（Travelling Salesman Problem, 简记TSP，亦称货郎担问题)：设 有n个城市和距离矩阵D=[dij]，其中dij表示 城市i到城市j的距离，i，j=1，2 … n，则问 题是要找出遍访每个城市恰好一次的一条 回路并使其路径长度为最短。

2.算法设计

旅行商问题是一个十分经典的NP难度问题，如果想找到真正的唯一最优的解复杂度是O(N!)的，所以求解这一类问题的策略就是找一个相对最优的解，也就是最优化问题。模拟退火算法就是一种启发式的组合优算法，通过不断迭代来寻找最优解。
旅行商问题的解可以表述为一个循环排列，我们可以把求解旅行商问题当做把n个城市的全排列进行对比，求最短的一条路径，但是这样复杂度太高，这里使用一个叫Mapkob链的东西，也就是城市的最排列次数不能超过Mapkob链长度（L = 20000）。对于当前已有的排列，产生新排列的过程使用2变换法，也就是在当前的队列中任选两个序号u和v（u<v），将u和v及其之间的顺序逆转。在逆转之后需要重新计算新路径长度，若新路径长度小于当前路径，就用新路径代替当前路径，否则使用模拟退火算法根据概率决定是否转移。

数据结构：

struct Path{   //路径，包括路径和路径长度
    int city[MaxCitylenth];
    double len;
};
double d[MaxCitylenth][MaxCitylenth];  //城市之间的距离
struct City{  //每个城市坐标
    double x;
    double y;
};   

1、新解产生：
新解的产生使用2变换法，即通过任选两个序号u和v（u<v），将u和v及其之间的顺序逆转：
（π1 …πu-1πuπu+1 …πv-1πvπv+1 … πn） 变为
（π1 …πu-1πvπv-1 …πu-1πuπv+1 … πn）
在逆转之后需要重新计算路径长度，但是通过观察，其实在u之前和v之后的长度已经计算过了，为了节约计算复杂度，这里我先将路径长度减去u和v之间的长度，然后加上逆转后u-v的路径长度能得到新的路径长度。
2、新解更新：
如果新解的路径长度小于当前解，就用新解替换当前解，否则使用模拟退火算法计算概率，通过概率确定是否接受这个新解。公式如下：

其中t我初试设置为100，每次更新是ti+1 = a*ti，其中a是温度下降控制参数，在代码里为0.9。对每个t计算L次路径，当对于临近的两个t的值，路径没有变化时，此时找到的路径就是这次优化的结果，退出程序。图1描述了新解产生与更新过程。

图1 新解产生更新流程图

###3.程序流程
1、按顺序1-n的顺序初始化path；
2、将s置为0，L置为20000，从0到L循环，由当前路径生成新路径；
3、如图一更新当前解，若新解路径长度小于当前解，就用新解替换当前解，否则用模拟退火算法判断是否接受新解；
4、若循环L次后解都没被替换，则s++，否则s = 0;
5、当s为2的时候，也就是连续两个Mapkob链没有变化时，退出程序，当前找到解就是最优解。

4.核心伪代码

Path GenerateNewpath(){  //生成新路径，使用2变换法
    Path newpath = path;
    int u = 0,v = 0;
    while(u == v){
        u = (int)(n * (rand() / (RAND_MAX + 1.0)));
        v = (int)(n * (rand() / (RAND_MAX + 1.0)));
    }
    int x = min(u,v);
    int y = max(u,v);
    if(x == 0 || y  == n-1){ //若x是起点或y是终点，则要减去终点到起点的路径长度
        newpath.len -= d[newpath.city[n-1]][newpath.city[0]];
    }
    if(x != 0){
        newpath.len -= d[newpath.city[x-1]][newpath.city[x]];
    }
    for(int i = x;i <= y;i++){  //减去u-v的路径长度
        if(i < n-1)  //如果v不是终点，则要减去v到v+1的长度
            newpath.len -= d[newpath.city[i]][newpath.city[i+1]];
    }
    for(int i = x,j = y;i < j;i++,j--){
        swap(newpath.city[i],newpath.city[j]);
    }
    if(x == 0||y == n-1){
        newpath.len += d[newpath.city[n-1]][newpath.city[0]];
    }
    if(x != 0){
        newpath.len += d[newpath.city[x-1]][newpath.city[x]];
    }
    for(int i = x;i <= y;i++){
        if(i < n-1){
            newpath.len += d[newpath.city[i]][newpath.city[i+1]];
        }
    }
    return newpath;
}
bool Accept(int l,double t){ //用模拟退火算法判断是否接受目前路径
    double rd = rand()/(RAND_MAX + 1.0);
    double e = exp(l/t);
    if(e > rd && e < 1){
        return true;
    }
    return false;
}

void SA(double a){   //模拟退火算法
    double t = T;
    int s = 0;
    while(s < 2){    //当s = 2时结束搜索
        //printf("s%d\n",s);
        int bchange = 0;
        int l = L;
        while(l--){  //对每个不同的t值遍历l次
            Path newpath = GenerateNewpath();
            if(newpath.len < path.len){
                path = newpath;
                bchange++;
            }
            else if(Accept(path.len-newpath.len,t)){
                path = newpath;
                bchange++;
            }
        }
        t = a*t;
        if(!bchange){
            s++;
        }
        else{
            s = 0;
        }
    }
    return ;

}

void Calculatedis(){ //利用欧式距离计算两个城市之间距离
    for(int i = 0;i < n;i++){
        for(int j = 0;j < n;j++){
     d[i][j] = sqrt((city[i].x-city[j].x)*(city[i].x-city[j].x)+(city[i].y-city[j].y)*(city[i].y-city[j].y));
        }
    }
    return ;
}

5.代码运行及测试
我在做的时候是不固定起点且默认所有城市都有一条路，若没有路就将路径长度置为无穷大。这里我的L设置为20000，t初试设置为100。数据第一行是表示城市个数，第二行表示衰减率。在网上找了一些测试样例，结果如下：
测试样例：
样例1：
5
0.9
100000 3 1 5 8
3 100000 6 7 9
1 6 100000 4 2
5 7 4 100000 3
8 9 2 3 100000
结果：

上面两行表示初试化的路径与长度，找到的最短路径是4 -> 2 -> 0 -> 1 -> 3 -> 4，路径长度为16。
当将衰减率为0.8时结果还是一样，当衰减率为0.5时结果还是一样，应该是数据比较简单。

样例2：
4
0.9
0 3 6 7
5 0 2 3
6 4 0 2
3 7 5 0
结果：

样例2：
这个样例是我在网上找到的，里面提供了各个城市的坐标，所以我用欧式距离算了一下城市距离矩阵。
27
0.9
41 94
37 84
53 67
25 62
7 64
2 99
68 58
71 44
54 62
83 69
64 60
18 54
22 60
83 46
91 38
25 38
24 42
58 69
71 71
74 78
87 76
18 40
13 40
82 7
62 32
58 35
45 21

结果：

这个跑出来最短距离是408.368。

样例3：
52
0.9
565
575
25
185
345
750
945
685
845
655
880
660
25
230
525
1000
580
1175
650
1130
1605
620
1220
580
1465
200
1530
5
845
680
725
370
145
665
415
635
510
875
560
365
300
465
520
585
480
415
835
625
975
580
1215
245
1320
315
1250
400
660
180
410
250
420
555
575
665
1150
1160
700
580
685
595
685
610
770
610
795
645
720
635
760
650
475
960
95
260
875
920
700
500
555
815
830
485
1170
65
830
610
605
625
595
360
1340
725
1740
245
结果：

一共52个城市，衰减率为0.9，最终结果为7753.772，运行时间0.881s。
同一个数据，当衰减率为0.7时，结果如下：

可以发现结果变差了，我们可以从中推出当衰减率为0.7时，算法陷入一个局部最优出不来了，而且此时花费的时间也比之前衰减率为0.9时快了将近一倍。

样例4：
34
0.9
103.73 36.03
101.74 36.56
104.06 30.67
114.48 38.03
102.73 25.04
106.71 26.57
114.31 30.52
113.65 34.76
117 36.65
118.78 32.04
117.27 31.86
120.19 30.26
115.89 28.68
119.3 26.08
113.23 23.16
113 28.21
110.35 20.02
123.38 41.8
125.35 43.88
126.63 45.75
112.53 37.87
108.95 34.27
121.3 25.03
116.46 39.92
121.48 31.22
106.54 29.59
117.2 39.13
111.65 40.82
108.33 22.84
91.11 29.97
106.27 38.47
87.68 43.77
114.17 22.28
113.54 22.19
这是中国34个城市的坐标，跑完结果如下:

最短的距离为157.64，运行时间0.988秒，算法效率还是挺高的。