Pyomoで魔法陣を解く

Pythonで整数計画問題

整数計画法で問題を解く必要が出てきたので、ライブラリを探してみた。バックエンドの最適化器はGLPKを使うとして、Pythonから使うフロントエンドには色々ある。あるのだが、なんだかうまくいかない。

• Python-MIP
  • M1 macbookではさっくり動かない。
• PulP
  • よさそうだし、小さいサンプルはさっくり動くのだけど、ちょっと大きめのサンプルを動かすと、GLPKでエラーがでて動かない。CPLEX形式のファイルを生成してglpsolを起動するのだけど、そのファイルの形式がGLPKが理解するフォーマットと微妙にズレているようで、エラーが出る

Pyomo

Pyomoは、線形計画問題を含む様々な種類の最適化問題を記述するためのライブラリで、バックエンドとしてGLPKを含むさまざまなパッケージを使用できる。githubを覗いてみると開発も続いているようなので、試してみた。

ConcreteモデルとAbstractモデル

PyomoにはConcreteモデルとAbstractモデルという概念がある。
前者はデータも埋め込んだコードで、後者はコードとデータを分離する。
後者のほうが一般には汎用性が高いのだが、機械的にコードを生成することを
目標としているので、Concreteモデルのほうを使う。

Pyomoプログラムの構成

• Set
  インデックスの空間を定義する。
• Param
  データファイルから読み込む値を定義する。
• Var
  最適化の対象となる変数。Setを用いて定義すると、Setで定義されたインデックスを用いた配列変数になる。
• Constraint
  制約を記述する。Setを指定することで、複数の制約式を同時に定義することができる。
• Objective
  最適化対象となる式の定義。最大化、最小化が選択できる。

魔法陣

魔法陣とは、NxNのグリッドに、1 - N*N の数をひとつづつ配置したもので、縦、横、対角線の数の和が常に等しくなるもの、と定義される。いかにN=3の魔法陣の例を示す。縦横対角線の和がすべて15になっている。

6 7 2 
1 5 9 
8 3 4

一般に Nの魔法陣内の数字すべての和は 1から$N^2$までを足し合わせたものであり、$(1+N^2)N^2/2$で表される。縦横はN 個あるので、各縦横の和は $(1+N^2)N^2/2N = (1+N^2)N/2$ で表される。N=3では15、N=4では34になる。

Pyomoによる魔法陣

2次元空間上に整数値を配置するという方法で良さそうな気がするが、それではうまくいかない。すべてのグリッドに別の値が入っている、という条件を表そうとすると != を用いて制約を書きたくなるのだが、線形計画問題ではこのような書き方はできないためだ。== が書けるのに、!= が書けないというのは直感的に理解できなくてびっくりした。
ではどうするかというと、深さ方向を追加して3次元空間で考える。深さはそのグリッドに置かれる値を示す。各グリッドには1/0のどちらかが入る。
たとえば x,y に d という値が置かれているということは array[x,y,p] == 1 where p == d
array[x,y,p] == 0 where p=!= d
で表現するということだ。この表現法は広く知られている方法らしく、こちらに書かれているのを参考に書いてみた。

from pyomo.environ import *

model = ConcreteModel()
N = 4
S = sum(range(1, N*N+1))//N

# model.Aの初期化に用いるジェネレータ
def init_idx2(_):
    for i in range(N):
        for j in range(N):
            yield (i,j)

# model.Bの初期化に用いるジェネレータ
def init_idx3(_):
    for i in range(N):
        for j in range(N):
            for k in range(N*N):
                yield (i,j,k)

#盤面のxもしくはyの1次元を表す
model.N = Set(initialize=range(N))
#盤面のz軸すなわち、置かれた数値を表す。
model.D = Set(initialize=range(N*N))
#盤面のxy平面2次元を表す。
model.A = Set(within=model.N*model.N, initialize=init_idx2)
#盤面のxyz 3次元全体を表す。
model.B = Set(within=model.N * model.N * model.D, initialize=init_idx3)

#盤面への配置を表す変数。0/1の2値で、1がある「深さ」で数値を表す
model.x = Var(model.B, within=Binary)

# X軸方向の和がSに等しいという制約
# 深さ方向に応じて重みを付けて足し合わせる。dではなくd+1にしているのは、深さが0オリジンであるため。 
def x_sum(M, p):
    return sum(M.x[(i,p,d)] * (d+1) for i in range(N) for d in range(N*N)) == S
model.x_sum_const = Constraint(model.N, rule = x_sum)

def y_sum(M, p):
    return sum(M.x[(p,i,d)] * (d+1) for i in range(N) for d in range(N*N)) == S
model.y_sum_const = Constraint(model.N, rule = y_sum)

# 各グリッドの値がそれぞれ異なることを示す制約
# 各深さでの平面上の値を加算して1であることを確認
def unique(M, d):
    return sum(M.x[i,j,d] for i in range(N) for j in range(N)) == 1
model.uqnique_const = Constraint(model.D, rule = unique)
        
# 各グリッドに1つだけしか値がないことを示す制約
# 深さ方向に加算した値が1であることで確認している
def onevalue(M, i, j):
    return sum(M.x[i,j,d] for d in range(N*N)) == 1
model.onevalue_const = Constraint(model.A, rule = onevalue)

# 右下がり対角に関する制約
model.diagonal = Constraint(expr=sum(model.x[(i,i,d)] * (d+1) for i in range(N) for d in range(N*N)) == S)

# 右上がり対角に関する制約
model.bdiagonal = Constraint(expr=sum(model.x[(i,N-1-i,d)] * (d+1) for i in range(N) for d in range(N*N)) == S)


opt = SolverFactory('glpk')
opt.solve(model, tee=False)
model.display()

解の表示

解は3次元の0/1配列なので解を見てもわからない。なので、以下の関数で表示する。

def display_board(M, N):
    for j in range(N):
        for i in range(N):
            for k in range(N*N):
                if M.x[i,j,k].value > 0:
                    print(k+1, end=" ")
                    break
        print("")

> display_board(model, N)
2 14 7 11 
3 15 6 10 
13 1 12 8 
16 4 9 5 

Constraintの書き方に関する注意点

Constraint の書き方が2種類ある。

model.diagonal = Constraint(expr=...)

のように1つの制約を書く方法と、下のようにインデックス空間を指定して一括して
多数の制約を記述する方法だ。

def onevalue(M, i, j):
    return sum(M.x[i,j,d] for d in range(N*N)) == 1
model.onevalue_const = Constraint(model.A, rule = onevalue)

rule で指定した関数で制約条件を記述する。この関数の引数には、モデルとインデックス空間のインデックス値が渡される。上の例では model.A が2次元空間であるため、インデックス i, j が与えられる。ここでは関数として書いたが、lambda式で書くことも可能。

Constraintに限らずさまざまな値にpprintが定義されていてこれを用いると定義式が確認できる。

> model.diagonal.pprint()
diagonal : Size=1, Index=None, Active=True
    Key  : Lower : Body                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      : Upper : Active
    None :  34.0 : x[0,0,0] + 2*x[0,0,1] + 3*x[0,0,2] + 4*x[0,0,3] + 5*x[0,0,4] + 6*x[0,0,5] + 7*x[0,0,6] + 8*x[0,0,7] + 9*x[0,0,8] + 10*x[0,0,9] + 11*x[0,0,10] + 12*x[0,0,11] + 13*x[0,0,12] + 14*x[0,0,13] + 15*x[0,0,14] + 16*x[0,0,15] + x[1,1,0] + 2*x[1,1,1] + 3*x[1,1,2] + 4*x[1,1,3] + 5*x[1,1,4] + 6*x[1,1,5] + 7*x[1,1,6] + 8*x[1,1,7] + 9*x[1,1,8] + 10*x[1,1,9] + 11*x[1,1,10] + 12*x[1,1,11] + 13*x[1,1,12] + 14*x[1,1,13] + 15*x[1,1,14] + 16*x[1,1,15] + x[2,2,0] + 2*x[2,2,1] + 3*x[2,2,2] + 4*x[2,2,3] + 5*x[2,2,4] + 6*x[2,2,5] + 7*x[2,2,6] + 8*x[2,2,7] + 9*x[2,2,8] + 10*x[2,2,9] + 11*x[2,2,10] + 12*x[2,2,11] + 13*x[2,2,12] + 14*x[2,2,13] + 15*x[2,2,14] + 16*x[2,2,15] + x[3,3,0] + 2*x[3,3,1] + 3*x[3,3,2] + 4*x[3,3,3] + 5*x[3,3,4] + 6*x[3,3,5] + 7*x[3,3,6] + 8*x[3,3,7] + 9*x[3,3,8] + 10*x[3,3,9] + 11*x[3,3,10] + 12*x[3,3,11] + 13*x[3,3,12] + 14*x[3,3,13] + 15*x[3,3,14] + 16*x[3,3,15] :  34.0 :   True

所感

いまのところ、PuLPのような理不尽な動作は見られず安定して使えている。GPLKのファイルを直接書くよりはずっと楽そうだ。
当面これでいこう。