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

データ分析（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

以上