[最优化]求解线性规划问题的单纯形算法

单纯形算法

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

极点的转移

线性规划的标准型为

minimizecTxsubject toAx=bx≥0$$\begin{gathered} minimize{c}^{T}x \\ subjecttoAx=b \\ x\ge 0 \end{gathered}$$

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

x1+y1 m+1xm+1+...+y1nxn=y10x2+y2 m+1xm+1+...+y2nxn=y20⋮xm+ym m+1xm+1+...+ymnxn=ym0$$\begin{gathered} x_1+y_{1m+1}{x}_{m+1}+...+y_{1n}{x}_n=y_{10} \\ {x}_2+y_{2m+1}{x}_{m+1}+...+y_{2n}{x}_n=y_{20} \\ ⋮ \\ {x}_m+y_{mm+1}{x}_{m+1}+...+y_{mn}{x}_n=y_{m0} \end{gathered}$$

可以将该方程组转换为[Im,Ym,n−m]x=y0$[I_m,Y_{m,n-m}]x=y_0$，这种形式的方程组Ax=b$Ax=b$称为典式，方程组的典式与原方程组的解是相同的，在典式表达式中，与基列向量对应的变量为基变量，其他的变量为非基变量，即在方程[Im,Ym,n−m]x=y0$[I_m,Y_{m,n-m}]x=y_0$中，x1,x2...,xm$x_1,x_2...,x_m$为基变量，其他的变量是非基变量。考虑增广矩阵规范型[Im,Ym,n−m,y0]$[I_m,Y_{m,n-m},y_0]$，其最后一列的各元素是向量b$b$关于基{a1,...,am$a_1,...,a_m$}的坐标。

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

aq=∑i=1myiqai=∑i=1i≠pmyiqai+ypqap$$\begin{gathered} 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 \end{gathered}$$

注意仅当ypq≠0$y_{pq}\ne 0$时，向量组{a1,...,ap−1,aq,ap+1,...,am$a_1,...,a_{p-1},a_q,a_{p+1},...,a_m$}才是线性无关的，才能形成一组基。
可以由上述方程求解出

ap=1ypqaq−∑i=1i≠pmyiqypqai$$\begin{gathered} a_p=\frac{1}{y_{pq}}{a}_q-\sum _{i=1 \\ i\ne p}^m\frac{y_{iq}}{y_{pq}}{a}_i \end{gathered}$$

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

aj=y1ja1+y2ja2+...+ymjam$$a_j=y_{1j}{a}_1+y_{2j}{a}_2+...+y_{mj}{a}_m$$

将ap$a_p$的表达式代入可以得到

aj=∑i=1i≠pm(yij−ypjypqyiq)ai+ypjypqaq$$\begin{gathered} 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 \end{gathered}$$

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

y′ij=yij−ypjypqyiq,i≠py′pj=ypjypq$$\begin{gathered} y_{ij}^{{}^{\prime }}=y_{ij}-\frac{y_{pj}}{y_{pq}}{y}_{iq},i\ne p \\ {y}_{pj}^{{}^{\prime }}=\frac{y_{pj}}{y_{pq}} \end{gathered}$$

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

确定转移的极点

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

y′ij=yij−ypjypqyiq,i≠py′pj=ypjypq$$\begin{gathered} y_{ij}^{{}^{\prime }}=y_{ij}-\frac{y_{pj}}{y_{pq}}{y}_{iq},i\ne p \\ {y}_{pj}^{{}^{\prime }}=\frac{y_{pj}}{y_{pq}} \end{gathered}$$

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

aq=y1qa1+y2qa2+...+ymqam$$a_q=y_{1q}{a}_1+y_{2q}{a}_2+...+y_{mq}{a}_m$$

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

ϵaq=ϵy1qa1+ϵy2qa2+...+ϵymqam$$ϵa_q=ϵy_{1q}{a}_1+ϵy_{2q}{a}_2+...+ϵy_{mq}{a}_m$$

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

y10a1+...+ym0am=b$$y_{10}{a}_1+...+y_{m0}{a}_m=b$$

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

(y10−ϵy1q)a1+(y20−ϵy2q)a2+...+(ym0−ϵymq)am+ϵaq=b$$(y_{10}-ϵy_{1q})a_1+(y_{20}-ϵy_{2q})a_2+...+(y_{m0}-ϵy_{mq})a_m+ϵa_q=b$$

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

z=c1(y10−y1qϵ)+...+cm(ym0−ymqϵ)+cqϵ=z0+[cq−(c1y1q+...+cmymq)]ϵ$$\begin{gathered} z=c_1(y_{10}-y_{1q}ϵ)+...+c_m(y_{m0}-y_{mq}ϵ)+c_qϵ \\ =z_0+[c_q-(c_1{y}_{1q}+...+c_m{y}_{mq})]ϵ \end{gathered}$$

其中，z0=c1y10+...+cmym0$z_0=c_1{y}_{10}+...+c_m{y}_{m0}$，令zq=c1y1q+...+cmymq$z_q=c_1{y}_{1q}+...+c_m{y}_{mq}$，可以得到z=z0+(cq−zq)ϵ$z=z_0+(c_q-z_q)ϵ$，当z−z0=(cq−zq)ϵ<0$z-z_0=(c_q-z_q)ϵ<0$时，说明基本可行解对应的目标函数值变小了，即认为这个向量加入基向量组后是有益的。

那么可以对每一个q,m<q≤n$q,m<q\le n$计算cq−zq$c_q-z_q$，令ri=0(i=1,...,m),ri=ci−zi(i=m+1,...,n)$r_i=0(i=1,...,m),r_i=c_i-z_i(i=m+1,...,n)$，称ri$r_i$为第i$i$个简化价值系数，也称为检验数。可以通过证明得到，当且仅当相应的检验数都是非负的时候，该基本可行解是最优解。

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

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

确定停止条件

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

单纯形算法

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

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

To be continue…