Merge branch 'stock-example' of github.com:leclere/StochDynamicProgramming.jl into stock-example

Author: frapac

examples/stock-example.jl

@@ -6,16 +6,16 @@
 # Compare different ways of solving a stock problem :
 # Min   E [\sum_{t=1}^TF c_t u_t]
 # s.t.    s_{t+1} = s_t + u_t - xi_t, s_0 given
-#         0 <= s_t <= 1    
+#         0 <= s_t <= 1
 #         u_min <= u_t <= u_max
-#         u_t choosen knowing xi_1 .. xi_t 
+#         u_t choosen knowing xi_1 .. xi_t
 #############################################################################
 
 push!(LOAD_PATH, "../src")
 using StochDynamicProgramming, JuMP, Clp, Distributions
 
 ######## Optimization parameters  ########
-# choose the LP solver used. 
+# choose the LP solver used.
 const SOLVER = ClpSolver()
 # const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0)
 # const SOLVER = GurobiSolver()
@@ -33,7 +33,7 @@ const CONTROL_MIN = 0
 const XI_MAX = 0.3
 const XI_MIN = 0
 const N_XI = 10
-
+# initial stock
 const S0 = 0.5
 
 # create law of noises
@@ -67,7 +67,7 @@ spmodel = LinearDynamicLinearCostSPmodel(N_STAGES,
                                                 xi_laws)
 
 set_state_bounds(spmodel, s_bounds)
-paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward path, stop at MAX_ITER 
+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
@@ -76,17 +76,17 @@ println("Lower bound obtained by SDDP: "*string(lb_sddp))
 
 ######### Solving the problem via Dynamic Programming
 
-stateSteps = [0.1]
-controlSteps = [0.1]
+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)
-V = sdp_optimize(spmodel,paramSDP)
-lb_sdp = StochDynamicProgramming.get_value(spmodel,paramSDP,V)
+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.generate_scenarios(xi_laws,1) 
+#scenarios = StochDynamicProgramming.generate_scenarios(xi_laws,1)
 #costs, stocks = forward_simulations(spmodel, params, V, pbs, scenarios)
 

src/SDDPoptimize.jl

@@ -124,7 +124,7 @@ function run_SDDP(model::SPModel,
         if (display > 0) && (iteration_count%display==0)
             println("Pass number ", iteration_count,
                     "\tLower-bound: ", round(get_bellman_value(model, param, 1, V[1], model.initialState),4),
-                    "\tTime: ", round(toq(),2))
+                    "\tTime: ", round(toq(),2),"s")
         end
 
     end

src/SDPoptimize.jl

@@ -140,34 +140,41 @@ Returns :
 function sdp_optimize(model::SPModel,
                   param::SDPparameters,
                   display=true::Bool)
-    
-    function true_fun(x...)
+
+    function true_fun(t,x,u,w)
         return true
     end
-    function zero_fun(x...)
+    function zero_fun(x)
         return 0
     end
-    
+
     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])
+            end
+            return test
+        end
         if in(:finalCostFunction,fieldnames(model))
-            SDPmodel = StochDynProgModel(model, model.finalCostFunction, true_fun)
+            SDPmodel = StochDynProgModel(model, model.finalCostFunction, cons_fun)
         else
-            SDPmodel = StochDynProgModel(model, zero_fun, true_fun)
+            SDPmodel = StochDynProgModel(model, zero_fun, cons_fun)
         end
     elseif isa(model,StochDynProgModel)
         SDPmodel = model
     else
-        error("cannot build StochDynProgModel from current SPmodel. You need to implement 
+        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
     if (param.infoStructure == "DH")
         V = sdp_solve_DH(SDPmodel, param, display)
     elseif (param.infoStructure == "HD")
         V = sdp_solve_HD(SDPmodel, param, display)
     else
-        error("param.infoStructure is neither 'DH' nor 'HD'") 
+        error("param.infoStructure is neither 'DH' nor 'HD'")
     end
 
     return V
@@ -198,7 +205,7 @@ Returns :
 function sdp_solve_DH(model::StochDynProgModel,
                   param::SDPparameters,
                   display=true::Bool)
-    
+
     TF = model.stageNumber
     next_state = zeros(Float64, model.dimStates)
     law = model.noises
@@ -264,10 +271,10 @@ function sdp_solve_DH(model::StochDynProgModel,
                     w_sample = samples[:, w]
                     proba = probas[w]
                     next_state = model.dynamics(t, x, u, w_sample)
-                    
-                    
+
+
                     if model.constraints(t, next_state, u, w_sample)
-                    
+
                         count_admissible_w = count_admissible_w + proba
                         ind_next_state = real_index_from_variable(next_state, x_bounds, x_steps)
                         next_V = Vitp[ind_next_state...]
@@ -383,7 +390,7 @@ function sdp_solve_HD(model::StochDynProgModel,
             #Compute expectation
             for w in 1:sampling_size
                 admissible_u_w_count = 0
-                best_V_x_w = 0.
+                best_V_x_w = Inf
                 next_V_x_w = Inf
                 w_sample = samples[:, w]
                 proba = probas[w]
@@ -410,7 +417,6 @@ function sdp_solve_HD(model::StochDynProgModel,
                 expected_V += proba*best_V_x_w
                 count_admissible_w += (admissible_u_w_count>0)*proba
             end
-
             if (count_admissible_w>0.)
                 expected_V = expected_V / count_admissible_w
             end
@@ -438,8 +444,9 @@ 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]
+    ind_x0 = real_index_from_variable(model.initialState, model.xlim, param.stateSteps)
+    Vi = value_function_interpolation(model, V, 1)
+    return Vi[ind_x0...,1]
 end
 
 

src/objects.jl

@@ -37,7 +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
@@ -103,7 +103,7 @@ type StochDynProgModel <: SPModel
     noises::Vector{NoiseLaw}
 
     function StochDynProgModel(model::LinearDynamicLinearCostSPmodel, final, cons)
-        return new(model.stageNumber-1, model.dimControls, model.dimStates,
+        return new(model.stageNumber, model.dimControls, model.dimStates,
                  model.dimNoises, model.xlim, model.ulim, model.initialState,
                  model.costFunctions, final, model.dynamics, cons,
                  model.noises)
@@ -122,7 +122,7 @@ type StochDynProgModel <: SPModel
             return saved_cost
         end
 
-        return new(model.stageNumber-1, model.dimControls, model.dimStates,
+        return new(model.stageNumber, model.dimControls, model.dimStates,
                  model.dimNoises, model.xlim, model.ulim, model.initialState,
                  cost, final, model.dynamics, cons,
                  model.noises)