Merge pull request #45 from leclere/correct-forwardpassnumber-use

[FIX] fix use of param.forwardPassNumber and extract_vector_from_3Dmatrix

Author: leclere

src/SDDPoptimize.jl

@@ -92,7 +92,7 @@ function run_SDDP(model::SPModel,
 
         # Build given number of scenarios according to distribution
         # law specified in model.noises:
-        aleas = simulate_scenarios(model.noises ,
+        noise_scenarios = simulate_scenarios(model.noises ,
                                     (model.stageNumber-1,
                                      param.forwardPassNumber,
                                      model.dimNoises))
@@ -102,7 +102,7 @@ function run_SDDP(model::SPModel,
                             param,
                             V,
                             problems,
-                            aleas)
+                            noise_scenarios)
 
         # Backward pass
         backward_pass!(model,
@@ -115,23 +115,25 @@ function run_SDDP(model::SPModel,
                       returnValueFunctions)
 
         iteration_count += 1
-        upb = upper_bound(costs)
 
-        V0 = get_bellman_value(model, param, 1, V[1], model.initialState)
 
-        time = toq()
+        
 
         if (display > 0) && (iteration_count%display==0)
             println("Pass number ", iteration_count,
-                    "\tEstimation of upper-bound: ", upb,
-                    "\tLower-bound: ", V0,
-                    "\tTime: ", time)
+                    "\tEstimation of upper-bound: ", upper_bound(costs),
+                    "\tLower-bound: ", get_bellman_value(model, param, 1, V[1], model.initialState),
+                    "\tTime: ", toq())
         end
 
     end
 
     # Estimate upper bound with a great number of simulations:
     if (display>0)
+        upb = upper_bound(costs)
+        V0 = get_bellman_value(model, param, 1, V[1], model.initialState)
+        time = toq()
+        
         println("Estimate upper-bound with Monte-Carlo ...")
         upb, costs = estimate_upper_bound(model, param, V, problems)
         println("Estimation of upper-bound: ", upb,
@@ -172,12 +174,9 @@ Float64 (estimation of the upper bound)
 """
 function estimate_upper_bound(model, param, V, problems, n_simulation=1000)
 
-    n_fpn = param.forwardPassNumber
-    param.forwardPassNumber = n_simulation
-
     aleas = simulate_scenarios(model.noises ,
                                     (model.stageNumber-1,
-                                     param.forwardPassNumber,
+                                     n_simulation,
                                      model.dimNoises))
 
     costs, stockTrajectories, _ = forward_simulations(model,
@@ -186,9 +185,6 @@ function estimate_upper_bound(model, param, V, problems, n_simulation=1000)
                                                         problems,
                                                         aleas)
 
-
-    param.forwardPassNumber = n_fpn
-
     return upper_bound(costs), costs
 end
 
@@ -463,4 +459,4 @@ function add_cuts_to_model!(model::SPModel, t::Int64, problem::JuMP.Model, V::Po
         lambda = vec(V.lambdas[i, :])
         @addConstraint(problem, V.betas[i] + dot(lambda, model.dynamics(t, x, u, w)) <= alpha)
     end
-end
+end

src/forwardBackwardIterations.jl

@@ -27,8 +27,8 @@ Parameters:
     Linear model used to approximate each value function
 
 - xi (Array{float})
-    the noise scenarios on which we simulate, each line being one scenario.
-    Generated if not given.
+    the noise scenarios on which we simulate, each column being one scenario :
+    xi[t,k,:] is the alea at time t of scenario k.
 
 - returnCosts (Bool)
     return the cost of each simulated scenario if true
@@ -47,12 +47,12 @@ Returns (according to the last parameters):
     If returnCosts=false, return nothing
 
 - stocks (Array{float})
-    the simulated stock trajectories. stocks(k,t,:) is the stock for
+    the simulated stock trajectories. stocks(t,k,:) is the stock for
     scenario k at time t.
 
 - controls (Array{Float64, 3})
-
-
+    the simulated controls trajectories. controls(t,k,:) is the control for
+    scenario k at time t.
 """
 function forward_simulations(model::SPModel,
                             param::SDDPparameters,
@@ -65,26 +65,28 @@ function forward_simulations(model::SPModel,
 
     # TODO: verify that loops are in the same order
     T = model.stageNumber
-    stocks = zeros(T, param.forwardPassNumber, model.dimStates)
+    nb_forward = size(xi)[2] 
+    
+    stocks = zeros(T, nb_forward, model.dimStates)
     # We got T - 1 control, as terminal state is included into the total number
     # of stages.
-    controls = zeros(T - 1, param.forwardPassNumber, model.dimControls)
+    controls = zeros(T - 1, nb_forward, model.dimControls)
 
     # Set first value of stocks equal to x0:
-    for i in 1:param.forwardPassNumber
-        stocks[1, i, :] = model.initialState
+    for k in 1:nb_forward
+        stocks[1, k, :] = model.initialState
     end
 
     costs = nothing
     if returnCosts
-        costs = zeros(param.forwardPassNumber)
+        costs = zeros(nb_forward)
     end
 
     for t=1:T-1
-        for k = 1:param.forwardPassNumber
+        for k = 1:nb_forward
 
-            state_t = extract_vector_from_3Dmatrix(stocks, k, t)
-            alea_t = extract_vector_from_3Dmatrix(xi, k, t)
+            state_t = extract_vector_from_3Dmatrix(stocks, t, k)
+            alea_t = extract_vector_from_3Dmatrix(xi, t, k)
 
             status, nextstep = solve_one_step_one_alea(
                                             model,
@@ -190,7 +192,7 @@ Parameters:
     Linear model used to approximate each value function
 
 - stockTrajectories (Array{Float64,3})
-    stockTrajectories[k,t,:] is the vector of stock where the cut is computed
+    stockTrajectories[t,k,:] is the vector of stock where the cut is computed
     for scenario k and time t.
 
 - law (Array{NoiseLaw})
@@ -218,21 +220,22 @@ function backward_pass!(model::SPModel,
                       updateV=false)
 
     T = model.stageNumber
-
+    nb_forward = size(stockTrajectories)[2]
+    
     # Estimation of initial cost:
     V0 = 0.
 
     costs::Vector{Float64} = zeros(1)
-    costs_npass = zeros(Float64, param.forwardPassNumber)
+    costs_npass = zeros(Float64, nb_forward)
     state_t = zeros(Float64, model.dimStates)
 
     for t = T-1:-1:1
         costs = zeros(law[t].supportSize)
 
-        for k = 1:param.forwardPassNumber
+        for k = 1:nb_forward
 
             subgradient_array = zeros(Float64, model.dimStates, law[t].supportSize)
-            state_t = extract_vector_from_3Dmatrix(stockTrajectories, k, t)
+            state_t = extract_vector_from_3Dmatrix(stockTrajectories, t, k)
 
             for w in 1:law[t].supportSize
 

src/utils.jl

@@ -96,7 +96,7 @@ function extract_vector_from_3Dmatrix(input_array::Array{Float64, 3},
                                       ny::Int64)
 
     state_dimension = size(input_array)[3]
-    return reshape(input_array[ny, nx, :], state_dimension)
+    return reshape(input_array[nx, ny, :], state_dimension)
 end
 
 

test/runtests.jl

@@ -45,9 +45,10 @@ end
 facts("Utils functions") do
     # Test extraction of vector in array:
     arr = rand(4, 4, 2)
-    vec = StochDynamicProgramming.extract_vector_from_3Dmatrix(arr, 2, 1)
-    @fact typeof(vec) --> Vector{Float64}
-    @fact size(vec) --> (2,)
+    v = StochDynamicProgramming.extract_vector_from_3Dmatrix(arr, 2, 1)
+    @fact typeof(v) --> Vector{Float64}
+    @fact size(v) --> (2,)
+    @fact v --> vec(arr[2, 1,:])
 
     # Test upper bound calculation:
     cost = rand(10)
@@ -108,7 +109,7 @@ facts("SDDP algorithm: 1D case") do
                                                 cost,
                                                 dynamic, laws)
     # Generate scenarios for forward simulations:
-    aleas = simulate_scenarios(model.noises,
+    noise_scenarios = simulate_scenarios(model.noises,
                           (model.stageNumber,
                            params.forwardPassNumber,
                            model.dimNoises))
@@ -135,7 +136,7 @@ facts("SDDP algorithm: 1D case") do
                                                            n_simulations)[1]
         @fact typeof(upb) --> Float64
 
-        sddp_costs, stocks = forward_simulations(model, params, V, pbs, aleas)
+        sddp_costs, stocks = forward_simulations(model, params, V, pbs, noise_scenarios)
 
         # Compare sddp cost with those given by extensive formulation:
         ef_cost = StochDynamicProgramming.extensive_formulation(model,params)
@@ -148,7 +149,7 @@ facts("SDDP algorithm: 1D case") do
         # Test hot start with previously computed value functions:
         V, pbs = solve_SDDP(model, params, 0, V)
         # Test if costs are roughly the same:
-        sddp_costs2, stocks = forward_simulations(model, params, V, pbs, aleas)
+        sddp_costs2, stocks = forward_simulations(model, params, V, pbs, noise_scenarios)
         @fact mean(sddp_costs) --> roughly(mean(sddp_costs2))
     end
 
@@ -244,13 +245,12 @@ facts("SDDP algorithm: 2D case") do
 
 
          # Test a simulation upon given scenarios:
-         params.forwardPassNumber = n_simulations
-        aleas = simulate_scenarios(model.noises,
+        noise_scenarios = simulate_scenarios(model.noises,
                               (model.stageNumber,
-                               params.forwardPassNumber,
+                               n_simulations,
                                model.dimNoises))
 
-        sddp_costs, stocks = forward_simulations(model, params, V, pbs, aleas)
+        sddp_costs, stocks = forward_simulations(model, params, V, pbs, noise_scenarios)
 
         # Compare sddp cost with those given by extensive formulation:
         ef_cost = StochDynamicProgramming.extensive_formulation(model,params)
@@ -495,4 +495,4 @@ facts("SDP algorithm") do
         @fact ((state_ref[1]<=state_neighbor[1])==(value_bar_ref[1]<=value_bar_neighbor[1])) --> true
     end
 
-end
+end