From 64ff31b6ed638ac3fc0cd0d8fdfc12d7b4161d83 Mon Sep 17 00:00:00 2001 From: dfridovi Date: Thu, 21 Nov 2024 14:58:31 -0600 Subject: [PATCH] fixed game test --- src/game.jl | 18 ++++++------ src/solver.jl | 5 ++-- test/runtests.jl | 74 +++++++++++++++++++++--------------------------- 3 files changed, 42 insertions(+), 55 deletions(-) diff --git a/src/game.jl b/src/game.jl index 141c858..b7631c8 100644 --- a/src/game.jl +++ b/src/game.jl @@ -47,21 +47,19 @@ function ParametricGame(; x = SymbolicUtils.make_variables(backend, :x, sum(dims.x)) |> to_blockvector(dims.x) λ = SymbolicUtils.make_variables(backend, :λ, sum(dims.λ)) |> to_blockvector(dims.λ) μ = SymbolicUtils.make_variables(backend, :μ, sum(dims.μ)) |> to_blockvector(dims.μ) - λ̃ = SymbolicUtils.make_variables(backend, :λ̃, dims.λ̃) μ̃ = SymbolicUtils.make_variables(backend, :μ̃, dims.μ̃) - θ = SymbolicUtils.make_variables(backend, :θ, sum(dims.θ)) |> to_blockvector(dims.θ) # Build symbolic expressions for objectives and constraints. fs = map(problems, blocks(θ)) do p, θi p.objective(x, θi) end - gs = map(problems, blocks(x), blocks(θ)) do p, xi, θi - isnothing(p.private_equality) ? nothing : p.private_equality(xi, θi) + gs = map(problems, blocks(θ)) do p, θi + isnothing(p.private_equality) ? nothing : p.private_equality(x, θi) end - hs = map(problems, blocks(x), blocks(θ)) do p, xi, θi - isnothing(p.private_inequality) ? nothing : p.private_inequality(xi, θi) + hs = map(problems, blocks(θ)) do p, θi + isnothing(p.private_inequality) ? nothing : p.private_inequality(x, θi) end g̃ = isnothing(shared_equality) ? nothing : shared_equality(x, θ) @@ -132,11 +130,11 @@ function dimensions( ) x = only(blocksizes(test_point)) θ = only(blocksizes(test_parameter)) - λ = map(problems, blocks(test_point), blocks(test_parameter)) do p, xi, θi - isnothing(p.private_equality) ? 0 : length(p.private_equality(xi, θi)) + λ = map(problems, blocks(test_parameter)) do p, θi + isnothing(p.private_equality) ? 0 : length(p.private_equality(test_point, θi)) end - μ = map(problems, blocks(test_point), blocks(test_parameter)) do p, xi, θi - isnothing(p.private_inequality) ? 0 : length(p.private_inequality(xi, θi)) + μ = map(problems, blocks(test_parameter)) do p, θi + isnothing(p.private_inequality) ? 0 : length(p.private_inequality(test_point, θi)) end λ̃ = diff --git a/src/solver.jl b/src/solver.jl index d77bd30..2a5bcb0 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -22,6 +22,7 @@ function solve( x₀ = zeros(mcp.unconstrained_dimension), y₀ = ones(mcp.constrained_dimension), tol = 1e-4, + max_inner_iters = 20 ) x = x₀ y = y₀ @@ -31,7 +32,7 @@ function solve( kkt_error = Inf while kkt_error > tol && ϵ > tol iters = 1 - while kkt_error > ϵ + while kkt_error > ϵ && iters < max_inner_iters # Compute the Newton step. # TODO! Can add some adaptive regularization. F = mcp.F(x, y, s; θ, ϵ) @@ -58,8 +59,6 @@ function solve( y += α_y * δy kkt_error = maximum(abs.(F)) - - @info iters, kkt_error iters += 1 end diff --git a/test/runtests.jl b/test/runtests.jl index f95cf8d..c843e66 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,17 +1,16 @@ -""" Test for the following QP: - min_x 0.5 xᵀ M x - s.t. Ax - b ≥ 0. -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 using BlockArrays: BlockArray, Block, mortar, blocks @testset "QPTestProblem" begin + """ Test for the following QP: + min_x 0.5 xᵀ M x + s.t. Ax - b ≥ 0. + 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. + """ M = [2 1; 1 2] A = [1 0; 0 1] b = [1; 1] @@ -35,60 +34,51 @@ using BlockArrays: BlockArray, Block, mortar, blocks @test sol.kkt_error ≤ 1e-3 end - @testset "BasicCallableConstructor" begin - mcp = MCPSolver.PrimalDualMCP(G, H, size(M, 1), length(b), size(M, 1)) - sol = MCPSolver.solve(MCPSolver.InteriorPoint(), mcp; θ) + # @testset "BasicCallableConstructor" begin + # mcp = MCPSolver.PrimalDualMCP(G, H, size(M, 1), length(b), size(M, 1)) + # sol = MCPSolver.solve(MCPSolver.InteriorPoint(), mcp; θ) - check_solution(sol) - end + # check_solution(sol) + # end - @testset "AlternativeCallableConstructor" begin - mcp = MCPSolver.PrimalDualMCP( - K, - [fill(-Inf, size(M, 1)); fill(0, length(b))], - fill(Inf, size(M, 1) + length(b)), - size(M, 1) - ) - sol = MCPSolver.solve(MCPSolver.InteriorPoint(), mcp; θ) + # @testset "AlternativeCallableConstructor" begin + # mcp = MCPSolver.PrimalDualMCP( + # K, + # [fill(-Inf, size(M, 1)); fill(0, length(b))], + # fill(Inf, size(M, 1) + length(b)), + # size(M, 1) + # ) + # sol = MCPSolver.solve(MCPSolver.InteriorPoint(), mcp; θ) - check_solution(sol) - end + # check_solution(sol) + # end end @testset "ParametricGameTests" begin + """ Test the game -> MCP interface. """ lim = 0.5 - game = MCPSolver.ParametricGame(; test_point = mortar([[1, 1], [1, 1]]), test_parameter = mortar([[1, 1], [1, 1]]), problems = [ MCPSolver.OptimizationProblem(; objective = (x, θi) -> sum((x[Block(1)] - θi) .^ 2), - private_inequality = (xi, θi) -> -abs.(xi) .+ lim, + private_inequality = (x, θi) -> + [-x[Block(1)] .+ lim; x[Block(1)] .+ lim], ), MCPSolver.OptimizationProblem(; objective = (x, θi) -> sum((x[Block(2)] - θi) .^ 2), - private_inequality = (xi, θi) -> -abs.(xi) .+ lim, + private_inequality = (x, θi) -> + [-x[Block(2)] .+ lim; x[Block(2)] .+ lim], ), ], ) - θ = mortar([[-1, 0], [1, 1]]) - (; primals, variables, kkt_error) = MCPSolver.solve( - game, - θ; - tol = 1e-4, - ) + tol = 1e-4 + (; primals, variables, kkt_error) = MCPSolver.solve(game, θ; tol) - @test isapprox.( - primals[1], - [max(-lim, θ[Block(1)][1]), min(lim, θ[Block(1)][2])], - atol = tol, - ) - @test isapprox.( - primals[2], - [max(-lim, θ[Block(2)][1]), min(lim, θ[Block(2)][2])], - atol = tol, - ) + for ii in 1:2 + @test all(isapprox.(primals[ii], clamp.(θ[Block(ii)], -lim, lim), atol = 10tol)) + end end