模板
教程
環境
性能
功能
管理
GNU線性編程中的中間問題
本文將繼續介紹有關 GNU 線性編程工具包（GNU Linear Programming Kit） 和 glpsol 客戶機工具以及 GNU MathProg 語言的使用。在本文中，我們將從一個日常飲食問題入手來介紹如何表述一個簡單的多類型變量，并聲明二元參數。然后通過一個郵局資源分配問題來介紹 MathProg 表達式和只使用整型的決策變量。
本文是有關 GNU Linear Programming Kit（GLPK） 的 3 部分系列文章 中的第二篇。有關 GLPK 的介紹，請閱讀本系列的第一部分 “GNU 線性編程工具包（GNU Linear Programming Kit），第 1 部分: 線性優化簡介”。
這個問題來自于Wayne L. Winston 所著的 Operations Research: Applications and Algorithms, 4th Edition（Thomson，2004 年）；參閱后面 參考資料 中的鏈接。
我的日常飲食要求吃的所有食物都必須屬于四種基本的食物：巧克力蛋糕、冰激凌、汽水飲料和芝士蛋糕?，F在有 4 種食品可以選擇：巧克力松糕、巧克力冰激凌、可樂和菠蘿芝士蛋糕。每個巧克力松糕價值 0.50 美元，每份巧克力冰激凌價值 0.20 美元，每瓶可樂價值 0.30 美元，每片菠蘿芝士蛋糕價值 0.80 美元。每天我必須攝取至少 500 卡路里能量、6 盎司巧克力、10 盎司糖，以及 8 盎司脂肪。每種食物的單位營養含量如下表所示?，F在我想以最小成本滿足自己的日常營養需求。
現在我們對這個問題的相關重要信息進行一下總結：
食物的單位成本
- 巧克力松糕： $0.50 （每片）
- 巧克力冰激凌：$0.20 （每份）
- 可樂： $0.30 （每瓶）
- 菠蘿芝士蛋糕：$0.80 （每片）
每天的食物需求
- 500 卡路里
- 6 盎司巧克力
- 10 盎司糖
- 8 盎司脂肪
食物營養含量
- 巧克力松糕：400 卡路里能量，3 盎司巧克力，2 盎司糖，2 盎司脂肪
- 巧克力冰激凌（每份）：200 卡路里能量，2 盎司巧克力，2 盎司糖，4 盎司脂肪
- 可樂（每瓶）：150 卡路里能量，0 盎司巧克力，4 盎司糖，1 盎司脂肪
- 菠蘿芝士蛋糕（1 片）：500 卡路里能量，0 盎司巧克力，4 盎司糖，5 盎司脂肪
對這個日常飲食問題進行建模在解決過 第 1 部分中的 Giapetto 問題 之后應該就很容易了。下面讓我們首先來確定一下決策變量。
此處的目標是以最小的成本滿足日常營養需要。因此，決策變量應該是需要購買的每種食品的數量：
食品變量：
x1：巧克力松糕片數x2：巧克力冰激凌的份數x3：可樂的瓶數x4：菠蘿芝士蛋糕的片數
現在我們就可以使用這些決策變量來編寫目標函數了。日常飲食的成本 z 需要最小化：
下一個步驟是為各個約束編寫不等式。注意日常飲食的食物提供的卡路里、巧克力、糖和脂肪數量都有限制。這 4 種需要每個都是一個限制，因此約束可以如下所示：
注意菠蘿芝士蛋糕和可樂中都不含巧克力。
最后，是一些符號約束：
試想一下這個問題的可行區域。這個問題有 4 個變量。因此，其解析域應該有 4 個軸。這很難繪制出來，因此我們需要使用自己的想象力。首先，解析域應該在一個 4 維空間中。使用每個約束，解析空間會逐漸收縮，最終看起來就像是一個多面體一樣。
注意：清單 1 中的行號僅僅是為了引用方便而給出的。
清單 1. 日常飲食問題的 GNU MathProg 解決方案
1 #
2 # Diet problem
3 #
4 # This finds the optimal solution for minimizing the cost of my diet
5 #
6
7 /* sets */
8 set FOOD;
9 set NEED;
10
11 /* parameters */
12 param NutrTable {i in FOOD, j in NEED};
13 param Cost {i in FOOD};
14 param Need {j in NEED};
15
16 /* decision variables: x1: Brownie, x2: Ice cream, x3: soda, x4: cake*/
17 var x {i in FOOD} >= 0;
18
19 /* objective function */
20 minimize z: sum{i in FOOD} Cost[i]*x[i];
21
22 /* Constraints */
23 s.t. const{j in NEED} : sum{i in FOOD} NutrTable[i,j]*x[i] >= Need[j];
24
25 /* data section */
26 data;
27
28 set FOOD := Brownie "Ice cream" soda cake;
29 set NEED := Calorie Chocolate Sugar Fat;
30
31 param NutrTable: Calorie Chocolate Sugar Fat:=
32 Brownie 400 3 2 2
33 "Ice cream" 200 2 2 4
34 soda 150 0 4 1
35 cake 500 0 4 5;
36
37 param Cost:=
38 Brownie 0.5
39 "Ice cream" 0.2
40 soda 0.3
41 cake 0.8;
42
43 param Need:=
44 Calorie 500
45 Chocolate 6
46 Sugar 10
47 Fat 8;
48
49 end;
|
第 8 行和第 9 行定義了兩個集合：FOOD 和 NEED，但是這兩個集合的元素都是在第 28 行和 29 行的數據部分聲明的。FOOD 集合中包含了前面給出的 4 種食物約束。
由于空格字符被用來分隔集合中的元素，因此 Ice cream 元素（這里沒有使用 icecream）需要在名字兩邊使用雙引號括起來。如果我們希望在 MathProg 名字中使用非 ASCII 字符，那么即使名字中間沒有空格，也應該使用雙引號將它們括起來。
現在我們回到模型部分上來。NEED 集合中保存了 4 種飲食的需要。第 12、13 和 14 行定義了問題的參數：Need、 Cost 和 NutrTable（營養表）。每個 FOOD 都有一個成本值。對于 NEED 集合中的每種營養都有一定的需求量。我們可以試圖為每個集合使用不同的命名索引變量；這是一個好主意，尤其是在進行調試時更為突出。在這種情況中，FOOD 集合使用 i，而 NEED 集合使用 j。數據部分中的 Cost 和 Need 參數是在 37 行到 47 行進行定義的。
在 12 行上定義的營養表和在 31 到 35 行給出的數據有兩個維度，這是因為每種食物都提供了多種營養價值。因此，營養表 NutrTable 就是 FOOD 和 NEED 集合之間的一個映射。注意行和列的次序與元素在每個集合中的順序是相同的，索引變量名在第 12、13 和 14 行是一致的。
在這個練習中，i 軸是行，j 軸是列，這符合大部分數學專才的習慣。這個兩維參數聲明（最多 N 列 M 行）的語法如下：
param label : J_SET_VAL_1 J_SET_VAL_2 ... J_SET_VAL_N :=
I_SET_VAL_1 param_1_1 param_1_2 ... param_1_N
I_SET_VAL_2 param_2_1 param_2_2 ... param_2_N
... ... ... ... ...
I_SET_VAL_M param_M_1 param_M_2 ... param_M_N;
|
不要忘記第一行末尾的 :=以及最后一行末尾的 ; 符號。
第 17 行聲明了決策變量，并聲明每個決策變量都不能是負數。這是一個非常簡單的定義，因為數組 x 是使用 FOOD 集合中的元素自動定義的。
第 20 行的目標函數要對食物的總體成本實現最小化，該值是每個決策變量（食物數量）乘以每單位食物成本的結果。注意 i 是 FOOD 集合的索引。
第 23 行聲明了所有這 4 種食品的約束。這是采用非常精簡的風格來編寫的，因此需要特別注意！冒號 : 左邊的定義告訴 GLPK 為 NEED 集合中的每種需要創建一個約束： const[Calorie]、 const[Chocolate]、 const[Sugar] 和 const[Fat]。每個約束都要有自己的營養需求，每種食品的數量乘以每單位該食品所提供的營養需要的總和就是這種營養的總量；這個總和應該大于或等于日常飲食對該種營養的最小需求。
展開之后，這個聲明應該如下所示（i 會遍歷整個 FOOD 集合）：
Problem: diet
Rows: 5
Columns: 4
Non-zeros: 18
Status: OPTIMAL
Objective: z = 0.9 (MINimum)
No. Row name St Activity Lower bound Upper bound Marginal
------ ------------ -- ------------- ------------- ------------- -------------
1 z B 0.9
2 const[Calorie]
B 750 500
3 const[Chocolate]
NL 6 6 0.025
4 const[Sugar] NL 10 10 0.075
5 const[Fat] B 13 8
No. Column name St Activity Lower bound Upper bound Marginal
------ ------------ -- ------------- ------------- ------------- -------------
1 x[Brownie] NL 0 0 0.275
2 x['Ice cream']
B 3 0
3 x[soda] B 1 0
4 x[cake] NL 0 0 0.5
Karush-Kuhn-Tucker optimality conditions:
KKT.PE: max.abs.err. = 1.82e-13 on row 2
max.rel.err. = 2.43e-16 on row 2
High quality
KKT.PB: max.abs.err. = 0.00e+00 on row 0
max.rel.err. = 0.00e+00 on row 0
High quality
KKT.DE: max.abs.err. = 5.55e-17 on column 3
max.rel.err. = 4.27e-17 on column 3
High quality
KKT.DB: max.abs.err. = 0.00e+00 on row 0
max.rel.err. = 0.00e+00 on row 0
High quality
End of output
|
這些結果說明日常飲食最低成本（優化值）是 0.90 美元。哪些約束共同決定了這個解決方案呢？
這個報告的第二部分表明巧克力和糖的約束都是有下界的，因此這份日常飲食使用了最少的巧克力和糖。這個臨界值告訴我們如果我們可以放松巧克力限制一個單位（變成 5 盎司，而不是 6 盎司），那么目標函數就可以改進 0.025（它會從 0.9 變成 0.875）。類似地，如果我們將糖約束放松一個單位，那么目標函數就會改進 0.075。道理是很顯然的：吃得越少，我們花的錢也就越少。重要的一點是要對它們進行邊界和數量的常規意義的考察。例如，如果我們被告知最好吃 200 磅的巧克力，但是不能攝取任何能量，那么我們就會對此表示懷疑（如果真能如此，我們倒是會感激不盡）。
報告的第三部分則給出了決策變量的優化值：3 份冰激凌和 1 瓶可樂。巧克力松糕和菠蘿芝士蛋糕都有一個臨界值，因為這些值受到了符號約束的限制。如果巧克力松糕變量的臨界值可以是 -1，那么目標函數就還可以改進 0.275，但是對于這個問題的具體情況來說，這當然沒什么用處。
這個問題引自于“Operations Research”：
一個郵局需要有不同數目的全職員工在每周的不同時間工作（總結如下）。聯盟規定每個全職員工必須連續工作 5 天，然后休息 2 天。例如，在周一到周五工作的員工必須在周六和周日休息。郵局希望只雇傭全職員工來滿足自己的日常需求，并且雇傭員工的人數要最少。
下面我們對這個問題的一些重要信息進行一下總結：
- 每個全職工人只能連續工作 5 天，然后就要休息 2 天
- 第一天（周一）：需要 17 個工人
- 第二天：需要 13 個工人
- 第三天：需要 15 個工人
- 第四天：需要 19 個工人
- 第五天：需要 14 個工人
- 第六天：需要 16 個工人
- 第七天（周日）：需要 11 個工人
郵局需要對滿足自己需要的雇員數目實現最小化。
下面讓我們開始分析這個問題的決策變量。我們應該使用 7 個變量，一周中的每天都要使用一個變量，其值等于在當天工作的員工數目。盡管乍看起來這已經解決了這個問題，但是這不能實現一個員工每周只能工作 5 天的約束，因為在員工某天工作并不能要求該員工在下一天也工作。
正確的途徑應該是確保在 i 天開始工作的員工在接下來的連續 4 天也會工作，因此正確的方法是使用 xi 表示從 i 天開始工作的員工數目。使用這種方法，強制這種連續約束就簡單多了。決策變量就變成了：
x1：從周一開始工作的員工數目x2：從周二開始工作的員工數目x3：從周三開始工作的員工數目x4：從周四開始工作的員工數目x5：從周五開始工作的員工數目x6：從周六開始工作的員工數目x7：從周日開始工作的員工數目
需要最小化的目標函數是所雇傭員工的數量，它可以這樣給出：
那么，約束都是什么呢？一周中的每天都有一個約束，這是為了確保當天的工人數量最少。讓我們以周一為例來看一下。哪些人在周一工作呢？在我腦海中浮現出來的第一個（片面）答案是 “那些在周一開始工作的人”。但是還有別人嗎？有，那些要連續工作 5 天的員工中，周日開始工作的員工在周一時應該也在工作（回想一下問題的定義）。同理，我們可以推論那些周六、周五、周四開始工作的員工在周一也都在工作。
這個約束確保周一至少有 17 名員工在工作。
類似地：
當然，不要忘記了符號約束：
注意：清單 4 中的行號僅僅是為了參考方便而給出的。
1 #
2 # Post office problem
3 #
4 # This finds the optimal solution for minimizing the number of full-time
5 # employees to the post office problem
6 #
7
8 /* sets */
9 set DAYS;
10
11 /* parameters */
12 param Need {i in DAYS};
13
14 /* Decision variables. x[i]: No. of workers starting at day i */
15 var x {i in DAYS} >= 0;
16
17 /* objective function */
18 minimize z: sum{i in DAYS} x[i];
19
20 /* Constraints */
21
22 s.t. mon: sum{i in DAYS: i<>'Tue' and i<>'Wed'} x[i] >= Need['Mon'];
23 s.t. tue: sum{i in DAYS: i<>'Wed' and i<>'Thu'} x[i] >= Need['Tue'];
24 s.t. wed: sum{i in DAYS: i<>'Thu' and i<>'Fri'} x[i] >= Need['Wed'];
25 s.t. thu: sum{i in DAYS: i<>'Fri' and i<>'Sat'} x[i] >= Need['Thu'];
26 s.t. fri: sum{i in DAYS: i<>'Sat' and i<>'Sun'} x[i] >= Need['Fri'];
27 s.t. sat: sum{i in DAYS: i<>'Sun' and i<>'Mon'} x[i] >= Need['Sat'];
28 s.t. sun: sum{i in DAYS: i<>'Mon' and i<>'Tue'} x[i] >= Need['Sun'];
29
30 data;
31
32 set DAYS:= Mon Tue Wed Thu Fri Sat Sun;
33
34 param Need:=
35 Mon 17
36 Tue 13
37 Wed 15
38 Thu 19
39 Fri 14
40 Sat 16
41 Sun 11;
42
43 end;
|
第 9 行聲明了一個名為 DAYS 的集合，其元素（本周中的天數，從周一開始）是在數據部分的第 32 行中聲明的。
第 12 行為 DAYS 中的每天都聲明了 Need 參數。第 34 行到 41 行定義了這個參數的值：每周的每天都所需要的最少員工數。
第 15 行將決策變量聲明為一個數組，它有七個變量，索引在 DAYS 集合中定義；分別表示從該天開始工作的員工數目。
第 22 行到 28 行為每周的每天定義了一個約束。注意編寫 7 個不等式作為 5 個不必有序的決策變量的和太過繁瑣了，因為在某些約束中，索引可能會覆蓋索引值 7。 GNU MathProg 為編程者提供了 表達式 來簡化這個問題。
每個約束都是所有這些決策變量中除去那兩天不工作的特殊日期之外（而不是直接包括工作的這 5 天）的綜合結果。這個表達式在花括號（{}）中使用，它定義了和的索引。表達式的語法如下：
{index_variable in your_set: your_expression}
這個表達式可以使用邏輯比較操作。在本例中，Monday 約束使用了：i<>'Tue' 和 i<>'Wed'，這表示 “當 i 不等于 Tue 并且 i 不等于 Wed 時”。對于其他約束來說全部類似。
== 邏輯比較符也可以在這些表達式中使用。
Problem: post
Rows: 8
Columns: 7
Non-zeros: 42
Status: OPTIMAL
Objective: z = 22.33333333 (MINimum)
No. Row name St Activity Lower bound Upper bound Marginal
------ ------------ -- ------------- ------------- ------------- -------------
1 z B 22.3333
2 mon NL 17 17 0.333333
3 tue B 15 13
4 wed NL 15 15 0.333333
5 thu NL 19 19 0.333333
6 fri NL 14 14 < eps
7 sat NL 16 16 0.333333
8 sun B 15.6667 11
No. Column name St Activity Lower bound Upper bound Marginal
------ ------------ -- ------------- ------------- ------------- -------------
1 x[Mon] B 1.33333 0
2 x[Tue] B 5.33333 0
3 x[Wed] NL 0 0 < eps
4 x[Thu] B 7.33333 0
5 x[Fri] NL 0 0 0.333333
6 x[Sat] B 3.33333 0
7 x[Sun] B 5 0
Karush-Kuhn-Tucker optimality conditions:
KKT.PE: max.abs.err. = 3.55e-15 on row 6
max.rel.err. = 2.37e-16 on row 6
High quality
KKT.PB: max.abs.err. = 0.00e+00 on row 0
max.rel.err. = 0.00e+00 on row 0
High quality
KKT.DE: max.abs.err. = 5.55e-17 on column 1
max.rel.err. = 2.78e-17 on column 1
High quality
KKT.DB: max.abs.err. = 5.55e-17 on row 6
max.rel.err. = 5.55e-17 on row 6
High quality
End of output
|
咳，等會兒！沒有人可以讓 1.33333 個員工在周一開始工作！記住前面說過的用常識檢查，這就是其中的一個例子。
GLPK 必須要將這些決策變量全部當作整型變量進行考慮。幸運的是，MathProg 有一種很好的方法來聲明整型變量。我們只需要將第 15 行修改成下面的形式：
var x {i in DAYS} >=0, integer;
這相當簡單。glpsol 的輸出對整型情況的結果稍有不同：
Reading model section from post-office-int.mod...
Reading data section from post-office-int.mod...
50 lines were read
Generating z...
Generating mon...
Generating tue...
Generating wed...
Generating thu...
Generating fri...
Generating sat...
Generating sun...
Model has been suclearcase/" target="_blank" >ccessfully generated
lpx_simplex: original LP has 8 rows, 7 columns, 42 non-zeros
lpx_simplex: presolved LP has 7 rows, 7 columns, 35 non-zeros
lpx_adv_basis: size of triangular part = 7
0: objval = 0.000000000e+00 infeas = 1.000000000e+00 (0)
7: objval = 2.600000000e+01 infeas = 0.000000000e+00 (0)
* 7: objval = 2.600000000e+01 infeas = 0.000000000e+00 (0)
* 10: objval = 2.233333333e+01 infeas = 0.000000000e+00 (0)
OPTIMAL SOLUTION FOUND
Integer optimization begins...
Objective function is integral
+ 10: mip = not found yet >= -inf (1; 0)
+ 19: mip = 2.300000000e+01 >= 2.300000000e+01 0.0% (9; 0)
+ 19: mip = 2.300000000e+01 >= tree is empty 0.0% (0; 17)
INTEGER OPTIMAL SOLUTION FOUND
Time used: 0.0 secs
Memory used: 0.2M (175512 bytes)
lpx_print_mip: writing MIP problem solution to `post-office-int.sol'...
|
注意輸出結果現在顯示已經找到一個整型優化解決方案，在找到這個解決方案之前，GLPK 已經對放松限制下的問題（不要求變量是整數）計算了優化解決方案。
Problem: post
Rows: 8
Columns: 7 (7 integer, 0 binary)
Non-zeros: 42
Status: INTEGER OPTIMAL
Objective: z = 23 (MINimum) 22.33333333 (LP)
No. Row name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 z 23
2 mon 18 17
3 tue 13 13
4 wed 15 15
5 thu 19 19
6 fri 14 14
7 sat 16 16
8 sun 20 11
No. Column name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 x[Mon] * 1 0
2 x[Tue] * 2 0
3 x[Wed] * 3 0
4 x[Thu] * 7 0
5 x[Fri] * 1 0
6 x[Sat] * 3 0
7 x[Sun] * 6 0
End of output
|
第一部分說明所找到的解決方案是整型優化的，為目標函數找到的最小值是 23。這個郵局需要雇傭 23 名全職員工才能滿足自己的需求。放松限制的目標函數（不將決策變量當作整數考慮）的優化值同時也給出來了。
讓我們暫時跳過這個報告的第二部分。第三部分給出了決策變量的值。這些值都是能使整個問題的目標函數最小化的整數值。
現在，讓我們來討論一下有關第二部分的問題。它給出了約束的行為。有些是有下界的，現在我們可能會期望一個臨界值或影子價格。不過，討論整型問題的臨界值并沒有意義。整型問題的可行域并不是一個連續區域。換而言之，可行域不是一個多面體；它是由這個多面體或放松限制問題的實際多面體邊界中的一些整數（x1、x2、...、 xn）對。這意味著可行域是由這個空間中的一些離散點構成的，因此放松某個限制可能會在新的多面體中獲得更好的解決方案，也可能并不能獲得更好的解決方案。
要更好地理解此處討論的問題，讓我們來了解一下這個簡單的可行域：
使用 x 表示的藍點是一些整數（x1,x2）對。這是一個 x1 X x2 兩維空間的整數可行域，其約束為 x1 >= 0、x2 >= 0 和 x1 + x2 <= 4。如果這種簡單情況的目標函數是 maximize z = x1 - x2，那么顯然優化解決方案就是 (2,0)，這個約束剛好是邊界約束（因為優化解決方案就在約束規則這條線上）。
如果這個約束放松一個單元，x1 + x2 <= 5，可行域現在就不同了。
優化解決方案仍然是整數點 (2,0)，不過現在可行域中有更多整數點了。
因此，整型問題的約束放松并不一定會改進解決方案，因為可行域是離散的，而不是連續的。
您可能會問的另外一個問題是：“整型問題的優化解決方案與放松限制問題的優化解決方案有關聯嗎？”答案要看 simplex 算法背后采用的代數關系，但是有關它是如何工作的解釋已經超出了本文的范圍。我們只需要知道非整型問題的優化解決方案通常都是一個多面體的頂點就好了。通常這都是正確的！
在上面的第一個可行域中，優化解決方案是解析空間中的右頂點，而這個解析空間是一個由所有約束構成的三角形。在目標函數延伸的這個方向上，目標函數會逐漸改進。在這個簡單的情況中，x1 - x2 的方向延伸是 (1,-1)。因此，優化解決方案是從軸的原點沿 (1,-1) 方向作為方向向量指向多面體邊界的最遠點。由于整型解決方案可能會在一個邊界上，也可能在多面體內部，因此最佳的情況是采用一個頂點上的整型解決方案。在這種特例中，整型和非整型的優化解決方案是相同的，但是并非所有情況都是如此。為什么呢？因為這個整數點距離原點的距離可能不如放松限制的解決方案點那樣遠。對于其他非最佳情況來說，最好的整數點就在多面體中，目標函數自然比放松限制的問題的的性能要差，這我們在郵局問題中就已經注意到了。
記住這個分析使用了一個簡單的兩變量最大化問題。對于最小化問題的相同分析會查找可以將目標函數最小化的點，是最靠近原點的地方（沿著目標函數延伸方向的反方向）。
在日常飲食問題中，我們看到了如何對一個簡單的多變量問題進行建模，如何在 GNU MathProg 中聲明二維參數，以及如何解釋最小化問題的結果。
郵局問題中引入了 MathProg 表達式和只使用整型的決策變量。我們學習了如何分析整型問題的 glpsol 輸出結果。
最后，我們討論了整型問題的可行域的可視化問題，以及整型和相關放松限制問題的目標函數的處理方法。
本系列文章的第 3 部分即最后一部分將討論一個讓一家香水制造商實現盈利最大化的問題，并使用籃球陣容的例子來解釋二元決策變量的問題。
