TWO-TIME-SCALE MARKOV CHAINS AND ... - CiteSeerX
c 2005 Society for Industrial and Applied Mathematics
SIAM J. APPL. MATH. Vol. 65, No. 2, pp. 567–586
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS TO QUASI-BIRTH-DEATH QUEUES∗ G. YIN† AND HANQIN ZHANG‡ Abstract. Aiming at reduction of complexity, this work is concerned with two-time-scale Markov chains and applications to quasi-birth-death queues. Asymptotic expansions of probability vectors are constructed and justified. Lumping all states of the Markov chain in each subspace into a single state, an aggregated process is shown to converge to a continuous-time Markov chain whose generator is an average with respect to the stationary measures. Then a suitably scaled sequence is shown to converge to a switching diffusion process. Extensions of the results are presented together with examples of quasi-birth-death queues. Key words. Markov chain, singular perturbation, countable state space, asymptotic expansion, occupation measure, aggregation, switching diffusion, quasi-birth-death queue AMS subject classifications. 34E05, 60J27, 60F05 DOI. 10.1137/S003613990139756X
1. Introduction. Much effort has gone into evaluating performance measures of queueing systems in the past two decades; see, for example, [3, 22, 24, 26] and the references therein. Since exact solutions are difficult to obtain in many queueing problems, we are content with approximate solutions. To treat large-scale Markovian queueing systems, one of the most popular methods is decomposition, which consists of breaking the underlying network into smaller pieces (e.g., one station in each piece); see [3, 24, 26] among others. Although one often uses time-homogeneous Markovian models for approximating the actual systems, many queueing systems in real life are nonstationary (time-dependent); for example, the arrival and service rates in the systems are time-varying. Recently, time-inhomogeneous Markovian queueing networks have been widely used to model telecommunication systems; see [5]. Developing computational methods and approximation techniques for these quantities involved in time-inhomogeneous queueing problems has long been regarded as a challenging task (see [14] and [20]); see also [12, 18] and the references therein for earlier effort in this direction. The motivation of our study stems from the recent advances in understanding asymptotic properties of singularly perturbed Markov chains aimed at reducing of complexity. We began our investigation in [13]. By combining matched asymptotic expansions and stochastic analysis, our effort was subsequently extended to treat more complex models and applications in control and optimization; see [29, 30]. These results are about Markov chains with finite-state spaces, and they are applicable to queues with a finite number of waiting rooms or finite capacity; see also [20] for a related work. Continuing our effort initiated in [28], using a model similar to that ∗ Received by the editors November 6, 2001; accepted for publication (in revised form) May 21, 2004; published electronically January 5, 2005. http://www.siam.org/journals/siap/65-2/39756.html † Department of Mathematics, Wayne State University, Detroit, MI 48202 (gyin@math. wayne.edu). The research of this author was supported in part by National Science Foundation grant DMS-0304928 and in part by the Wayne State University Research Enhancement Program. ‡ Institute of Applied Mathematics, Academy of Mathematics and Systems Sciences, Academia Sinica, Beijing, 100080, China ([email protected]). The research of this author was supported in part by a Distinguished Young Investigator grant from the National Natural Sciences Foundation of China and a grant from the Hundred Talents Program of the Chinese Academy of Sciences.
567
568
G. YIN AND HANQIN ZHANG
presented in [26] as our starting point, we consider queueing systems in which the generator of the queue length process includes both a fast varying part and a slowly changing part reflecting both strong and weak interactions of states belonging to different irreducible classes. Compared with the methods used in [3, 24, 26], we adopt a singular perturbation approach via time-scale separation to model the different transition intensities. In [28], treating countable-state-space Markov chains, one of the main assumptions used is that the underlying Markov chain is irreducible. In the context of queueing theory, it basically corresponds to the consideration of a single queue. The goal of this paper is to treat multistation queueing networks and to develop decomposition and aggregation methods for queueing network problems. Motivated by queueing networks involving quasi-birth-death processes, we start with decomposition by splitting the entire state space into a number of subspaces. Usually, the transitions within each subspace are much more intensive and frequent than those among different subspaces. We introduce a small parameter ε > 0 to highlight the different rates of transition intensities. Then we proceed to derive the asymptotic expansions of the probability distribution of the queue length. By lumping all the states in each subspace into a single state, we obtain a sequence of aggregated processes. We demonstrate that this sequence converges to a Markov chain. Furthermore, we obtain limit results for suitably scaled sequences. The rest of the paper is arranged as follows. Section 2 begins with the precise formulation of the problem. Section 3 provides asymptotic expansions. A sequence of occupation measures is defined in section 4; its probabilistic properties, including mean square estimate, aggregation, and switching diffusion limit of a suitably scaled sequence, are examined. Section 5 presents extensions of results and queueing examples. Finally, an appendix is provided to include some technical complements. 2. Formulation. Working with a finite time horizon t ∈ [0, T ] for some T > 0, our focus is on time-inhomogeneous Markov chains. Suppose that β(t) is a continuoustime Markov chain with countable state space N = {1, 2, . . .}. An infinite dimensional matrix-valued function Q(t) = (q ij (t)) defined on [0, T ] is a generator of the Markov chain β(t) if q ij (·) is Borel measurable and bounded for each i, j ∈ N, q ij (t) ≥ 0 for ∞ all i = j, j=1 q ij (t) = 0 for all i ∈ N, and for any bounded and Borel-measurable t function g(·) defined on N, g(β(t)) − 0 Q(s) g (·)(β(s))ds is a martingale, where (2.1)
Q(t) g (·)(i) =
∞
q ij (t) g (j) for each i ∈ N.
j=1
Throughout the paper, we use K to denote a generic positive constant. The conventions K + K = K and KK = K are used for simplicity. We use 11m0 and 11 to denote an m0 -dimensional and infinite dimensional column vector with all components being 1. For a vector z and a matrix H, we use z and H to denote their transposes, and we use z i and hij to denote the ith component of z and the ijth entry of H = (hij ), respectively. For a given matrix H = (hij )∞×∞ with infinite columns and infinite · rows, Ha is an augmented matrix given by Ha = (H : 11). In addition, we use a subscript to index a sequence. A Markov chain β(t) or its generator Q(t) is weakly irreducible, if the system · of equations g(t)Qa (t) = (0 : 1) has a unique solution g(t) = (g i (t)) with g i (t) ≥ 0 for each i ∈ N, where 0 = (0, 0, . . .) is an infinite dimensional 0 vector. The unique nonnegative solution is termed a quasi-stationary distribution. (An equivalent way to write the system of equations in the above definition is g(t)Q(t) = 0, g(t)11 = 1.) The
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
569
definition is an extension of the usual notion of irreducibility, and the weak irreducibility given in [13] for finite-state Markov chains. Compared with the usual definition of irreducibility, it deals with time-varying generators and allows some components of the quasi-stationary distribution to be 0. Suppose that the Markovian queueing network has n0 (n0 < ∞) interconnected stations. Consider a vector-valued queue-length process taking values in Nn0 = N × N × · · · × Nn0 (an n0 -fold product). We use α(·), a continuous-time Markov chain, to model the queue length of the queueing network. Suppose that Nn0 can be divided into l subsets. Within each subset, the transitions (such as arrivals and departures of customers, etc.) take place an order of magnitude more frequent than that of among different subsets. To highlight this contrast, we introduce a small parameter ε > 0. Thus the continuous-time Markov chain is ε-dependent, i.e., α(·) = αε (·). For t ∈ [0, T ], assume that the generator of the Markov chain αε (t) is given by (2.2)
Qε (t) =
A(t) + B(t), with A(t) = diag(A1 (t), . . . , Al (t)), ε
where A(t), B(t), and Ai (t) (for i = 1, . . . , l) are all generators of certain countablestate space Markov chains, and diag(A1 (t), . . . , Al (t)) denotes a block diagonal matrix, each of the entries of which has appropriate dimension. Note that B(t) is a generator and there is no need to assume it has the same diagonal matrix form as that of A(t). There is a certain hierarchy in the underlying network. Within each subset, one observes detailed variations of the networks such as arrivals and services of customers, etc., whereas at an upper system management level, instead of these variations, one observes the transitions among different subsets. Thus from an upper management point of view, by lumping all the states in each subspace into a single one, the system may be regarded as if it were a queue with finite waiting rooms (l rooms). Nevertheless, the aggregated process is generally non-Markovian. Fortunately, as will be shown, the aggregated process converges weakly to a limit process that is a finite-state Markovian queue. The significance of such a result is that in lieu of examining the detailed variations, one can study the aggregated process. To reflect the decomposition in the network, write the state space of the queueing network as M = M1 ∪ M2 · · · ∪ Ml , where Mi = {si1 , si2 , . . . , } for i = 1, . . . , l. Following our approach for singularly perturbed Markov chains with finite-state spaces, we construct matched asymptotic expansions of the probability vector (2.3)
21 l1 ij pε (t) = (p11 ε (t), . . . , pε (t), . . . , . . . , pε (t), . . .), pε (t) = P (αε (t) = sij ).
For future use, partition an infinite dimensional vector v as v = (v 1 , . . . , v l ) with v i = (v i1 , v i2 , . . .). That is, v i is an infinite dimensional vector corresponding to the subspace Mi . Since we are dealing with countable state space Markov chains, we work with an infinite dimensional vector space. It is naturalto consider the spaces ∞ l 1 = {(v 1 , . . . , v l ) : 1 ≤ i ≤ l, v ik ∈ R for each k ∈ N, and i=1 k=1 |v ik | < ∞}, ∞ = {(v 1 , . . . , v l ): 1 ≤ i ≤ l, v ik ∈ R for each k ∈ N, and sup1≤i≤l sup1≤k 0, exp(A(0)τ ) − diag(11ν 1 (0), . . . , 11ν l (0)) ≤ K exp(−κτ ), (3.1) ∞ where ν i (t) = (ν i1 (t), ν i2 (t), . . .) is the quasi-stationary distribution corresponding to the generator Ai (t). Remark 3.1. From (A1) and the definition of weak irreducibility, ν i (t)Aia (t) = · · (0 : 1) has a unique solution. Furthermore, ν i (t) = (0 : 1)(Aia (t)) (Aia (t)(Aia (t)) )−1 . (A2) is a Doeblin-type condition. A condition in a slightly different form is given in [7, p. 192]. We will obtain the asymptotic expansions in two steps. The first step is a formal construction in which we (see [4, 10, 25] and [11, 13] among others) find the outer expansions and the initial layer corrections. The second step involves validating the formal expansions and deriving the desired error estimates. 3.1. Formal expansions. Weseek matched asymptotic expansions of the form n n Oε,n (t) + Iε,n (t) = i=0 εi φi (t) + i=0 εi ψi (t/ε). For technical reasons, which will become clear in what follows, to justify the validity of the expansions, we also need to compute φn+1 (t) and ψn+1 (t/ε). To make sure that the matching condition is satisfied, we use Oε,n+1 (0) + Iε,n+1 (0) = p(0) or, more specifically, φ0 (0) + ψ0 (0) = p(0), φk (0) + ψk (0) = 0 for 1 ≤ k ≤ n + 1. Our approach is based on constructions of the sequences {φi (t)} and {ψi (t/ε)}. Substituting Oε,n+1 (t) into (2.4) and comparing coefficients of powers of εk , we obtain (3.2)
φ0 (t)A(t) = 0,
φk (t)A(t) = φ˙ k−1 (t) − φk−1 (t)B(t), 1 ≤ k ≤ n + 1.
The outer expansions give us satisfactory approximation for t > 0 away from an initial layer of the order O(ε), but it does not satisfy the initial condition and breaks down for sufficiently small t. To compensate, introduce a fast time variable τ = t/ε. By the Lipschitz continuity given in (A1), taking Taylor expansions of A(ετ ) and B(ετ ) i i i i n+1 n A(0) B(0) about 0 yields A(ετ ) = i=0 (ετi!) d dt + O((ετ )n+2 ), εB(ετ ) = i=0 ε (ετi!) d dt + i i n+2 O((ετ ) ). Substituting Iε,n+1 (t) into (2.4), using the above Taylor expansions of A(t) and B(t), and comparing powers of εi , we obtain the equations satisfied by the initial layer terms. We have ψ0 (0) = p(0) − φ0 (0), ψk (0) = −φk (0), and
(3.3)
dψ0 (τ ) dψk (τ ) = ψ0 (τ )A(0), = ψk (τ )A(0) + rk (τ ), 1 ≤ k ≤ n + 1, dτ dτ k−1 τ i+1 di+1 A(0) τ i di B(0) , 1 ≤ k ≤ n + 1. ψk−i−1 (τ ) + rk (τ ) = (i + 1)! dti+1 i! dti i=0
In [28], in which A(t) consists of only one block, the outer expansions and initial layers can be obtained separately. Here, the constructions of {φk (t)} have to be done
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
571
in conjunction with the initial layer corrections. The equations in (3.3) together with the initial data are known as abstract Cauchy problems; see [17, p. 21]. It is well known that {exp(A(0)τ )} is a uniformly continuous semigroup (see [17, 23]). With ψk (0) to be determined, the representation of the solutions of (3.3) is given by (3.4)
ψ0 (τ ) = (p(0) − φ0 (0)) exp(A(0)τ ), τ
rk (s) exp(A(0)(τ − s))ds, 1 ≤ k ≤ n + 1.
ψk (τ ) = ψk (0) exp(A(0)τ ) + 0
Step 1. Determine φ0 (t) and ψ0 (τ ). We begin with the first equation in (3.2). Using the partitioned vector form introduced right after (2.3), we obtain are not uniquely solvable φi0 (t)Ai (t) = 0 for each i = 1, . . . , l. These equations ∞ i since Ai (t) have a 0 eigenvalue. However, by attaching j=1 φij 0 (t) = θ0 (t) to the equations, the resulting system of equations has a unique solution thanks to the weak irreducibility of Ai (t). Thus, φi0 (t) must be proportional to ν i (t), the quasi-stationary distribution corresponding to Ai (t). That is, φi0 (t) = θ0i (t)ν i (t), where θ0i (t) ∈ R is to be determined. Define θ0 (t) = (θ01 (t), . . . , θ0l (t)) ∈ R1×l , 11 = diag(11, 11, . . . , 11), ν(t) = diag(ν 1 (t), . . . , ν l (t)). (3.5) It is immediate that A(t) is orthogonal to 11 (i.e., A(t)11 = 0) and φ0 (t) = θ0 (t)ν(t). For the equation with k = 1 in (3.2), multiplying from the right by 11 leads to (3.6)
θ˙0 (t) = θ0 (t)B(t),
B(t) = ν(t)B(t)11,
so B(t) is an average of B(t) with respect to ν 1 (t), . . . , ν l (t). Note that B(t) is an l × l matrix-valued function. Note also that (3.6) is a linear system of differential equations and that it has a unique solution for each initial condition. Choose θ0 (0) = p(0)11. Then the solution of (3.6) is uniquely determined. The φi0 (t) has the interpretation of total probability. Concerning the initial layer correction ψ0 (τ ), the solution is given by the first equation in (3.4) and is uniquely solved. We claim that ||ψ0 (τ )||∞ ≤ K exp(−κτ ) for some κ > 0 and K > 0. To prove this, note that 11ν(0) = diag(11ν 1 (0), . . . , 11ν l (0)) and (3.7)
ψ0 (0)11ν(0) = [p(0) − φ0 (0)]11ν(0) = (0, 0, . . .),
where φ0 (0) = θ0 (0)ν(0), φ0 (0)11 = θ0 (0), and θ0 (0) = p(0)11 are used. By virtue of the orthogonality and (A2), ||ψ0 (τ )||∞ = ||ψ0 (0)(exp(A(0)τ ) − 11ν(0))||∞ ≤ K exp(−κτ ). To determine φi (t) and ψi (t) for i ≥ 1, we need the following lemma. Lemma 3.2. Suppose that Q(t) is a generator of a countable-state-space Markov chain such that Q(t) is weakly irreducible for each t ∈ [0, T ] and (dn+1 /dtn+1 )Q(·) ∈ −1 (t) exists = Qa (t)Q (t). Then for k = 1, . . . , n+1, (dk /dtk )Q Lip[0, T ]. Denote Q(t) a n+1−k and belongs to C (the class of functions that are (n + 1 − k)-times continuously differentiable). Q −1 (t) = I, differentiating both sides of the above equation Proof. Since Q(t) leads to (3.8)
−1 (t) dQ −1 (t) dQ(t) Q −1 (t), = −Q dt dt
−1 (·) is differentiable. Repeatedly differentiating equation (3.8) yields the desired so Q result.
572
G. YIN AND HANQIN ZHANG
Step 2. Determine φ1 (t) and ψ1 (τ ). Lemma 3.2 and Remark 3.1 imply that φ0 (·) is differentiable and is a function of class C n+1 . We proceed to determine φ1 (t) and ψ1 (τ ). Note that the equation with k = 1 in (3.2) is a nonhomogeneous equation def whose right-hand side φ0 (t) = φ˙ 0 (t) − φ0 (t)B(t) is a known function since φ0 (t) has been found. Using (3.6), [φ˙ 0 (t)−φ0 (t)B(t)]11 = θ˙0 (t)−θ0 (t)ν(t)B(t)11 = (0, 0, . . .), and the Fredholm alternative yields that the equation with k = 1 in (3.2) has a particular solution φ1,p (t) being orthogonal to 11. Assume the solution of φ1 (t)A(t) = φ0 (t) to be of the form φ1 (t) = θ1 (t)ν(t)+φ1,p (t). Since φ1,p (t) is orthogonal to 11, postmultiplying the equation with k = 2 in (3.2) by 11 leads to θ˙1 (t) = θ1 (t)B(t) + φ1,p (t)B(t)11.
(3.9)
Once the initial condition is specified, (3.9) is uniquely solved. The initial condition θ1 (0) has to come from the initial layer correction term. With the selection of ψ1 (0) = −φ1 (0), by (3.4) with k = 1, the unique solution is given by τ ψ0 (s) exp(A(0)s)B(0) exp(A(0)(τ − s))ds ψ1 (τ ) = ψ1 (0) exp(A(0)τ ) + 0 τ (3.10) dA(0) exp(A(0)(τ − s))ds. + sψ0 (s) exp(A(0)s) dt 0 τ ∞ By the exponential decay of ψ0 (τ ), 0 ψ0 (s) exp(A(0)s)ds∞ = 0 ||ψ0 (s)||∞ ds < def ∞ ∞. Denote ψ 0 = 0 ψ0 (s) exp(A(0)s)dsB(0). Then from (A2),
lim ψ1 (0) exp(A(0)τ ) = lim ψ1 (0)[exp(A(0)τ ) − 11ν(0)] + ψ1 (0)11ν(0) τ →∞ τ →∞ = ψ1 (0)11ν(0), τ ψ0 (s) exp(A(0)s)B(0) exp(A(0)(τ − s))ds = ψ 0 11ν(0). lim τ →∞
0
Note that (3.11)
A(t)11 = 0,
dk A(t) dk A(t)11 = 0, 1 1 = dtk dtk
k = 1, . . . , n + 1.
Using the orthogonality (see (3.7) and (3.11)), we obtain τ dA(0) exp(A(0)(τ − s))ds ≤ Kτ 2 exp(−κτ ). sψ0 (s) exp(A(0)s) (3.12) dt 0 ∞ By demanding limτ →∞ ψ1 (τ ) = 0, taking the limit as τ → ∞ in (3.10), and using the estimates above, we arrive at (3.13)
ψ1 (0)11ν(0) + ψ 0 11ν(0) = 0.
An important observation indicates that there are only l unknowns in (3.13). Using the notation of partitioned vector given after (2.3), ψ1 (0) = (ψ11 (0), ψ12 (0), . . . , ψ1l (0)) 1 2 l i and ψ 0 = (ψ 0 , ψ 0 , . . . , ψ 0 ), the solution is given by ψ1i (0)11 = −ψ 0 11, i = 1, . . . , l. i Since ψ1 (0) must be chosen so that φ1 (0) + ψ1 (0) = 0, we obtain θ1i (0) = ψ 0 11, i = 1, . . . , l. Thus both φ1 (t) and ψ1 (τ ) have been determined. We claim that ψ1 (τ )
573
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
decays exponentially fast. By adding and subtracting appropriate terms, ∞ ψ0 (s) exp(A(0)s)dsB(0)11ν(0) ψ1 (τ ) = ψ1 (0)[exp(A(0)τ ) − 11ν(0)] − τ ∞ + ψ1 (0)11ν(0) + ψ0 (s) exp(A(0)s)dsB(0)11ν(0) 0 (3.14) τ + ψ0 (s) exp(A(0)s)B(0)[exp(A(0)(τ − s)) − 11ν(0)]ds 0 τ dA(0) exp(A(0)(τ − s))ds. + sψ0 (s) exp(A(0)s) dt 0 By (A2),
(3.15)
ψ1 (0)[exp(A(0)τ ) − 11ν(0)] ≤ K exp(−κτ ), and ∞ τ ψ0 (s) exp(A(0)s)B(0)[exp(A(0)(τ − s)) − 11ν(0)]ds 0
≤ K(1 + τ ) exp(−κτ ).
∞
It then follows from the above estimates, ||ψ1 (τ )||∞ ≤ K(1 + τ + τ 2 ) exp(−κτ ) ≤ K exp(−κ0 τ ), for some 0 < κ0 < κ. Step 3. Determine φk (t) and ψk (τ ), for 1 < k ≤ n + 1. We apply the same method for finding φ1 (t) and ψ1 (τ ) to determine φk (t) and ψk (τ ). For each i = 1, . . . , l, assume φik (t) is of the form φik (t) = θki (t)ν i (t) + φik,p (t), where φk,p (t) = (φ1k,p (t), . . . , φlk,p (t)) is a particular solution of the equation φk (t)A(t) = φ˙ k−1 (t) − φk−1 (t)B(t), and φk,p (t) is orthogonal to 11. We proceed to determine θk (t) = (θki (t)) ∈ R1×l . By substitution, similar to (3.9), it is easily seen that θk (t) satisfies the differential equation θ˙k (t) = θk (t)B(t) + φk,p (t)B(t)11.
(3.16)
To determine the initial condition θk (0), the definition of rk (τ ) in (3.4) gives us ψk (τ ) = ψk (0) exp(A(0)τ ) + 0
(3.17) ×
τ k−1
ψk−j−1 (s)
j=0 j j
τ j+1 dj+1 A(0) τ d B(0) exp(A(0)(τ − s))ds. + (j + 1)! dtj+1 j! dtj
Using the same techniques as that for θ1 (0), we obtain θk (0) and uniquely determine both φk (t) and ψk (τ ). We then verify the exponential decay property of ψk (τ ). The procedure is the same as for that of ψ1 (τ ); we record the result as follows. Theorem 3.3. Assume (A1) and (A2). Then the following assertions hold: (a) The sequence {φk (t)} can be constructed by solving the system of equations (3.18)
def φk (t)A(t) = φ˙ k−1 (t) − φk−1 (t)B(t) = φk−1 (t)
for k = 1, 2, . . . , n. The solution is φk (t) = θk (t)ν(t) + φk,p (t), where φk,p (t) is a particular solution of (3.18), which is orthogonal to 11, and θk (t) satisfies θ˙k (t) = θk (t)B(t) + φk,p (t)B(t)11. j k−1 ∞ j def A(0) )11ν(0) = (b) Find ψki (0)11 from ψk (0)11ν(0) = −( j=0 0 sj! ψk−j−1 (s)ds d dt j i −ψ k−1 11ν(0). Choose θki (0) = −ψki (0)11 = ψ k−1 11 for i = 1, . . . , l. Choose ψk (0) = −φk (0).
574
G. YIN AND HANQIN ZHANG
(c) For k = 0, 1, . . . , n+1, θk (·) and φk (·) are (n+1−k)-times continuously differentiable on [0, T ]; there exists a κ0 > 0 such that ||ψk (τ )||∞ ≤ K exp(−κ0 τ ). 3.2. Asymptotic justification. In this section, we obtain the desired error bounds and show that the asymptotic expansions hold uniformly in t ∈ [0, T ]. Define an operator Lε as (3.19)
Lε f (t) = ε
df (t) − f (t)A(t) − εf (t)B(t) dt
for a suitable smooth function f (·). First let us establish a lemma. Lemma 3.4. Consider Lε vε (t) = ∆ε,k (t), vε (0) = 0, with supt∈[0,T ] ||∆ε,k (t)||∞ = O(εk+1 ) for some k with 0 ≤ k ≤ n + 1. Then supt∈[0,T ] ||vε (t)||∞ = O(εk ). Proof. Note that the initial value problem given above is a time-dependent abstract Cauchy problem or an evolution equation. The solution is given by vε (t) = 1 t ε 0 ∆ε,k (s)Xε (t, s)ds, where Xε (t, s) is a fundamental solution or an evolution operator. In fact (see Ladas and Lakshmikantham [17, p. 56]), Xε (t, s) is an operator-valued function with values in L(1 ); the space of bounded linear operators, defined on 1 and strongly continuous in t, s for 0 ≤ s ≤ t ≤ T such that (∂/∂t)Xε (t, s) exists in strong topology of 1 ; (∂/∂t)Xε (t, s) ∈ L(1 ) for 0 ≤ s ≤ t ≤ T and (∂/∂t)Xε (t, s) is strongly continuous in t; the range of Xε (t, s) is in the domain of Qε (t); and it satisfies the homogeneous problem ε ∂Xε∂t(t,s) − Xε (t, s)A(t) − εXε (t, s)B(t) = 0, Xε (s, s) = I, where the initial value I is the infinite dimensional identity matrix. Since it represents transition probabilities, ||Xε (t, s)||∞ is bounded uniformly in ε for all t, s ∈ [0, T ]. t k Therefore, we have supt∈[0,T ] ||vε (t)||∞ ≤ K ε supt∈[0,T ] 0 ||∆ε,k (s)||∞ ds ≤ Kε . The lemma is obtained. For each k = 1, . . . , n + 1, define a vector-valued error eε,k (t) (3.20)
eε,k (t) = pε (t) −
k i=0
εi φi (t) −
k
εi ψi (t/ε),
i=0
where pε (·) is the solution of (2.4), φi (·) and ψi (·) are the outer expansions and initial layer corrections, respectively. We must show that eε,n (t) = O(εn+1 ). For this purpose, we first derive a lemma whose proof is similar in spirit to the corresponding results for weakly irreducible generators [28] and is thus omitted. Lemma 3.5. Under (A1) and (A2), for k = 1, . . . , n+1, supt∈[0,T ] ||Lε eε,k (t)||∞ = k+1 O(ε ), and hence supt∈[0,T ] ||eε,k (t)||∞ = O(εk ). Remark 3.6. Using Lemma 3.5, with k = 1, we obtain supt∈[0,T ] ||eε,1 (t)||∞ = O(ε). However, eε,1 (t) = eε,0 (t) − εφ1 (t) − εψ1 (t/ε). In view of the boundedness of φ1 (t) and ψ1 (t/ε), εφ1 (t) + εψ1 (t/ε) = O(ε). Thus, eε,0 (t) = O(ε) uniformly in t ∈ [0, T ]. We can proceed inductively. By virtue of Lemma 3.5, with k = n + 1, supt∈[0,T ] ||εε,n+1 ||∞ = O(εn+1 ). Going back one step and using eε,n+1 (t) = eε,n (t) + O(εn+1 ), similar to the case of k = 1, we obtain supt∈[0,T ] ||eε,n (t)||∞ = O(εn+1 ). We summarize the discussion thus far into the following theorem. This gives us the desired approximation results. Not only is the convergence of P (αε (t) = i) obtained, but also the rate of convergence is derived. Furthermore, a full asymptotic series is obtained. Theorem 3.7. Under (A1) and (A2), we can construct two sequences {φk (t)}nk=0 and {ψk (t/ε)}nk=0 by Theorem 3.3 such that φk (t) ∈ C n+1−k and ψk (t/ε) decay exponentially fast. Moreover with eε,n (t) defined by (3.20), supt∈[0,T ] ||eε,n (t)||∞ = O(εn+1 ).
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
575
Using the same techniques, we can also obtain asymptotic expansions of transition probability matrices. Since the proofs are essentially the same, we will state only the results and omit the detailed argument. Let Pε (t0 , t) be the transition matrix (pin,jk (t0 , t)) with pin,jk (t0 , t) = P (αε (t) = sjk |αε (t0 ) = sin ) for all sin , sjk ∈ M. ε ε Theorem 3.8. Assume (A1) and (A2) with n = 1. Then (3.21)
Pε (t0 , t) = Φ0 (t0 , t) + Ψ0 (t0 , (t − t0 )/ε) + O(ε + ε exp((t − t0 )/ε))
uniformly in (t0 , t), where 0 ≤ t0 ≤ t ≤ T , Φ0 (t0 , t) = 11Θ(t0 , t)ν(t),
dΨ0 (t0 , τ ) = Ψ0 (t0 , τ )A(t0 ), Ψ0 (t0 , t0 ) = I − Φ0 (t0 , t0 ), dτ
(3.22) where Θ(t0 , t) = (θij (t0 , t)) ∈ Rl×l is the solution of
d Θ(t0 , t) = Θ(t0 , t)B(t), Θ(t0 , t0 ) = I. dt
(3.23)
Remark 3.9. Owing to Theorems 3.7 and 3.8, for some κ0 > 0, pε (t) = ν(t) + O(ε + exp(−κ0 t/ε)), Pε (t0 , t) = Φ0 (t0 , t) + O(ε + exp(−κ0 (t − t0 )/ε)). 4. Occupation measures, aggregation, and switching diffusion. This section presents further asymptotic results that are of probabilistic feature. The statements of the results are given, and the proofs are relegated to the appendix. 4.1. Occupations measures. For any i = 1, . . . , l and j = 1, 2, . . ., define sequences of occupation measures as t ij (4.1) µε (t) = z ij (s, αε (s))ds, z ij (t, αε (t)) = I{αε (t)=sij } − ν ij (t)I{αε (t)∈Mi } , 0
where ν ij (t) is the jth component of the quasi-stationary distribution ν i (t) as defined in (A2). As a preparation, we first derive an order of magnitude estimate for µij ε (·) defined in (4.1). 2 Theorem 4.1. Under the conditions of Theorem 3.8, supt∈[0,T ] E(µij ε (t)) = O(ε). As alluded to in the introduction, we wish to reduce the complexity of the queueing network by aggregating states in each subspace as a single state. This leads to the definition of αε (t) = i if αε (t) ∈ Mi . Theorem 4.2. Under the conditions of Theorem 3.8, αε (·) converges weakly to α(·), a Markov chain generated by B(t) given in (3.6). 4.2. Switching diffusion limit. Let f (·) be a real-valued function defined on M satisfying that {f (sij ) : 1 ≤ i ≤ l, 1 ≤ j < ∞} ∈ 1 . For t ∈ [0, T ], define a sequence of real-valued functions by (4.2)
xε (t) =
t l ∞ 1 √ f (sij )[I{αε (u)=sij } − ν ij (u)I{αε (u)∈Mi } ]du. ε 0 i=1 j=1
Our interest lies in the asymptotic properties of xε (·). To obtain the desired limit property, we consider a pair of processes {Yε (·)} = {xε (·), αε (·)} and aim to show that Yε (·) converges weakly to Y (·) with a suitable generator.
576
G. YIN AND HANQIN ZHANG
Lemma 4.3. Under the conditions of Theorem 3.8, {Yε (·)} is tight in D2 [0, T ] (the space of R2 -valued functions that are right continuous and have left limits endowed with the Skorohod topology). Since {Yε (·)} is tight, by Prohorov’s theorem, we can extract a weakly convergent subsequence. Select such a subsequence and still use ε as its index for notational simplicity. Denote the limit by Y (·) = (x(·), α(·)). We proceed to characterize this limit process. Let CL2 be the collection of functions having bounded derivatives up to the second order with the second derivative being Lipschitz continuous. For each i = 1, . . . , l, g(·, i) ∈ CL2 , define an operator D(t) by (4.3)
D(t)g(x, i) =
∂2 1 2 σ (t, i) 2 g(x, i) + B(t)g(x, ·)(i), 2 ∂x
where σ 2 (t, i) > 0 is a smooth function. Theorem 4.4. Under the conditions of Theorem 3.8, Yε (·) converges weakly to Y (·), which is a solution of the martingale problem with operator D(·) defined by (4.3). Using the above theorem, we can also get a weak convergence of a reflected process. Define rε (t) = xε (t)−inf 0≤u≤t {xε (u)}. The weak convergence of xε (t) and the continuous mapping theorem (see [2, p. 30, Theorem 5.1]) yield that rε (·) converges weakly to r(·), a reflected switching diffusion process. 5. Extensions and examples. In this section we give the extensions of the results and present examples of queueing problems. 5.1. Infinitely many blocks. We consider the case (5.1)
Qε (t) =
A(t) + B(t) ε
with A(t) = diag(A1 (t), A2 (t), . . .).
That is, the matrix A(t) given by (2.2) is a diagonal block one with infinitely many blocks. Instead of condition (A2), we assume (A2 ). l (A2 ) There is a κl > 0 such that for exp(A (0)τ ) − 11ν l (0)∞ ≤ K exp(−κl τ ), any real number τ > 0, where ν l (t) = (ν l1 (t), ν l2 (t), . . .) is the quasi-stationary distribution corresponding to the generator Al (t), and inf l≥1 κl > 0. Theorem 5.1. Under (A1) and (A2 ), Theorem 3.7, Theorem 4.2, and Theorem 4.4 hold. The proof of this theorem follows the same line of argument as that of Theorems 3.7, 4.2, and 4.4 and is thus omitted. A special case of the theorem is that Qε (t) has the form (5.1) with A(t) having infinitely many blocks such that each Ai (t) is a generator of a finite-state Markov chain. Example 5.2. Consider a two-station queueing system, where each station has a single server and an exogenous Poisson arrival process. For the first station, having completed service there, customers either proceed to a queue in front of the second station with probability p1 or depart the system with probability (1 − p1 ). For the second station, after completing service there, they depart the system with probability (1 − p2 ) and go to a queue in front of the first station with probability p2 . Service is rendered in the order of the arrivals at each station. We assume that the first station can hold at most a total of m0 customers (including the customer in service) and any further arriving customers from outside will in fact be refused entry and hence will depart immediately without service at the first station, while the second station
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
577
has an unlimited waiting room. Furthermore, we assume that the system is a timedependent Markovian queueing network. Specifically, the rate of station i’s customer arriving from outside is λi (t), and the service rate is µi (t) (i = 1, 2). Assume that λi (t) and µi (t) are smooth, positive, real analytic functions of time t. We let Q(t) = (Q2 (t), Q1 (t)) be the Markovian queue length process that Qi (t) equals the number of customers at station i at time t. Then the state space of the queueing system can be represented by {(k, i) : 0 ≤ k < ∞ and 0 ≤ i ≤ m0 }. To determine the generator Q(t) = (q (j,n),(k,i) (t)), let π1 (t) = (1 − p1 )µ1 (t), π2 (t) = (1 − p2 )µ2 (t), π(t) = λ1 (t) + π1 (t), π ˜ (t) = λ2 (t) + p1 µ1 (t) + µ2 (t), ⎞ ⎛ −λ (t) λ (t) 1 1 ⎟ ⎜ π1 (t) −π(t) λ1 (t) ⎟ ⎜ 1 .. .. .. ⎟, A (t) = ⎜ . . . ⎟ ⎜ ⎠ ⎝ π1 (t) −π(t) λ1 (t) π1 (t) −π1 (t) B 1 (t) = diag(−λ2 (t), −[λ2 (t) + p1 µ1 (t)], . . . , −[λ2 (t) + p1 µ1 (t)]), B 2 (t) = diag(−(λ2 (t) + µ2 (t)), −˜ π (t), . . . , −˜ π (t)), ⎞ ⎛ λ (t) 2 ⎟ ⎜ p1 µ1 (t) λ2 (t) ⎟, C(t) = ⎜ .. .. ⎠ ⎝ . . ⎛ π (t) p µ (t) 2 2 2 .. ⎜ . D(t) = ⎜ ⎝
p1 µ1 (t) λ2 (t) ..
. π2 (t) p2 µ2 (t) µ2 (t)
⎞ ⎟ ⎟, ⎠
where each matrix is a (m0 + 1) × (m0 + 1) matrix. Then some tedious calculations yield ⎛ 1 ⎞ B (t) C(t) 2 ⎜ ⎟ (5.2) Q(t) = diag(A1 (t), A1 (t), . . .) + ⎝ D(t) B (t) C(t) ⎠. .. .. .. . . . For an initial time point t0 with t0 ∈ [0, T ], let P (t0 , t) with t > t0 be the transition matrix P (t0 , t) = (p(j,n),(k,i) (t0 , t)) with p(j,n),(k,i) (t0 , t) = P (Q(t) = (k, i)|Q(t0 ) = (j, n)) for all 0 ≤ j, k < ∞, 0 ≤ n, i ≤ m0 . Then we have the following system of equations for a quasi-birth-death process: d P (t0 , t) = P (t0 , t)Q(t). dt
(5.3) Assume that (5.4)
(1 − p1 ) p1 , λ1 (t) p1 µ1 (t), λ1 (t) λ2 (t), λ1 (t) µ2 (t). 1 inf t0 ≤t≤T λ1 (t) .
Then ν i (t) = (ν i0 (t), ν i1 (t), . . . , ν im0 (t)) with ⎤ ⎡ j j m0 λ1 (t) ⎦ λ1 (t) . ν 1j (t) = p0 (t), p0 (t) = 1 ⎣1 + π1 (t) π1 (t) j=1
Introduce ε =
578
G. YIN AND HANQIN ZHANG
Let π (t) = [λ2 (t) + p1 µ1 (t)] − p1 µ1 (t)p0 (t). It follows from some tedious calculations that ⎞ ⎛ − π (t) π (t) ⎟ ⎜ π (t) − µ2 (t) π (t) (5.5) B(t) = ⎝ µ2 (t) − ⎠. .. .. .. . . . By Theorem 5.1,
(5.6)
t − t0 Pε (t0 , t) = Φ0 (t0 , t) + [I − Φ(t0 , t0 )] exp A(t0 ) ε t − t0 + O ε + exp ε
uniformly in (t0 , t), where for 0 ≤ t0 ≤ t ≤ T , Φ0 (t0 , t) is given by (3.22) with Θ(t0 , t) satisfying (3.23). Using [23, Theorem 5.2, p. 128], the solution of (3.23) can be explicitly written as Θ(t0 , t) = U (t0 , t), where U (t0 , t) is a solution operator. Consequently Φ(t0 , t) given in (5.6) takes the form Φ(t0 , t) = 11U (t0 , t)ν. In particular, when the generators are time independent, U (t0 , t) becomes a semigroup [23, Chapter 1] and may be expressed as U (t0 , t) = exp((t − t0 )B). Next we consider the queue length process at the first station, Q1 (t). In general, Q1 (t) is not a Markov process. But from Theorem 5.1, Q1 (t) can be approximated very well by a Markov process α(t) with generator B(t) defined by (5.5). Furthermore, we consider a refined approximation of Q(t) = (Q2 (t), Q1 (t)). Define ∞ m0 t I{Q(u)=(i,j)} − ν (i,j) (u)I{Q1 (u)=j} du. It follows from Theoε (t) = i=0 j=0 0 √ rem 5.1 that (1/ ε)ε (t) can be well approximated by the switching diffusion prot cess 0 σ(u, α(u))dw(u), where w(·) is a standard one-dimensional Brownian motion, ∞ lk m0 m0 0k ∞ kl 0l σ 2 (t, i) = k=0 l=0 ν (t) 0 ψ0 (i, t, u)du + ν (t) 0 ψ0 (i, t, u)du for 0 ≤ i < ∞, and ψ0kl (i, s, t) is the (k, l)th entry of Ψ0 (i, t, u) given by ⎛ ⎛ 1 ⎞⎞ ν (t) · 1 ⎝ ⎝ Ψ0 (i, t, u) = I − : ⎠⎠ exp(A (t)u). 1 ν (t) Remark 5.3. References [12] and [18] establish asymptotic expansions for the queue length distribution of the time inhomogeneous single serve queue, a timedependent pure birth-death process. Here we consider a quasi-birth-death process. Furthermore, the queue length process at the first station (generally non-Markovian) can be approximated well by a Markov process with generator B(t). In [20], uniform acceleration expansions for time-varying generators were treated, and in [28] diffusion approximation was also considered. In these references, the Markov chains have one ergodic class, whereas in the current paper, multiergodic classes are considered. Compared with the diffusion approximations in [19] and [27], the switching diffusion approximation given here is related to the queue length process on the interval [0, t] and provides the evolution of the scaled sequence. The usual diffusion approximation leads to asymptotic normality, whereas switching diffusion limit yields Gaussian mixture distribution. The [λ1 (t) + p2 λ2 (t)]/[(1 − p1 p2 )µ1 (t)] and [λ2 (t) + p1 µ1 (t)]/[(1 − p1 p2 )µ2 (t)] are called the traffic intensities at station 1 and station 2, respectively. If [λ1 (t) + p2 λ2 (t)]/[(1 − p1 p2 )µ1 (t)] and [λ2 (t) + p1 µ1 (t)]/[(1 − p1 p2 )µ2 (t)] are less than but close
579
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
to one, the system is regarded as in heavy traffic, whereas if [λ1 (t) + p2 λ2 (t)]/[(1 − p1 p2 )µ1 (t)] and [λ2 (t) + p1 µ1 (t)]/[(1 − p1 p2 )µ2 (t)] are much less than one, the system is considered to be in light traffic. Our approximations are valid for either the heavy traffic case or the light traffic case as long as (5.4) holds. Using the Laplace transform technique, from the balance equation of the system, one may proceed as in [15] to carry out heavy traffic analysis for the sojourn time of time-homogeneous Markovian tandem queues with two servers. Comparing with (5.2), we can consider a general quasi-birth-death process with state space {(i, j), i ≥ 0, 1 ≤ j ≤ m0 } with a generator G(t) given by ⎛
⎞
A1 (t) − diag(C 1 (t)11) C 1 (t) 1 2 ⎜ B (t) A (t) − diag(B 1 (t)11 + C 2 (t)11) C 2 (t) (5.7) ⎝ .. .. . .
⎟ ⎠,
where i and j are called level and phase, respectively; see [22]. If the transitions between levels are much less frequent than the transitions between the phases inside the same level, then Q(t) given by (5.7) can be written as diag(A1⎛ (t), A2 (t), . . .) −diag(C 1 (t)11) C 1 (t) 1 1 ⎜ B (t) −diag(B (t)11 + C 2 (t)11) C 2 (t) +ε⎝ .. .. . .
⎞ ..
⎟ ⎠. .
Hence following Example 5.2, we can also get the asymptotic probability distribution of the process. 5.2. Asymptotic properties under γ-norm. For v = (v 1 , . . . , v l ) with v i = (v , v i2 , . . .), the γ-norm (which was named ν-norm in [21]; see also [1]) is defined as ||v||γ = max1≤i≤l sup1≤k 0, where ν l (t) = (ν l1 (t), ν l2 (t), . . .) is the quasi-stationary distribution corresponding to the generator Al (t), and inf l≥1 κl > 0. i1
Similar to Theorem 3.7, using essentially the same techniques, we obtain the following result. Theorem 5.4. Under (A1) and (A2 ), we construct two sequences {φk (t)}nk=0 and {ψk (t/ε)}nk=0 by Theorem 3.3 such that φk (t) ∈ C n+1−k and ψk (t/ε) decay exponentially fast. Moreover, with eε,n (t) defined by (3.20), supt∈[0,T ] ||eε,n (t)||γ = O(εn+1 ). Example 5.5. Consider the queueing system given by Example 5.2, but allow the first station to have unlimited waiting rooms. Then the state space of (Q2 (t), Q1 (t))
580
G. YIN AND HANQIN ZHANG
can be represented by {(k, i) : 0 ≤ k < ∞ and 0 ≤ i < ∞}. Let ⎞ ⎛ −λ1 (t) λ1 (t) ⎟ 1 (t) = ⎜ A ⎠, ⎝ (1 − p1 )µ1 (t) −π(t) λ1 (t) .. .. .. . . . 1 (t) = diag(−λ2 (t), −[λ2 (t) + p1 µ1 (t)], . . .), B 2 (t) = diag(−(λ2 (t) + µ2 (t)), −˜ π (t), B ⎞ . . .), ⎛ ⎛ π2 (t) p2 µ2 (t) λ2 (t) ⎟ ⎜ (t) λ (t) p π2 (t) p2 µ2 (t) µ =⎜ 1 2 1 , D(t) = C(t) ⎠ ⎝ ⎝ .. .. .. . . .
⎞ ..
⎟ ⎠. .
Then, some tedious calculations yield ⎛
(5.8)
1 (t) C(t) B ⎜ 1 1 (t), . . .) + ⎝ D(t) B 2 (t) C(t) = diag(A (t), A Q(t) .. .. . .
⎞ ..
⎟ ⎠. .
Letting γ 1i (0) = ((1−p1 )µ1 (0)/λ1 (0))i , then we have (A2 ) holds; see [21]. Therefore, by Theorem 5.4, we obtain the asymptotic probability distribution of (Q2 (t), Q1 (t)) under γ-norm. Appendix A. Proof of Theorem 4.1. Suppose that 0 ≤ s ≤ t. For each i = 1, . . . , l and j ∈ N, define ζε1,ij (s, t) = P (αε (s) = sij , αε (t) = sij ) − ν ij (t)P (αε (s) = sij , αε (t) ∈ Mi ), (A.1) ζε2,ij (s, t) = ν ij (s)ν ij (t)P (αε (s) ∈ Mi , αε (t) ∈ Mi ) −ν ij (s)P (αε (s) ∈ Mi , αε (t) = sij ). Theorem 3.8 and Remark 3.9 yield ζε1,ij (s, t) = O(ε + exp(−κ0 (t − s)/ε)). Similarly, ij 2 ζε2,ij (s, t) = O(ε + exp(−κ0 (t − s)/ε)). It follows from µij ε (0) = 0 that E(µε (t)) = t s 1,ij 2,ij 2 0 0 (ζε (r, s) + ζε (r, s))drds. The desired order estimates then follows. Proof of Theorem 4.2. The theorem will be proved in two steps. In the first step, we establish the tightness of {αε (·)}, and in the second step, we characterize the limit process. l l (1) Tightness. Note that αε (t) = i=1 iI{αε (t)=i} = i=1 iI{αε (t)∈Mi } . Define χε (t) = (χ1ε (t), . . . , χlε (t)), χiε (t) = (I{αε (t)=sij } ), χε (t) = (I{αε (t)=1} , . . . , I{αε (t)=l} ). In view of the definition of αε (t) and the Cram´er–Wold theorem [2, p. 49], to prove the tightness of {αε (·)}, it suffices to derive that of {χε (·)}. For any δ > 0, any t > 0, and any s > 0 with δ > s > 0 and t + s ∈ [0, T ], owing to the Markov property, E(χε (t + t+s s) − χε (t) − t χε (u)Qε (u)du|Ft,ε ) = 0, where Ft,ε denotes the σ-algebra generated by {αε (u) : u ≤ t}. Postmultiplying by 11 and noting A(t)11 = 0 lead to E(χε (t + s) − t+s χε (t) − t χε (u)B(u)11du|Ft,ε ) = 0. Thus by (A1), E(χε (t + s)|Ft,ε ) − χε (t) = O(s). t+s This implies that E(|χε (t+s)−χε (t)|2 Ft,ε ) = E(| t χε (u)B(u)11du|2 Ft,ε ) = O(s). Consequently, for 0 < δ ≤ s, limδ→0 lim supε→0 E|χε (t + s) − χε (t)|2 = 0. The desired tightness then follows from [16, p. 47]. (2) Characterization of the limit. Since {αε (·)} is tight, by Prohorov’s theorem, we can extract a weakly convergent subsequence. Select such a sequence and still use ε as its index for notational simplicity; denote its limit by α(·). We proceed to
581
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
characterize the limit process. We need to show that for any bounded function g(·) defined on {1, . . . , l},
t
g(α(t)) −
(A.2)
B(u)g(·)(α(u))du is a martingale. 0
To establish (A.2), it suffices that for any positive integers m, and k ≤ m, any bounded functions hk (·) defined on {1, . . . , l}, and any 0 < tk ≤ t ≤ t + s, (A.3) E
m
hk (α(tk )) g(α(t + s)) − g(α(t)) −
t+s
B(u)g(·)(α(u))du
= 0.
t
k=1
l To proceed, define g(αε (s)) = i=1 g(i)I{αε (s)∈Mi } . Then, g(αε (s)) = g(αε (s)). t Thus, g(αε (t)) − 0 Qε (u)g(·)(αε (u))du is a martingale. Consequently, (A.4)
m
t+s
Ehk (αε (tk ))[g(αε (t + s)) − g(αε (t)) −
Qε (u)g(·)(αε (u))du)] = 0. t
k=1
By virtue of the weak convergence and the Skorohod representation, as ε → 0, (A.5) m m E hk (α(tk ))[g(α(t + s)) − g(α(t))]. hk (αε (tk ))[g(αε (t + s) − g(αε (t))] → E k=1
k=1
The definition of g(·) leads to A(u)g(·)(αε (u)) = 0 for all u ∈ [0, T ]. As a result, (A.6) m E hk (αε (tk )) k=1 l ∞
=
+
t+s
Qε (u)g(·)(αε (u))du
t
E
i=1 j=1 l ∞ i=1 j=1
m
hk (αε (tk ))[
k=1 m
E
t+s
ν ij (u)I{αε (u)=i} B(u)g(·)(sij )du]
t
hk (αε (tk ))[
t+s
t
k=1
[I{αε (u)=sij } − ν ij (u)I{αε (u)=i} ]B(u)g(·)(sij )du].
Using Theorem 4.1, the last term in (A.6) goes to 0 as ε → 0. By virtue of the weak convergence of αε (·) and the Skorohod representation, (A.6) yields that E (A.7)
m
hk (αε (tk ))
k=1
→E
Qε (u)g(·)(αε (u))du t
m k=1
t+s
hk (α(tk ))
t+s
B(u)g(·)(α(u))du .
t
Combining (A.4), (A.5), and (A.7), (A.3) is verified. This completes the proof. Proof of Lemma 4.3. Note that Yε (·) is a sequence of vector-valued random processes with two components. Again, using the Cr´amer–Wold theorem, since {αε (·)} is tight, to prove the tightness of {Yε (·)}, it suffices to verify the tightness of {xε (·)}.
582
G. YIN AND HANQIN ZHANG
We claim that for any δ > 0, 0 ≤ s ≤ t, and t − s ≤ δ, sup E (|xε (t) − xε (s)|2 Fs,ε ) ≤ K(t − s). 0≤s≤t≤T
First, for u ∈ [s, t], E(z ij (u, αε (u))Fs,ε ) = O(ε + exp(−κ0 (t − s)/ε)) uniformly in i and j, by (4.2), the Markov property, and Theorem 3.8. Thus
(A.8)
t l ∞ 1 E(xε (t) − xε (s) Fs,ε ) = f (sij ) √ E(z ij (u, αε (u))Fs,ε )du ε s i=1 j=1 √ 1 = √ O(ε + exp(−κ0 (t − s)/ε)) = O( ε). ε
Denote ∆ε (s, t) = E((xε (t) − xε (s))2 Fs,ε ) ⎛⎛ ⎞ ⎞2 l ∞ t 1 ⎜ ⎟ = E ⎝⎝ f (sij ) zεij (αε (u))du⎠ Fs,ε ⎠ . ε s i=1 j=1 t l ∞ d It then follows dt ∆ε (s, t) = 2ε i=1 j=1 f (sij ) s E(ζε1,ij (u, t) + ζε2,ij (u, t)Fs,ε )du, where ζε1,ij (u, t) and ζε2,ij (u, t) are defined in (A.1). Thus (d/dt)∆ε (s, t) = O(1), ∆ε (s, s) = 0. An integration then leads to ∆ε (s, t) = O(t − s). Therefore, lim lim sup sup E(E(xε (t) − xε (s))2 Fs,ε ) = 0. δ→0
ε→0
0≤t−s≤δ
Moreover, for each δ > 0, and each rational t ≥ 0, using the Chebyshev’s inequality 2 2 2 inf ε P (|xε (t)| ≤ Kt,δ ) ≥ inf ε (1 − E|xε (t)| /(Kt,δ ) ) ≥ 1 − Kt/(Kt,δ ) . Select Kt,δ > KT /δ. Then inf ε P (|xε (t)| ≤ Kt,δ ) ≥ 1 − δ. This inequality, (A.8), and the criterion [16] then yield the desired tightness. Proof of Theorem 4.4. To prove the theorem, it suffices to show that t g(x(t), α(t)) − (A.9) D(u)g(x(u), α(u))du is a martingale. 0
To verify this martingale property, we prove that for any positive integer m, any k ≤ m, any bounded and continuous function hk (·, α) for each α, and any 0 < tk ≤ t ≤ t + s, m k=1
Ehk (x(tk ), α(tk ))
g(x(t + s), α(t + s)) − g(x(t), α(t)) t+s − D(u)g(x(u), α(u))du = 0. t
To accomplish this, we begin with the process indexed by ε. (1) Using an argument similar to [29, pp. 199–200], it can be verified that the martingale problem associated with operator D(t) has a unique (in the sense of distribution) solution. (2) In what follows, we often need to carry out estimates involving xε (t) and αε (t) intertwined. To untangle them, as a preparation, let η(t, x) be a real-valued function
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
583
that is Lipschitz continuous in (t, x) ∈ R2 . Using an argument similar to that of [29, Lemma 7.4, pp. 189–192], we can show t 2 (A.10) z ij (u, αε (u))η(u, xε (u))du = 0. lim sup E ε→0 0≤t≤T 0
(3) Define an operator Dε (t) by (A.11)
⎛ ⎞ l ∞ ∂ρ(t, x, α) ∂ρ(t, x, α) 1 Dε (t)ρ(t, x, α) = f (sij )z ij (t, α)⎠ +√ ⎝ ∂t ∂x ε i=1 j=1 + Qε (t)ρ(t, x, ·)(α) for α ∈ M, ρ(·, ·, α) ∈ C 1,1 .
t Then ρ(t, xε (t), αε (t)) − 0 Dε (u)ρ(u, xε (u), αε (u))du is a martingale; see [6, Chapter 2]. l Define g(x, α) = i=1 g(x, i)I{α∈Mi } . It is easily seen that A(u)g(x, ·)(αε (u)) = 0. Thus, (A.12) 1 ∂ Dε (u)g(xε (u), αε (u)) = √ b(u, αε (u)) g(xε (u), αε (u)) + B(u)g(xε (u), ·)(αε (u)), ∂x ε l ∞ where b(u, α) = i=1 j=1 f (sij )z ij (u, α). To obtain the desired limit, we use the techniques of perturbed test function methods [16], which require us to introduce a perturbation that is small in magnitude and that results in desired cancellations. To this end, define the perturbation by g(t, x, α) such that g(·, x, α) is uniform Lipschitz in t; both g(·) and (∂/∂x) g (·) are bounded; (∂/∂x) g (·, ·, α) is Lipschitz in (t, x) such that it is the solution of ∂g(x, α) A(t) g (t, x, ·)(α) = −b(t, α) . ∂x
(A.13)
It can be shown that such a function exists similar to [29, Remark 7.16]. Next, define the perturbed test function gε (·) by √ gε (t, x, α) = g(x, α) + ε (A.14) g (t, x, α). t Then using the definition of Dε (t), gε (t, xε (t), αε (t)) − 0 Dε (u)gε (u, xε (u), αε (u))du is a martingale. Consequently, m Ehk (xε (tk ), αε (tk )) gε (t + s, xε (t + s), αε (t + s)) − gε (t, xε (t), αε (t)) k=1 t+s
−
Dε (u)gε (u, xε (u), αε (u))du = 0. t
By the weak convergence of Yε (·) to Y (·), the Skorohod representation, and (A.14), m
lim
ε→0
Ehk (xε (tk ), αε (tk ))[gε (t + s, xε (t + s), αε (t + s)) − gε (t, xε (t), αε (t))]
k=1 m
=
k=1
Ehk (x(tk ), α(tk ))[g(x(t + s), α(t + s)) − g(x(t), α(t))].
584
G. YIN AND HANQIN ZHANG
Note that from (A.11), √
εDε (u) g (u, xε (u), αε (u)) =
l ∞
f (sij )z ij (u, αε (u))
i=1 j=1
(A.15)
∂ g (u, xε (u), αε (u)) ∂x
√ ∂ g (u, xε (u), αε (u)) + ε ∂u √ 1 + εB(u) g (u, xε (u), ·)(αε (u)). g (u, xε (u), ·)(αε (u)) + √ A(u) ε
It follows from (A.12), (A.13), and (A.15), upon cancellation, that (A.16) l ∞
Dε (u)gε (u, xε (u), αε (u)) =
f (sij )z ij (u, αε (u))
i=1 j=1
∂ g (u, xε (u), αε (u)) ∂x
+ B(u)g(xε (u), ·)(αε (u)) √ ∂ g (u, xε (u), αε (u)) √ + ε g (u, xε (u), ·)(αε (u)). + εB(u) ∂u Thus again, by the weak convergence and the Skorohod representation, t+s m lim Dε (u)gε (u, xε (u), αε (u))du Ehk (xε (tk ), αε (tk )) ε→0 t k=1 m l ∞ t+s ∂ g (u, xε (u), αε (u)) = lim Ehk (xε (tk ), αε (tk )) f (sij )z ij (u, αε (u)) du ε→0 ∂x t i=1 j=1 k=1 t+s B(u)g(xε (u), ·)(αε (u))du . + t
By virtue of the weak convergence of Yε (·) to Y (·), the Skorohod representation, t and the boundedness of g(·), (∂/∂x) g (·), we have 0 B(u)g(xε (u), ·)(αε (u))du → t B(u)g(x(u), ·)(α(u))du. 0 Define ∂ g (t, x, α) b(t, x, α) = b(t, α) . ∂x
(A.17) Then
t
b(u, xε (u), αε (u))du = (A.18)
0
+
t l ∞
t l ∞
b(u, xε (u), sij )ν ij (u)I{αε (u)=i} du
0 i=1 j=1
[I{αε (u)=sij } − ν ij (u)I{αε (u)=i} ]b(u, xε (u), sij )du.
0 i=1 j=1
t By (A.10), the last term in (A.18) goes to 0 in mean squares, so 0 b(u, xε (u), αε (u))du t ∞ converges to 0 b(u, x(u), α(u))du, where b(u, x, i) = j=1 ν ij (u)b(u, x, sij ). Combining the estimates obtained so far, we arrive at t+s Ehk (xε (tk ), αε (tk )) Dε (t)gε (u, xε (u), αε (u))du t t+s t+s b(u, x(u), α(u))du + B(u)g(x(u), ·)(α(u))du . → Ehk (x(tk ), α(tk )) t
t
585
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
(4) In view of (A.13), for i = 1, . . . , l, ⎞ ⎞ ⎛ b(t, si1 ) ∂g(x,i) g(t, x, si1 ) ∂x ⎟ ⎜ Ai (t) ⎝ g(t, x, si2 ) ⎠ = ⎝ b(t, si2 ) ∂g(x,i) ⎠, ∂x · · : : ⎛
which has a unique solution. This implies that g(t, x, sij ) is a function of (∂/∂x)g(x, i). g (t, x, sij ). Thus b(t, x, sij ) is a In view of (A.17), b(t, x, sij ) is a function of (∂/∂x) function of (∂ 2 /∂x2 )g(x, i). Denote b(t, x, i) = (1/2)a(t, i)(∂ 2 /∂x2 )g(x, i), where a(t, i) is an appropriate function. Using an argument similar to that [29, pp. 200–203], it can be shown that a(t, i) ≥ 0. Thus, a(t, i) can be written as a(t, i) = σ 2 (t, i), where for each i = 1, 2, . . . , l, ∞ ∞ ∞ ik 2 (A.19) f (sik )f (sim ) ν (t) ψ0km (i, u, t)du σ (t, i) = 0
k=1 m=1
+ν
im
(t)
∞
ψ0mk (i, u, t)du
,
0
and ψ0km (i, s, t) is the (k, m)th entry of Ψ0 (i, s, t) given by Ψ0 (i, s, t) = (I − (ν i (s), ν i (s), . . .) ) exp(Ai (s)t). This specifies the covariance structure of the limit process and establishes the desired martingale property and hence the theorem follows. Acknowledgment. We thank the editor and the reviewers for detailed comments and suggestions on an early version of the manuscript, which have led to much improvement of the paper. REFERENCES [1] E. Altman, K. E. Avrachenkov, and R. Nunez-Queija, Pertrubation analysis for denumerable Markov chains with applications to queueing models, Adv. Appl. Probab., 36 (2004), pp. 839–853. [2] P. Billingsley, Convergence of Probability Measures, Wiley, New York, 1968. [3] G. Bitran and D. Tirupati, Multiproduct queueing networks with deterministic routing: Decomposition approach and the notion of interference, Management Sci., 34 (1988), pp. 75– 100. [4] N. N. Bogoliubov and Y. A. Mitropolskii, Asymptotic Methods in the Theory of Nonlinear Oscillator, Gordon and Breach, New York, 1961. [5] J. Daigle and J. Langford, Models for analysis of packet voice communication systems, IEEE J. Sel. Areas Commun., SAC-4 (1986), pp. 847-855. [6] M. H. A. Davis, Markov Models and Optimization, Chapman & Hall, London, 1993. [7] J. L. Doob, Stochastic Processes, Wiley Classic Library Edition, Wiley, New York, 1990. [8] S. G. Eick, W. A. Massey, and W. Whitt, The physics of the Mt /G/∞ queue, Oper. Res., 41 (1993), pp. 731–742. [9] V. Hutson and J. S. Pym, Applications of Functional Analysis and Operator Theory, Academic Press, London, 1980. [10] A. M. Il’in, A. S. Kalashnikov, and O. A. Oleinik, Linear equations of the second order of parabolic type, Russian Math. Surveys, 17 (1962), pp. 1–143. [11] A. M. Il’in, R. Z. Khasminskii, and G. Yin, Asymptotic expansions of solutions of integrodifferential equations for transition densities of singularly perturbed switching diffusions: Rapid switchings, J. Math. Anal. Appl., 238 (1999), pp. 516–539.
586
G. YIN AND HANQIN ZHANG
[12] J. Keller, Time-dependent queues, SIAM Rev., 24 (1982), pp. 401–412. [13] R. Z. Khasminskii, G. Yin, and Q. Zhang, Asymptotic expansions of singularly perturbed systems involving rapidly fluctuating Markov chains, SIAM J. Appl. Math., 56 (1996), pp. 277–293. [14] C. Knessl, B. J. Matkowsky, Z. Schuss, and C. Tier, Asymptotic analysis of a statedependent M/G/1 queueing system, SIAM J. Appl. Math., 46 (1986), pp. 483–505. [15] C. Knessl and J. A. Morrison, Heavy traffic analysis of the sojourn time in tandem queues with overtaking, SIAM J. Appl. Math., 51 (1991), pp. 1740–1763. [16] H. J. Kushner, Approximation and Weak Convergence Methods for Random Processes, with Applications to Stochastic Systems Theory, MIT Press, Cambridge, MA, 1984. [17] G. E. Ladas and V. Lakshmikantham, Differential Equations in Abstract Spaces, Academic Press, New York, 1972. [18] W. A. Massey, Asymptotic analysis of the time dependent M/M/1 queue, Math. Oper. Res., 10 (1985), pp. 305–327. [19] W. A. Massey and W. Whitt, Unstable asymptotics for nonstationary queues, Math. Oper. Res., 19 (1994), pp. 267–291. [20] W. A. Massey and W. Whitt, Uniform acceleration expansions for Markov chains with timevarying rates, Ann. Appl. Probab., 8 (1998), pp. 1130–1155. [21] S. P. Meyn and R. L. Tweedie, Computable bounds for geometric convergence rates of Markov chains, Ann. Appl. Probab., 4 (1994), pp. 981–1021. [22] M. F. Neuts, Matrix-Geometric Solutions in Stochastic Models, Johns Hopkins University Press, Baltimore, 1981. [23] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Springer-Verlag, New York, 1983. [24] M. I. Reiman, Asymptotically exact decomposition approximations for open queueing networks, Oper. Res. Lett., 9 (1990), pp. 363–370. [25] A. B. Vasil’eava and V. F. Butuzov, Asymptotic Methods in Singular Perturbations Theory, Vysshaya Shkola, Moscow, 1990 (in Russian). [26] W. Whitt, The queueing network analyzer, Bell Syst. Tech. J., 62 (1983), pp. 2779–2815. [27] W. Whitt, Stochastic-Process Limits, Springer-Verlag, New York, 2001. [28] G. Yin and H. Zhang, Countable-state-space Markov chains with two-time scales and applications to queueing systems, Adv. Appl. Probab., 34 (2002), pp. 662–688. [29] G. Yin and Q. Zhang, Continuous-Time Markov Chains and Applications: A Singular Perturbation Approach, Springer-Verlag, New York, 1998. [30] G. Yin, Q. Zhang, and G. Badowski, Asymptotic properties of a singularly perturbed Markov chain with inclusion of transient states, Ann. Appl. Probab., 10 (2000), pp. 549–572.
SIAM J. APPL. MATH. Vol. 65, No. 2, pp. 567–586
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS TO QUASI-BIRTH-DEATH QUEUES∗ G. YIN† AND HANQIN ZHANG‡ Abstract. Aiming at reduction of complexity, this work is concerned with two-time-scale Markov chains and applications to quasi-birth-death queues. Asymptotic expansions of probability vectors are constructed and justified. Lumping all states of the Markov chain in each subspace into a single state, an aggregated process is shown to converge to a continuous-time Markov chain whose generator is an average with respect to the stationary measures. Then a suitably scaled sequence is shown to converge to a switching diffusion process. Extensions of the results are presented together with examples of quasi-birth-death queues. Key words. Markov chain, singular perturbation, countable state space, asymptotic expansion, occupation measure, aggregation, switching diffusion, quasi-birth-death queue AMS subject classifications. 34E05, 60J27, 60F05 DOI. 10.1137/S003613990139756X
1. Introduction. Much effort has gone into evaluating performance measures of queueing systems in the past two decades; see, for example, [3, 22, 24, 26] and the references therein. Since exact solutions are difficult to obtain in many queueing problems, we are content with approximate solutions. To treat large-scale Markovian queueing systems, one of the most popular methods is decomposition, which consists of breaking the underlying network into smaller pieces (e.g., one station in each piece); see [3, 24, 26] among others. Although one often uses time-homogeneous Markovian models for approximating the actual systems, many queueing systems in real life are nonstationary (time-dependent); for example, the arrival and service rates in the systems are time-varying. Recently, time-inhomogeneous Markovian queueing networks have been widely used to model telecommunication systems; see [5]. Developing computational methods and approximation techniques for these quantities involved in time-inhomogeneous queueing problems has long been regarded as a challenging task (see [14] and [20]); see also [12, 18] and the references therein for earlier effort in this direction. The motivation of our study stems from the recent advances in understanding asymptotic properties of singularly perturbed Markov chains aimed at reducing of complexity. We began our investigation in [13]. By combining matched asymptotic expansions and stochastic analysis, our effort was subsequently extended to treat more complex models and applications in control and optimization; see [29, 30]. These results are about Markov chains with finite-state spaces, and they are applicable to queues with a finite number of waiting rooms or finite capacity; see also [20] for a related work. Continuing our effort initiated in [28], using a model similar to that ∗ Received by the editors November 6, 2001; accepted for publication (in revised form) May 21, 2004; published electronically January 5, 2005. http://www.siam.org/journals/siap/65-2/39756.html † Department of Mathematics, Wayne State University, Detroit, MI 48202 (gyin@math. wayne.edu). The research of this author was supported in part by National Science Foundation grant DMS-0304928 and in part by the Wayne State University Research Enhancement Program. ‡ Institute of Applied Mathematics, Academy of Mathematics and Systems Sciences, Academia Sinica, Beijing, 100080, China ([email protected]). The research of this author was supported in part by a Distinguished Young Investigator grant from the National Natural Sciences Foundation of China and a grant from the Hundred Talents Program of the Chinese Academy of Sciences.
567
568
G. YIN AND HANQIN ZHANG
presented in [26] as our starting point, we consider queueing systems in which the generator of the queue length process includes both a fast varying part and a slowly changing part reflecting both strong and weak interactions of states belonging to different irreducible classes. Compared with the methods used in [3, 24, 26], we adopt a singular perturbation approach via time-scale separation to model the different transition intensities. In [28], treating countable-state-space Markov chains, one of the main assumptions used is that the underlying Markov chain is irreducible. In the context of queueing theory, it basically corresponds to the consideration of a single queue. The goal of this paper is to treat multistation queueing networks and to develop decomposition and aggregation methods for queueing network problems. Motivated by queueing networks involving quasi-birth-death processes, we start with decomposition by splitting the entire state space into a number of subspaces. Usually, the transitions within each subspace are much more intensive and frequent than those among different subspaces. We introduce a small parameter ε > 0 to highlight the different rates of transition intensities. Then we proceed to derive the asymptotic expansions of the probability distribution of the queue length. By lumping all the states in each subspace into a single state, we obtain a sequence of aggregated processes. We demonstrate that this sequence converges to a Markov chain. Furthermore, we obtain limit results for suitably scaled sequences. The rest of the paper is arranged as follows. Section 2 begins with the precise formulation of the problem. Section 3 provides asymptotic expansions. A sequence of occupation measures is defined in section 4; its probabilistic properties, including mean square estimate, aggregation, and switching diffusion limit of a suitably scaled sequence, are examined. Section 5 presents extensions of results and queueing examples. Finally, an appendix is provided to include some technical complements. 2. Formulation. Working with a finite time horizon t ∈ [0, T ] for some T > 0, our focus is on time-inhomogeneous Markov chains. Suppose that β(t) is a continuoustime Markov chain with countable state space N = {1, 2, . . .}. An infinite dimensional matrix-valued function Q(t) = (q ij (t)) defined on [0, T ] is a generator of the Markov chain β(t) if q ij (·) is Borel measurable and bounded for each i, j ∈ N, q ij (t) ≥ 0 for ∞ all i = j, j=1 q ij (t) = 0 for all i ∈ N, and for any bounded and Borel-measurable t function g(·) defined on N, g(β(t)) − 0 Q(s) g (·)(β(s))ds is a martingale, where (2.1)
Q(t) g (·)(i) =
∞
q ij (t) g (j) for each i ∈ N.
j=1
Throughout the paper, we use K to denote a generic positive constant. The conventions K + K = K and KK = K are used for simplicity. We use 11m0 and 11 to denote an m0 -dimensional and infinite dimensional column vector with all components being 1. For a vector z and a matrix H, we use z and H to denote their transposes, and we use z i and hij to denote the ith component of z and the ijth entry of H = (hij ), respectively. For a given matrix H = (hij )∞×∞ with infinite columns and infinite · rows, Ha is an augmented matrix given by Ha = (H : 11). In addition, we use a subscript to index a sequence. A Markov chain β(t) or its generator Q(t) is weakly irreducible, if the system · of equations g(t)Qa (t) = (0 : 1) has a unique solution g(t) = (g i (t)) with g i (t) ≥ 0 for each i ∈ N, where 0 = (0, 0, . . .) is an infinite dimensional 0 vector. The unique nonnegative solution is termed a quasi-stationary distribution. (An equivalent way to write the system of equations in the above definition is g(t)Q(t) = 0, g(t)11 = 1.) The
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
569
definition is an extension of the usual notion of irreducibility, and the weak irreducibility given in [13] for finite-state Markov chains. Compared with the usual definition of irreducibility, it deals with time-varying generators and allows some components of the quasi-stationary distribution to be 0. Suppose that the Markovian queueing network has n0 (n0 < ∞) interconnected stations. Consider a vector-valued queue-length process taking values in Nn0 = N × N × · · · × Nn0 (an n0 -fold product). We use α(·), a continuous-time Markov chain, to model the queue length of the queueing network. Suppose that Nn0 can be divided into l subsets. Within each subset, the transitions (such as arrivals and departures of customers, etc.) take place an order of magnitude more frequent than that of among different subsets. To highlight this contrast, we introduce a small parameter ε > 0. Thus the continuous-time Markov chain is ε-dependent, i.e., α(·) = αε (·). For t ∈ [0, T ], assume that the generator of the Markov chain αε (t) is given by (2.2)
Qε (t) =
A(t) + B(t), with A(t) = diag(A1 (t), . . . , Al (t)), ε
where A(t), B(t), and Ai (t) (for i = 1, . . . , l) are all generators of certain countablestate space Markov chains, and diag(A1 (t), . . . , Al (t)) denotes a block diagonal matrix, each of the entries of which has appropriate dimension. Note that B(t) is a generator and there is no need to assume it has the same diagonal matrix form as that of A(t). There is a certain hierarchy in the underlying network. Within each subset, one observes detailed variations of the networks such as arrivals and services of customers, etc., whereas at an upper system management level, instead of these variations, one observes the transitions among different subsets. Thus from an upper management point of view, by lumping all the states in each subspace into a single one, the system may be regarded as if it were a queue with finite waiting rooms (l rooms). Nevertheless, the aggregated process is generally non-Markovian. Fortunately, as will be shown, the aggregated process converges weakly to a limit process that is a finite-state Markovian queue. The significance of such a result is that in lieu of examining the detailed variations, one can study the aggregated process. To reflect the decomposition in the network, write the state space of the queueing network as M = M1 ∪ M2 · · · ∪ Ml , where Mi = {si1 , si2 , . . . , } for i = 1, . . . , l. Following our approach for singularly perturbed Markov chains with finite-state spaces, we construct matched asymptotic expansions of the probability vector (2.3)
21 l1 ij pε (t) = (p11 ε (t), . . . , pε (t), . . . , . . . , pε (t), . . .), pε (t) = P (αε (t) = sij ).
For future use, partition an infinite dimensional vector v as v = (v 1 , . . . , v l ) with v i = (v i1 , v i2 , . . .). That is, v i is an infinite dimensional vector corresponding to the subspace Mi . Since we are dealing with countable state space Markov chains, we work with an infinite dimensional vector space. It is naturalto consider the spaces ∞ l 1 = {(v 1 , . . . , v l ) : 1 ≤ i ≤ l, v ik ∈ R for each k ∈ N, and i=1 k=1 |v ik | < ∞}, ∞ = {(v 1 , . . . , v l ): 1 ≤ i ≤ l, v ik ∈ R for each k ∈ N, and sup1≤i≤l sup1≤k 0, exp(A(0)τ ) − diag(11ν 1 (0), . . . , 11ν l (0)) ≤ K exp(−κτ ), (3.1) ∞ where ν i (t) = (ν i1 (t), ν i2 (t), . . .) is the quasi-stationary distribution corresponding to the generator Ai (t). Remark 3.1. From (A1) and the definition of weak irreducibility, ν i (t)Aia (t) = · · (0 : 1) has a unique solution. Furthermore, ν i (t) = (0 : 1)(Aia (t)) (Aia (t)(Aia (t)) )−1 . (A2) is a Doeblin-type condition. A condition in a slightly different form is given in [7, p. 192]. We will obtain the asymptotic expansions in two steps. The first step is a formal construction in which we (see [4, 10, 25] and [11, 13] among others) find the outer expansions and the initial layer corrections. The second step involves validating the formal expansions and deriving the desired error estimates. 3.1. Formal expansions. Weseek matched asymptotic expansions of the form n n Oε,n (t) + Iε,n (t) = i=0 εi φi (t) + i=0 εi ψi (t/ε). For technical reasons, which will become clear in what follows, to justify the validity of the expansions, we also need to compute φn+1 (t) and ψn+1 (t/ε). To make sure that the matching condition is satisfied, we use Oε,n+1 (0) + Iε,n+1 (0) = p(0) or, more specifically, φ0 (0) + ψ0 (0) = p(0), φk (0) + ψk (0) = 0 for 1 ≤ k ≤ n + 1. Our approach is based on constructions of the sequences {φi (t)} and {ψi (t/ε)}. Substituting Oε,n+1 (t) into (2.4) and comparing coefficients of powers of εk , we obtain (3.2)
φ0 (t)A(t) = 0,
φk (t)A(t) = φ˙ k−1 (t) − φk−1 (t)B(t), 1 ≤ k ≤ n + 1.
The outer expansions give us satisfactory approximation for t > 0 away from an initial layer of the order O(ε), but it does not satisfy the initial condition and breaks down for sufficiently small t. To compensate, introduce a fast time variable τ = t/ε. By the Lipschitz continuity given in (A1), taking Taylor expansions of A(ετ ) and B(ετ ) i i i i n+1 n A(0) B(0) about 0 yields A(ετ ) = i=0 (ετi!) d dt + O((ετ )n+2 ), εB(ετ ) = i=0 ε (ετi!) d dt + i i n+2 O((ετ ) ). Substituting Iε,n+1 (t) into (2.4), using the above Taylor expansions of A(t) and B(t), and comparing powers of εi , we obtain the equations satisfied by the initial layer terms. We have ψ0 (0) = p(0) − φ0 (0), ψk (0) = −φk (0), and
(3.3)
dψ0 (τ ) dψk (τ ) = ψ0 (τ )A(0), = ψk (τ )A(0) + rk (τ ), 1 ≤ k ≤ n + 1, dτ dτ k−1 τ i+1 di+1 A(0) τ i di B(0) , 1 ≤ k ≤ n + 1. ψk−i−1 (τ ) + rk (τ ) = (i + 1)! dti+1 i! dti i=0
In [28], in which A(t) consists of only one block, the outer expansions and initial layers can be obtained separately. Here, the constructions of {φk (t)} have to be done
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
571
in conjunction with the initial layer corrections. The equations in (3.3) together with the initial data are known as abstract Cauchy problems; see [17, p. 21]. It is well known that {exp(A(0)τ )} is a uniformly continuous semigroup (see [17, 23]). With ψk (0) to be determined, the representation of the solutions of (3.3) is given by (3.4)
ψ0 (τ ) = (p(0) − φ0 (0)) exp(A(0)τ ), τ
rk (s) exp(A(0)(τ − s))ds, 1 ≤ k ≤ n + 1.
ψk (τ ) = ψk (0) exp(A(0)τ ) + 0
Step 1. Determine φ0 (t) and ψ0 (τ ). We begin with the first equation in (3.2). Using the partitioned vector form introduced right after (2.3), we obtain are not uniquely solvable φi0 (t)Ai (t) = 0 for each i = 1, . . . , l. These equations ∞ i since Ai (t) have a 0 eigenvalue. However, by attaching j=1 φij 0 (t) = θ0 (t) to the equations, the resulting system of equations has a unique solution thanks to the weak irreducibility of Ai (t). Thus, φi0 (t) must be proportional to ν i (t), the quasi-stationary distribution corresponding to Ai (t). That is, φi0 (t) = θ0i (t)ν i (t), where θ0i (t) ∈ R is to be determined. Define θ0 (t) = (θ01 (t), . . . , θ0l (t)) ∈ R1×l , 11 = diag(11, 11, . . . , 11), ν(t) = diag(ν 1 (t), . . . , ν l (t)). (3.5) It is immediate that A(t) is orthogonal to 11 (i.e., A(t)11 = 0) and φ0 (t) = θ0 (t)ν(t). For the equation with k = 1 in (3.2), multiplying from the right by 11 leads to (3.6)
θ˙0 (t) = θ0 (t)B(t),
B(t) = ν(t)B(t)11,
so B(t) is an average of B(t) with respect to ν 1 (t), . . . , ν l (t). Note that B(t) is an l × l matrix-valued function. Note also that (3.6) is a linear system of differential equations and that it has a unique solution for each initial condition. Choose θ0 (0) = p(0)11. Then the solution of (3.6) is uniquely determined. The φi0 (t) has the interpretation of total probability. Concerning the initial layer correction ψ0 (τ ), the solution is given by the first equation in (3.4) and is uniquely solved. We claim that ||ψ0 (τ )||∞ ≤ K exp(−κτ ) for some κ > 0 and K > 0. To prove this, note that 11ν(0) = diag(11ν 1 (0), . . . , 11ν l (0)) and (3.7)
ψ0 (0)11ν(0) = [p(0) − φ0 (0)]11ν(0) = (0, 0, . . .),
where φ0 (0) = θ0 (0)ν(0), φ0 (0)11 = θ0 (0), and θ0 (0) = p(0)11 are used. By virtue of the orthogonality and (A2), ||ψ0 (τ )||∞ = ||ψ0 (0)(exp(A(0)τ ) − 11ν(0))||∞ ≤ K exp(−κτ ). To determine φi (t) and ψi (t) for i ≥ 1, we need the following lemma. Lemma 3.2. Suppose that Q(t) is a generator of a countable-state-space Markov chain such that Q(t) is weakly irreducible for each t ∈ [0, T ] and (dn+1 /dtn+1 )Q(·) ∈ −1 (t) exists = Qa (t)Q (t). Then for k = 1, . . . , n+1, (dk /dtk )Q Lip[0, T ]. Denote Q(t) a n+1−k and belongs to C (the class of functions that are (n + 1 − k)-times continuously differentiable). Q −1 (t) = I, differentiating both sides of the above equation Proof. Since Q(t) leads to (3.8)
−1 (t) dQ −1 (t) dQ(t) Q −1 (t), = −Q dt dt
−1 (·) is differentiable. Repeatedly differentiating equation (3.8) yields the desired so Q result.
572
G. YIN AND HANQIN ZHANG
Step 2. Determine φ1 (t) and ψ1 (τ ). Lemma 3.2 and Remark 3.1 imply that φ0 (·) is differentiable and is a function of class C n+1 . We proceed to determine φ1 (t) and ψ1 (τ ). Note that the equation with k = 1 in (3.2) is a nonhomogeneous equation def whose right-hand side φ0 (t) = φ˙ 0 (t) − φ0 (t)B(t) is a known function since φ0 (t) has been found. Using (3.6), [φ˙ 0 (t)−φ0 (t)B(t)]11 = θ˙0 (t)−θ0 (t)ν(t)B(t)11 = (0, 0, . . .), and the Fredholm alternative yields that the equation with k = 1 in (3.2) has a particular solution φ1,p (t) being orthogonal to 11. Assume the solution of φ1 (t)A(t) = φ0 (t) to be of the form φ1 (t) = θ1 (t)ν(t)+φ1,p (t). Since φ1,p (t) is orthogonal to 11, postmultiplying the equation with k = 2 in (3.2) by 11 leads to θ˙1 (t) = θ1 (t)B(t) + φ1,p (t)B(t)11.
(3.9)
Once the initial condition is specified, (3.9) is uniquely solved. The initial condition θ1 (0) has to come from the initial layer correction term. With the selection of ψ1 (0) = −φ1 (0), by (3.4) with k = 1, the unique solution is given by τ ψ0 (s) exp(A(0)s)B(0) exp(A(0)(τ − s))ds ψ1 (τ ) = ψ1 (0) exp(A(0)τ ) + 0 τ (3.10) dA(0) exp(A(0)(τ − s))ds. + sψ0 (s) exp(A(0)s) dt 0 τ ∞ By the exponential decay of ψ0 (τ ), 0 ψ0 (s) exp(A(0)s)ds∞ = 0 ||ψ0 (s)||∞ ds < def ∞ ∞. Denote ψ 0 = 0 ψ0 (s) exp(A(0)s)dsB(0). Then from (A2),
lim ψ1 (0) exp(A(0)τ ) = lim ψ1 (0)[exp(A(0)τ ) − 11ν(0)] + ψ1 (0)11ν(0) τ →∞ τ →∞ = ψ1 (0)11ν(0), τ ψ0 (s) exp(A(0)s)B(0) exp(A(0)(τ − s))ds = ψ 0 11ν(0). lim τ →∞
0
Note that (3.11)
A(t)11 = 0,
dk A(t) dk A(t)11 = 0, 1 1 = dtk dtk
k = 1, . . . , n + 1.
Using the orthogonality (see (3.7) and (3.11)), we obtain τ dA(0) exp(A(0)(τ − s))ds ≤ Kτ 2 exp(−κτ ). sψ0 (s) exp(A(0)s) (3.12) dt 0 ∞ By demanding limτ →∞ ψ1 (τ ) = 0, taking the limit as τ → ∞ in (3.10), and using the estimates above, we arrive at (3.13)
ψ1 (0)11ν(0) + ψ 0 11ν(0) = 0.
An important observation indicates that there are only l unknowns in (3.13). Using the notation of partitioned vector given after (2.3), ψ1 (0) = (ψ11 (0), ψ12 (0), . . . , ψ1l (0)) 1 2 l i and ψ 0 = (ψ 0 , ψ 0 , . . . , ψ 0 ), the solution is given by ψ1i (0)11 = −ψ 0 11, i = 1, . . . , l. i Since ψ1 (0) must be chosen so that φ1 (0) + ψ1 (0) = 0, we obtain θ1i (0) = ψ 0 11, i = 1, . . . , l. Thus both φ1 (t) and ψ1 (τ ) have been determined. We claim that ψ1 (τ )
573
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
decays exponentially fast. By adding and subtracting appropriate terms, ∞ ψ0 (s) exp(A(0)s)dsB(0)11ν(0) ψ1 (τ ) = ψ1 (0)[exp(A(0)τ ) − 11ν(0)] − τ ∞ + ψ1 (0)11ν(0) + ψ0 (s) exp(A(0)s)dsB(0)11ν(0) 0 (3.14) τ + ψ0 (s) exp(A(0)s)B(0)[exp(A(0)(τ − s)) − 11ν(0)]ds 0 τ dA(0) exp(A(0)(τ − s))ds. + sψ0 (s) exp(A(0)s) dt 0 By (A2),
(3.15)
ψ1 (0)[exp(A(0)τ ) − 11ν(0)] ≤ K exp(−κτ ), and ∞ τ ψ0 (s) exp(A(0)s)B(0)[exp(A(0)(τ − s)) − 11ν(0)]ds 0
≤ K(1 + τ ) exp(−κτ ).
∞
It then follows from the above estimates, ||ψ1 (τ )||∞ ≤ K(1 + τ + τ 2 ) exp(−κτ ) ≤ K exp(−κ0 τ ), for some 0 < κ0 < κ. Step 3. Determine φk (t) and ψk (τ ), for 1 < k ≤ n + 1. We apply the same method for finding φ1 (t) and ψ1 (τ ) to determine φk (t) and ψk (τ ). For each i = 1, . . . , l, assume φik (t) is of the form φik (t) = θki (t)ν i (t) + φik,p (t), where φk,p (t) = (φ1k,p (t), . . . , φlk,p (t)) is a particular solution of the equation φk (t)A(t) = φ˙ k−1 (t) − φk−1 (t)B(t), and φk,p (t) is orthogonal to 11. We proceed to determine θk (t) = (θki (t)) ∈ R1×l . By substitution, similar to (3.9), it is easily seen that θk (t) satisfies the differential equation θ˙k (t) = θk (t)B(t) + φk,p (t)B(t)11.
(3.16)
To determine the initial condition θk (0), the definition of rk (τ ) in (3.4) gives us ψk (τ ) = ψk (0) exp(A(0)τ ) + 0
(3.17) ×
τ k−1
ψk−j−1 (s)
j=0 j j
τ j+1 dj+1 A(0) τ d B(0) exp(A(0)(τ − s))ds. + (j + 1)! dtj+1 j! dtj
Using the same techniques as that for θ1 (0), we obtain θk (0) and uniquely determine both φk (t) and ψk (τ ). We then verify the exponential decay property of ψk (τ ). The procedure is the same as for that of ψ1 (τ ); we record the result as follows. Theorem 3.3. Assume (A1) and (A2). Then the following assertions hold: (a) The sequence {φk (t)} can be constructed by solving the system of equations (3.18)
def φk (t)A(t) = φ˙ k−1 (t) − φk−1 (t)B(t) = φk−1 (t)
for k = 1, 2, . . . , n. The solution is φk (t) = θk (t)ν(t) + φk,p (t), where φk,p (t) is a particular solution of (3.18), which is orthogonal to 11, and θk (t) satisfies θ˙k (t) = θk (t)B(t) + φk,p (t)B(t)11. j k−1 ∞ j def A(0) )11ν(0) = (b) Find ψki (0)11 from ψk (0)11ν(0) = −( j=0 0 sj! ψk−j−1 (s)ds d dt j i −ψ k−1 11ν(0). Choose θki (0) = −ψki (0)11 = ψ k−1 11 for i = 1, . . . , l. Choose ψk (0) = −φk (0).
574
G. YIN AND HANQIN ZHANG
(c) For k = 0, 1, . . . , n+1, θk (·) and φk (·) are (n+1−k)-times continuously differentiable on [0, T ]; there exists a κ0 > 0 such that ||ψk (τ )||∞ ≤ K exp(−κ0 τ ). 3.2. Asymptotic justification. In this section, we obtain the desired error bounds and show that the asymptotic expansions hold uniformly in t ∈ [0, T ]. Define an operator Lε as (3.19)
Lε f (t) = ε
df (t) − f (t)A(t) − εf (t)B(t) dt
for a suitable smooth function f (·). First let us establish a lemma. Lemma 3.4. Consider Lε vε (t) = ∆ε,k (t), vε (0) = 0, with supt∈[0,T ] ||∆ε,k (t)||∞ = O(εk+1 ) for some k with 0 ≤ k ≤ n + 1. Then supt∈[0,T ] ||vε (t)||∞ = O(εk ). Proof. Note that the initial value problem given above is a time-dependent abstract Cauchy problem or an evolution equation. The solution is given by vε (t) = 1 t ε 0 ∆ε,k (s)Xε (t, s)ds, where Xε (t, s) is a fundamental solution or an evolution operator. In fact (see Ladas and Lakshmikantham [17, p. 56]), Xε (t, s) is an operator-valued function with values in L(1 ); the space of bounded linear operators, defined on 1 and strongly continuous in t, s for 0 ≤ s ≤ t ≤ T such that (∂/∂t)Xε (t, s) exists in strong topology of 1 ; (∂/∂t)Xε (t, s) ∈ L(1 ) for 0 ≤ s ≤ t ≤ T and (∂/∂t)Xε (t, s) is strongly continuous in t; the range of Xε (t, s) is in the domain of Qε (t); and it satisfies the homogeneous problem ε ∂Xε∂t(t,s) − Xε (t, s)A(t) − εXε (t, s)B(t) = 0, Xε (s, s) = I, where the initial value I is the infinite dimensional identity matrix. Since it represents transition probabilities, ||Xε (t, s)||∞ is bounded uniformly in ε for all t, s ∈ [0, T ]. t k Therefore, we have supt∈[0,T ] ||vε (t)||∞ ≤ K ε supt∈[0,T ] 0 ||∆ε,k (s)||∞ ds ≤ Kε . The lemma is obtained. For each k = 1, . . . , n + 1, define a vector-valued error eε,k (t) (3.20)
eε,k (t) = pε (t) −
k i=0
εi φi (t) −
k
εi ψi (t/ε),
i=0
where pε (·) is the solution of (2.4), φi (·) and ψi (·) are the outer expansions and initial layer corrections, respectively. We must show that eε,n (t) = O(εn+1 ). For this purpose, we first derive a lemma whose proof is similar in spirit to the corresponding results for weakly irreducible generators [28] and is thus omitted. Lemma 3.5. Under (A1) and (A2), for k = 1, . . . , n+1, supt∈[0,T ] ||Lε eε,k (t)||∞ = k+1 O(ε ), and hence supt∈[0,T ] ||eε,k (t)||∞ = O(εk ). Remark 3.6. Using Lemma 3.5, with k = 1, we obtain supt∈[0,T ] ||eε,1 (t)||∞ = O(ε). However, eε,1 (t) = eε,0 (t) − εφ1 (t) − εψ1 (t/ε). In view of the boundedness of φ1 (t) and ψ1 (t/ε), εφ1 (t) + εψ1 (t/ε) = O(ε). Thus, eε,0 (t) = O(ε) uniformly in t ∈ [0, T ]. We can proceed inductively. By virtue of Lemma 3.5, with k = n + 1, supt∈[0,T ] ||εε,n+1 ||∞ = O(εn+1 ). Going back one step and using eε,n+1 (t) = eε,n (t) + O(εn+1 ), similar to the case of k = 1, we obtain supt∈[0,T ] ||eε,n (t)||∞ = O(εn+1 ). We summarize the discussion thus far into the following theorem. This gives us the desired approximation results. Not only is the convergence of P (αε (t) = i) obtained, but also the rate of convergence is derived. Furthermore, a full asymptotic series is obtained. Theorem 3.7. Under (A1) and (A2), we can construct two sequences {φk (t)}nk=0 and {ψk (t/ε)}nk=0 by Theorem 3.3 such that φk (t) ∈ C n+1−k and ψk (t/ε) decay exponentially fast. Moreover with eε,n (t) defined by (3.20), supt∈[0,T ] ||eε,n (t)||∞ = O(εn+1 ).
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
575
Using the same techniques, we can also obtain asymptotic expansions of transition probability matrices. Since the proofs are essentially the same, we will state only the results and omit the detailed argument. Let Pε (t0 , t) be the transition matrix (pin,jk (t0 , t)) with pin,jk (t0 , t) = P (αε (t) = sjk |αε (t0 ) = sin ) for all sin , sjk ∈ M. ε ε Theorem 3.8. Assume (A1) and (A2) with n = 1. Then (3.21)
Pε (t0 , t) = Φ0 (t0 , t) + Ψ0 (t0 , (t − t0 )/ε) + O(ε + ε exp((t − t0 )/ε))
uniformly in (t0 , t), where 0 ≤ t0 ≤ t ≤ T , Φ0 (t0 , t) = 11Θ(t0 , t)ν(t),
dΨ0 (t0 , τ ) = Ψ0 (t0 , τ )A(t0 ), Ψ0 (t0 , t0 ) = I − Φ0 (t0 , t0 ), dτ
(3.22) where Θ(t0 , t) = (θij (t0 , t)) ∈ Rl×l is the solution of
d Θ(t0 , t) = Θ(t0 , t)B(t), Θ(t0 , t0 ) = I. dt
(3.23)
Remark 3.9. Owing to Theorems 3.7 and 3.8, for some κ0 > 0, pε (t) = ν(t) + O(ε + exp(−κ0 t/ε)), Pε (t0 , t) = Φ0 (t0 , t) + O(ε + exp(−κ0 (t − t0 )/ε)). 4. Occupation measures, aggregation, and switching diffusion. This section presents further asymptotic results that are of probabilistic feature. The statements of the results are given, and the proofs are relegated to the appendix. 4.1. Occupations measures. For any i = 1, . . . , l and j = 1, 2, . . ., define sequences of occupation measures as t ij (4.1) µε (t) = z ij (s, αε (s))ds, z ij (t, αε (t)) = I{αε (t)=sij } − ν ij (t)I{αε (t)∈Mi } , 0
where ν ij (t) is the jth component of the quasi-stationary distribution ν i (t) as defined in (A2). As a preparation, we first derive an order of magnitude estimate for µij ε (·) defined in (4.1). 2 Theorem 4.1. Under the conditions of Theorem 3.8, supt∈[0,T ] E(µij ε (t)) = O(ε). As alluded to in the introduction, we wish to reduce the complexity of the queueing network by aggregating states in each subspace as a single state. This leads to the definition of αε (t) = i if αε (t) ∈ Mi . Theorem 4.2. Under the conditions of Theorem 3.8, αε (·) converges weakly to α(·), a Markov chain generated by B(t) given in (3.6). 4.2. Switching diffusion limit. Let f (·) be a real-valued function defined on M satisfying that {f (sij ) : 1 ≤ i ≤ l, 1 ≤ j < ∞} ∈ 1 . For t ∈ [0, T ], define a sequence of real-valued functions by (4.2)
xε (t) =
t l ∞ 1 √ f (sij )[I{αε (u)=sij } − ν ij (u)I{αε (u)∈Mi } ]du. ε 0 i=1 j=1
Our interest lies in the asymptotic properties of xε (·). To obtain the desired limit property, we consider a pair of processes {Yε (·)} = {xε (·), αε (·)} and aim to show that Yε (·) converges weakly to Y (·) with a suitable generator.
576
G. YIN AND HANQIN ZHANG
Lemma 4.3. Under the conditions of Theorem 3.8, {Yε (·)} is tight in D2 [0, T ] (the space of R2 -valued functions that are right continuous and have left limits endowed with the Skorohod topology). Since {Yε (·)} is tight, by Prohorov’s theorem, we can extract a weakly convergent subsequence. Select such a subsequence and still use ε as its index for notational simplicity. Denote the limit by Y (·) = (x(·), α(·)). We proceed to characterize this limit process. Let CL2 be the collection of functions having bounded derivatives up to the second order with the second derivative being Lipschitz continuous. For each i = 1, . . . , l, g(·, i) ∈ CL2 , define an operator D(t) by (4.3)
D(t)g(x, i) =
∂2 1 2 σ (t, i) 2 g(x, i) + B(t)g(x, ·)(i), 2 ∂x
where σ 2 (t, i) > 0 is a smooth function. Theorem 4.4. Under the conditions of Theorem 3.8, Yε (·) converges weakly to Y (·), which is a solution of the martingale problem with operator D(·) defined by (4.3). Using the above theorem, we can also get a weak convergence of a reflected process. Define rε (t) = xε (t)−inf 0≤u≤t {xε (u)}. The weak convergence of xε (t) and the continuous mapping theorem (see [2, p. 30, Theorem 5.1]) yield that rε (·) converges weakly to r(·), a reflected switching diffusion process. 5. Extensions and examples. In this section we give the extensions of the results and present examples of queueing problems. 5.1. Infinitely many blocks. We consider the case (5.1)
Qε (t) =
A(t) + B(t) ε
with A(t) = diag(A1 (t), A2 (t), . . .).
That is, the matrix A(t) given by (2.2) is a diagonal block one with infinitely many blocks. Instead of condition (A2), we assume (A2 ). l (A2 ) There is a κl > 0 such that for exp(A (0)τ ) − 11ν l (0)∞ ≤ K exp(−κl τ ), any real number τ > 0, where ν l (t) = (ν l1 (t), ν l2 (t), . . .) is the quasi-stationary distribution corresponding to the generator Al (t), and inf l≥1 κl > 0. Theorem 5.1. Under (A1) and (A2 ), Theorem 3.7, Theorem 4.2, and Theorem 4.4 hold. The proof of this theorem follows the same line of argument as that of Theorems 3.7, 4.2, and 4.4 and is thus omitted. A special case of the theorem is that Qε (t) has the form (5.1) with A(t) having infinitely many blocks such that each Ai (t) is a generator of a finite-state Markov chain. Example 5.2. Consider a two-station queueing system, where each station has a single server and an exogenous Poisson arrival process. For the first station, having completed service there, customers either proceed to a queue in front of the second station with probability p1 or depart the system with probability (1 − p1 ). For the second station, after completing service there, they depart the system with probability (1 − p2 ) and go to a queue in front of the first station with probability p2 . Service is rendered in the order of the arrivals at each station. We assume that the first station can hold at most a total of m0 customers (including the customer in service) and any further arriving customers from outside will in fact be refused entry and hence will depart immediately without service at the first station, while the second station
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
577
has an unlimited waiting room. Furthermore, we assume that the system is a timedependent Markovian queueing network. Specifically, the rate of station i’s customer arriving from outside is λi (t), and the service rate is µi (t) (i = 1, 2). Assume that λi (t) and µi (t) are smooth, positive, real analytic functions of time t. We let Q(t) = (Q2 (t), Q1 (t)) be the Markovian queue length process that Qi (t) equals the number of customers at station i at time t. Then the state space of the queueing system can be represented by {(k, i) : 0 ≤ k < ∞ and 0 ≤ i ≤ m0 }. To determine the generator Q(t) = (q (j,n),(k,i) (t)), let π1 (t) = (1 − p1 )µ1 (t), π2 (t) = (1 − p2 )µ2 (t), π(t) = λ1 (t) + π1 (t), π ˜ (t) = λ2 (t) + p1 µ1 (t) + µ2 (t), ⎞ ⎛ −λ (t) λ (t) 1 1 ⎟ ⎜ π1 (t) −π(t) λ1 (t) ⎟ ⎜ 1 .. .. .. ⎟, A (t) = ⎜ . . . ⎟ ⎜ ⎠ ⎝ π1 (t) −π(t) λ1 (t) π1 (t) −π1 (t) B 1 (t) = diag(−λ2 (t), −[λ2 (t) + p1 µ1 (t)], . . . , −[λ2 (t) + p1 µ1 (t)]), B 2 (t) = diag(−(λ2 (t) + µ2 (t)), −˜ π (t), . . . , −˜ π (t)), ⎞ ⎛ λ (t) 2 ⎟ ⎜ p1 µ1 (t) λ2 (t) ⎟, C(t) = ⎜ .. .. ⎠ ⎝ . . ⎛ π (t) p µ (t) 2 2 2 .. ⎜ . D(t) = ⎜ ⎝
p1 µ1 (t) λ2 (t) ..
. π2 (t) p2 µ2 (t) µ2 (t)
⎞ ⎟ ⎟, ⎠
where each matrix is a (m0 + 1) × (m0 + 1) matrix. Then some tedious calculations yield ⎛ 1 ⎞ B (t) C(t) 2 ⎜ ⎟ (5.2) Q(t) = diag(A1 (t), A1 (t), . . .) + ⎝ D(t) B (t) C(t) ⎠. .. .. .. . . . For an initial time point t0 with t0 ∈ [0, T ], let P (t0 , t) with t > t0 be the transition matrix P (t0 , t) = (p(j,n),(k,i) (t0 , t)) with p(j,n),(k,i) (t0 , t) = P (Q(t) = (k, i)|Q(t0 ) = (j, n)) for all 0 ≤ j, k < ∞, 0 ≤ n, i ≤ m0 . Then we have the following system of equations for a quasi-birth-death process: d P (t0 , t) = P (t0 , t)Q(t). dt
(5.3) Assume that (5.4)
(1 − p1 ) p1 , λ1 (t) p1 µ1 (t), λ1 (t) λ2 (t), λ1 (t) µ2 (t). 1 inf t0 ≤t≤T λ1 (t) .
Then ν i (t) = (ν i0 (t), ν i1 (t), . . . , ν im0 (t)) with ⎤ ⎡ j j m0 λ1 (t) ⎦ λ1 (t) . ν 1j (t) = p0 (t), p0 (t) = 1 ⎣1 + π1 (t) π1 (t) j=1
Introduce ε =
578
G. YIN AND HANQIN ZHANG
Let π (t) = [λ2 (t) + p1 µ1 (t)] − p1 µ1 (t)p0 (t). It follows from some tedious calculations that ⎞ ⎛ − π (t) π (t) ⎟ ⎜ π (t) − µ2 (t) π (t) (5.5) B(t) = ⎝ µ2 (t) − ⎠. .. .. .. . . . By Theorem 5.1,
(5.6)
t − t0 Pε (t0 , t) = Φ0 (t0 , t) + [I − Φ(t0 , t0 )] exp A(t0 ) ε t − t0 + O ε + exp ε
uniformly in (t0 , t), where for 0 ≤ t0 ≤ t ≤ T , Φ0 (t0 , t) is given by (3.22) with Θ(t0 , t) satisfying (3.23). Using [23, Theorem 5.2, p. 128], the solution of (3.23) can be explicitly written as Θ(t0 , t) = U (t0 , t), where U (t0 , t) is a solution operator. Consequently Φ(t0 , t) given in (5.6) takes the form Φ(t0 , t) = 11U (t0 , t)ν. In particular, when the generators are time independent, U (t0 , t) becomes a semigroup [23, Chapter 1] and may be expressed as U (t0 , t) = exp((t − t0 )B). Next we consider the queue length process at the first station, Q1 (t). In general, Q1 (t) is not a Markov process. But from Theorem 5.1, Q1 (t) can be approximated very well by a Markov process α(t) with generator B(t) defined by (5.5). Furthermore, we consider a refined approximation of Q(t) = (Q2 (t), Q1 (t)). Define ∞ m0 t I{Q(u)=(i,j)} − ν (i,j) (u)I{Q1 (u)=j} du. It follows from Theoε (t) = i=0 j=0 0 √ rem 5.1 that (1/ ε)ε (t) can be well approximated by the switching diffusion prot cess 0 σ(u, α(u))dw(u), where w(·) is a standard one-dimensional Brownian motion, ∞ lk m0 m0 0k ∞ kl 0l σ 2 (t, i) = k=0 l=0 ν (t) 0 ψ0 (i, t, u)du + ν (t) 0 ψ0 (i, t, u)du for 0 ≤ i < ∞, and ψ0kl (i, s, t) is the (k, l)th entry of Ψ0 (i, t, u) given by ⎛ ⎛ 1 ⎞⎞ ν (t) · 1 ⎝ ⎝ Ψ0 (i, t, u) = I − : ⎠⎠ exp(A (t)u). 1 ν (t) Remark 5.3. References [12] and [18] establish asymptotic expansions for the queue length distribution of the time inhomogeneous single serve queue, a timedependent pure birth-death process. Here we consider a quasi-birth-death process. Furthermore, the queue length process at the first station (generally non-Markovian) can be approximated well by a Markov process with generator B(t). In [20], uniform acceleration expansions for time-varying generators were treated, and in [28] diffusion approximation was also considered. In these references, the Markov chains have one ergodic class, whereas in the current paper, multiergodic classes are considered. Compared with the diffusion approximations in [19] and [27], the switching diffusion approximation given here is related to the queue length process on the interval [0, t] and provides the evolution of the scaled sequence. The usual diffusion approximation leads to asymptotic normality, whereas switching diffusion limit yields Gaussian mixture distribution. The [λ1 (t) + p2 λ2 (t)]/[(1 − p1 p2 )µ1 (t)] and [λ2 (t) + p1 µ1 (t)]/[(1 − p1 p2 )µ2 (t)] are called the traffic intensities at station 1 and station 2, respectively. If [λ1 (t) + p2 λ2 (t)]/[(1 − p1 p2 )µ1 (t)] and [λ2 (t) + p1 µ1 (t)]/[(1 − p1 p2 )µ2 (t)] are less than but close
579
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
to one, the system is regarded as in heavy traffic, whereas if [λ1 (t) + p2 λ2 (t)]/[(1 − p1 p2 )µ1 (t)] and [λ2 (t) + p1 µ1 (t)]/[(1 − p1 p2 )µ2 (t)] are much less than one, the system is considered to be in light traffic. Our approximations are valid for either the heavy traffic case or the light traffic case as long as (5.4) holds. Using the Laplace transform technique, from the balance equation of the system, one may proceed as in [15] to carry out heavy traffic analysis for the sojourn time of time-homogeneous Markovian tandem queues with two servers. Comparing with (5.2), we can consider a general quasi-birth-death process with state space {(i, j), i ≥ 0, 1 ≤ j ≤ m0 } with a generator G(t) given by ⎛
⎞
A1 (t) − diag(C 1 (t)11) C 1 (t) 1 2 ⎜ B (t) A (t) − diag(B 1 (t)11 + C 2 (t)11) C 2 (t) (5.7) ⎝ .. .. . .
⎟ ⎠,
where i and j are called level and phase, respectively; see [22]. If the transitions between levels are much less frequent than the transitions between the phases inside the same level, then Q(t) given by (5.7) can be written as diag(A1⎛ (t), A2 (t), . . .) −diag(C 1 (t)11) C 1 (t) 1 1 ⎜ B (t) −diag(B (t)11 + C 2 (t)11) C 2 (t) +ε⎝ .. .. . .
⎞ ..
⎟ ⎠. .
Hence following Example 5.2, we can also get the asymptotic probability distribution of the process. 5.2. Asymptotic properties under γ-norm. For v = (v 1 , . . . , v l ) with v i = (v , v i2 , . . .), the γ-norm (which was named ν-norm in [21]; see also [1]) is defined as ||v||γ = max1≤i≤l sup1≤k 0, where ν l (t) = (ν l1 (t), ν l2 (t), . . .) is the quasi-stationary distribution corresponding to the generator Al (t), and inf l≥1 κl > 0. i1
Similar to Theorem 3.7, using essentially the same techniques, we obtain the following result. Theorem 5.4. Under (A1) and (A2 ), we construct two sequences {φk (t)}nk=0 and {ψk (t/ε)}nk=0 by Theorem 3.3 such that φk (t) ∈ C n+1−k and ψk (t/ε) decay exponentially fast. Moreover, with eε,n (t) defined by (3.20), supt∈[0,T ] ||eε,n (t)||γ = O(εn+1 ). Example 5.5. Consider the queueing system given by Example 5.2, but allow the first station to have unlimited waiting rooms. Then the state space of (Q2 (t), Q1 (t))
580
G. YIN AND HANQIN ZHANG
can be represented by {(k, i) : 0 ≤ k < ∞ and 0 ≤ i < ∞}. Let ⎞ ⎛ −λ1 (t) λ1 (t) ⎟ 1 (t) = ⎜ A ⎠, ⎝ (1 − p1 )µ1 (t) −π(t) λ1 (t) .. .. .. . . . 1 (t) = diag(−λ2 (t), −[λ2 (t) + p1 µ1 (t)], . . .), B 2 (t) = diag(−(λ2 (t) + µ2 (t)), −˜ π (t), B ⎞ . . .), ⎛ ⎛ π2 (t) p2 µ2 (t) λ2 (t) ⎟ ⎜ (t) λ (t) p π2 (t) p2 µ2 (t) µ =⎜ 1 2 1 , D(t) = C(t) ⎠ ⎝ ⎝ .. .. .. . . .
⎞ ..
⎟ ⎠. .
Then, some tedious calculations yield ⎛
(5.8)
1 (t) C(t) B ⎜ 1 1 (t), . . .) + ⎝ D(t) B 2 (t) C(t) = diag(A (t), A Q(t) .. .. . .
⎞ ..
⎟ ⎠. .
Letting γ 1i (0) = ((1−p1 )µ1 (0)/λ1 (0))i , then we have (A2 ) holds; see [21]. Therefore, by Theorem 5.4, we obtain the asymptotic probability distribution of (Q2 (t), Q1 (t)) under γ-norm. Appendix A. Proof of Theorem 4.1. Suppose that 0 ≤ s ≤ t. For each i = 1, . . . , l and j ∈ N, define ζε1,ij (s, t) = P (αε (s) = sij , αε (t) = sij ) − ν ij (t)P (αε (s) = sij , αε (t) ∈ Mi ), (A.1) ζε2,ij (s, t) = ν ij (s)ν ij (t)P (αε (s) ∈ Mi , αε (t) ∈ Mi ) −ν ij (s)P (αε (s) ∈ Mi , αε (t) = sij ). Theorem 3.8 and Remark 3.9 yield ζε1,ij (s, t) = O(ε + exp(−κ0 (t − s)/ε)). Similarly, ij 2 ζε2,ij (s, t) = O(ε + exp(−κ0 (t − s)/ε)). It follows from µij ε (0) = 0 that E(µε (t)) = t s 1,ij 2,ij 2 0 0 (ζε (r, s) + ζε (r, s))drds. The desired order estimates then follows. Proof of Theorem 4.2. The theorem will be proved in two steps. In the first step, we establish the tightness of {αε (·)}, and in the second step, we characterize the limit process. l l (1) Tightness. Note that αε (t) = i=1 iI{αε (t)=i} = i=1 iI{αε (t)∈Mi } . Define χε (t) = (χ1ε (t), . . . , χlε (t)), χiε (t) = (I{αε (t)=sij } ), χε (t) = (I{αε (t)=1} , . . . , I{αε (t)=l} ). In view of the definition of αε (t) and the Cram´er–Wold theorem [2, p. 49], to prove the tightness of {αε (·)}, it suffices to derive that of {χε (·)}. For any δ > 0, any t > 0, and any s > 0 with δ > s > 0 and t + s ∈ [0, T ], owing to the Markov property, E(χε (t + t+s s) − χε (t) − t χε (u)Qε (u)du|Ft,ε ) = 0, where Ft,ε denotes the σ-algebra generated by {αε (u) : u ≤ t}. Postmultiplying by 11 and noting A(t)11 = 0 lead to E(χε (t + s) − t+s χε (t) − t χε (u)B(u)11du|Ft,ε ) = 0. Thus by (A1), E(χε (t + s)|Ft,ε ) − χε (t) = O(s). t+s This implies that E(|χε (t+s)−χε (t)|2 Ft,ε ) = E(| t χε (u)B(u)11du|2 Ft,ε ) = O(s). Consequently, for 0 < δ ≤ s, limδ→0 lim supε→0 E|χε (t + s) − χε (t)|2 = 0. The desired tightness then follows from [16, p. 47]. (2) Characterization of the limit. Since {αε (·)} is tight, by Prohorov’s theorem, we can extract a weakly convergent subsequence. Select such a sequence and still use ε as its index for notational simplicity; denote its limit by α(·). We proceed to
581
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
characterize the limit process. We need to show that for any bounded function g(·) defined on {1, . . . , l},
t
g(α(t)) −
(A.2)
B(u)g(·)(α(u))du is a martingale. 0
To establish (A.2), it suffices that for any positive integers m, and k ≤ m, any bounded functions hk (·) defined on {1, . . . , l}, and any 0 < tk ≤ t ≤ t + s, (A.3) E
m
hk (α(tk )) g(α(t + s)) − g(α(t)) −
t+s
B(u)g(·)(α(u))du
= 0.
t
k=1
l To proceed, define g(αε (s)) = i=1 g(i)I{αε (s)∈Mi } . Then, g(αε (s)) = g(αε (s)). t Thus, g(αε (t)) − 0 Qε (u)g(·)(αε (u))du is a martingale. Consequently, (A.4)
m
t+s
Ehk (αε (tk ))[g(αε (t + s)) − g(αε (t)) −
Qε (u)g(·)(αε (u))du)] = 0. t
k=1
By virtue of the weak convergence and the Skorohod representation, as ε → 0, (A.5) m m E hk (α(tk ))[g(α(t + s)) − g(α(t))]. hk (αε (tk ))[g(αε (t + s) − g(αε (t))] → E k=1
k=1
The definition of g(·) leads to A(u)g(·)(αε (u)) = 0 for all u ∈ [0, T ]. As a result, (A.6) m E hk (αε (tk )) k=1 l ∞
=
+
t+s
Qε (u)g(·)(αε (u))du
t
E
i=1 j=1 l ∞ i=1 j=1
m
hk (αε (tk ))[
k=1 m
E
t+s
ν ij (u)I{αε (u)=i} B(u)g(·)(sij )du]
t
hk (αε (tk ))[
t+s
t
k=1
[I{αε (u)=sij } − ν ij (u)I{αε (u)=i} ]B(u)g(·)(sij )du].
Using Theorem 4.1, the last term in (A.6) goes to 0 as ε → 0. By virtue of the weak convergence of αε (·) and the Skorohod representation, (A.6) yields that E (A.7)
m
hk (αε (tk ))
k=1
→E
Qε (u)g(·)(αε (u))du t
m k=1
t+s
hk (α(tk ))
t+s
B(u)g(·)(α(u))du .
t
Combining (A.4), (A.5), and (A.7), (A.3) is verified. This completes the proof. Proof of Lemma 4.3. Note that Yε (·) is a sequence of vector-valued random processes with two components. Again, using the Cr´amer–Wold theorem, since {αε (·)} is tight, to prove the tightness of {Yε (·)}, it suffices to verify the tightness of {xε (·)}.
582
G. YIN AND HANQIN ZHANG
We claim that for any δ > 0, 0 ≤ s ≤ t, and t − s ≤ δ, sup E (|xε (t) − xε (s)|2 Fs,ε ) ≤ K(t − s). 0≤s≤t≤T
First, for u ∈ [s, t], E(z ij (u, αε (u))Fs,ε ) = O(ε + exp(−κ0 (t − s)/ε)) uniformly in i and j, by (4.2), the Markov property, and Theorem 3.8. Thus
(A.8)
t l ∞ 1 E(xε (t) − xε (s) Fs,ε ) = f (sij ) √ E(z ij (u, αε (u))Fs,ε )du ε s i=1 j=1 √ 1 = √ O(ε + exp(−κ0 (t − s)/ε)) = O( ε). ε
Denote ∆ε (s, t) = E((xε (t) − xε (s))2 Fs,ε ) ⎛⎛ ⎞ ⎞2 l ∞ t 1 ⎜ ⎟ = E ⎝⎝ f (sij ) zεij (αε (u))du⎠ Fs,ε ⎠ . ε s i=1 j=1 t l ∞ d It then follows dt ∆ε (s, t) = 2ε i=1 j=1 f (sij ) s E(ζε1,ij (u, t) + ζε2,ij (u, t)Fs,ε )du, where ζε1,ij (u, t) and ζε2,ij (u, t) are defined in (A.1). Thus (d/dt)∆ε (s, t) = O(1), ∆ε (s, s) = 0. An integration then leads to ∆ε (s, t) = O(t − s). Therefore, lim lim sup sup E(E(xε (t) − xε (s))2 Fs,ε ) = 0. δ→0
ε→0
0≤t−s≤δ
Moreover, for each δ > 0, and each rational t ≥ 0, using the Chebyshev’s inequality 2 2 2 inf ε P (|xε (t)| ≤ Kt,δ ) ≥ inf ε (1 − E|xε (t)| /(Kt,δ ) ) ≥ 1 − Kt/(Kt,δ ) . Select Kt,δ > KT /δ. Then inf ε P (|xε (t)| ≤ Kt,δ ) ≥ 1 − δ. This inequality, (A.8), and the criterion [16] then yield the desired tightness. Proof of Theorem 4.4. To prove the theorem, it suffices to show that t g(x(t), α(t)) − (A.9) D(u)g(x(u), α(u))du is a martingale. 0
To verify this martingale property, we prove that for any positive integer m, any k ≤ m, any bounded and continuous function hk (·, α) for each α, and any 0 < tk ≤ t ≤ t + s, m k=1
Ehk (x(tk ), α(tk ))
g(x(t + s), α(t + s)) − g(x(t), α(t)) t+s − D(u)g(x(u), α(u))du = 0. t
To accomplish this, we begin with the process indexed by ε. (1) Using an argument similar to [29, pp. 199–200], it can be verified that the martingale problem associated with operator D(t) has a unique (in the sense of distribution) solution. (2) In what follows, we often need to carry out estimates involving xε (t) and αε (t) intertwined. To untangle them, as a preparation, let η(t, x) be a real-valued function
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
583
that is Lipschitz continuous in (t, x) ∈ R2 . Using an argument similar to that of [29, Lemma 7.4, pp. 189–192], we can show t 2 (A.10) z ij (u, αε (u))η(u, xε (u))du = 0. lim sup E ε→0 0≤t≤T 0
(3) Define an operator Dε (t) by (A.11)
⎛ ⎞ l ∞ ∂ρ(t, x, α) ∂ρ(t, x, α) 1 Dε (t)ρ(t, x, α) = f (sij )z ij (t, α)⎠ +√ ⎝ ∂t ∂x ε i=1 j=1 + Qε (t)ρ(t, x, ·)(α) for α ∈ M, ρ(·, ·, α) ∈ C 1,1 .
t Then ρ(t, xε (t), αε (t)) − 0 Dε (u)ρ(u, xε (u), αε (u))du is a martingale; see [6, Chapter 2]. l Define g(x, α) = i=1 g(x, i)I{α∈Mi } . It is easily seen that A(u)g(x, ·)(αε (u)) = 0. Thus, (A.12) 1 ∂ Dε (u)g(xε (u), αε (u)) = √ b(u, αε (u)) g(xε (u), αε (u)) + B(u)g(xε (u), ·)(αε (u)), ∂x ε l ∞ where b(u, α) = i=1 j=1 f (sij )z ij (u, α). To obtain the desired limit, we use the techniques of perturbed test function methods [16], which require us to introduce a perturbation that is small in magnitude and that results in desired cancellations. To this end, define the perturbation by g(t, x, α) such that g(·, x, α) is uniform Lipschitz in t; both g(·) and (∂/∂x) g (·) are bounded; (∂/∂x) g (·, ·, α) is Lipschitz in (t, x) such that it is the solution of ∂g(x, α) A(t) g (t, x, ·)(α) = −b(t, α) . ∂x
(A.13)
It can be shown that such a function exists similar to [29, Remark 7.16]. Next, define the perturbed test function gε (·) by √ gε (t, x, α) = g(x, α) + ε (A.14) g (t, x, α). t Then using the definition of Dε (t), gε (t, xε (t), αε (t)) − 0 Dε (u)gε (u, xε (u), αε (u))du is a martingale. Consequently, m Ehk (xε (tk ), αε (tk )) gε (t + s, xε (t + s), αε (t + s)) − gε (t, xε (t), αε (t)) k=1 t+s
−
Dε (u)gε (u, xε (u), αε (u))du = 0. t
By the weak convergence of Yε (·) to Y (·), the Skorohod representation, and (A.14), m
lim
ε→0
Ehk (xε (tk ), αε (tk ))[gε (t + s, xε (t + s), αε (t + s)) − gε (t, xε (t), αε (t))]
k=1 m
=
k=1
Ehk (x(tk ), α(tk ))[g(x(t + s), α(t + s)) − g(x(t), α(t))].
584
G. YIN AND HANQIN ZHANG
Note that from (A.11), √
εDε (u) g (u, xε (u), αε (u)) =
l ∞
f (sij )z ij (u, αε (u))
i=1 j=1
(A.15)
∂ g (u, xε (u), αε (u)) ∂x
√ ∂ g (u, xε (u), αε (u)) + ε ∂u √ 1 + εB(u) g (u, xε (u), ·)(αε (u)). g (u, xε (u), ·)(αε (u)) + √ A(u) ε
It follows from (A.12), (A.13), and (A.15), upon cancellation, that (A.16) l ∞
Dε (u)gε (u, xε (u), αε (u)) =
f (sij )z ij (u, αε (u))
i=1 j=1
∂ g (u, xε (u), αε (u)) ∂x
+ B(u)g(xε (u), ·)(αε (u)) √ ∂ g (u, xε (u), αε (u)) √ + ε g (u, xε (u), ·)(αε (u)). + εB(u) ∂u Thus again, by the weak convergence and the Skorohod representation, t+s m lim Dε (u)gε (u, xε (u), αε (u))du Ehk (xε (tk ), αε (tk )) ε→0 t k=1 m l ∞ t+s ∂ g (u, xε (u), αε (u)) = lim Ehk (xε (tk ), αε (tk )) f (sij )z ij (u, αε (u)) du ε→0 ∂x t i=1 j=1 k=1 t+s B(u)g(xε (u), ·)(αε (u))du . + t
By virtue of the weak convergence of Yε (·) to Y (·), the Skorohod representation, t and the boundedness of g(·), (∂/∂x) g (·), we have 0 B(u)g(xε (u), ·)(αε (u))du → t B(u)g(x(u), ·)(α(u))du. 0 Define ∂ g (t, x, α) b(t, x, α) = b(t, α) . ∂x
(A.17) Then
t
b(u, xε (u), αε (u))du = (A.18)
0
+
t l ∞
t l ∞
b(u, xε (u), sij )ν ij (u)I{αε (u)=i} du
0 i=1 j=1
[I{αε (u)=sij } − ν ij (u)I{αε (u)=i} ]b(u, xε (u), sij )du.
0 i=1 j=1
t By (A.10), the last term in (A.18) goes to 0 in mean squares, so 0 b(u, xε (u), αε (u))du t ∞ converges to 0 b(u, x(u), α(u))du, where b(u, x, i) = j=1 ν ij (u)b(u, x, sij ). Combining the estimates obtained so far, we arrive at t+s Ehk (xε (tk ), αε (tk )) Dε (t)gε (u, xε (u), αε (u))du t t+s t+s b(u, x(u), α(u))du + B(u)g(x(u), ·)(α(u))du . → Ehk (x(tk ), α(tk )) t
t
585
TWO-TIME-SCALE MARKOV CHAINS AND APPLICATIONS
(4) In view of (A.13), for i = 1, . . . , l, ⎞ ⎞ ⎛ b(t, si1 ) ∂g(x,i) g(t, x, si1 ) ∂x ⎟ ⎜ Ai (t) ⎝ g(t, x, si2 ) ⎠ = ⎝ b(t, si2 ) ∂g(x,i) ⎠, ∂x · · : : ⎛
which has a unique solution. This implies that g(t, x, sij ) is a function of (∂/∂x)g(x, i). g (t, x, sij ). Thus b(t, x, sij ) is a In view of (A.17), b(t, x, sij ) is a function of (∂/∂x) function of (∂ 2 /∂x2 )g(x, i). Denote b(t, x, i) = (1/2)a(t, i)(∂ 2 /∂x2 )g(x, i), where a(t, i) is an appropriate function. Using an argument similar to that [29, pp. 200–203], it can be shown that a(t, i) ≥ 0. Thus, a(t, i) can be written as a(t, i) = σ 2 (t, i), where for each i = 1, 2, . . . , l, ∞ ∞ ∞ ik 2 (A.19) f (sik )f (sim ) ν (t) ψ0km (i, u, t)du σ (t, i) = 0
k=1 m=1
+ν
im
(t)
∞
ψ0mk (i, u, t)du
,
0
and ψ0km (i, s, t) is the (k, m)th entry of Ψ0 (i, s, t) given by Ψ0 (i, s, t) = (I − (ν i (s), ν i (s), . . .) ) exp(Ai (s)t). This specifies the covariance structure of the limit process and establishes the desired martingale property and hence the theorem follows. Acknowledgment. We thank the editor and the reviewers for detailed comments and suggestions on an early version of the manuscript, which have led to much improvement of the paper. REFERENCES [1] E. Altman, K. E. Avrachenkov, and R. Nunez-Queija, Pertrubation analysis for denumerable Markov chains with applications to queueing models, Adv. Appl. Probab., 36 (2004), pp. 839–853. [2] P. Billingsley, Convergence of Probability Measures, Wiley, New York, 1968. [3] G. Bitran and D. Tirupati, Multiproduct queueing networks with deterministic routing: Decomposition approach and the notion of interference, Management Sci., 34 (1988), pp. 75– 100. [4] N. N. Bogoliubov and Y. A. Mitropolskii, Asymptotic Methods in the Theory of Nonlinear Oscillator, Gordon and Breach, New York, 1961. [5] J. Daigle and J. Langford, Models for analysis of packet voice communication systems, IEEE J. Sel. Areas Commun., SAC-4 (1986), pp. 847-855. [6] M. H. A. Davis, Markov Models and Optimization, Chapman & Hall, London, 1993. [7] J. L. Doob, Stochastic Processes, Wiley Classic Library Edition, Wiley, New York, 1990. [8] S. G. Eick, W. A. Massey, and W. Whitt, The physics of the Mt /G/∞ queue, Oper. Res., 41 (1993), pp. 731–742. [9] V. Hutson and J. S. Pym, Applications of Functional Analysis and Operator Theory, Academic Press, London, 1980. [10] A. M. Il’in, A. S. Kalashnikov, and O. A. Oleinik, Linear equations of the second order of parabolic type, Russian Math. Surveys, 17 (1962), pp. 1–143. [11] A. M. Il’in, R. Z. Khasminskii, and G. Yin, Asymptotic expansions of solutions of integrodifferential equations for transition densities of singularly perturbed switching diffusions: Rapid switchings, J. Math. Anal. Appl., 238 (1999), pp. 516–539.
586
G. YIN AND HANQIN ZHANG
[12] J. Keller, Time-dependent queues, SIAM Rev., 24 (1982), pp. 401–412. [13] R. Z. Khasminskii, G. Yin, and Q. Zhang, Asymptotic expansions of singularly perturbed systems involving rapidly fluctuating Markov chains, SIAM J. Appl. Math., 56 (1996), pp. 277–293. [14] C. Knessl, B. J. Matkowsky, Z. Schuss, and C. Tier, Asymptotic analysis of a statedependent M/G/1 queueing system, SIAM J. Appl. Math., 46 (1986), pp. 483–505. [15] C. Knessl and J. A. Morrison, Heavy traffic analysis of the sojourn time in tandem queues with overtaking, SIAM J. Appl. Math., 51 (1991), pp. 1740–1763. [16] H. J. Kushner, Approximation and Weak Convergence Methods for Random Processes, with Applications to Stochastic Systems Theory, MIT Press, Cambridge, MA, 1984. [17] G. E. Ladas and V. Lakshmikantham, Differential Equations in Abstract Spaces, Academic Press, New York, 1972. [18] W. A. Massey, Asymptotic analysis of the time dependent M/M/1 queue, Math. Oper. Res., 10 (1985), pp. 305–327. [19] W. A. Massey and W. Whitt, Unstable asymptotics for nonstationary queues, Math. Oper. Res., 19 (1994), pp. 267–291. [20] W. A. Massey and W. Whitt, Uniform acceleration expansions for Markov chains with timevarying rates, Ann. Appl. Probab., 8 (1998), pp. 1130–1155. [21] S. P. Meyn and R. L. Tweedie, Computable bounds for geometric convergence rates of Markov chains, Ann. Appl. Probab., 4 (1994), pp. 981–1021. [22] M. F. Neuts, Matrix-Geometric Solutions in Stochastic Models, Johns Hopkins University Press, Baltimore, 1981. [23] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Springer-Verlag, New York, 1983. [24] M. I. Reiman, Asymptotically exact decomposition approximations for open queueing networks, Oper. Res. Lett., 9 (1990), pp. 363–370. [25] A. B. Vasil’eava and V. F. Butuzov, Asymptotic Methods in Singular Perturbations Theory, Vysshaya Shkola, Moscow, 1990 (in Russian). [26] W. Whitt, The queueing network analyzer, Bell Syst. Tech. J., 62 (1983), pp. 2779–2815. [27] W. Whitt, Stochastic-Process Limits, Springer-Verlag, New York, 2001. [28] G. Yin and H. Zhang, Countable-state-space Markov chains with two-time scales and applications to queueing systems, Adv. Appl. Probab., 34 (2002), pp. 662–688. [29] G. Yin and Q. Zhang, Continuous-Time Markov Chains and Applications: A Singular Perturbation Approach, Springer-Verlag, New York, 1998. [30] G. Yin, Q. Zhang, and G. Badowski, Asymptotic properties of a singularly perturbed Markov chain with inclusion of transient states, Ann. Appl. Probab., 10 (2000), pp. 549–572.
