#用SA算法求解旅行商问题

#问题描述

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

#算法设计

采用模拟退火算法计算旅行商问题, 模拟退火算法(Simulated Annealing Algorithm, 简称SA算法), 源于对固体退火过程的模拟, 采用Metropolis接受准则, 并用一组称为冷却表的参数控制算法进程, 使算法在多项式时间里给出一个近似最优解.

#Metropolis算法

模拟退火算法用Metropolis算法产生组合优化问题解的序列. 并由Metropolis准则对应的转移概率P.

$$ P(i=>j)= \begin{cases} 1,&{当f(j) \leq f(i)} \\ e^{\frac{f(i)−f(j)}{t}},&{否则} \end{cases} \tag{1} $$

#冷却进度表

冷却表(cooling schedule)是一组控制算法进程的参数, 用以逼近模拟退火算法的渐进收敛性态, 使算法在有限时间内执行迭代过程后返回一个近似最优解.
冷却表是影响模拟退火算法实验性能的重要因素, 其合理选取是算法应用的关键.
一个冷却表应当规定下述参数:

  • 控制参数t的初值t0
  • 控制参数t的衰减函数
  • 控制参数t的终值tf(停止准则)
  • Mapkob链长Lk

#初值t0

控制参数初值t0要足够大, 才能使算法进程在合理时间里搜索尽可能大的解空间范围.
本实验选取t0为100.

#衰减函数

在控制参数的衰减函数应使tk的衰减以小为宜.

$$ \alpha(t)=0.9*t \tag{2} $$

本实验选取的衰减函数如公式2所示.

#t的终值tf(停止准则)

控制参数终值tf通常由停止准则决定. 合理的停止准则既要确保算法收敛于某一近似解, 又要使最终解具有一定的质量.
本实验选取的停止准则为连续两个Mapkob链路径无变化则程序终止.

#Mapkob链长Lk

Mapkob链长Lk在控制参数t的衰减函数已选定的前提下, Lk应选得在控制参数的每一取值上都能恢复准平衡.
本实验选取的Mapkob链长为定长20000.

#新解的产生

采用2变换法. 任选序号u和v(u < v), 将 u和v及其之间的顺序逆转.

$$ \pi_{1}...\pi_{u-1}\color{red}{\pi_{u}\pi_{u+1}...\pi_{v-1}\pi_{v}}\pi_{v+1}...\pi_{n} \\ \downarrow{} \\ \pi_{1}...\pi_{u-1}\color{red}{\pi_{v}\pi_{v-1}...\pi_{u-1}\pi_{u}}\pi_{v+1}...\pi_{n} \tag{3} $$

#程序流程

程序流程如图1所示

图1 程序流程图

#核心伪代码

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
# 初始化路径
x1 = random.sample(range(citynum), citynum)

# 计算按此路径运动所花费的代价
f1 = calculate(x1)

# 初始化起始温度
T = T_max

flag = 0

# 当起始温度低于预设的最低温度时停止循环
while True:
    flag = flag + 1

    # Mapkob链长, 即迭代次数
    for i in range(iter):

        # 通过2变换法生成新路径
        x2 = generate(x1)

        # 计算新解所花费的代价
        f2 = calculate(x2)
        
        # 计算两个解的代价差
        df = f2 - f1

        # 用Metropolis算法判断是否接受新解
        if accept(df, T):
            x1 = x2
            f1 = calculate(x1)
            flag = 0
    
    # 模拟温度衰减
    T *= rate

    # 连续两次Mapkob链路径无变化则停止循环
    if flag >= 2:
        break

# 得到最终路径结果
x = x1

#代码运行及测试

代码的运行结果如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
time: 14.867s
result: 0 -> 3 -> 2 -> 1
cost: 108

time: 16.764s
result: 5 -> 6 -> 7 -> 1 -> 2 -> 9 -> 0 -> 8 -> 4 -> 3
cost: 229

time: 24.910s
result: 12 -> 15 -> 8 -> 0 -> 18 -> 3 -> 21 -> 16 -> 27 -> 20
 -> 29 -> 14 -> 10 -> 6 -> 28 -> 25 -> 17 -> 5 -> 4 -> 11
 -> 24 -> 1 -> 26 -> 9 -> 23 -> 7 -> 19 -> 13 -> 2 -> 22
cost: 576

由于数据计算结果过慢, 在代码中额外增加了终止条件: 当温度低于0.01时程序停止.
由于前两个测试数据城市数量较少, 已经得到最优解, 后续改进实验仅以30个城市的结果作为对比.
考虑到初始值为随机选取的, 且最终的路径为一个头尾相接的"环", 修改新解的产生方法

$$ \color{red}{\pi_{1}...\pi_{u-1}}\pi_{u}\pi_{u+1}...\pi_{v-1}\pi_{v}\color{red}{\pi_{v+1}...\pi_{n}} \\ \downarrow{} \\ \color{red}{\pi_{1}...\pi_{v+1}}\pi_{u}\pi_{u+1}...\pi_{v-1}\pi_{v}\color{red}{\pi_{u-1}...\pi_{n}} \tag{4} $$

代码的运行结果如下:

1
2
3
4
5
time: 22.672s
result: 9 -> 17 -> 5 -> 18 -> 16 -> 27 -> 20 -> 21 -> 4 -> 1
 -> 12 -> 15 -> 8 -> 3 -> 13 -> 25 -> 2 -> 29 -> 14 -> 10
 -> 6 -> 28 -> 22 -> 11 -> 24 -> 7 -> 19 -> 23 -> 26 -> 0
cost: 565

增大变换的自由度, 只交换u和v的两个序号

$$ \pi_{1}...\pi_{u-1}\color{red}{\pi_{u}}\pi_{u+1}...\pi_{v-1}\color{red}{\pi_{v}}\pi_{v+1}...\pi_{n} \\ \downarrow{} \\ \pi_{1}...\pi_{u-1}\color{red}{\pi_{v}}\pi_{u+1}...\pi_{v-1}\color{red}{\pi_{u}}\pi_{v+1}...\pi_{n} \tag{5} $$

代码的运行结果如下:

1
2
3
4
5
time: 11.749s
result: 0 -> 2 -> 22 -> 11 -> 7 -> 28 -> 6 -> 19 -> 23 -> 26
 -> 9 -> 21 -> 12 -> 13 -> 16 -> 27 -> 20 -> 15 -> 8 -> 24
 -> 4 -> 1 -> 5 -> 25 -> 17 -> 18 -> 3 -> 29 -> 14 -> 10
cost: 475

观察代码运行状态, 都是温度达到最小值以下后退出的程序, 可能未达到当前参数下的最佳结果, 故将循环终止条件改为只有连续两个Mapkob链路径无变化.
代码的运行结果如下:

1
2
3
4
5
time: 7.6435s
result: 8 -> 2 -> 21 -> 6 -> 28 -> 29 -> 14 -> 10 -> 17 -> 23
 -> 26 -> 9 -> 12 -> 15 -> 19 -> 25 -> 1 -> 5 -> 13 -> 16
 -> 27 -> 20 -> 18 -> 22 -> 11 -> 0 -> 24 -> 4 -> 7 -> 3
cost: 467

所以结果如表1所示:

表1 代码运行结果
代码版本 时间 花费
1 24.910s 576
2 22.672s 565
3 11.749s 475
4 7.6435s 467

#结论

算法在多项式时间里给出一个近似最优解.

#附录

#城市矩阵生成代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import random

def generate_city(num, min=10, max=100):
	data = []
	for i in range(num):
		temp_data = []
		for j in range(num):
			if(i == j):
				temp_data.append(0)
			else:
				temp_data.append(random.randint(min, max))
		data.append(temp_data)
	return data

print(generate_city(4))
print(generate_city(10))
print(generate_city(30))

#模拟退火算法代码

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
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import time
import random

class SA(object):
    def __init__(self, citymap, T_max=100, T_min=0.01, iter=20000, rate=0.90):
        self.citymap = citymap
        self.citynum = len(citymap)
        self.T_max = T_max
        self.T_min = T_min
        self.iter = iter
        self.rate = rate

        self.i = random.sample(range(self.citynum), self.citynum)

        start_time = time.time()
        self.solve()
        end_time = time.time()

        print('time: ' + str(end_time - start_time)[:6] + 's')
        print('result: ' + ' -> '.join(str(i) for i in self.x))
        print('cost: ' + str(self.calculate(self.x)))

    def calculate(self, city_list):
        cost = 0
        for i in range(self.citynum - 1):
            cost = cost + self.citymap[city_list[i]][city_list[i + 1]]
        cost = cost + self.citymap[city_list[i + 1]][city_list[0]]
        return cost

    # def generate(self, city_list):
    #     u, v = random.sample(range(self.citynum + 1), 2)
    #     if u > v:
    #         temp_list = city_list[u:] + city_list[:v]
    #         temp_list.reverse()
    #         return temp_list[self.citynum-u:] + city_list[v:u] + temp_list[:self.citynum-u]
    #     else:
    #         temp_list = city_list[u:v]
    #         temp_list.reverse()
    #         return city_list[:u] + temp_list + city_list[v:]

    # def generate(self, city_list):
    #     u, v = random.sample(range(self.citynum + 1), 2)
    #     if u > v:
    #         u, v = v, u
    #     temp_list = city_list[u:v]
    #     temp_list.reverse()
    #     return city_list[:u] + temp_list + city_list[v:]

    def generate(self, city_list):
        u, v = random.sample(range(self.citynum), 2)
        if u > v:
            u, v = v, u
        return city_list[:u] + city_list[v:v+1] + city_list[u+1:v] + city_list[u:u+1] + city_list[v+1:]

    def accept(self, df, T):
        if df < 0:
            return True
        else:
            if np.exp(-df / T) >= random.random():
                return True
            return False

    def solve(self):
        x1 = self.i
        f1 = self.calculate(x1)
        T = self.T_max
        flag = 0
        while True:
            flag = flag + 1
            for i in range(self.iter):
                x2 = self.generate(x1)
                f2 = self.calculate(x2)
                df = f2 - f1
                if self.accept(df, T):
                    x1 = x2
                    f1 = self.calculate(x1)
                    flag = 0
            T *= self.rate
            if flag >= 2:
                break
        self.x = x1

if __name__ == '__main__':
    SA([[0, 48, 46, 48], [1, 0, 63, 4], [39, 3, 0, 59], [67, 28, 56, 0]])
    SA([[0, 76, 62, 85, 81, 99, 78, 99, 37, 87],
        [31, 0, 24, 61, 56, 89, 84, 99, 48, 44],
        [12, 54, 0, 43, 53, 86, 73, 90, 37, 22],
        [77, 12, 23, 0, 64, 13, 61, 29, 94, 54],
        [17, 90, 81, 27, 0, 61, 94, 36, 79, 64],
        [47, 27, 66, 82, 87, 0, 17, 57, 91, 55],
        [98, 18, 93, 72, 89, 90, 0, 48, 47, 94],
        [88, 13, 97, 78, 80, 54, 21, 0, 10, 96],
        [32, 87, 75, 27, 10, 94, 36, 21, 0, 14],
        [18, 97, 75, 91, 57, 78, 76, 70, 13, 0]])
    SA([[0, 84, 19, 73, 83, 56, 55, 78, 68, 10, 33, 61, 23, 41, 17,
         81, 12, 22, 24, 60, 85, 53, 82, 97, 13, 44, 54, 96, 54, 18],
        [56, 0, 25, 76, 60, 11, 38, 35, 98, 76, 25, 67, 29, 51, 53,
         29, 94, 70, 76, 74, 83, 29, 54, 64, 75, 33, 25, 53, 23, 62],
        [43, 89, 0, 95, 100, 93, 38, 95, 77, 34, 33, 70, 90, 75, 66,
         75, 26, 81, 90, 19, 91, 27, 10, 25, 17, 37, 18, 95, 13, 27],
        [31, 26, 63, 0, 80, 71, 64, 70, 10, 54, 70, 47, 78, 12, 18,
         89, 69, 90, 59, 37, 24, 11, 90, 87, 71, 48, 38, 76, 84, 13],
        [24, 16, 72, 64, 0, 14, 64, 11, 82, 52, 87, 14, 43, 75, 32,
         79, 66, 54, 66, 50, 55, 13, 14, 85, 31, 58, 37, 60, 33, 83],
        [11, 15, 56, 46, 37, 0, 50, 86, 55, 76, 17, 59, 58, 15, 37,
         12, 22, 39, 33, 97, 65, 88, 57, 56, 64, 29, 69, 51, 66, 94],
        [63, 71, 50, 34, 97, 87, 0, 37, 75, 58, 68, 21, 57, 74, 75,
         84, 20, 90, 77, 18, 99, 59, 19, 22, 40, 33, 44, 81, 11, 34],
        [89, 79, 22, 15, 78, 55, 83, 0, 21, 74, 22, 69, 53, 25, 57,
         100, 88, 53, 68, 18, 99, 71, 81, 56, 72, 82, 100, 21, 14, 55],
        [19, 88, 13, 28, 61, 79, 73, 31, 0, 14, 85, 61, 36, 39, 62,
         44, 69, 63, 25, 55, 51, 100, 77, 36, 11, 78, 37, 76, 88, 66],
        [15, 99, 36, 61, 79, 85, 57, 26, 39, 0, 22, 54, 13, 78, 95,
         99, 94, 11, 60, 12, 43, 20, 60, 28, 24, 95, 94, 50, 51, 77],
        [34, 73, 96, 100, 81, 93, 21, 63, 34, 19, 0, 63, 47, 95, 36,
         20, 18, 13, 99, 58, 30, 98, 29, 75, 85, 58, 99, 39, 25, 79],
        [16, 90, 15, 85, 79, 53, 44, 13, 83, 86, 27, 0, 37, 41, 42,
         13, 50, 45, 97, 36, 96, 87, 59, 57, 42, 36, 85, 17, 55, 98],
        [99, 85, 69, 70, 44, 34, 39, 88, 35, 13, 26, 40, 0, 11, 19,
         10, 22, 28, 21, 52, 54, 91, 76, 78, 53, 66, 87, 17, 54, 87],
        [80, 66, 29, 77, 41, 92, 61, 25, 10, 100, 47, 39, 92, 0, 30,
         37, 10, 81, 59, 24, 75, 74, 34, 65, 82, 53, 100, 82, 69, 92],
        [87, 73, 97, 59, 72, 77, 34, 86, 57, 90, 10, 74, 38, 18, 0,
         90, 19, 55, 17, 71, 63, 95, 83, 88, 45, 32, 40, 51, 70, 71],
        [78, 51, 43, 92, 50, 81, 91, 59, 11, 48, 52, 58, 30, 79, 48,
         0, 26, 28, 52, 17, 24, 11, 63, 55, 46, 76, 59, 32, 43, 96],
        [77, 64, 11, 77, 89, 64, 23, 31, 81, 43, 55, 30, 38, 92, 63,
         67, 0, 100, 31, 40, 59, 39, 41, 35, 72, 19, 37, 12, 76, 65],
        [33, 80, 18, 86, 35, 18, 81, 22, 32, 47, 72, 43, 60, 73, 88,
         49, 30, 0, 13, 46, 84, 72, 25, 15, 88, 22, 69, 60, 97, 97],
        [57, 68, 69, 33, 74, 36, 24, 49, 32, 11, 48, 38, 56, 20, 60,
         24, 16, 69, 0, 60, 27, 98, 33, 26, 73, 37, 40, 83, 29, 89],
        [83, 79, 23, 58, 33, 65, 76, 97, 65, 17, 93, 87, 66, 25, 95,
         28, 90, 78, 33, 0, 64, 90, 41, 12, 74, 30, 100, 86, 73, 49],
        [57, 66, 50, 82, 42, 56, 54, 61, 55, 65, 23, 98, 78, 54, 44,
         17, 60, 89, 17, 44, 0, 14, 48, 40, 100, 44, 58, 11, 55, 25],
        [88, 18, 94, 98, 32, 55, 20, 40, 90, 52, 80, 30, 15, 89, 60,
         72, 20, 59, 77, 31, 83, 0, 83, 51, 90, 48, 63, 48, 42, 36],
        [98, 51, 76, 52, 42, 80, 57, 36, 98, 64, 28, 25, 18, 62, 22,
         59, 13, 95, 65, 69, 91, 28, 0, 88, 64, 38, 97, 66, 14, 96],
        [32, 20, 71, 94, 45, 93, 98, 10, 99, 78, 31, 35, 56, 18, 91,
         79, 23, 16, 55, 71, 53, 59, 75, 0, 58, 66, 11, 100, 81, 53],
        [68, 13, 14, 62, 23, 74, 55, 11, 83, 64, 73, 57, 73, 45, 28,
         81, 40, 93, 78, 29, 97, 97, 36, 17, 0, 70, 77, 65, 75, 31],
        [45, 12, 16, 84, 72, 22, 47, 93, 46, 37, 85, 89, 32, 92, 59,
         93, 62, 12, 80, 21, 91, 80, 49, 35, 90, 0, 68, 54, 17, 45],
        [21, 49, 13, 58, 100, 98, 56, 12, 74, 13, 87, 36, 48, 59, 55,
         71, 21, 70, 45, 38, 63, 91, 99, 100, 41, 38, 0, 91, 13, 48],
        [69, 59, 58, 28, 25, 98, 54, 39, 75, 33, 60, 52, 35, 43, 15,
         34, 90, 58, 62, 34, 11, 23, 41, 50, 81, 75, 83, 0, 60, 54],
        [36, 44, 92, 14, 41, 100, 17, 19, 89, 43, 30, 27, 65, 80, 65,
         81, 35, 90, 17, 75, 26, 89, 13, 63, 16, 35, 22, 86, 0, 19],
        [85, 23, 52, 95, 40, 77, 71, 81, 45, 66, 47, 52, 53, 69, 11,
         22, 22, 86, 90, 35, 33, 12, 36, 51, 86, 43, 59, 65, 64, 0]])