Extensions to MBM and cutting planes techniques.

Project.toml

@@ -20,7 +20,7 @@ Reexport = "1"
 julia = "1.10"
 Juniper = "0.9.3"
 Ipopt = "1.9.0"
-InfiniteOpt = "0.6"
+InfiniteOpt = "0.6.3"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

ext/InfiniteDisjunctiveProgramming.jl

@@ -9,8 +9,8 @@ import DisjunctiveProgramming as DP
 ################################################################################
 function DP.InfiniteGDPModel(args...; kwargs...)
     return DP.GDPModel{
-        InfiniteOpt.InfiniteModel, 
-        InfiniteOpt.GeneralVariableRef, 
+        InfiniteOpt.InfiniteModel,
+        InfiniteOpt.GeneralVariableRef,
         InfiniteOpt.InfOptConstraintRef
         }(args...; kwargs...)
 end
@@ -26,7 +26,7 @@ end
 ################################################################################
 DP.InfiniteLogical(prefs...) = DP.Logical(InfiniteOpt.Infinite(prefs...))
 
-_is_parameter(vref::InfiniteOpt.GeneralVariableRef) = 
+_is_parameter(vref::InfiniteOpt.GeneralVariableRef) =
     _is_parameter(InfiniteOpt.dispatch_variable_ref(vref))
 _is_parameter(::InfiniteOpt.DependentParameterRef) = true
 _is_parameter(::InfiniteOpt.IndependentParameterRef) = true
@@ -171,17 +171,207 @@ function DP.disaggregate_expression(
 end
 
 ################################################################################
-#                              ERROR MESSAGES
+#                          MBM FOR INFINITEMODEL
+################################################################################
+
+# Copy the InfiniteModel, strip everything but VariableInfo bounds,
+# add back the selected disjunct constraints, transcribe, and return
+# only the other disjunct's constraints plus variable bounds.
+function DP.copy_model_with_constraints(
+    model::InfiniteOpt.InfiniteModel,
+    constraints::Vector{<:DP.DisjunctConstraintRef},
+    method::DP._MBM
+    )
+    # Filter out every source constraint at copy time instead of
+    # copying then deleting. Equivalent end state, fewer allocations.
+    mini, ref_map = JuMP.copy_model(
+        model; filter_constraints = cref -> false
+        )
+
+    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)
+    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
+    # any term it sees.
+    decision_vars = DP.collect_all_vars(model)
+    fwd_map = Dict{InfiniteOpt.GeneralVariableRef,
+        Vector{InfiniteOpt.GeneralVariableRef}}()
+    for v in decision_vars
+        fwd_map[v] = [ref_map[v]]
+    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
+    return DP.GDPSubmodel(mini, decision_vars, fwd_map)
+end
+
+function DP.prepare_max_M_objective(
+    ::InfiniteOpt.InfiniteModel,
+    obj::JuMP.ScalarConstraint{T, S},
+    sub::DP.GDPSubmodel
+    ) where {T, S <: _MOI.LessThan}
+    flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map)
+    obj_func = DP.replace_variables_in_constraint(obj.func, flat_map)
+    return obj_func - obj.set.upper
+end
+
+function DP.prepare_max_M_objective(
+    ::InfiniteOpt.InfiniteModel,
+    obj::JuMP.ScalarConstraint{T, S},
+    sub::DP.GDPSubmodel
+    ) where {T, S <: _MOI.GreaterThan}
+    flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map)
+    obj_func = DP.replace_variables_in_constraint(obj.func, flat_map)
+    return obj.set.lower - obj_func
+end
+
+# Constant interpolation
+function _interpolate(
+    grids::NTuple{N, AbstractVector{<:Real}},
+    values::AbstractArray{<:Real, N}
+    ) where {N}
+    # mimic the call form of Interpolations.jl's interpolation
+    return (args...) -> _interpolate_at(grids, values, args)
+end
+
+function _interpolate_at(
+    grids::NTuple{N, AbstractVector{<:Real}},
+    values::AbstractArray{<:Real, N},
+    args::NTuple{N, <:Real}
+    ) where {N}
+    # lower-corner cell index per dimension
+    idx_lo = ntuple(d -> 
+        clamp(searchsortedlast(grids[d], args[d]),1, length(grids[d]) - 1), N
+    )
+    # max over the 2^N corners; bit d of k picks lower or upper
+    return maximum(
+        values[ntuple(d -> idx_lo[d] +((k >> (d - 1)) & 1), N)...]
+        for k in 0:(2^N - 1)
+        )
+end
+
+# Transcribe mini_expr, solve per support on the transcribed JuMP
+# model, and aggregate to a scalar if uniform, else to a parameter
+# function on main.
+function DP.raw_M(
+    sub::DP.GDPSubmodel{<:InfiniteOpt.InfiniteModel},
+    mini_expr::JuMP.AbstractJuMPScalar,
+    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}}()
+        )
+    M_vals = Array{typeof(method.default_M)}(undef, size(objectives))
+    for I in eachindex(objectives)
+        m = DP.raw_M(inner_sub, objectives[I], method)
+        m === nothing && return nothing
+        M_vals[I] = m
+    end
+    all(==(first(M_vals)), M_vals) && return first(M_vals)
+    mini_prefs = InfiniteOpt.parameter_refs(mini_expr)
+    reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map)
+    prefs = Tuple(reverse_map[p] for p in mini_prefs)
+    main = JuMP.owner_model(first(prefs))
+    grids = Tuple(InfiniteOpt.supports(p) for p in prefs)
+    param_func = InfiniteOpt.build_parameter_function(
+        error, _interpolate(grids, M_vals), prefs)
+    return InfiniteOpt.add_parameter_function(main, param_func)
+end
+
+################################################################################
+#                    CUTTING PLANES FOR INFINITEMODEL
 ################################################################################
-function DP.reformulate_model(::InfiniteOpt.InfiniteModel, ::DP.MBM)
-    error("The `MBM` method is not supported for `InfiniteModel`." *
-          "Please use `BigM`, `Hull`, `Indicator`, or `PSplit` instead.")
+
+# Build CP subproblem: reformulate the InfiniteModel in-place, transcribe,
+# copy, and wrap in GDPSubmodel with forward variable map.
+function DP.copy_and_reformulate(
+    model::InfiniteOpt.InfiniteModel,
+    decision_vars::Vector{InfiniteOpt.GeneralVariableRef},
+    reform_method::DP.AbstractReformulationMethod,
+    method::DP.CuttingPlanes
+    )
+    DP.reformulate_model(model, reform_method)
+    InfiniteOpt.build_transformation_backend!(model)
+    transcribed = InfiniteOpt.transformation_model(model)
+    transcription_fwd = Dict{InfiniteOpt.GeneralVariableRef,
+        Vector{JuMP.VariableRef}}()
+    for v in DP.collect_all_vars(model)
+        transcription_var = InfiniteOpt.transformation_variable(v)
+        var_prefs = InfiniteOpt.parameter_refs(v)
+        transcription_fwd[v] = isempty(var_prefs) ?
+            [transcription_var] : vec(transcription_var)
+    end
+    sub_copy, copy_map = JuMP.copy_model(transcribed)
+    fwd_map = Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.VariableRef}}()
+    for v in decision_vars
+        haskey(transcription_fwd, v) || continue
+        fwd_map[v] = [copy_map[transcribed_var] for transcribed_var in transcription_fwd[v]]
+    end
+    sub = DP.GDPSubmodel(sub_copy, decision_vars, fwd_map)
+    JuMP.set_optimizer(sub.model, method.optimizer)
+    JuMP.set_silent(sub.model)
+    return sub
 end
 
-function DP.reformulate_model(::InfiniteOpt.InfiniteModel, ::DP.CuttingPlanes)
-    error("The `CuttingPlanes` method is not supported for `InfiniteModel`." *
-          "Please use `BigM`, `Hull`, `Indicator`, or `PSplit` instead.")
+# Read per-support values from the transformation backend.
+function DP.extract_solution(model::InfiniteOpt.InfiniteModel)
+    dvars = DP.collect_cutting_planes_vars(model)
+    V = eltype(dvars)
+    T = JuMP.value_type(typeof(model))
+    sol = Dict{V, Vector{T}}()
+    for v in dvars
+        transcription_var = InfiniteOpt.transformation_variable(v)
+        var_prefs = InfiniteOpt.parameter_refs(v)
+        sol[v] = isempty(var_prefs) ? [JuMP.value(transcription_var)] :
+            JuMP.value.(vec(transcription_var))
+    end
+    return sol
 end
 
+# Add a pointwise-sum cut directly to the transformation backend and mark
+# it ready so the next optimize! doesn't re-transcribe and wipe the cut.
+function DP.add_cut(
+    model::InfiniteOpt.InfiniteModel,
+    decision_vars::Vector{InfiniteOpt.GeneralVariableRef},
+    rBM_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}},
+    sep_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}}
+    )
+    transcribed = InfiniteOpt.transformation_model(model)
+    cut_expr = zero(JuMP.GenericAffExpr{
+        JuMP.value_type(typeof(transcribed)),
+        JuMP.variable_ref_type(transcribed)})
+    for var in decision_vars
+        haskey(rBM_sol, var) || continue
+        haskey(sep_sol, var) || continue
+        rbm_vals = rBM_sol[var]
+        sep_vals = sep_sol[var]
+        transcription_var = InfiniteOpt.transformation_variable(var)
+        transcribed_vars = transcription_var isa AbstractArray ?
+            vec(transcription_var) : [transcription_var]
+        for k in eachindex(transcribed_vars)
+            xi = 2 * (sep_vals[k] - rbm_vals[k])
+            JuMP.add_to_expression!(cut_expr, xi, transcribed_vars[k])
+            JuMP.add_to_expression!(cut_expr, -xi * sep_vals[k])
+        end
+    end
+    JuMP.@constraint(transcribed, cut_expr >= 0)
+    InfiniteOpt.set_transformation_backend_ready(model, true)
+    return
 end
 
+end