Add Optimizer for automatic reformulation of nonlinear programs

.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: |

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"

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!
 ```

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
+```

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.

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)

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

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()

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")