运用python scipy来求解线性规划问题,很多新手对此不是很清楚,为了帮助大家解决这个难题,下面小编将为大家详细讲解,有这方面需求的人可以来学习下,希望你能有所收获。
最近项目中遇到了一个有意思的技术点,运用线性规划数学模型,求最优解的问题。
由于我后台是springcloud,所以我调查到两种实现方式。
第一种是运用python的scipy开源库写一个python脚本,然后java后台调用python脚本,求最优解,然后再将结果返回。
第二种是运用java中ibm组件Cplex直接求解,但分为收费和免费版,免费版决策变量限制为1000个,但也够了。找jar包和dll费劲点。
发现他们俩算出来的最优解是相同的,但各个决策变量不太相同。
这里记录一下运用python求解的方法。
项目中的不等式方程组是∑求和形式的,这里就手动打码先假设决策变量的数量为 5 * 3 个,经过化简后得到多元一次不等式组如下:
求解最大值 = AX + AY + AZ + BX + BY + BZ + CX + CY + CZ + DX + DY + DZ + EX + EY + EZ AX + AY + AZ <= 25400 BX + BY + BZ <= 18600 CX + CY + CZ <= 39800 DX + DY + DZ <= 53200 EX + EY + EZ <= 5900 AX + AY + AZ <= 10000 BX + BY + BZ <= 10000 CX + CY + CZ <= 10000 DX + DY + DZ <= 10000 EX + EY + EZ <= 10000 AX + BX + CX + DX + EX >= 15000 AY + BY + CY + DY + DY >= 5000 AZ + BZ + CZ + DZ + DZ >= 10000 (50.25-50)*AX + (49.86-50)*BX + (68.80-50)*CX + (49.79-50)*DX + (48.77-50)*EX >= 0 (50.25-60)*AY + (49.86-60)*BY + (68.80-60)*CY + (49.79-60)*DY + (48.77-60)*EY >= 0 (50.25-55)*AZ + (49.86-55)*BZ + (68.80-55)*CZ + (49.79-55)*DZ + (48.77-55)*DZ >= 0 (30.95*(1-2/100)-30)*AX + (31.52*(1-3/100)-30)*BX + (30.58*(1-1/100)-30)*CX + (30.17*(1-1/100)-30)*DX + (27.83*(1-1/100)-30)*EX >= 0 (30.95*(1-2/100)-30)*AY + (31.52*(1-3/100)-30)*BY + (30.58*(1-1/100)-30)*CY + (30.17*(1-1/100)-30)*DY + (27.83*(1-1/100)-30)*EY >= 0 (30.95*(1-2/100)-30)*AZ + (31.52*(1-3/100)-30)*BZ + (30.58*(1-1/100)-30)*CZ + (30.17*(1-1/100)-30)*DZ + (27.83*(1-1/100)-30)*EZ >= 0 (11.32*(1-2/100)-10)*AX + (12.83*(1-3/100)-10)*BX + (16.06*(1-1/100)-10)*CX + (5.68*(1-1/100)-10)*DX + (8.54*(1-1/100)-10)*EX >= 0 (11.32*(1-2/100)-10)*AY + (12.83*(1-3/100)-10)*BY + (16.06*(1-1/100)-10)*CY + (5.68*(1-1/100)-10)*DY + (8.54*(1-1/100)-10)*EY >= 0 (11.32*(1-2/100)-10)*AZ + (12.83*(1-3/100)-10)*BZ + (16.06*(1-1/100)-10)*CZ + (5.68*(1-1/100)-10)*DZ + (8.54*(1-1/100)-10)*EZ >= 0 (6*(1-2/100)-5)*AX + (4*(1-3/100)-5)*BX + (5*(1-1/100)-5)*CX + (2*(1-1/100)-5)*DX + (5*(1-1/100)-5)*EX <= 0 (6*(1-2/100)-5)*AY + (4*(1-3/100)-5)*BY + (5*(1-1/100)-5)*CY + (2*(1-1/100)-5)*DY + (5*(1-1/100)-5)*EY <= 0 (6*(1-2/100)-5)*AZ + (4*(1-3/100)-5)*BZ + (5*(1-1/100)-5)*CZ + (2*(1-1/100)-5)*DZ + (5*(1-1/100)-5)*EZ <= 0 非负约束: AX,AY,AZ,BX,BY,BZ,CX,CY,CZ,DX,DY,DZ,EX,EY,EZ >= 0
上面这些是根据我自己的项目得到的不等式组,可以根据自己的项目来做相应改动。
# coding=utf-8 import json import sys import numpy as np from scipy import optimize as op def solve(aaaaa, bbbbb, ccccc): # 定义目标函数系数矩阵 c = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # 定义不等式约束左端系数矩阵 a_ub = np.array([ [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1], [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1], [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1], [-1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0], [0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0], [0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1], [50 - 50.25, 0, 0, 50 - 49.86, 0, 0, 50 - 68.80, 0, 0, 50 - 49.79, 0, 0, 50 - 48.77, 0, 0], [0, 60 - 50.25, 0, 0, 60 - 49.86, 0, 0, 60 - 68.80, 0, 0, 60 - 49.79, 0, 0, 60 - 48.77, 0], [0, 0, 55 - 50.25, 0, 0, 55 - 49.86, 0, 0, 55 - 68.80, 0, 0, 55 - 49.79, 0, 0, 55 - 48.77], [30 - 30.95 * (1 - 2 / 100), 0, 0, 30 - 31.52 * (1 - 3 / 100), 0, 0, 30 - 30.58 * (1 - 1 / 100), 0, 0, 30 - 30.17 * (1 - 1 / 100), 0, 0, 30 - 27.83 * (1 - 1 / 100), 0, 0], [0, 30 - 30.95 * (1 - 2 / 100), 0, 0, 30 - 31.52 * (1 - 3 / 100), 0, 0, 30 - 30.58 * (1 - 1 / 100), 0, 0, 30 - 30.17 * (1 - 1 / 100), 0, 0, 30 - 27.83 * (1 - 1 / 100), 0], [0, 0, 30 - 30.95 * (1 - 2 / 100), 0, 0, 30 - 31.52 * (1 - 3 / 100), 0, 0, 30 - 30.58 * (1 - 1 / 100), 0, 0, 30 - 30.17 * (1 - 1 / 100), 0, 0, 30 - 27.83 * (1 - 1 / 100)], [10 - 11.32 * (1 - 2 / 100), 0, 0, 10 - 12.83 * (1 - 3 / 100), 0, 0, 10 - 16.06 * (1 - 1 / 100), 0, 0, 10 - 5.68 * (1 - 1 / 100), 0, 0, 10 - 8.54 * (1 - 1 / 100), 0, 0], [0, 10 - 11.32 * (1 - 2 / 100), 0, 0, 10 - 12.83 * (1 - 3 / 100), 0, 0, 10 - 16.06 * (1 - 1 / 100), 0, 0, 10 - 5.68 * (1 - 1 / 100), 0, 0, 10 - 8.54 * (1 - 1 / 100), 0], [0, 0, 10 - 11.32 * (1 - 2 / 100), 0, 0, 10 - 12.83 * (1 - 3 / 100), 0, 0, 10 - 16.06 * (1 - 1 / 100), 0, 0, 10 - 5.68 * (1 - 1 / 100), 0, 0, 10 - 8.54 * (1 - 1 / 100)], [6 * (1 - 2 / 100) - 5, 0, 0, 4 * (1 - 3 / 100) - 5, 0, 0, 5 * (1 - 1 / 100) - 5, 0, 0, 2 * (1 - 1 / 100) - 5, 0, 0, 5 * (1 - 1 / 100) - 5, 0, 0], [0, 6 * (1 - 2 / 100) - 5, 0, 0, 4 * (1 - 3 / 100) - 5, 0, 0, 5 * (1 - 1 / 100) - 5, 0, 0, 2 * (1 - 1 / 100) - 5, 0, 0, 5 * (1 - 1 / 100) - 5, 0], [0, 0, 6 * (1 - 2 / 100) - 5, 0, 0, 4 * (1 - 3 / 100) - 5, 0, 0, 5 * (1 - 1 / 100) - 5, 0, 0, 2 * (1 - 1 / 100) - 5, 0, 0, 5 * (1 - 1 / 100) - 5] ]) # 定义不等式约束右端项矩阵 # 非常接近0的值 e = 1e-10 b_ub = np.array( [25400, 18600, 39800, 53200, 5900, 10000, 10000, 10000, 10000, 10000, 20000, 20000, 20000, 20000, 20000, -15000, -5000, -10000, -e, -e, -e, -e, -e, -e, -e, -e, -e, e, e, e]) # 取值范围 ax = ay = az = bx = by = bz = cx = cy = cz = dx = dy = dz = ex = ey = ez = (0, 10000) bounds = (ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez) # 调用函数进行求解,默认求解最小值,求解最大值使用-c并取结果相反数,maxiter最大迭代次数 return op.linprog(-c, a_ub, b_ub, bounds=bounds, options={'maxiter': 50}) if __name__ == '__main__': aaaaa= """111""" bbbbb= """222""" ccccc= """333""" # sys.argv[i] res = solve(aaaaa, bbbbb, ccccc) print(res) exit() # 执行结果 # con: array([], dtype=float64) # fun: -44280.11226613013 # message: 'Optimization terminated successfully.' # nit: 8 # slack: array([1.54000001e+04, 8.60000008e+03, 2.98000001e+04, 4.32000001e+04, # 1.61988736e+03, 8.62767956e-05, 8.33909326e-05, 1.16868736e-04, # 8.36751933e-05, 5.71988736e+03, 1.00000001e+04, 1.00000001e+04, # 1.00000001e+04, 1.00000001e+04, 1.57198874e+04, 6.58111172e+03, # 3.95477828e+03, 3.74422226e+03, 1.21703633e+04, 4.95361919e+03, # 6.34258264e+03, 6.66164697e-05, 6.74106730e-05, 9.34142053e-07, # 7.20490633e+03, 2.40991082e+04, 1.26944988e+04, 1.52974587e+04, # 5.68045285e+03, 1.23360938e+04]) # status: 0 # success: True # x: array([6771.5673779 , 652.76382833, 2575.6687075 , 7211.86997104, # 773.63550554, 2014.49444003, 820.59747095, 5091.61932829, # 4087.78308389, 4309.51511365, 1687.54017157, 4002.9446311 , # 2467.561791 , 749.21944752, 1063.33139782])
需要注意的点:
一、在定义不等式左端和右端矩阵时,默认不等式都是小于等于的,可以先将不等式反转。
要不就直接取相反数(我就是)。类似这段中不等式左端和右端在定义矩阵时都取了相反数:
AX + BX + CX + DX + EX >= 15000
二、在定义不等式左端和右端矩阵时,先将决策变量位置对齐,要是对齐的位置上没有就用0替代。对齐后的不等式:
AX + 0AY + 0AZ + BX + 0BY + 0BZ + CX + 0CY + 0CZ + DX + 0DY + 0DZ + EX + 0EY + 0EZ >= 15000
三、在求解时默认计算最小值,如果想计算最大值就加一个负号,算出来的值把负号去掉就是最大值。类似这段中的 -c:
return op.linprog(-c, a_ub......................
看完上述内容是否对您有帮助呢?如果还想对相关知识有进一步的了解或阅读更多相关文章,请关注亿速云行业资讯频道,感谢您对亿速云的支持。
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。