线性规划

线性规划简介

简单来说线性规划（LP）是指在给出一系列约束条件后求解目标函数的极值问题，该类题型在中学数学中常见，这里不再赘述。一个简单的例子:

$$\max \ f(x)=2x_1 + 3x_2 - 5x3\\ s.t \begin{cases} x1 + 3x_2 + x_3 < 12\\ 2x1 - 5x_2 + x_3 \geq 10\\ x_1 + x_2 + x_3 = 7\\ x_1,x_2,x_3\leq 0 \end{cases}$$

$x_1$、$x_2$、$x_3$称为决策变量，$f(x)$称为目标函数，$s.t.$称为约束条件，这是线性规划问题的三个关键因素。在数学建模中，我们需要在具体问题中抽象出这三个关键因素，并分析它们的属性，步骤如下：

• 找到决策变量。判断它们是连续的还是离散的。
• 写出目标函数。这通常是问题所在，如最大利润，最小成本。目标函数的自变量为决策变量。
• 分析约束条件。约束条件指的是对单个决策变量的约束和对多个决策变量的关系进行约束。

PuLP库求解过程

python中的Pulp求解以上问题。以上述问题为例，详细介绍一下利用pulp库求解线性规划问题的详细步骤：

定义规划问题

LpProb = pulp.LpProblem(name = 'LpProb', sense = pulp.LpMaximize)

pulp.LpProblem()有两个参数：

• name：该问题的名字。
• sense：指定待解决的问题的类型。pulp.LpMaximise为求目标函数的最大值，pulp.LpMinimize为求目标函数的最小值。

定义决策变量

x1 = pulp.LpVariable('x1', lowBound = 0, upBound = 7, cat = 'Continuous')
x2 = pulp.LpVariable('x2', lowBound = 0, upBound = 7, cat = 'Continuous')
x3 = pulp.LpVariable('x3', lowBound = 0, upBound = 7, cat = 'Continuous')

利用pulp.LpVariable()函数创建问题的决策变量。它有五个参数：

• name：变量名，在.lp输出文件中将被用到
• lowBound：变量的下限，默认为none，代表负无穷
• upBound：变量的上限，默认为none，代表正无穷
• cat：设定决策变量的类型，默认为'Continuous'，表示变量是连续的。另外还有'Interger'表示变量是离散的，此时该变量的取值只能为整数；'Binary'表示这是一个二值变量，变量只能为$0/1$
• e：指明该决策变量是否在目标函数和约束条件中存在，这里暂时不需要设置该参数

因此，本题我们定义了三个决策变量$x_1$、$x_2$、$x_3$，它们下限均为$0$，上限均为$7$，且均为连续型变量。

定义目标函数

目标函数就是f(x)，添加目标函数的格式为：问题变量 += 目标函数表达式

LpProb += 2*x1 + 3*x2 - 5*x3

定义约束条件

定义约束条件的格式与定义目标函数的格式类似，但是约束条件式均为等式或不等式：问题变量 += 约束条件式

LpProb += 2*x1 - 5*x2 + x3 >= 10
LpProb += x1 + 3*x2 + x3 <= 12
LpProb += x1 + x2 + x3 == 7

求解

LpProb.solve()

solve()参数可以指定PuLP进行求解时使用的求解器，默认为CBC。也可以用其它求解器，如solve(CPLEX())

输出结果

仅仅调用solve()函数求解还不够，我们还要知道该问题是否有解，最后的结果是多少

print("status:", pulp.LpStatus[LpProb.status])
for v in LpProb.variables():
    print(v.name, "=", v.varValue)
print('F(x) = ', pulp.value(LpProb.objective))

首先请求解决方案的状态，可以是“Not Solved（未解决）”, “Infeasible（不可行）”, “Unbounded（无界）”, “Undefined（未定义）” 和 “Optimal（最优解）”。LpStatus是一个字典类型。

其次需要将每个决策变量的值都打印出来，variables()将返回一个列表，遍历打印每一个变量名与变量值。

最后打印目标函数的最优值。

算法步骤

综上，利用PuLP求解一般的线性规划问题步骤总结如下：

1. 使用pulp.LpProblem()创建问题实例
2. 使用pulp.LpVariable()定义决策变量
3. LpProb += 表达式设置目标函数
4. LpProb += 关系式设置约束条件
5. 调用solve()求解
6. 输出求解状态、变量取值、目标函数值

上述问题完整代码如下：

# 导入import库
import pulp
# 定义一个规划问题
LpProb = pulp.LpProblem("LPProbDemo1", sense = pulp.LpMaximize)
# 定义决策变量
x1 = pulp.LpVariable('x1', lowBound = 0, upBound = 7, cat = 'Continuous')
x2 = pulp.LpVariable('x2', lowBound = 0, upBound = 7, cat = 'Continuous')
x3 = pulp.LpVariable('x3', lowBound = 0, upBound = 7, cat = 'Continuous')
# 添加目标函数
LpProb += 2*x1 + 3*x2 - 5*x3
# 添加约束函数
LpProb += 2*x1 - 5*x2 + x3 >= 10
LpProb += x1 + 3*x2 + x3 <= 12
LpProb += x1 + x2 + x3 == 7
# 求解
LpProb.solve()
print("status:", pulp.LpStatus[LpProb.status])
for v in LpProb.variables():
    print(v.name, "=", v.varValue)
print('F(x) = ', pulp.value(LpProb.objective))

结果如下：

status: Optimal
x1 = 6.4285714
x2 = 0.57142857
x3 = 0.0
F(x) =  14.57142851

应用案例

以上介绍了PuLP求解线性规划问题的基本方法，下面介绍几个应用案例。

混合问题

问题描述

Uncle Ben’s希望尽可能便宜地生产他们的猫粮产品，同时确保它们符合罐头上显示的营养分析要求。因此，他们希望改变每种成分（主要成分是鸡肉，牛肉，羊肉，大米，小麦和凝胶）的数量，同时仍然符合他们的营养标准。鸡肉，牛肉和羊肉的成本分别为0.013美元，0.008美元和0.010美元，而大米，小麦和凝胶的成本分别为0.002美元，0.005美元和0.001美元。（所有费用均为每克。）

每种成分将组成最终产品中蛋白质，脂肪，纤维和盐的总重量。每克成分的贡献（以克为单位）如下表所示。

东西	蛋白	脂肪	纤维	盐
鸡	0.100	0.080	0.001	0.002
牛肉	0.200	0.100	0.005	0.005
羊肉	0.150	0.110	0.003	0.007
米饭	0.000	0.010	0.100	0.002
麦麸	0.040	0.010	0.150	0.008
凝胶	0.000	0.000	0.000	0.000

为了满足营养分析要求，我们需要每100克至少含有8克蛋白质，6克脂肪，但不超过2克纤维和0.4克盐，请问满足营养要求的最小成本开销是多少？

解题过程

首先要从问题中分析出决策变量、目标函数和约束条件。

设鸡肉、牛肉、羊肉、米饭、麦麸和凝胶的质量分别为$x_1$、$x_2$、$x_3$、$x_4$、$x_5$、$x_6$克，那么目标函数和约束条件如下：

$$\min \ f(x) = 0.013x_1 + 0.008x_2 + 0.010x_3 + 0.002x_4 + 0.005x_5 + 0.001x_6\\ s.t. \begin{cases} x_1 + x_2 + x_3 + x_4 + x_5 + x_6=100\\ 0.100x_1 + 0.200x_2 + 0.150x_3 + 0.040x_5 \geq 8\\ 0.080x_1 + 0.100x_2 + 0.110x_3 + 0.010x_4 + 0.010x_5 \geq 6\\ 0.001x_1 + 0.005x_2 + 0.003x_3 + 0.100x_4 + 0.150x_5 \leq 2\\ 0.002x_1 + 0.005x_2 + 0.007x_3 + 0.002x_4 + 0.008x_5 \leq 0.4\\ \end{cases}$$

from pulp import *
# 定义问题
blendLp = LpProblem('混合问题', LpMinimize)
# 定义约束变量
x1 = LpVariable('x1', 0, None, LpContinuous)
x2 = LpVariable('x2', 0, None, LpContinuous)
x3 = LpVariable('x3', 0, None, LpContinuous)
x4 = LpVariable('x4', 0, None, LpContinuous)
x5 = LpVariable('x5', 0, None, LpContinuous)
x6 = LpVariable('x6', 0, None, LpContinuous)
# 设置目标函数
blendLp += 0.013*x1 + 0.008*x2 + 0.01*x3 + 0.002*x4 + 0.005*x5 + 0.001*x6
# 设置约束条件
blendLp += x1 + x2 + x3 + x4 + x5 + x6 == 100
blendLp += 0.100*x1 + 0.200*x2 + 0.150*x3 + 0.040*x5 >= 8
blendLp += 0.080*x1 + 0.100*x2 + 0.110*x3 + 0.010*x4 + 0.010*x5 >= 6
blendLp += 0.001*x1 + 0.005*x2 + 0.003*x3 + 0.100*x4 + 0.150*x5 <= 2
blendLp += 0.002*x1 + 0.005*x2 + 0.007*x3 + 0.002*x4 + 0.008*x5 <=0.4
# 求解
blendLp.solve()
# 打印求解结果
print('status: ', LpStatus[blendLp.status])
for v in blendLp.variables():
    print(v.name, "=", v.varValue)
print('min = ', value(blendLp.objective))

当决策变量与约束条件较多时，使用字典更规范，也更易于修改

Ingerdient = ['chicken', 'beef', 'mutton', 'rice', 'wheat', 'gel']
    cost = {
    'chicken' : 0.013,
    'beef' : 0.008,
    'mutton' : 0.01,
    'rice' : 0.002,
    'wheat' : 0.04,
    'gel' : 0.001
    }
    protein = {
    'chicken' : 0.1,
    'beef' : 0.2,
    'mutton' : 0.15,
    'rice' : 0,
    'wheat' : 0.04,
    'gel' : 0
    }
    fat = {
    'chicken' : 0.08,
    'beef' : 0.1,
    'mutton' : 0.11,
    'rice' : 0.01,
    'wheat' : 0.01,
    'gel' : 0
    }
    fibre = {
    'chicken' : 0.001,
    'beef' : 0.005,
    'mutton' : 0.003,
    'rice' : 0.1,
    'wheat' : 0.15,
    'gel' : 0,
    }
    salt = {
    'chicken' : 0.002,
    'beef' : 0.005,
    'mutton' : 0.007,
    'rice' : 0.002,
    'wheat' : 0.008,
    'gel' : 0
    }
    blendLp = LpProblem("混合问题", LpMinimize)
    varsLp = LpVariable.dict('Ingerdient', Ingerdient, 0)
    blendLp += (lpSum([cost[i] * varsLp[i] for i in Ingerdient]), 
        "total cost of Ingerdient per can")
    blendLp += (lpSum([varsLp[i] for i in Ingerdient]) == 100, "PercentagesSum")
    blendLp += (lpSum([varsLp[i] * protein[i] for i in Ingerdient]) >= 8, "ProteinRequirement")
    blendLp += (lpSum([varsLp[i] * fat[i] for i in Ingerdient]) >= 6, "FatRequirement")
    blendLp += (lpSum([varsLp[i] * fibre[i] for i in Ingerdient]) <= 2, "FibreRequirement")
    blendLp += (lpSum([varsLp[i] * salt[i] for i in Ingerdient]) <=4, "saltRequirement" )
    blendLp.solve()

0-1规划

问题描述

公司有 5 个项目被列入投资计划，各项目的投资额和预期投资收益如下表所示（万元）：

项目	A	B	C	D	E
投资额	210	300	100	130	260
投资收益	150	210	60	80	180

公司只有 600万元资金可用于投资，综合考虑各方面因素，需要保证：

（1）项目 A、B、C 中必须且只能有一项被选中；
（2）项目 C、D 中最多只能选中一项；
（3）选择项目 E 的前提是项目 A被选中。

如何在上述条件下，进行投资决策，使收益最大。

解题过程

定义决策变量：定义五个二值变量$x_1$、$x_2$、$x_3$、$x_4$、$x_5$，表示每个项目是否被选中，被选中则为1，未被选中则为0.

目标函数和约束条件如下：

$$\max \ F=150x_1 + 210x_2 + 60x_3 + 80x_4 + 180x_5\\ s.t. \begin{cases} 210x_1 + 300x_2 + 100x_3 + 130x_4 + 260x_5 \leq600\\ x_1 + x_2 + x_3 = 1\\ x_3 + x_4 \leq1\\ x_5 \leq x_1 \end{cases}$$

代码如下

prob = LpProblem("0-1规划", LpMaximize)
x1 = LpVariable('A', cat = 'Binary')
x2 = LpVariable('B', cat = 'Binary')
x3 = LpVariable('C', cat = 'Binary')
x4 = LpVariable('D', cat = 'Binary')
x5 = LpVariable('E', cat = 'Binary')
prob += 150*x1 + 210*x2 + 60*x3 + 80*x4 + 180*x5
prob += 210*x1 + 300*x2 + 100*x3 + 130*x4 + 260*x5 <= 600
prob += x1 + x2 + x3 == 1
prob += x3 + x4 <= 1
prob += x1 - x5 <= 0
prob.solve()
print("status: ", LpStatus[prob.status])
for v in prob.variables():
    print(v.name, "=", v.varValue)
print("max = ", value(prob.objective))

得到如下结果

status:  Optimal
A = 1.0
B = 0.0
C = 0.0
D = 1.0
E = 1.0
max =  410.0

数独问题

问题描述

数独的规则为：

• 在9个单独的3x3盒子中的任何一个中，必须找到每个数字1到9。
• 在 9x9 网格的任何列中，必须找到每个数字 1 到 9。
• 在 9x9 网格的任何一行中，必须找到每个数字 1 到 9。

解题过程

问题的关键在于寻找决策变量。

最容易想到的方法是每一个单元设置一个变量，代表要填入的数，这样需要设置81个决策变量。但是这样无法设置约束条件，即同一行同一列和同一$3x3$的盒子中不能有相同值。

这里，我们创建 729 个单独的二进制 （0-1） 决策变量，表示9x9网格中每个单元的9种取值情况。每种取值有两种情况，即是（1）和不是（0）。我们用一个三元组表示决策变量，格式为$(vals, rows, cols)$，即第$rows$行第$cols$列值是否为$val$。例如$(5, 7, 3)$表示第7行第3列是否为5，如果是，该变量取值为$1$，否则取值为$0$。

另外注意，该题不需要设置目标函数，因为它并不是求最优解问题，只要满足约束条件，即为一个可行解。

理解了决策变量的含意后，就很容易设置约束条件了。

• 同一个单元中只能有一个值
• 同一个值在一行中有且只有一个
• 同一个值在一列中有且只有一个

from pulp import *
# 总共有9行9列，每一个单元中可以选择的值有9个
vals = rows = cols = range(1, 10)
#定义一个问题，不需要设置目标函数，也不需要有第二个参数
prob = LpProblem("数独问题")
# 设置决策变量
choices = LpVariable.dicts("choice", (vals, rows, cols), cat = 'Binary')
# 设置约束条件
# 同一个单元只有一个值
for r in rows:
    for c in cols:
        prob += lpSum([choices[v][r][c] for v in vals]) == 1
# 同一个值在一行中有且只有一个
for v in vals:
    for r in rows:
        prob += lpSum([choices[v][r][c] for c in cols]) == 1
# 同一个值在一列中有些只有一个
for v in vals:
    for c in cols:
        prob += lpSum([choices[v][r][c] for r in rows]) == 1
# 将已知的值填入
input_data = [
    (5, 1, 1),
    (6, 2, 1),
    (8, 4, 1),
    (4, 5, 1),
    (7, 6, 1),
    (3, 1, 2),
    (9, 3, 2),
    (6, 7, 2),
    (8, 3, 3),
    (1, 2, 4),
    (8, 5, 4),
    (4, 8, 4),
    (7, 1, 5),
    (9, 2, 5),
    (6, 4, 5),
    (2, 6, 5),
    (1, 8, 5),
    (8, 9, 5),
    (5, 2, 6),
    (3, 5, 6),
    (9, 8, 6),
    (2, 7, 7),
    (6, 3, 8),
    (8, 7, 8),
    (7, 9, 8),
    (3, 4, 9),
    (1, 5, 9),
    (6, 6, 9),
    (5, 8, 9),
]
for (v, r, c) in input_data:
    prob += choices[v][r][c] == 1
# 解决问题
prob.solve()
# 打印结果
print("status: ", LpStatus[prob.status])

for r in rows:
    if r in [1, 4, 7]:
        print("+-------+-------+-------+")
    for c in cols:
        for v in vals:
            if(value(choices[v][r][c])) == 1:
                if c in [1, 4, 7]:
                    print("| ", end = "")
                print(str(v) + " ", end="")
                if c == 9:
                    print("|")
print("+-------+-------+-------+")

结果如下，答案不止一种。

status:  Optimal
+-------+-------+-------+
| 5 3 9 | 6 7 2 | 8 1 4 |
| 6 7 3 | 1 9 5 | 4 2 8 |
| 3 9 8 | 7 4 1 | 5 6 2 |
+-------+-------+-------+
| 8 5 2 | 9 6 7 | 1 4 3 |
| 4 2 6 | 8 5 3 | 7 9 1 |
| 7 1 4 | 3 2 8 | 9 5 6 |
+-------+-------+-------+
| 9 6 1 | 5 3 4 | 2 8 7 |
| 2 8 7 | 4 1 9 | 6 3 5 |
| 1 4 5 | 2 8 6 | 3 7 9 |
+-------+-------+-------+

运输问题

问题描述

一家啤酒厂有两个仓库，从中将啤酒分发给五个酒吧。在每周开始时，每个酒吧都会向啤酒厂总部发送订单，订购很多箱啤酒，然后从相应的仓库发送到酒吧。啤酒厂希望有一个交互式计算机程序可以告诉他们哪个仓库应该供应哪个酒吧，以最大限度地降低整个运营的成本。例如，假设在给定的一周开始时，啤酒厂在仓库 A 有 1000 个箱子，在仓库 B 有 4000 个箱子，而条形图分别需要 500、900、1800、200 和 700 个箱子。哪个仓库应该供应哪个酒吧？

我们假设每个箱子都有固定的运输成本，成本如下（元/箱）

From Warehouse to Bar	A	B
1	2	3
2	4	1
3	5	3
4	2	2
5	1	3

解题过程

决策变量为每一个仓库送到每一个酒吧的箱子个数，$A_i$表示从仓库A运输到酒吧$i$的箱子个数，$B_i$表示从仓库B运输到酒吧$i$的箱子个数，$Wa_i$表示从仓库A运输到酒吧$i$的成本，$Wb_i$表示从仓库B运输到酒吧$i$的成本，$D_i$表示酒吧$i$的需求量

$$\min \ F = \sum_{i=1}^5 (A_iWa_i+B_iWb_i) \\ s.t:\begin{cases} \sum_{i = 1}^5 A_i \leq 1000 \\ \sum_{i = 1}^5 B_i \leq 4000\\ A_i + B_i \geq D_i &{i\in\{1, 2, 3, 4, 5\}}\\ \end{cases}$$

该题利用“分析决策变量，设置目标函数，设置约束条件”的步骤做，由于变量较多，利用字典来表示变量的值。

from pulp import *
warehouse = ['A', 'B']
supply = {'A': 1000, 'B': 4000}
bars = ['1', '2', '3', '4', '5']
demand = {
    '1': 500,
    '2': 900,
    '3': 1800,
    '4': 200,
    '5': 700
}
# 运输成本
cost = {
    'A': [2, 4, 5, 2, 1],
    'B': [3, 1, 3, 2, 3]
}
routes = [(w, b) for w in warehouse for b in bars]
# 定义问题
prob = LpProblem('运输问题', LpMinimize)
# 定义决策变量
vars = LpVariable.dicts("路线", (warehouse, bars), 0, None, LpInteger)
# 设置目标函数
prob += (lpSum([vars[w][b] * cost[w][int(b)-1] for (w, b) in routes]),
    '运输成本')
# 设置约束条件
for w in warehouse:
    prob += (lpSum([vars[w][b] for b in bars]) <= supply[w],
        "运出仓库%s" % w)
for b in bars:
    prob += (lpSum([vars[w][b] for w in warehouse]) >= demand[b], 
        "运进酒吧%s" % b)
# 求解问题
prob.solve()
print('status: ', LpStatus[prob.status])
# 输出结果
for v in prob.variables():
    print(v.name, v.varValue)
print('min = ', value(prob.objective))

输出结果如下：

status:  Optimal
路线_A_1 300.0
路线_A_2 0.0
路线_A_3 0.0
路线_A_4 0.0
路线_A_5 700.0
路线_B_1 200.0
路线_B_2 900.0
路线_B_3 1800.0
路线_B_4 200.0
路线_B_5 0.0
min =  8600.0

参考

• Python小白的数学建模课 - CSDN博客
• Optimization with PuLP — PuLP 2.6.0 documentation (coin-or.github.io)