一、基本了解

1、求解范围

1. 连续问题、整数问题、线性和二次凸问题、二次非凸问题、广义非线性问题等
2. 广义非线性问题——广义函数约束
   函数形式可以是高阶多项式、对数、指数、三角函数等非线性函数，那么Gurobi 会对这些函数自动分段线性化进行近似，用户可以通过参数来平衡近似的精度和速度。这样我们就允许在传统的线性和二次型模型中看到这些非线性函数，大大扩展了Gurobi求解器的适用范围，而且这些结果也具有全局最优性。

2、求解速度

线性规划速度 > 0-1 整数线性规划 > 整数线性规划 > 混合整数线性规划 > 非线性规划
在现实问题中，应用最广泛的数学规划问题类型是混合整数线性规划（MILP）。

3、学习资料

• 安装目录 examples, docs 目录下有完整的范例和使用手册、参考手册
• 官网 www.gurobi.com
  1. 在线手册：http://www.gurobi.com/documentation/
  2. 视频：http://www.gurobi.com/resources/seminars-andvideos/seminars-videos
  3. 在线课程：http://www.gurobi.com/academia/for-online-courses
• 中文网站 www.gurobi.cn

二、快速入门

1.常规步骤

MODEL = gurobipy.Model()#创建模型

x = MODEL.addVars(创建变量, vtype=gurobipy.GRB.CONTINUOUS, name="X")#创建变量

MODEL.update()#更新变量环境

MODEL.setObjective(创建目标函数, sense=gurobipy.GRB.MINIMIZE)#创建目标函数

MODEL.addConstrs(创建约束条件) #创建约束条件

MODEL.optimize() #执行最优化

2.基础模型

3.详细步骤

3.1 实例化求解器

MODEL = gurobipy.Model(name="")
name 为模块名默认为空

import gurobipy as GRB
MODEL = GRB.Model("Model_name")   #创建模型

3.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 = "" 	#变量名，默认为空

x1 = MODEL.addVar(lb=0, ub=1, name="x1")

(2) 多个变量

x = MODEL.addVars(*indexes, lb=0, ub=gurobipy.GRB.INFINITY, vtype=gurobipy.GRB.CONTINUOUS, name="")
# *indexes 不能含有中文字符 索引的集合

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

3.3 更新变量空间

MODEL.update()

3.4 添加目标函数

(1) 单目标函数

MODEL.setObjective(expression, sense=None)
expression #表达式，可以是一次或二次函数类型
sense #求解类型，可以是 GRB.MINIMIZE(最小值) 或 GRB.MAXIMIZE(最大值)
MODEL.setObjective(8 * x1 + 10 * x2 + 7 * x3 + 6 * x4 + 11 * x5 + 9 * x6, GRB.GRB.MINIMIZE)

(2) 多目标函数

MODEL.setObjectiveN(expression, index, priority=0, weight=1.0, abstol=0, reltol=0, name="")
expression 	# 表达式，可以是一次或二次函数类型
index	# 目标函数对应的序号 (默认 0，1，2，…), 以 index=0 作为目标函数的值, 其余值需要另外设置参数
priority	# 分层序列法多目标决策优先级(整数值), 值越大优先级越高
weight	#线性加权多目标决策权重(在优先级相同时发挥作用)
abstol 	# 分层序列法多目标决策时允许的目标函数值最大的降低量
reltol	#分层序列法多目标决策时允许的目标函数值最大的降低比率 reltol*|目标函数值|

注意： 所有的目标函数都为线性的,并且目标函数的优化方向一致(全部最大化或全部最小化)。可以通过乘以 -1 实现不同的优化方向。

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 的大小关系）

注意： Gurobi 按照优先级大小优化(先优化Obj1),若没有设定abstol或reltol,在优化低优先级目标时,不会改变高优先级的目标值。假设Obj1=10,在优化 Obj2 时只能在使得 Obj1=10 的所有解中挑选最优解

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 最后相加作为目标值, 由于优先级不同, weight不会起作用

4.获取目标函数值

通过参数 ObjNumber 选择特定的目标,进而获得对应的目标函数值。

MODEL.setParam(gurobipy.GRB.Param.ObjNumber, i) # 第 i 个目标
print(MODEL.ObjNVal)  # 打印第 i 个值

for i in range(model.NumObj):
model.setParam(GRB.Param.ObjNumber, i)
print('Obj%d = ' %(i+1), model.ObjNVal)

(3) 分段目标函数

setPWLObj(var,x,y)
• var	 #指定变量的目标函数是分段线性
• x 	#定义分段线性目标函数的点的横坐标值(非减序列)
• y 	#定义分段线性目标函数的点的纵坐标值

例1：

m.setPWLObj(x, [0, 2, 5], [0, 2, 3])
m.setAttr(GRB.Attr.ModelSense, -1)

注意： 对一些非线性模型，可以使用这一功能去线性逼近。
例2：

m.setPWLObj(x, [0,1/4, 1,7/4, 2], [3,41/9, 2,41/9,3])

3.5 添加约束条件

(1) 创建正常约束

MODEL.addConstr(expression, name="")
expression # 布尔表达式，可以是一次或二次函数类型
name	# 约束式的名称

# 单个约束
MODEL.addConstr(12 * x1 + 9 * x2 + 25 * x3 + 20 * x4 + 17 * x5 + 13 * x6 >= 60, "c0")

# 多个约束
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))
# 写法 2 的等价写法: 
MODEL.addConstrs(sum(x[i, j] for j in range(8) <= 1 for i in range(20))

(2) 创建广义约束

1. 最大值Max

addGenConstrMax(resvar, vars, constant, name )
• resvar 变量(x = max(x1,x2,10))
• vars 一组变量(可以包含常数)
• constant 常数
• name 广义约束名称

例如： z = max(x, y, 3)

m.addGenConstrMax(z, [x, y], 3, "maxconstr")
m.addGenConstrMax(z, [x, y, 3], name="maxconstr")

换成一般的约束表达方式：

m.addConstr(z == max_([x, y, 3]), "maxconstr")
m.addConstr(z == max_(x, y, 3), "maxconstr")

2. 最小值Min

addGenConstrMin( resvar, vars, constant, name ) 
• resvar 变量(x = min(x1,x2,10))
• vars 一组变量(可以包含常数)
• constant 常数
• name 广义约束名称

例如: z = min(x, y, 3)

m.addGenConstrMin(z, [x, y], 3, "minconstr")
m.addGenConstrMin(z, [x, y, 3], name="minconstr")

换成一般的约束表达方式：

m.addConstr(z == min_([x, y, 3]), "minconstr")
m.addConstr(z == min_(x, y, 3), "minconstr")

3. 绝对值Abs

addGenConstrAbs( resvar, argvars, name ) 
• resvar 变量
• argvar 变量
• name 广义约束名称

例如: x = |y|

m.addGenConstrAbs(x, y, "absconstr")

换成一般的约束表达方式：

m.addConstr(x == abs_(y), "absconstr")

4. 和

一组变量的值全等于1,则取1,否则取0。

addGenConstrAnd( resvar, vars, name ) 
• resvar 变量
• vars 一组变量
• name 广义约束名称

例如: x = 1且y = 1,那么z = 1,否则 z = 0

m.addGenConstrAnd(z, [x,y], "andconstr")

换成一般的约束表达方式：

m.addConstr(z == and_(x, y), "andconstr")

5. 与

一组变量的值有一个等于1,则取1,否则取0。

addGenConstrOr( resvar, vars, name ) 
• resvar 变量
• vars 一组变量
• name 广义约束名称

例如: x = 0且y = 0,那么z = 0,否则 z = 1

m.addGenConstrOr(z, [x,y], "orconstr")

换成一般的约束表达方式：

m.addConstr(z == or_(x, y), "orconstr")

6. 范围约束

MODEL.addRange(expression, min_value, max_value, name="")
# 表达式 min_value<=expression<=max_value 的简写, 但这里的 min_value, max_value 必须是具体的实数值, 不能含有变量

例如: 5 <= x + y + z <= 10

m.addRange(x+y+z, 5, 10, "c") 

换成一般的约束表达方式：

m.addConstr(x+y+z==[5,12])

7. 指示变量

MODEL.addGenConstrIndicator(binvar, binval, expression, name="")
# 指示变量 binvar 的值取 binval 时, 进行约束 expression
# 指示变量的值为1,约束成立,否则约束可以被违反。

# 案例: 如果生产某一类型汽车 (x[i] > 0), 则至少要生产 80 辆
# %% 方法一
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: 如果 z = 1,则 x+y <= 4

m.addGenConstrIndicator(z, True, x + y, GRB.LESS_EQUAL, 4, 'indicator')

换成一般的约束表达方式：

m.addConstr((z == 1) >>(x + y<= 4))

3.6 执行最优化

MODEL.Params.LogToConsole=True # 显示求解过程
MODEL.Params.MIPGap=0.0001 # 百分比界差
MODEL.Params.TimeLimit=100 # 限制求解时间为 100s
MODEL.Params.Presolve = -1 # 预处理程度, 0关闭,1保守,2激进
MODEL.Params.MIPFocus = 0 # 求解侧重点. 1快速找到可行解, 2证明最有, 3侧重边界提升, 0均衡搜索
MODEL.Params.SolutionLimit = inf # 求解数量, 默认求所有解, 比较出最优的结果, 只需要可行解时可以设置该参数为1
MODEL.Params.NonConvex = 1  # 默认求解器，改为 2 时可以解决非凸二次优化问题
​
MODEL.optimize()

整数规划/混合整数规划问题评价准则: 百分比界差

OPT 为当前的最优解，LP 为去掉整数约東后的松弛最优解百分比界差是数学规划领域常用的评价参数，具体可视为在去除问题参数整数取值的限制后，问题松弛解与最优解之间的绝对差距

3.7 查看模型结果

#单目标获得目标值和单个变量值
print("Obj = ",m.ObjVal)
for i in range(N):
print(y[i].VarName,' = ',y[i].x)
# 查看多个变量取值
for var in MODEL.getVars():
   print(f"{var.varName}: {round(var.X, 3)}")

# 查看多目标规划模型的目标函数值
for i in range(MODEL.NumObj):
  MODEL.setParam(gurobipy.GRB.Param.ObjNumber, i)
  print(f"Obj {i+1} = {MODEL.ObjNVal}")

三、常用的线性化方法

1. 最大值Max

2. 最小值Min

3. 绝对值Abs

4. Maxmin/Minmax目标函数

5. 带fixed cost目标函数

6. 分式目标函数

7. 逻辑或

8. 乘积式

9. 分段整数变量

四、进阶

(1) 列表、元组和字典建模

在建模过程中，经常要对带下标数据做挑选，不同下标的数据进行组合，这样面临着二个处理方法。

1. 全部循环。多维下标意味着多重循环+ if 条件
   这样的处理方法没有效率
2. 采用特殊的Gurobi 扩展对象 TupleList 和 TupleDict
   
   如果我们使用Gurobi对其进行建模，一般有下面3（或者说4）种方法：

• 按行建模：逐行进行建模。通过使用quicksum 或者创建表达式LinExpr对象结合addTerms，拼凑表达式，最后使用Model.addConstr完成添加约束的操作，从而完成建模；
• 按列建模：逐列进行建模。通过创建Column对象，使用addTerms函数拼凑好Column对象，最后使用Model.addVar函数直接完成建模。
• 按非零系数建模：类似于逐列进行建模，可以使用addTerms拼凑表达式。具体实现跟按行建模类似。
• 按矩阵方式建模：通过拼凑约束系数矩阵，然后将决策变量转化为MVar的对象，最后直接使用重载运算符@完成约束的添加，即如Model.addConstr(A @ x == b)，就可以完成建模。

列表、元组和字典建模优点

• 添加约束代码量少，只需要拼凑好系数矩阵，即可一行代码搞定约束的添加；
• 便于直接得到对偶问题的系数矩阵，从而快速得到对偶问题的具体形式。

列表、元组和字典建模缺点

• 对于约束系数稀疏的模型，会占用不必要的内存来存储那些系数为0的部分，模型较大时，有可能导致内存溢出；

• 虽然可以得到对偶问题的具体形式，但是对偶问题的公式形式却不能得出（这个严格来讲不算是这种建模方式的缺点）。

• 不方便对每个决策变量赋予独一无二的变量名。
  例题
  

  import gurobipy as grb
  import numpy as np
  
  if __name__ == "__main__":
  	m = grb.Model("LP")
  	m.setParam('OutputFlag', 1)
  	
  	x = m.addMVar(2, lb=0, ub=grb.GRB.INFINITY)
  	
  	c = np.array([[3, 5]])
  	A = np.array([[1, 0],
  				[0, 2],
  				[3, 2]])
  	b = np.array([4, 12, 18])
  	
  	m.addConstr(A @ x <= b)
  	
  	m.Params.timelimit = 9999999999999999999999
  	
  	m.setObjective(c @ x, grb.GRB.MAXIMIZE)
  	
  	m.update()
  	m.optimize()
  	print('x_1={}'.format(x[0].x))
  	print('x_2={}'.format(x[1].x))
  

(2) 初始解设置

求解器支持用户向模型输入一个或多个已知的全部或部分可行解（一个可行解包括了模型所有变量的取值）。例如，当前问题的可行方案，以及用户通过其他算法获得的可行解等，都可以直接输入给求解器以期加速或帮助求解器进行求解。

• 利用决策变量的Start属性直接设置变量的初始值
  

  # 将变量x_1_12的初值设置为0
  x[0,12].start = 1.0
  
  # 循环为x_ij赋值完整的一组初始解,data为存放初始解的矩阵
  for i in range(len(data)):
  	for j in range(len(data[i])):
  		x[i,j].start = data[i][j]
  

	# 为模型m2设置要输入的初始解个数为2个
	m2.NumStart = 2
	# 将该可行解设置为第一个输入的可行解
	m2.Params.StartNumber = 0
	# 将该可行解设置为第二个输入的可行解
	m2.Params.StartNumber = 1

详细参考：https://mp.weixin.qq.com/s/R76jmojKGqm6ZC1dy7hSQQ

(3) Solution Pool

3.1 什么是Solution Pool

在默认状态下，Gurobi的目标是为具有单一目标函数的MIP模型寻找一个最优解。然而在Gurobi的branch and cut算法迭代的过程中，可能会通过启发式算法(也就是求解日志行首标注为H的行)得到整数可行解，也有可能通过分支切割得到整数可行解(也就是求解日志行首标注*的行)。所以在求解MIP的过程中，求解器往往会得到多个次优的可行解(Sub-Optimal Solutions)。

Solution Pool(以下简称Pool)顾名思义就是将迭代过程中产生的整数可行解存储下来的集合。(solution pool存储的只是一部分整数可行解，并不是MIP的所有可行解)。许多情况下这些可行解能够为我们提供有价值的信息，因此为方便用户了解这些解的具体情况，Gurobi提供了Solution Pool一系列相关的操作方法。

3.2 Solution Pool参数与属性

Solution Pool相关Parameters:

具体可参考该文章：https://zhuanlan.zhihu.com/p/545591977

五、数学模型

1. 线性规划

• 例1. 最简单的线性规划模型
  

  from gurobipy import * # 在Python中调用gurobi求解包
  M_LP=Model("LP_Exam") # 创建模型
  
  # 变量声明
  OF =M_LP.addVar(lb=-GRB.INFINITY,ub=GRB.INFINITY, name="OF")
  x1 =M_LP.addVar(lb=-GRB.INFINITY,ub=GRB.INFINITY, name="x1")
  x2 =M_LP.addVar(lb=-GRB.INFINITY,ub=GRB.INFINITY, name="x2")
  x3 =M_LP.addVar(lb=-GRB.INFINITY,ub=GRB.INFINITY, name="x3")
  
  # 设置目标函数
  M_LP.setObjective(x1+3*x2+3*x3,GRB.MINIMIZE)
  
  # 添加约束
  M_LP.addConstr(x1+2*x2>=3,"Con1")
  M_LP.addConstr(x2+x3>=5,"Con2")
  M_LP.addConstr(x1+x3==4,"Con3")
  
  # Optimize model
  M_LP.optimize()
  # 输出名为‘LP_Expression’的 .lp文件
  M_LP.write("LP_Expression.lp")
  # 输出结果
  print('**************')
  print(' The optimal solution ')
  print('**************')
  print('OF is :',M_LP.ObjVal) # 输出目标值
  print('x1 is :',x1.x) # 输出 X1 的值
  print('x2 is :',x2.x)
  print('x3 is :',x3.x)
  

• 例2. 具有边界的变量

# 方法1：在变量声明时进行设定
		x1=M_LP.addVar(lb=0,ub=5, name="x1")
		# 方法2：将变量边界以约束的形式加入模型中
		M_LP.addConstr(x1>=0,"Bound_Con11")
		M_LP.addConstr(x1<=5,"Bound_Con12")

from gurobipy import * # 在Python中调用gurobi求解包

	Max_Obj=1  # 用于判定目标函数求最大值还是求最小值

	M_LP=Model("LP_Exam2")

	# 变量声明
	OF =M_LP.addVar(lb=-GRB.INFINITY,ub=GRB.INFINITY, name="OF")
	x1 =M_LP.addVar(lb=0,ub=5, name="x1")
	x2 =M_LP.addVar(lb=0,ub=3, name="x2")
	x3 =M_LP.addVar(lb=0,ub=2, name="x3")

	# 设置目标函数
	if Max_Obj==0:
		M_LP.setObjective(x1+2*x2-3*x3,GRB.MINIMIZE)
	else:
		M_LP.setObjective(x1+2*x2-3*x3,GRB.MAXIMIZE)

	# 添加约束
	M_LP.addConstr(x1+2*x2<=3,"Con1")
	M_LP.addConstr(x2+x3<=2,"Con2")
	M_LP.addConstr(x1+x2+x3==4,"Con3")

	# 对变量边界进行限定的第二种写法：以x1为例
	# M_LP.addConstr(x1>=0,"Bound_Con11")
	# M_LP.addConstr(x1<=5,"Bound_Con12")

	# Optimize model
	M_LP.optimize()

	M_LP.write("LP_Expression2.lp")

	print('**************')
	print(' The optimal solution ')
	print('**************')
	print('OF is :',M_LP.ObjVal) # 输出目标值
	print('x1 is :',x1.x) # .x 用于输出 X1 的值
	print('x2 is :',x2.x)
	print('x3 is :',x3.x)

2. 混合整数线性规划

• 例1.

  	M_LP=Model("LP_Exam2")
  
  	N_i=4
  	N_j=4
  	# 设置变量
  	x =M_MIP.addVars(N_i,N_j,vtype=GRB.BINARY, name="x")
  
  	# 设置约束
  	M_MIP.addConstrs((sum(x[i,j] for j in range(N_j))<=1 for i in range(N_i)),"Con1")
  	M_MIP.addConstrs((sum(x[i,j] for i in range(N_i))<=1 for j in range(N_j)),"Con2")
  
  	# 汇总解
  	x_c=np.zeros((N_i,N_j))
  	for i in range(N_i):
  		for j in range(N_j):
  			x_c[i,j]=x[i,j].x
  

  3. 非线性规划

  在新的Gurobi求解器（9.0以上版本）中，包含了非线性/非凸问题的求解。注意：

  1. 需要设置参数’NonConvex=2’；
  2. Gurobi接受的幂函数中，次数最大值为2.
  3. 虽然有解，但不一定是最优值。
  4. 求解时间看模型规模。
     例1. 多个线性变量相乘
     
     解决办法：添加中间变量

     from gurobipy import *
     
     # 创建模型
     M_NLP=Model("NLP")
     
     # 变量声明
     x1  =M_NLP.addVar(lb=1,ub=3, name="x1")
     x2  =M_NLP.addVar(lb=1,ub=3, name="x2")
     x3  =M_NLP.addVar(lb=1,ub=3, name="x3")
     x4  =M_NLP.addVar(lb=1,ub=3, name="x4")
     
     y14 =M_NLP.addVar(lb=0,ub=300, name="y14")
     y23 =M_NLP.addVar(lb=0,ub=300, name="y23")
     
     # 设置目标函数
     M_NLP.setObjective(y14*(x1+x2+x3)+x2,GRB.MAXIMIZE)
     
     # 添加约束
     M_NLP.addConstr(y14*y23>=20,"Con1")
     M_NLP.addConstr(x1*x1+x2*x2+x3*x3+x4*x4==30,"Con2")
     
     # 表示乘积项
     M_NLP.addConstr(y14==x1*x4,"Con_y14")
     M_NLP.addConstr(y23==x2*x3,"Con_y23")
     
     M_NLP.Params.NonConvex=2  ###求解重点
     
     # Optimize model
     M_NLP.optimize()
     
     M_NLP.write("NLP.lp")
     
     print('**************')
     print(' The optimal solution ')
     print('**************')
     print('Obj is :',M_NLP.ObjVal) # 输出目标值
     print('x1 is :',x1.x) # 输出 x1 的值
     print('x2 is :',x2.x) # 输出 x2 的值
     print('x3 is :',x3.x) # 输出 x3 的值
     print('x4 is :',x4.x) # 输出 x4 的值