Edited at

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

More than 1 year has passed since last update.


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

以上