diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 2b9863c..7ed9e29 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/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) diff --git a/src/hull.jl b/src/hull.jl index 779369e..8ec0949 100644 --- a/src/hull.jl +++ b/src/hull.jl @@ -281,6 +281,33 @@ 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}, @@ -288,9 +315,13 @@ function reformulate_disjunct_constraint( 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) diff --git a/test/constraints/hull.jl b/test/constraints/hull.jl index d01b30e..f599deb 100644 --- a/test/constraints/hull.jl +++ b/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) diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl index 36ddefe..749c471 100644 --- a/test/extensions/InfiniteDisjunctiveProgramming.jl +++ b/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