Merge branch 'sdp_getters3'

Author: trigaut

src/SDPoptimize.jl

@@ -117,12 +117,22 @@ function generate_grid(model::SPModel, param::SDPparameters)
 end
 
 """
-Try to construct a StochDynProgModel from an SPModel
+Transform a general SPmodel into a StochDynProgModel
+
+Parameters:
+- model (SPmodel)
+    the model of the problem
+
+- param (SDPparameters)
+    the parameters of the problem
+
+
+Returns :
+- sdpmodel : (StochDynProgModel)
+    the corresponding StochDynProgModel
 """
 function build_sdpmodel_from_spmodel(model::SPModel)
-    function true_fun(t,x,u,w)
-        return true
-    end
+
     function zero_fun(x)
         return 0
     end
@@ -147,6 +157,7 @@ function build_sdpmodel_from_spmodel(model::SPModel)
         error("cannot build StochDynProgModel from current SPmodel. You need to implement
         a new StochDynProgModel constructor.")
     end
+
     return SDPmodel
 end
 
@@ -276,7 +287,6 @@ function sdp_solve_DH(model::StochDynProgModel,
                 end
 
                 for w = 1:sampling_size
-
                     w_sample = samples[:, w]
                     proba = probas[w]
                     next_state = model.dynamics(t, x, u, w_sample)
@@ -495,16 +505,16 @@ function sdp_forward_simulation(model::SPModel,
                   scenarios::Array{Float64,3},
                   value::Array,
                   display=true::Bool)
-                  
-    SDPmodel = build_sdpmodel_from_spmodel(model)              
-    TF = SDPmodel.stageNumber 
-    nb_scenarios = size(scenarios)[2] 
-             
+
+    SDPmodel = build_sdpmodel_from_spmodel(model)
+    TF = SDPmodel.stageNumber
+    nb_scenarios = size(scenarios)[2]
+
     costs = zeros(nb_scenarios)
     states = zeros(TF,nb_scenarios)
     controls = zeros(TF-1,nb_scenarios)
-    
-    
+
+
     for k = 1:nb_scenarios
         #println(k)
         costs[k],states[:,k], controls[:,k] = sdp_forward_single_simulation(SDPmodel,
@@ -514,6 +524,148 @@ function sdp_forward_simulation(model::SPModel,
     return costs, controls, states
 end
 
+"""
+Get the optimal control at time t knowing the state of the system in the decision hazard case
+
+Parameters:
+- model (SPmodel)
+    the DPSPmodel of our problem
+
+- param (SDPparameters)
+    the parameters for the SDP algorithm
+
+- V (Array{Float64})
+    the Bellman Functions
+
+- t (int)
+    the time step
+
+- x (Array)
+    the state variable
+
+- w (Array)
+the alea realization
+
+Returns :
+- V(x0) (Float64)
+"""
+function get_control(model::SPModel,param::SDPparameters,V::Array{Float64}, t::Int64, x::Array)
+
+    SDPmodel = build_sdpmodel_from_spmodel(model)
+
+    product_controls = product([SDPmodel.ulim[i][1]:param.controlSteps[i]:SDPmodel.ulim[i][2] for i in 1:SDPmodel.dimControls]...)
+
+    law = SDPmodel.noises
+    best_control = tuple()
+    Vitp = value_function_interpolation(SDPmodel, V, t+1)
+
+    u_bounds = SDPmodel.ulim
+    x_bounds = SDPmodel.xlim
+    x_steps = param.stateSteps
+
+    best_V = Inf
+
+    for u in product_controls
+
+        count_admissible_w = 0.
+        current_V = 0.
+
+        if (param.expectation_computation=="MonteCarlo")
+            sampling_size = param.monteCarloSize
+            samples = [sampling(law,t) for i in 1:sampling_size]
+            probas = (1/sampling_size)
+        else
+            sampling_size = law[t].supportSize
+            samples = law[t].support
+            probas = law[t].proba
+        end
+
+        for w = 1:sampling_size
+
+            w_sample = samples[:, w]
+            proba = probas[w]
+
+            next_state = SDPmodel.dynamics(t, x, u, w_sample)
+
+            if SDPmodel.constraints(t, next_state, u, w_sample)
+                ind_next_state = real_index_from_variable(next_state, x_bounds, x_steps)
+                next_V = Vitp[ind_next_state...]
+                current_V += proba *(SDPmodel.costFunctions(t, x, u, w_sample) + next_V)
+                count_admissible_w = count_admissible_w + proba
+            end
+        end
+        current_V = current_V/count_admissible_w
+        if (current_V < best_V)&(count_admissible_w>0)
+            best_control = u
+            best_V = current_V
+        end
+    end
+
+    return best_control
+end
+
+"""
+Get the optimal control at time t knowing the state of the system and the alea in the hazard decision case
+
+Parameters:
+- model (SPmodel)
+    the DPSPmodel of our problem
+
+- param (SDPparameters)
+    the parameters for the SDP algorithm
+
+- V (Array{Float64})
+    the Bellman Functions
+
+- t (int)
+    the time step
+
+- x (Array)
+    the state variable
+
+- w (Array)
+the alea realization
+
+Returns :
+- V(x0) (Float64)
+"""
+function get_control(model::SPModel,param::SDPparameters,V::Array{Float64}, t::Int64, x::Array, w::Array)
+
+    SDPmodel = build_sdpmodel_from_spmodel(model)
+
+    product_controls = product([SDPmodel.ulim[i][1]:param.controlSteps[i]:SDPmodel.ulim[i][2] for i in 1:SDPmodel.dimControls]...)
+
+    law = SDPmodel.noises
+    best_control = tuple()
+    Vitp = value_function_interpolation(SDPmodel, V, t+1)
+
+    u_bounds = SDPmodel.ulim
+    x_bounds = SDPmodel.xlim
+    x_steps = param.stateSteps
+
+    best_V = Inf
+
+    for u = product_controls
+
+        next_state = SDPmodel.dynamics(t, x, u, w)
+
+        if SDPmodel.constraints(t, next_state, u, w)
+            ind_next_state = real_index_from_variable(next_state, x_bounds, x_steps)
+            next_V = Vitp[ind_next_state...]
+            current_V = SDPmodel.costFunctions(t, x, u, w) + next_V
+            if (current_V < best_V)
+                best_control = u
+                best_state = SDPmodel.dynamics(t, x, u, w)
+                best_V = current_V
+            end
+        end
+
+    end
+
+    return best_control
+
+end
+
 """
 Simulation of optimal control given an initial state and an alea scenario
 
@@ -599,21 +751,21 @@ function sdp_forward_single_simulation(model::StochDynProgModel,
                     probas = (1/sampling_size)
                 else
                     sampling_size = law[t].supportSize
-                    samples = law[t].support[:]
+                    samples = law[t].support
                     probas = law[t].proba
                 end
 
                 for w = 1:sampling_size
 
-                    w_sample = samples[w]
+                    w_sample = samples[:, w]
                     proba = probas[w]
 
                     next_state = model.dynamics(t, x, u, w_sample)
 
-                    if model.constraints(t, next_state, u, scenario[t])
+                    if model.constraints(t, next_state, u, w_sample)
                         ind_next_state = real_index_from_variable(next_state, x_bounds, x_steps)
                         next_V = Vitp[ind_next_state...]
-                        current_V += proba *(model.costFunctions(t, x, u, scenario[t]) + next_V)
+                        current_V += proba *(model.costFunctions(t, x, u, w_sample) + next_V)
                         count_admissible_w = count_admissible_w + proba
                     end
                 end

src/StochDynamicProgramming.jl

@@ -17,7 +17,7 @@ export solve_SDDP, NoiseLaw, simulate_scenarios,
         PiecewiseLinearCostSPmodel,
         PolyhedralFunction, NextStep, forward_simulations,
         StochDynProgModel, SDPparameters, sdp_optimize,
-        sdp_forward_simulation, sampling
+        sdp_forward_simulation, sampling, get_control, get_value
 
 include("objects.jl")
 include("utils.jl")

src/objects.jl

@@ -103,8 +103,7 @@ type StochDynProgModel <: SPModel
     noises::Vector{NoiseLaw}
 
     function StochDynProgModel(model::LinearDynamicLinearCostSPmodel, final, cons)
-        return new(model.stageNumber, model.dimControls, model.dimStates,
-                 model.dimNoises, model.xlim, model.ulim, model.initialState,
+        return StochDynProgModel(model.stageNumber, model.xlim, model.ulim, model.initialState,
                  model.costFunctions, final, model.dynamics, cons,
                  model.noises)
     end
@@ -122,16 +121,13 @@ type StochDynProgModel <: SPModel
             return saved_cost
         end
 
-        return new(model.stageNumber, model.dimControls, model.dimStates,
-                 model.dimNoises, model.xlim, model.ulim, model.initialState,
-                 cost, final, model.dynamics, cons,
-                 model.noises)
+        return StochDynProgModel(model.stageNumber, model.xlim, model.ulim, model.initialState,
+                 cost, final, model.dynamics, cons, model.noises)
     end
 
-    function StochDynProgModel(TF, N_CONTROLS, N_STATES, N_NOISES,
-                    x_bounds, u_bounds, x0, cost_t, finalCostFunction, dynamic,
-                    constraints, aleas)
-        return new(TF, N_CONTROLS, N_STATES, N_NOISES,
+    function StochDynProgModel(TF, x_bounds, u_bounds, x0, cost_t,
+                                finalCostFunction, dynamic, constraints, aleas)
+        return new(TF, length(u_bounds), length(x_bounds), length(aleas[1].support[:, 1]),
                     x_bounds, u_bounds, x0, cost_t, finalCostFunction, dynamic,
                     constraints, aleas)
     end