充足可能性問題と重み付き制約充足問題

重み付き制約充足問題

変数 \(x_j (j=1,2,\ldots,n)\) に対して有限集合から成る領域 (domain) \(D_j\) を定義する．

制約 \(C_i (i=1,2,\ldots,m)\) は変数の \(n\)　組の部分集合として定義される． 制約充足問題(constraint satisfaction problem)は，制約を満たすように，変数に領域の中の1つの値を割り当てることを目的とした問題である．

これの判定問題版（解があるか否かを判定する問題）が，最初の\(NP\)-完全問題として知られる制約充足問題 (satisfiability problem; SAT) である．

この問題は実行可能性判定問題の範疇に含まれるが，これに制約からの逸脱ペナルティを付加することによって，以下の組合せ最適化問題になる．

解ベクトル \(x\) が制約 \(C_i\) から逸脱している量を表すペナルティ関数を \(g_i (x)\) とする．各制約 \(C_i\) の重みを \(w_i\) としたとき，重み付き制約充足問題(weighted constraint satisfaction problem)は以下のように書ける．

\[ \begin{array}{lll} minimize & \sum_{i=1}^m w_i g_i(x) & \\ s.t. & x_j \in D_j & j=1,2,\ldots,n \end{array} \]

この問題に対しては，メタヒューリスティクスに基づく制約最適化ソルバーSCOP （Solver for COnstraint Programing）が開発されている． SCOPについての詳細は，付録1を参照されたい．

SCOPは実務的な問題を解くために開発されたものであるが，様々な基本的な組合せ最適化問題も簡単に解くことができる．例として，グラフ彩色問題のところで，SCOPを用いた数行の最適化コードを示している．

実務的な問題の例としては，シフトスケジューリングのことろで，看護婦スケジューリング問題のモデル化を示している．以下では，もう1つの実務的な例として，時間割作成問題への適用を考える．

時間割作成問題

通常の時間割作成の他に，大学では試験の時間割などにもニーズがある．

次の集合が与えられている．

授業（クラス）集合: 大学への応用の場合には，授業には担当教師が紐付けられている．
教室集合
学生集合
期集合：通常は1週間(5日もしくは6日）の各時限を考える．

以下の条件を満たす時間割を求める．

すべての授業をいずれか期へ割当
すべての授業をいずれか教室へ割当
各期、1つの教室では1つ以下の授業
同じ教員の受持ち講義は異なる期へ
割当教室は受講学生数以上の容量をもつ
同じ学生が受ける可能性がある授業の集合（カリキュラム）は，異なる期へ割り当てなければならない

考慮すべき付加条件には，以下のものがある．

1日の最後の期に割り当てられると履修学生数分のペナルティ
1人の学生が履修する授業が3連続するとペナルティ
1人の学生の授業数が、1日に1つならばペナルティ（\(0\) か \(2\) 以上が望ましい）

これは制約最適化問題として以下のように定式化できる．

授業 \(i\) の割り当て期を表す変数 \(X_i\) （領域はすべての期の集合）と割り当て教室を表す変数 \(Y_i\) （領域はすべての教室の集合）を用いる．

定式化では，これらの変数は値変数として表記する．それぞれ，授業 \(i\) を期 \(t\) に割り当てるとき \(1\) の \(0\)-\(1\) 変数 \(x_{it}\)，授業 \(i\) を教室 \(k\) に割り当てるとき \(1\) の \(0\)-\(1\) 変数 \(y_{ik}\) となる．

SCOPにおいては，以下の2つの制約は自動的に守られるので必要ない．

すべての授業をいずれか期へ割当
すべての授業をいずれか教室へ割当

他の制約は，以下のように記述できる．

各期 \(t\) の各教室 \(k\) への割り当て授業は1以下であることを表す制約： \[ \sum_{i} x_{it} y_{ik} \leq 1 \ \ \ \forall t,k \]
同じ教員の受持ち講義は異なる期へ割り当て：

教員 \(l\) の担当する授業の集合を \(E_l\) とすると，この制約は「\(X_{i} (i \in E_l)\) はすべて異なる値をもつ」と書くことができる． SCOPでは，このようなタイプの制約を相異制約とよび，そのまま記述できる．

割当教室は受講学生数以上の容量をもつ．

授業 \(i\) ができない（容量を超過する）教室の集合を \(K_i\) する．

\[ \sum_{i, k \in K_i} y_{ik} \leq 0 \]

同じ学生が受ける可能性がある授業の集合（カリキュラム）は，異なる期へ割り当てなければならない．

カリキュラム \(j\) に含まれる授業の集合を \(C_j\) としたとき，以下の制約として記述できる．

\[ \sum_{i \in C_j} x_{it} \leq 1 \ \ \ \forall j, t \]

1日の最後の期に割り当てられると履修学生数分のペナルティ

1日の最後の期の集合を \(L\)，授業 \(i\) の履修学生数を \(w_i\) とする． \[ \sum_{i, t \in L} w_i x_{it} \leq 0 \]

1人の学生が履修する授業が3連続すると1ペナルティ

\(T\) を1日のうちで最後の2時間でない期の集合とする． \[ \sum_{i \in C_j} x_{it} + x_{i,t+1} +x_{i,t+2} \leq 2 \ \ \ \forall j, t \in T \]

1人の学生の授業数が、1日に1つならば1ペナルティ（0か2以上が望ましい）

各日 \(d\) に含まれる期の集合を \(T_d\) とし，日 \(d\) におけるカリキュラム \(j\) に含まれる授業数が \(0\) か \(2\) 以上なのかを表す \(0\)-\(1\) 変数 \(z_{jd}\) を用いる．

\[ \sum_{i \in C_j} x_{it} \leq |T_d| z_{jd} \ \ \ \forall d, j \] \[ \sum_{i \in C_j} x_{it} \geq 2 z_{jd} \ \ \ \forall d, j \]

上の定式化をSCOPで解くことによって，付加的な制約を「できるだけ」満たす時間割を作成することができる．

時間割作成の国際コンペティション (ITC2007) では， SCOPは3つの異なるトラックですべてメダル（3位, 2位, 3位）を獲得している．

OR-Tools (cp-sat)

Googleが開発している制約最適化ソルバー OR-Tools (cp-sat)でも同様のことができる．cp_modelはSCOPとは異なる求解手法である制約伝播を用いているため，パズルのような問題を解くことが得意である． OR-Toolsについての詳細は，以下を参照されたい．

https://developers.google.com/optimization

以下に，簡単な最適化の例を示す．

変数

cp-satの変数はモデルインスタンスのNewIntVar(もしくはnew_int_var）メソッドで生成できる．引数は下限 lb，上限 ub，名前 name である．必ず下限と上限を設定するのが特徴で，変数も実数変数でなく，必ず整数値をとる変数を用いる必要がある．

NewBoolVarメソッドを使うとTrueもしくはFalseの値をとる論理変数が定義できる．以下で生成したbという論理変数に対して，b.Notメソッドで生成される変数は，bがTrueのときFalse，FalseのときTrueになる．

from ortools.sat.python import cp_model
model = cp_model.CpModel()
z = model.NewIntVar(-100, 100, "z")
# z = model.new_int_var(-100, 100, "z") #これでも同じだが，以下ではCamel Stypeに統一する．

# boolean variable b
b = model.NewBoolVar("b")
# implicitly available negation of b:
not_b = b.Not()  # will be 1 if b is 0 and 0 if b is 1
not_b

<ortools.sat.python.cp_model._NotBooleanVariable>

連続した整数値だけでなく，与えられた整数値のリストの値だけをとれる変数も定義できる．

それにはまずDomainクラスのインスタンスを生成し，それを引数としてモデルのNewIntVarFromDomainメソッドで生成する．

# Define a domain with selected values
domain = cp_model.Domain.FromValues([2, 5, 8, 10, 20, 50, 90])

x = model.NewIntVarFromDomain(domain, "x")
x

x(2, 5, 8, 10, 20, 50, 90)

目的関数

目的関数の最小化（最大化）はモデルの Minimize (Maximize)メソッドで定義する．以下の例では論理変数 \(x_i (i=0,1,\ldots,9)\) に対して，添字 \(i\) が偶数ならTrueのとき係数を \(i\) とし，奇数ならFalseのとき係数を \(i\) とした目的関数を最小化している．

最適化は， CpSolverのインスタンス solver を準備し，モデルインスタンス model をSolveメソッドの引数として与えることによって実行される．返値はstatusであり，定数OPTIMAL(4)のとき最適解が得られたことになる．最適目的関数値は，ソルバーインスタンスのObjectiveValueメソッドで，最適解はValueメソッドに変数インスタンスを渡すことによって得られる．

model = cp_model.CpModel()

x_vars = [model.NewBoolVar(f"x{i}") for i in range(10)]
model.Minimize(
    sum(i * x_vars[i] if i % 2 == 0 else i * x_vars[i].Not() for i in range(10))
)
solver = cp_model.CpSolver()
status = solver.Solve(model)
if status == cp_model.OPTIMAL:
    print('Optimal Value:', solver.ObjectiveValue())
    for i in range(10):
        print( solver.Value( x_vars[i]), end="  ")
else:
    print('Solver exited with nonoptimal status:', status)

Optimal Value: 0.0
0  1  0  1  0  1  0  1  0  1

制約

model = cp_model.CpModel()

num_vals = 3
x = model.new_int_var(0, num_vals - 1, "x")
y = model.new_int_var(0, num_vals - 1, "y")
z = model.new_int_var(0, num_vals - 1, "z")

#model.add(x + z == 3 * y)

model.add(x + y != z)

model.add(x < z)  

solver = cp_model.CpSolver()
status = solver.solve(model)

model.Minimize(x+y+z)

if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
    print(f"x = {solver.value(x)}")
    print(f"y = {solver.value(y)}")
    print(f"z = {solver.value(z)}")
else:
    print('Solver exited with nonoptimal status:', status)

x = 0
y = 0
z = 1

すべての解の列挙

cp-satでは，1つの最適解を出すのではなく，すべての実行可能解を列挙することもできる．

列挙した解はコールバック関数を渡すことによって得ることができる． cp_modelにはVarArraySolutionPrinterクラスが準備されているので，それを用いても良いが，以下では探索した解の個数をカウントする自前のコールバック関数を作成する．これは， CpSolverSolutionCallbackクラスから派生させ，コンストラクタで変数のリストを受け取り，新たな実行可能解がみつかると解の個数を増やし，solution_countメソッドで見つかった解の個数を返すように変更してある．

解の列挙はソルバーインスタンスのSearchForAllSolutionsメソッドで行う．引数はモデルインスタンス model とコールバック関数のインスタンスである．

以下の例では，変数 \(x,y,z\) がすべて異なる値をとる制約を課しているが，これは相異制約 AddAllDifferent を用いても記述できる．相異制約の引数は，変数のリストである．ただし，少数の変数に対しては相異制約は制約伝播を妨げるので，通常の異なることを表す制約 \(!=\) を用いた方が良い．

class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):
    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(v, "=", self.Value(v), end=" ")
        print()

    def solution_count(self):
        return self.__solution_count

model = cp_model.CpModel()

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

model.Add(x != y)
model.Add(y != z)
model.Add(x != z)
# model.AddAllDifferent([x,y,z])

solver = cp_model.CpSolver()
solution_printer = VarArraySolutionPrinter([x, y, z])
status = solver.SearchForAllSolutions(model, solution_printer)

print("Status =", solver.StatusName(status))
print("Number of solutions found", solution_printer.solution_count())

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

例題：覆面算

各文字に0から９の数字を入れて等式 \(SEND+MORE=MONEY\) を成立させたい．ただし，数字に重複があってはならず，先頭の文字に0を入れることはできない．

model = cp_model.CpModel()
S = model.NewIntVar(1,9,"S")
E = model.NewIntVar(0,9,"E")
N = model.NewIntVar(0,9,"N")
D = model.NewIntVar(0,9,"D")
M = model.NewIntVar(1,9,"M")
O = model.NewIntVar(0,9,"O")
R = model.NewIntVar(0,9,"R")
Y = model.NewIntVar(0,9,"Y")

model.Add(   1000*S + 100*E + 10*N + D + 
             1000*M + 100*O + 10*R + E 
 == 10000*M +1000*O + 100*N + 10*E + Y)

model.AddAllDifferent([S,E,N,D,M,O,R,Y])

solver = cp_model.CpSolver()
status = solver.Solve(model)

if status == cp_model.OPTIMAL:
    for v in [S,E,N,D,M,O,R,Y]:
        print(f"{v} = {solver.Value(v)}")

S = 9
E = 5
N = 6
D = 7
M = 1
O = 0
R = 8
Y = 2

#唯一の解であることの確認
solution_printer = VarArraySolutionPrinter([S,E,N,D,M,O,R,Y])
status = solver.SearchForAllSolutions(model, solution_printer)

print("Status =", solver.StatusName(status))
print("Number of solutions found:", solution_printer.solution_count())

S = 9 E = 5 N = 6 D = 7 M = 1 O = 0 R = 8 Y = 2 
Status = OPTIMAL
Number of solutions found: 1

練習問題：覆面算

https://www.math.uni-bielefeld.de/~sillke/PUZZLES/ALPHAMETIC/ にある覆面算の問題をどれか一つ解いてみよ．

例題：魔方陣

魔方陣とは，\(n \times n\) の正方形の方陣に1から\(n^2\)までの整数を1つずつ入れて，縦・横・対角線のいずれの列の和も同じになるものをいう．（一列の和は \(n(n^2+1)/2\) になる．）

注意： \(n=12\) より大きい値で試さないこと．

n = 10
model = cp_model.CpModel()
x = {}
for i in range(n):
    for j in range(n):
        x[i, j] = model.NewIntVar(1, n * n, f"x({i},{j})")

x_list = [x[i,j] for i in range(n) for j in range(n)]
model.AddAllDifferent(x_list)

s = n*(n**2+1)//2

for i in range(n):
    model.Add(sum([x[i, j] for j in range(n)]) == s) 
for j in range(n):
    model.Add(sum([x[i, j] for i in range(n)]) == s) 

model.Add(sum([x[i, i] for i in range(n)]) == s) 
model.Add(sum([x[i, n - i - 1] for i in range(n)]) == s)

solver = cp_model.CpSolver()
status = solver.Solve(model)

for i in range(n):
    for j in range(n):
        print(solver.Value(x[i,j]), end=" ")
    print()

65 1 60 39 51 90 73 33 78 15 
71 20 77 67 80 25 21 74 53 17 
69 59 12 8 64 88 79 76 40 10 
27 86 32 26 99 24 52 31 35 93 
92 83 41 56 19 61 34 23 46 50 
96 85 82 49 81 5 37 48 18 4 
6 58 62 87 14 84 68 42 54 30 
22 75 16 72 7 11 44 97 70 91 
2 9 57 38 43 28 94 36 98 100 
55 29 66 63 47 89 3 45 13 95

練習問題：完全方陣

\(n=4\) の場合の完全方陣（魔方陣でかつ四隅の合計と任意の\(2\times2\)の正方形の合計も同じ和になるもの）を求めよ．

例題：数独

数独は \(9 \times 9\) の正方形の枠内に1から\(n\)までの数字を入れるパズルである．初期配置に与えられた数字はそのままとし，縦・横の各列とブロックとよばれる \(3 \times 3\) の小正方形の枠内には，同じ数字を重複して入れてはいけないものとする．

以下のデータでは，数字の０が入っている枠は空白であるとする．

import math
problem = [
[1,0,0, 0,0,7, 0,9,0],
[0,3,0, 0,2,0, 0,0,8],
[0,0,9, 6,0,0, 5,0,0],

[0,0,5, 3,0,0, 9,0,0],
[0,1,0, 0,8,0, 0,0,2],
[6,0,0, 0,0,4, 0,0,0],

[3,0,0, 0,0,0, 0,1,0],
[0,4,0, 0,0,0, 0,0,7],
[0,0,7, 0,0,0, 3,0,0]]

model = cp_model.CpModel()  
n = len(problem)

cell_size = math.ceil(math.sqrt(n))
line_size = cell_size ** 2
line = range(0, line_size)
cell = range(0, cell_size)

x = {}
for i in line:
    for j in line:
        x[i, j] = model.NewIntVar(1, line_size, f"x({i},{j})")

for i in line:
    model.AddAllDifferent([x[i, j] for j in line])

for j in line:
    model.AddAllDifferent([x[i, j] for i in line])

for i in cell:
    for j in cell:
        one_cell = []
        for di in cell:
            for dj in cell:
                one_cell.append(x[(i * cell_size + di, j * cell_size + dj)])
        model.AddAllDifferent(one_cell)
    
for i in line:
    for j in line:
        if problem[i][j]:
            model.Add(x[i, j] == problem[i][j])

solver = cp_model.CpSolver()
status = solver.Solve(model)

for i in line:
    for j in line:
        print(solver.Value(x[i,j]), end=" ")
    print()

1 6 2 8 5 7 4 9 3 
5 3 4 1 2 9 6 7 8 
7 8 9 6 4 3 5 2 1 
4 7 5 3 1 2 9 8 6 
9 1 3 5 8 6 7 4 2 
6 2 8 7 9 4 1 3 5 
3 5 6 4 7 8 2 1 9 
2 4 1 9 3 5 8 6 7 
8 9 7 2 6 1 3 5 4

#唯一の解であることの確認
solution_printer = VarArraySolutionPrinter([x[i,j] for i in line for j in line])
status = solver.SearchForAllSolutions(model, solution_printer)

print("Status", solver.StatusName(status))
print("Number of solutions found",  solution_printer.solution_count())

x(0,0) = 1 x(0,1) = 6 x(0,2) = 2 x(0,3) = 8 x(0,4) = 5 x(0,5) = 7 x(0,6) = 4 x(0,7) = 9 x(0,8) = 3 x(1,0) = 5 x(1,1) = 3 x(1,2) = 4 x(1,3) = 1 x(1,4) = 2 x(1,5) = 9 x(1,6) = 6 x(1,7) = 7 x(1,8) = 8 x(2,0) = 7 x(2,1) = 8 x(2,2) = 9 x(2,3) = 6 x(2,4) = 4 x(2,5) = 3 x(2,6) = 5 x(2,7) = 2 x(2,8) = 1 x(3,0) = 4 x(3,1) = 7 x(3,2) = 5 x(3,3) = 3 x(3,4) = 1 x(3,5) = 2 x(3,6) = 9 x(3,7) = 8 x(3,8) = 6 x(4,0) = 9 x(4,1) = 1 x(4,2) = 3 x(4,3) = 5 x(4,4) = 8 x(4,5) = 6 x(4,6) = 7 x(4,7) = 4 x(4,8) = 2 x(5,0) = 6 x(5,1) = 2 x(5,2) = 8 x(5,3) = 7 x(5,4) = 9 x(5,5) = 4 x(5,6) = 1 x(5,7) = 3 x(5,8) = 5 x(6,0) = 3 x(6,1) = 5 x(6,2) = 6 x(6,3) = 4 x(6,4) = 7 x(6,5) = 8 x(6,6) = 2 x(6,7) = 1 x(6,8) = 9 x(7,0) = 2 x(7,1) = 4 x(7,2) = 1 x(7,3) = 9 x(7,4) = 3 x(7,5) = 5 x(7,6) = 8 x(7,7) = 6 x(7,8) = 7 x(8,0) = 8 x(8,1) = 9 x(8,2) = 7 x(8,3) = 2 x(8,4) = 6 x(8,5) = 1 x(8,6) = 3 x(8,7) = 5 x(8,8) = 4 
Status OPTIMAL
Number of solutions found 1

練習問題：数独

数独（ナンプレ）の問題を検索して探し，上のコードで解いてみよ．

練習問題：不等式パズル

\(5 \times 5\) の正方形の枠内に1から5の整数を入れる．初期配置に与えられた数字はそのままとし，さらに2つの格子間に不等式 \(<\) が付加されている．ただし，各行・各列ともに数字の重複があってはいけない．

初期配置を入れた行列（リストのリスト）においては，0は空白を意味する．

不等式を表すリストは，長さ4のリスト \([i,j,k,l]\) から構成され，これは\(i\)行\(j\)列の数字が\(k\)行\(l\)列の数字より小さいことを表す．

以下の例題は， https://en.wikipedia.org/wiki/Futoshiki からとったものである．

#初期配置
initial = [
    [0, 0, 0, 0, 0],
    [4, 0, 0, 0, 2],
    [0, 0, 4, 0, 0],
    [0, 0, 0, 0, 4],
    [0, 0, 0, 0, 0]]
#不等式（添字は1から始まるものとする）
inequality = [
    [1, 2, 1, 1],
    [1, 4, 1, 3],
    [1, 5, 1, 4],
    [4, 4, 4, 5],
    [5, 1, 5, 2],
    [5, 2, 5, 3]
]

例題：nクイーン

import time

class NQueenSolutionPrinter(cp_model.CpSolverSolutionCallback):
    def __init__(self, queens: list[cp_model.IntVar]):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__queens = queens
        self.__solution_count = 0
        self.__start_time = time.time()

    def solution_count(self) -> int:
        return self.__solution_count

    def on_solution_callback(self) -> None:
        current_time = time.time()
        print(
            "Solution %i, time = %f s"

        )
        self.__solution_count += 1

        all_queens = range(len(self.__queens))
        for i in all_queens:
            for j in all_queens:
                if self.value(self.__queens[j]) == i:
                    # There is a queen in column j, row i.
                    print("Q", end=" ")
                else:
                    print("_", end=" ")
            print()
        print()

board_size = 5
model = cp_model.CpModel()

queens = [
    model.new_int_var(0, board_size - 1, "x%i" % i) for i in range(board_size)
]

# All columns must be different because the indices of queens are all
# different, so we just add the all different constraint on the rows.
model.add_all_different(queens)

# No two queens can be on the same diagonal.
diag1 = []
diag2 = []
for i in range(board_size):
    q1 = model.new_int_var(0, 2 * board_size, "diag1_%i" % i)
    q2 = model.new_int_var(-board_size, board_size, "diag2_%i" % i)
    diag1.append(q1)
    diag2.append(q2)
    model.add(q1 == queens[i] + i)
    model.add(q2 == queens[i] - i)
model.add_all_different(diag1)
model.add_all_different(diag2)

solver = cp_model.CpSolver()
solution_printer = NQueenSolutionPrinter(queens)
solver.parameters.enumerate_all_solutions = True
status = solver.solve(model, solution_printer)

Solution 0, time = 0.003772 s
_ _ Q _ _ 
Q _ _ _ _ 
_ _ _ Q _ 
_ Q _ _ _ 
_ _ _ _ Q 

Solution 1, time = 0.004260 s
_ _ _ Q _ 
Q _ _ _ _ 
_ _ Q _ _ 
_ _ _ _ Q 
_ Q _ _ _ 

Solution 2, time = 0.004526 s
_ Q _ _ _ 
_ _ _ Q _ 
Q _ _ _ _ 
_ _ Q _ _ 
_ _ _ _ Q 

Solution 3, time = 0.004755 s
_ Q _ _ _ 
_ _ _ _ Q 
_ _ Q _ _ 
Q _ _ _ _ 
_ _ _ Q _ 

Solution 4, time = 0.004973 s
_ _ Q _ _ 
_ _ _ _ Q 
_ Q _ _ _ 
_ _ _ Q _ 
Q _ _ _ _ 

Solution 5, time = 0.005197 s
Q _ _ _ _ 
_ _ _ Q _ 
_ Q _ _ _ 
_ _ _ _ Q 
_ _ Q _ _ 

Solution 6, time = 0.005414 s
Q _ _ _ _ 
_ _ Q _ _ 
_ _ _ _ Q 
_ Q _ _ _ 
_ _ _ Q _ 

Solution 7, time = 0.005625 s
_ _ _ _ Q 
_ _ Q _ _ 
Q _ _ _ _ 
_ _ _ Q _ 
_ Q _ _ _ 

Solution 8, time = 0.005841 s
_ _ _ _ Q 
_ Q _ _ _ 
_ _ _ Q _ 
Q _ _ _ _ 
_ _ Q _ _ 

Solution 9, time = 0.006049 s
_ _ _ Q _ 
_ Q _ _ _ 
_ _ _ _ Q 
_ _ Q _ _ 
Q _ _ _ _

例題： ABCプレース

ABCプレースは正方形の枠内にA,B,Cの文字を入れるパズルである．ルールは以下の通り．

各行各列にA,B,Cの文字を1つずつ記入する（文字が入らないマスもある）。
枠外にある文字は、その列で最も辺に近い位置にある文字である．

以下の例は6行6列の問題である． A,B,Cは整数0,1,2で表し，以下の整数変数を用いる．

vh[i,k]: 各行i，各文字k (=0,1,2) に対して，文字が入る列の番号を変数とした整数変数
vv[i,k]: 各列i，各文字k (=0,1,2) に対して，文字が入る行の番号を変数とした整数変数

from ortools.sat.python import cp_model
import numpy as np
from itertools import product

data = """\
 B B  C 
A......B
B...... 
C...... 
 ......A
 ......C
 ...... 
  BC CA """.split('\n')
nn = len(data)-2

def chk(v, c, isless):
    k = ord(c) - 65
    if k >= 0:
        for j in range(3):
            if j != k:
                model.Add(v[k] < v[j] if isless else v[k] > v[j])
                
model = cp_model.CpModel()

vh = np.array([[model.NewIntVar(0, nn - 1, f'vh{i}{j}') for j in range(3)] for i in range(nn)])
vv = np.array([[model.NewIntVar(0, nn - 1, f'vv{i}{j}') for j in range(3)] for i in range(nn)])

for i in range(nn):
    model.AddAllDifferent(vh[i])
    model.AddAllDifferent(vv[i])
    chk(vv[i], data[0][i + 1], True)
    chk(vv[i], data[nn + 1][i + 1], False)
    chk(vh[i], data[i + 1][0], True)
    chk(vh[i], data[i + 1][nn + 1], False)
    
for i, j, k in product(range(nn), range(nn), range(3)):
    vb1, vb2 = model.NewBoolVar(f'vb1{i}{j}{k}'), model.NewBoolVar(f'vb2{i}{j}{k}')
    model.Add(vh[i][k] == j).OnlyEnforceIf(vb1)
    model.Add(vv[j][k] != i).OnlyEnforceIf(vb2)
    model.AddBoolXOr([vb1, vb2])
    
solver = cp_model.CpSolver()
status = solver.Solve(model)

res = np.vectorize(solver.Value)(vh)
cc = np.full((nn, nn), '.')
for i, j in product(range(nn), range(3)):
    cc[i][res[i][j]] = chr(65 + j)
print('\n'.join(' '.join(s) for s in cc))

. A . C B .
. . B . A C
. C A . . B
B . C . . A
A . . B C .
C B . A . .

例題：嘘つき島パズル

ある島には正直族と嘘つき族と呼ばれる2 種類の人たちが仲良く住んでいる．正直族は必ず本当のことを言い，嘘つき族は必ず嘘をつく．

あなたは，この島の人たちが正直族か嘘つき族なのかの調査を依頼された．

最初の家の旦那に聞いたところ「夫婦は両方とも嘘つき族だよ」という答えだった．
次の家に行って旦那に「ご夫婦は両方とも嘘つき族ですか？」と聞いたところ「少なくとも1人はね」という答えだった．
その次の家に行ったところ旦那が答えた「もし私が正直族なら，私の家内も正直族です」

さて，上の情報から，調査結果を報告せよ．

整数最適化でも定式化できるが，論理的な制約（AddBoolAnd, AddBoolOr, OnlyEnforceIf）を使うと，より簡潔に記述できる．

from ortools.sat.python import cp_model
# 1. 最初の家の旦那に聞いたところ「夫婦は両方とも嘘つき族だよ」という答えだった．
model = cp_model.CpModel()
x = model.NewBoolVar("x")  #旦那が正直ならTrueの論理変数
y = model.NewBoolVar("y")  #奥さんが正直ならTrueの論理変数
model.AddBoolAnd([x.Not(),y.Not()]).OnlyEnforceIf(x)  #旦那が正直なら「夫婦ともに嘘つき」
model.AddBoolOr([x,y]).OnlyEnforceIf(x.Not())         #旦那が嘘つきなら「夫婦ともに嘘つき」ではない => いずれかが正直
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.Value(x), solver.Value(y))  # 0  1 => 旦那は嘘つき，奥さんは正直 

#唯一解であることの確認
solution_printer = VarArraySolutionPrinter([x,y])
status = solver.SearchForAllSolutions(model, solution_printer)
print("Status =", solver.StatusName(status))
print("Number of solutions found:", solution_printer.solution_count())

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

#2. 旦那に「ご夫婦は両方とも嘘つき族ですか？」と聞いたところ「少なくとも1人はね」という答えだった．
model = cp_model.CpModel()
x = model.NewBoolVar("x")  #旦那が正直ならTrueの論理変数
y = model.NewBoolVar("y")  #奥さんが正直ならTrueの論理変数
model.AddBoolOr([x.Not(),y.Not()]).OnlyEnforceIf(x)  #旦那が正直なら「いずれかが嘘つき」
model.AddBoolAnd([x,y]).OnlyEnforceIf(x.Not())       #旦那が嘘つきなら「両方とも正直」
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.Value(x), solver.Value(y))  # 1  0 => 旦那は正直，奥さんは嘘つき 

#唯一解であることの確認
solution_printer = VarArraySolutionPrinter([x,y])
status = solver.SearchForAllSolutions(model, solution_printer)
print("Status =", solver.StatusName(status))
print("Number of solutions found:", solution_printer.solution_count())

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

# 3. 旦那が答えた「もし私が正直族なら，私の家内も正直族です」
model = cp_model.CpModel()
x = model.NewBoolVar("x")  #旦那が正直ならTrueの論理変数
y = model.NewBoolVar("y")  #奥さんが正直ならTrueの論理変数
model.Add(x<=y).OnlyEnforceIf(x)  #旦那が正直なら「旦那が正直なら奥さんも正直」 => x<=y 
model.Add(x>y).OnlyEnforceIf(x.Not()) #旦那が嘘つきなら「旦那が正直なら奥さんも正直」でない => x>y
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(solver.Value(x), solver.Value(y))  # 1  0 => 旦那は正直，奥さんは嘘つき 

#唯一解であることの確認
solution_printer = VarArraySolutionPrinter([x,y])
status = solver.SearchForAllSolutions(model, solution_printer)
print("Status =", solver.StatusName(status))
print("Number of solutions found:", solution_printer.solution_count())

1 1
x = 1 y = 1 
Status = OPTIMAL
Number of solutions found: 1

練習問題嘘つき島パズル (1)

ある島には正直族と嘘つき族と呼ばれる2 種類の人たちが仲良く住んでいる．正直族は必ず本当のことを言い，嘘つき族は必ず嘘をつく．

またこの島には，夜になると狼に変身して人を襲う狼男が紛れ込んでいる．狼男もこの島の住民なので，正直族か嘘つき族のいずれかに属する．あなたは，この島の人たちが狼男なのかの調査を依頼された．

3人のうち1人が狼男であることが分かっている． A,B,Cの3人組への証言は以下の通りである．

A: 「わたしは狼男です．」
B: 「わたしも狼男です．」
C: 「わたしたちの中の高々1人が正直族です．」

さて，狼男は誰か？また誰が正直族で誰が嘘つき族か？

練習問題嘘つき島パズル (2)

同じ島でまた別の3 人組 A,B,C の証言を得た．

A: 「わたしたちの中の少なくとも1 人が嘘つき族です．」
B: 「Cさんは正直族です．」

この3 人組も彼らのうちの1人が狼男で，彼は正直族であることが分かっているとき，狼男は誰か？また誰が正直族で誰が嘘つき族か？

在庫制約

モデルインスタンスの AddReservoirConstraint メソッドでタンクなどの貯蔵施設に対する在庫制約を定義できる．

引数は以下の通り．

times: 在庫を増加・減少させる時刻を表す非負の整数変数（もしくは整数）のリスト
demands: 在庫を増加・減少させる量を表す整数変数（もしくは整数）のリスト
min_level: 在庫の下限値
max_level: 在庫の上限値

以下の例では在庫の増加するための生産設備には，機械の制約の下で，需要（在庫の減少）はなるべく早くするものとする．

model = cp_model.CpModel()

#計画期間
horizon = 20

# 生産・消費量（生産時間は需要量に等しいものとする）
demands = [3, 5, 4, -4, -4, -4]
num_demands = len(demands) 

# 需要の発生時刻
times = [model.NewIntVar(0,10,f'times_{i}') for i in range(num_demands)]

J = [0,1,2] #生産するためのジョブ
end, interval = {}, {} #ジョブの終了時刻と区間を表す変数
for j in J:
    end[j] = model.NewIntVar(0, horizon, f"end({j})")
    interval[j] = model.NewIntervalVar(times[j], demands[j], end[j], f"job({j})")

#生産容量（機械）の制約
model.AddNoOverlap([interval[j] for j in J])

# タンクの在庫制約を追加（タンクの最小容量は0，最大容量は5）
model.AddReservoirConstraint(
    times,  
    demands,  
    0,  # 最小容量
    5  # 最大容量
)

#なるべく早く消費する
model.Minimize(sum(times[j] for j in [3,4,5]))

solver = cp_model.CpSolver()
status = solver.Solve(model)

print("times, deamnds")
for i in range(num_demands):
    print(solver.Value(times[i]), demands[i])

times, deamnds
9 3
4 5
0 4
0 -4
4 -4
9 -4

状態の推移を表すオートマトン制約

モデルインスタンスの AddAutomaton メソッドで有限状態オートマトンの状態推移を定義できる．

引数は以下の通り．

transition_variables: 変数のリスト；値はオートマトンの状態の推移（選択されたグラフの枝）を表す．
starting_state: オートマトンの初期状態
final_states: 最終状態のリスト
transition_triples: オートマトンの状態推移を表すリスト；現在の状態，変数の値，新しい状態のタプルを要素とする．

以下に2つの例を示す．最初の例は自明な状態推移である．もう1つの例は，固定された長さの文字列が “101” のパターンを含む最終状態に到達できるような状態遷移を求めるものである．

model = cp_model.CpModel()

# 変数の数を定義
num_vars = 3

# 0から2までの整数値を取ることができる変数を定義
vars = [model.NewIntVar(0, 2, f"var_{i}") for i in range(num_vars)]

# 初期状態
start_state = 0

# 最終状態
final_states = [0]

# 状態遷移
# (from_state, input, to_state) のリスト
transition_triples = [
    (0, 1, 1),  # 状態0から1を入力して状態1へ
    (1, 2, 2),  # 状態1から2を入力して状態2へ
    (2, 0, 0)   # 状態2から0を入力して状態0へ
]

# オートマトン制約を追加
model.AddAutomaton(vars, start_state, final_states, transition_triples)

solver = cp_model.CpSolver()
status = solver.Solve(model)

for i in range(num_vars):
    print(f"Variable {i}: {solver.Value(vars[i])}")

Variable 0: 1
Variable 1: 2
Variable 2: 0

model = cp_model.CpModel()

# 文字列の長さを設定
string_length = 4

# 0または1を取ることができる変数を定義
vars = [model.NewIntVar(0, 1, f"var_{i}") for i in range(string_length)]

# オートマトンの状態
# 状態0: パターン「101」をまだ見つけていない状態
# 状態1: 「1」を見つけた状態
# 状態2: 「10」を見つけた状態
# 状態3: 「101」を見つけた状態（最終状態）
start_state = 0
final_states = [3]

# 状態遷移
# (from_state, input, to_state) のリスト
transition_triples = [
    (0, 0, 0),  # 状態0から0を入力してそのまま状態0
    (0, 1, 1),  # 状態0から1を入力して状態1へ
    (1, 0, 2),  # 状態1から0を入力して状態2へ
    (2, 1, 3),  # 状態2から1を入力して状態3へ
    (1, 1, 1),  # 状態1から1を入力してそのまま状態1
    (2, 0, 0),  # 状態2から0を入力して状態0へ
    (3, 0, 3),  # 状態3は最終状態なので、どんな入力でもそのまま
    (3, 1, 3),  # 状態3は最終状態なので、どんな入力でもそのまま
]

# オートマトン制約を追加
model.AddAutomaton(vars, start_state, final_states, transition_triples)

solver = cp_model.CpSolver()

status = solver.Solve(model)

result = "".join(str(solver.Value(var)) for var in vars)
print("Generated binary string:", result)

#列挙
solution_printer = VarArraySolutionPrinter(vars)
status = solver.SearchForAllSolutions(model, solution_printer)

Generated binary string: 1011
var_0 = 1 var_1 = 0 var_2 = 1 var_3 = 0 
var_0 = 1 var_1 = 1 var_2 = 0 var_3 = 1 
var_0 = 1 var_1 = 0 var_2 = 1 var_3 = 1 
var_0 = 0 var_1 = 1 var_2 = 0 var_3 = 1

要素制約

モデルインスタンスの AddElement メソッドで，目標となる変数を整数変数で与えたインデックスに対応する変数（もしくは定数）配列に一致させることができる．

引数は以下の通り．

index: 整数変数で与えたインデックス
variables: 変数のリスト（定数を含んでも良い）
target: 目標となる変数

制約としては variables[index] == target であるが， indexは整数変数でも良い点が，通常の制約と異なる．

以下の例は，indexを0から4まで変化させたとき，変数の対応する値に目標変数が一致することを確認している．

from ortools.sat.python import cp_model

for i in range(5):
    model = cp_model.CpModel()
    
    # インデックスの値として0から4までの範囲
    index = model.NewIntVar(0, 4, 'index')
    
    # 参照される変数配列を定義する
    var_array = [model.NewIntVar(0,50,f"x({i})") for i in range(5)]

    # 変数の値を固定
    model.add( var_array[i]== i*10 )
        
    # 取得される値のための変数を定義する
    selected_value = model.NewIntVar(0, 50, 'selected_value')
    
    # AddElement制約を使用して、インデックスに基づいて選択された値を制約する
    model.AddElement(index, var_array, selected_value)
    
    solver = cp_model.CpSolver()

    # インデックスを特定の値に設定
    model.Add(index == i)

    status = solver.Solve(model)
    
    print(f'Index: {i}, Selected Value: {solver.Value(selected_value)}')

Index: 0, Selected Value: 0
Index: 1, Selected Value: 10
Index: 2, Selected Value: 20
Index: 3, Selected Value: 30
Index: 4, Selected Value: 40

インバース制約

通常の整数最適化では，\(n\)個の要素の完全マッチング（割当）をモデル化するには， \(n^2\)個の \(0\)-\(1\) 変数が必要となる．制約最適化では，マッチング制約を \(2n\)個の整数変数で記述できる．

モデルインスタンスの AddInverse メソッドで，2つの整数の変数配列が互いに割当になっていることを定義できる．

引数は以下の通り． - variables: 割当を表す整数変数の配列 - inverse_variables: 逆割当を表す整数変数の配列

制約として記述すると variables[i] = j \(\Leftrightarrow\) inverse_variables[j] = i となる．

以下の例では，5つの仕事を5人の作業員に割り振るときの最小費用の割当を求めている．仕事からみた作業員の割当と，作業員からみた仕事の割当が正しく定義されていることが確認できる．

from ortools.sat.python import cp_model

# タスクと労働者の数
num_tasks = 5
num_workers = 5

# 割当用の変数を定義
tasks = [model.NewIntVar(0, num_workers - 1, f'task_{i}') for i in range(num_tasks)]
workers = [model.NewIntVar(0, num_tasks - 1, f'worker_{i}') for i in range(num_workers)]

# インバース制約を追加して、タスクと労働者の間の一対一の関係を強制
model.AddInverse(tasks, workers)

# コスト行列を定義
cost_matrix = [
    [10, 19, 8, 15, 5],  # タスク0の各労働者への割当コスト
    [10, 18, 7, 17, 3],  # タスク1の各労働者への割当コスト
    [13, 16, 9, 14, 4],  # タスク2の各労働者への割当コスト
    [12, 17, 6, 10, 11], # タスク3の各労働者への割当コスト
    [14, 15, 9, 8, 2],   # タスク4の各労働者への割当コスト
]

# コストを計算する変数
total_cost = model.NewIntVar(0, 100, 'total_cost')

# タスク割り当てのコストを計算
costs = [model.NewIntVar(0, 20, f'cost_{i}') for i in range(num_tasks)]

# 各タスクの割当コストを定義
for i in range(num_tasks):
    model.AddElement(tasks[i], cost_matrix[i], costs[i])

model.Add(total_cost == sum(costs))

# 目的関数：合計コストを最小化
model.Minimize(total_cost)

solver = cp_model.CpSolver()
status = solver.Solve(model)

if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
    print("Total cost:", solver.Value(total_cost))

    print("Task assignments:")
    for i in range(num_tasks):
        print(f"  Task {i} -> Worker {solver.Value(tasks[i])}")

    print("Worker assignments:")
    for i in range(num_workers):
        print(f"  Worker {i} -> Task {solver.Value(workers[i])}")

Total cost: 43
Task assignments:
  Task 0 -> Worker 0
  Task 1 -> Worker 4
  Task 2 -> Worker 1
  Task 3 -> Worker 2
  Task 4 -> Worker 3
Worker assignments:
  Worker 0 -> Task 0
  Worker 1 -> Task 2
  Worker 2 -> Task 3
  Worker 3 -> Task 4
  Worker 4 -> Task 1

重み付き制約充足問題

時間割作成問題

OR-Tools (cp-sat)

変数

目的関数

制約

すべての解の列挙

例題：覆面算

練習問題：覆面算

例題：魔方陣

練習問題：完全方陣

例題：数独

練習問題：数独

練習問題：不等式パズル

例題：nクイーン

例題： ABCプレース

例題：嘘つき島パズル

練習問題 嘘つき島パズル (1)

練習問題 嘘つき島パズル (2)

在庫制約

状態の推移を表すオートマトン制約

要素制約

インバース制約

練習問題嘘つき島パズル (1)

練習問題嘘つき島パズル (2)