Vol. 64, No. 2, April 2021, pp. 61–86
NUMERICAL IMPLEMENTATION OF THE AUGMENTED TRUNCATION APPROXIMATION TO SINGLE-SERVER QUEUES
WITH LEVEL-DEPENDENT ARRIVALS AND DISASTERS
Masatoshi Kimura Tetsuya Takine Osaka University
(Received April 3, 2020; Revised September 3, 2020)
Abstract This paper considers the computation of the stationary queue length distribution in single-server queues with level-dependent arrivals and disasters. We assume that service times follow a general distribution and therefore, we consider the stationary queue length distribution via an imbedded Markov chain. Because this imbedded Markov chain has infinitely many states, level dependence, and bidirectional jumps of levels, it is hard to compute the solution of the global balance equation exactly. We thus consider the augmented truncation approximation. In particular, we focus on the computation of the truncated state transition probability matrix of the imbedded Markov chain, assuming that the underlying continuous-time absorbing Markov chain during a service time is not uniformizable. Under some stability conditions, we develop a computational procedure for the truncated transition probability matrix, where the upper bound of errors owing to truncation can be set in advance. We also provide some numerical examples and demonstrate that our procedure works well.
Keywords: Queue, disaster, level dependence, augmented truncation approximation, uniformization technique, error bound
1. Introduction
In queueing systems, a disaster is the event that all customers are removed from the system. For instance, resets in computer systems, which arise from overloads or virus infections, can be modeled by disasters. Such events are also called clearing [2], (total) catastrophe [4, 12], or mass exodus [5] in the literature. The stationary workload [2, 9, 10] and the stationary queue length [6, 7, 13, 18] in single-server queues with disasters and generally distributed service times have been studied in the literature.
Queueing systems whose arrival rates depend on the queue length are termed state-dependent queues or level-state-dependent queues. Such queueing systems naturally arise in practical modeling. In many situations, a long waiting line discourages prospective cus-tomers from joining the queue, so that the net arrival rate depends on the queue length. The level-dependent batch Markovian arrival process (LD-BMAP) is known as an generalization of the Markovian arrival process (MAP), and the stationary distribution in LD-BMAP/G/1 queues has been discussed in [8].
It would be natural that the occurrence of disasters also depends on the state of the system. For instance, a huge backlog may cause a fatal fault or a compulsory reset. This paper considers a single-server queue, where arrivals and disasters occur according to a level-dependent marked Markovian arrival process (LD-MMAP) described in later. We assume that customers are served on a FIFO basis and that their service times are independent and identically distributed (i.i.d.) according to a general distribution H(x) (x ≥ 0) with finite
mean E[H]. To the best of our knowledge, single-server queues with generally distributed service times and level-dependent disasters have not been studied so far.
We define L(t) (t ≥ 0) as the number of customers in the system at time t and S(t) (t ≥ 0) as the state of the auxiliary variable at time t. We call L(t) level at time t and S(t) phase at time t. We also define level k as the subset {k} × M (k = 0, 1, . . .) of the state space {0, 1, . . .} × M of {(L(t), S(t))}t≥0. Let N (s, t] (s < t) denote the number of customers entering the system in time interval (s, t] under the assumption that no service completions occur in (s, t]. We then define Ck, Dk, and Γk(k = 0, 1, . . .) as M×M matrices whose elements are given by
[Ck]i,j = lim ∆t→0 Pr N (t, t + ∆t] = 0, S(t + ∆t) = j | (L(t), S(t)) = (k, i) ∆t , i, j ∈ M, i ̸= j, [Dk]i,j = lim ∆t→0 Pr N (t, t + ∆t] = 1, S(t + ∆t) = j | (L(t), S(t)) = (k, i) ∆t , i, j ∈ M, [Γk]i,j = lim ∆t→0 Pr a disaster occurs in (t, t + ∆t], S(t + ∆t) = j | (L(t), S(t)) = (k, i) ∆t , i, j ∈ M, and [Ck]i,i =− X j∈M j̸=i [Ck]i,j+ X j∈M [Dk]i,j+ [Γk]i,j , i∈ M.
We assume Γ0 = O without loss of generality. We also assume Dk ̸= O (k = 0, 1, . . .) and 0 < [−Ck]i,i < ∞ (k = 0, 1, . . .; i ∈ M). By definition, {(Ck, Dk, Γk); k = 0, 1, . . .} describes the arrival/disaster process. Roughly speaking, when k (k = 0, 1, . . . ) customers stay in the system (i.e., L(t) = k), an arrival (resp. a disaster) occurs along with a transition driven by Dk (resp. Γk). By definition, (Ck + Dk + Γk) (k = 0, 1, . . . ) represents the infinitesimal generator of the underlying phase process {S(t)}t≥0 when L(t) = k. Note that (Ck+ Dk+ Γk)e = 0, k = 0, 1, . . . , (1.1) where e denotes a column vector with an appropriate dimension, whose elements are all equal to one. Note that when a service completes, L(t) decreases by one, while S(t) remains unchanged. Throughout this paper, matrices, row vectors, and column vectors are denoted by capital bold letters, lowercase bold Greek letters, and lowercase bold Roman letters, respectively.
We are interested in the stationary distribution of {(L(t), S(t))}t≥0. We thus assume that for any pair of states (k1, j1) and (k2, j2) (k1, k2 = 0, 1, . . .; j1, j2 ∈ M), there exists
a sample path from (k1, j1) to (k2, j2). Unfortunately, the necessary and sufficient stability
condition for queueing systems with level-dependent arrivals and disasters has not been clarified. In this paper, we assume the following sufficient condition for the stability. Assumption 1.1. {(L(t), S(t))}t≥0 satisfies either of the following two conditions. (i) There exists K1†∈ {0, 1, . . .} and α > 0 such that Dke− kΓke≤ −α (k ≥ K1†), or
(ii) there exists K2† ∈ {0, 1, . . .} and α > 0 such that Dke− µe ≤ −α (k ≥ K2†), where
µ = 1/E[H].
Assumption 1.1 (i) is a drift condition ignoring service completions, while Assumption 1.1 (ii) is a drift condition ignoring disasters. Note that under Assumption 1.1, the ar-rival/disaster process characterized by {(Ck, Dk, Γk); k = 0, 1, . . .} is non-explosive, i.e.,
for an arbitrary finite T > 0 and an arbitrary initial state (L(0), S(0)), the number of arrivals in time interval (0, T ] is finite with probability 1.
Remark 1.1. The stability under Assumption 1.1 (i) can be shown as follows. For a while, we assume that services never complete. In this case, {(L(t), S(t))}t≥0 forms a continuous-time Markov chain and it is readily shown by Pakes’s Lemma [3, Corollary 1.1 of Chap. 5] that {(L(t), S(t))}t≥0 is positive-recurrent under Assumption 1.1 (i). Therefore, in the model with a finite E[H] described above, the mean first passage time from any state to the empty system is finite, so that {(L(t), S(t))}t≥0 is stable. On the other hand, the stability under Assumption 1.1 (ii) can be shown by considering a modified model in which all cus-tomers remain in the system at the occurrence of disasters (i.e., the ordinary LD-MAP/G/1 queueing model without disasters), and applying Theorem 3.3-3 in [8].
Let π = (π0 π1 · · · ) denote the stationary distribution of {(L(t), S(t))}t≥0, where πℓ denotes a 1× M vector whose jth (j ∈ M) element is given by limt→∞Pr (L(t), S(t) = (ℓ, j). Since service times are generally distributed, we consider an imbedded Markov chain associated with {(L(t), S(t))}t≥0. Specifically, the standard approach to the analysis of the stationary distribution π is outlined as follows.
Step 1: Construct an imbedded Markov chain: We choose time instants at which services complete and disasters occur as imbedded Markov points. We define τn (n = 0, 1, . . .) as the nth imbedded Markov point, where τ0 = 0. Let (Ln, Sn) = (L(τn), S(τn)) (n = 0, 1, . . .). It is clear that {(Ln, Sn)}n=0,1,... forms a discrete-time Markov chain whose state transition probability matrix P takes the following form:
P = A0,0 A0,1 A0,2 A0,3 A0,4 A0,5 · · · A1,0 A1,1 A1,2 A1,3 A1,4 A1,5 · · · A2,0 A2,1 A2,2 A2,3 A2,4 A2,5 · · · A3,0 O A3,2 A3,3 A3,4 A3,5 · · · A4,0 O O A4,3 A4,4 A4,5 · · · A5,0 O O O A5,4 A5,5 · · · .. . ... ... ... ... ... . .. , (1.2)
where Ak,0 (k = 1, 2, . . .) and Ak,ℓ (k = 0, 1, . . ., ℓ = k− 1, k, . . .) are M × M matrices. Step 2: Compute the stationary distribution of the imbedded Markov chain: We
de-fine κ = (κ0 κ1 · · · ) as the stationary distribution of the imbedded Markov chain
{(Ln, Sn)}n=0,1,..., where κℓ denotes a 1× M vector whose jth (j ∈ M) element is given by limn→∞Pr((Ln, Sn) = (ℓ, j)). By definition, κ satisfies κ = κP and κe = 1.
Step 3: Compute the stationary distribution π = (π0 π1 · · · ) of {(L(t), S(t))}t≥0: π is given in terms of κk (k = 0, 1, . . .) as follows.
π0 = κ0(−C0)−1 E[T ] , πℓ = κ0(−C0)−1D0W1,ℓ E[T ] + ℓ X k=1 κkWk,ℓ E[T ] , ℓ = 1, 2, . . . , (1.3) where E[T ] denotes the average of τn+1 − τn (n = 0, 1, . . . ) and Wk,ℓ (k = 1, 2, . . . , ℓ = k, k + 1, . . . ) denotes an M× M matrix whose (i, j)th (i, j ∈ M) element represents the mean total sojourn time in state (ℓ, j) in (τn, τn+1], starting from (Ln, Sn) = (k, i). The purpose of this paper is to establish a method for numerical implementation of the above three steps in the LD-MMAP/G/1 queue with disasters. In what follows, we identify the problems to be resolved in each step.
In Step 1, we have to compute Ak,ℓ’s (k = 0, 1, . . ., ℓ = 0, k− 1, k, . . .) in P . As we will see in (2.1), Ak,ℓ’s are given in terms of M × M matrices BA,k,ℓ’s in BA:
BA= I O O O O · · ·
BA,1,0 BA,1,1 BA,1,2 BA,1,3 BA,1,4 · · ·
BA,2,0 O BA,2,2 BA,2,3 BA,3,4 · · ·
BA,3,0 O O BA,3,3 BA,3,4 · · ·
BA,4,0 O O O BA,4,4 · · · .. . ... ... ... ... . .. = Z ∞ 0 exp[U t]dH(t), (1.4)
where U denotes the infinitesimal generator of an absorbing Markov chain that represents the arrival process until a disaster occurs.
U = O O O O O · · · Γ1 C1 D1 O O · · · Γ2 O C2 D2 O · · · Γ3 O O C3 D3 · · · Γ4 O O O C4 · · · .. . ... ... ... ... . .. . (1.5)
In what follows, we use the following convention, as with BA in (1.4). For any matrix X
composed of M × M block matrices, let Xk,ℓ (k, ℓ≥ 0) denote the (k, ℓ)th block matrix of X. For example, for U in (1.5), Uk,k = Ck (k = 1, 2, . . .).
Remark 1.2. Note that exp[U t] (t≥ 0) can be regarded as a transition probability matrix, i.e., exp[U t] ≥ O and exp[Ut]e = e. Therefore, all elements in BA are finite because
R∞
0 exp[U t]dH(t)e =
R∞
0 dH(t)e = e.
The stable and accurate numerical implementation of block matrices BA,k,ℓ in BA is
crucial and it is one of the main topics of this paper. We define θ as θ = sup k=0,1,... max i∈M[−Ck]i,i .
If θ is finite, the uniformization technique is directly applicable and BA,k,ℓ (k = 0, 1, . . .,
ℓ = 0, k − 1, k, . . .) is given by the weighted sum of infinitely many M × M non-negative matrices, where the weights are given in terms of θ and the service time distribution H(x). The advantages of the uniformization technique are that (i) the procedure is numerically stable because it does not involve subtractions and that (ii) for a given error bound ϵ, we can set an appropriate truncation point n := n∗(ϵ) of the infinite sum representing BA,k,ℓ.
Note that θ can be infinite and in this case, there are no standard methods for computing block matrices BA,k,ℓ in BA in (1.4), which guarantee error bounds. In this paper, we
consider computing an approximation to BA,k,ℓ by using only the (m + 1)M × (m + 1)M
north-west corner submatrix U (m) of U because exp[U (m)t] is always uniformizable since [−Ck]i,i <∞ (i ∈ M) for all k = 1, 2, . . ..
The truncation of U at level m does not affect the accuracy of computed BA,k,ℓ (ℓ =
k, k + 1, . . . , m) because it can be represented completely in terms of U (m). Contrarily, the accuracy of computed BA,k,0is affected by the truncation of U by the following reason. Note
that BA,k,0 includes the probability that a disaster occurs before a service completion given
that the service starts with Ln= k ≤ m. We then consider the level L(τn+1−) immediately before the occurrence of a disaster. Because Pr(L(τn+1−) = ℓ | L(τn) = k) > 0 for all
ℓ > m, the exact expression of BA,k,0 is given in terms of the whole U . Therefore, if we
compute BA,k,0 using the truncated U (m) of U at level m, we will obtain an approximation
BA,k,0(m) to BA,k,0, which ignores all sample paths with L(τn+1−) > m. Consequently, in terms of an approximation Ak,0(m) to Ak,0, we have for k = 0, 1, . . . , m and i, j ∈ M,
[Ak,0(m)]i,j = Pr (Ln+1, Sn+1) = (0, j), L(τn+1−) ≤ m | (Ln, Sn) = (k, i)
.
The above discussion makes it clear that our approach has two sources of errors: one is the truncation of U at level m and the other is the truncation of the weighted infinite sum obtained by uniformization. While the latter is controllable, the error control of the former is not straightforward because it requires information beyond level m.
To control errors inherent in our approach, we adopt a cross-layer design of the compu-tational procedure. Specifically, we assume that the augmented truncation approximation (ATA) is employed in Step 2, as we will explain below, where the truncation level K in the ATA is assumed to be given in advance. We then utilize this information in Step 1, i.e., we consider a computational procedure for Ak,0 only for k = 0, 1, . . . , K, which satisfies a predefined error bound under Assumption 1.1.
In Step 2, we compute the stationary distribution κ of the imbedded Markov chain {(Ln, Sn)}n=0,1,.... Owing to the level dependence of arrivals and disasters, however, it is hard to obtain the analytical expression of the stationary distribution κ. We thus adopt the ATA, which attempts to compute an approximation κapprox to the stationary distribution κ as follows. We first choose an appropriate level K and partition the stationary distribution κ as κ = (κ(1)(K) κ(2)(K)), where
κ(1)(K) = (κ0 κ1 . . . κK), κ(2)(K) = (κK+1 κK+2 . . .).
In the ATA, we approximate κ(2)(K) to be 0 and attempt to compute κ(1)(K) based on the north-west corner submatrix P (K) of the transition probability matrix P in (1.2), where
P (K) = A0,0 A0,1 A0,2 A0,3 · · · A0,K−1 A0,K A1,0 A1,1 A1,2 A1,3 · · · A1,K−1 A1,K A2,0 A2,1 A2,2 A2,3 · · · A2,K−1 A2,K A3,0 O A3,2 A3,3 · · · A3,K−1 A3,K .. . ... ... ... . .. ... ... AK−1,0 O O O · · · AK−1,K−1 AK−1,K AK,0 O O O · · · AK,K−1 AK,K .
Specifically, we compute an approximation κapprox(K) to κ(1)(K) by solving
κapprox(K) = κapprox(K)[P (K) + PA(K)], κapprox(K)e = 1, (1.6)
where PA(K) denotes an augmentation matrix such that PA(K)≥ O and [P (K)+PA(K)]e
= e. We then adopt κapprox = (κapprox(K) 0) as an approximation to κ, which is known to converge to κ as K goes to infinity under some technical conditions [14, 16]. In the literature, upper bounds of the error∥κapprox−κ∥
1 for κapprox are also studied under various
drift conditions [14, 15].
Let κ(K) = (κ0(K) κ1(K) · · · κK(K)) denote the conditional stationary distribution given that the level is not greater than K, where κℓ(K) (ℓ = 0, 1, . . . , K) denotes a 1× M vector whose jth (j ∈ M) element is given by limn→∞Pr((Ln, Sn) = (ℓ, j) | Ln ≤ K).
Note here that κapprox(K) in (1.6) can be regarded as an approximation to the conditional
stationary distribution κ(K). If numerical errors in elements of P (K) are negligible, the error in an approximation κapprox = (κapprox(K) 0) to κ can be evaluated by the tail
probability and the difference between the ATA solution κapprox(K) and κ(K) [11, 19].
2ζ(K)≤ ∥κ − κapprox∥1 ≤ 2ζ(K) + ϵ(K), (1.7)
where ∥ · ∥1 stands for the total variation norm or ℓ1-norm and
ζ(K) = lim
n→∞Pr(Ln > K), ϵ(K) =∥κ
approx(K)− κ(K)∥ 1.
The tail probability ζ(K) decreases monotonically as K increases. Unfortunately, how-ever, it is hard to analyze the tail probability in our model due to the level-dependence. We thus evaluate it numerically after obtaining the ATA solution. On the other hand, the error ϵ(K) in the ATA solution κapprox(K) depends on the selection of the augmentation
matrix PA(K). It would be ideal (i.e., ϵ(K) = 0) if we could select it in such a way that
P (K) + PA(K) equals the state transition matrix of the censored Markov chain obtained
by observing{(Ln, Sn)}n=0,1,... only when Ln ≤ K [19]. In this paper, we discuss reasonable selections of the augmentation matrix PA(K), based on the results in [11, Section 2.4].
In Step 3, we compute the approximation πapprox = (πapprox
0 π
approx
1 · · · ) to the stationary
distribution π using κapprox = (κapprox(K) 0). It follows from (1.3) that πapprox
ℓ (ℓ = 0, 1, . . .) is given in terms of κapprox(K):
π0approx= κ approx 0 (K)(−C0)−1 c , (1.8) πℓapprox= κ approx 0 (K)(−C0)−1D0W1,ℓ c + min(ℓ,K)X k=1 κapproxk (K)Wk,ℓ c , ℓ = 1, 2, . . . , (1.9) where c denotes the normalizing constant so that πapproxe = 1. Note here that W
k,ℓ’s in (1.9) are given by block matrices in W defined as
W = I O O O O · · · W1,0 W1,1 W1,2 W1,3 W1,4 · · · W2,0 O W2,2 W2,3 W2,4 · · · W3,0 O O W3,3 W3,4 · · · W4,0 O O O W4,4 · · · .. . ... ... ... ... . .. = Z ∞ 0 exp[U t] 1− H(t)dt.
It then follows that
W = E[H]BW, where BW = Z ∞ 0 exp[U t]1− H(t) E[H] dt. (1.10)
Because (1− H(t))/E[H] is the probability density function of elapsed service times, BW
can be computed in the same way as BA in (1.4).
If the numerical errors in Wk,ℓ(k = 1, 2, . . . , K, ℓ = k, k + 1, . . .) are negligible, the error in πapprox can be evaluated as follows.
∥π − πapprox∥
where cp is a finite coefficient. (1.11) shows that we can control the error in πapprox through the error in κapprox.
In summary, to guarantee a sufficient accuracy of πapprox, we have to pay attention to
the followings.
1. K should be set appropriately so as to make the tail probability ζ(K) negligible, 2. the numerical errors in Ak,ℓ’s and Wk,ℓ’s should be sufficiently small, and
3. the error ϵ(K) in κapprox(K) due to the selection of PA(K) should be small as much as
possible.
In this paper, we mainly consider the last two points. Through numerical examples, we also discuss the upper bound of the tail probability and an appropriate truncation level K in the ATA.
The rest of this paper is organized as follows. In Section 2, we consider the numerical implementation of Ak,ℓ’s and Wk,ℓ’s. In Section 3, we briefly discuss the selection of the augmentation matrix and the error bound in the ATA solution. In section 4, we discuss the error in the approximation to the stationary distribution. In Section 5, we show some numerical examples. Finally, we conclude this paper in Section 6.
2. Computation of Ak,ℓ’s and Wk,ℓ’s
As stated in the preceding section, we will adopt the ATA in Step 2. We thus consider the computation of the north-west corner submatrix P (K) of P for a given K. By definition, block matrices Ak,ℓ’s (k, ℓ = 0, 1, . . . , K) in P (K) are given in terms of M × M block matrices in BA defined by (1.4). Ak,ℓ= (−C0)−1D0A1,ℓ, k = 0, ℓ = 0, 1, . . . , K, BA,1,0+ BA,1,1, k = 1, ℓ = 0, BA,1,ℓ+1, k = 1, ℓ = 1, 2, . . . , K, BA,k,0, k = 2, 3, . . . , K, ℓ = 0, BA,k,ℓ+1, k = 2, 3, . . . , K, ℓ = k− 1, k, . . . , K, O, otherwise. (2.1)
Therefore, we consider the computation of BA,k,ℓ’s below. We also discuss the computation
of Wk,ℓ in this section because W = E[H]BW and block matrices in BW can be computed
in the same way as those in BA.
2.1. The uniformizable case
We first consider the case that exp[U t] is uniformizable, i.e., there exists a finite θ such that [−Ck]i,i ≤ θ (k = 1, 2, . . . ; i ∈ M). In this case, by applying the uniformization technique to BA in (1.4) and BW in (1.10), we can rewrite them to be
BA= Z ∞ 0 e−θtexp[(I + θ−1U )θt]dH(t) = ∞ X n=0 γnVn, BW = Z ∞ 0 e−θtexp[(I + θ−1U )θt]1− H(t) E[H] dt = ∞ X n=0 ηnVn,
where γn and ηn (n = 0, 1, . . .) denote probability functions given by
γn= Z ∞ 0 e−θt(θt) n n! dH(t), ηn= Z ∞ 0 e−θt(θt) n n! · 1− H(t) E[H] dt,
and V = I + θ−1U ≥ O. Note here that V0 = I and V e = e. Block matrices V(n)
k,ℓ in Vn (n = 2, 3, . . .) can be computed recursively by Vk,ℓ(1) = Vk,ℓ and for n = 2, 3, . . .,
Vk,ℓ(n)= Vk,0(n−1)+ k+nX−1 r=k Vk,r(n−1)θ−1Γr, ℓ = 0, Vk,k(n−1)(I + θ−1Ck), ℓ = k, Vk,ℓ(n−1−1)θ−1Dℓ−1+ V (n−1) k,ℓ θ−1(I + θ−1Cℓ), ℓ = k + 1, k + 2, . . . , K + 1, O, otherwise.
In numerical implementation, we assume that the recursion of Vk,ℓ(n) for BA,k,ℓ stops at
n = nA and that for BW,k,ℓ at n = nW.
Bcomp,nA A = nA X n=0 γnVn, BWcomp,nW = nW X n=0 ηnVn.
We then adopt Bcomp,nA
A as an approximation to BA and Bcomp,nW W as an approximation to
BW, where O≤ BAcomp,nA ≤ BA and O ≤ BWcomp,nW ≤ BW. In what follows, we determine
nA = n∗A and nW = n∗W separately in such a way that for k = 0, 1, . . . , K,
K X ℓ=0 (Ak,ℓ− A comp,n∗A k,ℓ )e≤ ϵAe, ∞ X ℓ=k (Wk,ℓ− W comp,n∗W k,ℓ )e≤ ϵWe, (2.2)
where ϵA and ϵW are sufficiently small positive constants.
It follows from (2.1) that the inequality for Acomp,n∗A
k,ℓ (k = 0, 1, . . . , K) in (2.2) is equiv-alent to K+1X ℓ=0 (BA,k,ℓ− B comp,n∗A A,k,ℓ )e≤ ϵAe, k = 1, 2, . . . , K. (2.3)
Note here that K+1X
ℓ=0
(BA,k,ℓ− BA,k,ℓcomp,nA)e≤
∞ X
ℓ=0
(BA,k,ℓ− BA,k,ℓcomp,nA)e =
∞ X n=nA+1 γn ∞ X ℓ=0 Vk,ℓ(n)e = ∞ X n=nA+1 γne.
Therefore, we set nA= n∗A in such a way that
n∗A X n=0
γn≥ 1 − ϵA,
which guarantees (2.3).
Similarly, we have for Wcomp,nW
k,ℓ , ∞ X ℓ=k (Wk,ℓ− Wk,ℓcomp,nW)e≤ ∞ X ℓ=0 (Wk,ℓ− Wk,ℓcomp,nW)e = E[H] ∞ X n=nW+1 ηn ∞ X ℓ=0 Vk,ℓ(n)e = E[H] ∞ X n=nW+1 ηne.
Therefore, we set nW = n∗W in such a way that n∗W X n=0 ηn≥ 1 − ϵW E[H], which guarantees the inequality for Wcomp,n∗W
k,k+ℓ in (2.2). Note that for each k = 1, 2, . . . , K, the computation of Bcomp,n∗W
W,k,ℓ automatically stops at ℓ = k + n∗W and B
comp,n∗W
W,k,ℓ = O for all
ℓ > k + n∗W.
2.2. The non-uniformizable case
In this subsection, we consider the computation of approximations to Ak,ℓ’s and Wk,ℓ’s under the assumption that exp[U t] is non-uniformizable, i.e.,
sup k=0,1,... max i∈M[−Ck]i,i =∞. (2.4) We define U (m) (m = 1, 2, . . .) as U (m) = O O O O · · · O O Γ1 C1 D1 O · · · O O Γ2 O C2 D2 · · · O O Γ3 O O C3 · · · O O .. . ... ... ... . .. ... ... Γm−1 O O O · · · Cm−1 Dm−1 Γm O O O · · · O Cm , m = 1, 2, . . . .
By definition, U (m) (m = 1, 2, . . .) is a defective infinitesimal generator. Let θm (m = 1, 2, . . . ) denote the maximum of the absolute values of the diagonal elements of U (m).
θm = max k=1,2,...,m
i∈M
[−Ck]i,i. m = 1, 2, . . . . In what follows, we first consider Ak,ℓ’s and then consider Wk,ℓ’s.
We assume that U is truncated at level mA in computing Ak,ℓ’s, where mA ≥ K + 1.
We then define BA(mA) as BA(mA) = Z ∞ 0 exp[U (mA)t]dH(t) = ∞ X n=0 γn(mA)Vn(mA), (2.5)
where γn(mA) (n = 0, 1, . . .) denotes a probability function given by
γn(mA) = Z ∞ 0 e−θmAt(θmAt) n n! dH(t),
and V (mA) = I + θ−1mAU (mA). Note that V (mA) is substochastic, i.e., V (mA)e ≤ e.
For any m = 1, 2, . . ., block matrices Vk,ℓ(n)(m) in Vn(m) (n = 1, 2, . . .) can be computed recursively by Vk,ℓ(1)(m) = Vk,ℓ(m) and for n = 2, 3, . . .,
Vk,ℓ(n)(m) = Vk,0(n−1)(m) + min(m,k+nX −1) r=k Vk,r(n−1)(m)θ−1m Γr, ℓ = 0, Vk,k(n−1)(m)(I + θm−1Ck), ℓ = k, Vk,ℓ(n−1−1)(m)θm−1Dℓ−1+ V (n−1) k,ℓ (m)(I + θ−1m Cℓ), ℓ = k + 1, k + 2, . . . , min(m, k + n), O, otherwise. (2.6)
In numerical implementation, we have to truncate the infinite sum in (2.5), i.e., Bcomp,nA A (mA) = nA X n=0 γn(mA)Vn(mA). We adopt Bcomp,nA
A,k,ℓ (mA) as an approximation to BA,k,ℓ. Specifically, for k = 1, 2, . . . , K, an
approximation Bcomp,nA
A,k,ℓ (mA) to BA,k,ℓ is given by
Bcomp,nA A,k,ℓ (mA) = nA X n=1 γn(mA)V (n) k,0 (mA), ℓ = 0, nA X n=ℓ−k γn(mA)V (n) k,ℓ (mA), ℓ = k, k + 1, . . . , min(mA, k + nA), O, otherwise. (2.7)
Furthermore, it follows from (2.1) that an approximation Acomp,nA
k,ℓ (mA) to Ak,ℓ is given by Acomp,nA k,ℓ (mA) = (−C0)−1D0A1,ℓcomp,nA(mA), k = 0, ℓ = 0, 1, . . . , mA− 1, Bcomp,nA
A,1,0 (mA) + BA,1,1comp,nA(mA), k = 1, ℓ = 0,
Bcomp,nA A,1,ℓ+1 (mA), k = 1, ℓ = 1, 2, . . . , mA− 1, Bcomp,nA A,k,0 (mA), k = 2, 3, . . . , K, ℓ = 0, Bcomp,nA A,k,ℓ+1 (mA), k = 2, 3, . . . , K, ℓ = k− 1, k, . . . , mA− 1, O, otherwise. (2.8)
Note here that
Acomp,nA
k,ℓ (mA)≤ Ak,ℓ, k = 0, 1, . . . , K, ℓ = 0, 1, . . . , mA− 1. (2.9) For a given ϵA > 0, we would like to find mA= m∗A and nA= n∗A such that
K X ℓ=0 Ak,ℓ− A comp,n∗A k,ℓ (m∗A) e≤ ϵAe, k = 0, 1, . . . , K. (2.10)
Remark 2.1. It follows from (2.9) that for k = 0, 1, . . . , K, K X ℓ=0 Ak,ℓ− A comp,nA k,ℓ (mA) e≤ mXA−1 ℓ=0 Ak,ℓ− A comp,nA k,ℓ (mA) e ≤X∞ ℓ=0 Ak,ℓe− mXA−1 ℓ=0 Acomp,nA k,ℓ (mA)e = e− mXA−1 ℓ=0 Acomp,nA k,ℓ (mA)e. (2.11)
Therefore, if mA and nA are fixed, we can evaluate the upper bound of the left-hand side of
(2.10) by computing Acomp,nA
k,ℓ (mA) (k = 0, 1, . . . , K, ℓ = 0, k, k + 1, . . . , mA).
This approach has two potential sources of errors: the truncation of U at level mA and
the truncation of the infinite sum in (2.5) at n = nA. Note here that
mA
X ℓ=0
Therefore, if nA≤ mA− K, mXA−1 ℓ=0 Acomp,nA k,ℓ (mA)e = mA X ℓ=0 Bcomp,nA A,k,ℓ e = mA X ℓ=0 nA X n=0 γn(mA)V (n) k,ℓ (mA)e = nA X n=0 γn(mA)e,
for all k = 1, 2, . . . , K. It then follows from (2.11) that K X ℓ=0 Ak,ℓ− Acomp,nk,ℓ A(mA) e≤ 1− nA X n=0 γn(mA) e if nA≤ mA− K. (2.12)
Because the right-hand side of (2.12) takes the minimum value at nA = mA− K, (2.10)
would hold if we could find mA = m∗A such that
1− mX∗A−K n=0 γn(m∗A) e = Pr(Nθm∗ A (H) > m∗A− K) ≤ ϵA, (2.13) where Nθm∗ A
(H) denotes the number of Poisson arrivals with rate θm∗A during a randomly chosen service time H.
Note, however, that m∗Asatisfying (2.13) does not necessary exist in general. The reason is that (i) for a fixed x, Pr(NθmA(H) > x) is an increasing function of θmA and (ii) θmA is
a nondecreasing function of mA, which imply that Pr(NθmA(H) > mA− K) may or may
not decrease as mA increases. Besides, Pr(NθmA(H) > x) depends on the distribution of H. For example, if θmA/mA is equal to zero in the limit mA → ∞, it can be shown that
mA− K > E[NθmA(H)] = θmAE[H] for a sufficient large mA and in this case, the one-side
Chebyshev’s inequality (also called Cantelli’s inequality) [20, pp.198–199] yields
Pr(NθmA(H) > mA− K) ≤
θ2mAVar[H] + θmAE[H]
θ2
mAVar[H] + θmAE[H] + (mA− K − θmAE[H])
2, (2.14)
where the variance Var[H] of service times is assumed to be finite. Furthermore, it can be readily verified that the right-hand side of (2.14) is equal to zero in the limit mA → ∞ if
θmA/mA is equal to zero in the limit mA→ ∞. Another example is that if the service time
distribution belongs to a certain class of long-tailed distributions, we have [1, Eq. (1.1)]
lim x→∞
Pr(NθmA(H) > x) Pr(H > x/θmA)
= 1.
These examples suggest that for the existence of mA= m∗A satisfying (2.13), we need some
additional assumptions on the service time distribution and/or the asymptotic property of the sequence of θmA’s. In what follows, we develop a numerical procedure that always works
under Assumption 1.1 introduced in Section 1.
Theorem 2.1. Consider {(L(t), S(t))}t≥0 whose arrival/disaster process is non-explosive, i.e., P∞ℓ=0Ak,ℓe = e. We then have for k = 0, 1, . . . , K,
K X ℓ=0 Ak,ℓ− Acomp,nk,ℓ A(mA) e≤ max i∈M b + A,(K,i)(mA) e + 1− nA X n=0 γn(mA) e, (2.15)
where b+A,(k,i)(mA) (k = 0, 1, . . . , K; i∈ M) is defined as
b+A,(k,i)(mA) = Pr L(τn+1−) > mA| (Ln, Sn) = (k, i)
Proof. Note first that for an arbitrary ϵ > 0, K X ℓ=0 A1,ℓ− Acomp,n1,ℓ A(mA) e≤ ϵe ⇒ K X ℓ=0
(A0,ℓ− Acomp,n0,ℓ A(mA))e≤ (−C0)−1D0· ϵe
= ϵe. We thus consider (2.15) for k = 1, 2, . . . , K.
It follows from (2.1) and Bcomp,nA
A,k,ℓ (mA)≤ BA,k,ℓ (k = 1, 2, . . . , K, ℓ = 0, 1, . . . , mA) that
for k = 1, 2, . . . , K, K X ℓ=0 Ak,ℓ−Acomp,nk,ℓ A(mA) e = K+1X ℓ=0
BA,k,ℓ−BA,k,ℓcomp,nA(mA)
e≤
∞ X
ℓ=0
BA,k,ℓ−BA,k,ℓcomp,nA(mA)
e,
where Bcomp,nA
A,k,ℓ (mA) = O for ℓ > mA. In what follows, we will show for k = 1, 2, . . . , K,
∞ X
ℓ=0
BA,k,ℓ− BA,k,ℓcomp,nA(mA)
e≤ max i∈M b + A,(K,i)(mA) e + 1− nA X n=0 γn(mA) e. (2.17)
By definition, BA,k,ℓ(mA) = BA,k,ℓ (k = 1, 2, . . . , K, ℓ = 1, 2, . . . , mA), where mA ≥
K + 1. It then follows that for k = 1, 2, . . . , K, ∞ X ℓ=0 BA,k,ℓ− B comp,nA A,k,ℓ (mA) e = BA,k,0− BA,k,0(mA) e + mA X ℓ=0
BA,k,ℓ(mA)− BA,k,ℓcomp,nA(mA)
e + ∞ X ℓ=mA+1 BA,k,ℓe. (2.18) Note here that for k = 0, 1, . . . , K and j ∈ M,
h BA,k,0− BA,k,0(mA) e + ∞ X ℓ=mA+1 BA,k,ℓe i j = b + A,(k,j)(mA)≤ maxi∈M b + A,(K,i)(mA) . (2.19)
On the other hand, we have for k = 1, 2, . . . , K, mA
X ℓ=0
BA,k,ℓ(mA)− BA,k,ℓcomp,nA(mA)
e = mA X ℓ=0 ∞ X n=nA+1 γn(mA)V (n) k,ℓ (mA)e≤ ∞ X n=nA+1 γn(mA)e. (2.20) (2.17) now follows from (2.18), (2.19), and (2.20), which completes the proof.
By definition, maxi∈M b+A,(K,i)(mA)
on the right-hand side of (2.15) monotonically con-verges to zero as mA goes to infinity and for a fixed mA, the second term also converges to
zero monotonically as nAgoes to infinitely. Under Assumption 1.1 (i), we have the following
theorem, where
b+A,K(mA) = b+A,(K,1)(mA) b+A,(K,2)(mA) · · · b+A,(K,M )(mA)
T .
Theorem 2.2. We consider{(L(t), S(t))}t≥0, where (2.4) is assumed to hold. For arbitrary positive integers K and mA (mA≥ K + 1), we have
b+A,K(mA)≤ (−CK)−1DK(−CK+1)−1DK+1· · · (−CmA)
−1D
mAe. (2.21)
Furthermore, for arbitrary K and ϵ > 0, there exists mA= m∗A such that
(−CK)−1DK(−CK+1)−1DK+1· · · (−CmA)
−1Dm∗
Ae≤ ϵe, (2.22)
under Assumption 1.1 (i).
Proof. We first prove (2.21) by a probabilistic argument, even though it can also be shown algebraically. It follows from (2.16) that
b+A,(K,i)(mA) = Pr(L(τn+1−) > mA | (Ln, Sn) = (K, i))
= Pr(L(t) = mA+ 1 for some t∈ (τn, τn+1)| (Ln, Sn) = (K, i)).
Note here that τn+1 = min(dn+1, τn+ Hn), where dn+1 denotes the first occurrence time of a disaster after time τn and Hn denotes the service time starting at time τn. Let ˜N (s, t) (s < t) denote the number of arrivals during an interval (s, t) under the assumption that the service never completes. We then have
b+A,(K,i)(mA) = Pr(L(t) = mA+ 1 for some t ∈ (τn, τn+1)| (Ln, Sn) = (K, i)) ≤ Pr(L(t) = mA+ 1 for some t∈ (τn, dn+1)| (Ln, Sn) = (K, i)) = Pr( ˜N (τn, dn+1) > mA− K | (Ln, Sn) = (K, i))
= [(−CK)−1DK(−CK+1)−1DK+1· · · (−CmA)
−1D mAe]i,
from which (2.21) follows.
Next we consider (2.22) under Assumption 1.1 (i). It follows from (1.1) and Assumption 1.1 (i), i.e., Γke > {1/k} · Dke (k = K1†, K1†+ 1, . . . ), that (−Ck)−1Dke <{k/(k + 1)} · e (k = K1†, K1†+ 1, . . . ). We then have lim mA→∞ (−CK† 1) −1D K1†(−CK1†+1)−1DK1†+1· · · (−CmA) −1D mAe < lim mA→∞ K1† K1†+ 1 · K1†+ 1 K1†+ 2· · · mA mA+ 1 · e = limmA→∞ K1† mA+ 1 · e = 0,
so that for an arbitrary ϵ (0 < ϵ < 1), mA= m∗A satisfying (2.22) exists.
Next we consider Assumption 1.1 (ii). Preliminary to it, we show the following lemma, where for any vector x, we define diag(x) as a diagonal matrix whose ith diagonal element is given by the ith element of x.
Lemma 2.1. Let Ck♮ = Ck+ diag(Γke). (i) For any t ≥ 0, we have
exp[Ckt]≤ exp[Ck♮t], k = 1, 2, . . . . (2.23) (ii) For any M × 1 vector x ≥ 0, we have
(iii) For any M× 1 vector x ≥ 0 and τ ≥ 0, we have Z τ 0 exp[Ck♮t]Dkdt≤ Z τ 0
exp Ck♮− diag(x)t[Dk+ diag(x)(−Ck♮)−1Dk]dt. (2.25) Proof. We first consider (2.23). Since diag(Γke)≥ O, we have Ck ≤ Ck+ diag(Γe) = Ck♮. It then follows that [I + θ−1m
ACk] n≤ [I + θ−1 mAC ♮ k] n and therefore exp[Ckt] = e−θmAt ∞ X n=0 (θmAt) n n! [I + θ −1 mACk] n≤ e−θmAt ∞ X n=0 (θmAt) n n! [I + θ −1 mAC ♮ k] n= exp[C♮ kt]. Furthermore, (2.24) follows from− Ck♮−diag(x)(−Ck♮)−1Dk = Dk+diag(x)(−Ck♮)−1Dk. Finally, (2.25) follows from
Z τ
0
exp Ck♮− diag(x)t[Dk+ diag(x)(−Ck♮)−1Dk]dt
= I − exp Ck♮− diag(x)τ ]− Ck♮− diag(x)−1[Dk+ diag(x)(−Ck♮)−1Dk] = I − exp Ck♮− diag(x)τ ](−Ck♮)−1Dk ≥ (I − exp[C♮ kτ ])(−C ♮ k)−1Dk= Z τ 0 exp[Ck♮t]Dkdt,
where (2.24) is used in the second equality and the inequality can be shown in the same way as (2.23) because (−Ck♮)−1Dk≥ O.
Theorem 2.3. We consider {(L(t), S(t))}t≥0 under Assumption 1.1 (ii), where (2.4) is assumed to hold. For an arbitrary integer K > 0, we have
max i∈M b + A,(K,i)(mA) ≤ 1 − mA−max(K,KX †2) n=0 Z ∞ 0 e−µt(µt) n n! dH(t), mA≥ max(K, K † 2),
where µ = 1/E[H] and K2† is given in Assumption 1.1 (ii).
Because the proof of Theorem 2.3 is a bit lengthy, it is given in Appendix A. Theorems 2.1–2.3 enable us to compute the approximation Aapprox,nA
k,ℓ (mA) (k, ℓ ≤ K) to Ak,ℓ, which satisfies (2.10) under Assumption 1.1. We summarize the procedures in Figure 1.
Next, we consider the computation of an approximation to Wk,ℓ’s using U (mW). We
first define W (mW) (mW > K) as W (mW) = Z ∞ 0 exp[U (mW)t] 1− H(t) dt = E[H]· BW(mW), where BW(mW) = Z ∞ 0 exp[U (mW)t] 1− H(t) E[H] dt.
Note here that 1−H(t)/E[H] is the probability density function of the equilibrium random variable for service times and therefore, we have
BW(mW) =
∞ X n=0
Input: {(Ck, Dk, Γk); k = 0, 1, . . .}, K, and ϵA. Output: Acomp,n∗A k,ℓ (k, ℓ = 0, 1, . . . , K). Let m∗A:= min mA> K; (−CK)−1DK(−CK+1)−1DK+1· · · DmAe≤ ϵA 2 e . Let n∗A:= min nA≥ 0; nA X n=0 γn(m∗A)≥ 1 − ϵA 2 . Compute Vk,ℓ(n)(m∗A)’s (n = 0, 1, . . . , n∗A) by (2.6). Compute Acomp,n∗A k,ℓ (m∗A) (k, ℓ = 0, 1, . . . , K) by (2.8) with (2.7), where Acomp,n∗A k,ℓ (m∗A) = O (k = 3, 4, . . . , K, ℓ = 1, 2, , . . . , k− 2).
(a) Under Assumption 1.1 (i) Input: {(Ck, Dk, Γk); k = 0, 1, . . .}, K, and ϵA. Output: Acomp,n∗A k,ℓ (k, ℓ = 0, 1, . . . , K). Let m∗A:= min mA> K; mA−max(K,KX 2†) n=0 Z ∞ 0 e−µt(µt) n n! dH(t)≥ 1 − ϵA 2 . Let n∗A:= min nA≥ 0; nA X n=0 γn(m∗A)≥ 1 − ϵA 2 . Compute Vk,ℓ(n)(m∗A)’s (n = 0, 1, . . . , n∗A) by (2.6). Compute Acomp,n∗A k,ℓ (m∗A) (k, ℓ = 0, 1, . . . , K) by (2.8) with (2.7), where Acomp,n∗A k,ℓ (m∗A) = O (k = 3, 4, . . . , K, ℓ = 1, 2, , . . . , k− 2).
(b) Under Assumption 1.1 (ii)
Figure 1: Computational procedures for Ak,ℓ’s in the non-uniformizable case
where ηn(mW) (n = 0, 1, . . .) denotes a probability function given by
ηn(mW) = Z ∞ 0 e−θmWt(θmWt) n n! 1− H(t) E[H] dt,
and V (mW) = I + θm−1WU (mW). In numerical implementation, we truncate the infinite sum
in (2.26) at n = nW. Bcomp,nW W (mW) = nW X n=0 ηn(mW)Vn(mW).
It then follows that for k = 1, 2, . . . , K, an approximation Wcomp,nW
k,ℓ (mW) to Wk,ℓ is given by Wcomp,nW k,ℓ (mW) = E[H] nW X n=1 ηn(mW)V (n) k,0 (mW), ℓ = 0, E[H] nW X n=ℓ−k ηn(mW)V (n) k,ℓ (mW), ℓ = k, k + 1, . . . , min(mW, k + nW), O, otherwise. (2.27)
We would like to find mW = m∗W and nW = n∗W such that for a given ϵW > 0, ∞ X ℓ=k Wk,ℓ− Wcomp,nW k,ℓ (mW) e = E[H] ∞ X ℓ=k BW,k,ℓ− BW,k,ℓcomp,nW(mW) e≤ ϵWe, (2.28)
for all k = 1, 2, . . . , K. Note here that Bcomp,nW
W,k,ℓ (mW) = O for all ℓ > min(mW, k + nW).
Since Bcomp,nW W,k,ℓ (mW)≤ BW,k,ℓ, we have ∞ X ℓ=k BW,k,ℓ− B comp,nW W,k,ℓ (mW) e≤ ∞ X ℓ=0 BW,k,ℓ− B comp,nW W,k,ℓ (mW) e. Therefore, (2.28) holds if ∞ X ℓ=0 BW,k,ℓ− BW,k,ℓcomp,nW(mW) e≤ ϵW E[H]e, k = 1, 2, . . . , K. (2.29) Because the left-hand side of (2.29) takes the same form as that of (2.17), we can obtain the following corollary according to the same line of discussions as for Theorems 2.1–2.3. Corollary 2.1. For arbitrary K, mW (0 < K < mW) and nW > 0, we have
∞ X ℓ=0 BW,k,ℓ− BW,k,ℓcomp,nW(mW) e≤ max i∈M b + W,(K,i)(mA) e + 1− nW X n=0 ηn(mW) e, k = 1, 2, . . . , K, where b+W,K(mW) = b+W,(K,1)(mW) b+W,(K,2)(mW) · · · b+W,(K,M )(mW) T is given by b+W,K(mW) = BW,k,0− BW,k,0(mW) e + ∞ X ℓ=mW+1 BW,k,ℓe.
Under Assumption 1.1 (i),
b+W,K(mW)≤ ϵe, (2.30)
holds for a given ϵ > 0 if mW satisfies (−CK)−1DK(−CK+1)−1DK+1· · · (−CmW)−1DmWe≤
ϵe. On the other hand, under Assumption 1.1 (ii), (2.30) holds for a given ϵ > 0 if mW
satisfies mW ≥ max(K, K2†) and
1− mW−max(K,KX 2†) n=0 Z ∞ 0 e−µt(µt) n n! · 1− H(t) E[H] dt≤ ϵ, where µ = 1/E[H] and K2† is given in Assumption 1.1 (ii).
Figure 2 shows the computational procedures for Wk,ℓ (k = 1, 2, . . . , K, ℓ = 0, k, k + 1, . . . , min(mW, k + nW)) satisfying (2.28). Once we compute Wk,ℓcomp,nW(mW)’s, we can
estimate the error contained in those by ∞ X ℓ=k Wk,ℓ− Wcomp,nW k,ℓ (mW) e≤ E[H]e − mW X ℓ=0 Wcomp,nW k,ℓ (mW)e, (2.31)
Input: {(Ck, Dk, Γk); k = 0, 1, . . .}, K, and ϵW.
Output: m∗W and Wcomp,n∗W
k,ℓ (k = 1, 2, . . . , K, ℓ = 0, k, k + 1, . . . , m∗W). Let m∗W := min mW > K; (−CK)−1DK(−CK+1)−1DK+1· · · DmWe≤ ϵW 2E[H]e . Let n∗W := min nW ≥ 0; nW X n=0 ηn(m∗W)≥ 1 − ϵW 2E[H] . Compute Vk,ℓ(n)(m∗W)’s (n = 0, 1, . . . , n∗W) by (2.6). Compute Wcomp,n∗W k,ℓ (m∗W) (k = 1, 2, . . . , K, ℓ = 0, k, k + 1, . . . , min(m∗W, k + n∗W)).
(a) Under Assumption 1.1 (i) Input: {(Ck, Dk, Γk); k = 0, 1, . . .}, K, and ϵW.
Output: m∗W and Wcomp,n∗W
k,ℓ (k = 1, 2, . . . , K, ℓ = 0, k, k + 1, . . . , m∗W). Let m∗W := min mW> K; mW−max(K,KX 2†) n=0 Z ∞ 0 e−µt(µt) n n! 1− H(t) E[H] dt≥ 1 − ϵW 2E[H] . Let n∗W:= min nW≥ 0; nW X n=0 ηn(m∗W)≥ 1 − ϵW 2E[H] . Compute Vk,ℓ(n)(m∗W)’s (n = 0, 1, . . . , n∗W) by (2.6). Compute Wcomp,n∗W k,ℓ (m∗W) (k = 1, 2, . . . , K, ℓ = 0, k, k + 1, . . . , min(m∗W, k + n∗W)).
(b) Under Assumption 1.1 (ii)
Figure 2: Computational procedures for Wk,ℓ’s in the non-uniformizable case
3. The Augmented Truncation Approximation to the Imbedded Markov Chain In this section, we discuss the augmented truncation approximation (ATA) to the imbedded Markov chain with the truncated transition probability matrix P (K). In the ATA, the selection of the augmentation matrix PA(K) in (1.6) is crucial and we briefly summarize
how to manage it according to [11]. First of all, we restrict our attention to the linear augmentation, i.e., PA(K) = (I− P (K))eξ for some (K + 1)M × 1 probability vector ξ,
which does not lose generality because for any augmentation matrix PA(K), there exists
ξ := ξ(PA(K)) that yields the same ATA solution [11, Implication 1]. We then define
κapprox(K; ξ) as the linear ATA solution obtained by
κapprox(K; ξ) = κapprox(K; ξ)P (K) + (I− P (K))eξ, κapprox(K; ξ)e = 1. Let [ξ](k,i)(k = 0, 1, . . . , K; i∈ M) denote the (kM +i)th element of ξ and let [P ](k,i),(ℓ,j)
(k, ℓ = 0, 1, . . .; i, j ∈ M) denote the (kM + i, ℓM + j)th element of P . We then define Ξ(K) as
Ξ(K) =ξ∈ R(K+1)M; ξ≥ 0, ξe = 1,
[ξ](k,i) > 0 if (k, i)∈ J+(K), [ξ](k,i) = 0 if (k, i) /∈ J+(K)
,
Input: P (K) andJ+(K).
Output: κapprox = (κapprox(K) 0) and ε(K).
Compute κapprox(K; eT(k,i)) = x(k,i)/x(k,i)e for all (k, i)∈ J+(K),
where x(k,i) is the unique solution of x(k,i)(I− P (K)) = eT(k,i). Compute κapprox(K) by (3.3) and set κapprox= (κapprox(K) 0). Compute the error bound ε(K) by (3.4).
Figure 3: Computational procedure for κapprox = (κapprox(K) 0)
at least one state in levels K + 1 or higher.
J+(K) =n(0, i); i∈ M, ∞ X k=K+1 [eTΓk]i > 0 o ∪(K, i); i∈ M . (3.1) Note that there exists ξ∗ ∈ Ξ(K) such that κapprox(K; ξ∗) = κ(K) [11, Remark 2.2].
Furthermore, for an arbitrary ξ ∈ Ξ(K), κapprox(K; ξ) is given by a convex combination
of κapprox(K; eT(k,i))’s ((k, i) ∈ J+(K)) with positive weights [11, Eq. (2.33)], where eT(k,i) denotes the unit row vector whose (kM + i)th element is equal to one. Specifically,
κapprox(K; ξ) = X
(k,i)∈J+(K)
[α(K; ξ)](k,i)κapprox(K; eT(k,i)),
where α(K; ξ) ∈ Ξ(K) and there is a one-to-one correspondence between ξ and α(K; ξ) [11, Eq. (2.34)]. It then follows that [11, Eq. (2.37)]
∥κapprox(K; ξ)− κ(K)∥
1 ≤ ε(K) = max (k,i)∈J+(K)∥κ
approx(K; ξ)− κapprox(K; eT
(k,i))∥1. (3.2)
Based on this observation, the following approximation κapprox(K) to the conditional sta-tionary distribution κ(K) is suggested in [11, Implication 4]:
κapprox(K) = X (k,i)∈J+(K) 1 |J+(K)|κ approx(K; eT (k,i)), (3.3)
which presumes that the transition structure in levels higher than K is unavailable, and the upper bound ε(K) of the error ϵ(K) in (3.2) is given by
∥κ(K) − κapprox
(K)∥1 ≤ ε(K) = max (k,i)∈J+(K)∥κ
approx
(K)− κapprox(K; eT(k,i))∥1. (3.4)
In summary, we compute the approximation κapprox(K) and its error bound ε(K) by the
procedure in Figure 3.
Remark 3.1. Unfortunately, we cannot show limK→∞ε(K) = 0 in general. Note, however, that if
Pr((Ln+1, Sn+1) = (ℓ, j)| (Ln, Sn) = (k, i)) > 0, for all (k, i)∈ {0, 1, . . . , K1†} × M, (ℓ, j) ∈ J+(K†
1), we can show limK→∞ε(K) = 0 under Assumption 1.1 (i). Due to the shortage of space, we omit the details.
4. The Error Bound in the Approximation to the Stationary Distribution If we obtain an approximation κapprox to the stationary distribution κ of the imbedded
Markov chain {(Ln, Sn)}n=0,1,... and Wk,ℓcomp,nW(mW)’s, we can compute an approximation
πapprox = (πapprox
0 π
approx
1 · · · ) to the stationary distribution π by (1.8) and (1.9), where
Wk,ℓ is replaced by Wk,ℓcomp,nW(mW). Since Wk,ℓcomp,nW(mW) = O (ℓ > min(mW, k + nW)),
the computation of πkapprox’s stops at ℓ = min(mW, K + nW). Recall that for an arbitrary
ϵW > 0, Wk,ℓcomp,nW(mW)’s satisfy (2.28). We thus discuss the error in πapprox, assuming that
errors in Wcomp,nW
k,ℓ (mW)’s are negligible, i.e.,
π = κ cW κ cW e = κ cW E[T ], π approx= κapproxWc κapproxW ec , where c W = (−C0)−1 (−C0)−1D0W1,1 (−C0)−1D0W1,2 (−C0)−1D0W1,3 · · · O W1,1 W1,2 W1,3 · · · O O W2,2 W2,3 · · · O O O W3,3 · · · .. . ... ... ... . .. ,
and E[T ] = κ cW e denotes the average length of two consecutive imbedded points. Theorem 4.1. If errors in Wcomp,nW
k,ℓ (mW)’s are negligible, we have
∥π − πapprox∥ 1 ≤ Tmaxupper κapproxW ec ∥κ − κ approx∥ 1 1 + T upper max κapproxW ec − Tupper max ∥κ − κapprox∥1 . (4.1) where Tupper
max = E[H] + maxi∈M[(−C0)−1e]i. Remark 4.1. (1.11) comes from (4.1) and
Tmaxupper κapproxW ec ∥κ − κ approx∥ 1 1 + T upper max κapproxW ec − Tupper max ∥κ − κapprox∥1 = T upper max κapproxW ec 1 + T upper max κapproxW ec ∥κ − κapprox∥ 1+ o(∥κ − κapprox∥1).
Proof. Let Tmax= max(k,i)∈{0,1,...}×M[ cW e](k,i). Because diag−1( cW e) cW e = e, we have
∥π − πapprox∥ 1 = κdiag( cW e) κ cW e diag −1( cW e) cW − κapproxdiag( cW e) κapproxW ec diag −1( cW e) cW 1 ≤ κdiag( cW e) κ cW e − κapproxdiag( cW e) κapproxW ec 1 ≤ (κ− κapprox)diag( cW e) κapproxW ec 1+ 1 κ cW e − 1 κapproxW ec ∥κdiag(cW e)∥1 ≤ Tmax κapproxW ec ∥κ − κ approx∥ 1+ (κapprox− κ)cW e κ cW e· κapproxW ec · Tmax ≤ Tmax κapproxW ec ∥κ − κ approx∥ 1+ T2 max κ cW e· κapproxW ec ∥κ − κ approx∥ 1. (4.2)
