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

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

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

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

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

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

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

# Juliaについて

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

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

- [Julia/Benchmarks](https://julialang.org/benchmarks/)

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

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

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

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

[公式サイト](julialang.org)より各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`となります）．

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

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

```bash
$ 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
julia> exit()
$
```

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

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

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

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

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

```julia
(v1.2) pkg>
```

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

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

JuMPは線形計画問題

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

二次計画問題

```math
\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}
```

非線形計画問題

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

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

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

- [JuliaOpt](https://www.juliaopt.org/packages/)

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

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

```julia
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`とします．

```julia:hello.jl
println("Hello World")
```

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

```bash
$ julia hello.jl
Hello World
```

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

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

例えば，

- 『[Juliaという速くて書きやすい言語をちょっとだけ覗いてみたんだが、なにやらワクワクするものがあったので報告しようと思う](https://qiita.com/sadayuki-matsuno/items/fc5e9ec3894a4b7bfbfb)』
- [Juliaにおける最適化モデリングライブラリJuMP入門](https://myenigma.hatenablog.com/entry/2017/07/31/025341)

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

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

それでは，実際にJuliaを使って以下の簡単な線形計画問題を解いてみましょう．問題のソースは[JuliaOpt.org/Quick Start Guide](http://www.juliaopt.org/JuMP.jl/0.18/quickstart.html#simple-example)です．

```math
\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）を使いたいと思います[^GLPK]．

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

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

```julia
julia> Pkg.add("Clp")
```

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

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

残念ながら[JuliaOpt.org/Quick Start Guide](http://www.juliaopt.org/JuMP.jl/0.18/quickstart.html#simple-example)で公開されているコードも旧バージョン用で動作しないため，新しい表記に直したものを以下のように作成しました．これを`test.jl`という名前で保存します．

```julia: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))
```

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

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

```bash: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`と名付けました）．

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

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

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

結果：

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

# まとめ

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

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

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

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

参考にした文献
---
1. [Juliaにおける最適化モデリングライブラリJuMP入門](https://myenigma.hatenablog.com/entry/2017/07/31/025341)
2. [JuMP Julia for Mathematical Optimization/Quick Start Guide](http://www.juliaopt.org/JuMP.jl/0.18/quickstart.html)