[FIX] test and stock examples

Author: leclere

examples/multistock-example.jl

@@ -17,6 +17,7 @@ println("library loaded")
 
 run_sddp = true # false if you don't want to run sddp
 run_sdp  = true # false if you don't want to run sdp
+test_simulation = true # false if you don't want to test your strategies
 
 ######## Optimization parameters  ########
 # choose the LP solver used.
@@ -26,6 +27,8 @@ const SOLVER = ClpSolver() 			   # require "using Clp"
 # convergence test
 const MAX_ITER = 10 # number of iterations of SDDP
 
+const step = 0.1   # discretization step of SDP
+
 ######## Stochastic Model  Parameters  ########
 const N_STAGES = 6              # number of stages of the SP problem
 const N_STOCKS = 3              # number of stocks of the SP problem
@@ -40,9 +43,9 @@ const XI_MAX = 0.3              # bounds on the noise
 const XI_MIN = 0
 const N_XI = 3                 # discretization of the noise
 
-const r = 1                  # bound on cumulative control : \sum_{i=1}^N u_i < rN
+const r = 0.5                  # bound on cumulative control : \sum_{i=1}^N u_i < rN
 
-const S0 = [0.5*r for i=1:N_STOCKS]     # initial stock
+const S0 = [0.5 for i=1:N_STOCKS]     # initial stock
 
 # create law of noises
 proba = 1/N_XI*ones(N_XI) # uniform probabilities
@@ -79,7 +82,6 @@ if run_sddp
     # 10 forward pass, stop at MAX_ITER
     paramSDDP = SDDPparameters(SOLVER,
                                passnumber=10,
-                               gap=0,
                                max_iterations=MAX_ITER)
     V, pbs = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
     lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V)
@@ -90,7 +92,6 @@ end
 ######### Solving the problem via Dynamic Programming
 if run_sdp
     tic()
-    step = 0.1
     println("Starting resolution by SDP")
     stateSteps = [step for i=1:N_STOCKS] # discretization step of the state
     controlSteps = [step for i=1:N_STOCKS] # discretization step of the control
@@ -108,15 +109,11 @@ end
 
 ######### Comparing the solutions on simulated scenarios.
 #srand(1234) # to fix the random seed accross runs
-scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,1000)
-if run_sddp
-    costsddp, stocks_sddp = forward_simulations(spmodel, paramSDDP, pbs, scenarios)
-end
-
-if run_sdp
-    costsdp, stocks_sdp, controls = sdp_forward_simulation(spmodel_sdp,paramSDP,scenarios,Vs)
-end
-
-run_sddp && run_sdp && println("Simulated relative gain of sddp over sdp: "
+if run_sddp && run_sdp && test_simulation
+    scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,1000)
+    costsddp, stocks = forward_simulations(spmodel, paramSDDP, pbs, scenarios)
+    costsdp, states, controls = sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs)
+    println("Simulated relative gain of sddp over sdp: "
             *string(round(200*mean(costsdp-costsddp)/abs(mean(costsddp+costsdp)),3))*"%")
+end
 

examples/stock-example.jl

@@ -17,6 +17,7 @@ println("library loaded")
 run_sddp = true # false if you don't want to run sddp
 run_sdp  = true # false if you don't want to run sdp
 run_ef   = true # false if you don't want to run extensive formulation
+test_simulation = true # false if you don't want to test your strategies
 
 ######## Optimization parameters  ########
 # choose the LP solver used.
@@ -25,13 +26,13 @@ const SOLVER = ClpSolver() 			   # require "using Clp"
 
 # convergence test
 const MAX_ITER = 10 # number of iterations of SDDP
+const step = 0.1   # discretization step of SDP
 
 ######## Stochastic Model  Parameters  ########
 const N_STAGES = 6              # number of stages of the SP problem
 const COSTS = [sin(3*t)-1 for t in 1:N_STAGES]
 #const COSTS = rand(N_STAGES)    # randomly generating deterministic costs
 
-
 const CONTROL_MAX = 0.5         # bounds on the control
 const CONTROL_MIN = 0
 
@@ -70,9 +71,8 @@ if run_sddp
     println("Starting resolution by SDDP")
     # 10 forward pass, stop at MAX_ITER
     paramSDDP = SDDPparameters(SOLVER,
-                                                    passnumber=10,
-                                                    gap=0,
-                                                    max_iterations=MAX_ITER)
+                               passnumber=10,
+                               max_iterations=MAX_ITER)
     V, pbs = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
     lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V)
     println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
@@ -82,7 +82,6 @@ end
 ######### Solving the problem via Dynamic Programming
 if run_sdp
     tic()
-    step = 0.01
     println("Starting resolution by SDP")
     stateSteps = [step] # discretization step of the state
     controlSteps = [step] # discretization step of the control
@@ -100,23 +99,18 @@ if run_ef
     println("Starting resolution by Extensive Formulation")
     value_ef = extensive_formulation(spmodel, paramSDDP)[1]
     println("Value obtained by Extensive Formulation: "*string(round(value_ef,4)))
-    println("Relative error of SDP value: "*string(100*round(value_sdp/value_ef-1,4))*"%")
-    println("Relative error of SDDP lower bound: "*string(100*round(lb_sddp/value_ef-1,4))*"%")
+    println("Relative error of SDP value: "*string(100*round(abs(value_sdp/value_ef)-1,4))*"%")
+    println("Relative error of SDDP lower bound: "*string(100*round(abs(lb_sddp/value_ef)-1,4))*"%")
     toc(); println();
 end
 
 ######### Comparing the solutions on simulated scenarios.
 
 #srand(1234) # to fix the random seed accross runs
-scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,1000)
-if run_sddp
+if run_sddp && run_sdp && test_simulation
+    scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,1000)
     costsddp, stocks = forward_simulations(spmodel, paramSDDP, pbs, scenarios)
-end
-if run_sdp
-    costsdp, states, controls =sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs)
-end
-
-if run_sddp && run_sdp
+    costsdp, states, controls = sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs)
     println("Simulated relative gain of sddp over sdp: "
             *string(round(200*mean(costsdp-costsddp)/abs(mean(costsddp+costsdp)),3))*"%")
 end

src/SDDPoptimize.jl

@@ -10,10 +10,10 @@
 
 
 """
-Solve SDDP algorithm and return estimation of bellman functions.
+Solve spmodel using SDDP algorithm and return lower approximation of Bellman functions.
 
 # Description
-Alternate forward and backward phase till the stopping criterion is
+Alternate forward and backward phases untill the stopping criterion is
 fulfilled.
 
 # Arguments
@@ -32,7 +32,7 @@ fulfilled.
 * `problems::Array{JuMP.Model}`:
     the collection of linear problems used to approximate
     each value function
-* `sddp_stats::SDDPStat`:
+* `sddp_stats::SDDPStat`: statistics of the algorithm run
 
 """
 function solve_SDDP(model::SPModel, param::SDDPparameters, verbose=0::Int64)
@@ -46,10 +46,11 @@ function solve_SDDP(model::SPModel, param::SDDPparameters, verbose=0::Int64)
 end
 
 """
-Solve SDDP algorithm with hotstart and return estimation of bellman functions.
+Solve spmodel using SDDP algorithm and return lower approximation of Bellman functions.
+Use hotstart.
 
 # Description
-Alternate forward and backward phase till the stopping criterion is
+Alternate forward and backward phases untill the stopping criterion is
 fulfilled.
 
 # Arguments
@@ -70,7 +71,7 @@ fulfilled.
 * `problems::Array{JuMP.Model}`:
     the collection of linear problems used to approximate
     each value function
-* `sddp_stats::SDDPStat`:
+* `sddp_stats::SDDPStat`: statistics of the algorithm run
 """
 function solve_SDDP(model::SPModel, param::SDDPparameters, V::Vector{PolyhedralFunction}, verbose=0::Int64)
     check_SDDPparameters(model,param,verbose)

src/objects.jl

@@ -159,7 +159,7 @@ type SDDPparameters
     confidence_level::Float64
     # Maximum iterations of the SDDP algorithms:
     maxItNumber::Int64
-    # Prune cuts every %% iterations:
+    # Define the pruning method
     pruning::Dict{Symbol, Any}
     # Estimate upper-bound every %% iterations:
     compute_ub::Int64
@@ -175,7 +175,7 @@ type SDDPparameters
     function SDDPparameters(solver; passnumber=10, gap=0., confidence=.975,
                             max_iterations=20, prune_cuts=0,
                             pruning_algo="none",
-                            compute_ub=-1, montecarlo_final=-1, montecarlo_in_iter=100,
+                            compute_ub=-1, montecarlo_final=1000, montecarlo_in_iter=100,
                             mipsolver=nothing,
                             rho0=0., alpha=1.)