[DEV] Add SDDPStats Model to return statistics about iterations

Author: frapac

src/SDDPoptimize.jl

@@ -40,16 +40,16 @@ 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
 
 

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/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