Extracting minimizers

Adapted from: Example 6.23 of [L09]

[L09] Laurent, Monique. Sums of squares, moment matrices and optimization over polynomials. Emerging applications of algebraic geometry (2009): 157-270.

Introduction

Consider the polynomial optimization problem [L09, Example 6.23] of minimizing the linear function $-x_1 - x_2$ over the basic semialgebraic set defined by the inequalities $x_2 \le 2x_1^4 - 8x_1^3 + 8x_1^2 + 2$, $x_2 \le 4x_1^4 - 32x_1^3 + 88x_1^2 - 96x_1 + 36$ and the box constraints $0 \le x_1 \le 3$ and $0 \le x_2 \le 4$, World Scientific, 2009.

using DynamicPolynomials
@polyvar x[1:2]
p = -sum(x)
using SumOfSquares
K = @set x[1] >= 0 && x[1] <= 3 && x[2] >= 0 && x[2] <= 4 && x[2] <= 2x[1]^4 - 8x[1]^3 + 8x[1]^2 + 2 && x[2] <= 4x[1]^4 - 32x[1]^3 + 88x[1]^2 - 96x[1] + 36

Basic semialgebraic Set defined by no equality
6 inequalities
 x[1] ≥ 0
 3 - x[1] ≥ 0
 x[2] ≥ 0
 4 - x[2] ≥ 0
 2 - x[2] + 8*x[1]^2 - 8*x[1]^3 + 2*x[1]^4 ≥ 0
 36 - x[2] - 96*x[1] + 88*x[1]^2 - 32*x[1]^3 + 4*x[1]^4 ≥ 0

We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.

import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)

MathOptInterface.OptimizerWithAttributes(CSDP.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.Silent() => true])

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 function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.

function solve(d)
    model = SOSModel(solver)
    @variable(model, α)
    @objective(model, Max, α)
    @constraint(model, c, p >= α, domain = K, maxdegree = d)
    optimize!(model)
    println(solution_summary(model))
    return model
end

solve (generic function with 1 method)

The first level of the hierarchy gives a lower bound of -7`

model4 = solve(4)

A JuMP Model
Maximization problem with:
Variable: 1
Objective function type: VariableRef
`Vector{AffExpr}`-in-`SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Int64, DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}, FullSpace}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Putinar{SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}`: 1 constraint
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: CSDP
Names registered in the model: c, α

The second level improves the lower bound

model5 = solve(5)

A JuMP Model
Maximization problem with:
Variable: 1
Objective function type: VariableRef
`Vector{AffExpr}`-in-`SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Int64, DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}, FullSpace}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Putinar{SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}`: 1 constraint
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: CSDP
Names registered in the model: c, α

The third level finds the optimal objective value as lower bound...

model7 = solve(7)

A JuMP Model
Maximization problem with:
Variable: 1
Objective function type: VariableRef
`Vector{AffExpr}`-in-`SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Int64, DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}, FullSpace}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Putinar{SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}`: 1 constraint
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: CSDP
Names registered in the model: c, α

...and proves it by exhibiting the minimizer.

ν7 = moment_matrix(model7[:c])
η = atomic_measure(ν7, 1e-3) # Returns nothing as the dual is not atomic

Atomic measure on the variables x[1], x[2] with 1 atoms:
 at [2.32952019617943, 3.1784930775648848] with weight 0.9999999990880203

We can indeed verify that the objective value at x_opt is equal to the lower bound.

x_opt = η.atoms[1].center
p(x_opt)

-5.508013273744314

This page was generated using Literate.jl.