NPS55-82-Q15 WAL POSTGRADUATE SCHOOL Monterey, California FINITE MARKOV CHAIN MODELS SKIP-FREE IN ONE DIRECTION by G. Latouche D. P. Gaver P. A. Jacobs April 1982 FEDDOCS D 208.14/2: NPS-55-82-015 Approved for public release; distribution unlimited. Prepared for: Office of Naval Research n ^'ngton, VA 22217 DUDLEY KNOX LIBRARY NAVAL POSTGRADUATE SCHOOL MONTEREY, CA 93943-5101 NAVAL POSTGRADUATE SCHOOL MONTEREY, CALIFORNIA Rear Admiral J. J. Ekelund David A. Schrady Superintendent Acting Provost This report was prepared with the partial support of the Probability and Statistics Program of the Office of Naval Research, Arlington, VA. Professor D. P. Gaver also gratefully acknowledges the hospitality of the Universite Libre de Bruxelles, and the support of the IBM Company. The work of Associate Professor P. A. Jacobs was partially supported by The National Science Foundation, Grant Number ENG-79-10825. Reproduction of all or part of this report is authorized. UNCLASSIFIED SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) DUDLEY KNOX LIBRARY M^f T L ^° STGRADUATE SCHOOL READ INSTRUCTION^' BEFORE COMPLETING FORM REPORT DOCUMENTATION PAGE 1. REPORT NUMBER NPS55-82-015 2. GOVT ACCESSION NO 3. . RECIPIENT'S CATALOG NUMBER 4. TITLE (and Subtitle) FINITE MARKOV CHAIN MODELS SKIP-FREE IN ONE DIRECTION 5. TYPE OF REPORT ft PERIOD COVERED Technical Report 6. PERFORMING ORG. REPORT NUMBER 7. AUTHORfsJ D. P. Gaver P. A. Jacobs G. Latouche 8. CONTRACT OR GRANT NUMBERf*; 9. PERFORMING ORGANIZATION NAME AND ADDRESS Naval Postgraduate School Monterey, CA 93940 10. PROGRAM ELEMENT. PROJECT, TASK AREA ft WORK UNIT NUMBERS 61153N; RR014-05-GE N0001482WR20017 II. CONTROLLING OFFICE NAME AND ADDRESS Office of Naval Research Arlington, VA 22217 12. REPORT DATE April 1982 13. NUMBER OF PAGES 53 U. MONITORING AGENCY NAME A ADDRESSf// dlflerent from Controlling Oltlce) 15. SECURITY CLASS, (of thla report) UNCLASSIFIED 15a. DECLASSIFl CATION/ DOWN GRADING SCHEDULE 16. DISTRIBUTION ST ATEMEN T (of thla Report) Approved for public release; distribution unlimited 17. DISTRIBUTION STATEMENT (of tha abatract entered In Block 20, It different from Report) 18. SUPPLEMENTARY NOTES 19. KEY WORDS (Continue on tavaraa aide If neceeeary and Identify by block number) 10. ABSTRACT (Continue on reveree aide If neceaaary and Identify by block number) Finite Markov processes are considered, with bi-dimensional state space, such that transitions from state (n,i) to state (m,j) are possible only if m <_ n+1 . The analysis leads to efficient cimputational algorithms, to deter- mine the stationary probability distribution, and moments of first passage times. dd ,; FORM AN 7) 1473 COITION OF I NOV 08 IS OBSOLETE S/N 0102-014-A601 | UNCLASSIFIED SC CURITY CLASSIFICATION OF THIS FACE (When Data Entered) FINITE MARKOV CHAIN MODELS SKIP-FREE IN ONE DIRECTION by G. LATOUCHE P. A. JACOBS and D.P. GAVER Universite Libre de Naval Postgraduate School Bruxelles, Bruxelles Monterey, California Belgium U.S. A Technical Report N° 123 ABSTRACT Finite Markov processes are considered, with bi -dimensional state space, such that transitions from state (n,i) to state (m,j) are possible only if m < n+1. The analysis leads to efficient computational algorithms, to determine the stationary probability distribution, and moments of first passage times. 1. INTRODUCTION We consider irreducible, bi-dimensional Markov processes {(X(t),Y(t)), t > 0>, on a finite state space {(n,i), < n < N, 1 < i < K }, for which transitions from state (n,i) to state (m,j) are possible only if m<n+l. The infinitesimal generator for such a process has a lower block- Hessenberg form : Q ■ '0,0 '1,0 '2,0 '0,1 '1,1 '2,1 1,2 '2,2 2,3 Q N-1,0 Q N-1,1 Q N-1,2 Q N-1,3 Q N,N Q N,1 Q N,2 Q N,3 Qn- Qn-i N-l.N-1 V N-1,N Q N,N-1 Q N,N (1.1) The diagonal blocks Q are square matrices of order K , < n < N. Their diagonal elements are strictly negative, all other elements are non- negative. The blocks Q , n * m, are rectangular with appropriate dimen- sions; their entries are non-negative. The rowsums of Q are equal to zero. We denote by level n the set {(n,i), 1 < i < K } of states corres- ponding to the common value n for the first index. The structure (1.1) of Q permits the Markov process to move up by only one level at a time, and to move down by several levels at a time. The process is said to be skip-free to the right. 3. Many models in queueing theory lead to Markov processes (or Markov chains) for which the infinitesimal generator (respectively the transition probability matrix) is a (possibly infinite) block-Hessenberg matrix. The paradigm for the lower-Hessenberg case is the GI/M/1 queue, for the upper-Hessenberg case, it is the M/G/l queue. Neuts [7,9] has studied in great detail Markov processes and Markov chains with infinite block-Hessenberg matrices which exhibit a repetitive structure : for some n*, K n = K n* * Q„ m = Qr,_Li ^1 » f° r all n > n*, m < n+1, ^n,m ^n+l,m+l if the matrix is lower-Hessenberg, and a similar condition if the matrix is upper-Hessenberg. In particular, Neuts has shown that if a process with such a lower block-Hessenberg matrix is positive recurrent, its stationary probability vector has a modified matrix-geometric form. For processes with upper block-Hessenberg matrix, the quantities which may be most elegantly studied are first passage times to lower levels; these first passage times correspond to the busy period in queueing theory. There do not exist many published results for processes with finite block-Hessenberg matrices. Most of them deal with block-Jacobi matrices (the blocks Q are zero if either m > n+1 or m < n-1). Torrez [12] proposes to determine the stationary distribution for birth-and-death processes in a random environment, by using numerical procedures designed for band matrices. Neuts [10] and Carrol et al . [1] explicitely determine the stationary probability distribution for the PH/M/C and M/PH/1 queues with finite waiting room. Hajek [3] considers a model with repetitive structure : Q = Q, ,, 1 < n < N-1, Q ,i = i)i i and Q„, , = Q~ , , 1 < n < N-2. Processes with general block- ^n,n+l 1,2 ^n+l,n v,l 3 Jacobi matrices are analysed in Keilson et al . [4], and in Gaver et al . [2], where specific numerical examples are also given. Wikarski [13] determines, by purely analytic arguments, the stationary probability vector for processes such that K Q > K, > ... > K N , and such that the matrices Q n n+1 are of rank K- + i» < n < N-l. We do not need these restrictions. Moreover, our approach yields quantities which have a clear probabilistic interpretation. The purpose of the present paper is to show how the stationary probabi- lity distribution, and moments of first passage times from some level to another level, may be efficiently computed. The stationary probability distribution is analysed in the next section. In Section 3, we analyse the first passage time from level n to a higher level m, <n < m <N. The results in these two sections are related to those in [2]. The analysis of the first passage time from level n to a lower level k, < k < n < N, is more involved, and is presented in Section 4. We shall observe that it is possible to compute independently the stationary probability distribution, and the moments of first passage times to higher levels. On the other hand, it is necessary to compute both the stationary distribution, and moments of first passage times to higher levels, in order to compute the mo- ments of first passage times to lower levels. In Section 6, we give some numerical examples of common-cause of failure models. We consider systems with N machines and R repairmen. If n denotes the number of machines in working condition, the repair rate is ymin(R, N-n). The failures occur in this way : in an interval (t, t+dt), there is a jump from n to n-l with probability An dt + o(dt) (independent failures) or to n-k, 1 <k < n, with probability v(£)q (l-q) n ~ dt + o(dt) (common-cause failures). The results to follow are applicable, mutatis mutandis, to Markov pro- cesses with upper block-Hessenberg infinitesimal generator, and to finite Markov chains with block-Hessenberg transition probability matrix, as we briefly mention in Section 5. 5. Notational conventions . We define a number of vectors and matrices with indices. Occasionally, we refer to specific elements of those. In order to somewhat simplify the notations, we use Q n i(i»j) to represent the (i,j)th element of the matrix Vn" We systematically use "n" to denote a given level, "m" to denote a higher level, and "k" to denote a lower level. 2. THE STATIONARY PROBABILITY DISTRIBUTION We denote by p_ the vector of stationary probabilities associated with the Markov process Q : p_Q = , p_e = 1, (2.1) where e is a column vector with unit elements. We partition the vector p as p_ = (p~, p_i, ..., Pm)» where the subvectors p_^ have K elements and correspond to the states of level n, < n < N. In order to determine the vector p_, we need three preliminary lemmas. Lemma 1 Consider a Markov process on the state space {1,2,. . . ,r,r+l,r+2,. . . ,r+s}, with infinitesimal generator Q : Q = A A where A is a square matrix of order r, A is a rectangular r by s matrix. The states r+1 to r+s are all absorbing. a) The states 1 to r are all transient if and only if the matrix A is non- singular, in which case (-A~ ) > 0. b) The (i,j)th entry of (-A" ), for 1 < i , j < r, is the expected amount of time spent in the transient state j, starting from the transient state i, before absorption in any of the absorbing states. c) The (i,k)th entry of (-A~ A), for 1 < i < r, r+1 < k < r+s, is the probabi- lity that, starting from the transient state i, absorption occurs in the state k. d) The (i,k)th entry of (-A'^f-A" 1 A) , for 1 < i < r, r+1 < k < r+s, is the expectation of the time spent in the transient states, starting from the transient state i, before absorption, measured on those paths that lead to absorption in the state k. Proof . These results are well known, although we are not aware of publications where they are proved under the stated form. The first statement is proved in Neuts [9], Lemma 2.2.1; the second statement is a consequence of Neuts and Meier [11], Corollary 2; the third statement is proved in [2], Lemma 1. To prove the fourth statement, we must show that the (i,k)th entry of (-A )• (-A A) is equal to C t d v- k^' wnere v- k (t) = P [the process gets absorbed in k, before time t | the initial state is i]. Conditioning on the state of the process at time t > 0, one obtains, for h > v ijk (t+h) - (1 +Aii h)v ifk (t) + A ijk h + I Aij h v j>k (t) + o(h), hence v!^(t) = I A^ v j>k (t) +A i>k . or v'(t) = A v(t) + A , where v(t) is the matrix formed by the v i k (t) 's. We already know that v(0) = 0, and lim v(t) = (-A )A, therefore, we easily obtain that t-» • J~ t d v(t) = - [t((-A _1 )A - v(t))]~ + Jq ((-A _1 )A " V(t))dt ■ JJ ("A" 1 ) v'(t) dt - ("A" 1 ) [v(t)]J = (-AJ'^-A^A) We now define the processes S , < n < N. For 1 < n < N, S is the restriction of the original process Q, observed during those intervals of time spent at level n, before the process Q enters level n-1 for the first time. The state space of S is {(n,j), < j <K }. Clearly, all S , 1 < n < N, are transient Markov processes. The process Sq is the res- triction of the process Q observed at the lowest level 0; it is an ergodic Markov process. We denote by C the infinitesimal generator of the process S , < n <N. We also define the matrices n , , < k < n < N as follows : n n i,(i »j) is the probability that, starting from state (n,i), the first state visited at any level below level n is the state (k,j). Lemma 2 The matrices C , < n < N are recursively determined as follows C N = Q N> , . (2.2) e n - 1n.« + On.n+1 n n + l,n . < n < N-1. (2.3) Proof . The equation (2.2) is obvious : starting from level N, the process S N terminates as soon as the process Q enters any of the levels to N-1. Meanwhile, the transitions are governed by C\. ... Consider n fixed, 1 <n <N-1, and let Z (t), t > be defined as follows. Z n ("0 = i if the process S is in state (n,i) at time t. Assume that Z (x) = i and Z (x + dx) = j * i. Since j * i, a transition must have occured in the process Q. Z (x + dx) = j if and only if one of the following events have occured. a. The transition is from (n,i) to (n,j), this happens with probability Q n>n dx + o(dx). b. The transition is from (n,i) to (n+l,h) for some h, and the process Q returns to (n,j), after spending an unspecified amount of time at level n+1 or higher, before visiting any other state at level n or lower. This happens with probability Q n n+l (1»h) dx n n+1 n (h,j) + o(dx). It is now clear that P[Z n (x + dx) = j|Z n (x) = i] K n + 1 " Vn^'J) dT + h ^ ^n.n+l^ 1 ^) n n + l,n( h 'J) dT + o(dx), for 1 < i * j < K n , and similarly, that P[Z n (x + dx) = i|Z n (x) = i] K n + 1 - (1 ♦ Q Bi „<l.l) ^) ^ Q n>n+1 (i,h) n n+ltB (h.1) dt + o(dx), for 1< i < K n , which proves (2.3) for n > 1. The proof for n = is identical. Lemma 3 The matrices TT .,0<k<n<N are recursively determined as n j k follows Vk ■ (^N 1 ) Q N,k • forO<k<N-l, (2.4) "n.k •<<>««.* ♦«n.iHlW>' for < k < n * N-l (2.5) 9. Proof . The equation (2.4) is an immediate consequence of Lemma 1 and (2.2) A simple probabilistic argument leads to Vk = K!n> tQ n,k + Q n,n + 1 n n + l,k + Q n,n + 1 n n + l,n n n,k ] « which, together with (2.3), implies (2.5). □ Remark . The matrices C , < n < N and n.,0<k<n<N, must be numerically computed together, as indicated below. Algorithm A A.l C N := Q NjN ; A. 2 for k from to N-l do determine rL. . from (2.4) ; A. 3 for n from N-l to 1 step -1 do begin determine C from (2.3); for k from to n-l do determine rr . from (2.5) — — n,k v ' end ; A. 4 determine £ Q from (2.3)« We are now in a position to prove our first main result. Theorem 1. Th e vectors p. < n < N, are determined by the equations 2q C Q = 0, (2.6) fin-fin.lVl.nK 1 )' ft»rl<n<M» (2.7) N I fi.e=l. (2.8) n=0 ~" ~ 10. Proof . From the structure (1.1) of Q, it results that the system p Q = may be decomposed into N £0 %,0 + ^ fiv ^v,0 = °- ( 2 - 9 ) N ^-l\-l,k + Pk\,k + v ! k+1 ^ Q v,k = ^ forl<k<N-l, (2.10) fiN-1 Q N-1,N + ^N,N = °' (2.11) From (2.11) and Lemma 2, p^ = -p^ Q N _ M Q"j N = p^ Q N _ M (-Q-J M ), which proves (2.7) for n = N. The remainder of the proof uses an induction argument in two steps. Induction assumption 1 : for a given k, 1 < k < N-l, the equation (2.7) holds for k+1 < n < N. In order to prove that (2.7) holds for n = k, we first need to show that Z , £v ^,k " fit <- C t> n t,k ' for k+l< t < N. (2.12) Clearly, (2.12) holds for t = N, by (2.4). An induction argument will be used to prove (2.12) for arbitrary t > k+1. Induction assumption 2 : for a given a, k+1 < a+1 <N, the equation (2.12) holds for a+1 < t < N. We prove that it holds for t = a, by using successively the assump- tion 2, the assumption 1, and (2.5). N * J^A.k =£a Q a,k + Ha+l <" C a+l> "a+l.k ' v=a = K ^a,k + ( 5a,a+i n a+l,k)' 11. Finally, to show that (2.7) is satisfied for n = k, note that (2.10) can be written as N ° = J>k-i Q k-i,k + Mk.k + v * k+1 £v Q v,k = fik-1 Q k-l,k + h Q k,k + Hk+1 ( _e k + l) n k+l,k = 2k-lQk-l,k + £k (\,k + \,k + l n k + l,k)- The equation (2.6) is proved in the same way. Since C Q is the infinitesimal generator of a finite, irreducible Markov process, the equation (2.6) has a unique solution, up to a multiplicative constant, and that constant is determined by (2.8). a We have thus shown that the matrices C , < n < N, are the crucial quantities needed to determine the stationary probability distribution, We shall observe in Section 4 that the same matrices, together with the matrices n k , < k < n < N, play a central role in the evaluation of first passage times from a given level n to a lower level k. We must, however, first examine the first passage times to higher levels. 3. FIRST PASSAGE TIME TO HIGHER LEVELS We denote by i n n i the first passage time from level n to a different level n' : t^, = inf {t > : X(t) = n'|X(0) = n}, < n,n' <N, n * n'. We are concerned in this section with the first passage times t for J n,m m > n. We need to first establish some preliminary results We define the processes S n , < n < N in a similar fashion to the processes S n . For < n <N-1, the transient Markov process S is the restric- tion of the process Q, observed during those intervals of time spent at level n, before the process Q moves upward to level n+1 for the first time. The ergodic process S N is the restriction of the process Q observed at the highest level N. We denote by C n the infinitesimal generator of the process S , < n < N. 12, We also define the matrices r nm , < n < m < N as follows : r n m (i>j) is the probability that, starting from state (n,i), the first state visited at level m is the state (m,j). More formally, r n m (i»J) » p £ Y ( T n rJ = J I x (°) = n, Y(0) = 1]. Clearly, r e = e, for all n and m. (3.1) Lemma 4 . The matrices C , < n < N, are recursively determined as follows c o - V • < 3 - 2 > C n = "n,n + [fo "n.k r k,n- 'or 1< n < N . (3.3) a The proof is analogous to that of Lemma 2 and is omitted here. Lemma 5 follows The matrices r „,0<n<m<N are recursively determined as n,m r n,n + l " K 1 ) Vn + 1 • *r <n <M-1. ( 3 -4) r n m = r n m . r . (3.5) n,m n,m-l m-l,m ' = n r f , + , for < n < N-2, t=n+l w,t n+2 <m <N. (3.6) Proof . For n=0, the equation (3.4) is an immediate consequence of Lemma 1 and (3.2). For 1 <n < N-l, observing the system at the first departure from level n leads to -1 n_1 r n,n+l = ( " Q n,n )(Q n,n+l + k J Q Q n,k r k,n r n,n+l)» which, together with (3.3), implies (3.4). To prove (3.5), since Q has the structure (1.1), one merely decompose the probabilities in r according to 13. the first state visited at level m-1. Equation (3.6) is now obvious. The matrices C„, < n < N and r , < n < m < N, must be numeri n n,m cally computed together, as indicated below. Algorithm B B.l Cq := Qq^q ; determine Tq , from (3.4); B.2 for m from 1 to N-l do begin determine C from (3.3) ; determine r, ._, from (3.4); m,m+i v ' for n from to m-1 do determine r ml1 from (3.5) — — n,m+l v ' end ; B.3 determine C» from (3.3). We now define, for x > 0, < n,n' < N, n* n', 1 < i < K n , Kj<K n . , 9n,i;n',j< x > - P[T n,n' <x « Y < T n,n' + °) = ^^ = n ' Y (°) " l3 ' < 3 ' 7 ) we denote by G n .(0 the matrix such that G „ • (C) (1 » J ) equals the Laplace- Stieltjes transform of the possibly defective distribution g n . , .-(.). n ,i ,n ,j Clearly, since the Markov process may only move up by one level at a time, we have that G n >m («) " G „,n,-1<«> Vl. m <«) < 3 - 8 > m = n G t , t (0, forO<n<N-2, (3.9) t=n+l w »* n+2 < m < N. 14, (3.10) A simple argument conditioning on the state of the process at time h > shows that Vi;n + l,j( x+h ) = [1 + «n,n>U hl 9 n,i ;n + l,j( x ) + v J i (Q n.n)l,v h 9n,v;n + l,jW n-1 K k K n + k=0 v =l (Qn ' k)i * v H y=l 9k » v ' n ^ * 9n ^ n+1 'J (X) + (Qn,n + l>i,j h+0 < h >* forl<n<N-l, where * denotes the Stleltjes convolution, from which we conclude that 5 G n,n + 1< 5 > * «n,n G n,n+1< 5 > + «n.n+l + 3o Q ".l' G k.n (C)G n,n + l(«)' fo, - 1<n<N - 1 - and similarly, * 6 o,i^) ' ^0,0 %.i««) + %,v < 3 - u > This equations may be written as Wl< ?) • D n<5) "n,n + l • *r < n < N-1. (3.12) where D„(?) = (El - Q^,,)" 1 . (3.13) ° n M - (CI - Q n , n - V Q n>k \Jt))' 1 . for l< n < N-1. (3.14) Of course, we have G„ (0) = r m , and from (3.3) D n (0) = (-C' ), n,m v ' n,m v ' n x ' x n ' < n < m < N. 15. Let U n,n' -~ G "."' U) I** " (3 - 15) and Hn n ' = U n n' - * for < n,n' < N, n * n' . (3.16) We have U„ m (1»j) = E[x n m , Y(x n m + 0) = j|X(0) = n, Y(0) = i], n,m v J/ n,m v n,m ' -l v ' ■ x ' and H„, m (i) - «x I X(0) - n. Y(0) - 11. Theorem 2 . The expected values of the first passage times to higher levels are given by the following relations. "o.i ■ (-Co 1 ) £' < 3 - 17 > Hn,n + 1 " <<><& + "=! V Hfc.n*' for l * n < "- 1 ' (3 ' 18 > k=0 i^ m = u« m i + r nm1 u mlm ,forO<n< N-2, (3.19) — n,m — n,m-l n,m-l -m-l,m v n+2 < m < N . Proof . The equation (3.17) immediately results from Lemma 1. The equations (3.19) are directly proved by using the strong Markov property : W i)=E[ Vml X<0) - n. Y<0) - 1] = E [T n,m-l I X (°) = n ' Y <°> = i] + E[E[ Vl.ml Y ( T n,m-l + °)ll X (°) " "• Y < ) =i] ' Conditioning on the first level visited after level n, and using Lemma 1, we obtain 16, Vn + l = K!n>^ + [ z =0 ("O Vk Vn+1 = ( " Q "* n) i+ k=o (_g "' n) Qn « k (j ^ n + Fk » n ^* n+l) by (3.19). Hence, n-1 n-1 (Q n,n + k l Q Q n,k r k,n) Vn + 1 = " <• + £, Vk M|c.n>' which, together with (3.3), proves (3.18). Remarks . a. Another argument, using the equations (3.12), (3.15) and (3.16), requires the evaluation of the matrix 3 D (0/35 • This is done in Appendix A. Similarly, higher moments of first passage times may be obtained by purely routine, albeit tedious, manipulations. b. The vectors u^ may be computed together with the matrices r , by appropriately modifying Algorithm B. FIRST PASSAGE TIME TO LOWER LEVELS We denote by t k the first passage time from level n to a lower level k. In the preceding section, we have in essence decomposed the first passage time from level n to a higher level m, into a sum of first passage times from n to n+1, from n+1 to n+2,..., from m-1 to m. Since it is possible to move down several levels in one transition, no such simple decomposition is possible for t ■. Rather, we must explicitely consider the first level visited below n, starting from level n. In order to do tnis, we define the first passage time T n as follows T n = inf {t > : X(t) < n | X(0) = n} , = min {x n ^ k , < k < n-1}, for 1 < n < N. 17. We define, forx>0,0<k<n<N, l<i^K n , 1 < j < K k , g n>i;kjj (x) = P [T n <x, X(T n + 0) - k, Y(T n + 0) = J|X(0) = n, Y(0) = 1] we denote by G .(C)(i »j) the Laplace-Stieltjes transform of the possibly defective distribution g ... • (.), and b y G n u(0 tne matrix with entries given by S n>k (e)(1 »j). An argument similar to the one used to obtain (3.10) and (3.11) shows that * S N,k (s) = Vk + %,H § N,k< 5) • ^r < k < N-l (4.1) and * B n,k^> = Vk + Vn 5 n,k< 5 > < 4 ' 2 ) + <Wl C W (5) + W (0 W^' for < k < n < N-l. These equations may be written as 5 N,k^) = D N^ Q N,k * for 0<k<N-l. (4.3) B n,k^) = D n<0 W„.k + Q n,n + 1 W (5)] ' *>r < k < n < H-l. (4.4) where B N U) = (?I - Q N>N ) _1 , (4.5) *nM '- <* " Vn " <Wl W^^' for l < n < H ' 1 ' (4 ' 6) Obviously, we have that 5 . (0) = n . , and (Lemma 2) D n (0) = (-C" 1 ). Observe that n u e <e, for < k < n < N, n-l I n . e = e, for 1 < n < N. (4.7; k=0 n,K " ~ 18. Let D <^ = " IT 6 ".* (5) '5=0- (4 - 8) D nfk (i,j) - E[T n ,X(T n + 0) = k, Y(T n + 0) = j|X(0) = n, Y(0) = i]. Lemma 6 . The matrices k are recursively determined as follows °N,k = ^N^ n N,k • forO<k<N-l f (4.9) °n,k = K^Vk + Wl (0 n+l,k + °n + l,n n n,k> ] ' < 4 - 10 > for <k < n <N-1. Proof . It results from Lemma 1 that D N.k " (-tfH-tf «N.k>- forO<k<N-l, which together with (2.4) proves (4.9). For n <N-1, we condition on the first level visited after level n, and we obtain by using Lemma 1 and the strong Markov property that °n,k " Kln> Vk + Kin <Wl )(0 n + l,k + Vl.n n n,k + n n + l,n °n,k>* hence «n.n + Q n,n + 1 n n + l,n> °n,k * ' [n n,k + *n,n + l <°n + l,k + °n + l,n n n,k> ] ' which, together with (2.3), proves (4.10). a Corollary 1 . The vectors u„ defined as — n u. are recursively determined as follows % - K^ + Vn+1 W • for 1 < n < N-l. (4.13) u*. = ("Cm 1 ) £ , (4.12) a This corollary immediately results from Lemma 6 and (4.7). The Laplace-Stieltjes transforms of the distribution functions for the first passage times t . are given by the matrices G ,,(€). ^ one decomposes the trajectories from level n to level k, according to the first level v visited below n, one immediately obtains G n.k<5> ■ G n,k<«> + = S „,v<«> G v,k<«) v=U n-l + E (4.14) -1,4.1 g n (0 G .(c), for <k < n < N, v=k+l n,v x ' v,k v ' from which the moments of first passage times to specific lower levels may be derived. The first moments may be more easily obtained by using the strong Markov property. The next theorem is stated without proof. Theorem 3 . The expected values of the first passage times to lower levels are given by the following relations Hl.0"Hi- < 4 - 15 ) k-1 n-l iJn k = i, + E n n u , + I n u . , for < k < n < N, (4.16) — n,K — n _q n,v —v,k -k+i n » v — v»k v ' where the matrices n are determined by Lemma 3, and the vectors u . n,v J — v,k for v < k-1 are determined by theorem 2. a 20, We observe here a major difference between the first passage times to higher and to lower levels. We have observed in the preceding section (Remark b. after Theorem 2) that the vectors u for m > n may be computed in one loop, starting from n = 0, to n = N-l. The evaluation of the vectors u . for k < n requires two steps. One for computing the vectors u , starting from n = N, to n = 1. In the second step are computed the vectors u . , starting from n = 1, to n = N. We observe also that the vectors u^ for m > n have to be computed, even if one is only interested in the first passage times to lower levels. It is clear that one can combine the computation of the stationary probability distribution, and the expected first passage times, as is done in the algorithm C given in Appendix B. It is much easier to obtain moments for the first passage time from level n to any level at or below k < n-l. Let T n k = inf {t > : X(t) < k|X(0) = n}, = min {t : < v < k} , for < k < n < N. n,v We denote by G . (£) the matrices with entries equal to the Laplace-Stieltjes n jK transforms of the distribution functions for the first passage times T k> Also, let and u„ . = U . e , for < k < n < N. — n,k n,k — If we condition on the first level v visited below n, we obtain *n k< ? ) = * 5 (O + "i G n (0 G k (0. n ' K v =0 ' v=k+l ' ' from which the moments of T . may be obtained. Again, it is easier to deter- n j k mine the first moments by using the strong Markov property. The next theorem 21. is stated without proof. Theorem 4 . The expected values of the first passage time from level n to any level at or below k, are given by u„ „ i = iL , (4.17) — n,n-l — n v ' n-1 and Hn k = "n + z n n x, u x, k » for < k < n < N, (4.18) "»* " \i=k + 1 "» v v,l\. where u is defined in (4.11). d — n x # The vectors li . may be computed in one loop, beginning with k = N-1, n = N, and progressively decreasing k until k = 0; for each value of k, one computes the vectors Um k to u. + , . . In Section 5, we shall refer to Algorithm D (in Appendix B), which is designed to recursively compute the vectors lL . , N-1 > k > 0. 5. GENERAL REMARKS Remark a. We firstly comment on Equation (2.7). The (i,j)th entry of the matrix Q , (-C~ ) is equal to 1 K " 1 n L,n - -(-C^Xk.JHS.l) ' k=1 ("Vl.n-lH 1 ' 1 ) 22, The factor Q n _ 1>n (i » k y("Q n -l n-1^ 1 fi ) is the P robabilit y that » upon leaving the state (n-1,1), the process Q moves to the state (n,k). The (k,j)th entry of (-C~ ) is the expected time spent by the process in state (n,j), starting from (n,k), before hitting any state at a lower level, by Lemma 1 and the definition of £ . If that lower level is below level n-1, the process will be forced, by the structure (1.1) of Q, to pass through level n-1 before reaching again level n. Therefore, it results from (5.1) that [Q , _(-C~ )](i,j) is equal to (-Q j n _i)(i »i)times the expected time spent in the state (n,j), before the next return to level n-1, given that the process Q starts in the state (n-1,1). This is exactly the interpretation of Neuts' matrix R ([9], Theorem 1.7.2, and [8]), and shows that (2.7) is the counterpart to the matrix-geometric solutions. Remark b . If the matrix Q is upper block-Hessenberg, instead of having the form (1.1), then our results can be modified in an obvious way. In that case, the first passage times to lower levels may be directly computed, via an algorithm derived from Algorithm B, while the first passage times to higher levels are more involved. Also, the stationary probability vector is computed by determi- ning first the vector pj, (up to a multiplicative constant) and then recursively the vectors p^j, p^_ 2 j)q. This last observation points to the reason why it is difficult to eva- luate the stationary probability distribution for infinite upper block-Hessenberg matri ces . We briefly examine now the problem of obtaining a good approximation for the stationary distribution in such a case. Neuts [7] and Lucantoni and Neuts [6] consider infinite upper block-Hessenberg matrices of the following form : 23. Q = B B l B 2 b 3 B A l A 2 A 3 A Q A x A 2 A o A i and stationary probability vector x = ()u, x,, x,, ...). An explicit form is found in [6,7] for x« and x, , as well as complicated but tractable expres- sions for m, = I v x e , v=0 m = E v x e. 2 -n -v — v=(J From these, a heuristic method is proposed to truncate the matrix Q at some level N, and it is suggested that x,, x 3 , ..., Xm be determined by a block Gauss-Seidel iterative procedure. An alternative method would be to determine, as indicated in Section 1, the vector £ corresponding to the finite matrix, and to compare the vectors (pn»Pi) to (x_q,Xi), in order to check whether the level N is properly chosen. Also, we may truncate the matrix Q on the basis of first passage times argument, instead of using m, and nu. It is yery easy to adapt Algorithm D (in Appendix B), in order to compute successively the expected first passage time from level to any level at or above m, for m = 1,2,... . One may then choose N to be the smallest value of m for which the expected first passage time exceeds a given, large value. When N is determined, the vector p_ can be readily obtained, since most of the preliminary computations are already performed. 24. Remark c . We have not considered first passage times to specific states, but have concentrated our attention on first passage times to levels. The reason is that, in a number of applications, such as chains in a random environment, the levels themselves are of importance. Nevertheless, it is not very diffi- cult to extend our analysis to the study of x . . = inf{t > : X(t) = m, n , i ,m , j Y(t) = j | X(0) = n, Y(0) = i}. With the results of Sections 3 and 4, we can observe the process between successive visits to level m, until, for the first time, the state (m,j) is visited before the process leaves level m. Remark d. We have only considered continuous-parameter Markov processes. Our analysis may be immediately adapted to discrete-time Markov chains. The major difference is that, in all our equations, the matrices (-C~ ) and (-C~ ) will be replaced by fundamental matrices for absorbing Markov chains (Kemeny and Snell [5], Section 3.2). NUMERICAL EXAMPLES We have computed some characteristics for different machine-repairmen- like processes, using the concepts of common-cause failure (described in the introduction) and of randomly varying environments. We have considered four dif- ferent types of systems, described below in increasing order of their complexity. We have computed, for each system, the stationary distribution p and the vectors iL k - To do this, we have used Algorithm E (in Appendix B), which inclu- des parts of both Algorithms C and D. 25, We shall not attempt to answer any specific question about such systems. Rather, we shall describe the qualitative and quantitative diffe- rences in their behaviour. We shall use here, as the main performance measure, the expected times ijL . until there remain k or less operating machines, start- ing from a fully operational system. Type I . The first model is the classical machine-repairmen process. We con- sider a system of N machines and R repairmen, with failure and repair rates respectively given by x' and y. The infinitesimal generator is in this case a tri -diagonal matrix. Type II . This is a simple common-cause failure model. In superposition to a machine-repairmen system, with N machines, R repairmen, failure and repair rates given by x and y, we consider an independent point process with exponen- tial interevent time distribution, with parameter v : the process of common- cause failures. Whenever a common-cause failure occurs, each operating machine has a probability q of failing, independently of the other machines. Thus, if there are n operating machines, the number of machines that fail due to common- cause has a binomial distribution, with parameters n and q. The infinitesimal generator has the form (1.1) » with blocks of order 1. We define as system I I. A the common-cause failure system with N = 10, R = 1, X = .005, y = 1, v = .01, q = .2. If 1/y is equal to one hour, then 1/X «8 days, and 1/v »4 days. The system I. A is a machine-repairmen system with N = 10, R = 1, X' = .007, y = 1. That value for x' is chosen so that both systems are globally similar in the following sense : when n machines are ope- rating, they both have the same expected number of failures in an interval of length dt : this is n(X dt + q v dt) + o(dt) = .007 n dt for system II. A, and n X' dt + o(dt) = .007 n dt for system I. A. 26. We give in Table I the stationary probability distribution for these two systems- We observe that both systems have high, nearly equal, probabili- ties of having most machines in operation. Nevertheless, the two systems are different since the distribution II. A has a longer tail towards the states with few operating machines. This difference is strikingly apparent on Table II, where we give the expected times u 1Q k , < k < 9. For instance, if we compare u, Q 5 for both systems, we observe that the average time until 5 machines or more are under 3 repair, is smaller by a factor 10 for the common-cause failure system. We conclude that multiple failures have a profound effect on the dyna- mic behaviour of the system. In support of our conclusion that multiple failures are the main cause for difference, we observe that both systems have identical repair facilities (R = 1, y = 1), and that the stationary expected number of failures per unit of time (henceforth denote by cp) are nearly equal : they are respectively given by 10 ip T .= I n p n X' = .0694815, 1,A n=0 n 10 and ip n ,= Z n p(\ + vq) = .0693535. n. a n=0 n Type III. We assume now that the failure and repair parameters change according to an independent environmental process : they are respectively given by \. t u . , J J v. and q. when the environment state it j, 1 < j < K. We assume that the envi- J J ronmental process is an irreducible Markov process, with infinitesimal generator T. The system is now described by a pair of indices (n,j), where n denotes the number of operating machines, and j denotes the state of the environment. The infinitesimal generator has the form (1.1), with blocks of order K. A typical example, with N = 4, R = 2 and K = 3 is displayed in Appendix C. Note that the 27. off-diagonal blocks are diagonal matrices because the only possible transi- tions are from (n,j) to (n,j'), with j * j', or (n,j) to (n',j), with n * n ' . We firstly define a system with two environments. In environment 1 (the off-environment), no common-cause failure may occur. In environment 2 (the on-environment), common-cause failure may occur at the rate v*. Indivi- dual failures occur with the same rate X in both environments. In order to make meaningful comparisons with the system I I. A, we have chosen v* so that the expected time between two common-cause failures is equal to 1/v = (.01) , the same as for the system 1 1. A. The time between two succes- sive common-cause failure has a PH-distribution (Neuts [9], chapter 2), with representation (£, V), where 6 = [1 0], V = (T 21+ v*) T 21 12 -T 12 and its first moment is given by -£ V~ e = (Toi+ Ti 2 )(T, 2 v *)~ » hence T 21 U« = V ( 1 * ) . T 12 We have chosen T, 2 = .00015 and T 21 = .004, from which we obtain v* = .27666. If the unit of time is one hour, then l/Tjp « 9 months, 1/Toi « 10 days, 1/v* » 3.5 hours. We now define the systems III. A and III.B. For system III. A, we have N = 10, R = 1, K = 2, \ l = x 2 = .005, »i = v z = 1 * v i = °» v 2 = - 27666 » ^ = °» q 2 = .2. For system III.B, there is only one difference : u 2 = 2, the service rate is doubled in the on-environment. 28. In Table III, for each system, the first column represents tne stationary marginal distribution of the number of operating machines; the last row represents the stationary marginal probability of being in each environment; in the remaining two columns are given the stationary conditional distributions of the number of operating machines, given that the system is in the corresponding environment. The comparison of the column III.A.O in Table III and the column II. A in Table I shows that the two distributions are >jery similar, but the marginal distribution for the system in a random environment has a longer tail towards the states with few operating machines. We also observe that the marginal distribution is not a good descriptor for the system in a random environment : the two conditional distributions are quite different, which may have important implications, even if there is a small probability of being in the on-environment, In Table IV are given, for the same two systems, the values of u,^ . (1) and u, Q k (2), for < k < 9. The last row indicates the conditional probability of being in each environment, given that there are 10 machines in operation. The columns III.A.l and III. A. 2 are to be compared to the column II. A in Table II. Several differences may be observed. Firstly, the time u, Q ~ to com- 3 plete failure is shorter by a factor 10 for the system in a random environment. Secondly, the values of u, Q ,(1) to u,g 7 (1) are close to each other : they differ at most by a factor 2. This is due to the influence of the on-environment -1 3 it takes on the average (T,-) » 6.6 10 units of time to switch from environ- ment 1 to environment 2, at which time the common-cause failure process is turned on and brings more rapidly the system towards the states with few opera- ting machines. As u\ Q 6 (1) is so close to (T 12 )" » we have computed u 1Q 6 (1) + u 6 k (2), for < k < 5, and found that it differs from the corresponding u^ k (l) by at most 70 units of time. 29. We conclude that the systems 1 1. A and 1 1 1. A have markedly different dynamic behaviours, even when they have identical repair facilities and nearly equal stationary failure rates : 10 *III.A = n % n tp n.l x l + Pn,2 < x 2 + V^ 1 = - 06664 35, is close to <Pjt .. When examining the data for system III .B , we observe that the same comparisons can be made between III.B and 1 1. A as between 1 1 1. A and 1 1. A. 10 The stationary average service rate for system III.B is I (p n , m + P 9 u 9 ) = 1.036145, and ip m B = .0680697. The remaining examples will be derived from the system III.B. We shall next consider systems with three environments. Environment 1 is still the off-environment where no common-cause failure may occur. Environ- ments 2 and 3 together constitute the on-environment where common-cause failure may occur at a constant rate v*, but with different failure probabilities, which we denote respectively by f 2 and f 3 , in order to distinguish from the o probability q 2 in the preceding model. The infinitesimal generator T for envi- ronmental changes is chosen as follows : T = 12 21 21 '12 -(T 21 + T 23 ) T 32 T 23 -(T 21 + T 32 ) (6.1) thus the on-environment has exponential duration with parameter T 2 ,, just as for the systems I I I. A and III.B. We shall set T 32 > T 23 and f 2 < f 3 , so that the on-periods are comprised of intervals of relatively low failure probability f 2 , interrupted by shorter intervals of higher failure probability f 3 . 30. In order to make meaningful comparisons with the systems I I. A and III .B, we have set, as before, v* = v (1 + T 21 /T 12 ), so that the expected time between successive common-cause failures is the same for all systems. We also choose f 2 and f 3 such that q 2 = a f 2 + (1-a) f y where a = Ogl + T 32^ T 21 + T 32 + T 23^" » and is the conditional probability of being in environment 2, given that the system is either in environment 2 or 3. This choice of f 2 and f 3 ensures that, on the average, the probability that any machine fails when common-cause fai- lure occurs is the same for all systems. The system III .C is defined as follows : N = 10, R = 1, K = 3, X, = x 2 ■ X 3 = .005, u 1 = 1, m 2 = y 3 = 2, v x = 0, v 2 = v 3 = .27666, q 1 = 0, f 2 = .1504, f 3 = >36 » T 12 = - 00015 * T 21 = - 004 ' T 23 = '° 4, T 32 = ,125 ' For svstem m - D » the differences are as follows : f 2 = .1801, f 3 = .7, T 32 = 1. Thus, we have for system III.D yery short bursts of high failure probability. The stationary probability distributions are respectively given in Tables V and VI : in column the marginal distribution of the number of ope- rating machines; in columns 1 to 3 the conditional distributions, given the environment state. The last row of each table gives the marginal distribution of the environment state. We compare Tables V and VI with the columns III .B of Table III. The same observation that were made earlier are valid here : the marginal distribu- tion has a greater probability mass at the states with low operating machines. The probabilities corresponding to 6 to 10 operating machines given the off- environment (columns III.B.l, III.C.l, III.D.l) are nearly equal in all three systems. The conditional distributions given the on-environment are noticeably different for all systems; in particular, the distribution in column III.D. 3 is different from all the others, in that it is more or less uniform for n between 2 and 8. 31. In Table VII are given, for the same two systems, the values of u 10 k^) ' ^ or 1^J^3, 0<k<9. The last row indicates the conditional probability of being in each environment, given that there are 10 machines in operation. We observe that short bursts of higher failure probability are more detrimental to the performances of the system, since the values of u 10 . (j) for the system III.C are nearly consistantly greater than those for the system III.D, the exception being u, Q . (3) for high values of k. That exception is probably due to the very short duration of the environment 3. Finally, we have computed that ipjrj - = .0679624, and cpj , . Q = .0680954. Our general conclusion is as follows : if one considers the systems III. A (or III.B), II. A and I. A as increasingly simplified approximations of the system III.C (or the system III.D), we obtain increasingly optimistic measures of their performance, based on the expected first passage times u,q ■ , although the systems have nearly equal stationary expected number of failures per unit of time. Type IV . For our last numerical examples, we shall assume that the server does not change instantaneously the service rate when there is a change of environ- ment. More precisely, we assume that the environment evolves as in the system III.C and III.D, with environmental changes governed by a matrix T of the form (6.1). The "normal" service rate in the off-environment (environment 1) is equal to y,. When the system enters the on-environment (environment states 2 and 3), the server notices it when 2 or more machines fail simultaneously, or after a random interval of time of exponentially distributed duration, with parameter £,. The server then changes the service rate to y 2 - The service rate remains equal to p 2 until an interval of time has elapsed without any failure, that interval of time is exponentially distributed, with parameter £ 2 - In this case, we need to describe the state of the system by three indices (n,j,k), where n denotes the number of machines that are operating, j denotes the environment state, and k = 1 or 2 to indicate the service rate that 32. is in effect. For each given value of n, we order the states as follows (n.1,1), (n,2,l), (n,3,l), (n,l,2), (n,2,2), (n,3,2). The possible tran- sitions from (n,j,k) to (n,j',k') are indicated below, together with the corresponding transition rates. (n,2,l) (n.1.1) *- '2,1 2,3 (n,3,l) (n,l,2) The infinitesimal generator Q has the form (1.1), with blocks of order 6 given by yi .... Hi . ^2 • ^2 'n,n+l , 0<n<N-l, 33. 'n.n '1,2 * T * T '2,1 '2,3 T T * '2,1 '3,2 1,2 T 2,l * 2,3 T T * '2,1 '3,2 , < n <N, 'n,n-l nx 1 ' q n,l * • q ( 3 ) • • • nx, . • M n,l q (3 ! M n,l , 1 <n <N, *n,n-k ,(2) W q( 2 ) q n,k * «$ , 2 <n <N, 2 < k < n, where q$ = nx 1 + v^J) f* (l-f/-". 34. the entries represented by a "." are equal to zero, the entries represented by a "*" are such that the row-sums of Q are equal to zero. We define the following type IV system : N = 10, R = 1, K = 3, Xj = A 2 = X 3 = .005, uj = 1, y 2 = 2, vj = 0, v 2 = v 3 = v* = .27666, q. = 0, f 2 = .1801, f 3 = .7, T 12 = .00015, T 21 = .004, T^ = .04, T 32 = 1, q= ^ = v*/2 = .13833. The column of Table VIII gives the marginal distribution of the number of operating machines, in the remaining columns are given the condi- tional distributions, given (j,k) where j is the environment state, and k the index of the service rate. The last row gives the marginal stationary proba- bilities of the pairs (j,k). In Table IX are given the values of u\ Q (j,k), for 1 < j < 3, k = 1,2 , <n <9. The last row indicates the conditional probabilities of being in the state (10,j,k), given that there are 10 operating machines. We observe from the last row of Table VIII that the states in {(n,l,2), < n < 10} are essentially states of transitions towards the set {(n,l,l), <n < 10}. For both j = 2 and 3, the sets {(n,j,l), < n < 10} and {(n,j,2), < n < 10} have similar probabilities, but different distributions : the states in {(n,j,l), <n < 10} are more concentrated than those in {(n,j,2), < n < 10} towards large values of n, in other words, towards few machines under repair. This is seemingly in contradiction with the inequality y, < u 2 ; this is due to the fact that transitions from (j,l) to (j,2) occur whenever there are multiple failures. The expected first passage times are not much different from those of the system III.D (see Table VII). 35 Finally, we mention that the program implementing Algorithm E has been written in FORTRAN 77 and run on a CDC 170/750. Care has been taken not to waste memory space, by making sure that no space is reserved for matrices n. . with j > i, Q. . with j > n+1, or vectors u. . with j > i. l ,j i ,j i ,j In other respects, the implementation was straight forward. For systems with N = 10, the program uses on the average .2 seconds of central processor time for K = 3, and .6 seconds for K = 6. 7. ACKNOWLEDGEMENTS Work of Donald P. Gaver partly supported by a contract with the Office of Naval Research (NR 609-006). Work of P. A. Jacobs partly supported by National Science Foundation Grant Number ENG 79-10825. 36. REFERENCES [1] J.L. Carroll, A. Van de Liefvoort and L. Lipsky. Another look at the M/G/l and GI/M/1 queues. Operations Research, to appear. [2] D.P. Gaver, P. A. Jacobs and G. Latouche. Finite birth-and-death models in randomly changing environments. Tech. Report N°121. Laboratoire d' Informatique Theorique, Universite Libre de Bruxelles, Brussels, 1981. [3] B. Hajek. Birth-and-death processes on the integers with phases and general boun- daries. Journal of Applied Probability , to appear. [4] J. Keilson, U. Sumita and M. Zachman. Row-continuous finite Markov chains, structure and algorithms. Technical Report N°8115. Graduate School of Management, University of Rochester, (1981). [5] J.G. Kemeny and J.L. Snell. Finite Markov chains . Van Nostrand, Princeton, (1960). [6] D.M. Lucantoni and M.F. Neuts. Numerical methods for a class of Markov chains arising in queueing theory. Tech. Report 78/10, Applied Mathematics Institute, University of Delaware, Newark, 1978. [7] M.F. Neuts. Queues solvable without Rouche's theorem. Operations Research, 27 , (1979), 767-781. [8] M.F. Neuts. The probabilistic significance of the rate matrix in matrix-geometric invariant vectors. Journal of Applied Probability, 17 , (1980), 291-296. [9] M.F. Neuts. Matrix-Geometric Solutions in Stochastic Models - An Algorithmic Approach . The Johns Hopkins University Press, Baltimore, 1981. [10] M.F. Neuts. Explicit steady-state solutions to some elementary queueing models. Operations Research , to appear. 37. [11] M.F. Neuts and K.S. Meier. On the use of phase type distributions in reliability modelling of systems with a small number of components. OR Spektrum, 2 , (1981), 227-234. [12] W.C. Torrez. Calculating extinction probabilities for the birth-and-death chain in a random environment. Journal of Applied Probability, 16 , (1979), 709-720. [13] D. Wikarski. An algorithm for the solution of linear equation systems with block structure. Elektronische Informationsverarbeitung und Kybernetik, 16 , (1980), 615-620. 38. I. A II. A 1 2 .000002 3 .000017 4 .000100 5 .000441 6 .000011 .001472 7 .000230 .003731 8 .004104 .008752 9 .065136 .054839 10 .930519 .930647 Table I Stationary distribution for systems I. A and II. A Missing numbers are less than 10" . 39. I. A II .A 1 2 3 4 5 6 7 8 9 923 10 12 7 10 12 105 10 9 2 10 9 63 10 6 2 10 6 93 10 3 5 10 3 256.9 14.3 106 10 6 5 10 6 463 10 3 60 10 3 10 10 3 2 10 3 710.9 279.5 117.6 17.0 Table II Expected time, starting from 10, until there remain k or less operating machines, for systems I. A and II. A. 40. Ill .A.O III.A.l 1 1 1 . A . 2 III .B.O III.B.l III.B.2 .000005 .000141 .000010 1 .000037 .001015 .000003 .000094 2 .000141 .000001 .003871 .000018 .000487 3 .000382 .000002 .010500 .000065 .001815 4 .000827 •• .000006 .022717 .000189 .000001 .005199 5 .001524 .000013 .041815 .000466 .000003 .000466 6 .002480 .000027 .067911 .001002 .000010 .027451 7 .003667 .000125 .098122 .001923 .000102 .050497 8 .006577 .002196 .123402 .004864 .002165 .076818 9 .051019 .047572 .142937 .049665 .047543 .106236 10 .933342 .950059 .487564 .941804 .950174 .718606 Env . .963855 .036145 .963855 .036145 Table III Stationary distributions for systems III. A and III.B Missing numbers are less than 10" . 41 III.A.l III. A. 2 III.B.l III.B.2 1 2 3 4 5 6 7 8 9 218 10 3 39 10 3 16 10 3 10 10 3 8 10 3 8 10 3 7 10 3 4 10 3 454.6 20.0 211 10 3 32 10 3 9 10 3 4 10 3 2 10 3 841.7 443.9 163.7 15.1 3.6 2 10 6 166 10 3 39 10 3 16 10 3 10 10 3 8 10 3 7 10 3 4 10 3 454.7 20.0 1 10 6 159 10 3 32 10 3 9 10 3 3 10 3 1 10 3 548.8 177.4 15.4 3.6 .981118 .018882 .972421 .027579 Table IV Expected time, starting from (10, j), until there remain k or less operating machines, for systems III. A and III.B. 42, III. CO III.C.l III.C.2 III.C.3 .000004 .000014 .000468 1 .000023 .000087 .002443 2 .000071 .000348 .007189 3 -.000165 .000001 .001080 .015699 4 .000328 .000003 .002939 .028580 5 .000594 .000005 .007332 .045250 6 .001016 .000013 .017046 .062420 7 .001731 .000102 .036057 .074554 8 .004557 .002165 .065256 .078259 9 .049505 .047542 .0106656 .086266 10 Env. .942004 .950169 .763186 .598871 .963855 .027590 .008555 Table V Stationary distributions for system III.C Missing numbers are less than 10~ . 43. Ill .D.O III.D.l III .D.2 III.D.3 1 2 3 4 5 6 7 8 9 10 .000018 .000054 .000105 -.000175 .000284 .000495 .000921 .001739 .004677 .049594 .941938 .000001 .000002 .000003 .000005 .000012 .000101 .002165 .047542 .950169 .000193 .000716 .001700 .003398 .006519 .012724 .024937 .045872 .072797 .105755 .725389 .008376 .020724 .033000 .040077 .039427 .034300 .030762 .033505 .043182 .067661 .648986 Env. .963855 .034760 .001385 Table VI Stationary distributions for system III.D Missing numbers are less than 10 44, III.C.l III.C.2 III .C.3 III .D.l III .D.2 III.D.3 1 2 3 4 5 6 7 8 9 131 10 3 32 10 3 16 10 3 11 10 3 9 10 3 8 10 3 7 10 3 4 10 3 454.9 20.0 124 10 3 25 10 3 9 10 3 4 10 3 2 10 3 1 10 3 665.5 230.2 18.8 3.9 124 10 3 25 10 3 8 10 3 4 10 3 2 10 3 831.7 . 399.0 138.7 13.4 3.5 37 10 3 17 10 3 12 10 3 10 10 3 9 10 3 8 10 3 7 10 3 4 10 3 454.7 20.0 30 10 3 11 10 3 6 10 3 3 10 3 2 10 3 1 10 3 592.8 197.2 16.6 3.7 30 10 3 10 10 3 5 10 3 3 10 3 2 10 3 956.4 486.8 168.6 15.2 3.6 .972209 .022352 .005439 .972277 .026769 .000954 Table VII Expected time, starting from (10, j), until there remain k or less operating machines, for systems III.C and III.D 45, (1.1) (2,1) (3,1) (1.2) (2,2) (3,2) .000024 .000252 .004878 .000054 .000320 .013208 1 .000071 .001010 .015020 .000200 .001152 .029737 2 .000142 .002395 .027803 .000473 .002744 .043229 3 .000240 .000001 .004412 .036023 .000929 .005660 .049054 4 .000393 .000002 .007597 .035707 .001699 .011110 .046435 5 .000684 .000003 .014368 .031340 .003036 .020984 .040339 6 .001246 .000009 .030240 .031798 .005317 .037149 .037239 7 .002253 .000096 .061311 .044074 .008901 .059562 .040634 8 .005373 .002157 .104263 .068355 .014212 .082957 .049765 9 .050324 .047533 .145735 .107482 .041137 .111313 .071410 10 .939250 .950199 .628417 .597519 .924043 .667052 .578950 Env. .963268 .015221 .000580 .000588 .019539 .000804 Table VIII Stationary distribution for a Type IV system. Missing numbers are less than 10" . 46. (1.1) (2,1) (3,1) (1.2) (2,2) (3,2) 34 10 3 28 10 3 27 10 3 34 10 3 28 10 3 27 10 3 1 17 10 3 10 10 3 9 10 3 17 10 3 10 10 3 9 10 3 2 12 10 3 5 10 3 5 10 3 12 10 3 5 10 3 5 10 3 3 10 10 3 3 10 3 3 10 3 10 10 3 3 10 3 3 10 3 4 9 10 3 2 10 3 2 10 3 9 10 3 2 10 3 2 10 3 5 8 10 3 1 10 3 839.8 8 10 3 1 10 3 851.9 6 7 10 3 534.6 441.8 7 10 3 547.2 450.0 7 4 10 3 187.0 160.8 4 10 3 190.7 163.3 8 454.7 16.3 15.0 457.7 16.5 15.1 9 20.0 3.7 3.6 20.0 3.7 3.6 .974497 .010184 .000369 .000578 .013876 .000496 Table IX Expected time, starting from (10,i,k), until there remain n or less operating machines, for a type IV system. A.l APPENDIX A First order derivatives of the matrices D„(0 rr ' The equations (3.13) and (3.14) may be written as D K)( V - "0,0' " '• D n<5>( « -'Q„,„ " Jo Q n.k G k,n<«» " »• 1<n< "-*• By differentiating both sides of these equations with respect to 5, we obtain that DJ(0( 51 " Q , ) + D o^) = °» W«>( 5I " Q n,n " £ Qn.k G k,n<«>) + D n^) < ! " £ Q n .k G k,n<*» = °- 1 < n < N-l, or D£(e) = " (D fl (e)) 2 • (A.l) D;(0 - - D n (0(I - V Q njk G kfn (0) D n (0, for 1< n < N-l. (A.2) First order derivatives of the matrices D (£) The equations (4.5) and (4.6) may be written as D N (5)(EI - Q N , N ) ■ I. D n(0(d - Q„,„ - Q„,„ + i 5„ + i,„(0) - I. Kn<N-l. By differentiating both sides of these equations with respect to £, we obtain d n(0(«i - q n>n ) + D N (e) = 0, A. 2 or D'(c) = - (B N (€)r. 'N (A. 3) and B'(0(5I - Q n , n - Q n>n+1 5 n+1>n (0) ♦ D n (Od - Q n , n+1 6 n+1>n (0) - 0, or D;(0 « - D n U)d - Q n , n+1 8; +1 , n (0) D n (0. for 1 < n < H-l. (A.4) B.l APPENDIX B Algorithm C . Evaluation of the stationary probability vector p_ = (Pn>- • • >Pm) » and the expected first passage times ir , , for all pairs of levels n * n'. CI (Initial computations : the matrices n . , C and the vectors Oj . C N := Q N,N ; determine u». from (4.12); for k from to N-l do determine n,. . from (2.4); for n from N-l to 1 step -1 do begin determine C from (2.3); determine u from (4.13); for k from to n-l do determine n u from (2.5) end ; determi ne C Q from (2.3) j C2 (Determination of the vectors jr and u ,) Hq := solution of the system ttq C"q = 0, tTq e_ = 1; C := Q 0,0 ; Hi.o := ^1 ; for n from 1 to N do determine u^ Q from (4.16) ; determine u Q , from (3.17); determine r Q , from (3.4); for m from 1 to N-l do B.2 begi n determine it from (2.7); m re-normalize the vectors n n to n so that I n. e = 1: ~° Hn 1-0 - 1 " for n from m+1 to N do determine u from (4.16): — n,m v '' determine C from (3.3); determine u a1 from (3.18); m,m+l v ' determine r mJ , from (3.4); m,m+l x ' for n from to m-1 do begin determine u , from (3.19); determine r mJ _. from (3.5) n,m+l v ' end end ; determine it m from (2.7); n N renormalize the vectors n n to n M so that I n„ e = 1. "° -N n=0 -" " At the end of the algorithm, the vectors n are equal to the vectors p. n < N. We sys' overflow problems. <n <N. We systematically renormalize the vectors n in order to avoid Algorithm D. Recursive evaluation of the vectors Vj. defined by v. = u.. . , <k <N-1. Dl (Initial computations) £ N := Q N,N ; determine u N from (4.12); v^_ 1 := Ufti (equation (4.17)) B.3 D2 (Main loop, at the k th step, v^_, is determined) for k from N-l to 1 step -1 do begin for n from N to k+1 step -1 do determine TT . from (2.4) and (2.5) determine C. from (2.3); determine u. from (4.13); \ k-1 - := ^c' ( eQ . uation ( 4 - 17 )) for n from k+1 to N do determine ti . , from (4.18) ; end. Algorithm E . Evaluation of the stationary probability vector p_ = (Pq»...»p_m)» and the expected first passage times iT . , for < k < n < N. This algorithm ~™l I y l\ combines Algorithm D and part of Algorithm C. El (Initial computations) C N := Q N,N ; determine tL from (4.12); for k from to N-l do determine n., k from (2.4) ; E2 (Determination of the vectors if ,, and preparation to the evaluation of the ^1 y l\ vector p_) Hn n-1 := ^W ' ( ec l uation ( 4 - 17 )) for n from N-l to 1 step -1 do begin determine C from (2.3); determine u from (4.13); Vn-1 := % ; ( e q uation ( 4 - 17 )) for m from n+1 to N do determine u , from (4.18); — m,n-l x ' for k from to n-l do determine n ., from (2.5) end; 3.4 E3 (Determination of the vector p_) ttq := solution of the system ttq Cq ■ 0, tt, e_ ■ 1; for m from 1 to N-l do begin determine n from (2.7); m renormalize the vectors n n to n m so that I Tt. e = 1 -0 -m i=Q -i - end ; determine n^'from (2.7); N renormalize the vectors n n to n M so that I n. e = 1. "° -N i=0 - 1 " (The vectors jr are equal to the vectros t) , < n < N) C.l APPENDIX C Infinitesimal generator for Type III . We assume N = 4, R = 2 and K = 3. Entries represented by a "." are equal to zero, entries represented by a "*" are equal to minus the sum of the off-diagonal entries on the same row. Q = * T T '1,2 2,3 2y, T * T '2,1 '2,3 2^2 T T * '3,1 3,2 * 2u 3 q l,l * T T 1,2 1 ,3 2 "i • • • 43 • T * T 2,1 '2,3 • 2Vn • • • 43 T T * '3,1 3,2 • • 2u 3 43 • • 43 • • * T l,2 T l,3 2y. (2) q 2,2 • 43 • T 2,l * T 2,3 2u« (3) * * q 2,2 • • 43 T 3,l T 3,2 * • • 2 "3 qd) • • q 3,3 43 • • q (1) q 3,l • • * T T 1 1,2 1,3 llj • • • 43 ■ • 43 • • q( 2 ) q 3,l • T 2,l * T 2,3 • p 2 ' ■ ■ 43 • • 43 • • q {3) q 3,l T T „ * 3,1 3,2 * u 3 H 4,4 q 4,3 q (1) q 4,2 • • 43 • • * T T 1 1,2 1,3 • 43 ■ • 43 • • q(2) q 4,2 • • 43 • T * T '2,1 '2,3 ■ ■ 43 • • 43 • • q(3) q 4,2 • • 43 T T * 3,1 3,2 where q (k) _ n,i nA k + v k (J) q\ (l-q k ) n "\ 1< k < 3, 1< n < 4, 1 < i < n. DISTRIBUTION LIST NO. OF COPIES Library, Code 0142 4 Naval Postgraduate School Monterey, CA 93940 Dean of Research 1 Code 01 2A Naval Postgraduate School Monterey, CA 93940 Library, Code 55 2 Naval Postgraduate School Monterey, CA 93940 Professor D. P. Gaver 183 Code 55Gv Naval Postgraduate School Monterey, CA 93940 Defense Technical Information Center 2 ATTN: DTIC-DDR Cameron Station Alexandria, Virginia 22314 DUDLEY KNOX LIBRARY - RESEARCH REPORTS 5 6853 01068007 7