Question about optimizing simple polynomial program

[innuo]

I have been having trouble getting my SOS polynomial program to solve. I have reproduced the issue in the the following simple optimization problem which fails with a primal infeasibility.

It is very likely that I am doing something wrong. I would greatly appreciate it if someone could set me straight.

The problem that the following function implements is

min γ
s.t.
     γ \in [0, 1]
     γ >= x
     x^2 >= 0.2 
     x \in [0, 1]

function test_fun()
    solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => false)
    model = SOSModel(solver)

    @polyvar x

    @variable(model, 0 <= γ <= 1)
    @objective(model, Min, γ)

    S = @set x >= 0 && x <=1 

    @constraint(model, x*x >= 0.2, domain = S)
    @constraint(model, x <= γ, domain = S)

    JuMP.optimize!(model)
    @show JuMP.primal_status(model)
    @show objective_value(model)
end

Which produces:

julia> test_fun()
CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj: -1.0000000e+01 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 8.07e-01 Pobj: -4.2636829e+00 Ad: 8.60e-01 Dobj:  1.2715419e+01 
Iter:  2 Ap: 7.40e-01 Pobj: -1.6522138e+00 Ad: 1.00e+00 Dobj:  1.0174006e+01 
Iter:  3 Ap: 5.19e-01 Pobj: -1.2764407e+00 Ad: 1.00e+00 Dobj: -2.0333996e+01 
Iter:  4 Ap: 1.53e-01 Pobj: -1.2314787e+00 Ad: 4.68e-01 Dobj: -4.3146688e+02 
Iter:  5 Ap: 2.46e-02 Pobj: -1.2253435e+00 Ad: 1.00e+00 Dobj: -2.5078739e+05 
Iter:  6 Ap: 8.83e-04 Pobj: -1.2247766e+00 Ad: 5.87e-04 Dobj: -1.0852500e+06 
Iter:  7 Ap: 2.49e-04 Pobj: -1.2241718e+00 Ad: 9.06e-04 Dobj: -6.8809838e+06 
Iter:  8 Ap: 1.19e-03 Pobj: -1.2103102e+00 Ad: 1.00e+00 Dobj: -4.2286503e+10 
Declaring primal infeasibility.
Success: SDP is primal infeasible
Certificate of primal infeasibility: a'*y=-1.00000e+00, ||A'(y)-Z||=2.36482e-11
JuMP.primal_status(model) = MathOptInterface.INFEASIBLE_POINT
objective_value(model) = 1.2103101735382997
1.2103101735382997

[blegat]

The syntax >= is a bit confusing. It's maybe best that you don't use SOSModel but Model instead to distable this syntax. What you are doing is the following

model = Model(solver)
@polyvar x
@variable(model, 0 <= γ <= 1)
@objective(model, Min, γ)
S = @set x >= 0 && x <= 1 
@constraint(model, x*x - 0.2 in SOSCone(), domain = S)
@constraint(model, y - x in SOSCone(), domain = S)

This searches for a value of y such that y - x is SOS over S and x^2 - 0.2 is SOS over S. Note that at x = 0 which is in S, x^2 - 0.2 = -0.2 is negative so x^2 - 0.2 cannot be SOS over S hence this is infeasible.

What you can do is as follows:

model = Model(solver)
@polyvar x y
@variable(model, lb)
@objective(model, Max, lb)
S = @set x >= 0 && x <= 1 && x^2 >= 0.2 && x <= y && y >= 0 && y <= 1
@constraint(model, y - lb in SOSCone(), domain = S, maxdegree = d)

This will give you lower bounds to the optimization problem you want to solve. The lower bound will get better as you increase d. For large enough d, you will get the minimizer, see https://jump.dev/SumOfSquares.jl/stable/generated/Polynomial%20Optimization/polynomial_optimization/

[innuo]

@blegat Thanks for your response. Indeed the <= syntax was the cause for my confusion. I had assumed that the optimization is conducted over the intersection of the constraint and the domain. Your way of specifying the constraint as a polynomial in an SOSCone makes the infeasibility clear.

As a follow-up question, is there a way to create the domain sequentially? In other words, can the domain object be passed around between different julia functions to iteratively append constraints to it? Is there an example to which you could point me?

[blegat]

Yes, you can do as follows:

julia> using SemialgebraicSets

julia> using DynamicPolynomials

julia> @polyvar x y
(x, y)

julia> S = @set x >= 0 && x <= 1
Basic semialgebraic Set defined by no equality
2 inequalities
 x ≥ 0
 -x + 1 ≥ 0

julia> S = S ∩ @set(x^2 >= 0.2)
Basic semialgebraic Set defined by no equality
3 inequalities
 x ≥ 0
 -x + 1.0 ≥ 0
 x^2 - 0.2 ≥ 0

julia> S = S ∩ @set(x <= y && y >= 0 && y <= 1)
Basic semialgebraic Set defined by no equality
6 inequalities
 x ≥ 0
 -x + 1.0 ≥ 0
 x^2 - 0.2 ≥ 0
 -x + y ≥ 0
 y ≥ 0
 -y + 1.0 ≥ 0

[odow]

Closing because this seems resolved, and there doesn't seem to be anything actionable here.