使用Python解决地理旅行推销员问题

使用 pyconcorde 在实际路线问题中找到最优解决方案

著名的旅行推销员问题（TSP）是关于在一组节点（城市）之间找到最优路线并返回起点的问题。听起来很简单，但对于大量节点来说，通过蛮力法解决是不可能的，因为n个城市的可能排列数是n！。这意味着即使只有30个城市，您需要检查的行程数量也是265,252,859,812,191,058,636,308,480,000,000。对于大型TSP问题，通过蛮力法解决是不切实际的，即使是强大的计算机也无法解决。

幸运的是，已经开发出一些算法，可以大大减少解决大型TSP所需的计算量。其中一款软件Concorde几十年前就为学术界开发了。虽然作为独立软件使用起来相当复杂，且只面向专业人士，但pyconcorde已经开发为Concorde的Python包装器。本文不涉及Concorde使用的算法解释。但我们将深入探讨在Python中复现这些问题及其解决方案所需的代码。

现实世界的地理TSP

如何解决现实世界的地理旅行推销员问题？现实世界的点并不像上面的图片那样通过简单的2D线条相连。相反，地理特征通过各种可能的路线相连，并且这些路线会根据步行、骑行或驾车的方式而变化。

为什么数据科学家或软件工程师要解决现实世界的TSP问题？以下是一些使用案例示例：

1. 一家雇用快递员的公司需要一种计算在城市中最小化每个司机在道路上花费时间的最优路线的方法。
2. 旅游运营商需要找到在约束时间内连接一组目的地的最短路线。
3. 废物处理公司或地方政府需要安排资源，以确保拾取物品的顺序尽可能高效。

为了解决现实世界的TSP问题，可以使用routingpy库来找到地理点（以[经度，纬度]对）之间的路线、距离（以米为单位）和持续时间（以秒为单位）。在本文中，我们将描述解决此类问题的方法。

编码演示

这里概述了使用Python解决地理TSP的指南。解决问题的过程的大致结构如下：

1. 获取一个以[经度，纬度]对形式的n个坐标的列表。
2. 使用路由服务获取每个坐标之间的现实世界持续时间的矩阵（n x n），适用于相应的配置文件（步行、骑行、驾驶汽车、驾驶大型货车等）。该矩阵将是非对称的（从A到B的驾车时间不是从B到A的完全相反）。
3. 将（n x n）矩阵转换为对称矩阵（2n x 2n）。
4. 将该矩阵提供给Concorde求解器，以找到坐标的最优排序。
5. 使用路由服务创建现实世界的路线。
6. 在地图上可视化结果。
7. 可选地，创建最终路线的GPX文件。

每个步骤都将详细介绍。

步骤1：获取坐标

对于我们的例子，我们将考虑在英国的79个城市之间驾驶汽车的问题。下图显示了英国城市的地图（蓝色）。数据科学家可以通过多种方式找到坐标。如果需要，他们可以使用Google Maps或Google Earth手动查找坐标。

此示例中使用的代码结构和数据也可在此GitHub存储库中获得。

以下是包含城市坐标（在存储库中的gb_cities.csv文件）的CSV文件，以及使用pandas导入它所需的代码。

Place Name,Latitude,LongitudeAberdeen,57.149651,-2.099075Ayr,55.458565,-4.629179Basildon,51.572376,0.470009Bath,51.380001,-2.36Bedford,52.136436,-0.460739...

import pandas as pd
df = pd.read_csv('gb_cities.csv')
coordinates = df[['Longitude', 'Latitude']].values
names = df['Place Name'].values

步骤2：使用路由服务获取持续时间矩阵

通过routingpy库提供了几个路由服务。Graphhopper的API包括一个免费的层，允许限制使用。通过routingpy可用的其他路由器在文档中列出。

import routingpy as rp
import numpy as np
api_key = #在https://www.graphhopper.com/api获得免费密钥
api = rp.Graphhopper(api_key=api_key)
matrix = api.matrix(locations=coordinates, profile='car')
durations = np.matrix(matrix.durations)
print(durations)

这是durations，一个79 x 79的矩阵，表示坐标之间的驾驶时间（以秒为单位）：

matrix([[    0, 10902, 30375, ..., 23380, 25233, 19845],        [10901,     0, 23625, ..., 16458, 18312, 13095],        [30329, 23543,     0, ...,  8835,  9441, 12260],        ...,        [23397, 16446,  9007, ...,     0,  2789,  7924],        [25275, 18324,  9654, ...,  2857,     0,  9625],        [19857, 13071, 12340, ...,  8002,  9632,     0]])

可以通过以下方式确定城市之间的驾车时间：

1. 每行和每列对应一个城市：Aberdeen是第一行和第一列，Ayr是第二行和第二列，Basildon是第三行和第三列，依此类推。
2. 要找到Aberdeen和Ayr之间的时间，查看第1行，第2列：10,902秒。反向时间（Ayr到Aberdeen）为10,901秒。
3. 一般来说，第i个城市到第j个城市的时间位于第i行和第j列的交叉点。

请注意，矩阵中沿对角线有零值，因为每个点与自身的距离或持续时间为零。此外，矩阵不完全对称：由于不同的道路布局和交通热点，城市之间的驾车时间在相反方向上可能不相同。不过，它们大致相似，这是可以预期的。

步骤3：将非对称矩阵转换为对称矩阵

在使用此矩阵生成pyconcorde中的最佳顺序之前，我们需要使矩阵对称。Jonker和Volgenant（1983）描述了一种将非对称TSP转换为对称TSP的方法：Transforming asymmetric into symmetric traveling salesman problems, Operations Research Letters, 2(4), 161–163。以下是此转换的理论。如果需要，可以跳过此部分（滚动到名为Transforming the geographic asymmetric TSP的部分）。

Jonker/Volgenant非对称到对称转换

下面是一个有3个节点的非对称TSP以及它的距离矩阵的可视化。

matrix([[0, 5, 2],        [7, 0, 4],        [3, 4, 0]])

这是将其转换为对称TSP所使用的方法的草图：

1. 创建新的虚拟节点A’、B’和C’，用距离零将A连接到A’，B连接到B’，C连接到C’。
2. 将节点连接的权重表示如下：A到B现在由A’到B表示；B到A现在由B’到A表示。B到C现在由B’到C表示；C到B现在由C’到B表示。C到A现在由C’到A表示；A到C现在由A’到C表示。
3. 将所有其他边的权重设置为无穷大，以防止任何算法尝试在它们之间进行旅行。由于在使用pyconcorde时这将是不实际的，因此将所有其他权重设置为比我们拥有的最高权重更高的值。在这种情况下，我们将它们设置为99。

下面是结果距离矩阵。矩阵中节点的排序为：A、B、C、A’、B’、C’。

matrix([[ 0, 99, 99,  0,  7,  3],        [99,  0, 99,  5,  0,  4],        [99, 99,  0,  2,  4,  0],        [ 0,  5,  2,  0, 99, 99],        [ 7,  0,  4, 99,  0, 99],        [ 3,  4,  0, 99, 99,  0]])

请注意，对角线上的值仍然为零，这是符合预期的，而且矩阵现在是对称的。原始矩阵位于新矩阵的左下角，其转置位于右上角。与此同时，左上角和右下角的部分包含节点之间的非常大的权重。

A、B和C（左上角）不再相互连接（严格来说，它们是连接的，但权重非常大，而不是无穷大，为实际目的）。这意味着任何算法都不会寻求在这些节点之间找到路径。同样，A’、B’和C’（右下角）也不相互连接。相反，原始非对称网络的方向性在此处由原始节点A、B和C以及它们的虚拟节点A’、B’和C’上的权重表示。

原始非对称问题的解与新的对称TSP之间存在一对一的映射：

• A — B — C — A 对应于 A — A’ — B — B’ — C — C’ — A
• A — C — B — A 对应于 A — A’ — C — C’ — B — B’ — A

在每种情况下，虚拟节点A’、B’和C’与原始节点A、B和C交替出现，并且每个原始节点与其“伙伴”虚拟节点相邻（A与A’相邻，依此类推）。

转换地理非对称TSP

回到我们的实际例子。我们可以创建一个将非对称TSP矩阵转换为对称矩阵的函数：

def symmetricize(m, high_int=None):        # 如果未提供high_int，则将其设置为最大值的10倍：    if high_int is None:        high_int = round(10*m.max())            m_bar = m.copy()    np.fill_diagonal(m_bar, 0)    u = np.matrix(np.ones(m.shape) * high_int)    np.fill_diagonal(u, 0)    m_symm_top = np.concatenate((u, np.transpose(m_bar)), axis=1)    m_symm_bottom = np.concatenate((m_bar, u), axis=1)    m_symm = np.concatenate((m_symm_top, m_symm_bottom), axis=0)        return m_symm.astype(int) # Concorde需要整数权重

symmetricize(durations) 返回：

matrix([[     0, 461120, 461120, ...,  23397,  25275,  19857],        [461120,      0, 461120, ...,  16446,  18324,  13071],        [461120, 461120,      0, ...,   9007,   9654,  12340],        ...,        [ 23397,  16446,   9007, ...,      0, 461120, 461120],        [ 25275,  18324,   9654, ..., 461120,      0, 461120],        [ 19857,  13071,  12340, ..., 461120, 461120,      0]])

这个 158 x 158 的矩阵在左下角包含了 durations 的一个副本，并在右上角包含了一个转置的副本。461,120 的高值（durations 中最大值的十倍）表示，从实际意义上讲，持续时间为此值的节点是不连接的。

这个矩阵最终可以被传入 pyconcorde 以计算最优路径。

步骤 4：使用 Concorde 求解器

安装 pyconcorde

运行以下命令以安装 pyconcorde（安装在 Linux 或 Mac OS 上可行，但目前不支持 Windows）：

virtualenv venv                                  # 创建虚拟环境source venv/bin/activate                         # 激活虚拟环境git clone https://github.com/jvkersch/pyconcorde # 克隆 git 仓库cd pyconcorde                                    # 切换目录pip install -e .                                 # 安装 pyconcorde

在 Python 中求解 TSP

现在我们可以在 Python 脚本中从 concorde 中导入。

from concorde.problem import Problemfrom concorde.concorde import Concordedef solve_concorde(matrix):    problem = Problem.from_matrix(matrix)    solver = Concorde()    solution = solver.solve(problem)    print(f'Optimal tour: {solution.tour}')    return solution

我们可以将对称的持续时间矩阵传入 solve_concorde()。

durations_symm = symmetricize(durations)solution = solve_concorde(durations_symm)

这里是打印输出：

Optimal tour: [0, 79, 22, 101, 25, 104, 48, 127, 68, 147, 23, 102, 58, 137, 7, 86, 39, 118, 73, 152, 78, 157, 36, 115, 42, 121, 62, 141, 16, 95, 20, 99, 51, 130, 40, 119, 19, 98, 59, 138, 50, 129, 54, 133, 27, 106, 10, 89, 4, 83, 66, 145, 33, 112, 14, 93, 2, 81, 45, 124, 32, 111, 11, 90, 29, 108, 34, 113, 24, 103, 8, 87, 17, 96, 56, 135, 64, 143, 61, 140, 75, 154, 52, 131, 71, 150, 18, 97, 3, 82, 9, 88, 74, 153, 55, 134, 72, 151, 28, 107, 12, 91, 70, 149, 65, 144, 35, 114, 31, 110, 77, 156, 63, 142, 41, 120, 69, 148, 6, 85, 76, 155, 67, 146, 15, 94, 44, 123, 47, 126, 60, 139, 57, 136, 38, 117, 13, 92, 5, 84, 43, 122, 49, 128, 46, 125, 21, 100, 1, 80, 30, 109, 53, 132, 37, 116, 26, 105]

这个解决方案展示了最优路径中节点的排序。注意，这个解决方案如上所述，包含原始节点（编号从0到78）与它们的伙伴虚拟节点（79到157）交替出现：

• 0与79配对，
• 22与101配对，
• 25与104配对，等等…

这表明解决方案已经正确工作。

第五步：创建现实世界的路径

下一步是选择解决方案中的交替元素（对应于原始的79个城市的节点），然后根据顺序排列坐标。

# 选择交替元素：这些对应于原始的tour = solution.tour[::2]# 按照原始坐标和名称排序coords_ordered = [coordinates[i].tolist() for i in tour]names_ordered = [names[i] for i in tour]

下面是names_ordered中的前几个城市名称（最优路径中城市的真实排序）：

['阿伯丁', '邓迪', '爱丁堡', '纽卡斯尔', '桑德兰', '达勒姆', ...]

现在我们将第一个城市添加回来，以形成完整的循环路径，并最终使用Graphhopper方向API获得最终路线。

# 添加第一个城市以形成完整的循环coords_ordered_return = coords_ordered + [coords_ordered[0]]# 为排序后的循环获取完整的行车路线directions = api.directions(locations=coords_ordered_return, profile='car')

第六步：在地图上可视化

在地图上看到最终路线将使我们对结果有信心，同时也能让我们在实际环境中使用这个解决方案。以下代码将显示一个folium地图，可以保存为HTML文件。

import foliumdef generate_map(coordinates, names, directions):    # folium需要纬度和经度    coordinates = [(y, x) for (x, y) in coordinates]    route_points = [(y, x) for (x, y) in directions.geometry]    lat_centre = np.mean([x for (x, y) in coordinates])    lon_centre = np.mean([y for (x, y) in coordinates])    centre = lat_centre, lon_centre    m = folium.Map(location=centre, zoom_start=1, zoom_control=False)    # 绘制路径线    folium.PolyLine(route_points, color='red', weight=2).add_to(m)        # 绘制每个点并添加悬停提示      for i, (point, name) in enumerate(zip(coordinates, names)):        folium.CircleMarker(location=point,                      tooltip=f'{i}: {name}',                      radius=2).add_to(m)    custom_tile_layer = folium.TileLayer(        tiles='http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png',        attr='CartoDB Positron',        name='Positron',        overlay=True,        control=True,        opacity=0.7  # 调整不透明度以控制灰度级别    )    custom_tile_layer.add_to(m)    folium.LayerControl().add_to(m)    sw = (np.min([x for (x, y) in coordinates]), np.min([y for (x, y) in coordinates]))    ne = (np.max([x for (x, y) in coordinates]), np.max([y for (x, y) in coordinates]))    m.fit_bounds([sw, ne])    return mgenerate_map(coords_ordered, names_ordered, directions).save('gb_cities.html')

结果显示在本文顶部。点击此处以查看交互式地图。可以放大地图以查看更多细节，并悬停在各个城市上，可以看到它们在路径序列中的编号。下面是地图的一部分放大显示，显示路线经过谢菲尔德（位于林肯和切斯特菲尔德之间的最优路径上）。

步骤7：可选：创建GPX文件

如果需要在现实生活中遵循计算出的路线，例如在带有GPS的设备上（如手机或汽车导航系统），可以创建GPX文件。这不是优化问题的一部分，但如果您想将路线保存到文件中，则可以选择执行此额外步骤。GPX文件是从directions变量创建的：

def generate_gpx_file(directions, filename):    gpx_template = """<?xml version="1.0" encoding="UTF-8"?>    <gpx version="1.1" xmlns="http://www.topografix.com/GPX/1/1"        xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"        xsi:schemaLocation="http://www.topografix.com/GPX/1/1        http://www.topografix.com/GPX/1/1/gpx.xsd">        <trk>            <name>Track</name>            <trkseg>{}</trkseg>        </trk>    </gpx>    """    trkseg_template = """        <trkpt lat="{}" lon="{}"/>    """    trkseg_elements = ""    for point in directions.geometry:        trkseg_elements += trkseg_template.format(point[1], point[0])    gpx_data = gpx_template.format(trkseg_elements)    with open(filename, 'w') as file:        file.write(gpx_data)generate_gpx_file(directions, 'gb_cities.gpx')

可以从此处下载该问题的GPX文件。

结论

我们已经看到了如何将以下元素结合起来解决现实世界中的地理旅行推销员问题：

1. 来自routingpy库的方向和持续时间矩阵，指定适当的profile（交通方式）。
2. 通过pyconcorde封装器提供的高效而强大的Concorde求解器，以提供最佳路线。
3. 使用folium进行可视化，创建地图。

上面显示的驾驶旅行是解决79个城市旅行推销员问题的令人信服的解决方案，并且根据Concorde求解器是可证明的“最优解”。然而，由于我们使用的是真实世界的数据，最终结果只能与输入数据的质量相当。我们依赖于从routingpy获得的点对点持续时间矩阵是否代表真实世界。实际上，步行、骑自行车或驾车两点之间所需的时间将取决于一天中的时间或一周中的日期。这是我们使用的方法的局限性。增加对最终结果更有信心的一种方法是尝试使用其他路由服务执行相同的方法。每个路由服务（Graphhopper、ORS、Valhalla等）都有自己的API，可用于解决诸如此处所述的TSP问题，并且可以比较不同服务的结果。

尽管解决这样一个问题存在现实世界的限制，但上述方法提供了一个良好的起点，供销售人员或快递员在城市中以尽可能高效的方式前往目的地，或者供希望在旅行中尽可能多地参观景点的游客使用。通过在交互式地图上可视化结果并将路线存储为GPX文件，解决方案对最终用户有用，而不仅仅是实现代码的数据科学家。