1. 程式人生 > >運籌系列11:google的OR-Tools簡介

運籌系列11:google的OR-Tools簡介

1. 介紹

google的開源優化演算法包ortools,支援線性規劃、整數規劃,可以方便的求解Routing、Bin packing、Network flows、Assignment、Scheduling等問題。
官網地址為:https://developers.google.com/optimization
開原始碼地址為 https://github.com/google/or-tools
演算法包支援java、c++、c#、python。這裡使用python環境,安裝非常簡單:

pip install ortools

2. 線性規劃示例

ortools預設使用Glop(谷歌自己開發的線性規劃求解器)
這裡貼上一下官網的例子。相比其他專業的優化軟體,語法是有點冗長。

from ortools.linear_solver import pywraplp
solver = pywraplp.Solver('LinearExample',pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)

# Create the two variables and let them take on any value.
x = solver.NumVar(-solver.infinity(), solver.infinity(), 'x')
y = solver.NumVar(-solver.infinity(), solver.infinity(
), 'y') # Constraint 1: x + 2y <= 14. constraint1 = solver.Constraint(-solver.infinity(), 14) constraint1.SetCoefficient(x, 1) constraint1.SetCoefficient(y, 2) # Constraint 2: 3x - y >= 0. constraint2 = solver.Constraint(0, solver.infinity()) constraint2.SetCoefficient(x, 3) constraint2.SetCoefficient(
y, -1) # Constraint 3: x - y <= 2. constraint3 = solver.Constraint(-solver.infinity(), 2) constraint3.SetCoefficient(x, 1) constraint3.SetCoefficient(y, -1) # Objective function: 3x + 4y. objective = solver.Objective() objective.SetCoefficient(x, 3) objective.SetCoefficient(y, 4) objective.SetMaximization() # Solve the system. solver.Solve() opt_solution = 3 * x.solution_value() + 4 * y.solution_value() print('Number of variables =', solver.NumVariables()) print('Number of constraints =', solver.NumConstraints()) # 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)

輸出為:

Number of variables = 2
Number of constraints = 3
Solution:
x =  5.999999999999998
y =  3.9999999999999996
Optimal objective value = 33.99999999999999

3. 混合整數規劃

OR-Tools預設使用Coin-or branch and cut (CBC)求解整數規劃問題。此外,還支援SCIP、GLPK、Gurobi等。下面是官方例子,注意求解器和變數定義發生了變化。

from ortools.linear_solver import pywraplp
solver = pywraplp.Solver('SolveIntegerProblem',pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

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

# x + 7 * y <= 17.5
constraint1 = solver.Constraint(-solver.infinity(), 17.5)
constraint1.SetCoefficient(x, 1)
constraint1.SetCoefficient(y, 7)

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

# Maximize x + 10 * y.
objective = solver.Objective()
objective.SetCoefficient(x, 1)
objective.SetCoefficient(y, 10)
objective.SetMaximization()

"""Solve the problem and print the solution."""
result_status = solver.Solve()
# The problem has an optimal solution.
assert result_status == pywraplp.Solver.OPTIMAL

# The solution looks legit (when using solvers other than
# GLOP_LINEAR_PROGRAMMING, verifying the solution is highly recommended!).
assert solver.VerifySolution(1e-7, True)

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

# The objective value of the solution.
print('Optimal objective value = %d' % solver.Objective().Value())
print()
# The value of each variable in the solution.
variable_list = [x, y]

for variable in variable_list:
    print('%s = %d' % (variable.name(), variable.solution_value()))

輸出為

Number of variables = 2
Number of constraints = 2
Optimal objective value = 23

x = 3
y = 2

4. 約束優化

約束優化(constraint programming)和混合整數規劃的數學模型並無差異,但是它的求解方法是使用拉格朗日變換去掉約束條件,然後使用諸如牛頓法、梯度法、模擬退火、遺傳演算法、蟻群演算法等方法來快速得到一個可行解,而不是傳統MIP問題中尋找最優解
ortools可以使用原生CP求解器和CP-SAT求解器。根據官方文件說明,推薦使用CP-SAT求解器。求解status有如下幾種可能:OPTIMAL(4)、MODEL_SAT(2)、MODEL_UNSAT(3)、MODEL_INVALID(1)、UNKNOWN(0)。
另外注意CP-SAT目前不支援3.6版本,使用3.6會報錯:This program was compiled against version 3.5.1 of the Protocol Buffer runtime library, which is not compatible with the installed version (3.6.1). Contact the program author for an update.
下面舉例說明用法。假設有3個變數x, y, z,取值範圍為[0,1,2],約束條件為x≠y。

from ortools.sat.python import cp_model
# 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')

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

# Creates a solver and solves the model.
solver = cp_model.CpSolver()
status = solver.Solve(model)
print('status = %i' % status)
print('x = %i' % solver.Value(x))
print('y = %i' % solver.Value(y))
print('z = %i' % solver.Value(z))

求解結果為:

status = 2
x = 1
y = 0
z = 0

然後我們新增目標函式:

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

求解結果為:

status = 4
x = 1
y = 2
z = 2

假如想獲得中間結果,可以使用cp_model.CpSolver.SolveWithSolutionCallback(cp_model.CpModel, 回撥函式(變數list))。參考如下程式碼:

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from ortools.sat.python import cp_model
class VarArrayAndObjectiveSolutionPrinter(cp_model.CpSolverSolutionCallback):
    """Print intermediate solutions."""

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

    def OnSolutionCallback(self):
        print('Solution %i' % self.__solution_count)
        print('  objective value = %i' % self.ObjectiveValue())
        for v in self.__variables:
            print('  %s = %i' % (v, self.Value(v)), end=' ')
        print()
        self.__solution_count += 1

    def SolutionCount(self):
        return self.__solution_count

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

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

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

# Creates a solver and solves.
solver = cp_model.CpSolver()
solution_printer = VarArrayAndObjectiveSolutionPrinter([x, y, z])
status = solver.SolveWithSolutionCallback(model, solution_printer)

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

輸出結果為:

Solution 0
  objective value = 1
  x = 1   y = 0   z = 0
Solution 1
  objective value = 2
  x = 2   y = 0   z = 0
Solution 2
  objective value = 4
  x = 2   y = 1   z = 0
Solution 3
  objective value = 5
  x = 1   y = 2   z = 0
Solution 4
  objective value = 8
  x = 0   y = 1   z = 2
Solution 5
  objective value = 10
  x = 0   y = 2   z = 2
Solution 6
  objective value = 11
  x = 1   y = 2   z = 2
Status = OPTIMAL
Number of solutions found: 7