线性规划(一)——单纯形法

2018-07-14  本文已影响0人  心露边白

问题介绍

单纯形法(simplex method)是求解线性规划问题一种通用算法,在实际生产生活中有广泛的应用。有些教材在介绍单纯形法时使用了复杂的矩阵和下标运算,使得算法的脉络难以把握。本文试图从几何意义和简明的代数运算出发,解释单纯形法的运行逻辑。

标准形

标准形更加常用的是以下矩阵形式:

线性约束条件 Ax=bm 个线性方程组成,为了理论上的方便,一般要求这 m 个方程线性无关,即 \mathrm{rank}~A=m . 以下均假设该条件成立。


下面来讨论标准形的几何意义。众所周知,下面的方程
a_1x_1+\cdots+a_nx_n=b 表示 n 维欧氏空间 \mathbb{R}^n 中的一个 n-1超平面m 个线性无关方程构成的方程组 Ax=b 则表示 m 个超平面的交,也就是一个 n-m 维超平面。举例来说,若 (n,m)=(3,1), 则线性约束条件转化为 \mathbb{R}^3 中的一个平面(如下图所示)。

(n,m)=(3,1) 的线性规划问题

x\geqslant0 的非负约束条件下,n-m 维超平面 S 被坐标轴截为 n-m 维凸图形 S_0, 从而线性函数 z=c^\mathrm{T}x 的最小值一定在 S_0 的顶点上取。确切地说,这些顶点是 S 与某 m 条坐标轴张成的 m 维子空间的交点。在 (n,m)=(3,1) 的条件下,平面 S 被坐标轴截为三角形 S_0, 并且 S_0 的顶点恰好是 S 与三个一维坐标轴的交点。

S_0 的顶点由 Sm 维坐标轴子空间相交而成,而 m 条坐标轴又选自于 \mathbb{R}^nn 条坐标轴,因此 S_0 的顶点数不超过 C_n^m. (实际可能会少的多)


回到矩阵形式。记原始的线性规划问题为 \mathcal{P},则超平面 S 表示 Ax=b 的解,超平面凸图形 S_0 表示 \mathcal{P}可行域。由于线性函数 c^\mathrm{T}x 的最小值在 S_0 的顶点取到,因此 S_0 的顶点处的性质具有“重要的研究意义”。若 PS_0 的一个顶点,则 Pn 个分量中只有坐标轴对应的 m 个分量可能取非零值,其余 n-m 个分量均必取 0

对于“坐标轴”“顶点”的初步印象自然引出了“基向量组”“基本解”的概念:

基向量组、基变量、基本解

上述对于基本解的定义是合理的。基向量组 B=(P_{i_1},\cdots,P_{i_m}) 也可以理解为列向量的水平叠置,因而 B 也被称为基向量矩阵。在下面的叙述当中,基向量组和基向量矩阵的概念是等同的,不再加以区分。B 所对应的基本解仅可能在 i_1,\cdots,i_mm 个分量上取非零值,将这 m 个分量构成的向量记为 x_B~(x_B\in\mathbb{R}^m),则原 n 元线性约束条件可简化为 m 元的情形:
Ax=b\Rightarrow Bx_B=b由于 Bm 阶可逆矩阵,故 x_B=B^{-1}b。这说明基向量组 B 对应的基本解的确是唯一的。

更进一步,如果基向量组 B 对应的基本解 x 满足 x\geqslant0,则 x 位于可行域当中,因而 B 称为可行基向量组x 称为基本可行解。从几何意义来看,基本可行解就是可行域 S_0 的顶点。如果 \mathcal{P}最优解(使目标函数取最小值的解)存在,则一定存在一个基本可行解是最优的,换言之,最优解(目标函数最小值)一定能在基本可行解(S_0 的顶点)上取到

如上图所示,P_1,P_2,P_3 三个顶点都分别只有一个分量非零,因而它们都是基本解。由于它们都在各自坐标轴的正半轴上,因此都是基本可行解。如果移动 S 使得 Sx_1 轴的交点 P_1 在负半轴上,那么 P_1 就不是基本可行解了。


以上的讨论都是关于线性规划问题 \mathcal{P}静态性质。下面所要讨论的单纯形法,则是真正意义上给出了求出最优解的算法。尽管 \mathcal{P} 的几何意义对于单纯形法而言没有实质性的作用,但它提供了理解诸如“基本解”等概念的一种有效途径。

单纯形法的核心步骤是:(a) 最优性检验; (b) 基变换。简单来说,这两个步骤的主要作用是:

交替进行上述2个步骤,直至找到最优解或证明最优解不存在。

最优性检验

给定可行基向量组 B=(P_{\pi(1)},\cdots,P_{\pi(m)}),为了讨论基本解的最优性,需要计算目标函数 c^{\mathrm{T}}x 在基本解附近的值。

B\in\mathbb{R}^{m\times m}A\pi 分量部分构成的基向量组,而 N\in\mathbb{R}^{m\times(n-m)}A 的非 \pi 分量部分构成的非基向量组。设 x\in S_0 是可行域内任一点,x\pi 分量部分记为 x_B\in\mathbb{R}^m,非 \pi 分量部分记为 x_N\in\mathbb{R}^{n-m}c_B,c_N 的意义类似。于是
Ax=b\Rightarrow (B~N) \begin{pmatrix} x_B\\ x_N \end{pmatrix}=b \Rightarrow x_B=B^{-1}(b-Nx_N) 上式表明可行解的基分量部分 x_B 可由非基分量部分 x_N 唯一线性表示。这一结论从几何意义上是容易理解的:Sn-m 维超平面,含有 n-m 个自由变量,因此 S 上的任一点的坐标可以仅由其中的 n-m 个线性表示。

于是,z=c^{\mathrm{T}}x 也应当仅由 n-m 个分量表示:
\begin{aligned} z&=c^{\mathrm{T}}x\\ &=(c_B^{\mathrm{T}}~c_N^\mathrm{T})\begin{pmatrix}x_B\\ x_N\end{pmatrix}\\ &=c_B^{\mathrm{T}}x_B+c_N^{\mathrm{T}}x_N\\ &=c_B^{\mathrm{T}}B^{-1}(b-Nx_N)+c_N^{\mathrm{T}}x_N\\ &=c_B^{\mathrm{T}}B^{-1}b+(c_N^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}N)x_N \end{aligned}上式表明 z 可仅由 n-m 个非基分量的值 x_N 线性表示。

注意到目标函数在 B 对应的基本可行解 x^{(0)} 处的取值是 z_0=c_B^{\mathrm{T}}B^{-1}b,有
\begin{aligned} z&=z_0+(c_N^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}N)x_N\\ &=z_0+(c_B^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}B)x_B+(c_N^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}N)x_N\\ &=z_0+(c^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}A)x\\ :&=z_0+\lambda^{\mathrm{T}} x \end{aligned}其中 \lambda^{T}=c^{\mathrm{T}}-c_B^{T}B^{-1}A\in\mathbb{R}^n检验向量,其各个分量称为检验数。尽管我们最终将 z 重新写成了 x 的线性函数,但实际上 x_B 的系数(m 个检验数)均为 0,因此 z 一共只有 n-m 个自由变量(即 x_N)。基于此,目标函数 z=z_0+\lambda^{\mathrm{T}}x 的这种形式也被称为简化形式

在可行域的限制之下,x\geqslant0。如果 \lambda 的每一个分量均非负(检验数均非负),则 z\geqslant z_0 恒成立,且 x=x^{(0)} 时等号成立(注意 x_N^{(0)}=0)。\lambda 各个分量的正负性成为了判别最优解的基本依据。

在介绍下面的定理前,先引入如下记号:\alpha=(\alpha_{ij})_{m\times n}=B^{-1}A, P_j'=B^{-1}P_j=(\alpha_{1j},\cdots,\alpha_{mj})^{\mathrm{T}}, (1\hspace{-1pt}\leqslant \hspace{-1pt}j\hspace{-1pt}\leqslant \hspace{-1pt}n), \beta=B^{-1}b.

定理 (最优性检验):
设可行基向量组 B 对应的基本解为 x^{(0)},目标函数的简化形式为 z=z_0+\lambda^{\mathrm{T}}x

  • 若所有的检验数 \lambda_k\geqslant 0,则 x^{(0)} 是最优解,z_0 是目标函数的最小值。
  • 若存在某检验数 \lambda_k<0,且 \alpha_{ik}\leqslant0~~(1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m),则无最优解。

证明:第一部分的证明已经在上面给出。考虑第二部分:
因为 \lambda_k\neq0,所以 x_k 是非基分量。取的 xk 分量 x_k=m>0,而其它非基分量均为 0。此时 x_N=(0,\cdots,m,\cdots,0)^{\mathrm{T}},代入 x_B 的表达式可得
x_B=B^{-1}(b-Nx_N)=\beta-B^{-1}P_km=\beta-P_k'm 由于 P_k'\leqslant0,故 x_B\geqslant0,从而 x 的各个分量非负,是 \mathcal{P} 的可行解。但此时 z=z_0+\lambda_kmm 可以取充分大的正数,因此目标函数无下界,不存在最优解。


总而言之,最优性检验利用目标函数的简化形式检验数的正负来判断当前基本可行解的最优性。在最优性检验定理当中,还有一种情况未被涉及:存在某检验数 \lambda_k<0,且存在下标 i 使得 \alpha_{ik}>0。这正是下面的基变换所要解决的问题。

基变换

在最优性检验当中,我们利用可行基向量组 B 和对应的基本解 x^{(0)} 得到来目标函数的简化形式 z=z_0+\lambda^{\mathrm{T}}x 和检验向量 \lambda。在最优性检验定理的未涉及的情形中,存在检验数 \lambda_k<0,并且 \alpha 的第 kP_k' 至少有一个分量为正。于是,我们希望替换基向量组中某一个基向量,得到一个新的可行基向量组,并且保证目标函数不增(然后重新评估新基对应的基本解的最优性)。这一步骤就是所谓的基变换

在进一步叙述基变换的具体操作之前,我们需要补充一些线性代数的知识。事实上,我更希望基变换的计算最后是由“纯粹的”向量和矩阵运算完成的。我喜欢赋予向量和矩阵鲜明的代数意义,通过简洁而又富有想象力的符号运算来获取结果,而讨厌与令人头疼的下标和臃肿的 n 阶矩阵展开式打交道。在我使用的《算法分析与设计》一书中,基变换的核心步骤就是以下标运算和矩阵展开写就的,矩阵的一些优美性质被掩盖了。因此,我决定使用更加富有技巧和表现力的矩阵运算来展现这一部分。


Gauss-Jordan 矩阵

n 为给定的正整数。记 e_l 是第 l 个分量为 1,其余分量为 0n 维列向量,即
e_l=(0,\cdots,\underbrace{1}_{l\mathrm{-th}},\cdots,0)^{\mathrm{T}}考察 n 阶矩阵 H=I-ye_l^{\mathrm{T}},其中 y=(y_1,\cdots,y_n)\in\mathbb{R}^n,则
H= \begin{pmatrix} 1& & &{}-y_1& & &\\[4pt] &\ddots&&\vdots&&&\\ &&1&{}-y_{l-1}\\[2pt] &&&1-y_l\\[2pt] &&&{}-y_{l+1}&1\\[4pt] &&&\vdots&&\ddots\\ &&&{}-y_n&&&1 \end{pmatrix}具有以上形式的矩阵称为Gauss-Jordan 矩阵。矩阵 H 拥有一些优良的性质:


有了以上的准备工作后,就可以开始进行基变换了。设原来的可行基向量组是
B=(P_{\pi(1)},\cdots,P_{\pi(l)},\cdots,P_{\pi(m)}) 如果将第 l 个基向量 P_{\pi(l)} 替换为非基向量 P_k~(k 是使得 \lambda_k<0 的那个下标),则得到的新的向量组为
B'=(P_{\pi(1)},\cdots,P_{k},\cdots,P_{\pi(m)}) 由于 B' 仅在第 l 列上与 B 不同,故可考虑 m 阶 Gauss-Jordan 矩阵 H=I-ye_l^\mathrm{T} 使得 B'=BH。于是
\begin{aligned} &B'=B(I-ye_l^\mathrm T)\\ \Rightarrow~&B'=B-Bye_l^\mathrm T\\ \Rightarrow~& Bye_l^\mathrm T=B-B'=(0,\cdots,P_{\pi(l)}-P_k,\cdots,0)\\ \Rightarrow~& By=P_{\pi(l)}-P_k\\ \Rightarrow~& y=e_l-P_k' \end{aligned} 根据最优性检验的假设,P_k' 至少有一分量为正。注意到 y_l=1-\alpha_{lk},因此当 \alpha_{lk}>0 时,一定有 H 可逆,故 B'=BH 可逆,从而 B'A 的一个基向量组


关于基向量组 B',需要证明以下结论:


最后给出基变换之后的参数变化情况:
y=e_l-P_k'H=I-ye_l^\mathrm T,则:


总而言之,基变换从可行基 B 出发构造了一组新的可行基 B',并且保证目标函数在基本解处的值不增。每一轮基变换之后,原先的参数 (\alpha,\beta,\lambda,z_0) 就会更新为 (\alpha',\beta',\lambda',z_0')

Bland 规则

单纯形法交替执行最优性检测基变换,Bland 规则提供了一种简单的办法使算法不会陷入循环:


综合以上的讨论,单纯形法的伪代码如下:

算法(单纯形法):Ax=bx\geqslant0 的条件下,求 c^\mathrm{T}x 的最小值。

  • 输入:A\in\mathbb{R}^{m\times n},~b\in\mathbb{R}^{m\times1},~c\in\mathbb{R}^{n\times1},~\pi\in\mathbb{R}^m
    其中 \pi 指定了 m 个基分量的下标。
  • 输出:\pi\in\mathbb{R}^m,~\beta\in\mathbb{R}^m,~z_0\in\mathbb{R}
    其中 \beta 是基本可行解的 m 个基分量,z_0 是目标函数最小值。

1. 初始化参数
B=A(:,\pi) % 基向量组,m*m
\alpha=B^{-1}A % m*n
\beta=B^{-1}b % 约化的基本可行解,m*1
c_B=c(\pi) % m*1
\lambda=c-\alpha^\mathrm Tc_B % 检验向量,n*1
z_0=c_B^\mathrm T\beta % 目标函数最小值


2. 最优性检验/基变换
找到最小的 k 使得 \lambda(k)<0

k 不存在:返回 \pi,~\beta,~z_0

k 存在:检查 \alpha(:,k) 中是否有正数

不存在:返回错误信息

存在:选择 l 使得 \displaystyle\frac{\beta(l)}{\alpha(l,k)}=\min\bigg\{\frac{\beta(i)}{\alpha(i,k)}\bigg|~\alpha(i,k)>0,~1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m \bigg\}
\pi(l)=k % 更换基分量
y=e_l-\alpha(:,k) % 辅助计算

z_0=z_0+\lambda(k)\cdot\dfrac{\beta(l)}{\alpha(l,k)} % 更新z0

\lambda=\lambda-\lambda(k)\dfrac{\alpha(l,:)^\mathrm T}{\alpha(l,k)} % 更新lambda

\beta=\beta+\dfrac{y\beta(l)}{\alpha(l,k)} % 更新beta

\alpha=\alpha+\dfrac{y\alpha(l,:)}{\alpha(l,k)} % 更新alpha

goto 2

上一篇 下一篇

猜你喜欢

热点阅读