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,53 @@ 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
+    # no smaller than the exact M = max(4t - 1, 0) and max(1 - 4t, 0)
+    M_funcs = InfiniteOpt.all_parameter_functions(model)
+    @test all(value(M_funcs[1]) .>= [0, 0, 1, 2, 3] .- 1e-6)
+    @test all(value(M_funcs[2]) .>= [1, 0, 0, 0, 0] .- 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
+    # no smaller than the exact M = max(2xi - 2, 0) and max(2 - 2xi, 0)
+    M_funcs = InfiniteOpt.all_parameter_functions(model)
+    @test all(value(M_funcs[1]) .>= [0 0 1 2; 0 0 1 2] .- 1e-6)
+    @test all(value(M_funcs[2]) .>= [1 0 0 0; 1 0 0 0] .- 1e-6)
+end
+
 @testset "AbstractGPsDisjunctiveProgramming" begin
     test_gp_sampler_kwargs()
     test_gp_compute_M_scalar()
@@ -309,4 +356,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