Python
数学
最適化
パズル
組合せ最適化

魔方陣を通して組合せ最適化を学ぶ

データ分析(pandas)と最適化(pulp)を使ったモデリング

3×3の魔方陣を題材に、Pythonで組合せ最適化のモデルの作り方を学びましょう。

数理モデル

魔方陣を表す数理モデルは、言葉で表すと以下のようになります。

変数:行、列ごとの数字(1から9)
目的関数:なし
制約条件:各数字はちょうど1つ使う
     各行の和は15
     各列の和は15
     斜めの和は15

テクニック01)特定の数字に着目する場合、0-1変数を利用する
テクニック02)1変数が1行に対応する変数(Var)の表(df)を利用する

上記のテクニックを使って書き換えると、以下のようになります。

変数:Var(属性として、「行」、「列」、「数」を持つ0-1変数の表)
制約条件:1つのマス(行と列が同じ集まり)のVarの和=1
     数が同じ集まりのVarの和=1
     行が同じ集まりのVarと数の内積の和=15
     列が同じ集まりのVarと数の内積の和=15
     斜めの和=15

Pythonによる数理モデル

まずは、変数表(df)を作成します。

python3
import pandas as pd
from pulp import LpProblem, LpVariable, LpBinary, lpSum, lpDot, value

df = pd.DataFrame([(i,j,k,LpVariable(f'Var_{i}_{j}_{k}',cat=LpBinary))
                   for i in range(3)
                   for j in range(3)
                   for k in range(1,10)],
                  columns=['行','列','数','Var'])
print(df[:3])
print(df[-3:])
>>>
               Var
0  0  0  1  Var_0_0_1
1  0  0  2  Var_0_0_2
2  0  0  3  Var_0_0_3
                Var
78  2  2  7  Var_2_2_7
79  2  2  8  Var_2_2_8
80  2  2  9  Var_2_2_9

変数表は、3行×3列×9数字で、81レコードあります。

続いて、下記のテクニックを使って、Pythonの数理モデル(model)を作成してみましょう。

テクニック03)制約条件を指定するには、「数理モデル += 制約条件」とする
テクニック04)「XXXが同じ集まり」→ df.groupby('XXX') を利用する
テクニック05)「XXXの和」→ lpSum を利用する
テクニック06)「XXXの内積の和」→ lpDot を利用する

python3
model = LpProblem()
for key, group in df.groupby(('行','列')):
    model += lpSum(group.Var) == 1
for key, group in df.groupby('数'):
    model += lpSum(group.Var) == 1
for key, group in df.groupby('行'):
    model += lpDot(group., group.Var) == 15
for key, group in df.groupby('列'):
    model += lpDot(group., group.Var) == 15
model += lpDot(*df.query('行==列')[['数','Var']].T.values) == 15
model += lpDot(*df.query('行+列==2')[['数','Var']].T.values) == 15

斜めの和を表すコードは、一旦、こういうものだと思ってください。

求解と結果の表示

下記のテクニックを使って、ソルバー(cbc)で解を求めて、結果を表示してみましょう。
reshape(3,3)で3×3の形に変えています。

テクニック07)ソルバーで数理モデルを解くには、solveをよぶ
テクニック08)数理モデルのstatusが1ならば最適解が得られている
テクニック09)変数の表に結果の列(Val)を追加する
テクニック10)変数の値は、value を利用する(列ならapply(value))
テクニック10)結果が1のものだけ取り出す → df[df.Val>0.5]

python3
model.solve()
if model.status == 1:
    df['Val'] = df.Var.apply(value)
    print(df[df.Val>0.5]..values.reshape(3,3))
>>>
[[4 3 8]
 [9 5 1]
 [2 7 6]]

魔方陣が計算できました。

補足説明

「行,列」でgroupbyすると、各マスごとの表がgroupとして取り出せます。

python3
for key, group in df.groupby(('行','列')):
    print(key)
    print(group)
    break
>>>
(0, 0)
               Var
0  0  0  1  Var_0_0_1
1  0  0  2  Var_0_0_2
2  0  0  3  Var_0_0_3
3  0  0  4  Var_0_0_4
4  0  0  5  Var_0_0_5
5  0  0  6  Var_0_0_6
6  0  0  7  Var_0_0_7
7  0  0  8  Var_0_0_8
8  0  0  9  Var_0_0_9

「数」でgroupbyすると、各数ごとの3×3の表がgroupとして取り出せます。

python3
for key, group in df.groupby('数'):
    print(key)
    print(group)
    break
>>>
1
                Var
0   0  0  1  Var_0_0_1
9   0  1  1  Var_0_1_1
18  0  2  1  Var_0_2_1
27  1  0  1  Var_1_0_1
36  1  1  1  Var_1_1_1
45  1  2  1  Var_1_2_1
54  2  0  1  Var_2_0_1
63  2  1  1  Var_2_1_1
72  2  2  1  Var_2_2_1

df[df.Val>0.5]とすると、選ばれたレコードだけ取り出せます。

python3
print(df[df.Val>0.5])
>>>
                Var  Val
3   0  0  4  Var_0_0_4  1.0
11  0  1  3  Var_0_1_3  1.0
25  0  2  8  Var_0_2_8  1.0
35  1  0  9  Var_1_0_9  1.0
40  1  1  5  Var_1_1_5  1.0
45  1  2  1  Var_1_2_1  1.0
55  2  0  2  Var_2_0_2  1.0
69  2  1  7  Var_2_1_7  1.0
77  2  2  6  Var_2_2_6  1.0

以上