Corrected phasetype approximations of heavytailed queueing models in a Markovian environment
Abstract
Significant correlations between arrivals of loadgenerating events make the numerical evaluation of the workload of a system a challenging problem. In this paper, we construct highly accurate approximations of the workload distribution of the MAP/G/1 queue that capture the tail behavior of the exact workload distribution and provide a bounded relative error. Motivated by statistical analysis, we consider the service times as a mixture of a phasetype and a heavytailed distribution. With the aid of perturbation analysis, we derive our approximations as a sum of the workload distribution of the MAP/PH/1 queue and a heavytailed component that depends on the perturbation parameter. We refer to our approximations as corrected phasetype approximations, and we exhibit their performance with a numerical study.
definition \@definecounterremark \@definecounterapproximation
Corrected phasetype approximations of heavytailed queueing models in a Markovian environment



Eurandom and Department of Mathematics & Computer Science, Eindhoven University of Technology,  
P.O. Box 513, 5600 MB Eindhoven, The Netherlands  
Centrum Wiskunde & Informatica (CWI), P.O. Box 94079, 1090 GB Amsterdam, The Netherlands  
Department of Mechanical Engineering, Eindhoven University of Technology,  
P.O. Box 513, 5600 MB Eindhoven, The Netherlands 
\@float
copyrightbox[b]
\end@float
Markovian Arrival Process (MAP), Workload distribution, Heavytailed service times, Tail asymptotics, Perturbation analysis
The evaluation of performance measures in stochastic models is a key problem that has been widely studied in the literature [?, ?, ?, ?]. In this paper, we focus on the evaluation of the workload distribution of a single server queue where customers arrive according to a Markovian Arrival Process (MAP) [?, ?] and their service times follow some general distribution. Under the presence of heavytailed service times, such evaluations become more challenging and sometimes even problematic [?, ?]. In such cases, it is necessary to construct approximations. In this study, we propose to modify existing approximations by adding a small refinement term, which can serve two purposes. On the one hand, the refinement term helps in constructing approximations not only with a small absolute error, but also with a small relative error. On the other hand, it gives information on the accuracy of the approximation without the modification: the smaller the refinement term, the better the premodified approximation.
An important generalization of the Poisson point process is the MAP. In a MAP, the arrivals are not homogenous in time, but they are determined by a Markov process with a finite state space. The class of MAPs is a very rich class of point processes, containing many wellknown arrival processes as special cases. A special case of a MAP is the Markovmodulated Poisson process (MMPP), which is a popular model for bursty arrivals [?]. The class of MAPs contains also the class of phasetype renewal processes, i.e. renewal processes with phasetype interarrivals [?].
It has been shown that the Laplace transform of the workload of a MAP/G/1 queue has a matrix expression analogous to the PollazceckKhinchine equation of an M/G/1 queue [?, ?]. However, these closedform expressions are only practical in case of phasetype service times [?, ?], where the workload distribution has a phasetype representation [?] in a form which is explicit up to the solution of a matrix functional equation.
Since the class of phasetype distributions is dense in the class of all distributions on [?], a common approach to approximate the workload is by approximating the service time distribution with a phasetype one; see e.g. [?, ?]. We refer to these methods as phasetype approximations. There are many algorithms for phasetype approximations, which provide highly accurate approximations for the workload distribution when the service times are lighttailed. However, in many cases, a heavytailed distribution is most appropriate to model the service times [?, ?]. In these cases, the exponential decay of phasetype approximations gives a big relative error at the tail and the evaluation of the workload becomes more complicated. Since heavytailed distributions have cumbersome expressions for their Laplace transform, this prevents the usage of techniques that require transform expressions, such as [?].
In this paper, we develop approximations of the workload distribution for heavytailed service times that maintain the computational tractability of phasetype approximations, capture the correct tail behavior and provide small absolute and relative errors. In order to achieve these desirable characteristics, our key idea is to use a mixture model for the service times. The idea of our approach stems from fitting procedures of the service time distribution to data. Heavytailed statistical analysis suggests that only a small fraction of the upperorder statistics of a sample is relevant for estimating tail probabilities [?]. The remaining data set may be used to fit the bulk of the distribution, where a natural choice is to fit a phasetype distribution to the remaining data set [?]. As a result, a mixture model for the service times is a natural assumption.
We now briefly explain how to derive our approximations when the service time distribution is a mixture of a phasetype distribution and a heavytailed one. We show that if the service time distribution is such a mixture, then the workload can also be written as a mixture, in the sense that it involves the workload of a model with purely phasetype service times and some additional terms related to the heavytailed distribution of our mixture model. Consequently, we first need to compute the workload in a MAP/PH/1 queue and afterwards use this as a base to calculate the rest of the terms involving the heavytailed distribution.
As a first step to derive our approximations, we write the service time distribution as perturbation of the phasetype distribution by a function that contains the heavytailed component. By ignoring the perturbation term and by taking the service time distribution equal to the phasetype distribution, we find the workload of a resulting simpler MAP/PH/1 queue, which is a phasetype approximation of the workload. By applying perturbation analysis to all parameters that depend on the service time distribution, we can write the workload as a series expansion, where the constant term is the workload of the MAP/PH/1 queue used as base and all other terms contain the heavytailed component.
Large deviations theory suggests that a single catastrophic event, i.e. a stationary heavytailed service time, is sufficient to give a nonzero tail probability for the workload [?]. As we will see in Section LABEL:s3.workload_distribution_of_the_perturbed_model, the second term of the series expansion of the workload can be expressed in terms of such a catastrophic event. Thus, we define our approximations as the sum of the first two terms of the series expansion of the workload, and we show that the addition of the second term leads to improved approximations when compared to their phasetype counterparts. In other words, the second term makes the phasetype approximation more robust so that the relative error at the tail does not explode. Therefore, we call this term correction term, and inspired by the terminology corrected heavy traffic approximations [?] we refer to our approximations as corrected phasetype approximations. In a previous study [?], we applied this approach to Poisson arrivals.
The connection between the stationary workload distribution of a MAP/G/1 queue and ruin probabilities for a risk process in a Markovian environment, where the claim sizes in the risk model correspond to the service times and the arrival process of claims is the timereversed MAP of the queueing model, is well known [?, ?]. Thus, the corrected phasetype approximations can also be used to estimate the ruin probabilities of the above mentioned risk model. Finally, our technique can be applied to more general queueing models, i.e. queuing models with dependencies between interarrival and service times [?, ?], and also to models that allow for customers to arrive in batches (the arrival process is called Batch Markovial Arrival Process) [?, ?, ?].
A closely related work is Adan and Kulkarni [?]. They consider a single server queue, where the interarrival times and the service times depend on a common discrete Markov Chain. In addition, they assume that a customer arrives in each phase transition, and they find a closed form expression for the waiting time distribution under general service time distributions. However, when there exist also phase transitions not related to arrivals of customers, their results remain valid for the evaluation of the workload. This can be seen by using the standard technique of including dummy customers in the model; namely customers with zero service time.
The rest of the paper is organized as follows. In Section Corrected phasetype approximations of heavytailed queueing models in a Markovian environment, we introduce the model under consideration without assuming any special form for the service time distribution, and in Section Corrected phasetype approximations of heavytailed queueing models in a Markovian environment we find the general expressions for the Laplace transforms of the workload prior to a transition from each state. In Section LABEL:s3.Construction_of_the_corrected_phasetype_approximations, we consider service time distributions that are a mixture of a phasetype distribution and a heavytailed one, and we explain the idea to construct our approximations. Later in Section LABEL:s3.replace_base_model, we specialize the results of Section Corrected phasetype approximations of heavytailed queueing models in a Markovian environment for phasetype service times. We use as base model the phasetype model of Section LABEL:s3.replace_base_model, and we apply perturbation analysis to find in Section LABEL:s3.parameters_of_the_perturbed_system the perturbed parameters and in Section LABEL:s3.workload_distribution_of_the_perturbed_model the desired Laplace transforms of the workload in the mixture model. Using the latter results, we construct in Section LABEL:s3.properties_corrected_replace_approximation the approximations and we discuss their properties. Finally, in Section LABEL:s3.numerics, we use a specific mixture service time distribution for which the exact workload distribution can be calculated and we exhibit the accuracy of our approximations through numerical experiments. Finally, in the Appendix, we give the proofs of all theorems, the necessary theory on perturbation analysis, and other related results. Due to the complexity of the formulas, we use a simple running example in order to explain the idea behind the calculations.
We consider a single server queue with FIFO discipline, where customers arrive according to a Markovian Arrival Process (MAP). The arrivals are regulated by a Markov process with a finite state space , say with states. We assume that the service time distribution of a customer is independent of the state of upon his arrival. For this model, we are interested in finding accurate approximations for the workload distribution.
The intensity matrix governing is denoted by the decomposition , where the matrix is related to arrivals of dummy customers, while transitions in are related to arrivals of real customers. Note that the diagonal elements of the matrix may not be identically equal to zero. This means that if , then a real customer arrives with rate and we have a transition from state to itself. However, phase transitions not associated with arrivals (dummy customers) from any state to itself are not allowed. Since the matrix is an intensity matrix, its rows sum up to zero. Therefore, the diagonal elements of the matrix are negative and they are defined as .
In this paper, we are interested in modeling heavytailed service times. As stated earlier, motivated by statistical analysis, we assume that the service time distribution of a real customer is a mixture of a phasetype distribution, , and a heavytailed one, . Namely, the service time distribution of a real customer has the form
(1) where is typically small.
Our goal is to find the workload distribution for this mixture model. Towards this direction, we present in the next section existing results [?] for the evaluation of the workload distribution under the assumption of generally distributed service times. Ultimately, we wish to specialize these results to service times of the aforementioned form (Corrected phasetype approximations of heavytailed queueing models in a Markovian environment).
Since the results of this section are valid for any service time distribution, we suppress the index and we use the notation for the service time distribution of a real customer. We consider now the embedded Markov chain on the arrival epochs of customers (real and dummy) and we denote by the transition probability matrix of the regulating Markov chain , which we assume to be irreducible. If is the exponential exit rate from state , i.e.
(2) the transition probabilities can be calculated by
(3) where is the Kronecker delta ( when and when ). In addition, an arriving customer at a transition from state to state is tagged . If , then we define the probability
(4) which is the probability of an arriving customer to be dummy conditioned on the event that there is a phase transition from state to . Similarly, conditioned on the event that there is a phase transition from to , the arriving customer is real with probability
(5) If , then we define . Consequently, the conditional service time distribution of an arriving customer at a transition from to is , and its LaplaceStieltjes transform (LST) is , , where is the LST of the service time distribution of a real customer. In matrix form, the above quantities can be written as
(6) (7) (8) (9) Let now denote the Hadamard product between two matrices of same dimensions; i.e. if and are matrices, then the element of the matrix is equal to . We also define the matrix
(10) which we will need later. Finally, let be the stationary distribution of , and be the mean of the service time distribution . Then the system is stable if the mean service time of a customer is less than the mean interarrival times between two consecutive customers in steady state. Namely,
(11) where and is the column vector with appropriate dimensions and all elements equal to 1. Note that the element of the matrix is the unconditional probability that a real customer arrives at a transition from to .
From this point on, we use a simple running example so that we display the involved parameters and the derived formulas. The running example evolves progressively, which means that its parameters are introduced only once and the reader should consult a previous block of the example to recall the notation.
For our running example, we consider a MAP with Erlang2 distributed interarrival times, where the exponential phases have both rate (). Therefore, the matrices and are given as follows:
In this case, we have that , , , and all other elements of the matrices and are equal to zero. Observe that we only have transitions from state 1 to state 2 and from state 2 to state 1. Therefore, in state 1 we always have arrivals of dummy customers while in state 2 we only have arrivals of real customers. Thus, only the diagonal elements of the matrix are not equal to zero, so that and . Finally, the stability condition takes its known form .
Let now denote the steadystate workload of the system just prior to an arrival of a customer. If the arriving customer is real, then the workload just prior to its arrival equals the waiting time of the customer in the queue, which we denote by . In terms of Laplace transforms, the steadystate workload of the system just prior to an arrival of a customer in state is found as
where is the steadystate limit of . Gathering all the above Laplace transforms , , we construct the transform vector
(12) We first provide some general theorems for the transform vector , which we later on refine in order to provide more detailed information regarding the form of the elements , . In the following, stands for the identity matrix, with appropriate dimensions.
Theorem 2.1
Provided that the stability condition (Corrected phasetype approximations of heavytailed queueing models in a Markovian environment) is satisfied, the transform vector satisfies
(13) (14) where is a vector with unknown parameters that needs to be determined.
Note that the above theorem is similar to Theorem 3.1 in [?] and so does its proof. Therefore, we omit here the proof and we refer the reader to Theorem 3.1 of [?] for more details.

Let be a column vector of dimension , such that . Then, it can easily be verified that is the Laplace transform of the waiting time of a real customer. If, however, , then is the Laplace transform of the workload just prior to an arrival of a customer. Thus, for the study of our system it is sufficient to determine the transform vector . \@endtheorem
If denotes the determinant of the square matrix , then for the determination of the unknown vector , we have the following theorem.
Theorem 2.2
The next two statements hold:

The equation has exactly solutions , with and for .

Suppose that the stability condition (Corrected phasetype approximations of heavytailed queueing models in a Markovian environment) is satisfied and that the above mentioned solutions are distinct. Let be a nonzero column vector satisfying
Then is given by the unique solution to the following linear equations:
(15) (16)
Again, Theorem 2.2 is similar to Theorems 3.2 & 3.3 in [?], and therefore, its proof is omitted.
Theorem 2.2 on one hand provides us with an algorithm to calculate the vector and on the other hand it guarantees that all elements of the transform vector are welldefined on the positive halfplane. To understand the latter remark observe the following. For simplicity, we set
(17) Let be the adjoint matrix of , so . Postmultiplying Eq. (13) with , we have that , and consequently
(18) The first statement of Theorem 2.2 says that the determinant has the factors , , in its expression. This means that the transform vector has potential singularities on the positive half plane, as the determinant appears at the denominator. However, the second statement of Theorem 2.2 explains that the vector is such that these problematic factors are canceled out.
Observe that Theorem 2.2 does not give us any information about the form of the elements of the transform vector , which is the stepping stone for the construction of our approximations. For this reason, we proceed by finding an analytic expression for the aforementioned elements. It is apparent from Eq. (18) that for the evaluation of we only need and the adjoint matrix . For the determination of these quantities, we introduce the following notation:

As before, we denote the set of all states of the Markov process as .

If , for some set , then is the complementary set of with respect to . Observe that all subset relations will be used locally and that the symbol “” does not imply strict subsets. The number of elements in a set is denoted as .

For a subset of we define and . We also define .

Suppose that and that is a square matrix of dimension . Then is the submatrix of if we keep the rows in and the columns in . Whenever the notation becomes very complicated, to avoid any confusion with the indices, we will denote the th column and row of matrix with and , respectively. We also define .

Suppose that is a subset of , for some set , and that it follows some properties, i.e. “Property 1", etc. If we want to sum with respect to , then we write under the symbol of summation first , followed by the properties. Namely, we write . In some cases, to avoid lengthy expressions we will write instead of the double sum , where is a subset of , for some set . We apply the same rule also for multiple sums.

Suppose that and are two square matrices of dimension , and that and are two disjoint subsets of . For all , we use the notation for the matrix that has columns the union of the columns of matrix and the columns of matrix , ordered according to the index set ; e.g. if , , and , then .
Using the above notation, we proceed with refining the desired quantities. More precisely, we first find , then the adjoint matrix , and finally the vector that appears in the numerator of the transform vector (see Eq. (18)). Combining these results, one can easily derive . We start by finding the determinant of matrix (see Eq. (17)).
Theorem 2.3
The determinant of matrix can be explicitly calculated as follows:

See Appendix missingLABEL:appendix_B.
Observe that the determinant is an at most degree polynomial with respect to the LST of the service time distribution of a real customer. Moreover, the coefficients of this polynomial are all polynomials with respect to . Therefore, in case is a rational function in , then is also a rational function in and its eigenvalues can be easily calculated. Furthermore, the subset of that appears in the second summand has at least one element, thus in the formula of it always holds that .
The matrix has elements , , , and . We calculate its determinant using Theorem 2.3. It holds that for all subsets of , except for . Since , it is evident that only for and , because the st column of the matrix and the nd column of the matrix are zero. Combining all these we obtain
In a similar manner, we find the explicit form of the adjoint matrix in the following theorem.
Theorem 2.4
The adjoint matrix has elements
where , , and .

See Appendix missingLABEL:appendix_B.
The adjoint matrix is equal to the transpose of the cofactor matrix of . Therefore, similarly to , each element of is an at most degree polynomial with respect to . This observation explains also the similarity between the formula of and the diagonal elements of .
Using the same arguments as for the evaluation of the determinant, we have for the adjoint matrix
Observe that the elements of the transform vector are defined as (see Eq. (18)), where is a column vector with element equal to 1 in position and all other elements zero. The outcome of is the inner product of the vector with the th column of matrix . Therefore, we have the following theorem.
Theorem 2.5
The numerator of the th element of the transform vector takes the form

See Appendix missingLABEL:appendix_B.
Combining now the results of the Theorems 2.3 and 2.5 by using Eq. (18), one can find the transform vector .
For each state we have
The transform vector is then
The following remark connects the system of equations that is required for the evaluation of , which was introduced in Theorem 2.2, to the adjoint matrix .

