应用运筹学基础：线性规划 (2)

这一节课讲解了求解线性规划问题的方法：单纯形法（simplex method）。

基可行解和最优解的关系

接上一次课，首先是用反证法，证明一下最优解一定可以在基可行解处取到：

假设没有基可行解是最优解，设最优解为 $x = egin{bmatrix} x_1 & x_2 & dots & x_k & 0 & 0 & dots & 0 end{bmatrix}^T$ 不是基可行解，由上一次课的一个引理，$x_1$ 到 $x_k$ 对应的 $A$ 中 $k$ 列 $p_1, p_2, dots, p_k$ 是线性相关的。那么有不全为 0 的 $lambda$，使得 $sum_{i=1}^k lambda_i p_i = 0$。

记辅助向量 $v = egin{bmatrix} lambda_1 & lambda_2 & dots & lambda_k & 0 & 0 & dots & 0 end{bmatrix}^T$，我们仍然构造 $x' = x + epsilon v$ 和 $x'' = x - epsilon v$（熟悉的套路...）。显然有 $x = frac{x'+x''}{2}$，由于 $x$ 是最优解，当然 $x'$ 和 $x''$ 也是最优解。如果存在 $lambda_i < 0$，那么我们取 $epsilon = minlimits_{lambda_i < 0} -frac{x_i}{lambda_i}$，就能把 $x'$ 中 $x_1$ 到 $x_k$ 的某一个变成 0；如果不存在 $lambda_i < 0$，那么我们取 $epsilon = min frac{x_i}{lambda_i}$，就能把 $x''$ 中 $x_1$ 到 $x_k$ 的某一个变成 0。也就是说，非 0 元素少了一个。这样一直构造下去，我们就能不断减少非 0 元素，直到 $x_1$ 到 $x_k$ 在 $A$ 中对应的 $k$ 列线性无关，$x$ 就变成了一个基可行解。而且在构造过程中，最优性没有改变，$x$ 还是最优解。

单纯形法

接下来介绍用单纯形法（simplex method）求解线性规划问题的方式。单纯形法的每一步都在引入一个非基变量取代某一基变量，找出目标函数值更优的另一基本可行解，这样一步一步调整到最优解。

来看一个例题：$$egin{matrix} & maxlimits_{x} quad 3x_1 + 2x_2 \ ext{s.t.} & 2x_1 + x_2 le 12 \ & x_1 + 2x_2 le 9 \ & x_1, x_2 ge 0 end{matrix}$$ 引入松弛变量，将原问题变为 $$egin{matrix} & maxlimits_{x} quad z = 3x_1 + 2x_2 \ ext{s.t.} & 2x_1 + x_2 + x_3 = 12 \ & x_1 + 2x_2 + x_4 = 9 \ & x_1, x_2, x_3, x_4 ge 0 end{matrix}$$ 首先，我们有一个可行解：$x = egin{bmatrix} 0 & 0 & 12 & 9 end{bmatrix}^T$，当前目标函数值为 $z = 0$。我们把各个变量用非基变量进行表示，有 $$x_3 = 12 - 2x_1 - x_2 \ x_4 = 9 - x_1 - 2x_2 \ z = 3x_1 + 2x_2$$ 看起来 $x_1$ 对 $z$ 的贡献比 $x_2$ 来的大（$x_1$ 的系数比较大，这个系数称为“检验数”），我们考虑把 $x_1$ 变成基变量。要注意，$x_3 ge 0$ 和 $x_4 ge 0$ 的条件仍然需要成立。根据 $x_1$ 与 $x_3$ 和 $x_4$ 之间的表达式不难看出，当 $x_1 = 6$ 时，$x_3$ 最先变成 0，我们把它从基变量中去除。

现在，可行解变为 $x = egin{bmatrix}6 & 0 & 0 & 3end{bmatrix}^T$，当前目标函数值为 $z = 18$。继续把各个变量用非基变量进行表示，有 $$x_1 = 6 - x_2/2 - x_3/2 \ x_4 = 3 - 3x_2/2 + x_3/2 \ z = 18 + x_2/2 - 3x_3/2$$ 从系数可以看出，只有增加 $x_2$ 才能对目标函数值的增大有所贡献，考虑把 $x_2$ 变成基变量。从 $x_2$ 与 $x_1$ 和 $x_4$ 的关系式中也不难发现，当 $x_2 = 2$ 时，$x_4$ 最先变成 0，我们把它从基变量中去除。

现在，可行解变为 $x = egin{bmatrix} 5 & 2 & 0 & 0 end{bmatrix}^T$，当前目标函数值为 $z = 19$。继续把各个变量用非基变量进行表示，有 $$x_1 = 5 - 2x_3/3 + x_4/3 \ x_2 = 2 + x_3/3 - 2x_4/3 \ z = 19 - 4x_3/3 - x_4/3$$ 可以发现，$x_3$ 和 $x_4$ 与 $z$ 相关的系数都是负值，那么无论把哪个变量加入基变量，都只能让目标函数值变小了。所以此时我们就得到了线性规划问题的最优解：$x = egin{bmatrix} 5 & 2 & 0 & 0 end{bmatrix}^T$，目标函数值为 $z = 19$。

当然，完全有可能出现某一个变量可以无限增大，求不出最优解的情况。考虑以下线性规划问题：$$egin{matrix} & maxlimits_{x} quad 2x_1 + x_2 \ ext{s.t.} & -3x_1 + x_2 le 3 \ & -4x_1 + x_2 le 5 \ & x_1, x_2 ge 0 end{matrix}$$ 加入松弛变量后有 $$egin{matrix} & maxlimits_{x} quad 2x_1 + x_2 \ ext{s.t.} & -3x_1 + x_2 + x_3 = 3 \ & -4x_1 + x_2 + x_4 = 5 \ & x_1, x_2, x_3, x_4 ge 0 end{matrix}$$ 有初始可行解 $x = egin{bmatrix} 0 & 0 & 3 & 5 end{bmatrix}^T$，仍然用非基变量表示其它变量，有 $$x_3 = 3 + 3x_1 - x_2 \ x_4 = 5 + 4x_1 - x_2 \ z = 2x_1 + x_2$$ 从 $x_1$ 和 $x_3$ 与 $x_4$ 的关系式中可以看出，不管 $x_1$ 如何增大，两个基变量都不会小于 0，而且 $x_1$ 的检验数也是正数，那么这个优化问题没有最优解。

归纳起来，单纯形法的基本步骤如下（来自 wiki）：

1. 把线性规划问题的约束方程组表达成典范型方程组，找出基本可行解作为初始基可行解。
2. 若基本可行解不存在，即约束条件有矛盾，则问题无解。
3. 若基本可行解存在，从初始基可行解作为起点，根据最优性条件和可行性条件，引入非基变量取代某一基变量，找出目标函数值更优的另一基本可行解。
4. 按步骤3进行迭代，直到对应检验数满足最优性条件（这时目标函数值不能再改善），即得到问题的最优解。
5. 若迭代过程中发现问题的目标函数值无界，则终止迭代。

单纯形表

可以用矩阵的形式，把各个变量用非基变量进行表示。

设要优化的线性规划问题（标准形式）为 $$egin{matrix} & maxlimits_{x} quad z = c^Tx \ ext{s.t.} & Ax le b \ & x ge 0 end{matrix}$$ 我们可以把各个矩阵或向量根据基变量和非基变量分为两部分：令 $c^T = egin{bmatrix} c_B^T & c_N^T end{bmatrix}$，令 $A = egin{bmatrix} A_B & A_N end{bmatrix}$，令 $x = egin{bmatrix} x_B^T & x_N^T end{bmatrix}^T$，显然我们有 $$A_Bx_B + A_Nx_N = b \ z = c_B^Tx_B + c_N^Tx_N$$ 通过移项就能用 $x_N$ 表示其它变量：$$x_B = A_B^{-1}b - A_B^{-1}A_Nx_N \ z = c_B^TA_B^{-1}b + (c_N^T - c_B^TA_B^{-1}A_N)x_N$$ 当然啦，由于基可行解中 $x_N = 0$，我们有 $$x_B = A_B^{-1}b \ z = c_B^TA_B^{-1}b$$ 单纯形表可以帮助我们利用矩阵形式，通过列表的方式求解线性规划问题，也利于编程实现。

首先画一张这样的表：$$egin{array}{c|cc|c} & c_B^T & c_N^T & 0 \ hline x_B & A_B & A_N & bend{array}$$ 利用行变换将 $A_B$ 变为 $I$，根据线性代数的知识，这相当于左乘了一个 $A_B^{-1}$，那么表会变成这样：$$egin{array}{c|cc|c} & c_B^T & c_N^T & 0 \ hline x_B & I & A_B^{-1}A_N & A_B^{-1}bend{array}$$ 然后再把下面一行乘上 $c_B^T$，用上面一行减一减，得到 $$egin{array}{c|cc|c} & 0 & c_N^T-c_B^TA_B^{-1}A_N & -c_B^TA_B^{-1}b \ hline x_B & I & A_B^{-1}A_N & A_B^{-1}bend{array}$$ 这是一张非常神奇的表：第一行第二列（$x_B$ 那一列我们不看）就是检验数，第一行第三列就是 $-z$，而基变量和非基变量之间的系数也可以通过第二行的第二、三两列求出，计算起来非常方便。表格每这样计算一次，就是单纯形法里的一次迭代。

虽然这个表格看起来有点麻烦（比如，看起来每次迭代好像都要重新写上 $A_B$ 和 $A_N$ 什么的，从头开始变化？），但实际操作起来是非常简便的，每次迭代的变化不会很多。看一个例子就知道了：

$$egin{matrix} & maxlimits_{x} quad 3x_1 + 2x_2 \ ext{s.t.} & 2x_1 + x_2 le 12 \ & x_1 + 2x_2 le 9 \ & x_1, x_2 ge 0 end{matrix}$$ 这个例子和上一节的例子是一样的。仍然加入松弛变量变化为标准形式：$$egin{matrix} & maxlimits_{x} quad z = 3x_1 + 2x_2 \ ext{s.t.} & 2x_1 + x_2 + x_3 = 12 \ & x_1 + 2x_2 + x_4 = 9 \ & x_1, x_2, x_3, x_4 ge 0 end{matrix}$$ 初始的基可行解是 $x = egin{bmatrix} 0 & 0 & 12 & 9end{bmatrix}^T$，绘制初始的单纯形表：$$egin{array}{c|cccc|c} & 3 & 2 & 0 & 0 & 0 \ hline x_3 & 2 & 1 & 1 & 0 & 12 \ x_4 & 1 & 2 & 0 & 1 & 9 end{array}$$ 由于 $x_B = egin{bmatrix} x_3 & x_4 end{bmatrix} ^ T$，所以此时的 $A_B$ 由第二行和第三行的三、四两列组成，而且恰好是 $I$；$c_B^T$ 由第一行的第三、四两列组成，而且恰好是 0，我们直接来到了最后一步。

根据单纯形法，我们选择检验数中较大的那个（3，对应 $x_1$）。由于 $12/2 = 6 < 9/1 = 9$，所以 $x_1$ 成为基变量，$x_3$ 被移出基变量。修改表格为 $$egin{array}{c|cccc|c} & 3 & 2 & 0 & 0 & 0 \ hline x_1 & 2 & 1 & 1 & 0 & 12 \ x_4 & 1 & 2 & 0 & 1 & 9 end{array}$$ 由于 $x_B = egin{bmatrix} x_1 & x_4 end{bmatrix} ^ T$，所以此时的 $A_B$ 由第二行和第三行的一、四两列组成（有人可能会问：为什么不用原问题里的矩阵 $A$ 呢？因为矩阵进行行变化之后，不改变方程的解，而且每次都用修改后的 $A$，只要更改一列就可以把 $A_B$ 变成 $I$，较为方便）。我们发现，只要把第一列改成单位向量即可，通过行变换，变化后的表格为 $$egin{array}{c|cccc|c} & 3 & 2 & 0 & 0 & 0 \ hline x_1 & 1 & 1/2 & 1/2 & 0 & 6 \ x_4 & 0 & 3/2 & -1/2 & 1 & 3 end{array}$$ 接下来我们把下面的几列乘以 $c_B^T$，注意此时 $c_B^T$ 是由第一行第一、四列组成的，所以我们要计算的是 $$egin{bmatrix}3 & 0end{bmatrix}egin{bmatrix} 1 & 1/2 & 1/2 & 0 & 6 \ 0 & 3/2 & -1/2 & 1 & 3 end{bmatrix}$$ 把第一行减去矩阵计算的结果，我们就有了第二次迭代后的单纯形表：$$egin{array}{c|cccc|c} & 0 & 1/2 & -3/2 & 0 & -18 \ hline x_1 & 1 & 1/2 & 1/2 & 0 & 6 \ x_4 & 0 & 3/2 & -1/2 & 1 & 3 end{array}$$ 依然选择最大的检验数（$1/2$，对应 $x_2$），由于 $3/(3/2) = 2 < 6/(1/2) = 12$，所以 $x_2$ 成为基变量，$x_4$ 被移出基变量。进行和上面类似的过程。$$egin{array}{c|cccc|c} & 3 & 2 & 0 & 0 & 0 \ hline x_1 & 1 & 1/2 & 1/2 & 0 & 6 \ x_2 & 0 & 3/2 & -1/2 & 1 & 3 end{array}$$ $$egin{array}{c|cccc|c} & 3 & 2 & 0 & 0 & 0 \ hline x_1 & 1 & 0 & 2/3 & -1/3 & 5 \ x_2 & 0 & 1 & -1/3 & 2/3 & 2 end{array}$$ 第三次迭代后的单纯形表为：$$egin{array}{c|cccc|c} & 0 & 0 & -4/3 & -1/3 & -19 \ hline x_1 & 1 & 0 & 2/3 & -1/3 & 5 \ x_2 & 0 & 1 & -1/3 & 2/3 & 2 end{array}$$ 此时所有检验数都为非正数，那么单纯形法结束。单纯形表右上角的 $-z$ 就是最优目标函数的相反数，所以最优目标函数为 19。表格最右一列的 $A_B^{-1}b$ 就是 $x_B$ 的取值，所以最优解为 $x_1 = 5$，$x_2 = 2$。

退化

退化是指一个基可行解中，存在至少一个基变量为 0 的情况。也就是说，这个基变量可以和另一个非基变量任意互换，而不影响结果（反正两个变量在这个解里取值都是 0）。

退化会给单纯形法的求解带来一些麻烦：它可能会让单纯形法陷入循环，可能无法找到最优解。这个页面里提供了一个让单纯形法无限循环的例子。

不过，我们可以使用如下法则让单纯形法一定不会陷入循环（也就是一定能找到最优解）：在有多个可以入基的变量时，选择下标最小的一个；在有多个可以出基的变量时，也选择下标最小的一个。不过我不会证明...

单纯形法的证明

接下来证明单纯形法在非退化的情况下为什么可以取到最优解，以及为什么可以停止。

首先证明：若检验数向量 $hat{c} = c_N - c_B^TA_B^{-1}A_N$ 的每一维都小等于 0，那么此时基可行解 $x$ 是最优解。

反证：假设有一个更优的 $y$ 才是最优解。令 $d = x - y$，我们有 $$Ax - Ay = b - b = 0 \ = Ad = A_Bd_B + A_Nd_N$$ 那么 $$d_B = -A_B^{-1}A_Nd_N$$ 则 $$c^Td = c_B^Td_B + c_N^Td_N \ = (c_N^T-c_B^TA_B^{-1}A_N)d_N = hat{c}d_N$$ 要注意，这里的 B（基变量）和 N（非基变量）都是对于 $x$ 而言的。

显然根据基可行解的定义 $x_N = 0$，又 $y ge 0$，所以 $d_N ge 0$；又 $hat{c} le 0$，所以 $c^Td = hat{c}d_N le 0$，那么 $c^Ty = c^Tx + c^Td le c^Tx$，说明 $y$ 并没有比 $x$ 更优，矛盾。这就说明了为什么在检验数均小等于 0 时停止算法，就能得到最优解。

类似地，我们可以说明若基可行解 $x$ 是非退化的最优解，那么检验数向量的每一维都小等于 0。否则由于 $x$ 是非退化的最优解，如果有一个检验数大于 0，它对应的非基变量一定有增加的空间（而不会像退化的解一样，增加的空间为 0），那么就能构造一个更优的解。这就说明了，在非退化的情况下，肯定有解满足算法停止的情况。

最后简单说明非退化情况下的单纯形法一定可以停止。因为非退化的单纯形法的每一次迭代都会让答案更优一点，所以它访问的基可行解都是不会重复的。而基可行解的数量是有限的（根据上一节课，至多 $C_n^m$ 个），其中又存在着让算法停止的最优解，所以算法一定可以停止。单纯形法的最差复杂度是指数级的（这个最差情况由 Klee 和 Minty 提出，是一个高维立方体，详见维基百科），不过在大多数问题下，它的运行效率都还是比较快的。