【數學建模】線性規划各種問題的Python調包方法

本文轉載自查看原文 2019-07-12 12:24 1003 python學習/ 數學建模與數學

關鍵詞：Python、調包、線性規划、指派問題、運輸問題、pulp、混合整數線性規划（MILP）

注：此文章是線性規划的調包實現，具體步驟原理請搜索具體解法。

本文章的各個問題可能會采用多種調用方法，為什么？因為這些包各有特點，有些語法特別像matlab，只要稍稍改變即可達成代碼交換；而有些包利用了python本身的特性，在靈活度與代碼的可讀性上更高。我認為這些包各有優劣，各位各持所需吧。

看了本文章能做到什么？你可以在本文章內學到線性規划的幾個問題的求解方式，並學會如何用pulp包解決線性規划問題。無論是整數規划（Integer Program）、01規划（Binary Program）還是混合整數線性規划（MILP），你都可以得到很好的解題方法。

一、線性規划

該問題引用自《數學建模算法與應用-司守奎》第一章線性規划 3.線性規划
包的具體使用可參考scipy官網

求解最普通的線性規划問題： $\min{z} = 2x_1+3x_2+x_3 \\ \begin{equation} \left\{ \begin{array}{lr} x_1 + 4x_2+2x_3 \geq 8 \\ 3x_1 + 2x_2 \geq 6 \\ x_1,x_2,x_3 \geq 0 \end{array} \right. \end{equation}$

scipy調包代碼

import numpy as np z = np.array([2, 3, 1]) a = np.array([[1, 4, 2], [3, 2, 0]]) b = np.array([8, 6]) x1_bound = x2_bound = x3_bound =(0, None) from scipy import optimize res = optimize.linprog(z, A_ub=-a, b_ub=-b,bounds=(x1_bound, x2_bound, x3_bound)) print(res) #output: # fun: 7.0 # message: 'Optimization terminated successfully.' # nit: 2 # slack: array([0., 0.]) # status: 0 # success: True # x: array([0.8, 1.8, 0. ])

注意，此函數和matlab的linprog有很多相似之處。
首先默認就是求解最小值，如果想要求最大值，需要對目標函數的各參數取反（既對z取反），並在得出的結果（func）前取反。
並且所有約束條件都為≤，因此例子中傳入參數是前面都加了負號。
bound為邊界的二元元組，None時為無邊界。
如果存在類似 $x_1+2x_2+4x_3=101$ 這種情況，可以：

A_eq = [[1,2,4]] b_eq = [101]

並在linprog中傳入。
得出的結果為scipy.optimize.optimize.OptimizeResult（優化結果）類型，是類似字典的結構，例如想要優化結果代入目標函數的值，可以res.fun或res['fun']這樣取值。

pulp調包代碼

import pulp #目標函數的系數 z = [2, 3, 1] #約束 a = [[1, 4, 2], [3, 2, 0]] b = [8, 6] #確定最大化最小化問題，最大化只要把Min改成Max即可 m = pulp.LpProblem(sense=pulp.LpMinimize) #定義三個變量放到列表中 x = [pulp.LpVariable(f'x{i}', lowBound=0) for i in [1,2,3]] #定義目標函數，lpDot可以將兩個列表的對應位相乘再加和 #相當於z[0]*x[0]+z[1]*x[0]+z[2]*x[2] m += pulp.lpDot(z, x) #設置約束條件 for i in range(len(a)): m += (pulp.lpDot(a[i], x) >= b[i]) #求解 m.solve() #輸出結果 print(f'優化結果：{pulp.value(m.objective)}') print(f'參數取值：{[pulp.value(var) for var in x]}') #output: #優化結果：7.0 #參數取值：[2.0, 0.0, 3.0]

每一步的說明已經注釋在代碼中，可以看到輸出結果，兩者的變量取值並不一致，但代入目標函數的結果都是一樣的。
同樣的，如果存在類似 $x_1+2x_2+4x_3=101$ 這種情況，可以：

A_eq = [1,2,4] b_eq = 101 m += (pulp.lpDot(A_eq, x) == b_eq)

二、運輸問題

某商品有 $m$ 個產地、 $n$ 個銷地，各產地的產量分別為 $a_1, ..., a_m$ ，各銷地的需求量分別為 $b_1,...,b_n$ 。若該商品由 $i$ 產地運到 $j$ 銷地的單位運價為 $c_{ij}$ ，問應該如何調運才能使總運費最省？
引入變量 $x_{ij}$ ，其取值為由 $i$ 產地運往 $j$ 銷地的該商品數量，數學模型為 : $\Large \min {\sum_{i=1}^m\sum_{j=1}^{n}c_{ij}x_{ij}} \\ \left\{ \begin{array}{lr} \Large \sum_{j=1}^n x_{ij} = a_i,\quad i=1,...,m\\ \Large \sum_{i=1}^m x_{ij} = b_j, \quad j=1,2,..., n \\ \Large x_{ij} \geq 0 \end{array} \right.$
例題：

pulp調包代碼

import pulp import numpy as np from pprint import pprint def transportation_problem(costs, x_max, y_max): row = len(costs) col = len(costs[0]) prob = pulp.LpProblem('Transportation Problem', sense=pulp.LpMaximize) var = [[pulp.LpVariable(f'x{i}{j}', lowBound=0, cat=pulp.LpInteger) for j in range(col)] for i in range(row)] flatten = lambda x: [y for l in x for y in flatten(l)] if type(x) is list else [x] prob += pulp.lpDot(flatten(var), costs.flatten()) for i in range(row): prob += (pulp.lpSum(var[i]) <= x_max[i]) for j in range(col): prob += (pulp.lpSum([var[i][j] for i in range(row)]) <= y_max[j]) prob.solve() return {'objective':pulp.value(prob.objective), 'var': [[pulp.value(var[i][j]) for j in range(col)] for i in range(row)]}

然后構造參數傳遞進去：

if __name__ == '__main__': costs = np.array([[500, 550, 630, 1000, 800, 700], [800, 700, 600, 950, 900, 930], [1000, 960, 840, 650, 600, 700], [1200, 1040, 980, 860, 880, 780]]) max_plant = [76, 88, 96, 40] max_cultivation = [42, 56, 44, 39, 60, 59] res = transportation_problem(costs, max_plant, max_cultivation) print(f'最大值為{res["objective"]}') print('各變量的取值為：') pprint(res['var']) #output: #最大值為284230.0 #各變量的取值為： #[[0.0, 0.0, 6.0, 39.0, 31.0, 0.0], # [0.0, 0.0, 0.0, 0.0, 29.0, 59.0], # [2.0, 56.0, 38.0, 0.0, 0.0, 0.0], # [40.0, 0.0, 0.0, 0.0, 0.0, 0.0]]

三、指派問題

該問題引用自《數學建模算法與應用-司守奎》第一章線性規划 3.指派問題
調包解決方法參考https://blog.csdn.net/your_answer/article/details/79160045
可參考官網https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.linear_sum_assignment.html

擬分配n人去干 $n$ 項工作，沒人干且僅干一項工作，若分配第 $i$ 人去干第 $j$ 項工作，需花費 $c_{ij}$ 單位時間，問應如何分配工作才能使公認花費的總時間最少？
假設指派問題的系數矩陣為： $C= \left[ \begin{matrix} 12 & 7&9&7&9\\ 8&9&6&6&6 \\ 7&17&12&14&12\\ 15&14&6&6&10\\ 4&10&7&10&6 \end{matrix} \right]$ 引入變量 $x_{ij}$ ，若分配 $i$ 干 $j$ 工作，則取 $x_{ij}=1$ ，否則取 $x_{ij}=0$ ，上述指派問題的數學模型為 $\min {\sum_{i=1}^{n} {\sum_{j=1}^{n}{c_{ij}x_{ij}}}} \\ s.t. \sum_{j=1}^{n}x_{ij}=1 \\ \qquad \sum_{i=1}^{n}x_{ij}=1 \\ \qquad x_{ij}=0或1$ 指派問題的可行解用矩陣表示，每行每列有且只有一個元素為1，其余元素均為0。

調用scipy解決

原書使用匈牙利算法解決的，在這里我們用scipy的優化模塊解決

import numpy as np from scipy.optimize import linear_sum_assignment

引入包，linear_sum_assignment是專門用於解決指派問題的模塊。

efficiency_matrix = np.array([ [12,7,9,7,9], [8,9,6,6,6], [7,17,12,14,12], [15,14,6,6,10], [4,10,7,10,6] ]) row_index, col_index=linear_sum_assignment(efficiency_matrix) print(row_index+1) print(col_index+1) print(efficiency_matrix[row_index,col_index]) #output： #[1 2 3 4 5] #[2 3 1 4 5] #[7 6 7 6 6]

定義了開銷矩陣（指派問題的系數矩陣）efficiency_matrix，傳入linear_sum_assignment，結果返回的是最優指派的行和列，例如第一行選擇第二列，意為：將第一個人派往第二個工作。而根據numpy.array的性質，傳入行和列就會返回行列所對應的值，即為輸出的第三列

print(efficiency_matrix[row_index, col_index].sum())
#output： # 32

對其求和，即可得到指派問題的最小時間。

調用pulp解決

先定義通用解決方法，其中的flatten是遞歸展開列表用的。

def assignment_problem(efficiency_matrix): row = len(efficiency_matrix) col = len(efficiency_matrix[0]) flatten = lambda x: [y for l in x for y in flatten(l)] if type(x) is list else [x] m = pulp.LpProblem('assignment', sense=pulp.LpMinimize) var_x = [[pulp.LpVariable(f'x{i}{j}', cat=pulp.LpBinary) for j in range(col)] for i in range(row)] m += pulp.lpDot(efficiency_matrix.flatten(), flatten(var_x)) for i in range(row): m += (pulp.lpDot(var_x[i], [1]*col) == 1) for j in range(col): m += (pulp.lpDot([var_x[i][j] for i in range(row)], [1]*row) == 1) m.solve() print(m) return {'objective':pulp.value(m.objective), 'var': [[pulp.value(var_x[i][j]) for j in range(col)] for i in range(row)]}

然后定義矩陣，輸入求解

efficiency_matrix = np.array([
    [12, 7, 9, 7, 9], [8, 9, 6, 6, 6], [7, 17, 12, 14, 9], [15, 14, 6, 6, 10], [4, 10, 7, 10, 9] ]) res = assignment_problem(efficiency_matrix) print(f'最小值{res["objective"]}') print(res['var']) #output #最小值32.0 #[[0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0, 0.0]]

四、pulp的使用方式

基本使用

可以看出，pulp在線性規划這一部分，有更多的通用性，編寫程序更自由。前面的例子已經足夠豐富了，但是如果讀到這里還不去清楚pulp的使用方式的話……可以去百度一下，我這里也簡單講一講。

首先是定義一個問題，第一個參數為問題的名稱，不過並非必要的參數，而通過sense參數可決定優化的是最大值還是最小值

prob = pulp.LpProblem('problem name', sense=pulp.LpMinimize)

然后是定義變量：

x0 = pulp.LpVariable('x0', lowBound=0, upBound=None, cat=pulp.LpInteger) x1 = pulp.LpVariable('x1', lowBound=0, upBound=None, cat=pulp.LpInteger) x2 = pulp.LpVariable('x2', lowBound=0, upBound=None, cat=pulp.LpInteger)

這里定義了三個變量，第一個參數為變量名，lowBound 、upBound為上下限，cat為變量類型，默認為連續變量，還可以設為離散變量或二值變量，具體怎么設置？
如上述代碼所示，pulp.LpInteger為離散變量，pulp.LpBinary為二值變量，同時也可以傳入'Integer'字符串的方式指明變量類型。從上面幾個問題的代碼可以看出，我幾乎沒有單個定義變量，而是批量定義的。
然后是添加目標函數：

prob += 2*x0-5*x1+4*x2

只要讓問題（prob）直接加變量的表達式即可添加目標函數。
再之后是添加約束：

prob += (x0+x1-6*x2 <= 120)

只要讓問題（prob）直接加變量的判斷式即可添加約束

prob.solve()

調用solve方法解出答案，如果省去這一步，后面的變量和結果值都會顯示None。

print(pulp.value(prob.objective))
print(pulp.value(x0))

打印優化結果，並顯示當優化達到結果時x0的取值。

思考程序本質

problem對象是如何通過不斷加來獲得目標函數和約束的？熟悉python或者c++的朋友可能會想到一個詞：操作符重載。
沒錯，就是這么實現的，上述的對象幾乎都實現了不同的重載。
首先是Problem對象prob，全名pulp.pulp.LpProblem；當打印輸出（print）時，會打印出問題名，當不斷增加目標函數、約束時，也會隨着print輸出；而它的__add__一定是被定義過了，我們先說其他對象。
當我們定義一個變量時，它的類型是pulp.pulp.LpVariable，當我們對這些變量和其他變量做加法、和其他常數做乘法時，它會返回一個新的對象，經檢測，這個新對象的類名叫pulp.pulp.LpAffineExpression，顧名思義，叫做關系表達式；如果print，會打印這個關系表達式。
而如果對關系表達式做出：>=、<=、==時，會返回新的對象，類名叫做pulp.pulp.LpConstraint，即約束對象；如果print，會打印這個約束。
將關系表達式（pulp.pulp.LpAffineExpression）與問題（pulp.pulp.LpProblem）相加時，會返回新的問題對象，並且新的問題對象會將這個關系表達式作為目標函數。
將約束對象（pulp.pulp.LpConstraint）與問題（pulp.pulp.LpProblem）相加時，會返回新的問題對象，並且這個新的問題對象會多出約束對象所代表的約束條件。
調用問題對象的solve方法，解出線性規划的解。
訪問問題對象的objective成員變量，會得到目標函數（關系表達式對象）。
調用pulp的value方法，將獲得對變量代入值的結果，如果是關系表達式對象，將獲得優化結果；如果是變量對象，將獲得優化結果達到時的變量取值；如果是None，說明你忘調用solve了。

這個包的使用，我是在google中查到的，如果有興趣，可以去原地址看看，擁有更多的應用以及其他包的調用https://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721
也可以去bilibili看看相關的視頻https://www.bilibili.com/video/av18818604?from=search&seid=17463693264345909709

作者：crossous
鏈接：https://www.jianshu.com/p/9be417cbfebb
來源：簡書
簡書著作權歸作者所有，任何形式的轉載都請聯系作者獲得授權並注明出處。

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 Python數學建模系列（一）：規划問題之線性規划 Python數學建模系列（三）：規划問題之非線性規划數學建模（線性規划）python版【數學建模】基於問題的線性規划和混合整數規划求解 Python小白的數學建模課-03.線性規划【數學建模】day01-線性規划問題數學規划模型——線性規划問題 [數學建模(五)]線性規划，二次規划和非線性規划的MATLAB求解 python玄學建模(2)：非線性規划 Matlab數學建模學習筆記——非線性規划