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

心态与习惯

已于 2024-05-29 01:11:55 修改

阅读量2.8w

点赞数 49

分类专栏： Python Cplex/Gurobi 文章标签： python gurobi 数学规划混合整数规划线性规划

于 2020-09-03 11:27:45 首次发布

本文链接：https://blog.csdn.net/robert_chen1988/article/details/107085828

版权

Python 同时被 2 个专栏收录

58 篇文章 36 订阅

订阅专栏

Cplex/Gurobi

22 篇文章 40 订阅

订阅专栏

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

文章目录

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

$\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}$

# 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')

几点总结：

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

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 个约束条件

当删除一个约束条件，又增加一个新的约束条件时，这个新的约束条件的索引自动排到后面
用 getA() 得到线性约束条件的系数矩阵

A = model.getA()

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

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

m.Params.LogToConsole = 0

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

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

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

输出目标函数值用 .objVal
输出变量值有时用 .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 参数查看是否可行，无界，迭代超过限制的信息，官方的说明文档如下：
在这里插入图片描述
例如，下面通过访问模型模型不可行时，计算 IIS：

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

10. 用命令行运行

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

11. 几点心得

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

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

在这里插入图片描述

心态与习惯

关注

49
点赞
踩
307

收藏

觉得还不错? 一键收藏
打赏
41
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫

专栏目录