diff --git a/README.md b/README.md
index 52a725b..b86c53c 100644
--- a/README.md
+++ b/README.md
@@ -63,7 +63,7 @@ This allows mixing SAGE and SOS constraints in the same model.
 
 ### Polynomial optimization
 
-PolyJuMP also allows solving polynomial optimization problems using the `QCQP` and `KKT` solvers.
+PolyJuMP also allows solving polynomial optimization problems using the `QCQP`, `KKT` and `SAGE` solvers.
 Polynomial optimization problems do not involve any symbolic variables from DynamicPolynomials or TypedPolynomials,
 instead all variables are JuMP decision variables.
 
@@ -89,6 +89,53 @@ model = Model(optimizer_with_attributes(
 ))
 ```
 
+The `SAGE` solver computes a bound on the optimal objective value using the
+SAGE relaxation of the problem [CP16, MCW21]: it certifies the nonnegativity
+of the Lagrangian with the SAGE cone, using one SAGE multiplier for each
+inequality constraint and one free polynomial multiplier for each equality
+constraint.
+It is parametrized by an inner solver for the resulting relative entropy
+program. For instance, to compute a lower bound on the minimum of the Motzkin
+polynomial with `ECOS.Optimizer` as inner solver, use:
+```julia
+using JuMP, PolyJuMP, ECOS
+model = Model(() -> PolyJuMP.SAGE.Optimizer(ECOS.Optimizer))
+@variable(model, x)
+@variable(model, y)
+@objective(model, Min, x^4 * y^2 + x^2 * y^4 + 1 - 3 * x^2 * y^2)
+optimize!(model)
+objective_bound(model) # ≈ 0
+```
+The bound is returned as `MOI.ObjectiveBound`. In addition, candidate
+solutions are recovered from the dual of the SAGE constraint, which is a
+vector of pseudo-moments, following [MCW21, Section 4.2] (see also its
+reference implementation `poly_solrec` in
+[sageopt](https://github.com/rileyjmurray/sageopt)); `result_count(model)`
+gives the number of candidates found, sorted by feasibility and objective
+value. In the example above, the four minimizers `(±1, ±1)` are recovered:
+```julia
+value(x; result = 1), value(y; result = 1) # ≈ (1, 1)
+```
+The maximum degree of the multiplier of a constraint is chosen with the
+`PolyJuMP.MultiplierMaxdegree` constraint attribute:
+```julia
+@constraint(model, con, x^2 >= 1)
+MOI.set(model, PolyJuMP.MultiplierMaxdegree(), con, 2)
+```
+The `SumOfSquares.Optimizer` of [SumOfSquares.jl](https://github.com/jump-dev/SumOfSquares.jl)
+is the analogous solver certifying the nonnegativity of the Lagrangian with
+the SOS cone instead; increasing the `PolyJuMP.MultiplierMaxdegree` attributes
+then gives the higher levels of the Lasserre hierarchy.
+
+[CP16] Chandrasekaran, Venkat, and Parikshit Shah.
+*Relative entropy relaxations for signomial optimization.*
+SIAM Journal on Optimization 26.2 (2016): 1147-1173.
+
+[MCW21] Murray, Riley, Venkat Chandrasekaran, and Adam Wierman.
+*Signomials and polynomial optimization via relative entropy and partial
+dualization.* Mathematical Programming Computation 13 (2021): 257-295.
+[arXiv:1907.00814](https://arxiv.org/abs/1907.00814)
+
 ## Documentation
 
 Documentation for `PolyJuMP.jl` is included in the
diff --git a/src/KKT/KKT.jl b/src/KKT/KKT.jl
index bb36d6e..c973b24 100644
--- a/src/KKT/KKT.jl
+++ b/src/KKT/KKT.jl
@@ -13,7 +13,7 @@ Base.@kwdef mutable struct Options{T}
     feasibility_tolerance::T = Base.rtoldefault(T)
 end
 
-mutable struct Optimizer{T} <: MOI.AbstractOptimizer
+mutable struct Optimizer{T} <: PolyJuMP.AbstractPolynomialOptimizer{T}
     model::PolyJuMP.Model{T}
     options::Options{T}
     # Result
@@ -40,34 +40,6 @@ Optimizer() = Optimizer{Float64}()
 
 MOI.get(::Optimizer, ::MOI.SolverName) = "PolyJuMP.KKT"
 
-MOI.is_empty(model::Optimizer) = MOI.is_empty(model.model)
-
-function MOI.empty!(model::Optimizer)
-    MOI.empty!(model.model)
-    invalidate_solutions!(model)
-    return
-end
-
-function MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute)
-    return MOI.supports(model.model, attr)
-end
-
-function MOI.set(model::Optimizer, attr::MOI.AbstractModelAttribute, value)
-    MOI.set(model.model, attr, value)
-    invalidate_solutions!(model)
-    return
-end
-
-function MOI.get(
-    model::Optimizer,
-    attr::Union{
-        MOI.AbstractModelAttribute,
-        MOI.Bridges.ListOfNonstandardBridges,
-    },
-)
-    return MOI.get(model.model, attr)
-end
-
 function MOI.supports(::Optimizer{T}, attr::MOI.RawOptimizerAttribute) where {T}
     return hasfield(Options{T}, Symbol(attr.name))
 end
@@ -87,7 +59,7 @@ function MOI.get(model::Optimizer, attr::MOI.RawOptimizerAttribute)
     return getfield(model.options, Symbol(attr.name))
 end
 
-function invalidate_solutions!(model::Optimizer)
+function PolyJuMP._invalidate!(model::Optimizer)
     empty!(model.solutions)
     model.solve_time = NaN
     model.termination_status = MOI.OPTIMIZE_NOT_CALLED
@@ -95,36 +67,6 @@ function invalidate_solutions!(model::Optimizer)
     return
 end
 
-MOI.is_valid(model::Optimizer, i::MOI.Index) = MOI.is_valid(model.model, i)
-
-function MOI.add_variable(model::Optimizer)
-    invalidate_solutions!(model)
-    return MOI.add_variable(model.model)
-end
-
-function MOI.supports_constraint(
-    model::Optimizer,
-    ::Type{F},
-    ::Type{S},
-) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
-    return MOI.supports_constraint(model.model, F, S)
-end
-
-function MOI.add_constraint(
-    model::Optimizer,
-    func::MOI.AbstractFunction,
-    set::MOI.AbstractSet,
-)
-    ci = MOI.add_constraint(model.model, func, set)
-    invalidate_solutions!(model)
-    return ci
-end
-
-MOI.supports_incremental_interface(::Optimizer) = true
-function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
-    return MOI.Utilities.default_copy_to(dest, src)
-end
-
 function _add_to_system(system, lagrangian, ::SS.FullSpace, ::Bool)
     return lagrangian
 end
@@ -176,7 +118,7 @@ function _square(x::Vector{T}, n) where {T}
     return T[(i + n in eachindex(x)) ? x[i] : x[i]^2 for i in eachindex(x)]
 end
 
-function _optimize!(model::Optimizer{T}) where {T}
+function PolyJuMP._optimize!(model::Optimizer{T}) where {T}
     if isnothing(model.options.solver)
         system = SS.AlgebraicSet{T,PolyJuMP.PolyType{T}}()
     else
@@ -242,12 +184,6 @@ function _optimize!(model::Optimizer{T}) where {T}
     return
 end
 
-function MOI.optimize!(model::Optimizer)
-    return model.solve_time = @elapsed _optimize!(model)
-end
-
-MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time
-
 function MOI.get(model::Optimizer, ::MOI.RawStatusString)
     if model.termination_status === MOI.OPTIMIZE_NOT_CALLED
         return "`optimize!` has not yet been called"
diff --git a/src/PolyJuMP.jl b/src/PolyJuMP.jl
index 54ed872..054a668 100644
--- a/src/PolyJuMP.jl
+++ b/src/PolyJuMP.jl
@@ -33,6 +33,7 @@ include("data.jl")
 include("default.jl")
 
 include("model.jl")
+include("optimizer.jl")
 include("KKT/KKT.jl")
 include("QCQP/QCQP.jl")
 include("SAGE/SAGE.jl")
diff --git a/src/SAGE/SAGE.jl b/src/SAGE/SAGE.jl
index c4aaed7..00303e2 100644
--- a/src/SAGE/SAGE.jl
+++ b/src/SAGE/SAGE.jl
@@ -182,4 +182,6 @@ function PolyJuMP.bridges(
     )]
 end
 
+include("optimizer.jl")
+
 end
diff --git a/src/SAGE/bridges/age.jl b/src/SAGE/bridges/age.jl
index ce57448..edf4803 100644
--- a/src/SAGE/bridges/age.jl
+++ b/src/SAGE/bridges/age.jl
@@ -103,6 +103,32 @@ function MOI.Bridges.added_constraint_types(
     return [(F, MOI.EqualTo{T}), (G, MOI.RelativeEntropyCone)]
 end
 
+# The coefficients `c` of the AGE constraint only appear in the relative
+# entropy constraint `(c_k + ∑ν, c_{-k}, ν) ∈ RelativeEntropyCone` where `k`
+# is the index of the distinguished monomial. The dual is the adjoint of this
+# map applied to the dual `(u, v, w)` of the relative entropy constraint,
+# that is, `u` for the entry `k` and `v` for the other entries. This is the
+# dual AGE cone characterization given by the conic duality of the relative
+# entropy formulation [MCW21, (2)].
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDual,
+    bridge::AGEBridge,
+)
+    dual = MOI.get(model, attr, bridge.relative_entropy_constraint)
+    m = div(length(dual) + 1, 2)
+    v = Vector{eltype(dual)}(undef, m)
+    v[bridge.k] = dual[1]
+    j = 1
+    for i in 1:m
+        if i != bridge.k
+            j += 1
+            v[i] = dual[j]
+        end
+    end
+    return v
+end
+
 function MOI.Bridges.Constraint.concrete_bridge_type(
     ::Type{<:AGEBridge{T}},
     H::Type{<:MOI.AbstractVectorFunction},
diff --git a/src/SAGE/bridges/sage.jl b/src/SAGE/bridges/sage.jl
index dcadad1..6b5d42d 100644
--- a/src/SAGE/bridges/sage.jl
+++ b/src/SAGE/bridges/sage.jl
@@ -69,6 +69,21 @@ function MOI.Bridges.Constraint.concrete_bridge_type(
     return SAGEBridge{T,F,G}
 end
 
+# The signomial SAGE constraint `func ∈ SAGE` is reformulated into the
+# equality constraints `∑_k ν[k, i] - func_i = 0` in which `func` appears
+# with coefficient `-1` so the dual is `-μ` where `μ` is the dual of these
+# equality constraints; this is the adjoint of the reformulation map.
+# Since the SAGE cone is the sum of the AGE cones, its dual is the
+# intersection of the duals of the AGE cones and indeed, at the optimum, the
+# dual of each constraint `ν[k, :] ∈ AGE` also equals `-μ`.
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDual,
+    bridge::SAGEBridge,
+)
+    return [-MOI.get(model, attr, ci) for ci in bridge.equality_constraints]
+end
+
 function MOI.get(
     model::MOI.ModelLike,
     attr::DecompositionAttribute,
diff --git a/src/SAGE/bridges/signomial.jl b/src/SAGE/bridges/signomial.jl
index 90c47ba..57806c1 100644
--- a/src/SAGE/bridges/signomial.jl
+++ b/src/SAGE/bridges/signomial.jl
@@ -11,30 +11,46 @@ https://arxiv.org/abs/1810.01614
 Mathematical Programming Computation 13 (2021): 257-295.
 https://arxiv.org/pdf/1907.00814.pdf
 """
-struct SignomialsBridge{T,S,P,F} <: MOI.Bridges.Constraint.AbstractBridge
+struct SignomialsBridge{T,S,P,F,G} <: MOI.Bridges.Constraint.AbstractBridge
+    # Indices of the rows `i` of `set.α` that have an odd entry
+    odd::Vector{Int}
+    # For each odd row `i`, the constraint `vi - g[i] ≤ 0`
+    lower::Vector{MOI.ConstraintIndex{G,MOI.LessThan{T}}}
+    # For each odd row `i`, the constraint `vi + g[i] ≤ 0`
+    upper::Vector{MOI.ConstraintIndex{G,MOI.LessThan{T}}}
     constraint::MOI.ConstraintIndex{F,S}
 end
 
 function MOI.Bridges.Constraint.bridge_constraint(
-    ::Type{SignomialsBridge{T,S,P,F}},
+    ::Type{SignomialsBridge{T,S,P,F,G}},
     model,
     func::F,
     set,
-) where {T,S,P,F}
+) where {T,S,P,F,G}
     g = MOI.Utilities.scalarize(func)
+    odd = Int[]
+    lower = MOI.ConstraintIndex{G,MOI.LessThan{T}}[]
+    upper = MOI.ConstraintIndex{G,MOI.LessThan{T}}[]
     for i in eachindex(g)
         if any(isodd, set.α[i, :])
             vi = MOI.add_variable(model)
+            push!(odd, i)
             # vi ≤ -|g[i]|
-            MOI.Utilities.normalize_and_add_constraint(
-                model,
-                one(T) * vi - g[i],
-                MOI.LessThan(zero(T)),
+            push!(
+                lower,
+                MOI.Utilities.normalize_and_add_constraint(
+                    model,
+                    one(T) * vi - g[i],
+                    MOI.LessThan(zero(T)),
+                ),
             )
-            MOI.Utilities.normalize_and_add_constraint(
-                model,
-                one(T) * vi + g[i],
-                MOI.LessThan(zero(T)),
+            push!(
+                upper,
+                MOI.Utilities.normalize_and_add_constraint(
+                    model,
+                    one(T) * vi + g[i],
+                    MOI.LessThan(zero(T)),
+                ),
             )
             g[i] = vi
         end
@@ -44,7 +60,7 @@ function MOI.Bridges.Constraint.bridge_constraint(
         MOI.Utilities.vectorize(g),
         Cone(Signomials(set.cone.monomial), set.α),
     )
-    return SignomialsBridge{T,S,P,F}(constraint)
+    return SignomialsBridge{T,S,P,F,G}(odd, lower, upper, constraint)
 end
 
 function MOI.supports_constraint(
@@ -62,9 +78,9 @@ function MOI.Bridges.added_constrained_variable_types(
 end
 
 function MOI.Bridges.added_constraint_types(
-    ::Type{<:SignomialsBridge{T,S,P,F}},
-) where {T,S,P,F}
-    return [(F, S)]
+    ::Type{<:SignomialsBridge{T,S,P,F,G}},
+) where {T,S,P,F,G}
+    return [(F, S), (G, MOI.LessThan{T})]
 end
 
 function MOI.Bridges.Constraint.concrete_bridge_type(
@@ -72,7 +88,13 @@ function MOI.Bridges.Constraint.concrete_bridge_type(
     F::Type{<:MOI.AbstractVectorFunction},
     P::Type{Cone{Polynomials{M}}},
 ) where {T,M}
-    return SignomialsBridge{T,Cone{Signomials{M}},P,F}
+    G = MOI.Utilities.promote_operation(
+        -,
+        T,
+        MOI.ScalarAffineFunction{T},
+        MOI.Utilities.scalar_type(F),
+    )
+    return SignomialsBridge{T,Cone{Signomials{M}},P,F,G}
 end
 
 function MOI.get(
@@ -82,3 +104,26 @@ function MOI.get(
 )
     return MOI.get(model, attr, bridge.constraint)
 end
+
+# The dual of the polynomial SAGE constraint is the adjoint of the linear map
+# used in the reformulation, applied to the duals of the constraints created
+# by the bridge. For a row `i` with even exponents, `g[i]` only appears in row
+# `i` of the signomial constraint so the dual is the corresponding entry `w[i]`
+# of its dual `w`. For an odd row, `g[i]` appears with coefficient `-1` in
+# `lower[i]` and `+1` in `upper[i]` so the dual is the difference of the duals
+# of these two constraints. The result `v` satisfies `|v[i]| ≤ w[i]` for odd
+# rows, which matches the characterization of the dual of the polynomial SAGE
+# cone in terms of the dual of the signomial SAGE cone of [MCW20].
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDual,
+    bridge::SignomialsBridge,
+)
+    v = MOI.get(model, attr, bridge.constraint)
+    for (j, i) in enumerate(bridge.odd)
+        v[i] =
+            MOI.get(model, attr, bridge.upper[j]) -
+            MOI.get(model, attr, bridge.lower[j])
+    end
+    return v
+end
diff --git a/src/SAGE/optimizer.jl b/src/SAGE/optimizer.jl
new file mode 100644
index 0000000..6fffbd8
--- /dev/null
+++ b/src/SAGE/optimizer.jl
@@ -0,0 +1,185 @@
+"""
+    Optimizer{T}(solver; feasibility_tolerance = sqrt(Base.rtoldefault(T)))
+
+Optimizer computing a bound on the objective value of a polynomial
+optimization problem using its SAGE relaxation, see
+[`PolyJuMP.AbstractRelaxationOptimizer`](@ref) with the cone
+[`Polynomials`](@ref) as certificate of nonnegativity [CP16; MCW21].
+The relaxation is solved with `solver` and the bound can be queried with
+`MOI.ObjectiveBound()`. Candidate solutions are in addition recovered from
+the dual of the relaxation following [MCW21, Section 4.2]; the candidates
+are classified as feasible up to `feasibility_tolerance` (see
+`MOI.PrimalStatus`) and sorted by objective value.
+
+[CP16] Chandrasekaran, Venkat, and Parikshit Shah.
+"Relative entropy relaxations for signomial optimization."
+SIAM Journal on Optimization 26.2 (2016): 1147-1173.
+[MCW21] Murray, Riley, Venkat Chandrasekaran, and Adam Wierman.
+"Signomials and polynomial optimization via relative entropy and partial dualization."
+Mathematical Programming Computation 13 (2021): 257-295.
+https://arxiv.org/pdf/1907.00814.pdf
+"""
+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) = "PolyJuMP.SAGE"
+
+PolyJuMP.nonnegativity_cone(::Optimizer) = Polynomials()
+
+# Solve `A * z = b` over GF(2). Return `nothing` if the system is infeasible
+# and otherwise a particular solution and a basis of the nullspace of `A`.
+# This is used for the sign recovery of [MCW21, Section 4.2], see also
+# `sageopt.relaxations.symbolic_correspondences.mod2linsolve` in `sageopt`:
+# https://github.com/rileyjmurray/sageopt
+function _mod2_solve(A::Matrix{Bool}, b::Vector{Bool})
+    A = copy(A)
+    b = copy(b)
+    m, n = size(A)
+    pivot_cols = Int[]
+    for col in 1:n
+        r = length(pivot_cols) + 1
+        p = findfirst(i -> A[i, col], r:m)
+        if isnothing(p)
+            continue
+        end
+        p += r - 1
+        A[[r, p], :] = A[[p, r], :]
+        b[r], b[p] = b[p], b[r]
+        for i in 1:m
+            if i != r && A[i, col]
+                A[i, :] .⊻= A[r, :]
+                b[i] ⊻= b[r]
+            end
+        end
+        push!(pivot_cols, col)
+    end
+    if any(i -> b[i], (length(pivot_cols)+1):m)
+        return nothing
+    end
+    z = fill(false, n)
+    for (i, col) in enumerate(pivot_cols)
+        z[col] = b[i]
+    end
+    nullspace = map(setdiff(1:n, pivot_cols)) do col
+        w = fill(false, n)
+        w[col] = true
+        for (i, pivot) in enumerate(pivot_cols)
+            w[pivot] = A[i, col]
+        end
+        return w
+    end
+    return z, nullspace
+end
+
+"""
+    PolyJuMP.recover_solutions(model::Optimizer, relaxation, cref, lagrangian)
+
+Recover candidate solutions from the dual `v` of the SAGE constraint `cref`
+of the Lagrangian, following [MCW21, Section 4.2] whose reference
+implementation is `poly_solrec` in `sageopt`:
+https://github.com/rileyjmurray/sageopt/blob/master/sageopt/relaxations/poly_solution_recovery.py
+
+After normalizing `v` by its entry for the constant monomial, `v` is a
+pseudo-moment vector: at the optimum of a tight relaxation, `v[i]` is the
+monomial `monos[i]` evaluated at an optimal solution (or a convex combination
+of the moments of several optimal solutions). The magnitude of the variables
+is recovered by solving the least squares problem `α * y ≈ log.(abs.(v))` on
+the exponents `α` of the monomials with nonnegligible moment; variables
+appearing in no such monomial have all their moments negligible, hence
+magnitude zero. The signs are recovered from the linear system
+`α * z ≡ [v[i] < 0] (mod 2)` over GF(2); when the moments leave signs
+undetermined (e.g., for sign-symmetric problems), one candidate per element
+of the affine solution set is returned.
+"""
+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
+    v = MOI.get(
+        JuMP.backend(relaxation),
+        MOI.ConstraintDual(),
+        JuMP.index(cref),
+    )
+    monos = MP.monomials(lagrangian)
+    # The constant monomial is always present since `lagrangian` contains `t`
+    i0 = findfirst(iszero ∘ MP.degree, monos)
+    v /= v[i0]
+    ztol = sqrt(Base.rtoldefault(T))
+    vars = MP.variables(lagrangian)
+    rows = [i for i in eachindex(v) if i != i0 && abs(v[i]) > ztol]
+    cols = filter(eachindex(vars)) do j
+        return any(i -> !iszero(MP.degree(monos[i], vars[j])), rows)
+    end
+    F = float(T)
+    A = F[MP.degree(monos[i], vars[j]) for i in rows, j in cols]
+    b = F[log(abs(v[i])) for i in rows]
+    y = A \ b
+    mags = zeros(T, length(vars))
+    for (k, j) in enumerate(cols)
+        mags[j] = exp(y[k])
+    end
+    As = Bool[
+        isodd(MP.degree(monos[i], vars[j])) for i in rows, j in eachindex(vars)
+    ]
+    bs = Bool[v[i] < 0 for i in rows]
+    signs = _mod2_solve(As, bs)
+    if isnothing(signs)
+        # No consistent signs; only attempt nonnegative variable values,
+        # similar to the `heuristic_signs` option of `poly_solrec` in `sageopt`
+        patterns = [fill(false, length(vars))]
+    else
+        z, nullspace = signs
+        # Flipping the sign of a variable of magnitude zero gives the same
+        # solution so we do not enumerate these
+        filter!(
+            w -> any(j -> w[j] && !iszero(mags[j]), eachindex(mags)),
+            nullspace,
+        )
+        patterns = [z]
+        for w in nullspace
+            append!(patterns, [p .⊻ w for p in patterns])
+        end
+    end
+    x = MP.variables(model.model)
+    for pattern in patterns
+        values = zeros(T, length(x))
+        for (j, var) in enumerate(vars)
+            values[findfirst(isequal(var), x)] = pattern[j] ? -mags[j] : mags[j]
+        end
+        push!(
+            solutions,
+            PolyJuMP.Solution(values, model.model, model.feasibility_tolerance),
+        )
+    end
+    return solutions
+end
diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl
index bd8a942..1b57dc7 100644
--- a/src/nl_to_polynomial.jl
+++ b/src/nl_to_polynomial.jl
@@ -24,7 +24,7 @@ struct InvalidNLExpression <: Exception
 end
 
 function _to_polynomial!(d, ::Type, expr)
-    throw(
+    return throw(
         InvalidNLExpression(
             "Unexpected expression type `$(typeof(expr))` of `$expr`",
         ),
@@ -230,7 +230,7 @@ function _invalid_value(
     ::NLToPolynomial,
     attr::Union{MOI.VariablePrimal,MOI.ConstraintDual,MOI.ConstraintPrimal},
 )
-    throw(MOI.ResultIndexBoundsError(attr, 0))
+    return throw(MOI.ResultIndexBoundsError(attr, 0))
 end
 
 function MOI.get(model::NLToPolynomial, attr::MOI.AbstractModelAttribute)
diff --git a/src/optimizer.jl b/src/optimizer.jl
new file mode 100644
index 0000000..1217f7d
--- /dev/null
+++ b/src/optimizer.jl
@@ -0,0 +1,350 @@
+"""
+    abstract type AbstractPolynomialOptimizer{T} <: MOI.AbstractOptimizer end
+
+Optimizer for polynomial optimization problems storing the problem in a
+[`Model`](@ref).
+
+Subtypes should be mutable structs with the fields
+```julia
+model::PolyJuMP.Model{T}
+solve_time::Float64
+```
+and implement [`_invalidate!`](@ref) and [`_optimize!`](@ref).
+"""
+abstract type AbstractPolynomialOptimizer{T} <: MOI.AbstractOptimizer end
+
+"""
+    _invalidate!(model::AbstractPolynomialOptimizer)
+
+Invalidate the result of previous calls to `MOI.optimize!`; called when the
+polynomial optimization problem is modified.
+"""
+function _invalidate! end
+
+"""
+    _optimize!(model::AbstractPolynomialOptimizer)
+
+Solve the polynomial optimization problem stored in `model.model`; called by
+`MOI.optimize!` which records the elapsed time in `model.solve_time`.
+"""
+function _optimize! end
+
+MOI.is_empty(model::AbstractPolynomialOptimizer) = MOI.is_empty(model.model)
+
+function MOI.empty!(model::AbstractPolynomialOptimizer)
+    MOI.empty!(model.model)
+    _invalidate!(model)
+    return
+end
+
+function MOI.supports(
+    model::AbstractPolynomialOptimizer,
+    attr::MOI.AbstractModelAttribute,
+)
+    return MOI.supports(model.model, attr)
+end
+
+function MOI.set(
+    model::AbstractPolynomialOptimizer,
+    attr::MOI.AbstractModelAttribute,
+    value,
+)
+    MOI.set(model.model, attr, value)
+    _invalidate!(model)
+    return
+end
+
+function MOI.get(
+    model::AbstractPolynomialOptimizer,
+    attr::Union{
+        MOI.AbstractModelAttribute,
+        MOI.Bridges.ListOfNonstandardBridges,
+    },
+)
+    return MOI.get(model.model, attr)
+end
+
+function MOI.is_valid(model::AbstractPolynomialOptimizer, i::MOI.Index)
+    return MOI.is_valid(model.model, i)
+end
+
+function MOI.add_variable(model::AbstractPolynomialOptimizer)
+    _invalidate!(model)
+    return MOI.add_variable(model.model)
+end
+
+function MOI.supports_constraint(
+    model::AbstractPolynomialOptimizer,
+    ::Type{F},
+    ::Type{S},
+) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
+    return MOI.supports_constraint(model.model, F, S)
+end
+
+function MOI.add_constraint(
+    model::AbstractPolynomialOptimizer,
+    func::MOI.AbstractFunction,
+    set::MOI.AbstractSet,
+)
+    ci = MOI.add_constraint(model.model, func, set)
+    _invalidate!(model)
+    return ci
+end
+
+MOI.supports_incremental_interface(::AbstractPolynomialOptimizer) = true
+
+function MOI.copy_to(dest::AbstractPolynomialOptimizer, src::MOI.ModelLike)
+    return MOI.Utilities.default_copy_to(dest, src)
+end
+
+function MOI.optimize!(model::AbstractPolynomialOptimizer)
+    model.solve_time = @elapsed _optimize!(model)
+    return
+end
+
+function MOI.get(model::AbstractPolynomialOptimizer, ::MOI.SolveTimeSec)
+    return model.solve_time
+end
+
+"""
+    MultiplierMaxdegree()
+
+A constraint attribute for the maximum degree of the multiplier of the
+constraint in the certificate of nonnegativity of the Lagrangian used by an
+[`AbstractRelaxationOptimizer`](@ref).
+Increasing this degree gives a higher level of the hierarchy of relaxations,
+hence a possibly tighter objective bound at the price of a larger relaxation.
+By default, the degree `d - maxdegree(g)` is used for a constraint of
+polynomial `g` where `d` is the smallest even number larger than the maximum
+degree of the objective and constraint polynomials.
+"""
+struct MultiplierMaxdegree <: MOI.AbstractConstraintAttribute end
+
+function MOI.Bridges.Constraint.invariant_under_function_conversion(
+    ::MultiplierMaxdegree,
+)
+    return true
+end
+
+"""
+    abstract type AbstractRelaxationOptimizer{T} <: AbstractPolynomialOptimizer{T} end
+
+Optimizer computing a bound on the objective value of a polynomial
+optimization problem
+```
+min  f(x)
+s.t. g_i(x) ≥ 0
+     h_j(x) = 0
+```
+by solving the relaxation
+```
+max  t
+s.t. f - t - Σ_i σ_i * g_i - Σ_j μ_j * h_j ∈ C
+     σ_i ∈ C
+```
+where `μ_j` are free polynomials and `C` is the cone of certified nonnegative
+polynomials returned by [`nonnegativity_cone`](@ref)
+(and conversely for a `max` problem). The bound is the objective value of this
+relaxation and is returned as the `MOI.ObjectiveBound`. The degrees of the
+multipliers `σ_i` and `μ_j` are given by the [`MultiplierMaxdegree`](@ref)
+constraint attribute. Candidate primal solutions may in addition be recovered
+from the solution of the relaxation by implementing
+[`recover_solutions`](@ref); the `MOI.ResultCount` is the number of
+recovered candidates.
+
+Subtypes should be mutable structs with the fields
+```julia
+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
+```
+and implement [`nonnegativity_cone`](@ref) and [`recover_solutions`](@ref).
+"""
+abstract type AbstractRelaxationOptimizer{T} <: AbstractPolynomialOptimizer{T} end
+
+"""
+    nonnegativity_cone(model::AbstractRelaxationOptimizer)
+
+Return the set (e.g., `PolyJuMP.SAGE.Polynomials()` or
+`SumOfSquares.SOSCone()`) in which a polynomial is constrained to belong as a
+sufficient condition for its nonnegativity.
+"""
+function nonnegativity_cone end
+
+function _invalidate!(model::AbstractRelaxationOptimizer)
+    model.relaxation = nothing
+    empty!(model.solutions)
+    model.solve_time = NaN
+    return
+end
+
+"""
+    recover_solutions(
+        model::AbstractRelaxationOptimizer{T},
+        relaxation::JuMP.GenericModel{T},
+        cref::JuMP.ConstraintRef,
+        lagrangian,
+    ) where {T}
+
+Return a vector of candidate [`Solution`](@ref)s recovered from the solution
+of `relaxation`, where `cref` is the constraint of the `lagrangian` polynomial
+in the cone [`nonnegativity_cone`](@ref). Return an empty vector if no
+solution is recovered; `PolyJuMP.SAGE.Optimizer` implements the recovery from
+the dual of `cref` of [MCW21, Section 4.2].
+
+[MCW21] Murray, Riley, Venkat Chandrasekaran, and Adam Wierman.
+"Signomials and polynomial optimization via relative entropy and partial dualization."
+Mathematical Programming Computation 13 (2021): 257-295.
+https://arxiv.org/pdf/1907.00814.pdf
+"""
+function recover_solutions end
+
+function MOI.empty!(model::AbstractRelaxationOptimizer)
+    MOI.empty!(model.model)
+    empty!(model.multiplier_maxdegree)
+    _invalidate!(model)
+    return
+end
+
+function MOI.supports(
+    ::AbstractRelaxationOptimizer{T},
+    ::MultiplierMaxdegree,
+    ::Type{<:MOI.ConstraintIndex{<:ScalarPolynomialFunction{T}}},
+) where {T}
+    return true
+end
+
+function MOI.set(
+    model::AbstractRelaxationOptimizer{T},
+    ::MultiplierMaxdegree,
+    ci::MOI.ConstraintIndex{<:ScalarPolynomialFunction{T}},
+    degree::Integer,
+) where {T}
+    MOI.throw_if_not_valid(model, ci)
+    model.multiplier_maxdegree[ci] = degree
+    _invalidate!(model)
+    return
+end
+
+function MOI.get(
+    model::AbstractRelaxationOptimizer{T},
+    ::MultiplierMaxdegree,
+    ci::MOI.ConstraintIndex{<:ScalarPolynomialFunction{T}},
+) where {T}
+    return get(model.multiplier_maxdegree, ci, nothing)
+end
+
+_equalities(::SS.FullSpace) = []
+_equalities(set::SS.AbstractAlgebraicSet) = SS.equalities(set)
+_equalities(set::SS.BasicSemialgebraicSet) = _equalities(set.V)
+_inequalities(::SS.AbstractAlgebraicSet) = []
+_inequalities(set::SS.BasicSemialgebraicSet) = SS.inequalities(set)
+
+function _multiplier(relaxation::JuMP.GenericModel, x, degree)
+    basis = MB.SubBasis{MB.Monomial}(MP.monomials(x, 0:degree))
+    poly = JuMP.@variable(relaxation, variable_type = Poly(basis))
+    return MP.polynomial(poly)
+end
+
+function _optimize!(model::AbstractRelaxationOptimizer{T}) where {T}
+    pop = model.model
+    x = MP.variables(pop)
+    if pop.objective_sense == MOI.FEASIBILITY_SENSE ||
+       isnothing(pop.objective_function)
+        f = zero(PolyType{T})
+    else
+        f = pop.objective_function
+    end
+    eqs = _equalities(pop.set)
+    ineqs = _inequalities(pop.set)
+    maxdeg = max(
+        MP.maxdegree(f),
+        maximum(MP.maxdegree, eqs; init = 0),
+        maximum(MP.maxdegree, ineqs; init = 0),
+    )
+    maxdeg += isodd(maxdeg)
+    eq_degrees = Dict{Int,Int}()
+    ineq_degrees = Dict{Int,Int}()
+    for (ci, d) in model.multiplier_maxdegree
+        if ci isa MOI.ConstraintIndex{<:Any,<:MOI.EqualTo}
+            eq_degrees[ci.value] = d
+        else
+            ineq_degrees[ci.value] = d
+        end
+    end
+    relaxation = JuMP.GenericModel{T}(model.solver)
+    cone = nonnegativity_cone(model)
+    t = JuMP.@variable(relaxation, base_name = "t")
+    lagrangian = f - t
+    if pop.objective_sense == MOI.MAX_SENSE
+        lagrangian = -lagrangian
+    end
+    for (i, h) in enumerate(eqs)
+        d = get(eq_degrees, i, maxdeg - MP.maxdegree(h))
+        lagrangian -= _multiplier(relaxation, x, d) * h
+    end
+    for (i, g) in enumerate(ineqs)
+        d = get(ineq_degrees, i, maxdeg - MP.maxdegree(g))
+        σ = _multiplier(relaxation, x, d)
+        JuMP.@constraint(relaxation, σ in cone)
+        lagrangian -= σ * g
+    end
+    cref = JuMP.@constraint(relaxation, lagrangian in cone)
+    sense = pop.objective_sense == MOI.MAX_SENSE ? MOI.MIN_SENSE : MOI.MAX_SENSE
+    JuMP.set_objective(relaxation, sense, t)
+    JuMP.optimize!(relaxation)
+    model.relaxation = relaxation
+    model.solutions = recover_solutions(model, relaxation, cref, lagrangian)
+    postprocess!(model.solutions, pop, nothing)
+    return
+end
+
+function MOI.get(model::AbstractRelaxationOptimizer, ::MOI.TerminationStatus)
+    if isnothing(model.relaxation)
+        return MOI.OPTIMIZE_NOT_CALLED
+    end
+    return JuMP.termination_status(model.relaxation)
+end
+
+function MOI.get(model::AbstractRelaxationOptimizer, ::MOI.RawStatusString)
+    if isnothing(model.relaxation)
+        return "`optimize!` has not yet been called"
+    end
+    return JuMP.raw_status(model.relaxation)
+end
+
+function MOI.get(model::AbstractRelaxationOptimizer, ::MOI.ObjectiveBound)
+    return JuMP.objective_value(model.relaxation)
+end
+
+function MOI.get(model::AbstractRelaxationOptimizer, ::MOI.ResultCount)
+    return length(model.solutions)
+end
+
+function MOI.get(model::AbstractRelaxationOptimizer, attr::MOI.ObjectiveValue)
+    MOI.check_result_index_bounds(model, attr)
+    return model.solutions[attr.result_index].objective_value
+end
+
+function MOI.get(
+    model::AbstractRelaxationOptimizer,
+    attr::MOI.VariablePrimal,
+    vi::MOI.VariableIndex,
+)
+    MOI.throw_if_not_valid(model, vi)
+    MOI.check_result_index_bounds(model, attr)
+    return model.solutions[attr.result_index].values[vi.value]
+end
+
+function MOI.get(model::AbstractRelaxationOptimizer, attr::MOI.PrimalStatus)
+    if attr.result_index in 1:length(model.solutions)
+        return model.solutions[attr.result_index].status
+    end
+    return MOI.NO_SOLUTION
+end
+
+MOI.get(::AbstractRelaxationOptimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION
diff --git a/test/sage.jl b/test/sage.jl
index a73e716..9a36363 100644
--- a/test/sage.jl
+++ b/test/sage.jl
@@ -86,6 +86,196 @@ function test_domain(x, y, T, solver)
     @test_throws ErrorException @constraint(model, c3, p >= α, domain = S)
 end
 
+function test_optimizer_attributes(x, y, T, solver)
+    # We don't specify `T` to test the fallback
+    @test PolyJuMP.SAGE.Optimizer(solver) isa PolyJuMP.SAGE.Optimizer{Float64}
+    optimizer = PolyJuMP.SAGE.Optimizer{T}(solver)
+    @test MOI.get(optimizer, MOI.SolverName()) == "PolyJuMP.SAGE"
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
+    @test MOI.get(optimizer, MOI.RawStatusString()) ==
+          "`optimize!` has not yet been called"
+    @test MOI.get(optimizer, MOI.ResultCount()) == 0
+    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION
+    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.NO_SOLUTION
+    @test isnan(MOI.get(optimizer, MOI.SolveTimeSec()))
+    list = MOI.get(optimizer, MOI.Bridges.ListOfNonstandardBridges{T}())
+    @test PolyJuMP.Bridges.Constraint.ToPolynomialBridge{T} in list
+    @test PolyJuMP.Bridges.Objective.ToPolynomialBridge{T} in list
+    @test MOI.supports_incremental_interface(optimizer)
+    src = MOI.Utilities.Model{T}()
+    v = MOI.add_variable(src)
+    index_map = MOI.copy_to(optimizer, src)
+    @test MOI.is_valid(optimizer, index_map[v])
+    # The attribute getter of the optimizer is not covered through JuMP as
+    # the value is then cached by `MOI.Utilities.CachingOptimizer`
+    func = PolyJuMP.ScalarPolynomialFunction(
+        MP.polynomial(one(T) * x^2),
+        [index_map[v]],
+    )
+    ci = MOI.add_constraint(optimizer, func, MOI.GreaterThan(one(T)))
+    attr = PolyJuMP.MultiplierMaxdegree()
+    @test MOI.supports(optimizer, attr, typeof(ci))
+    @test isnothing(MOI.get(optimizer, attr, ci))
+    MOI.set(optimizer, attr, ci, 2)
+    @test MOI.get(optimizer, attr, ci) == 2
+    @test !MOI.is_empty(optimizer)
+    MOI.empty!(optimizer)
+    @test MOI.is_empty(optimizer)
+end
+
+function test_age_dual(x, y, T, solver)
+    model = Model(solver)
+    @variable(model, γ)
+    @objective(model, Max, γ)
+    # By the AM-GM inequality, `x^2 * y + x * y^2 + 1 >= 3 * x * y` for
+    # `x, y >= 0` with equality at `(1, 1)`
+    con = @constraint(
+        model,
+        x^2 * y + x * y^2 + 1 - γ * x * y in PolyJuMP.SAGE.Signomials(x * y)
+    )
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test value(γ) ≈ 3 rtol = 1e-3
+    # The dual is the vector of moments of `exp` evaluated at the optimal
+    # solution `(1, 1)`, up to the sign convention of duals of `Max` problems
+    v = MOI.get(backend(model), MOI.ConstraintDual(), JuMP.index(con))
+    @test abs.(v) ≈ ones(4) rtol = 1e-3
+end
+
+function test_optimizer_mod2_solve(x, y, T, solver)
+    # `0 * z = 1` is infeasible
+    @test isnothing(PolyJuMP.SAGE._mod2_solve(fill(false, 1, 1), [true]))
+    z, nullspace = PolyJuMP.SAGE._mod2_solve(fill(true, 1, 1), [true])
+    @test z == [true]
+    @test isempty(nullspace)
+end
+
+function test_optimizer_feasibility(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    @constraint(model, a^2 >= 1)
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_bound(model) ≈ 0 atol = 1e-4
+end
+
+function test_optimizer_unbounded(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    # Unbounded below so the relaxation is infeasible
+    @objective(model, Min, a)
+    optimize!(model)
+    @test termination_status(model) != MOI.OPTIMAL
+    @test result_count(model) == 0
+    @test primal_status(model) == MOI.NO_SOLUTION
+end
+
+function test_optimizer_not_tight(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    # `(a + 1)^2 * (a - 2)^2` whose signomial representative
+    # `a^4 - 2|a|^3 - 3a^2 - 4|a| + 4` is negative at `1` so the relaxation
+    # is not tight and the dual is not a pseudo-moment vector
+    @objective(model, Min, a^4 - 2a^3 - 3a^2 + 4a + 4)
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_bound(model) <= 1e-3
+end
+
+function test_optimizer_motzkin(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    @variable(model, b)
+    @objective(model, Min, a^4 * b^2 + a^2 * b^4 + 1 - 3 * a^2 * b^2)
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_bound(model) ≈ 0 atol = 1e-3
+    # The four minimizers `(±1, ±1)` are recovered from the dual
+    @test result_count(model) == 4
+    @test primal_status(model) == MOI.FEASIBLE_POINT
+    @test !isempty(raw_status(model))
+    @test solve_time(model) >= 0
+    for i in 1:4
+        @test abs(value(a; result = i)) ≈ 1 rtol = 1e-3
+        @test abs(value(b; result = i)) ≈ 1 rtol = 1e-3
+        @test objective_value(model; result = i) ≈ 0 atol = 1e-3
+    end
+    signs = [(value(a; result = i) > 0, value(b; result = i) > 0) for i in 1:4]
+    @test sort(signs) ==
+          [(false, false), (false, true), (true, false), (true, true)]
+end
+
+function test_optimizer_asymmetric(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    @objective(model, Min, a^2 - 2a + 3)
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_bound(model) ≈ 2 rtol = 1e-3
+    # The unique minimizer `1` is recovered from the dual; its sign
+    # tests the sign conventions of the `MOI.ConstraintDual` of the bridges
+    @test result_count(model) == 1
+    @test primal_status(model) == MOI.FEASIBLE_POINT
+    @test value(a) ≈ 1 rtol = 1e-3
+    @test objective_value(model) ≈ 2 rtol = 1e-3
+end
+
+function test_optimizer_zero_solution(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    @objective(model, Min, a^2 + 1)
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_bound(model) ≈ 1 rtol = 1e-3
+    # `a²` has zero moment in the dual so the magnitude of `a` is zero
+    @test result_count(model) == 1
+    @test primal_status(model) == MOI.FEASIBLE_POINT
+    @test value(a) ≈ 0 atol = 1e-3
+    @test objective_value(model) ≈ 1 rtol = 1e-3
+end
+
+function test_optimizer_constrained(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    @objective(model, Min, a^2)
+    @constraint(model, con, a^2 >= 1)
+    # Setting the attribute before `optimize!` covers its `MOI.supports`
+    # which is checked when the cache is copied to the optimizer
+    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 rtol = 1e-3
+    # Setting it after `optimize!` covers the direct forwarding to the
+    # attached optimizer
+    MOI.set(model, PolyJuMP.MultiplierMaxdegree(), con, 0)
+    @test MOI.get(model, PolyJuMP.MultiplierMaxdegree(), con) == 0
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_bound(model) ≈ 1 rtol = 1e-3
+    # The minimizers `±1` are recovered from the dual
+    @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
+end
+
+function test_optimizer_equality_max(x, y, T, solver)
+    model = Model(() -> PolyJuMP.SAGE.Optimizer{T}(solver))
+    @variable(model, a)
+    @objective(model, Max, 2 - a^2)
+    @constraint(model, eq, a^2 == 1)
+    # Covers the degree of a multiplier of an equality constraint
+    MOI.set(model, PolyJuMP.MultiplierMaxdegree(), eq, 0)
+    optimize!(model)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_bound(model) ≈ 1 rtol = 1e-3
+    # The maximizers `±1` are recovered from the dual
+    @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) ≈ 1 rtol = 1e-3
+end
+
 import ECOS
 const SOLVERS =
     [optimizer_with_attributes(ECOS.Optimizer, MOI.Silent() => true)]