运用python scipy来求解线性规划问题

运用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......................

看完上述内容是否对您有帮助呢？如果还想对相关知识有进一步的了解或阅读更多相关文章，请关注亿速云行业资讯频道，感谢您对亿速云的支持。