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