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


文章目录

  • 1. 创建模型
  • 2. 定义求解变量
  • 3. 定义目标函数
  • 4. 定义或删除约束条件
  • 5. 设置参数
  • 6. 输出变量值
  • 7. 输出约束条件的对偶值及右端项
  • 8. 检查约束条件
  • 9. 输出 lp, mps 等
  • 10. 模型状态
  • 10. 用命令行运行
  • 11. 几点心得


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

通过python实现混合整数规划模型 python求解混合整数规划_gurobi

# 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. 创建模型

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

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)

上面两个命令是等价的

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

3. 定义目标函数

  • 定义目标函数时，有时会用到 LinExpr（）定义一个表达式，然后用 setObjective() 设置，例如：
final_cash = LinExpr(X[0]+X[1]) # 直接 X[0]+X[1] 也可以
     
# Set objective
m.setObjective(final_cash, GRB.MAXIMIZE)

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

  • python 中的 gurobi 会自动将变量的运算视为一个线性表达式 LinExpr，但需要更新模型，线性表达式才能在模型中生成，用到 update 函数，即：
m.update()
  • LinExpr 也可以生成数组形式，例如：
I = [LinExpr() for i in range(3)]
  • setObjective() 中的表达式 LinExpr 一定要在出现 setObjective() 之前定义好，若之后变动，可能会计算出错
  • 目标函数中不能有包含求解变量的 min 或 max 表达式（此时可以让 min 或 max 表达式等于一个辅助变量，添加到约束条件中）
  • 变量 var 或者线性表达式 LinExpr 可以用 python 的 sum（）相加

4. 定义或删除约束条件

  • 使用 addConstr()，括号内跟约束条件即可。例如：
m.addConstr(x + y <=10)
  • 对于 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")

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

  • 也可以通过 addMVar 定义可以矩阵相乘的求解变量，这样就可以用矩阵乘法 (用到了python 的矩阵乘法符号 @) 定义形如 Ax >= b约束条件了
m.addConstr(A@x >= b)
  • 使用 MVar 定义约束条件 c(Ax+b) >= d 需要注意，此时需要拆开写，因为矩阵相乘符号 @ 不支持括号，要写成：
m.addConstr(c@A@x + c@b >= d)
  • 使用 getConstrs() 取出约束条件，并结合 RHS 可以更新约束条件右端的常数项，例如：
c=m.getConstrs()[0] #将第一个约束条件取出并传递给变量 c
c.RHS = 10  # 更新这个约束条件的常数项
  • 可以用 getConstrs() 函数得到所有的约束条件列表
  • 可以通过 remove() 函数删除某个约束条件，例如，下面的这个例子删除第一个约束条件：
model.remove(models.getConstrs()[0])
model.remove(model.getConstrs()[1:3]) # 移除第 2 与第 3 个约束条件
  • 当删除一个约束条件，又增加一个新的约束条件时，这个新的约束条件的索引自动排到后面

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

设置是否在控制台显示优化信息：

m.Params.LogToConsole = 0

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

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

6. 输出变量值

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

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

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

  • 约束条件的对偶值为
m.getAttr(GRB.Attr.Pi)
  • 约束条件的右端项（RHS，即约束条件中将参数移到右端时，参数的值）
m.getAttr(GRB.Attr.RHS)

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

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 参数查看是否可行，无界，迭代超过限制的信息，官方的说明文档如下：

通过python实现混合整数规划模型 python求解混合整数规划_gurobi_02


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

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

10. 用命令行运行

11. 几点心得

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