1. 介绍
CBC和Gurobi的一个python封装，支持cut generation, lazy constraints, MIP starts以及solution pools。
自带CBC，省去安装的麻烦。
python-mip使用cffi库调用C代码，使用Pypy进行编译，据称性能很好，比gurobi自带的python接口快25倍，比JuMP还要快（但总体上差的不多）。
2. 建模与求解
2.1 基本步骤
定义模型
m = Model() # default is CBC
m = Model(sense=MAXIMIZE, solver_name=GRB) # use GRB for Gurobi
添加目标函数：
m.objective = xsum(c[i]*x[i] for i in range(n)) # 默认是min
m.objective = maximize(xsum(c[i]*x[i] for i in range(n))) # 手动添加
添加变量：
x = m.add_var() # 默认是continues，此外还有binary和integer。此外还可以定义lb和ub，以及给变量命名
z = m.add_var(name='zCost', var_type=INTEGER, lb=-10, ub=10) # 默认是非负的，除非人工设置lb
y = [ m.add_var(var_type=BINARY) for i in range(10) ] # 变量数组
vz = m.var_by_name('zCost') # 获取变量的方法
vz.ub = 5
添加约束，非常独特，没有方法，直接用+
m += x + y <= 10
m.add_constr( x + y <= 1 ) # 也可以这么写
m += xsum(w[i]*x[i] for i in range(n) if i%2 == 0) <= c, 'even_sum' # 这里同时给约束进行了命名
constraint = m.constr_by_name('even_sum') # 获取constraint的方法
保存和读写，支持lp和mps
m.write('model.lp')
m.read('model.lp')
print('model has {} vars, {} constraints and {} nzs'.format(m.num_cols, m.num_rows, m.num_nz))
求解，调用optimize。除了下面的max_gap和max_seconds，还有max_nodes, max_solutions等
m.max_gap = 0.05
status = m.optimize(max_seconds=300)
if status == OptimizationStatus.OPTIMAL:
print('optimal solution cost {} found'.format(m.objective_value))
elif status == OptimizationStatus.FEASIBLE:
print('sol.cost {} found, best possible: {}'.format(m.objective_value, m.objective_bound))
elif status == OptimizationStatus.NO_SOLUTION_FOUND:
print('no feasible solution found, lower bound is: {}'.format(m.objective_bound))
if status == OptimizationStatus.OPTIMAL or status == OptimizationStatus.FEASIBLE:
print('solution:')
for v in m.vars:
if abs(v.x) > 1e-6: # only printing non-zeros
print('{} : {}'.format(v.name, v.x))
返回类型：
OPTIMAL
FEASIBLE
NO_SOLUTION_FOUND
INFEASIBLE or INT_INFEASIBLE
UNBOUNDED
ERROR
求解优先级：
default setting:
tries to balance between the search of improved feasible solutions and improved lower bounds;
feasibility:
focus on finding improved feasible solutions in the first moments of the search process, activates heuristics;
optimality:
activates procedures that produce improved lower bounds, focusing in pruning the search tree even if the production of the first feasible solutions is delayed.
2.2 不常见的步骤
- add_var时有个column选项，用于列生成算法
- x = m.add_var_tensor((3, 5), “x”)，生成矩阵变量
- optimize可以加上preprocess做预处理，也可以用feasibility pump启发式算法来找初始解。
- 在optimize时可以设置max_seconds_same_incumbent=inf, max_nodes_same_incumbent=1073741824，如果有找到可行解并且搜索没有改进，那么停止搜索。
- optimize可以指定threads，-1的话会用满所有的核。
- read除了lp和mps文件，还可以读入sol（初始整数可行解），bas（线性松弛问题的最优基）
- generate_constrs(model, depth=0, npass=0)用于割平面和lasy constraints
- IncumbentUpdater(model)可以用于接收是否有新的整数可行解生成
- update_incumbent(objective_value, best_bound, solution)可用于更新bound，其中solution类型为(List[Tuple[mip.Var,float]])
10.使用compute_features(model)和features可以快速得到模型画像，可以用于机器学习算法
2.3 常见cut清单
cut使用的是COIN-OR Cut Generation Library，可用的有：
CLIQUE = 12
Clique cuts [Padb73].
FLOW_COVER = 5
Lifted Simple Generalized Flow Cover Cut Generator.
GMI = 2
Gomory Mixed Integer cuts [Gomo69], as implemented by Giacomo Nannicini, focusing on numerically safer cuts.
GOMORY = 1
Gomory Mixed Integer cuts [Gomo69], as implemented by John Forrest.
KNAPSACK_COVER = 14
Knapsack cover cuts [Bala75].
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.
MIR = 6
Mixed-Integer Rounding cuts [Marc01].
ODD_WHEEL = 13
Lifted odd-hole inequalities.
PROBING = 0
Cuts generated evaluating the impact of fixing bounds for integer variables
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.
RESIDUAL_CAPACITY = 10
Residual capacity cuts [AtRa02], implemented by Francisco Barahona.
TWO_MIR = 7
Two-phase Mixed-integer rounding cuts.
ZERO_HALF = 11
Zero/Half cuts [Capr96].
2.4 求解状态清单
CUTOFF = 7
No feasible solution exists for the current cutoff
ERROR = -1
Solver returned an error
FEASIBLE = 3
An integer feasible solution was found during the search but the search was interrupted before concluding if this is the optimal solution or not.
INFEASIBLE = 1
The model is proven infeasible
INT_INFEASIBLE = 4
A feasible solution exist for the relaxed linear program but not for the problem with existing integer variables
LOADED = 6
The problem was loaded but no optimization was performed
NO_SOLUTION_FOUND = 5
A truncated search was executed and no integer feasible solution was found
OPTIMAL = 0
Optimal solution was computed
UNBOUNDED = 2
One or more variables that appear in the objective function are not included in binding constraints and the optimal objective value is infinity.
2.5 求解器清单
线性规划求解：
AUTO = 0
Let the solver decide which is the best method
BARRIER = 3
The barrier algorithm
DUAL = 1
The dual simplex algorithm
PRIMAL = 2
The primal simplex algorithm
2.6 SOS
SOS1：集合中只允许有一个非零变量
SOS2：集合中允许连续两个非零变量，用于piecewise linear approximations of non-linear functions，如下，使用add_sos()方法
下面是一个拟合log目标函数的例子：
3. 经典问题建模
3.1 knapsack：基本
p：profit；w：weight
p = [10, 13, 18, 31, 7, 15]
w = [11, 15, 20, 35, 10, 33]
c, I = 47, range(len(w))
m = Model("knapsack")
x = [m.add_var(var_type=BINARY) for i in I]
m.objective = maximize(xsum(p[i] * x[i] for i in I))
m += xsum(w[i] * x[i] for i in I) <= c
m.optimize()
selected = [i for i in I if x[i].x >= 0.99]
print("selected items: {}".format(selected))
3.2 TSP：引入cutting plane等技巧
c：cost，
第三行约束用来消除sub-tour，如图。增加中间变量y用于标识每个点是路径上的第几个点；当然也有其他的消除subtour的方式，即任意两个子集之间都需要由连接。
from itertools import product
from sys import stdout as out
from mip import Model, xsum, minimize, BINARY
# names of places to visit
places = ['Antwerp', 'Bruges', 'C-Mine', 'Dinant', 'Ghent',
'Grand-Place de Bruxelles', 'Hasselt', 'Leuven',
'Mechelen', 'Mons', 'Montagne de Bueren', 'Namur',
'Remouchamps', 'Waterloo']
# distances in an upper triangular matrix
dists = [[83, 81, 113, 52, 42, 73, 44, 23, 91, 105, 90, 124, 57],
[161, 160, 39, 89, 151, 110, 90, 99, 177, 143, 193, 100],
[90, 125, 82, 13, 57, 71, 123, 38, 72, 59, 82],
[123, 77, 81, 71, 91, 72, 64, 24, 62, 63],
[51, 114, 72, 54, 69, 139, 105, 155, 62],
[70, 25, 22, 52, 90, 56, 105, 16],
[45, 61, 111, 36, 61, 57, 70],
[23, 71, 67, 48, 85, 29],
[74, 89, 69, 107, 36],
[117, 65, 125, 43],
[54, 22, 84],
[60, 44],
[97],
[]]
# number of nodes and list of vertices
n, V = len(dists), set(range(len(dists)))
# distances matrix
c = [[0 if i == j
else dists[i][j-i-1] if j > i
else dists[j][i-j-1]
for j in V] for i in V]
model = Model()
# binary variables indicating if arc (i,j) is used on the route or not
x = [[model.add_var(var_type=BINARY) for j in V] for i in V]
# continuous variable to prevent subtours: each city will have a
# different sequential id in the planned route except the first one
y = [model.add_var() for i in V]
# objective function: minimize the distance
model.objective = minimize(xsum(c[i][j]*x[i][j] for i in V for j in V))
# constraint : leave each city only once
for i in V:
model += xsum(x[i][j] for j in V - {i}) == 1
# constraint : enter each city only once
for i in V:
model += xsum(x[j][i] for j in V - {i}) == 1
# subtour elimination
for (i, j) in product(V - {0}, V - {0}):
if i != j:
model += y[i] - (n+1)*x[i][j] >= y[j]-n
# optimizing
model.optimize()
# checking if a solution was found
if model.num_solutions:
out.write('route with total distance %g found: %s'
% (model.objective_value, places[0]))
nc = 0
while True:
nc = [i for i in V if x[nc][i].x >= 0.99][0]
out.write(' -> %s' % places[nc])
if nc == 0:
break
out.write('\n')
3.3 location：引入SOS变量
import matplotlib.pyplot as plt
from math import sqrt, log
from itertools import product
from mip import Model, xsum, minimize, OptimizationStatus
# possible plants
F = [1, 2, 3, 4, 5, 6]
# possible plant installation positions
pf = {1: (1, 38), 2: (31, 40), 3: (23, 59), 4: (76, 51), 5: (93, 51), 6: (63, 74)}
# maximum plant capacity
c = {1: 1955, 2: 1932, 3: 1987, 4: 1823, 5: 1718, 6: 1742}
# clients
C = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
# position of clients
pc = {1: (94, 10), 2: (57, 26), 3: (74, 44), 4: (27, 51), 5: (78, 30), 6: (23, 30),
7: (20, 72), 8: (3, 27), 9: (5, 39), 10: (51, 1)}
# demands
d = {1: 302, 2: 273, 3: 275, 4: 266, 5: 287, 6: 296, 7: 297, 8: 310, 9: 302, 10: 309}
# plotting possible plant locations
for i, p in pf.items():
plt.scatter((p[0]), (p[1]), marker="^", color="purple", s=50)
plt.text((p[0]), (p[1]), "$f_%d$" % i)
# plotting location of clients
for i, p in pc.items():
plt.scatter((p[0]), (p[1]), marker="o", color="black", s=15)
plt.text((p[0]), (p[1]), "$c_{%d}$" % i)
plt.text((20), (78), "Region 1")
plt.text((70), (78), "Region 2")
plt.plot((50, 50), (0, 80))
dist = {(f, c): round(sqrt((pf[f][0] - pc[c][0]) ** 2 + (pf[f][1] - pc[c][1]) ** 2), 1)
for (f, c) in product(F, C) }
m = Model()
z = {i: m.add_var(ub=c[i]) for i in F} # plant capacity
# Type 1 SOS: only one plant per region
for r in [0, 1]:
# set of plants in region r
Fr = [i for i in F if r * 50 <= pf[i][0] <= 50 + r * 50]
m.add_sos([(z[i], i - 1) for i in Fr], 1)
# amount that plant i will supply to client j
x = {(i, j): m.add_var() for (i, j) in product(F, C)}
# satisfy demand
for j in C:
m += xsum(x[(i, j)] for i in F) == d[j]
# SOS type 2 to model installation costs for each installed plant
y = {i: m.add_var() for i in F}
for f in F:
D = 6 # nr. of discretization points, increase for more precision
v = [c[f] * (v / (D - 1)) for v in range(D)] # points
# non-linear function values for points in v
vn = [0 if k == 0 else 1520 * log(v[k]) for k in range(D)]
# w variables
w = [m.add_var() for v in range(D)]
m += xsum(w) == 1 # convexification
# link to z vars
m += z[f] == xsum(v[k] * w[k] for k in range(D))
# link to y vars associated with non-linear cost
m += y[f] == xsum(vn[k] * w[k] for k in range(D))
m.add_sos([(w[k], v[k]) for k in range(D)], 2)
# plant capacity
for i in F:
m += z[i] >= xsum(x[(i, j)] for j in C)
# objective function
m.objective = minimize(
xsum(dist[i, j] * x[i, j] for (i, j) in product(F, C)) + xsum(y[i] for i in F) )
m.optimize()
plt.savefig("location.pdf")
if m.num_solutions:
print("Solution with cost {} found.".format(m.objective_value))
print("Facilities capacities: {} ".format([z[f].x for f in F]))
print("Facilities cost: {}".format([y[f].x for f in F]))
# plotting allocations
for (i, j) in [(i, j) for (i, j) in product(F, C) if x[(i, j)].x >= 1e-6]:
plt.plot(
(pf[i][0], pc[j][0]), (pf[i][1], pc[j][1]), linestyle="--", color="darkgray"
)
plt.savefig("location-sol.pdf")
4. 自定义约束
动机：放松约束容易求解。每次迭代时需要找到最违背的约束，这也是个优化问题，叫做separation problem。
4.1 Cutting plane，也叫行生成
TSP问题，subtour elimination用的如下形式：
距离矩阵如下：
将第三个约束分离（separate），得到如下解：
添加两条约束：𝑥(𝑑,𝑒)+𝑥(𝑒,𝑑)≤1 以及 𝑥(𝑐,𝑓)+𝑥(𝑓,𝑐)≤1，求解得到如下
增加约束：𝑥(𝑎,𝑔)+𝑥(𝑔,𝑎)+𝑥(𝑎,𝑏)+𝑥(𝑏,𝑎)+𝑥(𝑏,𝑔)+𝑥(𝑔,𝑏)≤2
求解得到结果为261。
接下来的问题是，每次如何获得这些约束？使用 Minimum cut problem in graphs，在python中可以使用networkx min-cut module
from itertools import product
from networkx import minimum_cut, DiGraph
from mip import Model, xsum, BINARY, OptimizationStatus, CutType
N = ["a", "b", "c", "d", "e", "f", "g"]
A = { ("a", "d"): 56, ("d", "a"): 67, ("a", "b"): 49, ("b", "a"): 50,
("f", "c"): 35, ("g", "b"): 35, ("g", "b"): 35, ("b", "g"): 25,
("a", "c"): 80, ("c", "a"): 99, ("e", "f"): 20, ("f", "e"): 20,
("g", "e"): 38, ("e", "g"): 49, ("g", "f"): 37, ("f", "g"): 32,
("b", "e"): 21, ("e", "b"): 30, ("a", "g"): 47, ("g", "a"): 68,
("d", "c"): 37, ("c", "d"): 52, ("d", "e"): 15, ("e", "d"): 20,
("d", "b"): 39, ("b", "d"): 37, ("c", "f"): 35, }
Aout = {n: [a for a in A if a[0] == n] for n in N}
Ain = {n: [a for a in A if a[1] == n] for n in N}
m = Model()
x = {a: m.add_var(name="x({},{})".format(a[0], a[1]), var_type=BINARY) for a in A}
m.objective = xsum(c * x[a] for a, c in A.items())
for n in N:
m += xsum(x[a] for a in Aout[n]) == 1, "out({})".format(n)
m += xsum(x[a] for a in Ain[n]) == 1, "in({})".format(n)
newConstraints = True
while newConstraints:
m.optimize(relax=True)
print("status: {} objective value : {}".format(m.status, m.objective_value))
G = 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 product(N, N) if i != j]:
cut_value, (S, NS) = minimum_cut(G, n1, n2)
if cut_value <= 0.99: # 如果任两点之间不连通，则最小割的值为0
m += (xsum(x[a] for a in A if (a[0] in S and a[1] in S)) <= len(S) - 1)
newConstraints = True
# 没有新的约束时，开始寻找整数约束
if not newConstraints and m.solver_name.lower() == "cbc":
cp = m.generate_cuts([CutType.GOMORY, CutType.MIR,
CutType.ZERO_HALF, CutType.KNAPSACK_COVER])
if cp.cuts:
m += cp
newConstraints = True
4.2 Branch&Cut
使用callback，将割平面与branch&bound结合，称为CGC，会被添加进入cut pool，与求解器自带的cut generator（主要是整数相关）一起添加到LP relaxation问题上。
注意在这个方法中，原始模型就是完备的，Cut的目的是快速改进LP relaxation。在TSP问题中，原始模型是3.2节的模型，cut可以是4.1节新增的cut，本质上是冗余的约束。
下面是操作步骤：
- 定义一个ConstrsGenerator类，实现generate_constrs方法。
- 由于模型有预处理步骤，因此要使用model.translate(self.x)获取新的对应变量
- model.cuts_generator = 新定义的ConstrsGenerator类
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
n = 30 # number of points
V = set(range(n))
seed(0)
p = [(randint(1, 100), randint(1, 100)) for i in V] # coordinates
Arcs = [(i, j) for (i, j) in product(V, V) if i != j]
# distance matrix
c = [[round(sqrt((p[i][0]-p[j][0])**2 + (p[i][1]-p[j][1])**2)) for j in V] for i in V]
model = Model()
# binary variables indicating if arc (i,j) is used on the route or not
x = [[model.add_var(var_type=BINARY) for j in V] for i in V]
# continuous variable to prevent subtours: each city will have a
# different sequential id in the planned route except the first one
y = [model.add_var() for i in V]
# objective function: minimize the distance
model.objective = minimize(xsum(c[i][j]*x[i][j] for (i, j) in Arcs))
# constraint : leave each city only once
for i in V:
model += xsum(x[i][j] for j in V - {i}) == 1
# constraint : enter each city only once
for i in V:
model += xsum(x[j][i] for j in V - {i}) == 1
# (weak) subtour elimination
# subtour elimination
for (i, j) in product(V - {0}, V - {0}):
if i != j:
model += y[i] - (n+1)*x[i][j] >= y[j]-n
# no subtours of size 2
for (i, j) in Arcs:
model += x[i][j] + x[j][i] <= 1
# computing farthest point for each point, these will be checked first for
# isolated subtours
F, G = [], nx.DiGraph()
for (i, j) in Arcs:
G.add_edge(i, j, weight=c[i][j])
for i in V:
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]))
class SubTourCutGenerator(ConstrsGenerator):
"""Class to generate cutting planes for the TSP"""
def __init__(self, Fl: List[Tuple[int, int]], x_):
self.F, self.x = Fl, x_
def generate_constrs(self, model: Model, depth: int = 0, npass: int = 0):
xf, cp, G = model.translate(self.x), CutPool(), nx.DiGraph()
for (u, v) in [(k, l) for (k, l) in product(V, V) 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(V, V) 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)-1
cp.add(cut)
if len(cp.cuts) > 256:
for cut in cp.cuts:
model += cut
return
for cut in cp.cuts:
model += cut
model.cuts_generator = SubTourCutGenerator(F, x)
model.optimize()
print('status: %s route length %g' % (model.status, model.objective_value))
4.3 Lasy Constraints
lazy constraint is a constraint that is only inserted into the model after the first integer solution that violates it is found.
Lasy constraints和cutting plane很类似，但是它是作用在integer solution上的，也就是说：lasy constraints=cutting plane+integer constraints/integer cuts.
下面的代码注释掉了weak subtour elimination和cuts_generator，增加了lazy_constrs_generator。
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
n = 30 # number of points
V = set(range(n))
seed(0)
p = [(randint(1, 100), randint(1, 100)) for i in V] # coordinates
Arcs = [(i, j) for (i, j) in product(V, V) if i != j]
# distance matrix
c = [[round(sqrt((p[i][0]-p[j][0])**2 + (p[i][1]-p[j][1])**2)) for j in V] for i in V]
model = Model()
# binary variables indicating if arc (i,j) is used on the route or not
x = [[model.add_var(var_type=BINARY) for j in V] for i in V]
# continuous variable to prevent subtours: each city will have a
# different sequential id in the planned route except the first one
y = [model.add_var() for i in V]
# objective function: minimize the distance
model.objective = minimize(xsum(c[i][j]*x[i][j] for (i, j) in Arcs))
# constraint : leave each city only once
for i in V:
model += xsum(x[i][j] for j in V - {i}) == 1
# constraint : enter each city only once
for i in V:
model += xsum(x[j][i] for j in V - {i}) == 1
# # (weak) subtour elimination
# # subtour elimination
# for (i, j) in product(V - {0}, V - {0}):
# if i != j:
# model += y[i] - (n+1)*x[i][j] >= y[j]-n
# # no subtours of size 2
# for (i, j) in Arcs:
# model += x[i][j] + x[j][i] <= 1
# computing farthest point for each point, these will be checked first for
# isolated subtours
F, G = [], nx.DiGraph()
for (i, j) in Arcs:
G.add_edge(i, j, weight=c[i][j])
for i in V:
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]))
class SubTourCutGenerator(ConstrsGenerator):
"""Class to generate cutting planes for the TSP"""
def __init__(self, Fl: List[Tuple[int, int]], x_):
self.F, self.x = Fl, x_
def generate_constrs(self, model: Model, depth: int = 0, npass: int = 0):
xf, cp, G = model.translate(self.x), CutPool(), nx.DiGraph()
for (u, v) in [(k, l) for (k, l) in product(V, V) 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(V, V) 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)-1
cp.add(cut)
if len(cp.cuts) > 256:
for cut in cp.cuts:
model += cut
return
for cut in cp.cuts:
model += cut
#model.cuts_generator = SubTourCutGenerator(F, x)
model.lazy_constrs_generator = SubTourCutGenerator(F, x)
model.optimize()
print('status: %s route length %g' % (model.status, model.objective_value))
5. 自定义初始解
一个优秀的初始解可以prune掉很多branch。使用start来定义一个初始解。在TSP案例中，只需要定义x就好了，y会自动计算出来。
from random import shuffle
S=[i for i in range(n)]
shuffle(S)
model.start = [(x[S[k-1]][S[k]], 1.0) for k in range(n)]
2529
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
