第11章整数計画法で数独を解く¶

整数計画法では、数独の候補を0-1変数、規則を線形等式で表します。たとえば左上のマスに5を置く候補は、置くなら1、置かないなら0です。同じマスについて1から9まで九つの変数を用意し、その合計を 1にすれば、どれか一つの数字だけが選ばれます。

整数計画法は、変数の一部または全部を整数に制限し、線形の制約を満たす解を探す方法です。通常は費用を最小にする、利益を最大にするといった目的関数も与えます。数独では目的値で盤面を選び分けず、制約をすべて満たす実行可能解を探します。

この章では、線形計画法と混合整数計画法のソルバー HiGHS をPython APIの highspy から使います。HiGHSはLP、MIP、QPを扱うオープンソースの最適化ソフトウェアです。[1] この章で使うのは、0-1変数を含む混合整数計画（MIP）の機能です。

二つの品物から選ぶ¶

0-1整数計画の小さな例として、品物AとBを袋へ入れるか決めます。品物Aは重さ2、価値3、品物Bは重さ3、価値4です。袋には重さ3まで入ります。二値変数 \(a,b\) を使うと、モデルは次のように書けます。

\[\begin{split}\begin{aligned} \text{maximize}\quad & 3a+4b \\ \text{subject to}\quad & 2a+3b \leq 3 \\ & a,b\in\{0,1\} \end{aligned}\end{split}\]

制約は袋の重さを表し、目的関数は価値です。AとBを同時には選べないので、価値4のBを選ぶのが最適です。highspy のモデリングAPIでは、二値変数、線形式、目的関数を式に近い形で書けます。 [2]

model = highspy.Highs()
model.setOptionValue("output_flag", False)

# 品物を選ぶなら1、選ばないなら0になる二値変数です。
item_a = model.addBinary(name="item_a")
item_b = model.addBinary(name="item_b")

# 重さの合計を3以下にし、価値の合計を最大化します。
model.addConstr(2 * item_a + 3 * item_b <= 3)
model.maximize(3 * item_a + 4 * item_b)

status = model.getModelStatus()
print(f"status: {model.modelStatusToString(status)}")
print(f"item_a: {round(model.val(item_a))}")
print(f"item_b: {round(model.val(item_b))}")
print(f"objective: {model.getObjectiveValue():.0f}")

status: Optimal
item_a: 0
item_b: 1
objective: 4

Optimal は、制約を満たす組合せが見つかり、それより良い目的値がないと確認できた状態です。数独でも同じ状態名が返りますが、目的関数の使い方が異なります。

次の表に出てくる LP緩和 は、0-1条件をいったん0以上1以下の実数まで広げた問題です。 分枝限定法 は、分数になった変数を0と1の側へ分け、不要な枝を除きながら整数解を探します。詳しい動きは後の補足で説明します。

人間とソルバーの役割分担¶
段階	人間（Pythonコード）が行うこと	ソルバー（HiGHS）が行うこと
入力	729個の0-1変数を作り、324本の線形等式と目的関数（定数0）を設定する	―
探索	―	LP緩和で上界を求め、分枝限定法（branch-and-cut）で整数解を探索する
出力	各変数の値を0.5で閾値判定し、盤面を復元・検証する	`Optimal`（変数値付き）または `Infeasible` を返す

729個の二値変数を作る¶

行 \(r\)、列 \(c\) のマスへ数字 \(d\) を置くかどうかを、二値変数 \(x_{r,c,d}\) で表します。

\[\begin{split}x_{r,c,d}=\begin{cases} 1 & \text{マス }(r,c)\text{ に }d\text{ を置くとき} \\ 0 & \text{それ以外} \end{cases}\end{split}\]

添字はそれぞれ1から9までなので、変数は \(9^3=729\) 個です。一つのマスを1〜9の整数変数一つで表す方法もありますが、整数計画法では線形等式を書きやすい二値変数を使います。

たとえば左上のマスについては、九候補を足した \(\sum_{d=1}^{9}x_{1,1,d}=1\) を置きます。これは「九候補のうち一つだけが1になる」という式です。 Pythonコードでは同じ候補を x[0, 0, digit] と書き、digit を0から8まで動かします。数式の添字は1から、Pythonの添字は0から始まる点だけが異なります。

数独の四種類の条件は、係数がすべて1の線形等式になります。

各マスに数字を一つ置く: マス \((r,c)\) ごとに \(\sum_{d=1}^{9}x_{r,c,d}=1\) とします。
各行で数字を一度ずつ使う: 行 \(r\) と数字 \(d\) ごとに \(\sum_{c=1}^{9}x_{r,c,d}=1\) とします。
各列で数字を一度ずつ使う: 列 \(c\) と数字 \(d\) ごとに \(\sum_{r=1}^{9}x_{r,c,d}=1\) とします。
各ブロックで数字を一度ずつ使う: ブロック \(B\) と数字 \(d\) ごとに \(\sum_{(r,c)\in B}x_{r,c,d}=1\) とします。

各種類に81本ずつあり、合わせて324本の等式です。初期配置が \(d\) のマスには、さらに \(x_{r,c,d}=1\) を加えます。

def build_model(board: Board):
    """数独を0-1整数計画モデルへ変換します。"""
    side = isqrt(len(board))
    box = isqrt(side)
    model = highspy.Highs()
    model.setOptionValue("output_flag", False)

    # x[r,c,d]は、マス(r,c)へ数字dを置くときだけ1になります。
    x = {
        (row, col, digit): model.addBinary(name=f"x_{row + 1}_{col + 1}_{digit + 1}")
        for row in range(side)
        for col in range(side)
        for digit in range(side)
    }

    # 各マスでは、候補となる数字をちょうど一つ選びます。
    for row in range(side):
        for col in range(side):
            model.addConstr(sum(x[row, col, digit] for digit in range(side)) == 1)

    # 各行と各列では、それぞれの数字をちょうど一回使います。
    for digit in range(side):
        for row in range(side):
            model.addConstr(sum(x[row, col, digit] for col in range(side)) == 1)
        for col in range(side):
            model.addConstr(sum(x[row, col, digit] for row in range(side)) == 1)

    # 各ブロックでも、それぞれの数字をちょうど一回使います。
    for digit in range(side):
        for box_row in range(0, side, box):
            for box_col in range(0, side, box):
                model.addConstr(
                    sum(
                        x[row, col, digit]
                        for row in range(box_row, box_row + box)
                        for col in range(box_col, box_col + box)
                    )
                    == 1
                )

    # 初期配置に対応する二値変数は1に固定します。
    for position, given in enumerate(board):
        if given:
            row, col = divmod(position, side)
            model.addConstr(x[row, col, given - 1] == 1)

    # 目的係数はすべて0です。今回は制約を満たす点だけを探します。
    return model, x

model.addBinary() が0-1変数を作り、model.addConstr() が線形等式を追加します。コードは盤面の一辺 side を使っているので、同じモデルで4×4の数独も表せます。

目的関数を0にする¶

小さな例では価値を最大化しました。数独では、どの完成盤面が得られてもよいため、729個の目的係数をすべて0にしています。すべての実行可能解の目的値は0です。HiGHSが一つ見つけ、ほかにより良い目的値がないと確認すると、モデル状態は Optimal になります。

\[\text{minimize}\quad 0\]

掲載コードでは addBinary() で作る全変数の目的係数を0のままにしており、この定数の目的関数に対応します。

ここで Optimal が表すのは、制約を満たす盤面が見つかり、目的値0について最適性の判定が終わったことです。実行可能性問題を、定数の目的関数を持つ最適化問題としてソルバーへ渡しています。

補足：LP緩和と分枝限定法

二値変数の条件を外し、\(0\leq x\leq1\) の実数を許した問題を LP緩和 と呼びます。品物の例を緩和すると、\(a=1,b=1/3\) も制約を満たし、目的値は \(13/3\) です。整数解の最適値4より大きいので、緩和問題の値は整数問題の上界として使えます。

分数になった変数について0の側と1の側へ問題を分け、各側の上界を調べる方法が 分枝限定法 です。実際のMIPソルバーは、不要な枝を除く前処理や不等式も組み合わせます。 HiGHSは分枝限定にカットを加えたbranch-and-cutソルバーを使います。[3]

解を盤面へ戻す¶

Optimal のとき、各二値変数の値を model.vals() で読みます。数値計算上の許容誤差があるため、コードでは0.5より大きい変数を1として扱います。各マスで一つだけ選ばれていることを確認して数字へ戻し、共通検証器で初期配置、行、列、ブロックをもう一度検査します。

一意性を調べるには、最初に得た盤面だけを除外します。完成盤面で1だった81個の変数を集めた集合を \(S\) とすると、追加する制約は次のとおりです。

\[\sum_{x\in S}x\leq80\]

最初の盤面では左辺が81になるので、この制約を満たしません。別の完成盤面では少なくとも一マスの数字が変わり、最初に選ばれた変数の一つ以上が0になります。このように、既知の整数解を除く制約を no-good制約 と呼びます。

def solve(board: Board, limit: int = 1) -> SolveResult:
    side = isqrt(len(board))
    model, x = build_model(board)
    solutions: list[Board] = []

    while len(solutions) < limit:
        model.run()
        model_status = model.getModelStatus()

        if model_status == highspy.HighsModelStatus.kOptimal:
            # 目的値0の最適解は、数独の全制約を満たす実行可能解です。
            solution = board_from_solution(model, x, side)
            valid, errors = validate_solution(solution, board)
            if not valid:
                raise RuntimeError("invalid HiGHS solution: " + "; ".join(errors))
            solutions.append(solution)

            # 完成盤面で選ばれた変数のうち、次の解では少なくとも一つを0にします。
            selected = [
                x[row, col, solution[row * side + col] - 1]
                for row in range(side)
                for col in range(side)
            ]
            model.addConstr(sum(selected) <= side * side - 1)

        elif model_status == highspy.HighsModelStatus.kInfeasible:
            status = "solved" if solutions else "unsat"
            return SolveResult(status, tuple(solutions), search_complete=True)
        else:
            # 制限到達や中断を、実行不能と取り違えないようにします。
            status = "solved" if solutions else "unknown"
            reason = model.modelStatusToString(model_status)
            return SolveResult(status, tuple(solutions), False, reason)

    return SolveResult("solved", tuple(solutions), search_complete=False)

no-good制約追加後のモデルが Infeasible と判定された状態は、初期解以外の解空間が空集合であることを意味します。Optimal が返された状態は二つ目の解が存在することを示します。時間制限や中断等で最終判定に到達しなかった状態は unknown として扱い、unsat とは明確に区別します。 HiGHSのモデル状態には最適、実行不能、制限到達などが独立して定義されています。[4]

実行と検証¶

実行手順とオプションは examples/08-integer-programming/README.md にまとめています。最初の解を得た後にno-good制約を加え、最大二解まで探した結果は次のとおりです。

HiGHSで三つの盤面を調べた結果¶
問題	`status`	解数	`search_complete`	`unique`
一意解	`solved`	1	`yes`	`yes`
矛盾	`unsat`	0	`yes`	対象外
複数解	`solved`	2	`no`	`no`

一意解問題では、最初の解を除外したモデルが Infeasible になり、二つ目の解がないと分かりました。矛盾問題は最初のモデルから Infeasible です。複数解問題ではno-good制約を加えた後も Optimal が返り、異なる二盤面が共通検証器を通りました。二解で探索を止めているため search_complete は no ですが、複数解であることは確定します。

解がない場合と、解が複数ある場合¶

矛盾問題では、初回のモデル構築・解探索の時点で実行不能（Infeasible）となります。初期配置を保ったまま、324本の線形等式を満たす0-1変数の割り当てが存在しないためです。

複数解問題では、no-good制約を加えた後もソルバーが再び最適解（Optimal）を返し、異なる二つの完成盤面が得られることから、一意解ではない（複数解を持つ）ことが判定されます。

この方法で分かること¶

数独の完成盤面からは、各マスで対応する変数を1にすれば線形等式をすべて満たす割り当てを作れます。反対に、モデルの実行可能解では各マス、行、列、ブロックの等式が成立するため、選ばれた変数を数独の完成盤面として読めます。この双方向の対応により、0-1整数計画への変換は数独の解を増減させません。

HiGHSが Optimal を返した場合は実行可能解を取得し、Infeasible を返した場合は解なしと判定します。最初の解を除外して Infeasible となった場合は一意解、二つ目の実行可能解が得られた場合は複数解と判定します。制限到達などの状態では結論を出さず、unknown として処理します。

本実装は決定論的に動作します。同一のモデルであっても、HiGHSのバージョンや内部設定により複数解の中で最初に得られる盤面が異なることがあります。前章までに扱ったSAT、SMT、制約プログラミングとの処理速度の比較は、変数数や一単一問題の結果のみから決定することはできません。

参考文献¶

[1]

HiGHS Development Team. HiGHS: High-Performance Parallel Linear Optimization Software. URL: https://highs.dev/ (visited on 2026-08-08).

[2]

HiGHS Development Team. HiGHS Python Interface: Modelling. URL: https://ergo-code.github.io/HiGHS/dev/interfaces/python/model-py/ (visited on 2026-08-08).

[3]

HiGHS Development Team. HiGHS Solvers. URL: https://ergo-code.github.io/HiGHS/dev/solvers/ (visited on 2026-08-08).

[4]

HiGHS Development Team. HiGHS Enums: HighsModelStatus. URL: https://ergo-code.github.io/HiGHS/dev/structures/enums/ (visited on 2026-08-08).

第11章 整数計画法で数独を解く¶