구속된 최단 경로 문제(1) LSP 확인

총결산

제약을 받는 가장 짧은 경로 문제에 대해 열계단법의 계산을 확인

참고 문헌: 궁본 유일랑, 수리 최적화 입문(3): 라그랑의 일완화와 열계단법(교정)https://ci.nii.ac.jp/naid/110009636306

예제

그림 등은 상술한 참고 문헌에 인용하여 아래 그림에서 보여준 가장 짧은 경로 문제를 고려한다.

가장 최적화된 문제 해결(JumP)

우선 위의 최적화 문제를 해결합시다. 이번 환경은 좀 낡았지만 Juria1.1시스템에 설치된 Cbc(v0.6.3)와Jump(v0.119.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

열계단법 실시(1)

개요

위의 참고 문헌에 근거하여 다음과 같은 방법을 실현할 수 있다.

Jump로 라그랑 데이 완화 문제 해결

LSP는 국경 통과 비용의 부등식 제한을 경감시켰다{e\in E} c(e) x_e\lequ의 형태이므로 목적 함수에 변수\sigma\left(\sum{e\in E}) c (x) xe-u\right)을 주고 반대로 제약을 없애고 완화된 문제를 해결한다. 구체적인 실시 방식은 라그랑일 상수\sigma를 받아들이고,완화 문제를 Jump로 다시 푸는 함수를 만듭니다. (예를 들어\sigma에 따라 도표의 무게를 다시 쓰는 최단 노선 문제의 해답기 등을 만들면 됩니다. 여기에는 Jump로 통일적으로 씁니다.)

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

이 함수는\sigma의 형식에 따라 완화 문제를 구체적으로 결정하고\sigma는 상수로 해답을 하며 라그랑일 완화 문제의 목적 함수 값과 최적 해답 x^\star를 구한다.

열계단법 실시(1)

위의 참고 문헌에 의하면 x^\star를 사용하는 열계단은
\partial f(\sigma) =\sum_{e\in E}c(e) x^\star_e - u
의 형식. 이 방법을 사용하여 위의 열계단법(1)을 실시하는 코드는 다음과 같다.\sigma의 초기값은 위의 참고 문헌에 따라 0으로 규정한다.iter 등은 debug에 사용되며, 마지막으로 계산된\sigma에 대해서는 relax를 다시 한 번 반복하십시오최적 값을 계산하기 위해 Solve를 호출합니다.

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)

실험

참고문헌에 따르면\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)을 실시해 테스트하고 싶다.