Merge pull request #108 from JuliaOpt/dev-nightly

Update stock-example

Author: leclere

examples/stock-example.jl

@@ -14,29 +14,30 @@
 using StochDynamicProgramming, JuMP, Clp, Distributions
 println("library loaded")
 
-run_sddp = true
-run_sdp = true
+run_sddp = true # false if you don't want to run sddp
+run_sdp  = true # false if you don't want to run sdp
+run_ef   = true # false if you don't want to run extensive formulation
 
 ######## Optimization parameters  ########
 # choose the LP solver used.
 const SOLVER = ClpSolver()
 # const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) # require "using CPLEX"
 
 # convergence test
-const MAX_ITER = 100 # maximum iteration of SDDP
+const MAX_ITER = 10 # number of iterations of SDDP
 
 ######## Stochastic Model  Parameters  ########
-const N_STAGES = 5
-const COSTS = rand(N_STAGES)
+const N_STAGES = 6              # number of stages of the SP problem
+const COSTS = rand(N_STAGES)    # generating deterministic costs
 
-const CONTROL_MAX = 0.5
+const CONTROL_MAX = 0.5         # bounds on the control
 const CONTROL_MIN = 0
 
-const XI_MAX = 0.3
+const XI_MAX = 0.3              # bounds on the noise
 const XI_MIN = 0
-const N_XI = 10
-# initial stock
-const S0 = 0.5
+const N_XI = 10                 # discretization of the noise
+
+const S0 = 0.5                  # initial stock
 
 # create law of noises
 proba = 1/N_XI*ones(N_XI) # uniform probabilities
@@ -55,33 +56,42 @@ function cost_t(t, x, u, w)
 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)
-
+s_bounds = [(0, 1)] # bounds on the state
+u_bounds = [(CONTROL_MIN, CONTROL_MAX)] # bounds on controls
+spmodel = LinearDynamicLinearCostSPmodel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws)
+set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
+println("Model set up")
 
 ######### Solving the problem via 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
+    println("Starting resolution by SDDP")
+    paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward pass, stop at MAX_ITER
+    V, pbs = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
     lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V)
-    println("Lower bound obtained by SDDP: "*string(lb_sddp))
+    println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
 end
 
 ######### Solving the problem via Dynamic Programming
 if run_sdp
-    stateSteps = [0.01]
-    controlSteps = [0.01]
+    println("Starting resolution by SDP")
+    stateSteps = [0.01] # discretization step of the state
+    controlSteps = [0.01] # discretization step of the control
     infoStruct = "HD" # noise at time t is known before taking the decision at time t
-
     paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, infoStruct)
     Vs = solve_DP(spmodel,paramSDP, 1)
-    lb_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs)
-    println("Value obtained by SDP: "*string(lb_sdp))
+    value_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs)
+    println("Value obtained by SDP: "*string(round(value_sdp,4)))
+end
+
+if run_ef
+    println("Starting resolution by Extensive Formulation")
+    value_ef = extensive_formulation(spmodel, paramSDDP)[1]
+    println("Value obtained by Extensive Formulation: "*string(round(value_ef,4)))
+    println("Relative error of SDP value: "*string(100*round(value_sdp/value_ef-1,4))*"%")
+    println("Relative error of SDDP lower bound: "*string(100*round(lb_sddp/value_ef-1,4))*"%")
 end
 
-######### Comparing the solution
+######### Comparing the solutions on simulated scenarios.
 scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,1000)
 if run_sddp
     costsddp, stocks = forward_simulations(spmodel, paramSDDP, pbs, scenarios)
@@ -90,5 +100,6 @@ if run_sdp
     costsdp, states, stocks =sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs)
 end
 if run_sddp && run_sdp
-    println(mean(costsddp-costsdp))
+    println("Simulated relative difference between sddp and sdp: "
+            *string(round(200*mean(costsddp-costsdp)/mean(costsddp+costsdp),3))*"%")
 end

src/SDDPoptimize.jl

@@ -39,6 +39,7 @@ fulfilled.
 function solve_SDDP(model::SPModel, param::SDDPparameters, display=0::Int64)
     # initialize value functions:
     V, problems = initialize_value_functions(model, param)
+    println("initial value function intialized")
     # Run SDDP upon example:
     sddp_stats = run_SDDP!(model, param, V, problems, display)
     return V, problems, sddp_stats
@@ -128,10 +129,16 @@ function run_SDDP!(model::SPModel,
         push!(stats.upper_bounds, upb)
 
         if (display > 0) && (iteration_count%display==0)
-            println("Pass number ", iteration_count,
+            if upb < Inf
+                println("Pass number ", iteration_count,
                     "\tUpper-bound: ", upb,
                     "\tLower-bound: ", round(stats.lower_bounds[end], 4),
                     "\tTime: ", round(stats.exectime[end], 2),"s")
+             else 
+                println("Pass number ", iteration_count,
+                    "\tLower-bound: ", round(stats.lower_bounds[end], 4),
+                    "\tTime: ", round(stats.exectime[end], 2),"s")
+             end
         end
 
     end
@@ -313,7 +320,7 @@ function initialize_value_functions(model::SPModel,
                                     param::SDDPparameters)
 
     solverProblems = build_models(model, param)
-
+    println("model builded")
     V = PolyhedralFunction[
                 PolyhedralFunction([], Array{Float64}(0, model.dimStates), 0) for i in 1:model.stageNumber]
 

src/extensiveFormulation.jl

@@ -13,6 +13,10 @@ measurability constraints.
 # Arguments:
 * `model::SPModel`
 * `param::SDDPparameters`
+# Returns
+* `objective value`
+* `first control` 
+* `status of optimization problem`
 """
 function extensive_formulation(model, param)
 
@@ -49,12 +53,16 @@ function extensive_formulation(model, param)
             end
         end
     end
-
-    #Instantiate the problem creating dynamic constraint at each node
-    for t = 1 : (T)
+    #Add state constraints
+    for t = 1 :(T+1)
         for n = 1 : N[t]
             @constraint(mod,[x[t,DIM_STATE*(n-1)+k] for k = 1:DIM_STATE] .>= [model.xlim[k][1] for k = 1:DIM_STATE])
             @constraint(mod,[x[t,DIM_STATE*(n-1)+k] for k = 1:DIM_STATE] .<= [model.xlim[k][2] for k = 1:DIM_STATE])
+        end
+    end
+    #Instantiate the problem creating dynamic constraint at each node
+    for t = 1 : (T)
+        for n = 1 : N[t]
             for xi = 1 : laws[t].supportSize
                 m = (n-1)*laws[t].supportSize+xi
 
@@ -83,12 +91,18 @@ function extensive_formulation(model, param)
     @constraint(mod, [x[1,k] for k = 1:DIM_STATE] .== X_init)
 
     #Define the objective of the function
-    @objective(mod, Min, sum{ sum{proba[t][laws[t].supportSize*(n-1)+k]*c[t,laws[t].supportSize*(n-1)+k],k = 1:laws[t].supportSize} , t = 1:T, n=1:div(N[t+1],laws[t].supportSize)})
+    @objective(mod, Min, 
+    sum{ 
+        sum{    proba[t][laws[t].supportSize*(n-1)+k]*c[t,laws[t].supportSize*(n-1)+k],
+            k = 1:laws[t].supportSize}, 
+        t = 1:T, n=1:div(N[t+1],laws[t].supportSize)}
+    )
 
     status = solve(mod)
     solved = (status == :Optimal)
 
     if solved
+        println("EF value: "*string(getobjectivevalue(mod)))
         firstControl = collect(values(getvalue(u)))[1:DIM_CONTROL*laws[1].supportSize]
         return getobjectivevalue(mod), firstControl, status
     else