最近になってプログラマーもすなる最適化というものに興味が出てきた。
今回はPythonのPuLPで簡単なカックロを解いてみたい。
記事中の追記の通り、この記事には誤りが含まれる
準備
バージョンは次の通りとする。
| バージョン | |
|---|---|
| Python | 3.9.2 |
| PuLP | 2.5.0 |
まずはPuLPのインストール
$ python --version 3.9.2 $ pip install pulp (中略) Successfully installed pulp-2.5.0
ドキュメンテーションは次のサイトを参照。
Optimization with PuLP — PuLP 2.7.0 documentation
PuLPの感覚の把握
まずは簡単な問題でデータ型の感覚をつかんでみる。
PuLPには大きく、変数のクラス（pulp.LpVariable）、式のクラス（pulp.LpAffineExpression）、制約のクラス（pulp.LpConstraint）、問題のクラス（pulp.LpProblem）、の4つのクラスが存在する。ドキュメントは下記サイトを参照。
pulp: Pulp classes — PuLP 2.7.0 documentation
PuLPでは、最初にLpProblemクラスのオブジェクトを作る。
import pulp prob = pulp.LpProblem('prob', sense=pulp.LpMaximize) print(type(prob)) # <class 'pulp.pulp.LpProblem'> print(prob) # prob: # MAXIMIZE # None # VARIABLES #
pulp.LpProblem()コンストラクタの引数は2つしかない。これほど引数が少ないのも珍しい話である。
第一引数nameには問題の名前（文字列）を記述する。ドキュメントによると.lpファイル（出力ファイル）の名前に使う。将来性を考えてまともな名前を付けておいた方が良いだろう。第二引数senseにはpulp.LpMinimize（デフォルト）かpulp.LpMaximizeを指定する。
できあがったLP問題クラスは、問題設定がされていない状態になっている。 ここに変数や目的関数、制約条件を追加していく。まずは変数を定義する。
v1 = pulp.LpVariable('v1', cat='Integer', lowBound=0, upBound=9) v2 = pulp.LpVariable('v2', cat='Integer', lowBound=0, upBound=9) print(type(v1)) # <class 'pulp.pulp.LpVariable'> print(v1) # v1
変数を定義するさいに、先ほどの問題クラスと紐づける必要はない。これは感覚としては少し奇妙だが、そういうものなのだろう。
pulp.LpVariable()コンストラクタの引数は5つである。nameは名前（文字列）を指定する。catは変数のカテゴリを指定する。カテゴリには'Continous'（デフォルト）、'Integer'、'Binary'の3種類を指定できる。ドキュメントには引数の型が書かれていないが、文字列で指定できる。しかしどうやら、文字列で指定する以外に、pulp.LpIntegerのように指定する方法もあるようだ。おそらく結果は同じである。lowBoundとupBoundは言うまでもなく、変数の下限上限を指定する。カテゴリがBinaryなら指定する意味はないと考えれば、指定するのは実数か整数であろう。変数の可動範囲は、境界を含むようだ。コンストラクタは5つ目の引数としてeという引数を持つ。これはcolumn based modelingという手法で用いるようだが、今は考えない。
print(<変数>)で変数名が出てくる、というのは面白い点で、変数や式を確認したいときに役立つだろう。
次に問題を定義する。
prob += (v1 + v2) # target for maximization prob += (v1 - v2 == 1) # constraint 1 prob += (v1 <= 5) # constraint 2
ここで何が起こっているのかは、一見してわかりにくい。1行目と2行目&3行目で用いている構文はほとんど同一だが、1行目では目的関数を登録、2行目では制約条件を登録している。
より具体的には、v1 + v2を実行することでpulp.LpAffineExpressionクラスのオブジェクトができる。+=の演算はLpProblemオブジェクトの__iadd__によって定義されているのだろう。対してv1 - v2 == 1やv1 <= 5の実行では、pulp.LpConstraintクラスのオブジェクトができる。同じ構文なのは一見して奇妙だが、おそらくLpProblemクラスの__iadd__は引数の型によって処理を変更しているのだろう（C++で言うところのオーバーロード）。可読性は落ちるので、コメントを入れておいた方が良さそうだ。
目的関数や制約を定義する際に、+や==や<=を当たり前のように使っているが、いかにも危なっかしい。ここで使える演算の種類も調べておいた方が良さそうだ。
ここまでで、問題が定義できていることを確認してみる。
print(prob) # prob: # MAXIMIZE # 1*v1 + 1*v2 + 0 # SUBJECT TO # _C1: v1 - v2 = 1 # # _C2: v1 <= 5 # # VARIABLES # 0 <= v1 <= 9 Integer # 0 <= v2 <= 9 Integer
出力を見ると、MAXIMIZEに登録した目的関数が、SUBJECT TOに登録した2つの制約が、VARIABLESに登録した2つの変数の型と範囲が、期待通りに表示されていることが確認できる。
目的関数と制約は、それぞれ、問題オブジェクトのobjective属性とconstraints属性に保存されている。さもありなん。
print(prob.objective) # v1 + v2 print(prob.constraints) # OrderedDict([('_C1', 1*v1 + -1*v2 + -1 = 0), ('_C2', 1*v1 + -5 <= 0)])
問題設定は完了した。いよいよ問題を解いてみる。
status = prob.solve() print(status) # 1 print(f'v1 = {v1.value()}') # v1 = 5.0 print(f'v2 = {v2.value()}') # v2 = 4.0 print(f'target value = {prob.objective.value()}') target value = 9.0
prob.solve()で問題を解く。引数solverでソルバーを指定する。デフォルトではCBCが使える。他のソルバーとしてGLPKなどがあるが、別途インストールが必要らしい。他に、おそらくソルバーに依存する引数が存在するが、今回はいったん置いておく。
ソルブの間は画面にメッセージが出力される。ソルブが完了すると、上のコード片のようにv1.value()で結果にアクセスできる。
statusが正常だったかを確認するには、pulp.LpStatusに照らし合わせる。
print(pulp.LpStatus) # {0: 'Not Solved', 1: 'Optimal', -1: 'Infeasible', -2: 'Unbounded', -3: 'Undefined'}
ここまでをまとめておく。
import pulp def test_pulp(): prob = pulp.LpProblem('prob', sense=pulp.LpMaximize) v1 = pulp.LpVariable('v1', cat='Integer', lowBound=0, upBound=9) v2 = pulp.LpVariable('v2', cat='Integer', lowBound=0, upBound=9) prob += (v1 + v2) # target for maximization prob += (v1 - v2 == 1) # constraint 1 prob += (v1 <= 5) # constraint 2 status = prob.solve() print('Optimization Status:', pulp.LpStatus[status]) print('Result:') print(f'\tv1 = {v1.value()}') print(f'\tv2 = {v2.value()}') print(f'\ttarget value = {prob.objective.value()}')
バイナリーの例
変数カテゴリとしてBinaryを使う例も載せておく。
def test_pulp_binary(): values = [1, 3, 4, 18, 20, 31, 35, 52, 58, 70, 78, 88, 127] prob = pulp.LpProblem('prob') bools = [pulp.LpVariable(f'value_{v}', cat='Binary') for v in values] prob += (pulp.lpDot(values, bools) == 200) status = prob.solve() if pulp.LpStatus[status] != 'Optimal': print('failed') return selected_values = [v for (v, b) in zip(values, bools) if (b.value())] print('+'.join([str(v) for v in selected_values]))
これは数の和が200となる組み合わせを見つける例で、部分和問題と呼ばれる。下のサイトに部分和問題の例と、応用例が載っている。
上の例では目的関数を設定していない。目的関数がなくても良いということである。また、上ではpulp.lpDot()という関数を用いている。これは1次元的な変数列がある場合の線形演算の表現に便利である。他にpulp.lpSum()というものもある。特にpulp.lpSum()は、下のサイトによるとリスト内包表記のような構文が使えるらしい。
pulpを使いこなすための備忘録 - 数学やプログラミングの備忘録
簡単なカックロへの適用
問題は簡単に、図1のものを使う。
2021/10/11追記：PuLPでは「!=」が使えないらしい。この例ではたまたまうまく動いているように見えたが、問題設定が正しくできていなかった。一般にはあるマスにある数が入るかどうかで一つのブール変数にする必要があるようだ。
def test_solve_kakuro_with_pulp(): prob = pulp.LpProblem('prob') v11, v12, v21, v22 = [ pulp.LpVariable(s, cat='Integer', lowBound=1, upBound=9) for s in 'v11,v12,v21,v22'.split(',')] # Constraints # ヨコのカギ prob += (v11 + v12 == 3) prob += (v21 + v22 == 8) # タテのカギ prob += (v11 + v21 == 4) prob += (v12 + v22 == 7) # 排他性 prob += (v11 != v12) prob += (v21 != v22) prob += (v11 != v21) prob += (v12 != v22) status = prob.solve() if pulp.LpStatus[status] != 'Optimal': print('failed') return print('Answer:') print(v11.value(), v12.value()) print(v21.value(), v22.value())
結果は次の通り
Answer: 1.0 2.0 3.0 5.0
なぜか結果がfloat型になっているが、答えは得ることができた。
今回の実装がどの程度実用に足りるものか不明だが、とりあえず一歩目は踏み出すことができた。
