diff --git a/README.md b/README.md index 5d68bae..f64fa6b 100644 --- a/README.md +++ b/README.md @@ -184,6 +184,15 @@ The following reformulation methods are currently supported: All variables must be included in exactly one partition. For manual partitioning, ensure each variable appears in exactly one group. For automatic partitioning, variables are divided as evenly as possible among the specified number of partitions. +6. [Cutting Planes](https://pubsonline.informs.org/doi/10.1287/ijoc.2015.0669): This method iteratively generates cutting planes using a separation problem and a relaxed Big-M formulation, then applies a final reformulation method. The `cutting_planes` struct is created with the following arguments: + + - `optimizer`: Optimizer to use when solving the separation and relaxed Big-M subproblems. This is a required value. + - `max_iter`: Maximum number of cutting plane iterations. Default: `3`. + - `seperation_tolerance`: Convergence tolerance for the separation problem objective. Default: `1e-6`. + - `final_reform_method`: Reformulation method to apply after cutting plane iterations. Default: `BigM()`. + - `M_value`: Big-M value to use in the relaxed Big-M reformulation during iterations. Default: `1e9`. + + ## Release Notes diff --git a/src/DisjunctiveProgramming.jl b/src/DisjunctiveProgramming.jl index 3bf91c8..5dd7e3b 100644 --- a/src/DisjunctiveProgramming.jl +++ b/src/DisjunctiveProgramming.jl @@ -20,6 +20,7 @@ include("constraints.jl") include("macros.jl") include("reformulate.jl") include("bigm.jl") +include("cuttingplanes.jl") include("hull.jl") include("mbm.jl") include("indicator.jl") diff --git a/src/cuttingplanes.jl b/src/cuttingplanes.jl new file mode 100644 index 0000000..5644a7b --- /dev/null +++ b/src/cuttingplanes.jl @@ -0,0 +1,149 @@ +function reformulate_model( + model::JuMP.AbstractModel, + method::cutting_planes + ) + _clear_reformulations(model) + var_type = JuMP.variable_ref_type(model) + obj = objective_function(model) + sense = objective_sense(model) + + #Creation of seperation (SEP) and relaxed big M model (rBM). + SEP, sep_ref_map, _ = copy_gdp_model(model) + rBM, rBM_ref_map, _ = copy_gdp_model(model) + reformulate_model(rBM, BigM(method.M_value)) + reformulate_model(SEP, Hull()) + main_to_SEP_map = Dict(v => sep_ref_map[v] for v in all_variables(model)) + main_to_rBM_map = Dict(v => rBM_ref_map[v] for v in all_variables(model)) + JuMP.set_optimizer(SEP, method.optimizer) + JuMP.set_optimizer(rBM, method.optimizer) + JuMP.set_silent(rBM) + JuMP.set_silent(SEP) + JuMP.relax_integrality(rBM) + JuMP.relax_integrality(SEP) + JuMP.@objective(rBM, sense, + _replace_variables_in_constraint(obj, main_to_rBM_map) + ) + #Mapping of variables between models. + rBM_to_SEP_map = Dict{var_type, var_type}() + SEP_to_rBM_map = Dict{var_type, var_type}() + for (var, rBM_var) in main_to_rBM_map + SEP_var = main_to_SEP_map[var] + rBM_to_SEP_map[rBM_var] = SEP_var + SEP_to_rBM_map[SEP_var] = rBM_var + end + + #Main cutting planes loop. + i = 1 + sep_obj = Inf + while i <= method.max_iter && sep_obj > method.seperation_tolerance + rBM_sol = _solve_rBM(rBM) + SEP_sol = _solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map, rBM_to_SEP_map) + sep_obj = objective_value(SEP) + _cutting_planes(model, rBM, main_to_rBM_map, + main_to_SEP_map, rBM_sol, SEP_sol + ) + i += 1 + end + + #Final reformulation with added cutting planes. + reformulate_model(model, method.final_reform_method) + return +end + +function _solve_rBM( + rBM::M, + ) where {M <: JuMP.AbstractModel} + T = JuMP.value_type(M) + optimize!(rBM, ignore_optimize_hook = true) + rBM_vars = JuMP.all_variables(rBM) + + #Solution to be passed to SEP model. + sol = Dict{JuMP.AbstractVariableRef,T}(var => zero(T) for var in rBM_vars) + for rBM_var in rBM_vars + sol[rBM_var] = JuMP.value(rBM_var) + end + return sol +end + +function _solve_SEP( + SEP::M, + rBM::M, + rBM_sol::Dict{<:JuMP.AbstractVariableRef,T}, + SEP_to_rBM_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef}, + rBM_to_SEP_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef} + ) where {M <: JuMP.AbstractModel, T <: Number} + + SEP_vars = [rBM_to_SEP_map[rBM_var] for rBM_var in JuMP.all_variables(rBM)] + + #Modified objective function for SEP. + obj_expr = sum( + (SEP_var - rBM_sol[SEP_to_rBM_map[SEP_var]])^2 for SEP_var in SEP_vars + ) + JuMP.@objective(SEP, Min, obj_expr) + optimize!(SEP, ignore_optimize_hook = true) + + #Solution to be used in cutting plane generation. + sol = Dict{JuMP.AbstractVariableRef, T}(var => zero(T) for var in SEP_vars) + for SEP_var in SEP_vars + sol[SEP_var] = JuMP.value(SEP_var) + end + return sol +end + +function _cutting_planes( + model::M, + rBM::M, + main_to_rBM_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef}, + main_to_SEP_map::Dict{<:JuMP.AbstractVariableRef,<:JuMP.AbstractVariableRef}, + rBM_sol::Dict{<:JuMP.AbstractVariableRef,T}, + SEP_sol::Dict{<:JuMP.AbstractVariableRef,T}, + ) where {M <: JuMP.AbstractModel, T <: Number} + main_vars = JuMP.all_variables(model) + + #Cutting plane generation + ξ_sep = Dict{JuMP.AbstractVariableRef,T}(var =>zero(T) for var in main_vars) + for var in main_vars + ξ_sep[var] = 2*(SEP_sol[main_to_SEP_map[var]] + -rBM_sol[main_to_rBM_map[var]] + ) + end + #Cutting plane added to main model. + main_cut = JuMP.@expression(model, + sum(ξ_sep[var]*(var - SEP_sol[main_to_SEP_map[var]]) + for var in main_vars + ) + ) + #Cutting plane added to rBM + rBM_cut = _replace_variables_in_constraint(main_cut, main_to_rBM_map) + JuMP.@constraint(model, main_cut >= 0.0) + JuMP.@constraint(rBM, rBM_cut >= 0.0) +end + +################################################################################ +# ERROR MESSAGES +################################################################################ + +function reformulate_model(::M, ::cutting_planes) where {M} + error("reformulate_model not implemented for model type `$(M)`.") +end + +function _solve_rBM(::M) where {M} + error("_solve_rBM not implemented for model type `$(M)`.") +end + +function _solve_SEP(::M, ::N, ::H, ::S, ::R) where {M, N, H, S, R} + error("_solve_SEP not implemented for argument types:\n + SEP: `$(M)`, rBM: `$(N)`,\n + rBM_sol: `$(H)`,\n + SEP_to_rBM_map: `$(S)`,\n + rBM_to_SEP_map: `$(R)`.") +end + +function _cutting_planes(::M, ::N, ::H, ::S, ::R, ::T) where {M, N, H, S, R, T} + error("_cutting_planes not implemented for argument types: \n + model: `$(M)`, rBM: `$(N)`,\n + main_to_rBM_map: `$(H)`, main_to_SEP_map: + `$(S)`,\n + rBM_sol: `$(R)`,\n + SEP_sol: `$(T)`.") +end diff --git a/src/datatypes.jl b/src/datatypes.jl index d46f4e0..3a5f09d 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -431,6 +431,36 @@ mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationM end end +""" + cutting_planes{O} <: AbstractReformulationMethod + +A type for using the cutting planes approach for disjunctive constraints. + +**Fields** +- `optimizer::O`: Optimizer to use when solving mini-models (required). +- `max_iter::Int`: Number of iterations (default = `3`). +- `seperation_tolerance::Float64`: Tolerance for the separation problem (default = `1e-6`). +- `final_reform_method::AbstractReformulationMethod`: Final reformulation +method to use after cutting planes (default = `BigM()`). +- `M_value::Float64`: Big-M value to use in the final reformulation (default = `1e9`). +""" +struct cutting_planes{O} <: AbstractReformulationMethod + optimizer::O; + max_iter::Int + seperation_tolerance::Float64 + final_reform_method::AbstractReformulationMethod + M_value::Float64 + function cutting_planes( + optimizer::O; + max_iter::Int = 3, + seperation_tolerance::Float64 = 1e-6, + final_reform_method = BigM(), + M_value::Float64 = 1e9 + ) where {O} + new{O}(optimizer, max_iter, seperation_tolerance, final_reform_method, M_value) + end +end + """ PSplit <: AbstractReformulationMethod diff --git a/test/constraints/cuttingplanes.jl b/test/constraints/cuttingplanes.jl new file mode 100644 index 0000000..3c5239e --- /dev/null +++ b/test/constraints/cuttingplanes.jl @@ -0,0 +1,155 @@ +using HiGHS + +function test_cutting_planes_datatype() + method = cutting_planes(HiGHS.Optimizer) + @test method.optimizer == HiGHS.Optimizer + @test method.max_iter == 3 + @test method.seperation_tolerance == 1e-6 + @test method.final_reform_method isa BigM + @test method.M_value == 1e9 + + method = cutting_planes(HiGHS.Optimizer;max_iter=10, + seperation_tolerance=1e-4, final_reform_method=Indicator(), M_value=1e6 + ) + @test method.max_iter == 10 + @test method.seperation_tolerance == 1e-4 + @test method.final_reform_method isa Indicator + @test method.M_value == 1e6 +end + +function test_solve_rBM() + rBM = JuMP.Model(HiGHS.Optimizer) + @variable(rBM, 0 <= x <= 100) + @variable(rBM, 0 <= y[1:2] <= 1) + @constraint(rBM, x <= 3 + 100(1 - y[1])) + @constraint(rBM, x <= 4 + 100(1 - y[2])) + @constraint(rBM, y[1] + y[2] == 1) + @objective(rBM, Max, x) + + solutions = DP._solve_rBM(rBM) + @test solutions[x] == 53.5 + @test solutions[y[1]] == 0.495 + @test solutions[y[2]] == 0.505 + + @test_throws ErrorException DP._solve_rBM(Dict()) +end + +function test_solve_SEP() + model = GDPModel() + @variable(model, 0 <= x <= 100) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 4, Disjunct(Y[2])) + @disjunction(model, [Y[1], Y[2]]) + @objective(model, Max, x) + var_type = JuMP.variable_ref_type(model) + method = cutting_planes(HiGHS.Optimizer) + obj = objective_function(model) + sense = objective_sense(model) + SEP, sep_ref_map, _ = DP.copy_gdp_model(model) + rBM, rBM_ref_map, _ = DP.copy_gdp_model(model) + DP.reformulate_model(rBM, DP.BigM(method.M_value)) + DP.reformulate_model(SEP, DP.Hull()) + main_to_SEP_map = Dict(v => sep_ref_map[v] for v in all_variables(model)) + main_to_rBM_map = Dict(v => rBM_ref_map[v] for v in all_variables(model)) + JuMP.set_optimizer(SEP, method.optimizer) + JuMP.set_optimizer(rBM, method.optimizer) + JuMP.set_silent(rBM) + JuMP.set_silent(SEP) + JuMP.relax_integrality(rBM) + JuMP.relax_integrality(SEP) + JuMP.@objective(rBM, sense, + DP._replace_variables_in_constraint(obj, main_to_rBM_map) + ) + rBM_to_SEP_map = Dict{var_type, var_type}() + SEP_to_rBM_map = Dict{var_type, var_type}() + for (var, rBM_var) in main_to_rBM_map + SEP_var = main_to_SEP_map[var] + rBM_to_SEP_map[rBM_var] = SEP_var + SEP_to_rBM_map[SEP_var] = rBM_var + end + rBM_sol = DP._solve_rBM(rBM) + SEP_sol = DP._solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map, rBM_to_SEP_map) + @test length(SEP_sol) == length(rBM_sol) + @test SEP_sol[rBM_to_SEP_map[main_to_rBM_map[x]]] ≈ 4.0 + + @test_throws ErrorException DP._solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map + , "not a dict" + ) +end + +function test_cutting_planes() + model = GDPModel() + @variable(model, 0 <= x <= 100) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 4, Disjunct(Y[2])) + @disjunction(model, [Y[1], Y[2]]) + @objective(model, Max, x) + var_type = JuMP.variable_ref_type(model) + method = cutting_planes(HiGHS.Optimizer) + obj = objective_function(model) + sense = objective_sense(model) + SEP, sep_ref_map, _ = DP.copy_gdp_model(model) + rBM, rBM_ref_map, _ = DP.copy_gdp_model(model) + DP.reformulate_model(rBM, DP.BigM(method.M_value)) + DP.reformulate_model(SEP, DP.Hull()) + main_to_SEP_map = Dict(v => sep_ref_map[v] for v in all_variables(model)) + main_to_rBM_map = Dict(v => rBM_ref_map[v] for v in all_variables(model)) + JuMP.set_optimizer(SEP, method.optimizer) + JuMP.set_optimizer(rBM, method.optimizer) + JuMP.set_silent(rBM) + JuMP.set_silent(SEP) + JuMP.relax_integrality(rBM) + JuMP.relax_integrality(SEP) + JuMP.@objective(rBM, sense, + DP._replace_variables_in_constraint(obj, main_to_rBM_map) + ) + rBM_to_SEP_map = Dict{var_type, var_type}() + SEP_to_rBM_map = Dict{var_type, var_type}() + for (var, rBM_var) in main_to_rBM_map + SEP_var = main_to_SEP_map[var] + rBM_to_SEP_map[rBM_var] = SEP_var + SEP_to_rBM_map[SEP_var] = rBM_var + end + rBM_sol = DP._solve_rBM(rBM) + SEP_sol = DP._solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map, rBM_to_SEP_map) + DP._cutting_planes(model, rBM, main_to_rBM_map, main_to_SEP_map, rBM_sol, SEP_sol) + + rBM_sol = DP._solve_rBM(rBM) + SEP_sol = DP._solve_SEP(SEP, rBM, rBM_sol, SEP_to_rBM_map, rBM_to_SEP_map) + + @test rBM_sol[main_to_rBM_map[x]] ≈ 4.0 + @test SEP_sol[rBM_to_SEP_map[main_to_rBM_map[x]]] ≈ 4.0 atol=1e-3 + + @test_throws ErrorException DP._cutting_planes(model, rBM, main_to_rBM_map, + main_to_SEP_map, rBM_sol, "not a dict" + ) +end + +function test_reformulate_model() + model = GDPModel() + @variable(model, 0 <= x[1:4] <= 100) + @variable(model, Y[1:2], Logical) + @constraint(model, x[1] + x[2] <= 3, Disjunct(Y[1])) + @constraint(model, x[3] + x[4] <= 4, Disjunct(Y[2])) + @disjunction(model, [Y[1], Y[2]]) + @objective(model, Max, x[1] + x[2]) + + method = cutting_planes(HiGHS.Optimizer) + DP.reformulate_model(model, method) + num_con = length( + JuMP.all_constraints(model; include_variable_in_set_constraints = false) + ) + @test num_con == 4 + @test_throws ErrorException DP.reformulate_model(42, method) +end + + +@testset "Cutting Planes" begin + test_cutting_planes_datatype() + test_solve_rBM() + test_solve_SEP() + test_cutting_planes() + test_reformulate_model() +end diff --git a/test/runtests.jl b/test/runtests.jl index 843a95f..283cbb2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,8 +17,9 @@ include("constraints/indicator.jl") include("constraints/mbm.jl") include("constraints/bigm.jl") include("constraints/psplit.jl") +include("constraints/cuttingplanes.jl") include("constraints/hull.jl") include("constraints/fallback.jl") include("constraints/disjunction.jl") include("print.jl") -include("solve.jl") \ No newline at end of file +include("solve.jl") diff --git a/test/solve.jl b/test/solve.jl index bb503bd..b8b8db3 100644 --- a/test/solve.jl +++ b/test/solve.jl @@ -54,15 +54,25 @@ function test_linear_gdp_example(m, use_complements = false) @test !value(W[1]) @test !value(W[2]) + @test optimize!(m, gdp_method = cutting_planes(HiGHS.Optimizer)) isa Nothing + @test termination_status(m) == MOI.OPTIMAL + @test objective_value(m) ≈ 11 atol=1e-3 + @test value.(x) ≈ [9,2] atol=1e-3 + @test !value(Y[1]) + @test value(Y[2]) + @test !value(W[1]) + @test !value(W[2]) + + m_copy, ref_map = JuMP.copy_model(m) lv_map = DP.copy_gdp_data(m, m_copy, ref_map) set_optimizer(m_copy, HiGHS.Optimizer) set_optimizer_attribute(m_copy, "output_flag", false) optimize!(m_copy, gdp_method = BigM()) @test termination_status(m_copy) == MOI.OPTIMAL - @test objective_value(m_copy) ≈ 11 - @test value.(x) == value.(ref_map[x]) - @test value.(Y[1]) == value.(lv_map[Y[1]]) + @test objective_value(m_copy) ≈ 11 atol=1e-3 + @test value.(x) ≈ value.(ref_map[x]) atol=1e-3 + @test value.(Y[1]) == value.(lv_map[Y[1]]) @test value.(Y[2]) == value.(lv_map[Y[2]]) @test !value(W[1]) @test !value(W[2])