1、工业场景
从很多种铁矿石中选出适合烧结配料的部分铁矿石及其比例,并使其成本最低。
2、数学模型
设Pi代表了第i种原料的成本,xi代表了第i种原料在总配料中的比例,其中i取值为1,2,…,n。计算1吨配料成本:
第种原料的成本是Yi=Pixi1,总成本是Y=∑Yi=∑Pixi,其中i取值为1,2,…,n。目标函数即为:
minY = ∑Yi=∑Pixi
3、使用函数
Scipy.optimize.linprog
def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None,bounds=None, method='interior-point', callback=None,options=None, x0=None):
Parameters----------c : 1-D arrayThe coefficients of the linear objective function to be minimized.A_ub : 2-D array, optionalThe inequality constraint matrix. Each row of ``A_ub`` specifies thecoefficients of a linear inequality constraint on ``x``.b_ub : 1-D array, optionalThe inequality constraint vector. Each element represents anupper bound on the corresponding value of ``A_ub @ x``.A_eq : 2-D array, optionalThe equality constraint matrix. Each row of ``A_eq`` specifies thecoefficients of a linear equality constraint on ``x``.b_eq : 1-D array, optionalThe equality constraint vector. Each element of ``A_eq @ x`` must equalthe corresponding element of ``b_eq``.bounds : sequence, optionalA sequence of ``(min, max)`` pairs for each element in ``x``, definingthe minimum and maximum values of that decision variable. Use ``None``to indicate that there is no bound. By default, bounds are``(0, None)`` (all decision variables are non-negative).If a single tuple ``(min, max)`` is provided, then ``min`` and``max`` will serve as bounds for all decision variables.
c 目标函数的决策变量对应的系数向量
在这里也就是由17种原料的每吨成本价格构成的向量
A_ub 不等式约束组成的决策变量系数矩阵
在这里也就是17种原料的每种化学成分构成的向量
b_ub 由A_ub对应不等式顺序的阈值向量
在这里也就是每种化学成分对应烧结矿的影响
A_eq 等式约束组成的决策变量系数矩阵
b_eq 由A_ub对应等式顺序的阈值向量
这一组是等式约束,在这里17种原料的比例之和等于1
bounds 表示决策变量x连续的定义域的n×2维矩阵,None表示无穷
这里表示的是对应17种原料的控制范围最小值最大值,是一个n×2维矩阵,n为原料的个数17,
method 调用的求解方法
callback 选择的回调函数
options 求解器选择的字典
x0 初始假设的决策变量向量,若可行linprog会对方法优化
4、代码
from scipy.optimize import linprog
import numpy as npc=np.array([946,1062,592,1155,1375,1248,847,519,848,600,300,650,447,62,193,1843,800]).transpose()
A_ub = np.array([[0.6137, 0.6663, 0.5647, 0.6451, 0.7210, 0.6887, 0.6115, 0.5384, 0.6141, 0.5510, 0.6230, 0.7000, 0.00, 0.00, 0.00, 0.00, 0.00],[-0.6137, -0.6663, -0.5647, -0.6451, -0.7210, -0.6887, -0.6115, -0.5384, -0.6141, -0.5510, -0.6230, -0.7000, -0.00, -0.00, -0.00, -0.00, -0.00],[0.0391,0.0507 ,0.0570,0.0210 ,0.0052 ,0.0238 ,0.0465 ,0.0828 ,0.0775,0.0525 ,0.0588 ,0.0255 ,0.0280 ,0.0250 ,0.0580 ,0.4200 ,0.4150 ],[-0.0391,-0.0507 ,-0.0570,-0.0210 ,-0.0052 ,-0.0238 ,-0.0465 ,-0.0828 ,-0.0775,-0.0525 ,-0.0588 ,-0.0255 ,-0.0280 ,-0.0250 ,-0.0580 ,-0.4200 ,-0.4150 ],[0.0003,0.0054 ,0.0005,0.0001 ,0.0007 ,0.0026,0.0002,0.0011 ,0.0020,0.0270 ,0.0134 ,0.0026 ,0.0360 ,0.0175 ,0.4290 ,0.00 ,0.00 ],[-0.0003,-0.0054 ,-0.0005,-0.0001 ,-0.0007 ,-0.0026,-0.0002,-0.0011 ,-0.0020,-0.0270 ,-0.0134 ,-0.0026 ,-0.0360 ,-0.0175 ,-0.4290 ,-0.00 ,-0.00 ],[0.0240,0.0070 ,0.0335,0.0177 ,0.0010 ,0.0066 ,0.0217 ,0.0379 ,0.0190,0.0230 ,0.0132 ,0.0020 ,0.0080 ,0.0080 ,0.0081 ,0.2000 ,0.2800 ],[-0.0240,-0.0070 ,-0.0335,-0.0177 ,-0.0010 ,-0.0066 ,-0.0217 ,-0.0379 ,-0.0190,-0.0230 ,-0.0132 ,-0.0020 ,-0.0080 ,-0.0080 ,-0.0081 ,-0.2000 ,-0.2800 ],[0.00096 ,0.00019 ,0.00041 ,0.00114 ,0.00021 ,0.00010 ,0.00055 ,0.00069 ,0.00051 ,0.00080 ,0.00334 ,0.00013 ,0.00 ,0.00 ,0.00 ,0.0003 ,0.0004 ],[-0.00096 ,-0.00019 ,-0.00041 ,-0.00114 ,-0.00021 ,-0.00010 ,-0.00055 ,-0.00069 ,-0.00051 ,-0.00080 ,-0.00334 ,-0.00013 ,-0.00 ,-0.00 ,-0.00 ,-0.0003 ,-0.0004 ],[0.00001 ,0.0011,0.01520 ,0.0002 ,0.00490 ,0.00060 ,0.00230 ,0.000 ,0.00420 ,0.00400 ,0.01680 ,0.00480 ,0.00 ,0.00 ,0.00 ,0.00 ,0.00 ],[-0.00001 ,-0.0011,-0.01520 ,-0.0002 ,-0.00490 ,-0.00060 ,-0.00230 ,-0.000 ,-0.00420 ,-0.00400 ,-0.01680 ,-0.00480 ,-0.00 ,-0.00 ,-0.00 ,-0.00 ,-0.00]])
b_ub = np.array([57,-55,7,-4,3.5,-1.5,3,-1.5,1.1,-0.04,0.7,-0.2])A_eq = np.array([[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]])
b_eq = np.array([100])x=linprog(c,A_ub,b_ub,A_eq,b_eq,method='highs',bounds=np.array([[10,20],[1,10],[1,10],[1,10],[1,10],[1,10],[1,10],[1,10],[10,20],[10,20],[0,1],[0,2],[5,8],[0,2],[3,5],[3,5],[0.5,2]]))
print(x)
5、计算结果
con: array([0.])crossover_nit: 0eqlin: marginals: array([-2092.4011976])residual: array([0.])fun: 86687.33832335332ineqlin: marginals: array([ -0. , -4850.2994012, -0. , -0. ,-0. , -0. , -0. , -0. ,-0. , -0. , -0. , -0. ])residual: array([2. , 0. , 0.7905684 , 2.2094316 , 1.62490589,0.37509411, 0.58011713, 0.91988287, 1.05100968, 0.00899032,0.31292236, 0.18707764])lower: marginals: <MemoryView of 'ndarray' at 0xb0686e8>residual: array([ 0. , 9. , 9. , 0. , 9. ,1.29607452, 9. , 1.20392548, 10. , 0. ,1. , 2. , 0. , 0. , 0. ,0. , 0. ])message: 'Optimization terminated successfully.'nit: 9slack: array([2. , 0. , 0.7905684 , 2.2094316 , 1.62490589,0.37509411, 0.58011713, 0.91988287, 1.05100968, 0.00899032,0.31292236, 0.18707764])status: 0success: Trueupper: marginals: <MemoryView of 'ndarray' at 0xb068540>residual: array([10. , 0. , 0. , 9. , 0. ,7.70392548, 0. , 7.79607452, 0. , 10. ,0. , 0. , 3. , 2. , 2. ,2. , 1.5 ])x: array([10. , 10. , 10. , 1. , 10. ,2.29607452, 10. , 2.20392548, 20. , 10. ,1. , 2. , 5. , 0. , 3. ,3. , 0.5 ])Process finished with exit code 0
最后的x: array就是我们想要的值
[10., 10., 10., 1., 10.,2.29607452, 10., 2.20392548, 20., 10.,1., 2., 5., 0., 3.,3., 0.5]
6、带入实际验证公式
铁元素含量
tensor(0.5879)
吨度价
tensor(17.1714)
7、同一批料对比人工测算
铁元素含量
tensor(0.5510)
吨度价
tensor(17.6456)