漢米爾頓路徑問題 (III) - 使用 Google OR-Tools 實作 === ###### tags: `TSP` `Hamiltonian Path` `Operations Research` `OR-Tools` `Python` ## 前言 [__上一篇文章__](https://hackmd.io/kdviO126QJq6B4YfxrGzwQ) 說明了最短漢米爾頓路徑問題的數學模式，這篇則將說明如何撰寫程式，把最短漢米爾頓路徑問題實作出來。 本文實作將使用 Python 程式語言，並使用 Google 開發的 OR-Tools 套件求解。本文將提供三種方式求解最短漢米爾頓路徑。 ## 問題描述 小乖是一個咖哩愛好者，有一天他搭高鐵到台北，決定要來品嘗一下台北好吃的咖哩。他在清單上列出了 6 家咖哩店，預計從台北車站出發，造訪完清單上所有咖哩店後，最後到大稻埕碼頭貨櫃市集喝酒，結束一日台北咖哩行。因為時間有限，小乖希望可以用最短旅行時間跑完上述的行程。 上述問題即是最短漢米爾頓路徑問題，給定起點 (台北車站) 和訖點 (大稻埕碼頭貨櫃市集)，求出一條經過所有節點且為距離或旅行時間為最短的路徑，也就是： ``` 台北車站 -> ? -> ? -> ? -> ? -> ? -> ? -> 大稻埕碼頭貨櫃市集 ``` 我們要做的就是把中間這些點做排序，使這條路徑為最短。 :::danger :warning: __聲明：本文絕對沒有業配！沒有業配！沒有業配！__ (但如果要找我業配也是可以^^) ::: 以下為起訖點及 6 家咖哩店的經緯度座標： |編號|地點|經度|緯度| |:--:|:--:|:--:|:--:| |0|台北車站 (起點)|121.51705498294422|25.04812699185495| |1|溫咖哩|121.52694765430479|25.04160674979448| |2|哞屋|121.52294339950605|25.064214083600284| |3|勝利洋食|121.54956207412458|25.043419518468635| |4|天菜咖哩|121.55331801759017|25.043348130259034| |5|島津和牛咖哩民生店|121.55186587878282|25.05846636112716| |6|新高軒信義店(A11-B2)|121.5672858678217|25.03670073623688| |7|大稻埕碼頭貨櫃市集 (訖點)|121.50742317199105|25.0561805065471| <br> 用地圖呈現點位分布如下圖，粉紅色表示起點，橘色表示終點： ![](https://hackmd.io/_uploads/HJNPrF19p.png) <br> ## 取得旅行時間矩陣 這裡我們透過 [Open Street Routing Machine (OSRM)](https://project-osrm.org/) 取得旅行時間矩陣，可以利用筆者開發的 [OSRM API](https://github.com/chhuang17/osrm_api) 取得相關資料。 ```python= from osrm_api import osrm startPoint = (121.51705498294422, 25.04812699185495) # 台北車站 endPoint = (121.50742317199105, 25.0561805065471) # 大稻埕碼頭貨櫃市集 # Curry shop curryShops = [ (121.52694765430479, 25.04160674979448), # 溫咖哩 (121.52294339950605, 25.064214083600284), # 哞屋Mon wo 焼き肉丼 、咖哩カレー (121.54956207412458, 25.043419518468635), # 勝利洋食 (121.55331801759017, 25.043348130259034), # 天菜咖哩 (121.55186587878282, 25.05846636112716), # 島津和牛咖哩民生店 (121.5672858678217, 25.03670073623688), # 新高軒信義店(A11-B2) ] # 點位經緯度列表 nodeList = [startPoint] + curryShops + [endPoint] # 點位編號 nodeNumbers = [i for i in range(len(nodeList))] # 透過 OSRM API 取得旅行時間矩陣 timeMatx = osrm.odMatrix(nodeList, get='duration', timeUnit='second') print(timeMatx) # Output # [[ 0. 188.7 330.1 330.3 375.1 510.2 492.9 261.1] # [171.2 0. 359.2 259.3 304.1 450.8 421.9 322.4] # [256.4 356.9 0. 494.6 539.4 417. 657.2 299.3] # [354.2 292.3 463.8 0. 122.9 248.8 272.3 496.4] # [368.6 306.7 478.2 61.3 0. 269.2 230. 451.9] # [484.4 465.4 376. 303.2 289.4 0. 408.3 534.3] # [491. 468.3 621.3 259.6 235. 435.9 0. 568.2] # [298.2 409. 575.9 558. 536.8 693.1 640. 0. ]] ``` ## 方法一：使用 MPSolver 求解 我們知道最短漢米爾頓路徑問題是一個整數規劃問題，而 OR-Tools 中的 MPSolver 可以用來解決線性規劃和混合整數規劃問題，因此當然可以用 MPSolver 求解最短漢米爾頓路徑。 ### 求解器和決策變數 首先建立求解器，這裡使用 `'SCIP'` 求解器，同時建立決策變數 $x_{00},\space x_{01},\cdots, x_{77}$。我們利用 `solver.IntVar()` 建立整數決策變數，第一個 args 為決策變數下界，第二個 args 為決策變數上界，第三個 args 為變數名稱。 ```python= from ortools.linear_solver import pywraplp # 使用 SCIP 求解器 solver = pywraplp.Solver.CreateSolver('SCIP') # 建立二元決策變數，xij = {0,1} variables = [solver.IntVar(0, 1, f"x{i:02d}_{j:02d}") for i in range(len(nodeNumbers)) for j in range(len(nodeNumbers))] ``` :::info :bulb: 為了方便 debug，這裡將變數 $x_{01}$ 命名為 `x_00_01`，變數 $x_{46}$ 命名為 `x_04_06`。 ::: 接著建立成本矩陣 $c_{ij}$，也就是旅行時間矩陣，並且對旅行時間矩陣稍作修改，在對角線的地方填入一個很大的數作為懲罰項，避免運算時選到自己。 ```python= c = timeMatx # Add penalty term at cii for i in range(c.shape[0]): c[i][i] = 1e7 ``` ### 目標式 目標式為最小化總成本，這裡的總成本就是總旅行時間，即： $$ \min \space\space \sum_{i}\sum_{j}c_{ij}x_{ij} $$ 首先建立 `objective` 物件，並針對每個變數設定其係數 (也就是 $c_{ij}$ )，最後設定求解目標為最小化。 :::info :bulb: 目標式係數一定要 `int`，不可以是 `float`, `double`, ...。 :bulb: `var` 的資料型態是 ortools 的 Variable, 訪問它的話要先轉換成 str。 ::: ```python= objective = solver.Objective() for var in variables: coeff = int(round(c[int(str(var)[1:3])][int(str(var)[4:6])], 0)) objective.SetCoefficient(var, coeff) objective.SetMinimization() ``` ### 流量限制式 1 - 除了終點，其餘節點流出的總和為 1 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{i}x_{ij} &= 1 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= nodeNumList = [f"{x:02d}" for x in range(len(nodeNumbers)-1)] for node in nodeNumList: constraint = solver.Constraint(1, 1) for var in variables: if str(var)[1:3] == node: constraint.SetCoefficient(var, 1) ``` ### 流量限制式 2 - 除了起點, 其餘節點流入的總和為 1 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{i}x_{ij} &= 1 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= nodeNumList = [f"{x:02d}" for x in range(1, len(nodeNumbers))] for node in nodeNumList: constraint = solver.Constraint(1, 1) for var in variables: if str(var)[4:6] == node: constraint.SetCoefficient(var, 1) ``` ### 流量限制式 3 - 流入起點的總和為 0 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{i}x_{i1} &= 0 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= constraint = solver.Constraint(0, 0) for var in variables: if str(var)[4:6] == "00": constraint.SetCoefficient(var, 1) ``` ### 流量限制式 4 - 流出終點的總和為 0 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{j}x_{nj} &= 0 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= constraint = solver.Constraint(0, 0) for var in variables: if str(var)[1:3] == f"{len(nodeNumbers)-1:02d}": constraint.SetCoefficient(var, 1) ``` ### 子迴圈限制式 考量實作的方便性，這裡使用 MTZ 限制式消除子迴圈。MTZ 限制式的數學式為： $$ \begin{align} &u_{i} - u_{j} + nx_{ij} \leq n - 1 \\ &u_{i}, u_{j} \geq 0, \space\space\forall i,j = 1,2,\cdots n, \space\space\forall i \neq j \\ &x_{ij} \in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= uVars = [solver.IntVar(0, solver.infinity(), f"u{i:02d}") for i in range(len(nodeNumbers))] for i, ui in enumerate(uVars): for j, uj in enumerate(uVars): for xij in variables: if i != j and str(ui)[1:3] == str(xij)[1:3] and str(uj)[1:3] == str(xij)[4:6]: constraint = solver.Constraint(-solver.infinity(), len(nodeNumbers)-1) constraint.SetCoefficient(ui, 1) constraint.SetCoefficient(uj, -1) constraint.SetCoefficient(xij, len(nodeNumbers)) ``` <br> :::info :bulb: 執行以下程式碼會將模式以數學規劃形式印出 ::: ```python= print(solver.ExportModelAsLpFormat(False).replace('\\', '').replace(',_', ','), sep='\n') ``` ### 模式求解 呼叫 `solver.Solve()` 函式求解，並選取求解後 >0 的決策變數，表示該兩點間存在一條可行路徑。 程式碼實作如下： ```python= solver.Solve() odPairs = [] for xij in variables: if xij.solution_value() > 0: odPairs.append((int(str(xij)[1:3]), int(str(xij)[4:6]))) ``` 接著要將 `odPairs` 轉換成路徑。`odPairs` 的長相大概會長這樣： ```python= # odPairs 長相, 意思是可行路徑分別為 0->4, 1->6, 2->5, ... , 3->7 odPairs = [(0, 4), (1, 6), (2, 5), ... , (3, 7)] ``` 我們要做的是將這些起訖對串成一條路徑 `route`，變成 $0-4-1-6-\cdots-3-7$，下面的程式碼即是做這樣的處理，最後將回傳 `orderList`。 ```python= route = [odPairs[0][0], odPairs[0][1]] odPairs.remove(odPairs[0]) while True: if len(odPairs) < 1: break else: for pair in odPairs: if pair[0] in route: route.append(pair[1]) odPairs.remove(pair) orderList = [nodeNumbers[idx] for idx in route] ``` ### 完整程式碼 ```python= from ortools.linear_solver import pywraplp def getHamtPath(timeMatx, nodeNumbers): """ 使用 SCIP 求解器求解最短漢米爾頓路徑 """ solver = pywraplp.Solver.CreateSolver('SCIP') variables = [solver.IntVar(0, 1, f"x{i:02d}_{j:02d}") for i in range(len(nodeNumbers)) for j in range(len(nodeNumbers))] c = timeMatx for i in range(c.shape[0]): c[i][i] = 1e7 # 目標式 - 最小化總旅行時間 objective = solver.Objective() for var in variables: coeff = int(round(c[int(str(var)[1:3])][int(str(var)[4:6])], 0)) objective.SetCoefficient(var, coeff) objective.SetMinimization() # 流量限制式 1 - 除了終點，其餘節點流出的總和為 1 nodeNumList = [f"{x:02d}" for x in range(len(nodeNumbers)-1)] for node in nodeNumList: constraint = solver.Constraint(1, 1) for var in variables: if str(var)[1:3] == node: constraint.SetCoefficient(var, 1) # 流量限制式 2 - 除了起點, 其餘節點流入的總和為 1 nodeNumList = [f"{x:02d}" for x in range(1, len(nodeNumbers))] for node in nodeNumList: constraint = solver.Constraint(1, 1) for var in variables: if str(var)[4:6] == node: constraint.SetCoefficient(var, 1) # 流量限制式 3 - 流入起點的總和為 0 constraint = solver.Constraint(0, 0) for var in variables: if str(var)[4:6] == "00": constraint.SetCoefficient(var, 1) # 流量限制式 4 - 流出終點的總和為 0 constraint = solver.Constraint(0, 0) for var in variables: if str(var)[1:3] == f"{len(nodeNumbers)-1:02d}": constraint.SetCoefficient(var, 1) # 子迴圈限制式: 使用 MTZ 限制式 uVars = [solver.IntVar(0, solver.infinity(), f"u{i:02d}") for i in range(len(nodeNumbers))] for i, ui in enumerate(uVars): for j, uj in enumerate(uVars): for xij in variables: if i != j and str(ui)[1:3] == str(xij)[1:3] and str(uj)[1:3] == str(xij)[4:6]: constraint = solver.Constraint(-solver.infinity(), len(nodeNumbers)-1) constraint.SetCoefficient(ui, 1) constraint.SetCoefficient(uj, -1) constraint.SetCoefficient(xij, len(nodeNumbers)) # # 執行以下程式碼會將模式以數學規劃形式印出 # print(solver.ExportModelAsLpFormat(False).replace('\\', '').replace(',_', ','), sep='\n') # 模式求解 solver.Solve() odPairs = [] for xij in variables: if xij.solution_value() > 0: odPairs.append((int(str(xij)[1:3]), int(str(xij)[4:6]))) # 路徑轉換 route = [odPairs[0][0], odPairs[0][1]] odPairs.remove(odPairs[0]) while True: if len(odPairs) < 1: break else: for pair in odPairs: if pair[0] in route: route.append(pair[1]) odPairs.remove(pair) orderList = [nodeNumbers[idx] for idx in route] return orderList ``` ## 方法二：使用 CP-SAT 求解 OR-Tools 提供另一種求解整數規劃的求解器，即 CP-SAT。據官方的說法，若所有的變數皆為整數，使用 CP-SAT 求解會更有效率。 ### 求解器和決策變數 首先宣告模型 `model`，並利用 `model.NewIntVar()` 宣告決策變數，函式 args 的填入規則與 MPSolver 求解器中的 `solver.IntVar()` 相同。 ```python= from ortools.sat.python import cp_model # 宣告模型 model = cp_model.CpModel() # 建立二元決策變數，xij = {0,1} variables = [model.NewIntVar(0, 1, f"x{i:02d}_{j:02d}") for i in range(len(nodeNumbers)) for j in range(len(nodeNumbers))] ``` :::info :bulb: 為了方便 debug，這裡將變數 $x_{01}$ 命名為 `x_00_01`，變數 $x_{46}$ 命名為 `x_04_06`。 ::: 接著建立成本矩陣 $c_{ij}$，作法跟上面相同： ```python= c = timeMatx # Add penalty term at cii for i in range(c.shape[0]): c[i][i] = 1e7 ``` ### 目標式 目標式為最小化總成本，這裡的總成本就是總旅行時間，即： $$ \min \space\space \sum_{i}\sum_{j}c_{ij}x_{ij} $$ 這裡會先把所有決策變數的係數先設定好，存到空列表 `objElements`，接著將 `objElements` 丟入 `cp_model.LinearExpr.Sum()` 加總後，再一起包進 `model.Minimize()`，宣告此模式為最小化問題。 :::info :bulb: 目標式係數一定要 `int`，不可以是 `float`, `double`, ...。 :bulb: `var` 的資料型態是 ortools 的 Variable, 訪問它的話要先轉換成 str。 ::: ```python= objElements = [] for var in variables: coeff = int(round(c[int(str(var)[1:3])][int(str(var)[4:6])], 0)) objElements.append(var * coeff) model.Minimize(cp_model.LinearExpr.Sum(objElements)) ``` ### 流量限制式 1 - 除了終點，其餘節點流出的總和為 1 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{i}x_{ij} &= 1 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= nodeNumList = [f"{x:02d}" for x in range(len(nodeNumbers)-1)] for node in nodeNumList: constrntElmts = [] for var in variables: if str(var)[1:3] == node: constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 1) ``` ### 流量限制式 2 - 除了起點, 其餘節點流入的總和為 1 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{i}x_{ij} &= 1 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= nodeNumList = [f"{x:02d}" for x in range(1, len(nodeNumbers))] for node in nodeNumList: constrntElmts = [] for var in variables: if str(var)[4:6] == node: constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 1) ``` ### 流量限制式 3 - 流入起點的總和為 0 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{i}x_{i1} &= 0 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= constrntElmts = [] for var in variables: if str(var)[4:6] == "00": constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 0) ``` ### 流量限制式 4 - 流出終點的總和為 0 假設共有 n 個點 (包括起點和終點)，起點編號為 $1$，終點的編號為 $n$，則將此流量限制式表示成數學式如下： $$ \begin{align} \sum_{j}x_{nj} &= 0 \\ x_{ij} &\in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= constrntElmts = [] for var in variables: if str(var)[1:3] == f"{len(nodeNumbers)-1:02d}": constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 0) ``` ### 子迴圈限制式 這裡同樣使用 MTZ 限制式消除子迴圈： $$ \begin{align} &u_{i} - u_{j} + nx_{ij} \leq n - 1 \\ &u_{i}, u_{j} \geq 0, \space\space\forall i,j = 1,2,\cdots n, \space\space\forall i \neq j \\ &x_{ij} \in {\{0, 1\}}, \space\space\forall i,j = 1,2,\cdots n \end{align} $$ 程式碼實作如下： ```python= uVars = [model.NewIntVar(0, 120, f"u{i:02d}") for i in range(len(nodeNumbers))] for i, ui in enumerate(uVars): for j, uj in enumerate(uVars): for xij in variables: if i != j and str(ui)[1:3] == str(xij)[1:3] and str(uj)[1:3] == str(xij)[4:6]: model.Add(1* ui - 1* uj + int(len(nodeNumbers))* xij <= int(len(nodeNumbers)-1)) ``` <br> :::info :bulb: 執行以下程式碼會將模式以數學規劃形式印出 ::: ```python= print(model.ExportModelAsLpFormat(False).replace('\\', '').replace(',_', ','), sep='\n') ``` ### 模式求解 建立求解器物件 `solver`，再呼叫 `solver.Solve()` 函式求解，並選取求解後 >0 的決策變數，表示該兩點間存在一條可行路徑。 程式碼實作如下： ```python= solver = cp_model.CpSolver() status = solver.Solve(model) odPairs = [] for xij in variables: if solver.Value(xij) > 0: odPairs.append((int(str(xij)[1:3]), int(str(xij)[4:6]))) ``` 接著要將 `odPairs` 轉換成路徑。`odPairs` 的長相大概會長這樣： ```python= # odPairs 長相, 意思是可行路徑分別為 0->4, 1->6, 2->5, ... , 3->7 odPairs = [(0, 4), (1, 6), (2, 5), ... , (3, 7)] ``` 我們要做的是將這些起訖對串成一條路徑 `route`，變成 $0-4-1-6-\cdots-3-7$，下面的程式碼即是做這樣的處理，最後將回傳 `orderList`。 ```python= route = [odPairs[0][0], odPairs[0][1]] odPairs.remove(odPairs[0]) while True: if len(odPairs) < 1: break else: for pair in odPairs: if pair[0] in route: route.append(pair[1]) odPairs.remove(pair) orderList = [nodeNumbers[idx] for idx in route] ``` ### 完整程式碼 ```python= from ortools.sat.python import cp_model def getHamtPath(timeMatx, nodeNumbers): """ 使用 CP-SAT 求解器求解最短漢米爾頓路徑 """ # 宣告模型 model = cp_model.CpModel() # 建立二元決策變數，xij = {0,1} variables = [model.NewIntVar(0, 1, f"x{i:02d}_{j:02d}") for i in range(len(nodeNumbers)) for j in range(len(nodeNumbers))] c = timeMatx for i in range(c.shape[0]): c[i][i] = 1e7 # 目標式 - 最小化總旅行時間 objElements = [] for var in variables: coeff = int(round(c[int(str(var)[1:3])][int(str(var)[4:6])], 0)) objElements.append(var * coeff) model.Minimize(cp_model.LinearExpr.Sum(objElements)) # 流量限制式 1 - 除了終點，其餘節點流出的總和為 1 nodeNumList = [f"{x:02d}" for x in range(len(nodeNumbers)-1)] for node in nodeNumList: constrntElmts = [] for var in variables: if str(var)[1:3] == node: constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 1) # 流量限制式 2 - 除了起點, 其餘節點流入的總和為 1 nodeNumList = [f"{x:02d}" for x in range(1, len(nodeNumbers))] for node in nodeNumList: constrntElmts = [] for var in variables: if str(var)[4:6] == node: constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 1) # 流量限制式 3 - 流入起點的總和為 0 constrntElmts = [] for var in variables: if str(var)[4:6] == "00": constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 0) # 流量限制式 4 - 流出終點的總和為 0 constrntElmts = [] for var in variables: if str(var)[1:3] == f"{len(nodeNumbers)-1:02d}": constrntElmts.append(var * int(1)) model.Add(cp_model.LinearExpr.Sum(constrntElmts) == 0) # 子迴圈限制式 uVars = [model.NewIntVar(0, 120, f"u{i:02d}") for i in range(len(nodeNumbers))] for i, ui in enumerate(uVars): for j, uj in enumerate(uVars): for xij in variables: if i != j and str(ui)[1:3] == str(xij)[1:3] and str(uj)[1:3] == str(xij)[4:6]: model.Add(1* ui - 1* uj + int(len(nodeNumbers))* xij <= int(len(nodeNumbers)-1)) # # 執行以下程式碼會將模式以數學規劃形式印出 # print(model.ExportModelAsLpFormat(False).replace('\\', '').replace(',_', ','), sep='\n') # 模式求解 solver = cp_model.CpSolver() status = solver.Solve(model) odPairs = [] for xij in variables: if solver.Value(xij) > 0: odPairs.append((int(str(xij)[1:3]), int(str(xij)[4:6]))) # 路徑轉換 route = [odPairs[0][0], odPairs[0][1]] odPairs.remove(odPairs[0]) while True: if len(odPairs) < 1: break else: for pair in odPairs: if pair[0] in route: route.append(pair[1]) odPairs.remove(pair) orderList = [nodeNumbers[idx] for idx in route] return orderList ``` ## 方法三：使用 RoutingModel 求解 OR-Tools 的 RoutingModel 可以更有效率地求解路徑規劃問題，依據筆者的經驗，當模式求解的點位 >10，建議用 RoutingModel 求解，求解效率會快上不少。 下面程式碼是根據[官方的範例程式碼](https://colab.research.google.com/github/google/or-tools/blob/stable/examples/notebook/constraint_solver/vrp_starts_ends.ipynb?hl=zh-tw#scrollTo=code)進行改寫。 ### 完整程式碼 ```python= from ortools.constraint_solver import routing_enums_pb2 from ortools.constraint_solver import pywrapcp def getHamtPath(timeMatx, nodeNumbers): """ 使用 RoutingModel 求解最短漢米爾頓路徑 """ data = {} data['distance_matrix'] = timeMatx data['num_vehicles'] = 1 # 這裡是將矩陣的第0列/第0行設為起點，可以根據需求進行更改 data['starts'] = [0] # 這裡是將矩陣的最後1列/最後1行設為終點，同樣可以根據需求進行更改 data['ends'] = [int(data['distance_matrix'].shape[0]-1)] # Create the routing index manager. manager = pywrapcp.RoutingIndexManager( len(data["distance_matrix"]), data["num_vehicles"], data["starts"], data["ends"] ) # Create Routing Model. routing = pywrapcp.RoutingModel(manager) # Create and register a transit callback. 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 = routing.RegisterTransitCallback(distance_callback) # Define cost of each arc. routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index) # Add Distance constraint. dimension_name = "Distance" routing.AddDimension( transit_callback_index, 0, # no slack 86400, # vehicle maximum travel distance True, # start cumul to zero dimension_name, ) distance_dimension = routing.GetDimensionOrDie(dimension_name) distance_dimension.SetGlobalSpanCostCoefficient(100) # Setting first solution heuristic. search_parameters = pywrapcp.DefaultRoutingSearchParameters() search_parameters.first_solution_strategy = ( routing_enums_pb2.FirstSolutionStrategy.SAVINGS ) # Solve the problem. solution = routing.SolveWithParameters(search_parameters) index = routing.Start(0) orderList = [nodeNumbers[index]] while not routing.IsEnd(index): index = solution.Value(routing.NextVar(index)) orderList.append(nodeNumbers[index]) return orderList ``` ## 最佳化演算結果 無論你是用上面三種方法的任意一種，基本上求出來的解都會相同，以下為求解的範例程式碼： ```python= path = getHamtPath(timeMatx, nodeNumbers) print(path) # Output # [0, 1, 6, 4, 3, 5, 2, 7] ``` 將結果繪製於地圖上，這裡使用 `folium` 套件進行地圖繪製： ```python= import folium numCurryShops = len(curryShops) myMap = folium.Map(location=[nodeList[0][1], nodeList[0][0]], zoom_start=13, tiles='CartoDB positron') for node in path: if node > numCurryShops: folium.Marker(location=[nodeList[node][1], nodeList[node][0]], icon=folium.Icon(color='orange', icon='cloud')).add_to(myMap) else: folium.Marker(location=[nodeList[node][1], nodeList[node][0]], icon=folium.Icon(color='blue', icon='cloud')).add_to(myMap) folium.Marker(location=[nodeList[0][1], nodeList[0][0]], icon=folium.Icon(color='pink', icon='cloud')).add_to(myMap) folium.PolyLine([(nodeList[node][1], nodeList[node][0]) for node in path], color='cyan', opacity=0.6).add_to(myMap) folium.TileLayer('cartodbdark_matter').add_to(myMap) display(myMap) ``` ![](https://hackmd.io/_uploads/B1u90s1cT.png) <br> ## 延伸閱讀 - [最短漢米爾頓路徑問題 (I) - 數學模式說明](https://hackmd.io/kdviO126QJq6B4YfxrGzwQ) - [OSRM Python API 使用說明](https://hackmd.io/1kdmHn2xTGS1jjLzC5kVTA)