[DEV] Add functions to prune cuts

Author: frapac

src/SDDPoptimize.jl

@@ -504,3 +504,77 @@ function add_cuts_to_model!(model::SPModel, t::Int64, problem::JuMP.Model, V::Po
     end
 end
 
+
+"""
+Prune all polyhedral functions in input array.
+
+Parameters:
+- model (SPModel)
+- params (SDDPparameters)
+- V Vector{PolyhedralFunction}
+    Polyhedral functions where cuts will be removed
+
+"""
+function prune_cuts!(model::SPModel, params::SDDPparameters, V::Vector{PolyhedralFunction})
+    for i in 1:length(V)
+        V[i] = prune_cuts(model, params, V[i])
+    end
+end
+
+
+"""
+Remove useless cuts in PolyhedralFunction.
+
+Parameters:
+- model (SPModel)
+- params (SDDPparameters)
+- V (PolyhedralFunction)
+    Polyhedral function where cuts will be removed
+
+Return:
+- PolyhedralFunction: pruned polyhedral function
+
+"""
+function prune_cuts(model::SPModel, params::SDDPparameters, V::PolyhedralFunction)
+    ncuts = V.numCuts
+    # Find all active cuts:
+    active_cuts = Bool[is_cut_active(model, i, V, params.solver) for i=1:ncuts]
+    return PolyhedralFunction(V.betas[active_cuts], V.lambdas[active_cuts, :], sum(active_cuts))
+end
+
+
+"""
+Test whether the cut number k is active in polyhedral function Vt.
+
+Parameters:
+- model (SPModel)
+- k (Int)
+    Position of cut to test in PolyhedralFunction object
+- Vt (PolyhedralFunction)
+    Object storing all cuts
+- solver
+    Solver to use to solve linear problem
+
+Return:
+- Bool: true if the cut is active, false otherwise
+
+"""
+function is_cut_active(model::SPModel, k::Int, Vt::PolyhedralFunction, solver)
+
+    m = Model(solver=solver)
+    @defVar(m, alpha)
+    @defVar(m, model.xlim[i][1] <= x[i=1:model.dimStates] <= model.xlim[i][2])
+
+    for i in 1:Vt.numCuts
+        if i!=k
+            lambda = vec(Vt.lambdas[i, :])
+            @addConstraint(m, Vt.betas[i] + dot(lambda, x) <= alpha)
+        end
+    end
+
+    λ_k = vec(Vt.lambdas[k, :])
+    @setObjective(m, Min, alpha - dot(λ_k, x) - Vt.betas[k])
+    solve(m)
+    return getObjectiveValue(m) < 0.
+end
+

src/utils.jl

@@ -78,7 +78,7 @@ Remove redundant cuts in Polyhedral Value functions
 
 """
 function remove_cuts(V::PolyhedralFunction)
-    Vf = hcat(V.lambdas, V.be   tas)
+    Vf = hcat(V.lambdas, V.betas)
     Vf = unique(Vf, 1)
     return PolyhedralFunction(Vf[:, end], Vf[:, 1:end-1], size(Vf)[1])
 end
@@ -89,7 +89,7 @@ Remove redundant cuts in a vector of Polyhedral Functions.
 
 """
 function remove_redundant_cuts!(Vts::Vector{PolyhedralFunction})
-    n_funct ions = length(Vts)
+    n_functions = length(Vts)
     for i in 1:n_functions
         Vts[i] = remove_cuts(Vts[i])
     end

test/runtests.jl

@@ -105,15 +105,12 @@ facts("SDDP algorithm: 1D case") do
     u_bounds, x0,
     cost,
     dynamic, laws)
+    set_state_bounds(model, x_bounds)
     # Generate scenarios for forward simulations:
     noise_scenarios = simulate_scenarios(model.noises,params.forwardPassNumber)
 
     sddp_costs = 0
     context("Linear cost") do
-        # Instantiate a SDDP linear model:
-        set_state_bounds(model, x_bounds)
-
-
         # Compute bellman functions with SDDP:
         V, pbs = solve_SDDP(model, params, 0)
         @fact typeof(V) --> Vector{StochDynamicProgramming.PolyhedralFunction}
@@ -152,6 +149,15 @@ facts("SDDP algorithm: 1D case") do
         @fact mean(sddp_costs) --> roughly(mean(sddp_costs2))
     end
 
+    context("Cuts pruning") do
+        v = V[1]
+        vt = PolyhedralFunction([v.betas[1]; v.betas[1] - 1.], v.lambdas[[1,1],:],  2)
+        StochDynamicProgramming.prune_cuts!(model, params, V)
+        isactive1 = StochDynamicProgramming.is_cut_active(model, 1, vt, params.solver)
+        isactive2 = StochDynamicProgramming.is_cut_active(model, 2, vt, params.solver)
+        @fact isactive1 --> true
+        @fact isactive2 --> false
+    end
 
     context("Piecewise linear cost") do
         # Test Piecewise linear costs: