🎒

[Package] JuMP+CbcとJuMP+SICPの比較: Pisingerのナップサック問題インスタンスを使って

2022/11/07に公開

背景

最適化業界にニュースがありました。

https://twitter.com/timoberthold/status/1588451599906918400

著名なソルバーの1つSCIPが、いつでもどこでも使えるようになりました。

これでライセンス(GurobiやCPLEX)のサーバーにアクセスできないような環境下で最適化のお仕事をしているとき、CBCよりも良い解が得られることが期待できます(アカデミックライセンスが使える方だと関係ないかもしれないですが…)。

使えるようになったことを信じて、早速JuMP+SCIPの実装を試してみます。最初はランダムなナップサック問題を使おうと思ったのですが…(参考)

https://zenn.dev/takilog/articles/c0c52c608be465

こういう簡単に作ったランダムインスタンスは、簡単過ぎてベンチマークに使えないことが知られています。有識者の方が詳しい情報を書いてくれていました。

https://zenn.dev/kounoike/articles/20210327-hard-knapsack

仕方ないので、上の記事に出てくる Pisinger さんのページからダウンロードできる生成器や、インスタンスを使おうと思います。

使ったデータ

Pisinger さんのページからダウンロードできるインスタンス hardinstances_pisinger.tgz を適当に解答し、knapPI_11_10000_1000.csv を見てみます(このインスタンスの作られ方については上の記事とREADME.mdを見て下さい)。

それぞれ適当な区切り(----と空行)でインスタンスと解・目的関数値・計算時間が保存されているので、最良の結果で計算時間が 0.5 以上な場合について、インスタンス情報を読み込み、CbcとSCIPで解きます。ソルバーですが、gap 0.0の設定だと計算が手元のM1 Macで終わらなくてダメだったので、ちょっとだけ設定してあります。Juliaの実装はこんな感じです。

using JuMP
using Cbc
using SCIP
using MathOptInterface

# ソルバー
@enum TargetSolver cbc scip

# ソルバーに応じて設定して解く
function solve(W, V, Capacity; ts=1, tsolver=cbc, level=1, gap=0.01)
    N = length(W)
    
    if tsolver == cbc
        model = Model(Cbc.Optimizer)
        set_optimizer_attribute(model, "log", level)
        set_optimizer_attribute(model, "seconds", ts)
        set_optimizer_attribute(model, "allowableGap", gap)
    elseif tsolver == scip
        model = Model(SCIP.Optimizer)
        set_optimizer_attribute(model, "display/verblevel", level)
        set_optimizer_attribute(model, "limits/time", ts)
        set_optimizer_attribute(model, "limits/gap", gap)
    end
    
    @variable(model, x[1:N], Bin)
    @objective(model, Max, sum(V[i] * x[i] for i in 1:N))
    @constraint(model, sum(W[i] * x[i] for i in 1:N) <= Capacity)
    optimize!(model)
    
    # 設定条件で最適値が得られたとき、解と目的関数値を返す
    # memo: もうちょっとちゃんと終了判定した方が良かった気がする
    if MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        obj = objective_value(model)
        resx = value.(x)
        sol = [i for i in 1:N if resx[i] > 0.5]
        return obj, sol
    else
        return nothing, nothing
    end
end

結果

knapPI_11_10000_1000.csv には、計算時間が 0.5 以上と記録されていたインスタンスが4つあったので、その結果についてCbcとSCIPを「gap 0.01、計算時間 30秒」のパラメータでソルバーを呼び出して解いた結果を表にまとめました。

ソルバーの計算時間は関数呼び出しの前後で普通にtimeで測定しているので、あまり正確じゃないです。

1つ目です(容量1,850,714のインスタンス)。

記録されている結果 Cbc SCIP
目的関数値 2,865,604 - 2,856,216 (99.67%)
計算時間 1.88 27.465 0.827

2つ目です(容量3,685,028のインスタンス)。

記録されている結果 Cbc SCIP
目的関数値 935,272 - 935,268 (99.99%)
計算時間 0.91 18.51 0.997

3つ目です(容量2,480,873のインスタンス)。

記録されている結果 Cbc SCIP
目的関数値 2,488,539 - 2,487,939 (99.97%)
計算時間 0.51 18.38 0.919

4つ目です(容量2,902,736のインスタンス)。

記録されている結果 Cbc SCIP
目的関数値 2,240,382 - 2,239,464 (99.95%)
計算時間 1.03 26.78 0.874

今までフリーでどこでも使えるソルバーで一番使われている(はず)Cbcと比べて、SCIPがかなりいい感じに計算してくれていることが分かります。これでみんな幸せになりますね(?)。

まとめ

我々はSCIPを手に入れた(真面目なベンチマークじゃなくてすいませんでした)。

(追記)ソースコードなどです。

https://github.com/cocomoff/CompareSolversRandomKnapsack

Discussion