StochDynamicProgramming: the stock 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 stock_example()
    model = SDDP.PolicyGraph(
        SDDP.LinearGraph(5);
        lower_bound = -2,
        optimizer = HiGHS.Optimizer,
    ) do sp, stage
        @variable(sp, 0 <= state <= 1, SDDP.State, initial_value = 0.5)
        @variable(sp, 0 <= control <= 0.5)
        @variable(sp, ξ)
        @constraint(sp, state.out == state.in - control + ξ)
        SDDP.parameterize(sp, 0.0:1/30:0.3) do ω
            return JuMP.fix(ξ, ω)
        end
        @stageobjective(sp, (sin(3 * stage) - 1) * control)
    end
    SDDP.train(model; log_frequency = 10)
    @test SDDP.calculate_bound(model) ≈ -1.471 atol = 0.001
    simulation_results = SDDP.simulate(model, 1_000)
    @test length(simulation_results) == 1_000
    μ = SDDP.Statistics.mean(
        sum(data[:stage_objective] for data in simulation) for
        simulation in simulation_results
    )
    @test μ ≈ -1.471 atol = 0.05
    return
end

stock_example()

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
  nodes           : 5
  state variables : 1
  scenarios       : 1.00000e+05
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [5, 5]
  AffExpr in MOI.EqualTo{Float64}         : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [3, 3]
  VariableRef in MOI.LessThan{Float64}    : [2, 3]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [3e-01, 2e+00]
  bounds range     [5e-01, 2e+00]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10  -1.573154e+00 -1.474247e+00  7.209086e-02      1050   1
        20  -1.346690e+00 -1.471483e+00  1.122549e-01      1600   1
        30  -1.308031e+00 -1.471307e+00  1.987529e-01      2650   1
        40  -1.401200e+00 -1.471167e+00  2.418110e-01      3200   1
        50  -1.557483e+00 -1.471097e+00  3.466070e-01      4250   1
        60  -1.534169e+00 -1.471075e+00  3.962181e-01      4800   1
        65  -1.689864e+00 -1.471075e+00  4.204600e-01      5075   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 4.204600e-01
total solves   : 5075
best bound     : -1.471075e+00
simulation ci  : -1.484094e+00 ± 4.058993e-02
numeric issues : 0
-------------------------------------------------------------------