运筹优化工具ortools解读与实践-ortools求解LP、IP、MIP问题

1.前言

在之前的文章中，运筹优化工具ortools解读与实践-概览 已经介绍过线性规划、整数规划、约束规划的概念。

针对以上问题，ortools提供了多种解决方案：

• MPSolver接口：ortools提供一个名为MPSolver的接口，实现了对众多LP、MIP第三方求解器的封装，其中也包含了google自研的线性规划求解器Glop。
• CP-SAT求解器：google自研解决CP问题的求解器。为了提高求解器计算速度，CP-SAT被设计成只支持整数问题。因此，当约束中含有非整数变量时，需要通过乘以一个较大的数转变为整数问题。
• Original CP Solver：解决CP问题的早期求解器，官方建议使用CP-SAT。

有意思的是，同样是自研求解器，Glop与第三方求解器并无差别，都需要MPSolver初始化之后使用；CP-SAT却是特殊的存在，不与MPSolver一起使用，编码风格也有较大变化。

这里以一个线性规划问题引出MPSolver及Glop求解器的使用方法（整数规划和混合整数规划同样适用）。

2.套路

针对线性规划问题，求解过程分为7步：

• 1.加载一个线性求解器的wrapper
• 2.实例化求解器
• 3.创建变量，为变量设置上下限
• 4.定义约束
• 5.定义目标函数
• 6.调用求解方法
• 7.输出求解结果

在ortools官方文档中，主推google自研的GLOP求解器，从使用体验来讲，GLOP开发接口非常友好，代码开发与数学建模完全一致，执行效率也非常高。

3.案例

按照运筹优化工具ortools解读与实践-概览中的“三段式”描述解决线性规划问题的过程。

定义问题

e.g.某企业有两种不同的投资渠道，A渠道有3倍潜在收益，B渠道有4倍的潜在收益。但考虑到企业发展和风险控制，制定了一定的投资规范：
1）企业账户总金额为14W
2）每向B渠道增加一笔投资，需要留出相同数额的准备金应对风险
3）投资B渠道的金额不应超过A渠道投资额的3倍
4）A渠道投资额与B渠道投资额之差，应该2W以内

定义问题需要将后续建模中优化目标、决策变量、约束条件以自然语言的形式明确清楚。其次，作为算法人员，需要根据个人知识储备及过往经验大体明确问题类型和可能的解决方案。

在工作当中，定义问题是最重要的环节。

建模问题

优化目标：Maximize 3x + 4y
决策变量： x、y
约束条件：

• x + 2y ≤ 14
• 3x – y ≥ 0
• x – y ≤ 2
• x ≥ 0
• y ≥ 0

求解问题

#1.加载一个线性求解器的wrapper
from ortools.linear_solver import pywraplp


#2.实例化求解器，这里使用的是Google自研的GLOP线性规划求解器
solver = pywraplp.Solver.CreateSolver('GLOP')

# 3.创建变量，为变量设置上下限
x = solver.NumVar(0, solver.infinity(), 'x')
y = solver.NumVar(0, solver.infinity(), 'y')

print('Number of variables =', solver.NumVariables())

#4.定义约束
# Constraint 0: x + 2y <= 14.
solver.Add(x + 2 * y <= 14.0)

# Constraint 1: 3x - y >= 0.
solver.Add(3 * x - y >= 0.0)

# Constraint 2: x - y <= 2.
solver.Add(x - y <= 2.0)

print('Number of constraints =', solver.NumConstraints())

#5.定义目标函数
solver.Maximize(3 * x + 4 * y)

#6.调用求解方法
status = solver.Solve()

#7.输出求解结果
if status == pywraplp.Solver.OPTIMAL:
    print('Solution:')
    print('Objective value =', solver.Objective().Value())
    print('x =', x.solution_value())
    print('y =', y.solution_value())
else:
    print('The problem does not have an optimal solution.')

print('\nAdvanced usage:')
print('Problem solved in %f milliseconds' % solver.wall_time())
print('Problem solved in %d iterations' % solver.iterations())

4.MPSolver结构

除了GLOP，ortools同样支持其他求解器。

• CLP_LINEAR_PROGRAMMING or CLP
• CBC_MIXED_INTEGER_PROGRAMMING or CBC
• GLOP_LINEAR_PROGRAMMING or GLOP
• BOP_INTEGER_PROGRAMMING or BOP
• SAT_INTEGER_PROGRAMMING or SAT or CP_SAT
• SCIP_MIXED_INTEGER_PROGRAMMING or SCIP
• GUROBI_LINEAR_PROGRAMMING or GUROBI_LP
• GUROBI_MIXED_INTEGER_PROGRAMMING or GUROBI or GUROBI_MIP
• CPLEX_LINEAR_PROGRAMMING or CPLEX_LP
• CPLEX_MIXED_INTEGER_PROGRAMMING or CPLEX or CPLEX_MIP
• XPRESS_LINEAR_PROGRAMMING or XPRESS_LP
• XPRESS_MIXED_INTEGER_PROGRAMMING or XPRESS or XPRESS_MIP
• GLPK_LINEAR_PROGRAMMING or GLPK_LP
• GLPK_MIXED_INTEGER_PROGRAMMING or GLPK or GLPK_MIP

例如，将GLOP变更为SCIP或CP-SAT，同样可以得到最优解。

# solver = pywraplp.Solver.CreateSolver('GLOP')
solver = pywraplp.Solver.CreateSolver('SCIP')

有几个重要的函数，分别是：

NumVar()：创建数值型变量，当然也可以创建其他变量，例如整型变量IntVar、布尔型变量BoolVar
solver.Add()：设置约束
solver.Objective().Value()：获取目标函数结果
solution_value：获取求解结果中变量的值
SetMaximization：指定求最大值还是求最小值

了解了以上主要内容，就可以照葫芦画瓢解决一些简单的线性规划问题了！遇到其他细节可以查阅ortools源代码，了解更多方法。

5.硬编码 VS 数组存储

可以看到，我们将模型要素硬编码到代码中，这种方式可读性非常强。同样也有缺点，如果我们的要素参数并不是固定的，会随着时间的推移发生变化，此时，硬编码的形式将不利于工程化实现。

为此，ortools提供了一种数组存储模型要素的方式。

我们将第3节案例，按照数组存储模型的方式进行编码，运行结果相同。

#1.加载一个线性求解器的wrapper
from ortools.linear_solver import pywraplp

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    #约束不等式左侧部分参数
    data['constraint_coeffs'] = [
        [1, 2],
        [-3, 1],
        [1, -1],
    ]
    #约束不等式右侧部分参数
    data['bounds'] = [14.0, 0.0, 2.0]
    #目标函数部分参数
    data['obj_coeffs'] = [3,4]
    data['num_vars'] = 2
    data['num_constraints'] = 3
    return data

data = create_data_model()

#2.实例化求解器，这里使用的是Google自研的GLOP线性规划求解器
solver = pywraplp.Solver.CreateSolver('GLOP')

infinity = solver.infinity()

# 3.创建变量，为变量设置上下限
x = {}
for j in range(data['num_vars']):
    x[j] = solver.NumVar(0, infinity, 'x[%i]' % j)

print('Number of variables =', solver.NumVariables())

#4.定义约束
for i in range(data['num_constraints']):
        constraint = solver.RowConstraint(-infinity, data['bounds'][i], '')
        for j in range(data['num_vars']):
            constraint.SetCoefficient(x[j], data['constraint_coeffs'][i][j])
print('Number of constraints =', solver.NumConstraints())


#5.定义目标函数
objective = solver.Objective()
for j in range(data['num_vars']):
    objective.SetCoefficient(x[j], data['obj_coeffs'][j])
objective.SetMaximization()

#6.调用求解方法
status = solver.Solve()

#7.输出求解结果
if status == pywraplp.Solver.OPTIMAL:
    print('Solution:')
    print('Objective value =', solver.Objective().Value())
    for j in range(data['num_vars']):
        print(x[j].name(), ' = ', x[j].solution_value())
else:
    print('The problem does not have an optimal solution.')

print('\nAdvanced usage:')
print('Problem solved in %f milliseconds' % solver.wall_time())
print('Problem solved in %d iterations' % solver.iterations())

特别注意，为了约束部分代码块编码方便，需要将约束统一，我这里统一改成了≤的形式。对应的，create_data_model()方法中的编码与建模时存在差异。

原约束表达形式:

• x + 2y ≤ 14
• -3x + y ≤ 0
• x – y ≤ 2
• x ≥ 0
• y ≥ 0
  引用
  调整后的表达形式:
• x + 2y ≤ 14
• 3x – y ≥ 0
• x – y ≤ 2
• x ≥ 0
• y ≥ 0

6.结语

本篇文章主要讲解了ortools利用MPSolver及其他求解器解决LP、IP、MIP问题的两种方法。这里引出两个问题：

• 不同求解器有哪些异同，实际项目中应该如何抉择？
• CP问题如何求解，google自研的CP-SAT有何特别之处？

以上两个问题将在后续文章中探讨。