システムとモデリング

エンジニアリング, DSM, Julia, SysML, Matlab/Simulink,他モデリング全般について。

Juliaで最適化計算

今回はJuliaのパッケージ"JuMP.jl"を使って最適化問題を解いてみます。

JuMP.jlについて

Juliaの最適化計算パッケージです。 www.juliaopt.org

特徴は下記のようですね。

  • 本来の数式に近い記述が可能
  • 速い
  • 記述がsolverに依存しない

例題を実行する。

下記ドキュメントにある例題を実行してみます。

JuMP — Julia for Mathematical Optimization — JuMP -- Julia for Mathematical Optimization 0.18 documentation

using JuMP
using Clp

m = Model(solver = ClpSolver())
@variable(m, 0 <= x <= 2 )
@variable(m, 0 <= y <= 30 )

@objective(m, Max, 5x + 3*y )
@constraint(m, 1x + 5y <= 3.0 )

print(m)

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))

これを実行すると下記が出力されます。

Max x1 + 2 x2 + 5 x3
Subject to
 -x1 + x2 + 3 x3 <= -5
 x1 + 3 x2 - 7 x3 <= 10
 0 <= x1 <= 10
 x2 >= 0
 x3 >= 0
Optimal Solutions:
x1 = 10.0
x2 = 2.1875
x3 = 0.9374999999999999

問題なく最適化計算できたようです。

練習問題を解く

筆者お気に入りのこの教科書の問題を解いてみます。

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

例題3.4

植え付け可能な625エーカーの土地を持つ家族農業経営を考えよう。植え付けを考慮されている作物はトウモロコシ、小麦、オート麦である。1,000エーカーフートの水が灌漑用に必要で、家族は1週あたり300時間仕事ができると見積もられる。その他のデータは下表にまとめられている。最大利益を得るために各作物の植え付け面積を決定しよう。

必要量 (エーカーあたり) トウモロコシ 小麦 オート麦
灌漑用水(エーカーフート) 3.0 1.0 1.5
労働力(人時/週) 0.8 0.2 0.3
収益(ドル) 400 200 250
変数 説明
x1 トウモロコシの植え付け面積 (>=0)
x2 小麦植え付け面積 (>=0)
x3 オート麦植え付け面積 (>=0)

これをJuliaで実装した場合下記のようなコードになります。

using JuMP
using Clp

m = Model(solver=ClpSolver())

@variable(m, x1 >=0)
@variable(m, x2 >=0)
@variable(m, x3 >=0)


@objective(m, Max, 400x1 + 200x2 + 250x3)

@constraint(m, 3.0x1 + 1.0x2 + 1.5x3 <=1000)
@constraint(m,  0.8x1 + 0.2x2 + 0.3x3 <=300)
@constraint(m,  x1 + x2 + x3 <= 625)

print(m)
status = solve(m)

println("Optimal Solutions:")
println("x1 = ", getvalue(x1))
println("x2 = ", getvalue(x2))
println("x3 = ", getvalue(x3)

このとき出力は以下のようになります。

Max 400 x1 + 200 x2 + 250 x3
Subject to
 3 x1 + x2 + 1.5 x3 <= 1000
 0.8 x1 + 0.2 x2 + 0.3 x3 <= 300
 x1 + x2 + x3 <= 625
 x1 >= 0
 x2 >= 0
 x3 >= 0
Optimal Solutions:
x1 = 187.49999999999997
x2 = 437.5
x3 = 0.0

教科書の解答では x1=187.5, x2=437.5, x3=0.0 なので、正解ですね。

今日はここまでにします