Python + Pyomoによる(非線形)数値最適化
TL; DR
- Pythonのデータマネジメント技術と数値最適化をスムーズに繋げたい
- Pyomoを使うことで自然な記法でモデルを組み立てることが出来る
- Webドキュメントは貧弱だが、コミュニティは活発!
概要
数値最適化は機械学習や数値モデリングの基礎も基礎ですが、単に役に立ちます。Pythonという観点を度外視すれば、いろいろな手段がありますが、基本的には共通するフローがあります。
- 現象から最適化したい目的関数を決定する
- 目的関数に付随する制約条件を立式化する
- データを用意し、モデリング言語を記述する
- モデルをソルバーへ投入し、解を得る
こういった作業をするための環境は商用でも色々と売られていまして、例えばAIMMS社とか、GAMS社、AMPL社といった会社は最適化の開発環境や言語を作って販売しています。
もちろん非商用でも、こういう団体は存在し、COIN-ORではライブラリやソルバーを開発して公開しています。商用のソフトウェアに比べれば少々パフォーマンスに劣りますが無料で使えるということを考えれば文句も出ません。
しかしながら、例えばAMPL言語の書き方を見てみると…ちょっとめんどくさそうです。Stanかよ。新しく言語を覚えるのは大体の場合面倒くさいので、実務でちゃちゃっと最適化したい場合には向きません。どうせだったらPythonで済ませられたらうれしいです。
Pythonで最適化?
Pythonで最適化する方法も色々あります。
Scipy.optimize
一番最初に思いつく方法だと思います。アルゴリズムも色々ありますし、導入も簡単です。ただし、簡単な最小化とか線形計画法(LP)に特化しています。もうちょっと痒いところに手が伸びるものが欲しいのです。
PuLP
日本語で調べるとこのライブラリについての記事が一番多く出てきます。PuLPは先程のCOIN-ORによって開発が進められているライブラリです。使いやすいですし、Sphinxでのドキュメントも有るのが個人的にはありがたいところです。
ただし自分が使う上では、
- 使えるソルバーが少ない
- 最近開発が停滞気味
といった2点が気になりました。
リポジトリに書いてあるように、PuLPはLP/MILP用ライブラリです。デフォルトソルバーはCOIN-OR CBCで、他に使えるのもCPLEX、GLPKなどLPソルバーです。非線形問題は解けません。
もう一つ、PuLPは2017年の中頃からレポジトリが更新されていません。ドキュメントも更新がされておらず、実際に使う上で良く分からないことも多いです。というか英語ドキュメントすらなさすぎて、日本語記事のほうが詳しいかも…。先行記事では
などがありました。他にもQiitaやBlog上に、沢山の詳しい記事があります。
CVXPY
こちらの記事がとても詳しいです*1。
Pyomo
PyomoはOpen-OR のウェブページにも載っていますが、メインで開発しているのは米サンディア国立研究所のようです。昔はCooprという名前でした。Pythonのオブジェクト指向な書き方をとても重視したライブラリです。しかし僕はPuLPの次にこれを使ってみたら、結構似ていると感じたので、移行も簡単でしょう。PyomoやPuLPは2008年ぐらいに沢山開発されたPythonライブラリで、今も生き残っているものたちです。なので、書きやすさは担保されていると思います(参考↓)。
Pyomoの嬉しいところは他に2つあって、
- 使えるソルバーが多い
- コミュニティがそこそこ活発
という点です。PyomoはNLフォーマットを利用するソルバー(AMPL Solver Library)に対応しています。PuLPで使えていたCBC、GLPK、CPLEXなども利用することが出来ますし、その他についてはAMPLのページを見たほうが早いですね。PuLPでは使えなかった非線形ソルバーで有名なSCIPが使えますし、Open-OR主導で開発されているBonminやCouenne、IPOPTなども利用できます。ありがたいものです。
もう一つ、PyomoはPuLPと同じぐらいドキュメントは弱いのですが*2、コミュニティは活発です。リポジトリは頻繁にコミットがありますし、チュートリアル用のリポジトリやモデル例もあります。個人的に一番ありがたかったのは、Google Groupが活発で、質問にすぐに答えてくれた点です。
そういうわけで、今回はPyomoを使って最適化を試してみます。先行記事には
がありました。他に日本語の記事があまり出てこないので、残しておきます。
導入
ローカルへ構築
ライブラリ自体は
$ pip install pyomo
で終わりですが、これだとソルバーが全然入っていません。ソルバーはコマンドラインから実行できるようにする必要があります。自分はOSXの端末に以下のソルバーを入れました。皆様はお好きにどうぞ。
- glpk
-
brew install glpk
で終了!
-
- CPLEX
- IPOPT
- SCIP
- そのまま導入しただけでは利用することが出来ません。
- Google Groupではnlオプションを付ければいいと書いてありましたが、僕は無理でした。
- SCIPをAMPLフォーマットで読み書きするようにコンパイルする必要があります。これはこちらをそのまま参考にしました。
- Bonmin、Couenne
- Open-ORのソルバーなので、IPOPTと同様の手順で入ります。
- IPOPTの時にOPENBLASやLAPACKを入れていれば、飛ばしてOKです。
Dockerを利用する(2019-03-29追記)
Dockerを使ったコンテナを作りました。以下のページに詳細は書きましたので、面倒くさい方はこちらをどうぞ。
非線形最適化問題を解く
どうせなので、非線形問題を解いてみましょう。例として、放送大学の資料にあった問題を解きます。
工数の最適配分の計算 * 工程i(i=1, 2, ..., n)にx_i(日)かけたときの仕上がりの程度がlog(a_i * x_i + 1) ただしa_i > 0 で決定する * 各工程には最低でt_i(>0)日必要である * 全行程をT日以内に終わらせる必要がある * 仕上がりを最大化するには、各工程にどのように割り振ればいいか
仕上がりの程度にlogが入っているので非線形最適化問題になります。細かく言うなら、変数であるxが整数しか取らないので、整数非線形計画法(Integer Non Linear Programming; INLP)というパラダイムで解かれる問題ですね。ただし整数計画問題(Integer Programming; IP)やより狭義の整数線形計画問題(Integer Linear Programming)に特化した研究は色々ありますが、INLPに特化したものは少なく、一部の変数だけが整数条件を取る混合整数非線形計画法(Mixed INLP; MINLP)の一部として取り組むことが多いです。ソルバーもMINLP用のものを探したほうがいいです。
以下Jupyter + Pythonでの解法です。
- 2, 26行目: これはJupyter環境でやっています。Pyomoの例題ではコマンドラインから動かしていますが、こうすればJupyterでも動きます。optの部分でsolverを指定して下さい。
- 7, 9行目: PyomoはList型のIndexingを自動でやってくれないので、Dict型に変換しています。
- 14行目: 変数についてはVarで宣言し、2つめの位置変数へデータ型を入れます。今回は日数なので、非負整数です。
- 16~20行目: 制約条件を書いていきます。最適化条件と順番が前後しても大丈夫です。
- Pyomoの変な特性ですが、グローバル変数とローカル変数がごっちゃになっても受け付けてしまいます。より良い書き方がありそうですが…現状知りません。
- 21行目: 制約条件時に、Iの各要素iで成立…としたい場合は、第一位置変数へ入力すれば受け付けてくれます。
- 22~24行目: 最適化条件を書きます。senceで最大化(maximize) or 最小化(minimize)を指定します。
- 27行目: opt.solveをする時に、tee=Trueとすればソルバーのメッセージが表示されます。デバッグ時に便利です。
- 28行目: モデルの結果を表示します。
以下出力です(デバッグ部分は除きました)。
Model Days allocation Variables: x : Size=4, Index=I Key : Lower : Value : Upper : Fixed : Stale : Domain 1 : 0 : 15.0 : None : False : False : NonNegativeIntegers 2 : 0 : 20.0 : None : False : False : NonNegativeIntegers 3 : 0 : 50.0 : None : False : False : NonNegativeIntegers 4 : 0 : 15.0 : None : False : False : NonNegativeIntegers Objectives: value : Size=1, Index=None, Active=True Key : Active : Value None : True : 19.26358596079164 Constraints: const1 : Size=1 Key : Lower : Body : Upper None : None : 100.0 : 100.0 const2 : Size=4 Key : Lower : Body : Upper 1 : 15.0 : 15.0 : None 2 : 20.0 : 20.0 : None 3 : 50.0 : 50.0 : None 4 : 5.0 : 15.0 : None
4つのタスクに対して、15:20:50:15と割り振ることでパフォーマンスが最大化され、パフォーマンス和は19.26です。さっさと終わる仕事の4番目ですが、プロトタイプが出来てからの伸びしろが良いのでこうなるわけです。
おわり
以上Pythonでの数値最適化でした。皆様も是非ガシガシ最適化して、ついでにブログもアップして、勉強させて下さい。
参考
*1:解説が薄いのは、このライブラリ(記事)を見つけたのが記事を書き終わってからというだけであり、他意はありません
*2:APIが纏められていない点が使いづらいです。出版されている本に載っているようですが、高い…
*3:昨年末に書かれた記事があったので、心配はしていませんでした