本稿では、最適化モデリングツールfloptの基本的な使い方やいくつかの機能の具体例を紹介します。
最適化モデリングツールとは、ユーザーが解きたい問題を表現、具現化する作業をサポートするソフトウェアです。
floptでは
- (制約付き) (非) 線形計画問題
- ブラックボックス最適化
- 充足最大化問題
- 巡回セールスマン問題のような順列を変数として持つ問題
- 二次計画問題
などの幅広い問題を表現できる柔軟性が見どころです。
本稿ではflopt:version 0.5.3を用いて説明を行います。
はじめに
最適化モデリングツールについて、線形計画法のモデリングであればPulpやScipy、二次計画であればCVXOPTやPyomoなどが知られています。
当然ですが、これらのツールはそれぞれでモデルリングのプログラムの書き方が異なります。例えばScipy.optimize.milpでMIP(混合整数線形計画法)を解こうとすると、ユーザーが不等式標準形で問題を表現し、そのときの制約行列とベクトルを関数に与える必要があり、これは非常に手間です。
また一般に、ある定式化は他の最適化問題での表現が可能である場合が多いです（例えばイジングモデルは二次計画やQUBOでの定式化へ書き換えができるます) 。そのため、求解に使用できるソルバはその定式化に応じて複数の選択肢があります。しかし、使うソルバに応じてモデリングプログラムをユーザーが(書き方を調べながら)作成する必要があり、これも大変手間でしょう。
floptはそのような問題を解決するために、様々な定式化について統一的にモデリングが行えるツールを目指し、Pythonライブラリとして開発されています。ユーザーはfloptで最適化問題をモデリングするだけで、他のモデリングツールへの変換や複数のソルバを簡単に使用することができます。
現在はScipy(linprog, milp, minimize)やOptuna、Hyperopt、Pulp, CVXOPTなどの外部ライブラリのソルバやランダムサンプリングや2-Optなどの内部実装のアルゴリズムを選択できます。
floptの特徴
- 非線形式、線形式、多項式、ブラックボックス関数をサポートする柔軟な数式表現能力
- 線形計画法において、(不)等式標準形への変換や変数積の自動線形化機能
- 他の問題表現形式への変換機能 (e.g. イジングモデル → 二次計画)
- 外部ライブラリを含む複数のソルバの呼び出しや、外部モデリングツールへの変換が可能
- Pulpライクな直感的なユーザーインターフェイス
他の問題表現形式へ変換したり、変数積を自動線形化する機能は他のモデリングツールにはみられないものです。
本稿ではfloptの基本的な使い方や上記の特徴の具体例を紹介します。
インストール
pipを用いてインストールできます。
pip install flopt
一般的な使い方
ここでは、下記の制約付き非線形問題のモデリング + 最適化をしてみます。
二つの変数a, bからなる非線形の目的関数と制約から構成されています。
minimize 2*(3*a+b**2)+3
s.t. a*b >= 2
a+b >= 2
0 <= a <= 1, a is integer
1 <= b <= 2, b is continuous
import flopt
# create variables
a = flopt.Variable("a", lowBound=0, upBound=1, cat="Continuous")
b = flopt.Variable("b", lowBound=1, upBound=2, cat="Continuous")
# create problem
prob = flopt.Problem(name="Test")
prob += 2*(3*a+b**2)+3 # set the objective function
prob += a*b >= 2 # add constraint
prob += a+b >= 2 # add constraint
# create solver
solver = flopt.Solver(algo="ScipySearch")
# run solver
prob.solve(solver, timelimit=10)
# get best solution
print("obj value =", prob.getObjectiveValue())
print("a =", a.value())
print("b =", b.value())
Pulpやその他のモデリングツールとは異なり、変数同士の掛け算や累乗計算を柔軟に表現できていますね。
Problemオブジェクトに+=で目的関数や制約を追加しています。(それぞれ .setObjective(), .addConstraint()関数を用いた操作も可能です）追加されるのが目的関数なのか制約なのかは式に(==, <=, >=)が含まれているかで判断されます。
また上の例ではScipySearchソルバを使用しています。内部ではscipy.optimize.minimizeソルバが実行されていて、これは制約付きの非線形問題を扱えるソルバです。
floptで使用できるソルバ
floptではいくつかの外部ソルバや内部で実装されたアルゴリズムをSolverとして実行できるようになっています。サポートしているSolverの一覧はここから確認できます。主にはScipyやOptuna、Hyperopt、Pulp, CVXOPTなどの外部ライブラリやランダムサンプリングや2-Optなどの内部実装のアルゴリズムを使用できます。さらに、
solver = flopt.Solver(algo="auto")
とalgoに"auto"を指定すると、おすすめのsolverを自動で選択してくれます。どのsolverが選択されているのかはsolver.select(prob).nameで確認できます。
選出されるsolverは問題と時間制限(timelimit)によって変わります。例えば下記の無制約問題の場合、
a = flopt.Variable("a", 0, 1, cat="Continuous")
b = flopt.Variable("b", 1, 2, cat="Continuous")
prob = flopt.Problem(name="Test")
prob += 2*a + 3*b # set the objective function
timelimit=1を設定した場合
solver = flopt.Solver(algo="auto")
solver.setParams(timelimit=1)
solver.select(prob).name
>>> 'RandomSearch'
timelimit=10を設定した場合
solver = flopt.Solver(algo="auto")
solver.setParams(timelimit=10)
solver.select(prob).name
>>> 'OptunaCmaEsSearch'
と選択されているソルバが変わっていますね。
ケーススタディー
ケーススタディーとしてドキュメントでは
についての実装方法を紹介しています。ぜひご覧ください。
問題の表現形式の変換
最適化問題ではさまざまな標準的な表現形式があります。 例えば下記はQUBOと呼ばれる問題形式です。
minimize x^T Q x + C
x_i is in {0, 1}
また量子アニーリング系だとイジング形式の定式化が使用されることが多いです。
minimize - x^T J x - h^T x + C
x_i is in {-1, 1}
そのほかにも Linear Programming (LP)やQuadratic Programming (QP)などの表現形式があります。floptではこれら任意の二つの表現形式間の変換を下図の変換グラフを用いて行います。
例えばLPを表現するために使用するLpStructureからQUBOを表現するのに使用するQuboStructunreへ変換するためには下図の順に変換が行われます。
このように連鎖的に問題形式を変換することで、任意の問題変換を実現しています。
ここではIsingからQUBOへの変換を見てみましょう。 次のようなIsing問題を例にします。
minimize 1-a*b-a
s.t. a, b is in {-1, 1}
floptでのモデリング。取る値が-1, 1の変数はcat="Spin"を指定して生成します、これをスピン変数と呼ぶことにします。
import flopt
# Variables
a = flopt.Variable("a", cat="Spin")
b = flopt.Variable("b", cat="Spin")
# Problem
prob = flopt.Problem()
prob += 1 - a * b - a
IsingStructureに変換したあと、.toQubo()でQuboStructureへ変換します。
from flopt.convert import IsingStructure
ising = IsingStructure.fromFlopt(prob)
qubo = ising.toQubo() # convert ising to QUBO
print(qubo.toFlopt().show())
>>> Name: None
>>> Type : Problem
>>> sense : minimize
>>> objective : -4.0*(a_b*b_b)+(2.0*b_b)+1.0
>>> #constraints : 0
>>> #variables : 2 (Binary 2)
変換結果を見やすくするために、あえて.toFlopt()でfloptのProblemオブジェクトに変換してから表示しています。
スピン変数a, bがそれぞれバイナリ変数a_b, b_bに変換されています。同様に線形計画問題やQPをIsingやQUBOもしくはその逆に変換したり自由に問題を変換することができます。
問題の線形化
floptには、目的関数や制約に変数積がある場合に、それらを線形に自動で変換する機能があるので紹介します。
:::note note
変数積を線形化するにはもちろん条件があって、基本的には変数に下限上限制約がついている必要があります。条件さえ揃えば整数変数同士や、整数変数と連続変数の変数積などを線形式に変換することができます。詳細はhttp://web.tuat.ac.jp/~miya/fujie_ORSJ.pdf の3章を参考にしてください。
:::
ここでは下記の問題を例にして線形化を行ってみます。
minimize a * b * c + 2
s.t. a * c == 2
0 <= a <= 1, a is integer
1 <= b <= 2, b is integer
1 <= c <= 3, c is continuous
まずは普通にfloptで問題を生成したのち
import flopt
a = flopt.Variable(name="a", lowBound=0, upBound=1, cat="Integer")
b = flopt.Variable(name="b", lowBound=1, upBound=2, cat="Integer")
c = flopt.Variable(name="c", lowBound=1, upBound=3, cat="Continuous")
prob = flopt.Problem()
prob += a * b + c + 2
prob += a * c == 2
flopt.convert.linearizeを実行するだけ。
flopt.convert.linearize(prob)
print(prob.show())
>>> Name: None
>>> Type : Problem
>>> sense : minimize
>>> objective : mul_0+2*mul_1+c+2
>>> #constraints : 15
>>> #variables : 10 (Continuous 2, Integer 2, Binary 6)
>>>
>>> C 0, name None, mul_2-2 == 0
>>> C 1, name for_bin_a_sum, bin_a_0+bin_a_1-1 == 0
>>> C 2, name for_bin_a_eq, a-bin_a_1 == 0
>>> C 3, name for_bin_b_sum, bin_b_0+bin_b_1-1 == 0
>>> C 4, name for_bin_b_eq, b-bin_b_0-(2*bin_b_1) == 0
>>> C 5, name for_mul_0_1, mul_0-bin_a_1 <= 0
>>> C 6, name for_mul_0_2, mul_0-bin_b_0 <= 0
>>> C 7, name for_mul_0_3, mul_0-(bin_a_1+bin_b_0-1) >= 0
>>> C 8, name for_mul_1_1, mul_1-bin_a_1 <= 0
>>> C 9, name for_mul_1_2, mul_1-bin_b_1 <= 0
>>> C 10, name for_mul_1_3, mul_1-(bin_a_1+bin_b_1-1) >= 0
>>> C 11, name for_mul_2_1, mul_2-bin_a_1 >= 0
>>> C 12, name for_mul_2_2, mul_2-(3*bin_a_1) <= 0
>>> C 13, name for_mul_2_3, mul_2-(c-(1-bin_a_1)) >= 0
>>> C 14, name for_mul_2_4, mul_2-(c-(3*(1-bin_a_1))) <= 0
新たに変数が追加されていることに気がつくと思います。内部では整数制約をバイナリ変数に分解したり、変数xとyの積を表す変数z=xyを追加したり、という作業を経て線形化を行っています。あとは適当にsolverを指定してsolve関数を実行すれば求解されます。
solver = flopt.Solver("auto")
prob.solve(solver, msg=True)
線形計画問題の標準形変換
さて、ここでは作成した線形計画問題を標準形に変換する機能を紹介します。
次のようなシンプルな線形計画問題を考えます。
minimize a + b + c + 2
s.t. a + b == 2
b - c <= 3
0 <= a <= 1, a is integer
1 <= b <= 2, b is continuous
1 <= c <= 3, c is continuous
floptでのモデリング。
import flopt
# variables
a = flopt.Variable(name="a", lowBound=0, upBound=1, cat="Integer")
b = flopt.Variable(name="b", lowBound=1, upBound=2, cat="Continuous")
c = flopt.Variable(name="c", lowBound=1, upBound=3, cat="Continuous")
# problem
prob = flopt.Problem(name="LP")
prob += a + b + c + 2
prob += a + b == 2
prob += b - c <= 3
これをLpStructureに変換します。
from flopt.convert import LpStructure
lp = LpStructure.fromFlopt(prob)
print(lp.show())
すると、probを下記の形で表現したときの行列G、Aやベクトルh, bなどがLpStructureオブジェクトに格納され、それらはLpStructureの属性として取り出すことができます(lp.A)など。
>>> LpStructure
>>> obj c.T.dot(x) + C
>>> s.t. Gx <= h
>>> Ax == b
>>> lb <= x <= ub
また
lp = lp.toAllEq()
を実行すると、全ての制約が等式になるように変換してくれます。いわゆる等式標準形です。上記の問題の例だと
print(lp.toAllEq().toFlopt().show())
>>> Name: None
>>> Type : Problem
>>> sense : minimize
>>> objective : a+b+c+2
>>> #constraints : 2
>>> #variables : 4 (Continuous 3, Integer 1)
>>>
>>> C 0, name None, a+b-2.0 == 0
>>> C 1, name None, -c+b+__s_0-3.0 == 0
変換結果を見やすくするために、あえて.toFlopt()でLpStructureをfloptのProblemオブジェクトに変換してから表示しています。
等式制約に変換するためのスラック変数__sがfloptにより追加されていますね。
また全ての制約を不等式形にするには.toAllNeq()を使用します。
print(lp.toAllNeq().toFlopt().show())
>>> Name: None
>>> Type : Problem
>>> sense : minimize
>>> objective : a+b+c+2
>>> #constraints : 3
>>> #variables : 3 (Continuous 2, Integer 1)
>>>
>>> C 0, name None, -c+b-3.0 <= 0
>>> C 1, name None, a+b-2.0 <= 0
>>> C 2, name None, -a-b+2.0 <= 0
a+b-2.0 == 0という制約がa+b-2.0 <= 0と-a-b+2.0 <= 0の二つの不等式制約に変換されています。
開発の経緯
もともとfloptはメタヒューリスティック手法を実装し、その評価を行うために開発していました。メタヒューリスティックは目立った無償のモデリングツールが存在していおらず、PuLPのようなユーザーインターフェイスのツールがあると便利そう、ということで開発を進めていきました。そのため、非線形の式も表現できるようにしたり、ベンチマーク問題を複数のソルバで解いてそのパフォーマンスの比較を行えるような機能も実装されています (ここを参照)。
その過程でScipyの非線形計画ソルバやOptunaソルバの性能比較もしてみたくなりソルバに追加されました。 そうするとどんどん他のソルバも使えるようにしたくなり、イジングモデルや線形計画ソルバも追加し、またそれらのモデル間の変換機能も実装していきました。変数積も自動で線形化できたら便利だよね、というのでその機能も追加されました。
floptの基本的な開発が終わったらメタヒューリスティックを実装しようと思っていましたが、最近は問題変換機能の方に力を入れており、現状ではそちらの方に見どころのあるツールになっているかなと思います。
まとめ
本稿ではfloptを用いたモデリングについて実装例、および、表現形式間の変換、線形計画法の自動変数積の線形化機能と標準形変換機能を紹介しました。floptはGithub上で公開されていますので、ご意見ご要望はissueや本稿のコメントに投稿もらえると嬉しいです！