线性规划、整数规划开源求解器——OR-Tools

longwilll

5901人浏览 · 2023-05-28 22:54:22

longwilll · 2023-05-28 22:54:22 发布

背景

在很多分配场景下，会有一个或多个总的待分配资源，目标是极大化回报。以折扣分配为例，公式化如下：
$\begin{aligned} \max _{x_{i, j}} & \sum_{i=1}^N \sum_{j=1}^M p_{i, j} x_{i, j} \\ \text { s.t. } & \sum_{i=1}^N \sum_{j=1}^M b_{i, j} x_{i, j} \leq B_j\\ & \sum_{j \in M} x_{i, j} \leq C_l, \forall l \in[L] \\ & x_{i, j} \in\{0,1\}, \forall i \in[N], \forall j \in[M] \end{aligned}$
具体例子
假设潜在的用户数为1000万，总共有5种券，每种券的数量为B_j，每个用户最多只能领一张券，那么上述公式则具体化为：

由于用户量极其庞大，潜在的解空间也非常大，为了加速求解，则需要专门的求解器进行求解。通过上式可以看到，目标函数和约束条件都为线性方程，我们称之为线性规划，如果求解空间属于整数，则称之为整数规划，如果求解空间是部分整数，则称之为混合整数规划。

求解器分类

根据是否收费，可以将求解器分为收费和开源两种

收费：国外的如IBM ILOG CPLEX Optimization、Gurobi等，国产的如杉数COPT、阿里MindOPT、华为OPTV等。

开源：国外的如谷歌 OR-Tools、德国的SCIP、中科院CMIP等。

详细调研可以参考: 市面上的数学规划求解器都有哪些?

开源求解器——OR-Tools

本着节约的精神，下面较为详细的介绍目前使用较广的开源求解器——OR-Tools 。这是一个用于优化的开源软件套件，专为解决世界上最棘手的车辆路线、流程、整数和线性编程以及约束编程问题而设计，具有如下特点：

 1. 它具有跨平台性。OR-Tools的核心算法是用C++进行编写的，这使其具有跨平台性。此外，它同样可以用于Python、Java或C#编译过程。

 2. 它是面向不同问题的优化工具套件。OR-Tools集合了各种先进的优化算法，它所包含的求解器主要分为约束规划、线性和整数规划、车辆路径规划以及图论算法这四个基本求解器，能够  按照优化问题的类型，提供相对应的不同类和接口。例如：对于最简单的线性规划问题，可以使用Linear Solver来解决。

它是开源且开放的。OR-Tools可以免费使用并且公开源代码。此外，OR-Tools还支持第三方求解器，可接入CPLEX等商用求解器以及SCIP等开源求解器。

安装or-tools

pip install ortools -i https://pypi.tuna.tsinghua.edu.cn/simple

MPSolver介绍

MPSolver是OR-Tools 中的一个接口，它是多个不同求解器的封装容器，可以调用收费的或者开源的求解器。根据要解决问题性质的不同，在创造求解器实例时，需要指定相应的后端(即求解器)。下面的代码都是基于python，对于其他语言可以参考代码块上方的链接。

线性规划求解器——GLOP

ortools内置的线性规划非商业求解器，一般用户求解纯线性问题，允许用户对求解器过程进行完全控制。

# Create the linear solver with the GLOP backend.
solver = pywraplp.Solver.CreateSolver('GLOP')
# 两种定义方式等价
solver = pywraplp.Solver('LinearProgrammingExample',
                             pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
if not solver:
    return
    
# 创建连续变量
infinity = solver.infinity()
# Create the variables x and y.
x = solver.NumVar(0.0, infinity, 'x')

代码示例

https://developers.google.cn/optimization/lp/lp_example?hl=zh-cn
from __future__ import print_function
from ortools.linear_solver import pywraplp
def LinearProgrammingExample():
    """Linear programming sample."""
    # Instantiate a Glop solver, naming it LinearExample.
    solver = pywraplp.Solver('LinearProgrammingExample',
                             pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
    # solver = pywraplp.Solver.CreateSolver('GLOP')

    # Create the two variables and let them take on any non-negative value.
    x = solver.NumVar(0, solver.infinity(), 'x')
    y = solver.NumVar(0, solver.infinity(), 'y')

    # Constraint 0: x + 2y <= 14.
    constraint0 = solver.Constraint(-solver.infinity(), 14)
    constraint0.SetCoefficient(x, 1)
    constraint0.SetCoefficient(y, 2)

    # Constraint 1: 3x - y >= 0.
    constraint1 = solver.Constraint(0, solver.infinity())
    constraint1.SetCoefficient(x, 3)
    constraint1.SetCoefficient(y, -1)

    # Constraint 2: x - y <= 2.
    constraint2 = solver.Constraint(-solver.infinity(), 2)
    constraint2.SetCoefficient(x, 1)
    constraint2.SetCoefficient(y, -1)

    # Objective function: 3x + 4y.
    objective = solver.Objective()
    objective.SetCoefficient(x, 3)
    objective.SetCoefficient(y, 4)
    objective.SetMaximization()

    # Solve the system.
    status = solver.Solve()
    opt_solution = 3 * x.solution_value() + 4 * y.solution_value()
    print('Number of variables =', solver.NumVariables())
    print('Number of constraints =', solver.NumConstraints())

    if status == pywraplp.Solver.OPTIMAL:
        # The value of each variable in the solution.
        print('Solution:')
        print('x = ', x.solution_value())
        print('y = ', y.solution_value())
        # The objective value of the solution.
        print('Optimal objective value =', opt_solution)
    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())

LinearProgrammingExample()
# 等价的方式
def LinearProgrammingExample():
    """Linear programming sample."""
    # Instantiate a Glop solver, naming it LinearExample.
    solver = pywraplp.Solver.CreateSolver('GLOP')
    if not solver:
        return

    # Create the two variables and let them take on any non-negative value.
    x = solver.NumVar(0, solver.infinity(), 'x')
    y = solver.NumVar(0, solver.infinity(), 'y')

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

    # 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)

    # solver.Add(x != y)

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

    # Objective function: 3x + 4y.
    solver.Maximize(3 * x + 4 * y)

    # Solve the system.
    status = solver.Solve()

    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())

混合整数规划求解器——SCIP

SCIP是目前混合整数规划（MIP）和混合整数非线性规划（MINLP）最快的非商业求解器之一。SCIP也是一个用于约束整数规划、分支定界以及分支定价的框架，主要由德国ZIB研究所开发。与大多数商业求解器不同，SCIP 允许用户对求解过程进行完全控制，并允许用户访问求解器内部的详细信息。

以下两个都是MPSolver，

Mixed-Integer Linear Programming（MIP）
 # Create the mip solver with the SCIP backend.
solver = pywraplp.Solver.CreateSolver('SCIP')
if not solver:
    return
 
# 创建整数变量
infinity = solver.infinity()
# x and y are integer non-negative variables.
x = solver.IntVar(0.0, infinity, 'x')

代码示例

https://developers.google.cn/optimization/mip/mip_example?hl=zh-cn
from ortools.linear_solver import pywraplp

def main():
    # Create the mip solver with the SCIP backend.
    solver = pywraplp.Solver.CreateSolver('SCIP')
    if not solver:
        return

    infinity = solver.infinity()
    # x and y are integer non-negative variables.
    x = solver.IntVar(0.0, infinity, 'x')
    y = solver.IntVar(0.0, infinity, 'y')

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

    # x + 7 * y <= 17.5.
    solver.Add(x + 7 * y <= 17.5)

    # x <= 3.5.
    solver.Add(x <= 3.5)

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

    # Maximize x + 10 * y.
    solver.Maximize(x + 10 * y)

    status = solver.Solve()

    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())
    print('Problem solved in %d branch-and-bound nodes' % solver.nodes())


if __name__ == '__main__':
    main()

以数组的形式快速的创建求解过程

from ortools.linear_solver import pywraplp
def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data['constraint_coeffs'] = [
        [5, 7, 9, 2, 1],
        [18, 4, -9, 10, 12],
        [4, 7, 3, 8, 5],
        [5, 13, 16, 3, -7],
    ]
    data['bounds'] = [250, 285, 211, 315]
    data['obj_coeffs'] = [7, 8, 2, 9, 6]
    data['num_vars'] = 5
    data['num_constraints'] = 4
    return data



def main():
    data = create_data_model()
    # Create the mip solver with the SCIP backend.
    solver = pywraplp.Solver.CreateSolver('SCIP')
    # 线性求解器 
    # x[j] = solver.IntVar(0, infinity, 'x[%i]' % j)
    # solver = pywraplp.Solver.CreateSolver('GLOP') 
    if not solver:
        return

    infinity = solver.infinity()
    x = {}
    for j in range(data['num_vars']):
        x[j] = solver.IntVar(0, infinity, 'x[%i]' % j)
        # 线性求解器
        # x[j] = solver.NumVar(0, infinity, 'x[%i]' % j)
    print('Number of variables =', solver.NumVariables())

    for i in range(data['num_constraints']):
        constraint = solver.RowConstraint(0, 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())
    # In Python, you can also set the constraints as follows.
    # for i in range(data['num_constraints']):
    #  constraint_expr = \
    # [data['constraint_coeffs'][i][j] * x[j] for j in range(data['num_vars'])]
    #  solver.Add(sum(constraint_expr) <= data['bounds'][i])

    objective = solver.Objective()
    for j in range(data['num_vars']):
        objective.SetCoefficient(x[j], data['obj_coeffs'][j])
    objective.SetMaximization()
    # In Python, you can also set the objective as follows.
    # obj_expr = [data['obj_coeffs'][j] * x[j] for j in range(data['num_vars'])]
    # solver.Maximize(solver.Sum(obj_expr))

    status = solver.Solve()

    if status == pywraplp.Solver.OPTIMAL:
        print('Objective value =', solver.Objective().Value())
        for j in range(data['num_vars']):
            print(x[j].name(), ' = ', x[j].solution_value())
        print()
        print('Problem solved in %f milliseconds' % solver.wall_time())
        print('Problem solved in %d iterations' % solver.iterations())
        print('Problem solved in %d branch-and-bound nodes' % solver.nodes())
    else:
        print('The problem does not have an optimal solution.')


if __name__ == '__main__':
    main()

OR-Tools 返回值及含义

状态及含义

状态	说明
OPTIMAL	已找到最佳可行解决方案。
FEASIBLE	找到了可行的解决方案，但我们不知道它是否最优。
INFEASIBLE	事实证明无法解决问题。
UNBOUNDED	解无界。
MODEL_INVALID	模型无效，如无效的系数。
NOT_SOLVED	未解决。
ABNORMAL	发生了某种异常。

整数规划求解器——CP-SAT

专门求解整数规划的求解器，功能强大，开发效率高，速度快，如果遇到约束条件为非整数的问题，则需要先将这些约束条件乘以足够大的整数，以便所有数都是整数。
详细的API介绍: https://developers.google.com/optimization/reference/python/sat/python/cp_model#newintvar

https://developers.google.cn/optimization/cp/cp_solver?hl=zh-cn
from ortools.sat.python import cp_model


class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):
    """Print intermediate solutions."""

    def __init__(self, variables):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variables = variables
        self.__solution_count = 0

    def on_solution_callback(self):
        self.__solution_count += 1
        for v in self.__variables:
            print('%s=%i' % (v, self.Value(v)), end=' ')
        print()

    def solution_count(self):
        return self.__solution_count


def SearchForAllSolutionsSampleSat():
    """Showcases calling the solver to search for all solutions."""
    # Creates the model.
    model = cp_model.CpModel()

    # Creates the variables.
    num_vals = 3
    x = model.NewIntVar(0, num_vals - 1, 'x')
    y = model.NewIntVar(0, num_vals - 1, 'y')
    z = model.NewIntVar(0, num_vals - 1, 'z')

    # Create the constraints.
    model.Add(x != y)

    # Create a solver and solve.
    solver = cp_model.CpSolver()
    solution_printer = VarArraySolutionPrinter([x, y, z])
    # Enumerate all solutions.
    solver.parameters.enumerate_all_solutions = True
    # Solve.
    status = solver.Solve(model, solution_printer)

    print('Status = %s' % solver.StatusName(status))
    print('Number of solutions found: %i' % solution_printer.solution_count())


SearchForAllSolutionsSampleSat()

https://developers.google.cn/optimization/cp/cp_example?hl=zh-cn
"""Simple solve."""
from ortools.sat.python import cp_model


def main():
    """Minimal CP-SAT example to showcase calling the solver."""
    # Creates the model.
    model = cp_model.CpModel()

    # Creates the variables.
    var_upper_bound = max(50, 45, 37)
    x = model.NewIntVar(0, var_upper_bound, 'x')
    y = model.NewIntVar(0, var_upper_bound, 'y')
    z = model.NewIntVar(0, var_upper_bound, 'z')

    # Creates the constraints.
    model.Add(2 * x + 7 * y + 3 * z <= 50)
    model.Add(3 * x - 5 * y + 7 * z <= 45)
    model.Add(5 * x + 2 * y - 6 * z <= 37)

    model.Maximize(2 * x + 2 * y + 3 * z)

    # Creates a solver and solves the model.
    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
        print(f'Maximum of objective function: {solver.ObjectiveValue()}\n')
        print(f'x = {solver.Value(x)}')
        print(f'y = {solver.Value(y)}')
        print(f'z = {solver.Value(z)}')
    else:
        print('No solution found.')

    # Statistics.
    print('\nStatistics')
    print(f'  status   : {solver.StatusName(status)}')
    print(f'  conflicts: {solver.NumConflicts()}')
    print(f'  branches : {solver.NumBranches()}')
    print(f'  wall time: {solver.WallTime()} s')


if __name__ == '__main__':
    main()

CP-SAT 返回值及含义

状态	说明
OPTIMAL	已找到最佳可行解决方案。
FEASIBLE	找到了可行的解决方案，但我们不知道它是否最优。
INFEASIBLE	事实证明无法解决问题。
MODEL_INVALID	给定的 CpModelProto 未通过验证步骤。您可以通过调用 ValidateCpModel(model_proto) 获取详细错误信息。
UNKNOWN	模型的状态未知，因为在导致求解器停止的原因（例如时间限制、内存限制或用户设置的自定义限制）之前，未找到解决方案（或者问题未证明无法解决）。