Python 中 TSP 的可扩展实现

Question

我有（比方说）100 分。我想在它们之间生成一条路径，越短越好。这是Travelling salesperson problem。自从TSP is NP-hard 之后，我对没有找到全局解决方案感到满意。 I 方法可以快速提供解决方案并很好地扩展。

生成示例点：

import numpy as np
points = np.random.RandomState(42).rand(100,2)

生成距离矩阵，其中i,j 条目包含point[i] 和point[j] 之间的距离。使用this 回答：

z = np.array([[complex(*c) for c in points]]) # notice the [[ ... ]]
distance_matrix = abs(z.T-z)

我该如何继续，找到一条至少具有局部最小路径长度的最短路径？

---

一些相关问题：

此线程：Travelling Salesman in scipy 提供了用于查找 TSP 解决方案的代码。但它以迭代方式工作，因此如果要访问的点数很大，脚本不会在合理的时间内产生解决方案。

此线程：How to solve the Cumulative Traveling Salesman Problem using or-tools in python? 没有代码答案，也不专注于经典 TSP。

这个帖子：Optimizing a Traveling Salesman Algorithm (Time Traveler Algorithm) 提供了问题的迭代解决方案（这意味着不好的缩放）

Answer 1

这是一个基于simulated annealing (pip install frigidum) 的示例。它很可能比其他解决方案慢。我希望发布的 Google OR 会更好，但由于它是一种不同的方法可能会引起人们的兴趣（尽管用于学习/教育）。

import numpy as np
import frigidum

from frigidum.examples import tsp

points = np.random.RandomState(42).rand(100,2)

tsp.nodes = points
tsp.nodes_count = points.shape[0]

z = np.array([[complex(*c) for c in points]]) # notice the [[ ... ]]
distance_matrix = abs(z.T-z)

tsp.dist_eu = distance_matrix

我们可以使用以下模拟退火方案（alpha 更接近 0.999，和/或以更长的计算成本创建 repeats）

local_opt = frigidum.sa(random_start=tsp.random_start,
           objective_function=tsp.objective_function,
           neighbours=[tsp.euclidian_bomb_and_fix, tsp.euclidian_nuke_and_fix, tsp.route_bomb_and_fix, tsp.route_nuke_and_fix, tsp.random_disconnect_vertices_and_fix],
           copy_state=frigidum.annealing.naked,
           T_start=5,
           alpha=.8,
           T_stop=0.001,
           repeats=10**2,
           post_annealing = tsp.local_search_2opt)

local_opt 是一个元组，它的第一个元素是解决方案，第二个元素是目标值（路线长度）。最后它将输出统计信息和找到的最小值； （这不是绘制的那个）。

(Local) Minimum Objective Value Found: 
   7.46016765

我使用 zabop 绘图代码来绘制解决方案；

import matplotlib.pyplot as plt

plt.scatter(points[:,0],points[:,1])

route = local_opt[0]
for a, b in zip(route[:-1],route[1:]):
    x = points[[a,b]].T[0]
    y = points[[a,b]].T[1]
    plt.plot(x, y,c='r',zorder=-1)
plt.gca().set_aspect('equal')

https://imgur.com/rg9Fgzk

Answer 2

此解决方案使用 Google 的 OR-Tools。 Installortools 通过python -m pip install --upgrade --user ortools。

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

生成点：

import numpy as np
points = np.random.RandomState(42).rand(100,2)

计算有问题的距离矩阵：

z = np.array([[complex(*c) for c in points]]) # notice the [[ ... ]]
distance_matrix_pre = abs(z.T-z)

自从docs说：

注意：由于路由求解器使用整数进行所有计算，因此距离回调必须返回整数距离 任意两个位置。如果data['distance_matrix']的任何条目 不是整数，您需要对矩阵条目或 将回调的值返回为整数。请参阅缩放距离 矩阵示例，显示如何避免由以下原因引起的问题 舍入误差。

我要将距离矩阵中的每个条目乘以 10000，然后向下舍入为 int。如果要使用更多小数位，请乘以更大的数字。

distance_matrix = np.floor(distance_matrix_pre*10000).astype(int).tolist()

继续关注Google's tutorial：

# Create the routing index manager.
manager = pywrapcp.RoutingIndexManager(len(points),1,0)

# Create Routing Model.
routing = pywrapcp.RoutingModel(manager)

def distance_callback(from_index, to_index):
    """Returns the distance between the two nodes."""
    # Convert from routing variable Index to distance matrix NodeIndex.
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    return distance_matrix[from_node][to_node]

transit_callback_index = routing.RegisterTransitCallback(distance_callback)

# Define cost of each arc.
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

输出取决于使用的初始路径选择策略（类似于数值解如何依赖于初始猜测）。这里我们指定初始策略：

# Setting first solution heuristic.
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

用谷歌的话来说：

代码将第一个解决策略设置为PATH_CHEAPEST_ARC， 它通过重复添加为求解器创建初始路径 不会导致先前访问过的权重最小的边 节点（仓库除外）。有关其他选项，请参阅First solution strategy。

使用上面定义的初始策略解决问题（即使用PATH_CHEAPEST_ARC）：

# Solve the problem.
solution = routing.SolveWithParameters(search_parameters)

获取route：

if solution:
    index = routing.Start(0)
    route = [manager.IndexToNode(index)]
    while not routing.IsEnd(index):
        index = solution.Value(routing.NextVar(index))
        route.append(manager.IndexToNode(index))

绘制结果：

plt.scatter(points[:,0],points[:,1])
for a, b in zip(route[:-1],route[1:]):
    x = points[[a,b]].T[0]
    y = points[[a,b]].T[1]
    plt.plot(x, y,c='r',zorder=-1)
plt.gca().set_aspect('equal')

---

其他可以使用的初始策略列表（source）：

initial_strategy_list = \
['AUTOMATIC','PATH_CHEAPEST_ARC','PATH_MOST_CONSTRAINED_ARC','EVALUATOR_STRATEGY','SAVINGS','SWEEP','CHRISTOFIDES','ALL_UNPERFORMED','BEST_INSERTION','PARALLEL_CHEAPEST_INSERTION','LOCAL_CHEAPEST_INSERTION','GLOBAL_CHEAPEST_ARC','LOCAL_CHEAPEST_ARC','FIRST_UNBOUND_MIN_VALUE']

如果想增加对解决方案的信心，可以遍历所有解决方案并找到最佳解决方案。