Abstract
For dealing numerically with the infinite-state-space Markov chains, a truncation of the state space is inevitable, that is, an approximation by a finite-state-space Markov chain has to be performed. In this paper, we consider level-dependent quasi-birth-death processes, and we focus on the computation of stationary expectations. In previous literature, efficient methods for computing approximations to these characteristics have been suggested and established. These methods rely on truncating the process at some level , and for , convergence of the approximation to the desired characteristic is guaranteed. This paper’s main goal is to quantify the speed of convergence. Under the assumption of an -modulated drift condition, we derive terms for a lower bound and an upper bound on stationary expectations which converge quickly to the same value and which can be efficiently computed.
1. Introduction
In many applications, discrete-time or continuous-time Markov chains are used to model real-life problems. For these models, characteristics of interest can be determined, enabling us to understand relationships between parameters and characteristics, solve optimization problems, and so on. When Markov models are used, the characterization of the dynamics by its probability transition matrix or generator matrix , respectively, is quite simple due to all theoretical problems being solved. However, in applications, the determination of interesting characteristics of the process is very important. In many models, there is no chance of finding explicit analytical representations of these characteristics, and hence, we have to use numerical calculations or simulation methods. In most practical situations, interesting characteristics refer to stationary probabilities or stationary expectations which can be derived from the invariant distribution. In queueing theory, elementary examples are given by moments of the stationary number of customers in a queueing system or queueing network, or by the stationary probability that the number of customers exceeds some threshold. Similarly, such characteristics are interesting in Markovian population models, epidemic models, etc. Formally, stationary expectations are given by , where is the state space of the Markov chain, is the invariant distribution, and is some cost or reward function.
The invariant distribution for an irreducible and positive recurrent Markov chain is characterized by the unique probability vector which solves a certain homogeneous system of linear equations. The number of states coincides with the number of equations and the number of variables. Unfortunately, in realistic models, we have a very large number of states or even infinitely many states. In these situations, numerically computing stationary probabilities or stationary expectations requires a truncation of the system of equations (which means a truncation of the state space), and this truncation results in inevitable errors.
Hence, numerical calculations only compute some approximation to the desired characteristic , and from both a mathematical and a practical point of views, bounds on the truncation error are interesting. In this paper, we will derive lower and upper bounds on which is equivalent to computing an approximation to and bounds on . Since all values are computed numerically, our results and methods provide a posteriori error estimates. Competing methods for finding a posteriori error estimates were given in [1–3], and an approach for finding a priori estimates can be found in [4, 5]. In Section 8, we will briefly compare these methods to the method suggested in this paper.
Although our theoretical results will be quite general, we will only discuss the efficient implementation of the computation of the bounds for level-dependent quasi-birth-death processes (LDQBDs). These processes are characterized by the block-tridiagonal structure of the transition probability matrix (discrete-time) or generator matrix (continuous time). This structure is very typical for a large class of queueing models, population and epidemic models, etc.
The rest of the paper is organized as follows: After introducing some notation in Section 2, we will present our theoretical results on lower and upper bounds on stationary expectations in Section 3. In Section 4, we will derive a method for computing these bounds efficiently for LDQBDs. Afterwards, we will apply our method to the elementary queue (Section 5), a variant of the queue (Section 6), and a popular retrial queueing model (Section 7). Whereas we have an explicit analytical representation of the invariant distribution and of stationary characteristics in the first example, numerical computations are the only chance to find stationary expectations in the latter two examples. Finally, we give a comparison of the suggested method to competing ones (Section 8).
2. Preliminaries
2.1. Notations and an Auxiliary Result
In order to introduce our goals formally, we review some well-known facts on the limit behaviour of Markov chains. In this context, and in the whole paper, we will use the following notations: (i) is a finite or infinite identity matrix. The dimension will become clear from the context(ii)1 denotes a finite or infinite (column) vector with all entries being 1. Alternatively, 1 refers to a function with constant value 1. The meaning will be clear from the context(iii) denotes the indicator function of set , that is, if , and otherwise. With slight abuse of notation, we use for Boolean expressions . If is true, is 1, and otherwise, is 0(iv)id is the identity function(v)We use the notation and .
Furthermore, we will use a probabilistic proof of invertibility frequently throughout the paper: If is finite and transient, is invertible. We give some details: A finite substochastic matrix is referred to as transient if (component-wise). This is true if and only if contains no stochastic submatrix for some . Hence, is transient if and only if, for all , there are with for , and .
Due to being finite, entails for all eigenvalues of P, and hence, converges, and by a standard argument, we obtain which means that is invertible.
2.2. Results on the Limit Behaviour of the Markov Chains
We briefly summarize some results on the limit behaviour of Markov chains since these results are the reason why it is important to be able to compute (approximations to) stationary expectations .
We consider Markov chains in discrete time or Markov chains in continuous time with countable state space . Throughout the paper, we assume that is irreducible and recurrent with transition probability matrix or that is nonexplosive, irreducible and recurrent with generator matrix . Then, there is an invariant measure , that is, a positive vector with or , respectively. In both cases, is unique up to constant multiples. For , we introduce the notation . If and converge absolutely with for functions , a fundamental result on the limit behaviour of Markov chains is in the sense of almost sure convergence (see [6], pp. 85-86, 203-209). If we have the stronger assumption of positive recurrence, the sum of the entries of is finite. Hence, by multiplying with an appropriate constant, we obtain the uniquely determined invariant distribution which still is an invariant measure, but additionally satisfies . Then, we can set and constant, and obtain almost surely if converges absolutely. If measures costs or rewards that are associated with the time spent in state , is the long-run average of costs or rewards per time unit. In particular, for , we can specify and is the long-run average of the proportion of time spent in set . For , is the long-run average of the proportion of time spent in state . Note that there is another interpretation of the entries of the invariant distribution : For irreducible and positive recurrent discrete-time Markov chains which are additionally aperiodic, converges to some random variable in distribution, and for nonexplosive, irreducible, and positive recurrent continuous-time Markov chains, converges to some random variable in distribution. Then, is the distribution of or (see any textbook on Markov chains, e.g., [6, 7]), that is, and . Hence, if converges absolutely for some function , we have or , respectively.
Due to these results on the limit behaviour, in many applications, all characteristics of interest can be deduced from for functions , where we will set in case of positive recurrence. Furthermore, by splitting functions into positive and negative part, we can assume for all without restriction. We remark that is vector valued, .
2.3. Notations for Block-Structured Markov Chains
Throughout this paper, let the state space be partitioned into finite and pairwise disjoint sets which will be referred to as levels. We introduce some notation: (i)We fix . will be chosen such that the -modulated drift condition (see below) holds on . Furthermore, throughout this paper, is some arbitrary, but fixed state within level (ii) denotes the invariant measure with . In case of positive recurrence, is the invariant distribution(iii)For , or , respectively, shall denote the transition probabilities or transition rates, respectively, for transitions from states to states .(iv)For , we write and for functions .
Later on, multidimensional functions will be useful for computing approximations to several stationary expectations at the same time. With this notation, we have , and we are interested in computing approximations to and use them for approximating , or we want to find an approximation to directly by other means. As in other approaches for finding error bounds on the approximation of stationary characteristics, we assume that an -modulated drift condition is satisfied, that is,
In scalar notation, (6) reads as for . -modulated drift conditions are a popular tool to prove positive recurrence and convergence of or . In our examples in Sections 5–7, we will use this tool, too. For details, we refer to Theorem 1 and Theorem 2 in [8]. Note that the idea of making use of -modulated drift conditions for this purpose is much older (see, e.g., [9–13], but some of these classical references formally require , whereas the results in [8] only require .
Finally, throughout the paper, we will require that
We will comment this structural requirement in Section 8.
3. Upper and Lower Bounds on : Theoretical Results
Our goal is to find a lower bound and an upper bound on such that (i)Both bounds converge to as the truncation level tends to (ii)Both bounds can be computed numerically in an efficient manner.
In this section, we focus on the “mathematical” criteria, that is, we derive exact terms for the lower and upper bounds in which both converge to . The latter goal of computational efficiency will be considered in Section 4 for the special case of QBDs.
First, we focus on the discrete-time setting. We use the notation for block-structured Markov chains as introduced above, and we begin with a representation of invariant measures and their approximations. For this purpose, we introduce the hitting times , , and .
Lemma 1. Let be an irreducible recurrent discrete-time Markov chain with block-structured transition probability matrix as in the preliminaries, let with and . Let and . (a) is well defined, and we have and all increase monotonically in . (b) exists for all with and exists with For the limits, we have for all and . (c) is stochastic and irreducible. Let be the invariant measure for subject to , and let for . Then is the unique invariant measure for subject to .(d)We have where .
Proof. We refer to the representations of , , , and as “probabilistic interpretations.” (a)The definition of implies Since is finite and transient, is invertible, and hence, is well defined for . Afterwards, the definition of for is explicit. We omit the proof of the probabilistic interpretation since such considerations are standard in the context of finding probabilistic interpretations for invariant measures (see, e.g., [7, 14]), and since these considerations are similar to those in the proof of Lemma 3.2b) or Lemma 4 below. Note that the probabilistic interpretation directly implies the monotonicity since increases monotonically. (b)More precisely, converges to almost surely. By monotone convergence, the probabilistic interpretation of and follows. Note that these expectations are finite due to recurrence. for and can be proved by means of the probabilistic interpretations. Again, we omit all further details.(c)The probabilistic interpretation of can be rewritten as Due to recurrence, we have almost surely, and stochasticity and irreducibility follow easily. Since is finite, it admits an invariant measure which becomes unique by requiring . Multiplying by immediately leads to for , and analogously, we find . Hence, is invariant for . Note that we can find the probabilistic interpretation by standard considerations (again, we omit details). This is the standard representation of invariant measures with (see, e.g., [7, 14]). (d)The statement is a direct consequence of and the representation of found in (c).Both and depend on for all . Our goal is to find bounds both on and which can be computed from , since these matrices only depend on the finite matrix . We begin with bounds on .
Lemma 2. Let all requirements of Lemma 1 hold and choose such that is an invertible diagonal matrix. (a)Such a matrix exists, and for two such matrices we have with an invertible diagonal matrix D.(b)For any such choice of , we have (c)Let and be defined as in (b). Then increases monotonically in N, decreases monotonically in , and both converge to for all .
Proof. We adapt ideas for proving classical bounds on quotients of entries of the invariant measure (see, e.g., [15]).
(a)The probabilistic interpretation of implies that is transient. Hence, is invertible and we could choose . If both and are invertible diagonal matrices, so is .(b)First, we note that and do not depend on the choice of since a left-hand multiplication by some diagonal matrix has no impact on the quotients under consideration. Now, we will directly construct a matrix with for all . With and , this works as follows: We set
Additionally, we set . Then for all, is the standard construction (see, e.g., [7, 14]) of the minimal subinvariant measure for the substochastic matrix subject to . Indeed,
for . Hence, is a diagonal matrix. Analogously, the rows of are subinvariant measures for . This latter matrix is finite and stochastic, and hence, every subinvariant measure is invariant. In particular, every row of is a constant multiple of , that is,
for all . Due to , we have by a trivial induction, and thus, . From , we obtain
for all . By setting , we obtain (due to ), and by setting , we obtain .(c)We have , implying , and by monotone convergence, we obtain
By setting or , respectively, the statement follows.
Lemma 2 gives a lower and an upper bound on all entries of ψK. It remains to find bounds on . A lower bound is a direct consequence of Lemma 1.
Lemma 3. Let all requirements of Lemma 1 hold, let , and let (a)We have .(b) increases monotonically in with limit .
As pointed out above, for finding an upper bound, we have to use some information on with or . This is done by using the -modulated drift condition (5). Note that the function is given, and the function solving the drift condition can often be found by “educationally guessing” even if there is no chance to find explicit analytical terms for or .
Lemma 4. Let all requirements of Lemma 1 hold, let the -modulated drift condition (5) hold, and let the structural requirement (7) be fulfilled. Furthermore, define . (a)For , we have .(b)Additionally, let . Then converges to .
Proof. (a) First, remember
where for implying
For , we have
and together, we obtain
Analogously, we find
Hence, it remains to show for all
. For this purpose, we first write
Let . By means of the tower property of conditional expectations, we find
Together, we obtain
for the entries of . Below, we will show for , and due to and , the desired statement follows.
For proving , we define for . Then for , and by induction, we obtain
for all and all , since transitions from have probability . Due to recurrence, is finite almost surely, implying almost surely, and by monotone convergence, , implying for all as desired.
(b) Since converges to , it suffices to prove that converges to . For this purpose, we write
Again, let . For , we have (due to the respective probabilistic interpretation, see proof of Lemma 1 for the interpretation of ), and since we assume that converges, we have convergence to for . Since the choice of is arbitrary, this is true for any .
By summarizing all previous results, we obtain our main result for the discrete-time setting where we will omit the dependency of on . Furthermore, we will phrase the result directly for both discrete-time and continuous-time setting. All previous results may be transferred into the context of continuous-time Markov chains easily by means of the embedded jump chain. Since similar considerations have been used in the literature frequently (see, e.g., [4, 16, 17]), we omit further details.
Theorem 5. Let be an irreducible and recurrent discrete-time Markov chain, or let be a nonexplosive, irreducible, and recurrent continuous-time Markov chain, respectively. Let the transition probability matrix or , respectively, be block-structured as introduced in the preliminaries. Let the drift condition (5) or (6), respectively, and the structural requirement (7) hold. Furthermore, let (i) be defined by and or (ii) be chosen such that or AT, respectively, is a diagonal matrix, let and for , and let and as row vectors(iii) or , resp., for (iv) or , resp., be component-wise finite for all (v)and (a) is uniquely defined for , and or T, resp., is invertible, in particular, there is a matrix A with the desired properties(b)We have (c)If , we have
Remark 6. Theorem 5 is phrased for Markov chains with an infinite state space, and the goal is to find approximations to which rely on the transition probabilities for transitions within a finite subset of the state space. Since even finite state spaces can be so large that an exact computation of cannot be performed, we might want to apply Theorem 5 to Markov chains with a finite state space in order to reduce the state space to a smaller one.
Indeed, for , , for , and , (a) and (b) in Theorem 5 remain true. Note that formally, we cannot allow since would become stochastic, and would not be invertible anymore. However, means that we have no reduction of the state space, and hence, this setting makes no sense. Part (c) becomes obsolete.
4. Efficiently Computing the Bounds for Quasi-Birth-Death Process (Discrete-Time Setting)
We turn now to developing methods for computing both bounds on efficiently (and simultaneously). As pointed out in the introduction, we focus on quasi-birth-death processes (QBDs). QBDs are characterized by or , resp., for with , that is, each jump changes the level at most by 1.
Level-independent QBDs were analysed by Neuts [18] by means of matrix-geometric methods. First algorithmic approaches for level-dependent QBDs (LDQBDs) are due to Bright and Taylor [17] and Hanschke [19]. The approach in [17] generalizes Neuts’ probabilistic interpretations of the matrices which arise in the matrix-geometric method, whereas the approach in [19] is motivated by the relationship between second-order vector-matrix difference equations and matrix-valued continued fractions. Remarkably, both methods are equivalent (up to suggestions regarding some initializations). More details and comparisons can be found in [20, 21]. All these methods intend to find (approximations to) an invariant measure or the invariant distribution. In [22], an algorithm was suggested which allows to compute stationary expectations directly. In the continuous-time setting, this method reads as follows: (i)Choose large , set and (ii)For , replace by and by (iii)Determine as an (approximative) solution to (iv)Return .
Since the memory requirement does not depend directly on (with , it depends on max{}), the truncation level can be chosen very large. Despite the possibility to choose large , results on the truncation error still remained desirable from both a mathematical and a practical points of view.
It is not difficult to prove that (with the notation of the present paper) the method from [22] computes the lower bound on for ( has no impact on the lower bound). The goal of this section is to generalize this method in such a way that any value is allowed, and that the upper bound will be computed, too. For means of conciseness, we focus on the continous-time setting.
Note that the QBD property directly implies for , that is, (7) holds. The drift condition (6) is met if and the definition of simplifies to
Most importantly, the computation of the matrices according to the system becomes much easier (the matrix-analytic methods in, e.g., [17, 19–21] use this fact). Many of the following identities are consequences of the probabilistic interpretation of the involved matrices, and we omit these proofs since the considerations are similar to those in [17] anyway. (i)Let be defined by and Then, the inverses in (39) exist, and we have (ii)Similarly, for defined by and we have (iii)The structure of allows to write .
If the level sizes are relatively small, the computational effort induced by (39) and (41) is acceptable. Then, it seems natural to compute all and , then all and finally and . However, in [22], it was pointed out that for , computing can be performed in a much more efficient way (in particular with respect to memory requirement) by using a Horner-type scheme. Here, we generalize this procedure slightly with respect to two issues: We consider an arbitrary , and we consider the simultaneous computation of and . Introduce
With this notation, we have
The recurrence schemes for , , and all start at , and are used in the sequel for computing the other values for . Similarly, the computation of and starts at , and then the values for can be computed. For finding , we only need the values for . Hence, for , is only used for finding , , and . Similarly, the other matrices are only used in a single step. Hence, there is no need to store these matrices. Instead, we suggest to use the recurrence scheme for all these matrices as an “update” procedure. In total, we suggest the following method. (i)Choose, set , , and .(ii)For , update (iii)Set , and .(iv)For , update (v)Set and determine such that is a diagonal matrix, define by , and define by .(vi)Return as a lower bound and as an upper bound.
Usually the matrices and can be generated when they are needed. Up to four of these matrices are needed at the same time, and we need memory for saving . In total, a small number of finite matrices have to be stored at the same time. In particular, if is bounded by , the memory requirement is , and this bound does not depend on or . Note that the “price” for this low memory requirement is that we only calculate the values for cost/reward functions which have to be specified before the computation procedure begins. Since we do not store the values or the matrices , adding any “new interesting” function requires a complete restart of the method.
A discussion of all numerical details of the algorithm is beyond the scope of this paper. For two specific aspects (avoiding ill-conditioned problems when computing and avoiding instabilities in the update step for ), we refer to the Appendix.
5. Application to the Queue
We first consider an example where we have explicit terms for . Of course, we would never use numerical methods for finding bounds on in such a situation, but the following considerations illustrate how the method works, and they show how sharp the upper bound can be. Numerical examples for situations in which we do not have an explicit analytical representation for will follow in Section 6 and Section 7.
Consider the queue, that is, we have and
is referred to as arrival rate, and is called service rate. is bounded, and therefore, nonexplosive. For any , is irreducible, and we will assume positive recurrence which is equivalent to .
Here, we have , and we will first consider choices of which allow choosing .
The invariant measure with is given by , and the -matrices shall solve for which reads as
This system of equations can be solved explicitly (e.g., by using standard results concerning linear difference equations with constant coefficients), and we find
Due to , we do not have to find the upper and lower bounds on , and instead, we focus on bounds on .
First, consider . As pointed out above, we have positive recurrence, and therefore, should be finite. From the explicit term, we obtain . To show the finiteness of by means of the drift criteria, set . Then, for all , we have
By Theorem 2 in [8] (or older references, see remark in the preliminaries), we obtain positive recurrence. Furthermore, the requirements of Theorem 5 are met. Therefore, we have , where
Hence, the lower bound converges to , and the upper bound coincides with which is the best possible result. It is clear that this can only happen if we have (instead of mere ≤ for all . Consider now where , that is, for and otherwise. Obviously, we can still use the same function , and with the same computation as before, we obtain
Let , and let . If we are interested in computing the bounds on , we combine the previous bounds, and the numerical algorithm would return
For any choice of , both bounds will converge to which is the stationary probability for set . As an easy choice, let for some . Then, we know that the stationary probability is given by , and for , we obtain the bounds
Next, let , that is, for all . Note that is the expected stationary number of customers in the system which is one of the most important performance measures in queueing theory. We set , and then we have
The exact value is , the computed lower bound is given by and the upper bound is given by
Hence again, we find the best possible upper bound . Appropriately combining the bounds on and leads to the upper and lower bounds on the expected stationary number of customers in the system.
At a first glance, the choice might be “difficult to guess,” but note that we have used a function of the form for dealing with , and it is quite natural to choose a function of the form for dealing with . Afterwards, it is not difficult to determine such that the -modulated drift condition is met. Nevertheless, in other situations, we will not be able to find such an optimal function . Therefore, we demonstrate next that the method still works if we use some function which is far from being optimal.
Let us directly consider , that is, , and . Then, is the stationary probability for set , and is the expected stationary number of customers in the system. Let us choose for and . Then for and all for some sufficiently large if we choose . Due to , the values will change. Now, we have
The solution to this system is given by
We omit an explicit representation of and (which is strictly larger than here), but we remark that and due to , this difference tends to 0 as . Note that the speed of convergence of to 0 is still exponentially fast, but slower than for the optimal choices of .
6. Application to a Variant of the Queue
Next, we consider an application of our method to a queueing model where we do not know the exact invariant distribution and have no chance but to use numerical methods. We consider a variant of the queue where arriving customers decide whether to join the queue or not depending on the number of customers they find in the queue. Precisely, we assume that (i)Customers arrive according to the Poisson process with intensity (ii)An arriving customer joins the queue with probability if there are other customers in the system upon the arrival(iii)The service time is phase-type distributed with parameters and , that is, the cumulative distribution function of a service time is given by and its expectation by (iv)There is one server.
This queueing system can be modelled as the Markov chain , where is two-dimensional: is the number of customers in the system at time , and for is the current service phase. By setting the service phase to 1 for , we obtain levels with and for . The generator matrix is given by
For , we obtain a level-independent QBD (the classical queue) for which an explicit analytical representation of the solution exists (see, e.g., [23]). For nonconstant , we obtain an LDQBD, and we have to use numerical methods. In what follows, we assume that which is equivalent to the existence of some with for all with some . Without restriction, we assume . Condition (63) guarantees that we have positive recurrence and that the stationary number of customers in the system has finite expectation as we will prove below. Hence, under condition (63), it makes sense to compute (i)the stationary probability that an arriving customer joins the queue where for all (ii)the stationary probability that at least customers are in the system where for and (iii)the expectation of the stationary number of customers in the system where for all .
In order to use our algorithmic method, we set additionally. For applying our method and for proving for under condition (63), we check some -modulated drift conditions. Set for , and then we have for and all choices of for which (component-wise). For we have , and we set for . Due to , we obtain for . From the special case , we obtain positive recurrence from Theorem 2 in [8], and by the same result, we have for For , we have , and we set for (note that is a nonnegative matrix, and hence, we have , implying for . Hence,
Choose such that
Then, we have for . The existence of such a value proves , that is, the stationary number of customers has finite expectation.
For applying our algorithmic method, we have to specify . In our numerical examples, we set , , . This allows to choose , and it is easy to derive that for . Furthermore, we consider the special case where the service times are Erlang-2-distributed, that is, and implying and
Then, (68) is guaranteed if . For and , this results in . Since we have to choose , we set . Note that simplifies to
In Table 1, some numerically computed results are listed. The parameters are chosen as specified above, and additionally, we set for . The results clearly illustrate the convergence of both bounds to the same value. For obtaining precise results for , the truncation level has to be chosen a bit larger than for obtaining quite precise bounds on and . This is not surprising since due to the fact that is the stationary probability for at least customers, the truncation level should be chosen significantly larger than 30. However, the results demonstrate that even for the computation of , the truncation level 55 yields very precise results.
7. Application to a Retrial Queueing Model
As a final example, we consider the queue with retrials which can be seen as some kind of “prototype” for LDQBDs. (i)Customers arrive according to the Poisson process with intensity (ii)The service times are independent and identically exponentially distributed with parameter (iii)There are servers(iv)The system capacity is , that is, customers can wait at the same time(v)Customers who cannot enter the system due to lack of waiting capacity enter the “orbit” of retrials(vi)Each customer in the orbit will retry to enter the system after a time which is exponentially distributed with parameter and independent of all other random variables(vii)Retrying customers who cannot enter the system due to lack of waiting capacity stay in the orbit of retrying customers.
Let , where is the number of customers in the orbit of retrials, and is the number of customers in the queue (including service). Then, is a continuous-time Markov chain with state space , where . Obviously, a transition with a positive rate will change the level at most by 1, and hence, we have a quasi-birth-death process with the following state transitions from state (): (i)Arrivals occur with rate . For , the arriving customer enters the queueing system, and the new state is . For , the arriving customer enters the orbit of retrials, and the new state is ()(ii)For , service completions occur with rate , and the new state is (iii)For and , successful retrials occur with rate , and the new state is .
In matrix notation, we have for , for , and for .
Note that retrial queues have been discussed intensively in the literature. We refer to the textbook [24] which gives an overview on retrial queueing models and computational methods for determining invariant distributions, stationary expectations, etc. Computational methods are very important in this context since even for , there are only explicit analytical representations of the invariant distribution if (see [24]).
Of course, there are many interesting characteristics such as the mean number of customers in the orbit and the mean number of customers in the queueing system. Our purpose is to demonstrate that our method provides reliable results for interesting characteristics, but we do not want to discuss the retrial queueing model in details. Therefore, we focus on a single characteristic, the stationary probability for arriving customers to find the queuing system full and being forced to join the orbit. If denotes the invariant distribution (if there is one), this probability is given by
The invariant distribution exists if and only if we have positive recurrence. Below, we will restate a proof for the fact that this is the case if . Note that is also a necessary requirement for stability/positive recurrence; we refer to [24] for more details.
We demonstrate how to prove that is sufficient for positive recurrence by means of drift criteria. The function which we use for this purpose can also be used for finding the upper bounds on and for any . In particular, by setting , we obtain bounds on
In any case (both for proving recurrence if and for applying our truncation method), we want to find a function such that we have for with some sufficiently large , where . We use the approach . For and , the entries of compute as
If we choose , we surely have for with some sufficiently large (due to ). Additionally, we have to take into account. Hence, we have to guarantee . Due to , such a choice with is possible. For example, we can choose and , since then we have
For this choice, we have for all .
As pointed out above, we have to guarantee that for . Hence, we choose
In total, we set where and , and we choose as above. Then, we can apply our algorithm which returns bounds on and where is the invariant measure with . By dividing lower and upper bounds appropriately, we find bounds on
In Tables 2 and 3, some numerically computed values for lower and upper bound are listed for different choices of the arrival rate and fixed values . In both cases, we see that lower and upper bounds converge to the same limit (which coincides with the value computed in Table 3.8 in the textbook [24] by other means). In Table 2, we have low traffic due to . Then, the blocking probability is quite low, and—more important for the evaluation of our method—we obtain precise results for low truncation levels. In Table 3, we have more traffic due to . Hence, the blocking probability is much larger, and we need higher truncation levels to obtain precise results.
The fact that we have to choose higher truncation levels in the situation of Table 3 is anything but surprising since due to more traffic, the Markov chain will assume states within higher levels much more often. However, it is important to remark that it is not clear how large should be chosen. The suggested method can be iteratively applied with increased truncation levels until a prescribed accuracy is achieved. Unfortunately, when increasing , we have to restart most of the calculations ( and do not change, but , , and do). However, for each such calculation, the memory requirement is the same since it does not depend on . In particular, in all situations considered in Tables 2 and 3, the memory requirement coincides.
As pointed out above, is an important performance measure, but there are more interesting characteristics for this model. Note that with the same choice for , we can find bounds on the stationary probability for any finite set of (finite or infinite) sets. In particular, if denotes the stationary probabilities for customers in the queueing system (waiting or in service), we could determine bounds on for simultaneously with little additional effort. Other performance measures (e.g., moments of the numbers of retrying customers) will require other choices for .
We conclude this discussion of the basic retrial queueing model by remarking that there are many extensions of the retrial queueing model which improve the applicability. For example, we could consider Markovian arrival processes, phase-type distributed service times, impatient customers who leave the queue and enter the orbit, impatient customers who leave the orbit, queueing networks where declined external arrivals join the orbit of retrying customers, etc. For all such models, the suggested algorithm allows to compute important performance measures precisely and efficiently. Note that such retrial queues provide good models for many problems in telecommunications and computer networks, e.g., in the design of wireless networks. However, a deep discussion of a realistic (and hence complex) a model justifies a separate publication and is far beyond the scope of this methodological paper.
8. Evaluation of the Method
8.1. The -Modulated Drift Condition
The major requirement for finding the bounds in Theorem 5 and for the resulting algorithmic procedure is the -modulated drift condition (5) or (6). We give some comments on this condition. (i)As pointed out above, -modulated drift conditions are very popular for proving convergence/finiteness of or . For , we remark that the existence of an -modulated drift condition is even equivalent to the convergence of . Hence, from a purely theoretical point of view, the requirement that an -modulated drift condition is satisfied is no restriction(ii)Most often, finding a function which satisfies the -modulated drift condition is the easiest way (or even the only way) to prove positive recurrence or to prove . If convergence of is not guaranteed, it makes no sense to compute approximations to (which are always finite) by numerical means. Hence, also from a practical point of view, finding a function which satisfies the appropriate -modulated drift condition does not require any additional effort in many situations(iii)Every function which can be used for proving can be used in Theorem 5. Note that this is not true for all truncation bounds which rely on -modulated drift conditions, see Section 8.6.
8.2. The Structural Requirement (7)
We give some comments on condition (7). Remember that is determined by the drift condition. If (7) is not met, we can always restructure the levels by setting and for , and set . Then (7) is fulfilled. If we have the QBD structure for the levels , this is preserved for the levels . Consequently, we have . This is unfortunate since (which becomes after the restructuring) should be quite small due to the operations involving matrices of dimension .
In total, from a purely mathematical point of view, can be guaranteed without restriction, and then the structural requirement (7) can be omitted. From a practical point of view, it is advantageous to keep and in particular, small, and in this context, we benefit from allowing . Then, the price is the requirement (7).
8.3. Asymptotic Properties of the Bounds
In order to compare our bounds to competing methods, we consider the asymptotic behaviour of the difference between upper and lower bounds.
In order to keep things simple and concise, we focus on the approximation of (where ) for positive recurrent continuous-time Markov chains, and we assume (which is a clearly restrictive condition), implying . The bounds on are given by and , and for , the difference behaves asymptotically like where we have used that converges to . In particular, if (e.g., if is an indicator function), we have the asymptotic expansion
These considerations clearly motivate to choose as small as possible.
8.4. Augmented Matrices versus Nonaugmented Matrices
Theorem 5 and the resulting method for computing bounds on or rely on considering the northwest-corner truncation which is substochastic. From , we obtain an approximation to , where is uniquely determined by and by satisfying up to the scalar equation corresponding to state .
There are many authors who prefer to “repair” the row sums of , that is, they deal with an augmented stochastic matrix . Due to finiteness and stochasticity, it is possible to determine an invariant distribution of directly. Under some mild constraints, converges to (see, e.g., [3]). The fact that admits an invariant distribution (in contrast to ) can be interpreted as an advantage of methods relying on augmented matrices. There are more theoretical issues why considering augmented matrices is popular, and additionally, the augmented matrices often allow meaningful probabilistic interpretations. A more detailed list of advantages can be found in [3] (for the continuous-time setting).
Despite these undeniable advantages, the approach relying on nonaugmented matrices was chosen in this paper. The simple reason is that converges (pointwise) monotonically to which makes lower bounds on trivial. As we will see in the next lines, this fact is more important than all the advantages of augmentation procedures.
8.5. Comparison to Other A Posteriori Bounds on Truncation Errors
The most recent results for truncation bounds on stationary expectations can be found in [1–3] where the authors consider continuous-time Markov chains. They consider northwest-corner truncations and their augmented versions (with row sums 0). With our notation and the invariant distribution of , Theorem 2.1 in [2] reads as
The authors require that (i)the Markov chain is irreducible and positive recurrent(ii)(iii) is constructed such that arises from truncation and augmentation and the other entries are arbitrary such that is conservative(iv)e.g., all other entries are 0(v) satisfies the -modulated drift condition (vi) is arbitrary(vii), where with the transition function arising from .
The right-hand side in (86) has some similarities to the asymptotic expansion (85). Asymptotically, the entries of behave like . Hence, (86) mainly differs by the factor and the summand . Furthermore, note that the right-hand side of (86) provides a bound on the difference between and its approximation. Hence, “” is larger by factor 2.
Finally, it is important to remark that cannot be determined by means of numerical computations. Therefore, it is suggested in [2] to use that increases monotonically in with limit . Hence, the entries of can be used to find computable bounds on , and thus, on .
Replacing in this manner has two effects: (i)There is no evidence that the bound in (86) is sharp. Even if it was, the bound which results from this replacement cannot be sharp anymore(ii)Evaluating the computable bound which arises from (86) requires computing at least some rows of which can be costly in terms of computational effort.
Together, we see that the algorithm suggested in this paper will provide slightly better bounds which are significantly easier to compute at the same time. Finally, in [1–3] is clearly restrictive. In particular, this requirement prevents us from using indicator functions in order to obtain stationary probabilities.
8.6. Comparison to A Priori Bounds on Truncation Errors
Theorem 5 provides a posteriori truncation bounds (similar to (86) from [2]), since both the approximation to and the bounds are computed numerically (simultaneously).
A priori bounds do not rely on numerical computations. Hence, before any computation scheme starts, a priori bounds would allow to determine the truncation level in such a way that the (relative) error of the approximation is smaller than a prescribed bound.
An approach in this direction was introduced in [5] and generalized in [4] where a method was developed which determines a level such that for some prescribed and some scalar function . The method works as follows: (i)Choosesome function and determine (continuous-time setting) .(ii)Set and .(iii)Set and .(iv)If is finite, we have . In particular, we can find such that and obtain (88).(v)If is not finite, we have to change our choice of .
At a first glance, this method provides a priori error bounds, but there are some problems. In the original papers [4, 5], it was already pointed out that it is not easy to find appropriate functions , and that very often, the bounds are quite pessimistic. In order to illustrate these effects, we consider the queue from Section 5 and set . Whereas Theorem 5 can be applied with for an appropriate constant , this choice will fail here: We obtain and for , that is, is constant for . In particular, if is very small (and hence, is large), we have . Note that the choice of is unimportant since it cancels out.
The next intuitive choice is , resulting in and . If , we have , , and
For , we obtain where . Hence, the method guarantees that whereas the exact value is . For and , the latter term converges to 1 much faster than . For example, for , we have , and hence, .
Summarizing this brief numerical example, the most intuitive choice for does not work, and the next choice provides very pessimistic estimates. It might be possible to find functions which provide sharp estimates, but it does not seem likely that such functions can be “guessed.”
There is an even more serious problem when interpreting the results from [4, 5] as a priori estimates for the truncation error: The numerator in (88) refers to the exact invariant measure which cannot be computed by numerical means (up to some trivial exceptions like the queue). For honest a priori bounds, we should be interested in guaranteeing with a numerically calculated approximation to . Unfortunately, the methods in [4, 5] cannot be modified in this direction. One might argue that this slight modification is negligible in view of the pessimistic bounds which are obtained. However, if we guessed an optimal function by accident, we would not have a bound on anymore. In total, the results from [4, 5] may be exploited to find theoretical estimates on stationary probabilities or expectations. With respect to numerical means, they can only be used for rough preliminary considerations (how large should the truncation level be chosen) before applying the method suggested in this paper. Even that might be hard due to the increased difficulty to find appropriate functions (in comparison to finding appropriate functions for applying Theorem 5).
9. Conclusion and Further Research
In this paper, we have developed an exact and efficient method for numerically computing lower and upper bounds on where is an invariant measure for an irreducible recurrent discrete-time or continuous-time QBD and is a multidimensional nonnegative cost/reward function. By combining the lower and upper bounds appropriately, we obtain lower and upper bounds on in the case of positive recurrence, where is the unique invariant distribution and is an arbitrary cost/reward function. The term can be interpreted as the stationary expectation of the Markov chain measured by , or as the long-run average of cost/rewards measured by . If is chosen to be an indicator function, we obtain stationary probabilities. The computation of in Section 6 illustrates that even very small probabilities can be computed with high relative precision.
The algorithm was developed as an extension of previous matrix-analytic methods which also dealt algorithmically with level-dependent QBDs, but for which no error bounds were known up to now. Due to the fact that we compute lower and upper error bounds, we obtain automatically an approximation to and (a posteriori) error bounds. In our examples, we have seen that the difference between lower and upper bounds converges to 0 very fast, that is, we obtain very precise results. For optimal choices of , we have observed that the upper bound coincides with the true value of .
Obvious directions of future publications concern detailed applications of our algorithm to specific Markov models, e.g., we could consider enhanced retrial queuing models with Markovian arrival processes, phase-type distributed service times, impatient customers, etc., but there is a huge variety of other practical processes which are modelled by infinite LDQBDs. If we focus on the analysis of a specific model, we should spend effort into finding very good functions which may result in obtaining very precise results for small truncation levels.
In this context, it may be worth to find not only an upper truncation level but also a lower one. Then, the numerical computation would focus on a very small region of the state space, and the method would become even more efficient.
Note that the main results of this paper are not restricted to QBDs. Theorem 5 only requires that no transitions from to occur with positive probability/rate if . Hence, these results also apply if or , respectively, has upper block-Hessenberg structure (also referred to as the block- structure). We have focused on QBDs here since this specialized structure allows a very efficient computation of lower and upper bounds. Future research could deal with developing efficient algorithms for more general structures.
We remark that some of the key results also transfer to nondiscrete state spaces: Let , and let the levels be petite in the sense of [10] (this property is satisfied by compact sets for many Markov chains). Then becomes a (sub-Markovian) kernel, and does so, too. The multiplication is defined straightforward by . We still have , where solves , and where . Under the structural assumptions (no transitions from to for ), the results on the bounds on remain unchanged. However, for finding by means of numerical calculations (and for finding a lower and an upper bound on ), it seems that some sort of discretization is inevitable.
We conclude with a remark concerning the relationship to nonprobabilistic topics: As mentioned above, for , the basic algorithm for computing approximations to the subvectors of the invariant measure can be deduced from general considerations concerning (matrix-valued) continued fractions. Hence, in some sense, our computation of bounds on is related to the speed-of-convergence statements for matrix-valued continued fractions. Note that the probabilistic interpretation of continued fractions (arising in nonprobabilistic contexts) and their convergents has led to new convergence criteria for continued fractions and their generalizations (see [25, 26]), and with further research, considerations similar to those in this paper might provide speed-of-convergence results.
Appendix
A. Alternative Algebraic Proof of the Upper Bound on
Additionally, we give a less stochastic, more algebraic proof of the key result, that is, the upper bound in Lemma 4, where we concentrate on discrete setting. Let be blocks with transition probabilities where the index (i) corresponds to states within (ii) corresponds to states within (iii) corresponds to states within (iv) corresponds to states within .
Then, the structural requirement (7) on the transitions reads as and . The -modulated drift condition (5) can be rewritten as
Let us set where is stochastic, and the property of irreducibility is inherited to from . Let the blocks of an invariant measure for be denoted by . Then and is an invariant measure for . These formulas can be proven in a purely algebraic manner although they have a stochastic interpretation: If the original Markov chain with transition probability matrix is only observed at return times to set , we obtain a new Markov chain with transition probability matrix (see [27]). The above representation of allows to write and due to , has the same meaning as before. As pointed out above, we do not want to give alternative proofs for all results here. Instead, we focus on the upper bound on given in Lemma 4.
From (A.1), we find and with (A.2), we obtain
Therefore, we can write
In the next lines, we will prove that the middle-term matrix is which obviously yields the desired statement . We have where we used the assumptions .
B. Remarks Concerning Avoidance of Numerical Problems
B.1. On Finding
In the discrete-time setting of Theorem 5, the matrix converges to a stochastic matrix as . Hence, for relatively small truncation levels , we can choose and compute this matrix inverse by standard methods. For large (when we want to obtain precise results), the inversion of is ill-conditioned. This problem can be avoided easily with a little additional effort: Remember that for an appropriate choice of , the th row of is the minimal subinvariant measure for the transient substochastic matrix subject to . Hence, we fix and solve the inhomogeneous system of equations , where . With and the row vectors , this system reads as .
As , converges to an irreducible stochastic matrix, and hence, converges to a transient substochastic matrix. Hence, is invertible, and this new system of equations is not susceptible to numerical problems. Of course, we have to find all rows of A separately, that is, we solve systems of linear equations of the dimension . However, if is relatively small, this effort can be neglected in comparison to the other computational steps. These considerations directly transfer to the continuous-time setting.
B.2. Exploiting the GTH Advantage in the Update Step for the Matrix
For many systems of linear equations arising in the context of the Markov chains, the coefficient matrix satisfies some row sum conditions which is preserved during all steps of solving the linear system of equations. This fact can be used for implementing these steps in a way which is not affected by numerical instability. Such procedures were introduced by Grassmann et al. [28] and therefore referred to as GTH advantage. In the present paper, this advantage should be used in the update step for the matrix . We give some more details.
First, note that the probabilistic interpretation guarantees that all entries of are nonnegative for all . Hence, computing the diagonal entries of results in computing differences of positive numbers. Furthermore, when performing any elimination procedure, the updates of the diagonal entries result in computing differences of positive numbers again. All other operations only compute sums of nonnegative numbers (which cannot cause numerical problems). For avoiding the need of subtractions, we remark that the row sums of and of add up to 0 (this statement can be deduced from the probabilistic interpretation). Let be the column vector of row sums of . Then the -dimensional matrix () has row sums of 0 where only the diagonal entries can be positive. An important finding on the application of (scaled) Gaussian elimination to systems with coefficient matrices of this type guarantees that these properties are preserved in each step of the elimination procedure. Hence, the pivot element in each step (which are the entries which might be affected by numerical instabilities) can be “corrected” by computing the negative sum of the other entries in the same rows. We demonstrate this by means of a simple example. (i)Let We want to compute .(ii)First, write
By performing the steps of Gaussian elimination, we convert the left-hand side to the unity matrix, and the right-hand side becomes the desired inverse matrix. The column in between acts as a control column. Before starting the Gaussian elimination, we can correct the first pivot element 4 by recomputing it as −(−2 −1 −1). (iv)The first row is scaled by . Afterwards a constant multiple (with factor 2) of the first row is added to the second row. This results in
The submatrix still has row sums of 0. The upper left entry 2 can be corrected by replacing it by before applying the next step. (v)After the next step, we obtain
Again, has a row sum of 0 which can be used for correcting the pivot element which is used for the last step of the elimination procedure.
Due to the need of using the control vector , and due to computing the inverse first instead of solving a matrix-valued system of linear equations directly, we have slightly more computational effort. The effects which are caused by instability justify this additional effort.
As an example for instability, consider the update procedure of in the retrial queueing model in Section 7. In fact, note that for this specific model, we could even simplify our algorithmic approach since we are able to determine explicitly. Hence, there is no need to compute recursively. This explicit representation is given by which can be proven by induction. The key ideas for the induction step are the recurrence relation and the observation that the last column of is given by ; we omit further details.
Instead, we use this explicit representation to evaluate numerically computed values of and and to demonstrate the effect of the GTH advantage. Let , , and . If we apply some ordinary elimination procedures to compute the values of (and give out the values of as a test), and if we use double precision in C++, the first four decimal places of the entries of are “correct” for . The first four entries in the last column of are −0.8501 instead of −0.8500, the first four entries in the last column of are −0.9003 instead of −0.9000, and the first four entries in the last column of are −0.9517 instead of −0.9500. Obviously, the error of these values increases drastically. The corresponding entries in have values between −8.9 and −9.4 (instead of −1.15), and they differ. For larger indices, some of these entries assume even positive values. Unsurprisingly, the resulting values for lower and upper bound on are completely wrong. When using the GTH advantage, this effect does not occur, and all matrices coincide with the values which can be obtained from the explicit representation.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The author acknowledges support by Open Access Publishing Fund of Clausthal University of Technology.