机器学习-线性规划(LP)

 

线性规划问题

首先引入如下的问题：

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

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的钙，怎样买最省钱？

设买燕麦、全奶、草莓派、猪肉为x₁,x₂,x₃,x₄

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

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

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

可行域

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

max x1+x2

s.t 4x1–x2≤8

    2x1+x2≤10

    5x1–2x2≥−2

    x1,x2≥0

 

我们可以画出下面的图：

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

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

 

线性规划标准形式

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

min c^Tx

s.t  Ax≤b

     x≥0

就是

• 求的是min
• 所有的约束为<=的形式
• 所有的变量均 >=0

如何变为标准形式？

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

线性规划松弛形式

松弛形式为：

 min c^Tx

 s.t. Ax=b

      x≥0 

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

如： 

x1+x2≤2

x1+x2≥1

x1,x2≥0

 

写为松弛形式就是

x1+x2+x3=2

x1+x2+x4=1

x1,x2,x3,x4≥0

 <= vs <

为什么我们的线性规划的形式都是可以 <= 或者 >=的形式的？把等号去掉可以么？不可以

举个例子

 max x

 s.t. x≤1

 max x

 s.t. x<1

 显然第二个是无解的。

 

单纯形算法的思想与例子

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

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

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

它主要就三个步骤

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

以下面的线性规划为例:

min −x1−14x2−6x3

s.t. x1+x2+x3≤4

     x1≤2

     x3≤3

     3x2+x3≤6

     x1,x2,x3≥0

 

将其写为松弛的形式：

 

min −x1−14x2–6x3

s.t. x1+x2+x3+x4=4

     x1+x5=2

     x3++x6=3

     3x2+x3+x7=6

     x1,x2,x3,x4,x5,x6,x7≥0

 

其实，就是等价于（仍然要求 x₁,x₂,x₃,x₄,x₅,x₆,x₇ >=0）：

 

z=−x1−14x2–6x3

x4=4−x1–x2−x3

x5=2–x1

x6=3–x3

x7=6–3x2–x3

 

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

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

这里，容易得到基本解为：(x₁,x₂….x₇) = (0,0,0,4,2,3,6)，而目标值z = 0，其实就是把基本变量x_i设置为b_i。

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

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

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

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

z=−2−14x2–6x3+x5

x4=2–x2−x3+x5

x1=2–x5

x6=3–x3

x7=6–3x2–x3

 

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

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

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

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

z=−30+8x3+14x4−13x5

x2=2−x3−x4+x5

x1=2–x5

x6=3–x3

x7=2x3+3x4–3x5

 此时，我们的基本解变为(x₁,x₂….x₇) = (2,2,0,0,0,3,0)， Z = -30

我们可以继续的选择增大x₅，第4个等式具有最严格的限制（0 – 3x₅ >=0），我们有x₅ = 2/3 x₃ + x₄ – 1/3 x₇

带入得

 

z=−30–23x3+x4+133x7

x2=2–13x3−13x7

x1=2–23x3–x4+13x7

x6=3–x3

x5=23x3+x4–13x7

 

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

z=−32+x1+2x4+4x7

x2=1+12x1+12x4–12x7

x3=3–32x1–32x4+12x7

x6=32x1+32x4–12x7

x5=2–x1

 

现在，已经没有可以继续增大的值了，停止转动，z=-32就是我们的解，而此时，基本解为：(x₁,x₂….x₇) = (0,1,3,0,2,0,0)，看看最开始的目标函数：z = -x₁ -14x₂ – 6x₃ ,我们将x₂=1,x₃=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规则：

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

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

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

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

 

无界(unbounded)的情况

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

min−x1–x2

s.t. x1–x2≤1

     −x1+x2≤1

     x1,x2≥0

 

画出区域为：

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

上述的写成松弛形式为：

min −x1–x2

s.t.  x1–x2+x3=1

      −x1+x2+x4=1

      x1,x2,x3,x4≥0

 

也就是，

 

z=−x1–x2

x3=1–x1+x2

x4=1+x1–x2

 

选择x₁ 为替入变量，x₃为替出变量，有：

 

z=−1–2x2+x3

x1=1+x2–x3

x4=2–x3

 

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

 

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

 

从几何角度看单纯形算法

上面我们介绍单纯形算法的时候，是通过最直观的等式变换（就是旋转操作）介绍的。

我们知道，线性规划就是在可行域围成的多胞形中求解，现在从几何的视图来看看单纯形算法。

 

只需考虑顶点

让我再次召唤之前的图：

直观上看，最优解就在顶点上，不需要考虑内部点。

一个引入的证明

我们假设x⁽⁰⁾ 是最优解，连接x⁽¹⁾和x⁽⁰⁾ 与 x⁽²⁾和x⁽³⁾相交于点x’

我们可以把x⁽⁰⁾ 分解，x⁽⁰⁾ = λ₁ x⁽¹⁾ + (1 – λ₁)x’ 其中λ₁ = p / (p + q)

同样的把x‘ 分解，x’ = λ₂ x⁽²⁾ + (1 – λ₂)x⁽³⁾ 其中λ₂ = r / (r + s)

因此有：x⁽⁰⁾ = λ₁ x⁽¹⁾ + (1 – λ₁)λ₂ x⁽²⁾ + (1 – λ₁) (1 – λ₂)x⁽³⁾，而λ₁ + (1 – λ₁)λ₂ + (1 – λ₁) (1 – λ₂) = 1

设 c^T x⁽¹⁾ 小于等于 c^T x⁽²⁾， c^T x⁽³⁾，因此有：

 

CTx(0)=λ1CTx(1)+(1–λ1)λ2CTx(2)+(1–λ1)(1–λ2)CTx(3)

         ≥λ1CTx(1)+(1–λ1)λ2CTx(1)+(1–λ1)(1–λ2)CTx(1)

         =CTx(1)

 

因此，x⁽¹⁾ 并不比x⁽⁰⁾ 差。

小结

用几何的角度看待单纯形算法，主要有几点：

1. 最优解可以在顶点上找到，不许考虑内部点
2. 一次旋转就是一个顶点沿着一条边λ走θ倍到另一个顶点的过程

当然也需要注意初始化单纯形算法，比如之前的情况：

我们的顶点要在可行域才行，而不要跑到(0,0)去了。初始方法和之前的一样。

 

单纯形算法的调用(Python内置工具包)

python真的是非常强大。scipy包里面包含了很多科学计算相关的模块方法。

'''
原题目：
有2000元经费，需要采购单价为50元的若干桌子和单价为20元的若干椅子，你希望桌椅的总数尽可能的多，但要求椅子数量不少于桌子数量，且不多于桌子数量的1.5倍，那你需要怎样的一个采购方案呢？
解：要采购x1张桌子，x2把椅子

max z= x1 + x2
s.t. x1 - x2 <= 0
1.5x1 >= x2
50x1 + 20x2 <= 2000
x1, x2 >=0

在python中此类线性规划问题可用lp solver解决
scipy.optimize._linprog def linprog(c: int,
            A_ub: Optional[int] = None,
            b_ub: Optional[int] = None,
            A_eq: Optional[int] = None,
            b_eq: Optional[int] = None,
            bounds: Optional[Iterable] = None,
            method: Optional[str] = 'simplex',
            callback: Optional[Callable] = None,
            options: Optional[dict] = None) -> OptimizeResult

矩阵A：就是约束条件的系数（等号左边的系数）
矩阵B：就是约束条件的值（等号右边）
矩阵C：目标函数的系数值
'''

from scipy import  optimize as opt
import numpy as np
#参数
#c是目标函数里变量的系数
c=np.array([1,1])
#a是不等式条件的变量系数
a=np.array([[1,-1],[-1.5,1],[50,20]])
#b是是不等式条件的常数项
b=np.array([0,0,2000])
#a1，b1是等式条件的变量系数和常数项，这个例子里无等式条件,不要这两项
#a1=np.array([[1,1,1]])
#b1=np.array([7])
#限制
lim1=(0,None) #(0,None)->(0,+无穷)
lim2=(0,None)
#调用函数
ans=opt.linprog(-c,a,b,bounds=(lim1,lim2))
#输出结果
print(ans)

#注意：我们这里的应用问题，椅子不能是0.5把，所以最后应该采购37把椅子

手动自己实现一个简单的单纯性算法

import numpy as np

class Simplex(object):
    def __init__(self, obj, max_mode=False):
        self.max_mode = max_mode  # 默认是求min，如果是max目标函数要乘-1
        self.mat = np.array([[0] + obj]) * (-1 if max_mode else 1)      #矩阵先加入目标函数

    def add_constraint(self, a, b):
        self.mat = np.vstack([self.mat, [b] + a])      #矩阵加入约束

    def solve(self):
        m, n = self.mat.shape  # 矩阵里有1行目标函数，m - 1行约束，应加入m-1个松弛变量
        temp, B = np.vstack([np.zeros((1, m - 1)), np.eye(m - 1)]), list(range(n - 1, n + m - 1))  # temp是一个对角矩阵，B是个递增序列
        mat = self.mat = np.hstack([self.mat, temp])  # 横向拼接
        while mat[0, 1:].min() < 0:   #判断目标函数里是否还有负系数项
            col = np.where(mat[0, 1:] < 0)[0][0] + 1  # 在目标函数里找到第一个负系数的变量，找到替入变量
            row = np.array([mat[i][0] / mat[i][col] if mat[i][col] > 0 else 0x7fffffff for i in
                            range(1, mat.shape[0])]).argmin() + 1  # 找到最严格约束的行，也就找到替出变量
            if mat[row][col] <= 0: return None  # 若无替出变量，原问题无界
            mat[row] /= mat[row][col]    #替入变量和替出变量互换
            ids = np.arange(mat.shape[0]) != row
            mat[ids] -= mat[row] * mat[ids, col:col + 1]  # 对i!= row的每一行约束条件，做替入和替出变量的替换
            B[row] = col  #用B数组记录替入的替入变量
        return mat[0][0] * (1 if self.max_mode else -1), {B[i]: mat[i, 0] for i in range(1, m) if B[i] < n} #返回目标值，对应x的解就是基本变量为对应的bi，非基本变量为0


"""
       minimize -x1 - 14x2 - 6x3
       st
        x1 + x2 + x3 <=4
        x1 <= 2
        x3 <= 3
        3x2 + x3 <= 6
        x1 ,x2 ,x3 >= 0
       answer :-32
    """
t = Simplex([-1, -14, -6])
t.add_constraint([1, 1, 1], 4)
t.add_constraint([1, 0, 0], 2)
t.add_constraint([0, 0, 1], 3)
t.add_constraint([0, 3, 1], 6)
print(t.solve())
print(t.mat)

　　

小结

给定一个线性规划L，就只有如下三种情形：

1. 有一个有限目标值的最优解
2. 不可行
3. 无界

 

转载：细语呢喃 > 线性规划-单纯形算法详解
本文在转载的过程中加入了一些自己代码的实现，想看原作者详情的请移步到：https://www.hrwhisper.me/introduction-to-simplex-algorithm/