GNU 线性编程工具包

GNUGNU 线性编程工具包 线性编程工具包 GNUGNU LinearLinear ProgrammingProgramming KitKit 第 第 1 1 部分部分 线性优化简介线性优化简介 为复杂数学问题寻找最佳解决方案 级别 中级 Rodrigo Ceron rceron 资深软件工程师 IBM 2006 年 9 月 04 日 GNU Linear Programming Kit 对于解决具有多种约束的数学问题来说是一个功能非常强大的工具 本文简要 介绍了如何使用 GLPK glpsol 客户机工具 和 GNU MathProg 语言来解决 Giapetto 的 Woodcarving 公司 一家玩具制造商 的作业优化问题 简介 线性编程是一个用来解决优化问题的工具 在 1947 年 George Dantzig 开发了一种效率方法 simplex 算法 来解决线性编程的问题 由于 simplex 算法的出现 线性编程已经在工业界 银行界 教育界 林业 石油行业以及 运输业界中广泛地用来解决优化问题 在对财富 500 强公司的调查中 85 的被调查者都说他们已经使用了线性编程 引自 Operations Research Applications and Algorithms 4th Edition Wayne L Winston 著 Thomson 2004 请参看下面 参考资料 中的链接 有很多工具都可以用来解决线性编程的问题 尽管有些专用工具都非常出名 但是开放源码社区中的很多人可能并不了解免费的 GLPK 工具 本系列文章一共 3 篇 将逐渐介绍 GLPK 的功能和用法 本文是其中的第 1 篇 在本文中 我们将对 GLPK 简要进行介绍 然后展示并应用 GLPK 中的 文档选项 未显示需 要 JavaScrip t 的文档选 项 GNU MathProg Language 如果我们仅仅从运筹学理论入手 希望学习进行建模和解决线性编程的问题 那么本文就是一个很好的指南 GNU Linear Programming Kit GNU Linear Programming Kit GLPK 是一个使用了众所周知的运筹学算法来解决线性问题的程序库 这些程序实现了simplex 算法 branch and bound 算法 primal dual interior point 算法以及很多其他算法 请查看 GLPK 下载包中所包含的 GLPK 手册以便了解更多的内容 要下载 GLPK 请参看 参考资料 一节中给出的 gnu org 上面 GLPK 页面的链接 GLPK 不是一个程序 它无法运行 也没有 main 函数 实际上 客户机需要将问题数据通过 GLPK API 提供给算法程序 并接收结果 GLPK 有一个默认的客户机 即 glpsol 程序 它可以与这个 API 进行交互 通常诸如 glpsol 之类的程序都被称为 solver 而不是客户机 因此从现在开始我们就会看到这个术语 GNU MathProg 建模语言 GNU MathProg 建模语言非常简单 但却可以很好地声明线性问题 通常 问题的声明包括 问题决策变量 目标函数 注意目标 objective 是一个名词 而不是一个形容词 这个名字是运筹学理论中的标准术语 问题约束 问题参数 数据 让我们从一个简单的两变量的例子入手 Giapetto 的 Woodcarving 公司 Giapetto 的 Woodcarving 公司 这个问题引自于 Operations Research Giapetto 的 Woodcarving 公司生产两种木头制作的玩具 士兵和火车 一个士兵的销售价格为 27 回页首回页首 回页首回页首 美元 需要耗费价值 10 美元的原料 制造每个士兵需要耗费 Giapetto 的可变人力成本和间接成本一共 14 美元 一辆火车的销售价格为 21 美元 需要耗费价值 9 美元的原料 制造每辆火车需要耗费 Giapetto 的可变人力成本和间接成本一共 10 美元 这家木头士兵和火车的制造商需要两类熟练工人 木工和修整工 一个 士兵需要 2 小时的修整工劳动和 1 小时的木工劳动 一辆火车需要 1 小时的修整工劳动和 1 小时的木工劳动 每周 Giapetto 可以获得所有必需的原料 但是只能提供 100 个修整工时和 80 个木工工时 市场对于火车的需求是无限的 但是每周最多可以销售 40 个士兵 Giapetto 希望能够使每周的收益 收入 成本 最大化 下面我们对这个问题的重要信息和假设小结一下 1 有两种木制玩具 士兵和火车 2 一个士兵的销售价格为 27 美元 需要耗费价值 10 美元的原料 另外需要耗费可变人力成本和间接成本一共 14 美元 3 一辆火车的销售价格为 21 美元 需要耗费价值 9 美元的原料 另外需要耗费可变人力成本和间接成本一共 10 美元 4 一个士兵需要 2 小时的修整工劳动和 1 小时的木工劳动 5 一辆火车需要 1 小时的修整工劳动和 1 小时的木工劳动 6 每周最多可以获得 100 个修整工时和 80 个木工工时 7 每周对于火车的需求是无限的 但是最多可以销售 40 个士兵 我们的目标是确定每周生产士兵和火车的数量 从而可以使每周的收益最大化 建模 要对一个线性问题进行建模 必须首先确定决策变量 这是因为这些变量在每 个 simplex 算法循环中都会变化 并决定了目标函数的值 这样才会产生优化解决方案 在 Giapetto 的商店中 目标函数是收益 这是一个每周生产的士兵和火车的函数 因此 这个问题中的两个决策变量是 x1 每周生产的士兵数量 x2 每周生产的火车数量 确定了决策变量之后 这个问题的目标函数就简化为每个玩具的收入减去成本 这是 x1和 x2的一个函数 回页首回页首 请注意 收益线性依赖于 x1和 x2 这是一个线性问题 乍一看 收益可以通过简单地增大 x1和 x2来实现最大化 然而 如果生活是如此简单 那就让我们都开始生产火车和 士兵 并迁居到 Caribbean 好了 不幸的是 有很多约束会对可以选择的决策变量造成影响 否则这个模 型很可能就是错的 回想一下我们对这个问题的 假设 前三条确定了决策变量和目标函数 第 4 个假设和第 6 个假设说完成士兵的制造需要耗费一定的木工和修整工工时 此处的约束是 Giapetto 并没有无限的木工和修整工工时 这就是约束 下面我们来分析并澄清这个问 题 一个士兵需要 2 小时的修整工时劳动 Giapetto 每周最多有 100 小时的修整工劳动力 因此每周就不能生产超过 50 个士兵 类似地 木工工时的约束也使它不能每周生产超过 80 个士兵 注意第一个约束比第二个约束更为严格 第一个约束实际上是第二个 约束的一个子集 因此第二个约束实际上是冗余的 上面这段描述展示了如何对优化问题进行建模 但这只是不完全的一个分析 因为所有必需的变量都还没有充分考虑 这尚不是 Giapetto 问题的彻底解决方案 那么我们应该如何解决这个问题呢 我们首先通过分析约束因素来发现所有的约束条件 首先 到底是什么约束的 修整工时 由于士兵和火车都需要修整时间 因此它们都需要考虑在 内 假设要制造 10 个士兵和 20 辆火车 这所需要的修整工时是 10 乘以 2 小时 用来生产士兵 加上 20 乘以 1 小时 用来生产火车 总共是 40 小时的修整劳动力 按照决策变量的术语来说 通用的约束为 有很多 x1 x2 对能够满足这个不等式 因此它无法确定 Giapetto 商店的最佳组合 现在我们已经知道了修整工时的约束 同样也可以得出木工工时的约束 非常不错 这个问题只剩下一个约束了 记得每周对士兵的需求么 根据问题 描述 每周最多可以生产 40 个士兵 对于火车的需求是无限的 因此对它就没有什么约束了 这个模型就完成了 它包括以下等式 注意最后一个约束 它确保决策变量的值都是正数 问题并没有显式地说明这 一点 但它却非常重要 也很显然 现在 GLPK 可以解析这个模型了 由于 GLPK 很擅长解决线性优化问题 一点理论知识 下面我们来了解以下问题的解析空间 有了这两个决策变量 它就有了两个维 度 图图 1 1 GiapettoGiapetto 的极大域的极大域 回页首回页首 x1 x2 超出第一象限 其中所有值都为正数 的结果都已经被舍弃了 然而 我们需 要注意这个解析空间仍然是无限的 这依然可能 成为我希望移居 Caribbean 的一种情况 正如我们给出的约束一样 这个无限的解析空间达到了边界 有了上面的 不等式 6 结果就更加有趣了 图图 2 2 GiapettoGiapetto 考虑修整约束时的解析域考虑修整约束时的解析域 这个解析空间包含了 x1 x2 在第一象限的所有可能值 它们可以满足修整工时约束 在加上 不等式 7 之后 结果就收缩了 图图 3 3 GiapettoGiapetto 考虑修整约束和木工约束时的解析域考虑修整约束和木工约束时的解析域 注意 这个解析空间更小 这意味着其中只有更少的 x1 x2 值 在应用不等式 8 之后 结果还会进一步变小 图图 4 4 GiapettoGiapetto 的可行域的可行域 这样解析空间更小了 这个可以满足所有约束的解析空间称为可行域 feasib le region 图 4 给出了 Giapetto 商店的可行域 这个区域中任何 x1 x2 对都是这个问题的一个可能解决方案 现在的问题是 哪个值能够使得 Giapetto 的收益最大化呢 使用 GLPK 来解析模型 GLPK 对于解析这个问题来说是一个很好的工具 Giapetto 问题中的数学公式需要使用 GNU MathProg 语言进行编写 需要声明的关键内容有 决策变量 目标函数 约束 问题数据集 下面的代码显示了如何使用 MathProg 来解决 Giapetto 回页首回页首 问题 代码中的行号都不是代码本身的一部分 添加这些行号只是为了能够方 便地对代码进行引用 清单清单 1 1 GiapettoGiapetto 问题的第一个解决方案 问题的第一个解决方案 giapetto solgiapetto sol 1 2 Giapetto s problem 3 4 This finds the optimal solution for maximizing Giapetto s profit 5 6 7 Decision variables 8 var x1 0 soldier 9 var x2 0 train 10 11 Objective function 12 maximize z 3 x1 2 x2 13 14 Constraints 15 s t Finishing 2 x1 x2 100 16 s t Carpentry x1 x2 80 17 s t Demand x1 0约束 正如我们在第 8 行和第 9 行中看到的一样 GNU MathProg 中的每条语句都必须以分号 结束 记住 x1表示 soldier的数量 x2表示 train的数量 这些变量应该分别称为 soldiers和 trains 但是这可能会把读者中的数学界人士搞得混乱不堪 第 12 行声明了目标函数 线性问题可以是最大化问题 也可以是最小化问题 记住 Giapetto 的数学模型是一个最大化问题 因此我们应该使用 maximize作为关键字 而不是相反的 minimize关键字 目标函数称为 z 它等于 3x1 2x2 请注意 冒号 字符分隔开了目标函数及其定义 星号 字符表示乘法 类似地 加号 减号 和斜线 字符分别表示加法 减法和除法 第 15 16 17 行定义了约束 尽管 s t 对于声明一个约束来说并不需要在行首 但是这样可以提高代码的可读 性 这三个 Giapetto 约束分别被标记成 Finishing Carpentry 和 Demand 每个都是作为一个数学模型来声明的 符号 表示不等式 不要忘记每个声明末尾的 每个 GNU MathProg 文件必须要以 end 结束 正如我们在 19 行看到的一样 现在 glpsol 可以使用这个文件作为输入 但是请等一下 这个问题的数据部分在什么地方 呢 噢 这个问题是如此简单 问题数据都作为声明中决策变量的系数直接包 含在了目标函数和约束声明中 举例来说 在目标函数中 系数 3和 1就是问题数据集的一部分 当我使用数据集重写这个问题时 大家就可以非 常清楚这到底是如何工作的 现在 我们不用关心这个问题 使用 mod作为 MathProg 输入文件的扩展名并将答案重定向到一个扩展名为 sol的文件中是一个好习惯 不过这并不是必需的要求 您可以使用任何自 己喜欢的文件名和扩展名 这个例子中 Giapetto 的 MathProg 文件是 giapetto mod 输出结果在 giapetto sol中 现在 请在自己喜欢的终端中运行 glpsol glpsol m giapetto mod o giapetto sol 这个命令使用了两个 glpsol 选项 m选项告诉 glpsol 输入是使用 GNU MathProg 编写的 o选项告诉 glpsol 将自己的输出结果发送到 giapetto sol中 解答报告现在就在 giapetto sol中了 但是有关 GLPK 所消耗的时间和内存信息会显示在标准输出上 清单清单 2 2 glpsolglpsol 的输出的输出 ceron curly glpsol m giapetto mod o giapetto sol Reading model section from giapetto real mod 19 lines were read Generating z Generating Finishing Generating Carpentry Generating Demand Model has been successfully generated lpx simplex original LP has 4 rows 2 columns 7 non zeros lpx simplex presolved LP has 2 rows 2 columns 4 non zeros lpx adv basis size of triangular part 2 0 objval 0 000000000e 00 infeas 0 000000000e 00 0 2 objval 1 400000000e 02 infeas 0 000000000e 00 0 OPTIMAL SOLUTION FOUND Time used 0 0 secs Memory used 0 1M 151326 bytes lpx print sol writing LP problem solution to giapetto sol 所生成的报告显示 glpsol 会读取这个模型 调用一个 GLPK API 函数来生成目标函数 然后调用另外一个 API 函数来生成约束 在生成模型之后 glpsol 简要地解释了 GLPK 在内部是如何对这个问题进行处理的 最后是有关解答和 GLPK 为解答这个问题所消耗的资源的信息 这个解答被明确说明是优化的 这非常棒 但是决策变量的实际优化值是什么呢 它们都保存在 giapetto sol文件中 清单清单 3 3 GiapettoGiapetto 问题的解答 问题的解答 giapetto solgiapetto sol Problem giapetto Rows 4 Columns 2 Non zeros 7 Status OPTIMAL Objective z 180 MAXimum No Row name St Activity Lower bound Upper bound Marginal 1 z B 180 2 Finishing NU 100 100 1 3 Carpentry NU 80 80 1 4 Demand B 20 40 No Column name St Activity Lower bound Upper bound Marginal 1 x1 B 20 0 2 x2 B 60 0 Karush Kuhn Tucker optimality conditions KKT PE max abs err 0 00e 00 on row 0 max rel err 0 00e 00 on row 0 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 0 00e 00 on column 0 max rel err 0 00e 00 on column 0 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 解答被划分成 4 部分 有关问题和目标函数优化值的信息 有关目标函数状态和约束的确切信息 有关决策变量的优化值的确切信息 有关优化条件的信息 如果有的话 对于这个特定的问题来说 我们可以看到解决方案是 OPTIMAL的 Giapetto 每周最大收益是 180 美元 Finishing 约束的状态是 NU St 列 NU表示这个约束的上界有一个非基本变量 如果我们了解一些基本的 运筹学理论 就可以构建 simplex 场景并将之提取出来 如果不了解运筹学理论 下面是一个实际的解释 不论何时约束达到自己的上界或下界时 我们就称之为是一个边界约束 boun ded constraint 边界约束阻碍了目标函数达到更为优化的值 例如我们可以认 为它是一个音量调节旋钮 现在它已经无法再进行调节了 当这种情况发生时 glpsol 就会将约束状态显示为 NU或 NL 分别对应上界和下界的情况 它还会显示边界 marginal 值 也称 为假定价格 shadow price 边界值是约束如果可以放松一个单位 如果音量调节旋钮可以再调节一点点 目标函数可以改进的值 注意这种改进依赖于我们的目标是要对目标函数进行 最大化还 是最小化 举例来说 Giapetto 问题所寻求的就是最大化 边界值为 1 就表示如果我们可以增加一个小时的修整工时 目标函数就可以增大 1 我们知道这是要多 一个小时 而不是少一个小时 这是因为修整工时约束是一个上界 木工和士兵的需求约束是类似的 对于木工约束来说 注意它也有一个上界 因此 这个约束中放松一个单位 增加一个小时 可以使目标函数的优化值增 加边界值 1 成为 181 然而 士兵需求是没有边界的 因此其状态是 B 这个约束中的放松不会改变目标函数的优化值 一次只尝试放松每个边界约束一个单位 解决修改后的问题 看一下目标函数 的优化值发生了什么变化 还要检查修改无边界约束的值不会对解答造成任何 变化 这正是我们期望的 最后 glpsol 的报告给出了决策变量的值 x1 20和 x2 60 这就告诉 Giapetto 它应该生产 20 个士兵和 60 辆火车才能实现每周收益的最大化 Giapetto 的问题很小 我们可能会纳闷 在有更多决策变量和约束的问题中 我们只能 分别逐一声明每个变量和每个约束吗 如果希望调节问题中的数据 例如士兵 的销售价格 应该怎样做呢 我们只能逐一修改这些值出现的地方吗 下一节 将讨论这个问题 在 Giapetto 问题中使用模型和数据段 MathProg 模型通常都有一个模型段和一个数据段 有时可以在两个不同的文件中 这样 glpsol 就可以简单地使用不同的数据集来解析某个模型 从而找到对这些新数据应该 采用哪种解决方案 下面的列表以更优雅的方式说明了 Giapetto 的问题 清单清单 4 4 使用一个模型段和一个数据段来分析使用一个模型段和一个数据段来分析 GiapettoGiapetto 问题 问题 giapetto2 modgiapetto2 mod 1 2 Giapetto s problem 3 4 This finds the optimal solution for maximizing Giapetto s profit 5 6 7 Set of toys 8 set TOY 9 10 Parameters 11 param Finishing hours i in TOY 12 param Carpentry hours i in TOY 13 param Demand toys i in TOY 14 param Profit toys i in TOY 15 16 Decision variables 17 var x i in TOY 0 18 19 Objective function 20 maximize z sum i in TOY Profit toys i x i 21 22 Constraints 23 s t Fin hours sum i in TOY Finishing hours i x i 100 24 s t Carp hours sum i in TOY Carpentry hours i x i 80 25 s t Dem i in TOY x i Demand toys i 26 27 回页首回页首 28 data 29 data section 30 31 set TOY soldier train 32 33 param Finishing hours 34 soldier 2 35 train 1 36 37 param Carpentry hours 38 soldier 1 39 train 1 40 41 param Demand toys 42 soldier 40 43 train 6 02E 23 44 45 param Profit toys 46 soldier 3 47 train 2 48 49 end 我们并没有使用两个单独的文件 而是在一个具有模型段 第 1 行 到第 27 行 和一个数据段 第 28 行到第 49 行 的文件中声明了这个问题 第 8 行定义了一个 SET SET是一个元素的值域 例如 如果我声明数学变量 xi for all i in 1 2 3 4 就说明 x是一个包含 4 个元素的数组 因此我们就可以使用 x1 x2 x3 x4了 在这个例子中 1 2 3 4 就是这个集合 因此 在 Giapetto 问题中 有一个集合 TOY 这个集合的实际值是什么呢 在这个文件的数据段中 请查看第 31 行 它定义了 TOY集合包含 soldier和 train 哦 这个集合不是一个数字类型的范围 那么这怎样实现呢 GLPK 是以一种非常有趣的方法来处理这个问题的 稍后我们就会看到如何使用这种 方法 现在 请习惯数据段中 SET声明所使用的语法 setlabel value1 value2 valueN 第 11 12 和 13 行定义了这个问题的参数 一共有三个参数 Finishing hours Carpentry hours和 Demand toys 这些参数构成了这个问题的数据矩阵 并用来计算约束 正如 我们稍后会看到的一样 我们以 Finishing hours参数为例来看一下这个问题 它是在 TOY集合上定义的 因此 TOY集合中的每种玩具都有一个 Finishing hours值 记住每个士兵需要 2 小时的修整工作 每辆火车需要 1 小时的修整时间 现在请看一下第 33 34 和 35 行的内容 这是对每种玩具的修整工时的定义 实际上 这些行的内容声明 Finishing hours soldier 2 Finishing hours train 1 因此 Finishing hours就是一个 1 行 2 列的矩阵 木工工时和需求参数都是类似地进行声明的 注意对于火车的需求是无限的 因此在 43 行上我们使用了一个很大的值作为上界 这个值看起来够大么 第 17 行声明了一个变量 x 对于 TOY中的每个 i 即 x soldier 和 x train 我们都要限制它们大于或等于 0 有了这些集合之后 声明变量就变得非常简单了 不是么 第 20 行声明了对象函数是 z的最大值 这是每种玩具的总体收益 一共有两种玩具 火车和士兵 例 如 对于士兵来说 收益是士兵个数乘以每个士兵的收益 第 23 24 和 25 行上的约束也是类似地进行声明的 以修整工时约束为例 它是每种玩具的修 整工时乘以所生产的这种玩具的数量之积的总和 对于这两种玩具 火车和士 兵 来说 这必须小于或等于 100 类似地 木工工时综合也应该小于或等于 80 需求约束与上面两个约束有所不同 因为它是为每种玩具而定义的 而不是为 所有玩具类型总和进行定义的 因此 我们需要两个约束 一个用于火车 另 外一个用于士兵 正如我们在第 25 行中看到的一样 注意索引变量 i in TOY 都是在 之前出现的 这告诉 GLPK 为 TOY 中的每种玩具类型都创建一个约束 每个约束的等式都在 之后出现 在本例中 GLPK 会创建 Dem soldier x soldier Demand soldier Dem train x train 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 行 的语法如下 清单清单 2 2 两维参数声明的语法两维参数声明的语法 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集合 清单清单 3 3 这个问题中这个问题中 glpsolglpsol 的输出的输出 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 名员工在工作 类似地 回页首回页首 当然 不要忘记了符号约束 邮局问题的 GNU MathProg 解决方案 注意 清单 4 中的行号仅仅是为了参考方便而给出的 清单清单 4 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时 对于其他约束来说全部类似 逻辑比较符也可以在这些表达式中使用 清单清单 5 5 这个问题在这个问题在 glpsolglpsol 中的解答中的解答 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 0 integer 这相当简单 glpsol的输出对整型情况的结果稍有不同 清单清单 6 6 glpsolglpsol 对整型约束邮局问题的输出结果对整型约束邮局问题的输出结果 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 successfully 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 已经对放松限制下的问题 不要求变量是整数 计算了优化解决方案 清单清单 7 7 邮局问题的整型解决方案邮局问题的整型解决方案 Problem post Rows 8 Columns 7 7 integer 0 binary Non zeros 42 Status INTEGER OPTIMAL Objective z 23 MINimum 22 33333333 LP N