Two approaches to the construction of perturbation bounds for continuous-time Markov chains

The paper is largely of a review nature. It considers two main methods used to study stability and obtain appropriate quantitative estimates of perturbations of (inhomogeneous) Markov chains with continuous time and a finite or countable state space. An approach is described to the construction of perturbation estimates for the main five classes of such chains associated with queuing models. Several specific models are considered for which the limit characteristics and perturbation bounds for admissible"perturbed"processes are calculated.


Introduction
In the paper, some topics are considered that are related with the stability of both homogeneous and non-homogeneous continuous-time Markov chains with respect to the perturbation of their intensities (infinitesimal characteristics). It is assumed that the evolution of the system under consideration is described by a Markov chain with the known state space and it is the infinitesimal matrix that is given inexactly. Different classes of admissible perturbations can be considered. The "perturbed" infinitesimal matrix can be arbitrary, and the small deviation of its norm from that of the original matrix is assumed or it can be assumed that the structure of the infinitesimal matrix is known and only its elements are "perturbed" within the same structure. Below we will give a detailed description of these cases. In some papers it is assumed that the perturbations have a special form, for example, are expanded in a power series of a small parameter. This assumption seems to be too restrictive and unrealistic.
The study of stability of characteristics of stochastic models has been actively developing since the 1970-s [19,43,71]. At that time V. M. Zolotarev proposed an approach within which limit theorems of probability theory were treated as stability theorems. This approach was cultivated in the works of V. M. Zolotarev, V. V. Kalashnikov, V. M. Kruglov, V. Yu. Korolev and their colleagues in the framework of international seminars on stability problems for stochastic models founded by V. M. Zolotarev (see the series of the proceedings of the seminar published as Springer Lecture Notes starting from [20] or as issues of the Journal of Mathematical Sciences). This approach proved to be very productive for the study of random sums in queueing theory, renewal theory and theory of branching processes [15].
Since 1980-s the problems related to the estimation of stability of Markov chains with respect to perturbations of their characteristics have been thoroughly studied by N. V. Kartashov for homogeneous discrete-time chains with general state space and, in parallel, by A. I. Zeifman for inhomogeneous continuous-time chains within the seminar mentioned above, see [21,22,47]. In particular, a general approach for inhomogeneous continuous-time chains was developed in [47]. In that paper both uniform and strong cases were considered. Later birth-death processes were considered in [48] and general properties and estimates for inhomogeneous finite chains were considered in [49]. The paper [51] was specially devoted to estimates for general birth-death processes with the queueing system M t |M t |N considered as an example. It should be mentioned that these papers were not noticed in western papers. For example, in [3] it was stated that there were NO papers on stability of the (simplest stationary!!) system M |M |1. For the first time we used the term 'perturbation bounds' instead of 'stability' in the paper [57] on the referee's prompt. The same situation takes place with the Kartashov's papers cited above. The methods proposed in those papers seem to be used by most authors of subsequent studies in estimation of perturbations of discrete-time chains. Possibly, poor acquaintance with the early papers of Kartashov and Zeifman can be explained by the differences in terminology mentioned above: in the original (and foundational) papers the term 'stability' was used (in the proceedings of the seminar with the consonant appellation 'Stability Problems for Stochastic Models').
The present paper deals only with continuous-time chains, therefore the subsequent remarks mainly regard such a case.
Note that to obtain explicit and exact estimates of perturbation bounds of a chain, it is required to have estimates of the rate of convergence of the chain to its limit characteristics in the form of explicit inequalities. Moreover, the sharper convergence rate estimates, the more accurate perturbation bounds. These bounds can be more easily obtained for finite homogeneous Markov chains. Therefore, most publications concern just this situation, see, e. g., [2,7,30,31,32,33]. As this is so, two main approaches can be highlighted.
The first of them can be used for the case of weak ergodicity of a chain in the uniform operator topology. The first bounds in this direction were obtained in [47]. The principal progress related to the replacement of the constant S with log S in the bound was implemented in [30] and continued in Mitrophanov's papers [31,32,33] for the case of homogeneous chains and then in [55,57] and in the subsequent papers of these authors for the inhomogeneous chains. The contemporary state of affairs in this field and new applied problems related to the link between convergence rate and perturbation bounds in the 'uniform' case were described in [34]. In some recent papers uniform perturbation bounds of homogeneous Markov chains were studied by the techniques of stochastic differential equations, see for instance [42] and the references therein.
The second approach is used in the case where the uniform ergodicity is not assured, that is typical for the processes most interesting from the practical viewpoint. For example, birth-death processes used for modeling queueing systems, and real processes in biology, chemistry, physics, as a rule, are not uniformly ergodic.
Following the ideas of N. V. Kartashov (see a detailed description in [23]), most authors use the probability methods to study ergodicity and perturbation bounds of stationary chains (with a finite, countable or general state space) in various norms [3,12,35]. For a wide class of (mainly) stationary discrete-time chains a close approach was considered in [29] and more recent papers [1,18,25,26,28,36,38,44,45,46,70].
In the works of the authors of the present paper perturbation bounds for non-stationary finite or infinite continuous-time chains were studied by other methods.
The first papers dealing with non-statonary queueing models appeared in the 1970-s (see [14,16], and the more recent paper [27]). Moreover, as far back as in [13] it was noted that it is principally possible to use the logarithmic matrix norm for the study of convergence rate of continuous-time Markov chains. The corresponding general approach employing the theory of differential equations in Banach spaces was developed in a series of papers by the authors of the present paper, see a detailed description in [59,60]. In [47] (also see [48,49]) a method for the study of perturbation bounds for the vector of state probabilities of a continuous-time Markov chain with respect to the perturbations of infinitesimal characteristics of the chain in the total variation norm (l 1 -norm) was proposed. The paper [51] contained a detailed study of the stability of essentially non-stationary birth-death processes with respect to conditionally small perturbations. Convergence rate estimates in terms of weight norms and hence, the corresponding bounds for new classes of Markov chains were considered in [61,66,41,68].
In the present paper both approaches are considered as well as the classes of inhomogeneous Markov chains for which at least one of these approaches yields reasonable perturbation bounds for basic probability characteristics.
The paper is organized as follows. In Section 2 basic notions and preliminary results are introduced. In Section 3 general theorems on perturbation bounds are considered. Section 4 contains convergence rate estimates and perturbation bounds for basic classes of the chains under consideration. Finally, in Section 5 some special queueing models are studied.

Basic notions and preliminaries
Let X = X(t), t ≥ 0, be, in general, inhomogeneous continuous-time Markov chain with a finite or countable state space E S = 0, 1, . . . , S, S ≤ ∞. The transition probabilities for X = X(t) will be denoted p ij (s, t) = Pr {X(t) = j |X(s) = i }, i, j ≥ 0, 0 ≤ s ≤ t. Let p i (t) = Pr {X(t) = i} be the state probabilities of the chain and p(t) = (p 0 (t), p 1 (t), . . . ) T be the corresponding vector of state probabilities. In what follows it is assumed that As usual, we assume that if a chain is inhomogeneous, then all the infinitesimal characteristics (intensity functions) Further, to provide the possibility to obtain more evident estimates we will assume that |a ii (t)| ≤ L < ∞ for almost all t ≥ 0.
Then the state probabilities satisfy the forward Kolmogorov system dp dt where A(t) = Q T (t), and Q(t) is the infinitesimal matrix of the process. Let · be the usual l 1 -norm, i.e. x = |x i |, and B = sup Then for almost all t ≥ 0, and we can apply all results of [5] to equation (3) in the space l 1 . Let p 0 (t) = 1− i≥1 p i (t). Then from (3) we obtain the following equation (for a detailed discussion, see, e. g., [17,52]): where f (t) = (a 10 , a 20 , · · · ) T , ByX =X(t) we will denote the 'perturbed' Markov chain with the same state space, state probabilitiesp i (t), transposed infinitesimal matrix i,j=0 and so on, and the 'perturbations' themselves, that is, the differences between the corresponding 'perturbed' and original characteristics will be denoted byâ ij (t),Â(t). Let Now briefly describe the main classes of the chains under consideration. The details concerning the first four classes can be found in [65,66].
. This is an inhomogeneous birth-death process (BDP) with the intensities λ i (t) (of birth) and µ i+1 (t) (of death) correspondingly.
. This chain describes, for instance, the number of customers in a queueing system in which the customers arrive in groups, but are served one by one (in this case a k (t) is the arrival intensity of a group of k customers, and µ i (t) is the service intensity of the ith customer). The simplest models of this type were considered in [37], also see [65,66].
Class III. Let a ij (t) = 0 for i > j + 1, a i,i+k (t) = b k (t), k ≥ 1, and a i+1,i (t) = λ i (t). This situation occurs in modelling queueing systems with arrivals of single customers and group service.
. This process appears in the description of a system with group arrival and group service, for earlier studies see [39,40,61].
Class V. Consider a Markov chain with 'catastrophes' used for modelling of some queueing systems, see, e. g., [10,11,6,24,69,57]. Here the intensities have a general form whereas a single (although substantial) restriction consists in that the zero state is attainable from any other state and the corresponding intensities q k,0 (t) = a 0,k (t) for k ≥ 1 are called the intensities of catastrophes. Now consider the following example illustrating some specific features of the problem under consideration.
Example [57]. Consider a homogeneous BDP (class I) with the intensities λ k (t) = 1, µ k (t) = 4 for all t, k and denote by A the corresponding transposed intensity matrix. Then, as is known (see, e. g., [50]), the BDP is strongly ergodic and stable in the corresponding norm. On the other hand, take a perturbed process with the transposed infinitesimal matrixĀ = A+Â, where â 00 = −ε,â k0 = ε k(k+1) for k ≥ 1 andâ ij = 0 for the other i, j. Then the perturbed Markov chainX(t) (describing the 'M|M|c queue with mass arrivals when empty', see [4,24,69]) is not ergodic, since from the condition Ap = 0 it follows that the coordinates of the stationary distribution (if it exists) must satisfy the condition 4p k+1 =p k +p 0 ε k+1 ≥p 0 ε k+1 , which is impossible.
As it has already been noted, the (upper) bounds of perturbations are closely connected with the (correspondingly, upper) estimates for the convergence rate (also see the two next sections). On the other side, it is also possible to construct important lower estimates of the rate of convergence providing that the influence of the initial conditions cannot fade too rapidly, see [67]. It turns out that it is principally impossible to construct lower bounds for perturbations. Indeed, if we consider the same BDP and as a perturbed BDP choose a BDP with the intensitiesλ k (t) = 1 + ε,μ k (t) = 4(1 + ε), then the stationary distribution for the perturbed process will be the same as for the original BDP for any positive ε.

General theorems concerning perturbation bounds
First consider uniform bounds.
Theorem 1. Let the Markov chain X(t) be exponentially weakly ergodic, that is, for any initial conditions p * (s) ∈ Ω, p * * (s) ∈ Ω and any s ≥ 0, t ≥ s there holds the inequality Then, for the perturbations small enough (Â(t) ≤ ε for almost all t ≥ 0), the perturbed chainX(t) is also exponentially weakly ergodic and the following perturbation bound takes place: For the proof we will use the approach proposed in [30] and modified in [55] for the inhomogeneous case, also see [57]. Let Then moreover, whence, as t → ∞, we obtain (7).
If under the conditions of theorem 1 the Markov chain X(t) has a finite state space, then both Markov chains X(t) andX(t) have limit expectations and Now turn to the weighted bounds. Here we use the approach proposed in [47], also see the detailed description in [59,60]. Let Below we will assume that the following conditions hold: for almost all t ≥ 0.
Theorem 2. If a Markov chain X(t) is 1D-exponentially weakly ergodic, thenX(t) is also 1D-exponentially weakly ergodic and the following perturbation estimate in the 1D-norm holds: Proof. The detailed consideration can be found in [60]. Here we only outline the scheme of reasoning. Let V (t, s) andV (t, s) be the Cauchy operators for equation (4) and for the corresponding 'perturbed' equation, respectively. Then for all t ≥ s ≥ 0. Then, rewriting equation (4) as after some algebra we obtain the following inequality in the 1D-norm: On the other hand, z(t) 1D ≤ M e −at z(0) 1D + M a f, for any 0 ≤ s ≤ t. Hence, under any initial condition z(0) ∈ l 1D we obtain the following inequalities for the 1D-norm: So, the first assertion of the theorem is proved. Then second assertion follows from the inequality z 1E ≤ W −1 z 1D (see, e. g., [52]) and estimate (22) Remark 1. A number of consequences of this statement can be formulated, for example, this follows from The respective perturbation bounds can be formulated for strongly ergodic (for instance, homogeneous) Markov chains, see [60].

Remark 2.
As it was shown in [60], the bounds presented in theorem 2 and its corollaries are sufficiently sharp.

Convergence rate estimates and perturbation bounds for main classes
For Markov chains of classes I -IV an important role is played by the matrix B * * (t) = DB(t)D −1 . To begin with, write out this matrix for each of these classes.
For class I this matrix has the form in the case of a countable state space (S = ∞); in the case of a finite state space (S < ∞).
For class II this matrix has the form in the case of a countable state space (S = ∞); in the case of a finite state space (S < ∞).
For class III this matrix has the form in the case of a countable state space (S = ∞); in the case of a finite state space (S < ∞).
In the proofs of the following theorems we use the notion of the logarithmic norm of a linear operator function and related estimates of the norm of the Cauchy operator of a linear differential equation. The corresponding results were described in detail in our preceding works, see [17,52,66]. Here we restrict ourselves only to the necessary minimum.
Recall that the logarithmic norm of an operator function B * * (t) is defined as the number Let V (t, s) = V (t)V −1 (s) be the Cauchy operator of the differential equation Then the estimate holds. Moreover, if for each t ≥ 0 B * * (t) maps l 1 into itself, then the logarithmic norm can be calculated by the formula Now let Also note that if in classes II-IV the intensities a k (t) and b k (t) do not increase in k for each t, then in all the cases the matrix B * * (t) is essentially nonnegative (that is, its non-diagonal elements are nonnegative), then in (33) and (34) the signs of the absolute value can be omitted.
The following statement ([66, theorem 1]) is given here for convenience. hold. Then the Markov chain X(t) is weakly ergodic and for any initial condition s ≥ 0, w(s) and for all t ≥ s there holds the estimate Now let instead of (35), for all 0 ≤ s ≤ t a stronger condition holds.  (35) and (37) are equivalent.
Theorem 5. Let the conditions of theorem 4 hold. Then the Markov chain X(t) is 1D-exponentially weakly ergodic, under perturbations small enough (16) the perturbed chainX(t) is also 1D-exponentially weakly ergodic and perturbation bound (38) in the 1D-norm holds. If, moreover, W = inf i≥1 d i i > 0, then both chains X(t) andX(t) have limit expectations and estimate (18) holds for the perturbation of the mathematical expectation.
To obtain perturbation bounds in the natural norm it suffices to use inequality (24) mentioned above.

Corollary 2.
Under the conditions of theorem 5 the following perturbation bound in the natural l 1 -(total variation) norm holds: Note that it is convenient to use the results formulated above for the construction of perturbation bounds for Markov chains of the first four classes, see, e. g., [51,60,61,66].
For chains of the fifth class, as a rule, it is convenient to use the approach based on uniform bounds as will be shown below. These models were considered, e. g., in [57,62,63]. Let Theorem 6. Let the intensities of catastrophes be essential, that is Then the chain X (t) is weakly ergodic in the uniform operator topology and for any initial conditions p * (0) , p * * (0) and any 0 ≤ s ≤ t there holds the following convergence rate estimate: To prove this theorem we will use the same technique as in [57]. Rewrite the direct Kolmogorov system (3) in the form Here The solution to this equation can be written as where U * (t, s) is the Cauchy operator of the differential equation Note that the matrix A * (t) is essentially nonnegative for all t ≥ 0. Its logarithmic norm is equal to and hence, hold. Then the chain X (t) is weakly exponentially ergodic in the uniform operator topology, and if the perturbations are small enough, that is, Â (t) ≤ ε for almost all t ≥ 0, then the perturbed chainX(t) is also exponentially weakly ergodic and the perturbation bound (7) holds with c = c * , b = b * .
Here, to compare both approaches, we will mostly deal with the queueing system M t |M t |N |N with losses and 1-periodic intensities. In the preceding papers on this model, other problems were considered. For example, in [8] the asymptotics of the rate of convergence to the stationary mode as N → ∞, was studied, whereas the paper [9] dealt with the asymptotics of the convergence parameter under various limit relations between the intensities and the dimensionality of the model. In [54,56] perturbation bounds were considered under additional assumptions.
Let N ≥ 1 be the number of servers in the system. Assume that the customers arrival intensity λ(t) and the service intensity of a server µ(t) are 1-periodic nonnegative functions integrable on the interval [0, 1]. Then the number of customers in the system (queue length) X(t) is a finite Markov chain of class I, that is, a BDP with the intensities λ k−1 (t) = λ(t), µ k (t) = kµ(t) for k = 1, . . . , N .
It should be especially noted that the process X(t) is weakly ergodic (obviously, exponentially and uniformly ergodic, since the intensities are periodic and the state space is finite), if and only if see, e. g., [53]. For definiteness, assume that 1 0 µ(t) dt > 0. Apply the approach described in theorems 3 and 4.
Let all d k = 1. Then and in (34) we have α i (t) = µ(t) for all i, hence, α (t) = µ(t). Therefore, theorem 3 yields the estimate To find the constants in the estimates, let µ * = 1 0 µ(τ ) dτ and consider Find the bound for the second summand in (52). Assuming u = {t}, we obtain Then Therefore, for the queueing system M t |M t |N |N the conditions of theorem 5 and corollary 5 hold with These statements imply the following perturbation bounds: for the vector od=f state probabilities and for limit expectations.
Moreover, for these bounds to be consistent, the additional in=formation is required concerning the form of the perturbed intensity matrix. The simplest bounds can be obtained, if it is assumed that the perturbed Markov chain is also a BDP with the same state space and birth and death intensities λ k−1 (t) and µ k (t) respectively. Then if the birth and death intensities themselves do not exceed ε for almost all t ≥ 0, then |f −f| ≤ ε and |B −B| ≤ 5ε, so that the bounds (56) and (57) have the form for the vectors of state probabilities and for the limit expectations.
On the other hand, theorem 7 can be applied as well. To construct the bounds for the corresponding parameters, use (24) and the fact that D 1 = N . Then theorem 7 is valid for the queueing system M t |M t |N |N with the following values of the parameters: According to this theorem we obtain the estimate Moreover, the Markov chains X(t) andX(t) have limit expectations and It is worth noting that for estimates (61) and (62) to hold, only the condition of the smallness of perturbations is required and no additional information concerning the structure of the intensity matrix is required.
Thus, in the example with the finite state space uns=der consideration, uniform bounds turn out to be more exact.
On Fig. 1-5 there are the plots of the expected number of customers in the system for some of most probable states with ω = 1; on Fig. 6-7 there are the plots of of the expected number of customers with ω = 0.5.
On the other hand, as it has already been noted, for the Markov chains of classes I-IV with countable state space no uniform bounds could be constructed.
Consider the construction of bounds on the example of a rather simple model, which, however, does not belong to the most well-studied class I (that is, which is not a BDP). Let a queueing system be given in which the customers can appear separately or in pairs with the corresponding intensities a 1 (t) = λ(t) and a 2 (t) = 0.5λ(t), but are served one by one on one of two servers with constant intensities µ k (t) = min(k, 2)µ, where λ(t) is a 1-periodic function integrable on the interval [0, 1]. Then the number of customers in this system belongs to class II and the corresponding matrix B * * (t) has the form where a 11 (t) = − (1.5λ(t) + µ), a kk (t) = − (1.5λ(t) + 2µ), if k ≥ 2. This matrix is essentially nonnegative, so that in the expression for the logarithmic norm the signs of the absolute value can be omitted. Let d 1 = 1, d k+1 = δd k and choose δ > 1. For this purpose consider the expressions from (34). We have Then for δ ≤ 2 we obtain and the condition will a fortiori hold, if µ > λ * with a corresponding choice of δ ∈ (1, 2]. The further reasoning is almost the same as in the preceding example: instead of (54) we obtain where now K * = sup Hence, the conditions of theorem 5 and corollary 5 for the number of customers in the system under consideration hold for To construct meaningful perturbation bounds, it is necessarily required to have additional information concerning the form of the perturbed intensity matrix. So, example 1 in Section 2 shows that if a possibility of the arrival of an arbitrary number of customers ('mass arrival' in the terminology of [69]) to an empty queue is assumed, then an arbitrarily small (in the uniform norm) perturbation of the intensity matrix can 'spoil' all the characteristics of the process. For example, satisfactory bounds can be constructed, if we know that the intensity matrix of the perturbed system has the same form, that is, the customers can appear either separately or in pairs and are served one by one. Then, if the perturbations of the intensities themselves do not exceed ε for almost all t ≥ 0, then |f −f| ≤ 5ε and |B −B| ≤ 5ε, so that instead of (56) and (57) we obtain for the vectors of state probabilities and lim sup for the limit expectations. The particular example: let λ(t) = 1 + sin 2πt, µ(t) = 3. Choose δ = 2. Then we have α(t) = µ − 2.5λ(t), α * = 0.5, W = 1.
Further we follow the method that was described in [58,64] in detail. Namely, we choose the dimensionality of the truncated process (300 in our case), the interval on which the desired accuracy is achieved ([0, 100]) in the example under consideration) and the limit interval itself (here it is [100, 101]).     Fig. 5: Example 1. The perturbation bounds for the 'limit' probability Pr(X(t) = 210), t ∈ [19,20], ω = 1.