基于Gurobi軟件Callback功能的旅行商問(wèn)題求解

度巍 陳昊澤

摘要：作為經(jīng)典組合優(yōu)化問(wèn)題，旅行商問(wèn)題（Traveling Salesman Problem簡(jiǎn)稱(chēng)TSP） 一直是大學(xué)交通運(yùn)輸與應(yīng)用數(shù)學(xué)等專(zhuān)業(yè)的教學(xué)與科研熱點(diǎn)。在基于混合整數(shù)規(guī)劃模型的TSP求解中，需要解決如何避免出現(xiàn)子環(huán)路問(wèn)題，Gurobi作為當(dāng)前最先進(jìn)的運(yùn)籌優(yōu)化軟件，其具有的Callback功能使模型在求解過(guò)程中，動(dòng)態(tài)地添加子環(huán)路約束成為可能。文章針對(duì)當(dāng)前相關(guān)網(wǎng)絡(luò)資源存在的問(wèn)題，構(gòu)建了用Python編寫(xiě)的基于Callback功能動(dòng)態(tài)添加子環(huán)路消除約束的TSP求解代碼，通過(guò)多個(gè)算例驗(yàn)證了代碼的求解可行性，為逐步將Gurobi引入課堂教學(xué)提供了素材。

關(guān)鍵詞：旅行商問(wèn)題;子環(huán)路消除;Gurobi;Callback功能

中圖分類(lèi)號(hào)：G642? ? ? ? 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2022）25-0009-02

開(kāi)放科學(xué)（資源服務(wù)） 標(biāo)識(shí)碼（OSID） ：

1 引言

新一代大規(guī)模優(yōu)化軟件Gurobi，因其在求解數(shù)學(xué)規(guī)劃問(wèn)題的卓越性能，逐漸成為高校師生運(yùn)籌學(xué)教學(xué)與科研首選軟件。在Gurobi的高級(jí)功能中，Callback函數(shù)可獲取求解過(guò)程中的信息，動(dòng)態(tài)加入新的約束條件或者其他算法，為實(shí)現(xiàn)各種復(fù)雜問(wèn)題的求解創(chuàng)造了條件。本文通過(guò)Callback功能，在求解旅行商問(wèn)題過(guò)程中，動(dòng)態(tài)添加子環(huán)路消除約束條件，實(shí)現(xiàn)了一種新的求解旅行商問(wèn)題方法，并通過(guò)與現(xiàn)有模型在不同問(wèn)題規(guī)模求解時(shí)間的對(duì)比，幫助學(xué)生加深對(duì)旅行商問(wèn)題的理解，為開(kāi)展Gurobi軟件的教學(xué)與相關(guān)科研提供了素材。

2 旅行商問(wèn)題的數(shù)學(xué)規(guī)劃模型

組合優(yōu)化領(lǐng)域中的旅行商問(wèn)題（Traveling Salesman Problem簡(jiǎn)稱(chēng)TSP） ，廣泛應(yīng)用物流配送、電纜和光纜布線等領(lǐng)域[1]，如何找到大規(guī)模TSP的最優(yōu)解一直是國(guó)內(nèi)外學(xué)者的研究熱點(diǎn)[2]，近幾十年來(lái)涌現(xiàn)出各種求解算法。旅行商問(wèn)題一般表述為：某旅行推銷(xiāo)員要去若干個(gè)城市推銷(xiāo)商品，然后回到其出發(fā)城市，已知任意兩個(gè)城市之間的行走距離，求出該推銷(xiāo)員經(jīng)過(guò)每個(gè)城市有且一次同時(shí)總距離最短的巡回路線。若城市個(gè)數(shù)為[n]，城市集合為[V={1，2，…，n}]，從城市[i]到城市[j]的距離為[cij]，若在巡回路線中，從城市[i]直接走到城市[j]，0-1變量[xij]取值為1，否則取值為0，則求解旅行商問(wèn)題涉及如下整數(shù)規(guī)劃模型：[3]

[min Z=i=1nj=1ncijxij? ?i，j∈V，i≠js.t.? ? ? ? ? ?i=1nxij=1，? ??j∈V，i≠j? ? ? ? ? ? ? ? j=1nxij=1，? ??i∈V，i≠j? ? ? ? ? ? ? ? xij∈{0，1}， ?i，j∈V，i≠j]? ? ? ? （1）

模型（1）中的目標(biāo)函數(shù)表示巡回路線總的長(zhǎng)度，第1類(lèi)與第2類(lèi)約束條件分別表示每個(gè)城市節(jié)點(diǎn)只能進(jìn)出一次。然而如果僅僅是求解模型（1），獲得的解有可能出現(xiàn)子環(huán)路情形，即解所對(duì)應(yīng)的路線中存在少于[n]的若干城市節(jié)點(diǎn)首尾相連的情形。為了避免該子環(huán)路問(wèn)題，有學(xué)者做過(guò)大量的研究，如在模型（1）上增加如下MTZ約束：

[ui-uj+nxij≤n-1， ?i，j∈V，i≠j]? ? ? ? （2）

（2）式中，實(shí)數(shù)決策變量[ui]表示城市[i]在巡回路線中的序號(hào)。添加MTZ約束是目前大多數(shù)運(yùn)籌優(yōu)化軟件在避免出現(xiàn)子環(huán)路時(shí)采取的方法。如高校運(yùn)籌學(xué)教學(xué)應(yīng)用廣泛的LINGO軟件，在其模型庫(kù)中，求解TSP采取的是文獻(xiàn)[4]對(duì)應(yīng)的增強(qiáng)MTZ條件：

[uj≥ui+xij-（n-2）（1-xij）+（n-3）xji? ?i，j∈{2…n}，i≠jui≤n-1-（n-2）x1i? ? ? ? ? ? ? ? ? ? ? ? ? ? ??i∈{2…n}ui≥1+（n-2）x1i? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??i∈{2…n}]? ? ? （3）

解決避免出現(xiàn)子環(huán)路的另一個(gè)思路是在模型（1）上添加如下子環(huán)路消除約束：

[i，j∈Sxij≤S-1， 2≤S≤n-1，S?V]? ? ? （4）

其中[S]表示所有的子環(huán)路集合，[S]表示[S]包含的城市節(jié)點(diǎn)個(gè)數(shù)，（4）式的理解非常直觀，即讓[S]中依次相連的節(jié)點(diǎn)之間對(duì)應(yīng)的[xij]之和少于[S]，從而把子環(huán)路出現(xiàn)時(shí)對(duì)應(yīng)的解排除掉。然而，由于該約束數(shù)量與[2n]同階，過(guò)去很難應(yīng)用于實(shí)際問(wèn)題的求解中。目前Gurobi軟件具有的Callback功能[5]，可以動(dòng)態(tài)地將求解過(guò)程中出現(xiàn)的每個(gè)子環(huán)路所對(duì)應(yīng)的消除約束（4）式依次添加到模型中，基于Callback功能的完整求解流程如圖1所示：

求解TSP流程圖

3 模型求解的關(guān)鍵代碼

基于上述途徑求解TSP的代碼用Python語(yǔ)言實(shí)現(xiàn)，通過(guò)導(dǎo)入gurobipy庫(kù)實(shí)現(xiàn)調(diào)用Gurobi的Callback功能，部分代碼參考了網(wǎng)絡(luò)資源[6]。筆者在對(duì)文獻(xiàn)[6]的考察中，發(fā)現(xiàn)其相關(guān)代碼并不能求出最優(yōu)解，進(jìn)一步對(duì)代碼做了改進(jìn)完善。由于代碼較長(zhǎng)，只給出關(guān)鍵代碼部分。

檢查當(dāng)前解所對(duì)應(yīng)的巡游路線函數(shù)段：

def subtour（graph）：

unvisited = range（0，n）? #構(gòu)造探尋可能構(gòu)成子環(huán)路城市節(jié)點(diǎn)序列

cycle = range（0， n）? # cycle用來(lái)返回找到的子環(huán)路

edges = findEdges（graph）

edges = tuplelist（edges）

tempcycle=[]

thiscycle1 = [unvisited[0]] #從unvisited的第一個(gè)節(jié)點(diǎn)可以尋找子環(huán)路

unvisited.remove（0）

while len（thiscycle1）

point1=thiscycle1[-1]

neighbors1 = [j for i， j in edges.select（point1， '*'） if j in unvisited]

if? len（neighbors1） >= 1：

thiscycle1.extend（neighbors1）

if （thiscycle1[-1]， thiscycle1[0]） in edges：

tempcycle=thiscycle1

break

cycle = tempcycle

return cycle

對(duì)當(dāng)前解是否存在子環(huán)路進(jìn)行檢查，如果存在子環(huán)路，構(gòu)造出對(duì)應(yīng)的約束條件（4）函數(shù)段：

def subtourelim（model， where）：

if （where == GRB.Callback.MIPSOL）：

x_value = np.zeros（[nodeNum ， nodeNum ]）

for m in model.getVars（）：

if （m.varName.startswith（'x'））：

a = （int）（m.varName.split（'_'）[1]）

b = （int）（m.varName.split（'_'）[2]）

x_value[a][b] = model.cbGetSolution（m）

tour = subtour（x_value）? #從subtour函數(shù)獲取當(dāng)前解對(duì)應(yīng)的環(huán)路

if （len（tour） < n）：? #如果得到的環(huán)路包含節(jié)點(diǎn)個(gè)數(shù)小于城市總數(shù)，即為子環(huán)路

# 將該子環(huán)路所對(duì)應(yīng)的子環(huán)路消除約束添加到模型中

tt1 = LinExpr（0）

for i in range（0，len（tour））：

if （i < len（tour）-1）：

tt1.addTerms（1， X[tour[i]，tour[i+1]]）

elif （i == len（tour）-1）：

tt1.addTerms（1， X[tour[i]，tour[0]]）

model.cbLazy（tt1<= len（tour） - 1）

不同于文獻(xiàn)[6]在構(gòu)建基本模型時(shí)額外添加一個(gè)虛擬起始點(diǎn)的做法，本文代碼直接構(gòu)建模型（1）中的決策變量與目標(biāo)函數(shù)：

model = Model（'TSP'）

X = {}? #構(gòu)建決策變量

mu = {}

for j in range（nodeNum）：

if （i ！= j）：

X[i， j] = model.addVar（vtype=GRB.BINARY， name='x_' + str（i） + '_' + str（j））

obj = LinExpr（0）? # 構(gòu)建目標(biāo)函數(shù)

for key in X.keys（）：

i = key[0]

j = key[1]

if （i < nodeNum and j < nodeNum）： obj.addTerms（cost[key[0]][key[1]]， X[key]）

4 實(shí)例求解分析

結(jié)合上面三個(gè)關(guān)鍵段以及導(dǎo)入數(shù)據(jù)，調(diào)用Gurobi求解模型等代碼部分，能實(shí)現(xiàn)對(duì)TSP的求解。在Windows10操作系統(tǒng)、8G內(nèi)存、Inter Core2.5GHz實(shí)驗(yàn)平臺(tái)上，做了一系列不同規(guī)模TSP的求解實(shí)驗(yàn)，對(duì)較小規(guī)模的TSP，本文構(gòu)造的模型可以較快地求出最優(yōu)巡回路徑，如針對(duì)文獻(xiàn)[6]提到的Solomon算例集中的C201節(jié)點(diǎn)數(shù)據(jù)，取前11個(gè)節(jié)點(diǎn)，能在幾秒內(nèi)計(jì)算出最優(yōu)巡游序列為1→7→9→10→11→6→3→2→8→4→5→1，進(jìn)一步將本文的模型與基于MTZ約束的模型在不同節(jié)點(diǎn)數(shù)量下的相同問(wèn)題求解時(shí)間進(jìn)行了對(duì)比，結(jié)果如表1所示：

可以看到，在實(shí)例規(guī)模較小時(shí)，兩種模型求解時(shí)間差別不大，但隨著問(wèn)題規(guī)模的增加，本文構(gòu)建的Callback動(dòng)態(tài)添加子環(huán)路消除約束模型在求解時(shí)間上迅速增加。這是因?yàn)樽迎h(huán)路個(gè)數(shù)隨著問(wèn)題規(guī)模的增加呈現(xiàn)指數(shù)增加，模型在每一次添加新的子環(huán)路消除約束后又要重新計(jì)算新的解，而相比較下，基于MTZ約束的模型一次性添加好所有約束，在求解時(shí)間上明顯更有優(yōu)勢(shì)，這很好解釋了當(dāng)前類(lèi)似LINGO、GAMS等軟件在求解TSP時(shí)基于MTZ約束的原因。

5 結(jié)束語(yǔ)

旅行商問(wèn)題是高校交通運(yùn)輸本科專(zhuān)業(yè)開(kāi)設(shè)的諸如《交通運(yùn)籌學(xué)》《交通系統(tǒng)分析》等課程中重要知識(shí)點(diǎn)，目前課堂的教學(xué)往往是對(duì)相關(guān)模型與求解算法進(jìn)行介紹，這很難使學(xué)生對(duì)該問(wèn)題形成深刻印象，更不能體會(huì)到求解旅行商問(wèn)題的難度與應(yīng)用價(jià)值。本文基于Gurobi軟件的求解優(yōu)化問(wèn)題強(qiáng)大功能，實(shí)現(xiàn)了過(guò)去無(wú)法做到的通過(guò)添加子環(huán)路消除約束求解旅行商問(wèn)題途徑，讓學(xué)生直觀感受到不同的求解模型在計(jì)算時(shí)間上的差別，加深了對(duì)該問(wèn)題的理解，直觀體會(huì)到問(wèn)題的計(jì)算復(fù)雜程度。同時(shí)也將本科生在計(jì)算機(jī)課程學(xué)到的Python語(yǔ)言應(yīng)用于專(zhuān)業(yè)知識(shí)實(shí)踐，實(shí)現(xiàn)了大學(xué)各個(gè)課程間的相互融入，為學(xué)生進(jìn)一步學(xué)習(xí)車(chē)輛路徑問(wèn)題、交通分配模型打下基礎(chǔ)。

參考文獻(xiàn)：

[1] William J.Cook.迷茫的旅行商一個(gè)無(wú)處不在的計(jì)算機(jī)算法問(wèn)題[M]. 隋春寧，譯.北京：人民郵電出版社，2013.

[2] 李汝佳.基于蟻群算法的旅游路線規(guī)劃問(wèn)題研究[J].電腦知識(shí)與技術(shù)，2019，15（8）：137-140.

[3] 韓中庚.運(yùn)籌學(xué)及其工程應(yīng)用[M].北京：清華大學(xué)出版社，2014.

[4] Desrochers M，Laporte G.Improvements and extensions to the Miller-Tucker-Zemlin subtour elimination constraints[J].Operations Research Letters，1991，10（1）：27-36.

[5] Gurobi OptimizationCompany. Gurobi documents[EB/OL].[2021-10-01] https：//www.gurobi.com/documentation/9.5/refman/index.html

[6] 劉興祿.TSP中兩種不同消除子環(huán)路的方法及callback實(shí)現(xiàn)[EB/OL].[2021-01-16]. https：//blog.csdn.net/weixin_398332 90/article/details/113452735？utm_medium=distribute.pc_relevant.none-task-blog-2～default～baidujs_title～default-0.no_ search_link&spm=1001.2101.3001.4242.1.

【通聯(lián)編輯：王力】