この記事では、Python の強力なライブラリである PuLP を使用して、線形プログラミングで数独の問題を解く方法についてご紹介してきます。

PuLP とは

PuLP は、線形計画問題（Linear Programming, LP）や混合整数線形計画問題（Mixed-Integer Linear Programming, MILP）など、様々な最適化問題を解くための Python ライブラリです。主に最適化問題のモデリングと解決のために設計されており、Python で書かれた使いやすいインターフェースを提供します。

PuLP の主な特徴は以下の通りです。

使いやすい: PuLP は Python 言語を使って直感的に数理最適化問題を記述できるように設計されており、Python から簡単に利用することができます。
柔軟性: 線形プログラミング、整数線形プログラミング、混合整数線形プログラミングなど、さまざまな種類の最適化問題を扱うことができます。
多様なソルバーのサポート: PuLP は、CBC、GLPK、CPLEX、Gurobi など、さまざまなソルバーと連携して動作します。これにより、ユーザは様々なソルバーを試すことができ、特定の問題に最適なソルバーを選択できるという柔軟性があります。
- CBC (Coin-or branch and cut): CBC は、整数プログラミング問題を解くためのオープンソースソルバーで、ブランチアンドカットアルゴリズムを使用しており、主に大規模な整数線形問題に適しています。
- GLPK (GNU Linear Programming Kit): GLPK は、線形プログラミング（LP）や混合整数プログラミング（MIP）の問題を解くための無料のオープンソースソルバーで、シンプレックス法とブランチアンドバウンドアルゴリズムを提供しています。主に、線形プログラミング問題や混合整数問題向けに設計されています。
- CPLEX: CPLEX は、IBM によって開発された商用製品で、高性能な数学的最適化ソルバーです。線形プログラミング、混合整数プログラミング、二次プログラミングなど幅広い問題を解くことができます。高度なアルゴリズムと平行処理能力を持ち合わせているため、特に複雑で大規模な最適化問題に適しています。
- Gurobi: Gurobi は、業界最高水準の性能を誇る商用の最適化ソルバーで、線形プログラミング、混合整数プログラミング、二次プログラミングなどの問題を高速に解くことができます。シンプレックス法、内点法、ブランチアンドカットなどのアルゴリズムを使用しています。
オープンソース: PuLP はオープンソースとなっており、無料で利用、改変、配布することができます。

PuLP は主に、リソース割り当て、スケジューリング、生産計画、輸送問題などの最適化問題を解決するために使用されます。数独のようなパズルを解くためにも使用することができ、その柔軟性とアクセシビリティから、商用利用や学術利用の両方で広く使われています。

数独とは

数独は、9x9 のグリッドに数字を埋めるパズルゲームです。各行、各列、そして９つの 3x3 のブロックそれぞれに 1 から 9 までの数字が重複なく配置される必要があります。

数独の主な制約は以下の通りです。

行の制約：各行に1から9までの数字が一度だけ現れる。
列の制約：各列に1から9までの数字が一度だけ現れる。
ブロックの制約：各3x3のブロックに1から9までの数字が一度だけ現れる。

これらの制約を満たす解が、数独の正しい解となります。

実装

それでは、Python を使って、PuLP を使って数独パズルをモデル化し、解いていきましょう。

PuLP ライブラリのインストール

まず Python の線形プログラミングライブラリである PuLP をインストールする必要があります。そのために以下をコマンド実行します。

$ pip install pulp

Python での実装

以下に PuLP を使用して数独を解くための Python コードを示します。

数独を表示するための関数定義

PuLP で解いていく数独の問題、および解いた後の結果を見やすく表示するための関数を作成します。

def show_sudoku(puzzle):
    # 数独の盤面を視覚的に表示する関数
    for r in range(len(puzzle)):
        # 3行ごとに区切り線を表示して盤面を区切る
        if r in [0, 3, 6]:
            print("+-------+-------+-------+")
        for c in range(len(puzzle[r])):
            # 各3列ごとに左境界を表示
            if c in [0, 3, 6]:
                print("| ", end="")
            # セルに数字がある場合は数字を、ない場合は空白を表示
            if puzzle[r][c] != 0:
                print(puzzle[r][c], end=" ")
            else:
                print(end="  ")
            # 行の最後に右境界を表示
            if c == 8:
                print("|")
    # 最後に底の境界線を表示
    print("+-------+-------+-------+")

数独を解くための関数定義

数独の初期設定をした後、制約処理を定義し、prob.solve() で線形プログラミングで問題を解いていきます。
pulp.lpSum() は、和や内積を計算する際に使用するもので、numpy で np.sum() するよりも高速で処理することができます。

import pulp

def solve_sudoku(puzzle):
    # 数独のグリッドを定義（1から9までの数字）
    VALS = ROWS = COLS = range(1, 10)

    # 数独問題を解くための線形プログラミングの問題インスタンスを作成
    prob = pulp.LpProblem("Sudoku Problem")

    # 各セルに入る可能性のある数字に対してバイナリ変数を定義（0または1）
    choices = pulp.LpVariable.dicts("Choice", [(val, row, col) for val in VALS for row in ROWS for col in COLS], cat="Binary")
    
    # 数独パズルの初期状態を設定（与えられた数値を固定）
    for row in ROWS:
        for col in COLS:
            if puzzle[row-1][col-1] != 0:  # 0は空欄を意味する
                for val in VALS:
                    if val == puzzle[row-1][col-1]:
                        prob += choices[val, row, col] == 1  # 与えられた数値を固定
                    else:
                        prob += choices[val, row, col] == 0  # 他の数値を排除

    # 数独のルールに従った制約を設定
    for val in VALS:
        # 各行において、各数字が1回だけ現れるように制約を設定
        for row in ROWS:
            prob += pulp.lpSum([choices[val, row, col] for col in COLS]) == 1
        # 各列において、各数字が1回だけ現れるように制約を設定
        for col in COLS:
            prob += pulp.lpSum([choices[val, row, col] for row in ROWS]) == 1
        # 各3x3ボックスにおいて、各数字が1回だけ現れるように制約を設定
        for box_row in range(3):
            for box_col in range(3):
                prob += pulp.lpSum([choices[val, 3*box_row + i, 3*box_col + j] for i in range(1, 4) for j in range(1, 4)]) == 1

    # 各セルには1つの数字のみが入るという制約を設定
    for row in ROWS:
        for col in COLS:
            prob += pulp.lpSum([choices[val, row, col] for val in VALS]) == 1

    # 線形プログラミング問題を解く
    prob.solve()

    # 解決策を取得して整理（数独盤面を作成）
    solution = [[0 for _ in range(9)] for _ in range(9)]
    for val in VALS:
        for row in ROWS:
            for col in COLS:
                if pulp.value(choices[val, row, col]) == 1:
                    solution[row-1][col-1] = val  # 解決策を盤面に反映
    return solution

問題の定義

実際に解いていく数独パズルを定義します。

# 数独パズルの例（0は空欄を意味する）
puzzle = [
    # ここに数独パズルの初期状態が格納される
    [7, 2, 0, 0, 0, 0, 0, 0, 0],
    [4, 0, 9, 5, 3, 7, 0, 8, 2],
    [5, 0, 0, 0, 1, 0, 0, 0, 0],
    [0, 0, 1, 0, 2, 0, 5, 6, 3],
    [3, 5, 0, 6, 8, 0, 1, 0, 9],
    [9, 7, 6, 0, 0, 0, 4, 2, 0],
    [6, 0, 7, 2, 4, 1, 8, 3, 0],
    [0, 0, 5, 9, 7, 8, 0, 4, 0],
    [2, 8, 0, 3, 6, 5, 9, 1, 0]
]

数独問題を解く

solve_sudoku() 関数を呼び出し、定義した数独問題を解きます。

# 数独パズルを解く
solution = solve_sudoku(puzzle)

# 問題と解答を表示

show_sudoku() 関数を呼び出し、数独の問題と問題を解いた後の解答を見やすく表示します。

print("問題:")
show_sudoku(puzzle)
print("")
print("解答:")
show_sudoku(solution)

実行結果

実装した Python コードを実行すると、以下のように解答が出力されます。

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /home/stack/.local/lib/python3.10/site-packages/pulp/solverdir/cbc/linux/64/cbc /tmp/d5ab89859a8e4b3ca16b8c6f6f72568a-pulp.mps timeMode elapsed branch printingOptions all solution /tmp/d5ab89859a8e4b3ca16b8c6f6f72568a-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 743 COLUMNS
At line 5533 RHS
At line 6272 BOUNDS
At line 7003 ENDATA
Problem MODEL has 738 rows, 730 columns and 3330 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 0 - 0.00 seconds
Cgl0002I 368 variables fixed
Cgl0004I processed model has 0 rows, 0 columns (0 integer (0 of which binary)) and 0 elements
Cbc3007W No integer variables - nothing to do
Cuts at root node changed objective from 0 to -1.79769e+308
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                0.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

問題:
+-------+-------+-------+
| 7 2   |       |       |
| 4   9 | 5 3 7 |   8 2 |
| 5     |   1   |       |
+-------+-------+-------+
|     1 |   2   | 5 6 3 |
| 3 5   | 6 8   | 1   9 |
| 9 7 6 |       | 4 2   |
+-------+-------+-------+
| 6   7 | 2 4 1 | 8 3   |
|     5 | 9 7 8 |   4   |
| 2 8   | 3 6 5 | 9 1   |
+-------+-------+-------+

解答:
+-------+-------+-------+
| 7 2 8 | 4 9 6 | 3 5 1 |
| 4 1 9 | 5 3 7 | 6 8 2 |
| 5 6 3 | 8 1 2 | 7 9 4 |
+-------+-------+-------+
| 8 4 1 | 7 2 9 | 5 6 3 |
| 3 5 2 | 6 8 4 | 1 7 9 |
| 9 7 6 | 1 5 3 | 4 2 8 |
+-------+-------+-------+
| 6 9 7 | 2 4 1 | 8 3 5 |
| 1 3 5 | 9 7 8 | 2 4 6 |
| 2 8 4 | 3 6 5 | 9 1 7 |
+-------+-------+-------+

最後に

今回は、PuLP を使って、無事に数独の正しい答えを導き出すことができました。
PuLP を使うと様々な最適化問題を解くことができるので、今後いろいろ挑戦していこうと思います。