# Computation of Invariant Measures and Stationary Expectations for Markov Chains with Block-Band Transition Matrix

Authors: Hendrik Baumann; Thomas Hanschke
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

This paper deals with the computation of invariant measures and stationary expectations for discrete-time Markov chains governed by a block-structured one-step transition probability matrix. The method generalizes in some respect Neuts’ matrix-geometric approach to vector-state Markov chains. The method reveals a strong relationship between Markov chains and matrix continued fractions which can provide valuable information for mastering the growing complexity of real-world applications of large-scale grid systems and multidimensional level-dependent Markov models. The results obtained are extended to continuous-time Markov chains.


---

## Body

## 1. Introduction

This paper is concerned with the computation of invariant measures of discrete-time Markov chainsX=Xkk≥0 having a block-partitioned one-step transition matrix P of the form
(1)P=P00P0100⋯⋯⋯⋯P10P11P120⋯⋯⋯⋯⋮⋮⋱⋱Pn0Pn1⋯PnnPn,n+100⋯0Pn+1,1⋯Pn+1,nPn+1,n+1Pn+1,n+20⋯⋱⋱⋱⋱⋱⋱,where all elements are nonnegative matrices of suitable dimension, that is, Pin∈ℝdi×dn for d0,d1,⋯∈ℕ=1,2,3,⋯. The states of the Markov chain are denoted by i,j, i≥0 and 1≤j≤di, and are ordered in lexicographic order, that is, 0,1, …, 0,d0, 1,1, …, 1,d1, 2,1, …. Following Neuts [1], the set of states i,1,⋯,i,di,i≥0, will be called the leveli. We refer to the first component of Xk (the level) as Xk1, and by Xk2, we denote the second component of Xk. Throughout this paper, it is assumed that n≥1 and that the Markov chain is irreducible and recurrent. Note that irreducibility implies Pi,i+1≠0 for all i∈0,1,2,⋯. From the general theory of Markov chains, we then know that the system
(2)xP=x,x≥0,has a unique solution (up to a constant factor) which is strictly positive (see Karlin and Taylor ([2], chapter 11)). For the solution vector x, we use the notation x=x0,x1,x2,⋯ and xi=xi,1,xi,2,⋯,xi,di,i≥0. Every constant positive multiple of x is called an invariant measure of X. The class of irreducible recurrent chains can be further subdivided into positive and null-recurrent chains according to whether ∑i=0∞xieT<∞ or =∞, with e being a vector with all its components equal to one.Vector-state Markov chains of the form (1) generalize the theory of matrix-geometric processes which was pioneered by M. Neuts [1]. As he points out in his monograph, the matrix-geometric approach is applicable to a wide class of stochastic models. The characteristic feature of these models is that the blocks within the transition probability matrix (or the generator matrix in case of continuous-time models) only depend on their relative position to the main diagonal, that is, Pij only depends on j−i. Otherwise, Neuts’ method fails.A general strategy for solving systems of the form (1) is the following:
(i)
Truncate the system at levelN, that is, consider the northwest corner truncation NP of P with N+1×N+1 blocks(ii)
Derive a stochastic matrixNP~ by appropriately augmenting NP(iii)
SolvexN·NP~=xN.ThenxN might converge to the desired solution x under some appropriate conditions (a natural condition is some normalization, e.g., x0,1N=1).In case of (1) with n=1 and Pij only depending on j−i, level-independent quasi-birth-death processes arise. For n=1 and arbitrary Pij, the quasi-birth-death process is level-dependent. Several publications have dealt with the problem of computing the invariant measures for this class of stochastic processes numerically: In [3] and in [4], algorithms based on taboo probabilities were developed. In [5], a method relying on matrix continued fractions was discussed. In principle, these three algorithms differ with respect to the question how to compensate the truncation of the infinite state space which is inevitable for algorithmic procedures, that is, how to augment NP. The major part of the algorithmic procedure is similar in all algorithms, and as pointed out in [6], for large truncation levels, it can be used without any augmentation. Since NP is not stochastic anymore, the last step of the above general method has to be interpreted appropriately: xN shall solve xN·NP=xN up to the first scalar equation. This simplified matrix-analytic algorithm becomes even more powerful when only stationary expectations have to be computed since then, the memory-efficient variant suggested in [7] can be applied.Recently, Masuyama published a technical report [8] in which Markov chains (with a focus on continuous-time Markov chains) with a general block structure are discussed. Under somewhat restrictive conditions, it is shown that the above method (truncation, augmentation, finding a solution of the truncated and augmented system) converges to the desired invariant measure if the truncation level tends to ∞. Since Markov chains with general block structure are discussed, transition matrices of the form (1) or its continuous analogue ((20) below) can be interpreted as special cases. However, the focus in [8] is the discussion of augmentation procedures.We will treat the problem of finding an invariant measure numerically with a different perspective: After the truncation, we perform no augmentation but findxN in a similar manner as in the method for quasi-birth-death processes described in [7]. In some situations, not using an augmentation might result in requiring a higher truncation level than with an appropriate augmentation (although this is not clear). However, we discuss not only a basic algorithm for computing invariant measures but also the direct computation of stationary characteristics for which we introduce a memory-efficient implementation. Due to the memory-efficiency, our method can be applied with a very high truncation level.In total, in some way, we combine the ideas of Neuts’ original matrix-geometric approach with considerations for level-dependent quasi-birth-death processes (although the resulting transition structures are not as general as in [8]), and with memory-efficient algorithms as developed in [7] for quasi-birth-death processes.The rest of the paper is organized as follows: We begin with a derivation of our results and methods based on taboo probabilities in Section2 and extend the results to continuous-time Markov chains in Section 3. The resulting algorithm for computing invariant measures is presented in Section 4 whereby we also derive a memory-efficient implementation for computing stationary expectations. Section 5 is dedicated to an alternate approach. This derivation is based on generalized continued fractions or matrix continued fractions, respectively. Although these considerations do not result in new computational methods, they give us a better insight in the algebraic and numerical aspects of the approach. Finally, in Section 6, we demonstrate the applicability of our methods by briefly discussing a queueing model with batch service and server breakdowns.


## 2. A Derivation Based on Taboo Probabilities

According to [2, 9, 10], an invariant measure x=xℓℓ∈E of an irreducible recurrent Markov chain X with state space E and transition probability matrix P=Piji,j∈E is given by xℓ0=1 and
(3)xℓ=∑k=1∞ℙXk=ℓ,Xk−1,⋯,X1≠ℓ0∣X0=ℓ0,for ℓ≠ℓ0, where ℓ0 is some fixed state and ℙA denotes the probability of event A. The summands on the right-hand side are taboo probabilities. Irreducibility of P guarantees that the series converges for all ℓ∈E. xℓ can be interpreted as the expected number of visits to state ℓ before the Markov chain returns to its initial state ℓ0. Note that x is the unique invariant measure of P with xℓ0=1.In the present case, we chooseℓ0=0,1 and obtain
(4)xi,j=∑k=1∞ℙXk=i,j,Xk−1,⋯,X1≠0,1∣X0=0,1,for i,j≠0,1. By σN=infk∈1,2,3,⋯:Xk1≥N, we denote the first hitting time of X of some level ≥N.Following [9, 10], we have
(5)xi=limN→∞1h0,1NhiNi=0,1,2,⋯,where hN=h0N,⋯,hNN solves
(6)hNI−NP=e1.NP denotes the northwest corner truncation of P at level N, I is an identity matrix of appropriate size, and e1 is the vector with unity in the first position and zeroes elsewhere. The solution is unique, and with xiN=1/h0,1NhiN, we have x0,1N=1 and
(7)xi,jN=∑k=1∞ℙXk=i,j,Xk−1,⋯,X1≠0,1,σN+1>k∣X0=0,1,for i≤N with i,j≠0,1. Obviously, xi,jN increases monotically in N, and by definition, we have xi=limN→∞xiN.Furthermore, we setτi=inf{k∈1,2,3,⋯:Xk1≤i to be the first positive hitting time of some level ≤i (note that due to irreducibility and recurrence, τi<∞ almost surely), and we define the matrices
(8)Ri=ri,j,j′1≤j≤di1≤j′≤di+1=∑k=1∞ℙXk=i+1,j′,τi>k∣X0=i,j1≤j≤di1≤j′≤di+1,RiN=∑k=1∞ℙXk=i+1,j′,τi>k,σN+1>k∣X0=i,j1≤j≤di1≤j′≤di+1.Again, the entries ofRi can be interpreted as expectations. In this case, conditioned on the initial state i,j, the entry at position j,j′ is the expected number of visits in i+1,j′ before returning to some level ≤i. Similarly, the entries of RiN are expectations for the expected number of visits in i+1,j′ before returning to some level ≤i or hitting some level >N. Obviously, RiN=0 for i≥N. Again, componentwise convergence of RiN⟶Ri is a direct consequence of the monotone convergence theorem. We are now ready to find the representation of the invariant measures which will be exploited numerically later on.Theorem 1.
With the structure (1) ofPand the notations introduced above, we have(9)xi+1=xiRii=0,1,2,⋯,(10)xi+1N=xiNRiNi=0,1,2,⋯,N−1,(11)x0=x0∑i=0n∏r=0i−1RrPi0,(12)Ri=Pi,i+1I−∑j=i+1i+1+n∏r=i+1j−1RrPj,i+1−1i=0,1,2,⋯,(13)RiN=Pi,i+1I−∑j=i+1i+1+n∏r=i+1j−1RrNPj,i+1−1i=0,1,2,⋯,N−1.
In (12) and (13), the inverse matrices exist. In (11), the solution of the system of linear equations is unique up to constant multiples, and additionally,(14)x0N=x0N∑i=0n∏r=0i−1RrNPi0holds up to the first scalar equation.Proof.
The interpretations ofRi are very similar to the level-independent case (see [1]) and to quasi-birth-death processes (see [4, 5]). Therefore, we do not state all proofs in detail.
(i)
Fori≥1, we have
(15)xi+1,j′=∑k=2∞ℙXk=i+1,j′,Xk−1,⋯,X1≠0,1∣X0=0,1=∑k=2∞∑ℓ=0k−1ℙXk=i+1,j′,Xk−11,⋯,Xℓ+11>iXℓ,⋯,X1≠0,1X0=0,1=∑k=2∞∑ℓ=1k−1∑j=1mℙXk=i+1,j′,Xk−11,⋯,Xℓ+11>iXℓ=i,j,Xℓ−1,⋯,X1≠0,1X0=0,1=∑j=1m∑k=2∞∑ℓ=1k−1ℙXℓ=i,j,Xℓ−1,⋯,X1≠0,1∣X0=0,1·ℙXk=i+1,j′,Xk−11,⋯,Xℓ+11>i∣Xℓ=i,j=∑j=1m∑ℓ=1∞∑k=ℓ+1∞ℙXℓ=i,j,Xℓ−1,⋯,X1≠0,1∣X0=0,1·ℙXk−ℓ=i+1,j′,Xk−ℓ−11,⋯,X11>i∣X0=i,j=∑j=1m∑ℓ=1∞ℙXℓ=i,j,Xℓ−1,⋯,X1≠0,1∣X0=0,1·∑k=1∞ℙXk=i+1,j′,τi>k∣X0=i,j=∑j=1mxi,jri,j,j′.This proves (9) for i≥1. For i=0, we have an additional summand for k=1, but on the other hand, for k≥2, the sum regarding j starts with j=2. It turns out that (9) still holds. (10) is proven analogously.(ii)
Using the structure (1) of P and (9) repetitively, it becomes clear that x0 satisfies equation (11) which can be rewritten as x0=x0K0, where
(16)Ki=∑j=ii+n∏r=ij−1RrPjii=0,1,2,⋯.By a calculation similar to that for proving (9), we see that the entry of Ki at position j,j′ is given by the probability
(17)∑k=1∞ℙXk=i,j′,τi≥k∣X0=i,j=ℙτi<∞,τi<τi−1,Xτi=i,j′∣X0=i,j,that the Markov chain returns to level i before hitting a level <i and that this return is to i,j′. In particular, the entries of K0 represent the probabilities that the Markov chain eventually returns to level 0 and that this return occurs in 0,j′. Irreducibility and recurrence of X imply stochasticity and irreducibility (and due to K0 being finite, recurrence) of K0. Hence, x0 is an invariant measure for K0 which is uniquely defined up to a constant factor.(iii)
Alternatively, (11) can be shown by arguments similar to those leading to (9). These arguments also confirm that x0N satisfies the corresponding equation up to the first scalar equation(iv)
Additionally, define(18)Ji=∑k=1∞ℙXk=i,j′,τi−1>k∣X0=i,jj,j′∈E,as the matrix of the expectations of the number of visits in i,j′ before hitting some level <i. Then considerations similar to those yielding (9) lead to I+JiI−Ki=I−KiI+Ji=I, and hence, I+Ji=I−Ki−1 for i≥1 and Ri=Pi,i+1I+Ji+1 for i≥0.
The fact thatRiN converges monotonically increasing to Ri entails that the Ri are minimal nonnegative solutions to the equations
(19)RiI−∑j=i+1i+1+n∏r=i+1j−1RrPj,i+1=Pi,i+1,which corresponds to similar results in [1]. However, we will not make use of this fact in Algorithms 1 and 2 which we suggest for computing invariant measures and stationary expectations.


## 3. Continuous-Time Markov Chains

For continuous-time Markov chainsY=Ytt≥0 with block-structured generator matrix
(20)Q=Q00Q0100⋯⋯⋯⋯Q10Q11Q120⋯⋯⋯⋯⋮⋮⋱⋱Qn0Qn1⋯QnnQn,n+100⋯0Qn+1,1⋯Qn+1,nQn+1,n+1Qn+1,n+20⋯⋱⋱⋱⋱⋱⋱,invariant measures satisfy yQ=0. For the construction of a continuous-time Markov chain from its infinitesimal generator as well as for precise definitions of properties such as regularity, recurrence, and ergodicity, the reader is referred to [11, 12]. Assuming regularity, irreducibility, and recurrence, the solution to yQ=0 is unique up to constant multiples and strictly positive. Furthermore, these assumptions allow constructing the embedded jump chain which is a discrete-time Markov chain X=Xnn=0∞ with transition probability matrix with structure (1) and Pij=Di−1Qij+δijI where Di is the diagonal matrix with entries −qi,j,i,j at position j,j. Irreducibility of Y implies −qi,j,i,j>0 for all i,j∈E, and hence, Di is indeed invertible for all i∈0,1,2,⋯. Furthermore, irreducibility and recurrence of Y imply irreducibility and recurrence of X and vice versa (note that this is not true for ergodicity). Finally, if x is invariant for P, an invariant measure y for Q can be constructed by setting yi=xiDi−1.Renaming the matrices from Section2 to R~i, we have xiR~i=xi+1 for the invariant measure of the embedded jump chain. Hence, we obtain
(21)yi+1=xi+1Di+1−1=xiR~iDi+1−1=yiDiR~iDi+1−1=yiRi,with Ri=DiR~iDi+1−1. From (11), we obtain
(22)y0=x0D0−1=x0∑i=0n∏r=0i−1R~rPi0D0−1=y0D0∑i=0n∏r=0i−1Dr−1RrDr+1Di−1Qi0+δi0ID0−1=y0∑i=0n∏r=0i−1RrQi0+y0.Taking into account thatR~i satisfies (12), we arrive at
(23)Ri=DiR~iDi+1−1=DiPi,i+1I−∑j=i+1i+1+n∏r=i+1j−1R~rPj,i+1−1Di+1−1=DiDi−1Qi,i+1Di+1−Di+1∑j=i+1i+1+n∏r=i+1j−1Dr−1RrDr+1Dj−1Qj,i+1+δj,i+1I−1=−Qi,i+1∑j=i+1i+1+n∏r=i+1j−1RrQj,i+1−1i=0,1,2,⋯.Since thej′th diagonal entry of Di+1−1 can be interpreted as the mean sojourn time in i+1,j′, the entries of R~iDi+1−1 are the expected time intervals spent in state i+1,j′ before returning to some level ≤i, conditioned on starting in i,j. Therefore, the entries of Ri represent the proportion of the expected times spent in i+1,j′ before returning to some level ≤i, and the time spent in the initial state i,j.Summarizing, we getTheorem 2.
With the above assumptions onQ, we have(24)yi+1=yiRii=0,1,2,⋯,(25)yi+1N=yiNRiNi=0,1,2,⋯,N−1,(26)0=y0∑i=0n∏r=0i−1RrQi0,where the solution of the system of linear equations in (26) is unique up to constant multiples. The Ri’s satisfy
(27)Ri=−Qi,i+1∑j=i+1i+1+n∏r=i+1j−1RrQj,i+1−1i=0,1,2,⋯,and can be uniquely constructed by Ri=limN→∞RiN for i=0,1,2,⋯, where RiN=0 for i≥N and
(28)RiN=−Qi,i+1∑j=i+1i+1+n∏r=i+1j−1RrNQj,i+1−1i=0,1,2,⋯,N−1.y0N satisfies
(29)0=y0N∑i=0n∏r=0i−1RrNQi0,up to the first scalar equation.


## 4. Algorithms

Next, we demonstrate how to use the results numerically. We focus on continuous-time Markov chains since the methods for discrete-time Markov chains follow an analogous scheme.Basically, we have to use (28) for computing approximations to RN−1,⋯,R0, then we determine an approximation for y0 by solving (26) up to the first scalar equation, and finally, we use (24) for computing approximations to y1,⋯,yN. In order to save operations, we suggest the following algorithm for computing an approximation to the invariant distribution.Algorithm 1: Computation of an invariant measure.
ChooseN largeInitializeHj=−QNj for j=N−n,⋯,Nfori=N−1,N−2,⋯,0 doComputeRi=Qi,i+1Hi+1−1Free memory space forHi+1UpdateHj=−Qij+RiHj for j=maxi−n+1,0,⋯,iInitializeHi−n=−Qi,i−n if i−n>0end forSety0,1=1 and determine all other entries of y0 by solving y0H0 up to the first scalar equation.Fori=0,1,⋯,N−1, compute yi+1=yiRireturn yIn case of positive recurrence, an invariant distribution can be obtained afterwards by renormalizingy. We give some remarks:
(i)
Algorithm 3 in [8] is quite similar. While our algorithm relies on the nonstochastic matrix P=Pini,n=1N or the nonconservative matrix Q=Qini,n=1N, respectively, in [8], an augmentation is performed first in order to obtain a stochastic or conservative matrix. Note that Theorem 1 guarantees that for N⟶∞, we have convergence of the approximations xN or yN to the desired invariant measure without any augmentation. The only use of augmentations is that the truncation error might be decreased a bit for small values of the truncation levelN. Instead of considering an augmentation, we suggest to choose N very large. If such a choice results in computational problems, we suggest to use the memory-efficient algorithm for computing stationary expectations (see below). Furthermore, note that the results in [8] were derived under different conditions: While Pin≠0 is allowed for all i≤n+1 (lower-block Hessenberg matrices), some restrictive conditions (Assumptions 2.1 and 2.2 in [8]) are assumed(ii)
A major problem is the choice of the truncation levelN. There are some approaches concerning a priori estimates of the truncation error (e.g., [13]), but still the problem is not completely solved. Instead of using such error estimates, alternatively, as in [8], we might increase the truncation level iteratively, that is, we start with some small N and increase N by 1 in each step until yN−yN−1<ε for some fixed norm · and some fixed small ε>0. Note that for each new value of N, we have to restart Algorithm 1. Furthermore, note that a small difference yN−yN−1 does not imply that yN−y is small where y is the true invariant measure(iii)
In many applications, the state space is finite with a transition structure corresponding to the northwest corner of (1) or (20), respectively. Then Algorithm 1 is exact up to numerical deviations(iv)
Instead of computing the inverse matrixHi+1−1 first and then perform the matrix multiplication for determining Ri=Qi,i+1Hi+1−1, it is recommended to solve the system RiHi+1=Qi,i+1 directly(v)
The matricesQij can be generated when they are used for computation. Hence, there is no need to store these matrices(vi)
We have to store aboutN matrices at the same time. If the dimensions di of the matrices are constant, say di=d, the memory requirement is approximately Nd2 real numbers(vii)
Concerning the computational effort, we have to solve approximatelyN linear systems of equations of the form AX=B where the coefficients A and B and the solution X are matrices of dimension d×d (for constant di=d). In addition, we have to perform N·n matrix multiplications. The computational effort is approximately γ1Nd3+γ2Nndα with some constants γ1,γ2 and α∈2,3 (the exponent α and the constants γ1,γ2 depend on the method of choice, we have α=3 for “conventional matrix multiplication” and α=2.807 for Strassen’s algorithm, see [14])(viii)
Forn=1, Algorithm 1 more or less coincides with the methods proposed in [3–5]. One might expect that also in general, a discrete-time Markov chain with transition probability matrix (1) and bandwidth >3 or a continuous-time Markov chain with generator matrix (20) and bandwidth >3, respectively, may be reduced to the block-tridiagonal case by appropriately rearranging the states, i.e., define a new level 0 as the union of the original levels 0,⋯,n−1and define a new level 1 as the union of the original levels n,⋯,2n−1, …. But the increased block size heavily impacts the performance of the algorithm. If we apply the same truncation procedure, the truncation level for the modified matrix structure is ≈N/n, and for di=d, all matrices involved have the dimension nd×nd. This results in a computational effort of approximately γ1N/nnd3+γ2N/nndα=γ1Nn2d3+γ2Nnα−1dα. Since α>2, the computational effort for the transformed system becomes significantly larger than the effort for the original system. Similarly, the memory requirement of Nd2 is replaced by N/n·nd2=Nnd2. For these reasons, one should discard this supposedly simple “trick” (compare ([15], pp.~73ff)).In applications, often only stationary expectations have to be considered: Consider the case of positive recurrence, and denote byπ the invariant distribution. Then the ergodic theorem for discrete-time Markov chains guarantees that
(30)limK→∞1K∑k=0K−1fXk=πf=xfx1,for any function f:E⟶ℝ for which πf≔∑ℓ∈Eπℓfℓ converges absolutely. Again, by 1, we refer to a vector of appropriate dimension consisting of 1s. Analogously, for continuous-time Markov chains, we have
(31)limT→∞1T∫0TfXtdt=πf=yfy1,if πf converges absolutely. Now, let f=f1,⋯,fL,fL+1:E⟶ℝL+1 be a vector valued function with fL+1ℓ=1 for all ℓ∈E. Then for continuous-time Markov chains, we have πfℓ=yfℓ/yfL+1 for ℓ=1,⋯,L. Hence, computing yf allows to determine the stationary characteristics πfℓ which coincide with the long-run costs or rewards measured by fℓ.When computingyf, we can apply a Horner-type scheme which allows us not to store the matrices Ri and therefore save quite a lot memory requirement. This idea is very similar to that in [7]. With
(32)fi=f1i,1⋯fL+1i,1⋮⋱⋮f1i,di⋯fL+1i,di,we write
(33)yf≈∑j=0NyjNfj=y0N∑j=0n∏r=0j−1RrNfj=y0Nf0+R0Nf1+R1Nf2+R2N⋯+RN−2NfN−1+RN−1NfN=y0NZ0N,where
(34)ZiN=∑j=in∏r=ij−1RrNfj=fi+RiNfi+1+Ri+1N⋯+RN−2NfN−1+RN−1NfN.Obviously, we haveZNN=fN and
(35)ZiN=fi+RiNZi+1Ni=0,⋯,N−1.The matricesRiN are computed by a backward scheme, and (35) is a backward scheme, too. Since we are only interested in computing Z0N, there is no need to store the matrices RiN or ZiN. We suggestAlgorithm 2: Computation of stationary expectations.
ChooseN largeInitializeHj=−QNj for j=N−n,⋯,N and Z=fNfori=N−1,N−2,⋯,0 doUpdateR=Qi,i+1Hi+1−1UpdateZ=fi+RZFree memory space forHi+1UpdateHj=−Qij+RHj for j=maxi−n+1,0,⋯,iInitializeHi−n=−Qi,i−n if i−n≥0end forSety0,1=1 and determine all other entries of y0 by solving y0H0 up to the first scalar equationUpdateZ=y0Z and renormalize Z=z1,⋯,zL+1 such that zL+1=1return ZAgain, we give some comments:(i)
The problem of finding an appropriate truncation levelN to keep the error in prescribed limits is still to be solved. In [16], a generalization of the idea from [13] is discussed, and ∑i=0Nyifi and ∑i=0∞yifi are compared. Unfortunately, the method is quite cumbersome to apply, especially because it would be more desirable to compare ∑i=0∞yifi with ∑i=0NyiNfi. Again, an alternative consists in iteratively increasing N until ∑i=0NyiNfi−∑i=0N−1yiN−1fi becomes small. Note that again, this method is somehow problematic since the value of this difference does not tell anything about the true error ∑i=0NyiNfi−∑i=0∞yifi(ii)
For finite state spaces with the transition structure under consideration, Algorithm2 is exact up to numerical deviations(iii)
Of course, it is still recommended to solve the systemRHi+1=Qi,i+1 directly, and still, the matrices Qij should be generated when they are needed(iv)
In total, we have to store the matricesHi and Z. Since we need n matrices Hi at the same time, the total amount of memory requirement is Ond2 if the matrix dimensions di=d are constant. In particular, the memory requirement does not depend directly on the truncation level N(v)
As pointed out above, in ([8], section 3.2), lower-block Hessenberg matrices (block GI/M/1-type matrices) are considered, that is, Pij=0 for j>i+1, but Pij≠0 is allowed for all j<i. All our theoretical considerations carry over (we only have to replace n by ∞ or n by N at some points). Unfortunately, when applying Algorithm 2 to such block matrices, we still have to store N matrices, and hence, it is not memory-efficient anymore. Since the memory-efficient implementation makes the theoretical considerations applicable to a wide range of stochastic models, we have restricted ourselves to considering band-block-structured matrices(vi)
The considerations concerning the computational effort are similar to those for Algorithm1. In particular, we do not recommend to reduce the structure to the case n=1(vii)
Note that often, “natural” state spaces are multidimensional. However, if there is a partitionE=E0∪E1∪E2∪⋯ into finite levels Ei such that transitions with positive probability/transition rate increase the level Ei at most by 1 and decrease it at most by n, an appropriate renumbering of the states results in a structure given by (1) or (20). We refer to the end of Section 6 for an example and to [7, 13, 17] for the corresponding discussion concerning quasi-birth-death processes (that is, n=1).


## 5. Relationship to Matrix Continued Fractions

This section is devoted to the interdependence between Markov chains and matrix continued fractions. The results obtained are of importance not only from an inner mathematical but also from a numerical point of view. Matrix continued fractions are generalizations of ordinary continued fractions and generalized continued fractions (orn-fractions, see [18, 19]). We will show that with each transition matrix P of the form (1), one can associate a matrix continued fraction (or a matrix-driven Jacobi-Perron algorithm) being convergent if X is irreducible and recurrent. In addition, it is proven that its approximants and its limits coincide with the coefficients of the recurrence relations (10) and (9), respectively.A milestone in pure mathematics has been the arithmetic characterization of quadratic irrationalities which has been discovered by J. L. Lagrange (see [20]). Lagrange’s theorem states that the expansion of a real number as a simple continued fraction is periodic if and only if it is the real root of a quadratic equation in one variable with rational coefficients. Since then, mathematicians like C. G. J. Jacobi [21], O. Perron [22], P. Bachmann [23], L. Bernstein [24], and F. Schweiger [25] tried to expand Lagrange’s approach to higher degree irrationalities. G. C. F. Jacobi and O. Perron developed a tool for approximating real algebraic irrationalities of degree n≥3. To the honor of Jacobi and Perron, this algorithm has been named as the Jacobi-Perron algorithm. Although the Jacobi-Perron algorithm does not lead to the desired results, it has been successfully addressed to a lots of other unsolved problems in approximation theory and numerical analysis. In particular, in connection with modern computer technology, the Jacobi-Perron algorithm turns out to be an efficient tool for identifying and calculating subdominant solutions of higher order linear vector recursions.In the sequel, we will show that the schemes (10) and (13) are equivalent to a specific matrix continued fraction or a matrix-driven Jacobi-Perron algorithm associated with the transition probability matrix (1). For this purpose, we assume that di=d for all i=0,1,2,⋯ and that Pi,i+1 is invertible for all i=0,1,2,⋯. Define
(36)Ai+n−kk=−Pi+n−k,iPi−1,i−1fork=0,1,⋯,n−1;i=1,2,⋯,−Pii−IPi−1,i−1fork=n;i=1,2,⋯,where I denotes the unit matrix. Then (36) and (2) yield the n+1th-order homogeneous linear vector recursion
(37)xi−1=xiAin+xi+1Ai+1n−1+⋯+xi+nAi+n0i=1,2,⋯,subject to the boundary condition
(38)x0P00−I+x1P10+⋯+xnPn0=0.In this perspective, the problem at hand is to find the “right” initial values, sayxl,xl+1,⋯,xl+n,l sufficiently large, such that xi>0 for all i. In order to find a solution to this problem based on generalized continued fractions, we generalize and combine the results from [5, 26].Remember that(39)xi=limN→∞Nh0,1N−1hii=0,1,2,⋯,where Nh=Nh0,⋯,NhN solves the truncated system (6). (36) and (6) lead to the finite difference scheme
(40)Nhi−1=NhiAin+Nhi+1Ai+1n−1+⋯+Nhi+nAi+n0i=1,2,⋯,N,coupled to the boundary conditions:
(41)e1=Nh0P00−I+Nh1P10+⋯+NhnPn0,(42)NhN+1=NhN+2=⋯=NhN+n=0,where e1 is a 1×m vector with unity in the first position and zeroes elsewhere. For simplification we transform the system (40) into a first order scheme of the form
(43)Nhi−1=Nhi·Aii=1,⋯,N,by putting
(44)Nhi=Nhi,Nhi+1,⋯,Nhi+ni=0,1,2,⋯,Ai=AinI0⋯00Ai+1n−10I⋯00⋮⋮⋮⋮⋮Ai+n−1100⋯0IAi+n000⋯00i=1,2,⋯.Now, letℓ≤N be fixed. By repeated substitutions, we obtain from (43)
(45)Nhℓ=NhNΦℓN−ℓ,with
(46)ΦℓN−ℓ=ANAN−1·⋯·Aℓ+1.For the remainder, it will be convenient to partitionΦℓi as
(47)Φℓi=Φℓ,i00Φℓ,i01⋯Φℓ,i0nΦℓ,i10Φℓ,i11⋯Φℓ,i1n⋮⋮⋮Φℓ,in0Φℓ,in1⋯Φℓ,inni=0,1,2,⋯,where the matrices Φℓ,irs are of dimension d×d. In view of the boundary condition (42), we infer from (45) that
(48)Nhℓ+k=NhNΦℓ,N−ℓ0kk=0,1,⋯,n.For computing thed×d-matrices Φℓ,N−ℓ0k, the following result reduces the computational effort.Theorem 3.
LetAi+n−k+ℓk,k=0,1,⋯,n;i=1,2,⋯, be defined as in (36) and let, for convenience,(49)Ai+ℓk=Ifork=0;i=1,⋯,n;0fork=1,⋯,n−1;i=1,⋯,n−k.
Then(50)Φℓ,i0k=Ai+ℓnΦℓ,i−10k+Ai+ℓn−1Φℓ,i−20k+⋯+Ai+ℓ0Φℓ,i−n−10kk=0,1,⋯,n;i=1,2,⋯,where
(51)Φℓ,−j00=Iforj=0,0forj=1,⋯,nΦℓ,−j0k=Iforj=n−k+10elsewherek=1,⋯,n;j=0,1,⋯,nProof.
From (46), it is clear that the matrices Φℓi can be recursively computed by
(52)Φℓi=Ai+ℓΦℓi−1i=1,2,⋯,with initial condition
(53)Φℓ0=I.
Rewriting (52) in expanded form we get
(54)Φℓ,i00=Ai+ℓnΦℓ,i−100+Φℓ,i−110Φℓ,i01=Ai+ℓnΦℓ,i−101+Φℓ,i−111⋮Φℓ,i10=Ai+ℓ+1n−1Φℓ,i−100+Φℓ,i−120Φℓ,i11=Ai+ℓ+1n−1Φℓ,i−101+Φℓ,i−121⋮Φℓ,in0=Ai+ℓ+n0Φℓ,i−100Φℓ,in1=Ai+ℓ+n0Φℓ,i−101⋮i=1,2,⋯.
By repeated substitutions, it is readily seen that the matricesΦℓ,i0k satisfy the recurrence relation (50) at least for i=n+1,n+2,⋯. To proceed in this manner we assume, for the moment, that the recurrence relation (50) may be continued for i=0,−1,⋯,−n+1, without affecting the initial condition (53) and the scheme (54). Since the matrices Aℓ,⋯,A−n+1+ℓ do not appear in formula (46), they can be chosen arbitrarily. By putting
(55)A−j+ℓ=A−j+ℓnI0⋯0A−j+1+ℓn−10I⋯0⋮⋮⋮⋮A−j+n−1+ℓ100⋯IA−j+n+ℓ000⋯0=M=0I0⋯000I⋯0⋮⋮⋮⋮000⋯II00⋯0j=0,⋯,n−1,and Φℓ−n=M we get
(56)Φℓ−n+1=M2=00I0⋯0000I⋯0⋮⋮⋮⋮⋮I000⋯00I00⋯0,⋯,Φℓ−1=Mn=00⋯00II0⋯000⋮⋮⋮⋮⋮00⋯I0000⋯0I0,which explains (49) and (51). Since Φℓ0=Mn+1=I, (55) becomes compatible with (52) and (53), which proves our result for i=1,⋯,n.To make (48) congruent with (9) and (10), respectively, we first prove that detΦi−1,N−i−100≠0 for almost all N under a mild condition on P.Consider the recurrence relation (40) for i=ℓ+1,ℓ+2,⋯,N. Assembling the unknowns Nhℓ,Nhℓ+1,⋯,NhN into an 1×m·N−ℓ+1 vector and setting the boundaries NhN=v,NhN+1=NhN+2=⋯=NhN+n=0, (40) permits a presentation of the form
(57)Nhℓ,Nhℓ+1,⋯,NhNW=0,⋯,0,v,where
(58)Since we assumeddetPi,i+1≠0 for all i, W is a regular matrix. For its inverse, we choose the partition
(59)W−1=XYUV,where U is a square matrix of order d. From (48), we know that
(60)Nhℓ=NhNΦℓ,N−ℓ00.Solving (57) successively with v replaced by 1,0,⋯,0, 0,1,⋯,0, …, 0,⋯,0,1, we get
(61)Φℓ,N−ℓ00=U.The matrixU is called the algebraic complement of the matrix C. A theorem which is due to Jacobi (see ([27], pp. 123-124, formula (4.45))) states that detU≠0 iff detC≠0. Since we assumed detPN,N+1≠0, C is weakly diagonally dominant. To guarantee that C is regular, it is therefore sufficient to assume that C is irreducible (see ([28], p. 531, Theorem 3)).From (48), we conclude that
(62)Nhℓ+k=NhNΦℓ,N−ℓ0kk=1,⋯,nandNhN=NhℓΦℓ,N−ℓ00−1,implying
(63)Nhℓ+k=NhℓRℓ,ℓ+kNk=1,⋯,n;ℓ=0,1,⋯,N,where
(64)Rℓ,ℓ+kN=Φℓ,N−ℓ00−1Φℓ,N−ℓ0kk=1,⋯,n;ℓ=0,1,⋯,N.From (63) and (41), we obtain the decoupled recursion
(65)Nhℓ+1=NhℓRℓ,ℓ+1Nℓ=0,1,⋯,N,and the boundary condition
(66)Nh0P00−I+R01NP10+⋯+R0NNPn0=e1,which enable us to compute Nhℓ for ℓ=0,1,⋯,N. In view of the representations (49), (50), and (51), we make use of the notation
(67)Rℓ,ℓ+1NRℓ,ℓ+2N⋮Rℓ,ℓ+nN=Aℓ+10Aℓ+20Aℓ+30⋯AN0Aℓ+11Aℓ+21Aℓ+31⋯AN1⋯⋯⋯⋯⋯Aℓ+1n−1Aℓ+2n−1Aℓ+3n−1⋯ANn−1Aℓ+1nAℓ+2nAℓ+3n⋯ANn=II⋯IAℓ+n+10⋯AN000⋯Aℓ+n1Aℓ+n+11⋯AN1⋯⋯⋯⋯⋯⋯⋯0Aℓ+2n−1⋯Aℓ+nn−1Aℓ+n+1n−1⋯ANn−1Aℓ+1nAℓ+2n⋯Aℓ+nnAℓ+n+1n⋯ANn.(67) can be interpreted as the Nth convergent of the generalized matrix continued fraction (GCF) or the matrix-driven Jacobi-Perron algorithm associated with the limiting process
(68)Rℓ,ℓ+k=limN→∞Rℓ,ℓ+kNk=1,⋯,n,compare with [22, 29, 30]. If the limits (68) do exist, we obtain
(69)xℓ+1=xℓRℓ,ℓ+1,ℓ=0,1,2,⋯,x0P00−I+R01P10+⋯+R0nPn0=0,as a consequence of (65), (39), and (38).For the calculation ofRℓ,ℓ+kN and Rℓ,ℓ+k, we recommend the following algorithms:
(a)
Forward algorithm (FA): this is nothing else than a straightforward application of the relations (64), (49), (50), and (51).(b)
Backward algorithm (BA): since a matrix-driven Jacobi-Perron algorithm may be interpreted as a special case of a matrix continued fraction, its Nth convergent Rℓ,ℓ+kN=Φℓ,N−ℓ00−1Φℓ,N−ℓ0k, k=1,⋯,n, can be calculated from a linear fractional transformation (see ([29], Theorem 2.2)). Restricted to the matrix-driven Jacobi-Perron algorithm, this algorithm runs as follows. For N sufficiently large, define
(70)RN+1,N+1+kN=0k=1,⋯,n.Then(71)Rℓ,ℓ+1N=Bℓ+Rℓ+1,ℓ+1+nN−1Aℓ0ℓ=N,N−1,⋯,0,(72)Rℓ,ℓ+kN=Bℓ+Rℓ+1,ℓ+1+nN−1Aℓk−1+Rℓ+1,ℓ+kNk=2,⋯,n;ℓ=N,N−1,⋯,0,(73)Rℓ,ℓ+k=limN→∞Rℓ,ℓ+kNk=1,⋯,n;ℓ=0,1,2,⋯.From a numerical point of view,BA requires less computational effort than FA. The only disadvantage of BA is that for each convergent, the whole algorithm has to be repeated which can be time-consuming if the algorithm converges slowly. Since we are not interested in a particular Ri,i+k but in the family of matrices Ri,i+k, we take BA for the preferable algorithm in all situations.It is easy to see that the recurrence relations (71)–(73) are equivalent to (12) and (13) with Rℓ,ℓ+k=Rℓ·Rℓ+1·⋯·Rℓ+k−1. Hence, the derivation based on continued fractions leads to the same results as the derivation based on taboo probabilities.Observe that the derivation based on taboo probabilities still works ifdetPi,i+1≠0 for some i, whereas the derivation based on matrix continued fractions fails. Nevertheless, the derivation based on matrix continued fractions enables us to detect that (65) is a decoupled version of the difference equation (40) which is important with respect to numerical applications. In [31], it was shown that for n=1 and d=1, the difference equation (37) is subject to numerical instability, because the invariant measure x is asymptotically dominated by other solutions of the difference equation (see [32]). In addition, it was shown that the continued fraction algorithms (65) and (71) result in a decoupled recursion for the computation of x in which the dominant solutions of (37) do not appear. But however, we cannot assure that in general, the decoupled recursion (65) is kept free from all sources of numerical instability. To take this issue forward, we need tailored Pringsheim- and Pincherle-type convergence criteria for matrix continued fractions, compare with [29].Remark 4.
The existence of the limits (68) has been shown in Theorem 1 by probabilistic arguments. To prove (68) in a purely analytic manner, further research is required. The reason is the following: Obviously, the convergence of the matrix-driven Jacobi-Perron algorithm is implied by the diagonal dominance of the transition probability matrix P. In general, there exists a great variety of convergence criteria for ordinary continued fractions and their generalizations which appeal to the diagonal dominance of their coefficients (the so-called Pringsheim criteria). But in the case of the matrix-driven Jacobi-Perron algorithm, these results are only available in terms of norm-based criteria, which are too restrictive to imply convergence in the present case. In so far, it may be worth to study the relationship between Markov chains and approximation theory in more detail in order to discover alternate techniques for detecting convergence of the matrix-driven Jacobi-Perron algorithm. Notice that the approximants of the invariant measure x of (1) are entirely determined by the finite difference schemes (65)–(67). In so far, the analytic proof of the existence of the limits (68) is mainly of theoretical interest.


## 6. Numerical Example

Consider a queueing system consisting of a single server which is subject to breakdowns. The up time of the system is exponentially distributed with parameterα>0. The down time is assumed to be Erlang distributed with shape parameter K∈ℕ and rate parameter Kβ for some β>0 (implying that the mean down time is 1/β). Customers arrive according to a Poisson stream with rate λ>0. By κi>0, we denote the probability that an arriving customer joins the queue when i customers are present and the server is active. Analogously, θi≥0 shall denote the probability that an arriving customer joins the queue when i customers are present and the server is down. Notice that it could make sense to assume θi=0, that is, customers will not join the queue when the system is interrupted. When less than n customers are present, no service is provided. Otherwise, n customers are served at a time. The service time of each batch is assumed to be exponentially distributed with parameter μ>0. In addition, we assume that the arrival process, the service process, and the breakdown mechanism are stochastically independent. The model under consideration is an extension of the usual M/Mn,n/1 queue, which is described in more detail in [33].Remember that an Erlang distributed random variable may be interpreted as a phase-type distributed random variable [1] with distribution function Fx=1−αexpTx1 for x≥0, where
(74)T=−KβKβ−KβKβ⋱⋱−KβKβ−Kβandα=1,0,⋯,0.Erlang distributions are frequently used when modelling random variables with a small variance or even approximating deterministic behaviour. This is due to the fact that the variance of the Erlang distribution with parametersK and Kβ is 1/Kβ2, and this is indeed the smallest variance which can occur for a phase-type distribution with order K and expectation 1/β (see [34] for more details).Due to the lack-of-memory property of the exponential distribution and the independence assumptions, the system under consideration is characterized by a two-dimensional continuous-time Markov chainY=Ytt≥0 with Yt=Nt,Ut where
(i)
Nt is the number of customers in the system at time t(ii)
Ut=j∈1,⋯,K indicates that the server is down and the down time is in phase j at time t(iii)
Ut=K+1 indicates that the server is up at time t.Hence, the underlying state space isE=0,1,2,⋯×1,2,⋯,K+1, and we obtain a generator matrix of structure (20) with
(75)Qi,i+1=λθi0⋯000λθi⋱00⋮⋱⋱⋮⋮00⋯λθi00⋯⋯0λκii=0,1,2,⋯,Qi,i−n=0⋯00⋮⋱⋮⋮0⋯000⋯0μi=n,n+1,⋯,Qi,i=−λθi−KβKβ−λθi−KβKβ⋱⋱−λθi−KβKβα0⋯0−λκi−αi=0,1,⋯,n−1,Qi,i=−λθi−KβKβ−λθi−KβKβ⋱⋱−λθi−KβKβα0⋯0−λκi−α−μi=n,n+1,⋯.A sufficient condition forY to be regular and ergodic is
(76)limsupi+j→∞λ·maxθi,κi<nμ·βα+β.This statement can be deduced from Tweedie’s criterion ([35], Theorem 2.3). The condition is to find a nonnegative sequence ξi,j and an integer N>0 satisfying ξi,j⟶∞ for i+j⟶∞ and
(77)λθiξi,j−ξi+1,j+Kβξi,j−ξi,j+1≥ε,j=1,⋯,K,λκiξi,K+1−ξi+1,K+1+αξi,K+1−ξi,0+μξi,K+1−ξi−n,K+1≥ε,i≥n, i+j≥N, and some ε>0. A simple check shows that
(78)ξi,j=i+K+1−jK·nμα+βi=0,1,2,⋯;j=1,2,⋯,K+1,satisfies this set of inequalities as long as (76) holds.Next, we focus on the calculation of the following performance values:(i)
The expected stationary numberEN∞ of customers in the system is given by
(79)EN∞=∑i=0∞πi·ii⋮ii,where πi=πi,1,⋯,πi,K+1(ii)
The corresponding variance is calculated according toVARN∞=EN∞2−EN∞2 where
(80)EN∞2=∑i=0∞πi·i2i2⋮i2i2(iii)
The stationary probability that no service is provided despite≥n customers waiting is given by
(81)p~≔ℙU∞≤K∣N∞≥n=ℙN∞≥n,U∞≤KℙN∞≥n,where
(82)ℙN∞≥n,U∞≤K=∑i=n∞∑j=1Kπi,j=∑i=n∞πi·11⋮10,ℙN∞≥n=∑i=n∞∑j=1K+1πi,j=∑i=n∞πi·11⋮11(iv)
The stationary probability of the system being down is(83)ℙU∞≤K=∑i=0∞∑j=1Kπi,j=∑i=0∞πi·11⋮10.The latter probability can be computed exactly as1/β/1/α+1/β=α/α+β. Therefore, the numerical computation of this probability only helps us to check the implementation.For calculating all these values, we apply Algorithm2 with L=5 and
(84)fi=ii21n,n+1,∞i1n,n+1,∞i11ii21n,n+1,∞i1n,n+1,∞i11⋮⋮⋮ii21n,n+1,∞i1n,n+1,∞i11ii201n,n+1,∞i01,where 1Ax=1 for x∈A and 1Ax=0 for x∉A. After applying the algorithm, the variance can be deduced from the first two entries of the resulting vector. Similarly, the conditional probability p~=ℙU∞≤K∣N∞≥n can be obtained from the third and the forth entry.In Tables1, 2, and 3, we present some numerical results for n=3, K=10, μ=1.0, λ=1.0 or λ=4.0, κi=1/i/n+1, and θi=κi/10 for i=0,1,2,⋯, and various values of α and β.Table 1
EN∞ for λ=1.0 (a) and λ=4.0 (b), respectively.
(a)α\β210.50.10.51.8251.8551.9252.5370.11.8061.8131.8322.0960.051.8031.8071.8171.9850.011.8011.8021.8041.848(b)α\β210.50.10.54.4134.5314.7876.9120.14.3324.3584.4235.1920.054.3224.3354.3694.8440.014.3144.3164.3234.445Table 2
VARN∞ for λ=1.0 (a) and λ=4.0 (b), respectively.
(a)α\β210.50.10.51.8531.8921.9772.7980.11.8261.8341.8542.1320.051.8221.8261.8372.0140.011.8191.8201.8221.871(b)α\β210.50.10.56.9697.1557.54411.6360.16.8186.8546.9277.9410.056.8006.8176.8537.4620.016.7856.7896.7966.959Table 3
p~ for λ=1.0 (a) and λ=4.0 (b), respectively.
(a)α\β210.50.10.50.2040.3440.5230.8680.10.04880.09520.1810.5920.050.02500.05000.09940.4240.010.005110.01040.02160.130(b)α\β210.50.10.50.2020.3400.5120.8460.10.04830.09350.1750.5390.050.02480.04900.09580.3720.010.005060.01020.02080.107In Table4, we fix λ=μ=1.0 and α=β=0.1 and vary the coefficient of variation of the down time (i.e., CV=1/K) from 1 to 1/20.Table 4
EN∞, VARN∞, and p~ as a function of the coefficient of variation of the down time (i.e., CV=1/K).
K1351020EN∞2.2442.1362.1132.0962.087VARN∞2.4312.2002.1602.1322.118pu0.6310.6040.5980.5920.589The usualM/G/1 queue with service interruptions has been investigated in great detail by Gaver [36] and Federgruen and Green [37]. Their explicit formulas tell us that the mean number of customers in the system and the variance of the number of customers in the system are strongly affected by the arrival rate (i.e., λ), the ratio of the mean repair time and the mean service time (i.e., μ/β), and the availability of the system (i.e., ℙU∞=K+1=β/α+β). Obviously, the same is true for the state dependent bulk service queue with unreliable server (see Tables 1 and 2).The conditional probabilityp~ may be interpreted as a certain measure for the non-value added time of the server. The data in Table 3 suggest that for λ⟶∞, convergence is to α/α+β=ℙU∞≤K. This can be deduced from the law of total probability which states
(85)αα+β=ℙU∞≤K=p~·ℙN∞≥n+ℙU∞≤K∣N∞<n·1−ℙN∞≥n.SinceℙN∞≥n⟶1 for λ⟶∞, we conclude from (85) that p~⟶ℙU∞≤K=α/α+β as λ⟶∞. So it is not surprising that α and β are the main influencing factors, whereas the arrival rate λ has only little effect on p~ which may be explained as follows: Due to its state dependence, the arrival stream in our example is massively dampened when the number of customers grows. In other words, in spite of an excessive arrival rate (i.e., λ>μ), the system always reaches a comparatively moderate and balanced state.Finally, Table4 reveals that the coefficient of variation of the down time has a relatively small influence on the main characteristics of the system.It should be mentioned that our model allows the specification of more complex interarrival, service, and up and down times by fitting phase-type distributions to empirical data or by approximating general distributions by phase-type distributions. Since it is only our intention to demonstrate the practical suitability of our algorithm, we do not deepen the discussion of the numerical results.Obviously, a lot of other queueing situations can be taken into account within our model. In [17], a tandem queue was analysed in which
(i)
customers arrive at the first station according to a Markovian arrival process(ii)
customers served at the first station are transferred to the second station(iii)
customers served there leave the system(iv)
there arec1 and c2 servers at both stations, respectively(v)
at both stations, the service time is phase-type distributed.As pointed out in [17], the state space of the underlying queueing process consists of the number of waiting customers (at the first station and at the second station), the number of servers serving in the various phases of the phase-type distribution (at the first station and at the second station), and the environmental state for the Markovian arrival process. But by appropriately renumbering the states, it is seen that the process may be transformed into a quasi-birth-death process having a generator matrix of the form (20) with n=1.Next, consider a batch service at both stations with a maximum batch sizen at the second station. Then departures at the second station decrease the total number of customers in the system by n. Hence, a structure of type (20) with nonconstant block-dimensions comes into play. Since the state space remains finite, algorithms 1 and 2 are exact up to numerical errors. Similarly, many other queueing models may be solved in this manner.


## 7. Conclusion and Further Research

In this paper, we have suggested methods for computing invariant measures and stationary expectations for Markov chains with a band-block-structured transition probability matrix or generator matrix, respectively. We have presented a probabilistic derivation as well as a derivation based on generalized matrix-valued continued fractions. Although the latter derivation only works under some restrictive assumptions, we have gained some additional numerical and algebraic insights. From a practical point of view, we want to emphasize the benefits of Algorithm2. This method generalizes the established algorithm for computing stationary expectations for quasi-birth-death processes [7], and in situations where we are only interested in stationary expectations (what should be the case in most practical applications), it is extremely efficient with respect to memory requirement.As pointed out in the remarks concerning the Algorithms1 and 2, all considerations but the memory-efficient computation of stationary expectations can be easily extended to lower-block Hessenberg matrices. Unfortunately, for general lower-block Hessenberg matrices (and large truncation level N), the memory requirement will be crucial. Hence, further research could concern memory-efficient use of the lower-block Hessenberg transition structure (skip-free to the right). Other structures being of interest in practical applications are upper-block Hessenberg matrices or general block-band matrices with upper bandwith >1. Another interesting point concerns the choice of the truncation level N and the derivation of precise error bounds.



---

*Source: 4318906-2020-07-08.xml*