Solving a polynomial optimization problem
using DynamicPolynomials, MomentToolsX = @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; optimizer = CSDP.Optimizer
v, M = maximize(x1, [e1, e2], [p1,p2], X, 4, optimizer)
vCSDP 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.3773383e+01 Ad: 6.75e-01 Dobj: -3.1219527e-01
Iter: 3 Ap: 1.46e-01 Pobj: 1.6205657e+02 Ad: 1.91e-01 Dobj: -7.3806838e-01
Iter: 4 Ap: 7.10e-01 Pobj: -7.0126415e+00 Ad: 6.50e-01 Dobj: -1.5156365e+00
Iter: 5 Ap: 8.65e-01 Pobj: -7.5100797e+00 Ad: 8.90e-01 Dobj: -1.3959313e+00
Iter: 6 Ap: 7.93e-01 Pobj: -3.6855309e+00 Ad: 7.29e-01 Dobj: -1.4143253e+00
Iter: 7 Ap: 7.37e-01 Pobj: -2.5848103e+00 Ad: 7.43e-01 Dobj: -1.4136246e+00
Iter: 8 Ap: 5.35e-01 Pobj: -2.1333028e+00 Ad: 7.75e-01 Dobj: -1.4141761e+00
Iter: 9 Ap: 6.34e-01 Pobj: -1.7466820e+00 Ad: 7.25e-01 Dobj: -1.4141777e+00
Iter: 10 Ap: 7.08e-01 Pobj: -1.5388419e+00 Ad: 7.76e-01 Dobj: -1.4142100e+00
Iter: 11 Ap: 7.56e-01 Pobj: -1.4544602e+00 Ad: 7.16e-01 Dobj: -1.4142117e+00
Iter: 12 Ap: 6.99e-01 Pobj: -1.4305242e+00 Ad: 7.56e-01 Dobj: -1.4142132e+00
Iter: 13 Ap: 5.90e-01 Pobj: -1.4222141e+00 Ad: 7.46e-01 Dobj: -1.4142134e+00
Iter: 14 Ap: 6.66e-01 Pobj: -1.4173958e+00 Ad: 9.44e-01 Dobj: -1.4142135e+00
Iter: 15 Ap: 8.91e-01 Pobj: -1.4146634e+00 Ad: 1.00e+00 Dobj: -1.4142135e+00
Iter: 16 Ap: 9.38e-01 Pobj: -1.4142490e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 17 Ap: 9.85e-01 Pobj: -1.4142149e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 18 Ap: 1.00e+00 Pobj: -1.4142138e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 19 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 20 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 9.37e-01 Dobj: -1.4142136e+00
Iter: 21 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 8.42e-01 Dobj: -1.4142136e+00
Success: SDP solved
Primal objective value: -1.4142136e+00
Dual objective value: -1.4142136e+00
Relative primal infeasibility: 6.98e-11
Relative dual infeasibility: 7.51e-10
Real Relative Gap: 4.20e-10
XZ Relative Gap: 7.21e-09
DIMACS error measures: 6.98e-11 0.00e+00 5.92e-09 0.00e+00 4.20e-10 7.21e-09
1.414213562364079The 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 Array{Float64,2}:
1.41421 1.41421 1.41421
-1.73205 1.41416 1.73198Each 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.22888491281955842, 0.4931176123147139, 0.277997474859004], [1.414213562391005 1.4142135621730105 1.4142135621419918; -1.7320521660841413 1.4141606157582571 1.7319755718907128])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 Array{Tuple{AbstractPolynomialLike{Int64},String},1}:
(x1, "sup")
(x1² - 2, "=0")
(x1x2³ - 3x1x2 - 2x2² + 6, "=0")
(x1, ">=0")
(-x2 + 2, ">=0")v, M = optimize(pop, X, 4, 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.3773383e+01 Ad: 6.75e-01 Dobj: -3.1219527e-01
Iter: 3 Ap: 1.46e-01 Pobj: 1.6205657e+02 Ad: 1.91e-01 Dobj: -7.3806838e-01
Iter: 4 Ap: 7.10e-01 Pobj: -7.0126412e+00 Ad: 6.50e-01 Dobj: -1.5156365e+00
Iter: 5 Ap: 8.65e-01 Pobj: -7.5100797e+00 Ad: 8.90e-01 Dobj: -1.3959313e+00
Iter: 6 Ap: 7.93e-01 Pobj: -3.6855309e+00 Ad: 7.29e-01 Dobj: -1.4143253e+00
Iter: 7 Ap: 7.37e-01 Pobj: -2.5848117e+00 Ad: 7.43e-01 Dobj: -1.4136246e+00
Iter: 8 Ap: 5.35e-01 Pobj: -2.1333020e+00 Ad: 7.75e-01 Dobj: -1.4141762e+00
Iter: 9 Ap: 6.34e-01 Pobj: -1.7466786e+00 Ad: 7.25e-01 Dobj: -1.4141777e+00
Iter: 10 Ap: 7.08e-01 Pobj: -1.5388415e+00 Ad: 7.76e-01 Dobj: -1.4142100e+00
Iter: 11 Ap: 7.56e-01 Pobj: -1.4544598e+00 Ad: 7.16e-01 Dobj: -1.4142117e+00
Iter: 12 Ap: 6.99e-01 Pobj: -1.4305245e+00 Ad: 7.56e-01 Dobj: -1.4142132e+00
Iter: 13 Ap: 5.90e-01 Pobj: -1.4222141e+00 Ad: 7.46e-01 Dobj: -1.4142134e+00
Iter: 14 Ap: 6.66e-01 Pobj: -1.4173956e+00 Ad: 9.44e-01 Dobj: -1.4142135e+00
Iter: 15 Ap: 8.91e-01 Pobj: -1.4146603e+00 Ad: 1.00e+00 Dobj: -1.4142135e+00
Iter: 16 Ap: 9.38e-01 Pobj: -1.4142481e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 17 Ap: 9.85e-01 Pobj: -1.4142148e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 18 Ap: 1.00e+00 Pobj: -1.4142138e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 19 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 1.00e+00 Dobj: -1.4142136e+00
Iter: 20 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 9.37e-01 Dobj: -1.4142136e+00
Iter: 21 Ap: 1.00e+00 Pobj: -1.4142136e+00 Ad: 8.41e-01 Dobj: -1.4142136e+00
Success: SDP solved
Primal objective value: -1.4142136e+00
Dual objective value: -1.4142136e+00
Relative primal infeasibility: 7.60e-11
Relative dual infeasibility: 7.53e-10
Real Relative Gap: 4.04e-10
XZ Relative Gap: 7.23e-09
DIMACS error measures: 7.60e-11 0.00e+00 5.94e-09 0.00e+00 4.04e-10 7.23e-09
(1.4142135623642407,
A Moment program with:
A JuMP Model
Maximization problem with:
Variables: 45
Objective function type: JuMP.GenericAffExpr{Float64,JuMP.VariableRef}
`JuMP.GenericAffExpr{Float64,JuMP.VariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 44 constraints
`Array{JuMP.VariableRef,1}`-in-`MathOptInterface.PositiveSemidefiniteConeTriangle`: 1 constraint
`Array{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},1}`-in-`MathOptInterface.PositiveSemidefiniteConeSquare`: 2 constraints
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: Dual model with CSDP attached
Names registered in the model: basis, degree, dual, index, moments, monomials, nu, variables, y)get_minimizers(M)2×3 Array{Float64,2}:
1.41421 1.41421 1.41421
-1.73205 1.41416 1.73198