Journal Search Engine

ISSN : 2005-0461(Print)
ISSN : 2287-7975(Online)

Journal of Society of Korea Industrial and Systems Engineering Vol.40 No.4 pp.137-146
DOI : https://doi.org/10.11627/jkise.2017.40.4.137

Approximate Calculation of Order Fill Rate under Purchase Dependence

Changkyu Park†

College of Business, University of Ulsan

Corresponding Author : ckparkuou@ulsan.ac.kr

Received October 10, 2017 Finally Revised November 8, 2017 Accepted November 9, 2017

Abstract

For the competitive business environment under purchase dependence, this paper proposes a new approximate calculation of order fill rate which is a probability of satisfying a customer order immediately using the existing inventory. 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. Traditional computational methods have a difficulty of the curse of dimensionality for the large cases, when deriving the stationary joint distribution which is utilized to calculate the order fill rate. In order to escape the curse of dimensionality and protect the solution from diverging for the large cases, we develop a greedy iterative search algorithm based on the Gauss-Seidel method. We show that the greedy iterative search algorithm is a dependable algorithm to derive the stationary joint distribution of on-hand inventories in the retailer system by conducting a comparison analysis of a greedy iterative search algorithm with the simulation. In addition, we present some managerial insights such as : (1) The upper bound of order fill rate can be calculated by the one-item 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.

Key Words : Approximate Calculation , Order Fill Rate , Purchase Dependence , Quasi-Birth-and-Death Process

구매종속성을 고려한 주문충족률의 근사적 계산

박창규^†

울산대학교 경영대학

초록

키워드 :

This article has been cited by 0 article in crossref

Cited-By

Funding:

the 2016 Research Fund of University of Ulsan.

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 order-based 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 item-based approach. In the order-based 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 stocked-out 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 quasi-birth-and-death (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 order-based approach. She derived an exact solution of the order fill rate for a multiitem, multi-product 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 multi-component, multi-product system, where each of the components is managed by the make-to-stock 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 matrix-geometric 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 well-known Gauss-Seidel method. The Gauss-Seidel method is modified from the Jacobi method. Thus it needs less iteration to make the same degree of accuracy. In order to speed-up 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 Gauss-Seidel 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. $Ω$ = {1, 2, …, J} is the set of all item indices. For any $K \subseteq Ω$ , 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 Ω$ , 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|} 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 First-Come-First-Serve (FCFS) basis. When an order comes and a retailer has some items in out-of-stock, 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 $μ_{i}$ . The supplier processes the orders on a FCFS basis.

3.1. The Order Fill Rate

For any i ∈ Ω and K ∈ T, let $λ$ 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 $λ$ ), q_i be the probability that an order requires item i (i.e., $q_{i} = \sum_{K : i \in K} q^{K}$ ), and λ_i be the aggregate demand rate of item i (i.e., λ_i = q_i $λ$ ).

We denote that

I_i = the on-hand inventory of item i in steady-state, 0 ≤ I_i ≤ s_i,

the fill rate of type-K order is denoted by

F^K = the joint probability that all items in a type-K order are delivered immediately

$= P (I_{i} > 0; i \in K)$

Then, the order fill rate of a retailer system F is obtained by(1)

F = \sum_{K \in T} q^{K} F^{K}

(1)

3.2. The Stationary Joint Distribution

Due to purchase dependence, a type-K order is cancelled if at least one item i (i ∈ K) is out of stock. Since the fill rate of type-K order is the joint probability that all items in the type-K order are delivered immediately, the fill rate of type-K 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} (s_{i} + 1)$ states. A state transition can only take place if a demand or an item on order arrives. In other word, with transition rate $λ^{K}$ , the state n transits to the state $n^{'} = ({n^{'}}_{i} : i \in K)$ , where

{n^{'}}_{j} = {\begin{matrix} n_{j} - 1 & if j \in K and 0 < n_{i} for all i \in K \\ n_{j} & o t h e r w i s e \end{matrix}

(2)

and with transition rate $μ_{i}$ (i ∈ K), the state n transits to the state $n^{'} = ({n^{'}}_{i} : i \in K)$ , where

{n^{″}}_{j} = {\begin{matrix} n_{j} + 1 & if j = i and n_{i} < s_{i} \\ n_{j} & o t h e r w i s e \end{matrix}

(3)

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 type-1 order requires an item 1, the type-2 order requires an item 2, and the type-3 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

(\sum_{i \in K} δ_{i} μ_{i} + \sum_{K \in T} ξ_{K} λ^{K}) p_{I} = \sum_{i \in K} μ_{i} p_{I - e_{i}} + \sum_{K \in T} λ^{K} p_{I + \sum_{i \in K} e_{i}}

(4)

where $δ_{i}$ = 1 for I_i < s_i or 0 otherwise, $ξ_{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, matrix-geometric 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 Gauss-Seidel 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 right-hand side vector, x is the vector of unknowns (i.e., pI, the equilibrium probabilities), and the order $m = \prod_{i \in K} (s_{i} + 1)$ . 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⁽⁰⁾, and generate a sequence of solution vectors x⁽¹⁾ → x⁽²⁾ → ⋯ according to the algorithm. 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, an initial solution x⁽⁰⁾ set up and we compute the maximum absolute residual r_max and the summation of absolute residual r_sum corresponding to x⁽⁰⁾. 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}^{(k)}$ 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 right-hand 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}^{(k)}$ 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⁽⁰⁾ 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 on-hand 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¹ = q² = q³ = 0.05, q⁴ = q⁵ = q⁶ = 0.07, and q⁷ = 0.64 for case 1, and q¹ = 0.04, q² = 0.05, q³ = 0.06, q⁴ = 0.08, q⁵ = 0.06, q⁶ = 0.07, and q⁷ = 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¹ = q² = q³ = q⁴ = 0.05 and q⁵ = 0.85 for case 3, and q¹ = q² = q³ = q⁴ = q⁵ = 0.05 and q⁶ = 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} = \dots = s_{J}$ . For these base stock levels, we assigned three values of 5, 10, and 15. We set that all replenishment rates $μ_{i}$ and the overall order rate for the retailer system $λ$ 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 warm-up 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 on-hand 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 $λ$ = 1.5 and $μ_{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 well-known that the convergence of Gauss-Seidel 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

x_{i}^{(k)} = \frac{1}{a_{i i}} (b_{i} - \sum_{j = 1}^{i - 1} a_{i} x_{j}^{(k)} - \sum_{j = i + 1}^{n} a_{i} x_{j}^{(k - 1)})

(5)

\begin{matrix} x_{i}^{(k)} = (1 - ω) x_{i}^{(k - 1)} + \\ \frac{ω}{a_{i i}} (b_{i} - \sum_{j = 1}^{i - 1} a_{i} x_{j}^{(k)} - \sum_{j = i + 1}^{n} a_{i} x_{j}^{(k - 1)}) \end{matrix}

(6)

For choices of ω, the procedure is called as either an under- relaxation method with 0 < ω < 1 or an over-relaxation 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 $λ$ = 1.5 and $μ_{i}$ = 1.0 (i ∈ K).) We show only the results of under-relaxation methods because the results of over-relaxation 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 out-of-stock. 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 $μ_{i}$ and the overall order rate for the retailer system $λ$ are equal to 1.0.

The degree of purchase dependence (dp) is defined as

d p = \sum_{K \in T} \frac{| K | - 1}{J - 1} q^{K}

(7)

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 type-K pure system.) The fill rate of type-K 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 type-K pure system, the overall order rate is approximated 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

$= \frac{1}{| K |} \sum_{i \in K} λ_{i}$ .

Then, the fill rate of type-K pure system is calculated by

$F^{p K}$ = the joint probability of filling immediately all items in the type-K pure system

$= P (I_{i} > 0 : i \in K)$

We calculate the fill rates of type-K pure systems using the GISA after decomposing customer orders into item demands.

<Table 5> and <Figure 5> show the fill rates of both type-K 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 one-item 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 Gauss-Seidel 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 on-hand 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 one-item 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.

Acknowledgement

This work was supported by the 2016 Research Fund of University of Ulsan.

Figure

<Figure 1> .The Transition Rate Diagram and Balance Equations

<Figure 2> .The Pseudo-Code of GISA

<Figure 3> .The Convergent Speeds Depending on Initial Solutions

<Figure 4> .The Convergent Speeds of Under-Relaxation Method

<Figure 5> .Fill Rates Against the dp

Table

<Table 1>.

A Summary of Preconditioned Iterative Methods

Research	Preconditioner	Remark
Milaszewicz [8]	I +C	C eliminates the elements of the first column below the diagonal of A.
Gunawardena et al. [2]	I +S	S transforms the first upper codiagonal to zero.
Usui et al. [18]	I +U	U is a strictly upper triangular part of-A.
Kotakemori et al. [7]	I+βU	β≥ 1
Kohno et al. [5]	I+αS	α>1
Kotakemori et al. [6]	I+Smax	Smax is constructed by only the largest element at each row of the upper triangular part of A.
Hadjidimos et al.[3]	I+αC
Morimoto et al.[9]	I+S+Smax
Noutsos and Tzoumas[12]	Generalized preconditioner	Eliminate two or more off-diagonal elements in each row.
Niki et al. [11]	I+S+R	R is referred to as the nth preconditioner.
Zheng and Miao [21]	I+ Smax+R
I+ Smax+Rmax	Rmaxis constructed by the largest element of R.

<Table 2>.

Comparison of GISA with the Gaussian Elimination and Computer Simulation

Case	si	Gaussian	GISA	Simulation	Absolute
1	5	0.8476	0.8476	0.0000	0.0000	0.8448	0.0028
10	-	0.9474	0.0000	0.0000	0.9483	0.0009
15	-	0.9794	0.0001	0.0005	0.9796	0.0003
2	5	0.8475	0.8474	0.0000	0.0001	0.8462	0.0012
10	-	0.9471	0.0001	0.0008	0.9468	0.0003
15	-	0.9789	0.0001	0.0008	0.9790	0.0002
3	5	0.8098	0.8098	0.0000	0.0000	0.8142	0.0044
10	-	0.9155	0.0001	0.0006	0.9221	0.0066
15	-	0.9556	0.0000	0.0005	0.9616	0.0060
4	5	-	0.8407	0.0001	0.0006	0.8477	0.0070
10	-	0.9478	0.0001	0.0005	0.9526	0.0048
15	-	0.9852	0.0000	0.0003	0.9838	0.0014

<Table 3>.

The Fill Rates Generated by the GISA with Different Initial Solutions

Initial solution	F	rmax	rsum
10	0.7704	0.0007	0.0065
50	0.7670	0.0000	0.0001
100	0.7665	0.0002	0.0008
150	0.7663	0.0003	0.0012

<Table 4>.

Case 1 Values of dp and Associated q^K

K	Item	dp	qK
1	1	0	0		0.33	0.20	0.15	0.10	0.10	0.05	0.00	0.00	0.00
2	0	1	0	0.33	0.20	0.15	0.10	0.10	0.05	0.00	0.00	0.00
3	0	0	1	0.33	0.20	0.15	0.10	0.10	0.05	0.00	0.00	0.00
4	1	1	0	0.00	0.10	0.15	0.15	0.10	0.07	0.10	0.05	0.00
5	1	0	1	0.00	0.10	0.15	0.15	0.10	0.07	0.10	0.05	0.00
6	0	1	1	0.00	0.10	0.15	0.15	0.10	0.07	0.10	0.05	0.00
7	1	1	1	0.00	0.10	0.10	0.25	0.40	0.64	0.70	0.85	1.00

<Table 5>.

Fill Rates of Type-K Pure and Whole Systems

dp	Fp7	Fp4=Fp5= Fp6	Fp1= Fp2= Fp3	F
0.000	0.9930	0.9950	0.9973	0.9973
0.250	0.9644	0.9729	0.9841	0.9781
0.325	0.9498	0.9613	0.9767	0.9658
0.475	0.9130	0.9328	0.9561	0.9318
0.550	0.8915	0.9147	0.9429	0.9130
0.745	0.8297	0.8569	0.9005	0.8476
0.850	0.7950	0.8246	0.8740	0.8046
0.925	0.7703	0.8011	0.8540	0.7757
1.000	0.7460	0.7776	0.83333	0.7460

Reference

C. Gao , H. Shen , T.C. Cheng (2010) Order-fulfillment performance analysis of an assemble-to-order system with unreliable machines., Int. J. Prod. Econ., Vol.126 (2) ; pp.341-349
A.D. Gunawardena , S.K. Jain , L. Snyder (1991) Modified iterative methods for consistent linear systems., Linear Algebra Appl., Vol.154-156 ; pp.123-143
A. Hadjidimos , D. Noutsos , M. Tzoumas (2003) More on modifications and improvements of classical iterative schemes for M-matrices., Linear Algebra Appl., Vol.364 ; pp.253-279
S. Iravani , L. Luangkesorn , D. Simchi-Levi (2003) On assemble-to-order systems with flexible customers., IIE Trans., Vol.35 (5) ; pp.389-403
T. Kohno , H. Kotakemori , H. Niki (1997) Improving the modified Gauss-Seidel method for Z-matrices., Linear Algebra and its Application, Vol.267 ; pp.113-123
H. Kotakemori , K. Harada , M. Morimoto , H. Niki (2002) A comparison theorem for the iterative method with the preconditioner (I+Smax)., J. Comput. Appl. Math., Vol.145 (2) ; pp.373-378
H. Kotakemori , H. Niki , N. Okamoto (1996) Accelerated iterative method for Z-matrices., J. Comput. Appl. Math., Vol.75 (1) ; pp.87-97
J.P. Milaszewicz (1987) Improving Jacobi and Gauss-Seidel iterations., Linear Algebra Appl., Vol.93 ; pp.161-170
M. Morimoto , K. Harada , M. Sakakihara , H. Sawami (2004) The Gauss-Seidel iterative method with the preconditioning matrix (I+S+Sm)., Jpn. J. Ind. Appl. Math., Vol.21 (1) ; pp.25-34
W. Niethammer (2001) The successive order relaxation method (SOR) and Markov chains., Ann. Oper. Res., Vol.103 (1) ; pp.351-358
H. Niki , T. Kohno , M. Morimoto (2008) The preconditioned Gauss-Seidel method faster than the SOR method., J. Comput. Appl. Math., Vol.219 (1) ; pp.59-71
D. Noutsos , M. Tzoumas (2006) On optimal improvements of classical iterative schemes for Z-matrices., J. Comput. Appl. Math., Vol.188 (1) ; pp.89-106
C. Park , J. Seo (2013) Consideration of purchase dependence in inventory management., Comput. Ind. Eng., Vol.66 (2) ; pp.274-285
K.J. Park (2009) Effect of customer demand on total inventory cost and order fill rate in a supply chain., Journal of Society of Korea Industrial and Systems Engineering, Vol.32 (3) ; pp.93-98
B. Philippe , Y. Saad , W.J. Stewart (1992) Numerical methods in Markov chain modeling., Oper. Res., Vol.40 (6) ; pp.1156-1179
J.S. Song (1998) On the order fill rate in a multi-item, basestock inventory system., Oper. Res., Vol.46 (6) ; pp.831-845
J.S. Song , S. Xu , B. Liu (1999) Order-fulfillment performance measures in an assemble-to-order system with stochastic lead times., Oper. Res., Vol.47 (1) ; pp.131-149
M. Usui , H. Niki , T. Kohno (1994) Adaptive Gauss- Seidel method for linear systems., Int. J. Comput. Math., Vol.51 (1-2) ; pp.119-125
J. Ye , S. Li (1994) Folding algorithm : a computational method for finite QBD processes with level-dependent transitions., IEEE Trans. Commun., Vol.42 (234) ; pp.625-639
S.C. Yoon (2015) A study on inventory control method for an item with stock keeping units., Journal of Society of Korea Industrial and Systems Engineering, Vol.38 (1) ; pp.124-130
B. Zheng , S.X. Miao (2009) Two new modified Gauss- Seidel methods for linear system with M-matrices., J. Comput. Appl. Math., Vol.233 (4) ; pp.922-930