Hull: evaulate_at_zero function

ext/InfiniteDisjunctiveProgramming.jl

@@ -38,6 +38,33 @@ function DP.requires_disaggregation(vref::InfiniteOpt.GeneralVariableRef)
     return !_is_parameter(vref)
 end
 
+# Value of a reference with every decision variable at zero. Infinite
+# parameters and parameter functions stay symbolic, so transcription
+# evaluates the perspective constant f(0) at each support.
+DP.evaluate_at_zero(vref::InfiniteOpt.GeneralVariableRef) =
+    _evaluate_at_zero(vref, InfiniteOpt.dispatch_variable_ref(vref))
+_evaluate_at_zero(vref, ::InfiniteOpt.DecisionVariableRef) = 0.0
+_evaluate_at_zero(vref, pref::InfiniteOpt.FiniteParameterRef) =
+    convert(Float64, JuMP.parameter_value(pref))
+_evaluate_at_zero(vref, ::Union{InfiniteOpt.IndependentParameterRef,
+    InfiniteOpt.DependentParameterRef, InfiniteOpt.ParameterFunctionRef}) =
+    vref
+function _evaluate_at_zero(vref, ::InfiniteOpt.MeasureRef)
+    error("Hull cannot reformulate a nonlinear disjunct constraint " *
+        "containing the measure `$vref`: its argument cannot be " *
+        "evaluated with the decision variables at zero. Use BigM or " *
+        "MBM for this constraint.")
+end
+# Replace every reference by its value at zero; the result is a number
+# unless infinite parameters remain, then it is an expression in them
+function DP.evaluate_at_zero(
+    expr::JuMP.GenericNonlinearExpr{InfiniteOpt.GeneralVariableRef}
+    )
+    expr = InfiniteOpt.map_expression(DP.evaluate_at_zero, expr)
+    isempty(InfiniteOpt.parameter_refs(expr)) || return expr
+    return JuMP.value(DP.evaluate_at_zero, expr)
+end
+
 function DP.VariableProperties(vref::InfiniteOpt.GeneralVariableRef)
     info = DP.get_variable_info(vref)
     name = JuMP.name(vref)

src/hull.jl

@@ -281,16 +281,47 @@ function reformulate_disjunct_constraint(
     reform_con = JuMP.build_constraint(error, new_func, con.set)
     return [reform_con]
 end
+################################################################################
+#                       NONLINEAR EXPRESSION AT ZERO
+################################################################################
+
+"""
+    evaluate_at_zero(
+        expr::Union{Number, JuMP.AbstractJuMPScalar}
+        )::Union{Float64, JuMP.AbstractJuMPScalar}
+
+Return `expr` with every decision variable set to zero, the constant
+\$f(0)\$ of the perspective function in the hull reformulation of a
+nonlinear disjunct constraint. Expressions are evaluated by
+`JuMP.value(evaluate_at_zero, expr)`, so every variable leaf comes back
+through this function. Extensions with their own
+`JuMP.AbstractVariableRef` subtype add a method for it, since what
+"zero" means depends on the reference kind, and a method for nonlinear
+expressions over their subtype, since JuMP's nonlinear walk only
+evaluates its own variable type. An extension may return an expression
+in quantities that are not decision variables, such as infinite
+parameters; it is then substituted into the perspective constraint
+as is.
+"""
+evaluate_at_zero(c::Number) = convert(Float64, c)
+evaluate_at_zero(::JuMP.GenericVariableRef) = 0.0
+evaluate_at_zero(expr::JuMP.AbstractJuMPScalar) =
+    JuMP.value(evaluate_at_zero, expr)
+
 function reformulate_disjunct_constraint(
     model::JuMP.AbstractModel, 
     con::JuMP.ScalarConstraint{T, S}, 
     bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
     method::_Hull
 ) where {T <: JuMP.GenericNonlinearExpr, S <: Union{_MOI.LessThan, _MOI.GreaterThan, _MOI.EqualTo}}
     con_func = _disaggregate_nl_expression(model, con.func, bvref, method)
-    con_func0 = JuMP.value(v -> 0.0, con.func)
-    if isinf(con_func0)
-        error("Operator `$(con.func.head)` is not defined at 0, causing the perspective function on the Hull reformulation to fail.")
+    con_func0 = evaluate_at_zero(con.func)
+    if con_func0 isa Number && isinf(con_func0)
+        error("Cannot apply Hull to the nonlinear disjunct constraint on " *
+            "`$(con.func)`: it evaluates to $(con_func0) at zero, so the " *
+            "perspective constant f(0) is not finite (for example " *
+            "`log(x)` or `1/x`). Shift the variables so the function is " *
+            "finite at zero, or use BigM or MBM.")
     end
     ϵ = method.value
     set_value = _set_value(con.set)
@@ -307,9 +338,14 @@ function reformulate_disjunct_constraint(
     con_func = JuMP.@expression(model, [i=1:con.set.dimension],
         _disaggregate_nl_expression(model, con.func[i], bvref, method)
     )
-    con_func0 = JuMP.value.(v -> 0.0, con.func)
-    if any(isinf.(con_func0))
-        error("At least of of the operators `$([func.head for func in con.func])` is not defined at 0, causing the perspective function on the Hull reformulation to fail.")
+    con_func0 = evaluate_at_zero.(con.func)
+    rows = findall(v -> v isa Number && isinf(v), con_func0)
+    if !isempty(rows)
+        error("Cannot apply Hull to the vector disjunct constraint on " *
+            "`$(con.func)`: rows $(rows) evaluate to $(con_func0[rows]) " *
+            "at zero, so the perspective constant f(0) is not finite (for " *
+            "example `log(x)` or `1/x`). Shift the variables so the " *
+            "function is finite at zero, or use BigM or MBM.")
     end
     ϵ = method.value
     new_func = JuMP.@expression(model, [i=1:con.set.dimension], 
@@ -338,9 +374,13 @@ function reformulate_disjunct_constraint(
     method::_Hull
 ) where {T <: JuMP.GenericNonlinearExpr, S <: _MOI.Interval}
     con_func = _disaggregate_nl_expression(model, con.func, bvref, method)
-    con_func0 = JuMP.value(v -> 0.0, con.func)
-    if isinf(con_func0)
-        error("Operator `$(con.func.head)` is not defined at 0, causing the perspective function on the Hull reformulation to fail.")
+    con_func0 = evaluate_at_zero(con.func)
+    if con_func0 isa Number && isinf(con_func0)
+        error("Cannot apply Hull to the nonlinear disjunct constraint on " *
+            "`$(con.func)`: it evaluates to $(con_func0) at zero, so the " *
+            "perspective constant f(0) is not finite (for example " *
+            "`log(x)` or `1/x`). Shift the variables so the function is " *
+            "finite at zero, or use BigM or MBM.")
     end
     ϵ = method.value
     new_func = JuMP.@expression(model, ((1-ϵ)*bvref+ϵ) * con_func - ϵ*(1-bvref)*con_func0)

test/constraints/hull.jl

@@ -422,6 +422,27 @@ function test_vector_quadratic_hull_1sided(moiset)
     @test ref[1].set == moiset(2)
 end
 #less than, greater than, equalto
+function test_evaluate_at_zero()
+    model = GDPModel()
+    @variable(model, 1 <= x <= 10)
+    @variable(model, 1 <= y <= 10)
+    @test evaluate_at_zero(3) === 3.0
+    @test evaluate_at_zero(x) === 0.0
+    @test evaluate_at_zero(2x + 1.5) === 1.5
+    @test evaluate_at_zero(x * y + 2x + 1.5) === 1.5
+    @test evaluate_at_zero(exp(x) - 1) === 0.0
+    @test evaluate_at_zero(log(x + 1) * (y + 2)) === 0.0
+    @test evaluate_at_zero((x - 1)^2 / (y + 4)) === 0.25
+    branch = NonlinearExpr(:ifelse,
+        Any[NonlinearExpr(:>=, Any[x, 1]), x, y + 3])
+    @test evaluate_at_zero(branch) === 3.0
+    @test evaluate_at_zero(log(x)) == -Inf
+    @operator(model, op_f, 1, v -> v^2 + 7)
+    @test evaluate_at_zero(op_f(x)) === 7.0
+    @test_throws ErrorException evaluate_at_zero(
+        NonlinearExpr(:op_missing, Any[x]))
+end
+
 function test_scalar_nonlinear_hull_1sided_error()
     model = GDPModel()
     @variable(model, 10 <= x <= 100)
@@ -678,6 +699,7 @@ end
         test_scalar_quadratic_hull_1sided(s)
         test_scalar_nonlinear_hull_1sided(s)
     end
+    test_evaluate_at_zero()
     test_scalar_nonlinear_hull_1sided_error()
     for s in (MOI.Nonpositives, MOI.Nonnegatives, MOI.Zeros)
         test_vector_var_hull_1sided(s)

test/extensions/InfiniteDisjunctiveProgramming.jl

@@ -755,6 +755,79 @@ function test_mbm_with_derivatives()
          MOI.ALMOST_LOCALLY_SOLVED]
 end
 
+# evaluate_at_zero on each kind of InfiniteOpt reference
+function test_evaluate_at_zero_infinite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 3)
+    @variable(model, x, Infinite(t))
+    @variable(model, z)
+    @finite_parameter(model, p == 2.0)
+    @parameter_function(model, pf == t -> 2 * t)
+    @test evaluate_at_zero(x) === 0.0
+    @test evaluate_at_zero(z) === 0.0
+    @test evaluate_at_zero(deriv(x, t)) === 0.0
+    @test evaluate_at_zero(x(0.5)) === 0.0
+    @test evaluate_at_zero(p) === 2.0
+    @test evaluate_at_zero(exp(x) - 1 + p * z) === 0.0
+    @test evaluate_at_zero(log(x + p) * (z + 1)) ≈ log(2.0)
+    @test evaluate_at_zero(t) === t
+    @test evaluate_at_zero(pf) === pf
+    f0 = evaluate_at_zero(exp(x) - t)
+    @test f0 isa JuMP.GenericNonlinearExpr && parameter_refs(f0) == (t,)
+    @test_throws ErrorException evaluate_at_zero(∫(x, t))
+end
+
+# Hull on a non-quadratic nonlinear disjunct constraint of an infinite
+# variable: the perspective constant f(0) is evaluated without a
+# JuMP.value method for GeneralVariableRef. Optimum x = 0 in disjunct 1.
+function test_hull_infinite_nonlinear()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, exp(x) - 1 <= 2, Disjunct(Y[1]))
+    @constraint(model, x >= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Min, ∫(x, t))
+    @test DP.reformulate_model(model, Hull()) isa Nothing
+    juniper = JuMP.optimizer_with_attributes(
+        Juniper.Optimizer,
+        "nl_solver" => JuMP.optimizer_with_attributes(
+            Ipopt.Optimizer, "print_level" => 0),
+        "log_levels" => Symbol[])
+    set_optimizer(model, juniper)
+    @test optimize!(model, gdp_method = Hull()) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    @test objective_value(model) ≈ 0.0 atol = 1e-4
+    @test all(value(Y[1]))
+end
+
+# Hull on a nonlinear disjunct constraint that depends on the infinite
+# parameter: f(0) = -t is an expression, substituted per support by the
+# transcription. Only disjunct 1 is feasible, so x(t) = log(1 + t).
+function test_hull_infinite_nonlinear_parameter()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, exp(x) - 1 <= t, Disjunct(Y[1]))
+    @constraint(model, x <= -1, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x, t))
+    juniper = JuMP.optimizer_with_attributes(
+        Juniper.Optimizer,
+        "nl_solver" => JuMP.optimizer_with_attributes(
+            Ipopt.Optimizer, "print_level" => 0),
+        "log_levels" => Symbol[])
+    set_optimizer(model, juniper)
+    @test optimize!(model, gdp_method = Hull()) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    @test value(x) ≈ log.(1 .+ supports(t)) atol = 1e-4
+    @test all(value(Y[1]))
+end
+
 function test_CuttingPlanes_infinite_simple()
     model = InfiniteGDPModel(HiGHS.Optimizer)
     set_silent(model)
@@ -1203,6 +1276,9 @@ end
     @testset "Integration" begin
         test_infiniteopt_extension()
         test_methods()
+        test_evaluate_at_zero_infinite()
+        test_hull_infinite_nonlinear()
+        test_hull_infinite_nonlinear_parameter()
     end
 
     @testset "Cutting Planes" begin