1 Introduction↩︎

Many real-time decision systems share a common abstraction: at each epoch the decision maker observes contextual covariates, commits to an action under uncertainty, and carries the consequences into future decisions. The same structure appears across many real-time operations. In real-time electricity dispatch, system operators clear security-constrained economic dispatch every five minutes to rebalance generation against fluctuating load and prices [1]. In ride-hail matching, platforms assign drivers to riders as requests arrive under demand and travel-time uncertainty [2]. In dynamic pricing, sellers revise prices as demand and changing circumstances are observed over time, where digital technologies enabled continuously adjusting prices in many business contexts [3]. These systems run at large scale under tight latency budgets. For instance, Uber recorded roughly \(11.3\) billion trips in 2024 [4], and Alibaba reported a peak of \(583{,}000\) orders created per second during its 2020 Singles’ Day event [5]. When a high-quality decision cannot be produced in time, operators fall back on simpler heuristics that forgo the value of optimization, and at these volumes even small per-order losses accumulate rapidly.

Sequential contextual stochastic programs formalize this setting by representing the underlying uncertainty through a conditional distribution based on the covariate information. For the solution of these problems, a standard approach approximates this distribution by a finite scenario set and solves the resulting stochastic program under these scenarios at every epoch by utilizing an off-the-shelf optimization solver. While this delivers high decision quality offline, the scenario-based reformulation typically runs in seconds to minutes per instance depending on its size, falling short of the sub-second responses that peak periods demand, which require the development of novel computationally efficient solution techniques.

Omnichannel order fulfillment in large e-commerce networks is a representative problem setting where the fast response times to satisfy the order assignments are critical for achieving effective operations of these systems [6]. Specifically, each arriving order must be assigned to a distribution center and carrier service that together balance shipping cost, delivery-time risk, multi-unit consolidation discounts, and the opportunity cost of consuming inventory needed by future orders. Leveraging context information which can include factors such as order characteristics, state of the distribution and carrier network, weather information and current demand, can be critical in this regard to represent the underlying uncertainties such as delivery times accurately. To this end, although contextual stochastic optimization (CSO) has lowered the realized cost of these decisions by integrating distributional forecasts of delivery times [7], extending CSO to sequential online environments remains computationally challenging [8].

Concretely, in the JD.com transactional setting used in this paper for the order fulfillment problem [9], peak hours of a typical day see hundreds of order arrivals per minute, with even higher rates during holiday seasons, while the contextual sample average approximation (C-SAA) reformulation [10] typically takes seconds to minutes per order to solve to a tight optimality gap [7]. Solver-based C-SAA is therefore unavailable when peak periods demand sub-second responses, and operational systems fall back on heuristics that ignore the modeled uncertainty.

This paper develops a learning-based optimization proxy for sequential contextual stochastic programs and specializes it to omnichannel order fulfillment. The proxy is a scenario-embedded hierarchical neural network trained offline on solver-based labels, paired at deployment with a dedicated decoder that returns a feasible plan; together they replace the per-epoch solve with a single forward pass. Unlike prior optimization proxies that do not explicitly condition on the contextual uncertainty distribution [11]–[15], the proposed architecture is explicitly scenario-aware, capturing both immediate stochastic costs and the sequential impact of resource consumption.

The contributions of this paper are four-fold:

1. A scenario-conditioned optimization proxy for sequential contextual stochastic programs: This model maps contextual features and scenarios to decisions in a single forward pass, ensuring feasibility via a lightweight decoder and utilizing a multi-term loss function that integrates solver imitation, constraint penalties, and self-supervised cost-alignment.

2. A two-stage contextual stochastic model for sequential omnichannel fulfillment: This formulation explicitly ties each arriving order’s assignment to the inventory available for future orders. To address this problem, a tailored proxy is introduced, featuring a hierarchical neural network over distribution centers and carrier services that generates scenario-conditioned scores. An inventory-weighted decoder then converts these scores into a feasible plan in a single pass.

3. Theoretical guarantees for the fulfillment proxy: The proposed inventory-weighted decoder is shown to yield a feasible decision in time log-linear in the size of the fulfillment network (the number of distribution centers and carrier services) per product line. Furthermore, its expected cost upper-bounds the optimal two-stage objective, and its out-of-sample value admits unbiased, consistent sample-mean and sample-variance estimators.

4. Empirical evaluation via an industry-calibrated simulator: Using a simulation framework built from JD.com transactional records and proprietary-calibrated carrier services, the operational impact of the proxy is demonstrated. Compared to an online solver-based baseline utilizing the same fitted forecaster, the proposed approach reduces per-order decision latency by approximately \(2{,}800\times\) while lowering realized fulfillment costs by \(3.3\%\). Beyond computational efficiency and better out-of-sample performance, the proxy improves operational outcomes. Relative to established fulfillment policies, it lowers total realized fulfillment cost by at least \(10.7\%\) and roughly halves the late-delivery rate of the strongest classical baselines. Fewer late deliveries raise on-time service, the dimension that drives customer satisfaction.

The remainder of the paper is organized as follows. Section 2 reviews related work. Section 3 formalizes the two-stage contextual stochastic model and its C-SAA reformulation. Section 4 develops the scenario-embedded proxy and instantiates it for omnichannel fulfillment. Section 5 reports the computational study on JD.com data. Section 6 concludes the paper and discusses future directions.

2 Literature Review↩︎

This section surveys three relevant domains which this paper provides contributions: contextual stochastic optimization, machine learning for stochastic and combinatorial optimization, and dynamic e-commerce fulfillment.

2.1 Contextual Stochastic Optimization↩︎

Classical stochastic programming typically represents uncertainty with an empirical distribution estimated from historical data, ignoring any auxiliary information that might be available. In contrast, contextual stochastic optimization (CSO) uses concurrent auxiliary information, or covariates, to construct conditional distributions and thereby produce context-dependent decisions [8]. CSO methodologies fall broadly into three streams: sequential predict-and-optimize, integrated learning and optimization, and direct decision-rule learning. Sequential predict-and-optimize methods train a predictive model and pass its output to a downstream solver [16]. Integrated approaches, or decision-focused learning, instead train predictive models directly with respect to downstream decision quality [17], [18]. However, in settings with non-linear and asymmetric costs, such as e-commerce fulfillment with severe late-delivery penalties, relying only on point forecasts may still be insufficient.

To address this limitation, recent CSO methods increasingly use conditional scenario generation to represent predictive uncertainty [19]–[21]. In omnichannel fulfillment, [7] established the benefit of distribution-aware C-SAA relative to point-forecast and empirical-SAA baselines when delivery uncertainty is modeled explicitly. The present paper extends this from a static offline setting to an online two-stage contextual stochastic program that also captures the future opportunity cost of inventory consumption. To satisfy tight latency requirements, the proposed framework adopts direct decision-rule learning, training a parameterized policy that maps contextual states directly to fulfillment decisions. While related to sequential control, direct reinforcement learning is less attractive here due to hard operational constraints and the availability of offline optimization labels.

2.2 Machine Learning for Stochastic and Combinatorial Optimization↩︎

The intersection of machine learning and operations research has yielded diverse methodologies for combinatorial and stochastic optimization, ranging from accelerating exact heuristics [22] to differentiable categorical relaxations [23]. In multi-stage stochastic programming, learning-based proxies address intractability primarily through direct solution prediction, value function approximation, and optimization-integrated learning [24]. For instance, Neur2SP [25] approximates the expected second-stage value function by pre-training a neural network that is subsequently encoded as fixed constraints in a downstream mixed-integer program. [26] avoids pre-training by embedding the second-stage cost approximation directly into the optimization formulation alongside the first-stage decisions, yielding new theoretical optimality bounds.

A directly relevant thread is constrained optimization learning [15], which imitates complex optimization oracles for sub-second inference in autonomous fleet control and vehicle relocation [11], [27], recurrent dual-sourcing inventory control [28], and end-to-end load planning [14]. However, these works primarily study deterministic settings or rely on point-forecast approximations, leaving scenario-aware policies for sequential contextual stochastic programs largely unexplored. The proposed proxy differs from these threads on three axes. First, a permutation-invariant DeepSets-style encoder consumes a sampled scenario set as a first-class input rather than a point forecast. Furthermore, solver-label imitation is augmented by a constraint-violation penalty and a self-supervised cost-alignment term, rather than the decision-focused gradients of the respective learning setting [17], [18]. Additionally, decisions are first-stage and per-epoch in a sequential inventory-coupled setting, rather than the static second-stage value-function approximations of [25] and [26].

2.3 Dynamic E-Commerce Fulfillment↩︎

Omnichannel fulfillment requires assigning orders to facilities and carrier services to balance costs, consolidation, and service levels [6], [29]. While rolling-horizon policies and dynamic inventory heuristics improve myopic decisions by addressing future demand uncertainty and capacity limits [30]–[32], they typically assume deterministic delivery times. Conversely, recent data-driven studies emphasize timeliness by explicitly modeling delivery uncertainty [7], [33]–[35]. Because existing models have yet to integrate both uncertainty streams simultaneously, formulating a sequential model that jointly captures future demand and delivery uncertainties constitutes a primary contribution of this paper. Beyond these modeling contributions, the paper develops a learning-based solution approach for this problem. A scenario-aware optimization proxy returns fulfillment decisions in a single forward pass and runs several orders of magnitude faster than solver-based stochastic optimization. Existing order-fulfillment methods do not develop such a proxy, so the need for fast, scenario-aware solutions in this setting remains as an open challenge.

3 Problem Formulation↩︎

Many e-commerce firms face an online order-fulfillment problem in their daily operations. Customers place orders continuously throughout the day, and each order must be served quickly under tight delivery promises. An order consists of one or more products, each requested in some quantity, together with the customer’s destination and order attributes such as product categories, brands, order type, and applicable discounts. When an order arrives, the firm observes its contents, the on-hand inventory at each distribution center, the available carrier services and their shipping costs, and contextual signals predictive of delivery times and future demand. The firm must then immediately decide, for each product, which distribution center fulfills it and which carrier service ships it. Because these decisions consume inventory that later orders also draw on, each assignment shapes the cost and feasibility of fulfilling future orders. The remainder of this section formalizes this setting.

The online order fulfillment is formulated over a planning horizon \(\mathcal{T} = \{1, 2, \ldots, T\}\), where each epoch \(t \in \mathcal{T}\) corresponds to the arrival of a customer order. The fulfillment network consists of a set of locations/distribution centers (DCs) \(\mathcal{J}\), a set of products \(\mathcal{I}\), and a set of carrier-services \(\mathcal{K}\). Each order requests quantities \(\mathbf{q}^t = (q^t_i)_{i \in \mathcal{I}}\). The current on-hand inventory levels are denoted by \(\mathrm{Inv}^t_{i,j}\) for each product \(i \in \mathcal{I}\) at each location \(j \in \mathcal{J}\). The core task is to create a dynamic fulfillment plan for each arriving order, assigning products to locations and carrier-services to minimize total expected costs while respecting operational constraints (see App. 7 for a full nomenclature).

For each order, the primary operational decision is the fulfillment plan \(\mathbf{z}^{t} = (z^t_{i,j,k})_{i \in \mathcal{I},\, j \in \mathcal{J},\, k \in \mathcal{K}}\), where \(z^t_{i,j,k}\) is the integer quantity of product \(i\) sourced from location \(j\) and shipped using carrier-service \(k\). An auxiliary slack variable \(\mathbf{u}^t = (u^t_i)_{i \in \mathcal{I}}\) represents unfulfilled demand for each product, ensuring feasibility when the requested quantity exceeds globally available inventory. The joint decision at epoch \(t\) then belongs to the feasible set \(\mathcal{Z}^t\) defined by fulfillment and inventory restrictions: \[\mathcal{Z}^t = \left\{ (\mathbf{z}^t, \mathbf{u}^t) \in \mathbb{Z}^{|\mathcal{I}||\mathcal{J}||\mathcal{K}|}_{+} \times \mathbb{Z}^{|\mathcal{I}|}_{+} : \sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} z^t_{i,j,k} + u^t_i = q^t_i, \;\forall i \in \mathcal{I}; \; \sum_{k \in \mathcal{K}} z^t_{i,j,k} \le \mathrm{Inv}^t_{i,j}, \;\forall i \in \mathcal{I},\, j \in \mathcal{J} \right\}.\]

The decision-making process is subject to two primary sources of uncertainty. The first is delivery timeliness. Let \(\tilde{\boldsymbol{\Delta}}^t = (\tilde{\Delta}^t_{j,k})_{j \in \mathcal{J},\, k \in \mathcal{K}}\) be a random vector, where \(\tilde{\Delta}^t_{j,k}\) represents the deviation in days from the promised delivery date when an item is sourced from location \(j\) and shipped using carrier-service \(k\). A positive realization \(\Delta^t_{j,k}\) indicates a late delivery, while a negative value indicates an early one. The realization of \(\Delta^t_{j,k}\) for the chosen location-carrier pairs is observed only after the decision \(\mathbf{z}^t\) is executed, analogous to a partial feedback setting [7]. The second source of uncertainty is future demand, represented by a random vector \(\tilde{\mathbf{D}}^t = (\tilde{D}^t_i)_{i \in \mathcal{I}}\) denoting the cumulative demand for each product from \(t+1\) to \(T\).

The per-epoch cost function \(C(\mathbf{z}^t, \mathbf{u}^t, \tilde{\boldsymbol{\Delta}}^t)\) combines deterministic shipping costs, penalties for stochastic delivery deviations, and a unit stockout penalty \(\rho\) applied to unfulfilled demand. To encourage consolidation, a discount \(\beta \in (0,1)\) is applied to the shipping cost \(c_{i,j,k}^{\mathrm{ship}}\) for any location-carrier pair \((j,k)\) that sources two or more units for the order. Asymmetric penalties, \(\gamma^+ \gg \gamma^- > 0\), are applied for late (\(\tilde{\Delta}_{j,k}^{\,t+} = \max\{0, \tilde{\Delta}^t_{j,k}\}\)) and early (\(\tilde{\Delta}_{j,k}^{\,t-} = \max\{0, -\tilde{\Delta}^t_{j,k}\}\)) deliveries. The cost \(C(\mathbf{z}^t, \mathbf{u}^t, \tilde{\boldsymbol{\Delta}}^t)\) is given by \[\sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} \left[ \left( \sum_{i \in \mathcal{I}} c_{i,j,k}^{\mathrm{ship}} z^t_{i,j,k} \right) \left( 1 - \beta \, \mathbf{1}\Bigl\{\sum_{i \in \mathcal{I}} z^t_{i,j,k} \ge 2\Bigr\} \right) + \left( \gamma^+ \tilde{\Delta}_{j,k}^{\,t+} + \gamma^- \tilde{\Delta}_{j,k}^{\,t-} \right) \sum_{i \in \mathcal{I}} z^t_{i,j,k} \right] + \rho \sum_{i \in \mathcal{I}} u^t_i.\] As this cost function induces non-linear terms, in the computational implementation, the non-linear indicator function governing the consolidation discount is linearized using standard big-\(M\) constraints and auxiliary binary variables.

At the start of epoch \(t\), a contextual vector \(\mathbf{X}^t \in \mathcal{X}\) is observed, containing features of the order, customer, and current network state. Because fulfillment decisions consume inventory and thereby affect future feasible actions, the problem is inherently sequential. The system state therefore evolves over time as inventory is depleted and new exogenous information arrives.

3.1 Two-Stage Contextual Stochastic Model↩︎

The sequential nature of online fulfillment gives rise to a multi-stage stochastic program that is typically intractable. A rolling-horizon framework with stage aggregation is therefore employed. At each epoch \(t\), a two-stage contextual stochastic program is solved for the immediate decision \((\mathbf{z}^t, \mathbf{u}^t)\), collapsing future epochs \(t+1, \dots, T\) into a single second-stage problem [36]. The first-stage decision is committed before any current-period or future uncertainty is realized. The second-stage variables \(\mathbf{v}^t = (v^t_{i,j,k})_{i, j, k}\) and \(\mathbf{w}^t = (w^t_i)_{i \in \mathcal{I}}\) represent the aggregated fulfillment plan and slack stockout for orders from \(t+1\) to \(T\). The problem solved at each epoch is \[\label{eq:full} \min_{\;(\mathbf{z}^t, \mathbf{u}^t) \in \mathcal{Z}^t} \; \mathbb{E}_{\tilde{\boldsymbol{\Delta}}^t,\, \tilde{\mathbf{D}}^t \sim \mathbb{P}(\cdot \mid \mathbf{X}^t)} \Big[ C\!\big(\mathbf{z}^t, \mathbf{u}^t, \tilde{\boldsymbol{\Delta}}^t\big) + Q\!\big(\mathbf{z}^t, \tilde{\mathbf{D}}^t\big) \Big],\tag{1}\] where \(\mathbb{P}(\cdot \mid \mathbf{X}^t)\) represents the conditional distribution given the context \(\mathbf{X}^t\) and \(Q(\mathbf{z}^t, \tilde{\mathbf{D}}^t)\) is the second-stage cost function representing the expected future shipping and stockout costs: \[Q(\mathbf{z}^t, \tilde{\mathbf{D}}^t) = \min_{\mathbf{v}^{t} \ge 0,\, \mathbf{w}^t \ge 0} \left\{ \sum_{i \in \mathcal{I}} \sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} c^{\mathrm{ship}}_{i,j,k} \, v^t_{i,j,k} + \rho \sum_{i \in \mathcal{I}} w^t_i \;\middle|\; \begin{align} &\sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} v^t_{i,j,k} + w^t_i = \tilde{D}^t_i, && \forall i \in \mathcal{I}, \\ &\sum_{k \in \mathcal{K}} \bigl(z^t_{i,j,k} + v^t_{i,j,k}\bigr) \le \mathrm{Inv}^t_{i,j}, && \forall i \in \mathcal{I},\, j \in \mathcal{J}. \end{align} \right\}.\] The first-stage decision \(\mathbf{z}^t\) enters the second stage through the post-allocation inventory state, which constrains the feasible recourse and thereby determines \(Q\). The second-stage variable \(\mathbf{v}^t\) is declared continuous: with integer right-hand sides \(\mathbf{Inv}^t\) and \(\tilde{\mathbf{D}}^t\), the transportation-network constraint matrix is totally unimodular and an integer optimum is recovered without imposing integrality. The second-stage cost includes only base shipping and stockout costs on the forecast future demand; it omits future delivery-time penalties and consequently considers an approximation of the true future cost, since the aggregate delivery-time deviation across the remaining horizon cannot be forecast at decision time.

3.2 Contextual Sample Average Approximation (C-SAA)↩︎

To solve the stochastic program 1 , the conditional distribution \(\mathbb{P}(\cdot \mid \mathbf{X}^t)\) is approximated by a finite sampled scenario set \(\tilde{\Omega}^t = \{\omega_1, \dots, \omega_E\}\) given the context \(\mathbf{X}^t\), with each scenario \(\omega\) providing a delivery-time realization \(\boldsymbol{\Delta}^t_{\omega} = (\Delta^t_{j,k,\omega})_{j,k}\) and per-product cumulative demands \(D^t_{i,\omega}\). The details of the scenario generation process and integration of the context information are discussed in Section 4.2.1. By introducing per-scenario second-stage variables \(\mathbf{v}^t_{\omega} = (v^t_{i,j,k,\omega})\) and \(\mathbf{w}^t_{\omega} = (w^t_{i,\omega})\), the C-SAA extensive form can be written as follows [7], [10]: \[\label{eq:csaa95extensive} \begin{align} \min_{\substack{(\mathbf{z}^t, \mathbf{u}^t) \in \mathcal{Z}^t \\ \mathbf{v}^t_{\omega} \ge 0,\; \mathbf{w}^t_{\omega} \ge 0}} \quad & \frac{1}{E} \sum_{\omega \in \tilde{\Omega}^t} \Bigg( C(\mathbf{z}^t, \mathbf{u}^t, \boldsymbol{\Delta}^t_{\omega}) + \sum_{i \in \mathcal{I}} \sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} c^{\mathrm{ship}}_{i,j,k} \, v^t_{i,j,k,\omega} + \rho \sum_{i \in \mathcal{I}} w^t_{i,\omega} \Bigg) \\ \text{s.t.} \quad & \sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} v^t_{i,j,k,\omega} + w^t_{i,\omega} = D^t_{i,\omega}, && \forall i \in \mathcal{I},\, \omega \in \tilde{\Omega}^t, \\ & \sum_{k \in \mathcal{K}} \bigl(z^t_{i,j,k} + v^t_{i,j,k,\omega}\bigr) \le \mathrm{Inv}^{t}_{i,j}, && \forall i \in \mathcal{I},\, j \in \mathcal{J},\, \omega \in \tilde{\Omega}^t . \end{align}\tag{2}\]

As a single finite-sample C-SAA instance is sensitive to scenario sampling variability, an iterative evaluation procedure is utilized to generate multiple candidate solutions under different scenario sets. The proposed procedure follows the standard multiple-replication strategy of the SAA algorithms [10], where \(S\) candidate first-stage decisions are generated by solving independent instances of 2 with scenario size \(N_1\), then these solutions are evaluated on a larger common scenario set of size \(N_2 \gg N_1\) to estimate their out-of-sample performances to select the best candidate solution. Details of this procedure is provided in Appendix 8.

4 A Scenario-Embedded Optimization Proxy for Sequential Contextual Stochastic Program↩︎

Solving a solver-based scenario approximation of the two-stage contextual stochastic program could still be computationally prohibitive in real time during peak operating periods. This motivates the low-latency optimization proxy pipeline developed in this section to address these operational challenges. Section 4.1 states this design as a general proxy framework for sequential contextual stochastic programs. Section 4.2 instantiates the framework for the fulfillment setting of Section 3.

4.1 General Proxy Framework↩︎

Consider a sequential contextual stochastic program in which, at each epoch \(t\), the decision maker observes a state or context \(\mathbf{X}^t\), samples a scenario set \(\mathbf{\Omega}^t\) from a fitted conditional forecaster, and chooses a first-stage decision in a feasible set \(\mathcal{Z}^t\) before the next epoch. The computational bottleneck arises when a solver-based approximation, such as SAA or another scenario-based reformulation, can produce high-quality first-stage decisions offline but is too slow to run at every online epoch. The proxy framework replaces this online solve with a learned decision rule trained on offline solver labels.

The framework has four ingredients. First, an offline oracle maps historical contexts and scenario sets to first-stage action labels \(\mathbf{a}^{t,*}\). Second, a scenario-embedded model \(f_\theta(\mathbf{X}^t,\mathbf{\Omega}^t)\) maps the same information to decision scores or logits. Third, a decoder \(D\) maps these scores into a feasible decision in \(\mathcal{Z}^t\). Fourth, training combines label imitation with problem-aware regularization, such as a soft feasibility penalty and a scenario-based cost-alignment term. Online, the solver call is replaced by one forward pass through \(f_\theta\) followed by the decoder \(D\).

This framework is intended for latency-critical settings where high-quality offline labels are available and feasibility can be enforced by a decoder, projection, or efficient repair rule. The quality of the solutions obtained through this approach depends on the offline label generator, the fitted scenario model, and the expressiveness of the proxy class, which can be specialized depending on the problem setting considered. To this end, the fulfillment application below is one instantiation of the proposed proxy framework, where the offline oracle is finite-sample C-SAA, the generic action label \(\mathbf{a}^{t,*}\) corresponds to the fulfillment decision \((\mathbf{z}^{t,*},\mathbf{u}^{t,*})\), the decoder is inventory-weighted feasible inference, and the scenario set contains demand and delivery-time samples.

4.2 Fulfillment Instantiation↩︎

The remainder of this section instantiates the framework for omnichannel fulfillment. Figure 1 schematizes the pipeline: demand and delivery-time forecasters feed finite-sample C-SAA on historical orders offline, and the trained proxy serves each arriving order online.

4.2.1 Contextual Inputs and Scenario Generation↩︎

The fulfillment context decomposes as \(\mathbf{X}^t = \bigl(\mathbf{X}_{\mathrm{ord}}^t,\; \{\mathbf{X}_{\mathrm{dc},j}^t\}_{j \in \mathcal{J}},\; \{\mathbf{X}_{\mathrm{cost},j,k}^t\}_{j \in \mathcal{J}, k \in \mathcal{K}}\bigr).\) The order-context vector \(\mathbf{X}_{\mathrm{ord}}^t\) concatenates SKU and brand categorical embeddings with numeric order fields (number of SKUs in the order, total quantity, order type, bundle and coupon discount flags, gift-item flag, average discount rate), calendar features (e.g.hour and day of week), and customer demographics; the per-DC vector \(\mathbf{X}_{\mathrm{dc},j}^t\) collects the on-hand inventory \(\mathrm{Inv}_{j}^t\) and a DC categorical embedding for each candidate \(j \in \mathcal{J}\); and the per-(DC, carrier) vector \(\mathbf{X}_{\mathrm{cost},j,k}^t\) collects the base shipping cost \(c^{\mathrm{ship}}_{j,k}\) and a carrier categorical embedding for each option \((j,k)\).

Two fitted forecasters approximate the conditional distribution \(\mathbb{P}(\cdot \mid \mathbf{X}^t)\) through predicted quantiles, generating a paired scenario set \(\mathbf{\Omega}^t = (\mathbf{\Omega}^t_{\mathrm{dem}},\,\mathbf{\Omega}^t_{\mathrm{time}})\) of \(E\) Monte Carlo scenarios per random vector by inverse-transform sampling: for each \(\omega \in \{1, \ldots, E\}\), a uniform draw \(u_\omega \sim U(0,1)\) is mapped through the linearly interpolated predicted CDF to a realization. A multi-horizon quantile recurrent network (MQRNN) [37] consumes order-context features and historical demand to produce non-crossing per-period demand quantiles, which are inverse-sampled per period and summed across the horizon to yield \(\mathbf{\Omega}^t_{\mathrm{dem}} = (\tilde{\Omega}^t_{i, \mathrm{dem}})_{i \in \mathcal{I}}\) with \(\tilde{\Omega}^t_{i, \mathrm{dem}} = \{D_{i,1}^t, \ldots, D_{i,E}^t\}\). A multi-quantile multilayer perceptron (MQMLP) consumes order-context and per-(DC, carrier) cost features to produce delivery-time quantiles, sampled the same way to yield \(\mathbf{\Omega}^t_{\mathrm{time}} = (\tilde{\Omega}^t_{j,k, \mathrm{time}})_{j \in \mathcal{J}, k \in \mathcal{K}}\) with \(\tilde{\Omega}^t_{j,k, \mathrm{time}} = \{\Delta_{j,k,1}^{\,t}, \ldots, \Delta_{j,k,E}^{\,t}\}\). The online supplement reports the full feature lists, network architectures, training objectives, and held-out forecasting performance.

4.2.2 Scenario-Embedded Hierarchical Proxy↩︎

Because the constraints in \(\mathcal{Z}^t\) are separable across products, a single shared network \(f_\theta\) processes one product line at a time and the order-level decision concatenates per-product outputs. The consolidation discount couples the multiple product lines of the same order that ship together but is not separable; the per-product decomposition therefore sees this coupling only indirectly through the jointly optimized C-SAA labels. Each training instance \(n \in \mathcal{N}\) is one historical (order, product-line) record, paired with its contextual features, sampled scenarios, and C-SAA label.

4.2.2.1 Architecture Walkthrough.

Figure 2 schematizes the proxy. The upper panel shows the architecture: four input modules (global, scenario, DC, cost) produce embeddings that are concatenated per DC and passed to a DC head and a DC-conditioned carrier head, which together return \((\boldsymbol{\ell}^{\mathrm{dc}}, \boldsymbol{\ell}^{\mathrm{carrier}})\). The lower panels summarize the composite training loss applied at training time and the inventory-weighted decoder applied at inference (developed in Sections 4.2.3 and 4.2.4). Figure 3 formalizes the forward pass of the optimization proxy. Each module \(\phi\) and head \(g\) is a multilayer perceptron (MLP) whose hidden layers apply a fully connected linear projection, batch normalization, dropout regularization, and a ReLU activation in sequence; hidden widths and layer counts are tuned and reported in Appendix 13.

None

4.2.2.2 Global and DC modules.

A shared MLP \(\phi_{\mathrm{ord}}\) encodes the order context \(\mathbf{X}_{\mathrm{ord}}\) into a global embedding; a second shared MLP \(\phi_{\mathrm{dc}}\) encodes each per-DC feature vector \(\mathbf{X}_{\mathrm{dc},j}\) into a per-DC embedding, parameter-shared across candidates \(j\) so the embedding size is independent of \(|\mathcal{J}|\): \[\begin{align} \mathbf{e}_{\mathrm{global}} &= \phi_{\mathrm{ord}}(\mathbf{X}_{\mathrm{ord}}), \quad \mathbf{e}_{\mathrm{dc},j} = \phi_{\mathrm{dc}}(\mathbf{X}_{\mathrm{dc},j}), \qquad j \in \mathcal{J}. \end{align}\]

4.2.2.3 Scenario module.

The sampled scenario sets \(\mathbf{\Omega}_{\mathrm{dem}}\) and \(\mathbf{\Omega}_{\mathrm{time}}\) are unordered and may vary in size at deployment. For the current product line, identify \(\mathbf{\Omega}_{\mathrm{dem}} = \{D_\omega\}_{\omega=1}^E\) as the scalar demand scenarios and \(\mathbf{\Omega}_{\mathrm{time}} = \{\boldsymbol{\Delta}_\omega\}_{\omega=1}^E\) as the deviation matrices \(\boldsymbol{\Delta}_\omega = (\Delta_{j,k,\omega})_{j \in \mathcal{J}, k \in \mathcal{K}}\) under a fixed \((j,k)\) ordering. Following the DeepSets template [25], [38], each scenario is encoded individually and the resulting embeddings are averaged. The encoder \(\phi_{\mathrm{scen}}^{\mathrm{dem}}: \mathbb{R} \to \mathbb{R}^h\) (with \(h\) the shared embedding dimension across the proxy’s modules, a tuned hyperparameter reported in Appendix 13) maps \(D_\omega\) to an embedding capturing where the scenario sits in the demand distribution. The encoder \(\phi_{\mathrm{scen}}^{\mathrm{time}}: \mathbb{R}^{|\mathcal{J}||\mathcal{K}|} \to \mathbb{R}^h\) maps the flattened matrix \(\operatorname{vec}(\boldsymbol{\Delta}_\omega)\) to an embedding capturing which DC–carrier options are jointly stressed in that sampled future. Averaging the per-scenario embeddings is the network’s analogue of an expectation and matches the arithmetic-mean structure of the C-SAA objective: \[\begin{align} \mathbf{m}^{\mathrm{dem}} &= \frac{1}{E}\sum_{D_\omega \in \mathbf{\Omega}_{\mathrm{dem}}} \phi_{\mathrm{scen}}^{\mathrm{dem}}(D_\omega), \quad \mathbf{m}^{\mathrm{time}} = \frac{1}{E}\sum_{\boldsymbol{\Delta}_\omega \in \mathbf{\Omega}_{\mathrm{time}}} \phi_{\mathrm{scen}}^{\mathrm{time}}(\operatorname{vec}(\boldsymbol{\Delta}_\omega)), \\ \mathbf{e}_{\mathrm{scen}} &= \mathbf{m}^{\mathrm{dem}} + \mathbf{m}^{\mathrm{time}}. \end{align}\] The two branches use independent weights at the same hidden width, so addition preserves disentangled demand and delivery-time signals without a learned fusion layer.

4.2.2.4 Cost module.

The candidate option set \(\{(j,k)\}_{j,k}\) has \(|\mathcal{J}|\,|\mathcal{K}|\) entries, so feeding all per-option features into the DC head directly would either inflate parameter count or require a second permutation-invariant pool over carriers. Instead, the proxy summarizes the per-DC cost and timeliness features through \(\mathrm{stats}(\mathbf{c}_j, \boldsymbol{\Delta}_j) \in \mathbb{R}^8\), with \(\mathbf{c}_j = (c^{\mathrm{ship}}_{j,k})_{k \in \mathcal{K}}\) and \(\boldsymbol{\Delta}_j = (\Delta_{j,k,\omega})_{k \in \mathcal{K},\, \omega \in [E]}\). The first five outputs are carrier-axis statistics of \(\mathbf{c}_j\) (\(\overline{\mathbf{c}_j}, \min \mathbf{c}_j, \sigma(\mathbf{c}_j), p_{90}(\mathbf{c}_j), \mathrm{gap}_2(\mathbf{c}_j) = c^{(2)}_j - c^{(1)}_j\)) corresponding to mean, minimum, standard deviation, 90th percentile, the gap between the two smallest entries, respectively, and the last three are (carrier, scenario)-axis statistics of \(\boldsymbol{\Delta}_j\) (\(\overline{\boldsymbol{\Delta}_j}, \sigma(\boldsymbol{\Delta}_j), p_{90}(\boldsymbol{\Delta}_j)\)). The cost embedding is \(\mathbf{e}_{\mathrm{cost},j} = \phi_{\mathrm{cost}}(\mathrm{stats}(\mathbf{c}_j, \boldsymbol{\Delta}_j))\). These dispersion statistics expose per-carrier variability to the DC head without committing to a specific carrier.

4.2.2.5 Joint embedding and hierarchical heads.

The four embeddings concatenate per DC, \(\mathbf{h}_j = [\mathbf{e}_{\mathrm{global}} \parallel \mathbf{e}_{\mathrm{scen}} \parallel \mathbf{e}_{\mathrm{dc},j} \parallel \mathbf{e}_{\mathrm{cost},j}]\) for \(j \in \mathcal{J}\). The DC head, MLP \(g_{\mathrm{dc}}\), maps \(\mathbf{h}_j\) to a scalar logit \(\ell^{\mathrm{dc}}_j = g_{\mathrm{dc}}(\mathbf{h}_j)\). The carrier head \(g_{\mathrm{carrier}}\) then conditions the carrier choice on the joint embedding, the raw DC logit, and the per-option features: \[\label{eq:proxy95ell95carrier} \ell^{\mathrm{carrier}}_{j,k} \;=\; g_{\mathrm{carrier}}\!\left(\mathbf{h}_j,\; \ell^{\mathrm{dc}}_j,\; \mathbf{X}_{\mathrm{cost},j,k}\right), \qquad j \in \mathcal{J},\, k \in \mathcal{K}.\tag{3}\] Two design choices matter in 3 . First, the carrier head receives the raw, pre-softmax DC logit, so carrier selection can depend on the strength of the learned DC ranking rather than only on the relative ordering of DCs. Second, the two heads imply a hierarchical factorization \(p(j) \cdot p(k \mid j)\). The DC-head softmax realizes \(p(j)\) over distribution centers, and the carrier-head softmax realizes \(p(k \mid j)\) over each DC’s eligible carrier services. This factorization matches the structure of the fulfillment decision: an order line first selects a DC, then a carrier-service from that DC’s eligible set. It reduces the complexity of learning the full joint \(|\mathcal{J}|\,|\mathcal{K}|\) assignment.

4.2.3 Composite Training Loss↩︎

Figure 4 states the training optimization. It takes the \(|\mathcal{N}|\) training instances as input and minimizes the average per-instance composite loss \(\mathcal{L}^n\), which combines three terms targeting C-SAA-label imitation, operational feasibility, and self-supervised cost alignment.

The offline C-SAA oracle returns a quantity decision \((\mathbf{z}^{n,*},u^{n,*})\) for the product line, which may distribute the requested quantity across multiple \((j,k)\) options. The learning task is framed as predicting a single primary option per product line rather than the full multi-DC split, which encourages shipping consolidation; if a shortage forces a split, the inventory-weighted decoder of Section 4.2.4 recovers it at inference. The primary location and carrier-service labels are extracted as the most-served option of the C-SAA decision: \[j^{n,*} \in \arg\max_{j\in\mathcal{J}}\sum_{k\in\mathcal{K}} z^{n,*}_{j,k}, \qquad k^{n,*} \in \arg\max_{k\in\mathcal{K}_{j^{n,*}}^n} z^{n,*}_{j^{n,*},k},\] with deterministic tie breaking.

None

4.2.3.1 Selection Loss.

Let \(\mathbf{p}^{n,\mathrm{dc}} = \mathrm{softmax}(\boldsymbol{\ell}^{n,\mathrm{dc}})\) and \(\mathbf{p}^{n,\mathrm{carrier}} = \mathrm{softmax}_k(\boldsymbol{\ell}^{n,\mathrm{carrier}})\) denote the head probabilities, with the carrier softmax taken along the carrier-service axis at each fixed DC. The selection loss \(\mathcal{L}_{\mathrm{sel}}\) supervises both heads against the primary labels \((j^{n,*}, k^{n,*})\) via the negative log-likelihood of the labeled pair under the hierarchical predicted distribution (equivalently, the categorical cross-entropy on the unnormalized logits), with separate weights on the DC and conditional-carrier terms: \[\mathcal{L}^n_{\mathrm{sel}} \;=\; -\lambda_{\mathrm{sel}} \Bigl(\log p^{n,\mathrm{dc}}_{j^{n,*}} \;+\; \lambda_{\mathrm{carrier}}\, \log p^{n,\mathrm{carrier}}_{j^{n,*},\, k^{n,*}}\Bigr).\] Here \(\lambda_{\mathrm{sel}}\) scales the overall selection weight and \(\lambda_{\mathrm{carrier}}\) scales the relative importance of the conditional carrier signal.

4.2.3.2 Constraint Loss.

The one-hot DC assignment \(\mathbf{a}^n \in \{0,1\}^{|\mathcal{J}|}\) is drawn from the DC logits via the Straight-Through Gumbel-Softmax (ST-GumbelSoftmax) estimator of [23]: \(\mathbf{a}^n\) is the one-hot encoding of \(\arg\max_j (\ell^{n,\mathrm{dc}}_j + g_j)\) with independent and identically distributed (i.i.d.) standard Gumbel noise \(g_j\), but its gradient is computed as if \(\mathbf{a}^n\) equaled the temperature-\(\tau\) softmax \(\mathrm{softmax}(\boldsymbol{\ell}^{n,\mathrm{dc}}/\tau)\), so a discrete DC choice can pass to the loss while gradients still flow back through \(g_{\mathrm{dc}}\). The discrete sample is essential here: evaluating the inventory cap on the soft mass \(\sum_j p^{n,\mathrm{dc}}_j \min\{q^n, \mathrm{Inv}^n_j\}\) instead would let a hedge between an infeasible primary and a feasible secondary DC escape penalty.

Two feasibility requirements in \(\mathcal{Z}^t\) bind on the DC selected by \(\mathbf{a}^n\): the per-DC inventory cap (\(\sum_k z_{j,k} \le \mathrm{Inv}^n_j\)) and the demand-satisfaction requirement that the requested quantity \(q^n\) be served. The constraint penalty \(\mathcal{L}_{\mathrm{const}}\) encodes both as the demand shortfall after each DC’s allocation is capped by its on-hand inventory: \[\mathcal{L}^n_{\mathrm{const}} \;=\; \lambda_{\mathrm{const}} \left[\, q^n \;-\; \sum_{j \in \mathcal{J}} \min\!\bigl\{q^n a^n_j,\; \mathrm{Inv}^n_{j}\bigr\} \,\right]^+,\] where \([x]^+ = \max(0, x)\) and \(\lambda_{\mathrm{const}}\) scales the penalty. Because \(\mathbf{a}^n\) is one-hot, the served quantity reduces to \(\min\{q^n, \mathrm{Inv}^n_{j^*}\}\) at the chosen primary DC \(j^*\), so the penalty fires precisely when that DC’s inventory is insufficient.

4.2.3.3 Cost-Alignment Loss.

The cost-alignment penalty \(\mathcal{L}_{\mathrm{cost}}\) is a self-supervised regularizer that pulls probability mass toward low-cost (DC, carrier-service) pairs. It operates on the head probabilities \((\mathbf{p}^{n,\mathrm{dc}}, \mathbf{p}^{n,\mathrm{carrier}})\) defined above, so the gradient flows through a smooth surface. Letting \(c^n_{j,k,\omega}\) be the realized immediate fulfillment cost for option \((j,k)\) under scenario \(\omega\) (shipping plus delivery-time penalty, no second-stage recourse) drawn from \(\mathbf{\Omega}^n_{\mathrm{cost}}\), the term computes the expected immediate cost under the predicted joint \(p^{n,\mathrm{dc}}_j p^{n,\mathrm{carrier}}_{j,k}\): \[\mathcal{L}^n_{\mathrm{cost}} \;=\; \lambda_{\mathrm{cost}}\, \frac{1}{E} \sum_{\omega=1}^E \sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} p^{n,\mathrm{dc}}_j\, p^{n,\mathrm{carrier}}_{j,k}\, c^n_{j,k,\omega}.\]

4.2.4 Feasible Inference via Inventory-Weighted Decoding↩︎

At an online epoch \(t\), the trained proxy first maps the current context and sampled scenarios to the same hierarchical logits used during training, \[(\boldsymbol{\ell}^{\mathrm{dc}}, \boldsymbol{\ell}^{\mathrm{carrier}}) = f_\theta\!\left(\mathbf{X}_{\mathrm{ord}}^t,\{\mathbf{X}_{\mathrm{dc},j}^t\}_{j\in\mathcal{J}}, \{\mathbf{X}_{\mathrm{cost},j,k}^t\}_{j,k},\mathbf{\Omega}^t_{\mathrm{dem}},\mathbf{\Omega}^t_{\mathrm{time}}\right),\] and converts them to probabilities \(\mathbf{p}^{\mathrm{dc}} = \mathrm{softmax}(\boldsymbol{\ell}^{\mathrm{dc}})\) and \(\mathbf{p}^{\mathrm{carrier}} = \mathrm{softmax}_k(\boldsymbol{\ell}^{\mathrm{carrier}})\). These probabilities are not used directly as fractional decisions. Instead, the decoder reweights each DC’s predicted probability by the fractional coverage that the DC can offer for the current line: for a product line with requested quantity \(q^t\), \[r_j^t=\min\{\mathrm{Inv}_j^t/q^t,1\},\qquad s_j^t=p_j^{\mathrm{dc}} r_j^t,\] where \(r_j^t \in [0,1]\) leaves a fully-covering DC undisturbed and discounts a partially-covering one in proportion to its coverage. The reweighting approximates the second-stage recourse cost \(Q\) that the training loss does not directly model: depleting a DC whose stock barely covers the current line raises the likelihood of a costly split or stockout on a future order that the same DC could have served, so steering volume toward DCs with healthier inventory buffers preserves serviceability over the horizon. The score \(s_j^t\) also softens the single-DC preference learned through \(\mathcal{L}_{\mathrm{const}}\), allowing multi-DC splits when shortages require them.

The decoder then sorts eligible DCs by \(s_j^t\) and scans them greedily. For each visited DC, it chooses the highest-probability eligible carrier service and assigns as much residual demand as inventory permits, where \(\mathcal{K}_j^t\) denotes the set of carrier services eligible at DC \(j\) at epoch \(t\): \[k^*(j)\in\arg\max_{k\in\mathcal{K}_j^t}p_{j,k}^{\mathrm{carrier}},\qquad x_j^t=\min\{q_{\mathrm{rem}}^t,\mathrm{Inv}_j^t\},\qquad z_{j,k^*(j)}^{\mathrm{proxy},t}=x_j^t.\] The residual \(q_{\mathrm{rem}}^t\) (the quantity of the line not yet assigned, initialized to the requested quantity \(q^t\)) is updated after each assignment, and any remaining quantity becomes unmet demand \(u^{\mathrm{proxy},t}\). Thus the learned probabilities determine the ordering and carrier choice, while the decoder enforces the demand-allocation, inventory, and eligibility constraints in \(\mathcal{Z}^t_i\) corresponding to the product-line \(i\) restriction of \(\mathcal{Z}^t\), which can be defined as \(\mathcal{Z}^t_i =\{ (\mathbf{z}^t, \mathbf{u}^t) \in \mathbb{Z}^{|\mathcal{J}||\mathcal{K}|}_{+} \times \mathbb{Z}_{+}: \sum_{j \in \mathcal{J}} \sum_{k \in \mathcal{K}} z^t_{j,k} + u^t = q^t; \; \sum_{k \in \mathcal{K}} z^t_{j,k} \le \mathrm{Inv}^t_{j},\, j \in \mathcal{J}\}\). Algorithm 5 summarizes the end-to-end online inference procedure for one product line.

Proposition 1 establishes feasibility, runtime, and policy-value estimation properties of the proposed framework.

Proof.

Since the constraints in \(\mathcal{Z}^t\) are separable across products with \(\mathcal{Z}^t = \prod_{i \in \mathcal{I}} \mathcal{Z}^t_i\), and the order-level decision concatenates per-product outputs of Algorithm 5, it suffices to verify \((\mathbf{z}^{\mathrm{proxy},t}_i, u^{\mathrm{proxy},t}_i) \in \mathcal{Z}^t_i\) for the product \(i\). Each loop iteration assigns an integer \(x_j = \min\{q_{\mathrm{rem}}, \mathrm{Inv}_j^t\} \le \mathrm{Inv}_j^t\) to a single carrier at a previously unvisited \(j\), so \(\sum_k z^{\mathrm{proxy},t}_{j,k} \le \mathrm{Inv}_j^t\), and the loop invariant \(q^t = q_{\mathrm{rem}} + \sum_j x_j\) gives \(\sum_{j,k} z^{\mathrm{proxy},t}_{j,k} + u^{\mathrm{proxy},t} = q^t\) on termination. Hence \((\mathbf{z}^{\mathrm{proxy},t},\mathbf{u}^{\mathrm{proxy},t}) \in \mathcal{Z}^t\). Then, by utilizing this feasible solution, an upper bound can be constructed as \(f_{\mathrm{proxy}}^t = \mathbb{E}[\,C(\mathbf{z}^{\mathrm{proxy},t},\mathbf{u}^{\mathrm{proxy},t},\tilde{\boldsymbol{\Delta}}^t) + Q(\mathbf{z}^{\mathrm{proxy},t},\tilde{\mathbf{D}}^t) \mid \mathbf{X}^t\,] \ge \min_{(\mathbf{z},\mathbf{u}) \in \mathcal{Z}^t} \mathbb{E}[\,C(\mathbf{z},\mathbf{u},\tilde{\boldsymbol{\Delta}}^t) + Q(\mathbf{z},\tilde{\mathbf{D}}^t)\mid \mathbf{X}^t\,]\).

For runtime, sorting eligible DCs by \(s_j^t\) costs \(O(|\mathcal{J}|\log|\mathcal{J}|)\). The allocation scan over sorted DCs costs \(O(|\mathcal{J}|)\), and selecting \(k^*(j)\) scans the eligible carrier sets at total cost \(\sum_j |\mathcal{K}_j^t|\). Thus, the decoder costs \(O(|\mathcal{J}|\log|\mathcal{J}| + \sum_j |\mathcal{K}_j^t|)\), showing the desired result. When all carriers are eligible at all DCs, \(\sum_j |\mathcal{K}_j^t|=|\mathcal{J}||\mathcal{K}|\). 0◻

To evaluate the policy obtained in practice, an independent out-of-sample scenario set can be utilized. The following remark clarifies the statistical interpretation of this evaluation procedure.

4.3 Training Data Collection and Model Selection↩︎

This Section describes the experimental protocol used to generate proxy-training labels and evaluate out-of-sample policy performance. The training horizon spans three consecutive weeks, partitioned chronologically: the initial two weeks form the forecast training split \(\mathcal{D}_{\mathrm{fcast}}\), the third week constitutes the proxy training split \(\mathcal{D}_{\mathrm{proxy}}\), and a held-out fourth week is used for online evaluation. Forecasters are trained on \(\mathcal{D}_{\mathrm{fcast}}\) alone, with no exposure to \(\mathcal{D}_{\mathrm{proxy}}\) or the evaluation week to prevent leakage into labels or test; they then generate context-conditional scenarios for each order in \(\mathcal{D}_{\mathrm{proxy}}\), and C-SAA solves the resulting instances against the inventory state at each order’s arrival epoch to produce the C-SAA-optimal label \(\mathbf{z}^{n,*}\) (optimal over the finite scenario sample, not the true conditional distribution). The same scenario tensors feed both C-SAA (for labels) and the proxy (as inputs), so the proxy is trained on the uncertainty distribution that produced its labels.

The proxy is then trained on \(\mathcal{D}_{\mathrm{proxy}}\) with hierarchical location-carrier labels \((j^{n,*}, k^{n,*})\) extracted from \(\mathbf{z}^{n,*}\). Before the final evaluation, the forecasters are refitted on \(\mathcal{D}_{\mathrm{train}} = \mathcal{D}_{\mathrm{fcast}} \cup \mathcal{D}_{\mathrm{proxy}}\) for deployment fidelity; the proxy itself is not retrained. A temporal subset of \(\mathcal{D}_{\mathrm{proxy}}\) is held out as a validation set, and the proxy checkpoint is selected by Hit@5, defined as the rate at which the C-SAA-chosen (DC, carrier) pair appears in the proxy’s top-five predicted options. Top-five is used because multiple near-optimal pairs typically exist per instance, and at deployment the inventory-weighted decoder reshuffles top-ranked options against real-time inventory.

5 Computational Study↩︎

The proxy is evaluated in a simulated online environment built from JD.com transaction records [9], with independent daily instances spanning the 06:00–18:00 peak operational window. The planning horizon \(\mathcal{T}\) corresponds to this window, and second-stage demand \(\tilde{\mathbf{D}}^t\) aggregates remaining intra-day arrivals from the demand forecaster’s 24-hour rolling output. Inventory and queue states are reinitialized each day from historical demand parameters with no intra-day replenishment. The primary metric is cumulative realized objective value, comprising base shipping cost, delivery-time penalties, multi-unit consolidation discounts, and stockout penalties. Late delivery percentage and cumulative lateness are reported as service-level metrics.

Because the public JD.com release is dominated by a single carrier and lacks exact facility coordinates, the dataset is augmented with a synthetic multi-carrier-service layer and synthetic geocoordinates, calibrated against a separate proprietary multi-carrier dataset. Preprocessing, calibration, and instance construction are detailed in Appendices 9 and 10. The first three weeks supply training data and the final week is held out for online evaluation. A counterfactual simulator generates calibrated delivery-time realizations for unobserved fulfillment decisions [7], and all reported metrics average 50 independent simulation runs.

Evaluation proceeds order by order in rolling-horizon fashion, with the inventory state updated after each decision. The C-SAA solver uses Gurobi 13.0.1 with four threads per instance, \(S=10\) candidate solutions at \(N_1=50\) scenarios each, common-set evaluation on \(N_2=500\) scenarios, and a 0.01% relative MIP gap. Experiments ran on a shared HPC cluster (RHEL 9.6, dual Intel Xeon Gold 6226 at 2.7 GHz, NVIDIA RTX 6000 GPUs).

5.1 Comparative Analysis↩︎

The proxy and the online C-SAA solver are compared against five baselines covering heuristic, point-forecast, and benchmark fulfillment approaches. The first three baselines are extended to the sequential setting from [7].

• Greedy: minimizes the item-specific base shipping cost \(c^{\mathrm{ship}}_{i,j,k}\) for each product \(i\) and ignores stochastic delivery-time deviations. It mirrors the common practice of always sourcing from the cheapest feasible option.

• Predict-Then-Optimize (PTO): substitutes \(\mathbb{P}(\cdot \mid \mathbf{X}^t)\) with its conditional mean and solves the deterministic equivalent. It is a forecast-driven policy that acts on expected values and ignores the spread of uncertainty.

• Empirical-SAA: samples \(\mathbf{\Omega}^t\) uniformly from the historical unconditional empirical distributions, ignoring contextual signals. It accounts for uncertainty but not for the context of the specific order.

• Deterministic LP (DTLP) [30]: solves a linear program on a demand forecast to obtain dual prices that value the future opportunity cost of consuming inventory, then fulfills each arriving order against those prices.

• Primal-Dual [39]: maintains dual prices on inventory and updates them online as orders arrive, without relying on a demand forecast.

Table ¿tbl:tab:sim95summary95test? reports aggregate out-of-sample performance over the test week: the total realized objective (base shipping cost, delivery-time penalties, consolidation discounts, and stockout penalties), the average late-delivery rate and cumulative lateness, the per-order runtime, and the \(q_{50}\)–\(q_{95}\) quantiles of each. In the calibrated simulation environment, the proxy attains a total objective of \(829{,}805\), improving on the online finite-sample C-SAA reference (\(857{,}892\)) by 3.3%. The average late delivery rate falls from 5.71% to 4.93%, and the daily gap is stable. The proxy’s edge over its own labels reflects scenario-noise averaging. The C-SAA reference at \((S,N_1,N_2)=(10,50,500)\) is itself exposed to scenario generation error, finite-\(N_1\) sampling variance, and two-stage aggregation bias, while the proxy averages over the full corpus of these noisy labels with \(\mathcal{L}_{\mathrm{cost}}\) as a structural regularizer. The proxy improves on C-SAA in aggregate but can lag where C-SAA’s finite sample lands near the true optimum, as on the largest tier in Table ¿tbl:tab:runtime95gap95size95test?. A higher-budget C-SAA may narrow the aggregate gap. Deterministic and point-prediction baselines lag substantially: Greedy reaches \(928{,}970\) with a 9.80% late rate, and PTO and Empirical-SAA both exceed \(1{,}000{,}000\) in objective and 11% in late rate, confirming the value of scenario-aware contextual information [7].

The proxy executes each decision in 0.006 seconds, an acceleration above \(2{,}800\times\) over C-SAA’s 17.084 seconds. It is also faster than the simpler Greedy (0.104s) and PTO (0.223s) baselines, fitting within the sub-second latency budgets of high-throughput fulfillment systems.

Relative to DTLP, the strongest classical baseline on objective, the proxy lowers total cost by \(10.7\%\) (\(829{,}805\) versus \(928{,}889\)) and roughly halves the late-delivery rate (\(4.93\%\) versus \(9.77\%\)); the cost gap widens to \(48.7\%\) against the primal-dual algorithm. Both run within the same sub-second budget as the proxy, yet neither matches it on cost or service. The primal-dual algorithm is the fastest method (\(0.004\)s) but the most expensive (objective \(1{,}618{,}102\)): it updates inventory dual prices online without a demand forecast, so they stay far from what the forecast-based and scenario-aware policies achieve. DTLP performs about as well as greedy on both cost and service (objective \(928{,}889\); late rate \(9.77\%\)) and runs in \(0.101\)s. Both price inventory through dual variables but ignore the contextual delivery-time uncertainty that the proxy and C-SAA capture through scenario forecasts, so neither avoids the options likely to ship late, and both leave late rates well above the proxy’s \(4.93\%\). The proxy outperforms both on every cost and service metric at comparable per-order speed.

width=,center

Figure 6 breaks the realized cost down by date, summing over the orders served each day and averaging across the 50 replications. The proxy stays at or below every baseline on each evaluation date and keeps its margin over the baseline policies across the week.

5.2 Computational Scalability↩︎

This Section highlights the computational performance of the proposed proxy against the baseline C-SAA approach by examining its scalability in terms of number of scenarios and size of the instances considered.

5.2.1 Scenario Size Scalability↩︎

Figure 7 reports solution quality and latency as the scenario sample size \(N_1\) varies. For the C-SAA solver, total objective value decreases monotonically with \(N_1\), but average runtime grows roughly linearly to about 30 seconds at \(N_1=90\). The proxy, by contrast, attains its minimum objective at the training value \(N_1=50\) and degrades little for larger \(N_1\), with per-decision latency remaining near 0.006 seconds across the tested \(N_1\) range.

5.2.2 Instance Size Scalability↩︎

Table ¿tbl:tab:runtime95gap95size95test? stratifies the test evaluation set by combinatorial complexity, defined as the product of the order quantity and the total number of eligible fulfillment options. As instance size grows, the C-SAA solver’s runtime increases from roughly 15 to 27 seconds, while the proxy’s per-decision time stays near 0.006 seconds across the tested tiers. On the largest instances the C-SAA point estimate is marginally lower than the proxy’s, but the 95% confidence intervals overlap, so the two are statistically indistinguishable on that tier. The proxy maintains a lower late delivery rate across all tiers.

width=,center

5.3 Proxy Ablations↩︎

An ablation study quantifies the contribution of individual framework components across the neural architecture, loss formulation, and inference strategy. Aggregate impacts on objective cost and service reliability are summarized in Tables ¿tbl:tab:proxy95architecture95ablation95summary?, ¿tbl:tab:proxy95loss95ablation95summary?, and ¿tbl:tab:proxy95inference95strategy95ablation95summary?. These ablations also stand in for natural learning-based baselines: no_cost_loss (imitation-heavy, drops \(\mathcal{L}_{\mathrm{cost}}\)), no_selection_loss (cost-only, drops solver-label supervision), and no_scenario_module (no-scenario proxy, same decoder). All are inferior to the full model on objective cost.

5.3.0.1 Architecture and Features:

Table ¿tbl:tab:proxy95architecture95ablation95summary? confirms the value of the hierarchical, multimodal design. Removing the DC module (no_dc_module) strips location-specific features and inventory state, costing 11.7% in objective and pushing the late rate to 5.45%. Flattening the heads (single_tower) costs 11.2% in objective with a comparable late rate, indicating that joint heads sacrifice cost economies for speed. Removing the scenario module (no_scenario_module) forces the proxy onto deterministic features and produces a consistent cost increase.

width=,center

5.3.0.2 Loss Formulation:

Table ¿tbl:tab:proxy95loss95ablation95summary? isolates the contribution of each composite-loss term. Removing the selection loss (no_selection_loss) costs 3.7% in objective, showing that direct cost-based training alone, without C-SAA label supervision, is insufficient. Removing the cost-alignment term \(\mathcal{L}_{\mathrm{cost}}\) costs 3.2% in objective and pushes the late rate to 5.63%, isolating its role in smoothing expected scenario costs and consistent with the proxy’s improvement over the finite-sample C-SAA solver. Removing the constraint loss \(\mathcal{L}_{\mathrm{const}}\) costs 1.8% in objective with comparable service-level degradation, confirming its role in shaping operationally feasible primary-DC commitments before decoding.

width=,center

5.3.0.3 Inference Strategy:

Table ¿tbl:tab:proxy95inference95strategy95ablation95summary? compares the proposed Inv-Weighted decoder against four alternative decoders specified in Appendix 13.4. Inv-Weighted and Greedy-Prob are within 0.1% of each other in objective, but Greedy-Prob lacks the explicit bias toward fully feasible first choices that Inv-Weighted’s inventory-adequacy ratio \(r_j^t\) provides. Hard-feasibility filters (Top-K-Feasible-DC, Feasible-DC-First) lose ground because they discard high-probability DCs that fall marginally short of full coverage. The joint-space search (Top-K-Feasible-Joint) breaks the factorized hierarchy and produces a 30.1% cost explosion with a late rate above 11%.

width=,center

5.4 Managerial Insights↩︎

The computational study suggests four practical implications for large omnichannel fulfillment networks. First, the main value of the proposed approach is not just faster execution. By shifting solver compute offline, stochastic fulfillment can be deployed in real time: offline solver-based solutions supervise a learned policy that then runs online in milliseconds. Second, scenario-aware methods materially outperformed point-estimate baselines on realized cost, so scenario generation is a primary modeling lever rather than a peripheral forecasting step. Third, the inference-strategy ablation shows that how predictive outputs are translated into feasible actions has a measurable effect on realized cost and service, so inventory-aware feasible inference belongs in the policy rather than in post-processing. Finally, the results support a hybrid operating model. The learned proxy serves high-throughput real-time traffic where latency precludes online stochastic optimization, and online C-SAA remains useful for offline policy generation and selected high-stakes orders.

6 Conclusion↩︎

This paper formulated online omnichannel order fulfillment as a two-stage contextual stochastic program and addressed the latency of its mixed-integer C-SAA reformulation through a scenario-embedded hierarchical optimization proxy paired with an inventory-weighted feasibility decoder. The network maps context and scenarios to distribution-center and carrier-service logits in a single forward pass, and the decoder turns them into a feasible per-order plan. The composite training loss augments solver-label imitation with a constraint-violation penalty and a self-supervised cost-alignment term, so the trained proxy can match or improve upon its finite-sample C-SAA reference under the same scenario forecaster. Theoretical analysis establishes that the decoder returns a feasible decision in \(O(|\mathcal{J}|\log|\mathcal{J}|+|\mathcal{J}||\mathcal{K}|)\) time per product line, that its expected cost upper-bounds the two-stage stochastic optimum, and that its out-of-sample value admits unbiased, consistent sample-mean and sample-variance estimators.

The study on a simulator built from JD.com transactional records evaluated the proxy against online C-SAA and heuristic baselines under a common fitted forecaster. The proxy reduced per-order decision latency by roughly \(2{,}800\times\) relative to online finite-sample C-SAA and ran in constant time as scenario density grew, whereas C-SAA latency scaled linearly. On realized fulfillment cost the proxy improved over the same reference by \(3.3\%\), with the strongest gains on small-to-medium instances and statistical parity on the largest tier.

Several limitations point to natural extensions. The per-product-line decomposition keeps inference fast but captures the multi-line consolidation discount only indirectly through the jointly optimized C-SAA labels. An end-to-end architecture would broaden applicability to multi-item baskets. The second-stage recourse aggregates intra-day demand into a single approximation and would benefit from finer multi-period coupling. Methodologically, extending the uncertainty quantification of the demand and delivery-time forecasters, for instance through conformal prediction or distributionally robust forecasting, would supply coverage guarantees on the sampled scenarios under distribution shift.

Acknowledgments↩︎

This research was supported by the NSF AI Institute for Advances in Optimization (Award 2112533).

7 Nomenclature↩︎

Table 1 summarizes the notation, sets, and mathematical operators employed throughout this paper.

8 C-SAA Evaluation Procedure↩︎

Algorithm 8 summarizes the resulting candidate-generation and evaluation procedure, where SolveSAA\((\Omega)\) denotes solving (2) with scenario set \(\Omega\) and returns the first-stage solution and the resulting optimal objective function value.

9 Fulfillment Instance Setup↩︎

This section details the construction of the daily fulfillment instances, including the generation of feasible options, state initialization, and the composite realized cost structure used in the simulation environment.

9.1 Simulation Scope and State Initialization↩︎

The simulation environment models independent daily instances rather than a continuous multi-day rolling horizon. At the start of each peak operational window, the system state is freshly initialized, meaning no inventory or queue backlogs carry over between simulation dates.

• Inventory Initialization: Initial on-hand inventory is computed separately for each simulation date using historical demand strictly preceding that date. For a given facility and item, the target stock level follows the standard safety-stock formulation \(\text{target} = \mu_{\text{daily}} + z \cdot \sigma_{\text{daily}}\), where \(\mu_{\text{daily}}\) and \(\sigma_{\text{daily}}\) are the mean and standard deviation of historical daily demand. For local distribution centers, the \(z\)-factor corresponds to a 50% cycle service level, and items are additionally filtered through a deterministic 75% stocking probability hash. Central warehouses aggregate demand regionally and use a higher 80% cycle service level. The simulator does not model intra-day replenishment; inventory strictly depletes as sequential fulfillment decisions allocate units.

• Queue State Warm-starting: Before processing the first order of the day, the simulator warm-starts the operational queue (e.g., waiting orders, waiting units, and pending shipments) using a historical snapshot from the preceding hours. This ensures the simulator begins from a realistic, non-empty operational state. Subsequent processing delays are estimated dynamically based on base unit processing times and queue backlogs, capped at a maximum of 720 minutes.

9.2 Feasible Options and Base Costs↩︎

For each arriving order, feasible fulfillment options are restricted to distribution centers and carrier services capable of delivering to the customer’s geographic state. The shipping distance from the origin DC to the synthetic customer location is computed using a vectorized Haversine formula based on hub-based pseudo-coordinates. The deterministic base shipping cost for an eligible (DC, carrier) pair is \(\text{Base Cost} = (\text{Distance in km}) \cdot \alpha_{\text{carrier}} + \text{Fixed DC Cost}\), where \(\alpha_{\text{carrier}}\) is the externally calibrated carrier distance coefficient and the fixed facility operational cost is 2.0 for local DCs and 4.0 for central warehouses. The remaining cost-function parameters of the main-paper definition take the values \(\beta = 0.5\) (consolidation discount), \(\gamma^+ = 40.0\) and \(\gamma^- = 0.2\) (per-unit-per-day late and early penalties), and \(\rho = 200.0\) (per-unit lost-sales penalty).

9.3 Delivery Time Simulator Construction↩︎

Because realized delivery-time deviations are only observed for the location–carrier options selected in the historical system, counterfactual delivery outcomes are unavailable for alternative fulfillment decisions. Following [7], a separate counterfactual simulator is therefore constructed to support online policy evaluation. The simulator uses the same multi-quantile MLP (MQMLP) class and feature pipeline as the delivery-time forecaster (defined in Appendix 13), but is invoked at a different point in the loop and on disjoint data: the simulator is fit on the historical realized-delivery records and is queried only after a policy decision is made (to score that decision), whereas the forecaster’s role at decision time is to produce the scenario set \(\mathbf{\Omega}^t_{\mathrm{time}}\) that policies see as input. The simulator therefore acts as a shared counterfactual evaluator: every policy compared in Section 5.3 consumes the same scenario inputs from the forecaster and is scored by the same simulator-sampled realizations, so the comparison among policies is consistent within the simulator, while not eliminating action-dependent bias of the simulator itself; calibration diagnostics by carrier and policy-relevant action region are an item for future work. For a realized context \(\mathbf{X}^t\) and a selected location–carrier pair \((j,k)\), the simulator produces conditional quantile estimates that define an approximate distribution for the delivery-time deviation, from which a realization is sampled to evaluate the chosen action. Because the simulator is a fitted predictive model rather than the true delivery process, it provides an approximate counterfactual evaluation environment for comparing policies.

10 Data Preprocessing and Carrier-Service Augmentation↩︎

10.1 Raw Data Cleaning↩︎

The foundational dataset is JD.com e-commerce transaction records [9]. Several preprocessing filters were applied to clean the inputs for the simulation and forecasting models. Entries with missing delivery promise days and duplicate (order_ID, sku_ID) pairs were removed. Records with missing brand_ID, single-item gift orders, and multi-package orders (identifiable by more than one unique package_ID per order) were excluded to focus the evaluation on standard multi-item consolidation dynamics.

Timestamps were parsed to compute the end-to-end delivery time (from order_time to arr_time) and the DC-to-DC travel time (from ship_out_time to arr_station_time). Records exhibiting negative durations due to logging errors were discarded. Delivery hours were converted to discrete days via \(\lceil \text{hours}/24 \rceil\) (counting same-day delivery as 1 day). The delivery deviation was explicitly defined as delivery_time_days minus promise_delivery_days. Finally, to remove extreme outliers, records with a delivery_time_days exceeding 5 days or a deviation falling outside the \([-5, 5]\) day window were filtered out.

10.2 Carrier-Service Augmentation↩︎

The baseline JD.com dataset is predominantly serviced by a single internal carrier and lacks precise geographic coordinates for its fulfillment facilities. To train and evaluate the hierarchical proxy and carrier-specific simulators, a calibrated transformation layer is applied to the cleaned data to generate a synthetic, carrier-aware spatial dataset. The two augmentation parameter families described below (the per-distance-bin carrier-eligibility shares and the per-carrier delivery-time multipliers) are calibrated from a separate proprietary multi-carrier fulfillment dataset rather than chosen ad hoc, so that the augmented benchmark preserves the joint distribution of distance, carrier, and delivery-time variability observed in real operations even though the underlying transaction stream remains the public JD.com release. Aggregate calibration footprint (no proprietary records released): 17 carrier services, 5 origin-destination distance bins, an eligibility share per (distance bin, carrier) pair, and a median delivery-time ratio per (carrier, distance bin) pair (so \(17 \times 5 = 85\) shares and \(85\) ratios).

The augmentation proceeds sequentially. First, distribution centers are assigned synthetic coordinates based on regional metadata. Synthetic customer locations are similarly generated by mapping the original destination facility’s three-digit (ZIP3) postal codes to consistent five-digit (ZIP5) codes, enabling precise geocoding into latitude and longitude coordinates. Second, the shipping distance from the origin distribution center to the synthetic customer location is computed utilizing a vectorized Haversine formula. Third, each order is categorized into a discrete distance bin. Fourth, a carrier-service identifier is sampled from the proprietary-calibrated eligibility pool for the corresponding distance bin. Finally, a carrier-level delivery-time perturbation is applied: the original observed delivery time is multiplied by the proprietary-calibrated median delivery-time ratio for the sampled (carrier, distance) pair.

This process yields a deterministic hybrid dataset (controlled via a fixed random seed) that preserves the anchoring of the original transactional observations while introducing spatial coordinates and re-scaling delivery times to reflect the realistic stochastic variations induced by different carrier selections.

11 Forecasting Model Feature Details↩︎

11.1 Demand Forecasting Features↩︎

The demand pipeline constructs rolling-window tensors using a historical lookback of 36 hours and a forecast horizon of 24 hours. The historical covariate tensor consists of 23 features per time step: 4 cyclical calendar variables (sine and cosine transformations of hour and day of the week) and 19 aggregated order covariates (including SKU counts, bundle and coupon discount indicators, total quantity, and customer demographic encodings). The historical target tensor contains the realized demand quantity per time step. The future covariate tensor is strictly limited to the 4 cyclical calendar variables over the prediction horizon. Missing numerical values are imputed with a constant (\(-1\)), and discrete covariates are factorized into integer indices. Baseline tree-based models flatten these temporal tensors into summary statistics (mean and standard deviation) or use the fully flattened temporal sequences.

11.2 Delivery-Time Forecasting Features↩︎

The delivery-time models use a shared feature set constructed from preprocessed order data. The target variable is the realized delivery duration in days. Numerical features include temporal indicators, order complexity metrics, and DC-operations snapshots (specifically, the volume of shipped orders and SKUs in the preceding two hours, alongside current waiting queues). Categorical features are limited to the origin and destination facility identifiers, which are label-encoded and embedded in the deep learning pipelines.

12 Forecasting Model Performance↩︎

Quantile regression models were trained and evaluated for both delivery-time and demand forecasting. The candidate model families are XGBoost’s XGBRegressor (XGBoost), scikit-learn’s QuantileRegressor (QR), quantile-forest’s RandomForestQuantileRegressor (QRF), an MLP-based MQMLP for delivery time, and a Multi-Quantile Recurrent Neural Network (MQRNN) for demand. All models predict 19 quantiles (0.05 to 0.95 in increments of 0.05); predictions are rounded to integer values before scoring against the ground truth. Reported metrics are pinball loss and Continuous Ranked Probability Score (CRPS), averaged across the 17 carriers and 19 quantiles for delivery time (Table 2) and across 24 future horizons and 19 quantiles for demand (Table 3). For demand, MQRNN processes full temporal sequences through its recurrent architecture, whereas XGBoost, QR, and QRF use summary statistics (mean and standard deviation) of the temporal features.

13 Neural Network Details↩︎

13.1 Delivery Time Forecasting Model (MQMLP)↩︎

The multi-quantile MLP (MQMLP) delivery time model is a feed-forward MLP that predicts multiple delivery-time quantiles. Inputs consist of numerical features and origin/destination DC IDs; the DC IDs are mapped through learned embeddings and concatenated to the numerical features before entering the network. The output layer returns one value per requested quantile.

13.1.0.1 Non-crossing quantiles.

To avoid quantile crossing, the head is built around the median (\(0.5\)) as an anchor. The model predicts (i) a median value and (ii) nonnegative “gaps” for higher and lower quantiles via a softplus transform. Higher quantiles are the median plus cumulative upward gaps; lower quantiles are the median minus cumulative downward gaps. This construction produces nondecreasing quantiles by design.

13.1.0.2 Bounded support.

Predictions are mapped to the fixed interval \([t_{\min}, t_{\max}] = [1, 5]\) via \(\hat{q} = t_{\min} + (t_{\max}-t_{\min}) \cdot \sigma(\text{raw})\) with a monotone sigmoid \(\sigma\).

13.1.0.3 Training objective.

Training uses the pinball (quantile) loss summed over all predicted quantiles. Early stopping is applied on a validation split.

13.1.0.4 Hyperparameters.

Because the models are trained separately for each carrier, the selected hyperparameter configuration may vary by carrier. Table 4 reports the hyperparameter search space. For each model, 500 unique hyperparameter combinations are sampled, with learning rate and weight decay drawn on a log scale. For each carrier, the final hyperparameter setting is selected as the configuration attaining the lowest validation loss.

13.2 Demand Forecasting (MQRNN)↩︎

The MQRNN demand model follows an encoder–decoder design. Two input encoders map (i) the history window (calendar features, order-statistics, and SKU/brand IDs via learned embeddings) and (ii) the prediction window (calendar features and the same IDs) into compact contexts. An LSTM then processes the history context concatenated with the observed target series. A two-level decoder produces multi-horizon quantiles: a global MLP builds horizon-agnostic and horizon-specific contexts, and a local head refines each horizon.

13.2.0.1 Non-crossing & non-negative quantiles.

The quantile head is centered at the median as an anchor. It predicts nonnegative “gaps” for higher and lower quantiles using a softplus transform and accumulates them outwards from the median, which ensures ordered (non-crossing) quantiles. A final ReLU enforces non-negative forecasts.

13.2.0.2 Training objective.

Training uses the pinball (quantile) loss summed over all horizons and quantiles, with early stopping on a validation split.

13.2.0.3 Hyperparameters.

Table 5 summarizes the hyperparameter search space and the selected best configuration for the MQRNN model. A total of 50 unique hyperparameter combinations are sampled, with learning rate and weight decay drawn on a log scale, consistent with the tuning procedure used for the delivery model. Training is run for up to 200 epochs, and early stopping is applied using the patience value selected from the search space.

13.3 The Optimization Proxy↩︎

This appendix provides concrete dimensions and feature inventories for the symbolic specification of \(f_\theta\) given in Section 4.2.1 (Figure 3, with notation defined in Table 1). Table 6 lists the constituent fields of each input embedding source: categorical fields use the lookup widths shown in the table and continuous fields are min–max scaled prior to entering the network. Table 7 reports the per-module input source, hidden-layer structure, and output dimension. Table 8 reports the hyperparameter search space and the selected production configuration; the hidden width \(h\), dropout probability, and number of head hidden layers are tuned, while the structural choices in Table 7 are fixed.

13.4 Decoder Variants for the Inference Strategy Ablation↩︎

The four alternatives to Algorithm 1 compared in Table 5 consume the same logits \((\boldsymbol{\ell}^{\mathrm{dc}}, \boldsymbol{\ell}^{\mathrm{carrier}})\) and return a feasible decision in \(\mathcal{Z}^t\); they differ only in the DC-ranking score and eligibility filter, with the per-DC carrier choice \(k^*(j) = \arg\max_{k \in \mathcal{K}_j^t} p_{j,k}^{\mathrm{carrier}}\) shared and \(K\) denoting the top-\(K\) rank cutoff.

• Inv-Weighted (proposed, Algorithm 1): \(s_j^t = p_j^{\mathrm{dc}}\, r_j^t\) with \(r_j^t = \min\{\mathrm{Inv}_j^t/q^t,\, 1\}\); greedy fill on \(\mathcal{J}^{\mathrm{elig}} = \{j : \mathrm{Inv}_j^t > 0,\, \mathcal{K}_j^t \neq \emptyset\}\) in decreasing \(s_j^t\).

• Greedy-Prob: identical to Inv-Weighted but with \(s_j^t = p_j^{\mathrm{dc}}\) (the multiplier \(r_j^t\) is dropped).

• Top-K-Feasible-DC: restricts the eligible set to the top-\(K\) DCs by \(p_j^{\mathrm{dc}}\) that also satisfy \(\mathrm{Inv}_j^t \ge q^t\), then greedy-fills on this restricted set; relaxes to \(\mathcal{J}^{\mathrm{elig}}\) if no DC qualifies.

• Feasible-DC-First: first selects the highest-\(p_j^{\mathrm{dc}}\) DC satisfying \(\mathrm{Inv}_j^t \ge q^t\); falls back to the Inv-Weighted fill if none meets full coverage.

• Top-K-Feasible-Joint: ranks \((j,k) \in \mathcal{J} \times \mathcal{K}\) pairs by \(p_j^{\mathrm{dc}}\, p_{j,k}^{\mathrm{carrier}}\), retains the top \(K\), and assigns greedily under \(\mathrm{Inv}_j^t > 0\) and \(k \in \mathcal{K}_j^t\); multiplying \(p_j^{\mathrm{dc}}\) with \(p_{j,k}^{\mathrm{carrier}} = p(k \mid j)\) as if independent breaks the hierarchical factorization \(f_\theta\) is trained against.

14 Reproducibility↩︎

The computational pipeline runs in six stages: preprocessing the public JD.com transaction records; the carrier-service augmentation of Appendix 10; training the demand and delivery-time scenario generators; offline C-SAA label generation; proxy training; and online simulation. Fixed random seeds control augmentation, scenario generation, proxy training, and simulation replications; C-SAA scenario streams use order-indexed seeds derived from the global seed. After preprocessing, the study contains \(237{,}407\) orders, \(1{,}488\) SKUs, \(55\) distribution centers, and \(17\) carrier services. The proxy is trained on \(31{,}907\) peak-hour order-SKU examples with \(500\) scenarios per example, and the online test set contains \(38{,}416\) peak-hour orders over March 26–31, 2018. Reported runtimes measure policy execution time and exclude scenario generation and post-decision simulator sampling.

References↩︎

---

1. Corresponding author: joe.ye@gatech.edu↩︎