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

公開:2020/12/08
更新:2020/12/08
4 min読了の目安(約4100字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では通過費用に関する不等式制約を緩和しています.制約項は eEc(e)xeu\sum_{e\in E} c(e) x_e \leq u の形をしていますので,変数σ\sigmaを追加して目的関数に σ(eEc(x)xeu)\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を定数と思って解くことで,ラグランジュ緩和問題の目的関数値と最適解xx^\starが求まります.

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

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

f(σ)=eEc(e)xeu\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)

実験

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

ϵ\epsilon σ\sigma^\starの値 f(σ)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)を実装してテストしたいと思います.