diff --git a/docs/src/api.md b/docs/src/api.md index 873c3ab3..260cb49e 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -86,6 +86,13 @@ Modules = [DecisionFocusedLearningBenchmarks.StochasticVehicleScheduling] Private = false ``` +### Dynamic Replenishment + +```@autodocs +Modules = [DecisionFocusedLearningBenchmarks.DynamicReplenishment] +Private = false +``` + ### Warcraft ```@autodocs @@ -103,6 +110,7 @@ Modules = [ DecisionFocusedLearningBenchmarks.ContextualStochasticArgmax, DecisionFocusedLearningBenchmarks.DynamicVehicleScheduling, DecisionFocusedLearningBenchmarks.DynamicAssortment, + DecisionFocusedLearningBenchmarks.DynamicReplenishment, DecisionFocusedLearningBenchmarks.FixedSizeShortestPath, DecisionFocusedLearningBenchmarks.Maintenance, DecisionFocusedLearningBenchmarks.PortfolioOptimization, diff --git a/docs/src/benchmarks/dynamic/04_replenishment.jl b/docs/src/benchmarks/dynamic/04_replenishment.jl new file mode 100644 index 00000000..dbdc95e6 --- /dev/null +++ b/docs/src/benchmarks/dynamic/04_replenishment.jl @@ -0,0 +1,230 @@ +# # Dynamic Replenishment +# A retailer must decide how many units of each item to reorder at each time step. Demand +# follows an endogenous customer choice model, in which the purchase probability depends on +# the items currently available and on the stock levels. +# The objective is to maximize total revenue over a finite horizon, defined as the sales +# margin minus the stock costs. +# Replenished items take a certain amount of time to reach the store, but they can already be +# sold while in transit. The physical inventory has soft lower and upper bounds, and violating +# them incurs a large penalty. +# Items are also subject to coupling production constraints with random quotas. + +using DecisionFocusedLearningBenchmarks +using Plots + +b = DynamicReplenishmentBenchmark() + +# ## Observable input +# +# Generate one environment and roll it out with the random policy to collect a sample +# trajectory. At each step the agent observes item prices and features, the current virtual +# and physical stock levels, the sales, replenishment and stock histories, and the remaining +# quotas: +policies = generate_baseline_policies(b) +env = generate_environments(b, 1)[1] +_, trajectory = evaluate_policy!(policies.random, env) + +# The observable state at step 1: stock levels (virtual stock includes the units still in +# transit, physical stock only the units already in the store) together with the static +# utility of each item, which drives the customer choice model: +plot_context(b, trajectory[1]) + +# ## A training sample +# +# Each step in a trajectory is a labeled tuple `(x, (θ, η), y)` plus state and reward: +# - `x`: `(d+19) × ∑ᵢ ubᵢ` feature matrix per step, with one column per candidate stock level +# of each item (`ubᵢ` is the replenishment upper bound of item ``i``). The rows hold the +# static features (price and item features), the dynamic item features (current stock, mean +# sales, mean stock, mean number of customers, days on lot) and the stock-level features +# (deviations from the stock bounds, the mean stock and the quota). The last row is the item +# identifier. +# - `(θ, η)`: predicted utility scores. `θ` is the predicted utility of each item, `η` the +# predicted marginal cost of each additional unit held in stock. +# - `y`: replenishment decision at this step (vector of length ``N``) +# - `instance`: the state, containing the physical and virtual stock levels together with the +# sales, replenishment and stock histories +# - `reward`: sales margin minus stock costs at time step ``t`` +# +# One step in the trajectory: +plot_sample(b, trajectory[1]) + +# A few steps side by side: +plot_trajectory(b, trajectory[1:min(4, length(trajectory))]) + +# ## DFL pipeline components + +# The DFL agent chains two components: a neural network predicting utility scores per item: +model = generate_statistical_model(b) # state features → predicted utility scores (θ, η) +# and a maximizer choosing the best replenishment decision based on the predicted scores, +# the stock levels and the production constraints: +maximizer = generate_maximizer(b) + +# At each step, the model maps the current state (prices, features, stock levels, quotas) to +# the utility scores ``(\theta, \eta)``. The maximizer then selects the feasible replenishment +# decision that maximizes the total predicted utility. + +# --- +# ## Problem Description +# +# ### Overview +# +# In the **Dynamic Replenishment problem**, a retailer has a catalog of ``N`` items, of which +# only a subset is present in its inventory at any given time. The inventory splits into units +# that are physically in the store and units that are still in transit; the items offered to +# customers are all those with a positive virtual or physical inventory. +# The retailer pays a holding cost for every unit in stock and for every unit in transit, and +# it faces soft lower and upper bounds on the physical stock level: violating either bound +# incurs a penalty. +# At each time step the retailer decides how many units of each item to replenish, subject to +# coupling production quotas that limit how many units can be reordered. +# Customer demand is stochastic and follows a multinomial logit choice model. +# +# The problem is characterized by: +# - **Endogenous noise**: what customers buy depends on which items are actually available, hence on the past replenishment decisions +# - **Combinatorial action space**: the number of feasible replenishment decisions is exponential in the number of items +# +# ### Mathematical Formulation +# +# **State** ``s_t = (p, f, vs_t, ps_t, t, \mathcal{H}_t^s, \mathcal{H}_t^r, \mathcal{H}_t^p, \mathcal{H}_t^c)`` where: +# - ``p``: fixed item prices +# - ``f``: static item features +# - ``vs_t``: current virtual stock levels +# - ``ps_t``: current physical stock levels +# - ``\mathcal{H}_t^s``: stock history +# - ``\mathcal{H}_t^r``: revenue history +# - ``\mathcal{H}_t^p``: purchase history +# - ``\mathcal{H}_t^c``: customer history +# - ``t``: current time step +# +# **Action:** ``a_t \in \mathbb{N}^N`` with ``A a_t \leq b_t``, where ``A`` is the (coupling) production constraint matrix and ``b_t`` the quotas of each constraint at time step ``t``. +# +# **Customer choice** (multinomial logit): each item is assigned a static utility score ``v_i`` +# given by the customer choice model. The default model is linear, so the utility of an item is +# a linear combination of its features. Only the items that are actually in the inventory can +# be bought, that is the offer set ``\mathcal{O}_t = \{i : vs_t^i > 0\}``. The probability that +# a customer purchases item ``i`` at time step ``t`` is: +# ```math +# \mathbb{P}(i \mid s_t) = \frac{\exp(v_i)}{\sum_{j \in \mathcal{O}_t} \exp(v_j) + 1} +# ``` +# where the ``+1`` in the denominator accounts for the no-purchase option. +# The Gumbel-max trick is used to sample the purchased item from this distribution: customer +# ``k`` purchases the item ``i^\star`` such that +# ```math +# i^\star = \operatorname*{argmax}_{i \in \mathcal{O}_t \cup \{0\}} \left(v_i + \epsilon_i^k\right) +# ``` +# where ``\epsilon_i^k`` is a Gumbel random variable and ``0`` denotes the no-purchase option. +# At each time step, a random number of customers (drawn from a Poisson distribution of rate +# ``\lambda``) arrives and makes purchases. We write ``q_t^i`` for the number of units of item +# ``i`` purchased at time step ``t``. +# +# **Transition dynamics:** for each item ``i``, with ``\tau`` the delivery delay: +# - ``vs_{t+1}^i = vs_t^i + a_t^i - q_t^i`` +# - ``ps_{t+1}^i = \max(0, ps_t^i + a_{t-\tau}^i - q_t^i)`` +# +# **Reward:** for each item ``i``: +# - the sales margin is ``q_t^i m_i`` +# - the virtual stock cost is ``c_{vs}^i \, vs_t^i`` +# - the physical stock cost is ``c_{ps}^i \, ps_t^i`` +# - the penalty for violating the soft stock bounds is ``c_{lb}^i \max(0, lb_i - ps_t^i) + c_{ub}^i \max(0, ps_t^i - ub_i)`` +# +# The total reward at time step ``t`` is therefore: +# ```math +# r(s_t, a_t) = \sum_{i=1}^N q_t^i m_i - c_{vs}^i \, vs_t^i - c_{ps}^i \, ps_t^i - c_{lb}^i \max(0, lb_i - ps_t^i) - c_{ub}^i \max(0, ps_t^i - ub_i) +# ``` +# +# **Objective:** +# ```math +# \max_\pi \; \mathbb{E}\!\left[\sum_{t=1}^T r(s_t, \pi(s_t))\right] +# ``` +# +# ## Key Components +# +# ### [`DynamicReplenishmentBenchmark`](@ref) +# +# | Parameter | Description | Default | +# |-----------|-------------|---------| +# | `N` | Number of items in the catalog | 10 | +# | `λ` | Poisson arrival rate of customers per step | 15 | +# | `d` | Static feature dimension per item (in addition to price) | 5 | +# | `nb_constraints` | Number of coupling production constraints | 2 | +# | `constraints_matrix` | Coupling matrix ``A`` (`nb_constraints × N`) | random 0/1 matrix | +# | `quotas` | Quotas ``b_t`` per constraint and per step (`max_steps × nb_constraints`) | random in ``[10, 30]`` | +# | `stock_inf` | Soft lower bound on the physical stock | 0 | +# | `stock_sup` | Soft upper bound on the physical stock | 30 | +# | `ub_same_item` | Upper bound on the number of units of the same item | 30 | +# | `delivery_delay` | Delivery delay ``\tau``, in time steps | 3 | +# | `max_steps` | Steps per episode | 10 | +# | `customer_choice_model` | Model mapping item features to static utilities | random linear model | +# +# Prices are drawn uniformly in ``[1, 10]`` and item features uniformly in ``[-10, 10]``. The +# stock costs derive from the prices (``c_{vs}^i = p_i / 10T`` and ``c_{ps}^i = p_i / 5T``), +# and the bound violation cost is ``\max_i p_i``. The static utilities are obtained by applying +# the customer choice model to the standardized features, with a ``0`` appended for the +# no-purchase option. +# +# Either provide only `nb_constraints`, in which case a random constraints matrix and random +# quotas are generated, or provide both `constraints_matrix` and `quotas` and they are used +# as is. +# +# ### State Observation +# +# Agents observe a ``(d+19) \times \sum_i ub_i`` feature matrix, with one column per candidate +# stock level of each item. Each column concatenates: +# - the ``d+1`` static features of the item (price and item features) +# - 9 dynamic item features: current stock, mean sales, mean stock, mean number of past customers, and mean days on lot (each also scaled by the price) +# - 8 stock-level features: the deviations from `stock_inf`, from `stock_sup`, from the mean stock and from the quota of the step (each also scaled by the price) +# - the item identifier, used by the statistical model to group the columns per item +# +# ### Environment Generation +# +# Each environment starts from an initial stock drawn uniformly in ``\{0,\ldots,5\}`` per item, +# from which the virtual and physical stocks, the histories and the per-item replenishment +# upper bounds are initialized. A scenario is sampled at the same time and fixes, for the whole +# episode, the number of customers arriving at each step (Poisson with rate ``\lambda``) and +# the Gumbel perturbation of the utilities of each customer. Resetting the environment restores +# the initial stock and, by default, resamples the scenario. + +# ## Baseline Policies +# +# | Policy | Description | +# |--------|-------------| +# | Greedy | Solves the replenishment problem with the prices as item utilities and no stock penalization, which favors the most expensive items | +# | Random | Goes through the items in a random order and replenishes a random feasible quantity of each | +# | Lazy | Never replenishes anything | +# | SAA | Solves a multi-stage sample average approximation of the problem over the remaining horizon, on a set of sampled scenarios, and applies the first-stage decision | +# +# ## DFL Policy +# +# ```math +# \xrightarrow[\text{State}]{s_t} +# \fbox{Neural network $\varphi_w$} +# \xrightarrow[\text{Utilities}]{(\theta, \eta) \in \mathbb{R}^N \times \mathbb{R}^{\sum_i ub_i}} +# \fbox{Replenishment maximizer} +# \xrightarrow[\text{Replenishment}]{a_t} +# ``` +# +# **Model:** two heads sharing the same feature matrix: +# - `θ_model = Chain(Dense(d+10 => 1))`: one replenishment utility ``\theta_i`` per item, from the item features only +# - `η_model = Chain(Dense(d+18 => 1), softplus)`: one nonnegative marginal stock cost ``\eta_{i,j}`` per candidate stock level ``j`` of item ``i`` +# +# **Maximizer:** the predicted scores parametrize the linear objective of a MILP solved at each +# step, where ``y_i`` is the replenishment of item ``i``, ``s_i`` its current stock and +# ``z_{i,j} = 1`` if the post-replenishment stock of item ``i`` reaches at least ``j`` units: +# ```math +# \begin{aligned} +# \max_{y, z} \quad & \sum_{i = 1}^N \overbrace{\theta_i y_i}^{\text{replenishment revenue}} +# + \overbrace{\underbrace{\eta_{i,1} z_{i,1}}_{\text{intercept}} +# - \sum_{j = 2}^{ub_i} \underbrace{z_{i,j} \sum_{k = 2}^{j} \eta_{i,k}}_{\text{decreasing slope}}}^{\text{stock cost}} \\ +# \text{s.t.} \quad +# & y_i + s_i = \sum_{j = 1}^{ub_i} z_{i,j}, \quad \forall i \in [N] \\ +# & z_{i,j} \geq z_{i,j+1}, \quad \forall i \in [N], \ j \in [ub_i - 1] \\ +# & A y \leq b_t \\ +# & y_i \in \{0, \ldots, ub_i\}, \quad z_{i,j} \in \{0, 1\} +# \end{aligned} +# ``` +# The first constraint links the ``z`` variables to the post-replenishment stock, the second +# enforces that they are non-increasing in ``j`` (so they encode a stock level rather than an +# arbitrary set), and the third is the coupling quota constraint of the current step. Since +# ``\eta \geq 0``, the stock cost is a concave piecewise-linear function of the stock level, +# which lets the model penalize large inventories without making the problem nonlinear. +# diff --git a/ext/DFLBenchmarksPlotsExt.jl b/ext/DFLBenchmarksPlotsExt.jl index 117d1747..3c8b4161 100644 --- a/ext/DFLBenchmarksPlotsExt.jl +++ b/ext/DFLBenchmarksPlotsExt.jl @@ -23,5 +23,6 @@ include("plots/svs_plots.jl") include("plots/dvs_plots.jl") include("plots/dynamic_assortment_plots.jl") include("plots/maintenance_plots.jl") +include("plots/dynamic_replenishment_plots.jl") end diff --git a/ext/plots/dynamic_replenishment_plots.jl b/ext/plots/dynamic_replenishment_plots.jl new file mode 100644 index 00000000..5286b324 --- /dev/null +++ b/ext/plots/dynamic_replenishment_plots.jl @@ -0,0 +1,179 @@ +has_visualization(::DynamicReplenishmentBenchmark) = true + +function plot_context(bench::DynamicReplenishmentBenchmark, sample::DataSample; kwargs...) + static_u = bench.static_utilities[1:(end - 1)] # drop the "no purchase" option + state = sample.state + stock = Float64.(state.stock) + stock_p = Float64.(state.physical_stock) + N = length(stock) + + p1 = bar( + 1:N, + stock; + label="Virtual stock", + color="#5fd6a8", + ylabel="Count", + title="Stock levels", + xticks=(1:N, fill("", N)), + ) + bar!(p1, 1:N, stock_p; label="Physical stock", color="#0f7d52") + + p2 = bar( + 1:N, + static_u; + legend=false, + xlabel="Item", + ylabel="Utility", + title="Static utilities", + color="#2a78d6", + ) + + l = Plots.@layout [a{0.6h}; b{0.4h}] + return Plots.plot(p1, p2; layout=l, size=(800, 600), kwargs...) +end + +function bar_plot_stock_repl_sales( + stock, + stock_p, + repl, + sales=nothing, + nb_customers=nothing; + xlabel::String="Item", + ylabel::String="Count", + title::String="Stock, replenishment and sales", + legend, + kwargs..., +) + xmax = length(stock) + if sales !== nothing || nb_customers !== nothing + w = 0.5 + xs_left = (1:xmax) .- w/2 + xs_right = (1:xmax) .+ w/2 + else + w = 1 + xs_left = 1:xmax + end + + p = bar( + xs_left, + stock .+ repl; + bar_width=w, + label="Replenishment", + color="#2a78d6", # blue + xlabel=xlabel, + ylabel=ylabel, + title=title, + legend=legend, + xticks=1:xmax, + size=(800, 500), + ) + bar!(p, xs_left, stock; bar_width=w, label="Virtual stock", color="#5fd6a8") + bar!(p, xs_left, stock_p; bar_width=w, label="Physical stock", color="#0f7d52") + + if nb_customers !== nothing + bar!(p, xs_right, -nb_customers; bar_width=w, label="No buy", color="#9a9a9a") + end + if sales !== nothing + bar!(p, xs_right, -sales; bar_width=w, label="Sales", color="#e34948") + end + return p +end + +""" +Bar plot of stock level of each items. +""" +function plot_sample( + b::DynamicReplenishmentBenchmark, + sample::DataSample; + with_legend=true, + with_title=true, + kwargs..., +) + state = sample.state + stock = Float64.(state.stock) + stock_p = Float64.(state.physical_stock) + repl = Float64.(sample.y) + sales = Float64.(sample.next_sales) + return p = bar_plot_stock_repl_sales( + stock, + stock_p, + repl, + sales, + nothing; + xlabel="Item", + ylabel="Count", + title=with_title ? "Stock, replenishment and sales" : "", + legend=with_legend ? :topright : false, + kwargs..., + ) +end + +""" +Plot a full episode. + +With `aggregated=true`, quantities are summed over items and shown time step by time +step; the shelf bounds `[stock_inf, stock_sup]` are drawn as two dashed lines, since they +constrain exactly that total and are what the over/under stock penalty is computed on. +Note that the physical stock of a step is the one observed *before* that step's sales, so +the bar to compare against the bounds for the penalty incurred at step `t` is the one at +`t + 1`. + +With `aggregated=false`, each time step gets its own item-by-item subplot; the bounds are +not drawn there, as they bear on the total rather than on any single item. +""" +function plot_trajectory( + bench::DynamicReplenishmentBenchmark, + trajectory::Vector{<:DataSample}; + max_steps=10, + cols=3, + aggregated::Bool=false, + kwargs..., +) + n = min(length(trajectory), max_steps) + rows = ceil(Int, n / cols) + steps = round.(Int, range(1, length(trajectory); length=n)) + upper_middle = div(cols, 2) + 1 + if aggregated + states = [sample.state for sample in trajectory[steps]] + stocks = [sum(state.stock) for state in states] + stocks_p = [sum(state.physical_stock) for state in states] + repls = [sum(sample.y) for sample in trajectory[steps]] + sales = [sum(sample.next_sales) for sample in trajectory[steps]] + nb_customers = [sample.customers for sample in trajectory[steps]] + p = bar_plot_stock_repl_sales( + stocks, + stocks_p, + repls, + sales, + nb_customers; + xlabel="Time step", + ylabel="Count", + title="Total stock, replenishment and sales over time", + legend=:topright, + kwargs..., + ) + # Single `hline!` call for both bounds: two lines, one legend entry. + hline!( + p, + Float64[bench.stock_inf, bench.stock_sup]; + linestyle=:dash, + color="#b23a48", + lw=2, + label="Shelf bounds", + ) + return p + else + plots = [ + plot_sample( + bench, + trajectory[t]; + with_legend=(t == 1), + with_title=(t == upper_middle), + kwargs..., + ) for t in steps + ] + return Plots.plot( + plots...; layout=(rows, cols), size=(cols * 300, rows * 250), kwargs... + ) + end +end diff --git a/src/DecisionFocusedLearningBenchmarks.jl b/src/DecisionFocusedLearningBenchmarks.jl index 7a68f6d2..ef7aba3c 100644 --- a/src/DecisionFocusedLearningBenchmarks.jl +++ b/src/DecisionFocusedLearningBenchmarks.jl @@ -59,6 +59,7 @@ include("ContextualStochasticArgmax/ContextualStochasticArgmax.jl") include("DynamicVehicleScheduling/DynamicVehicleScheduling.jl") include("DynamicAssortment/DynamicAssortment.jl") include("Maintenance/Maintenance.jl") +include("DynamicReplenishment/DynamicReplenishment.jl") using .Utils @@ -109,6 +110,7 @@ using .ContextualStochasticArgmax using .DynamicVehicleScheduling using .DynamicAssortment using .Maintenance +using .DynamicReplenishment export Argmax2DBenchmark export ArgmaxBenchmark @@ -122,5 +124,6 @@ export SubsetSelectionBenchmark export WarcraftBenchmark export MaintenanceBenchmark export ContextualStochasticArgmaxBenchmark +export DynamicReplenishmentBenchmark end # module DecisionFocusedLearningBenchmarks diff --git a/src/DynamicReplenishment/DynamicReplenishment.jl b/src/DynamicReplenishment/DynamicReplenishment.jl new file mode 100644 index 00000000..6ec84c3b --- /dev/null +++ b/src/DynamicReplenishment/DynamicReplenishment.jl @@ -0,0 +1,468 @@ +module DynamicReplenishment + +using ..Utils + +using JuMP: + Model, + @variable, + @objective, + @constraint, + optimize!, + value, + fix, + primal_status, + termination_status, + objective_value, + set_silent, + MOI, + AffExpr, + set_attribute, + set_start_value +using Random: Random, AbstractRNG, seed!, randperm, Xoshiro +using Distributions: Poisson, Uniform, Gumbel +using Flux: Chain, Dense, @layer, softplus, relu +using InferOpt: LinearMaximizer +using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES +using LinearAlgebra: dot, I +using Statistics: mean, quantile, std +using StatsBase: ZScoreTransform, fit, transform +""" +$TYPEDEF + +Benchmark for a replenishment (retail) problem with production constraints. +Items are chosen according to a a given customer choice model which is endogenous. + +# Fields +$TYPEDFIELDS +""" +struct DynamicReplenishmentBenchmark{M} <: AbstractDynamicBenchmark{true} + "customer choice model (price, features)" + customer_choice_model::M + "Poisson arrival rate of customers" + λ::Float64 + "number of items" + N::Int + "dimension of feature vectors (in addition to price: number of objects)" + d::Int + "Coupling matrix for quota constraints (nb_constraints x N)" + constraints_matrix::Matrix{Int} + "quotas for each constraint at each time step (max_steps x nb_constraints)" + quotas::Matrix{Int} + "Lower stock bound" + stock_inf::Int + "Upper stock bound" + stock_sup::Int + "upper bound of same item in stock" + ub_same_item::Int + "delivery delay in days" + delivery_delay::Int + "prices of the items (N)" + prices::Vector{Float64} + "items' features (d x N)" + features::Matrix{Float64} + "price and features, centered and scaled across items ((d+1) x N)" + scaled_features::Matrix{Float64} + "the transform mapping `[prices'; features]` to `scaled_features`" + feature_transform::ZScoreTransform{Float64,Vector{Float64}} + "cost of virtual stock (N)" + virtual_stock_cost::Vector{Float64} + "cost of physical stock (N)" + physical_stock_cost::Vector{Float64} + "over stock bound cost" + over_stock_bound_cost::Float64 + "number of steps per episode" + max_steps::Int + "max quota per time step per item (max_steps x N)" + max_quotas::Matrix{Int} + "static utilities of the items from the customer choice model (N)" + static_utilities::Vector{Float64} +end + +""" + DynamicReplenishmentBenchmark(; + N=10, + λ=15, + d=5, + nb_constraints=2, + constraints_matrix=nothing, + quotas=nothing, + stock_inf=0, + stock_sup=30, + ub_same_item=30, + delivery_delay=3, + max_steps=10, + over_stock_bound_cost=nothing + ) + +Constructor for [`DynamicReplenishmentBenchmark`](@ref). +By default, the benchmark has 10 items, feature dimension 5 (+1 for price), 10 steps per +episode, a simple linear customer choice model (all weights are negative), a poisson arrival of 15 customer per time step, and is endogenous. It generates +- random prices uniformly in [1, 10] +- random features uniformly in [-10, 10] +- stock costs are dependant on the price +The user can choose between +- only providing a number of constraints, in which case the constructor generates a random constraints matrix and random quotas +- providing both a constraints matrix and quotas, in which case the constructor uses them as is. +For quotas, the user can choose between fixed quotas (same for all time steps) or random quotas (different for each time step). + +`over_stock_bound_cost` sets the unit penalty applied to a physical stock outside +`[stock_inf, stock_sup]` (see [`compute_total_cost`](@ref)). It was originally introduced +to force the anticipative baseline to order a little, so it would have something to imitate. +`nothing` (default) keeps that historical behaviour (`maximum(prices)`); pass `0` to disable +the surcost entirely and rely only on `stock_inf`/`stock_sup` as soft targets driving +`virtual_stock_cost`/`physical_stock_cost`. +""" +function DynamicReplenishmentBenchmark(; + N::Int=10, + λ::Int=15, + d::Int=5, + nb_constraints::Int=2, + constraints_matrix=nothing, + quotas=nothing, + stock_inf::Int=0, + stock_sup::Int=30, + ub_same_item::Int=30, + delivery_delay::Int=3, + max_steps::Int=10, + over_stock_bound_cost::Union{Real,Nothing}=nothing, + customer_choice_model=nothing, + prices=nothing, + features=nothing, + seed=nothing, + rng=Xoshiro(seed), +) + if isnothing(constraints_matrix) || isnothing(quotas) + if !isnothing(constraints_matrix) || !isnothing(quotas) + @warn "If either constraints_matrix or quotas is provided, both must be provided. Generating random constraints and quotas." + end + constraints_matrix = rand(rng, 0:1, nb_constraints, N) + quotas = rand(rng, 10:30, max_steps, nb_constraints) + else + @assert size(constraints_matrix, 1) == size(quotas, 2) "The number of constraints in the constraints matrix must match the number of columns in the quotas matrix." + @assert size(constraints_matrix, 2) == N "The number of items ($N) must match the number of columns ($(size(constraints_matrix, 2))) in the constraints matrix." + @assert size(quotas, 1) == max_steps "The number of steps ($max_steps) must match the number of rows ($(size(quotas, 1))) in the quotas matrix." + nb_constraints = size(constraints_matrix, 1) + end + + constraints_matrix = vcat(constraints_matrix, I) + quotas = hcat(quotas, fill(ub_same_item, max_steps, N)) + + if isnothing(prices) + prices = rand(rng, Uniform(1.0, 10.0), N) + else + @assert length(prices) == N "`prices` must have length N=$N, got $(length(prices))." + prices = collect(Float64, prices) + end + if isnothing(features) + features = rand(rng, Uniform(-10.0, 10.0), (d, N)) + else + @assert size(features) == (d, N) "`features` must be of size (d, N)=($d, $N), got $(size(features))." + features = Matrix{Float64}(features) + end + if customer_choice_model === nothing + price_w = rand(rng, Uniform(-1.0, -0.7), 1) + features_w = rand(rng, Uniform(-0.8, -0.1), d) + customer_choice_model = Chain(Dense(reshape(vcat(price_w, features_w), 1, :)), vec) + else + try + customer_choice_model(rand(rng, d + 1, N)) + catch e + throw( + ArgumentError( + "customer_choice_model is incompatible with d=$d (expected input of size (d+1, N)): $e", + ), + ) + end + end + full_features = vcat(prices', features) # (d+1, N) + dt = fit(ZScoreTransform, full_features; dims=2) + scaled_features = transform(dt, full_features) + full_features = scaled_features + static_utilities = customer_choice_model(full_features) + # add no purchase option + static_utilities = vcat(static_utilities, 0.0) + + virtual_stock_cost = prices ./ (max_steps * 10) + physical_stock_cost = prices ./ (max_steps * 5) + over_stock_bound_cost = if isnothing(over_stock_bound_cost) + maximum(prices) + else + Float64(over_stock_bound_cost) + end + max_quotas = Matrix{Int}(undef, max_steps, N) + for i in 1:N, t in 1:max_steps + max_quotas[t, i] = minimum( + quotas[t, c] for + c in axes(constraints_matrix, 1) if constraints_matrix[c, i] == 1 + ) + end + + return DynamicReplenishmentBenchmark{typeof(customer_choice_model)}( + customer_choice_model, + λ, + N, + d, + constraints_matrix, + quotas, + stock_inf, + stock_sup, + ub_same_item, + delivery_delay, + prices, + features, + scaled_features, + dt, + virtual_stock_cost, + physical_stock_cost, + over_stock_bound_cost, + max_steps, + max_quotas, + static_utilities, + ) +end + +# Accessor functions +customer_choice_model(b::DynamicReplenishmentBenchmark) = b.customer_choice_model +poisson_arrival_rate(b::DynamicReplenishmentBenchmark) = b.λ +item_count(b::DynamicReplenishmentBenchmark) = b.N +feature_count(b::DynamicReplenishmentBenchmark) = b.d +max_steps(b::DynamicReplenishmentBenchmark) = b.max_steps +constraints_matrix(b::DynamicReplenishmentBenchmark) = b.constraints_matrix +quotas(b::DynamicReplenishmentBenchmark) = b.quotas +stock_inf(b::DynamicReplenishmentBenchmark) = b.stock_inf +stock_sup(b::DynamicReplenishmentBenchmark) = b.stock_sup +ub_same_item(b::DynamicReplenishmentBenchmark) = b.ub_same_item +delivery_delay(b::DynamicReplenishmentBenchmark) = b.delivery_delay +prices(b::DynamicReplenishmentBenchmark) = b.prices +features(b::DynamicReplenishmentBenchmark) = b.features +scaled_features(b::DynamicReplenishmentBenchmark) = b.scaled_features +feature_transform(b::DynamicReplenishmentBenchmark) = b.feature_transform +virtual_stock_cost(b::DynamicReplenishmentBenchmark) = b.virtual_stock_cost +physical_stock_cost(b::DynamicReplenishmentBenchmark) = b.physical_stock_cost +over_stock_bound_cost(b::DynamicReplenishmentBenchmark) = b.over_stock_bound_cost +nb_constraints(b::DynamicReplenishmentBenchmark) = size(b.constraints_matrix, 1) +max_quotas(b::DynamicReplenishmentBenchmark) = b.max_quotas + +# The objective is a margin net of stock costs, maximized. +Utils.is_minimization_problem(::DynamicReplenishmentBenchmark) = false + +include("utils.jl") + +include("state.jl") +include("scenario.jl") +include("environment.jl") +include("statistical_model.jl") +include("policies.jl") +include("maximizer.jl") +include("anticipative_solver.jl") +include("features.jl") + +""" +$TYPEDSIGNATURES + +Creates a random environment for the dynamic replenishment benchmark using the provided random number generator. +""" +function Utils.build_environment( + b::DynamicReplenishmentBenchmark, + rng::AbstractRNG; + stock_ini_fill_rate=nothing, + kwargs..., +) + return if isnothing(stock_ini_fill_rate) + Environment(b, rng) + else + Environment(b, rng; stock_ini_fill_rate) + end +end + +""" +$TYPEDSIGNATURES + +Rebuild an environment sitting at `sample.state` and running on `scenario`, so that the +callable returned by [`Utils.generate_parametric_anticipative_solver`](@ref) can be applied +at any epoch of a stored trajectory. + +The state is shared with `sample`, not copied: the solvers reading it must not mutate it. +""" +function Utils.build_environment( + ::DynamicReplenishmentBenchmark, sample::DataSample, scenario::Scenario +) + hasproperty(sample, :state) || error( + "`build_environment` needs the epoch state in `sample.state`, got a DataSample " * + "with $(propertynames(sample)).", + ) + state = sample.state + state isa DRPState || + error("`sample.state` should be a `DRPState`, got a $(typeof(state)) instead.") + return Environment(; + config=state.config, state=state, scenario=scenario, stock_ini=stock_ini(state) + ) +end + +""" +$TYPEDEF + +Callable wrapping [`replenishment_problem`](@ref) with a fixed `model_builder`, so it can +be passed to `LinearMaximizer` without a closure. +""" +struct MaximizerProblem{M} + model_builder::M +end +function (p::MaximizerProblem)(Θ; kwargs...) + return replenishment_problem(Θ; kwargs..., model_builder=p.model_builder) +end + +function Utils.generate_maximizer( + ::DynamicReplenishmentBenchmark; model_builder=highs_model +) + return LinearMaximizer(MaximizerProblem(model_builder); g) +end + +""" +$TYPEDEF + +Callable wrapping [`anticipative_solver`](@ref) with a fixed `model_builder`, returned by +[`Utils.generate_anticipative_solver`](@ref). +""" +struct AnticipativeSolverCall{M} + model_builder::M + "relative MIP gap tolerance" + mip_gap::Float64 + "solver time limit in seconds, `nothing` to solve to optimality" + time_limit::Union{Float64,Nothing} +end + +function AnticipativeSolverCall( + model_builder::M; mip_gap::Real=0.0, time_limit::Union{Real,Nothing}=nothing +) where {M} + return AnticipativeSolverCall{M}( + model_builder, + Float64(mip_gap), + isnothing(time_limit) ? nothing : Float64(time_limit), + ) +end + +function (s::AnticipativeSolverCall)( + env::Utils.SeededEnvironment; reset_env=false, kwargs... +) + _, trajectory = anticipative_solver( + env.env, + env.rng; + reset_env, + mip_gap=s.mip_gap, + time_limit=s.time_limit, + kwargs..., + model_builder=s.model_builder, + ) + return trajectory +end + +function Utils.generate_anticipative_solver( + ::DynamicReplenishmentBenchmark; + model_builder=highs_model, + mip_gap::Real=0.0, + time_limit::Union{Real,Nothing}=nothing, +) + return AnticipativeSolverCall(model_builder; mip_gap, time_limit) +end + +""" +$TYPEDEF + +Callable wrapping [`anticipative_solver`](@ref) (scenario-conditioned) with a fixed +`model_builder`, returned by [`Utils.generate_parametric_anticipative_solver`](@ref). +""" +struct ParametricAnticipativeSolverCall{M} + model_builder::M + "relative MIP gap tolerance" + mip_gap::Float64 + "solver time limit in seconds, `nothing` to solve to optimality" + time_limit::Union{Float64,Nothing} +end + +function ParametricAnticipativeSolverCall( + model_builder::M; mip_gap::Real=0.0, time_limit::Union{Real,Nothing}=nothing +) where {M} + return ParametricAnticipativeSolverCall{M}( + model_builder, + Float64(mip_gap), + isnothing(time_limit) ? nothing : Float64(time_limit), + ) +end + +function (s::ParametricAnticipativeSolverCall)( + θ, scenario::Scenario, env::Utils.SeededEnvironment; reset_env=true, kwargs... +) + # reset_env && Utils.reset_to_initial!(env) + _, trajectory = anticipative_solver( + env.env, + env.rng, + scenario; + reset_env=false, + θ, + mip_gap=s.mip_gap, + time_limit=s.time_limit, + kwargs..., + model_builder=s.model_builder, + ) + return trajectory +end + +function Utils.generate_parametric_anticipative_solver( + ::DynamicReplenishmentBenchmark; + model_builder=highs_model, + mip_gap::Real=0.0, + time_limit::Union{Real,Nothing}=nothing, +) + return ParametricAnticipativeSolverCall(model_builder; mip_gap, time_limit) +end + +""" +$TYPEDSIGNATURES + +Returns baseline policies for the dynamic replenishment benchmark: `Greedy`, `Random`, +`Lazy`, `MeanAnticipative`, `MeanAnticipativePerEpoch` and `SAA`. + +Both mean policies need `anticipative_results` (anticipative decisions), if empty, they fall back to `Lazy`. +Remaining keyword arguments go to the SAA +policy. +""" +function Utils.generate_baseline_policies( + ::DynamicReplenishmentBenchmark; + model_builder=highs_model, + anticipative_results::AbstractVector{<:DataSample}=DataSample[], + order_item::Function=mean_feature_order, + kwargs..., +) + greedy = Policy{DynamicReplenishmentBenchmark}( + "Greedy", "policy that replenishes items in decreasing price order", greedy_policy + ) + random = Policy{DynamicReplenishmentBenchmark}( + "Random", + "Policy that replenishes items in a random order with random quantities", + random_policy, + ) + lazy = Policy{DynamicReplenishmentBenchmark}( + "Lazy", "Policy that replenishes nothing", lazy_policy + ) + mean_anticipative = Policy{DynamicReplenishmentBenchmark}( + "MeanAnticipative", + "Policy that replenishes items in increasing mean feature order, with quantities equal to the mean of the anticipative results over the whole horizon", + MeanAnticipativePolicyCall(anticipative_results, order_item), + ) + mean_anticipative_per_epoch = Policy{DynamicReplenishmentBenchmark}( + "MeanAnticipativePerEpoch", + "Policy that replenishes items in increasing mean feature order, with quantities equal to the mean of the anticipative results at the current epoch", + MeanAnticipativePolicyCall(anticipative_results, order_item; per_epoch=true), + ) + saa = Policy{DynamicReplenishmentBenchmark}( + "SAA", + "Policy that solves a sample average approximation problem.", + SAAPolicyCall(model_builder; kwargs...), + ) + return (; greedy, random, lazy, mean_anticipative, mean_anticipative_per_epoch, saa) +end + +export DynamicReplenishmentBenchmark + +end diff --git a/src/DynamicReplenishment/anticipative_solver.jl b/src/DynamicReplenishment/anticipative_solver.jl new file mode 100644 index 00000000..f7c60b6d --- /dev/null +++ b/src/DynamicReplenishment/anticipative_solver.jl @@ -0,0 +1,434 @@ + +""" +$TYPEDSIGNATURES + +Compute big M values for a scenario of a specific environment. +""" +function compute_bigM_sales(env::Environment, scenario::Scenario) + T = max_steps(env.config) + N = item_count(env.config) + max_q = max_quotas(env.config) + s0 = stock(env) + n_customers = nb_customers(scenario) + utilities = scenario.utilities + big_M = Vector{Vector{Vector{Int}}}(undef, T - current_epoch(env) + 1) + for (t_m, t) in enumerate(current_epoch(env):T) + big_M[t_m] = Vector{Vector{Int}}(undef, n_customers[t]) + for k in 1:n_customers[t] + big_M[t_m][k] = zeros(Int, N + 1) + sorted_indices = sortperm(utilities[t][k]) # ascending order + no_buy_index = findfirst(==(N + 1), sorted_indices) + for (index, i_1) in enumerate(sorted_indices[1:(end - 1)]) + if index >= no_buy_index + higher_items = [ + i_2 for i_2 in sorted_indices[(index + 1):end] if i_2 <= N + ] + ini_stock_sum = sum(s0[i_2] for i_2 in higher_items) + quota_sum = sum(max_q[τ, i_2] for τ in 1:t for i_2 in higher_items) + # M = ∑_τ=1^t ∑_{i_2: u_{i_2} > u_{i_1}} max_q[τ][i_2] + stock_ini[i_2] + 1 + big_M[t_m][k][i_1] = quota_sum + ini_stock_sum + 1 + end + end + end + end + return big_M +end + +function compute_bigM_physical_stock(env::Environment, scenario::Scenario) + T = max_steps(env.config) - current_epoch(env) + 1 + N = item_count(env.config) + delay = delivery_delay(env.config) + max_q = max_quotas(env.config)[current_epoch(env):end, :] + s0 = stock(env) + n_customer = nb_customers(scenario)[current_epoch(env):end] # borne des ventes + + big_M = zeros(Int, N, T + 1) + for i in 1:N, t in 2:(T + 1) + t_arrived = max(0, t - delay) + pos = s0[i] + sum(max_q[τ, i] for τ in 1:t_arrived; init=0) # borne x_hi + neg = sum(n_customer[1:(t - 1)]) # borne -x_lo (≤ ventes cumulées) + big_M[i, t] = max(pos, neg) + end + return big_M +end + +function stock_constraints!(m, y, s, α, T, N, nb_customers, stock_ini) + # Initial stock + @constraint(m, [i in 1:N], s[1, i] == stock_ini[i]) + # Stock dynamics + @constraint( + m, + [i in 1:N, t in 1:T], + s[t + 1, i] == s[t, i] + y[t, i] - sum(α[i, t, k] for k in 1:nb_customers[t]) + ) + return nothing +end + +function customer_constraints!(m, α, T, N, nb_customers) + # Each customer buys at most one vehicle (no purchase option included) + @constraint( + m, [t in 1:T, k in 1:nb_customers[t]], sum(α[i, t, k] for i in 1:(N + 1)) == 1 + ) + return nothing +end + +function sales_order_constraints!(m, y, s, α, T, N, nb_customers, utilities, bigM_s) + for t in 1:T + for k in 1:nb_customers[t] + sorted_indices = sortperm(utilities[t][k]) # ascending order + no_buy_index = findfirst(==(N+1), sorted_indices) + for (index, i_1) in enumerate(sorted_indices[1:(end - 1)]) + # no-buy case + if index < no_buy_index + @constraint(m, α[i_1, t, k] == 0) + continue + else + # don't sell i_1 if ∃ i_2 in stock s.t. u_{i_2} > u_{i_1} + if k == 1 + @constraint( + m, + α[i_1, t, k] <= ( + 1 - + sum( + s[t, i_2] + y[t, i_2] for + i_2 in sorted_indices[(index + 1):end] if i_2 <= N + ) / bigM_s[t][k][i_1] + ), + ) + else + @constraint( + m, + α[i_1, t, k] <= ( + 1 - + sum( + s[t, i_2] + y[t, i_2] - + sum(α[i_2, t, j] for j in 1:(k - 1)) for + i_2 in sorted_indices[(index + 1):end] if i_2 <= N + ) / bigM_s[t][k][i_1] + ), + ) + end + end + end + end + end + return nothing +end + +""" +$TYPEDSIGNATURES + +Add quota constraints to the model. +""" +function quota_constraints!(m, y, T, N, constraints_matrix, quotas) + nb_cons = size(constraints_matrix, 1) + @constraint( + m, + [c in 1:(nb_cons), t in 1:T], + sum(constraints_matrix[c, i] * y[t, i] for i in 1:N) <= quotas[t, c] + ) + return nothing +end +""" +$TYPEDSIGNATURES + +Add physical stock constraints (exact linearization of v = max(0, x) via indicator binaries). +""" +function physical_stock_constraints!( + m, y, α, v, z, T, N, delivery_delay, stock_ini, nb_customers, bigM_ps +) + @constraint(m, [i in 1:N], v[1, i] == stock_ini[i]) + + # v = max(0, x_phys) via indicatrice z (z=1 ⟺ x_phys ≥ 0) + @constraint( + m, + [i in 1:N, t in 2:(T + 1)], + v[t, i] >= + stock_ini[i] + sum(y[τ, i] for τ in 1:(t - delivery_delay); init=zero(AffExpr)) - + sum(α[i, τ, k] for τ in 1:(t - 1) for k in 1:nb_customers[τ]; init=zero(AffExpr)) + ) + @constraint( + m, + [i in 1:N, t in 2:(T + 1)], + v[t, i] <= + stock_ini[i] + sum(y[τ, i] for τ in 1:(t - delivery_delay); init=zero(AffExpr)) - + sum(α[i, τ, k] for τ in 1:(t - 1) for k in 1:nb_customers[τ]; init=zero(AffExpr)) + + bigM_ps[i, t] * (1 - z[i, t]) + ) + @constraint(m, [i in 1:N, t in 2:(T + 1)], v[t, i] <= bigM_ps[i, t] * z[i, t]) + return nothing +end + +""" +$TYPEDSIGNATURES + +Add stock bounds constraints. +""" +function stock_bounds_constraints!(m, v, T, N, s_min, s_sup, stock_inf, stock_sup) + @constraint(m, [t in 1:T], s_min[t] >= stock_inf - sum(v[t + 1, i] for i in 1:N)) + @constraint(m, [t in 1:T], s_sup[t] >= sum(v[t + 1, i] for i in 1:N) - stock_sup) + return nothing +end + +""" +$TYPEDSIGNATURES + +Compute the base objective function. +""" +function compute_objective(y, s, α, v, T, s_min, s_sup, env, nb_customers) + N = item_count(env) + # margin + # `init`: a scenario can have no customer at all over the horizon (small λ or short + # horizon), in which case the inner sum is empty and would throw without it. + margin = sum( + prices(env)[i] * + sum(α[i, t, k] for t in 1:T for k in 1:nb_customers[t]; init=zero(AffExpr)) for + i in 1:N; + init=zero(AffExpr), + ) + # virtual stock cost + virtual_stock = sum(virtual_stock_cost(env)[i] * s[t + 1, i] for t in 1:T for i in 1:N) + # physical stock cost + physical_stock = sum( + physical_stock_cost(env)[i] * v[t + 1, i] for t in 1:T for i in 1:N + ) + # over bound stock + under_stock_min = sum(s_min) + over_stock_sup = sum(s_sup) + + objective = + margin - virtual_stock - physical_stock - + over_stock_bound_cost(env) * (under_stock_min + over_stock_sup) + + return objective +end + +function solver_variable_to_dataset( + env::Environment, + scenario::Scenario, + s_val, + y_val, + α_val, + v_val, + obj_val; + θ=nothing, + κ=1.0, + state::DRPState=env.state, + mip_gap=nothing, +) + s_val = Int.(round.(s_val)) # (T+1, N) + y_val = Int.(round.(y_val)) # (T, N) + α_val = Int.(round.(α_val)) # (N+1, T, k) + v_val = Int.(round.(v_val)) # (T+1, N) + + config = env.config + T = max_steps(config) - current_epoch(env) + 1 + N = item_count(config) + n_customers = nb_customers(scenario)[current_epoch(env):end] + max_q = max_quotas(config)[current_epoch(env):end, :] + # sales_full[t, i] = total units of item i sold at epoch t + sales_full = zeros(Int, T, N) + for t in 1:T, i in 1:N + sales_full[t, i] = sum( + round(Int, value(α_val[i, t, k])) for k in 1:n_customers[t]; init=0 + ) + end + dataset = Vector{DataSample}(undef, T) + + # initial state, before any replenishment/sales (epoch 0 / pre-action). + init_state = deepcopy(state) + x_init = compute_features(init_state) + y_init = y_val[1, :] + dataset[1] = DataSample(; + y=y_init, + x=x_init, + state=init_state, + extra=(; next_sales=sales_full[1, :], customers=n_customers[1]), + ) + for t in 2:T + state_t = DRPState(; + config=config, + current_epoch=t, + stock=s_val[t, :], + stock_history=s_val[1:t, :], + replenishment_history=y_val[1:(t - 1), :], + sales_history=sales_full[1:(t - 1), :], + customer_history=n_customers[1:(t - 1)], + ub_per_item=s_val[t, :] .+ max_q[t, :], + ) + y_true = y_val[t, :] + x = compute_features(state_t) + dataset[t] = DataSample(; + y=y_true, + x, + state=state_t, + extra=(; next_sales=sales_full[t, :], customers=n_customers[t]), + ) + end + final_state = DRPState(; + config=config, + current_epoch=T + 1, + stock=s_val[T + 1, :], + stock_history=s_val[1:(T + 1), :], + replenishment_history=y_val[1:T, :], + sales_history=sales_full[1:T, :], + customer_history=n_customers[1:T], + ub_per_item=s_val[T + 1, :] .+ max_q[end, :], + ) + final_obj_val = total_cost(final_state) + + if !isnothing(θ) + g_y = g(dataset[1].y; state=dataset[1].state) + @assert length(θ) == N + sum(ub_per_item(dataset[1].state)) + final_obj_val += κ * dot(θ, g_y) + end + if !isapprox(obj_val, final_obj_val, atol=1e-3, rtol=1e-3) + # Écart entre l'objectif rapporté par le solveur et l'objectif recalculé + # depuis la trajectoire arrondie : signe d'instabilité numérique du solveur + # (le même MILP produit les warnings SCIP "LP solution value is above SCIP's + # infinity value" même sans θ). On logue au lieu de planter — la trajectoire + # elle-même reste utilisable — pour avoir enfin les chiffres la prochaine + # fois que ça se produit. + abs_diff = abs(obj_val - final_obj_val) + rel_diff = abs_diff / max(abs(final_obj_val), 1.0) + epoch = current_epoch(env) + max_abs_theta = isnothing(θ) ? nothing : maximum(abs, θ) + nonfinite_theta = isnothing(θ) ? nothing : count(!isfinite, θ) + @warn "Anticipatif paramétrique : objectif solveur et objectif recalculé en désaccord (instabilité numérique probable) — trajectoire conservée telle quelle" obj_val final_obj_val abs_diff rel_diff mip_gap epoch max_abs_theta nonfinite_theta maxlog = + 20 + end + return dataset +end + +""" +$TYPEDSIGNATURES + +Construct yη vector for +""" +function g_model(m, N, ub, y, s) + # Same encoding as `replenishment_problem`: z is the staircase indicator of the + # stock level, z[i, j] = 1 iff j <= s[i] + y[i]. + @variable(m, z_eta[i in 1:N, j in 1:ub[i]], Bin) + + @constraint(m, [i in 1:N], sum(z_eta[i, j] for j in 1:ub[i]) == s[i] + y[i]) + @constraint(m, [i in 1:N, j in 1:(ub[i] - 1)], z_eta[i, j] >= z_eta[i, j + 1]) + + y_eta_vec = Vector{AffExpr}(undef, sum(ub)) + row = 1 + for i in 1:N + y_eta_vec[row] = 1 * z_eta[i, 1] + for k in 2:ub[i] + # max(0, s[i] + y[i] - (k - 1)) = number of levels j >= k that are filled + y_eta_vec[row + k - 1] = -sum(z_eta[i, j] for j in k:ub[i]) + end + row += ub[i] + end + return vcat(vec(y), y_eta_vec) +end + +""" +$TYPEDSIGNATURES + +Solve the anticipative problem for a given instance and scenario. +""" +function anticipative_solver( + env::Environment, + rng::AbstractRNG, + scenario::Scenario=env.scenario; + model_builder=highs_model, + reset_env::Bool=true, + verbose::Bool=false, + bigM_s=nothing, + bigM_ps=nothing, + θ=nothing, + state::DRPState=env.state, + κ::Float64=1.0, + mip_gap::Float64=0.0, + time_limit::Union{Real,Nothing}=nothing, +) + if reset_env + reset!(env, rng) + scenario = env.scenario + state = env.state + end + + if isnothing(bigM_s) + bigM_s = compute_bigM_sales(env, scenario) + end + if isnothing(bigM_ps) + bigM_ps = compute_bigM_physical_stock(env, scenario) + end + + @assert !is_terminated(env) + + m = model_builder() + verbose || set_silent(m) + set_attribute(m, MOI.RelativeGapTolerance(), mip_gap) + isnothing(time_limit) || set_attribute(m, MOI.TimeLimitSec(), Float64(time_limit)) + N = item_count(env) + T = max_steps(env) - current_epoch(env) + 1 + n_customers = nb_customers(scenario)[current_epoch(env):end] + q = quotas(env)[current_epoch(env):end, :] + s0 = stock(env) + ## Variables + @variable(m, y[1:T, 1:N] >= 0, Int) # replenishments + @variable(m, s[1:(T + 1), 1:N] >= 0, Int) # stock + @variable(m, α[i in 1:(N + 1), t in 1:T, k in 1:n_customers[t]], Bin) # sales + @variable(m, v[1:(T + 1), 1:N] >= 0, Int) # physical stock + @variable(m, z[i in 1:N, t in 2:(T + 1)], Bin) # auxiliary binary for physical stock linearization + @variable(m, s_min[1:T] >= 0, Int) # stock under min + @variable(m, s_sup[1:T] >= 0, Int) # stock over max + + ## Constraints + stock_constraints!(m, y, s, α, T, N, n_customers, s0) + customer_constraints!(m, α, T, N, n_customers) + sales_order_constraints!( + m, y, s, α, T, N, n_customers, scenario.utilities[current_epoch(env):end], bigM_s + ) + quota_constraints!(m, y, T, N, constraints_matrix(env), q) + physical_stock_constraints!( + m, y, α, v, z, T, N, delivery_delay(env), s0, n_customers, bigM_ps + ) + stock_bounds_constraints!(m, v, T, N, s_min, s_sup, stock_inf(env), stock_sup(env)) + + ## Objective + objective = compute_objective(y, s, α, v, T, s_min, s_sup, env, n_customers) + if !isnothing(θ) + if !all(isfinite, θ) || maximum(abs, θ) > 1e12 + @warn "Solveur anticipatif paramétré : θ hors échelle en entrée" max_abs_theta = maximum( + abs, θ + ) nonfinite_theta = count(!isfinite, θ) epoch = current_epoch(env) maxlog = 10 + end + g_y = g_model(m, N, ub_per_item(state), y[1, :], s[1, :]) + @assert length(θ) == N + sum(ub_per_item(state)) + objective += κ * dot(θ, g_y) + end + @objective(m, Max, objective) + + optimize!(m) + if primal_status(m) == MOI.FEASIBLE_POINT + obj_val = objective_value(m) + if !isnothing(θ) && (!isfinite(obj_val) || abs(obj_val) >= 1e20) + @warn "Parametric Anticipatif Solver: infinite objective" objective = obj_val max_abs_theta = maximum( + abs, θ + ) nonfinite_theta = count(!isfinite, θ) epoch = current_epoch(env) maxlog = 10 + end + dataset = solver_variable_to_dataset( + env, + scenario, + value.(s), + value.(y), + value.(α), + value.(v), + obj_val; + θ=θ, + κ=κ, + state=state, + mip_gap=mip_gap, + ) + return obj_val, dataset + else + @warn("No feasible points found.") + return nothing, nothing + end +end diff --git a/src/DynamicReplenishment/environment.jl b/src/DynamicReplenishment/environment.jl new file mode 100644 index 00000000..5d2fafd8 --- /dev/null +++ b/src/DynamicReplenishment/environment.jl @@ -0,0 +1,147 @@ +""" +$TYPEDEF + +Environment for the Dynamic Replenishment problem. + +# Fields +$TYPEDFIELDS +""" +@kwdef mutable struct Environment{B<:DynamicReplenishmentBenchmark,S<:DRPState} <: + Utils.AbstractEnvironment + "associated benchmark" + config::B + "current state" + state::S + "scenario the environment will use when not given a specific one" + scenario::Scenario + "initial stock" + stock_ini::Vector{Int} +end + +# Accessor functions +customer_choice_model(env::Environment) = customer_choice_model(env.config) +poisson_arrival_rate(env::Environment) = poisson_arrival_rate(env.config) +item_count(env::Environment) = item_count(env.config) +feature_count(env::Environment) = feature_count(env.config) +max_steps(env::Environment) = max_steps(env.config) +constraints_matrix(env::Environment) = constraints_matrix(env.config) +quotas(env::Environment) = quotas(env.config) +stock_inf(env::Environment) = stock_inf(env.config) +stock_sup(env::Environment) = stock_sup(env.config) +ub_same_item(env::Environment) = ub_same_item(env.config) +delivery_delay(env::Environment) = delivery_delay(env.config) +prices(env::Environment) = prices(env.config) +features(env::Environment) = features(env.config) +virtual_stock_cost(env::Environment) = virtual_stock_cost(env.config) +physical_stock_cost(env::Environment) = physical_stock_cost(env.config) +over_stock_bound_cost(env::Environment) = over_stock_bound_cost(env.config) +max_quotas(env::Environment) = max_quotas(env.config) + +current_epoch(env::Environment) = current_epoch(env.state) +stock_ini(env::Environment) = env.stock_ini +stock(env::Environment) = stock(env.state) +ub_per_item(env::Environment) = ub_per_item(env.state) + +""" +$TYPEDSIGNATURES + +Draw an initial stock with a fill_rate relative to the stock_sup, spread +uniformly at random over the items. +""" +function draw_stock_ini(rng::AbstractRNG, N::Int, stock_sup::Int, fill_rate::Real) + 0 <= fill_rate <= 1 || + throw(ArgumentError("`stock_ini_fill_rate` must be in [0, 1], got $fill_rate.")) + stock_ini = zeros(Int, N) + for _ in 1:round(Int, fill_rate * stock_sup) + stock_ini[rand(rng, 1:N)] += 1 + end + return stock_ini +end + +""" +$TYPEDSIGNATURES + +Creates an [`Environment`](@ref) from an instance of the dynamic replenishment benchmark. +The initial stock fills half of `stock_sup` by default, see [`draw_stock_ini`](@ref). Pass +`stock_ini` directly to bypass the random draw entirely (e.g. all-zero stock). +""" +function Environment( + config::DynamicReplenishmentBenchmark, + rng::AbstractRNG; + stock_ini_fill_rate::Real=0.5, + stock_ini=draw_stock_ini( + rng, item_count(config), stock_sup(config), stock_ini_fill_rate + ), +) + scenario = Utils.generate_scenario(config; rng=rng) + initial_state = DRPState(config, stock_ini) + return Environment(; config, state=initial_state, scenario, stock_ini) +end + +function Environment( + config::DynamicReplenishmentBenchmark, + scenario::Scenario, + rng::AbstractRNG; + stock_ini_fill_rate::Real=0.5, + stock_ini=draw_stock_ini( + rng, item_count(config), stock_sup(config), stock_ini_fill_rate + ), +) + initial_state = DRPState(config, stock_ini) + return Environment(; config, state=initial_state, scenario, stock_ini) +end + +""" +$TYPEDSIGNATURES + +Get the current state of the environment. +""" +function Utils.observe(env::Environment) + return compute_features(env.state), env.state +end + +""" +$TYPEDSIGNATURES + +Check if the episode is terminated, i.e. if the current epoch is the last one. +""" +Utils.is_terminated(env::Environment) = current_epoch(env) > max_steps(env) + +""" +$TYPEDSIGNATURES + +Reset the environment to its initial state. +Also reset the rng to `seed` if `reset_rng` is set to true. +""" +function Utils.reset!(env::Environment, rng::AbstractRNG; reset_scenario::Bool=true) + if reset_scenario + env.scenario = Utils.generate_scenario(env.config; rng) + end + reset_state!(env.state, rng; reset_stock_ini=false) + return nothing +end + +""" +$TYPEDSIGNATURES + +Apply the replenishment to the stock, apply the sales and increase time. +""" +function Utils.step!(env::Environment, replenishment, rng::AbstractRNG) + @assert !Utils.is_terminated(env) "Environment is terminated, cannot act!" + apply_replenishment!(env.state, replenishment) + delta_cost = apply_sales!(env.state; env.scenario[current_epoch(env)]...) + add_customers!(env.state; env.scenario[current_epoch(env)]...) + env.state.current_epoch += 1 + if !is_terminated(env) + env.state.ub_per_item = + env.state.stock .+ max_quotas(env.config)[current_epoch(env), :] + env.state.physical_stock = compute_physical_stock( + env.config, + current_epoch(env), + stock_ini(env.state), + env.state.replenishment_history, + env.state.sales_history, + ) + end + return delta_cost +end diff --git a/src/DynamicReplenishment/features.jl b/src/DynamicReplenishment/features.jl new file mode 100644 index 00000000..cf098c1a --- /dev/null +++ b/src/DynamicReplenishment/features.jl @@ -0,0 +1,258 @@ +function mean_feature_matrix(state::DRPState, feature_matrix) + N = item_count(state.config) + @assert size(feature_matrix, 2) == N + isempty(feature_matrix) && return zeros(Float64, N) + return vec(mean(feature_matrix; dims=1)) +end + +mean_sales_history(state::DRPState) = mean_feature_matrix(state, sales_history(state)) +function mean_replenishment_history(state::DRPState) + return mean_feature_matrix(state, replenishment_history(state)) +end +mean_stock_history(state::DRPState) = mean_feature_matrix(state, stock_history(state)) + +function items_with_positive_stock(state::DRPState) + isempty(stock_history(state)) && return Int[] + return [ + i for i in 1:item_count(state.config) if sum(view(stock_history(state), :, i)) > 0 + ] +end + +""" +Return mean features of item if it was ever in stock, +return mean features across all items that have had positive stock otherwise. +""" +function item_mean_feature( + feature_row_means::Vector{Float64}, item::Int, positive_stock_items::Vector{Int} +) + isempty(positive_stock_items) && return 0.0 + if item ∈ positive_stock_items + return feature_row_means[item] + else + return mean_or_zero(feature_row_means[positive_stock_items]) + end +end + +""" +Compute the Days On Lot of the items. +It corresponds to the number of time step an item stays physically in stock before being sold. +It can be negative if the item is sold before arriving. +""" +function compute_dol_item(state::DRPState, item::Int) + replenishments = replenishment_history(state) # (t, N) + sales = sales_history(state) # (t, N) + s0 = stock_ini(state)[item] + t_now = current_epoch(state) + + repl_item = isempty(replenishments) ? Int[] : replenishments[:, item] + sales_item = isempty(sales) ? Int[] : sales[:, item] + + cum_arrivals = s0 .+ cumsum(repl_item) # cum_arrivals[t] = total arrived by end of epoch t + total_nb_item = s0 + sum(repl_item) + + total_nb_item == 0 && return Float64[] + + cum_sales = cumsum(sales_item) # cum_sales[t] = total sold by end of epoch t + total_nb_sales = isempty(cum_sales) ? 0 : cum_sales[end] + + if total_nb_sales == 0 + return fill(Float64(t_now), total_nb_item) + end + + dols = zeros(Float64, total_nb_item) + for j in 1:total_nb_item + date_repl_j = if j <= s0 + 1 + else + t = findfirst(>=(j), cum_arrivals) + t === nothing ? t_now : t + 1 # arrives the epoch *after* the replenishment lands + end + + end_date_j = if j <= total_nb_sales + something(findfirst(>=(j), cum_sales), t_now) + else + t_now + end + + dols[j] = end_date_j - date_repl_j + 1 + end + return dols +end + +""" +Number of dynamic (state-dependent) columns appended per item. +""" +nb_dynamic_item_features(config) = 17 + +""" +Number of rows of the item block, i.e. the input size of the `θ` model. +""" +item_features_size(config) = feature_count(config) + 1 + nb_dynamic_item_features(config) + +""" +Number of stock-level columns appended by [`create_stock_features`](@ref). +""" +const NB_STOCK_FEATURES = 12 + +""" +Number of rows of the stock block, i.e. the input size of the `η` model. +Be careful: the feature matrix has one additional row, the item identifier. +""" +stock_features_size(config) = item_features_size(config) + NB_STOCK_FEATURES + +""" +$TYPEDSIGNATURES + +Create features per item. +The first `feature_count(config) + 1` columns correspond to static features (scaled price +and item features). The remaining [`nb_dynamic_item_features`](@ref) columns are dynamic: + +- current *virtual* stock, and scaled with price +- mean sales, and scaled with price +- mean stock, and scaled with price +- mean number of customers in the past +- mean days on lot, and scaled with price +- current *physical* stock, and scaled with price +- stock in transit (`stock - physical_stock`), and scaled with price +- four state-level columns, identical for every item: the total physical stock, the + slack to `stock_inf` and to `stock_sup`, and the remaining horizon +""" +function create_items_features(state::DRPState) + config = state.config + N = item_count(config) + nb_static = feature_count(config) + 1 # replaces instance.nb_features + nb_features = item_features_size(config) + item_features = zeros(Float32, N, nb_features) + + # precompute once + pos_items = items_with_positive_stock(state) + mean_sales = mean_sales_history(state) + mean_stock = mean_stock_history(state) + current_stock = stock(state) + phys_stock = physical_stock(state) + static_features = scaled_features(config) + + # state-level quantities, shared by every item + total_physical = sum(phys_stock) + inf_slack = total_physical - stock_inf(config) + sup_slack = stock_sup(config) - total_physical + remaining_horizon = max_steps(config) - current_epoch(state) + + for i in 1:N + p = prices(config)[i] + ## static features + item_features[i, 1:nb_static] = static_features[:, i] + ## current total stock (virtual + physical) + item_features[i, nb_static + 1] = current_stock[i] + item_features[i, nb_static + 2] = current_stock[i] * p + ## mean sales + ms = item_mean_feature(mean_sales, i, pos_items) + item_features[i, nb_static + 3] = ms + item_features[i, nb_static + 4] = ms * p + ## mean stock + mst = item_mean_feature(mean_stock, i, pos_items) + item_features[i, nb_static + 5] = mst + item_features[i, nb_static + 6] = mst * p + ## mean customers in the past + item_features[i, nb_static + 7] = mean_or_zero(state.customer_history) + ## dol item_features + dols = compute_dol_item(state, i) + item_features[i, nb_static + 8] = mean_or_zero(dols) + item_features[i, nb_static + 9] = mean_or_zero(dols) * p + ## physical stock + item_features[i, nb_static + 10] = phys_stock[i] + item_features[i, nb_static + 11] = phys_stock[i] * p + ## stock in transit + in_transit = current_stock[i] - phys_stock[i] + item_features[i, nb_static + 12] = in_transit + item_features[i, nb_static + 13] = in_transit * p + ## state-level features (same for all items) + item_features[i, nb_static + 14] = total_physical + item_features[i, nb_static + 15] = inf_slack + item_features[i, nb_static + 16] = sup_slack + item_features[i, nb_static + 17] = remaining_horizon + end + return item_features +end + +""" +$TYPEDSIGNATURES + +Create features per stock level per archetype. Row `(i, j)` describes the candidate +post-decision (virtual) stock level `j` for item `i`. + +The first `item_features_size(config)` columns repeat the item features, the next +[`NB_STOCK_FEATURES`](@ref) are dynamic stock features, and the last one is the item +identifier used by [`StatisticalModel`](@ref): + +- deviation from `stock_inf` and scaled with price +- deviation from `stock_sup` and scaled with price +- deviation from mean stock and scaled with price +- deviation from `max_quotas` and scaled with price +- *coupled* deviation from `stock_inf` and `stock_sup`, i.e. what the total physical stock would be if item + `i` were at level `j`, and scaled with price +""" +function create_stock_features(state::DRPState, item_features::Matrix{Float32}) + config = state.config + N = item_count(config) + ub = ub_per_item(state) + nb_fi = size(item_features, 2) + total_rows = sum(ub) + stock_features = zeros(Float32, total_rows, nb_fi + NB_STOCK_FEATURES + 1) # +1 for unique index + t = current_epoch(state) + + pos_items = items_with_positive_stock(state) + mean_stock = mean_stock_history(state) + phys_stock = physical_stock(state) + total_physical = sum(phys_stock) + + stock_inf = config.stock_inf + stock_sup = config.stock_sup + + starts = [1; cumsum(ub)[1:(end - 1)] .+ 1] + ends = cumsum(ub) + + for i in 1:N + rows = starts[i]:ends[i] + stock_features[rows, 1:nb_fi] .= item_features[i:i, :] + + p = prices(config)[i] + js = 1:ub[i] + stock_inf_dev = js .- stock_inf + stock_sup_dev = stock_sup .- js + stock_mean_dev = js .- item_mean_feature(mean_stock, i, pos_items) + max_quotas_dev = max_quotas(state.config)[t, i] .- js + # Total stock if item i ended at level j + others_physical = total_physical - phys_stock[i] + total_inf_dev = (others_physical .+ js) .- stock_inf + total_sup_dev = stock_sup .- (others_physical .+ js) + + stock_features[rows, nb_fi + 1] = stock_inf_dev + stock_features[rows, nb_fi + 2] = stock_inf_dev .* p + stock_features[rows, nb_fi + 3] = stock_sup_dev + stock_features[rows, nb_fi + 4] = stock_sup_dev .* p + stock_features[rows, nb_fi + 5] = stock_mean_dev + stock_features[rows, nb_fi + 6] = stock_mean_dev .* p + stock_features[rows, nb_fi + 7] = max_quotas_dev + stock_features[rows, nb_fi + 8] = max_quotas_dev .* p + stock_features[rows, nb_fi + 9] = total_inf_dev + stock_features[rows, nb_fi + 10] = total_inf_dev .* p + stock_features[rows, nb_fi + 11] = total_sup_dev + stock_features[rows, nb_fi + 12] = total_sup_dev .* p + stock_features[rows, nb_fi + NB_STOCK_FEATURES + 1] .= i # identifier for the item for the statistical model + end + return stock_features +end + +""" +$TYPEDSIGNATURES + +Create features from state. +""" +function compute_features(state::DRPState) + # archetype features + item_features = create_items_features(state) + # stock features + stock_features = create_stock_features(state, item_features) + return stock_features' +end diff --git a/src/DynamicReplenishment/maximizer.jl b/src/DynamicReplenishment/maximizer.jl new file mode 100644 index 00000000..7aec0b3e --- /dev/null +++ b/src/DynamicReplenishment/maximizer.jl @@ -0,0 +1,106 @@ +function _obj_function(N::Int, ub::Vector{Int}, Θ, y, z) + θ = Θ[1:N] + η = Vector{Vector{Float64}}(undef, N) + offset = N + for i in 1:N + η[i] = Θ[(offset + 1):(offset + ub[i])] + offset += ub[i] + end + utility_reward = sum(θ[i] * y[i] for i in 1:N) + stock_penalization = sum( + η[i][1] * z[i, 1] - sum(z[i, j] * sum(η[i][k] for k in 2:j) for j in 2:ub[i]) for + i in 1:N + ) + return utility_reward + stock_penalization +end + +""" +$TYPEDSIGNATURES + +Solve the Replenishment Problem defined by the config and cost vectors θ and η. +""" +function replenishment_problem( + Θ; state::DRPState, y_true=nothing, model_builder=highs_model +) + config = state.config + N = item_count(config) + ub = ub_per_item(state) + t = current_epoch(state) + + if !all(isfinite, Θ) || maximum(abs, Θ) > 1e12 + @warn "Maximiser: Θ infinite" max_abs_theta = maximum(abs, Θ) nonfinite_theta = count( + !isfinite, Θ + ) maxlog = 10 + end + + m = model_builder() + set_silent(m) + # Variables + @variable(m, 0 <= y[i in 1:N] <= ub[i], Int) + # penalization + @variable(m, z[i in 1:N, j in 1:ub[i]], Bin) + + # Objective function + @objective(m, Max, _obj_function(N, ub, Θ, y, z)) + # Constraints + ## penalization constraints + @constraint(m, [i in 1:N], y[i] + state.stock[i] == sum(z[i, j] for j in 1:ub[i])) + ## quota constraints + @constraint( + m, + [c in 1:nb_constraints(config)], + sum(config.constraints_matrix[c, i] * y[i] for i in 1:N) <= config.quotas[t, c] + ) + ## structural constraints + @constraint(m, [i in 1:N, j in 1:(ub[i] - 1)], z[i, j] >= z[i, j + 1]) + + if !isnothing(y_true) + z_true = get_z_from_y(y_true, state) + for i in 1:N + fix(y[i], y_true[i]; force=true) + for j in 1:ub[i] + fix(z[i, j], z_true[i, j]; force=true) + end + end + end + + optimize!(m) + + if primal_status(m) == MOI.FEASIBLE_POINT + obj = objective_value(m) + if !isfinite(obj) || abs(obj) >= 1e20 + @warn "Maximiser: infinite objective" objective = obj max_abs_theta = maximum( + abs, Θ + ) nonfinite_theta = count(!isfinite, Θ) maxlog = 10 + end + end + + return Int.(round.(value.(y))) +end + +function g(y; state::DRPState, kwargs...) + N = item_count(state.config) + ub = ub_per_item(state) + stock_and_replenishment = round.(Int, state.stock .+ y) + yη = Vector{Float64}(undef, sum(ub)) + row = 1 + for i in 1:N + yη[row] = stock_and_replenishment[i] > 0 ? 1 : 0 + for k in 2:ub[i] + yη[row + k - 1] = -max(0, stock_and_replenishment[i] - (k - 1)) + end + row += ub[i] + end + return vcat(vec(y), vec(yη)) +end + +function get_z_from_y(y_true::Vector{Int}, state::DRPState) + N = length(y_true) + ub = ub_per_item(state) + stock_and_replenishment = round.(Int, state.stock .+ y_true) + z_true = zeros(Int, N, maximum(ub)) + for i in 1:N + z_true[i, 1:min(ub[i], stock_and_replenishment[i])] .= 1 + end + return z_true +end diff --git a/src/DynamicReplenishment/policies.jl b/src/DynamicReplenishment/policies.jl new file mode 100644 index 00000000..9bfb75ac --- /dev/null +++ b/src/DynamicReplenishment/policies.jl @@ -0,0 +1,422 @@ +# policy utils + +function mean_feature_order(env::Environment; rng::AbstractRNG=Xoshiro(nothing)) + item_features = features(env) # (d, N), price is a separate field + # if no features other than price, return a random order + if size(item_features, 1) == 0 + return randperm(rng, item_count(env)) + end + # else order items by the mean of their features + return sortperm(vec(mean(item_features; dims=1))) +end + +function max_quotas_item_after_repl( + item_idx::Int, + replenishment::Vector{Int}, + time_idx::Int, + quotas::Matrix{Int}, + cons_mat::Matrix{Int}, +) + N = size(cons_mat, 2) + nb_constraints = size(cons_mat, 1) + return max( + 0, + minimum([ + quotas[time_idx, c] - sum(replenishment[j] * cons_mat[c, j] for j in 1:N) for + c in 1:nb_constraints if cons_mat[c, item_idx] == 1 + ]), + ) +end + +# Policies + +function greedy_policy(env::Environment; model_builder=highs_model) + _, state = observe(env) + N = item_count(env) + ub = ub_per_item(env) + Θ = zeros(N + sum(ub)) + Θ[1:N] .= prices(env) + return (replenishment_problem(Θ; state, model_builder=model_builder)) +end + +function lazy_policy(env::Environment) + N = item_count(env) + return zeros(Int, N) +end + +function random_policy(env::Environment, rng::AbstractRNG=Xoshiro(nothing)) + N = item_count(env) + cons_mat = constraints_matrix(env) + q = quotas(env) + replenishment = zeros(Int, N) + order_item = randperm(rng, N) + t = current_epoch(env) + for item in order_item + max_quota_item = max_quotas_item_after_repl(item, replenishment, t, q, cons_mat) + replenishment[item] = rand(rng, 0:max_quota_item) + end + return replenishment +end + +""" +$TYPEDSIGNATURES + +Mean replenishment over `anticipative_results`. +When `per_epoch` is `true`, only the samples decided at epoch `t` are averaged. +""" +function mean_replenishment( + anticipative_results::AbstractVector{<:DataSample}, N::Int, t::Int; per_epoch::Bool +) + total = zeros(Float64, N) + count = 0 + for sample in anticipative_results + per_epoch && + hasproperty(sample, :state) && + current_epoch(sample.state) != t && + continue + total .+= sample.y + count += 1 + end + count == 0 && return nothing + return total ./ count +end + +""" +$TYPEDSIGNATURES + +Replenish the mean quantities of a set of anticipative decisions, item by item in `order_item` order. +- `per_epoch=false` averages decisions of the complete dataset +- `per_epoch=true` averages only the decisions taken at the current +epoch +""" +function mean_anticipative_policy( + env::Environment; + rng::AbstractRNG=Xoshiro(nothing), + anticipative_results::AbstractVector{<:DataSample}=DataSample[], + order_item::Function=mean_feature_order, + per_epoch::Bool=false, +) + if isempty(anticipative_results) + @warn "mean_anticipative_policy: no anticipative results provided, falling back to lazy policy" + return lazy_policy(env) + end + N = item_count(env) + t = current_epoch(env) + anticipative_repl = mean_replenishment(anticipative_results, N, t; per_epoch) + if isnothing(anticipative_repl) + @warn "mean_anticipative_policy: no demonstration at epoch $t, averaging over the whole horizon instead" maxlog = + 1 + anticipative_repl = mean_replenishment(anticipative_results, N, t; per_epoch=false) + end + replenishment = zeros(Int, N) + q = quotas(env) + cons_mat = constraints_matrix(env) + for item in order_item(env; rng=rng) + max_quota_item = max_quotas_item_after_repl(item, replenishment, t, q, cons_mat) + replenishment[item] = min(round(Int, anticipative_repl[item]), max_quota_item) + end + return replenishment +end + +""" +$TYPEDEF + +Callable wrapping [`mean_anticipative_policy`](@ref) with a dataset of anticipative results. +""" +struct MeanAnticipativePolicyCall + anticipative_results::Vector{DataSample} + order_item::Function + "average only the decisions decided at the epoch of the decision" + per_epoch::Bool +end + +function MeanAnticipativePolicyCall(anticipative_results, order_item; per_epoch::Bool=false) + return MeanAnticipativePolicyCall(anticipative_results, order_item, per_epoch) +end + +function (p::MeanAnticipativePolicyCall)(env::Environment; kwargs...) + return mean_anticipative_policy( + env; + kwargs..., + anticipative_results=p.anticipative_results, + order_item=p.order_item, + per_epoch=p.per_epoch, + ) +end + +""" +$TYPEDSIGNATURES + +Solve a sample average approximation of the replenishment problem over the remaining horizon +on `nb_scenarios` sampled scenarios, and return the first-stage replenishment decision. + +The first replenishment is constrained to be identical across scenarios, so the returned +decision is implementable without knowing the realized demand. When `θ` is given, the +objective is augmented with `κ * dot(θ, g(y))` to bias the decision towards the predicted +utilities. + +The solver is stopped after `time_limit` seconds (10 minutes by default), returning the +best feasible solution found so far; pass `time_limit=nothing` to disable it. If no feasible +point was found at all, the policy falls back to replenishing nothing. + +`warm_start` hands the solver an initial all-zero replenishment. +""" +function saa_policy( + env::Environment; + nb_scenarios::Int=5, + rng::AbstractRNG=Xoshiro(nothing), + model_builder=highs_model, + verbose::Bool=false, + mip_gap::Float64=0.0, + time_limit::Union{Real,Nothing}=600.0, + warm_start::Bool=true, + θ=nothing, + state::DRPState=env.state, + κ::Float64=1.0, +) + scenarios = [generate_scenario(env.config; rng=rng) for _ in 1:nb_scenarios] + bigM_s = [compute_bigM_sales(env, scenario) for scenario in scenarios] + bigM_ps = [compute_bigM_physical_stock(env, scenario) for scenario in scenarios] + + @assert !is_terminated(env) + + m = model_builder() + verbose || set_silent(m) + set_attribute(m, MOI.RelativeGapTolerance(), mip_gap) + isnothing(time_limit) || set_attribute(m, MOI.TimeLimitSec(), Float64(time_limit)) + N = item_count(env) + T = max_steps(env) - current_epoch(env) + 1 + n_customers = [nb_customers(scenario)[current_epoch(env):end] for scenario in scenarios] + q = quotas(env)[current_epoch(env):end, :] + s0 = stock(env) + ## Variables + @variable(m, y[1:nb_scenarios, 1:T, 1:N] >= 0, Int) # replenishments + @variable(m, s[1:nb_scenarios, 1:(T + 1), 1:N] >= 0, Int) # stock + @variable( + m, + α[s_idx in 1:nb_scenarios, i in 1:(N + 1), t in 1:T, k in 1:n_customers[s_idx][t]], + Bin + ) # sales + @variable(m, v[1:nb_scenarios, 1:(T + 1), 1:N] >= 0, Int) # physical stock + @variable(m, z[1:nb_scenarios, 1:N, 2:(T + 1)], Bin) # auxiliary binary for physical stock linearization + @variable(m, s_min[1:nb_scenarios, 1:T] >= 0, Int) # stock under min + @variable(m, s_sup[1:nb_scenarios, 1:T] >= 0, Int) # stock over max + + ## Constraints + @constraint(m, [s_idx in 2:nb_scenarios, i in 1:N], y[1, 1, i] == y[s_idx, 1, i]) # first replenishment is the same for all scenarios + objective = zero(AffExpr) + for s_idx in 1:nb_scenarios + stock_constraints!( + m, + y[s_idx, :, :], + s[s_idx, :, :], + α[s_idx, :, :, :], + T, + N, + n_customers[s_idx], + s0, + ) + customer_constraints!(m, α[s_idx, :, :, :], T, N, n_customers[s_idx]) + sales_order_constraints!( + m, + y[s_idx, :, :], + s[s_idx, :, :], + α[s_idx, :, :, :], + T, + N, + n_customers[s_idx], + scenarios[s_idx].utilities[current_epoch(env):end], + bigM_s[s_idx], + ) + quota_constraints!(m, y[s_idx, :, :], T, N, constraints_matrix(env), q) + physical_stock_constraints!( + m, + y[s_idx, :, :], + α[s_idx, :, :, :], + v[s_idx, :, :], + z[s_idx, :, :], + T, + N, + delivery_delay(env), + s0, + n_customers[s_idx], + bigM_ps[s_idx], + ) + stock_bounds_constraints!( + m, + v[s_idx, :, :], + T, + N, + s_min[s_idx, :], + s_sup[s_idx, :], + stock_inf(env), + stock_sup(env), + ) + ## Objective + objective += compute_objective( + y[s_idx, :, :], + s[s_idx, :, :], + α[s_idx, :, :, :], + v[s_idx, :, :], + T, + s_min[s_idx, :], + s_sup[s_idx, :], + env, + n_customers[s_idx], + ) + end + + if θ !== nothing + g_y = g_model(m, N, ub_per_item(state), y[1, 1, :], s[1, 1, :]) + @assert length(θ) == N + sum(ub_per_item(state)) + objective += κ * dot(θ, g_y) + end + @objective(m, Max, objective) + + # Warm start: replenishing nothing is always feasible + if warm_start + for s_idx in 1:nb_scenarios, t in 1:T, i in 1:N + set_start_value(y[s_idx, t, i], 0) + end + end + + optimize!(m) + if primal_status(m) == MOI.FEASIBLE_POINT + return round.(Int, value.(y[1, 1, :])) + else + # Fall back to a lazy decision instead of `nothing` + @warn "SAA: no feasible point found, falling back to lazy replenishment" + return zeros(Int, N) + end +end + +""" +$TYPEDEF + +Callable wrapping [`saa_policy`](@ref) with a fixed `model_builder`. +""" +struct SAAPolicyCall{M} + model_builder::M + verbose::Bool + mip_gap::Float64 + nb_scenarios::Int + "solver time limit in seconds, `nothing` to disable" + time_limit::Union{Float64,Nothing} + "hand the solver an all-zero replenishment as initial solution" + warm_start::Bool +end + +function SAAPolicyCall( + model_builder::M; + verbose::Bool=false, + mip_gap::Float64=1e-2, + nb_scenarios::Int=1, + time_limit::Union{Real,Nothing}=600.0, + warm_start::Bool=true, +) where {M} + return SAAPolicyCall{M}( + model_builder, + verbose, + mip_gap, + nb_scenarios, + isnothing(time_limit) ? nothing : Float64(time_limit), + warm_start, + ) +end + +function (p::SAAPolicyCall)(env::Environment; kwargs...) + return saa_policy( + env; + kwargs..., + model_builder=p.model_builder, + verbose=p.verbose, + mip_gap=p.mip_gap, + nb_scenarios=p.nb_scenarios, + time_limit=p.time_limit, + warm_start=p.warm_start, + ) +end + +""" +$TYPEDEF + +Full-horizon anticipative policy for the dynamic replenishment benchmark. + +# Fields +$TYPEDFIELDS +""" +struct AnticipativePolicy{M} <: + Utils.AbstractTrajectoryPolicy{DynamicReplenishmentBenchmark} + "JuMP model builder handed to [`anticipative_solver`](@ref)" + model_builder::M + "relative MIP gap tolerance" + mip_gap::Float64 + "solver time limit in seconds, `nothing` to solve to optimality" + time_limit::Union{Float64,Nothing} +end + +function AnticipativePolicy( + model_builder::M=highs_model; mip_gap::Real=0.0, time_limit::Union{Real,Nothing}=nothing +) where {M} + return AnticipativePolicy{M}( + model_builder, + Float64(mip_gap), + isnothing(time_limit) ? nothing : Float64(time_limit), + ) +end + +""" +$TYPEDSIGNATURES + +Solve the whole remaining horizon at once via [`anticipative_solver`](@ref) and return +`(total_reward, dataset)`. +""" +function (p::AnticipativePolicy)(env::Environment, rng::AbstractRNG; kwargs...) + total_reward, dataset = anticipative_solver( + env, + rng; + reset_env=false, + mip_gap=p.mip_gap, + time_limit=p.time_limit, + kwargs..., + model_builder=p.model_builder, + ) + if isnothing(total_reward) + # `anticipative_solver` already warned; keep the `(reward, dataset)` contract so a + # single infeasible episode does not blow up a whole benchmark run. + return 0.0, DataSample[] + end + return total_reward, dataset +end + +""" +$TYPEDSIGNATURES + +Record a dynamic replenishment rollout step with the fields: +- the pre-decision `state` in the context, +- `next_sales` and `customers` in `extra`: the sales and the number of customers of time t +- the reward of the step in `extra`. +""" +function Utils.rollout_step!( + policy::Utils.AbstractPolicy{DynamicReplenishmentBenchmark}, + env::Environment, + rng::AbstractRNG; + kwargs..., +) + y = policy(env; kwargs...) + features, state = observe(env) + state_copy = deepcopy(state) + t = current_epoch(env) + reward = step!(env, y, rng) + return reward, + DataSample(; + x=features, + y=y, + state=state_copy, + extra=(; + next_sales=sales_history(state)[t, :], + customers=customer_history(state)[t], + reward, + ), + ) +end diff --git a/src/DynamicReplenishment/scenario.jl b/src/DynamicReplenishment/scenario.jl new file mode 100644 index 00000000..5fe188c5 --- /dev/null +++ b/src/DynamicReplenishment/scenario.jl @@ -0,0 +1,43 @@ +""" +$TYPEDEF + +# Fields +$TYPEDFIELDS +""" +@kwdef struct Scenario + "Perturbed utilities for each customers: utilities[t][k][i] = static_utilities[i] + ε[t][k] where ε[t][k] ~ Gumbel(0, 1): utility of archetype i for customer k at time t" + utilities::Vector{Vector{Vector{Float64}}} +end + +function nb_customers(scenario::Scenario) + return [length(scenario.utilities[t]) for t in 1:length(scenario.utilities)] +end + +function Base.getindex(scenario::Scenario, idx::Integer) + return (; utilities=scenario.utilities[idx]) +end + +""" +$TYPEDSIGNATURES + +Sample a scenario given the customer choice model and static utilities. +""" +function Utils.generate_scenario( + config::DynamicReplenishmentBenchmark; + seed=nothing, + rng::AbstractRNG=Xoshiro(seed), + temp=1.0, + random_utility_model=Gumbel(0.0, 1.0), +) + N = item_count(config) + T = max_steps(config) + λ = poisson_arrival_rate(config) + nb_customers = rand(rng, Poisson(λ), T) + utilities = [ + [ + config.static_utilities .+ temp * rand(rng, random_utility_model, N+1) for + _ in 1:nb_customers[t] + ] for t in 1:T + ] + return Scenario(; utilities=utilities) +end diff --git a/src/DynamicReplenishment/state.jl b/src/DynamicReplenishment/state.jl new file mode 100644 index 00000000..ec4b33a3 --- /dev/null +++ b/src/DynamicReplenishment/state.jl @@ -0,0 +1,301 @@ +""" +$TYPEDEF + +State data structure for the Dynamic Replenishment Problem. +Convention: all history matrices are (time, item), i.e. `history[t, i]`. + +# Fields +$TYPEDFIELDS +""" +mutable struct DRPState{B<:DynamicReplenishmentBenchmark} + "The benchmark configuration." + config::B + "Current time step (epoch) of the simulation." + current_epoch::Int + "Current stock levels for each item (N)" + stock::Vector{Int} + "Current physical stock levels for each item (N)" + physical_stock::Vector{Int} + "History of stock levels for each item (current_epoch, N)" + stock_history::Matrix{Int} + "History of replenishments for each item (current_epoch-1, N)" + replenishment_history::Matrix{Int} + "History of sales for each item (current_epoch-1, N)" + sales_history::Matrix{Int} + "History of number of customers per time step (current_epoch-1)" + customer_history::Vector{Int} + "Upper bound of replenishment per item (N)" + ub_per_item::Vector{Int} +end + +""" +$TYPEDSIGNATURES + +Compute physical stock at time `t` from raw historical data. +""" +function compute_physical_stock( + config, + t::Int, + s0::Vector{Int}, + replenishment_history::AbstractMatrix{Int}, + sales_history::AbstractMatrix{Int}, +) + N = item_count(config) + t_repl = t - delivery_delay(config) + t_sales = t - 1 + sales_sum = if t_sales <= 0 + zeros(Int, N) + else + vec(sum(view(sales_history, 1:t_sales, :); dims=1)) + end + if t <= delivery_delay(config) + return max.(0, s0 .- sales_sum) + else + repl_sum = if t_repl <= 0 + zeros(Int, N) + else + vec(sum(view(replenishment_history, 1:t_repl, :); dims=1)) + end + return max.(0, s0 .+ repl_sum .- sales_sum) + end +end + +""" +$TYPEDSIGNATURES + +Compute the cumulative cost of a complete history (epochs 1:T where T = size(sales_history,1)), +from raw data. +""" +function compute_total_cost( + config, + s0::Vector{Int}, + stock_history::AbstractMatrix{Int}, + replenishment_history::AbstractMatrix{Int}, + sales_history::AbstractMatrix{Int}, +) + T = size(sales_history, 1) + T == 0 && return 0.0 + + margin_sales = sum(sum(prices(config) .* sales_history[t, :]) for t in 1:T) + v_stock_cost = sum( + sum(virtual_stock_cost(config) .* stock_history[t + 1, :] for t in 1:T) + ) + phys_stocks = [ + compute_physical_stock(config, t + 1, s0, replenishment_history, sales_history) for + t in 1:T + ] + p_stock_cost = sum(sum(physical_stock_cost(config) .* phys_stocks[t]) for t in 1:T) + stock_bound_cost = sum( + over_stock_bound_cost(config) * ( + max(0, stock_inf(config) - sum(phys_stocks[t])) + + max(0, sum(phys_stocks[t]) - stock_sup(config)) + ) for t in 1:T + ) + return margin_sales - v_stock_cost - p_stock_cost - stock_bound_cost +end + +""" +$TYPEDSIGNATURES + +Construct a `DRPState`. +By convention, the initial physical stock is equal to the initial stock. +""" + +function DRPState(; + config::B, + current_epoch::Int, + stock::Vector{Int}, + stock_history::Matrix{Int}, + replenishment_history::Matrix{Int}, + sales_history::Matrix{Int}, + customer_history::Vector{Int}, + ub_per_item::Vector{Int}, +) where {B<:DynamicReplenishmentBenchmark} + s0 = stock_history[1, :] + physical_stock = compute_physical_stock( + config, current_epoch, s0, replenishment_history, sales_history + ) + return DRPState{B}( + config, + current_epoch, + stock, + physical_stock, + stock_history, + replenishment_history, + sales_history, + customer_history, + ub_per_item, + ) +end + +function DRPState( + config::B, stock_ini::Vector{Int} +) where {B<:DynamicReplenishmentBenchmark} + N = length(stock_ini) + return DRPState(; + config, + current_epoch=1, + stock=copy(stock_ini), + stock_history=reshape(copy(stock_ini), 1, N), + replenishment_history=zeros(Int, 0, N), + sales_history=zeros(Int, 0, N), + customer_history=Int[], + ub_per_item=stock_ini .+ max_quotas(config)[1, :], + ) +end + +""" +$TYPEDSIGNATURES + +Deep-copy a `DRPState` without copying `config`: the benchmark configuration is never +mutated after construction and is shared across every state of an episode, so cloning it +on every `deepcopy` (e.g. once per epoch in `rollout!`) would be wasted work (it can hold +heavy fields such as a `customer_choice_model` neural network). +""" +function Base.deepcopy_internal(state::DRPState{B}, stackdict::IdDict) where {B} + haskey(stackdict, state) && return stackdict[state] + new_state = DRPState{B}( + state.config, + state.current_epoch, + Base.deepcopy_internal(state.stock, stackdict), + Base.deepcopy_internal(state.physical_stock, stackdict), + Base.deepcopy_internal(state.stock_history, stackdict), + Base.deepcopy_internal(state.replenishment_history, stackdict), + Base.deepcopy_internal(state.sales_history, stackdict), + Base.deepcopy_internal(state.customer_history, stackdict), + Base.deepcopy_internal(state.ub_per_item, stackdict), + ) + stackdict[state] = new_state + return new_state +end + +current_epoch(state::DRPState) = state.current_epoch +stock(state::DRPState) = state.stock +physical_stock(state::DRPState) = state.physical_stock +total_stock(state::DRPState) = sum(state.stock) +stock_history(state::DRPState) = state.stock_history +replenishment_history(state::DRPState) = state.replenishment_history +sales_history(state::DRPState) = state.sales_history +customer_history(state::DRPState) = state.customer_history +stock_ini(state::DRPState) = stock_history(state)[1, :] +ub_per_item(state::DRPState) = state.ub_per_item + +""" +$TYPEDSIGNATURES + +Compute the cumulative cost of the state's history so far, from raw data (see +[`compute_total_cost`](@ref)). +""" +function total_cost(state::DRPState) + return compute_total_cost( + state.config, + stock_ini(state), + stock_history(state), + replenishment_history(state), + sales_history(state), + ) +end + +function reset_state!( + state::DRPState, rng::AbstractRNG; reset_stock_ini=false, stock_ini_fill_rate::Real=0.5 +) + N = item_count(state.config) + if reset_stock_ini + s0 = draw_stock_ini(rng, N, stock_sup(state.config), stock_ini_fill_rate) + else + s0 = stock_ini(state) + end + state.current_epoch = 1 + state.stock = copy(s0) + state.physical_stock = copy(s0) + state.stock_history = reshape(copy(s0), 1, N) + state.replenishment_history = zeros(Int, 0, N) + state.sales_history = zeros(Int, 0, N) + state.customer_history = Int[] + state.ub_per_item = s0 .+ max_quotas(state.config)[1, :] + return state +end + +function is_feasible(state::DRPState, replenishment::Vector{Int}; verbose=false) + config = state.config + cons_mat = constraints_matrix(config) + q = quotas(config) + for c in 1:nb_constraints(config) + if sum(cons_mat[c, :] .* replenishment) > q[state.current_epoch, c] + verbose && + @warn "Replenishment violates quota constraint $c at epoch $(state.current_epoch) : $(sum(cons_mat[c, :] .* replenishment)) > $(q[state.current_epoch, c])" + return false + end + end + return true +end + +""" +$TYPEDSIGNATURES + +Compute the cost delta incurred at the current epoch (sales margin minus virtual and +physical stock costs minus over/under stock bound penalties), matching +[`compute_total_cost`](@ref) exactly. Must be called after `state.sales_history` and +`state.stock_history` have been updated for the current epoch (i.e. once this epoch's +sales are known). Does not read or write `state.physical_stock`: the physical stock used +here already reflects this epoch's own (now realized) sales, whereas `state.physical_stock` +must keep representing the physical stock as observed *before* those sales, so that the +state used to make the next decision never leaks this epoch's outcome. +""" +function update_cost!(state::DRPState) + config = state.config + t = current_epoch(state) + # sales reward + sales_t = view(state.sales_history, t, :) + margin = sum(prices(config) .* sales_t) + # physical stock after this epoch's replenishment delivery and sales are accounted for + post_sale_physical_stock = compute_physical_stock( + config, t + 1, stock_ini(state), state.replenishment_history, state.sales_history + ) + physical_cost = sum(physical_stock_cost(config) .* post_sale_physical_stock) + # virtual stock cost + virtual_stock = view(state.stock_history, t + 1, :) + virtual_cost = sum(virtual_stock_cost(config) .* virtual_stock) + # over / under stock costs (based on physical stock, matching the anticipative solver) + total_physical = sum(post_sale_physical_stock) + under = max(0, stock_inf(config) - total_physical) + over = max(0, total_physical - stock_sup(config)) + penalty = over_stock_bound_cost(config) * (under + over) + + return margin - virtual_cost - physical_cost - penalty +end + +function apply_replenishment!(state::DRPState, replenishment::Vector{Int}) + state.stock .+= replenishment + state.replenishment_history = vcat(state.replenishment_history, replenishment') + return nothing +end + +function apply_sales!(state::DRPState; utilities::Vector{Vector{Float64}}) + N = length(state.stock) + sales = zeros(Int, N) + nb_customers = length(utilities) + for k in 1:nb_customers + order_of_sales = sortperm(utilities[k]; rev=true) + for item_index in order_of_sales + if item_index == N + 1 + break + end + if state.stock[item_index] > 0 + sales[item_index] += 1 + state.stock[item_index] -= 1 + break + end + end + end + state.sales_history = vcat(state.sales_history, sales') + state.stock_history = vcat(state.stock_history, state.stock') + delta_cost = update_cost!(state) + return delta_cost +end + +function add_customers!(state::DRPState; utilities::Vector{Vector{Float64}}) + nb_customers = length(utilities) + state.customer_history = push!(state.customer_history, nb_customers) + return state +end diff --git a/src/DynamicReplenishment/statistical_model.jl b/src/DynamicReplenishment/statistical_model.jl new file mode 100644 index 00000000..d28e20af --- /dev/null +++ b/src/DynamicReplenishment/statistical_model.jl @@ -0,0 +1,40 @@ +""" +$TYPEDEF + +# Fields +$TYPEDFIELDS +""" +@kwdef struct StatisticalModel{L1,L2} + "replenishment reward" + θ_model::L1 + "stock penalization" + η_model::L2 +end + +@layer StatisticalModel + +""" +$TYPEDSIGNATURES + +""" +function Utils.generate_statistical_model( + b::DynamicReplenishmentBenchmark; seed=nothing, kwargs... +) + isnothing(seed) || seed!(seed) + θ_model = Chain(Dense(item_features_size(b) => 1)) + η_model = Chain(Dense(stock_features_size(b) => 1), softplus) + return StatisticalModel(; θ_model, η_model) +end + +function (m::StatisticalModel)(x) + item_ids = @view x[end, :] + starts = [findfirst(==(i), item_ids) for i in 1:maximum(Int, item_ids)] + + # the stock block and the trailing item identifier sit on top of the item block + nb_item_features = size(x, 1) - (NB_STOCK_FEATURES + 1) + x_features = @view x[1:(end - 1), :] + x_item = x_features[1:(nb_item_features), starts] + θ = m.θ_model(x_item) + η = m.η_model(x_features) + return vcat(vec(θ), vec(η)) +end diff --git a/src/DynamicReplenishment/utils.jl b/src/DynamicReplenishment/utils.jl new file mode 100644 index 00000000..8f71ed76 --- /dev/null +++ b/src/DynamicReplenishment/utils.jl @@ -0,0 +1,3 @@ +mean_or_zero(x) = isempty(x) ? 0.0 : mean(x) +max_or_zero(x) = isempty(x) ? 0.0 : maximum(x) +min_or_zero(x) = isempty(x) ? 0.0 : minimum(x) diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index b0b0f3f8..35810d7a 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -3,7 +3,7 @@ module Utils using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES using Flux: softplus using HiGHS: HiGHS -using JuMP: Model +using JuMP: Model, set_attribute using LinearAlgebra: dot using Random: Random, Xoshiro, AbstractRNG using SCIP: SCIP diff --git a/src/Utils/interface/dynamic_benchmark.jl b/src/Utils/interface/dynamic_benchmark.jl index 76f596fc..cb13d815 100644 --- a/src/Utils/interface/dynamic_benchmark.jl +++ b/src/Utils/interface/dynamic_benchmark.jl @@ -17,6 +17,9 @@ of decisions) as in [`AbstractStochasticBenchmark`](@ref). environment must not manage its own seed/rng: draw randomness from the passed `rng`. Users obtain wrapped environments via [`generate_environment`](@ref) (one) or [`generate_environments`](@ref) (many). +- [`build_environment`](@ref)`(bench, sample, scenario)`: rebuild a bare environment sitting + at the state carried by `sample` and running on `scenario`. Lets solvers that take an + environment be evaluated at any state of a stored trajectory. - [`generate_baseline_policies`](@ref)`(bench)`: returns named baseline callables of signature `(env) -> Vector{DataSample}` (full trajectory rollout). - [`generate_anticipative_solver`](@ref)`(bench)`: returns a callable @@ -70,6 +73,21 @@ environments cannot be drawn independently (e.g. loaded from files), override """ function build_environment end +""" + build_environment(::AbstractDynamicBenchmark, sample::DataSample, scenario) + -> AbstractEnvironment + +**Optional.** Rebuild a bare environment positioned at the state carried by `sample` and +running on `scenario`. Useful for solvers that takes an environment as argument. +Useful for algorithms that need to run an anticipative or parametric solver from a state of +a stored trajectory. +""" +function build_environment(b::AbstractDynamicBenchmark, ::DataSample, scenario) + return error( + "build_environment(::$(typeof(b)), ::DataSample, scenario) is not implemented" + ) +end + """ $TYPEDSIGNATURES diff --git a/src/Utils/model_builders.jl b/src/Utils/model_builders.jl index 4f0c838b..ed340a40 100644 --- a/src/Utils/model_builders.jl +++ b/src/Utils/model_builders.jl @@ -15,5 +15,7 @@ Initialize a SCIP model (with disabled logging). """ function scip_model() model = Model(SCIP.Optimizer) + # Accept partial primal starts however few variables they fix (default 0.85) + set_attribute(model, "heuristics/completesol/maxunknownrate", 1.0) return model end diff --git a/test/replenishment.jl b/test/replenishment.jl new file mode 100644 index 00000000..c01cc939 --- /dev/null +++ b/test/replenishment.jl @@ -0,0 +1,627 @@ +using Statistics: mean + +const DR = DecisionFocusedLearningBenchmarks.DynamicReplenishment + +@testset "DynamicReplenishment - Benchmark Construction" begin + b = DynamicReplenishmentBenchmark() + @test b.N == 10 + @test b.λ == 15 + @test b.d == 5 + @test b.stock_inf == 0 + @test b.stock_sup == 30 + @test b.ub_same_item == 30 + @test b.delivery_delay == 3 + @test b.max_steps == 10 + @test size(b.constraints_matrix) == (12, 10) + @test size(b.quotas) == (10, 12) + @test size(b.max_quotas) == (b.max_steps, b.N) + @test all(b.max_quotas .≥ 0) + @test all(b.max_quotas .≤ b.ub_same_item) + + b_custom = DynamicReplenishmentBenchmark(; + N=5, + λ=10, + constraints_matrix=[1 1 1 0 0; 0 0 0 1 1; 0 0 0 0 1], + quotas=[20 15 5; 10 20 5], + d=3, + stock_inf=2, + stock_sup=10, + ub_same_item=17, + delivery_delay=1, + max_steps=2, + ) + @test b_custom.N == 5 + @test b_custom.λ == 10 + @test b_custom.d == 3 + @test b_custom.stock_inf == 2 + @test b_custom.stock_sup == 10 + @test b_custom.ub_same_item == 17 + @test b_custom.delivery_delay == 1 + @test b_custom.max_steps == 2 + @test size(b_custom.constraints_matrix) == (8, 5) + @test b_custom.quotas[1, :] == [20, 15, 5, 17, 17, 17, 17, 17] + @test size(b_custom.quotas) == (2, 8) + + @test b_custom.max_quotas[1, :] == [17, 17, 17, 15, 5] + @test b_custom.max_quotas[2, :] == [10, 10, 10, 17, 5] + + @test DR.item_count(b) == 10 + @test DR.feature_count(b) == 5 + @test DR.max_steps(b) == 10 + @test DR.stock_inf(b) == 0 + @test DR.stock_sup(b) == 30 + @test DR.ub_same_item(b) == 30 + @test DR.delivery_delay(b) == 3 + @test DR.poisson_arrival_rate(b) == 15.0 + @test length(DR.prices(b)) == 10 + @test all(1.0 .≤ DR.prices(b) .≤ 10.0) + @test size(DR.features(b)) == (5, 10) + @test length(DR.virtual_stock_cost(b)) == 10 + @test length(DR.physical_stock_cost(b)) == 10 + @test DR.nb_constraints(b) == 12 + + @test size(DR.scaled_features(b)) == (6, 10) + @test all(isapprox.(vec(sum(DR.scaled_features(b); dims=2)), 0.0; atol=1e-10)) + + # Two instances sharing prices and customer weights but differing in features: + # same catalogue, different stores. + shared = (; + prices=DR.prices(b_custom), + customer_choice_model=DR.customer_choice_model(b_custom), + constraints_matrix=[1 1 1 0 0; 0 0 0 1 1; 0 0 0 0 1], + quotas=[20 15 5; 10 20 5], + N=5, + d=3, + max_steps=2, + ) + b_same = DynamicReplenishmentBenchmark(; features=DR.features(b_custom), shared...) + @test DR.prices(b_same) == DR.prices(b_custom) + @test b_same.static_utilities == b_custom.static_utilities + + b_other = DynamicReplenishmentBenchmark(; features=zeros(3, 5) .+ (1:5)', shared...) + @test DR.prices(b_other) == DR.prices(b_custom) + @test b_other.static_utilities != b_custom.static_utilities + + @test_throws AssertionError DynamicReplenishmentBenchmark(; N=4, prices=[1.0, 2.0]) + @test_throws AssertionError DynamicReplenishmentBenchmark(; + N=4, d=2, features=zeros(3, 4) + ) +end + +@testset "DynamicReplenishment - Environment Initialization" begin + b = DynamicReplenishmentBenchmark() + rng = Xoshiro(42) + env1 = DR.Environment(b, rng) + @test !is_terminated(env1) + @test DR.item_count(env1) == 10 + @test DR.max_steps(env1) == 10 + @test length(DR.stock_ini(env1)) == 10 + # The total is a fixed fraction of stock_sup, spread over the items: it does not + # scale with N the way a per-item draw would. + @test sum(DR.stock_ini(env1)) == round(Int, 0.5 * DR.stock_sup(b)) + @test all(DR.stock_ini(env1) .≥ 0) + + env_full = DR.Environment(b, rng; stock_ini_fill_rate=1.0) + @test sum(DR.stock_ini(env_full)) == DR.stock_sup(b) + env_empty = DR.Environment(b, rng; stock_ini_fill_rate=0.0) + @test all(iszero, DR.stock_ini(env_empty)) + + # Same fill rate, twice the items: the shelf still starts at the same level. + b_wide = DynamicReplenishmentBenchmark(; N=20) + @test sum(DR.stock_ini(DR.Environment(b_wide, rng))) == + sum(DR.stock_ini(DR.Environment(b, rng))) + + @test_throws ArgumentError DR.Environment(b, rng; stock_ini_fill_rate=1.5) + + @test DR.current_epoch(env1) == 1 + @test env1.stock_ini == DR.stock_ini(env1) + @test DR.stock(env1) == DR.stock_ini(env1) + + state_ini = env1.state + @test state_ini.current_epoch == 1 + @test state_ini.stock == DR.stock_ini(env1) + @test size(state_ini.stock_history) == (1, 10) + @test size(state_ini.replenishment_history) == (0, 10) + @test size(state_ini.sales_history) == (0, 10) + @test length(state_ini.customer_history) == 0 + @test DR.total_cost(state_ini) == 0.0 + + # custom environment + env2 = DR.Environment(b, rng; stock_ini=fill(5, 10)) + @test DR.stock_ini(env2) == fill(5, 10) + @test DR.stock(env2) == fill(5, 10) +end + +@testset "DynamicReplenishment - Environment Reset" begin + b = DynamicReplenishmentBenchmark() + rng = Xoshiro(42) + env = DR.Environment(b, rng) + + s0 = copy(DR.stock_ini(env)) + N = DR.item_count(b) + repl = zeros(Int, N) + step!(env, repl, rng) + reset!(env, rng) + + @test !is_terminated(env) + @test DR.stock(env) == s0 + @test DR.current_epoch(env) == 1 +end + +@testset "DynamicReplenishment - Environment Step" begin + b = DynamicReplenishmentBenchmark() + rng = Xoshiro(42) + env = DR.Environment(b, rng) + N = DR.item_count(b) + + action = zeros(Int, N) + reward = step!(env, action, rng) + @test reward isa Float64 + @test DR.current_epoch(env) == 2 + + # run to termination + while !is_terminated(env) + repl = zeros(Int, N) + @test DR.is_feasible(env.state, repl) + step!(env, repl, rng) + end + @test is_terminated(env) + @test_throws AssertionError step!(env, zeros(Int, N), rng) +end + +@testset "DynamicReplenishment - Feasibility" begin + b = DynamicReplenishmentBenchmark( + N=2, max_steps=2, constraints_matrix=[1 1], quotas=[1; 1] + ) + N = DR.item_count(b) + stock_ini = [0, 0] + rng = Xoshiro(42) + env = DR.Environment(b, rng; stock_ini=stock_ini) + # zero replenishment is always feasible + @test DR.is_feasible(env.state, zeros(Int, N)) + @test DR.is_feasible(env.state, [0, 1]) + @test DR.is_feasible(env.state, [1, 0]) + @test DR.is_feasible(env.state, [0, 0]) + + # replenishment exceeding quota is infeasible + @test !DR.is_feasible(env.state, [2, 0]) + @test !DR.is_feasible(env.state, [0, 2]) + @test !DR.is_feasible(env.state, [1, 1]) +end + +@testset "DynamicReplenishment - State" begin + b = DynamicReplenishmentBenchmark() + rng = Xoshiro(42) + env = DR.Environment(b, rng) + N = DR.item_count(b) + state = env.state + + @test DR.current_epoch(state) == 1 + @test length(DR.stock(state)) == N + @test size(DR.stock_history(state)) == (1, N) + @test size(DR.replenishment_history(state)) == (0, N) + @test size(DR.sales_history(state)) == (0, N) + @test length(DR.customer_history(state)) == 0 + @test DR.total_cost(state) == 0.0 + @test DR.stock_ini(state) == DR.stock_ini(env) + + # after one step + reward = step!(env, zeros(Int, N), rng) + @test DR.current_epoch(state) == 2 + @test size(DR.stock_history(state)) == (2, N) + @test size(DR.replenishment_history(state)) == (1, N) + @test size(DR.sales_history(state)) == (1, N) + @test length(DR.customer_history(state)) == 1 + @test isapprox(DR.total_cost(state), reward; atol=1e-8) +end + +@testset "DynamicReplenishment - Observe" begin + b = DynamicReplenishmentBenchmark() + rng = Xoshiro(42) + env = DR.Environment(b, rng) + N = DR.item_count(b) + ub = DR.ub_per_item(env) + + x, state = observe(env) + + @test !any(isnan.(x)) + + # x is stock_features' : (nb_features, sum(UB)) + @test size(x, 2) == sum(ub) + # the item block, the stock block and the trailing item identifier + @test size(x, 1) == DR.stock_features_size(b) + 1 + static_features = x[1:(DR.feature_count(b) + 1), :] + + # The static block is the *scaled* price/features matrix, the one + # `static_utilities` was built from — not the raw values. + @test DR.create_items_features(state)[:, 1:(DR.feature_count(b) + 1)]' ≈ + Float32.(DR.scaled_features(b)) + + starts = [1; cumsum(ub)[1:(end - 1)] .+ 1] + ends = cumsum(ub) + for i in 1:N + rows = starts[i]:ends[i] + ref = static_features[:, starts[i]] + block = static_features[:, rows] + @test all(block .≈ ref) + end +end + +@testset "DynamicReplenishment - Feature block sizes" begin + # The statistical model slices `x` by hand, so the declared sizes and the actual + # feature matrix must stay in lockstep: a feature added on one side only would + # silently shift every column of the other block. + for b in ( + DynamicReplenishmentBenchmark(; seed=0), + DynamicReplenishmentBenchmark(; N=4, d=2, delivery_delay=1, max_steps=4, seed=1), + DynamicReplenishmentBenchmark(; N=3, d=7, delivery_delay=5, max_steps=4, seed=2), + ) + rng = Xoshiro(0) + env = DR.Environment(b, rng) + state = env.state + + item_features = DR.create_items_features(state) + @test size(item_features) == (DR.item_count(b), DR.item_features_size(b)) + + x, _ = observe(env) + @test size(x, 1) == DR.stock_features_size(b) + 1 + @test size(x, 1) - (DR.NB_STOCK_FEATURES + 1) == DR.item_features_size(b) + + model = generate_statistical_model(b) + @test size(model.θ_model[1].weight, 2) == DR.item_features_size(b) + @test size(model.η_model[1].weight, 2) == DR.stock_features_size(b) + @test all(isfinite.(model(x))) + end +end + +@testset "DynamicReplenishment - Physical stock features" begin + # λ=1 keeps the demand far below the initial stock, so the physical stock never hits + # the `max(0, ...)` floor and the transit identity below is exact + b = DynamicReplenishmentBenchmark(; + N=3, d=2, λ=1, delivery_delay=3, max_steps=6, stock_inf=5, stock_sup=20, seed=0 + ) + rng = Xoshiro(1) + env = DR.Environment(b, rng; stock_ini=[6, 5, 4]) + step!(env, [2, 0, 1], rng) + step!(env, [0, 3, 0], rng) + + state = env.state + @test DR.current_epoch(state) == 3 + item_features = DR.create_items_features(state) + + nb_static = DR.feature_count(b) + 1 + phys = DR.physical_stock(state) + virt = DR.stock(state) + + # physical stock and its price-scaled version + @test item_features[:, nb_static + 10] ≈ Float32.(phys) + @test item_features[:, nb_static + 11] ≈ Float32.(phys .* DR.prices(b)) + # stock in transit: the two orders placed so far, none of them delivered yet since + # delivery_delay = 3 and we are at epoch 3 + @test item_features[:, nb_static + 12] ≈ Float32.(virt .- phys) + @test item_features[:, nb_static + 12] ≈ Float32.([2, 3, 1]) + @test item_features[:, nb_static + 13] ≈ Float32.((virt .- phys) .* DR.prices(b)) + + # state-level features are identical for every item + total_physical = sum(phys) + @test all(item_features[:, nb_static + 14] .≈ Float32(total_physical)) + @test all(item_features[:, nb_static + 15] .≈ Float32(total_physical - DR.stock_inf(b))) + @test all(item_features[:, nb_static + 16] .≈ Float32(DR.stock_sup(b) - total_physical)) + @test all( + item_features[:, nb_static + 17] .≈ + Float32(DR.max_steps(b) - DR.current_epoch(state)), + ) + # the state-level block ends there: the last dynamic column is the horizon + @test nb_static + 17 == DR.item_features_size(b) +end + +@testset "DynamicReplenishment - Stock level features stay linear in j" begin + # The bound penalty is paid on the physical stock, which the replenishment contained in + # a virtual level `j` only reaches `delivery_delay` epochs later: a hinge in `j` would + # place its kink at a threshold the sales of the lead time will have moved. The stock + # level block is therefore linear in `j`, bounds included. + b = DynamicReplenishmentBenchmark(; + N=2, d=2, λ=1, delivery_delay=2, max_steps=4, stock_inf=10, stock_sup=12, seed=4 + ) + state = DR.Environment(b, Xoshiro(0); stock_ini=[2, 2]).state + ub = DR.ub_per_item(state) + stock_features = DR.create_stock_features(state, DR.create_items_features(state)) + nb_fi = DR.item_features_size(b) + rows_1 = 1:ub[1] + # levels straddling both bounds, so a hinge would be visible if one had been added + @test ub[1] >= 12 + @test stock_features[rows_1, nb_fi + 9] ≈ Float32.((2 .+ (1:ub[1])) .- 10) + @test stock_features[rows_1, nb_fi + 11] ≈ Float32.(12 .- (2 .+ (1:ub[1]))) + @test any(stock_features[rows_1, nb_fi + 9] .< 0) + @test any(stock_features[rows_1, nb_fi + 11] .< 0) +end + +@testset "DynamicReplenishment - Coupled stock bound features" begin + b = DynamicReplenishmentBenchmark(; + N=3, d=2, delivery_delay=2, max_steps=5, stock_inf=8, stock_sup=25, seed=3 + ) + rng = Xoshiro(2) + env = DR.Environment(b, rng; stock_ini=[5, 4, 3]) + step!(env, [1, 2, 0], rng) + + state = env.state + ub = DR.ub_per_item(state) + item_features = DR.create_items_features(state) + stock_features = DR.create_stock_features(state, item_features) + nb_fi = DR.item_features_size(b) + + phys = DR.physical_stock(state) + total_physical = sum(phys) + starts = [1; cumsum(ub)[1:(end - 1)] .+ 1] + ends = cumsum(ub) + + for i in 1:DR.item_count(b) + rows = starts[i]:ends[i] + js = 1:ub[i] + p = DR.prices(b)[i] + others = total_physical - phys[i] + expected_inf = (others .+ js) .- DR.stock_inf(b) + expected_sup = DR.stock_sup(b) .- (others .+ js) + + @test stock_features[rows, nb_fi + 9] ≈ Float32.(expected_inf) + @test stock_features[rows, nb_fi + 10] ≈ Float32.(expected_inf .* p) + @test stock_features[rows, nb_fi + 11] ≈ Float32.(expected_sup) + @test stock_features[rows, nb_fi + 12] ≈ Float32.(expected_sup .* p) + # the coupled deviations stay linear in j: no hinge at this level, since the + # replenishment j contains only reaches the physical stock delivery_delay epochs + # later, so a kink here would sit at a threshold the future does not respect + @test all(diff(stock_features[rows, nb_fi + 9]) .≈ 1) + @test all(diff(stock_features[rows, nb_fi + 11]) .≈ -1) + # the item identifier stays the very last column + @test all(stock_features[rows, end] .≈ i) + end + + # `stock_inf`/`stock_sup` are global bounds: comparing them to the level of a single + # item ignores the stock the other items hold, which is exactly `others_physical` + @test stock_features[:, nb_fi + 1] != stock_features[:, nb_fi + 9] + + # the uncoupled and coupled columns agree only for an item whose siblings are empty + env_solo = DR.Environment(b, Xoshiro(2); stock_ini=[5, 0, 0]) + state_solo = env_solo.state + sf_solo = DR.create_stock_features(state_solo, DR.create_items_features(state_solo)) + ub_solo = DR.ub_per_item(state_solo) + starts_solo = [1; cumsum(ub_solo)[1:(end - 1)] .+ 1] + ends_solo = cumsum(ub_solo) + rows_1 = starts_solo[1]:ends_solo[1] + @test sf_solo[rows_1, nb_fi + 1] ≈ sf_solo[rows_1, nb_fi + 9] + @test sf_solo[rows_1, nb_fi + 3] ≈ sf_solo[rows_1, nb_fi + 11] + rows_2 = starts_solo[2]:ends_solo[2] + @test sf_solo[rows_2, nb_fi + 9] ≈ sf_solo[rows_2, nb_fi + 1] .+ 5 + @test sf_solo[rows_2, nb_fi + 11] ≈ sf_solo[rows_2, nb_fi + 3] .- 5 +end + +@testset "DynamicReplenishment - Statistical Model" begin + b = DynamicReplenishmentBenchmark() + N = DR.item_count(b) + + model = generate_statistical_model(b) + @test model isa DR.StatisticalModel + + rng = Xoshiro(42) + env = DR.Environment(b, rng) + x, _ = observe(env) + @test !any(isnan.(x)) + + ub = DR.ub_per_item(env) + θη = model(x) + @test length(θη) == N + sum(ub) + @test all(isfinite.(θη)) +end + +@testset "DynamicReplenishment - Anticipative trajectory is feasible for the CO layer" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=3) + rng = Xoshiro(0) + env = generate_environments(b, 1; seed=0) + ant_solver = generate_anticipative_solver(b) + ant_traj = ant_solver(env[1]) + model = generate_statistical_model(b) + maximizer = generate_maximizer(b) + while !is_terminated(env[1]) + x, _ = observe(env[1]) + θ = model(x) + y_true = ant_traj[DR.current_epoch(env[1].env)].y + y_hat = maximizer(θ; state=env[1].env.state, y_true=y_true) + @test y_true == y_hat + step!(env[1].env, y_true, rng) + end +end + +@testset "DynamicReplenishment - Parametric Anticipative trajectory is feasible for the CO layer" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=3) + rng = Xoshiro(0) + env = generate_environments(b, 1; seed=0) + model = generate_statistical_model(b) + param_ant_solver = DR.generate_parametric_anticipative_solver(b) + while !is_terminated(env[1]) + x, _ = observe(env[1]) + @test !any(isnan.(x)) + θ = model(x) + y_true = param_ant_solver(θ, env[1].env.scenario, env[1])[1].y + y_hat = DR.generate_maximizer(b)(θ; state=env[1].env.state, y_true=y_true) + @test y_true == y_hat + step!(env[1].env, y_true, rng) + end +end +@testset "DynamicReplenishment - Policies" begin + b = DynamicReplenishmentBenchmark() + environments = generate_environments(b, 5; seed=0) + policies = generate_baseline_policies(b) + + @test policies.greedy.name == "Greedy" + @test policies.random.name == "Random" + @test policies.lazy.name == "Lazy" + + r_greedy, greedy_traj = evaluate_policy!(policies.greedy, environments, 5) + @test length(r_greedy) == length(environments) + env = environments[1] + reset!(env) + greedy_action = policies.greedy(env.env) + @test DR.is_feasible(env.env.state, greedy_action) + random_action = policies.random(env.env) + @test DR.is_feasible(env.env.state, random_action) + lazy_action = policies.lazy(env.env) + @test DR.is_feasible(env.env.state, lazy_action) + @test all(lazy_action .== 0) + mean_ant_action = DR.mean_anticipative_policy(env.env; anticipative_results=greedy_traj) + @test DR.is_feasible(env.env.state, mean_ant_action) + @test sort(DR.mean_feature_order(env.env)) == 1:DR.item_count(b) +end + +@testset "DynamicReplenishment - MeanAnticipative per epoch" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=3) + env = generate_environments(b, 1; seed=0)[1] + rng = Xoshiro(0) + _, ant_traj = DR.anticipative_solver(env.env, rng) + N = DR.item_count(b) + + policies = generate_baseline_policies(b; anticipative_results=ant_traj) + @test policies.mean_anticipative.name == "MeanAnticipative" + @test policies.mean_anticipative_per_epoch.name == "MeanAnticipativePerEpoch" + + # Only one expert trajectory here: the per-epoch mean is that epoch's decision, + # whereas the whole-horizon mean mixes every epoch together. + for (t, sample) in enumerate(ant_traj) + @test DR.mean_replenishment(ant_traj, N, t; per_epoch=true) ≈ sample.y + end + horizon_mean = sum(sample.y for sample in ant_traj) ./ length(ant_traj) + @test DR.mean_replenishment(ant_traj, N, 1; per_epoch=false) ≈ horizon_mean + + # At the first epoch the per-epoch variant replays the expert decision exactly: + # it satisfies the quotas, so nothing is clipped. + @test DR.current_epoch(env.env) == 1 + per_epoch_action = policies.mean_anticipative_per_epoch(env.env) + @test per_epoch_action == ant_traj[1].y + @test DR.is_feasible(env.env.state, per_epoch_action) + + # An epoch no demonstration covers: the mean is empty, and the policy falls back + # to the whole-horizon mean instead of failing. + @test isnothing(DR.mean_replenishment(ant_traj, N, DR.max_steps(b) + 1; per_epoch=true)) + + # Without demonstrations, both variants fall back to Lazy. + empty_policies = generate_baseline_policies(b) + @test all(iszero, empty_policies.mean_anticipative_per_epoch(env.env)) +end + +@testset "DynamicReplenishment - Anticipative Solver" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=3) + rng = Xoshiro(42) + env = generate_environments(b, 1; seed=0) + scenario = DR.generate_scenario(b; rng=rng) + env[1].env.scenario = scenario + ant_obj, ant_traj = DR.anticipative_solver(env[1].env, rng, scenario) + policies = generate_baseline_policies(b) + r_greedy, greedy_traj = evaluate_policy!(policies.greedy, env) + @test length(ant_traj) == DR.max_steps(b) == length(greedy_traj) + for g_sample in greedy_traj + @test DR.is_feasible(g_sample.state, g_sample.y) + end + for ant_sample in ant_traj + @test DR.is_feasible(ant_sample.state, ant_sample.y) + end + @test r_greedy[1] <= ant_obj +end + +@testset "DynamicReplenishment - rollout cost matches anticipative solver objective" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=4) + env = generate_environment(b; seed=0) + # `AnticipativePolicy` is a trajectory policy: `evaluate_policy!` resets the wrapper to + # its initial seed, then returns the MILP objective and the planned trajectory in one + # call, without stepping the environment. + ant_obj, ant_traj = evaluate_policy!(DR.AnticipativePolicy(), env) + @test length(ant_traj) == DR.max_steps(b) + + # Replaying those decisions through the simulator must reproduce the MILP objective: + # this is what makes `ant_obj` usable as the reference bound for optimality gaps. + # Resetting to the same seed replays the very same scenario. + reset_to_initial!(env) + total_reward = sum(step!(env.env, sample.y, env.rng) for sample in ant_traj) + + @test isapprox(total_reward, ant_obj; rtol=1e-5) + @test isapprox(DR.total_cost(env.env.state), ant_obj; rtol=1e-5) +end + +@testset "DynamicReplenishment - rollout logs the anticipative solver's fields" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=4) + env = generate_environment(b; seed=0) + + ant_obj, ant_traj = evaluate_policy!(DR.AnticipativePolicy(), env) + total_reward, traj = evaluate_policy!(generate_baseline_policies(b).greedy, env) + + # A rollout sample carries everything an anticipative sample does, plus the step reward. + for field in propertynames(ant_traj[1]) + @test hasproperty(traj[1], field) + end + @test hasproperty(traj[1], :reward) + # `state` (not `instance`) is what both dataset kinds use, so `per_epoch` filtering works + @test all(DR.current_epoch(s.state) == t for (t, s) in enumerate(traj)) + @test DR.mean_replenishment(traj, DR.item_count(b), 2; per_epoch=true) ≈ traj[2].y + + # `next_sales` / `customers` describe the epoch the decision was taken at + final_state = env.env.state + for (t, s) in enumerate(traj) + @test s.next_sales == DR.sales_history(final_state)[t, :] + @test s.customers == DR.customer_history(final_state)[t] + end + @test isapprox(sum(s.reward for s in traj), total_reward; rtol=1e-5) + @test total_reward <= ant_obj +end + +@testset "DynamicReplenishment - Parametric Anticipative Solver" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=3) + rng = Xoshiro(42) + env = generate_environments(b, 1; seed=0) + model = generate_statistical_model(b) + x, _ = observe(env[1]) + θ = model(x) + param_ant_solver = DR.generate_parametric_anticipative_solver(b) + ant_traj = param_ant_solver(θ, env[1].env.scenario, env[1]) + policies = generate_baseline_policies(b) + r_greedy, greedy_traj = evaluate_policy!(policies.greedy, env) + @test length(ant_traj) == DR.max_steps(b) == length(greedy_traj) + for g_sample in greedy_traj + @test DR.is_feasible(g_sample.state, g_sample.y) + end + for ant_sample in ant_traj + @test DR.is_feasible(ant_sample.state, ant_sample.y) + end +end + +@testset "DynamicReplenishment - Θ ̇g(y) = obj(co_layer)" begin + b = DynamicReplenishmentBenchmark(; N=5, max_steps=5) + rng = Xoshiro(42) + env = generate_environments(b, 1) + policies = generate_baseline_policies(b) + _, traj = evaluate_policy!(policies.greedy, env) + maximizer = generate_maximizer(b) + model = generate_statistical_model(b) + x_1, state_1 = traj[1].x, traj[1].state + N = DR.item_count(b) + ub = DR.ub_per_item(state_1) + + Θ = model(x_1) + Y_oracle = maximizer(Θ; state=state_1) + Z_oracle = DR.get_z_from_y(Y_oracle, state_1) + maximizer_obj = DR._obj_function(N, ub, Θ, Y_oracle, Z_oracle) + Θ_dot_g_y = DR.dot(Θ, DR.g(Y_oracle; state=state_1)) + @test isapprox(maximizer_obj, Θ_dot_g_y; rtol=1e-5) +end + +@testset "DynamicReplenishment - Plots" begin + using Plots + + b = DynamicReplenishmentBenchmark() + envs = generate_environments(b, 5) + policies = generate_baseline_policies(b) + _, traj = evaluate_policy!(policies.greedy, envs) + + @test has_visualization(b) + fig1 = plot_sample(b, traj[1]) + @test fig1 isa Plots.Plot + fig2 = plot_trajectory(b, traj) + @test fig2 isa Plots.Plot +end diff --git a/test/runtests.jl b/test/runtests.jl index b2f5072b..c2120199 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,4 +27,5 @@ using Random include("dynamic_vsp_plots.jl") end include("dynamic_assortment.jl") + include("replenishment.jl") end