diff --git a/src/datatypes.jl b/src/datatypes.jl index 6ba1315..0f9141c 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -388,15 +388,21 @@ 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}} + 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, - Dict{LogicalVariableRef{M}, T}(), + 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}}(), + Dict{LogicalVariableRef{M}, Tuple{M, Dict}}() ) end end diff --git a/src/mbm.jl b/src/mbm.jl index e9b1f13..a515251 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -1,54 +1,106 @@ +################################################################################ +# HELPER FUNCTIONS +################################################################################ +# Check if M result (scalar or vector) contains only zeros +_is_all_zeros(M::Number) = iszero(M) +_is_all_zeros(M::AbstractVector) = all(iszero, M) + ################################################################################ # 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}() + # 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 + ) + disjunct_cons[d] = Vector{JuMP.AbstractConstraint}() + _reformulate_disjunct(model, disjunct_cons[d], d, mbm) + end + # Collect constraints from non-deactivated disjuncts + # it needs to be in a separate loop because disjuncts are only deactivated + # by looking reforming other disjuncts (subproblem infeasibility) + ref_cons = Vector{JuMP.AbstractConstraint}() for d in disj.indicators - mbm.conlvref = filter(x -> x != d, disj.indicators) - _reformulate_disjunct(model, ref_cons, d, mbm) + #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 -#Reformualates a disjunct the disjunct of interest -#represented by lvref and the other indicators in conlvref + +# Reformulates a disjunct represented by lvref using per-constraint M values. +# Gets its own set of M_{ie,i'} values for each other 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) - + # Filter out deactivated disjuncts from binary variable mapping + # in the event we've identified some infeasible disjuncts already + active_conlvref = filter(d -> !(d in method.deactivated), method.conlvref) + bconref = Dict(d => binary_variable(d) for d in active_conlvref) + constraints = _indicator_to_constraints(model)[lvref] - filtered_constraints = [c for c in constraints if c isa DisjunctConstraintRef] - - for d in method.conlvref - d_constraints = _indicator_to_constraints(model)[d] - disjunct_constraints = [c for c in d_constraints if c isa DisjunctConstraintRef] - if !isempty(disjunct_constraints) - method.M[d] = maximum( - _maximize_M( - model, - JuMP.constraint_object(cref), + filtered_constraints = [ + c for c in constraints if c isa DisjunctConstraintRef + ] + + # For each constraint, compute its own set of M values + for cref in filtered_constraints + empty!(method.M) + + for d in method.conlvref + # Skip already-deactivated disjuncts + d in method.deactivated && continue + + d_constraints = _indicator_to_constraints(model)[d] + disjunct_constraints = [ + c for c in d_constraints if c isa DisjunctConstraintRef + ] + if !isempty(disjunct_constraints) + M_result = _maximize_M( + model, + JuMP.constraint_object(cref), disjunct_constraints, method - ) for cref in filtered_constraints - ) + ) + # 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) + # 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), " * + "adding without Big-M relaxation." + push!(ref_cons, con) + else + append!(ref_cons, reformulate_disjunct_constraint(model, con, + bconref, method)) end - end - for cref in filtered_constraints - con = JuMP.constraint_object(cref) - append!(ref_cons, reformulate_disjunct_constraint(model, con, - bconref, method)) end return ref_cons end @@ -64,13 +116,14 @@ 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 +# 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}, @@ -78,15 +131,14 @@ 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 + 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}, @@ -94,27 +146,26 @@ 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 + 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} - 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 + con.func[i] + sum(method.M[d][i] * bconref[d] for d in keys(method.M)) ) lower_expr = JuMP.@expression(model, [i=1:con.set.dimension], - con.func[i] - m_sum + con.func[i] - sum(method.M[d][i] * bconref[d] for d in keys(method.M)) ) upper_con = JuMP.build_constraint(error, upper_expr, MOI.Nonnegatives(con.set.dimension) @@ -153,6 +204,7 @@ function reformulate_disjunct_constraint( return [reform_con] end +# Uses per-bound M values: method.M[d][1] for lower, method.M[d][2] for upper function reformulate_disjunct_constraint( model::JuMP.AbstractModel, con::JuMP.ScalarConstraint{T, S}, @@ -160,21 +212,24 @@ 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)) + # 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)) ) - 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 +# 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}, @@ -182,21 +237,21 @@ 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) + # 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)) ) - 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 @@ -216,65 +271,80 @@ end ################################################################################ # Dispatches over constraint types to reformulate into >= or <= # in order to solve the mini-model +# 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} +) where {T, S <: _MOI.Nonpositives, R} val_type = JuMP.value_type(typeof(model)) - return maximum( + results = [ _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 - ) + ] + any(r === nothing for r in results) && return nothing + return results end +# 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} +) where {T, S <: _MOI.Nonnegatives, R} val_type = JuMP.value_type(typeof(model)) - return maximum( + results = [ _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 - ) + ] + any(r === nothing for r in results) && return nothing + return results end +# 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} +) where {T, S <: _MOI.Zeros, R} val_type = JuMP.value_type(typeof(model)) - return max( - maximum( - _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i],MOI.GreaterThan(zero(val_type))), - constraints, - method - ) for i in 1:objective.set.dimension - ), - maximum( - _maximize_M( - model, - JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), - constraints, - method - ) for i in 1:objective.set.dimension + results = [] + for i in 1:objective.set.dimension + M_ge = _maximize_M( + model, + JuMP.ScalarConstraint( + objective.func[i], MOI.GreaterThan(zero(val_type)) + ), + 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( @@ -286,50 +356,56 @@ function _maximize_M( return _mini_model(model, objective, constraints, method) end +# 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 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 ) + (M_lower === nothing || M_upper === nothing) && return nothing + return [M_lower, M_upper] end +# 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 max( - _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( @@ -342,36 +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) - if JuMP.termination_status(sub_model) != MOI.OPTIMAL || - !JuMP.has_values(sub_model) || - JuMP.primal_status(sub_model) != MOI.FEASIBLE_POINT - M = method.default_M - else - M = JuMP.objective_value(sub_model) - end - return M + + return (sub_model, var_map) end ################################################################################ diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index ed74a5a..14cd808 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -1,9 +1,19 @@ using HiGHS function test_mbm() + @test DP._MBM( + DP.MBM(HiGHS.Optimizer), JuMP.Model() + ).optimizer == HiGHS.Optimizer - @test DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()).optimizer == HiGHS.Optimizer - + # Test _is_all_zeros + @test DP._is_all_zeros(0) + @test DP._is_all_zeros(0.0) + @test !DP._is_all_zeros(1) + @test !DP._is_all_zeros(5.0) + @test DP._is_all_zeros([0, 0, 0]) + @test DP._is_all_zeros([0.0, 0.0]) + @test !DP._is_all_zeros([0, 1, 0]) + @test !DP._is_all_zeros([5.0, 0.0]) end function test__replace_variables_in_constraint() @@ -65,76 +75,161 @@ 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)== -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), + @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 + # Create fresh _MBM after changing bounds (as reformulate_disjunction does) JuMP.fix(y, 5; force=true) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con], mbm)== 10 + mbm2 = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + @test DP._mini_model(model, constraint_object(con2), + 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) - @test DP._mini_model(model, constraint_object(con2), - DisjunctConstraintRef[con2], mbm) == 1.0e9 + mbm3 = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + @test DP._mini_model(model, constraint_object(con2), + 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], + mbm4 + ) == nothing end function test_maximize_M() model = GDPModel() - @variable(model, 0 <= x[1:2] <= 50) + # Different bounds for x[1] and x[2] to demonstrate per-row M values + @variable(model, x[1:2]) + set_lower_bound(x[1], 0); set_upper_bound(x[1], 10) + set_lower_bound(x[2], 0); set_upper_bound(x[2], 5) @variable(model, Y[1:6], Logical) @constraint(model, lessthan, x[1] <= 1, Disjunct(Y[1])) @constraint(model, greaterthan, x[1] >= 1, Disjunct(Y[1])) @constraint(model, interval, 0 <= x[1] <= 55, Disjunct(Y[2])) @constraint(model, equalto, x[1] == 1, Disjunct(Y[3])) - @constraint(model, nonpositives, -x in MOI.Nonpositives(2), + # Vector constraints: x >= 0 (both rows) + @constraint(model, nonpositives, -x in MOI.Nonpositives(2), Disjunct(Y[4])) - @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), + @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), Disjunct(Y[5])) + # Vector equality: x == 1 (both rows) @constraint(model, zeros, -x .+ 1 in MOI.Zeros(2), Disjunct(Y[6])) mbm = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) - @test DP._maximize_M(model, constraint_object(interval), + + # Interval returns [M_lower, M_upper] + # M_lower = max(0 - x[1]) s.t. 0<=x[1]<=10 = 0 at x[1]=0 + # M_upper = max(x[1] - 55) s.t. 0<=x[1]<=10 = -45 at x[1]=10, clamped to 0 + @test DP._maximize_M(model, constraint_object(interval), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 0.0 - @test DP._maximize_M(model, constraint_object(lessthan), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [0.0, 0.0] + + # Scalar LessThan/GreaterThan still return scalars + # lessthan: x[1] <= 1 vs interval 0 <= x[1] <= 55 + # max(x[1] - 1) s.t. 0<=x[1]<=10 = 9 at x[1]=10 + @test DP._maximize_M(model, constraint_object(lessthan), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 49 - @test DP._maximize_M(model, constraint_object(greaterthan), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == 9.0 + + # greaterthan: x[1] >= 1 vs interval 0 <= x[1] <= 55 + # max(1 - x[1]) s.t. 0<=x[1]<=10 = 1 at x[1]=0 + @test DP._maximize_M(model, constraint_object(greaterthan), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), + DP._indicator_to_constraints(model)[Y[2]]), mbm) == 1.0 - @test DP._maximize_M(model, constraint_object(equalto), + + # EqualTo returns [M_lower, M_upper] + @test DP._maximize_M(model, constraint_object(equalto), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[3]]), - mbm) == 0 - @test DP._maximize_M(model, constraint_object(nonpositives), + DP._indicator_to_constraints(model)[Y[3]]), + mbm) == [0.0, 0.0] + + # Vector constraints: per-row M values + # nonpositives: x >= 0 against Y[2] (interval only on x[1]) + # Row 1: max(0 - x[1]) s.t. 0<=x[1]<=10 = 0 + # Row 2: max(0 - x[2]) s.t. 0<=x[2]<=5 = 0 + @test DP._maximize_M(model, constraint_object(nonpositives), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 0 - @test DP._maximize_M(model, constraint_object(nonnegatives), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [0.0, 0.0] + + @test DP._maximize_M(model, constraint_object(nonnegatives), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 0 - @test DP._maximize_M(model, constraint_object(zeros), + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [0.0, 0.0] + + # Row 1: max(|x[1] - 1|) s.t. 0<=x[1]<=10 = max(9, 1) = 9 + # Row 2: max(|x[2] - 1|) s.t. 0<=x[2]<=5 = max(4, 1) = 4 + @test DP._maximize_M(model, constraint_object(zeros), Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), - mbm) == 49 - @test_throws ErrorException DP._maximize_M(model, "odd", + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == [9.0, 4.0] + + @test_throws ErrorException DP._maximize_M(model, "odd", Vector{DisjunctConstraintRef}( - DP._indicator_to_constraints(model)[Y[2]]), + DP._indicator_to_constraints(model)[Y[2]]), mbm) + + # Add an infeasible disjunct (x >= 100 but bounds are 0-10) + @variable(model, Y_infeas, Logical) + @constraint(model, infeas_con, x[1] >= 100, Disjunct(Y_infeas)) + + # Scalar constraint against infeasible disjunct -> nothing + @test DP._maximize_M(model, constraint_object(lessthan), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Vector constraint (Nonpositives) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(nonpositives), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Vector constraint (Nonnegatives) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(nonnegatives), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Vector constraint (Zeros) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(zeros), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Bidirectional scalar (EqualTo) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(equalto), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing + + # Bidirectional scalar (Interval) against infeasible -> nothing + @test DP._maximize_M(model, constraint_object(interval), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y_infeas]), + mbm) == nothing end function test_reformulate_disjunct_constraint() @@ -144,57 +239,97 @@ function test_reformulate_disjunct_constraint() @constraint(model, lessthan, x[1] <= 1, Disjunct(Y[1])) @constraint(model, greaterthan, x[1] >= 1, Disjunct(Y[1])) @constraint(model, equalto, x[1] == 1, Disjunct(Y[2])) - @constraint(model, nonpositives, -x in MOI.Nonpositives(2), + @constraint(model, nonpositives, -x in MOI.Nonpositives(2), Disjunct(Y[3])) - @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), + @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), Disjunct(Y[4])) @constraint(model, zeros, -x .+ 1 in MOI.Zeros(2), Disjunct(Y[5])) @disjunction(model, disjunction,[Y[1], Y[2], Y[3], Y[4], Y[5]]) - method = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + bconref = Dict(Y[i] => binary_variable(Y[i]) for i in 1:5) + + # Test scalar constraints (LessThan, GreaterThan) with scalar M values + method_scalar = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) for i in 1:5 - method.M[Y[i]] = Float64(i) + method_scalar.M[Y[i]] = Float64(i) end - bconref = Dict(Y[i] => binary_variable(Y[i]) for i in 1:5) - reformulated_constraints = [reformulate_disjunct_constraint(model, - constraint_object(constraints), bconref, method) - for constraints in [lessthan, greaterthan, equalto, nonpositives, - nonnegatives, zeros, disjunction]] - @test reformulated_constraints[1][1].func == JuMP.@expression(model, - x[1] - sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[1][1].set == MOI.LessThan(1.0) - @test reformulated_constraints[2][1].func == JuMP.@expression(model, - x[1] + sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[2][1].set == MOI.GreaterThan(1.0) - @test reformulated_constraints[3][1].func == JuMP.@expression(model, - x[1] + sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[3][1].set == MOI.GreaterThan(1.0) - @test reformulated_constraints[3][2].func == JuMP.@expression(model, - x[1] - sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[3][2].set == MOI.LessThan(1.0) - @test reformulated_constraints[4][1].func == JuMP.@expression(model, - -x .- sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[4][1].set == MOI.Nonpositives(2) - @test reformulated_constraints[5][1].func == JuMP.@expression(model, - x .+ sum(method.M[i] * bconref[i] for i in keys(method.M))) && - reformulated_constraints[5][1].set == MOI.Nonnegatives(2) - @test reformulated_constraints[6][1].func == JuMP.@expression(model, - -x .+(1 + sum(method.M[i] * bconref[i] for i in keys(method.M)))) && - reformulated_constraints[6][1].set == MOI.Nonnegatives(2) - @test reformulated_constraints[6][2].func == JuMP.@expression(model, - -x .+(1 - sum(method.M[i] * bconref[i] for i in keys(method.M)))) && - reformulated_constraints[6][2].set == MOI.Nonpositives(2) - @test reformulated_constraints[7][1].func == JuMP.@expression(model, - x[1] - 52*bconref[Y[3]] - 53*bconref[Y[4]] - bconref[Y[1]] - - 5*bconref[Y[5]] - 2*bconref[Y[2]]) && - reformulated_constraints[7][1].set == MOI.LessThan(1.0) - @test reformulated_constraints[7][2].func == JuMP.@expression(model, - x[1] + 52*bconref[Y[3]] + 53*bconref[Y[4]] + bconref[Y[1]] - + 5*bconref[Y[5]] + 2*bconref[Y[2]]) && - reformulated_constraints[7][2].set == MOI.GreaterThan(1.0) - - @test_throws ErrorException reformulate_disjunct_constraint(model, - "odd", bconref, method) + ref_lessthan = reformulate_disjunct_constraint( + model, constraint_object(lessthan), bconref, method_scalar) + ref_greaterthan = reformulate_disjunct_constraint( + model, constraint_object(greaterthan), bconref, method_scalar) + @test ref_lessthan[1].func == JuMP.@expression(model, + x[1] - sum(method_scalar.M[i] * bconref[i] + for i in keys(method_scalar.M))) && + ref_lessthan[1].set == MOI.LessThan(1.0) + @test ref_greaterthan[1].func == JuMP.@expression(model, + x[1] + sum(method_scalar.M[i] * bconref[i] + for i in keys(method_scalar.M))) && + ref_greaterthan[1].set == MOI.GreaterThan(1.0) + + # Test bidirectional constraint (EqualTo) with [M_lower, M_upper] values + method_equalto = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + for i in 1:5 + method_equalto.M[Y[i]] = [Float64(i), Float64(i)] + end + ref_equalto = reformulate_disjunct_constraint( + model, constraint_object(equalto), bconref, method_equalto) + @test ref_equalto[1].func == JuMP.@expression(model, + x[1] + sum(method_equalto.M[i][1] * bconref[i] + for i in keys(method_equalto.M))) && + ref_equalto[1].set == MOI.GreaterThan(1.0) + @test ref_equalto[2].func == JuMP.@expression(model, + x[1] - sum(method_equalto.M[i][2] * bconref[i] + for i in keys(method_equalto.M))) && + ref_equalto[2].set == MOI.LessThan(1.0) + # Test vector constraints with per-row M values [M_row1, M_row2] + method_vector = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + for i in 1:5 + method_vector.M[Y[i]] = [Float64(i), Float64(i)] + end + ref_nonpositives = reformulate_disjunct_constraint( + model, constraint_object(nonpositives), bconref, method_vector) + ref_nonnegatives = reformulate_disjunct_constraint( + model, constraint_object(nonnegatives), bconref, method_vector) + ref_zeros = reformulate_disjunct_constraint( + model, constraint_object(zeros), bconref, method_vector) + @test ref_nonpositives[1].func == JuMP.@expression(model, [j=1:2], + -x[j] - sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_nonpositives[1].set == MOI.Nonpositives(2) + @test ref_nonnegatives[1].func == JuMP.@expression(model, [j=1:2], + x[j] + sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_nonnegatives[1].set == MOI.Nonnegatives(2) + @test ref_zeros[1].func == JuMP.@expression(model, [j=1:2], + -x[j] + 1 + sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_zeros[1].set == MOI.Nonnegatives(2) + @test ref_zeros[2].func == JuMP.@expression(model, [j=1:2], + -x[j] + 1 - sum(method_vector.M[i][j] * bconref[i] + for i in keys(method_vector.M))) && + ref_zeros[2].set == MOI.Nonpositives(2) + + # Test nested disjunction reformulation with proper nested structure + # Create outer disjunct with inner disjunction + model2 = GDPModel() + @variable(model2, 0 <= z <= 50) + @variable(model2, Outer[1:2], Logical) + @variable(model2, Inner[1:2], Logical) + @constraint(model2, inner_lt, z <= 1, Disjunct(Inner[1])) + @constraint(model2, inner_gt, z >= 1, Disjunct(Inner[2])) + @disjunction(model2, inner_disj, Inner) + method_nested = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + bconref2 = Dict(Outer[2] => binary_variable(Outer[2])) + method_nested.M[Outer[2]] = 10.0 #Dummy M value for testing. + #Normally _reformulate_disjunct will this without having to assign a value + ref_disjunction = reformulate_disjunct_constraint( + model2, constraint_object(inner_disj), bconref2, method_nested) + @test length(ref_disjunction) >= 2 + @test JuMP.coefficient(ref_disjunction[1].func, z) == 1.0 + @test JuMP.coefficient(ref_disjunction[2].func, z) == 1.0 + + @test_throws ErrorException reformulate_disjunct_constraint(model, + "odd", bconref, method_scalar) end function test_reformulate_disjunct() @@ -221,13 +356,14 @@ function test_reformulate_disjunct() func_3 = reformulated_disjunct[3].func @test JuMP.coefficient(func_1, x[1]) == 1.0 - @test JuMP.coefficient(func_1, binary_variable(Y[2])) == -1.5 + @test JuMP.coefficient(func_1, binary_variable(Y[2])) == 0.0 + # Per-bound M: lower bound uses M_lower=1.5, upper bound uses M_upper=2.5 @test JuMP.coefficient(func_2, x[1]) == 1.0 - @test JuMP.coefficient(func_2, binary_variable(Y[1])) == 2.5 + @test JuMP.coefficient(func_2, binary_variable(Y[1])) == 1.5 # M_lower @test JuMP.coefficient(func_3, x[1]) == 1.0 - @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -2.5 + @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -2.5 # -M_upper end function test_reformulate_disjunction() @@ -238,39 +374,361 @@ function test_reformulate_disjunction() @constraint(model, greaterthan, x >= 1, Disjunct(Y[1])) @constraint(model, interval, 0 <= x <= 55, Disjunct(Y[2])) disj = disjunction(model, [Y[1], Y[2]]) - + method = DP.MBM(HiGHS.Optimizer) ref_cons = reformulate_disjunction(model, constraint_object(disj), method) - @test length(ref_cons) == 4 + # 3 constraints: lessthan, greaterthan (with Big-M), interval (global) + @test length(ref_cons) == 3 @test ref_cons[1].set == MOI.LessThan(2.0) - @test ref_cons[2].set == MOI.GreaterThan(1.0) - - @test ref_cons[3].set == MOI.GreaterThan(0.0) - - @test ref_cons[4].set == MOI.LessThan(55.0) + # Interval is global (M=0 for both bounds in Y[1]'s region 1<=x<=2) + @test ref_cons[3].set == MOI.Interval(0.0, 55.0) - func_1 = ref_cons[1].func # x - 53 Y[2] <= 2.0 - func_2 = ref_cons[2].func # x + 53 Y[2] >= 1.0 - func_3 = ref_cons[3].func # x - Y[1] >= 0.0 - func_4 = ref_cons[4].func # x + Y[1] <= 55.0 + # Per-constraint, per-bound M values: + # - lessthan (x <= 2) in Y[2] region (0 <= x <= 55): max(x-2) at + # x=55 -> M=53 + # - greaterthan (x >= 1) in Y[2] region: max(1-x) at x=0 -> M=1 + # - interval in Y[1] region (1 <= x <= 2): + # - M_lower (x >= 0): max(0-x) at x=1 -> M_lower=-1, clamped to 0 + # - M_upper (x <= 55): max(x-55) at x=2 -> M_upper=-53, clamped to 0 + # - Both M=0 -> detected as global, added without Big-M + func_1 = ref_cons[1].func # x - 53*Y[2] <= 2.0 + func_2 = ref_cons[2].func # x + 1*Y[2] >= 1.0 + func_3 = ref_cons[3].func # x (global, no binary variables) @test JuMP.coefficient(func_1, x) == 1.0 @test JuMP.coefficient(func_1, binary_variable(Y[2])) == -53.0 @test JuMP.coefficient(func_2, x) == 1.0 - @test JuMP.coefficient(func_2, binary_variable(Y[2])) == 53.0 + @test JuMP.coefficient(func_2, binary_variable(Y[2])) == 1.0 + # Global constraint has just x, no binary variables @test JuMP.coefficient(func_3, x) == 1.0 - @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -1.0 - @test JuMP.coefficient(func_4, x) == 1.0 - @test JuMP.coefficient(func_4, binary_variable(Y[1])) == 1.0 + #Test infeasible disjunct detection and deactivation + model2 = GDPModel() + @variable(model2, 0 <= z <= 1) + @variable(model2, W[1:2], Logical) + # W[1]: z >= 5 is infeasible (z has upper bound 1) + @constraint(model2, z >= 5, Disjunct(W[1])) + # W[2]: z <= 0.5 is feasible + @constraint(model2, z <= 0.5, Disjunct(W[2])) + disj2 = disjunction(model2, [W[1], W[2]]) + + method2 = DP.MBM(HiGHS.Optimizer) + ref_cons2 = @test_logs (:warn, r"infeasible, deactivating") begin + reformulate_disjunction(model2, constraint_object(disj2), method2) + end + + @test length(ref_cons2) == 1 + @test ref_cons2[1].set == MOI.LessThan(0.5) + + #Test multiple infeasible disjuncts + model3 = GDPModel() + @variable(model3, 0 <= w <= 1) + @variable(model3, V[1:3], Logical) + @constraint(model3, w >= 5, Disjunct(V[1])) # infeasible + @constraint(model3, w >= 10, Disjunct(V[2])) # infeasible + @constraint(model3, w <= 0.5, Disjunct(V[3])) # feasible + disj3 = disjunction(model3, [V[1], V[2], V[3]]) + + method3 = DP.MBM(HiGHS.Optimizer) + #warn about V[1] and V[2] being infeasible + ref_cons3 = @test_logs (:warn,) (:warn,) begin + reformulate_disjunction(model3, constraint_object(disj3), method3) + end + + # Only V[3]'s constraint should be reformulated + @test length(ref_cons3) == 1 + @test ref_cons3[1].set == MOI.LessThan(0.5) + + model4 = GDPModel() + @variable(model4, 0 <= u <= 10) + @variable(model4, U[1:2], Logical) + @constraint(model4, u <= 3, Disjunct(U[1])) + @constraint(model4, u >= 5, Disjunct(U[2])) + disj4 = disjunction(model4, [U[1], U[2]]) + + method4 = DP.MBM(HiGHS.Optimizer) + ref_cons4 = @test_nowarn begin + reformulate_disjunction(model4, constraint_object(disj4), method4) + end + @test length(ref_cons4) == 2 + + # Disjunct 1: x <= 10 (will be global because D2's region has x <= 5) + # Disjunct 2: x <= 5 + # max(x - 10) s.t. x <= 5 = 5 - 10 = -5, clamped to 0 + # So M = 0 for D1's constraint -> it's global + model5 = GDPModel() + @variable(model5, 0 <= g <= 10) + @variable(model5, G[1:2], Logical) + @constraint(model5, g <= 10, Disjunct(G[1])) # global: other region is g<=5 + @constraint(model5, g <= 5, Disjunct(G[2])) + disj5 = disjunction(model5, [G[1], G[2]]) + + method5 = DP.MBM(HiGHS.Optimizer) + ref_cons5 = reformulate_disjunction(model5, constraint_object(disj5), method5) + # G[1]'s constraint is global (added without Big-M) + # G[2]'s constraint has M > 0 (needs Big-M): max(g - 5) s.t. g<=10 = 5 + @test length(ref_cons5) == 2 + # Check that one constraint has binary variable coefficient (Big-M term) + # and one doesn't (global) + global_con = ref_cons5[1] # g <= 10 (global) + bigm_con = ref_cons5[2] # g <= 5 with Big-M + @test global_con.set == MOI.LessThan(10.0) + @test bigm_con.set == MOI.LessThan(5.0) + # Global constraint should have no binary variable coefficient + @test JuMP.coefficient(global_con.func, binary_variable(G[2])) == 0.0 + # Big-M constraint should have binary variable coefficient = -5 + @test JuMP.coefficient(bigm_con.func, binary_variable(G[1])) == -5.0 +end + +################################################################################ +# LOW-LEVEL UNIT TESTS FOR MBM INFRASTRUCTURE +################################################################################ + +# Test _copy_model directly +function test__copy_model() + # Test with standard JuMP Model + model = Model() + @variable(model, x >= 0) + @variable(model, y <= 10) + copied = DP._copy_model(model) + @test copied isa Model + @test num_variables(copied) == 0 # _copy_model creates empty model + + # Test with GDPModel + gdp_model = GDPModel() + @variable(gdp_model, z) + copied_gdp = DP._copy_model(gdp_model) + @test copied_gdp isa Model + @test num_variables(copied_gdp) == 0 +end + +# Test VariableProperties struct and constructors +function test_variable_properties() + # Test VariableProperties(vref::GenericVariableRef) - standard JuMP variable + model = GDPModel() + @variable(model, 0 <= x <= 10, start = 5) + props_x = DP.VariableProperties(x) + @test props_x.info.has_lb == true + @test props_x.info.lower_bound == 0 + @test props_x.info.has_ub == true + @test props_x.info.upper_bound == 10 + @test props_x.info.has_start == true + @test props_x.info.start == 5 + @test props_x.name == "x" + @test props_x.variable_type === nothing + + # Test with binary variable + @variable(model, y, Bin) + props_y = DP.VariableProperties(y) + @test props_y.info.binary == true + @test props_y.info.integer == false + + # Test with integer variable + @variable(model, z, Int) + props_z = DP.VariableProperties(z) + @test props_z.info.binary == false + @test props_z.info.integer == true + + # Test with fixed variable + @variable(model, w == 42) + props_w = DP.VariableProperties(w) + @test props_w.info.has_fix == true + @test props_w.info.fixed_value == 42 + + # Test VariableProperties(expr) - blank info constructor + expr = 2*x + 3*y + props_expr = DP.VariableProperties(expr) + @test props_expr.info.has_lb == false + @test props_expr.info.has_ub == false + @test props_expr.info.has_fix == false + @test props_expr.info.has_start == false + @test props_expr.info.binary == false + @test props_expr.info.integer == false + @test props_expr.name == "" +end + +# Test _make_variable_object +function test__make_variable_object() + model = GDPModel() + @variable(model, 0 <= x <= 10, start = 5) + + # Create VariableProperties and then make variable object + props = DP.VariableProperties(x) + var_obj = DP._make_variable_object(props) + @test var_obj isa JuMP.ScalarVariable + @test var_obj.info.has_lb == true + @test var_obj.info.lower_bound == 0 + @test var_obj.info.has_ub == true + @test var_obj.info.upper_bound == 10 + + # Test with blank properties (from expression) + props_blank = DP.VariableProperties(2*x) + var_obj_blank = DP._make_variable_object(props_blank) + @test var_obj_blank isa JuMP.ScalarVariable + @test var_obj_blank.info.has_lb == false + @test var_obj_blank.info.has_ub == false +end + +# Test create_variable +function test_create_variable() + model = GDPModel() + @variable(model, 0 <= x <= 10, start = 5) + + # Create a new variable from properties + props = DP.VariableProperties(x) + new_var = DP.create_variable(model, props) + @test new_var isa VariableRef + @test has_lower_bound(new_var) + @test lower_bound(new_var) == 0 + @test has_upper_bound(new_var) + @test upper_bound(new_var) == 10 + @test start_value(new_var) == 5 + @test name(new_var) == "x" + + # Test with binary variable + @variable(model, y, Bin) + props_bin = DP.VariableProperties(y) + new_bin = DP.create_variable(model, props_bin) + @test is_binary(new_bin) + + # Test with integer variable + @variable(model, z, Int) + props_int = DP.VariableProperties(z) + new_int = DP.create_variable(model, props_int) + @test is_integer(new_int) + + # Test with fixed variable + @variable(model, w == 42) + props_fix = DP.VariableProperties(w) + new_fix = DP.create_variable(model, props_fix) + @test is_fixed(new_fix) + @test fix_value(new_fix) == 42 +end + +# Test variable_copy +function test_variable_copy() + source_model = GDPModel() + @variable(source_model, 0 <= x <= 10, start = 5) + @variable(source_model, y, Bin) + @variable(source_model, z, Int) + @variable(source_model, w == 42) + + target_model = GDPModel() + + # Copy bounded variable + x_copy = DP.variable_copy(target_model, x) + @test x_copy isa VariableRef + @test owner_model(x_copy) === target_model + @test has_lower_bound(x_copy) + @test lower_bound(x_copy) == 0 + @test has_upper_bound(x_copy) + @test upper_bound(x_copy) == 10 + @test start_value(x_copy) == 5 + @test name(x_copy) == "x" + + # Copy binary variable + y_copy = DP.variable_copy(target_model, y) + @test is_binary(y_copy) + @test owner_model(y_copy) === target_model + + # Copy integer variable + z_copy = DP.variable_copy(target_model, z) + @test is_integer(z_copy) + @test owner_model(z_copy) === target_model + + # Copy fixed variable + w_copy = DP.variable_copy(target_model, w) + @test is_fixed(w_copy) + @test fix_value(w_copy) == 42 + @test owner_model(w_copy) === target_model +end + +# Test _create_submodel (combines _copy_model + variable_copy + constraint handling) +function test__create_submodel() + model = GDPModel() + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 5) + @variable(model, Y[1:2], Logical) + @constraint(model, con1, x + y <= 8, Disjunct(Y[1])) + @constraint(model, con2, x - y >= 2, Disjunct(Y[2])) + @disjunction(model, [Y[1], Y[2]]) + + mbm = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + constraints = Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y[1]] + ) + + sub_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() 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]))