機器學習-線性規划(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/