Support cutting planes technique.

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
 

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")

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

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
 

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

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")
+include("solve.jl")