- Email: [email protected]
Sensitivity formulas of performance in two-server cyclic queueing networks with phase-type distributed service times
Intl. Trans. in Op. Res. 6 (1999) 649±663
www.elsevier.com/locate/orms
Sensitivity formulas of performance in two-server cyclic queueing networks with phase-type distributed service p times Baoqun Yin*, Yaping Zhou, Hongsheng Xi, Demin Sun Department of Automation, China University of Science and Technology, Hefei 230026, People's Republic of China Received 9 May 1997; received in revised form 2 October 1998; accepted 26 May 1999
Abstract With the approach of the in®nitesimal generator perturbation, sensitivities of the steady-state performance are studied in two-server cyclic queueing networks with phase-type distributed service times. Sensitivity formulas expressed by the potentials of the networks are given for parameterdependent performance functions. An algorithm to compute the potentials and the performance derivatives of the networks is given, and a numerical example is provided. # 1999 IFORS. Published by Elsevier Science Ltd. All rights reserved. Keywords: Sensitivity formulas; Phase-type distributions; Potentials
1. Introduction The approach of perturbation analysis (PA) of closed queueing networks (CQNs) has become an increasingly important research area. The main objective of the PA research is to obtain performance derivatives with respect to the parameters of systems. The approach is widely used in studying practical engineering systems, such as computer systems, communication networks and integrated manufacturing systems (see e.g., Caramanis and Liberopoulos, 1992; Haurie et al., 1994). The traditional PA (IPA) is very ecient for some systems (see e.g., Ho and Cao, 1991; Glasserman, 1991; Cao, 1994; Chong and Ramadge, Supported by Huawei research foundation and China high performance computing foundation. * Corresponding author. E-mail address: [email protected] (B. Yin) p
0969-6016/99/$ - see front matter # 1999 IFORS. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 9 - 6 0 1 6 ( 9 9 ) 0 0 0 1 6 - 7
650
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
1993; Suri and Zazanis, 1988; Tang and Chen, 1994), but for others, it is not (see e.g., Cao, 1987a; Heidelberger et al., 1988). After the hard works of many researchers, many useful processing methods that can be applied to dierent cases where IPA fails have been proposed (see e.g., Cassandras and Strickland, 1989; Ho et al., 1983; Vakili, 1991). The concept of perturbation realization in the PA approach has been successfully applied in the performance sensitivity analysis of queueing networks. With realization factors, the problem of sensitivity analysis of the steady-state performance with respect to the perturbation of a service rate has been solved (see e.g., Cao, 1994). However, it is dicult to use the approach in sensitivity analysis with respect to the perturbation of a routing probability (see e.g., Cao, 1995). It was often assumed that the performance function was independent on the parameters of networks in previous discussions (see e.g., Cao, 1994; Cao and Chen, 1997; Yin et al., 1996). Another problem worth noting is that IPA can only provide unbiased and strongly consistent estimates of the performance derivatives for some CQNs (see e.g., Cao, 1987a; Heidelberger et al., 1988). The available conditions to ensure unbiasedness and strong consistency are not easy (see e.g., Cao, 1987b; Heidelberger et al., 1988; Glasserman et al., 1991; Konstantopoulos and Zazanis, 1992). Therefore, the PA approach needs to be improved. The approach of the perturbation of its in®nitesimal generator for a Markov process, which was recently proposed by Cao and Chen (1997), can be used to solve the above problems. The approach provides new derivative formulas of the steady-state performance of Markov processes with respect to their in®nitesimal generators. The quantity, performance potential, involved in the derivative formulas can be easily estimated by analyzing a single sample path of a Markov process. And the derivative estimates obtained by using the approach are unbiased, and converge a.s. (almost surely) to accurate values of the steadystate performance derivatives. By using this approach, sensitivities of the steady-state performance of a class of CQNs with respect to their parameters (service rate and routing probability) were studied in Yin et al. (1996, 1997), and sensitivity formulas were given for parameter-dependent performance functions. This approach provides a united framework for both IPA and non-IPA approaches. In this paper, we generalize the corresponding results in Yin et al. (1997) to the case of the Markov-type queueing networks whose service distributions are not exponential. For general Markov-type queueing networks with non-exponential distributions (for example, PHdistribution), their state spaces, obtained by using the method of enlarging state space, are very complicated. We only consider a simple case, two-server cyclic queueing network. The processing method for a multi-server closed queueing network is similar. We ®rst state some results for performance sensitivity analysis of Markov processes in Ref. (Cao and Chen, 1997) and give sensitivity formulas of performance with respect to the parameters of the systems for parameter-dependent performance functions. We then brie¯y review the concept of phase-type distribution, and give the expression of the in®nitesimal generator of a two-server cyclic queueing system. Finally, we derive sensitivity formulas of performance which can be expressed by the potentials of the systems. An algorithm to calculate potentials and performance derivatives based on a single sample path is given, and a numerical example is also provided to show the application of the algorithm.
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
651
2. Sensitivity formulas for a class of Markov processes In this section, we review two fundamental concepts and quantities, realization matrix and performance potentials. It is shown that the sensitivity of the steady-state performance can be easily calculated by using the potentials of the systems. 2.1. The performance Consider a positive recurrent, irreducible Markov process Y fYt ;tr0g. It has a ®nite state space F f1,2, . . . ,K g with an in®nitesimal generator A aij , i,j 2 F. By the theory of Markov processes, Y has a unique steady-state probability measure. Denote by P p 1, 2, . . . ,p K the row vector representing the probability. Let Rk be the space of the kdimensional real numbers, and simply denote R 1 by R; let e 1,1, . . . ,1 0 , where the prime represents the transpose. Then, we have Pe 1, and Ae 0, PA 0:
1
Let f:F4R be a performance function. The performance measure Zf of Y is de®ned as its expected value with respect to the steady-state probability measure P, i.e., Zf
K X p i f i Pf:
2
i1
where the last f f 1,f 2, . . . ,f K 0 is a column vector. f represents both a function and a vector. Now, suppose that the in®nitesimal generator A is a function of the parameter y de®ned on the interval J R. We also assume that the performance function f depends on y. The quantities Zf , P, then, are also functions of the parameter y. We assume that all aij 's and f are dierentiable functions on J. 2.2. The realization matrix and potentials In this subsection, we brie¯y review the concepts of group inverse of the in®nitesimal generator, realization matrix and potential in Cao and Chen (1997). Since the in®nitesimal generator A in Eq. (1) is not of full rank, its regular inverse does not exist. For every P satisfying Eq. (1), the corresponding A is not unique. De®ne 3 AP A£jA 0,Ae 0 It is easy to verify that AP is a group with identity element I ÿ eP on matrix multiplication. De®nition 1. The inverse of A in the group AP is called the group inverse of the in®nitesimal generator A, denoted as A #. Lemma 1. The matrix A eP is nonsingular.
652
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
It is easy to verify A A ePÿ1 ÿ eP A ePÿ1 ÿePA I ÿ eP. Therefore, we have A# A eP ÿ1 ÿeP:
4
For any i,j 2 F, let Y fi g fYt :Y0 i; tr0g and Y fj g fYt :Y0 j; tr0g be two Markov processes with the initial states i and j, respectively. De®ne Z fi,j g Y fi g ,Y fj g . Then, Z fi,j g is a Markov process with state space F F. Z fi,j g is positive recurrent, irreducible since Y fi g and Y fj g are. Thus, the ®rst passage time from any state to any other state has a ®nite mean. Let g 2 Sg. S f k,k:k 2 Fg. De®ne T fi,j g infft:tr0,Z fi,j t De®nition 2. ( dij E
) i ÿ fjg f Y t ÿ f Y fti g dt ,
T fi,jg h
0
i,j 2 F,
is called a perturbation realization factor of the Markov process Y with respect to f; the matrix D f dij is called a realization matrix. From the de®nition, we have D f 0 ÿD f . Lemma 2.
# " ( " #) T ÿ T dt ÿ E dij lim E f Y fjg f Y fti g dt , t T41
0
0
i,j 2 F:
5
From Lemma 2, we have dij dik dkj ,
i,j,k 2 F:
6
This property is similar to that of the potential energy in physics. In view of this, we may pick up any integer k 2 F and any real number c and de®ne a quantity, x i , for any i 2 F, i 6 k , as follows: x k c,
x i x k dk i :
7
De®nition 3. x i , i 2 F is called the potential of the Markov process Y at state i with respect to the performance function f; the column vector x f x 1 ,x 2 , . . . ,x k 0 is called the potential vector, or simply the potential, of the Markov process Y with respect to f.
From Eq. (6), we can know dij x j ÿ x i , i,j 2 F, or equivalently,
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663 0
D f e x f ÿx f e 0 :
Lemma 3. The realization matrix D f Lyapunov equation
653
8
is the unique solution in V fex 0 ÿ xe 0 : x 2 RK g to the
AD DA 0 ÿF:
9
where F ef 0 ÿ fe 0 :
The proofs of Lemma 1±3 can be found in Cao and Chen (1997). Lemma 4. The set of all potential vectors x f is the solution space fx f :x f ÿA# f be,b 2 Rg to the linear equation Ax ÿf Zf e:
10
Proof. Substituting Eq. (8) into AD f D f A 0 ÿF, and noting that F ef 0 ÿ fe 0 , we have 0 ÿ ÿ 11 e Ax f f Ax f f e 0 : Because the solution to equation ez 0 ze 0 is z be, b 2 R, Eq. (11) is equivalent to equation Ax f f be, b 2 R. Note that PA 0, Pf Zf , and Pe 1, we must have b Zf in equation Ax f f be, b 2 R. Therefore, x f is the solution to Eq. (10).
Because P satisfying Pe 1, PA 0 is unique, the rank of A is K ÿ 1. It is easy to know that the solution to equation Ax 0 is x be, b 2 R. Note that x ÿA# f is a solution to Eq. (10), therefore, all solutions to Eq. (10) are x ÿA# f be, b 2 R. In addition, it is easy to verify that ÿA# f be is the potential of Y with respect to f for each b. From this, we get Lemma 4q. From Eq. (10), we have @A f @x f @f @Zf x A ÿ e: @y @y @y @y Multiplying both sides on the left by P, since PA 0 and Pe 1, we get @Zf @A @f P x f P : @y @y @y
12
654
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
3. Phase-type distributions In this section, we brie¯y review the concept of phase-type distribution (see e.g., Neuts, 1981). Consider a Markov process on the states f1,2, . . . ,m 1g with an in®nitesimal generator T t0 ^ : 13 Q 0 0 where the m m matrix T satis®es Tii <0, for 1RiRm, and Tij r0, for i 6 j. Also the mdimensional vector t0 satis®es Te t0 0, with e 1,1, . . . ,1 0 . The initial probability vector of Q^ is given by (a,am1 ), with ae am1 1. We assume that states 1,2, . . . ,m are all transient, so that absorption into the state m 1, from any initial state, is certain. De®nition 4. A probability distribution F x on 0,1 is an m-order phase-type distribution (PH-distribution) if and only if it is the distribution of the time until absorption in a ®nite Markov process of the type de®ned in Eq. (13) and with initial probability vector (a,am1 ). The pair (a,T) is called a representation of F x.
A useful equivalent condition that the states 1,2, . . . ,m are all transient is given by the following lemma (see e.g., Neuts, 1981). Lemma 5. The states 1,2, . . . ,m are transient if and only if the matrix T is nonsingular. Lemma 6. The probability distribution F x in De®nition 4 is given by F x 1 ÿ a exp Tx e,
for xr0:
14
Now, we consider an example of PH-distribution, which has an explicit physical meaning. Consider a service station consisting of k phase, each having an exponential distribution with service rate miP , i 1,2, . . . ,k. The customer arriving at the station enters phase i with k the customer goes to probability ai , i1 ai 1. After the completion of its service at phase i, P k phase j with probability qi,j and leaves the station with probability qi , j1 qi,j qi 1. Let Q qi,j . We assume that the matrix Ik ÿ Q is nonsingular so that the probability of a customer leaving the station is one. Where Ik represents the k-order identity matrix. Then, the service time of the customer at the station has a k-order phase-type distribution with representation (a,T), where a a1 ,a2 , . . . ,ak , and T diag m1 ,m2 , . . . ,mk Q ÿ Ik .
4. System models and sensitivity formulas Consider a cyclic queueing network consisting of two single FCFS server nodes and N single-class customer; each server has a buer of an in®nite capacity. The customer service
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
655
times at two servers are PH-distributions with m1-order representation (a,T) and m2-order representation (b,S), respectively, i.e., two distributed functions of the service times are given by F1 x 1 ÿ a exp Tx e,
15
F2 x 1 ÿ b exp Sx e,
16
and
respectively. Where a and b are probability vectors of m1-dimensional and m2-dimensional, respectively; T and S are two nonsingular matrices, and satisfy Tii <0, Sii <0 and Tij r0, Sij r0 for i 6 j. Let (O,F,P) be the probability space on which ransom variables of service times are de®ned. Then, the state process of such a network can be described by a QBD process Y fYt ; tr0g with a ®nite state space: E 0,j: 1RjRm1 [ i,j,k: 1RiRN ÿ 1,1RjRm1 ,1RkRm2 [ N,k : 1RkRm2 :
17
The index i, 1RiRN ÿ 1, denotes the number of customers in server 2; the index j, 1RjRm1 , represents the phase of service in server 1; while the index k, 1RkRm2 indicates the phase of service in server 2. We assume that Y is irreducible. The states are labeled in the lexicographic order, that is, 0,1, . . . , 0,m1 , 1,1,1, . . . , 1,1,m2 , 1,2,1, . . . , 1,2,m2 , . . . , 1,m1 ,m2 , . . . , N,1, . . . , N,m2 . Then, the in®nitesimal generator of the QBD process Y is given by 2 3 T t0 a b 0 0 0 0 0 6 Im1 S 0 A1 7 A0 0 0 0 0 6 7 60 7 A2 A1 A0 0 0 0 6 7 18 A6 7 6 7 40 t0 Im2 5 0 0 0 A2 A1 0 0 0 0 0 0 a s b S where t0 ÿTe, s0 ÿSe; t0 a represents the multiplication of m1 1 matrix t 0 and 1 m1 matrix a, and s0 b is similar. represents the Kronecker product of matrices. The matrices A0 ,A1 ,A2 are given by A0 t0 a Im2 ,
19
A1 T Im2 Im1 S,
20
A2 Im1 s0 b:
21
and
656
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
Since the state space is ®nite, Y is positive recurrent. By the theory of Markov process, Y has a unique steady-state probability vector P P0 ,P1 , . . . ,PNÿ1 ,PN , where P0 is an m1-dimensional row vector; P1 , . . . ,PNÿ1 are all of m1 m2 -dimensional; PN is of m2-dimensional. Let F fn: 0RnRN g be the discrete state space of the QBD process Y, where n represents the number of customers in server 2. Assume that f:F4R is a performance function. The performance measure Zf of the network is de®ned as its expected value with respect to the steady-state measure P, i.e., Zf
N X f np n Pf:
22
n0
where p n jPn j, n 0,1, . . . ,N, with jPn j representing the sum of all components of the vector Pn ; f f0 ,f1 , . . . ,fN 0 is an m1 N ÿ 1m1 m2 m2 -dimensional column vector; fn is the same dimensional vector as Pn , and its components are all f n. Now, we suppose that vectors a,b and matrices T,S are all dierentiable functions of parameter y on the interval J R. Then, the matrix A is also dierentiable on J. We also assume that f depends on y, and is a dierentiable function of y.
5. Some special cases In this section, we discuss sensitivity formulas of the networks with phase-type distributions of special forms, which have explicit physical meaning. 5.1. Phase-type distributions We ®rst consider the models of the phase-type distributions given in the end of section Phase-Type Distributions. We assume that two servers have service distributions of such a form with representation a a1 ,a2 , . . . ,am1 , T diag m11 ,m12 , . . . ,m1m1 Q1 ÿ Im1 and b b1 ,b2 , . . . ,bm2 , S diag m21 ,m22 , . . . ,m2m2 , respectively. Where Ql q li,j , l 1,2. Now, we only consider sensitivities with respect to service rates m1j , j 1,2, . . . ,m1 . Denote by Ti the ith row vector of the matrix T. It is easy to get that 2 3 0 6. 7 6 .. 7 6 7 6 7 60 7 @T 6 7 23 m1i 6 Ti 7, @m1i 6 7 60 7 6 7 6 .. 7 4. 5 0
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
2
3
0 7 6. 7 6 .. 7 6 7 6 0 7 6 @ 0 7 1 6 t a b m1i qi 6 a b 7, m1i 7 6 @m1i 7 60 7 6 7 6 .. 5 4. 0 2
0 6. 6 .. 6 6 60 @A 6 m1i 6 Ti Im2 @m1i 6 60 6 6 .. 4. 0 and
2
657
24
3 7 7 7 7 7 7 7, 7 7 7 7 5
0 6. 6 .. 6 6 60 @A0 6 m1i m1i q i 1 6 a Im2 6 @m1i 60 6 6 .. 4. 0
25
3 7 7 7 7 7 7 7: 7 7 7 7 5
26
where q i 1 1 ÿ
m1 X j1
q i,j1 ,
i 1,2, . . . ,m1 . Pk pk1 ,pk2 , . . . ,pkm1 for k 1,2, . . . ,N ÿ 1 with Let P0 p01 ,p02 , . . . ,p0m1 , pij pij1 ,pij2 , . . . ,pijm2 , j 1,2, . . . ,m1 and PN pN1 ,pN2 , . . . ,pNm2 . From Eqs. (12), (23)±(26), we then have m1i where
Nÿ1 X @Zf @f Pk B ki x f m1i P , i 1,2, . . . ,m1 , @m1i @m1i k0
27
658
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
h
B 0i Ti
m1i q i 1 a b
i 0 ,
0
28
" B ki
#
and
Ti Im2
m1i q i 1 a
Im2
" i B Nÿ1
0
0 , k 1,2, . . . ,N ÿ 2,
29
# Ti Im2
m1i q i 1 Im2
:
30
@Z
Similarly, the expression for m2i @ m f can be derived. 2i
5.2. Coxian distributions For the above phase-type distribution, let a1 1, ai 0, i 2,3, . . . ,k, and qi,i1 ai , qi,j 0, j 6 i 1, i 1,2, . . . ,k ÿ 1, qk 1, then we get the Coxian distribution Ck . It is a phase-type distribution with representation a 1,0, . . . ,0, and 2 3 ÿ1 a1 0 6 7 ÿ1 6 7 6 7 6 7: T diag m1 ,m2 , . . . ,mk 6 7 6 7 4 akÿ1 5 0 ÿ1 Now suppose that the service distributions of two servers are the Coxian distributions Cm1 and Cm2 with representations (a,T) and (b,S), respectively. Where a 1,0, . . . ,0, b 1,0, . . . ,0, and 2 3 ÿ1 a1 0 6 7 ÿ1 6 7 6 7 7, T diag m11 ,m12 , . . . ,m1m1 6 6 7 6 7 4 am1 ÿ1 5 0 ÿ1 2 6 6 6 S diag m21 ,m22 , . . . ,m2m2 6 6 6 4
ÿ1
0 In this case, Eqs. (28)±(30) become
b1 ÿ1
0
3
7 7 7 7: 7 7 bm2 ÿ1 5 ÿ1
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
" B 0i
ÿm1i
m1i 1 ÿ ai
m1i ai
0
659
#
0 ,
"
#
B ki
ÿm1i Im2
m1i ai Im2
m1i 1 ÿ ai Im2
0
0 ,
k 1,2, . . . ,N ÿ 2, and
" i B Nÿ1
#
ÿm1i Im2
m1i ai Im2
m1i 1 ÿ ai Im2 :
where am1 0. 5.3. Erlang distributions In the k-order Coxian distribution Ck , let q1 q2 qkÿ1 0, m1 m2 mk km, we then get the k-order Erlang distribution Ek . Its representation is given by a 1,0, . . . ,0 and 3 2 ÿkm km 0 7 6 ÿkm 7 6 7 6 7 kmJk , T6 7 6 7 6 4 km 5 0 ÿkm where
2
ÿ1 1 6 ÿ1 6 6 Jk 6 6 6 4 0
0
3
7 7 7 7: 7 7 1 5 ÿ1
Now, we suppose that the service distributions of two servers are Erlang distributions Em1 and Em2 . Their representations are given by (a,T) and (b,S), respectively. Where a 1,0, . . . ,0, b 1,0, . . . ,0, T m1 m1 Jm1 and S m2 m2 Jm2 . Then, we have t0 ÿTe 0, . . . ,0,m1 m1 0 . From this, we can obtain m1
@T T m1 m1 Jm1 , @m1
660
2
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
0 6 @ 0 m1 t a b 6 40 @m1 m1 m 1 m1 and
0 0 0
3 0 7 7, 05 0
@A1 T Im2 m1 m1 Jm1 Im2 , @m1 2
0 @A0 6 6 m1 4 0 @m1 m1 m1 Im2
0 0 0
3 0 7 7: 05 0
Substituting the above equations into Eq. (12), we can get the sensitivity formulas for the case of Erlang distribution. 5.4. Hyperexponential distributions For the phase-type distribution of Section 5.1, let qi,j 0, for all i, j, i.e., Q 0, we then get the k-order Hyperexponential distribution Hk . Its representation is given by a a1 ,a2 , . . . ,ak , and T diag ÿm1 , ÿ m2 , . . . , ÿ mk . Now we assume that the service distributions of two servers are Hm1 and Hm2 , whose T diag ÿm11 , ÿ m12 , . . . , ÿ m1m1 , and representations are a a1 ,a2 , . . . ,am1 , b b1 ,b2 , . . . ,bm2 , S diag ÿm21 , ÿ m22 , . . . , ÿ m2m2 . In this case, substituting Ti 1 1 0, . . . ,0, ÿ m1i ,0, . . . ,0, and q 1 1 q2 qm1 1 into Eqs. (28)±(30), we can obtain the i expressions of B k , k 0,1, . . . ,N ÿ 1. Finally the sensitivity formulas can be obtained. 6. An algorithm and an example In this section, we provide an algorithm to compute the potential x f and the performance derivative @Zf =@mij of the network. This algorithm is directly applicable to the controlling and optimization of CQNs, since it is based on analyzing a single sample path. A numerical example is provided to illustrate the application of the algorithm. Let S fj g i infft:tr0,Y fjt g ig be the ®rst passage time from state j to state i for the Markov process Y. By Theorem 2 in Cao and Chen (1997), we have ( ) i S fjg i h ÿ ÿ Zf dt : f Y fjg 31 dij E t 0
It is easy to verify that x f ÿD f P 0
32
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
661
is a potential vector of the network.From Eq. (27), we have Nÿ1 X
@Zf @m1i
k0
Pk B ki x f m1i
P
@f , i 1,2, . . . ,m1 : @m1i
33
6.1. Algorithm 1. 2. 3. 4.
S fj g i Estimate E 0 f Y fjt g dt, ES fj g i and Zf on a single sample path. Calculate dij 's and the realization matrix D f by using Eq. (31). Set x f ÿD f P 0 . Calculate @Zf =@mij by using Eq. (33).
6.2. Example Consider a two-server cyclic queueing networks with N 3, each having the two-order Coxian distribution C2. In this case, F f0,1,2,3g, and there are 12 states in E. The parameters of the network are given by m11 1, m12 2,
m21 2, m22 3, a1 b1 0:5:
The performance function is f n n, for n 2 F; i.e., Zf is the steady-state mean response time of server 2. The means of the potentials for simulation run of 1,000,000 times are as follows: ÿ1:1219 2:1111
ÿ0:3779 0:8534
ÿ0:0014 3:3763
ÿ0:7298 2:1343
1:2028 0:1791 3:9805 2:8676
The theoretical values are given by ÿ1:1089 2:1061
ÿ0:3671 0:8466
0:0038 ÿ0:7246 1:2026 0:1700 3:3668 2:1167 3:9879 2:8588
The mean of the performance derivative @Zf =@m11 for same run is given by 0.8292, while its theoretical value is 0.8259. The standard derivations of our estimates are reasonably small.
7. Conclusions In this paper, with the approach of the in®nitesimal generator perturbation, we studied sensitivities of the steady-state performance for two-server cyclic queueing networks. Each server had the service time of the phase-type distribution. With the potentials of the networks, sensitivity formulas of the performance with respect to the parameters of the systems were given for parameter-dependent performance functions. The main contribution of this paper is that we apply the approach of the in®nitesimal generator perturbation to the queueing
662
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
networks with nonexponential distributions, and provide a method to solve the problems of sensitivity analysis for such a class of networks. Especially, we derived sensitivity formulas of performance with PH-distributions of special forms for two-server cyclic queueing networks. These formulas provide a theoretical basis with the single sample path-based performance sensitivity analysis of CQNs. The system models with PH-distributions of special forms can be directly applied to many practical engineering systems. The accurate estimations of the potentials in these formulas can be obtained from a single sample path, therefore, these formulas can be used directly in optimization problems for these systems. Since a Markov system is the fundamental model for many CQNs, the above approach is expected to apply to a wide class of systems. In future, we shall consider multi-server closed queueing networks with PH-distributions. The key point may be the expression of the in®nitesimal generator of the QBD process describing this system. We shall extend the above results to other system models such as generalized semi-Markov processes (GSMP). The key point of extension is that we need to ®nd a matrix for a GSMP, it corresponds with the in®nitesimal generator in a Markov process. Another future research directions is the application of the potential in stochastic optimization.
Acknowledgements The authors would like to thank the referees for their valuable comments and suggestions which helped to improve the paper.
References Cao, X.R., 1987a. First-order perturbation analysis of a single multi-class ®nite source queue. Performance Evaluation 7, 31±41. Cao, X.R., 1987b. Realization probability in closed Jacksson networks and its application. Advanced in Applied Probability 19, 708±738. Cao, X.R., 1994. Realization Probabilities: The Dynamics of Queueing Systems. New York, Springer-Verlag. Cao, X.R., 1995. An overview of perturbation analysis. In: Proceedings of Chinese Control Conference, pp. 22±39. Cao, X.R., Chen, H.F., 1997. Perturbation realization, potentials, and sensitivity analysis of Markov processes. IEEE Transactions on Automatic Control 42, 1382±1393. Caramanis, M., Liberopoulos, G., 1992. Perturbation analysis for the design of ¯exible manufacturing system ¯ow controllers. Operations Research 40, 1107±1125. Cassandras, C.G., Strickland, S.G., 1989. On-line sensitivity analysis of Markov chains. IEEE Transaction on Automation and Control 34, 76±86. Chong, E.K.P., Ramadge, P.J., 1993. Optimization of queues using in®nitesimal perturbation analysis-based algorithms with general update times. SIAM J. of Control and Optimization 31, 698±732. Glasserman, P., 1991. Gradient Estimates via Perturbation Analysis. Kluwer, Boston. Glasserman, P., Hu, J.Q., Strickland, S.G., 1991. Strongly consistent steady-state derivative estimates. Probability in the Engineering and Informational Science 5, 391±413. Haurie, A., L'Ecuyer, P., Delft, C., 1994. Convergence of stochastic approximation couple with perturbation analysis in a class of manufacturing ¯ow control models. Discrete Event Dynamic Systems: Theory and Application 4, 87±111.
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
663
Heidelberger, P., Cao, X.R., Zazanis, M.A., Suri, R., 1988. Convergence properties of in®nitesimal perturbation analysis estimates. Management Science 34, 1281±1302. Ho, Y.C., Cao, X.R., Cassandras, C.G., 1983. In®nitesimal and ®nite perturbation analysis for queueing networks. Automatica 19, 439±445. Ho, Y.C., Cao, X.R., 1991. Perturbation Analysis of Discrete-Event Dynamic Systems. Boston, Kluwer. Konstantopoulos, P., Zazanis, M., 1992. Sensitivity analysis for stationary and ergodic queues. Advanced in Applied Probability 24, 738±750. Neuts, M.F., 1981. Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach. The Johns Hopkins University Press, Baltimore. Suri, R., Zazanis, M.A., 1988. Perturbation analysis gives strongly consistent sensitivity estimates for the M/G/1 queue. Management Science 34, 39±64. Tang, Q.Y., Chen, H.F., 1994. Convergence of perturbation analysis based optimization algorithm with ®xed number of customers period. Discrete Event Dynamic Systems: Theory and Application 4, 359±375. Vakili, P., 1991. A standard clock technique for ecient simulation. Operations Research Letters 10, 445±452. Yin, B.Q., Zhou, Y.P., Yang, X.X., Xi, H.S., Sun, D.M., 1996. Sensitivity formulas of performance in closed statedependent queueing networks. In: Proceedings of Chinese Control Conference, pp. 671±677. Yin, B.Q., Zhou, Y.P., Xi, H.S., Sun, D.M., 1997. Sensitivity formulas of performance and realization factors with parameter-dependent performance functions in closed queueing networks. In: Proceedings of The Second Chinese World Congress on Intelligent Control and Intelligent Automation, Xi'an, P. R. China, pp. 1884±1889.
www.elsevier.com/locate/orms
Sensitivity formulas of performance in two-server cyclic queueing networks with phase-type distributed service p times Baoqun Yin*, Yaping Zhou, Hongsheng Xi, Demin Sun Department of Automation, China University of Science and Technology, Hefei 230026, People's Republic of China Received 9 May 1997; received in revised form 2 October 1998; accepted 26 May 1999
Abstract With the approach of the in®nitesimal generator perturbation, sensitivities of the steady-state performance are studied in two-server cyclic queueing networks with phase-type distributed service times. Sensitivity formulas expressed by the potentials of the networks are given for parameterdependent performance functions. An algorithm to compute the potentials and the performance derivatives of the networks is given, and a numerical example is provided. # 1999 IFORS. Published by Elsevier Science Ltd. All rights reserved. Keywords: Sensitivity formulas; Phase-type distributions; Potentials
1. Introduction The approach of perturbation analysis (PA) of closed queueing networks (CQNs) has become an increasingly important research area. The main objective of the PA research is to obtain performance derivatives with respect to the parameters of systems. The approach is widely used in studying practical engineering systems, such as computer systems, communication networks and integrated manufacturing systems (see e.g., Caramanis and Liberopoulos, 1992; Haurie et al., 1994). The traditional PA (IPA) is very ecient for some systems (see e.g., Ho and Cao, 1991; Glasserman, 1991; Cao, 1994; Chong and Ramadge, Supported by Huawei research foundation and China high performance computing foundation. * Corresponding author. E-mail address: [email protected] (B. Yin) p
0969-6016/99/$ - see front matter # 1999 IFORS. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 9 - 6 0 1 6 ( 9 9 ) 0 0 0 1 6 - 7
650
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
1993; Suri and Zazanis, 1988; Tang and Chen, 1994), but for others, it is not (see e.g., Cao, 1987a; Heidelberger et al., 1988). After the hard works of many researchers, many useful processing methods that can be applied to dierent cases where IPA fails have been proposed (see e.g., Cassandras and Strickland, 1989; Ho et al., 1983; Vakili, 1991). The concept of perturbation realization in the PA approach has been successfully applied in the performance sensitivity analysis of queueing networks. With realization factors, the problem of sensitivity analysis of the steady-state performance with respect to the perturbation of a service rate has been solved (see e.g., Cao, 1994). However, it is dicult to use the approach in sensitivity analysis with respect to the perturbation of a routing probability (see e.g., Cao, 1995). It was often assumed that the performance function was independent on the parameters of networks in previous discussions (see e.g., Cao, 1994; Cao and Chen, 1997; Yin et al., 1996). Another problem worth noting is that IPA can only provide unbiased and strongly consistent estimates of the performance derivatives for some CQNs (see e.g., Cao, 1987a; Heidelberger et al., 1988). The available conditions to ensure unbiasedness and strong consistency are not easy (see e.g., Cao, 1987b; Heidelberger et al., 1988; Glasserman et al., 1991; Konstantopoulos and Zazanis, 1992). Therefore, the PA approach needs to be improved. The approach of the perturbation of its in®nitesimal generator for a Markov process, which was recently proposed by Cao and Chen (1997), can be used to solve the above problems. The approach provides new derivative formulas of the steady-state performance of Markov processes with respect to their in®nitesimal generators. The quantity, performance potential, involved in the derivative formulas can be easily estimated by analyzing a single sample path of a Markov process. And the derivative estimates obtained by using the approach are unbiased, and converge a.s. (almost surely) to accurate values of the steadystate performance derivatives. By using this approach, sensitivities of the steady-state performance of a class of CQNs with respect to their parameters (service rate and routing probability) were studied in Yin et al. (1996, 1997), and sensitivity formulas were given for parameter-dependent performance functions. This approach provides a united framework for both IPA and non-IPA approaches. In this paper, we generalize the corresponding results in Yin et al. (1997) to the case of the Markov-type queueing networks whose service distributions are not exponential. For general Markov-type queueing networks with non-exponential distributions (for example, PHdistribution), their state spaces, obtained by using the method of enlarging state space, are very complicated. We only consider a simple case, two-server cyclic queueing network. The processing method for a multi-server closed queueing network is similar. We ®rst state some results for performance sensitivity analysis of Markov processes in Ref. (Cao and Chen, 1997) and give sensitivity formulas of performance with respect to the parameters of the systems for parameter-dependent performance functions. We then brie¯y review the concept of phase-type distribution, and give the expression of the in®nitesimal generator of a two-server cyclic queueing system. Finally, we derive sensitivity formulas of performance which can be expressed by the potentials of the systems. An algorithm to calculate potentials and performance derivatives based on a single sample path is given, and a numerical example is also provided to show the application of the algorithm.
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
651
2. Sensitivity formulas for a class of Markov processes In this section, we review two fundamental concepts and quantities, realization matrix and performance potentials. It is shown that the sensitivity of the steady-state performance can be easily calculated by using the potentials of the systems. 2.1. The performance Consider a positive recurrent, irreducible Markov process Y fYt ;tr0g. It has a ®nite state space F f1,2, . . . ,K g with an in®nitesimal generator A aij , i,j 2 F. By the theory of Markov processes, Y has a unique steady-state probability measure. Denote by P p 1, 2, . . . ,p K the row vector representing the probability. Let Rk be the space of the kdimensional real numbers, and simply denote R 1 by R; let e 1,1, . . . ,1 0 , where the prime represents the transpose. Then, we have Pe 1, and Ae 0, PA 0:
1
Let f:F4R be a performance function. The performance measure Zf of Y is de®ned as its expected value with respect to the steady-state probability measure P, i.e., Zf
K X p i f i Pf:
2
i1
where the last f f 1,f 2, . . . ,f K 0 is a column vector. f represents both a function and a vector. Now, suppose that the in®nitesimal generator A is a function of the parameter y de®ned on the interval J R. We also assume that the performance function f depends on y. The quantities Zf , P, then, are also functions of the parameter y. We assume that all aij 's and f are dierentiable functions on J. 2.2. The realization matrix and potentials In this subsection, we brie¯y review the concepts of group inverse of the in®nitesimal generator, realization matrix and potential in Cao and Chen (1997). Since the in®nitesimal generator A in Eq. (1) is not of full rank, its regular inverse does not exist. For every P satisfying Eq. (1), the corresponding A is not unique. De®ne 3 AP A£jA 0,Ae 0 It is easy to verify that AP is a group with identity element I ÿ eP on matrix multiplication. De®nition 1. The inverse of A in the group AP is called the group inverse of the in®nitesimal generator A, denoted as A #. Lemma 1. The matrix A eP is nonsingular.
652
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
It is easy to verify A A ePÿ1 ÿ eP A ePÿ1 ÿePA I ÿ eP. Therefore, we have A# A eP ÿ1 ÿeP:
4
For any i,j 2 F, let Y fi g fYt :Y0 i; tr0g and Y fj g fYt :Y0 j; tr0g be two Markov processes with the initial states i and j, respectively. De®ne Z fi,j g Y fi g ,Y fj g . Then, Z fi,j g is a Markov process with state space F F. Z fi,j g is positive recurrent, irreducible since Y fi g and Y fj g are. Thus, the ®rst passage time from any state to any other state has a ®nite mean. Let g 2 Sg. S f k,k:k 2 Fg. De®ne T fi,j g infft:tr0,Z fi,j t De®nition 2. ( dij E
) i ÿ fjg f Y t ÿ f Y fti g dt ,
T fi,jg h
0
i,j 2 F,
is called a perturbation realization factor of the Markov process Y with respect to f; the matrix D f dij is called a realization matrix. From the de®nition, we have D f 0 ÿD f . Lemma 2.
# " ( " #) T ÿ T dt ÿ E dij lim E f Y fjg f Y fti g dt , t T41
0
0
i,j 2 F:
5
From Lemma 2, we have dij dik dkj ,
i,j,k 2 F:
6
This property is similar to that of the potential energy in physics. In view of this, we may pick up any integer k 2 F and any real number c and de®ne a quantity, x i , for any i 2 F, i 6 k , as follows: x k c,
x i x k dk i :
7
De®nition 3. x i , i 2 F is called the potential of the Markov process Y at state i with respect to the performance function f; the column vector x f x 1 ,x 2 , . . . ,x k 0 is called the potential vector, or simply the potential, of the Markov process Y with respect to f.
From Eq. (6), we can know dij x j ÿ x i , i,j 2 F, or equivalently,
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663 0
D f e x f ÿx f e 0 :
Lemma 3. The realization matrix D f Lyapunov equation
653
8
is the unique solution in V fex 0 ÿ xe 0 : x 2 RK g to the
AD DA 0 ÿF:
9
where F ef 0 ÿ fe 0 :
The proofs of Lemma 1±3 can be found in Cao and Chen (1997). Lemma 4. The set of all potential vectors x f is the solution space fx f :x f ÿA# f be,b 2 Rg to the linear equation Ax ÿf Zf e:
10
Proof. Substituting Eq. (8) into AD f D f A 0 ÿF, and noting that F ef 0 ÿ fe 0 , we have 0 ÿ ÿ 11 e Ax f f Ax f f e 0 : Because the solution to equation ez 0 ze 0 is z be, b 2 R, Eq. (11) is equivalent to equation Ax f f be, b 2 R. Note that PA 0, Pf Zf , and Pe 1, we must have b Zf in equation Ax f f be, b 2 R. Therefore, x f is the solution to Eq. (10).
Because P satisfying Pe 1, PA 0 is unique, the rank of A is K ÿ 1. It is easy to know that the solution to equation Ax 0 is x be, b 2 R. Note that x ÿA# f is a solution to Eq. (10), therefore, all solutions to Eq. (10) are x ÿA# f be, b 2 R. In addition, it is easy to verify that ÿA# f be is the potential of Y with respect to f for each b. From this, we get Lemma 4q. From Eq. (10), we have @A f @x f @f @Zf x A ÿ e: @y @y @y @y Multiplying both sides on the left by P, since PA 0 and Pe 1, we get @Zf @A @f P x f P : @y @y @y
12
654
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
3. Phase-type distributions In this section, we brie¯y review the concept of phase-type distribution (see e.g., Neuts, 1981). Consider a Markov process on the states f1,2, . . . ,m 1g with an in®nitesimal generator T t0 ^ : 13 Q 0 0 where the m m matrix T satis®es Tii <0, for 1RiRm, and Tij r0, for i 6 j. Also the mdimensional vector t0 satis®es Te t0 0, with e 1,1, . . . ,1 0 . The initial probability vector of Q^ is given by (a,am1 ), with ae am1 1. We assume that states 1,2, . . . ,m are all transient, so that absorption into the state m 1, from any initial state, is certain. De®nition 4. A probability distribution F x on 0,1 is an m-order phase-type distribution (PH-distribution) if and only if it is the distribution of the time until absorption in a ®nite Markov process of the type de®ned in Eq. (13) and with initial probability vector (a,am1 ). The pair (a,T) is called a representation of F x.
A useful equivalent condition that the states 1,2, . . . ,m are all transient is given by the following lemma (see e.g., Neuts, 1981). Lemma 5. The states 1,2, . . . ,m are transient if and only if the matrix T is nonsingular. Lemma 6. The probability distribution F x in De®nition 4 is given by F x 1 ÿ a exp Tx e,
for xr0:
14
Now, we consider an example of PH-distribution, which has an explicit physical meaning. Consider a service station consisting of k phase, each having an exponential distribution with service rate miP , i 1,2, . . . ,k. The customer arriving at the station enters phase i with k the customer goes to probability ai , i1 ai 1. After the completion of its service at phase i, P k phase j with probability qi,j and leaves the station with probability qi , j1 qi,j qi 1. Let Q qi,j . We assume that the matrix Ik ÿ Q is nonsingular so that the probability of a customer leaving the station is one. Where Ik represents the k-order identity matrix. Then, the service time of the customer at the station has a k-order phase-type distribution with representation (a,T), where a a1 ,a2 , . . . ,ak , and T diag m1 ,m2 , . . . ,mk Q ÿ Ik .
4. System models and sensitivity formulas Consider a cyclic queueing network consisting of two single FCFS server nodes and N single-class customer; each server has a buer of an in®nite capacity. The customer service
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
655
times at two servers are PH-distributions with m1-order representation (a,T) and m2-order representation (b,S), respectively, i.e., two distributed functions of the service times are given by F1 x 1 ÿ a exp Tx e,
15
F2 x 1 ÿ b exp Sx e,
16
and
respectively. Where a and b are probability vectors of m1-dimensional and m2-dimensional, respectively; T and S are two nonsingular matrices, and satisfy Tii <0, Sii <0 and Tij r0, Sij r0 for i 6 j. Let (O,F,P) be the probability space on which ransom variables of service times are de®ned. Then, the state process of such a network can be described by a QBD process Y fYt ; tr0g with a ®nite state space: E 0,j: 1RjRm1 [ i,j,k: 1RiRN ÿ 1,1RjRm1 ,1RkRm2 [ N,k : 1RkRm2 :
17
The index i, 1RiRN ÿ 1, denotes the number of customers in server 2; the index j, 1RjRm1 , represents the phase of service in server 1; while the index k, 1RkRm2 indicates the phase of service in server 2. We assume that Y is irreducible. The states are labeled in the lexicographic order, that is, 0,1, . . . , 0,m1 , 1,1,1, . . . , 1,1,m2 , 1,2,1, . . . , 1,2,m2 , . . . , 1,m1 ,m2 , . . . , N,1, . . . , N,m2 . Then, the in®nitesimal generator of the QBD process Y is given by 2 3 T t0 a b 0 0 0 0 0 6 Im1 S 0 A1 7 A0 0 0 0 0 6 7 60 7 A2 A1 A0 0 0 0 6 7 18 A6 7 6 7 40 t0 Im2 5 0 0 0 A2 A1 0 0 0 0 0 0 a s b S where t0 ÿTe, s0 ÿSe; t0 a represents the multiplication of m1 1 matrix t 0 and 1 m1 matrix a, and s0 b is similar. represents the Kronecker product of matrices. The matrices A0 ,A1 ,A2 are given by A0 t0 a Im2 ,
19
A1 T Im2 Im1 S,
20
A2 Im1 s0 b:
21
and
656
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
Since the state space is ®nite, Y is positive recurrent. By the theory of Markov process, Y has a unique steady-state probability vector P P0 ,P1 , . . . ,PNÿ1 ,PN , where P0 is an m1-dimensional row vector; P1 , . . . ,PNÿ1 are all of m1 m2 -dimensional; PN is of m2-dimensional. Let F fn: 0RnRN g be the discrete state space of the QBD process Y, where n represents the number of customers in server 2. Assume that f:F4R is a performance function. The performance measure Zf of the network is de®ned as its expected value with respect to the steady-state measure P, i.e., Zf
N X f np n Pf:
22
n0
where p n jPn j, n 0,1, . . . ,N, with jPn j representing the sum of all components of the vector Pn ; f f0 ,f1 , . . . ,fN 0 is an m1 N ÿ 1m1 m2 m2 -dimensional column vector; fn is the same dimensional vector as Pn , and its components are all f n. Now, we suppose that vectors a,b and matrices T,S are all dierentiable functions of parameter y on the interval J R. Then, the matrix A is also dierentiable on J. We also assume that f depends on y, and is a dierentiable function of y.
5. Some special cases In this section, we discuss sensitivity formulas of the networks with phase-type distributions of special forms, which have explicit physical meaning. 5.1. Phase-type distributions We ®rst consider the models of the phase-type distributions given in the end of section Phase-Type Distributions. We assume that two servers have service distributions of such a form with representation a a1 ,a2 , . . . ,am1 , T diag m11 ,m12 , . . . ,m1m1 Q1 ÿ Im1 and b b1 ,b2 , . . . ,bm2 , S diag m21 ,m22 , . . . ,m2m2 , respectively. Where Ql q li,j , l 1,2. Now, we only consider sensitivities with respect to service rates m1j , j 1,2, . . . ,m1 . Denote by Ti the ith row vector of the matrix T. It is easy to get that 2 3 0 6. 7 6 .. 7 6 7 6 7 60 7 @T 6 7 23 m1i 6 Ti 7, @m1i 6 7 60 7 6 7 6 .. 7 4. 5 0
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
2
3
0 7 6. 7 6 .. 7 6 7 6 0 7 6 @ 0 7 1 6 t a b m1i qi 6 a b 7, m1i 7 6 @m1i 7 60 7 6 7 6 .. 5 4. 0 2
0 6. 6 .. 6 6 60 @A 6 m1i 6 Ti Im2 @m1i 6 60 6 6 .. 4. 0 and
2
657
24
3 7 7 7 7 7 7 7, 7 7 7 7 5
0 6. 6 .. 6 6 60 @A0 6 m1i m1i q i 1 6 a Im2 6 @m1i 60 6 6 .. 4. 0
25
3 7 7 7 7 7 7 7: 7 7 7 7 5
26
where q i 1 1 ÿ
m1 X j1
q i,j1 ,
i 1,2, . . . ,m1 . Pk pk1 ,pk2 , . . . ,pkm1 for k 1,2, . . . ,N ÿ 1 with Let P0 p01 ,p02 , . . . ,p0m1 , pij pij1 ,pij2 , . . . ,pijm2 , j 1,2, . . . ,m1 and PN pN1 ,pN2 , . . . ,pNm2 . From Eqs. (12), (23)±(26), we then have m1i where
Nÿ1 X @Zf @f Pk B ki x f m1i P , i 1,2, . . . ,m1 , @m1i @m1i k0
27
658
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
h
B 0i Ti
m1i q i 1 a b
i 0 ,
0
28
" B ki
#
and
Ti Im2
m1i q i 1 a
Im2
" i B Nÿ1
0
0 , k 1,2, . . . ,N ÿ 2,
29
# Ti Im2
m1i q i 1 Im2
:
30
@Z
Similarly, the expression for m2i @ m f can be derived. 2i
5.2. Coxian distributions For the above phase-type distribution, let a1 1, ai 0, i 2,3, . . . ,k, and qi,i1 ai , qi,j 0, j 6 i 1, i 1,2, . . . ,k ÿ 1, qk 1, then we get the Coxian distribution Ck . It is a phase-type distribution with representation a 1,0, . . . ,0, and 2 3 ÿ1 a1 0 6 7 ÿ1 6 7 6 7 6 7: T diag m1 ,m2 , . . . ,mk 6 7 6 7 4 akÿ1 5 0 ÿ1 Now suppose that the service distributions of two servers are the Coxian distributions Cm1 and Cm2 with representations (a,T) and (b,S), respectively. Where a 1,0, . . . ,0, b 1,0, . . . ,0, and 2 3 ÿ1 a1 0 6 7 ÿ1 6 7 6 7 7, T diag m11 ,m12 , . . . ,m1m1 6 6 7 6 7 4 am1 ÿ1 5 0 ÿ1 2 6 6 6 S diag m21 ,m22 , . . . ,m2m2 6 6 6 4
ÿ1
0 In this case, Eqs. (28)±(30) become
b1 ÿ1
0
3
7 7 7 7: 7 7 bm2 ÿ1 5 ÿ1
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
" B 0i
ÿm1i
m1i 1 ÿ ai
m1i ai
0
659
#
0 ,
"
#
B ki
ÿm1i Im2
m1i ai Im2
m1i 1 ÿ ai Im2
0
0 ,
k 1,2, . . . ,N ÿ 2, and
" i B Nÿ1
#
ÿm1i Im2
m1i ai Im2
m1i 1 ÿ ai Im2 :
where am1 0. 5.3. Erlang distributions In the k-order Coxian distribution Ck , let q1 q2 qkÿ1 0, m1 m2 mk km, we then get the k-order Erlang distribution Ek . Its representation is given by a 1,0, . . . ,0 and 3 2 ÿkm km 0 7 6 ÿkm 7 6 7 6 7 kmJk , T6 7 6 7 6 4 km 5 0 ÿkm where
2
ÿ1 1 6 ÿ1 6 6 Jk 6 6 6 4 0
0
3
7 7 7 7: 7 7 1 5 ÿ1
Now, we suppose that the service distributions of two servers are Erlang distributions Em1 and Em2 . Their representations are given by (a,T) and (b,S), respectively. Where a 1,0, . . . ,0, b 1,0, . . . ,0, T m1 m1 Jm1 and S m2 m2 Jm2 . Then, we have t0 ÿTe 0, . . . ,0,m1 m1 0 . From this, we can obtain m1
@T T m1 m1 Jm1 , @m1
660
2
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
0 6 @ 0 m1 t a b 6 40 @m1 m1 m 1 m1 and
0 0 0
3 0 7 7, 05 0
@A1 T Im2 m1 m1 Jm1 Im2 , @m1 2
0 @A0 6 6 m1 4 0 @m1 m1 m1 Im2
0 0 0
3 0 7 7: 05 0
Substituting the above equations into Eq. (12), we can get the sensitivity formulas for the case of Erlang distribution. 5.4. Hyperexponential distributions For the phase-type distribution of Section 5.1, let qi,j 0, for all i, j, i.e., Q 0, we then get the k-order Hyperexponential distribution Hk . Its representation is given by a a1 ,a2 , . . . ,ak , and T diag ÿm1 , ÿ m2 , . . . , ÿ mk . Now we assume that the service distributions of two servers are Hm1 and Hm2 , whose T diag ÿm11 , ÿ m12 , . . . , ÿ m1m1 , and representations are a a1 ,a2 , . . . ,am1 , b b1 ,b2 , . . . ,bm2 , S diag ÿm21 , ÿ m22 , . . . , ÿ m2m2 . In this case, substituting Ti 1 1 0, . . . ,0, ÿ m1i ,0, . . . ,0, and q 1 1 q2 qm1 1 into Eqs. (28)±(30), we can obtain the i expressions of B k , k 0,1, . . . ,N ÿ 1. Finally the sensitivity formulas can be obtained. 6. An algorithm and an example In this section, we provide an algorithm to compute the potential x f and the performance derivative @Zf =@mij of the network. This algorithm is directly applicable to the controlling and optimization of CQNs, since it is based on analyzing a single sample path. A numerical example is provided to illustrate the application of the algorithm. Let S fj g i infft:tr0,Y fjt g ig be the ®rst passage time from state j to state i for the Markov process Y. By Theorem 2 in Cao and Chen (1997), we have ( ) i S fjg i h ÿ ÿ Zf dt : f Y fjg 31 dij E t 0
It is easy to verify that x f ÿD f P 0
32
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
661
is a potential vector of the network.From Eq. (27), we have Nÿ1 X
@Zf @m1i
k0
Pk B ki x f m1i
P
@f , i 1,2, . . . ,m1 : @m1i
33
6.1. Algorithm 1. 2. 3. 4.
S fj g i Estimate E 0 f Y fjt g dt, ES fj g i and Zf on a single sample path. Calculate dij 's and the realization matrix D f by using Eq. (31). Set x f ÿD f P 0 . Calculate @Zf =@mij by using Eq. (33).
6.2. Example Consider a two-server cyclic queueing networks with N 3, each having the two-order Coxian distribution C2. In this case, F f0,1,2,3g, and there are 12 states in E. The parameters of the network are given by m11 1, m12 2,
m21 2, m22 3, a1 b1 0:5:
The performance function is f n n, for n 2 F; i.e., Zf is the steady-state mean response time of server 2. The means of the potentials for simulation run of 1,000,000 times are as follows: ÿ1:1219 2:1111
ÿ0:3779 0:8534
ÿ0:0014 3:3763
ÿ0:7298 2:1343
1:2028 0:1791 3:9805 2:8676
The theoretical values are given by ÿ1:1089 2:1061
ÿ0:3671 0:8466
0:0038 ÿ0:7246 1:2026 0:1700 3:3668 2:1167 3:9879 2:8588
The mean of the performance derivative @Zf =@m11 for same run is given by 0.8292, while its theoretical value is 0.8259. The standard derivations of our estimates are reasonably small.
7. Conclusions In this paper, with the approach of the in®nitesimal generator perturbation, we studied sensitivities of the steady-state performance for two-server cyclic queueing networks. Each server had the service time of the phase-type distribution. With the potentials of the networks, sensitivity formulas of the performance with respect to the parameters of the systems were given for parameter-dependent performance functions. The main contribution of this paper is that we apply the approach of the in®nitesimal generator perturbation to the queueing
662
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
networks with nonexponential distributions, and provide a method to solve the problems of sensitivity analysis for such a class of networks. Especially, we derived sensitivity formulas of performance with PH-distributions of special forms for two-server cyclic queueing networks. These formulas provide a theoretical basis with the single sample path-based performance sensitivity analysis of CQNs. The system models with PH-distributions of special forms can be directly applied to many practical engineering systems. The accurate estimations of the potentials in these formulas can be obtained from a single sample path, therefore, these formulas can be used directly in optimization problems for these systems. Since a Markov system is the fundamental model for many CQNs, the above approach is expected to apply to a wide class of systems. In future, we shall consider multi-server closed queueing networks with PH-distributions. The key point may be the expression of the in®nitesimal generator of the QBD process describing this system. We shall extend the above results to other system models such as generalized semi-Markov processes (GSMP). The key point of extension is that we need to ®nd a matrix for a GSMP, it corresponds with the in®nitesimal generator in a Markov process. Another future research directions is the application of the potential in stochastic optimization.
Acknowledgements The authors would like to thank the referees for their valuable comments and suggestions which helped to improve the paper.
References Cao, X.R., 1987a. First-order perturbation analysis of a single multi-class ®nite source queue. Performance Evaluation 7, 31±41. Cao, X.R., 1987b. Realization probability in closed Jacksson networks and its application. Advanced in Applied Probability 19, 708±738. Cao, X.R., 1994. Realization Probabilities: The Dynamics of Queueing Systems. New York, Springer-Verlag. Cao, X.R., 1995. An overview of perturbation analysis. In: Proceedings of Chinese Control Conference, pp. 22±39. Cao, X.R., Chen, H.F., 1997. Perturbation realization, potentials, and sensitivity analysis of Markov processes. IEEE Transactions on Automatic Control 42, 1382±1393. Caramanis, M., Liberopoulos, G., 1992. Perturbation analysis for the design of ¯exible manufacturing system ¯ow controllers. Operations Research 40, 1107±1125. Cassandras, C.G., Strickland, S.G., 1989. On-line sensitivity analysis of Markov chains. IEEE Transaction on Automation and Control 34, 76±86. Chong, E.K.P., Ramadge, P.J., 1993. Optimization of queues using in®nitesimal perturbation analysis-based algorithms with general update times. SIAM J. of Control and Optimization 31, 698±732. Glasserman, P., 1991. Gradient Estimates via Perturbation Analysis. Kluwer, Boston. Glasserman, P., Hu, J.Q., Strickland, S.G., 1991. Strongly consistent steady-state derivative estimates. Probability in the Engineering and Informational Science 5, 391±413. Haurie, A., L'Ecuyer, P., Delft, C., 1994. Convergence of stochastic approximation couple with perturbation analysis in a class of manufacturing ¯ow control models. Discrete Event Dynamic Systems: Theory and Application 4, 87±111.
B. Yin et al. / Intl. Trans. in Op. Res. 6 (1999) 649±663
663
Heidelberger, P., Cao, X.R., Zazanis, M.A., Suri, R., 1988. Convergence properties of in®nitesimal perturbation analysis estimates. Management Science 34, 1281±1302. Ho, Y.C., Cao, X.R., Cassandras, C.G., 1983. In®nitesimal and ®nite perturbation analysis for queueing networks. Automatica 19, 439±445. Ho, Y.C., Cao, X.R., 1991. Perturbation Analysis of Discrete-Event Dynamic Systems. Boston, Kluwer. Konstantopoulos, P., Zazanis, M., 1992. Sensitivity analysis for stationary and ergodic queues. Advanced in Applied Probability 24, 738±750. Neuts, M.F., 1981. Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach. The Johns Hopkins University Press, Baltimore. Suri, R., Zazanis, M.A., 1988. Perturbation analysis gives strongly consistent sensitivity estimates for the M/G/1 queue. Management Science 34, 39±64. Tang, Q.Y., Chen, H.F., 1994. Convergence of perturbation analysis based optimization algorithm with ®xed number of customers period. Discrete Event Dynamic Systems: Theory and Application 4, 359±375. Vakili, P., 1991. A standard clock technique for ecient simulation. Operations Research Letters 10, 445±452. Yin, B.Q., Zhou, Y.P., Yang, X.X., Xi, H.S., Sun, D.M., 1996. Sensitivity formulas of performance in closed statedependent queueing networks. In: Proceedings of Chinese Control Conference, pp. 671±677. Yin, B.Q., Zhou, Y.P., Xi, H.S., Sun, D.M., 1997. Sensitivity formulas of performance and realization factors with parameter-dependent performance functions in closed queueing networks. In: Proceedings of The Second Chinese World Congress on Intelligent Control and Intelligent Automation, Xi'an, P. R. China, pp. 1884±1889.