Uniformization for nonhomogeneous Markov chains

Uniformization for nonhomogeneous Markov chains

Operations Research Letters 12 (1992) 283-291 North-Holland November 1992 Uniformization for nonhomogeneous Markov chains N i c o M. v a n D i j k U...

444KB Sizes 1 Downloads 74 Views

Operations Research Letters 12 (1992) 283-291 North-Holland

November 1992

Uniformization for nonhomogeneous Markov chains N i c o M. v a n D i j k Unicersity of Amsterdam, Amsterdam, The Netherlands Received February 1990 Revised July 1992

The discrete Poissonian representation for transition probabilities of homogeneous continuous-time Markov chains, known as uniformization or randomization, is extended to time-inhomogeneous chains. For computational purposes a discrete-time approximation is also provided.

Continuous-time Markov chains; uniformization; time-discretization

1. Introduction

The technique of uniformization, introduced by Jensen [8] and also known as randomization, is used to transform continuous-term Markov chains into discrete-time Markov chains. It is widely known as a powerful tool for evaluating continuous-time Markov chain applications such as those arising in computer performance evaluation, telecommunications and reliability analysis. Especially over the last decade it has been extensively employed for computation and sensitivity analysis of both steady-state and transient measures (e.g. [5-7,10,14,15,19,24]). Next to the essential requirement of uniformly bounded transition rates, the uniformization method seems to have a second restriction. It seems to be limited to time-homogeneous Markov chains. This latter restriction is reasonably justified for steady-state analysis but far less realistic for transient analysis. For example, the availability or reliability of a system which is subject to breakdowns will generally deteriorate by age. A time-inhomogeneous version of uniformization, however, has not been reported. This note provides such a version. Roughly speaking, this version is similar in spirit to the homogeneous case in that it comes down to a random sampling of Poissonian-event epochs at which transitions may take place according to conditional jump probabilities. This result can be of practical interest for the following reasons: (i) It suggests a computational tool of combining iterative numerical computations with random sampling or Monte-Carlo simulation. (ii) Truncation error bounds can be derived more easily when truncation procedures are employed for numerical computation. For the purpose of practical computation, or more precisely of avoiding the necessity of storing a continuum of transition matrices, a finite-grid approximation is also provided along with an error bound for its accuracy. Correspondence to: Professor Nico M. van Dijk, Department of Econometrics, University of Amsterdam, 1018 NB Amsterdam, The Netherlands. 0167-6377/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved

283

Volume 12, Number 5

OPERATIONS RESEARCH LETI'ERS

November 1992

2. Preliminaries Consider a continuous-time inhomogeneous Markov chain (or Markov jump process) with state space S = {1, 2 . . . . } and time-dependent jump or transition rates at time t:

qt(i, J) for a jump or transition from state i into a state j at time t. Here, j = i is included, where for j = i, the jump represents some event at which the state remains unchanged or at which a dummy transition from a state into itself occurs. R e m a r k 2.1. The main reason for including these dummy transitions is purely technical to establish the desired result in the next section. But also from a practical point of view it can be natural. For example, in an M / M / s / s - l o s s system with arrival rate A, in state s in which all servers are busy a ' d u m m y transition' occurs at rate A which represents that a customer arrives but is rejected. Assume that for some /3 < ~, all i and t < Z: (2.1)

qt(i) = Y'~ qt(i, j)
Further, the transition rates are assumed to be continuous in t, for t < Z and all i, j. R e m a r k 2.2. Again for technical reasons later on, we include j = i in the summation of (2.1). The technical L e m m a 2.1 below will remain valid also if we exclude j = i in (2.1). Definition 2.1. A family of transition probabilities {Ps,t 10 < s _< t < Z} is said to be Markov, hereafter called Markov semigroup, if for all i, j and s, t, v with s < s +. t + v < Z:

P,,,+,+.(i,j) = ~ P,,,+,(i,k)Ps+t,s+t+.(k,j).

(2.2)

k~S

Based on this Markov (or semigroup) property, the following Lemma, which is a minor modification of well-known results in the literature in that dummy transitions are included, shows how a unique Markov semigroup can be constructed from the jump rates {qt(i, j)}. Since this construction will be used later on, it is included in the proof. Lemma 2.1. T h e r e exists a unique Markov semigroup {Ps,t 10 _< s < t < Z} such that for all t, t + At < Z and all i, j ~ S:

Pt,t+At( i, j) / A t --->q,( i, j), [Pt,t+at(i, i ) / A t - 1] --* k

j ~ i,

Eq,(i, j),

j=i,

(2.3)

j~i

as At ~ O. H e r e the convergence is to be read in weak convergent sense. That is, the left hand sides in (2.3) need to be uniformly bounded for all t, t + At < Z and all i and i, j. Proof. U n d e r the assumption that qt(i, i) = 0 ~or all t a n d i and with the .convergence in (2.3) uniform in all i, j and t, this can be found at various places in the literature (e.g. [3], T h e o r e m 4 at p. 364). Also the extension to weak convergence is concludable from literature in combination with results for weak infinitesimal generators (e.g. [2]). Let us just give the minor differences due to the inclusion of dummy transitions and the construction of the transition probabilities which will be used later on. As

qt(i,i)-

Eqt( i,j)=j~S

284

~-,qt(i,j), j4-i

(2.4)

Volume 12, Number 5

OPERATIONS RESEARCH LETTERS

November 1992

condition (2.3) is equivalent to:

[Pt,t+,at(i, j)

j)-- l{j=iJ[k~Sqt(i, k ) ]

- l{j=i}](At)-l-'->qt(i,

for all t, t + At < Z and all i, j (j = i included), as At l~m = 1 if event A is satisfied and 0 otherwise. With define transition probabilities by the construction:

(2.5)

~ O. Here l~a I denotes an indicator of event A, i.e. {qt(i)} as per (2.1), that is with j = i included, now

Ps,t( i, J) = E P~t( i, J), k=O

P° is,t,, j ) = l{j=i}exp[- fstqu(i ) d u ] ,

JslZ

Pk+lti,., , , J') = texp -

q,(i) du

(2.6)

]I

]

~., q~.(i, m)PL~t(m, j) dr,

m~S

k_> O.

Then, in analogy with [3, p. 347-353, p. 364-366,] the convergence condition (2.5) and thus (2.3) is verified, also for j = i, and uniqueness of this semigroup is concluded based on its infinitesimal conditions (2.1) or (2.3). As only difference with [3], dominated rather than uniform convergence arguments are to be used. Here in fact, the uniqueness is determined only by the transition rates qt(i, j), j 4: i, for leaving a state, so that the result is not affected by including positive rates qt(i, i). (This later uniqueness can also be concluded using [2, p. 22-34].) Remark 2.3. Theorem 4 on p. 364 of [3] also presents the construction of a Markov process {X, It _> 0} with transition probabilities Ps,~ given by (2.6) and sample paths which are right-continuous. As the finite-dimensional distributions are uniquely determined by (2.6), the corresponding Markov process is unique within the space of right-continuous sample paths, more precisely D[0, ~) (for the precise definition see [1]) by virtue of Theorem 14.5 of Billingsley [1].

3. Uniformization

Now consider the following two Markov chains in the setting of Section 2: (i) A CTMC with transition rates {qt(i, J)l t < Z} as described in Section 2, that is satisfying condition (2.1) for some/3 < ~ and possibly with qt(i, i) > 0. (Note that we thus consider a new CTMC). Let {Ps,t I s, t < Z} be the corresponding Markov semigroup as per construction (2.6). (ii) A CTMC with transition rates J) l t < Z} given by:

{qt(i,

~t(i, j) =/3Pt(i, j),

(3.1)

Pt(i' J) =

(3.2)

where

Y'.qt(i, j)¢1-1 , j = i .

j÷i

In the description of Section 2 that is, whenever the system is in state i, a jump will occur after an exponential time with one and the same parameter /3 independent of that state i as following from:

qt(i) = Y'. qt(i, j) =/3 Y'. Pt(i, j) =/3. j~S

(3.3)

j~S

285

Volume 12, Number 5

OPERATIONS RESEARCH LETTERS

November 1992

Given that this jump occurs at time t, the state will change from state i into state j with probability Pt(i, j). Particularly, a 'dummy transition' into state i will occur with probability P,(i,i). By virtue of Lemma 2.1, where we recall that dummy transitions were allowed, also for this second model (ii) there thus exists a unique corresponding Markov semigroup which satisfies condition (2.3) with qt(i) and qt(i, j) replaced by ~t(i) and ~t(i, j). By replacing q,(i) and q,(i, j) by ~,(i) and ~t(i, j) in (2.6) and substituting (3.1) and (3.3), this semigroup for model (ii) is given by:

Us,,(i, j ) :

E Usk,t(i, j ) , k=0

Us°,(i, :) = l{;=qexp[- ( t - s)[3],

(3.4)

J ) = fs/eXp[-(/: -s)[3][3[ ~mP,(i, m ) U ~ t ( m , j)] dr,

u ki ,+s ,l t(

k>_ 0.

The following theorem will show that the transition semigroups of model (i) and model (ii) are identical and given by a Poissonian power series. Theorem 3.1. For all 0 <_s <_t <_Z:

oo

: k=0

is ,

k!

""

{ t l < t 2 < . . .
where dH(t 1. . . . . t k) is the density of the order statistics t 1 <_t 2 <_ • • • < t k of a k-dimensional uniform distribution at [s, t] × . . . × Is, t] c R k and where PtlPt2 . . . Ptk is the standard matrix product of the transition matrices at times t~ < t 2 <_ • • • <_ t k. Proof. We use the notation o ( A s ) when o ( A s ) / A s -~ 0 as As ~ 0 and o(1) for o(1) -, 0 as As ~ O. Now

first note that: o

Y'. Us ,ks

+ As

<

--

e -("-')t3

[3 e -("-")t3 du

dv < [32

du dv

k=2

(3.6)

= l[3Z[As] 2 = o ( A s ) .

Further, as for all v ~ [s, s + As]: e -'("-S)=[1+o(1)],

and

P,(i,j)=Ps(i,j)[l+o(1)],

while Uslt can be written explicitly as: Usl, t ( i ,

j)

=

fst[3 e-/3("-')P,(i, j )

(3.7)

e -13(t-v) dr,

for any t ~ [s, s + As), we can write: [U~°~+as(i, J) + usl,s+as(i, J)] = 1~_;~[1 - [3 As + o(Zls)] + [3[as][1 + = [[3P,(i, j ) [ A s ] [ 1 + o(1)] ~ l - [ 3 As +[3 A s P s ( i , i ) [ l 286

o(1)]Ps(i,

+o(1)] +o(As),

j ) [ 1 + o(1)] [1 + o(1)]

j-~i, j=i.

(3.8)

Volume 12, Number 5

OPERATIONS RESEARCH LETTERS

November 1992

Substituting (3.2) and letting As ~ 0 thus shows:

qs(i, j), 1 - Eqs(i,

[Us°s+a~(i, j) + U~s+as(i, j ) ] / d s ~

j ~ i,

J), J =i.

(3.9)

j4-i

Combining (3.4), (3.6) and (3.9) verifies condition (2.3) with Us,t for Ps,r By Lemma 2.1 the semigroups Us,, and Ps,, thus coincide for 0 < s _< t < Z. To prove the specific form given in (3.5), we conclude from (3.4) where we restrict to a matrix notation: U " = f t e - I 3 ( V l - S ) / 3 e c evk~t 1 d v 1 S,t

JS

=

~1

I'

fst-fl(vl-s)l'4p[(te-fl(v2-vl)lqpljk-2du2]due ,-

Pl/.J.tU1

P"

t)2v/)2,I

I

=fstf[....fT,._l/3k. o e-~(~.~-s, e-~(v2-~.0 . . . e--~(,k--,k 0 e l~(t--,'O ×P,.P,, "'" P,,~ d r , " - dv 2 dv 1 /3k e_~(t_S) fsstfvt . --ft pviPv2 . . . pvk dVk . . . dv2 dvl 1

=/3k e-~(t-s) (t --S) k

k!

Vk-]

fstfs t

"'"

is tPtIPt2

"'" P'k d / ~ ( t , . . . . . tk),

(3.10)

{s
as the density d H ( - ) of the order statistics X~ _<.g2 -< "'" < X, of k independent uniformly distributed random variables X1, X 2 , . . . , X k at [s, t] is given by k ! / ( t - s)' (e.g. [9, p. 100-103] or [17, p. 240, exercise 29]). [] Remark 3.1. (Interpretation). Clearly, by substituting d H k ( t i . . . . . t k) = k ! / ( t - s ) k, we can rewrite (3.5) by: at)

Ps.,=u,,=

f... f

PtPt2...Ptkdt dt2...dt,.

(3.11)

{S
k=0

However, the form (3.5) is given because of its more natural interpretation as an extension of the standard Poissonian expansion in the homogeneous case (see remark 3.2). Furthermore, it also directly leads to a truncation error bound (see remark 3.3). The interpretation is as follows: During the time interval Is, t] jumps occur according to a Poisson process with parameter/3. Given that k jumps have occurred the jump epochs are uniformly distributed and have a density d H k for the ordered times. At these times the state changes according to the transition probabilities (3.2). Remark 3.2 (Homogeneous case). The standard expression for the homogeneous case: c~

Psi,

=

U~,,

=

E f l k ( t - - s ) k e-~(t-s) pk k=0

(3.12)

k!

is included by setting P = Pt for all t, and noting that:

f,

fs

d H ( q . . . . . tk)

(t

. s). k. .

dt <

(3.13)

{s~tj<_ ...
Volume 12, Number 5

OPERATIONS RESEARCH LETTERS

November 1992

R e m a r k 3.3 ( T r u n c a t i o n error b o u n d ) . With fiL S~t the truncated version of (3.5) at level k = L, and with [[" [I denoting the standard s u p r e m u m norm, one directly derives E k=L+l

flk(t--S)ke-¢('-')< k!

-s)

(3.14)

L!

R e m a r k 3.4 ( C o m p u t a t i o n by simulation). As expression (3.5) involves integrals with respect to the density of a uniform distribution an approximate computation can be obtained by Monte-Carlo simulation as follows. Truncate the poisson distribution at some level k = L as in (3.14) Next compute the Poisson probabilities with p a r a m e t e r [13(t - s)] for any k < L. Then for fixed k do the following: Step 1. G e n e r a t e k random numbers x i ~ [s, t], i = 1 , . . . , k. Step 2. Take their order statistics t I =~1 < t2 =ff2 < " " " < tk =£k" Step 3. Compute (approximate or simulate) the product PtlPt: . . . Ptc Step 4. R e p e a t these steps a large number of times and average.

4. Discrete-time approximation A major complication of computing the non-homogeneous expression (3.5) is the fact that a continuum of transition matrices is required. In this section we study a finite-grid approximation to overcome this complication. Such approximations are well-known in numerical analysis (cf. [13]) and have also been investigated for stochastic processes (cf. [11,12,20,24]). However, as the present application can only be concluded from this literature by combination of several steps, we prefer to give a short selfcontained treatment. For simplicity (see R e m a r k 4.2 below for a possible relaxation) we assume that the transition rates are Lipschitz. That is, for some constant K and all t, t + A t < Z: 1[Qt+At -- Qt 11 -< A t . K ,

(4.1)

where Qt denotes the matrix of transition rates qt(i, j ) for j ~ i and with the diagonal elements given by - q , ( i ) as per (2.1). Here, I1" II denotes the usual s u p r e m u m norm II A II = s u p i Y ' v l a i j l for a matrix A. Let h > 0 be some step-size where we assume h
Lemma 4.1. F o r arbitrary n < m , s = nh a n d t = m h < Z define:

u

=E

[

(t-s)lk k!

"''

m,m e-t~(t-s)

k=0

1

Z

Z

n

n

"'"

m,

_

E

P,,,hP.:h " " P , , ~ h H ( n l , n : . . . . . n k )

,

n

{nl <_nz< - ..- <_nk}

(4.2) where H is the probability m a s s distribution o f the order statistics n I < n 2 <_ • • • <_ n k o f a k - d i m e n s i o n a l u n i f o r m distribution at {n, n + 1 . . . . . m - 1}k. Then, Uh II.,m

(4.3)

.mhll _
Proof. As for all v < Z: (4.4)

PL.=[I+hQ~],

first note that by virtue of (4.1): [[ P,: - P[~,h 'in 288

II = h [I Q~, -

Q[~,h-']h 1[ < h 2 K

(4.5)

Volume 12, Number 5

OPERATIONS RESEARCH LETTERS

November 1992

As a consequence, for arbitrary t I < t 2 < • • • < t~ and by using that transition matrices have an o p e r a t o r n o r m 1 with respect to the s u p r e m u m norm, we obtain by iteratively using (4.5) for v = t~, t 2 . . . . , t~:

Il PttPt2 " " " et~ -- Pn~hPn2h " " " Pnkh II

-IIe,,[P,2

p, -eo2h

-< II P,= " ' " Pt~--Pnzh

e n k h ] ~- [ P t l - P n l h ] [ e n 2

" ' " Pnkh

h " ' " Pnkh][[

II +heK < "'" < kh2K.

(4.6)

Further, with

H ( nl, n 2 , . . . , nk) = k ! / ( m

- n) ~,

(4.7a)

d H ( t l , t 2 .... , tk) = k ! / ( t - s) ~,

(4.7b)

note that for arbitrary function g(t 1, t 2 . . . . . tk):

f, ft.. f,

g ( n , h , n2h . . . . . n~h) d H ( t l , t2 . . . . . tk )

{t~
= ~.,

m--1

~_, "'"

n

n

{hi
...

~_, g ( n l h , n2h . . . . . n k h ) H ( n 1, n 2 .... , n k ) .

(4.8)

n
Applying (4.8) with g(tl, t 2 . . . . . tk) = Pt,Pt2 " "" Ptk, and using (4.6) and (4.7) yields:

ft fst'"

fst P,P,2 " " P, k d H ( t l , t 2 .... ,tk)

{tl<~t2<_ ""


rn--I m - - 1 -

E

E

n

n

~1 '

.....

n

{nl <_n2 < . . .
f f f s t ' ' ' fs t [P',Pt2"'' Ptk-P~,hP~2h ' ' ' Pnkh I d H ( t , , t 2 . . . . . tk)) I

{tl <--t2 < "'" <-0,}


d H ( t l , t2 . . . . . t k ) = k h 2 K "

(4.9)

{t l
Combining (3.5), (4.2) and (4.9) and recalling that h
h z g k=0 ~-~ k e -t~('-s) [/3( t k!- s)] ~

T h e o r e m 3.1 completes the proof.

h(t-s)K.

(4.10)

[]

R e m a r k 4.1. As s t a n d a r d (cf. [20]), u n d e r condition (2.1) it can also be proven that for some constant C:

Ilenh.mh--es.,ll <_hC,

n=[sh-1],

m=[th-l].

L e m m a 4.1 thus leads to an error b o u n d of o r d e r h for the discrete-time approximation applied to arbitrary intervals [s, t]. 289

Volume 12, Number 5

OPERATIONS RESEARCH LETTERS

November 1992

Remark 4.2. The Lipschitz condition (4.1) can be too strong in a practical situation. However, in the proof of Lemma 4.1 we only used that: I1Q,, - Q[t,h-ll h II -< hK for all v < Z. As a consequence, with {Q~.Iv ~ [0, Z]} piecewise constant at intervals [rh, rh + h ) ; r = O, 1,... , [ Z h - l ] , Lemma 4.1 even applies with K = 0. Other relaxations are clearly possible. Remark 4.3. Another possible complication is the size of the state space when very large or infinite. In that case, matrix truncation methods (such as the North-West Corner Rule, e.g. [18]) can be used. Error bounds for the effects of these truncations may be concluded by results from [23]. Remark 4.4. Lemma 4.1 can be seen as a time-inhomogeneous extension of the standard Euler approximation (cf. [13,21]): [I+hQ]"~T

t,

n=[th-'],

for homogeneous semigroups {Ttlt > 0} which satisfy the relation d / d t T~ = Q T t = TtQ. As such, the probabilistic and direct proof by means of uniformization is of interest in itself. Remark 4.5. Similarly to [22] also an approximate uniformized expression of the form (4.2) can be proposed when the boundedness condition (2.1) is violated. Remark 4.6. Recently, in [16] an elegant method is proposed to approximate transition probabilities of homogeneous continuous-time Markov chains by inspecting the process at exponential times. This method resembles the uniformization method but is not the same. It turned out to be amazingly efficient. In [24] further extensions of this method were given. By the results of this paper an extension of this method to nonhomogeneous Markov processes also seems possible.

Acknowledgement I am grateful to professor W. Grassman for addressing the nonhomogeneous case as an open problem during his talk at the first international workshop on numerical solutions of Markov chains, Raleigh, January 8-10, 1990.

References [1] [2] [3] [4] [5] [6] [7]

[8] [9] [10] [11] [12]

290

P. Billingsley, 1968, Convergence of Probability Measures, Wiley, New York, 1968. C.B. Dynkin, Markov Processes I, Springer, Berlin, 1965. I.I. Gihman and A.V. Skorohod, Introduction to the Theory of Random Processes, W.B. Saunders, Philadelphia, PA, 1969. I.I. Gihman and A.V. Skorohod, Controlled Stochastic Processes, Springer, Berlin, 1979. W. Grassman (1990), "Finding transient solutions in Markovian event systems through randomization", in: Proceedings of." "The First International Conference on the Numerical Solution of Markov Chains", January 8-10, Raleigh, NC. 375-385 (1990). D. Gross and D.R. Miller, "The randomization technique as a modeling tool and solution procedure for transient Markov processes", Oper. Res. 32, 343-361 (1984). P. Heidelberger, and A. Goyal, "Sensitivity analysis of continuous time Markov chains using uniformization", in: G. Iazeolla, P.J. Courtois and O.J. Boxma (eds.), Computer Performance and Reliability, Elsevier Science Publishers, Amsterdam, 1988, 93-104. A. JerJsen, "Markov chains as an aid in the study of Markov processes, Skan. Aktuarietisdskr. 36, 87-91 (1953). S. Karlin and H.M. Taylor, A Second Course in Stochastic Processes, Academic Press, New York, 1981. J. Kohlas, Stochastic Methods of Operations Research, Cambridge Press, 1982. T.G. Kurtz, Extensions of Trotter's operator semigroup approximations theorem, J. Funct. Anal. 3, 111-132 (1969). H.J. Kushner, Probability Methods for Approximations in Stochastic Control and for Elliptic Equations, Academic Press, New York, 1977.

Volume 12, Number 5

OPERATIONS RESEARCH LETTERS

November 1992

[13] T. Meis and U. Marcowitz, Numerical Solution of Partial Differential Equations, Springer-Verlag, Berlin, 1981. [14] B. Melamed and M. Yadin, "Randomization procedures in the computation of cumulative-time distributions over discrete state Markov processes", Oper. Res. 32, 926-943 (1984). [15] A. Reibman and K. Trivedi, Numerical transient analysis of Markov models, Comput. Oper. Res. 15, 19-36 (1988). [16] S.M. Ross, "Approximating transition probabilities and mean occupation times in continuous-time Markov chains", Probab. Engrg. Inform. Sciences 1,251-264 (1987). [17] S.M. Ross, A First Course in Probability, 3rd ed., Collier MacMillan, New York, 1988. [18] E. Seneta, Non-negative Matrices and MarkoL~ Chains, Springer-Verlag, New York, 1980. [19] K.S. Trivedi, Probability and Statistics with Reliability, Queueing and Computer Science Applications, Prentice-Hall, Englewood Cliffs, NJ, 1982. [20] N.M. Van Dijk, "Controlled Markov processes: time-discretization", CWI Tract 11, Amsterdam. [21] N.M. Van Dijk, "A note on error bounds for approximating transition probabilities in continuous-time Markov chains". Probab. Engrg. Inform. Sciences 2, 471-474 (1988). [22] N.M. Van Dijk, "Truncation of Markov chains with applications to queueing", Oper. Res. 39, 1018-1026 (1991). [23] N.M. Van Dijk, "Approximation uniformization for continuous-time Markov chains with an application to performability analysis", Stochastic Process AppL 40, 339-357 (1992). [24] B.S. Yoon and J.G. Shanthikumar, (1989), "Bounds and approximations for the transient behavior of continuous-time Markov chains", Probab. Engrg. Inform. Sciences 3, 175-198 (1989).

291