优化理论之单纯形法

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

单纯行法

核心思想

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

单纯形法核心思路：

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

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

单纯形方法原理

考虑线性规划问题 \[ \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} \] 其中，\(\boldsymbol{A}\)是\(m×n\)矩阵，秩为\(m\)。\(\boldsymbol{c}^T\)是n维行向量，\(\boldsymbol{b}≥0\)是m维列向量(每个分量都大于等于0)。我们以列向量来表示\(\boldsymbol{A}\) \[ \boldsymbol{A}=(\boldsymbol{p_1},\boldsymbol{p_2},\dotsb,\boldsymbol{p_n}) \] 显然\(\boldsymbol{A}\boldsymbol{x}\)可表示为： \[ \boldsymbol{A}\boldsymbol{x}=\sum_{j=1}^n x_j\boldsymbol{p_j} \] \(x_j\)可以看成是列向量\(\boldsymbol{p_j}\)线性组合的系数。现将\(\boldsymbol{A}\)分解成\((\boldsymbol{B},\boldsymbol{N})\)（可能经过列调换），使得其中\(\boldsymbol{B}\)是基矩阵，\(\boldsymbol{N}\)是非基矩阵。在此分解下，我们假设能得到一个基本可行解： \[ \boldsymbol{x^{(0)}}=\begin{bmatrix} \boldsymbol{B^{-1}b}\\ \boldsymbol{0}\\ \end{bmatrix} \] 这个基本可行解对应了线性规划的某个极点，在点\(\boldsymbol{x^{(0)}}\)处的目标函数值是 \[ 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} \] 其中，\(\boldsymbol{c^T_B}\)是\(\boldsymbol{c}^T\)中与基变量对应的分量组成的m维行向量，\(\boldsymbol{c^T_N}\)是\(\boldsymbol{c}^T\)中与非基变量对应的分量组成的n-m维行向量（因为为了让\(\boldsymbol{B}\)的秩为m，可能对\(\boldsymbol{A}\)中的列变量进行过变换，对应的\(\boldsymbol{x,c}\)都要变换）。

现在我们从这个基本可行解\(\boldsymbol{x^{(0)}}\)出发，利用非基变量得出一个改进的基本可行解。设 \[ \boldsymbol{x}=\begin{bmatrix} \boldsymbol{x_B}\\ \boldsymbol{x_N} \end{bmatrix} \] 是任一个可行解，则由\(\boldsymbol{Ax=b}\)可得： \[ \boldsymbol{x_B}=\boldsymbol{B^{-1}b-B^{-1}Nx_N}\tag{1.3} \] 在此点处的目标函数值是 \[ \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} \] 由于\(\boldsymbol{c^T_B}\)是m维行向量，\(\boldsymbol{B^{-1}}\)是\(m×m\)维矩阵，\(\boldsymbol{p_j}\)是m维列向量，因此它们的积是一个数字。由\((1.4)\)可知，适当的选择非基变量的值\(x_j\)，可以使得 \[ \sum_{j ∈ R}(z_j-c_j)x_j>0\tag{1.5} \] 例如，当\(z_j-c_j<0\Rightarrow x_j=0\)；当\(z_j-c_j≥0\Rightarrow x_j>0\)。从而得到使目标函数的值减少的新的基本可行解。为了方便，我们不妨令\(n-m-1\)个非基变量取0，一个非基变量，比如\(x_k\)大于0。需要注意的是，这个\(x_k\)的系数\(z_j-c_j\)应该是大于0的。

根据\((1.4)\)，正系数\(z_j-c_j\)越大，目标函数下降越快，我们不妨用贪心法，选择\(x_k\)，使 \[ z_k-c_k=\max_{j\in R}\{z_j-c_j\}>0\tag{1.6} \] \(x_k\)由0变成正数后，得到方程组\(\boldsymbol{Ax=b}\)的另一个解： \[ \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} \] 其中，\(\boldsymbol{\bar b, y_k}\)都是m维列向量，\(x_k\)是标量。我们知道标准形式下的可行解要求\(\boldsymbol{x_B}≥0\)，若我们把新得到的解\(\boldsymbol{x_B}\)按分量写出，即 \[ \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} \] \[ \boldsymbol{x_N}=(0,0,\dotsb,0,x_k,0,\dotsb,0)^T\tag{1.9} \] 在新得到的点，目标函数的函数值为 \[ f=f_0-(z_k-c_k)x_k\tag{1.10} \] 再来，我们分析怎样确定\(x_k\)的取值。根据\((1.8)\)的可行性限制条件(非负)，我们不难得出，对\(\forall i ∈ \{1,2,\dotsb,m\}且y_{k_i}>0\).： \[ x_{B_i}=\bar b_i -y_{k_i}x_k≥0\Rightarrow x_k≤\frac{\bar b_i}{y_{k_i}} \] 因为若\(y_{k_i}\)为负，\(x_k\)为正，这一分量不会成为破坏非负限制条件的分量。为了保证非负限制条件，我们必须有 \[ x_k≤\min\big\{\frac{\bar b_i}{y_{k_i}}\big |y_{k_i}>0\big\} \] 同时，为了让\(f\)尽可能减少，我们希望\(x_k\)尽量取大的值，所以有 \[ 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} \] 认为第\(r\)项即为最小值。回顾一下刚刚的步骤，在目标函数中，因为\(x_i\)是未定变量，因此我们先简单的使用贪心法，选择正系数最大的\((z_k-c_k)\)；接下来，再判断\((z_k-c_k)\)对应的变量\(x_k\)的取值范围，得到其最大值。

根据\((1.8)\)和\((1.11)\)，我们得到了一个改进可行解，其\(x_r\rightarrow 0, x_k\rightarrow\frac{\bar b_r}{y_{k_r}}\)，即 \[ \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 \]

这个可行解一定是基本可行解。这是因为原来的基矩中， \[ \boldsymbol{B}=(\boldsymbol{p_{B_1},p_{B_2},\dotsb,p_{B_m}}) \] 中的m个列是线性无关的，其中不包含原属于\(\boldsymbol{N}\)中的列\(\boldsymbol{p_{k}}\)。根据\((1.7)\)，有\(\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}}\)，即\(\boldsymbol{p_k}\)是线性无关向量组\((\boldsymbol{p_{B_1},p_{B_2},\dotsb,p_{B_m}})\)的线性组合，且至少存在一个非0系数\(y_{k_r}\).根据替换定理，用\(\boldsymbol{p_{k}}\)来替代\(\boldsymbol{p_{B_r}}\)得到如下向量组依然是线性无关的。 \[ \boldsymbol{B}=(\boldsymbol{p_{B_1},p_{B_2},\dotsb,p_{B_r}\rightarrow p_{k},\dotsb,p_{B_m}}) \] 这样可以组成新的基矩阵\(\boldsymbol{B'}\)，并得到新的基解\(\boldsymbol{x_B'}\)，而从\((1.8)和(1.11)\)，我们确定\(\boldsymbol{x_B'}\)满足每一个分量非负的条件，因此其必然是一个基本可行解。

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

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

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

其中，\(z_j-c_j=\boldsymbol{c_B^T B^{-1} p_j}-c_j,\ j=1,\dotsb,n\)

在线性规划中，通常称\(z_j-c_j\)为判别数或检验数。

一般计算步骤

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

收敛性

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

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

证明：

单纯形法收敛性非退化

单纯形法收敛性非退化

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

使用表格形式的单纯形方法

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

退化情形

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