吃遍中国最地道的100道美食（权威榜单）的最佳路线（TSP问题）

时间：2023-04-04 02:20:42来源：food栏目：餐饮美食新闻阅读：

前言

大家好，我是吃货，会写代码的那种。

中国最地道的美食，我相信榜单很多。

为什么中餐在国际上不够有名，因为美食太多了，很难选出能让大家信服的最最代表性的，以至于最后胜出并盛行的竟是“炒面”。

但是，刀削面，小面，兰州拉面，担担面，油泼面。。。哪个会服？

中国农业部推荐了中国最地道的100道乡间美食，我是信服的。

美食链接：http://www.moa.gov.cn/ztzl/zgnmfsj/xwzx/201809/t20180913_6157258.htm

除港澳台的31省市的地道美食地图如下：

（图片来自人民日报）

作为一个吃货，肯定要想吃完这些美食。

作为一个会写代码的吃货，就要思考一个问题: 我们要游历全国，吃遍这些美食，最优的路线是什么？

这就是一个（吃货）旅行商问题TSP。

TSP问题

TSP问题是英文缩写是Travelling Tasting Salesman Problem（旅行商吃货问题）。

旅行商问题就是“给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路”。

已知TSP算法最坏情况下的时间复杂度随着城市数量的增多而成超多项式（可能是指数）级别增长。

一般来说，对于上万个城市数量的TSP问题，常见的求解器都是可以搞定的。

TSP应该在工业上除了明显的物流场景，还有一些有趣的场景：比如PCB电路板的焊接点顺序。

我们需要分解一下这个有趣问题的任务（源码和数据见文末）：

如何获取城市以及美食清单？如何获取城市距离？这里需要分情况，坐飞机去还是开车？如何求解这个旅行商问题，得到最后的路线图？如何画出路线图？

获取美食清单

获取数据，大家想到的第一方法一定是爬虫，农业部的页面是htm，属于入门爬虫级别。但是不要过于依赖爬虫技术。因为有些时候写代码是效率最差的，尤其是当你爬美食网址，容易走神。

htm网址，其实只要右键保存，所有的图片自动就会给你保存到文件夹了。还需要爬吗？

无非是文字需要解析而已，我们要从网页里面提取到美食的清单和对应的城市以及图片。

准备城市数据

人民网推荐的榜单很全面的覆盖了除港澳台的所有省市，我们可以使用分享的数据集来获取经纬度信息。

file = ../data/china_cities.csv all_cities_df = pd.read_csv(file,encoding=gb18030) main_cities = all_cities_df[all_cities_df.sort_values(by=行政代码)[行政代码].astype(str).str.find(0100)>0] main_cities.columns = [id,name,lon,lat] main_cities[code] = main_cities[name].str[0:2] meal_city_df = pd.merge(left = meal_df,right=main_cities,on=code,how=left) meal_city_df.head()

我们需要对数据进行聚合，因为每个城市可能有好几道名菜，我们需要对城市进行group，然后针对不同的列进行不同的统计聚合。

比如：经纬度信息，取第一个值即可id：我们可以计数，用于统计每个城市的榜上有名的菜个数。

另外，我们可以将每个城市的菜清单整理成List，这样可以一起显示。

meal_city_df_count = meal_city_df.groupby(code).agg({name:first,lon:first,lat:first,id:count}) meal_list = meal_city_df.groupby(name)[meal].apply(list).to_frame().reset_index() meal_list.columns = [name,meal_list] meal_list[meal_list] = meal_list[meal_list].str.join(",") meal_city_df_count = pd.merge(left = meal_city_df_count,right = meal_list,on=name) meal_city_df_count.tail()

我们可以初步看一下这100道名菜的分布。

求解TSP问题

TSP作为经典的优化问题，算法已经很成熟了，这里我们直接调用OR-Tools进行求解。

OR-Tools求解TSP问题需要两步。

第一步：准备好城市距离数据

我们提供每两个城市的距离，这里有31个省市，最后就形成[31x31]的矩阵。

我们已经获取了每个城市经纬度信息，计算飞行距离就比较容易了。经纬度能确定一个点，两个点之间的距离是一个很好的量化指标。距离其实有很多计算方式，我们默认的“距离”一般指的是欧式距离，对于小范围的距离计算问题，比如一个城市内，可能欧式距离也可以求得比较理想的结果。但是对于大范围距离，就会造成很大的误差。因为：地球是圆的，球面上的两点不适用欧式距离。飞行距离就近似于地球大圆距离，可以采用haversine_distances进行计算。

还有一个数据就是定义出发城市，这里我随手定义一个内蒙古。

为什么可以随手呢？

因为最后解决的路径是一个闭环，从闭环的任意一点出发都可以。

from sklearn.metrics.pairwise import haversine_distances from math import radians def geo_distance(p1,p2): dis = haversine_distances([[radians(_) for _ in p1], [radians(_) for _ in p2]])[0][1]* 6371000/1000 # multiply by Earth radius to get kilometers return dis from ortools.constraint_solver import routing_enums_pb2 from ortools.constraint_solver import pywrapcp meal_city_df_count[loc] = meal_city_df_count.apply(lambda x: list([x[lon], x[lat]]),axis=1) xv,yv= np.meshgrid(meal_city_df_count[loc], meal_city_df_count[loc],indexing=ij) data = {} data[distance_matrix] = [[round(geo_distance(xv[i][j],yv[i][j])) for i in range(xv.shape[0])] for j in range(xv.shape[1])] data[num_vehicles] = 1 data[depot] = 2

第二步：原封不动的复制官方定义的函数

OR-Tools的案例很好，以至于这些函数都不用改一行代码，只需要传入data，得到solution变量即可。

这里对函数进行说明：

RoutingIndexManager 函数就是创建路径索引，无非就是给每个城市按照顺序编个号码而已。RoutingModel：创建路由模型，也就是实例化一个模型distance_callback：告诉模型如何计算使用城市的距离信息，采用RegisterTransitCallback函数传入距离信息到模型。这样模型就可以计算总距离了SetArcCostEvaluatorOfAllVehicles：计算其它成本，这里我们默认距离就是成本SolveWithParameters：开始计算，很快就出结果了，毕竟才31个城市。# Create the routing index manager. manager = pywrapcp.RoutingIndexManager(len(data[distance_matrix]), data[num_vehicles], data[depot]) # Create Routing Model. my_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 data[distance_matrix][from_node][to_node] transit_callback_index = my_routing.RegisterTransitCallback(distance_callback) # Define cost of each arc. my_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) # Solve the problem. solution = my_routing.SolveWithParameters(search_parameters) # Print solution on console. if solution: print_solution(manager, my_routing, solution)

可以看到总的飞行距离10954公里，路线也已经规划好了。

Objective: 10954 km Route for vehicle 0: 2 -> 10 -> 17 -> 16 -> 3 -> 6 -> 9 -> 8 -> 14 -> 26 -> 4 -> 30 -> 0 -> 18 -> 23 -> 15 -> 20 -> 21 -> 11 -> 19 -> 12 -> 25 -> 27 -> 5 -> 1 -> 13 -> 24 -> 29 -> 22 -> 7 -> 28 -> 2

展示路线

这里我采用的plotly，画出路线的基本思路就是将相邻的城市依次连接起来。从图中我们可以看出：

多数情况下，路径就是从一个省到邻省。从江苏开始后直飞东三省，之后再飞回上海。另外一个出人意料的地方就是从云南直飞新疆，而不是西藏，这点与我们的直觉不太吻合。真相就是：我们采用的是大圆距离，而且地图上会有麦卡托投影的视觉误差。

自驾版本

如果想要自驾过去，我们同样可以采用这样的方式进行求解。

自驾版本的计算要稍微复杂一些，因为需要获取城市之间的实时路线规划，然后进行求解。

一般我会采用openstreetmap来进行规划，因为它的API是免费的。

我们这次再看一下自驾版本结果：

吃完老北京炸酱面（北京）之后，我们应该直奔东北了，取尝尝杀猪菜了（黑龙江）吃饭盐水鸭（江苏）之后，之后的确是应该去尝蟹壳黄了(上海)这次担担面（四川）离风干牛羊肉（西藏）更近了。bjective: 19258 km Route for vehicle 0: 2 -> 3 -> 30 -> 4 -> 26 -> 6 -> 9 -> 8 -> 14 -> 0 -> 18 -> 23 -> 15 -> 20 -> 21 -> 11 -> 19 -> 12 -> 1 -> 25 -> 27 -> 5 -> 24 -> 13 -> 29 -> 22 -> 7 -> 28 -> 17 -> 16 -> 10 -> 2

可以看到总的飞行距离19258 公里，接近飞行距离的两倍，不过自驾的风景应该是飞行风景的N倍吧。