交通 | Python-MIP库的基础应用——以旅行商问题为例

推文作者：薛桂琴，大连海事大学，研究方向: 即时配送、交通优化

编者按：

本文介绍了python的MIP库的安装以及应用，以经典的旅行商问题为例介绍了MIP库求解的实现方式、对偶变量获取以及如何获得多个可行解。

0 引言

目前，python语言随着各种开源框架的应用而得到越来越多的使用，而运筹学的MIP库在python中扮演着重要的角色。本文介绍了MIP库的安装以及基础的使用，包括约束、变量的添加以及不等式应用等，并以经典的旅行商问题为例，介绍了python的MIP的使用情况。

1 工具包的安装

正常情况下，使用 pip install mip即可完成这个包的安装。

pip install mip

但是在我安装的使用报了下面这样一个错误：

ERROR: Could not install packages due to an EnvironmentError: [WinError 5] 拒绝访问。: ‘d:\…\anaconda3\lib\site-packages\_cffi_backend.cp38-win_amd64.pyd’

Consider using the --useroption or check the permissions.

如果你遇见了这个错误，可以使用在安装指令中添加 --user 选项，最终命令是：

pip install mip --user

2 基础API的使用

2.1 模型类 Model

这是mip库的主类，通过该类可以完成模型的构建、决策变量的添加、约束的添加、目标函数的设置、模型的求解以及解的状态查询。

在构造Model对象时，可以不传递任何参数，也可以指定相关参数：

# 方式1
mp = Model()
# 方式2:这种方式指定了模型的名称、模型的目标函数类型、求解模型要调用的求解器；这里指定为Gurobi，还可以使用开源的 CBC
mp = Model(name='TSP',sense=MINIMIZE,solver_name='GRB')

2.1.1 向模型中加入约束

添加约束涉及的api：

•add_constr(expr, name)：包含约束表达式和约束的名称两个参数

•+= ： 这是mip库特有的添加约束的方式

•如果你的约束中包含求和项，则需调用 xsum(expr) 方法，这个方法类似于gurobipy中的 quicksum() 方法

•add_cut() ：向模型中添加有效不等式时，使用该方法

•add_lazy_constr()：当初始模型不完整且发现整数解时，使用该方法添加规避掉当前整数解的约束

# 添加单个约束
model.add_constr(x1 + x2 <= 9, name='约束1')
model += x1 + x2 <= 9, "约束1"
# 添加求和约束
model += xsum(x[i] for i in range(2)) <= 9, '约束1'

2.1.2 向模型中添加变量

添加决策变量涉及的API：

•add_var(name=‘’, lb=0.0, ub=inf, obj=0.0, var_type=‘C’, column=None)：决策变量的名称、下界、上界、目标系数、类型；Column在列生成时使用

•add_var_tensor(shape, name, **kwargs)：一次性添加多个决策变量，shape属性表示决策变量的维度

x = model.add_var(name='x', lb=0.0, ub=10.0, var_type='C') # 创建一个连续型决策变量
x = model.add_var_tensor(name='x',shape=(3,3),var_type='B') # 创建3*3个二进制决策变量
x = [[ model.add_var(name='x',var_type='I') for i in range(0,3)] for j in range(0,3)] # 创建3*3个整数决策变量

2.1.3 模型状态检查

状态说明：

•ERROR = -1：求解器返回了一个错误

•OPTIMAL = 0：得到的最优解

•INFEASIBLE = 1：建立的模型是一个不可行模型

•UNBOUNDED = 2：目标函数值是无穷的，且存在一个或多个决策变量不受约束条件限制

•FEASIBLE = 3：在搜索过程中找到了整数解，但在证明他是最优解之前，求解被中断

•INT_INFEASIBLE = 4：问题的线性松弛存在可行解，但是不存在整数解

•NO_SOLUTION_FOUND = 5：执行了截断(truncated)搜索，没找到可行解

•LOADED = 6：问题被加载了，但是还有执行优化过程

•CUTOFF = 7：在当前cutoff下没有找到可行的解决方案

代码示例：

status = lp.optimize()
if status == OptimizationStatus.OPTIMAL:print('模型已求解到最优')

2.1.4 几个高效率的方法

下面介绍几个提升编码效率的有用的方法：

•xum()：该方法用于构造带求和形式的线性表达式

•minimize()：该方法用于构造最小化目标函数

•maximize()：该方法用于构造最大化目标函数

•start()：该方法可以为模型设置初始解

•model.threads：该属性可以配置求解模型用到的线程数目，默认使用求解器自身的线程配置；设置的线程多了可能会加速但是会增加内存消耗。

•model.read(filePath)：从指定的路径下读取lp或mps模型文件

2.2 有效不等式的使用

2.2.1 生割或添加lazy_cut

要用的类是：ConstrsGenerator，该类的标准定义为：

class myConstrsGen:def _nint__():# 类的构造函数def generate_constrs(model, depth=0, npass=0):# 编写你的生割代码# 求解器引擎将自动调用该方法，因此这个方法是生割类必须有的方法，不能少

2.2.2 调用coin-or库中的有效不等式

这些有效不等式是定义在 CutType类中的，主要包含以下几种：

•PROBING = 0：Cuts generated evaluating the impact of fixing bounds for integer variables

•GOMORY = 1：Gomory Mixed Integer cuts, as implemented by John Forrest.

•GMI = 2：Gomory Mixed Integer cuts, as implemented by Giacomo Nannicini, focusing on numerically safer cuts.

•RED_SPLIT = 3：Reduce and split cuts [AGY05], implemented by Francois Margot.

•RED_SPLIT_G = 4：Reduce and split cuts [AGY05], implemented by Giacomo Nannicini.

•FLOW_COVER = 5：Lifted Simple Generalized Flow Cover Cut Generator.

•MIR = 6：Mixed-Integer Rounding cuts

•TWO_MIR = 7：Two-phase Mixed-integer rounding cuts.

•LATWO_MIR = 8：Lagrangean relaxation for two-phase Mixed-integer rounding cuts, as in LAGomory

•LIFT_AND_PROJECT = 9：Lift-and-project cuts [BCC93], implemented by Pierre Bonami.

•RESIDUAL_CAPACITY = 10：Residual capacity cuts [AtRa02], implemented by Francisco Barahona.

•ZERO_HALF = 11：Zero/Half cuts

•CLIQUE=12：Clique cuts

•ODD_WHEEL = 13：Lifted odd-hole inequalities.

•KNAPSACK_COVER = 14：Knapsack cover cuts

使用方法是：

 cp = model.generate_cuts([CutType.GOMORY, CutType.MIR, CutType.ZERO_HALF, CutType.KNAPSACK_COVER])

2.2.3 割池 CutPool

CutPool类中仅有一个方法， add()；可以使用该方法将构造的一个或多个cut存储起来，最后统一添加到数学模型中。其典型的用法是：

cp = CutPool()
cp.add(x1 + x2 <= 12)
cp.add(x1 + x2 >= 5)
for cut in cp.cuts:model += cut

2.3 获取多个可行解

解释：

•num_solutions 属性，记录了可行解的个数

•objective_values 属性，记录了可行解的目标函数值

•xi 方法，记录了决策变量的值

for idx in range(multiSol.num_solutions):print(f'第{idx}个可行解的目标函数值为{multiSol.objective_values[idx]}')print(f"x1={x1.xi(idx)},x2={x2.xi(idx)}")

3 应用示例

3.1 旅行商问题

3.1.1 问题描述与模型定义

给定一组城市坐标，旅行商需要从其中的某一个城市出发，每个城市访问一次，最终返回到他出发的城市。

•集合定义：

• $N=0,1,2,...,n-1$：标识城市列表索引
• $A=(i,j),\forall i\in N,j\in N$：表示城市间的连接弧段集合

•参数：

• $c_{ij}$：表示城市$i$到城市$j$之间的距离

•决策变量：

• $x_{ij}$：二进制决策变量，标识城市$i$到城市$j$之间的弧段被旅行商使用时，取值为1；否则取值为0
• $u_i$：整数决策变量，标识城市$i$在旅行商路径中的序号；例如，城市$i$是旅行商需要访问的第5个城市，则对应的决策变量取值为$u_i=4$，注意通常情况下，我们将起点序号标记为0。

目标函数是：最小化旅行商访问所有城市的总距离
$$min:\sum _{i\in N}\sum _{j\in N}{x}_{ij}{c}_{ij}$$

约束条件1：所有节点的出度为1
$$\sum _{j\in N}{x}_{ij}=1,\forall i\in N$$

约束条件2：所有节点的入度为1
$$\sum _{i\in N}{x}_{ij}=1,\forall j\in N$$

约束条件3：MTZ子回路消除约束
$$u_i+1\le u_j+(n-1)\ast (1-x_{ij}),\forall i,j\in N,i\ne 0$$

注：MTZ形式的子回路消除约束中，编码时务必要注意 i 的取值不能为0，否则这个约束会在出发点处，产生矛盾；举个例子，假如我们有3个城市，分别记为0，1，2，得到的一条路径为 0-1-2-0；那么对应的 u 的取值为$u_0=0$
, $u_1=1$,$u_2=2$；由于2还要回到0，这就要求$u_0=3$，就跟前面的$u_0=0$向矛盾了。所以这里要对出发的城市做一个特殊处理。

约束条件4：决策变量的定义域约束
$$u_{i} \in [0,n-1],\forall i \in N \ x_{ij} \in {0,1},\forall i,j \in A$$

3.1.2 实现代码1

代码的编写思路：

•我们在[0, 200 ]之间随机生成节点的坐标，然后以节点间的欧式距离构造节点间的距离矩阵

•利用 mip 库的 Model 类创建tsp的模型

•利用 mip 库的 add_var_tensor() 方法构造决策变量矩阵

•利用 mip 库的 add_constr() 方法向模型中添加约束条件

•利用 mip 库的 xsum() 方法构造具有求和关系的约束表达式

•利用 mip 库的 write() 方法将建立的数学模型导出为 lp 格式的模型文件

•利用 mip 库的 optimize() 方法求解建立的数学模型

•最终我们使用 matplotlib.pyplot 将我们的路径打印出来

from mip import *
import numpy as np
import matplotlib.pyplot as plt
# 定义城市列表、城市坐标
num_of_city = 30
coors = np.random.randint(low=0, high=200, size=(num_of_city,2))
# 定义建模需要的节点集合和弧段集合
N = [i for i in range(0, num_of_city)]
A = [(i,j) for i in N for j in N if i!=j]
# 定义成本矩阵
cost = {(i,j):np.sqrt((coors[i,0]-coors[j,0])**2 + (coors[i,1]-coors[j,1])**2) for i,j in A}
# 构建模型
tsp = Model(name='TSP',sense=MINIMIZE,solver_name='GRB')
# 添加决策变量
x = tsp.add_var_tensor(shape=(num_of_city,num_of_city),var_type='B',name='x')
u = tsp.add_var_tensor(shape=(num_of_city,), name='u')
# 设置目标函数
tsp.objective = xsum(x[i,j]*cost[i,j] for i,j in A)
# 添加节点的入度约束
for j in N:tsp.add_constr(xsum(x[i,j] for i in N if i!=j) == 1)
# 添加节点的出度约束
for i in N:tsp.add_constr(xsum(x[i,j] for j in N if i!=j) == 1)
# 添加MTZ形式的子回路消除约束
for i,j in A:if i!=0:tsp.add_constr(u[i] + 1 <= u[j] + num_of_city*(1-x[i,j]))
tsp.optimize()
tsp.write('tsp.lp')
active_x = [(i,j) for i,j in A if x[i,j].x >= 0.99]
print(active_x)
for i,j in active_x:plt.plot([coors[i,0],coors[j,0]],[coors[i,1],coors[j,1]], c='c')
plt.show()

3.1.3 实现方式2

在上一个实现方式中，我们追求代码简洁，同时倾向于通过方法调用的方式构建模型；在这个实现这个实现方式中，我们中规中矩的使用官方推荐的方式来构建我们的TSP模型。

from mip import *
import numpy as np
import matplotlib.pyplot as plt
# 定义城市列表、城市坐标
num_of_city = 30
coors = np.random.randint(low=0, high=200, size=(num_of_city,2))
# 定义建模需要的节点集合和弧段集合
N = [i for i in range(0, num_of_city)]
A = [(i,j) for i in N for j in N if i!=j]
# 定义成本矩阵
cost = {(i,j):np.sqrt((coors[i,0]-coors[j,0])**2 + (coors[i,1]-coors[j,1])**2) for i,j in A}
# 构建模型
tsp = Model(name='TSP',sense=MINIMIZE,solver_name='GRB')
# 添加决策变量
x = [[tsp.add_var(var_type=BINARY) for i in N] for j in N]
u = [tsp.add_var(var_type='I') for i in N]
print(x)
# 设置目标函数
tsp.objective = xsum(x[i][j]*cost[i,j] for i,j in A if i!=j)
# 添加节点的入度约束
for j in N:tsp += xsum(x[i][j] for i in N if i!=j) == 1, f'indegree_{j}'
# 添加节点的出度约束
for i in N:tsp += xsum(x[i][j] for j in N if i!=j) == 1, f'oudegree_{i}'
# 添加MTZ形式的子回路消除约束
for i,j in A:if i!=0:tsp += u[i] + 1 <= u[j] + num_of_city*(1-x[i][j]),  f'sec_[{i},{j}]'
tsp.optimize()
tsp.write('tsp.lp')
active_x = [(i,j) for i,j in A if x[i][j].x >= 0.99]
print(active_x)
for i,j in active_x:plt.plot([coors[i,0],coors[j,0]],[coors[i,1],coors[j,1]], c='c')
plt.show()

两种实现方式的区别在于：

• 在添加决策变量时，使用 mip 库的 add_var() 方法，以单个变量的形式添加到模型中的；由于列表定义的差异，在使用决策变量时，应该用$x[i][j]$而非第一种方式中的$x[i,j]$

• 在添加约束条件时，使用 mip 库提供的 += 操作符，而不是第一种方式的 add_constr() 方法

3.1.4 实现方式3

由于子回路消除约束是一类特殊的约束，可以在发现子回路时，再将它加入到模型中也是可以的。因此，我们可以首先求解TSP的线性松弛模型，然后根据松弛解构造图，并利用最小割算法生成子回路。最小割算法在networkX中提供有直接可用的方法 minimum_cut。

from mip import *
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx# 定义城市列表、城市坐标
num_of_city = 30
coors = np.random.randint(low=0, high=200, size=(num_of_city,2))
# 定义建模需要的节点集合和弧段集合
N = [i for i in range(0, num_of_city)]
A = [(i,j) for i in N for j in N if i!=j]
# 定义成本矩阵
cost = {(i,j):np.sqrt((coors[i,0]-coors[j,0])**2 + (coors[i,1]-coors[j,1])**2) for i,j in A}
# 构建模型
tsp = Model(name='TSP',sense=MINIMIZE,solver_name='CBC')
# 添加决策变量
x = tsp.add_var_tensor(shape=(num_of_city,num_of_city),var_type='B',name='x')
u = tsp.add_var_tensor(shape=(num_of_city,), name='u')
# 设置目标函数
tsp.objective = xsum(x[i,j]*cost[i,j] for i,j in A)
# 添加节点的入度约束
for j in N:tsp.add_constr(xsum(x[i,j] for i in N if  i!=j) == 1)
# 添加节点的出度约束
for i in N:tsp.add_constr(xsum(x[i,j] for j in N if i!=j) == 1)newConstraints = Truewhile newConstraints:tsp.optimize(relax=True)print("status: {} objective value : {}".format(tsp.status, tsp.objective_value))G = nx.DiGraph()for a in A:G.add_edge(a[0], a[1], capacity=x[a].x)newConstraints = False# 使用最小割方法找到子回路时，将子回路加入到模型当中for (n1, n2) in [(i, j) for (i, j) in A if i != j]:cut_value, (S, NS) = nx.minimum_cut(G, n1, n2)if cut_value <= 0.99:tsp += (xsum(x[a] for a in A if (a[0] in S and a[1] in S)) <= len(S) - 1)newConstraints = True# 如果子回路消除约束不产生cut时，我们使用其他的割平面方法产生有效不等式if not newConstraints and tsp.solver_name.lower() == "cbc":cp = tsp.generate_cuts([CutType.GOMORY, CutType.MIR, CutType.ZERO_HALF, CutType.KNAPSACK_COVER])if cp.cuts:tsp += cpnewConstraints = Truetsp.write('tsp.lp')
active_x = [(i,j) for i,j in A if x[i,j].x >= 0.99]
print(active_x)
for i,j in active_x:plt.plot([coors[i,0],coors[j,0]],[coors[i,1],coors[j,1]], c='c')
plt.show()

3.1.5 实现方式4

最后我们讨论一种lazy cut的实现方式，在这种方式建立在初始模型是不完整的假设上的，也就是当模型捕获到一个整数解时，他会自动检查这个解是否时可行的；如果得到的整数解时不可行的，那么我们会添加一个或几个约束来确保未来这个不可行的整数解不会再出现。

from typing import List, Tuple
from random import seed, randint
from itertools import product
from math import sqrt
import networkx as nx
from mip import Model, xsum, BINARY, minimize, ConstrsGenerator, CutPool
import matplotlib.pyplot as pltclass SubTourCutGenerator(ConstrsGenerator):"""Class to generate cutting planes for the TSP"""def __init__(self, Fl: List[Tuple[int, int]], x_, N_):self.F, self.x, self.N = Fl, x_, N_def generate_constrs(self, model: Model, depth: int = 0, npass: int = 0):xf, N_, cp, G = tsp.translate(self.x), self.N, CutPool(), nx.DiGraph()for (u, v) in [(k, l) for (k, l) in product(N_, N_) if k != l and xf[k][l]]:G.add_edge(u, v, capacity=xf[u][v].x)for (u, v) in F:val, (S, NS) = nx.minimum_cut(G, u, v)if val <= 0.99:aInS = [(xf[i][j], xf[i][j].x) for (i, j) in product(N_, N_) if i != j and xf[i][j] and i in S and j in S]if sum(f for v, f in aInS) >= (len(S)-1)+1e-4:cut = xsum(1.0*v for v, fm in aInS) <= len(S)-1cp.add(cut)if len(cp.cuts) > 256:for cut in cp.cuts:model += cutreturnfor cut in cp.cuts:model += cut
if __name__ == '__main__':num_of_city = 50  # 城市数目N = set(range(num_of_city))seed(0)coors = [(randint(1, 100), randint(1, 100)) for i in N]  # coordinatesArcs = [(i, j) for (i, j) in product(N, N) if i != j]# 距离矩阵c = [[round(sqrt((coors[i][0]-coors[j][0])**2 + (coors[i][1]-coors[j][1])**2)) for j in N] for i in N]tsp = Model()# 二进制决策变量表示弧段是否被旅行商使用vx = [[tsp.add_var(var_type=BINARY) for j in N] for i in N]# 设置最小化距离目标函数tsp.objective = minimize(xsum(c[i][j]*x[i][j] for (i, j) in Arcs))# 出度一次约束for i in N:tsp += xsum(x[i][j] for j in N - {i}) == 1# 入度一次约束for i in N:tsp += xsum(x[j][i] for j in N - {i}) == 1# 移除仅包含两个点的子回路for (i, j) in Arcs:tsp += x[i][j] + x[j][i] <= 1# 对于点集中的每个点，我们找到距离其最远的点，然后首先检查他们是否构成了子回路F, G = [], nx.DiGraph()for (i, j) in Arcs:G.add_edge(i, j, weight=c[i][j])for i in N:P, D = nx.dijkstra_predecessor_and_distance(G, source=i)DS = list(D.items())DS.sort(key=lambda x: x[1])F.append((i, DS[-1][0]))# tsp.cuts_generator = SubTourCutGenerator(F, x, N)tsp.lazy_constrs_generator = SubTourCutGenerator(F, x, N)tsp.optimize()print('解的状态: %s 路径长度 %g' % (tsp.status, tsp.objective_value))tsp.write('tsp.lp')active_x = [(i,j) for i,j in Arcs if x[i][j].x >= 0.99]print(active_x)for i,j in active_x:plt.plot([coors[i][0],coors[j][0]],[coors[i][1],coors[j][1]], c='c')for i,txt in enumerate(N):plt.annotate(txt,(coors[i][0]+0.5, coors[i][1]+0.5))plt.show()

3.2获取线性规划问题的对偶变量

3.2.1 一个简单的线性规划问题

下面给出一个简单的线性规划问题，其对应的表达式为：

$$max:2x_1+3x_2-5x_3$$
$$s.t$$
$$\x_{1} + x_{2} + x_{3} = 7 \ 2x_{1} - 5x_{2} + x_{3} \geq 10 \ x_{1} + 3x_{2} + x_{3} \leq 12 \ x_{1},x_{2},x_{3} \geq 0$$

求解盖模型的代码为：

from mip import *
lp = Model()
x = lp.add_var_tensor(shape=(3,),var_type='C',name='x')
lp.objective = maximize(2*x[0] + 3*x[1] - 5*x[2])
lp += x[0] + x[1] + x[2] == 7, '约束1'
lp += 2*x[0] - 5*x[1] + x[2] >= 10, '约束2'
lp += x[0] + 3*x[1] + x[2] <= 12, '约束3'
lp.optimize()
active_x = [x[i].x for i in range(3)]
print(active_x)

模型的最优解为：[6.428571428571429, 0.5714285714285712, 0.0]，最优值为14.571

3.2.2 获取感兴趣的模型信息

• 获取模型对应的对偶变量

active_pi = [lp.constr_by_name(f'约束{i}').pi for i in range(1,4)]
print(active_pi)

模型的对偶变量取值为：[2.2857142857142856, -0.14285714285714285, 0.0]

• 获取模型的目标函数值

print(f'目标函数值为{lp.objective_value}')
print(f'目标函数值的界为{lp.objective_bound}')

3.3 获取优化模型的多个可行解​

​

3.3.1 简单的整数规划示例​

​下面给出一个简单的整数规划问题，其对应的表达式为：

$$min:5x1+5x2$$
$$s.t$$
$$x1+x2\ge 6$$
$$x1+x2\ge 4$$
$$x1+3x2\ge 6$$
$$x1,x2\in Z^{+}$$

求解盖模型的代码为：

3.3.2 获取多个可行解的验证代码

from mip import *
multiSol = Model()
x1 = multiSol.add_var(var_type='I',name='x1')
x2 = multiSol.add_var(var_type='I',name='x1')
multiSol.objective = 5*x1 + 5*x2
multiSol += 3*x1 + x2 >= 6
multiSol += x1 + x2 >= 4
multiSol += x1 + 3*x2 >= 6
status = multiSol.optimize(max_solutions=10)
if status == OptimizationStatus.OPTIMAL:print('模型已求解到最优')for idx in range(multiSol.num_solutions):print(f'第{idx}个可行解的目标函数值为{multiSol.objective_values[idx]}')print(f"x1={x1.xi(idx)},x2={x2.xi(idx)}")

打印结果：

模型已求解到最优
第0个可行解的目标函数值为20.0
x1=3.0,x2=1.0
第1个可行解的目标函数值为30.0
x1=0.0,x2=6.0

4.总结与展望

本文以旅行商问题为例，介绍了python的MIP库的应用，实现了在求解问题之外的多种应用，这对未来算法优化打下了坚实的基础。

参考文献：

https://buildmedia.readthedocs.org/media/pdf/python-mip/latest/python-mip.pdf