Zenn
🎢

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

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

まとめ

  • 制約付き最短路問題に対して，劣勾配法の計算を確認してみます
  • 参考文献:宮本裕一郎, 数理最適化入門(3): ラグランジュ緩和と劣勾配法(チュートリアル) https://ci.nii.ac.jp/naid/110009636306

例題

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

最適化問題を解く(JuMP)

まずは上の最適化問題を解いてみましょう．今回の環境は少し古いですが，Julia 1.1系の上にインストールしたCbc(v0.6.3)とJuMP(v0.19.0)を利用しています．グラフは上の図のグラフ構造を直接埋め込んでいます．

using JuMP
using Cbc

MN = 4
Nodes = [1, 2, 3, 4];  # s, a, b, t
Edges = [(1, 2), (1, 3), (2, 3), (2, 4), (3, 4)];
Ng_plus, Ng_minus = Dict(k => Int[] for k in 1:MN), Dict(k => Int[] for k in 1:MN)
for (u, v) in Edges
    push!(Ng_plus[u], v)
    push!(Ng_minus[v], u)
end
C = Dict((1, 2) => 2,(1, 3) => 6,(2, 3) => 1,(2, 4) => 5,(3, 4) => 3)
U = Dict((1, 2) => 0,(1, 3) => 2,(2, 3) => 9,(2, 4) => 9,(3, 4) => 5)
maxU = 10;

# 最適化問題のhard coding
model = Model(with_optimizer(Cbc.Optimizer))
@variable(model, x[u in Nodes, v in Nodes; (u, v) in Edges], Bin)
@objective(model, Min, sum(x[u, v] * C[u, v] for (u, v) in Edges))
@constraint(model, x[1, 2] + x[1, 3] == 1)
@constraint(model, x[1, 2] == x[2, 3] + x[2, 4])
@constraint(model, x[1, 3] + x[2, 3] == x[3, 4])
@constraint(model, x[2, 4] + x[3, 4] == 1)
@constraint(model, sum(x[u, v] * U[u, v] for (u, v) in Edges) ≤ maxU)

# 解く
optimize!(model)

# 出力する
println(objective_value(model))
for (u, v) in Edges
    println("x($u,$v)=$(Int(value(x[u, v])))")
end

結果は次のようになります．

7.0
x(1,2)=1
x(1,3)=0
x(2,3)=0
x(2,4)=1
x(3,4)=0

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

概要

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

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

LSPでは通過費用に関する不等式制約を緩和しています．制約項は \sum_{e\in E} c(e) x_e \leq u の形をしていますので，変数\sigmaを追加して目的関数に \sigma\left(\sum_{e\in E}c(x)x_e - u\right) を付与し，逆に制約を取り除いて緩和した問題を解きます．具体的な実装としてはラグランジュ定数\sigmaを受け取って，緩和問題をJuMPで解き直すような関数を作ります (ここは例えば，グラフの重みを\sigmaに応じて書き換えた最短路問題のソルバなどで良いですが，ここではまとめてJuMPで書いています)．

function relax_solve(;σ = 1, 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)

    obj = sum(x[u, v] * C[u, v] for (u, v) in Edges)
    obj += σ * (sum(x[u, v] * U[u, v] for (u, v) in Edges) - maxU)  # !!追加!!

    @objective(model, Min, obj)
    
    @constraint(model, sum(x[s, u] for u in Ng_plus[s]) == 1)
    @constraint(model, sum(x[u, t] for u in Ng_minus[t]) == 1)
    
    for n in 1:MN
        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])
        @constraint(model, iterm - oterm == 0)
    end
    
    optimize!(model)
    pv = objective_value(model)
    ∂v = sum(value(x[u, v]) * U[u, v] for (u, v) in Edges) - maxU  # 劣勾配
    return pv, ∂v
end

obj, ∂obj = relax_solve(σ=1)
println(obj)
println(∂obj)

この関数は\sigmaの形に応じて1つの緩和問題が具体的に決まり，\sigmaを定数と思って解くことで，ラグランジュ緩和問題の目的関数値と最適解x^\starが求まります．

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

上の参考文献によれば，x^\starを用いて劣勾配は

\partial f(\sigma) = \sum_{e\in E}c(e) x^\star_e - u

の形です．これを用いて上の劣勾配法(1)を実装したコードが以下のものです．\sigmaの初期値については，上の参考文献に従って0にしています．max_iterとかはdebug用に使ってます．最後に計算で求めた\sigmaに対して，もう一度relax_solveを呼び出して最適値を計算しています．

σ, i, ϵ, max_iter = 0, 1, 0.0001, 2000

while i <= max_iter
    _, ∂fσ = relax_solve(σ=σ)
    (i % 25 == 0) && println("$i\t $σ")

    # step 3
    if abs(∂fσ / i) < ϵ
        println("$i\t $σ")
        break
    end
    
    # step 4
    σ = max(0, σ + ∂fσ / i)
    i += 1
end

opt_val, _ = relax_solve(σ=σ)
println(opt_val)

実験

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

\epsilon \sigma^\starの値 f(\sigma^\star)の値
0.1 0.45992063492063495 6.5400793650793645
0.01 0.21177310907693775 6.788226890923065
0.001 0.2034766797336718 6.796523320266328

計算した結果を参考文献の表と比較してみます．良さそうですね．また最適解が7.0ですので，最後に求まった6.79はそこそこ（？）の数値になっていそうですね（6.79より小さくすることはできません）．

次回はフロー制約を緩和したLKNAPについての劣勾配法(2)を実装してテストしたいと思います．

Discussion