Hi all. I’m seeking help to diagnose a simple inventory problem over a finite horizon I’ve set up with SDDP. Currently, my model is such that the stochasticity is demand, for which we assume we know the intensity of the Poisson distributed demand at each time point. The decision variable is whether or not to place an order, and if so, the order size; orders arrive after a lead time. We may then dispense stock from inventory to fulfill demand to avoid incurring stockout costs, or keep inventory, incurring a mild holding cost. At the end of the time horizon, leftover stock must be disponsed with a cost.
I’m running into occasional issues in training or (strangely enough), training will complete succesfully but out of sample MC simulation using the same distributional assumptions will fail. The numerical stability report does not seem too bad, see:
julia> SDDP.numerical_stability_report(model)
numerical stability report
matrix range [1e+00, 1e+03]
objective range [1e+00, 2e+02]
bounds range [0e+00, 0e+00]
rhs range [0e+00, 0e+00]
The errors I get (which, again, are sometimes during training, sometimes during simulation if the training completes its alloted iterations) are of the form:
caused by: Unable to retrieve solution from node 22.
Termination status : OTHER_ERROR
Primal status : NO_SOLUTION
Dual status : NO_SOLUTION.
The current subproblem was written to `subproblem_22.mof.json`.
I’ve attached the model below. I do not think I am violating relatively complete recourse. If anyone has any suggestions for what may be causing these occasional errors, I would be grateful for help. Thanks. I’d like to understand what is going wrong in this simple case before in the future I end up adding complications such as stochastic lead times, and stochastic inventory perishing.
using SDDP, HiGHS, Distributions, Plots
time_horizon = 50
# ---------- parameters ----------
initial_inventory = 0
shipping_delay = 3
storage_cost = 1 # storage cost per unit per day
order_cost = 10 # cost to place an order (regardless of size)
order_size_cost = 2 # cost per unit to order
wastage_cost = 50 # cost per unit of disposing of excess inventory at end of time horizon
stockout_cost = 200 # cost per unit of demand that cannot be fulfilled
bigM = 1_000 # a number bigger than any reasonable order size (order sizes will be capped at this val)
# demand (patient arrivals) is a Poisson process, but with a known intensity function, given here
# it linearly increases to a peak, then linearly decreases to zero
demand_λ = zeros(time_horizon)
demand_midpoint = Int(floor((time_horizon-shipping_delay)/2))
demand_λ[shipping_delay+2:demand_midpoint] = range(0.1, 10, length(shipping_delay+2:demand_midpoint))
demand_λ[demand_midpoint+1:end] = range(9.9, 0, length(demand_midpoint+1:time_horizon))
p = plot(demand_λ, title="demand intensity", legend=false)
for _ in 1:10
plot!(p, rand.(Poisson.(demand_λ)), legend=false, linealpha=0.2, color="red")
end
p
model = SDDP.LinearPolicyGraph(
stages = time_horizon,
sense = :Min,
lower_bound = 0,
optimizer = HiGHS.Optimizer,
) do sp, t
# state variables
@variables(sp, begin
x_inventory >= 0, SDDP.State, (initial_value=initial_inventory)
x_transit[1:shipping_delay] >= 0, SDDP.State, (initial_value=0)
end)
# decision variables
@variables(sp, begin
u_order, Bin # order or not
u_order_size >= 0 # order size
u_dispense >= 0
end)
# helper variables
@variables(sp, begin
h_lost >= 0
end)
# random variables (demand)
@variable(sp, ω_demand)
Ω = range(quantile.(Poisson(demand_λ[t]), [0.05, 0.95])...) # get most of the probability mass
P = pdf.(Poisson(demand_λ[t]), Ω)
P = P ./ sum(P)
SDDP.parameterize(sp, Ω, P) do ω
return JuMP.fix(ω_demand, ω)
end
# constraints for ordering and shipping
@constraint(sp, u_order_size <= bigM * u_order)
@constraint(sp, x_transit[1].out == u_order_size)
@constraint(sp, c_transit[k=2:shipping_delay], x_transit[k].out == x_transit[k-1].in)
# constraints for inventory
@constraint(sp, u_dispense <= ω_demand) # dispensing more than demand does not make sense
@constraint(sp, u_dispense <= x_inventory.in) # dispensing more than inventory does not make sense
@constraint(sp, x_inventory.out == x_inventory.in - u_dispense + x_transit[end].in)
# stock outs for patients
@constraint(sp, h_lost == ω_demand - u_dispense)
# costs that are generic for all time steps
@expression(sp, ordinary_costs, (h_lost * stockout_cost) + (u_order * order_cost) + (u_order_size * order_size_cost))
if t == time_horizon
# on final time point, we have to throw out all the unused inventory
@stageobjective(sp, ordinary_costs + (x_inventory.out * wastage_cost))
else
# on intermediate times, we incur storage costs
@stageobjective(sp, ordinary_costs + (x_inventory.out * storage_cost))
end
end
SDDP.numerical_stability_report(model)
SDDP.train(
model;
iteration_limit = 100,
risk_measure = SDDP.Expectation()
# risk_measure = SDDP.AVaR(0.5)
)
# out of sample MC simulation
sampler = SDDP.OutOfSampleMonteCarlo(
model,
use_insample_transition=true
) do t
Ω = range(quantile.(Poisson(demand_λ[t]), [0.01, 0.99])...) # get most of the probability mass
P = pdf.(Poisson(demand_λ[t]), Ω)
P = P ./ sum(P)
return [SDDP.Noise(ω, p) for (ω, p) in zip(Ω, P)]
end
n_replications = 100
simulations = SDDP.simulate(model, n_replications, [:x_inventory, :x_transit, :u_dispense, :u_order_size, :h_lost], sampling_scheme = sampler);
# ribbon plot
p1 = Plots.plot(
SDDP.publication_plot(simulations; title = "inventory") do sim
return sim[:x_inventory].in
end;
xlabel = "Time"
)
p2 = Plots.plot(
SDDP.publication_plot(simulations; title = "dispense") do sim
return sim[:u_dispense]
end;
xlabel = "Time"
)
p3 = Plots.plot(
SDDP.publication_plot(simulations; title = "order size") do sim
return sim[:u_order_size]
end;
xlabel = "Time"
)
p4 = Plots.plot(
SDDP.publication_plot(simulations; title = "stockouts") do sim
return sim[:h_lost]
end;
xlabel = "Time"
)
plot(p1,p2,p3,p4, layout=(2,2), size=(1024,768))