線性規劃-單純形算法詳解

本博客若無特殊說明則由 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