• Editorial Board +
• For Contributors +
• Journal Search +
Journal Search Engine
ISSN : 1598-7248 (Print)
ISSN : 2234-6473 (Online)
Industrial Engineering & Management Systems Vol.16 No.3 pp.363-374
DOI : https://doi.org/10.7232/iems.2017.16.3.363

# Iterative Approach to Calculating the Order Fill Rate Under Purchase Dependence

Changkyu Park*
College of Business, University of Ulsan, Republic of Korea
Corresponding Author, ckparkuou@ulsan.ac.kr
February 14, 2017 May 22, 2017 June 12, 2017

## ABSTRACT

For the business environment under purchase dependence that characterizes customer purchase patterns observed in such areas as marketing, manufacturing systems, and distribution systems, this paper proposes a 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. In order to avoid the curse of dimensionality and prevent the solution from diverging for the large instances, we develop the greed iterative search algorithm (GISA) based on the Gauss-Seidel method. Through the computational analysis that compares the GISA with the simulation, we demonstrate that the GISA is a dependable algorithm to derive the stationary joint distribution of on-hand inventories in the type-K pure system. In addition, we present some managerial insights.

## 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 order-based 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 item-based approach. In the order-based 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 purchase-order-request (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 stocked-out 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 out-of-stock. 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 multi-component production-inventory system, such as an assemble-to-order (ATO) system, purchase dependence for different components may be induced by a bill-of-materials 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. Multi-item distribution systems in which customer orders consist of a set of different items, such as online retailing and mail-order 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 on-hand inventories modelled by a quasi-birth-and-death (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 matrix-geometric techniques, there still exists the curse of dimensionality. Thus we attempt to use an iterative method, such as the Gauss-Seidel method.

The Gauss-Seidel method is a well-known method used to solve a linear system of equations: Ax = b, A = (aij) ∈ Rm×m, x, b ∈ Rm. 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 Gauss-Seidel method, several preconditioned iterative methods have been reported in the literature. Table 1 summarizes these studies. Unfortunately, however, the Gauss-Seidel 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 Gauss-Seidel 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 on-hand 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 First-Come-First-Serve (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 order-based approach. She presented an exact evaluation of the order fill rate for a multi-item, multi-product 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 one-dimensional distributions. On the other hand, the other studies for the exact modelling approach formulated the system using a continuous-time 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 matrix-geometric techniques.

Song et al. (1999) studied a multi-component, 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, lost-sales, 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 easy-to-compute performance bounds and used them as surrogates for the performance measures in several optimization problems that seek the best trade-off between inventory and customer service. Greedy- type algorithms were developed to solve the surrogate problems. Lu et al. (2003) extended the single-product 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 Stein-Chen approximation as well as its error-bound of the order fill rate for an ATO system. Their approximation gave an expression for the order fill rate in terms of item-based 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| ≤ 2J-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 si 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 si, then order up to si. Without loss of generality, we assume that at time 0 item i is stocked at level si 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, si 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 type-K pure system.) The order fill rate of a retailer system is approximated by combining all fill rates of type-K pure systems.

For any iΩ and KT, let λ be the overall order rate, qK be the probability that an order is of type K, qi be the probability that an order requires item i (i.e., $q i = ∑ K : i ∈ K q K$), and λi be the aggregate demand rate of item i (i.e., λi = qiλ). Then we approximate the overall order rate for the type-K pure system by averaging the demand rates of items in the type-K order. Let

$λ p K$ = the overall order rate for the type-K pure system

$= 1 | K | ∑ i ∈ K λ i$
(1)

For example, consider the case 1: J = 3 with |T| = 7 in Table 5. In the case 1, λp1 = λ1, λp2 = λ2, λp3 = λ3, λp4 = (λ1 + λ2)/2, λp5 = (λ1 + λ3)/2, λp6 = ( λ2 + λ3)/2, and λp7 = ( λ1 + λ2 + λ3)/3.

By denoting that

• Ii = the inventory on hand of item i in steady-state, 0 ≤ Iisi, the fill rate of type-K pure system is defined by

• FpK = the joint probability that all items in the type-K pure system are filled immediately

$= P ( I i > 0 : i ∈ K )$
(2)

Then the order fill rate of a retailer system F is calculated approximately by

$F = ∑ K ∈ T q K F p K$
(3)

The order fill rate, F expressed in Eq. (3) is a weighted average of fill rates of type-K pure systems. For example, the case 1 in Table 5 has 7 pure systems. Then the order fill rate is $F = q 1 F p 1 + q 2 F p 2 + q 3 F p 3 + q 4 F p 4 + q 5 F p 5 + q 6 F p 6 + q 7 F p 7$.

### 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 type-K order is entirely lost if, upon its arrival, at least one item i (iK) is out of stock. Since the fill rate of type-K pure system is defined by the joint probability that all items in the type-K pure system are filled immediately, the fill rate of type-K pure system is determined by the joint distribution of (Ii : iK). This section illustrates the procedure of deriving the joint distribution of (Ii : iK).

It is evident that a stochastic process {(Ii(t) : iK), t ≥ 0}, where Ii(t) is the inventory on hand of item i at time t, is a continuous-time Markov chain with finite state space ${ n = ( n i : i ∈ K ) | 0 ≤ n i ≤ s i , i ∈ K }$. There are $∏ i ∈ K ( s i + 1 )$ 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 : iK), where

$n j ′ = { n j − 1 if j ∈ K and 0 < n i for all i ∈ K n j otherwise$
(4)

and with transition rate μi (iK), the state n transits to the state n″ = (n ″ i : iK), where

$n j ′ ​ ′ = { n j + 1 if j = i and n i < S i n j otherwise$
(5)

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 type-3 pure system in which the type-3 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 (Ii : iK) and pI be the equilibrium probabilities. Then

$( ∑ i ∈ K δ i μ i + ξ λ p K ) p I = ∑ i ∈ K μ i p I − e i + λ p K p I + 1$
(6)

where δi = 1 for Ii < si or 0 otherwise, ξ = 1 for ∀(Ii > 0) or 0 otherwise, ei is a unit vector having 1 at the position of Ii, 1 is an all-ones vector, and pI = 0 for any Ii ∉ [0, si].

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, matrix-geometric 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 Gauss-Seidel 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 = (aij) ∈ Rm×m, x, b ∈ Rm. A is the coefficient matrix, b is the right-hand side vector, x is the vector of unknowns (i.e., pI, the equilibrium probabilities), and the order $m = ∏ i ∈ K ( s i + 1 )$. 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 fill-in 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 pseudo-code of GISA. In the initialization step, the GISA sets up an initial solution x(0) and calculates the maximum absolute residual rmax 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 ( k )$ denote the ith component of the k th iterate solution vector x(k) (i.e., [x(k)]i) and bi the ith component of the right-hand 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 ( k )$ 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 rsum.

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 rmax 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 on-hand inventories in the type-K pure system, so eventually the fill rate of the pure system. The fill rate of type-K pure system is defined as the joint probability that all items in the type-K pure system are filled immediately. We measure the accuracy of GISA with the fill rate of type-K 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., s1 = s2 = ⋅⋅⋅ = s|K|. For these base stock levels, we chose three values of 5, 10, and 15. By letting all replenishment rates μi (iK) be equal to 1.0, we chose three values of overall order rate for the type-K pure system so that the first value leads to a low fill rate of type-K 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 warm-up. 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 type-K 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 rmax and rsum are near to zeros. For large instances, even though rmax and rsum 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 on-hand inventories in the type-K 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 type-K 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 si = 5 and si = 10. (We set that λpk = 1.0 and μi = 1.0 (iK)) 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 type-K pure system (i.e., |K| = 5) calculated by the GISA with five different initial solution vectors (We set that si = 10, λpk = 1.0 and μi = 1.0 (iK)). 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 rsum.

### 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 type-K 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., s1 = s2 = … = sJ. 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)

$d p = ∑ K ∈ T | K | − 1 J − 1 q K$
(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 qK for cases 1&2 and cases 3&4, respectively

For each case, we calculated the fill rates of type-K 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 type-K 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 warm-up. 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 qK for each case.

Tables 7 through 10 illustrates the fill rates of type-K 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 one-item 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 in-depth 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 type-K pure systems with qK 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 type-K 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 Gauss-Seidel 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 on-hand inventories in the type-K 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 one-item 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 in-depth 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 weight-averaging the fill rates of type-K 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 type-K pure systems to derive the order fill rate.

## ACKNOWLEDGEMENTS

This work was supported by the Ministry of Education of the Republic of Korea and the National Research Foundation of Korea (NRF-2016S1A5A2A01022007).

## Figure

The transition rate diagram and balance equations.

The pseudo-code of GISA.

The convergent speeds.

## Table

A summary of preconditioned iterative methods

A summary of the studies that focused on the performance evaluation

Comparison of GISA with the Gaussian elimination and computer simulation

The fill rates generated by the GISA with different initial solutions

Cases 1 and 2: values of dp and associated qK

Cases 3 and 4: values of dp and associated qK

The fill rates of case 1

The fill rates of case 2

The fill rates of case 3

The fill rates of case 4

## REFERENCES

1. Athappilly K. , Razi M.A. , Tarn J.M. (2010) A multi-technique data mining approach to exploring consumer behaviors , Hum. Syst. Manag, Vol.29 (3) ; pp.153-163
2. Benjaafar S. , Hafsi EI (2006) Production and inventory control of a single product assemble-toorder system with multiple customer classes , Manage. Sci, Vol.52 (12) ; pp.1896-1912
3. Dayanik S. , Song J.S. , Xu S.H. (2003) The effectiveness of several performance bounds for capacitated production, partial-order-service, assemble-toorder systems , Manuf. Serv. Oper. Manag, Vol.5 (3) ; pp.230-251
4. Gao C. , Shen H. , Cheng T.C. (2010) Orderfulfillment performance analysis of an assemble-toorder system with unreliable machines , Int. J. Prod. Econ, Vol.126 (2) ; pp.341-349
5. Gunawardena A.D. , Jain S.K. , Snyder L. (1991) Modified iterative methods for consistent linear systems , Linear Algebra Appl, Vol.154-156 ; pp.123-143
6. Hadjidimos A. , Noutsos D. , Tzoumas M. (2003) More on modifications and improvements of classical iterative schemes for M-matrices , Linear Algebra Appl, Vol.364 ; pp.253-279
7. Hoen K.M. , Gullu R. , van Houtum G.J. , Vliegen I.M. (2011) A simple and accurate approximation for the order fill rates in lost-sales assembleto- order systems , Int. J. Prod. Econ, Vol.133 (1) ; pp.95-104
8. Iravani S. , Luangkesorn L. , Simchi-Levi D. (2003) On assemble-to-order systems with flexible customers , IIE Trans, Vol.35 (5) ; pp.389-403
9. Iravani S. , Luangkesorn L. , Simchi-Levi D. (2004) A general decomposition algorithm for parallel queues with correlated arrivals , Queueing Syst, Vol.47 (4) ; pp.313-344
10. Kohno T. , Kotakemori H. , Niki H. (1997) Improving the modified Gauss-Seidel method for Z-matrices , Linear Algebra and its Application, Vol.267 ; pp.113-123
11. Kotakemori H. , Harada K. , Morimoto M. , Niki H. (2002) A comparison theorem for the iterative method with the preconditioner (I + Smax) , J. Comput. Appl. Math, Vol.145 (2) ; pp.373-378
12. Kotakemori H. , Niki H. , Okamoto N. (1996) Accelerated iterative method for Z-matrices , J. Comput. Appl. Math, Vol.75 (1) ; pp.87-97
13. Lu Y. , Song J.S. , Yao D.D. (2003) Order fill rate, lead time variability, and advance demand information in an assemble-to-order system , Oper. Res, Vol.51 (2) ; pp.292-308
14. Milaszewicz J.P. (1987) Improving Jacobi and Gauss– Seidel iterations , Linear Algebra Appl, Vol.93 ; pp.161-170
15. Morimoto M. , Harada K. , Sakakihara M. , Sawami H. (2004) The Gauss–Seidel iterative method with the preconditioning matrix (I +S +Sm) , Jpn. J. Ind. Appl. Math, Vol.21 (1) ; pp.25-34
16. Niki H. , Kohno T. , Morimoto M. (2008) The preconditioned Gauss-Seidel method faster than the SOR method , J. Comput. Appl. Math, Vol.219 (1) ; pp.59-71
17. Noutsos D. , Tzoumas M. (2006) On optimal improvements of classical iterative schemes for Zmatrices , J. Comput. Appl. Math, Vol.188 (1) ; pp.89-106
18. Park C. , Seo J. (2013) Consideration of purchase dependence in inventory management , Comput. Ind. Eng, Vol.66 (2) ; pp.274-285
19. Song J.S. (1998) On the order fill rate in a multi-item, base-stock inventory system , Oper. Res, Vol.46 (6) ; pp.831-845
20. Song J.S. (2002) Order-based backorders and their implications in multi-item inventory systems , Manage. Sci, Vol.48 (4) ; pp.499-516
21. Song J.S. , Yao D.D. (2002) Performance analysis and optimization of assemble-to-order systems with random lead times , Oper. Res, Vol.50 (5) ; pp.889-903
22. Song J.S. , Zipkin P. de Kok T. , Graves S. (2003) Handbooks in Operations Research and Management Science, North-Holland, ; pp.561-596
23. Song J.S. , Xu S. , Liu B. (1999) Order-fulfillment performance measures in an assemble-to-order system with stochastic lead times , Oper. Res, Vol.47 (1) ; pp.131-149
24. Usui M. , Niki H. , Kohno T. (1994) Adaptive Gauss-Seidel method for linear systems , Int. J. Comput. Math, Vol.51 (1-2) ; pp.119-125
25. Vliegen I.M. , van Houtum G.J. (2009) Approximate evaluation of order fill rates for an inventory system of service tools , Int. J. Prod. Econ, Vol.118 (1) ; pp.339-351
26. Yang Y. , Liu H. , Cai Y. (2013) Discovery of online shopping patterns across websites , INFORMS J. Comput, Vol.25 (1) ; pp.161-176
27. Ye J. , Li S. (1994) Folding algorithm: A computational method for finite QBD processes with leveldependent transitions , IEEE Trans. Commun, Vol.42 (234) ; pp.625-639
28. Zheng B. , Miao S.X. (2009) Two new modified Gauss-Seidel methods for linear system with Mmatrices , J. Comput. Appl. Math, Vol.233 (4) ; pp.922-930
29. Zhou W. , Chao X. (2012) Stein-Chen approximation and error bounds for order fill rates in assembleto- order systems , Nav. Res. Logist, Vol.59 (8) ; pp.643-655