From ae550d1d87a94c985645e6759cbace832bfc3970 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Sun, 22 Mar 2026 20:39:03 -0400 Subject: [PATCH 01/13] inital split --- src/datatypes.jl | 60 ++- src/mbm.jl | 732 ++++++++++++++++++++---------------- src/utilities.jl | 2 +- src/variables.jl | 3 +- test/constraints/mbm.jl | 804 +++++++++++++++++++++++++++++++++------- 5 files changed, 1131 insertions(+), 470 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index 6ba13159..f37c8e2d 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -388,15 +388,22 @@ end mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMethod optimizer::O - M::Dict{LogicalVariableRef{M}, T} - default_M::T - conlvref::Vector{LogicalVariableRef{M}} + M::Dict{LogicalVariableRef{M}, Any} + default_M::T + conlvref::Vector{LogicalVariableRef{M}} + deactivated::Set{LogicalVariableRef{M}} + # Stored submodels: indicator => GDPSubmodel. + # Typed Any so extensions can store different types. + store::Dict{LogicalVariableRef{M}, Any} function _MBM(method::MBM{O, T}, model::M) where {O, T, M <: JuMP.AbstractModel} - new{O, T, M}(method.optimizer, - Dict{LogicalVariableRef{M}, T}(), + new{O, T, M}( + method.optimizer, + Dict{LogicalVariableRef{M}, Any}(), method.default_M, - Vector{LogicalVariableRef{M}}() + Vector{LogicalVariableRef{M}}(), + Set{LogicalVariableRef{M}}(), + Dict{LogicalVariableRef{M}, Any}() ) end end @@ -444,23 +451,38 @@ A type for using the cutting planes approach for disjunctive constraints. method to use after cutting planes (default = `BigM()`). - `M_value::Float64`: Big-M value to use in the final reformulation (default = `1e9`). """ -struct cutting_planes{O} <: AbstractReformulationMethod +struct cutting_planes{O, T} <: AbstractReformulationMethod optimizer::O; max_iter::Int - seperation_tolerance::Float64 + seperation_tolerance::T final_reform_method::AbstractReformulationMethod - M_value::Float64 + M_value::T function cutting_planes( - optimizer::O; - max_iter::Int = 3, - seperation_tolerance::Float64 = 1e-6, - final_reform_method = BigM(), - M_value::Float64 = 1e9 - ) where {O} - new{O}(optimizer, max_iter, seperation_tolerance, final_reform_method, M_value) + optimizer::O; + max_iter::Int = 3, + seperation_tolerance::T = 1e-6, + final_reform_method = BigM(), + M_value::T = 1e9 + ) where {O, T} + new{O, T}(optimizer, max_iter, seperation_tolerance, final_reform_method, M_value) end end +################################################################################ +# GDP SUBMODEL +################################################################################ + +# Unified submodel wrapper for MBM and cutting planes. +# Holds a flat JuMP model, ordered decision variables, +# and a forward map (orig var → submodel vars). +struct GDPSubmodel{M <: JuMP.AbstractModel, + V <: JuMP.AbstractVariableRef, + W <: JuMP.AbstractVariableRef} + model::M + dec_vars::Vector{V} + fwd::Dict{V, Vector{W}} +end + """ PSplit <: AbstractReformulationMethod @@ -635,9 +657,9 @@ end """ VariableProperties(expr)::VariableProperties -Creates a `VariableProperties` object with blank variable info (no bounds, not fixed, -not binary/integer) from an expression. The `expr` argument is provided for -extensions to infer additional properties (e.g., parameter dependencies in InfiniteOpt). +Creates a `VariableProperties` object with blank variable info (no bounds, not fixed, +not binary/integer) from an expression. The `expr` argument is provided for +extensions to infer additional properties. ## Arguments - `expr`: Expression for extensions to extract metadata from diff --git a/src/mbm.jl b/src/mbm.jl index e9b1f133..9d731740 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -1,140 +1,178 @@ +################################################################################ +# HELPER FUNCTIONS +################################################################################ +# Check if M result contains only zeros. +_is_all_zeros(M::Number) = iszero(M) +_is_all_zeros(M::AbstractVector) = all(iszero, M) +_is_all_zeros(::Any) = false + ################################################################################ # CONSTRAINT, DISJUNCTION, DISJUNCT REFORMULATION ################################################################################ -#Reformulates the disjunction using multiple big-M values +# Reformulates the disjunction using multiple big-M values per constraint function reformulate_disjunction( - model::JuMP.AbstractModel, - disj::Disjunction, + model::JuMP.AbstractModel, + disj::Disjunction, method::MBM -) + ) mbm = _MBM(method, model) - ref_cons = Vector{JuMP.AbstractConstraint}() + disjunct_cons = Dict{LogicalVariableRef, Vector{JuMP.AbstractConstraint}}() + for d in disj.indicators + d in mbm.deactivated && continue + mbm.conlvref = filter( + x -> x != d && !(x in mbm.deactivated), disj.indicators) + disjunct_cons[d] = Vector{JuMP.AbstractConstraint}() + _reformulate_disjunct(model, disjunct_cons[d], d, mbm) + end + # Collect constraints from non-deactivated disjuncts. It needs to be + # in a separate loop because disjuncts are only deactivated by looking + # at reforming other disjuncts (subproblem infeasibility). + ref_cons = Vector{JuMP.AbstractConstraint}() for d in disj.indicators - mbm.conlvref = filter(x -> x != d, disj.indicators) - _reformulate_disjunct(model, ref_cons, d, mbm) + d in mbm.deactivated && continue + haskey(disjunct_cons, d) && append!(ref_cons, disjunct_cons[d]) end return ref_cons end -#Reformualates a disjunct the disjunct of interest -#represented by lvref and the other indicators in conlvref + +# Reformulates a disjunct represented by lvref using per-constraint M values. +# Gets its own set of M_{ie,i'} values for each other disjunct term i'. function _reformulate_disjunct( - model::JuMP.AbstractModel, - ref_cons::Vector{JuMP.AbstractConstraint}, - lvref::LogicalVariableRef, - method::_MBM -) - - empty!(method.M) + model::JuMP.AbstractModel, + ref_cons::Vector{JuMP.AbstractConstraint}, + lvref::LogicalVariableRef, method::_MBM + ) !haskey(_indicator_to_constraints(model), lvref) && return - bconref = Dict(d => binary_variable(d) for d in method.conlvref) - + # Filter out deactivated disjuncts from binary variable mapping in + # the event we've identified some infeasible disjuncts already + active_conlvref = filter(d -> !(d in method.deactivated), method.conlvref) + bconref = Dict(d => binary_variable(d) for d in active_conlvref) + constraints = _indicator_to_constraints(model)[lvref] - filtered_constraints = [c for c in constraints if c isa DisjunctConstraintRef] - - for d in method.conlvref - d_constraints = _indicator_to_constraints(model)[d] - disjunct_constraints = [c for c in d_constraints if c isa DisjunctConstraintRef] - if !isempty(disjunct_constraints) - method.M[d] = maximum( - _maximize_M( - model, - JuMP.constraint_object(cref), - disjunct_constraints, - method - ) for cref in filtered_constraints - ) + filtered_constraints = [ + c for c in constraints if c isa DisjunctConstraintRef] + + # For each constraint, compute its own set of M values + for cref in filtered_constraints + empty!(method.M) + + for d in method.conlvref + # Skip already-deactivated disjuncts + d in method.deactivated && continue + + d_constraints = _indicator_to_constraints(model)[d] + disjunct_constraints = [ + c for c in d_constraints if c isa DisjunctConstraintRef] + if !isempty(disjunct_constraints) + M_result = _maximize_M(model, JuMP.constraint_object(cref), + disjunct_constraints, method) + # Check for infeasibility: disjunct d + # has empty feasible region + if M_result === nothing + push!(method.deactivated, d) + @warn "Disjunct $(d) is infeasible, deactivating." + delete!(bconref, d) + else + method.M[d] = M_result + end + end + end + + con = JuMP.constraint_object(cref) + # Check if all M values are zero for that constraint. If so, it + # should be enforced globally (no reformulation with binaries). + if !isempty(method.M) && all( + _is_all_zeros(method.M[d]) for d in keys(method.M)) + push!(ref_cons, con) + else + append!(ref_cons, + reformulate_disjunct_constraint(model, con, bconref, method)) end - end - for cref in filtered_constraints - con = JuMP.constraint_object(cref) - append!(ref_cons, reformulate_disjunct_constraint(model, con, - bconref, method)) end return ref_cons end function reformulate_disjunct_constraint( - model::JuMP.AbstractModel, - con::Disjunction, - bconref::Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + model::JuMP.AbstractModel, con::Disjunction, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) - + ) ref_cons = reformulate_disjunction(model, con, MBM(method.optimizer)) new_ref_cons = Vector{JuMP.AbstractConstraint}() for ref_con in ref_cons - append!(new_ref_cons, - reformulate_disjunct_constraint(model, ref_con, bconref, method) - ) + append!(new_ref_cons, + reformulate_disjunct_constraint(model, ref_con, bconref, method)) end return new_ref_cons end +# Uses per-row M values: method.M[d][row] for each disjunct d function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.VectorConstraint{T, S, R}, - bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) where {T, S <: _MOI.Nonpositives, R} - m_sum = sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) where {T, S <: _MOI.Nonpositives, R} new_func = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] - m_sum - ) + con.func[i] - sum(method.M[d][i] * bconref[d] for d in keys(method.M))) reform_con = JuMP.build_constraint(error, new_func, con.set) return [reform_con] end - +# Uses per-row M values: method.M[d][row] for each disjunct d function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.VectorConstraint{T, S, R}, - bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) where {T, S <: _MOI.Nonnegatives, R} - m_sum = sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) where {T, S <: _MOI.Nonnegatives, R} new_func = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] + m_sum - ) + con.func[i] + sum(method.M[d][i] * bconref[d] for d in keys(method.M))) reform_con = JuMP.build_constraint(error, new_func, con.set) return [reform_con] end +# Uses per-row M values: method.M[d][row] for each disjunct d function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.VectorConstraint{T, S, R}, - bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) where {T, S <: _MOI.Zeros, R} - m_sum = sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) where {T, S <: _MOI.Zeros, R} upper_expr = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] + m_sum - ) + con.func[i] + sum(method.M[d][i] * bconref[d] for d in keys(method.M))) lower_expr = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] - m_sum - ) - upper_con = JuMP.build_constraint(error, upper_expr, - MOI.Nonnegatives(con.set.dimension) - ) - lower_con = JuMP.build_constraint(error, lower_expr, - MOI.Nonpositives(con.set.dimension) - ) + con.func[i] - sum(method.M[d][i] * bconref[d] for d in keys(method.M))) + upper_con = JuMP.build_constraint( + error, upper_expr, MOI.Nonnegatives(con.set.dimension)) + lower_con = JuMP.build_constraint( + error, lower_expr, MOI.Nonpositives(con.set.dimension)) return [upper_con, lower_con] end - function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.ScalarConstraint{T, S}, - bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) where {T, S <: _MOI.LessThan} - new_func = JuMP.@expression(model, - con.func - sum(method.M[i] * bconref[i] for i in keys(method.M))) + ) where {T, S <: _MOI.LessThan} + new_func = JuMP.@expression(model, con.func - sum( + method.M[i] * bconref[i] for i in keys(method.M))) reform_con = JuMP.build_constraint(error, new_func, con.set) return [reform_con] end @@ -142,336 +180,386 @@ end function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.ScalarConstraint{T, S}, - bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) where {T, S <: _MOI.GreaterThan} - new_func = JuMP.@expression(model, - con.func + sum(method.M[i] * bconref[i] for i in keys(method.M)) - ) + ) where {T, S <: _MOI.GreaterThan} + new_func = JuMP.@expression(model, con.func + sum( + method.M[i] * bconref[i] for i in keys(method.M))) reform_con = JuMP.build_constraint(error, new_func, con.set) return [reform_con] end +# Uses per-bound M values: method.M[d][1] for lower, method.M[d][2] for upper function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.ScalarConstraint{T, S}, - bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) where {T, S <: _MOI.EqualTo} - upper_func = JuMP.@expression(model, - con.func - sum(method.M[i] * bconref[i] for i in keys(method.M)) - ) - lower_func = JuMP.@expression(model, - con.func + sum(method.M[i] * bconref[i] for i in keys(method.M)) - ) - upper_con = JuMP.build_constraint(error, upper_func, - MOI.LessThan(con.set.value) - ) - lower_con = JuMP.build_constraint(error, lower_func, - MOI.GreaterThan(con.set.value) - ) + ) where {T, S <: _MOI.EqualTo} + # M[d][1] = M for GreaterThan (lower), M[d][2] = M for LessThan (upper) + lower_func = JuMP.@expression(model, con.func + sum( + method.M[d][1] * bconref[d] for d in keys(method.M))) + upper_func = JuMP.@expression(model, con.func - sum( + method.M[d][2] * bconref[d] for d in keys(method.M))) + lower_con = JuMP.build_constraint( + error, lower_func, MOI.GreaterThan(con.set.value)) + upper_con = JuMP.build_constraint( + error, upper_func, MOI.LessThan(con.set.value)) return [lower_con, upper_con] end +# Uses per-bound M values: method.M[d][1] for lower, method.M[d][2] for upper function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.ScalarConstraint{T, S}, - bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + bconref::Union{ + Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} + }, method::_MBM -) where {T, S <: _MOI.Interval} - set_values = _set_values(con.set) - upper_func = JuMP.@expression(model, - con.func - sum(method.M[i] * bconref[i] for i in keys(method.M)) - ) - upper_con = JuMP.build_constraint(error, upper_func, - MOI.LessThan(set_values[2]) - ) - - lower_func = JuMP.@expression(model, - con.func + sum(method.M[i] * bconref[i] for i in keys(method.M)) - ) - lower_con = JuMP.build_constraint(error, lower_func, - MOI.GreaterThan(set_values[1]) - ) - + ) where {T, S <: _MOI.Interval} + set_values = _set_values(con.set) + # M[d][1] = M for GreaterThan (lower), M[d][2] = M for LessThan (upper) + lower_func = JuMP.@expression(model, con.func + sum( + method.M[d][1] * bconref[d] for d in keys(method.M))) + upper_func = JuMP.@expression(model, con.func - sum( + method.M[d][2] * bconref[d] for d in keys(method.M))) + lower_con = JuMP.build_constraint( + error, lower_func, MOI.GreaterThan(set_values[1])) + upper_con = JuMP.build_constraint( + error, upper_func, MOI.LessThan(set_values[2])) return [lower_con, upper_con] end -function reformulate_disjunct_constraint( - ::JuMP.AbstractModel, - ::F, - ::Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, - Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, +function reformulate_disjunct_constraint(::JuMP.AbstractModel, + ::F, + ::Union{Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr}}, ::_MBM -) where {F} - error("Constraint type $(typeof(F)) is not supported by the " * - "Multiple Big-M reformulation method.") + ) where {F} + error("Constraint type $(typeof(F)) is not supported by " * + "the Multiple Big-M reformulation method.") end ################################################################################ -# MULTIPLE BIG-M REFORMULATION +# MULTIPLE BIG-M REFORMULATION ################################################################################ -# Dispatches over constraint types to reformulate into >= or <= -# in order to solve the mini-model -function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + +# Prepare flat objectives for _raw_M. Returns a vector of objective expressions +# ready to maximize. Base: single flat constraint via fwd[v][1]. +function _prepare_objectives( + ::JuMP.AbstractModel, + obj::JuMP.ScalarConstraint{T, S}, + sub::GDPSubmodel + ) where {T, S <: _MOI.LessThan} + flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd) + expr = -obj.set.upper +_replace_variables_in_constraint(obj.func, flat_map) + return [expr] +end + +function _prepare_objectives( + ::JuMP.AbstractModel, + obj::JuMP.ScalarConstraint{T, S}, + sub::GDPSubmodel + ) where {T, S <: _MOI.GreaterThan} + flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd) + expr = obj.set.lower -_replace_variables_in_constraint(obj.func, flat_map) + return [expr] +end + +# Solve the submodel for each objective and return raw M values as a vector, or +# nothing if infeasible. Order of results matches order of objectives. +function _raw_M( + sub::GDPSubmodel, + objectives::Vector{<:JuMP.AbstractJuMPScalar}, method::_MBM -) where { T, S <: _MOI.Nonpositives, R} - val_type = JuMP.value_type(typeof(model)) - return maximum( - _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), - constraints, - method - ) for i in 1:objective.set.dimension ) + M_vals = typeof(method.default_M)[] + for obj_expr in objectives + JuMP.@objective(sub.model, Max, obj_expr) + JuMP.optimize!(sub.model) + if JuMP.termination_status(sub.model) == _MOI.INFEASIBLE + return nothing + elseif !JuMP.is_solved_and_feasible(sub.model) + push!(M_vals, method.default_M) + # Clear NaN start values from non-feasible solve + # so the next objective doesn't inherit them + for v in JuMP.all_variables(sub.model) + JuMP.set_start_value(v, nothing) + end + else + push!(M_vals, max( + JuMP.objective_value(sub.model), + zero(method.default_M)) + ) + end + end + return M_vals end +# Condense flat per-support values to final form. Base: return +# scalar from single-element vector. Extensions may override. +function condense_values( + ::JuMP.AbstractModel, + vals::AbstractVector + ) + return vals[1] +end + +# Dispatch over constraint types to compute M values. Scalar +# LE/GE: prepare objectives, solve, finalize. function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.ScalarConstraint{T, S}, + constraints::Vector{<:DisjunctConstraintRef}, + method::_MBM + ) where {T, S <: Union{_MOI.LessThan, _MOI.GreaterThan}} + sub = _get_submodel(model, constraints, method) + objectives = _prepare_objectives(model, objective, sub) + raw = _raw_M(sub, objectives, method) + raw === nothing && return nothing + return condense_values(model, raw) +end + +# Helper: get or create the submodel for a set of constraints. +function _get_submodel( + model::JuMP.AbstractModel, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where { T, S <: _MOI.Nonnegatives, R} - val_type = JuMP.value_type(typeof(model)) - return maximum( - _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(zero(val_type))), - constraints, - method - ) for i in 1:objective.set.dimension ) + indicator = _constraint_to_indicator( + model)[first(constraints)] + if !haskey(method.store, indicator) + method.store[indicator] = create_submodel( + model, constraints, method) + end + return method.store[indicator] end +# EqualTo: solve both GreaterThan and LessThan directions, finalize each. function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.ScalarConstraint{T, S}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where { T, S <: _MOI.Zeros, R} - val_type = JuMP.value_type(typeof(model)) - return max( - maximum( - _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i],MOI.GreaterThan(zero(val_type))), - constraints, - method - ) for i in 1:objective.set.dimension - ), - maximum( - _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), - constraints, - method - ) for i in 1:objective.set.dimension - ) - ) + ) where {T, S <: _MOI.EqualTo} + sub = _get_submodel(model, constraints, method) + set_value = objective.set.value + ge_obj = JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)) + le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)) + raw_lower = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) + raw_upper = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + (raw_lower === nothing || raw_upper === nothing) && + return nothing + return [condense_values(model, raw_lower),condense_values(model, raw_upper)] end +# Interval: solve both lower and upper bound directions, finalize each. function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.ScalarConstraint{T, S}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.ScalarConstraint{T, S}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where {T, S <: Union{_MOI.LessThan, _MOI.GreaterThan}} - return _mini_model(model, objective, constraints, method) + ) where {T, S <: _MOI.Interval} + sub = _get_submodel(model, constraints, method) + set_values = _set_values(objective.set) + ge_obj = JuMP.ScalarConstraint(objective.func, + MOI.GreaterThan(set_values[1])) + le_obj = JuMP.ScalarConstraint(objective.func, + MOI.LessThan(set_values[2])) + raw_lower = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) + raw_upper = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + (raw_lower === nothing || raw_upper === nothing) && + return nothing + return [condense_values(model, raw_lower),condense_values(model, raw_upper)] end +# Nonpositives: per-row LessThan solves for each dimension of the vector. function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.ScalarConstraint{T, S}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where {T, S <: _MOI.EqualTo} - set_value = objective.set.value - return max( - _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)), - constraints, - method - ), - _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)), - constraints, - method - ) - ) + ) where {T, S <: _MOI.Nonpositives, R} + sub = _get_submodel(model, constraints, method) + val_type = JuMP.value_type(typeof(model)) + results = Any[] + for i in 1:objective.set.dimension + le_obj = JuMP.ScalarConstraint( + objective.func[i], MOI.LessThan(zero(val_type))) + raw = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + raw === nothing && return nothing + push!(results, condense_values(model, raw)) + end + return results end +# Nonnegatives: per-row GreaterThan solves for each dimension of the vector. function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.ScalarConstraint{T, S}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where {T, S <: _MOI.Interval} - set_values = _set_values(objective.set) # Returns (lower, upper) - return max( - _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_values[1])), - constraints, - method - ), - _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_values[2])), - constraints, - method - ) - ) + ) where {T, S <: _MOI.Nonnegatives, R} + sub = _get_submodel(model, constraints, method) + val_type = JuMP.value_type(typeof(model)) + results = Any[] + for i in 1:objective.set.dimension + ge_obj = JuMP.ScalarConstraint( + objective.func[i], MOI.GreaterThan(zero(val_type))) + raw = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) + raw === nothing && return nothing + push!(results, condense_values(model, raw)) + end + return results end +# Zeros: per-row element-wise max of GE and LE raw values, then finalize. function _maximize_M( - ::JuMP.AbstractModel, - ::F, - ::Vector{<:DisjunctConstraintRef}, - ::_MBM -) where {F} - error("This type of constraints and objective constraint has " * - "not been implemented for MBM subproblems\nF: $(F)") + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, + method::_MBM + ) where {T, S <: _MOI.Zeros, R} + sub = _get_submodel(model, constraints, method) + val_type = JuMP.value_type(typeof(model)) + results = Any[] + for i in 1:objective.set.dimension + ge_obj = JuMP.ScalarConstraint( + objective.func[i], MOI.GreaterThan(zero(val_type))) + le_obj = JuMP.ScalarConstraint( + objective.func[i], MOI.LessThan(zero(val_type))) + raw_ge = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) + raw_le = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + (raw_ge === nothing || raw_le === nothing) && + return nothing + push!(results, condense_values(model, max.(raw_ge, raw_le))) + end + return results +end + +function _maximize_M( + ::JuMP.AbstractModel, ::F, + ::Vector{<:DisjunctConstraintRef}, ::_MBM + ) where {F} + error("This type of constraints and objective constraint " * + "has not been implemented for MBM subproblems\nF: $(F)") end -# Solve a mini-model to find the maximum value of the objective -# function for M value -function _mini_model( - model::JuMP.AbstractModel, - objective::JuMP.ScalarConstraint{T,S}, - constraints::Vector{<:DisjunctConstraintRef}, +# Create a submodel for a disjunct's feasible region. Returns +# GDPSubmodel. Extensions may override for custom construction. +function create_submodel( + model::JuMP.AbstractModel, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where {T,S <: Union{_MOI.LessThan, _MOI.GreaterThan}} + ) var_type = JuMP.variable_ref_type(model) sub_model = _copy_model(model) - new_vars = Dict{var_type, var_type}() - for var in collect_all_vars(model) - new_vars[var] = variable_copy(sub_model, var) + dec_vars = collect_all_vars(model) + fwd = Dict{var_type, Vector{var_type}}() + + for var in dec_vars + copy_var = variable_copy(sub_model, var) + fwd[var] = [copy_var] end - for con in [JuMP.constraint_object(con) for con in constraints] - expr = _replace_variables_in_constraint(con.func, new_vars) - JuMP.@constraint(sub_model, expr * 1.0 in con.set) + + for cref in constraints + con = JuMP.constraint_object(cref) + flat_map = Dict(v => ws[1] for (v, ws) in fwd) + expr = _replace_variables_in_constraint( + con.func, flat_map) + T = one(JuMP.value_type(typeof(sub_model))) + JuMP.@constraint(sub_model, expr * T in con.set) end - _constraint_to_objective(sub_model, objective, new_vars) + JuMP.set_optimizer(sub_model, method.optimizer) JuMP.set_silent(sub_model) - JuMP.optimize!(sub_model) - if JuMP.termination_status(sub_model) != MOI.OPTIMAL || - !JuMP.has_values(sub_model) || - JuMP.primal_status(sub_model) != MOI.FEASIBLE_POINT - M = method.default_M - else - M = JuMP.objective_value(sub_model) - end - return M -end - -################################################################################ -# CONSTRAINT TO OBJECTIVE -################################################################################ -function _constraint_to_objective( - sub_model::JuMP.AbstractModel, - obj::JuMP.ScalarConstraint{<:JuMP.AbstractJuMPScalar, MOI.LessThan{T}}, - new_vars::Dict{V,K} -) where {T,V <: JuMP.AbstractVariableRef, K <: JuMP.AbstractVariableRef} - JuMP.@objective(sub_model, Max, - - obj.set.upper + _replace_variables_in_constraint(obj.func, new_vars) - ) -end -function _constraint_to_objective( - sub_model::JuMP.AbstractModel, - obj::JuMP.ScalarConstraint{<:JuMP.AbstractJuMPScalar, MOI.GreaterThan{T}}, - new_vars::Dict{V,K} -) where {T,V <: JuMP.AbstractVariableRef, K <: JuMP.AbstractVariableRef} - JuMP.@objective(sub_model, Max, - - _replace_variables_in_constraint(obj.func, new_vars) + obj.set.lower - ) -end -function _constraint_to_objective( - sub_model::JuMP.AbstractModel, - obj::JuMP.ScalarConstraint, - new_vars::Dict{V,K} -) where {V, K} - error("This type of constraint is not supported, only greater " * - "than and less than constraints are supported with " * - "intervals and equalities being converted.") + return GDPSubmodel(sub_model, dec_vars, fwd) end ################################################################################ -# REPLACE VARIABLES IN CONSTRAINT +# REPLACE VARIABLES IN CONSTRAINT ################################################################################ +# Replace variable refs in an expression using a map. Uses AbstractDict +# because the InfiniteModel MBM path maps decision vars to VariableRefs +# and parameter functions to Numbers in the same dict (via _build_flat_map). function _replace_variables_in_constraint( - fun:: JuMP.AbstractVariableRef, - var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} -) + fun::JuMP.AbstractVariableRef, + var_map::AbstractDict + ) return var_map[fun] end +# Infer the variable reference type from the map values, falling back to the +# expression's own type. +function _var_ref_type( + ::Type{JuMP.GenericAffExpr{C, V}}, + var_map::AbstractDict + ) where {C, V} + for val in values(var_map) + if val isa JuMP.AbstractVariableRef + return typeof(val) + end + end + return V +end + +# Dispatch for affine/quadratic term addition when var_map values may be Numbers +# (parameter functions evaluated at supports). +_add_aff_term(aff, c, r::Number) = aff.constant += c * r +_add_aff_term(aff, c, r) = JuMP.add_to_expression!(aff, c, r) + function _replace_variables_in_constraint( - fun::T, - var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} -) where {T <: JuMP.GenericAffExpr} - new_aff = JuMP.zero(T) + fun::T, var_map::AbstractDict + ) where {T <: JuMP.GenericAffExpr} + C = JuMP.value_type(T) + W = _var_ref_type(T, var_map) + new_aff = zero(JuMP.GenericAffExpr{C, W}) for (var, coef) in fun.terms - new_var = var_map[var] - JuMP.add_to_expression!(new_aff, coef, new_var) + _add_aff_term(new_aff, coef, var_map[var]) end - new_aff.constant = fun.constant + new_aff.constant = new_aff.constant + fun.constant return new_aff end +_add_quad_term(q, c, ra::Number, rb::Number) = q.aff.constant += c * ra * rb +_add_quad_term(q, c, ra::Number, rb) = JuMP.add_to_expression!(q.aff, c * ra, rb) +_add_quad_term(q, c, ra, rb::Number) = JuMP.add_to_expression!(q.aff, c * rb, ra) +_add_quad_term(q, c, ra, rb) = JuMP.add_to_expression!(q, c, ra, rb) function _replace_variables_in_constraint( - fun::T, - var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} -) where {T <: JuMP.GenericQuadExpr} - new_quad = JuMP.zero(T) + fun::T, var_map::AbstractDict + ) where {T <: JuMP.GenericQuadExpr} + C = JuMP.value_type(T) + W = _var_ref_type(typeof(fun.aff), var_map) + new_quad = zero(JuMP.GenericQuadExpr{C, W}) for (vars, coef) in fun.terms - JuMP.add_to_expression!(new_quad, coef, - var_map[vars.a], var_map[vars.b]) + _add_quad_term(new_quad, coef, var_map[vars.a], var_map[vars.b]) end new_aff = _replace_variables_in_constraint(fun.aff, var_map) JuMP.add_to_expression!(new_quad, new_aff) return new_quad end -function _replace_variables_in_constraint( - fun::Number, - var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} -) +function _replace_variables_in_constraint(fun::Number, var_map::AbstractDict) return fun end -function _replace_variables_in_constraint( - fun::T, - var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} -) where {T <: JuMP.GenericNonlinearExpr} - new_args = Any[_replace_variables_in_constraint(arg, var_map) - for arg in fun.args] +function _replace_variables_in_constraint(fun::T, + var_map::AbstractDict) where {T <: JuMP.GenericNonlinearExpr} + new_args = Any[_replace_variables_in_constraint( + arg, var_map) for arg in fun.args] return T(fun.head, new_args) end -function _replace_variables_in_constraint( - fun::Vector, - var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} -) - return [_replace_variables_in_constraint(expr, var_map) - for expr in fun] -end - -function _replace_variables_in_constraint( - ::F, - ::S -) where {F, S} - error("_replace_variables_in_constraint not implemented for " * - "$(typeof(F)) and $(typeof(S))") +function _replace_variables_in_constraint(fun::Vector, var_map::AbstractDict) + return [_replace_variables_in_constraint(expr, + var_map) for expr in fun] end diff --git a/src/utilities.jl b/src/utilities.jl index 0f211b4c..e9850246 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,7 +1,7 @@ ################################################################################ # MODEL COPYING ################################################################################ -# extentsion point for model copying +# Extension point for model copying (creates empty model). function _copy_model( model::M ) where {M <: JuMP.AbstractModel} diff --git a/src/variables.jl b/src/variables.jl index f96cfe02..e1d2d2b3 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -424,8 +424,7 @@ function _interrogate_variables(interrogator::Function, nlp::JuMP.GenericNonline for arg in nlp.args _interrogate_variables(interrogator, arg) end - # TODO avoid recursion. See InfiniteOpt.jl for alternate method that avoids stackoverflow errors with deeply nested expressions: - # https://github.com/infiniteopt/InfiniteOpt.jl/blob/cb6dd6ae40fe0144b1dd75da0739ea6e305d5357/src/expressions.jl#L520-L534 + return end diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index ed74a5a5..f57e406b 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -1,9 +1,55 @@ using HiGHS function test_mbm() + @test DP._MBM( + DP.MBM(HiGHS.Optimizer), JuMP.Model() + ).optimizer == HiGHS.Optimizer - @test DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()).optimizer == HiGHS.Optimizer + # Test _is_all_zeros + @test DP._is_all_zeros(0) + @test DP._is_all_zeros(0.0) + @test !DP._is_all_zeros(1) + @test !DP._is_all_zeros(5.0) + @test DP._is_all_zeros([0, 0, 0]) + @test DP._is_all_zeros([0.0, 0.0]) + @test !DP._is_all_zeros([0, 1, 0]) + @test !DP._is_all_zeros([5.0, 0.0]) + @test !DP._is_all_zeros("not a number") # non-numeric fallback +end +# _var_ref_type when all map values are Numbers (no VariableRef +# found → returns V). Covers line 537. +function test__var_ref_type_numeric_map() + model = Model() + @variable(model, x[1:2]) + aff = @expression(model, 2*x[1] + 3*x[2]) + var_map = Dict{VariableRef, Any}(x[1] => 5.0, x[2] => 3.0) + @test DP._var_ref_type(typeof(aff), var_map) == VariableRef +end + +# _replace_variables_in_constraint with QuadExpr where var_map +# maps some vars to Numbers. Covers lines 569, 571, 574. +function test__replace_variables_quad_numeric_map() + model = Model() + sub = Model() + @variable(model, x[1:3]) + @variable(sub, y) + quad1 = @expression(model, x[1] * x[2]) + + # both Number (line 569) + map1 = Dict{VariableRef, Any}(x[1] => 2.0, x[2] => 3.0) + result1 = DP._replace_variables_in_constraint(quad1, map1) + @test result1.aff.constant ≈ 6.0 + + # ra Number, rb VariableRef (line 571) + map2 = Dict{VariableRef, Any}(x[1] => 2.0, x[2] => y) + result2 = DP._replace_variables_in_constraint(quad1, map2) + @test result2.aff.terms[y] ≈ 2.0 + + # rb Number, ra VariableRef (line 574) + map3 = Dict{VariableRef, Any}(x[1] => y, x[2] => 3.0) + result3 = DP._replace_variables_in_constraint(quad1, map3) + @test result3.aff.terms[y] ≈ 3.0 end function test__replace_variables_in_constraint() @@ -34,107 +80,241 @@ function test__replace_variables_in_constraint() expected = JuMP.@expression(sub_model, sin(new_vars[x[3]]) - 0.0) @test JuMP.isequal_canonical(expr3, expected) @test expr4 == [new_vars[x[i]] for i in 1:3] - @test_throws ErrorException DP._replace_variables_in_constraint( + @test_throws MethodError DP._replace_variables_in_constraint( "String", new_vars) end -function test__constraint_to_objective() +function test__prepare_objectives() model = Model() sub_model = Model() @variable(model, x[1:2]) - @constraint(model, lessthan, x[1] <= 1) + @constraint(model, lessthan, x[1] <= 1) @constraint(model, greaterthan, x[2] >= 1) - @constraint(model, interval, 0 <= x[1] <= 55) - @constraint(model, equalto, x[1] == 1) - new_vars = Dict{AbstractVariableRef, AbstractVariableRef}() - [new_vars[x[i]] = @variable(sub_model) for i in 1:2] - DP._constraint_to_objective(sub_model, constraint_object(lessthan), - new_vars) - @test objective_function(sub_model) == JuMP.@expression(sub_model, - new_vars[x[1]] - 1) - DP._constraint_to_objective(sub_model, constraint_object(greaterthan), - new_vars) - @test objective_function(sub_model) == JuMP.@expression(sub_model, - 1 - new_vars[x[2]]) - @test_throws ErrorException DP._constraint_to_objective(sub_model, - constraint_object(interval), new_vars) + new_vars = Dict{VariableRef, Vector{VariableRef}}() + for i in 1:2 + new_vars[x[i]] = [@variable(sub_model)] + end + sub = DP.GDPSubmodel(sub_model, + collect(keys(new_vars)), new_vars) + + # LessThan: max(f - upper) = max(x[1] - 1) + objs_le = DP._prepare_objectives( + model, constraint_object(lessthan), sub) + @test length(objs_le) == 1 + @test objs_le[1] == JuMP.@expression(sub_model, + new_vars[x[1]][1] - 1) + + # GreaterThan: max(lower - f) = max(1 - x[2]) + objs_ge = DP._prepare_objectives( + model, constraint_object(greaterthan), sub) + @test length(objs_ge) == 1 + @test objs_ge[1] == JuMP.@expression(sub_model, + 1 - new_vars[x[2]][1]) end -function test_mini_model() +function test_raw_M() model = GDPModel() @variable(model, 0 <= x, start = 1) @variable(model, 0 <= y) - @variable(model, Y[1:4], Logical) - @constraint(model, con, 3*-x <= 4, Disjunct(Y[1])) - @constraint(model, con2, 3*x + y >= 15, Disjunct(Y[2])) - @constraint(model, infeasiblecon, 3*x + y == 15, Disjunct(Y[3])) - @constraint(model, intervalcon, 0 <= x <= 55, Disjunct(Y[4])) - @disjunction(model, [Y[1], Y[2], Y[3], Y[4]]) - mbm = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) - @test DP._mini_model(model, constraint_object(con), - DisjunctConstraintRef[con2], mbm)== -4 + @variable(model, Y[1:5], Logical) + @constraint(model, con, 3*-x <= 4, + Disjunct(Y[1])) + @constraint(model, con2, 3*x + y >= 15, + Disjunct(Y[2])) + @constraint(model, infeasiblecon, + 3*x + y == 15, Disjunct(Y[3])) + @constraint(model, intervalcon, + 0 <= x <= 55, Disjunct(Y[4])) + @constraint(model, truly_infeasible, + x >= 100, Disjunct(Y[5])) + @disjunction(model, + [Y[1], Y[2], Y[3], Y[4], Y[5]]) + mbm = DP._MBM( + DP.MBM(HiGHS.Optimizer), JuMP.Model()) + sub = DP.create_submodel(model, + DisjunctConstraintRef[con2], mbm) + objs = DP._prepare_objectives(model, + constraint_object(con), sub) + raw = DP._raw_M(sub, objs, mbm) + @test DP.condense_values(model, raw) == 0.0 set_upper_bound(x, 1) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con], mbm)== 15 + sub2 = DP.create_submodel(model, + DisjunctConstraintRef[con], mbm) + objs2 = DP._prepare_objectives(model, + constraint_object(con2), sub2) + raw = DP._raw_M(sub2, objs2, mbm) + @test DP.condense_values(model, raw) == 15 set_integer(y) - @constraint(model, con3, y*x == 15, Disjunct(Y[1])) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con], mbm)== 15 + @constraint(model, con3, y*x == 15, + Disjunct(Y[1])) + objs3 = DP._prepare_objectives(model, + constraint_object(con2), sub2) + raw = DP._raw_M(sub2, objs3, mbm) + @test DP.condense_values(model, raw) == 15 + # Fresh _MBM after changing bounds JuMP.fix(y, 5; force=true) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con], mbm)== 10 + mbm2 = DP._MBM( + DP.MBM(HiGHS.Optimizer), JuMP.Model()) + sub3 = DP.create_submodel(model, + DisjunctConstraintRef[con], mbm2) + objs4 = DP._prepare_objectives(model, + constraint_object(con2), sub3) + raw = DP._raw_M(sub3, objs4, mbm2) + @test DP.condense_values(model, raw) == 10 + # Infeasible region → nothing delete_lower_bound(x) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con2], mbm) == 1.0e9 + mbm3 = DP._MBM( + DP.MBM(HiGHS.Optimizer), JuMP.Model()) + sub4 = DP.create_submodel(model, + DisjunctConstraintRef[con2], mbm3) + objs5 = DP._prepare_objectives(model, + constraint_object(con2), sub4) + @test DP._raw_M(sub4, objs5, mbm3) == nothing + + # infeasible (x >= 100 but x <= 1) + set_upper_bound(x, 1) + mbm4 = DP._MBM( + DP.MBM(HiGHS.Optimizer), JuMP.Model()) + sub5 = DP.create_submodel(model, + DisjunctConstraintRef[truly_infeasible], + mbm4) + objs6 = DP._prepare_objectives(model, + constraint_object(con), sub5) + @test DP._raw_M(sub5, objs6, mbm4) == nothing + + # Unbounded subproblem → default_M fallback. + # No lower bound on x means max(5 - x) s.t. x <= 3 + # is unbounded (DUAL_INFEASIBLE). + model_ub = GDPModel() + @variable(model_ub, xu) # no bounds + @variable(model_ub, Yu[1:2], Logical) + @constraint(model_ub, ub_con1, xu <= 3, Disjunct(Yu[1])) + @constraint(model_ub, ub_con2, xu >= 5, Disjunct(Yu[2])) + @disjunction(model_ub, [Yu[1], Yu[2]]) + mbm_ub = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + sub_ub = DP.create_submodel(model_ub, + DisjunctConstraintRef[ub_con1], mbm_ub) + objs_ub = DP._prepare_objectives(model_ub, + constraint_object(ub_con2), sub_ub) + raw_ub = DP._raw_M(sub_ub, objs_ub, mbm_ub) + @test raw_ub == [mbm_ub.default_M] end function test_maximize_M() model = GDPModel() - @variable(model, 0 <= x[1:2] <= 50) + # Different bounds for x[1] and x[2] to demonstrate per-row M values + @variable(model, x[1:2]) + set_lower_bound(x[1], 0); set_upper_bound(x[1], 10) + set_lower_bound(x[2], 0); set_upper_bound(x[2], 5) @variable(model, Y[1:6], Logical) @constraint(model, lessthan, x[1] <= 1, Disjunct(Y[1])) @constraint(model, greaterthan, x[1] >= 1, Disjunct(Y[1])) @constraint(model, interval, 0 <= x[1] <= 55, Disjunct(Y[2])) @constraint(model, equalto, x[1] == 1, Disjunct(Y[3])) - @constraint(model, nonpositives, -x in MOI.Nonpositives(2), + # Vector constraints: x >= 0 (both rows) + @constraint(model, nonpositives, -x in MOI.Nonpositives(2), Disjunct(Y[4])) - @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), + @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), Disjunct(Y[5])) + # Vector equality: x == 1 (both rows) @constraint(model, zeros, -x .+ 1 in MOI.Zeros(2), Disjunct(Y[6])) mbm = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) - @test DP._maximize_M(model, constraint_object(interval), + + # Interval returns [M_lower, M_upper] + # M_lower = max(0 - x[1]) s.t. 0<=x[1]<=10 = 0 at x[1]=0 + # M_upper = max(x[1] - 55) s.t. 0<=x[1]<=10 = -45 at x[1]=10, clamped to 0 + @test DP._maximize_M(model, constraint_object(interval), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 0.0 - @test DP._maximize_M(model, constraint_object(lessthan), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [0.0, 0.0] + + # Scalar LessThan/GreaterThan still return scalars + # lessthan: x[1] <= 1 vs interval 0 <= x[1] <= 55 + # max(x[1] - 1) s.t. 0<=x[1]<=10 = 9 at x[1]=10 + @test DP._maximize_M(model, constraint_object(lessthan), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 49 - @test DP._maximize_M(model, constraint_object(greaterthan), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == 9.0 + + # greaterthan: x[1] >= 1 vs interval 0 <= x[1] <= 55 + # max(1 - x[1]) s.t. 0<=x[1]<=10 = 1 at x[1]=0 + @test DP._maximize_M(model, constraint_object(greaterthan), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), + DP._indicator_to_constraints(model)[Y[2]]), mbm) == 1.0 - @test DP._maximize_M(model, constraint_object(equalto), + + # EqualTo returns [M_lower, M_upper] + @test DP._maximize_M(model, constraint_object(equalto), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[3]]), - mbm) == 0 - @test DP._maximize_M(model, constraint_object(nonpositives), + DP._indicator_to_constraints(model)[Y[3]]), + mbm) == [0.0, 0.0] + + # Vector constraints: per-row M values + # nonpositives: x >= 0 against Y[2] (interval only on x[1]) + # Row 1: max(0 - x[1]) s.t. 0<=x[1]<=10 = 0 + # Row 2: max(0 - x[2]) s.t. 0<=x[2]<=5 = 0 + @test DP._maximize_M(model, constraint_object(nonpositives), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 0 - @test DP._maximize_M(model, constraint_object(nonnegatives), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [0.0, 0.0] + + @test DP._maximize_M(model, constraint_object(nonnegatives), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 0 - @test DP._maximize_M(model, constraint_object(zeros), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [0.0, 0.0] + + # Row 1: max(|x[1] - 1|) s.t. 0<=x[1]<=10 = max(9, 1) = 9 + # Row 2: max(|x[2] - 1|) s.t. 0<=x[2]<=5 = max(4, 1) = 4 + @test DP._maximize_M(model, constraint_object(zeros), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 49 - @test_throws ErrorException DP._maximize_M(model, "odd", + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [9.0, 4.0] + + @test_throws ErrorException DP._maximize_M(model, "odd", Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), + DP._indicator_to_constraints(model)[Y[2]]), mbm) + + # Add an infeasible disjunct (x >= 100 but bounds are 0-10) + @variable(model, Y_infeas, Logical) + @constraint(model, infeas_con, x[1] >= 100, Disjunct(Y_infeas)) + + # Scalar constraint against infeasible disjunct -> nothing + @test DP._maximize_M(model, constraint_object(lessthan), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Vector constraint (Nonpositives) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(nonpositives), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Vector constraint (Nonnegatives) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(nonnegatives), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Vector constraint (Zeros) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(zeros), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Bidirectional scalar (EqualTo) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(equalto), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Bidirectional scalar (Interval) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(interval), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing end function test_reformulate_disjunct_constraint() @@ -144,57 +324,97 @@ function test_reformulate_disjunct_constraint() @constraint(model, lessthan, x[1] <= 1, Disjunct(Y[1])) @constraint(model, greaterthan, x[1] >= 1, Disjunct(Y[1])) @constraint(model, equalto, x[1] == 1, Disjunct(Y[2])) - @constraint(model, nonpositives, -x in MOI.Nonpositives(2), + @constraint(model, nonpositives, -x in MOI.Nonpositives(2), Disjunct(Y[3])) - @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), + @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), Disjunct(Y[4])) @constraint(model, zeros, -x .+ 1 in MOI.Zeros(2), Disjunct(Y[5])) @disjunction(model, disjunction,[Y[1], Y[2], Y[3], Y[4], Y[5]]) - method = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + bconref = Dict(Y[i] => binary_variable(Y[i]) for i in 1:5) + + # Test scalar constraints (LessThan, GreaterThan) with scalar M values + method_scalar = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) for i in 1:5 - method.M[Y[i]] = Float64(i) + method_scalar.M[Y[i]] = Float64(i) end - bconref = Dict(Y[i] => binary_variable(Y[i]) for i in 1:5) - reformulated_constraints = [reformulate_disjunct_constraint(model, - constraint_object(constraints), bconref, method) - for constraints in [lessthan, greaterthan, equalto, nonpositives, - nonnegatives, zeros, disjunction]] - @test reformulated_constraints[1][1].func == JuMP.@expression(model, - x[1] - sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[1][1].set == MOI.LessThan(1.0) - @test reformulated_constraints[2][1].func == JuMP.@expression(model, - x[1] + sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[2][1].set == MOI.GreaterThan(1.0) - @test reformulated_constraints[3][1].func == JuMP.@expression(model, - x[1] + sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[3][1].set == MOI.GreaterThan(1.0) - @test reformulated_constraints[3][2].func == JuMP.@expression(model, - x[1] - sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[3][2].set == MOI.LessThan(1.0) - @test reformulated_constraints[4][1].func == JuMP.@expression(model, - -x .- sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[4][1].set == MOI.Nonpositives(2) - @test reformulated_constraints[5][1].func == JuMP.@expression(model, - x .+ sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[5][1].set == MOI.Nonnegatives(2) - @test reformulated_constraints[6][1].func == JuMP.@expression(model, - -x .+(1 + sum(method.M[i] * bconref[i] for i in keys(method.M)))) && - reformulated_constraints[6][1].set == MOI.Nonnegatives(2) - @test reformulated_constraints[6][2].func == JuMP.@expression(model, - -x .+(1 - sum(method.M[i] * bconref[i] for i in keys(method.M)))) && - reformulated_constraints[6][2].set == MOI.Nonpositives(2) - @test reformulated_constraints[7][1].func == JuMP.@expression(model, - x[1] - 52*bconref[Y[3]] - 53*bconref[Y[4]] - bconref[Y[1]] - - 5*bconref[Y[5]] - 2*bconref[Y[2]]) && - reformulated_constraints[7][1].set == MOI.LessThan(1.0) - @test reformulated_constraints[7][2].func == JuMP.@expression(model, - x[1] + 52*bconref[Y[3]] + 53*bconref[Y[4]] + bconref[Y[1]] - + 5*bconref[Y[5]] + 2*bconref[Y[2]]) && - reformulated_constraints[7][2].set == MOI.GreaterThan(1.0) - - @test_throws ErrorException reformulate_disjunct_constraint(model, - "odd", bconref, method) + ref_lessthan = reformulate_disjunct_constraint( + model, constraint_object(lessthan), bconref, method_scalar) + ref_greaterthan = reformulate_disjunct_constraint( + model, constraint_object(greaterthan), bconref, method_scalar) + @test ref_lessthan[1].func == JuMP.@expression(model, + x[1] - sum(method_scalar.M[i] * bconref[i] + for i in keys(method_scalar.M))) && + ref_lessthan[1].set == MOI.LessThan(1.0) + @test ref_greaterthan[1].func == JuMP.@expression(model, + x[1] + sum(method_scalar.M[i] * bconref[i] + for i in keys(method_scalar.M))) && + ref_greaterthan[1].set == MOI.GreaterThan(1.0) + # Test bidirectional constraint (EqualTo) with [M_lower, M_upper] values + method_equalto = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + for i in 1:5 + method_equalto.M[Y[i]] = [Float64(i), Float64(i)] + end + ref_equalto = reformulate_disjunct_constraint( + model, constraint_object(equalto), bconref, method_equalto) + @test ref_equalto[1].func == JuMP.@expression(model, + x[1] + sum(method_equalto.M[i][1] * bconref[i] + for i in keys(method_equalto.M))) && + ref_equalto[1].set == MOI.GreaterThan(1.0) + @test ref_equalto[2].func == JuMP.@expression(model, + x[1] - sum(method_equalto.M[i][2] * bconref[i] + for i in keys(method_equalto.M))) && + ref_equalto[2].set == MOI.LessThan(1.0) + + # Test vector constraints with per-row M values [M_row1, M_row2] + method_vector = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + for i in 1:5 + method_vector.M[Y[i]] = [Float64(i), Float64(i)] + end + ref_nonpositives = reformulate_disjunct_constraint( + model, constraint_object(nonpositives), bconref, method_vector) + ref_nonnegatives = reformulate_disjunct_constraint( + model, constraint_object(nonnegatives), bconref, method_vector) + ref_zeros = reformulate_disjunct_constraint( + model, constraint_object(zeros), bconref, method_vector) + @test ref_nonpositives[1].func == JuMP.@expression(model, [j=1:2], + -x[j] - sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_nonpositives[1].set == MOI.Nonpositives(2) + @test ref_nonnegatives[1].func == JuMP.@expression(model, [j=1:2], + x[j] + sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_nonnegatives[1].set == MOI.Nonnegatives(2) + @test ref_zeros[1].func == JuMP.@expression(model, [j=1:2], + -x[j] + 1 + sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_zeros[1].set == MOI.Nonnegatives(2) + @test ref_zeros[2].func == JuMP.@expression(model, [j=1:2], + -x[j] + 1 - sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_zeros[2].set == MOI.Nonpositives(2) + + # Test nested disjunction reformulation with proper nested structure + # Create outer disjunct with inner disjunction + model2 = GDPModel() + @variable(model2, 0 <= z <= 50) + @variable(model2, Outer[1:2], Logical) + @variable(model2, Inner[1:2], Logical) + @constraint(model2, inner_lt, z <= 1, Disjunct(Inner[1])) + @constraint(model2, inner_gt, z >= 1, Disjunct(Inner[2])) + @disjunction(model2, inner_disj, Inner) + method_nested = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + bconref2 = Dict(Outer[2] => binary_variable(Outer[2])) + method_nested.M[Outer[2]] = 10.0 #Dummy M value for testing. + #Normally _reformulate_disjunct will this without having to assign a value + ref_disjunction = reformulate_disjunct_constraint( + model2, constraint_object(inner_disj), bconref2, method_nested) + @test length(ref_disjunction) >= 2 + @test JuMP.coefficient(ref_disjunction[1].func, z) == 1.0 + @test JuMP.coefficient(ref_disjunction[2].func, z) == 1.0 + + @test_throws ErrorException reformulate_disjunct_constraint(model, + "odd", bconref, method_scalar) end function test_reformulate_disjunct() @@ -221,13 +441,14 @@ function test_reformulate_disjunct() func_3 = reformulated_disjunct[3].func @test JuMP.coefficient(func_1, x[1]) == 1.0 - @test JuMP.coefficient(func_1, binary_variable(Y[2])) == -1.5 + @test JuMP.coefficient(func_1, binary_variable(Y[2])) == 0.0 + # Per-bound M: lower bound uses M_lower=1.5, upper bound uses M_upper=2.5 @test JuMP.coefficient(func_2, x[1]) == 1.0 - @test JuMP.coefficient(func_2, binary_variable(Y[1])) == 2.5 + @test JuMP.coefficient(func_2, binary_variable(Y[1])) == 1.5 # M_lower @test JuMP.coefficient(func_3, x[1]) == 1.0 - @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -2.5 + @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -2.5 # -M_upper end function test_reformulate_disjunction() @@ -238,43 +459,374 @@ function test_reformulate_disjunction() @constraint(model, greaterthan, x >= 1, Disjunct(Y[1])) @constraint(model, interval, 0 <= x <= 55, Disjunct(Y[2])) disj = disjunction(model, [Y[1], Y[2]]) - + method = DP.MBM(HiGHS.Optimizer) ref_cons = reformulate_disjunction(model, constraint_object(disj), method) - @test length(ref_cons) == 4 + # 3 constraints: lessthan, greaterthan (with Big-M), interval (global) + @test length(ref_cons) == 3 @test ref_cons[1].set == MOI.LessThan(2.0) - @test ref_cons[2].set == MOI.GreaterThan(1.0) - - @test ref_cons[3].set == MOI.GreaterThan(0.0) - - @test ref_cons[4].set == MOI.LessThan(55.0) + # Interval is global (M=0 for both bounds in Y[1]'s region 1<=x<=2) + @test ref_cons[3].set == MOI.Interval(0.0, 55.0) - func_1 = ref_cons[1].func # x - 53 Y[2] <= 2.0 - func_2 = ref_cons[2].func # x + 53 Y[2] >= 1.0 - func_3 = ref_cons[3].func # x - Y[1] >= 0.0 - func_4 = ref_cons[4].func # x + Y[1] <= 55.0 + # Per-constraint, per-bound M values: + # - lessthan (x <= 2) in Y[2] region (0 <= x <= 55): max(x-2) at + # x=55 -> M=53 + # - greaterthan (x >= 1) in Y[2] region: max(1-x) at x=0 -> M=1 + # - interval in Y[1] region (1 <= x <= 2): + # - M_lower (x >= 0): max(0-x) at x=1 -> M_lower=-1, clamped to 0 + # - M_upper (x <= 55): max(x-55) at x=2 -> M_upper=-53, clamped to 0 + # - Both M=0 -> detected as global, added without Big-M + func_1 = ref_cons[1].func # x - 53*Y[2] <= 2.0 + func_2 = ref_cons[2].func # x + 1*Y[2] >= 1.0 + func_3 = ref_cons[3].func # x (global, no binary variables) @test JuMP.coefficient(func_1, x) == 1.0 @test JuMP.coefficient(func_1, binary_variable(Y[2])) == -53.0 @test JuMP.coefficient(func_2, x) == 1.0 - @test JuMP.coefficient(func_2, binary_variable(Y[2])) == 53.0 + @test JuMP.coefficient(func_2, binary_variable(Y[2])) == 1.0 + # Global constraint has just x, no binary variables @test JuMP.coefficient(func_3, x) == 1.0 - @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -1.0 - @test JuMP.coefficient(func_4, x) == 1.0 - @test JuMP.coefficient(func_4, binary_variable(Y[1])) == 1.0 + #Test infeasible disjunct detection and deactivation + model2 = GDPModel() + @variable(model2, 0 <= z <= 1) + @variable(model2, W[1:2], Logical) + # W[1]: z >= 5 is infeasible (z has upper bound 1) + @constraint(model2, z >= 5, Disjunct(W[1])) + # W[2]: z <= 0.5 is feasible + @constraint(model2, z <= 0.5, Disjunct(W[2])) + disj2 = disjunction(model2, [W[1], W[2]]) + + method2 = DP.MBM(HiGHS.Optimizer) + ref_cons2 = @test_logs (:warn, r"infeasible, deactivating") begin + reformulate_disjunction(model2, constraint_object(disj2), method2) + end + + @test length(ref_cons2) == 1 + @test ref_cons2[1].set == MOI.LessThan(0.5) + + #Test multiple infeasible disjuncts + model3 = GDPModel() + @variable(model3, 0 <= w <= 1) + @variable(model3, V[1:3], Logical) + @constraint(model3, w >= 5, Disjunct(V[1])) # infeasible + @constraint(model3, w >= 10, Disjunct(V[2])) # infeasible + @constraint(model3, w <= 0.5, Disjunct(V[3])) # feasible + disj3 = disjunction(model3, [V[1], V[2], V[3]]) + + method3 = DP.MBM(HiGHS.Optimizer) + #warn about V[1] and V[2] being infeasible + ref_cons3 = @test_logs (:warn,) (:warn,) begin + reformulate_disjunction(model3, constraint_object(disj3), method3) + end + + # Only V[3]'s constraint should be reformulated + @test length(ref_cons3) == 1 + @test ref_cons3[1].set == MOI.LessThan(0.5) + + model4 = GDPModel() + @variable(model4, 0 <= u <= 10) + @variable(model4, U[1:2], Logical) + @constraint(model4, u <= 3, Disjunct(U[1])) + @constraint(model4, u >= 5, Disjunct(U[2])) + disj4 = disjunction(model4, [U[1], U[2]]) + + method4 = DP.MBM(HiGHS.Optimizer) + ref_cons4 = @test_nowarn begin + reformulate_disjunction(model4, constraint_object(disj4), method4) + end + @test length(ref_cons4) == 2 + + # Disjunct 1: x <= 10 (will be global because D2's region has x <= 5) + # Disjunct 2: x <= 5 + # max(x - 10) s.t. x <= 5 = 5 - 10 = -5, clamped to 0 + # So M = 0 for D1's constraint -> it's global + model5 = GDPModel() + @variable(model5, 0 <= g <= 10) + @variable(model5, G[1:2], Logical) + @constraint(model5, g <= 10, Disjunct(G[1])) # global: other region is g<=5 + @constraint(model5, g <= 5, Disjunct(G[2])) + disj5 = disjunction(model5, [G[1], G[2]]) + + method5 = DP.MBM(HiGHS.Optimizer) + ref_cons5 = reformulate_disjunction(model5, constraint_object(disj5), method5) + # G[1]'s constraint is global (added without Big-M) + # G[2]'s constraint has M > 0 (needs Big-M): max(g - 5) s.t. g<=10 = 5 + @test length(ref_cons5) == 2 + # Check that one constraint has binary variable coefficient (Big-M term) + # and one doesn't (global) + global_con = ref_cons5[1] # g <= 10 (global) + bigm_con = ref_cons5[2] # g <= 5 with Big-M + @test global_con.set == MOI.LessThan(10.0) + @test bigm_con.set == MOI.LessThan(5.0) + # Global constraint should have no binary variable coefficient + @test JuMP.coefficient(global_con.func, binary_variable(G[2])) == 0.0 + # Big-M constraint should have binary variable coefficient = -5 + @test JuMP.coefficient(bigm_con.func, binary_variable(G[1])) == -5.0 +end + +################################################################################ +# LOW-LEVEL UNIT TESTS FOR MBM INFRASTRUCTURE +################################################################################ + +# Test _copy_model directly +function test__copy_model() + # Test with standard JuMP Model + model = Model() + @variable(model, x >= 0) + @variable(model, y <= 10) + copied = DP._copy_model(model) + @test copied isa Model + @test num_variables(copied) == 0 # _copy_model creates empty model + + # Test with GDPModel + gdp_model = GDPModel() + @variable(gdp_model, z) + copied_gdp = DP._copy_model(gdp_model) + @test copied_gdp isa Model + @test num_variables(copied_gdp) == 0 +end + +# Test VariableProperties struct and constructors +function test_variable_properties() + # Test VariableProperties(vref::GenericVariableRef) - standard JuMP variable + model = GDPModel() + @variable(model, 0 <= x <= 10, start = 5) + props_x = DP.VariableProperties(x) + @test props_x.info.has_lb == true + @test props_x.info.lower_bound == 0 + @test props_x.info.has_ub == true + @test props_x.info.upper_bound == 10 + @test props_x.info.has_start == true + @test props_x.info.start == 5 + @test props_x.name == "x" + @test props_x.variable_type === nothing + + # Test with binary variable + @variable(model, y, Bin) + props_y = DP.VariableProperties(y) + @test props_y.info.binary == true + @test props_y.info.integer == false + + # Test with integer variable + @variable(model, z, Int) + props_z = DP.VariableProperties(z) + @test props_z.info.binary == false + @test props_z.info.integer == true + + # Test with fixed variable + @variable(model, w == 42) + props_w = DP.VariableProperties(w) + @test props_w.info.has_fix == true + @test props_w.info.fixed_value == 42 + + # Test VariableProperties(expr) - blank info constructor + expr = 2*x + 3*y + props_expr = DP.VariableProperties(expr) + @test props_expr.info.has_lb == false + @test props_expr.info.has_ub == false + @test props_expr.info.has_fix == false + @test props_expr.info.has_start == false + @test props_expr.info.binary == false + @test props_expr.info.integer == false + @test props_expr.name == "" +end + +# Test _make_variable_object +function test__make_variable_object() + model = GDPModel() + @variable(model, 0 <= x <= 10, start = 5) + + # Create VariableProperties and then make variable object + props = DP.VariableProperties(x) + var_obj = DP._make_variable_object(props) + @test var_obj isa JuMP.ScalarVariable + @test var_obj.info.has_lb == true + @test var_obj.info.lower_bound == 0 + @test var_obj.info.has_ub == true + @test var_obj.info.upper_bound == 10 + + # Test with blank properties (from expression) + props_blank = DP.VariableProperties(2*x) + var_obj_blank = DP._make_variable_object(props_blank) + @test var_obj_blank isa JuMP.ScalarVariable + @test var_obj_blank.info.has_lb == false + @test var_obj_blank.info.has_ub == false +end + +# Test create_variable +function test_create_variable() + model = GDPModel() + @variable(model, 0 <= x <= 10, start = 5) + + # Create a new variable from properties + props = DP.VariableProperties(x) + new_var = DP.create_variable(model, props) + @test new_var isa VariableRef + @test has_lower_bound(new_var) + @test lower_bound(new_var) == 0 + @test has_upper_bound(new_var) + @test upper_bound(new_var) == 10 + @test start_value(new_var) == 5 + @test name(new_var) == "x" + + # Test with binary variable + @variable(model, y, Bin) + props_bin = DP.VariableProperties(y) + new_bin = DP.create_variable(model, props_bin) + @test is_binary(new_bin) + + # Test with integer variable + @variable(model, z, Int) + props_int = DP.VariableProperties(z) + new_int = DP.create_variable(model, props_int) + @test is_integer(new_int) + + # Test with fixed variable + @variable(model, w == 42) + props_fix = DP.VariableProperties(w) + new_fix = DP.create_variable(model, props_fix) + @test is_fixed(new_fix) + @test fix_value(new_fix) == 42 +end + +# Test variable_copy +function test_variable_copy() + source_model = GDPModel() + @variable(source_model, 0 <= x <= 10, start = 5) + @variable(source_model, y, Bin) + @variable(source_model, z, Int) + @variable(source_model, w == 42) + + target_model = GDPModel() + + # Copy bounded variable + x_copy = DP.variable_copy(target_model, x) + @test x_copy isa VariableRef + @test owner_model(x_copy) === target_model + @test has_lower_bound(x_copy) + @test lower_bound(x_copy) == 0 + @test has_upper_bound(x_copy) + @test upper_bound(x_copy) == 10 + @test start_value(x_copy) == 5 + @test name(x_copy) == "x" + + # Copy binary variable + y_copy = DP.variable_copy(target_model, y) + @test is_binary(y_copy) + @test owner_model(y_copy) === target_model + + # Copy integer variable + z_copy = DP.variable_copy(target_model, z) + @test is_integer(z_copy) + @test owner_model(z_copy) === target_model + + # Copy fixed variable + w_copy = DP.variable_copy(target_model, w) + @test is_fixed(w_copy) + @test fix_value(w_copy) == 42 + @test owner_model(w_copy) === target_model +end + +# Test _create_submodel (combines _copy_model + variable_copy + constraint handling) +function test__create_submodel() + model = GDPModel() + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 5) + @variable(model, Y[1:2], Logical) + @constraint(model, con1, x + y <= 8, Disjunct(Y[1])) + @constraint(model, con2, x - y >= 2, Disjunct(Y[2])) + @disjunction(model, [Y[1], Y[2]]) + + mbm = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + constraints = Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y[1]] + ) + + sub = DP.create_submodel(model, constraints, mbm) + + # Check submodel struct + @test sub isa DP.GDPSubmodel + @test sub.model isa Model + # Check variable mapping exists + @test haskey(sub.fwd, x) + @test haskey(sub.fwd, y) + + # Check mapped variables (length-1 vectors) + @test length(sub.fwd[x]) == 1 + @test length(sub.fwd[y]) == 1 + + # Check bounds in submodel + xm = sub.fwd[x][1] + ym = sub.fwd[y][1] + @test has_lower_bound(xm) + @test lower_bound(xm) == 0 + @test has_upper_bound(xm) + @test upper_bound(xm) == 10 + + @test has_lower_bound(ym) + @test lower_bound(ym) == 0 + @test has_upper_bound(ym) + @test upper_bound(ym) == 5 + + # Check constraint was added to submodel + @test num_constraints(sub.model, AffExpr, + MOI.LessThan{Float64}) == 1 +end + +# Test get_variable_info +function test_get_variable_info() + model = GDPModel() + @variable(model, 0 <= x <= 10, start = 5) + @variable(model, y, Bin) + @variable(model, z == 42) + + # Test bounded variable + info_x = DP.get_variable_info(x) + @test info_x.has_lb == true + @test info_x.lower_bound == 0 + @test info_x.has_ub == true + @test info_x.upper_bound == 10 + @test info_x.has_start == true + @test info_x.start == 5 + @test info_x.binary == false + @test info_x.integer == false + + # Test binary variable + info_y = DP.get_variable_info(y) + @test info_y.binary == true + @test info_y.integer == false + + # Test fixed variable + info_z = DP.get_variable_info(z) + @test info_z.has_fix == true + @test info_z.fixed_value == 42 + + # Test with overridden kwargs + info_custom = DP.get_variable_info(x; has_lb = false, has_ub = false) + @test info_custom.has_lb == false + @test info_custom.has_ub == false end @testset "MBM" begin + test__copy_model() + test_variable_properties() + test__make_variable_object() + test_create_variable() + test_variable_copy() + test__create_submodel() + test_get_variable_info() test_mbm() + test__var_ref_type_numeric_map() + test__replace_variables_quad_numeric_map() test__replace_variables_in_constraint() - test__constraint_to_objective() - test_mini_model() + test__prepare_objectives() + test_raw_M() test_maximize_M() test_reformulate_disjunct_constraint() test_reformulate_disjunct() From 2270d22fe21d92ff3b2959b6d7435147aba1a51c Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Mon, 23 Mar 2026 00:10:59 -0400 Subject: [PATCH 02/13] . --- src/mbm.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 9d731740..77d0a2a2 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -512,11 +512,6 @@ function _var_ref_type( return V end -# Dispatch for affine/quadratic term addition when var_map values may be Numbers -# (parameter functions evaluated at supports). -_add_aff_term(aff, c, r::Number) = aff.constant += c * r -_add_aff_term(aff, c, r) = JuMP.add_to_expression!(aff, c, r) - function _replace_variables_in_constraint( fun::T, var_map::AbstractDict ) where {T <: JuMP.GenericAffExpr} @@ -524,11 +519,15 @@ function _replace_variables_in_constraint( W = _var_ref_type(T, var_map) new_aff = zero(JuMP.GenericAffExpr{C, W}) for (var, coef) in fun.terms - _add_aff_term(new_aff, coef, var_map[var]) + JuMP.add_to_expression!(new_aff, coef, var_map[var]) end new_aff.constant = new_aff.constant + fun.constant return new_aff end + +# Dispatch for quadratic term addition when var_map values may be Numbers +# (parameter functions evaluated at supports). JuMP's 3-arg +# add_to_expression!(quad, c, ra, rb) doesn't support Number×Number. _add_quad_term(q, c, ra::Number, rb::Number) = q.aff.constant += c * ra * rb _add_quad_term(q, c, ra::Number, rb) = JuMP.add_to_expression!(q.aff, c * ra, rb) _add_quad_term(q, c, ra, rb::Number) = JuMP.add_to_expression!(q.aff, c * rb, ra) From 47a73045a46c2f0274e7c5ca4777f5629c99039f Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Mon, 23 Mar 2026 00:15:28 -0400 Subject: [PATCH 03/13] . --- src/mbm.jl | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 77d0a2a2..a2305ef5 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -525,14 +525,6 @@ function _replace_variables_in_constraint( return new_aff end -# Dispatch for quadratic term addition when var_map values may be Numbers -# (parameter functions evaluated at supports). JuMP's 3-arg -# add_to_expression!(quad, c, ra, rb) doesn't support Number×Number. -_add_quad_term(q, c, ra::Number, rb::Number) = q.aff.constant += c * ra * rb -_add_quad_term(q, c, ra::Number, rb) = JuMP.add_to_expression!(q.aff, c * ra, rb) -_add_quad_term(q, c, ra, rb::Number) = JuMP.add_to_expression!(q.aff, c * rb, ra) -_add_quad_term(q, c, ra, rb) = JuMP.add_to_expression!(q, c, ra, rb) - function _replace_variables_in_constraint( fun::T, var_map::AbstractDict ) where {T <: JuMP.GenericQuadExpr} @@ -540,7 +532,8 @@ function _replace_variables_in_constraint( W = _var_ref_type(typeof(fun.aff), var_map) new_quad = zero(JuMP.GenericQuadExpr{C, W}) for (vars, coef) in fun.terms - _add_quad_term(new_quad, coef, var_map[vars.a], var_map[vars.b]) + JuMP.add_to_expression!(new_quad, + coef * var_map[vars.a] * var_map[vars.b]) end new_aff = _replace_variables_in_constraint(fun.aff, var_map) JuMP.add_to_expression!(new_quad, new_aff) From 2da191fcd0afe0be118a5e4e7332b46fc3638332 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Mon, 23 Mar 2026 00:31:56 -0400 Subject: [PATCH 04/13] . --- src/mbm.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index a2305ef5..2b845d9e 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -283,16 +283,16 @@ function _raw_M( M_vals = typeof(method.default_M)[] for obj_expr in objectives JuMP.@objective(sub.model, Max, obj_expr) + # Clear start values before each solve to prevent NaN + # residue from a previous non-feasible solve + for v in JuMP.all_variables(sub.model) + JuMP.set_start_value(v, nothing) + end JuMP.optimize!(sub.model) if JuMP.termination_status(sub.model) == _MOI.INFEASIBLE return nothing elseif !JuMP.is_solved_and_feasible(sub.model) push!(M_vals, method.default_M) - # Clear NaN start values from non-feasible solve - # so the next objective doesn't inherit them - for v in JuMP.all_variables(sub.model) - JuMP.set_start_value(v, nothing) - end else push!(M_vals, max( JuMP.objective_value(sub.model), From 01807a913d507458a8cd5ae120f5d48ed492fc7f Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Mon, 23 Mar 2026 09:09:23 -0400 Subject: [PATCH 05/13] . --- src/variables.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/variables.jl b/src/variables.jl index e1d2d2b3..f96cfe02 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -424,7 +424,8 @@ function _interrogate_variables(interrogator::Function, nlp::JuMP.GenericNonline for arg in nlp.args _interrogate_variables(interrogator, arg) end - + # TODO avoid recursion. See InfiniteOpt.jl for alternate method that avoids stackoverflow errors with deeply nested expressions: + # https://github.com/infiniteopt/InfiniteOpt.jl/blob/cb6dd6ae40fe0144b1dd75da0739ea6e305d5357/src/expressions.jl#L520-L534 return end From 0403594714a7ebc0037c28181b56697048f676d6 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Mon, 23 Mar 2026 09:28:15 -0400 Subject: [PATCH 06/13] . --- src/datatypes.jl | 19 +++++++++++--- src/mbm.jl | 56 +++++++++++++++++++++++++++-------------- src/variables.jl | 3 +-- test/constraints/mbm.jl | 22 ++++++++-------- 4 files changed, 65 insertions(+), 35 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index f37c8e2d..c75e54f1 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -472,9 +472,22 @@ end # GDP SUBMODEL ################################################################################ -# Unified submodel wrapper for MBM and cutting planes. -# Holds a flat JuMP model, ordered decision variables, -# and a forward map (orig var → submodel vars). +""" + GDPSubmodel{M, V, W} + +A unified submodel wrapper used by MBM and cutting plane +reformulations. It encapsulates a flat JuMP optimization +submodel built from a single disjunct's feasible region, +along with mappings back to the original model's variables. + +## Fields +- `model::M`: The JuMP submodel representing a disjunct's + feasible region (constraints and variable bounds). +- `dec_vars::Vector{V}`: Ordered decision variables in + the submodel, matching the original model's ordering. +- `fwd::Dict{V, Vector{W}}`: Forward map from original + model variables to their submodel counterparts. +""" struct GDPSubmodel{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, W <: JuMP.AbstractVariableRef} diff --git a/src/mbm.jl b/src/mbm.jl index 2b845d9e..025c1427 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -9,10 +9,10 @@ _is_all_zeros(::Any) = false ################################################################################ # CONSTRAINT, DISJUNCTION, DISJUNCT REFORMULATION ################################################################################ -# Reformulates the disjunction using multiple big-M values per constraint + function reformulate_disjunction( model::JuMP.AbstractModel, - disj::Disjunction, + disj::Disjunction, method::MBM ) mbm = _MBM(method, model) @@ -251,9 +251,14 @@ end # MULTIPLE BIG-M REFORMULATION ################################################################################ -# Prepare flat objectives for _raw_M. Returns a vector of objective expressions -# ready to maximize. Base: single flat constraint via fwd[v][1]. -function _prepare_objectives( +""" + prepare_objectives(model, obj::ScalarConstraint, sub::GDPSubmodel) + +Convert a constraint into objective expressions for M-value maximization. +Returns a vector of JuMP expressions to pass to `_raw_M`. The base method +produces a single-element vector by mapping variables through `sub.fwd`. +""" +function prepare_objectives( ::JuMP.AbstractModel, obj::JuMP.ScalarConstraint{T, S}, sub::GDPSubmodel @@ -263,7 +268,7 @@ function _prepare_objectives( return [expr] end -function _prepare_objectives( +function prepare_objectives( ::JuMP.AbstractModel, obj::JuMP.ScalarConstraint{T, S}, sub::GDPSubmodel @@ -303,8 +308,15 @@ function _raw_M( return M_vals end -# Condense flat per-support values to final form. Base: return -# scalar from single-element vector. Extensions may override. +""" + condense_values(model, vals::AbstractVector) + +Reduce a vector of raw M values from `_raw_M` to the final form used +during constraint reformulation. The base method returns the single +element from a length-1 vector. Extensions (e.g., InfiniteOpt) override +to aggregate across multiple support points (e.g., interpolating K +per-support M values into a parameter function). +""" function condense_values( ::JuMP.AbstractModel, vals::AbstractVector @@ -321,7 +333,7 @@ function _maximize_M( method::_MBM ) where {T, S <: Union{_MOI.LessThan, _MOI.GreaterThan}} sub = _get_submodel(model, constraints, method) - objectives = _prepare_objectives(model, objective, sub) + objectives = prepare_objectives(model, objective, sub) raw = _raw_M(sub, objectives, method) raw === nothing && return nothing return condense_values(model, raw) @@ -353,8 +365,8 @@ function _maximize_M( set_value = objective.set.value ge_obj = JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)) le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)) - raw_lower = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) - raw_upper = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + raw_lower = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) + raw_upper = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) (raw_lower === nothing || raw_upper === nothing) && return nothing return [condense_values(model, raw_lower),condense_values(model, raw_upper)] @@ -373,8 +385,8 @@ function _maximize_M( MOI.GreaterThan(set_values[1])) le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_values[2])) - raw_lower = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) - raw_upper = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + raw_lower = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) + raw_upper = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) (raw_lower === nothing || raw_upper === nothing) && return nothing return [condense_values(model, raw_lower),condense_values(model, raw_upper)] @@ -393,7 +405,7 @@ function _maximize_M( for i in 1:objective.set.dimension le_obj = JuMP.ScalarConstraint( objective.func[i], MOI.LessThan(zero(val_type))) - raw = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + raw = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) raw === nothing && return nothing push!(results, condense_values(model, raw)) end @@ -413,7 +425,7 @@ function _maximize_M( for i in 1:objective.set.dimension ge_obj = JuMP.ScalarConstraint( objective.func[i], MOI.GreaterThan(zero(val_type))) - raw = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) + raw = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) raw === nothing && return nothing push!(results, condense_values(model, raw)) end @@ -435,8 +447,8 @@ function _maximize_M( objective.func[i], MOI.GreaterThan(zero(val_type))) le_obj = JuMP.ScalarConstraint( objective.func[i], MOI.LessThan(zero(val_type))) - raw_ge = _raw_M(sub,_prepare_objectives(model, ge_obj, sub),method) - raw_le = _raw_M(sub,_prepare_objectives(model, le_obj, sub),method) + raw_ge = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) + raw_le = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) (raw_ge === nothing || raw_le === nothing) && return nothing push!(results, condense_values(model, max.(raw_ge, raw_le))) @@ -452,8 +464,14 @@ function _maximize_M( "has not been implemented for MBM subproblems\nF: $(F)") end -# Create a submodel for a disjunct's feasible region. Returns -# GDPSubmodel. Extensions may override for custom construction. +""" + create_submodel(model, constraints, method::_MBM) + +Build a `GDPSubmodel` representing a disjunct's feasible region for +MBM subproblem solves. Copies the model's decision variables and adds +the given disjunct constraints. Submodels are cached in `method.store` +by indicator. +""" function create_submodel( model::JuMP.AbstractModel, constraints::Vector{<:DisjunctConstraintRef}, diff --git a/src/variables.jl b/src/variables.jl index f96cfe02..e1d2d2b3 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -424,8 +424,7 @@ function _interrogate_variables(interrogator::Function, nlp::JuMP.GenericNonline for arg in nlp.args _interrogate_variables(interrogator, arg) end - # TODO avoid recursion. See InfiniteOpt.jl for alternate method that avoids stackoverflow errors with deeply nested expressions: - # https://github.com/infiniteopt/InfiniteOpt.jl/blob/cb6dd6ae40fe0144b1dd75da0739ea6e305d5357/src/expressions.jl#L520-L534 + return end diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index f57e406b..b061473c 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -84,7 +84,7 @@ function test__replace_variables_in_constraint() "String", new_vars) end -function test__prepare_objectives() +function test_prepare_objectives() model = Model() sub_model = Model() @@ -99,14 +99,14 @@ function test__prepare_objectives() collect(keys(new_vars)), new_vars) # LessThan: max(f - upper) = max(x[1] - 1) - objs_le = DP._prepare_objectives( + objs_le = DP.prepare_objectives( model, constraint_object(lessthan), sub) @test length(objs_le) == 1 @test objs_le[1] == JuMP.@expression(sub_model, new_vars[x[1]][1] - 1) # GreaterThan: max(lower - f) = max(1 - x[2]) - objs_ge = DP._prepare_objectives( + objs_ge = DP.prepare_objectives( model, constraint_object(greaterthan), sub) @test length(objs_ge) == 1 @test objs_ge[1] == JuMP.@expression(sub_model, @@ -134,21 +134,21 @@ function test_raw_M() DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub = DP.create_submodel(model, DisjunctConstraintRef[con2], mbm) - objs = DP._prepare_objectives(model, + objs = DP.prepare_objectives(model, constraint_object(con), sub) raw = DP._raw_M(sub, objs, mbm) @test DP.condense_values(model, raw) == 0.0 set_upper_bound(x, 1) sub2 = DP.create_submodel(model, DisjunctConstraintRef[con], mbm) - objs2 = DP._prepare_objectives(model, + objs2 = DP.prepare_objectives(model, constraint_object(con2), sub2) raw = DP._raw_M(sub2, objs2, mbm) @test DP.condense_values(model, raw) == 15 set_integer(y) @constraint(model, con3, y*x == 15, Disjunct(Y[1])) - objs3 = DP._prepare_objectives(model, + objs3 = DP.prepare_objectives(model, constraint_object(con2), sub2) raw = DP._raw_M(sub2, objs3, mbm) @test DP.condense_values(model, raw) == 15 @@ -158,7 +158,7 @@ function test_raw_M() DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub3 = DP.create_submodel(model, DisjunctConstraintRef[con], mbm2) - objs4 = DP._prepare_objectives(model, + objs4 = DP.prepare_objectives(model, constraint_object(con2), sub3) raw = DP._raw_M(sub3, objs4, mbm2) @test DP.condense_values(model, raw) == 10 @@ -168,7 +168,7 @@ function test_raw_M() DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub4 = DP.create_submodel(model, DisjunctConstraintRef[con2], mbm3) - objs5 = DP._prepare_objectives(model, + objs5 = DP.prepare_objectives(model, constraint_object(con2), sub4) @test DP._raw_M(sub4, objs5, mbm3) == nothing @@ -179,7 +179,7 @@ function test_raw_M() sub5 = DP.create_submodel(model, DisjunctConstraintRef[truly_infeasible], mbm4) - objs6 = DP._prepare_objectives(model, + objs6 = DP.prepare_objectives(model, constraint_object(con), sub5) @test DP._raw_M(sub5, objs6, mbm4) == nothing @@ -195,7 +195,7 @@ function test_raw_M() mbm_ub = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub_ub = DP.create_submodel(model_ub, DisjunctConstraintRef[ub_con1], mbm_ub) - objs_ub = DP._prepare_objectives(model_ub, + objs_ub = DP.prepare_objectives(model_ub, constraint_object(ub_con2), sub_ub) raw_ub = DP._raw_M(sub_ub, objs_ub, mbm_ub) @test raw_ub == [mbm_ub.default_M] @@ -825,7 +825,7 @@ end test__var_ref_type_numeric_map() test__replace_variables_quad_numeric_map() test__replace_variables_in_constraint() - test__prepare_objectives() + test_prepare_objectives() test_raw_M() test_maximize_M() test_reformulate_disjunct_constraint() From 387257dcb2a3d5a35add6fd1585d25ec493afcd7 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Mon, 23 Mar 2026 09:34:01 -0400 Subject: [PATCH 07/13] . --- src/datatypes.jl | 6 +++--- src/mbm.jl | 2 +- src/variables.jl | 3 ++- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index c75e54f1..5bfe7658 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -439,17 +439,17 @@ mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationM end """ - cutting_planes{O} <: AbstractReformulationMethod + cutting_planes{O,T} <: AbstractReformulationMethod A type for using the cutting planes approach for disjunctive constraints. **Fields** - `optimizer::O`: Optimizer to use when solving mini-models (required). - `max_iter::Int`: Number of iterations (default = `3`). -- `seperation_tolerance::Float64`: Tolerance for the separation problem (default = `1e-6`). +- `seperation_tolerance::T`: Tolerance for the separation problem (default = `1e-6`). - `final_reform_method::AbstractReformulationMethod`: Final reformulation method to use after cutting planes (default = `BigM()`). -- `M_value::Float64`: Big-M value to use in the final reformulation (default = `1e9`). +- `M_value::T`: Big-M value to use in the final reformulation (default = `1e9`). """ struct cutting_planes{O, T} <: AbstractReformulationMethod optimizer::O; diff --git a/src/mbm.jl b/src/mbm.jl index 025c1427..25c62e5a 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -9,7 +9,7 @@ _is_all_zeros(::Any) = false ################################################################################ # CONSTRAINT, DISJUNCTION, DISJUNCT REFORMULATION ################################################################################ - +#Reformulates the disjunction using multiple big-M values function reformulate_disjunction( model::JuMP.AbstractModel, disj::Disjunction, diff --git a/src/variables.jl b/src/variables.jl index e1d2d2b3..f96cfe02 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -424,7 +424,8 @@ function _interrogate_variables(interrogator::Function, nlp::JuMP.GenericNonline for arg in nlp.args _interrogate_variables(interrogator, arg) end - + # TODO avoid recursion. See InfiniteOpt.jl for alternate method that avoids stackoverflow errors with deeply nested expressions: + # https://github.com/infiniteopt/InfiniteOpt.jl/blob/cb6dd6ae40fe0144b1dd75da0739ea6e305d5357/src/expressions.jl#L520-L534 return end From c52d97442c4f2f7b78f10078cfaea50e2c0935f6 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Tue, 31 Mar 2026 11:20:33 -0400 Subject: [PATCH 08/13] . --- README.md | 2 +- ext/InfiniteDisjunctiveProgramming.jl | 4 +- src/cuttingplanes.jl | 12 +++--- src/datatypes.jl | 12 +++--- src/mbm.jl | 40 ++++++++++--------- test/constraints/cuttingplanes.jl | 22 +++++----- test/constraints/mbm.jl | 8 ++-- .../InfiniteDisjunctiveProgramming.jl | 2 +- test/solve.jl | 2 +- 9 files changed, 53 insertions(+), 51 deletions(-) diff --git a/README.md b/README.md index fce0aa0f..543d1f7a 100644 --- a/README.md +++ b/README.md @@ -184,7 +184,7 @@ The following reformulation methods are currently supported: All variables must be included in exactly one partition. For manual partitioning, ensure each variable appears in exactly one group. For automatic partitioning, variables are divided as evenly as possible among the specified number of partitions. -6. [Cutting Planes](https://pubsonline.informs.org/doi/10.1287/ijoc.2015.0669): This method iteratively generates cutting planes using a separation problem and a relaxed Big-M formulation, then applies a final reformulation method. The `cutting_planes` struct is created with the following arguments: +6. [Cutting Planes](https://pubsonline.informs.org/doi/10.1287/ijoc.2015.0669): This method iteratively generates cutting planes using a separation problem and a relaxed Big-M formulation, then applies a final reformulation method. The `CuttingPlanes` struct is created with the following arguments: - `optimizer`: Optimizer to use when solving the separation and relaxed Big-M subproblems. This is a required value. - `max_iter`: Maximum number of cutting plane iterations. Default: `3`. diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index fad26baa..426f7ac6 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -178,8 +178,8 @@ function DP.reformulate_model(::InfiniteOpt.InfiniteModel, ::DP.MBM) "Please use `BigM`, `Hull`, `Indicator`, or `PSplit` instead.") end -function DP.reformulate_model(::InfiniteOpt.InfiniteModel, ::DP.cutting_planes) - error("The `cutting_planes` method is not supported for `InfiniteModel`." * +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.") end diff --git a/src/cuttingplanes.jl b/src/cuttingplanes.jl index 6d2c387f..2020a3cf 100644 --- a/src/cuttingplanes.jl +++ b/src/cuttingplanes.jl @@ -1,6 +1,6 @@ function reformulate_model( model::JuMP.AbstractModel, - method::cutting_planes + method::CuttingPlanes ) _clear_reformulations(model) var_type = JuMP.variable_ref_type(model) @@ -39,7 +39,7 @@ function reformulate_model( rBM_sol = _solve_rBM(rBM) SEP_sol = _solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map, rBM_to_SEP_map) sep_obj = objective_value(SEP) - _cutting_planes(model, rBM, main_to_rBM_map, + _CuttingPlanes(model, rBM, main_to_rBM_map, main_to_SEP_map, rBM_sol, SEP_sol ) i += 1 @@ -90,7 +90,7 @@ function _solve_SEP( return sol end -function _cutting_planes( +function _CuttingPlanes( model::M, rBM::M, main_to_rBM_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef}, @@ -123,7 +123,7 @@ end # ERROR MESSAGES ################################################################################ -function reformulate_model(::M, ::cutting_planes) where {M} +function reformulate_model(::M, ::CuttingPlanes) where {M} error("reformulate_model not implemented for model type `$(M)`.") end @@ -139,8 +139,8 @@ function _solve_SEP(::M, ::N, ::H, ::S, ::R) where {M, N, H, S, R} rBM_to_SEP_map: `$(R)`.") end -function _cutting_planes(::M, ::N, ::H, ::S, ::R, ::T) where {M, N, H, S, R, T} - error("_cutting_planes not implemented for argument types: \n +function _CuttingPlanes(::M, ::N, ::H, ::S, ::R, ::T) where {M, N, H, S, R, T} + error("_CuttingPlanes not implemented for argument types: \n model: `$(M)`, rBM: `$(N)`,\n main_to_rBM_map: `$(H)`, main_to_SEP_map: `$(S)`,\n diff --git a/src/datatypes.jl b/src/datatypes.jl index 5bfe7658..fe28697e 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -390,11 +390,11 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho optimizer::O M::Dict{LogicalVariableRef{M}, Any} default_M::T - conlvref::Vector{LogicalVariableRef{M}} + subproblem_indicators::Vector{LogicalVariableRef{M}} deactivated::Set{LogicalVariableRef{M}} - # Stored submodels: indicator => GDPSubmodel. + # Cached submodels: indicator => GDPSubmodel. # Typed Any so extensions can store different types. - store::Dict{LogicalVariableRef{M}, Any} + model_cache::Dict{LogicalVariableRef{M}, Any} function _MBM(method::MBM{O, T}, model::M) where {O, T, M <: JuMP.AbstractModel} new{O, T, M}( @@ -439,7 +439,7 @@ mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationM end """ - cutting_planes{O,T} <: AbstractReformulationMethod + CuttingPlanes{O,T} <: AbstractReformulationMethod A type for using the cutting planes approach for disjunctive constraints. @@ -451,13 +451,13 @@ A type for using the cutting planes approach for disjunctive constraints. method to use after cutting planes (default = `BigM()`). - `M_value::T`: Big-M value to use in the final reformulation (default = `1e9`). """ -struct cutting_planes{O, T} <: AbstractReformulationMethod +struct CuttingPlanes{O, T} <: AbstractReformulationMethod optimizer::O; max_iter::Int seperation_tolerance::T final_reform_method::AbstractReformulationMethod M_value::T - function cutting_planes( + function CuttingPlanes( optimizer::O; max_iter::Int = 3, seperation_tolerance::T = 1e-6, diff --git a/src/mbm.jl b/src/mbm.jl index 25c62e5a..762f7695 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -19,7 +19,7 @@ function reformulate_disjunction( disjunct_cons = Dict{LogicalVariableRef, Vector{JuMP.AbstractConstraint}}() for d in disj.indicators d in mbm.deactivated && continue - mbm.conlvref = filter( + mbm.subproblem_indicators = filter( x -> x != d && !(x in mbm.deactivated), disj.indicators) disjunct_cons[d] = Vector{JuMP.AbstractConstraint}() _reformulate_disjunct(model, disjunct_cons[d], d, mbm) @@ -40,13 +40,14 @@ end function _reformulate_disjunct( model::JuMP.AbstractModel, ref_cons::Vector{JuMP.AbstractConstraint}, - lvref::LogicalVariableRef, method::_MBM + lvref::LogicalVariableRef, + method::_MBM ) !haskey(_indicator_to_constraints(model), lvref) && return # Filter out deactivated disjuncts from binary variable mapping in # the event we've identified some infeasible disjuncts already - active_conlvref = filter(d -> !(d in method.deactivated), method.conlvref) - bconref = Dict(d => binary_variable(d) for d in active_conlvref) + active_subproblem_indicators = filter(d -> !(d in method.deactivated), method.subproblem_indicators) + bconref = Dict(d => binary_variable(d) for d in active_subproblem_indicators) constraints = _indicator_to_constraints(model)[lvref] filtered_constraints = [ @@ -56,7 +57,7 @@ function _reformulate_disjunct( for cref in filtered_constraints empty!(method.M) - for d in method.conlvref + for d in method.subproblem_indicators # Skip already-deactivated disjuncts d in method.deactivated && continue @@ -93,7 +94,8 @@ function _reformulate_disjunct( end function reformulate_disjunct_constraint( - model::JuMP.AbstractModel, con::Disjunction, + model::JuMP.AbstractModel, + con::Disjunction, bconref::Union{ Dict{<:LogicalVariableRef, <:JuMP.AbstractVariableRef}, Dict{<:LogicalVariableRef, <:JuMP.GenericAffExpr} @@ -309,7 +311,7 @@ function _raw_M( end """ - condense_values(model, vals::AbstractVector) + aggregate_M_values(model, vals::AbstractVector) Reduce a vector of raw M values from `_raw_M` to the final form used during constraint reformulation. The base method returns the single @@ -317,11 +319,11 @@ element from a length-1 vector. Extensions (e.g., InfiniteOpt) override to aggregate across multiple support points (e.g., interpolating K per-support M values into a parameter function). """ -function condense_values( +function aggregate_M_values( ::JuMP.AbstractModel, vals::AbstractVector ) - return vals[1] + return only(vals) end # Dispatch over constraint types to compute M values. Scalar @@ -336,7 +338,7 @@ function _maximize_M( objectives = prepare_objectives(model, objective, sub) raw = _raw_M(sub, objectives, method) raw === nothing && return nothing - return condense_values(model, raw) + return aggregate_M_values(model, raw) end # Helper: get or create the submodel for a set of constraints. @@ -347,11 +349,11 @@ function _get_submodel( ) indicator = _constraint_to_indicator( model)[first(constraints)] - if !haskey(method.store, indicator) - method.store[indicator] = create_submodel( + if !haskey(method.model_cache, indicator) + method.model_cache[indicator] = create_submodel( model, constraints, method) end - return method.store[indicator] + return method.model_cache[indicator] end # EqualTo: solve both GreaterThan and LessThan directions, finalize each. @@ -369,7 +371,7 @@ function _maximize_M( raw_upper = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) (raw_lower === nothing || raw_upper === nothing) && return nothing - return [condense_values(model, raw_lower),condense_values(model, raw_upper)] + return [aggregate_M_values(model, raw_lower),aggregate_M_values(model, raw_upper)] end # Interval: solve both lower and upper bound directions, finalize each. @@ -389,7 +391,7 @@ function _maximize_M( raw_upper = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) (raw_lower === nothing || raw_upper === nothing) && return nothing - return [condense_values(model, raw_lower),condense_values(model, raw_upper)] + return [aggregate_M_values(model, raw_lower),aggregate_M_values(model, raw_upper)] end # Nonpositives: per-row LessThan solves for each dimension of the vector. @@ -407,7 +409,7 @@ function _maximize_M( objective.func[i], MOI.LessThan(zero(val_type))) raw = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) raw === nothing && return nothing - push!(results, condense_values(model, raw)) + push!(results, aggregate_M_values(model, raw)) end return results end @@ -427,7 +429,7 @@ function _maximize_M( objective.func[i], MOI.GreaterThan(zero(val_type))) raw = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) raw === nothing && return nothing - push!(results, condense_values(model, raw)) + push!(results, aggregate_M_values(model, raw)) end return results end @@ -451,7 +453,7 @@ function _maximize_M( raw_le = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) (raw_ge === nothing || raw_le === nothing) && return nothing - push!(results, condense_values(model, max.(raw_ge, raw_le))) + push!(results, aggregate_M_values(model, max.(raw_ge, raw_le))) end return results end @@ -469,7 +471,7 @@ end Build a `GDPSubmodel` representing a disjunct's feasible region for MBM subproblem solves. Copies the model's decision variables and adds -the given disjunct constraints. Submodels are cached in `method.store` +the given disjunct constraints. Submodels are cached in `method.model_cache` by indicator. """ function create_submodel( diff --git a/test/constraints/cuttingplanes.jl b/test/constraints/cuttingplanes.jl index 3c5239e9..1264e79a 100644 --- a/test/constraints/cuttingplanes.jl +++ b/test/constraints/cuttingplanes.jl @@ -1,14 +1,14 @@ using HiGHS -function test_cutting_planes_datatype() - method = cutting_planes(HiGHS.Optimizer) +function test_CuttingPlanes_datatype() + method = CuttingPlanes(HiGHS.Optimizer) @test method.optimizer == HiGHS.Optimizer @test method.max_iter == 3 @test method.seperation_tolerance == 1e-6 @test method.final_reform_method isa BigM @test method.M_value == 1e9 - method = cutting_planes(HiGHS.Optimizer;max_iter=10, + method = CuttingPlanes(HiGHS.Optimizer;max_iter=10, seperation_tolerance=1e-4, final_reform_method=Indicator(), M_value=1e6 ) @test method.max_iter == 10 @@ -43,7 +43,7 @@ function test_solve_SEP() @disjunction(model, [Y[1], Y[2]]) @objective(model, Max, x) var_type = JuMP.variable_ref_type(model) - method = cutting_planes(HiGHS.Optimizer) + method = CuttingPlanes(HiGHS.Optimizer) obj = objective_function(model) sense = objective_sense(model) SEP, sep_ref_map, _ = DP.copy_gdp_model(model) @@ -78,7 +78,7 @@ function test_solve_SEP() ) end -function test_cutting_planes() +function test_CuttingPlanes() model = GDPModel() @variable(model, 0 <= x <= 100) @variable(model, Y[1:2], Logical) @@ -87,7 +87,7 @@ function test_cutting_planes() @disjunction(model, [Y[1], Y[2]]) @objective(model, Max, x) var_type = JuMP.variable_ref_type(model) - method = cutting_planes(HiGHS.Optimizer) + method = CuttingPlanes(HiGHS.Optimizer) obj = objective_function(model) sense = objective_sense(model) SEP, sep_ref_map, _ = DP.copy_gdp_model(model) @@ -114,7 +114,7 @@ function test_cutting_planes() end rBM_sol = DP._solve_rBM(rBM) SEP_sol = DP._solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map, rBM_to_SEP_map) - DP._cutting_planes(model, rBM, main_to_rBM_map, main_to_SEP_map, rBM_sol, SEP_sol) + DP._CuttingPlanes(model, rBM, main_to_rBM_map, main_to_SEP_map, rBM_sol, SEP_sol) rBM_sol = DP._solve_rBM(rBM) SEP_sol = DP._solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map, rBM_to_SEP_map) @@ -122,7 +122,7 @@ function test_cutting_planes() @test rBM_sol[main_to_rBM_map[x]] ≈ 4.0 @test SEP_sol[rBM_to_SEP_map[main_to_rBM_map[x]]] ≈ 4.0 atol=1e-3 - @test_throws ErrorException DP._cutting_planes(model, rBM, main_to_rBM_map, + @test_throws ErrorException DP._CuttingPlanes(model, rBM, main_to_rBM_map, main_to_SEP_map, rBM_sol, "not a dict" ) end @@ -136,7 +136,7 @@ function test_reformulate_model() @disjunction(model, [Y[1], Y[2]]) @objective(model, Max, x[1] + x[2]) - method = cutting_planes(HiGHS.Optimizer) + method = CuttingPlanes(HiGHS.Optimizer) DP.reformulate_model(model, method) num_con = length( JuMP.all_constraints(model; include_variable_in_set_constraints = false) @@ -147,9 +147,9 @@ end @testset "Cutting Planes" begin - test_cutting_planes_datatype() + test_CuttingPlanes_datatype() test_solve_rBM() test_solve_SEP() - test_cutting_planes() + test_CuttingPlanes() test_reformulate_model() end diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index b061473c..ee4d7e0c 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -137,21 +137,21 @@ function test_raw_M() objs = DP.prepare_objectives(model, constraint_object(con), sub) raw = DP._raw_M(sub, objs, mbm) - @test DP.condense_values(model, raw) == 0.0 + @test DP.aggregate_M_values(model, raw) == 0.0 set_upper_bound(x, 1) sub2 = DP.create_submodel(model, DisjunctConstraintRef[con], mbm) objs2 = DP.prepare_objectives(model, constraint_object(con2), sub2) raw = DP._raw_M(sub2, objs2, mbm) - @test DP.condense_values(model, raw) == 15 + @test DP.aggregate_M_values(model, raw) == 15 set_integer(y) @constraint(model, con3, y*x == 15, Disjunct(Y[1])) objs3 = DP.prepare_objectives(model, constraint_object(con2), sub2) raw = DP._raw_M(sub2, objs3, mbm) - @test DP.condense_values(model, raw) == 15 + @test DP.aggregate_M_values(model, raw) == 15 # Fresh _MBM after changing bounds JuMP.fix(y, 5; force=true) mbm2 = DP._MBM( @@ -161,7 +161,7 @@ function test_raw_M() objs4 = DP.prepare_objectives(model, constraint_object(con2), sub3) raw = DP._raw_M(sub3, objs4, mbm2) - @test DP.condense_values(model, raw) == 10 + @test DP.aggregate_M_values(model, raw) == 10 # Infeasible region → nothing delete_lower_bound(x) mbm3 = DP._MBM( diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl index b206d8e0..8ba1e5c2 100644 --- a/test/extensions/InfiniteDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteDisjunctiveProgramming.jl @@ -339,7 +339,7 @@ end function test_unsupported_methods_error() model = InfiniteGDPModel(HiGHS.Optimizer) @test_throws ErrorException DP.reformulate_model(model, MBM(HiGHS.Optimizer)) - @test_throws ErrorException DP.reformulate_model(model, cutting_planes(HiGHS.Optimizer)) + @test_throws ErrorException DP.reformulate_model(model, CuttingPlanes(HiGHS.Optimizer)) end function test_methods() diff --git a/test/solve.jl b/test/solve.jl index b8b8db32..c3ca9441 100644 --- a/test/solve.jl +++ b/test/solve.jl @@ -54,7 +54,7 @@ function test_linear_gdp_example(m, use_complements = false) @test !value(W[1]) @test !value(W[2]) - @test optimize!(m, gdp_method = cutting_planes(HiGHS.Optimizer)) isa Nothing + @test optimize!(m, gdp_method = CuttingPlanes(HiGHS.Optimizer)) isa Nothing @test termination_status(m) == MOI.OPTIMAL @test objective_value(m) ≈ 11 atol=1e-3 @test value.(x) ≈ [9,2] atol=1e-3 From 527bd541603f9cbfd2bb91bbba3e1b5f9093ed3e Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Tue, 31 Mar 2026 13:48:48 -0400 Subject: [PATCH 09/13] renaming and spaces --- src/mbm.jl | 38 +++++++++++++++++++++----------------- 1 file changed, 21 insertions(+), 17 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 762f7695..6ef2369d 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -46,7 +46,10 @@ function _reformulate_disjunct( !haskey(_indicator_to_constraints(model), lvref) && return # Filter out deactivated disjuncts from binary variable mapping in # the event we've identified some infeasible disjuncts already - active_subproblem_indicators = filter(d -> !(d in method.deactivated), method.subproblem_indicators) + active_subproblem_indicators = filter( + d -> !(d in method.deactivated), + method.subproblem_indicators + ) bconref = Dict(d => binary_variable(d) for d in active_subproblem_indicators) constraints = _indicator_to_constraints(model)[lvref] @@ -83,7 +86,8 @@ function _reformulate_disjunct( # Check if all M values are zero for that constraint. If so, it # should be enforced globally (no reformulation with binaries). if !isempty(method.M) && all( - _is_all_zeros(method.M[d]) for d in keys(method.M)) + _is_all_zeros(method.M[d]) for d in keys(method.M) + ) push!(ref_cons, con) else append!(ref_cons, @@ -290,11 +294,6 @@ function _raw_M( M_vals = typeof(method.default_M)[] for obj_expr in objectives JuMP.@objective(sub.model, Max, obj_expr) - # Clear start values before each solve to prevent NaN - # residue from a previous non-feasible solve - for v in JuMP.all_variables(sub.model) - JuMP.set_start_value(v, nothing) - end JuMP.optimize!(sub.model) if JuMP.termination_status(sub.model) == _MOI.INFEASIBLE return nothing @@ -367,11 +366,13 @@ function _maximize_M( set_value = objective.set.value ge_obj = JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)) le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)) - raw_lower = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) - raw_upper = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) + raw_lower = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) + raw_upper = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) (raw_lower === nothing || raw_upper === nothing) && return nothing - return [aggregate_M_values(model, raw_lower),aggregate_M_values(model, raw_upper)] + return [aggregate_M_values(model, raw_lower), + aggregate_M_values(model, raw_upper) + ] end # Interval: solve both lower and upper bound directions, finalize each. @@ -387,11 +388,14 @@ function _maximize_M( MOI.GreaterThan(set_values[1])) le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_values[2])) - raw_lower = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) - raw_upper = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) + raw_lower = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) + raw_upper = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) (raw_lower === nothing || raw_upper === nothing) && return nothing - return [aggregate_M_values(model, raw_lower),aggregate_M_values(model, raw_upper)] + return [ + aggregate_M_values(model, raw_lower), + aggregate_M_values(model, raw_upper) + ] end # Nonpositives: per-row LessThan solves for each dimension of the vector. @@ -407,7 +411,7 @@ function _maximize_M( for i in 1:objective.set.dimension le_obj = JuMP.ScalarConstraint( objective.func[i], MOI.LessThan(zero(val_type))) - raw = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) + raw = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) raw === nothing && return nothing push!(results, aggregate_M_values(model, raw)) end @@ -427,7 +431,7 @@ function _maximize_M( for i in 1:objective.set.dimension ge_obj = JuMP.ScalarConstraint( objective.func[i], MOI.GreaterThan(zero(val_type))) - raw = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) + raw = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) raw === nothing && return nothing push!(results, aggregate_M_values(model, raw)) end @@ -449,8 +453,8 @@ function _maximize_M( objective.func[i], MOI.GreaterThan(zero(val_type))) le_obj = JuMP.ScalarConstraint( objective.func[i], MOI.LessThan(zero(val_type))) - raw_ge = _raw_M(sub,prepare_objectives(model, ge_obj, sub),method) - raw_le = _raw_M(sub,prepare_objectives(model, le_obj, sub),method) + raw_ge = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) + raw_le = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) (raw_ge === nothing || raw_le === nothing) && return nothing push!(results, aggregate_M_values(model, max.(raw_ge, raw_le))) From 603c5523a7fdd363dd51f0a896fa881e75edf5cb Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Tue, 31 Mar 2026 15:44:57 -0400 Subject: [PATCH 10/13] . --- src/mbm.jl | 24 ++++++++++++------------ test/constraints/mbm.jl | 20 ++++++++++---------- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 6ef2369d..efb223ae 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -295,15 +295,15 @@ function _raw_M( for obj_expr in objectives JuMP.@objective(sub.model, Max, obj_expr) JuMP.optimize!(sub.model) - if JuMP.termination_status(sub.model) == _MOI.INFEASIBLE - return nothing - elseif !JuMP.is_solved_and_feasible(sub.model) - push!(M_vals, method.default_M) - else + if JuMP.is_solved_and_feasible(sub.model) push!(M_vals, max( JuMP.objective_value(sub.model), zero(method.default_M)) ) + elseif JuMP.termination_status(sub.model) == _MOI.INFEASIBLE + return nothing + else + push!(M_vals, method.default_M) end end return M_vals @@ -349,7 +349,7 @@ function _get_submodel( indicator = _constraint_to_indicator( model)[first(constraints)] if !haskey(method.model_cache, indicator) - method.model_cache[indicator] = create_submodel( + method.model_cache[indicator] = copy_model_with_constraints( model, constraints, method) end return method.model_cache[indicator] @@ -471,14 +471,14 @@ function _maximize_M( end """ - create_submodel(model, constraints, method::_MBM) + copy_model_with_constraints(model, constraints, method) -Build a `GDPSubmodel` representing a disjunct's feasible region for -MBM subproblem solves. Copies the model's decision variables and adds -the given disjunct constraints. Submodels are cached in `method.model_cache` -by indicator. +Build a `GDPSubmodel` with disjunct constraints passed. +This builds a model seperate from the original model with copied constraints +and variables, and maps between the original model's variables and the +submodel's variables. """ -function create_submodel( +function copy_model_with_constraints( model::JuMP.AbstractModel, constraints::Vector{<:DisjunctConstraintRef}, method::_MBM diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index ee4d7e0c..b6bffe82 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -132,14 +132,14 @@ function test_raw_M() [Y[1], Y[2], Y[3], Y[4], Y[5]]) mbm = DP._MBM( DP.MBM(HiGHS.Optimizer), JuMP.Model()) - sub = DP.create_submodel(model, + sub = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con2], mbm) objs = DP.prepare_objectives(model, constraint_object(con), sub) raw = DP._raw_M(sub, objs, mbm) @test DP.aggregate_M_values(model, raw) == 0.0 set_upper_bound(x, 1) - sub2 = DP.create_submodel(model, + sub2 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con], mbm) objs2 = DP.prepare_objectives(model, constraint_object(con2), sub2) @@ -156,7 +156,7 @@ function test_raw_M() JuMP.fix(y, 5; force=true) mbm2 = DP._MBM( DP.MBM(HiGHS.Optimizer), JuMP.Model()) - sub3 = DP.create_submodel(model, + sub3 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con], mbm2) objs4 = DP.prepare_objectives(model, constraint_object(con2), sub3) @@ -166,7 +166,7 @@ function test_raw_M() delete_lower_bound(x) mbm3 = DP._MBM( DP.MBM(HiGHS.Optimizer), JuMP.Model()) - sub4 = DP.create_submodel(model, + sub4 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con2], mbm3) objs5 = DP.prepare_objectives(model, constraint_object(con2), sub4) @@ -176,7 +176,7 @@ function test_raw_M() set_upper_bound(x, 1) mbm4 = DP._MBM( DP.MBM(HiGHS.Optimizer), JuMP.Model()) - sub5 = DP.create_submodel(model, + sub5 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[truly_infeasible], mbm4) objs6 = DP.prepare_objectives(model, @@ -193,7 +193,7 @@ function test_raw_M() @constraint(model_ub, ub_con2, xu >= 5, Disjunct(Yu[2])) @disjunction(model_ub, [Yu[1], Yu[2]]) mbm_ub = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) - sub_ub = DP.create_submodel(model_ub, + sub_ub = DP.copy_model_with_constraints(model_ub, DisjunctConstraintRef[ub_con1], mbm_ub) objs_ub = DP.prepare_objectives(model_ub, constraint_object(ub_con2), sub_ub) @@ -733,8 +733,8 @@ function test_variable_copy() @test owner_model(w_copy) === target_model end -# Test _create_submodel (combines _copy_model + variable_copy + constraint handling) -function test__create_submodel() +# Test _copy_model_with_constraints (combines _copy_model + variable_copy + constraint handling) +function test__copy_model_with_constraints() model = GDPModel() @variable(model, 0 <= x <= 10) @variable(model, 0 <= y <= 5) @@ -748,7 +748,7 @@ function test__create_submodel() DP._indicator_to_constraints(model)[Y[1]] ) - sub = DP.create_submodel(model, constraints, mbm) + sub = DP.copy_model_with_constraints(model, constraints, mbm) # Check submodel struct @test sub isa DP.GDPSubmodel @@ -819,7 +819,7 @@ end test__make_variable_object() test_create_variable() test_variable_copy() - test__create_submodel() + test__copy_model_with_constraints() test_get_variable_info() test_mbm() test__var_ref_type_numeric_map() From 82dc6aea36db6e20c3521cdadab2bddb16e8a184 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Tue, 31 Mar 2026 15:55:50 -0400 Subject: [PATCH 11/13] removal of deactivation logic --- src/datatypes.jl | 2 -- src/mbm.jl | 53 +++++++++++++---------------------------- test/constraints/mbm.jl | 29 +++------------------- 3 files changed, 20 insertions(+), 64 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index fe28697e..d948b7cf 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -391,7 +391,6 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho M::Dict{LogicalVariableRef{M}, Any} default_M::T subproblem_indicators::Vector{LogicalVariableRef{M}} - deactivated::Set{LogicalVariableRef{M}} # Cached submodels: indicator => GDPSubmodel. # Typed Any so extensions can store different types. model_cache::Dict{LogicalVariableRef{M}, Any} @@ -402,7 +401,6 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho Dict{LogicalVariableRef{M}, Any}(), method.default_M, Vector{LogicalVariableRef{M}}(), - Set{LogicalVariableRef{M}}(), Dict{LogicalVariableRef{M}, Any}() ) end diff --git a/src/mbm.jl b/src/mbm.jl index efb223ae..7e9c49f3 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -16,21 +16,11 @@ function reformulate_disjunction( method::MBM ) mbm = _MBM(method, model) - disjunct_cons = Dict{LogicalVariableRef, Vector{JuMP.AbstractConstraint}}() - for d in disj.indicators - d in mbm.deactivated && continue - mbm.subproblem_indicators = filter( - x -> x != d && !(x in mbm.deactivated), disj.indicators) - disjunct_cons[d] = Vector{JuMP.AbstractConstraint}() - _reformulate_disjunct(model, disjunct_cons[d], d, mbm) - end - # Collect constraints from non-deactivated disjuncts. It needs to be - # in a separate loop because disjuncts are only deactivated by looking - # at reforming other disjuncts (subproblem infeasibility). ref_cons = Vector{JuMP.AbstractConstraint}() for d in disj.indicators - d in mbm.deactivated && continue - haskey(disjunct_cons, d) && append!(ref_cons, disjunct_cons[d]) + mbm.subproblem_indicators = filter( + x -> x != d, disj.indicators) + _reformulate_disjunct(model, ref_cons, d, mbm) end return ref_cons end @@ -44,13 +34,8 @@ function _reformulate_disjunct( method::_MBM ) !haskey(_indicator_to_constraints(model), lvref) && return - # Filter out deactivated disjuncts from binary variable mapping in - # the event we've identified some infeasible disjuncts already - active_subproblem_indicators = filter( - d -> !(d in method.deactivated), - method.subproblem_indicators - ) - bconref = Dict(d => binary_variable(d) for d in active_subproblem_indicators) + bconref = Dict( + d => binary_variable(d) for d in method.subproblem_indicators) constraints = _indicator_to_constraints(model)[lvref] filtered_constraints = [ @@ -61,37 +46,33 @@ function _reformulate_disjunct( empty!(method.M) for d in method.subproblem_indicators - # Skip already-deactivated disjuncts - d in method.deactivated && continue - d_constraints = _indicator_to_constraints(model)[d] disjunct_constraints = [ c for c in d_constraints if c isa DisjunctConstraintRef] if !isempty(disjunct_constraints) - M_result = _maximize_M(model, JuMP.constraint_object(cref), + M_result = _maximize_M(model, + JuMP.constraint_object(cref), disjunct_constraints, method) - # Check for infeasibility: disjunct d - # has empty feasible region if M_result === nothing - push!(method.deactivated, d) - @warn "Disjunct $(d) is infeasible, deactivating." - delete!(bconref, d) - else - method.M[d] = M_result + error("Disjunct $(d) has an infeasible feasible region. + Check the disjunct constraints and variable bounds." + ) end + method.M[d] = M_result end end con = JuMP.constraint_object(cref) - # Check if all M values are zero for that constraint. If so, it - # should be enforced globally (no reformulation with binaries). + # Check if all M values are zero for that constraint. + # If so, enforce globally (no reformulation needed). if !isempty(method.M) && all( - _is_all_zeros(method.M[d]) for d in keys(method.M) - ) + _is_all_zeros(method.M[d]) + for d in keys(method.M)) push!(ref_cons, con) else append!(ref_cons, - reformulate_disjunct_constraint(model, con, bconref, method)) + reformulate_disjunct_constraint( + model, con, bconref, method)) end end return ref_cons diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index b6bffe82..a4296e18 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -492,7 +492,7 @@ function test_reformulate_disjunction() # Global constraint has just x, no binary variables @test JuMP.coefficient(func_3, x) == 1.0 - #Test infeasible disjunct detection and deactivation + # Test infeasible disjunct throws error model2 = GDPModel() @variable(model2, 0 <= z <= 1) @variable(model2, W[1:2], Logical) @@ -503,31 +503,8 @@ function test_reformulate_disjunction() disj2 = disjunction(model2, [W[1], W[2]]) method2 = DP.MBM(HiGHS.Optimizer) - ref_cons2 = @test_logs (:warn, r"infeasible, deactivating") begin - reformulate_disjunction(model2, constraint_object(disj2), method2) - end - - @test length(ref_cons2) == 1 - @test ref_cons2[1].set == MOI.LessThan(0.5) - - #Test multiple infeasible disjuncts - model3 = GDPModel() - @variable(model3, 0 <= w <= 1) - @variable(model3, V[1:3], Logical) - @constraint(model3, w >= 5, Disjunct(V[1])) # infeasible - @constraint(model3, w >= 10, Disjunct(V[2])) # infeasible - @constraint(model3, w <= 0.5, Disjunct(V[3])) # feasible - disj3 = disjunction(model3, [V[1], V[2], V[3]]) - - method3 = DP.MBM(HiGHS.Optimizer) - #warn about V[1] and V[2] being infeasible - ref_cons3 = @test_logs (:warn,) (:warn,) begin - reformulate_disjunction(model3, constraint_object(disj3), method3) - end - - # Only V[3]'s constraint should be reformulated - @test length(ref_cons3) == 1 - @test ref_cons3[1].set == MOI.LessThan(0.5) + @test_throws ErrorException reformulate_disjunction( + model2, constraint_object(disj2), method2) model4 = GDPModel() @variable(model4, 0 <= u <= 10) From 9c5235b923f381dbe3620dd9c341f91020756d12 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Wed, 1 Apr 2026 12:08:28 -0400 Subject: [PATCH 12/13] . --- src/datatypes.jl | 4 ++-- src/mbm.jl | 14 +++++++------- test/constraints/mbm.jl | 12 ++++++------ 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index d948b7cf..866dbd2b 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -483,7 +483,7 @@ along with mappings back to the original model's variables. feasible region (constraints and variable bounds). - `dec_vars::Vector{V}`: Ordered decision variables in the submodel, matching the original model's ordering. -- `fwd::Dict{V, Vector{W}}`: Forward map from original +- `fwd_map::Dict{V, Vector{W}}`: Forward map from original model variables to their submodel counterparts. """ struct GDPSubmodel{M <: JuMP.AbstractModel, @@ -491,7 +491,7 @@ struct GDPSubmodel{M <: JuMP.AbstractModel, W <: JuMP.AbstractVariableRef} model::M dec_vars::Vector{V} - fwd::Dict{V, Vector{W}} + fwd_map::Dict{V, Vector{W}} end """ diff --git a/src/mbm.jl b/src/mbm.jl index 7e9c49f3..79839edf 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -243,14 +243,14 @@ end Convert a constraint into objective expressions for M-value maximization. Returns a vector of JuMP expressions to pass to `_raw_M`. The base method -produces a single-element vector by mapping variables through `sub.fwd`. +produces a single-element vector by mapping variables through `sub.fwd_map`. """ function prepare_objectives( ::JuMP.AbstractModel, obj::JuMP.ScalarConstraint{T, S}, sub::GDPSubmodel ) where {T, S <: _MOI.LessThan} - flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd) + flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map) expr = -obj.set.upper +_replace_variables_in_constraint(obj.func, flat_map) return [expr] end @@ -260,7 +260,7 @@ function prepare_objectives( obj::JuMP.ScalarConstraint{T, S}, sub::GDPSubmodel ) where {T, S <: _MOI.GreaterThan} - flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd) + flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map) expr = obj.set.lower -_replace_variables_in_constraint(obj.func, flat_map) return [expr] end @@ -467,16 +467,16 @@ function copy_model_with_constraints( var_type = JuMP.variable_ref_type(model) sub_model = _copy_model(model) dec_vars = collect_all_vars(model) - fwd = Dict{var_type, Vector{var_type}}() + fwd_map = Dict{var_type, Vector{var_type}}() for var in dec_vars copy_var = variable_copy(sub_model, var) - fwd[var] = [copy_var] + fwd_map[var] = [copy_var] end for cref in constraints con = JuMP.constraint_object(cref) - flat_map = Dict(v => ws[1] for (v, ws) in fwd) + flat_map = Dict(v => ws[1] for (v, ws) in fwd_map) expr = _replace_variables_in_constraint( con.func, flat_map) T = one(JuMP.value_type(typeof(sub_model))) @@ -486,7 +486,7 @@ function copy_model_with_constraints( JuMP.set_optimizer(sub_model, method.optimizer) JuMP.set_silent(sub_model) - return GDPSubmodel(sub_model, dec_vars, fwd) + return GDPSubmodel(sub_model, dec_vars, fwd_map) end ################################################################################ diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index a4296e18..ac04f217 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -731,16 +731,16 @@ function test__copy_model_with_constraints() @test sub isa DP.GDPSubmodel @test sub.model isa Model # Check variable mapping exists - @test haskey(sub.fwd, x) - @test haskey(sub.fwd, y) + @test haskey(sub.fwd_map, x) + @test haskey(sub.fwd_map, y) # Check mapped variables (length-1 vectors) - @test length(sub.fwd[x]) == 1 - @test length(sub.fwd[y]) == 1 + @test length(sub.fwd_map[x]) == 1 + @test length(sub.fwd_map[y]) == 1 # Check bounds in submodel - xm = sub.fwd[x][1] - ym = sub.fwd[y][1] + xm = sub.fwd_map[x][1] + ym = sub.fwd_map[y][1] @test has_lower_bound(xm) @test lower_bound(xm) == 0 @test has_upper_bound(xm) From 560bdf6935b5c8d749699bad9c75b99c5eba9943 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Wed, 1 Apr 2026 12:39:48 -0400 Subject: [PATCH 13/13] singular prepare_max_M_objective --- src/mbm.jl | 130 +++++++++++++++++----------------------- test/constraints/mbm.jl | 51 +++++++--------- 2 files changed, 77 insertions(+), 104 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 79839edf..91a2f873 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -239,75 +239,51 @@ end ################################################################################ """ - prepare_objectives(model, obj::ScalarConstraint, sub::GDPSubmodel) + prepare_max_M_objective(model, obj::ScalarConstraint, sub::GDPSubmodel) -Convert a constraint into objective expressions for M-value maximization. -Returns a vector of JuMP expressions to pass to `_raw_M`. The base method -produces a single-element vector by mapping variables through `sub.fwd_map`. +Convert a constraint into an objective expression for M-value +maximization. Returns a single JuMP expression to pass to `_raw_M`. """ -function prepare_objectives( +function prepare_max_M_objective( ::JuMP.AbstractModel, obj::JuMP.ScalarConstraint{T, S}, sub::GDPSubmodel ) where {T, S <: _MOI.LessThan} flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map) - expr = -obj.set.upper +_replace_variables_in_constraint(obj.func, flat_map) - return [expr] + expr = -obj.set.upper + _replace_variables_in_constraint(obj.func, flat_map) + return expr end -function prepare_objectives( +function prepare_max_M_objective( ::JuMP.AbstractModel, obj::JuMP.ScalarConstraint{T, S}, sub::GDPSubmodel ) where {T, S <: _MOI.GreaterThan} flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map) - expr = obj.set.lower -_replace_variables_in_constraint(obj.func, flat_map) - return [expr] + expr = obj.set.lower - _replace_variables_in_constraint(obj.func, flat_map) + return expr end -# Solve the submodel for each objective and return raw M values as a vector, or -# nothing if infeasible. Order of results matches order of objectives. +# Solve the submodel for a single objective expression. +# Returns a scalar M value, or nothing if infeasible. function _raw_M( sub::GDPSubmodel, - objectives::Vector{<:JuMP.AbstractJuMPScalar}, + objective::JuMP.AbstractJuMPScalar, method::_MBM ) - M_vals = typeof(method.default_M)[] - for obj_expr in objectives - JuMP.@objective(sub.model, Max, obj_expr) - JuMP.optimize!(sub.model) - if JuMP.is_solved_and_feasible(sub.model) - push!(M_vals, max( - JuMP.objective_value(sub.model), - zero(method.default_M)) - ) - elseif JuMP.termination_status(sub.model) == _MOI.INFEASIBLE - return nothing - else - push!(M_vals, method.default_M) - end + JuMP.@objective(sub.model, Max, objective) + JuMP.optimize!(sub.model) + if JuMP.is_solved_and_feasible(sub.model) + return max(JuMP.objective_value(sub.model), zero(method.default_M)) + elseif JuMP.termination_status(sub.model) == _MOI.INFEASIBLE + return nothing + else + return method.default_M end - return M_vals -end - -""" - aggregate_M_values(model, vals::AbstractVector) - -Reduce a vector of raw M values from `_raw_M` to the final form used -during constraint reformulation. The base method returns the single -element from a length-1 vector. Extensions (e.g., InfiniteOpt) override -to aggregate across multiple support points (e.g., interpolating K -per-support M values into a parameter function). -""" -function aggregate_M_values( - ::JuMP.AbstractModel, - vals::AbstractVector - ) - return only(vals) end # Dispatch over constraint types to compute M values. Scalar -# LE/GE: prepare objectives, solve, finalize. +# LE/GE: prepare objective, solve, return scalar. function _maximize_M( model::JuMP.AbstractModel, objective::JuMP.ScalarConstraint{T, S}, @@ -315,10 +291,8 @@ function _maximize_M( method::_MBM ) where {T, S <: Union{_MOI.LessThan, _MOI.GreaterThan}} sub = _get_submodel(model, constraints, method) - objectives = prepare_objectives(model, objective, sub) - raw = _raw_M(sub, objectives, method) - raw === nothing && return nothing - return aggregate_M_values(model, raw) + return _raw_M(sub, + prepare_max_M_objective(model, objective, sub), method) end # Helper: get or create the submodel for a set of constraints. @@ -336,7 +310,7 @@ function _get_submodel( return method.model_cache[indicator] end -# EqualTo: solve both GreaterThan and LessThan directions, finalize each. +# EqualTo: solve both GreaterThan and LessThan directions. function _maximize_M( model::JuMP.AbstractModel, objective::JuMP.ScalarConstraint{T, S}, @@ -347,16 +321,14 @@ function _maximize_M( set_value = objective.set.value ge_obj = JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)) le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)) - raw_lower = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) - raw_upper = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) + raw_lower = _raw_M(sub, prepare_max_M_objective(model, ge_obj, sub), method) + raw_upper = _raw_M(sub, prepare_max_M_objective(model, le_obj, sub), method) (raw_lower === nothing || raw_upper === nothing) && return nothing - return [aggregate_M_values(model, raw_lower), - aggregate_M_values(model, raw_upper) - ] + return [raw_lower, raw_upper] end -# Interval: solve both lower and upper bound directions, finalize each. +# Interval: solve both lower and upper bound directions. function _maximize_M( model::JuMP.AbstractModel, objective::JuMP.ScalarConstraint{T, S}, @@ -369,17 +341,14 @@ function _maximize_M( MOI.GreaterThan(set_values[1])) le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_values[2])) - raw_lower = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) - raw_upper = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) + raw_lower = _raw_M(sub, prepare_max_M_objective(model, ge_obj, sub), method) + raw_upper = _raw_M(sub, prepare_max_M_objective(model, le_obj, sub), method) (raw_lower === nothing || raw_upper === nothing) && return nothing - return [ - aggregate_M_values(model, raw_lower), - aggregate_M_values(model, raw_upper) - ] + return [raw_lower, raw_upper] end -# Nonpositives: per-row LessThan solves for each dimension of the vector. +# Nonpositives: per-row LessThan solves. function _maximize_M( model::JuMP.AbstractModel, objective::JuMP.VectorConstraint{T, S, R}, @@ -392,14 +361,16 @@ function _maximize_M( for i in 1:objective.set.dimension le_obj = JuMP.ScalarConstraint( objective.func[i], MOI.LessThan(zero(val_type))) - raw = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) + raw = _raw_M(sub, + prepare_max_M_objective(model, le_obj, sub), + method) raw === nothing && return nothing - push!(results, aggregate_M_values(model, raw)) + push!(results, raw) end return results end -# Nonnegatives: per-row GreaterThan solves for each dimension of the vector. +# Nonnegatives: per-row GreaterThan solves. function _maximize_M( model::JuMP.AbstractModel, objective::JuMP.VectorConstraint{T, S, R}, @@ -411,15 +382,18 @@ function _maximize_M( results = Any[] for i in 1:objective.set.dimension ge_obj = JuMP.ScalarConstraint( - objective.func[i], MOI.GreaterThan(zero(val_type))) - raw = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) + objective.func[i], + MOI.GreaterThan(zero(val_type))) + raw = _raw_M(sub, + prepare_max_M_objective(model, ge_obj, sub), + method) raw === nothing && return nothing - push!(results, aggregate_M_values(model, raw)) + push!(results, raw) end return results end -# Zeros: per-row element-wise max of GE and LE raw values, then finalize. +# Zeros: per-row element-wise max of GE and LE values. function _maximize_M( model::JuMP.AbstractModel, objective::JuMP.VectorConstraint{T, S, R}, @@ -431,14 +405,20 @@ function _maximize_M( results = Any[] for i in 1:objective.set.dimension ge_obj = JuMP.ScalarConstraint( - objective.func[i], MOI.GreaterThan(zero(val_type))) + objective.func[i], + MOI.GreaterThan(zero(val_type))) le_obj = JuMP.ScalarConstraint( - objective.func[i], MOI.LessThan(zero(val_type))) - raw_ge = _raw_M(sub, prepare_objectives(model, ge_obj, sub), method) - raw_le = _raw_M(sub, prepare_objectives(model, le_obj, sub), method) + objective.func[i], + MOI.LessThan(zero(val_type))) + raw_ge = _raw_M(sub, + prepare_max_M_objective(model, ge_obj, sub), + method) + raw_le = _raw_M(sub, + prepare_max_M_objective(model, le_obj, sub), + method) (raw_ge === nothing || raw_le === nothing) && return nothing - push!(results, aggregate_M_values(model, max.(raw_ge, raw_le))) + push!(results, max(raw_ge, raw_le)) end return results end diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index ac04f217..d181c082 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -52,7 +52,7 @@ function test__replace_variables_quad_numeric_map() @test result3.aff.terms[y] ≈ 3.0 end -function test__replace_variables_in_constraint() +function test_replace_variables_in_constraint() model = Model() sub_model = Model() @variable(model, x[1:3]) @@ -84,7 +84,7 @@ function test__replace_variables_in_constraint() "String", new_vars) end -function test_prepare_objectives() +function test_prepare_max_M_objective() model = Model() sub_model = Model() @@ -99,17 +99,15 @@ function test_prepare_objectives() collect(keys(new_vars)), new_vars) # LessThan: max(f - upper) = max(x[1] - 1) - objs_le = DP.prepare_objectives( + obj_le = DP.prepare_max_M_objective( model, constraint_object(lessthan), sub) - @test length(objs_le) == 1 - @test objs_le[1] == JuMP.@expression(sub_model, + @test obj_le == JuMP.@expression(sub_model, new_vars[x[1]][1] - 1) # GreaterThan: max(lower - f) = max(1 - x[2]) - objs_ge = DP.prepare_objectives( + obj_ge = DP.prepare_max_M_objective( model, constraint_object(greaterthan), sub) - @test length(objs_ge) == 1 - @test objs_ge[1] == JuMP.@expression(sub_model, + @test obj_ge == JuMP.@expression(sub_model, 1 - new_vars[x[2]][1]) end @@ -134,43 +132,39 @@ function test_raw_M() DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con2], mbm) - objs = DP.prepare_objectives(model, + obj = DP.prepare_max_M_objective(model, constraint_object(con), sub) - raw = DP._raw_M(sub, objs, mbm) - @test DP.aggregate_M_values(model, raw) == 0.0 + @test DP._raw_M(sub, obj, mbm) == 0.0 set_upper_bound(x, 1) sub2 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con], mbm) - objs2 = DP.prepare_objectives(model, + obj2 = DP.prepare_max_M_objective(model, constraint_object(con2), sub2) - raw = DP._raw_M(sub2, objs2, mbm) - @test DP.aggregate_M_values(model, raw) == 15 + @test DP._raw_M(sub2, obj2, mbm) == 15 set_integer(y) @constraint(model, con3, y*x == 15, Disjunct(Y[1])) - objs3 = DP.prepare_objectives(model, + obj3 = DP.prepare_max_M_objective(model, constraint_object(con2), sub2) - raw = DP._raw_M(sub2, objs3, mbm) - @test DP.aggregate_M_values(model, raw) == 15 + @test DP._raw_M(sub2, obj3, mbm) == 15 # Fresh _MBM after changing bounds JuMP.fix(y, 5; force=true) mbm2 = DP._MBM( DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub3 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con], mbm2) - objs4 = DP.prepare_objectives(model, + obj4 = DP.prepare_max_M_objective(model, constraint_object(con2), sub3) - raw = DP._raw_M(sub3, objs4, mbm2) - @test DP.aggregate_M_values(model, raw) == 10 + @test DP._raw_M(sub3, obj4, mbm2) == 10 # Infeasible region → nothing delete_lower_bound(x) mbm3 = DP._MBM( DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub4 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[con2], mbm3) - objs5 = DP.prepare_objectives(model, + obj5 = DP.prepare_max_M_objective(model, constraint_object(con2), sub4) - @test DP._raw_M(sub4, objs5, mbm3) == nothing + @test DP._raw_M(sub4, obj5, mbm3) == nothing # infeasible (x >= 100 but x <= 1) set_upper_bound(x, 1) @@ -179,9 +173,9 @@ function test_raw_M() sub5 = DP.copy_model_with_constraints(model, DisjunctConstraintRef[truly_infeasible], mbm4) - objs6 = DP.prepare_objectives(model, + obj6 = DP.prepare_max_M_objective(model, constraint_object(con), sub5) - @test DP._raw_M(sub5, objs6, mbm4) == nothing + @test DP._raw_M(sub5, obj6, mbm4) == nothing # Unbounded subproblem → default_M fallback. # No lower bound on x means max(5 - x) s.t. x <= 3 @@ -195,10 +189,9 @@ function test_raw_M() mbm_ub = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) sub_ub = DP.copy_model_with_constraints(model_ub, DisjunctConstraintRef[ub_con1], mbm_ub) - objs_ub = DP.prepare_objectives(model_ub, + obj_ub = DP.prepare_max_M_objective(model_ub, constraint_object(ub_con2), sub_ub) - raw_ub = DP._raw_M(sub_ub, objs_ub, mbm_ub) - @test raw_ub == [mbm_ub.default_M] + @test DP._raw_M(sub_ub, obj_ub, mbm_ub) == mbm_ub.default_M end function test_maximize_M() @@ -801,8 +794,8 @@ end test_mbm() test__var_ref_type_numeric_map() test__replace_variables_quad_numeric_map() - test__replace_variables_in_constraint() - test_prepare_objectives() + test_replace_variables_in_constraint() + test_prepare_max_M_objective() test_raw_M() test_maximize_M() test_reformulate_disjunct_constraint()