线性规划：单纯形算法

考虑线性规划的标准形式（更多介绍参考《线性规划的标准形》）：
$$\begin{array}{rl}\min & c^{T}x\\ \text{s.t.}& Ax=b\\ & x\ge 0\end{array}$$
其中 $c,x\in {\mathbb{R}}^n$，$A\in {\mathbb{R}}^{m\times n}$，$b\in {\mathbb{R}}^m\ge 0$。

单纯形算法的思路是从多面体的一个顶点出发，然后沿着降低目标的方向，迭代到另一个顶点，直到目标值无法降低，于是得到最优解。

基本可行解

把 $A$ 拆成两个部分：
$$A=[BN]$$
其中 $B\in {\mathbb{R}}^{m\times m}$， $N\in {\mathbb{R}}^{m\times (n?m)}$。

相应地，把 $c$ 和 $x$ 分别拆成两个部分，即
$$c=\left[\begin{array}{c}{c}_B\\ c_N\end{array}\right]x=\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right].$$
注意：$c_B,x_B$ 的下标 与 $B$ 的列分别对应。

假设 $A$ 满秩，我们可以交换列的位置保证 $B$ 可逆，所以 $B$ 的列向量是一组 基。$B$ 称为 基矩阵，$N$ 称为 非基矩阵。 把 $x_B$ 对应的变量称为 基变量， $x_N$ 对应的变量称为 非基变量。

这样一来，约束条件可以写成
$$[BN]\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right]=b?B{x}_B+N{x}_N=b.$$
我们有
$$x_B=B^{?1}b?B^{?1}N{x}_N.$$
令$x_N=0$，得到 基本解
$$x_0=\left[\begin{array}{c}{B}^{?1}b\\ 0\end{array}\right].$$
换句话说，基本解 $x_0$ 满足约束条件 $A{x}_0=b$，但有可能不可行。如果 $x_0\ge 0$，那么称它为 基本可行解。换句话说，基本可行解是多面体的一个顶点（证明略）。

假设已知一个初始的基本可行解，接下来我们要想办法迭代到另一个基本可行解，并且使得目标函数值下降。

迭代条件

要想目标函数降低，我们自然想到梯度下降法。令 $f(x)$ 代表目标函数，梯度 $?f(x)$ 就是函数值增加的方向。所以应该沿着负梯度的方向移动 $x$。

把 $x_B$ 代入目标函数 $f(x):=c^{T}x$ 可得
$$\begin{array}{rl}f(x)& =c^{T}x\\ & =c_B^T{x}_B+c_N^T{x}_N\\ & =c_B^T{B}^{?1}b+(c_N^T?c_B^T{B}^{?1}N)x_N.\end{array}$$
对 $f(x)$ 关于 $b$ 求导， 得到 Shadow Price （也称为 对偶变量）：
$${\lambda }^T:=c_B^T{B}^{?1}.$$
对 $f(x)$ 关于 $x_N$ 求导（求梯度），得到 Reduced Cost ：
$${\mu }^T:=c_N^T?c_B^T{B}^{?1}N=c_N^T?{\lambda }^{T}N.$$
令 $J$ 代表非基变量的下标。$?j\in J$，$x_j$ 每增加一个单位，目标函数就增加 ${\mu }_j$。

如果存在 ${\mu }_j<0$，我们只要把 $x_j$ 从 0 增加 $\delta$，那么目标函数就可以降低 $?{\mu }_j\delta$。

当 ${\mu }_j\ge 0$，$?j\in J$，这意味着目标函数值无法降低，此时达到最优解（证明略）。

如何迭代

假设存在 $j\in J$ 使得 ${\mu }_j<0$，于是我们想增加 $x_j$ 使得目标函数降低。但是问题来了，增加多少可以保证 $x$ 可行？

回顾基本可行解的定义
$$x=\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right]=\left[\begin{array}{c}{B}^{?1}b?B^{?1}N{x}_N\\ x_N\end{array}\right].$$
注意到 $x_N=0$，我们只要保证 $x_B\ge 0$。接下来把 $x_B$ 换一种写法，令 $\tilde{b}=B^{?1}b$，${\tilde{a}}_j=B^{?1}{a}_j$，我们有
$$x_B=?\begin{array}{c}{x}_{B_1}\\ x_{B_2}\\ ?\\ x_{B_m}\end{array}?=?\begin{array}{c}{\tilde{b}}_1\\ {\tilde{b}}_1\\ ?\\ {\tilde{b}}_m\end{array}???\begin{array}{c}{\tilde{a}}_{j_1}\\ {\tilde{a}}_{j_2}\\ ?\\ {\tilde{a}}_{j_m}\end{array}?x_j=?\begin{array}{c}{\tilde{b}}_1?{\tilde{a}}_{j_1}{x}_j\\ {\tilde{b}}_2?{\tilde{a}}_{j_2}{x}_j\\ ?\\ {\tilde{b}}_m?{\tilde{a}}_{j_m}{x}_j\end{array}?.$$
现在要增加 $x_j$， 那么只要保证 $x_B$ 的每个分量非负即可。

注意到 $x_0$ 是基本可行解，所以 $\tilde{b}=B^{?1}b\ge 0$。容易验证，下面的取值（称之为 Minimum Ratio Test）可以保证 $x_B\ge 0$。
$$x_j:=\min \left\{\frac{{\tilde{b}}_i}{{\tilde{a}}_{j_i}}\text{ and }{\tilde{a}}_{j_i}>0,i=1,2,\dots ,m\right\}.$$
这样一来，$x_j\ge 0$，相应地，另一个变量 $x_{B_i}={\tilde{b}}_i?{\tilde{a}}_{j_i}=0$， 其中 $i$ 是上式达到最小值的下标。换句话说，此时 $x_j$ 成为基变量（入基），而 $x_{B_i}$ 成为非基变量（出基）。可以证明，新得到的解仍然是一个基本可行解。

注意：如果 ${\tilde{a}}_{j_i}\le 0$，$?i$ ，这意味着 $x_j$ 可以无限大，最优目标函数值为 $?\infty$，于是最优解不存在。

算法描述

第0步：输入基本可行解对应的基矩阵 $B$。

第1步：判断当前的解是否最优。如果 $\mu \ge 0$，当前是最优解，算法停止。

第2步：计算入基变量和出基变量。如果存在 $j\in J$ 使得 ${\mu }_j<0$，那么 $x_j$ 是入基变量。执行 Minimum Ratio Test， 找到出基变量 $x_i$。

第3步：判断问题是否无界。如果 ${\tilde{a}}_{j_i}\le 0$，$?i$，则代表无界，算法停止。

第4步：执行出入基操作，更新基矩阵 $B$，然后执行第1步。

算法实现

下面我们用Python来实现单纯形算法。

先定义算法的输入和输出。

class SimplexA(object):"""单纯形算法（基本版）。Note:1、系数矩阵满秩。2、未处理退化情形。3、输入基本可行解（对应的列）。"""def __init__(self, c, A, b, v0):""":param c: n * 1 vector:param A: m * n matrix:param b: m * 1 vector:param v0: basic variables, list of variable indices注意：v0是 B 的列下标。x0 = B^{-1}b 即为基本可行解（需要保证x0非负）。"""# 输入self._c = np.array(c)self._A = np.array(A)self._b = np.array(b)self._basic_vars = v0  # basic variablesself._m = len(A)self._n = len(c)self._non_basic_vars = self._init_non_basic_vars()  # non basic variables# 辅助变量self._iter_num = 0self._B_inv = None  # inverse of Bself._lambda = None  # shadow priceself._mu = None  # reduced cost# 输出self._obj = None  # objective function valueself._sol = None  # solutionself._status = None

接下来要实现单纯形算法 SimplexA.solve()，思路如下。

class SimplexA(object):# ...# 其它函数省略……def solve(self):self._iter_num = 0  # 记录迭代次数self._check_init_solution()  # 检查初始基本解是否可行self._update_reduced_cost()self._update_obj()self._update_solution()while not self._is_optimal():  # 判断是否最优或者无界if self._status == "UNBOUNDED":breakself._pivot()  # 迭代（选主元入基，执行Minimum Ratio Test，然后出基）self._update_reduced_cost()  # 更新Reduced Cost: mu = c_N - lambda * Nself._update_obj()self._update_solution()self._iter_num += 1if self._status != "UNBOUNDED":self._status = 'OPTIMAL'print("Done >> status: {}".format(self._status))return self

目标函数值和Reduced Cost的计算可以直接套公式。关键是实现每次迭代的出入基操作，即 SimplexA._pivot()。

class SimplexA(object):# ...# 其它函数省略……def _pivot(self):""" 选主元，入基和出基。"""j_ind = np.argmin(self._mu)# 入基变量 x_jj = self._non_basic_vars[j_ind]# 出基变量 x_ii = self._minimum_ratio_test(j)if i is None:self._status = 'UNBOUNDED'return# update basic varsfor k in range(self._m):if self._basic_vars[k] == i:self._basic_vars[k] = jbreak# update non basic varsself._non_basic_vars[j_ind] = idef _minimum_ratio_test(self, j):""" Minimum Ratio Test.给定入基的非基变量，返回出基的基变量。:param j: 入基变量 x_j 的下标 j:return: 出基变量 x_i 的下标 i"""b_bar = np.dot(self._B_inv, self._b)a_in = np.dot(self._B_inv, self._A[:, j])ratios = list(map(lambda b, a: b/a if a > 1e-6 else np.infty, b_bar, a_in))i_ind = np.argmin(ratios)if ratios[i_ind] != np.infty:return self._basic_vars[i_ind]else:return None

完整代码

遗留问题

需要注意的是，前面描述的算法并没有完全解决问题。比如：初始的基本可行解如何找到？矩阵 $A$ 不是满秩如何处理？

此外，迭代过程有没有可能出现死循环？比如初始解为 $x_0$， 经过一系列迭代又回到 $x_0$。遗憾的是，这种情况有可能发生，称为 退化情形。

矩阵 $A$ 不满秩时意味着有冗余的约束。如果已知初始的基本可行解，于是已知 $A$ 的基，那么去掉冗余的约束即可。综上所述，还有两个问题有待解决：

1. 计算一个初始基本可行解；
2. 处理退化情形。

我们将在后续的教程中介绍处理的方法。

参考文献

[1] Christopher Griffin. Linear Programming: Penn State Math 484 Lecture Notes. Version 1.8.3. Chapter 5 （点此下载）