はじめに

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

• 整数計画問題の定式化と解法/大阪大学 梅谷先生
• 線形計画問題/筑波大学 久野先生

環境

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

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

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

栄養問題

栄養問題の例をモデリングしてみます．だいたい式のまま書くだけです．

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)")

Obj 197.5 (x₁=2.0, x₂=1.5)

輸送問題

次は こちら(NTTデータ数理システムのページ)に出てくる輸送問題の例 をコーディングしてみます．大体上の栄養問題と同じような感じが，可読性のため変数を配列化しています．

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)

Welcome to the CBC MILP Solver Version: 2.10.3 Build Date: Oct 7 2019

command line - Cbc_C_Interface -solve -quit (default strategy 1) Presolve 4 (-1) rows, 2 (0) columns and 8 (-1) elements 0 Obj 101 Primal inf 2.6833303 (3) 2 Obj 197.5 Optimal - objective value 197.5 After Postsolve, objective 197.5, infeasibilities - dual 0 (0), primal 0 (0 ) Optimal objective 197.5 - 2 iterations time 0.002, Presolve 0.00 Total time (CPU seconds): 0.00 (Wallclock seconds): 0.00

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

Factory 1 -> Store 1: 0.0 Factory 1 -> Store 2: 200.0 Factory 1 -> Store 3: 0.0 Factory 2 -> Store 1: 200.0 Factory 2 -> Store 2: 0.0 Factory 2 -> Store 3: 200.0

println(objective_value(model))

1620.0

最適輸送問題（参照: 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))

 

線形回帰問題

次はスライドに出てくる線形回帰問題(係数\(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)

1.5980308730324215

println(opt_b)

0.23304672266423465

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

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)

 

整数計画

場合によりますが，例えば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

Welcome to the CBC MILP Solver Version: 2.10.3 Build Date: Oct 7 2019

command line - Cbc_C_Interface -solve -quit (default strategy 1) Presolve 40 (0) rows, 22 (0) columns and 120 (0) elements Perturbing problem by 0.001% of 0.97479864 - largest nonzero change 9.55386 93e-06 ( 0.0025402793%) - largest zero change 0 0 Obj 0 Primal inf 29.82845 (20) 21 Obj 1.4828617 Optimal - objective value 1.4828451 Optimal objective 1.482845069 - 21 iterations time 0.002 Total time (CPU seconds): 0.00 (Wallclock seconds): 0.00

write_model (generic function with 1 method)

対応している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（や他のソルバーのシェル）から呼び出せば良いでしょう．

リンク

• この記事を踏まえた例題レポジトリ https://github.com/cocomoff/JuliaAC2019
• JuliaOpt https://github.com/JuliaOpt
• JuMP.jl https://github.com/JuliaOpt/JuMP.jl
• Gurobi.jl https://github.com/JuliaOpt/Gurobi.jl
• MathOptFormat.jl https://github.com/odow/MathOptFormat.jl

---