当前位置：物联沃-IOTWORD物联网 > 技术教程 > Gurobi教程-从入门到入土-一篇顶万篇
代码收藏家代码收藏家 技术教程

Gurobi教程-从入门到入土-一篇顶万篇

再详尽的帮助文档也不如举几个示例能让人看得明白，于是想通过一个帖子涵盖Gurobi的所有操作。

安装与激活

软件下载、免学术ip申请学术许可，参见 Gurobi中国官方网站。具体过程简单且网上很多教程。值得注意的是当运行pycharm时发现Gurobi的license过期，在重新安装 gurobi.lic 时若选择的地址不是系统环境变量，需要在系统环境变量中添加该.lic，变量名为“GRB_LICENSE_FILE”，并重启电脑！

版权

以下引自官网：

GUROBI 为中国学校教师、学生提供半年免费使用版本（可延续），功能没有限制，需要符合以下规定。同时激活时的机器IP地址将会被记录，如果存在违反以下规定的情况，将有可能面临不利后果。
（1）在中国高等教育机构的教师或者在籍学生（研究所不在此列。中科院系统需要具备中科院大学证件）；
（2）用于非商业用途；
（3）学术许可只能申请本人使用，为其他人申请将会造成申请人和使用人同时失去使用学术许可的永久资格；
（4）学术许可严禁安装和使用在非学术科研的场合和设备上；
（5）学生和教师严禁在企业和不符申请资格的科研机构上安装和使用学术许可；
（6）企业和科研机构违规安装和使用学术许可，将会被追究经济和法律责任；
（7）使用Gurobi 的科研论文的第一作者需要具有 Gurobi 正版许可（合规学术许可或者正版商业许可）；
（8）遵守其他版权规定 http://www.gurobi.com/pdf/eula/eula.pdf

其中第（7）条需要格外注意。

MIP模型

举个例子

例如求解VRP问题。

import numpy as np
import matplotlib.pyplot as plt
from gurobipy import *

rnd = np.random
rnd.seed(1) # 随机种子，如有本行，则程序每次运行结果一样。可任意赋值

n = 10 # 一共几个客户点/城市/需求点
xc = rnd.rand(n+1)*200 # 随机生成每个城市的横坐标，范围[0,200]
yc = rnd.rand(n+1)*100 # 随机生成每个城市的纵坐标，范围[0,100]

# 可以画图看一眼生成的城市什么样子
# plt.plot(xc[0], yc[0], c='r',marker='s' ) # 索引为0的点，即depot/仓库/出发点
# plt.scatter(xc, yc, c='b') # 客户点

N = list(range(1,n+1)) # 客户点集合
V = list(range(0,n+1)) # 所有点集合（仓库+客户点）
A = [(i,j) for i in V for j in V if i !=j] # 城市之间有哪些连线/路段/弧段
C = {(i,j): np.hypot(xc[i]-xc[j], yc[i]-yc[j]) for i,j in A} # np.hypot：二范数=求平方和；计算弧段的长度
Q = 20 # 车最大载重
q = {i:rnd.randint(1,10) for i in N} # 随机生成客户点的需求量，范围[1,10]

mdl = Model('CVRP') # 起名字

x = mdl.addVars(A, vtype=GRB.BINARY) # 增加变量xij，表示是否链接ij
u = mdl.addVars(N, vtype=GRB.CONTINUOUS) # 增加变量ui，表示车在该点处累计载货量

mdl.modelSense = GRB.MINIMIZE # 目标为最小化
mdl.setObjective(quicksum(x[i,j]*C[i,j] for i,j in A )) # 目标函数为总距离

# 添加所有约束
mdl.addConstrs(quicksum(x[i,j] for j in V if i != j)==1 for i in N) 
mdl.addConstrs(quicksum(x[i,j] for i in V if i!=j)==1 for j in N)
mdl.addConstrs((x[i,j]==1) >> 
               (u[i] + q[j] == u[j]) for i,j in A if i!=0 and j!=0)
mdl.addConstrs(u[i] >= q[i] for i in N)
mdl.addConstrs(u[i] <= Q for i in N)

mdl.optimize() # 优化

#优化完成，下面输出结果
active_arts = [a for a in A if x[a].x > 0.9] # 输出最优解的所有连线，即xij中是1的(i,j)
# 由于存在误差，xij可能为0.999999999，因此不要用==1
print(active_arts)

# 画图
for index, (i,j) in enumerate(active_arts):
    plt.plot([xc[i],xc[j]],[yc[i],yc[j]],c='r')
plt.plot(xc[0], yc[0], c='g',marker='s' )
plt.scatter(xc, yc, c='b')

输出的啥

求解某VRP模型后，输出结果如下：

第一部分

可以不看。

Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 40 rows, 120 columns and 220 nonzeros
Model fingerprint: 0x56854507

第二部分

总结了模型的规模。包括有几个约束条件，几个变量。Coefficient statistics暂时不知道什么意思。

Model has 90 general constraints
Variable types: 10 continuous, 110 integer (110 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+01, 2e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 2e+01]

第三部分

Gurobi求解MIP的presolve阶段。

Presolve added 171 rows and 8 columns
Presolve time: 0.01s
Presolved: 211 rows, 128 columns, 1318 nonzeros
Variable types: 38 continuous, 90 integer (90 binary)

第四部分

展示求解过程，其中1：

Root relaxation: objective 3.612672e+02, 40 iterations, 0.00 seconds
    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  361.26720    0   21          -  361.26720      -     -    0s
H    0     0                    1130.6043611  361.26720  68.0%     -    0s
H    0     0                    1108.0315407  361.26720  67.4%     -    0s
H    0     0                     899.3261060  361.26720  59.8%     -    0s
H    0     0                     773.2318022  361.26720  53.3%     -    0s
H    0     0                     748.6243052  361.26720  51.7%     -    0s
H    0     0                     696.3684124  361.26720  48.1%     -    0s
H    0     0                     689.4096665  366.86305  46.8%     -    0s
     0     0  370.85380    0   22  689.40967  370.85380  46.2%     -    0s
H    0     0                     663.9889500  370.85380  44.1%     -    0s
     0     0  375.47206    0   22  663.98895  375.47206  43.5%     -    0s
     0     0  375.47206    0   22  663.98895  375.47206  43.5%     -    0s
     0     0  378.67819    0   22  663.98895  378.67819  43.0%     -    0s
     0     0  378.67819    0   22  663.98895  378.67819  43.0%     -    0s
     0     2  378.67819    0   21  663.98895  378.67819  43.0%     -    0s
H   31    40                     611.6880568  395.45082  35.4%  20.0    0s

Nodes（1-2列）：Gurobi在利用 branch-and-cut（分支剪界法）求解MIP时的过程。第1列：已搜索过的节点个数。第2列：还没有被搜索的叶子节点。第1列会越来越多，第2列数字可能偶尔减小。已搜索节点数（第1列）总是为0的话，表示Gurobi MIP solver正在处理根节点。在根节点需要花时间产生cutting plane，以及尝试多种heuristic方法来减小后面branch-and-cut tree的规模。当某行前面标有 H 或 *，表示找到了一个新的可行解。H 和 *分别表示该可行解是通过 heuristic 或 branching 方法找到的。

Current Node（3-5列）：branch-and-cut tree中被搜索过的特定节点的信息。第3列：相应relaxation（松弛）的目标值。第4列：branch-and-cut tree中节点的深度。第5列：非整数值的整数变量的个数。

Objective Bounds（6-8列）：关于可行解的区间信息。对最小化问题来说，第6列：上界（目前找到的最优解）。第7列：下界（线性松弛后的最优解，由于某些把某些变量松弛为连续值，因此一定是不可行解，但全局的可行的最优解不会比该松弛解更优）。第8列：Gap（(上界-下界)/上界）。

Work（9-10列）：求解运算的工作量。第9列：branch-and-cut tree中每个节点simplex迭代的平均次数，Number of iterations/node。第10列：模型求解时间。

第五部分

Solution count：共找到几个可行解，分别是什么

Cutting planes:
  Learned: 4
  Gomory: 7
  Cover: 17
  Implied bound: 18
  MIR: 5
  GUB cover: 2
  Zero half: 2
  Relax-and-lift: 4
  
Explored 2398 nodes (36083 simplex iterations) in 0.23 seconds
Thread count was 16 (of 16 available processors)

第六部分

最重要最关心的部分。
Solution count：共找到几个可行解，分别是什么

Solution count 9: 611.688 663.989 689.41 ... 1130.6
Optimal solution found (tolerance 1.00e-04)
Best objective 6.116880568072e+02, best bound 6.116880568072e+02, gap 0.0000%

隐藏的坑

python数据存储误差

有时为了能够将模型线性化，需要借助大M法。通常将M设为一个很大量级的数，如1e+6等。但在有些情况下，M过大将导致一些错误。python存储数据时可能存在误差，即使是0-1变量，优化结束后可能并不是整数0或1，而是3.3306690738754696e-6或0.99999879454，此时乘以一个较大的M时，量级将回到一个不可忽略的范围，导致约束条件判断错误。
如在路径规划问题中消除子环时，M本应可以是越大越好的正数，但实际操作时仍推荐找到M的比较紧的上界2。同理，代码中的M0也是如此。

mdl.addConstrs(miu[d, i, j, m] - miu[d, i, j, n] + M * y[d,i,j,m,n] 
			   <= M - 1 + M0 * x[j,n]
			   for d in D for i in y_i_set for j in y_j_set if j > i for m in V for n in V)

同时这个特性导致我们想查找结果为1的变量时，最好不要写x[i]==1，而是x[i]>0.9。

额外技巧

设置%gap或求解时间限制

面对大规模问题，不指望精确算法能够给出全局最优解，因此设置参数以控制运算时间。

  • %gap=|BP−BF|/|BP|， 其中 BP 为整数连续化后得到的最优解，BF 则为当前求解器求得的整数时的解。就是输出结果中 Objective Bounds 的 Gap 值。3
  • TimeLimit，求解时限。

    • 设置方法4：

    # 在调用 model.optimize() 之前进行如下设置。 MIPGap 和 TimeLimit 的单位分别是分数和秒。 
model.Params.MIPGap = 0.05    # %gap = 5%; 或 model.setParam('MIPGap', 0.05)
model.Params.TimeLimit = 300  # TimeLimit = 5 minutes; 或 model.setParam('Timelimit', 300)

    1. https://blog.csdn.net/weixin_33467739/article/details/112110951 ↩︎

    2. https://www.zhihu.com/question/42175410 ↩︎

    3. https://blog.csdn.net/qinghai5068/article/details/51443851 ↩︎

    4. https://stackoom.com/question/4LZjl ↩︎

    您可能感兴趣的内容:

    物联沃分享整理
    物联沃-IOTWORD物联网 » Gurobi教程-从入门到入土-一篇顶万篇

    代码收藏家 普通

    分享到：

    相关推荐

    发表评论