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 6ba13159..866dbd2b 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -388,15 +388,20 @@ 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 + subproblem_indicators::Vector{LogicalVariableRef{M}} + # Cached submodels: indicator => GDPSubmodel. + # Typed Any so extensions can store different types. + model_cache::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}}(), + Dict{LogicalVariableRef{M}, Any}() ) end end @@ -432,35 +437,63 @@ mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationM end """ - cutting_planes{O} <: AbstractReformulationMethod + CuttingPlanes{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} <: AbstractReformulationMethod +struct CuttingPlanes{O, T} <: AbstractReformulationMethod optimizer::O; max_iter::Int - seperation_tolerance::Float64 + seperation_tolerance::T final_reform_method::AbstractReformulationMethod - M_value::Float64 - 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) + M_value::T + function CuttingPlanes( + 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 +################################################################################ + +""" + 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_map::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} + model::M + dec_vars::Vector{V} + fwd_map::Dict{V, Vector{W}} +end + """ PSplit <: AbstractReformulationMethod @@ -635,9 +668,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..91a2f873 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -1,140 +1,165 @@ +################################################################################ +# 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 function reformulate_disjunction( - model::JuMP.AbstractModel, + model::JuMP.AbstractModel, disj::Disjunction, method::MBM -) + ) mbm = _MBM(method, model) - ref_cons = Vector{JuMP.AbstractConstraint}() + ref_cons = Vector{JuMP.AbstractConstraint}() for d in disj.indicators - mbm.conlvref = filter(x -> x != d, disj.indicators) + mbm.subproblem_indicators = filter( + x -> x != d, disj.indicators) _reformulate_disjunct(model, ref_cons, d, mbm) 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, + model::JuMP.AbstractModel, + ref_cons::Vector{JuMP.AbstractConstraint}, + lvref::LogicalVariableRef, method::_MBM -) - - empty!(method.M) + ) !haskey(_indicator_to_constraints(model), lvref) && return - bconref = Dict(d => binary_variable(d) for d in method.conlvref) - + bconref = Dict( + d => binary_variable(d) for d in method.subproblem_indicators) + 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.subproblem_indicators + 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) + if M_result === nothing + 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, enforce globally (no reformulation needed). + 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 +167,376 @@ 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_max_M_objective(model, obj::ScalarConstraint, sub::GDPSubmodel) + +Convert a constraint into an objective expression for M-value +maximization. Returns a single JuMP expression to pass to `_raw_M`. +""" +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 +end + +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 +end + +# Solve the submodel for a single objective expression. +# Returns a scalar M value, or nothing if infeasible. +function _raw_M( + sub::GDPSubmodel, + objective::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 ) + 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 end +# Dispatch over constraint types to compute M values. Scalar +# LE/GE: prepare objective, solve, return scalar. 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) + return _raw_M(sub, + prepare_max_M_objective(model, objective, sub), method) +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.model_cache, indicator) + method.model_cache[indicator] = copy_model_with_constraints( + model, constraints, method) + end + return method.model_cache[indicator] end +# EqualTo: solve both GreaterThan and LessThan directions. 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_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 [raw_lower, raw_upper] end +# Interval: solve both lower and upper bound directions. 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_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 [raw_lower, raw_upper] end +# Nonpositives: per-row LessThan solves. 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_max_M_objective(model, le_obj, sub), + method) + raw === nothing && return nothing + push!(results, raw) + end + return results end +# Nonnegatives: per-row GreaterThan solves. 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_max_M_objective(model, ge_obj, sub), + method) + raw === nothing && return nothing + push!(results, raw) + end + return results end +# Zeros: per-row element-wise max of GE and LE values. 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_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, max(raw_ge, raw_le)) + end + return results 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}, +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 + +""" + copy_model_with_constraints(model, constraints, method) + +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 copy_model_with_constraints( + 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_map = Dict{var_type, Vector{var_type}}() + + for var in dec_vars + copy_var = variable_copy(sub_model, var) + fwd_map[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_map) + 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_map) 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 + 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) + JuMP.add_to_expression!(new_aff, coef, var_map[var]) end - new_aff.constant = fun.constant + new_aff.constant = new_aff.constant + fun.constant return new_aff end 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]) + 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) 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/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 ed74a5a5..d181c082 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -1,12 +1,58 @@ 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 -function test__replace_variables_in_constraint() +# _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() model = Model() sub_model = Model() @variable(model, x[1:3]) @@ -34,107 +80,234 @@ 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_max_M_objective() 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) + obj_le = DP.prepare_max_M_objective( + model, constraint_object(lessthan), sub) + @test obj_le == JuMP.@expression(sub_model, + new_vars[x[1]][1] - 1) + + # GreaterThan: max(lower - f) = max(1 - x[2]) + obj_ge = DP.prepare_max_M_objective( + model, constraint_object(greaterthan), sub) + @test obj_ge == 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.copy_model_with_constraints(model, + DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective(model, + constraint_object(con), sub) + @test DP._raw_M(sub, obj, mbm) == 0.0 set_upper_bound(x, 1) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con], mbm)== 15 + sub2 = DP.copy_model_with_constraints(model, + DisjunctConstraintRef[con], mbm) + obj2 = DP.prepare_max_M_objective(model, + constraint_object(con2), sub2) + @test DP._raw_M(sub2, obj2, mbm) == 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])) + obj3 = DP.prepare_max_M_objective(model, + constraint_object(con2), sub2) + @test DP._raw_M(sub2, obj3, mbm) == 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.copy_model_with_constraints(model, + DisjunctConstraintRef[con], mbm2) + obj4 = DP.prepare_max_M_objective(model, + constraint_object(con2), sub3) + @test DP._raw_M(sub3, obj4, mbm2) == 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.copy_model_with_constraints(model, + DisjunctConstraintRef[con2], mbm3) + obj5 = DP.prepare_max_M_objective(model, + constraint_object(con2), sub4) + @test DP._raw_M(sub4, obj5, mbm3) == nothing + + # infeasible (x >= 100 but x <= 1) + set_upper_bound(x, 1) + mbm4 = DP._MBM( + DP.MBM(HiGHS.Optimizer), JuMP.Model()) + sub5 = DP.copy_model_with_constraints(model, + DisjunctConstraintRef[truly_infeasible], + mbm4) + obj6 = DP.prepare_max_M_objective(model, + constraint_object(con), sub5) + @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 + # 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.copy_model_with_constraints(model_ub, + DisjunctConstraintRef[ub_con1], mbm_ub) + obj_ub = DP.prepare_max_M_objective(model_ub, + constraint_object(ub_con2), sub_ub) + @test DP._raw_M(sub_ub, obj_ub, mbm_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 +317,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 +434,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 +452,351 @@ 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 throws error + 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) + @test_throws ErrorException reformulate_disjunction( + model2, constraint_object(disj2), method2) + + 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 _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) + @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.copy_model_with_constraints(model, constraints, mbm) + + # Check submodel struct + @test sub isa DP.GDPSubmodel + @test sub.model isa Model + # Check variable mapping exists + @test haskey(sub.fwd_map, x) + @test haskey(sub.fwd_map, y) + + # Check mapped variables (length-1 vectors) + @test length(sub.fwd_map[x]) == 1 + @test length(sub.fwd_map[y]) == 1 + + # Check bounds in submodel + 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) + @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__copy_model_with_constraints() + test_get_variable_info() test_mbm() - test__replace_variables_in_constraint() - test__constraint_to_objective() - test_mini_model() + test__var_ref_type_numeric_map() + test__replace_variables_quad_numeric_map() + test_replace_variables_in_constraint() + test_prepare_max_M_objective() + test_raw_M() test_maximize_M() test_reformulate_disjunct_constraint() test_reformulate_disjunct() 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