基于Gurobi软件Callback功能的旅行商问题求解

度巍 陈昊泽

摘要：作为经典组合优化问题，旅行商问题（Traveling Salesman Problem简称TSP） 一直是大学交通运输与应用数学等专业的教学与科研热点。在基于混合整数规划模型的TSP求解中，需要解决如何避免出现子环路问题，Gurobi作为当前最先进的运筹优化软件，其具有的Callback功能使模型在求解过程中，动态地添加子环路约束成为可能。文章针对当前相关网络资源存在的问题，构建了用Python编写的基于Callback功能动态添加子环路消除约束的TSP求解代码，通过多个算例验证了代码的求解可行性，为逐步将Gurobi引入课堂教学提供了素材。

关键词：旅行商问题;子环路消除;Gurobi;Callback功能

中图分类号：G642        文献标识码：A

文章编号：1009-3044（2022）25-0009-02

开放科学（资源服务） 标识码（OSID） ：

1 引言

新一代大规模优化软件Gurobi，因其在求解数学规划问题的卓越性能，逐渐成为高校师生运筹学教学与科研首选软件。在Gurobi的高级功能中，Callback函数可获取求解过程中的信息，动态加入新的约束条件或者其他算法，为实现各种复杂问题的求解创造了条件。本文通过Callback功能，在求解旅行商问题过程中，动态添加子环路消除约束条件，实现了一种新的求解旅行商问题方法，并通过与现有模型在不同问题规模求解时间的对比，帮助学生加深对旅行商问题的理解，为开展Gurobi软件的教学与相关科研提供了素材。

2 旅行商问题的数学规划模型

组合优化领域中的旅行商问题（Traveling Salesman Problem簡称TSP） ，广泛应用物流配送、电缆和光缆布线等领域[1]，如何找到大规模TSP的最优解一直是国内外学者的研究热点[2]，近几十年来涌现出各种求解算法。旅行商问题一般表述为：某旅行推销员要去若干个城市推销商品，然后回到其出发城市，已知任意两个城市之间的行走距离，求出该推销员经过每个城市有且一次同时总距离最短的巡回路线。若城市个数为[n]，城市集合为[V={1，2，…，n}]，从城市[i]到城市[j]的距离为[cij]，若在巡回路线中，从城市[i]直接走到城市[j]，0-1变量[xij]取值为1，否则取值为0，则求解旅行商问题涉及如下整数规划模型：[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）中的目标函数表示巡回路线总的长度，第1类与第2类约束条件分别表示每个城市节点只能进出一次。然而如果仅仅是求解模型（1），获得的解有可能出现子环路情形，即解所对应的路线中存在少于[n]的若干城市节点首尾相连的情形。为了避免该子环路问题，有学者做过大量的研究，如在模型（1）上增加如下MTZ约束：

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

（2）式中，实数决策变量[ui]表示城市[i]在巡回路线中的序号。添加MTZ约束是目前大多数运筹优化软件在避免出现子环路时采取的方法。如高校运筹学教学应用广泛的LINGO软件，在其模型库中，求解TSP采取的是文献[4]对应的增强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）

解决避免出现子环路的另一个思路是在模型（1）上添加如下子环路消除约束：

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

其中[S]表示所有的子环路集合，[S]表示[S]包含的城市节点个数，（4）式的理解非常直观，即让[S]中依次相连的节点之间对应的[xij]之和少于[S]，从而把子环路出现时对应的解排除掉。然而，由于该约束数量与[2n]同阶，过去很难应用于实际问题的求解中。目前Gurobi软件具有的Callback功能[5]，可以动态地将求解过程中出现的每个子环路所对应的消除约束（4）式依次添加到模型中，基于Callback功能的完整求解流程如图1所示：

求解TSP流程图

3 模型求解的关键代码

基于上述途径求解TSP的代码用Python语言实现，通过导入gurobipy库实现调用Gurobi的Callback功能，部分代码参考了网络资源[6]。笔者在对文献[6]的考察中，发现其相关代码并不能求出最优解，进一步对代码做了改进完善。由于代码较长，只给出关键代码部分。

检查当前解所对应的巡游路线函数段：

def subtour（graph）：

unvisited = range（0，n）  #构造探寻可能构成子环路城市节点序列

cycle = range（0， n）  # cycle用来返回找到的子环路

edges = findEdges（graph）

edges = tuplelist（edges）

tempcycle=[]

thiscycle1 = [unvisited[0]] #从unvisited的第一个节点可以寻找子环路

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

对当前解是否存在子环路进行检查，如果存在子环路，构造出对应的约束条件（4）函数段：

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函数获取当前解对应的环路

if （len（tour） < 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）

不同于文献[6]在构建基本模型时额外添加一个虚拟起始点的做法，本文代码直接构建模型（1）中的决策变量与目标函数：

model = Model（'TSP'）

X = {}  #構建决策变量

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）  # 构建目标函数

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 实例求解分析

结合上面三个关键段以及导入数据，调用Gurobi求解模型等代码部分，能实现对TSP的求解。在Windows10操作系统、8G内存、Inter Core2.5GHz实验平台上，做了一系列不同规模TSP的求解实验，对较小规模的TSP，本文构造的模型可以较快地求出最优巡回路径，如针对文献[6]提到的Solomon算例集中的C201节点数据，取前11个节点，能在几秒内计算出最优巡游序列为1→7→9→10→11→6→3→2→8→4→5→1，进一步将本文的模型与基于MTZ约束的模型在不同节点数量下的相同问题求解时间进行了对比，结果如表1所示：

可以看到，在实例规模较小时，两种模型求解时间差别不大，但随着问题规模的增加，本文构建的Callback动态添加子环路消除约束模型在求解时间上迅速增加。这是因为子环路个数随着问题规模的增加呈现指数增加，模型在每一次添加新的子环路消除约束后又要重新计算新的解，而相比较下，基于MTZ约束的模型一次性添加好所有约束，在求解时间上明显更有优势，这很好解释了当前类似LINGO、GAMS等软件在求解TSP时基于MTZ约束的原因。

5 结束语

旅行商問题是高校交通运输本科专业开设的诸如《交通运筹学》《交通系统分析》等课程中重要知识点，目前课堂的教学往往是对相关模型与求解算法进行介绍，这很难使学生对该问题形成深刻印象，更不能体会到求解旅行商问题的难度与应用价值。本文基于Gurobi软件的求解优化问题强大功能，实现了过去无法做到的通过添加子环路消除约束求解旅行商问题途径，让学生直观感受到不同的求解模型在计算时间上的差别，加深了对该问题的理解，直观体会到问题的计算复杂程度。同时也将本科生在计算机课程学到的Python语言应用于专业知识实践，实现了大学各个课程间的相互融入，为学生进一步学习车辆路径问题、交通分配模型打下基础。

参考文献：

[1] William J.Cook.迷茫的旅行商一个无处不在的计算机算法问题[M]. 隋春宁，译.北京：人民邮电出版社，2013.

[2] 李汝佳.基于蚁群算法的旅游路线规划问题研究[J].电脑知识与技术，2019，15（8）：137-140.

[3] 韩中庚.运筹学及其工程应用[M].北京：清华大学出版社，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中两种不同消除子环路的方法及callback实现[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.

【通联编辑：王力】