MyEnigma

とある自律移動システムエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

Juliaにおける最適化モデリングライブラリJuMP入門

目次

 

はじめに

以前、Python製の最適化モデリングライブラリ

cvxpyを紹介しましたが、

myenigma.hatenablog.com

 

今回は、高速科学演算向け言語であるJuliaの

最適化モデリングツールであるJuMP (Julia for Mathematical Optimization)

の紹介をしようと思います。

JuMP — Julia for Mathematical Optimization — JuMP -- Julia for Mathematical Optimization 0.18 documentation

github.com

 

JuMPの詳細に関しては、下記の論文を参照ください。

 

プログラミング言語Juliaに関しては下記の記事を参照ください。

myenigma.hatenablog.com

 

JuMPの特徵

下記が公式ドキュメントに書かれている

JuMPの特徴です。

使いやすさ

自然な数式に似たシンタックスで

最適化問題を記述できるため簡単に利用できる。

またコードがドキュメントとしても機能する。

 

計算速度

ベンチマーク結果によると、

JuMPの最適化モデリングの速度は、

最適化モデリング専用の言語であるAMPLなどと

同じ位の速度を有しています。

また、JuMPは各ソルバーとのデータのやり取りを

ファイル作成ではなく、メモリ上でやり取りすることにより、

高速に計算できるようになっています。

 

また、一般的なモデリングツールは、

プログラミング言語のオペレータオーバロード機能を使って、

モデリングを実施しますが、

JuMPはJuliaのマクロ機能を使っているため、

高速に最適化モデリングを実現できます。

 

ソルバーの独立利用のためのインターフェース

JuMPは一般的なソルバーを独立的に利用できるように

インターフェースが設計されているため、

オープンソースのソルバーから、

商用ソルバーまで様々なソルバーを簡単に切り替えて利用することができます。

現在使用できるソルバーとしては、

上記のソルバーの中には、

特別な入力を必要とするソルバーもありますが

JuMPはその違いを吸収できます。

 

それぞれのソルバーがどのような最適化問題に使えるのかは、

下記のリンクを参照ください。

http://jump.readthedocs.io/en/latest/installation.html

 

また、この設計により、

混合整数問題でしばしば利用されるbranch-and-bound法

などで必要な最適化中のソルバーとの情報のやりとりが可能になり、

最適化中に動的に最適化問題を変更することが可能であることも、

JuMPの特徴です。

 

システムに組み込みやすい

JuMPはJuliaで書かれており、

ソルバーはバイナリのものを使うように設計されています。

したがって、他の言語からJuMPのコードを呼んだり、

大きなソフトウェア(シミュレーションやWebサーバなど)

にJuMPのコードを組み込むのも簡単です。

 

ライセンスがMPL

JuMPはライセンスがMozilla Public License(MPL)であるため、

商用利用も可能です。

 

インストール方法

JuliaのREPLで、

julia> Pkg.add("JuMP")

とすることでJuMPがインストールされます。

 

続いて、それぞれのソルバーをインストールします。

GLPKのインストール

julia> Pkg.add("GLPKMathProgInterface")

Clpのインストール

julia> Pkg.add("Clp")

Cbcのインストール

julia> Pkg.add("Cbc")

Culpのインストール

julia> Pkg.add("Culp")

ECOSのインストール

julia> Pkg.add("ECOS")

GitHub - JuliaOpt/ECOS.jl: Julia wrapper for the ECOS conic optimization solver

 

Ipoptのインストール

julia> Pkg.add("Ipopt")

 

NLoptのインストール

julia> Pkg.add("NLopt")

 

Mosekのインストール

Mosekをインストールした後、

julia> Pkg.add("Mosek")

julia> Pkg.build("Mosek")

でOKです。

 

CPLEXのインストール

CPLEXのインストールは下記の記事を参考にして、

CPLEXをインストールし、

github.com

 

bashrcに下記のようにCPLEXのパスを指定したあと、

export LD_LIBRARY_PATH="/Users/hoge/Applications/IBM/ILOG/CPLEX_Studio1271/cplex/bin/x86-64_osx":$LD_LIBRARY_PATH

下記のコマンドでインストールします。

$CPLEX_STUDIO_BINARIES=/path/to/cplex/bin/x86-64_linux julia -e 'Pkg.add("CPLEX"); Pkg.build("CPLEX")'

 

各ソルバーのパラメータ設定

各ソルバーで様々なパラメータを設定できますが、

よく使うものをまとめておきます。

最適化の標準出力を止める

CPLEXの場合

using JuMP, CPLEX
m = Model(with_optimizer(Cplex.Optimizer, CPX_PARAM_SCRIND=0))

Gurobiの場合

using JuMP, Gurobi
m = Model(with_optimizer(Gurobi.Optimizer, OutputFlag=0))

Ipoptの場合

using JuMP, Ipopt
m = Model(with_optimizer(Ipopt.Optimizer, print_level=0))

  

最適化計算時間に上限を設定する

それぞれ1時間(3600秒)を上限としてみます。

 

CPLEXの場合

using JuMP, CPLEX
m = Model(with_optimizer(Cplex.Optimizer, CPX_PARAM_TILIM=3600))

Gurobiの場合

using JuMP, Gurobi
m = Model(with_optimizer(Gurobi.Optimizer, TimeLimit=0))

その他のパラメータ

IBM Knowledge Center

www.gurobi.com

  

JuMPによるシンプルな線形計画問題の解法サンプルコード

下記の記事を参照ください。

myenigma.hatenablog.com

  

非線形最適化のサンプルコード

cvxpyは凸最適化のみに対応しているため、

非線形最適問題は解くことができませんが、

JuMPの場合は、非線形最適化のモデリングも可能です。

 

具体的にJuMPを使って、

下記の非線形最適化問題を解いてみたいと思います。

f:id:meison_amsl:20170802081157p:plain

 

下記が実際のJulia + JuMPを使ったサンプルコードです。

今回は非線形ソルバーであるIpoptを使用しています。

非線形最適化問題であるため、

初期値を変えると、最適値が異なっていることが

実行時の出力から確認できます。

 

using JuMP
using Ipopt

function main_ipopt()
    println("JuMP sample")
    model = Model(with_optimizer(Ipopt.Optimizer))

    @variable(model, z1, start= -0.3)
    @variable(model, z2, start= 0.5)
    @NLobjective(model, Min, 3.0*sin(-2*pi*z1)+2*z1+4+cos(2*pi*z2)+z2)
    @constraint(model, -1<=z1<=1)
    @constraint(model, -1<=z2<=1)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    optimize!(model)
    println("Objective value: ", objective_value(model)) 
    println("z1:",value(z1))
    println("z2:",value(z2))
end

function main_ipopt2()
    println("JuMP sample")
    model = Model(with_optimizer(Ipopt.Optimizer))

    @variable(model, z1)
    @variable(model, z2)
    @NLobjective(model, Min, log(1+z1^2)-z2)
    @NLconstraint(model, -(1+z1^2)^2+z2^2==4)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    optimize!(model)
    println("Objective value: ", objective_value(model)) 
    println("z1:",value(z1))
    println("z2:",value(z2))
end


main_ipopt()
main_ipopt2()

github.com

 

JuMPの最適化コードを静的コンパイルして高速化する方法

下記の記事の方法を使うと、

シンプルな最適化コードであれば

Juliaのコードをコンパイルして高速実行できます。

myenigma.hatenablog.com

 

参考資料

github.com

Jupyter Notebook Viewer

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com