线性规划-单纯形算法详解

本博客若无特殊说明则由 hrwhisper 原创发布
转载请点名出处：细语呢喃 > 线性规划-单纯形算法详解
本文地址：https://www.hrwhisper.me/introduction-to-simplex-algorithm

线性规划-单纯形算法详解

 2016年11月24日   hrwhisper  7 Comments 15,246 views

本文将详细的介绍单纯形算法，包括但不限于

• LP问题
• 单纯形算法原理
• 无界、无解、循环等情况
• python代码实现

线性规划问题

首先引入如下的问题：

假设食物的各种营养成分、价格如下表：

Food	Energy（能量）	Protein（蛋白质）	Calcium（钙）	Price
Oatmeal（燕麦）	110	4	2	3
Whole milk（全奶）	160	8	285	9
Cherry pie（草莓派）	420	4	22	20
Pork with beans（猪肉）	260	14	80	19

要求我们买的食物中，至少要有2000的能量，55的蛋白质，800的钙，怎样买最省钱？

设买燕麦、全奶、草莓派、猪肉为x1,x2,x3,x4

于是我们可以写出如下的不等式组

其实这些不等式组就是线性规划方程（Linear programming formulation）：

简单的说，线性规划就是在给定限制的情况下，求解目标。

可行域

来看一个算法导论中的例子，考虑如下的线性规划：

maxs.t.x1+x24x1–x22x1+x25x1–2x2x1,x2≤≤≥≥810−20

我们可以画出下面的图：

看图a，灰色的区域就是这几个约束条件要求x1,x2所在的区域，而我们最后的解x1,x2也要在这里面。我们把这个区域称为可行域（feasible region）

图b可以直观的看出，最优解为8, 而 x1= 2 , x2=6

 

线性规划标准形式

线性规划的标准形式如下：

mins.t.AxxcTx≤b≥0

 

就是

• 求的是min（算法导论的是max，本文为min）
• 所有的约束为<=的形式
• 所有的变量均 >=0

如何变为标准形式？

• 原来是max, 直接*-1求min
• 若原来约束为=，转为 >= 和<=
• 约束原来为 >= 同样的*-1，就改变了<=
• 若有变量 xi < 0 ，那么用 x‘ – x”来替代，其中 x’>=0 x”>=0

 

线性规划松弛形式

松弛形式为：

mins.t.AxxcTx=b≥0

就是通过引入变量把原来的 <= ，变为=的松弛形式.

如：

x1+x2x1+x2x1,x2≤2≥1≥0

写为松弛形式就是

x1+x2+x3x1+x2+x4x1,x2,x3,x4=2=1≥0

<= vs <

有砸场子的同学会问(╯‵□′)╯︵┻━┻，为什么我们的线性规划的形式都是可以 <= 或者 >=的形式的？把等号去掉可以么？

就是不可以(￣ε(#￣)

举个例子

maxs.t.xx≤1

maxs.t.xx<1

显然第二个是无解的。

 

单纯形算法的思想与例子

如何求解线性规划问题呢？

有一些工具如GLPK，Gurobi 等，不在本文的介绍范围内。

本文要介绍的是单纯形算法，它是求解线性规划的经典方法，虽然它的执行时间在最坏的情况下是非多项式的（指数时间复杂度），但是，在绝大部分情况下或者说实际运行过程中却是多项式时间。

它主要就三个步骤

1. 找到一个初始的基本可行解
2. 不断的进行旋转（pivot）操作
3. 重复2直到结果不能改进为止

以下面的线性规划为例:

mins.t.−x1x1x1x1−+,14x2x23x2x2−++,6x3x3x3x3x3≤≤≤≤≥42360

将其写为松弛的形式：

mins.t.−x1−14x2–6x3x1+x2+x3+x4x1+x5x3++x63x2+x3+x7x1,x2,x3,x4,x5,x6,x7====≥04236

其实，就是等价于（仍然要求 x1,x2,x3,x4,x5,x6,x7 >=0）：

z=−x1−14x2–6x3x4=4−x1–x2−x3x5=2–x1x6=3–x3x7=6–3x2–x3

在上述的等式的左边称为基本变量，而右边称为非基本变量。

现在来考虑基本解就是把等式右边的所有非基本变量设为0，然后计算左边基本变量的值。

这里，容易得到基本解为：(x1,x2….x7) = (0,0,0,4,2,3,6)，而目标值z = 0，其实就是把基本变量xi设置为bi。

一般而言，基本解是可行的，我们称其为基本可行解。初始的基本解不可行的情况见后面的讨论，这里假设初始的基本解就是基本可行解，因此三个步骤中第一步完成了。

现在开始，来讨论上面的第二个步骤，就是旋转的操作。

我们每次选择一个在目标函数中的系数为负的非基本变量xe，然后尽可能的增加xe而不违反约束，并将xe用基本变量xl表示， 然后把xe变为基本变量，xl变为非基本变量。

这里，假设我们选择增加x1，那么在上述的等式（不包括目标函数z那行）中，第1个等式限制了x1 <=4（因为x4>=0），第2个等式有最严格的限制，它限制了x1 <=2，因此我们最多只能将x1增加到2，根据上面的第二个等式，我们有： x1 = 2 – x5，带入上面的等式就实现了xe和xl的替换：

z=−2−14x2–6x3+x5x4=2–x2−x3+x5x1=2–x5x6=3–x3x7=6–3x2–x3

这样其实就是一个转动(pivot)的过程，一次转动选取一个非基本变量（也叫替入变量）xe 和一个基本变量（也叫替出变量） xl ，然后替换二者的角色。执行一次转动的过程与之前所描述的线性规划是等价的。

同样的，将非基本变量设为0，于是得到：(x1,x2….x7) = (2,0,0,2,0,3,6)， Z = -2，说明我们的目标减少到了-2

接下来是单纯形算法的第三步，就是不断的进行转动，直到无法进行改进为止，继续看看刚才的例子：

我们接着再执行一次转动，这次我们可以选择增大x2或者x3，而不能选择x5，因为增大x5之后，z也增大，而我们要求的是最小化z。假设选择了x2，那么第1个等式限制了x2 <=2 , 第4个等式限制了x2 <= 2，假设我们选择x4为替出变量，于是有： x2 = 2 – x3 – x4 + x5 ，带入得：

z=−30+8x3+14x4−13x5x2=2−x3−x4+x5x1=2–x5x6=3–x3x7=2x3+3x4–3x5

此时，我们的基本解变为(x1,x2….x7) = (2,2,0,0,0,3,0)， Z = -30

我们可以继续的选择增大x5，第4个等式具有最严格的限制（0 – 3x5 >=0），我们有x5 = 2/3 x3 + x4 – 1/3 x7

带入得

z=−30–23x3+x4+133x7x2=2–13x3−13x7x1=2–23x3–x4+13x7x6=3–x3x5=23x3+x4–13x7

此时，我们的基本解变为(x1,x2….x7) = (2,2,0,0,0,3,0)， Z = -30，这时候并没有增加，但是下一步，我们可以选择增加 x3。第2个和第3个有最严格的限制，我们选第2个的话，得：x3 = 3 – 3/2 x1 – 3/2 x4 + 1/2 x7，然后老样子，继续带入：

z=−32+x1+2x4+4x7x2=1+12x1+12x4–12x7x3=3–32x1–32x4+12x7x6=32x1+32x4–12x7x5=2–x1

现在，已经没有可以继续增大的值了，停止转动，z=-32就是我们的解，而此时，基本解为：(x1,x2….x7) = (0,1,3,0,2,0,0)，看看最开始的目标函数：z = -x1 -14x2 – 6x3 ,我们将x2=1,x3=3带入得，z=-32，说明我们经过一系列的旋转，最后得到了目标值。

 

退化(Degeneracy)

在旋转的过程中，可能会存在保持目标值不变的情况，这种现象称为退化。比如上面的例子中，两次等于-30.

可以说退化可能会导致循环（cycling）的情况，这是使得单纯形算法不会终止的唯一原因。还好上面的例子中，我们没有产生循环的情况，再次旋转，目标值继续降低。

《算法导论》是这样介绍退化产生循环的：

Degeneracy can prevent the simplex algorithm from terminating, because it can lead to a phenomenon known as cycling: the slack forms at two different iterations of SIMPLEX are identical. Because of degeneracy, SIMPLEX could choose a sequence of pivot operations that leave the objective value unchanged but repeat a slack form within the sequence. Since SIMPLEX is a deterministic algorithm, if it cycles, then it will cycle through the same series of slack forms forever, never terminating.

如何避免退化？一个方法就是使用Bland规则：

在选择替入变量和替出变量的时候，我们总是选择满足条件的下标最小值。

• 替入变量xe：目标条件中，系数为负数的第一个作为替入变量
• 替出变量xl：对所有的约束条件中，选择对xe约束最紧的第一个

在上面的例子中，我也是这么做的。^ ^

另一个方法是加入随机扰动。

 

无界(unbounded)的情况

有的线性规划问题是无界的，举个栗子

对于下面的线性规划

mins.t.−−x1–x2x1–x2x1+x2x1,x2≤1≤1≥0

画出区域为：

显然可以不断的增大。让我们来看看单纯形算法是如何应对的：

上述的写成松弛形式为：

mins.t.−−x1–x2x1–x2+x1+x2x1,x2,x3,x4≥0x3+x4=1=1

也就是，

zx3x4===−x1–x21–x1+x21+x1–x2

选择x1 为替入变量，x3为替出变量，有：

zx1x4===−1–2x2+x31+x2–x32–x3

这时候我们只能选择x2 为替入变量,才能使得目标值变小，但是我们发现，对于x2没有任何的约束，也就是说，x2可以无限大，所以这是没有边界的情况。

这个情况是我们有一个替入变量，但是找不到一个替出变量导致的，这时候就是无界的情况了，写算法的时候注意判断一下即可。

 

单纯形算法的具体实现

说了那么多，代码怎么写呢？

看一下最开始的线性规划的问题（已经是松弛形式）：

mins.t.−x1−14x2–6x3x1+x2+x3+x4x1+x5x3++x63x2+x3+x7x1,x2,x3,x4,x5,x6,x7====≥04236

我们可以得到下面的矩阵：

C=(−1−14−60000)

我们可以得到下面的矩阵：

C=(−1−14−60000)B=⎛⎝⎜⎜000)B=⎛⎝⎜⎜⎜4236⎞⎠⎟⎟⎟A=42