Stochastic Lot-Sizing with Fixed Ordering Costs

This example shows how to train target-state decision rules for a stochastic inventory problem with ex-ante ordering decisions.

The example has two purposes:

show the complete optimization model before discussing implementation details; and
show the code in the same order a reader would run it.

using DecisionRules
using Flux
using HiGHS
using JuMP
using Random
using Statistics

The runnable experiment lives outside the documentation tree. The file defines the demand process, JuMP builders, and policy architecture used below.

include(joinpath(@__DIR__, "..", "..", "..", "examples", "inventory_control",
    "build_inventory_problem.jl"))

Main.build_lstm_exante_policy

Information Pattern

At the beginning of period t, the controller knows

\[x_t = (s_{t-1}, d_{t-1}, d_{t-2}),\]

where s is net inventory and d is realized demand. The controller chooses the order quantity before seeing current demand d_t. This is an ex-ante decision.

The neural policy receives [d_t, x_t...] during training because DecisionRules policies output target states after the stage uncertainty is sampled. The implementation below uses that target only to guide the optimization model; the actual order still respects the model's information pattern.

Complete Stage Model

For each period t = 1, ..., T, the stage model is

\[\begin{aligned} \min_{q_t,z_t,s_t^{mid},s_t,h_t,b_t} \quad & K z_t + c q_t + h h_t + p b_t + \lambda |s_t^{mid} - \hat{s}_t| \\ \text{s.t.}\quad & 0 \le q_t \le Q_{\max} z_t, && \text{(1) order capacity} \\ & z_t \in \{0,1\}, && \text{(2) setup decision} \\ & s_t^{mid} = s_{t-1} + q_t, && \text{(3) order arrives} \\ & s_t = s_t^{mid} - d_t, && \text{(4) demand realizes} \\ & h_t - b_t = s_t, && \text{(5) inventory split} \\ & h_t \ge 0,\; b_t \ge 0. && \text{(6) split bounds} \end{aligned}\]

The relaxed model removes (2) and replaces (1) by $0 \le q_t \le Q_{\max}$; it also removes the fixed cost K z_t from the objective.

The target \hat{s}_t is not an operational requirement. It is the state target produced by the neural decision rule, and the penalty term gives the policy a gradient signal.

Parameters

inventory_parameters = (
    T = INVENTORY_T,
    setup_cost = INVENTORY_K,
    unit_order_cost = INVENTORY_C,
    holding_cost = INVENTORY_H,
    backlog_cost = INVENTORY_P,
    order_capacity = INVENTORY_Q_MAX,
    initial_inventory = INVENTORY_I0,
    target_penalty = INVENTORY_PENALTY,
)

(T = 12, setup_cost = 500.0, unit_order_cost = 2.0, holding_cost = 1.0, backlog_cost = 25.0, order_capacity = 350.0, initial_inventory = 30.0, target_penalty = 75.0)

Demand Process

Demand has a hidden seasonal phase, a persistent hidden regime, and an AR(1) shock:

\[\epsilon_t = 0.92 \epsilon_{t-1} + 0.35 \eta_t,\]

\[d_t = \operatorname{clip}\!\left( m_{\kappa_t} + w_{\kappa_t}(0.85 r_t + 0.42 \epsilon_t + 0.12 \eta'_t) \right),\]

where r_t is the hidden regime and $\kappa_t = 1 + ((t + \phi - 1) \bmod T)$ is the hidden seasonal index.

Random.seed!(11)
demand_paths = [sample_inventory_demand_path() for _ in 1:3]

3-element Vector{Vector{Float64}}:
 [140.27957223310318, 120.13602692784835, 84.52413828360736, 66.42497332612164, 30.203507099786783, 21.532366274129416, 23.569729350226368, 25.382477943381183, 36.62016501123337, 59.941965455171044, 84.87346850226275, 116.74446297187289]
 [93.00459235715715, 79.48901751183757, 66.51138359599099, 47.56538250301797, 30.26403947318441, 22.000849292036136, 18.0501383475498, 21.10920899665569, 27.249268559662696, 38.384409214789926, 59.21646965840418, 77.88399163890075]
 [34.98148478302812, 28.301769921307635, 31.919096823374215, 43.22633221808005, 66.37498702882962, 92.31135599583436, 114.47850802136863, 134.47476479238463, 109.50631135153354, 87.58163704187803, 64.33436178525484, 60.45803932952831]

Build the Continuous and Integer Models

The builders return the JuMP model(s), input-state parameters, output-target parameters, an uncertainty sampler, and the initial state.

relaxed_subproblems,
relaxed_state_in,
relaxed_state_out,
relaxed_sampler,
initial_state = build_inventory_subproblems(;
    num_scenarios = 100,
    integer = false,
)

integer_subproblems,
integer_state_in,
integer_state_out,
integer_sampler,
_ = build_inventory_subproblems(;
    num_scenarios = 100,
    integer = true,
)

(JuMP.Model[A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z], Vector{Any}[[s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in]], Vector{Tuple{Any, JuMP.VariableRef}}[[(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)]], Main.InventoryProcessSampler(JuMP.VariableRef[demand, demand, demand, demand, demand, demand, demand, demand, demand, demand, demand, demand]), [30.0, 0.0, 0.0])

The deterministic equivalent is the full-horizon model used by direct transcription training.

integer_det_equivalent,
integer_det_state_in,
integer_det_state_out,
integer_det_sampler,
_ = build_inventory_det_equivalent(;
    num_scenarios = 50,
    integer = true,
)

(A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 171
├ num_constraints: 219
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 72
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 12
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 12
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 48
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 12
│ ├ JuMP.VariableRef in MOI.ZeroOne: 12
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 51
└ Names registered in the model
  └ :back, :deficit_arr, :demand, :inv_hold, :last_demand_init, :last_demand_out, :last_demand_target, :norm_deficit_arr, :prev_demand_init, :prev_demand_out, :prev_demand_target, :q, :s_init, :s_mid, :s_out, :s_target, :z, Vector{Any}[[s_init, last_demand_init, prev_demand_init], [s_out[1], last_demand_out[1], prev_demand_out[1]], [s_out[2], last_demand_out[2], prev_demand_out[2]], [s_out[3], last_demand_out[3], prev_demand_out[3]], [s_out[4], last_demand_out[4], prev_demand_out[4]], [s_out[5], last_demand_out[5], prev_demand_out[5]], [s_out[6], last_demand_out[6], prev_demand_out[6]], [s_out[7], last_demand_out[7], prev_demand_out[7]], [s_out[8], last_demand_out[8], prev_demand_out[8]], [s_out[9], last_demand_out[9], prev_demand_out[9]], [s_out[10], last_demand_out[10], prev_demand_out[10]], [s_out[11], last_demand_out[11], prev_demand_out[11]]], Vector{Tuple{Any, JuMP.VariableRef}}[[(s_target[1], s_out[1]), (last_demand_target[1], last_demand_out[1]), (prev_demand_target[1], prev_demand_out[1])], [(s_target[2], s_out[2]), (last_demand_target[2], last_demand_out[2]), (prev_demand_target[2], prev_demand_out[2])], [(s_target[3], s_out[3]), (last_demand_target[3], last_demand_out[3]), (prev_demand_target[3], prev_demand_out[3])], [(s_target[4], s_out[4]), (last_demand_target[4], last_demand_out[4]), (prev_demand_target[4], prev_demand_out[4])], [(s_target[5], s_out[5]), (last_demand_target[5], last_demand_out[5]), (prev_demand_target[5], prev_demand_out[5])], [(s_target[6], s_out[6]), (last_demand_target[6], last_demand_out[6]), (prev_demand_target[6], prev_demand_out[6])], [(s_target[7], s_out[7]), (last_demand_target[7], last_demand_out[7]), (prev_demand_target[7], prev_demand_out[7])], [(s_target[8], s_out[8]), (last_demand_target[8], last_demand_out[8]), (prev_demand_target[8], prev_demand_out[8])], [(s_target[9], s_out[9]), (last_demand_target[9], last_demand_out[9]), (prev_demand_target[9], prev_demand_out[9])], [(s_target[10], s_out[10]), (last_demand_target[10], last_demand_out[10]), (prev_demand_target[10], prev_demand_out[10])], [(s_target[11], s_out[11]), (last_demand_target[11], last_demand_out[11]), (prev_demand_target[11], prev_demand_out[11])], [(s_target[12], s_out[12]), (last_demand_target[12], last_demand_out[12]), (prev_demand_target[12], prev_demand_out[12])]], Main.InventoryProcessSampler(JuMP.VariableRef[demand[1], demand[2], demand[3], demand[4], demand[5], demand[6], demand[7], demand[8], demand[9], demand[10], demand[11], demand[12]]), [30.0, 0.0, 0.0])

Integer Sensitivity Strategies

Mixed-integer models do not have ordinary LP duals. DecisionRules therefore makes the chosen postprocessing strategy explicit.

fixed_discrete = FixedDiscreteIntegerStrategy()
continuous_relaxation = ContinuousRelaxationIntegerStrategy()

ContinuousRelaxationIntegerStrategy()

FixedDiscreteIntegerStrategy solves the MIP, fixes the incumbent integer variables, re-solves the fixed LP, and reads local dual information.

ContinuousRelaxationIntegerStrategy relaxes integer variables first and reads duals from the relaxed LP. This is smoother and faster, but the gradient is for the relaxation, not for an integer-feasible decision.

Score-Function Correction

Local LP duals do not see a discrete switch such as "open the setup variable". A score-function correction estimates the effect of target changes by solving perturbed integer rollouts:

\[\nabla L = \alpha \nabla L_{\mathrm{dual}} + (1-\alpha) \frac{1}{M} \sum_{m=1}^{M} (R_m - b) \nabla_\theta \sum_{t=1}^{T} \left\langle \delta_{m,t}/\sigma^2, \hat{x}_{t+1}(\theta) \right\rangle .\]

There are two different solves in the mixed-gradient training loop:

train_multistage(...; integer_strategy = fixed_discrete) controls the deterministic-equivalent solve used for the dual-gradient term $\nabla L_{\mathrm{dual}}$. This solve needs a postprocessing strategy because duals are not directly defined for a MIP.
ScoreFunctionConfig(integer_subproblems, ...) controls the Monte Carlo rollout term. These rollout models are solved exactly as they are built. Because integer_subproblems contain binary setup variables, the rollout costs R_m are true MIP rollout costs.

In short: integer_strategy is for reading local duals; score-function rollouts are for measuring realized costs.

score_function = ScoreFunctionConfig(
    integer_subproblems,
    integer_state_in,
    integer_state_out;
    dual_weight = 0.5,
    perturbation_std = 1.0,
    num_rollouts = 8,
)

score_schedule = ScoreFunctionSchedule(
    score_function;
    sf_start = 200,
    ramp_batches = 300,
    perturbation_std_initial = 0.1,
    num_rollouts_initial = 2,
)

ScoreFunctionSchedule(ScoreFunctionConfig(JuMP.Model[A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z, A JuMP Model
├ solver: HiGHS
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 17
├ num_constraints: 21
│ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 6
│ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 1
│ ├ Vector{JuMP.VariableRef} in MOI.NormOneCone: 1
│ ├ JuMP.VariableRef in MOI.GreaterThan{Float64}: 4
│ ├ JuMP.VariableRef in MOI.LessThan{Float64}: 1
│ ├ JuMP.VariableRef in MOI.ZeroOne: 1
│ └ JuMP.VariableRef in MOI.Parameter{Float64}: 7
└ Names registered in the model
  └ :_deficit, :back, :demand, :inv_hold, :last_demand_in, :last_demand_out, :last_demand_target, :norm_deficit, :prev_demand_in, :prev_demand_out, :prev_demand_target, :q, :s_in, :s_mid, :s_out, :s_target, :z], Vector{Any}[[s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in], [s_in, last_demand_in, prev_demand_in]], Vector{Tuple{Any, JuMP.VariableRef}}[[(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)], [(s_target, s_out), (last_demand_target, last_demand_out), (prev_demand_target, prev_demand_out)]], 0.5, 1.0, 8, :mean), 200, 300, 0.5, 0.1, 1.0, 2, 8)

Policy

A DecisionRules policy is any callable π(x) -> target where x is the concatenation [uncertainty..., state...] and target is the desired next state. The only requirement is that it is differentiable via Zygote.gradient and registered with Functors.@functor so that Flux.loadmodel! can checkpoint its parameters.

Feedforward policy

The simplest architecture is a feedforward MLP. This policy is ex-ante: it ignores the current demand d_t (index 1) and uses only the state entries [inventory, d_{t-1}, d_{t-2}].

using Functors: @functor

struct ExAntePolicy{N}
    net::N
end

@functor ExAntePolicy (net,)

The callable normalizes features to ≈[0,1] and maps through the network. The sigmoid output bounds the target to [0, 500].

function (p::ExAntePolicy)(x)
    inventory = Float32(x[2])
    d_prev    = Float32(x[3])
    d_prev2   = Float32(x[4])
    features  = Float32[inventory / 100, d_prev / 100, d_prev2 / 100]
    target    = 500f0 .* Flux.sigmoid.(p.net(features))
    return Float32[target[1], x[1], d_prev]
end

Random.seed!(2024)
policy = ExAntePolicy(Chain(Dense(3, 32, relu), Dense(32, 24, relu), Dense(24, 1)))

Main.ExAntePolicy{Flux.Chain{Tuple{Flux.Dense{typeof(NNlib.relu), Matrix{Float32}, Vector{Float32}}, Flux.Dense{typeof(NNlib.relu), Matrix{Float32}, Vector{Float32}}, Flux.Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}}}}(Chain(Dense(3 => 32, relu), Dense(32 => 24, relu), Dense(24 => 1)))

Recurrent (LSTM) policy

When the uncertainty process has temporal structure (regimes, trends, seasonality), a recurrent encoder can learn patterns that a feedforward MLP cannot detect from a fixed-length window.

The design below uses Flux.LSTMCell to process one lagged demand value per stage. The LSTM hidden state accumulates across stages within a scenario, then resets between scenarios via Flux.reset!.

The affine output raw × 200 + 150 avoids sigmoid saturation and centers the target on typical inventory levels.

mutable struct RecurrentExAntePolicy{E,C,S}
    encoder::E
    combiner::C
    state::S
end

@functor RecurrentExAntePolicy (encoder, combiner)

function (p::RecurrentExAntePolicy)(x)
    d_prev    = Float32(x[3])
    inventory = Float32(x[2])
    d_prev2   = Float32(x[4])
    T = eltype(first(p.state))
    encoded, new_state = p.encoder(T[d_prev / 100], p.state)
    p.state = new_state
    raw = p.combiner(vcat(encoded, T[inventory / 100, d_prev2 / 100]))
    target = raw[1] * 200f0 + 150f0
    return Float32[target, x[1], d_prev]
end

function Flux.reset!(p::RecurrentExAntePolicy)
    p.state = Flux.initialstates(p.encoder)
    return nothing
end

Random.seed!(2024)
lstm_encoder = Flux.LSTMCell(1 => 16)
lstm_policy = RecurrentExAntePolicy(
    lstm_encoder,
    Dense(16 + 2, 1),
    Flux.initialstates(lstm_encoder),
)

Main.RecurrentExAntePolicy{Flux.LSTMCell{Matrix{Float32}, Matrix{Float32}, Vector{Float32}}, Flux.Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}, Tuple{Vector{Float32}, Vector{Float32}}}(LSTMCell(1 => 16), Dense(18 => 1), (Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]))

Training Calls

The continuous problem uses ordinary dual information.

train_multistage(
    policy,
    initial_state,
    relaxed_subproblems,
    relaxed_state_in,
    relaxed_state_out,
    relaxed_sampler;
    num_batches = 400,
    num_train_per_batch = 5,
    optimizer = Flux.Adam(0.0015),
    integer_strategy = NoIntegerStrategy(),
    penalty_schedule = [(1, 80, 0.4), (81, 400, 1.0)],
)

The integer deterministic-equivalent run uses the fixed-discrete local dual path plus the scheduled score-function correction.

train_multistage(
    policy,
    initial_state,
    integer_det_equivalent,
    integer_det_state_in,
    integer_det_state_out,
    integer_det_sampler;
    num_batches = 800,
    num_train_per_batch = 10,
    optimizer = Flux.Adam(0.0008),
    integer_strategy = fixed_discrete,
    penalty_schedule = [(1, 120, 0.4), (121, 800, 1.0)],
    score_function = score_schedule,
)

Evaluation

A trained policy should be evaluated by stage-wise rollout, because that is the deployment semantics: solve one period, observe the realized next state, then solve the next period.

uncertainty_sample = sample(integer_sampler)
rollout_cost = simulate_multistage(
    integer_subproblems,
    integer_state_in,
    integer_state_out,
    initial_state,
    uncertainty_sample,
    policy;
    integer_strategy = fixed_discrete,
)

10134.684587793752

Experiment Scripts

Each variant can be trained independently via SLURM or directly:

# Single variant
julia --project=. train_dr_inventory.jl integer_lstm

# All variants in parallel via SLURM
cd examples/inventory_control && bash launch_all.sh

Available variant tags: relaxed, relaxed_lstm, relaxed_hp, relaxed_lstm_hp, integer, integer_cr, integer_sf, integer_hp, integer_lstm, integer_lstm_sf.

After training, run the comparison script to regenerate tables and figures:

julia --project=. evaluate_inventory.jl
julia --project=. solve_sddp.jl
julia --project=. compare_results.jl

The figures used by this page are generated by compare_results.jl.

Demand process

Relaxed results

Integer results

Relaxed (continuous) results

SDDP uses a PAR(1) approximation of the true latent demand process, which is not exact for this problem. Despite this advantage for TS-DDR, the gap between the best TS-DDR variant and SDDP is ~7%.

The LSTM encoder closes ~25% of the gap versus the feedforward baseline by learning temporal demand patterns from lagged observations.

Method	N	Mean cost	Std	vs SDDP
SDDP (PAR)	300	2434.0	—	0.0%
TS-DDR (LSTM)	300	2610.6	540.3	+7.3%
TS-DDR (feedforward)	300	2667.3	593.5	+9.6%
TS-DDR (HighPenalty)	300	2677.5	547.0	+10.0%
TS-DDR (LSTM+HP)	300	2712.0	554.6	+11.4%

Integer (MIP) results

SDDP uses an AlternativeForwardPass: MIP in the forward pass, LP relaxation in the backward pass for valid cuts. The TS-DDR gap is ~36%.

Method	N	Mean cost	Std	vs SDDP
SDDP (MIP fwd)	300	5871.6	1087.4	0.0%
TS-DDR (FixedDiscrete)	300	8015.8	718.3	+36.5%
TS-DDR (MixedGrad)	300	8268.0	715.3	+40.8%
TS-DDR (ContRelax)	300	8318.1	718.8	+41.7%
TS-DDR (HighPenalty)	300	8388.4	615.9	+42.8%
SDDP (LP relax)	300	8274.2	912.5	+40.9%
Base-stock (S*=160)	300	9035.6	506.8	+53.9%
Random (untrained)	300	9594.6	361.1	+63.4%