[DEV] Integrate territory to SDDP interface

Author: frapac

src/SDDPoptimize.jl

@@ -69,7 +69,7 @@ function run_SDDP!(model::SPModel,
 
     #Initialization of the counter
     stats = SDDPStat(0, [], [], [], 0)
-    territory = [Territories(model.dimStates) for i in 1:model.stageNumber-1]
+    territory = (param.pruning[:type]=="territory")? [Territories(model.dimStates) for i in 1:model.stageNumber-1]: nothing
 
     (verbose > 0) && println("Initialize cuts")
 
@@ -114,24 +114,29 @@ function run_SDDP!(model::SPModel,
         iteration_count += 1
         stats.niterations += 1
 
-        ####################
-        # If specified, prune cuts:
-        if (param.pruning[:period] > 0) && (iteration_count%param.pruning[:period]==0)
-            (verbose > 0) && println("Prune cuts ...")
-            remove_redundant_cuts!(V)
-            prune_cuts!(model, param, V)
-            problems = hotstart_SDDP(model, param, V)
-        end
 
         # Update estimation of lower bound:
         lwb = get_bellman_value(model, param, 1, V[1], model.initialState)
         push!(stats.lower_bounds, lwb)
+        push!(stats.exectime, toq())
+        push!(stats.upper_bounds, upb)
 
-        for t in 1:model.stageNumber-1
-            states = reshape(stockTrajectories[t, :, :], param.forwardPassNumber, model.dimStates)
-            find_territory!(territory[t], V[t], states)
+        if (verbose > 0) && (iteration_count%verbose==0)
+            print("Pass number ", iteration_count)
+            (upb < Inf) && print("\tUpper-bound: ", upb)
+            println("\tLower-bound: ", round(stats.lower_bounds[end], 4),
+                "\tTime: ", round(stats.exectime[end], 2),"s")
         end
 
+        ####################
+        # Cut pruning
+        prune_cuts!(model, param, V, stockTrajectories, territory, iteration_count, verbose)
+        if (param.pruning[:period] > 0) && (iteration_count%param.pruning[:period]==0)
+            problems = hotstart_SDDP(model, param, V)
+        end
+
+
+
         ####################
         # If specified, compute upper-bound:
         if (param.compute_upper_bound > 0) && (iteration_count%param.compute_upper_bound==0)
@@ -144,18 +149,8 @@ function run_SDDP!(model::SPModel,
             end
         end
 
-        push!(stats.exectime, toq())
-        push!(stats.upper_bounds, upb)
-
-        if (verbose > 0) && (iteration_count%verbose==0)
-            print("Pass number ", iteration_count)
-            (upb < Inf) && print("\tUpper-bound: ", upb)
-            println("\tLower-bound: ", round(stats.lower_bounds[end], 4),
-                "\tTime: ", round(stats.exectime[end], 2),"s")
-        end
 
     end
-    println([length(v) for v in territory[9].territories])
 
     ##########
     # Estimate final upper bound with param.monteCarloSize simulations:

src/objects.jl

@@ -154,7 +154,7 @@ type SDDPparameters
                             max_iterations=20, prune_cuts=0,
                             pruning_algo="exact",
                             compute_ub=-1, montecarlo=10000, mipsolver=nothing)
-        prune_cuts = Dict(:period=>prune_cuts, :type=>pruning_algo)
+        prune_cuts = Dict(:pruning=>prune_cuts>0, :period=>prune_cuts, :type=>pruning_algo)
         return new(solver, mipsolver, passnumber, gap, max_iterations, prune_cuts, compute_ub, montecarlo)
     end
 end

src/pruning.jl

@@ -5,16 +5,13 @@
 #  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 #############################################################################
 
-
 type Territories
     ncuts::Int
     territories::Array{Array}
     nstates::Int
     states::Array{Float64, 2}
 end
 
-Territories(ndim) = Territories(0, [], 0, Array{Float64}(0, ndim))
-
 
 """
 Exact pruning of all polyhedral functions in input array.
@@ -24,10 +21,53 @@ Exact pruning of all polyhedral functions in input array.
 * `params::SDDPparameters`:
 * `Vector{PolyhedralFunction}`:
     Polyhedral functions where cuts will be removed
+* `trajectories::Array{Float64, 3}`
+    Previous trajectories
+* `territory::Array{Territories}`
+    Container storing the territory for each cuts
+* `it::Int64`:
+    current iteration number
+* `verbose::Int64`
 """
-function prune_cuts!(model::SPModel, params::SDDPparameters, V::Vector{PolyhedralFunction})
-    for i in 1:length(V)-1
-        V[i] = exact_prune_cuts(model, params, V[i])
+function prune_cuts!(model::SPModel,
+                    param::SDDPparameters,
+                    V::Vector{PolyhedralFunction},
+                    trajectories::Array{Float64, 3},
+                    territory::Union{Void, Array{Territories}},
+                    it::Int64,
+                    verbose::Int64)
+    # Basic pruning: remove redundant cuts
+    remove_redundant_cuts!(V)
+
+    # If pruning is performed with territory heuristic, update territory
+    # at given iteration:
+    if param.pruning[:type] == "territory"
+        for t in 1:model.stageNumber-1
+            states = reshape(trajectories[t, :, :], param.forwardPassNumber, model.dimStates)
+            find_territory!(territory[t], V[t], states)
+        end
+    end
+
+    # If specified to prune cuts at this iteration, do it:
+    if param.pruning[:pruning] && (it%param.pruning[:period]==0)
+        # initial number of cuts:
+        ncuts_initial = get_total_number_cuts(V)
+        (verbose > 0) && print("Prune cuts ...")
+
+        for i in 1:length(V)-1
+            if param.pruning[:type] == "exact"
+                # apply exact cuts pruning:
+                V[i] = exact_prune_cuts(model, param, V[i])
+            elseif param.pruning[:type] == "territory"
+                # apply heuristic to prune cuts:
+                V[i] = remove_empty_cuts!(territory[i], V[i])
+            end
+        end
+
+        # final number of cuts:
+        ncuts_final = get_total_number_cuts(V)
+
+        (verbose > 0) && println(" Deflation: ", ncuts_final/ncuts_initial)
     end
 end
 
@@ -96,12 +136,8 @@ end
 ########################################
 # Territory algorithm
 ########################################
-function territory_prune_cuts!(territory::Territories,
-                               V::PolyhedralFunction,
-                               states::Array{Float64, 2})
-    find_territory!(territory, V, states)
-    return remove_empty_cuts(territory, V)
-end
+
+Territories(ndim) = Territories(0, [], 0, Array{Float64}(0, ndim))
 
 
 """ Update territories with cuts previously computed during backward pass.  """
@@ -112,9 +148,11 @@ function find_territory!(territory, V, states)
     nt = nc - territory.ncuts
 
     for i in 1:nt
-        add_cut!(territory, V)
+        add_cut!(territory)
         update_territory!(territory, V, nc - nt + i)
     end
+
+    # ensure that territory has the same number of cuts as V!
     assert(territory.ncuts == V.numCuts)
 
     for i in 1:nx
@@ -125,15 +163,14 @@ function find_territory!(territory, V, states)
 end
 
 
-""" Update territories with new cuts added during previous backward pass.  """
+"""Update territories considering new cut with index `indcut`."""
 function update_territory!(territory, V, indcut)
     for k in 1:territory.ncuts
         if k == indcut
             continue
         end
-        terr = copy(territory.territories[k])
         todelete = []
-        for (num, (ix, cost)) in enumerate(terr)
+        for (num, (ix, cost)) in enumerate(territory.territories[k])
             x = collect(territory.states[ix, :])
 
             costnewcut = cutvalue(V, indcut, x)
@@ -148,31 +185,32 @@ function update_territory!(territory, V, indcut)
 end
 
 
-"""Add a new cut to the"""
-function add_cut!(territory, V)
+"""Add cut to `territory`."""
+function add_cut!(territory)
     push!(territory.territories, [])
     territory.ncuts += 1
 end
 
 
-function add_state!(territory, V, x)
+"""Add a new state and update territories."""
+function add_state!(territory::Territories, V::PolyhedralFunction, x::Array{Float64})
+    # Get cut which is the supremum at point `x`:
     bcost, bcuts = optimalcut(x, V)
 
-    indx = territory.nstates + 1
-    push!(territory.territories[bcuts], (indx, bcost))
+    # Add `x` to the territory of cut `bcuts`:
+    territory.nstates += 1
+    push!(territory.territories[bcuts], (territory.nstates, bcost))
 
+    # Add `x` to the list of visited state:
     territory.states = vcat(territory.states, x')
-    territory.nstates += 1
 end
 
 
-"""
-Remove empty cuts in PolyhedralFunction
-"""
-function remove_empty_cuts!(territory, V)
-    assert(length(territory) == V.numCuts)
+"""Remove empty cuts in PolyhedralFunction"""
+function remove_empty_cuts!(territory::Territories, V::PolyhedralFunction)
+    assert(territory.ncuts == V.numCuts)
 
-    nstates = [length(cont) for cont in territory.territories]
+    nstates = [length(terr) for terr in territory.territories]
     active_cuts = nstates .> 0
 
     territory.territories = territory.territories[active_cuts]
@@ -183,9 +221,7 @@ function remove_empty_cuts!(territory, V)
 end
 
 
-"""
-Get cut which approximate the best value function at point `x`.
-"""
+"""Get cut which approximate the best value function at point `x`."""
 function optimalcut(xf::Vector{Float64}, V::PolyhedralFunction)
     bestcost = -Inf::Float64
     bestcut = -1
@@ -206,12 +242,29 @@ function optimalcut(xf::Vector{Float64}, V::PolyhedralFunction)
 end
 
 
-""" Get approximation of value function at given point `x`.  """
-function cutvalue(V, indc, x)
+"""
+Get approximation of value function at given point `x`.
+
+# Arguments
+- `V::PolyhedralFunction`
+    Approximation of the value function as linear cuts
+- `indc::Int64`
+    Index of cut to consider
+- `x::Array{Float64}`
+    Coordinates of state
+
+# Return
+`cost::Float64`
+    Value of cut `indc` at point `x`
+"""
+function cutvalue(V::PolyhedralFunction, indc::Int, x::Array{Float64})
     cost = V.betas[indc]
     for j in 1:size(V.lambdas, 2)
         cost += V.lambdas[indc, j]*x[j]
     end
     return cost
 end
 
+# Return total number of cuts in PolyhedralFunction array:
+get_total_number_cuts(V::Array{PolyhedralFunction}) = sum([v.numCuts for v in V])
+

test/sddp.jl

@@ -70,8 +70,8 @@ facts("SDDP algorithm: 1D case") do
 
         # Test if the first subgradient has the same dimension as state:
         @fact size(V[1].lambdas, 2) --> model.dimStates
-        @fact V[1].numCuts --> n_scenarios*max_iterations + n_scenarios
-        @fact size(V[1].lambdas, 1) --> n_scenarios*max_iterations + n_scenarios
+        @fact V[1].numCuts < n_scenarios*max_iterations + n_scenarios --> true
+        @fact size(V[1].lambdas, 1) --> V[1].numCuts
 
         # Test upper bounds estimation with Monte-Carlo:
         n_simulations = 100
@@ -108,7 +108,6 @@ facts("SDDP algorithm: 1D case") do
     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_relevant(model, 1, vt, params.SOLVER)
         isactive2 = StochDynamicProgramming.is_cut_relevant(model, 2, vt, params.SOLVER)
         @fact isactive1 --> true
@@ -153,7 +152,7 @@ facts("SDDP algorithm: 1D case") do
     context("Stopping criterion") do
         # Compute upper bound every %% iterations:
         params.compute_upper_bound = 1
-        params.pruning = Dict(:period=>1, :type=>"exact")
+        params.pruning = Dict(:pruning=>true, :period=>1, :type=>"exact")
         params.maxItNumber = 30
         V, pbs = solve_SDDP(model, params, V, 0)
         V0 = StochDynamicProgramming.get_lower_bound(model, params, V)