[UPD] update example and adapt SDPoptimize accordingly

Author: Vincent Leclere

examples/stock-example.jl

@@ -10,15 +10,17 @@
 #         u_min <= u_t <= u_max
 #         u_t choosen knowing xi_1 .. xi_t
 #############################################################################
-
 push!(LOAD_PATH, "../src")
 using StochDynamicProgramming, JuMP, Clp, Distributions
+println("library loaded")
+
+run_sddp = false
+run_sdp = true
 
 ######## Optimization parameters  ########
 # choose the LP solver used.
 const SOLVER = ClpSolver()
-# const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0)
-# const SOLVER = GurobiSolver()
+# const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) # require "using CPLEX"
 
 # convergence test
 const MAX_ITER = 100 # maximum iteration of SDDP
@@ -40,10 +42,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-1)
-for t=1:N_STAGES-1
-   xi_laws[t] = xi_law
-end
+xi_laws = NoiseLaw[xi_law for t in 1:N_STAGES-1] 
 
 # Define dynamic of the stock:
 function dynamic(t, x, u, xi)
@@ -55,38 +54,42 @@ 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)
-paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward path, stop at MAX_ITER
+    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)
+    
 
 ######### 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))
+if run_sddp
+    paramSDDP = SDDPparameters(SOLVER, 2, 0, MAX_ITER) # 10 forward pass, stop at MAX_ITER
+    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))
+end
 
 ######### Solving the problem via Dynamic Programming
+if run_sdp
+    stateSteps = [0.01]
+    controlSteps = [0.01]
+    infoStruct = "HD" # noise at time t is known before taking the decision at time t
+
+    paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, infoStruct)
+    Vs = sdp_optimize(spmodel,paramSDP)
+    lb_sdp = StochDynamicProgramming.get_value(spmodel,paramSDP,Vs)
+    println("Value obtained by SDP: "*string(lb_sdp))
+end
 
-stateSteps = [0.01]
-controlSteps = [0.01]
-infoStruct = "HD" # noise at time t is known before taking the decision at time t
-
-paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, infoStruct)
-Vs = sdp_optimize(spmodel,paramSDP)
-lb_sdp = StochDynamicProgramming.get_value(spmodel,paramSDP,Vs)
-println("Lower bound obtained by SDP: "*string(lb_sdp))
+######### Comparing the solution
+scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,2)
+if run_sddp
+    costs, stocks = forward_simulations(spmodel, paramSDDP, V, pbs, scenarios)
+end
+if run_sdp
+    costs, states, stocks =sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs)
+end
+println(costs)
 
 
-######### Comparing the solution
-#scenarios = StochDynamicProgramming.generate_scenarios(xi_laws,1)
-#costs, stocks = forward_simulations(spmodel, params, V, pbs, scenarios)
 

src/SDDPoptimize.jl

@@ -94,10 +94,7 @@ function run_SDDP(model::SPModel,
 
         # Build given number of scenarios according to distribution
         # law specified in model.noises:
-        noise_scenarios = simulate_scenarios(model.noises ,
-                                    (model.stageNumber-1,
-                                     param.forwardPassNumber,
-                                     model.dimNoises))
+        noise_scenarios = simulate_scenarios(model.noises, param.forwardPassNumber)
 
         # Forward pass
         costs, stockTrajectories, _ = forward_simulations(model,

src/SDPoptimize.jl

@@ -117,30 +117,9 @@ function generate_grid(model::SPModel, param::SDPparameters)
 end
 
 """
-Value iteration algorithm to compute optimal value functions in
-the Decision Hazard (DH) as well as the Hazard Decision (HD) case
-
-Parameters:
-- model (SPmodel)
-    the DPSPmodel of our problem
-
-- param (SDPparameters)
-    the parameters for the SDP algorithm
-
-- display (Bool)
-    the output display or verbosity parameter
-
-
-Returns :
-- value_functions (Array)
-    the vector representing the value functions as functions of the state
-    of the system at each time step
-
+Try to construct a StochDynProgModel from an SPModel
 """
-function sdp_optimize(model::SPModel,
-                  param::SDPparameters,
-                  display=true::Bool)
-
+function build_sdpmodel_from_spmodel(model::SPModel)
     function true_fun(t,x,u,w)
         return true
     end
@@ -150,11 +129,12 @@ function sdp_optimize(model::SPModel,
 
     if isa(model,PiecewiseLinearCostSPmodel)||isa(model,LinearDynamicLinearCostSPmodel)
         function cons_fun(t,x,u,w)
-            test = true
             for i in 1:model.dimStates
-                test &= (x[i]>=model.xlim[i][1])&(x[i]<=model.xlim[i][2])
+                if (x[i]<=model.xlim[i][1]) || (x[i]>=model.xlim[i][2])
+                    return false
+                end
             end
-            return test
+            return true
         end
         if in(:finalCostFunction,fieldnames(model))
             SDPmodel = StochDynProgModel(model, model.finalCostFunction, cons_fun)
@@ -167,6 +147,35 @@ function sdp_optimize(model::SPModel,
         error("cannot build StochDynProgModel from current SPmodel. You need to implement
         a new StochDynProgModel constructor.")
     end
+    return SDPmodel
+end
+
+"""
+Value iteration algorithm to compute optimal value functions in
+the Decision Hazard (DH) as well as the Hazard Decision (HD) case
+
+Parameters:
+- model (SPmodel)
+    the DPSPmodel of our problem
+
+- param (SDPparameters)
+    the parameters for the SDP algorithm
+
+- display (Bool)
+    the output display or verbosity parameter
+
+
+Returns :
+- value_functions (Array)
+    the vector representing the value functions as functions of the state
+    of the system at each time step
+
+"""
+function sdp_optimize(model::SPModel,
+                  param::SDPparameters,
+                  display=true::Bool)
+
+    SDPmodel = build_sdpmodel_from_spmodel(model::SPModel)
 
     #Display start of the algorithm in DH and HD cases
     if (param.infoStructure == "DH")
@@ -449,6 +458,58 @@ function get_value(model::SPModel,param::SDPparameters,V::Array{Float64})
     return Vi[ind_x0...,1]
 end
 
+"""
+Simulation of optimal trajectories given model and Bellman functions
+
+Parameters:
+- model (SPmodel)
+    the DPSPmodel of our problem
+
+- param (SDPparameters)
+    the parameters for the SDP algorithm
+
+- scenarios (Array)
+    the scenarios of uncertainties realizations we want to simulate on
+    scenarios[t,k,:] is the alea at time t for scenario k
+
+- value_functions (Array)
+    the vector representing the value functions as functions of the state
+    of the system at each time step
+
+- display (Bool)
+    the output display or verbosity parameter
+
+Returns :
+
+- costs (Vector{Float64})
+    the cost of the optimal control over the scenario provided
+
+- stocks (Array{Float64})
+    the state of the controlled system at each time step
+
+- controls (Array{Float64})
+    the controls applied to the system at each time step
+"""
+function sdp_forward_simulation(model::SPModel,
+                  param::SDPparameters,
+                  scenarios::Array{Float64,3},
+                  value::Array,
+                  display=true::Bool)
+                  
+    SDPmodel = build_sdpmodel_from_spmodel(model)              
+    TF = SDPmodel.stageNumber 
+    nb_scenarios = size(scenarios)[2] 
+             
+    costs = zeros(nb_scenarios)
+    controls = zeros(TF,nb_scenarios)
+    states = zeros(TF-1,nb_scenarios)
+    
+    for k = 1:nb_scenarios
+        costs[k],states[:,k], controls[:,k]= sdp_forward_simulation(SDPmodel,
+                  param,scenarios[:,k],model.initialState,value,display)
+    end
+    return costs, controls, states
+end
 
 """
 Simulation of optimal control given an initial state and an alea scenario
@@ -485,7 +546,7 @@ Returns :
     the controls applied to the system at each time step
 
 """
-function sdp_forward_simulation(model::SPModel,
+function sdp_forward_simulation(model::StochDynProgModel,
                   param::SDPparameters,
                   scenario::Array,
                   X0::Array,

src/forwardBackwardIterations.jl

@@ -58,15 +58,23 @@ function forward_simulations(model::SPModel,
                             param::SDDPparameters,
                             V::Vector{PolyhedralFunction},
                             solverProblems::Vector{JuMP.Model},
-                            xi::Array{Float64, 3},
+                            xi::Array{Float64},
                             returnCosts=true::Bool,
                             init=false::Bool,
                             display=false::Bool)
-
-    # TODO: verify that loops are in the same order
+    
     T = model.stageNumber
     nb_forward = size(xi)[2] 
 
+    if ndims(xi)!=3
+        if ndims(xi)==2
+            warn("noise scenario are not given in the right shape. Assumed to be real valued noise.")
+            xi = reshape(xi,(T,nb_forward,1)) 
+        else
+            error("wrong dimension of noise scenarios")
+        end
+     end
+
     stocks = zeros(T, nb_forward, model.dimStates)
     # We got T - 1 control, as terminal state is included into the total number
     # of stages.

src/noises.jl

@@ -155,6 +155,7 @@ end
 
 
 """
+DEPRECATED
 Simulate n scenarios according to a given NoiseLaw
 
 Parameters:
@@ -170,6 +171,7 @@ Returns :
     an Array of scenarios, scenarios[i,:] being the ith noise scenario
 """
 function generate_scenarios(laws::Vector{NoiseLaw}, n::Int64)
+    warn("deprecated generate_scenarios use simulate_scenarios")
     if n <= 0
         error("negative number of simulations")
     end
@@ -187,8 +189,43 @@ function generate_scenarios(laws::Vector{NoiseLaw}, n::Int64)
     return scenarios
 end
 
+"""
+Simulate n scenarios and return a 3D array
+
+Parameters:
+- laws (Vector{NoiseLaw})
+    Distribution laws corresponding to each timestep
+
+- n (Int64)
+    number of scenarios to simulate
 
+Return:
+- scenarios Array{Float64, 3}
+    scenarios[t,k,:] is the noise at time t for scenario k
 """
+function simulate_scenarios(laws, n::Int64)
+    T = length(laws)
+    dimAlea = size(laws[1].support)[1]
+    dims =(T,n,dimAlea)
+    if typeof(laws) == Distributions.Normal
+        scenarios = rand(laws, dims)
+    else
+        scenarios = zeros(dims)
+
+        for k=1:dims[2]
+            for t=1:dims[1]
+                gen = Categorical(laws[t].proba)
+                scenarios[t, k, :] = laws[t].support[:, rand(gen)]
+            end
+
+        end
+    end
+
+    return scenarios
+end
+
+"""
+DEPRECATED
 Simulate n scenarios and return a 3D array
 
 
@@ -205,7 +242,7 @@ Return:
 
 """
 function simulate_scenarios(laws, dims::Tuple)
-
+    warn("decrecated call to simulate_scenarios")
     if typeof(laws) == Distributions.Normal
         scenarios = rand(laws, dims)
     else