模板
教程
环境
性能
功能
管理
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 年）；参阅后面 参考资料 中的链接。
我的日常饮食要求吃的所有食物都必须属于四种基本的食物：巧克力蛋糕、冰激凌、汽水饮料和芝士蛋糕。现在有 4 种食品可以选择：巧克力松糕、巧克力冰激凌、可乐和菠萝芝士蛋糕。每个巧克力松糕价值 0.50 美元，每份巧克力冰激凌价值 0.20 美元，每瓶可乐价值 0.30 美元，每片菠萝芝士蛋糕价值 0.80 美元。每天我必须摄取至少 500 卡路里能量、6 盎司巧克力、10 盎司糖，以及 8 盎司脂肪。每种食物的单位营养含量如下表所示。现在我想以最小成本满足自己的日常营养需求。
现在我们对这个问题的相关重要信息进行一下总结：
食物的单位成本
- 巧克力松糕： $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 部分即最后一部分将讨论一个让一家香水制造商实现盈利最大化的问题，并使用篮球阵容的例子来解释二元决策变量的问题。
