diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 08f30aca4..69cbd76d1 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -17,6 +17,13 @@ jobs: with: # Build documentation on the latest Julia 1.x version: '1' + - name: PolyJuMP + shell: julia --project=docs/ {0} + run: | + using Pkg + Pkg.add([ + PackageSpec(name="PolyJuMP", rev="master"), + ]) - name: Install dependencies shell: julia --project=docs/ {0} run: | diff --git a/Project.toml b/Project.toml index ca1301d3f..7aff00f15 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,7 @@ MultivariateBases = "0.3.4" MultivariateMoments = "0.5" MultivariatePolynomials = "0.5.19" MutableArithmetics = "1" -PolyJuMP = "0.8" +PolyJuMP = "0.8.2" Reexport = "1" SemialgebraicSets = "0.3" StarAlgebras = "0.3" diff --git a/docs/src/reference/internal.md b/docs/src/reference/internal.md index 474ecfe8b..4a90b6193 100644 --- a/docs/src/reference/internal.md +++ b/docs/src/reference/internal.md @@ -5,4 +5,10 @@ SumOfSquares.Certificate.Symmetry.orthogonal_transformation_to SumOfSquares.Certificate.Symmetry._reorder! SumOfSquares.Certificate.Symmetry._rotate_complex PolyJuMP.QCQP._subs_ensure_moi_order +PolyJuMP.Model +PolyJuMP.Solution +PolyJuMP._optimize! +PolyJuMP.recover_solutions +PolyJuMP.nonnegativity_cone +PolyJuMP._invalidate! ``` diff --git a/docs/src/reference/optimizer.md b/docs/src/reference/optimizer.md new file mode 100644 index 000000000..b4b7999e1 --- /dev/null +++ b/docs/src/reference/optimizer.md @@ -0,0 +1,16 @@ +# Optimizer + +Optimizers provide allows solving Polynomial Optimization programs using a hierarchy of relaxations. + +```@docs +PolyJuMP.SAGE.Optimizer +SumOfSquares.Optimizer +PolyJuMP.AbstractRelaxationOptimizer +PolyJuMP.AbstractPolynomialOptimizer +``` + +The degree of the hierarchy is controlled by the following attribute: + +```@docs +PolyJuMP.MultiplierMaxdegree +``` diff --git a/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl b/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl index c24fcf8cd..575b70f65 100644 --- a/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl +++ b/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl @@ -31,7 +31,7 @@ model = Model(Ipopt.Optimizer) @variable(model, a >= 0) @variable(model, b >= 0) @constraint(model, a + b >= 1) -@NLobjective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3) +@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3) optimize!(model) # As we can see below, the termination status is `LOCALLY_SOLVED` and not of `OPTIMAL` @@ -68,11 +68,11 @@ function ∇²f(H, a, b) end using Ipopt gmodel = Model(Ipopt.Optimizer) -@variable(gmodel, a >= 0) -@variable(gmodel, b >= 0) -@constraint(gmodel, a + b >= 1) +@variable(gmodel, α >= 0) +@variable(gmodel, β >= 0) +@constraint(gmodel, α + β >= 1) register(gmodel, :f, 2, f, ∇f, ∇²f) -@NLobjective(gmodel, Min, f(a, b)) +@NLobjective(gmodel, Min, f(α, β)) optimize!(gmodel) # Even if we have the algebraic expressions of gradient and hessian, @@ -85,9 +85,9 @@ solution_summary(gmodel) # and the same solution is found: -@test value(a) ≈ 0.5 rtol=1e-5 #src -@test value(b) ≈ 0.5 rtol=1e-5 #src -value(a), value(b) +@test value(α) ≈ 0.5 rtol=1e-5 #src +@test value(β) ≈ 0.5 rtol=1e-5 #src +value(α), value(β) # ## Sum-of-Squares approach @@ -99,6 +99,56 @@ scs = SCS.Optimizer import Dualization dual_scs = Dualization.dual_optimizer(scs) +# The Sum-of-Squares approach can be applied to the same `model` by simply +# changing its optimizer to `SumOfSquares.Optimizer`. +# This optimizer computes a lower bound using the Sum-of-Squares relaxation +# detailed in the section "How it works" below. +# The SDP solver solving this relaxation is given as argument to its constructor: + +set_optimizer(model, () -> SumOfSquares.Optimizer(dual_scs)) +optimize!(model) + +# The termination status is now `OPTIMAL`: the relaxation was solved to +# optimality so its objective value, queried with `objective_bound`, is a +# **global** lower bound of `0` for the polynomial problem. +# The `result_count` is however zero: no candidate solution could be recovered +# from the relaxation (we detail why at the end of the section "How it works"). +# Combining this lower bound with the solution found by Ipopt, we know at this +# point that the optimal value is in the interval $[0, 1/4]$. + +@test termination_status(model) == MOI.OPTIMAL #src +@test objective_bound(model) ≈ 0 atol = 1e-3 #src +@test result_count(model) == 0 #src +solution_summary(model) + +# ### SAGE approach + +# The Sum-of-Squares certificate is not the only nonnegativity certificate that +# can be used in such relaxation. The SAGE certificate leads to a relative +# entropy program, solved with `PolyJuMP.SAGE.Optimizer` as follows: + +set_optimizer(model, () -> PolyJuMP.SAGE.Optimizer(dual_scs)) +optimize!(model) + +# The SAGE relaxation also certifies (up to the tolerance of the solver) the +# lower bound `0`. This time, a candidate solution is recovered from the dual +# of the relaxation. It is feasible but its objective value $16/27 \approx 0.59$ +# does not close the gap with the lower bound; the interval is still $[0, 1/4]$. + +@test termination_status(model) == MOI.OPTIMAL #src +@test objective_bound(model) ≈ 0 atol = 1e-2 #src +@test result_count(model) == 1 #src +@test primal_status(model) == MOI.FEASIBLE_POINT #src +@test objective_value(model) ≈ 16/27 rtol = 1e-2 #src +solution_summary(model) + +# The recovered candidate is the following: + +@test value(a) ≈ 2/3 rtol = 1e-2 #src +@test value(b) ≈ 2/3 rtol = 1e-2 #src +value(a), value(b) + +# ### How it works # A Sum-of-Squares certificate that $p \ge \alpha$ over the domain `S`, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. # The following program searches for the largest lower bound and finds zero. diff --git a/src/SumOfSquares.jl b/src/SumOfSquares.jl index 2eb5d4900..248259f0d 100644 --- a/src/SumOfSquares.jl +++ b/src/SumOfSquares.jl @@ -65,6 +65,7 @@ Reexport.@reexport using JuMP include("utilities.jl") include("constraints.jl") include("variables.jl") +include("optimizer.jl") function setdefaults!(data::PolyJuMP.Data) PolyJuMP.setdefault!(data, PolyJuMP.NonNegPoly, SOSCone) diff --git a/src/optimizer.jl b/src/optimizer.jl new file mode 100644 index 000000000..868ddf144 --- /dev/null +++ b/src/optimizer.jl @@ -0,0 +1,85 @@ +""" + Optimizer{T}(solver) + +Optimizer computing a bound on the objective value of a polynomial +optimization problem using its Lasserre / Sum-of-Squares relaxation, see +[`PolyJuMP.AbstractRelaxationOptimizer`](@ref) with the cone [`SOSCone`](@ref) +as certificate of nonnegativity. +The level of the hierarchy is chosen for each constraint with the +`PolyJuMP.MultiplierMaxdegree` constraint attribute. +The relaxation is solved with `solver` and the bound can be queried with +`MOI.ObjectiveBound()`. +""" +mutable struct Optimizer{T} <: PolyJuMP.AbstractRelaxationOptimizer{T} + model::PolyJuMP.Model{T} + multiplier_maxdegree::Dict{MOI.ConstraintIndex,Int} + solver::Any + relaxation::Union{Nothing,JuMP.GenericModel{T}} + solutions::Vector{PolyJuMP.Solution{T}} + feasibility_tolerance::T + solve_time::Float64 +end + +function Optimizer{T}( + solver; + feasibility_tolerance = sqrt(Base.rtoldefault(T)), +) where {T} + return Optimizer{T}( + PolyJuMP.Model{T}(), + Dict{MOI.ConstraintIndex,Int}(), + solver, + nothing, + PolyJuMP.Solution{T}[], + feasibility_tolerance, + NaN, + ) +end + +Optimizer(solver; kws...) = Optimizer{Float64}(solver; kws...) + +MOI.get(::Optimizer, ::MOI.SolverName) = "SumOfSquares" + +PolyJuMP.nonnegativity_cone(::Optimizer) = SOSCone() + +""" + PolyJuMP.recover_solutions(model::Optimizer, relaxation, cref, lagrangian) + +Recover candidate solutions from the dual of the SOS constraint `cref` of the +Lagrangian. This dual is a moment matrix and, if the relaxation is tight and +the moments correspond to an atomic measure, the atoms of this measure are +optimal solutions [HL05]. The extraction of the atoms is implemented by +`MultivariateMoments.atomic_measure`; an empty vector of solutions is +returned when it detects that the moment matrix is not atomic. + +[HL05] Henrion, Didier, and Jean-Bernard Lasserre. +"Detecting global optimality and extracting solutions in GloptiPoly." +Positive polynomials in control. Springer (2005): 293-310. +""" +function PolyJuMP.recover_solutions( + model::Optimizer{T}, + relaxation::JuMP.GenericModel{T}, + cref::JuMP.ConstraintRef, + lagrangian, +) where {T} + solutions = PolyJuMP.Solution{T}[] + if JuMP.termination_status(relaxation) != MOI.OPTIMAL + return solutions + end + ν = MultivariateMoments.moment_matrix(cref) + measure = MultivariateMoments.atomic_measure(ν, sqrt(Base.rtoldefault(T))) + if isnothing(measure) + return solutions + end + x = MP.variables(model.model) + for atom in measure.atoms + values = zeros(T, length(x)) + for (j, var) in enumerate(measure.variables) + values[findfirst(isequal(var), x)] = atom.center[j] + end + push!( + solutions, + PolyJuMP.Solution(values, model.model, model.feasibility_tolerance), + ) + end + return solutions +end diff --git a/test/optimizer.jl b/test/optimizer.jl new file mode 100644 index 000000000..c1bb3c3cb --- /dev/null +++ b/test/optimizer.jl @@ -0,0 +1,90 @@ +module TestOptimizer + +using Test + +import MathOptInterface as MOI +using JuMP +using SumOfSquares + +import Clarabel +const SOLVER = + optimizer_with_attributes(Clarabel.Optimizer, MOI.Silent() => true) + +function test_optimizer_attributes() + optimizer = SumOfSquares.Optimizer(SOLVER) + @test optimizer isa SumOfSquares.Optimizer{Float64} + @test MOI.get(optimizer, MOI.SolverName()) == "SumOfSquares" + @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED + @test MOI.get(optimizer, MOI.ResultCount()) == 0 + list = MOI.get(optimizer, MOI.Bridges.ListOfNonstandardBridges{Float64}()) + @test PolyJuMP.Bridges.Constraint.ToPolynomialBridge{Float64} in list + @test PolyJuMP.Bridges.Objective.ToPolynomialBridge{Float64} in list +end + +function test_optimizer_unconstrained() + model = Model(() -> SumOfSquares.Optimizer(SOLVER)) + @variable(model, a) + @objective(model, Min, a^4 - 2a^2 + 1) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test objective_bound(model) ≈ 0 atol = 1e-6 + # The two minimizers `±1` are the atoms of the moment matrix + @test result_count(model) == 2 + @test primal_status(model) == MOI.FEASIBLE_POINT + @test sort([value(a; result = i) for i in 1:2]) ≈ [-1, 1] rtol = 1e-3 + @test objective_value(model) ≈ 0 atol = 1e-6 +end + +function test_optimizer_multiplier_maxdegree() + model = Model(() -> SumOfSquares.Optimizer(SOLVER)) + @variable(model, a) + @objective(model, Min, a) + @constraint(model, con, 1 - a^2 >= 0) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test objective_bound(model) ≈ -1 atol = 1e-6 + # `t` and the constant multiplier of `con` + @test num_variables(unsafe_backend(model).relaxation) == 2 + # The minimizer `-1` is the atom of the moment matrix + @test result_count(model) == 1 + @test primal_status(model) == MOI.FEASIBLE_POINT + @test value(a) ≈ -1 rtol = 1e-3 + MOI.set(model, PolyJuMP.MultiplierMaxdegree(), con, 2) + @test MOI.get(model, PolyJuMP.MultiplierMaxdegree(), con) == 2 + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test objective_bound(model) ≈ -1 atol = 1e-6 + # `t` and the quadratic multiplier of `con` + @test num_variables(unsafe_backend(model).relaxation) == 4 + @test result_count(model) == 1 + @test value(a) ≈ -1 rtol = 1e-3 +end + +function test_optimizer_equality() + model = Model(() -> SumOfSquares.Optimizer(SOLVER)) + @variable(model, a) + @objective(model, Min, a) + @constraint(model, a^2 == 1) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test objective_bound(model) ≈ -1 atol = 1e-6 + # The minimizer `-1` is the atom of the moment matrix + @test result_count(model) == 1 + @test primal_status(model) == MOI.FEASIBLE_POINT + @test value(a) ≈ -1 rtol = 1e-3 + @test objective_value(model) ≈ -1 rtol = 1e-3 +end + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$name" begin + getfield(@__MODULE__, name)() + end + end + end +end + +end # module + +TestOptimizer.runtests() diff --git a/test/runtests.jl b/test/runtests.jl index 9486b290b..0effe0ea2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,7 @@ include("Mock/mock_tests.jl") # Tests needing a solver # FIXME these tests should be converted to Literate and moved to `examples` or # converted to be used with `MockOptimizer` and moved to `test/Tests` +include("optimizer.jl") include("solvers.jl") include("sospoly.jl") include("sosquartic.jl")