线性规划单纯形求解理论
线性规划（Linear Programming, LP）是优化理论中用于在给定约束条件下最大化或最小化线性目标函数的一种数学方法。线性规划的最优解总是出现在可行域的顶点上，这是因为目标函数在可行域内的变化是线性的，因此在顶点处函数的值可能达到极值（最大或最小）。求解线性规划问题的常用方法之一是单纯形法（Simplex Method），其主要思路是从一个可行解开始，通过沿着可行域边界在顶点之间移动，以找到一个最优解。每一步移动都使得目标函数的值得以改进，直到无法进一步改进为止，即达到最优解。
| 唯一最优解 | 无穷多最优解 |
|---|---|
 |
 |
一、线性规划原型
Wyndor Glass公司生产高质量的玻璃产品。该公司拥有 3 个工厂。铝框架和硬件在工厂 1 制造，木质框架在工厂 2 生产，玻璃生产和产品组装在工厂 3 完成。由于收人下滑, 高层管理者决定调整公司的生产线。不赢利的产品将被停止生产, 而将生产能力转移到有较大销售潜力的两种新产品。
产品 1：具有铝框架的 8 英尺玻璃门；产品 2: 具有双悬木质框架的 \(4 \times 6\) 英尺窗
产品 1 需要工厂 1 和工厂 3 的生产能力, 而不需要工厂 2 的生产能力。产品 2 需要工厂 2 和工厂 3 的生产能力。进行市场细分研究后得出结论: 公司能够销售工厂所能生产的全部产品。然而, 由于两种产品都为使用工厂 3 的生产能力而竞争, 不清楚如何确定两种产品的组合可以实现利润最大化, 因此组织了一个运筹小组来研究这个问题。
运筹小组首先与高层管理者讨论, 以明确该研究的管理目标。经过讨论之后明确了下面的问题。决定两种产品生产率的依据应该是总利润最大化。限制条件是三个工厂有限的生产能力(每种产品将以 20 个为一批进行生产, 因此生产率定义为每周生产该产品的批数)。满足这些限定条件的任何生产率的组合都是可行的, 包括一种产品产量为零, 其他产品尽量多的生产情况，见下面生产计划表。
| 工厂 | 产品 1 | 产品 2 | 每周可用的生产时间/小时 |
|---|---|---|---|
| 1 | 1 | 0 | 4 |
| 2 | 0 | 2 | 12 |
| 3 | 3 | 2 | 18 |
| 每批的利润/美元 | 3000 | 5000 |
目标是通过生产两种产品来最大化总利润。设：\(x_1\): 产品 1 的生产批数（具有铝框架的 8 英尺玻璃门）；\(x_2\): 产品 2 的生产批数（具有双悬木质框架的 \(4 \times 6\)英尺窗) ，该问题的线性规划模型为
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
# 定义线性规划模型的系数
c = [-3000, -5000] # 利润函数的系数，负号表示求最大值问题
A = [
[1, 0], # 工厂 1 约束
[0, 2], # 工厂 2 约束
[3, 2] # 工厂 3 约束
]
b = [4, 12, 18] # 约束条件右边的常数
# 使用 scipy.optimize.linprog 求解线性规划问题
result = linprog(c, A_ub=A, b_ub=b, method='highs')
# 定义x1的范围
x1 = np.linspace(0, 10, 400)
# 各个约束条件的线性函数
x2_3 = (18 - 3 * x1) / 2 # 工厂 3 约束
# 设置图形
plt.figure(figsize=(8, 6))
# 画出约束线
plt.axvline(x=4, color='blue', linestyle='--', label='Factory 1 Constraint: x1 <= 4') # 垂直线表示x1 <= 4
plt.axhline(y=6, color='green', linestyle='--', label='Factory 2 Constraint: x2 <= 6') # 水平线表示x2 <= 6
plt.plot(x1, x2_3, label='Factory 3 Constraint: 3x1 + 2x2 <= 18', color='orange')
# 填充可行域
y1 = np.maximum(0, np.minimum(6, x2_3))
plt.fill_between(x1, 0, y1, where=(x1 <= 4) & (y1 >= 0), color='grey', alpha=0.5)
# 标出最优解点
if result.success:
x_opt = result.x
plt.scatter(x_opt[0], x_opt[1], color='red', label=f'Optimal Solution: ({x_opt[0]:.2f}, {x_opt[1]:.2f})')
plt.annotate(f'({x_opt[0]:.2f}, {x_opt[1]:.2f})', (x_opt[0], x_opt[1]), textcoords="offset points", xytext=(-10, -10), ha='center', fontsize=10, color='red')
# 图形设置
plt.xlim(0, 10)
plt.ylim(0, 10)
plt.xlabel('Number of batches of Product 1')
plt.ylabel('Number of batches of Product 2')
plt.title('Feasible Region and Optimal Solution for Wyndor Glass Company')
plt.legend()
plt.grid()
plt.show()
二、线性规划的基本术语和概念
线性规划的标准型为：
其中
- \(X\)是一个\(n\)维变量向量，表示决策变量：\[X^T = (x_1,x_2, \ldots ,x_n ) \]
- \(C\)是一个\(n\)维系数向量，表示目标函数中的系数：\[C^T= (c_1 ,c_2,\ldots,c_n ) \]
- \(A\)是一个\(m \times n\)的矩阵，表示约束条件中的系数：\[A = \begin{bmatrix} a_{11} & a_{12} & \ldots & a_{1n} \\ a_{21} & a_{22} & \ldots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \ldots & a_{mn} \end{bmatrix} = \begin{bmatrix} P_{1} & P_{2} & \ldots & P_{n} \end{bmatrix} \]
| 术语 | 解释 |
|---|---|
| 可行解 (或容许解) | 满足(2.1) 约束条件的解 \(x=\left(x_1, x_2, \cdots, x_n\right)^T\)称为线性规划问题的可行解 |
| 可行域 | 所有可行解的集合称为可行解集或可行域 |
| 最优解 | 使得目标函数取到最小值的可行解称为线性规划问题的最优可行解，简称为最优解或者解 |
| 基 | 假设 \(A\) 是约束方程组的系数矩阵, 其秩数为 \(m, B\)是矩阵 \(\boldsymbol{A}\) 中由 \(\boldsymbol{m}\) 列构成的非奇异子矩阵, 则称 \(B\) 是线性规划问题的一组基 |
| 基本解 | 若令 (2.1) 中的非基变量 \(x_{m+1}=\cdots=x_n=0\), 得到解 \(x=\) \(\left(x_1, x_2, \cdots, x_m, 0, \cdots, 0\right)^T\) ，称 \(x\) 为基本解（这个解的非0分量的数目不大于方程个数 \(m\) ) |
| 基本可行解 | 满足非负约束条件的基本解称为基本可行解。基本可行解的非 0 分量的数目不大于 \(\boldsymbol{m}\) ，都是非负的 |
| 可行基 | 对应于基本可行解的基称为可行基 |
| 最优基 | 对应于基本最优解的基称为最优基 |
矩阵 \(\boldsymbol{B}\) 是由 \(\boldsymbol{A}\)的\(\boldsymbol{m}\) 个线性无关的列向量组成，不失一般性，可假设:
称 \(P_j(j=1,2, \cdots, m)\) 为基向量，与基向量 \(P_j\) 相对应的变量 \(x_j\) \((j=1,2, \cdots, m)\) 为基变量，否则称为非基变量。
这时线性规划(1)的一个基是 \(B=\left(P_1, P_2, \cdots, P_m\right)\) ，
设 \(x_B\) 是对应于这个基的基向量，即:
约束方程组 \(A x=b\) 基本解的数目至多是 \(C_n{ }^m\) 个。基本可行解是同时满足约束方程和变量非负约束的解。变量一般要求都是非负的，所以基本可行解的要求非 0 变量也是要非负的。
三、线性规划的求解理论及几何意义
3.1 求解理论
基本解是对于部分约束等式的解，几何意义为图上的一个顶点。
基本可行解是对于在满足所有约束等式的解，几何意义为图上的可行的顶点。
顶点：设 \(S\) 是凸集，\(x \in S\)，不存在 \(S\) 中的另外两个点 \(x_1\) 和 \(x_2\)，及 \(\lambda \in(0,1)\)，使 \(x=\lambda x_1+(1-\lambda) x_2\), 则称 x 是凸集 S 的极点，或称顶点。
定理1：设 \(A=\left(a_{i j}\right)_{m \times n}\), 其秩为 \(m\)且 \(m<n，x=\left(x_1, x_2, \ldots, x_n\right)^T\)，\(b=\left(b_1, b_2, \ldots, b_m\right)^T\)，则 \(x\) 为凸集 \(R=\{x \mid A x=b, x \geq 0\}\)一个顶点的充要条件是\(x\)为
的一个基本可行解。阐明了基本可行解与可行域的顶点之间一一对应的关系。
可行域这个凸多面体(单纯形)的任何一个顶点对应了矩阵A的一组基。
假设\(X^*\)是一个基本可行解，如果它是可行域\(D\)内的顶点，那么\(X^*\)不能表示为两个不同可行解\(X^1\)和\(X^2\)的凸组合：\[X^* = \lambda X^1 + (1 - \lambda) X^2, \quad 0 < \lambda < 1 \]若\(X^*\)是基本可行解，存在\(m\)个变量使得其为非零，其余为零（定义中提到的子矩阵\(B\)非奇异）。则\(X^1\)和\(X^2\)在相应位置上的变量也需为零，否则\(X^*\)无法表示为\(X^1\)和$X^2的凸组合。因此，基本可行解是可行域的顶点。
这样我们就把凸多面体(单纯形)可行域这一具象的概念和矩阵的一组基这一抽象概念建立起了联系，而后者则是利于编程实现的。只要能找到矩阵的一组基，就意味找找到了可行域中的一个顶点。
定理2：（线性规划问题的基本定理）
(1）若线性规划问题存在可行解，则一定存在基本可行解；
(2) 若线性规划问题存在最优解, 则一定存在最优的基本可行解。
线性规划问题若有最优解,必定在某顶点取到。
假设我们选择\(x_j\)作为入基变量，并通过步长\(\theta\)来控制它的变化。我们考虑新的基向量\(X_B'\)和目标函数值的变化。
目标函数在当前解时的值为：\[Z = C_B^T X_B \]入基后，新的目标函数值为：
\[Z' = C_B^T X_B' + c_j \theta \]我们希望通过调节\(\theta\)增加\(Z'\)值。
由于\(X_B' = X_B - \theta d\)，新的目标函数值可以表示为：\[Z' = C_B^T (X_B - \theta d) + c_j \theta\\ Z' = C_B^T X_B - \theta C_B^T d + c_j \theta\\ Z' = Z + \theta (c_j - C_B^T d)\]其中\(d = B^{-1}a_j\)：
\[Z' = Z + \theta (c_j - C_B^T B^{-1}a_j)\\ Z' = Z + \theta \sigma_j \]
- 当 \(\sigma_j > 0\)时，通过适当选择\(\theta >0\)，目标函数值\(Z'\)会增加。
- 选择合适的\(x_j\)入基（即使其检验数为正的非基变量入基）可以提升目标函数的值。
- 为使得目标函数值改变更大，选取\(\sigma_j > 0\)最大的那个非基变量进基。
- 为保持可行解的非负性，我们还需要确保新的基变量\(X_B'\geq 0\)，这要求：
\[\theta \leq \min \left\{ \frac{X_{B_i}}{d_i} \mid d_i > 0 \right\} \]通过选择合适的步长\(\theta\)，我们既能提升目标函数值又能保持解的可行性。
\(d = B^{-1}a_j\)的推导：
考虑约束条件\(A X = b\)，可以将其分为基变量和非基变量两部分：\[B X_B + N X_N = b \]其中 $ B $ 是当前基矩阵，$ N $ 是非基矩阵。
当我们选择 $ x_j $（属于 $ X_N $）作为入基变量时，我们有 $ x_j = \theta $（其中 $ \theta $ 是我们调节的步长）。
那么新的变量表示为：\[X_B' = X_B - \theta d \quad \text{和} \quad X_N' = X_N + \theta e_j \]其中 $ e_j $ 是一个单位向量，只有第 $ j $ 个位置为 1。代入约束条件中：
\[B (X_B - \theta d) + N (X_N + \theta e_j) = b \]这化简为：
\[B X_B - \theta B d + N X_N + \theta N e_j = b \]考虑到当前解 $ (X_B, X_N) $ 满足 $ B X_B + N X_N = b $，我们得到：
\[-\theta B d + \theta a_j = 0 \]其中 \(a_j\) 是 \(A\) 矩阵中对应于 \(x_j\) 的列向量（即 \(N e_j = a_j\)）。
我们得到：$ B d = a_j \quad d = B^{-1} a_j $
推论1: 若线性规划的可行域 \(R=\{x \mid A x=b, x \geq 0\}\) 是非空的, 则它至少有一个顶点。基本可行解 \(\Leftrightarrow R\) 的顶点。
推论2: 若线性规划存在有限的最优解, 则至少存在一个 \(R\) 的顶点是有限最优解。
推论3：线性规划的可行域 \(R\) 有有限多个顶点。基本可行解最多 \(C_n^m\) 个。
为了选择入基变量$x_j $，我们要求其检验数为正：
\[\sigma_j = c_j - C_B^T d > 0 \]这表明增加\(x_j\)可以使目标函数\(Z\)增加。因此，目标函数的增量为：
\[\Delta Z = \sigma_j \theta > 0 \quad \text{当} \quad \theta > 0 \]假设$d \leq 0 $，即方向向量中所有元素都小于等于零，即单纯形迭代时有一列元素都非正，这意味着：
\[X_B' = X_B - \theta d \geq 0 \]对于任意大的\(\theta\)都成立，因为$\theta d \leq 0 $，所以：
\[X_B \geq \theta d \geq 0 \]这意味着我们可以将\(x_j\)增加到任意大而不违反非负约束$X \geq 0 $。
- 由于\(\theta\)可以无限大，并且每增加一个单位的\(x_j\)都能使目标函数值增加，因此目标函数可以增大到无穷大。
- 这表明，如果方向向量\(d \leq 0\)且对应检验数$\sigma_j > 0 $，则线性规划问题的目标函数是无界的。这种情况仅仅是理论上的论述，在经济社会实践活动中永远不会发生。
3.2几何意义
这里以上面 Wyndor Glass 问题为例阐述线性规划求解几何意义。该例题的模型和图形见上图所示。其中，每个约束边界是一条直线，约束边界的交点是这个问题的顶点解（Corner-Point Solution），在可行域上的顶点解被称为顶点可行解（CPF solution）。在 Wyndor Glass 公司的例子中，两条约束边界产生一个 CPF 解，因此，每个 CPF 解都有两个相邻的 CPF 解（每个都在两条边之一的另一端，正如下表所列举的那样）。
| CPF 解 | 相邻 CPF 解 |
|---|---|
| (0, 0) | (0, 6) 和 (4, 0) |
| (0, 6) | (2, 6) 和 (0, 0) |
| (2, 6) | (4, 3) 和 (0, 6) |
| (4, 3) | (4, 0) 和 (2, 6) |
| (4, 0) | (0, 0) 和 (4, 3) |
我们之所以对相邻的 CPF 解感兴趣，是因为的这些解的基本特性能够为我们判断某个 CPF 解是否最优提供有效的思路。考虑任一个至少拥有一个最优解的线性规划问题。如果一个 CPF 解没有比它更好（以 Z 来衡量）的相邻 CPF 解，那么它就是最优解。例如根据上表，由于 CPF 解 (2,6) 的 Z=36，大于 (0,6) 对应的 Z=30 和 (4,3) 对应的 Z=27，因此，它必定是最优的。
import numpy as np
import matplotlib.pyplot as plt
# 定义 x1 和 x2 的范围
x1 = np.linspace(0, 5, 400)
x2 = np.linspace(0, 7, 400)
# 定义约束条件
x2_3 = (18 - 3 * x1) / 2 # 工厂 3 的约束
# 定义各个顶点
vertices = [(0, 0), (0, 6), (2, 6), (4, 3), (4, 0)]
# 计算顶点的 Z 值
z_values = [3000 * v[0] + 5000 * v[1] for v in vertices]
# 设置图形
plt.figure(figsize=(8, 6))
# 画出约束线
plt.axvline(x=4, color='blue', linestyle='--', label='$x_1 \\leq 4$') # 垂直线表示x1 <= 4
plt.axhline(y=6, color='green', linestyle='--', label='$x_2 \\leq 6$') # 水平线表示x2 <= 6
plt.plot(x1, x2_3, label='$3x_1 + 2x_2 \\leq 18$', color='orange')
# 填充可行域
y1 = np.minimum(6, np.maximum(0, (18 - 3 * x1) / 2))
plt.fill_between(x1, 0, y1, where=(x1 <= 4) & (y1 >= 0), color='grey', alpha=0.5)
# 标出顶点和 Z 值
for i, v in enumerate(vertices):
plt.scatter(*v, color='red')
plt.text(v[0], v[1], f'Z={z_values[i]}', fontsize=10, ha='center')
# 图形设置
plt.xlim(0, 5)
plt.ylim(0, 7)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.title('Feasible Region for the Linear Programming Problem')
plt.legend()
plt.grid()
plt.show()
四、示例
根据线性规划基本定理，我们可以通过穷举所有顶点来求解线性规划的最优解。下面以一个标准形式的线性规划问题为例，展示求解过程。
4.1 求解顶点步骤
Step 1: 令 \(x_3 = 0, x_4 = 0\)，代入(4.2-4.4)，构成方程组：
解得顶点 \(x^{(1)}\) 。
Step 2: 令\(x_2 = 0, x_3 = 0\)，代入(4.2-4.4)，构成方程组：
解得顶点\(x^{(2)}\) 。
Step 3: 令\(x_1 = 0, x_4 = 0\)，代入(4.2-4.4)，构成方程组：
解得顶点\(x^{(3)}\)。
Step 4: 令\(x_1 = 0, x_2 = 0\)，代入(4.2-4.4)，构成方程组：
解得顶点\(x^{(4)}\) 。
Step 5: 令\(x_1 = 0, x_3 = 0\)，代入(4.2-4.4)，构成方程组：
解得顶点\(x^{(5)}\) 。
Step 6: 令\(x_2 = 0, x_4 = 0\)，代入(4.2-4.4)，构成方程组：
解得顶点\(x^{(6)}\)。
4.2 顶点坐标
- $x^{(1)} = (0.5, 2.5, 0, 0) $
- $x^{(2)} = (3, 0, 0, 5) $
- \(x^{(3)} = (0, 2, 1, 0)\)
- $ x^{(4)} = (0, 0, 3, 2) $
- $ x^{(5)} = (0, 3, 0, -1)$ （舍弃，因为包含负数分量）
- $ x^{(6)} = (-2, 0, 5, 0)$ （舍弃，因为包含负数分量）
舍弃包含负数分量的顶点\(x^{(5)}\)和\(x^{(6)}\)，只保留顶点 \(x^{(1)}\)，\(x^{(2)}\)，\(x^{(3)}\)和\(x^{(4)}\)。
import numpy as np
import matplotlib.pyplot as plt
# 定义x1的范围
x1 = np.linspace(0, 4, 400)
# 根据约束条件定义x2
y1 = 3 - x1 # x1 + x2 <= 3
y2 = x1 + 2 # -x1 + x2 <= 2, 等价于 x2 <= x1 + 2
# 绘制图形
plt.figure(figsize=(8, 6))
plt.plot(x1, y1, label='$x_1 + x_2 = 3$', color='r')
plt.plot(x1, y2, label='$-x_1 + x_2 = 2$', color='b')
# 填充可行域
plt.fill_between(x1, 0, np.minimum(y1, y2), where=(np.minimum(y1, y2) >= 0), color='gray', alpha=0.5)
# 标记顶点和最优解
vertices = [(0, 2), (3, 0), (0.5, 2.5), (0, 0)] # 顶点坐标
optimal_vertex = (0.5, 2.5) # 最优解
for vertex in vertices:
plt.plot(vertex[0], vertex[1], 'ko') # 黑色圆点
plt.text(vertex[0], vertex[1], f' {vertex}', fontsize=12)
# 标记最优解
plt.plot(optimal_vertex[0], optimal_vertex[1], 'go', markersize=10, label='Optimal Solution')
plt.text(optimal_vertex[0], optimal_vertex[1], f' Optimal {optimal_vertex}', fontsize=12, color='green')
# 绘制目标函数方向
c = np.array([1, 2]) # 目标函数的系数
x1_opt = np.linspace(0, 4, 400)
x2_opt = (5.5 - c[0] * x1_opt) / c[1] # 使用目标值5.5绘制等值线
plt.plot(x1_opt, x2_opt, '--', label='Objective Function: $x_1 + 2x_2 = 5.5$', color='purple')
# 设置图形标签
plt.xlim(0, 4)
plt.ylim(0, 4)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.legend()
plt.title('Feasible Region and Optimal Solution of Linear Programming (Maximization)')
# 显示图形
plt.grid(True)
plt.show()
总结
线性规划的求解理论基于可行域的凸性和目标函数的线性特性，通过探索可行域的顶点或在内部逐渐接近最优点，最终找到使目标函数达到极值的解。无论是单纯形法还是内点法，都利用了线性结构的优势，确保能够高效地找到最优解。线性规划的求解过程主要包括下面一些步骤。
建立模型：根据实际问题的需求，构造线性目标函数和一组线性约束条件。目标函数通常表示我们希望优化的量，如利润或成本。约束条件表示资源的限制，如时间、材料等。
识别可行域： 可行域是所有满足线性约束条件的点的集合。由于约束条件是线性的，因此可行域是一个凸多边形或多面体。线性规划的一个重要性质是：如果有解，最优解必定位于可行域的边界上，且通常是在顶点处。
寻找最优解： 线性规划的核心在于寻找可行域上的一个点，使目标函数达到最优值。求解线性规划问题的经典方法是单纯形法（Simplex Method），它的基本思想是：从可行域的一个顶点开始，沿着边移动到相邻的顶点，不断改善目标函数的值，直到达到最优解或确定无界解。
验证最优性和可行性： 在找到一个候选解后，需要验证其是否是可行域中的点，以及是否能实现目标函数的最优值。通过检查约束条件是否被满足，以及目标函数在该点的值，来验证这一点。
敏感性分析： 在得到最优解之后，通常还会进行敏感性分析，以评估解对输入参数变化的敏感度。敏感性分析有助于理解模型的稳定性和可靠性，并指导决策者在实际操作中进行调整。
单纯形法之所以有效，是因为线性目标函数在可行域内的变化是单调的。因此，在移动过程中，一旦到达一个顶点，若在该点处目标函数不能进一步改进，则该点即为最优解。单纯形法通过系统性地探索顶点，能够保证找到最优解。
