Merge pull request #120 from JuliaOpt/naive-parallel-sddp

Naive parallelization

Author: leclere

examples/parallel_sddp.jl

@@ -0,0 +1,58 @@
+
+import StochDynamicProgramming
+
+
+"""
+Solve SDDP in parallel, dispatching both forward and backward passes to process, 
+which is not the most standard parallelization of SDDP.
+
+# Arguments
+* `model::SPmodel`:
+    the stochastic problem we want to optimize
+* `param::SDDPparameters`:
+    the parameters of the SDDP algorithm
+* `V::Array{PolyhedralFunction}`:
+    the current estimation of Bellman's functions
+* `n_parallel_pass::Int`: default is 4
+    Number of parallel pass to compute
+* `synchronize::Int`: default is 5
+    Synchronize the cuts between the different processes every "synchronise" iterations
+* `display::Int`: default is 0
+    Says whether to display results or not
+
+# Return
+* `V::Array{PolyhedralFunction}`:
+    the collection of approximation of the bellman functions
+"""
+function psolve_sddp(model, params, V; n_parallel_pass=4,
+                     synchronize=5, display=0)
+    # Redefine seeds in every processes to maximize randomness:
+    @everywhere srand()
+
+    mitn = params.maxItNumber
+    params.maxItNumber = synchronize
+
+    # Count number of available CPU:
+    ncpu = nprocs() - 1
+    (display > 0) && println("\nLaunch simulation on ", ncpu, " processes")
+    workers = procs()[2:end]
+
+    fpn = params.forwardPassNumber
+    # As we distribute computation in n process, we perform forward pass in parallel:
+    params.forwardPassNumber = max(1, round(Int, params.forwardPassNumber/ncpu))
+
+    # Start parallel computation:
+    for i in 1:n_parallel_pass
+        # Distribute computation of SDDP in each process:
+        refs = [@spawnat w StochDynamicProgramming.solve_SDDP(model, params, V, display)[1] for w in workers]
+        # Catch the result in the main process:
+        V = StochDynamicProgramming.catcutsarray([fetch(r) for r in refs]...)
+        # We clean the resultant cuts:
+        StochDynamicProgramming.remove_redundant_cuts!(V)
+        (display > 0) && println("Lower bound at pass ", i, ": ", StochDynamicProgramming.get_lower_bound(model, params, V))
+    end
+
+    params.forwardPassNumber = fpn
+    params.maxItNumber = mitn
+    return V
+end

src/SDDPoptimize.jl

@@ -311,6 +311,8 @@ function initialize_value_functions(model::SPModel,
         V[end] = model.finalCost
         build_terminal_cost!(model, solverProblems[end], V[end])
     elseif isa(model.finalCost, Function)
+        # In this case, define a trivial value functions for final cost to avoid problem:
+        V[end] = PolyhedralFunction(zeros(1), zeros(1, model.dimStates), 1)
         model.finalCost(model, solverProblems[end])
     end
 

src/utils.jl

@@ -72,6 +72,24 @@ function read_polyhedral_functions(dump::AbstractString)
 end
 
 
+"""Concatenate collection of arrays of PolyhedralFunction."""
+function catcutsarray(polyfunarray::Vector{StochDynamicProgramming.PolyhedralFunction}...)
+    assert(length(polyfunarray) > 0)
+    ntimes = length(polyfunarray[1])
+    # Concatenate cuts in polyfunarray, and discard final time as we do not add cuts at final time:
+    concatcuts = StochDynamicProgramming.PolyhedralFunction[catcuts([V[t] for V in polyfunarray]...) for t in 1:ntimes-1]
+    return vcat(concatcuts, polyfunarray[1][end])
+end
+
+
+"""Concatenate collection of PolyhedralFunction."""
+function catcuts(Vts::StochDynamicProgramming.PolyhedralFunction...)
+    betas = vcat([V.betas for V in Vts]...)
+    lambdas = vcat([V.lambdas for V in Vts]...)
+    numcuts = sum([V.numCuts for V in Vts])
+    return StochDynamicProgramming.PolyhedralFunction(betas, lambdas, numcuts)
+end
+
 """
 Extract a vector stored in a 3D Array