Merge pull request #36 from hdavid16:split_interval_constraints

allow interval constraints (lb <= expr <= ub)

Project.toml

@@ -1,7 +1,7 @@
 name = "DisjunctiveProgramming"
 uuid = "0d27d021-0159-4c7d-b4a7-9ccb5d9366cf"
 authors = ["hdavid16 <[email protected]>"]
-version = "0.1.5"
+version = "0.1.6"
 
 [deps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"

README.md

@@ -10,7 +10,7 @@ Pkg.add("DisjunctiveProgramming")
 
 ## Disjunctions
 
-After defining a JuMP model, disjunctions can be added to the model by specifying which of the original JuMP model constraints should be assigned to each disjunction. The constraints that are assigned to the disjunctions will no longer be general model constraints, but will belong to the disjunction that they are assigned to. These constraints must be either `GreaterThan`, `LessThan`, or `EqualTo` constraints. Constraints that are of `Interval` type are currently not supported. It is assumed that the disjuncts belonging to a disjunction are proper disjunctions (mutually exclussive) and only one of them will be selected.
+After defining a JuMP model, disjunctions can be added to the model by specifying which of the original JuMP model constraints should be assigned to each disjunction. The constraints that are assigned to the disjunctions will no longer be general model constraints, but will belong to the disjunction that they are assigned to. These constraints must be either `GreaterThan`, `LessThan`, `EqualTo`, or `Interval` constraints. 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`).
 
 When a disjunction is defined using the `@disjunction` macro, the disjunctions are reformulated to algebraic constraints via either,
 - The Big-M method (when `reformulation = :BMR` in the `@disjunction` macro) 
@@ -48,34 +48,29 @@ using JuMP
 using DisjunctiveProgramming
 
 m = Model()
-@variable(m, 0<=x[1:2]<=10)
+@variable(m, -1<=x<=10)
 
-@constraint(m, con1[i=1:2], x[i] <= [3,4][i])
-@constraint(m, con2[i=1:2], zeros(2)[i] <= x[i])
-@constraint(m, con3[i=1:2], [5,4][i] <= x[i])
-@constraint(m, con4[i=1:2], x[i] <= [9,6][i])
+@constraint(m, con1, 0 <= x <= 3)
+@constraint(m, con2, 5 <= x <= 9)
 
-@disjunction(m,(con1,con2),(con3,con4), reformulation=:BMR, name = :y)
-@proposition(m, y[1] ∨ y[2])
+@disjunction(m,con1,con2,reformulation=:BMR,name=:y)
+@proposition(m, y[1] ∨ y[2]) #this is a redundant proposition
 
 print(m)
 
+┌ Warning: con1 : x in [0.0, 3.0] uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound.
+┌ Warning: con2 : x in [5.0, 9.0] uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound.
+
 Feasibility
 Subject to
- y[1] + y[2] == 1.0 #XOR constraint
- y[1] + y[2] >= 1.0 #reformulated logical proposition (redundant here)
- con1[1] : x[1] + 7 y[1] <= 10.0
- con1[2] : x[2] + 6 y[1] <= 10.0
- con2[1] : -x[1] <= 0.0
- con2[2] : -x[2] <= 0.0
- con3[1] : -x[1] + 5 y[2] <= 0.0
- con3[2] : -x[2] + 4 y[2] <= 0.0
- con4[1] : x[1] + y[2] <= 10.0
- con4[2] : x[2] + 4 y[2] <= 10.0
- x[1] >= 0.0
- x[2] >= 0.0
- x[1] <= 10.0
- x[2] <= 10.0
+ XOR(y) : y[1] + y[2] == 1.0
+ y[1] ∨ y[2] : y[1] + y[2] >= 1.0
+ con1[lb] : -x + y[1] <= 1.0
+ con1[ub] : x + 7 y[1] <= 10.0
+ con2[lb] : -x + 6 y[2] <= 1.0
+ con2[ub] : x + y[2] <= 10.0
+ x >= -1.0
+ x <= 10.0
  y[1] binary
  y[2] binary
 ```

examples/ex1.jl

@@ -4,28 +4,10 @@ using DisjunctiveProgramming
 m = Model()
 @variable(m, -1<=x<=10)
 
-@constraint(m, con1, x<=3)
-@constraint(m, con2, 0<=x)
-@constraint(m, con3, x<=9)
-@constraint(m, con4, 5<=x)
+@constraint(m, con1, 0 <= x <= 3)
+@constraint(m, con2, 5 <= x <= 9)
 
-@disjunction(m,(con1,con2),con3,con4,reformulation=:CHR,name=:y)
-# @proposition(m, y[1] ∨ y[2])
-# @proposition(m, 
-#     (
-#         ((y[1] ∨ y[2]) ∧ (¬y[1] ∨ ¬y[2]))
-#         ∨
-#         y[3]
-#     )
-#     ∧
-#     (
-#         ¬((y[1] ∨ y[2]) ∧ (¬y[1] ∨ ¬y[2]))
-#         ∨
-#         ¬y[3]
-#     )
-#     ∧
-#     (¬y[1] ∨ ¬y[2] ∨ ¬y[3])
-# )
-# @proposition(m, y[1] ⊻ y[2] ⊻ y[3])
+@disjunction(m,con1,con2,reformulation=:CHR,name=:y)
+@proposition(m, y[1] ∨ y[2]) #this is a redundant proposition
 
 print(m)

examples/ex2.jl

@@ -5,11 +5,9 @@ using DisjunctiveProgramming
 m = Model()
 @variable(m, 0<=x[1:2]<=10)
 
-@constraint(m, con1[i=1:2], x[i]<=[3,4][i])
-@constraint(m, con2, 0<=x[1])
-@constraint(m, con3, 0<=x[2])
-@constraint(m, con4[i=1:2], [5,4][i]<=x[i])
-@constraint(m, con5, x[1] <= 9)
-@constraint(m, con6, x[2] <= 6)
-
-@disjunction(m,(con1,con2,con3),(con4,con5,con6),reformulation=:BMR,name=:y)
+@constraint(m, con1[i=1:2], 0 <= x[i]<=[3,4][i])
+@constraint(m, con2[i=1:2], [5,4][i] <= x[i] <= [9,6][i])
+
+@disjunction(m,con1,con2,reformulation=:BMR,name=:y)
+
+print(m)

examples/ex3.jl

@@ -2,20 +2,11 @@ using JuMP
 using DisjunctiveProgramming
 
 m = Model()
-@variable(m, -1<=x<=10)
+@variable(m, -10 ≤ x[1:2] ≤ 10)
 
-@constraint(m, con1, x<=3)
-@constraint(m, con2, 0<=x)
-@constraint(m, con3, x<=9)
-@constraint(m, con4, 5<=x)
+nl_con1 = @NLconstraint(m, exp(x[1]) >= 1)
+nl_con2 = @NLconstraint(m, exp(x[2]) <= 2)
 
-M = 10
-@disjunction(
-    m,
-    (con1,con2),
-    con3,
-    con4,
-    reformulation=:BMR,
-    M=M,
-    name=:y
-)
+@disjunction(m, nl_con1, nl_con2, reformulation=:CHR, name=:z)
+
+print(m)

src/DisjunctiveProgramming.jl

@@ -5,6 +5,7 @@ using JuMP, IntervalArithmetic, Symbolics, Suppressor
 export add_disjunction, add_proposition
 export @disjunction, @proposition
 
+include("constraint.jl")
 include("logic.jl")
 include("utils.jl")
 include("big_M.jl")

src/big_M.jl

@@ -1,9 +1,7 @@
 function BMR!(m, constr, bin_var, i, k, M)
     if ismissing(k)
-        @assert is_valid(m,constr) "$constr is not a valid constraint in the model."
         ref = constr
     else
-        @assert is_valid(m,constr[k...]) "$constr is not a valid constraint in the model."
         ref = constr[k...]
     end
     if ismissing(M)

src/constraint.jl

@@ -0,0 +1,131 @@
+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
+JuMP.name(con_ref::NonlinearConstraintRef) = ""
+
+function check_disjunction(m, disj)
+    disj_new = [] #create a new array where the disjunction will be copied to so that we can split constraints that use an Interval set
+    for constr in disj
+        if constr isa Tuple #NOTE: Make it so that it must be bundled in a Tuple (not Array), to avoid confusing it with a Variable Array
+            constr_list = []
+            for constr_j in constr
+                push!(constr_list, check_constraint!(m, constr_j))
+            end
+            push!(disj_new, Tuple(constr_list))
+        elseif constr isa Union{ConstraintRef, Array, Containers.DenseAxisArray, Containers.SparseAxisArray}
+            push!(disj_new, check_constraint!(m, constr))
+        end
+    end
+
+    return disj_new
+end
+
+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)
+    if constr isa ConstraintRef
+        new_constr = split_interval_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))
+            new_constr = constr
+        else
+            split_flag = true
+            if constr isa Union{Array, Containers.DenseAxisArray}
+                idxs = Iterators.product(axes(constr)...)
+            elseif constr isa Containers.SparseAxisArray
+                idxs = keys(constr.data)
+            end
+            constr_dict = Dict(union(
+                [
+                    split_interval_constraint(m, constr[idx...]) |>
+                        i -> isnothing(i) ? 
+                            (idx...,"") => constr[idx...] : 
+                            [(idx...,"lb") => i[1], (idx...,"ub") => i[2]]
+                    for idx in idxs
+                ]...
+            ))
+            new_constr = Containers.SparseAxisArray(constr_dict)
+            m[constr_name] = new_constr
+        end
+    end
+
+    split_flag && @warn "$constr uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound."
+    delete_original_constraint!(m, constr)
+
+    return new_constr
+end
+
+function gen_constraint_name(constr)
+    constr_name = name.(constr)
+    if any(isempty.(constr_name))
+        constr_name = gensym("constraint")
+    elseif !isa(constr_name, String)
+        c_names = union(first.(split.(constr_name,"[")))
+        if length(c_names) == 1
+            constr_name = c_names[1]
+        else
+            constr_name = gensym("constraint")
+        end
+    end
+
+    return Symbol("$(constr_name)_split")
+end
+
+function split_interval_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
+            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]
+        end
+    elseif constr isa ConstraintRef
+        constr_obj = constraint_object(constr)
+        if constr_obj.set isa MOI.Interval
+            lb = constr_obj.set.lower
+            ub = constr_obj.set.upper
+            ex = constr_obj.func
+            return [
+                @constraint(m, lb <= ex, base_name = "$(constr_name)[lb]"),
+                @constraint(m, ex <= ub, base_name = "$(constr_name)[ub]")
+            ]
+        end
+    end
+    return nothing
+end
+
+function delete_original_constraint!(m, constr)
+    if constr isa ConstraintRef
+        if !isa(constr, NonlinearConstraintRef)
+            delete(m, constr)
+            # unregister(m, constr)
+        end
+    elseif constr isa Union{Array, Containers.DenseAxisArray, Containers.SparseAxisArray}
+        if !isa(first(constr), NonlinearConstraintRef)
+            delete.(m, constr)
+            # unregister(m, constr)
+        end
+    end
+end

src/convex_hull.jl

@@ -1,6 +1,5 @@
 function CHR!(m, constr, bin_var, i, k, eps)
     ref = ismissing(k) ? constr : constr[k...] #get constraint
-    @assert is_valid(m,ref) "$constr is not a valid constraint in the model."
     #create convex hull constraint
     if ref isa NonlinearConstraintRef || constraint_object(ref).func isa QuadExpr
         nonlinear_perspective_function(ref, bin_var, i, eps)
@@ -66,8 +65,10 @@ function disaggregate_variables(m, disj, bin_var)
                 m[var_name_i] = Containers.SparseAxisArray(var_i_dict)
             end
             #apply bounding constraints on disaggregated variable
-            m[Symbol(var_name_i,"_lb")] = @constraint(m, LB * m[bin_var][i] .- m[var_name_i] .<= 0)
-            m[Symbol(var_name_i,"_ub")] = @constraint(m, m[var_name_i] .- UB * m[bin_var][i] .<= 0)
+            var_i_lb = "$(var_name_i)_lb" 
+            var_i_ub = "$(var_name_i)_ub" 
+            m[Symbol(var_i_lb)] = @constraint(m, LB * m[bin_var][i] .- m[var_name_i] .<= 0, base_name = var_i_lb)
+            m[Symbol(var_i_ub)] = @constraint(m, m[var_name_i] .- UB * m[bin_var][i] .<= 0, base_name = var_i_ub)
         end
     end
 end
@@ -80,7 +81,8 @@ function sum_disaggregated_variables(m, disj, bin_var)
             var_name_i = Symbol("$(var_name)_$bin_var$i")
             push!(dis_vars, m[var_name_i])
         end
-        m[Symbol(var_name,"_aggregation")] = @constraint(m, var .== sum(dis_vars))
+        aggr_con = "$(var_name)_aggregation"
+        m[Symbol(aggr_con)] = @constraint(m, var .== sum(dis_vars), base_name = aggr_con)
     end
 end
 
@@ -97,8 +99,6 @@ end
 
 function linear_perspective_function(ref, bin_var, i)
     #check constraint type
-    ref_obj = constraint_object(ref)
-    @assert ref_obj.set isa MOI.LessThan || ref_obj.set isa MOI.GreaterThan || ref_obj.set isa MOI.EqualTo "$ref must be one the following: GreaterThan, LessThan, or EqualTo."
     bin_var_ref = ref.model[bin_var][i]
     #replace each variable with its disaggregated version
     for var_ref in ref.model[:gdp_variable_refs]

src/logic.jl

@@ -18,9 +18,9 @@ function to_cnf!(m::Model, expr::Expr)
     unique!(lhs)
     #generate JuMP constraints for the logical proposition
     if length(lhs) == 1
-        m[expr_name] = @constraint(m, lhs[1] >= 1)
+        m[expr_name] = @constraint(m, lhs[1] >= 1, base_name = string(expr_name))
     else
-        m[expr_name] = @constraint(m, [i = eachindex(lhs)], lhs[i] >= 1)
+        m[expr_name] = @constraint(m, [i = eachindex(lhs)], lhs[i] >= 1, base_name = string(expr_name))
     end
 end
 

src/macro.jl

@@ -33,7 +33,8 @@ function add_disjunction(m::Model,disj...;reformulation::Symbol,M=missing,ϵ=1e-
     #create variable if it doesn't exist
     m[disj_name] = @variable(m, [eachindex(disj)], Bin, base_name = string(disj_name))
     #add xor constraint on binary variable
-    m[Symbol("XOR(",disj_name,")")] = @constraint(m, sum(m[disj_name][i] for i in eachindex(disj)) == 1)
+    xor_con = "XOR($disj_name)"
+    m[Symbol(xor_con)] = @constraint(m, sum(m[disj_name][i] for i in eachindex(disj)) == 1, base_name = xor_con)
     #apply reformulation
     param = reformulation == :BMR ? M : ϵ
     reformulate_disjunction(m, disj, disj_name, reformulation, param)