From e2457a452ebf354b14c88b85198789405e4061d6 Mon Sep 17 00:00:00 2001 From: dfridovi Date: Sat, 16 Nov 2024 19:11:05 -0600 Subject: [PATCH] test runs but getting numerical issues --- Project.toml | 4 ++++ src/MCPSolver.jl | 3 +++ src/mcp.jl | 2 +- src/solver.jl | 37 ++++++++++++++++++++++++------------- test/runtests.jl | 19 +++++++++++-------- 5 files changed, 43 insertions(+), 22 deletions(-) diff --git a/Project.toml b/Project.toml index e5eb7a2..8e604bb 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,8 @@ version = "0.1.0" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be" +Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TrajectoryGamesBase = "ac1ac542-73eb-4349-ae1b-660ab3609574" @@ -15,6 +17,8 @@ TrajectoryGamesBase = "ac1ac542-73eb-4349-ae1b-660ab3609574" BlockArrays = "0.16.43" DataStructures = "0.18.20" FastDifferentiation = "0.4.2" +Infiltrator = "1.8.3" +LinearAlgebra = "1.11.0" SparseArrays = "1.11.0" Symbolics = "6.19.0" TrajectoryGamesBase = "0.3.10" diff --git a/src/MCPSolver.jl b/src/MCPSolver.jl index 740f4f7..93937dd 100644 --- a/src/MCPSolver.jl +++ b/src/MCPSolver.jl @@ -3,6 +3,9 @@ module MCPSolver using SparseArrays: SparseArrays using FastDifferentiation: FastDifferentiation as FD using Symbolics: Symbolics +using LinearAlgebra: I, norm + +using Infiltrator include("SymbolicUtils.jl") include("sparse_utils.jl") diff --git a/src/mcp.jl b/src/mcp.jl index 500842f..af1b8d5 100644 --- a/src/mcp.jl +++ b/src/mcp.jl @@ -51,7 +51,7 @@ function PrimalDualMCP( ) where {T<:Union{FD.Node,Symbolics.Num}} # Create symbolic slack variable `s` and parameter `ϵ`. s_symbolic = SymbolicUtils.make_variables(backend, :s, length(y_symbolic)) - ϵ_symbolic = SymbolicUtils.make_variables(backend, :ϵ, 1) + ϵ_symbolic = only(SymbolicUtils.make_variables(backend, :ϵ, 1)) z_symbolic = [x_symbolic; y_symbolic; s_symbolic] F_symbolic = [ diff --git a/src/solver.jl b/src/solver.jl index a4bfc91..cd01cb7 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -1,6 +1,8 @@ abstract type SolverType end struct InteriorPoint <: SolverType end +using Infiltrator: @infiltrate + """ Basic interior point solver, based on Nocedal & Wright, ch. 19. Computes step directions `δz` by solving the relaxed primal-dual system, i.e. ∇F(z; ϵ) δz = -F(z; ϵ). @@ -19,39 +21,49 @@ function solve( ::InteriorPoint, mcp::PrimalDualMCP; x₀ = zeros(mcp.unconstrained_dimension), - y₀ = zeros(mcp.constrained_dimension), + y₀ = ones(mcp.constrained_dimension), tol = 1e-4 ) x = x₀ y = y₀ s = ones(length(y)) - ϵ = 1.0 + ϵ = 10.0 kkt_error = Inf - while kkt_error > tol + while kkt_error > tol && ϵ > tol iters = 1 while kkt_error > ϵ # Compute the Newton step. - F = mcp.F(x, y, s, ϵ) - δz = -mcp.∇F(x, y, s, ϵ) \ F + F = mcp.F(x, y, s; ϵ) + δz = -(mcp.∇F(x, y, s; ϵ) + ϵ*I) \ F # Fraction to the boundary linesearch. δx = @view δz[1:mcp.unconstrained_dimension] δy = @view δz[(mcp.unconstrained_dimension + 1):(mcp.unconstrained_dimension + mcp.constrained_dimension)] δs = @view δz[(mcp.unconstrained_dimension + mcp.constrained_dimension + 1):end] - α_s = fraction_to_the_boundary_linesearch(s, δs) - α_y = fraction_to_the_boundary_linesearch(y, δy) + α_s = fraction_to_the_boundary_linesearch(s, δs; tol) + α_y = fraction_to_the_boundary_linesearch(y, δy; tol) + + if isnan(α_s) || isnan(α_y) + @warn "Linesearch failed. Exiting prematurely." + break + end # Update variables accordingly. x += α_s * δx s += α_s * δs y += α_y * δy - kkt_error = norm(F) + kkt_error = maximum(abs.(F)) iters += 1 + @info iters, s end + @info x + @info y + @info s + @info ϵ ϵ *= 1 - exp(-iters) end @@ -61,12 +73,11 @@ end """Helper function to compute the step size `α` which solves: α* = max(α ∈ [0, 1] : v + α δ ≥ (1 - τ) v). """ -function fraction_to_the_boundary_linesearch(v, δ; τ = 0.995, decay = 0.5) +function fraction_to_the_boundary_linesearch(v, δ; τ = 0.995, decay = 0.5, tol = 1e-4) α = 1.0 - while any(v + α * δ .≥ (1 - τ) * v) - if α < 1e-2 - @warn "Fraction to the boundary linesarch failed." - break + while any(v + α * δ .< (1 - τ) * v) + if α < tol + return NaN end α *= decay diff --git a/test/runtests.jl b/test/runtests.jl index 9b372a5..50c370a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,3 @@ -using Test: @testset, @test -using MCPSolver - - """ Test for the following QP: min_x 0.5 xᵀ M x s.t. Ax - b ≥ 0. @@ -9,17 +5,24 @@ Taking `y ≥ 0` as a Lagrange multiplier, we obtain the KKT conditions: G(x, y) = Mx - Aᵀy = 0 0 ≤ y ⟂ H(x, y) = Ax - b ≥ 0. """ + +using Test: @testset, @test +using MCPSolver + + @testset "QPTestProblem" begin - M = [2, 1; 1, 2] - A = [1, 0.3] - b = 1.1 + M = [2 1; 1 2] + A = [1.0 0.0; 0.0 1.0] + b = [1.0; 1.0] G(x, y) = M * x - A' * y H(x, y) = A * x - b - mcp = MCPSolver.to_symbolic_mcp(G, H, 2, 1) + mcp = MCPSolver.to_symbolic_mcp(G, H, size(M, 1), length(b)) sol = MCPSolver.solve(MCPSolver.InteriorPoint(), mcp) + println(sol) + @test all(abs.(G(sol.x, sol.y)) .≤ 1e-3) @test all(H(sol.x, sol.y) .≥ 0) @test all(sol.y .≥ 0)