Juliaによる数理最適化入門：簡単な線形計画問題を解いてみる

はじめに

本記事では数理最適化問題をJuliaを用いて解くことが目的です．Juliaの基本的な使い方に関しては後に紹介する記事などを参考にしていただければ幸いです．

本記事の構成は以下のようになります．

• Juliaとは？
• Julia本体と追加で使用するJuliaパッケージのインストール方法
  • モデリング言語パッケージJuMPと線形計画問題ソルバClpの導入
• 簡単な線形計画問題を解いてみる

なお，使用するJuliaバージョンは1.2.0です．

今回私が使用した環境はmacOS Mojave (v10.14)になります．

あと，プログラミングスキルは低めですので，プログラムに関するご指摘はもちろん，用語の使い方などたくさんのコメントをお待ちしております！

Juliaについて

Juliaは2012年にリリースされた比較的新しいプログラミング言語で，LLVMを用いたJITコンパイラによる高速な演算処理を実現します．

他言語と比較してどれほど速いのか具体的に詳しい情報が知りたい方は以下のリンクにある画像に注目していただけるとわかると思います．

• Julia/Benchmarks

まだ数値計算ではバックグラウンドのあるPythonが主流ですが，Juliaは実行処理の速さで近年注目され始めています．

Juliaのインストールとコマンドの設定

Homebrewが利用できる方はコンソールでbrew cask install juliaと入力するとインストールできます．

以下はHomebrewを使わない方法です．一応紹介しておきます．

公式サイトより各OSに対応したパッケージをダウンロードしてインストールしてください（2019年8月現在では最新の安定版はv1.2.0です）．

コマンドの設定（Homebrew使用者はスルー）

インストールした時に生成されるアプリケーション実行ファイルJulia-x.x.x.app（macOSの場合）を起動すればJuliaコンソールが立ち上がりすぐに利用できます．しかし，ターミナルからjuliaと入力することでJuliaコンソールを立ち上げられるようにしておくと，より便利かつ快適に利用できます．

以下のリンク作成のコマンドをターミナルで実行してください（Julia-x.x.appのx.xにはバージョンの数字が入ります．1.2の場合はJulia-1.2.appとなります）．

$ ln -fs "/Applications/Julia-x.x.app/Contents/Resources/julia/bin/julia" /usr/local/bin/julia

以下はJuliaのREPLに入るかどうかの確認となります．

$ which julia
/usr/local/bin/julia
$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.2.0 (2019-08-20)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia>

exit()と入力すればJuliaを終了します．

julia> exit()
$

モデリング言語JuMPのインストール

最適化問題を記述する際に便利となるモデリング言語を導入します．これにより，最適化問題の記述を簡易化します．

先ほどインストールしたJuliaのコンソール（Julia-x.x.app）を起動して以下のコマンドを入力します．先ほどjuliaコマンドの設定をしたので，ターミナルでjuliaと入力すると同じく起動されます．

julia> using Pkg # パッケージ管理の呼び出し
julia> Pkg.add("JuMP") # JuMPの追加（ダウンロード＆インストール）

またはjulia>表示時に]を入力すると，パッケージ管理モードに移行できます．

(v1.2) pkg>

直接gitリンクを指定することで導入することもできます．

(v1.2) pkg> add https://github.com/JuliaOpt/JuMP.jl

JuMPは線形計画問題

\begin{align}
\begin{array}{cl}
\mbox{minimize} & c^\top x\\
\mbox{subject to} & Ax=b,\ x\geq 0
\end{array}
\end{align}

二次計画問題

\begin{align}
\begin{array}{cl}
\mbox{minimize} & \frac{1}{2}x^\top Mx+q^\top x\\
\mbox{subject to} & Ax\leq b
\end{array}
\end{align}

非線形計画問題

\begin{align}
\begin{array}{cl}
\mbox{minimize} & f(x)\\
\mbox{subject to} & g(x)=0,\ h(x)\leq 0
\end{array}
\end{align}

の記述に対応しています．

ロバスト最適化問題や多目的最適化問題などのモデリング言語はそれぞれ別の追加パッケージが必要なので，適宜下記のサイトを参考に探してみてください．

• JuliaOpt

インストールしたパッケージのリストを表示

インストールしたパッケージはJuliaのREPLで以下のように入力することで確認できます．

julia> using Pkg

julia> Pkg.installed()
Dict{String,Union{Nothing, VersionNumber}} with 2 entries:
  "Clp"  => v"0.6.2"
  "JuMP" => v"0.19.2"

本記事では基本的にREPLで使うことはなく，"Hello World"の例のようにスクリプトファイルを作成して，コマンドにより実行させる形を取ることにします．

jlファイルの実行

Juliaのスクリプトの拡張子はjlとなります．試しに"Hello World"をコンソールに出力してみましょう．ファイル名をhello.jlとします．

hello.jl

println("Hello World")

これをjuliaコマンドで実行します．

$ julia hello.jl
Hello World

println文は正常に動作することが確認できました．

さて，『はじめに』でも述べましたが，本記事はJuliaを使って最適化計算を行うことに主眼を置いているため，以降の基本的な処理の文法（if文，forループなど）は別の記事，もしくはJuliaのドキュメントを参考にしてください．

例えば，

• 『Juliaという速くて書きやすい言語をちょっとだけ覗いてみたんだが、なにやらワクワクするものがあったので報告しようと思う』
• Juliaにおける最適化モデリングライブラリJuMP入門

はJulia特有の演算子をコンパクトにわかりやすく説明しています．

試しに線形計画問題を解いてみる

それでは，実際にJuliaを使って以下の簡単な線形計画問題を解いてみましょう．問題のソースはJuliaOpt.org/Quick Start Guideです．

\begin{align}
\begin{array}{cl}
\mbox{Maximize} & 5x+3y\\
\mbox{subject to} & x+5y\leq 3\\
& 0\leq x\leq 2\\
& 0\leq y\leq 30
\end{array}
\end{align}

Clpのインストール

線形計画問題を解くためのソルバーはたくさんありますが，今回はCoin-OR Linear Programming（CLP）を使いたいと思います¹．

Julia REPLで以下のようにコマンドを入力してください．

julia> Pkg.add("Clp")

Julia/JuMPによるコーディング

現在ウェブでみられるJulia用の最適化プログラムはJulia v1.0未満しか対応しておらず，本記事でインストールしたv1.2ではシンタックスエラーとなり動作しません．

残念ながらJuliaOpt.org/Quick Start Guideで公開されているコードも旧バージョン用で動作しないため，新しい表記に直したものを以下のように作成しました．これをtest.jlという名前で保存します．

test.jl

using JuMP
using Clp

# 最適化問題mを定義し，使用するソルバ（ここではOptimizerと呼んでいる）を設定
m = Model(with_optimizer(Clp.Optimizer))

# 区間制約は変数定義と同時に設定可能
@variable(m, 0 <= x <= 2 )
@variable(m, 0 <= y <= 30 )

# 目的関数は最大化，最小化を選べる(乗算演算子*は省略可能)
@objective(m, Max, 5x + 3*y )

# 不等式制約
@constraint(m, 1x + 5y <= 3.0 )

# 最適化問題の表示
print(m)

# 最適化計算の開始
JuMP.optimize!(m)

# 結果の出力（目的関数値，最適解）
println("Objective value: ", JuMP.objective_value(m))
println("x = ", JuMP.value(x))
println("y = ", JuMP.value(y))

ご覧のように，可読性が高く，視認性も高いため，どのような最適化問題を考えているのかコードを読むだけですぐに読み取れるかと思います．

実行結果は以下のように出力されました．

result

Max 5 x + 3 y
Subject to
 x ≥ 0.0
 y ≥ 0.0
 x ≤ 2.0
 y ≤ 30.0
 x + 5 y ≤ 3.0
Coin0506I Presolve 0 (-1) rows, 0 (-2) columns and 0 (-2) elements
Clp3002W Empty problem - 0 rows, 0 columns and 0 elements
Clp0000I Optimal - objective value 10.6
Coin0511I After Postsolve, objective 10.6, infeasibilities - dual 0 (0), primal 0 (0)
Clp0032I Optimal objective 10.6 - 0 iterations time 0.002, Presolve 0.00
Objective value: 10.6
x = 2.0
y = 0.2
dualvar(1x+5y<=3) = -0.6

この結果を読み取ると，最適解$(x^*,y^*)=(2.0,0.2)$が得られ，最適値は$10.6$となりました．

また，制約条件に名前をつけておくことも可能です（以下では制約条件$x+5y\leq 3$に対してconstrと名付けました）．

@constraint(m, constr, 1x + 5y <= 3.0 )

これによって，制約条件constrに対する双対変数（Lagrange乗数）をJuMP.dualを用いることで調べることができます．

println("dualvar(1x+5y<=3) = ", JuMP.dual(constr))

結果：

dualvar(1x+5y<=3) = -0.6

まとめ

本記事ではJuliaのインストール方法と基本操作をわずかですが，紹介しました．また，簡単な線形計画問題を例にモデリング言語JuMPの記法を紹介するとともに，実際にJuliaで最適解を計算しました．

使ってみた感想として，Juliaの公式サイトでは「高速」を一番の売りにしていましたが，どうも本記事の方法ではむしろPythonより遅いのではないかと思うほど，JITコンパイルに時間がかかっていました（実行自体は速いと思いますが，行単位でのJITコンパイルに多くの時間を要している？）．

jlファイルは事前コンパイルを施すことによってC言語のように実行ファイルを作成して実行できるそうですが，如何せん私自身がJulia入門者のため，その方法を見つけ出すことができませんでした．

何か進展があれば，続編を投稿してみたいと思います．

参考にした文献

1. Juliaにおける最適化モデリングライブラリJuMP入門
2. JuMP Julia for Mathematical Optimization/Quick Start Guide

1. 同じくGNU Linear Programming Kit（GLPK）を試しましたが，julia v1.2.0ではパッケージのビルドに失敗するので今回は使わないことにしました． ↩