ラベル Julia の投稿を表示しています。 すべての投稿を表示
ラベル Julia の投稿を表示しています。 すべての投稿を表示

2020年9月3日木曜日

JuMP でビンパッキング問題おためし

Julia の数理最適化パッケージ JuMP を用いて混合整数線型計画問題を試してみる。
JuMP で通常の線形計画問題を解くのは JuMP で線型計画問題おためし を参照。


1. 環境の準備


今回用いるパッケージは下記。
  • JuMP
  • CBC (MILP:混合整数計画問題 のソルバ)
※環境構築に関しては JuMP で線型計画問題おためし を参照
※各 Solver に関しては Getting Solvers を参照


2. ビンパッキング問題


ビンパッキング問題の詳細は Wikipedia を参照してもらうとして、 今回は下記の状況を扱う。

  • それぞれ重さの異なる荷物を大きめのビンに詰める事を考える
  • 使用するビンの数を最も少なくする格納方法を考えたい
  • 荷物 10 個
    • 重さ 2 kg: 4 個
    • 重さ 3 kg: 3 個
    • 重さ 4 kg: 2 個
    • 重さ 5 kg: 1 個
  • ビンごとの最大容量 5 kg
※ちなみに 「ビン」 は 「瓶」 ではなく 「Bin(容器)」 の事


3. 定式化


荷物の数を n とすると必要なビンの数は最大でも n 個となるので、使われないビンの存在も許容してビンの数を(多めに) n と設定する。
このとき各ビンに荷物が 1 つでも格納されたかどうかを判定する変数 y[1〜n] を、各成分に 0 か 1 が格納されるベクトルで表現する。

\begin{eqnarray} \forall j \in \{ 1, 2, \dots, n \}, \ \ y_{j} \in \{ 0, 1 \} \nonumber . \end{eqnarray} 次に行列 $X = \{ x_{i, j} \}$ を下記で定義する。
  • 行(i): 各荷物(1〜n)に対応
  • 列(j): 各ビン(1〜n)に対応
  • 各 $x_{i, j}$ には 0 か 1 が格納される

最後にベクトル $s$ 及びスカラ $B$ を下記で定義する。
  • $s$: 各荷物の重さ
    • 今回の場合は $s = [ 2, 2, 2, 2, 3, 3, 3, 4, 4, 5 ]$
  • $B$: 各ビンごとの最大容量
    • 今回の場合は $B = 5$

3.1 目的関数


出来るだけビンの数を少なくしたいので、各 $y_{j}$ の和を最小とする事を考える。 \begin{eqnarray} minimize \sum_{j} y_{j} \end{eqnarray}

3.2 制約条件(1)


まず、各荷物 1〜n がそれぞれ 1 つのビンに必ず格納される事を表現する。
\begin{eqnarray} \forall i, \ \ \sum_{j} x_{i, j} = 1 \nonumber \end{eqnarray}
この条件は n 個の成分が全て 1 であるベクトル $e (\forall i, \ e_{i} = 1)$ を用いて下記のように書き換え可能。 \begin{eqnarray} X \cdot e = 1 \end{eqnarray}

3.3 制約条件(2)


次に、各ビンに格納された荷物の重さの合計が $B$ を以下となる事を表現する。 \begin{eqnarray} \forall j, \ \ \sum_{i} s_{i} \cdot x_{i, j} \leqq B \cdot y_{j} \nonumber \end{eqnarray}
この条件は行列の転置を用いて下記のように書き換え可能。 \begin{eqnarray} {}^t\!X \cdot s \leqq B \cdot y \end{eqnarray}

3.4 制約条件(3)


最後に、行列 $X$ およびベクトル $y$ の各成分に関する条件を表現する。
\begin{eqnarray} \forall i, j, \ \ x_{i, j} \in \{ 0, 1 \} \\ \forall j, \ \ y_{j} \in \{ 0, 1 \} \end{eqnarray}

3.5 まとめ


上記 (1)〜(5) を用いて下記のようにまとめられる。
\begin{eqnarray} minimize & & \sum_{j} y_{j} \nonumber \\ subject \ to & & X \cdot e = 1 \nonumber \\ & & {}^t\!X \cdot s \leqq B \cdot y \nonumber \\ & & \forall i, j, \ \ x_{i, j} \in \{ 0, 1 \} \nonumber \\ & & \forall j, \ \ y_{j} \in \{ 0, 1 \} \nonumber \end{eqnarray}

4. JuMP


JuMP を用いて最適解を求める。
MILP(混合整数線型問題) を扱う事のできる CBC をソルバーに用いる。
using JuMP
using Cbc
using LinearAlgebra

# 荷物の数 => 最大のビンの数
n = 10

# ビンごとの最大容量
B = 5

# 荷物ごとの重さ
s = [2, 2, 2, 2, 3, 3, 3, 4, 4, 5]


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

# 変数を追加
# X: 荷物とビンの関係
# y: ビン使用の有無
# binary = true で 2 値(0 or 1)変数である事を指定
@variable(model, X[i=1:n, j=1:n], binary = true) # (4)
@variable(model, y[1:n], binary = true)          # (5)

# 目的関数を最小化
@objective(model, Min, sum(y))                   # (1)

# 制約条件を追加
@constraint(model, X * ones(n) .== 1)            # (2)
@constraint(model, transpose(X) * s .<= B .* y)  # (3)

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


### 実行結果 ###

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: May 23 2020 

command line - Cbc_C_Interface -solve -quit (default strategy 1)
Continuous objective value is 6 - 0.00 seconds
〜中略〜

Result - Optimal solution found

Objective value:                7.00000000
Enumerated nodes:               0
Total iterations:               83
Time (CPU seconds):             0.37
Time (Wallclock seconds):       0.12

Total time (CPU seconds):       0.37   (Wallclock seconds):       0.12

4.1 実行結果


正常に終了し、各荷物の各ビンへの格納が正しく行われている事を確認できる。
println("TerminationStatus: ", termination_status(model))
# => TerminationStatus: OPTIMAL

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

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

# 荷物の格納されたビンの情報のみを表示
for i = 1:n, j = 1:n
    x = value(X[i,j])
    if x > 0
        println("X[", i, "," , j, "]: ", x)
    end
end
X[1,2]: 1.0
X[2,6]: 1.0
X[3,7]: 1.0
X[4,9]: 1.0
X[5,7]: 1.0
X[6,9]: 1.0
X[7,2]: 1.0
X[8,10]: 1.0
X[9,8]: 1.0
X[10,4]: 1.0

ビンと荷物の対応
  • ビン 2:
    • 荷物1(2kg)
    • 荷物7(3kg)
  • ビン 4:
    • 荷物10(5kg)
  • ビン 6:
    • 荷物2(2kg)
  • ビン 7:
    • 荷物3(2kg)
    • 荷物5(3kg)
  • ビン 8:
    • 荷物9(4kg)
  • ビン 9:
    • 荷物4(2kg)
    • 荷物6(3kg)
  • ビン 10:
    • 荷物8(4kg)


参考

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 を入力してリターンキーを押下