CP-SAT（Google OR-Tools）入門

1.概要

1.1 CP-SATとは

CP-SATは、Googleが開発したオープンソースの最適化計算ツール「OR-Tools」に含まれているエンジンです。プログラム上で条件を柔軟に書ける「制約プログラミング（CP）」の表現力と、論理パズルのような問題を高速で解く「SATソルバー」の探索技術が組み合わさっています。

1.2 MIP（混合整数最適化） vs メタヒューリスティクス vs CP-SAT

解きたい問題の性質に合わせて、計算手法を選ぶ必要があります。

特徴	MIP（例: Gurobi, CBC, SCIP）	メタヒューリスティクス（GA, 局所探索）	CP-SAT
得意な問題	コスト最小化、配送計画、連続変数が混ざる問題	超大規模（数万変数）の近似解探索	複雑な論理条件・時間枠・順序・離散組み合わせ問題
条件の表現	すべてを「線形不等式」で表す必要あり	目的関数に評価式を書けば何でも動く	不等式に加え、AND/OR/IF文、全体制約を直接表現可能
最適性の保証	あり（下界・上界で証明可能）	なし（良い解は見つかるが最善かは不明）	あり（最適解であること、あるいは解なしを証明可能）
小数の扱い	連続変数・小数係数を自然に扱える	制限なし	整数のみ（小数は整数へ変換する必要あり）

1.3 制約伝播（Constraint Propagation）の直感的理解

CP-SATが計算を速く終わらせることができるのは、すべてのパターンを総当たりする前に「制約伝播」という仕組みで、あり得ない選択肢（ドメイン）を一気に削り落とすからです。

• 例：$x$と$y$がどちらも1から9の整数で、$x + y = 3$という条件がある場合、計算を始める前に「$x$も$y$も1か2のどちらか」以外はあり得ないと自動的に候補を絞り込みます。

1.4 インストール

pip install ortoolsまたはuv add ortoolsでインストールしてください。

以降では、Python 3.14、ortools 9.15.6755で動作確認しています。

2. 変数型とドメインの扱い

CP-SATの世界では、すべての事象を以下の3つの変数型で表現します。

1. IntVar（整数変数）
2. BoolVar（真偽値変数 / 0-1変数）
3. IntervalVar（区間変数）

2.1 IntVar（整数変数）

個数や金額、時刻など、とびとびの数（離散的な値）を表します。

from ortools.sat.python import cp_model

model = cp_model.CpModel()

# 書式: model.new_int_var(下限, 上限, "名前")
x = model.new_int_var(0, 100, "x")

# 離散的な飛び飛びの値を指定したい場合（例: 1, 3, 5〜10）
domain = cp_model.Domain.from_values([1, 3])
domain = domain.union_with(cp_model.Domain(5, 10))
y = model.new_int_var_from_domain(domain, "y")

2.2 BoolVar（真偽値変数 / 0-1変数）

「採用するか（1）、しないか（0）」や「あるイベントが起きたかどうか」を表す、最もよく使う変数です。

# 内部的には new_int_var(0, 1, "is_selected") と同等
is_selected = model.new_bool_var("is_selected")

# 否定（NOT）の取得: 1 - is_selected と同じ意味
not_selected = ~is_selected

2.3 IntervalVar（区間変数）

開始時刻（start）、所要時間（size/duration）、終了時刻（end）を1セットにしたスケジューリング専用変数です。常に $\text{start} + \text{size} = \text{end}$ が保たれます。

start = model.new_int_var(0, 100, "start")
size = 8  # 固定値でも IntVar でも可
end = model.new_int_var(0, 100, "end")

# 区間変数の作成
task_interval = model.new_interval_var(
    start, size, end, "task_interval"
)

【重要】小数の扱い方について
CP-SATでは model.add(x * 1.5 <= y) のような小数は書けません。
小数第2位まで扱う場合は、すべて100倍して整数化します。

• 例: $1.5x \le y \implies 150x \le 100y$ （または $3x \le 2y$）

2.4 Optional IntervalVar（オプショナル区間変数）

実行するかどうか（採用されるかどうか）が確定していないタスクのスケジューリングには、オプショナル区間変数を使用します。

この変数は、指定した真偽値変数（BoolVar）が 1 のときは通常の区間変数として機能し、0 のときは「開始時刻＝終了時刻（所要時間0）、かつ他のタスクと重複しない」ものとして無視されるという非常に便利な性質を持っています。

start = model.new_int_var(0, 100, "start")
size = 8
end = model.new_int_var(0, 100, "end")
is_performed = model.new_bool_var("is_performed") # 実行するかどうか

# オプショナル区間変数の作成
optional_interval = model.new_optional_interval_var(
    start, size, end, is_performed, "optional_interval"
)

3. よく使うモデリングパターン

現実の条件をCP-SATのプログラムに翻訳するためのパターン集です。

3.1 論理結合（AND / OR / XOR）

a = model.new_bool_var("a")
b = model.new_bool_var("b")
target = model.new_bool_var("target")

# 1. target <=> a AND b
model.add_bool_and([a, b]).only_enforce_if(target)
model.add_bool_or([~a, ~b]).only_enforce_if(~target)

# 2. target <=> a OR b
model.add_bool_or([a, b]).only_enforce_if(target)
model.add_bool_and([~a, ~b]).only_enforce_if(~target)

# 3. 排他的論理和（XOR: どちらか一方のみ1）
model.add(a + b == 1)

3.2 条件付き制約（only_enforce_if）

「もし〜ならば」という条件を表すときによく使います。

flag = model.new_bool_var("flag")
work_hours = model.new_int_var(0, 24, "work_hours")

# flag が 1 ならば work_hours >= 8
model.add(work_hours >= 8).only_enforce_if(flag)

# flag が 0 ならば work_hours == 0
model.add(work_hours == 0).only_enforce_if(~flag)

3.3 最大・最小・絶対値の等価制約

最大値などを求める関数も、直接書くことができます。

x = model.new_int_var(-10, 10, "x")
y = model.new_int_var(0, 10, "y")
max_val = model.new_int_var(0, 20, "max_val")
abs_x = model.new_int_var(0, 10, "abs_x")

# max_val = max(x, y, 5)
model.add_max_equality(max_val, [x, y, 5])

# abs_x = |x|
model.add_abs_equality(abs_x, x)

3.4 グローバル制約（重複禁止・要素参照）

ソルバーの内部アルゴリズムと直結しているため、不等式を並べるより遥かに高速です。

vars_list = [model.new_int_var(1, 9, f"v{i}") for i in range(5)]

# 1. 全て異なる値をとる（数独、順列）
model.add_all_different(vars_list)

# 2. 変数インデックスによる配列アクセス（Element制約）
# target == costs[selected_index]
costs = [100, 250, 400, 700]
selected_index = model.new_int_var(0, 3, "idx")
target_cost = model.new_int_var(0, 1000, "cost")
model.add_element(selected_index, costs, target_cost)

3.5 スケジューリング向けのグローバル制約（Cumulative）

スケジューリングにおいて、「NoOverlap（重複禁止）」の次に重要になるのが「Cumulative（累積制約）」です。これは、各タスクがリソース（人員、電力、機械のキャパシティなど）を消費する場合に、どの瞬間においてもリソース消費の合計が上限を超えないことを保証する強力な制約です。

# リソースの上限（例：同時に作業できる人員は最大5人）
capacity = 5

intervals = []
demands = []

# タスク1: 所要時間3, 人員2名を消費
s1 = model.new_int_var(0, 10, "s1")
e1 = model.new_int_var(0, 10, "e1")
intervals.append(model.new_interval_var(s1, 3, e1, "task1"))
demands.append(2)

# タスク2: 所要時間4, 人員3名を消費
s2 = model.new_int_var(0, 10, "s2")
e2 = model.new_int_var(0, 10, "e2")
intervals.append(model.new_interval_var(s2, 4, e2, "task2"))
demands.append(3)

# タスク3: 所要時間2, 人員2名を消費
s3 = model.new_int_var(0, 10, "s3")
e3 = model.new_int_var(0, 10, "e3")
intervals.append(model.new_interval_var(s3, 2, e3, "task3"))
demands.append(2)

# どの時刻でも、実行中のタスクの demands の合計が capacity(5) 以下になるようにする
model.add_cumulative(intervals, demands, capacity)

4. ソフト制約と目的関数

問題を解くときには、「違反しても計算はできるが、なるべく避けたい条件（ソフト制約）」がたくさん出てきます。

4.1 ハード制約とソフト制約の分離

• ハード制約: 破綻を許さない物理・法的条件（例: 定員超過の禁止） $\rightarrow$ model.add()
• ソフト制約: 努力目標・希望条件 $\rightarrow$ 違反量を変数化し目的関数でペナルティ化

4.2 スラック変数（緩衝変数）によるペナルティモデリング

「できれば40時間以内にしたいが、超えたら1時間あたり1000のペナルティをつける」という場合の書き方です。

actual_hours = model.new_int_var(0, 60, "hours")
overtime = model.new_int_var(0, 20, "overtime")

# overtime >= actual_hours - 40 かつ overtime >= 0
model.add_max_equality(overtime, [actual_hours - 40, 0])

# コスト最小化
penalty_cost = overtime * 1000
model.minimize(penalty_cost)

4.3 辞書式多目的最適化（Lexicographical Optimization）

「第1目標を絶対に優先し、その次に第2目標をよくしたい」といったように優先順位がはっきりしている場合、1つの式にまとめてしまうとうまくいきません。この場合は2段階に分けて計算するのが一般的です。

unassigned_count = model.new_int_var(0, 100, "unassigned_count")
total_cost = model.new_int_var(0, 100, "total_cost")

# --- フェーズ1: 第1目標（未割当数）の最小化 ---
model.minimize(unassigned_count)
solver = cp_model.CpSolver()
solver.solve(model)
best_unassigned = round(solver.objective_value)

# --- フェーズ2: 第1目標を制約として固定し、第2目標（総コスト）を最小化 ---
model.add(unassigned_count == best_unassigned)
model.minimize(total_cost)
solver.solve(model)

5. 実践

5.1 ビンパッキング問題（トラック積載最適化）

荷物をできるだけ少ない台数のトラックに詰め込む問題です。

from ortools.sat.python import cp_model

weights = [48, 30, 19, 36, 36, 27, 42, 18, 24, 30]
bin_capacity = 100
num_items = len(weights)
num_bins = len(weights)  # 最悪でも荷物と同じ台数

model = cp_model.CpModel()

# x[i, j] = 荷物 i が ビン j に入るなら 1
x = {}
for i in range(num_items):
    for j in range(num_bins):
        x[i, j] = model.new_bool_var(f"x_{i}_{j}")

# y[j] = ビン j が使用されるなら 1
y = [model.new_bool_var(f"y_{j}") for j in range(num_bins)]

# 制約1: 各荷物は必ず1つのビンに入る
for i in range(num_items):
    model.add(sum(x[i, j] for j in range(num_bins)) == 1)

# 制約2: ビンの容量制約（ビンが使われていなければ上限0）
for j in range(num_bins):
    model.add(
        sum(
            weights[i] * x[i, j]
            for i in range(num_items)
        ) <= bin_capacity * y[j]
    )

# 目的関数: 使用するビン数の最小化
model.minimize(sum(y))

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

if status in {cp_model.OPTIMAL, cp_model.FEASIBLE}:
    print(f"最小必要ビン数: {round(solver.objective_value)}")
    for j in range(num_bins):
        if solver.value(y[j]) == 1:
            items = [
                i for i in range(num_items)
                if solver.value(x[i, j]) == 1
            ]
            total_w = sum(weights[i] for i in items)
            print(f"ビン {j}: 荷物 {items}")
            print(f"  合計重量: {total_w}/{bin_capacity}")

5.2 ナーススケジューリング（実務レベルのシフト編成）

「連続勤務の制限」や「夜勤の翌日は休み」といった、日付をまたぐ条件をプログラムにします。

from ortools.sat.python import cp_model

model = cp_model.CpModel()

num_nurses = 4
num_days = 8
shifts = [0, 1, 2]  # 0: 休, 1: 日勤, 2: 夜勤

# 変数定義: shift[(n, d, s)]
work = {}
for n in range(num_nurses):
    for d in range(num_days):
        for s in shifts:
            work[(n, d, s)] = model.new_bool_var(f"w_{n}_{d}_{s}")

# 1. 各看護師は1日1シフト
for n in range(num_nurses):
    for d in range(num_days):
        model.add(sum(work[(n, d, s)] for s in shifts) == 1)

# 2. 各シフトの必要人数（日勤:2人, 夜勤:1人）
for d in range(num_days):
    model.add(sum(work[(n, d, 1)] for n in range(num_nurses)) == 2)
    model.add(sum(work[(n, d, 2)] for n in range(num_nurses)) == 1)

# 3. 夜勤の翌日は必ず「休み」(夜勤の翌日に日勤や夜勤は不可)
for n in range(num_nurses):
    for d in range(num_days - 1):
        # d日目が夜勤(2)ならば、d+1日目は休み(0)
        model.add(
            work[(n, d + 1, 0)] == 1
        ).only_enforce_if(work[(n, d, 2)])

# 4. 4連勤以上を禁止（どの連続する4日間でも、休みの合計が1日以上）
for n in range(num_nurses):
    for d in range(num_days - 3):
        model.add(sum(work[(n, d + i, 0)] for i in range(4)) >= 1)

# 充足可能性を解く
solver = cp_model.CpSolver()
status = solver.solve(model)

if status in {cp_model.OPTIMAL, cp_model.FEASIBLE}:
    symbols = {0: " 休 ", 1: " 日 ", 2: " 夜 "}
    print("\n=== シフトスケジュール結果 ===")
    print("         " + " ".join(
        [f" D{d+1}" for d in range(num_days)]
    ))
    for n in range(num_nurses):
        row = [
            symbols[s]
            for d in range(num_days)
            for s in shifts
            if solver.value(work[(n, d, s)]) == 1
        ]
        print(f"看護師{n}: " + "".join(row))

5.3 ジョブショップスケジューリング（JSP）

複数の製品を、決まった手順で機械に通して作る際の日程計画です。add_no_overlap という機能が活躍します。

from ortools.sat.python import cp_model

model = cp_model.CpModel()

# ジョブ定義: 各タスク (使用する機械ID, 処理時間)
# ジョブ0: 機械0(3h) -> 機械1(2h) -> 機械2(2h)
jobs_data = [
    [(0, 3), (1, 2), (2, 2)],
    [(0, 2), (2, 1), (1, 4)],
    [(1, 4), (2, 3)]
]

num_machines = 3
all_machines = range(num_machines)
horizon = sum(task[1] for job in jobs_data for task in job)

all_tasks = {}
machine_to_intervals = {m: [] for m in all_machines}

for j, job in enumerate(jobs_data):
    for t, (machine, duration) in enumerate(job):
        start = model.new_int_var(0, horizon, f"start_{j}_{t}")
        end = model.new_int_var(0, horizon, f"end_{j}_{t}")
        interval = model.new_interval_var(
            start, duration, end, f"interval_{j}_{t}"
        )
        all_tasks[(j, t)] = (start, end, interval)
        machine_to_intervals[machine].append(interval)

# 制約1: 同一ジョブ内の工程順序（先行タスクのend <= 後続タスクのstart）
for j, job in enumerate(jobs_data):
    for t in range(len(job) - 1):
        model.add(all_tasks[(j, t + 1)][0] >= all_tasks[(j, t)][1])

# 制約2: 機械の重複稼働禁止（NoOverlap）
for machine in all_machines:
    model.add_no_overlap(machine_to_intervals[machine])

# 目的: メイクスパン（全作業完了時刻）の最小化
makespan = model.new_int_var(0, horizon, "makespan")
model.add_max_equality(
    makespan,
    [
        all_tasks[(j, len(job) - 1)][1]
        for j, job in enumerate(jobs_data)
    ]
)
model.minimize(makespan)

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

if status in {cp_model.OPTIMAL, cp_model.FEASIBLE}:
    print(f"\n総完了時間 (Makespan): {round(solver.objective_value)}")
    for j, job in enumerate(jobs_data):
        for t, (machine, duration) in enumerate(job):
            st = solver.value(all_tasks[(j, t)][0])
            en = solver.value(all_tasks[(j, t)][1])
            print(f"Job {j} Task {t} (機械{machine}): {st} -> {en}")

6. デバッグ・性能チューニング

プログラムを書く上で一番の壁になる「答えが見つからない（INFEASIBLE）」「計算が終わらない」ときの対処法です。

6.1 「INFEASIBLE（解なし）」の体系的デバッグ法

条件を厳しくしすぎて答えが存在しなくなった場合、CP-SATはどこが悪いのか直接は教えてくれません。以下の3つの方法で原因を探します。

1. 二分探索法（コメントアウト法）: 追加した条件を半分ずつコメントアウトして無効にし、どの条件のせいで答えが消えるのかを突き止めます。

2. ソフト制約化:
   怪しい条件にペナルティ用の変数をつけ、「ペナルティ最小化」で計算させます。

   # model.add(x + y <= 10) の代わりにスラック変数を追加
   slack = model.new_int_var(0, 1000, "slack_violation")
   model.add(x + y <= 10 + slack)
   model.minimize(slack)
   # 解いた結果、slack > 0 になっている箇所が矛盾の原因
   

3. Assumptions（仮定）機能による診断:
   OR-Toolsには、どの制約が矛盾しているかをソルバーに抽出させる機能（sufficient_assumptions_for_infeasibility）があります。

   # 制約ごとに診断用フラグを用意
   c1 = model.new_bool_var("c1")
   c2 = model.new_bool_var("c2")
   
   model.add(x >= 5).only_enforce_if(c1)
   model.add(x <= 3).only_enforce_if(c2)
   
   # 両方の制約を前提（仮定）としてセット
   model.add_assumptions([c1, c2])
   
   solver = cp_model.CpSolver()
   status = solver.solve(model)
   
   if status == cp_model.INFEASIBLE:
       # 互いに矛盾している仮定のリスト（極小不整合部分集合: MUS）を出力
       conflict = solver.sufficient_assumptions_for_infeasibility()
       print("矛盾している制約フラグ:", conflict)
   

6.2 探索ログ（Log Search Progress）の完全解読

ソルバーが裏でどう計算しているかを把握するために、ログを出力させます。

solver.parameters.log_search_progress = True

【出力ログの代表例】

#1       0.01s best:32   next:[] 
#2       0.03s best:25   next:[] 
*        0.05s best:18   best_bound:12  no_lp:[...]
===== Solutions =====
Status: OPTIMAL
Objective: 18
Best bound: 18
Booleans: 124
Branches: 820

• best: 現在見つかっているとりあえずの答えのスコア。計算が進むにつれて良くなっていきます。
• best_bound: 数学的に計算された「理論上の限界値」。本当の最適解であっても、この値より良くなることは絶対にありません。
• ギャップの消滅: best == best_bound になった瞬間に、ソルバーは「最適解であることが100%証明された」と判断して計算を終えます。
• Branches: パターンを分岐した回数。この回数が増え続けているのにbestやboundが変わらない場合は、変数の範囲が広すぎることなどが疑われます。

6.3 大規模問題のチューニング・パラメータ設計

複雑な問題を解くときのおすすめ設定です。

solver = cp_model.CpSolver()

# 1. タイムアウト設定（必須）
solver.parameters.max_time_in_seconds = 60.0

# 2. 並列ワーカー数（CPUコア数に合わせる）
# CP-SATはワーカーごとに異なる探索戦略
# （LNS, LP緩和, SAT探索など）を分散協調実行します
solver.parameters.num_search_workers = 8

# 3. 許容ギャップ（早期終了）
# 「理論限界値から 2% 以内の誤差なら十分」として計算を打ち切る
solver.parameters.relative_gap_limit = 0.02

6.4 高速化のためのプログラムのコツ

1. 対称性の排除（同じ意味のパターンの省略）

• 例: 「同じ性能のトラック3台」がある場合、トラック1, 2, 3 の順番を入れ替えただけの同じ意味の答えが6通り存在してしまい、無駄な計算が増えます。
• 対策: トラックを使う順番を固定する条件（トラック0の使用 >= トラック1の使用 >= トラック2の使用）を追加して、無駄な探索を省きます。

1. 変数の範囲（ドメイン）は可能な限り小さく設定する

• 最大値がある程度わかっているなら new_int_var(0, 1000000) のように大きく書かず、理論上の上限を計算してできるだけ狭い範囲を指定してください。

7. 演習

7.1 整数最適化

• 変数 $x, y$ はそれぞれ $0$ 以上 $10$ 以下の整数
• 制約 1: $x + 2y \le 14$
• 制約 2: $3x - y \ge 0$
• 目的: $2x + y$ を最大化する

from ortools.sat.python import cp_model

7.2 ナップサック問題

重さ制限 $15\text{kg}$ のリュックに、価値の合計が最大になるように荷物を詰める。

荷物	価値 (円)	重さ (kg)
A	1200	3
B	1500	5
C	2000	8
D	3000	10

from ortools.sat.python import cp_model

7.3 解答例

整数最適化

from ortools.sat.python import cp_model

model = cp_model.CpModel()
x = model.new_int_var(0, 10, "x")
y = model.new_int_var(0, 10, "y")
model.add(x + 2 * y <= 14)
model.add(3 * x - y >= 0)
model.maximize(2 * x + y)
solver = cp_model.CpSolver()
status = solver.solve(model)
if status in {cp_model.OPTIMAL, cp_model.FEASIBLE}:
    print("最適解:", solver.objective_value)
    print("x =", solver.value(x), "y =", solver.value(y))

ナップサック問題

from ortools.sat.python import cp_model

model = cp_model.CpModel()
x = [model.new_bool_var(c) for c in ["A", "B", "C", "D"]]

model.add(sum(v * c for v, c in zip(x, [3, 5, 8, 10], strict=True)) <= 15)
model.maximize(sum(v * c for v, c in zip(x, [1200, 1500, 2000, 3000], strict=True)))
solver = cp_model.CpSolver()
status = solver.solve(model)
if status in {cp_model.OPTIMAL, cp_model.FEASIBLE}:
    print("最適解:", solver.objective_value)
    print("x =", [solver.value(v) for v in x])