线性规划问题与单纯形法

---

线性规划是运筹学中一个很重要的分支，本文通过一个实例简单的介绍一下什么是线性规划问题，以及单纯形法的计算步骤。

---

问题 ： 在一个工厂中，有两种产品A与B，它们的单价分别为6元与7元。A产品需要2千克原料P与4千克原料Q，而B产品需要3千克P与1千克Q，且P原料总数为16千克，Q原料总数为12千克，设计一种方案，即需要分别生产多少单位A和B可以使得利润最大化？
问题出处 ：《运筹学基础及其MATLAB应用》 – 李工农

---

我们分析上述问题，假设要生产x1单位A，x2单位B，那我们可以得到z = 6 * x1 + 7 * x2，所以我们的问题就转化为了怎样使得z最大。如果只按照上面这一式子求最大值，显然求不出来，故我们需要增加一些约束条件。因为P与Q原料是有限的，所以我们很容易得出2 * x1 + 3 * x2应该小与等于16，同时4 * x1 + 1 * x2应该小与等于12。此外，由于生产单位数不可能是负数，所以还应增加一约束条件x1 >= 0，x2 >= 0。我们将上面几个式子汇总一下，max表示求最大，s.t.(subject to)表示约束条件：
max z = 6 * x1 + 7 * x2
s.t.
2 * x1 + 3 * x2 <= 16
4 * x1 + 1 * x2 <= 12
x1,x2 >= 0
此时，我们再引入两个剩余参数x3与x4，目的是使得约束条件的前两个式子能够左右相等，由于前两个不等式左边恒小与等于右边，所以我们将x3与x4分别放到这两个式子左边，为了能够左右两边值相同，x3与x4必大于等于0。
max z = 6 * x1 + 7 * x2
s.t.
2 * x1 + 3 * x2 + x3 = 16
4 * x1 + 1 * x2 + x4 = 12
x1,x2,x3,x4 >= 0
上式就是一个标准的线性规划模型，现在给出线性规划的通俗定义，形如
max z = c^TX
s.t.
$\Sigma {}_{j=1}^n{P}_j{x}_j=b$
$x_j>=0(j=1,2,...,n)$
这类模型的问题都可以划为线性规划的范畴。其中，c^T等价于上面问题中的[6,7,0,0]，称作价值系数，X表示[x1,x2,x3,x4]^T，称作决策变量，b表示[16,12]^T，称作资源系数，最后P表示[2,3,1,0 ;4,1,0,1]，称作技术系数。另外我们称z为目标函数，我们需要求的便是目标函数的最优解。

---

到现在，我们已经大致明白了线性规划问题究竟是什么样的一个模型。然后，我们还需要知道一些线性规划中一些名词：

• 可行解：满足全部约束条件的解
• 可行域：可行解所构成的集合
• 最优解：所有可行解对应的目标函数值中，最大的一个值所对应的解
• 基：假设技术系数P是一个m x n的一个矩阵，其秩rank ( P ) = m，n >= m我们可以知道P中至少有m个线性无关的列向量，我们将这m个列向量组成一个新的矩阵，记作B，这就是一组基。例如问题中P = [2,3,1,0; 4,1,0,1]，则我们可以挑出一组线性无关的列向量（基）[1,0; 0,1]
• 基变量：基对应的X称作基变量$X_B$，例如[1,0; 0,1]对应的基变量就是[x3, x4]
• 非基：顾名思义P矩阵除去基后剩余的矩阵称作非基N，所以P = [B,N]
• 非基变量 ：X除去基变量后剩余的元素称作非基变量$X_N$，X = [$X_B,X_N$]
• 基解：我们设PX=b(即约束条件中的第一个式子)，其等价于(B,N)($X_B^T,X_N^T$)^T = b，化简等于$B{X}_B+N{X}_N=b$。此时我们令非基变量$X_N=0$，可得到$X_B=B^{-1}b$，这是我们称由基变量与非基
  变量合成的X，即($b^T{B}^{-T},0$)是问题的一个基解，显然取基不同，基解也不尽相同
• 基可行解：在基解中若满足X >= 0，称该基解为基可行解
• 可行基：基可行解对应的基称作可行基

---

了解了这些名词后，我们正式引入单纯形法解决这类问题。以下列出了单纯形法的基本步骤：

1. 先找出一个基可行解，称作初始基可行解
2. 判断该基可行解是否是最优解
3. 若判断为最优解算法终止，否则更换基可行解，返回第二步

我们通过解决最初提出的问题更详细的描述一下单纯形法。

1. 初始条件为P = [2,3,1,0 ; 4,1,0,1]，c = (6,7,0,0)^T，b = [16,12]^T，X = [x1,x2,x3,x4]^T
2. 先找基B，我们可以很容易的发现一对线性无关的单位向量[1,0 ; 0,1]，令其为B，N为[2,3 ; 4,1]，由于$B{X}_B+N{X}_N=b$，其中，$X_B=[x_3,x_4{]}^T,X_N=[x_1,x_2{]}^T$，令$X_N=0$，可以得到$X_B=[16,12{]}^T$，X = [0,0,16,12]^T，由于X >= 0，故该解是一个基可行解
3. 我们由约束方程组可以得到：
   $(1)x_3=16-2\ast x_1-3\ast x_2$
   $(2)x_4=12-4\ast x_1-x_2$
   将(1)式与(2)式代入目标函数可得:
   $(3)z=6\ast x_1+7\ast x_2+0\ast x_3+0\ast x_4=6\ast x_1+7\ast x_2$
   当$x_1与x_2$不在为0，且逐渐增大时，z也会随之增加，故不是最优解。其实，我们判断是否是最优解就是通过看目标函数中非基变量的系数是否有大于0的数来区分的，只要有大于0的数，增加相应的变量肯定会增加目标函数的值。
4. 更换基，此时我们需要确定换入元素与换出元素。换入元素我们选择目标函数中系数最大的对应元素，即$x_2$，由于$x_1$仍是非基变量，所以令其为0，我们可将(1)(2)式变为：
   $(4)x_3=16-3\ast x_2>=0$
   $(5)x_4=12-x_2>=0$
   我们选择令$x_3与x_4$等于0时所解得的$x_2$最小的那个变量作为换出变量，即$x_3$。由(1)式可得$(6)x_2=\frac{16}{3}-\frac{2}{3}\ast x_1-\frac{1}{3}\ast x_3$，将(6)式代入(2)式可得：
   $(6)x_2=\frac{16}{3}-\frac{2}{3}\ast x_1-\frac{1}{3}\ast x_3$
   $(7)x_4=\frac{20}{3}-\frac{10}{3}\ast x_1+\frac{1}{3}\ast x_3$
5. 将(6)(7)式代入目标函数中可得：
   $(8)z=6\ast x_1+7\ast (\frac{16}{3}-\frac{2}{3}\ast x_1-\frac{1}{3}\ast x_3)=\frac{112}{3}+\frac{4}{3}\ast x_1-\frac{7}{3}\ast x_3$，令$x_1与x_3$等于0可以得到一组解，按照第3步的判定规则我们认为该解也不是最优解
6. 参照第四步的换基规则确定换出元素与换入元素分别是$x_4与x_1$，且通过(6)(7)式可以得到：
   $(9)x_2=4+\frac{1}{5}\ast x_4-\frac{2}{5}\ast x_3$
   $(10)x_1=2-\frac{3}{10}\ast x_4+\frac{1}{10}\ast x_3$
7. 将(9)(10)式代入目标函数中可得：
   $(11)z=40-\frac{11}{5}\ast x_3-\frac{2}{5}\ast x_4$，令$x_3与x_4$等于0可以得到一组解，按照第3步的判定规则我们认为该解是最优解。最优解为[2,4,0,0]，值为40。

---

为了能更直观，更快速的计算，可以将上述内容做成如下一张表格（单纯形表）：

	$c_j$		6	7	0	0	$\theta$
$c_B$	$X_B$	b	x1	x2	x3	x4
0	$x_3$	16	2	[3]	1	0	$\frac{16}{3}$
0	$x_4$	12	4	1	0	1	12
${\sigma }_j=c_j-{\Sigma }_{i=1}^m{c}_i{a}_{ij}$			6	7	0	0
7	$x_2$	$\frac{16}{3}$	$\frac{2}{3}$	1	$\frac{1}{3}$	0	8
0	$x_4$	$\frac{20}{3}$	[$\frac{10}{3}$]	0	-$\frac{1}{3}$	1	2
${\sigma }_j=c_j-{\Sigma }_{i=1}^m{c}_i{a}_{ij}$			$\frac{4}{3}$	0	-$\frac{7}{3}$	0
7	$x_2$	4	0	1	$\frac{2}{5}$	- $\frac{1}{5}$
6	$x_1$	2	1	0	- $\frac{1}{10}$	$\frac{3}{10}$
${\sigma }_j=c_j-{\Sigma }_{i=1}^m{c}_i{a}_{ij}$			0	0	- $\frac{11}{5}$	- $\frac{2}{5}$

在这张表中，[]中的值称作旋转主元，其横向对应的变量为换出变量，纵向对应的为换入变量。在每次迭代中，最小的$\theta 与最大的\sigma$指向一个旋转主元。$\sigma$的计算是由$c_j-{\Sigma }_{i=1}^m{c}_i{a}_{ij}$计算得来的，举个例子，在第二次迭代中，${\sigma }_1=6-7\ast \frac{2}{3}-0\ast \frac{10}{3}$。而$\theta$的计算是b中元素分别除以最大$\sigma$所在列的元素。对于每次基变换其实只需要对增广矩阵[b,P]进行行列变换即可，其目的是使得换入变量所在列与之前其他单位向量继续能够构成单位矩阵。最后评判是否是最优解只需要看非基变量的$\theta$是否有大与零。

---

MATLAB源码点击此处下载