Fix monomial_variable_index in QCQP optimizer (#108)

* Fix monomial_variable_index in QCQP optimizer

* Add test

* Use IntervalArithmetics

---------

Co-authored-by: Benoît Legat <[email protected]>

Author: odow

Project.toml

@@ -6,6 +6,7 @@ version = "0.7.1"
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
@@ -18,6 +19,7 @@ SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204"
 [compat]
 DataStructures = "0.18"
 DynamicPolynomials = "0.5"
+IntervalArithmetic = "0.22"
 JuMP = "1"
 MathOptInterface = "1"
 MultivariateBases = "0.2"

src/QCQP/MOI_wrapper.jl

@@ -112,6 +112,7 @@ function MOI.supports_constraint(
 ) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
     return MOI.supports_constraint(model.model, F, S)
 end
+
 function MOI.supports_constraint(
     model::Optimizer{T},
     ::Type{<:PolyJuMP.ScalarPolynomialFunction{T}},
@@ -131,6 +132,7 @@ function MOI.add_constraint(
 )
     return MOI.add_constraint(model.model, func, set)
 end
+
 function MOI.add_constraint(
     model::Optimizer{T},
     func::PolyJuMP.ScalarPolynomialFunction{T,P},
@@ -250,31 +252,65 @@ function _add_variables!(
     return d
 end
 
+import IntervalArithmetic
+
 function monomial_variable_index(
     model::Optimizer{T},
     d::Dict,
     div,
     mono::MP.AbstractMonomialLike,
 ) where {T}
     if !haskey(d, mono)
+        # If we don't have a variable for `mono` yet,
+        # we create one now by equal to `x * y`.
+        mono_bounds = IntervalArithmetic.interval(one(T))
+        for var in MP.variables(mono)
+            deg = MP.degree(mono, var)
+            if deg == 0
+                continue
+            end
+            vi = monomial_variable_index(model, d, div, var)
+            lb_var, ub_var = MOI.Utilities.get_bounds(model, T, vi)
+            F = float(T)
+            var_bounds = IntervalArithmetic.interval(
+                lb_var == typemin(lb_var) ? typemin(F) : float(lb_var),
+                ub_var == typemax(ub_var) ? typemax(F) : float(ub_var),
+            )
+            mono_bounds *= var_bounds^deg
+        end
         x = div[mono]
         vx = monomial_variable_index(model, d, div, x)
         y = MP.div_multiple(mono, x)
         vy = monomial_variable_index(model, d, div, y)
-        lx, ux = MOI.Utilities.get_bounds(model, T, vx)
-        ly, uy = MOI.Utilities.get_bounds(model, T, vy)
-        bounds = (lx * ly, lx * uy, ux * ly, ux * uy)
-        l = min(bounds...)
-        if vx == vy
-            l = max(l, zero(T))
+        lb = IntervalArithmetic.inf(mono_bounds)
+        ub = IntervalArithmetic.sup(mono_bounds)
+        if isfinite(lb)
+            if isfinite(ub)
+                if lb == ub
+                    set = MOI.EqualTo(T(lb))
+                else
+                    set = MOI.Interval(T(lb), T(ub))
+                end
+            else
+                set = MOI.GreaterThan(T(lb))
+            end
+        else
+            if isfinite(ub)
+                set = MOI.LessThan(T(ub))
+            else
+                set = nothing
+            end
         end
-        u = max(bounds...)
-        d[mono], _ =
-            MOI.add_constrained_variable(model.model, MOI.Interval(l, u))
-        MOI.add_constraint(
+        if isnothing(set)
+            d[mono] = MOI.add_variable(model.model)
+        else
+            d[mono], _ = MOI.add_constrained_variable(model.model, set)
+        end
+        MOI.Utilities.normalize_and_add_constraint(
             model,
             MA.@rewrite(one(T) * d[mono] - one(T) * vx * vy),
-            MOI.EqualTo(zero(T)),
+            MOI.EqualTo(zero(T));
+            allow_modify_function = true,
         )
     end
     return d[mono]
@@ -286,8 +322,9 @@ function _add_constraints(model::Optimizer, cis, index_to_var, d, div)
         set = MOI.get(model, MOI.ConstraintSet(), ci)
         func, index_to_var = _subs!(func, index_to_var)
         quad = _quad_convert(func.polynomial, d, div)
-        MOI.add_constraint(model, quad, set)
+        MOI.Utilities.normalize_and_add_constraint(model, quad, set)
     end
+    return
 end
 
 function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T}

test/qcqp.jl

@@ -157,6 +157,48 @@ function test_no_monomials(x, y, T)
     return
 end
 
+function test_scalar_constant_not_zero(x, y, T)
+    inner = Model{T}()
+    model = PolyJuMP.JuMP.GenericModel{T}() do
+        return PolyJuMP.QCQP.Optimizer{T}(MOI.Utilities.MockOptimizer(inner))
+    end
+    PolyJuMP.@variable(model, -1 <= x[1:3] <= 2)
+    PolyJuMP.@constraint(model, 4 * x[1] * x[2] * x[3] == 1)
+    PolyJuMP.@objective(model, Min, sum(x))
+    PolyJuMP.optimize!(model)
+    for (F, S) in MOI.get(inner, MOI.ListOfConstraintTypesPresent())
+        for ci in MOI.get(inner, MOI.ListOfConstraintIndices{F,S}())
+            if F isa MOI.ScalarQuadraticFunction
+                f = MOI.get(inner, MOI.ConstraintFunction(), ci)
+                @test iszero(f.constant)
+            end
+        end
+    end
+    return
+end
+
+function test_unbound_polynomial(x, y, T)
+    inner = Model{T}()
+    model = PolyJuMP.JuMP.GenericModel{T}() do
+        return PolyJuMP.QCQP.Optimizer{T}(MOI.Utilities.MockOptimizer(inner))
+    end
+    PolyJuMP.@variable(model, x >= 0)
+    PolyJuMP.@objective(model, Min, x^3)
+    PolyJuMP.optimize!(model)
+    F = MOI.VariableIndex
+    S = MOI.Interval{T}
+    @test MOI.get(inner, MOI.NumberOfConstraints{F,S}()) == 0
+    S = MOI.LessThan{T}
+    @test MOI.get(inner, MOI.NumberOfConstraints{F,S}()) == 0
+    S = MOI.GreaterThan{T}
+    @test MOI.get(inner, MOI.NumberOfConstraints{F,S}()) == 2
+    for ci in MOI.get(inner, MOI.ListOfConstraintIndices{F,S}())
+        set = MOI.get(inner, MOI.ConstraintSet(), ci)
+        @test set.lower == zero(T)
+    end
+    return
+end
+
 function runtests(x, y)
     for name in names(@__MODULE__; all = true)
         if startswith("$name", "test_")