MBM/Hull on InfiniteGDP disjunctions where constraints carry different infinite parameters (or none)

ext/InfiniteDisjunctiveProgramming.jl

@@ -170,6 +170,37 @@ function DP.disaggregate_expression(
     return JuMP.@expression(model, sum(terms))
 end
 
+# Return a copy that also contains the indicator's parameters
+function _copy_over_indicator(
+    vref::InfiniteOpt.GeneralVariableRef,
+    bvref::JuMP.AbstractJuMPScalar,
+    indicator_copies::Dict
+    )
+    # parameters and variables already over the indicator's groups stay
+    _is_parameter(vref) && return vref
+    issubset(InfiniteOpt.parameter_group_int_indices(bvref),
+        InfiniteOpt.parameter_group_int_indices(vref)) && return vref
+    # one copy per variable, shared by every row and objective it is in
+    haskey(indicator_copies, vref) && return indicator_copies[vref]
+    return indicator_copies[vref] = _copy_over_indicator(
+        InfiniteOpt.dispatch_variable_ref(vref), vref, bvref,
+        indicator_copies)
+end
+# a derivative follows its argument's copy, keeping the finite differences
+function _copy_over_indicator(
+    ::InfiniteOpt.DerivativeRef, vref, bvref, indicator_copies)
+    argument = _copy_over_indicator(InfiniteOpt.derivative_argument(vref),
+        bvref, indicator_copies)
+    return InfiniteOpt.deriv(argument, InfiniteOpt.operator_parameter(vref))
+end
+# the copy keeps the bounds and spans the union of both parameter sets
+function _copy_over_indicator(::Any, vref, bvref, indicator_copies)
+    prefs = InfiniteOpt.parameter_refs(vref + bvref)
+    properties = DP.VariableProperties(DP.get_variable_info(vref), "",
+        nothing, InfiniteOpt.Infinite(prefs...))
+    return DP.create_variable(JuMP.owner_model(vref), properties)
+end
+
 ################################################################################
 #                          MBM FOR INFINITEMODEL
 ################################################################################
@@ -188,17 +219,19 @@ function DP.copy_model_with_constraints(
         model; filter_constraints = cref -> false
         )
 
+    indicator = DP._constraint_to_indicator(model)[first(constraints)]
+    bvref = ref_map[DP.binary_variable(indicator)]
+    indicator_copies = Dict{InfiniteOpt.GeneralVariableRef,
+        InfiniteOpt.GeneralVariableRef}()
     for cref in constraints
         con = JuMP.constraint_object(cref)
         T = one(JuMP.value_type(typeof(mini)))
-        JuMP.@constraint(mini, ref_map[con.func] * T in con.set)
+        func = InfiniteOpt.map_expression.(
+            v -> _copy_over_indicator(v, bvref, indicator_copies),
+            ref_map[con.func])
+        JuMP.@constraint(mini, func * T in con.set)
     end
 
-    InfiniteOpt.build_transformation_backend!(mini)
-    transcribed = InfiniteOpt.transformation_model(mini)
-    JuMP.set_optimizer(transcribed, method.optimizer)
-    JuMP.set_silent(transcribed)
-
     # fwd_map needs every ref reachable from disjunct constraints —
     # decision vars + parameters + parameter functions so the
     # objective substitution in `prepare_max_M_objective` can look up
@@ -207,14 +240,20 @@ function DP.copy_model_with_constraints(
     fwd_map = Dict{InfiniteOpt.GeneralVariableRef,
         Vector{InfiniteOpt.GeneralVariableRef}}()
     for v in decision_vars
-        fwd_map[v] = [ref_map[v]]
+        fwd_map[v] = [_copy_over_indicator(ref_map[v], bvref,
+            indicator_copies)]
     end
     for p in InfiniteOpt.all_parameters(model)
         fwd_map[p] = [ref_map[p]]
     end
     for pf in InfiniteOpt.all_parameter_functions(model)
         fwd_map[pf] = [ref_map[pf]]
     end
+
+    InfiniteOpt.build_transformation_backend!(mini)
+    transcribed = InfiniteOpt.transformation_model(mini)
+    JuMP.set_optimizer(transcribed, method.optimizer)
+    JuMP.set_silent(transcribed)
     return DP.GDPSubmodel(mini, decision_vars, fwd_map)
 end
 
@@ -300,15 +339,18 @@ function DP.compute_M(
     method::DP._MBM
     )
     objectives = InfiniteOpt.transformation_expression(mini_expr)
+    transcribed = InfiniteOpt.transformation_model(sub.model)
+    inner_sub = DP.GDPSubmodel(transcribed, JuMP.VariableRef[],
+        Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}())
+    if objectives isa JuMP.AbstractJuMPScalar
+        return DP.compute_M(inner_sub, objectives, method)
+    end
     # transcription orders the dimensions by parameter group, which is
     # not the ascending order `parameter_refs` gives the grids below
     group_idxs = InfiniteOpt.parameter_group_int_indices(mini_expr)
     if length(group_idxs) > 1 && ndims(objectives) == length(group_idxs)
         objectives = permutedims(objectives, sortperm(group_idxs))
     end
-    transcribed = InfiniteOpt.transformation_model(sub.model)
-    inner_sub = DP.GDPSubmodel(transcribed, JuMP.VariableRef[],
-        Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}())
     prefs, grids = _support_grids(sub, mini_expr)
     M_vals = DP.sample_M_values(method.sampler, objectives,
         inner_sub, method, grids)

src/hull.jl

@@ -32,13 +32,14 @@ function _disaggregate_variable(
     lb, ub = variable_bound_info(vref)
     info = get_variable_info(vref; has_lb = true, has_ub = true, 
                              lower_bound = lb, upper_bound = ub)
+    #get binary indicator variable
+    bvref = binary_variable(lvref)
+    #the copy lives where its bound rows do, lb*bvref <= dvref <= ub*bvref
     old_props = VariableProperties(vref)
-    properties = VariableProperties(info, "$(vref)_$(lvref)", 
-                                    old_props.set, old_props.variable_type)
+    properties = VariableProperties(info, "$(vref)_$(lvref)", old_props.set,
+        VariableProperties(vref + bvref).variable_type)
     dvref = create_variable(model, properties)
     push!(_reformulation_variables(model), dvref)
-    #get binary indicator variable
-    bvref = binary_variable(lvref)
     #temp storage
     push!(method.disjunction_variables[vref], dvref)
     method.disjunct_variables[vref, bvref] = dvref

test/extensions/AbstractGPsDisjunctiveProgramming.jl

@@ -297,6 +297,45 @@ function test_gp_unknown_sampler_error()
             sampler = _UnimplementedSampler()))
 end
 
+# MBM-GP on shared variables (see the InfiniteOpt extension tests), with
+# every support solved so the result is exact
+
+# finite d under Y(t): d = min_t max(2 - 2t, 1 + 2t) = 1.5
+function test_gp_shared_finite_variable()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1],
+        supports = [0.0, 0.25, 0.5, 0.75, 1.0])
+    @variable(model, 0 <= d <= 3)
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, d <= 2 - 2t, Disjunct(Y[1]))
+    @constraint(model, d <= 1 + 2t, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, d)
+    method = MBM(HiGHS.Optimizer; sampler = GPSampler(frac_supports = 1.0))
+    optimize!(model, gdp_method = method)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+end
+
+# z(t) under Y(t, xi): z = min_xi max(2 - xi, xi) = 1 at each of two t
+function test_gp_shared_infinite_variable()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 1.0])
+    @infinite_parameter(model, xi ∈ [0, 2], supports = [0.5, 1.0, 1.5, 2.0])
+    @variable(model, 0 <= z <= 3, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t, xi))
+    @constraint(model, z <= 2 - xi, Disjunct(Y[1]))
+    @constraint(model, z <= xi, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, support_sum(z, t))
+    method = MBM(HiGHS.Optimizer; sampler = GPSampler(frac_supports = 1.0))
+    optimize!(model, gdp_method = method)
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 2.0 atol = 1e-6
+end
+
 @testset "AbstractGPsDisjunctiveProgramming" begin
     test_gp_sampler_kwargs()
     test_gp_compute_M_scalar()
@@ -309,4 +348,6 @@ end
     test_gp_detect_uniform_M_off_dependent()
     test_gp_unknown_sampler_error()
     test_gp_infeasible_disjunct()
+    test_gp_shared_finite_variable()
+    test_gp_shared_infinite_variable()
 end

test/extensions/InfiniteDisjunctiveProgramming.jl

@@ -577,6 +577,68 @@ function test_add_cut_infinite()
     @test InfiniteOpt.transformation_backend_ready(model)
 end
 
+# compute_M on a finite disjunct constraint: no supports to sample, one
+# M subproblem. Slack r(d) = 5 - d maximized over d <= 3 gives 5.
+function test_compute_M_infinite_finite_constraint()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, 0 <= d <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, con, d >= 5, Disjunct(Y[1]))
+    @constraint(model, con2, d <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(MBM(HiGHS.Optimizer), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    @test isempty(InfiniteOpt.parameter_refs(obj))
+    @test DP.compute_M(sub, obj, mbm) == 5.0
+end
+
+# MBM with a disjunction made only of finite constraints in an
+# InfiniteModel: d = 3 and x = 1 - t, objective 3 + 1/2.
+function test_mbm_finite_disjunction()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 11)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, 0 <= d <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x >= 1 - t)
+    @constraint(model, d >= 3, Disjunct(Y[1]))
+    @constraint(model, d >= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Min, d + ∫(x, t))
+    @test optimize!(model, gdp_method = MBM(HiGHS.Optimizer)) isa Nothing
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 3.5 atol = 1e-6
+    @test value(Y[1])
+end
+
+# MBM with finite and infinite constraints in the same disjunct:
+# d = 1 with x = 1 - t <= 2d, objective 1 + 1/2.
+function test_mbm_mixed_finite_disjunct()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 11)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, 0 <= d <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x >= 1 - t)
+    @constraint(model, x <= d, Disjunct(Y[1]))
+    @constraint(model, d >= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 2d, Disjunct(Y[2]))
+    @constraint(model, d >= 1, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Min, d + ∫(x, t))
+    @test optimize!(model, gdp_method = MBM(HiGHS.Optimizer)) isa Nothing
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+    @test value(Y[2])
+end
+
 # MBM with finite + integer variables in an InfiniteModel.
 function test_mbm_finite_and_integer_var()
     model = InfiniteGDPModel(HiGHS.Optimizer)
@@ -1147,6 +1209,123 @@ function test_add_cut_weighted_coefficients()
     end
 end
 
+# Shared variables: a variable whose infinite parameters do not contain
+# its indicator's, so one copy serves indicator supports that may pick
+# different disjuncts. Optima are by brute force over the disjunct
+# choices at each support.
+
+# finite d under Y(t): per t, d <= 2 - 2t or d <= 1 + 2t; maximize d.
+# Best choice per t is the larger bound, so d = min_t max(...) = 1.5.
+function test_shared_finite_variable()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1],
+        supports = [0.0, 0.25, 0.5, 0.75, 1.0])
+    @variable(model, 0 <= d <= 3)
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, d <= 2 - 2t, Disjunct(Y[1]))
+    @constraint(model, d <= 1 + 2t, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, d)
+
+    optimize!(model, gdp_method = BigM(10.0))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+
+    optimize!(model, gdp_method = MBM(HiGHS.Optimizer))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+
+    optimize!(model, gdp_method = Hull())
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+end
+
+# z(t) under Y(t, xi): per (t, xi), z <= 2 - xi or z <= xi; maximize
+# the sum over t of z, so z = min_xi max(2 - xi, xi) = 1 at each t.
+function test_shared_infinite_variable()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 1.0])
+    @infinite_parameter(model, xi ∈ [0, 2], supports = [0.5, 1.0, 1.5, 2.0])
+    @variable(model, 0 <= z <= 3, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t, xi))
+    @constraint(model, z <= 2 - xi, Disjunct(Y[1]))
+    @constraint(model, z <= xi, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, support_sum(z, t))
+
+    optimize!(model, gdp_method = BigM(10.0))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 2.0 atol = 1e-6
+
+    optimize!(model, gdp_method = MBM(HiGHS.Optimizer))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 2.0 atol = 1e-6
+
+    optimize!(model, gdp_method = Hull())
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 2.0 atol = 1e-6
+end
+
+# Y(t) over x(t, xi): no shared variable, the indicator only has fewer
+# parameters than the constraint. Per t, x <= t + xi (sum 3t + 1.5) or
+# x <= 1.2 at a cost of 1.5 (value 2.1): 2.1 + 3.0 + 4.5 = 9.6.
+function test_fewer_indicator_parameters()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @infinite_parameter(model, xi ∈ [0, 1], supports = [0.2, 0.5, 0.8])
+    @variable(model, 0 <= x <= 3, Infinite(t, xi))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x <= t + xi, Disjunct(Y[1]))
+    @constraint(model, x <= 1.2, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max,
+        support_sum(support_sum(x, xi) - 1.5 * binary_variable(Y[2]), t))
+
+    optimize!(model, gdp_method = BigM(10.0))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 9.6 atol = 1e-6
+
+    optimize!(model, gdp_method = MBM(HiGHS.Optimizer))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 9.6 atol = 1e-6
+
+    optimize!(model, gdp_method = Hull())
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 9.6 atol = 1e-6
+end
+
+# derivative of a shared z(t) under Y(t, xi): per (t, xi), dz <= 2 - xi or
+# dz <= xi, so dz = max(2 - xi, xi) = 1.5 at both xi and z(1) - z(0) = 1.5.
+function test_shared_variable_derivative()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1],
+        supports = [0.0, 0.25, 0.5, 0.75, 1.0])
+    @infinite_parameter(model, xi ∈ [0, 2], supports = [0.5, 1.5])
+    @variable(model, -5 <= z <= 5, Infinite(t))
+    @variable(model, -10 <= dz <= 10, Deriv(z, t))
+    @variable(model, Y[1:2], InfiniteLogical(t, xi))
+    @constraint(model, dz <= 2 - xi, Disjunct(Y[1]))
+    @constraint(model, dz <= xi, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, z(1) - z(0))
+
+    optimize!(model, gdp_method = BigM(100.0))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+
+    optimize!(model, gdp_method = MBM(HiGHS.Optimizer))
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+
+    optimize!(model, gdp_method = Hull())
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 1.5 atol = 1e-6
+end
+
 @testset "InfiniteDisjunctiveProgramming" begin
 
     @testset "Model" begin
@@ -1193,13 +1372,22 @@ end
         test_compute_M_infinite_two_params()
         test_compute_M_infinite_dependent_params()
         test_compute_M_infinite_dependent_varying()
+        test_compute_M_infinite_finite_constraint()
+        test_mbm_finite_disjunction()
+        test_mbm_mixed_finite_disjunct()
         test_mbm_finite_and_integer_var()
         test_mbm_infinite_simple()
         test_mbm_infinite_param_dependent()
         test_mbm_vs_bigm_infinite()
         test_mbm_with_derivatives()
     end
 
+    @testset "Shared variables" begin
+        test_shared_finite_variable()
+        test_shared_infinite_variable()
+        test_fewer_indicator_parameters()
+        test_shared_variable_derivative()
+    end
     @testset "Integration" begin
         test_infiniteopt_extension()
         test_methods()