优化理论之单纯形法

Nov 2, 2019 · 优化理论  ·
分享到:

由优化理论可知，如果线性规划中的最优解，那么这个最优解必是最优基本可行解。因此线性规划的核心就变成了求最优基本可行解。这构成了单纯形法的基础。

单纯形法核心思路：

  • 找到一个基本可行解
  • 求一个使目标函数值有所改善的基本可行解
  • 不断改进的基本可行解直到最优解

因此，大方向可以分为两个：1.找一个初时基本可行；2.改进基本可行解。对于第一个问题，用二阶段法或者大M法，第二个问题则需要考虑退化情形。

考虑线性规划问题 min⁡ f=defcTxs.t. Ax=b,x≥0(1.1) \begin{aligned} \min\ &f \overset{\mathop{def}}{=} \boldsymbol{c}^T\boldsymbol{x}\\ \mathop{s.t.}\ &\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b},\\ &\boldsymbol{x}\geq 0 \end{aligned}\tag{1.1} 其中，A\boldsymbol{A}是m×nm×n矩阵，秩为mm。cT\boldsymbol{c}^T是n维行向量，b≥0\boldsymbol{b}≥0是m维列向量(每个分量都大于等于0)。我们以列向量来表示A\boldsymbol{A} A=(p1,p2,⋯ ,pn) \boldsymbol{A}=(\boldsymbol{p_1},\boldsymbol{p_2},\dotsb,\boldsymbol{p_n}) 显然Ax\boldsymbol{A}\boldsymbol{x}可表示为： Ax=∑j=1nxjpj \boldsymbol{A}\boldsymbol{x}=\sum_{j=1}^n x_j\boldsymbol{p_j} xjx_j可以看成是列向量pj\boldsymbol{p_j}线性组合的系数。现将A\boldsymbol{A}分解成(B,N)(\boldsymbol{B},\boldsymbol{N})（可能经过列调换），使得其中B\boldsymbol{B}是基矩阵，N\boldsymbol{N}是非基矩阵。在此分解下，我们假设能得到一个基本可行解： x(0)=[B−1b0] \boldsymbol{x^{(0)}}=\begin{bmatrix} \boldsymbol{B^{-1}b}\\ \boldsymbol{0}\\ \end{bmatrix} 这个基本可行解对应了线性规划的某个极点，在点x(0)\boldsymbol{x^{(0)}}处的目标函数值是 f0=cTx(0)=(cBT,cNT)[B−1b0]=cBTB−1b(1.2) f_0=\boldsymbol{c}^T\boldsymbol{x^{(0)}}=(\boldsymbol{c_B^T,c_N^T})\begin{bmatrix} \boldsymbol{B^{-1}b}\\ \boldsymbol{0}\\ \end{bmatrix}\\ =\boldsymbol{c_B^T B^{-1}b}\tag{1.2} 其中，cBT\boldsymbol{c^T_B}是cT\boldsymbol{c}^T中与基变量对应的分量组成的m维行向量，cNT\boldsymbol{c^T_N}是cT\boldsymbol{c}^T中与非基变量对应的分量组成的n-m维行向量（因为为了让B\boldsymbol{B}的秩为m，可能对A\boldsymbol{A}中的列变量进行过变换，对应的x,c\boldsymbol{x,c}都要变换）。

现在我们从这个基本可行解x(0)\boldsymbol{x^{(0)}}出发，利用非基变量得出一个改进的基本可行解。设 x=[xBxN] \boldsymbol{x}=\begin{bmatrix} \boldsymbol{x_B}\\ \boldsymbol{x_N} \end{bmatrix} 是任一个可行解，则由Ax=b\boldsymbol{Ax=b}可得： xB=B−1b−B−1NxN(1.3) \boldsymbol{x_B}=\boldsymbol{B^{-1}b-B^{-1}Nx_N}\tag{1.3} 在此点处的目标函数值是 f=cTx=(cBT,cNT)[xBxN]=cBTxB+cNTxN=cBT(B−1b−B−1NxN)+cNTxN=cBTB−1b⏟基本可行解函数值−(cBTB−1N−cNT)xN=f0−∑j∈R(cBTB−1pj⏟zj−cj)xj, R为非基变量下标集=f0−∑j∈R(zj−cj)xj(1.4) \begin{aligned} f&=\boldsymbol{c^Tx}=(\boldsymbol{c^T_B,c^T_N})\begin{bmatrix}\boldsymbol{x_B}\\\boldsymbol{x_N}\end{bmatrix}\\ &=\boldsymbol{c^T_B x_B}+\boldsymbol{c^T_N x_N}\\ &=\boldsymbol{c^T_B(B^{-1}b-B^{-1}Nx_N)}+\boldsymbol{c^T_N x_N}\\ &=\underbrace{\boldsymbol{c^T_B B^{-1}b}}_{基本可行解函数值}-\boldsymbol{(c^T_B B^{-1}N-c^T_N)x_N}\\ &=f_0 - \sum_{j\in R}(\underbrace{\boldsymbol{c^T_B B^{-1}p_j}}_{z_j}-c_j)x_j,\ R为非基变量下标集\\ &=f_0 - \sum_{j ∈ R}(z_j-c_j)x_j \end{aligned}\tag{1.4} 由于cBT\boldsymbol{c^T_B}是m维行向量，B−1\boldsymbol{B^{-1}}是m×mm×m维矩阵，pj\boldsymbol{p_j}是m维列向量，因此它们的积是一个数字。由(1.4)(1.4)可知，适当的选择非基变量的值xjx_j，可以使得 ∑j∈R(zj−cj)xj>0(1.5) \sum_{j ∈ R}(z_j-c_j)x_j>0\tag{1.5} 例如，当zj−cj<0⇒xj=0z_j-c_j<0\Rightarrow x_j=0；当zj−cj≥0⇒xj>0z_j-c_j≥0\Rightarrow x_j>0。从而得到使目标函数的值减少的新的基本可行解。为了方便，我们不妨令n−m−1n-m-1个非基变量取0，一个非基变量，比如xkx_k大于0。需要注意的是，这个xkx_k的系数zj−cjz_j-c_j应该是大于0的。

根据(1.4)(1.4)，正系数zj−cjz_j-c_j越大，目标函数下降越快，我们不妨用贪心法，选择xkx_k，使 zk−ck=max⁡j∈R{zj−cj}>0(1.6) z_k-c_k=\max_{j\in R}\{z_j-c_j\}>0\tag{1.6} xkx_k由0变成正数后，得到方程组Ax=b\boldsymbol{Ax=b}的另一个解： xB=B−1b−B−1NxN=B−1b⏟bˉ−B−1pk⏟ykxk=bˉ−ykxk(1.7) \begin{aligned} \boldsymbol{x_B}&=\boldsymbol{B^{-1}b-B^{-1}Nx_N}\\ &=\boldsymbol{\underbrace{B^{-1}b}_{\bar b}-\underbrace{B^{-1}p_k}_{y_k}}x_k\\ &=\boldsymbol{\bar b}- \boldsymbol{y_k}x_k \end{aligned}\tag{1.7} 其中，bˉ,yk\boldsymbol{\bar b, y_k}都是m维列向量，xkx_k是标量。我们知道标准形式下的可行解要求xB≥0\boldsymbol{x_B}≥0，若我们把新得到的解xB\boldsymbol{x_B}按分量写出，即 xB=[xB1xB2⋮xBm]=[bˉ1bˉ2⋮bˉm]−[yk1yk2⋮ykm]xk≥0(1.8) \boldsymbol{x_B}=\begin{bmatrix}x_{B_1}\\x_{B_2}\\ \vdots \\x_{B_m}\end{bmatrix}=\begin{bmatrix}\bar{b}_1\\\bar{b}_2\\ \vdots \\\bar{b}_m\end{bmatrix}-\begin{bmatrix}y_{k_1}\\y_{k_2}\\ \vdots \\y_{k_m}\end{bmatrix}x_k ≥0 \tag{1.8} xN=(0,0,⋯ ,0,xk,0,⋯ ,0)T(1.9) \boldsymbol{x_N}=(0,0,\dotsb,0,x_k,0,\dotsb,0)^T\tag{1.9} 在新得到的点，目标函数的函数值为 f=f0−(zk−ck)xk(1.10) f=f_0-(z_k-c_k)x_k\tag{1.10} 再来，我们分析怎样确定xkx_k的取值。根据(1.8)(1.8)的可行性限制条件(非负)，我们不难得出，对∀i∈{1,2,⋯ ,m}且yki>0\forall i ∈ \{1,2,\dotsb,m\}且y_{k_i}>0.： xBi=bˉi−ykixk≥0⇒xk≤bˉiyki x_{B_i}=\bar b_i -y_{k_i}x_k≥0\Rightarrow x_k≤\frac{\bar b_i}{y_{k_i}} 因为若ykiy_{k_i}为负，xkx_k为正，这一分量不会成为破坏非负限制条件的分量。为了保证非负限制条件，我们必须有 xk≤min⁡{bˉiyki∣yki>0} x_k≤\min\big\{\frac{\bar b_i}{y_{k_i}}\big |y_{k_i}>0\big\} 同时，为了让ff尽可能减少，我们希望xkx_k尽量取大的值，所以有 xk=min⁡{bˉiyki∣yki>0}=bˉrykr(1.11) x_k=\min\big\{\frac{\bar b_i}{y_{k_i}}\big |y_{k_i}>0\big\}=\frac{\bar b_r}{y_{k_r}}\tag{1.11} 认为第rr项即为最小值。回顾一下刚刚的步骤，在目标函数中，因为xix_i是未定变量，因此我们先简单的使用贪心法，选择正系数最大的(zk−ck)(z_k-c_k)；接下来，再判断(zk−ck)(z_k-c_k)对应的变量xkx_k的取值范围，得到其最大值。

根据(1.8)(1.8)和(1.11)(1.11)，我们得到了一个改进可行解，其xr→0,xk→bˉrykrx_r\rightarrow 0, x_k\rightarrow\frac{\bar b_r}{y_{k_r}}，即 xB′=(xB1,xB2,⋯ ,xBr−1,0,xBr+1,⋯ ,xBm,0,⋯ ,xk,⋯ ,0)T \boldsymbol{x_B'}=(x_{B_1},x_{B_2},\dotsb,x_{B_{r-1}},0,x_{B_{r+1}},\dotsb,x_{B_m},0,\dotsb,x_k,\dotsb,0)^T

这个可行解一定是基本可行解。这是因为原来的基矩中， B=(pB1,pB2,⋯ ,pBm) \boldsymbol{B}=(\boldsymbol{p_{B_1},p_{B_2},\dotsb,p_{B_m}}) 中的m个列是线性无关的，其中不包含原属于N\boldsymbol{N}中的列pk\boldsymbol{p_{k}}。根据(1.7)(1.7)，有yk=B−1pk⇒pk=Byk=∑i=1mykipBi\boldsymbol{y_{k}}=\boldsymbol{B^{-1}p_k}\Rightarrow \boldsymbol{p_k}=\boldsymbol{By_k}=\sum_{i=1}^m y_{k_i}\boldsymbol{p_{B_i}}，即pk\boldsymbol{p_k}是线性无关向量组(pB1,pB2,⋯ ,pBm)(\boldsymbol{p_{B_1},p_{B_2},\dotsb,p_{B_m}})的线性组合，且至少存在一个非0系数ykry_{k_r}.根据替换定理，用pk\boldsymbol{p_{k}}来替代pBr\boldsymbol{p_{B_r}}得到如下向量组依然是线性无关的。 B=(pB1,pB2,⋯ ,pBr→pk,⋯ ,pBm) \boldsymbol{B}=(\boldsymbol{p_{B_1},p_{B_2},\dotsb,p_{B_r}\rightarrow p_{k},\dotsb,p_{B_m}}) 这样可以组成新的基矩阵B′\boldsymbol{B'}，并得到新的基解xB′\boldsymbol{x_B'}，而从(1.8)和(1.11)(1.8)和(1.11)，我们确定xB′\boldsymbol{x_B'}满足每一个分量非负的条件，因此其必然是一个基本可行解。

经过上述转换，xkx_k由原来的非基变量变成了基变量，而原来的基变量xBrx_{B_r}变成了非基变量，同时目标函数值比原来减少了(zk−ck)xk>0(z_k-c_k)x_k>0。重复以上过程，进一步改进基本可行解，直到式(1.6)(1.6)的结果只能为非正数，以致任何一个非基变量取正值都无法再使得目标函数值降低，此时即为最优解。

总结定理1：若在极小化问题中，对于某个基本可行解，所有(zj−cj)≤0(z_j-c_j)≤0，则这个这个基本可行解是最优解；

若在极大化问题中，对于某个基本可行解，所有(zj−cj)≥0(z_j-c_j)≥0，则这个这个基本可行解是最优解；

其中，zj−cj=cBTB−1pj−cj, j=1,⋯ ,nz_j-c_j=\boldsymbol{c_B^T B^{-1} p_j}-c_j,\ j=1,\dotsb,n

在线性规划中，通常称zj−cjz_j-c_j为判别数或检验数。

  1. 首先，我们要获得一个初始基本可行解（可通过直接计算，大M法，两阶段法获得），设初始基矩阵为B\boldsymbol{B}
  2. 在限制条件Ax=[B N][xBxN]=b\boldsymbol{Ax}=\bigl[\boldsymbol{B}\ \boldsymbol{N}\bigr] \bigl[ \begin{smallmatrix} \boldsymbol{x_B} \\ \boldsymbol{x_N} \end{smallmatrix} \bigr]=\boldsymbol{b}中，令非基变量xN=0\boldsymbol{x_N=0}，则BxB=b⇒xB=B−1b=bˉ\boldsymbol{Bx_B=b}\Rightarrow \boldsymbol{x_B}=\boldsymbol{B^{-1}b}=\boldsymbol{\bar b}，计算目标函数值f=cBTxBf=\boldsymbol{c^T_Bx_B}
  3. 对于此基矩阵的解集，求单纯形乘子wT=cBTB−1\boldsymbol{w^T=c_B^T B^{-1}}（m维行向量）。对于所有非基变量，计算判别数cBTB−1pj−cj=wTpj−cj=zj−cj\boldsymbol{c_B^T B^{-1}p_j}-c_j=\boldsymbol{w^T p_j}-c_j=z_j-c_j，求得zk−ck=max⁡j∈R{zj−cj}z_k-c_k=\max\limits_{j\in R}\{z_j-c_j\}。若zk−ck≤0z_k-c_k≤0，停止计算，目前基本可行解是最优解。否则，进行步骤（3）
  4. xB′=B−1b⏟bˉ−B−1pk⏟ykxk\boldsymbol{x_B'}=\boldsymbol{\underbrace{B^{-1}b}_{\bar b}-\underbrace{B^{-1}p_k}_{y_k}}x_k。若yk≤0\boldsymbol{y_k}≤0，即每一个分量都不大于0，则停止计算，问题不存在有限最优解；否则，进行步骤（4）
  5. 确定下标rr，使得xk=bˉrykr=min⁡{bˉiyki∣yki>0}x_k=\frac{\bar b_r}{y_{k_r}}=\min\big\{\frac{\bar b_i}{y_{k_i}}\big |y_{k_i}>0\big\}。xBrx_{B_r}为离基变量，xkx_k为进基变量，用pk\boldsymbol{p_k}替代变量pBrp_{B_r}得到新的基矩阵B′→B\boldsymbol{B'\rightarrow B}，返回步骤（1）

对于非退化情形，单纯形法有优秀的结论：

定理2：对于非退化问题，单纯形方法经过有限次迭代后，能达到最优解或得出无界的结论。在此条件下，单纯形法是收敛的。

证明：

单纯形法收敛性非退化

Figure 1: 单纯形法收敛性非退化

单纯形法收敛性非退化

对于非退化情形，每次迭代都有 xB=B−1b=bˉ>0 \boldsymbol{x_B=B^{-1}b=\bar b}>0 此时能保证xr=brykr>0x_r=\frac{b_r}{y_{k_r}}>0。因此经迭代，目标函数值减小，并且由此可知，各次迭代得到的基本可行解互不相同。由于基本可行解的个数有限，因此经有限次迭代必能达到最优解。对于退化情形，我们在后面将要证明，如果最优解存在，只要采取一定的措施，也能做到有限步收敛。

用单纯形方法求解线性规划问题的过程，实际上就是解线性方程组。只是在每次迭代中，要按一定规则选择自由未知量，以便能够得到改进的基本可行解。这个求解过程可以通过变换单纯形表来实现。具体做法可见《最优化理论与算法（第2版）》第3.1.4节。

对于退化情形，即存在基变量为0的情形（基变量为0是否一定退化有待验证），经过单纯形法多次迭代后，可能出现循环解现象（退化常见而循环不常见）。因此，我们可以采用摄动法、Bland法，字典序法。