Pythonで使える最適化モジュールPuLPを紹介します. なお,ここではPython3を用いています. PuLPのインストールのやり方などは,適当にWebで検索してみてください.
まず,以下の線形最適化問題 $$ \text{max. } 3 x + 4 y + 2 z $$ $$ \text{s. t. } 2 x \le 4, $$ $$ x + 2 z \le 8, $$ $$ 3 y + z \le 6, $$ $$ x \ge 0, \ y \ge 0, \ z \ge 0 $$ を解いてみます.
最初にPuLPをimportします.
import pulp
PuLPを使う準備ができたので,PuLPのメソッドLpProblemを用いて,線形最適化問題のinstanceを作ります. ここではinstanceの名前をLPとします.
LP = pulp.LpProblem(name='LP_example', sense=pulp.LpMaximize)
最大化問題か最小化問題かはsenseというパラメータで指定します. 最大化問題にしたい場合はsense=pulp.LpMaximize,最小化問題にしたい場合はsense=pulp.LpMinimizeとします. senseの指定を省略するとデフォルト値のpulp.LpMinimizeが適用されるみたいです.
次に,PuLPのメソッドLpVariableを用いて,線形最適化問題の変数のinstaceを作ります.
x = pulp.LpVariable(name='x', lowBound=0)
y = pulp.LpVariable(name='y', lowBound=0)
z = pulp.LpVariable(name='z', lowBound=0)
パラメータlowBoundは変数の下限です. 線形最適化問題を記述する際は通常,変数の下限・上限も制約式として記述します. しかし多くの最適化ソフトウェアでは,変数を設定する際に下限・上限も設定できます. PuLPも例外ではないようです.
次に,問題に目的関数を加えます.
LP += 3 * x + 4 * y + 2 * z
問題に線形式(の表現)を加えると,最初のものが目的関数として設定されるようです.
次に,問題に制約を加えます.
LP += 2 * x <= 4
LP += x + 2 * z <= 8
LP += 3 * y + z <= 6
続けて問題に線形不等式(あるいは等式)を加えると,制約として追加できるようです. これで問題の設定は完了です.
PuLPの問題のinstanceを,Pythonのprintで表示してみます.
print(LP)
LP_example: MAXIMIZE 3*x + 4*y + 2*z + 0 SUBJECT TO _C1: 2 x <= 4 _C2: x + 2 z <= 8 _C3: 3 y + z <= 6 VARIABLES x Continuous y Continuous z Continuous
どうやらこういう形式で表示されるようです.
問題を設定(あるいは入力)できたので,最適解を求めてみます. 最適化のメソッドはsolve()です.
LP.solve()
1
正常に解けると,1が返ってくるみたいです. 問題の設定によっては「解がない」「解はあるけれど際限なく最大化(あるいは最小化)できる」ということもあるので,その場合には別の値が返ってくることでしょう.
最後に,最適解を確認します.
print(x.value(), y.value(), z.value())
2.0 1.0 3.0
地味ですが,最適解を確認できました. 最適値は最適解から簡単に得られます.
次は以下の整数最適化問題 $$ \text{max. } 16 x_a + 19 x_b + 23 x_c + 28 x_d $$ $$ \text{s. t. } 2 x_a + 3 x_b + 4 x_c + 5 x_d \le 7, $$ $$ x_a, \ x_b, \ x_c, \ x_d \in \{0, 1\} $$ を解いてみる. この問題は俗に0-1ナップサック問題とよばれるものの一例である.
今度は問題のinstanceをknapという名前にする. 問題と変数のinstanceを作る.
knap = pulp.LpProblem(name='01knapsack', sense=pulp.LpMaximize)
x_a = pulp.LpVariable(name='x_a', cat='Binary')
x_b = pulp.LpVariable(name='x_b', cat='Binary')
x_c = pulp.LpVariable(name='x_c', cat='Binary')
x_d = pulp.LpVariable(name='x_d', cat='Binary')
メソッドLpVariableのパラメータcatは変数種類である. 連続変数はContinuous,整数変数はInteger,Binary変数はBinaryである. 省略するとデフォルト値のContinuousが適用されるみたいである. よって先ほどの線形最適化問題の場合には上手くいっていた.
次に制約を加えて問題を表示する.
knap += 16 * x_a + 19 * x_b + 23 * x_c + 28 * x_d
knap += 2 * x_a + 3 * x_b + 4 * x_c + 5 * x_d <= 7
print(knap)
01knapsack: MAXIMIZE 16*x_a + 19*x_b + 23*x_c + 28*x_d + 0 SUBJECT TO _C1: 2 x_a + 3 x_b + 4 x_c + 5 x_d <= 7 VARIABLES 0 <= x_a <= 1 Integer 0 <= x_b <= 1 Integer 0 <= x_c <= 1 Integer 0 <= x_d <= 1 Integer
さらに最適化して,得られた最適解も確認する.
knap.solve()
print(x_a.value(), x_b.value(), x_c.value(), x_d.value())
1.0 0.0 0.0 1.0
問題Technobubbleの説明は https://code.google.com/codejam/contest/11254486/dashboard#s=p2 にあるので,説明は省略する. ものすごくざっくり言うと,2単語からなるトピックがたくさん与えられるので偽物の最大値を答えよ,という問題である.
問題のSampleの1つ目は
HYDROCARBON COMBUSTION,
QUAIL BEHAVIOR,
QUAIL COMBUSTION,
という3つのトピックである.
$p$番目のトピックが偽物であるとき1,そうでないとき0となる変数$x_p$を用意する. すると,問題の目的はすべてのトピック$p$に関する$x_p$の和の最大化である. 問題の設定から,1つ目の単語のそれぞれに関して,それを使っているトピックのうち少なくとも1つは本物でなければならない. これを変数$x_p$を用いて表すと,1つ目の単語$w$それぞれに関して,
($w$を使っているトピック$p$に関する$x_p$の和)$\le$($w$を使っているトピックの数)$-1$
となる. 2つ目の単語に関しても同様である. これらが,変数$x_p$が満たすべき制約となる.
上記の具体的な問題例を上記と同様の整数最適化問題として定式化すると, $$ \text{max. } x_1 + x_2 + x_3 $$ $$ \text{s. t. } x_1 \le 0, \quad \text{HYDROCARBONは,1つ目の単語として,トピック1のみで使われているから} $$ $$ x_2 + x_3 \le 1, \quad \text{QUAILは,1つ目の単語として,トピック2と3で使われているから} $$ $$ x_1 + x_3 \le 1, \quad \text{COMBUSTIONは,2つ目の単語として,トピック1と3で使われているから}$$ $$ x_2 \le 0, \quad \text{BEHAVIORは,2つ目の単語として,トピック1のみで使われているから}$$ $$ x_1, \ x_2, \ x_3 \in \{0,1\}$$ となる.
さらに,一般的な定式化を示すためにまず集合の記号を用意する.
1つ目の単語の集合を$W_1 = \{w_{11},w_{12},\dots,w_{1L}\}$とし,2つ目の単語の集合を$W_2 = \{w_{21},w_{22},\dots,w_{2M}\}$とする. ここで$L, M$はそれぞれ1つ目の単語の数,2つ目の単語の数である.
トピックの集合を$P = \{p_1,p_2,\dots,p_N\}$とする. ここで$N$はトピックの数である. そして$p_n \in W_1 \times W_2 \ (n \in \{1,\dots,N\})$である.
さらに,単語$w \in W_1$が使われているトピックの集合を$Q_1(w) \subset P$とする. 同様に,単語$w \in W_2$が使われているトピックの集合を$Q_2(w) \subset P$とする.
以上の記号の準備のもとに,問題Technobubbleを整数最適化問題として記述すると $$ \text{max. } \sum_{p \in P} x_p$$ $$ \text{s. t.} \ \sum_{p \in Q_1(w)} x_p \le |Q_1(w)| - 1 \quad (w \in W_1) $$ $$\sum_{p \in Q_2(w)} x_p \le |Q_2(w)| - 1 \quad (w \in W_2) $$ $$x_n \in \{0,1\} \quad (n \in \{1,\dots,N\})$$ となる.
では,この定式化に基づいて,PuLPでTechnobubbleを解いてみる.
まず,Sampleの1つ目をtupleのlistで用意する.
topics1 = [('HYDROCARBON', 'COMBUSTION'),
('QUAIL', 'BEHAVIOR'),
('QUAIL', 'COMBUSTION')]
次に,トピックの集合を入力として,上記の定式化に基づいて解を出力する関数solveを以下に定義する.
def solve(topic):
'''
まず,定式化におけるW_1, W_2, Q_1, Q_2を作る.
'''
first_word = {} # 1つ目の単語の集合W_1に対応する辞書
second_word = {} # 2つ目の単語の集合W_2に対応する辞書
topic_using_first_word = {} # 1つ目の単語wを使っているトピックを引き出すための辞書,Q_1に対応
topic_using_second_word = {} # 2つ目の単語wを使っているトピックを引き出すための辞書,Q_2に対応
for i in range(len(topic)):
(w1, w2) = topic[i]
first_word[w1] = True # トピックに現れたら辞書first_wordのkeyとして登録
second_word[w2] = True # トピックに現れたら辞書second_wordのkeyとして登録
if w1 not in topic_using_first_word:
topic_using_first_word[w1] = [i] # 1つ目の単語w1が使われているトピックを覚える
else:
topic_using_first_word[w1].append(i) # 1つ目の単語w1が使われているトピックを覚える
if w2 not in topic_using_second_word:
topic_using_second_word[w2] = [i] # 2つ目の単語w1が使われているトピックを覚える
else:
topic_using_second_word[w2].append(i) # 2つ目の単語w1が使われているトピックを覚える
'''
ここまでで,定式化に必要な集合を作り終えたので,以下でPuLPを用いて解く.
'''
problem = pulp.LpProblem(name='max fake', sense=pulp.LpMaximize) # 最大化問題として設定
x = {}
for i in range(len(topic)):
x[i] = pulp.LpVariable(name='t'+str(i), cat='Binary') # それぞれの変数x_p(ここでの名前は'tp')をBinary変数として設定
problem += pulp.LpAffineExpression([(x[i], 1) for i in range(len(topic))]) # ここで目的関数を設定,メソッドLpAffineExpressionの使い方の説明は割愛
for w1 in first_word: # 1つ目の単語のそれぞれに関して制約を加える
problem += pulp.LpAffineExpression([(x[i], 1) for i in topic_using_first_word[w1]]) <= len(topic_using_first_word[w1]) - 1
for w2 in second_word: # 2つ目の単語のそれぞれに関して制約を加える
problem += pulp.LpAffineExpression([(x[i], 1) for i in topic_using_second_word[w2]]) <= len(topic_using_second_word[w2]) - 1
problem.solve() # 最適化
return sum([x[i].value() for i in range(len(topic))]) # 最適解から最適値を計算,リスト内包表記はこういうときに便利
関数を定義したら,先ほどのサンプルで解いてみる.
int(solve(topics1)) # 最後の出力はGoogle code jamの形式に合わせて整数で
1
ここではSampleの1つ目を解くだけにとどめるが,もちろんLargeも実用時間内に解ける. 宮本の手元のパソコンで解いたら計算時間は2.5秒くらいであった.
一般には,大規模な整数最適化問題の最適解の発見(保証)には多大な時間がかかる. しかし,この問題Technobubbleを上記の方法で定式化したものは完全ユニモジュラー(Totally Unimodular)という性質を満たすので,高速に解ける. 完全ユニモジュラー性の証明は自分で考えてみよう!(わかっている人には自明かもしれない.) この性質から,実質的には線形最適化で最適解が見つかっている.
数理最適化モジュールをただ使えるというだけでなく,「上手く使える」ようになるためには完全ユニモジュラーのような概念の修得・理解が大事である. 完全ユニモジュラーを含む数理最適化の理論的側面は,大学院生向けの講義「数理最適化特論」で扱っている. 興味があったらいずれ履修してみて欲しい.
最後に,今更言うことでもないが,Google code jamで数理最適化ソルバーを用いるのは反則である.