From b89f03277bbf958239274860ffa70e0af3801790 Mon Sep 17 00:00:00 2001 From: mtanneau Date: Fri, 4 Mar 2022 20:33:27 -0500 Subject: [PATCH] Fix residual computation --- src/IPM/HSD/HSD.jl | 2 +- src/IPM/MPC/MPC.jl | 2 +- test/IPM/HSD.jl | 75 ++++++++++++++++++++++++++++++++++---------- test/IPM/MPC.jl | 77 ++++++++++++++++++++++++++++++++++++---------- 4 files changed, 121 insertions(+), 35 deletions(-) diff --git a/src/IPM/HSD/HSD.jl b/src/IPM/HSD/HSD.jl index b50f1149..37efe7a4 100644 --- a/src/IPM/HSD/HSD.jl +++ b/src/IPM/HSD/HSD.jl @@ -113,7 +113,7 @@ function compute_residuals!(hsd::HSD{T} # Residuals norm res.rp_nrm = norm(res.rp, Inf) - res.rl_nrm = norm(res.ru, Inf) + res.rl_nrm = norm(res.rl, Inf) res.ru_nrm = norm(res.ru, Inf) res.rd_nrm = norm(res.rd, Inf) res.rg_nrm = norm(res.rg, Inf) diff --git a/src/IPM/MPC/MPC.jl b/src/IPM/MPC/MPC.jl index 76b37782..b2ccefd5 100644 --- a/src/IPM/MPC/MPC.jl +++ b/src/IPM/MPC/MPC.jl @@ -128,7 +128,7 @@ function compute_residuals!(mpc::MPC{T}) where{T} # Residuals norm res.rp_nrm = norm(res.rp, Inf) - res.rl_nrm = norm(res.ru, Inf) + res.rl_nrm = norm(res.rl, Inf) res.ru_nrm = norm(res.ru, Inf) res.rd_nrm = norm(res.rd, Inf) diff --git a/test/IPM/HSD.jl b/test/IPM/HSD.jl index 062ec070..563a1316 100644 --- a/test/IPM/HSD.jl +++ b/test/IPM/HSD.jl @@ -73,26 +73,15 @@ function run_tests_hsd(T::Type) hsd.pt.y .= T.([0, 1]) hsd.pt.zl .= T.([0, 0]) hsd.pt.zu .= T.([0, 0]) - + hsd.pt.τ = 1 hsd.pt.κ = 0 hsd.pt.μ = 0 ϵ = sqrt(eps(T)) - - @testset "Residuals" begin - @inferred TLP.compute_residuals!(hsd) - TLP.compute_residuals!(hsd) - - @test isapprox(hsd.res.rp_nrm, zero(T); atol=ϵ, rtol=ϵ) - @test isapprox(hsd.res.ru_nrm, zero(T); atol=ϵ, rtol=ϵ) - @test isapprox(hsd.res.rd_nrm, zero(T); atol=ϵ, rtol=ϵ) - @test isapprox(hsd.res.rg_nrm, zero(T); atol=ϵ, rtol=ϵ) - - end + TLP.compute_residuals!(hsd) @testset "Convergence" begin - hsd.solver_status = TLP.Trm_Unknown TLP.update_solver_status!(hsd, ϵ, ϵ, ϵ, ϵ) @test hsd.solver_status == TLP.Trm_Optimal @@ -106,8 +95,62 @@ function run_tests_hsd(T::Type) end end +function test_hsd_residuals(T::Type) + # Simple example: + #= + min x1 - x2 + s.τ. x1 + x2 = 1 + x1 - x2 = 0 + 0 <= x1 <= 2 + 0 <= x2 <= 2 + =# + kkt_options = TLP.KKTOptions{T}() + A = Matrix{T}([ + [1 1]; + [1 -1] + ]) + b = Vector{T}([1, 0]) + c = Vector{T}([1, -1]) + c0 = zero(T) + l = Vector{T}([0, 0]) + u = Vector{T}([2, 2]) + dat = Tulip.IPMData(A, b, true, c, c0, l, u) + + hsd = TLP.HSD(dat, kkt_options) + pt = hsd.pt + res = hsd.res + + # Primal-dual solution + x = pt.x .= T[3, 5] + xl = pt.xl .= T[1, 8] + xu = pt.xu .= T[2, 1] + y = pt.y .= T[10, -2] + zl = pt.zl .= T[2, 1] + zu = pt.zu .= T[5, 7] + τ = pt.τ = T(1//2) + κ = pt.κ = T(1//10) + μ = pt.μ = 0 + + TLP.compute_residuals!(hsd) + + @test res.rp ≈ (τ .* b) - A * x + @test res.rl ≈ (τ .* l) - (x - xl) + @test res.ru ≈ (τ .* u) - (x + xu) + @test res.rd ≈ (τ .* c) - A' * y - zl + zu + @test res.rg ≈ c'x - (b'y + l'zl - u'zu) + κ + + @test res.rp_nrm == norm(res.rp, Inf) + @test res.rl_nrm == norm(res.rl, Inf) + @test res.ru_nrm == norm(res.ru, Inf) + @test res.rd_nrm == norm(res.rd, Inf) + @test res.rg_nrm == norm(res.rg, Inf) + + return nothing +end + @testset "HSD" begin - for T in TvTYPES - @testset "$T" begin run_tests_hsd(T) end + @testset "$T" for T in TvTYPES + run_tests_hsd(T) + test_hsd_residuals(T) end -end \ No newline at end of file +end diff --git a/test/IPM/MPC.jl b/test/IPM/MPC.jl index 2818a271..db87b1d4 100644 --- a/test/IPM/MPC.jl +++ b/test/IPM/MPC.jl @@ -73,27 +73,16 @@ function run_tests_mpc(T::Type) ipm.pt.y .= T.([0, 1]) ipm.pt.zl .= T.([0, 0]) ipm.pt.zu .= T.([0, 0]) - + ipm.pt.τ = 1 ipm.pt.κ = 0 ipm.pt.μ = 0 ϵ = sqrt(eps(T)) - - @testset "Residuals" begin - @inferred TLP.compute_residuals!(ipm) - TLP.compute_residuals!(ipm) - - @test isapprox(ipm.res.rp_nrm, zero(T); atol=ϵ, rtol=ϵ) - @test isapprox(ipm.res.ru_nrm, zero(T); atol=ϵ, rtol=ϵ) - @test isapprox(ipm.res.rd_nrm, zero(T); atol=ϵ, rtol=ϵ) - @test isapprox(ipm.res.rg_nrm, zero(T); atol=ϵ, rtol=ϵ) - - end + TLP.compute_residuals!(ipm) @testset "Convergence" begin - - ipm.solver_status = TLP.Trm_Unknown + ipm.solver_status = TLP.Trm_Unknown TLP.update_solver_status!(ipm, ϵ, ϵ, ϵ, ϵ) @test ipm.solver_status == TLP.Trm_Optimal @@ -106,8 +95,62 @@ function run_tests_mpc(T::Type) end end +function test_mpc_residuals(T::Type) + # Simple example: + #= + min x1 - x2 + s.τ. x1 + x2 = 1 + x1 - x2 = 0 + 0 <= x1 <= 2 + 0 <= x2 <= 2 + =# + kkt_options = TLP.KKTOptions{T}() + A = Matrix{T}([ + [1 1]; + [1 -1] + ]) + b = Vector{T}([1, 0]) + c = Vector{T}([1, -1]) + c0 = zero(T) + l = Vector{T}([0, 0]) + u = Vector{T}([2, 2]) + dat = Tulip.IPMData(A, b, true, c, c0, l, u) + + ipm = TLP.MPC(dat, kkt_options) + pt = ipm.pt + res = ipm.res + + # Primal-dual solution + x = pt.x .= T[3, 5] + xl = pt.xl .= T[1, 8] + xu = pt.xu .= T[2, 1] + y = pt.y .= T[10, -2] + zl = pt.zl .= T[2, 1] + zu = pt.zu .= T[5, 7] + pt.τ = 1 + pt.κ = 0 + pt.μ = 0 + + TLP.compute_residuals!(ipm) + + @test res.rp ≈ b - A*x + @test res.rl ≈ l - (x - xl) + @test res.ru ≈ u - (x + xu) + @test res.rd ≈ c - A' * y - zl + zu + @test res.rg == 0 + + @test res.rp_nrm == norm(res.rp, Inf) + @test res.rl_nrm == norm(res.rl, Inf) + @test res.ru_nrm == norm(res.ru, Inf) + @test res.rd_nrm == norm(res.rd, Inf) + @test res.rg_nrm == norm(res.rg, Inf) + + return nothing +end + @testset "MPC" begin - for T in TvTYPES - @testset "$T" begin run_tests_mpc(T) end + @testset "$T" for T in TvTYPES + run_tests_mpc(T) + test_mpc_residuals(T) end -end \ No newline at end of file +end