Merge branch 'hdavid16:master' into master

Author: shivankj11

README.md

@@ -19,7 +19,9 @@ After defining a JuMP model, disjunctions can be added to the model by using the
 
 NOTES: 
 - Vectorized constraints (using `.` notation) are not currently supported. The current workarround is to first creating the constraint outside of the `@disjunction` macro and then passing the reference to the constraint to the `@disjunction` macro.
-- Any constraints that are of `Interval` type are split into two constraints (one for each bound). It is assumed that the disjuncts belonging to a disjunction are proper disjunctions (mutually exclussive) and only one of them will be selected (`XOR`).
+- Any constraints that are of `EqualTo` type are split into two constraints (e.g., `f(x) == 0` -> `0 <= f(x) <= 0`). This is necessary only for the Big-M reformulation of equality constraints, but is currently applied regardless of the reformulation technique.
+- Any constraints that are of `Interval` type are split into two constraints (one for each bound).
+- It is assumed that the disjuncts belonging to a disjunction are proper disjunctions (mutually exclussive) and only one of them will be selected (`XOR`).
 
 The valid key-word arguments for the `@disjunction` macro are:
 - `reformulation::Symbol`: `:big_m` for [Big-M Reformulation](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities#Big-M_Reformulation), `:convex_hull` for [Convex-Hull Reformulation](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities#Convex-Hull_Reformulation)

examples/ex3.jl

@@ -5,14 +5,10 @@ m = Model()
 @variable(m, -5 ≤ x ≤ 10)
 @disjunction(
     m,
-    # begin
-    #     exp(x) ≤ 2
-    #     -3 ≤ x
-    # end,
-    (
-        exp(x) ≤ 2,
+    begin
+        exp(x) ≤ 2
         -3 ≤ x
-    ),
+    end,
     begin
         3 ≤ exp(x)
         5 ≤ x
@@ -46,4 +42,5 @@ print(m)
 #  Perspective Functions:
 #  (-1.0e-6 + -1.9999989999999999 * z[1]) + (1.0e-6 + 0.999999 * z[1]) * exp(x_z1 / (1.0e-6 + 0.999999 * z[1])) <= 0
 #  (1.0000000000000002e-6 + 2.999999 * z[2]) + (-1.0e-6 + -0.999999 * z[2]) * exp(x_z2 / (1.0e-6 + 0.999999 * z[2])) <= 0
+#  -1.0 * x_z1 + -3.0 * z[1] <= 0
 #  -1.0 * x_z2 + 5.0 * z[2] <= 0

src/big_M.jl

@@ -1,4 +1,4 @@
-function big_m_reformulation!!(constr, bin_var, i, k, M)
+function big_m_reformulation!(constr, bin_var, i, k, M)
     if ismissing(k)
         ref = constr
     else
@@ -7,6 +7,8 @@ function big_m_reformulation!!(constr, bin_var, i, k, M)
     if ismissing(M)
         M = apply_interval_arithmetic(ref)
         # @warn "No M value passed for $ref. M = $M was inferred from the variable bounds."
+    elseif !ismissing(k)
+        M = M[k]
     end
     if ref isa NonlinearConstraintRef
         nonlinear_bigM(ref, bin_var, M, i)

src/constraint.jl

@@ -1,21 +1,24 @@
 is_interval_constraint(con_ref::ConstraintRef{<:AbstractModel}) = constraint_object(con_ref).set isa MOI.Interval
 is_interval_constraint(con_ref::NonlinearConstraintRef) = count(i -> i == :(<=), Meta.parse(string(con_ref)).args) == 2
+is_equality_constraint(con_ref::ConstraintRef{<:AbstractModel}) = constraint_object(con_ref).set isa MOI.EqualTo
+is_equality_constraint(con_ref::NonlinearConstraintRef) = Meta.parse(string(con_ref)).args[1] == :(==)
 JuMP.name(con_ref::NonlinearConstraintRef) = ""
 
 function check_constraint!(m, constr)
     @assert all(is_valid.(m, constr)) "$constr is not a valid constraint."
     split_flag = false
     constr_name = gen_constraint_name(constr)
+    constr_str = split(string(constr),"}")[end]
     if constr isa ConstraintRef
-        new_constr = split_interval_constraint(m, constr)
+        new_constr = split_constraint(m, constr)
         if isnothing(new_constr)
             new_constr = constr
         else
             split_flag = true
             m[constr_name] = new_constr
         end
     elseif constr isa Union{Array, Containers.DenseAxisArray, Containers.SparseAxisArray}
-        if !any(is_interval_constraint.(constr))
+        if !any(is_interval_constraint.(constr)) && !any(is_equality_constraint.(constr))
             new_constr = constr
         else
             split_flag = true
@@ -26,7 +29,7 @@ function check_constraint!(m, constr)
             end
             constr_dict = Dict(union(
                 [
-                    split_interval_constraint(m, constr[idx...]) |>
+                    split_constraint(m, constr[idx...]) |>
                         i -> isnothing(i) ? 
                             (idx...,"") => constr[idx...] : 
                             [(idx...,"lb") => i[1], (idx...,"ub") => i[2]]
@@ -39,7 +42,7 @@ function check_constraint!(m, constr)
     end
 
     if split_flag
-        @warn "$(split(string(constr),"}")[end]) uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound."
+        @warn "$constr_str uses the `MOI.Interval` or `MOI.EqualTo` set. Each instance of the interval set has been split into two constraints, one for each bound."
         delete_original_constraint!(m, constr)
     end
 
@@ -62,47 +65,67 @@ function gen_constraint_name(constr)
     return Symbol("$(constr_name)_split")
 end
 
-function split_interval_constraint(m, constr, constr_name = name(constr))
+function split_constraint(m, constr, constr_name = name(constr))
     if isempty(constr_name)
         constr_name = "[$constr]"
     end
     if constr isa NonlinearConstraintRef
         constr_expr = Meta.parse(string(constr))
-        if count(x -> x == :(<=), constr_expr.args) == 2
+        if constr_expr.args[1] == :(==) #replace == for lb <= expr <= ub and split
+            lb, ub = 0, 0#rhs is always 0, but could get obtained from: constr_expr.args[3]
+            constr_expr_func = copy(constr_expr.args[2])
+            new_constraints = split_nonlinear_constraint(m, constr, constr_expr_func, lb, ub)
+            return new_constraints
+        elseif count(x -> x == :(<=), constr_expr.args) == 2 #split lb <= expr <= ub
             lb = constr_expr.args[1]
             ub = constr_expr.args[5]
             constr_expr_func = copy(constr_expr.args[3]) #get func part of constraint
-            replace_JuMPvars!(constr_expr_func, m) #replace Expr with JuMP vars
-            #replace original constraint with lb <= func
-            lb_constr = JuMP._NonlinearConstraint(
-                JuMP._NonlinearExprData(m, constr_expr_func), 
-                lb, 
-                Inf
-            )
-            m.nlp_data.nlconstr[constr.index.value] = lb_constr
-            #create new constraint for func <= ub
-            constr_expr_ub = Expr(:call, :(<=), constr_expr_func, ub)
-            ub_constr = add_nonlinear_constraint(m, constr_expr_ub)
-            #return split constraint
-            return [constr, ub_constr]
+            new_constraints = split_nonlinear_constraint(m, constr, constr_expr_func, lb, ub)
+            return new_constraints
         end
     elseif constr isa ConstraintRef
         constr_obj = constraint_object(constr)
-        if constr_obj.set isa MOI.Interval
+        if constr_obj.set isa MOI.EqualTo
+            lb = constr_obj.set.value
+            ub = lb
+            new_constraints = split_linear_constraint(m, constr_obj, constr_name, lb, ub)
+            return new_constraints
+        elseif constr_obj.set isa MOI.Interval
             lb = constr_obj.set.lower
             ub = constr_obj.set.upper
-            ex = constr_obj.func
-            lb_name = name_split_constraint(constr_name, :lb)
-            ub_name = name_split_constraint(constr_name, :ub)
-            return [
-                @constraint(m, lb <= ex, base_name = lb_name),
-                @constraint(m, ex <= ub, base_name = ub_name)
-            ]
+            new_constraints = split_linear_constraint(m, constr_obj, constr_name, lb, ub)
+            return new_constraints
         end
     end
     return nothing
 end
 
+function split_linear_constraint(m, constr_obj, constr_name, lb, ub)
+    ex = constr_obj.func
+    lb_name = name_split_constraint(constr_name, :lb)
+    ub_name = name_split_constraint(constr_name, :ub)
+    return [
+        @constraint(m, lb <= ex, base_name = lb_name),
+        @constraint(m, ex <= ub, base_name = ub_name)
+    ]
+end
+
+function split_nonlinear_constraint(m, constr, constr_expr_func, lb, ub)
+    replace_JuMPvars!(constr_expr_func, m) #replace Expr with JuMP vars
+    #replace original constraint with lb <= func
+    lb_constr = JuMP._NonlinearConstraint(
+        JuMP._NonlinearExprData(m, constr_expr_func), 
+        lb, 
+        Inf
+    )
+    m.nlp_data.nlconstr[constr.index.value] = lb_constr
+    #create new constraint for func <= ub
+    constr_expr_ub = Expr(:call, :(<=), constr_expr_func, ub)
+    ub_constr = add_nonlinear_constraint(m, constr_expr_ub)
+
+    return [constr, ub_constr]
+end
+
 function delete_original_constraint!(m, constr)
     if constr isa ConstraintRef
         if !isa(constr, NonlinearConstraintRef)

src/convex_hull.jl

@@ -1,4 +1,4 @@
-function convex_hull_reformulation!!(constr, bin_var, i, k, eps)
+function convex_hull_reformulation!(constr, bin_var, i, k, eps)
     ref = ismissing(k) ? constr : constr[k...] #get constraint
     #create convex hull constraint
     if ref isa NonlinearConstraintRef || constraint_object(ref).func isa QuadExpr
@@ -13,11 +13,17 @@ function disaggregate_variables(m, disj, bin_var)
     var_refs = m[:gdp_variable_refs]
     @assert all((has_upper_bound.(var_refs) .&& has_lower_bound.(var_refs)) .|| is_binary.(var_refs)) "All variables must be bounded to perform the Convex-Hull reformulation."
     #reformulate variables
+    obj_dict = object_dictionary(m)
+    bounds_dict = :variable_bounds_dict in keys(obj_dict) ? obj_dict[:variable_bounds_dict] : Dict() #NOTE: should pass as an keyword argument
     for var_name in m[:gdp_variable_names]
         var = m[var_name]
         #define UB and LB
         if var isa VariableRef
-            LB, UB = get_bounds(var)
+            if string(var) in keys(bounds_dict)
+                LB, UB = bounds_dict[string(var)]
+            else
+                LB, UB = get_bounds(var)
+            end
         elseif var isa Array{VariableRef} || var isa Containers.DenseAxisArray || var isa Containers.SparseAxisArray
             #initialize UB and LB with same container type as variable
             if var isa Array{VariableRef} || var isa Containers.DenseAxisArray
@@ -34,7 +40,11 @@ function disaggregate_variables(m, disj, bin_var)
             end
             #populate UB and LB
             for idx in eachindex(var)
-                LB[idx], UB[idx] = get_bounds(var[idx])
+                if string(var[idx]) in keys(bounds_dict)
+                    LB[idx], UB[idx] = bounds_dict[string(var[idx])]
+                else
+                    LB[idx], UB[idx] = get_bounds(var[idx])
+                end
             end
         end
         #disaggregate variable and add bounding constraints
@@ -99,8 +109,8 @@ end
 function add_disaggregated_variable(m, LB, UB, var, base_name)
     @variable(
         m, 
-        lower_bound = LB, 
-        upper_bound = UB, 
+        lower_bound = min(LB,0), 
+        upper_bound = max(UB,0), 
         binary = is_binary(var), 
         integer = is_integer(var),
         base_name = base_name

src/reformulate.jl

@@ -22,11 +22,10 @@ function reformulate_disjunction(m::Model, disj...; bin_var, reformulation, para
         Symbol(bin_var,"[$i]") => disj[i] for i in eachindex(disj)
     )
     new_constraints[Symbol(bin_var,"_XOR")] = constraint_by_name(m, "XOR(disj_$bin_var)")
-    for var in vars
-        agg_con_name = "$(var)_$(bin_var)_aggregation"
-        agg_con = constraint_by_name(m, agg_con_name)
-        if !isnothing(agg_con)
-            new_constraints[Symbol(agg_con_name)] = agg_con
+    if reformulation == :convex_hull
+        for var in m[:gdp_variable_refs]
+            agg_con_name = "$(var)_$(bin_var)_aggregation"
+            new_constraints[Symbol(agg_con_name)] = constraint_by_name(m, agg_con_name)
         end
     end
     return new_constraints
@@ -53,6 +52,8 @@ function check_disjunction!(m, disj)
             push!(disj_new, Tuple(constr_list))
         elseif constr isa Union{ConstraintRef, Array, Containers.DenseAxisArray, Containers.SparseAxisArray}
             push!(disj_new, check_constraint!(m, constr))
+        elseif isnothing(constr)
+            push!(disj_new, constr)
         end
     end
 
@@ -84,9 +85,9 @@ end
 
 function call_reformulation(reformulation, constr, bin_var, i, k, param)
     if reformulation == :big_m
-        big_m_reformulation!!(constr, bin_var, i, k, param)
+        big_m_reformulation!(constr, bin_var, i, k, param)
     elseif reformulation == :convex_hull
-        convex_hull_reformulation!!(constr, bin_var, i, k, param)
+        convex_hull_reformulation!(constr, bin_var, i, k, param)
     end
 end
 

src/utils.jl

@@ -10,11 +10,18 @@ function apply_interval_arithmetic(ref)
     ref_func_expr, ref_type, ref_rhs = parse_constraint(ref)
     #create a map of variables to their bounds
     interval_map = Dict()
-    vars = ref.model[:gdp_variable_refs]
+    vars = all_variables(ref.model)#ref.model[:gdp_variable_refs]
+    obj_dict = object_dictionary(ref.model)
+    bounds_dict = :variable_bounds_dict in keys(obj_dict) ? obj_dict[:variable_bounds_dict] : Dict() #NOTE: should pass as an keyword argument
     for var in vars
-        UB = has_upper_bound(var) ? upper_bound(var) : (is_binary(var) ? 1 : Inf)
-        LB = has_lower_bound(var) ? lower_bound(var) : (is_binary(var) ? 0 : -Inf)
-        interval_map[string(var)] = LB..UB
+        if string(var) in keys(bounds_dict)
+            bnds = bounds_dict[string(var)]
+            interval_map[string(var)] = bnds[1]..bnds[2]
+        else
+            UB = has_upper_bound(var) ? upper_bound(var) : (is_binary(var) ? 1 : Inf)
+            LB = has_lower_bound(var) ? lower_bound(var) : (is_binary(var) ? 0 : -Inf)
+            interval_map[string(var)] = LB..UB
+        end
     end
     ref_func_expr = replace_intevals!(ref_func_expr, interval_map)
     #get bounds on the entire expression
@@ -34,12 +41,13 @@ end
 function parse_constraint(ref)
     if ref isa NonlinearConstraintRef
         ref_str = string(ref)
-        ref_func = replace(split(ref_str, r"[=<>]")[1], " " => "")
+        ref_func = replace(split(ref_str, r"[=<>]")[1], " " => "") #remove operator and spaces
         ref_type = occursin(">", ref_str) ? :lower : :upper
         ref_rhs = 0 #Could be calculated with: parse(Float64,split(ref_str, " ")[end]). NOTE: @NLconstraint will always have a 0 RHS.
     elseif ref isa ConstraintRef
         ref_obj = constraint_object(ref)
-        ref_func = string(ref_obj.func)
+        ref_str = string(ref_obj.func)
+        ref_func = replace(ref_str, " " => "") #remove spaces
         ref_type = fieldnames(typeof(ref_obj.set))[1]
         ref_rhs = normalized_rhs(ref)
     end