Merge pull request #99 from JuliaOpt/refactoring

Refactoring

Author: leclere

src/SDDPoptimize.jl

@@ -34,22 +34,22 @@ fulfilled.
     each value function
 * `count_callsolver::Int64`:
     number of times the solver has been called
-    
+
 """
 function solve_SDDP(model::SPModel, param::SDDPparameters, display=0::Int64)
     # initialize value functions:
     V, problems = initialize_value_functions(model, param)
     # Run SDDP upon example:
-    count_callsolver = run_SDDP!(model, param, V, problems, display)
-    return V, problems, count_callsolver
+    sddp_stats = run_SDDP!(model, param, V, problems, display)
+    return V, problems, sddp_stats
 end
 
 
 function solve_SDDP(model::SPModel, param::SDDPparameters, V::Vector{PolyhedralFunction}, display=0::Int64)
     # First step: process value functions if hotstart is called
     problems = hotstart_SDDP(model, param, V)
-    count_callsolver = run_SDDP!(model, param, V, problems, display)
-    return V, problems, count_callsolver
+    sddp_stats = run_SDDP!(model, param, V, problems, display)
+    return V, problems, sddp_stats
 end
 
 
@@ -61,7 +61,7 @@ function run_SDDP!(model::SPModel,
                     display=0::Int64)
 
     #Initialization of the counter
-    count_callsolver::Int = 0
+    stats = SDDPStat(0, [], [], [], 0)
 
     if display > 0
       println("Initialize cuts")
@@ -77,10 +77,9 @@ function run_SDDP!(model::SPModel,
     iteration_count::Int64 = 0
 
     while (iteration_count < param.maxItNumber) & (~stopping_test)
+
         # Time execution of current pass:
-        if display > 0
-            tic()
-        end
+        tic()
 
         # Build given number of scenarios according to distribution
         # law specified in model.noises:
@@ -103,32 +102,37 @@ function run_SDDP!(model::SPModel,
                       false)
 
         #Update the number of call
-        count_callsolver += callsolver_forward + callsolver_backward
+        stats.ncallsolver += callsolver_forward + callsolver_backward
 
         iteration_count += 1
+        stats.niterations += 1
 
         if (param.compute_cuts_pruning > 0) && (iteration_count%param.compute_cuts_pruning==0)
             (display > 0) && println("Prune cuts ...")
             remove_redundant_cuts!(V)
             prune_cuts!(model, param, V)
             problems = hotstart_SDDP(model, param, V)
         end
+        lwb = get_bellman_value(model, param, 1, V[1], model.initialState)
+        push!(stats.lower_bounds, lwb)
 
         if (param.compute_upper_bound > 0) && (iteration_count%param.compute_upper_bound==0)
             (display > 0) && println("Compute upper-bound with ",
                                       param.monteCarloSize, " scenarios...")
             upb, costs = estimate_upper_bound(model, param, upperbound_scenarios, problems)
             if param.gap > 0.
-                lwb = get_bellman_value(model, param, 1, V[1], model.initialState)
                 stopping_test = test_stopping_criterion(lwb, upb, param.gap)
             end
         end
 
+        push!(stats.exectime, toq())
+        push!(stats.upper_bounds, upb)
+
         if (display > 0) && (iteration_count%display==0)
             println("Pass number ", iteration_count,
                     "\tUpper-bound: ", upb,
-                    "\tLower-bound: ", round(get_bellman_value(model, param, 1, V[1], model.initialState),4),
-                    "\tTime: ", round(toq(),2),"s")
+                    "\tLower-bound: ", round(stats.lower_bounds[end], 4),
+                    "\tTime: ", round(stats.exectime[end], 2),"s")
         end
 
     end
@@ -149,7 +153,7 @@ function run_SDDP!(model::SPModel,
                  round(mean(costs),4), " +/- ", round(1.96*std(costs)/sqrt(length(costs)),4))
     end
 
-    return count_callsolver
+    return stats
 end
 
 
@@ -180,15 +184,11 @@ function estimate_upper_bound(model::SPModel, param::SDDPparameters,
                                 n_simulation=1000::Int)
 
     aleas = simulate_scenarios(model.noises, n_simulation)
-
-    callsolver::Int = 0
-    costs, stockTrajectories, _ = forward_simulations(model,
-                                                        param,
-                                                        problem,
-                                                        aleas)
-
+    costs, stockTrajectories, _ = forward_simulations(model, param, problem, aleas)
     return upper_bound(costs), costs
 end
+
+
 function estimate_upper_bound(model::SPModel, param::SDDPparameters,
                                 aleas::Array{Float64, 3},
                                 problem::Vector{JuMP.Model})
@@ -199,13 +199,7 @@ end
 
 """Build a collection of cuts initialized at 0"""
 function get_null_value_functions_array(model::SPModel)
-
-    V = Vector{PolyhedralFunction}(model.stageNumber)
-    for t = 1:model.stageNumber
-        V[t] = PolyhedralFunction(zeros(1), zeros(1, model.dimStates), 1)
-    end
-
-    return V
+    return [PolyhedralFunction(zeros(1), zeros(1, model.dimStates), 1) for i in 1:model.stageNumber]
 end
 
 
@@ -298,7 +292,6 @@ function build_models(model::SPModel, param::SDDPparameters)
 end
 
 
-
 """
 Initialize value functions along a given trajectory
 
@@ -461,7 +454,8 @@ Compute optimal control at point xt and time t.
 # Return
     `Vector{Float64}`: optimal control at time t
 """
-function get_control(model::SPModel, param::SDDPparameters, lpproblem::Vector{JuMP.Model}, t::Int, xt::Vector{Float64}, xi::Vector{Float64})
+function get_control(model::SPModel, param::SDDPparameters, lpproblem::Vector{JuMP.Model},
+                     t::Int, xt::Vector{Float64}, xi::Vector{Float64})
     return solve_one_step_one_alea(model, param, lpproblem[t], t, xt, xi)[2].optimal_control
 end
 
@@ -565,3 +559,4 @@ function is_cut_relevant(model::SPModel, k::Int, Vt::PolyhedralFunction, solver;
     sol = getobjectivevalue(m)
     return sol < epsilon
 end
+

src/compare.jl

@@ -12,7 +12,7 @@ Compare different sets of parameters to solve an instance of SDDP
 
 # Description
 Take a collection of SDDP parameters and compare the time of execution,
-the memory used, an estimation of the gap to optimality and the number 
+the memory used, an estimation of the gap to optimality and the number
 of calls to the solver
 
 # Arguments
@@ -24,7 +24,7 @@ of calls to the solver
     Set of scenarios used to calculate costs
 * `seeds::Int`
     The random number generator seeds
-    
+
 # Output
 * `Display in the terminal`
     Print information in the terminal
@@ -44,21 +44,21 @@ function benchmark_parameters(model,
 
         srand(seeds)
 
-        (V, pbs, callsolver), t1, m1 = @timed solve_SDDP(model, sddpparams, 0)
+        (V, pbs, sddpstats), t1, m1 = @timed solve_SDDP(model, sddpparams, 0)
         V0, t2, m2 = @timed get_bellman_value(model, sddpparams, 1, V[1], model.initialState)
         (upb, costs), t3, m3 = @timed estimate_upper_bound(model, sddpparams, scenarios, pbs)
 
         solvingtime = t1
         simulationtime = t2+t3
         solvingmemory = m1
         simulationmemory = m2+m3
-        
+
         print("Instance \t")
         print("Solving time = ",round(solvingtime,4),"\t")
         print("Solving memory = ", solvingmemory,"\t")
         print("Simulation time = ",round(simulationtime,4),"\t")
         print("Simulation memory = ", simulationmemory,"\t")
         print("Gap < ", round(100*(upb-V0)/V0),"% with prob 97.5%\t")
-        println("number external solver call = ", callsolver)
+        println("number external solver call = ", sddpstats.ncallsolver)
     end
 end

src/noises.jl

@@ -129,40 +129,6 @@ function sampling(law::Vector{NoiseLaw}, t::Int64)
 end
 
 
-"""
-DEPRECATED
-Simulate n scenarios according to a given NoiseLaw
-
-Parameters:
-- law::Vector{NoiseLaw}
-    Vector of discrete independent random variables
-- n::Int
-    number of simulations to compute
-
-Returns :
-- scenarios Array(Float64,n,T)
-    an Array of scenarios, scenarios[i,:] being the ith noise scenario
-"""
-function generate_scenarios(laws::Vector{NoiseLaw}, n::Int64)
-    warn("deprecated generate_scenarios use simulate_scenarios")
-    if n <= 0
-        error("negative number of simulations")
-    end
-    Tf = length(laws)
-    scenarios = Array{Vector{Float64}}(n,Tf)
-    for i = 1:n#TODO can be parallelized
-        scenario = []
-        for t=1:Tf
-            new_val = laws[t].support[:, rand(Categorical(laws[t].proba))]
-            push!(scenario, new_val)
-        end
-        scenarios[i,:]=scenario
-    end
-
-    return scenarios
-end
-
-
 """
 Simulate n scenarios and return a 3D array
 

src/objects.jl

@@ -236,7 +236,20 @@ function set_max_iterations(param::SDDPparameters, n_iter::Int)
     param.maxItNumber = n_iter
 end
 
-
+# Define an object to store evolution of solution
+# along iterations:
+type SDDPStat
+    # Number of iterations:
+    niterations::Int64
+    # evolution of lower bound:
+    lower_bounds::Vector{Float64}
+    # evolution of upper bound:
+    upper_bounds::Vector{Float64}
+    # evolution of execution time:
+    exectime::Vector{Float64}
+    # number of calls to solver:
+    ncallsolver::Int64
+end
 
 
 type NextStep