【线性规划(三)】单纯型法（上）

1 最直观的一种方法：穷举顶点法

首先来复习一下线性规划基本定理:

对于标准形式的线性规划问题，如果该问题存在有界的最优解，那么至少有一个最优解在顶点上。

根据线性规划基本定理，我们不难想到一个求解线性规划最优解的方法就是穷举所有的顶点，然后找出目标函数最优的那个顶点就是最优解了。 例如考虑如下的标准形式的线性规划问题：

\begin{aligned} \text{minimize } x_1+2x_2 \;\;\;\;(1.1) \\ x_1+x_2+x_3 = 3 \;\;\;\;(1.2) \\ -x_1 + x_2 +x_4 = 2 \;\;\;\;(1.3) \\ x_1,x_2,x_3,x_4 \ge0\;\;\;\;(1.4) \\ \end{aligned} \\

画出该线性规划问题的可行域，并在图中标出所有顶点的坐标：

这里我们简单回顾一下在上一节笔记中我们求所有顶点的过程：

Step 1: 令 x_3=0, x_4=0 代入(1.2-1.4)，构成方程组 x_1+x_2=3, -x_1+x_2=2 ，解上述方程组即可得到顶点 v_1 ;

Step 2: 令 x_2=0, x_3=0 代入(1.2-1.4)，构成方程组 x_1=3, -x_1+x_4=2 ，解上述方程组即可得到顶点 v_2 ;

Step 3: 令 x_1=0, x_4=0 代入(1.2-1.4)，构成方程组 x_2+x_3=3, -x_2=2 ，解上述方程组即可得到顶点 v_3 ;

Step 4: 令 x_1=0, x_2=0 代入(1.2-1.4)，构成方程组 x_3=3, x_4=2 ，解上述方程组即可得到顶点 v_4 ;

Step 5: 令 x_1=0, x_3=0 代入(1.2-1.4)，构成方程组 x_2=3, x_2+x_3=2 ，解上述方程组即可得到顶点 v_5 ;

Step 6: 令 x_2=0, x_4=0 代入(1.2-1.4)，构成方程组 x_1+x_3=3, -x_1=2 ，解上述方程组即可得到顶点 v_6 ;

v_1=\left( \begin{array}{c} 0.5\\ 2.5\\ 0\\ 0\\ \end{array} \right) ,v_2=\left( \begin{array}{c} 3\\ 0\\ 0\\ 5\\ \end{array} \right) ,v_3=\left( \begin{array}{c} 0\\ 2\\ 1\\ 0\\ \end{array} \right) ,v_4=\left( \begin{array}{c} 0\\ 0\\ 3\\ 2\\ \end{array} \right) ,v_5=\left( \begin{array}{c} 0\\ 3\\ 0\\ -1\\ \end{array} \right) ,v_6=\left( \begin{array}{c} -2\\ 0\\ 5\\ 0\\ \end{array} \right) \\

v_5 和 v_6 顶点包含负数的分量，因此 v_5 和 v_6 不是可行的顶点，故此舍弃掉 v_5 和 v_6 顶点，只保留顶点 v_1, v_2, v_3 和 v_4

接下来我们将上述求解过程推广到一般的标准形式的线性规划问题，考虑如下的标准形式线性规划问题。

\begin{aligned} \underset{x\in R^n}{\min} \; c^Tx \;\;\;\;\;\;\;\;\;\;\;\;\; (1.5) \\ \text{s.t.} \ Ax=b,\ x\ge 0 \;\;\;\; (1.6) \end{aligned} \\

其中 A\in R^{m \times n} ，并且 n \ge m 。

从这里我们不难总结出对于一般的标准形式的线性规划问题的穷举顶点算法为：

• Step 1: 从 n 个变量中任意抽取其中 m 个，然后将剩余的 n-m 个变量赋值为0。
• Step 2: 抽取得到的 m 个变量构成 m 乘 m 的方程组，求解这个方程组即可获得顶点。
• Step 3: 判断这个顶点是否满足约束 x \ge 0 ，若不满足则舍弃掉，若满足则保存该顶点。

那么从上述算法流程不难看出：

1. 上述算法要循环 C^{n}_{m} 次必然可以求解到最优解，也就是说需要求解 C^{n}_{m} 次 规模大小为 m 个未知数的线性代数方程组。
   
2. 上述算法所需的计算量是非常恐怖的，尤其是在 m 比较大的情况下。
   

到这里我们会进一步思考：穷举所有顶点是不是显得比较笨一些，那有没有一种更加效率的方法来求解线性规划问题呢？

答案就是 我们线性规划最重要的算法：单纯型法。

2 从穷举顶点法过渡到单纯型法

下面我们直接给出单纯型算法的大致框架：

• Step 1：从一个初始顶点出发；
• Step 2：检查是否是最优解（最优解条件），若是最优解则停止，否则进入下一步；
• Step 3：从当前顶点移动到更好的顶点，然后跳转到 Step 2；

仔细分析一下单纯型算法的每一步。Step1 是从一个初始顶点出发，这是很自然的想法，因为根据线性规划基本定理 若线性规划存在最优解，顶点中肯定至少有一个是最优解。Step 2 是判断算法是否终止，如果找到最优解的就停止，如果没有就接下去再找。具体怎么判断当前解是不是最优解接下来我们会详细介绍。Step 3 进入到这一步说明当前的解还不是最优的解，那我们自然会想办法让改变当前的解，让它往更接近最优解的方向上去靠拢。

对比穷举算法和单纯形算法可以发现： 在第1节中的穷举算法中，我们是把所以的顶点都拿出来比较一番，然后就可以找出最优解了。单纯型法和穷举算法的主要区别在于 单纯型法是一个迭代的方法。单纯型法是通过从一个可能不是很优的可行解出发，然后逐步逐步改进这个可行解，直至达到最优解。

上述算法的缺陷：顶点可能压根就不是一个可行解。根据基可行解的定义，进一步改进上述算法可得：

• Step 1：从一个初始的基可行解出发；
• Step 2：检查是否是最优解（最优解条件），若是最优解则停止，否则进入下一步；
• Step 3：从当前基可行解移动到更好的基可行解，然后跳转到 Step 2；

上述算法即为单纯型算法的最基本的步骤。显然，在上述算法过程中我们需要解决三个主要的问题：

1. 如何找到一个初始的基可行解；
2. 如何检查当前解是否是最优解；
3. 如何从当前的基可行解移动到更好的基可行解；

首先问题1和问题2暂时先放一放留在后边再解答，我们先来考虑问题3。为了解决问题3，我们先定义相邻顶点的概念：

我们还是把上述线性规划构成的可行域画出来，如下图所示：

在上图中我们标出了可行域的所有的顶点为

v_1=\left( \begin{array}{c} 0.5\\ 2.5\\ 0\\ 0\\ \end{array} \right) ,v_2=\left( \begin{array}{c} 3\\ 0\\ 0\\ 5\\ \end{array} \right) ,v_3=\left( \begin{array}{c} 0\\ 2\\ 1\\ 0\\ \end{array} \right) ,v_4=\left( \begin{array}{c} 0\\ 0\\ 3\\ 2\\ \end{array} \right) \\

从图中不难看出 v_1 和 v_2 是相邻节点, v_2 和 v_4 是相邻节点, v_3 和 v_4 是相邻节点。 v_1 的基变量为 x_1,x_2 非基变量为 x_3,x_4 ， v_2 的基变量为 x_1,x_4 非基变量为 x_2,x_3 ，观察可知将 v_1 的非基变量 x_4 转为 基变量。由于 x_4 从等于0，变化成大于等于0，为了保持约束的平衡还需要 v_1 的基变量 x_2 转为 非基变量。这样就完成了从 v_1 变化为 v_2

总结一下一般的情况，顶点 v_i 和顶点 v_j 是相邻节点的定义为 v_i 和 v_j 的非基变量只有一个不一样，基变量也只有一个不一样，其余变量均相等。

进一步将这一过程用代数方式严谨表述出来，其中 B \in R^{m \times m} , N \in R^{m \times(n-m)}

\begin{aligned} x=\left[ \begin{array}{c} x_B\\ x_N\\ \end{array} \right] \end{aligned} \\

其中 x_B \in R^m 是基变量， x_N \in R^{(n-m)} 。

由此可以将线性规划(1.5-1.6)改写为如下形式：

\begin{aligned} \underset{x\in R^n}{\min} \;c_B^Tx_B+c_N^Tx_N \;\;\;\;\;\;\;\;\;\;\;\;\; (2.1) \\ \text{s.t.} \ Bx_B + Nx_N=b,\ x_B \ge 0, x_N \ge 0 \;\;\;\; (2.2) \end{aligned} \\

由式(2.2)可得：

\begin{aligned} x_B=B^{-1}b-B^{-1}Nx_N\;\;\;\; (2.3) \end{aligned} \\

在非基变量中我们选择其中一个分量 x_q 进入到 基变量中

\begin{aligned} x_B=B^{-1}b-x_qB^{-1}A_q\;\;\;\; (2.4) \end{aligned} \\

A_q 是与 x_q 所对应的 A 的列。式(2.4)反映了基变量和非基变量之前的关系。

设当前的解为 x ，我们想要得到下一步的解为 x(\lambda)

\begin{array}{c} x\left( \lambda \right) :=x+\lambda d_q\;\;\;\;\left( 2.5 \right)\\ d_q=\left[ \begin{array}{c} -B^{-1}A_q\\ e_q\\ \end{array} \right] \;\;\;\;\left( 2.6 \right)\\ \end{array}

其中 d_q 为搜索方向， \lambda \ge 0 为搜索步长。 d_q \in R^n ， e_q \in R^{n-m} 其在 x_q 对应的位置为 1，其余位置为 0。

验证 d_q 是一个可行方向吗？

Ad_q=\left[ \left. B \right|N \right] \left( \frac{-B^{-1}A_q}{e_q} \right) =A_q-N_q=0 \;\;\;\; (2.7) \\

由此说明 d_q 在 A 矩阵的零空间里。

不光要找到一个相邻的节点，还要找到一个充分好的相邻节点： 将式(2.5)带入到目标函数中有：

\begin{aligned} c^Tx\left( \lambda \right) &=c^T\left( x+\lambda d_q \right) =c^Tx+\lambda \left[ \left. c_{B}^{T} \right|c_{N}^{T} \right] \left( \frac{-B^{-1}A_q}{e_q} \right) \\ &=c^Tx+\lambda \left( c_q-c_{B}^{T}B^{-1}A_q \right) \\ & =c^Tx+\lambda r_q \end{aligned} \\

如果 r_q < 0 ，那么这就是一个好的移动方向。其中

\begin{aligned} r_q=c_q-c_{B}^{T}B^{-1}A_q \;\;\;\; (2.8) \end{aligned} \\

r_q 表示的含义为非基变量 q 加入进来可以让目标函数下降的值，我们将 r_q 称之为 reduced cost

有可能有多个非基变量都可以使得 r_q<0 ，那么选择哪个非基变量进基呢？

1. 如果选择一个非基变量进基，显然有 r_q = c_q - c^T_BB^{-1}A_q = 0
2. 选择一个好的方向 d_q 使得 r_q < 0

目前比较常用的有两种选择非基变量进基的方法：

1. 最小索引法：采用 索引下标 q 最小的 r_q<0 的非基变量；
2. 最速下降法：采用目标函数最小的优化问题来选择出非基变量：

\begin{aligned} \underset{j\in N}{\min}\frac{c^Td_j}{\lVert d_j \rVert} \;\;\;\; (2.9) \end{aligned} \\

上述目标函数中 c^Td_j 表示让下一步搜索的目标函数最小，除以 \lVert d_j \rVert 代表着要做归一化。看起来最速下降法选择出来的非基变量可能会更好一些。事实上，目前没有明显的证据表明哪种选择非基变量的方法就一定好。

通过上面的方法确定了搜索方向 d_q 之后，接下来就是步长 \lambda \ge 0 的选择。

根据定理3可知，若 d_q \ge 0 是一个可行方向，那么表示此线性规划问题是没有下界的。若 d_q 存在某些分量是小于0的，之前我们验证过 Ad_q = 0 ，为了保证下一步到达的位置依然是可行解，只需要保证 x+\lambda d_q \ge 0 即可，根据这一原则有如下步长表达式：

\begin{aligned} \lambda =\underset{i\in B}{\min}\left\{ \left. \frac{x_i}{-d_{i}^{q}} \right|d_{i}^{q}<0 \right\} \;\;\;\; (2.10) \end{aligned} \\

其中 x_i 为当前搜索点的第 i 个分量， d^q_i 为 d_q 的第 i 个分量。

上式的含义就是在 d^q_i<0 的那些分量中，选择让 \frac{x_i}{-d_{i}^{q}} 最小的作为步长。

3 单纯型算法流程

• Step 1: 找到基可行解 x , 设其基矩阵为 B 和非基矩阵 N ;
• Step 2: 根据式(2.8)计算所有非基矩阵的 r_q ，若所有的 r_q \ge 0 则停止输出最优解；否则 跳转到 Step 3;
• Step 3: 根据式(2.5)计算移动方向 d_q ，若 d_q \ge 0 , 则表明线性规划是无界的。 否则，根据式(2.10)计算出步长 \lambda 。更新 x\gets x+\lambda d_q ，然后更新基矩阵，并跳转到 Step 2;

下面通过一个具体的实例来展示单纯型法的计算流程

\begin{aligned} \min -x_1-x_2-x_3 \\ \text{s.t.} \;\; 2x_1 + x_4 =1 \\ 2x_2 + x_5 = 1 \\ 2x_3 + x_6 =1 \\ x_1,x_2,x_3,x_4,x_5,x_6 \ge 0 \end{aligned} \\

由此可知约束矩阵

A=\left( \begin{matrix} 2& 0& 0& 1& 0& 0\\ 0& 2& 0& 0& 1& 0\\ 0& 0& 2& 0& 0& 1\\ \end{matrix} \right) \\ b=\left( \begin{matrix} 1& 1& 1\\ \end{matrix} \right) ^T \\ c=\left( \begin{matrix} -1& -1& -1& 0& 0& 0\\ \end{matrix} \right) ^T \\

Step 1: 取 x_4,x_5,x_6 作为基变量，有

\left[ \frac{x_B}{x_N} \right] =\left( \left. x_4\ x_5\ x_6 \right|x_1\ x_2\ x_3 \right) ^T=\left( \left. 1\ 1\ 1 \right|0\ 0\ 0 \right) ^T \\

可以得到基矩阵和非基矩阵：

B=\left[ \begin{matrix} 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ \end{matrix} \right] ,\ \ N=\left[ \begin{matrix} 2& 0& 1\\ 0& 2& 0\\ 0& 0& 2\\ \end{matrix} \right] \\ \tilde{B}=\left\{ 4,5,6 \right\} ,\ \ \tilde{N}=\left\{ 1,2,3 \right\} \\ c_B=\left( 0\ 0\ 0 \right) ^T,\ \ c_N=\left( -1\ -1\ -1 \right) ^T \\

Step 2: 根据式(2.8)计算所有非基矩阵的 reduced cost:

r_4=c_4-c_{B}^{T}B^{-1}A_4=-1-\left( 0\ 0\ 0 \right) \left[ \begin{matrix} 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ \end{matrix} \right] ^T\left( \begin{array}{c} 2\\ 0\\ 0\\ \end{array} \right) =-1 \\

r_5, r_6 的计算方式和 r_4 类似就不在这里展开写了，最终整理可得：

r=\left( 0\ 0\ 0\ -1\ -1\ -1 \right) ^T \\

由于 r 不满足所有分量都大于等于0，因此跳转到 Step 3

Step 3：

采用最小索引法来选择非基变量进基，因此非基变量 x_1 被作为进基变量，由此知 q=1 。根据式(2.6)可以算出迭代方向：

d_1=\left( \frac{-B^{-1}A_1}{e_1} \right) =\left( -2\ 0\ 0\ 1\ 0\ 0 \right) ^T \\

进一步根据式(2.10)计算出步长：

\begin{aligned} \lambda =\underset{i\in B}{\min}\left\{ \left. \frac{x_i}{-d_{i}^{q}} \right|d_{i}^{q}<0 \right\} =\min \left\{ \frac{1}{2} \right\} =0.5 \end{aligned} \\

更新变量：

\begin{aligned} x&:=x+\lambda d_q=\left( 1\ 1\ 1\ 0\ 0\ 0 \right) ^T+0.5\left( -2\ 0\ 0\ 1\ 0\ 0 \right) ^T\\ &=\left( 0\ 1\ 1\ 0.5\ 0\ 0 \right) ^T \end{aligned} \\

基变量变为：

x_B=\left( x_5,x_6,x_1 \right) ^T=\left( 1\ 1\ 0.5 \right) ^T \\

非基变量为：

x_N=\left( x_2,x_3,x_4 \right) ^T=\left( 0\ 0\ 0 \right) ^T \\

4 附录（理论证明）

为了保证前面算法过程描述的完整性，我们将单纯型算法的理论保证部分单独抽出来放到这一节来集中讲述。由于这些定理的证明过程可能会涉及到一些线性规划课程以外的知识，所以我们对于证明过程采取一些简化，重点会放在这些定理的意义上。那么整个单纯型算法有以下五个比较重要的定理。

定理 1 ：若满足

x^*=\left[ \frac{x_{B}^{*}}{x_{N}^{*}} \right] =\left[ \frac{B^{-1}b}{0} \right] \ge 0 \\

是基可行解，并且不存在 r_q < 0 , 则表明 x^* 是线性规划的最优解。

该定理的证明在这里不予给出了。想要证明该定理需要用到凸优化的知识，大致思路是线性规划就是一个凸优化问题，若 不存在 r_q < 0 则表明是当前顶点是一个局部最优点，同时在凸优化问题中局部最优就等价于全局最优，由此定理1可以得到证明。

定理 2 ： 若

x =\left[ \frac{B^{-1}b}{0} \right] \\

是基可行解， B 是基矩阵。若对于某些非基变量 x_q ，使得 r_q < 0 , 那么式(2.5)中的 d_q 可以使得目标函数得到改进。

定理2 在前面推导出方向 d_q 的过程中已经被证明过了，这里就不再赘述了。

定理 3 ： x 是线性规划的基可行解。对于一些非基变量 x_q ，若存在可行方向 d_q \ge 0 , 使得 reduced cost r_q < 0 , 那么线性规划是无下界的。

定理3 是用作判定线性规划是否存在有界的最优解的定理。显然若 若存在可行方向 d_q \ge 0 , 使得 reduced cost r_q < 0 ，那么我们只需要让步长取无穷大就可以一直让目标函数保持下降，这就说明了最优解会出现在无穷远处，定理3可以被证明。

定理 4 若

x =\left[ \frac{B^{-1}b}{0} \right] \\

是基可行解， B 是基矩阵。若对于某些非基变量 x_q ，使得 reduced cost r_q < 0 , 那么采用式(2.5)计算 d_q ，采用式(2.10)计算步长 \lambda , 可以使得目标函数得到改进。

定理4 在前面推导出方向 d_q 的过程中已经被证明过了，这里就不再赘述了。

定理 5 : 在非退化情况下，线性规划单纯性算法可以保证在有限步终止。

关于退化的情况我们会在后边介绍，这里就把非退化的情况理解为“正常情况”即可。退化的情况有可能让单纯型法重复/循环遍历某些顶点，这样的情况是非常不好的一种情况。若没有退化情况发生，则由于顶点的数量是有限的，单纯型法最多会在穷举所有顶点后停止，因此定理5是成立的。

上面5条定理看似非常多非常混乱，其实每一条都是非常有存在的必要的，最后我们再简单梳理一下：

1. 定理1 是判断当前顶点是不是最优解
2. 定理2和定理4 是说明合理地选择迭代方向和迭代步长可以保证让每次迭代后目标函数有所改进
3. 定理3 是判断当前线性规划是不是无界
4. 定理5 是对单纯型法收敛性的一个基本保障

有了如上的梳理相信你对这5条定理会有一个基本的认识，哪怕你对这些定理的推导过程丝毫不了解，这也不是很要紧。只要你知道每条定理的出发点，那么如果需要用到这个定理的证明去查对应的教科书即可。

参考文献

【1】Padberg M. Linear optimization and extensions[M]. Springer Science & Business Media, 2013.