Zenn
🎢

制約付き最短路問題で劣勾配法を確認する (2)LKNAP

2020/12/11に公開
Julia
数理最適化
tech

まとめ

例題 (続)

図などは上記の参考文献から引用しています．次の図のような制約付き最短路問題を考えます．

制約は2つあり，1つはネットワークフローが保存される制約(等式制約)で，もう1つは移動に関わる費用コストが与えられた上限uを超えないという制約(不等式制約)でした．前回は不等式制約を緩和して劣勾配を実装し，参考文献に書かれた実験と似た結果が得られることを確認しました．

今回は等式制約を緩和した問題に対する劣勾配法 (劣勾配法2) を実装して確認してみます．参考文献の5章に相当します．

劣勾配法(2)を実装する

概要

上の参考文献によれば，次のようなアプローチを実装すると良いです．


前回は1本の不等式制約(\sum_{e\in E}c(e)x_e \leq u)を緩和していたため，ラグランジュ定数(もしかしてラグランジュ乗数が正しいでしょうか?) \sigmaに関する劣勾配を求めていました．今回はグラフの各ノードが等式制約(ネットワークフロー)を満たすというノード数|V|本文の等式制約があり，ラグランジュ乗数もベクトル\mathbf{\lambda}になります．

劣勾配はベクトル\mathbf{\lambda}に対して得られた解\mathbf{x}^\starから上の式(1式目はフローが出発するソースs，2式目はフローが到着するシンクt，3式目はそれ以外のノードでフローが保存される式であり具体的には頂点aとbに対応します)．

具体的な緩和問題の形は参考文献にある通り

のものですが，これはよく見るとナップサック問題的な最適化問題になっています (そのためL-KNAP => LKNAP，たぶん)．つまり費用制約u分だけ辺を選択することができて，求めたいものは固定されたラグランジュ緩和ベクトル\mathbf{\lambda}に対した解です．ここを高速なソルバで解いて使うと良いですが，ここでは前回と同じく緩和問題自体もJuMPでモデリングした上で解くことにします．

ラグランジュ緩和問題LKNAPをJuMPで解く

前回記事と同じくJulia(w/ JuMP, Cbc Solver)の環境です．前回の実装は雑に隣接行列などをhard codingしていたのですが，少しだけ直しました．後は上の緩和問題をそのままモデリングして解いています．

function relax_solve(;λ = zeros(MN), s = 1, t = 4)
    model = Model(with_optimizer(Cbc.Optimizer, logLevel=0))
    @variable(model, x[u in Nodes, v in Nodes; (u, v) in Edges], Bin)

    # LKNAPの目的関数
    obj = sum(x[u, v] * C[u, v] for (u, v) in Edges)
    obj += λ[s] * (1 - sum(x[s, u] for u in Ng_plus[s]))  # 1式目
    obj += λ[t] * (-1 + sum(x[u, t] for u in Ng_minus[t]))  # 2式目
    for n in 1:MN  # 3式目(s/t以外)
        if n == s || n == t
            continue
        end
        iterm = sum(x[u, n] for u in Ng_minus[n])
        oterm = sum(x[n, u] for u in Ng_plus[n])
        obj += λ[n] * (-oterm + iterm)
    end

    @objective(model, Min, obj)
    @constraint(model, sum(x[u, v] * U[u, v] for (u, v) in Edges) ≤ maxU)
    
    optimize!(model)
    
    # 劣勾配ベクトルの計算
    ∂vg = zeros(MN)
    ∂vg[s] = 1 - sum(value(x[s, u]) for u in Ng_plus[s])  # 1式目
    ∂vg[t] = -1 + sum(value(x[u, t]) for u in Ng_minus[t]) # 2式目
    for n in 1:MN # 3式目(s/t以外)
        if n == s || n == t
            continue
        end
        iterm = sum(value(x[u, n]) for u in Ng_minus[n])
        oterm = sum(value(x[n, u]) for u in Ng_plus[n])
        ∂vg[n] = iterm - oterm
    end

    pv = objective_value(model)
    return pv, ∂vg
end

試しにrelax_solveを呼び出して，最適値と劣勾配を求めてみます．

println(relax_solve(λ=[0.0, 0.0, 0.0, 0.0], s=1, t=4))
# > (0.0, [1.0, 0.0, 0.0, -1.0])

良さそうですね (NaNとかが返ってきていなかったので，とりあえず動作を信じておくことにしましょう)．

劣勾配法(2)を実装する

前回の記事で劣勾配法(1)を実装し，上に緩和問題を解いて劣勾配を計算する関数を用意したので，あとは2つをガッチャンコすると劣勾配法(2)が実装できます．劣勾配ベクトルの初期値は参考文献の通り0ベクトルにしました．

λ = [0.0, 0.0, 0.0, 0.0]
i, ϵ, max_iter = 1, 0.1, 2000

while i <= max_iter
    _, ∂f = relax_solve(λ=λ, s=1, t=4)
    (i % 25 == 0) && println("$i\t $λ")
    
    flag = true
    for v in 1:MN
        if abs(∂f[v] / i) > ϵ
            flag = false
        end
    end
    if flag
        println("end $(λ)")
        break
    end
    
    λ .+= (∂f ./ i)  # ベクトルの更新
    i += 1
end

obj, _ = relax_solve(λ=λ, s=1, t=4)
println(obj)

ほとんど劣勾配法(1)と同じ実装です．

実験

参考文献に従って，\epsilon=0.1,0.01,0.0001に設定した場合の計算結果を表に示します．

\epsilon \mathbf{\lambda}^\star f(\mathbf{\lambda}^\star)
0.1 [2.375, 0.453968, 0.0, -2.8297] 5.203968253968253
0.01 [3.20074, 1.19816, -0.398126, -4.00077] 6.596281283998133
0.001 [3.50634, 1.50048, -0.502052, -4.50477] 7.0

前回の記事の記憶では最適値は7.0でしたので，良いですね．計算した結果を参考文献の表と比較してみます．

\epsilon=0.1のときは少し値が違う結果になっていますが，\epsilonを小さくした結果，s=1,t=4で実装しているため，\lambda_s=3.50634, \lambda_t=-4.50477, \lambda_a=1.50048,\lambda_b=-0.502052となり，ほぼ参考文献に書かれているベクトルの値と最適解に一致しました．良さそうですね．

次回はUme本こと梅谷先生の本で触れられている被覆問題に関する劣勾配法を実装して挙動を確認してみたいと思います (書けたら書くフラグ)．

Discussion