重み付き制約充足問題
変数 $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
以下に，簡単な最適化の例を示す．
from ortools.sat.python import cp_model
model = cp_model.CpModel()
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")
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)
solver = cp_model.CpSolver()
status = solver.Solve(model)
if status == cp_model.OPTIMAL:
print("Maximum of objective function: %i" % solver.ObjectiveValue())
print()
print("x value: ", solver.Value(x))
print("y value: ", solver.Value(y))
print("z value: ", solver.Value(z))
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
# 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())
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)}")
solution_printer = VarArraySolutionPrinter([S,E,N,D,M,O,R,Y])
status = solver.SearchForAllSolutions(model, solution_printer)
print("Status = %s" % solver.StatusName(status))
print("Number of solutions found: %i" % solution_printer.solution_count())
練習問題：覆面算
https://www.math.uni-bielefeld.de/~sillke/PUZZLES/ALPHAMETIC/ にある覆面算の問題をどれか一つ解いてみよ．
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)
status = solver.Solve(model)
for i in range(n):
for j in range(n):
print(solver.Value(x[i,j]), end=" ")
print()
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])
status = solver.Solve(model)
for i in line:
for j in line:
print(solver.Value(x[i,j]), end=" ")
print()
solution_printer = VarArraySolutionPrinter([x[i,j] for i in line for j in line])
status = solver.SearchForAllSolutions(model, solution_printer)
print("Status = %s" % solver.StatusName(status))
print("Number of solutions found: %i" % solution_printer.solution_count())
練習問題：不等式パズル
$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]
]