Python SciPy optimize.linprog用法及代碼示例

本文簡要介紹 python 語言中 scipy.optimize.linprog 的用法。

用法:

scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='highs', callback=None, options=None, x0=None, integrality=None)#

線性規劃：最小化受線性等式和不等式約束的線性目標函數。

線性規劃解決以下形式的問題：

\[\begin{split}\min_x \ & c^T x \\ \mbox{such that} \ & A_{ub} x \leq b_{ub},\\ & A_{eq} x = b_{eq},\\ & l \leq x \leq u ,\end{split}\]

其中\(x\)是決策變量的向量； \(c\) 、 \(b_{ub}\) 、 \(b_{eq}\) 、 \(l\) 和 \(u\) 是向量；和\(A_{ub}\) 和\(A_{eq}\) 是矩陣。

或者，那是：

• minimize

  c @ x

• such that

  A_ub @ x <= b_ub
  A_eq @ x == b_eq
  lb <= x <= ub

請注意，默認情況下為 lb = 0 和 ub = None 。其他邊界可以使用 bounds 指定。

參數 ：：

c： 一維數組

要最小化的線性目標函數的係數。

A_ub： 二維陣列，可選

不等式約束矩陣。 A_ub 的每一行指定 x 上的線性不等式約束的係數。

b_ub： 一維數組，可選

不等式約束向量。每個元素代表 A_ub @ x 對應值的上限。

A_eq： 二維陣列，可選

等式約束矩陣。 A_eq 的每一行指定 x 上的線性等式約束的係數。

b_eq： 一維數組，可選

等式約束向量。 A_eq @ x 的每個元素必須等於 b_eq 的相應元素。

bounds： 順序，可選

x 中每個元素的 (min, max) 對序列，定義該決策變量的最小值和最大值。如果提供單個元組 (min, max)，則 min 和 max 將充當所有決策變量的邊界。使用None表示沒有界限。例如，默認界限(0, None)意味著所有決策變量都是非負的，而(None, None)對意味著根本沒有界限，即所有變量都允許為任何實數。

method： str，可選

用於解決標準形式問題的算法。支持‘highs’(默認)、“highs-ds”、“highs-ipm”、“interior-point”(舊版)、“修訂版單工”(舊版)和‘simplex’(舊版)。舊方法已棄用，並將在 SciPy 1.11.0 中刪除。

callback： 可調用的，可選的

如果提供了回調函數，則每次算法迭代將至少調用一次。回調函數必須接受由以下字段組成的單個 scipy.optimize.OptimizeResult ：

x 一維數組

當前解向量。

樂趣 浮點數

目標函數 c @ x 的當前值。

成功 bool

True 當算法成功完成時。

鬆弛 一維數組

鬆弛的(名義上為正)值 b_ub - A_ub @ x 。

騙局 一維數組

等式約束的(名義上為零)殘差 b_eq - A_eq @ x 。

階段 int

正在執行的算法的階段。

狀態 int

表示算法狀態的整數。

0：名義上進行優化。

1：達到迭代限製。

2：問題似乎不可行。

3：問題似乎是無限的。

4：遇到數字困難。

尼特 int

當前迭代次數。

信息 str

算法狀態的字符串說明符。

HiGHS 方法當前不支持回調函數。

options： 字典，可選

求解器選項字典。所有方法都接受以下選項：

馬克西特 int

要執行的最大迭代次數。默認值：參見method-specific 文檔。

顯示 bool

設置為True 以打印收斂消息。默認值：False。

預解決 bool

設置為 False 以禁用自動預求解。默認值：True。

除了 HiGHS 求解器之外的所有方法也接受：

tol 浮點數

一個容差，它確定殘差何時為 “close enough” 到零被視為完全為零。

自動縮放 bool

設置為True 以自動執行平衡。如果約束中的數值相隔幾個數量級，請考慮使用此選項。默認值：False。

rr bool

設置為False 以禁用自動冗餘刪除。默認值：True。

rr_method string

用於在預求解後從等式約束矩陣中識別和刪除冗餘行的方法。對於密集輸入的問題，可用的冗餘去除方法有：

“SVD”：

對矩陣重複執行奇異值分解，根據左奇異向量中與零奇異值對應的非零值檢測冗餘行。當矩陣接近滿秩時可能會很快。

“pivot”：

使用 [5] 中提出的算法來識別冗餘行。

“ID”：

使用隨機插值分解。標識未在矩陣的全秩插值分解中使用的矩陣轉置列。

None:

如果矩陣接近滿秩，即矩陣秩與行數之差小於五，則使用“svd”。如果沒有，請使用“pivot”。此默認行為如有更改，恕不另行通知。

默認值：無。對於稀疏輸入的問題，忽略此選項，並使用[5]中提出的基於主元的算法。

對於 method-specific 選項，請參閱 show_options('linprog') 。

x0： 一維數組，可選

猜測決策變量的值，這些值將通過優化算法進行細化。此參數目前僅由“修訂單純形”方法使用，並且隻能在 x0 表示基本可行解時使用。

integrality： 一維數組或整數，可選

指示每個決策變量的完整性約束類型。

0：連續變量；沒有完整性約束。

1：整數變量；決策變量必須是整數界限.

2：Semi-continuous變量；決策變量必須在界限或取值0.

3：Semi-integer變量；決策變量必須是整數界限或取值0.

默認情況下，所有變量都是連續的。

對於混合完整性約束，提供形狀為 c.shape 的數組。為了從較短的輸入推斷每個決策變量的約束，參數將使用 np.broadcast_to 廣播到 c.shape。

該參數當前僅由 'highs' 方法使用，否則將被忽略。

返回 ：：

res： OptimizeResult

A scipy.optimize.OptimizeResult由以下字段組成。請注意，字段的返回類型可能取決於優化是否成功，因此建議檢查OptimizeResult.status在依賴其他字段之前：

x 一維數組

在滿足約束條件的同時最小化目標函數的決策變量的值。

樂趣 浮點數

目標函數c @ x的最優值。

鬆弛 一維數組

鬆弛變量的(名義上為正)值 b_ub - A_ub @ x 。

騙局 一維數組

等式約束的(名義上為零)殘差 b_eq - A_eq @ x 。

成功 bool

True 當算法成功找到最優解時。

狀態 int

一個整數，表示算法的退出狀態。

0：優化成功終止。

1：達到迭代限製。

2：問題似乎不可行。

3：問題似乎是無限的。

4：遇到數字困難。

尼特 int

在所有階段執行的迭代總數。

信息 str

算法退出狀態的字符串說明符。

注意：

本節介紹可通過 ‘method’ 參數選擇的可用求解器。

‘highs-ds’和‘highs-ipm’是HiGHS simplex 和interior-point 方法求解器的接口[13]， 分別。‘highs’(默認)自動在兩者之間進行選擇。這些是 SciPy 中最快的線性規劃求解器，特別是對於大型稀疏問題；這兩個中哪個更快是problem-dependent。其他求解器 (‘interior-point’,‘修訂單工’， 和‘simplex’) 是遺留方法，將在 SciPy 1.11.0 中刪除。

方法highs-ds是 C++ 高性能對偶修正單純形實現 (HSOL) 的包裝器[13],[14].方法highs-ipm是一個 C++ 實現的包裝器i內部-p軟膏m方法[13];它具有交叉例程，因此與單純形求解器一樣準確。方法高點自動在兩者之間進行選擇。對於涉及的新代碼linprog，我們建議明確選擇這三個方法值之一。

方法interior-point使用 primal-dual 路徑跟隨算法，如[4].該算法支持稀疏約束矩陣，通常比單純形法更快，尤其是對於大型稀疏問題。但是請注意，返回的解可能比單純形方法的精度稍差，並且通常不會與約束定義的多麵體的頂點相對應。

方法修正單純形使用修改後的單純形法，如[9], 除了分解[11]在算法的每次迭代中，有效地維護和用於求解線性係統的基本矩陣，而不是其逆矩陣。

方法單純形使用 Dantzig 單純形算法的傳統 full-tableau 實現[1],[2](不是Nelder-Mead 單純形)。包含此算法是為了向後兼容和教育目的。

申請前interior-point,修正單純形， 或者單純形, 一個基於預求解過程[8]嘗試識別瑣碎的不可行性、瑣碎的無界性和潛在的問題簡化。具體來說，它檢查：

• A_eq 或 A_ub 中的零行，表示微不足道的約束；

• 中的零列A_eq 和 A_ub，代表無約束變量；

• A_eq 中的列單例，表示固定變量；和

• A_ub 中的列單例，表示簡單邊界。

如果 presolve 顯示問題是無界的(例如，無約束和無界的變量具有負成本)或不可行(例如，A_eq 中的一行零對應於 b_eq 中的非零)，求解器以適當的狀態代碼終止.請注意，一旦檢測到任何無界跡象，presolve 就會終止；因此，當問題實際上是不可行的(但尚未檢測到不可行)時，問題可能會被報告為無界的。因此，如果知道問題是否實際上不可行很重要，請使用選項 presolve=False 再次解決問題。

如果在單次預求解中既沒有檢測到不可行性也沒有檢測到無界，則在可能的情況下收緊邊界，並從問題中刪除固定變量。然後，刪除A_eq 矩陣的線性相關行(除非它們表示不可行)以避免主求解例程中的數值困難。請注意，也可能會刪除幾乎線性相關的行(在規定的容差內)，這可能會在極少數情況下改變最佳解決方案。如果這是一個問題，請消除問題公式中的冗餘並使用選項 rr=False 或 presolve=False 運行。

這裏可以進行一些潛在的改進：應該實施[8]中概述的額外預求解檢查，預求解例程應該運行多次(直到無法進行進一步的簡化)，並且應該從[5]中獲得更多的效率改進在冗餘去除例程中實現。

在預求解之後，通過將(緊縮的)簡單邊界轉換為上限約束，為不等式約束引入非負鬆弛變量，並將無界變量表示為兩個非負變量之間的差異，將問題轉換為標準形式。可選地，問題通過平衡[12]自動縮放。所選算法解決標準形式問題，後處理例程將結果轉換為原始問題的解決方案。

參考：

[1]

Dantzig, George B.，線性規劃和擴展。蘭德公司研究研究普林斯頓大學。新聞，普林斯頓，新澤西州，1963

[2]

希利爾，S.H.和利伯曼，G.J. (1995)，“數學編程導論”，McGraw-Hill，第 4 章。

[3]

Bland, Robert G. 單純形法的新有限旋轉規則。運籌學數學 (2)，1977：第 103-107 頁。

[4]

Andersen、Erling D. 和 Knud D. Andersen。 “線性規劃的 MOSEK 內點優化器：齊次算法的實現。”高性能優化。施普林格美國，2000。197-232。

[5] (1,2,3)

Andersen, Erling D. “在 large-scale 線性規劃中查找所有線性相關的行。”優化方法和軟件 6.3 (1995): 219-227。

[6]

Freund, Robert M. “Primal-Dual Interior-Point 基於牛頓法的線性規劃方法。”未發表的課程筆記，2004 年 3 月。可於 2017 年 2 月 25 日在https://ocw.mit.edu/courses/sloan-school-of-management/15-084j-nonlinear-programming-spring-2004/lecture-notes/lec14_int_pt_mthd.pdf

[7]

福勒，羅伯特。 “通過Interior-Point 方法求解線性規劃。”未發表的課程筆記，2005 年 8 月 26 日。可於 2017 年 2 月 25 日在 http://www.4er.org/CourseNotes/Book%20B/B-III.pdf

[8] (1,2)

Andersen、Erling D. 和 Knud D. Andersen。 “在線性規劃中進行預求解。”數學編程 71.2(1995)：221-245。

[9]

Bertsimas、Dimitris 和 J. Tsitsiklis。 “線性規劃導論。”雅典娜科學 1 (1997): 997。

[10]

安徒生，Erling D.，等人。大規模線性規劃的內點方法的實現。 HEC/日內瓦大學，1996 年。

[11]

Bartels, Richard H. “單純形法的穩定性”。數值數學雜誌 16.5 (1971): 414-434。

[12]

Tomlin, J. A. “關於縮放線性規劃問題”。數學規劃研究 4(1975)：146-166。

[13] (1,2,3)

Huangfu, Q.、Galabova, I.、Feldmeier, M. 和 Hall, J. A. J.“HiGHS - 用於線性優化的高性能軟件。” https://highs.dev/

[14]

Huangfu, Q. 和 Hall, J. A. J. “並行化對偶修正單純形法”。數學規劃計算, 10 (1), 119-142, 2018. DOI: 10.1007/s12532-017-0130-5

例子：

考慮以下問題：

\[\begin{split}\min_{x_0, x_1} \ -x_0 + 4x_1 & \\ \mbox{such that} \ -3x_0 + x_1 & \leq 6,\\ -x_0 - 2x_1 & \geq -4,\\ x_1 & \geq -3.\end{split}\]

該問題並未以 linprog 接受的形式呈現。通過將 “greater than” 不等式約束轉換為 “less than” 不等式約束(將兩邊乘以因子 \(-1\) )可以輕鬆解決此問題。另請注意，最後一個約束實際上是簡單邊界 \(-3 \leq x_1 \leq \infty\) 。最後，由於 \(x_0\) 沒有界限，我們必須顯式指定界限 \(-\infty \leq x_0 \leq \infty\) ，因為默認情況下變量為非負數。將係數收集到數組和元組中後，該問題的輸入為：

>>> from scipy.optimize import linprog
>>> c = [-1, 4]
>>> A = [[-3, 1], [1, 2]]
>>> b = [6, 4]
>>> x0_bounds = (None, None)
>>> x1_bounds = (-3, None)
>>> res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])
>>> res.fun
-22.0
>>> res.x
array([10., -3.])
>>> res.message
'Optimization terminated successfully. (HiGHS Status 7: Optimal)'

邊際值(又稱對偶值/影子價格/拉格朗日乘數)和殘差(鬆弛)也可用。

>>> res.ineqlin
  residual: [ 3.900e+01  0.000e+00]
 marginals: [-0.000e+00 -1.000e+00]

例如，由於與第二個不等式約束相關的邊際為 -1，因此如果我們在第二個不等式約束的右側添加少量 eps，我們預計目標函數的最優值將減少 eps :

>>> eps = 0.05
>>> b[1] += eps
>>> linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds]).fun
-22.05

另外，由於第一個不等式約束的殘差為 39，因此我們可以將第一個約束的右側減少 39，而不會影響最優解。

>>> b = [6, 4]  # reset to original values
>>> b[0] -= 39
>>> linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds]).fun
-22.0