From 9be42ec5d2729f604c424d167eaea0be07991e15 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 8 Jan 2026 14:37:32 -0500 Subject: [PATCH 01/13] Refactor disjunct reformulation to use per-constraint M values and improve test coverage --- src/mbm.jl | 47 +++++++++++++++++++------------------ test/constraints/mbm.jl | 51 ++++++++++++++++++++++------------------- 2 files changed, 53 insertions(+), 45 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index e9b1f13..c836053 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -15,39 +15,42 @@ function reformulate_disjunction( 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. +# Per Trespalacios & Grossmann (2015) Eq. (9), each constraint e in term i +# gets its own set of M_{ie,i'} values for each other term i'. function _reformulate_disjunct( - model::JuMP.AbstractModel, - ref_cons::Vector{JuMP.AbstractConstraint}, + 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) - + 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), + # For each constraint, compute its own set of M values + for cref in filtered_constraints + empty!(method.M) # Clear M for each constraint + + 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] = _maximize_M( + model, + JuMP.constraint_object(cref), disjunct_constraints, method - ) for cref in filtered_constraints - ) + ) + end end - end - for cref in filtered_constraints - con = JuMP.constraint_object(cref) - append!(ref_cons, reformulate_disjunct_constraint(model, con, + + con = JuMP.constraint_object(cref) + append!(ref_cons, reformulate_disjunct_constraint(model, con, bconref, method)) end return ref_cons diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index ed74a5a..2b1c901 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -177,22 +177,21 @@ function test_reformulate_disjunct_constraint() @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)))) && + @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)))) && + @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, + + @test length(reformulated_constraints[7]) >= 2 + @test reformulated_constraints[7][1].set == MOI.LessThan(1.0) + @test reformulated_constraints[7][2].set == MOI.GreaterThan(1.0) + # Verify x[1] has coefficient 1.0 in both constraints + @test JuMP.coefficient(reformulated_constraints[7][1].func, x[1]) == 1.0 + @test JuMP.coefficient(reformulated_constraints[7][2].func, x[1]) == 1.0 + + @test_throws ErrorException reformulate_disjunct_constraint(model, "odd", bconref, method) end @@ -238,30 +237,36 @@ 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 @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) - 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 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]: _maximize_M for Interval takes max(M_lower, M_upper) + # - M_lower (x >= 0): max(-x) s.t. 1<=x<=2 → -1 + # - M_upper (x <= 55): max(x-55) s.t. 1<=x<=2 → -53 + func_1 = ref_cons[1].func # x - 53*Y[2] <= 2.0 + func_2 = ref_cons[2].func # x + 1*Y[2] >= 1.0 (per-constraint M=1) + func_3 = ref_cons[3].func # x + (-1)*Y[1] >= 0.0 + func_4 = ref_cons[4].func # x - (-1)*Y[1] <= 55.0 → x + Y[1] <= 55.0 @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 @test JuMP.coefficient(func_3, x) == 1.0 @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -1.0 From 0d66676fc249796286f9e37d88f4f878cef5ab83 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 8 Jan 2026 14:54:00 -0500 Subject: [PATCH 02/13] removing extra autocomplete coments --- src/mbm.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index c836053..1eb0e96 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -15,8 +15,8 @@ function reformulate_disjunction( end return ref_cons end + # Reformulates a disjunct represented by lvref using per-constraint M values. -# Per Trespalacios & Grossmann (2015) Eq. (9), each constraint e in term i # gets its own set of M_{ie,i'} values for each other term i'. function _reformulate_disjunct( model::JuMP.AbstractModel, @@ -32,7 +32,7 @@ function _reformulate_disjunct( # For each constraint, compute its own set of M values for cref in filtered_constraints - empty!(method.M) # Clear M for each constraint + empty!(method.M) for d in method.conlvref d_constraints = _indicator_to_constraints(model)[d] From d21c4e286822f5cd1a29d9e044ca25feb751115f Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 8 Jan 2026 17:30:40 -0500 Subject: [PATCH 03/13] Enhance MBM structure to support per-row and per-bound M values in constraints --- src/datatypes.jl | 10 +- src/mbm.jl | 216 +++++++++++++++++++------------------ test/constraints/mbm.jl | 230 ++++++++++++++++++++++++++-------------- 3 files changed, 272 insertions(+), 184 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index 6ba1315..387e49b 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -388,15 +388,15 @@ 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}, Union{T, Vector{T}}} + default_M::T + conlvref::Vector{LogicalVariableRef{M}} function _MBM(method::MBM{O, T}, model::M) where {O, T, M <: JuMP.AbstractModel} new{O, T, M}(method.optimizer, - Dict{LogicalVariableRef{M}, T}(), + Dict{LogicalVariableRef{M}, Union{T, Vector{T}}}(), method.default_M, - Vector{LogicalVariableRef{M}}() + Vector{LogicalVariableRef{M}}() ) end end diff --git a/src/mbm.jl b/src/mbm.jl index 1eb0e96..af1adac 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -74,6 +74,7 @@ function reformulate_disjunct_constraint( return new_ref_cons end +# Per-row M values: method.M[d] is a Vector, index with [j] for row j function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.VectorConstraint{T, S, R}, @@ -81,15 +82,13 @@ function reformulate_disjunct_constraint( 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)) - new_func = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] - m_sum + new_func = JuMP.@expression(model, [j=1:con.set.dimension], + con.func[j] - sum(method.M[d][j] * bconref[d] for d in keys(method.M)) ) reform_con = JuMP.build_constraint(error, new_func, con.set) return [reform_con] end - function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.VectorConstraint{T, S, R}, @@ -97,9 +96,8 @@ function reformulate_disjunct_constraint( 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)) - new_func = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] + m_sum + new_func = JuMP.@expression(model, [j=1:con.set.dimension], + con.func[j] + sum(method.M[d][j] * bconref[d] for d in keys(method.M)) ) reform_con = JuMP.build_constraint(error, new_func, con.set) return [reform_con] @@ -112,17 +110,16 @@ function reformulate_disjunct_constraint( 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)) - upper_expr = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] + m_sum + upper_expr = JuMP.@expression(model, [j=1:con.set.dimension], + con.func[j] + sum(method.M[d][j] * bconref[d] for d in keys(method.M)) ) - lower_expr = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] - m_sum + lower_expr = JuMP.@expression(model, [j=1:con.set.dimension], + con.func[j] - sum(method.M[d][j] * bconref[d] for d in keys(method.M)) ) - upper_con = JuMP.build_constraint(error, upper_expr, + upper_con = JuMP.build_constraint(error, upper_expr, MOI.Nonnegatives(con.set.dimension) ) - lower_con = JuMP.build_constraint(error, lower_expr, + lower_con = JuMP.build_constraint(error, lower_expr, MOI.Nonpositives(con.set.dimension) ) return [upper_con, lower_con] @@ -156,6 +153,7 @@ function reformulate_disjunct_constraint( return [reform_con] end +# Per-bound M values: method.M[d] = [M_lower, M_upper] function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.ScalarConstraint{T, S}, @@ -163,18 +161,18 @@ function reformulate_disjunct_constraint( 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)) + lower_func = JuMP.@expression(model, + con.func + sum(method.M[d][1] * bconref[d] for d in keys(method.M)) ) - upper_con = JuMP.build_constraint(error, upper_func, - MOI.LessThan(con.set.value) + 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, + 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 @@ -185,21 +183,19 @@ function reformulate_disjunct_constraint( 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)) + set_values = _set_values(con.set) + lower_func = JuMP.@expression(model, + con.func + sum(method.M[d][1] * bconref[d] for d 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)) + 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, + 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 @@ -219,65 +215,74 @@ end ################################################################################ # Dispatches over constraint types to reformulate into >= or <= # in order to solve the mini-model +# Per-row M values for vector constraints: returns Vector{T} function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM ) where { T, S <: _MOI.Nonpositives, R} val_type = JuMP.value_type(typeof(model)) - return maximum( + return [ _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), - constraints, + model, + JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), + constraints, method ) for i in 1:objective.set.dimension - ) + ] end function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM ) where { T, S <: _MOI.Nonnegatives, R} val_type = JuMP.value_type(typeof(model)) - return maximum( + return [ _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(zero(val_type))), - constraints, + model, + JuMP.ScalarConstraint( + objective.func[i], + MOI.GreaterThan(zero(val_type)) + ), + constraints, method ) for i in 1:objective.set.dimension - ) + ] end +# For Zeros, each row is an equality: return per-row M = max(M_lower, M_upper) function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM ) where { T, S <: _MOI.Zeros, R} val_type = JuMP.value_type(typeof(model)) - return max( - maximum( + return [ + max( _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i],MOI.GreaterThan(zero(val_type))), - constraints, + 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, + model, + JuMP.ScalarConstraint( + objective.func[i], + MOI.LessThan(zero(val_type)) + ), + constraints, method - ) for i in 1:objective.set.dimension - ) - ) + ) + ) for i in 1:objective.set.dimension + ] end function _maximize_M( @@ -289,56 +294,65 @@ function _maximize_M( return _mini_model(model, objective, constraints, method) end +# Per-bound M values for bidirectional constraints: returns [M_lower, M_upper] 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 <: _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 - ) + M_lower = _mini_model( + model, + JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)), + constraints, + method + ) + M_upper = _mini_model( + model, + JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)), + constraints, + method ) + return [M_lower, M_upper] end 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 <: _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 - ) + M_lower = _mini_model( + model, + JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_values[1])), + constraints, + method ) + M_upper = _mini_model( + model, + JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_values[2])), + constraints, + method + ) + return [M_lower, M_upper] end +# Nested Disjunctions: M computation skipped, handler creates fresh MBM function _maximize_M( - ::JuMP.AbstractModel, - ::F, - ::Vector{<:DisjunctConstraintRef}, + ::JuMP.AbstractModel, + ::Disjunction, + ::Vector{<:DisjunctConstraintRef}, + ::_MBM +) + return nothing +end + +function _maximize_M( + ::JuMP.AbstractModel, + ::F, + ::Vector{<:DisjunctConstraintRef}, ::_MBM ) where {F} error("This type of constraints and objective constraint has " * diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index 2b1c901..2298d8a 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -91,49 +91,78 @@ 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 + @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, -45.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] + + # zeros: x == 1 against Y[2] - per-row M values differ! + # 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) end @@ -144,56 +173,100 @@ 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 length(reformulated_constraints[7]) >= 2 - @test reformulated_constraints[7][1].set == MOI.LessThan(1.0) - @test reformulated_constraints[7][2].set == MOI.GreaterThan(1.0) - # Verify x[1] has coefficient 1.0 in both constraints - @test JuMP.coefficient(reformulated_constraints[7][1].func, x[1]) == 1.0 - @test JuMP.coefficient(reformulated_constraints[7][2].func, x[1]) == 1.0 + 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) # Outer indicators + @variable(model2, Inner[1:2], Logical) # Inner disjunction indicators + # Inner disjunction constraints (will be nested inside Outer[1]) + @constraint(model2, inner_lt, z <= 1, Disjunct(Inner[1])) + @constraint(model2, inner_gt, z >= 1, Disjunct(Inner[2])) + @disjunction(model2, inner_disj, Inner) + # bconref contains ONLY the OTHER outer indicator (Outer[2]) + # This is what happens when reformulating Outer[1]'s constraints + 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 # Scalar M for outer reformulation + ref_disjunction = reformulate_disjunct_constraint( + model2, constraint_object(inner_disj), bconref2, method_nested) + @test length(ref_disjunction) >= 2 + # Verify structure: inner reformulation + outer Big-M terms + @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) - + "odd", bconref, method_scalar) end function test_reformulate_disjunct() @@ -222,11 +295,12 @@ function test_reformulate_disjunct() @test JuMP.coefficient(func_1, x[1]) == 1.0 @test JuMP.coefficient(func_1, binary_variable(Y[2])) == -1.5 + # 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() @@ -251,28 +325,28 @@ function test_reformulate_disjunction() @test ref_cons[4].set == MOI.LessThan(55.0) - # Per-constraint M values: + # 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]: _maximize_M for Interval takes max(M_lower, M_upper) - # - M_lower (x >= 0): max(-x) s.t. 1<=x<=2 → -1 - # - M_upper (x <= 55): max(x-55) s.t. 1<=x<=2 → -53 + # - interval in Y[1] region (1 <= x <= 2): + # - M_lower (x >= 0): max(0-x) at x=1 → M_lower=-1 + # - M_upper (x <= 55): max(x-55) at x=2 → M_upper=-53 func_1 = ref_cons[1].func # x - 53*Y[2] <= 2.0 - func_2 = ref_cons[2].func # x + 1*Y[2] >= 1.0 (per-constraint M=1) - func_3 = ref_cons[3].func # x + (-1)*Y[1] >= 0.0 - func_4 = ref_cons[4].func # x - (-1)*Y[1] <= 55.0 → x + Y[1] <= 55.0 + func_2 = ref_cons[2].func # x + 1*Y[2] >= 1.0 + func_3 = ref_cons[3].func # x + M_lower*Y[1] >= 0.0 → x + (-1)*Y[1] >= 0 + func_4 = ref_cons[4].func # x - M_upper*Y[1] <= 55 → x - (-53)*Y[1] <= 55 @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])) == 1.0 + @test JuMP.coefficient(func_2, binary_variable(Y[2])) == 1.0 @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 JuMP.coefficient(func_4, binary_variable(Y[1])) == 53.0 # -M_upper = -(-53) = 53 end @testset "MBM" begin From 13156b62d6aa139e98c4d968f68db89c1ef8491e Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 8 Jan 2026 17:32:24 -0500 Subject: [PATCH 04/13] Refactor reformulate_disjunct_constraint functions to clarify per-row and per-bound M value usage --- src/mbm.jl | 160 +++++++++++++++++++++++++---------------------------- 1 file changed, 75 insertions(+), 85 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index af1adac..624f064 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -74,7 +74,7 @@ function reformulate_disjunct_constraint( return new_ref_cons end -# Per-row M values: method.M[d] is a Vector, index with [j] for row j +# 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}, @@ -82,13 +82,14 @@ function reformulate_disjunct_constraint( Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, method::_MBM ) where {T, S <: _MOI.Nonpositives, R} - new_func = JuMP.@expression(model, [j=1:con.set.dimension], - con.func[j] - sum(method.M[d][j] * bconref[d] for d in keys(method.M)) + new_func = JuMP.@expression(model, [i=1:con.set.dimension], + 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}, @@ -96,30 +97,31 @@ function reformulate_disjunct_constraint( Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, method::_MBM ) where {T, S <: _MOI.Nonnegatives, R} - new_func = JuMP.@expression(model, [j=1:con.set.dimension], - con.func[j] + sum(method.M[d][j] * bconref[d] for d in keys(method.M)) + new_func = JuMP.@expression(model, [i=1:con.set.dimension], + 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}, + bconref::Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, method::_MBM ) where {T, S <: _MOI.Zeros, R} - upper_expr = JuMP.@expression(model, [j=1:con.set.dimension], - con.func[j] + sum(method.M[d][j] * bconref[d] for d in keys(method.M)) + upper_expr = JuMP.@expression(model, [i=1:con.set.dimension], + con.func[i] + sum(method.M[d][i] * bconref[d] for d in keys(method.M)) ) - lower_expr = JuMP.@expression(model, [j=1:con.set.dimension], - con.func[j] - sum(method.M[d][j] * bconref[d] for d in keys(method.M)) + lower_expr = JuMP.@expression(model, [i=1: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, + upper_con = JuMP.build_constraint(error, upper_expr, MOI.Nonnegatives(con.set.dimension) ) - lower_con = JuMP.build_constraint(error, lower_expr, + lower_con = JuMP.build_constraint(error, lower_expr, MOI.Nonpositives(con.set.dimension) ) return [upper_con, lower_con] @@ -153,7 +155,7 @@ function reformulate_disjunct_constraint( return [reform_con] end -# Per-bound M values: method.M[d] = [M_lower, M_upper] +# 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}, @@ -161,6 +163,7 @@ function reformulate_disjunct_constraint( Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, method::_MBM ) where {T, S <: _MOI.EqualTo} + # M[d][1] = M for GreaterThan (lower bound), 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)) ) @@ -176,6 +179,7 @@ function reformulate_disjunct_constraint( 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}, @@ -184,6 +188,7 @@ function reformulate_disjunct_constraint( method::_MBM ) where {T, S <: _MOI.Interval} set_values = _set_values(con.set) + # M[d][1] = M for GreaterThan (lower bound), 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)) ) @@ -215,13 +220,13 @@ end ################################################################################ # Dispatches over constraint types to reformulate into >= or <= # in order to solve the mini-model -# Per-row M values for vector constraints: returns Vector{T} +# Returns Vector{T} - one M per row for tighter per-row relaxations function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where { T, S <: _MOI.Nonpositives, R} +) where {T, S <: _MOI.Nonpositives, R} val_type = JuMP.value_type(typeof(model)) return [ _maximize_M( @@ -233,51 +238,43 @@ function _maximize_M( ] end +# Returns Vector{T} - one M per row for tighter per-row relaxations function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where { T, S <: _MOI.Nonnegatives, R} +) where {T, S <: _MOI.Nonnegatives, R} val_type = JuMP.value_type(typeof(model)) return [ _maximize_M( model, - JuMP.ScalarConstraint( - objective.func[i], - MOI.GreaterThan(zero(val_type)) - ), + JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(zero(val_type))), constraints, method ) for i in 1:objective.set.dimension ] end -# For Zeros, each row is an equality: return per-row M = max(M_lower, M_upper) +# Returns Vector{T} - one M per row, each is max(M_ge, M_le) for that row function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM -) where { T, S <: _MOI.Zeros, R} +) where {T, S <: _MOI.Zeros, R} val_type = JuMP.value_type(typeof(model)) return [ max( _maximize_M( model, - JuMP.ScalarConstraint( - objective.func[i], - MOI.GreaterThan(zero(val_type)) - ), + JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(zero(val_type))), constraints, method ), _maximize_M( model, - JuMP.ScalarConstraint( - objective.func[i], - MOI.LessThan(zero(val_type)) - ), + JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), constraints, method ) @@ -294,65 +291,58 @@ function _maximize_M( return _mini_model(model, objective, constraints, method) end -# Per-bound M values for bidirectional constraints: returns [M_lower, M_upper] +# Returns [M_lower, M_upper] for per-bound relaxations 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 <: _MOI.EqualTo} set_value = objective.set.value - M_lower = _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)), - constraints, - method - ) - M_upper = _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)), - constraints, - method - ) - return [M_lower, M_upper] + return [ + _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 + ) + ] end +# Returns [M_lower, M_upper] for per-bound relaxations 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 <: _MOI.Interval} set_values = _set_values(objective.set) # Returns (lower, upper) - M_lower = _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_values[1])), - constraints, - method - ) - M_upper = _mini_model( - model, - JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_values[2])), - constraints, - method - ) - return [M_lower, M_upper] -end - -# Nested Disjunctions: M computation skipped, handler creates fresh MBM -function _maximize_M( - ::JuMP.AbstractModel, - ::Disjunction, - ::Vector{<:DisjunctConstraintRef}, - ::_MBM -) - return nothing + return [ + _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 + ) + ] end function _maximize_M( - ::JuMP.AbstractModel, - ::F, - ::Vector{<:DisjunctConstraintRef}, + ::JuMP.AbstractModel, + ::F, + ::Vector{<:DisjunctConstraintRef}, ::_MBM ) where {F} error("This type of constraints and objective constraint has " * From 74fc0943da9fe09c8e95cf9729bdbbda6e68aeac Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Thu, 8 Jan 2026 22:31:27 -0500 Subject: [PATCH 05/13] Refactor code formatting for improved readability in MBM functions and tests --- src/mbm.jl | 34 +++++++++++++++++++++++++--------- test/constraints/mbm.jl | 10 +++++++--- 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 624f064..871b06c 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -28,7 +28,9 @@ function _reformulate_disjunct( bconref = Dict(d => binary_variable(d) for d in method.conlvref) constraints = _indicator_to_constraints(model)[lvref] - filtered_constraints = [c for c in constraints if c isa DisjunctConstraintRef] + 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 @@ -163,7 +165,8 @@ function reformulate_disjunct_constraint( Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, method::_MBM ) where {T, S <: _MOI.EqualTo} - # M[d][1] = M for GreaterThan (lower bound), M[d][2] = M for LessThan (upper) + # M[d][1] = M for GreaterThan (lower bound), M[d][2] = M for LessThan + # (upper bound) lower_func = JuMP.@expression(model, con.func + sum(method.M[d][1] * bconref[d] for d in keys(method.M)) ) @@ -188,7 +191,8 @@ function reformulate_disjunct_constraint( method::_MBM ) where {T, S <: _MOI.Interval} set_values = _set_values(con.set) - # M[d][1] = M for GreaterThan (lower bound), M[d][2] = M for LessThan (upper) + # M[d][1] = M for GreaterThan (lower bound), M[d][2] = M for LessThan + # (upper bound) lower_func = JuMP.@expression(model, con.func + sum(method.M[d][1] * bconref[d] for d in keys(method.M)) ) @@ -231,7 +235,9 @@ function _maximize_M( return [ _maximize_M( model, - JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), + JuMP.ScalarConstraint( + objective.func[i], MOI.LessThan(zero(val_type)) + ), constraints, method ) for i in 1:objective.set.dimension @@ -249,7 +255,9 @@ function _maximize_M( return [ _maximize_M( model, - JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(zero(val_type))), + JuMP.ScalarConstraint( + objective.func[i], MOI.GreaterThan(zero(val_type)) + ), constraints, method ) for i in 1:objective.set.dimension @@ -268,13 +276,17 @@ function _maximize_M( max( _maximize_M( model, - JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(zero(val_type))), + JuMP.ScalarConstraint( + objective.func[i], MOI.GreaterThan(zero(val_type)) + ), constraints, method ), _maximize_M( model, - JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), + JuMP.ScalarConstraint( + objective.func[i], MOI.LessThan(zero(val_type)) + ), constraints, method ) @@ -326,13 +338,17 @@ function _maximize_M( return [ _mini_model( model, - JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_values[1])), + JuMP.ScalarConstraint( + objective.func, MOI.GreaterThan(set_values[1]) + ), constraints, method ), _mini_model( model, - JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_values[2])), + JuMP.ScalarConstraint( + objective.func, MOI.LessThan(set_values[2]) + ), constraints, method ) diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index 2298d8a..74ce474 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -2,7 +2,9 @@ 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 end @@ -326,7 +328,8 @@ function test_reformulate_disjunction() @test ref_cons[4].set == MOI.LessThan(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 + # - 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 @@ -346,7 +349,8 @@ function test_reformulate_disjunction() @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])) == 53.0 # -M_upper = -(-53) = 53 + # -M_upper = -(-53) = 53 + @test JuMP.coefficient(func_4, binary_variable(Y[1])) == 53.0 end @testset "MBM" begin From d1c9fb608225e34d942f92abdb33407facd2dbd6 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Thu, 8 Jan 2026 22:37:24 -0500 Subject: [PATCH 06/13] comments for test --- test/constraints/mbm.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index 74ce474..e79573d 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -249,21 +249,18 @@ function test_reformulate_disjunct_constraint() # Create outer disjunct with inner disjunction model2 = GDPModel() @variable(model2, 0 <= z <= 50) - @variable(model2, Outer[1:2], Logical) # Outer indicators - @variable(model2, Inner[1:2], Logical) # Inner disjunction indicators - # Inner disjunction constraints (will be nested inside Outer[1]) + @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) - # bconref contains ONLY the OTHER outer indicator (Outer[2]) - # This is what happens when reformulating Outer[1]'s constraints 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 # Scalar M for outer reformulation + 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 - # Verify structure: inner reformulation + outer Big-M terms @test JuMP.coefficient(ref_disjunction[1].func, z) == 1.0 @test JuMP.coefficient(ref_disjunction[2].func, z) == 1.0 From 84b7e8287753f48dac5e02d8d88afb48c2bb7662 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Sun, 11 Jan 2026 18:51:22 -0500 Subject: [PATCH 07/13] max(M,0) --- src/mbm.jl | 2 +- test/constraints/mbm.jl | 24 ++++++++++++------------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 871b06c..b2cef00 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -394,7 +394,7 @@ function _mini_model( else M = JuMP.objective_value(sub_model) end - return M + return max(M, 0) end ################################################################################ diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index e79573d..5e4e492 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -74,8 +74,8 @@ function test_mini_model() @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 + @test DP._mini_model(model, constraint_object(con), + DisjunctConstraintRef[con2], mbm)== 0.0 set_upper_bound(x, 1) @test DP._mini_model(model, constraint_object(con2), DisjunctConstraintRef[con], mbm)== 15 @@ -113,11 +113,11 @@ function test_maximize_M() # 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 + # 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, -45.0] + mbm) == [0.0, 0.0] # Scalar LessThan/GreaterThan still return scalars # lessthan: x[1] <= 1 vs interval 0 <= x[1] <= 55 @@ -292,7 +292,7 @@ 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 @@ -329,12 +329,12 @@ function test_reformulate_disjunction() # 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 - # - M_upper (x <= 55): max(x-55) at x=2 → M_upper=-53 + # - 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 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 + M_lower*Y[1] >= 0.0 → x + (-1)*Y[1] >= 0 - func_4 = ref_cons[4].func # x - M_upper*Y[1] <= 55 → x - (-53)*Y[1] <= 55 + func_3 = ref_cons[3].func # x + M_lower*Y[1] >= 0.0 → x + 0*Y[1] >= 0 + func_4 = ref_cons[4].func # x - M_upper*Y[1] <= 55 → x - 0*Y[1] <= 55 @test JuMP.coefficient(func_1, x) == 1.0 @test JuMP.coefficient(func_1, binary_variable(Y[2])) == -53.0 @@ -343,11 +343,11 @@ function test_reformulate_disjunction() @test JuMP.coefficient(func_2, binary_variable(Y[2])) == 1.0 @test JuMP.coefficient(func_3, x) == 1.0 - @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -1.0 + @test JuMP.coefficient(func_3, binary_variable(Y[1])) == 0.0 @test JuMP.coefficient(func_4, x) == 1.0 - # -M_upper = -(-53) = 53 - @test JuMP.coefficient(func_4, binary_variable(Y[1])) == 53.0 + # -M_upper = -0 = 0 + @test JuMP.coefficient(func_4, binary_variable(Y[1])) == 0.0 end @testset "MBM" begin From c984a7a5bf730f89aa0f286601c1537063595722 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 13 Jan 2026 13:19:55 -0500 Subject: [PATCH 08/13] Infeasible and M = 0 enhancements --- src/datatypes.jl | 4 +- src/mbm.jl | 236 ++++++++++++++++++++++++++++++----------------- test/solve.jl | 12 ++- 3 files changed, 165 insertions(+), 87 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index 387e49b..5247d9a 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -391,12 +391,14 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho M::Dict{LogicalVariableRef{M}, Union{T, Vector{T}}} default_M::T conlvref::Vector{LogicalVariableRef{M}} + deactivated::Set{LogicalVariableRef{M}} function _MBM(method::MBM{O, T}, model::M) where {O, T, M <: JuMP.AbstractModel} new{O, T, M}(method.optimizer, Dict{LogicalVariableRef{M}, Union{T, Vector{T}}}(), method.default_M, - Vector{LogicalVariableRef{M}}() + Vector{LogicalVariableRef{M}}(), + Set{LogicalVariableRef{M}}() ) end end diff --git a/src/mbm.jl b/src/mbm.jl index b2cef00..4d2dff4 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -1,17 +1,53 @@ +################################################################################ +# HELPER FUNCTIONS +################################################################################ +# Check if M result (scalar or vector) contains only zeros +function _is_all_zeros(M) + if M isa Number + return M == 0 + elseif M isa AbstractVector + return all(m == 0 for m in M) + end + return false +end + ################################################################################ # CONSTRAINT, DISJUNCTION, DISJUNCT REFORMULATION ################################################################################ -#Reformulates the disjunction using multiple big-M values +# 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}() + # Track index ranges for each disjunct's constraints: (start, end) + disjunct_ranges = Dict{LogicalVariableRef, Tuple{Int, Int}}() for d in disj.indicators - mbm.conlvref = filter(x -> x != d, disj.indicators) + d in mbm.deactivated && continue + mbm.conlvref = filter( + x -> x != d && !(x in mbm.deactivated), + disj.indicators + ) + start_idx = length(ref_cons) + 1 _reformulate_disjunct(model, ref_cons, d, mbm) + end_idx = length(ref_cons) + if end_idx >= start_idx + disjunct_ranges[d] = (start_idx, end_idx) + end + end + # Remove constraints from deactivated disjuncts (in reverse order) + indices_to_remove = Int[] + for deact in mbm.deactivated + if haskey(disjunct_ranges, deact) + start_idx, end_idx = disjunct_ranges[deact] + append!(indices_to_remove, start_idx:end_idx) + end + end + sort!(indices_to_remove, rev=true) + for idx in indices_to_remove + deleteat!(ref_cons, idx) end return ref_cons end @@ -25,7 +61,9 @@ function _reformulate_disjunct( 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 + 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 = [ @@ -34,26 +72,46 @@ function _reformulate_disjunct( # For each constraint, compute its own set of M values for cref in filtered_constraints - empty!(method.M) + 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) - method.M[d] = _maximize_M( + 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." + # Remove from bconref since it's now deactivated + delete!(bconref, d) + else + method.M[d] = M_result + end end end con = JuMP.constraint_object(cref) - append!(ref_cons, reformulate_disjunct_constraint(model, con, - bconref, method)) + # Check if all M values are zero (constraint is global) + if !isempty(method.M) && + all(_is_all_zeros(method.M[d]) for d in keys(method.M)) + @info "Constraint is global (M ≤ 0 for all disjuncts), " * + "adding without Big-M relaxation." + push!(ref_cons, con) + else + append!(ref_cons, reformulate_disjunct_constraint(model, con, + bconref, method)) + end end return ref_cons end @@ -224,15 +282,15 @@ end ################################################################################ # Dispatches over constraint types to reformulate into >= or <= # in order to solve the mini-model -# Returns Vector{T} - one M per row for tighter per-row relaxations +# Returns Vector{T} or nothing - one M per row for tighter per-row relaxations function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM ) where {T, S <: _MOI.Nonpositives, R} val_type = JuMP.value_type(typeof(model)) - return [ + results = [ _maximize_M( model, JuMP.ScalarConstraint( @@ -242,17 +300,19 @@ function _maximize_M( method ) for i in 1:objective.set.dimension ] + any(r === nothing for r in results) && return nothing + return results end -# Returns Vector{T} - one M per row for tighter per-row relaxations +# Returns Vector{T} or nothing - one M per row for tighter per-row relaxations function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM ) where {T, S <: _MOI.Nonnegatives, R} val_type = JuMP.value_type(typeof(model)) - return [ + results = [ _maximize_M( model, JuMP.ScalarConstraint( @@ -262,36 +322,40 @@ function _maximize_M( method ) for i in 1:objective.set.dimension ] + any(r === nothing for r in results) && return nothing + return results end -# Returns Vector{T} - one M per row, each is max(M_ge, M_le) for that row +# Returns Vector{T} or nothing - one M per row, each is max(M_ge, M_le) function _maximize_M( - model::JuMP.AbstractModel, - objective::JuMP.VectorConstraint{T, S, R}, - constraints::Vector{<:DisjunctConstraintRef}, + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, method::_MBM ) where {T, S <: _MOI.Zeros, R} val_type = JuMP.value_type(typeof(model)) - return [ - max( - _maximize_M( - model, - JuMP.ScalarConstraint( - objective.func[i], MOI.GreaterThan(zero(val_type)) - ), - constraints, - method + results = [] + for i in 1:objective.set.dimension + M_ge = _maximize_M( + model, + JuMP.ScalarConstraint( + objective.func[i], MOI.GreaterThan(zero(val_type)) ), - _maximize_M( - model, - JuMP.ScalarConstraint( - objective.func[i], MOI.LessThan(zero(val_type)) - ), - constraints, - method - ) - ) for i in 1:objective.set.dimension - ] + constraints, + method + ) + M_le = _maximize_M( + model, + JuMP.ScalarConstraint( + objective.func[i], MOI.LessThan(zero(val_type)) + ), + constraints, + method + ) + (M_ge === nothing || M_le === nothing) && return nothing + push!(results, max(M_ge, M_le)) + end + return results end function _maximize_M( @@ -303,56 +367,56 @@ function _maximize_M( return _mini_model(model, objective, constraints, method) end -# Returns [M_lower, M_upper] for per-bound relaxations +# Returns [M_lower, M_upper] or nothing for per-bound relaxations 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 <: _MOI.EqualTo} set_value = objective.set.value - return [ - _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 - ) - ] + M_lower = _mini_model( + model, + JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value)), + constraints, + method + ) + M_upper = _mini_model( + model, + JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value)), + constraints, + method + ) + (M_lower === nothing || M_upper === nothing) && return nothing + return [M_lower, M_upper] end -# Returns [M_lower, M_upper] for per-bound relaxations +# Returns [M_lower, M_upper] or nothing for per-bound relaxations 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 <: _MOI.Interval} set_values = _set_values(objective.set) # Returns (lower, upper) - return [ - _mini_model( - model, - JuMP.ScalarConstraint( - objective.func, MOI.GreaterThan(set_values[1]) - ), - constraints, - method + M_lower = _mini_model( + model, + JuMP.ScalarConstraint( + objective.func, MOI.GreaterThan(set_values[1]) ), - _mini_model( - model, - JuMP.ScalarConstraint( - objective.func, MOI.LessThan(set_values[2]) - ), - constraints, - method - ) - ] + constraints, + method + ) + M_upper = _mini_model( + model, + JuMP.ScalarConstraint( + objective.func, MOI.LessThan(set_values[2]) + ), + constraints, + method + ) + (M_lower === nothing || M_upper === nothing) && return nothing + return [M_lower, M_upper] end function _maximize_M( @@ -387,9 +451,13 @@ function _mini_model( 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 + status = JuMP.termination_status(sub_model) + # Detect infeasibility: other disjunct has empty feasible region + if status == MOI.INFEASIBLE + return nothing + elseif status != 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) diff --git a/test/solve.jl b/test/solve.jl index b8b8db3..c78b032 100644 --- a/test/solve.jl +++ b/test/solve.jl @@ -7,7 +7,7 @@ function test_linear_gdp_example(m, use_complements = false) @variable(m, Y2, Logical, logical_complement = Y1) Y = [Y1, Y2] else - @variable(m, Y[1:2], Logical) + @variable(m, Y[1:3], Logical) end @variable(m, W[1:2], Logical) @objective(m, Max, sum(x)) @@ -16,7 +16,13 @@ function test_linear_gdp_example(m, use_complements = false) @constraint(m, w2[i=1:2], [2,4][i] ≤ x[i] ≤ [3,5][i], Disjunct(W[2])) @constraint(m, y2[i=1:2], [8,1][i] ≤ x[i] ≤ [9,2][i], Disjunct(Y[2])) @disjunction(m, inner, [W[1], W[2]], Disjunct(Y[1])) - @disjunction(m, outer, [Y[1], Y[2]]) + if use_complements + @disjunction(m, outer, [Y[1], Y[2]]) + else + #Infeasible disjunct + @constraint(m, y3[i=1:2], x[i] ≤ [-44,44][i], Disjunct(Y[3])) + @disjunction(m, outer, [Y[1], Y[2], Y[3]]) + end @test optimize!(m, gdp_method = BigM()) isa Nothing @test termination_status(m) == MOI.OPTIMAL @@ -98,6 +104,8 @@ function test_quadratic_gdp_example(use_complements = false) #psplit does not wo @objective(m, Max, sum(x)) @constraint(m, y1_quad, x[1]^2 + x[2]^2 ≤ 16, Disjunct(Y[1])) + # This constraint is always satisfied + @constraint(m, y1_global, x[1] + x[2] ≤ 20, Disjunct(Y[1])) @constraint(m, w1[i=1:2], [1, 2][i] ≤ x[i] ≤ [3, 4][i], Disjunct(W[1])) @constraint(m, w1_quad, x[1]^2 ≥ 2, Disjunct(W[1])) From d6afb11af1b446d88f127071bcaaf8b981e604cb Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 13 Jan 2026 13:20:08 -0500 Subject: [PATCH 09/13] tests for enhancements --- test/constraints/mbm.jl | 171 ++++++++++++++++++++++++++++++++++------ 1 file changed, 147 insertions(+), 24 deletions(-) diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index 5e4e492..2c4bcc9 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -67,28 +67,41 @@ function test_mini_model() model = GDPModel() @variable(model, 0 <= x, start = 1) @variable(model, 0 <= y) - @variable(model, Y[1:4], Logical) + @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])) - @disjunction(model, [Y[1], Y[2], Y[3], Y[4]]) + #infeasible constraint (x >= 100 but x <= 1 after bounds) + @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()) @test DP._mini_model(model, constraint_object(con), DisjunctConstraintRef[con2], mbm)== 0.0 set_upper_bound(x, 1) - @test DP._mini_model(model, constraint_object(con2), + @test DP._mini_model(model, constraint_object(con2), DisjunctConstraintRef[con], mbm)== 15 set_integer(y) @constraint(model, con3, y*x == 15, Disjunct(Y[1])) - @test DP._mini_model(model, constraint_object(con2), + @test DP._mini_model(model, constraint_object(con2), DisjunctConstraintRef[con], mbm)== 15 JuMP.fix(y, 5; force=true) - @test DP._mini_model(model, constraint_object(con2), + @test DP._mini_model(model, constraint_object(con2), DisjunctConstraintRef[con], mbm)== 10 + # With x <= 1 and y = 5, con2's region (3x + y >= 15) is infeasible: + # 3(1) + 5 = 8 < 15, so returns nothing (detected infeasibility) delete_lower_bound(x) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con2], mbm) == 1.0e9 + @test DP._mini_model(model, constraint_object(con2), + DisjunctConstraintRef[con2], mbm) == nothing + + # infeasible (x >= 100 but x has upper bound 1) + set_upper_bound(x, 1) + @test DP._mini_model( + model, + constraint_object(con), + DisjunctConstraintRef[truly_infeasible], + mbm + ) == nothing end function test_maximize_M() @@ -154,7 +167,6 @@ function test_maximize_M() DP._indicator_to_constraints(model)[Y[2]]), mbm) == [0.0, 0.0] - # zeros: x == 1 against Y[2] - per-row M values differ! # 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), @@ -166,6 +178,46 @@ function test_maximize_M() Vector{DisjunctConstraintRef}( 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() @@ -314,27 +366,25 @@ function test_reformulate_disjunction() 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) # 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 + # 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 + # - 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 + M_lower*Y[1] >= 0.0 → x + 0*Y[1] >= 0 - func_4 = ref_cons[4].func # x - M_upper*Y[1] <= 55 → x - 0*Y[1] <= 55 + 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 @@ -342,12 +392,85 @@ function test_reformulate_disjunction() @test JuMP.coefficient(func_2, x) == 1.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])) == 0.0 - @test JuMP.coefficient(func_4, x) == 1.0 - # -M_upper = -0 = 0 - @test JuMP.coefficient(func_4, binary_variable(Y[1])) == 0.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 @testset "MBM" begin From 0b5118f1f1a8d9e4ae8879216811a8b2c6392364 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 13 Jan 2026 13:32:52 -0500 Subject: [PATCH 10/13] change to typing --- src/mbm.jl | 12 +++--------- test/constraints/mbm.jl | 10 +++++++++- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 4d2dff4..3a64256 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -2,14 +2,8 @@ # HELPER FUNCTIONS ################################################################################ # Check if M result (scalar or vector) contains only zeros -function _is_all_zeros(M) - if M isa Number - return M == 0 - elseif M isa AbstractVector - return all(m == 0 for m in M) - end - return false -end +_is_all_zeros(M::Number) = iszero(M) +_is_all_zeros(M::AbstractVector) = all(iszero, M) ################################################################################ # CONSTRAINT, DISJUNCTION, DISJUNCT REFORMULATION @@ -462,7 +456,7 @@ function _mini_model( else M = JuMP.objective_value(sub_model) end - return max(M, 0) + return max(M, zero(M)) end ################################################################################ diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index 2c4bcc9..945b554 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -1,11 +1,19 @@ using HiGHS function test_mbm() - @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]) end function test__replace_variables_in_constraint() From 0fdac4824c584f020a5bac8ebc442b86ac5bf3a3 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Wed, 14 Jan 2026 09:47:32 -0500 Subject: [PATCH 11/13] clean up --- src/mbm.jl | 31 +++++++++++-------------------- 1 file changed, 11 insertions(+), 20 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 3a64256..0c9541c 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -15,33 +15,24 @@ function reformulate_disjunction( method::MBM ) mbm = _MBM(method, model) - ref_cons = Vector{JuMP.AbstractConstraint}() - # Track index ranges for each disjunct's constraints: (start, end) - disjunct_ranges = Dict{LogicalVariableRef, Tuple{Int, Int}}() + # Store constraints per disjunct, then flatten excluding deactivated + 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 ) - start_idx = length(ref_cons) + 1 - _reformulate_disjunct(model, ref_cons, d, mbm) - end_idx = length(ref_cons) - if end_idx >= start_idx - disjunct_ranges[d] = (start_idx, end_idx) - end - end - # Remove constraints from deactivated disjuncts (in reverse order) - indices_to_remove = Int[] - for deact in mbm.deactivated - if haskey(disjunct_ranges, deact) - start_idx, end_idx = disjunct_ranges[deact] - append!(indices_to_remove, start_idx:end_idx) - end + disjunct_cons[d] = Vector{JuMP.AbstractConstraint}() + _reformulate_disjunct(model, disjunct_cons[d], d, mbm) end - sort!(indices_to_remove, rev=true) - for idx in indices_to_remove - deleteat!(ref_cons, idx) + # Collect constraints from non-deactivated disjuncts + # it needs to be in a separate loop because disjuncts are only deactivated + # by looking 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]) end return ref_cons end From e60206353cd8e28136d3af1b545d18a192880257 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Wed, 14 Jan 2026 15:57:14 -0500 Subject: [PATCH 12/13] comments --- src/mbm.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/mbm.jl b/src/mbm.jl index 0c9541c..601271b 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -31,14 +31,16 @@ function reformulate_disjunction( # by looking reforming other disjuncts (subproblem infeasibility) ref_cons = Vector{JuMP.AbstractConstraint}() for d in disj.indicators + #Skip deactivated disjuncts d in mbm.deactivated && continue + #Skip disjuncts with no constraints haskey(disjunct_cons, d) && append!(ref_cons, disjunct_cons[d]) end return ref_cons end # Reformulates a disjunct represented by lvref using per-constraint M values. -# gets its own set of M_{ie,i'} values for each other term i'. +# Gets its own set of M_{ie,i'} values for each other term i'. function _reformulate_disjunct( model::JuMP.AbstractModel, ref_cons::Vector{JuMP.AbstractConstraint}, @@ -47,6 +49,7 @@ 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_conlvref = filter(d -> !(d in method.deactivated), method.conlvref) bconref = Dict(d => binary_variable(d) for d in active_conlvref) @@ -87,7 +90,8 @@ function _reformulate_disjunct( end con = JuMP.constraint_object(cref) - # Check if all M values are zero (constraint is global) + # 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)) @info "Constraint is global (M ≤ 0 for all disjuncts), " * @@ -112,9 +116,9 @@ function reformulate_disjunct_constraint( 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 From 4242a9d468d034582dbbdc78d6724edea02495d0 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Tue, 20 Jan 2026 13:11:57 -0500 Subject: [PATCH 13/13] store mbm problems --- src/datatypes.jl | 8 +- src/mbm.jl | 68 +++++++---- test/constraints/mbm.jl | 253 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 301 insertions(+), 28 deletions(-) diff --git a/src/datatypes.jl b/src/datatypes.jl index 5247d9a..0f9141c 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -392,13 +392,17 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho default_M::T conlvref::Vector{LogicalVariableRef{M}} deactivated::Set{LogicalVariableRef{M}} + # Stored submodels: indicator => (submodel, var_map) + store::Dict{LogicalVariableRef{M}, Tuple{M, Dict}} function _MBM(method::MBM{O, T}, model::M) where {O, T, M <: JuMP.AbstractModel} - new{O, T, M}(method.optimizer, + new{O, T, M}( + method.optimizer, Dict{LogicalVariableRef{M}, Union{T, Vector{T}}}(), method.default_M, Vector{LogicalVariableRef{M}}(), - Set{LogicalVariableRef{M}}() + Set{LogicalVariableRef{M}}(), + Dict{LogicalVariableRef{M}, Tuple{M, Dict}}() ) end end diff --git a/src/mbm.jl b/src/mbm.jl index 601271b..a515251 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -418,40 +418,62 @@ function _maximize_M( "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 +# Solve a mini-model to find the maximum value of the objective function for M. +# Uses stored submodels to avoid rebuilding from scratch each time. function _mini_model( - 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}} + indicator = _constraint_to_indicator(model)[first(constraints)] + + # Get or create stored submodel for this disjunct's feasible region + if !haskey(method.store, indicator) + method.store[indicator] = _create_submodel(model, constraints, method) + end + sub_model, var_map = method.store[indicator] + + # Set objective and solve + _constraint_to_objective(sub_model, objective, var_map) + JuMP.optimize!(sub_model) + status = JuMP.termination_status(sub_model) + + # Return M value + if status == _MOI.INFEASIBLE + return nothing + elseif status != _MOI.OPTIMAL || + !JuMP.has_values(sub_model) || + JuMP.primal_status(sub_model) != _MOI.FEASIBLE_POINT + return method.default_M + end + return max(JuMP.objective_value(sub_model), zero(method.default_M)) +end + +# Create a submodel for a disjunct's feasible region +function _create_submodel( + model::JuMP.AbstractModel, + constraints::Vector{<:DisjunctConstraintRef}, + method::_MBM +) var_type = JuMP.variable_ref_type(model) sub_model = _copy_model(model) - new_vars = Dict{var_type, var_type}() + var_map = Dict{var_type, var_type}() + for var in collect_all_vars(model) - new_vars[var] = variable_copy(sub_model, var) + var_map[var] = variable_copy(sub_model, var) end - for con in [JuMP.constraint_object(con) for con in constraints] - expr = _replace_variables_in_constraint(con.func, new_vars) + + for cref in constraints + con = JuMP.constraint_object(cref) + expr = _replace_variables_in_constraint(con.func, var_map) JuMP.@constraint(sub_model, expr * 1.0 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) - status = JuMP.termination_status(sub_model) - # Detect infeasibility: other disjunct has empty feasible region - if status == MOI.INFEASIBLE - return nothing - elseif status != 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 max(M, zero(M)) + + return (sub_model, var_map) end ################################################################################ diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index 945b554..14cd808 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -93,22 +93,26 @@ function test_mini_model() @constraint(model, con3, y*x == 15, Disjunct(Y[1])) @test DP._mini_model(model, constraint_object(con2), DisjunctConstraintRef[con], mbm)== 15 + # Create fresh _MBM after changing bounds (as reformulate_disjunction does) JuMP.fix(y, 5; force=true) + mbm2 = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con], mbm)== 10 + DisjunctConstraintRef[con], mbm2)== 10 # With x <= 1 and y = 5, con2's region (3x + y >= 15) is infeasible: # 3(1) + 5 = 8 < 15, so returns nothing (detected infeasibility) delete_lower_bound(x) + mbm3 = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con2], mbm) == nothing + DisjunctConstraintRef[con2], mbm3) == nothing # infeasible (x >= 100 but x has upper bound 1) set_upper_bound(x, 1) + mbm4 = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) @test DP._mini_model( model, constraint_object(con), DisjunctConstraintRef[truly_infeasible], - mbm + mbm4 ) == nothing end @@ -481,7 +485,250 @@ function test_reformulate_disjunction() @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_model, var_map = DP._create_submodel(model, constraints, mbm) + + # Check submodel was created + @test sub_model isa Model + + # Check variable mapping exists + @test haskey(var_map, x) + @test haskey(var_map, y) + + # Check mapped variables have correct bounds in submodel + @test has_lower_bound(var_map[x]) + @test lower_bound(var_map[x]) == 0 + @test has_upper_bound(var_map[x]) + @test upper_bound(var_map[x]) == 10 + + @test has_lower_bound(var_map[y]) + @test lower_bound(var_map[y]) == 0 + @test has_upper_bound(var_map[y]) + @test upper_bound(var_map[y]) == 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__replace_variables_in_constraint() test__constraint_to_objective()