We consider a simple stock management problem over a finite horizon {1,…,TF}. The stockfollow the dynamic
whereStis the stock at the beginning of time-step t,Utis the amount of products boughtduring [t,t + 1[ for a cost ct, andξtis the amount of product used during [t,t + 1[. We assumethat the demandξtis a random exogenous noise with known law over a finitesupport. Furthermore we assume a Hazard-Decision setting, where the controlUtis choosen knowing the consumptionξt. Finally, we assume bound constraintover the level of stock at each time step and over the control. Thus, the problemreads
1.2 Setting up the StochDynamicProgramming Library
In order to solve the stock management problem using various stochastic optimizationapproaches we are going to use the StochDynamicProgramming library available in Julia.Julia is an efficient programming language designed for scientific computing. Moreinformation and installation instructions are available at url http://www.julialang.comInstalling the library once Julia is installed is easy. Just launch a Julia terminal andtype
The first line update the list of packages available with their relevant version. The second lineadd “StochDynamicProgramming” as a required package, and install it as well as all itsdependencies packages. All installations are done using git protocol and calling githubrepositories. In case of trouble you can replace git protocol by https protocolusing
In order to solve the extended formulation of our example, as well as in order to use theSDDP algorithm we will require an external linear solver. If you already have CPLEX orGurobi installed just type
julia> Pkg.add("CPLEX")
or
julia> Pkg.add("Gurobi")
Else you can intall a free powerful linear solver with
julia> Pkg.add("Clp")
1.3 Setting up the problem
Once this is done, open (with a text editor) stock-example.jl in the file examples of your Juliainstallation
~/.julia/v0.5/StochDynamicProgramming/examples
You can run it with julia stock-example.jl from a terminal launched in the samedirectory
Let’s comment a few lines of the example.
julia> using StochDynamicProgramming, Clp
This line tell Julia that you are going to use the library StochDynamicProgramming.jl and Clp.jl.The first library is a library to model and solve stochastic optimization problems. The secondlibrary is a wrapper to the linear solver Clp. It can be replaced by Gurobi or CPLEX if youhave 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
These lines define constant parameters of the problem. You can modify them to test variations ofthe 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]
The first line define the probabilities ofξt, the second the support ofξt. The third line create theactual object representing the random variable. The fourth line state that create the vector(ξt)t. The noises are assumed to be independent.
function dynamic(t, x, u, xi) return [x[1] + u[1] - xi[1]] end
This define the function linking st+1and st(see constraint(1b)).
function cost_t(t, x, u, w) return COSTS[t] * u[1] end
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
Finally, these lines define the spmodel object, which represent our stochastic model including,dynamics, initial state, costs, bounds, and laws of the noises.
2 Solving the problem
Now that the problem is set-up it can be solved through the different methods.
The first line define the parameters of the algorithm SDDP. The first argument specify the solverrequired, the second argument the number of forward pass done at each iterations and thethird the number of iteration of the SDDP algorithm before stopping. The second lineactually solve the problem.
Question 1Open a Julia terminal. Type
Julia> using StochDynamicProgramming
(auto completion should be working). Then type ’?’ to access the help, and type
help?> solve_SDDP
to obtain a description of the function used here.
Question 2In the file 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
julia stock-example
from a bash terminal, or by typing in a julia terminal
julia> include("stock-example.jl")
Remark: in order to navigate from a Julia terminal just type ’;’ and you can use any bashcommand.
Compare computation time and bounds obtained for multiple values of number of iterationsand forward pass.
2.2 Discretized Dynamic Programming
Question 3Set ’run_sdp’ to ’true’. Run the example to compare the lower bound ofsddp and the approximated value of the problem given by sdp.
Test different discretization step size (e.g. [0.2,0.1,0.05,0.01]), and compare thecomputation time and quality of approximation. In theory, how should evolve thecomputation time with the discretization step ? Do you observe it in practice ?
Question 4We have solved the problem in Hazard-Decision framework. If theproblem was in Decision-Hazard framework, would the value be lower or higher ? Modifythe code to check. Don’t forget to revert to hazard-decision before continuing.
2.3 Extensive Formulation
Question 5Can you determine inequalities between the lower bound of sddp, thevalue computed by sdp and the value computed by extensive formulation ?
Set ’run_ef’ to ’true’. Numerically check your intuition.
Question 6Theoretically, how should evolve the computing time of the three methodswith N_XIand N_STAGES? Numerically test your intuition by increasing / decreasingthe parameters N_XIand N_STAGES.
2.4 Evaluating the solution
Once we have obtained the solution to our problem we would like to simulate the strategyobtained and evaluate the costs.
Question 7Simulate the strategies obtained by sddp over 1000 scenarios. Print themean cost obtained by each strategy, and the relative difference. Compare to the lowerbound of sddp, the value given by sdp and the value given by extensive formulation.Repeat 10 times and comment.
Question 8Fix the random seed, compare the quality of the solution obtained withvarious parameters of sddp and sdp on the same scenarios.
3 Extension
We present here a few extensions to manipulate the package.
3.1 Multi-stock problem
The example ’multistock-example.jl’ is extremely close to ’stock-example.jl’ but incorporatemultiple stocks linked through a constraint on the controls.
Question 9Theoretically 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 10Theoretically what is the computing time dependence in the parameterN_STAGES? Test it.
Question 11Theoretically what is the computing time dependence in the parameterN_STOCKS? Test it.
Question 12Theoretically what is the computing time dependence in the parameterN_XI? Test it.
The parameter pruning_algocan be set to the string value "exact"or to to the stringvalue "level1". The prune_cutsparameter is used to define at which periodicity the pruningmethod must be applied.
Question 15test the pruning on the stock example and make comments on thecomputing time efficiency.