diff --git a/Project.toml b/Project.toml index 7bbcd68..88df6ec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DisjunctiveProgramming" uuid = "0d27d021-0159-4c7d-b4a7-9ccb5d9366cf" authors = ["hdavid16 "] -version = "0.1.6" +version = "0.1.7" [deps] IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" diff --git a/README.md b/README.md index c395910..32f7aa6 100644 --- a/README.md +++ b/README.md @@ -10,18 +10,24 @@ 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`, `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`). +After defining a JuMP model, disjunctions can be added to the model by using the `@disjunction` macro. This macro is called by `@disjunction(m, disjuncts...; kwargs...), where `disjuncts...` is a list of at least two expressions of the form: +1. A valid expression accepted by [JuMP.@constraint](https://jump.dev/JuMP.jl/stable/reference/constraints/#JuMP.@constraint). Names for the constraints or containers of constraints cannot be passed (use option 2). +2. A valid expression accepted by [JuMP.@constraints](https://jump.dev/JuMP.jl/stable/reference/constraints/#JuMP.@constraints) (using `begin...end) +3. A valid expression accepted by [JuMP.@NLconstraint](https://jump.dev/JuMP.jl/stable/reference/nlp/#JuMP.@NLconstraint). Containers of constraints cannot be passed (use option 4). Naming of non-linear constraints is not currently supported. +4. A valid expression accepted by [JuMP.@NLconstraints](https://jump.dev/JuMP.jl/stable/reference/nlp/#JuMP.@NLconstraints) (using `begin...end) +5. `Tuple` of expressions accepted by options 1 and/or 3. -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) -- The Convex-Hull (when `reformulation = :CHR` in the `@disjunction` macro) -These approaches are described [here](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities). 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("$_$name$i_lb")]` and `model[Symbol("$_$name$i_ub")]` for the lower bound and upper bound constraints, respectively. The aggregation constraint can be accessed with `model[Symbol("$_aggregation")]` When the Convex-Hull reformulation is applied to a nonlinear model, the perspective function proposed in [Furman, et al. [2020]](https://link.springer.com/article/10.1007/s10589-020-00176-0) is used. +NOTE: 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`). -When calling the `@disjunction` macro, a `name::Symbol` keyword argument can be specified to define the name of the binary indicator variable to be used for that disjunction. Otherwise, (`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($name)")]`. +The valid key-word arguments for the `@disjunction` macro are: +- `reformulation::Symbol`: `:BMR` for [Big-M Reformulation](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities#Big-M_Reformulation), `:CHR` for [Convex-Hull Reformulation](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities#Convex-Hull_Reformulation) +- `M`: Big-M value used when `reformulation = :BMR`. +- `ϵ`: 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 = :CHR`. +- `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)")]`. -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 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("$_$name$i_lb")]` and `model[Symbol("$_$name$i_ub")]` for the lower bound and upper bound constraints, respectively. The aggregation constraint can be accessed with `model[Symbol("$_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. -For empty disjuncts, use `nothing` for their positional argument (e.g., `@disjunction(m, con1, nothing, reformulation = :BMR)`). +For empty disjuncts, use `nothing` for their positional argument (e.g., `@disjunction(m, x <= 1, nothing, reformulation = :BMR)`). NOTE: `:gdp_variable_refs` and `:gdp_variable_names` are forbidden JuMP model object names when using *DisjunctiveProgramming.jl*. They are used to store the variable names and variable references in the original model. @@ -39,7 +45,7 @@ The logical proposition is then internally reformulated to an algebraic constrai The example below is from the [Northwestern University Process Optimization Open Textbook](https://optimization.mccormick.northwestern.edu/index.php/Disjunctive_inequalities). -To perform the Big-M reformulation, `:BMR` 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),(10,10))`). +To perform the Big-M reformulation, `:BMR` 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 = :CHR`. Variables must have bounds for the reformulation to work. @@ -48,29 +54,30 @@ using JuMP using DisjunctiveProgramming m = Model() -@variable(m, -1<=x<=10) - -@constraint(m, con1, 0 <= x <= 3) -@constraint(m, con2, 5 <= x <= 9) - -@disjunction(m,con1,con2,reformulation=:BMR,name=:y) +@variable(m, -5 ≤ x ≤ 10) +@disjunction( + m, + 0 ≤ x ≤ 3, + 5 ≤ x ≤ 9, + 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. - +┌ Warning: disj_y[1] : 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: disj_y[2] : 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 - 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 + XOR(disj_y) : y[1] + y[2] == 1.0 <- XOR constraint + y[1] ∨ y[2] : y[1] + y[2] >= 1.0 <- reformulated logical proposition (name is the proposition) + disj_y[1][lb] : -x + 5 y[1] <= 5.0 <- left-side of constraint in 1st disjunct (name is assigned to disj_y[1][lb]) + disj_y[1][ub] : x + 7 y[1] <= 10.0 <- right-side of constraint in 1st disjunct (name is assigned to disj_y[1][ub]) + disj_y[2][lb] : -x + 10 y[2] <= 5.0 <- left-side of constraint in 2nd disjunct (name is assigned to disj_y[2][lb]) + disj_y[2][ub] : x + y[2] <= 10.0 <- right-side of constraint in 2nd disjunct (name is assigned to disj_y[2][ub]) + x >= -5.0 <- variable lower bound + x <= 10.0 <- variable upper bound + y[1] binary <- indicator variable (1st disjunct) is binary + y[2] binary <- indicator variable (2nd disjunct) is binary ``` diff --git a/examples/ex1.jl b/examples/ex1.jl index d113ad2..0ba9a82 100644 --- a/examples/ex1.jl +++ b/examples/ex1.jl @@ -2,12 +2,29 @@ using JuMP using DisjunctiveProgramming m = Model() -@variable(m, -1<=x<=10) - -@constraint(m, con1, 0 <= x <= 3) -@constraint(m, con2, 5 <= x <= 9) - -@disjunction(m,con1,con2,reformulation=:CHR,name=:y) +@variable(m, -5 ≤ x ≤ 10) +@disjunction( + m, + 0 ≤ x ≤ 3, + 5 ≤ x ≤ 9, + reformulation=:BMR, + name=:y +) @proposition(m, y[1] ∨ y[2]) #this is a redundant proposition -print(m) \ No newline at end of file +print(m) + +# ┌ Warning: disj_y[1] : 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: disj_y[2] : 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 +# XOR(disj_y) : y[1] + y[2] == 1.0 <- XOR constraint +# y[1] ∨ y[2] : y[1] + y[2] >= 1.0 <- reformulated logical proposition (name is the proposition) +# disj_y[1][lb] : -x + 5 y[1] <= 5.0 <- left-side of constraint in 1st disjunct (name is assigned to disj_y[1][lb]) +# disj_y[1][ub] : x + 7 y[1] <= 10.0 <- right-side of constraint in 1st disjunct (name is assigned to disj_y[1][ub]) +# disj_y[2][lb] : -x + 10 y[2] <= 5.0 <- left-side of constraint in 2nd disjunct (name is assigned to disj_y[2][lb]) +# disj_y[2][ub] : x + y[2] <= 10.0 <- right-side of constraint in 2nd disjunct (name is assigned to disj_y[2][ub]) +# x >= -5.0 <- variable lower bound +# x <= 10.0 <- variable upper bound +# y[1] binary <- indicator variable (1st disjunct) is binary +# y[2] binary <- indicator variable (2nd disjunct) is binary \ No newline at end of file diff --git a/examples/ex2.jl b/examples/ex2.jl index 9ed2683..9e84021 100644 --- a/examples/ex2.jl +++ b/examples/ex2.jl @@ -3,11 +3,36 @@ using JuMP using DisjunctiveProgramming m = Model() -@variable(m, 0<=x[1:2]<=10) +@variable(m, -5 ≤ x[1:2] ≤ 10) +@disjunction( + m, + begin + con1[i=1:2], 0 ≤ x[i] ≤ [3,4][i] + end, + begin + con2[i=1:2], [5,4][i] ≤ x[i] ≤ [9,6][i] + end, + reformulation = :BMR, + name = :y +) +print(m) -@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) \ No newline at end of file +# ┌ Warning: [con1[1] : x[1] in [0.0, 3.0], con1[2] : x[2] in [0.0, 4.0]] uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound. +# ┌ Warning: [con2[1] : x[1] in [5.0, 9.0], con2[2] : x[2] in [4.0, 6.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 +# XOR(disj_y) : y[1] + y[2] == 1.0 <- XOR constraint +# con1[1][lb] : -x[1] + 5 y[1] <= 5.0 <- left-side of con1[1] +# con1[1][ub] : x[1] + 7 y[1] <= 10.0 <- right-side of con1[1] +# con1[2][lb] : -x[2] + 5 y[1] <= 5.0 <- left-side of con1[2] +# con1[2][ub] : x[2] + 6 y[1] <= 10.0 <- right-side of con1[2] +# con2[1][lb] : -x[1] + 10 y[2] <= 5.0 <- left-side of con2[1] +# con2[1][ub] : x[1] + y[2] <= 10.0 <- right-side of con2[1] +# con2[2][lb] : -x[2] + 9 y[2] <= 5.0 <- left-side of con2[2] +# con2[2][ub] : x[2] + 4 y[2] <= 10.0 <- right-side of con2[2] +# x[1] >= -5.0 <- varaible bounds +# x[2] >= -5.0 <- variable bounds +# x[1] <= 10.0 <- variable bounds +# x[2] <= 10.0 <- variable bounds +# y[1] binary <- indicator variable (1st disjunct) is binary +# y[2] binary <- indicator variable (2nd disjunct) is binary \ No newline at end of file diff --git a/examples/ex3.jl b/examples/ex3.jl index f74032c..0862a51 100644 --- a/examples/ex3.jl +++ b/examples/ex3.jl @@ -2,11 +2,40 @@ using JuMP using DisjunctiveProgramming m = Model() -@variable(m, -10 ≤ x[1:2] ≤ 10) +@variable(m, -5 ≤ x ≤ 10) +@disjunction( + m, + ( + exp(x) ≤ 2, + -3 ≤ x + ), + begin + 3 ≤ exp(x) + 5 ≤ x + end, + reformulation=:CHR, + name=:z +) +print(m) -nl_con1 = @NLconstraint(m, exp(x[1]) >= 1) -nl_con2 = @NLconstraint(m, exp(x[2]) <= 2) - -@disjunction(m, nl_con1, nl_con2, reformulation=:CHR, name=:z) - -print(m) \ No newline at end of file +# Feasibility +# 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) +# x_z2_ub : -10 z[2] + x_z2 <= 0.0 <- upper-bound constraint on disaggregated variable x_z2 (x in 2nd disjunct) +# x >= -5.0 <- lower-bound on x +# x_z1 >= -5.0 <- lower-bound on x_z1 (disaggregated x in 1st disjunct) +# x_z2 >= -5.0 <- lower-bound on x_z2 (disaggregated x in 2nd disjunct) +# x <= 10.0 <- upper-bound on x +# x_z1 <= 10.0 <- upper-bound on x_z1 (disaggregated x in 1st disjunct) +# x_z2 <= 10.0 <- upper-bound on x_z2 (disaggregated x in 2nd disjunct) +# z[1] binary <- indicator variable (1st disjunct) is binary +# z[2] binary <- indicator variable (2nd disjunct) is binary +# 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_z2 + 5.0 * z[2] <= 0 \ No newline at end of file diff --git a/src/DisjunctiveProgramming.jl b/src/DisjunctiveProgramming.jl index da89e24..1d94395 100644 --- a/src/DisjunctiveProgramming.jl +++ b/src/DisjunctiveProgramming.jl @@ -2,7 +2,7 @@ module DisjunctiveProgramming using JuMP, IntervalArithmetic, Symbolics, Suppressor -export add_disjunction, add_proposition +export add_disjunction!, add_proposition!, reformulate_disjunction export @disjunction, @proposition include("constraint.jl") @@ -11,6 +11,6 @@ include("utils.jl") include("big_M.jl") include("convex_hull.jl") include("reformulate.jl") -include("macro.jl") +include("macros.jl") end # module diff --git a/src/big_M.jl b/src/big_M.jl index 7e5b096..0c0c469 100644 --- a/src/big_M.jl +++ b/src/big_M.jl @@ -1,4 +1,4 @@ -function BMR!(m, constr, bin_var, i, k, M) +function BMR!(constr, bin_var, i, k, M) if ismissing(k) ref = constr else diff --git a/src/constraint.jl b/src/constraint.jl index d292ae3..63406bb 100644 --- a/src/constraint.jl +++ b/src/constraint.jl @@ -2,23 +2,6 @@ is_interval_constraint(con_ref::ConstraintRef{<:AbstractModel}) = constraint_obj 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 @@ -55,8 +38,10 @@ function check_constraint!(m, 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) + 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." + delete_original_constraint!(m, constr) + end return new_constr end diff --git a/src/convex_hull.jl b/src/convex_hull.jl index 9e64a4c..bfb468f 100644 --- a/src/convex_hull.jl +++ b/src/convex_hull.jl @@ -1,4 +1,4 @@ -function CHR!(m, constr, bin_var, i, k, eps) +function CHR!(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 @@ -103,11 +103,7 @@ function linear_perspective_function(ref, bin_var, i) #replace each variable with its disaggregated version for var_ref in ref.model[:gdp_variable_refs] #get disaggregated variable reference - if occursin("[", string(var_ref)) - var_name_i = replace(string(var_ref), "[" => "_$bin_var$i[") - else - var_name_i = "$(var_ref)_$bin_var$i" - end + var_name_i = name_disaggregated(var_ref, bin_var, i) var_i_ref = variable_by_name(ref.model, var_name_i) #check var_ref is present in the constraint coeff = normalized_coefficient(ref, var_ref) @@ -124,14 +120,10 @@ end function nonlinear_perspective_function(ref, bin_var, i, eps) #create symbolic variables (using Symbolics.jl) - sym_vars = Dict() - for var_ref in ref.model[:gdp_variable_refs] - #get disaggregated variable reference - var_name_i = replace(string(var_ref), "[" => "_$bin_var$i[") - var_sym = Symbol(var_ref) - var_i_sym = Symbol(var_name_i) - sym_vars[eval(:(Symbolics.@variables($var_sym)[1]))] = eval(:(Symbolics.@variables($var_i_sym)[1])) - end + sym_vars = Dict( + symbolic_variable(var_ref) => symbolic_variable(name_disaggregated(var_ref, bin_var, i)) + for var_ref in ref.model[:gdp_variable_refs] + ) ϵ = eps #epsilon parameter for perspective function (See Furman, Sawaya, Grossmann [2020] perspecive function) bin_var_sym = Symbol("$bin_var[$i]") λ = Num(Symbolics.Sym{Float64}(bin_var_sym)) @@ -147,7 +139,9 @@ function nonlinear_perspective_function(ref, bin_var, i, eps) #first term g1 = FSG1*substitute(gx, Dict(var => var_i/FSG1 for (var,var_i) in sym_vars)) #second term - g2 = FSG2*substitute(gx, Dict(var => 0 for var in keys(sym_vars))) + g0 = substitute(gx, Dict(var => 0 for var in keys(sym_vars))) + @assert !isinf(g0.val) "Convex-hull reformulation has failed for non-linear constraint $ref: $gx is not defined at 0. Perspective function is undetermined." + g2 = FSG2*g0 #create perspective function and simplify pers_func = simplify(g1 - g2, expand = true) #replace FSG expressions & simplify diff --git a/src/gdp.jl b/src/gdp.jl new file mode 100644 index 0000000..384ce3f --- /dev/null +++ b/src/gdp.jl @@ -0,0 +1,9 @@ +abstract type :ModelVariables end + +function GDPModel() + model = Model() + model.ext[:gdp_variable_refs] = [] + model.ext[:gdp_variable_names] = [] + + return model +end \ No newline at end of file diff --git a/src/macro.jl b/src/macro.jl deleted file mode 100644 index 1a5aca6..0000000 --- a/src/macro.jl +++ /dev/null @@ -1,53 +0,0 @@ -macro disjunction(args...) - pos_args, kw_args, _ = Containers._extract_kw_args(args) - - #get kw_args - reformulation = filter(i -> i.args[1] == :reformulation, kw_args) - if !isempty(reformulation) - reformulation = reformulation[1].args[2] - else - throw(UndefKeywordError(:reformulation)) - end - M = filter(i -> i.args[1] == :M, kw_args) - M = !isempty(M) ? esc(M[1].args[2]) : :(missing) - eps = filter(i -> i.args[1] == :ϵ, kw_args) - eps = !isempty(eps) ? esc(eps[1].args[2]) : :(1e-6) - name = filter(i -> i.args[1] == :name, kw_args) - name = !isempty(name) ? esc(name[1].args[2]) : :(missing) - - #get args - m = esc(pos_args[1]) - disj = [esc(a) for a in pos_args[2:end]] - - #build disjunction - :(add_disjunction($m, $(disj...), reformulation = $reformulation, M = $M, ϵ = $eps, name = $name)) -end - -function add_disjunction(m::Model,disj...;reformulation::Symbol,M=missing,ϵ=1e-6,name=missing) - @assert m isa Model "A valid JuMP Model must be provided." - @assert reformulation in [:BMR, :CHR] "Invalid reformulation method passed to keyword argument `:reformulation`. Valid options are :BMR (Big-M Reformulation) and :CHR (Convex-Hull Reormulation)." - @assert length(disj) > 1 "At least 2 disjuncts must be included." - #create binary indicator variables for each disjunction - disj_name = ismissing(name) ? Symbol("disj",gensym()) : name - @assert !in(disj_name, keys(object_dictionary(m))) "The disjunction name $disj_name already exists as a model object. Specify a name that is not present in the model's object_dictionary." - #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 - 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) -end - -macro proposition(m, expr) - #get args - m = esc(m) - expr = Expr(:quote, expr) - :(add_proposition($m, $expr)) -end - -function add_proposition(m::Model, expr::Expr) - @assert m isa Model "A valid JuMP Model must be provided." - to_cnf!(m, expr) -end \ No newline at end of file diff --git a/src/macros.jl b/src/macros.jl new file mode 100644 index 0000000..746ca30 --- /dev/null +++ b/src/macros.jl @@ -0,0 +1,124 @@ +macro disjunction(m, args...) + #get disjunction (pos_args) and keyword arguments + pos_args, kw_args, _ = Containers._extract_kw_args(args) + @assert length(pos_args) > 1 "At least 2 disjuncts must be included. If there is an empty disjunct, use `nothing`." + + #get kw_args and set defaults if missing + reformulation = filter(i -> i.args[1] == :reformulation, kw_args) + if !isempty(reformulation) + reformulation = reformulation[1].args[2] + reformulation_kind = eval(reformulation) + @assert reformulation_kind in [:BMR, :CHR] "Invalid reformulation method passed to keyword argument `:reformulation`. Valid options are :BMR (Big-M Reformulation) and :CHR (Convex-Hull Reormulation)." + if reformulation_kind == :BMR + M = filter(i -> i.args[1] == :M, kw_args) + param = !isempty(M) ? M[1].args[2] : :(missing) + elseif reformulation_kind == :CHR + ϵ = filter(i -> i.args[1] == :ϵ, kw_args) + param = !isempty(ϵ) ? ϵ[1].args[2] : :(1e-6) + else + end + else + throw(UndefKeywordError(:reformulation)) + end + name = filter(i -> i.args[1] == :name, kw_args) + if !isempty(name) + name = name[1].args[2] + disj_name = Symbol("disj_",eval(name)) + else + name = :(Symbol("disj",gensym())) + disj_name = eval(name) + end + + #create constraints for each disjunction + disj_names = [Symbol("$(disj_name)[$i]") for i in eachindex(pos_args)] + disjunction = [] + for (d,dname) in zip(pos_args,disj_names) + if Meta.isexpr(d, :tuple) + for (j,di) in enumerate(d.args) + i = findfirst(x -> x == d, pos_args) + dname_j = Symbol("$(disj_name)[$i,$j]") + d.args[j] = add_disjunction_constraint(m, di, dname_j) + end + push!(disjunction, d) + else + push!(disjunction, add_disjunction_constraint(m, d, dname)) + end + end + + #XOR constraint name + xor_con = Symbol("XOR($disj_name)") + + #build disjunction + code = quote + @assert !in($name, keys(object_dictionary($m))) "The disjunction name $name is already registered in the model. Specify new name." + $m[$name] = @variable($m, [$eachindex($disjunction)], Bin, base_name = string($name)) + @constraint($m, $xor_con, sum($m[$name]) == 1) + reformulate_disjunction($m, $(disjunction...); bin_var = $name, reformulation = $reformulation, param = $param) + end + + return esc(code) +end + +function add_disjunction_constraint(m, d, dname) + if Meta.isexpr(d, :block) + d = quote + try + @constraints($m,$d) + catch e + if e isa ErrorException + @NLconstraints($m,$d) + else + throw(e) + end + end + end + elseif Meta.isexpr(d, (:call, :comparison)) + d = quote + try + @constraint($m,$dname,$d) + catch e + if e isa ErrorException + @NLconstraint($m,$d) + else + throw(e) + end + end + end + end + + return d +end + +function add_disjunction!(m::Model,disj...;reformulation::Symbol,M=missing,ϵ=1e-6,name=missing) + #run checks + @assert reformulation in [:BMR, :CHR] "Invalid reformulation method passed to keyword argument `:reformulation`. Valid options are :BMR (Big-M Reformulation) and :CHR (Convex-Hull Reormulation)." + @assert length(disj) > 1 "At least 2 disjuncts must be included. If there is an empty disjunct, use `nothing`." + + #create binary indicator variables for each disjunction + bin_var = ismissing(name) ? Symbol("disj",gensym()) : name + @assert !in(bin_var, keys(object_dictionary(m))) "The disjunction name $bin_var is already registered in the model. Specify new name." + + #create indicator variable + m[bin_var] = @variable(m, [eachindex(disj)], Bin, base_name = string(bin_var)) + + #add xor constraint on binary variable + xor_con = "XOR($bin_var)" + m[Symbol(xor_con)] = @constraint(m, sum(m[bin_var]) == 1, base_name = xor_con) + + #apply reformulation + param = reformulation == :BMR ? M : ϵ + reformulate_disjunction(m, disj; bin_var, reformulation, param) +end + +macro proposition(m, expr) + #get args + expr = QuoteNode(expr) + code = :(add_proposition!($m, $expr)) + + return esc(code) +end + +function add_proposition!(m::Model, expr::Expr) + @assert m isa Model "A valid JuMP Model must be provided." + to_cnf!(m, expr) +end diff --git a/src/reformulate.jl b/src/reformulate.jl index 42218c8..d76e286 100644 --- a/src/reformulate.jl +++ b/src/reformulate.jl @@ -1,6 +1,6 @@ -function reformulate_disjunction(m, disj, bin_var, reformulation, param) +function reformulate_disjunction(m::Model, disj...; bin_var, reformulation, param) #check disj - disj = check_disjunction(m, disj) + disj = check_disjunction!(m, disj) #get original variable refs and variable names vars = setdiff(all_variables(m), m[bin_var]) var_names = unique(Symbol.([split("$var","[")[1] for var in vars])) @@ -15,39 +15,56 @@ function reformulate_disjunction(m, disj, bin_var, reformulation, param) disaggregate_variables(m, disj, bin_var) sum_disaggregated_variables(m, disj, bin_var) end - reformulate(m, disj, bin_var, reformulation, param) + reformulate(disj, bin_var, reformulation, param) # return m[bin_var] end -function reformulate(m, disj, bin_var, reformulation, param) +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 reformulate(disj, bin_var, reformulation, param) for (i,constr) in enumerate(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 for (j,constr_j) in enumerate(constr) - apply_reformulation(m, constr_j, bin_var, reformulation, param, i, j) + apply_reformulation(constr_j, bin_var, reformulation, param, i, j) end elseif constr isa Union{ConstraintRef, Array, Containers.DenseAxisArray, Containers.SparseAxisArray} - apply_reformulation(m, constr, bin_var, reformulation, param, i) + apply_reformulation(constr, bin_var, reformulation, param, i) end end end -function apply_reformulation(m, constr, bin_var, reformulation, param, i, j = missing) +function apply_reformulation(constr, bin_var, reformulation, param, i, j = missing) param = get_reform_param(param, i, j) #M or eps if constr isa ConstraintRef - call_reformulation(reformulation, m, constr, bin_var, i, missing, param) + call_reformulation(reformulation, constr, bin_var, i, missing, param) elseif constr isa Union{Array, Containers.DenseAxisArray, Containers.SparseAxisArray} for k in eachindex(constr) - call_reformulation(reformulation, m, constr, bin_var, i, k, param) + call_reformulation(reformulation, constr, bin_var, i, k, param) end end end -function call_reformulation(reformulation, m, constr, bin_var, i, k, param) +function call_reformulation(reformulation, constr, bin_var, i, k, param) if reformulation == :BMR - BMR!(m, constr, bin_var, i, k, param) + BMR!(constr, bin_var, i, k, param) elseif reformulation == :CHR - CHR!(m, constr, bin_var, i, k, param) + CHR!(constr, bin_var, i, k, param) end end diff --git a/src/utils.jl b/src/utils.jl index 064d380..dbb7319 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -70,6 +70,7 @@ function replace_NLconstraint(ref, sym_expr, op, rhs) #operator (not the symbol). #This is done to later replace the symbolic variables with JuMP variables, #without messing with the math operators. + #NOTE: Maybe not necessary anymore due to implementation of replace_JuMPvars if ref isa NonlinearConstraintRef expr = Base.remove_linenums!(build_function(sym_expr)).args[2].args[1] elseif ref isa ConstraintRef @@ -98,16 +99,15 @@ end function replace_Symvars!(expr, model; logical_proposition = false) #replace JuMP variables with symbolic variables + name = join(split(string(expr)," ")) + var = variable_by_name(model, name) + if !isnothing(var) + logical_proposition && @assert is_binary(var) "Only binary variables are allowed in $expr." + return Symbol(name) + end if expr isa Expr - name = join(split(string(expr)," ")) - var = variable_by_name(model, name) - if !isnothing(var) - logical_proposition && @assert is_binary(var) "Only binary variables are allowed in $expr." - expr = Symbol(name) - else - for i in eachindex(expr.args) - expr.args[i] = replace_Symvars!(expr.args[i], model) - end + for i in eachindex(expr.args) + expr.args[i] = replace_Symvars!(expr.args[i], model) end end @@ -155,4 +155,41 @@ function replace_intevals!(expr, intervals) end return expr -end \ No newline at end of file +end + +function symbolic_variable(var_ref) + var_sym = Symbol(var_ref) + return eval(:(Symbolics.@variables($var_sym)[1])) +end + +function name_disaggregated(var_ref, bin_var, i) + #get disaggregated variable reference + if occursin("[", string(var_ref)) + var_name_i = replace(string(var_ref), "[" => "_$bin_var$i[") + else + var_name_i = "$(var_ref)_$bin_var$i" + end + + return var_name_i +end + +# function is_linear_func(expr::Expr, m) +# m = eval(m) +# if expr.head == :call +# func = expr.args[2] isa Number ? expr.args[3] : expr.args[2] +# elseif expr.head == :comparison +# func = expr.args[3] +# end +# sym_vars = symbolic_variable.(all_variables(m)) +# func = eval(replace_Symvars!(func, m)) +# try +# # eval(func) +# replace_JuMPvars!(func, m) +# return true +# catch e +# if e isa ErrorException +# return false +# end +# end +# # Symbolics.islinear(func, sym_vars) +# end \ No newline at end of file