Two perishable inventory systems with one-way substitution

Motivated by the ABO issue of the blood bank system, in which the portions stored have constant shelf life, we consider two subsystems of perishable inventory. The two Perishable Inventory Subsystems-PIS A and PIS B, are correlated to each other through a one-way substitution of demands. Specifically, the input streams and the demand streams applied to each subsystem are four Poisson processes, which are independent of one another. However, if the shelf of PIS A (blood of type O) is empty of items, an arriving demand of type A is unsatisfied, since demand of type A cannot be satisfied by an item of type B (blood portions of type AB), but if the shelf of PIS B is empty of items, an arriving demand of type B is applied to PIS A, since demands of type B can be satisfied by both types. This one-way substitution of the issuing policy generates for PIS A a modulated Poisson demand process operating in a two-state non-Markovian environment. The performance analysis of PIS B is known from previous work. Thus, in this study we focus on the marginal performance analysis of PIS A. Based on a fluid formulation and a Markovian approximation for the one-way substitution demand process, we develop a unified approach to efficiently and accurately approximate the performance of the PIS A. The effectiveness of the approach is investigated by extensive numerical experiments.


Introduction
We consider a stochastic input-output inventory system composed of two correlated Perishable Inventory Systems-PIS A and PIS B. Items of type A arrive in the shelf of PIS A according to a Poisson process with rate λ A and items of type B arrive in the shelf of PIS B according to a Poisson process with rate λ B . The shelf life of all items is a constant that without loss of generality is equal to 1. Demands of type A, that arrive according to a Poisson process with rate μ A apply to PIS A and demands of type B, that arrive according to a Poisson process with rate μ B , apply first to PIS B, but if PIS B is empty the application is endorsed to PIS A. The four Poisson streams are independent. However, since demands of type A can be satisfied only by items of type A, a demand of type A leaves unsatisfied if shelf A is empty. Unsymmetrically, demands of type B can be satisfied by either items of type B or by items of type A. We assume that the issuing policy in both subsystems is first-in-first-out (FIFO).
We are interested in the marginal performance analysis of each subsystem. The performance measures are the long-run averages for the number of items on shelf s A (s B ), the rate of item loss due to perishing A ( B ) and the rate of overall demand loss m. In fact, several versions of the marginal performance analysis of PIS B have already been carried out in previous work [for the analysis of the basic PIS B as described above, see Kaspi and Perry (1983)]. However, the performance analysis of PIS A appears to be new. To see the intricateness of a rigorous analysis, note that while the arrival process of items into PIS B is Poisson with rate λ B and the demand process is Poisson with rate μ B , the arrival process into PIS A is Poisson with rate λ A , but the demand process applied to PIS A is a modulated Poisson process operating in a two-state random environment, which is determined according to the environment status of PIS B. Namely, when shelf B is not empty, the demand process applied to PIS A is a Poisson process with rate μ A , but when shelf B is empty, the demand process applied to PIS A is a Poisson process with rate μ A + μ B . Also, while the time periods in which shelf B is empty are exponentially distributed (with parameter λ B ), the time periods in which shelf B is not empty are not exponential (but the law of these time periods can be computed), so that the random environment associated with the demand arrival process into PIS A is not Markovian. The latter fact makes a rigorous analysis of PIS A too complicated, if possible at all. In light of the intricateness of the input process into PIS A, which is a non-renewal process, it seems that it is not likely to be able to perform an exact analysis of the relevant performance measures of PIS A. Accordingly, for the analysis of PIS A we apply an accurate approximation based on the following approach (Osogami and Harchol-Balter 2006): Let U n be the length of the nth non-emptiness period in PIS B, which we call ON period. Clearly U 1 , U 2 , . . . are independent and identically distributed (i.i.d.) random variables and let U be the generic random variable of this sequence. Since we know how to compute the law of U we also know to compute its moments. Now take a random variablê U with a known phase-type distribution such that EU k = EÛ k for all k = 1, 2, . . . , n, for some predetermined n. Then intuitively, the model in which the original U is replaced bŷ U will be a good approximation to the original model and the approximation is expected to improve as n increases, see Osogami and Harchol-Balter (2006).
Our model is motivated by the ABO issue associated with the blood bank system. In practice, there are four types of blood-O, A, B and AB (in this study we assume a generic model of only two types). The blood portions arrive in accordance with four independent Poisson streams and are classified into the four categories (shelves)-O, A, B and AB. According to the formal statistics, more than 40 % of the population belong to category O and less than 5 % belong to category AB. Demands of type O can be satisfied by only blood portions of type O, but demands of type AB can be satisfied by either blood portions of type O or blood portions of type AB. It turns out that a shortage in blood portions of type O might have been a disaster, but a shortage in blood portions of type AB is important only from a managerial point of view. According to the formal standards and regulations, the maximum shelf life of all blood portions is 21 days. That means that after 21 days any blood portion cannot be used for transfusion. Then, after 1 unit of time (21 days), the blood portion is removed from the shelf where the removed items generate the so-called outdating process. Except for extreme urgency, it is natural to believe that the issuing policy of blood portions is FIFO. Namely, if the shelf is not empty, any arriving demand is satisfied by the oldest portion on the shelf. From a modeling point of view this means that the FIFO issuing policy can be used as a good approximation to reality.
Over the last three decades several review papers on inventory models with deterioration in the utility of the items have been published in the operations research literature. We indicate here only the comprehensive overviews (Nahmias 1982(Nahmias , 2011Karaesmen et al. 2009); the other reviews focus on other operations research issues and as a result, are not relevant to this study. The three aforementioned monographs refer to more than 200 works about perishable inventory models. Undoubtedly, the general literature about perishable inventory models is very rich, but the dominant component of the field aims at models that apply optimization and/or control; namely, models that study the behavior of the optimal values of certain decision variables. Apparently, in most of the operations research models there exists a controller who faces the problems of 'when to place an order' and 'how much to order'. However, there is a very large class of perishable inventory systems for which this type of well motivated problems is of irrelevance. These include models operating in the absence of ordering policies; that is, models that are run without controllers who face the above traditional problems of 'when to order' and 'how much to order'. For example, some blood bank systems with stochastic input (the arrival process of items) and stochastic output (the demand for those items and also outdatings of the items) might be regarded as appropriate applications of such models. Astonishingly, this natural stream of problems received only a sparse attention in the operations research literature. It appears that the inventory literature on perishable inventories with random input (and with the absence of control ordering policies) is not rich at all. The reason for that is probably the fact that such models are closely related to stochastic queueing models. As a result, these models focus generally on performance analysis, not on optimization. More precisely, these models may comprise two separated phases; the first phase is that of performance analysis. Then, the relevant measures and functionals that are found in the first phase lay the groundwork for the second phase, which is optimization. However, the type of control in the second phase is completely different from the traditional problems of when and how much to order, since naturally, the decision factors are different. Due to the intricateness of the stochastic model introduced in this study, we restrict our attention to the first phase of a prototype problem and focus on the performance analysis of a certain stochastic perishable inventory system. In the language of operations research applied to health care, this problem is called the ABO transfusion of blood problem (as mentioned above).
Every natural objective function will take into account at least three factors (or measures): the unsatisfied demands, the outdatings and the number of items on the shelf for the holding costs. In particular, in the blood bank application, unsatisfied demand might be a disaster. On the other hand, too many outdatings show that the system is run inefficiently and the controller will have to pay high holding costs if too many items wait on the shelf. Our analysis is based on a one-dimensional process, which is the age of the oldest item on the shelf, and once we know the steady state law of the age of the oldest item on the shelf, we are able to compute the rate of the unsatisfied demand process, the rate of the outdating process and the law of the number of items on the shelf in steady state (recall that the latter process is not a Markov process).
As indicated above, tractable PISs are mostly restricted to renewal arrival processes. Hence the current PIS with a non-renewal demand process motivates us to develop a new approach to analyze a single PIS under a pure Markovian setting. Accordingly, we relax the renewal property and exploit the Markov property by taking approximations when necessary, for tractability. The approach is a combination of two simple ideas of practical importance: (i) Approximate the non-renewal arrival process by a Markovian arrival process (cf. Neuts 1979); specifically, a Markov-modulated Poisson process, and (ii) Use a fluid formulation for the virtual waiting time process (e.g. Asmussen 2003, p. 308).
The rest of the paper is organized as follows. In Sect. 2, we present a fluid formulation for the PIS. We apply the fluid formulation to PIS B and immediately obtain a fluid model driven by a continuous time Markov chain (CTMC). In Sect. 3, we address PIS B and illustrate the proposed approach in detail. We re-derive several known results. We introduce explicit formulas for quantities of our interest, in particular, for the ON period. A direct application of the fluid formulation to PIS A results in a non-Markovian model. In Sect. 4, we propose to approximate the ON period by certain phase type distributed random variables by using the same moments of the ON period. We extend the state space of the driving CTMC. Then the approximate evaluation of PIS A is similar to that of PIS B. In Sect. 5, we present numerical experiments to investigate approximation errors. Finally, in the last section, we summarize our results and indicate possible directions for future research.

Fluid formulation
Consider A(t), the age of the oldest stock in a PIS at time t ≥ 0. By convention, let A(t) = 0 if the system is empty (out-of-stock) at time t. Illustrated in upper Fig. 1, the sample path of A(t) has a linear growth of rate 1 when the inventory is in stock. The linear growth is interpreted as the aging of the oldest stock. A downward jump occurs when the oldest stock is removed from the system. A removal of the oldest stock happens due to either (a) an arrival of demand (recall that in Sect. 1 we assume FIFO, namely the oldest stock is always used to satisfy a demand); or (b) perishing.
Let 0 ≤ S 1 ≤ S 2 ≤ . . . be all epochs that an item is removed. We assume that the sequence {S i } i=1,2,··· corresponds to a locally finite counting process, i.e. P{lim i→∞ S i < ∞} = 0. Clearly the size of the ith jump, i = 1, 2, . . ., is the inter-arrival time of the ith supply and the next, possibly truncated due to nonnegativity of A(t). We assume A(t) is right continuous. We shall analyze the age process {A(t)} t≥0 and show that the performance measures of our interest can be derived from the stationary distribution of the age process. Our approach is to couple the age process with a family of bivariate processes {(I r , A r )} parameterized by a single positive constant r . The motivation is easy to see by thinking of vertical jumps in a sample path as segments of infinite slope (r = ∞). We construct A r in such a way that the sample paths of A r converge to the sample paths of A as r → ∞. Under a Markovian setting, it turns out that the (I r , A r ) process is a piecewise deterministic Markov process (PDMP, cf. Davis 1984) with continuous (linear) sample paths, also known as fluid models. Thus we formulate and solve the PIS problem as a fluid model. A similar approach is also used in Liu and Kulkarni (2008) to analyze the busy period of a M/P H/1 queueing model with impatient customers. Although the probability laws can be derived directly for the A process within the PDMP framework, we introduce such a construction because it is a unified approach, which is directly amenable to computations. The construction, and our rationale to introduce it, shall become clear in the following sections. The numerical experiments in Sect. 5 also serve as a convincing support of the approach.
The idea is explained in Fig. 1, where A r is constructed from A by (i) extending instants when the oldest stock in the system perishes (so A reaches 1) with exponential times with mean 1/r (during which A stays 1) and (ii) replacing vertical jumps by jumps of the same size, but with slope r . Formally, for any r > 0, let us define process {A r (t)} t≥0 as the following transform of the A process.
where (a) X 0 , X 1 , . . . are i.i.d. exponential random variables with mean 1; (b) 1 {x} denotes the indicator function that equals 1 if condition x is true and 0 otherwise.
Clearly X i is an exponential time with mean 1/r when the ith supply perishes, and X i equals 0 when the ith supply is removed due to an arrival of demand. Y i is the duration of the ith downward jump with slope r . The instants T i are obtained from S i by inserting all pieces X j and Y j with j < i. Let So I r = + when A r increases or A r = 1, and I r = − otherwise (when A r decreases or A r = 0). Figure 1 is an illustration of the construction in (1) and (2), with r = 1. The construction, in the limit, gives an equivalent representation of the original A process. More precisely, we have the following proposition. Let τ be any stopping time of the A process, and τ r be the corresponding one for the A r process, τ r := τ + i:S i <τ ( X i + Y i ). Further let the random variable A(∞) and A r (∞), respectively, be distributed as the limiting distribution of the A and A r processes.
Proof (a) and (b) immediately follow from the construction of the A r process. (c) and (d) follow from dominated, respectively monotone convergence. To prove (e), note that A is regenerative. Let C be the first cycle in the A process, and C r be the corresponding one in the A r process. By (a) and (b), with probability 1. Then, by dominated convergence The implication of Proposition 1 is that it is sufficient to perform our probabilistic computations, for both transient and limiting analysis, on the (I r , A r ) process, then take the limit with respect to r . With this in mind we now turn our attention to the analysis of the (I r , A r ) process. We treat the PIS B case in the next section before we go into a slightly more general setting for PIS A.

PIS B, Poisson demands
In this section we apply the fluid formulation to the A process in PIS B, where both demand and supply are independent Poisson processes of rate μ B and λ B respectively. We give results for the stationary distributions of the (I r , A r ) process (Proposition 2) and the A process (Theorem 3). We derive the performance measures from the stationary distribution (Sect. 3.2). We also study a first passage time of the age process in PIS B, which is useful when we proceed to PIS A. It is known from general PDMP studies that the stationary distribution, or the transform of a first passage time of the (I r , A r ) process, is the solution to a certain boundary value problem of linear ordinary differential equations (ODE). We exploit this result and give our solutions explicitly.
For conciseness, we omit the subscripts in λ B and μ B when handling solely PIS B. Since the inter-supply times, and the inter-demand times, for PIS B are i.i.d. exponential random variables, The (I r , A r ) process evolves as follows. When I r = + (resp. I r = −), A r increases (decreases) with rate 1 (r ), unless A r = 1 (A r = 0). In the latter case A r stays flat until I r switches to the other state. When I r = + (resp. I r = −), I r stays for an exponentially distributed time with mean μ −1 (r −1 λ −1 ) in that state, then switches to the other, unless A r = 1 (A r = 0). In the latter case I r remains + (resp. −) for an exponentially distributed time with mean r −1 (λ −1 ) and then switches to the other, which in consequence takes A r away from the boundaries and the evolution again is driven according to the rules for A r ∈ (0, 1). The description above is exactly the so-called fluid model driven by a CTMC with a state space {+, −} (cf. Kulkarni 1997), with special behavior on the boundaries. It is easy to see that the generator of I r is as follows:

Stationary distribution
Let diag( v) denote the diagonal matrix of a vector v and It can be proved rigorously within the PDMP framework that the stationary distribution of the (I r , A r ) process has a density f r (i, x) at (i, x) ∈ {+, −} × (0, 1) and atoms p r (+), p r (−) at (+, 1) and (−, 0) respectively, which satisfy the following system of equations: where p r and f r (x) are row vectors as follows Intuitively, the ith equation of (5a)-(5c) are the global balance equations for the state sets Thus (5) becomes a standard boundary value problem of linear ODE. The explicit solution is given in the following proposition. (3). Let R be as in (4). The (I r , A r ) process has a unique stationary distribution, which has a density f r of the form

Proposition 2 Let Q, Q and Q be as in
and atoms p r (+), p r (−) at (+, 1) and (−, 0) respectively. The atoms are uniquely determined by where 0 is the zero row vector and 1 is the column vector with all entries being 1.
Proof From (5a) and (7) we get (8). Evaluating f r (1) using (8) and substituting in (5c) we get (9a). Since Q + Q is irreducible, the null space of Q exp(R −1 Q) + Q has exactly one dimension. Then p r is completely determined by the normalization equation (9b).
Remark 3.1 The matrix R −1 Q in the exponent does not depend on r and is diagonalizable.
The stationary distribution of the A process can be obtained by taking the limit r → ∞ (Proposition 1). Let where p is the proportion of time the system is out of stock. Alternatively it is simply the conditional distribution given that I r = + and A r < 1, or I r = − and A r = 0 (e.g. Asmussen 2003, Proposition 1.12, p. 309). Hence, for any given r > 0, where σ r is the probability that I r = + and A r < 1, or I r = − and A r = 0, Either way we give the explicit result in the following theorem and omit the proof.
Theorem 3 (Stationary distribution) The age of the oldest stock in PIS B has a unique stationary distribution, which has a density and an atom Remark 3.2 Theorem 3 is well known in the context of the finite dam model and coincides with Corollary 2.7 of Perry and Asmussen (1995).
We are interested in the distribution of the stock level in steady state, denoted by N . The generating function of N can be obtained by conditioning on the age of the oldest stock in steady state, namely E(z N ) = E(E(z N |A)), |z| ≤ 1. Notice that if A = a > 0, then N − 1, the number of supplies since the arrival of the oldest stock, has a Poisson distribution with mean λa, i.e., E(z N −1 |A = a > 0) = e −λa(1−z) . The probability distribution can be obtained from the coefficients of the power series of the generating function. We give the result in the following corollary and omit the proof.

Corollary 4 (Stock Level) Let
The moment generating function of the stock level in steady state is Remark 3.3 A direct calculation of the long-run average stock level is as follows.

Performance measures
Recall that the performance measures we consider for PIS B are long-run average stock level (s B ) and perishing rate ( B ). Also the substitution demand rate (m B ) is of interest, since the performance of PIS A depends on m B . Clearly m B = μp. Since every supply is either perished or issued to a demand, and every demand is either satisfied by a supply or lost, in the long run, we have λ − B = μ − m B . We use this conservation law to compute Corollary 5 (PIS B Performance) When λ = μ, When λ → μ,

First passage times
We start by identifying the ON and OFF periods in PIS B, which are closely related to the performance measures of our interest. The ON (resp. OFF) period is the period during which A(t) > 0 (A(t) = 0) and demands are satisfied (unsatisfied and routed to PIS A). Let U (resp. D) be the generic random variables for the duration of the ON (OFF) period. The ON and OFF periods affect the performance of PIS A by modulating its demand process.
Obviously D is exponentially distributed with mean λ −1 . Now we define, for the (I r , A r ) process, the first passage time, τ r , and the Laplace-Stieltjes transform (LST), φ α,r , which are related to U . Let Let φ(α) = E(e −αU ) be the LST of the ON period. From Proposition 1 we have An immediate result from the PDMP theory is that the column vector φ α,r (x) = [φ α,r (+, x) φ α,r (−, x)] T satisfies the following differential equation: where I is the identity matrix. Next we specify the boundary conditions. Clearly Recall that the trajectory of A r , given I r (0) = + and A r (0) = 1, stays on the boundary for an exponentially distributed time with mean r −1 , and then leaves upon I r switching from + to −. Therefore we have the following factorization, With these two boundary conditions, the solution is uniquely determined as given in the following proposition. (3). Let R be as in (4). The LST of the first passage time defined in (14) is given by

Proposition 6 Let Q, Q and Q be as in
where Proof From (17) and (16) Solving the equation above yields (20).
The following theorem is a direct result from Proposition 6 and (15).

Theorem 7 (ON Period, LST) The LST of the ON period of PIS B is given by
Remark 3.4 Theorem 7 is in agreement with Perry and Asmussen (1995, Corollary 3.1, Model II).
Remark 3.5 An alternative way to determine p, instead of using (10) is to use the mean of U and D. Notice that the process {1 {A(t)>0} } t≥0 is an alternating renewal process. Then The ith moment of the ON period can be obtained from Theorem 7 by We list the explicit formulas for the first three moments in: These moments are useful when we deal with PIS A in the next section. The squared coefficient of variation (SCV) of the ON period, defined as c 2 U (λ, μ) = E(U 2 )/E 2 (U ) − 1, is given as follows.
The observations above provide certain heuristic information to explore the parameter space in the numerical experiments, presented in Sect. 5.

PIS A, modulated Poisson demands
As mentioned in Sect. 1, the demand process of PIS A is a modulated Poisson process. For PIS A, the resulting (I r , A r ) process by the fluid formulation of (1) and (2) is not Markovian any more. We adopt an approximation as follows. First we use a phase type (PH) distribution (cf. Neuts 1975) with an irreducible representation ( γ , T ) to approximate the distribution of the ON period U in PIS B, i.e. To specify the representation ( γ , T ), a probably over-simplified solution is to take γ = 1 and T = −1/E(U ), i.e., approximate U by an exponential random variable with mean E(U ). It is worth noting that closed form solutions are developed in Osogami and Harchol-Balter (2006) for mapping a general distribution (on the positive half-line) to a PH distribution, which matches the first three moments. In this section let us assume ( γ , T ) is given.
The PH approximation enables us to enlarge the state space of the I r process in order to render the (I r , A r ) process Markovian. Let n − 1 ≥ 1 be the number of phases in the PH random variable that approximates the ON period. Let J (t) = i, i = 1, 2, . . . , n − 1 if PIS B is in phase i of an ON period and J (t) = n if PIS B is in an OFF period at time t. Then the process {J (t), t ≥ 0} is a CTMC with a generator M (of size n) as follows.
The steady-state analysis now proceeds as a straightforward extension of our treatment of PIS B. Recall that the fluid formulation approach is to construct a fluid model driven by a CTMC of finite state space, the I r process. The additional modulating CTMC J introduced here induces an extension of the state space of I r from {+, −} to {+1, +2, . . . , +n, −1, −2, . . . , −n}. This extended state now captures information about the J process as well, as the numeric part in our notation. For matrix notations, we index the states by 1 through 2n in the order they are listed above. The essence of the construction in (1) is to insert pieces (parameterized by r > 0) into the original A process so that, as r → ∞, these pieces vanish and the two processes coincide with probability 1. During the inserted periods, the J process is not allowed to make a transition, i.e., its state will be suspended. This implies that the conditional (I r , A r ) process is identical to the original (J, A) process, and hence, the same approach as (11) can be employed to obtain the stationary distribution of (J, A). Specifically we extend the matrices in (3) as follows.
where I is the identity matrix of size n. Notice the diagonal matrices r I in Q, and r λ A I in Q, effectively suspend the J process by disallowing phase transitions. The entries of the row vector μ are the demand rates of PIS A modulated by J , i.e., all entries of μ are μ A , except the last one, which is μ A + μ B . The dimension of R is also extended as For the stationary distribution, Proposition 2 is readily extendable. First we introduce the following vector notations.
Then we extend the notation introduced in (6) as Equipped with these notations, we extend Proposition 2 as follows.
Proposition 9 Let Q, Q and Q be as in (21). Let R be as in (22). The (I r , A r ) process has a unique stationary distribution, which has a density f r of the form (8) and atoms p r (+ j), p r (− j) at (+ j, 1) and (− j, 0) respectively, for j = 1, 2, . . . , n. The atoms are uniquely determined by (9).
For the stationary distribution of the (J, A) process, f (x) and p, we extend (10)-(12) as follows. (23) However we can no longer get explicit expressions of such a simple form as in Theorem 3. We can again compute the performance measures from the stationary distribution of the (J, A) process. Clearly m = μ p T . By the conservation law for supply and demand, we get (13), we get the average stock level Although first passage times in PIS A are irrelevant to the performance measures in our current consideration, we note that Proposition 6 is readily extendable as well. Denote Then we have: Proposition 10 Let Q, Q and Q be as in (21). Let R be as in (22). The LST of the first passage time defined in (14) is given by To this end, the approach outlined above provides a unified treatment for similar models where the supply and demand are two independent Markovian arrival processes. Certainly it becomes delicate to construct the Q t matrix and even more so to specify the Markovian arrival process of the unsatisfied demands. The size of the state space of the driving Markov chain may grow rapidly, which is a common limitation for approaches based on Markovianization with supplementary variables.

Numerical experiments
To validate our approach, in this section we conduct experiments which focus on approximation errors in comparison to discrete event simulation of the inventory systems. Three approximations are in consideration: • Poisson approximation (named PA): isolated PIS A has demand stream as Poisson process (PP) with rate μ A + m B . • Exponential approximation (named EA): isolated PIS A has demand stream as ON-OFF Markov-modulated PP with rates μ A and μ A + μ B respectively. The ON period is approximated by an exponential random variable, the mean of which coincides with the mean duration of the ON period in PIS B. • Three-moment approximation (named M3A): Same as EA, except that the ON period is approximated by a PH random variable that matches the first three moments of the ON period, using the algorithm developed by Osogami and Harchol-Balter (2006, Fig . 8).
Clearly, PA is the most straightforward one to use and it is usually seen in the literature (e.g. Zhao et al. 2006;Reijnen et al. 2009). The other two belong to the PH approximation discussed in Sect. 4. Compared to EA, M3A uses a more refined approximation for the ON period, thus better approximates the exact superposition demand process of PIS A. In principle we can continue to refine such an approximation to be as precise as desired, attributed to the denseness of the class of PH distributions, which is a well-known fact stated as follows (e.g. Wolff 1989, p. 271). For any non-negative random variable, the ON period U in our case, there exists a sequence of PH random variables that converges to U in distribution. The main difficulty in practice is to find the PH approximation. A fruitful approach in this area is the so-called moment matching algorithm, for which we refer the reader to Osogami and Harchol-Balter (2006) and the references therein.
Although it would be of practical interest to bound the error for each approximation (so that one can choose an approximation of the lowest refinement level from all approximations that meet a given precision requirement), we are not going to pursue the matter here. However we are interested to see whether it is possible to draw qualitative conclusions at this stage. Intuitively, the significance of an exact description of the total demand process diminishes, if the substitution demands have a relatively negligible contribution. We may use a ratio as follows to roughly quantify the impact of PIS B on PIS A, Then one may think of η as an "amplifier" of the approximation error and expect that, for any of the three approximations, the accuracy decays as η increases, given other possible factors remain the same. Heuristically, another factor of importance seems the variability of the ON period in PIS B in the sense that the larger the normalized variance, the more significant an exact description of the total demand process will be. Hence the advantage of using a more refined approximation should be more prominent when the ON period is of higher variability. When we consider these factors simultaneously, it is plausible to conjecture that the normalized error is expected to increase in η, in c 2 U , and to decrease in the refinement level of the approximations. Therefore, at least, one would be somewhat assured to use PA when both η and c 2 U are small, otherwise be alerted about the potential pitfall. The main purpose of the experiments in this section is to evaluate the three approximations. Meanwhile we try to seek numerical evidence for the above discussion on the choice of an economically adequate approximation. The following experiments are carried out in three settings to generate test cases. We start with a setting to get a general review of the approximations for a fairly wide range of system parameters. This setting also conveniently serves as a cross validation for our implementations of the analytical computation and the simulation. Then we proceed to an interesting extreme setting for which we are able give a high contrast demonstration of the approximation quality. Finally we return to our motivating application and test the approximations for several series of realistic system parameters.

The wide setting
The test cases in this setting are generated as follows. We fix μ A = 1 (so that η = m B ) and λ A = μ A . 1 Then we vary both μ B and the supply-demand ratio of PIS Thus we obtain 25 test cases in total.
The side-by-side comparison of M3A and simulation is reported in Table 1. We do not list the perishing rate ( ) since it is computed in terms of the demand lost rate (m).
We make the following observations from a close examination of Table 1. First, the wide coverage of these 25 test cases is evident by the ranges of η (from nearly 0 to 3.04), c 2 U (from 0.1 to 2.67) and the size of the matrix T (from 2 to 20). Second, the precision is extraordinarily high. The absolute differences between the values by M3A and by simulation are in the scale of 10 −4 . Hence M3A is remarkably accurate in this setting. Third, performance evaluation by M3A is efficient. For all cases, it takes merely several milliseconds on an ordinary desktop computer, which makes M3A accessible for evaluation-intensive optimization procedures.
On the other hand, EA (even PA) also performs reasonably well in this setting. Table 2 illustrates the relative errors. The outcome can be mostly explained by our intuitive rationale about the relation between η, c 2 U and the approximation error. For example, a comparison between Case λ B = 2, μ B = 4 and Case λ B = 4, μ B = 2 reveals the influence of η; a comparison between Case λ B = 0.25, μ B = 1, Case λ B = 0.5, μ B = 1 and Case λ B = 4, μ B = 4 reveals the influence of c 2 U . The top three errors of PA indeed involve large η and/or c 2 U . For the remaining 22 cases, the error of PA is less than 4 %. This observation motivates the next setting where we shall see the effort for a refined approximation is well paid off.

The extreme setting
Here we consider an interpretation of our PIS A/B model as follows. Let us think of A and B as two quality grades of a product with shelf life, say, one month. A customer who demands a grade B product will always be willing to accept a grade A (which is a higher grade) The length of the 99 % confidence intervals (CI-99) of all simulation estimates are in the scale of 10 −4 The M3A computation times are several milliseconds (λ A = μ A = 1) product but never the other way around. For certain reasons, we regard all customers equally important and decide to satisfy any demand whenever it is possible. Now suppose grade B product is a "fast mover", a product of high supply and demand rates, say, thousands of units per month; grade A product is a "slow mover", a product of low supply and demand rates, say, tens of units per month. An interesting case arises if PIS B is balanced, i.e., the supply rate equals the demand rate. In this case a shortage of grade B product can be quite unlikely. For example, if λ B = μ B = 1000, then for a long term the probability of shortage is less than 0.1 %. However, such an event is not too unlikely to be negligible for PIS A. For example, if μ A = 9, then the substitution demand amounts to almost 10 % of the total demand of PIS A. Whenever a shortage of grade B product happens, PIS A bears a billowing surge of demands. This hints at a high variability in the superposition demand process of PIS A, for which we would expect a high contrast in accuracies of the three approximations.   We generate 10 test cases as follows. We fix λ A = μ A = 1, balance supply and demand for PIS B, then vary μ B in {2 i ; i = −2, . . . , 7}.
The side-by-side comparison of M3A and simulation is reported in Table 3. Figure 3 illustrates the relative errors for all three approximations. As we expect, the errors increase as η and c 2 U increase simultaneously. We can clearly see that among the three approximations, M3A is the most accurate while PA is the least. The experiment in this setting also reveals a limitation of our approach. Here we record an error slightly higher than 5 % for the average stock level evaluated by M3A. If we keep increasing μ B , then the accuracy of M3A may eventually become insufficient. This observation suggests a direction of further investigation in heavy-tailed traffic queueing systems for alternatives.
Another observation is that our approximations apparently tend to under-estimate the average stock level and the demand lost rate.

The realistic setting
In this setting, we start from (λ A = 25, μ A = 20, μ B = 30, λ B = 40), which is supposedly close to the reality of our motivating blood bank application. We try to see the influence of substitution, as well as the relation between system parameters and performance. We generate five series of parameters as follows.
(a) Vary each parameter of the four, so that the supply-demand ratio for the corresponding system (PIS A or B) varies within the interval [0.5, 1.5]. This results in four series, named λ A -series, μ A -series, etc. For example, the λ A -series is generated by varying λ A from 10 to 30. (b) Vary the shelf life parameter (b) within the interval [1,2]. This results in one series, named b-series.
For λ A -, μ A -and b-series, we expect hardly any accuracy difference between the three approximations. All of them should be quite accurate, since the influence of PIS B on PIS A is negligible in all cases (the largest η being in the scale of 10 −5 ). Our expectation is indeed

Summary and further research
In this study we introduce a prototype PIS with random input. As far as is known to the authors, such inventory models that are subject to one way substitution have never been introduced in the operations research literature. Accordingly, due to the intricateness of the model we have focused mainly on the first part of the problem, namely, on the performance analysis part of the stochastic model. It appears that even if the demand processes of type A and of type B are independent Poisson processes the total demand arrival process into PIS A is neither a Poisson process nor a renewal process. We suggest a methodology based on a certain approximation for the analysis of PIS A. The Laplace transform of the ON period can be found. Thus, the moments of the ON periods can be obtained by taking derivatives at 0. Our approach is based on the idea that the law of the ON period can be approximated by phase-type (PH) distributions that have the same moments as those of the original distribution of the ON period. Intuitively, as the fitness between the moments increases, the approximation improves. The drawback of our approach is the fact that, theoretically, moments do not determine the distribution. For that reason, we are unable to compute bounds to that approximation. However, practically and intuitively, the higher fitness among the moments the better approximation is obtained. The model can be extended in several directions, which are left for further research: • In the current study we focus only on the marginal analysis of PIS A. The joint analysis of PIS A and PIS B is completely forsaken. • In practice, there are more than two types of blood. As a result, more complicated relations than a one way substitution exist. • Both the item and the demand arrival processes are assumed to be Poisson processes. This assumption can be extended to arrival processes in a random environment. Namely, it is known that demands for blood portions are subject to changes due to disasters such as tsunami, earthquake, terror attack, and so on. • The demand arrival rate can be controlled by increasing the publicity according to the state of the content level. • Considering optimization, where the decision variables are the arrival rate, the demand rate and the shelf life of the items.