PuLPクイック入門

PuLPは, (混合) 整数計画問題を解くことができるPythonライブラリです.

問題の定義: m = pulp.LpProblem(name, sense)
- name: モデルの名前. 任意の名前でOK.
- sense: pulp.LpMinimize(最小化問題)かpulp.LpMaximize(最大化問題)
変数の定義: pulp.LpVariable(name, lowBound, upBound, cat)
- name: 変数の名前. 任意の名前でOKだが, 他の変数と重複してはいけない.
- lowBound: 下限値. 省略すると下限なし.
- upBound: 上限値. 省略すると上限なし.
- cat: 変数の種類. 文字列で指定する.
  - "Binary": 0, 1の二値変数.
  - "Integer": 整数変数.
  - "Continuous": 連続変数.
目的関数の設定: m += 式
- 例: m += x+y
  - +=という表記がややこしいが, 目的関数を何度も設定すると上書きされる
  - m.setObjective(式)だと思うとわかりやすい
制約の書き方: m += 条件式
- 制約が追加される
- 条件式で使える種類 (!=, >, <は使えない)
  - 式1 >= 式2
  - 式1 == 式2
  - 式1 <= 式2
- 例: m += x>=y
  - x>=yの制約が追加されます
問題を解く: m.solve()
解いたあとのステータスを得る: pulp.LpStatus[m.status]
- "Undefined"
- "Unbounded"
- "Infeasible"
- "Not Solved"
- "Optimal"
変数や目的関数の値を得る: pulp.value(・)
- 変数の値: pulp.value(変数)
- 目的関数の値: pulp.value(m.objective)
便利な関数
- pulp.lpSum(・): 合計をとる
- pulp.lpDot(係数, 変数): 内積をとる
- pulp.LpVariable.dicts(name, (iter1, iter2, ...), lowBound, upBound, cat): 変数を一気に作る
  - 例: vars = pulp.LpVariable.dicts("X", (range(3), range(4)))
    - f"X_{i}_{j}"の名前のついた変数が生成される(i=0,1,2, j=0,1,2,3)
    - vars[i][j]のように各変数にアクセスできる
      - numpyのようにvars[:, j]こういうアクセスができないので, このような表現を多用したいなら, vars = np.array([[pulp.LpVariable(f"X_{i}_{j}") for j in range(4)] for i in range(3)])で生成したほうがいい

数独

数独はニコリの登録商標です (ナンプレとも呼ばれてます).
「各行, 各列, 3x3の枠内の数字は重複してはいけない」という条件のもと, 以下のような9x9のマスに1-9の数字を埋めていく問題です.

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

今回はこの3x3のブロックを$m \times n$として一般化します.
盤面のサイズは$mn \times mn$とします. ($m=n=3$が普通の数独です)
たとえば, $m=2, n=3$ならこうなります.

---------------------------
|  3       1 |  6         |
|      5     |  1   2     |
---------------------------
|          6 |  3         |
|  4       2 |            |
---------------------------
|  2       3 |            |
|            |            |
---------------------------

ソルバー

定式化

なにも考えずにやろうとすると, 盤面の変数行列を$X\in\{1,...,mn\}^{mn \times mn}$として, たとえば行の制約を$X_{ij_1}\ne X_{ij_2} (j_1\ne j_2)$としたくなります.
しかし, $\ne$は使えないので別の表現をする必要があります.

変数: $X\in\{0,1\}^{mn \times mn \times mn}$
- $X_{ijk}$: 盤面$(i,j)$の値が $k+1$ であるなら $1$
  - 添字は0始まりとしてます
定数: $P\in\{0,1,...,mn\}^{mn \times mn}$
- 与えられる問題を表す行列
- 0は空欄を表す
制約:
- $\displaystyle\sum_k X_{ijk}=1, \forall i,j$
  - 1マスには1つの数字しか入らない
- $\displaystyle\sum_j X_{ijk}=1, \forall i,k$
  - 行の数字は重複してはいけない
- $\displaystyle\sum_i X_{ijk}=1, \forall j,k$
  - 列の数字は重複してはいけない
- $\displaystyle\sum_{i\le I<i+m,~j \le J<j+n}X_{IJk}=1, \forall i\in\{0,m,2m,...,mn\},j\in\{0,n,2n,...,mn\},k$
  - $m \times n$ブロック内は重複してはいけない
- $X_{ij,P_{ij}-1}=1, \forall i,j, P_{ij}\ne 0$
  - 与えられた数字を満たさないといけない
目的関数: なし
- あえて書くなら, 0 (定数)

コード

class Sudoku:
    """数独クラス"""
    def __init__(self, m, n):
        """考える小さいブロックのサイズ(m, n): 普通の数独ならm=3, n=3"""
        self.m = m
        self.n = n
        # 入る数字の最大数
        self.max_num = m * n
        # ループ用
        self.r = range(self.max_num)

    def init(self, problem):
        """与えられた問題を解くための初期化用関数
         problem: 0が空欄をあらわすnp.ndarray(int)
        """
        # 最適化問題
        self.sudoku = pulp.LpProblem("sudoku", pulp.LpMinimize)
        # max_num x max_numの盤面の数を表す変数
        self.board = np.array(
            [
                [
                    [
                        pulp.LpVariable(f"board_{i}_{j}_{k}", cat="Binary")
                        for k in self.r
                    ]
                    for j in self.r
                ]
                for i in self.r
            ]
        )
        # 数字はどれか一つ
        for i in self.r:
            for j in self.r:
                self.sudoku += pulp.lpSum(self.board[i,j,:]) == 1
        # 行に同じ数字が入ってはいけない
        for i in self.r:
            for k in self.r:
                self.sudoku += pulp.lpSum(self.board[i,:,k]) == 1
        # 列に同じ数字が入ってはいけない
        for j in self.r:
            for k in self.r:
                self.sudoku += pulp.lpSum(self.board[:,j,k]) == 1
        # 正方形の枠に同じ数字が入ってはいけない
        for i in range(0, self.max_num, self.m):
            for j in range(0, self.max_num, self.n):
                for k in self.r:
                    self.sudoku += pulp.lpSum(self.board[i:i+self.m, j:j+self.n, k]) == 1
        # 初期値をいれる
        for i in self.r:
            for j in self.r:
                if problem[i, j]:# 空欄でないとき
                    self.sudoku += self.board[i, j, problem[i, j]-1] == 1

定式化した表現を素直に書いただけです.
コード内では$X$をself.boardで表しています.

すべての解を列挙する

先程の数独クラスで, self.init(problem)してself.sudoku.solve()すれば解が得られますが, 1つしか得られません.
すべての解を列挙するためには, 一度出た解を再び出力しないような制約を加えて, 再び最適化する必要があります.

得られた1つの解を$\tilde{X}$とします.
各マスには1つの数字しか入ってないので, 当然$\tilde{X}$の合計値は, $\displaystyle\sum_{ijk}\tilde{X}_{ijk}=(mn)^2$です.
別の解を得るためには, 「割り当てられた箇所の合計 $\le (mn)^2-1$」 という制約を追加します.
式で書くと,

\sum_{ijk: \tilde{X}_{ijl}=1 } X_{ijl} \le (mn)^2-1

この制約を加えることで, $\tilde{X}$の割り当てから, 最低でも1つは割り当てを変えなければいけなくなります.

    def solve(self, problem):
        """ソルバー. 全ての解が順番に得られる
         problem: 0が空欄をあらわすnp.ndarray(int)
        """
        self.init(problem)
        while True:
            self.sudoku.solve()
            if pulp.LpStatus[self.sudoku.status] != "Optimal":
                break
            answer = np.vectorize(pulp.value)(self.board)
            answer = np.where(answer)[2].reshape(self.max_num,self.max_num)+1
            yield answer
            # 割り当てられた変数から最低でも1つ外すようにすることで別の解を探すようにする
            self.sudoku += pulp.lpSum([
                self.board[i, j, k]
                for i in self.r
                for j in self.r
                for k in self.r
                if pulp.value(self.board[i, j, k]) == 1
            ]) <= self.max_num * self.max_num-1

これですべての解を列挙できるようになりました.

問題作成

作ったソルバーを使って数独の問題を作ります.
ここで, 問題とは, 解が存在し, 一意であるものと定義します.

作った問題の解が一意であるかの判定は, すべての解を列挙するself.solve(・)を使って表現します.
解の個数が2以上になった時点で一意ではありません.
解の個数が0なら解は存在しません.

    def is_unique(self, problem):
        """解が一意かどうか判定
         True: 一意
         False: 一意でない
         None: 解なし
        """
        cnt = 0
        for ans in self.solve(problem):
            cnt+=1
            # 一意でない
            if cnt>=2:
                return False
        # 一意
        if cnt==1:
            return True
        # 解なし
        return None

単純な方法

空の盤面にランダムに数字を入れていき, ソルバーで解が一意かどうか判定しながら, 一意になるまで数字を追加します.

    def create_problem_from_empty(self, try_num):
        """問題作成 (空の盤面にランダムに数字を入れていくバージョン)
         実行速度遅い.
        """
        # 空の盤面作成
        problem = np.zeros((self.max_num, self.max_num), dtype=int)
        # 埋めていく
        for _ in range(try_num):
            print(">", end="")
            while True:
                i, j = np.random.randint(self.max_num, size=2)
                if not problem[i, j]:
                    break
            # 盤面(i,j)を適当に埋める
            problem[i, j] = np.random.randint(self.max_num) + 1
            is_unique = self.is_unique(problem)
            # 一意でないなら盤面を空白に戻す
            if is_unique is None:
                problem[i, j] = 0
                continue
            # 一意なら終了
            if is_unique:
                print("")
                return problem
        return None

try_numは埋める試行をする回数です. 十分大きな値を設定する必要があります.
単純で実装しやすいですが, はじめのうちはほとんど空欄の状態の数独を何度も解くことになるので, とても時間がかかります.
$𝑚𝑛>10$からきつくなってきました.

埋まった状態の盤面から穴を開けていく方法

この方法では, 最初に答えの盤面を用意してから空欄にしていきます.
ソルバーにとっては, ほとんど埋まった状態から解くことになり, 探索空間が少なく楽になります.
一方で, すべて埋まった盤面を用意するのが必要になります.
真面目にやると計算時間がやばいので, ヒューリスティックに埋めていきます.


    def __fill_rows(self, problem, i, try_num):
        """i行以降をランダムに埋める再帰関数 (埋まった状態の盤面作成用ヒューリスティック)"""
        if i >= self.max_num:
            return True
        for _ in range(try_num):
            if not self.__fill_vals(problem, i, 0, try_num):
                continue
            if self.__fill_rows(problem, i+1, try_num):
                return True
        return False

    def __fill_vals(self, problem, i, j, try_num):
        """i行目のj列以降をランダムに埋める再帰関数 (埋まった状態の盤面作成用ヒューリスティック)"""
        if j >= self.max_num: # 最後の列まで埋まったら
            return True
        for _ in range(try_num):
            problem[i, j] = np.random.randint(self.max_num)+1
            if self.check_duplicate(problem, i, j):
                continue # 重複したら別の数字をトライ
            if self.__fill_vals(problem, i, j+1, try_num): #次の値を埋めに行く
                return True
        problem[i, j] = 0
        return False

    def check_duplicate(self, problem, i, j):
        """(i,j)が属する行、列、ブロックで、0以外の数字が重複してたらTrue"""
        row = problem[i, :]
        zero_row = (row==0).any()
        col = problem[:, j]
        zero_col = (col==0).any()
        blk = problem[i//self.m*self.m:(i//self.m+1)*self.m, j//self.n*self.n:(j//self.n+1)*self.n]
        zero_blk = (blk == 0).any()
        return (np.unique(row).shape[0]-zero_row != (row!=0).sum() or \
            np.unique(col).shape[0]-zero_col != (col!=0).sum() or \
            np.unique(blk).shape[0]-zero_blk != (blk != 0).sum()
        )

0行目から順番に埋めていきます.
「 $i$ 行以降を埋める」という動作は, 「 $i$ 行目を埋める」+「 $i+1$ 行目以降を埋める」という表現ができます. (再帰)
「 $i$ 行目の $j$ 列以降を埋める」という動作は, 「 $i$ 行目の $j$ 列目を埋める」+「 $i$ 行目の $j+1$ 列以降を埋める」という表現ができます. (再帰)

try_numは, その行を作る試行を行う回数, および, その値を埋める試行を行う回数です.
再帰で書くことで, 「この行何回試行してもうまくいかないな」となったら前の行に戻り前の行を疑って修正します.
ヒューリスティックなので, 失敗することもあります.

次に, この埋められた盤面を使って問題を作ります.

    def create_problem_from_fill(self, try_num, level=1):
        """問題作成 (埋まった状態から穴をあけていくバージョン)
         最初に答えの盤面を作る必要がある. 
         try_num: 答えの盤面作成で試行する回数. m*nくらいが良さそう. 
         level: 穴あけ失敗許容回数
        """
        # 空の盤面作成
        problem = np.zeros((self.max_num, self.max_num), dtype=int)
        # 盤面を全て埋める
        if not self.__fill_rows(problem, 0, try_num):
            return None

        cnt = 0 # 穴あけに失敗した回数
        while True:
            print(">", end="")
            while True:
                i, j = np.random.randint(self.max_num, size=2)
                if problem[i, j]:
                    break
            tmp = problem[i, j]
            # 盤面に穴をあける
            problem[i, j] = 0
            if not self.is_unique(problem):
                # 一意でなくなったら値をもとに戻す
                problem[i, j] = tmp
                cnt += 1
                if cnt >= level:
                    break
        print("")
        return problem

盤面に穴をあけて, 解いた結果, 一意ならOK.
一意でないなら値をもとに戻して終了します. (穴あけ失敗)
この穴あけ失敗回数の許容回数をlevelとしてます.
levelが大きくなると, なるべく空欄の多い問題を作ろうとするので, 問題の難易度が高くなります.

実行結果

ちゃんと計算時間の平均とるの面倒なので, だいたいで書きます.
Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz 1.80 GHzにて.

問題作成時間:

m	n	単純	穴あけ(level=1)	穴あけ(level=2)
2	3	2s-3s	1s-2s	2s-3s
3	3	15s-30s	3s-7s	6s-8s
3	4	3min	13s-15s	45s-52s
4	4	未実行	50s-60s	2min40s

prologで実装したときはもう少し早かったのでかなり大きいサイズまでやれました.
prologソースのコメントに詳しく書いてます. (そもそもprolog知っている人どれくらいいるだろうか. )
こういう整数計画はprologの方が得意な気がする.
バックトラックで自動的にすべての解を出してくれるし.

最後に自動生成した問題をいくつかおいておきます. やりたい方はどうぞ.

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

----------------------------------------------------
|          7   4 |     10   9     |      3   1     |
|                |  7       1   2 |             10 |
|          1     |  3           6 |                |
----------------------------------------------------
|  5             |          2     |      9  10  11 |
|     10       1 |      7         |  6   5      12 |
|     12   8     |      6       9 |      4         |
----------------------------------------------------
|                |  9             |  2       3     |
|          3   8 |     11  12  10 |  5             |
|                |      3       5 |                |
----------------------------------------------------
|      8   9     |  5       6     |              7 |
|          5  12 |                |  3           4 |
|      3         |  8      10     |                |
----------------------------------------------------

おまけ

JavaScriptで書いてブラウザで動かせるようにしておきました。