备忘：MIP算法

来源：https://docs.python-mip.com/en/latest/quickstart.html#creating-models

快速开始

本章介绍了使用 Python-MIP 构建和优化模型所需的主要组件。方法及其参数的完整描述可在第 4 章中找到。

在 Python 代码中启用 Python-MIP 的第一步是添加：

from mip import *

加载后，Python-MIP 将显示其安装的版本：

Using Python-MIP package version 1.6.2

创建模型

模型类代表优化模型。下面的代码使用默认设置创建一个空的混合整数线性规划问题。

m = Model()

默认情况下，优化意义设置为最小化，选定的求解器设置为 CBC。如果 Gurobi 已安装并配置，则将使用它。您可以使用构造函数的附加参数来更改模型的客观意义或强制选择特定的求解器引擎：

m = Model(sense=MAXIMIZE, solver_name=CBC) # use GRB for Gurobi

创建模型后，您应该包括决策变量、目标函数和约束。这些任务将在下一节中讨论。

变量

使用该方法将决策变量添加到模型中add_var()。没有参数，一个带域的变量R+R+被创建并返回它的引用：

x = m.add_var()

通过使用 Python 列表初始化语法，您可以轻松创建变量向量。假设您的模型将有n 个二元决策变量（在下面的示例中 n=10），指示是否选择了 10 个项目中的每一个。下面的代码创建了 10 个二进制变量y[0]……，y[n-1]并将它们的引用存储在一个列表中。

n = 10
y = [ m.add_var(var_type=BINARY) for i in range(n) ]

其他变量类型是CONTINUOUS（默认）和INTEGER. 可以为变量指定的一些附加属性是它们的下限和上限（分别为 和 ）以及名称（属性lb）。命名变量是可选的，如果您计划以.LP 或 .MPS 文件格式保存模型（请参阅保存、加载和检查模型属性），这将特别有用。以下代码创建了一个名为的整数变量，该变量被限制在范围内ubnamezCost{−10,…,10}{−10,…,10}. 请注意，变量的引用存储在名为z.

z = m.add_var(name='zCost', var_type=INTEGER, lb=-10, ub=10)

您不需要存储变量的引用，尽管编写约束通常更容易这样做。如果您不存储这些引用，您可以在之后使用 Model 函数获取它们var_by_name()。以下代码检索命名变量的引用zCost并将其上限设置为 5：

vz = m.var_by_name('zCost')
vz.ub = 5

约束

约束是线性表达式，分别涉及变量、等于或等于、小于或等于和大于或等于的含义，以及一个常数==。约束<=>=x+y≤10x+y≤10可以很容易地包含在模型中m：

m += x + y <= 10

求和表达式可以用函数来实现xsum()。如果对于背包问题nn物品，每件都有重量wiwi，我们想包括一个约束来选择具有二进制变量的项目xixi尊重背包容量cc，则可以使用以下代码在模型中包含此约束m：

m += xsum(w[i]*x[i] for i in range(n)) <= c

在求和中条件包含变量也很容易。假设只有偶数索引项受到容量限制：

m += xsum(w[i]*x[i] for i in range(n) if i%2 == 0) <= c

最后，命名约束可能很有用。这样做很简单：在线性表达式之后包含约束的名称，用逗号分隔。下面给出一个例子：

m += xsum(w[i]*x[i] for i in range(n) if i%2 == 0) <= c, 'even_sum'

与变量一样，约束的引用可以通过它们的名称来检索。模型函数constr_by_name()对此负责：

constraint = m.constr_by_name('even_sum')

目标函数

默认情况下，使用最小化意义创建模型。以下代码将目标函数更改为∑n−1i=0cixi∑i=0n−1cixi通过设置objective我们示例模型的属性m：

m.objective = xsum(c[i]*x[i] for i in range(n))

为了指定目标是最小化还是最大化目标函数，包括了两个有用的函数：minimize()和maximize()。下面是两个使用示例：

m.objective = minimize(xsum(c[i]*x[i] for i in range(n)))

m.objective = maximize(xsum(c[i]*x[i] for i in range(n)))

您还可以通过将sense模型属性设置为MINIMIZE或来更改优化方向MAXIMIZE。

保存、加载和检查模型属性

模型方法write()和read()可用于分别保存和加载 MIP 模型。模型支持的文件格式是LP 文件格式，它更易读，更适合调试，而MPS 文件格式，推荐用于扩展兼容性，因为它是一种较旧且更广泛采用的格式。调用该write()方法时，文件扩展名（.lp 或 .mps）用于定义文件格式。因此，要将m使用 lp 文件格式的模型保存到文件 model.lp 我们可以使用：

m.write('model.lp')

同样，我们可以读取模型，从而从读取的 LP 或 MPS 文件中创建变量和约束。读取模型后，其所有属性都可用，例如约束矩阵中的变量、约束和非零的数量：

m.read('model.lp')
print('model has {} vars, {} constraints and {} nzs'.format(m.num_cols, m.num_rows, m.num_nz))

优化和查询优化结果

MIP 求解器执行分支切割 (BC) 算法，该算法将在有限时间内提供最佳解决方案。在许多情况下，这个时间可能对于您的需要来说太大了。幸运的是，即使完整的树搜索过于昂贵，通常在搜索开始时就可以获得结果。有时在处理第一个树节点时会产生一个可行的解决方案，并花费大量额外的精力来改进对偶边界，这是对最优解决方案成本的有效估计。当这个估计（最小化的下限）与找到的最佳解决方案（上限）的成本完全匹配时，搜索就结束了。对于实际应用，通常执行截断搜索。这optimize()执行配方优化的方法，接受可选的处理限制作为参数。以下代码执行分支和切割算法以求解模型m长达 300 秒。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
 
m.max_gap = 0.05
status = m.optimize(max_seconds=300)
if status == OptimizationStatus.OPTIMAL:
    print('optimal solution cost {} found'.format(m.objective_value))
elif status == OptimizationStatus.FEASIBLE:
    print('sol.cost {} found, best possible: {}'.format(m.objective_value, m.objective_bound))
elif status == OptimizationStatus.NO_SOLUTION_FOUND:
    print('no feasible solution found, lower bound is: {}'.format(m.objective_bound))
if status == OptimizationStatus.OPTIMAL or status == OptimizationStatus.FEASIBLE:
    print('solution:')
    for v in m.vars:
       if abs(v.x) > 1e-6: # only printing non-zeros
          print('{} : {}'.format(v.name, v.x))

可以使用额外的处理限制：max_nodes限制搜索树中探索节点的最大数量，并max_solutions在获得多个可行解后完成BC算法。指定结束搜索的边界应该有多紧也是明智的。模型属性max_mip_gap指定结束搜索的上限与下限的允许百分比偏差。在我们的示例中，只要上下界的距离小于或等于 5%（参见第 1 行），就可以完成搜索。

该optimize()方法返回OptimizationStatusBC搜索的状态（ ）： OPTIMAL如果搜索结束并找到最优解；FEASIBLE如果找到了可行的解决方案，但没有时间证明该解决方案是否最优； NO_SOLUTION_FOUND如果在截断搜索中没有找到解决方案；INFEASIBLE或者INT_INFEASIBLE如果模型不存在可行的解决方案； UNBOUNDED如果缺少约束或ERROR如果在优化过程中发生了一些错误。在上面的示例中，如果可行的解决方案可用（第 8 行），则打印值不为零的变量。还要注意，即使没有可行的解决方案可用，双重界限（最小化情况下的下限）也可用（第 8 行）：如果执行了截断执行，即求解器由于时间限制而停止，您可以检查这个双重界限是对检查gap 属性的解决方案质量的估计。

在树搜索过程中，通常会找到许多不同的可行解。求解器引擎将此解决方案存储在解决方案池中。以下代码打印优化旅行商问题时找到的所有路线。

for k in range(model.num_solutions):
    print('route {} with length {}'.format(k, model.objective_values[k]))
    for (i, j) in product(range(n), range(n)):
        if x[i][j].xi(k) >= 0.98:
            print('\tarc ({},{})'.format(i,j))

性能调优

MIP 求解器的树搜索算法提供了一组改进的可行解和下界。根据您的应用程序，您将对快速生成可行的解决方案更感兴趣，而不是对可能需要昂贵计算的改进下限更感兴趣，即使从长远来看，这些计算证明值得证明找到的解决方案的最优性。模型属性 emphasis提供了三种不同的设置：

  1. 默认设置：

    试图在寻找改进的可行解决方案和改进的下限之间取得平衡；

  2. 可行性：

    专注于在搜索过程的最初时刻寻找改进的可行解决方案，激活启发式；

  3. 最优性：

    激活产生改进下界的过程，即使第一个可行解决方案的产生被延迟，也专注于修剪搜索树。

将此设置更改为 1 或 2 会触发在搜索树的每个节点处处理的几种算法的激活/停用，这些算法会影响求解器的性能。即使平均而言，这些设置会如前所述改变求解器的性能，但根据您的公式，这些更改的影响可能会非常不同，通常值得在您的应用程序中使用这些不同的设置检查求解器的行为。

另一个可能值得调整的参数是cuts 属性，它控制在生成切割平面时应该花费多少计算工作量。

posted @ 2022-04-10 14:30  博二爷  阅读(656)  评论(0)  编辑  收藏  举报