ナップサック問題の分枝限定法-likeな最適化操作をJuMPを使って動作確認する

9 min読了の目安(約8900字TECH技術記事 2

劣勾配法以来の数理最適化ネタになります.

この記事では「しっかり学ぶ数理最適化 (Ume本)」の§4.4あたりにあるナップサック問題の分枝限定法について,動作確認していきたいと思います.なお本記事のプログラムについては以下のバージョンで動作確認しています.

  • Julia 1.5.3
  • JuMP 0.21.5
  • Cbc 0.7.1
  • MathOptInterface 0.9.19

ナップサック問題のインスタンス

書籍の例を最初に試してみます.以下のものです.

\min z(x_1,x_2,x_3,x_4) = 3x_1 + 4x_2 + x_3 + 2x_4
s.t. 2x_1 + 3x_2 + x_3 + 3x_4 \leq 4
all variables are integer

こちらをJuMPで書いてみます.今回は制約不等式が1本ですが,念の為行列っぽい形にしておきます (この記事では特に利用していないですが).

m = Model(with_optimizer(Cbc.Optimizer))
set_optimizer_attribute(m, "logLevel", 0)
c = [3 4 1 2]
A = [2 3 1 3]
b = [4]
@variable(m, x[1:4], Bin)
@objective(m, Max, sum(c[i] * x[i] for i in 1:4))
for i in 1:size(A)[1]
    @constraint(m, sum(A[i, j] * x[j] for j in 1:4) <= b[i])
end

解くと最適解が得られます.ここで得られた最適解z^\star=5.0を分枝限定法で求めていくのが本稿の目的です.

optimize!(m)
println(objective_value(m))
> 5.0
println(getvalue.(x))
> [0.0, 1.0, 1.0, 0.0]

緩和問題 P_0 の実装

上の定式化では,決定変数x_iを二値変数としていましたが,ここでは[0, 1]の値域へ緩和した問題もソルバで解いてみることにします.これは変数宣言をBinから[0, 1]へ変更すれば良いです.

# 変数宣言部
@variable(m, 0 <= x[1:4] <= 1)

緩和問題P_0を解くと,以下の結果を得ることができます.

optimize!(m)
println(objective_value(m))
> 5.666666666666666
println(getvalue.(x))
> [1.0, 0.6666666666666666, 0.0, 0.0]

実際には解がOPTIMALでは求まらない (例えばinfeasible) 場合もあるので,MathOptInterfaceを使って以下のようなブロックで処理することになります.

if termination_status(m) == MathOptInterface.OPTIMAL
    # 目的関数値や変数を取得できる
    obj0 = objective_value(m);
    X0 = value.(x);
    println(obj0, " ", X0)
end

分枝限定法っぽい操作の手作業実装

緩和問題P_0の解が得られたので,解の確認と子問題の生成をやります.

解の確認

得られた解のうち,整数値になっていない要素を探します.Julia的にどのように数値が整数かどかを判定するか難しい (JuMPがFloatにしてしまうので) ですが,ここではceilとfloorをInt的に計算した結果が一致するかどうかで確認することにします.例えばP_0において

X = value.(x)
for i in 1:4
    xif, xic = floor(Int, X[i]), ceil(Int, X[i])
    (xif != xic) && println("x$i is non-integer")
end

と判定します.この出力は

5.666666666666666 [1.0, 0.6666666666666666, 0.0, 0.0]
x2 is non-integer

となって,x_2についてx_2=0またはx_2=1と制限した子問題を作ることになります.

子問題P_1

すでに緩和問題P_0を上で実装しているため,緩和問題は次の制約を付け加えるだけで実装できます.たとえばx_2=0を制限した子問題P_1は以下の一行を追加することで計算でき,出力としては緩和解(x_2=0)として 4.66666 と,更にこの際にx_4は非整数であったことが分かります.

# 問題定義はP0と同じ 
@constraint(m, x[2] == 0)

# 出力部分
4.666666666666666 [1.0, 0.0, 1.0, 0.33333333333333326]
4 is non-integer

子問題P_5

以上の操作を繰り返していくと,最終的にこの問題が得られます.

# 問題
m = Model(with_optimizer(Cbc.Optimizer))
set_optimizer_attribute(m, "logLevel", 0)
c = [3 4 1 2]
A = [2 3 1 3]
b = [4]
@variable(m, 0 <= x[1:4] <= 1)
@objective(m, Max, sum(c[i] * x[i] for i in 1:4))
for i in 1:size(A)[1]
    @constraint(m, sum(A[i, j] * x[j] for j in 1:4) <= b[i])
end
@constraint(m, x[2] == 1)
@constraint(m, x[1] == 0)
optimize!(m)

# 出力
if termination_status(m) == MathOptInterface.OPTIMAL
    obj5 = objective_value(m);
    X5 = value.(x);
    println(obj5, " ", X5)
    for i in 1:4
        xif, xic = floor(Int, X5[i]), ceil(Int, X5[i])
        (xif != xic) && println("$i $xif $xic")
    end
end

# 出力 (Julia)
> 5.0 [0.0, 1.0, 1.0, 0.0]

こうして最適化z^\star=5.0が求められました.

幅優先探索を用いた探索

制約を手作業で加えていく操作をJupyter notebookで試して雰囲気を掴めたので,この操作を自動化します.

解の計算と制限の計算

制約はx_i=jという形なので,これをJuliaの辞書型で持つことにします.実装としては

  • 最適解が緩和問題について求まらないとき,適当な結果を返す (第一引数がfalse)
  • 最適解が緩和問題について求まる場合,以下を返す
    • 最適解が求まったフラグ: true
    • 目的関数値 obj
    • 解 X
    • Xの中で非整数である要素のリスト non_integer_list

という雑な関数programを作ります.得られた結果のnon_integer_listが空であるとき,得られるobjの値は暫定の解として利用することができます.

function problem(; constraints=Dict())
    # 問題
    m = Model(with_optimizer(Cbc.Optimizer))
    set_optimizer_attribute(m, "logLevel", 0)
    
    # 例題 (page 250)
    N = 4
    c = [3 4 1 2]
    A = [2 3 1 3]
    b = [4]
    @variable(m, 0 <= x[1:N] <= 1)
    @objective(m, Max, sum(c[i] * x[i] for i in 1:N))
    for i in 1:size(A)[1]
        @constraint(m, sum(A[i, j] * x[j] for j in 1:N) <= b[i])
    end
    for (k, v) in constraints
        @constraint(m, x[k] == v)
    end
    optimize!(m)
    
    # 実行可能なら解析
    if termination_status(m) == MathOptInterface.OPTIMAL
        obj = objective_value(m)
        X = value.(x)
        non_integer = false
        non_integer_list = Int[]
        for i in 1:N
            xif, xic = floor(Int, X[i]), ceil(Int, X[i])
            if (xif != xic)
                non_integer = true
                push!(non_integer_list, i)
            end
        end
        return true, obj, X, non_integer_list
    else
        return false, nothing, nothing, nothing, nothing
    end
end

作成した情報を探索のためのQueueに乗せる

データ構造を用いて分枝限定法の木探索的なことをやりつつ,最適解を探索していきます.ここで次のように名前を付けておきます.もっと他のデータをもたせてもいいですが,とりあえずはこの形にします.

struct State
    constraints::Dict
end

幅優先探索 (コスト監視付き) を実装する

探索の部分です.最良優先探索をやるべきだった気がしますが,とりあえず幅優先探索をしました.

using DataStructures
queue = Queue{State}()
best_obj = 0.0
best_sol = nothing
root = State(Dict()) # 制限はない
enqueue!(queue, root)

while length(queue) > 0
    state = dequeue!(queue)    
    is_feasible, obj, x, list_non_integer = problem(constraints=state.constraints)
    println(state, " | ", list_non_integer)
    !is_feasible && continue
    
    # 整数解であれば (feasible),暫定最良値を更新
    if length(list_non_integer) == 0
        if best_obj < obj
            best_obj = obj
            best_sol = x
        end
    end

    # もし暫定最良値よりも良い可能性があれば分枝する
    if obj > best_obj
        for n in list_non_integer
            p_zero = copy(state.constraints); p_zero[n] = 0
            p_one = copy(state.constraints); p_one[n] = 1
            enqueue!(queue, State(p_zero))
            enqueue!(queue, State(p_one))
        end
    end
end

# 解を出力する
println()
println("best obj: $(best_obj)")
println("best sol: $(Int.(best_sol))")

適用例

さて,先程の上の例題 (Ume本の例題) について,用意した探索コードを実行した結果を以下に示します.

State(Dict{Any,Any}()) | [2]  # P0です
State(Dict{Any,Any}(2 => 0)) | [4]  # P1です
State(Dict{Any,Any}(2 => 1)) | [1]  # P2です
State(Dict{Any,Any}(4 => 0,2 => 0)) | Int64[]  # P3です
State(Dict{Any,Any}(4 => 1,2 => 0)) | [1]  # P4です
State(Dict{Any,Any}(2 => 1,1 => 0)) | Int64[]  # P5です
State(Dict{Any,Any}(2 => 1,1 => 1)) | nothing  # P6です

best obj: 5.0
best sol: [0, 1, 1, 0]

無事最適解が得られました.

ランダムな例題を試す

上で作成した関数を次の形に直します.

function problem(; constraints=Dict(), c=[3 4 1 2], A=[2 3 1 3], b=[4])
    # 問題
    m = Model(with_optimizer(Cbc.Optimizer))
    set_optimizer_attribute(m, "logLevel", 0)
    
    # 例題 (page 250)
    N = length(c)
    @variable(m, 0 <= x[1:N] <= 1)
    @objective(m, Max, sum(c[i] * x[i] for i in 1:N))
    for i in 1:size(A)[1]
        @constraint(m, sum(A[i, j] * x[j] for j in 1:N) <= b[i])
    end
    for (k, v) in constraints
        @constraint(m, x[k] == v)
    end
    optimize!(m)
    
    # 実行可能なら解析
    if termination_status(m) == MathOptInterface.OPTIMAL
        obj = objective_value(m)
        X = value.(x)
        non_integer_list = Int[]
        for i in 1:N
            xif, xic = floor(Int, X[i]), ceil(Int, X[i])
            if (xif != xic)
                non_integer = true
                push!(non_integer_list, i)
            end
        end
        return true, obj, X, non_integer_list
    else
        return false, nothing, nothing, nothing, nothing
    end
end

また次のように雑なランダム問題を作成する環境を用意します.

M = 10
c = rand(1:10, M);
A = reshape(rand(1:10, M), (1, M));
b = [rand(20:30)];
problem(c=c, A=A, b=b)

例えばN=10の場合に,作成した問題インスタンス(c, A, b)とJuMPで解いた厳密解は次のようになります.

  • c: [8, 9, 6, 5, 3, 10, 7, 9, 8, 7]
  • A: [9 3 3 2 5 1 7 6 5 2]
  • b: [26]
  • obj: 55.0
  • X: [0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0]

これを探索ベースの実装で求められるかどうかを確認してみます.出力はこのようになります.

State(Dict{Any,Any}()) | [7]
State(Dict{Any,Any}(7 => 0)) | [1]
State(Dict{Any,Any}(7 => 1)) | [8]
State(Dict{Any,Any}(7 => 0,1 => 0)) | [5]
State(Dict{Any,Any}(7 => 0,1 => 1)) | [8]
State(Dict{Any,Any}(7 => 1,8 => 0)) | [1]
State(Dict{Any,Any}(7 => 1,8 => 1)) | [9]
State(Dict{Any,Any}(7 => 0,5 => 0,1 => 0)) | Int64[]
State(Dict{Any,Any}(7 => 0,5 => 1,1 => 0)) | [8]
State(Dict{Any,Any}(7 => 0,8 => 0,1 => 1)) | [5]
State(Dict{Any,Any}(7 => 0,8 => 1,1 => 1)) | Int64[]
State(Dict{Any,Any}(7 => 1,8 => 0,1 => 0)) | [5]
State(Dict{Any,Any}(7 => 1,8 => 0,1 => 1)) | [3]
State(Dict{Any,Any}(7 => 1,9 => 0,8 => 1)) | [1]
State(Dict{Any,Any}(7 => 1,9 => 1,8 => 1)) | Int64[]
State(Dict{Any,Any}(7 => 0,5 => 1,8 => 0,1 => 0)) | Int64[]
State(Dict{Any,Any}(7 => 0,5 => 1,8 => 1,1 => 0)) | [9]
State(Dict{Any,Any}(7 => 1,9 => 0,8 => 1,1 => 0)) | [5]
State(Dict{Any,Any}(7 => 1,9 => 0,8 => 1,1 => 1)) | [2]
State(Dict{Any,Any}(7 => 0,9 => 0,5 => 1,8 => 1,1 => 0)) | Int64[]
State(Dict{Any,Any}(7 => 0,9 => 1,5 => 1,8 => 1,1 => 0)) | [3]

best obj: 55.0
best sol: [0, 1, 0, 1, 0, 1, 1, 1, 1, 1]

無事にQueueを利用した探索によって,JuMPによる最適解と同じ結果を求めることができました.