Partially observable inventory management

This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.

using SDDP, HiGHS, Random, Statistics, Testfunction inventory_management_problem()    demand_values = [1.0, 2.0]    demand_prob = Dict(:Ah => [0.2, 0.8], :Bh => [0.8, 0.2])    graph = SDDP.Graph(        :root_node,        [:Ad, :Ah, :Bd, :Bh],        [            (:root_node => :Ad, 0.5),            (:root_node => :Bd, 0.5),            (:Ad => :Ah, 1.0),            (:Ah => :Ad, 0.8),            (:Ah => :Bd, 0.1),            (:Bd => :Bh, 1.0),            (:Bh => :Bd, 0.8),            (:Bh => :Ad, 0.1),        ],    )    SDDP.add_ambiguity_set(graph, [:Ad, :Bd], 1e2)    SDDP.add_ambiguity_set(graph, [:Ah, :Bh], 1e2)    model = SDDP.PolicyGraph(        graph;        lower_bound = 0.0,        optimizer = HiGHS.Optimizer,    ) do subproblem, node        @variables(            subproblem,            begin                0 <= inventory <= 2, (SDDP.State, initial_value = 0.0)                buy >= 0                demand            end        )        @constraint(subproblem, demand == inventory.in - inventory.out + buy)        if node == :Ad || node == :Bd || node == :D            fix(demand, 0)            @stageobjective(subproblem, buy)        else            SDDP.parameterize(subproblem, demand_values, demand_prob[node]) do ω                fix(demand, ω)                return            end            @stageobjective(subproblem, 2 * buy + inventory.out)        end    end    # Train the policy.    Random.seed!(123)    SDDP.train(        model;        iteration_limit = 100,        cut_type = SDDP.SINGLE_CUT,        log_frequency = 10,        parallel_scheme = SDDP.Serial(),    )    results = SDDP.simulate(model, 500; parallel_scheme = SDDP.Serial())    objectives =        [sum(s[:stage_objective] for s in simulation) for simulation in results]    sample_mean = round(Statistics.mean(objectives); digits = 2)    sample_ci = round(1.96 * Statistics.std(objectives) / sqrt(500); digits = 2)    @test SDDP.calculate_bound(model)  sample_mean atol = sample_ci    returnendinventory_management_problem()
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-26
-------------------------------------------------------------------
problem
  nodes           : 4
  state variables : 1
  scenarios       : Inf
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [7, 7]
  AffExpr in MOI.EqualTo{Float64}         : [1, 1]
  AffExpr in MOI.GreaterThan{Float64}     : [2, 2]
  VariableRef in MOI.EqualTo{Float64}     : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [5, 5]
  VariableRef in MOI.LessThan{Float64}    : [3, 3]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [1e+00, 2e+00]
  bounds range     [2e+00, 1e+02]
  rhs range        [0e+00, 0e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10   2.694235e+00  1.060052e+01  8.913491e-01       900   1
        20   5.498088e+00  1.456766e+01  9.949920e-01      1720   1
        30   2.912197e+01  1.665921e+01  1.220594e+00      3036   1
        40   1.608515e+01  1.792397e+01  1.498812e+00      4192   1
        50   3.997965e+00  1.830624e+01  1.706787e+00      5020   1
        60   1.102045e+01  1.840685e+01  1.938382e+00      5808   1
        70   8.045596e+00  1.846781e+01  2.170205e+00      6540   1
        80   4.710242e+01  1.851858e+01  2.348811e+00      7088   1
        90   3.901806e+01  1.865685e+01  2.765516e+00      8180   1
       100   8.003646e+00  1.869095e+01  2.946005e+00      8664   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 2.946005e+00
total solves   : 8664
best bound     :  1.869095e+01
numeric issues : 0
-------------------------------------------------------------------