Constraint Optimization(借助ortools)

注：中文非直接翻译英文，而是理解加工后的笔记，记录英文仅为学其专业表述。

目录

• 概述
• SP-SAT Solver
  • 加目标函数，寻找最优解
  • 加回调函数，展示所有可行解
  • 展示 intermediate solutions
• Original CP Solver
• Cryptarithmetic Puzzles
  • 建模
• The N-queens Problem
  • 传播和回溯
• Setting solver limits
  • 时间限制
  • 找到指定数量的解后停止搜索

概述

Constraint optimization, or constraint programming(CP)，约束优化(规划)，用于在一个非常大的候选集合中找到可行解，其中的问题可以用任意约束来建模。

CP基于可行性(寻找可行解)而不是优化(寻找最优解)，它更关注约束和变量，而不是目标函数。

SP-SAT Solver

CP-SAT求解器技术上优于传统CP求解器。

The CP-SAT solver is technologically superior to the original CP solver and should be preferred in almost all situations.

为了增加运算速度，CP求解器处理的都是整数。

如果有非整数项约束，可以先将其乘一个整数，使其变成整数项。

If you begin with a problem that has constraints with non-integer terms, you need to first multiply those constraints by a sufficiently large integer so that all terms are integers

来看一个简单例子：寻找可行解。

• 变量x，y，z，每个只能取值0，1，2。
  • eg: x = model.NewIntVar(0, num_vals - 1, 'x')
• 约束条件：x ≠ y
  • eg: model.Add(x != y)

核心步骤：

• 声明模型
• 创建变量
• 创建约束条件
• 调用求解器
• 展示结果

from ortools.sat.python import cp_model

def SimpleSatProgram():
    """Minimal CP-SAT example to showcase calling the solver."""
    # 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)

    if status == cp_model.FEASIBLE:
        print('x = %i' % solver.Value(x))
        print('y = %i' % solver.Value(y))
        print('z = %i' % solver.Value(z))

SimpleSatProgram()

运行得

x = 1
y = 0
z = 0

status = solver.Solve(model)返回值状态含义:

• OPTIMAL 找到了最优解。
• FEASIBLE 找到了一个可行解，不过不一定是最优解。
• INFEASIBLE 无可行解。
• MODEL_INVALID 给定的CpModelProto没有通过验证步骤。可以通过调用ValidateCpModel(model_proto)获得详细的错误。
• UNKNOWN 由于达到了搜索限制，模型的状态未知。

加目标函数，寻找最优解

假设目标函数是求x + 2y + 3z的最大值，在求解器前加

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

并替换展示条件

if status == cp_model.OPTIMAL:

运行得

x = 1
y = 2
z = 2

加回调函数，展示所有可行解

去掉目标函数，加上打印回调函数

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])
    status = solver.SearchForAllSolutions(model, solution_printer)

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

SearchForAllSolutionsSampleSat()

运行得

x=0 y=1 z=0
x=1 y=2 z=0
x=1 y=2 z=1
x=1 y=2 z=2
x=1 y=0 z=2
x=1 y=0 z=1
x=2 y=0 z=1
x=2 y=1 z=1
x=2 y=1 z=2
x=2 y=0 z=2
x=0 y=1 z=2
x=0 y=1 z=1
x=0 y=2 z=1
x=0 y=2 z=2
x=1 y=0 z=0
x=2 y=0 z=0
x=2 y=1 z=0
x=0 y=2 z=0
Status = OPTIMAL
Number of solutions found: 18

展示 intermediate solutions

与展示所有可行解的异同主要在:

• 加 model.Maximize(x + 2 * y + 3 * z)
• 换 status = solver.SolveWithSolutionCallback(model, solution_printer)

输出结果为:

x=0 y=1 z=0
x=0 y=2 z=0
x=0 y=2 z=1
x=0 y=2 z=2
x=1 y=2 z=2
Status = OPTIMAL
Number of solutions found: 5

与官网结果不一致。https://developers.google.cn/optimization/cp/cp_solver?hl=es-419#run_intermediate_sol

应该与python和or-tools版本有关。

Original CP Solver

一些特殊的小问题场景，传统CP求解器性能优于CP-SAT。

区别于CP-SAT的包，传统CP求解器用的是from ortools.constraint_solver import pywrapcp

打印方式也有区别，传统CP求解器打印方式是用while solver.NextSolution():来遍历打印。

构造器下例子中只用了一个phase，复杂问题可以用多个phase，以便于在不同阶段用不同的搜索策略。Phase()方法中的三个参数:

• vars 包含了该问题的变量，下例是 [x, y, z]
• IntVarStrategy 选择下一个unbound变量的规则。默认值 CHOOSE_FIRST_UNBOUND，表示每个步骤中，求解器选择变量出现顺序中的第一个unbound变量。
• IntValueStrategy 选择下一个值的规则。默认值 ASSIGN_MIN_VALUE，表示选择还没有试的变量中的最小值。另一个选项 ASSIGN_MAX_VALUE 反之。

import sys
from ortools.constraint_solver import pywrapcp

def main():
  solver = pywrapcp.Solver("simple_example")
  # Create the variables.
  num_vals = 3
  x = solver.IntVar(0, num_vals - 1, "x")
  y = solver.IntVar(0, num_vals - 1, "y")
  z = solver.IntVar(0, num_vals - 1, "z")
  # Create the constraints.
  solver.Add(x != y)
  # Call the solver.
  db = solver.Phase([x, y, z], solver.CHOOSE_FIRST_UNBOUND, solver.ASSIGN_MIN_VALUE)
  solver.Solve(db)
  print_solution(solver, x, y, z)

def print_solution(solver, x, y, z):
  count = 0

  while solver.NextSolution():
    count += 1
    print("x =", x.Value(), "y =", y.Value(), "z =", z.Value())
  print("
Number of solutions found:", count)

if __name__ == "__main__":
  main()

打印结果:

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

Number of solutions found: 18

Cryptarithmetic Puzzles

Cryptarithmetic Puzzles 是一种数学问题，每个字母代表唯一数字，目标是找到这些数字，使给定方程成立。
由

      CP
+     IS
+    FUN
--------
=   TRUE

找到

      23
+     74
+    968
--------
=   1065

这是其中一个答案，还可以有其他解。

用 CP-SAT 或者 传统CP 求解器均可。

建模

约束构建如下:

• 等式 CP + IS + FUN = TRUE
• 每个字母代表一个数字
• C、I、F、T不能为0，因为开头不写0

这个问题官方解法代码有误，回头研究。

The N-queens Problem

n皇后问题是个很经典的问题。如果用深度优先搜索的解法见 link

In chess, a queen can attack horizontally, vertically, and diagonally. The N-queens problem asks:

How can N queens be placed on an NxN chessboard so that no two of them attack each other?

在国际象棋中，皇后可以水平、垂直和对角攻击。N-queens问题即: 如何把N个皇后放在一个NxN棋盘上，使它们不会互相攻击?

这并不是最优化问题，我们是要找到所有可行解。

CP求解器是尝试所有的可能来找到可行解。

传播和回溯

约束优化有两个关键元素:

• Propagation，传播可以加速搜索过程，因为减少了求解器继续做无谓的尝试。每次求解器给变量赋值时，约束条件都对为赋值的变量增加限制。
• Backtracking，回溯发生在继续向下搜索也肯定无解的情况下，故需要回到上一步尝试未试过的值。

构造约束:

• 约束禁止皇后区位于同一行、列或对角线上
• queens[i] = j means there is a queen in column i and row j.
• 设置在不同行列，model.AddAllDifferent(queens)
• 设置两个皇后不能在同一对角线，更为复杂些。
  • 如果对角线是从左到右下坡的，则斜率为-1（等同于有直线 y = -x + b），故满足 queens(i) + i 具有相同值的所有 queens(i) 在一条对角线上。
  • 如果对角线是从左到右上坡的，则斜率为 1（等同于有直线 y = x + b），故满足 queens(i) - i 具有相同值的所有 queens(i) 在一条对角线上。

import sys
from ortools.sat.python import cp_model

def main(board_size):
  model = cp_model.CpModel()
  # Creates the variables.
  # The array index is the column, and the value is the row.
  queens = [model.NewIntVar(0, board_size - 1, 'x%i' % i)
            for i in range(board_size)]
  # Creates the constraints.
  # The following sets the constraint that all queens are in different rows.
  model.AddAllDifferent(queens)

  # Note: all queens must be in different columns because the indices of queens are all different.

  # The following sets the constraint that no two queens can be on the same diagonal.
  for i in range(board_size):
    # Note: is not used in the inner loop.
    diag1 = []
    diag2 = []
    for j in range(board_size):
      # Create variable array for queens(j) + j.
      q1 = model.NewIntVar(0, 2 * board_size, 'diag1_%i' % i)
      diag1.append(q1)
      model.Add(q1 == queens[j] + j)
      # Create variable array for queens(j) - j.
      q2 = model.NewIntVar(-board_size, board_size, 'diag2_%i' % i)
      diag2.append(q2)
      model.Add(q2 == queens[j] - j)
    model.AddAllDifferent(diag1)
    model.AddAllDifferent(diag2)
  ### Solve model.
  solver = cp_model.CpSolver()
  solution_printer = SolutionPrinter(queens)
  status = solver.SearchForAllSolutions(model, solution_printer)
  print()
  print('Solutions found : %i' % solution_printer.SolutionCount())


class SolutionPrinter(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):
    self.__solution_count += 1
    for v in self.__variables:
      print('%s = %i' % (v, self.Value(v)), end = ' ')
    print()

  def SolutionCount(self):
    return self.__solution_count


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

  def OnSolutionCallback(self):
    self.__solution_count += 1

    for v in self.__variables:
      queen_col = int(self.Value(v))
      board_size = len(self.__variables)
      # Print row with queen.
      for j in range(board_size):
        if j == queen_col:
          # There is a queen in column j, row i.
          print("Q", end=" ")
        else:
          print("_", end=" ")
      print()
    print()

  def SolutionCount(self):
    return self.__solution_count


if __name__ == '__main__':
  # By default, solve the 8x8 problem.
  board_size = 8
  if len(sys.argv) > 1:
    board_size = int(sys.argv[1])
  main(board_size)

欲图形化打印，调用SolutionPrinter处换成DiagramPrinter。

Setting solver limits

时间限制

eg:

solver.parameters.max_time_in_seconds = 10.0

找到指定数量的解后停止搜索

eg:

from ortools.sat.python import cp_model


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

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

    def on_solution_callback(self):
        self.__solution_count += 1
        for v in self.__variables:
            print('%s=%i' % (v, self.Value(v)), end=' ')
        print()
        if self.__solution_count >= self.__solution_limit:
            print('Stop search after %i solutions' % self.__solution_limit)
            self.StopSearch()

    def solution_count(self):
        return self.__solution_count


def StopAfterNSolutionsSampleSat():
    """Showcases calling the solver to search for small number of 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 a solver and solve.
    solver = cp_model.CpSolver()
    solution_printer = VarArraySolutionPrinterWithLimit([x, y, z], 5)
    status = solver.SearchForAllSolutions(model, solution_printer)
    print('Status = %s' % solver.StatusName(status))
    print('Number of solutions found: %i' % solution_printer.solution_count())
    assert solution_printer.solution_count() == 5


StopAfterNSolutionsSampleSat()