Merge branch 'territory-heuristic' into cuts-removal

Conflicts:
	src/SDDPoptimize.jl
	src/objects.jl
	test/sddp.jl

Author: frapac

src/SDDPoptimize.jl

@@ -108,7 +108,7 @@ function run_SDDP!(model::SPModel,
 
     #Initialization of the counter
     stats = SDDPStat()
-
+    territory = (param.pruning[:type]=="territory")? [Territories(model.dimStates) for i in 1:model.stageNumber-1]: nothing
     (verbose > 0) && println("Initialize cuts")
 
     # If computation of upper-bound is needed, a set of scenarios is built
@@ -136,7 +136,15 @@ function run_SDDP!(model::SPModel,
 
         ####################
         # cut pruning
-        prune_cuts!(model,param,V,stats.niterations,verbose)
+        prune_cuts!(model, param, V, stockTrajectories, territory, stats.niterations, verbose)
+        ####################
+        # 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 =#
+
+
 
         ####################
         # In iteration upper bound estimation
@@ -149,7 +157,6 @@ function run_SDDP!(model::SPModel,
         updateSDDPStat!(stats, callsolver_forward+callsolver_backward, lwb, upb, toq())
 
         print_current_stats(stats,verbose)
-
         ####################
         # Stopping test
         stopping_test = test_stopping_criterion(param,stats)
@@ -463,3 +470,4 @@ function add_cuts_to_model!(model::SPModel, t::Int64, problem::JuMP.Model, V::Po
         @constraint(problem, V.betas[i] + dot(lambda, xf) <= alpha)
     end
 end
+

src/cutpruning.jl

@@ -1,10 +1,33 @@
+#############################################################################
 #  Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard
 #  This Source Code Form is subject to the terms of the Mozilla Public
 #  License, v. 2.0. If a copy of the MPL was not distributed with this
 #  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 #############################################################################
-# Define the Forward / Backward iterations of the SDDP algorithm
-#############################################################################
+
+type Territories
+    ncuts::Int
+    territories::Array{Array}
+    nstates::Int
+    states::Array{Float64, 2}
+end
+
+
+""" Remove redundant cuts in Polyhedral Value function `V`"""
+function remove_cuts(V::PolyhedralFunction)
+    Vf = hcat(V.lambdas, V.betas)
+    Vf = unique(Vf, 1)
+    return PolyhedralFunction(Vf[:, end], Vf[:, 1:end-1], size(Vf)[1])
+end
+
+
+""" Remove redundant cuts in a vector of Polyhedral Functions `Vts`."""
+function remove_redundant_cuts!(Vts::Vector{PolyhedralFunction})
+    n_functions = length(Vts)-1
+    for i in 1:n_functions
+        Vts[i] = remove_cuts(Vts[i])
+    end
+end
 
 
 """
@@ -15,39 +38,53 @@ Exact pruning of all polyhedral functions in input array.
 * `params::SDDPparameters`:
 * `Vector{PolyhedralFunction}`:
     Polyhedral functions where cuts will be removed
-* `iteration_count::Int64`:
+* `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,
                     param::SDDPparameters,
                     V::Vector{PolyhedralFunction},
-                    iteration_count::Int64,
+                    trajectories::Array{Float64, 3},
+                    territory::Union{Void, Array{Territories}},
+                    it::Int64,
                     verbose::Int64)
-    # If specified, prune cuts:
-    if (param.compute_cuts_pruning > 0) && (iteration_count%param.compute_cuts_pruning==0)
-        (verbose > 0) && println("Prune cuts ...")
-        remove_redundant_cuts!(V)
-        for i in 1:length(V)-1
-            V[i] = exact_prune_cuts(model, param, V[i])
+    # 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
-        problems = hotstart_SDDP(model, param, V)
     end
-end
 
-""" Remove redundant cuts in Polyhedral Value function `V`"""
-function remove_cuts(V::PolyhedralFunction)
-    Vf = hcat(V.lambdas, V.betas)
-    Vf = unique(Vf, 1)
-    return PolyhedralFunction(Vf[:, end], Vf[:, 1:end-1], size(Vf)[1])
-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
 
-""" Remove redundant cuts in a vector of Polyhedral Functions `Vts`."""
-function remove_redundant_cuts!(Vts::Vector{PolyhedralFunction})
-    n_functions = length(Vts)-1
-    for i in 1:n_functions
-        Vts[i] = remove_cuts(Vts[i])
+        # final number of cuts:
+        ncuts_final = get_total_number_cuts(V)
+
+        (verbose > 0) && println(" Deflation: ", ncuts_final/ncuts_initial)
     end
 end
 
@@ -111,3 +148,140 @@ function is_cut_relevant(model::SPModel, k::Int, Vt::PolyhedralFunction, solver;
     sol = getobjectivevalue(m)
     return sol < epsilon
 end
+
+
+########################################
+# Territory algorithm
+########################################
+
+Territories(ndim) = Territories(0, [], 0, Array{Float64}(0, ndim))
+
+
+""" Update territories with cuts previously computed during backward pass.  """
+function find_territory!(territory, V, states)
+    nc = V.numCuts
+    # get number of new positions to analyse:
+    nx = size(states, 1)
+    nt = nc - territory.ncuts
+
+    for i in 1:nt
+        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
+        x = collect(states[i, :])
+        add_state!(territory, V, x)
+    end
+
+end
+
+
+"""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
+        todelete = []
+        for (num, (ix, cost)) in enumerate(territory.territories[k])
+            x = collect(territory.states[ix, :])
+
+            costnewcut = cutvalue(V, indcut, x)
+
+            if costnewcut > cost
+                push!(todelete, num)
+                push!(territory.territories[indcut], (ix, costnewcut))
+            end
+        end
+        deleteat!(territory.territories[k], todelete)
+    end
+end
+
+
+"""Add cut to `territory`."""
+function add_cut!(territory)
+    push!(territory.territories, [])
+    territory.ncuts += 1
+end
+
+
+"""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)
+
+    # 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')
+end
+
+
+"""Remove empty cuts in PolyhedralFunction"""
+function remove_empty_cuts!(territory::Territories, V::PolyhedralFunction)
+    assert(territory.ncuts == V.numCuts)
+
+    nstates = [length(terr) for terr in territory.territories]
+    active_cuts = nstates .> 0
+
+    territory.territories = territory.territories[active_cuts]
+    territory.ncuts = sum(active_cuts)
+    return PolyhedralFunction(V.betas[active_cuts],
+                              V.lambdas[active_cuts, :],
+                              sum(active_cuts))
+end
+
+
+"""Get cut which approximate the best value function at point `x`."""
+function optimalcut(xf::Vector{Float64}, V::PolyhedralFunction)
+    bestcost = -Inf::Float64
+    bestcut = -1
+    nstates = size(V.lambdas, 2)
+    ncuts = size(V.lambdas, 1)
+
+    @inbounds for i in 1:ncuts
+        cost = V.betas[i]
+        for j in 1:nstates
+            cost += V.lambdas[i, j]*xf[j]
+        end
+        if cost > bestcost
+            bestcost = cost
+            bestcut = i
+        end
+    end
+    return bestcost, bestcut
+end
+
+
+"""
+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])
+

src/objects.jl

@@ -19,7 +19,7 @@ type PolyhedralFunction
     numCuts::Int64
 end
 
-PolyhedralFunction(ndim) =  PolyhedralFunction([], Array{Float64}(0, ndim), 0)
+PolyhedralFunction(ndim) = PolyhedralFunction([], Array{Float64}(0, ndim), 0)
 
 
 type LinearSPModel <: SPModel
@@ -152,7 +152,7 @@ type SDDPparameters
     # Maximum iterations of the SDDP algorithms:
     maxItNumber::Int64
     # Prune cuts every %% iterations:
-    compute_cuts_pruning::Int64
+    pruning::Dict{Symbol, Any}
     # Estimate upper-bound every %% iterations:
     compute_ub::Int64
     # Number of MonteCarlo simulation to perform to estimate upper-bound:
@@ -165,15 +165,18 @@ type SDDPparameters
     acceleration::Dict{Symbol, Float64}
 
     function SDDPparameters(solver; passnumber=10, gap=0.,
-                            max_iterations=20, compute_cuts_pruning=0,
+                            max_iterations=20, prune_cuts=0,
+                            pruning_algo="exact",
                             compute_ub=-1, montecarlo_final=10000, montecarlo_in_iter = 100,
                             mipsolver=nothing,
                             rho0=0., alpha=1.)
         is_acc = (rho0 > 0.)
         accparams = is_acc? Dict(:ρ0=>rho0, :alpha=>alpha, :rho=>rho0): Dict()
 
+        prune_cuts = Dict(:pruning=>prune_cuts>0, :period=>prune_cuts, :type=>pruning_algo)
         return new(solver, mipsolver, passnumber, gap,
-                   max_iterations, compute_cuts_pruning, compute_ub, montecarlo_final,montecarlo_in_iter, is_acc, accparams)
+                   max_iterations, prune_cuts, compute_ub,
+                   montecarlo_final, montecarlo_in_iter, is_acc, accparams)
     end
 end
 

test/sddp.jl

@@ -46,7 +46,7 @@ facts("SDDP algorithm: 1D case") do
                                                     passnumber=n_scenarios,
                                                     gap=epsilon,
                                                     max_iterations=max_iterations,
-                                                    compute_cuts_pruning=1)
+                                                    prune_cuts=0)
 
     V = nothing
     model = StochDynamicProgramming.LinearSPModel(n_stages, u_bounds,
@@ -71,8 +71,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
@@ -109,7 +109,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, param, V,1,0)
         isactive1 = StochDynamicProgramming.is_cut_relevant(model, 1, vt, param.SOLVER)
         isactive2 = StochDynamicProgramming.is_cut_relevant(model, 2, vt, param.SOLVER)
         @fact isactive1 --> true
@@ -163,7 +162,7 @@ facts("SDDP algorithm: 1D case") do
     context("Stopping criterion") do
         # Compute upper bound every %% iterations:
         param.compute_ub = 1
-        param.compute_cuts_pruning = 1
+        param.pruning = Dict(:pruning=>true, :period=>1, :type=>"exact")
         param.maxItNumber = 30
         V, pbs = solve_SDDP(model, param, V, 0)
         V0 = StochDynamicProgramming.get_lower_bound(model, param, V)