Constraints on Positive Moment Sequences

using DynamicPolynomials, MultivariateSeries, MomentPolynomialOpt
using JuMP, MosekTools; mpo_optimizer(Mosek.Optimizer, "QUIET" => true);

We define a moment model in variables $x,y$:

X = @polyvar x y

M = MOM.Model()
d = 10

10

We define now a moment sequence $\mu_1$ with is Positive Semi-Definite (.i.e non-negative on squares) and non-negative on the unit ball, truncated in degree $2d$:

mu1 = moments(M,X,2*d,:PSD)

g1 = 1-x^2-y^2

# p1 * mu >= 0
MOM.add_constraint_nneg(M, mu1, g1);

We consider a second moment sequence $\mu_2$ which is PSD and non-negative on the bux $\mathbb{B}=[-1,1]^2$:

mu2 = moments(M, X, 2*d, :PSD)

q1 = 1-x^2
q2 = 1-y^2

# q1 * mu_2 >= 0, q2 * mu_2 >=0
MOM.add_constraint_nneg(M, mu2, q1)
MOM.add_constraint_nneg(M, mu2, q2);

Now, we impose that $\mu_1+\mu_2= \lambda$ where $\lambda$ is the Lebesgue measure on the box $\mathbb{B}$. It translates into moment constraintes:

function lebesgue(m)
    e = exponents(m); 
    return ((1-(-1)^(e[1]+1))/(e[1]+1))*((1-(-1)^(e[2]+1))/(e[2]+1))
end

L = monomials(X, 0:2*d) 
for m in L
    @constraint(M, mmt(mu1, m) + mmt(mu2,m) - lebesgue(m) == 0)
end

Finally, we optimize the total mass of $\mu_1$ to get the volume of $S$.

# sup  <mu_1,1>  
@objective(M, Max, mmt(mu1,1) )

v, M = optimize(M)

(3.5833175500047845, A JuMP Model
Maximization problem with:
Variables: 462
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 231 constraints
`Vector{AffExpr}`-in-`MathOptInterface.PositiveSemidefiniteConeTriangle`: 5 constraints
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: Dual model with Mosek attached
Names registered in the model: mu, type)