Add add_all_cuts (#422)

docs/src/apireference.md

@@ -89,6 +89,7 @@ SDDP.SDDiP
 ```@docs
 SDDP.simulate
 SDDP.calculate_bound
+SDDP.add_all_cuts
 SDDP.Historical
 ```
 

docs/src/examples/the_farmers_problem.jl

@@ -256,4 +256,4 @@ SDDP.train(model; iteration_limit = 20)
 # Birge and Louveaux report that the optimal objective value is \\\$108,390.
 # Check that we got the correct solution using [`SDDP.calculate_bound`](@ref):
 
-@assert SDDP.calculate_bound(model) == 108_390.0
+@assert isapprox(SDDP.calculate_bound(model), 108_390.0, atol = 0.1)

docs/src/examples/vehicle_location.jl

@@ -102,7 +102,12 @@ function vehicle_location_model(integrality_handler)
             )
         end
     end
-    SDDP.train(model, iteration_limit = 50, log_frequency = 10)
+    SDDP.train(
+        model;
+        iteration_limit = 50,
+        log_frequency = 10,
+        cut_deletion_minimum = 100,
+    )
     if integrality_handler == SDDP.ContinuousRelaxation()
         @test isapprox(SDDP.calculate_bound(model), 1700, atol = 5)
     end

src/plugins/bellman_functions.jl

@@ -140,7 +140,7 @@ end
 
 # Internal function: check if the candidate point dominates the incumbent.
 function _dominates(candidate, incumbent, minimization::Bool)
-    return minimization ? candidate > incumbent : candidate < incumbent
+    return minimization ? candidate >= incumbent : candidate <= incumbent
 end
 
 # Internal function: update the Level-One datastructures inside `bellman_function`.
@@ -702,3 +702,38 @@ function read_cuts_from_file(
     end
     return
 end
+
+"""
+    add_all_cuts(model::PolicyGraph)
+
+Add all cuts that may have been deleted back into the model.
+
+## Explanation
+
+During the solve, SDDP.jl may decide to remove cuts for a variety of reasons.
+
+These can include cuts that define the optimal value function, particularly
+around the extremes of the state-space (e.g., reservoirs empty).
+
+This function ensures that all cuts discovered are added back into the model.
+
+You should call this after [`train`](@ref) and before [`simulate`](@ref).
+"""
+function add_all_cuts(model::PolicyGraph)
+    for node in values(model.nodes)
+        global_theta = node.bellman_function.global_theta
+        for cut in global_theta.cut_oracle.cuts
+            if cut.constraint_ref === nothing
+                add_cut_constraint_to_model(global_theta, cut)
+            end
+        end
+        for approximation in node.bellman_function.local_thetas
+            for cut in approximation.cut_oracle.cuts
+                if cut.constraint_ref === nothing
+                    add_cut_constraint_to_model(approximation, cut)
+                end
+            end
+        end
+    end
+    return
+end

test/plugins/bellman_functions.jl

@@ -114,3 +114,71 @@ using Test
         rm("model.cuts.json")
     end
 end
+
+@testset "add_all_cuts-SINGLE_CUT" begin
+    model = SDDP.LinearPolicyGraph(
+        stages = 3,
+        lower_bound = 0.0,
+        optimizer = GLPK.Optimizer,
+    ) do sp, t
+        @variable(sp, 5 <= x <= 15, SDDP.State, initial_value = 10)
+        @variable(sp, g >= 0)
+        @variable(sp, h >= 0)
+        @variable(sp, u >= 0)
+        @constraint(sp, inflow, x.out == x.in - h - u)
+        @constraint(sp, demand, h + g == 0)
+        @stageobjective(sp, 10 * u + g)
+        SDDP.parameterize(sp, [(0, 7.5), (5, 5.0), (10, 2.5)]) do ω
+            set_normalized_rhs(inflow, ω[1])
+            set_normalized_rhs(demand, ω[2])
+            return
+        end
+    end
+    SDDP.train(model; iteration_limit = 10)
+    for (t, node) in model.nodes
+        @test num_constraints(
+            node.subproblem,
+            AffExpr,
+            MOI.GreaterThan{Float64},
+        ) < 10
+    end
+    SDDP.add_all_cuts(model)
+    for (t, node) in model.nodes
+        n = num_constraints(node.subproblem, AffExpr, MOI.GreaterThan{Float64})
+        @test t == 3 || n == 10
+    end
+end
+
+@testset "add_all_cuts-MULTI_CUT" begin
+    model = SDDP.LinearPolicyGraph(
+        stages = 3,
+        lower_bound = 0.0,
+        optimizer = GLPK.Optimizer,
+    ) do sp, t
+        @variable(sp, 5 <= x <= 15, SDDP.State, initial_value = 10)
+        @variable(sp, g >= 0)
+        @variable(sp, h >= 0)
+        @variable(sp, u >= 0)
+        @constraint(sp, inflow, x.out == x.in - h - u)
+        @constraint(sp, demand, h + g == 0)
+        @stageobjective(sp, 10 * u + g)
+        SDDP.parameterize(sp, [(0, 7.5), (5, 5.0), (10, 2.5)]) do ω
+            set_normalized_rhs(inflow, ω[1])
+            set_normalized_rhs(demand, ω[2])
+            return
+        end
+    end
+    SDDP.train(model; iteration_limit = 10, cut_type = SDDP.MULTI_CUT)
+    for (t, node) in model.nodes
+        @test num_constraints(
+            node.subproblem,
+            AffExpr,
+            MOI.GreaterThan{Float64},
+        ) < 31
+    end
+    SDDP.add_all_cuts(model)
+    for (t, node) in model.nodes
+        n = num_constraints(node.subproblem, AffExpr, MOI.GreaterThan{Float64})
+        @test t == 3 || n == 31
+    end
+end