前言

大家好，我是吃貨，會寫代碼的那種。

中國最地道的美食，我相信榜單很多。

為什么中餐在國際上不夠有名，因為美食太多了，很難選出能讓大家信服的最最代表性的，以至于最后勝出并盛行的竟是“炒面”。

但是，刀削面，小面，蘭州拉面，擔(dān)擔(dān)面，油潑面。。。哪個會服？

中國農(nóng)業(yè)部推薦了中國最地道的100道鄉(xiāng)間美食，我是信服的。

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

除港澳臺的31省市的地道美食地圖如下：

（圖片來自人民日報）

作為一個吃貨，肯定要想吃完這些美食。

作為一個會寫代碼的吃貨，就要思考一個問題: 我們要游歷全國，吃遍這些美食，最優(yōu)的路線是什么？

這就是一個（吃貨）旅行商問題TSP。

TSP問題

TSP問題是英文縮寫是Travelling Tasting Salesman Problem（旅行商吃貨問題）。

旅行商問題就是“給定一系列城市和每對城市之間的距離，求解訪問每一座城市一次并回到起始城市的最短回路”。

已知TSP算法最壞情況下的時間復(fù)雜度隨著城市數(shù)量的增多而成超多項式（可能是指數(shù)）級別增長。

一般來說，對于上萬個城市數(shù)量的TSP問題，常見的求解器都是可以搞定的。

TSP應(yīng)該在工業(yè)上除了明顯的物流場景，還有一些有趣的場景：比如PCB電路板的焊接點順序。

我們需要分解一下這個有趣問題的任務(wù) （源碼和數(shù)據(jù)見文末）：

如何獲取城市以及美食清單？如何獲取城市距離？這里需要分情況，坐飛機(jī)去還是開車？如何求解這個旅行商問題，得到最后的路線圖？如何畫出路線圖？

獲取美食清單

獲取數(shù)據(jù)，大家想到的第一方法一定是爬蟲，農(nóng)業(yè)部的頁面是htm，屬于入門爬蟲級別。但是不要過于依賴爬蟲技術(shù)。因為有些時候?qū)懘a是效率最差的，尤其是當(dāng)你爬美食網(wǎng)址，容易走神。

htm網(wǎng)址，其實只要右鍵保存，所有的圖片自動就會給你保存到文件夾了。還需要爬嗎？

無非是文字需要解析而已，我們要從網(wǎng)頁里面提取到美食的清單和對應(yīng)的城市以及圖片。

準(zhǔn)備城市數(shù)據(jù)

人民網(wǎng)推薦的榜單很全面的覆蓋了除港澳臺的所有省市，我們可以使用分享的數(shù)據(jù)集來獲取經(jīng)緯度信息。

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

我們需要對數(shù)據(jù)進(jìn)行聚合，因為每個城市可能有好幾道名菜，我們需要對城市進(jìn)行g(shù)roup，然后針對不同的列進(jìn)行不同的統(tǒng)計聚合。

比如：經(jīng)緯度信息，取第一個值即可id：我們可以計數(shù)，用于統(tǒng)計每個城市的榜上有名的菜個數(shù)。

另外，我們可以將每個城市的菜清單整理成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作為經(jīng)典的優(yōu)化問題，算法已經(jīng)很成熟了，這里我們直接調(diào)用OR-Tools進(jìn)行求解。

OR-Tools求解TSP問題需要兩步。

第一步：準(zhǔn)備好城市距離數(shù)據(jù)

我們提供每兩個城市的距離，這里有31個省市，最后就形成[31x31]的矩陣。

我們已經(jīng)獲取了每個城市經(jīng)緯度信息，計算飛行距離就比較容易了。經(jīng)緯度能確定一個點，兩個點之間的距離是一個很好的量化指標(biāo)。距離其實有很多計算方式，我們默認(rèn)的“距離”一般指的是歐式距離，對于小范圍的距離計算問題，比如一個城市內(nèi)，可能歐式距離也可以求得比較理想的結(jié)果。但是對于大范圍距離，就會造成很大的誤差。因為：地球是圓的，球面上的兩點不適用歐式距離。飛行距離就近似于地球大圓距離，可以采用haversine_distances進(jìn)行計算。

還有一個數(shù)據(jù)就是定義出發(fā)城市，這里我隨手定義一個內(nèi)蒙古。

為什么可以隨手呢？

因為最后解決的路徑是一個閉環(huán)，從閉環(huán)的任意一點出發(fā)都可以。

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

第二步：原封不動的復(fù)制官方定義的函數(shù)

OR-Tools的案例很好，以至于這些函數(shù)都不用改一行代碼，只需要傳入data，得到solution變量即可。

這里對函數(shù)進(jìn)行說明：

RoutingIndexManager 函數(shù)就是創(chuàng)建路徑索引，無非就是給每個城市按照順序編個號碼而已。RoutingModel：創(chuàng)建路由模型，也就是實例化一個模型distance_callback：告訴模型如何計算使用城市的距離信息，采用RegisterTransitCallback函數(shù)傳入距離信息到模型。這樣模型就可以計算總距離了SetArcCostEvaluatorOfAllVehicles：計算其它成本，這里我們默認(rèn)距離就是成本SolveWithParameters：開始計算，很快就出結(jié)果了，畢竟才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公里，路線也已經(jīng)規(guī)劃好了。

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，畫出路線的基本思路就是將相鄰的城市依次連接起來。從圖中我們可以看出：

多數(shù)情況下，路徑就是從一個省到鄰省。從江蘇開始后直飛東三省，之后再飛回上海。另外一個出人意料的地方就是從云南直飛新疆，而不是西藏，這點與我們的直覺不太吻合。真相就是：我們采用的是大圓距離，而且地圖上會有麥卡托投影的視覺誤差。

自駕版本

如果想要自駕過去，我們同樣可以采用這樣的方式進(jìn)行求解。

自駕版本的計算要稍微復(fù)雜一些，因為需要獲取城市之間的實時路線規(guī)劃，然后進(jìn)行求解。

一般我會采用openstreetmap來進(jìn)行規(guī)劃，因為它的API是免費的。

我們這次再看一下自駕版本結(jié)果：

吃完老北京炸醬面（北京）之后，我們應(yīng)該直奔東北了，取嘗嘗殺豬菜了（黑龍江）吃飯鹽水鴨（江蘇）之后，之后的確是應(yīng)該去嘗蟹殼黃了(上海)這次擔(dān)擔(dān)面（四川）離風(fēng)干牛羊肉（西藏）更近了。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 公里，接近飛行距離的兩倍，不過自駕的風(fēng)景應(yīng)該是飛行風(fēng)景的N倍吧。

總結(jié)

目前，我在刀削面，肉夾饃，火鍋，松鼠桂魚之間游蕩，希望有朝一日自駕嘗遍中國美食，Bite of China。

冷鏈服務(wù)業(yè)務(wù)聯(lián)系電話：13613841283

標(biāo)簽:

食品安全網(wǎng) ：https://www.food12331.com

上一篇：一百種家常菜做法

下一篇：會被秒贊的100條美食文案