线性规划：单纯形算法之处理退化

上一篇文章介绍了单纯形算法（《线性规划：单纯形算法》），但是还有一些遗留问题没有解决，比如退化情形（Degeneracy）。

本文介绍如何处理退化情形。

退化情形

从几何上看上看，造成退化情形的原因是顶点重合。下图是三维空间中非退化情形下的多面体，每个顶点由三个面交汇。

当我们把线段 a b ab ab缩短至一点，于是得到退化情形。

如上图所示， a a a和 b b b重合时四个面交于一点，于是迭代过程可能在 a a a 和 b b b 之间转圈（Cycling）。

处理方法

回顾单纯形算法，根据 Minimum Ratio Test 计算入基变量 x j x_j xj 的值。

x j : = min ⁡ { b ~ i a ~ j i and a ~ j i > 0 , i = 1 , 2 , … , m } . x_j := \min \left\{ \frac{\tilde{b}_i}{\tilde{a}_{j_i}} \text{ and } \tilde{a}_{j_i}>0, \quad i=1,2,\ldots, m\right\}. xj:=min{a~jib~i and a~ji>0,i=1,2,…,m}.
其中 b ~ = B − 1 b \tilde{b} = B^{-1}b b~=B−1b， a ~ j = B − 1 a j \tilde{a}_j = B^{-1}a_j a~j=B−1aj， a j a_j aj 代表矩阵 A A A 的列向量。

令 I 0 I_0 I0 代表上式达到最小值的下标集合。
I 0 = { r : b ~ r a ~ j r = min ⁡ [ b ~ i a ~ j i and a ~ j i > 0 , i = 1 , 2 , … , m ] } . I_0 = \left\{r:\frac{\tilde{b}_r}{\tilde{a}_{j_r}} = \min \left[ \frac{\tilde{b}_i}{\tilde{a}_{j_i}} \text{ and } \tilde{a}_{j_i}>0, \quad i=1,2,\ldots, m\right]\right\}. I0={r:a~jrb~r=min[a~jib~i and a~ji>0,i=1,2,…,m]}.

如果 I 0 I_0 I0 只包含一个元素 r r r，那么出基变量为 x B r x_{B_r} xBr，不会出现退化情形。否则可能出现退化情形，即变量 x B r x_{B_r} xBr 先出基，迭代数次又入基，这样一直循环下去导致算法死循环。

解决思路是防止同一个变量反复出基和入基。当 Minimum Ratio Test 发现 I 0 I_0 I0 包含多个出基变量时，用一些规则防止它们反复出基和入基。

下面介绍 Lexicographic Minimum Ratio Test。

在退化的情形下， I 0 I_0 I0 包含了多个元素，这此时我们定义新的下标集合 I 1 I_1 I1：
I 1 = { r : a ~ 1 r a ~ j r = min ⁡ [ a ~ 1 i a ~ j i : i ∈ I 0 ] } . I_1 =\left\{r:\frac{\tilde{a}_{1_r}}{\tilde{a}_{j_r}} = \min \left[ \frac{\tilde{a}_{1_i}}{\tilde{a}_{j_i}}: i \in I_0\right]\right\}. I1={r:a~jra~1r=min[a~jia~1i:i∈I0]}.
注意：上述分式中的分子 a ~ 1 = B − 1 a 1 \tilde{a}_{1} = B^{-1} a_1 a~1=B−1a1，而 a ~ 1 i \tilde{a}_{1_i} a~1i 代表 a ~ 1 \tilde{a}_{1} a~1 中下标为 i i i 的分量。

如果 I 1 I_1 I1 只包含一个元素 r r r，那么交换 x j x_j xj 与 x B r x_{B_r} xBr 从而得到新的基。否则，我们定义类似的集合 I 2 , I 3 I_2, I_3 I2,I3，直到 I k I_k Ik，使得 I k I_k Ik 只包含一个元素。
I k = { r : a ~ k r a ~ j r = min ⁡ [ a ~ k i a ~ j i : i ∈ I k − 1 ] } . I_k =\left\{r:\frac{\tilde{a}_{k_r}}{\tilde{a}_{j_r}} = \min \left[ \frac{\tilde{a}_{k_i}}{\tilde{a}_{j_i}}: i \in I_{k-1}\right]\right\}. Ik={r:a~jra~kr=min[a~jia~ki:i∈Ik−1]}.
可以证明，这样的 I k I_k Ik 始终存在。

设 r r r 是 I k I_k Ik 中唯一的元素，那么交换 x j x_j xj 与 x B r x_{B_r} xBr 从而得到新的基。

通过这种方式，单纯形算法可以在有限步内返回最优解（证明略）。

算法实现

下面用Python实现 Lexicographic Minimum Ratio Test。

单纯形算法的实现这里不重复展示，可以参考上一篇文章。只需要继承之前的类SimplexA，然后重写方法 _minimum_ratio_test(self, j) 即可。

class SimplexAD(SimplexA):"""单纯形算法：处理退化情形。"""def _minimum_ratio_test(self, j):""" Lexicographically minimum ratio test.:param j: 入基变量 x_j 的下标 j:return: 出基变量 x_i 的下标 i"""a_in = np.dot(self._B_inv, self._A[:, j])d0 = self._get_d0(a_in)  # 计算 I_0 的数组下标if d0 is None:return None# The index of "basic_vars".# The associated basic variable will be moved out.i_ind = self._get_out_ind(d0, 0, a_in)return self._basic_vars[i_ind]def _get_d0(self, a_in):""" 根据 Minimum Ratio Test 计算 I_0.实际上计算 d0 即可（定义如下）：1、计算 ratios = [r_1, r_2, ..., r_m].2、计算 d0 = arg min(ratios)，即 ratios 中最小值的下标3、I_0 = [self._basic_vars[i] for i in d0]"""b_bar = np.dot(self._B_inv, self._b)ratios = list(map(lambda b, a: b / a if a > 1e-6 else np.infty, b_bar, a_in))d0 = arg_min(ratios)if ratios[d0[0]] == np.infty:# a_in 的分量 <= 0# 最优目标函数值无界return Nonereturn d0def _get_out_ind(self, d0, it, a_in):""" 用递归的方法计算出基变量（在self._basic_vars[]中的index）.:param d0: I_0 （在self._basic_vars[]的indices）.:param it: iteration number, equals 0 at the beginning.:param a_in: B_inv * A[:, j]:return: the index to "basic_vars" (to be moved out)"""if len(d0) == 1:return d0[0]a_it = np.dot(self._B_inv, self._A[:, it])ratios = [a_it[k] / a_in[k] for k in d0]indices = arg_min(ratios)d1 = [d0[k] for k in indices]return self._get_out_ind(d1, it + 1, a_in)

完整代码

参考文献

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