单纯形法求解线性规划

线性规划问题

线性规划问题是在一些线性约束下求解目标函数极值的问题.

在二维的平面上, 每个线性约束确定了一条直线的某一侧. 一些线性约束最终形成了一个区域. 从几何直观上很容易看出这个区域一定是凸的, 也就是任意两个可行点的连线上都是可行点.
目标函数对应了一组平行的直线, 移动这条直线使得其与可行范围有交点, 则直线对应的目标函数值 $z$ 就是可达到的. 所以在高中里求解二维线性规划问题就是画图后求解可行区域的角点对应的目标函数值, 从中选出最大或者最小的作为极值(在线性规划有解的情形下).

但是高维的线性规划问题无法直观画图, 也就不能用直观的办法解决. 一般形式的线性规划问题都可以转化为以下标准形式:

$$\begin{align*} \max & \quad z = c^T x \\ \text{s.t. }& \left\{ \begin{array}{ll} A x = b \ge 0, \\ x \ge 0, \end{array} \right. \end{align*}$$

其中 $c = (c_1, c_2, \cdots, c_n)^T$, $X = (x_1, x_2, \cdots, x_n)^T$, $b = (b_1, b_2, \cdots, b_m)^T$, $A = (a_{ij})_{mn}$ 且行向量线性无关.

将满足约束条件 $Ax=b$ 的集合称为 $S$,也称 $S$ 为可行区域, 对于其中任意 $x ,y \in S$ 有:

$$Ax=b, Ay=b$$

$x$ 和 $y$之间任意一个元素可以表示为 $v= \alpha x+(1-\alpha)y$,计算 $Av=A[\alpha x+(1-\alpha)y]=\alpha b+b-\alpha b=b$, 显然有 $v \in S$, 所以任意维的线性规划中的约束条件描述的都是一个凸集.

标准形式的转化

任意一个线性规划问题可以通过以下几个步骤转化为等价的标准形式.

目标函数的转化

统一求目标函数的极大值. 若是求 $z$ 的极小值, 则等价于求解 $-z$ 的极大值. 即：

$$\max z' = \Sigma_{j} c_j x_j$$

变量的转换

• (1) 本来就是要求 $x_j \ge 0$ 的不做变换;
• (2) 对于小于等于零的变量 $x_j$, 取负号令其变为大于等于零的变量, 即若 $x_j \le 0$, 则定义新变量 $x_j' = -x_j$, 约束条件为 $x_j' \ge 0$;
• (3) 如果约束条件为 $x_j \ge a$, 则定义新的变量 $x_j' = x_j - a$ 满足约束条件 $x_j' \ge 0$.
• (4) 若 $x_j$ 取值无约束,可令两个新的非负变量 $x_j', x_j''$, 然后用 $x_j = x_j' - x_j''$ 替换原问题中的 $x_j$.

约束条件的转换

引入松弛变量, 将所有不等式全部转换为等式:

• 对于 $\ge$ 的约束加入一个变量 $x_s$, $x_s \ge 0$;
• 对于 $\le$ 的约束加入一个变量 $-x_s$, $x_s \ge 0$;

在实际问题中它表示未被充分利用的资源或缺少的资源, 所以在引入模型后它们在目标函数中的系数均为零.

约束条件中如果有严格的不等号($<$ 或者 $>$), 则在优化过程中可以直接不考虑这些约束, 因为这些条件永远不可能是勒紧的, 而只是检验优化的得到的结果是否满足这些条件.

约束条件的右端项常数的转换

$b_i < 0$ 时, 将等式或不等式两端同乘$(-1)$.

人工变量的加入

为了构造出初始基变量, 约束条件还可能需要加上人工变量. 人工变量最终必须等于 $0$ 才能保持原问题性质不变. 为保证人工变量为0, 在目标函数中令其系数为 $M$. $M$ 为无限大的正数, 这是一个惩罚项, 倘若人工变量不为零,则目标函数就永远达不到最优, 所以必须将人工变量逐步从基变量中替换出去. 如若到最终表中人工变量仍没有置换出去, 那么这个问题就没有可行解, 当然亦无最优解.

解的存在性

假设不存在冗余的约束条件, 转换出的问题对应的矩阵 $A$ 的秩为 $m$, 且 $m \le n$, 这样方程才能有解.

单纯形法的思路

基变量与非基变量(自由变元)

n个变元满足m个约束条件, 则可以有 n-m 个自由变元. 选取 m 个基变量 $x_j'(j = 1,2,\cdots,m)$, 其对应约束系数矩阵的列向量线性无关. 通过线性行变换, 基变量可由非基变量表示:

$$x_i' = C_i + \Sigma_{j = m+1}^n m_{ij}x_j'\quad (i = 1,2,\cdots,m)$$

令非基变量等于 $0$, 可求得基变量的值:

$$x_i' = C_i$$

而目标函数也可以由非基变量表示:

$$z = z_0 + \Sigma_{j = m+1}^n \sigma_j x_j'$$

当达到最优解时, 所有的 $\sigma_j$ 应小于等于 $0$. 当存在 $\sigma_j > 0$时, 当前解不是最优解.

$\sigma_j$ 可以通过把每个基变量用非基变量和常数表示的形式代入 $z$ 的原始形式得出. 对于基变量, 不会出现在表达式中, 对应的 $\sigma_j = 0$. 对于非基变量, $\sigma_j$ 由两部分组成. 一部分是本身在 $z$ 中的系数, 直接由 $C$ 给出; 另一部分是通过基变量带来的系数, 也就是负的 基变量的系数乘以对应行中该非基变量参数的加和. 这部分文字描述比较困难, 在下面的例子中再解释.

若 $\sigma_j > 0$, $x_j' = 0$ 代表可行域的某个边界, 当脱离此边界变大时, 目标函数的值也会增大, 因此 $z_0$ 不是最优解.

如果存在 $\sigma_j > 0$, 考虑将其对应的变量换入基变量. 当存在多个 $\sigma > 0$ 时, 选择最大的一个 $\sigma$ 对应的变量作为基变量, 这表明目标函数随着 $x_j'$ 的增加而增长得最快.

随着 $x_j'$ 的增加其他 $x$ 可能会小于 $0$, 所以应该增大 $x_j'$ 直到有一个 $x$ 刚好等于 $0$ 为止, 然后将这个 $x$ 换出基变量. 具体的计算方法就是找到使得 $\frac{b_j}{a_ij}$ 最小的 $j$.

在线性行变换中, 把第 $j$ 个基变量对应的列向量变换成只有第 $j$ 个元素是 $1$, 其余都是 $0$ 的向量. 这样使用非基变量表示基变量时相当于得到一个单位矩阵的方程, 就不用再求解了(相当于之前做了高斯消元). 由于我们总是令非基变量为 $0$, 所以第 $j$ 个基变量的值就是第 $j$ 行对应的常数项.

选择被替换的基变量

假设选择非基变量 $x_s'$ 作为下一轮的基变量, 替换原先的基变量 $x_j'$. 选择 $x_j'$ 的原则是让替换后的 $x_s'$ 尽可能大.

终止条件

目标函数用非基变量表示时各个非基变量的系数 $\sigma_j$ 都小于等于 $0$, 则目标函数达到最优.

最优时如果存在非基变量的系数为 $0$, 则说明目标函数的梯度与某一边界正交, 所以目标函数的最优解有无穷多个.

---

上面这个思路很清晰, 但是细节上总是说不清, 再搬运一个…

单纯形法的算法步骤

• (1) 确定初始可行基和初始基可行解

  建立初始单纯形表;

• (2) 最优性检验

  若在当前表的目标函数对应的行中, 所有非基变量的系数非正, 则可判断得到最优解 (目标值不会再继续增大, 不会出现更优解), 可停止计算. 否则转入下一步;

• (3) 若单纯形表中 $1$ 至 $m$ 列构成单位矩阵, 在 $j=m+1$ 至 $n$ 列中, 若有某个对应 $x_k$ 的系数列向量 $P_k \le 0$, 则此问题是无界, 停止计算. 否则转入下一步;

• (4) 挑选目标函数对应行中系数最大的非基变量作为进基变量. 假设 $x_k$ 为进基变量, 按 $\theta$ 规则 [1] 计算, 可确定 $x_l$ 为出基变量, 转下一步;

$$\theta = \min \left( \frac{b_i}{a_{ik}} | a_{ik} >0 = \frac{b_l}{a_{lk}} \right) \tag{1}$$

• (5) 以 $a_lk$ 为主元素进行迭代(即用高斯消去法或称为旋转运算, 也就是初等行变换), 把 $x_k$ 所对应的列向量进行变换[2];

$$P_k = \begin{pmatrix} a_{1k} \\ a_{2k} \\ \vdots \\ a_{lk} \\ \vdots \\ a_{mk} \end{pmatrix} \Rightarrow \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 1 \\ \vdots \\ 0 \end{pmatrix} \tag{2}$$

也就是把 $x_k$ 用新的非基变量表出, 关键是要不含有其他基变量, 且其他基变量的表示中也不含有 $x_k$. 或者说基变量的列向量构成单位阵的一个重排.

• (6) 重复2-5步, 直到所有检验数非正后终止, 得到最优解.

注意每次迭代我们都直接让基变量的值取了系数向量. 这就要求系数向量永远都要是大于等于零的. 第一次选定初始可行基需要单独注意这一点. 在之后的迭代中由于我们总是控制选择最小的 $\frac{b_i}{a_{ik}}$, 这一点总是可以做到.

一个算例

求解线性规划问题

$$\begin{align*} \max & \quad z = x_1 + 14 x_2 + 6 x_3 \\ \text{s.t. }& \left\{ \begin{array}{ll} x_1 + x_2 + x_3 \le 4 \\ x_1 \le 2 \\ x_3 \le 3 \\ 3x_2 + x_3 \le 6 \\ x_1 , x_2 , x_3 \ge 0, \end{array} \right. \end{align*}$$

引入松弛变量, 转化为标准形式:

$$\begin{align*} \max & \quad z = x_1 + 14 x_2 + 6 x_3 \\ \text{s.t. }& \left\{ \begin{array}{ll} x_1 + x_2 + x_3 + x_4 = 4 \\ x_1 + x_5 = 2 \\ x_3 + x_6 = 3 \\ 3x_2 + x_3 + x_7 = 6 \\ x_1, x_2, x_3, x_4, x_5, x_6, x_7 \ge 0, \end{array} \right. \end{align*}$$

构造单纯形表如下:

		$C_j$	1	14	6	0	0	0	0
C	X	b	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$
0	$x_4$	4	1	1	1	1	0	0	0
0	$x_5$	2	1	0	0	0	1	0	0
0	$x_6$	3	0	0	1	0	0	1	0
0	$x_7$	6	0	3	1	0	0	0	1
	z		1	14	6	0	0	0	0

忽略前两列, 这张表展示的就是约束条件的 $m$ 个方程和最顶上一行目标函数. 找 $m$ 个变量作为一组基, 想通过初等行变换把它们用其他非基变量表示, 则矩阵中这 $m$ 个变量应该没有互相干扰, 其对应的列向量是单位矩阵的一个重排. 每一行的等式其实是一个基变量的表示方程. 在第二列记录各行对应的基变量, 第一列记录这个基变量在目标函数中的系数. 这里由于每个约束都引入了松弛变量, 直接就构成了单位阵, 直接就可以拿来做可行基.

非基变量 $(x_1, x_2, x_3)$ 的系数为 $(1,14,6)$, 因此可以通过增加这些变量来增长目标函数.

选取目标函数对应非基变量系数最大的 $x_2$ 换入基变量. 按照 $\theta$ 规则, 求得各个基变量对应的比值是 $(4, \infty, \infty, 2)$, 所以 $x_7$ 出基. 把第 $4$ 个式子整体除以 $3$, 得到 $x_2$ 用其他非基变量的表示形式, 然后使用初等行变换消去其他行中的 $x_2$, 得到新的单纯形表:

		$C_j$	1	14	6	0	0	0	0
C	X	b	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$
0	$x_4$	2	1	0	0.67	1	0	0	-0.33
0	$x_5$	2	1	0	0	0	1	0	0
0	$x_6$	3	0	0	1	0	0	1	0
14	$x_2$	2	0	1	0.33	0	0	0	0.33
	z		1	0	1.33	0	0	0	-4.67

这里 $z$ 用非基变量的表示可以如下计算. 对于 $x_1$, 原先的系数为 $1$, 减去通过 $x_4$ 带来的影响 $1 * 0$ 以及通过 $x_5$ 带来的影响 $1 * 0$, 得到 $1$. 对于 $x_3$, 原先的系数为 $6$, 减去通过 $x_4$ 带来的影响 $14 * 0.33$, 得到 $1.33$. 对于 $x_7$, 原先的系数为 $0$, 减去通过 $x_4$ 带来的影响 $14 * 0.33$, 得到 $-4.67$.

$x_3$ 进基, $x_4$ 出基, 再次计算, 得到:

		$C_j$	1	14	6	0	0	0	0
C	X	b	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$
6	$x_3$	3	1.5	0	1	1.5	0	0	-0.5
0	$x_5$	2	1	0	0	0	1	0	0
0	$x_6$	0	-1.5	0	0	-1.5	0	1	0.5
14	$x_2$	1	-0.5	1	0	-0.5	0	0	0.5
	z	-32	-1	0	0	-2	0	0	-4

此时我们发现所有非基变量的系数全部非正, 即增大任何基变量的值并不能使得目标函数增大. 于是我们可以断定该问题的最优解是 $z = 32, X = (0, 1, 3, 0, 2, 0, 0)$.