【Python】在Python中使用NetworkX解决德国的旅行推销员问题

使用Christofides算法，在访问每个城市一次的同时，可以发现穿越德国16个联邦州首府的最短路线。

我已经在德国居住六年了。德国由16个联邦州组成，到目前为止，我已经访问了五个州的首府。最近，一个想法打动了我：

我想穿越所有16个联邦州的首府，从柏林出发，到柏林结束，每座城市只访问一次。这样做的最短可能路线是什么？

这类问题称为旅行推销员问题（TSP）。顾名思义，旅行推销员会挨家挨户，或在特定地区内的地区之间，尝试销售其业务的产品或服务。

TSP处理的是一次找到覆盖所有点的最短路线，因为这是成本、时间和能量方面最有效的路线。

TSP是一个NP难问题

决策问题是一个有“是”或“否”答案的问题。可能有不同的复杂性类来解决这个问题。例如：

• P是一个复杂性类，其中决策问题可以用确定性图灵机（DTM）在多项式时间内求解。对于一些非负整数k，多项式时间算法的执行时间为O(nᵏ)，其中n是输入的复杂度。

可以用多项式时间算法解决的问题称为可处理问题。此类算法的示例包括线性搜索 (O(n))、二分搜索(O(log n))、插入排序(O(n²))、归并排序 (O(n log n))和矩阵乘法 (O(n³))。

• NP（Non-deterministic Polynomial time，非确定性多项式时间）是指一组可以用非确定性图灵机（NTM）在多项式时间内求解的决策问题。NTM可以有一组动作来解决问题，但每次执行的转换都是不确定的（随机的）（Tutorialspoint，2022）。一旦为NP问题提供了潜在解决方案，DTM就可以在多项式时间内验证其正确性（Viswarupan，2016）。

如果明天有人找到它的解决方案，今天不确定的算法也可能是确定的（Bari，2018年和Maru，2020年）。这意味着DTM在多项式时间内可以解决的问题也可以在多项式时间里由NTM解决。因此，P是NP的子集。

• NP-hard是“至少与NP中最难的问题一样难”的决策问题的复杂性类别。当NP中的每个问题L都可以在多项式时间内简化为H时，问题H是NP难的；也就是假设H的解需要1个单位时间，可在多项式时间内求解L。

TSP是一个NP难问题。这意味着，如果有大量城市，在“合理”的时间内评估每个可能的解决方案是不可行的（Hayes，2019）。

哈密顿循环

要解决旅行推销员问题（TSP），首先需要了解哈密顿循环（也称为哈密顿循环）的概念。

汉密尔顿循环是一个通过图形的图形循环（闭环），每个节点只访问一次（Mathworld，2022a）。这是以威廉·罗文·汉密顿爵士的名字命名的，他通过一个名为哈密顿之谜的游戏（也称为伊科西亚游戏）引入了这个概念。

在一个有n个节点的图中，可能的哈密顿循环总数由（n-1）给出！如果我们把顺时针和逆时针路径看作两条不同的路径。如果我们认为顺时针和逆时针路径相同，则存在（n-1）/2个可能的哈密顿循环。

让我们以四个城市为例：柏林、杜塞尔多夫、汉堡和慕尼黑。我的目标是穿越这四座城市，从柏林出发，到柏林结束，同时穿越这四个城市之间的一次。本例中可能的唯一路径数（哈密顿循环）为（4–1）！=3! = 6.

在下面的代码中，城市是这四个城市的列表，但起点和终点都是柏林。其他三个城市被定义为*rest。这三个城市的排列同时占据了所有三个城市，导致了所有可能的汉密尔顿循环。

这在下面的Python代码中显示。

from itertools import permutations

cities = ["Berlin","Dusseldorf","Hamburg","Munich","Berlin"]
start, *rest, end = cities

#start and end are both Berlin.
#We need to get permutation of rest of the cities considering all of them
#3P3 = 3!/0! = 6 = (4-1)!
paths = [(start, *path, end) for path in permutations(rest,len(rest))]
paths

[('Berlin', 'Dusseldorf', 'Hamburg', 'Munich', 'Berlin'),
 ('Berlin', 'Dusseldorf', 'Munich', 'Hamburg', 'Berlin'),
 ('Berlin', 'Hamburg'


    
, 'Dusseldorf', 'Munich', 'Berlin'),
 ('Berlin', 'Hamburg', 'Munich', 'Dusseldorf', 'Berlin'),
 ('Berlin', 'Munich', 'Dusseldorf', 'Hamburg', 'Berlin'),
 ('Berlin', 'Munich', 'Hamburg', 'Dusseldorf', 'Berlin')]

在下面的要点中，make_tsp_tree函数首先创建一个哈密顿路径列表，然后从这些路径列表中创建一个定向前缀树，然后通过删除根节点和nil节点返回图形对象G。

import networkx as nx
import matplotlib.pyplot as plt

cities = ["Berlin","Dusseldorf","Hamburg","Munich","Berlin"]

def make_tsp_tree(cities):
    """
    Create all 哈密顿 paths from start to end city from a list of cities.
    Creates a directed prefix tree from a list of the created paths. 
    Remove the root node and nil node.
    """
    start, *rest, end = cities
    paths = [(start, *path, end) for path in permutations(rest)]

    #Creates a directed prefix from a list of paths
    G = nx.prefix_tree(paths)

    #remove synthetic root node (None) and nil node (NIL)
    G.remove_nodes_from([0, -1])

    return G

#Get graph object with all 哈密顿 paths
G = make_tsp_tree(cities)


#Create node positions for G using Graphviz.
#Available layouts: https://graphviz.org/docs/layouts/
#dot gives hierarchical or layered drawings of directed graphs
#possible layouts: dot, neato, twopi, circo, fdp, osage, patchwork, sfdp
positions = nx.nx_agraph.graphviz_layout(G, "dot")

plt.figure(figsize = (12,8))


#Draw networkx Graph object with labels off
nx.draw_networkx(G,
                 pos = positions,
                 node_color = "mistyrose",
                 with_labels= False)


#Draw NetworkX Graph object with defined labels for each node
nx.draw_networkx_labels(G,
                        pos = positions,
                        labels = dict(G.nodes(data = "source")))

plt.title(


    
"哈密顿 paths for Berlin, Hamburg, Dusseldorf and Munich.\n Start and end at Berlin.")
plt.axis("off")
plt.show()

这四个城市列表的汉密尔顿路径如下图所示：

在德国地图上对16个州首府进行地理编码和绘图

我将德国16个州首府的名单定义为首府。使用一个称为地理编码的过程，我可以得到所有16个城市的坐标。

本文详细描述了使用geopy包进行地理编码的过程。

capitals = [‘Berlin’, ‘Bremen’, ‘Dresden’, ‘Dusseldorf’,
            ‘Erfurt’, ‘Hamburg’, ‘Hannover’, ‘Kiel’,
            ‘Magdeburg’, ‘Mainz’, ‘Munich’, ‘Potsdam’,   ‘Saarbrucken’, ‘Schwerin’, ‘Stuttgart’, ‘Wiesbaden’]

使用GeoPandas包，我在名为德国的GeoDataFrame中提取了德国地图。

要绘制德国首都，首先我使用NetworkX创建一个图形对象G。然后我创建了16个节点，每个城市一个。

接下来，我为每个城市的位置创建一个节点和坐标字典。我还为标签创建了另一个节点和首都名称字典。代码和生成的地图在下面的截图中给出。

G = nx.Graph()
#Create a graph object with number of nodes same as number of cities
nodes = np.arange(0, len(capitals))
G.add_nodes_from(nodes)

#Create a dictionary of node and coordinate of each state for positions
positions = {node:coordinate for node, coordinate in zip(nodes, coordinates)}

#Create a dictionary of node and capital for labels
labels = {node:capital for node, capital in zip(nodes, capitals)}

fig, ax = plt.subplots(figsize = (10, 7))
germany.plot(color = "whitesmoke", edgecolor = "black", ax = ax)
nx.draw_networkx(G, pos = positions,
                 labels = labels, ax = ax,
                 bbox = dict(facecolor = "skyblue", boxstyle = "round",
                             ec = "black", pad = 0.3),)

plt.title("Map of Germany with capital cities of 16 federal states")
plt.axis("off")
plt.show()

德国任何两个首都之间的所有可能路径

下一步是确定并绘制德国任何两个城市之间的所有可能路径。

由于我们的列表中有16个城市，因此任意两个城市之间可能的路径数如下所示ⁿCᵣ = ¹⁶C₂ = 120.在现有图形对象G中，我添加任意两个节点之间的边，但其本身除外，如下所示：

for i in nodes:
     for j in nodes:
          if i!=j:
               G.add_edge(i, j)

使用德国最近添加的边绘制更新的图形对象G，如下所示：

我创建了H作为G的备份副本，供以后使用。

接下来，我计算德国所有两个城市之间的欧氏距离，该信息存储为G中边的权重。因此，G现在是一个完整的加权图，其中一条边连接着每对图节点，每条边都有一个与其相关联的权重。

    
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)):

        #Multidimensional Euclidean distan
          dist = math.hypot(pos[i][0] - pos[j][0], pos[i][1] - pos[j][1])
          dist = dist
          G.add_edge(i, j, weight=dist)

下表以矩阵形式显示了城市之间的距离。在下表中，如主对角线所示，城市之间的距离为零。主对角线上方的所有条目都反映在对角线下方的相等条目中。因此，它是一个空心（零对角）对称矩阵的示例。

求解TSP的算法

要查找解决方案循环，NetworkX对无向图使用一种默认算法，称为Christofides算法，作为解决方案循环一部分的边列表存储为edge_list。这在下面的代码中实现：

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

Christofides算法

Christofides算法在距离形成度量空间的情况下（它们是对称的，并遵循三角形不等式，即∆ABC，a+b≥c）（古德里奇和塔马西亚，2014年）。该算法以塞浦路斯数学家Nicos Christofides命名，他于1976年设计了该算法。到目前为止，该算法提供了求解TSP的最佳近似比（Christofites，2022）。

该算法的基本步骤如下。这些步骤已在此处的NetworkX源代码中实现。我还举例说明了我们问题的算法的逐步实现。

1. 求G的最小（加权）生成树T。

加权、连通、无向图（此处为G）的最小生成树是一个由边的子集组成的图，这些边将所有连接节点（此处为16）连接在一起，同时最小化边上的权重总和。

NetworkX使用Kruskal算法来寻找最小生成树（NetworkX，2015）。德国16个联邦州的首府城市情况如下：

2. 在T中生成一组奇数度为O的节点。

在下一步中，将创建一组名为O的节点，其中每个节点的阶数都是奇数。从上面的树T中，我们可以看到柏林、波茨坦、斯图加特和施韦林都有2级，即偶数。因此，将在新的集合O中删除这些节点，如下所示：

3. 在由O的顶点给出的诱导子图中找到一个最小权完美匹配M。

要实现完美匹配，图中需要有偶数个节点，每个节点都与另一个节点相连。因此，完美匹配是包含n/2条边（最大可能）的匹配（Mathworld，2022b）。

最小权重完美匹配M计算完美匹配，以最小化匹配边的权重，如下图所示。

4. 合并M和T的边，形成连通多重图H。

在下一步中，该算法将步骤1中T的边与步骤3中M的边合并，形成一个连通多重图H。

5. 使用M和T的边缘构建欧拉回路。

在下一步中，Christofides算法使用M和T的边构建欧拉回路。欧拉回路是有限图中的一条轨迹，每条边只访问一次。欧拉回路中的每个节点都必须具有偶数度。哈密顿循环和欧拉回路之间的区别在于，哈密顿周期只通过每个节点一次，并在初始节点处结束，而欧拉回路只通过每个边一次，然后在初始节点结束。

6. 通过跳过重复的节点，将欧拉回路转换为汉密尔顿循环。

就TSP而言，可能有一些节点在欧拉回路中重复出现。此类节点的度超过2，即节点访问次数超过一次。Christofides算法利用一个快捷函数从欧拉回路中删除重复的节点并创建一个循环。因此，TSP的解决方案是通过使用循环内节点之间的连续边列表来实现的。

Christofides算法产生的解的权重在最优解的2/3以内。

TSP解决方案直接使用NetworkX实施

使用如上所述的Christofides算法，按周期给出TSP的解。由于cycle只包含首都城市的指数，因此我得到了表示城市顺序的解决方案路径tsp_cycle，如下面的截图所示。

如下面截图所示，我绘制了德国地图，其中任何两个城市之间的所有可能路径都用蓝线表示。红线表示edge_list。

#Create a dictionary of node and capital for labels
labels = {node:capital for node, capital in zip(nodes, capitals)}
tsp_cycle = [labels[value] for value in cycle]
print (tsp_cycle)

['Berlin', 'Dresden', 'Munich', 'Stuttgart', 'Mainz', 'Wiesbaden', 'Saarbrucken', 'Dusseldorf', 'Hannover', 'Bremen', 'Hamburg', 'Kiel', 'Schwerin', 'Magdeburg', 'Erfurt', 'Potsdam', 'Berlin']

fig, ax = plt.subplots(figsize = (12, 8))
germany.plot(color = "whitesmoke", edgecolor = "black", ax = ax)

#All possible routes
nx.draw_networkx_edges(H, pos = positions, edge_color="blue", width=0.5, ax = ax)

#TSP solution route
nx.draw_networkx(G, pos = positions, labels = labels, node_size=0,
                 edgelist=edge_list, edge_color="red", width = 3, ax = ax)

plt.axis("off")
plt.show()

下图通过删除两个城市之间的蓝色（所有可能的）边缘，并仅以红色绘制德国TSP的解决方案，提供了一个更干净的解决方案。

为每条边绘制具有唯一颜色的解决方案

我想通过为城市之间的每条边缘绘制独特的颜色，使我的解决方案更具吸引力。

计算机显示器上使用的任何颜色都可以表示为可见光谱中RGB（红色、绿色和蓝色）的比例（Mathisfun，2021）。因此，在Python中，任何颜色都可以表示为#RRGGBB，其中RR、GG和BB具有从00到ff的十六进制值。#000000表示白色，#ffffff表示黑色。

我创建了一个名为get_colors（n）的函数，它基于十六进制系统中#RRGGBB的随机值创建一个随机颜色列表。

import random
get_colors = lambda n: list(map(lambda i: “#” +”%06x” % random.randint(0, 0xffffff),range(n)))

在下面的代码中，我传递get_colors（16）作为edge_color。

fig, ax = plt.subplots(figsize = (20, 12))

germany.plot(ax = ax, color = “whitesmoke”, edgecolor = “black”)

# Draw the route
nx.draw_networkx(G, pos = positions, labels = labels, 
 edgelist=edge_list, edge_color=get_colors(16), width=3,
 node_color = “snow”, node_shape = “s”, node_size = 300,
 bbox = dict(facecolor = “#ffffff”, boxstyle = “round”,ec = “silver”),
 ax = ax)

plt.title(“Travelling Salesman Problem Solution for Germany”, fontsize = 15)
plt.axis(“off”)

因此，我得到了下面的图，它比前一个图色彩丰富，更有吸引力，并且解决方案中每个边都有独特的颜色。

使用folium绘制解决方案

上面的图也可以使用folium包绘制。

为此，我将坐标修改为folium_coordinates，循环为route。这些是相同的数据，但采用了新的格式：一个与folium兼容的GPS坐标（纬度和经度）列表。

结果图如下所示：

import folium

folium_coordinates = []
for x,y in coordinates:
    folium_coordinates.append([y,x])

route = []
for stop in cycle:
    route.append(folium_coordinates[stop])

m1 = folium.Map(location = [51, 10],   #latitude (N), longitude (E)
                tiles = "OpenStreetMap", 
                zoom_start= 6
                )

for coordinate, capital in zip(folium_coordinates, capitals):
    folium.Marker(location = coordinate,
                 popup = capital).add_to(m1)

folium.PolyLine(route).add_to(m1)

m1.save("germany tsp.html")

m1

结论

这篇帖子的问题陈述是寻找穿越德国16个联邦州首府城市的最短路线，从柏林开始，到柏林结束，其间每座城市只访问一次。

我首先描述了解决任何决策问题的不同复杂性类：P（多项式时间）、NP（非确定性多项式时间）和NP难。接下来，我讨论了哈密顿循环的概念。

我使用地理编码找到了德国16个联邦州首府的坐标，并使用GeoPandas软件包在德国地图上绘制了它们。我在地图上添加了任意两个城市之间的所有边缘。

接下来，我演示了Christofides算法如何通过逐步实现为旅行推销员问题提供解决方案。最后，我使用NetworkX、GeoPandas和Matplotlib包为这个问题绘制了一个干净的解决方案。我还为每一条边绘制了一种独特的颜色的解决方案，并使用folium进行绘图。

感谢阅读！

参考引用

• Bari, 2018. NP-Hard and NP-Complete Problems.

• Christofides, 2022. Worst-Case Analysis of a New Heuristics for the Travelling Salesman Problem.

• Goodrich and Tamassia, 2014. Algorithm Design and Applications | Wiley

• Hayes, 2019. Solving Travelling Salesperson Problems with Python.

• Maru, 2020. P, NP, NP Hard and NP Complete Problem | Reduction | NP Hard and NP Compete | Polynomial Class.

• Matuszek, 1996. Polynomial-time algorithm.

• Mathisfun.com, 2021. Hexadecimal Colours.

• Mathworld, 2022a. 哈密顿ian Cycle.

• Mathworld, 2022b. Perfect Matching.

• MIT OpenCourseWare, 2016. R9. Approximation Algorithms: Traveling Salesman Problem

• NetworkX, 2015. minimum_spanning_tree.

• Tutorialspoint, 2022. Non-Deterministic Turing Machine.

• Viswarupan, 2016. P vs NP Problem.

    
往期


    
精彩回顾






    

适合初学者入门人工智能的路线及资料下载
(图文+视频)机器学习入门系列下载
机器学习及深度学习笔记等资料打印
《统计学习方法》的代码复现专辑
机器学习交流qq群955171419，加入微信群请扫码