Computers & Operations Research 38 (2011) 1584–1595
Contents lists available at ScienceDirect
Computers & Operations Research journal homepage: www.elsevier.com/locate/caor
A multi-class closed queueing maintenance network model with a parts inventory system Chan-Woo Park a,, Hyo-Seong Lee b a b
National Railway Safety R&D Programme Office, Korea Railroad Research Institute, 360-1 Woram-dong, Uiwang-si, Gyeonggi-do 437-757, Korea Department of Industrial & Management Engineering, College of Engineering, Kyung Hee University, Korea
a r t i c l e i n f o
a b s t r a c t
Available online 3 March 2011
The system we address is a maintenance network of repairable items where a set of bases is supported by a centrally located repair depot and a consumable replacement parts inventory system. If an item fails, a replacement part must be obtained at the parts inventory system before the failed item enters the repair depot. The ordering policy for the parts is the (S,Q) inventory policy. An approximation method for this system is developed to obtain performance measures such as steady-state probabilities of the number of items at each site and the expected backorders at the parts inventory system. The proposed system is modelled as a multi-class closed queueing network with a synchronization station and analyzed using a product-form approximation method. Particularly, the product-form approximation method is adapted so that the computational effort on estimating the parameters of the equivalent multi-class network is minimized. In analyzing a sub-network, a recursive method is used to solve balance equations by exploiting the special structure of the Markov chain. Numerical tests show that the approximation method provides fairly good estimation of the performance measures of interests. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Spares provisioning problem Multi-class closed queueing network Product-form approximation method Recursive technique Performance analysis
1. Introduction During the past decades, a substantial number of mathematical models have been developed related to various aspects of managing inventory systems. One of these models, a spares provisioning problem, was designed to determine the stock level of repairable-items such as expensive vehicles, electronic devices, and airplane engines. These repairable items are typically expensive and their individual failure rates usually relatively low. Their proper management, however, is imperative since the investment cost in spares is considerable. Since Sherbrooke [1] developed a multi-echelon technique for repairable item control (METRIC), the spares provisioning problem for repairable items has been of interests. Muckstadt [2], Graves [3] and Schultz [4] are only a few among the researchers that have investigated the spares provisioning problem for repairable items (see [5 and 6] for surveys). Most of the previous studies are based on the assumption that repair capacity is infinite. Lau et al. [7] considered a multi-echelon repairable item inventory system with infinite repair capacity and passivation. De Smidt-Destombes et al. [8] presented a spare parts model with cold-standby redundancy on system level under the
Corresponding author. Tel.: + 82 31 460 5545.
E-mail addresses:
[email protected] (C.-W. Park),
[email protected] (H.-S. Lee). 0305-0548/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2011.01.020
assumption that a failure is immediately restored. Kranenburg and Van Houtum [9] studied a multi-item, single-stage spare parts inventory model with multiple customer classes. Hill et al. [10] studied a single-item, two-echelon inventory with lost sales and a batch ordering policy. Al-Rifaia and Rossettib [11] and Topan et al. [12] considered two-echelon inventory systems in which the central warehouse operates under a batch ordering policy and the local warehouses implement basestock policy. Due to substantial computing time required for exact analysis, most research has focused on developing approximation methods using queueing theory or Markov models. For finite repair capacity models, Gross et al. [13] and Madu [14] proposed mathematical models using a closed queueing network which consists of one base and a two-stage repair depot. They proposed a heuristic method to decide the number of repairable items and the capacity of repair depot. Albright and Gupta [15] analyzed a multi-echelon, multi-indentured model that consists of a finite number of identical repairmen. Sleptchenko et al. [16] analyzed the trade-off between inventory levels and repair capacity in a multi-echelon, multi-indentured model. Dı´az and Fu [17], Spanjers et al. [18], Jung et al. [19] and Kim et al. [20] developed an approximation method for a multi-item, multi-echelon model with finite repair capacity. Rappold and Van Roo [21] presented an approach to model and solve the joint problem of facility location, inventory allocation, and capacity investment when demand is stochastic. De Smidt-Destombes
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
et al. [22,23] considered k-out-of-N systems under block replacement sharing limited spares and repair capacity. In the above-mentioned studies, it is assumed that repair of the failed items does not require parts. Therefore, when an item fails, it is directly sent to a repair depot. While this is a realistic assumption in some cases, in many real cases, parts are needed to repair failed items. Abboud and Daigle [24] introduced a spares provisioning problem incorporating a parts inventory system. They assumed that the system consists of a base, a parts inventory system and a repair depot. When an item in the base fails, the failed item is supplied with a part at the parts inventory system and sent to the repair depot. They assumed that (S 1,S) inventory policy is used at the parts inventory system. For this system, they developed a method using Little’s result for the case that the value of S is zero or infinite. In this study, the model of Abboud and Daigle is extended by relaxing the following assumptions: First, while there is only one base in the model of Abboud and Daigle, our model assumes that there are multiple bases and the mean time before failure for items owned by one base differs from that for items owned by another base. Second, in our model we use a continuous review (S,Q) inventory policy. In the (S,Q) inventory policy, each time the accumulated demand for parts reaches level Q, a batch order of size Q is placed so that sum of the on-hand inventory and the onorder amount becomes a target number, S. The (S1,S) policy is a special case of the (S,Q) inventory policy with Q¼1. The advantage of the (S,Q) inventory policy over the (S 1,S) inventory policy is that the ordering cost can be significantly reduced by allowing batch orders. However, due to the characteristic of batch orders, the analysis becomes much more complicated. The purpose of this paper is to develop a new approximate algorithm for the analysis of the proposed system based on a product-form approximation method [25]. In Section 2, the proposed system is defined and is modelled as a multi-class closed queueing network. The product-form approximation method is adapted to allow for the efficient analysis of the multi-class closed queueing network. Sections 3 and 4 explain the algorithm in
1585
detail. Performance measures are calculated in Section 5. In Section 6, the results obtained from the algorithm are compared with those obtained by simulation and the cost effectiveness of the (S,Q) inventory policy is evaluated with numerical results. Finally, Section 7 concludes the paper.
2. Model description and analysis overview 2.1. Model description A maintenance network consists of bases, a parts inventory system, and a repair depot (see Fig. 1 for a pictorial illustration). It is assumed that the total number of bases is R and the total number of items owned by base r is Nr, r ¼ 1, . . . R. At base r there can be as many as Dr items in operation at any time (Dr rNr). Since items in operation are subject to failure, each base maintains (Nr Dr) spare items to increase system availability. Operating times of the item before failure at base r follow an exponential distribution with rate lr. When an item fails, it always requires a part, which must be obtained from the parts inventory system. If a part is available at the parts inventory system at the instant an item fails, the failed item is immediately sent to the repair depot. On the other hand, if a part is not available when an item fails, the failed item should wait at the parts inventory system until a part becomes available before it enters the repair depot. Therefore, at the parts inventory system synchronization occurs between failed items and parts. In the repair depot, M identical repair channels are retained to service failed items. Repair times are assumed to follow an exponential distribution with rate m. Once the repair of an item owned by base r is completed, it is sent to base r immediately. At the parts inventory system, each time the accumulated demands for parts reach Q, an order of size Q is placed so that the sum of the on-hand inventory and the on-order amount is equal to S. Order lead times are assumed to follow an exponential distribution with rate t.
Fig. 1. The multi-class maintenance network with a part inventory system.
1586
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
In this study, items circulate permanently within the network in two alternating modes: operative mode and failure mode. Items in the failure mode should be repaired before they go into operative mode. Hence, if we consider items owned by each base as each separate class of customers and the parts as external resources, the maintenance network can be modelled as a multi-class closed queueing network. If base r is expressed as node r, r ¼ 1, . . . R, the parts inventory system as node (R+1), and the repair depot as node (R+2), then the network consists of (R+2) nodes. Since node (R+1) is a synchronization station, the multi-class closed queueing network is not a separable network [26]. Therefore, to analyze the network, we have to rely on an approximation method.
2.2. Analysis overview To analyze the multi-class closed queueing network, we will devise a product-form approximation method, based on Marie’s method [27]. The product-form approximation method has mainly been developed in the context of a single-job class closed queueing network, where the equilibrium distribution of jobs in the network does not have a product-form solution. The principle of the method is to approximate the performance of the original network by that of a product-form queueing network with loaddependent service rates. Details of the product-form approximation method can be found in Baynat and Dallery [25]. Very few results have been obtained on the extension of product-form approximation methods to multi-class networks. A general extension of product-form approximation methods was attempted by Neuse and Chandy [28]. Their method is an extension of Marie’s method, based on Norton’s aggregation theorem. Baynat and Dallery [29] reported that the method is computationally very demanding because it involves repetitive solution of a closed multi-class queueing network. Baynat and Dallery proposed an alternative method of extending Marie’s method to multi-class networks. The idea is to associate with the original network, R equivalent single-class product-form networks, each of them corresponding to a particular class of customers. Then, the method estimates the parameters of these R equivalent networks. Although their method is fairly general and has definite advantage over that by Neurse and Chandy, it needs to solve R equivalent single-class product-from networks and bear heavy computational cost in the analysis of the sub-networks in isolation. In this study, Marie’s method is adapted to make the analysis especially suitable for our multi-class closed queueing maintenance network. Note that each item owned by base r can visit only three nodes among the (R+ 2) nodes in the maintenance network: base r, the parts inventory system and the repair depot. The basic idea is to simplify the estimation of the parameters of the multi-class network using the routing information of the maintenance network. The main advantage of this approach is that the multi-class network can be treated as if it were a singleclass network, thus reducing the computational complexity significantly. Our general approach is presented below. The algorithm is presented in detail in Section 3. In Marie’s method, the closed queueing network is first partitioned into a set of sub-networks. Each sub-network is then approximated to an exponential station with state-dependent service rates. The resulting network, which is called an equivalent network, is then a product-form network which can be solved quickly using any computational algorithm if the state-dependent services are known. To obtain the state-dependent service rates, each sub-network is analyzed in isolation as an open queueing system fed by a state-dependent Poisson arrival process. According to Baynat and Dallery [25], for the product-form approximation method to be accurate, certain conditions must be
satisfied. One of these conditions is that the original network should be partitioned so that transitions between two sub-networks involve only a single customer. In our system, owing to the (S,Q) policy at the parts inventory system, more than one item can be released from the parts inventory system to the repair depot. This implies that these two stations should be in the same sub-network. In this paper, the network is partitioned into (R+1) sub-networks. Sub-network R0 consists of the parts inventory system and the repair depot and sub-network Rr, r ¼ 1, . . . ,R, consists of single station, which is base r in the original network. We denote station i in the product-form network corresponding to sub-network Ri by Si, and the state-dependent service rate of Si per unit time by mi(n). Let pi,j(r) denote the probability that a class r customer that departs from station Si will proceed next to station Sj. Then the maintenance network can be approximated to a BCMP network with R distinct customer classes and (R+1) stations, as shown in Fig. 2 (see [26] for further details of the BCMP theorem). Let the BCMP network be in state n ¼ ðn 0 , . . . ,n i , . . . ,n R Þ, where n i ¼ ðni1 ,ni2 , . . . ,niR Þ, implying that at station Si there are nir P customers of class r for 1 rrrR. Note that Nr ¼ Ri¼ 0 nir for r ¼ 1,2, . . . R. If we let N represent the total number of customers P in the network, then N ¼ Rr ¼ 1 Nr . Now the set of all feasible states can be explicitly defined by ( u SðNR ,R þ1Þ ¼ n 0 r nir r Nr ,0 ri r R,1 r r rR and ) R X nir ¼ Nr ,1r r rR : i¼0
NRu
¼ ðN1 ,N2 , . . . ,NR Þ. where Let ni denote the total number of customers at station Si, i.e., P ni ¼ Rr ¼ 1 nir . Then, according to the BCMP theorem, the steadystate probability of state n ¼ ðn 0 , . . . ,n R Þ is given by Pðn Þ ¼
R 1 Y f ðn Þ, u GðNR Þ i ¼ 0 i i
ð1Þ
where R Y 1 nir e : n ! ir m ðkÞ k¼1 i r ¼ 1 ir
fi ðn i Þ ¼ Qni
ni !
ð2Þ
In Eq. (1), GðNRu Þ is a normalization constant chosen so that all Pðn Þ sum to one as follows: GðNRu Þ ¼
X
R Y
fi ðn i Þ:
n A SðN u R ,R þ 1Þ i ¼ 0
Fig. 2. The BCMP network with distinct customer classes and stations.
ð3Þ
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
In Eq. (2), eir, the relative frequency of visits to station Si by class r customers, is the solutions of R sets of linear equations eir ¼
R X
ejr pji ðrÞ,
0 ri r R,
1r r r R:
ð4Þ
j¼0
Using the routing information of the network as shown in Fig. 2, Eq. (4) of station Si except for station S0 can be transformed as follows (1ri, r rR): ( e0r , r ¼ i, ð5Þ eir ¼ 0, r ai: In this study, e0r is set to 1 for r ¼ 1,2, . . . R. From Eqs. (2) and (5), fi ðn i Þ, 1rirR, is simply expressed as follows: ( hi ðnii Þ, ni1 ¼ . . . ¼ niði1Þ ¼ niði þ 1Þ ¼ niR ¼ 0, fi ðn i Þ ¼ ð6Þ 0, o=w: where hi ðnÞ ¼
8 > < 1,
n ¼ 0,
> : Qn
enii
k¼1
mi ðkÞ
,
o=w:
From Eq. (6), it is clear that if the state-dependent service rates of each station are known, every station except for station S0 can be treated as if it is Type 1 (FCFS) node in a single-class separable network [26]. If we take advantage of this characteristic and the relation n0r ¼Nr nrr for r ¼ 1, . . . R, the customer distributions in station S0 can then be obtained. Thus, the first problem is to obtain the state-dependent service rates mi(n) for station Si in the productform network. Though the state-dependent service rates mi(n) for station Si are known, efficient algorithms for computing the normalization constant GðNuR Þ as well as the customer distributions in each station are essential. These problems are dealt with in Section 3. Let li(n) denote the state-dependent arrival rate of sub-network Ri when there are n customers in sub-network Ri, i ¼ 0, . . . ,R. State-dependent service rate mi(n) then can be obtained by analyzing sub-network Ri. Sub-network Ri can be analyzed as an open queueing system with arrival rate li(n) when there are n items in the system. Suppose we know state-dependent arrival rates li(n) for sub-network Ri. Let us define the probability that there are n items at sub-network Ri as P i ðnÞ.
1587
Since sub-network Ri(i¼1,y,R) can be modelled as li(ni)/M/Di/Ni queue as shown in Fig. 3, the steady-state probabilities P i ðnÞ can simply be obtained. Sub-network R0 is complicated to analyze because it consists of the parts inventory system and the repair depot as shown in Fig. 4. The detailed analysis method of sub-network R0 is presented in Section 4.
3. Maintenance network analysis algorithm This section explains in detail how the product-form approximation method is applied to the maintenance network. Take the state-dependent service rates mi(n) for the product-form network as shown in Fig. 2 and define the steady-state probability that there are n customers at station Si of the product-form network as Pi(n). Then the conditional throughputs of station Si, 0 rirR, can be computed as follows: 8 > < m ðn þ 1ÞPi ðn þ 1Þ, n ¼ 0,1, . . . ,N 1, i i Pi ðnÞ ð7Þ Xi ðnÞ ¼ > : 0, n ¼ Ni : where N0 ¼N. Once the conditional throughputs Xi(n) are computed, the state-dependent arrival rates of each sub-network can be obtained by setting them equal to the conditional throughputs of the corresponding station in the product-form network
li ðnÞ ¼ Xi ðnÞ, n ¼ 0, . . . ,Ni :
ð8Þ
Thus, if we know steady-state probability Pi(n) for station Si, the state-dependent arrival rates of each sub-network can be obtained by Eqs. (7) and (8). Thus, sub-network Ri corresponding to station Si can be analyzed in isolation as an open queueing system with arrival rate li(n). The steady-state probability Pi(n), however, is not known and needs to be determined. In this study, the convolution algorithm [30] is used to compute Pi(n). From the classical formulas of the convolution algorithm and the routing information, steady-state probability Pi(n) for station Si except for S0 is expressed as follows: Pi ðnÞ ¼
hi ðnÞ G ðN nÞ, GðNRu Þ i=ðR þ 1Þ i
0 rn rNi ,
1 r ir R:
ð9Þ
In Eq. (9), Gi/(R + 1)(n) is the normalization constant for the product-form network with station i removed and population ðN1 , . . . ,Ni1 ,n,Ni þ 1 , . . . ,NR Þ. More formally X
Gi=ðR þ 1Þ ðnÞ ¼
R Y
fj ðn j Þ
j ¼ 0 n A SðNuR ,R þ 1Þ jai n i ¼ ð0,0, . . . ,Ni n, . . . ,0Þ
¼ GðN1 , . . . ,Ni1 ,n,Ni þ 1 , . . . ,NR Þ
n X
hi ðkÞGi=ðR þ 1Þ ðnkÞ,
k¼1
Fig. 3. Sub-network Rr, r ¼ 1, . . . R.
1 r n r Ni ,
Fig. 4. Sub-network R0.
1 ri rR:
ð10Þ
1588
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
Here ð0,0, . . . ,Ni n, . . . ,0Þ is (Ni n) in the ith component and (R 1) zeros. Values for normalization constant Gi/(R + 1)(n) are computed iteratively starting with initial condition Gi=ðR þ 1Þ ð0Þ ¼ GðN1 , . . . , Ni1 ,0,Ni þ 1 , . . . ,NR Þ. In Eq. (9), to compute normalization constant GðN uR Þ, let us define nuR ¼ ðn1 ,n2 , . . . ,nR Þ and the following auxiliary function: X
gi ðnuR Þ ¼
i Y
fj ðn j Þ,
0 ri r R:
ð11Þ
n A Sðnu R ,iÞ j ¼ 0
Note that GðnuR Þ is equal to gR ðnuR Þ. Now, for i40, Eq. (11) can be transformed as follows: ni X
gi ðnuR Þ ¼
0
1 r ir R:
hi ðkÞgi1 ðnuR k i Þ,
ð12Þ
¼
aX 2 ðr,nÞ k ¼ a1
1 k e0r hr ðNr kÞP0r1 ðnkÞ: k! ðr,nÞ
ð17Þ
P where a1 ðr,nÞ ¼ maxð0,n jr1 ¼ 1 Nj Þ and a2 ðr,nÞ ¼ minðn,Nr Þ. As stated above, if we obtain all steady-state probabilities Pi(n) of the BCMP network, the state-dependent arrival rates of each sub-network can be obtained by Eqs. (7) and (8), and each subnetwork can be analyzed in isolation as an open queueing system. If the steady-state probability P i ðnÞ for sub-network Ri is known, the conditional throughput of sub-network Ri when there are n items in sub-network Ri can be obtained as follows: 8 n ¼ 0, > < 0, ð18Þ ni ðnÞ ¼ l ðn1ÞP i ðn1Þ, n ¼ 1,2, . . . ,N : > i : i P i ðnÞ
k¼0 0
where k i ¼ ð0,0, . . . ,k, . . . ,0Þ (i.e. k in the ith component and (R1) zeros). Using Eq. (12), normalization constant GðnuR Þ can be iteratively computed with initial condition g0 ðnuR Þ ¼ f0 ðnuR Þ. Let G0=ðR þ 1Þ ðn 0 Þ be the normalization constant for the productform network with station 0 removed and population n 0 . From Q the definition of G0=ðR þ 1Þ ðn 0 Þ, G0=ðR þ 1Þ ðn 0 Þ is equal to Rj¼ 1 hj ðn0j Þ. Then, steady-state probability P0 ðn 0 Þ is expressed as follows: P0 ðn 0 Þ ¼
R fi ðn 0 Þ f0 ðn 0 Þ Y G0=ðR þ 1Þ ðNRu n 0 Þ ¼ h ðN n Þ: u u GðNR Þ GðNR Þ j ¼ 1 j j 0j
ð13Þ
Let n r0 ¼ ðn01 ,n02 , . . . ,n0r Þ and S0(r,n) denote the set of the feasible states for station S0 when there are R ¼r and n customers at station S0 as follows: 8 9 < = r X n0j ¼ n , S0 ðr,nÞ ¼ n r0 0 r n0j r Nj for 1 r j r r and : ; j¼1
Once the conditional throughputs have been obtained, the state-dependent service rates of the product-form network are set equal to the conditional throughputs of the corresponding subnetwork in isolation, i.e.,
mi ðnÞ ¼ ni ðnÞ, n ¼ 0, . . . ,Ni :
ð19Þ i
It should be noted that for sub-network R , i ¼ 1, . . . ,R, since the operating times of the item before failure at base i are exponential, the conditional throughput can be more easily obtained by ni(n)¼ min(Di,n)li, n ¼ 1, . . . ,Ni . Therefore, we do not have to compute mi(n) for i ¼ 1, . . . ,R. To summarize, if mi(ni) is known, li(ni) can be obtained by Eqs. (7) and (8). Conversely, if li(ni) is known, mi(ni) can be obtained by Eqs. (18) and (19). Since none of these values is known, determining mi(ni) and li(ni) is a fixed point problem and an iterative procedure is required to compute them. The iterative algorithm is as follows: 3.1. Maintenance network analysis (MNA) algorithm
r ¼ 1, . . . R, n ¼ 0, . . . ,N: Then, steady-state probability P0(n) for station S0 is defined as follows: X
P0 ðnÞ ¼
n R0 A S0 ðR,nÞ
R f0 ðn 0 Þ Y h ðN n Þ: u GðN R Þ j ¼ 1 j j 0j
ð14Þ
To compute steady-state probability P0(n), let us define the following auxiliary function: X
P0r ðnÞ ¼ n r0
r Y 1 n0j e h ðN n Þ, n ! 0j j j 0j 0 A S ðr,nÞ j ¼ 1 0j
r ¼ 1, . . . R, n ¼ 0, . . . ,N: ð15Þ
Then, from Eqs. (2) and (15), steady-state probability P0(n) in Eq. (14) is expressed as follows: P0 ðnÞ ¼
n! P R ðnÞ: GðNRu ÞPnk ¼ 1 m0 ðkÞ 0
ð16Þ
Thus, if we know the value of the auxiliary function P0R ðnÞ, steady-state probability P0(n) can be obtained by Eq. (16). For r 41, the auxiliary function P0R ðnÞ of Eq. (15) can be iteratively computed with initial condition P01 ðnÞ ¼ ½en01 h1 ðN1 nÞ=n! as follows: P0r ðnÞ ¼
aX 2 ðr,nÞ
X
k ¼ a1 ðr,nÞ n r0 A S0 ðr,nÞ and n0r ¼ k
¼
aX 2 ðr,nÞ k ¼ a1
1 k e0r hr ðNr kÞ k! ðr,nÞ
r Y 1 n0j e h ðN n Þ n ! 0j j j 0j j ¼ 1 0j
X n r1 A S0 ðr1,nkÞ 0
r1 Y 1 n0j e h ðN n Þ n ! 0j j j 0j j ¼ 1 0j
Step 0: Initialization. (a) Set mi(ni)¼ min(Di,ni)li for ni ¼ 1, . . . ,Ni , i ¼ 1, . . . ,R. (b) Set m0(n0) to some initial value for n0 ¼ 1, . . . ,N. Step 1: Iterative Procedure. Step 1.1: Analysis of the BCMP network. (a) Compute normalization constant GðNuR Þ using Eq. (12). (b) Compute Pi(ni) using Eq. (9) for ni ¼ 0, . . . ,Ni , i ¼ 1, . . . ,R. (c) Compute P0(n0) using Eq. (16) for n0 ¼ 0, . . . ,N. (d) Compute li(ni) using Eqs. (7) and (8) for ni ¼ 0, . . . ,Ni , i ¼ 0, . . . ,R. Step 1.2: Analysis of sub-network R0. (a) Compute P 0 ðn0 Þ using the methods proposed in Section 4 for n0 ¼ 0, . . . ,N. (b) Compute m0(n0) using Eqs. (18) and (19) for n0 ¼ 0, . . . ,N. Step 2: Go to Step 1 until convergence of state-dependent service rates m0(n) occurs. Step 3: Compute performance measures using Eqs. (31)–(40) in Section 5.
4. Analysis of sub-network R0 This section presents a method for analysis of sub-network R0. This section assumes that items from outside arrive according to a Poisson process with arrival rate l0(n0). Then, sub-network R0 can be modelled as a continuous time Markov chain because all events causing state transitions occur according to exponential
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
distributions. Let (i,k) represent the state that there are i parts in the parts inventory system (if negative, it indicates a state in which there are i items waiting for the parts) and k items are in the repair depot. Thus, the state of the Markov chain can be represented by a two-dimensional vector (i,k) and the state space can be defined as follows: SIR ¼ fði,kÞj N r i rS, 0 rk rN þ gi g,
OIR ði,kÞP IR ði,kÞ ¼ I1IR ði,kÞ þI2IR ði,kÞ þ I3IR ði,kÞ:
ð21Þ
The left-hand side of Eq. (21) implies the flow rate out of state (i,k). The term OIR(i,k) in Eq. (21) can be expressed as OIR ði,kÞ ¼ l0 ðjgi j þkÞ þ mI ðiÞ þ mR ðkÞ:
ii) I2IR ði,kÞ: the flow rate into state (i,k) due to the arrival of ordered parts 8 I m ðiQ ÞPIR ðiQ ,kÞ, Q r ir S, > > > > < mI ðiQ ÞP IR ðiQ ,kQ þ iÞ, 1 r i rQ 1, k ZQ i, I2IR ði,kÞ ¼ > mI ðiQ ÞPIR ðiQ ,kQ Þ, Q N ri r0, kZ Q , > > > : 0, o=w:
ð20Þ
where gi ¼min(i,0). Let mI ðiÞ ¼ bðSiÞ=Q ct represent the state-dependent orderfulfillment rate when there are i parts in the parts inventory system and mR(k)¼min(M,k)m be the state-dependent repair rate when there are k items in the repair depot. Clearly, this Markov chain is irreducible and positive recurrent. Thus, its stationary probabilities exist uniquely. Therefore, the stationary probabilities can be computed by solving balance equations. The state transition-rate diagram of the continuous time Markov chain when S ¼2, Q¼2 and N ¼4 is depicted in Fig. 5. From the state transition-rate diagram, the balance equation of state (i, k) can be obtained as follows:
ð22Þ
The right-hand side of Eq. (21) implies the flow rate into state (i,k), which is composed of the following three terms. i) I1IR ði,kÞ: the flow rate into state (i,k) due to the arrival of a failed item 8 IR 0 ri rS1, k Z1, > < l0 ðk1ÞP ði þ1,k1Þ, IR ð23Þ I1 ði,kÞ ¼ l0 ðjiþ 1j þ kÞP IR ðiþ 1,kÞ, N r i r1, > : 0, o=w:
1589
ð24Þ iii) I3IR ði,kÞ: the flow rate into state (i,k) due to the completion of repair( mR ðk þ1ÞPIR ði,k þ 1Þ, k r N1 þ gi , I3IR ði,kÞ ¼ : ð25Þ 0, o=w: Although steady-state probabilities can be obtained by solving the balance equations, the steady-state probabilities are not easy to compute if S and N are large because the number of balance equations increases very rapidly. Therefore, we need an efficient method which exploits the special structure of the Markov chain to obtain the steady-state probabilities when S and N is large. In this study, we will use a recursive technique, which was developed by Herzog et al. [31] and was used by Buzacott and Kostelski [32], among others. The recursive technique can be applied when all the states can be expressed as a linear combination of a subset of states which are called the boundary states. Thanks to the special structure of the Markov chain depicted in Fig. 5, it can be checked that if we know the steady-state probabilities of the state in TBIR ¼{(n,0)| N +1 rnrS}, we can compute the steady-state probabilities of all the other states recursively. Hence, instead of solving the balance equations associated with all states, we have only to solve the balance equations associated with (S+ N) boundary states in TBIR (|TBIR|¼ S+ N). For instance, when S¼5 and N ¼20, we have only to solve 25 balance equations associated with the states in TBIR instead of solving 336 balance equations associated with all states in SIR.
Fig. 5. The state transition-rate diagram for sub-network R0 (S¼ 2, Q ¼2, N ¼4).
1590
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
Let us express the steady-state probability of state (i, k) as a linear combination of the steady-state probabilities of the boundary states as follows: S X
P IR ði,kÞ ¼
CnIR ði,kÞP IR ðn,0Þ:
where 1R ¼ ð0,0, . . . ,1, . . . 0Þ is a vector in which there is 1 in the rth component and zero elsewhere. The mean number of items in operation at base r, D r , can be computed as
ð26Þ
n ¼ N þ 1
Dr ¼
where CnIR ði,kÞ is the coefficient of the probability PIR(n,0) for state (i,k). Therefore, if we know the values of CnIR ði,kÞ as well as the probabilities of the boundary states, the probabilities of all states can be computed using Eq. (26). To find CnIR ði,kÞ, we execute the following recursive scheme for state (n,0) for all (n,0) A TBIR.
Ni X
1r r rR:
½minðn,Dr ÞP r ðnÞ,
ð33Þ
n¼1
In the parts inventory system, the mean number of on-hand parts, I, and the mean number of items waiting for parts, N I , are calculated as follows: I¼
S N X X i P IR ði,kÞ:
ð34Þ
i¼1 k¼0
4.1. Recursive scheme of sub-network R0 IR
Step 1: Set the steady-state probability P (n,0) of boundary state (n,0) equal to 1 and the probabilities of the other boundary states equal to 0. From Eq. (26), we have the following relationship: CnIR ði,kÞ ¼ PIR ði,kÞ,
8ði,kÞ A SIR :
NI ¼
The mean number of items at the repair depot, N R , and the mean number of busy repair channels, M, can be computed as follows: NR ¼
N X
k
k¼1
M¼
N X
S X
Apply Eqs. (29) and (27) successively to find CnIR ði,k þ1Þ for 0rk rN 1, k N + 1rirS. There are (S + N) balance equations associated with the states in the set NBIR ¼{(i,k)| N + 1rirS, k¼N + gi} which are not used in calculating CnIR ði,kÞ. To determine the steady-state probabilities of the boundary states, we solve a set of simultaneous equations which consist of balance equations associated with (S +N 1) states in NBIR and a normalization equation. Once the steadystate probabilities of the boundary states are computed, all the other probabilities are computed using Eq. (26). If the steadystate probabilities of each state in sub-network R0 are computed, we can compute the probability that there are n items at subnetwork R0 as follows: P 0 ðnÞ ¼
S X
P IR ði,n þ gi Þ:
ð30Þ
i ¼ n
5. Computing performance measures Once the MNA algorithm has converged, we can compute the various performance measures of interests. Let the probability that a failed item will wait at the parts inventory system be denoted by Pb/k. Since a part is not available at the moment an item fails with probability Pb/k, Pb/k can be computed as Pb=k ¼ 1
S X N X
P IR ði,kÞ:
ð31Þ
i¼1k¼0
Using the normalization constant, the throughput of class r customers at station Si is calculated as Xir ¼ eir GðNRu 1 R Þ=GðNRu Þ,
0 ri r R,1 r r rR:
ð32Þ
P IR ði,kÞ:
ð36Þ
i ¼ N þ k
minðk,MÞ
k¼1
Apply Eqs. (28) and (27) successively to find CnIR ðN,0Þ. Step 3: For 0 rkrN 1, k N + 1rirS, the balance equation associated with state (i,k) can be transformed into the following equation: PIR ði,k þ1Þ ¼ OIR ði,kÞP IR ði,kÞI1IR ði,kÞI2IR ði,kÞ =mR ðk þ1Þ: ð29Þ
ð35Þ
i¼1 k¼0
ð27Þ
Step 2: The balance equation associated with state ( N,0) can be transformed into the following equation: PIR ðN,0Þ ¼ I1IR ðN,0Þ þ I2IR ðN,0Þ þI3IR ðN,0Þ =OIR ðN,0Þ: ð28Þ
N N i X X i P IR ði,kÞ:
S X
PIR ði,kÞ:
ð37Þ
i ¼ N þ k
Using Little’s law, M can be computed as X0/m where P X0 ¼ Rr ¼ 1 X0r . If we choose an item randomly at sub-network R0, the probability that this item belongs to base r is p0 ðrÞ ¼ X0r =X0 ,
1 r r rR:
ð38Þ
Using Eq. (38), the mean number of items owned by base r at parts inventory system, N I ðrÞ, and the mean number of items owned by base r at repair depot, N R ðrÞ, can be computed as follows: N I ðrÞ ¼ p0 ðrÞN I , N R ðrÞ ¼ p0 ðrÞN R ,
1 rr r R: 1 r r rR:
ð39Þ ð40Þ
6. Numerical results This section reports the numerical results of the method proposed in this paper. The algorithm developed in this study was implemented on an IBM PC and tested using various scenario sets. The results of the approximation method were compared with those of simulation obtained using the SIMLIB programme [33] on the same PC. For the simulation of each problem, 10 replications were made and each replication was run for 105 simulation times. The convergence criterion was that the values of m0(n) in successive iterations of the algorithm should be within 10 5. Four scenario sets, (A, B, C, D), are presented in this section. In scenario set A, the number of bases is changed. In scenario sets B and C, the effects of changing the values of S and Q in the (S,Q) inventory policy are investigated, respectively. In scenario set D, the number of repair channels at the repair depot is changed. The experimental data are presented in Table 1. The experimental results are summarized in Tables 2–5. In Tables 2–5, we also provide the absolute relative errors and 95% confidence intervals with respect to the simulation results. As shown in Tables 2–5, the proposed approximation method, in general, yields fairly accurate results. Most of the estimated values have a relative error of less than 0.3% and are within the
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
1591
Table 1 Experiment data for 20 scenarios. Scenario set
A B C D
R
Base i (1r irR)
1,2,3,4,5 3 3 4
Parts inventory system
Repair center
Ni
Di
li
S
Q
t
M
m
5 12 2i 12 2i 9i
3 10 2i 10 2i 8i
2 2 +0.2i 2 +0.2i 2 +0.1i
3 4,5,6,7,8 8 10
2 3 4,5,6,7,8 4
3 2 2 2
3 7 7 6,7,8,9,10
4 2 2 2
Table 2 Results of scenario set A. Performance measures
Simulation ( 795% confidence interval)
Approximation (relative error %) R ¼1
R¼ 2
R¼ 3
R ¼4
R¼ 5
R ¼1
R ¼2
R ¼3
R ¼4
R ¼5
I
0.4076 (0.0980%) 1.0711
0.6849 (0.0438%) 0.5152
0.7694 (0.0000%) 0.3673
0.7813 (0.0128%) 0.3476
0.7818 (0.0384%) 0.3467
0.4080 ( 7 0.0005) 1.0701
0.6852 (7 0.0002) 0.5147
0.7694 (7 0.0007) 0.3672
0.7814 (7 0.0005) 0.3473
0.7815 (7 0.0003) 0.3473
NI
(0.0934%) 0.3317
(0.0971%) 1.1378
(0.0272%) 1.6983
(0.0864%) 1.8360
(0.1728%) 1.8464
( 7 0.0012) 0.3323
(7 0.0004) 1.1383
(7 0.0012) 1.6953
(7 0.0009) 1.8371
(7 0.0006) 1.8446
NR
(0.1806%) 1.3679
(0.0439%) 3.3526
(0.1770%) 7.0010
(0.0599%) 11.9147
(0.0976%) 17.0228
( 7 0.0011) 1.3675
(7 0.0013) 3.3544
(7 0.0030) 7.0014
(7 0.0029) 11.9170
(7 0.0029) 17.0241
M
(0.0293%) 1.3208
(0.0537%) 2.3419
(0.0057%) 2.8732
(0.0193%) 2.9913
(0.0076%) 2.9998
( 7 0.0018) 1.3203
(7 0.0042) 2.3430
(7 0.0088) 2.8734
(7 0.0068) 2.9912
(7 0.0071) 2.9998
(0.0379%) 5.2830 (0.0114%)
(0.0469%) 9.3678 (0.0032%)
(0.0070%) 11.4929 (0.0217%)
(0.0033%) 11.9650 (0.0217%)
(0.0000%) 11.9991 (0.0333%)
( 7 0.0016) 5.2836 ( 7 0.0028)
(7 0.0016) 9.3675 (7 0.0030)
(7 0.0009) 11.4954 (7 0.0065)
(7 0.0001) 11.9624 (7 0.0059)
(7 0.0000) 12.0031 (7 0.0080)
2.6415 (0.0038%)
2.3419 (0.0299%)
1.9155 (0.0052%)
1.4956 (0.1138%)
1.1999 (0.0083%)
2.6416 ( 7 0.0011)
2.3412 (7 0.0013)
1.9156 (7 0.0029)
1.4939 (7 0.0020)
1.2000 (7 0.0013)
Pb/k
R P
X0r
r¼1
DR
Table 3 Results of scenario set B. Performance measures
Simulation ( 7 95% confidence interval)
Approximation (relative error %) S ¼4
S ¼5
S¼6
S¼ 7
S ¼8
S¼ 4
S¼ 5
S ¼6
S ¼7
S ¼8
I
0.8911 (0.0000%) 0.1983
0.8191 (0.0000%) 0.3782
0.7300 (0.0274%) 0.6470
0.6291 (0.0318%) 1.0165
0.5234 (0.1336%) 1.4918
0.8911 ( 70.0006) 0.1981
0.8191 (70.0004) 0.3779
0.7298 (7 0.0006) 0.6472
0.6289 (7 0.0007) 1.0178
0.5241 (7 0.0006) 1.4894
N I ð1Þ
(0.1010%) 1.6824
(0.0794%) 1.3587
(0.0309%) 1.0677
(0.1277%) 0.8154
(0.1611%) 0.6044
( 70.0017) 1.6831
(70.0012) 1.3589
(7 0.0022) 1.0676
(7 0.0029) 0.8161
(7 0.0027) 0.6063
N I ð2Þ
(0.0416%) 1.3764
(0.0147%) 1.1117
(0.0094%) 0.8737
(0.0858%) 0.6672
(0.3134%) 0.4946
( 70.0028) 1.3773
(70.0023) 1.1116
(7 0.0015) 0.8741
(7 0.0022) 0.6671
(7 0.0024) 0.4959
N I ð3Þ
(0.0653%) 1.0515
(0.0090%) 0.8493
(0.0458%) 0.6675
(0.0150%) 0.5098
(0.2621%) 0.3779
( 70.0023) 1.0517
(70.0022) 0.8490
(7 0.0011) 0.6673
(7 0.0019) 0.5095
(7 0.0019) 0.3786
N R ð1Þ
(0.0190%) 5.7394
(0.0353%) 6.0522
(0.0300%) 6.3354
(0.0589%) 6.5826
(0.1849%) 6.7901
( 70.0011) 5.7441
(70.0013) 6.0580
(7 0.0010) 6.3411
(7 0.0016) 6.5903
(7 0.0011) 6.7941
N R ð2Þ
(0.0818%) 4.6954
(0.0957%) 4.9519
(0.0899%) 5.1842
(0.1168%) 5.3867
(0.0589%) 5.5568
( 70.0038) 4.6938
(70.0032) 4.9496
(7 0.0039) 5.1825
(7 0.0038) 5.3858
(7 0.0035) 5.5534
N R ð3Þ
(0.0341%) 3.5871
(0.0465%) 3.7833
(0.0328%) 3.9609
(0.0167%) 4.1157
(0.0612%) 4.2458
( 70.0033) 3.5832
(70.0032) 3.7794
(7 0.0028) 3.9553
(7 0.0032) 4.1114
(7 0.0031) 4.2393
M
(0.1088%) 6.9121
(0.1032%) 6.9415
(0.1416%) 6.9620
(0.1046%) 6.9760
(0.1533%) 6.9851
( 70.0024) 6.9122
(70.0015) 6.9418
(7 0.0026) 6.9619
(7 0.0029) 6.9761
(7 0.0021) 6.9851
D1
(0.0014%) 5.6584 (0.0389%) 4.6292 (0.0173%) 3.5365 (0.0198%) 2.5720
(0.0043%) 5.6820 (0.0123%) 4.6491 (0.0430%) 3.5519 (0.0028%) 2.5827
(0.0014%) 5.6985 (0.0176%) 4.6629 (0.0172%) 3.5626 (0.0253%) 2.5902
(0.0014%) 5.7096 (0.0228%) 4.6723 (0.0343%) 3.5699 (0.0364%) 2.5953
(0.0000%) 5.7170 (0.0332%) 4.6785 (0.0427%) 3.5747 (0.0419%) 2.5986
( 70.0008) 5.6562 ( 70.0034) 4.6300 ( 70.0023) 3.5372 ( 70.0018) 2.5727
(70.0006) 5.6827 (70.0019) 4.6511 (70.0028) 3.5518 (70.0013) 2.5829
(7 0.0005) 5.6975 (7 0.0062) 4.6637 (7 0.0043) 3.5617 (7 0.0032) 2.5911
(7 0.0004) 5.7083 (7 0.0035) 4.6707 (7 0.0025) 3.5686 (7 0.0020) 2.5935
(7 0.0003) 5.7189 (7 0.0032) 4.6805 (7 0.0030) 3.5762 (7 0.0011) 2.5994
D2
(0.0272%) 1.9288
(0.0077%) 1.9371
(0.0347%) 1.9429
(0.0694%) 1.9468
(0.0308%) 1.9494
( 70.0023) 1.9281
(70.0020) 1.9380
(7 0.0027) 1.9426
(7 0.0022) 1.9461
(7 0.0026) 1.9498
D3
(0.0363%) 1.3602
(0.0464%) 1.3661
(0.0154%) 1.3702
(0.0360%) 1.3730
(0.0205%) 1.3749
( 70.0020) 1.3594
(70.0018) 1.3657
(7 0.0024) 1.3712
(7 0.0021) 1.3729
(7 0.0015) 1.3757
(0.0588%)
(0.0293%)
(0.0729%)
(0.0073%)
(0.0582%)
( 70.0017)
(70.0013)
(7 0.0018)
(7 0.0021)
(7 0.0014)
Pb/k
X01 X02 X03
1592
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
Table 4 Results of scenario set C. Performance measures
Approximation (relative error %)
Simulation (7 95% confidence interval)
Q¼ 4
Q¼ 5
Q¼ 6
Q ¼7
Q ¼8
Q¼ 4
Q ¼5
Q¼ 6
Q ¼7
Q¼ 8
I
0.5610 (0.0178%) 1.4223
0.5953 (0.0841%) 1.3535
0.6207 (0.0644%) 1.2903
0.6419 (0.0624%) 1.2257
0.6623 (0.0756%) 1.1584
0.5609 ( 70.0009) 1.4231
0.5948 (7 0.0009) 1.3553
0.6211 (7 0.0010) 1.2884
0.6415 (7 0.0009) 1.2273
0.6618 (7 0.0007) 1.1609
N I ð1Þ
(0.0562%) 0.7737
(0.1328%) 0.9398
(0.1475%) 1.1048
(0.1304%) 1.2656
(0.2154%) 1.4220
( 70.0037) 0.7739
(7 0.0040) 0.9375
(7 0.0044) 1.1061
(7 0.0039) 1.2649
(7 0.0029) 1.4203
N I ð2Þ
(0.0258%) 0.6331
(0.2453%) 0.7690
(0.1175%) 0.9039
(0.0553%) 1.0353
(0.1197%) 1.1630
( 70.0028) 0.6329
(7 0.0021) 0.7667
(7 0.0042) 0.9041
(7 0.0031) 1.0346
(7 0.0043) 1.1604
N I ð3Þ
(0.0316%) 0.4837
(0.3000%) 0.5875
(0.0221%) 0.6905
(0.0677%) 0.7909
(0.2241%) 0.8883
( 70.0021) 0.4841
(7 0.0021) 0.5861
(7 0.0033) 0.6908
(7 0.0029) 0.7898
(7 0.0033) 0.8851
N R ð1Þ
(0.0826%) 6.6271
(0.2389%) 6.4704
(0.0434%) 6.3180
(0.1393%) 6.1728
(0.3615%) 6.0351
( 70.0016) 6.6328
(7 0.0016) 6.4796
(7 0.0023) 6.3226
(7 0.0027) 6.1814
(7 0.0026) 6.0439
N R ð2Þ
(0.0859%) 5.4230
(0.1420%) 5.2942
(0.0728%) 5.1689
(0.1391%) 5.0495
(0.1456%) 4.9360
( 70.0049) 5.4209
(7 0.0034) 5.2955
(7 0.0045) 5.1664
(7 0.0039) 5.0504
(7 0.0040) 4.9381
N R ð3Þ
(0.0387%) 4.1434
(0.0245%) 4.0448
(0.0484%) 3.9488
(0.0178%) 3.8571
(0.0425%) 3.7701
( 70.0027) 4.1388
(7 0.0030) 4.0414
(7 0.0047) 3.9439
(7 0.0034) 3.8540
(7 0.0038) 3.7695
M
(0.1111%) 6.9682
(0.0841%) 6.9429
(0.1242%) 6.9089
(0.0804%) 6.8661
(0.0159%) 6.8153
( 70.0022) 6.9681
(7 0.0013) 6.9430
(7 0.0034) 6.9091
(7 0.0032) 6.8667
(7 0.0028) 6.8169
D1
(0.0014%) 5.7034 (0.0105%) 4.6671 (0.0193%) 3.5659 (0.0084%) 2.5925
(0.0014%) 5.6831 (0.0176%) 4.6501 (0.0215%) 3.5526 0.0479 2.5832
(0.0029%) 5.6557 (0.0283%) 4.6271 (0.0022%) 3.5349 (0.0339%) 2.5708
(0.0087%) 5.6214 (0.0071%) 4.5984 (0.0239%) 3.5126 (0.0456%) 2.5552
(0.0235%) 5.5804 (0.0197%) 4.5642 (0.0044%) 3.4861 (0.0258%) 2.5366
( 70.0003) 5.7028 ( 70.0029) 4.6662 ( 70.0028) 3.5662 ( 70.0029) 2.5931
(7 0.0007) 5.6821 (7 0.0041) 4.6491 (7 0.0039) 3.5509 (7 0.0028) 2.5827
(7 0.0010) 5.6573 (7 0.0041) 4.6272 (7 0.0035) 3.5361 (7 0.0030) 2.5712
(7 0.0012) 5.6210 (7 0.0030) 4.5995 (7 0.0033) 3.5110 (7 0.0024) 2.5535
(7 0.0015) 5.5793 (7 0.0027) 4.5644 (7 0.0026) 3.4852 (7 0.0031) 2.5357
D2
(0.0231%) 1.9446
(0.0194%) 1.9375
(0.0156%) 1.9280
(0.0666%) 1.9160
(0.0355%) 1.9017
( 70.0040) 1.9453
(7 0.0033) 1.9370
(7 0.0017) 1.9286
(7 0.0024) 1.9142
(7 0.0028) 1.9007
D3
(0.0360%) 1.3715
(0.0258%) 1.3664
(0.0311%) 1.3596
(0.0940%) 1.3510
(0.0526%) 1.3408
( 70.0016) 1.3708
(7 0.0024) 1.3664
(7 0.0025) 1.3595
(7 0.0017) 1.3505
(7 0.0015) 1.3399
(0.0511%)
(0.0000%)
(0.0074%)
(0.0370%)
(0.0672%)
( 70.0016)
(7 0.0010)
(7 0.0017)
(7 0.0016)
(7 0.0010)
Pb/k
X01 X02 X03
95% confidence intervals of the simulation results. The maximum relative error observed in the four scenario sets is very small (0.36%). As in the previous algorithms based on Mare’s method, we could not prove the convergence of our algorithm. However, all of the problems tested to date have always converged. Typically, the number of iterations required in the MNA algorithm is less than 5. Although the results of only four scenario sets are presented here, the other problems that were tested showed similar results. As noted above, for a given configuration of the system (values of R, Ni,Di, li, t, M and m are given), if the (S,Q) inventory policy is specified, we can evaluate the performances of the system. It is then interesting to use the proposed approximation method to obtain the optimal (S,Q) inventory policy for a given configuration of the system. To this end, the following cost function to operate the (S,Q) inventory policy has been chosen CðS,Q Þ ¼
R X
½CD ðrÞðDr D r Þ þ CH I þ CQ ðX0 =Q Þ þ CP X0
ð41Þ
r¼1
where C(S,Q) is the expected cost per unit time with control values S and Q, CD(r) is the penalty cost per unit time per each shortage of item at base r when the number of items in operation at base r is less than Dr, r ¼1,y,R, CH is the unit holding cost of parts per unit time at the parts inventory system, CQ is the ordering cost of parts at the parts inventory system, and CP is the unit purchase cost of parts at the parts inventory system. We have tested numerous problem sets of the (S,Q) inventory policy for the cost function detailed above. It is observed from extensive numerical tests that the cost function C(S,Q) is convex with respect to Q for a given value of S. Let us denote the optimal
Q value which minimizes C(S,Q) for a given S ¼s by Qn(s). It is also observed that the cost function C(s,Qn(s)) is unimodal with respect to s. In Table 6, an example is presented that shows the behaviour of the cost function C(S,Q). From Table 6, we can check the unimodality property of the cost function C(s,Qn(s)) as well as the convexity property of the cost function C(S,Q) with respect to Q for a given value of S. Although we could not prove these properties, in all the problems we have tested, both convexity and unimodality properties are observed to hold. If we assume that these two properties of the cost function hold, we can obtain the optimal (S,Q) inventory policy which minimizes the expected cost per unit time using the following simple procedure. For a given value of S, we search for Qn(S) in the range of 1rQrN + S. We start with S¼0 and find Qn(0) to obtain C(0,Qn(0)). To find Qn(0), we simply compute C(0,k) from k ¼1 up to the point k¼Qn(0) + 1, at which C(0,k) is first increased. Then, due to the convexity of C(0,k), the minimum value of C(0,k) is found at C(0,Qn(0)). Once C(0,Qn(0)) is obtained, we find C(S,Qn(S)) sequentially in the order of S ¼ 1,2,3 . . .. Our procedure repeats this process (increase S by 1 and compute C(S,Qn(S)) for a new value of S) until C(S,Qn(S)) is first increased. Then, by the unimodality property, the local minimum point encountered is a global minimum and the optimal control values Sn and Qn are obtained. In order to verify the properties of the cost function C(S, Q), we have made extensive numerical tests. Among these, we present three problem sets, (E, F, G). For these problem sets, using a procedure developed in Section 6, we have done a simple sensitivity analysis to measure the influence of cost parameters in Eq. (41) on the optimal (S,Q) inventory policy. In problem set E, the ordering cost, CQ, is changed. In problem set F, the effects of changing the holding cost, CH, is investigated. In problem set G, the penalty cost due to shortage of operative items at each base,
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
1593
Table 5 Results of scenario set D. Performance measures
Simulation ( 795% confidence interval)
Approximation (relative error %) M ¼6
M¼7
M¼8
M¼ 9
M ¼10
M ¼6
M ¼7
M¼ 8
M¼9
M ¼10
I
0.2858 (0.1049%) 3.1709
0.3797 (0.0791%) 2.5700
0.4697 (0.0213%) 2.0685
0.5502 (0.0364%) 1.6642
0.6175 (0.0324%) 1.3526
0.2861 ( 7 0.0009) 3.1699
0.3794 (7 0.0011) 2.5723
0.4698 (70.0008) 2.0688
0.5500 (7 0.0005) 1.6651
0.6173 (7 0.0006) 1.3531
N I ð1Þ
(0.0315%) 0.2036
(0.0894%) 0.3224
(0.0145%) 0.4643
(0.0541%) 0.6169
(0.0370%) 0.7635
( 7 0.0067) 0.2041
(7 0.0063) 0.3222
(70.0047) 0.4643
(7 0.0024) 0.6169
(7 0.0039) 0.7631
N I ð2Þ
(0.2450%) 0.1800
(0.0621%) 0.2855
(0.0000%) 0.4118
(0.0000%) 0.5480
(0.0524%) 0.6791
( 7 0.0009) 0.1803
(7 0.0017) 0.2853
(70.0014) 0.4120
(7 0.0017) 0.5476
(7 0.0015) 0.6784
N I ð3Þ
(0.1664%) 0.1557
(0.0701%) 0.2475
(0.0485%) 0.3574
(0.0730%) 0.4763
(0.1032%) 0.5908
( 7 0.0007) 0.1559
(7 0.0017) 0.2475
(70.0011) 0.3575
(7 0.0014) 0.4764
(7 0.0013) 0.5900
N I ð4Þ
(0.1283%) 0.1309
(0.0000%) 0.2083
(0.0280%) 0.3013
(0.0210%) 0.4019
(0.1356%) 0.4990
( 7 0.0007) 0.1311
(7 0.0013) 0.2082
(70.0011) 0.3011
(7 0.0013) 0.4015
(7 0.0015) 0.4977
N R ð1Þ
(0.1526%) 6.0570
(0.0480%) 5.6541
(0.0664%) 5.2356
(0.0996%) 4.8245
(0.2612%) 4.4518
( 7 0.0006) 6.0593
(7 0.0010) 5.6571
(70.0011) 5.2397
(7 0.0010) 4.8295
(7 0.0013) 4.4585
N R ð2Þ
(0.0380%) 5.3547
(0.0530%) 5.0069
(0.0782%) 4.6438
(0.1035%) 4.2855
(0.1503%) 3.9594
( 7 0.0041) 5.3545
(7 0.0029) 5.0082
(70.0023) 4.6434
(7 0.0034) 4.2867
(7 0.0027) 3.9615
N R ð3Þ
(0.0037%) 4.6335
(0.0260%) 4.3392
(0.0086%) 4.0305
(0.0280%) 3.7245
(0.0530%) 3.4450
( 7 0.0018) 4.6326
(7 0.0030) 4.3386
(70.0019) 4.0280
(7 0.0026) 3.7222
(7 0.0029) 3.4453
N R ð4Þ
(0.0194%) 3.8950
(0.0138%) 3.6526
(0.0621%) 3.3970
(0.0618%) 3.1426
(0.0087%) 2.9094
( 7 0.0021) 3.8907
(7 0.0027) 3.6499
(70.0023) 3.3930
(7 0.0026) 3.1388
(7 0.0030) 2.9066
M
(0.1105%) 5.9993
(0.0740%) 6.9938
(0.1179%) 7.9663
(0.1211%) 8.8788
(0.0963%) 9.6798
( 7 0.0022) 5.9993
(7 0.0022) 6.9937
(70.0020) 7.9665
(7 0.0020) 8.8785
(7 0.0016) 9.6812
D1
(0.0000%) 3.6447 (0.0027%) 3.2221 (0.0124%) 2.7881 (0.0179%) 2.3437 (0.0043%) 1.7356
(0.0014%) 4.2399 (0.0496%) 3.7546 (0.0160%) 3.2539 (0.0031%) 2.7390 (0.0073%) 2.0190
(0.0025%) 4.8199 (0.0166%) 4.2751 (0.0094%) 3.7104 (0.0081%) 3.1273 (0.0352%) 2.2952
(0.0034%) 5.3622 (0.0354%) 4.7631 (0.0168%) 4.1396 (0.0386%) 3.4928 (0.0115%) 2.5534
(0.0145%) 5.8369 (0.0086%) 5.1913 (0.0231%) 4.5168 (0.0465%) 3.8146 (0.0420%) 2.7795
( 7 0.0000) 3.6448 ( 7 0.0028) 3.2225 ( 7 0.0026) 2.7886 ( 7 0.0023) 2.3436 ( 7 0.0023) 1.7366
(7 0.0002) 4.2378 (7 0.0026) 3.7552 (7 0.0031) 3.2540 (7 0.0024) 2.7388 (7 0.0014) 2.0207
(70.0006) 4.8191 (70.0030) 4.2747 (70.0026) 3.7107 (70.0019) 3.1284 (70.0023) 2.2958
(7 0.0010) 5.3603 (7 0.0021) 4.7639 (7 0.0026) 4.1412 (7 0.0037) 3.4924 (7 0.0019) 2.5532
(7 0.0014) 5.8364 (7 0.0030) 5.1901 (7 0.0039) 4.5147 (7 0.0028) 3.8130 (7 0.0024) 2.7778
D2
(0.0576%) 1.4646
(0.0841%) 1.7066
(0.0261%) 1.9432
(0.0078%) 2.1650
(0.0612%) 2.3597
( 7 0.0033) 1.4652
(7 0.0021) 1.7064
(70.0020) 1.9443
(7 0.0025) 2.1651
(7 0.0024) 2.3590
D3
(0.0410%) 1.2122
(0.0117%) 1.4148
(0.0566%) 1.6132
(0.0046%) 1.7998
(0.0297%) 1.9638
( 7 0.0015) 1.2113
(7 0.0018) 1.4136
(70.0016) 1.6137
(7 0.0025) 1.7999
(7 0.0027) 1.9624
D4
(0.0743%) 0.9766
(0.0849%) 1.1413
(0.0310%) 1.3030
(0.0056%) 1.4553
(0.0713%) 1.5894
( 7 0.0015) 0.9776
(7 0.0016) 1.1408
(70.0018) 1.3036
(7 0.0016) 1.4560
(7 0.0022) 1.5899
(0.1023%)
(0.0438%)
(0.0460%)
(0.0481%)
(0.0314%)
( 7 0.0019)
(7 0.0014)
(70.0012)
(7 0.0013)
(7 0.0015)
Pb/k
X01 X02 X03 X04
Table 6 Example of the expected cost per unit time with control values S and Q. R ¼3, (N1,N2,N3) ¼ (9,8,7), (D1,D2,D3) ¼ (8,7,6), (l1,l2,l3)¼ (0.3,0.4,0.5), t ¼ 0.5, M ¼5, m ¼ 0.5, CD(1) ¼ CD(2) ¼CD(3) ¼140, CQ ¼ 8, CH ¼4, CP ¼ 40 Q
1 2 3 4 5 6 7 8 9 10 Q*(S) Cost(S,Q*(S))
S 0
1
2
3
4
5
6
7
2140.88 2134.57 2136.85 2142.65 2150.79 2160.78 2172.31 2185.10 2198.66 2213.36 2 2134.57
2139.62 2132.02 2132.64 2136.59 2142.82 2150.91 2160.55 2171.60 2183.59 2196.61 2 2132.02
2139.08 2130.61 2129.98 2132.40 2136.95 2143.28 2151.14 2160.45 2170.89 2182.29 3 2129.98
2139.20 2130.17 2128.60 2129.80 2132.90 2137.66 2143.89 2151.52 2160.40 2170.28 3 2128.60
2140.08 2130.60 2128.32 2128.53 2130.43 2133.81 2138.57 2144.64 2151.99 2160.40 3 2128.32
2141.75 2131.84 2128.96 2128.39 2129.29 2131.49 2134.96 2139.62 2145.51 2152.51 4 2128.39
2144.17 2133.80 2130.42 2129.18 2129.28 2130.50 2132.82 2136.24 2140.78 2146.44 4 2129.18
2147.20 2136.38 2132.54 2130.75 2130.18 2130.60 2131.96 2134.28 2137.62 2142.01 5 2130.18
CD(r), is changed. The experimental results together with the experimental data are summarized in Tables 7–9. In Tables 7–9, the optimal (S,Q) inventory policy is compared with the (S¼0, Q¼1) inventory policy, which was used in the study of Abboud and Daigle [24]. We can verify from Tables 7–9, the optimal (S,Q)
8 2150.66 2139.46 2135.21 2132.97 2131.83 2131.58 2132.15 2133.56 2135.84 2139.08 6 2131.58
Optimal inventory policy when S ¼4, Q¼ 3, with cost 2128.32
policy saves a cost considerably compared to the (S¼0, Q¼1) policy. From Table 7, we can check that both the optimal S value and the optimal Q value increase as the ordering cost, CQ, increases. Thus, Table 7 shows that the batch ordering policy is better than
1594
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
Table 7 Sensitivity analysis of the ordering cost (scenario set E). R ¼3, (N1,N2,N3) ¼ (9,8,7), (D1,D2,D3) ¼ (8,7,6), (l1,l2,l3)¼ (0.3,0.4,0.5), t ¼ 0.5, M ¼5, m ¼0.5 CD(1) ¼ CD(2) ¼CD(3)¼ 50, CH ¼4, CP ¼ 40 CQ
4 8 12 16 20
Optimal (S,Q) inventory policy
(S¼ 0,Q¼ 1) Cost
Optimal S
Optimal Q
Optimal Cost
3 4 4 4 5
4 5 5 5 7
825.7143 828.1854 830.1688 831.8388 833.3581
831.5085 841.4738 851.4392 861.4045 871.3699
Table 8 Sensitivity analysis of the holding cost (scenario set F). R ¼3, (N1,N2,N3) ¼ (9,8,7), (D1,D2,D3) ¼ (8,7,6), (l1,l2,l3)¼ (0.3,0.4,0.5), t ¼ 0.5, M ¼5, m ¼0.5 CD(1) ¼ CD(2) ¼CD(3)¼ 50, CQ ¼ 7, CP ¼40 CH
3 5 7 9 11
Optimal (S,Q) inventory policy
(S¼ 0, Q¼ 1) Cost
Optimal S
Optimal Q
Optimal cost
4 3 2 2 2
5 4 4 4 4
827.2936 827.7866 828.1702 828.3375 828.5047
838.9825 838.9825 838.9825 838.9825 838.9825
Table 9 Sensitivity analysis of the shortage penalty cost (scenario set G). R ¼3, (N1,N2,N3) ¼ (9,8,7), (D1,D2,D3) ¼ (8,7,6), (l1,l2,l3)¼ (0.3,0.4,0.5), t ¼ 0.5, M ¼5, m ¼0.5 CQ ¼7, CH ¼4, CP ¼40 CD(r) (r¼ 1,2,3)
40 60 80 100 120
Optimal (S,Q) inventory policy Optimal S
Optimal Q
Optimal cost
3 3 4 4 4
5 4 4 4 3
682.8447 972.1981 1261.2003 1550.1026 1838.8963
(S¼ 0,Q¼ 1) Cost
694.6046 983.3604 1272.1163 1560.8722 1849.6280
the one-by-one ordering policy if the ordering cost, CQ, is relatively high. As shown in Table 8, the optimal S value and the optimal Q value decrease as the holding cost, CH, increases. Finally, from Table 9, we can observe that the optimal S value increases while the optimal Q value decreases as the penalty cost, CD(r), increases.
7. Conclusions In this study we considered a multi-echelon repair system with multiple bases in which a parts inventory system is incorporated. While the previous studies assumed that either parts for repair always exist or are ordered whenever an item fails, in this study, the (S,Q) inventory policy is assumed. In addition, to build the model in more realistic settings, it is assumed that there are multiple bases. We modelled the proposed system as a multiclass closed queueing network with a synchronized station and developed an approximation method based on Marie’s method. By taking advantage of the special structure of the routing probability, we developed an algorithm in which the computational effort required to estimate the parameters of the equivalent network is comparable to that for a single class network. In analyzing sub-network R0, a recursive technique was used.
Comparisons with simulation showed that the approximation method provided fairly good estimation of the desired performance measures. Also, we evaluated the cost effectiveness of the (S,Q) inventory policy with numerical results. Then we proposed a procedure to minimize the costs relevant to the (S,Q) inventory policy. For future studies, we can consider an optimal design of the multi-echelon repair system to determine the stock level of the repairable items and the repair capacity. References [1] Sherbrooke CC. METRIC: a multi-echelon technique for recoverable item control. Operations Research 1968;16(1):122–41. [2] Muckstadt JA. A model for a multi-item, multi-echelon, multi-indenture inventory system. Management Science 1973;20(4):472–81. [3] Graves SC. A multi-echelon inventory model for repairable items with one-for one replenishment. Management Science 1985;31(10):1247–56. [4] Schultz CR. Spares parts inventory and cycle time reduction. International Journal of Production Research 2004;42(4):759–76. [5] Dı´az A, Fu MC. Multi-echelon models for repairable item: a review. Working paper, University of Maryland, 1995. [6] Kennedy WJ, Patterson JW, Fredendall LD. An overview of recent literature on spare parts inventories. International Journal of Production Economics 2002;76(2):201–15. [7] Lau HC, Song HW, See CT, Cheng SY. Evaluation of time-varying availability in multi-echelon spare parts systems with passivation. European Journal of Operational Research 2006;170(1):91–105. [8] De Smidt-Destombes KS, Van Elst NP, Barros AL, Mulder H. A spare parts model with cold-standby redundancy on system level. Computers & Operations Research 2011;38(7):985–91. [9] Kranenburg AA, Van Houtum GJ. A multi-item spare parts inventory model with customer differentiation. Working paper, Technische Universiteit Eindhoven, 2004. [10] Hill RM, Seifbarghy M, Smith DK. A two-echelon inventory model with lost sales. European Journal of Operational Research 2007;181(2):753–66. [11] Al-Rifaia MH, Rossettib MD. An efficient heuristic optimization algorithm for a two-echelon (R, Q) inventory system. International Journal of Production Economics 2007;109(1–2):195–213. [12] Topan E, Bayindir ZP, Tan T. An exact solution procedure for multi-item twoechelon spare parts inventory control problem with batch ordering in the central warehouse. Operations Research Letters 2010;38(5):454–61. [13] Gross D, Miller DR, Soland RM. A closed queueing network model for multiechelon repairable item provisioning. IIE Transactions 1983;15(4):344–52. [14] Madu CN. A closed queueing maintenance network with two repair centres. Journal of the Operational Research Society 1988;39(10):959–67. [15] Albright SC, Gupta A. Steady-state approximation of a multiechelon multiindentured repairable-item inventory system with a single repair facility. Naval Research Logistics 1993;40(4):479–93. [16] Sleptchenko A, Van Der Heijden MC, Van Harten A. Trade-off between inventory and repair capacity in spare part networks. Journal of the Operational Research Society 2003;54(3):263–72. [17] Dı´az A, Fu MC. Models for multi-echelon repairable item inventory systems with limited repair capacity. European Journal of Operational Research 1997; 97(3):392–480. [18] Spanjers L, Van Ommeren JCW, Zijm WHM. Closed loop two-echelon repairable item systems. OR Spectrum 2005;23(2-3):369–98. [19] Jung BR, Sun BG, Kim JS, Ahn SE. Modeling lateral transshipments in multiechelon repairable item inventory systems with finite repair channels. Computers & Operations Research 2003;30(9):1401–17. [20] Kim JS, Kim TY, Hur S. An algorithm for repairable item inventory with depot spares and general repair time distribution. Applied Mathematical Modelling 2007;31(5):795–804. [21] Rappold JA, Van Roo BD. Designing multi-echelon service parts networks with finite repair capacity. European Journal of Operational Research 2009; 199(3):781–92. [22] De Smidt-Destombes KS, Van Der Heijden MC, Van Harten A. Availability of k-out-of-N systems under block replacement sharing limited spares and repair capacity. International Journal of Production Economics 2007;107(2): 404–21. [23] De Smidt-Destombes KS, Van Der Heijden MC, Van Harten A. Joint optimization of spare part inventory, maintenance frequency and repair capacity for k-out-of-N systems. International Journal of Production Economics 2009; 108(1):260–8. [24] Abboud NE, Daigle JN. A little’s result approach to the service constrained spares provisioning problem for repairable items. Operations Research 1997;45(4):577–83. [25] Baynat B, Dallery Y. A unified view of product-form approximation techniques for general closed queueing networks. Performance Evaluation 1993; 18(3):205–24. [26] Baskett F, Chandy KM, Muntz RR, Palacios FG. Open, closed and mixed networks of queues with different classes of customers. Journal of the ACM 1975;22(2):248–60.
C.-W. Park, H.-S. Lee / Computers & Operations Research 38 (2011) 1584–1595
[27] Marie RA. An approximate analytical method for general queueing networks. IEEE Transactions on Software Engineering 1979;5(5):530–8. [28] Neuse D, Chandy KM. HAM: the heuristic aggregation method for solving general closed queueing network models of computer systems. ACM SIGMETRICS Performance Evaluation Review 1982;11(4):195–212. [29] Baynat B, Dallery Y. A product-form approximation method for general closed queueing networks with several classes of customers. Performance Evaluation 1996;24(3):165–88.
1595
[30] Bruell SC, Balbo G. Computational algorithms for closed queueing networks. New York: Elsevier North-Holland; 1980. [31] Herzog U, Woo LS, Chandy KM. Solution of queuing problems by a recursive technique. IBM Journal of Research and Development 1975;19(3):295–300. [32] Buzacott JA, Kostelski D. Matrix-geometric and recursive algorithm solution of a two-stage unreliable flow line. IIE Transactions 1987;19(4):429–38. [33] Law AM, Kelton WD. Simulation modeling and analysis. New York: McGrawHill; 1999.
