python 调用 gurobi 求解数学规划的几点总结，输出 mps，lp_gurobi mps

gurobi 对 python 支持的不错，我已经编写了几个规划求解的例子。每次重新编程时，之前例子里的一些知识点又忘记了，觉得很有必要总结一下。

例如，下面的 python 代码调用 gurobi 求解一个简单的混合整数规划问题：

max ⁡ x + y + 2 z s . t . x + 2 y + 3 z ≤ 4 x + y ≥ 1 x , y , z ∈ { 0 , 1 }

$$\begin{aligned} &\max\quad & x +y +2z\\ &s.t. &\\ &&x + 2 y + 3 z \leq 4\\ && x + y \geq 1\\ && x,y, z\in\{0, 1\} \end{aligned}$$ ​maxs.t.​x+y+2zx+2y+3z≤4x+y≥1x,y,z∈{0,1}​

# This example formulates and solves the following simple MIP model:
#  maximize
#        x +   y + 2 z
#  subject to
#        x + 2 y + 3 z <= 4
#        x +   y       >= 1
#        x, y, z binary

from gurobipy import *

try:

    # Create a new model
    m = Model("mip1")

    # Create variables
    x = m.addVar(vtype=GRB.BINARY, name="x") # default bounds for continuous type is [0, infinite]
    y = m.addVar(vtype=GRB.BINARY, name="y")
    z = m.addVar(vtype=GRB.BINARY, name="z")

    # Set objective
    m.setObjective(x + y + 2 * z, GRB.MAXIMIZE)

    # Add constraint: x + 2 y + 3 z <= 4
    m.addConstr(x + 2 * y + 3 * z <= 4, "c0")

    # Add constraint: x + y >= 1
    m.addConstr(x + y >= 1, "c1")

    m.optimize()

    for v in m.getVars():
        print('%s %g' % (v.varName, v.x))

    print('Obj: %g' % m.objVal)

except GurobiError as e:
    print('Error code ' + str(e.errno) + ": " + str(e))

except AttributeError:
    print('Encountered an attribute error')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

几点总结：

gurobi 的申请方式参看：http://www.gurobi.cn/NewsView1.Asp?id=4

1. 创建模型

• m = Model(‘model_name’) 引号里面是自己起的模型名
• 另一种方法，直接用 read() 读取 mps 或 lp 格式的文件。例如：

model = read('test.mps')
1

2. 定义求解变量

• 使用 m.addVar(vtype = GRB.CONTINUOUS)，小括号里面是求解变量的类型，gurobi 中 GRB.CONTINUOUS 表示大于等于零的连续型实数，GRB.INTEGER 表示整数型，GRB.BINARY 表示0-1型。
• 小括号还可以跟 lb, ub 表示变量的上下界， 跟 obj 表示目标函数中该变量的系数（此时就不用专门再定义目标函数了）。
• 小括号里还可以用 name 参数跟一个字符串给变量命名，这样输出模型时更加清晰
• 批量定义变量可以用 m.addVars()，或者像定义数组那样使用 m.addVar()，例如：

X = [m.addVar(vtype = GRB.CONTINUOUS) for t in range(3)]

X = m.addVars(3, vtype = GRB.CONTINUOUS)
1
2
3

上面两个命令是等价的

• 还可以通过 addMVar 定义一个变量的数组，这样在添加约束条件时，就可以通过矩阵乘法添加了，不用一个一个条件添进去

3. 定义目标函数

• 定义目标函数时，有时会用到 LinExpr（）定义一个表达式，然后用 setObjective() 设置，例如：

final_cash = LinExpr(X[0]+X[1]) # 直接 X[0]+X[1] 也可以
     
# Set objective
m.setObjective(final_cash, GRB.MAXIMIZE)
1
2
3
4

若目标函数的表达式比较简单，也可以直接放到 setObjective 里。

• python 中的 gurobi 会自动将变量的运算视为一个线性表达式 LinExpr，但需要更新模型，线性表达式才能在模型中生成，用到 update 函数，即：

m.update()
1

• LinExpr 也可以生成数组形式，例如：

I = [LinExpr() for i in range(3)]
1

• setObjective() 中的表达式 LinExpr 一定要在出现 setObjective() 之前定义好，若之后变动，可能会计算出错

• 目标函数中不能有包含求解变量的 min 或 max 表达式（此时可以让 min 或 max 表达式等于一个辅助变量，添加到约束条件中）

• 变量 var 或者线性表达式 LinExpr 可以用 python 的 sum（）相加

4. 定义或删除约束条件，约束条件系数

• 使用 addConstr()，括号内跟约束条件即可。例如：

m.addConstr(x + y <=10)
1

• 对于 min 或 max 表达式的约束条件，可以用 addGenConstrMax()、addGenConstrMin()，或者在 addConstr() 里面使用 max_()、min_()。例如：

  # x5 = max(x1, x3, x4, 2.0)
  m.addGenConstrMax(x5, [x1, x3, x4], 2.0, "maxconstr")

  # alternative form
  m.addGenConstrMax(x5, [x1, x3, x4, 2.0], name="maxconstr")

  # overloaded forms
  m.addConstr(x5 == max_([x1, x3, x4, 2.0]), name="maxconstr")
  m.addConstr(x5 == max_(x1, x3, x4, 2.0), name="maxconstr")
1
2
3
4
5
6
7
8
9

不过我发现，当遇到数组变量时， max_ 或 min_ 可能会出错，所以一般还是用 addGenConstrMax()、addGenConstrMin() 来处理 min 或 max 表达式

• 也可以通过 addMVar 定义可以矩阵相乘的求解变量，这样就可以用矩阵乘法 (用到了python 的矩阵乘法符号 @) 定义形如 Ax >= b约束条件了

m.addConstr(A@x >= b)
1

• 使用 MVar 定义约束条件 c(Ax+b) >= d 需要注意，此时需要拆开写，因为矩阵相乘符号 @ 不支持括号，要写成：

m.addConstr(c@A@x + c@b >= d)
1

• 使用 getConstrs() 取出约束条件，并结合 RHS 可以更新约束条件右端的常数项，例如：

c=m.getConstrs()[0] #将第一个约束条件取出并传递给变量 c
c.RHS = 10  # 更新这个约束条件的常数项
1
2

• 可以用 getConstrs() 函数得到所有的约束条件
• 可以通过 remove() 函数删除某个约束条件，例如，下面的这个例子删除第一个约束条件：

model.remove(models.getConstrs()[0])
model.remove(model.getConstrs()[1:3]) # 移除第 2 与第 3 个约束条件
1
2

• 当删除一个约束条件，又增加一个新的约束条件时，这个新的约束条件的索引自动排到后面

• 用 getA() 得到线性约束条件的系数矩阵

A = model.getA()
1

5. 设置参数

• 可以通过 Params 设置，例如设置求解时间上限等：

m.Params.timeLimit = 100.0 # 等价于 m.setParam('timeLimit', 100)
m.Params.InfUnbdInfo = 1  # Determines whether simplex (and crossover) will compute additional information when a model is determined to be infeasible or unbounded
1
2

• 设置是否在运行过程中显示优化信息 (不在 ipython 的控制台显示信息，能够节约运行时间；但如果通过命令行运行程序，速度差别不大)：

m.Params.LogToConsole = 0
1

• 设置线性规划的求解方法：

 m.params.Method = 1  # 使用对偶单纯形法
 m.params.Method = 0  # 使用原始单纯形法(迭代慢，但占内存小）
 m.params.Method = 2  # 使用内点法（gurobi称作barrier法）
1
2
3

6. 输出目标函数值和变量值

• 输出目标函数值用 .objVal
• 输出变量值有时用 .getValue()，有时用 .X，
• 若其中一个求解变量本质上是其他变量的表达式(LinExpr)，用 getValue()，否则用 X。例如：Q.X，其中 Q 为模型中的求解变量； I.getValue()，其中 I 为求解变量 Q 的表达式。

7. 输出约束条件的对偶值及右端项

对于 gurobi 的一个线性规划模型 m

• 约束条件的对偶值为

m.getAttr(GRB.Attr.Pi)
1

• 约束条件的右端项（RHS，即约束条件中将参数移到右端时，参数的值）

m.getAttr(GRB.Attr.RHS)
1

若约束条件有多个，上面的两个命令得出的都是列表值

8. 检查约束条件

• 使用 m.computeIIS() 检查不可行的约束条件（模型得不到可行解时，才能用这个函数），然后通过命令 m.write(“model.ilp”)，可以进一步将不可行的详细信息输出到 model.ilp 这个文件中
• 使用 m.feasRelax() 通过松弛最少的不可行约束条件，得到一个可行解（模型得不到可行解时，才能用这个函数）
• 这两个约束条件检查一般放在求解 m.optimize() 之前，或者模型m.optimize() 后不可行时也可以用

9. 输出 lp, mps 等

• 可以用 m.write() 输出模型，括号内可以跟多种文件格式，例如：lp，mps，或者 ilp 来输出 IIS 等。
• 格式 sol 可以将线性规划的求解结果输出到一个 sol 文件中（要在模型求解后）
• 格式 dlp 可以输出纯线性规划模型的对偶模型（输出的对偶模型都为标准型，但 gurobi 求解出的约束条件的对偶值 PI 却是原模型对应的 PI 值，并不是标准型约束条件的对偶值）

10. 模型状态

gurobi 的 模型可以通过访问 Status 参数查看是否可行，无界，迭代超过限制的信息，官方的说明文档如下：

例如，下面通过访问模型模型不可行时，计算 IIS：

if m.Status == 3: # model is infeasible 
    print("Model is infeasible")
    m.computeIIS()
    m.write("model.ilp")
1
2
3
4

10. 用命令行运行

• gurobi 也支持用命令行运行，调出命令行窗口，使用 gurobi_cl 跟上不同的指令，具体可参看：
  https://www.gurobi.com/documentation/9.0/refman/grb_command_line_tool.html

11. 几点心得

• gurobi 会自动调用多线程进行并行计算，所以对于大规模问题，电脑内存一定要大
• gurobi 官方说对于线性规划，单纯形法的计算效果一般最好。我之前一直以为内点法最好
• 在一个循环中不断构建新的 model， 与在循环外构造好 model并更新约束条件，计算速度差别不大

转载于个人公众号：Python 数据科学与数学建模