Solving a polynomial optimization problem

using DynamicPolynomials, MomentPolynomialOpt

X  = @polyvar x1 x2

e1 = x1^2-2
e2 = (x2^2-3)*(x1*x2-2)

p1 = x1
p2 = 2-x2;

We are looking for the points with maximal $x_1$ in the set $e_{1}=e_{2}=0$ such that $p_1\geq 0$, $p_2\geq 0$.

We solve a SDP relaxation of order $d=4$, where the variables of the underlying convex optimization problem are the moments of order $\le 2d$ in the variables $x_1, x_2$.

using CSDP;
v, M = maximize(x1, [e1, e2], [p1,p2], X, 4, CSDP.Optimizer)
v

CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 7.96e-01 Pobj:  1.3512389e+00 Ad: 7.77e-01 Dobj: -1.6779112e-01 
Iter:  2 Ap: 5.89e-01 Pobj:  1.3774675e+01 Ad: 6.75e-01 Dobj: -3.1220844e-01 
Iter:  3 Ap: 1.45e-01 Pobj:  1.6171643e+02 Ad: 1.90e-01 Dobj: -7.3834519e-01 
Iter:  4 Ap: 7.11e-01 Pobj: -7.5743600e+00 Ad: 6.50e-01 Dobj: -1.5160115e+00 
Iter:  5 Ap: 8.65e-01 Pobj: -7.6704896e+00 Ad: 8.84e-01 Dobj: -1.3959650e+00 
Iter:  6 Ap: 7.98e-01 Pobj: -3.7100039e+00 Ad: 7.34e-01 Dobj: -1.4143632e+00 
Iter:  7 Ap: 7.35e-01 Pobj: -2.6064472e+00 Ad: 7.48e-01 Dobj: -1.4135877e+00 
Iter:  8 Ap: 5.35e-01 Pobj: -2.1446579e+00 Ad: 8.01e-01 Dobj: -1.4142031e+00 
Iter:  9 Ap: 6.65e-01 Pobj: -1.7242201e+00 Ad: 7.02e-01 Dobj: -1.4141766e+00 
Iter: 10 Ap: 6.91e-01 Pobj: -1.5362595e+00 Ad: 8.31e-01 Dobj: -1.4142130e+00 
Iter: 11 Ap: 7.66e-01 Pobj: -1.4508474e+00 Ad: 6.82e-01 Dobj: -1.4142118e+00 
Iter: 12 Ap: 6.46e-01 Pobj: -1.4304738e+00 Ad: 7.98e-01 Dobj: -1.4142134e+00 
Iter: 13 Ap: 6.36e-01 Pobj: -1.4212252e+00 Ad: 7.70e-01 Dobj: -1.4142135e+00 
Iter: 14 Ap: 7.07e-01 Pobj: -1.4166180e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 15 Ap: 9.27e-01 Pobj: -1.4144197e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 16 Ap: 9.49e-01 Pobj: -1.4142278e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 17 Ap: 9.44e-01 Pobj: -1.4142151e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 18 Ap: 1.00e+00 Pobj: -1.4142137e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 19 Ap: 9.95e-01 Pobj: -1.4142136e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 20 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 9.25e-01 Dobj: -1.4142136e+00 
Iter: 21 Ap: 6.57e-01 Pobj: -1.4142136e+00 Ad: 6.03e-01 Dobj: -1.4142136e+00 
Success: SDP solved
Primal objective value: -1.4142136e+00 
Dual objective value: -1.4142136e+00 
Relative primal infeasibility: 1.17e-09 
Relative dual infeasibility: 6.11e-10 
Real Relative Gap: 4.77e-09 
XZ Relative Gap: 4.20e-09 
DIMACS error measures: 1.17e-09 0.00e+00 4.81e-09 0.00e+00 4.77e-09 4.20e-09





1.4142135623764274

The output of the function maximize is the optimal value v and the optimization model M.

The points which reach the optimal value, can be obtained as follows:

Xi = get_minimizers(M)

2×3 Matrix{Float64}:
 1.41421  1.41421   1.41421
 1.73196  1.41415  -1.73205

Each column of this matrix represents a point. It is an $n\times r$ matrix, where $n$ is the number of coordinates in X and $r$ is the number of points.

The weighted sum of Dirac measures associated to the optimal moment sequence can be obtained as follows:

w, Xi = get_measure(M)

([0.2846422969212594, 0.48611125587805604, 0.22924644720334292], [1.4142135622523302 1.4142135622777432 1.414213562369705; 1.7319614907125995 1.414148281759339 -1.7320524561055801])

w is the vector of weights and Xi is the matrix of points, that is support of the measure $\mu=\sum_i \omega_i \delta_{\Xi_i}$.

Here is another way to solve it. We describe it as a Polynomial Optimization Problem and use the function optimize:

pop = [(x1, "sup"), (e1,"=0"),(e2 ,"=0"),(p1,">=0"),(p2,">=0")]

5-element Vector{Tuple{AbstractPolynomialLike{Int64}, String}}:
 (x1, "sup")
 (-2 + x1², "=0")
 (6 - 2x2² - 3x1x2 + x1x2³, "=0")
 (x1, ">=0")
 (2 - x2, ">=0")

v, M = optimize(pop, X, 4, CSDP.Optimizer)

CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 7.96e-01 Pobj:  1.3512389e+00 Ad: 7.77e-01 Dobj: -1.6779112e-01 
Iter:  2 Ap: 5.89e-01 Pobj:  1.3774675e+01 Ad: 6.75e-01 Dobj: -3.1220844e-01 
Iter:  3 Ap: 1.45e-01 Pobj:  1.6171643e+02 Ad: 1.90e-01 Dobj: -7.3834519e-01 
Iter:  4 Ap: 7.11e-01 Pobj: -7.5743600e+00 Ad: 6.50e-01 Dobj: -1.5160115e+00 
Iter:  5 Ap: 8.65e-01 Pobj: -7.6704896e+00 Ad: 8.84e-01 Dobj: -1.3959650e+00 
Iter:  6 Ap: 7.98e-01 Pobj: -3.7100039e+00 Ad: 7.34e-01 Dobj: -1.4143632e+00 
Iter:  7 Ap: 7.35e-01 Pobj: -2.6064472e+00 Ad: 7.48e-01 Dobj: -1.4135877e+00 
Iter:  8 Ap: 5.35e-01 Pobj: -2.1446579e+00 Ad: 8.01e-01 Dobj: -1.4142031e+00 
Iter:  9 Ap: 6.65e-01 Pobj: -1.7242201e+00 Ad: 7.02e-01 Dobj: -1.4141766e+00 
Iter: 10 Ap: 6.91e-01 Pobj: -1.5362595e+00 Ad: 8.31e-01 Dobj: -1.4142130e+00 
Iter: 11 Ap: 7.66e-01 Pobj: -1.4508474e+00 Ad: 6.82e-01 Dobj: -1.4142118e+00 
Iter: 12 Ap: 6.46e-01 Pobj: -1.4304738e+00 Ad: 7.98e-01 Dobj: -1.4142134e+00 
Iter: 13 Ap: 6.36e-01 Pobj: -1.4212252e+00 Ad: 7.70e-01 Dobj: -1.4142135e+00 
Iter: 14 Ap: 7.07e-01 Pobj: -1.4166180e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 15 Ap: 9.27e-01 Pobj: -1.4144197e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 16 Ap: 9.49e-01 Pobj: -1.4142278e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 17 Ap: 9.44e-01 Pobj: -1.4142151e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 18 Ap: 1.00e+00 Pobj: -1.4142137e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 19 Ap: 9.95e-01 Pobj: -1.4142136e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00 
Iter: 20 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 9.25e-01 Dobj: -1.4142136e+00 
Iter: 21 Ap: 6.57e-01 Pobj: -1.4142136e+00 Ad: 6.03e-01 Dobj: -1.4142136e+00 
Success: SDP solved
Primal objective value: -1.4142136e+00 
Dual objective value: -1.4142136e+00 
Relative primal infeasibility: 1.17e-09 
Relative dual infeasibility: 6.11e-10 
Real Relative Gap: 4.77e-09 
XZ Relative Gap: 4.20e-09 
DIMACS error measures: 1.17e-09 0.00e+00 4.81e-09 0.00e+00 4.77e-09 4.20e-09





(1.4142135623764274, A JuMP Model
Maximization problem with:
Variables: 45
Objective function type: JuMP.VariableRef
`JuMP.AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 44 constraints
`Vector{JuMP.AffExpr}`-in-`MathOptInterface.PositiveSemidefiniteConeTriangle`: 3 constraints
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: Dual model with CSDP attached
Names registered in the model: mu, type)

get_minimizers(M)

2×3 Matrix{Float64}:
  1.41421  1.41421  1.41421
 -1.73205  1.41415  1.73196