Merge branch 'hdavid16:master' into master

Author: shivankj11

Project.toml

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

README.md

@@ -18,18 +18,18 @@ After defining a JuMP model, disjunctions can be added to the model by using the
 5. `Tuple` of expressions accepted by options 1 and/or 3.
 
 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.
+- Vectorized constraints (using `.` notation) are not currently supported. The current workarround is to first create the constraint with the `@constraint` macro and then use the `add_disjunction!`, instead of the `@disjunction` macro. The `add_disjunction!` function receives the same arguments as the `@disjunction` macro, with the exception that instead of creating the constraints in the disjunctions, references to previously created constraints are used for the disjuncts.
 - 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)
+- `reformulation::Symbol`: `:big_m` for [Big-M Reformulation](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities#Big-M_Reformulation), `:hull` for [Hull Reformulation](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities#Convex-Hull_Reformulation)
 - `M`: Big-M value used when `reformulation = :big_m`.
-- `ϵ`: epsilon tolerance for the perspective function proposed by [Furman, et al. [2020]](https://link.springer.com/article/10.1007/s10589-020-00176-0). Only used when `reformulation = :convex_hull`.
+- `ϵ`: epsilon tolerance for the perspective function proposed by [Furman, et al. [2020]](https://link.springer.com/article/10.1007/s10589-020-00176-0). Only used when `reformulation = :hull`.
 - `name::Symbol`: Name for the disjunction (also name for indicator variable used on that disjunction). If not passed (`name = missing`), a symbolic name will be generated with the prefix `disj`. The mutual exclussion constraint on the binary indicator variables can be accessed with `model[Symbol("XOR(disj_$name)")]`.
 
-When a disjunction is defined using the `@disjunction` macro, the disjunctions are reformulated to algebraic constraints via either Big-M or Convex-Hull reformulations. For the Convex-Hull reformulation, disaggregated variables are generated by adding the suffix `_$name$i` to the original variables, where `i` is the index of the disjunct in that disjunction. Bounding constraints are applied to the disaggregated variables and can be accessed with `model[Symbol("$<original var>_$name$i_lb")]` and `model[Symbol("$<original var>_$name$i_ub")]` for the lower bound and upper bound constraints, respectively. The aggregation constraint can be accessed with `model[Symbol("$<original var>_aggregation")]`. For Big-M reformulations, the user may provide an `M` object that represents the BigM value(s). The `M` object can be a `Number` that is applied to all constraints in the disjuncts, or a `Vector`/`Tuple` of values that are used for each of the disjuncts. For Convex-Hull reformulations, the user may provide an `ϵ` value for the perspective function (default is `ϵ = 1e-6`). The `ϵ` object can be a `Number` that is applied to all perspective functions, or a `Vector`/`Tuple` of values that are used for each of the disjuncts.
+When a disjunction is defined using the `@disjunction` macro, the disjunctions are reformulated to algebraic constraints via either Big-M or Hull reformulations. For the Hull reformulation, disaggregated variables are generated by adding the suffix `_$name$i` to the original variables, where `i` is the index of the disjunct in that disjunction. Bounding constraints are applied to the disaggregated variables and can be accessed with `model[Symbol("$<original var>_$name$i_lb")]` and `model[Symbol("$<original var>_$name$i_ub")]` for the lower bound and upper bound constraints, respectively. The aggregation constraint can be accessed with `model[Symbol("$<original var>_aggregation")]`. For Big-M reformulations, the user may provide an `M` object that represents the BigM value(s). The `M` object can be a `Number` that is applied to all constraints in the disjuncts, or a `Vector`/`Tuple` of values that are used for each of the disjuncts. For Hull reformulations, the user may provide an `ϵ` value for the perspective function (default is `ϵ = 1e-6`). The `ϵ` object can be a `Number` that is applied to all perspective functions, or a `Vector`/`Tuple` of values that are used for each of the disjuncts.
 
 For empty disjuncts, use `nothing` for their positional argument (e.g., `@disjunction(m, x <= 1, nothing, reformulation = :big_m)`).
 
@@ -51,7 +51,7 @@ The example below is from the [Northwestern University Process Optimization Open
 
 To perform the Big-M reformulation, `:big_m` is passed to the `reformulation` keyword argument. If nothing is passed to the keyword argument `M`, tight Big-M values will be inferred from the variable bounds using IntervalArithmetic.jl. If `x` is not bounded, Big-M values must be provided for either the whole system (e.g., `M = 10`) or for each of the constraint arrays in the example (e.g., `M = (10,10)`).
 
-To perform the Convex-Hull reformulation, `reformulation = :convex_hull`. Variables must have bounds for the reformulation to work.
+To perform the Hull reformulation, `reformulation = :hull`. Variables must have bounds for the reformulation to work.
 
 ```julia
 using JuMP

examples/ex3.jl

@@ -13,7 +13,7 @@ m = Model()
         3 ≤ exp(x)
         5 ≤ x
     end,
-    reformulation=:convex_hull,
+    reformulation=:hull,
     name=:z
 )
 print(m)
@@ -22,7 +22,6 @@ print(m)
 # Subject to
 #  XOR(disj_z) : z[1] + z[2] == 1.0             <- XOR constraint
 #  x_aggregation : x - x_z1 - x_z2 == 0.0       <- aggregation of disaggregated variables
-#  disj_z[1,2] : -x_z1 - 3 z[1] <= 0.0          <- convex-hull reformulation of 2nd constraint if 1st disjunct (named disj_z[1,2] to indicate 1st disjunct, 2nd constraint)
 #  x_z1_lb : -5 z[1] - x_z1 <= 0.0              <- lower-bound constraint on disaggregated variable x_z1 (x in 1st disjunct)
 #  x_z1_ub : -10 z[1] + x_z1 <= 0.0             <- upper-bound constraint on disaggregated variable x_z1 (x in 1st disjunct)
 #  x_z2_lb : -5 z[2] - x_z2 <= 0.0              <- lower-bound constraint on disaggregated variable x_z2 (x in 2nd disjunct)

src/DisjunctiveProgramming.jl

@@ -4,12 +4,14 @@ using JuMP, IntervalArithmetic, Symbolics, Suppressor
 
 export add_disjunction!, add_proposition!, reformulate_disjunction
 export @disjunction, @proposition
+export choose!
 
 include("constraint.jl")
 include("logic.jl")
+include("bounds.jl")
 include("utils.jl")
-include("big_M.jl")
-include("convex_hull.jl")
+include("bigm.jl")
+include("hull.jl")
 include("reformulate.jl")
 include("macros.jl")
 

src/big_M.jl

@@ -0,0 +1,40 @@
+"""
+    big_m_reformulation!(constr::ConstraintRef, bin_var, M, i, j, k)
+
+Perform Big-M reformulation on a linear or quadratic constraint at index k of constraint j in disjunct i.
+
+    big_m_reformulation!(constr::NonlinearConstraintRef, bin_var, M, i, j, k)
+
+Perform Big-M reformulaiton on a nonlinear constraint at index k of constraint j in disjunct i.
+
+    big_m_reformulation!(constr::AbstractArray{<:ConstraintRef}, bin_var, M, i, j, k)
+
+Perform Big-M reformulation on a constraint at index k of constraint j in disjunct i.
+"""
+function big_m_reformulation!(constr::ConstraintRef, bin_var, M, i, j, k)
+    M = get_reform_param(M, i, j, k; constr)
+    bin_var_ref = constr.model[bin_var][i]
+    add_to_function_constant(constr, -M)
+    set_normalized_coefficient(constr, bin_var_ref , M)
+end
+function big_m_reformulation!(constr::NonlinearConstraintRef, bin_var, M, i, j, k)
+    M = get_reform_param(M, i, j, k; constr)
+    #create symbolic variables (using Symbolics.jl)
+    for var_ref in constr.model[:gdp_variable_refs]
+        var_sym = Symbol(var_ref)
+        eval(:(Symbolics.@variables($var_sym)[1]))
+    end
+    bin_var_sym = Symbol("$bin_var[$i]")
+    λ = Num(Symbolics.Sym{Float64}(bin_var_sym))
+    
+    #parse constr
+    op, lhs, rhs = parse_constraint(constr)
+    replace_Symvars!(lhs, constr.model) #convert JuMP variables into Symbolic variables
+    gx = eval(lhs) #convert the LHS of the constraint into a Symbolic expression
+    gx = gx - M*(1-λ) #add bigM
+    
+    #update constraint
+    replace_constraint(constr, gx, op, rhs)
+end
+big_m_reformulation!(constr::AbstractArray{<:ConstraintRef}, bin_var, M, i, j, k) =
+    big_m_reformulation(constr[k], bin_var, M, i, j, k)

src/bigm.jl

@@ -0,0 +1,77 @@
+"""
+    apply_interval_arithmetic(constr)
+
+Apply interval arithmetic on a constraint to find the bounds on the constraint.
+"""
+function apply_interval_arithmetic(constr)
+    #convert constraints into Expr to replace variables with interval sets and determine bounds
+    constr_type, constr_func_expr, constr_rhs = parse_constraint(constr)
+    #create a map of variables to their bounds
+    interval_map = Dict()
+    vars = all_variables(constr.model)#constr.model[:gdp_variable_constrs]
+    obj_dict = object_dictionary(constr.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
+        LB, UB = get_bounds(var, bounds_dict)
+        interval_map[string(var)] = LB..UB
+    end
+    constr_func_expr = replace_intevals!(constr_func_expr, interval_map)
+    #get bounds on the entire expression
+    func_bounds = eval(constr_func_expr)
+    Mlo = func_bounds.lo - constr_rhs
+    Mhi = func_bounds.hi - constr_rhs
+    M = constr_type == :(<=) ? Mhi : Mlo
+    isinf(M) && error("M parameter for $constr cannot be infered due to lack of variable bounds.")
+    return M
+end
+
+"""
+    get_bounds(var::VariableRef)
+
+Get bounds on a variable.
+
+    get_bounds(var, bounds_dict::Dict)
+
+Get bounds on a variable. Check if a bounds dictionary has been provided with bounds for that value.
+
+    get_bounds(var::AbstractArray{VariableRef}, bounds_dict::Dict, LB, UB)
+    
+Update lower bound `LB` and upper bound `UB` on a variable container.
+"""
+function get_bounds(var::VariableRef)
+    LB = has_lower_bound(var) ? lower_bound(var) : (is_binary(var) ? 0 : -Inf)
+    UB = has_upper_bound(var) ? upper_bound(var) : (is_binary(var) ? 1 : Inf)
+    return LB, UB
+end
+function get_bounds(var::VariableRef, bounds_dict::Dict)
+    if string(var) in keys(bounds_dict)
+        return bounds_dict[string(var)]
+    else
+        return get_bounds(var)
+    end
+end
+function get_bounds(var::AbstractArray{VariableRef}, bounds_dict::Dict, LB, UB)
+    #populate UB and LB
+    for idx in eachindex(var)
+        LB[idx], UB[idx] = get_bounds(var[idx], bounds_dict)
+    end
+    return LB, UB
+end
+function get_bounds(var::Array{VariableRef}, bounds_dict::Dict)
+    #initialize
+    LB, UB = zeros(size(var)), zeros(size(var))
+    return get_bounds(var, bounds_dict, LB, UB)
+end
+function get_bounds(var::Containers.DenseAxisArray, bounds_dict::Dict)
+    #initialize
+    LB = Containers.DenseAxisArray(zeros(size(var)), axes(var)...)
+    UB = Containers.DenseAxisArray(zeros(size(var)), axes(var)...)
+    return get_bounds(var, bounds_dict, LB, UB)
+end
+function get_gounds(var::Containers.SparseAxisArray, bounds_dict::Dict)
+    #initialize
+    idxs = keys(var.data)
+    LB = Containers.SparseAxisArray(Dict(idx => 0. for idx in idxs))
+    UB = Containers.SparseAxisArray(Dict(idx => 0. for idx in idxs))
+    return get_bounds(var, bounds_dict, LB, UB)
+end