Gurobi +Python 高效数学规划及建模实例

Gurobi 可以解决的数学问题：

线性问题（Linear problems）

二次型目标问题（Quadratic problems）

混合整数线性和二次型问题（Mixed integer linear and quadratic problems）

等等

2 快速入门

在利用 Python+Gurobi 建立数学规划模型时，通常会按照设置变量、更新变量空间、设置目标函数、设置约束条件、执行最优化的顺序进行。

python复制代码import gurobipy

# 创建模型

MODEL = gurobipy.Model()

# 创建变量

X = MODEL.addVar(vtype=gurobipy.GRB.INTEGER,name="X")

# 更新变量环境

MODEL.update()

# 创建目标函数

MODEL.setObjective('目标函数表达式', gurobipy.GRB.MINIMIZE)

# 创建约束条件

MODEL.addConstr('约束表达式，逻辑运算')

# 执行线性规划模型

MODEL.optimize()

# 输出模型结果

print("Obj:", MODEL.objVal)

for x in X:

    print(f"{x.varName}：{round(x.X,3)}")

2.1 辅助函数

2.1.1 列表推导式/列表解析式

列表推导式基于迭代器进行构造，效率高，常用于遍历、过滤、重复计算，快速生成满足要求的序列。

形式上使用中括号作为定界符

列表推导式以 if - for 为主，主要有以下两种嵌套方式：

python复制代码# if 在后 —— 过滤元素

[expression for exp1 in sequence1 if condition1

            for exp2 in sequence2 if condition2

            ...

            for exp2 in sequence2 if conditionn]

# 符合 condition 条件的值会被储存在列表中，

# 不符合的不保存 (类似于 pass 和 filter 函数)

python复制代码# if 在前 —— 分段映射

[expression1 if condition1 else expression2 for exp1 in sequence1

                                            for exp2 in sequence2

                                            ...

                                            for exp2 in sequence2]

python复制代码# 嵌套循环

[f'{cl}班学生{num}' for cl in list('ABCDE') for num in range(1,6)] # for 从左到右循环

>>> ['A班学生1', 'A班学生2', 'A班学生3', 'A班学生4', 'A班学生5', 'B班学生1', 'B班学生2', 'B班学生3', 'B班学生4', 'B班学生5', 'C班学生1', 'C班学生2', 'C班学生3', 'C班学生4', 'C班学生5', 'D班学生1', 'D班学生2', 'D班学生3', 'D班学生4', 'D班学生5', 'E班学生1', 'E班学生2', 'E班学生3', 'E班学生4', 'E班学生5']

以上几种可以混合使用，十分高效。

举例：

列出当前文件夹下的所有 Python 源文件

python复制代码[filename for filename in os.listdir() if filename.endswith(('.py', '.pyw'))]

# for 循环遍历 os.listdir() 中的所有的文件，if 过滤以 .py, .pyw 结尾的文件

2.1.2 quicksum()

在 gurobipy 模块中，quicksum 相当于 sum 函数及求和符号：

∑j∈Jxi,j≤5∀i∈I\sum_{j \in J} x_{i, j} \leq 5 \quad \forall i \in Ij∈J∑xi,j≤5∀i∈I

python复制代码for i in I:

    MODEL.addConstr(quicksum(x[i,j] for j in J)<=5)

## 用 sum 时

for i in I:

    MODEL.addConstr(sum(x[i,j] for j in J)<=5)

要注意: quicksum() 是 gurobi 推荐的用法, 它与 sum() 是一致的，但是出现字典时例外。

2.1.3 multidict()

扩展字典，便于处理同一个对象的不同属性约束。

python复制代码EMPLOYEE, MIN, MAX, COST, START, END = gurobipy.multidict({

"SMITH"   : [6, 8, 30, 6, 20], "JOHNSON": [6, 8, 50, 0, 24], 'WILLIAMS': [6, 8, 30, 0, 24],

'JONES'   : [6, 8, 30, 0, 24], 'BROWN': [6, 8, 40, 0, 24], 'DAVIS': [6, 8, 50, 0, 24],

'MILLER'  : [6, 8, 45, 6, 18], 'WILSON': [6, 8, 30, 0, 24], 'MOORE': [6, 8, 35, 0, 24],

'TAYLOR'  : [6, 8, 40, 0, 24], 'ANDERSON': [2, 3, 60, 0, 6], 'THOMAS': [2, 4, 40, 0, 24],

'JACKSON' : [2, 4, 60, 8, 16], 'WHITE': [2, 6, 55, 0, 24], 'HARRIS': [2, 6, 45, 0, 24],

'MARTIN'  : [2, 3, 40, 0, 24], 'THOMPSON': [2, 5, 50, 12, 24], 'GARCIA': [2, 4, 50, 0, 24],

'MARTINEZ': [2, 4, 40, 0, 24], 'ROBINSON': [2, 5, 50, 0, 24]})

EMPLOYEE 变成了list，而 MIN、MAX 等等则变成了字典

EMPLOYEE = [“SMITH”, "JOHNSON",...]

MIN = {“SMITH”: 6, "JOHNSON": 6,...}

2.1.4 tuplelist()

扩展列表元组，可以使用通配符筛选变量组。若使用 tuplelist 创建变量：

python复制代码Cities= [('A','B'), ('A','C'), ('B','C'),('B','D'),('C','D')]

Routes = gurobipy.tuplelist(Cities)

Routes.select('A','*') # 选出所有第一个, 素为 "A" 的, 组

# 此外, 对于addVars创建的变量, x.select("*", i, j) 也可以进行筛选

# 可以将 '*' 理解为切片操作(列表的 list[:])

2.1.5 prod() 和 sum() 下标聚合

∑i∑jxijcij\sum_{i} \sum_{j} x_{i j} c_{i j}i∑j∑xijcij

python复制代码gurobipy.quicksum(cost[i,j] * x[i,j] for i,j in arcs)

# 等价于 x.prod(cost)

# 效果为对应元素相乘再相加

2.2 添加决策变量 Model.addVar() 和 Model.addVars()

2.2.1 创建一个变量

x = MODEL.addVar(lb=0.0, ub=gurobipy.GRB.INFINITY, vtype=gurobipy.GRB.CONTINUOUS, name="")

lb = 0.0：变量的下界，默认为 0

ub = GRB.INFINITY：变量的上界，默认为无穷大

vtype = GRB.CONTINUOUS ：变量的类型，默认为连续型，可改为 GRB.BINARY 0-1变量，GRB.INTEGER 整型

name = ""：变量名，默认为空

使用方法：

python复制代码x1 = MODEL.addVar(lb=0, ub=1, name="x1")

2.2.2 创建多个变量

x = MODEL.addVars(*indexes, lb=0, ub=gurobipy.GRB.INFINITY, vtype=gurobipy.GRB.CONTINUOUS, name="")

python复制代码x = MODEL.addVars(3, 4, 5, vtype=gurobipy.GRB.BINARY, name="C") # 创建 3*4*5 个变量，使用 x[1,2,3] 进行访问

# lb，ub，vtype 可以单独设置(同样维度数据)，也可以全部设置(单个值)

# 此外，若使用 m.addVars(iterable1,iterable2,iterable3) ，相当于笛卡尔积，参考案例 4.2

# 这种创建方法可以使用通配符命令，能简化代码

2.2.3 添加完变量后需要 MODEL.update() 更新变量空间

2.3 添加目标函数 Model.setObjective() 和 Model.setObjectiveN()

2.3.1 单目标优化

MODEL.setObjective(expression, sense=None)

expression： 表达式，可以是一次或二次函数类型

sense：求解类型，可以是GRB.MINIMIZE 或 GRB.MAXIMIZE

python复制代码MODEL.setObjective(8 * x1 + 10 * x2 + 7 * x3 + 6 * x4 + 11 * x5 + 9 * x6, gurobipy.GRB.MINIMIZE)

# 执行完成最优化后，可以通过 MODEL.objVal 获取目标函数值(如果存在的话)

2.3.2 多目标优化(默认最小值)

expression: 表达式，可以是一次或二次函数类型

index: 目标函数对应的序号 (默认 0，1，2，…), 以 index=0 作为目标函数的值, 其余值需要另外设置参数

priority: 分层序列法多目标决策优先级(整数值), 值越大优先级越高

weight: 线性加权多目标决策权重(在优先级相同时发挥作用)

abstol: 分层序列法多目标决策时允许的目标函数值最大的降低量

reltol: 分层序列法多目标决策时允许的目标函数值最大的降低比率 reltol*|目标函数值|

多目标优化案例

python复制代码# 1. 合成型 ######################################################

# Obj1 = x + y          weight = 1

# Obj2 = x - 5 * y      weight = -2

# MODEL.setObjectiveN(x + y, index=0, weight=1, name='obj1')

# MODEL.setObjectiveN(x -5 * y, index=1, weight=-2, name='obj2')

# 即转化为：(x + y) - 2 * (x - 5 * y) = - x + 11 * y

# 2. 分层优化型(目标约束时使用) #####################################

# Obj1 = x + y          priority = 5

# Obj2 = x - 5 * y      priority = 1

# MODEL.setObjectiveN(x + y, index=0, priority=5, name='obj1')

# MODEL.setObjectiveN(x -5 * y, index=1, priority=1, name='obj2')

# 即转化为：先优化 Obj1，再优化 Obj2（按照 priority 的大小关系）

 

# 3. 混合优化型 ##################################################

# MODEL.setObjectiveN(x + y, index=0, weight=1, priority=5, name='obj1')

# MODEL.setObjectiveN(x -5 * y, index=1, weight=-2, priority=1, name='obj2')

# 则 先优化 Obj1 再优化 Obj2 最后相加作为目标值

 

# 4. 执行最优化结束后，获取目标函数值  ###############################

# MODEL.setParam(gurobipy.GRB.Param.ObjNumber, i)   # 第 i 个目标

# print(MODEL.ObjNVal)  # 打印第 i 个值

2.4 添加约束条件 Model.addConstr() 和 Model.addConstrs()

2.4.1 创建一个常规一次/二次/等式约束

MODEL.addConstr(expression, name="")

expression: 布尔表达式，可以是一次或二次函数类型

name: 约束式的名称

使用方法：

python复制代码MODEL.addConstr(12 * x1 + 9 * x2 + 25 * x3 + 20 * x4 + 17 * x5 + 13 * x6 >= 60, "c0")

# 注意，此处是调用 c 库实现的，若有 a < b < c，只能拆开写成 a < b，和 b < c

2.4.2 创建多个常规一次/二次/等式约束

MODEL.addConstrs(expressions, name="")

addConstrs 用于添加多个约束条件，参数与 addConstr 类似，最大的优势是可以将外层的 for 循环转化为内部迭代器形式

∑j=07xij≤1,∀i=0,1,⋯ ,19xij=0 or 1\begin{aligned}

\sum_{j=0}^{7} x_{i j} \leq 1, \forall i=0,1, \cdots, 19 \\

x_{i j}=0 & \text { or } 1

\end{aligned}j=0∑7xij≤1,∀i=0,1,⋯,19xij=0 or 1

python复制代码x = MODEL.addVars(20, 8, vtype=gurobipy.GRB.BINARY)

# 写法 1

for i in range(20):

    MODEL.addConstr(x.sum(i, "*") <= 1)

# 写法 2

MODEL.addConstrs(x.sum(i, "*") <= 1 for i in range(20))

# 可以将 '*' 理解为切片操作(列表的 list[:])

2.4.3 创建一个范围约束

MODEL.addRange(expression, min_value, max_value, name="")

表达式 min_value<=expression<=max_value 的简写, 但这里的 min_value, max_value 必须是具体的实数值, 不能含有变量

2.4.4 创建一个指示变量约束

MODEL.addGenConstrIndicator(binvar, binval, expression, name="")

指示变量 binvar 的值取 binval 时, 进行约束 expression

案例:

x1,x2,x3=0 或 ≥80x_{1}, x_{2}, x_{3}=0 \text { 或 } \geq 80x1,x2,x3=0 或 ≥80

方法一: 构造指示变量 , 则上述约束转化为

80⋅yi≤xi≤M⋅yi(M 是一个很大的数, 可取 1000)80 \cdot y_{i} \leq x_{i} \leq M \cdot y_{i} \quad(M \text { 是一个很大的数, 可取 } 1000)80⋅yi≤xi≤M⋅yi(M 是一个很大的数, 可取 1000)

方法二: 转为二次约束, 若矩阵非正定矩阵, 则无法求解

xi(xi−80)≥0 x_{i}( x_{i}-80) \geq 0xi(xi−80)≥0

案例: 如果生产某一类型汽车 (x[i] > 0), 则至少要生产 80 辆

python复制代码# %% 方法一

for i in range(3):

MODEL.addGenConstrIndicator(y[i + 1], 0, x[i + 1] >= 80)

MODEL.addGenConstrIndicator(y[i + 1], 1, x[i + 1] == 0)

    

# 以上代码等价于

for i in range(3):

MODEL.addConstr(x[i + 1] >= 80 * y[i + 1])

MODEL.addConstr(x[i + 1] <= 1000 * y[i + 1])

2.5 执行最优化 Model.optimize()

python复制代码MODEL.Params.LogToConsole=True # 显示求解过程

MODEL.Params.MIPGap=0.0001 # 百分比界差

MODEL.Params.TimeLimit=100 # 限制求解时间为 100s

MODEL.optimize()

2.6查看模型结果

2.6.1 模型是否取得最优解

MODEL.status == gurobipy.GRB.Status.OPTIMAL

True 最优

False 非最优, 在TimeLimit模式下, 不能用这个方法判断最优解

2.6.2 单目标优化 —— 查看目标函数值

python复制代码# 查看单目标规划模型的目标函数值

print("Optimal Objective Value", MODEL.objVal)

2.6.3 多目标优化 —— 查看目标函数值

python复制代码# 查看多目标规划模型的目标函数值

for i in range(MODEL.NumObj):

MODEL.setParam(gurobipy.GRB.Param.ObjNumber, i)

print(f"Obj {i+1} = {MODEL.ObjNVal}")

2.6.4 查看变量取值

python复制代码# 查看变量取值，这个方法用的很少，请看第 4 部分案例

for var in MODEL.getVars():

   print(f"{var.varName}: {round(var.X, 3)}")

2.6.5 LP 问题灵敏度分析

通过 MODEL.getVars() 得到的变量 var

使用 var.X 可以获取变量值, var.RC 可以获取 Reduced Cost；

var.Obj   可以得到在目标函数中的系数

var.SAObjLow 可以得到 Allowable Minimize

var.SAObjUp 可以得到 Allowable Maximize

通过 MODEL.getConstrs() 得到的约束式 Constr

Constr.Slack 可以得到 Slack or Surplus

Constr.pi 可以得到 Dual Price

Constr.RHS 可以得到约束式右侧常值

Constr.SARHSLow 可以得到 Allowable Minimize

Constr.SARHSUp 可以得到  Allowable Maximize

3 实例：一般线性规划问题

min⁡Z=8x1+10x2+7x3+6x4+11x5+9x6 s.t. 12x1+9x+25x3+20x4+17x5+13x6≥6035x1+42x2+18x3+31x4+56x5+49x6≥15037x1+53x2+28x3+24x4+29x5+20x6≥1250≤xj≤1,j=1,2,⋯ ,6\begin{array}{ll}

\min & Z=8 x_{1}+10 x_{2}+7 x_{3}+6 x_{4}+11 x_{5}+9 x_{6} \\

\text { s.t. } & 12 x_{1}+9 x+25 x_{3}+20 x_{4}+17 x_{5}+13 x_{6} \geq 60 \\

& 35 x_{1}+42 x_{2}+18 x_{3}+31 x_{4}+56 x_{5}+49 x_{6} \geq 150 \\

& 37 x_{1}+53 x_{2}+28 x_{3}+24 x_{4}+29 x_{5}+20 x_{6} \geq 125 \\

& 0 \leq x_{j} \leq 1, j=1,2, \cdots, 6

\end{array}min s.t. Z=8x1+10x2+7x3+6x4+11x5+9x612x1+9x+25x3+20x4+17x5+13x6≥6035x1+42x2+18x3+31x4+56x5+49x6≥15037x1+53x2+28x3+24x4+29x5+20x6≥1250≤xj≤1,j=1,2,⋯,6

python复制代码import gurobipy

# 创建模型

c = [8, 10, 7, 6, 11, 9]

p = [[12, 9, 25, 20, 17, 13],

[35, 42, 18, 31, 56, 49],

[37, 53, 28, 24, 29, 20]]

r = [60, 150, 125]

MODEL = gurobipy.Model("Example")

# 创建变量

x = MODEL.addVars(6, lb=0, ub=1, name='x')

# 更新变量环境

MODEL.update()

# 创建目标函数

MODEL.setObjective(x.prod(c), gurobipy.GRB.MINIMIZE)

# 创建约束条件

MODEL.addConstrs(x.prod(p[i]) >= r[i] for i in range(3))

# 执行线性规划模型

MODEL.optimize()

print("Obj:", MODEL.objVal)

for v in MODEL.getVars():

print(f"{v.varName}：{round(v.x,3)}")

上一条：MAC安装gurobi

下一条：Gurobi求解器使用实例