1.INTRODUCTION
This paper concerns the evaluation and analysis of customer service level in the situation that a customer order typically consists of several items and customers are satisfied only if their requests are met immediately. As a measure of customer satisfaction, service levels have become increasingly important for companies under competitive pressures in today’s marketplace in order to provide quicker response to customer needs. One of the most commonly used service level measures in such situation is the order fill rate, probability of filling an entire customer order immediately from the shelf.
Song (1998) called the order fill rate as the orderbased performance measure. Contrast to the orderbased approach, however, most standard inventory models do not take into account connections between items. They assume that demands for each item are independent of the others. Thus Song (1998) referred to this as the itembased approach. In the orderbased approach, since each customer order requires the simultaneous availability of several items, the order fill rate askes the following question. For any given safety stock level of each item, what is the probability that a customer order can be satisfied immediately?
The issue arose when Park and Seo (2013) analyzed the inventory operations practice of a spare parts distributor for ship engines and generators. After analyzing the historical purchaseorderrequest (POR) data, they recognized that a considerable proportion of previous orders had been cancelled due to the lack of a few spare parts, which were out of stock. In other words, if an order included at least one spare part that was out of stock, the order itself was cancelled, even though the other spare parts had large quantities of stocks. The historical POR data showed that the portion of orders cancelled due to stockedout parts was over 30%. In many orders, shipping companies buy spare parts only when a spare parts distributor keeps all of the spare parts in stock, because the shipping companies want to purchase all of the spare parts together. Under this purchase pattern, if one of the spare parts is not in stock, even though all of the other spare parts are in stock, the situation is the same as if all spare parts were outofstock. Park and Seo (2013) referred to this type of dependence as ‘purchase dependence.’
Purchase dependence is different to demand dependence. Purchase dependence deals with the purchase behavior of customers, whereas demand dependence deals with demand correlation between items, between regions, or over time. Purchase dependence can be observed in several areas such as marketing, manufacturing systems, and distribution systems. With the perspective of marketing, customer purchase patterns have been discovered by extracting associations among items from purchase transactions, so called market basket analysis. It discovers rules showing that customers are more or less likely to buy certain items given the purchase of other items. Many applications of market basket analysis are available in the literature for the obvious reason that information obtained from the analysis can be used to help store managers make better decisions about shelf allocation, product promotion, and so on. See Athappilly et al. (2010), Yang et al. (2013), and the references therein.
In a multicomponent productioninventory system, such as an assembletoorder (ATO) system, purchase dependence for different components may be induced by a billofmaterials or product options. There are multiple components and multiple products. Inventories are kept only at the component level; final products are assembled in response to customer orders. A customer order for an end product usually requires a number of different components. The order is satisfied only if all the required components are available; otherwise, it is lost or becomes a backorder. Multiitem distribution systems in which customer orders consist of a set of different items, such as online retailing and mailorder businesses, have a similar issue.
This paper proposes a new approximate approach that calculates the order fill rate under purchase dependence. The idea of this paper’s approach was motivated by the decomposition approach proposed by Iravani et al. (2004). Like Iravani et al. (2004), we divide customer orders into item demands, that is, construct new item demand processes that have the same marginal distributions as the original order processes. However, we calculate fill rates of all order types to approximate the order fill rate, whereas Iravani et al. (2004) use the component fill rate of each individual item to approximate the order fill rate.
In addition, we derive the fill rates of all order types from the stationary joint distribution of onhand inventories modelled by a quasibirthanddeath (QBD) process. The stationary joint distribution is the solution of the balance equations plus a normalization equation. However, conventional computational methods for solving the linear equations suffer from the curse of dimensionality for the large number of equations and unknowns. Even though Song et al. (1999), Iravani et al. (2003) and Gao et al. (2010) utilize matrixgeometric techniques, there still exists the curse of dimensionality. Thus we attempt to use an iterative method, such as the GaussSeidel method.
The GaussSeidel method is a wellknown method used to solve a linear system of equations: Ax = b, A = (a_{ij}) ∈ R^{m×m}, x, b ∈ R^{m}. It is a modification of the Jacobi method, so requires less iteration to produce the same degree of accuracy. In order to improve the convergence rate of the GaussSeidel method, several preconditioned iterative methods have been reported in the literature. Table 1 summarizes these studies. Unfortunately, however, the GaussSeidel method and several preconditioned iterative methods cannot use to solve our linear equations system because these methods diverge. Thus, in order to avoid the curse of dimensionality and prevent the solution from diverging, we develop the greedy iterative search algorithm (GISA) based on the GaussSeidel method.
The remainder of this paper is organized as follows. In the following section, we review the related literature. In Section 3, we explain the model formulation that calculates the order fill rate under purchase dependence through the stationary joint distribution of onhand inventories modelled by a QBD process. Section 4 describes the procedure of GISA that calculates the stationary joint distribution. In Section 5, we present the computational analysis conducted to verify the accuracy of GISA, investigate the impacts of an initial solution, and obtain the managerial insights. In Section 6, we present our conclusions.
2.LITERATURE REVIEW
In the ATO system, considerable research has been devoted to the performance evaluation: for an overview of modelling and analyzing the ATO system, see Song and Zipkin (2003) and Benjaafar and EI Hafsi (2006). We only focus on the studies on the performance evaluation (e.g., the order fill rate) of ATO systems. Most of the existing research on the performance evaluation focused on the base stock policy and the FirstComeFirstServe (FCFS) service discipline, but they differed in their model assumptions and solution approaches. These studies can be divided into two groups: the exact modelling approach and the approximate modelling approach. Table 2 provides a summary of the studies that focused on the performance evaluation.
As the exact modelling approach, Song (1998) took a first step in investigating the orderbased approach. She presented an exact evaluation of the order fill rate for a multiitem, multiproduct ATO system with backlogging and constant lead times. A customer order consisted of several items in different amounts. The computation of the order fill rate was accomplished through a series of convolution of onedimensional distributions. On the other hand, the other studies for the exact modelling approach formulated the system using a continuoustime Markov chain and derived performance measures from the exact stationary joint distribution of outstanding item orders which can be determined by modelling the system as a QBD process. An exact performance analysis was carried out using matrixgeometric techniques.
Song et al. (1999) studied a multicomponent, multiproduct ATO system, where each of the components is made to stock and replenished by its dedicated machine that has exponential processing times. They designed a generalized model that has both complete backlogging and lost sales as a special case. In addition, they distinguished total order service (TOS), which means that an order is fulfilled completely or rejected as a whole and partial order service (POS), which means that partial fulfilment occurs. Iravani et al. (2003) and Gao et al. (2010) analyzed a similar situation but took a somewhat different perspective. Iravani et al. (2003) incorporated customer flexibility in ATO systems. In their model customers might be willing to change some of their requested components that are not in stock (i.e., flexible customers). Meanwhile, in the model of Gao et al. (2010) each item was replenished by an independent unreliable machine.
As mentioned by Iravani et al. (2004), one major limitation with the exact modelling approach for ATO systems is its inability to analyze large ATO systems. Current exact modelling approaches are computationally intensive mainly due to the need for incorporating correlations among demands for different items. The half part of Table 2 shows studies based on the approximate modelling approach.
Vliegen and van Houtum (2009) and Hoen et al. (2011) assumed constant lead times. In the study on service tools by Vliegen and van Houtum (2009), both demands and returns were coupled. The approximation of Vliegen and van Houtum (2009) defined two adjusted systems with minimum and maximum coupled returns and combined them with the coupling factor. However, Hoen et al. (2011) combined adjusted systems with minimum and maximum coupled demands for ATO systems.
Dayanik et al. (2003) considered an ATO system with backlogging, lostsales, POS, exponential processing times and single servers for the production of components, whereas Iravani et al. (2004) employed a batch replenishment policy. Dayanik et al. (2003) formulated and examined a number of performance bounds for an ATO system with Markovian distributions. Iravani et al. (2004) introduced a decomposition method that estimates the performance measures of each individual queue by decomposing the large system of dependent parallel queues into individual independent queues and then uses the performance measures of each individual queue to approximate the performance measures of the system by utilizing the concept of chaining conditional probabilities.
Song and Yao (2002) considered only one endproduct. It simplified the analysis considerably and allowed them to obtain an easier evaluation of the order fill rate. They derived easytocompute performance bounds and used them as surrogates for the performance measures in several optimization problems that seek the best tradeoff between inventory and customer service. Greedy type algorithms were developed to solve the surrogate problems. Lu et al. (2003) extended the singleproduct ATO system studied by Song and Yao (2002). Based on the joint distribution and the first two moments, they derived approximations and bounds for the order fulfillment performance measures. Zhou and Chao (2012) developed an extremely simple SteinChen approximation as well as its errorbound of the order fill rate for an ATO system. Their approximation gave an expression for the order fill rate in terms of itembased fill rates.
3.MODEL FORMULATION
3.1.Model Description
Under the situation that purchase dependence exits, we consider a retailer system that sells J different items with Ω = {1, 2, …, J} being the set of all item indices. For any K ⊆ Ω, let K denote the number of elements in K. The overall order process is stationary in time and forms a Poisson process. Each customer requests at most one unit of each item but may require several items simultaneously. For any K ⊆ Ω, we say that an order is of type K if it requires one and only one unit of each item in K and 0 units in Ω  K. Let T = {1, 2, …, T} denote the set of order types, where 1 ≤ T ≤ 2^{J}1. We assume that there is a fixed probability that an order is of type K. Each order type is independent of the other order types. This implies that the demand process for each item is also a Poisson process. This paper follows the basic model assumptions and notations made in Song et al. (1999).
Orders are filled on a FCFS basis. When an order arrives and a retailer has some, but not all, of its items in stock, the order is lost entirely due to purchase dependence. This situation is known as TOS in Song et al. (1999). Each item is controlled by an independent base stock policy. Let s_{i} be the base stock level for item i. That is, at each demand epoch, if the inventory position (i.e., the inventory on hand plus inventory on order) of item i is less than s_{i}, then order up to s_{i}. Without loss of generality, we assume that at time 0 item i is stocked at level s_{i} for all i. Then each demand for item i triggers an order for that item if the item i is delivered to the customer. Hence there can be, at most, s_{i} outstanding orders of item i. Orders for item i are replenished by a facility i and the lead times are exponentially distributed with the rate μ_{i}. The facility processes the orders on a FCFS basis.
3.2.The Order Fill Rate
This section explains the new approximate approach devised by this paper to calculate the order fill rate under purchase dependence. As mentioned in Section 1, the approximate approach utilizes the idea of decomposition approach. The approximate approach is supported by Song (2002). She presented that a dominating factor for the computation time of order fill rate is not how many items of which an order consists, but rather the number of order types that have overlapping items with this order. Thus the easiest case is a pure system, where there is only one order type. In this approximate approach, we decompose customer orders into item demands, that is, construct new item demand processes that have the same marginal distributions as the original order processes. Then we calculate the fill rate of a pure system that has only one order type K. (We refer to the system as the typeK pure system.) The order fill rate of a retailer system is approximated by combining all fill rates of typeK pure systems.
For any i ∈ Ω and K ∈ T, let λ be the overall order rate, q^{K} be the probability that an order is of type K, q_{i} be the probability that an order requires item i (i.e., ${q}_{i}={\displaystyle \sum _{K:i\in K}{q}^{K}}$), and λ_{i} be the aggregate demand rate of item i (i.e., λ_{i} = q_{i}λ). Then we approximate the overall order rate for the typeK pure system by averaging the demand rates of items in the typeK order. Let
${\lambda}^{pK}$ = the overall order rate for the typeK pure system
For example, consider the case 1: J = 3 with T = 7 in Table 5. In the case 1, λ^{p}^{1} = λ_{1}, λ^{p}^{2} = λ_{2}, λ^{p}^{3} = λ_{3}, λ^{p}^{4} = (λ_{1} + λ_{2})/2, λ^{p}^{5} = (λ_{1} + λ_{3})/2, λ^{p}^{6} = ( λ_{2} + λ_{3})/2, and λp7 = ( λ_{1} + λ_{2} + λ_{3})/3.
By denoting that

I_{i} = the inventory on hand of item i in steadystate, 0 ≤ I_{i} ≤ s_{i}, the fill rate of typeK pure system is defined by

F^{pK} = the joint probability that all items in the typeK pure system are filled immediately
Then the order fill rate of a retailer system F is calculated approximately by
The order fill rate, F expressed in Eq. (3) is a weighted average of fill rates of typeK pure systems. For example, the case 1 in Table 5 has 7 pure systems. Then the order fill rate is $F={q}^{1}{F}^{p1}+{q}^{2}{F}^{p2}+{q}^{3}{F}^{p3}+{q}^{4}{F}^{p4}+{q}^{5}{F}^{p5}+{q}^{6}{F}^{p6}+{q}^{7}{F}^{p7}$.
3.3.The Stationary Joint Distribution
In this section, we consider a pure system that has only one order type K. Due to purchase dependence, a typeK order is entirely lost if, upon its arrival, at least one item i (i ∈ K) is out of stock. Since the fill rate of typeK pure system is defined by the joint probability that all items in the typeK pure system are filled immediately, the fill rate of typeK pure system is determined by the joint distribution of (I_{i} : i ∈ K). This section illustrates the procedure of deriving the joint distribution of (I_{i} : i ∈ K).
It is evident that a stochastic process {(I_{i}(t) : i ∈ K), t ≥ 0}, where I_{i}(t) is the inventory on hand of item i at time t, is a continuoustime Markov chain with finite state space $\left\{\text{n}=\left({n}_{i}:\text{\hspace{0.17em}}i\in K\right)0\le {n}_{i}\le {s}_{i},\text{\hspace{0.17em}}i\in K\right\}$. There are $\prod _{i\in K}\left({s}_{i}+1\right)$ states. A state transition can only occur if a demand or an item on order arrives. That is, with transition rate λ^{pK}, the state n transits to the state n′ = (n ′ _{i} : i ∈ K), where
and with transition rate μ_{i} (i ∈ K), the state n transits to the state n″ = (n ″ _{i} : i ∈ K), where
It is easy to see that this Markov chain is irreducible, so its stationary distribution exists uniquely (Song et al., 1999). Using Eqs. (4) and (5), we can establish the balance equation for each state. For example, Figure 1 shows the transition rate diagram and balance equations for the type3 pure system in which the type3 order requires both items 1 and 2. In the equilibrium situation the rate of leaving any node must be equal to the rate of entering that node. The ‘rate in = rate out’ principle allows for the following generalization.
Let I be a vector of (I_{i} : i ∈ K) and p_{I} be the equilibrium probabilities. Then
where δ_{i} = 1 for I_{i} < s_{i} or 0 otherwise, ξ = 1 for ∀(I_{i} > 0) or 0 otherwise, e_{i} is a unit vector having 1 at the position of I_{i}, 1 is an allones vector, and p_{I} = 0 for any I_{i} ∉ [0, s_{i}].
The stationary joint distribution is then the solution of balance equations plus a normalization equation. Conventional computational algorithms, such as Gaussian elimination, for solving the linear equations can be employed. However, the direct approach can suffer from severe computation and space complexities for the large number of equations and unknowns. In order to reduce the computation and space complexities, Song et al. (1999), Iravani et al. (2003) and Gao et al. (2010) utilized the special structure of the QBD process that allows obtaining a unified, matrixgeometric solution of the stationary joint distribution. (Ye and Li, 1994) presented the folding algorithm as a new computational method for steady state analysis of finite QBD process) Still, there exists the curse of dimensionality. Thus we develop the GISA such that it avoids the curse of dimensionality using the special structure of the QBD process depicted by Eq. (6) and prevents the solution from diverging by adding the newly invented steps into the GaussSeidel method.
4.GREEDY ITERATIVE SEARCH ALGORITHM
The balance equations generated by Eq. (6) and the normalization equation form a linear equations system: Ax = b, A = (a_{ij}) ∈ R^{m×m}, x, b ∈ R^{m}. A is the coefficient matrix, b is the righthand side vector, x is the vector of unknowns (i.e., p_{I}, the equilibrium probabilities), and the order $m={\displaystyle \prod _{i\in K}\left({s}_{i}+1\right)}$. The linear equations system can be large and sparse. The large number of equations and unknowns bring a serious storage concern to a computer computation. When it is possible, Gaussian elimination remains a very economical, accurate, and useful algorithm. However, when the order m is so large that it is impossible to store the fillin resulted from the Gaussian elimination, it is desirable to solve the linear equations system by methods that never alter the matrix A and never require storing more than a few vectors of length m. Iterative methods are especially suitable for this purpose.
In this section, we describe the procedure of GISA. In the GISA, beginning with an initial solution vector x^{(0)}, we generate a sequence of solution vectors x^{(1)} → x^{(2)} → … according to the iteration scheme. We hope that as k → ∞, x^{(k)} will converge to the exact solution. Figure 2 shows the pseudocode of GISA. In the initialization step, the GISA sets up an initial solution x^{(0)} and calculates the maximum absolute residual r_{max} corresponding to x^{(0)}. Let [y]_{i} be the ith component of vector y.
Then the GISA generates a sequence of solution vectors x^{(k)} according to the following scheme. Let ${x}_{i}^{\left(k\right)}$ denote the ith component of the k th iterate solution vector x^{(k)} (i.e., [x^{(k)}]i) and b_{i} the ith component of the righthand side b (i.e., [b]_{i}). The GISA corrects the ith component of the current solution, in the order i = 1, 2, …, m, to annihilate the ith component of the residual. The solution is updated immediately if the newly computed component (i.e., probability) ${x}_{i}^{\left(k\right)}$ is nonnegative and the newly computed maximum absolute residual does not increase. These steps are invented in this paper to prevent the solution from diverging. At the end of iteration k, the GISA calculates the summation of absolute residuals r_{sum}.
The GISA continues to generate a sequence of solution vectors x^{(k)} until a convergence is reached. In the GISA, we regard the following conditions as a convergence: (i) the iteration k reaches to the maximum number of iterations; (ii) the maximum absolute residual r_{max} is less than the maximum value of tolerance; or (iii) there is no change in the solution vector x^{(k)}. As the order m becomes so large, the sequence of solution vectors x^{(k)} would be influenced by an initial solution x^{(0)}. The impacts of initial solution are examined in Section 5.2.
5.COMPUTATIONAL ANALYSIS
This section describes numerical and simulation experiments conducted to verify the accuracy of GISA, investigate the impacts of an initial solution, and obtain the managerial insights through the analysis of solution procedure and the order fill rates.
5.1.The Accuracy of GISA
In this paper, we developed the GISA to derive the stationary joint distribution of onhand inventories in the typeK pure system, so eventually the fill rate of the pure system. The fill rate of typeK pure system is defined as the joint probability that all items in the typeK pure system are filled immediately. We measure the accuracy of GISA with the fill rate of typeK pure system.
In order to test the accuracy of GISA, we designed the test bed as follows. We considered two pure systems, i.e., K = 3 or 5. Within both pure systems, equal base stock levels were chosen, i.e., s_{1} = s_{2} = ⋅⋅⋅ = s_{K}. For these base stock levels, we chose three values of 5, 10, and 15. By letting all replenishment rates μ_{i} (i ∈ K) be equal to 1.0, we chose three values of overall order rate for the typeK pure system so that the first value leads to a low fill rate of typeK pure system, the second value to a medium fill rate, and the last value to a high fill rate of approximately 95%.
For comparison, we calculated the exact fill rate using the Gaussian elimination. However, for the large instances which require manipulation of large matrices, the Gaussian elimination was not practical. Thus we developed a simulation program to obtain the fill rate for the large instances. For comparison purpose, we conducted the simulation run for five times per each instance. For each time, the simulation was run for 10,000 customer orders with 1,000 orders of warmup. Then the fill rate was calculated by averaging the five fill rates resulted from simulations.
The results of this comparison are presented in Table 3. For small instances, the fill rates of typeK pure system calculated by the GISA are almost same as those by the Gaussian elimination. This demonstrates the accuracy of GISA. We can verify the fact by reviewing that r_{max} and r_{sum} are near to zeros. For large instances, even though r_{max} and r_{sum} illustrate the accuracy of GISA, Table 3 shows the comparison between the GISA and computer simulation at the last column. The average value of absolute difference is 0.0028. Based on the results in Table 3, we can conclude that the GISA is a dependable algorithm to derive the stationary joint distribution of onhand inventories in the typeK pure system.
5.2.The Effects of an Initial Solution
As mentioned in Section 4, the sequence of solution vectors generated by the GISA would be influenced by an initial solution. This section investigates what impacts an initial solution gives to the sequence of solution vectors.
We differentiated the initial solutions by controlling the number of components having a positive value in the initial solution vector. In other words, we assigned a positive value to the only controlled components in the initial solution vector but zeros to the rest. In this section, we assigned an equal positive value obtained by dividing 1.0 by the number of components having a positive value.
In order to examine the impact of initial solutions on the convergent speed of the GISA, we calculated the fill rates of typeK pure system (i.e., K = 5) as increasing the number of components having a positive value in the initial solution vector (e.g., 10, 100, and 1,000). Figure 3 shows the convergent speeds for the two cases of s_{i} = 5 and s_{i} = 10. (We set that λ^{pk} = 1.0 and μ_{i} = 1.0 (i ∈ K)) From Figure 3, we can observe that, in general, the convergent speed becomes faster as the number of components having a positive value decreases. In addition, we can observe that it takes more iteration to reach to the convergence as the base stock level increases.
On the other hand, we can notice that the sequence of solution vectors will converge to a certain solution because the GISA has the steps that prevent the solution from diverging. However, we recognized that for the large instances, the GISA could converge to somewhat different solutions as shown in Table 4, depending on the initial solutions. Table 4 presents the fill rates of typeK pure system (i.e., K = 5) calculated by the GISA with five different initial solution vectors (We set that s_{i} = 10, λ^{pk} = 1.0 and μ_{i} = 1.0 (i ∈ K)). Among the five trials, it will be a wise choice to select the result of trial 4 because the result of trial 4 has the highest degree of accuracy, i.e., a minimum value of r_{sum}.
5.3.Analysis on the Order Fill Rate
This paper obtains the order fill rate of a retailer system under purchase dependence by combining all fill rates of typeK pure systems. In this section, we analyze the solution procedure and the order fill rate, thus deduce the managerial insights.
Like Section 5.1, we designed the test bed as follows. For the combination of the number of items and the number of order types, we considered four cases such that (i) case 1: J = 3 with T = 7, (ii) case 2: J = 3 with T = 7 (asymmetric), (iii) case 3: J = 3 with T = 4, and (iv) case 4: J = 5 with T = 6. In cases 1 and 2, all possible order types occur and thus T = 7. For cases 3 and 4, we take order types for individual items and one order type for the set of all items, thus T = J + 1. We mainly kept order types symmetric for the order rate, but for case 2, we considered asymmetric order types. Within all retailer systems, equal base stock levels were chosen, i.e., s_{1} = s_{2} = … = s_{J}. For these base stock levels, we chose a value of 5. By letting all replenishment rates μ_{i} be equal to 1.0, we chose three values of overall order rate for the retailer system, i.e., λ = 1.0, 1.25, and 1.50.
We defined the degree of purchase dependence (dp) as(7)
This means that dp = 1 if all orders come only from the one order type for the set of all items; and dp = 0 if all orders come from order types for individual items only. The degree of purchase dependence measures the extent of how much item demands are conjoint by customer orders. Table 5 and Table 6 shows the values of dp and associated q^{K} for cases 1&2 and cases 3&4, respectively
For each case, we calculated the fill rates of typeK pure systems using the GISA after decomposing customer orders into item demands. Then the order fill rate of the retailer system was calculated by combining all fill rates of typeK pure systems. Like Section 5.1, for comparison purpose, we conducted the simulation run for five times per each instance. For each time, the simulation was run for 10,000 customer orders with 1,000 orders of warmup. Then the order fill rate was calculated by averaging the five order fill rates. At this time, the customer orders were composed of all order types in the set T with the proportion of q^{K} for each case.
Tables 7 through 10 illustrates the fill rates of typeK pure systems and the order fill rate of the retailer system obtained by the GISA and the simulation for each case, respectively. From these results, we can observe several managerial insights. First, we can observe that the order fill rate obtained by the simulation lies between the fill rates of pure systems that consist of all items and only one item. In order words, we can conclude that these two fill rates present the lower and upper bounds for the order fill rate. The lower bound can be obtained by the fill rate of pure system that consists of all items, while the upper bound by the oneitem pure system.Table 8Table 9Table 10
Secondly, as the degree of purchase dependence decreases but other conditions remain same, it is noticed that (i) the gap between the lower and upper bounds diminishes, (ii) the order fill rate increases, and (iii) the order fill rate gets closer to the upper bound. This observation implies that when purchase dependence exists, treating item demands as independent may result in either unnecessary overstocking of some items or unsatisfactory service at the order level. Thus we can conclude that it is important to consider knowledge of purchase dependence in designing an inventory replenishment policy for better profitability and lower inventory costs. However, it re quires the further indepth future studies to derive a conclusion that when it is tolerable of treating item demands as independent in terms of the degree of purchase dependence.
Lastly, we can see the comparison of order fill rates obtained by Eq. (3) and the simulation at the last column. The average value of absolute difference is 0.017. Based on this observation, we can acknowledge that for the large instances, the order fill rates calculated by Eq. (3) can be used as an approximate alternative. However, it is needed to improve the accuracy of order fill rate because Eq. (3) simply approximates the order fill rate by weightaveraging the fill rates of typeK pure systems with q^{K} probabilities that an order is of type K. Thus it is required that the future studies devise an elaborate method of combining the fill rates of typeK pure systems to derive the order fill rate.
6.CONCLUSION
For the retailer system in which purchase patterns are characterized by purchase dependence, this paper proposed the new approximate approach to calculate the order fill rate, probability of filling an entire customer order immediately from the shelf. The new approximate approach divides customer orders into item orders and calculates fill rates of all order types to approximate the order fill rate. On the process of deriving the stationary joint distribution that is used to calculate fill rates of all order types, conventional computational methods suffer from the curse of dimensionality for the large instances.
In order to avoid the curse of dimensionality and prevent the solution from diverging, we developed the GISA based on the GaussSeidel method. Through the computational analysis that compared the GISA with the simulation, we demonstrated that the GISA is a dependable algorithm to derive the stationary joint distribution of onhand inventories in the typeK pure system. In addition, this paper obtained some managerial insights as follows. (1) The lower bound of order fill rate can be obtained by the fill rate of pure system that consists of all items, while the upper bound by the oneitem pure system. (2) As the degree of purchase dependence decreases but other conditions remain same, it is noticed that (i) the gap between the lower and upper bounds diminishes, (ii) the order fill rate increases, and (iii) the order fill rate gets closer to the upper bound. (3) Based on the comparison of order fill rates obtained by Eq. (3) and the simulation, we can acknowledge that it would be recommendable for the large instances to use the order fill rates calculated by Eq. (3) as an approximate alternative.
As the future studies, we leave two themes. First, we observed that the order fill rate gets closer to the upper bound generated by the fill rate of pure system as the degree of purchase dependence decreases. This requires the further indepth analysis on when it is tolerable of treating item demands as independent in terms of the degree of purchase dependence. Secondly, the order fill rate calculated by Eq. (3) was simply approximated by weightaveraging the fill rates of typeK pure systems with probabilities that an order is of type K. In order to improve the accuracy of order fill rate, it is required that the future studies devise an elaborate method of combining the fill rates of typeK pure systems to derive the order fill rate.