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

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…

技术

2018 · 05 · 04