1 Introduction↩︎

The joint replenishment problem (JRP) concerns the management of multi-item inventory systems, in which a fixed replenishment cost is incurred jointly by items procured together. A JRP is called stochastic (SJRP) when demands are random. A classical problem in operations management, the JRP has been extensively studied over several decades with wide-ranging applications; see, for example, the literature surveys in [1], [2], [3] and [4]. In particular, the SJRP is relevant to sectors involving large-scale operations with highly variable demand and vast product assortments, such as modern retail, where joint replenishment can amortize fixed costs such as transportation, and yield substantial economies of scale.

While coordinated replenishment offers attractive cost-saving opportunities, effective deployment in practice is challenging, especially when the number of distinct items, i.e., the problem dimension, is large. A fundamental reason is that item-level inventory decisions become interdependent under nonzero fixed replenishment costs. This interdependence creates two difficulties: (i) the computational effort required to solve the associated Markov decision process (MDP) scales with the sizes of the inventory state and decision spaces, both of which grow exponentially with dimension, leading to the “curse of dimensionality”; and (ii) fixed costs induce large multi-item orders, producing jumps in the inventory state and amplifying dimensionality effects. As a result, there has been limited progress in the literature on obtaining optimal replenishment policies in the multidimensional setting; see, for example, Section 6 of [5] for a review of related literature. Instead, a substantial body of literature has focused on designing effective heuristics for the SJRP; see, for example, [3]. These heuristics typically propose intuitively appealing classes of parameterized inventory control policies, reducing the control problem to the optimization of the policy parameters.

In this paper, we revisit the classical SJRP [1], which is formulated in discrete time, and develop an effective computational method for the high-dimensional setting. Toward this end, we consider a closely related continuous-time impulse control problem, in which instantaneous jumps in the system state can be effected at a fixed cost, and whose solution naturally leads to implementable policies for the original inventory control problem. Traditional numerical approaches to impulse control problems use grid-based methods, such as finite difference and finite element schemes, to solve the associated Hamilton-Jacobi-Bellman (HJB) partial differential equations (PDEs), but they also suffer from the curse of dimensionality. Compared to the discrete-time SJRP, the continuous-time impulse control formulation has a key advantage: It admits an equivalent probabilistic representation of the associated value function [6]. This representation allows us to employ the deep backward stochastic differential equation (BSDE) method developed in the seminal work of [7] for tackling high-dimensional problems. Since [7], a large literature has used deep neural networks to solve high-dimensional PDEs following the deep BSDE framework. We build on and contribute to this literature by developing a computational method for solving the high-dimensional impulse control problem. More importantly, we leverage that solution to propose effective control policies for the classical SJRP.

To put our method in perspective, we briefly review the deep BSDE framework, highlight our main challenges, and discuss the innovations we introduce to address them. Deep BSDE is a simulation-based computational approach for solving high-dimensional PDEs using deep neural networks. It starts from a probabilistic (BSDE) representation of the target function, typically a stochastic equality satisfied almost surely (a.s.) along sample paths of a reference process, which is a simulatable random process chosen by the algorithm designer. Based on this representation, one constructs a loss function and trains neural network approximations for the target function (and its gradient) by minimizing the loss, using samples generated by simulating the reference process. When the representation is a stochastic equality, the loss can be taken as the standard mean squared error between the left- and right-hand sides, with the function and gradient replaced by their neural network approximations, respectively. Under appropriate choices of the loss and the reference process, the success of the deep BSDE method then hinges primarily on extensive engineering efforts in hyperparameter tuning and neural network training.

For our impulse control problem, two immediate challenges arise in designing the loss and the reference process within the deep BSDE framework. First, while an equivalent probabilistic representation exists, it takes the form of a family of almost-sure stochastic inequalities rather than an equality [6]. To address this more complex case, we design a novel, nonstandard loss that dualizes inequality constraints by incorporating a loss term for violations of the stochastic inequality and a tunable penalty parameter. The penalty parameter introduces considerable complications for neural network training, as its tuning schedule must be carefully coupled with the learning rate schedule to ensure training stability. Second, since the system state can be instantaneously displaced, deriving an effective control policy from the trained neural networks demands more than accurate local approximations of the value function and its gradient, as is common in prior literature (see, e.g., [7]–[9]); rather, the value function and its gradient must be well approximated across an expansive region of the state space. This, in turn, requires the reference process to explore the state space broadly. Indiscriminate exploration, however, yields irrelevant samples that hinder training. Thus, the reference process must strike an effective balance between exploration and exploitation; although this trade-off has been considered in prior work, it is substantially more prominent and nuanced in the context of impulse control. Toward this end, we design several sophisticated reference processes with carefully designed state-dependent jumps to ensure that relevant samples are generated for efficient training.

As mentioned earlier, our neural network solution to the impulse control problem naturally yields implementable inventory control policies for the original (discrete-time) formulation of the SJRP. We demonstrate the effectiveness of these policies through a series of test problems in the original SJRP setting. For low-dimensional test problems, we compare our policies primarily to the optimal ones derived by solving the MDPs associated with the original inventory control problem. We also compare our policy with the best available heuristics. Across settings, the performance of the proposed policy remains consistently within 1% of the optimal policy. It also matches or outperforms the best-performing heuristic. Moreover, the policy itself closely tracks the optimal policy (see Figures 4 and 5). For the high-dimensional test problems, however, solving the underlying MDP is computationally intractable. Thus, we benchmark our approach against the best available heuristic policies. For the high-dimensional test problems considered, our policy matches or beats the best performing heuristic. The rest of the paper is structured as follows. Section 2 provides an overview of related literature. Section 3 presents our model of the SJRP and its impulse control approximation. Section 4 introduces the HJB equation, and the aforementioned stochastic identity, Eq. 13 . Section 5 describes our algorithm. Section 6 presents our test problems and the benchmark policies. Computational results are described in Section 7. We conclude in Section 8 with some discussion on possible future research directions. Appendices 9–12 contain further technical details, which include a formal derivation of the key stochastic identity, how the hyperparameters of our neural networks are tuned, a complete description of the benchmark policies, and the validation of our computational approach in a one-dimensional inventory model with fixed cost.

2 Literature Review↩︎

We review the following three relevant streams of literature: (i) Joint replenishment problems, (ii) impulse control problems, and (iii) deep learning methods for solving high-dimensional PDEs.

Joint replenishment problems. The JRP is a classical problem in inventory management that has received considerable research attention over the past few decades. Broadly speaking, the JRP literature can be divided into problems with deterministic demand and those with stochastic demand [1]. The deterministic JRP can be viewed as the most basic multi-item extension of the classical economic order quantity (EOQ) model [10]. It assumes constant demand rates over an infinite planning horizon, with replenishment occurring at regular intervals. The objective is to determine the optimal frequency for replenishing the different items (see, e.g., [11]–[13]). One class of effective policies is characterized by the power-of-two rule under which orders are placed at power-of-two multiples of a base period. [14]–[16] have shown that the best power-of-two policies can perform within two to six percent of optimality. Much of the (deterministic) JRP literature focuses on the following cost structure: There is a group replenishment fixed cost, which is incurred every time an order is placed, and there are individual fixed costs for each item. By contrast, the generalized JRP (GJRP) allows the shared fixed costs to depend on the given collection of items included in the replenishment [17]. [18] prove the existence of optimal policies for the GJRP with capacity constraints using duality theory. [19] use approximate dynamic programming to determine near-optimal policies in this setting.

Our paper, however, concerns the SJRP; we now review this literature. For the single-item problem with procurement fixed costs, [20] introduces the notion of \(K\)-convexity to show that the optimal policy has the following simple structure: If the inventory level drops to or below a fixed level \(s\), an order is placed to raise it back to another fixed level \(S\). Efficient algorithms exist to find optimal parameters for such \((s,S)\) policies; see, for example, [21], [22]. In a single-item variant involving advance demand information, [23] give a geometric characterization of a convexity notion that generalizes \(K\)-convexity and simplifies the classical optimality proofs of \((s,S)\) policies and their state-dependent generalizations. Since Scarf’s seminal work, there have been several attempts to identify effective structural properties of optimal policies for multi-item inventory models. We refer the readers to [24] for a detailed review of these works, and discuss only the most relevant literature. [25] considers an infinite-horizon, \(d\)-item SJRP in discrete time, with integer-valued demand and the inventory process living in the state space \(\mathbb{Z}^d\). The demand in each time period is allowed to depend on the initial inventory levels in that period. The cost structure is identical to ours, i.e., there is a single fixed cost and linear, item-specific variable costs, whenever an order is placed. He considers both the total discounted cost and long-run average cost objectives and establishes the optimality of a so-called \((\sigma,S)\) policy under both objectives, where \(\sigma\) is a subset of \(\mathbb{Z}^d\), and \(S\) is a vector in \(\mathbb{Z}^d\). Under this policy, if the inventory state \(x\) is in \(\sigma\) and \(x\leq S\) component-wise, then one orders up to the levels \(S\), and does not order if \(x\notin \sigma\). Note that the policy is not specified when \(x \in \sigma\) and \(x \not \leq S\), but if the initial inventory state \(x\leq S\), then under a \((\sigma,S)\) policy, the inventory state always stays below \(S\) (component-wise), so a \((\sigma,S)\) policy is optimal in this case. Furthermore, for the long-run average cost objective, even if the initial state \(x\not \leq S\), an optimal policy can simply first wait for the inventory levels to deplete and drop below \(S\), and then follow a \((\sigma, S)\) policy. [26] proves the optimality of \((\sigma,S)\) policies under conditions that are different from those in [25], and characterizes the complement of \(\sigma\), a set where no order is placed, as being increasing under a certain partial ordering. His proof makes use of the notion of \((K,\eta)\)-quasiconvexity in \(\mathbb{R}\), a generalization of the concept of quasiconvexity used in [27] and in [28], for various single-item inventory problems. [29] recover results in [26] using generalized notions of \(K\)-convexity in \(\mathbb{R}^d\). Even though the optimal policy structure is known when the initial inventory state \(x\leq S\), determining the ordering region \(\sigma\) and the order-up-to level vector \(S\) is extremely challenging computationally, especially in high dimensions.

Due to the computational intractability, research on SJRP often focuses on classes of parameterized control policies that are intuitively appealing. Once attention is restricted to a particular class of control policies, the problem reduces to optimizing the policy parameters. This is often accomplished by numerical evaluations of the policy performances through simulation. See, for example, [30]–[36]. For a comprehensive review of the various heuristic control policies, we refer to [3].

Impulse control problems. One year after Scarf established the optimality of the \((s,S)\) policy in discrete time, [37] demonstrates that this policy is also optimal in continuous time with Poisson demand. To prove this, [37] uses what is known as uniformization. This approach restricts ordering decisions to occur only at the jumps of the Poissonian demand, effectively reducing the problem to its discrete-time counterpart. To address more general demand models such as diffusion or mixed-diffusion processes, impulse control problems offer a flexible framework that can effectively capture the lump-sum ordering nature of continuous-time inventory control problems with fixed costs. These stochastic control problems allow the system controller to intervene through jumps in the underlying state space [38]. Several authors have used the associated HJB equations to prove the optimality of \((s,S)\)–type policies for various one-dimensional inventory models; see, for example,
[39]–[43]. For applications of the one-dimensional impulse control problems in economics, we refer the reader to [44].

While numerous studies have investigated the continuous-time inventory model with fixed costs for single-item problems, only a handful of papers have addressed multi-item problems ([24]). A pioneering attempt to solve the multi-item setting is made by [45], who considers a two-item inventory system with constant, deterministic demand rates. There are individual fixed costs for ordering a single item, and a fixed cost of ordering both items, which is lower than the sum of individual fixed costs. The variable costs are assumed to be zero. Using the associated HJB equation, she establishes the optimality of a \((\gamma, \Gamma)\) policy, where \(\gamma\) represents the ordering boundary, and \(\Gamma\) denotes the order-up-to level boundary, similar to the \((\sigma, S)\) policy in the discrete-time setting. [46] also consider two-item inventory systems, and extend [45] in several ways. First, they allow variable costs, and establish the optimality of a \((\gamma, \Gamma)\) policy in the deterministic setting with constant demand rates. They provide characterizations of the boundaries \(\gamma\) and \(\Gamma\) in this case, facilitating their numerical constructions. Second, they establish various theoretical properties of a stochastic model with diffusion demand, such as the \((K, \eta)\)-quasiconvexity of the value function, and used the finite-difference scheme to obtain numerical solutions. Their numerical study suggests that the optimal policy in the stochastic case is also of the \((\gamma, \Gamma)\) type. Note that their finite-difference scheme is not computationally feasible in high-dimensional settings, i.e., when the number of distinct items is large.

Other common numerical methods for solving impulse control problems include iterated optimal stopping and policy iteration; see, e.g., [47], [48], and Chapter 12 in [49]. Unfortunately, these methods face the same computational tractability issues.

Deep learning methods for solving high-dimensional PDEs. Deep learning has gained significant traction as a method for solving high-dimensional PDEs since the seminal papers [7], [50]. Their method uses the probabilistic characterization of semilinear PDEs as backward stochastic differential equations (BSDEs); see, for example, [51]. The method outlined in [7] proposes to discretize the BSDE forward in time using a standard Euler scheme, and approximate the value function and its gradient with neural networks. It then learns the value function and its gradient by optimizing the neural network parameters through minimization of the expected squared distance between the known terminal value and the discretized BSDE’s terminal value. We refer to [52] for a detailed convergence study. With the advent of high-performance GPU computing, these deep-learning methods have been successfully applied to solve high-dimensional PDEs; see [53]–[55] for extensive literature reviews on the subject.

[9] and [8], [56] introduced these techniques to the operations management domain. The former studies a semilinear parabolic PDE linked to a drift control problem arising in dynamic call center scheduling (also see [57]), while the latter two develop computational methods for solving control problems in stochastic processing networks in the heavy-traffic regime. [58] solves a sequential vehicle routing problem motivated by an eviction enforcement application using a similar approach. [59] studies dynamic control of a make-to-order manufacturing system with throughput time constraints in high dimensions using deep learning. [60] solves a high-dimensional dynamic stochastic matching problem building on the deep BSDE method. The preceding literature predominantly addresses drift-rate control problems within the canonical deep BSDE setting, where the value function satisfies an almost-sure equality and the neural network training relies on a standard mean-squared loss. In contrast, while our impulse control problem is also tackled within the broad deep BSDE framework, it differs fundamentally: Our approach is based on the framework of BSDEs with constrained jumps [6], which provides a probabilistic representation of the value function of impulse control problems as the minimal solution to a constrained BSDE — a substantive methodological departure from the unconstrained BSDEs that underlie the canonical deep BSDE setting. Consequently, the stochastic identity satisfied by the value function involves an almost-sure inequality rather than an equality. This renders the standard mean-squared loss unsuitable, and necessitates a novel loss function penalizing violations of the inequality, deliberately crafted reference processes, and careful tuning of joint penalty-rate and learning-rate schedules; see also the discussion in Section 5. Relatedly, our deep-learning framework is inspired in part by the mathematical finance literature [61], [62], where deep learning is employed to identify optimal hedging strategies for pricing financial derivatives, problems often modeled by nonlinear PDEs. In particular, the stochastic identity (Eq. 13 ) on which our approach is based exploits a connection between impulse control problems and the so-called stochastic target problems, the latter of which have been studied extensively in the mathematical finance literature (see, e.g., [63, Ch. 7]). Finally, [64] studies an optimal switching problem — a problem class related to impulse control [65] — in energy markets, building on the framework of [66]. Unlike our setting, theirs involves a finite set of modes and exogenous, rather than controller-induced, jumps in the state. They further establish error bounds for their approach, whereas we do not pursue such bounds here; see also [67] and [52] for error analyses of other deep learning methods for solving PDEs.

3 Model↩︎

Consider a warehouse that stocks and sells \(d\) different items. The warehouse manager may place orders for any subset of these items. Define the order vector \(y=(y_1,\ldots,y_d)^\top\), where \(y_i\geq0\) denotes the number of units of item \(i\) ordered and the superscript \(^\top\) denotes transpose, and let \(c(y)\) denote the total ordering cost, where \[c(y) = \begin{cases} c_0 + \sum_{i=1}^d c_i y_i,\; &\text{if } \sum_{i=1}^d y_i > 0,\\ 0, &\text{ otherwise.} \\ \end{cases}\] That is, the warehouse manager incurs a fixed cost \(c_0\) for placing the order—for example, \(c_0\) may represent the fixed cost of using a shared resource associated with an order such as the transportation cost incurred from using a truck—and \(c_iy_i\) denotes the associated variable cost of ordering \(y_i\) units of item \(i\) for \(i=1,\ldots,d\).

We let \(x=(x_1,\ldots,x_d)^\top\) denote a generic system state, where \(x_i\) represents the inventory level of item \(i\). The warehouse manager reviews the system state periodically at times \(n=1,2,\ldots\). At the beginning of each period, upon observing the system state, but before observing demand, he makes an ordering decision; it is possible that the manager decides not to place an order. If an order is placed, it is received immediately, i.e., the delivery lead time is zero. Then, at the end of the time period, demand is realized, and any unmet demand is backlogged. In addition to ordering costs, there are linear holding and backlog costs. To be specific, a penalty cost \(p_i\) is incurred for each unit of backlogged demand of item \(i\), and a holding cost \(h_i\) is charged for each unit of inventory of item \(i\) carried over from one period to the next. Thus, if we let \[f_i(x_i)=\begin{cases} h_ix_i, \;& x_i\geq0, \\ -p_ix_i, \;& x_i<0, \end{cases} \label{eq:inv95cost}\tag{1}\] then we have that \[f(x) = \sum_{i=1}^d f_i(x_i)\] is the per-period state cost borne by the system manager when the system state is \(x\). All costs are incurred at the end of the period.

We now describe the system dynamics mathematically. First, let \(X(n)\) denote the system state (i.e., the vector of inventory levels for items \(1,\ldots,d\)) at the end of period \(n\). Second, let \(\xi(n)\) denote the \(d\)–dimensional demand vector in period \(n\). We assume that \(\{\xi(n)\}_{n \geq 1}\) are independent and identically distributed nonnegative random vectors with mean vector \(\mu \in \mathbb{R}^d\) and positive-definite covariance matrix \(\Sigma \in \mathbb{R}^{d \times d}\). To facilitate future analysis, we define the cumulative demand during the first \(n\) periods as \[\label{eq:cumdemand} D(n) = \sum_{k=1}^n \xi(k), \;n\geq1.\tag{2}\] Third, let \(y(n)\geq0\) denote the order vector in period \(n\). Then, the evolution of the system state can be expressed as \[\label{eq:inventorydynamics} X(n) = X(n-1) + y(n)- \xi(n),\;n\geq1,\tag{3}\] where \(X(0)\) denotes the initial system state.

Focusing on the infinite-horizon discounted cost criterion, the warehouse manager’s problem is to choose a non-anticipating control (policy), denoted by \(u=\{y(n), n\geq1\}\), consisting of order vectors for each period. The associated expected total discounted cost can be written as follows: \[\label{eq:MDPobjective1} {\tilde{J}(x, u)}=\E_x\Big[\sum_{n=1}^{\infty} \gamma^{n-1} \left(f(X(n)))+c(y(n))\right)\Big],\tag{4}\] where \(\gamma\in(0,1)\) is the discount factor, and \(\E_x[\cdot]\) denotes the expectation conditional on the initial state \(X(0)=x\).

For every initial state \(x\), we seek to find the optimal ordering policy that minimizes the total expected discounted cost 4 . That is, we want to determine \[\label{eq:MDPvaluefunction} {\tilde{V}(x)\mathrel{\vcenter{:}}= \inf_u \tilde{J}(x,u)},\tag{5}\] where the infimum is taken over the class of non-anticipating controls. As mentioned earlier, it is computationally prohibitive to solve 5 because the standard MDP-based techniques suffer from the curse of dimensionality. With this in mind, we next describe a continuous-time approximation of the foregoing problem.

Impulse Control Approximation. We now formally derive an impulse control approximation to the preceding inventory control problem. One can rewrite the stochastic dynamics of the inventory state process 3 as \[\label{eq:discrete95process} X(N) = x + \sum_{n=1}^N y(n)- D(N),\;N\geq1.\tag{6}\] The key simplification is to view the warehouse manager’s problem in continuous time and allow him to place orders at any point in time. With slight abuse of notation, we define a new continuous-time process \[\label{eq:continuous95process} X(t) = x + \sum_{\{j : \tau_j \leq t\}} y^{(j)} - D(t),\;t\geq 0,\tag{7}\] where the time \(t\) replaces \(N\) in Eq. 6 and \(y^{(j)}\) represents the \(j\)th order placed at time \(\tau_j\). In addition, we approximate the cumulative demand process \(D(\cdot)\) by a Gaussian process. To be specific, we assume \[\label{eq:browniandemand} D(t)= \mu t + \sigma B(t),\;t\geq0,\tag{8}\] where \(B(\cdot)\) is a \(d\)–dimensional standard Brownian motion, \(\mu \in \mathbb{R}^d\) is the drift vector, and \(\sigma \in \mathbb{R}^{d\times d}\) is chosen so that \(\sigma \sigma^{\top} = \Sigma\), where \(\Sigma\) is the covariance matrix of the per-period demand vector \(\xi(n)\). This Gaussian approximation of the original demand process is obtained by matching its first two moments. Such Gaussian approximations are ubiquitous in the inventory management literature, e.g., the newsvendor model and its extensions, because they enhance tractability. Despite its benefits, this approach introduces a potential approximation error. Therefore, once we derive the proposed policy under this approximation, we test its performance in the context of the discrete-time formulation using the original demand distribution. As the computational results in Section 7 demonstrate, our approach yields policies for the original discrete-time problem formulation that match or outperform the best available benchmark, confirming the practical effectiveness of this approximation.

The warehouse manager’s goal is to minimize the total expected discounted cost. More precisely, let \(0 < \tau_1 < \tau_2 < \ldots\) denote the sequence of strictly increasing, non-anticipating stopping times at which he places orders. Then, under the control policy \(u=\{(y^{(j)},\tau_j),j\geq1\}\), the system state evolves according to the \(d\)-dimensional stochastic process \[X^{u}(t) = x - \mu t - \sigma B(t) + \sum_{\{j : \tau_j \leq t\}} y^{(j)}, {\;t\geq 0.}\]

Letting the annual interest rate be \(r > 0\), the expected total discounted cost under policy \(u\) is given by \[{J(x, u)} = \E_x\Big[\int_0^\infty e^{-r t} f(X^u(t)) dt + \sum_j e^{-r \tau_j} c(y^{(j)})\Big].\] We then define the optimal value function as follows: \[\label{eq:impulsevaluefunction} V(x)=\inf_{u\in\mathcal{U}} {J(x,u)}\;\text{for each } x\in\mathbb{R}^d,\tag{9}\] where the infimum is taken over the set \(\mathcal{U}\) of non-anticipating policies.

4 The HJB Equation and an Equivalent Stochastic Identity↩︎

In this section, we describe an identity, Eq. 13 below, that is similar in spirit to the results used by [7] to justify their deep BSDE method for solving semilinear PDEs. That earlier work provided the motivation for our approach, though the specifics of our approach to solving the impulse control problem differ significantly.

The value function \(V\) we seek in 9 can be characterized as the solution to a particular kind of highly nonlinear PDE, namely, the HJB equation associated with the impulse control problem. Using the infinitesimal dynamic programming principle, one can show formally that the value function \(V\) should satisfy the following HJB equation: \[\label{eq:qvi} \left(\left(\mathcal{L}+r\right) V - f \right) (x) \vee \left( V(x) - \inf_{y \in \mathbb{R}_+^d} \left(V(x+y) + c(y)\right)\right) = 0, {\;x\in\mathbb{R}^d,}\tag{10}\] where \(\mathcal{L}:=-\frac{1}{2}\cdot \text{Tr}\left(\sigma \sigma^\top \text{Hess}(V) \right) + \mu^\top \nabla V\) is the diffusion operator; \(\nabla V\) and \(\text{Hess}(V)\) denote the gradient and Hessian of \(V\), respectively; and \(a \vee b\mathrel{\vcenter{:}}=\max\{a,b\}\) for scalars \(a\) and \(b\). To elaborate, the first term on the left-hand side of Eq. 10 , which we will refer to as the “no action” term, represents taking no action, which intuitively means that the warehouse manager opts to do nothing when the system state is \(x\). The second term, which we will refer to as the “intervention” term, corresponds to placing an order when the system state is \(x\), with the warehouse manager determining the order size vector \(y^*\) that minimizes future expected discounted costs.

The HJB equation 10 will not be used directly to learn the value function \(V\) in our approach. Instead, we rely on a related probabilistic characterization of \(V\), to be described shortly in Eq. 13 . To this end, we begin by specifying what we call a reference policy that is used to generate sample paths of the inventory process. That is, it facilitates our sampling of the paths living in the state space. We refer to the state process under the reference policy as the reference process, denoted by \(\tilde{X}(\cdot)\). Roughly speaking, we choose a reference policy so that the reference process tends to occupy the parts of the state space that we expect the optimal policy to visit frequently, so as to enable efficient and stable training of the neural networks. Note that the reference process plays the same role as the forward diffusion process in standard deep BSDE methods; see, e.g., Eq. (16) in [7].

As mentioned earlier, [25] characterized the optimal policy for the SJRP as a \((\sigma, S)\) policy. Also, for our approximating impulse control problem, a careful examination of the HJB equation 10 suggests that whenever an order is placed, the inventory vector is brought up to an order-up-to vector that is independent of the current inventory vector, provided the latter does not exceed the order-up-to vector (more will be said about this in Section 5). Thus, we build on a (randomized) order-up-to policy for choosing a reference policy. Additionally, we set the order times under the reference policy as the event times of a Poisson process with rate \(\lambda\). We let \(T_j\) denote the \(j\)th event time of the Poisson process for \(j=1,2,\dots\). Given these event times, we let \(Z_j \in \mathbb{R}^d,\) for \(j=1,2,\dots,\) be random vectors drawn independently from a common distribution \(\Phi\). In particular, \(Z_j\) will serve as the (random) order-up-to vector for the \(j\)th order under the reference policy. That is, \((Z_j - \tilde{X}(T_j-))^+\) denotes the order quantity associated with the order-up-to vector \(Z_j\); here, \((x)^+\) denotes the componentwise positive parts of \(x\in \mathbb{R}^d\), and \(\tilde{X}(T_j-)\) denotes the inventory vector just before time \(T_j\).

Randomizing the order-up-to vectors \(Z_j\) facilitates greater exploration of the state space and improves the training process. To further enhance state space exploration, we let \(\left(\zeta_j\right)_{j=1}^\infty\) denote a sequence of i.i.d random vectors in \(\mathbb{R}^d\), whose components are independent exponential random variables with mean \(\alpha_i = \alpha \mu_i/\lambda\) for \(i=1, \dots, d\). Then, we let \[Y_j = \max\left\{ Z_j - \tilde{X}(T_j-), \zeta_j \right\}\] denote the vector of order quantities of various items for the \(j\)th order, where the \(\max\{\cdot, \cdot\}\) is applied componentwise. This adjustment guarantees positive order quantities even if \(Z_j \ngeq \tilde{X}(T_j-)\), i.e., if \(Z_j\) does not dominate \(\tilde{X}(T_j-)\) componentwise, in which case some components of \((Z_j - \tilde{X}(T_j-))^+\) are zero. To ensure the stability of the reference process, we set \(\lambda \mathbb{E}[\zeta_j] < \mu_i\) for \(i=1,\ldots,d\), i.e., \(\alpha < 1\). Further details of the reference process parameters will be described in Section 5. Given \(\lambda, \Phi\) and \(\alpha\), we let \[\label{eq:uncontrolled95control} U(t; \alpha, \lambda, \Phi) = \sum_{j: T_j \le t} Y_j, \quad t \ge 0,\tag{11}\] and refer to the process \(U(\cdot; \alpha, \lambda, \Phi)\) as the reference policy. When the meaning is clear from the context, we suppress dependence on \(\alpha\), \(\lambda\), and \(\Phi\), and simply write \(U(\cdot)\). The corresponding reference process is given by \[\label{eq:uncontrolled95reference} \tilde{X}(t) = x - \mu t - \sigma B(t) + U(t), \quad t \ge 0.\tag{12}\]

Stochastic identity. We now present the equivalent probabilistic representation of the value function \(V\). Fix a reference policy \(U(\cdot)\). Then, for any \(T>0\) and \(\tilde{X}(0)=x\in\mathbb{R}^d\), \(V(x)\) can be formally characterized as follows: \[\label{eq:keyidentityprop} \begin{align} V(x) = \sup \Big\{ v \in \mathbb{R}\;:\; & v - \int_0^T e^{-r t} \nabla V(\tilde{X}(t))^\top \sigma dB(t) \\ & \leq e^{-r T} V(\tilde{X}(T)) + \int_0^T e^{-r t} f(\tilde{X}(t))dt + \sum_{j : T_j \leq T} c(Y_j)e^{-rT_j},\;\textrm{a.s.} \Big\}. \end{align}\tag{13}\]

The right-hand side of Eq. 13 defines a stochastic target problem, which, roughly speaking, is a control problem that aims to steer the system state so as to hit a target set at the terminal time \(T\). The connection between impulse control and stochastic target problems was first established in [65], for the case of discrete jump sizes. [6] later extended this connection for more general impulse control problems. For some intuition on the stochastic identity 13 , note that under suitable integrability conditions, the almost-sure inequality implies that \(V(x)\) satisfies \[V(x) \leq \E\left[\int_0^T e^{-r t} f(\tilde{X}(t))dt + \sum_{j : T_j \leq T} c(Y_j)e^{-rT_j} + e^{-r T} V(\tilde{X}(T)) \right].\] The right-hand side is the expected discounted cost under a policy that follows the reference policy on \([0,T]\) and acts optimally from time \(T\) onwards. Since \(V(x)\) is the optimal cost over all policies, any such policy must cost at least \(V(x)\). The surprising part is the exact characterization of \(V(x)\) as the largest value for which the almost-sure inequality holds; we provide a more substantive formal derivation in Appendix 9.

5 Computational Method↩︎

The identity 13 serves as the basis of our computational method. Seeking an approximate solution to it, we define the discretized loss function displayed in Eq. 15 . Recall the reference process defined in Eq. 12 . As a preliminary to defining the loss function, we fix a partition \(0=t_0<t_1<\ldots<t_N=T\) of the time horizon and simulate \(K\) discretized sample paths of the reference process at times \(t_0,t_1,\ldots,t_N\); see Subroutine 1.

Reference policy. As reviewed in Section 2, decades of research activity has culminated in effective heuristic policies for the SJRP. These include the \((R,S)\), \((Q,S)\), \((R,Q,S)\), and can-order policies, which are described in Section 6 and Appendix C.2. We consider their performances as benchmarks for our proposed policy. Crucially, they each have an order-up-to vector \(S\) as a policy parameter. In what follows, we set \[\label{eq:reference-S} \mathbb{E}[Z_j] = {\psi}S \quad \text{for } j \ge 1,\tag{14}\] where \(S\) is derived for one of these three benchmark policies, and \(\psi\geq1\) is a scaling factor. More specifically, assuming the optimal order-up-to levels are positive,¹ we set \(Z_j\) to be a vector of lognormal random variables, where its \(i\)th component has mean \(\psi S_i\) and variance \((\nu\psi S_i)^2\). Here, \(\nu > 0\) is a tuning parameter and is used to control the degree of exploration of the state space. We introduce the tuning parameter \(\psi\) to control the degree of exploration during the training process. By default, we set \(\psi=1\), unless the order-up-to vector \(S\) is too close to the origin, in which case we consider \(\psi>1\) to increase exploration. For virtually all test problems, we use \(\psi=1\) and report its value only when it differs from one; see Tables 8–13 in Appendix 10.2. Similarly, the parameter \(\alpha\) that determines the means of the exponential random vectors \(\zeta_j\), for \(j \ge 1\), is also viewed as a tuning parameter.

Loss function. Our method, described in Algorithm 2, computes the loss function defined in Eq. 15 by summing over the sample paths of the reference process and over discrete time steps to approximate the various integrals over \([0,T]\) appearing in the identity 13 . We approximate the value function \(V(\cdot)\) and its gradient \(\nabla V(\cdot)\) by deep neural networks, \(H_\theta\) and \(G_\vartheta\), respectively, with associated parameter vectors \(\theta\) and \(\vartheta\). Using the sample paths of the time-discretized reference process, the loss function \(L_\beta\) is defined as \[\label{eq:lossfunction} \begin{align} L_\beta & (\theta,\vartheta):= \frac{1}{K}\sum_{k=1}^{K}\Big[-H_\theta(x^{(k)}_0)+ \beta\cdot l\Big(H_\theta(x^{(k)}_0) - e^{-r T} H_\theta(\tilde{X}^{(k)}(T))\\ & - \sum_{n=0}^{N-1} e^{-r t_n} G_\vartheta(\tilde{X}^{(k)}(t_n))^\top \sigma\Delta B^{(k)}(t_n)-\sum_{n=0}^{N-1} e^{-r t_n} f(\tilde{X}^{(k)}(t_n)) \Delta t_n - \sum_{n=0}^{N-1} e^{-r t_n} c(\Delta U^{(k)}(t_n))\Big)\Big], \end{align}\tag{15}\] where the tuning parameter \(\beta\gg1\) can be viewed as a Lagrange multiplier. The penalty function \(l(\cdot)\) is given by \[\label{eq:squaredrectloss} {l(x)=\left[(x)^+\right]^2.}\tag{16}\] The function \(L_\beta\) consists of (the negative of) the value function estimate, \(-H_\theta(x_0^{(k)})\), and a penalty term, \(\beta l(\cdot)\), averaged over \(K\) initial states \(x_0^{(k)},\;k=1,\ldots,K\). As mentioned in Step 6 of Algorithm 2, the process \(\tilde{X}(\cdot)\) is simulated continuously, meaning that the terminal state of an iteration becomes the initial state of the next iteration. Loosely speaking, this approach can be seen as an attempt to capture the long-run behavior of the reference process. The penalty term penalizes violations of the almost sure inequality that appears in Eq. 13 . The higher \(\beta\) is, the lower the violation probability is. Note from Eq. 13 that we seek to choose the largest value \(v\) while satisfying the inequality there almost surely. Thus, the loss function \(L_\beta\) is designed so that as we minimize it, we increase the value function estimate \(H_\theta\), while simultaneously decreasing the violation probability. Also, we average over \(K\) initial states to avoid overfitting at any particular state. The tuning parameter \(\beta\) seeks to balance these two different objectives.

The penalty rate and learning rate schedules. Our computational method uses a stochastic gradient descent (SGD) type method to minimize the loss function. This involves choosing a learning rate schedule as done for deep learning algorithms. Additionally, we also need to choose the parameter \(\beta\) that is used to penalize the violations of the stochastic inequality on the right-hand side of Eq. 13 . As one would expect, the higher the value of \(\beta\), the lower the violation probability is. However, we observed empirically that the training is less stable with a higher value of \(\beta\). To deal with this instability, we adopt the approach of gradually increasing the penalty parameter \(\beta\) during training, while decreasing the learning rate as needed. More specifically, we begin with a relatively low value of the penalty parameter \(\beta\), at which we observe that a significant fraction of simulated sample paths violate the stochastic inequality. At the same time, we set the learning rate parameter relatively large and iterate until the loss appears to converge, lowering the learning rate if the training becomes unstable. Then, we increase the penalty parameter (e.g., by a factor of ten) and repeat the same procedure: We iterate until the loss appears to converge, and decrease the learning rate as needed if we observe instability. We keep increasing the penalty parameter in this manner until the violation probability is sufficiently low, e.g., less than 1%. For this final value of the penalty parameter, we choose a learning rate schedule manually by observing the evolution of the loss as we iterate, and lowering the learning rate as needed. See Appendix 10.1 for more details.

Using the approach just described, and through extensive experimentation, we identified a set of a priori fixed penalty rate and learning rate schedules — \(\{\beta_m\}\) and \(\{\eta_m\}\), where \(m\) indexes the iterations — to be considered in our experiments, as reflected in Algorithm 2. In particular, both schedules are treated as inputs to the algorithm and held fixed across all iterations of a given training run. As the reader will see, while Algorithm 2 formally allows \(\{\beta_m\}\) and \(\{\eta_m\}\) to vary with \(m\), in practice we restrict to piecewise-constant schedules, with values changing only at pre-selected small set of iterations. See Appendix 10.2 for additional details on the specific penalty parameter and learning rate schedules used for our test problems. Further implementation details of our computational method are provided in Appendix 10.

The proposed policy. After obtaining the neural network approximations \(H_\theta\) and \(G_\vartheta\) of the value function and its gradient, respectively, we propose a policy exploiting the structure of the HJB equation 10 . To be specific, Eq. 10 implies that either the no action term is zero, i.e., \[\label{eq:qvi95diffusion} \mathcal{L}V(x)+r V(x) - f(x) = 0,\tag{17}\] or the intervention term is zero, i.e., \[\label{eq:qvi95obstacle} V(x) = \inf_{y \in \mathbb{R}_+^d} \left\{V(x+y) + c(y)\right\}.\tag{18}\] In the former case, it is optimal to take no action. In the latter case, it is optimal to place an order, and the order vector is given by the minimizer of the right-hand side of Eq. 18 .

One way to operationalize this is to check if condition 17 holds. If 17 holds for some state \(x\), no action is required; otherwise, intervention through impulse control is warranted. That is, if for our neural network estimates of \(V\) and \(\nabla V\), we have \[\label{eq:qvi95diffusion95approx} \mathcal{L}V(x)+r V(x) - f(x) \leq -\epsilon < 0,\tag{19}\] then we intervene through impulse control. Here \(\epsilon\) is a suitably chosen numerical tolerance to account for optimization and sampling errors. To obtain an estimate of the value of \(\mathcal{L}V\), note that estimates of the value function \(V\) and its gradient \(\nabla V\) are already available from the deep neural network approximations, \(H_\theta\) and \(G_\vartheta\), respectively. It then remains to estimate the Hessian of \(V\) through, for example, automatic differentiation of \(G_\vartheta\) [66], [68].

If one decides to intervene, the impulse control \(y^*\) is determined by selecting \[\label{eq:argminimpulse} y^*\in \operatorname{argmin}_{y \in \mathbb{R}_+^d} \left\{H_\theta(x+y) + c(y)\right\}.\tag{20}\] To solve the minimization problem in 20 , one can use a variety of numerical optimization techniques. Further details are provided in Appendix 10.

We interpret this policy in the context of the original (discrete-time) formulation of the SJRP as follows: Given a system state \(x\), the warehouse manager checks if an intervention is warranted. If so, he places an order, and his order vector is \(y^*\) given in Eq. 20 ; see Algorithm 3 for a detailed description.

Simulation and optimization. A priori, under the proposed policy, every time that the warehouse manager decides to place an order, he needs to solve the minimization problem in Eq. 20 to determine the order size vector \(y^*\). This is computationally burdensome for performance simulations. However, the following observation allows us to compute a single order-up-to vector only once at the beginning of a simulation, resulting in a substantial speedup. Consider the minimization problem in Eq. 20 , and assume that the optimal solution \(y^* \neq 0\). Because the inventory state \(x\) is not a decision variable, the problem is in fact equivalent to \[\label{eq:argminimpulse95alt} \min_{y \in \mathbb{R}_+^d} \left\{H_\theta(x+y) + c(y) + \sum_{i=1}^d c_i x_i \right\},\tag{21}\] which can be further simplified and bounded as follows: \[\begin{align} \min_{y \in \mathbb{R}_+^d} \left\{H_\theta(x+y) + c(y) + \sum_{i=1}^d c_i x_i \right\} &~= \min_{y \in \mathbb{R}_+^d} \left\{H_\theta(x+y) + c(x+y) \right\} \nonumber\\ &~\geq \min_{z \in \mathbb{R}^d} \left\{H_\theta(z) + c(z) \right\}. \label{eq:single95order95up95to95level} \end{align}\tag{22}\] Here, the equality uses the fact that \(y^* \neq 0\), which ensures that \(c(y^*) + \sum_{i=1}^d c_i x_i = c(x+y^*)\). Let \(z^*\) be the optimal solution to the minimization problem in Eq. 22 . Then, as long as \(x \le z^*\) componentwise, the minimization problems in Eqs. 20 and 22 are equivalent. This in turn implies that under our proposed policy, the order-up-to vector is simply \(z^*\), irrespective of the initial state \(x_0\), as long as \(x_0\leq z^*\). Thus, if we assume that \(z^* \geq 0\), which holds in all our test problems, then, since we start our simulations from the initial state \(x_0=0\), the order-up-to vector \(z^*\) does not need to be recomputed at every ordering decision.

An alternative approach for deriving the optimal order quantity. One can alternatively use the neural network approximation of the gradient, \(G_{\vartheta}\), for solving 20 . In particular, the first-order optimality condition yields \(\nabla_z V(z) = -c\), where \(c\) is the vector of variable ordering costs. Thus, substituting \(G_\vartheta\) for \(\nabla V\) then yields \[\label{eq:stationary95conditionZ} G_\vartheta(z^*) = -c.\tag{23}\] Solving this system gives the optimal order-up-to vector. One can compare these approaches in two ways: (i) by fixing the same hyperparameter configuration for both methods a priori, or (ii) by considering the best-performing hyperparameter configuration for each method separately. Under the first comparison, the results of the two methods can differ substantially. Under the second comparison, however, the results are close, except for some test problems in the high-variability case. As mentioned earlier, since the order-up-to vectors only need to be computed once at the beginning of a simulation, running both optimization routines yields negligible overhead. Thus, we compute both vectors and select the one with the lower simulated cost. In practice, both optimization problems are solved over a compact box \([\underline{z},\overline{z}] := \{z\in\mathbb{R}_+^d:\underline{z}\leq z\leq\overline{z}\}\). The box is chosen large enough to contain a near-optimal order-up-to vector, while remaining small enough to restrict the search to regions that were thoroughly explored during training, where the neural network approximations are considered reliable; see Appendix 10 for further details.

To repeat, under natural conditions one needs to solve 20 or 23 only once as argued above. This is done in Step 1 of Algorithm 3 as the preprocessing step.

6 Test Problems and Benchmark Policies↩︎

We begin by describing the test problems. We consider a total of 43 test problems² that can be divided into three groups. First, we consider 27 twelve-dimensional problems based on an example studied in [30]. These problems vary across three attributes: Demand variability, fixed ordering cost, and backlog penalty costs. For each attribute, we consider low, medium, and high values, yielding 27 test problems.

Computing the optimal policy for these problems is intractable due to the curse of dimensionality. Thus, to assess the performance of our proposed policy, we use the state-of-the-art heuristic policies drawn from the literature as benchmarks; see below for more on those benchmark policies. However, when the number of items is small, e.g., \(d=2\), one can compute the optimal policy using standard dynamic programming methods. This, of course, helps assess how close our proposed policy is to the optimal one for low-dimensional problems. Therefore, we next consider 7 test problems of dimension two. To do so, we consider only two products of the 12-dimensional problem. For our base case of the two-dimensional test problems, we choose the medium values of demand variability, ordering cost, and backlog penalty costs. Additionally, while fixing all else we consider setting each of the three attributes, one at a time, to low and high values, resulting in a total of seven test problems of dimension two. Lastly, to assess the scalability of our method, we consider nine test problems of dimension 50 that are designed by building on the earlier test problems. More will be said about them below.

For all test problems, we set the time unit as a year, the initial inventory to zero, and the annual interest rate to \(r=5\%\), and allow the warehouse manager to make ordering decisions weekly.

12-dimensional test problems. We build on the 12-item problem studied in [30]. We set the holding cost parameters \(h_i=2\) for all items (\(i=1,\dots,12\)). To capture the three levels of demand variability, we proceed as follows: For low variability, we model the item-level demand as independent Poisson processes. In this case, the coefficient of variation (CV) of item-level annual demand ranges from \(0.158\) (for items with annual demand rate \(40\)) to \(0.224\) (for items with annual demand rate \(20\)); cf. Table 1. For medium and high variability cases, we model each item-level annual demand as an independent negative binomial random variable, and set their parameters so that the CV of the annual demand is \(0.5\) per item for the medium variability case, and it is \(1\) per item for the high-variability case. We choose the negative binomial distribution because it can effectively model both high and low demand rates, and its dispersion parameter can be used to fit a wide range of empirical retail demand distributions [69], [70]. Table 1 shows the annual demand rates and variable costs for each item. As mentioned earlier, we consider three possible values for the fixed cost of ordering and the backlog penalty costs. More specifically, we set \(p_i=p\) for all \(i\) and let \(p \in \{10, 50, 100\}\) and \(c_0 \in \{20, 100, 200\}\).

Low-dimensional test problems. Letting \(d=2\), the parameters of our test problem are taken from the preceding one. Namely, we focus on items 1 and 7 there. For the base case, we set \(c_0=50\) and \(p=50\), and model item-level annual demand as a negative binomial random variable with annual CV \(=0.5\), i.e., the medium variability case. The other six 2-dimensional test problems follow from varying one of the following at a time: The demand variability, the fixed cost of ordering, and the backlog penalty parameter. More specifically, these alternative parameters correspond to demand variability \(\in \{\text{low, high}\}\), \(c_0 \in \{20, 100\}\), and \(p \in \{10, 100\}\). For the low variability case, we model the item-level demands as independent Poisson processes as done earlier. For the high-variability case, item-level demands follow negative binomial distributions with annual CV \(=1.0\).

High-dimensional test problems. As mentioned earlier, we also consider nine 50-dimensional test problems. To set the problem parameters, we partition the items into three groups. To be specific, items 1-15 form the first group, items 16-30 form the second group, and items 31-50 form the third group. Within each group, items share the same annual demand rate, per-unit holding cost, backlog penalty, and variable cost; these parameters differ across groups. We vary the fixed cost \(c_{0}\in\{50,150,250\}\), and again model demand as Poisson (low variability) or negative binomial with an item-level CV of 0.5 (medium variability) or 1.0 (high variability). These lead to nine additional test problems. Table 2 summarizes the key parameters for each group.

Next, we describe the benchmark policies we use to assess the performance of our proposed policy. Crucially, we make all performance comparisons in the context of the original discrete-time formulation, cf. Eqs. 1 –5 , as opposed to the impulse control approximation 7 –9 , even though our proposed policy is derived using the latter.

Benchmark policies. As mentioned earlier, it is computationally feasible to solve the low-dimensional test problems. Thus, we compute the optimal policy using policy iteration and use it as the benchmark for those cases. For the other test problems, we rely on effective heuristic policies drawn from the extensive literature on the SJRP. More specifically, we consider the following benchmark policies:

Under the \((R,S)\) policy, the inventory levels of all items are replenished every \(R\) time periods, with each item \(i = 1, \dots, d\) having its inventory level brought back to its target \(S_i\). The \((Q, S)\) policy, studied in [71], places an order whenever the total aggregate demand since the previous replenishment exceeds a specified value \(Q\), raising the inventory, for each item \(i = 1, \dots, d\), back to its order-up-to level \(S_i\). The \((R, Q, S)\) policy, proposed by [72], builds on the \((R, S)\) and \((Q, S)\) policies, and can be described as follows: It places an order whenever (a) the total aggregate demand since the previous replenishment reaches \(Q\), or (b) \(R\) time units have elapsed since the last decision epoch, whichever occurs first, raising each item \(i = 1, \dots, d\) to its order-up-to level \(S_i\). Under the can-order policy, each item \(i\) is managed with three control parameters: A reorder point \(s_i\), a can-order level \(c_i\), and an order-up-to level \(S_i\). For each item \(i = 1, \dots, d\), a replenishment is triggered when its inventory level drops to or below \(s_i\), bringing it back up to \(S_i\). If at this point the inventory level of any other item \(j \neq i\) is at or below its can-order level \(c_j\), this item \(j\) is also included in the replenishment, and its inventory is raised to \(S_j\). The independent \((s,S)\) policy assumes the setting without any benefits from joint ordering and calculates the individual optimal \((s,S)\) levels as if there were no opportunity for group savings on the fixed cost. Further implementation details of these benchmark policies are provided in Appendix 11.2.

7 Computational Results↩︎

For the test problems introduced in Section 6, we now compare the performance of our proposed policy, henceforth referred to as the neural network (NN) policy, to the benchmarks. For the benchmark policies, i.e., \((R,S)\), \((Q,S)\), \((R,Q,S)\), can-order, and independent \((s,S)\), we compute their performance estimates by averaging over 100,000 sample paths of length 10,000 time periods each. The resulting standard errors are negligible (\(\leq0.01\%\) of the respective estimated performance), so we do not report them. For the NN policy, depending on the problem dimension, we estimate its performance by averaging over 100 or 1,000 sample paths of length 10,000 periods each; further details are provided below. The numbers of simulated sample paths under the NN policy are chosen to strike a balance between computational tractability and statistical accuracy: On the one hand, checking the no-action condition Ineq. 19 is computationally intensive, as it involves evaluating the approximate Hessian of the value function through automatic differentiation; on the other hand, we simulate enough paths to obtain standard errors that are less than \(1\%\) of the respective estimated performance. For a test problem, we say that the NN policy matches a benchmark if the estimated benchmark cost is within \(1\%\) of the estimated cost under the NN policy, i.e., \[\frac{\left|~\text{Estimated benchmark cost} - \text{Estimated NN cost}~\right|}{\text{Estimated NN cost}} \leq 1\%,\] and that the NN policy beats a benchmark if the estimated benchmark cost is at least \(1\%\) more than the estimated NN cost, i.e., \[\frac{\text{Estimated benchmark cost} - \text{Estimated NN cost}}{\text{Estimated NN cost}} > 1\%.\]

On our hardware, each neural network training takes roughly 2.5 hours in 12 dimensions and 4 hours in 50 dimensions. On a single CPU, the performance simulations take roughly 6 hours in 12 dimensions and 10 hours in 50 dimensions. Computation time thus grows only modestly with the problem dimension and does not exhibit exponential growth in \(d\). This is also because we use neural networks whose capacity scales modestly with the problem dimension (e.g., three hidden layers for the 12-dimensional instances and four for the 50-dimensional instances). Further implementation details, such as machine configurations, programming language, and software packages, are described in Appendix 10.1. Precise training and simulation times for each instance are reported in Appendix 10.2. Our code is available at https://github.com/wvaneeke/DeepSJRP.

Before presenting the computational results in detail, we first summarize the main findings. Overall, the NN policy performs strongly and often surpasses the benchmarks by a significant margin. Across the 36 high-dimensional test problems (27 in dimension twelve and 9 in dimension fifty), it matches or beats the best benchmark, and does so by a clear margin in most medium- to high-variability cases. In all seven two-dimensional test problems, the performance gap between the NN policy and the MDP solution is within \(1\%\), and the NN policy structure closely aligns with that of the MDP, suggesting that the NN policy learns the optimal policy in low dimensions. We now proceed to the detailed results.

7.0.0.1 Results in dimension 12.

Tables 3, 4, and 5 report computational results in the low, medium, and high variability settings, respectively, for the main 12-dimensional test problems. Each table reports the simulated performances of our NN policy, together with relative gaps between the various benchmark policies and our policy. The performance of the NN policy is estimated by averaging over 100 sample paths. The standard errors are shown as percentages of the estimated averages. Reported standard errors are small (\(<1\%\)), indicating statistically reliable estimates.

For the nine test problems with low variability in Table 3, our NN policy performs on par with the best benchmark policy, which turn out to be \((Q,S)\) and \((R,Q,S)\) in all cases (these two policies in fact attain identical costs to within simulation error across all instances in our study). In each case, the optimal \(R\) is large enough that orders are nearly always triggered by the \(Q\) threshold rather than by the periodic review, so \((R,Q,S)\) effectively reduces to \((Q,S)\). The same is true for all subsequent test problems as well. Thus, in the discussions that follow, we do not discuss \((R,Q,S)\) explicitly. The NN policy outperforms the can-order policy by a significant margin in most cases, but the performance gap shrinks as the backlog penalty increases. Intuitively, the shrinking performance gaps can be explained by the can-order policy’s tendency to order more proactively than other benchmark policies, suggesting that it is better at avoiding states with large backlogs and performs better as the backlog penalty increases. The NN policy also significantly outperforms the independent \((s,S)\) policy. This is expected, as the independent \((s,S)\) policy does not exploit cost-saving opportunities of joint replenishment and instead triggers an order each time an item reaches its reorder point individually.

We now focus on the performance comparisons among NN, \((R,S)\), and \((Q,S)\) in Table 3. For all values of the fixed cost, the NN policy tends to perform better when compared to \((R,S)\) and \((Q,S)\), as the backlog penalty increases. A possible explanation is that the neural networks are flexible enough to learn the shape of the ordering boundary well, so under high backlog penalties, the NN policy can better avoid costly states with large backlogs. When the fixed cost is high, it is optimal to place orders less frequently, allowing demand to pool between replenishments and thereby reducing effective variability. In this regime, the system behaves nearly deterministically, and we expect \((Q,S)\) and \((R,S)\) to be close to optimal. The intuition that \((R,S)\) is near-optimal is threefold: (i) Infrequent orders imply that aggregate inter-order demand is less variable, bringing the system closer to being deterministic; (ii) in a deterministic setting, it is optimal to place orders at regularly spaced times; and (iii) if optimal order times are approximately regular, a policy that optimizes order-up-to levels and (deterministic) review times, namely \((R,S)\), is near-optimal. The slight advantage of \((Q,S)\) over \((R,S)\) likely stems from its additional adaptivity in timing orders. Consequently, in a system with low effective variability (low demand variability and large fixed cost), there is little room for the NN policy to improve upon these benchmarks.

In the medium-variability setting (Table 4), the NN policy often delivers sizable gains over \((Q,S)\) and \((R,S)\), with the advantage growing as the backlog penalty increases. This aligns with the intuitive explanation discussed in the low-variability setting: The NN policy appears to learn the ordering boundary more accurately. The NN policy also dominates can-order in all cases, though the performance gap narrows as backlog penalties rise; again this observation is consistent with our earlier explanation on can-order’s proactive replenishments. Finally, as expected, the NN policy substantially outperforms the independent \((s,S)\) policy. In the high-variability setting (Table 5), the qualitative pattern is similar to the medium-variability cases, with the NN policy matching or beating the best benchmark in all cases, and often with high margins.

7.0.0.2 Low-dimensional results.

Table 6 and Figures 4 and 5 summarize our computational results in the two-dimensional test problems. To estimate performance of the NN policy, we average over 1,000 sample paths, each of length 10,000 periods. As mentioned earlier, the NN performance estimates are within \(1\%\) of the MDP solutions. Across the base case (Figure 4) and the six variants (Figure 5), the NN policy closely tracks the MDP solution: The decision boundaries largely overlap (dense purple regions), with minor local deviations near extremely-valued states. These visual observations are consistent with the small optimality gaps in Table 6 (\(\approx\) \(-0.01\%\) to \(-0.69\%\) relative gap between the MDP solution and the NN policy across different cases), indicating that not only does the NN policy perform near-optimally, but it also captures the geometry of the optimal policy. Among the benchmark policies, the can-order policy performs the best and closely matches the MDP solution in all seven test problems. Recall the intuition that a good reference policy should closely track the optimal policy. Indeed, the can-order policy is a strong candidate, and we used it as the reference policy for training in most instances. One exception is the high-variance case (\(\mathrm{CV} = 1.0\), see the 3rd row of Table 6). In this case, we use the optimal policy derived from the MDP. That is, we use its order-up-to vector as the mean of \(Z\) (see Eq. 14 ), and compute the average number of orders per year and set it as the arrival rate \(\lambda\) of the Poisson process that underlies our reference process. Among the reference policies we considered, this led to the best performance, supporting the intuition that under a good reference policy, sample paths of the state process should occupy parts of the state space most frequently visited by an optimal policy. Still, the can-order policy is slightly cheaper than the NN policy in this case (\(-0.09\%\)). Intuitively, because high demand variability makes large backlogs more likely, the can-order policy’s proactive replenishment is particularly effective in this regime.

7.0.0.3 Results in dimension 50.

Finally, we turn to the 50-dimensional test problems. We estimate NN performance by averaging over 100 sample paths per test problem. Standard errors are small (\(\approx 0.05\% - 0.20\%\)), indicating reliable estimates. The NN policy beats the best benchmark in all cases, with performance gaps up to 5.59%. Under low variability, NN performs comparably to \((R,S)\) and \((Q,S)\), while clearly outperforming can-order and independent \((s,S)\) policies. As variability increases, NN’s advantage over \((R,S)\) and \((Q,S)\) generally widens, whereas gaps versus can-order and independent \((s,S)\) remain relatively stable.

8 Concluding Remarks↩︎

In this paper, we tackle a classical inventory control problem in high dimensions by formulating it as an impulse control problem and approximating its solution via deep-learning techniques. Through extensive hyperparameter tuning, a tailored training algorithm, and careful engineering of policy extraction, we were able to match or beat the best available benchmarks. We note that the vast array of SJRP policies surveyed in [3] is largely tailored to settings that also feature item-specific minor fixed costs, with the modified periodic review policy of [30] serving as a representative example. Nonetheless, under the single joint fixed cost structure we consider, these policies either reduce to or are dominated by one of our main benchmarks: For example, the modified periodic review policy of [30] and the \((R, s, S)\) policy of [35] are both dominated by the optimized \((R, S)\) policy, and the \((Q, s, S)\) policy of [73] is dominated by the optimized \((Q, S)\) policy.

While we focus on the SJRP with linear holding and backlogging costs and unconstrained order vectors (subject only to non-negativity constraints), the methodology extends naturally in several directions. Within the present framework, it already accommodates general Lipschitz continuous inventory cost functions \(f\) (cf., e.g., growth conditions (2.4) and (2.5) in [6]) as well as capacity constraints on order vectors, for example, requiring the total order size not to exceed a given capacity. Another natural extension is to the lost sales case, where unmet demand is forfeited rather than backlogged. In this setting, the inventory state is constrained to the nonnegative orthant by minimally pushing it back into the state space whenever it reaches the boundary. The associated pushing process on each face of the orthant represents the cumulative lost demand for the corresponding product. Under this formulation, one accounts for lost-demand costs rather than the backlogging costs considered here. The resulting HJB equation is therefore posed on the nonnegative orthant and is accompanied by Neumann-type boundary conditions arising from the pushing, or reflection, at the boundary. Despite these differences, we expect that our approach can be extended to address this problem, although doing so would require resolving a number of technical subtleties. We leave this extension to future work.

Beyond these extensions, incorporating lead times represents another important direction for future research. In this regard, [74] study an impulse control problem with fixed delivery lags. They show that the optimal policy depends only on the aggregate state, defined as the sum of on-hand inventory and outstanding orders. They also show that the HJB equation takes a modified form with an intervention operator that averages holding and backlog costs over demand during the lead time. These structural results suggest that our framework may extend to this setting, though the requisite modifications to the stochastic identity 13 remain to be developed. Looking further ahead, it seems natural to extend our methodology to more complex supply-chain models, such as those with demand spillovers between facilities.

Complementary to these problem extensions, a natural direction concerns our use of the Brownian approximation for cumulative demand: Reflected Brownian motions or subordinator-like processes are potential alternatives. Beyond the inventory setting itself, our approach not only provides effective solutions to supply chain and inventory management problems, but also opens up avenues for applying these numerical methods to broader classes of impulse control problems. In particular, the techniques developed here may be used to tackle impulse control problems in other fields, most notably in financial engineering — for example, cash flow management and investment decisions [75], [76]. We leave these explorations for future work.

References↩︎

---

1. This is a reasonable assumption for our test problems because the backlog costs are considerably larger than the holding cost.↩︎

2. Additionally, Appendix D considers a one-dimensional impulse control problem studied in [41] and verifies the validity of our solution method by comparing it to the one derived in [41]; see Figure 7 in Appendix 12.↩︎

3. For clarity within this appendix, \(\beta_1, \ldots, \beta_M\) and \(\eta_1, \ldots, \eta_K\) denote the distinct values taken by the piecewise-constant penalty rate and learning rate schedules, so \(M\) and \(K\) here count distinct levels rather than the total number of iterations and the batch size as in Algorithm 2.↩︎