单纯形法（Simplex Method）是解决线性规划问题的一种高效且广泛使用的算法。由乔治·丹齐克（George Dantzig）在20世纪40年代提出，这一方法通过系统地检查可行解空间的极点，从而找到最优解。由于其计算效率高，单纯形法迅速成为线性规划问题中最重要和最常用的算法之一。它的应用范围广泛，能够有效解决实际中的大规模优化问题，因此在现代工业和经济管理中扮演着关键角色。
| 初始单纯形表 | 迭代后单纯形表 |
|---|---|
 |
 |
一、单纯形表的消元法建构
线性规划的标准型为
其中： \(\boldsymbol{A}\) 是 \(m \times n\) 矩阵， \(m \leqslant n\) 且秩 \(r(\boldsymbol{A})=m\) ，即 \(\boldsymbol{A}\) 中至少有一个 \(m \times m\) 满秩子矩阵。
1.1 基本可行解的代数消元(最优解判断条件)
不妨设 \(\boldsymbol{B}=\left(\begin{array}{llll}\boldsymbol{P}_1 & \boldsymbol{P}_2 & \cdots & \boldsymbol{P}_m\end{array}\right)\) 是线性规划的一个基, 将有关矩阵和向量分块，记 \(\boldsymbol{A}=(\boldsymbol{B}, N), C=\left(\boldsymbol{C}_B, \boldsymbol{C}_N\right), X=\left(\boldsymbol{X}_B, \boldsymbol{X}_N\right)^{\mathrm{T}}\) ，为求线性规划的基本最优解，先要求线性规划的基本可行解，这样就需将约束中的基变量用非基变量表示出来。
用 \(\boldsymbol{B}^{-1}\) 左乘约束方程 \(\boldsymbol{A} \boldsymbol{X}=\boldsymbol{b}\) 的两端，得
即
其中 \(\boldsymbol{E}\) 是单位矩阵。整理, 得
将其代入目标函数中，得
即
非基变量 \(x_j\) 前面的系数 \(c_j-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j\) 称为变量 \(x_j\) 的检验数，表示该变量增加或减少一个单位所引起目标函数值的变化，它可以用来判断将 \(x_j\) 变成基变量后能否改进目标函数值，以后记为 \(\sigma_j=c_j-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j\) 。
最优解判定定理：对某基本可行解 \(\boldsymbol{X}_B=\boldsymbol{B}^{-1} \boldsymbol{b}\) 其他 \(\boldsymbol{x}_N=0\), 若所有的 \(\sigma_j=\boldsymbol{c}_{\boldsymbol{j}}\) \(-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j \leqslant 0\) ，则该解为最优解。
1.2 基本可行解的矩阵消元
考虑线性方程组 \(\left\{\begin{array}{l}\boldsymbol{A} \boldsymbol{X}=\boldsymbol{b} \\ \boldsymbol{z}=\boldsymbol{C} \boldsymbol{X}\end{array}\right.\), 其变量为 \(\left[\begin{array}{l}\boldsymbol{X} \\ \boldsymbol{z}\end{array}\right]\), 为便于求解, 整理得方程组 \(\left\{\begin{array}{c}0 \cdot z+\boldsymbol{A} \boldsymbol{X}=\boldsymbol{b} \\ -z+\boldsymbol{C X}=\boldsymbol{0}\end{array}\right.\), 其增广矩阵见下表, 应用高斯消元法, 求解线性方程组的解。
| 常数项 | z | XB | XN | |
|---|---|---|---|---|
| b | 0 | B | N | （1-1） |
| 0 | -1 | \(C_B\) | \(C_N\) | （1-2） |
用 \(\boldsymbol{B}^{-1}\) 左乘方程组(1-1)的两端, 将 \(\boldsymbol{X}_B\) 的系数化为单位矩阵, 得下表。
| 常数项 | z | $$X_B$$ | $$X_N$$ | |
|---|---|---|---|---|
| \(B^{-1} b\) | 0 | \(B^{-1}B=E\) | \(B^{-1}N\) | （1-3） |
| 0 | -1 | \(C_B\) | \(C_N\) |
将方程组（1-3）左乘 \(-\boldsymbol{C}_B\) 加到方程（1-2）两边，消元化简得下表。
| 常数项 | z | \(X_B\) | \(X_N\) | |
|---|---|---|---|---|
| \(B^{-1} b\) | 0 | $$B^{-1}B=E$$ | $$B^{-1}N$$ | |
| $$-C_BB^{-1}b$$ | -1 | $$C_B - C_BB^{-1}B = 0$$ | $$C_N - C_BB^{-1}N $$ | （1-4） |
注意式(1-4)中基向量 \(\boldsymbol{X}_B\) 、非基向量 \(\boldsymbol{X}_N\) 的系数，它们形式相似，当前者 \(\boldsymbol{C}_B-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{B}=0\)时，后者即为检验数 \(\boldsymbol{C}_{\boldsymbol{N}}-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{N}\) 的矩阵形式。也就是说，用消元法将目标函数中基变量的系数化为零的同时，就会得到非基变量的检验数。事实上， \(\boldsymbol{C}_B-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{B}\) 也可看成基向量的检验数, 应用消元法后, 当基向量的检验数化为零时, 非基变量的系数 \(\boldsymbol{C}_{\boldsymbol{N}}-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{N}\) 就是其检验数。
从上述论述可知, 计算检验数的方法有两种：一是通过公式 \(\sigma_j=c_j-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j\) 计算得出；二是应用消元法，将目标函数中基向量的系数都化为零时，非基向量的系数就是检验数。这也是单纯表结构的出处。
1.3 单纯形法=逆矩阵法
单纯形法的核心思想是通过变换基本变量和非基本变量，找到一个新的基本可行解，并通过比较检验数判断该解是否为最优解。在单纯形法中，每一个基本可行解都对应着一个基矩阵的逆矩阵。基矩阵是从约束条件系数矩阵中选取的列所构成的矩阵，它的逆矩阵在每次迭代中都要更新。
这个逆矩阵的求解是单纯形法消元过程的核心，因为它直接影响到新解的生成和目标函数的优化。具体来说，当我们选择一个进入基的变量并排除一个离开基的变量时，这相当于对基矩阵进行了一次列替换操作。为了保持计算的有效性，我们需要快速更新逆矩阵。这通常通过一个修正的高斯消元法来完成，使得新基的逆矩阵可以通过旧基的逆矩阵迅速计算出来。
单纯形的每次迭代，都可以使用逆矩阵来计算检验数，从而判断是否可以进一步优化目标函数。检验数的计算过程依赖于逆矩阵，因为检验数用于评估当前解的可行性及其对目标函数的影响。如果所有检验数都满足最优性条件，那么当前解就是最优解；否则，我们需要通过逆矩阵的进一步操作继续迭代。
例1：已知初始单纯形表和最终单纯形表, 试求解以下问题。
（1）在初始单纯形表中找出最优基 \(\boldsymbol{B}\) ，在最终单纯形表里找出 \(\boldsymbol{B}^{-1}\) 。
（2）完成最终单纯形表。(3) 给出最优解与最优值。
| $$C_B$$ | $$X_B$$ | $$B^{-1}b$$ | $$x_1$$ | $$x_2$$ | $$x_3$$ | $$x_4$$ | $$x_5$$ | $$x_6$$ | |
|---|---|---|---|---|---|---|---|---|---|
| 初始表 | 0 | \(x_4\) | 60 | 3 | 1 | 1 | 1 | 0 | 0 |
| 0 | \(x_5\) | 10 | 1 | -1 | 2 | 0 | 1 | 0 | |
| 0 | \(x_6\) | 20 | 1 | 1 | -1 | 0 | 0 | 1 | |
| \(σ_j\) | 2 | -1 | 1 | 0 | 0 | 0 | |||
| 最终表 | \(x_4\) | -1 | -2 | ||||||
| \(x_1\) | 1/2 | 1/2 | |||||||
| \(x_2\) | -1/2 | 1/2 | |||||||
| \(σ_j\) |
- 问题描述
已知初始单纯形表和最终单纯形表，求解以下问题：
(1)在初始单纯形表中找出最优基\(\boldsymbol{B}\)，在最终单纯形表里找出\(\boldsymbol{B}^{-1}\)。
(2)完成最终单纯形表。
(3)给出最优解与最优值。
初始单纯形表和最终单纯形表如下：
| $$C_B$$ | $$X_B$$ | $$B^{-1}b$$ | $$x_1$$ | $$x_2$$ | $$x_3$$ | $$x_4$$ | $$x_5$$ | $$x_6$$ | |
|---|---|---|---|---|---|---|---|---|---|
| 初始表 | 0 | \(x_4\) | 60 | 3 | 1 | 1 | 1 | 0 | 0 |
| 0 | \(x_5\) | 10 | 1 | -1 | 2 | 0 | 1 | 0 | |
| 0 | \(x_6\) | 20 | 1 | 1 | -1 | 0 | 0 | 1 | |
| \(σ_j\) | 2 | -1 | 1 | 0 | 0 | 0 | |||
| 最终表 | \(x_4\) | -1 | -2 | ||||||
| \(x_1\) | 1/2 | 1/2 | |||||||
| \(x_2\) | -1/2 | 1/2 | |||||||
| \(σ_j\) |
确定最优基和\(B^{-1}\)
- 最优基 \(\boldsymbol{B}\)：从初始基变量\(x_4\)、\(x_5\)、\(x_6\) 通过单纯形法操作，最终基变量变为\(x_4\)、\(x_1\)、\(x_2\)
- \(B^{-1}\)：通过最终单纯形表的\(x_4\)、\(x_5\)、\(x_6\)列得到 \(B^{-1}\)：
计算\(B^{-1}\)
- 右端项\(b = \begin{bmatrix} 60 \quad 10 \quad 20 \end{bmatrix}\)
计算最终单纯形表中的列向量
-\(x_3\)列向量计算：
- 初始表中\(x_3\)列为\([1, 2, -1]^T\)
因此\(x_3\)列为\(\begin{bmatrix} 1 \quad 0.5 \quad -1.5 \end{bmatrix}\)。
计算检验数\(σ_j\)
目标函数为\(z = 2x_1 -x_2+x_3\)
- 基变量的成本系数向量\(C_B = [0, 2, -1]\)
| 变量 | $$P_j$$ | $$σ_j = c_j-C_B B^{-1}P_j - $$ | 计算过程 |
|---|---|---|---|
| $$x_1 $$ | $$ \begin{bmatrix} 0 \quad 1 \quad 0 \end{bmatrix} $$ | 0 | $$ 2-[0, 2, -1] \cdot \begin{bmatrix} 0 \ 1 \ 0 \end{bmatrix} = 0 $$ |
| $$x_2$$ | $$ \begin{bmatrix} 0 \quad 0 \quad 1 \end{bmatrix} $$ | 0 | $$-1-[0, 2, -1] \cdot \begin{bmatrix} 0 \ 0 \ 1 \end{bmatrix} = 0 $$ |
| $$ x_3 $$ | $$ \begin{bmatrix} 1 \quad 0.5 \quad -1.5 \end{bmatrix} $$ | -1.5 | $$1-[0, 2, -1] \cdot \begin{bmatrix} 1 \ 0.5 \ -1.5 \end{bmatrix} = -1.5 $$ |
| $$x_4 $$ | $$ \begin{bmatrix} 1 \quad 0 \quad 0 \end{bmatrix} $$ | 0 | $$ 0-[0, 2, -1] \cdot \begin{bmatrix} 1 \ 0 \ 0 \end{bmatrix} = 0 $$ |
| $$ x_5$$ | $$ \begin{bmatrix} -1 \quad 0.5 \quad -0.5 \end{bmatrix} $$ | -1.5 | $$0-[0, 2, -1] \cdot \begin{bmatrix} -1 \ 0.5 \ -0.5 \end{bmatrix} = -1.5 $$ |
| $$ x_6 $$ | $$ \begin{bmatrix} -2 \quad 0.5 \quad 0.5 \end{bmatrix} $$ | -0.5 | $$0- [0, 2, -1] \cdot \begin{bmatrix} -2 \ 0.5 \ 0.5 \end{bmatrix} = -0.5 $$ |
最终单纯形表
| $$C_B$$ | $$X_B$$ | $$B^{-1}b$$ | $$x_1$$ | $$x_2$$ | $$x_3$$ | $$x_4$$ | $$x_5$$ | $$x_6$$ | |
|---|---|---|---|---|---|---|---|---|---|
| 最终表 | 0 | $$x_4$$ | 10 | 0 | 0 | 1 | 1 | -1 | -2 |
| 3 | $$x_1$$ | 15 | 1 | 0 | 0.5 | 0 | 1/2 | 1/2 | |
| 2 | $$x_2$$ | 5 | 0 | 1 | -1.5 | 0 | -1/2 | 1/2 | |
| $$σ_j$$ | 0 | 0 | -1.5 | 0 | -1.5 | -0.5 |
最优解与最优值
- 最优解：^*=\left[\begin{array}{llllll}15 & 5 & 0 & 10 & 0 & 0\end{array}\right]^T。
- 最优值：\(X, z^*=25\)
二、单纯形的迭代步骤
例1：对于线性规划问题
其中， \(x_3\) 和 \(x_4\) 是松驰变量。
决策变量：\(\mathbf{x}=\left[\begin{array}{l}x_1& x_2 & x_3 & x_4\end{array}\right]^T\) 是决策变量向量。
目标函数: Max \(\mathbf{c}^{\top} \mathbf{x}\) 其中， \(\mathbf{c}=\left[\begin{array}{c}2\quad 3\quad 0 \quad0\end{array}\right]^T\) 为目标函数的系数向量。
约束条件: \(\mathbf{A} \mathbf{x}=\mathbf{b}\) ，其中 \(\mathbf{A}=\left[\begin{array}{cccc}4 & 1 & 1 & 0 \\ -1 & 1 & 0 & 1\end{array}\right]\) 是约束系数矩阵。
\(\mathbf{b}=\left[\begin{array}{c}16 & 6\end{array}\right]^T\) 是约束右端项向量。因此，标准型的线性规划问题可以表示为：
初始时令 \(x_N=[0,0]\)
**第一轮迭代: **此时 \(c_N-c_B B^{-1} N=[2,3]\) 因此选择 \(x_2\) 作为入基变量更为高效，且 \(b-N x_N=\left[\begin{array}{c}16-4 x_1-x_2 \\ 6+x_1-x_2\end{array}\right]>=0\) 故 \(x_2=6\) 并令 \(x_4=0, x_2\) 入基, \(x_4\) 出基。经过此轮迭代后，各变量如下
**第二轮迭代: **此时 \(c_N-c_B B^{-1} N=[5,-3]\) ，选择 \(x_1\) 作为入基变量更为高效，且
\(b-N x_N=\left[\begin{array}{c}16-4 x_1 \\ 6+x_1\end{array}\right]>=0\) 故 \(x_1=4\) 并令 \(x_3=0, x_1\) 入基， \(x_3\) 出基。经过此轮迭代后，各变量如下
**第三轮迭代: **此时 \(c_N-c_B B^{-1} N=[-1,-2]\) 都小于 0 ，达到收敛条件，此时令 \(x_N=[0,0]\) 解得 \(x_B=[2,8]\) 最小值为28。
**例2：*求解下面线性模型
初始单纯形表
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(x_5\) | \(b\) | \(\theta\) | |
|---|---|---|---|---|---|---|---|
| 目标函数 | 70 | 30 | 0 | 0 | 0 | ||
| 约束 1 | 3 | 9 | 1 | 0 | 0 | 540 | |
| 约束 2 | 5 | 5 | 0 | 1 | 0 | 450 | |
| 约束 3 | 9 | 3 | 0 | 0 | 1 | 720 |
判断当前顶点是否是最优解
对于最大化问题,若当前目标函数中的非基变量的系数小于等于0时，则所得的解为最优解。而在当前的例子中，非基变量的系数分别为70和30，意味着在可行域内随着非基变量\(x_1\)和\(x_2\)的增大，目标函数就会继续增大，因此当前的解不是最优解。故判断当前所得的解是否为最优解时，只需判断目标函数中非基变量的系数是否小于等于0。
进基和出基变量
变量的出基与入基，在几何图像上表现为顶点的变化。入基的规则为选择使目标函数z变化最快的非基变量入基，即选择系数最大且为正数的非基变量入基，故在本例中选择\(x_1\)入基。出基的规则则需要引入一个新的量\(\theta\)，\(\theta=b/a_i\)（\(a_i\)为非基变量系数，\(a_i\)的选择的是刚刚入基的非基变量的系数），选择最小的\(\theta\)出基。
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(x_5\) | \(b\) | \(\theta\) | |
|---|---|---|---|---|---|---|---|
| 目标函数 | 70 | 30 | 0 | 0 | 0 | ||
| 约束 1 | 3 | 9 | 1 | 0 | 0 | 540 | 180 |
| 约束 2 | 5 | 5 | 0 | 1 | 0 | 450 | 90 |
| 约束 3 | 9 | 3 | 0 | 0 | 1 | 720 | 80 |
第一次迭代，经过\(x_1\)入基与\(x_5\)出基的运算后，得到的结果如下所示。
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(x_5\) | \(b\) | \(\theta\) | |
|---|---|---|---|---|---|---|---|
| \(Z\) | 0 | 20/3 | 0 | 0 | -70/9 | 5600 | |
| \(x_3\) | 0 | 8 | 1 | 0 | -1/3 | 300 | |
| \(x_4\) | 0 | 10/3 | 0 | 1 | -5/9 | 50 | |
| \(x_1\) | 1 | 1/3 | 0 | 0 | 1/9 | 80 |
令非基变量的值为0，则得到第一次迭代的解，如下所示
判断该解是否为最优解(重复第二个步骤，直到是最优解为止)，显然由于非基变量\(x_2\)的系数大于0，故当前位置还不是最优解，再次重复上面步骤。
三、练习
| Food | Energy（能量） | Protein（蛋白质） | Calcium（钙） | Price |
|---|---|---|---|---|
| Oatmeal（燕麦） | 110 | 4 | 2 | 3 |
| Whole milk（全奶） | 160 | 8 | 285 | 9 |
| Cherry pie（草莓派） | 420 | 4 | 22 | 20 |
| Pork with beans（猪肉） | 260 | 14 | 80 | 19 |
决策变量：
- \(x_1\)：购买燕麦（Oatmeal）的数量
- \(x_2\)：购买全奶（Whole milk）的数量
- \(x_3\)：购买草莓派（Cherry pie）的数量
- \(x_4\)：购买猪肉（Pork with beans）的数量
目标函数：
最小化总价格：
约束条件：
- 能量需求：至少2000单位能量
- 蛋白质需求：至少55单位蛋白质
- 钙需求：至少800单位钙
- 非负性约束
from pulp import LpProblem, LpMinimize, LpVariable, lpSum, value, LpStatus# 定义问题
problem = LpProblem("Minimize Cost", LpMinimize)# 定义变量
x1 = LpVariable('x1', lowBound=0, cat='Continuous')
x2 = LpVariable('x2', lowBound=0, cat='Continuous')
x3 = LpVariable('x3', lowBound=0, cat='Continuous')
x4 = LpVariable('x4', lowBound=0, cat='Continuous')# 目标函数
problem += 3*x1 + 9*x2 + 20*x3 + 19*x4, "Total Cost"# 约束条件
problem += 110*x1 + 160*x2 + 420*x3 + 260*x4 >= 2000, "Energy"
problem += 4*x1 + 8*x2 + 4*x3 + 14*x4 >= 55, "Protein"
problem += 2*x1 + 285*x2 + 22*x3 + 80*x4 >= 800, "Calcium"# 求解
problem.solve()# 输出结果
print(f"Status: {LpStatus[problem.status]}")
print(f"Oatmeal (燕麦) units: {value(x1)}")
print(f"Whole milk (全奶) units: {value(x2)}")
print(f"Cherry pie (草莓派) units: {value(x3)}")
print(f"Pork with beans (猪肉) units: {value(x4)}")
print(f"Total minimum cost: {value(problem.objective)}")
Status: Optimal
Oatmeal (燕麦) units: 14.24428
Whole milk (全奶) units: 2.7070577
Cherry pie (草莓派) units: 0.0
Pork with beans (猪肉) units: 0.0
Total minimum cost: 67.09635929999999
总结
单纯形法是线性规划问题中的一种经典算法，它通过逐步优化，找到能够最大化或最小化目标函数的可行解。尽管单纯形法在最坏情况下可能需要指数级时间，但在实际应用中，单纯形法通常能高效地找到最优解，因此它被广泛应用于各种线性规划问题中，例如生产计划、资源分配、运输优化等。
单纯形法的直观性在于从一个顶点开始沿边界移动到另一个顶点，不断提高目标函数值，直到达到最优解。这种方法简单且易于理解，对大多数线性规划问题都能有效求解。然而，单纯形法也存在一些局限性，例如在处理退化问题时可能会出现循环，导致算法陷入无限循环的困境。为了克服这些局限性，研究者们提出了多种改进方案，如反周期规则来防止循环、对偶单纯形法来处理不可行的初始解等。此外，内点法作为一种替代算法，通过从可行域的内部逼近最优解，提供了不同于单纯形法的求解思路，并且在某些情况下表现出更好的最坏情况性能。通过这些改进和新方法的引入，线性规划问题的求解在理论和实践上都得到了极大的丰富和发展，使得我们能够解决更复杂、更大规模的问题。这些进步不仅提升了算法的效率，也拓展了线性规划在各个领域的应用范围。单纯形法及其改进方法的持续发展，确保了线性规划在优化和决策问题中的核心地位。
参考文献
- https://zhuanlan.zhihu.com/p/672071565
- python最优化算法实战---线性规划之单纯形法
在给定的营业时间内（14小时），啤酒厂需要制定一个生产计划，以最大化其收益。生产生啤需要1小时，生产黑啤需要2小时。生啤的售价为20美元/箱，黑啤的售价为30美元/箱。每日销售上限为100箱。
