[UPD] update stock-example and SDDPoptimize to have a simple SDDP working

Author: Vincent Leclere

examples/stock-example.jl

@@ -72,16 +72,16 @@ paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward path, stop at M
 
 ######### Solving the problem via SDDP
 #V, pbs = solve_SDDP(spmodel, paramSDDP, 10) # display information every 10 iterations
-#lb = StochDynamicProgramming.get_lower_bound(spmodel, params, V)
+#lb = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V)
 #println("Lower bound obtained by SDDP: "*string(lb))
 
 ######### Solving the problem via Dynamic Programming
 
 stateSteps = [0.1]
 controlSteps = [0.1]
 infoStruct = "HD" # noise at time t is known before taking the decision at time t
-monteCarloSize = 1000
-paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, monteCarloSize, infoStruct)
+
+paramSDP = SDPparameters(spmodel, stateSteps, controlSteps,  infoStruct)
 V = sdp_optimize(spmodel,paramSDP)
 
 

src/SDDPoptimize.jl

@@ -123,9 +123,8 @@ function run_SDDP(model::SPModel,
 
         if (display > 0) && (iteration_count%display==0)
             println("Pass number ", iteration_count,
-                    "\tEstimation of upper-bound: ", upper_bound(costs),
-                    "\tLower-bound: ", get_bellman_value(model, param, 1, V[1], model.initialState),
-                    "\tTime: ", toq())
+                    "\tLower-bound: ", round(get_bellman_value(model, param, 1, V[1], model.initialState),4),
+                    "\tTime: ", round(toq(),2))
         end
 
     end
@@ -137,11 +136,11 @@ function run_SDDP(model::SPModel,
 
         println("Estimate upper-bound with Monte-Carlo ...")
         upb, costs = estimate_upper_bound(model, param, V, problems)
-        println("Estimation of upper-bound: ", upb,
-                "\tExact lower bound: ", V0,
-                "\t Gap (\%) <  ", (V0-upb)/V0 , " with prob. > 97.5 \%")
+        println("Estimation of upper-bound: ", round(upb,4),
+                "\tExact lower bound: ", round(V0,4),
+                "\t Gap <  ", round(100*(upb-V0)/V0) , "\%  with prob. > 97.5 \%")
         println("Estimation of cost of the solution (fiability 95\%):",
-                 mean(costs), " +/- ", 1.96*std(costs)/sqrt(length(costs)))
+                 round(mean(costs),4), " +/- ", round(1.96*std(costs)/sqrt(length(costs)),4))
     end
 
     return V, problems
@@ -432,6 +431,28 @@ function get_bellman_value(model::SPModel, param::SDDPparameters,
 end
 
 
+"""
+Compute value of Bellman function at point xt. Return V_t(xt)
+
+Parameters:
+- model (SPModel)
+    Parametrization of the problem
+
+- param (SDDPparameters)
+    Parameters of SDDP
+
+- V (Vector{Polyhedral function})
+    Estimation of bellman function as Polyhedral function
+
+Return:
+current lower bound of the problem (Float64)
+"""
+function get_lower_bound(model::SPModel, param::SDDPparameters,
+                            V::Vector{PolyhedralFunction})
+    return get_bellman_value(model, param, 1, V[1], model.initialState)
+end
+
+
 """
 Add several cuts to JuMP.Model from a PolyhedralFunction
 

src/SDPoptimize.jl

@@ -288,8 +288,8 @@ function sdp_solve_HD(model::SPModel,
     end
 
     if display
-            println("Starting stochastic dynamic programming hazard decision computation")
-        end
+        println("Starting stochastic dynamic programming hazard decision computation")
+    end
 
     #Loop over time
     for t = (TF-1):-1:1

src/objects.jl

@@ -37,6 +37,7 @@ type LinearDynamicLinearCostSPmodel <: SPModel
         for i = 1:dimStates
             push!(xbounds, (-Inf, Inf))
         end
+       
         return new(nstage, dimControls, dimStates, dimNoises, xbounds, ubounds, x0, cost, dynamic, aleas)
     end
 end
@@ -109,7 +110,6 @@ type StochDynProgModel <: SPModel
     end
 
     function StochDynProgModel(model::PiecewiseLinearCostSPmodel, final, cons)
-
         function cost(t,x,u,w)
             saved_cost = -Inf
             current_cost = 0