·完整·单纯形算法(Simplex Algorithm)，附C源码

前段时间参加了华为的2017软件精英挑战赛，用到了单纯形算法求解线性规划问题，学习了正单纯形，对偶单纯形以及割平面法并用C语言实现了完整的simplex算法。

单纯形算法是用来求解线性规划问题的，其被用在了众多SMT Solver中，如Yices, Z3等，都是基于的单纯形算法，其算法效率在最坏情况下为指数级别，不过实现证明其效率在大多数情况下都是令人满意的。

标准形式

$\max{\{b+\sum_{1\le k \le n}c_kx_k\}}$

$\sum_{1\le k \le n}A_1x_k\le b_1$
…
$\sum_{1\le k \le n}A_mx_k\le b_m$
$x_1,...,x_n\ge 0$

对于$b_1,...,b_m\ge0$的情况，可使用正单纯形算法求解,对于存在$b_i\lt0$的情况，则需要使用对偶单纯形求解。这里求出来的解都是非整数解，若要得到目标式的整数最优解，在前面的基础上还需要使用割平面法，即引入Gomory(高莫雷)约束。

一些特殊情况的处理

1.对于不满足$x_i\ge0$的情况，可将$x_i$化为两个非负数的差值，即
$x_i=x_j-x_k$, 其中$x_j,x_k\ge0$
2.如果限制式子为等式，可将其化为一个大于等于和小于等于，如
$a+b=c \Leftrightarrow (a+b\le c \wedge a+b\ge c)$
3.如果要求的目标式$Z$的最小值，则可转化为先求$-Z$的最大值

正单纯形算法

首先将标准形式写为如下形式(Slack Form)
$Z=b+\sum_{1\le k \le n}c_kx_k$

$y_1=b_1-\sum_{1\le k \le n}A_1x_k$
…
$y_m=b_m-\sum_{1\le k \le n}A_mx_k$

这里引入了新的松弛变量$y_1,...,y_m\ge0$,不难证明，新的形式跟原不等式等价。$x_k$称为basic variable(基变量)， $y_m$称为non-basic variable(非基变量)。因为$b\ge0$,所以所有的非基变量为0为一组可行解，我们就是从这个初始可行解触发，一步一步找到最优解。下面是具体的步骤：
1.查找$c_k$且$c_k>0$,若不存在(目标式已经无法再继续增大)，跳转到第4步
2.找到$\min{\dfrac{b_i}{A_i}}$,即限制条件最强的一组式子
3.交换$x_k$和$y_i$, $x_k$变为新的基变量，$y_i$变为新的非基变量，将其他式子中的$x_k$替换为新的值，然后回到第1步
4.将所有的非基变量取0值，即可算出所有的基变量的值以及目标式的最大值，算法结束

单纯形的物理意义

如上图所示，为7个限制不等式所构成的凸多面体，每个限制条件可看成正凸多面体的一个面，凸多面体的内部(包括表面)即为可行解区域。可以证明，最优解如果存在，则一定在凸多面体的顶点上。对于三维空间来说(更高维空间可类似推广)，某个顶点比为三个或以上的平面的交点，这个交点满足对应相交平面的限制不等式。如点A，同时满足式子②③⑦。
算法中第3步中的交换因子(pivot)，可以看成是从凸多面体的一个顶点走向另一个顶点。

对偶单纯形

在标准形式(Slack Form)中，如果$b_m<0$，则不能直接用正单纯形求解，需要使用对偶单纯形，先将标准型化为正单纯形的标准形式，即$b_m\ge0$，然后再使用正单纯形进行求解。
对偶单纯形的求解步骤如下：
1.引入新的松弛变量$z\ge0$，得到如下新的标准型
$Z=b+\sum_{1\le k \le n}c_kx_k-z$
$y_1=b_1-\sum_{1\le k \le n}A_1x_k+z$
…
$y_m=b_m-\sum_{1\le k \le n}A_mx_k+z$
2.取$\min{\{b_i\}}, 0\le i\le m$，交换$z$与$y_i$，$z$变为基变量，$y_i$变为非基变量。执行完此步骤后，新的所有常数$b$都将变为非负数
3.为了使得加入$z$后的标准型跟原标准型等价，需要满足$z=0$的条件。此时$z$已经是基变量，所以我们需要先将其从基变量中置换出来，即变为非基变量。如果无法将$z$变为非基变量，则原式无解
4.当$z$再次从基变量置换为非基变量后，此时$z$可以取最小值0。然后我们可将$z$忽略，使用正单纯形算法求解目标值的最大值

割平面法求整数最优解

以上求得的都是非整数最优解，即LP(Linear Programing)，如果要求整数最优解(Interger Linear Programing)，则需要在得到非整数最优解的基础上，对其进行切割，以得到整数解，这里采用割平面法。
设求出的最优解:
$X=(b_1,...,b_m)$
设$x_i$不为整数，
$x_i=b_i-\sum_{1\le k \le n}A_kx_k$, $x_i$为基变量，$x_k$为非基变量
将$b_i$与$A_k$分离成一个整数与真分数之和，即
$b_i=[b_i]+f_i$, $A_k=[A_k]+f_k$
其中$[b_i],[A_k]$为整数，$f_i,f_k$为真分数，则有
$x_i=[b_i]+f_i-\sum_{1\le k \le n}[A_k]x_k-\sum_{1\le k \le n}f_kx_k$
把整数项移到左边，分数项移到右边，得
$x_i-[b_i]+\sum_{1\le k \le n}[A_k]x_k=f_i-\sum_{1\le k \le n}f_kx_k$
因为等式两边都为整数，且
$f_i-\sum_{1\le k \le n}f_kx_k\le f_i<1$
则 $f_i-\sum_{1\le k \le n}f_kx_k\le0$, 这就是高莫雷(Gomory)约束。
对上式引入一个新的松弛变量$g_i$得
$g_i=-f_i+\sum_{1\le k \le n}f_kx_k$
将其添加到原式中，使用对偶单纯形进行求解。如果算出解仍不为整数解，则选取新的Gomory约束，继续进行切割

用C实现的源码地址如下:
(https://github.com/zjc666/Simplex_Algorithm)