# 2-approximation algorithms for two variants of the static stochastic joint replenishment problem

Guillaume Massonnet, Gautier Stauffer

► To cite this version:

Guillaume Massonnet, Gautier Stauffer. 2-approximation algorithms for two variants of the static stochastic joint replenishment problem. Operations Research Letters, 2025, 62, pp.107321. 10.1016/j.orl.2025.107321. hal-05509909

HAL Id: hal-05509909

https://hal.science/hal-05509909v1

Submitted on 13 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# 2-approximation algorithms for two variants of the static stochastic joint replenishment problem

Guillaume Massonnet<sup>a, ID</sup>, Gautier Stauffer<sup>b, ID, *</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

**Keywords:**

Approximation algorithm

Stochastic JRP

Static strategy

## ABSTRACT

In this work, it is demonstrated that the technique developed by Gayon et al. (2017) [11] for designing 2-approximation algorithms for the one-warehouse multi-retailer (OWMR) problem is highly versatile and can be applied to stochastic variants of the original OWMR problem. The mechanics of this approach are illustrated on two (static) stochastic variants of the Joint Replenishment Problem (JRP). Additionally, the simplicity and versatility of this technique suggest its potential for broader applications in the stochastic setting.

## 1. Introduction

The joint replenishment problem is a classic model in inventory management. It captures the economies of scale achieved when the orders for different products are synchronized, allowing to mutualize some of the costs incurred by common operations such as transportation or handling. Standard models include a fixed cost each time an order is placed (regardless of the product mix) and additional fixed costs for each product involved in the corresponding order (regardless of the quantity). We consider the discrete-time version of the problem with time-varying demand-requests that are released in a finite set of periods when orders can also be placed. We assume linear costs for holding inventory and linear penalties for backorders. The deterministic version of the problem has been widely studied in the literature and is known to be NP-hard [2,6], but there exists exact mixed integer programming formulations suitable to solve medium size problems [9,23,21], efficient approximation algorithms [15,19,10,24] and good heuristics [22]. All the methods mentioned above assume perfect information about future customers demands, which is unrealistic in most practical settings where this knowledge often takes the form of probabilistic forecasts. Considering demand uncertainty, on the other hand, calls for solution methods that can tackle a stochastic version of the JRP, which still remains a major challenge. Until now, the literature focuses generally on optimizing parameters of pragmatic policies (e.g. base-stock, can-order) through stochastic programming formulations, or on solving special cases when assumptions are made on the demand structure (Poisson distribution, balance assumption, etc.) [7]. In fact, the structure of the optimal poli-

cies is unknown in general [7] under static, static-dynamic and dynamic uncertainty strategies [4]. For a comprehensive survey on the JRP and variants, see [18,12,3,26,20,7].

In order to derive good solutions within a limited computational time, approximation algorithms have proven to be viable alternatives in the realm of stochastic inventory management and have gained significant traction over the past decades following the pioneering work of [14,13,17,16]. However, unlike the deterministic setting, no approximations are known for any variant of the stochastic joint replenishment problem. This paper demonstrates that the 2-approximation algorithm based on the *decomposition/uncrossing/recombination* technique introduced in [11] for the deterministic JRP, naturally extends to variants of the *Static Stochastic Joint Replenishment Problem (SSJRP)*. The method comprises three main steps:

**Step 1: Decomposition.** Several subproblems, a global one and one for each item, whose sum of optimal values provides a lower bound on the value of the optimum, are solved independently.

**Step 2: Uncrossing** The ordering periods in the solutions to the single-item subproblems are adjusted and synchronized with a set of candidate orders from the remaining (global) subproblem, while keeping ordering costs under control.

**Step 3: Recombination** A feasible solution to the original problem is built by recombining elements of solutions from each subproblem, while keeping holding and backorder costs under control.

<footnote><sup>a</sup> IMT Atlantique, LS2N, 44300 Nantes, France</footnote>

<footnote><sup>b</sup> Faculty of Business and Economics (HEC Lausanne), Department of Operations, University of Lausanne, Quartier Unil-Chamberonne, 1015 Lausanne, Switzerland</footnote>

<footnote>* Corresponding author.
E-mail addresses: guillaume.massonnet@imt-atlantique.fr (G. Massonnet), gautier.stauffer@unil.ch (G. Stauffer).</footnote>
The stochastic variants considered in the remainder of this paper are based on two natural models for demand stochasticity that are presented in [25] and [5,1] for single-echelon problems, respectively. The former addresses a single-item dynamic lot-sizing problem with non-negative and independent stochastic demands with a known distribution. It computes an optimal production plan by extending the Wagner-Whitin's shortest path approach, where the costs of the arcs in the underlying graph are determined by computing the convolution of demands over time. This requires a precomputing phase that involves one-dimensional numerical integration to evaluate the integral of the corresponding convolution functions at specific points, before solving a one-dimensional convex optimization problem to compute a unique fractile of a weighted sum of these functions and obtain the optimal inventory level in each period. These tasks can be performed efficiently and with high accuracy using simple techniques such as rectangle rules for numerical integration or the bisection method for convex optimization. In certain cases, closed-form formulas or tabular reference functions may be used, as pointed in [25], e.g. for normally distributed demands. The second type of demand uncertainty, introduced in [5] and later in [1] for a single-echelon dynamic lot-sizing problem, is related to the timing of the demand rather than the quantity. In short, every demand occurs entirely during a single period within a specific time interval according to a given distribution and must be met completely before the end of the interval (no partial delivery is allowed). Under some practical assumptions, this problem can again be solved efficiently by a reduction to a shortest path problem in the vein of the Wagner-Whitin's approach.

## 1.1. The setting

We define $[k] = \{1, \dots, k\}$ for any integer $k \in \mathbb{N}$. We also introduce the indicator function $\delta: \mathbb{R} \to \{0, 1\}$, such that $\delta(x) = 1$ for $x > 0$, and $\delta(x) = 0$, otherwise. In the following we assume that we have a finite horizon of $T$ periods of time and $n$ products labeled $1, \dots, n$ with (stochastic) demands $D_i^t$, for $i \in [n]$ and $t \in [T]$. In order to fulfill the demands, we need to order the items from a supplier. Each time we place an order, at the beginning of a period, we pay a fixed cost $K_0$ independently of the subset of items and a fixed cost $K_i$ if we order some items of product $i$, for $i \in [n]$ (whatever the quantity). $K_0$ is called the joint set-up cost while $K_i$, $i \in [n]$ are called the individual set-up costs. Hence if we order a subset $S$ of items with $\emptyset \neq S \subseteq [n]$, we pay $K_0 + \sum_{i \in S} K_i$ in ordering costs. We assume that there is no lead time: orders are received at the beginning of the period where they are ordered. There might be on hand inventory $s_i \ge 0$ at the beginning of the horizon, for each product $i \in [n]$. Demand that cannot be satisfied from stock is backlogged. We also assume that at the end of period $t$, we pay a holding cost $h_t^i \ge 0$ per unit of item $i$ in stock and a per unit penalty cost of $\pi_t^i \ge 0$ for each unit of item $i$ that is backlogged. We assume additionally that $\pi_t^i = \alpha^i h_t^i$ for some $\alpha^i \ge 0$. We want to find an optimal static replenishment strategy, that is, we want to decide, before the first period, which quantity $q_t^i$ of each item $i \in [n]$ to order/receive at the beginning of each period $t \in [T]$, to minimize the sum of the ordering costs plus the expected holding and backorder costs over the different realizations of the demand. The replenishment strategy is said to be static as it does not depend on the realization of the demand. Although this type of policy might be suboptimal with respect to more dynamic/adaptative strategies, it has some virtue in practice as it allows to anticipate capacity requirements for instance (especially when the model is used in production settings as discussed in [25]).

## 1.2. A useful uncrossing procedure

Before closing this section we introduce a simple function $U: 2^{[T]} \times 2^{[T]} \to 2^{[T]}$ that uncrosses two sets of ordering periods contained in $[T]$ and that will be used for each of the two demand models considered later. We first need to introduce two auxiliary functions $\tau, \bar{\tau}$:

$[T] \times 2^{[T]} \to [T] \cup \{0, T+1\}$ defined as follows. For $t \in [T]$ and $\Theta_0 \subseteq [T]$, we define:

$$ \underline{\tau}(t, \Theta_0) = \max\{\theta \in \{0\} \cup \Theta_0 : \theta \le t\} \quad (1) $$

$$ \bar{\tau}(t, \Theta_0) = \min\{\theta \in \Theta_0 \cup \{T+1\} : \theta \ge t\} \quad (2) $$

$\underline{\tau}(t, \Theta_0)$ returns the largest period in $0 \cup \Theta_0$ that is less than or equal to $t$, while $\bar{\tau}(t, \Theta_0)$ returns the smallest period in $\Theta_0 \cup T+1$ that is greater than or equal to $t$. We call $\Theta_0$ the reference set. In the remainder, we sometimes omit the reference set $\Theta_0$ when it is clear from the context and simply refer to these functions as $\underline{\tau}(t)$ and $\bar{\tau}(t)$. For two subsets of periods $\Theta$ and $\Theta_0$, $U(\Theta, \Theta_0)$ is defined as follows:

$$ U(\Theta, \Theta_0) = \left( \bigcup_{\theta \in \Theta} \{\tau(\theta, \Theta_0), \bar{\tau}(\theta, \Theta_0)\} \right) \setminus \{0, T+1\} \quad (3) $$

**Claim 1.** Let $T \in \mathbb{N} \setminus \{0\}$. Let $\Theta = \{\theta_1, \dots, \theta_m\} \subseteq [T]$ for some $m \in \mathbb{N}$ with $\theta_1 < \dots < \theta_m$ and $\Theta_0 = \{\theta_1^0, \dots, \theta_{m_0}^0\} \subseteq [T]$ for some $m_0 \in \mathbb{N}$ with $\theta_1^0 < \dots < \theta_{m_0}^0$. Suppose $U(\Theta, \Theta_0) = \{\rho_1, \dots, \rho_r\}$ for some $r \in \mathbb{N}$ with $\rho_1 < \dots < \rho_r$. For convenience we introduce $\theta_0 = \theta_0^0 = \rho_0 = 0$ and $\theta_{m+1} = \theta_{m_0+1}^0 = \rho_{r+1} = 0$. Then for all $k \in \{0, \dots, r\}$, we have either:

(i) $\rho_k$ and $\rho_{k+1}$ coincide with two periods $\theta_0^0$ and $\theta_{l+1}^0$ for some $l \in \{0, \dots, m_0\}$, or

(ii) $\rho_k$ and $\rho_{k+1}$ are contained in an interval $[\theta_l, \theta_{l+1}]$ for some $l \in \{0, \dots, m\}$.

Moreover $U(\Theta, \Theta_0) \subseteq \Theta_0$ and $|U(\Theta, \Theta_0)| \le 2 \cdot |\Theta|$ by construction.

**Proof.** Let $k \in \{0, \dots, r\}$. We are in one of two cases: either there exists $\theta \in \Theta$ in the open interval $(\rho_k, \rho_{k+1})$ or there is none. In the first case, it implies that $\rho_k = \tau(\theta)$ and $\rho_{k+1} = \bar{\tau}(\theta)$ coincide with two periods $\theta_0^0$ and $\theta_{l+1}^0$ for some $l \in \{0, \dots, m_0\}$. In the second case, there exists $l \in \{0, \dots, m\}$ such that $\rho_k = \bar{\tau}(\theta_l)$ (with the convention $\bar{\tau}(0) = 0$). Since no element of $\Theta$ is contained in $(\rho_k, \rho_{k+1})$, we have that $\theta_{l+1} \ge \rho_{k+1}$ and the result follows. $\square$

Fig. 1 gives an example.

# 2. Stochastic volume model

## 2.1. The model

We assume here, following [25], that the demands $D_i^t$, for $i \in [n]$ and $t \in [T]$, are mutually independent random variables and follow known probability density functions $g_i^t$. At the end of the horizon, remaining on-hand inventory can be sold at production cost (inducing no additional costs) and backlogged demands (if any) are left unfilled.<fnref n="1" /> Following [25] for the single product setting, we consider equivalently the problem of deciding the cumulative quantity $S_i^t := \sum_{l' \in [l]} q_{l'}^t$ to order in or before period $t$ for each product $i$. For all $i \in [n]$, let $S_i^t = s^i$ be the inventory of product $i$ available at the beginning of the horizon. Our goal is to find a static replenishment strategy that minimizes the expected cost incurred over the $T$ planning periods. That is, we want to minimize the following expression over all $S_i^t$ such that $s_i \le S_i^t \le \dots \le S_T^t$ for all $i \in [n]$ and $t \in [T]$:

<footnote n="1">¹ Note that in a static framework we cannot add an order in the last period in scenarios where there is backorder as this would be a dynamic approach. However we could still consider such an extension in practice.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Example of output $U(\Theta, \Theta_0) = \{\rho_1, ..., \rho_5\}$ for $\Theta_0 = \{\theta_1^0, ..., \theta_6^0\}$ and $\Theta = \{\theta_1, \theta_2, \theta_3\}$. Case (i) of Claim 1 applies to $k \in \{1, 3, 4\}$, while case (ii) applies to $k \in \{0, 2, 5\}$.</figcaption></figure>

$$ \sum_{t \in [T]} K_0 \delta \left( \sum_{i \in [n]} (S_t^i - S_{t-1}^i) \right) + \sum_{i \in [n]} \left( K_i \delta(S_t^i - S_{t-1}^i) + h_t^i \int_{q=0}^{S_t^i} (S_t^i - q) f_t^i(q) dq + \pi_t^i \int_{q=S_t^i}^{\infty} (q - S_t^i) f_t^i(q) dq \right) $$

where $f_t^i(q)$ is the convolution of the demand of product $i$ for periods 1 to $t$. We call this problem (P).

## 2.2. Solving the single-product version using Vargas' approach (slightly revisited)

We shortly recall the solution approach proposed in [25] when there is only one product, say $i$, and no joint-ordering cost $K_0$. The single-product problem ($P_i$) consists in minimizing

$$ \sum_{t \in [T]} K_i \delta(S_t^i - S_{t-1}^i) + h_t^i \int_{q=0}^{S_t^i} (S_t^i - q) f_t^i(q) dq + \pi_t^i \int_{q=S_t^i}^{\infty} (q - S_t^i) f_t^i(q) dq $$

over all $S_t^i$ such that $s_i \le S_t^i \le \dots \le S_T^i$ for all $i \in [n]$ and $t \in [T]$.

($P_i$) reduces to a shortest path problem in a directed acyclic graph as follows. The idea is to build a directed graph $G = (V, A)$ (à la Wagner-Within) with node set $V = \{0, 1, ..., T + 1\}$ (nodes 0 and $T + 1$ are artificial nodes representing the beginning and the end of the horizon, respectively) and arc set $A = \{(u, v) \in V^2 : u < v\}$. We let

$$ \lambda^i(S^i, u, v) := \sum_{t=u}^{v-1} (h_t^i \int_{q=0}^{S_t^i} (S^i - q) f_t^i(q) dq + \pi_t^i \int_{q=S_t^i}^{\infty} (q - S^i) f_t^i(q) dq) \quad (4) $$

$\lambda^i(S^i, u, v)$ represents the holding and backlogging penalty costs incurred over periods $\{u, ..., v-1\}$ when (i) the cumulative quantity ordered in or before period $u$ is equal to $S^i$ and (ii) no order is placed between periods $u$ and $v-1$. We express the cost of an arc $(u, v) \in A$ as

$$ c^i(u, v) = \begin{cases} \lambda^i(s_i, u, v) & \text{if } u = 0 \\ \min_{S^i \ge s_i} K_i + \lambda^i(S^i, u, v) & \text{otherwise} \end{cases} \quad (5) $$

For all $0 < u < v \le T + 1$, $\lambda^i(S^i, u, v)$ is the sum of newsvendor formulas and thus inherits convexity and differentiability. Moreover, it also has a unique minimum as shown by Vargas in [25]. We denote by $S_*^i(u, v)$ the unique minimizer, with the convention that for all $v > 0$, $S_*^i(0, v) = s_i$. $S_*^i(u, v)$ is the solution to a one-dimensional convex optimization problem and can be computed with arbitrary precision with a classical numerical method (e.g. the bisection method).

Proposition 1 in [25] shows<fnref n="2" /> that $S_*^i(u, v) \le S^i(u, v)$ for all $0 < u < v < w \le T + 1$. Letting $\bar{S}^i(u, v) := \max\{s_i, S_*^i(u, v)\}$ for all $0 \le u < v \le T + 1$,<fnref n="3" /> we have:

**Claim 2.** For all $0 \le u < v < w \le T + 1$, $\bar{S}^i(u, v) \le \bar{S}^i(v, w)$.

Consider a path $P$ from 0 to $T + 1$ in $G$ and let $0 = \theta_0^i < \theta_1^i < \dots < \theta_m^i < \theta_{m+1}^i = T + 1$ be the nodes on the path. By Claim 2, the solution defined by $S_{\theta_l}^i = S_{\theta_{l+1}}^i = \dots = S_{\theta_{l+1}-1}^i = \bar{S}^i(\theta_l, \theta_{l+1})$ for all $l = 0, \dots, m$ is a solution of same cost as $P$ to ($P_i$). Conversely, consider an optimal solution $(S_t^i)_{t \in [T]}$ to ($P_i$) and let $\theta_1^i < \dots < \theta_m^i$ be the periods $t \in [T]$ for which $S_t^i > S_{t-1}^i$ (with the convention $S_0^i = s_i$). We can assume w.l.o.g. that $S_{\theta_l}^i = S_{\theta_{l+1}}^i = \dots = S_{\theta_{l+1}-1}^i = \bar{S}^i(\theta_l, \theta_{l+1})$ for all $l = 0, \dots, m$ as the cost can only decrease by this substitution and Claim 2 guarantees that the solution is still feasible i.e. $S_T^i \ge \dots \ge S_1^i \ge s_i$. It follows that one can build a path in $G$ of same cost as the corresponding solution to ($P_i$) using the node set $0 = \theta_0^i, \theta_1^i, \dots, \theta_m^i, \theta_{m+1}^i = T + 1$. Hence a shortest path from 0 to $T + 1$ in $G$ can be turned into an optimal solution to the original single-product problem ($P_i$).

## 2.3. The 2-approximation algorithm

We now show how to combine this solution with the decomposition, uncrossing and recombination technique presented in [11] for the deterministic JRP, to derive a 2-approximation algorithm.

**Decomposition** We consider a (multi-product) relaxation ($\tilde{P}_0$) of problem (P) in which for all $i \in [n]$, the individual fixed costs $K_i$ are set to zero and the holding costs and penalty costs are halved. We can again reduce problem ($\tilde{P}_0$) to a shortest path problem in the same graph $G$, except that the cost $\tilde{c}^0(u, v)$ of an arc $(u, v)$ is now defined as follows:

$$ \tilde{c}^0(u, v) = \begin{cases} \sum_{i \in [n]} \frac{1}{2} \lambda^i(s_i, u, v) & \text{if } u = 0 \\ \min_{S^1 \ge s_1, \dots, S^n \ge s_n} K_0 + \sum_{i \in [n]} \frac{1}{2} \lambda^i(S^i, u, v) & \text{otherwise} \end{cases} $$

It is clear that, for $u \ne 0$, $\tilde{c}^0(u, v)$ is separable and that the minimization can be performed independently over $S^1, \dots, S^n$, hence we have

$$ \tilde{c}^0(u, v) = K_0 \delta(u) + \sum_{i \in [n]} \frac{1}{2} \lambda^i(\bar{S}^i(u, v), u, v) \quad (6) $$

For all $i \in [n]$, we also consider single-product relaxations ($\tilde{P}_i$) of ($P_i$) where holding costs and penalty costs are halved by setting the cost of an arc $(u, v)$ to $\tilde{c}^i(u, v) = K_i \cdot \delta(u) + \frac{1}{2} \lambda^i(\bar{S}^i(u, v), u, v)$. Following the same discussion as above, the multi-product relaxation ($\tilde{P}_0$) and the single-product relaxations ($\tilde{P}_i$) can be solved through shortest path. Let

<footnote n="2">² Note that it uses the hypothesis that $\pi_t^i = \alpha^i \cdot h_t^i$.</footnote>

<footnote n="3">³ Observe that $\bar{S}^i(u, v)$ minimizes $\lambda^i(S^i, u, v)$ over $S_i \ge s_i$.</footnote>
us denote by OPT the optimal value to (P), by $\overline{OPT}_0$ the optimal value to $(\tilde{P}_0)$ and by $\overline{OPT}_i$ the optimal value to $(\tilde{P}_i)$. We have:

**Claim 3.** $OPT \ge \overline{OPT}_0 + \sum_{i \in [n]} \overline{OPT}_i$.

**Proof.** We can build feasible solutions to $(\tilde{P}_0)$ and all problems $(\tilde{P}_i)$ from a solution to OPT by simply ordering product $i$ exactly in the same periods as OPT, and in the same quantity, in all problems. As the holding and penalty costs are halved in $(\tilde{P}_0)$ and $(\tilde{P}_i)$, and since the joint set-up costs are accounted only in $(\tilde{P}_0)$ and the individual set-up costs are accounted only in $(\tilde{P}_i)$, the sum of the costs of the corresponding solutions is no more than OPT. The result follows by bounding from below the costs of these feasible solutions by $\overline{OPT}_0$ and $\overline{OPT}_i$, respectively. $\square$

**Uncrossing** Let $\tilde{\Theta}_0 = \{\theta_1^0, \dots, \theta_{m_0}^0\}$ and $\tilde{\Theta}_i = \{\theta_1^i, \dots, \theta_{m_i}^i\}$ be the (ordered) set of ordering periods in an optimal solution to problem $(\tilde{P}_0)$ and $(\tilde{P}_i)$, respectively. We let

$$R_i = U(\tilde{\Theta}_i, \tilde{\Theta}_0) = \{\rho_1^i, \dots, \rho_{k_i}^i\} \subseteq \tilde{\Theta}_0$$

where $U(\cdot, \cdot)$ is defined by (3). $R_i$ represents a (sub)set of candidate ordering periods for product $i$. For convenience, we define $\rho_0^i = 0$ and $\rho_{k_i+1}^i = T+1$.

**Recombination** We build a solution $\mathcal{R}$ to problem (P) by ordering $\tilde{S}^i(\rho_k^i, \rho_{k+1}^i) - \tilde{S}^i(\rho_{k-1}^i, \rho_k^i)$ units of product $i$ in periods $\rho_k^i$ for $k=1, \dots, k_i$. Notice that the resulting solution is feasible by Claim 2.

**Theorem 1.** The algorithm above is a 2-approximation for the SSJRP under the stochastic volume model introduced in [25]. That is, the expected cost of solution $\mathcal{R}$ is at most twice the expected cost of the optimal solution OPT.

**Proof.** It suffices to prove that the uncrossing and recombination operations of the algorithm produce a solution with a cost of at most twice $\overline{OPT}_0 + \sum_{i=1}^n \overline{OPT}_i$ in combination with Claim 3 to prove the result.

From Claim 1, we know that $R_i \subseteq \tilde{\Theta}_0$, and thus $|U_{i \in [n]} R_i| \le |\tilde{\Theta}_0|$, and $|R_i| \le 2 \cdot |\tilde{\Theta}_i|$ for all $i \in [n]$. Hence the sum of the setup costs incurred in $\mathcal{R}$ is at most twice the sum of the setup costs in the optimal solutions to $(\tilde{P}_0)$ and $(\tilde{P}_i)$ associated with $\tilde{\Theta}_0$ and $\tilde{\Theta}_i$.

Now regarding the holding and penalty costs, there are two types of situations to analyze from Claim 1 for any $i$ and any $k \in \{0, \dots, k_i\}$. First suppose that $\rho_k^i$ and $\rho_{k+1}^i$ coincide with two periods $\theta_i^0, \theta_{i+1}^0$ for some $l \in \{0, \dots, m_0\}$. Then the holding and penalty costs of product $i$ for solution $\mathcal{R}$ in problem (P), over the periods $[\rho_k^i, \rho_{k+1}^i)$, is at most twice those of product $i$ for the solution to sub-problem $(\tilde{P}_0)$ over the same periods (as the corresponding costs are halved in $(\tilde{P}_0)$). Assume now that $\rho_k^i$ and $\rho_{k+1}^i$ are contained in an interval $[\theta_i^i, \theta_{i+1}^i]$ for some $l \in \{0, \dots, m_i\}$. Defining $\phi_i^l(x) := h_i^l \int_{q=0}^x (x-q)f_i^l(q) dq + \pi_i^l \int_{q=x}^\infty (q-x)f_i^l(q) dq \ge 0$, we have in this case:

$$\begin{aligned} & \lambda^i(\tilde{S}^i(\theta_i^i, \theta_{i+1}^i), \theta_i^i, \theta_{i+1}^i) \\ &= \sum_{t=\theta_i^i}^{\theta_{i+1}^i-1} \phi_t^i(\tilde{S}^i(\theta_i^i, \theta_{i+1}^i)) \ge \sum_{t=\rho_k}^{\rho_{k+1}^i-1} \phi_t^i(\tilde{S}^i(\theta_i^i, \theta_{i+1}^i)) \\ &\ge \sum_{t=\rho_k}^{\rho_{k+1}^i-1} \phi_t^i(\tilde{S}^i(\rho_k^i, \rho_{k+1}^i)) = \lambda^i(\tilde{S}^i(\rho_k^i, \rho_{k+1}^i), \rho_k^i, \rho_{k+1}^i) \end{aligned}$$

The first inequality comes from the non-negativity of $\phi_i^l(\cdot)$, while the second comes from the fact that $\tilde{S}^i(u, v)$ is a minimizer of $\lambda^i(\tilde{S}^i(u, v), u, v)$ over $S^i \ge s_i$ for $u > 0$ and $\tilde{S}^i(u, v) = s_i$ if $u = 0$. The holding and penalty costs of product $i$ for solution $\mathcal{R}$ in problem (P) over the corresponding

interval is thus also at most twice the corresponding costs for the optimal solution of problem $(\tilde{P}_i)$. $\square$

We conclude this section with one final remark regarding the numerical efficiency of our approach in practice. To guarantee computational tractability in this setting, we assume access to efficient methods for the precomputing phase in Vargas's approach. Specifically, this means either having exact analytical formulae, tabular references with a known precision, or the ability to approximate arc costs to arbitrary precision using numerical methods. The accuracy of Vargas's approach (and our extension to the joint replenishment problem) directly depends on the fidelity of this precomputing step: if arc costs are computed within an error of at most $\epsilon$, then the deviation in the optimal solution is bounded by $n \cdot \epsilon$, where $n$ is the number of vertices in the graph.

# 3. Stochastic timing model

## 3.1. The model

Following [5,1], we now assume that the volume of the demands is deterministic but their timing is stochastic. More precisely for each product $i \in [n]$, we assume that there are $n_i \in \mathbb{N}$ demand requests to fulfill over the planning horizon. For a given request $j \in \{1, \dots, n_i\}$, the quantity $d_j^i > 0$ of the $j$th request is known but the period in which it occurs is random and follows a discrete probability distribution $p_j^i(\cdot)$ over a subset of periods $\{l_j^i, \dots, u_j^i\} \subseteq [T]$. All demands $d_j^i$ must be served in a single period (once fully available, no partial delivery allowed) and can be served from inventory (carried possibly from before $l_j^i$) or backlogged, but $u_j^i$ is a strict deadline that cannot be exceeded. In other words, any unit that satisfies demand $j$ must be ordered before $u_j^i$. Note that deterministic demands are a special case in which every demand request has a deterministic distribution over a single period $l_j^i = u_j^i$. For a discussion of the relevance of the corresponding model, see [1]. We assume here for simplicity that the initial stock $s_i$ are zero, although the approach can be extended seamlessly to any starting inventory/backlog value.

The expected holding and backlogging cost for producing one unit of product $i$ in period $t$ to satisfy demand $d_j^i$ is denoted by $EC_j^i(t)$, which is defined as follows for any $t \le u_j^i$:

$$EC_j^i(t) = \sum_{l=t}^{u_j^i} h_l^i \sum_{k=l+1}^{u_j^i} p_j^i(k) + \sum_{l=l_j^i}^{t-1} \pi_l^i \sum_{k=l_j^i}^{l} p_j^i(k) \quad (7)$$

By convention, we define $EC_j^i(t) = +\infty$ for all $t \in \{0\} \cup \{u_j^i + 1, \dots, T+1\}$. The goal is to find quantities $q_j^i(t)$ for each $i \in [n]$, $j \in [n_i]$ and $t \le u_j^i$ with $\sum_{i=1}^{u_j^i} q_j^i(t) = d_j^i$ that minimizes the total cost

$$\sum_{t \in [T]} \left[ K_0 \delta \left( \sum_{i \in [n]} q_j^i(t) \right) + \sum_{i \in [n]} K_i \delta \left( \sum_{j \in [n_i]} q_j^i(t) \right) + \sum_{i \in [n]} \sum_{j \in [n_i]} EC_j^i(t) q_j^i(t) \right]$$

We call this problem (P').

## 3.2. Solving the single item problem (in a similar spirit as in [1])

When there is no joint setup cost $K_0$ and only one product, say $i$, the single-product problem $(P'_i)$ consists in determining $q_j^i(t)$ for each $j \in [n_i]$ and $t \le u_j^i$, with $\sum_{i=1}^{u_j^i} q_j^i(t) = d_j^i$, that minimizes the cost

$$\sum_{t \in [T]} K_i \delta \left( \sum_{j \in [n_i]} q_j^i(t) \right) + \sum_{j \in [n_i]} EC_j^i(t) q_j^i(t)$$

For any $i \in [n]$ and $j \in [n_i]$, we define $\tilde{t}_j^i := \min\{t: 1 \le t \le T \text{ and } t = \arg \min_{t \in [T]} EC_j^i(t)\}$. We also call regeneration interval any interval of the
form $[t_1, t_2]$, such that orders are placed in periods $t_1 < t_2$ and no order occurs in periods $t$, for $t_1 < t < t_2$. Recall that we assume that the ratio between inventory holding and backlogging costs is time-independent ($\pi_i^t = \alpha^t h_i^t$ for all $t \in [T]$). Proposition 1, Theorem 3 and Proposition 2 in [1] derive useful dominance properties that allow to solve $(P'_i)$ in polynomial time through dynamic programming / shortest path. We present these results below for completeness.

**Proposition 1.** For all product $i \in [n]$ and all $j \in [n_i]$, we have $\vec{i}_j^i \in [l_j^i, u_j^i]$.

**Theorem 2.** $\forall t \in [T]$, $EC_j^t(u)$ is strictly decreasing on $\{1, \dots, \vec{i}_j^t\}$ and non-decreasing on $\{\vec{i}_j^t, \dots, u_j^t\}$.

**Proposition 2.** There is an optimal solution such that for any regeneration interval $[t_1, t_2]$ involved in the solution and for every $i \in [n]$ and $j \in [n_i]$ such that $t_1 \le \vec{i}_j^t \le t_2$, $d_j^t$ will be produced either in $t_1$ or in $t_2$.

Consider again the directed graph $\mathcal{G} = (\mathcal{V}, \mathcal{A})$ with node set $\mathcal{V} = \{0, 1, \dots, T+1\}$ (nodes 0 and $T+1$ are artificial nodes representing the beginning and the end of the horizon, respectively) and with arc set $\mathcal{A} = \{(u, v) \in \mathcal{V}^2: u < v\}$. Recall that for all $j \in [n_i]$ and $t > u_j^t$, $EC_j^t(u) = EC_j^t(0) = +\infty$. For an arc $(u, v)$ and $i \in [n]$, $j \in [n_i]$ such that $u \le \vec{i}_j^t < v$, we let

$$ \gamma_j^t(u, v) = \begin{cases} EC_j^t(u) \cdot d_j^i & \text{if } EC_j^t(u) < +\infty \text{ and } EC_j^t(u) \le EC_j^t(v) \\ EC_j^t(v) \cdot d_j^i & \text{otherwise} \end{cases} $$

Note that $\gamma_j^t(u, v) = +\infty$ if and only if $u = 0$ and $v > u_j^t$. We define in addition the cost of an arc $(u, v)$ as

$$ c^i(u, v) := K^i \delta(u) + \sum_{j \in [n_i]: u \le \vec{i}_j^t < v} \gamma_j^i(u, v) $$

A (finite cost) path from 0 to $T+1$ in $\mathcal{G}$ corresponds to a sequence of periods $\Theta_i = \{\theta_0^i, \theta_1^i, \dots, \theta_{m_i}^i, \theta_{m_i+1}^i\}$, with $0 = \theta_0^i < \theta_1^i < \dots < \theta_{m_i}^i < \theta_{m_i+1}^i = T+1$, that can be interpreted as a solution to problem $(P'_i)$ of same cost as follows. We consider the nodes $\theta_1^i, \dots, \theta_{m_i}^i$ in the path as the set of ordering periods. Since $\Theta_i$ partitions the planning horizon in disjoint intervals, there is a unique $l \in \{0, \dots, m_i\}$ such that $\theta_l^i \le \vec{i}_j^t < \theta_{l+1}^i$ for any given demand request $j \in [n_i]$. Demand $d_j^i$ is ordered in period $\theta_l^i$ if $\gamma_j^t(\theta_l^i, \theta_{l+1}^i) = EC_j^t(\theta_l^i) \cdot d_j^i$ and in period $\theta_{l+1}^i$ otherwise. It follows that any path with finite cost corresponds to a feasible solution to problem $(P'_i)$ of same cost. Conversely, any optimal solution to $(P'_i)$ with a sequence of ordering periods $\theta_1^i, \dots, \theta_{m_i}^i$ can be mapped to a path in $\mathcal{G}$ whose cost is no greater by serving demand request $d_j^i$ from the cheapest ordering period between $\theta_l^i$ and $\theta_{l+1}^i$ (from Proposition 2) where $l \in \{0, \dots, m_i\}$ is the unique value such that $\theta_l^i \le \vec{i}_j^t < \theta_{l+1}^i$. Hence, solving the shortest path problem from 0 to $T+1$ solves the single-item problem $(P'_i)$. (Note that since we assume $s_i = 0$, the first order happens in a period $\theta_l^i \le \min_{j \in [n_i]} u_j^t$.)

### 3.3. The 2-approximation algorithm

As in the stochastic demand volume case, it is possible to combine the solution of single-echelon problems with the decomposition, un-crossing and recombination technique [11] to derive a 2-approximation algorithm for $(P')$. We describe the approach in the remainder of this section.

**Decomposition** We again consider a (multi-product) relaxation $(\tilde{P}'_0)$ of problem $(P')$ with fixed costs $K_i$ set to zero and where the holding costs and penalty costs are halved. As in the stochastic volume case, the problem reduces to a shortest path problem in the same graph $\mathcal{G} = (\mathcal{V}, \mathcal{A})$, but the cost of an arc $(u, v)$ is now

$$ \tilde{c}^0(u, v) = K^0 \delta(u) + \frac{1}{2} \sum_{i \in [n]} \sum_{j \in [n_i]: u \le \vec{i}_j^t < v} \gamma_j^i(u, v) $$

We can also easily solve the single-product relaxation $(\tilde{P}'_i)$ of $(P'_i)$ where holding costs and penalty costs are halved by setting the cost of arc $(u, v)$ to $\tilde{c}^i(u, v) = K^i \delta(u) + \frac{1}{2} \sum_{j \in [n_i]: u \le \vec{i}_j^t < v} \gamma_j^i(u, v)$. Following the same discussion as above, the multi-product relaxation $(\tilde{P}'_0)$ and the single-product relaxations $(\tilde{P}'_i)$ can be solved through shortest path in $\mathcal{G}$, each with their respective arc costs. Let us denote by $OPT'$ the optimal value to $(P')$, by $OPT'_0$ the optimal value to $(\tilde{P}'_0)$ and by $OPT'_i$ the optimal value to $(\tilde{P}'_i)$. Following the same proof as for Claim 3, we have:

**Claim 4.** $OPT' \ge OPT'_0 + \sum_{i \in [n]} OPT'_i$.

**Uncrossing** Let $\tilde{\Theta}_0 = \{\theta_0^0, \dots, \theta_{m_0}^0\}$ and $\tilde{\Theta}_i = \{\theta_i^0, \dots, \theta_{m_i}^0\}$ be the (ordered) set of ordering periods in an optimal solution to problem $(\tilde{P}'_0)$ and $(\tilde{P}'_i)$, respectively. We let

$$ R_i = U(\tilde{\Theta}_i, \tilde{\Theta}_0) = \{\rho_i^0, \dots, \rho_{k_i}^0\} \subseteq \tilde{\Theta}_0 $$

where $U(\cdot, \cdot)$ is defined by (3). $R_i$ represents a (sub)set of candidate ordering periods for product $i$. For convenience, we define $\rho_i^0 = 0$ and $\rho_{k_i+1}^0 = T+1$.

**Recombination** We now build a solution $\mathcal{R}$ to problem $(P')$ by ordering in periods $\rho_k^i$ for $k=1, \dots, k_i$ as follows. For any $i \in [n]$ and $j \in [n_i]$, let $k \in \{0, \dots, k_i\}$ be the index such that $\rho_k^i \le \vec{i}_j^t < \rho_{k+1}^i$. By Claim 1, $\{\rho_k^i, \rho_{k+1}^i\}$ are either two consecutive ordering periods $\tilde{\Theta}_0$, i.e. there exists $l \in \{0, \dots, m_0\}$ such that $\theta_l^0 = \rho_k^i < \rho_{k+1}^i = \theta_{l+1}^0$, or are included between two consecutive orders of $\tilde{\Theta}_i$, i.e. $\theta_l^i \le \rho_k^i < \rho_{k+1}^i \le \theta_{l+1}^i$ for some $l \in \{0, \dots, m_i\}$. In the first case, we order demand $d_j^i$ as in the optimal solution to problem $(\tilde{P}'_0)$, i.e. in period $\theta_l^0$ or $\theta_{l+1}^0$. In the second case ($\theta_l^i \le \rho_k^i \le \vec{i}_j^t < \rho_{k+1}^i \le \theta_{l+1}^i$), we order demand $d_j^i$ in $\rho_k^i$ if $d_j^i$ was ordered in period $\theta_l^i$ in the optimal solution to problem $(\tilde{P}'_i)$, and in $\rho_{k+1}^i$ otherwise.

**Theorem 3.** The algorithm above is a 2-approximation for SSJRP under the stochastic demand timing model of [5,1]. That is, the expected cost of $\mathcal{R}$ is at most twice the expected cost of an optimal solution to $(P')$.

**Proof.** We need to prove that the expected cost of solution $\mathcal{R}$ in problem $(P')$ is at most twice $OPT'$. Using Claim 4, it is again enough to bound the cost of solution $\mathcal{R}$ by twice the value of $OPT'_0 + \sum_{i \in [n]} OPT'_i$.

From Claim 1, we know that $R_i \subseteq \tilde{\Theta}_0$, and thus $|\cup_{i \in [n]} R_i| \le |\tilde{\Theta}_0|$, and $|R_i| \le 2 \cdot |\tilde{\Theta}_i|$ for all $i \in [n]$. Hence the sum of the setup costs incurred in $\mathcal{R}$ is at most twice the sum of the setup costs in $(\tilde{P}'_0)$ and $(\tilde{P}'_i)$.

Now regarding the holding and backorder costs, consider $i \in [n]$, $j \in [n_i]$ and let $k \in \{0, \dots, k_i\}$ such that $\rho_k^i \le \vec{i}_j^t < \rho_{k+1}^i$. We distinguish the two different situations that are considered in the recombination phase of the algorithm (see again Claim 1). Assume first that $\{\rho_k^i, \rho_{k+1}^i\}$ coincide with two consecutive ordering periods in $\tilde{\Theta}_0$, i.e. $\theta_l^0 = \rho_k^i < \rho_{k+1}^i = \theta_{l+1}^0$ for some $l \in \{0, \dots, m_0\}$. In this case, the expected holding and backorder costs associated with demand $d_j^i$ are considered explicitly in the optimal solution associated with $\tilde{\Theta}_0$ for problem $(\tilde{P}'_0)$. Since the corresponding linear costs are halved in this subproblem, the final holding and backlogging costs incurred to serve demand $d_j^i$ in $\mathcal{R}$ at most double the ones accounted for the same demand request in $OPT'_0$. Assume now that $\{\rho_k^i, \rho_{k+1}^i\}$ is included into an interval defined by two consecutive orders of $\tilde{\Theta}_i$, i.e. $\theta_l^i \le \rho_k^i \le \vec{i}_j^t < \rho_{k+1}^i \le \theta_{l+1}^i$ for some $l \in \{0, \dots, m_i\}$. If $d_j^i$ is ordered in period $\theta_l^i$ in the optimal solution to $(\tilde{P}'_i)$, it is then ordered in period $\rho_k^i$ in the final solution $\mathcal{R}$. By
Theorem 2, $EC_j^i(\theta_i^i) \ge EC_j^i(\rho_k^i) \ge EC_j^i(\vec{r}_j^i)$. Conversely, if $d_j^i$ is ordered in period $\theta_{i+1}^i \le u_j^i$ in the optimal solution to $(\vec{P}'_i)$, then its ordering period in $\mathcal{R}$ is set to $\rho_{k+1}^i$ and, using again Theorem 2, we have that $EC_j^i(\theta_{i+1}^i) \ge EC_j^i(\rho_{k+1}^i) \ge EC_j^i(\vec{r}_j^i)$. In both cases, the expected holding and backorder costs incurred by demand $d_j^i$ is at most twice the holding and backorder costs of that same demand considered in $\overline{OPT}'_j$ (as the holding and backorder per-unit cost are halved in subproblem $(\vec{P}'_i)$). $\square$

## 4. Perspectives

We emphasize that we have kept models and settings simple for clarity, though it is likely that the mechanics illustrated here extend to richer SSJRP settings where we have for instance time-independent per-unit production/purchasing costs for each product, the joint ordering costs are time-dependent, or other natural variants already considered in [11] in the deterministic setting.

We also strongly believe that the techniques from [11] can also be used in more dynamic stochastic JRP settings. In fact, we have already used similar ideas to design efficient 2-competitive online algorithms for opportunistic maintenance problems emerging from humanitarian logistics [8].

### CRediT authorship contribution statement

**Guillaume Massonnet:** Writing – review & editing, Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Conceptualization.
**Gautier Stauffer:** Writing – review & editing, Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Conceptualization.

### References

[1] K. Akartunalı, S. Dauzère-Pérès, Dynamic lot sizing with stochastic demand timing, Eur. J. Oper. Res. 302 (1) (2022) 221–229.

[2] E. Arkin, D. Joneja, R. Roundy, Computational complexity of uncapacitated multi-echelon production planning problems, Oper. Res. Lett. 8 (2) (1989) 61–66.

[3] L.d.S.L. Bastos, M.L. Mendes, D.R.d.L. Nunes, A.C.S. Melo, M.P. Carneiro, A systematic literature review on the joint replenishment problem solutions: 2006-2015, Production 27 (2017) e20162229.

[4] J.H. Bookbinder, J.-Y. Tan, Strategies for the probabilistic lot-sizing problem with service-level constraints, Manag. Sci. 34 (9) (1988) 1096–1108.

[5] M.C. Burstein, C.H. Nevison, R.C. Carlson, Dynamic lot-sizing when demand timing is uncertain, Oper. Res. 32 (2) (1984) 362–379.

[6] T. Cohen-Hillel, L. Yedidsion, The periodic joint replenishment problem is strongly NP-hard, Math. Oper. Res. 43 (4) (2018) 1269–1289.

[7] T. de Kok, C. Grob, M. Laumanns, S. Minner, J. Rambau, K. Schade, A typology and literature review on stochastic multi-echelon inventory models, Eur. J. Oper. Res. 269 (3) (2018) 955–983.

[8] B. Detienne, M. Gaury, G. Stauffer, An efficient 2-competitive online algorithm for kit update at MSF logistique, in: Journées Polyèdres et Optimisation Combinatoire - JPOC'13, Clermont-Ferrand, 2023.

[9] A. Federgruen, M. Tzur, The joint replenishment problem with time-varying costs and demands: efficient, asymptotic and $\epsilon$-optimal solutions, Oper. Res. 42 (6) (1994) 1067–1086.

[10] J.P. Gayon, G. Massonnet, C. Rapine, G. Stauffer, Constant approximation algorithms for the one warehouse multiple retailers problem with backlog or lost-sales, Eur. J. Oper. Res. 250 (1) (2016) 155–163.

[11] J.-P. Gayon, G. Massonnet, C. Rapine, G. Stauffer, Fast approximation algorithms for the one-warehouse multi-retailer problem under general cost structures and capacity constraints, Math. Oper. Res. 42 (3) (2017) 854–875.

[12] M. Khouja, S. Goyal, A review of the joint replenishment problem literature: 1989–2005, Eur. J. Oper. Res. 186 (1) (2008) 1–16.

[13] R. Levi, G. Janakiraman, M. Nagarajan, A 2-approximation algorithm for stochastic inventory control models with lost sales, Math. Oper. Res. 33 (2) (2008) 351–374.

[14] R. Levi, M. Pál, R.O. Roundy, D.B. Shmoys, Approximation algorithms for stochastic inventory control models, Math. Oper. Res. 32 (2) (2007) 284–302.

[15] R. Levi, R. Roundy, D. Shmoys, M. Sviridenko, A constant approximation algorithm for the one-warehouse multiretailer problem, Manag. Sci. 54 (4) (2008) 763–776.

[16] R. Levi, R. Roundy, V.A. Truong, X. Wang, Provably near-optimal balancing policies for multi-echelon stochastic inventory control models, Math. Oper. Res. 42 (1) (2016) 256–276.

[17] R. Levi, C. Shi, Approximation algorithms for the stochastic lot-sizing problem with order lead times, Oper. Res. 61 (3) (2013) 593–602.

[18] J. Muckstadt, A. Sapra, Principles of Inventory Management: When You Are Down to Four, Order More, 2009.

[19] T. Nonner, M. Sviridenko, An efficient polynomial-time approximation scheme for the joint replenishment problem, in: M. Goemans, J. Correa (Eds.), Integer Programming and Combinatorial Optimization, Springer, Berlin, Heidelberg, 2013, pp. 314–323.

[20] L. Peng, L. Wang, S. Wang, A review of the joint replenishment problem from 2006 to 2022, Manag. Syst. Eng. 1 (1) (2022) 9.

[21] E.P. Robinson, L.-L. Gao, A dual ascent procedure for multiproduct dynamic demand coordinated replenishment with backlogging, Manag. Sci. 42 (11) (1996) 1556–1564.

[22] E.P. Robinson, A. Narayanan, L.-L. Gao, Effective heuristics for the dynamic demand joint replenishment problem, J. Oper. Res. Soc. 58 (6) (2007) 808–815.

[23] O. Solyali, H. Süral, The one-warehouse multi-retailer problem: reformulation, classification, and computational results, Ann. Oper. Res. 196 (1) (2012) 517–541.

[24] G. Stauffer, Approximation algorithms for k-echelon extensions of the one warehouse multi-retailer problem, Math. Methods Oper. Res. 88 (3) (2018) 445–473.

[25] V. Vargas, An optimal solution for the stochastic version of the Wagner–Whitin dynamic lot-size model, Eur. J. Oper. Res. 198 (2) (2009) 447–451.

[26] P.H. Zipkin, Foundations of Inventory Management, 1 edition, McGraw-Hill/Irwin, 2000.