2020年8月29日土曜日

JuMP で線型計画問題おためし

Julia の数理最適化パッケージ JuMP を用いて線型計画問題のお試しをしてみる。
サンプル問題は ここ の 『1.1 学生宿舎の朝食』 から拝借した。

1. 環境の準備


今回用いる下記パッケージをインストールする。
※プロジェクトの作成および Jupyter を用いる場合は Jupyter でプロジェクトを指定して Julia カーネル追加 を参照 ※各 Solver に関しては Getting Solvers を参照

今回は TestOR プロジェクトを作成済みとして以下を進めていく
# REPL を起動
$ julia

# パッケージ管理モードに移行
julia> ]

# 作成済の TestOR プロジェクトへ移行
(@v1.5) pkg> activate TestOR

# 必要なパッケージを追加
(TestOR) pkg> add JuMP
(TestOR) pkg> add Clp
(TestOR) pkg> add Cbc

2. 線型計画問題として解く


2.1 モデルの作成


変数として下記 2 つを考える。
  • 牛乳 1 単位(1/2 カップ): m
  • シリアル 1 単位(1/4 袋): s
using JuMP
using Clp

# CLP を用いたモデルを作成
model = Model(Clp.Optimizer)

# 変数を追加
# m: 牛乳
# s: シリアル
@variable(model, 0 <= m)
@variable(model, 0 <= s)

# 制約条件を追加
@constraint(model, 9 <= 3 * m + 2 * s)
@constraint(model, 1/3 <= m/15 + (2/15) * s)
@constraint(model, 1/4 <= m/6)
@constraint(model, m/3 <= s)
@constraint(model, s <= 2m)

# 目的関数を最小化
@objective(model, Min, 50 * m + 65 * s)

2.2 モデルの確認


model オブジェクトを確認すると目的関数および制約条件の一覧が表示される(Jupyter の場合)。



2.3 最適化の実施


optimize! コマンドで最適化を実施。
optimize!(model)

最適化の結果を参照するには下記関数の実行が必要。
  • termination_status (終了ステータス)
  • primal_status (Solution Status)
  • objective_value (目的関数の値)
  • value(variable) (変数の値)
※termination_status の結果については Termination statuses を、primal_status の結果については Solution Statuses を参照
println("TerminationStatus: ", termination_status(model))
# => OPTIMAL

println("PrimalStatus: ", primal_status(model))
# => FEASIBLE_POINT

println("ObjectiveValue: ", objective_value(model))
# => 197.5

println("milk: ", value(m))
# => 2.0

println("cerial: ", value(s))
# => 1.5
正しく回答が得られた。

3. 行列を用いて問題を表現する


上記 2. で記載した各制約条件は、行列 $A$ およびベクトル $b, x$ を用いて $b <= Mx $ と表現可能。
行列 $A$ およびベクトル $b$ を適切に表現する事により下記のように表現可能となる。
using JuMP
using Clp
using LinearAlgebra

# 行列 A
A = [
    3 2
    1/15 2/15
    1/6 0
    -1/3 1
    2 -1
]

# ベクトル b
b = [
    9
    1/3
    1/4
    0
    0
]

# CLP を用いたモデルを作成
model = Model(Clp.Optimizer)

# 変数を追加
# x[1]: 牛乳
# x[2]: シリアル
@variable(model, x[1:2] >= 0)

# 制約条件を追加
@constraint(model, b .<= A * x)

# 目的関数を最小化
# dot 関数で内積を計算: 50 * x[1] + 65 * x[2]
@objective(model, Min, LinearAlgebra.dot([50 65], x))

# 最適化の実施
optimize!(model)

結果を表示。
println("TerminationStatus: ", termination_status(model))
# => OPTIMAL

println("PrimalStatus: ", primal_status(model))
# => FEASIBLE_POINT

println("ObjectiveValue: ", objective_value(model))
# => 197.5

println("milk: ", value(x[1]))
# => 2.0

println("cerial: ", value(x[2]))
# => 1.5


4. 混合整数計画問題として解く


上記では牛乳およびシリアルの使用単位に制限が存在しなかったが、下記では牛乳 1/2 カップおよびシリアル 1/4 袋ごとにしか使用できないという前提で問題を考えてみる。
ソルバーとして CBC を指定し、変数を整数に制限する事により整数計画問題(MILP=Mixed-Integer Linear Programming)として問題を解く。
using JuMP
using Cbc # Clp から変更
using LinearAlgebra

# 行列 A
A = [
    3 2
    1/15 2/15
    1/6 0
    -1/3 1
    2 -1
]

# ベクトル b
b = [
    9
    1/3
    1/4
    0
    0
]

# CBC を用いたモデルを作成
# CLP から変更している事に注意
model = Model(Cbc.Optimizer)

# 変数を整数として追加
# x[1]: 牛乳
# x[2]: シリアル
@variable(model, x[i=1:2] >= 0, integer = true) # integer=true で整数を指定("Int" でも可)

# 制約条件を追加
@constraint(model, b .<= A * x)

# 目的関数を最小化
# dot 関数で内積を計算: 50 * x[1] + 65 * x[2]
@objective(model, Min, LinearAlgebra.dot([50 65], x))

# 最適化の実施
optimize!(model)

結果を表示
println("TerminationStatus: ", termination_status(model))
# => OPTIMAL

println("PrimalStatus: ", primal_status(model))
# => FEASIBLE_POINT

println("ObjectiveValue: ", objective_value(model))
# => 215.0

println("milk: ", value(x[1]))
# => 3.0

println("cerial: ", value(x[2]))
# => 1.0
牛乳とシリアルの最適解がそれぞれ整数値として得られている。

参考
  • 講義資料: https://www.cs.tsukuba.ac.jp/~takahito/ucourse/sys_math/part1.pdf

2020年8月27日木曜日

Jupyter でプロジェクトを指定して Julia カーネル追加

以前に試した内容が頭から抜けてたので備忘録として残しておく。

Julia でプロジェクトを作成し、そのプロジェクトを指定して jupyter に kernel を追加する。

下記の手順で進めていく。
  1. Julia で新規 project を作成する
  2. 適当なパッケージを追加する
  3. Jupyter に新規 kernel を追加する
  4. 追加した kernel を試す

1. Project の作成


適当なディレクトリに移動して下記を実施。
# REPL を起動
$ julia

# "]" コマンドでパッケージモードに移行
julia> ]

# "Test" プロジェクトを作成
(@v1.5) pkg> generate Test
 Generating  project Test:
    Test/Project.toml
    Test/src/Test.jl
REPL を抜けてカレントディレクトリを確認すると Test ディレクトリが作成されており、下記のファイルおよびディレクトリが格納されている事が確認できる。
  • Project.toml
  • src

2. パッケージの追加


テスト用に下記パッケージを追加する。
  • IJulia - Jupyter 用パッケージ(必須)
  • Plots
# REPL を起動
$ julia

# "]" コマンドでパッケージモードに移行
julia> ]

# Test プロジェクト環境へ移行
(@v1.5) pkg> activate Test
 Activating environment at `~/prog/julia/Test/Project.toml`

# Test プロジェクト環境で IJulia パッケージ(Jupyter用)をインストール
(Test) pkg> add IJulia

# Plots パッケージ(テスト実行用)をインストール
# GR も同時にインストールされる模様
(Test) pkg> add Plots
下記にてインストール済みパッケージの一覧を確認できる。
(Test) pkg> status
Project Test v0.1.0
Status `~/prog/julia/Test/Project.toml`
  [28b8d3ca] GR v0.51.0
  [7073ff75] IJulia v1.21.3
  [91a5bcdd] Plots v1.6.0

3. Jupyter に Kernel を追加


上記で作成した Test プロジェクト環境を Jupyter 上の Kernel として追加する。
# REPL を起動
$ julia

# "]" コマンドでパッケージモードに移行
julia> ]

# Test プロジェクト環境へ移行
(@v1.5) pkg> activate Test
 Activating environment at `~/prog/julia/Test/Project.toml`

# BackSpace キーでパッケージモードから抜ける
# Test 環境は継続する事に注意
(Test) pkg> (BackSpace キーを押下)

# IJulia.installkernel コマンドで新規 kernel を追加
# "--project" オプションでプロジェクト(Test)を指定している事に注意
julia> using IJulia
julia> IJulia.installkernel("JuliaTest", "--project=/各々の環境でフルパス指定/Test")
[ Info: Installing JuliaTest kernelspec in /各々の環境に依存/Jupyter/kernels/juliatest-1.5
"/各々の環境に依存/Jupyter/kernels/juliatest-1.5"

3.1 設定ファイルでプロジェクトを指定


上記の IJulia.installkernel における "--project" 指定は各 kernel 用の設定ファイルを編集する事でも代用可能となる。
下記にその方法(Mac のみ)を記載する。

Julia のバージョンとプロジェクト名は各々の環境に合わせて変更する。
下記の場合だと "argv" の 5 個目の引数が当該オプションとなるのでこれを編集もしくは追加する事でプロジェクトを指定可能。
$ cat ~/Library/Jupyter/kernels/juliatest-1.5/kernel.json
{
  "display_name": "JuliaTest 1.5.1",
  "argv": [
    "/Applications/Julia-1.5.app/Contents/Resources/julia/bin/julia",
    "-i",
    "--startup-file=yes",
    "--color=yes",
    "--project=/各々の環境でフルパス指定/Test",
    "/各々の環境に依存/.julia/packages/IJulia/tOM8L/src/kernel.jl",
    "{connection_file}"
  ],
  "language": "julia",
  "env": {},
  "interrupt_mode": "signal"
}

3.2 kernel 一覧の確認


下記コマンドで kernel 一覧に今回作成した新規 kernel(juliatest-1.5) が追加されている事が確認できる。
$ jupyter kernelspec list
Available kernels:
  ...その他の kernel...
  juliatest-1.5  /ホームディレクトリのパス/Library/Jupyter/kernels/juliatest-1.5

4. 追加した kernel を試す


Jupyter を起動して kernel の一覧に今回追加した kernel が存在する事を確認。



"JuliaTest 1.5.1" kernel を指定したファイルを開いて下記コードを試す事により、以前に追加したパッケージ Plots を呼び出し可能である事が確認できる。
using Plots
gr()

x = randn(10, 3)
plot(x)


4.1 追加したカーネルの削除


今回のお試し用の環境を削除するにはコンソール上から対象カーネルを指定して下記を実行する。
確認のメッセージが表示されるので "y" を入力してリターンキーを押下するとカーネルの削除が実施される。
$ jupyter kernelspec remove juliatest-1.5
Kernel specs to remove:
  juliatest-1.5  /ホームディレクトリのパス/Library/Jupyter/kernels/juliatest-1.5
Remove 1 kernel specs [y/N]:y # y を入力してリターンキーを押下