MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

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

Introduction to Non-Linear Optimization (Computer Science)

Introduction to Non-Linear Optimization (Computer Science)

目次

 

はじめに

以前、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は一般的なソルバーを独立的に利用できるように

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

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

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

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

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

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

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

 

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

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

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

 

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

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”)’

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

以前、cvxpyで解いたシンプルな線形計画問題をJumpで解いてみました。

myenigma.hatenablog.com

   

using JuMP
using Mosek

function main()
    println("JuMP sample2")
    model = Model(solver=MosekSolver())

    @variable(model, 2.5 <= z1 <= 5.0)
    @variable(model, -1.0 <= z2 <= 1.0)
    @NLobjective(model, Min, abs(z1+5.0) + abs(z2-3.0))

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

main()

github.com

 

上記を実行すると下記のように、

cvxpyで同じ問題を解いた時と同じ結果が得られます。

f:id:meison_amsl:20170731012256p:plain

 

一つ注意点として、abs関数をコスト関数の中で利用する場合は、

@objectiveではなく、@NLobjectiveを使う必要があるみたいです。

stackoverflow.com

 

下記はその他のサンプルです。

f:id:meison_amsl:20170801064857p:plain

using JuMP
using Mosek

function main()
    println("JuMP sample3")
    model = Model(solver=MosekSolver())

    @variable(model, z1>=0.0)
    @variable(model, z2>=0.0)
    @objective(model, Min, 3.0*z1 + 2.0*z2)

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

main()

 

f:id:meison_amsl:20170801070633p:plain

function main3()
    println("JuMP sample5")
    model = Model(solver=MosekSolver())

    @variable(model, z1>=0.0)
    @variable(model, z2>=0.0)
    @objective(model, Min, z1-z2)
    @constraint(model, z1-z2>=30)
    @constraint(model, 2.0*z1+z2>=1)

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

main3()

 

f:id:meison_amsl:20170801071107p:plain

function main5()
    println("JuMP sample7")
    model = Model(solver=MosekSolver())

    @variable(model, z1>=1.0)
    @variable(model, z2>=1.0)
    @objective(model, Min, z1^2+z2^2)

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

 

f:id:meison_amsl:20170801072757p:plain

function main8()
    println("JuMP sample8")
    model = Model(solver=MosekSolver())

    @variable(model, z1>=0.0)
    @variable(model, z2>=0.0)
    @variable(model, z3>=0.0)
    @objective(model, Min, 0.5*(z1^2+z2^2+0.1*z3^2)+0.55*z3)
    @constraint(model, z1+z2+z3==1)

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
    println("z3:",getvalue(z3))
end

 

JuMPによる最小二乗法解法サンプルコード

以前、cvxpyで解いた最小二乗法問題をJumpで解いてみました。

myenigma.hatenablog.com

 

using JuMP
using Mosek

function main()
    println("JuMP sample1")
    model = Model(solver=MosekSolver())
    m = 10
    n = 5
    A = rand(m, n)
    b = rand(m, 1)
    println(A)
    println(b)

    @variable(model, 0.0<=x[1:size(A,2)]<=1.0)
    @objective(model, Min, sum((A*x-b).^2))

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("x:",getvalue(x))
end

main()

github.com

すると下記のように最小二乗法の計算結果が得られます。

modelをprintすると、

最適化問題が数式みたいに表示されるのはいい感じですね。

ただ行列をprintした時に、

行と列で綺麗に表示されないのが悲しい。。。

f:id:meison_amsl:20170730084633p:plain

f:id:meison_amsl:20170730084746p:plain

 

JuMPによる混合整数最適化問題の解法

続いて、

下記の2つの混合整数(Mixed integer)最適化問題をJuMPを使って解いてみます。

f:id:meison_amsl:20170803100420p:plain

f:id:meison_amsl:20170803100435p:plain

 

下記がJuliaのサンプルコードです。

今回はMixed integerを解くことができる。

MosekとCPLEXをソルバーとして使っています。

using JuMP
using Mosek
using CPLEX

function main1()
    println("JuMP mix integer sample1")
    model = Model(solver=MosekSolver())

    @variable(model, z1, Int)
    @variable(model, z2, Int)
    @objective(model, Min, -6*z1-5*z2)
    @constraint(model, z1+4*z2<=16)
    @constraint(model, 6*z1+4*z2<=28)
    @constraint(model, 2*z1-5*z2<=6)
    @constraint(model, 0<=z1<=10)
    @constraint(model, 0<=z2<=10)

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

function main2()
    println("JuMP mix integer sample2")
    model = Model(solver=CplexSolver())

    @variable(model, z1>=0)
    @variable(model, z2>=0)
    @variable(model, z3, Bin)
    @objective(model, Min, -z1-2*z2)
    @constraint(model, 3*z1+4*z2-12<=4*z3)
    @constraint(model, 4*z1+3*z2-12<=4*(1-z3))

    println("The optimization problem to be solved is:")
    println(model)

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
    println("z3:",getvalue(z3))
end

main1()
main2()

github.com

 

JuMPによるナップサック問題の解法サンプルコード

以前cvxpyで解いたナップサック問題を

JuMPで解いてみました。

myenigma.hatenablog.com

 

using JuMP
using Mosek

function main()
    println("JuMP knapsack")
    model = Model(solver=MosekSolver())

    size = [21, 11, 15, 9, 34, 25, 41, 52]
    weight = [22, 12, 16, 10, 35, 26, 42, 53]
    capacity = 100

    @variable(model, x[1:length(size)], Bin)
    @objective(model, Max, dot(weight,x))
    @constraint(model, dot(size, x) <= capacity)

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("Solution is:")
    for i = 1:length(size)
        println("x[$i] = ", getvalue(x[i]))
    end
end

main()

github.com

 

下記のようにちゃんとcvxpyの時と同じ結果が得られています。

f:id:meison_amsl:20170731015249p:plain

 

注意点としては、最適化変数をboolenにしたい時は、

Boolを与えるのではなく、JuMP専用のbinを型として与える必要があります。

 

続いて、同じくcvxpyでも実施した

一つの商品を複数回選んで良い場合のコードはこちらです。

using JuMP
using CPLEX

function main()
    println("JuMP knapsack")
    model = Model(solver=CplexSolver())

    size = [21, 11, 15, 9, 34, 25, 41, 52]
    weight = [22, 12, 16, 10, 35, 26, 42, 53]
    capacity = 100

    @variable(model, x[1:length(size)], Int)
    @constraint(model, x .>= 0)
    @objective(model, Max, dot(weight,x))
    @constraint(model, dot(size, x) <= capacity)

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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("Solution is:")
    for i = 1:length(size)
        println("x[$i] = ", getvalue(x[i]))
    end
end

main()

github.com

今回はなんとなくソルバーとしてCPLEXを使ってみました。

ソルバーオブジェクトの名前がCplexSolverであることに注意が必要です。

以前の結果と同じ結果が得られています。

f:id:meison_amsl:20170731021343p:plain

 

ナップサック問題におけるcvxpyとJuMPの計算速度比較

Juliaは計算が早いことで有名な言語なので、

cvxpyとJuMPで計算速度のベンチマークを取ってみました。

下記がcvxpyのコード

#
# A sample code to solve a knapsack_problem with cvxpy
#
# Author: Atsushi Sakai (@Atsushi_twi)
#

import cvxpy
import numpy as np
import random
import time

N_item = 30000

size = np.array([random.randint(1, 50) for i in range(N_item)])
weight = np.array([random.randint(1, 50) for i in range(N_item)])
capacity = 1000

x = cvxpy.Bool(size.shape[0])
objective = cvxpy.Maximize(weight * x)
constraints = [capacity >= size * x]

start = time.time()
prob = cvxpy.Problem(objective, constraints)
prob.solve(solver=cvxpy.ECOS_BB)
elapsed_time = time.time() - start
print("elapsed_time:{0}".format(elapsed_time) + "[sec]")
result = [round(ix[0, 0]) for ix in x.value]

print("status:", prob.status)
print("optimal value", prob.value)

下記がjuliaのコードです。

using JuMP
using CPLEX

function main()
    println("JuMP knapsack")
    model = Model(solver=CplexSolver())

    nitem = 30000
    capacity = 1000

    size = rand(1:50, nitem)
    weight = rand(1:50, nitem)

    @variable(model, x[1:length(size)], Bin)
    @objective(model, Max, dot(weight,x))
    @constraint(model, dot(size, x) <= capacity)

    @time status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
end

main()

github.com

 

ナップサック問題の商品の数を30000個にした場合の比較です。

cvxpy JuMP
計算時間 14.2[s] 1.4[s]

十倍以上早いですね。。。

Juliaのコンパイル時間をいれても3秒ほどだったので

かなり早くなっています。

CPLEXがすごいという噂もありますが。

 

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

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

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

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

 

具体的にJuMPを使って、

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

f:id:meison_amsl:20170802081157p:plain

 

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

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

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

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

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

 

using JuMP
using Ipopt
using NLopt

function main_ipopt()
    println("JuMP sample")
    model = Model(solver=IpoptSolver())

    @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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

function main_nlopt()
    println("JuMP and NL opt sample")
    model = Model(solver=NLoptSolver(algorithm=:LD_MMA))

    @variable(model, z1, start= -0.3)
    @variable(model, z2, start= 0.0)
    @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

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

main_ipopt()
main_nlopt()

github.com

 

参考資料

github.com

Jupyter Notebook Viewer

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

Introduction to Non-Linear Optimization (Computer Science)

Introduction to Non-Linear Optimization (Computer Science)

 

MyEnigma Supporters

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

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

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

myenigma.hatenablog.com