这篇文章你将知道什么是LP，LP的建模过程，最后给出一个LP简单实例，并用scipy和ortools两个工具求解。
一般形如
min ? z = ∑ j = 1 n c j x j s . t . { ∑ j = 1 n a i j x j ≤ b i , i = 1 , 2 , ? ? , m x j ≥ 0 , j = 1 , 2 , ? ? , n \begin{equation*} \begin{split} &\min {z= \sum_{j=1}^nc_j x_j}\\ &s.t. \quad \left\{ \begin{array}{lc} \sum\limits_{j=1}^n a_{ij}x_j \leq b_i, i=1,2,\cdots,m\\ \\ x_j\geq 0, j = 1, 2, \cdots, n \end{array} \right. \end{split} \end{equation*} ?minz=j=1∑n?cj?xj?s.t.? ? ??j=1∑n?aij?xj?≤bi?,i=1,2,?,mxj?≥0,j=1,2,?,n???
的数学模型叫LP模型。
LP建模过程遵循一般的最优化建模过程，与我们小学5年级列方程解应用题有点类似，总结如下
任何建模工作都是基于具体的业务问题，首先要对业务有充分的理解，足够的熟稔，最好能够做到驾轻就熟，信手拈来，弄清业务流转和信息节点，在吃透业务逻辑的基础上聚焦要解决的问题，可以是具体业绩问题，也可以是形而上的问题，利用极限思想对问题取极限，让问题动起来变成了目标，有了目标就可以有的放矢，去挖掘哪些因素促进了目标极限过程，是有利的还是有弊的，将这些信息有条理的罗列出来；
我们在做业务数据分析过程中就会发现有些是恒定不变的信息，称为常量，有些是随时间或者业务形态而发生变化的信息，称为变量，但是有些常量和变量的出现他是有条件的，这个时候就可以做出一些适当的假设，以保证这些常量和变量按照预想的方式出现；
在第一步的时候，我就大体知道目标是什么了，又在第二步知道了哪些决策变量，于是可以利用这些决策变量把目标函数的表达式写出来，由于是LP问题，那么目标函数应具备线性形式，如果最初不具备线性形式也没关系，看后期能不能通过一些线性或者非线性变换转为线性形式；
第一步罗列出来的有利信息，在有了决策变量的前提下，可以再进一步加工，通过决策变量写出具体的表达式来作为约束条件；
某家化工厂要生产一种由甲、乙两种原料混合配置而成的产品，而甲、乙两种原料均包含A、B、C三种有效化学成分，甲、乙两种原料中的化学成分含量以及客户需要的最终产品中化学成分配方要求如下表所示
| 化学成分 | 甲原料化学成分含量(%) | 甲原料化学成分含量(%) | 最终产品中化学成分最低含量(%) |
|---|---|---|---|
| A | 12 | 3 | 4 |
| B | 2 | 3 | 2 |
| C | 3 | 15 | 5 |
已知甲，乙两种原料成本分别是5元/千克，4元/千克，请问如何配置原料可使得总成本最小？
设1千克最终产品中有甲原料 x 1 x_1 x1?千克，乙原料 x 2 x_2 x2?千克，那么， x 1 + x 2 = 1 x_1+x_2 =1 x1?+x2?=1，此时的成本是
5 x 1 + 4 x 2 5x_1 +4x_2 5x1?+4x2?
同时，还需要满足各种化学成分的最低含量约束，
A成分的
12
%
x
1
+
3
%
x
2
1
×
100
%
≥
4
%
\frac{12\%x_1+3\%x_2}{1}\times 100\% \geq 4\%
112%x1?+3%x2??×100%≥4%
B成分的
2
%
x
1
+
3
%
x
2
1
×
100
%
≥
2
%
\frac{2\%x_1+3\%x_2}{1}\times 100\% \geq 2\%
12%x1?+3%x2??×100%≥2%
C成分的
3
%
x
1
+
15
%
x
2
1
×
100
%
≥
5
%
\frac{3\%x_1+15\%x_2}{1}\times 100\% \geq 5\%
13%x1?+15%x2??×100%≥5%
整理一下，可以建立如下LP数学模型
min
?
z
=
5
x
1
+
4
x
2
s
.
t
.
{
x
1
+
x
2
=
1
12
x
1
+
3
x
2
≥
4
2
x
1
+
3
x
2
≥
2
3
x
1
+
15
x
2
≥
5
x
j
≥
0
,
j
=
1
,
2
\begin{equation*} \begin{split} &\min {z= 5x_1+4x_2}\\ &s.t. \quad \left\{ \begin{array}{lc} x_1+x_2 =1\\ 12x_1+3x_2\geq 4\\ 2x_1+3x_2\geq 2\\ 3x_1+15x_2\geq 5\\ x_j\geq 0, j = 1, 2 \end{array} \right. \end{split} \end{equation*}
?minz=5x1?+4x2?s.t.?
?
??x1?+x2?=112x1?+3x2?≥42x1?+3x2?≥23x1?+15x2?≥5xj?≥0,j=1,2???
有了数学模型后接下来就是求解的问题了,下面给出2种求解方法
为了利用scipy模块的optimize来进行求解还需要将数学模型化成标准形式
min
?
z
=
5
x
1
+
4
x
2
s
.
t
.
{
x
1
+
x
2
=
1
?
12
x
1
?
3
x
2
≤
?
4
?
2
x
1
?
3
x
2
≤
?
2
?
3
x
1
?
15
x
2
≤
?
5
x
j
≥
0
,
j
=
1
,
2
\begin{equation*} \begin{split} &\min {z= 5x_1+4x_2}\\ &s.t. \quad \left\{ \begin{array}{lc} x_1+x_2 =1\\ -12x_1-3x_2\leq -4\\ -2x_1-3x_2\leq -2\\ -3x_1-15x_2\leq -5\\ x_j\geq 0, j = 1, 2 \end{array} \right. \end{split} \end{equation*}
?minz=5x1?+4x2?s.t.?
?
??x1?+x2?=1?12x1??3x2?≤?4?2x1??3x2?≤?2?3x1??15x2?≤?5xj?≥0,j=1,2???
核心代码
# -*- encoding: utf-8 -*-
'''
@Project : lp
@Desc : 线性规划
@Time : 2024/01/20 14:07:47
@Author : 帅帅de三叔,zengbowengood@163.com
'''
import time
import numpy as np
from scipy.optimize import linprog
start_time = time.time()
c = np.array([5, 4])
A_ub = np.array([[-12, -3], [-2, -3], [-3, -15]])
b_ub = np.array([-4, -2, -5])
A_eq = np.array([[1,1]])
b_eq = np.array([1])
x1_bounds = (0, None)
x2_bounds = (0, None)
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,bounds=[x1_bounds, x2_bounds])
end_time = time.time()
cost_time = end_time - start_time
print("每公斤产品的原料成本：",result.fun)
print("每公斤成品需要甲，乙各多少千克 ", result.x)
print("优化算法信息：", result.message)
print("耗时：", cost_time)
最后输出
每公斤产品的原料成本： 4.111111111111112
每公斤成品需要甲，乙各多少千克 [0.11111111 0.88888889]
优化算法信息： Optimization terminated successfully. (HiGHS Status 7: Optimal)
耗时： 0.004002571105957031
# -*- encoding: utf-8 -*-
'''
@Project : lp
@Desc : 线性规划
@Time : 2024/01/20 14:07:47
@Author : 帅帅de三叔,zengbowengood@163.com
'''
import time
start_time = time.time()
from ortools.linear_solver import pywraplp
solver = pywraplp.Solver.CreateSolver('GLOP')
x1 = solver.NumVar(0, solver.infinity(), 'x1')
x2 = solver.NumVar(0, solver.infinity(), 'x2')
solver.Add(x1+x2 == 1)
solver.Add(12*x1+3*x2>=3)
solver.Add(2*x1+3*x2>=2)
solver.Add(3*x1+15*x2>=5)
solver.Minimize(5*x1+4*x2)
status = solver.Solve()
if status== pywraplp.Solver.OPTIMAL:
print("x1 = ", x1.solution_value())
print("x2 = ", x2.solution_value())
print("最优化结果：",solver.Objective().Value())
else:
print("no optimal solution")
print('\nAdvanced usage:')
print('Problem solved in %f milliseconds' % solver.wall_time())
print('Problem solved in %d iterations' % solver.iterations())
end_time = time.time()
cost_time = end_time - start_time
print("耗时：", cost_time)
最后输出结果
x1 = 0.0
x2 = 1.0
最优化结果： 4.0
Advanced usage:
Problem solved in 64.000000 milliseconds
Problem solved in 0 iterations
耗时： 0.07652115821838379
1，常见的运筹优化类问题及常用的优化算法
https://blog.csdn.net/qq_38384924/article/details/120849695
2，一步一步走向锥规划 - LP
https://zhuanlan.zhihu.com/p/541266572
3，运筹说 第13期 | 线性规划硬核知识点梳理—数学模型&图解法
https://zhuanlan.zhihu.com/p/382644742
4，1-3 LP问题的解
https://www.bilibili.com/video/BV15g411w7ox/
5，运筹优化工具ortools解读与实践-概览
https://segmentfault.com/a/1190000042004080
6，scipy.optimize.linprog
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html
7,运筹优化工具ortools解读与实践-ortools求解LP、IP、MIP问题
https://segmentfault.com/a/1190000042014666?sort=votes