[UPD] Improve initialization of value functions

Author: frapac

src/SDDPoptimize.jl

@@ -98,8 +98,7 @@ function run_SDDP!(model::SPModel,
                       V,
                       problems,
                       stockTrajectories,
-                      model.noises,
-                      false)
+                      model.noises)
 
         #Update the number of call
         stats.ncallsolver += callsolver_forward + callsolver_backward
@@ -315,7 +314,8 @@ function initialize_value_functions(model::SPModel,
 
     solverProblems = build_models(model, param)
 
-    V = Array{PolyhedralFunction}(model.stageNumber)
+    V = PolyhedralFunction[
+                PolyhedralFunction([], Array{Float64}(0, model.dimStates), 0) for i in 1:model.stageNumber]
 
     # Build scenarios according to distribution laws:
     aleas = simulate_scenarios(model.noises, param.forwardPassNumber)
@@ -339,8 +339,7 @@ function initialize_value_functions(model::SPModel,
                   V,
                   solverProblems,
                   stockTrajectories,
-                  model.noises,
-                  true)
+                  model.noises)
 
     return V, solverProblems
 end

src/forwardBackwardIterations.jl

@@ -169,16 +169,13 @@ the current estimation of Vt.
     for scenario k and time t.
 * `law::Array{NoiseLaw}`:
     Conditionnal distributions of perturbation, for each timestep
-* `init::Bool`:
-    If specified, then init PolyhedralFunction
 """
 function backward_pass!(model::SPModel,
                         param::SDDPparameters,
                         V::Vector{PolyhedralFunction},
                         solverProblems::Vector{JuMP.Model},
                         stockTrajectories::Array{Float64, 3},
-                        law,
-                        init=false::Bool)
+                        law)
 
     callsolver::Int = 0
 
@@ -222,11 +219,7 @@ function backward_pass!(model::SPModel,
                 beta = costs_npass - dot(subgradient, state_t)
 
                 # Add cut to polyhedral function and JuMP model:
-                if init
-                    V[t] = PolyhedralFunction([beta], reshape(subgradient, 1, model.dimStates), 1)
-                else
-                    add_cut!(model, t, V[t], beta, subgradient)
-                end
+                add_cut!(model, t, V[t], beta, subgradient)
                 if t > 1
                     add_cut_to_model!(model, solverProblems[t-1], t, beta, subgradient)
                 end

test/sddp.jl

@@ -75,9 +75,9 @@ facts("SDDP algorithm: 1D case") do
         @fact length(V) --> n_stages
 
         # Test if the first subgradient has the same dimension as state:
-        @fact length(V[1].lambdas[1, :]) --> model.dimStates
-        @fact V[1].numCuts --> n_scenarios*max_iterations + 1
-        @fact length(V[1].lambdas[:, 1]) --> n_scenarios*max_iterations + 1
+        @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
 
         # Test upper bounds estimation with Monte-Carlo:
         n_simulations = 100