线性规划单纯形法精解
单纯形法（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\) |
确定最优基和\(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}^T\)。
计算检验数\(σ_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 |
最优解与最优值
- 最优解：\(X^*=\left[\begin{array}{llllll}15 & 5 & 0 & 10 & 0 & 0\end{array}\right]^T\)。
- 最优值：\(z^*=25\)
二、单纯形的迭代步骤
例2：对于线性规划问题
其中， \(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}^T \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\) 作为入基变量更为高效，且 $$\theta =\min [\frac{16}{1},\frac{6}{1}]=6 $$ 故\(x_2\) 入基, \(x_4\) 出基。经过此轮迭代后，各变量如下
第二轮迭代:此时 \(c_N-c_B B^{-1} N=[5,-3]\) ，选择 \(x_1\) 作为入基变量更为高效，且
故\(x_2\) 故 \(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。
| 可行域顶点 | 迭代路径 |
|---|---|
 |
 |
例3：求解下面线性模型
初始单纯形表
| $$x_1$$ | $$x_2$$ | $$x_3$$ | $$x_4$$ | $$x_5$$ | $$b$$ | $$\theta$$ | |
|---|---|---|---|---|---|---|---|
| 约束 1 | 3 | 9 | 1 | 0 | 0 | 540 | |
| 约束 2 | 5 | 5 | 0 | 1 | 0 | 450 | |
| 约束 3 | 9 | 3 | 0 | 0 | 1 | 720 | |
| 目标函数 | 70 | 30 | 0 | 0 | 0 |
判断当前顶点是否是最优解
对于最大化问题,若当前目标函数中的非基变量的系数小于等于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$$ | |
|---|---|---|---|---|---|---|---|
| 约束 1 | 3 | 9 | 1 | 0 | 0 | 540 | 180 |
| 约束 2 | 5 | 5 | 0 | 1 | 0 | 450 | 90 |
| 约束 3 | 9 | 3 | 0 | 0 | 1 | 720 | 80 |
| 目标函数 | 70 | 30 | 0 | 0 | 0 |
第一次迭代，经过\(x_1\)入基与\(x_5\)出基的运算后，得到的结果如下所示。
| $$x_1$$ | $$x_2$$ | $$x_3$$ | $$x_4$$ | $$x_5$$ | $$b$$ | $$\theta$$ | |
|---|---|---|---|---|---|---|---|
| \(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 | |
| \(Z\) | 0 | 20/3 | 0 | 0 | -70/9 | 5600 |
令非基变量的值为0，则得到第一次迭代的解，如下所示
判断该解是否为最优解(重复第二个步骤，直到是最优解为止)，显然由于非基变量\(x_2\)的检验数大于0，故当前位置还不是最优解，再次重复上面步骤。
三、练习
例4：用单纯表求解下面线性规划
import pandas as pd
import numpy as np
"""
系数矩阵的形式:
b x1 x2 x3 x4 x5
obj 0 70 30 0 0 0
x3 540 3 9 1 0 0
x4 450 5 5 0 1 0
x5 720 9 3 0 0 1
①第一行是目标函数的系数；亦即各变量对应的检验数；第2~4行是约束条件的系数
②第一列是约束方程的常数项
③对于目标函数的更新我们同样采用矩阵的变换,所以obj对应的第一列表示的是目标函数的相反数
"""
"""
运行如下代码后得到的结果如下所示
最终的最优单纯性法是:
b x1 x2 x3 x4 x5
obj -5700 0 0.0 0 -2.0 -6.666667
x3 180 0 0.0 1 -2.4 1.000000
x2 15 0 1.0 0 0.3 -0.166667
x1 75 1 0.0 0 -0.1 0.166667
目标函数的值: 5700
最优决策变量是:
x1 = 75
x2 = 15
"""
# Initial Simplex tableau, set dtype to float to handle division correctly
matrix = pd.DataFrame(
data=np.array([
[0, 2, -1, 1, 0, 0, 0],
[60, 3, 1, 1, 1, 0, 0],
[10, 1, -1, 2, 0, 1, 0],
[20, 1, 1, -1, 0, 0, 1]
], dtype=float), # Explicitly set dtype to float
index=['obj', 'x4', 'x5', 'x6'],
columns=['b', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6']
)
# Print the initial tableau
iteration = 1
print(f"第 {iteration} 个单纯性表是:")
print(matrix)
# Continue the iteration until no positive coefficients are in the objective function row
while True:
# Check the coefficients of the objective function row (检验数)
c = matrix.iloc[0, 1:]
# Stop if all reduced costs are less than or equal to zero (optimal solution found)
if c.max() <= 0:
break
iteration += 1
# Choose entering variable (the one with the highest positive coefficient)
in_x = c.idxmax()
in_x_v = c[in_x] # Coefficient of the entering variable
# Choose leaving variable by calculating θ and finding the minimum
b = matrix.iloc[1:, 0]
in_x_a = matrix.iloc[1:][in_x]
# Ensure no division by zero, and select the minimum positive theta value
theta = b / in_x_a
theta = theta[theta > 0] # Consider only positive θ values
if theta.empty: # If theta is empty, break to avoid errors (degenerate case)
print("No feasible solution found.")
break
out_x = theta.idxmin()
# Pivoting: adjust the rows for the entering and leaving variables
matrix.loc[out_x, :] = matrix.loc[out_x, :] / matrix.loc[out_x, in_x]
for idx in matrix.index:
if idx != out_x:
matrix.loc[idx, :] -= matrix.loc[out_x, :] * matrix.loc[idx, in_x]
# Update the basis
basis = matrix.index.tolist()
basis[basis.index(out_x)] = in_x
matrix.index = basis
# Print current tableau
print(f"第 {iteration} 个单纯性表是:")
print(matrix)
# Output the final results
print("最终的最优单纯形表是:")
print(matrix)
print("目标函数的值:", -matrix.iloc[0, 0])
# Determine the optimal decision variables
x_count = (matrix.shape[1] - 1) - (matrix.shape[0] - 1)
X = matrix.iloc[0, 1:].index.tolist()[:x_count]
print("最优决策变量是:")
for xi in X:
if xi in matrix.index: # Check if the variable is in the current basis
print(f"{xi} = {matrix.loc[xi, 'b']}")
else:
print(f"{xi} = 0") # If not in basis, the value is zero
第 1 个单纯性表是:
b x1 x2 x3 x4 x5 x6
obj 0.0 2.0 -1.0 1.0 0.0 0.0 0.0
x4 60.0 3.0 1.0 1.0 1.0 0.0 0.0
x5 10.0 1.0 -1.0 2.0 0.0 1.0 0.0
x6 20.0 1.0 1.0 -1.0 0.0 0.0 1.0
第 2 个单纯性表是:
b x1 x2 x3 x4 x5 x6
obj -20.0 0.0 1.0 -3.0 0.0 -2.0 0.0
x4 30.0 0.0 4.0 -5.0 1.0 -3.0 0.0
x1 10.0 1.0 -1.0 2.0 0.0 1.0 0.0
x6 10.0 0.0 2.0 -3.0 0.0 -1.0 1.0
最终的最优单纯形表是:
b x1 x2 x3 x4 x5 x6
obj -25.0 0.0 0.0 -1.5 0.0 -1.5 -0.5
x4 10.0 0.0 0.0 1.0 1.0 -1.0 -2.0
x1 15.0 1.0 0.0 0.5 0.0 0.5 0.5
x2 5.0 0.0 1.0 -1.5 0.0 -0.5 0.5
目标函数的值: 25.0
最优决策变量是:
x1 = 15.0
x2 = 5.0
x3 = 0
总结
单纯形法是线性规划问题中的一种经典算法，它通过逐步优化，找到能够最大化或最小化目标函数的可行解。尽管单纯形法在最坏情况下可能需要指数级时间，但在实际应用中，单纯形法通常能高效地找到最优解，因此它被广泛应用于各种线性规划问题中，例如生产计划、资源分配、运输优化等。
单纯形法的直观性在于从一个顶点开始沿边界移动到另一个顶点，不断提高目标函数值，直到达到最优解。这种方法简单且易于理解，对大多数线性规划问题都能有效求解。然而，单纯形法也存在一些局限性，例如在处理退化问题时可能会出现循环，导致算法陷入无限循环的困境。为了克服这些局限性，研究者们提出了多种改进方案，如反周期规则来防止循环、对偶单纯形法来处理不可行的初始解等。此外，内点法作为一种替代算法，通过从可行域的内部逼近最优解，提供了不同于单纯形法的求解思路，并且在某些情况下表现出更好的最坏情况性能。通过这些改进和新方法的引入，线性规划问题的求解在理论和实践上都得到了极大的丰富和发展，使得我们能够解决更复杂、更大规模的问题。这些进步不仅提升了算法的效率，也拓展了线性规划在各个领域的应用范围。单纯形法及其改进方法的持续发展，确保了线性规划在优化和决策问题中的核心地位。
