1 Set up

1.1 Problem Presentation

1.2 Setting up the StochDynamicProgramming Library

1.3 Setting up the problem

2 Solving the problem

2.1 Stochastic Dual Dynamic Programming

2.2 Discretized Dynamic Programming

2.3 Extensive Formulation

2.4 Evaluating the solution

3 Extension

3.1 Multi-stock problem

3.2 Display

3.3 Cut Pruning

1.1 Problem Presentation

1.2 Setting up the StochDynamicProgramming Library

1.3 Setting up the problem

2 Solving the problem

2.1 Stochastic Dual Dynamic Programming

2.2 Discretized Dynamic Programming

2.3 Extensive Formulation

2.4 Evaluating the solution

3 Extension

3.1 Multi-stock problem

3.2 Display

3.3 Cut Pruning

Abstract

We consider a simple stock management problem over a ﬁnite horizon {1,…,T_{F }}. The stock
follow the dynamic

where S_{t} is the stock at the beginning of time-step t, U_{t} is the amount of products bought
during [t,t + 1[ for a cost c_{t}, and ξ_{t} is the amount of product used during [t,t + 1[. We assume
that the demand ξ_{t} is a random exogenous noise with known law over a ﬁnite
support. Furthermore we assume a Hazard-Decision setting, where the control U_{t}
is choosen knowing the consumption ξ_{t}. Finally, we assume bound constraint
over the level of stock at each time step and over the control. Thus, the problem
reads

In order to solve the stock management problem using various stochastic optimization approaches we are going to use the StochDynamicProgramming library available in Julia. Julia is an eﬃcient programming language designed for scientiﬁc computing. More information and installation instructions are available at url http://www.julialang.com Installing the library once Julia is installed is easy. Just launch a Julia terminal and type

The ﬁrst line update the list of packages available with their relevant version. The second line add “StochDynamicProgramming” as a required package, and install it as well as all its dependencies packages. All installations are done using git protocol and calling github repositories. In case of trouble you can replace git protocol by https protocol using

In order to solve the extended formulation of our example, as well as in order to use the SDDP algorithm we will require an external linear solver. If you already have CPLEX or Gurobi installed just type

or

Else you can intall a free powerful linear solver with

Once this is done, open (with a text editor) stock-example.jl in the ﬁle examples of your Julia installation

You can run it with julia stock-example.jl from a terminal launched in the same directory

Let’s comment a few lines of the example.

This line tell Julia that you are going to use the library StochDynamicProgramming.jl and Clp.jl. The ﬁrst library is a library to model and solve stochastic optimization problems. The second library is a wrapper to the linear solver Clp. It can be replaced by Gurobi or CPLEX if you have them installed on your computer.

const SOLVER = ClpSolver() # require "using Clp"

const MAX_ITER = 10 # number of iterations of SDDP

const N_STAGES = 6 # number of stages of the SP problem

const COSTS = [sin(3*t)-1 for t in 1:N_STAGES]

const CONTROL_MAX = 0.5 # bounds on the control

const CONTROL_MIN = 0

const XI_MAX = 0.3 # bounds on the noise

const XI_MIN = 0

const N_XI = 10 # discretization of the noise

const S0 = 0.5 # initial stock

const MAX_ITER = 10 # number of iterations of SDDP

const N_STAGES = 6 # number of stages of the SP problem

const COSTS = [sin(3*t)-1 for t in 1:N_STAGES]

const CONTROL_MAX = 0.5 # bounds on the control

const CONTROL_MIN = 0

const XI_MAX = 0.3 # bounds on the noise

const XI_MIN = 0

const N_XI = 10 # discretization of the noise

const S0 = 0.5 # initial stock

These lines deﬁne constant parameters of the problem. You can modify them to test variations of the same problem.

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 = NoiseLaw[xi_law for t in 1:N_STAGES-1]

xi_support = collect(linspace(XI_MIN,XI_MAX,N_XI))

xi_law = NoiseLaw(xi_support, proba)

xi_laws = NoiseLaw[xi_law for t in 1:N_STAGES-1]

The ﬁrst line deﬁne the probabilities of ξ_{t}, the second the support of ξ_{t}. The third line create the
actual object representing the random variable. The fourth line state that create the vector
(ξ_{t})_{t}. The noises are assumed to be independent.

This deﬁne the function linking s_{t+1} and s_{t} (see constraint (1b)).

This function is the current cost to minimize.

s_bounds = [(0, 1)] # bounds on the state

u_bounds = [(CONTROL_MIN, CONTROL_MAX)] # bounds on controls

spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws)

set_state_bounds(spmodel, s_bounds) # adding the bounds to the model

u_bounds = [(CONTROL_MIN, CONTROL_MAX)] # bounds on controls

spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws)

set_state_bounds(spmodel, s_bounds) # adding the bounds to the model

Finally, these lines deﬁne the spmodel object, which represent our stochastic model including, dynamics, initial state, costs, bounds, and laws of the noises.

Now that the problem is set-up it can be solved through the diﬀerent methods.

paramSDDP = SDDPparameters(SOLVER,

passnumber=10,

max_iterations=MAX_ITER)

V, pbs = solve_SDDP(spmodel, paramSDDP, 2)

passnumber=10,

max_iterations=MAX_ITER)

V, pbs = solve_SDDP(spmodel, paramSDDP, 2)

The ﬁrst line deﬁne the parameters of the algorithm SDDP. The ﬁrst argument specify the solver required, the second argument the number of forward pass done at each iterations and the third the number of iteration of the SDDP algorithm before stopping. The second line actually solve the problem.

Question 2 In the ﬁle stock-example.jl, set run_sdp’, ’run_ef’ and ’test_simulation’ to ’false’ ’.

Modify the example such that:

- SDDP run for 20 iterations
- SDDP make 5 forward pass

Then run the example either by

from a bash terminal, or by typing in a julia terminal

Remark: in order to navigate from a Julia terminal just type ’;’ and you can use any bash command.

Compare computation time and bounds obtained for multiple values of number of iterations and forward pass.

Question 3 Set ’run_sdp’ to ’true’. Run the example to compare the lower bound of sddp and the approximated value of the problem given by sdp.

Test diﬀerent discretization step size (e.g. [0.2,0.1,0.05,0.01]), and compare the computation time and quality of approximation. In theory, how should evolve the computation time with the discretization step ? Do you observe it in practice ?

Question 4 We have solved the problem in Hazard-Decision framework. If the problem was in Decision-Hazard framework, would the value be lower or higher ? Modify the code to check. Don’t forget to revert to hazard-decision before continuing.

Question 5 Can you determine inequalities between the lower bound of sddp, the value computed by sdp and the value computed by extensive formulation ?

Set ’run_ef’ to ’true’. Numerically check your intuition.

Question 6 Theoretically, how should evolve the computing time of the three methods with N_XI and N_STAGES ? Numerically test your intuition by increasing / decreasing the parameters N_XI and N_STAGES.

Once we have obtained the solution to our problem we would like to simulate the strategy obtained and evaluate the costs.

Question 7 Simulate the strategies obtained by sddp over 1000 scenarios. Print the mean cost obtained by each strategy, and the relative diﬀerence. Compare to the lower bound of sddp, the value given by sdp and the value given by extensive formulation. Repeat 10 times and comment.

Question 8 Fix the random seed, compare the quality of the solution obtained with various parameters of sddp and sdp on the same scenarios.

We present here a few extensions to manipulate the package.

The example ’multistock-example.jl’ is extremely close to ’stock-example.jl’ but incorporate multiple stocks linked through a constraint on the controls.

Question 9 Theoretically what is the computing time dependence in the parameter ’step’ ? Test it by setting step to 0.05 and 0.2. Comment on the quality of the solution.

Question 10 Theoretically what is the computing time dependence in the parameter N_STAGES ? Test it.

Question 11 Theoretically what is the computing time dependence in the parameter N_STOCKS ? Test it.

Question 12 Theoretically what is the computing time dependence in the parameter N_XI ? Test it.

In a Julia terminal in the example directory.

Question 13 Plot the cost distribution of SDDP and SDP.

Question 14 Plot some optimal trajectories of stocks.

A cut pruning method is available in the package and can be used as follows.

paramSDDP_pruned = SDDPparameters(SOLVER,

passnumber=10,

max_iterations=MAX_ITER,

pruning_algo="exact",

prune_cuts = 1)

passnumber=10,

max_iterations=MAX_ITER,

pruning_algo="exact",

prune_cuts = 1)

The parameter pruning_algo can be set to the string value "exact" or to to the string value "level1". The prune_cuts parameter is used to deﬁne at which periodicity the pruning method must be applied.

Question 15 test the pruning on the stock example and make comments on the computing time eﬃciency.