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

看完上述內容是否對您有幫助呢？如果還想對相關知識有進一步的了解或閱讀更多相關文章，請關注億速云行業資訊頻道，感謝您對億速云的支持。