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("\nNumber 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()