用Python手搓单纯形法求解器：从理论到代码的深度实践

为什么我们需要自己实现单纯形法？

在运筹学和优化领域，单纯形法就像是一把瑞士军刀——它可能不是最快的工具，但却是最可靠、最通用的线性规划求解方法之一。许多商业软件如Excel Solver、MATLAB都内置了这个算法，但作为学习者和实践者，亲自动手实现它能带来三大独特价值：

1. 透视算法本质 ：通过代码实现，你能真正理解检验数计算、基变换等关键步骤的数学含义
2. 调试思维训练 ：处理退化、循环等边界情况能大幅提升你的算法调试能力
3. 定制化扩展 ：自己实现的求解器可以轻松集成到特定业务场景中，不受商业软件限制

下面我们从一个实际的生产优化问题出发，逐步构建完整的单纯形法求解器。

1. 问题建模与标准化

案例：家具厂生产优化

假设一家家具厂生产桌子和椅子，每张桌子利润为3元，每把椅子利润为4元。生产受以下限制：

• 木材供应：2单位/桌 + 1单位/椅 ≤ 40
• 人工工时：1单位/桌 + 3单位/椅 ≤ 30

数学建模 ：

max Z = 3x₁ + 4x₂
s.t.:
    2x₁ + x₂ ≤ 40
    x₁ + 3x₂ ≤ 30
    x₁, x₂ ≥ 0

标准化转换

单纯形法要求约束条件为等式形式，我们引入松弛变量：

# 原始约束
2x₁ + x₂ + s₁ = 40  
x₁ + 3x₂ + s₂ = 30

# 目标函数调整
max Z = 3x₁ + 4x₂ + 0s₁ + 0s₂

Python实现标准化 ：

def standard_form(c, A, b):
    """将线性规划转化为标准形
    参数：
        c: 目标函数系数 [n,]
        A: 约束矩阵 [m,n]
        b: 约束右侧常数 [m,]
    返回：
        标准形系数 (c', A', b')
    """
    m, n = A.shape
    # 添加松弛变量
    c_slack = np.zeros(m)
    A_slack = np.eye(m)
    return (
        np.hstack([c, c_slack]),  # 新目标系数
        np.hstack([A, A_slack]),  # 新约束矩阵
        b.copy()                  # 约束右侧不变
    )

2. 单纯形表数据结构设计

高效的单纯形表实现是算法核心。我们使用NumPy数组存储，并设计专门的类来管理：

class SimplexTable:
    def __init__(self, c, A, b):
        self.c = c          # 目标函数系数
        self.A = A          # 约束矩阵
        self.b = b          # 右侧常数
        self.basis = [...]  # 当前基变量索引
        self.z = 0          # 当前目标值
        
    def compute_reduced_cost(self):
        """计算所有非基变量的检验数"""
        pass
    
    def pivot(self, entering, leaving):
        """执行旋转操作"""
        pass
    
    def get_solution(self):
        """提取当前解"""
        pass

关键数据结构示例 ：

CB	基变量	b	x₁	x₂	s₁	s₂	θ
0	s₁	40	2	1	1	0	40
0	s₂	30	1	3	0	1	10
-	σ	-	3	4	0	0	-

3. 算法核心步骤实现

3.1 初始可行解确定

def initialize(self):
    # 检查是否需要两阶段法
    if np.all(self.b >= 0):
        # 简单情况：松弛变量直接形成可行基
        self.basis = list(range(len(self.c)-self.A.shape[0], len(self.c)))
    else:
        # 需要两阶段法处理
        self.phase1()

3.2 检验数计算与入基选择

def select_entering(self):
    """选择入基变量——最大检验数规则"""
    reduced_costs = self.compute_reduced_cost()
    entering = np.argmax(reduced_costs)
    return entering if reduced_costs[entering] > 0 else None

3.3 θ比值测试与出基选择

def select_leaving(self, entering):
    """通过最小比值测试选择出基变量"""
    ratios = []
    for i in range(len(self.b)):
        if self.A[i, entering] > 0:
            ratios.append(self.b[i] / self.A[i, entering])
        else:
            ratios.append(np.inf)
    
    leaving = np.argmin(ratios)
    return leaving if ratios[leaving] < np.inf else None

3.4 基变换（旋转操作）

def pivot(self, entering, leaving):
    """执行旋转操作"""
    pivot_val = self.A[leaving, entering]
    
    # 更新主元行
    self.A[leaving] /= pivot_val
    self.b[leaving] /= pivot_val
    
    # 更新其他行
    for i in range(len(self.b)):
        if i != leaving and self.A[i, entering] != 0:
            factor = self.A[i, entering]
            self.A[i] -= factor * self.A[leaving]
            self.b[i] -= factor * self.b[leaving]
    
    # 更新基变量记录
    self.basis[leaving] = entering

4. 完整算法流程

将上述步骤组合成完整算法：

def solve(self):
    self.initialize()
    
    while True:
        entering = self.select_entering()
        if entering is None:  # 最优解条件
            break
            
        leaving = self.select_leaving(entering)
        if leaving is None:   # 无界情况
            raise ValueError("Problem is unbounded")
            
        self.pivot(entering, leaving)
    
    return self.get_solution()

5. 可视化迭代过程

理解单纯形法的几何意义至关重要。我们可以用Matplotlib展示解在可行域顶点间的移动：

import matplotlib.pyplot as plt

def plot_iteration(path, A, b):
    """绘制单纯形法的迭代路径"""
    fig, ax = plt.subplots()
    
    # 绘制约束条件
    x = np.linspace(0, 25, 100)
    for i in range(A.shape[0]):
        y = (b[i] - A[i,0]*x) / A[i,1]
        ax.plot(x, y, label=f'约束{i+1}')
    
    # 绘制可行域顶点
    vertices = compute_vertices(A, b)
    ax.scatter(*zip(*vertices), c='red')
    
    # 绘制迭代路径
    path = np.array(path)
    ax.plot(path[:,0], path[:,1], 'bo-')
    
    ax.set_xlim(0, 25)
    ax.set_ylim(0, 15)
    ax.legend()
    plt.show()

6. 处理特殊情况的技巧

6.1 退化与循环

当基变量取值为0时可能出现退化，导致算法循环。我们采用Bland规则：

def select_entering_bland(self):
    """Bland规则：选择下标最小的正检验数变量"""
    reduced_costs = self.compute_reduced_cost()
    for j in range(len(reduced_costs)):
        if reduced_costs[j] > 0:
            return j
    return None

6.2 初始可行解获取

当原点不是可行解时，需要使用两阶段法：

def phase1(self):
    """第一阶段：构造辅助问题寻找初始可行解"""
    # 添加人工变量构建辅助问题
    aux_c = np.zeros(self.A.shape[1])  # 原变量系数为0
    aux_c = np.hstack([aux_c, np.ones(self.A.shape[0])])  # 人工变量系数为1
    
    aux_A = np.hstack([self.A, np.eye(self.A.shape[0])])
    aux_table = SimplexTable(aux_c, aux_A, self.b.copy())
    
    # 解辅助问题
    aux_table.solve()
    
    if not np.isclose(aux_table.z, 0):
        raise ValueError("Problem is infeasible")
    
    # 转换到原问题
    self.basis = ...

7. 性能优化技巧

工业级实现需要考虑数值稳定性与计算效率：

def revised_simplex(self):
    """修正单纯形法：只存储基逆矩阵"""
    Binv = np.eye(len(self.basis))  # 基逆矩阵
    
    while True:
        # 计算对偶变量
        y = self.c[self.basis] @ Binv
        
        # 选择入基变量
        entering = self.select_entering_revised(y)
        
        # 计算入基方向
        d = Binv @ self.A[:, entering]
        
        # 选择出基变量
        leaving = self.select_leaving_revised(d)
        
        # 更新基逆矩阵
        Binv = self.update_basis_inverse(Binv, entering, leaving)

8. 完整代码架构

最终我们的求解器包含以下模块：

simplex_solver/
│── core/
│   ├── tableau.py      # 单纯形表实现
│   ├── phase1.py       # 第一阶段处理
│   └── revised.py      # 修正单纯形法
│── utils/
│   ├── io.py           # 模型读取/输出
│   └── visualization.py# 求解过程可视化
│── examples/           # 示例问题
└── tests/              # 单元测试

使用示例：

from simplex_solver import SimplexSolver

# 构建模型
c = np.array([3, 4])          # 目标函数
A = np.array([[2, 1], [1, 3]]) # 约束矩阵
b = np.array([40, 30])        # 右侧常数

# 求解
solver = SimplexSolver(c, A, b)
solution = solver.solve()

print(f"最优解: x1={solution.x[0]:.2f}, x2={solution.x[1]:.2f}")
print(f"最优值: {solution.z:.2f}")

9. 验证与测试

确保算法正确性的关键测试案例：

import unittest

class TestSimplex(unittest.TestCase):
    def test_unbounded(self):
        c = np.array([1, 1])
        A = np.array([[-1, 1], [-1, -1]])
        b = np.array([-1, -2])
        with self.assertRaises(ValueError):
            SimplexSolver(c, A, b).solve()
    
    def test_optimal(self):
        c = np.array([3, 2])
        A = np.array([[1, 2], [1, 1]])
        b = np.array([6, 4])
        solver = SimplexSolver(c, A, b)
        self.assertAlmostEqual(solver.solve().z, 8)

10. 扩展应用

单纯形法不仅适用于标准线性规划，还可扩展应用于：

1. 整数规划松弛 ：作为分支定界法的底层求解器
2. 灵敏度分析 ：研究参数变化对解的影响
3. 对偶理论 ：通过单纯形乘子获取对偶解

class SensitivityAnalyzer:
    def __init__(self, solver):
        self.solver = solver
        
    def analyze(self, param_range):
        """分析目标系数变化对解的影响"""
        results = []
        for delta in param_range:
            modified_c = self.solver.c + delta
            results.append(SimplexSolver(modified_c, self.solver.A, 
                                       self.solver.b).solve())
        return results

通过这个完整的实现过程，你不仅掌握了单纯形法的数学原理，还获得了将其转化为高效代码的实践经验。这种"通过编码理解算法"的方式，往往比单纯的理论学习更能带来深刻认知。