Go to Qiita Advent Calendar 2024 Top
LoginSignup
11
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

JuliaAdvent Calendar 2019

Day 15
@takilog(たきたき)

JuMPで数理最適化を楽しむ

Posted at

はじめに

本記事はJulia Advent Calendar 2019の15日目の記事として書きました．コンセプトは「数理計画問題をJulia上のJuMPでサクっと楽しむ」です．書く上で次の資料を参考に書いていたり，例題を借りてきたりしています．著者はこれが終わったら，JuMPの日本語訳をやるつもりです(こういうのってどうやってやりはじめたらいいのでしょう？何も分からない…)．

この記事では以下のような数理最適化の解説資料に沿って，JuMPの使い方を見ていきます．

環境

執筆時点で最新だったJulia1.3の上で実行しました．

  • Julia 1.3.0 (2019-11-26)
  • Cbc v0.6.6
  • JuMP v0.20.1
  • Weave v0.9.1 (執筆用)

線形計画問題のモデリング

栄養問題

栄養問題の例をモデリングしてみます．だいたい式のまま書くだけです．JuMPを使用する場合，変数はvariable，制約はconstraint，目的関数はobjectiveのマクロで記述します．このあたりの具体的な書き方はいろいろ例を見るか，ドキュメント(英語)を見ると分かると思います．

例えば変数ですが

を確認してみると，変数の名前や値域も設定できるようになっているとわかります．ここではx1, x2を定義していますが，arrayでx[i=1:2]という形で持っても大丈夫です．問題を解くときはoptimize!です（モデルの中身に計算がいろいろと記録されるので，!がついてます）．

using JuMP
using Cbc

model = Model(with_optimizer(Cbc.Optimizer))

# model
@variable(model, x₁ ≥ 0)
@variable(model, x₂ ≥ 0)
@objective(model, Min, 50x₁ + 65x₂)
@constraint(model, 3x₁ + 2x₂ ≥ 9)
@constraint(model, x₁/15 + 2x₂/15 ≥ 1/3)
@constraint(model, x₁/6 ≥ 1/4)
@constraint(model, 3x₂ ≥ x₁)
@constraint(model, x₂ ≤ 2x₁)

# solve
optimize!(model)

# get objective value and variables
obj = objective_value(model)
vx1 = value(x₁)
vx2 = value(x₂)
println("Obj $obj (x₁=$vx1, x₂=$vx2)")

輸送問題

次は輸送問題をコーディングしてみます．大体上の栄養問題と同じような感じが，可読性のため変数を配列化しています．

using JuMP
using Cbc

# 例題インスタンスの定数
Nf = 2;
Ns = 3;
C = [[3.4, 2.2, 2.9], [3.4, 2.4, 2.5]];
Supply = [250, 450];
Demand = [200, 200, 200];

# モデルと計算
model = Model(with_optimizer(Cbc.Optimizer, logLevel=0))

@variable(model, x[i=1:Nf, j=1:Ns] ≥ 0)

for factory in 1:Nf
    @constraint(model, sum(x[factory, :]) <= Supply[factory])
end

for store in 1:Ns
    @constraint(model, sum(x[:, store]) == Demand[store])
end

@objective(model, Min, sum(C[i][j] * x[i, j] for i in 1:Nf, j in 1:Ns))
optimize!(model)

for f in 1:Nf, s in 1:Ns
    println("Factory $f -> Store $s: $(value(x[f, s]))")
end
println(objective_value(model))

他にも最適輸送問題（参照: Computational Optimal Transport 精読会 記録 (#1: 2.1-2.3)）におけるKantrovich緩和問題も似た形で解くことができます．全て動くデータ例はgithubのレポジトリに置いてありますので興味がある人は見てみてください．

using Plots
gr()

model = Model(with_optimizer(Cbc.Optimizer, logLevel=0))
N = 100
M = 100
x = rand(100)
y = rand(100)
x ./= sum(x)
y ./= sum(y)

@variable(model, 0.0 <= P[1:N, 1:M] <= 1.0)

# x(row) and y(column) margin
for i in 1:M
    @constraint(model, sum(P[i, j] for j in 1:M) == x[i])
end
for i in 1:N
    @constraint(model, sum(P[j, i] for j in 1:N) == y[i])
end

@objective(model, Min, sum(abs(i - j) * P[i, j] for i in 1:N, j in 1:M));

optimize!(model)

valP = zeros(100, 100)
for i in 1:100, j in 1:100
    valP[i, j] = value(P[i, j])
end

heatmap(1:100, 1:100, valP, colorbar=false, size=(400, 400))

fig1.png

線形回帰問題

次はスライドに出てくる線形回帰問題(係数$a$と切片$b$を求める)を説いてみます．スライドにあるように，誤差(絶対値)を開いて両方向を$z$で抑え，この和を最小化することで切片$a$と$b$を求めます．わかりやすいようにPlots.jlを利用した図を添えました．

まず，疑似データを$a=1.5, b=0.3$として，適当なノイズを乗せて作成します．

using Random
using Plots

# 疑似データの作成
a = 1.5
b = 0.3
n = 20
ϵ = 0.1
x = rand(n)
y = a * x .+ b .+ ϵ * randn(n);
xx = 0:0.01:1;
yy = a * xx .+ b;
scatter(x, y, mark=:circle)
plot!(xx, yy, color=:red)

次に$a$と$b$をデータ$x, y$から求めてみます

using JuMP
using Cbc

model = Model(with_optimizer(Cbc.Optimizer))

@variable(model, a)
@variable(model, b)
@variable(model, z[1:n] >= 0)

for i in 1:n
    @constraint(model, y[i] - (a * x[i] + b) ≥ - z[i])
    @constraint(model, y[i] - (a * x[i] + b) ≤ z[i])
end

@objective(model, Min, sum(z[:]))
optimize!(model)
opt_a = value(a)
opt_b = value(b)
println(opt_a)
println(opt_b)

結果を重ねてプロットしてみると，そこそこ良い回帰直線が得られました．

xxx = 0:0.01:1;
yyy = opt_a * xx .+ opt_b;
scatter(x, y, mark=:circle)
plot!(xx, yy, color=:red)
plot!(xxx, yyy, color=:blue)

fig2.png

整数計画

場合によりますが，例えば2値(0と1)を取るような整数変数は，変数制限時の引数にBinを与えることで実現できます．一般の整数の場合はIntを与えます．ただし整数計画を利用する場合には，整数計画に対応したソルバー(with_optimizerのところ)を用意する必要があります．

小規模な問題であればCbcで大丈夫ですが，大規模な問題を解くためには，おそらく，CPLEX.jlやGurobi.jlを利用する必要があります．私は仕事ではGurobi.jlを通してGurobi 7系もしくは8系を使っています．利用しているGurobiが今後9系になったとき，Gurobi.jlが互換性があるかどうかは分からないです．

JuMPを利用していると，Gurobi.jlのmodelとJuMPの想定しているmodelが別の型になるので注意が必要です．例えばGurobi.jlにはモデルをファイルに出力する機能があります（https://github.com/JuliaOpt/Gurobi.jl/blob/master/src/grb_model.jl#L144-L152 参照）が，JuMPを通している場合はこの関数を利用できません．そのためGurobiが備えている解のファイル出力など（拡張子を.solファイルにする）などを使えません．実装を見てみるとgrb_ccallというマクロ経由で関数を利用しています．実際に実装を見てみると以下のような形式です．

macro grb_ccall(func, args...)
    f = "GRB$(func)"
    args = map(esc,args)
    if Sys.isunix()
        return quote
            ccall(($f,libgurobi), $(args...))
        end
    elseif Sys.iswindows()
        return quote
            ccall(($f,libgurobi), $(esc(:stdcall)), $(args...))
        end
    end
    error("System not recognised.s")
end

function write_model(model::Model, filename::String)
    @assert isascii(filename) # TODO: support non-ascii file names
    ret = @grb_ccall(write, Cint, (Ptr{Cvoid}, Ptr{UInt8}),
        model.ptr_model, filename)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

対応しているGurobiの関数をサイト（https://www.google.co.jp/search?q=GRBwrite&ie=utf-8&oe=utf-8&hl=ja) から確認してみると，いろいろな機能がついています．

filename: The name of the file to be written. The file type is encoded in the file name suffix. Valid suffixes are .mps, .rew, .lp, or .rlp for writing the model itself, .ilp for writing just the IIS associated with an infeasible model (see GRBcomputeIIS for further information), .sol for writing the current solution, .mst for writing a start vector, .hnt for writing a hint file, .bas for writing an LP basis, or .prm for writing modified parameter settings. If your system provides compressing utilities (e.g., 7z or zip for Windows, and gzip, bzip2, or unzip for Linux or Mac OS); then the files can be compressed, so additional suffixes of .gz, .bz2, or .7z are accepted.

そのため，この関数を使いたい場合には，JuMPのモデルから裏側で利用しているbackendのモデルを取ってくるか，MathOptFormat.jlを経由してファイル書き出しをする必要があります．

このあたりの細かい制約についてはまたどこかでまとめたいと思います(フラグ)．

適当にbackendのモデルを取ってくると，せっかくJuMP経由で読みやすい問題を書いた変数名などが飛んでしまうので，MathOptFormat.jlが使える場合は使ったほうがいいのではないかと思います（もしかしたら回避できるかもしれません）．

実際，細かいGurobiのチューニングをしたい場合には，Gurobiが提供しているPython形式の対話式シェルで実行した方が良いことは多いかもしれません．こういう場合には，MathOptFormat.jlを経由してモデルファイルを書き出し，それをGurobi（や他のソルバーのシェル）から呼び出せば良いでしょう．

解説資料

リンク

11
8
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
11
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?