1.Introduction
For companies in today’s competitive marketplace, order fulfillment performance has become more and more imperative in order to respond quickly to customer needs [14, 20]. This paper considers the calculation of order fill rate when a customer order is composed of several items and the customer order is met only if all the items can be deliverable at once. The order fill rate is defined as the probability of satisfying a customer order immediately using the existing inventory.
Song [16] stated the order fill rate as the orderbased performance measure. However, most standard inventory models do not consider connections between items. The models assume that each item demand is independent. This case is referred as the itembased approach. In the orderbased ap proach, the following question is important, because customer order requires several items at the same time: What is the probability that a customer order can be met immediately, given any safety stock level?
Park and Seo [13] faced the same question when analyzing the inventory of a spare parts distributor for ship engines and generators. After investigating the historical purchaseorder request (POR) data, Park and Seo [13] identified that a many of orders had been cancelled because of the lack of a few spare parts. That is, if an order contained at least one item that was out of stock, the whole order was cancelled, even though the other items had large quantities of stocks. They demonstrated that over 30% was explained by the portion of orders cancelled because of stockedout parts. In most of cases, shipping companies purchase spare parts only when all of the spare parts are available immediately, because they want to buy all of the spare parts together. In this case, if one of the spare parts is out of stock, while all of the other parts are in stock, the situation is the same as if all items were out of stock. Park and Seo [13] called this type of dependence as ‘purchase dependence.’
Purchase dependence is different to demand dependence. Purchase dependence treats the purchase behavior of customers, while demand dependence considers demand correlation between items, between regions, or over time. Purchase dependence can be observed in such areas as marketing, manufacturing systems, and distribution systems.
This paper proposes a new approximate calculation of order fill rate under purchase dependence. We calculate the order fill rate using the stationary joint distribution of onhand inventories. The stationary joint distribution is modelled by a quasibirthanddeath (QBD) process and it consists of the balance equations and a normalization equation. However, traditional computational methods for solving the linear equations face with a difficulty of the curse of dimensionality for the large number of equations and unknowns. Thus we develop the greedy iterative search algorithm (GISA).
This paper is organized as follows. In the following section, we review the related literature. Section 3 explains the model that calculates the order fill rate. In Section 4, we depict the procedure of GISA. Section 5 conducts the computational analysis. In Section 6, we finalize this paper by providing our conclusions.
2.Literature Review
Firstly, Song [16] examined the orderbased approach. She derived an exact solution of the order fill rate for a multiitem, multiproduct system with backlogging and constant lead times. A customer order was composed of several items in different amounts. Using a series of convolution of onedimensional distributions, she calculated the order fill rate.
On the other hand, Song, et al. [17] studied a multicomponent, multiproduct system, where each of the components is managed by the maketostock and refilled by the exponential processing machine. They formulated a generalized model having both complete backlogging and lost sales as a special case. In addition, they differentiated total order service and partial order service. The total order service is that an order is fulfilled completely or rejected as a whole, while the partial order service is that partial fulfilment occurs. Iravani et al. [4] and Gao et al. [1] examined a similar situation. However, they considered a little different perspective. Iravani et al. [4] considered customer flexibility. The customer flexibility means that customers change some of their requested components not in stock. On the other hand, Gao et al. [1] considered that each item was refilled by an independent unreliable machine.
One major limitation with the above studies is their inability to analyze large systems. Their modelling approaches are difficult to calculate, because they require incorporating correlations among demands for different items. Even though the matrixgeometric techniques are used by Song et al. [17], Iravani et al. [4] and Gao et al. [1], the curse of dimensionality is still remained.
In order to solve a linear system of equations : Ax = b, A = (a_{ij}) ∈ R^{mxm}, x, b ∈ R^{m}, we usually use a wellknown GaussSeidel method. The GaussSeidel method is modified from the Jacobi method. Thus it needs less iteration to make the same degree of accuracy. In order to speedup the convergence rate, several preconditioned iterative methods have been developed. These studies are summarized in <Table 1>. However, our linear equations system is not solved by the GaussSeidel method and several preconditioned iterative methods, because these methods diverge.
3.Model Description
When purchase dependence exists, we consider a retailer system that has J different items. $\Omega $ = {1, 2, …, J} is the set of all item indices. For any $K\subseteq \Omega $, let K denote the number of elements in K. The overall order process is stationary in time and follows a Poisson process. Each customer requests at most one unit of each item but may require several items simultaneously. For any $K\subseteq \Omega $, an order is of type K if it requires one and only one unit of each item in K and 0 units in $\Omega $  K. Let T = {1, 2, …, T} be the set of order types (1 ≤ T ≤ 2^{J}1). We assume that the probability of type K is constant. Every order type is independent each other. This means that the demand process for each item follows also a Poisson process. This paper uses the same assumptions and notations defined in Song et al. [17].
Orders are satisfied by a FirstComeFirstServe (FCFS) basis. When an order comes and a retailer has some items in outofstock, the order is cancelled. Song et al. [17] called this situation as a total order service. Each item is managed independently by an base stock policy. Let s_{i} be the base stock level for item i. This means that if the inventory position of item i is less than si, then order up to s_{i} at each demand epoch. The inventory position is calculated by the inventory on hand plus inventory on order. It is assumed that item i has an inventory at level s_{i} for all i at time 0. Then, if the item i is transported to the customer, each demand for item i generates an order for that item. Hence there can exist (at most) s_{i} outstanding orders of item i. Orders for item i are refilled by a supplier i and the lead times are distributed exponentially with the rate ${\mu}_{i}$. The supplier processes the orders on a FCFS basis.
3.1.The Order Fill Rate
For any i ∈ Ω and K ∈ T, let $\lambda $ be the overall order rate, q^{K} be the probability that an order is of type K, λ^{K} be the order rate of demand type K (i.e., λ^{K} = q^{K}$\lambda $), 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}$\lambda $).
We denote that
the fill rate of typeK order is denoted by

F^{K} = the joint probability that all items in a typeK order are delivered immediately
$=P\left({I}_{i}>0;\hspace{0.17em}i\in K\right)$
Then, the order fill rate of a retailer system F is obtained by(1)
3.2.The Stationary Joint Distribution
Due to purchase dependence, a typeK order is cancelled if at least one item i (i ∈ K) is out of stock. Since the fill rate of typeK order is the joint probability that all items in the typeK order are delivered immediately, the fill rate of typeK order is decided by the joint distribution of (I_{i} : i ∈ K). In this section, we demonstrate how to derive the joint distribution of (I_{i} : i ∈ K).
A stochastic process {(I_{i}(t) : i ∈ K), t ≥ 0}, is a continuous time Markov chain with finite state space {n = (n_{i} : i ∈ K)0 ≤ n_{i} ≤ s_{i}, i ∈ K}, where I_{i}(t) is the inventory on hand of item i at time t. There are $\prod _{i\in K}\left({s}_{i}+1\right)$ states. A state transition can only take place if a demand or an item on order arrives. In other word, with transition rate ${\lambda}^{K}$, the state n transits to the state ${\mathbf{n}}^{\prime}=\left({{n}^{\prime}}_{i}:\hspace{0.17em}i\in K\right)$, where
and with transition rate ${\mu}_{i}$ (i ∈ K), the state n transits to the state ${n}^{\prime}=\left({{n}^{\prime}}_{i}:\hspace{0.17em}i\in K\right)$, where
It is easy to see that this Markov chain is irreducible, so its stationary distribution exists uniquely [17]. The balance equation for each state can be formulated by utilizing equations (2) and (3). <Figure 1> shows the transition rate diagram and balance equations for the retailer system in which the type1 order requires an item 1, the type2 order requires an item 2, and the type3 order requires both items 1 and 2. In the equilibrium situation, it should be that the rate of leaving any node equals to the rate of entering that node. Using the ‘rate in = rate out’ principle, we can make the following generalization.
Let I be a vector of (I_{i} : i ∈ K) and p_{I} be the equilibrium probabilities. Then
where ${\delta}_{i}$ = 1 for I_{i} < s_{i} or 0 otherwise, ${\xi}_{K}$ = 1 for ∀I_{i} > 0 (i ∈ K) or 0 otherwise, ei is a unit vector having 1 at the position of I_{i}, and p_{I} = 0 for any I_{i} ∉ [0,s_{i}].
Then, the stationary joint distribution can be obtained by solving the balance equations and a normalization equation. In order to solve the linear equations, we can use traditional computational algorithms such as Gaussian elimination. However, the direct approach can face with severe computation and space complexities for the large size of problem. Voiding the computation and space complexities, Song et al. [17], Iravani et al. [4] and Gao et al. [1] used the special structure of the QBD process. They obtained a unified, matrixgeometric solution of the stationary joint distribution. On the other hand, Ye and Li [19] developed the folding algorithm to conduct a steady state analysis of finite QBD process. The curse of dimensionality is still remained. Thus we develop the greedy iterative search algorithm to avoid the curse of dimensionality. The algorithm utilizes the special structure of the QBD process shown in equation (4) and stops the solution from diverging by adding the newly invented steps into the GaussSeidel method.
4.Greedy Iterative Search Algorithm
A linear equations system is established by the balance equations depicted by equation (4) and the normalization equation such that Ax = b, A = (a_{ij}) ∈ R^{mxm}, x, b ∈ R^{m}. Here, we say that Ax is the coefficient matrix, b is the righthand side vector, x is the vector of unknowns (i.e., pI, the equilibrium probabilities), and the order $m={\displaystyle \prod _{i\in K}\left({s}_{i}+1\right)}$. The linear equations system, sometimes, can be large and sparse. Then, a serious storage problem can be brought by the large number of equations and unknowns. Gaussian elimination is a very accurate, economical, and useful algorithm, if possible. But if the order m is so large and the results of Gaussian elimination cannot be stored, it would be better to solve the linear equations system by other ways that never change the matrix A and never require storing more than a few vectors of length m. For this purpose, iterative methods are recommended.
In this section, we describe the procedure of GISA. We begin with an initial solution vector x^{(0)}, and generate a sequence of solution vectors x^{(1)} → x^{(2)} → ⋯ according to the algorithm. We hope that as k → ∞, x^{(k)} will converge to the exact solution. <Figure 2> shows the pseudocode of GISA. In the initialization step, an initial solution x^{(0)} set up and we compute the maximum absolute residual r_{max} and the summation of absolute residual r_{sum} corresponding to x^{(0)}. Let [y]_{i} be the ith element of vector y.
Then we create the solution vectors x^{(k)} according to the following rule. Let ${x}_{i}^{\left(k\right)}$ denote the ith element of the kth iterate solution vector x^{(k)} (i.e., [x^{(k)}]_{i}) and b_{i} the ith element of the righthand side b (i.e., [b]_{i}). The algorithm adjusts the ith element of the current solution, in the order i = 1, 2, …, m, to abolish the ith element of the residual. The solution is updated immediately if the newly computed element ${x}_{i}^{\left(k\right)}$ is nonnegative, and the newly computed maximum absolute residual or the newly computed summation of absolute residual does not increase. These steps are invented in this paper to prevent the solution from diverging.
The algorithm continues to create solution vectors x^{(k)} until a convergence is achieved. A convergence is reached when the iteration k becomes to the maximum number of iterations; the summation of absolute residual r_{sum} is less than the maximum value of tolerance; or no longer change is made in the solution vector x^{(k)}. The solution vectors x^{(k)} would be affected by an initial solution x^{(0)} as the order m becomes so large. Section 5.2 describes the impacts of initial solution.
5.Computational Analysis
In this section, we illustrate numerical and simulation experiments that are conducted to demonstrate the GISA performance, examine the influences of an initial solution and a relaxation method, and obtain the managerial insights.
5.1.The Performance of GISA
The GISA was developed to derive the stationary joint distribution of onhand inventories in the retailer system and calculate the order fill rate of retailer system. The order fill rate is calculated by the probability that a customer order can be met immediately. We measure the performance of GISA by the order fill rate of retailer system.
For the purpose of testing the GISA performance, the test bed was designed as follows. Combining the number of items and the number of order types, four cases are chosen. Case 1 is that J = 3 with T = 7. Case 2 is that J = 3 with T = 7 (asymmetric). Case 3 is that J = 3 with T = 4. Case 4 is that J = 5 with T = 6. All possible order types occur in cases 1 and 2, and thus T = 7. Type 1 – 3 requires each item, respectively, type 4 – 6 requires combination of two items, respectively, and type 7 requires all items. We set that q^{1} = q^{2} = q^{3} = 0.05, q^{4} = q^{5} = q^{6} = 0.07, and q^{7} = 0.64 for case 1, and q^{1} = 0.04, q^{2} = 0.05, q^{3} = 0.06, q^{4} = 0.08, q^{5} = 0.06, q^{6} = 0.07, and q^{7} = 0.64 for case 2. Cases 3 and 4 consider order types for individual items and one order type for all items, so T = J + 1. We set that q^{1} = q^{2} = q^{3} = q^{4} = 0.05 and q^{5} = 0.85 for case 3, and q^{1} = q^{2} = q^{3} = q^{4} = q^{5} = 0.05 and q^{6} = 0.75 for case 4. We treat order types symmetric for the order rate, however, we dealt with asymmetric order types for case 2.
Within all retailer systems, equal base stock levels were chosen, i.e., ${s}_{1}={s}_{2}=\cdots ={s}_{J}$. For these base stock levels, we assigned three values of 5, 10, and 15. We set that all replenishment rates ${\mu}_{i}$ and the overall order rate for the retailer system $\lambda $ are equal to 1.0.
For the purpose of comparison, we obtained the exact order fill rate by utilizing the Gaussian elimination. However, the Gaussian elimination was not suitable for the large cases requiring the large matrix manipulation. In this case, a simulation program was developed to calculate the order fill rate. We implemented the simulation run for 5 times per each case. For each time, the simulation was implemented for 10,000 customer orders with the warmup of 1,000 orders. Then the order fill rate was averaged the 5 order fill rates resulted from simulations.
The results of this comparison are presented in <Table 2>. For small cases, the order fill rates of retailer system calculated by the algorithm are almost same as those by the Gaussian elimination. This demonstrates the GISA performance. We can notice the fact by just reviewing that r_{max} and r_{sum} are near to zeros. For large cases, even though r_{max} and r_{sum} demonstrate the performance of GISA, <Table 2> shows the comparison between the algorithm and computer simulation. The average value of absolute difference is 0.0030. Based on the results in <Table 2>, the GISA can be considered as a reliable algorithm to obtain the stationary joint distribution of onhand inventories in the retailer system.
5.2.The Impacts of an Initial Solution and a Relation Method
The initial solutions are distinguished by assigning a different number of elements having a positive value in the initial solution vector. That is, we allocated a positive value to the only controlled elements in the initial solution vector, but zeros to the rest. We calculated an equal positive value as dividing 1.0 by the number of positive elements.
In order to investigate the effect of initial solutions on the convergent speed of the GISA, we computed the order fill rates of case 4 in Section 5.1 as increasing the number of positive elements in the initial solution vector (e.g., 10, 50, 100, and 150). For the two cases of s_{i} = 5 and s_{i} = 10, <Figure 3> illustrates the convergent speeds. (This section sets that $\lambda $ = 1.5 and ${\mu}_{i}$ = 1.0 (i ∈ K).) From <Figure 3>, we can observe that, in general, there is no significant difference in the convergent speed depending on the initial solutions. (In the case of s_{i} = 10, even though the initial solution having 10 positive values reaches to the convergence faster than others, the accuracy is very poor.) However, it is observed that the more number of the base stock level requires the more iteration to reach to the convergence.
Meanwhile, we can recognize that the solution vectors will converge to a certain solution, because the GISA algorithm provides the steps of stopping the solution from diverging. However, it is observed that for the large cases, the GISA algorithm could converge to different solutions depending on the initial solutions, which is shown in <Table 3>. <Table 3> shows the order fill rates of the case of s_{i} = 10 in <Figure 3>. Among the 4 cases, the order fill rate of the initial solution having 50 positive values results in the highest degree of accuracy, i.e., a minimum value of r_{sum}.
It is wellknown that the convergence of GaussSeidel method may sometimes be improved by combining with relaxation [10, 15]. The relaxation method modifies the Gauss Seidel procedure in <Figure 2>, as given by
For choices of ω, the procedure is called as either an under relaxation method with 0 < ω < 1 or an overrelaxation method with 1 < ω.
Like to the examination on the convergent speed of the GISA algorithm, we computed the order fill rates of case 4 for various values of ω. For the two cases of s_{i} = 5 and s_{i} = 10, <Figure 4> illustrates the convergent speeds. (This section sets that $\lambda $ = 1.5 and ${\mu}_{i}$ = 1.0 (i ∈ K).) We show only the results of underrelaxation methods because the results of overrelaxation methods are too poor. From <Figure 4>, we can observe that the relaxation methods do not improve the convergence speed of the GISA.
5.3.Analysis of Purchase Dependence
This paper calculates the order fill rate of a retailer system considering purchase dependence. Purchase dependence implies that the order is cancelled when an order comes and a retailer has some items in outofstock. In this section, we analyze the order fill rate as changing the degree of purchase dependence, thus deduce the managerial insights.
For analysis purpose, we utilized the case 1 in Section 5.1 as the test bed. In case 1, all base stock levels were assigned equally, i.e., ${s}_{1}={s}_{2}={s}_{3}=5$. We set that all replenishment rates ${\mu}_{i}$ and the overall order rate for the retailer system $\lambda $ are equal to 1.0.
The degree of purchase dependence (dp) is defined as
The degree of purchase dependence dp = 1 if all orders are the one order type for the set of all items. Inversely, if all orders are different order types for individual items only, dp = 0. The dp assesses how much demands are connected each other. <Table 4> shows the values of dp and associated q^{K} for case 1.
In order to investigate the attributes of order fill rate, we decomposed the order fill rate into the fill rate of a pure system. The pure system has one order type K. (We called the system as the typeK pure system.) The fill rate of typeK pure system is obtained by decomposing customer orders into item demands, in other words, by building new item demand processes. The new demand processes follow the same marginal distributions as the original order processes. For the typeK pure system, the overall order rate is approximated by averaging the demand rates of items in the typeK order. Let

${\lambda}^{pK}$ = the overall order rate for the typeK pure system
$=\frac{1}{\leftK\right}{\displaystyle \sum _{i\in K}{\lambda}_{i}}$.
Then, the fill rate of typeK pure system is calculated by

${F}^{pK}$ = the joint probability of filling immediately all items in the typeK pure system
$=P\left({I}_{i}>0\hspace{0.17em}:\hspace{0.17em}i\in K\right)$
We calculate the fill rates of typeK pure systems using the GISA after decomposing customer orders into item demands.
<Table 5> and <Figure 5> show the fill rates of both typeK pure systems and the retailer system obtained by the GISA. We can recognize some managerial insights. First, the value of order fill rate poses between values of tow pure systems. One pure system consists of all items. The other consists of only one item. That is, these two fill rates are the lower and upper bounds for the order fill rate. The lower bound can be calculated by the fill rate of pure system having all items. The upper bound can be calculated by the oneitem pure system.
Secondly, as the degree of purchase dependence declines while other conditions remain same, it is observed that the difference between the lower and upper bounds reduces, the order fill rate increases, and the order fill rate gets closer to the upper bound. From these results, it can be inferable that if we deal with purchase dependent item demands independently, we would face either unnecessary overstocking of some items or unsatisfactory service at the order level. For better profit and lower inventory costs, treating knowledge of purchase dependence is important in designing an inventory replenishment policy.
6.Conclusion
For the retailer system having purchase dependence, we developed a new approximate calculation of the order fill rate. The order fill rate is the probability of meeting an customer order immediately from existing inventory. During the derivation of the stationary joint distribution, traditional computational methods face with the curse of dimensionality for the large problems.
This paper developed the GISA algorithm to elude the curse of dimensionality and keep the solution from diverging. The GISA algorithm is grounded on the GaussSeidel method. From the comparison analysis of the GISA and the simulation, this paper illustrated that the GISA is a reliable algorithm to gain the stationary joint distribution of onhand inventories in the retailer system.
In addition, we observed some managerial insights. (1) The upper bound of order fill rate can be calculated by the oneitem pure system, while the lower bound can be provided by the pure system that consists of all items. (2) As the degree of purchase dependence declines while other conditions remain same, it is observed that the difference between the lower and upper bounds reduces, the order fill rate increases, and the order fill rate gets closer to the upper bound.
As the degree of purchase dependence declines, the order fill rate gets closer to the upper bound created by the pure system. It is wondered when treating purchase dependent demands as independent is tolerable. We leave this theme as the future studies.