Merge pull request #97 from JuliaOpt/add-benchmark

Add benchmark

Author: leclere

examples/testbenchmark.jl

@@ -0,0 +1,75 @@
+#  Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard
+#  This Source Code Form is subject to the terms of the Mozilla Public
+#  License, v. 2.0. If a copy of the MPL was not distributed with this
+#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#############################################################################
+# Compare different ways of solving a stock problem :
+# Min   E [\sum_{t=1}^TF c_t u_t]
+# s.t.    s_{t+1} = s_t + u_t - xi_t, s_0 given
+#         0 <= s_t <= 1
+#         u_min <= u_t <= u_max
+#         u_t choosen knowing xi_1 .. xi_t
+#############################################################################
+
+using StochDynamicProgramming, JuMP, Clp, Distributions, CPLEX
+println("library loaded")
+
+run_sddp = true
+
+######## Optimization parameters  ########
+# choose the LP solver used.
+#const SOLVER = ClpSolver()
+const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0)
+
+
+# convergence test
+const MAX_ITER = 1 # maximum iteration of SDDP
+
+######## Stochastic Model  Parameters  ########
+const N_STAGES = 5
+const COSTS = rand(N_STAGES)
+
+const CONTROL_MAX = 0.5
+const CONTROL_MIN = 0
+
+const XI_MAX = 0.3
+const XI_MIN = 0
+const N_XI = 10
+# initial stock
+const S0 = 0.5
+
+# create law of noises
+proba = 1/N_XI*ones(N_XI) # uniform probabilities
+xi_support = collect(linspace(XI_MIN,XI_MAX,N_XI))
+xi_law = NoiseLaw(xi_support, proba)
+xi_laws = NoiseLaw[xi_law for t in 1:N_STAGES-1]
+
+# Define dynamic of the stock:
+function dynamic(t, x, u, xi)
+    return [x[1] + u[1] - xi[1]]
+end
+
+# Define cost corresponding to each timestep:
+function cost_t(t, x, u, w)
+    return COSTS[t] * u[1]
+end
+
+######## Setting up the SPmodel
+s_bounds = [(0, 1)]
+u_bounds = [(CONTROL_MIN, CONTROL_MAX)]
+spmodel = LinearDynamicLinearCostSPmodel(N_STAGES, u_bounds, [S0], cost_t, dynamic, xi_laws)
+set_state_bounds(spmodel, s_bounds)
+
+
+######### Define scenarios
+scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws, 1000)
+
+######## Define different scenarios
+paramSDDP1 = SDDPparameters(SOLVER, 4, 0, MAX_ITER) #forwardpassnumber, sensibility
+paramSDDP2 = SDDPparameters(SOLVER, 10, 0, 2)
+
+######## Define parameters collection
+paramSDDP = [paramSDDP1 for i in 1:4]
+
+#Benchmark the collection of parameters
+benchmark_parameters(spmodel, paramSDDP, scenarios, 12)

src/SDDPoptimize.jl

@@ -21,7 +21,7 @@ fulfilled.
     the stochastic problem we want to optimize
 * `param::SDDPparameters`:
     the parameters of the SDDP algorithm
-* `display::Int64)`:
+* `display::Int64`:
     Default is `0`
     If non null, display progression in terminal every
     `n` iterations, where `n` is number specified by display.
@@ -32,20 +32,24 @@ fulfilled.
 * `problems::Array{JuMP.Model}`:
     the collection of linear problems used to approximate
     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:
-    run_SDDP!(model, param, V, problems, display)
-    return V, problems
+    count_callsolver = run_SDDP!(model, param, V, problems, display)
+    return V, problems, count_callsolver
 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)
-    run_SDDP!(model, param, V, problems, display)
-    return V, problems
+    count_callsolver = run_SDDP!(model, param, V, problems, display)
+    return V, problems, count_callsolver
 end
 
 
@@ -56,6 +60,9 @@ function run_SDDP!(model::SPModel,
                     problems::Vector{JuMP.Model},
                     display=0::Int64)
 
+    #Initialization of the counter
+    count_callsolver::Int = 0
+
     if display > 0
       println("Initialize cuts")
     end
@@ -80,20 +87,23 @@ function run_SDDP!(model::SPModel,
         noise_scenarios = simulate_scenarios(model.noises, param.forwardPassNumber)
 
         # Forward pass
-        stockTrajectories = forward_simulations(model,
+        _, stockTrajectories,_,callsolver_forward = forward_simulations(model,
                             param,
                             problems,
-                            noise_scenarios)[2]
+                            noise_scenarios)
+
 
         # Backward pass
-        backward_pass!(model,
+        callsolver_backward = backward_pass!(model,
                       param,
                       V,
                       problems,
                       stockTrajectories,
                       model.noises,
                       false)
 
+        #Update the number of call
+        count_callsolver += callsolver_forward + callsolver_backward
 
         iteration_count += 1
 
@@ -138,6 +148,8 @@ function run_SDDP!(model::SPModel,
         println("Estimation of cost of the solution (fiability 95\%):",
                  round(mean(costs),4), " +/- ", round(1.96*std(costs)/sqrt(length(costs)),4))
     end
+
+    return count_callsolver
 end
 
 
@@ -169,6 +181,7 @@ function estimate_upper_bound(model::SPModel, param::SDDPparameters,
 
     aleas = simulate_scenarios(model.noises, n_simulation)
 
+    callsolver::Int = 0
     costs, stockTrajectories, _ = forward_simulations(model,
                                                         param,
                                                         problem,
@@ -328,7 +341,7 @@ function initialize_value_functions(model::SPModel,
     end
 
 
-    backward_pass!(model,
+    callsolver = backward_pass!(model,
                   param,
                   V,
                   solverProblems,
@@ -552,4 +565,3 @@ function is_cut_relevant(model::SPModel, k::Int, Vt::PolyhedralFunction, solver;
     sol = getobjectivevalue(m)
     return sol < epsilon
 end
-

src/StochDynamicProgramming.jl

@@ -18,7 +18,8 @@ export solve_SDDP,
         PiecewiseLinearCostSPmodel, extensive_formulation,
         PolyhedralFunction, NextStep, forward_simulations,
         StochDynProgModel, SDPparameters, solve_DP,
-        sdp_forward_simulation, sampling, get_control, get_bellman_value
+        sdp_forward_simulation, sampling, get_control, get_bellman_value,
+        benchmark_parameters
 
 include("objects.jl")
 include("utils.jl")
@@ -27,4 +28,5 @@ include("forwardBackwardIterations.jl")
 include("SDDPoptimize.jl")
 include("extensiveFormulation.jl")
 include("SDPoptimize.jl")
+include("compare.jl")
 end

src/compare.jl

@@ -0,0 +1,64 @@
+#  Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard
+#  This Source Code Form is subject to the terms of the Mozilla Public
+#  License, v. 2.0. If a copy of the MPL was not distributed with this
+#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#############################################################################
+#  Compare the optimal values and control returned by different instances
+#  of SDDP on the same problem
+#############################################################################
+
+"""
+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 
+of calls to the solver
+
+# Arguments
+* `model::SPmodel`:
+    The stochastic problem we want to benchmark
+* `SDDParametersCollection::Array{Any,1}`
+    Collection of SDDPparameters
+* `scenarios::Array{Float64,3}`
+    Set of scenarios used to calculate costs
+* `seeds::Int`
+    The random number generator seeds
+    
+# Output
+* `Display in the terminal`
+    Print information in the terminal
+"""
+function benchmark_parameters(model,
+                               SDDParametersCollection,
+                               scenarios::Array{Float64,3},
+                               seeds::Int)
+
+    #Execute a first time each function to compile them
+    (V, pbs, callsolver), t1, m1 = @timed solve_SDDP(model, SDDParametersCollection[1], 0)
+     V0, t2, m2 = @timed get_bellman_value(model, SDDParametersCollection[1], 1, V[1], model.initialState)
+    (upb, costs), t3, m3 = @timed estimate_upper_bound(model, SDDParametersCollection[1], scenarios, pbs)
+
+
+    for sddpparams in SDDParametersCollection
+
+        srand(seeds)
+
+        (V, pbs, callsolver), 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)
+    end
+end

src/forwardBackwardIterations.jl

@@ -33,12 +33,17 @@ scenario according to the current value functions.
 * `controls::Array{Float64, 3}`:
     the simulated controls trajectories. controls(t,k,:) is the control for
     scenario k at time t.
+* `callsolver::Int64`:
+    the number of solver's call'
+
 """
 function forward_simulations(model::SPModel,
                             param::SDDPparameters,
                             solverProblems::Vector{JuMP.Model},
                             xi::Array{Float64})
 
+    callsolver::Int = 0
+
     T = model.stageNumber
     nb_forward = size(xi)[2]
 
@@ -69,6 +74,7 @@ function forward_simulations(model::SPModel,
             state_t = collect(stocks[t, k, :])
             alea_t = collect(xi[t, k, :])
 
+            callsolver += 1
             status, nextstep = solve_one_step_one_alea(
                                         model,
                                         param,
@@ -89,7 +95,7 @@ function forward_simulations(model::SPModel,
             end
         end
     end
-    return costs, stocks, controls
+    return costs, stocks, controls, callsolver
 end
 
 
@@ -174,6 +180,8 @@ function backward_pass!(model::SPModel,
                         law,
                         init=false::Bool)
 
+    callsolver::Int = 0
+
     T = model.stageNumber
     nb_forward = size(stockTrajectories)[2]
 
@@ -193,6 +201,7 @@ function backward_pass!(model::SPModel,
 
                 alea_t  = collect(law[t].support[:, w])
 
+                callsolver += 1
                 solved, nextstep = solve_one_step_one_alea(model, param, solverProblems[t], t, state_t, alea_t)
                 if solved
                     subgradient_array[:, w] = nextstep.sub_gradient
@@ -225,4 +234,5 @@ function backward_pass!(model::SPModel,
 
         end
     end
+    return callsolver
 end

test/runtests.jl

@@ -20,4 +20,3 @@ include("sddp.jl")
 
 # Test DP:
 include("sdp.jl")
-

test/sddp.jl

@@ -173,6 +173,12 @@ facts("SDDP algorithm: 1D case") do
         @fact V[1].betas --> Vdump[1].betas
         @fact V[1].lambdas --> Vdump[1].lambdas
     end
+
+    context("Compare parameters") do
+        paramSDDP = [params for i in 1:3]
+        scenarios = StochDynamicProgramming.simulate_scenarios(laws, 1000)
+        benchmark_parameters(model, paramSDDP, scenarios, 12)
+    end
 end