From 79e8a4893a3ab1702f62b01a4457317e47ab5104 Mon Sep 17 00:00:00 2001 From: David Fridovich-Keil Date: Wed, 20 Nov 2024 09:03:19 -0600 Subject: [PATCH] [skip ci] working on custom game solver --- src/game.jl | 187 +++++++++++++++++++++++++++++++++++++++++++++++ src/solver.jl | 2 - test/runtests.jl | 40 ++++++++++ 3 files changed, 227 insertions(+), 2 deletions(-) create mode 100644 src/game.jl diff --git a/src/game.jl b/src/game.jl new file mode 100644 index 0000000..7adf9fa --- /dev/null +++ b/src/game.jl @@ -0,0 +1,187 @@ +"Utility to represent a parameterized optimization problem." +Base.@kwdef struct OptimizationProblem{T1,T2,T3} + objective::T1 + private_equality::T2 = nothing + private_inequality::T3 = nothing +end + +"""A structure to represent a game with objectives and constraints parameterized by +a vector θ ∈ Rⁿ. + +We will assume that the players' primal variables are stacked into a BlockVector +x with the i-th block corresponding to the i-th player's decision variable. +Similarly, we will assume that the parameter θ has a block structure so that the +i-th player's problem depends only upon the i-th block of θ. +""" +struct ParametricGame{T1,T2,T3} + problems::Vector{<:OptimizationProblem} + shared_equality::T1 + shared_inequality::T2 + dims::T3 + + "PrimalDualMCP representation." + mcp::PrimalDualMCP +end + +function ParametricGame(; + test_point, + test_parameter, + problems, + shared_equality = nothing, + shared_inequality = nothing, + parametric_mcp_options = (;), +) + N = length(problems) + @assert N == length(blocks(test_point)) + + dims = dimensions( + test_point, + test_parameter, + problems, + shared_equality, + shared_inequality, + ) + + # Define primal and dual variables for the game, and game parameters.. + # Note that BlockArrays can handle blocks of zero size. + backend = SymbolicUtils.SymbolicsBackend() + 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) + end + hs = map(problems, blocks(x), blocks(θ)) do p, xi, θi + isnothing(p.private_inequality) ? nothing : p.private_inequality(xi, θi) + end + + g̃ = isnothing(shared_equality) ? nothing : shared_equality(x, θ) + h̃ = isnothing(shared_inequality) ? nothing : shared_inequality(x, θ) + + # Build gradient of each player's Lagrangian. + ∇Ls = map(fs, gs, hs, blocks(x), blocks(λ), blocks(μ)) do f, g, h, xi, λi, μi + L = + f - (isnothing(g) ? 0 : sum(λi .* g)) - (isnothing(h) ? 0 : sum(μi .* h)) - (isnothing(g̃) ? 0 : sum(λ̃ .* g̃)) - + (isnothing(h̃) ? 0 : sum(μ̃ .* h̃)) + SymbolicUtils.gradient(L, xi) + end + + # Build MCP representation. + F = Vector{Symbolics.Num}( + filter!( + !isnothing, + [ + reduce(vcat, ∇Ls) + reduce(vcat, gs) + reduce(vcat, hs) + g̃ + h̃ + ], + ), + ) + + z = Vector{Symbolics.Num}( + filter!( + !isnothing, + [ + x + mapreduce(b -> length(b) == 0 ? nothing : b, vcat, blocks(λ)) + mapreduce(b -> length(b) == 0 ? nothing : b, vcat, blocks(μ)) + length(λ̃) == 0 ? nothing : λ̃ + length(μ̃) == 0 ? nothing : μ̃ + ], + ), + ) + + z̲ = [ + fill(-Inf, length(x)) + fill(-Inf, length(λ)) + fill(0, length(μ)) + fill(-Inf, length(λ̃)) + fill(0, length(μ̃)) + ] + + z̅ = [ + fill(Inf, length(x)) + fill(Inf, length(λ)) + fill(Inf, length(μ)) + fill(Inf, length(λ̃)) + fill(Inf, length(μ̃)) + ] + + mcp = PrimalDualMCP(F |> collect, z |> collect, θ |> collect, z̲, z̅) + + ParametricGame(problems, shared_equality, shared_inequality, dims, mcp) +end + +function dimensions( + test_point, + test_parameter, + problems, + shared_equality, + shared_inequality, +) + 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)) + end + μ = map(problems, blocks(test_point), blocks(test_parameter)) do p, xi, θi + isnothing(p.private_inequality) ? 0 : length(p.private_inequality(xi, θi)) + end + + λ̃ = + isnothing(shared_equality) ? 0 : + length(shared_equality(test_point, test_parameter)) + μ̃ = + isnothing(shared_inequality) ? 0 : + length(shared_inequality(test_point, test_parameter)) + + (; x, θ, λ, μ, λ̃, μ̃) +end + +"Solve a parametric game." +function solve( + game::ParametricGame, + θ; + solver_type = InteriorPoint(), + x₀ = nothing, + y₀ = nothing, + tol = 1e-4 +) + initial_x = + !isnothing(x₀) ? x₀ : zeros(sum(game.dims.x) + sum(game.dims.λ) + game.dims.λ̃) + initial_y = !isnothing(y₀) ? y₀ : zeros(sum(game.dims.μ) + game.dims.μ̃) + + (; x, y, s, kkt_error) = solve( + solver_type, + game.mcp; + θ, + x₀ = initial_x, + y₀ = initial_y, + tol + ) + + # Unpack primals per-player for ease of access later. + end_dims = cumsum(game.dims.x) + primals = map(1:num_players(game)) do ii + (ii == 1) ? x[1:end_dims[ii]] : x[(end_dims[ii - 1] + 1):end_dims[ii]] + end + + (; primals, variables = (; x, y), kkt_error) +end + +"Return the number of players in this game." +function num_players(game::ParametricGame) + length(game.problems) +end diff --git a/src/solver.jl b/src/solver.jl index ac6014f..65d650f 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -1,8 +1,6 @@ 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; ϵ). diff --git a/test/runtests.jl b/test/runtests.jl index 08b69de..7e9ba71 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,9 @@ Taking `y ≥ 0` as a Lagrange multiplier, we obtain the KKT conditions: """ using Test: @testset, @test + using MCPSolver +using BlockArrays: BlockArray, Block, mortar, blocks @testset "QPTestProblem" begin M = [2 1; 1 2] @@ -52,3 +54,41 @@ using MCPSolver check_solution(sol) end end + +@testset "ParametricGameTests" begin + lim = 0.5 + + game = ParametricGame(; + test_point = mortar([[1, 1], [1, 1]]), + test_parameter = mortar([[1, 1], [1, 1]]), + problems = [ + OptimizationProblem(; + objective = (x, θi) -> sum((x[Block(1)] - θi) .^ 2), + private_inequality = (xi, θi) -> -abs.(xi) .+ lim, + ), + OptimizationProblem(; + objective = (x, θi) -> sum((x[Block(2)] - θi) .^ 2), + private_inequality = (xi, θi) -> -abs.(xi) .+ lim, + ), + ], + ) + + + θ = mortar([[-1, 0], [1, 1]]) + (; primals, variables, kkt_error) = solve( + game; + θ, + tol = 1e-4, + ) + + @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, + ) +end