Merge pull request #122 from JuliaOpt/dev-nightly

v0.1.4

Author: leclere

doc/quickstart.rst

@@ -1,7 +1,7 @@
 .. _quickstart:
 
 ====================
-Step-by-step example
+SDDP: Step-by-step example
 ====================
 
 This page gives a short introduction to the interface of this package. It explains the resolution with SDDP of a classical example: the management of a dam over one year with random inflow.
@@ -52,7 +52,7 @@ Dynamic
 We write the dynamic (which return a vector)::
 
     function dynamic(t, x, u, xi)
-        return [x[1] + u[1] - xi[1]] 
+        return [x[1] + u[1] - xi[1]]
     end
 
 
@@ -105,6 +105,10 @@ As our problem is purely linear, we instantiate::
 
     spmodel = LinearDynamicLinearCostSPmodel(N_STAGES,u_bounds,X0,cost_t,dynamic,xi_laws)
 
+We add the state bounds to the model afterward::
+
+    set_state_bounds(spmodel, s_bounds)
+
 
 Solver
 ^^^^^^

doc/quickstart_sdp.rst

@@ -0,0 +1,164 @@
+.. _quickstart_sdp:
+
+====================
+SDP: Step-by-step example
+====================
+
+This page gives a short introduction to the interface of this package. It explains the resolution with Stochastic Dynamic Programming of a classical example: the management of a dam over one year with random inflow.
+
+Use case
+========
+In the following, :math:`x_t` will denote the state and :math:`u_t` the control at time :math:`t`.
+We will consider a dam, whose dynamic is:
+
+.. math::
+   x_{t+1} = x_t - u_t + w_t
+
+At time :math:`t`, we have a random inflow :math:`w_t` and we choose to turbine a quantity :math:`u_t` of water.
+
+The turbined water is used to produce electricity, which is being sold at a price :math:`c_t`. At time :math:`t` we gain:
+
+.. math::
+    C(x_t, u_t, w_t) = c_t \times u_t
+
+We want to minimize the following criterion:
+
+.. math::
+    J = \underset{x, u}{\min} \sum_{t=0}^{T-1} C(x_t, u_t, w_t)
+
+We will assume that both states and controls are bounded:
+
+.. math::
+    x_t \in [0, 100], \qquad u_t \in [0, 7]
+
+
+Problem definition in Julia
+===========================
+
+We will consider 52 time steps as we want to find optimal value functions for one year::
+
+    N_STAGES = 52
+
+
+and we consider the following initial position::
+
+    X0 = [50]
+
+Note that X0 is a vector.
+
+Dynamic
+^^^^^^^
+
+We write the dynamic (which return a vector)::
+
+    function dynamic(t, x, u, xi)
+        return [x[1] + u[1] - xi[1]]
+    end
+
+
+Cost
+^^^^
+
+we store evolution of costs :math:`c_t` in an array `COSTS`, and we define the cost function (which return a float)::
+
+    function cost_t(t, x, u, w)
+        return COSTS[t] * u[1]
+    end
+
+Noises
+^^^^^^
+
+Noises are defined in an array of Noiselaw. This type defines a discrete probability.
+
+
+For instance, if we want to define a uniform probability with size :math:`N= 10`, such that:
+
+.. math::
+    \mathbb{P} \left(X_i = i \right) = \dfrac{1}{N} \qquad \forall i \in 1 .. N
+
+we write::
+
+    N = 10
+    proba = 1/N*ones(N) # uniform probabilities
+    xi_support = collect(linspace(1,N,N))
+    xi_law = NoiseLaw(xi_support, proba)
+
+
+Thus, we could define a different probability laws for each time :math:`t`. Here, we suppose that the probability is constant over time, so we could build the following vector::
+
+    xi_laws = NoiseLaw[xi_law for t in 1:N_STAGES-1]
+
+
+Bounds
+^^^^^^
+
+We add bounds over the state and the control::
+
+    s_bounds = [(0, 100)]
+    u_bounds = [(0, 7)]
+
+
+Problem definition
+^^^^^^^^^^^^^^^^^^
+
+We have two options to contruct a model that can be solved by the SDP algorithm.
+We can instantiate a model that can be solved by SDDP as well::
+
+    spmodel = LinearDynamicLinearCostSPmodel(N_STAGES,u_bounds,X0,cost_t,dynamic,xi_laws)
+
+    set_state_bounds(spmodel, s_bounds)
+
+Or we can instantiate a StochDynProgModel that can be solved only by SDP but we
+need to define the constraint function and the final cost function::
+
+    function constraints(t, x, u, xi) # return true when there is no constraints ecept state and control bounds
+        return true
+    end
+
+    function final_cost_function(x)
+        return 0
+    end
+
+    spmodel = StochDynProgModel(N_STAGES, s_bounds, u_bounds, X0, cost_t,
+                                final_cost_function, dynamic, constraints,
+                                xi_laws)
+
+
+Solver
+^^^^^^
+
+It remains to define SDP algorithm parameters::
+
+    stateSteps = [1] # discretization steps of the state space
+    controlSteps = [0.1] # discretization steps of the control space
+    infoStruct = "HD" # noise at time t is known before taking the decision at time t
+    paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, infoStruct)
+
+
+Now, we solve the problem by computing Bellman values::
+
+    Vs = solve_DP(spmodel,paramSDP, 1)
+
+:code:`V` is an array storing the value functions
+
+We have an exact lower bound given by :code:`V` with the function::
+
+    value_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs)
+
+
+Find optimal controls
+=====================
+
+Once Bellman functions are computed, we can control our system over assessments scenarios, without assuming knowledge of the future.
+
+We build 1000 scenarios according to the laws stored in :code:`xi_laws`::
+
+    scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,1000)
+
+We compute 1000 simulations of the system over these scenarios::
+
+    costsdp, states, controls =sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs)
+
+:code:`costsdp` returns the costs for each scenario, :code:`states` the simulation of each state variable along time, for each scenario, and
+:code:`controls` returns the optimal controls for each scenario
+

examples/battery_storage_parallel.jl

@@ -77,7 +77,7 @@ println("library loaded")
     end
 
     function constraint(t, x, u, xi)
-    	return( (x[1] <= s_bounds[1][2] )&(x[1] >= s_bounds[1][1]))
+    	return true
     end
 
     function finalCostFunction(x)

examples/benchmark.jl

@@ -167,20 +167,23 @@ function init_problem()
     x_bounds = [(VOLUME_MIN, VOLUME_MAX), (VOLUME_MIN, VOLUME_MAX)]
     u_bounds = [(CONTROL_MIN, CONTROL_MAX), (CONTROL_MIN, CONTROL_MAX), (0, Inf), (0, Inf)]
 
-    model = LinearDynamicLinearCostSPmodel(N_STAGES,
-                                                u_bounds,
-                                                x0,
-                                                cost_t,
-                                                dynamic,
-                                                aleas)
+    model = LinearSPModel(N_STAGES,
+                            u_bounds,
+                            x0,
+                            cost_t,
+                            dynamic,
+                            aleas)
 
     set_state_bounds(model, x_bounds)
 
 
     EPSILON = .05
     MAX_ITER = 20
     solver = SOLVER
-    params = SDDPparameters(solver, N_SCENARIOS, EPSILON, MAX_ITER)
+    params = SDDPparameters(solver,
+                            passnumber=N_SCENARIOS,
+                            gap=EPSILON,
+                            max_iterations=MAX_ITER)
 
     return model, params
 end
@@ -259,7 +262,7 @@ function benchmark_sdp()
     end
 
     function constraints(t, x, u, w)
-        return (VOLUME_MIN<=x[1]<=VOLUME_MAX)&(VOLUME_MIN<=x[2]<=VOLUME_MAX)
+        return true
     end
 
     function finalCostFunction(x)

examples/dam.jl

@@ -142,16 +142,19 @@ function init_problem()
 
     x_bounds = [(0, 100)]
     u_bounds = [(0, 7), (0, 7)]
-    model = StochDynamicProgramming.LinearDynamicLinearCostSPmodel(N_STAGES,
-                                                u_bounds,
-                                                x0,
-                                                cost_t,
-                                                dynamic, aleas)
+    model = StochDynamicProgramming.LinearSPModel(N_STAGES,
+                                                  u_bounds,
+                                                  x0,
+                                                  cost_t,
+                                                  dynamic, aleas)
 
     set_state_bounds(model, x_bounds)
 
     solver = SOLVER
-    params = StochDynamicProgramming.SDDPparameters(solver, N_SCENARIOS, EPSILON, MAX_ITER)
+    params = StochDynamicProgramming.SDDPparameters(solver,
+                                                    passnumber=N_SCENARIOS,
+                                                    gap=EPSILON,
+                                                    max_iterations=MAX_ITER)
 
     return model, params
 end

examples/damsvalley.jl

@@ -84,7 +84,7 @@ end
 const FORWARD_PASS = 10.
 const EPSILON = .05
 # Maximum number of iterations
-const MAX_ITER = 50
+const MAX_ITER = 10
 ##################################################
 
 """Build probability distribution at each timestep.
@@ -107,25 +107,25 @@ function init_problem()
 
     x_bounds = [(VOLUME_MIN, VOLUME_MAX) for i in 1:N_DAMS]
     u_bounds = vcat([(CONTROL_MIN, CONTROL_MAX) for i in 1:N_DAMS], [(0., 200) for i in 1:N_DAMS]);
-    model = LinearDynamicLinearCostSPmodel(N_STAGES,
-                                                u_bounds,
-                                                X0,
-                                                cost_t,
-                                                dynamic,
-                                                aleas,
-                                                final_cost_dams)
+    model = LinearSPModel(N_STAGES, u_bounds,
+                          X0, cost_t,
+                          dynamic, aleas,
+                          Vfinal=final_cost_dams)
 
     # Add bounds for stocks:
     set_state_bounds(model, x_bounds)
 
     # We need to use CPLEX to solve QP at final stages:
     solver = CPLEX.CplexSolver(CPX_PARAM_SIMDISPLAY=0, CPX_PARAM_BARDISPLAY=0)
-    params = SDDPparameters(solver, FORWARD_PASS, EPSILON, MAX_ITER)
 
+    params = SDDPparameters(solver,
+                            passnumber=FORWARD_PASS,
+                            gap=EPSILON,
+                            max_iterations=MAX_ITER)
     return model, params
 end
 
 # Solve the problem:
 model, params = init_problem()
-V, pbs = solve_SDDP(model, params, 1)
+V, pbs = @time solve_SDDP(model, params, 1)
 

examples/parallel_sddp.jl

@@ -0,0 +1,58 @@
+
+import StochDynamicProgramming
+
+
+"""
+Solve SDDP in parallel, dispatching both forward and backward passes to process, 
+which is not the most standard parallelization of SDDP.
+
+# Arguments
+* `model::SPmodel`:
+    the stochastic problem we want to optimize
+* `param::SDDPparameters`:
+    the parameters of the SDDP algorithm
+* `V::Array{PolyhedralFunction}`:
+    the current estimation of Bellman's functions
+* `n_parallel_pass::Int`: default is 4
+    Number of parallel pass to compute
+* `synchronize::Int`: default is 5
+    Synchronize the cuts between the different processes every "synchronise" iterations
+* `display::Int`: default is 0
+    Says whether to display results or not
+
+# Return
+* `V::Array{PolyhedralFunction}`:
+    the collection of approximation of the bellman functions
+"""
+function psolve_sddp(model, params, V; n_parallel_pass=4,
+                     synchronize=5, display=0)
+    # Redefine seeds in every processes to maximize randomness:
+    @everywhere srand()
+
+    mitn = params.maxItNumber
+    params.maxItNumber = synchronize
+
+    # Count number of available CPU:
+    ncpu = nprocs() - 1
+    (display > 0) && println("\nLaunch simulation on ", ncpu, " processes")
+    workers = procs()[2:end]
+
+    fpn = params.forwardPassNumber
+    # As we distribute computation in n process, we perform forward pass in parallel:
+    params.forwardPassNumber = max(1, round(Int, params.forwardPassNumber/ncpu))
+
+    # Start parallel computation:
+    for i in 1:n_parallel_pass
+        # Distribute computation of SDDP in each process:
+        refs = [@spawnat w StochDynamicProgramming.solve_SDDP(model, params, V, display)[1] for w in workers]
+        # Catch the result in the main process:
+        V = StochDynamicProgramming.catcutsarray([fetch(r) for r in refs]...)
+        # We clean the resultant cuts:
+        StochDynamicProgramming.remove_redundant_cuts!(V)
+        (display > 0) && println("Lower bound at pass ", i, ": ", StochDynamicProgramming.get_lower_bound(model, params, V))
+    end
+
+    params.forwardPassNumber = fpn
+    params.maxItNumber = mitn
+    return V
+end

examples/stock-example.jl

@@ -58,14 +58,18 @@ end
 ######## Setting up the SPmodel
 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)
+spmodel = LinearSPModel(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
     println("Starting resolution by SDDP")
-    paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward pass, stop at MAX_ITER
+    # 10 forward pass, stop at MAX_ITER
+    paramSDDP = StochDynamicProgramming.SDDPparameters(SOLVER,
+                                                    passnumber=10,
+                                                    gap=0,
+                                                    max_iterations=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(round(lb_sddp,4)))