- Email: [email protected]
Spatially-coherent uniformization of a stochastic fluid model to a Quasi-Birth-and-Death process
Performance Evaluation 70 (2013) 578–592
Contents lists available at ScienceDirect
Performance Evaluation journal homepage: www.elsevier.com/locate/peva
Spatially-coherent uniformization of a stochastic fluid model to a Quasi-Birth-and-Death process Nigel G. Bean a , Małgorzata M. O’Reilly b,∗ a
Applied Mathematics, University of Adelaide, SA 5005, Australia
b
School of Mathematics, University of Tasmania, Tas 7001, Australia
article
info
Article history: Available online 24 June 2013 Keywords: Stochastic fluid model Quasi-Birth-and-Death process Markov chain Uniformization
abstract We derive a uniformization of a stochastic fluid model (SFM) to a Quasi-Birth-and-Death process (QBD) that is spatially-coherent since the continuous level in the SFM has a natural correspondence to the discrete level in the QBD. As a consequence of this, the QBD can be used as a direct approximation of the original SFM, in those situations where a discrete state space is an advantage. We treat the unbounded as well as the bounded cases and illustrate the theory with a numerical example. The key fluid generator, Q, and matrix Ψ for the SFMs emerge from the QBD calculations in the natural limit. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Let {(ϕ(t )), t ≥ 0} be an irreducible continuous-time Markov Chain (CTMC) with a finite state space S = {1, 2, . . . , n} and infinitesimal generator T = [Ti,j ]. Let {(X (t ), ϕ(t )), t ≥ 0} be a Markovian stochastic fluid model (SFM) with phase variable ϕ(t ) and unbounded level variable X (t ) ∈ (−∞, +∞). There are real rates ci for all i ∈ S , such that when ϕ(t ) = i, then the rate at which the level is changing is ci . The buffer collecting fluid in this process is denoted by X . The CTMC is referred to as the driving process. The literature on SFMs has focused on the stationary and transient analysis, as well as the development of efficient algorithms [1–9] and performance-driven applications [10–14]. The first paper that used the class of effective methods known as matrix-analytic methods [15] to consider a Markovian SFM was the seminal paper by Ramaswami [9]. He proposed one way of mapping an SFM to a Quasi-Birth-and-Death process (QBD). However, that mapping is not spatially-coherent, in the sense that the level of the QBD does not correspond to the level of the SFM in any meaningful way. Instead, it was derived with the specific purpose of allowing exact algorithmic evaluation of the key matrix 9 for the SFMs, using known QBD techniques. Bean, O’Reilly and Taylor [6] developed a direct method for the analysis of SFMs which uses the key fluid generators Q(s) and does not require the construction of a QBD. In the work by Ahn and Ramaswami [2], Ramaswami’s mapping was further developed in such a way that the workload in the queue (not the level of the queue) represented the level of the SFM in an appropriate limit. In this paper we derive a spatially-coherent uniformization of an SFM to a QBD. This preserves a natural mapping between the continuous level in the SFM and the discrete level in the QBD, and so the QBD can be used as a direct approximation of the original SFM. We first construct a spatially-coherent discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} through a uniformization of the SFM. This construction yields a discrete-time QBD, but by incorporating the underlying Poisson process directly back into it, we develop the spatially-coherent continuous-time QBD {(X1x (t ), ϕ1x (t )), t ≥ 0}. A new class of important models, referred to as multi-dimensional SFMs [16–18,14], has recently emerged which are very hard to analyze. In particular, we are interested in the stochastic fluid–fluid model constructed in [18]. One attractive
∗
Corresponding author. E-mail addresses: [email protected] (N.G. Bean), [email protected] (M.M. O’Reilly).
0166-5316/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.peva.2013.05.006
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
579
approach to the analysis of these models is to discretize the first fluid and then treat the second fluid as being driven by a QBD. Therefore, since the QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} can be used as an approximation of the SFM, the significance of these results goes beyond the development of new ways of approximating the performance measures of an SFM. Instead it is possible to use this direct connection to translate known results, and methods of analysis, for QBDs to the field of SFMs. In particular, we expect to exploit this crucial feature in future analyses of multi-dimensional SFMs and thus provide a matrix-analytic methods based analysis of this important, but difficult to treat, class of models. We formalize the spatially-coherent uniformization in our paper and show that the transition function of the driving process behind the SFM is the same as that of the continuous-time QBD, following an argument similar to the usual uniformization argument. Hence, this QBD should be a good approximation for the SFM as it captures the stochastic phase-dynamics correctly; all we then need is to capture the level-dynamics up to the resolution allowed by the particular choice of 1x. By using a fairly natural choice of the uniformization parameter we are able to show that the level-dynamics converge in probability to those of the fluid. We then show that, in the limit as 1x → 0+ , the QBDs converge to the SFM (weak convergence of the processes), and that the key fluid generator Q and matrix 9 also emerge from the QBD calculations in this limit. We treat the unbounded as well as the bounded cases and illustrate the theory with a numerical example. The structure of the paper is as follows. First, we describe the uniformization method with general parameters ϑi (1x) > 0, i ∈ S . In Section 2 we construct the discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} and in Section 3 the continuous-time QBD {(X1x (t ), ϕ1x (t )), t ≥ 0}, for the unbounded SFM with X (t ) ∈ (−∞, +∞). Next, we focus on the appropriate choice of parameters ϑi (1x) in Section 4, and analyze the corresponding properties. The proofs of the main convergence results are given in Section 5. In Section 6 we show the connection with the fluid generator Q and matrix 9. The uniformization of the bounded SFM with X (t ) ∈ [0, B] is presented in Section 7. A numerical example is given in Section 8, followed by concluding remarks in Section 9.
¯ 1x (n)), n ≥ 0} 2. The spatially-coherent discrete-time QBD {(X¯ 1x (n), ϕ We partition the set of all phases as S = S1 ∪ S2 ∪ S0 , where S1 = {i ∈ S : ci > 0}, S2 = {i ∈ S : ci < 0}, S0 = {i ∈ S : ci = 0}, and the generator as T11 T21 T01
T=
T12 T22 T02
T10 T20 , T00
according to the partitioning of S . Further, let C1 = diag(ci )i∈S1 and C2 = diag(|ci |)i∈S2 . For a given positive 1x > 0, and for all i ∈ S , let
ϑi (1x) > 0
(1)
be some function of 1x. For the moment, we do not assume any particular form of ϑi (1x), except the positivity assumption. Later, in Section 4 we define a suitable form. We construct {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0}, the spatially-coherent discrete-time QBD with state space G¯ = {(k, i) : k ∈ Z, i ∈ S }, by connecting the phase and level of this QBD to the state of the SFM at arrival times of a Poisson process with an appropriately chosen parameter, following the standard uniformization approach [19]. Let
γi (1x) = λi + ϑi (1x),
(2)
where λi = Ti,i , for all i ∈ S . Choose
γ (1x) ≥ max {γi (1x)} ,
(3)
i
and consider the Poisson process {N (t ), t ≥ 0} with parameter γ (1x). X (0) The initial condition is (X¯ 1x (0), ϕ¯ 1x (0)) = (⌊ 1x ⌋, ϕ(0)). Then, by the standard theory of Markov Chains, for example
[19], assuming that an arrival in N (·) occurs at time t and that (X¯ 1x (N (t −)), ϕ¯ 1x (N (t −))) = (k, i),
• with probability p1 = λi /γ (1x) it is marked as a type-1 arrival and the process moves to some phase j ̸= i without a change in the discrete level so that a transition from (k, i) to (k, j) occurs (where j ̸= i is chosen with probability Ti,j /λi ), • with probability p2 = ϑi (1x)/γ (1x) it is marked as a type-2 arrival and the process does not change phase but changes level so that a transition from (k, i) to (k + 1, i) occurs if ci > 0, or a transition from (k, i) to (k − 1, i) occurs if ci < 0, or • with probability p3 = 1 − p1 − p2 it is marked as a type-3 arrival and the process remains in state (k, i). ¯ (1x) = [P¯ (1x)(k,i)(m,j) ] is given by Consequently, the corresponding one-step transition probability matrix P 1 − γ (1x)/γ (1x) m = k, j = i, c ̸= 0 i i 1 − λi /γ (1x) m = k, j = i, ci = 0 Ti,j /γ (1x) m = k, j ̸= i P¯ (1x)(k,i)(m,j) = m = k + 1, j = i, ci > 0 ϑi (1x)/γ (1x) m = k − 1, j = i, ci < 0 ϑi (1x)/γ (1x) 0
otherwise,
580
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
and can be written in the form
..
..
(−1) . ¯ ¯ (1x) = A (1x) P ¯ (0)
.
¯ (0) (1x) A ¯A(−1) (1x) ..
..
.
¯ (+1) (1x) A ¯ (0) (1x) A ..
.
¯ (0)
.
¯ (−1)
¯ (+1) (1x) A ..
.
,
¯ (−1)
where matrices A (1x) = [A (1x)i,j ]i,j∈S , A (1x) = [A (1x)i,j ]i,j∈S and A¯ (+1) (1x) = [A¯ (+1) (1x)i,j ]i,j∈S corresponding to the transitions to the same level, one level down and one level up, respectively, are given by
1 − γi (1x)/γ (1x) j = i, ci ̸= 0 ¯A(0) (1x)i,j = 1 − λi /γ (1x) j = i, ci = 0 Ti,j /γ (1x) j ̸= i ϑi (1x)/γ (1x) j = i, ci < 0 A¯ (−1) (1x)i,j = 0 otherwise ϑi (1x)/γ (1x) j = i, ci > 0 A¯ (+1) (1x)i,j = 0 otherwise. 3. The spatially-coherent continuous-time QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} In this section, using the discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} and standard techniques [19], we essentially reverse the uniformization to develop a continuous-time QBD. Specifically, we fix 1x > 0 and embed the events of the discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} at the time points of a Poisson process of rate γ (1x) to create the continuoustime Markov chain {(X1x (t ), ϕ1x (t )), t ≥ 0}. This CTMC then has state space G = {(k, i) : k ∈ Z, i ∈ S }, levels k ∈ Z, phases i ∈ S , and generator T(1x) = [T (1x)(k,i)(m,j) ] given by
T (1x)(k,i)(m,j) =
−γ (1x) i −λi
m = k, j = i, ci ̸= 0 m = k, j = i, ci = 0 m = k, j ̸= i m = k + 1, j = i, ci > 0 m = k − 1, j = i, ci < 0 otherwise.
Ti,j
ϑi (1x) ϑi (1x) 0
Clearly T(1x) is well defined. It has negative on-diagonals, nonnegative off-diagonals, and all rows sums equal to zero. We note that the above CTMC has the structure of an unbounded continuous-time QBD with
..
..
(−1) . A (1x) T(1x) = (0)
. (0) A (1x) A(−1) (1x) .. .
..
. (+1) A (1x) A(0) (1x) .. .
(0)
A(+1) (1x)
..
.
,
where matrices A (1x) = [A (1x)i,j ]i,j∈S , A (1x) = [A (1x)i,j ]i,j∈S and A(+1) (1x) = [A(+1) (1x)i,j ]i,j∈S corresponding to the transitions to the same level, one level down and one level up, respectively, are given by
−γi (1x) A (1x)i,j = −λi (0)
Ti,j A(−1) (1x)i,j =
A(+1) (1x)i,j =
ϑi (1x) 0
ϑi (1x) 0
(−1)
(−1)
j = i, ci ̸= 0 j = i, ci = 0 j ̸= i j = i, ci < 0 otherwise j = i, ci > 0 otherwise.
Now, assume (X1x (0), ϕ 1x (0)) = (k, i) for some k ∈ Z and i ∈ S . Let τi be the random variable that records the total time until a transition to some (m, j) with j ̸= i occurs. That is, τi is the time until the process {(X1x (t ), ϕ1x (t )), t ≥ 0} changes phase (regardless of level). Further, for all i, j ∈ S and t, define P (t )ij = P (ϕ(t ) = j|ϕ(0) = i), P1x (t )ij = P (ϕ1x (t ) = j|ϕ1x (0) = i). Below we establish a result which holds for any parameters ϑi (1x) > 0, i ∈ S , and guarantees the equivalence of the phase process in the uniformized QBD {(X1x (t ), ϕ1x (t )), t ≥ 0}, to the phase process in the original SFM, {(X (t ), ϕ(t )), t ≥ 0}.
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
581
Theorem 1. For any 1x > 0, the process {(X1x (t ), ϕ1x (t )), t ≥ 0} satisfies the following properties: (i) The random variable τi follows the exponential distribution with parameter λi . (ii) The probability that the process moves to a state with phase j, given that it leaves phase i (regardless of level), is Tij /λi . Therefore, for all i, j ∈ S and t, P (t )ij = P1x (t )ij .
(4)
Proof. We apply classic arguments from the theory of Markov Chains [19]. For i with ci = 0 the result follows trivially, since then, for any k ∈ Z, we have T(k,i)(k,i) = −λi and the only possible way for the process to leave phase i is to move from (k, i) to (k, j) for some j ̸= i. Suppose that ci ̸= 0. Without loss of generality assume ci > 0. First, note that the distribution of time spent in (m, i) is exponential with parameter γi (1x). An arrival of type 1, that is from (k, i) to (m, j) with j ̸= i can only occur if m = k. Also, note that for i with ci > 0 the only possible transitions during which a phase does not change are arrivals of type 2, that is from (m, i) to (m + 1, i) for some m. Now, given that the process makes n − 1 transitions before changing phase i for the first time at transition n, the distribution of τi is Erlang of order n and with parameter γi (1x). It therefore has density function given by
(γi (1x)t )n−1 γi (1x)e−γi (1x)t . (n − 1)! The probability that the process makes n − 1 transitions before changing phase i for the first time at transition n is equivalent to the probability that the process makes n − 1 transitions of type 1 (from (m, i) to (m + 1, i) for some m), and then one transition of type 2 (from (m, i) to (m, j) for some j ̸= i), which is equal to
λi 1− γi (1x)
n −1
λi . γi (1x)
Hence, by the Law of Total Probability, the density function of τi is given by ∞ (γi (1x)t )n−1 n =1
(n − 1)!
γi (1x)e−γi (1x)t 1 −
λi γi (1x)
n−1
λi γi (1x)
= λi e−γi (1x)t eγi (1x)t (1−λi /γi (1x)) = λi e−λi t , and so property (i) follows. Next, the probability that the process moves to a state with phase j ̸= i, given that a transition from (k, i) to a state with phase different than i occurs, is equal to Tij /γi (1x)
λi /γi (1x)
=
Tij
λi
,
which shows property (ii). Properties (i) and (ii) show that a CTMC obtained by observing only the phase in the QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} has the same probability law as the CTMC driving the original SFM {(X (t ), ϕ(t )), t ≥ 0}, before uniformization. Hence (4) follows. We are now interested in a suitable choice of parameters ϑi (1x) so as to establish a meaningful physical connection between the level of the QBD and the level of the SFM. Assume that the QBD and the SFM both start from level zero and in the same phase at time zero, and that ϕ(·) and ϕ1x (·) are coupled, in the sense that they share the same sample path. We desire the property such that at all times t ≥ 0, if the QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} is on level k, then the SFM {(X (t ), ϕ(t )), t ≥ 0} is near the level k1x with probability nearly 1. We formalize this idea in Section 4, followed by the proof of weak convergence of the processes in Section 5. 4. Level properties In this section, we assume for all i ∈ S that the parameters ϑi (1x) are of the form
ϑi (1x) =
|ci | , 1x
∀1x > 0,
(5)
and study the level-based properties of the process {(X1x (t ), ϕ1x (t )), t ≥ 0}. Consider the Poisson process {(Ni,1x (t )), t ≥ 0} with rate ϑi (1x) for every i ∈ S . Let Li,1x (t ) be the random variable defined by Li,1x (t ) = Ni,1x (t )1x. The following results are immediate from the properties of the underlying Poisson process.
582
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Lemma 1. For all 1x > 0 and t ≥ 0, E (Li,1x (t )) = |ci |t , and, for all T ≥ 0, lim Var(Li,1x (t )) = 0,
1x→0+
where the convergence is uniform on the interval t ∈ [0, T ]. Let A(t ) be the event that a change in phase does not occur during the time interval [0, t ] (but any number of changes in level may occur). Corollary 1. For all 1x > 0 and t ≥ 0, E (X1x (t )1x|(X1x (0), ϕ1x (0)) = (0, i), A(t )) = ci t , and, for all t ≥ 0, lim Var(X1x (t )1x|(X1x (t ), ϕ1x (0)) = (0, i), A(t )) = 0.
1x→0+
Theorem 2. For any ϵ > 0, i ∈ S , and T > 0, the family {Fi,1x (t )} of probabilities Fi,1x (t ) = P (Li,1x (t ) − |ci |t ≥ ϵ),
converges uniformly to 0, as 1x → 0+ , on the interval t ∈ [0, T ]. Proof. The result follows immediately from Lemma 1, since Var(Li,1x (t )) = E ((Li,1x (t ) − |ci |t )2 ), upon use of Markov’s inequality. 5. Weak convergence For all n = 1, 2, . . . , define Yn (·) = X1x (·)1x
and
ϕn (·) = ϕ1x (·) where 1x =
1 n
.
(6)
Assume without loss of generality that Yn (0) = 0 and X (0) = 0. Also, assume that parameters ϑi (1x) = ϑi (1/n), i ∈ S , are of the form (5). Consider the phase process {ϕ(·)} and {ϕn (·), n ≥ 1}. Let Ft be the σ -algebra for the process ϕ(·) up to and including time t and Ftn be the σ -algebra for the process ϕn (·) up to and including time t, for all n ≥ 1. Below, we establish the weak convergence of {Yn (·), ϕn (·)} ⇒ {X (·), ϕ(·)} as n → ∞. We need to work with two spaces R and S , and with the two associated Euclidean metrics r and s, making (R, r ) and (S , s) both metric spaces. In fact, they are both separable and complete metric spaces and, therefore, so is the product space (R, r ) × (S , s). Finally, both the SFM and the QBD have sample paths in DR×S [0, ∞) and we equip this space with the usual Skorokhod J1 -topology. For the technical details see Billingsley [20, Chapter 3] and Ethier and Kurtz [21, Chapter 3]. The structure of the argument in this section is as follows. Theorem 4 (which relies on Theorems 1 and 3) proves that for any individual t ≥ 0 that the distribution {Yn (t ), ϕn (t )} converges weakly to the distribution {Y (t ), ϕ(t )}; Theorem 5 proves that every sub-sequence of the sequence of processes {Yn (·), ϕn (·)} contains a further sub-sequence that is weakly convergent; finally, Theorem 6 shows that the finite-dimensional distributions converge weakly, which, in essence, shows that the process limit that was shown to exist in Theorem 5 must be unique since the sequence of processes is relative compact. This is sufficient to prove that the sequence of processes {Yn (·), ϕn (·)} weakly converges to the process {Y (·), ϕ(·)}. The theorems cited in the proofs are given in the Appendix, for completeness. Theorem 1 shows that ϕn (·) obeys the same probability law as ϕ(·). We shall therefore couple the processes by assuming that ϕn (·) and ϕ(·) share the same sample path, and so ϕn (t ) = ϕ(t ) for all t ≥ 0. Hence Ft ≡ Ftn . P
Theorem 3. {Yn (T )|FT } ⇒ X (T )|FT as n → ∞ and for all T ≥ 0. Proof. Given FT , define ui as the total time spent in phase i during the time interval [0, T ], for all i ∈ S . Therefore,
ui = T
i∈S
and X (T ) =
i∈S
Further, given FT , it is also clear that Yn (T ) =
i∈S
sign(ci )Li,n (ui ).
ci ui .
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
583
In Theorem 2 we showed that P (Li,n (t ) − |ci |t ≥ ϵ) → 0 uniformly in t ∈ [0, T ]. Therefore, since there are only finitely many phases |S |, ϵ P Li,n (ui ) − |ci |ui ≥ P Yn (T ) − X (T ) ≥ ϵ|FT ≤ |S | i∈S which converges to 0. Hence we have shown convergence in probability for the conditional probabilities.
D
Theorem 4. {(Yn (t ), ϕn (t ))} ⇒ (X (t ), ϕ(t )) as n → ∞ and for all t ≥ 0. Proof. For any measurable set A ⊂ R and j ∈ S , Theorem 1 shows that P (Yn (t ) ∈ A, ϕn (t ) = j) = P (Yn (t ) ∈ A, ϕ(t ) = j). Therefore, lim P (Yn (t ) ∈ A, ϕn (t ) = j) = lim
n→∞
n→∞
=
P (Yn (t ) ∈ A, ϕ(t ) = j|Ft )dP (Ft )
P (X (t ) ∈ A, ϕ(t ) = j|Ft )dP (Ft ),
by Theorem 3 since we have shown convergence in probability. The result then follows immediately.
Theorem 5. The sequence of processes {Yn (·), ϕn (·)} is relatively compact in DR×S [0, ∞). Proof. To prove relative compactness of {Yn (·), ϕn (·)} it suffices to prove relative compactness for both of the sequences {Yn (·)} and {ϕn (·)} separately [21, Theorem 3.2.2, p. 104 and Proposition 3.2.4, p. 107]. Now ϕn (t ) ∈ {1, 2, . . . , m}, for all t and so lives in a compact space. Prohorov’s Theorem [21, Theorem 3.2.2, Page 104] immediately implies that {ϕn (·)} is relatively compact. We are thus left to show that {Yn (·)} is relatively compact. We use Theorem 3.8.6 [21, p. 137] which first requires us to show that condition (a) of Theorem 3.7.2 [21, p. 128] holds. Condition (a) of Theorem 3.7.2 is that for every η > 0 and rational t ≥ 0, there exists a compact set Γη,t ∈ R such that η
inf P {Yn (t ) ∈ Γη,t } ≥ 1 − η,
(7)
n ≥1
η
where Γη,t = x ∈ S : infy∈Γη,t d(x, y) < η . Note that Eq. (7) is essentially equivalent to asserting tightness of the random variable Yn (t ), but only for rational t. To show that this holds, we shall define Γη,t = [−Z , Z ], Z > 0, and then prove the existence of a suitable Z . If Z is finite, then Γη,t is a compact set. Now, in order for Yn (t ) to escape the interval [−Z , Z ] there has to be at least Zn changes of level in the QBD. Let k = argmaxi∈S {|ci |} and so phase k represents the phase with the highest rate of level changes. Therefore,
P (Yn (t ) ̸∈ [−Z , Z ]) ≤ P (QBD has at least Zn changes of level in time t )
≤ P (Nk,n (t ) ≥ Zn). Recall that Nk,n (·) is the underlying Poisson process corresponding to the fastest rate of level changes in the QBD, and any move into a different phase will only lead to a slower Poisson process or possibly a phase that moves in the opposite direction. Now, we show that a suitable Z exists for all t and all η > 0. The Chernoff bound for the Poisson distribution implies that P (Nk,n (t ) ≥ Zn) ≤
e−β (eβ)Zn
(Zn)Zn
,
where β = ϑk (1/n)t = |ck |nt. We therefore need to choose Z such that fn (Z ) =
e−|ck |nt (e|ck |nt )Zn
(Zn)Zn
≤ η,
∀n ≥ 1.
If Z ≥ e|ck |t, then fn (Z ) is decreasing in n and so the supremum will be realized at n = 1. We therefore only need to consider f1 (Z ). As f1 (Z ) is decreasing in Z ≥ e|ck |t and has limZ →∞ f1 (Z ) = 0, it is always possible to choose a finite Zη,t such that f1 (Zη,t ) ≤ η. Therefore, we can find a Zη,t such that P (Nk,n (t ) ≥ Zη,t n) ≤ η for all n ≥ 1, which implies that with Γη,t = [−Zη,t , Zη,t ], inf P (Yn (t ) ∈ Γη,t ) ≥ 1 − η,
n ≥1
and hence η
inf P (Yn (t ) ∈ Γη,t ) ≥ 1 − η,
n ≥1
as required.
584
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Theorem 3.8.6 [21, p. 137] now requires us, for T > 0 to find β > 0 and a sequence of nonnegative random variables {νn (δ) : 0 < δ < 1} satisfying certain conditions. After noting Remark 3.8.7 [21, p. 138], it suffices to show that E [qβ (Yn (t + u), Yn (t ))|Ftn ] ≤ E [νn (δ)|Ftn ], and lim sup E [νn (δ)] = 0,
δ→0
n
where q is the metric r ∧ 1. Let β = 1 and, as before, let k = argmaxi∈S |ci | and νn (u) = Lk,n (u) =
Nk,n (u) . n
Therefore,
E [q(Yn (t + u), Yn (t ))|Ft ] ≤ E [r (Yn (t + u), Yn (t ))|Ft ] n
n
≤ = = ≤
E [r (Lk,n (t + u), Lk,n (t ))] E [r (Lk,n (u), Lk,n (0))],
(time-homogeneity)
E [νn (u)] E [νn (δ)],
since νn (δ) ≥ νn (u) for all u ≤ δ . Now E [νn (δ)] = |ck |δ , independent of n (see Lemma 1) and so lim sup E [νn (δ)] = lim |ck |δ = 0.
δ→0
δ→0
n
Therefore, all the conditions of Theorem 3.8.6 [21, p. 137] hold and hence {Yn (·)} is relatively compact.
Theorem 6.
{(Yn (·), ϕn (·))} ⇒ (X (·), ϕ(·)) in DR×S [0, ∞). Proof. Theorem 3.7.8(b) [21, p. 131] shows that we only need to prove that the finite dimensional distributions of {(Yn (·), ϕn (·))} converge in distribution to those of (X (·), ϕ(·)), since we have already shown that {(Yn (·), ϕn (·))} is relatively compact in DR×S [0, ∞). Proof of weak convergence of the finite dimensional distributions follows immediately from Theorem 4 and the fact that the processes on the relevant intervals are conditionally independent. 6. Connection with the fluid generator Q and matrix 9 In this section we demonstrate how the fundamental matrices of the SFM and the QBD are related. Assume that the parameters ϑi (1x), i ∈ S , are of the form (5). Consider the process {(X (t ), ϕ(t )), t ≥ 0}. Let θ (x) = inf{t > 0 : X (t ) = x}. For all i ∈ S1 and j ∈ S2 , let
Ψi,j = P (θ (0) < ∞, ϕ(θ (0)) = j|X (0) = 0, ϕ(0) = i)
(8)
and define the matrix 9 = [Ψi,j ]. The matrix 9 is probably the key performance measure in the theory of SFMs [5,6,9] and is the minimum nonnegative solution of the Riccati equation [6] Q12 + Q11 X + XQ22 + XQ21 X = 0,
(9)
where 1 −1 Q11 = C− T01 ], 1 [T11 − T10 (T00 ) 1 −1 Q22 = C− T02 ], 2 [T22 − T20 (T00 ) 1 −1 Q12 = C− T02 ], 1 [T12 − T10 (T00 ) 1 −1 Q21 = C− T01 ], 2 [T21 − T20 (T00 )
and
Q=
Q11 Q21
Q12 Q22
is the key fluid generator in the theory of the SFMs [6]. Further, consider the process {(X1x (t ), ϕ1x (t )), t ≥ 0}. Let τ (k) = inf{t > 0 : X1x (t ) = k}. For all i, j ∈ S , let G(1x)i,j = P (τ (k − 1) < ∞, ϕ1x (τ (k − 1)) = j|X1x (0) = k, ϕ1x (t ) = i),
(10)
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
585
and define the matrix G(1x) = [G(1x)i,j ], the key G-matrix in the theory of QBDs [15]. We partition G(1x) according to S1 ∪ S2 ∪ S0 as follows:
G(1x) =
0 0 0
G(1x)12 G(1x)22 G(1x)02
0 0 . 0
Below, we establish a result that connects the matrices G(1x) and 9. Theorem 7.
9
0 lim G(1x) = 0 1x→0+ 0
I
(−T00 )−1 (T01 9 + T02 )
0 0 . 0
Proof. For notational convenience, let A(−1) = A(−1) (1x), A(0) = A(0) (1x) and A(+1) = A(+1) (1x). Partition A(−1) , A(0) and A(+1) according to S1 ∪ S2 ∪ S0 as
0 A(−1) = 0 0
0
(−1)
A22 0
0 0 , 0
(0)
A11 0) A(0) = A(21 (0) A01
(0)
(0)
A12 (0) A22 (0) A02
A10 (0) A20 , (0) A00
(+1)
A11 A(+1) = 0 0
0 0 0
0 0 , 0
and define the matrices (0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
˜ 11 = A11 − A10 (A00 )−1 A01 , A ˜ 22 = A22 − A20 (A00 )−1 A02 , A ˜ 12 = A12 − A10 (A00 )−1 A02 , A ˜ 21 = A21 − A20 (A00 )−1 A01 . A
(11)
Matrix A(0) is clearly an infinitesimal generator of an irreducible non-conservative CTMC, and so by arguments similar to (0) (0) ˜ (110) and A˜ (220) , are also such generators and therefore are invertible. those used in [6], matrices A00 and A22 , as well as matrices A By [15], G(1x) is the minimum nonnegative solution of the equation 0 = A(−1) + A(0) G(1x) + A(+1) G(1x)2 , which gives the set of equations (0)
(0)
(0)
(+1)
(0)
(0)
(0)
(−1)
(0)
(0)
(0)
0 = A11 G(1x)12 + A12 G(1x)22 + A10 G(1x)02 + A11 G(1x)12 G(1x)22 , 0 = A21 G(1x)12 + A22 G(1x)22 + A20 G(1x)02 + A22 ,
(12) (13)
0 = A01 G(1x)12 + A02 G(1x)22 + A00 G(1x)02 .
(14)
By (14), we have (0)
(0)
(0)
(0)
(0)
(−1)
G(1x)02 = (−A00 )−1 (A01 G(1x)12 + A02 G(1x)22 ),
(15)
and so by (13),
˜ 22 )−1 (A˜ 21 G(1x)12 + A22 ). G(1x)22 = (−A
(16)
Note that (0)
(0)
A11 = T11 − V11 ,
A12 = T12 ,
(0)
A10 = T10 ,
(+1)
= V11 ,
(−1)
= V22 ,
A11
where V11 = diag(ϑi (1x))i∈S1 , (0)
(0)
A22 = T22 − V22 ,
A21 = T21 ,
(0)
A20 = T20 ,
where V22 = diag(ϑi (1x))i∈S2 , and (0)
A00 = T00 ,
(0)
A01 = T01 ,
Consequently, since
1xV11 = C1 ,
1xV22 = C2 ,
(0)
A02 = T02 .
A22
586
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
1 by left-multiplying both sides of Eq. (12) by C− 1 and using (15) and (16), it follows that 1 −1 −1 −1 0 = C− T01 G(1x)12 1 (T11 − V11 )G(1x)12 + C1 T12 G(1x)22 + C1 T10 (−T00 )
1) 1 1 −1 ˜ (0) −1 A˜ (210) G(1x)12 + A(− + C− T02 G(1x)22 + C− , 1 T10 (−T00 ) 1 V11 G(1x)12 (−A22 ) 22 and so
0 = Q12 G(1x)22 + Q11 G(1x)12 + G(1x)12
1) (−A˜ (220) )−1 A(− −I 22 1x
+ G(1x)12
(−A˜ (220) )−1 A˜ (210) 1x
G(1x)12 .
(17)
At this point we could observe that lim1x→0+ G(1x)22 = I and then note that lim1x→0+ G(1x)12 satisfies Eq. (9). However, we have no way of knowing which solution it might be to Eq. (9) (beyond being nonnegative) and so we cannot connect it to 9. Instead, we shall connect G(1x)12 directly to the minimal nonnegative solution to a perturbed form of Eq. (9) and hence show convergence to 9. By (16), we can write (17) in the equivalent form
˜ (1x)12 + Q˜ (1x)11 G(1x)12 + G(1x)12 Q˜ (1x)22 + G(1x)12 Q˜ (1x)21 G(1x)12 , 0=Q
(18)
where
˜ (1x)12 = Q12 I + Q˜ (1x)22 1x , Q ˜ (1x)11 = Q11 + Q12 Q˜ (1x)21 1x, Q 1) (−A˜ (220) )−1 A(− −I 22 , 1x (−A˜ (220) )−1 A˜ (210) = . 1x
˜ (1x)22 = Q ˜ (1x)21 Q
We now take the limit as 1x → 0+ of these coefficient matrices in (18). We have
˜ (1x)21 = lim Q
1x→0+
(0)
˜ 22 )−1 C2 Q21 lim (−1xA
1x→0+
= lim (C2 − 1xC2 Q22 )−1 C2 Q21 1x→0+
= Q21 and
1) (−1xA˜ (220) )−1 (A(− + A˜ (220) ) 22 1x→0+ = lim (C2 − 1xC2 Q22 )−1 (C2 Q22 )
˜ (1x)22 = lim Q
1x→0+
lim
1x→0+
= Q22 , and so
˜ (1x)12 = Q12 , lim Q
1x→0+
˜ (1x)11 = Q11 . lim Q
1x→0+
Clearly, the convergence of the above matrices is uniform, as they are finite-dimensional. ˜ (1x), of a SFM with fluid generator Now, consider a 9-matrix, denoted as 9
˜ (1x) = Q
˜ (1x)11 Q ˜ (1x)21 Q
˜ (1x)12 Q ˜ (1x)22 . Q
˜ (1x) is a generator, simply verify that it has nonnegative off-diagonals and all row sums equal to zero. Since, To see that Q ˜ (1x) is also the minimal nonnegative solution by above, G(1x)12 is the minimal nonnegative solution to (18), and, by [6], 9 ˜ ˜ to (18), it follows that 9(1x) = G(1x)12 . Further, since Q(1x) converges uniformly to Q, then by an argument analogous ˜ (1x) = 9. Specifically, we can derive integral expressions for 9 ˜ (i, 1x) to that in [22, Lemma 1] it follows that lim1x→0+ 9 ˜ (1x) in which exactly i peaks occur, and then let Sm1x = m ˜ corresponding to paths contributing to 9 9 ( i , 1 x ) and argue i=1
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
587
that 1x lim lim Sm ∆→0+ m→∞
˜ (1x) = lim 9 ∆→0+
= lim lim Sm1x m→∞ ∆→0+
= 9.
Consequently, lim1x→0+ G(1x)12 = 9. Also, by (16), lim G(1x)22 =
1x→0+
(0)
(0)
(−1)
˜ 22 )−1 A˜ 21 G(1x)12 + A22 lim (−A
1x→0+
1) = lim (−1xA˜ (220) )−1 1xA˜ (210) G(1x)12 + 1xA(− 22 1x→0+ −1
= C2 (0 + C2 ) = I, and so by (15), lim G(1x)02 = (−T00 )−1 (T01 9 + T02 ).
1x→0+
˜ (1x) converges to Q, it is not the most physically intuitive such perturbation. In fact, It is worth noting that although Q ˜ we regard Q(1x) as a construction to prove the above theorem, not a useful or insightful matrix in its own right. In the following result, we present a different matrix that converges to Q and has a much more appealing physical interpretation that strengthens the connection between the QBD and the SFM. Corollary 2.
Q = lim
1x→0+
1) (−A˜ (110) )−1 A(+ −I 11 ( 0 ) −1 ˜ (0) ˜ 1x (−A22 ) A21
1
(−A˜ (110) )−1 A˜ (120) (−1) (−A˜ (220) )−1 A22 −I
.
(19)
Proof. The expressions for the block matrices in the bottom row of (19) have been derived in the Proof of Theorem 7. The expressions for the block matrices in the top row follow by symmetry. Note that by (11), 1) • (−A˜ (110) )−1 A(+ is the probability matrix of the first transition from level k and the set of phases S1 ∪ S0 (with phases in 11 S2 taboo) to level (k + 1) in some phase in S1 , given the process starts from level k in some phase in S1 , • (−A˜ (110) )−1 A˜ (120) is the probability matrix of the first transition from level k and the set of phases S1 ∪ S0 (with phases in S2
taboo) to level k in some phase in S2 , given the process starts from level k in some phase in S1 ,
• (−A˜ (220) )−1 A˜ (210) is the probability matrix of the first transition from level k and the set of phases S2 ∪ S0 (with phases in S1 taboo) to level k in some phase in S1 , given the process starts from level k in some phase in S2 , and
1) • (−A˜ (220) )−1 A(− is the probability matrix of the first transition from level k and the set of phases S2 ∪ S0 (with phases in 22 S1 taboo) to level (k − 1) in some phase in S2 , given the process starts from level k in some phase in S2 .
Therefore, the right-hand side of (19) can be interpreted as the corresponding rate matrix, with respect to level, of these transitions. We note that the extension to the bounded SFMs follows by using 9, 4 [6] (a concept symmetrical to 9), fluid generator Q, and level-forward and level-reversed arguments. Here, we focused on establishing the results for matrix 9 and the fluid generator Q mainly, since these two are the key matrices in the theory of SFMs from which the other quantities can be derived. 7. Uniformization of a bounded SFM Now, let {(X (t ), ϕ(t )), t ≥ 0} be a Markovian SFM with phase variable ϕ(t ), level variable X (t ), a lower boundary X (t ) ≥ 0, upper boundary X (t ) ≤ B, and real rates ci for all i ∈ S , such that
• when 0 < X (t ) < B and ϕ(t ) = i, then the rate at which the level is changing is ci ; • when X (t ) = 0 and ϕ(t ) = i, then the rate at which the level is changing is max{ci , 0}; and • when X (t ) = B and ϕ(t ) = i, then the rate at which the level is changing is min{ci , 0}.
588
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
That is, when the SFM hits level zero, which must occur in some phase j with cj < 0, it will stay on level zero until a transition to some phase k with ck > 0 occurs. Similarly, when the SFM hits level B, which must occur in some phase j with cj > 0, it will stay on level B until a transition to some phase k with ck < 0 occurs. We choose some large positive integer K , let 1x = B/K and, as before, construct the spatially-coherent continuous-time QBD with state space G = {(k, i) : k = 0, 1, 2, . . . , K , i ∈ S } and generator T(1x) given by
T(1x) =
B(0) (1x) A(−1) (1x)
A(+1) (1x) A(0) (1x)
. ···
. ···
..
..
··· A(+1) (1x) .. .
A(−1) (1x)
··· B(K ) (1x)
,
where B(0) (1x) = [B(0) (1x)i,j ], B(K ) (1x) = [B(K ) (1x)i,j ], and
−γi (1x) B (1x)i,j = −λi
j = i, ci > 0 j = i, ci ≤ 0 j ̸= i
(0)
Ti,j
−γi (1x) (K ) B (1x)i,j = −λi
j = i, ci < 0 j = i, ci ≥ 0 j ̸= i,
Ti,j
and with the remaining block matrices as defined earlier in Section 3. We note that when the bounded QBD level is between levels 0 and K , on level 0 in phase i with ci ≥ 0, or on level K in phase i with ci ≤ 0, the process behaves just like the unbounded QBD of the earlier sections. When the QBD is on level 0 in phase i with ci < 0, or on level K in phase i with ci > 0, the process behaves just like the unbounded QBD in phase i with ci = 0. Under this interpretation it is clear that Theorem 1, Lemma 1 and Theorem 2 must still apply. Further, the weak convergence arguments in Section 5 can all be reworked in this case. Alternatively, we can apply the Continuous Mapping Theorem [21, Corollary 3.1.9, p. 103] on noting that the transformation from the unbounded SFM to the bounded SFM (and similarly for the spatially-coherent QBD approximations) is continuous. 8. Numerical example Consider an SFM {(X (t ), ϕ(t )), t ≥ 0}, X (t ) ∈ [0, +∞), with S = {1, 2, 3, 4, 5}, rates vector c = [ci ] given by
c= 1
4
−2
0 ,
−1
and generator
−8 1 T= 1 1 1
1 −6 1 1 1
4 2 −5 2 1
2 2 2 −5 1
1 1 1 . 1 −4
Calculations reveal that the drift for this process is equal to µ = −0.1524, and so the process is positive recurrent. Also, we have
−7.7500 0.3125 Q= 0.6250 1.2500
1.2500 −1.4375 0.6250 1.2500
4.2500 0.5625 −2.3750 2.2500
2.2500 0.5625 1.1250 −4.7500
(20)
and
9=
0.6854 0.6850
0.3146 . 0.3150
(21)
Consider {(X1x (t ), ϕ1x (t )), t ≥ 0} with ϑi (1x) given by (5) and 1x = 1/n. Denote by 9(n) the approximation of 9 based on Theorem 7, calculated using MATLAB. We illustrate the convergence of 9(n) to 9 in Figs. 1–2 and note, because both 9 and 9(n) are conservative generators, that the errors must necessarily sum to zero across each row. The convergence is such that a log–log plot is almost linear with asymptotic slope of −1, suggesting that the error has the form of O(n−1 ). We also present simulations to enable comparison between the fluid model (X (t ), ϕ(t )) and the approximating scaled QBD (Yn (t ), ϕ(t )), for various values of n. In order to simulate these processes in a fully-coupled manner, the phase transitions (i.e. the process ϕ(t )) were sampled once and kept the same for all the processes. Of course, that is the only source of randomness in the SFM and so the process (X (t ), ϕ(t )) was calculated directly from this. For the QBD approximations,
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
589
Fig. 1. 91j − 9(n)1j with scale 10−4 (y-axis).
Fig. 2. 92j − 9(n)2j with scale 10−5 (y-axis).
once the residence time is known for state i, say τ , the level change in the QBD (for each value of n) during this residence time is then given by a Poisson distribution with parameter ϑi (1/n)τ . We, therefore, sampled from the necessary Poisson distributions and built the QBD level accordingly. This is the most efficient way to simulate the QBD, and also enables the required coupling. In Fig. 3 we show the first 30 phases in the given sample path. On the graph we give the value of the fluid X (t ) and for two (scaled) QBD approximations Yn (t ) with n = 102 and 104 . As can clearly be seen, with n = 102 the approximation is very good, and once we make n = 104 it is essentially indistinguishable from the fluid. To explore this more fully, we present Figs. 4 and 5 that show the difference between the processes X (t ) and Yn (t ) for values of n ranging through 10, 102 , 103 , 104 , 105 and 106 . It is clear that each time the value of n is increased by an order of magnitude, there is a significant reduction in the magnitude of this difference. 9. Conclusion In this paper we have developed a natural spatial-uniformization of the Stochastic Fluid Model (SFM) into a QBD. The approach taken follows the standard uniformization argument, with a further discretization of the continuous fluid space. By reinstating continuous time, we developed a continuous-time QBD whose level naturally approximates the level of the SFM. This differs from the two previous QBDs derived from an SFM, where the level had no direct relationship to the SFM [9], or which required a further calculation involving the workload remaining in the queue of the QBD [2], and where the level represented, instead, the number of customers in the queue. We were able to show formal weak convergence of the appropriately scaled approximating QBDs to the SFM and that the generator of the fluid Q and the key matrix 9 are the limits of the natural QBD quantities. More importantly, this leaves us with a discrete-space approximation to the SFM, which we can use in the development of numerical schemes for the analysis of the stochastic fluid–fluid model, constructed in [18].
590
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Fig. 3. The fluid level X (t ) and two approximating scaled QBD levels Yn (t ) for n = 102 and 104 .
Fig. 4. X (t ) − Yn (t ) for n = 10, 102 , 103 .
Fig. 5. X (t ) − Yn (t ) for n = 104 , 105 , 106 .
We emphasize that the aim of the paper was to construct a meaningful uniformization in which the level of the SFM is directly linked with the level of the QBD so that we can use this theoretical connection between the two models and derive
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
591
more information about the SFMs. In particular, the spatially-coherent uniformization derived in the paper opens the door to the treatment of the multi-dimensional SFMs, since discretizing one of the SFMs in the two-dimensional model results in a one-dimensional SFM which we can treat using the existing methods. This work, which requires rigorous construction and validation, is under preparation and will be reported in our upcoming paper. Another important advantage achieved, since the theory of the QBDs is well understood, is that we can use the results and insights for the QBDs in order to derive further properties of the SFM, which were not considered previously. For example, our work on the quasi-stationary distribution of the SFMs is under preparation, and we will report on this in due course. We believe that more such opportunities will arise from this construction. Acknowledgment The authors would like to thank the Australian Research Council for funding this research through Discovery Project DP110101663. Appendix. Statement of necessary theorems from Ethier and Kurtz Theorem (3.2.2, p. 104 [21]). Let (S , d) be complete and separable and let M ⊂ P (S ) (where P (S ) is the family of Borel probability measures on S). Then the following are equivalent: (a) M is tight. (b) For each ϵ > 0, there exists a compact set K ⊂ S such that inf P (K ϵ ) ≥ 1 − ϵ,
P ∈M
where K ϵ = x ∈ S : infy∈K d(x, y) < ϵ . (c) M is relatively compact.
Proposition (3.2.4, p. 107 [21]). ∞ ∞Let −(kSk , dk ), k = 1, 2, . . . , be metric spaces, and define the metric space (S , d) by letting S = S and d ( x , y ) = (dk (xk , yk ) ∧ 1) for all x, y ∈ S. Let {Pα } ⊂ P (S ) (where α ranges over some index k k=1 k=1 2 set), and for k = 1, 2, . . . and each α , define Pαk ∈ P (S ) to be the kth marginal distribution of Pα (i.e., Pαk = Pα πk−1 , where the projection πk : S −→ Sk is given by πk (x) = xk ). Then {Pα } is tight if and only if {Pαk } is tight for k = 1, 2, . . . . Theorem (3.7.2, p. 128 [21]). Let (E , r ) be complete and separable, and let {Xα } be a family of processes with sample paths in DE [0, ∞). Then {Xα } is relatively compact if and only if the following two conditions hold: (a) For every η > 0 and rational t ≥ 0, there exists a compact set Γη,t ⊂ E such that η
inf P {Xα (t ) ∈ Γη,t } ≥ 1 − η, α η
where Γη,t = x ∈ S : infy∈Γη,t d(x, y) < η . (b) For every η > 0 and T > 0, there exists δ > 0 such that
sup P {w ′ (Xα , δ), T ≥ η} ≤ η. α
Theorem (3.7.8(b), p. 131 [21]). Let E be separable and let Xn , n = 1, 2, . . . , and X be processes with sample paths in DE [0, ∞). (b) If {Xn } is relatively compact and there exists a dense set D ⊂ [0, ∞) such that
(Xn (t1 ), . . . , Xn (tk )) ⇒ (X (t1 ), . . . , X (tk )) holds for every finite set {t1 , . . . , tk } ⊂ D, then Xn ⇒ X . Theorem (3.8.6, p. 137 [21]). Let (E , r ) be complete and separable, and let {Xα } be a family of processes with sample paths in DE [0, ∞). Suppose that condition (a) of Theorem 3.7.2 holds. Then the following are equivalent: (a) {Xα } is relatively compact. (b) For each T > 0, there exist β > 0 and a family {γα (δ) : 0 < δ < 1 all α} of nonnegative random variables satisfying E qβ (Xα (t + u), Xα (t )) |Ftα qβ (Xα (t ), Xα (t − v)) ≤ E γα (δ)|Ftα
for 0 ≤ t ≤ T , 0 ≤ u ≤ δ , and 0 ≤ v ≤ δ ∧ t, where Ftα ≡ Ft α ; in addition, X
lim sup E [γα (δ)] = 0
δ→0 α
and lim sup E qβ (Xα (δ), Xα (0)) = 0.
δ→0 α
592
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Remark (3.8.7, p. 138 [21]). (a) If, as will typically be the case, E qβ (Xα (t + u), Xα (t )) |Ftα ≤ E γα (δ)|Ftα
instead of E qβ (Xα (t + u), Xα (t )) |Ftα qβ (Xα (t ), Xα (t − v)) ≤ E γα (δ)|Ftα ,
then E qβ (Xα (δ), Xα (0)) ≤ E [γα (δ)] and we need only verify
lim sup E [γα (δ)] = 0
δ→0 α
in condition (b) of Theorem 3.8.6. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
S. Ahn, V. Ramaswami, Fluid flow models and queues—a connection by stochastic coupling, Stochastic Models 19 (3) (2003) 325–348. S. Ahn, V. Ramaswami, Transient analysis of fluid flow models via stochastic coupling to a queue, Stochastic Models 20 (1) (2004) 71–101. S. Ahn, V. Ramaswami, Efficient algorithms for transient analysis of stochastic fluid flow models, Journal of Applied Probability 42 (2) (2005) 531–549. S. Asmussen, Stationary distributions for fluid flow models with or without Brownian noise, Stochastic Models 11 (1) (1995) 21–49. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Algorithms for return probabilities for stochastic fluid flows, Stochastic Models 21 (1) (2005) 149–184. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Hitting probabilities and hitting times for stochastic fluid flows, Stochastic Processes and their Applications 115 (9) (2005) 1530–1556. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Algorithms for the Laplace–Stieltjes transforms of first return times for stochastic fluid flows, Methodology and Computing in Applied Probability 10 (3) (2008) 381–408. V. Ramaswami, Matrix analytic methods: a tutorial overview with some extensions and new results, in: Matrix-Analytic Methods in Stochastic Models (Flint, MI), in: Lecture Notes in Pure and Appl. Math., vol. 183, Dekker, New York, 1997, pp. 261–296. V. Ramaswami, Matrix analytic methods for stochastic fluid flows, in: Proceedings of the 16th International Teletraffic Congress, Edinburgh, 7–11 June 1999, pp. 1019–1030. D. Anick, D. Mitra, M.M. Sondhi, Stochastic theory of a data-handling system with multiple sources, Bell System Technical Journal 61 (8) (1982) 1871–1894. A. Badescu, L. Breuer, A. Da Silva Soares, G. Latouche, M.-A. Remiche, D. Stanford, Risk processes analyzed as fluid queues, Scandinavian Actuarial Journal 105 (2) (2005) 127–141. N.G. Bean, M.M. O’Reilly, J.E. Sargison, A stochastic fluid flow model of the operation and maintenance of power generation systems, IEEE Transactions on Power Systems 25 (3) (2010) 1361–1374. N.G. Bean, M.M. O’Reilly, Performance measures of a multi-layer Markovian fluid model, Annals of Operations Research 160 (2008) 99–120. G. Latouche, P.G. Taylor, A stochastic fluid model for an ad hoc mobile network, Queueing Systems 63 (1–4) (2009) 109–129. G. Latouche, V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, in: ASA-SIAM Series on Statistics and Applied Probability, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1999. N.G. Bean, G.G. Latouche, G.T. Nguyen, The stochastic fluid–fluid model applied to a two-server queue with assistance, 2013, in preparation. N.G. Bean, M.M. O’Reilly, A stochastic two-dimensional fluid model, Stochastic Models 29 (1) (2013) 31–63. N.G. Bean, M.M. O’Reilly, The stochastic fluid–fluid model: a stochastic fluid model driven by an uncountable-state process, which is a stochastic fluid model itself, 2013, submitted for publication. S.M. Ross, Introduction to Probability Models, Elsevier, New York, 2007. P. Billingsley, Convergence of Probability Measures, in: Wiley Series in Probability and Mathematical Statistics: Tracts on Probability and Statistics, Wiley, 1968. S.N. Ethier, T.G. Kurtz, Markov Processes Characterization and Convergence, in: Wiley Series in Probability and Mathematical Statistics, John Wiley & Sons, 1986. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Hitting probabilities and hitting times for stochastic fluid flows: the bounded model, Probability in the Engineering and Informational Sciences 23 (1) (2009) 121–147.
Nigel G. Bean received the B.Sc. degree in mathematics and the Honours degree in applied mathematics from the University of Adelaide, Adelaide, Australia, in 1988 and 1989, respectively. He received the Ph.D. degree from the University of Cambridge, Cambridge, U.K., in 1993, where he studied under Prof. F. Kelly, FRS. Since then he has been employed at the University of Adelaide and was appointed Chair of Applied Mathematics in 2004. His main research interests are in stochastic modeling, particularly Markov chains, and applied operations research. Prof. Bean is the recipient of the 2001 J.H. Michell Medal, awarded by the Australian and New Zealand Industrial and Applied Mathematics (ANZIAM) division of the Australian Mathematical Society, and the 2003 P.A.P. Moran Medal, awarded by the Australian Academy of Science.
Małgorzata M. O’Reilly received the M.Sc. degree in mathematics and education from Wroclaw University, Wroclaw, Poland, in 1987, and the Ph.D. degree in applied probability from the University of Adelaide, Adelaide, Australia, in 2002. She has been a lecturer in probability and operations research at the University of Tasmania, Hobart, Tasmania, Australia, since 2005. Her research is in the area of stochastic modeling, Markov chains, and operations research.
Contents lists available at ScienceDirect
Performance Evaluation journal homepage: www.elsevier.com/locate/peva
Spatially-coherent uniformization of a stochastic fluid model to a Quasi-Birth-and-Death process Nigel G. Bean a , Małgorzata M. O’Reilly b,∗ a
Applied Mathematics, University of Adelaide, SA 5005, Australia
b
School of Mathematics, University of Tasmania, Tas 7001, Australia
article
info
Article history: Available online 24 June 2013 Keywords: Stochastic fluid model Quasi-Birth-and-Death process Markov chain Uniformization
abstract We derive a uniformization of a stochastic fluid model (SFM) to a Quasi-Birth-and-Death process (QBD) that is spatially-coherent since the continuous level in the SFM has a natural correspondence to the discrete level in the QBD. As a consequence of this, the QBD can be used as a direct approximation of the original SFM, in those situations where a discrete state space is an advantage. We treat the unbounded as well as the bounded cases and illustrate the theory with a numerical example. The key fluid generator, Q, and matrix Ψ for the SFMs emerge from the QBD calculations in the natural limit. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Let {(ϕ(t )), t ≥ 0} be an irreducible continuous-time Markov Chain (CTMC) with a finite state space S = {1, 2, . . . , n} and infinitesimal generator T = [Ti,j ]. Let {(X (t ), ϕ(t )), t ≥ 0} be a Markovian stochastic fluid model (SFM) with phase variable ϕ(t ) and unbounded level variable X (t ) ∈ (−∞, +∞). There are real rates ci for all i ∈ S , such that when ϕ(t ) = i, then the rate at which the level is changing is ci . The buffer collecting fluid in this process is denoted by X . The CTMC is referred to as the driving process. The literature on SFMs has focused on the stationary and transient analysis, as well as the development of efficient algorithms [1–9] and performance-driven applications [10–14]. The first paper that used the class of effective methods known as matrix-analytic methods [15] to consider a Markovian SFM was the seminal paper by Ramaswami [9]. He proposed one way of mapping an SFM to a Quasi-Birth-and-Death process (QBD). However, that mapping is not spatially-coherent, in the sense that the level of the QBD does not correspond to the level of the SFM in any meaningful way. Instead, it was derived with the specific purpose of allowing exact algorithmic evaluation of the key matrix 9 for the SFMs, using known QBD techniques. Bean, O’Reilly and Taylor [6] developed a direct method for the analysis of SFMs which uses the key fluid generators Q(s) and does not require the construction of a QBD. In the work by Ahn and Ramaswami [2], Ramaswami’s mapping was further developed in such a way that the workload in the queue (not the level of the queue) represented the level of the SFM in an appropriate limit. In this paper we derive a spatially-coherent uniformization of an SFM to a QBD. This preserves a natural mapping between the continuous level in the SFM and the discrete level in the QBD, and so the QBD can be used as a direct approximation of the original SFM. We first construct a spatially-coherent discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} through a uniformization of the SFM. This construction yields a discrete-time QBD, but by incorporating the underlying Poisson process directly back into it, we develop the spatially-coherent continuous-time QBD {(X1x (t ), ϕ1x (t )), t ≥ 0}. A new class of important models, referred to as multi-dimensional SFMs [16–18,14], has recently emerged which are very hard to analyze. In particular, we are interested in the stochastic fluid–fluid model constructed in [18]. One attractive
∗
Corresponding author. E-mail addresses: [email protected] (N.G. Bean), [email protected] (M.M. O’Reilly).
0166-5316/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.peva.2013.05.006
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
579
approach to the analysis of these models is to discretize the first fluid and then treat the second fluid as being driven by a QBD. Therefore, since the QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} can be used as an approximation of the SFM, the significance of these results goes beyond the development of new ways of approximating the performance measures of an SFM. Instead it is possible to use this direct connection to translate known results, and methods of analysis, for QBDs to the field of SFMs. In particular, we expect to exploit this crucial feature in future analyses of multi-dimensional SFMs and thus provide a matrix-analytic methods based analysis of this important, but difficult to treat, class of models. We formalize the spatially-coherent uniformization in our paper and show that the transition function of the driving process behind the SFM is the same as that of the continuous-time QBD, following an argument similar to the usual uniformization argument. Hence, this QBD should be a good approximation for the SFM as it captures the stochastic phase-dynamics correctly; all we then need is to capture the level-dynamics up to the resolution allowed by the particular choice of 1x. By using a fairly natural choice of the uniformization parameter we are able to show that the level-dynamics converge in probability to those of the fluid. We then show that, in the limit as 1x → 0+ , the QBDs converge to the SFM (weak convergence of the processes), and that the key fluid generator Q and matrix 9 also emerge from the QBD calculations in this limit. We treat the unbounded as well as the bounded cases and illustrate the theory with a numerical example. The structure of the paper is as follows. First, we describe the uniformization method with general parameters ϑi (1x) > 0, i ∈ S . In Section 2 we construct the discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} and in Section 3 the continuous-time QBD {(X1x (t ), ϕ1x (t )), t ≥ 0}, for the unbounded SFM with X (t ) ∈ (−∞, +∞). Next, we focus on the appropriate choice of parameters ϑi (1x) in Section 4, and analyze the corresponding properties. The proofs of the main convergence results are given in Section 5. In Section 6 we show the connection with the fluid generator Q and matrix 9. The uniformization of the bounded SFM with X (t ) ∈ [0, B] is presented in Section 7. A numerical example is given in Section 8, followed by concluding remarks in Section 9.
¯ 1x (n)), n ≥ 0} 2. The spatially-coherent discrete-time QBD {(X¯ 1x (n), ϕ We partition the set of all phases as S = S1 ∪ S2 ∪ S0 , where S1 = {i ∈ S : ci > 0}, S2 = {i ∈ S : ci < 0}, S0 = {i ∈ S : ci = 0}, and the generator as T11 T21 T01
T=
T12 T22 T02
T10 T20 , T00
according to the partitioning of S . Further, let C1 = diag(ci )i∈S1 and C2 = diag(|ci |)i∈S2 . For a given positive 1x > 0, and for all i ∈ S , let
ϑi (1x) > 0
(1)
be some function of 1x. For the moment, we do not assume any particular form of ϑi (1x), except the positivity assumption. Later, in Section 4 we define a suitable form. We construct {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0}, the spatially-coherent discrete-time QBD with state space G¯ = {(k, i) : k ∈ Z, i ∈ S }, by connecting the phase and level of this QBD to the state of the SFM at arrival times of a Poisson process with an appropriately chosen parameter, following the standard uniformization approach [19]. Let
γi (1x) = λi + ϑi (1x),
(2)
where λi = Ti,i , for all i ∈ S . Choose
γ (1x) ≥ max {γi (1x)} ,
(3)
i
and consider the Poisson process {N (t ), t ≥ 0} with parameter γ (1x). X (0) The initial condition is (X¯ 1x (0), ϕ¯ 1x (0)) = (⌊ 1x ⌋, ϕ(0)). Then, by the standard theory of Markov Chains, for example
[19], assuming that an arrival in N (·) occurs at time t and that (X¯ 1x (N (t −)), ϕ¯ 1x (N (t −))) = (k, i),
• with probability p1 = λi /γ (1x) it is marked as a type-1 arrival and the process moves to some phase j ̸= i without a change in the discrete level so that a transition from (k, i) to (k, j) occurs (where j ̸= i is chosen with probability Ti,j /λi ), • with probability p2 = ϑi (1x)/γ (1x) it is marked as a type-2 arrival and the process does not change phase but changes level so that a transition from (k, i) to (k + 1, i) occurs if ci > 0, or a transition from (k, i) to (k − 1, i) occurs if ci < 0, or • with probability p3 = 1 − p1 − p2 it is marked as a type-3 arrival and the process remains in state (k, i). ¯ (1x) = [P¯ (1x)(k,i)(m,j) ] is given by Consequently, the corresponding one-step transition probability matrix P 1 − γ (1x)/γ (1x) m = k, j = i, c ̸= 0 i i 1 − λi /γ (1x) m = k, j = i, ci = 0 Ti,j /γ (1x) m = k, j ̸= i P¯ (1x)(k,i)(m,j) = m = k + 1, j = i, ci > 0 ϑi (1x)/γ (1x) m = k − 1, j = i, ci < 0 ϑi (1x)/γ (1x) 0
otherwise,
580
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
and can be written in the form
..
..
(−1) . ¯ ¯ (1x) = A (1x) P ¯ (0)
.
¯ (0) (1x) A ¯A(−1) (1x) ..
..
.
¯ (+1) (1x) A ¯ (0) (1x) A ..
.
¯ (0)
.
¯ (−1)
¯ (+1) (1x) A ..
.
,
¯ (−1)
where matrices A (1x) = [A (1x)i,j ]i,j∈S , A (1x) = [A (1x)i,j ]i,j∈S and A¯ (+1) (1x) = [A¯ (+1) (1x)i,j ]i,j∈S corresponding to the transitions to the same level, one level down and one level up, respectively, are given by
1 − γi (1x)/γ (1x) j = i, ci ̸= 0 ¯A(0) (1x)i,j = 1 − λi /γ (1x) j = i, ci = 0 Ti,j /γ (1x) j ̸= i ϑi (1x)/γ (1x) j = i, ci < 0 A¯ (−1) (1x)i,j = 0 otherwise ϑi (1x)/γ (1x) j = i, ci > 0 A¯ (+1) (1x)i,j = 0 otherwise. 3. The spatially-coherent continuous-time QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} In this section, using the discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} and standard techniques [19], we essentially reverse the uniformization to develop a continuous-time QBD. Specifically, we fix 1x > 0 and embed the events of the discrete-time QBD {(X¯ 1x (n), ϕ¯ 1x (n)), n ≥ 0} at the time points of a Poisson process of rate γ (1x) to create the continuoustime Markov chain {(X1x (t ), ϕ1x (t )), t ≥ 0}. This CTMC then has state space G = {(k, i) : k ∈ Z, i ∈ S }, levels k ∈ Z, phases i ∈ S , and generator T(1x) = [T (1x)(k,i)(m,j) ] given by
T (1x)(k,i)(m,j) =
−γ (1x) i −λi
m = k, j = i, ci ̸= 0 m = k, j = i, ci = 0 m = k, j ̸= i m = k + 1, j = i, ci > 0 m = k − 1, j = i, ci < 0 otherwise.
Ti,j
ϑi (1x) ϑi (1x) 0
Clearly T(1x) is well defined. It has negative on-diagonals, nonnegative off-diagonals, and all rows sums equal to zero. We note that the above CTMC has the structure of an unbounded continuous-time QBD with
..
..
(−1) . A (1x) T(1x) = (0)
. (0) A (1x) A(−1) (1x) .. .
..
. (+1) A (1x) A(0) (1x) .. .
(0)
A(+1) (1x)
..
.
,
where matrices A (1x) = [A (1x)i,j ]i,j∈S , A (1x) = [A (1x)i,j ]i,j∈S and A(+1) (1x) = [A(+1) (1x)i,j ]i,j∈S corresponding to the transitions to the same level, one level down and one level up, respectively, are given by
−γi (1x) A (1x)i,j = −λi (0)
Ti,j A(−1) (1x)i,j =
A(+1) (1x)i,j =
ϑi (1x) 0
ϑi (1x) 0
(−1)
(−1)
j = i, ci ̸= 0 j = i, ci = 0 j ̸= i j = i, ci < 0 otherwise j = i, ci > 0 otherwise.
Now, assume (X1x (0), ϕ 1x (0)) = (k, i) for some k ∈ Z and i ∈ S . Let τi be the random variable that records the total time until a transition to some (m, j) with j ̸= i occurs. That is, τi is the time until the process {(X1x (t ), ϕ1x (t )), t ≥ 0} changes phase (regardless of level). Further, for all i, j ∈ S and t, define P (t )ij = P (ϕ(t ) = j|ϕ(0) = i), P1x (t )ij = P (ϕ1x (t ) = j|ϕ1x (0) = i). Below we establish a result which holds for any parameters ϑi (1x) > 0, i ∈ S , and guarantees the equivalence of the phase process in the uniformized QBD {(X1x (t ), ϕ1x (t )), t ≥ 0}, to the phase process in the original SFM, {(X (t ), ϕ(t )), t ≥ 0}.
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
581
Theorem 1. For any 1x > 0, the process {(X1x (t ), ϕ1x (t )), t ≥ 0} satisfies the following properties: (i) The random variable τi follows the exponential distribution with parameter λi . (ii) The probability that the process moves to a state with phase j, given that it leaves phase i (regardless of level), is Tij /λi . Therefore, for all i, j ∈ S and t, P (t )ij = P1x (t )ij .
(4)
Proof. We apply classic arguments from the theory of Markov Chains [19]. For i with ci = 0 the result follows trivially, since then, for any k ∈ Z, we have T(k,i)(k,i) = −λi and the only possible way for the process to leave phase i is to move from (k, i) to (k, j) for some j ̸= i. Suppose that ci ̸= 0. Without loss of generality assume ci > 0. First, note that the distribution of time spent in (m, i) is exponential with parameter γi (1x). An arrival of type 1, that is from (k, i) to (m, j) with j ̸= i can only occur if m = k. Also, note that for i with ci > 0 the only possible transitions during which a phase does not change are arrivals of type 2, that is from (m, i) to (m + 1, i) for some m. Now, given that the process makes n − 1 transitions before changing phase i for the first time at transition n, the distribution of τi is Erlang of order n and with parameter γi (1x). It therefore has density function given by
(γi (1x)t )n−1 γi (1x)e−γi (1x)t . (n − 1)! The probability that the process makes n − 1 transitions before changing phase i for the first time at transition n is equivalent to the probability that the process makes n − 1 transitions of type 1 (from (m, i) to (m + 1, i) for some m), and then one transition of type 2 (from (m, i) to (m, j) for some j ̸= i), which is equal to
λi 1− γi (1x)
n −1
λi . γi (1x)
Hence, by the Law of Total Probability, the density function of τi is given by ∞ (γi (1x)t )n−1 n =1
(n − 1)!
γi (1x)e−γi (1x)t 1 −
λi γi (1x)
n−1
λi γi (1x)
= λi e−γi (1x)t eγi (1x)t (1−λi /γi (1x)) = λi e−λi t , and so property (i) follows. Next, the probability that the process moves to a state with phase j ̸= i, given that a transition from (k, i) to a state with phase different than i occurs, is equal to Tij /γi (1x)
λi /γi (1x)
=
Tij
λi
,
which shows property (ii). Properties (i) and (ii) show that a CTMC obtained by observing only the phase in the QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} has the same probability law as the CTMC driving the original SFM {(X (t ), ϕ(t )), t ≥ 0}, before uniformization. Hence (4) follows. We are now interested in a suitable choice of parameters ϑi (1x) so as to establish a meaningful physical connection between the level of the QBD and the level of the SFM. Assume that the QBD and the SFM both start from level zero and in the same phase at time zero, and that ϕ(·) and ϕ1x (·) are coupled, in the sense that they share the same sample path. We desire the property such that at all times t ≥ 0, if the QBD {(X1x (t ), ϕ1x (t )), t ≥ 0} is on level k, then the SFM {(X (t ), ϕ(t )), t ≥ 0} is near the level k1x with probability nearly 1. We formalize this idea in Section 4, followed by the proof of weak convergence of the processes in Section 5. 4. Level properties In this section, we assume for all i ∈ S that the parameters ϑi (1x) are of the form
ϑi (1x) =
|ci | , 1x
∀1x > 0,
(5)
and study the level-based properties of the process {(X1x (t ), ϕ1x (t )), t ≥ 0}. Consider the Poisson process {(Ni,1x (t )), t ≥ 0} with rate ϑi (1x) for every i ∈ S . Let Li,1x (t ) be the random variable defined by Li,1x (t ) = Ni,1x (t )1x. The following results are immediate from the properties of the underlying Poisson process.
582
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Lemma 1. For all 1x > 0 and t ≥ 0, E (Li,1x (t )) = |ci |t , and, for all T ≥ 0, lim Var(Li,1x (t )) = 0,
1x→0+
where the convergence is uniform on the interval t ∈ [0, T ]. Let A(t ) be the event that a change in phase does not occur during the time interval [0, t ] (but any number of changes in level may occur). Corollary 1. For all 1x > 0 and t ≥ 0, E (X1x (t )1x|(X1x (0), ϕ1x (0)) = (0, i), A(t )) = ci t , and, for all t ≥ 0, lim Var(X1x (t )1x|(X1x (t ), ϕ1x (0)) = (0, i), A(t )) = 0.
1x→0+
Theorem 2. For any ϵ > 0, i ∈ S , and T > 0, the family {Fi,1x (t )} of probabilities Fi,1x (t ) = P (Li,1x (t ) − |ci |t ≥ ϵ),
converges uniformly to 0, as 1x → 0+ , on the interval t ∈ [0, T ]. Proof. The result follows immediately from Lemma 1, since Var(Li,1x (t )) = E ((Li,1x (t ) − |ci |t )2 ), upon use of Markov’s inequality. 5. Weak convergence For all n = 1, 2, . . . , define Yn (·) = X1x (·)1x
and
ϕn (·) = ϕ1x (·) where 1x =
1 n
.
(6)
Assume without loss of generality that Yn (0) = 0 and X (0) = 0. Also, assume that parameters ϑi (1x) = ϑi (1/n), i ∈ S , are of the form (5). Consider the phase process {ϕ(·)} and {ϕn (·), n ≥ 1}. Let Ft be the σ -algebra for the process ϕ(·) up to and including time t and Ftn be the σ -algebra for the process ϕn (·) up to and including time t, for all n ≥ 1. Below, we establish the weak convergence of {Yn (·), ϕn (·)} ⇒ {X (·), ϕ(·)} as n → ∞. We need to work with two spaces R and S , and with the two associated Euclidean metrics r and s, making (R, r ) and (S , s) both metric spaces. In fact, they are both separable and complete metric spaces and, therefore, so is the product space (R, r ) × (S , s). Finally, both the SFM and the QBD have sample paths in DR×S [0, ∞) and we equip this space with the usual Skorokhod J1 -topology. For the technical details see Billingsley [20, Chapter 3] and Ethier and Kurtz [21, Chapter 3]. The structure of the argument in this section is as follows. Theorem 4 (which relies on Theorems 1 and 3) proves that for any individual t ≥ 0 that the distribution {Yn (t ), ϕn (t )} converges weakly to the distribution {Y (t ), ϕ(t )}; Theorem 5 proves that every sub-sequence of the sequence of processes {Yn (·), ϕn (·)} contains a further sub-sequence that is weakly convergent; finally, Theorem 6 shows that the finite-dimensional distributions converge weakly, which, in essence, shows that the process limit that was shown to exist in Theorem 5 must be unique since the sequence of processes is relative compact. This is sufficient to prove that the sequence of processes {Yn (·), ϕn (·)} weakly converges to the process {Y (·), ϕ(·)}. The theorems cited in the proofs are given in the Appendix, for completeness. Theorem 1 shows that ϕn (·) obeys the same probability law as ϕ(·). We shall therefore couple the processes by assuming that ϕn (·) and ϕ(·) share the same sample path, and so ϕn (t ) = ϕ(t ) for all t ≥ 0. Hence Ft ≡ Ftn . P
Theorem 3. {Yn (T )|FT } ⇒ X (T )|FT as n → ∞ and for all T ≥ 0. Proof. Given FT , define ui as the total time spent in phase i during the time interval [0, T ], for all i ∈ S . Therefore,
ui = T
i∈S
and X (T ) =
i∈S
Further, given FT , it is also clear that Yn (T ) =
i∈S
sign(ci )Li,n (ui ).
ci ui .
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
583
In Theorem 2 we showed that P (Li,n (t ) − |ci |t ≥ ϵ) → 0 uniformly in t ∈ [0, T ]. Therefore, since there are only finitely many phases |S |, ϵ P Li,n (ui ) − |ci |ui ≥ P Yn (T ) − X (T ) ≥ ϵ|FT ≤ |S | i∈S which converges to 0. Hence we have shown convergence in probability for the conditional probabilities.
D
Theorem 4. {(Yn (t ), ϕn (t ))} ⇒ (X (t ), ϕ(t )) as n → ∞ and for all t ≥ 0. Proof. For any measurable set A ⊂ R and j ∈ S , Theorem 1 shows that P (Yn (t ) ∈ A, ϕn (t ) = j) = P (Yn (t ) ∈ A, ϕ(t ) = j). Therefore, lim P (Yn (t ) ∈ A, ϕn (t ) = j) = lim
n→∞
n→∞
=
P (Yn (t ) ∈ A, ϕ(t ) = j|Ft )dP (Ft )
P (X (t ) ∈ A, ϕ(t ) = j|Ft )dP (Ft ),
by Theorem 3 since we have shown convergence in probability. The result then follows immediately.
Theorem 5. The sequence of processes {Yn (·), ϕn (·)} is relatively compact in DR×S [0, ∞). Proof. To prove relative compactness of {Yn (·), ϕn (·)} it suffices to prove relative compactness for both of the sequences {Yn (·)} and {ϕn (·)} separately [21, Theorem 3.2.2, p. 104 and Proposition 3.2.4, p. 107]. Now ϕn (t ) ∈ {1, 2, . . . , m}, for all t and so lives in a compact space. Prohorov’s Theorem [21, Theorem 3.2.2, Page 104] immediately implies that {ϕn (·)} is relatively compact. We are thus left to show that {Yn (·)} is relatively compact. We use Theorem 3.8.6 [21, p. 137] which first requires us to show that condition (a) of Theorem 3.7.2 [21, p. 128] holds. Condition (a) of Theorem 3.7.2 is that for every η > 0 and rational t ≥ 0, there exists a compact set Γη,t ∈ R such that η
inf P {Yn (t ) ∈ Γη,t } ≥ 1 − η,
(7)
n ≥1
η
where Γη,t = x ∈ S : infy∈Γη,t d(x, y) < η . Note that Eq. (7) is essentially equivalent to asserting tightness of the random variable Yn (t ), but only for rational t. To show that this holds, we shall define Γη,t = [−Z , Z ], Z > 0, and then prove the existence of a suitable Z . If Z is finite, then Γη,t is a compact set. Now, in order for Yn (t ) to escape the interval [−Z , Z ] there has to be at least Zn changes of level in the QBD. Let k = argmaxi∈S {|ci |} and so phase k represents the phase with the highest rate of level changes. Therefore,
P (Yn (t ) ̸∈ [−Z , Z ]) ≤ P (QBD has at least Zn changes of level in time t )
≤ P (Nk,n (t ) ≥ Zn). Recall that Nk,n (·) is the underlying Poisson process corresponding to the fastest rate of level changes in the QBD, and any move into a different phase will only lead to a slower Poisson process or possibly a phase that moves in the opposite direction. Now, we show that a suitable Z exists for all t and all η > 0. The Chernoff bound for the Poisson distribution implies that P (Nk,n (t ) ≥ Zn) ≤
e−β (eβ)Zn
(Zn)Zn
,
where β = ϑk (1/n)t = |ck |nt. We therefore need to choose Z such that fn (Z ) =
e−|ck |nt (e|ck |nt )Zn
(Zn)Zn
≤ η,
∀n ≥ 1.
If Z ≥ e|ck |t, then fn (Z ) is decreasing in n and so the supremum will be realized at n = 1. We therefore only need to consider f1 (Z ). As f1 (Z ) is decreasing in Z ≥ e|ck |t and has limZ →∞ f1 (Z ) = 0, it is always possible to choose a finite Zη,t such that f1 (Zη,t ) ≤ η. Therefore, we can find a Zη,t such that P (Nk,n (t ) ≥ Zη,t n) ≤ η for all n ≥ 1, which implies that with Γη,t = [−Zη,t , Zη,t ], inf P (Yn (t ) ∈ Γη,t ) ≥ 1 − η,
n ≥1
and hence η
inf P (Yn (t ) ∈ Γη,t ) ≥ 1 − η,
n ≥1
as required.
584
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Theorem 3.8.6 [21, p. 137] now requires us, for T > 0 to find β > 0 and a sequence of nonnegative random variables {νn (δ) : 0 < δ < 1} satisfying certain conditions. After noting Remark 3.8.7 [21, p. 138], it suffices to show that E [qβ (Yn (t + u), Yn (t ))|Ftn ] ≤ E [νn (δ)|Ftn ], and lim sup E [νn (δ)] = 0,
δ→0
n
where q is the metric r ∧ 1. Let β = 1 and, as before, let k = argmaxi∈S |ci | and νn (u) = Lk,n (u) =
Nk,n (u) . n
Therefore,
E [q(Yn (t + u), Yn (t ))|Ft ] ≤ E [r (Yn (t + u), Yn (t ))|Ft ] n
n
≤ = = ≤
E [r (Lk,n (t + u), Lk,n (t ))] E [r (Lk,n (u), Lk,n (0))],
(time-homogeneity)
E [νn (u)] E [νn (δ)],
since νn (δ) ≥ νn (u) for all u ≤ δ . Now E [νn (δ)] = |ck |δ , independent of n (see Lemma 1) and so lim sup E [νn (δ)] = lim |ck |δ = 0.
δ→0
δ→0
n
Therefore, all the conditions of Theorem 3.8.6 [21, p. 137] hold and hence {Yn (·)} is relatively compact.
Theorem 6.
{(Yn (·), ϕn (·))} ⇒ (X (·), ϕ(·)) in DR×S [0, ∞). Proof. Theorem 3.7.8(b) [21, p. 131] shows that we only need to prove that the finite dimensional distributions of {(Yn (·), ϕn (·))} converge in distribution to those of (X (·), ϕ(·)), since we have already shown that {(Yn (·), ϕn (·))} is relatively compact in DR×S [0, ∞). Proof of weak convergence of the finite dimensional distributions follows immediately from Theorem 4 and the fact that the processes on the relevant intervals are conditionally independent. 6. Connection with the fluid generator Q and matrix 9 In this section we demonstrate how the fundamental matrices of the SFM and the QBD are related. Assume that the parameters ϑi (1x), i ∈ S , are of the form (5). Consider the process {(X (t ), ϕ(t )), t ≥ 0}. Let θ (x) = inf{t > 0 : X (t ) = x}. For all i ∈ S1 and j ∈ S2 , let
Ψi,j = P (θ (0) < ∞, ϕ(θ (0)) = j|X (0) = 0, ϕ(0) = i)
(8)
and define the matrix 9 = [Ψi,j ]. The matrix 9 is probably the key performance measure in the theory of SFMs [5,6,9] and is the minimum nonnegative solution of the Riccati equation [6] Q12 + Q11 X + XQ22 + XQ21 X = 0,
(9)
where 1 −1 Q11 = C− T01 ], 1 [T11 − T10 (T00 ) 1 −1 Q22 = C− T02 ], 2 [T22 − T20 (T00 ) 1 −1 Q12 = C− T02 ], 1 [T12 − T10 (T00 ) 1 −1 Q21 = C− T01 ], 2 [T21 − T20 (T00 )
and
Q=
Q11 Q21
Q12 Q22
is the key fluid generator in the theory of the SFMs [6]. Further, consider the process {(X1x (t ), ϕ1x (t )), t ≥ 0}. Let τ (k) = inf{t > 0 : X1x (t ) = k}. For all i, j ∈ S , let G(1x)i,j = P (τ (k − 1) < ∞, ϕ1x (τ (k − 1)) = j|X1x (0) = k, ϕ1x (t ) = i),
(10)
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
585
and define the matrix G(1x) = [G(1x)i,j ], the key G-matrix in the theory of QBDs [15]. We partition G(1x) according to S1 ∪ S2 ∪ S0 as follows:
G(1x) =
0 0 0
G(1x)12 G(1x)22 G(1x)02
0 0 . 0
Below, we establish a result that connects the matrices G(1x) and 9. Theorem 7.
9
0 lim G(1x) = 0 1x→0+ 0
I
(−T00 )−1 (T01 9 + T02 )
0 0 . 0
Proof. For notational convenience, let A(−1) = A(−1) (1x), A(0) = A(0) (1x) and A(+1) = A(+1) (1x). Partition A(−1) , A(0) and A(+1) according to S1 ∪ S2 ∪ S0 as
0 A(−1) = 0 0
0
(−1)
A22 0
0 0 , 0
(0)
A11 0) A(0) = A(21 (0) A01
(0)
(0)
A12 (0) A22 (0) A02
A10 (0) A20 , (0) A00
(+1)
A11 A(+1) = 0 0
0 0 0
0 0 , 0
and define the matrices (0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
(0)
˜ 11 = A11 − A10 (A00 )−1 A01 , A ˜ 22 = A22 − A20 (A00 )−1 A02 , A ˜ 12 = A12 − A10 (A00 )−1 A02 , A ˜ 21 = A21 − A20 (A00 )−1 A01 . A
(11)
Matrix A(0) is clearly an infinitesimal generator of an irreducible non-conservative CTMC, and so by arguments similar to (0) (0) ˜ (110) and A˜ (220) , are also such generators and therefore are invertible. those used in [6], matrices A00 and A22 , as well as matrices A By [15], G(1x) is the minimum nonnegative solution of the equation 0 = A(−1) + A(0) G(1x) + A(+1) G(1x)2 , which gives the set of equations (0)
(0)
(0)
(+1)
(0)
(0)
(0)
(−1)
(0)
(0)
(0)
0 = A11 G(1x)12 + A12 G(1x)22 + A10 G(1x)02 + A11 G(1x)12 G(1x)22 , 0 = A21 G(1x)12 + A22 G(1x)22 + A20 G(1x)02 + A22 ,
(12) (13)
0 = A01 G(1x)12 + A02 G(1x)22 + A00 G(1x)02 .
(14)
By (14), we have (0)
(0)
(0)
(0)
(0)
(−1)
G(1x)02 = (−A00 )−1 (A01 G(1x)12 + A02 G(1x)22 ),
(15)
and so by (13),
˜ 22 )−1 (A˜ 21 G(1x)12 + A22 ). G(1x)22 = (−A
(16)
Note that (0)
(0)
A11 = T11 − V11 ,
A12 = T12 ,
(0)
A10 = T10 ,
(+1)
= V11 ,
(−1)
= V22 ,
A11
where V11 = diag(ϑi (1x))i∈S1 , (0)
(0)
A22 = T22 − V22 ,
A21 = T21 ,
(0)
A20 = T20 ,
where V22 = diag(ϑi (1x))i∈S2 , and (0)
A00 = T00 ,
(0)
A01 = T01 ,
Consequently, since
1xV11 = C1 ,
1xV22 = C2 ,
(0)
A02 = T02 .
A22
586
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
1 by left-multiplying both sides of Eq. (12) by C− 1 and using (15) and (16), it follows that 1 −1 −1 −1 0 = C− T01 G(1x)12 1 (T11 − V11 )G(1x)12 + C1 T12 G(1x)22 + C1 T10 (−T00 )
1) 1 1 −1 ˜ (0) −1 A˜ (210) G(1x)12 + A(− + C− T02 G(1x)22 + C− , 1 T10 (−T00 ) 1 V11 G(1x)12 (−A22 ) 22 and so
0 = Q12 G(1x)22 + Q11 G(1x)12 + G(1x)12
1) (−A˜ (220) )−1 A(− −I 22 1x
+ G(1x)12
(−A˜ (220) )−1 A˜ (210) 1x
G(1x)12 .
(17)
At this point we could observe that lim1x→0+ G(1x)22 = I and then note that lim1x→0+ G(1x)12 satisfies Eq. (9). However, we have no way of knowing which solution it might be to Eq. (9) (beyond being nonnegative) and so we cannot connect it to 9. Instead, we shall connect G(1x)12 directly to the minimal nonnegative solution to a perturbed form of Eq. (9) and hence show convergence to 9. By (16), we can write (17) in the equivalent form
˜ (1x)12 + Q˜ (1x)11 G(1x)12 + G(1x)12 Q˜ (1x)22 + G(1x)12 Q˜ (1x)21 G(1x)12 , 0=Q
(18)
where
˜ (1x)12 = Q12 I + Q˜ (1x)22 1x , Q ˜ (1x)11 = Q11 + Q12 Q˜ (1x)21 1x, Q 1) (−A˜ (220) )−1 A(− −I 22 , 1x (−A˜ (220) )−1 A˜ (210) = . 1x
˜ (1x)22 = Q ˜ (1x)21 Q
We now take the limit as 1x → 0+ of these coefficient matrices in (18). We have
˜ (1x)21 = lim Q
1x→0+
(0)
˜ 22 )−1 C2 Q21 lim (−1xA
1x→0+
= lim (C2 − 1xC2 Q22 )−1 C2 Q21 1x→0+
= Q21 and
1) (−1xA˜ (220) )−1 (A(− + A˜ (220) ) 22 1x→0+ = lim (C2 − 1xC2 Q22 )−1 (C2 Q22 )
˜ (1x)22 = lim Q
1x→0+
lim
1x→0+
= Q22 , and so
˜ (1x)12 = Q12 , lim Q
1x→0+
˜ (1x)11 = Q11 . lim Q
1x→0+
Clearly, the convergence of the above matrices is uniform, as they are finite-dimensional. ˜ (1x), of a SFM with fluid generator Now, consider a 9-matrix, denoted as 9
˜ (1x) = Q
˜ (1x)11 Q ˜ (1x)21 Q
˜ (1x)12 Q ˜ (1x)22 . Q
˜ (1x) is a generator, simply verify that it has nonnegative off-diagonals and all row sums equal to zero. Since, To see that Q ˜ (1x) is also the minimal nonnegative solution by above, G(1x)12 is the minimal nonnegative solution to (18), and, by [6], 9 ˜ ˜ to (18), it follows that 9(1x) = G(1x)12 . Further, since Q(1x) converges uniformly to Q, then by an argument analogous ˜ (1x) = 9. Specifically, we can derive integral expressions for 9 ˜ (i, 1x) to that in [22, Lemma 1] it follows that lim1x→0+ 9 ˜ (1x) in which exactly i peaks occur, and then let Sm1x = m ˜ corresponding to paths contributing to 9 9 ( i , 1 x ) and argue i=1
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
587
that 1x lim lim Sm ∆→0+ m→∞
˜ (1x) = lim 9 ∆→0+
= lim lim Sm1x m→∞ ∆→0+
= 9.
Consequently, lim1x→0+ G(1x)12 = 9. Also, by (16), lim G(1x)22 =
1x→0+
(0)
(0)
(−1)
˜ 22 )−1 A˜ 21 G(1x)12 + A22 lim (−A
1x→0+
1) = lim (−1xA˜ (220) )−1 1xA˜ (210) G(1x)12 + 1xA(− 22 1x→0+ −1
= C2 (0 + C2 ) = I, and so by (15), lim G(1x)02 = (−T00 )−1 (T01 9 + T02 ).
1x→0+
˜ (1x) converges to Q, it is not the most physically intuitive such perturbation. In fact, It is worth noting that although Q ˜ we regard Q(1x) as a construction to prove the above theorem, not a useful or insightful matrix in its own right. In the following result, we present a different matrix that converges to Q and has a much more appealing physical interpretation that strengthens the connection between the QBD and the SFM. Corollary 2.
Q = lim
1x→0+
1) (−A˜ (110) )−1 A(+ −I 11 ( 0 ) −1 ˜ (0) ˜ 1x (−A22 ) A21
1
(−A˜ (110) )−1 A˜ (120) (−1) (−A˜ (220) )−1 A22 −I
.
(19)
Proof. The expressions for the block matrices in the bottom row of (19) have been derived in the Proof of Theorem 7. The expressions for the block matrices in the top row follow by symmetry. Note that by (11), 1) • (−A˜ (110) )−1 A(+ is the probability matrix of the first transition from level k and the set of phases S1 ∪ S0 (with phases in 11 S2 taboo) to level (k + 1) in some phase in S1 , given the process starts from level k in some phase in S1 , • (−A˜ (110) )−1 A˜ (120) is the probability matrix of the first transition from level k and the set of phases S1 ∪ S0 (with phases in S2
taboo) to level k in some phase in S2 , given the process starts from level k in some phase in S1 ,
• (−A˜ (220) )−1 A˜ (210) is the probability matrix of the first transition from level k and the set of phases S2 ∪ S0 (with phases in S1 taboo) to level k in some phase in S1 , given the process starts from level k in some phase in S2 , and
1) • (−A˜ (220) )−1 A(− is the probability matrix of the first transition from level k and the set of phases S2 ∪ S0 (with phases in 22 S1 taboo) to level (k − 1) in some phase in S2 , given the process starts from level k in some phase in S2 .
Therefore, the right-hand side of (19) can be interpreted as the corresponding rate matrix, with respect to level, of these transitions. We note that the extension to the bounded SFMs follows by using 9, 4 [6] (a concept symmetrical to 9), fluid generator Q, and level-forward and level-reversed arguments. Here, we focused on establishing the results for matrix 9 and the fluid generator Q mainly, since these two are the key matrices in the theory of SFMs from which the other quantities can be derived. 7. Uniformization of a bounded SFM Now, let {(X (t ), ϕ(t )), t ≥ 0} be a Markovian SFM with phase variable ϕ(t ), level variable X (t ), a lower boundary X (t ) ≥ 0, upper boundary X (t ) ≤ B, and real rates ci for all i ∈ S , such that
• when 0 < X (t ) < B and ϕ(t ) = i, then the rate at which the level is changing is ci ; • when X (t ) = 0 and ϕ(t ) = i, then the rate at which the level is changing is max{ci , 0}; and • when X (t ) = B and ϕ(t ) = i, then the rate at which the level is changing is min{ci , 0}.
588
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
That is, when the SFM hits level zero, which must occur in some phase j with cj < 0, it will stay on level zero until a transition to some phase k with ck > 0 occurs. Similarly, when the SFM hits level B, which must occur in some phase j with cj > 0, it will stay on level B until a transition to some phase k with ck < 0 occurs. We choose some large positive integer K , let 1x = B/K and, as before, construct the spatially-coherent continuous-time QBD with state space G = {(k, i) : k = 0, 1, 2, . . . , K , i ∈ S } and generator T(1x) given by
T(1x) =
B(0) (1x) A(−1) (1x)
A(+1) (1x) A(0) (1x)
. ···
. ···
..
..
··· A(+1) (1x) .. .
A(−1) (1x)
··· B(K ) (1x)
,
where B(0) (1x) = [B(0) (1x)i,j ], B(K ) (1x) = [B(K ) (1x)i,j ], and
−γi (1x) B (1x)i,j = −λi
j = i, ci > 0 j = i, ci ≤ 0 j ̸= i
(0)
Ti,j
−γi (1x) (K ) B (1x)i,j = −λi
j = i, ci < 0 j = i, ci ≥ 0 j ̸= i,
Ti,j
and with the remaining block matrices as defined earlier in Section 3. We note that when the bounded QBD level is between levels 0 and K , on level 0 in phase i with ci ≥ 0, or on level K in phase i with ci ≤ 0, the process behaves just like the unbounded QBD of the earlier sections. When the QBD is on level 0 in phase i with ci < 0, or on level K in phase i with ci > 0, the process behaves just like the unbounded QBD in phase i with ci = 0. Under this interpretation it is clear that Theorem 1, Lemma 1 and Theorem 2 must still apply. Further, the weak convergence arguments in Section 5 can all be reworked in this case. Alternatively, we can apply the Continuous Mapping Theorem [21, Corollary 3.1.9, p. 103] on noting that the transformation from the unbounded SFM to the bounded SFM (and similarly for the spatially-coherent QBD approximations) is continuous. 8. Numerical example Consider an SFM {(X (t ), ϕ(t )), t ≥ 0}, X (t ) ∈ [0, +∞), with S = {1, 2, 3, 4, 5}, rates vector c = [ci ] given by
c= 1
4
−2
0 ,
−1
and generator
−8 1 T= 1 1 1
1 −6 1 1 1
4 2 −5 2 1
2 2 2 −5 1
1 1 1 . 1 −4
Calculations reveal that the drift for this process is equal to µ = −0.1524, and so the process is positive recurrent. Also, we have
−7.7500 0.3125 Q= 0.6250 1.2500
1.2500 −1.4375 0.6250 1.2500
4.2500 0.5625 −2.3750 2.2500
2.2500 0.5625 1.1250 −4.7500
(20)
and
9=
0.6854 0.6850
0.3146 . 0.3150
(21)
Consider {(X1x (t ), ϕ1x (t )), t ≥ 0} with ϑi (1x) given by (5) and 1x = 1/n. Denote by 9(n) the approximation of 9 based on Theorem 7, calculated using MATLAB. We illustrate the convergence of 9(n) to 9 in Figs. 1–2 and note, because both 9 and 9(n) are conservative generators, that the errors must necessarily sum to zero across each row. The convergence is such that a log–log plot is almost linear with asymptotic slope of −1, suggesting that the error has the form of O(n−1 ). We also present simulations to enable comparison between the fluid model (X (t ), ϕ(t )) and the approximating scaled QBD (Yn (t ), ϕ(t )), for various values of n. In order to simulate these processes in a fully-coupled manner, the phase transitions (i.e. the process ϕ(t )) were sampled once and kept the same for all the processes. Of course, that is the only source of randomness in the SFM and so the process (X (t ), ϕ(t )) was calculated directly from this. For the QBD approximations,
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
589
Fig. 1. 91j − 9(n)1j with scale 10−4 (y-axis).
Fig. 2. 92j − 9(n)2j with scale 10−5 (y-axis).
once the residence time is known for state i, say τ , the level change in the QBD (for each value of n) during this residence time is then given by a Poisson distribution with parameter ϑi (1/n)τ . We, therefore, sampled from the necessary Poisson distributions and built the QBD level accordingly. This is the most efficient way to simulate the QBD, and also enables the required coupling. In Fig. 3 we show the first 30 phases in the given sample path. On the graph we give the value of the fluid X (t ) and for two (scaled) QBD approximations Yn (t ) with n = 102 and 104 . As can clearly be seen, with n = 102 the approximation is very good, and once we make n = 104 it is essentially indistinguishable from the fluid. To explore this more fully, we present Figs. 4 and 5 that show the difference between the processes X (t ) and Yn (t ) for values of n ranging through 10, 102 , 103 , 104 , 105 and 106 . It is clear that each time the value of n is increased by an order of magnitude, there is a significant reduction in the magnitude of this difference. 9. Conclusion In this paper we have developed a natural spatial-uniformization of the Stochastic Fluid Model (SFM) into a QBD. The approach taken follows the standard uniformization argument, with a further discretization of the continuous fluid space. By reinstating continuous time, we developed a continuous-time QBD whose level naturally approximates the level of the SFM. This differs from the two previous QBDs derived from an SFM, where the level had no direct relationship to the SFM [9], or which required a further calculation involving the workload remaining in the queue of the QBD [2], and where the level represented, instead, the number of customers in the queue. We were able to show formal weak convergence of the appropriately scaled approximating QBDs to the SFM and that the generator of the fluid Q and the key matrix 9 are the limits of the natural QBD quantities. More importantly, this leaves us with a discrete-space approximation to the SFM, which we can use in the development of numerical schemes for the analysis of the stochastic fluid–fluid model, constructed in [18].
590
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Fig. 3. The fluid level X (t ) and two approximating scaled QBD levels Yn (t ) for n = 102 and 104 .
Fig. 4. X (t ) − Yn (t ) for n = 10, 102 , 103 .
Fig. 5. X (t ) − Yn (t ) for n = 104 , 105 , 106 .
We emphasize that the aim of the paper was to construct a meaningful uniformization in which the level of the SFM is directly linked with the level of the QBD so that we can use this theoretical connection between the two models and derive
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
591
more information about the SFMs. In particular, the spatially-coherent uniformization derived in the paper opens the door to the treatment of the multi-dimensional SFMs, since discretizing one of the SFMs in the two-dimensional model results in a one-dimensional SFM which we can treat using the existing methods. This work, which requires rigorous construction and validation, is under preparation and will be reported in our upcoming paper. Another important advantage achieved, since the theory of the QBDs is well understood, is that we can use the results and insights for the QBDs in order to derive further properties of the SFM, which were not considered previously. For example, our work on the quasi-stationary distribution of the SFMs is under preparation, and we will report on this in due course. We believe that more such opportunities will arise from this construction. Acknowledgment The authors would like to thank the Australian Research Council for funding this research through Discovery Project DP110101663. Appendix. Statement of necessary theorems from Ethier and Kurtz Theorem (3.2.2, p. 104 [21]). Let (S , d) be complete and separable and let M ⊂ P (S ) (where P (S ) is the family of Borel probability measures on S). Then the following are equivalent: (a) M is tight. (b) For each ϵ > 0, there exists a compact set K ⊂ S such that inf P (K ϵ ) ≥ 1 − ϵ,
P ∈M
where K ϵ = x ∈ S : infy∈K d(x, y) < ϵ . (c) M is relatively compact.
Proposition (3.2.4, p. 107 [21]). ∞ ∞Let −(kSk , dk ), k = 1, 2, . . . , be metric spaces, and define the metric space (S , d) by letting S = S and d ( x , y ) = (dk (xk , yk ) ∧ 1) for all x, y ∈ S. Let {Pα } ⊂ P (S ) (where α ranges over some index k k=1 k=1 2 set), and for k = 1, 2, . . . and each α , define Pαk ∈ P (S ) to be the kth marginal distribution of Pα (i.e., Pαk = Pα πk−1 , where the projection πk : S −→ Sk is given by πk (x) = xk ). Then {Pα } is tight if and only if {Pαk } is tight for k = 1, 2, . . . . Theorem (3.7.2, p. 128 [21]). Let (E , r ) be complete and separable, and let {Xα } be a family of processes with sample paths in DE [0, ∞). Then {Xα } is relatively compact if and only if the following two conditions hold: (a) For every η > 0 and rational t ≥ 0, there exists a compact set Γη,t ⊂ E such that η
inf P {Xα (t ) ∈ Γη,t } ≥ 1 − η, α η
where Γη,t = x ∈ S : infy∈Γη,t d(x, y) < η . (b) For every η > 0 and T > 0, there exists δ > 0 such that
sup P {w ′ (Xα , δ), T ≥ η} ≤ η. α
Theorem (3.7.8(b), p. 131 [21]). Let E be separable and let Xn , n = 1, 2, . . . , and X be processes with sample paths in DE [0, ∞). (b) If {Xn } is relatively compact and there exists a dense set D ⊂ [0, ∞) such that
(Xn (t1 ), . . . , Xn (tk )) ⇒ (X (t1 ), . . . , X (tk )) holds for every finite set {t1 , . . . , tk } ⊂ D, then Xn ⇒ X . Theorem (3.8.6, p. 137 [21]). Let (E , r ) be complete and separable, and let {Xα } be a family of processes with sample paths in DE [0, ∞). Suppose that condition (a) of Theorem 3.7.2 holds. Then the following are equivalent: (a) {Xα } is relatively compact. (b) For each T > 0, there exist β > 0 and a family {γα (δ) : 0 < δ < 1 all α} of nonnegative random variables satisfying E qβ (Xα (t + u), Xα (t )) |Ftα qβ (Xα (t ), Xα (t − v)) ≤ E γα (δ)|Ftα
for 0 ≤ t ≤ T , 0 ≤ u ≤ δ , and 0 ≤ v ≤ δ ∧ t, where Ftα ≡ Ft α ; in addition, X
lim sup E [γα (δ)] = 0
δ→0 α
and lim sup E qβ (Xα (δ), Xα (0)) = 0.
δ→0 α
592
N.G. Bean, M.M. O’Reilly / Performance Evaluation 70 (2013) 578–592
Remark (3.8.7, p. 138 [21]). (a) If, as will typically be the case, E qβ (Xα (t + u), Xα (t )) |Ftα ≤ E γα (δ)|Ftα
instead of E qβ (Xα (t + u), Xα (t )) |Ftα qβ (Xα (t ), Xα (t − v)) ≤ E γα (δ)|Ftα ,
then E qβ (Xα (δ), Xα (0)) ≤ E [γα (δ)] and we need only verify
lim sup E [γα (δ)] = 0
δ→0 α
in condition (b) of Theorem 3.8.6. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
S. Ahn, V. Ramaswami, Fluid flow models and queues—a connection by stochastic coupling, Stochastic Models 19 (3) (2003) 325–348. S. Ahn, V. Ramaswami, Transient analysis of fluid flow models via stochastic coupling to a queue, Stochastic Models 20 (1) (2004) 71–101. S. Ahn, V. Ramaswami, Efficient algorithms for transient analysis of stochastic fluid flow models, Journal of Applied Probability 42 (2) (2005) 531–549. S. Asmussen, Stationary distributions for fluid flow models with or without Brownian noise, Stochastic Models 11 (1) (1995) 21–49. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Algorithms for return probabilities for stochastic fluid flows, Stochastic Models 21 (1) (2005) 149–184. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Hitting probabilities and hitting times for stochastic fluid flows, Stochastic Processes and their Applications 115 (9) (2005) 1530–1556. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Algorithms for the Laplace–Stieltjes transforms of first return times for stochastic fluid flows, Methodology and Computing in Applied Probability 10 (3) (2008) 381–408. V. Ramaswami, Matrix analytic methods: a tutorial overview with some extensions and new results, in: Matrix-Analytic Methods in Stochastic Models (Flint, MI), in: Lecture Notes in Pure and Appl. Math., vol. 183, Dekker, New York, 1997, pp. 261–296. V. Ramaswami, Matrix analytic methods for stochastic fluid flows, in: Proceedings of the 16th International Teletraffic Congress, Edinburgh, 7–11 June 1999, pp. 1019–1030. D. Anick, D. Mitra, M.M. Sondhi, Stochastic theory of a data-handling system with multiple sources, Bell System Technical Journal 61 (8) (1982) 1871–1894. A. Badescu, L. Breuer, A. Da Silva Soares, G. Latouche, M.-A. Remiche, D. Stanford, Risk processes analyzed as fluid queues, Scandinavian Actuarial Journal 105 (2) (2005) 127–141. N.G. Bean, M.M. O’Reilly, J.E. Sargison, A stochastic fluid flow model of the operation and maintenance of power generation systems, IEEE Transactions on Power Systems 25 (3) (2010) 1361–1374. N.G. Bean, M.M. O’Reilly, Performance measures of a multi-layer Markovian fluid model, Annals of Operations Research 160 (2008) 99–120. G. Latouche, P.G. Taylor, A stochastic fluid model for an ad hoc mobile network, Queueing Systems 63 (1–4) (2009) 109–129. G. Latouche, V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, in: ASA-SIAM Series on Statistics and Applied Probability, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1999. N.G. Bean, G.G. Latouche, G.T. Nguyen, The stochastic fluid–fluid model applied to a two-server queue with assistance, 2013, in preparation. N.G. Bean, M.M. O’Reilly, A stochastic two-dimensional fluid model, Stochastic Models 29 (1) (2013) 31–63. N.G. Bean, M.M. O’Reilly, The stochastic fluid–fluid model: a stochastic fluid model driven by an uncountable-state process, which is a stochastic fluid model itself, 2013, submitted for publication. S.M. Ross, Introduction to Probability Models, Elsevier, New York, 2007. P. Billingsley, Convergence of Probability Measures, in: Wiley Series in Probability and Mathematical Statistics: Tracts on Probability and Statistics, Wiley, 1968. S.N. Ethier, T.G. Kurtz, Markov Processes Characterization and Convergence, in: Wiley Series in Probability and Mathematical Statistics, John Wiley & Sons, 1986. N.G. Bean, M.M. O’Reilly, P.G. Taylor, Hitting probabilities and hitting times for stochastic fluid flows: the bounded model, Probability in the Engineering and Informational Sciences 23 (1) (2009) 121–147.
Nigel G. Bean received the B.Sc. degree in mathematics and the Honours degree in applied mathematics from the University of Adelaide, Adelaide, Australia, in 1988 and 1989, respectively. He received the Ph.D. degree from the University of Cambridge, Cambridge, U.K., in 1993, where he studied under Prof. F. Kelly, FRS. Since then he has been employed at the University of Adelaide and was appointed Chair of Applied Mathematics in 2004. His main research interests are in stochastic modeling, particularly Markov chains, and applied operations research. Prof. Bean is the recipient of the 2001 J.H. Michell Medal, awarded by the Australian and New Zealand Industrial and Applied Mathematics (ANZIAM) division of the Australian Mathematical Society, and the 2003 P.A.P. Moran Medal, awarded by the Australian Academy of Science.
Małgorzata M. O’Reilly received the M.Sc. degree in mathematics and education from Wroclaw University, Wroclaw, Poland, in 1987, and the Ph.D. degree in applied probability from the University of Adelaide, Adelaide, Australia, in 2002. She has been a lecturer in probability and operations research at the University of Tasmania, Hobart, Tasmania, Australia, since 2005. Her research is in the area of stochastic modeling, Markov chains, and operations research.