吃遍中國最地道的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倍吧。