针对中国31个省会城市的TSP问题:
通过模拟退火算法解决,坐标输入:
[[1304, 2312], [3639, 1315], [4177, 2244], [3712, 1399], [3488, 1535], [3326, 1556], [3238, 1229], [4196, 1004], [4312, 790], [4386, 570], [3007, 1970], [2562, 1756], [2788, 1491], [2381, 1676], [1332, 695], [3715, 1678], [3918, 2179], [4061, 2370], [3780, 2212], [3676, 2578], [4029, 2838], [4263, 2931], [3429, 1908], [3507, 2367], [3394, 2643], [3439, 3201], [2935, 3240], [3140, 3550], [2545, 2357], [2778, 2826], [2370, 2975]]
这里有3种扰动操作:1.两元素互换 2.任意长度片段进行逆序操作 3.选取任意元素进行插入操作
具体实现:
import sys
import math
import os
import random
from matplotlib import pyplot as plt
END_T = 1 # 结束温度
T = 100 # 温度
L = 100 * 31 # 马尔科夫链长
K = 0.9 # 降温系数
dist = [] # 邻接矩阵
city = [[1304, 2312], [3639, 1315], [4177, 2244], [3712, 1399], [3488, 1535], [3326, 1556],
[3238, 1229], [4196, 1004], [4312, 790], [4386, 570], [3007, 1970], [2562, 1756],
[2788, 1491], [2381, 1676], [1332, 695], [3715, 1678], [3918, 2179], [4061, 2370],
[3780, 2212], [3676, 2578], [4029, 2838], [4263, 2931], [3429, 1908], [3507, 2367],
[3394, 2643], [3439, 3201], [2935, 3240], [3140, 3550], [2545, 2357], [2778, 2826],
[2370, 2975]] # 各城市坐标
initial_path = list(range(0, 31))
losses = []
def __initial__():
for i in city:
dist.append([math.sqrt(math.pow(i[0] - j[0], 2) + math.pow(i[1] - j[1], 2)) for j in city])
def E(x: list):
length = dist[x[0]][x[30]]
i = 1
while i < 31:
length += dist[x[i]][x[i - 1]]
i += 1
return length
def accept_probability(xj, xi):
return math.exp((E(xi) - E(xj)) / T)
def random_operation(path):
""" 互换操作 """
a, b = random.randint(1, 30), random.randint(1, 30)
while a == b: # 不能自己和自己交换,且从点0出发
b = random.randint(1, 30)
path[a], path[b] = path[b], path[a]
return path
def random_operation2(path):
""" 逆序操作 """
a, b = random.randint(1, 30), random.randint(1, 30)
if a > b:
a, b = b, a
while a < b:
path[a], path[b] = path[b], path[a]
a += 1
b -= 1
return path
def random_operation3(path):
""" 插入操作 """
a, b = random.randint(1, 30), random.randint(1, 30)
if a == b: # 没有操作
return path
if a < b: # 后插
i, temp = a, path[a]
while i < b:
path[i] = path[i + 1]
i += 1
path[b] = temp
else: # 前插
i, temp = a, path[a]
while i > b:
path[i] = path[i - 1]
i -= 1
path[b] = temp
return path
def 模拟退火(operator=None):
if operator is None:
operator = random_operation
losses = []
xi, xj = initial_path.copy(), None
x_best = xi
temperature = T
while temperature > END_T:
temperature *= K
loop = L
while loop > 0:
loop -= 1
xj = xi.copy()
losses.append(E(x_best))
operator(xj) # 扰动操作
if E(xj) < E(xi):
xi = xj.copy() # 距离更短则无条件接受
if E(xj) < E(x_best):
x_best = xj.copy() # 全局最小
continue
probablity = accept_probability(xj, xi)
# 距离更长,按概率接受
if random.random() < probablity: # 小于概率,接受
xi = xj.copy()
# print('Dist:', E(xi), '\tTemperature=', temperature)
return [x_best, losses]
def get_union_path_pic(ele1, ele2, ele3):
eles = [ele1, ele2, ele3]
titles = ['互换操作', '逆序操作', '插入操作']
counter = 0
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
while counter < 3:
plt.subplot(3, 2, 2 * counter + 1)
plt.title(titles[counter])
i = 1
while i < 31:
plt.plot([city[eles[counter][0][i]][0], city[eles[counter][0][i - 1]][0]],
[city[eles[counter][0][i]][1], city[eles[counter][0][i - 1]][1]], color='blue')
i += 1
plt.plot([city[eles[counter][0][0]][0], city[eles[counter][0][30]][0]],
[city[eles[counter][0][0]][1], city[eles[counter][0][30]][1]], color='red')
# line
plt.subplot(3, 2, 2 * counter + 2)
plt.plot(eles[counter][1])
plt.title('Distance=' + str(E(eles[counter][0])))
counter += 1
plt.show()
def get_path_pic(path, losses):
plt.subplot(1, 2, 1)
i = 1
while i < 31:
plt.plot([city[path[i]][0], city[path[i - 1]][0]], [city[path[i]][1], city[path[i - 1]][1]], color='blue')
i += 1
plt.plot([city[path[0]][0], city[path[30]][0]], [city[path[0]][1], city[path[30]][1]], color='red')
# line
plt.subplot(1, 2, 2)
plt.plot(losses)
plt.title('Distance=' + str(E(path)))
plt.show()
if __name__ == '__main__':
__initial__()
ele1 = 模拟退火()
print('Dist:', E(ele1[0]), 'Path', ele1[0])
ele2 = 模拟退火(operator=random_operation2)
print('Dist:', E(ele2[0]), 'Path', ele2[0])
ele3 = 模拟退火(operator=random_operation3)
print('Dist:', E(ele3[0]), 'Path', ele3[0])
get_union_path_pic(ele1, ele2, ele3)
运行效果: