Zenn
🛤️

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

2021/02/19に公開2件
Julia
数理最適化
tech

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

この記事では「しっかり学ぶ数理最適化 (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による最適解と同じ結果を求めることができました．

Discussion

こーのいけこーのいけ

分枝限定法の重要な性質の一つに「子問題の緩和問題の最適化は、親問題の緩和問題の最適化結果を使うと効率的に計算可能である」ことがあると思います。
この性質を無視した実装をしてしまうと、分枝限定法は遅いという印象にしかならないでしょう。

整数線形計画問題としてみるならばきちんと単体法の辞書を保存しておき、双対単体法を使うことですね。
ナップサック問題に特化するのであれば、品物を 価値 / 重さ でソートしておき、貪欲法を使って緩和問題を解いていき、0に固定したときには追加で入れるものを探す、1に固定したときは制約を満たせるまで逆順に探す、という工夫により、全体を見ずに子問題の緩和問題の解が得られます。

たきろぐたきろぐ

コメントありがとうございます！それはまったくその通りです．この記事の内容は「分枝限定法をまともに実装する」話ではなく，その動作を確認するためのプログラムを書いたもよです．