Newsvendor

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

This example is based on the classical newsvendor problem, but features an AR(1) spot-price.

   V(x[t-1], ω[t]) =         max p[t] × u[t]
                      subject to x[t] = x[t-1] - u[t] + ω[t]
                                 u[t] ∈ [0, 1]
                                 x[t] ≥ 0
                                 p[t] = p[t-1] + ϕ[t]

The initial conditions are

x[0] = 2.0
p[0] = 1.5
ω[t] ~ {0, 0.05, 0.10, ..., 0.45, 0.5} with uniform probability.
ϕ[t] ~ {-0.25, -0.125, 0.125, 0.25} with uniform probability.
using SDDP, HiGHS, Statistics, Test

function joint_distribution(; kwargs...)
    names = tuple([first(kw) for kw in kwargs]...)
    values = tuple([last(kw) for kw in kwargs]...)
    output_type = NamedTuple{names,Tuple{eltype.(values)...}}
    distribution = map(output_type, Base.product(values...))
    return distribution[:]
end

function newsvendor_example(; cut_type)
    model = SDDP.PolicyGraph(
        SDDP.LinearGraph(3);
        sense = :Max,
        upper_bound = 50.0,
        optimizer = HiGHS.Optimizer,
    ) do subproblem, stage
        @variables(subproblem, begin
            x >= 0, (SDDP.State, initial_value = 2)
            0 <= u <= 1
            w
        end)
        @constraint(subproblem, x.out == x.in - u + w)
        SDDP.add_objective_state(
            subproblem;
            initial_value = 1.5,
            lower_bound = 0.75,
            upper_bound = 2.25,
            lipschitz = 100.0,
        ) do y, ω
            return y + ω.price_noise
        end
        noise_terms = joint_distribution(;
            demand = 0:0.05:0.5,
            price_noise = [-0.25, -0.125, 0.125, 0.25],
        )
        SDDP.parameterize(subproblem, noise_terms) do ω
            fix(w, ω.demand)
            price = SDDP.objective_state(subproblem)
            @stageobjective(subproblem, price * u)
        end
    end
    SDDP.train(
        model;
        log_frequency = 10,
        time_limit = 20.0,
        cut_type = cut_type,
    )
    @test SDDP.calculate_bound(model) ≈ 4.04 atol = 0.05
    results = SDDP.simulate(model, 500)
    objectives =
        [sum(s[:stage_objective] for s in simulation) for simulation in results]
    @test round(Statistics.mean(objectives); digits = 2) ≈ 4.04 atol = 0.1
    return
end

newsvendor_example(; cut_type = SDDP.SINGLE_CUT)
newsvendor_example(; cut_type = SDDP.MULTI_CUT)
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-26
-------------------------------------------------------------------
problem
  nodes           : 3
  state variables : 1
  scenarios       : 8.51840e+04
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [6, 6]
  AffExpr in MOI.EqualTo{Float64}         : [1, 3]
  AffExpr in MOI.LessThan{Float64}        : [2, 2]
  VariableRef in MOI.GreaterThan{Float64} : [3, 4]
  VariableRef in MOI.LessThan{Float64}    : [3, 3]
numerical stability report
  matrix range     [8e-01, 2e+00]
  objective range  [1e+00, 2e+00]
  bounds range     [1e+00, 1e+02]
  rhs range        [5e+01, 5e+01]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10   5.500000e+00  5.473308e+00  1.803751e-01      1350   1
        20   4.062500e+00  4.451907e+00  2.823040e-01      2700   1
        30   2.993750e+00  4.101372e+00  3.903320e-01      4050   1
        40   5.750000e+00  4.095757e+00  5.059659e-01      5400   1
        50   5.125000e+00  4.093536e+00  6.259680e-01      6750   1
        60   3.737500e+00  4.089225e+00  7.485571e-01      8100   1
        70   4.500000e+00  4.088152e+00  8.732200e-01      9450   1
        80   4.950000e+00  4.087904e+00  1.002296e+00     10800   1
        90   4.125000e+00  4.087756e+00  1.134150e+00     12150   1
       100   3.275000e+00  4.086300e+00  1.275548e+00     13500   1
       110   5.250000e+00  4.086125e+00  1.411344e+00     14850   1
       120   4.975000e+00  4.086054e+00  1.558461e+00     16200   1
       130   3.500000e+00  4.085422e+00  1.700193e+00     17550   1
       140   4.593750e+00  4.085327e+00  1.842921e+00     18900   1
       150   3.000000e+00  4.085258e+00  1.986608e+00     20250   1
       160   2.700000e+00  4.085247e+00  2.136295e+00     21600   1
       170   3.812500e+00  4.085151e+00  2.286552e+00     22950   1
       180   3.875000e+00  4.085121e+00  2.475532e+00     24300   1
       190   4.737500e+00  4.085102e+00  2.626780e+00     25650   1
       200   2.906250e+00  4.085073e+00  2.778281e+00     27000   1
       210   3.750000e+00  4.085050e+00  2.935613e+00     28350   1
       220   5.050000e+00  4.085037e+00  3.091011e+00     29700   1
       230   2.925000e+00  4.085012e+00  3.241679e+00     31050   1
       240   4.500000e+00  4.084970e+00  3.397754e+00     32400   1
       250   4.875000e+00  4.084908e+00  3.558243e+00     33750   1
       260   3.675000e+00  4.084905e+00  3.707877e+00     35100   1
       270   4.725000e+00  4.084903e+00  3.871622e+00     36450   1
       280   3.437500e+00  4.084900e+00  4.031053e+00     37800   1
       290   3.750000e+00  4.084879e+00  4.197158e+00     39150   1
       300   4.125000e+00  4.084879e+00  4.365886e+00     40500   1
       310   4.875000e+00  4.084803e+00  4.532583e+00     41850   1
       320   5.625000e+00  4.084803e+00  4.693489e+00     43200   1
       330   3.825000e+00  4.084800e+00  4.857059e+00     44550   1
       340   3.300000e+00  4.084796e+00  5.023379e+00     45900   1
       350   5.887500e+00  4.084796e+00  5.194910e+00     47250   1
       360   3.600000e+00  4.084786e+00  5.390445e+00     48600   1
       370   3.618750e+00  4.084786e+00  5.554714e+00     49950   1
       380   4.968750e+00  4.084782e+00  5.718562e+00     51300   1
       390   3.300000e+00  4.084780e+00  5.884730e+00     52650   1
       400   4.006250e+00  4.084779e+00  6.061958e+00     54000   1
       410   4.050000e+00  4.084779e+00  6.240460e+00     55350   1
       420   5.125000e+00  4.084776e+00  6.442415e+00     56700   1
       430   4.000000e+00  4.084776e+00  6.617793e+00     58050   1
       440   4.125000e+00  4.084776e+00  6.798742e+00     59400   1
       450   3.112500e+00  4.084776e+00  6.978961e+00     60750   1
       460   3.750000e+00  4.084771e+00  7.159292e+00     62100   1
       470   3.187500e+00  4.084767e+00  7.353450e+00     63450   1
       480   3.031250e+00  4.084757e+00  7.541337e+00     64800   1
       490   4.181250e+00  4.084753e+00  7.726033e+00     66150   1
       500   3.187500e+00  4.084746e+00  7.918957e+00     67500   1
       510   4.875000e+00  4.084741e+00  8.102117e+00     68850   1
       520   3.900000e+00  4.084737e+00  8.282420e+00     70200   1
       530   4.350000e+00  4.084737e+00  8.465570e+00     71550   1
       540   4.212500e+00  4.084737e+00  8.648547e+00     72900   1
       550   4.250000e+00  4.084734e+00  8.834632e+00     74250   1
       560   3.750000e+00  4.084734e+00  9.017368e+00     75600   1
       570   5.306250e+00  4.084730e+00  9.199379e+00     76950   1
       580   4.725000e+00  4.084730e+00  9.378212e+00     78300   1
       590   4.250000e+00  4.084730e+00  9.551538e+00     79650   1
       600   4.000000e+00  4.084730e+00  9.736022e+00     81000   1
       610   4.600000e+00  4.084730e+00  9.904554e+00     82350   1
       620   3.375000e+00  4.084730e+00  1.008167e+01     83700   1
       630   3.981250e+00  4.084725e+00  1.026663e+01     85050   1
       640   3.250000e+00  4.084725e+00  1.044389e+01     86400   1
       650   3.625000e+00  4.084725e+00  1.062162e+01     87750   1
       660   4.781250e+00  4.084725e+00  1.081924e+01     89100   1
       670   4.275000e+00  4.084725e+00  1.100683e+01     90450   1
       680   2.731250e+00  4.084725e+00  1.121048e+01     91800   1
       690   5.237500e+00  4.084725e+00  1.141237e+01     93150   1
       700   3.325000e+00  4.084725e+00  1.160249e+01     94500   1
       710   4.750000e+00  4.084725e+00  1.180494e+01     95850   1
       720   4.537500e+00  4.084725e+00  1.200195e+01     97200   1
       730   4.725000e+00  4.084725e+00  1.220351e+01     98550   1
       740   4.475000e+00  4.084725e+00  1.239690e+01     99900   1
       750   2.893750e+00  4.084725e+00  1.258976e+01    101250   1
       760   3.525000e+00  4.084725e+00  1.278545e+01    102600   1
       770   3.525000e+00  4.084725e+00  1.298323e+01    103950   1
       780   3.262500e+00  4.084725e+00  1.318031e+01    105300   1
       790   4.918750e+00  4.084725e+00  1.338157e+01    106650   1
       800   3.750000e+00  4.084725e+00  1.358311e+01    108000   1
       810   4.687500e+00  4.084725e+00  1.379702e+01    109350   1
       820   4.018750e+00  4.084725e+00  1.400796e+01    110700   1
       830   4.725000e+00  4.084725e+00  1.421774e+01    112050   1
       840   4.268750e+00  4.084725e+00  1.442755e+01    113400   1
       850   5.175000e+00  4.084725e+00  1.462805e+01    114750   1
       860   3.125000e+00  4.084725e+00  1.483835e+01    116100   1
       870   2.762500e+00  4.084725e+00  1.503277e+01    117450   1
       880   3.375000e+00  4.084725e+00  1.524285e+01    118800   1
       890   3.875000e+00  4.084725e+00  1.544789e+01    120150   1
       900   3.093750e+00  4.084725e+00  1.564575e+01    121500   1
       910   4.125000e+00  4.084724e+00  1.586248e+01    122850   1
       920   4.750000e+00  4.084724e+00  1.608119e+01    124200   1
       930   5.431250e+00  4.084724e+00  1.629172e+01    125550   1
       940   3.000000e+00  4.084724e+00  1.650036e+01    126900   1
       950   4.243750e+00  4.084724e+00  1.670391e+01    128250   1
       960   3.756250e+00  4.084724e+00  1.690768e+01    129600   1
       970   3.750000e+00  4.084724e+00  1.712042e+01    130950   1
       980   4.350000e+00  4.084724e+00  1.732700e+01    132300   1
       990   3.375000e+00  4.084722e+00  1.754584e+01    133650   1
      1000   4.750000e+00  4.084722e+00  1.776167e+01    135000   1
      1010   3.825000e+00  4.084722e+00  1.797122e+01    136350   1
      1020   4.000000e+00  4.084722e+00  1.817086e+01    137700   1
      1030   4.500000e+00  4.084722e+00  1.838402e+01    139050   1
      1040   3.637500e+00  4.084722e+00  1.859575e+01    140400   1
      1050   4.050000e+00  4.084722e+00  1.879971e+01    141750   1
      1060   4.925000e+00  4.084722e+00  1.901227e+01    143100   1
      1070   4.500000e+00  4.084722e+00  1.924959e+01    144450   1
      1080   5.000000e+00  4.084722e+00  1.946797e+01    145800   1
      1090   4.393750e+00  4.084722e+00  1.968023e+01    147150   1
      1100   3.750000e+00  4.084722e+00  1.990651e+01    148500   1
      1105   2.531250e+00  4.084722e+00  2.000622e+01    149175   1
-------------------------------------------------------------------
status         : time_limit
total time (s) : 2.000622e+01
total solves   : 149175
best bound     :  4.084722e+00
numeric issues : 0
-------------------------------------------------------------------

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-26
-------------------------------------------------------------------
problem
  nodes           : 3
  state variables : 1
  scenarios       : 8.51840e+04
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [6, 6]
  AffExpr in MOI.EqualTo{Float64}         : [1, 3]
  AffExpr in MOI.LessThan{Float64}        : [2, 2]
  VariableRef in MOI.GreaterThan{Float64} : [3, 4]
  VariableRef in MOI.LessThan{Float64}    : [3, 3]
numerical stability report
  matrix range     [8e-01, 2e+00]
  objective range  [1e+00, 2e+00]
  bounds range     [1e+00, 1e+02]
  rhs range        [5e+01, 5e+01]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10   4.637500e+00  4.048083e+00  2.578578e-01      1350   1
        20   4.475000e+00  4.040445e+00  6.921079e-01      2700   1
        30   3.262500e+00  4.038663e+00  1.345685e+00      4050   1
        40   4.125000e+00  4.038482e+00  2.160043e+00      5400   1
        50   3.518750e+00  4.038076e+00  3.067011e+00      6750   1
        60   4.556250e+00  4.037709e+00  4.099762e+00      8100   1
        70   4.012500e+00  4.037637e+00  5.331666e+00      9450   1
        80   4.468750e+00  4.037465e+00  6.671495e+00     10800   1
        90   3.562500e+00  4.037463e+00  8.097321e+00     12150   1
       100   2.775000e+00  4.037463e+00  9.685515e+00     13500   1
       110   5.000000e+00  4.037463e+00  1.139315e+01     14850   1
       120   3.625000e+00  4.037417e+00  1.329704e+01     16200   1
       130   5.343750e+00  4.037466e+00  1.513318e+01     17550   1
       140   3.037500e+00  4.037417e+00  1.706756e+01     18900   1
       150   3.825000e+00  4.037414e+00  1.929656e+01     20250   1
       154   3.962500e+00  4.037414e+00  2.018312e+01     20790   1
-------------------------------------------------------------------
status         : time_limit
total time (s) : 2.018312e+01
total solves   : 20790
best bound     :  4.037414e+00
numeric issues : 0
-------------------------------------------------------------------