求解线性规划问题的单纯形算法

1947年，丹齐格提出了一种求解线性规划问题的方法，即今天所称的单纯形法，这是一种简洁且高效的算法，被誉为20世纪对科学发展和工程实践影响最大的十大算法之一。 上文提到线性规划问题的最优解一定是基本可行解，单纯形法的思路即在不同的基向量下求不同的基本可行解，然后找到最优的解。从几何的角度来看，也就是从一个极点转换到另一个极点，直至找到最优极点的过程。 那么这样的话可以把算法分成三个子问题，1.如何从一个基本可行解转换到另一个基本可行解。2.如何确定应该转移到哪个极点。3.什么时候停止转移操作，即如何判断当前基本可行解是否为最优解。

极点的转移

线性规划的标准型为

$$\begin{align*} minimize\quad c^{T}x\\ subject\ to\quad Ax=b\\ x\ge0 \end{align*}$$

考虑到方程$Ax=b$，展开成如下规范型的方程组

$$\begin{align*} x_{1}+y_{1\ m+1}x_{m+1}+...+y_{1n}x_{n}=y_{10}\\ x_{2}+y_{2\ m+1}x_{m+1}+...+y_{2n}x_{n}=y_{20}\\ \vdots\\ x_{m}+y_{m\ m+1}x_{m+1}+...+y_{mn}x_{n}=y_{m0} \end{align*}$$

可以将该方程组转换为$[I_{m},Y_{m,n-m}]x=y_{0}$，这种形式的方程组$Ax=b$称为典式，方程组的典式与原方程组的解是相同的，在典式表达式中，与基列向量对应的变量为基变量，其他的变量为非基变量，即在方程$[I_{m},Y_{m,n-m}]x=y_{0}$中，$x_{1},x_{2}...,x_{m}$为基变量，其他的变量是非基变量。考虑增广矩阵规范型$[I_{m},Y_{m,n-m},y_{0}]$，其最后一列的各元素是向量$b$关于基{$a_{1},...,a_{m}$}的坐标。

现在考虑增广矩阵的更新，即用某个非基变量替换某个基变量，求新的基变量对应的典式表达式，比如用非基变量$a_{q},m<q\le n$替换基变量$a_{p},1\le p\le m$。在原矩阵上，$a_{q}$可以表示为

$$a_{q}=\sum_{i=1}^{m}y_{iq}a_{i}=\sum_{i=1\\i\ne p}^{m}y_{iq}a_{i}+y_{pq}a_{p}$$

注意仅当$y_{pq}\ne 0$时，向量组{$a_{1},...,a_{p-1},a_{q},a_{p+1},...,a_{m}$}才是线性无关的，才能形成一组基。 可以由上述方程求解出

$$a_{p}=\frac{1}{y_{pq}}a_{q}-\sum_{i=1\\i\ne p}^{m}\frac{y_{iq}}{y_{pq}}a_{i}$$

利用原来的增广矩阵，可以将任意向量$a_{j}(m<j\le n)$表示为

$$a_{j}=y_{1j}a_{1}+y_{2j}a_{2}+...+y_{mj}a_{m}$$

将$a_{p}$的表达式代入可以得到

$$a_{j}=\sum_{i=1\\i\ne p}^{m}(y_{ij}-\frac {y_{pj}}{y_{pq}}y_{iq})a_{i}+\frac{y_{pj}}{y_{pq}}a_{q}$$

利用$y_{ij}^{'}$表示新的增广矩阵中的元素，由上式可得

$$\begin{align*} y_{ij}^{'}=y_{ij}-\frac{y_{pj}}{y_{pq}}y_{iq},i\ne p\\ y_{pj}^{'}=\frac{y_{pj}}{y_{pq}} \end{align*}$$

利用以上公式对矩阵的变换称为关于元素$(p,q)$的枢轴变换，以上即为从一个极点转换到另一个极点的公式。

确定转移的极点

确定转移的极点即将向量组转换到另一组坐标系上，详细来说就是将向量组中m个向量之一取出，从另外n-m个向量中选择一个向量加入到基向量组形成新的基。我们记为将$a_{p}$向量取出（出基），将$a_{q}$向量加入（入基），又可以将这个问题分为两部分，即分别确定$p$和$q$。 对于$q$的确定，先看枢轴变换的公式如下

$$\begin{align*} y_{ij}^{'}=y_{ij}-\frac{y_{pj}}{y_{pq}}y_{iq},i\ne p\\ y_{pj}^{'}=\frac{y_{pj}}{y_{pq}} \end{align*}$$

再看另一种方法将一个非基变量$a_{q}$变成基向量 利用当前的基向量来表示向量$a_{q}$

$$a_{q}=y_{1q}a_{1}+y_{2q}a_{2}+...+y_{mq}a_{m}$$

上式两边同时乘上$\epsilon>0$可以得到

$$\epsilon a_{q}=\epsilon y_{1q}a_{1}+\epsilon y_{2q}a_{2}+...+\epsilon y_{mq}a_{m}$$

将基本解$x=[y_{10},...,y_{m0},0,...,0]^T$代入方程$Ax=b$可以得到

$$y_{10}a_{1}+...+y_{m0}a_{m}=b$$

联立上面两个等式，可以得到

$$(y_{10}-\epsilon y_{1q})a_{1}+(y_{20}-\epsilon y_{2q})a_{2}+...+(y_{m0}-\epsilon y_{mq})a_{m}+\epsilon a_{q}=b$$

即向量$[y_{10}-\epsilon y_{1q},...,y_{m0}-\epsilon y_{mq},0,...,\epsilon,...,0]$是方程$Ax=b$的一个解，当$\epsilon$不断增大，向量的前m个元素中首次出现0，此时可以得到一个基本可行解，变换过后的基本可行解对应的目标函数值为

$$\begin{align*} z=c_{1}(y_{10}-y_{1q}\epsilon)+...+c_{m}(y_{m0}-y_{mq}\epsilon)+c_{q}\epsilon\\ =z_{0}+[c_{q}-(c_{1}y_{1q}+...+c_{m}y_{mq})]\epsilon \end{align*}$$

其中，$z_{0}=c_{1}y_{10}+...+c_{m}y_{m0}$，令$z_{q}=c_{1}y_{1q}+...+c_{m}y_{mq}$，可以得到$z=z_{0}+(c_{q}-z_{q})\epsilon$，当$z-z_{0}=(c_{q}-z_{q})\epsilon<0$时，说明基本可行解对应的目标函数值变小了，即认为这个向量加入基向量组后是有益的。

那么可以对每一个$q,m<q\le n$计算$c_{q}-z_{q}$，令$r_{i}=0(i=1,...,m),r_{i}=c_{i}-z_{i}(i=m+1,...,n)$，称$r_{i}$为第$i$个简化价值系数，也称为检验数。可以通过证明得到，当且仅当相应的检验数都是非负的时候，该基本可行解是最优解。

即可以通过计算检验数的方法来判断当前解是否为最优解，当检验数中有负数的时候，则说明当前解不是最优解，可以通过将第$q$个向量加入基向量组的方法减小目标函数值，如果有多个负数，则选择检验数最小的第$q$个向量加入基向量组。

上面确定了$q$的选择方法，下面来介绍一下如何确定出基向量$a_{p}$，由上文可知，对于向量$[y_{10}-\epsilon y_{1q},...,y_{m0}-\epsilon y_{mq},0,...,\epsilon,...,0]$，当$\epsilon$不断增大，在前$m$个元素中会有第$i$个元素$y_{i0}-\epsilon y_{iq}$最先等于0，那么我们则选择$a_{i}$作为出基向量$a_{p}$，即$p=arg\ min_{i}\{y_{i0}/y_{iq}:y_{iq}>0\}$，需要注意的是，如果不存在$y_{iq}>0$，则问题有无界解，停止运算，另外如果同时出现多个0元素，那么得到退化的基本可行解，此时选择最小的$p$值。

##确定停止条件 上文介绍了如何选择转移的极点数，在确定$q$的时候，计算检验数，若所有的检验数都是非负，那么此时的基本可行解是最优解，停止计算。另外在确定$p$的时候，若随着$\epsilon$不断增大，向量的前m个元素也随之增大，即对于$0<i\le m,y_{iq}<0$，那么该问题有无界解，也停止计算。

##单纯形算法 上文介绍了单纯形算法几个子问题的处理方法，下面来总结一下单纯形算法 **1.**根据初始基本可行解构造增广矩阵规范型； **2.**计算非基变量的检验数； **3.**如果对于所有$j$都有$r_{j}\le 0$，则停止计算，当前基本可行解是最优解，否则进入下一步 **4.**在小于0的检验数中选择检验数最小的$q$ **5.**如果不存在$y_{iq}>0$，则停止计算，问题有无界解；否则计算$p=arg\ min_{i}\{y_{i0}/y_{iq}:y_{iq>0}\}$，如果得到多个满足条件的下标$i$，令$p$为最小的下标值。 **6.**以元素$(p,q)$作为枢轴元素，进行枢轴变换，更新增广矩阵规范型。 **7.**转到步骤2

另外对于单纯形算法，还有其矩阵表达式，对于寻找初始基本可行解所遇到的困难，又提出了两阶段单纯形法，为了减小计算量，又提出了修正单纯形法等等。

To be continue…