Order Fill Rate, Leadtime Variability, and Advance ... - Semantic Scholar

Order Fill Rate, Leadtime Variability, and Advance Demand Information in an Assemble-to-Order System Yingdong Lu ∗ IBM Research Division, T. J. Watson Research Center Yorktown Height, NY 10598 Tel. (914) 945-3738, e-mail: [email protected] Jing-Sheng Song† Graduate School of Management University of California, Irvine, CA 92697 Tel. (949) 824-2482, e-mail: [email protected] David D. Yao‡ Department of Industrial Engineering and Operations Research Columbia University, New York, NY 10027 Tel. (212) 854-2934, e-mail: [email protected] April 2001

Abstract We study an assemble-to-order system with stochastic leadtimes for component replenishment. There are multiple product types, of which orders arrive at the system following batch Poisson processes. Base-stock policies are used to control component inventories. We analyze the system as a set of queues driven by a common, multiclass batch Poisson input, and derive the joint queue-length distribution. The result leads to simple, closed-form expressions of the first two moments, in particular the covariances, which capture the dependence structure of the system. Based on the joint distribution and the moments, we derive easy-to-compute approximations and bounds for the order fulfillment performance measures. We also examine the impact of demand and leadtime variability, and investigate the value of advance demand information. Subject classifications: Inventory/production: multi-item; operating characteristics; stochastic models. Queues: approximations; batch; multichannel; correlation.

∗

Most of this research was carried out while the author was a post doctoral researcher at the Graduate School of Management of UC Irvine, supported in part by NSF grants DMI-9896339 and DMI-0084922. † Supported in part by NSF grants DMI-9896339 and DMI-0084922. ‡ Supported in part by NSF grants ECS-9705392 and DMI-0085124. Part of this author’s research was undertaken while he was on leave at the Chinese University of Hong Kong, Dept. of Systems Engineering and Engineering Management, and supported by a Direct Grant from CUHK and by HK/RGC grant CUHK4376/99E.

1

Introduction

In the movement towards improved supply chain management, more and more enterprises have adopted the assemble-to-order (ATO) system, which is a hybrid between make-to-stock at the component (subassembly) level and assemble-to-order at the end-product (final assembly) level. With today’s advanced information technology to easily obtain and exchange information, the ATO system appears to be an ideal business process that enables mass customization and quick response. The usual inventory-service tradeoff, a key factor in supply chain management, becomes even more prominent in the ATO system, as each customer order typically involves a large number of components, and the stockout of any component will cause a delay in supplying the order. Hence, it is critically important to characterize the order fill rate and the necessary inventory investment. The analysis, however, is difficult as the production-inventory dynamics among the components are highly correlated, driven by a common demand stream. It is the objective of this paper to develop models and tools that facilitate this analysis and generate insights to important system design issues. We consider an ATO system supporting multiple types of demand, which arrive at the system following batch Poisson processes. The inventory of each component is controlled by a base-stock policy, and the replenishment leadtimes are i.i.d. (independent and identically distributed) random variables. This results in a set of M X /G/∞ queues, driven by a common multi-class arrival stream. This model extends the single-product ATO system studied in Song and Yao [23]. In the multiproduct setting, the key challenge is that in addition to the assembly feature, which requires the simultaneous availability of several components in order to fulfill one demand, there is the distributional characteristic of each component appearing on the bill-of-material (BOM) of several product types – the so called component commonality. Hence, not surprisingly, the approaches used here are mostly different from those in [23]. In addition, we study a few new topics, such as the response-time based order fill rate and related insights. Specifically, we derive the joint queue-length distribution in terms of its generating function. The exact analysis sheds light on the dependence structure of the system, which facilitates the understanding of several managerial issues, such as the effect of product structures, the effect of demand and leadtime variablity, and the value of advance demand information. It also allows us to obtain various moments — mean, variance and covariance — in simple, explicit formulas, which, in turn, lead to several approximation schemes for key performance measures that are otherwise computationally intractable. A brief overview of the related literature is in order. Studies of ATO systems differ quite

1

substantially in the detailed modeling assumptions and approaches. For example, Hausman, Lee and Zhang [9] and Zhang [26] study periodic-review (discrete-time) models with multivariate normal demand and constant component replenishment leadtimes. Song [20, 21] studies continuous-review models with multivariate compound Poisson demand and deterministic leadtimes. Song et al. [22], Glasserman and Wang [7] and Wang [24] also consider multivariate (compound) Poisson demand but the supply process for each component is capacitated and modeled as a single-server queue. In contrast, in our study the supply submodel is a parallel processing system, as reflected in the infinite-server queues. This is appropriate when components are procured from exogenous suppliers, whose capacity is deemed ample relative to the orders placed from the ATO system in question. Gallien and Wein [8] use the same leadtime model as ours here, but focus on a single demand stream and assume order synchronization. Cheung and Hausman [4] also assume i.i.d. leadtimes in the context of a repair shop. Their multivariate Poisson demand model is a special case of ours (the unit demand case). They use a combination of order synchronization and disaggregation in their analysis. The rest of the paper is organized as follows. The formal model description is presented in §2, followed by the derivation in §3 of the generating function for the joint queue-length distribution, the moments, and the analysis on leadtime and batch-size variability. Several approximation schemes and bounds are developed in §4. Response-time based order fill rates are studied §5, and the value of advance demand information discussed in §6. Numerical results are presented in §7, and, finally, brief concluding remarks is summarized in §8.

2

Model Description

We consider an inventory system of m different components from which multiple types of products can be assembled upon customers’ requests. See Figure 1. [INSERT FIGURE 1 HERE] Let I = {1, 2, · · · , m} denote the set of all component indices. Customer orders arrive at the system following a stationary Poisson process, denoted {A(t), t ≥ 0}, with rate λ. Each order may require several components in different amounts simultaneously. For any subset of components K ⊆ I, we say an order is of type K if it consists positive units of component in K and 0 units in I \K. We assume that there is a fixed probability q K that an order is of type K,




K

q K = 1. Thus,

the type-K order stream forms a compound Poisson process with rate λK = q K λ. A type-K order K requests QK j units of component j ∈ K, and QK = (Qj , j ∈ K) follows a known discrete probability

distribution. We assume that each order’s type is independent of the other orders’ types and of all

2

other events. We denote K to be the set of all demand types, that is, K = {K ⊂ I : q K > 0}. Note that K is not necessarily the set of all possible subsets of I. This demand model is flexible in modeling a variety of situations. One special case is a system with only one demand “type”: All demands request m components in stock, however, the units of each component requested can be different with Qi ≥ 1 being the batch size for i, i = 1, ..., m. Here, if Qi are constant, the system is equivalent to an assembly system with a single final product, where (Qi ) represent the BOM. If Qi are random variables, the model represents a system with a variety of products, in which the variety is characterized by the number of units of each component contained in the product. Each combination of Qi ’s corresponds to the BOM of a specific product. Another special case is a system in which each customer demand requires one and only one component. In this case, the system reduces to m independent single-component systems. Below we shall use the results from this simple system to bound the performance of the original, general system. For each component i, let Ki denote the family of subsets of K that contain i. It is clear that the demand process for component i forms a compound Poisson process with rate λi =




K∈Ki

λK = qi λ

and batch size Qi , the mixture of QK i for all K ∈ Ki . (Throughout the paper we use subscripts to indicate the component types and superscripts the demand types.) Demands are filled on a first-come-first-served (FCFS) basis. If there is enough on-hand inventory for all the components required by a demand upon its arrival, the demand is filled immediately. In other words, we assume that the time to assemble the components into the end-product is negligible. We also assume complete backlogging for demands that cannot be filled immediately. When a demand arrives and some of its required components are in stock but others are not, we can either ship the in-stock components or put them aside as committed inventory. However, a demand is considered backlogged unless it can be satisfied completely. When there are backorders, they are also filled on a FCFS basis. The inventory of each component is controlled by an independent base-stock policy, with si := the base-stock level for component i. That is, upon each demand arrival, if the inventory position (i.e., the on-hand inventory plus onorder position minus backorders) of component i is less than si , then order up to si ; otherwise, do not order. This type of policy is known to be optimal when there are no economies of scale in replenishment and when the component demands are independent. When component demands are correlated as in our case, the optimal policy should be a coordinated policy. Unfortunately, the precise form of the optimal policy is still unknown. Hence, the independent base-stock policies are widely used 3

in practice, partly because of its simplicity, partly because it provides a benchmark on how much inventory is needed in order to provide a certain service level. Since we follow a base-stock policy and the demand arrives in batches, each component replenishment order comprises several units. We assume that, for each unit of component i, the replenishment leadtimes are i.i.d. random variables with a common cumulative distribution Gi . Let Li denote the generic random variable with distribution Gi and mean E[Li ] = i . Denote Gci = 1 − Gi . Assume the leadtimes are independent among the components; that is, Li is independent of Lj for any i = j. (Alternatively, one can assume that each order experiences Li units of leadtime, which is independent of the order size. We will discuss this assumption in the Concluding Remarks.) The performance measure of primary interest is the immediate availability of all components needed for an arriving demand, or the “off-shelf” fill rate. (We will study a more general performance measure – the response-time based order fill rate – in §5.) It turns out that this measure can be obtained directly from the steady-state joint distribution of outstanding orders, as explained in the following. For any time t, let Ii (t) be the net inventory of component i at t, and Xi (t) be the number of outstanding orders of component i at t. Then, by the nature of the base-stock control, we have: Ii (t) = si − Xi (t),

i = 1, ..., m.

(1)

In the next section, we will show that the joint distribution of (X1 (t), ..., Xm (t)) has a steady-state limit (X1 , ..., Xm ). Therefore, we have the steady-state net inventory Ii = si − Xi ,

i = 1, ..., m.

From the property of Poisson Arrival See Time Average, we have fi := off-shelf fill rate of component i = P[Xi + Qi ≤ si ], fK

:=

f

= =

off-shelf fill rate of type-K demand = P[Xi + QK i ≤ si , ∀i ∈ K], average (over all demand types) off-shelf fill rate



qK f K .

K

Thus, to study the performance of the ATO system, the key is to obtain the joint distribution of (X1 , ..., Xm ).

4

3

Performance Analysis

3.1

Generating Functions and Moments

We now derive the joint distribution of the outstanding orders vector X(t) = (X1 (t), ..., Xm (t)), and its steady-state limit. By the nature of the base-stock policy, each type-K order triggers replenishment orders of all the components in subset K. Thus, the arrival processes to the supply subsystem are the same as the external demand processes to the ATO system. Also, for each component i, the number of outstanding orders Xi (t) is exactly the number of jobs in service in an MiQi /Gi /∞ queue with Poisson arrival rate λi and batch size Qi . Since each order generates simultaneous arrivals at several component queues, the m queues are not independent. See Figure 2. [INSERT FIGURE 2 HERE] On the other hand, the m queues are dependent only through the common arrival streams; their operations are otherwise independent, in particular the processing times are independent across the queues. In other words, given the number of demand arrivals up to t, Xi (t)’s are independent of one another. This fact is crucial, as it allows us to proceed, after conditioning upon the number of arrivals, as if we are analyzing independent systems. This way, we can derive (see details in the Appendix) the joint distribution of the outstanding orders, in terms of its generating function, as summarized in the following proposition. Proposition 1 Let X = (X1 , ..., Xm ) follows the limiting distribution of X(t) = (X1 (t), ..., Xm (t)) as t → ∞. Then, the generating function of X is: ψ(z1 , ..., zm ) := E[ 

= exp

 K∈K

K

λ

 ∞ 0

m  Xj

j=1

(ψ

zj ]

QK



(G1 (u) +

z1 Gc1 (u), ..., Gm (u)

+

zm Gcm (u))

− 1)du .

(2)

Below, we shall focus on the random vector X. In the special case of unit arrivals, i.e., Qi ≡ 1, for all i = 1, ..., m, the generating function takes the following form: ψ(z1 , ..., zm ) = exp[



K∈K

K

λ

 ∞  0

(

j∈K

[Gj (u) + zj Gcj (u)] − 1)du],

which corresponds to a multivariate Poisson distribution. In particular, ψ(1, ..., zi , ..., 1) = exp(



K∈Ki

5

λK i (zi − 1)).

(3)




That is, for each i, Xi is a Poisson random variable with parameter λi i = (

K∈Ki

λK )i , as

expected. The multivariate Poisson distribution implies that each Xi can be expressed as a finite sum of independent univariate Poisson distributions, see, e.g., [12]. Specifically, letting N (a) denote a Poisson random variable with parameter (mean) a, we can write 

Xi =

N (λK θSK )

Si,S⊆K

(4)

where all the Poisson variables under the summation are independent, and θSK =

 ∞ 

j∈S

Gcj (u)



Gj (u)du.

(5)

j∈K\S

In principle, we can compute the joint distribution of X = (X1 , ..., Xm ) through the above expressions. This, however, involves 2m − 1 (independent) Poisson variates. Therefore, even in the case of unit arrivals, the exact evaluation of the joint distribution would still be impractical even for moderately large m. Next, we derive the means, variances and covariances, making use of the generating function in Proposition 1. Let Kij denote the family of subsets that contain both i and j. Then, the first two moments of Xj can be derived as follows: µj := E[Xj ] =

 ∂ |zj =1 ψ(z1 , ..., zm ) = j λK E(QK j ), ∂zj K∈K

(6)

j

E[Xj2 ] = (

∂2 ∂ + )|z =1 ψ(z1 , ..., zm ) 2 ∂zj ∂zj j 2

= E(Xj ) + E (Xj ) +



K

λ

K∈Kj

2 [E((QK j ) )

−

E(QK j )]

 ∞ 0

[Gcj (u)]2 du,

(7)

and for i = j, E[Xi Xj ] =

∂2 |z =1 ψ(z1 , ..., zm ) ∂zj ∂zi j

= E(Xi )E(Xj ) +



K∈Kij

K λK E(QK i Qj )

 ∞ 0

Gci (u)Gcj (u)du.

(8)

Hence, from (7), we have σj2 := Var[Xj ] =

E(Xj ) +

 K∈Kj

2 K λK [E((QK j ) ) − E(Qj )]

6

 ∞ 0

[Gcj (u)]2 du.

(9)

Note that since E(QK j ) ≥ 1 for all K ∈ Kj , we have 2 K 2 K E[(QK j ) ] ≥ [E(Qj )] ≥ E(Qj ),

and hence σj2 ≥ µj . Note that σj2 = µj if and only if we have unit arrivals, i.e., Qi ≡ 1, which is consistent with the fact that Xj follows a Poisson distribution. From (8), we have, for i = j, 

σij := Cov[Xi , Xj ] =

K

λ

K∈Kij

K E(QK i Qj )

 ∞ 0

Gci (u)Gcj (u)du.

(10)

Note that σij = 0 if Kij is empty. That is, Xi and Xj are independent if there is no demand types that require both i and j. In other words, the correlation of the queues is solely induced by the common arrivals. Indeed, (10) also indicates that, as long as there exists common arrivals to both queue i and queue j, Xi and Xj are always positively correlated. With (9) and (10), it is easy to calculate the correlation coefficients ρij =

σij . σi σj

It is insightful to examine the explicit formulas in some special cases. In particular, assuming unit demand, we obtain


ρij Here λij =




=

K∈Kij

K∈Kij




i

λK

K∈Ki

∞ 0

λK

Gci (u)Gcj (u)du



j




K∈Kj

λK



=

λij

λi λj

  ∞ 0

Gci (u)Gcj (u)du

i j



.

(11)

λK is the aggregate demand rate for the component pair (i, j), and λi and λj

are the aggregate demand rate for components i and j, respectively. Notice that

∞ 0



∞ c c Gci (u)Gcj (u)du 0 Gi (u)Gj (u)du

≤ 1. =   ∞ c ∞ c i j G (u)du G (u)du i

0

0

j

So, λij . ρij ≤

λi λj If the proportion of the demand types that require both i and j is very small, so that λij  λi and λij  λj . Then the correlation between Xi and Xj is negligible, and the two queues can hence be treated as independent. On the other hand, if the components i and j are always demanded together, such as in the pure assembly (i.e., single-product) system, then λi = λj = λij . In this case, ∞

ρij =

0

Gci (u)Gcj (u)du

. i j 7

Thus, although the correlation is induced by the common arrival, the level of correlation is independent of the demand rate. The correlation is imperfect; its level depends on the leadtime distributions. Suppose, in addition, that leadtimes Li and Lj are exponentially distributed, then  ∞ 0

Gci (u)Gcj (u)du =

and

i j , i + j

(12)



ρij =

i j 1  = . i + j i /j + j /i

If both leadtimes have the same distribution, so that i = j , then ρij = 0.5. Thus, the imperfect correlation is a consequence of the statistical variabilities of the leadtimes. Now, let us examine what happens if there is no statistical variability in leadtimes. That is, suppose both Li and Lj are deterministic. In this case,  ∞ 0

Gci (u)Gcj (u)du = min{i , j },

and ρij =

(13)

min{i , j } . max{i , j }

Thus, the difference between the lengths of the leadtimes determines the level of correlation. The bigger the difference is, the smaller the correlation becomes. If i = j , then the correlation is perfect, as expected, because the two queues are completely synchronized due to the common arrival stream and the same fixed service times. Finally, we remark that the generating function in Proposition 1 can also be used to obtain all other higher moments of X. In particular, for any subset {i1 , ..., ik } of {1, 2, ..., m}, let U be any ¯ = {i1 , ..., ik } \ U . We have, subset of {i1 , ..., ik } with cardinality |U | and denote U E[Xi1 ...Xik ] =



k 



(λK )k−j+1

¯ h∈U

K:{i1 ,...,ik }⊆K j=1 U ⊂{i1 ,...,ik },|U |=j

3.2



h E[



h∈U

QK h ]

 ∞  0

h∈U

Gch (u)du.

(14)

The Impact of Variability

We now compare the original system with another system, where the processing times of queue i, are i.i.d. random variables with mean i (same as in the original system) and a distribution function ˜ i satisfying G  ∞ x

Gci (u)du

≥

 ∞ x

˜ ci (u)du, G 8

i = 1, ..., m;

(15)

˜c = 1 − G ˜ i . All other aspects of the two systems are identical. for any x ≥ 0. Like Gc , here G i ˜ i be the random variables associated with the distributions Gi and G ˜ i , respectively. Let Li and L ˜ i in the sense of the convex Then, the above condition means that Li is more “variable” than L ˜ i ; i.e., Ef (Li ) ≥ Ef (L ˜ i ) for all convex function f (·) (see, e.g., [17]), and ordering, denoted Li ≥cx L hence in particular, ˜ i − i ]2 = Var[L ˜ i ]. Var[Li ] = E[Li − i ]2 ≥ E[L Similarly, we can address the issue of variability associated with the batch-size distribution. ˜ K ), for component i of type-K product. Assume QK Suppose the new system has a batch size of (Q i

i

˜ K in the same sense of convex ordering as above, whereas all other aspects is more variable than Q i of the new system are the same as the original system. Proposition 2 Use “tilde” to denote the new system that has reduced variability in either leadtimes or batch sizes, in the sense of convex ordering as detailed above. We have

(i) E[

j∈T



Xj ] ≤ E[

j∈T

˜ j ], for any T ⊆ {1, ..., m}. X

˜ and X ≤uo X, ˜ where ≥lo (≤uo ) denotes lower-orthant (upper-orthant) ordering. (ii) X ≥lo X That is, for any (x1 , ..., xm ), the following hold: ˜ 1 ≤ x1 , ..., X ˜m ≤ xm ], P[X1 ≤ x1 , ..., Xm ≤ xm ] ≤ P[X and ˜ 1 ≥ x1 , ..., X ˜m ≥ xm ], P[X1 ≥ x1 , ..., Xm ≥ xm ] ≤ P[X respectively, for lower-orthant and upper-orthant orderings. (iii) f K ≤ f˜K for all K ∈ K and f ≤ f˜. The proof of the above proposition is presented in the Appendix. Part (iii) of the proposition indicates that variability degrades order fill rates. Part (i) of the proposition implies, in particular, ˜i, X ˜j ]. Cov[Xi , Xj ] ≤ Cov[X That is, reducing the variability of leadtimes or batch sizes will result in a higher correlation among the queue lengths of outstanding jobs. Combining the above inequality with the covariance expression in (10), and letting the leadtimes in the modified system be (deterministic) constants, we have, Cov[Xi , Xj ] ≤ min{i , j }

 K∈Kij

9

K λK E(QK i Qj ).

Next, for any pair of queues i and j, suppose in the original system the leadtimes Li and Lj are exponentially distributed, whereas the modified system has constant leadtimes as above, with i = E(Li ) and j = E(Lj ); and assume, without loss of generality, i ≤ j . Then, applying (12) and (13), we have

j Cov[Xi , Xj ] = , ˜i , X ˜j ] i + j Cov[X

which ranges from 0.5 to 1, depending on how close is i to j .

4

Approximations

In this section, we develop computational efficient approximations and bounds for the joint distribution of X, which, in turn, form easy-to-compute approximations for the order fill rates.

4.1

Factorized Normal Approximation

It is quite natural to use multivariate Normal distribution N (µ, Σ) to approximate (X1 , ..., Xm ), since the first two moments, including the covariances, are explicitly derived in §3.1. Here, µ := (µ1 , ..., µm ) with µj given in (6), and the covariance matrix Σ := (σij ) with σij given in (10). More specifically, P[X1 = x1 , ..., Xm = xm ] ≈

1 (2π)m/2 |Σ|1/2

1 exp{− (x − µ)Σ−1 (x − µ) }, 2

where x = (x1 , ..., xm ) (a row vector), and (x − µ) is the transpose of (x − µ). To compute the cumulative distribution, we have, P[X1 ≤ x1 , ..., Xm ≤ xm ] ≈

 x1 −∞

...

 xm −∞

1 (2π)m/2 |Σ|1/2

1 exp{− (y − µ)Σ−1 (y − µ) }dy. 2

(16)

Unfortunately, however, because of the multiple integrals involved, even the state-of-the-art computational packages of multivariate normal distributions are still limited to rather low dimensions, e.g., m ≤ 7; refer to, for example, Drezner [5]. For higher dimensions, we propose a further approximation on the multivariate normal distribution using the following approach. Applying H¨ older’s inequality, we have

1/2

K K 2 E[QK i Qj ] ≤ E[(Qi ) ]

and

 ∞ 0

Gci (u)Gcj (u)du

≤

 ∞ 0

1/2

2 E[(QK j ) ]

(Gci (u))2 du

1/2  ∞

·

0

,

(Gcj (u))2 du

1/2

.

Hence, denote 

ηi := 



K∈Ki

1/2 K

λ

2  E((QK i ) )

·

10

 ∞ 0

(Gci (u))2 du

1/2

,

(17)

and similarly denote ηj . Then, we can obtain a factorized (between i and j) upper bound on the covariance: σij =



K

λ

K∈Ki,j

≤



K

λ

K∈Ki,j

K E(QK i Qj )

 ∞ 0

Gci (u)Gcj (u)du

 2 K 2 1/2 E[(QK ) ]E[(Q ) ] i j

·

 ∞ 0

(Gci (u))2 du

 ∞ 0

(Gcj (u))2 du

1/2

≤ ηi ηj . The second inequality above takes into account Kij ⊆ Ki ∩ Kj . Note that we have ηi < σi , which follows from comparing (17) against (9), taking into account the following: 

E(Xi ) =

K

λ

K∈Ki



>

K∈Ki

E(QK i )

λK E(QK i )

 ∞

Gci (u)du

0

 ∞ 0

(Gci (u))2 du.

Similarly, we have ηj < σj . Now, let V := (V1 , .., Vm ) be a multivariate normal distribution with the same mean µ = (µ1 , ..., µm ) as before, but with the covariance σij replaced by ηi ηj . Then, we can write 

Vi = µi + ηi Z0 +

σi2 − ηi2 Zi ,

i = 1, ..., m;

(18)

where Z0 , Z1 , ..., Zm are i.i.d. standard normal variates. Clearly, V has the same marginal distributions as N (µ, Σ). But due to its special form in (18), we can avoid evaluating the multiple integrals in (16); instead, we only need to evaluate the following single integral: P[X1 ≤ x1 , ..., Xm ≤ xm ] ≈ P[V1 ≤ x1 , ..., Vm ≤ xm ]  ∞

=

−∞

φ(u)

m 

xi − µi − ηi u Φ(  )du, σi2 − ηi2 i=1

(19)

where φ and Φ are the normal density and distribution functions, respectively. We shall refer to the above as the factorized Normal approximation.

4.2

Pairwise Approximation

Next, we develop a second type of approximation, which we call pairwise approximation. It is based on the following result (proved in the Appendix):

11

Proposition 3 For any (x1 , ..., xm ) ≥ 0, we have P[X1 ≤ x1 , ..., Xm ≤ xm ] ≥

m 

P[Xj ≤ xj ].

(20)

j=1

More generally, letting (S1 , ...Sk ) be a partition of I = {1, ..., m}, we have P[X1 ≤ x1 , ..., Xm ≤ xm ] ≥

k 

P[Xj ≤ xj , j ∈ S ].

(21)

=1

In particular, Proposition 3 implies the following: fK ≥

 i∈K

P[Xi + QK i ≤ si ] ≥



P[Xi + Qi ≤ si ] =

i∈K



fi .

(22)

i∈K

Thus, the product of the component fill rates is a lower bound of the order fill rate. Note that in two earlier works, Song [21] and Song and Yao [23], the inequality in (20) has been established for special cases of the system considered here (specifically, those with constant leadtimes [21] and with single-class demand [23]). Our second approximation scheme is simply to apply the relation in (21) to a pairwise partition of I, and use the lower bound as an approximation for the joint distribution. (In the case of m being odd, let the last subset Sk be a singleton set.) To evaluate this lower bound, we can first replace each pair by a bivariate normal distribution, and then evaluate the product of k bivariate normal distributions. The only remaining question is how do we form the k pairs. In our numerical experiments, we have tried to take the best (largest) lower bound among all possible ways (≤ m2 /2) of forming pairs. The difference (among these different ways), however, appears to be insignificant. Hence, our recommendation is to form pairs in decreasing order of the covariances. That is, (i) Set I0 = I; set  = 1. (ii) Set S = {i , j } such that (i , j ) = arg max{σij , (i, j) ∈ I0 }. (iii) Stop, if I0 is either a singleton set or a null set. Otherwise, set I0 ← I0 \ {i , j },  ←  + 1; and go to (ii). Finally, in the case of unit arrivals, X = (X1 , ..., Xm ) follows a multivariate Poisson distribution, as discussed in §3.1. The pairwise approximation described above applies readily here, the only modification is: instead of evaluating bivariate normal distributions, here we evaluate bivariate Poisson distributions. In particular, we still form the pairs in decreasing order of their covariances. 12

The factorized approach could, in principle, also be applied here. However, our numerical experiments have shown that the quality of this approximation is not quite as good as in the case of normal distributions. The reason is because while we can always write a normal variate X as X = µ+σZ with Z being the standard normal, there is no analogous expression for a Poisson variate √ N . (If anything, we can write N = λ + λZ; but this goes back to the normal approximation.)

5

Response-Time Based Order Fill Rate

A key service measure is order response time: the time length between the point when the order is placed and the point when the product is received by the customer. This response time usually consists of two parts — the time needed to have all the components ready for assembly, denoted W ; and the time for outbound logistics, denoted W o , which includes time for processing orders, assembly time, and transportation time for delivering the order. What affects service is mainly the W part — when any component inventory has a stockout, the delay could be substantial, whereas W o is often nearly deterministic. Therefore, we shall focus on W , and denote by f K (w) the probability of having all the components i ∈ K ready within w units of time (i.e., W ≤ w). (In this sense, w is the allowed “time window” to fill the order.) In other words, f K (w) is the probability that the overall response time to type-K orders is under w + W o units of time. Suppose a type-K order arrives at time τ . Consider component i ∈ K. Define Di (t) to be the cumulative demand for component i in [0, t]. Let Di (t, t + u] = Di (t + u) − Di (t) for any t, u ≥ 0. Then, Xi (τ ) + Di (τ, τ + w] − Xi (τ + w) is the total number of departures from queue i in (τ, τ + w]. Since this is the total addition to the component i’s inventory between τ and τ + w, the order that arrives at τ can be filled by τ + w if and only if the net inventory at τ plus this addition is nonnegative. That is, Ii (τ ) + [Xi (τ ) + Di (τ, τ + w] − Xi (τ + w)] ≥ 0. Using (1), the above inequality simplifies to: Xi (τ + w) − Di (τ, τ + w] ≤ si ,

i ∈ K.

(23)

Observe that Xi (τ + w), the number of outstanding orders of component i at time τ + w, has the following decomposition: K

Xi (τ + w) =

Xiw (τ )

+

Qi  n=1

1{Lni > w} + Xi (τ, τ + w].

13

(24)

The three terms on the right hand side represent the number of units still in queue i (at τ + w) that correspond to the arrivals, respectively, in (0, τ ), at τ , and in (τ, τ + w]. 1{·} denotes the indicator function, and Lni is an independent copy of Li representing the leadtime for the nth unit of QK i . Note that Xi (τ, τ + w] := Xi (τ + w) − Xi (τ ) and Xiw is different from Xi . Moreover, the three terms on the right hand side of (24) are independent of one another, due to Poisson arrivals and the infinite servers in each queue. Thus, combining (23)and (24), we know the demand at τ can be supplied by τ + w if and only if K

Xiw (τ )

+ Xi (τ, τ + w] − Di (τ, τ + w] ≤ si −

Qi  n=1

1{Lni > w},

i ∈ K.

(25)

When τ → ∞, X w (τ ) has a limit X w , which has the following generating function (see derivation in the Appendix): exp{



(λS

 ∞ 0

S∈K

S

[ψ Q (G1 (u + w) + z1 Gc1 (u + w), ..., Gm (u + w) + zm Gcm (u + w)) − 1]du}.

(26)

Denote: Yiw := Di (τ, τ + w] − Xi (τ, τ + w], which represents the number of jobs that arrived in (τ, τ + w] but also departed by τ + w. Denote Y w := (Yiw ), which has the following generating function (derived in the Appendix): exp{



S

 w

(λ (

S∈K

0

S

[ψ Q (Gc1 (u) + z1 G1 (u), ..., Gcm (u) + zm Gm (u)) − 1]du)}.

(27)

Moreover, X w and Y w are independent, because of the independent-increment property of the Poisson arrival process. Denote Yi := Xiw − Yiw ,

i = 1, ..., m.

Then, Y = (Y1 , ..., Ym ) has the following probability generating function: ψY (z1 , ..., zm ) =

exp{



(λS

S∈K  w

−λS (

0

 ∞ 0

S

[ψ Q (G1 (u + w) + z1 Gc1 (u + w), ..., Gm (u + w) + zm Gcm (u + w)) − 1]du

S

[ψ Q (Gc1 (u) + G1 (u)/z1 , ..., Gcm (u) + Gm (u)/zm ) − 1]du)}.

(28)

Applying (23), we have f K (w) := order fill rate of type-K demand within time window w K

=

P[Yi +

Qi  n=1

1{Lni > w} ≤ si , 14

∀i ∈ K].

(29)

(Kruse [14] obtained a similar result in a single-component, unit-demand, base-stock inventory sys
QK i

n n=1 1{Li

tem.) The random variable

> w} is a mixture of geometric distributions. Conditioning

c on QK i = k, it has the geometric distribution with parameters k (the number of trials) and Gi (w)

(the success rate). Note that when w = 0, the second term in (28) vanishes, reducing the generating function to the one in (2). Thus, when w = 0, Y reduces to X, and f K (0) coincides with the off-shelf fill rate f K , as expected. Since Yi is independent of QK i and Li as explained earlier, to derive the above fill rate, the key is the distribution function of Yi , which follows from its generating function derived above. While this can only be done numerically in general, we can easily derive the first two moments in closed form from the generating function in (28). It is straightforward to see that all the approximation schemes and bounds developed in the last section automatically apply here. In particular, we have 





∂ |z =1 ψY (z1 , ..., zm ) =  λS E(QSi ) (i − w) , ∂zi i S∈K

E[Yi ] =

i



Var[Yi ] =



 S∈Ki



+



λS E(QSi )

 S∈Ki

 w 0

 ∞

G(u)du + 

λS (E(QSi )2 − E(QSi ))

w

 w 0

Gci (u)du



[Gi (u)]2 du +

 ∞ w



[Gci (u)]2 du ,

and for i = j, 

Cov[Yi , Yj ] = 



 S∈Kij

λS E[QSi QSj ]

 w 0

 ∞

Gi (u)Gj (u)du +

w



Gci (u)Gcj (u)du .

Similar to the properties of X discussed in §3, here, too, we have Var[Yi ] ≥ E[Yi ] for all i. Also, for any pair i and j, Yi and Yj are positively correlated as long as there exists a demand type that requires both i and j. Otherwise, they are independent. However, in the unit demand case, E[Yi ] = λi (i − w) = λi and Var[Yi ] = λi

 ∞ w

Gci (u)du −

 w 0

 ∞

G(u)du +

w



 w 0

Gi (u)du ,

Gci (u)du



.

So, for w > 0, Var[Yi ] > E[Yi ], i.e., Yi does not have a Poisson distribution.

15

6

The Value of Advance Demand Information

With the availability of advanced information technology, there is an increasing trend in industry to try to induce early customer order information so to help streamline operations. This strategy is believed to be effective in delivering what customers want, and when and where they want, in a cost efficient manner. For one thing, knowing early what customers want helps inventory planning by ordering and stocking the right items in right amount. In this section, we apply the analytical approach developed in the last section to gain insight into the value of early demand information on order fulfillment performance. Suppose a customer sends a signal of his/her order w periods earlier than when he/she actually places the order. That is, w is the amount of time in advance that we obtain the information on what a customer wants to order. As long as we have all the items requested by the order available during the window w, we will be able to fill the order when it actually arrives at the system. We adopt the following revised base-stock policy: we keep the original base-stock levels. However, we place a replenishment order of a component whenever there is a demand signal that requests this item. This revised base-stock policy is shown to be optimal in various single-item systems by Hariharan and Zipkin [10], who refer to w as demand leadtime. For any demand type K, note that f K (w), the fill rate with a time window w studies in the last section, measures exactly the probability we can satisfy the order at the moment the demand occurs, provided that we obtained the demand information w periods earlier and we use the revised base-stock policy. On the other hand, the probability of satisfying an order at the moment the order occurs without any advance demand information is the immediate fill rate f K (0). It is trivially held that f K (w) > f K (0). That is, for the same level of inventory investment, the advance demand information definitely benefits the order fulfillment measure. To assess the value of advance demand information more precisely, it is interesting to compare f K (w)

with the immediate fill rate of an modified system, in which the leadtimes are truncated by

w. This allows us to relate the impact of the advance demand information to the effect of leadtime reduction. The importance of this kind of comparison has been eloquently argued in [10]. Specifically, consider a modified system that differs from the original one only in the leadtime ˆ i = [Li − w]+ . Note that for any u ≥ 0, distributions, with component i’s leadtime being L ˆ i > u) = P([Li − w]+ > u) = P(Li > u + w) = Gci (u + w). P(L ˆm ) be the steady-state vector of outstanding orders in the new system. Then, ˆ = (X ˆ 1 , ..., X Let X ˆ has the same distribution as X w = (X w , ..., X w ). examining (2) and (26) we conclude that X 1

16

2

Thus, the immediate order fill rate for the new system is: K

fˆK (0) = P[Xiw +

Qi  n=0

ˆ ni > 0} ≤ si , ∀i ∈ K]. 1{L

(30)

ˆ i . On the other hand, for the original system, from (29), we ˆ n is an independent copy of L Here, L i have K

K

f (w) =

P[Xiw

−

Yiw

+

Qi  n=1

1{Lni > w} ≤ si ,

∀i ∈ K].

(31)

ˆ i > 0, and Y w ≥ 0, we have Since Li > w if and only if L i fˆK (0) ≤ f K (w).

(32)

Proposition 4 Consider a multi-product ATO system with stochastic leadtimes, the original system, along with a modified system where the leadtime for each component is reduced by a common amount w (or down to zero if the leadtime is less than w). The immediate (“off-shelf”) fill rate in the modified system is no greater than the fill rate within a response-time window of w in the original system. Thus, inequality in (32) indicates that knowing demand in advance (by w time units) is more effective, in terms of order fill rate, than reducing the supply leadtimes of components (by the same amount w). In fact, advance demand information can do even better. The inequality in (32) holds when both systems maintain the same base-stock levels (for component inventory). Suppose fˆK (0) in (32) corresponds to the optimal base-stock levels, in terms of maximizing fˆK (0), in the modified system, and let the original system follow the same base-stock levels. Then, (32) holds. We can then optimize the base-stock levels in the original system to further increase the right hand side of (32). It is interesting to note that Song [21] shows that if all leadtimes are deterministic, then fˆK (0) = f K (w). The equation suggests that obtaining customer order information w periods earlier has precisely the same effect on order fulfillment performance as reducing the component replenishment leadtimes by w. This is consistent with the discussions in [10]. Thus, the inequality (32) implies that the impact of advance demand information in systems with random leadtimes is in general greater than that in systems with deterministic leadtimes. One plausible interpretation is that reducing leadtimes by certain amount as suggested in the modified system does not really reduce the leadtime variability, whereas the advance demand information provides the same degree of flexibility in order fulfillment as in systems with deterministic leadtimes. 17

7

Numerical Results

In this section we illustrate the effectiveness of the approximations through numerical results. We focus on an ATO system of desktop computers studied in Cheng et al. [3], using a slightly simplified bill of materials (BOM) (through combining certain common components of different products). Specifically, a family of product consists 6 different types of products, which are assembled from a set of 14 different components. The components are obtained from external suppliers with different supply leadtimes. Table 1 provides the list of all the components and their average leadtimes, along with the BOM of each product. i 1 2 3 4 5 6 7 8 9 10 11 12 13 14

Description Shell Processor: 450MHz Processor: 500MHz Processor: 550MHz Processor: 600MHz 64 MB Memory 6.8GB Hard Drive 13.5GB Hard Drive Hard Drive (Common Parts) Software 1.0 Software 2.0 CDROM VGA Card Ethernet Card

i 7 12 12 12 12 15 18 18 8 4 4 10 6 10

K:

1 1 1 1 1 1 1 1 -

2 1 1 1 1 1 1 1 -

3 1 1 1 1 1 1 1 1 1

4 1 1 1 1 1 1 1 -

5 1 1 1 1 1 1 1 -

6 1 1 1 1 1 -

Table 1: Average component leadtimes and the BOM of each product. The aggregate demand process for the six products forms a Poisson process with mean λ = 50. The proportions of the six products are: q 1 = 0.10, q 2 = 0.15, q 3 = 0.4, q 4 = 0.2, q 5 = 0.1, q 6 = 0.05. We consider four different sets of base-stock levels at the components as summarized in Table 2. The data (base-stock levels) in Table 2 are chosen such that the off-shelf fill rates (lower bounds) are achieved at about 98%, 95%, 90% and 80% (respectively for the four groups) for all 6 products. To demonstrate the impact of leadtime variability on order fill rates, we consider both exponential distribution with mean i and Erlang distribution with parameters (i /2, 2).

18

Component Level 1 Level 2 Level 3 Level 4

1 52 50 48 47

2 13 12 12 10

3 17 16 16 14

4 67 62 60 60

5 9 8 7 6

6 96 94 90 86

7 52 44 42 46

8 87 83 83 78

9 67 59 59 60

10 30 29 27 25

11 10 10 10 9

12 36 35 34 30

13 54 46 42 44

14 36 35 34 30

Table 2: Four different groups of base-stock levels Table 3 reports the numerical results from four different approximations along with simulation. The results are also plotted in Figure 3. The four approximation schemes are: (i) Lower Bound (Marginal), or LB1, which is the lower bound in (20); (ii) Lower Bound (Pairwise), or LB2, which is the pairwise lower bound in (21), using bivariate Poisson distributions; (iii) Normal Approximation, or NA1, which is the factorized normal approximation given in (19); and (iv) Normal Approximation 2, or NA2, which uses bivariate normal distribution in the pairwise lower bound. [INSERT TABLE 3 AND FIGURE 3 HERE] From these results, we observe the following: • The factorized normal approximation overestimates the fill rates in all cases. • The pairwise normal approximation also overestimates the fill rate in most (but not all) cases. • Not surprisingly, the pairwise lower bound outperforms the marginal lower bound in all cases. • The pairwise normal approximation also outperforms the factored normal approximation in all cases. • For high fill rates (0.95 and above), the pairwise lower bound is better than the pairwise normal approximation. In fact, in this range, it is the best approximation among four. In this case, the marginal lower bound is also very reliable. • For fill rates around the range of 0.80-0.85, the pairwise normal approximation appears to be the best among all four. In the range of 0.90-0.95, the two pairwise bounds are the best. Comparing the results with different leadtime distributions, we also observe that the pairwise approximations definitely captures more of the probabilistic natures of the system such as the effect of leadtime distributions, while the marginal lower bound does not. Since the factorized normal approximation and the marginal lower bound require much less computational effort, and their performance is reasonably good, they can be used in practice as quick estimates. Also, one of them is a (proven) lower bound, while the other offers (in most cases) an upper bound. 19

[INSERT FIGURE 4 HERE] In Figure 4 we plot the simulation results reported in Table 3. The purpose is to illustrate the effect of leadtime variability. Clearly, the leadtime variability has more pronounced impact on fill rates in the lower ranges.

8

Concluding Remarks

We have studied a multiproduct ATO system with i.i.d. leadtimes and independent base-stock controls at the components. Modeling the demand process as a multivariate compound Poisson process, we have derived the joint distribution of outstanding orders in terms of generating function, along with close-form formulas for the first and second moments of the joint distribution. Examples were given that shed light on how product structures and leadtime distributions affect correlation. We have also established several comparison results, based on stochastic order relations, which yield insight to the impact of variability on system performance. In addition, we discussed the value of advance demand information. Finally, based on the exact analysis, we have developed several approximation schemes and bounds for the joint distribution and related service measures such as order fill rates. As mentioned in §2, sometimes the multiple units in a replenishment order can be treated as a single, indivisible unit at the supply system. Therefore, the entire batch experiences the same leadtime. In this case, the arrival process to the supply subsystem is a multivariate Poisson process (not compounded). Our analysis in §3 is still valid, with batch sizes Qi = 1 for all i. As shown in (3), the joint congestion level at the supply subsystem, X, has a multivariate Poisson distribution; and the net inventory of component i can be expressed as si −

Xi  n=1

Qni ,

where Qni is an independent copy of Qi . Thus, with the joint distributions of X and Q, one can compute the type-K order fill rate, P[

Xi 

n=1

Qni + QK i ≤ si , ∀i ∈ K];

and bounds based on similar ideas can also be developed. Further work is needed to develop computational efficient means to solve problems that seek the optimal inventory-service tradeoff in the type of multiproduct ATO systems studied here. Along this line, several structural properties of the system need to be further explored. These are the focus of our follow-up studies. 20

Appendix Proof of Proposition 1 Consider any given subset of components K ∈ K. Let AK (t) be the cumulative number of typeK arrivals up to time t. Then AK (t) follows a Poisson distribution with parameter λK , and its generating function is: E[z A

K (t)

K t(z−1)

] = eλ

.

At any time t, the number of jobs in queue j associated with type-K arrivals can be expressed as follows: K,n

XjK (t)

AK (t) Qj

=

 

n=1 k=1

1kjn (t),

(33)

k denotes the nth independent copy of QK where QK,n j , and 1jn (t) is the indicator function that the j

arrival to queue j generated by the nth demand is still in service by time t, n = 1, ..., AK (t). Let τnK be the arrival epoch of the nth demand, n = 1, ..., AK (t). The event that the arrival to queue j generated by nth arrival remains in service by time t is equivalent to that the processing time for this job is longer than t − τn . That is, 

1kjn (t)

1 w.p. Gcj (t − τnK ) 0 w.p. Gj (t − τnK )

=

k = 1, ..., QK,n . j

Hence, for any zj ≥ 0, j = 1, ..., m, we have 

E

m  XjK (t)

zj

j=1 AK (t)

 

=

n=1 j∈K



|AK (t), τ1K , τ2K , ..., τAKK (t) , QK,1 , ..., QK,A

K (t)



K,n

[Gj (t − τn ) + zj Gcj (t − τnK )]Qj

Observe that given AK (t), τnK follows independent uniform distribution over [0, t]. Hence, 

E

m  XjK (t)

j=1

=

  t 1

t

0

zj

ψ

QK



|AK (t)

(G1 (u) +

z1 Gc1 (u), ..., Gm (u)

+

zm Gcm (u))du

AK (t)

,

(34)

K

where ψ Q (z1 , ..., zm ) denotes the generating function of QK . Taking expectations on both sides, with respect to AK (t), we obtain the generating function for X K (t): ψtK (z1 , ..., zm ) 

= exp λK

 t 0



K

[ψ Q (G1 (u) + z1 Gc1 (u), ..., Gm (u) + zm Gcm (u)) − 1)]du . 21

Letting t → ∞, we have ψ K (z1 , ..., zm ) 

= exp λK

 ∞ 0



K

[ψ Q (G1 (u) + z1 Gc1 (u), ..., Gm (u) + zm Gcm (u))] − 1)]du .

(35)

K ) be the random vector corresponding to the above generating function ψ K . Let X K = (X1K , ..., Xm

We have X K (t) → X K in distribution. As the arrival processes, {AK (t), K ∈ K}, are independent across K, the generating function of X(t) = (X1 (t), ..., Xm (t)) is: ψt (z1 , ..., zm ) = exp[



λK

 t

K∈K

0

K

(ψ Q (G1 (u) + z1 Gc1 (u), ..., Gm (u) + zm Gcm (u)) − 1)du].

(36)

Letting t → ∞, we obtain the desired expression in (2).

Proof of Proposition 2 To show (i), note that the inequality in (15) implies  ∞ 0

j∈T

Gcj (u)du

≤

 ∞ 0

j∈T

˜ cj (u)du. G

This can be directly verified via integration by parts. (Also see Theorem 7.3 of [2]). From (14), we

observe that E[

j∈T

Xj ] is a linear combination, with positive coefficients, of

∞ 0

j∈T 

Gcj (u)du,

for T  ⊆ T ; hence, the desired inequality. To show the orderings in (ii) we first consider unit arrivals. In this case, X follows a multivariate Poisson distribution, which is known ([12, 13]) to be the limiting distribution of a sequence of multivariate binomial variates {Yn }. The generating function of Yn is: 

ψn (z1 , ..., zm ) = 

1 

···

a1 =0

1  am =0

n am  pa1 ...am z1a1 ...zm

with npa1 ...am → wa1 ...am , as n → ∞, where the limits wa1 ...am relate to Xi ’s via the following: E[Xi ] =

1 

...

a1 =0

E[Xi Xj ] =

1  a1 =0

...

1  ai =1

1 

...

ai =1

...

1  aj =1

1  am =0

...

1  am =0

22

wa1 ...ai−1 1ai+1 ...am ,

wa1 ...ai−1 1ai+1 ...aj−1 1aj+1 ...am .

(37)

In general, for any T ⊆ {1, ..., m}, we have, E[



1 

Xj ] =

j∈T

a1 =δ1,T

where

1 

...



am =δm,T

wa1 ...am ,

0 j∈T . 1 j ∈ Tc

δj,T =

In particular, we can choose pa1 ...am = wa1 ...am /n for each Yn . Also, from (37), we see that Yn is the summation of independent and identically distributed multivariate Bernoulli variates Zin , which have the following common probability generating function, ψZin (z1 , ..., zm ) =

1 

···

a1 =0

1  am =0

am pa1 ...am z1a1 ...zm

The above argument applies to the “tilde” system as well. Hence, similarly, we have multivariate binomial variates Y˜n , and the corresponding Z˜in . From the inequality established in Part (i), we can see that for each T , 1  a1 =δ1,T

1 

...

1 

pa1 ...am ≥

am =δm,T

...

a1 =δ1,T

1 

p˜a1 ..am .

am =δm,T

Notice that for the multivariate Bernoulli variates, Zin and Z˜in , the above implies that for any vector u = (u1 , ..., um ), with ui ∈ {0, 1}, we have P[Zin ≥ ui , ∀i] ≥ P[Z˜in ≥ ui , ∀i]. That is, Zin ≥uo Z˜in . Since this ordering is preserved under i.i.d. summation (e.g., Theorem 4.B.10 ˜ in [17]), we have Yn ≥uo Y˜n . Letting n → ∞, we obtain the desired X ≥uo X. K Next, consider batch arrivals. To be specific, write P[QK j = rj , j = 1, ..., m} = pr1 ...rm . Each

vector (r1 , ..., rm ) corresponds to a fixed bill of materials. We can decompose the arrival process for type-K demands into several independent Poisson arrival processes, one for each configuration (r1 , .., rm ). In other words, X =




r1 ,...,rm

X (r1 ,...,rm) , where X (r1 ,...rm) is the number of outstanding

orders generated by corresponding demand and has the following generating function: 

ψ (r1 ,...,rm) (z1 , ..., zm ) = E  

= exp λpr1 ...rm

m (r ,...,rm )  Xj 1

j=1

 ∞  m 0

[

zj

  

(Gj (u) + zj Gcj (u))rj − 1]du .

j=1

23

On the other hand, consider another generating function: 

rj  ∞  m 

exp λpr1 ...rm

0

[



(Gj (u) + zjk Gcj (u)) − 1]du .

(38)

j=1 k=1

It is not difficult to see that the random vector that corresponds to the above generating function is the limit of a sequence of multivariate binomial variates. This implies that the generating functions of the latter will converge to the generating function in (38). Letting zjk = zj for all k = 1, ..., rj , we obtained the desired result. As for the lower-orthant order, notice that for any T ⊂ {1, ..., m},  

δ1,T

δm,T





...

a1 =0

am =0

n

pa1 ...am  = ψn (δ1,T , ..., δm,T ),

meanwhile, ψn (δ1,T , ..., δm,T ) → ψ(δ1,T , ..., δm,T ),

n → ∞.

Now, the convex order implies that, for any U ⊂ {1, ..., m}, the following holds:  ∞  0



Gi (u) − 1 du ≤

i∈U

 ∞  0



˜ i (u) − 1 du. G

i∈U

(The above can again be directly verified via integration by parts; or see Theorem 7.4 of [2].) Hence, we have, for any T ,



ψ(δ1,T , ..., δm,T ) = exp

˜ 1,T , ..., δm,T ) = exp ≤ ψ(δ



K

λ

KK



K

λ

 ∞  0

i∈T ∩K

0

i∈T ∩K

 ∞ 

KK





Gi (u) − 1 du 



˜ i (u) − 1 du . G

Consequently, for n large enough, we have, ψn (δ1,T , ..., δm,T ) ≤ ψ˜n (δ1,T , ..., δm,T ); or equivalently,  

δ1,T



a1 =0

δm,T

...



am =0

n



pa1 ...am  ≤ 

δ1,T



a1 =0

δm,T

...



am =0

n

p˜a1 ...am  ,

˜ by the same argument as in the case of which implies that Z n ≤lo Z˜ n , and hence X ≤lo X upper-orthant ordering. Finally, for the comparison result with respect to the variability of batch sizes, the above argument also applies, once we observe that a convex order on the batch size, along with (14), will ˜ j ], which is the starting point of the above argument. guarantee E[ j∈T Xj ] ≤ E[ j∈T X 24

Proof of Proposition 3 Consider demand type K ∈ K, with the corresponding arrival process being AK (t); the associated outstanding orders of component j, XjK (t) follows the expression in (33). Now, let AK j (t), j = 1, ..., m, be i.i.d. copies of AK (t), and write ˆ K (t) = X j

K,n AK j (t) Qj

 

n=1 k=1

1kjn (t).

K (t)). Note the following: for any (a , ..., a ), we have ˆ K (t), ..., X ˆm ˆ K (t) = (X Write X 1 m 1 m 

P[AK (t) ≤ a1 , ..., AK (t) ≤ am ] = P[AK (t) ≤ min{aj }] ≥ j

j=1

P[AK j (t) ≤ aj ].

That is, K (AK (t), ..., AK (t)) ≤lo (AK 1 (t), ..., Am (t)).

(39)

For each j = 1, ..., m, define a function QK,n

hj (y) := 1[

y  j 

n=1 k=1

1kjn (t) ≤ ai ].

Obviously, hj (y) is non-negative and decreasing (i.e., non-increasing) in y. Hence, following Theorem 4.G.1 of [17], we know the lower-orthant ordering in (39) implies 

E

m 





hj (AK (t)) ≥ E 

j=1

m 

j=1

  hi (AK j (t)) =

m  j=1

E[hi (AK j (t))].

ˆ K (t). Since the lower-orthant ordering is preserved under i.i.d. summation That is, X K (t) ≤lo X ˆ which is the desired inequality in (20). (Here X ˆ (resp. X) is and taking limit, we have X ≤lo X, ˆ K is ˆ K (t) (resp. X K (t)) over K ∈ K and with t → ∞. Also note that X the i.i.d. summation of X j equal in distribution to XjK , and independent across j = 1, ..., m.) The inequality in (21) is similarly proved, with k i.i.d. copies of the arrival process, one for each partition.

Derivation of (26) and (27) To start with, consider the simple case of a single demand type which requests one unit from each of the m components. In this case, Di (t) = A(t), the number of arrivals during the interval [0, t), for all j and t. Conditioning on A(τ ), i.e., the number of arrivals during the interval [0, τ ) and

25

A(τ, τ + w], the number of arrivals during (τ, τ + w], we have 

E

 m  Xiw (τ )+Xiw (τ,τ +w]−A(τ,τ +w]   j=1

zj

 

 m  Xiw (τ )+Xiw (τ,τ +w]−A(τ,τ +w] = E E  zj |A(τ ), A(τ, τ + w) j=1

 

= E E 

m  Xjw (τ )

j=1

zj





|A(τ ) E

m  Xiw (τ,τ +w]−A(τ,τ +w]  j=1

zj



|A(τ, τ + w] .

Similar to (34), we have 

E

m  Xjw (τ )

j=1

zj





|A(τ ) = 

 τ  m 0

[

A(τ )

[Gj (u + w) + zj Gcj (u + w)] − 1]du

.

j=1

and 

E

m  A(τ,τ +w]−Xi (τ,τ +w]  j=1

zj



|A(τ, τ + w]

 A(τ,τ +w]  w  m c = [ [Gj (u) + zj Gj (u)] − 1]du 0

j=1

Also, from the independent increment property of Poisson process, we know that these two terms are independent. Thus, 

E

 m  w w Xi (τ )+Xi (τ,τ +w]−A(τ,τ +w]   j=1

zj

= exp(λ

 τ  m 0

[

[Gj (u + w) +

j=1

zj Gcj (u

+ w)] − 1]du − λ

 w  m 0

[

[Gcj (u) + Gj (u)/zj ] − 1]du).

j=1

Letting τ → ∞, we obtain 

E

 m  Xiw −Yiw   j=1

zj

= exp(λ

 ∞  m 0

[

[Gj (u + w) +

j=1

zj Gcj (u

+ w)] − 1]du − λ

 w  m 0

[

[Gcj (u) + Gj (u)/zj ] − 1]du).

j=1

We can now follow the steps in deriving (36) to extend the above expression to the multiproduct case, namely to (26) and (27).

26

References [1] Bahadur, R. R., A Representation of the Joint Distribution of Response of n Dichotomous Items, in Studies in Item Analysis and Predictions, Solomon, J. ed., 158-168, Stanford University Press, Stanford, CA, 1961. [2] Barlow, R.E. and Proschan, F., Statistical Theory of Reliability and Life Testing: Probability Models, Hol, Rinehart and Winston, Inc, 1975. [3] Cheng, F., Ettl, M., Lin, G. and Yao, D.D., Inventory-Service Optimization in Configure-to-Order Systems: From Machine-Type Models to Building Blocks, IBM Research Report, RC 21781(98077), June 2000. [4] Cheung, K.L. and Hausman, W., Multiple Failures in a Multi-Item Spare Inventory Model, IIE Transactions, 27 (1995), 171-180. [5] Drezner, Z. Computation of the Multivariate Normal Integral, ACM Transactions on Mathematical Software, 18, (1992)470-480. [6] Fendick, K. W., An Asymptotic Exact Decomposition of Coupled Brownian Motion, J. Appl. Prob., 30 (1993), 819-834. [7] Glasserman, P. and Wang, Y., Leadtime-Inventory tradeoffs in Assemble-to-Order Systems, Operations Research, 46 (1998), 858-871. [8] Gallien, J. and Wein, L., A Simple and Effective Component Procurement Policy for Stochastic Assembly Systems, working paper, MIT Sloan School/Operations Research Center, 1998. [9] Hausman, W.H., Lee, H.L. and Zhang, A.X., Order Response Time Reliability in a Multi-Item Inventory System, European J. of Operational Research, 109 (1998), 646-659. [10] Hariharan, R. and Zipkin, P., Customer-order Information, Leadtimes, and Inventories, Management Science, 41 (1995), 1599-1607. [11] Johnson, N. L. & Kotz, K. Distributions in Statistics: Continuous Multivariate Distributions. John Wiley & Sons, Inc. (1972) [12] Johnson, N. L. Kotz, S. and Balakrishnan, N., Discrete Multivariate Distribution, Wiley, New York, 1996. [13] Kawamura, K., The Structure of Multivariate Poisson Distribution, Kodai Mathematical Journal, 2 (1979) 337-345. [14] Kruse, W.K., Waiting Time in an S − 1, S Inventory System with Arbitrarily Distributed Lead Times, Operations Research, 28 (1980), 348-352. [15] Matsunawa, T., Poisson Distribution, Encyclopedia pf Statistical Sciences. 7 20-25. New York, Wiley. [16] Ross, S.M., Stochastic Processes, 2nd ed., Wiley, New York, 1996. [17] Shaked, M. and Shanthikumar, J.G., Stochastic Orders and Their Applications, Academic Press, New York, 1994. [18] Shanthikumar, J.G. and Yao, D.D., Strong Stochastic Convexity and Its Applications, J. Appl. Prob. 28 (1991), 131-145. 27

[19] Shanthikumar, J.G. and Yao, D.D., Bivariate Characterization of Some Stochastic Order Relations, Adv. Appl. Prob. 23 (1991), 642-659. [20] Song, J.S., Order-Based Backorders and Their Implications in Multi-Item Inventory Systems, preprint. [21] Song, J.S., On the Order Fill Rate in a Multi-Item, Base-Stock Inventory System, Operations Research, 46 (1998), 831-845. [22] Song, J.S., Xu, S. and Liu, B. Order Fulfillment Performance Measures in an Assemblyto-Order System with Stochastic Leadtimes, Operations Research, 47 (1999), 131-149. [23] Song, J.S. and Yao, D.D., Performance Analysis and Optimization in Assemble-to-Order Systems with Random Leadtimes, Operations Research, forthcoming. [24] Wang, Y. Near-Optimal Base-Stock Policies in Assemble-to-Order Systems under Service Levels Requirements, working paper, MIT Sloan School, 1999. [25] Wolff, R., Stochastic Modeling and the Theory of Queues, Prentice Hall, Englewood Cliffs, New Jersey, 1989. [26] Zhang, A.X. Demand Fulfillment Rates in an Assemble-to-Order System with Multiple Products and Dependent Demands, Production and Operations Management, 6 (1997), 309-324.

28

Figure 1. The Assemble-to-Order System Suppliers

Components (Items) L1

L2

1

Products (Customer demands)

Q121 Q122

2 QK2 QK m

Lm

λ12

m

QK1

λΚ

Backorders

Figure 2. The Supply Subsystem Suppliers X1

Xm

X2

Q

Q

12 1

λ12

12

2

QK2

QK m

QK1

λΚ Arrivals (Replenishment Orders)

Figure 3. Approximations vs. Simulation The Desktop Example Exponential, Level 1

Erlang, Level1

0.9950 0.9900 0.9850 0.9800 0.9750 0.9700 0.9650 0.9600

Sim LB1 LB2 NA1 NA2 1

2

3

4

5

0.9900 0.9850 0.9800 0.9750 0.9700 0.9650

6

Sim LB1 LB2 NA1 NA2 1

2

Exponential, Level 2

Sim LB1 LB2 NA1 NA2 2

3

4

5

6

LB1 LB2 NA1 NA2 2

3

Sim

5

6

0.9400

LB1

0.9200

0.9000

LB2

0.9000

0.8800

NA1

0.8800

NA2

0.8600 4

4

Erlang, Level 3

0.9200

3

6

Sim

1

0.9400

2

5

0.9900 0.9800 0.9700 0.9600 0.9500 0.9400 0.9300

Exponential, Level 3

1

4

Erlang, Level 2

0.9900 0.9800 0.9700 0.9600 0.9500 0.9400 0.9300 1

3

5

Sim LB1 LB2 NA1

0.8600

6

NA2 1

2

Expone ential, Leve l 4

3

4

5

6

Erlang, Lev el 4

0.9000

0.8600 0.8400 0.8200 0.8000 0.7800 0.7600 0.7400 0.7200 1

2

3

4

5

6

Sim LB1 LB2 NA1

0.8500

Sim

0.8000

LB1

NA2

0.7000

LB2

0.7500

NA1 NA2

1

2

3

4

5

6

Figure 4. Effect of Leadtime Variability

1.00

Level 1, Exp

0.95

Fill rate

Level 2, Exp Level 3, Exp

0.90

Level 4, Exp Level 1, Erl 0.85

Level 2, Erl Level 3, Erl

0.80

Level 4, Erl

0.75 1

2

3

4

Product

5

6