使用SciPy.optimize.linprog函数解决线性规划问题指南_scipy linprog

1、线性规划问题定义

linprog函数可以求解线性等式约束和不等式约束下的最小目标值问题。这类问题就是众所周知的线性规划。线性规划可以解决以下形式的问题：
$$\begin{gathered} \underset{x}{\min }{c}^{T}x \\ \mathrm{such}\mathrm{that}{A}_{ub}x\le b_{ub}, \\ {A}_{eq}x=b_{eq}, \\ l\le x\le u, \end{gathered}$$
其中$x$是一个决策变量向量；$c,b_{ub},b_{eq},l,u$是向量；$A_{ub},A_{eq}$是矩阵。在本教程中，我们将尝试使用 linprog 解决一个典型的线性规划问题。

2、线性规划问题示例

考虑下面这个简单的线性规划问题：
$$\begin{gathered} \underset{x_1,x_2,x_3,x_4}{\max }29x_1+45x_2 \\ \mathrm{such}\mathrm{that}{x}_1-x_2-3x_3\le 5 \\ 2x_1-3x_2-7x_3+3x_4\ge 10 \\ 2x_1+8x_2+x_3=60 \\ 4x_1+4x_2+x_4=60 \\ 0\le x_0 \\ 0\le x_1\le 5 \\ {x}_2\le 0.5 \\ -3\le x_3 \end{gathered}$$
我们需要进行一些数学操作，将目标问题转换成 linprog 可以接受的形式。

首先，我们来看看目标函数。我们希望最大化目标函数，但 linprog 只能接受最小化问题。将最大化$29x_1+45x_2$转换为最小化$-29x_1-45x_2$，就很容易解决这个问题。另外，目标函数中没有显示$x_3,x_4$。这意味着与$x_3,x_4$对应的权重为零。因此，目标函数可以转换为：
$$\underset{x_1,x_2,x_3,x_4}{\min }-29x_1-45x_2+0x_3+0x_4$$

如果我们定义决策变量$ x = [x_{1},x_{2},x_{3},x_{4}]^{T} $的向量，那么在这个问题中，linprog的目标权重向量$c$应该是
$$c=[-29,-45,0,0{]}^T$$

接下来，我们来看看两个不等式约束。第一个不等式是“小于”不等式，所以它已经是 linprog 接受的形式。第二个不等式是“大于”不等式，因此我们需要将两边乘以-1，将其转换为“小于”不等式。明确显示零系数后，我们得到：

$$\begin{gathered} x_1-x_2-3x_3+0x_4\le 5 \\ -2x_1+3x_2+7x_3-3x_4\le -10 \end{gathered}$$

这些方程可以转换成矩阵形式：

$$A_{ub}x\le b_{ub}$$

其中

A u b = [ 1 − 1 − 3 0 − 2 3 7 − 3 ] b u b = [ 5 − 10 ] A_{ub} =

$$\begin{bmatrix} 1 & -1 & -3 & 0 \\ -2 & 3 & 7 & -3 \\ \end{bmatrix}$$ \\ b_{ub} = $$\begin{bmatrix} 5 \\ -10 \\ \end{bmatrix}$$ Aub​=[1−2​−13​−37​0−3​]bub​=[5−10​]

接下来，我们来看看两个相等约束条件。明确显示权重为零的是：

$$\begin{gathered} 2x_1+8x_2+1x_3+0x_4=60 \\ 4x_1+4x_2+0x_3+1x_4=60 \end{gathered}$$

这些方程可以转换成矩阵形式：

$$A_{eq}x=b_{eq}$$

其中

A e q = [ 2 8 1 0 4 4 0 1 ] b e q = [ 60 60 ] A_{eq} =

$$\begin{bmatrix} 2 & 8 & 1 & 0 \\ 4 & 4 & 0 & 1 \\ \end{bmatrix}$$ \\ b_{eq} = $$\begin{bmatrix} 60 \\ 60 \\ \end{bmatrix}$$ Aeq​=[24​84​10​01​]beq​=[6060​]

最后，让我们考虑对单个决策变量的单独不等式约束，即所谓的“边界约束”或“简单约束”。这些约束可以通过linprog的界限参数来应用。正如linprog文档中所指出的，边界的默认值是0或无，这意味着每个决策变量的下限都是0，而每个决策变量的上限都是无穷大：所有的决策变量都是非负的。我们的界限是不同的，因此需要将每个决策变量的下限和上限指定为一个元组，并将这些元组组合成一个列表。
最后，我们可以使用 linprog 解决转换后的问题。

import numpy as np
from scipy.optimize import linprog
c = np.array([-29.0, -45.0, 0.0, 0.0])
A_ub = np.array([[1.0, -1.0, -3.0, 0.0],
                [-2.0, 3.0, 7.0, -3.0]])
b_ub = np.array([5.0, -10.0])
A_eq = np.array([[2.0, 8.0, 1.0, 0.0],
                [4.0, 4.0, 0.0, 1.0]])
b_eq = np.array([60.0, 60.0])
x0_bounds = (0, None)
x1_bounds = (0, 5.0)
x2_bounds = (-np.inf, 0.5)  # +/- np.inf can be used instead of None
x3_bounds = (-3.0, None)
bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
print(result.message)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

The problem is infeasible. (HiGHS Status 8: model_status is Infeasible; primal_status is At lower/fixed bound)
1

结果表明，我们的问题是不可行的，即没有满足所有约束条件的解向量。这并不一定意味着我们做错了什么，有些问题确实是不可行的。
然而，假设我们决定对$x_1$的边界约束太紧，可以放宽到$0\le x_1\le 6$。接下来调整我们的代码 x1_bounds = (0,6）以反应更改并再次执行它：

x1_bounds = (0, 6)
bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
print(result.message)
1
2
3
4

Optimization terminated successfully. (HiGHS Status 7: Optimal)
1

结果显示优化成功。我们可以检查目标值（result.fun）是否与$c^{T}x$相同：

x = np.array(result.x)
obj = result.fun
print(c @ x)
1
2
3

-505.974358974359
1

print(obj)
1

-505.974358974359
1

我们还可以检查所有约束条件是否在合理的公差范围内得到满足：

print(b_ub - (A_ub @ x).flatten())  # this is equivalent to result.slack
1

[ 0.00000000e+00 -7.10542736e-15]
1

print(b_eq - (A_eq @ x).flatten())  # this is equivalent to result.con
1

[-7.10542736e-15  0.00000000e+00]
1

print([0 <= result.x[0], 0 <= result.x[1] <= 6.0, result.x[2] <= 0.5, -3.0 <= result.x[3]])
1

[True, True, True, True]
1