线性规划（混合整数）

本文通过举例说明线性规划（LP）和混合整数线性规划（MILP）可以解决的几个问题，说明在Sage中使用线性规划（LP）和混合整数线性规划（MILP）。给出的大多数例子都是出于图论的考虑，在没有任何特定知识的情况下应该是可以理解的。作为组合数学中的一种工具，使用线性规划就等于理解了如何通过线性约束来重新表述一个优化（或存在）问题。

这是这本书一章的翻译 Calcul mathematique avec Sage .

定义

这里我们给出线性规划的一般定义：它是由矩阵定义的 A: mathbb{{R}}^m mapsto mathbb{{R}}^n ，以及两个矢量 b,c in mathbb{{R}}^n . 求解线性规划就是寻找一个向量 x 最大化 客观的 函数和满足一组约束条件，即。

\[c^t x=max{x'text{这样}Ax'leq b}c^t x'\]

在哪里下单 u leq u' 两个向量之间意味着 u' 成对大于 u . 我们还写道：

\[text{Max:}&c^t x\text{这样：}&Axleq b\]

等价地，我们也可以说求解一个线性规划等于最大化一个定义在多面体上的线性函数（preimage或 A^{{-1}} (leq b) ). 然而，这些定义并没有告诉我们如何在组合数学中使用线性规划。在下面，我们将展示如何解决诸如背包问题、最大匹配问题和流程问题等优化问题。

混合整数线性规划

线性规划的这个定义带来了坏消息：线性规划可以在多项式时间内求解。这确实是个坏消息，因为这意味着除非我们定义指数级的线性规划，否则我们不能期望线性规划求解NP完全问题，这将是令人失望的。从一个好的方面来说，如果我们允许指定不同种类的约束，那么求解线性规划就变成了NP完全的：要求某些变量是整数而不是实值。这种线性规划实际上被称为“混合整数线性规划”（有些变量可以是整数，有些则是实数）。因此，我们可以在MILP框架中找到 wide 表现力的范围。

实用的

这个 MILP 类

这个 MILP Sage中的类代表一个MILP！它也可用于求解正则线性规划问题。它有非常少的方法，意味着定义我们的约束和变量集，然后读取解算器在计算后找到的解。也可以将用Sage定义的MILP导出到 .lp 或 .mps 文件，大多数解算器都能理解。

让我们请Sage解决以下LP：

\[text{Max:}&x+y+3z\text{例如：}&x+2yleq 4\text{}&5z-yleq 8\text{&x，y，zgeq 0\\]

为了实现这个目标，我们需要定义一个对应的 MILP 对象，以及3个变量 x ， y 和 z ：：

sage: p = MixedIntegerLinearProgram()
sage: v = p.new_variable(real=True, nonnegative=True)
sage: x, y, z = v['x'], v['y'], v['z']

接下来，我们设置目标函数

sage: p.set_objective(x + y + 3*z)

最后我们设置了限制条件

sage: p.add_constraint(x + 2*y <= 4)
sage: p.add_constraint(5*z - y <= 8)

这个 solve 方法默认返回目标函数达到的最佳值

sage: round(p.solve(), 2)
8.8

我们可以读取解算器找到的最佳分配 x ， y 和 z 通过 get_values 方法

sage: round(p.get_values(x), 2)
4.0
sage: round(p.get_values(y), 2)
0.0
sage: round(p.get_values(z), 2)
1.6

变量

在前面的例子中，我们通过 v['x'] ， v['y'] 和 v['z'] . 这就是说，更大的LP/MILP将要求我们将LP变量与许多Sage对象相关联，这些对象可以是整数、字符串，甚至是图的顶点和边。例如：

sage: x = p.new_variable(real=True, nonnegative=True)

带这个新对象 x 现在我们可以使用 x[1],...,x[15] .

sage: p.add_constraint(x[1] + x[12] - x[14] >= 8)

注意，我们不需要定义 x . 事实上， x 会像字典一样接受任何不可变的对象作为键。我们现在可以写了

sage: p.add_constraint(x["I am a valid key"] +
....:                  x[("a",pi)] <= 3)

因为任何不可变的对象都可以用作键，双索引变量 x^{{1,1}}, ..., x^{{1,15}}, x^{{2,1}}, ..., x^{{15,15}} 可被引用 x[1,1],...,x[1,15],x[2,1],...,x[15,15]

sage: p.add_constraint(x[3,2] + x[5] == 6)

类型化变量和边界

类型： 如果希望变量只采用整数或二进制值，请使用 integer=True 或 binary=True 的参数 new_variable 方法。或者，拨打 set_integer 和 set_binary 方法。

界限： 如果希望变量只取非负值，可以在调用时这样说 new_variable 有了争论 nonnegative=True . 如果要为变量设置不同的上限/下限，请添加约束或使用 set_min ， set_max 方法。

基本线性规划

背包

这个 背包 问题是：给定一组同时具有权重和 有用性 ，我们想填充一个容量受到限制的袋子，同时最大限度地发挥包中物品的效用（我们将考虑物品有用性的总和）。在本教程中，我们设置了一个限制，即袋子只能承载一定的总重量。

为了实现这一点，我们必须与每个对象相关联 o 我们收藏的 C 二进制变量 taken[o] ，当对象在包中时设置为1，否则设置为0。我们正在努力解决以下问题

\[text{Max:}&sum{oin L}text{usability}}u otimestext{taken}}u o\text{这样：}&sum{oin L}text{weight}u otimestext{taken}u oleq C\\]

使用Sage，我们将为我们的物品随机分配权重：

sage: C = 1

sage: L = ["pan", "book", "knife", "gourd", "flashlight"]

sage: L.extend(["random_stuff_" + str(i) for i in range(20)])

sage: weight = {}
sage: usefulness = {}

sage: set_random_seed(685474)
sage: for o in L:
....:     weight[o] = random()
....:     usefulness[o] = random()

我们现在可以定义MILP本身了

sage: p = MixedIntegerLinearProgram()
sage: taken = p.new_variable(binary=True)

sage: p.add_constraint(p.sum(weight[o] * taken[o] for o in L) <= C)

sage: p.set_objective(p.sum(usefulness[o] * taken[o] for o in L))

sage: p.solve() # abs tol 1e-6
3.1502766806530307
sage: taken = p.get_values(taken)

找到的解决办法（当然）是可以接受的

sage: sum(weight[o] * taken[o] for o in L)  # abs tol 1e-6
0.6964959796619171

我们要不要拿个手电筒？

sage: taken["flashlight"] # abs tol 1e-6
1.0

明智的建议。基于纯粹的随机考虑。

匹配

给定一张图 G 匹配是一组成对不相交的边。空集是一个微不足道的匹配。所以我们把注意力集中在最大匹配上：我们想在图中找到一个基数最大的匹配。计算图的最大匹配是一个多项式问题，这是Edmonds的著名结果。Edmonds的算法是基于局部改进的，并且证明了给定的匹配在不能改进的情况下是最大的。这个算法并不是图论所能提供的最难实现的，尽管这个问题可以用一个非常简单的MILP来建模。

为此，我们需要——如前所述——将一个二进制变量与我们的每一个对象相关联：图的边（值为1表示对应的边包含在最大匹配中）。我们对边缘的限制是它们是不相交的，这就足够了， x 和 y 有两面性 m_x, m_y 它们的相关变量，不平等 m_x + m_y leq 1 是满意的，因为我们确信他们两个不可能都属于匹配。因此，我们能够写出我们想要的MILP。然而，当且仅当两条边有一个共同的端点时，可以很容易地减少不等式的数量 v . 我们可以要求最多一个边缘事件到 v 在匹配的内部，这是一个线性约束。我们将解决：

\[文本{Max:}&sum{e在e（G）}m_e\text{例如：}&for all v，sum{substack{e在e（G）\vsim e}}m_eleq 1\]

让我们写下这个弥勒的圣典：

sage: g = graphs.PetersenGraph()
sage: p = MixedIntegerLinearProgram()
sage: matching = p.new_variable(binary=True)

sage: p.set_objective(p.sum(matching[e] for e in g.edges(labels=False)))

sage: for v in g:
....:     p.add_constraint(p.sum(matching[e]
....:         for e in g.edges_incident(v, labels=False)) <= 1)

sage: p.solve()
5.0

sage: matching = p.get_values(matching)
sage: [e for e, b in matching.items() if b == 1]  # not tested
[(0, 1), (6, 9), (2, 7), (3, 4), (5, 8)]

流

图论中的另一个基本算法：最大流！它由一个有向图和两个顶点组成 s, t ，在发送最大 flow 从 s 到 t 利用 G ，每个人都有一个最大的容量。

thematic_tutorials/media/lp_flot1.png

这个问题的定义几乎就是它的线性规划公式。我们正在寻找与每条边相关的实际值，该值将表示通过这些边的流的强度，在两种类型的约束下：

  • 到达顶点的流量（不同于 s 或 t )等于流出的流量。

  • 通过一条边的流量受这条边的容量的限制。

这就是说，我们必须使流出的流量最大化 s ：所有的一切都将在 t ，因为其他顶点发送的和接收的一样多。我们可以用下面的线性规划来模拟流动问题

\[文本{Max:}&sum{sv在G}f{sv}\text{这样：}&forall ving，{substack{vneq s\vneq t}，sum{vuin G}f{vu}-sum{uvin G}f{uv}=0\&for all uvin G，f{uv}leq 1\\]

我们将在Chvatal图的一个方向上求解流问题，其中所有边的容量为1:：

sage: g = graphs.ChvatalGraph()
sage: g = g.minimum_outdegree_orientation()

sage: p = MixedIntegerLinearProgram()
sage: f = p.new_variable(real=True, nonnegative=True)
sage: s, t = 0, 2

sage: for v in g:
....:     if v != s and v != t:
....:         p.add_constraint(
....:             p.sum(f[v,u] for u in g.neighbors_out(v))
....:             - p.sum(f[u,v] for u in g.neighbors_in(v)) == 0)

sage: for e in g.edges(labels=False):
....:     p.add_constraint(f[e] <= 1)

sage: p.set_objective(p.sum(f[s,u] for u in g.neighbors_out(s)))

sage: p.solve()  # rel tol 2e-11
2.0
thematic_tutorials/media/lp_flot2.png

解算器（后端）

Sage通过调用特定的库来解决线性程序。当前支持以下库：

  • CBC: 一个来自 COIN-OR ，在Eclipse公共许可证（EPL）下提供，EPL是一个开源许可证，但与GPL不兼容。可以使用shell命令安装CBC和Sage CBC后端：

    $ sage -i -c sage_numerical_backends_coin

  • CPLEX: 专利，但通过IBM的学术计划，研究人员和学生可以免费获得。因为 :trac:`27790` ，仅支持12.8及更高版本。

    根据网站上的说明安装CPLEX，其中包括获取许可证密钥。

    然后找到ilogcplexstudio安装的安装目录，其中包含子目录 cplex ， doc ， opl 等。设置环境变量 CPLEX_HOME 例如，使用以下shell命令（在macOS上）：：

    $ export CPLEX_HOME=/Applications/CPLEX_Studio1210

    或者（在Linux上）：

    $ export CPLEX_HOME=/opt/ibm/ILOG/CPLEX_Studio1210

    现在验证您将在子目录中找到的CPLEX二进制文件 cplex/bin/ARCH-OS 正确启动，例如：

    $ $CPLEX_HOME/cplex/bin/x86-64_osx/cplex
Welcome to IBM(R) ILOG(R) CPLEX(R) Interactive Optimizer...

    此环境变量只需在以下步骤中设置：使用shell命令安装Sage CPLEX后端：

    $ sage -i -c sage_numerical_backends_cplex

  • CVXOPT: Python凸优化软件的LP解算器使用内部点法，始终安装在Sage中。

    根据GPL授权。

  • GLPK: A solver from GNU

    根据GPLv3授权。此解算器始终作为默认解算器安装在Sage中。

  • Gurobi: 专利，但研究人员和学生可以通过古洛比的学术项目免费获得。

    根据网站上的说明安装Gurobi，其中包括获取许可证密钥。安装程序应该使交互式Gurobi shell gurobi.sh 在您的 PATH . 通过键入shell命令来验证这一点 gurobi.sh ：：

    $ gurobi.sh
Python 3.7.4 (default, Aug 27 2019, 11:27:39)
...
Gurobi Interactive Shell (mac64), Version 9.0.0
Copyright (c) 2019, Gurobi Optimization, LLC
Type "help()" for help
gurobi>

    如果这不起作用，请调整 PATH 或者创建符号链接以便 gurobi.sh 被发现。

    现在使用shell命令安装Sage Gurobi后端：

    $ sage -i -c sage_numerical_backends_gurobi

  • PPL: 一个来自布森的解决方案。

    此解算器提供精确（任意精度）计算，始终安装在Sage中。

    根据GPLv3授权。