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

1. 環境の準備

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

JuMP
CLP (LP:線形計画問題のソルバ)

CBC (MILP:混合整数計画問題のソルバ)

※各 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月29日土曜日

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