Ana Bušić
email: ana.busic@inria.fr

Alan Scheller

Comparing Markov chains: Aggregation and precedence relations applied to sets of states, with applications to assemble-to-order systems

Keywords: Markov chains, precedence relations, bounding models, aggregation, supermodularity, assemble-to-order (ATO) systems Mathematics Subject Classification (MSC2010): 60E15, 60J10, Secondary: 90B25, 90B05

Solving Markov chains is in general difficult if the state space of the chain is very large (or infinite) and lacking a simple repeating structure. One alternative to solving such chains is to construct models that are simple to analyze and provide bounds for a reward function of interest. We present a new bounding method for Markov chains inspired by Markov reward theory: Our method constructs bounds by redirecting selected sets of transitions, facilitating an intuitive interpretation of the modifications of the original system. We show that our method is compatible with strong aggregation of Markov chains; thus we can obtain bounds for an initial chain by analyzing a much smaller chain. We illustrate our method by using it to prove monotonicity results and bounds for assemble-to-order systems.

1 Introduction.

In Markov chain modeling, one often faces the problem of combinatorial state space explosion: modeling a system completely requires an unmanageable -combinatorial -number of states. Many high-level formalisms, such as queueing networks or stochastic Petri nets, have been developed to simplify the specification and storage of the Markov chain. However, these models rarely have closed-form solutions, thus one must resort to numerical methods. Unfortunately, numerical methods are inefficient when the size of the state space becomes very large or for models with infinite state space that do not exhibit a special repeating structure that admits a matrix analytic approach (Neuts [31]). Typically, the latter approach is quite limited if the state space is infinite in more than one dimension. An alternative approach to cope with state space explosion is to design new models that (i) provide bounds for a specific measure of interest (for instance the probability of a failure in a complex system); and (ii) are simpler to analyze than the original system.

Establishing point (i): the relationship between the original and new (bounding) systems may be based on different arguments. Potentially the most general way of obtaining bounds is by stochastic comparison, which gives bounds for a whole family of reward functions (for instance increasing or convex functions). Furthermore, stochastic comparison provides bounds for both the steady-state and transient behavior of the studied model. Many results have been obtained using strong stochastic order (i.e. generated by increasing functions) and coupling arguments (Lindvall [START_REF] Lindvall | Lectures on the coupling method[END_REF]). Recently, an algorithmic approach has been proposed (Fourneau and Pekergin [START_REF] Fourneau | An algorithmic approach to stochastic bounds[END_REF]) to construct stochastic bounds, based on stochastic monotonicity; this stochastic monotonicity provides simple algebraic sufficient conditions for the stochastic comparison of Markov chains. Ben Mamoun et al. [START_REF] Ben Mamoun | Increasing convex monotone Markov chains: Theory, algorithm and applications[END_REF] showed that an algorithmic approach is also possible using increasing convex ordering that allows one to compare variability. The clear advantage of stochastic comparison is its generality: it provides bounds for a whole family of rewards, both for the steady-state and transient behavior of the studied system. Its drawback is that, due to its generality, it requires strong constraints that may not apply to the system of interest. For more details on the theoretical aspects of stochastic comparison we refer the reader to Müller and Stoyan [START_REF] Müller | Comparison methods for stochastic models and risks[END_REF], and Shaked and Shantikumar [START_REF] Shaked | Stochastic orders and their applications[END_REF].

For this reason more specialized methods than stochastic comparison have also been developed, which apply only to one specific function, and only in the steady-state. Muntz et al. [START_REF] Muntz | Bounding availability of repairable computer systems[END_REF] proposed an algebraic approach to derive bounds of steady-state rewards without computing the steady-state distribution of the chain, founded on results of Courtois and Semal [START_REF] Courtois | Bounds for the positive eigenvectors of nonnegative matrices and for their approximations by decomposition[END_REF][START_REF]On polyhedra of Perron-Frobenius eigenvectors[END_REF] on eigenvectors of nonnegative matrices. This approach was specially designed for reliability analysis of highly reliable systems, and requires special constraints on the structure of the underlying Markov chain. This approach was further improved and generalized by various authors (Lui and Muntz [START_REF] Lui | Computing bounds on steady state availability of repairable computer systems[END_REF], Semal [START_REF]Refinable bounds for large Markov chains[END_REF], Carrasco [START_REF] Carrasco | Bounding steady-state availability models with group repair and phase type repair distributions[END_REF], Mahevas and Rubino [START_REF] Mahevas | Bound computation of dependability and performance measures[END_REF]), but the primary assumption for its applicability is still that there is a very small portion of the state space that has a very high probability of occurrence, while the other states are almost never visited.

Similarly, Van Dijk [START_REF] Van Dijk | Bounds and error bounds for queueuing networks[END_REF] (see also Van Dijk and Van der Wal [START_REF] Van Dijk | Simple bounds and monotonicity results for finite multiserver exponential tandem queues[END_REF]) proposed a different method for comparing two chains in terms of a particular reward function, often referred to as the Markov reward approach. This method allows the comparison of mean cumulated and stationary rewards for two given chains. A simplified version of the Markov reward approach, called the precedence relation method, was proposed by Van Houtum et al. [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF]. The origin of the precedence relation method dates back to Adan et al. [START_REF] Adan | Upper and lower bounds for the waiting time in the symmetric shortest queue system[END_REF], and it has been successfully applied to various problems (Van Houtum et al. [START_REF] Van Houtum | Performance analysis of parallel identical machines with a generalized shortest queue arrival mechanism[END_REF], Tandra et al. [START_REF] Tandra | Join minimum cost queue for multiclass customers: Stability and performance bounds[END_REF], Leemans [START_REF] Leemans | Provable bounds for the mean queue lengths in a heterogeneous priority queue[END_REF]). The advantage of this method is its straightforward description of the modifications of the initial model.

The precedence relation method consists of two steps. Precedence relations are first established on the states of the system, based on the reward function (or family of functions) one wants to study. Then an upper (resp. lower) bound for the initial model can be obtained simply by redirecting the transitions to greater (resp. smaller) states with respect to the precedence relations established in the first step. This can also be thought of as replacing transitions in the original model with ones to greater (resp. smaller) states. A significant drawback of the precedence relation method is that transitions can be redirected only to greater states (in the case of an upper bound) or only to smaller states (lower bound): This method does not allow one to obtain bounding models by changing the probability of transitions, or by redirecting some transitions to greater states and others to smaller states, which is typically needed, for example, if one wants to keep the mean behavior of a part of the system (for instance arrivals to a queue), but change its variability. One small example of such a system is given in Section 3.

We propose a generalization of precedence relations to sets of states. By establishing precedence relations on sets of states, instead of individual states, one can avoid the above problems as long as the new destination states considered together are greater (resp. smaller) than the original set of destinations. As the modification of the probability of a transition can also be seen as replacement of one transition by two new transitions, one of which is a loop, this significantly increases the applicability of the precedence relation method, by allowing one to compare systems with different variability.

We now discuss point (ii): how to derive models that are simpler to solve. In the context of stochastic comparison, different types of bounding models have been used: models having closed form solutions, models that are easier to analyze using numerical methods, and aggregation (Fourneau and Pekergin [START_REF] Fourneau | An algorithmic approach to stochastic bounds[END_REF], Fourneau et al. [START_REF] Fourneau | Algorithms for an irreducible and lumpable strong stochastic bound[END_REF]). To our knowledge, the precedence relation method has been combined only with the first two simplifications. We show here that it is also compatible with aggregation. Thus we prove the validity of applying the precedence relation method on sets of states, in concert with the simplifying technique of aggregation.

We illustrate our new technique to prove new results for Assemble-to-Order (ATO) systems. To facilitate this, we first present a general framework for ATO systems, and then a detailed description of a class within this framework that we focus on: ATO models with Partial Order Service (ATO-POS). We then prove a monotonicity property with respect to supermodularity for ATO-POS systems, also discussing how this monotonicity extends, or fails to extend, to other classes within our ATO framework. We then apply precedence relations on sets of states to two specific problems in the ATO-POS class: ATO-POS models with individual returns (see [START_REF] Dayanik | The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-to-order systems[END_REF] for example), and the service tools problem. This latter problem, introduced by Vliegen and Van Houtum [45], models a single location multi-item inventory system in which customers demand different sets of service tools, needed for a specific maintenance action, from a stock point.

Considering the ATO models with POS, we show that in addition to deriving bounding models, our new method can be used to establish monotonicity results with respect to model parameters. We illustrate this in Section 6.1 to prove monotonicity of order fill rates with respect to coupling in arrivals for ATO-POS systems with state dependent replenishment times. This extends results in Dayanik et al. [START_REF] Dayanik | The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-to-order systems[END_REF] to the case of uncapacitated ATO systems. Then in 6.2 we apply our technique on the service tool problem, proving that the approximation proposed by Vliegen and Van Houtum [45] provides a bound for their original model. Our bounding model has a state space that is highly reduced compared to Vliegen and Van Houtum's original system: its dimension is equal to the number of different types of tools (I), while the original model has dimension 2 I . This paper is organized as follows. In Section 2 we give an overview of the precedence relation method proposed by Van Houtum et al. [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF]. In Section 3 we show the limits of this method and we propose and prove the validity of our generalization. Section 4 is devoted to aggregation and its connections with our method. In Section 5 we give an overview of known results for our general class of ATO systems and we establish a monotonicity property for supermodular rewards for a class of state dependent ATO models. This property is then used in Section 6 to obtain our results for ATO-POS and service tools systems. Finally, Appendix A contains a supermodularity characterization on a discrete lattice, necessary for our proofs in Sections 5 and 6.

Precedence relations.

Let {X n } n≥0 be a discrete time Markov chain (DTMC) on a countable state space X and denote by P its transition matrix. A probability measure π on X is stationary if πP = π. It is attractive if for any probability measure ν on X , the sequence of Cesàro averages of νP n converges weakly to π, i.e.

π(x) = lim t→∞ 1 t t-1 n=0 (νP n )(x), x ∈ X .
A DTMC is said to be stable if it has a unique and attractive stationary probability measure. This implies in particular that the graph of the Markov chain has a unique strongly connected component (recurrent class) that will be reached with probability one from any other state. We consider in the following only stable DTMC.

For a given reward (or cost) function r : X → R, the mean stationary reward is given by:

a = x∈X r(x)π(x).
Directly computing the stationary distribution π is often difficult if, for instance, the state space is infinite in many dimensions or finite, but prohibitively large. The main idea of the precedence relation method proposed by Van Houtum et al. [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF] is to obtain upper or lower bounds for a without explicitly computing π. By redirecting selected transitions of the original model, the graph of the chain is modified to obtain a new chain that is significantly easier to analyze. For example, one might truncate the chain by blocking the outgoing transitions from a subset of states.

Some special care needs to be taken in order to ensure that the reward function of the new chain provides a bound on the reward function of the original chain. We denote by v t (x) (resp. by v t (x)) the expected cumulated reward over the first t periods for the original (resp. modified) chain when starting in state x ∈ X :

v t+1 (x) = r(x) + y∈X P [x, y]v t (y), t ≥ 0, (1)
where v 0 (x) := 0, ∀x ∈ X . Then for any state x that belongs to the unique recurrent class, the mean stationary reward satisfies:

a = lim t→∞ v t (x) t .
If we can show that

v t (x) ≤ v t (x), ∀x, ∀t ≥ 0, (2)
then we have also the comparison of mean stationary rewards:

a = lim t→∞ v t (x) t ≤ lim t→∞ v t (x ) t = a,
where x and x belong respectively to the unique recurrent class of the initial and the modified chain.

In order to establish (2), a precedence relation is defined on state space X as follows:

x y if v t (x) ≤ v t (y), ∀t ≥ 0.
When talking about rewards, we will often say that a state x is less attractive than y (for the reward function r) if x y.

The following theorem states that redirecting (i.e. replacing with) transitions to less (more) attractive states results in an lower (upper) bound for mean stationary reward (Van Houtum et al. [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF]Theorem 4.1]): Theorem 1. Let {X n } be a DTMC and let {Y n } be a chain obtained from {X n } by replacing some transitions (x, y) with transitions (x, y ) such that y y . Then:

v t (x) ≤ v t (x), ∀x, ∀t ≥ 0.
If both chains are stable, then the mean stationary rewards satisfy a ≤ a.

The above theorem allows one to easily construct bounding models by redirecting possibly only a few transitions. Van Houtum et al. [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF] illustrated their approach on the example of a system with the Join the Shortest Queue routing. In the following section we illustrate some of the limits of the precedence relation approach before proposing its generalization.

3 Precedence relations on sets of states.

The precedence relation method allows one to redirect transitions: the destination of the transition is modified, while its probability remains the same. Redirecting transitions to less (more) attractive states results in an lower (upper) bound. However, the following simple example shows that we cannot use the precedence relation method to compare models with the same average arrival rate, but different variabilities. In that case, the modified model is obtained by redirecting simultaneously some transitions to more and some to less attractive states, so Theorem 1 no longer applies.

Example 1 (Single queue with batch arrivals). We consider a single queue with two types of jobs:

• Class 1 jobs arrive individually following a Poisson process with rate λ 1 .

• Class 2 jobs arrive by batches of size 2, following a Poisson process with rate λ 2 .

We assume a single exponential server with rate µ, and let x denote the number of jobs in the system. Then the following events can occur in the system: event rate transition condition type 1 arrival

λ 1 x + 1 - type 2 arrival λ 2 x + 2 - service µ x -1 x > 0
Without loss of generality, we assume that λ 1 + λ 2 + µ = 1. Thus we can consider λ 1 , λ 2 and µ as the probabilities of the events in the corresponding discrete time (uniformized) chain. Suppose that we are interested in the mean number of jobs. The appropriate reward function is thus r(x) = x, ∀x. The corresponding t-period rewards satisfy:

v t+1 (x) = r(x) + λ 1 v t (x + 1) + λ 2 v t (x + 2) + µv t (x -1)1 {x>0} + µv t (x)1 {x=0} , x ≥ 0, t ≥ 0,
with v 0 (x) := 0, ∀x ≥ 0. Denote respectively by A 1 , A 2 and S the t-period rewards in new states after an arrival of type 1, an arrival of type 2 and a service in state x:

• A 1 (x, t) = v t (x + 1), • A 2 (x, t) = v t (x + 2), • S(x, t) = v t (x -1)1 {x>0} + v t (x)1 {x=0} .
Then:

v t+1 (x) = r(x) + λ 1 A 1 (x, t) + λ 2 A 2 (x, t) + µS(x, t), x ≥ 0, t ≥ 0. (3)
For all t ≥ 0:

v t (x) ≤ v t (x + 1), x ≥ 0, (4)
i.e. x x + 1, x ≥ 0. This can be easily shown by induction on t. Assume (4) holds for a given t ≥ 0 (for t = 0, v 0 := 0 is obviously non-decreasing). We can consider separately one period rewards, arrivals of each type, and service.

• One period rewards. v 1 = r is obviously non-decreasing.

• Type 1 arrivals.

For x ≥ 0, A 1 (x + 1, t) -A 1 (x, t) = v t (x + 2) -v t (x + 1) ≥ 0, since v t is non-decreasing. • Type 2 arrivals. For x ≥ 0, A 2 (x + 1, t) -A 2 (x, t) = v t (x + 3) -v t (x + 2) ≥ 0, since v t is non-decreasing. • Service. For x ≥ 0, S(x + 1, t) -S(x, t) = v t (x) -v t (x -1)1 {x>0} -v t (x)1 {x=0} = (v t (x) - v t (x -1))1 {x>0} ≥ 0, since v t is non-decreasing.
It follows from relation (3) that v t+1 is also non-decreasing, so (4) holds for all t ≥ 0. Now, suppose some class 1 jobs become class 2 jobs, keeping the total arrival rate constant. This means that these jobs arrive less often (only half of the previous rate), but they arrive in batches of size 2 (Figure 1). Then:

λ 1 = λ 1 -and λ 2 = λ 2 + 2 ,
where 0 < ≤ λ 1 . The total arrival rate is the same in both models, but the arrival process of the second system is more variable. Different transitions for both models are shown in Figure 2. Note that a part of the transition rate that corresponds to the arrivals of type 1 is replaced by a new transition that corresponds to the arrivals of type 2, but the rate is divided by two. This can be also seen as replacing one transition with rate by two transitions, each with rate /2: a transition to state x + 2 and a loop (see Figure 2). Thus we cannot apply Theorem 1, since it allows neither replacing only a part of a transition, nor replacing one transition with two new ones: one going to a more and one to a less attractive state. In the following we propose a more general method that allows us to replace a transition, or more generally a set of transitions, by another set of transitions having the same aggregate rate. Also, only a portion of some transitions may be redirected.

Generalization of precedence relations.

To aid intuition, we will introduce the main ideas by considering a single state x ∈ X . (The general result will be proved later, in Theorem 2.) Assume we want to replace (redirect) the outgoing transitions from x to a subset A with transitions to another subset B. For instance, in Figure 3 we want to replace transitions to A = {a 1 , a 2 } (on the left) with transitions to B = {b 1 , b 2 , b 3 } (on the right). We might also have some transitions from x to states that are not in A and that we do not want to redirect (transitions to states u and v in Figure 3).

Furthermore, we might want to redirect transitions only partially: in Figure 4 only the half of probability of transitions to A is replaced by transitions to B. Thus in order to describe the redirection of a set of transitions we will need to provide:

• the set A (resp. B) and the probabilities of transitions to each state in A (resp. B); • the weight factor ∆ (the amount of each transition to be redirected; the same scalar ∆ is applied to all transitions to states in A). Because we can include only part of a transition in set A, this common ∆ is not restrictive. Information on sets A and B and the corresponding transition probabilities will be given by two vectors α = (α(z)) z∈X and β = (β(z)) z∈X . Since the information on the amount to be redirected will be given by a weight factor, we only need the relative probabilities of transitions to the respective sets A and B. Therefore it is convenient to renormalize the vectors α and β. The modifications in Figures 3 and4 can now be described by vectors α = 0.5δ a1 + 0.5δ a2 and β = 0.2δ b1 + 0.2δ b2 + 0.6δ b3 , where for any y ∈ X , δ y denotes the Dirac measure in y:

δ y (z) = 1, z = y, 0, z = y , z ∈ X .
The weight factor ∆ is equal to 0.5 in Figure 3 and to 0.25 in Figure 4. We now formally generalize the previous example. Let α and β be two stochastic vectors: α(z) ≥ 0, β(z) ≥ 0, ∀z ∈ X and ||α|| 1 = ||β|| 1 = 1 (where ||α|| 1 := z∈X α(z) is the usual 1norm). Let {X n } be a DTMC with transition probability matrix P and t-period reward functions v t , satisfying the following relation:

z∈X α(z)v t (z) ≤ z∈X β(z)v t (z), ∀t ≥ 0. ( 5
)
Let A and B denote the supports of vectors α and β respectively:

A = supp(α) = {z ∈ X : α(z) > 0}, B = supp(β) = {z ∈ X : β(z) > 0}.
For the example in Figures 3 and4, A = {a 1 , a 2 } and B = {b 1 , b 2 , b 3 }. If (5) holds, we will say that the set of states A is less attractive than the set B with respect to probability vectors α and β, and we will denote this:

A α,β B.
We will show that if A α,β B, replacing the outgoing transitions to A (with probabilities α) by the outgoing transitions to B (with probabilities β) leads to an upper bound for t-period rewards (and thus also for the mean stationary reward, when it exists). Before giving this result in Theorem 2, note that relation ( 5) is indeed a generalization of precedence relations of states: Remark 1. Suppose x y, for some x, y ∈ X . Set α = δ x and β = δ y . Then (5) becomes:

v t (x) ≤ v t (y), t ≥ 0,
which is equivalent to x y by definition.

To see that [START_REF] Carrasco | Bounding steady-state availability models with group repair and phase type repair distributions[END_REF] indeed is more general than the precedence relation method (Van Houtum et al. [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF]), let α = δ x and β = 1 2 δ y + 1 2 δ z , x, y, z ∈ X . Then (5) becomes:

v t (x) ≤ 1 2 v t (y) + 1 2 v t (z), t ≥ 0.
We can write this {x} δx, 1 2 δy+ 1 2 δz {y, z}. By taking y = x + 1 and z = x -1 this is exactly the relation we need to prove in Example 1. The proof of this relation for Example 1 will be given in Section 3.2.

Replacing the outgoing transitions that correspond to the set A and probabilities α by the transitions that correspond to the set B and probabilities β (called (α, β)-redirection in the following), can be also represented in matrix form. The matrix T α,β (x) defined as:

T α,β (x)[w, z] = β(z) -α(z), w = x, 0, w = x,
describes the (α, β)-redirection of the outgoing transitions from state x ∈ X . The transition matrix of the modified chain after (α, β)-redirection of the outgoing transitions from state x ∈ X is then given by:

P = P + ∆ α,β (x)T α,β (x),
with the weight factor ∆ α,β (x), 0 ≤ ∆ α,β (x) ≤ 1. (Note that if ∆ α,β (x) = 0, we do not modify the chain.) In order for P to be a stochastic matrix, the weight factor ∆ α,β (x) must satisfy:

0 ≤ P [x, y] + ∆ α,β (x)(β(y) -α(y)) ≤ 1, y ∈ X , (6)
which can be also written as:

∆ α,β (x) ≤ min min y : α(y)>β(y) P [x, y] α(y) -β(y)
, min y : α(y)<β(y)

1 -P [x, y] β(y) -α(y) .
Without loss of generality, we can assume that the supports of α and β are disjoint: A ∩ B = ∅. Indeed, if there exists y ∈ X such that α(y) > 0 and β(y) > 0, then we can define new vectors α = 1 1-c (α -cδ y ) and β = 1 1-c (β -cδ y ), where c = min{α(y), β(y)}. Relation ( 5) is then equivalent to:

z∈X α (z)v t (z) ≤ z∈X β (z)v t (z), t ≥ 0.
Assuming A and B are disjoint, relation [START_REF] Courtois | Bounds for the positive eigenvectors of nonnegative matrices and for their approximations by decomposition[END_REF] has an intuitive interpretation given as Proposition 1: one can only redirect existing transitions. Proposition 1. For vectors α and β with supports A ∩ B = ∅ condition (6) is equivalent to:

α(y)∆ α, β (x) ≤ P [x, y], y ∈ X . (7)
Proof. Relation (7) follows trivially from (6) as α(y) > 0 and A ∩ B = ∅ imply β(y) = 0. In order to see that ( 7) implies ( 6), we will consider the following cases for an arbitrary y ∈ X :

• α(y) > 0. Then β(y) = 0 so relation [START_REF] Courtois | Bounds for the positive eigenvectors of nonnegative matrices and for their approximations by decomposition[END_REF] becomes: 0

≤ P [x, y] -∆ α,β (x)α(y) ≤ 1.
As we assumed that 0 ≤ ∆ α,β (x) ≤ 1, and 0 ≤ α(y) ≤ 1, the right inequality is trivial and the left one is simply relation [START_REF]On polyhedra of Perron-Frobenius eigenvectors[END_REF].

• β(y) > 0. Then α(y) = 0 so relation [START_REF] Courtois | Bounds for the positive eigenvectors of nonnegative matrices and for their approximations by decomposition[END_REF] becomes: 0

≤ P [x, y] + ∆ α,β (x)β(y) ≤ 1.
The left inequality is trivial. For the right one we have, using the fact that P is a stochastic matrix and β(y) ≤ 1:

P [x, y] + ∆ α,β (x)β(y) ≤ 1 - z =y P [x, z] + ∆ α,β (x) ≤ 1 - z =y ∆ α,β (x)α(z) + ∆ α,β (x) = 1 + ∆ α,β (x)   1 - z =y α(z)   = 1 + ∆ α,β (x)α(y) = 1,
where the second inequality follows from [START_REF]On polyhedra of Perron-Frobenius eigenvectors[END_REF].

• Finally the case α(y) = β(y) = 0 is trivial.

Until now we have considered only one state x and only one relation (α, β). Typically, we will redirect outgoing transitions for a subset of the state space, and we may need more than one relation. Let R be a set of relations that are valid for our model (i.e. they satisfy ( 5)). We will denote by R x ⊂ R the set of all relations that will be used for a state x ∈ X . (If the outgoing transitions for a state x ∈ X do not change, we will set R x := ∅.) Note that the same relation (α, β) could be applied to different states x and x (i.e. (α, β) ∈ R x ∩ R x ). In that case the corresponding weight factors ∆ (α,β) (x) and ∆ (α,β) (x ) need not to be equal. Also, there could be different relations (α, β) and (α , β ) that have the same supports A and B; we might even have that A (α,β) B, but B (α ,β ) A (if supp(α) = supp(β ) = A and supp(β) = supp(α ) = B). Thus our method can be made arbitrarily general. The following theorem states that under conditions similar to (5) and ( 6) (for all states and for a family of precedence relations), the t-period rewards v t of the modified chain satisfy:

v t (x) ≤ v t (x), x ∈ X , t ≥ 0.
Theorem 2. Let {X n } be a DTMC with transition probability matrix P and a reward r that is bounded from below:

∃m ∈ R, r(x) ≥ m, ∀x ∈ X . (8)
Denote by v t , t ≥ 0, the corresponding t-period rewards. Let R be a set of pairs of stochastic vectors such that for all (α, β) ∈ R:

y∈X α(y)v t (y) ≤ y∈X β(y)v t (y), t ≥ 0. ( 9
)
Let R x ⊂ R, x ∈ X (the precedence relations that will be applied to a state x ∈ X ). Let {Y n } be a DTMC with transition probability matrix P given by:

P = P + x∈X (α,β)∈Rx ∆ α,β (x)T α,β (x),
where the factors ∆ α,β (x), x ∈ X , (α, β) ∈ R x , satisfy:

0 ≤ P [x, y] + (α,β)∈Rx ∆ α,β (x)(β(y) -α(y)) ≤ 1, x, y ∈ X ( 10
) (i.e. 0 ≤ P [x, y] ≤ 1, x, y ∈ X ).
Then the t-period rewards v t of the modified chain satisfy:

v t (x) ≤ v t (x), x ∈ X , t ≥ 0. ( 11
)
Symmetrically, if ( 9) is replaced by:

y∈X α(y)v t (y) ≥ y∈X β(y)v t (y), t ≥ 0, (12)
for all (α, β) ∈ R, then the t-period rewards v t of the modified chain satisfy:

v t (x) ≥ v t (x), x ∈ X , t ≥ 0. ( 13
)
Proof. We will prove [START_REF] Fourneau | An algorithmic approach to stochastic bounds[END_REF] by induction on t. For t = 0 we have v 0 (x) = v 0 (x) := 0, x ∈ X , so ( 11) is trivially satisfied. Suppose ( 11) is satisfied for t ≥ 0. Then for t + 1 we have:

v t+1 (x) = r(x) + y∈X P [x, y] v t (y) ≥ r(x) + y∈X P [x, y]v t (y) = r(x) + y∈X   P [x, y] + (α,β)∈Rx ∆ α,β (x)T α,β (x)[x, y]   v t (y).
Relation [START_REF] Fourneau | Algorithms for an irreducible and lumpable strong stochastic bound[END_REF] implies that for all y ∈ X the series absolutely converges (in R + ∪ {+∞}) if v t is bounded from below. (Note that if r is bounded from below then v t is bounded from below for each t.) Thus, simplifying and interchanging summation:

v t+1 (x) ≥ v t+1 (x) + (α,β)∈Rx ∆ α,β (x) y∈X (β(y) -α(y)) v t (y) ≥ v t+1 (x). Remark 2. If for all x ∈ X |R x | < ∞ and for all (α, β) ∈ R x , |supp(α) ∪ supp(β)| < ∞, (14)
then the result of Theorem 2 holds without condition [START_REF] Dayanik | The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-to-order systems[END_REF]. In other words, if for all x ∈ X , the initial and the modified chain differ only in a finite number of the outgoing transitions, then the reward function r can be unbounded.

Corollary 1. Under the same conditions as in Theorem 2, if both the initial and the modified chain are stable, then the mean stationary reward of the modified chain {Y n } is an upper bound (a lower bound in case of ( 12)) for the mean stationary reward of the original chain {X n }.

Proving the relations.

Thus the steps in order to prove a bound are to first identify the set R, and then to prove the corresponding relations for the t-period rewards. We will illustrate this steps on our simple example of a queue with batch arrivals, discussed in Example 1.

Example 2. Consider again the two models from Example 1. We will show that:

v t (x) ≤ v t (x), x ≥ 0, t ≥ 0, (15)
where v t denotes the t-period rewards for the original chain and v t for the modified chain for reward function r(x) = x, ∀x. For each x ≥ 0, we want to replace a part of the transition that goes to x + 1 by two new transitions that go to x and x + 2. We will define vectors α x and β x as follows:

α x = δ x+1 , β x = 1 2 (δ x + δ x+2 ) . Let R x = {(α x , β x )}, x ∈ X . Then R = ∪ x∈X R x = {(α x , β x ) : x ∈ X }. Furthermore, let: ∆(x) := ∆ αx,βx (x) = , x ∈ X , with 0 < ≤ λ 1 .
Let P be the transition probability matrix of the original discrete time model. The transition matrix of the modified chain is then given by:

P = P + x∈X (α,β)∈Rx ∆ α,β (x)T α,β (x) = P + x∈X T αx,βx (x).
Relation ( 9) for (α x , β x ), x ≥ 0, is equivalent to convexity of functions v t , t ≥ 0:

v t (x + 1) ≤ 1 2 v t (x + 2) + 1 2 v t (x), x ≥ 0, t ≥ 0.
Thus, if we prove that v t , t ≥ 0 are convex, then Theorem 2 implies (15), as our reward function r is positive, and ( 10) holds from the definition of .

In the proof of convexity of v t , t ≥ 0, we will also use the fact that for each t ≥ 0, the function v t is non-decreasing, which was shown in Example 1, relation [START_REF] Bijvank | Service inventory management: Solution techniques for inventory systems without backorders[END_REF].

Assume that v t is convex for a given t ≥ 0 (for t = 0, v 0 := 0 is obviously convex). Then for t + 1 we have (see [START_REF] Ben Mamoun | Increasing convex monotone Markov chains: Theory, algorithm and applications[END_REF] in Example 1):

v t+1 (x) = r(x) + λ 1 A 1 (x, t) + λ 2 A 2 (x, t) + µS(x, t), x ≥ 0.
We consider separately one period rewards, arrivals of each type, and service.

• One period rewards. v 1 = r is obviously convex.

• Type 1 arrivals. For x ≥ 0, A 1 (x + 2, t) + A 1 (x, t) -2A 1 (x + 1, t) = v t (x + 3) + v t (x + 1) - 2v t (x + 2) ≥ 0, since v t is convex. • Type 2 arrivals. For x ≥ 0, A 2 (x + 2, t) + A 2 (x, t) -2A 2 (x + 1, t) = v t (x + 4) + v t (x + 2) - 2v t (x + 3) ≥ 0, since v t is convex. • Service. For x ≥ 0, S(x + 2, t) + S(x, t) -2S(x + 1, t) = v t (x + 1) + v t (x -1)1 {x>0} + v t (x)1 {x=0} -2v t (x) = 1 {x>0} (v t (x + 1) + v t (x -1) -2v t (x)) + 1 {x=0} (v t (x + 1) -v t (x)) ≥ 0, since v t is non-decreasing and convex.
Thus v t+1 is convex, so by induction on t, v t is convex for all t ≥ 0. Applying Theorem 2 we have v t (x) ≤ v t (x), x ≥ 0, t ≥ 0, that is the mean number of jobs in the system increases if we have more variable arrivals.

Remark 3. The goal of this example is only to illustrate the generalization of the precedence relation method. Note that the above result can be also obtained by using stochastic recurrences and icx-order (see for instance Baccelli and Brémaud [START_REF] Baccelli | Elements of queuing theory[END_REF]) or using event based dynamic programming (Koole [START_REF] Koole | Structural results for the control of queueing systems using event-based dynamic programming[END_REF], Koole [START_REF]Monotonicity in Markov reward and decision chains: theory and applications[END_REF]).

A primary use of precedence relations is to enable bounds to be established by analyzing simpler (smaller) systems. As mentioned in the introduction, one common simplification of a chain is to reduce its state space using aggregation. We examine this technique next.

Aggregation.

In this and the following sections we assume that the state space of the chain is finite, as we will use results in Kemeny and Snell [START_REF] Kemeny | Finite Markov chains[END_REF] on the aggregation of finite Markov chains. Let C = {C k } k∈K be a partition of the state space X into macro-states:

∪ k∈K C k = X , C i ∩ C j = ∅, ∀i = j.
Definition 1. (Kemeny and Snell [START_REF] Kemeny | Finite Markov chains[END_REF]) A Markov chain X = {X n } n≥0 is strongly aggregable (or lumpable) with respect to partition C if the process obtained by merging the states that belong to the same set into one state is still a Markov chain, for all initial distributions of X 0 .

There are necessary and sufficient conditions for a chain to be strongly aggregable: Theorem 3. (Matrix characterization, Kemeny and Snell [16, Theorem 6.3.2]) A DTMC X = {X n } n≥0 with probability transition matrix P is strongly aggregable with respect to C if and only if:

∀i ∈ K, ∀j ∈ K, y∈Cj P [x, y] is constant for all x ∈ C i . ( 16
)
Then we can define a new (aggregated) chain Y = {Y n } n≥0 with transition matrix Q. For all i, j ∈ K:

Q[i, j] = y∈Cj P [x, y], x ∈ C i .
There are many results on aggregation of Markov chains, however they primarily consider the stationary distribution. Surprisingly, we were not able to find the following simple property, so we provide it here with a proof. Proposition 2. Let X = {X n } n≥0 be a Markov chain satisfying [START_REF] Kemeny | Finite Markov chains[END_REF] and Y = {Y n } n≥0 the aggregated chain. Let r : X → R be a reward function that is constant within each macro-state, i.e. there exist r k ∈ R, k ∈ K such that for all k ∈ K:

r(x) = r k , ∀x ∈ C k .
Denote by v t and w t the t-period rewards for chains X and Y . Then for all k ∈ K:

v t (x) = w t (k), x ∈ C k , t ≥ 0. ( 17
)
Proof. We will show [START_REF] Kök | Assortment planning: Review of literature and industry practice, Retail Supply Chain Management[END_REF] by induction on t. Suppose that ( 17) is satisfied for t ≥ 0 (for t = 0 this is trivially satisfied). Then for t + 1 and

k ∈ K: v t+1 (x) = r(x) + y∈X P [x, y]v t (y) = r k + j∈K y∈Cj P [x, y]v t (y).
By the induction hypothesis v t (y) = w t (j), j ∈ K, y ∈ C j , and from Theorem 3, for

x ∈ C k , y∈Cj P [x, y] = Q[k, j], j ∈ K. Thus: v t+1 (x) = r k + j∈K Q[k, j]w t (j) = w t+1 (k).
By taking the limit, the above result also gives us the equality of mean stationary rewards.

Corollary 2. Let X and Y be two Markov chains satisfying the assumptions of Proposition 2.

If both chains are stable (which is equivalent to having only one recurrent class in the case of a finite chain), then the mean stationary rewards satisfy a = a (where a = lim t→∞ vt(x) t , for any x from the unique recurrent class of the chain X, and a = lim t→∞ wt(k) t

for any k from the unique recurrent class of the aggregated chain Y ).

5 Assemble-to-order systems.

Although assemble-to-order (ATO) systems typically refer to systems in which components are kept to stock and the end-product is only assembled when a customer demand occurs, the same model can also be used in much more general settings. The ATO model can be used to model processes at, for instance, supermarkets, Dell (Kapuscinski et al. [START_REF] Kapuscinski | Inventory decisions in dell's supply chain[END_REF]), Amazon (Xu et al. [46]), maintenance companies (Vliegen and van Houtum [45]) and for assortment planning (Kök et al. [START_REF] Kök | Assortment planning: Review of literature and industry practice, Retail Supply Chain Management[END_REF]).

In this section, we first present a general framework for ATO systems, in which these different applications can be included (Section 5.1). Then, in Section 5.2, we give a detailed model description for one instance of our framework, the so-called ATO models with partial order service (POS). For this specific model we prove a monotonicity property with respect to supermodularity (Section 5.3). Afterwards, in Section 5.4, we discuss the other instances within the framework, and whether or not the same property holds. In Section 6, we apply this monotonicity property to prove results for two specific examples from our framework.

Framework.

ATO systems have received a lot of attention in literature over the last decades. For an overview of research in this area, see Song and Zipkin [START_REF] Song | Supply chain operations: Assemble-to-order systems, Supply Chain Management: Design, Coordination and Operation[END_REF] and, more recently, Bijvank [START_REF] Bijvank | Service inventory management: Solution techniques for inventory systems without backorders[END_REF].

In these overview papers (Song and Zipkin [START_REF] Song | Supply chain operations: Assemble-to-order systems, Supply Chain Management: Design, Coordination and Operation[END_REF], and Bijvank [START_REF] Bijvank | Service inventory management: Solution techniques for inventory systems without backorders[END_REF]), several distinguishing characteristics are mentioned. A first distinction is between one-period models, periodic review models, and continuous review ATO models. We focus on continuous-review models. A second characteristic is the replenishment of the components, which is influenced by the service capacity and the lead time distribution of the replenishment times. We focus on exponential replenishment times. Third, the way an out-of-stock situation is handled is important. Three options are possible: A demand can be backordered; fulfilled partly (just the components that are available), referred to as partial order service (POS); or lost fully, referred to as total order service (TOS). We discuss all three of these options below. Finally, Song and Zipkin [START_REF] Song | Supply chain operations: Assemble-to-order systems, Supply Chain Management: Design, Coordination and Operation[END_REF] distinguish between studies on optimal policies and studies on performance evaluation for a given policy. In this work, we focus on the performance evaluation of ATO systems. As is common in the ATO literature, we assume at most one component of each type will be demanded for any item (the unit demand cases).

One characteristic not yet studied extensively in the field of ATO systems is the effect of coupling in returns/replenishments of components used in the assembly of the end-product. It is usually assumed that components have independent replenishment times. However, when service tools for maintenance actions are used, they are returned together (Vliegen and van Houtum [45]). This is also true for other products that are borrowed or rented, i.e., a rental car together with the child seat or the GPS system. But even for replenishment systems, a coupling between the replenishment moments of components might exist: a supplier might use the same shipment for sending all components. Therefore, we distinguish between two different situations for the replenishment/return of components: components are either returned individually (and independently) or jointly. This distinction can be seen on the vertical axis of our framework (see Table 1). On the horizontal axis of the framework, we distinguish the different policies for an out-of-stock situation, as also seen in Song and Zipkin [START_REF] Song | Supply chain operations: Assemble-to-order systems, Supply Chain Management: Design, Coordination and Operation[END_REF]: backordering, total order service, and partial order service. Thus six different returns/out of stock combinations can be distinguished. We now shortly discuss the related literature of all these pairs. A1. In the upper left corner (A1), ATO systems with individual returns and total order service are considered, for instance, Dell computers (Kapuscinski et al. [START_REF] Kapuscinski | Inventory decisions in dell's supply chain[END_REF]). A customer wants to have a computer with some specific elements. If one of these parts is not available, the customer does not want to have the computer at all: the demand for all parts is lost. After parts are used in the assembly of a computer, they are replenished with part specific replenishment rates. This kind of model has also been studied by Song et al. [START_REF] Song | Order-fulfillment performance measures in an assemble-toorder system with stochastic leadtimes[END_REF]. They present an exact evaluation of ATO systems with backordering (backlog queue is infinite) and lost sales (with both POS and TOS; backlog queue equal to zero) by using a matrix geometric approach. This exact method, however, needs a B1. In the lower left corner (B1), models are considered that have joint returns of resources used plus total order service. The best known example of this is telecommunication systems, or specifically loss networks (Kelly [START_REF] Kelly | Loss networks[END_REF]). In these networks, demands arrive, for example a phone call, that need several links to be simultaneously available. If all links are available, the call is completed. After the call is finished, all links simultaneously become available again. When one or more of the links is not available, the call is not completes, . and the demand for all links is lost (total order service). Although loss networks have a product-form solution (Ross [START_REF] Ross | Multiservice loss models for broadband telecommunication networks[END_REF]), exactly computing the blocking probabilities for this system is known to be a #P-complete problem (Louth et al. [START_REF] Louth | Computational complexity of loss networks[END_REF]) due to the normalizing constant.

A2. In the upper middle section (A2), ATO systems with individual returns and partial order service are considered. A example from practice is grocery shopping: Customer come to a supermarket with a list of things they want to buy. When one of the items is not available, the customer still takes the rest of the items home. He or she, however, will be less satisfied after leaving the supermarket. Furthermore, there is no direct relationship between the items that were demanded together and their replenishment time; especially when items have different suppliers, there is no coupling in returns. In the literature, this kind of model has been studied by, among others, Song et al. [START_REF] Song | Order-fulfillment performance measures in an assemble-toorder system with stochastic leadtimes[END_REF]. As already discussed, Song et al. [START_REF] Song | Order-fulfillment performance measures in an assemble-toorder system with stochastic leadtimes[END_REF] studied a model that could apply to both POS and TOS, and presented an exact, but time-consuming, evaluation for this model. Dayanik et al. [START_REF] Dayanik | The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-to-order systems[END_REF] studied bounds on the performance of the system. Furthermore, the ATO-POS model is another special case of the model discussed in Iravani et al. ([13]).

B2. In the lower middle section (B2), ATO systems with joint returns and partial order service are considered. An example of this element of the framework is the service tool problem as considered in Vliegen and van Houtum [45]. In this problem, to perform a maintenance action, several service tools are needed at the same time. Whenever one or more tools are not present, they are sent by an emergency shipment to be able to start the maintenance action as soon as possible. For the supply location under consideration the demand for these emergency shipped tools is lost. Furthermore, after usage all tools return to the location they were sent from together. So, there is also a coupling in returns. Vliegen and van Houtum [45] developed three different approximations, of which the third one (a weighted average of approximation 1 and 2) leads to accurate and, for small instances, efficient results. For larger instances, however, these approximations are still time consuming.

A3. In the upper right corner (A3), ATO systems with individual returns and backordering are considered. Within this portion of the framework (and also in B3), a distinction needs to be made with regard to the order in which backorders are filled. While most papers assume the First-Come-First-Served (FCFS) rule, Lu et al. [START_REF] Lu | No-holdback allocation rules for continuous-time assemble-toorder systems[END_REF] analyze the class of no-holdback (NHB) rules. Under NHB rules, a demand is backordered if and only if at least one of its components is out of stock. In this case, components are not put aside as committed stock as is done in FCFS. In our framework, we consider the FCFS rule only. A real-life example for section A3 is a mail-order or e-commerce retailer. Often customer orders consisting of multiple items are shipped only when all items are present; the whole order therefore might be delayed/backordered. This model is among others studied by Song et al. [START_REF] Song | Order-fulfillment performance measures in an assemble-toorder system with stochastic leadtimes[END_REF] in the case where the backlog queue is infinitely long. Song and Yao [START_REF] Song | Performance analysis and optimisation of assemle-to-order systems with random lead times[END_REF] study a single-item model for which they derive an exact evaluation as well as easy-to-computer bounds. Lu et al. [START_REF] Lu | Order fill rate, leadtime variability, and advance demand information in an assemble-to-order system[END_REF] extend this model to a multi-item setting and derive the joint queue-length distribution, and Lu et al. [START_REF]Backorder minimization in multiproduct assemble-to-order systems[END_REF] study the optimal component inventories under a given budget so at to minimize the weighted average backorders. Furthermore, they derive bounds and approximations for the expected number of backorders, formulate surrogate optimization problems and develop algorithms to solve these problems. Lu [START_REF] Lu | Performance analysis for assemble-to-order systems with general renewal arrivals and random batch demands[END_REF] generalizes the model of Lu et al. [START_REF]Backorder minimization in multiproduct assemble-to-order systems[END_REF] to a situation with general assumptions of demand patterns and replenishment leadtime distributions.

B3. In the lower right corner (B3), ATO systems with joint returns and backordering are considered. This situation can be seen in online retail settings, where items are supplied from the same supplier. In this case, the customer order, which is either fulfilled or backlogged, triggers a replenishment order for the items demanded. If these items are delivered by the same supplier, a coupling in the return epoch might exist. We are not aware of any papers in this part of the framework that assume exponential lead times. So, the model of this part of our framework has not been applied so far. However, for a slightly different situation, where equal deterministic replenishment times are assumed (in special cases of Song [START_REF] Song | On the order fill rate in a multi-item, base-stock inventory system[END_REF], Song [START_REF]Order-based backorders and their implications in multi-iyem inventory systems[END_REF]), joint returns in an ATO system with backordering do exist.

We can conclude that for all sections of our framework bounding models and/or approximations have been developed, or might be useful. For loss networks (B1), however, a product-form solution exists, which makes the exact calculation somewhat easier. We therefore focus on the other sections of the framework.

In the next subsection, we will give a detailed model description for the ATO model with POS and individual returns (A2). For this specific model we prove a monotonicity property with respect to supermodularity in Section 5.3. Afterwards, in Section 5.4, we discuss the other categories of the framework, and whether or not the same property holds.

Model description ATO-POS (A2).

We consider an assemble-to-order system with POS and individual returns, that we will call the ATO-POS model. We discuss some possible extensions to other variants (backordering and TOS) in Section 5.4.

The set of items is denoted by I and the number of item types I = |I|. We assume a base stock replenishment policy and denote by S = {S 1 , . . . , S I } the base stock vector. The state of the system is given by vector x = (x 1 , . . . , x I ), where x i is the number of items of type i in replenishment. We assume exponential demands and replenishments. There is a server for each item type and the replenishment rate of item i may be state dependent: the service rate for item type i in state x ∈ X is equal to µ i (x). Furthermore, we assume that replenishment rates for type i items depend only on x i , i.e.: for all x, y ∈ X such that x i = y i , µ i (x) = µ i (y), and we will write µ i (x i ) for µ i (x). We will assume further that

i∈I µ i (x i ) > 0, ∀x = (0, . . . , 0). ( 18
)
We can have demands for sets of items: λ A denotes the demand rate of subset ∅ = A ⊂ I. Then λ = ∅ =A⊂I λ A denotes the total demand rate. We assume partial order fulfillment: if some of demanded items are not on stock, then the present items are sent from stock and the missing items are considered as lost sales. A system with 2 items is given in Figure 5.

µ 2 (x 2 ) µ 1 (x 1 ) S 1 -x 1 λ 2 λ 1,2 λ 1 S 2 -x 2
X = {x = (x 1 , . . . , x I ) | 0 ≤ x i ≤ S i , i ∈ I} ,
where x i denotes the number of items of type i in replenishment. In the following, for i ∈ I, we denote by e i the state with all the components equal to 0 except the component i that is equal to 1.

There are two types of transitions:

• Demands. For each subset ∅ = A ⊂ I, a demand of set A occurs with rate λ A . The new state is described by operator d A : X → X :

d A (x) = x + i∈A 1 {xi<Si} e i .
• Replenishment. For each item type i ∈ I, the new state after replenishment of item i is described by operator r i : X → X :

r i (x) = x -1 {xi>0} e i .
Denote by γ i the maximal service rate for item type i:

µ i (x i ) ≤ γ i , x ∈ X .
Then the outgoing rate for each state is smaller than ∅ =A⊂I λ A + i∈I γ i . We will take the somewhat unconventional uniformization constant:

ν = ∅ =A⊂I |A|λ A + i∈I γ i .
The choice of this uniformization constant will be very useful in Section 6, where we will establish a monotonicity property for the ATO-POS model (A2) with respect to the amount of coupling in the system. We consider a Discrete Time Markov Chain (DTMC) after uniformization with parameter ν. To simplify notation and without loss of generality, we assume that ν = 1. Then a demand of item set A in the uniformized chain occurs with probability λ A and the replenishment probability of item i is state-dependent and equal to µ i (x i ) in state x ∈ X . Figure 6 illustrates the state space and transitions for four different states in the 2-item case (the loops are omitted). Relation 18 implies that the Markov chain has only one recurrent class (and it is finite so this implies stability), so that the mean stationary rewards will be always well defined for our ATO-POS model.

x 1 S 1 S 2 x 2 µ 1 (2) µ 1 (1) λ 2 λ 2 λ 1 λ 1 µ 2 (4) λ 1 + λ 1,2 µ 1 (5) λ 2 + λ 1,2 λ 1,2 λ 1,2 µ 1 (3) µ 2 (2) µ 2 (2)
Figure 6: State space for the 2-item case.

ATO-POS model (A2) and supermodular rewards.

Let r : X → R be any reward function and denote by v t : X → R, t ≥ 0, its cumulative t-period reward:

v t+1 (x) = r(x)+ ∅ =A⊂I λ A v t (d A (x))+ i∈I µ i (x i )v t (r i (x))+   ν - ∅ =A⊂I λ A - i∈I µ i (x i )   v t (x), t ≥ 0, (19)
with v 0 (x) := 0, x ∈ X . To simplify the discussion, denote by D A (x, t) : X × N 0 → R and R i (x, t) : X × N 0 → R respectively the operators for demand and replenishment terms:

D A (x, t) = λ A v t (d A (x)), ∅ = A ⊂ I R i (x, t) = µ i (x i )v t (r i (x)), i ∈ I.
Furthermore, we define the uniformization operator U : X × N 0 → R as:

U (x, t) =   ν - ∅ =A⊂I λ A - i∈I µ i (x i )   v t (x) =   ∅ =A⊂I (|A| -1)λ A + i∈I (γ i -µ i (x i ))   v t (x). Then: v t+1 (x) = r(x) + ∅ =A⊂I D A (x, t) + i∈I R i (x, t) + U (x, t), x ∈ X , t ≥ 0, (20)
with v 0 (x) := 0, x ∈ X . We now prove the monotonicity property of the ATO-POS model with respect to supermodularity: Proposition 3. If the reward function r is supermodular, then for all t ≥ 0, t-period reward v t is also a supermodular function.

Proof. After Proposition 7 in Appendix A, proving supermodularity of v t is equivalent to showing that for all i, j ∈ I, i = j, and for all x ∈ X such that x + e i + e j ∈ X :

v t (x + e i + e j ) + v t (x) ≥ v t (x + e i ) + v t (x + e j ), ∀t ≥ 0. ( 21
)
We will show relation ( 21) by induction on t. Suppose that relation [START_REF] Lindvall | Lectures on the coupling method[END_REF] holds for a given t ≥ 0 (the case t = 0 is trivial). We will show that then it also holds for t + 1. Let i, j ∈ I, i = j, be arbitrary and fixed. We need to show that for all x ∈ X such that x + e i + e j ∈ X :

v t+1 (x + e i + e j ) + v t+1 (x) ≥ v t+1 (x + e i ) + v t+1 (x + e j ). (22)
We will prove this by considering demands separate from replenishments and uniformization terms. For demands we will prove the following relations for all x ∈ X such that x + e i + e j ∈ X :

D A (x + e i + e j , t) + D A (x, t) ≥ D A (x + e i , t) + D A (x + e j , t), ∅ = A ⊂ I. (23)
Due to the state-dependent replenishment rates, we need to consider replenishments and uniformization terms together. The uniformization term U (x, t) can be decomposed into the following I + 1 parts:

U (x, t) = U D (x, t) + k∈I U k (x, t), (24)
where:

U D (x, t) =   ∅ =A⊂I (|A| -1)λ A   v t (x), U k (x, t) = (γ k -µ k (x k ))v t (x), k ∈ I.
We will prove the following relations involving replenishments and uniformization terms for all x ∈ X such that x + e i + e j ∈ X :

R k (x + e i + e j , t) + U k (x + e i + e j , t) + R k (x, t) + U k (x, t) ≥ R k (x + e i , t) + U k (x + e i , t) + R k (x + e j , t) + U k (x + e j , t), k ∈ I, (25)
and

U D (x + e i + e j , t) + U D (x, t) ≥ U D (x + e i , t) + U D (x + e j , t). (26)
Then [START_REF] Louth | Computational complexity of loss networks[END_REF] follows from ( 20), ( 23), ( 24), ( 25), [START_REF]Backorder minimization in multiproduct assemble-to-order systems[END_REF], and supermodularity of the reward function r.

Demands. Consider a demand of an arbitrary and fixed subset A ⊂ I, A = ∅, and assume that λ A > 0 (the claim is trivial otherwise). Denote by x = x + a∈A 1 {xa<Sa} e a . We have 3 different cases:

• Case 1: i ∈ A and j ∈ A. Then ( 23) is equivalent to: v t (x + e i + e j ) + v t (x ) ≥ v t (x + e i ) + v t (x + e j ), which follows by induction hypothesis for state x .

• Case 2: i ∈ A and j ∈ A (case i ∈ A and j ∈ A is symmetric). Then we have:

v t (x + 1 {xi<Si-1} e i + e j ) + v t (x ) ≥ v t (x + 1 {xi<Si-1} e i ) + v t (x + e j ).
If x i = S i -1, then this trivially holds. If x i < S i -1, then the above relation follows by the induction hypothesis for state x .

• Case 3: i ∈ A and j ∈ A. Then we have:

v t (x + 1 {xi<Si-1} e i + 1 {xj <Sj -1} e j ) + v t (x ) ≥ v t (x + 1 {xi<Si-1} e i ) + v t (x + 1 {xj <Sj -1} e j ).
If at least one of x i = S i -1 or x j = S j -1 holds, then this is trivially satisfied. If x i < S i -1 and x j < S j -1, then the above relation follows by the induction hypothesis for state x .

Replenishments. Consider item k ∈ I. We have 2 different cases:

• Case 1: k ∈ {i, j}. Then ( 25) is equivalent to:

µ k (x k )v t (x + e i + e j -1 {x k >0} e k ) + (γ k -µ k (x k ))v t (x + e i + e j ) + µ k (x k )v t (x -1 {x k >0} e k ) + (γ k -µ k (x k ))v t (x) ≥ µ k (x k )v t (x + e i -1 {x k >0} e k ) + (γ k -µ k (x k ))v t (x + e i ) + µ k (x k )v t (x + e j -1 {x k >0} e k ) + (γ k -µ k (x k ))v t (x + e j ). (27)
By the induction hypothesis for state x -1 {x k >0} e k :

v t (x + e i + e j -1 {x k >0} e k ) + v t (x -1 {x k >0} e k ) ≥ v t (x + e i -1 {x k >0} e k ) + v t (x + e j -1 {x k >0} e k ),
and by the induction hypothesis for state x: v t (x + e i + e j ) + v t (x) ≥ v t (x + e i ) + v t (x + e j ), so ( 27) holds.

• Case 2: k = i (k = j is symmetric). Then [START_REF] Lu | Order fill rate, leadtime variability, and advance demand information in an assemble-to-order system[END_REF] becomes:

µ i (x i + 1)v t (x + e j ) + (γ i -µ i (x i + 1))v t (x + e i + e j ) + µ i (x i )v t (x -1 {xi>0} e i ) + (γ i -µ i (x i ))v t (x) ≥ µ i (x i + 1)v t (x) + (γ i -µ i (x i + 1))v t (x + e i ) + µ i (x i )v t (x + e j -1 {xi>0} e i ) + (γ i -µ i (x i ))v t (x + e j ).
This is equivalent to:

(γ i -µ i (x i + 1)) (v t (x + e i + e j ) + v t (x) -v t (x + e i ) -v t (x + e j ) I ) + µ i (x i )(v t (x + e j ) + v t (x -1 {xi>0} e i ) -v t (x) -v t (x + e j -1 {xi>0} e i ) II ) ≥ 0.
For part I, by the induction hypothesis for state x: v t (x + e i + e j ) + v t (x) -v t (x + e i ) -v t (x + e j ) ≥ 0.

For part II:

-If x i > 0, then by the induction hypothesis for state x -e i :

v t (x + e j ) + v t (x -e i ) -v t (x) -v t (x + e j -e i ) ≥ 0.
-If x i = 0, then:

v t (x + e j ) + v t (x) -v t (x) -v t (x + e j ) = 0.
Uniformization terms. Assume ∅ =A⊂I (|A| -1)λ A > 0 (the claim is trivial otherwise). Then ( 26) is equivalent to:

v t (x + e i + e j ) + v t (x) ≥ v t (x + e i ) + v t (x + e j ), (28)
which is satisfied by the induction hypothesis for state x.

It is easy to check that the proof remains valid if we invert the inequalities. Thus, the ATO-POS model is also monotone with respect to submodularity: Proposition 4. If the reward function r is submodular, then for all t ≥ 0, t-period reward v t is also a submodular function.

Examples. We give now some examples of supermodular reward functions.

Order fill rates for individual demand streams. Denote by OFR A : X → {0, 1} the order fill rate for the demand stream for set A ⊂ I, A = ∅:

OFR A (x) = i∈A 1 {xi<Si} , x ∈ X . ( 29
)
Lemma 1. Reward function OFR A is supermodular for all A ⊂ I, A = ∅.

Proof. Let A ⊂ I, A = ∅, and i, j ∈ I, i = j. We need to show that for all x ∈ X such that

x + e i + e j ∈ X (see Proposition 8):

OFR A (x + e i + e j ) + OFR A (x) ≥ OFR A (x + e i ) + OFR A (x + e j ). ( 30
)
Suppose first that i, j ∈ A. We have 4 different cases:

• There is a k ∈ A\{i, j} such that x k = S k . Then OFR A (x + e i + e j ) = OFR A (x) = OFR A (x + e i ) = OFR A (x + e j ) = 0, and relation [START_REF] Muntz | Bounding availability of repairable computer systems[END_REF] clearly holds.

• For all k ∈ A\{i, j}, x k < S k , x i = S i -1 and x j = S j -1. Then OFR A (x) = 1 and all the other terms are equal to 0, thus relation ( 30) holds.

• For all k ∈ A\{i, j}, x k < S k , x i = S i -1 and x j < S j -1 (x i < S i -1 and x j = S j -1 is symmetrical). Then OFR A (x) = OFR A (x + e j ) = 1 and the other two terms are equal to 0, so the both sides of relation ( 30) are equal to 1.

• Finally, if x k < S k , ∀k ∈ A\{i, j}, x i < S i -1, and x j < S j -1, then the both sides of ( 30) are equal to 2.

If i ∈ A, then OFR A (x + e i + e j ) = OFR A (x + e j ) and OFR A (x) = OFR A (x + e i ) so [START_REF] Muntz | Bounding availability of repairable computer systems[END_REF] clearly holds.

Aggregated order fill rate. Denote by OFR : X → {0, 1} the aggregate order fill rate:

OFR(x) = 1 λ ∅ =A⊂I λ A OFR A (x), x ∈ X . ( 31
)
Lemma 2. The aggregated order fill rate OFR is supermodular.

Proof. A linear combination of supermodular functions is supermodular.

Since the proofs for other reward functions are straightforward analogs to the proofs of Lemma 1 and 2, further proofs are omitted. Other supermodular reward functions, interesting in the framework of ATO systems, are:

• Item fill rates for individual items. Denote by IFR i : X → {0, 1} the item fill rate for item i ∈ I:

IFR i (x) = 1 xi<Si , x ∈ X (32)
• Aggregated item fill rate. Denote by IFR : X → {0, 1} the aggregated item fill rate:

IFR(x) = i∈I f i 1 xi<Si , x ∈ X, (33)
where f i ≥ 0 is the weight factor for item i, i f i = 1. This would be an interesting reward function when the availability of the items is not equally important, and the experience of the customers is not solely determined by the availability of all items.

• Average amount of backorders/queue length for individual items. Denote by BO i : X → R + the average amount of backorders for item i ∈ I:

BO i (x) = (x i -S i ) + , x ∈ X ( 34
)
This would be an interesting reward function for the right portion of the framework, ATO systems with backordering (A3 and B3).

• Aggregate amount of backorders/queue length. Denote by BO : X → R + the aggregate average amount of backorders:

BO(x) = i∈I f i (x i -S i ) + , x ∈ X, (35)
where f i ≥ 0 is the weight factor for item i, i f i = 1. Also this reward function is interesting for the right portion of the framework, ATO systems with backordering (A3 and B3).

Further discussion and extensions.

In Section 5.3, we proved a monotonicity property with respect to supermodularity (Proposition 3) or submodularity (Proposition 4) for the ATO-POS model (A2). We consider here the possible extensions of these results to the individual return model with backordering (A3) and TOS (A1).

For the lower part (B) of the framework, supermodularity is not well-defined. The monotonicity property shown for the ATO-POS model therefore does not make sense for the models in this part. However, the precedence relation method combined with aggregation can be used to prove bounds for the performance of the models, as is shown in Section 6.2.

Backordering (A3)

. Assume now that demands that cannot be immediately fulfilled are backordered. More precisely, the items that are available are sent immediately to the customer (or kept apart as reserved) and the missing items are sent as soon as they become available. The state-space of this model, referred to as ATO-B in the following, is:

X = {x = (x 1 , . . . , x I ) | x i ≥ 0, i ∈ I} ,
where x i denotes the number of items of type i in replenishment, as before.

Returns remain the same as in the ATO-POS model, while demands become simpler as we no longer have the boundary effect at x i = S i for i ∈ I. The operator d A becomes:

d A (x) = x + i∈A e i .
In this case, similar results as in propositions 3 and 4 hold for the ATO-B model: Proposition 5. Let X be a Markov chain of the ATO-B model. If the reward function r is supermodular (resp. submodular), then for all t ≥ 0, t-period reward v t is also a supermodular (resp. submodular) function.

Proof. The proof is similar to the proof of Propositions 3. The models differ only in demands, so we only need to prove the analog to relation [START_REF] Lu | Performance analysis for assemble-to-order systems with general renewal arrivals and random batch demands[END_REF], for any subset A ⊂ I, A = ∅, such that λ A > 0. As there is no boundary effect for demands, we no longer need to consider three different cases. Denote x = x + a∈A e a . Then for the ATO-B model relation ( 23) is equivalent to:

v t (x + e i + e j ) + v t (x ) ≥ v t (x + e i ) + v t (x + e j ),
which follows by the induction hypothesis for state x . TOS (A1). Unfortunately, supermodularity does not hold for ATO systems with TOS. The primary reason for this is that whenever one item is out-of-stock this influences the demand for other items.

To illustrate this, we now present a small example of this fact; see Table 5.4. In this example, we compare two situations, both including three items with S = {2, 2, 2}. Situation 1 has more coupling in demands than situation 2, while the aggregate demand for each item remains the same. Specifically in situation 2, the demand for set {1, 3} is partly split into separate demands for item {1} and item {3}. This decoupling leads to decreased service for sets {1, 3} and {1, 2, 3} (see Table 5.4), since the stock levels for items 1 and 3 are typically lower in situation 2. The reason for this is that the unavailability of one item does not lead to a lost demand for the other item as frequently. However, since the service level for set {1, 2, 3} is also lower in situation 2, the availability of item 2 will be higher. Therefore, the fill rate for set {2} will be higher (see Table 5.4). While in the other models (ATO-POS and ATO-B) the decoupling of demand leads to a situation where either the service levels decrease or stay equal, decoupling in the ATO-TOS model can also lead to an increase in service for some demand streams. 6 Applications.

Input Situation 1 Situation 2 µ 1 1 λ 1 0 0.3 λ 2 1 1 λ 3 0 0.3 λ 1,2 0 0 λ 1,3 0.5 0.2 λ 2,3 0 0 λ 1,
In this section we illustrate our analytical technique -applying precedence relations on sets of states -on two applications. In Section 6.1, we consider ATO-POS model and we prove a monotonicity property with respect to the amount of coupled demands, valid for any supermodular reward function. This result has been previously established in Dayanik et al. [START_REF] Dayanik | The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-to-order systems[END_REF] for capacitated (i.e. single server) ATO-POS model and for order fill rates. We extend here their result to the more general case of state dependent replenishment times (this covers many standard assumptions, such as capacitated server, ample server, or many server case) and for general supermodular rewards. In Section 6.2, we consider the service tool problem, introduced by Vliegen and Van Houtum [45]. We apply our generalized precedence relation method together with aggregation to establish a formal proof of bounds conjectured in [45] based on numerical experiments.

ATO-POS and coupling in demands.

Let X be a Markov chain of the ATO-POS model as described in Section 5.2. Consider a new ATO-POS model Y that is obtained from X by modifying its demand rates in order to decrease the coupling in demands. More formally, for each A ⊂ I such that |A| > 1 let A be a constant 0 ≤ A ≤ λ A . We can define the chain Y by modifying the demand rates to:

λ A = λ A -A , λ i = λ i + {A⊂I | |A|>1, i∈A} A .
We have the following monotonicity property of ATO-POS model with respect to the amount of coupling in the system: Proposition 6. For any supermodular reward function, the chain Y gives lower bounds for mean cumulated and stationary rewards of the chain X.

Proof. Chain Y is obtained from X by replacing the demands of sets of items by individual demands. To describe this transformation formally, recall that for each state x ∈ X and each A ⊂ I, A = ∅, d A (x) = x + i∈A 1 {xi<Si} e i denotes the state after a demand of set A. For each x ∈ X and each A ⊂ I such that |A| > 1, we will define probability vectors α x,A and β x,A as follows:

α x,A (z) = 1 |A| δ d A (x) (z) + (|A| -1)δ x (z) , β x,A (z) = 1 |A| i∈A δ d {i} (x) (z), z ∈ X ,
and the weight factor ∆ x,A = A . Let P be the transition matrix of the original chain X. The transition matrix P of the chain Y is then given by:

P = P + x∈X A⊂I, |A|>1 |A|T α x,A ,β x,A (x).
Let r : X → R be any supermodular reward function and denote by v t its cumulative t-period reward, defined by relation [START_REF]Monotonicity in Markov reward and decision chains: theory and applications[END_REF]. Similarly, denote by w t the t-period rewards for the chain Y . If we show that for all x ∈ X and for all A ⊂ I such that |A| > 1, functions v t satisfy:

z∈X α x,A (z)v t (z) ≥ z∈X β x,A (z)v t (z), t ≥ 0, (36)
then by Theorem 2 it follows that:

v t (x) ≥ w t (x), x ∈ X , t ≥ 0. (37)
For x ∈ X and A ⊂ I such that |A| > 1, relation ( 36) is equivalent to:

v t (x + i∈A 1 {xi<Si} e i ) + (|A| -1)v t (x) ≥ i∈A v t (x + 1 {xi<Si} e i ), t ≥ 0. ( 38
)
By Proposition 3, for any t ≥ 0, v t is a supermodular function. From Proposition 7 3 in Appendix A, it follows that for all ∅ = A ⊂ I and x ∈ X such that x + i∈A e i ∈ X :

v t (x + i∈A e i ) + (|A| -1)v t (x) ≥ i∈A v t (x + e i ), t ≥ 0. ( 39
)
This now trivially implies relation [START_REF] Song | Performance analysis and optimisation of assemle-to-order systems with random lead times[END_REF], so we can apply Theorem 2.

x 1 S 1 S 2 x 2 µ 1 (2) µ 1 (1) λ 2 + λ 1,2 λ 2 + λ 1,2 λ 1 + λ 1,2 µ 2 (4) λ 1 + λ 1,2 µ 1 (5) λ 2 + λ 1,2 µ 1 (3) µ 2 (2)
µ 2 (2) In the extreme case, if we split all the joint demands into individual demands (i.e. if we take

λ 1 + λ 1,2
A = λ A , for all A ⊂ I such that |A| > 1)
, we obtain the independent stock assumption (ISA) model. In Figure 7 we give the state space and transitions for four different states in the 2-item case ISA model. We can also split a demand of set A into demands for sets A 1 , . . . , A K such that A = ∪ K k=1 A k and A i ∩ A j = ∅, i = j. Again using supermodularity characterization (53) (Proposition 7 in Appendix A) and the general precedence relation method, it can be easily proved that this is also a lower bound for the original model, which is tighter than the ISA model. This lower bound is also easier to solve than the original model as it allows us to decompose the problem into smaller subproblems.

More precisely, let I 1 , . . . , I K be any decomposition of the item set I (i.e. I = ∪ K k=1 I k and I i ∩I j = ∅, i = j). Then for any subset ∅ = A ⊂ I, we can define the sets

A k := A∩I k , 1 ≤ k ≤ K.
In the bounding model, a demand of set A is decomposed into demands of sets A 1 , . . . , A K . (If some of the sets A k , 1 ≤ k ≤ K are empty sets, we can ignore them.) This bounding model can then be decomposed into K subproblems where subproblem k involves only items in set I k . These bounds are generalization to the state dependent replenishment case of setwise bounds, derived in Dayanik et al. [START_REF] Dayanik | The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-to-order systems[END_REF] for the capacitated replenishment case.

Remark 4. The results remain valid for the backordering case (A3) of our framework, by using Proposition 5 for the ATO-B model (instead of Proposition 3) in the proof of Proposition 6. We only modify a finite number of transitions per state, so the result holds even for unbounded supermodular functions. However, in the case of the ATO-B model, stability is no longer guaranteed by relation [START_REF] Koole | Structural results for the control of queueing systems using event-based dynamic programming[END_REF], so special care should be taken to check whether the mean stationary rewards are well defined for given model parameters, which will typically be implied by

µ i (x i ) > i∈A λ A , x i > k i , ∀i ∈ I,
for some k i ≥ 0, i ∈ I.

Service tools.

The service tool model was introduced by Vliegen and Van Houtum [45]. In our framework in Table 1, this model corresponds to the B2 case.

As the items that are delivered together return together, we need to keep track of the sets of items delivered together. Consider a very simple case of two different item types. Then the state of the system is given by a vector (n {1} , n {2} , n {1,2} ) where n {1} (resp. n {2} ) is the number of items of type 1 (resp. 2) at the customer that were delivered individually, and n {1,2} is the number of sets {1, 2} at the customers that were delivered together and that will return together. Given, for example stock levels S 1 = 1, S 2 = 2, all possible states of the system are: (0, 0, 0), (1, 0, 0), (0, 1, 0), (1, 1, 0), (0, 0, 1), (0, 2, 0), (1, 2, 0) and (0, 1, 1). Note that if set {1, 2} is demanded, and item type 2 is out of stock, this becomes a demand for item type 1 (and similarly if item 1 is out of stock). The Markov chain for this case is given in Figure 8. Markov chain description. For the general I-item case, the state of the system is given by a vector n = (n A ) ∅ =A⊂I where n A ≥ 0 is the number of sets A at the customer that were delivered together. For each i ∈ I, we denote by ξ i (n) the total number of items of type i at the customer:

ξ i (n) = A⊂I 1 {i∈A} n A . ( 40
)
The state space of the model is then:

E = n = (n A ) ∅ =A⊂I : ξ i (n) ≤ S i , ∀i ∈ I .
For each A ⊂ I, A = ∅, we will denote by e A the state in which all the components are equal to 0, except the component A that is equal to 1.

We have the following transitions:

• Demands. For each subset ∅ = A ⊂ I, a demand of set A occurs with rate λ A . For all i ∈ A the amount of items of type i at the customer stays the same. In state n ∈ E, we can deliver an item of type i only if ξ i (n) < S i . The new state is described by operator φ A : E → E:

φ A (n) = n + e {i∈A : ξi(n)<Si} .
• Returns. In state n ∈ E, for each subset ∅ = A ⊂ I such that n A > 0, a return of set A occurs with rate µn A . The new state is described by operator ψ A : E → E:

ψ A (n) = n -1 {n A >0} e A .
We will consider the uniformized chain with uniformization constant

ν = ∅ =A⊂I |A|λ A + I i=1 S i µ.
Without loss of generality, we assume that ν = 1.

We will refer to this (uniformized) model as M 0 . Note that M 0 has only one recurrent class, as µ > 0. It is also aperiodic (note the loop at state (0, . . . , 0)). As the state space is finite, the stationary distribution always exists. Though its state space is finite, its dimension is equal to 2 I -1, thus the Markov chain becomes numerically intractable even for small values of I and S i , i ∈ I. For example, for I = 5 and S i = 5, ∀i, the cardinality of the state space is |E| = 210 832 854.

Bounds. The complexity of the original service tool model (Vliegen and Van Houtum [45]) comes from the need to track which items were delivered together. If we assume that all the items return individually, then there is no need to track which items were delivered together. We thus obtain the ATO-POS model discussed in Section 5.2, with the state space:

X = {x : 0 ≤ x i ≤ S i , i ∈ I},
where the state x = (x 1 , . . . , x I ) of the system represents the number of items at the customer for each item type, i.e. x i = ξ i (n), ∀i. Due to the fixed return rate for all item sets, we have µ i (x) = µx i , ∀x ∈ X , ∀i ∈ I. We will denote this ATO-POS model by M 1 .

The cardinality of X , |X | = I i=1 (S i + 1), is considerably lower than for the original model. For example, for I = 5 and S i = 5, ∀i, we have |X | = 7776 (compared to |E| = 210 832 854).

Note that model M 1 can be obtained from the original model M 0 in two steps:

1. By redirecting transitions that correspond to returns. We will denote by M 1 the model obtained from model M 0 by replacing all joint returns by individual returns (see Figure 9). For example, in model M 0 in state (0, 1, 1) we have one joint return of set {1, 2} and a uniformization loop; these are replaced by two new transitions: one to state (0, 2, 0) (corresponding to an individual return of item 1) and one to state (1, 1, 0) (an individual return of item 2).

2. Notice that the obtained model M 1 is lumpable with respect to the partition of the state space induced by function ξ = (ξ i ) i∈I , see [START_REF] Tandra | Join minimum cost queue for multiclass customers: Stability and performance bounds[END_REF], and model M 1 is the lumped version of M 1 . We now need not track the history of joint demands, only the total number of items of each type at the customer. Proof of the bounds. In the following, we will show that model M 1 gives a lower bound for the aggregate order fill rate of the original model. Model M 1 is obtained from the original model by replacing the returns of sets of items by individual returns. In order to describe this transformation formally, for each n ∈ E and each A ⊂ I such that |A| > 1 and n A > 0 we will define probability vectors α n,A and β n,A as follows:

α n,A (m) = 1 |A| (δ n-e A (m) + (|A| -1)δ n (m)) , β n,A (m) = 1 |A| i∈A δ n-e A +e A\{i} (m), m ∈ E,
and the weight factor ∆ n,A : ∆ n,A = µn A |A|.

Let P be the transition matrix of the original chain. The transition matrix P 1 of the chain M 1 is then given by:

P 1 = P + n∈E A⊂I, |A|>1 : n A >0 µn A |A|T α n,A ,β n,A (n).
For example, if I = 2 and S 1 = 1, S 2 = 2 (as in Figure 9, on the left), then for state n = (0, 0, 1) and A = {1, 2}: α (0,0,1),{1,2} = 1 2 δ (0,0,1)-e {1,2} + δ (0,0,1) = 1 2 δ (0,0,0) + δ (0,0,1) , β (0,0,1),{1,2} = 1 2 δ (0,0,1)-e {1,2} +e {2} + δ (0,0,1)-e {1,2} +e {1} = 1 2 δ (0,1,0) + δ (1,0,0) , and ∆ (0,0,1),{1,2} = 2µ. Vectors α (0,0,1),{1,2} and β (0,0,1),{1,2} formally describe the redirection of outgoing transitions from state (0, 0, 1) (see Figure 9, on the left): Old (blue/dashed) transitions to states (0, 0, 0) (joint return of set {1, 2}) and (0, 0, 1) (uniformization loop) are replaced by new (red/bold) transitions to states (0, 1, 0) (individual return of item type 1) and (1, 0, 0) (individual return of item type 2). The corresponding weight factor ∆ (0,0,1),{1,2} = 2µ states that the redirected amount of the transitions should be equal to α (0,0,1),{1,2} (0, 0, 1)∆ (0,0,1),{1,2} = P (0,0,1),(0,0,0) . Intuitively, we redirect entirely the transitions corresponding to joint returns. Note that after applying the above modification we still have a loop at state (0, 0, 1), but with a modified probability: the new probability of the loop transition is now equal to λ {1,2} + µ, compared to λ {1,2} + 2µ in the original chain (for this example, ν = λ {1} + λ {2} + 2λ {1,2} + 3µ). Let r : E → R be any reward function and let w t , t ≥ 0, denote the t-period rewards for model M 0 :

w t+1 (n) = r(n) + ∅ =A⊂I λ A w t (φ A (n)) + B⊂I : n B >0 µn B w t (ψ B (n)) +   ν - ∅ =A⊂I λ A -µ B⊂I : n B >0 µn B   w t (n), ∀n ∈ E,
where w 0 (n) := 0, n ∈ E. Similarly, denote by w t the t-period rewards for model M 1 . If we show that for all n ∈ E and for all A ⊂ I such that n A > 0, functions w t satisfy:

k∈E α n,A (k)w t (k) ≥ k∈E β n,A (k)w t (k), t ≥ 0, (41)
then by Theorem 2 it follows that:

w t (n) ≥ w t (n), n ∈ E, t ≥ 0. ( 42
)
For n ∈ E and A ⊂ I such that n A > 0, relation ( 41) is equivalent to:

w t (n -e A ) + (|A| -1)w t (n) ≥ i∈A w t (n -e A + e A\{i} ), t ≥ 0. ( 43
)
Due to the complex structure of the state space E, relation ( 43) is difficult to check (and might even not hold). However, ( 43) is only a sufficient condition for [START_REF] Van Dijk | Simple bounds and monotonicity results for finite multiserver exponential tandem queues[END_REF]. The "dual" sufficient condition for [START_REF] Van Dijk | Simple bounds and monotonicity results for finite multiserver exponential tandem queues[END_REF] is to show that for all n ∈ E and for all A ⊂ I, A = ∅, functions w t satisfy:

w t (n -e A ) + (|A| -1)w t (n) ≥ i∈A w t (n -e A + e A\{i} ), t ≥ 0. ( 44
)
Intuitively, instead of starting with model M 0 as the original model, we can start with model M 1 .

Then the transformation of model M 1 to model M 0 can be described using probability vectors α n,A = β n,A and β n,A = α n,A , and the weight factor ∆ n,A = ∆ n,A , for each n ∈ E and each A ⊂ I such that |A| > 1 and n A > 0. Transition matrices P (model M 0 ) and P 1 (model M 1 ) clearly satisfy:

P = P 1 + n∈E A⊂I : n A >0 µn A |A|T α n,A ,β n,A (n).
Furthermore, relation ( 44) is clearly equivalent to:

k∈E α n,A (k)w t (k) ≤ k∈E β n,A (k)w t (k), t ≥ 0. ( 45
)
So proving [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF], and using Theorem 2 and Corollary 1, we will show that model M 0 gives an upper bound for model M 1 , which is equivalent to showing that M 1 gives a lower bound for model M 0 . The advantage of ( 44) is the lumpability of model M 1 . Let r be a reward that is constant within every macro-state C x = ξ -1 (x), x ∈ X , and denote this common value by r(x), x ∈ X :

r(n) = r(x), ∀n ∈ C x .
Let w t and w t be as before, the t-period rewards for original model and model M 1 , and let v t be the t-period reward for model M 1 :

v t+1 (x) = r(x)+ ∅ =A⊂I λ A v t (d A (x))+ I k=1 µx k v t (r k (x))+   ν - ∅ =A⊂I λ A -µ I k=1 x k   v t (x), (46)
for all x ∈ X , t ≥ 0, where v 0 (x) := 0, x ∈ X . Now by Proposition 2, for all x ∈ X :

w t (n) = v t (x), ∀n ∈ C x . (47)
This property allows us to consider the relations for model M 1 , instead of M 1 : relations ( 44) and (47) imply that to show [START_REF] Van Dijk | Simple bounds and monotonicity results for finite multiserver exponential tandem queues[END_REF], it is sufficient to show that for all x ∈ X and for all A ⊂ I, A = ∅, functions v t satisfy:

v t (x - i∈A e i ) + (|A| -1)v t (x) ≥ i∈A v t (x -e i ), t ≥ 0. (48)
Using Proposition 8 in Appendix A, relation ( 48) is equivalent to the supermodularity of functions v t , t ≥ 0.

Theorem 4. Let r : E → R be a reward function that is constant within every macro-state C x = ξ -1 (x), x ∈ X , and r : X → R the corresponding aggregated reward function:

r(x) = r(n), x ∈ X , n ∈ C x .
Let w t : E → R and v t : X → R, t ≥ 0 be the t-period rewards respectively for model M 0 and model M 1 . If the reward function r is supermodular, then:

w t (n) ≥ v t (x), n ∈ C x , t ≥ 0,
and the mean stationary rewards satisfy:

a = lim t→∞ w t (n) t ≥ lim t→∞ v t (x) t = a.
Proof. Let r : E → R be a reward function that satisfies the assumptions of the theorem. Then Proposition 3 and Proposition 8 (in Appendix A) imply that the t-period rewards v t , t ≥ 0, for model M 1 satisfy (48). By (47), this is equivalent to [START_REF] Van Houtum | Bounds for performance characteristics: a systematic approach via cost structures[END_REF] and to (45). The result now follows from Theorem 2 and its Corollary 1.

It remains for us to show that the aggregate order fill rate is indeed a function satisfying the conditions of Theorem 4. The aggregated order fill rate is a linear combination of order fill rates for individual demand streams. Denote by OFR A : E → {0, 1} the order fill rate for the demand stream for set A ⊂ I, A = ∅:

OFR A (n) = i∈A 1 {ξi(n)<Si} , n ∈ E.
Reward function OFR A is clearly constant within every macro-state and we will denote also by OFR A : X → {0, 1} its aggregated version:

OFR A (x) = i∈A 1 {xi<Si} , x ∈ X . (49)
We showed in Lemma 1, that this function is supermodular. The aggregated order fill rate is supermodular as an linear combination of supermodular functions. Thus Theorem 4 implies: Corollary 3. The order fill rates for any individual demand stream and the aggregate order fill rate in model M 1 give the lower bounds for the individual and aggregate order fill rates in model M 0 .

Conclusions.

We have established a new method to compare Markov chains: A generalization of the precedence relation method to sets of states, which we have shown to be compatible with aggregation. Our methodology preserves the core advantage of precedence relations: the modifications of the original model are easy to understand and allow intuitive interpretation.

Precedence relations for sets of states enable the construction of bounding chains by replacing one set of transitions with another set. As illustrated on the example of a queue with batch arrivals, this can be used to compare systems with arrivals that have the same mean but different variability. This is one illustration of how our technique can be applied to derive bounds by replacing a part of a system with a simplified version having the same mean behavior. We expect that there are many other examples. Note that this is not typically possible using some classical methods for Markov chain comparison, for example strong stochastic ordering or the classical precedence relation method.

By showing that our new method is compatible with strong aggregation, we facilitate the construction of bounding chains with a state space of significantly reduced cardinality. This much smaller chain can then be used to derive bounds on the reward function of interest. We have also shown how precedence relations can be established both on the original or the aggregated bounding chain; in some cases the latter may be much easier. For example, in the service tools problem, we have shown that the precedence relations we need to compare the two chains are equivalent to the supermodularity property of cumulated rewards for the aggregated chain.

We apply our new method to the general family of assemble-to-order systems in order to obtain two different types of new analytical results:

• Monotonicity of supermodular rewards (e.g. item fill rates, order fill rates) with respect to coupling in arrivals for ATO systems with state dependent replenishment times. This extends results in Dayanik et al. [START_REF] Dayanik | The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-to-order systems[END_REF] to more general ATO systems.

• Proving bounds for the service tools problem, conjectured by Vliegen and Van Houtum [45].

In addition, due to aggregation, our bounding model has a state space that is highly reduced compared to the original system: Its dimension is reduced from exponential to linear, with respect to the number of different types of tools.

Having established the validity of generalized precedence relations combined with aggregation, there are several promising directions for extensions. We list three possibilities.

Our generalization of precedence relations to sets of states can be, to some extent, compared with the generalization of strong stochastic order to other integral orders, such as convex or increasing convex order. One possible future research direction is to compare the method presented here with comparison techniques based on stochastic monotonicity and different integral stochastic orders. One could expect, for selected families of functions and under certain conditions, that the two methods would be equivalent. If true, such an equivalence could allow on one hand the definition of a model-driven family of functions for which the precedence relations hold, and on the other hand enable the use of arguments of integral stochastic orders that allow both steady-state and transient comparison of Markov chains.

We studied the service tools model in which all the returns have equal rates, which is a natural assumption in our application. It would be interesting to study a generalization of this model to allow different return rates. Note that the monotonicity property with respect to supermodularity (Section 5.3), used to establish bounds for the service tools model, holds for ATO-POS with general replenishment rates. On the other hand, aggregation of states requires equal return rates for all the states that are merged together. If we allow set-dependent return rates in the service tools model M 0 (denote by ν A the return rate for set A), then a lower bound is an ATO-POS model with return rates obtained by taking the minimal return rates over all subsets that contain a given item (i.e. µ i = min{ν A : i ∈ A}, ∀i). The proof of this statement involves both generalized and classical precedence relations and it is out of the scope of this paper. These bounds may still be used if the differences in rates are not too high. Otherwise, this will result in very loose bounds and it is reasonable to expect that more accurate bounding models need to be found.

Finally, our framework in Section 5 only considers ATO systems with exponential replenishment times. Another stream of literature includes ATO systems with deterministic replenishment times (Song [START_REF] Song | On the order fill rate in a multi-item, base-stock inventory system[END_REF], Song [START_REF]Order-based backorders and their implications in multi-iyem inventory systems[END_REF], Hoen et al. [START_REF] Hoen | A simple and accurate approximation for the order fill rates in lost-sales assemble-to-order systems[END_REF], Dogru et al. [START_REF] Dogru | A stochastic programming based inventory policy for assemble-to-order systems with applications to the w model[END_REF]). To cover this our framework needs to be extended, possibly to Erlang lead times, which can approach the deterministic in the limit.

[45] I. M. H. Vliegen and G. J. van Houtum, Approximate evaluation of order fill rates for an inventory system of service tools, International Journal of Production Economics 118 (2009), no. 1, 339-351.

[46] P. J. Xu, R. Allgor, and S. C. Graves, Benefits of reevaluating real-time order fulfillment decisions, Manufacturing & Service Operations Management 11 (2009), no. 2, 340-355.

A Supermodularity and its characterization.

Definition 2. Let (X , ) be a lattice and f a real function on X . Then f is said to be supermodular if f (x ∧ y) + f (x ∨ y) ≥ f (x) + f (y), ∀x, y ∈ X .

In Proposition 7 we will give a characterization of supermodularity for the case of a finitedimensional lattice. Without loss of generality, we will assume the set (X , ) to be a subset of (N n 0 , ≤), where ≤ denotes the usual componentwise partial order. In that case, the ∧ (meet) and ∨ (join) operators are defined componentwise: (x 1 , x 2 , . . . , x n ) ∧ (y 1 , y 2 , . . . , y n ) = (x 1 ∧ y 1 , x 2 ∧ y 2 , . . . , x n ∧ y n ), (x 1 , x 2 , . . . , x n ) ∨ (y 1 , y 2 , . . . , y n ) = (x 1 ∨ y 1 , x 2 ∨ y 2 , . . . , x n ∨ y n ), where x i ∧ y i = min{x i , y i } and x i ∨ y i = max{x i , y i }, for all i.

Before stating the proposition, we introduce some additional notation that will be used. We recall that e i , 1 ≤ i ≤ n, denotes the vector with all the coordinates equal to 0, except the coordinate i that is equal to 1. Similarly, we will denote by e A the vector with all the coordinates equal to 0, except the coordinates that belong to the set A: e A := i∈A e i , A ⊂ {1, . . . , n}, A = ∅, and e ∅ := (0, . . . , 0). Proposition 7. Let (X , ) be a subspace of (N n 0 , ≤) and f : X → R. The following statements are equivalent:

1. f is supermodular.

2. For all i, j ∈ {1, . . . , n}, i = j, and for all x ∈ X such that x + e i + e j ∈ X : f (x + e i + e j ) + f (x) ≥ f (x + e i ) + f (x + e j ).

(51)

3. For all A ⊂ {1, . . . , n}, A = ∅, and for all x ∈ X such that x + i∈A e i ∈ X :

f (x + i∈A e i ) + (|A| -1)f (x) ≥ i∈A f (x + e i ). (52)
4. For all K ∈ N, A 1 , . . . , A K ⊂ {1, . . . , n}, A k = ∅, 1 ≤ k ≤ K, A i ∩ A j = ∅, i = j, and x ∈ X such that x + K k=1 e A k ∈ X :

f (x + K k=1 e A k ) + (K -1)f (x) ≥ K k=1 f (x + e A k ). ( 53
)
Proof.

• 1 ⇒ 4. We will show this implication by induction on K. For K = 1 relation (53) is trivially satisfied: f (x + e A1 ) ≥ f (x + e A1 ). In order to better understand relation (53), we will write it explicitly also for K = 2: f (x + e A1 + e A2 ) + f (x) ≥ f (x + e A1 ) + f (x + e A2 ). This follows trivially from supermodularity of f , since: Using induction hypothesis for K, we obtain:

f (x + K+1 k=1 e A k ) + f (x) ≥ K k=1 f (x + e A k ) -(K -1)f (x) + f (x + e A K+1 ),
which is equivalent to (53) for K + 1:

f (x + K+1 k=1 e A k ) + Kf (x) ≥ K+1 k=1
f (x + e A k ).

• 4 ⇒ 3. Consider an arbitrary but fixed subset A ⊂ {1, . . . , n}, A = ∅, and a state x ∈ X such that x + i∈A e i ∈ X . Let K = |A|, and denote by i 1 , . . . , i K the elements of A. Define A k = {i k }, k = 1, . . . , K. Then A i ∩ A j = ∅, i = j. For this choice of sets, (53) becomes:

f (x + K k=1 e i k ) + (K -1)f (x) ≥ K k=1 f (x + e i k ),
which is precisely (52).

• 3 ⇒ 2 follows directly by taking A = {i, j}.

• 2 ⇒ 1. Consider two arbitrary states x, y ∈ X . Then y can be written as:

y = x + i∈A α i e i - i∈B β i e i ,
where A, B ⊂ {1, . . . , n} and A ∩ B = ∅. Denote by a = i∈A α i e i and b = i∈B β i e i . Then y = x + a -b and:

x ∧ y = x -b, and x ∨ y = x + a.

If B = ∅, then x ≤ y, which implies x ∧ y = x and x ∨ y = y. Therefore relation ( 50) is trivially satisfied. The case A = ∅ is similar. We consider now the non-trivial case where both A = ∅ and B = ∅. We will first show the following relation for arbitrary i ∈ B and j ∈ A: f (x + e j ) + f (x -β i e i ) ≥ f (x) + f (x + e j -β i e i ).

(54) Indeed, the above relation can be obtained by adding the following relations (relation (51) for states x -me i , 1 ≤ m ≤ β i ):

f (x + e j ) + f (x -e i ) ≥ f (x) + f (x + e j -e i )

f (x + e j -e i ) + f (x -2e i ) ≥ f (x -e i ) + f (x + e j -2e i ) . . . what we needed to show.

Proposition 7 remains valid if we replace the sign + by -in 2, 3, and 4, and the proof is similar. Thus we have: Proposition 8. Let (X , ) be a subspace of (N n 0 , ≤) and f : X → R. The following statements are equivalent:

1. f is supermodular.

2. For all i, j ∈ {1, . . . , n}, i = j, and for all x ∈ X such that x -e i -e j ∈ X : f (x -e i -e j ) + f (x) ≥ f (x -e i ) + f (x -e j ).

(57)

Figure 1 :

 1 Figure 1: Batch arrivals.

Figure 2 :

 2 Figure 2: Batch arrivals: redirecting transitions.

Figure 3 :

 3 Figure 3: Redirecting the sets of transitions.

Figure 4 :

 4 Figure 4: Partial redirection of transitions.

Figure 5 :

 5 Figure 5: A system with 2 item types.

Figure 7 :

 7 Figure 7: State space of ISA model.

Figure 8 :

 8 Figure 8: Markov chain for the service tools model for I = 2, S1 = 1 and S2 = 2.

1 Figure 9 :

 19 Figure 9: Markov chain M 1 for I = 2, S1 = 1 and S2 = 2. The blue (dashed) transitions represent the transitions of the original model that have been replaced by red (bold) transitions. To simplify the figures, the (uniformization) loops are not shown when they are not part of redirected transitions.

  (x + e A1 ) ∧ (x + e A2 ) = x and (x + e A1 ) ∨ (x + e A2 ) = x + e A1 + e A2 .Assume now that relation (53) holds for some K ≥ 2. For states x+ K k=1 e A k and x+e A K+1 we have:x + K k=1 e A k ∧ x + e A K+1 = x and x + K k=1 e A k ∨ x + e A K+1 = x + K+1 k=1 e A k .Thus supermodularity of f implies:f (x + K+1 k=1 e A k ) + f (x) ≥ f (x + K k=1 e A k ) + f (x + e A K+1).

f 1 l=1β 1 l=1β 2 l=1α

 112 (x + e j -(β i -1)e i ) + f (x -β i e i ) ≥ f (x -(β i -1)e i ) + f (x + e j -β i e i ) Denote the elements of B by B = {b 1 , b 2 , . . . , b |B| }. Then adding the following relations (obtained by applying relation (54) for k = 1, . . . , |B| to state x -k-1 l=1 β b l e b l , with i = b k and j ∈ A):f (x + e j ) + f (x -β b1 e b1 ) ≥ f (x) + f (x + e j -β b1 e b1 ) f (x + e j -β b1 e b1 ) + f (x -β b1 e b1 -β b2 e b2 ) ≥ f (x -β b1 e b1 ) + f (x + e j -β b1 e b1 -β b2 e b2 ). . .f (x + e j -|B|b l e b l ) + f (x -|B| l=1 β b l e b l ) ≥ f (x -|B|b l e b l ) + f (x + e j -|B| l=1 β b l e b l )and using b = |B| l=1 β b l e b l gives:f (x + e j ) + f (x -b) ≥ f (x) + f (x + e j -b), j ∈ A. (55)By adding the following equations (obtained by applying relation (55) for k = 0, . . . , α j -1 to state x + ke j , j ∈ A):f (x + e j ) + f (x -b) ≥ f (x) + f (x + e j -b) f (x + 2e j ) + f (x + e j -b) ≥ f (x + e j ) + f (x + 2e j -b) . . . f (x + α j e j ) + f (x + (α j -1)e j -b) ≥ f (x + (α j -1)e j ) + f (x + α j e j -b)we obtain:f (x + α j e j ) + f (x -b) ≥ f (x) + f (x + α j e j -b), j ∈ A.(56)Denote the elements of A by A = {a 1 , a 2 , . . . , a |A| }. Then adding the following relations (obtained by applying relation (56) for k = 1, . . . , |A| to state x + k-1 l=1 α a l e a l , with j = a k ): f (x + α a1 e a1 ) + f (x -b) ≥ f (x) + f (x + α a1 e a1 -b) f (x + 2 l=1 +f (x + α a1 e a1 -b)α a l e a l ) ≥ f (x + a l e a l -b) + f (x + α a1 e a1 ) . . . f (x + |A| l=1 α a l e a l ) + f (x + |A|-1 l=1 α a l e a l -b) ≥ f (x + |A|-1 l=1 α a l e a l ) + f (x + |A| l=1 α a l e a l -b) and using a = |A| l=1 α a l e a l gives: f (x + a) + f (x -b) ≥ f (x) + f (x + a -b),

Table 1 :

 1 A framework for ATO systems. considerable amount of time for large problem instances. Iravani et al. ([START_REF] Iravani | On assemble-to-order systems with flexible customers[END_REF]) extended the model of Song et al. by allowing customers to be selective and/or flexible. Selective customers are not going to buy the final product unless all key items are available. If all key-items are available, the product is bought, even if other non-key-items are missing. Flexible customers are willing to substitute one or more non-available items by an item that is available. The ATO-TOS model is a special case of the model ofIravani et al. ([13]).

	Returns/Out-of-stock	1. TOS	2. POS	3. Backordering
	A. Individual	Song et al. [37]	Song et al. [37]	Song et al. [37]
		Kapuscinski et al. [14]	Dayanik et al. [8]	Song and Yao [38]
		Iravani et al. [13]	Iravani et al. [13]	Lu et al. [25]
				Lu et al. [26]
				Lu [23]
	B. Joint	Kelly [15]	Vliegen and Van Houtum [45]

Table 2 :

 2 Counterexample for ATO-TOS model.

	2,3	0.5	0.5
		Output
		Situation 1 Situation 2
	OFR 1	-	0.8549
	OFR 2	0.7132	0.7159
	OFR 3	-	0.8549
	OFR 1,3	0.7662	0.7143
	OFR 1,2,3	0.5521	0.5362
	OFR	0.6862	0.7153

Acknowledgments

This work is partially supported by ANR project Blanc SMS, Prins Bernhard Cultuurfonds, and Beta Research School for Operations Management and Logistics.

3. For all A ⊂ {1, . . . , n}, A = ∅, and for all x ∈ X such that x -i∈A e i ∈ X :

(58)