线性规划

线性规划（Linear Programming, LP）是最简单的优化问题类型，指目标函数与约束条件均是关于变量的线性函数的优化问题。

建立线性规划模型

对于任何线性规划问题，我们通过下列步骤建立模型：

定义问题变量
给出目标函数
提供约束条件

考虑以下工厂生产问题，我们将通过它了解线性规划的建模过程。

工厂现有的某机器能生产甲、乙、丙三种产品，它们的具体生产信息如下表所示：

每种产品的单件利润
每种产品在生产单件时所需的人工工时
每种产品在生产单件时所需的机器工时
每种产品的数量

产品	利润	人工工时	机器工时	产品数量
A	10	2	4
B	7	0.5	3
C	3	1	1
约束		= 45	<= 100	<= 60

假设该机器平均每周固定分配人工 45 小时、机器工作不超过 100 小时，且不应生产超过 60 件产品；生产的产品全部可以卖出。求如何分配三种产品的生产数量（允许是小数），使得平均每周的总利润最大。

定义问题变量

该问题需要求解即每种产品的生产数量。因此，我们将产品 A、B、C 的生产数量作为变量，分别定义为 \(x_1\), \(x_2\) 与 \(x_3\).
给出目标函数

该问题的目标是最大化总利润，而每种产品的单位利润已经给出。将生产产品的数量与其利润相乘，最后求和，就得到总利润：

\[\max \, 10x_1 + 7x_2 + 3x_3\]

为了统一形式，本节将一律采用最小化问题的形式。我们可以将总利润乘以 \(-1\) 并求其最小值（记为最小化目标函数 \(z\)）：

\[\min \, z = -10x_1 - 7x_2 - 3x_3\]
提供约束条件

根据生产信息表格，我们用不等式约束机器工时、产品数量，用等式约束人工工时：
- 机器工时约束：\(4x_1 + 3x_2 + x_3 \leq 100\)
- 产品数量约束：\(x_1 + x_2 + x_3 \leq 60\)
- 人工工时约束：\(2x_1 + 0.5x_2 + x_3 = 45\)
最后，每种产品的生产数量不能为负数：
- 非负约束：\(x_1, x_2, x_3 \geq 0\)

总结以上步骤，我们列出优化模型如下。由于目标函数与所有约束条件都是线性的，因此它是一个线性规划问题：

\begin{align*} \min\ & z(x_1, x_2, x_3) = -10x_1 - 7x_2 - 3x_3 \\ \text{s.t.}\quad & 4x_1 + 3x_2 + x_3 \leq 100 \\ & x_1 + x_2 + x_3 \leq 60 \\ & 2x_1 + 0.5x_2 + x_3 = 45 \\ & x_1, x_2, x_3 \geq 0 \end{align*}

或者其矩阵形式：

\begin{align*} \min\ & z(x_1, x_2, x_3) = \begin{pmatrix} -10 \\ -7 \\ -3 \end{pmatrix}^T \begin{pmatrix} x_1 \\ x_2 \\ x_3\end{pmatrix} \\ \text{s.t.}\quad & \begin{pmatrix} 4 & 3 & 1 \\ 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3\end{pmatrix} \leq \begin{pmatrix} 100 \\ 60 \end{pmatrix} \\ & \begin{pmatrix} 2 & 0.5 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = 45\\ & \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} \geq \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \end{align*}

该矩阵形式可以进一步变形为一般形式，以减少矩阵的数量。这里我们将约束条件分为小于等于约束、等于约束与变量上下限约束三个部分，仅仅是便于后续的代码使用。

线性规划的一般形式

通过上述例题的模型建立，我们总结：对于拥有 \(n\) 个变量、\(m\) 个约束的线性规划，其模型的一般形式如下

\begin{align*} \min \quad & z(x) = c_1 x_1 + c_2 x_2 + \dotsm + c_n x_n + d \\ \text{s.t.} \quad & a_{11}x_1 + a_{12}x_2 + \dotsm + a_{1n}x_n \leq b_1 \\ & a_{21}x_1 + a_{22}x_2 + \dotsm + a_{2n}x_n \leq b_2 \\ & \vdots \\ & a_{m1}x_1 + a_{m2}x_2 + \dotsm + a_{mn}x_n \leq b_m \end{align*}

或者，我们使用其矩阵标准形式：

\begin{align*} \min \quad & z(x) = c^T x \\ \text{s.t.} \quad & \boldsymbol{A}x \leq b \end{align*}

其中变量 \(x \in \mathbb{R}^n\) 是 n 维向量，\(c \in \mathbb{R}^n\) 描述了组成目标函数的参数，\(\boldsymbol{A} \in \mathbb{R^{m\times n}}\) 与 \(b \in \mathbb{R^m}\) 描述了关于变量 \(x\) 的 m 个不等式约束条件。

对于一个不符合标准形式的线性规划模型，我们可以通过以下技巧（在引言中已有所介绍），将其等价变形为标准形式：

目标函数中的常数项：由于目标函数中的常数 \(d \in \mathbb{R}\) 的取值不影响最优解 \(x^*\) 的求解结果，因此我们可以将它省略。
最小化与最大化目标函数：对于求函数 \(z(x)\) 最大值的目标函数，我们可以将其转为求 \(-z(x)\) 的最小值。

\[\max\,z(x) \iff \min\,-z(x)\]
大于等于形式的约束：对于大于等于形式的不等式约束 \(\boldsymbol{A}_{geq} x \geq b_{geq}\)，我们可以将两边乘以 \(-1\) 来将其转变为小于形式的约束：

\[\boldsymbol{A}_{geq} x \geq b_{geq} \iff -\boldsymbol{A}_{geq} x \leq -b_{geq}\]
等于形式的约束：对于等于形式的约束 \(\boldsymbol{A}_{eq} x = b_{eq}\)，我们可以用一个大于等于约束与一个小于等于约束来描述，然后进一步转为两个不等式描述：

\[\begin{split}\boldsymbol{A}_{eq} x = b_{eq} \iff \begin{cases} \boldsymbol{A}_{eq} x \leq b_{eq} \\ \boldsymbol{A}_{eq} x \geq b_{eq} \end{cases} \iff \begin{cases} \boldsymbol{A}_{eq} x \leq b_{eq} \\ -\boldsymbol{A}_{eq} x \leq -b_{eq} \end{cases}\end{split}\]

线性规划的代码求解

关于线性规划的代码求解，我们有很多工具可以使用。本节将会使用：

[1]:

import scipy
import cvxpy

for lib in (scipy, cvxpy):
    print(lib.__name__, lib.__version__, sep='\t')

scipy   1.14.1
cvxpy   1.6.4

scipy.optimize.linprog

Scipy 可能是许多人第一次在 Python 上解线性规划问题时使用的工具。scipy.optimize.linprog 允许我们使用在矩阵形式中使用小于等于约束、等式约束，以及变量上下限约束：

\begin{align*} \min\ & z(x) = c^T x \\ \text{s.t.}\quad & \boldsymbol{A}_{leq} x \leq b_{leq} \\ & \boldsymbol{A}_{eq} x = b_{eq} \\ & \ell \leq x \leq u \end{align*}

[2]:

import numpy as np
from scipy.optimize import linprog

c = np.array([-10, -7, -3])
A_leq = np.array([[4, 3, 1],
                  [1, 1, 1]])
b_leq = np.array([100, 60])
A_eq = np.array([[2, 0.5, 1]])
b_eq = np.array([45])
x_bounds = [(0, None) for _ in range(3)]

result = linprog(c, A_leq, b_leq, A_eq, b_eq, bounds=x_bounds)
print(result.message)

Optimization terminated successfully. (HiGHS Status 7: Optimal)

以上 Optimal 状态信息表示 LP 问题成功解得最优值。我们利用 result.x 查看取得最优值时，变量的取值：

[3]:

x = np.array(result.x)
x

[3]:

array([ 0., 22., 34.])

查看目标函数最优值 result.fun ，或者使用 c @ x 手动计算：

[4]:

print("Optimal:\t", result.fun)
print("Manual compute:\t", c @ x)

Optimal:         -256.0
Manual compute:  -256.0

还可以用 result.slack 查看不等式约束的满足情况，或用 result.con 查看等式约束的满足情况：

[5]:

print(f"Leq constraints: {b_leq - (A_leq @ x).flatten()}",
      f"\t{result.slack = }",
      f"Eq constraints: {b_eq - (A_eq @ x).flatten()}",
      f"\t{result.con = }",
      sep='\n')

Leq constraints: [0. 4.]
        result.slack = array([0., 4.])
Eq constraints: [0.]
        result.con = array([0.])

cvxpy

cvxpy 的用法也比较类似，不过需要用 cvxpy.Variable() 方法声明变量。

cvxpy 可以使用 <= 或者 >= 来表示带等的不等于号；但是，不能像 1<= x <=2 这样连写。

[6]:

import numpy as np
import cvxpy as cp

c = np.array([-10, -7, -3])
A_leq = np.array([[4, 3, 1],
                  [1, 1, 1]])
b_leq = np.array([100, 60])
A_eq = np.array([[2, 0.5, 1]])
b_eq = np.array([45])

x = cp.Variable(3)
x_lower = np.zeros(3)
problem = cp.Problem(cp.Minimize(c @ x),
                     [A_leq @ x <= b_leq, A_eq @ x <= b_eq, x >= x_lower])
problem.solve()

[6]:

np.float64(-255.99999999194907)

LP 问题的解结果（例如，求得最优 cvxpy.OPTIMAL 、不可行 cvxpy.INFEASIBLE 或者未约束 cvxpy.UNBOUNDED 等）：

[7]:

print(f"Status: {problem.status}")
print(problem.status == cp.OPTIMAL)

Status: optimal
True

LP 问题的解存放在 x.value 与 problem.value 中：

[8]:

print("A feasible solution:\t", x.value)
print("Optimal value:\t\t", problem.value)

A feasible solution:     [9.27094654e-09 2.20000000e+01 3.40000000e+01]
Optimal value:           -255.99999999194907

对偶解（Dual Solution）中符号为正的第 i 个值表示在最优解时约束 \(a_i x \leq b_i\) 恰好取得等号。这也意味着改变 \(b_i\) 会引起最优解的改变。

[9]:

constraint_labels = ["A_leq", "A_eq", "x_bounds"]
print("Dual solution:")
for cgroup, label in zip(problem.constraints, constraint_labels):
     print(f"    {label:<10}", cgroup.dual_value)

Dual solution:
    A_leq      [2.20000000e+00 8.68986474e-10]
    A_eq       [0.8]
    x_bounds   [3.99999999e-01 8.81435962e-11 5.94961809e-11]

我们传入了三组约束。在上述对偶解中，小于号约束组 A_leq 的第二值为 0，表示原问题的 \(x_1 + x_2 + x_3 \leq 60\) 是松弛的（在最优解时未取得等号）。同理，变量约束组 x_bounds 的第二、三值也是松弛的，表示 \(x_2, x_3\) 的最优解值不在边界上（而显然 \(x_1^* = 0\) 在边界上）。