Support Multiple Big-M reformulation.

README.md

@@ -172,6 +172,10 @@ 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.
+
 ## 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).

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,6 +21,7 @@ include("macros.jl")
 include("reformulate.jl")
 include("bigm.jl")
 include("hull.jl")
+include("mbm.jl")
 include("indicator.jl")
 include("print.jl")
 

src/datatypes.jl

@@ -367,6 +367,25 @@ struct BigM{T} <: AbstractReformulationMethod
     end
 end
 
+"""
+    MBM{O} <: 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).
+"""
+struct MBM{O} <: AbstractReformulationMethod
+    optimizer::O
+
+    # Constructor with optimizer (required)
+    function MBM(optimizer::O) where {O}
+        new{O}(optimizer)
+    end
+end
+
+
+
 """
     Hull{T} <: AbstractReformulationMethod
 

src/mbm.jl

@@ -0,0 +1,332 @@
+################################################################################
+#               CONSTRAINT, DISJUNCTION, DISJUNCT REFORMULATION
+################################################################################
+#Reformulates the disjunction using multiple big-M values
+function reformulate_disjunction(model::JuMP.AbstractModel, 
+    disj::Disjunction,
+    method::MBM
+    )
+    ref_cons = Vector{JuMP.AbstractConstraint}() 
+    for d in disj.indicators
+        constraints = copy(_indicator_to_constraints(model)[d])  # Make a copy to avoid modifying during iteration
+        for cref in constraints
+            con = JuMP.constraint_object(cref)
+            if con isa JuMP.ScalarConstraint{<:Any, <:MOI.Interval}
+                    # Create lower and upper bound constraints
+                    lower_con = JuMP.build_constraint(error, con.func, MOI.GreaterThan(con.set.lower))
+                    upper_con = JuMP.build_constraint(error, con.func, MOI.LessThan(con.set.upper))
+                    # Create new disjunct constraints
+                    JuMP.add_constraint(model, _DisjunctConstraint(lower_con, d))
+                    JuMP.add_constraint(model, _DisjunctConstraint(upper_con, d))
+                    JuMP.delete(model, cref)
+                end 
+            end
+        end
+    for d in disj.indicators
+        _reformulate_disjunct(model, ref_cons, d, LogicalVariableRef[x for x in disj.indicators if x != d], method)
+    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,
+    conlvref::Vector{LogicalVariableRef},
+    method::AbstractReformulationMethod
+    )
+
+    bconref = Dict(d => binary_variable(d) for d in conlvref)
+    !haskey(_indicator_to_constraints(model), lvref) && return
+    constraints = _indicator_to_constraints(model)[lvref]
+    M = Dict{LogicalVariableRef,Float64}()
+    sizehint!(M, length(conlvref))
+    for d in conlvref 
+        M[d] = maximum(_maximize_M(model, JuMP.constraint_object(cref), 
+        Vector{DisjunctConstraintRef}(_indicator_to_constraints(model)[d]), method) for cref in constraints)
+    end
+
+    isempty(constraints) || sizehint!(ref_cons, length(constraints))
+    for cref in constraints
+        con = JuMP.constraint_object(cref)
+        append!(ref_cons, reformulate_disjunct_constraint(model, con, Dict{LogicalVariableRef,VariableRef}(bconref), M, method))
+    end
+    return ref_cons
+end
+
+function reformulate_disjunct_constraint(
+    model::JuMP.AbstractModel,
+    con::JuMP.VectorConstraint{T, S, R},
+    bconref:: Dict{LogicalVariableRef,JuMP.VariableRef},
+    M::Dict{LogicalVariableRef,Float64},
+    method::MBM
+) where {T, S <: _MOI.Nonpositives, R}
+    m_sum = sum(M[i] * bconref[i] for i in keys(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:: Dict{LogicalVariableRef,JuMP.VariableRef},
+    M::Dict{LogicalVariableRef,Float64},
+    method::MBM
+) where {T, S <: _MOI.Nonnegatives, R}
+    m_sum = sum(M[i] * bconref[i] for i in keys(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:: Dict{LogicalVariableRef,JuMP.VariableRef},
+    M::Dict{LogicalVariableRef,Float64},
+    method::MBM
+) where {T, S <: _MOI.Zeros, R}
+    m_sum = sum(M[i] * bconref[i] for i in keys(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:: Dict{LogicalVariableRef,JuMP.VariableRef},
+    M::Dict{LogicalVariableRef,Float64},
+    method::MBM
+) where {T, S <: _MOI.LessThan}
+    new_func = JuMP.@expression(model, con.func - sum(M[i] * bconref[i] for i in keys(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:: Dict{LogicalVariableRef,JuMP.VariableRef},
+    M::Dict{LogicalVariableRef,Float64},
+    method::MBM
+) where {T, S <: _MOI.GreaterThan}
+    new_func = JuMP.@expression(model, con.func + sum(M[i] * bconref[i] for i in keys(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:: Dict{LogicalVariableRef,JuMP.VariableRef},
+    M::Dict{LogicalVariableRef,Float64},
+    method::MBM
+) where {T, S <: _MOI.EqualTo}
+    upper_func = JuMP.@expression(model, con.func - sum(M[i] * bconref[i] for i in keys(M)))
+    lower_func = JuMP.@expression(model, con.func + sum(M[i] * bconref[i] for i in keys(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::Any,
+    bconref:: Dict{LogicalVariableRef,JuMP.VariableRef},
+    M::Dict{LogicalVariableRef,Float64},
+    method::MBM
+) 
+    error("Constraint type $(typeof(con)) is not supported by the Multiple Big-M reformulation method.")
+end
+
+################################################################################
+#                          MULTIPLE BIG-M REFORMULATION
+################################################################################
+#Dispatches over constaint 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}
+    return maximum(_maximize_M(model, JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(0.0)), 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}
+    return maximum(_maximize_M(model, JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(0.0)), 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}
+    return max(
+        maximum(_maximize_M(model, JuMP.ScalarConstraint(objective.func[i], MOI.GreaterThan(0.0)), constraints, method) for i in 1:objective.set.dimension),
+        maximum(_maximize_M(model, JuMP.ScalarConstraint(objective.func[i], MOI.LessThan(0.0)), 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}
+    return max(
+        _mini_model(model, JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(objective.set.value)), constraints, method),
+        _mini_model(model, JuMP.ScalarConstraint(objective.func, MOI.LessThan(objective.set.value)), 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}}
+    sub_model = JuMP.Model()
+    new_vars = Dict{JuMP.VariableRef, JuMP.VariableRef}()
+    for var in JuMP.all_variables(model)
+        new_vars[var] = JuMP.@variable(sub_model, base_name= "sub_model_$(JuMP.name(var))")
+        if JuMP.is_fixed(var)
+            JuMP.fix(new_vars[var], JuMP.fix_value(var); force=true)
+        end
+        if JuMP.is_integer(var)
+            JuMP.set_integer(new_vars[var])
+        end
+        if JuMP.has_upper_bound(var)
+            JuMP.set_upper_bound(new_vars[var], JuMP.upper_bound(var))
+        end
+        if JuMP.has_lower_bound(var)
+            JuMP.set_lower_bound(new_vars[var], JuMP.lower_bound(var))
+        end
+        if JuMP.has_start_value(var)
+            JuMP.set_start_value(new_vars[var], JuMP.start_value(var))
+        end
+    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 = 1e9
+    else
+        M = JuMP.objective_value(sub_model)
+    end
+    return M
+end
+#Catches any constraints that were not reformulated in _maximize_M
+#_mini_model requires objective to be >= or <= in order to run
+function _mini_model(
+    model::JuMP.AbstractModel, 
+    objective::JuMP.ScalarConstraint{T,S}, 
+    constraints::Vector{DisjunctConstraintRef}, 
+    method::MBM
+) where {T,S <: Union{_MOI.Nonpositives, _MOI.Nonnegatives, _MOI.Zeros, MOI.EqualTo, MOI.Interval}}
+    error("This type of constraints and objective constraint has not been implemented for MBM subproblems")
+end
+
+################################################################################
+#                          CONSTRAINT TO OBJECTIVE
+################################################################################
+function constraint_to_objective(sub_model::JuMP.AbstractModel,obj::JuMP.ScalarConstraint{<:JuMP.AbstractJuMPScalar, MOI.LessThan{Float64}}, new_vars::Dict{JuMP.VariableRef, JuMP.VariableRef})
+    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{Float64}}, new_vars::Dict{JuMP.VariableRef, JuMP.VariableRef})
+    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{JuMP.VariableRef, JuMP.VariableRef})
+    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.GenericVariableRef, var_map::Dict{JuMP.VariableRef, JuMP.VariableRef})
+    return var_map[fun]
+end
+
+function replace_variables_in_constraint(fun::JuMP.AffExpr, var_map::Dict{JuMP.VariableRef, JuMP.VariableRef})
+    new_aff = JuMP.zero(JuMP.AffExpr)
+    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::JuMP.QuadExpr, var_map::Dict{JuMP.VariableRef, JuMP.VariableRef})
+    new_quad = JuMP.zero(JuMP.QuadExpr)
+    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.VariableRef, JuMP.VariableRef})
+    return fun
+end
+
+function replace_variables_in_constraint(fun::JuMP.NonlinearExpr, var_map::Dict{JuMP.VariableRef, JuMP.VariableRef})
+    new_args = Any[replace_variables_in_constraint(arg, var_map) for arg in fun.args]
+    return JuMP.NonlinearExpr(fun.head, new_args)
+end
+
+function replace_variables_in_constraint(fun::Vector{T}, var_map::Dict{JuMP.VariableRef, JuMP.VariableRef}) where T
+    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