Unclear README.md

[votroto]

Hello, I was trying to use the new functionality of PolyJuMP, but I have trouble getting anything to work.

I tried to guess how to use the module based on what the readme says and this was the my best guess

using JuMP, PolyJuMP, Gurobi, DynamicPolynomials

model = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))
setpolymodule!(model, PolyJuMP.SAGE)
@polyvar x y
@variable(model, a)
@constraint(model, a * x * y^2 + y^3 >= a * x)

optimize!(model)

Apart from some typos in the snippets, it also fails with

Constraints of type MathOptInterface.VectorAffineFunction{Float64}-in-PolyJuMP.SAGE.Cone{PolyJuMP.SAGE.Polynomials{Nothing}} are not supported by the solver.

So I thought, maybe the QCQP optimizer was not tested yet and tried the same with the HC using KKT.

using HomotopyContinuation
model = Model(optimizer_with_attributes(
    PolyJuMP.KKT.Optimizer,
    "solver" => HomotopyContinuation.SemialgebraicSetsHCSolver(),
))

but this fails with the same error.

So I tried switching from the PolyJuMP.SAGE to SumOfSquares

using JuMP, PolyJuMP, Gurobi, DynamicPolynomials, SumOfSquares, HomotopyContinuation

model = Model(optimizer_with_attributes(
    PolyJuMP.KKT.Optimizer,
    "solver" => HomotopyContinuation.SemialgebraicSetsHCSolver(),
))
setpolymodule!(model, SumOfSquares)
@polyvar x y
@variable(model, a)
@constraint(model, a * x * y^2 + y^3 >= a * x)

optimize!(model)

but this also fails with pretty much the same error

Constraints of type MathOptInterface.VectorAffineFunction{Float64}-in-SumOfSquares.SOSPolynomialSet{FullSpace, Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}} are not supported by the solver.

I tried to switching to the inline constraint set

@constraint(model, a * x * y^2 + y^3 - a * x in SumOfSquares.SOSCone())

same thing...

So, maybe I misunderstood what the package does and it only handles quadratic things without reformulations. Then I tried a sanity check, which Gurobi can handle on its own.

using JuMP, PolyJuMP, Gurobi

model = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))
@variable(model, a)
@constraint(model, a^2 <= 1)

optimize!(model)

but that also fails

LoadError: MethodError: no method matching decompose(::Nothing)

So maybe it strictly requires something non-convex, or maybe an objective function?

using JuMP, PolyJuMP, Gurobi

model = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))
@variable(model, a)
@variable(model, b)
@constraint(model, a * b <= 1)
@objective(model, Max, a + b)
optimize!(model)

nope!

LoadError: MethodError: no method matching decompose(::Nothing)

Ok, I'm out of ideas. Is it possible to give a minimal working example in the readme of how the package is to be used, please?

Thank you.

Versions

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.10.0 (2023-12-25)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

(nonconvex) pkg> st
Project nonconvex v0.1.0
Status `~/projects/nonconvex/Project.toml`
  [7c1d4256] DynamicPolynomials v0.5.3
  [2e9cd046] Gurobi v1.2.1
  [f213a82b] HomotopyContinuation v2.9.2
  [4076af6c] JuMP v1.18.1
  [ddf597a6] PolyJuMP v0.7.1
  [4b9e565b] SumOfSquares v0.7.3

[blegat]

Thanks for trying out these new features.
The QCQP and KKT optimizers, do not handle inequalities with @polyvar (which means the polynomial should be nonnegative for all values of these @polyvar). Instead, all variables should be JuMP variables.
So try

using JuMP, PolyJuMP, Gurobi, DynamicPolynomials

model = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))
a = 2
@variable(model, x)
@variable(model, y)
@constraint(model, a * x * y^2 + y^3 >= a * x)
@objective(model, Min, x + y)
optimize!(model)

or try this one with KKT.

For SAGE, it's indeed constraints with @polyvar but you should use a conic solver like ECOS, not QCQP or KKT

[votroto]

I tried the snippet, but it does not work work for me.

PolyJuMP.InvalidNLExpression("Unexpected expression type `MathOptInterface.ScalarAffineFunction{Float64}` of `0.0 + 2.0 MOI.VariableIndex(1)`")

But this does, which is really promising:

using JuMP, PolyJuMP, Gurobi, DynamicPolynomials

model = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))
@variable(model, x)
@variable(model, y)
@constraint(model, x + y == 1)
@constraint(model, x >= 0)
@constraint(model, y >= 0)
@objective(model, Min, x^3 + y)
optimize!(model)

Unfortunately changing the objective to x^2+y breaks it again.

[votroto]

Another weird thing is the interaction with variables constrained on creation.

If I formulate the problem manually, everything works.

model = Model(Gurobi.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
@variable(model, xx)

@constraint(model, xx == x*x)
@constraint(model, x + y == 1)

@objective(model, Min, xx*x + y)
optimize!(model)

but the equivalent via QCQP, it fails with LoadError: Gurobi Error 10003: Element 0 of a double array is Nan.

model = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))
@variable(model, x >= 0)
@variable(model, y >= 0)

@constraint(model, x + y == 1)

@objective(model, Min, x^3 + y)
optimize!(model)

My gurobi version

Gurobi command line optimizer, version 11.0.0 build v11.0.0rc2 (linux64 - "Arch Linux")

[votroto]

I had no idea that it is now possible to formulate expressions such as x^3 in JuMP variables. That's brilliant!

Since I'm already out of topic for this entire issue... sorry... here goes:

The following example works as expected for SCS and SumOfSquares, but reports infeasible for ECOS with PolyJuMP.SAGE.

using JuMP, PolyJuMP, DynamicPolynomials, ECOS, SumOfSquares

@polyvar x y
p = x^3 - x^2 + 2x*y -y^2 + y^3
S = @set x >= 0 && y >= 0 && x + y >= 1

model = Model(ECOS.Optimizer)
setpolymodule!(model, PolyJuMP.SAGE)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c3, p >= α, domain = S)
optimize!(model)

[blegat]

I had no idea that it is now possible to formulate expressions such as x^3 in JuMP variables. That's brilliant!

Thanks! The goal now is to also add these solvers for SumOfSquares and SAGE that automatically rewrite it by moving the JuMP decision variables into DynamicPolynomials variables and then add the α JuMP decision variables for which the value will be reported as objective_bound

Thanks for reporting the bug, the domain keyword is not supported, this should be an error