Gurobi 求解线性规划问题¶
本文针对钢铁企业对供应商的选择问题，构建线性规划模型，并使用 Gurobi 求解器进行求解。
题目背景¶
考虑一家小型的钢铁企业，该企业炼钢时使用的主要原材料是炼焦煤，每年的需求量在 100 到 150 万吨。现在需要帮助该公司安排明年的生产，选择原料的供应商。目前他们收到了 8 家供应商的报价，如下表。表格中的信息包括了每家供应商的最大供应量、是否为工会的公司、运输的方式、炼焦煤的可燃率、单位价格。
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |
|---|---|---|---|---|---|---|---|---|
| 供应量 （千吨/年） | 300 | 600 | 510 | 655 | 575 | 680 | 450 | 490 |
| 工会U/ 非工会 N | U | U | N | U | N | U | N | N |
| 卡车T/ 铁路 R | R | T | R | T | T | T | R | R |
| 可燃率（%） | 15 | 16 | 18 | 20 | 21 | 22 | 23 | 25 |
| 价格 (￥/吨) | 49.5 | 50.0 | 61.0 | 63.5 | 66.5 | 71.0 | 72.5 | 80.0 |
现在该公司有如下考量：
-
根据对未来的形势估计和生产现状的考察，该公司计划明年要订购 122.5 万吨的炼焦煤。
-
要求这些炼焦煤的平均可燃率至少要达到 19%。
-
为了不影响劳工关系，该公司决定至少 50% 的炼焦煤要来自于工会煤矿。
-
每年由火车运输的炼焦煤量至多不超过 650 千吨，而由卡车运输的炼焦煤量至多不超过 720 千吨。
现在要使该公司原材料的成本最小，应该从各个供应商处订购多少吨炼焦煤？
-
建立模型并求解（软件求解），帮助该公司制定采购方案。
-
该公司原材料的边际成本是多少？（即多订购一吨炼焦煤的成本）
-
该公司会考虑扩展他们的卡车或者火车运输吗？如果会的话，他们会愿意为扩展相应的运输业务支付多少呢？
-
该公司现有的对于工会或者非工会煤矿的订购限制要求是否增加了原材料的订购成本？是否需要修订该政策呢？
-
目前的情况来看，从供应商 6 处订购是否经济呢？如果不经济，那么供应商 6 需要如何降低自身的定价，从而使得该公司会向他们订购？
-
该公司会为从供应商 2 处多订购炼焦煤而提高定价吗？如果会的话，那么该公司会愿意以更高的价格增加多少的订购量？同样的，对于供应商 7 呢？
建模与分析¶
使用 Gurobi 求解¶
导入包¶
Python
import numpy as np
import gurobipy as gp
from gurobipy import GRB
np.set_printoptions(linewidth=np.inf)
设定线性规划问题的参数¶
Python
c = gp.tuplelist(
[
-49.5,
-50,
-61,
-63.5,
-66.5,
-71,
-72.5,
-80,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
]
)
b = gp.tuplelist(
[
1225000,
1225000 * 0.19,
1225000 * 0.5,
650000,
720000,
300000,
600000,
510000,
655000,
575000,
680000,
450000,
490000,
]
)
A = gp.tuplelist(
[
[1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[
0.15,
0.16,
0.18,
0.2,
0.21,
0.22,
0.23,
0.25,
-1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
],
[1, 1, 0, 1, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
],
)
建立 Gurobi 模型¶
Python
m = gp.Model("Exercise")
x = m.addVars(
range(20),
vtype=GRB.CONTINUOUS,
name="x",
)
cons = m.addConstrs(
(x.prod(A[i]) == b[i] for i in range(len(b))), name="constrs" # <= >=
)
m.setObjective(x.prod(c), GRB.MAXIMIZE)
求最优解¶
Text Only
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (mac64[arm])
CPU model: Apple M1 Pro
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 13 rows, 20 columns and 48 nonzeros
Model fingerprint: 0x22742dfc
Coefficient statistics:
Matrix range [1e-01, 1e+00]
Objective range [5e+01, 8e+01]
Bounds range [0e+00, 0e+00]
RHS range [2e+05, 1e+06]
Presolve removed 8 rows and 12 columns
Presolve time: 0.00s
Presolved: 5 rows, 8 columns, 28 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 -6.0637500e+07 2.120000e+05 0.000000e+00 0s
5 -7.3267500e+07 0.000000e+00 0.000000e+00 0s
Solved in 5 iterations and 0.01 seconds (0.00 work units)
Optimal objective -7.326750000e+07
获取原问题的最优解¶
Text Only
x[0] = 55000.000000
x[1] = 600000.000000
x[2] = 0.000000
x[3] = 20000.000000
x[4] = 100000.000000
x[5] = 0.000000
x[6] = 450000.000000
x[7] = 0.000000
x[8] = 0.000000
x[9] = 62500.000000
x[10] = 145000.000000
x[11] = 0.000000
x[12] = 245000.000000
x[13] = 0.000000
x[14] = 510000.000000
x[15] = 635000.000000
x[16] = 475000.000000
x[17] = 680000.000000
x[18] = 0.000000
x[19] = 490000.000000
获取对偶变量¶
Python
p_list = []
for constr in m.getConstrs():
print(constr.constrName, "=", constr.Pi)
p_list.append(constr.Pi)
Text Only
constrs[0] = -4.499999999999915
constrs[1] = -300.00000000000057
constrs[2] = -0.0
constrs[3] = -0.0
constrs[4] = 1.0000000000000338
constrs[5] = -0.0
constrs[6] = 1.4999999999999734
constrs[7] = -0.0
constrs[8] = -0.0
constrs[9] = -0.0
constrs[10] = -0.0
constrs[11] = 1.0000000000000426
constrs[12] = -0.0
计算其他量¶
Python
# 求 C_B 矩阵，它是 C 中基变量对应的列向量组成的矩阵，维数是 13*1
C_B = np.array([c[i] for i in range(A.shape[1]) if x[i].x != 0])
Text Only
array([-49.5, -50. , -63.5, -66.5, -72.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ])
Python
# 求 B 矩阵，它是 A 的子矩阵，也就是 A 中基变量对应的列向量组成的矩阵，维数是 13*13
B = np.array([A[:, i] for i in range(A.shape[1]) if x[i].x != 0]).T
Text Only
array([1. , 0.16, 1. , 0. , 1. , 0. , 1. , 0. , 0. , 0. , 0. , 0. , 0. ])
Text Only
-7.105427357601002e-15
Text Only
array([[ 1., 0., 0., 0., -1., 0., 0., 0., 0., 0., 0., -1., 0.],
[ -0., -0., -0., -0., -0., -0., 1., -0., -0., -0., -0., -0., -0.],
[ 15., -100., -0., -0., 6., -0., -5., -0., -0., -0., -0., 8., -0.],
[ -15., 100., -0., -0., -5., -0., 4., -0., -0., -0., -0., -8., -0.],
[ -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., 1., -0.],
[ 16., -100., -1., 0., 5., 0., -4., 0., 0., 0., 0., 7., 0.],
[ -1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
[ -1., 0., 0., 0., 1., 1., 0., 0., 0., 0., 0., 1., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[ -15., 100., 0., 0., -6., 0., 5., 0., 1., 0., 0., -8., 0.],
[ 15., -100., -0., -0., 5., -0., -4., -0., -0., 1., -0., 8., -0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]])
Text Only
[1225000, 232750.0, 612500.0, 650000, 720000, 300000, 600000, 510000, 655000, 575000, 680000, 450000, 490000]
Text Only
array([ 55000., 600000., 20000., 100000., 450000., 62500., 145000., 245000., 510000., 635000., 475000., 680000., 490000.])