scipy中的旅行商问题

Question

我该如何在Python中解决旅行商问题？我没有找到任何库，应该有办法使用scipy的优化函数或其他库来解决。

我这个懒得要命的暴力破解方案是：

tsp_solution = min( (sum( Dist[i] for i in izip(per, per[1:])), n, per) for n, per in enumerate(i for i in permutations(xrange(Dist.shape[0]), Dist.shape[0])) )[2]

这里的Dist（numpy.array）是距离矩阵。如果Dist太大，这个方法会花费很长时间。

Answer 1

在Python中有一个叫做“graph”的库，专门用来解决网络问题（比如理论计算机科学和运筹学中的顶点覆盖优化问题），这个库叫做NetworkX。在它的文档中（参考/算法/近似和启发式方法），你可以找到一些例子（对于这类问题，它可能比scipy更方便，包括最短路径问题、最小成本流问题，等等）。还有一个叫做旅行商问题：

NP难的图形旅行商问题（GTSP）是要找到一个总权重最小的闭合路径，这条路径需要访问无向边加权且不一定是完全图中的每个顶点。

"""
==========================
Traveling Salesman Problem
==========================

This is an example of a drawing solution of the traveling salesman problem

The function used to produce the solution is `christofides <networkx.algorithms.approximation.traveling_salesman.christofides>`,
where given a set of nodes, it calculates the route of the nodes
that the traveler has to follow in order to minimize the total cost.
"""

import matplotlib.pyplot as plt
import networkx as nx
import networkx.algorithms.approximation as nx_app
import math

G = nx.random_geometric_graph(20, radius=0.4, seed=3)
pos = nx.get_node_attributes(G, "pos")

# Depot should be at (0,0)
pos[0] = (0.5, 0.5)

H = G.copy()


# Calculating the distances between the nodes as edge's weight.
for i in range(len(pos)):
    for j in range(i + 1, len(pos)):
        dist = math.hypot(pos[i][0] - pos[j][0], pos[i][1] - pos[j][1])
        dist = dist
        G.add_edge(i, j, weight=dist)

cycle = nx_app.christofides(G, weight="weight")
edge_list = list(nx.utils.pairwise(cycle))

# Draw closest edges on each node only
nx.draw_networkx_edges(H, pos, edge_color="blue", width=0.5)

# Draw the route
nx.draw_networkx(
    G,
    pos,
    with_labels=True,
    edgelist=edge_list,
    node_color = 'r',
    edge_color="red",
    node_size=200,
    width=3,
)

print("The route of the traveller is:", cycle)
plt.show()

注意：在networkX的文档中，有不同的算法可以用来解决你的问题。

公式化：覆盖旅行商问题的分支与切割算法

Answer 2

scipy.optimize 里的函数并不是特别适合直接用来解决旅行商问题（TSP）。如果你想要一个简单的解决方案，我推荐使用 2-opt 算法，这是一种被广泛接受的解决 TSP 的算法，而且实现起来相对简单。下面是我实现这个算法的代码：

import numpy as np

# Calculate the euclidian distance in n-space of the route r traversing cities c, ending at the path start.
path_distance = lambda r,c: np.sum([np.linalg.norm(c[r[p]]-c[r[p-1]]) for p in range(len(r))])
# Reverse the order of all elements from element i to element k in array r.
two_opt_swap = lambda r,i,k: np.concatenate((r[0:i],r[k:-len(r)+i-1:-1],r[k+1:len(r)]))

def two_opt(cities,improvement_threshold): # 2-opt Algorithm adapted from https://en.wikipedia.org/wiki/2-opt
    route = np.arange(cities.shape[0]) # Make an array of row numbers corresponding to cities.
    improvement_factor = 1 # Initialize the improvement factor.
    best_distance = path_distance(route,cities) # Calculate the distance of the initial path.
    while improvement_factor > improvement_threshold: # If the route is still improving, keep going!
        distance_to_beat = best_distance # Record the distance at the beginning of the loop.
        for swap_first in range(1,len(route)-2): # From each city except the first and last,
            for swap_last in range(swap_first+1,len(route)): # to each of the cities following,
                new_route = two_opt_swap(route,swap_first,swap_last) # try reversing the order of these cities
                new_distance = path_distance(new_route,cities) # and check the total distance with this modification.
                if new_distance < best_distance: # If the path distance is an improvement,
                    route = new_route # make this the accepted best route
                    best_distance = new_distance # and update the distance corresponding to this route.
        improvement_factor = 1 - best_distance/distance_to_beat # Calculate how much the route has improved.
    return route # When the route is no longer improving substantially, stop searching and return the route.

这里有一个使用这个函数的例子：

# Create a matrix of cities, with each row being a location in 2-space (function works in n-dimensions).
cities = np.random.RandomState(42).rand(70,2)
# Find a good route with 2-opt ("route" gives the order in which to travel to each city by row number.)
route = two_opt(cities,0.001)

还有这是在图上显示的近似解决路径：

import matplotlib.pyplot as plt
# Reorder the cities matrix by route order in a new matrix for plotting.
new_cities_order = np.concatenate((np.array([cities[route[i]] for i in range(len(route))]),np.array([cities[0]])))
# Plot the cities.
plt.scatter(cities[:,0],cities[:,1])
# Plot the path.
plt.plot(new_cities_order[:,0],new_cities_order[:,1])
plt.show()
# Print the route as row numbers and the total distance travelled by the path.
print("Route: " + str(route) + "\n\nDistance: " + str(path_distance(route,cities)))

如果算法的速度对你很重要，我建议你提前计算好城市之间的距离，并把这些距离存储在一个矩阵里。这样可以大大减少计算时间。

编辑：自定义起点和终点

如果你想要一条非循环的路径（也就是说，路径的结束地点和起始地点不同），你需要修改路径距离的公式为：

path_distance = lambda r,c: np.sum([np.linalg.norm(c[r[p+1]]-c[r[p]]) for p in range(len(r)-1)])

然后使用下面的代码重新排列城市，以便绘图：

new_cities_order = np.array([cities[route[i]] for i in range(len(route))])

按照现在的代码，起始城市是固定为 cities 中的第一个城市，而结束城市是可以变化的。

如果你想让结束城市是 cities 中的最后一个城市，可以通过修改 two_opt() 中的 swap_first 和 swap_last 的范围来限制可交换城市的范围，代码如下：

for swap_first in range(1,len(route)-3):
    for swap_last in range(swap_first+1,len(route)-1):

如果你想让起始城市和结束城市都可以变化，可以通过扩展 swap_first 和 swap_last 的范围来实现，代码如下：

for swap_first in range(0,len(route)-2):
    for swap_last in range(swap_first+1,len(route)):