运筹优化-LP

原创已于 2024-01-21 15:51:24 修改 · 1.1k 阅读

21 ·

CC 4.0 BY-SA版权

文章标签：

#python #线性规划 #运筹优化 #ortools #scipy

于 2024-01-20 15:41:18 首次发布

LP的定义

一般形如

$\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*}$

的数学模型叫LP模型。

LP的建模过程

LP建模过程遵循一般的最优化建模过程，与我们小学5年级列方程解应用题有点类似，总结如下

分析业务，理解要义，区分目标和条件

任何建模工作都是基于具体的业务问题，首先要对业务有充分的理解，足够的熟稔，最好能够做到驾轻就熟，信手拈来，弄清业务流转和信息节点，在吃透业务逻辑的基础上聚焦要解决的问题，可以是具体业绩问题，也可以是形而上的问题，利用极限思想对问题取极限，让问题动起来变成了目标，有了目标就可以有的放矢，去挖掘哪些因素促进了目标极限过程，是有利的还是有弊的，将这些信息有条理的罗列出来；

定义决策变量，适当假设

我们在做业务数据分析过程中就会发现有些是恒定不变的信息，称为常量，有些是随时间或者业务形态而发生变化的信息，称为变量，但是有些常量和变量的出现他是有条件的，这个时候就可以做出一些适当的假设，以保证这些常量和变量按照预想的方式出现；

写出目标函数

在第一步的时候，我就大体知道目标是什么了，又在第二步知道了哪些决策变量，于是可以利用这些决策变量把目标函数的表达式写出来，由于是LP问题，那么目标函数应具备线性形式，如果最初不具备线性形式也没关系，看后期能不能通过一些线性或者非线性变换转为线性形式；

列出约束条件

第一步罗列出来的有利信息，在有了决策变量的前提下，可以再进一步加工，通过决策变量写出具体的表达式来作为约束条件；

标准型

LP的标准型应该满足(1)目标函数最小，（2）约束条件为等式；（3）右端常数非负；（4）决策变量非负，按照标准型的要求，我们将定义里面的形式改写成标准形式如下
$\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 +x_{j+i}= b_i, i=1,2,\cdots,m\\ \\ x_j\geq 0, j = 1, 2, \cdots, n \end{array} \right. \end{split} \end{equation*}$

单纯形算法

单纯形算法和IPM算法是LP求解最常用的两类算法，这里先给出单纯形算法，单纯形算法的一般步骤

1，将LP的数学模型转为标准型；

2，找到一个基本可行解, 寻找目标函数的替入变量, 约束条件的替出变量, 替换目标函数, 不断旋转；

3，重复步骤2, 直到目标函数的系数中没有负系数为止，算法收敛，结束。

对偶性

后续补充

实例

某家化工厂要生产一种由甲、乙两种原料混合配置而成的产品，而甲、乙两种原料均包含A、B、C三种有效化学成分，甲、乙两种原料中的化学成分含量以及客户需要的最终产品中化学成分配方要求如下表所示

化学成分	甲原料化学成分含量(%)	甲原料化学成分含量(%)	最终产品中化学成分最低含量(%)
A	12	3	4
B	2	3	2
C	3	15	5

已知甲，乙两种原料成本分别是5元/千克，4元/千克，请问如何配置原料可使得总成本最小？

设1千克最终产品中有甲原料 $x_1$ 千克，乙原料 $x_2$ 千克，那么， $x_1+x_2 =1$ ，此时的成本是

$5x_1 +4x_2$

同时，还需要满足各种化学成分的最低含量约束，

A成分的
$\frac{12\%x_1+3\%x_2}{1}\times 100\% \geq 4\%$

B成分的
$\frac{2\%x_1+3\%x_2}{1}\times 100\% \geq 2\%$

C成分的
$\frac{3\%x_1+15\%x_2}{1}\times 100\% \geq 5\%$

整理一下，可以建立如下LP数学模型

$\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*}$
有了数学模型后接下来就是求解的问题了,下面给出2种求解方法

scipy.optimize

为了利用scipy模块的optimize来进行求解还需要将数学模型化成标准形式
$\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*}$
核心代码

# -*- 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

ortools

# -*- 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/83327579
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
8，LP线性规划求解之单纯形算法
https://www.cnblogs.com/chenjieyouge/p/11901150.html
9，运筹总结
https://zhuanlan.zhihu.com/p/480995829?utm_id=0