使用Python和SciPy解决线性规划问题

发表于 更新于 分类于

Python是一种很流行的编程语言，由于其语法简单而功能强大，已被广泛使用。目前的最新版本是3.10。

SciPy是一个Python的软件包，其基于NumPy软件包扩展了大量功能，包括统计方法、信号处理、线性代数以及规划问题。最新版本为1.9.1。

本文将基于Python这一编程语言以及Scipy这一软件包，说明如何使用这二者进行线性规划问题的求解。

0、准备（Python和SciPy的安装）

(1) Python的安装

有多种方法可以安装Python，这里列举两个方法。

  1. 从Python官方网站安装。点击网页上的Download按钮，下载安装包，然后直接执行安装即可。（比较麻烦，不推荐）
  2. 安装Anaconda。Anaconda是一个优秀的Python发行版，内置了大量科学计算和数据分析的软件包，并带有conda这一包管理器和Jupyter这一IDE，是数据分析的首选。从官方网站下载安装包，然后直接执行安装即可。（推荐）

(2) SciPy的安装

如果已经安装了Anaconda，则无需进行这一步，Anaconda已经打包了SciPy在自己的发行版中。

如果安装的是Python，可以使用下面的方法安装。首先，确定Python和pip已经被加载到环境变量中；然后，打开一个命令行窗口（Windows的“命令提示符”和macOS的“终端”），输入下面的指令：

1
2
pip install scipy  -i https://pypi.tuna.tsinghua.edu.cn/simple
# 其中，-i 参数指明使用清华大学软件源进行安装，以提高下载速度

即可执行安装。

1、快速上手

SciPy提供的线性规划计算工具为scipy.optimize.linprog函数。因此要想计算线性规划问题，需要按照要求导入scipy这一软件包

1
2
import numpy as np
from scipy.optimize import linprog

与我们在教科书中学习到的形式略有不同，SciPy中的线性规划问题形式如下：

$$
\begin{align}
min\ z =c^{T}x ; \\
s.t.\ ,\\
A_{ub}x \le b_{ub}\ , \\
A_{eq}x = b_{eq}\ , \\
l \le x \le u
\end{align}
$$

其中，$x$是待求解的决策变量向量；$c,\ b_{ub},\ b_{eq},\ l,\ u$均为向量，而$A_{ub},\ A_{eq}$均为矩阵。在SciPy中，线性规划的目标是最小化目标函数的值，这与我们在教材上学到的最大化刚好相反，因此在构造线性规划问题的参数向量时，目标函数的系数应该取实际问题的相反数。

以一个实际问题为例（参考SciPy官方文档）：

$$
\begin{align}
max\ z=29x_{1}+45x_{2}; \\
s.t.\ , \\
x_1-x_2-3x_3 \le 5 \ ,\\
2x_2-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_1 \ ,\\
0 \le x_2 \le 5 \ ,\\
x_3 \le 0.5 \ ,\\
-3 \le x_4
\end{align}
$$

在上面这个问题中，我们的目标是最大化目标函数$z$，因此为了使用SciPy求解，我们必须对目标函数的系数向量取相反数

1
c = np.array([-29.0, -45.0, 0.0, 0.0])

在约束条件中，有几个不等式，我们将其全部化成小于等于的形式，得到约束条件的系数矩阵A_ub和列向量b_ub

1
2
3
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和列向量b_eq

1
2
3
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])

此外，我们注意到，变量 $x_1$ 到 $x_4$ 还有取值范围的限制。因此这一项也应该被我们考虑进去。我们使用一个一维数组存储 $x_1$ 到 $x_4$ 的取值范围，其中每个变量的取值范围用一个元组（Tuple）进行存储

1
2
3
4
5
x1_bounds = (0, None)       # None代表可以取到无穷大
x2_bounds = (0, 5.0)
x3_bounds = (-np.inf, 0.5)  # +/- np.inf 可以代替None用以表示无穷大或者无穷小
x4_bounds = (-3.0, None)
bounds = [x1_bounds, x2_bounds, x3_bounds, x4_bounds]  # 这个数组存储了变量x1到x4的取值范围。

之后，用一行代码即可完成线性规划的求解

1
2
3
4
# 别忘了   import numpy as np
# 也别忘了 from scipy.optimize import linprog
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
print(result) # 输出求解结果

SciPy给出的结果如下所示：

1
2
3
4
5
6
7
8
9
    con: array([15.53611731, 16.61287584])
    fun: -370.23223964124924
message: 'The algorithm terminated successfully and determined that the problem
          is infeasible.'
    nit: 6
  slack: array([ 0.79315063, -1.76308707])
 status: 2
success: False
      x: array([ 6.60059411,  3.9736669 , -0.52664072,  1.09008012])

看到那个success: False了吗？这意味着我们的问题没有可行解。这并不一定意味着我们做错了，因为一些问题确实不存在可行解。但是，也许我们将条件放宽一点，例如令 0 <= x_2 <= 6 ，就能得到可行解。

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

SciPy给出的结果如下所示：

1
2
3
4
5
6
7
8
    con: array([9.79743930e-09, 1.04769242e-08])
    fun: -505.97435889005345
message: 'Optimization terminated successfully.'
    nit: 4
  slack: array([ 6.53052723e-10, -2.26972219e-09])
 status: 0
success: True
      x: array([ 9.41025641,  5.17948718, -0.25641026,  1.64102564])

这回确实得到了可行解。我们看一下目标函数的值是多少：

1
2
3
x = np.array(result.x)
print(c @ x) # @运算符是NumPy软件包中的点乘运算符
# out: -505.97435889005345

2、函数原型与参数含义

我们现在研究一下scipy.optimize.linprog这个函数，以便学习其用法。

参考SciPy文档，这一函数的原型如下：

1
2
3
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)

接下来，我们逐一分析相关参数。首先是与线性规划问题紧密相关的5个参数：

  • c：一个一维数组（向量），代表目标函数的系数向量

  • A_ub：一个二维数组（矩阵），是不等式约束条件的系数矩阵

  • b_ub：一个一维数组（向量），是不等式约束条件的资源向量

  • A_eq：一个二维数组（矩阵），是等式约束条件的系数矩阵

  • b_eq：一个一维数组（向量），是等式约束条件的资源向量

  • bounds：一个元组类型的一维数组，每一个元组都定义了一个决策变量 $x_i$ 的取值范围(min, max)，可以在元组中使用None代表正负无穷（例如，(None, 0)代表取值范围为负无穷大到0）。这是一个可选参数，如果不给出，则默认值为(0, None)，即0到正无穷大。

  • method：字符串，指定求解方法，默认值为'highs'，可选择的其他方法包括'highs-ds'，'highs-ipm'，'interior-point'，'revised simplex'以及'simplex'。这一参数将在SciPy的1.11.0以后的版本中废除。

  • x0：一个一维数组（向量），代表一个初始基可行解。这个参数目前只在'revised simplex'方法中被用到，因此可以不给出。

  • integrality：一个一维数组，指定每个决策变量的取整约束条件，默认是不存在取整约束。这个参数仅在highs方法中被用到。对每一个决策变量，不同的取整约束如下：

    • 0：连续变量，可以取到取值范围内的所有实数，不存在取整约束。
    • 1：整数变量，只能取到取值范围内的所有整数。
    • 2：半连续变量，可以取到取值范围内的所有实数，以及0。
    • 3：半整数变量，只能取到取值范围内的所有整数，以及0。

此外，还有两个可选参数callback和options，这两个参数包含一些与程序运行以及返回值有关的选项，而与线性规划问题本身无关，因此本文略过。

3、返回值

线性规划求解函数scipy.optimize.linprog的返回值是一个scipy.optimize.optimize.OptimizeResult类型的对象，其包含下列属性：

  • x：决策变量的值。如果问题存在可行解的话，x就是一个可行解。
  • fun：目标函数的值，也就是c@x（向量CCC点乘向量XXX）的结果。
  • slack：不等式约束中，松弛变量的值，原则上是正的。其等于b_ub-A_ub@x。
  • con：等式约束中，残差的值，原则上等于0。其等于b_eq-A_eq@x
  • success：一个布尔变量，代表是否找到了最优解。
  • status：状态码。含义如下：
    • 0：成功找到最优解
    • 1：达到了迭代次数的最大限制
    • 2：问题似乎无解
    • 3：问题似乎存在无界解
    • 4：遇到了数值困难
  • nit：迭代次数
  • message：对问题求解结果的文字描述。

4、例题

某工厂在计划期内要安排生产I、II两种产品，已知生产单位产品所需设备台时及原材料消耗如下表所示

资源 产品I 产品II 现有条件
设备 1台时/件 2台时/件 8台时
原材料A 4kg/件 0 16kg
原材料B 0 4kg/件 12kg

该工厂每生产一件产品I可获利2元，每生产一件产品II可获利3元，如何安排生产计划使得该工厂获利最多？

解：

我们建立上述问题的线性规划数学模型，编写程序如下：

1
2
3
4
5
6
7
8
9
import numpy as np
from scipy.optimize import linprog
c = np.array([-2.0, -3.0])
A_ub = np.array([[1.0, 2.0],
                [4.0, 0.0],
                [0.0, 4.0]])
b_ub = np.array([8.0, 16.0, 12.0])
result = linprog(c, A_ub=A_ub, b_ub=b_ub)
print(result)

程序的运行结果如下所示：

1
2
3
4
5
6
7
8
    con: array([], dtype=float64)
    fun: -13.999999982532964
message: 'Optimization terminated successfully.'
    nit: 4
  slack: array([1.02599227e-08, 1.66172107e-08, 4.00000001e+00])
 status: 0
success: True
      x: array([4., 2.])

我们可以看到，SciPy成功帮助我们求出了最优解，即x1=4,x2=2，此时目标函数值 z=−z1=13.999≈14 ，也就是说，原问题的最大利润值是14元。