【建模算法】Python调用scikit-opt工具箱中的模拟退火算法求解TSP问题

TSP (traveling salesman problem，旅行商问题)是典型的NP完全问题，即其最坏情况下的时间复杂度随着问题规模的增大按指数方式增长，到目前为止还未找到一个多项式时间的有效算法。本文探讨了Python调用scikit-opt工具箱中的模拟退火算法求解TSP问题。

​ 本案例以31个城市为例，假定31个城市的位置坐标如表1所列。寻找出一条最短的遍历31个城市的路径。

scikit-opt工具箱是一个封装了7种启发式算法的 Python 代码库，包含差分进化算法、遗传算法、粒子群算法、模拟退火算法、蚁群算法、鱼群算法、免疫优化算法。官网链接：https://scikit-opt.github.io/scikit-opt/#/zh/

pip install scikit-opt

主要是将自己问题的数据输入到封装的代码中，操作简单，修改码量很少
其中只有两行代码需要修改

num_points = 31                #num_points指需要遍历的节点个数
points_coordinate = np.random.rand(num_points, 2)    #points_coordinate是自定义问题中的节点坐标信息

其中的best_points表示最短距离的num_points点的序列（即points_coordinate矩阵中第几个点，但没有回到起点）

from sko.SA import SA_TSP

sa_tsp = SA_TSP(func=cal_total_distance, x0=range(num_points), T_max=100, T_min=1, L=10 * num_points)

best_points, best_distance = sa_tsp.run()
print(best_points, best_distance, cal_total_distance(best_points))

其中的best_points_ 表示从起点出发并回到起点的路径。
（其实与best_points差别就是在best_points后加入起点（best_points的第一个点））

from matplotlib.ticker import FormatStrFormatter

fig, ax = plt.subplots(1, 2)

best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(sa_tsp.best_y_history)
ax[0].set_xlabel("Iteration")
ax[0].set_ylabel("Distance")
ax[1].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1],
           marker='o', markerfacecolor='b', color='c', linestyle='-')
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].set_xlabel("Longitude")
ax[1].set_ylabel("Latitude")
plt.show()

城市序列与距离：

TSP图和Loss图：

动画演示：

python源码：

#模拟退火算法求解31座城市的TSP问题完整代码：
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
import sys
from sko.SA import SA_TSP
from matplotlib.ticker import FormatStrFormatter
from time import perf_counter

def cal_total_distance(routine):
    '''The objective function. input routine, return total distance.
    cal_total_distance(np.arange(num_points))
    '''    
    num_points, = routine.shape
    return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])

def print_route(best_points):
    result_cur_best=[]
    for i in best_points:
        result_cur_best+=[i]
    for i in range(len(result_cur_best)):
        result_cur_best[i] += 1
    result_path = result_cur_best
    result_path.append(result_path[0])
    return result_path

if __name__=="__main__":
    file_name = 'data.csv'
    points_coordinate = np.loadtxt(file_name, delimiter=',')
    num_points = points_coordinate.shape[0]
    distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
    distance_matrix = distance_matrix * 111000  # 1 degree of lat/lon ~ = 111000m，原模型将坐标数据视为经纬度数据，要进行数据转换后使用
    
    start=perf_counter()       #计时开始
    # 执行模拟退火(SA)算法
    sa_tsp = SA_TSP(func=cal_total_distance, x0=range(num_points), T_max=100, T_min=1, L=10 * num_points)  #调用工具箱
    
    # 结果输出
    best_points, best_distance = sa_tsp.run()
    print("运行时间是: {:.5f}s".format(perf_counter()-start))   #计时结束
    print("最优路线：", print_route(best_points))
    print("最优值：", cal_total_distance(best_points)/111000)   #数据还原

    best_points_ = np.concatenate([best_points, [best_points[0]]])
    best_points_coordinate = points_coordinate[best_points_, :]
    sa_tsp.best_y_history=np.array(sa_tsp.best_y_history)/111000   #数据还原

    #绘图
    #fig, ax = plt.subplots(1, 2)
    plt.rcParams['font.sans-serif'] = 'SimHei'  # 设置中文显示
    plt.rcParams['axes.unicode_minus'] = False
    fig1, ax1 = plt.subplots(1, 1)
    ax1.set_title('优化过程', loc='center')
    ax1.plot(sa_tsp.best_y_history)
    ax1.set_xlabel("迭代次数")
    ax1.set_ylabel("最优值")

    fig2, ax2 = plt.subplots(1, 1)
    ax2.set_title('轨迹图', loc='center')
    line=ax2.plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], marker='>', mec='r', mfc='w',label=u'路线')
    ax2.legend()  # 让图例生效
    for i in range(num_points):
        plt.text(best_points_coordinate[:, 0][i] + 0.05, best_points_coordinate[:, 1][i] + 0.05, str(best_points[i]+1), color='red')
    ax2.xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
    ax2.yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
    ax2.set_xlabel("横坐标")
    ax2.set_ylabel("纵坐标")
    plt.show()