[UPD] stock-example run for both sdp and sddp. compare only lower bounds.

Author: Vincent Leclere

examples/stock-example.jl

@@ -41,7 +41,7 @@ const S0 = 0.5
 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 = Vector{NoiseLaw}(N_STAGES)
+xi_laws = Vector{NoiseLaw}(N_STAGES-1)
 for t=1:N_STAGES-1
    xi_laws[t] = xi_law
 end
@@ -71,9 +71,9 @@ set_state_bounds(spmodel, s_bounds)
 paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward path, stop at MAX_ITER 
 
 ######### Solving the problem via SDDP
-#V, pbs = solve_SDDP(spmodel, paramSDDP, 10) # display information every 10 iterations
-#lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V)
-#println("Lower bound obtained by SDDP: "*string(lb_sddp))
+V, pbs = solve_SDDP(spmodel, paramSDDP, 10) # display information every 10 iterations
+lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V)
+println("Lower bound obtained by SDDP: "*string(lb_sddp))
 
 ######### Solving the problem via Dynamic Programming
 
@@ -83,10 +83,11 @@ infoStruct = "HD" # noise at time t is known before taking the decision at time
 
 paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, infoStruct)
 V = sdp_optimize(spmodel,paramSDP)
-
+lb_sdp = StochDynamicProgramming.get_value(spmodel,paramSDP,V)
+println("Lower bound obtained by SDP: "*string(lb_sdp))
 
 
 ######### Comparing the solution
-#scenarios = generate_scenarios(xi_laws,1000) 
+#scenarios = StochDynamicProgramming.generate_scenarios(xi_laws,1) 
 #costs, stocks = forward_simulations(spmodel, params, V, pbs, scenarios)
 

src/SDPoptimize.jl

@@ -157,8 +157,8 @@ function sdp_optimize(model::SPModel,
     elseif isa(model,StochDynProgModel)
         SDPmodel = model
     else
-        error("cannot build StochDynProgModel from current SPmodel. Implement new
-        StochDynProgModel constructor.")
+        error("cannot build StochDynProgModel from current SPmodel. You need to implement 
+        a new StochDynProgModel constructor.")
     end
 
     #Display start of the algorithm in DH and HD cases
@@ -421,7 +421,26 @@ function sdp_solve_HD(model::StochDynProgModel,
     return V
 end
 
+"""
+Get the optimal value of the problem from the optimal Bellman Function
+
+Parameters:
+- model (SPmodel)
+    the DPSPmodel of our problem
+
+- param (SDPparameters)
+    the parameters for the SDP algorithm
+
+- V (Array{Float64})
+    the Bellman Functions
 
+Returns :
+- V(x0) (Float64)
+"""
+function get_value(model::SPModel,param::SDPparameters,V::Array{Float64})
+    ind_x0 = index_from_variable(model.initialState, model.xlim, param.stateSteps)
+    return V[ind_x0...,1]
+end
 
 
 """

src/noises.jl

@@ -173,12 +173,12 @@ function generate_scenarios(laws::Vector{NoiseLaw}, n::Int64)
     if n <= 0
         error("negative number of simulations")
     end
-    Tf = length(law)
+    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(law[t].proba))]
+            new_val = laws[t].support[:, rand(Categorical(laws[t].proba))]
             push!(scenario, new_val)
         end
         scenarios[i,:]=scenario