StochDynamicProgramming: the multistock problem
This tutorial was generated using Literate.jl. Download the source as a .jl
file. Download the source as a .ipynb
file.
This example comes from StochDynamicProgramming.jl.
using SDDP, HiGHS, Test
function test_multistock_example()
model = SDDP.LinearPolicyGraph(;
stages = 5,
lower_bound = -5.0,
optimizer = HiGHS.Optimizer,
) do subproblem, stage
@variable(
subproblem,
0 <= stock[i = 1:3] <= 1,
SDDP.State,
initial_value = 0.5
)
@variables(subproblem, begin
0 <= control[i = 1:3] <= 0.5
ξ[i = 1:3] # Dummy for RHS noise.
end)
@constraints(
subproblem,
begin
sum(control) - 0.5 * 3 <= 0
[i = 1:3], stock[i].out == stock[i].in + control[i] - ξ[i]
end
)
Ξ = collect(
Base.product((0.0, 0.15, 0.3), (0.0, 0.15, 0.3), (0.0, 0.15, 0.3)),
)[:]
SDDP.parameterize(subproblem, Ξ) do ω
return JuMP.fix.(ξ, ω)
end
@stageobjective(subproblem, (sin(3 * stage) - 1) * sum(control))
end
SDDP.train(
model;
iteration_limit = 100,
cut_type = SDDP.SINGLE_CUT,
log_frequency = 10,
)
@test SDDP.calculate_bound(model) ≈ -4.349 atol = 0.01
simulation_results = SDDP.simulate(model, 5000)
@test length(simulation_results) == 5000
μ = SDDP.Statistics.mean(
sum(data[:stage_objective] for data in simulation) for
simulation in simulation_results
)
@test μ ≈ -4.349 atol = 0.1
return
end
test_multistock_example()
-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
nodes : 5
state variables : 3
scenarios : 1.43489e+07
existing cuts : false
options
solver : serial mode
risk measure : SDDP.Expectation()
sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
VariableRef : [13, 13]
AffExpr in MOI.EqualTo{Float64} : [3, 3]
AffExpr in MOI.LessThan{Float64} : [1, 1]
VariableRef in MOI.GreaterThan{Float64} : [7, 7]
VariableRef in MOI.LessThan{Float64} : [6, 7]
numerical stability report
matrix range [1e+00, 1e+00]
objective range [3e-01, 2e+00]
bounds range [5e-01, 5e+00]
rhs range [2e+00, 2e+00]
-------------------------------------------------------------------
iteration simulation bound time (s) solves pid
-------------------------------------------------------------------
10 -3.878303e+00 -4.434982e+00 1.888111e-01 1400 1
20 -4.262885e+00 -4.399265e+00 3.079290e-01 2800 1
30 -3.075162e+00 -4.382527e+00 4.321101e-01 4200 1
40 -3.761147e+00 -4.369587e+00 6.219380e-01 5600 1
50 -4.323162e+00 -4.362199e+00 7.603130e-01 7000 1
60 -3.654943e+00 -4.358401e+00 9.001839e-01 8400 1
70 -4.010883e+00 -4.357368e+00 1.041047e+00 9800 1
80 -4.314412e+00 -4.355714e+00 1.186634e+00 11200 1
90 -4.542422e+00 -4.353708e+00 1.337672e+00 12600 1
100 -4.178952e+00 -4.351685e+00 1.482508e+00 14000 1
-------------------------------------------------------------------
status : iteration_limit
total time (s) : 1.482508e+00
total solves : 14000
best bound : -4.351685e+00
simulation ci : -4.246786e+00 ± 8.703997e-02
numeric issues : 0
-------------------------------------------------------------------