diff --git a/src/cuttingplanes.jl b/src/cuttingplanes.jl index 2020a3c..dbc818e 100644 --- a/src/cuttingplanes.jl +++ b/src/cuttingplanes.jl @@ -1,149 +1,138 @@ -function reformulate_model( - model::JuMP.AbstractModel, - method::CuttingPlanes - ) - _clear_reformulations(model) - var_type = JuMP.variable_ref_type(model) - obj = objective_function(model) - sense = objective_sense(model) +################################################################################ +# FUNCTIONS FOR LOOP +################################################################################ - #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 collect_all_vars(model)) - main_to_rBM_map = Dict(v => rBM_ref_map[v] for v in collect_all_vars(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 +# Collect decision variables for cutting planes. Extensions may override to +# customize variable collection. +function collect_cutting_planes_vars(model::JuMP.AbstractModel) + return collect_all_vars(model) +end + +# Extract solution from a solved model (in-place). Extensions +# override for models where values live on a backend. +function extract_solution(model::JuMP.AbstractModel) + dvars = collect_cutting_planes_vars(model) + V = eltype(dvars) + T = JuMP.value_type(typeof(model)) + return Dict{V, Vector{T}}( + v => [JuMP.value(v)] for v in dvars) +end + +# Extract solution from a GDPSubmodel (SEP path). +function extract_solution(sub::GDPSubmodel) + V = eltype(sub.decision_vars) + T = JuMP.value_type(typeof(sub.model)) + sol = Dict{V, Vector{T}}() + for var in sub.decision_vars + sol[var] = JuMP.value.(sub.fwd_map[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) - _CuttingPlanes(model, rBM, main_to_rBM_map, - main_to_SEP_map, rBM_sol, SEP_sol + return sol +end + +# Set quadratic separation objective: min Σ (x_k - rBM_k)². +function _set_separation_objective( + sub::GDPSubmodel, + rBM_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}} + ) + obj_expr = zero(JuMP.GenericQuadExpr{ + JuMP.value_type(typeof(sub.model)), + JuMP.variable_ref_type(sub.model)} ) - i += 1 + for var in sub.decision_vars + sub_vars = sub.fwd_map[var] + vals = rBM_sol[var] + for k in 1:length(sub_vars) + JuMP.add_to_expression!(obj_expr, + (sub_vars[k] - vals[k]) * + (sub_vars[k] - vals[k]) + ) + end end - - #Final reformulation with added cutting planes. - reformulate_model(model, method.final_reform_method) + JuMP.@objective(sub.model, Min, obj_expr) return end -function _solve_rBM( - rBM::M, - ) where {M <: JuMP.AbstractModel} - T = JuMP.value_type(M) - optimize!(rBM, ignore_optimize_hook = true) - rBM_vars = collect_all_vars(rBM) +# Solve the separation problem. Returns (separation_obj, separation_sol). +function _solve_separation( + separation::GDPSubmodel, + rBM_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}} + ) + _set_separation_objective(separation, rBM_sol) + JuMP.optimize!(separation.model, ignore_optimize_hook = true) + separation_obj = JuMP.objective_value(separation.model) + separation_sol = extract_solution(separation) + return separation_obj, separation_sol +end - #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) +# Add cut: Σ_var Σ_k 2*(sep_k - rBM_k)*(x_k - sep_k) ≥ 0 +function add_cut( + model::JuMP.AbstractModel, + decision_vars::Vector{<:JuMP.AbstractVariableRef}, + rBM_sol::Dict{<:JuMP.AbstractVariableRef,<:Vector{<:Number}}, + separation_sol::Dict{<:JuMP.AbstractVariableRef,<:Vector{<:Number}} + ) + cut_expr = zero(JuMP.GenericAffExpr{ + JuMP.value_type(typeof(model)), + JuMP.variable_ref_type(model)}) + for var in decision_vars + rbm_vals = rBM_sol[var] + sep_vals = separation_sol[var] + for k in 1:length(rbm_vals) + xi = 2 * (sep_vals[k] - rbm_vals[k]) + JuMP.add_to_expression!(cut_expr, xi, var) + JuMP.add_to_expression!( + cut_expr, -xi * sep_vals[k]) + end end - return sol + cref = JuMP.@constraint(model, cut_expr >= 0) + push!(_reformulation_constraints(model), cref) + return 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 collect_all_vars(rBM)] +################################################################################ +# UNIFIED CUTTING PLANES LOOP +################################################################################ - #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 +function reformulate_model( + model::JuMP.AbstractModel, + method::CuttingPlanes ) - JuMP.@objective(SEP, Min, obj_expr) - optimize!(SEP, ignore_optimize_hook = true) + _clear_reformulations(model) + decision_vars = collect_cutting_planes_vars(model) - #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 + # Build separation subproblem from the clean (unreformulated) model + separation = copy_and_reformulate(model, decision_vars, Hull(), method) + JuMP.relax_integrality(separation.model) -function _CuttingPlanes( - 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 = collect_all_vars(model) + # rBM: BigM in-place, relax logical vars + reformulate_model(model, BigM(method.M_value)) + JuMP.set_optimizer(model, method.optimizer) + JuMP.set_silent(model) + relaxed_vars = relax_logical_vars(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]] - ) + # Cutting plane loop: rBM <-> SEP until convergence + for iter in 1:method.max_iter + JuMP.optimize!(model, ignore_optimize_hook = true) + rBM_sol = extract_solution(model) + separation_obj, separation_sol = _solve_separation(separation, rBM_sol) + if separation_obj <= method.seperation_tolerance + break + end + add_cut(model, decision_vars, rBM_sol, separation_sol) 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) + + unrelax_logical_vars(relaxed_vars) + _set_solution_method(model, method) + _set_ready_to_optimize(model, true) + return end ################################################################################ # ERROR MESSAGES ################################################################################ -function reformulate_model(::M, ::CuttingPlanes) where {M} +function reformulate_model( + ::M, ::CuttingPlanes + ) 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 _CuttingPlanes(::M, ::N, ::H, ::S, ::R, ::T) where {M, N, H, S, R, T} - error("_CuttingPlanes 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 866dbd2..bdcd4c2 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -481,7 +481,7 @@ along with mappings back to the original model's variables. ## Fields - `model::M`: The JuMP submodel representing a disjunct's feasible region (constraints and variable bounds). -- `dec_vars::Vector{V}`: Ordered decision variables in +- `decision_vars::Vector{V}`: Ordered decision variables in the submodel, matching the original model's ordering. - `fwd_map::Dict{V, Vector{W}}`: Forward map from original model variables to their submodel counterparts. @@ -490,7 +490,7 @@ struct GDPSubmodel{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, W <: JuMP.AbstractVariableRef} model::M - dec_vars::Vector{V} + decision_vars::Vector{V} fwd_map::Dict{V, Vector{W}} end diff --git a/src/mbm.jl b/src/mbm.jl index 91a2f87..973c051 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -434,10 +434,9 @@ end """ copy_model_with_constraints(model, constraints, method) -Build a `GDPSubmodel` with disjunct constraints passed. -This builds a model seperate from the original model with copied constraints -and variables, and maps between the original model's variables and the -submodel's variables. +Create a new model with only the variables (and their bounds) +from `model` and the selected `constraints`. Returns a +`GDPSubmodel` wrapping the new model. """ function copy_model_with_constraints( model::JuMP.AbstractModel, @@ -446,19 +445,18 @@ function copy_model_with_constraints( ) var_type = JuMP.variable_ref_type(model) sub_model = _copy_model(model) - dec_vars = collect_all_vars(model) + decision_vars = collect_all_vars(model) fwd_map = Dict{var_type, Vector{var_type}}() - for var in dec_vars + for var in decision_vars copy_var = variable_copy(sub_model, var) fwd_map[var] = [copy_var] end for cref in constraints con = JuMP.constraint_object(cref) - flat_map = Dict(v => ws[1] for (v, ws) in fwd_map) - expr = _replace_variables_in_constraint( - con.func, flat_map) + flat_map = Dict(v => only(ws) for (v, ws) in fwd_map) + expr = _replace_variables_in_constraint(con.func, flat_map) T = one(JuMP.value_type(typeof(sub_model))) JuMP.@constraint(sub_model, expr * T in con.set) end @@ -466,7 +464,7 @@ function copy_model_with_constraints( JuMP.set_optimizer(sub_model, method.optimizer) JuMP.set_silent(sub_model) - return GDPSubmodel(sub_model, dec_vars, fwd_map) + return GDPSubmodel(sub_model, decision_vars, fwd_map) end ################################################################################ diff --git a/src/utilities.jl b/src/utilities.jl index e985024..b0203d7 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -8,6 +8,100 @@ function _copy_model( return M() end +""" + copy_and_reformulate(model, decision_vars, reform_method, method) + +Copy the GDP model, reformulate the copy with `reform_method`, +and wrap in a `GDPSubmodel`. The original model is not +modified. The copy's objective is rewritten in terms of the +copied variables. +""" +function copy_and_reformulate( + model::JuMP.AbstractModel, + decision_vars::Vector{<:JuMP.AbstractVariableRef}, + reform_method::AbstractReformulationMethod, + method::AbstractReformulationMethod + ) + copy, ref_map, _ = copy_gdp_model(model) + reformulate_model(copy, reform_method) + obj = JuMP.objective_function(model) + sense = JuMP.objective_sense(model) + V = JuMP.variable_ref_type(model) + orig_to_copy = Dict{V, V}( + v => ref_map[v] for v in decision_vars) + JuMP.@objective(copy, sense, + _replace_variables_in_constraint(obj, orig_to_copy) + ) + fwd_map = Dict{V, Vector{V}}(v => [ref_map[v]] for v in decision_vars) + sub = GDPSubmodel(copy, decision_vars, fwd_map) + JuMP.set_optimizer(sub.model, method.optimizer) + JuMP.set_silent(sub.model) + return sub +end + +""" + reformulate_and_relax(model, decision_vars, reform_method, method) + +Reformulate the model in-place with `reform_method` and relax +integrality. Returns `(GDPSubmodel, undo_fn)` where `undo_fn` +restores integrality. +""" +function reformulate_and_relax( + model::JuMP.AbstractModel, + decision_vars::Vector{<:JuMP.AbstractVariableRef}, + reform_method::AbstractReformulationMethod, + method::AbstractReformulationMethod + ) + reformulate_model(model, reform_method) + V = JuMP.variable_ref_type(model) + fwd_map = Dict{V, Vector{V}}(v => [v] for v in decision_vars) + sub = GDPSubmodel(model, decision_vars, fwd_map) + JuMP.set_optimizer(sub.model, method.optimizer) + JuMP.set_silent(sub.model) + undo_relax = JuMP.relax_integrality(model) + return sub, undo_relax +end + +################################################################################ +# LOGICAL VARIABLE RELAXATION +################################################################################ +""" + relax_logical_vars(model::JuMP.AbstractModel) + +Relax the binary variables associated with logical indicators +to continuous variables in `[0, 1]`. Returns a vector of the +relaxed variable references, which can be passed to +[`unrelax_logical_vars`](@ref) to restore integrality. +""" +function relax_logical_vars(model::JuMP.AbstractModel) + V = JuMP.variable_ref_type(model) + binary_refs = V[] + for (_, bvar) in _indicator_to_binary(model) + bvar isa V || continue + push!(binary_refs, bvar) + JuMP.unset_binary(bvar) + JuMP.set_lower_bound(bvar, 0.0) + JuMP.set_upper_bound(bvar, 1.0) + end + return binary_refs +end + +""" + unrelax_logical_vars( + binary_refs::Vector{<:JuMP.AbstractVariableRef} + ) + +Restore the binary constraint on variables previously relaxed +by [`relax_logical_vars`](@ref). +""" +function unrelax_logical_vars( + binary_refs::Vector{<:JuMP.AbstractVariableRef} + ) + for v in binary_refs + JuMP.set_binary(v) + end +end + ################################################################################ # ALL VARIABLES ################################################################################ diff --git a/test/constraints/cuttingplanes.jl b/test/constraints/cuttingplanes.jl index 1264e79..9f3a793 100644 --- a/test/constraints/cuttingplanes.jl +++ b/test/constraints/cuttingplanes.jl @@ -7,8 +7,8 @@ function test_CuttingPlanes_datatype() @test method.seperation_tolerance == 1e-6 @test method.final_reform_method isa BigM @test method.M_value == 1e9 - - method = CuttingPlanes(HiGHS.Optimizer;max_iter=10, + + method = CuttingPlanes(HiGHS.Optimizer;max_iter=10, seperation_tolerance=1e-4, final_reform_method=Indicator(), M_value=1e6 ) @test method.max_iter == 10 @@ -17,24 +17,40 @@ function test_CuttingPlanes_datatype() @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()) +function test_copy_and_reformulate() + 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) + + method = CuttingPlanes(HiGHS.Optimizer) + decision_vars = DP.collect_cutting_planes_vars(model) + + # Build SEP subproblem (copy-based) + separation = DP.copy_and_reformulate(model, decision_vars, + Hull(), method) + + # GDPSubmodel with decision_vars and fwd_map map + @test separation isa DP.GDPSubmodel + @test length(separation.decision_vars) == length(decision_vars) + @test length(separation.fwd_map) == length(decision_vars) + + # Each fwd_map value is a length-1 vector + for (var, sub_vars) in separation.fwd_map + @test length(sub_vars) == 1 + @test sub_vars[1] isa JuMP.VariableRef + end + + # Subproblem is solvable + JuMP.relax_integrality(separation.model) + optimize!(separation.model, ignore_optimize_hook = true) + @test termination_status(separation.model) == MOI.OPTIMAL end -function test_solve_SEP() +function test_reformulate_and_relax() model = GDPModel() @variable(model, 0 <= x <= 100) @variable(model, Y[1:2], Logical) @@ -42,43 +58,33 @@ function test_solve_SEP() @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 = CuttingPlanes(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 + decision_vars = DP.collect_cutting_planes_vars(model) + + # Setup rBM on original model (no copy) + rBM, undo = DP.reformulate_and_relax(model, decision_vars, BigM(method.M_value), method) + + # rBM wraps the original model + @test rBM isa DP.GDPSubmodel + @test rBM.model === model + @test undo !== nothing + + # Identity forward map + @test length(rBM.fwd_map) == length(decision_vars) + for v in decision_vars + @test rBM.fwd_map[v] == [v] 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" - ) + + # Solvable with relaxed integrality + optimize!(model, ignore_optimize_hook = true) + @test termination_status(model) == MOI.OPTIMAL + + # Restore integrality + undo() end -function test_CuttingPlanes() +function test_cp_loop_helpers() model = GDPModel() @variable(model, 0 <= x <= 100) @variable(model, Y[1:2], Logical) @@ -86,45 +92,85 @@ function test_CuttingPlanes() @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 = CuttingPlanes(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._CuttingPlanes(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._CuttingPlanes(model, rBM, main_to_rBM_map, - main_to_SEP_map, rBM_sol, "not a dict" - ) + decision_vars = DP.collect_cutting_planes_vars(model) + + # Build SEP first (from clean model) + separation = DP.copy_and_reformulate(model, decision_vars, + Hull(), method) + JuMP.relax_integrality(separation.model) + + # Setup rBM on original model + rBM, undo = DP.reformulate_and_relax(model, decision_vars, BigM(method.M_value), method) + optimize!(model, ignore_optimize_hook = true) + + # Extract solution + rBM_sol = DP.extract_solution(rBM) + @test haskey(rBM_sol, x) + @test length(rBM_sol[x]) == 1 + + # Set SEP objective and solve + DP._set_separation_objective(separation, rBM_sol) + optimize!(separation.model, ignore_optimize_hook = true) + @test termination_status(separation.model) == MOI.OPTIMAL + + # SEP solution extraction + separation_sol = DP.extract_solution(separation) + @test haskey(separation_sol, x) + @test separation_sol[x][1] ≈ 4.0 atol = 0.1 + + undo() +end + +function test_cp_cut_generation() + 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) + + method = CuttingPlanes(HiGHS.Optimizer) + decision_vars = DP.collect_cutting_planes_vars(model) + + # Build SEP first (from clean model) + separation = DP.copy_and_reformulate(model, decision_vars, + Hull(), method) + JuMP.relax_integrality(separation.model) + + # Setup rBM on original model, solve + DP.reformulate_model(model, BigM(method.M_value)) + JuMP.set_optimizer(model, HiGHS.Optimizer) + JuMP.set_silent(model) + relaxed = DP.relax_logical_vars(model) + optimize!(model, ignore_optimize_hook = true) + rBM_sol = DP.extract_solution(model) + + # Solve SEP + DP._set_separation_objective(separation, rBM_sol) + optimize!(separation.model, ignore_optimize_hook = true) + separation_sol = DP.extract_solution(separation) + + # Add cut to original model + num_con_before = length(JuMP.all_constraints( + model; + include_variable_in_set_constraints = false + )) + DP.add_cut(model, decision_vars, rBM_sol, separation_sol) + num_con_after = length(JuMP.all_constraints( + model; + include_variable_in_set_constraints = false + )) + @test num_con_after == num_con_before + 1 + + # Re-solve with cut → should tighten + optimize!(model, ignore_optimize_hook = true) + rBM_sol2 = DP.extract_solution(model) + @test rBM_sol2[x][1] ≈ 4.0 atol = 0.1 + + DP.unrelax_logical_vars(relaxed) end function test_reformulate_model() @@ -141,15 +187,38 @@ function test_reformulate_model() num_con = length( JuMP.all_constraints(model; include_variable_in_set_constraints = false) ) - @test num_con == 4 + # 3 BigM constraints + 0-3 cuts depending on + # convergence (2 disjunct constraints + 1 exactly-one) + @test num_con >= 3 @test_throws ErrorException DP.reformulate_model(42, method) end +# Maximization where Hull is strictly tighter than BigM, +# forcing many CP iterations with a tight tolerance. +function test_cp_many_iterations() + model = GDPModel(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 10) + @variable(model, Y[1:2], Logical) + @constraint(model, x + y <= 5, Disjunct(Y[1])) + @constraint(model, x + y <= 8, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x + y) + cutting_planes = CuttingPlanes(HiGHS.Optimizer; + max_iter = 50, seperation_tolerance = 1e-10) + @test optimize!(model, gdp_method = cutting_planes) isa Nothing + @test termination_status(model) in + [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] +end + @testset "Cutting Planes" begin test_CuttingPlanes_datatype() - test_solve_rBM() - test_solve_SEP() - test_CuttingPlanes() + test_copy_and_reformulate() + test_reformulate_and_relax() + test_cp_loop_helpers() + test_cp_cut_generation() test_reformulate_model() + test_cp_many_iterations() end