diff --git a/README.md b/README.md index 80868bd1..09ffc22c 100644 --- a/README.md +++ b/README.md @@ -172,6 +172,11 @@ The following reformulation methods are currently supported: 3. [Indicator](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints): This method reformulates each disjunct constraint into an indicator constraint with the Boolean reformulation counterpart of the Logical variable used to define the disjunct constraint. +4. [MBM](https://doi.org/10.1016/j.compchemeng.2015.02.013): The multiple big-m method creates multiple M values for each disjunct constraint. The 'MBM' struct is created with the following required argument: + + - `optimizer`: Optimizer to use when solving subproblems to determine M values. This is a required value. + - `default_M`: Default big-M value to use if no big-M is specified for a logical variable (1e9). + ## Release Notes Prior to `v0.4.0`, the package did not leverage the JuMP extension capabilities and was not as robust. For these earlier releases, refer to [Perez, Joshi, and Grossmann, 2023](https://arxiv.org/abs/2304.10492v1) and the following [JuliaCon 2022 Talk](https://www.youtube.com/watch?v=AMIrgTTfUkI). diff --git a/src/DisjunctiveProgramming.jl b/src/DisjunctiveProgramming.jl index ae9892e0..5e3b1623 100644 --- a/src/DisjunctiveProgramming.jl +++ b/src/DisjunctiveProgramming.jl @@ -6,7 +6,6 @@ Reexport.@reexport using JuMP # Use Meta for metaprogramming using Base.Meta - # Create aliases import JuMP.MOI as _MOI import JuMP.MOIU.CleverDicts as _MOIUC @@ -22,8 +21,10 @@ include("macros.jl") include("reformulate.jl") include("bigm.jl") include("hull.jl") +include("mbm.jl") include("indicator.jl") include("print.jl") +include("utilities.jl") # Define additional stuff that should not be exported const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include] diff --git a/src/datatypes.jl b/src/datatypes.jl index 25291142..22480181 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -367,6 +367,40 @@ struct BigM{T} <: AbstractReformulationMethod end end +""" + MBM{O, T, L <: LogicalVariableRef} <: AbstractReformulationMethod + +A type for using the multiple big-M reformulation approach for disjunctive constraints. + +**Fields** +- `optimizer::O`: Optimizer to use when solving mini-models (required). +- `default_M::T`: Default big-M value to use if no big-M is specified for a logical variable (1e9). +""" +mutable struct MBM{O, T} <: AbstractReformulationMethod + optimizer::O + default_M::T + + # Constructor with optimizer (required) and optional default_M + function MBM(optimizer::O, default_M::T = 1e9) where {O, T} + new{O, T}(optimizer, default_M) + end +end + +mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMethod + optimizer::O + M::Dict{LogicalVariableRef{M}, 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}(), + method.default_M, + Vector{LogicalVariableRef{M}}() + ) + end +end + """ Hull{T} <: AbstractReformulationMethod diff --git a/src/mbm.jl b/src/mbm.jl new file mode 100644 index 00000000..58120a8f --- /dev/null +++ b/src/mbm.jl @@ -0,0 +1,477 @@ +################################################################################ +# CONSTRAINT, DISJUNCTION, DISJUNCT REFORMULATION +################################################################################ +#Reformulates the disjunction using multiple big-M values +function reformulate_disjunction( + model::JuMP.AbstractModel, + disj::Disjunction, + method::MBM +) + mbm = _MBM(method, model) + 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) + end + return ref_cons +end +#Reformualates a disjunct the disjunct of interest +#represented by lvref and the other indicators in conlvref +function _reformulate_disjunct( + 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), + disjunct_constraints, + method + ) for cref in filtered_constraints + ) + end + end + for cref in filtered_constraints + con = JuMP.constraint_object(cref) + append!(ref_cons, reformulate_disjunct_constraint(model, con, + bconref, method)) + end + return ref_cons +end + +function reformulate_disjunct_constraint( + model::JuMP.AbstractModel, + con::Disjunction, + bconref::Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + method::_MBM +) + + ref_cons = reformulate_disjunction(model, con, MBM(method.optimizer)) + new_ref_cons = Vector{JuMP.AbstractConstraint}() + for ref_con in ref_cons + append!(new_ref_cons, + reformulate_disjunct_constraint(model, ref_con, bconref, method) + ) + end + return new_ref_cons +end + +function reformulate_disjunct_constraint( + model::JuMP.AbstractModel, + con::JuMP.VectorConstraint{T, S, R}, + bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + method::_MBM +) where {T, S <: _MOI.Nonpositives, R} + m_sum = sum(method.M[i] * bconref[i] for i in keys(method.M)) + new_func = JuMP.@expression(model, [i=1:con.set.dimension], + con.func[i] - m_sum + ) + 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}, + bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + method::_MBM +) where {T, S <: _MOI.Nonnegatives, R} + m_sum = sum(method.M[i] * bconref[i] for i in keys(method.M)) + new_func = JuMP.@expression(model, [i=1:con.set.dimension], + con.func[i] + m_sum + ) + 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}, + 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 + ) + lower_expr = JuMP.@expression(model, [i=1:con.set.dimension], + con.func[i] - m_sum + ) + upper_con = JuMP.build_constraint(error, upper_expr, + MOI.Nonnegatives(con.set.dimension) + ) + lower_con = JuMP.build_constraint(error, lower_expr, + MOI.Nonpositives(con.set.dimension) + ) + return [upper_con, lower_con] +end + + +function reformulate_disjunct_constraint( + model::JuMP.AbstractModel, + con::JuMP.ScalarConstraint{T, S}, + bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + method::_MBM +) where {T, S <: _MOI.LessThan} + new_func = JuMP.@expression(model, + con.func - sum(method.M[i] * bconref[i] for i in keys(method.M))) + reform_con = JuMP.build_constraint(error, new_func, con.set) + return [reform_con] +end + +function reformulate_disjunct_constraint( + model::JuMP.AbstractModel, + con::JuMP.ScalarConstraint{T, S}, + bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + method::_MBM +) where {T, S <: _MOI.GreaterThan} + new_func = JuMP.@expression(model, + con.func + sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) + reform_con = JuMP.build_constraint(error, new_func, con.set) + return [reform_con] +end + +function reformulate_disjunct_constraint( + model::JuMP.AbstractModel, + con::JuMP.ScalarConstraint{T, S}, + bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + method::_MBM +) where {T, S <: _MOI.EqualTo} + upper_func = JuMP.@expression(model, + con.func - sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) + lower_func = JuMP.@expression(model, + con.func + sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) + upper_con = JuMP.build_constraint(error, upper_func, + MOI.LessThan(con.set.value) + ) + lower_con = JuMP.build_constraint(error, lower_func, + MOI.GreaterThan(con.set.value) + ) + return [lower_con, upper_con] +end + +function reformulate_disjunct_constraint( + model::JuMP.AbstractModel, + con::JuMP.ScalarConstraint{T, S}, + bconref:: Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + method::_MBM +) where {T, S <: _MOI.Interval} + set_values = _set_values(con.set) + upper_func = JuMP.@expression(model, + con.func - sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) + upper_con = JuMP.build_constraint(error, upper_func, + MOI.LessThan(set_values[2]) + ) + + lower_func = JuMP.@expression(model, + con.func + sum(method.M[i] * bconref[i] for i in keys(method.M)) + ) + lower_con = JuMP.build_constraint(error, lower_func, + MOI.GreaterThan(set_values[1]) + ) + + return [lower_con, upper_con] +end + +function reformulate_disjunct_constraint( + ::JuMP.AbstractModel, + ::F, + ::Union{Dict{<:LogicalVariableRef,<:JuMP.AbstractVariableRef}, + Dict{<:LogicalVariableRef,<:JuMP.GenericAffExpr}}, + ::_MBM +) where {F} + error("Constraint type $(typeof(F)) is not supported by the " * + "Multiple Big-M reformulation method.") +end + +################################################################################ +# MULTIPLE BIG-M REFORMULATION +################################################################################ +# Dispatches over constraint types to reformulate into >= or <= +# in order to solve the mini-model +function _maximize_M( + model::JuMP.AbstractModel, + objective::JuMP.VectorConstraint{T, S, R}, + constraints::Vector{<:DisjunctConstraintRef}, + method::_MBM +) where { T, S <: _MOI.Nonpositives, R} + val_type = JuMP.value_type(typeof(model)) + return maximum( + _maximize_M( + model, + JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), + constraints, + method + ) for i in 1:objective.set.dimension + ) +end + +function _maximize_M( + 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( + _maximize_M( + model, + JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(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}, + method::_MBM +) where { T, S <: _MOI.Zeros, R} + val_type = JuMP.value_type(typeof(model)) + return max( + maximum( + _maximize_M( + model, + JuMP.ScalarConstraint(objective.func[i],MOI.GreaterThan(zero(val_type))), + constraints, + method + ) for i in 1:objective.set.dimension + ), + maximum( + _maximize_M( + model, + JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(zero(val_type))), + constraints, + method + ) for i in 1:objective.set.dimension + ) + ) +end + +function _maximize_M( + model::JuMP.AbstractModel, + objective::JuMP.ScalarConstraint{T, S}, + constraints::Vector{<:DisjunctConstraintRef}, + method::_MBM +) where {T, S <: Union{_MOI.LessThan, _MOI.GreaterThan}} + return _mini_model(model, objective, constraints, method) +end + +function _maximize_M( + 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 + ) + ) +end + +function _maximize_M( + 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 + ) + ) +end + +function _maximize_M( + ::JuMP.AbstractModel, + ::F, + ::Vector{<:DisjunctConstraintRef}, + ::_MBM +) where {F} + error("This type of constraints and objective constraint has " * + "not been implemented for MBM subproblems\nF: $(F)") +end + +# Solve a mini-model to find the maximum value of the objective +# function for M value +function _mini_model( + model::JuMP.AbstractModel, + objective::JuMP.ScalarConstraint{T,S}, + constraints::Vector{<:DisjunctConstraintRef}, + method::_MBM +) where {T,S <: Union{_MOI.LessThan, _MOI.GreaterThan}} + var_type = JuMP.variable_ref_type(model) + sub_model = _copy_model(model) + new_vars = Dict{var_type, var_type}() + for var in JuMP.all_variables(model) + new_vars[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) + 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 +end + +################################################################################ +# CONSTRAINT TO OBJECTIVE +################################################################################ +function _constraint_to_objective( + sub_model::JuMP.AbstractModel, + obj::JuMP.ScalarConstraint{<:JuMP.AbstractJuMPScalar, MOI.LessThan{T}}, + new_vars::Dict{V,K} +) where {T,V <: JuMP.AbstractVariableRef, K <: JuMP.AbstractVariableRef} + JuMP.@objective(sub_model, Max, + - obj.set.upper + _replace_variables_in_constraint(obj.func, new_vars) + ) +end +function _constraint_to_objective( + sub_model::JuMP.AbstractModel, + obj::JuMP.ScalarConstraint{<:JuMP.AbstractJuMPScalar, MOI.GreaterThan{T}}, + new_vars::Dict{V,K} +) where {T,V <: JuMP.AbstractVariableRef, K <: JuMP.AbstractVariableRef} + JuMP.@objective(sub_model, Max, + - _replace_variables_in_constraint(obj.func, new_vars) + obj.set.lower + ) +end + +function _constraint_to_objective( + sub_model::JuMP.AbstractModel, + obj::JuMP.ScalarConstraint, + new_vars::Dict{V,K} +) where {V, K} + error("This type of constraint is not supported, only greater " * + "than and less than constraints are supported with " * + "intervals and equalities being converted.") +end + +################################################################################ +# REPLACE VARIABLES IN CONSTRAINT +################################################################################ + +function _replace_variables_in_constraint( + fun:: JuMP.AbstractVariableRef, + var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} +) + return var_map[fun] +end + +function _replace_variables_in_constraint( + fun::T, + var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} +) where {T <: JuMP.GenericAffExpr} + new_aff = JuMP.zero(T) + for (var, coef) in fun.terms + new_var = var_map[var] + JuMP.add_to_expression!(new_aff, coef, new_var) + end + new_aff.constant = fun.constant + return new_aff +end + +function _replace_variables_in_constraint( + fun::T, + var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} +) where {T <: JuMP.GenericQuadExpr} + new_quad = JuMP.zero(T) + for (vars, coef) in fun.terms + JuMP.add_to_expression!(new_quad, coef, + var_map[vars.a], var_map[vars.b]) + end + new_aff = _replace_variables_in_constraint(fun.aff, var_map) + JuMP.add_to_expression!(new_quad, new_aff) + return new_quad +end + +function _replace_variables_in_constraint( + fun::Number, + var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} +) + return fun +end + +function _replace_variables_in_constraint( + fun::T, + var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} +) where {T <: JuMP.GenericNonlinearExpr} + new_args = Any[_replace_variables_in_constraint(arg, var_map) + for arg in fun.args] + return T(fun.head, new_args) +end + +function _replace_variables_in_constraint( + fun::Vector, + var_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} +) + return [_replace_variables_in_constraint(expr, var_map) + for expr in fun] +end + +function _replace_variables_in_constraint( + ::F, + ::S +) where {F, S} + error("_replace_variables_in_constraint not implemented for " * + "$(typeof(F)) and $(typeof(S))") +end diff --git a/src/utilities.jl b/src/utilities.jl new file mode 100644 index 00000000..c3aa6f16 --- /dev/null +++ b/src/utilities.jl @@ -0,0 +1,9 @@ +################################################################################ +# MODEL COPYING +################################################################################ +# extentsion point for model copying +function _copy_model( + model::M + ) where {M <: JuMP.AbstractModel} + return M() +end \ No newline at end of file diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl new file mode 100644 index 00000000..ed74a5a5 --- /dev/null +++ b/test/constraints/mbm.jl @@ -0,0 +1,282 @@ +using HiGHS + +function test_mbm() + + @test DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()).optimizer == HiGHS.Optimizer + +end + +function test__replace_variables_in_constraint() + model = Model() + sub_model = Model() + @variable(model, x[1:3]) + @constraint(model, con1,x[1] <= 1) + @constraint(model, con2,x[2]*x[1] <= 1) + @constraint(model, con3,sin(x[3]) <= 0) + @constraint(model, con4, [x[1],x[2],x[3]] in MOI.Zeros(3)) + + #Test GenericVariableRef + new_vars = Dict{AbstractVariableRef, AbstractVariableRef}() + [new_vars[x[i]] = @variable(sub_model) for i in 1:3] + varref = DP._replace_variables_in_constraint(x[1], new_vars) + expr1 = DP._replace_variables_in_constraint( + constraint_object(con1).func, new_vars) + expr2 = DP._replace_variables_in_constraint( + constraint_object(con2).func, new_vars) + expr3 = DP._replace_variables_in_constraint( + constraint_object(con3).func, new_vars) + expr4 = DP._replace_variables_in_constraint( + constraint_object(con4).func, new_vars) + @test expr1 == JuMP.@expression(sub_model, new_vars[x[1]] + 1 - 1) + @test expr2 == JuMP.@expression(sub_model, new_vars[x[2]]*new_vars[x[1]]) + @test varref == new_vars[x[1]] + @test expr3 isa JuMP.NonlinearExpr + expected = JuMP.@expression(sub_model, sin(new_vars[x[3]]) - 0.0) + @test JuMP.isequal_canonical(expr3, expected) + @test expr4 == [new_vars[x[i]] for i in 1:3] + @test_throws ErrorException DP._replace_variables_in_constraint( + "String", new_vars) +end + +function test__constraint_to_objective() + model = Model() + sub_model = Model() + + @variable(model, x[1:2]) + @constraint(model, lessthan, x[1] <= 1) + @constraint(model, greaterthan, x[2] >= 1) + @constraint(model, interval, 0 <= x[1] <= 55) + @constraint(model, equalto, x[1] == 1) + new_vars = Dict{AbstractVariableRef, AbstractVariableRef}() + [new_vars[x[i]] = @variable(sub_model) for i in 1:2] + DP._constraint_to_objective(sub_model, constraint_object(lessthan), + new_vars) + @test objective_function(sub_model) == JuMP.@expression(sub_model, + new_vars[x[1]] - 1) + DP._constraint_to_objective(sub_model, constraint_object(greaterthan), + new_vars) + @test objective_function(sub_model) == JuMP.@expression(sub_model, + 1 - new_vars[x[2]]) + @test_throws ErrorException DP._constraint_to_objective(sub_model, + constraint_object(interval), new_vars) +end + +function test_mini_model() + model = GDPModel() + @variable(model, 0 <= x, start = 1) + @variable(model, 0 <= y) + @variable(model, Y[1:4], Logical) + @constraint(model, con, 3*-x <= 4, Disjunct(Y[1])) + @constraint(model, con2, 3*x + y >= 15, Disjunct(Y[2])) + @constraint(model, infeasiblecon, 3*x + y == 15, Disjunct(Y[3])) + @constraint(model, intervalcon, 0 <= x <= 55, Disjunct(Y[4])) + @disjunction(model, [Y[1], Y[2], Y[3], Y[4]]) + mbm = DP._MBM(DP.MBM(HiGHS.Optimizer), JuMP.Model()) + @test DP._mini_model(model, constraint_object(con), + DisjunctConstraintRef[con2], mbm)== -4 + set_upper_bound(x, 1) + @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), + DisjunctConstraintRef[con], mbm)== 15 + JuMP.fix(y, 5; force=true) + @test DP._mini_model(model, constraint_object(con2), + DisjunctConstraintRef[con], mbm)== 10 + delete_lower_bound(x) + @test DP._mini_model(model, constraint_object(con2), + DisjunctConstraintRef[con2], mbm) == 1.0e9 +end + +function test_maximize_M() + model = GDPModel() + @variable(model, 0 <= x[1:2] <= 50) + @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), + Disjunct(Y[4])) + @constraint(model, nonnegatives, x in MOI.Nonnegatives(2), + Disjunct(Y[5])) + @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), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == 0.0 + @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), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == 1.0 + @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), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == 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), + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y[2]]), + mbm) == 49 + @test_throws ErrorException DP._maximize_M(model, "odd", + Vector{DisjunctConstraintRef}( + DP._indicator_to_constraints(model)[Y[2]]), + mbm) +end + +function test_reformulate_disjunct_constraint() + model = GDPModel() + @variable(model, 0 <= x[1:2] <= 50) + @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, equalto, x[1] == 1, Disjunct(Y[2])) + @constraint(model, nonpositives, -x in MOI.Nonpositives(2), + Disjunct(Y[3])) + @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()) + for i in 1:5 + method.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) + +end + +function test_reformulate_disjunct() + model = GDPModel() + @variable(model, 1 <= x[1:2] <= 5) + @variable(model, Y[1:2], Logical) + @constraint(model, greaterthan, x[1] >= 1, Disjunct(Y[1])) + @constraint(model, interval, x[1] == 2.5, Disjunct(Y[2])) + + method = DP.MBM(HiGHS.Optimizer) + disj = constraint_object(disjunction(model, [Y[1], Y[2]])) + reformulated_disjunct = reformulate_disjunction(model, disj, method) + + @test length(reformulated_disjunct) == 3 + + @test reformulated_disjunct[1].set == MOI.GreaterThan(1.0) + @test reformulated_disjunct[2].set == MOI.GreaterThan(2.5) + @test reformulated_disjunct[3].set == MOI.LessThan(2.5) + + # Test that the expressions have the right structure + # Check coefficients and variables in the affine expressions + func_1 = reformulated_disjunct[1].func + func_2 = reformulated_disjunct[2].func + 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_2, x[1]) == 1.0 + @test JuMP.coefficient(func_2, binary_variable(Y[1])) == 2.5 + + @test JuMP.coefficient(func_3, x[1]) == 1.0 + @test JuMP.coefficient(func_3, binary_variable(Y[1])) == -2.5 +end + +function test_reformulate_disjunction() + model = GDPModel() + @variable(model, x) + @variable(model, Y[1:2], Logical) + @constraint(model, lessthan, x <= 2, Disjunct(Y[1])) + @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 + + @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_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 +end + +@testset "MBM" begin + test_mbm() + test__replace_variables_in_constraint() + test__constraint_to_objective() + test_mini_model() + test_maximize_M() + test_reformulate_disjunct_constraint() + test_reformulate_disjunct() + test_reformulate_disjunction() +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index d985839a..c653a3fe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,6 @@ import DisjunctiveProgramming as DP using DisjunctiveProgramming using Test - include("utilities.jl") # RUN ALL THE TESTS @@ -15,6 +14,7 @@ include("constraints/selector.jl") include("constraints/proposition.jl") include("constraints/disjunct.jl") include("constraints/indicator.jl") +include("constraints/mbm.jl") include("constraints/bigm.jl") include("constraints/hull.jl") include("constraints/fallback.jl") diff --git a/test/solve.jl b/test/solve.jl index a24f63d8..9f79f8e1 100644 --- a/test/solve.jl +++ b/test/solve.jl @@ -27,7 +27,6 @@ function test_linear_gdp_example(m, use_complements = false) @test value(Y[2]) @test !value(W[1]) @test !value(W[2]) - @test optimize!(m, gdp_method = Hull()) isa Nothing @test termination_status(m) == MOI.OPTIMAL @test objective_value(m) ≈ 11 @@ -46,6 +45,16 @@ function test_linear_gdp_example(m, use_complements = false) @test value(variable_by_name(m, "x[2]_Y[1]")) ≈ 0 @test value(variable_by_name(m, "x[2]_Y[2]")) ≈ 2 end + + @test optimize!(m, gdp_method = MBM(HiGHS.Optimizer)) isa Nothing + @test termination_status(m) == MOI.OPTIMAL + @test objective_value(m) ≈ 11 + @test value.(x) ≈ [9,2] + @test !value(Y[1]) + @test value(Y[2]) + @test !value(W[1]) + @test !value(W[2]) + end function test_generic_model(m)