Advent Calendar 2日目の記事 組合せ最適化でカックロを解く
数独を通して組合せ最適化を学ぼう
i. 目的
「様々な課題を組合せ最適化を使って解決できたら」と思った人に向けて、数独1を題材にして説明していきます。
一般に、パズルは解く過程を楽しむものです。
組合せ最適化を使って解くパズルは、モデル化の過程を楽しむものです。
ii. 組合せ最適化とは
数理モデルを使った最適化で、下記のような、いろいろな問題があります。
- 参考：典型問題と実行方法
もう少し詳しい説明は、下記をご覧ください。
- 参考：組合せ最適化を使おう
組合せ最適化は、難しい問題が多いですが、ソルバーを使うと数理モデルを作成するだけで、ソルバーが答えを出してくれます。
最近では、ソルバーの性能が上がってきて、過去、解けなかったものも解けるようになってきています。
iii. 手順について
課題を解くための手順です。
- 問題をまとめる
- ヒアリングしたり、頭を整理したりして、必要な情報をまとめていきます。
- 数理モデルを考える
- 変数、目的関数、制約条件を考えます。たくさんのアプローチ方法を考えることが重要です。この能力は、経験によって鍛えられます。
- 実装して解く
- ソルバーで実行できるモデルを記述します。お勧めするのは Pythonで記述する方法です。解くのはソルバーを呼ぶだけです。エラーを取り除くデバッグも含まれます。時間がかかるなどして解けない場合は、2.に戻って別のアプローチを探しましょう。
- 結果を確認する
- 結果を見てモデルの通りに実装されているか確認しましょう。テスト結果がおかしければ、3.に戻って直してください。モデル通りにできていたら、次にその課題の専門家に結果を見てもらいましょう。ここでは、図示が有効です。専門家に おかしいと言われたら、1.に戻って見直します。
1. 問題をまとめる
今回は、下図のような数独を解きます。
- 9x9の全マスに、1～9の数字を必ず埋めます。(1)
- どの1行(2)、どの1列(3)、どの3x3(4)において、同じ数字は1回だけ使えます。
- 数字が埋まっている場所は、その数字を使います。(5)
2. 数理モデルを考える
下図のような、9x9x9の729個の箱を考えましょう。この3軸は、行、列、数字に対応します。
1つの箱には、選ばれているか／選ばれていないか のどちらかの状態を持ちます。
これを1と0の数字で表しましょう。行i 列j 数字k の箱が1の場合、i行j列のマスの数字がkであることを意味するものとします。この729個の箱が0-1変数になります。
制約条件は、下記の通りです。
- 1つのマスに同じ数字は1つ。1x1x9の9箱の合計が1。 …(1)
- 1行のマスに同じ数字は1つ。1x9x1の9箱の合計が1。 …(2)
- 1列のマスに同じ数字は1つ。9x1x1の9箱の合計が1。 …(3)
- 3x3のマスに同じ数字は1つ。3x3x1の9箱の合計が1。 …(4)
- 数字指定の場合、その数字。1x1x1の1箱の合計が1。 …(5)
今回、目的関数はありません（便宜上、式 0 になります）。
3. 実装して解く
Pythonで数理モデルを作成する方法は、いろいろありますが、よく使われる PuLPとpandasを利用します。
PuLPについては、下記を参照ください。
インストール
Pythonはバージョン3.6をインストールしてください。
追加ライブラリのために、下記を実行してください(処理系によっては、pipの代わりにpip3になります)。
pip install pandas pulp ortoolpy
変数表の作成
pandasで下記のような表を作成します。1レコード1変数(Var)に対応させ729レコードになります（1レコードは表の1行に対応します）。
| 行 | 列 | 数 | _3x3 | 固 | Var | |
|---|---|---|---|---|---|---|
| 0 | 0 | 0 | 1 | 0 | False | v000001 |
| 0 | 0 | 0 | 2 | 0 | False | v000002 |
| ... |
この表を使うと制約条件は、以下のように表せます（カッコ内の数字は、2. 数理モデルを考えると対応します）。
- 1つのマスに同じ数字は1つ。「行と列」が同じレコード集合の変数の和は1。 …(1)
- 1行のマスに同じ数字は1つ。「行と数」が同じレコード集合の変数の和は1。 …(2)
- 1列のマスに同じ数字は1つ。「列と数」が同じレコード集合の変数の和は1。 …(3)
- 3x3のマスに同じ数字は1つ。「3x3の種別と数」が同じレコード集合の変数の和は1。 …(4)
- 数字指定の場合、その数字。「固」がTrueなら変数は1。 …(5)
キーが同じ集合は pandasの groupで、できます。
プログラム
import re, pandas as pd
from itertools import product
from pulp import LpProblem, lpSum, value
from ortoolpy import addbinvars
s = ('. . 6 |. . . |. . 1 '
'. 7 . |. 6 . |. 5 . '
'8 . . |1 . 3 |2 . . '
'------+------+------'
'. . 5 |. 4 . |8 . . '
'. 4 . |7 . 2 |. 9 . '
'. . 8 |. 1 . |7 . . '
'------+------+------'
'. . 1 |2 . 5 |. . 3 '
'. 6 . |. 7 . |. 8 . '
'2 . . |. . . |4 . . ')
data = re.sub(r'[^\d.]','',s) # 数字とドット以外を削除
r = range(9)
a = pd.DataFrame([(i,j,(i//3)*3+j//3,k+1,c==str(k+1))
for (i,j),c in zip(product(r,r),data) for k in r],
columns=['行','列','_3x3','数','固'])
a['Var'] = addbinvars(len(a))
m = LpProblem()
for cl in [('行','列'),('行','数'),('列','数'),('_3x3','数')]:
for _,v in a.groupby(cl):
m += lpSum(v.Var) == 1
for _,r in a[a.固==True].iterrows():
m += r.Var == 1
m.solve() # ソルバーで求解
4. 結果を確認する
結果を列として表に追加しましょう。
結果が1（＝選ばれた）のレコードを抜き出して、9x9に整形して表示します。
a['Val'] = a.Var.apply(value)
print(a[a.Val>0.5].数.values.reshape(9,9))
>>>
[[5 3 6 8 2 7 9 4 1]
[1 7 2 9 6 4 3 5 8]
[8 9 4 1 5 3 2 6 7]
[7 1 5 3 4 9 8 2 6]
[6 4 3 7 8 2 1 9 5]
[9 2 8 5 1 6 7 3 4]
[4 8 1 2 9 5 6 7 3]
[3 6 9 4 7 1 5 8 2]
[2 5 7 6 3 8 4 1 9]]
A. Q&A
Q: 4. 結果を確認するで、a.Val>0.5は、a.Val==1ではないのですか？
- A: ソルバーでは、整数変数も実数として計算しています。そのため、ごくたまに、0.99999999などの出力になることもあります。従って、等号ではなく不等号で判断しています。
Q: lpSumとは何ですか？sumではダメですか？
-
A: 機能的には同じですが、
lpSumの方が効率的なので、lpSumを使いましょう。