JuMP & SumOfSquares generate SDPs with ~18x more constraints than YALMIP

Specific Domains Optimization (Mathematical)
15

I checked the code in your top post. I detailed the issue in https://github.com/jump-dev/SumOfSquares.jl/issues/205.
In short, SumOfSquares uses the standard conic form for SumOfSquares constraints.
In your case, you should use the geometric conic form as discussed in https://github.com/jump-dev/SumOfSquares.jl/issues/205.
As your polynomials are quadratic, you can in fact use JuMP directly, e.g.,

model = Model(dual_optimizer(Mosek.Optimizer))

# Variables
@variable(model,H[1:n,1:n],PSD); # Matrix parametrizing V
@polyvar(x[1:n]) # State
@polyvar(u[1:m]) # Control

# Lyapunov candidate V
V = x'*H*x;

# constraints u.*(u-W1*x) == 0
h = u .* (u - W1 * x);

Vp = subs(V, x => A*x + B*u)

# Constraints

using LinearAlgebra
function matrix(poly)
    M = Matrix{AffExpr}(undef, n + m, n + m)
    v = [x; u]
    for i in 1:(n + m)
        for j in 1:(n + m)
            coef = DynamicPolynomials.coefficient(poly, v[i] * v[j])
            if i == j
                M[i, j] = coef
            else
                M[j, i] = M[i, j] = coef / 2
            end
        end
    end
    return Symmetric(M)
end

# (V - Vp - x'*x) on set (h == 0)
@variable(model, R[1:m])
@constraint(model, matrix(V - Vp - x'*x - sum(R[i] * h[i] for i in 1:m)) in PSDCone());

@variable(model, t)
@constraint(model, [t; vec(H)] in SecondOrderCone())
@objective(model, Min, t)

JuMP.optimize!(model)
@show termination_status(model)
@show objective_value(model)
5 Likes