基于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）<n：

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约束的原因。