Waiting time and queue length analysis of Markov-modulated fluid priority queues

This paper considers a multi-type fluid queue with priority service. The input fluid rates are modulated by a Markov chain, which is common for all fluid types. The service rate of the queue is constant. Various performance measures are derived, including the Laplace–Stieltjes transform and the moments of the stationary waiting time of the fluid drops and the queue length distributions. An Erlangization-based numerical method is also provided to approximate the waiting time and the queue length distributions up to arbitrary precision. All performance measures are formulated as reward accumulation problems during busy periods of simple Markovian fluid flow models, for which efficient matrix-analytic solutions are provided, enabling us to solve large models with several hundred states.

conceptual simplicity makes priority systems easy to implement both in hardware and in software.
Queueing systems with priority service have been studied for a very long time. The first results for priority queues with exponentially distributed inter-arrival and service times have been considered in [1] more than 50 years ago. Nowadays, priority systems can be analyzed efficiently with complex arrival processes (like marked Markovian arrival processes) and more general service time distributions (see [2] or [3]).
There are queueing applications where the demands are easier to model by continuous fluid flows rather than discrete jobs. In these systems, the jobs are considered to be infinitesimally small, like fluid drops, and the queue storing the jobs is a fluid storage. The analysis of continuous (fluid based) queues has a long history as well; however, for some reason, the results for continuous queues have always been published years after their discrete counterparts. For instance, the matrix-analytic methods enabling the numerically efficient analysis of many discrete-state (non-fluid) queues have been adapted to continuous queues only in 1999 (see [4]).
Not many results exist for fluid queues with multi-type jobs. Among the existing results, the fluid queues with priority service have been studied the most, due to their practical relevance in the performance evaluation of various telecommunication systems. The authors in [5] and [6] assume simple on-off modulated rate for the high priority demands and constant fluid rate for the low priority demands, and obtain the Laplace-Stieltjes transform (LST) of the joint queue length distribution for this two-class system. The procedure proposed in [7] assumes simple on-off sources with exponentially distributed on and off periods for both the high and the low priority traffic. The novelty of that approach is that the authors could avoid the characterization of the joint queue length distribution. The idea is to create a simplified process to obtain the marginal distribution of the buffer content. In the simplified low priority queue length process, however, there are non-Markovian periods, representing the periods when the high priority class is busy. In the solution proposed in this paper, a somewhat similar reduction is introduced, too, but the reduction results in a completely Markovian system. The assumptions on the input process in the aforementioned papers, however, can be too restrictive for practical applications.
In this paper, a more general system is considered, where the fluid input rates of the arrival process are modulated by an arbitrary Markov chain. The same model has been studied in [8], where the partial differential equations for the joint distribution of the queue lengths are derived. However, these differential equations, especially the boundary functions, are difficult to solve. [8] is able to provide the mean, the variance and the covariance of the queue lengths only in case of two-state Markov chains. In [9], the LST of the stationary joint distribution of the buffer contents is obtained in a closed form based on a spectral decomposition approach, from which the tail distributions and the queue length moments are derived. Using Little's law, the mean waiting times are also presented. [10] follows a different approach to obtain the LST of the joint distribution, based on the analysis of the idle and busy periods of the high priority queue.
The (matrix-analytic) approach presented in this paper is based on the workload process analysis, just like in [3] for the discrete case. While the main steps of the analysis are the same, adapting the method of [3] to continuous queues is not straightforward at all. The specialties of continuous systems require different solutions in many steps of the procedure. Additionally, several ingredients of the procedure, which have been available for discrete systems for a long time (like the relation between the queue length distributions at departure instants and at random points in time) are not yet available for the continuous case.
The drawback of the presented approach is that the joint behavior of the queues cannot be characterized this way. (Note that the characterization of the joint queue length was the starting point of most related past papers.) The advantage of the approach presented in this paper is that it can provide the LST and the moments for both the individual queue lengths and the waiting times of fluid drops in a nice and elegant form. (Note that the waiting time of fluid drops in a fluid priority queue has never been considered in the past, apart from the expected value.) Furthermore, numerical methods for the waiting time and the queue length distributions based on Erlangization are also provided, which are created with probabilistic interpretations, and have better numerical behavior than a generic inverse Laplace transform algorithm.
It is particularly important for us to develop algorithms that behave well numerically, even for large models. All performance measures are formulated as reward accumulation problems during busy periods of simple Markovian fluid flow models, for which matrix-analytic solutions are provided. The computation bottlenecks are the solutions of Riccati-and Sylvester-type equations, for which efficient implementations exist. As a result, the presented procedure can solve large models up to many hundreds of phases, while the past solutions mentioned above are less tractable, they cannot go above a couple of tens of states (except for [9], where an efficient solution is provided if the background Markov chain is constructed by a Kronecker composition).
The rest of the paper is organized as follows: Section 2 describes the system, introduces the notation and lists the main steps of the computation of the performance measures. Markov-modulated fluid flow models are considered in Sect. 3. The first part of this section gives an overview of known results. In the second part of this section, the reward accumulated during the busy period is characterized. The waiting time of the fluid drops is derived in Sect. 4. The queue length distribution of a given fluid class is studied in Sect. 5. Finally, some numerical examples conclude the paper in Sect. 6.

The description of the system
In the system considered in the paper, K different (fluid) job types are distinguished, where class 1 has the lowest, and class K the highest priority. The priority of the fluid can be interpreted as "weight," and can correspond to the importance or the urgency of the jobs. The rates at which the fluid belonging to various job types enters the system are modulated by a continuous-time Markov chain (CTMC, also referred to as the background process) with generator denoted by Q. (The irreducibility of the background process is assumed throughout the paper.) The stationary probability vector of the CTMC is π ; hence, π Q = 0, π1 = 1 hold. (1 denotes the column vector of ones of appropriate size.) The rate at which class k fluid flows into the queue in state i is r . For simplicity, we denote the input rates of classes having priority equal to or higher than k by r i , and the corresponding matrix by . The rate at which fluid can leave the queue is denoted by d.
Using the notation introduced above, the mean arrival rate of class k fluid can be obtained by λ (k) = π R (k) 1, and the total input rate is λ = K k=1 λ (k) . The system is assumed to be stable; hence, λ < d.
In those states of the system where the fluid arrival rate is greater than the service rate available, the fluid is directed into a buffer. The buffer space is assumed to be infinite. Due to the priority service, lower priority fluid cannot be drained when any of the higher priority queues contain fluid. When the higher priority queues are empty and the input fluid rates of the higher priority classes are less than d, the lower priority queues are drained with the remaining service capacity. (Hence, the system is work conserving.)

Concept of the solution
In this paper, all performance measures will be derived from the analysis of the workload process. The workload, also called the unfinished work, or the virtual waiting time, is the amount of work present in the system, not in fluid level, but in time [11]. The fluid levels of different priorities are correlated; their joint analysis is rather complex. The workload process, however, makes it possible to obtain several performance measures from a Markov chain having only a single continuous variable, which is much more tractable numerically.
The solution is based on the "tagged customer" approach. When a class k fluid drop arrives into the queue, it finds a certain amount of class k+ workload in the system. This amount of workload has to be served before the tagged fluid drop can leave the system. The workload of lower priority classes can be neglected. While the tagged fluid drop waits for its service in the queue, further class (k + 1)+ fluid can arrive, which has to be served before the tagged one. After a certain amount of time, the tagged fluid drop reaches the head of the queue and can leave the system. Since each step of the solution relies on Markovian fluid flow models, the following section gives an overview of their stationary distribution and the behavior of the busy period.

Markovian fluid flow models
In this section, we recap the definition of the "standard" Markovian fluid flow models (without priorities) and those properties that will be necessary for the performance analysis of the fluid priority queue. The main purpose of this section is to introduce the notation and to review the related matrix-analytic solution methodology. The second part of this section presents some new contributions as well, related to the busy period analysis.
Markovian fluid flow models (or, in short, fluid models) consist of a continuoustime background Markov chain and a fluid storage (buffer). The rate at which fluid flows into or out from the buffer depends on the state of the background Markov chain. The buffer has a lower boundary (i.e., the minimal fluid level is 0), and its capacity is considered to be infinite.
Formally, fluid models are two-dimensional Markov processes {A(t), J (t), t > 0}, where A(t) is the fluid level and J (t) is the state of a CTMC (the background process) at time t.
The Markov chain is assumed to be irreducible with state space S; the number of states is N = |S|. Its generator matrix is denoted by F = [ f i j , i, j ∈ S]; the stationary distribution vector is ϑ, which satisfies ϑF = 0, ϑ1 = 1. A fluid rate is associated with each state of the background process, c i , i ∈ S, that determines the rate at which the level of the fluid buffer changes.
Since the analysis of the fluid priority queue involves the solution of various fluid flow models, we use different notation to make it easier to distinguish between the two systems. The generator matrix and the fluid input rates of the fluid priority queue are characterized by the matrices Q and R, while the generator and the fluid rates of fluid flow models are denoted by F and C throughout the paper. The evolution of the fluid level A(t) can be described by The size N row vector π(x) denotes the stationary density of the fluid level, whose ith entry is defined by At level 0, probability mass can accumulate as well. The ith entry of the row vector p represents the probability mass at level 0, p i = lim t→∞ P(A(t) = 0, J (t) = i). The vector π(x) is the solution of the matrix differential equation [12] d dt where C is a diagonal matrix composed of the fluid rates c i . The boundary equations providing π(0) and p are π(0)C = p F,

The stationary distribution of the fluid level
Several numerical methods are available to solve (2) and (3), including the eigenvalue decomposition-based method [12], the Schur decomposition-based method [13] and the matrix-analytic method [4]. The subsequent sections of this paper rely heavily on the matrix-analytic solution of fluid models; thus, the rest of the section provides a short summary of this approach. For the matrix-analytic method, the state space of the Markov chain has to be partitioned into three sets, S = S + ∪ S − ∪ S 0 , according to the sign of the fluid rates, i.e., . From now on, it is assumed that the matrices F and C are partitioned, and thus, where all diagonal entries of C + = diag c i i∈S + are positive and all diagonal entries of C − = diag c i i∈S − are negative. Furthermore, let us introduce the matrix F • as the generator of J (t) restricted to nonzero states, and hence, The stationary solution of the differential equation (2) and (3) has a matrixexponential form, as proven in [4]. If the drift of the queue is negative, ϑC1 < 0, the density function of the stationary fluid level can be expressed as: where the matrix is the minimal nonnegative solution of the non-symmetric algebraic Riccati equation (NARE) and the matrix K is obtained by The probability mass at level zero is p = 0 p − p 0 . The vectors κ, p − and p 0 are the solutions of the linear equations with the normalization condition

Reward accumulation during the busy period
The duration of the busy period of fluid flow models, τ = inf(t > 0 : A(t) = 0), has been studied extensively in the literature: in [14] the Laplace-Stieltjes transform (LST) of its distribution is provided, and in [15] an Erlangization-based algorithm is described to approximate its distribution in the time domain. Both of these papers are based on matrix-analytic methods. (We note that the busy period of Markovian fluid flows can be analyzed with other methodologies as well, like with the spectral decomposition method, as in [16].) In this section, we consider a more general problem. We assign a reward rate to each state of the background process and characterize the accumulated reward over the busy period of the fluid model. It turns out that the main performance measures of the fluid priority queue can be related to the accumulated reward over the busy period.
The introduction of reward accumulation makes the system three-dimensional, represents the reward accumulated up to time t. The reward accumulation is linear; thus, B(t) = t 0 d J (u) du, where d i denotes the reward rate associated with state i. Let us also introduce the diagonal matrix of reward rates D = diag d i . The joint distribution of the accumulated reward and the state transitions over the busy period is given by the matrix (y), defined as for i ∈ S + , j ∈ S − . Observe that in the special case of D = I, the matrix (y) equals the distribution of the busy period duration starting with zero initial fluid level.
The following theorem, from [17], provides the LST of (y), defined by

Theorem 1 [17, Theorem 3]
The matrix * (v), describing the LST of the accumulated reward over a busy period starting from an empty buffer, satisfies the NARE where the matrices F * with H * 00 (v) = (νD 0 − F 00 ) −1 being the LST of the reward accumulated in the zero states.
To obtain the moments of the performance measures of the fluid priority queue, the derivatives of the matrix * (v) will be required at v → 0. Let us introduce the matrices F The next theorem provides the derivatives of * (v).
for m > 0, and for m = 0 we have that (0) = = . The matrix U m−1 depends only on the derivatives of * (v) up to order m − 1: Proof Taking the mth derivative of (12) and letting v → 0 leads to which, after some manipulations, provides the theorem.
Appendix A describes how to obtain the matrix derivatives F The rest of this section focuses on the distribution of the reward accumulated over the busy period. Unfortunately, there are no exact solutions available to obtain (t). One can either rely on a generic numerical inverse Laplace transform (ILT) method to calculate it from * (v), or, one can apply the Erlangization method to approximate it. Most numerical ILT algorithms have several well-known drawbacks: they require complex and high precision arithmetic, have convergence problems and are prone to over-and undershooting, meaning that they can result in negative or higher than 1 values instead of valid probabilities. The Erlangization method can be viewed as a special kind of ILT method as well. Compared to generic ILT methods, Erlangization is slower, in the sense that it needs more evaluations of * (v) to achieve reasonable accuracy. On the other hand, it can be implemented using real, machine precision arithmetic and has no over-and undershooting problems; the result is always a valid probability. Additionally, Erlangization has an elegant stochastic interpretation, which can be exploited when tailored to specific problems. Below we present an Erlangization-based numerical method to approximate (t), following the approach of [15]. The basic algorithm in [15] is extended in two aspects: first, that paper studies the duration of the busy period, while we need the accumulated reward in this paper; second, zero rates are excluded in [15], while we need them here.
According to the Erlangization algorithm, the order-n approximation of the matrix (t) is where f E(n,n/t) (u) is an order-n Erlang density with rate parameter ν = n/t. As n tends to infinity, f E(n,n/t) (u) tends to a Dirac delta function δ(u − t); hence, n (t) tends to (t). The key idea behind Erlangization is that, by probabilistic interpretation, the integral in (15) is often easier to compute than (t) itself. In our case, the integral can be interpreted as the probability that the reward accumulated during the busy period is less than an Erlang(n, ν) variable.

Theorem 2 The order-n approximation of the joint distribution of the accumulated reward and the state transitions over the busy period is
where we haveˆ 0 (ν) = * (ν), and for > 0 the matricesˆ (ν) are obtained recursively by the solutions of the Sylvester equations where the matrices F * Proof The proof of the theorem is presented in Appendix B.

Waiting time analysis of the fluid priority queue
This section presents the LST, the moments and the distribution function of the waiting time of class k fluid drops, for k = 1, . . . , K .

The distribution of the workload at fluid drop arrivals
As the first step of the analysis, the distribution of the total workload (also called backlog, unfinished work, or virtual waiting time in the literature) of classes ≥ k is determined at class k fluid drop arrival instants. This is the amount of workload (and possibly more, due to further higher priority arrivals) that has to be accomplished before the tagged fluid drop can leave the system. In non-fluid queues, where the jobs are discrete, the difference between the discretestate queue length process and the continuous-state workload process is more obvious, but in the case of fluid queues both processes have continuous state space, and we will show that their dynamics are rather similar, too. Still, for the analysis of fluid priority queues, using the workload process is essential: it allows us to derive the performance measures from a process having only one continuous variable. Basing the analysis on the queue length process, we would need to solve a system with two continuous variables (the lengths of priority k and the priority (k + 1)+ queues-since they are correlated), which is much more involved.
Let us denote the workload of classes ≥ k at time t by V (k+) (t). It is important to note that the workload is measured in time, and not in fluid level. V (k+) (t) is the time needed by the system to accomplish all the work and become idle, given that no class ≥ k fluid enters the queue after time t. V (k+) (t) = 0 holds if and only if the class ≥ k fluid level is zero, but in general the evolution of the workload process differs from the evolution of the fluid level process.
The joint density function of the workload and the state of the background process at time t is denoted by the row vector v (k+) (t, x), whose elements are v The joint probability that the workload of the system is 0 and the state of the background process is i is stored by the row vector α (k+) with the boundary condition and α Proof where the first term corresponds to the case when there is no state transition in (t −Δ, t), and the second one when there is a state transition from j to i. Over a Δ long interval, the amount of class ≥ k workload brought into the system is Δr Rearranging the terms of (20) and dividing both sides by Δ leads to thus, either the workload was 0 in t − Δ, or it was positive, but close enough to zero such that it became zero in (t − Δ, t). Rearranging the terms, dividing by Δ and tending Δ → 0 gives (19). Lemma 1 has an important corollary: the workload process behaves like an ordinary Markovian fluid flow model with parameters thus, its stationary density, v (k+) (x) = lim t→∞ v (k+) (t, x), and probability mass at zero, α (k+) = lim t→∞ α (k+) (t), are given by the matrix-exponential form where the vectors κ (k+) , p (k+) and the matrices K (k+) , A (k+) are given by (6)- (10) with parameters (21). Hence, the density of the workload and the probability mass at zero embedded at class k fluid drop arrival instants arê

Obtaining the properties of the waiting time
As described in Sect. 2, the waiting time of a fluid drop is equal to the workload found in the system at its arrival, increased by the higher priority workload brought to the system while waiting. Figure 1 shows an example for the evolution of the remaining waiting time process of a class k fluid drop, W (k) (t, x), where the arriving fluid drop found x amount of class k+ workload in the system. In the time periods where the curve is drawn with a solid line, the background process is in a state where r ((k+1)+) = 0 and there is no class (k + 1)+ workload in the system either; hence, all the server capacity is spent on class k, and the remaining waiting time decreases with slope −1. In the intervals where the curve is dashed, class (k + 1)+ fluid is also present, and the service capacity is either fully taken by the higher priority classes, or, if their incoming rate is less than d and there is no high priority fluid accumulated in the buffer, it is shared between class k and the higher priority classes. While Fig. 1 depicts the remaining waiting time of a class k fluid drop, we can observe that it follows the same trajectory as the class (k + 1)+ workload process, V ((k+1)+) (t), given that V ((k+1)+) (0) = x.
Let T denote the random variable representing the waiting time of class k fluid drops. (For the sake of simplicity, we neglect the superscript (k) for T in this section.) t W (k) (t, x) Extra class (k + 1)+ work is served Extra class (k + 1)+ work is served x τ ((k+1)+) Fig. 1 The evolution of the remaining waiting time process t V ′((k+1)+) (t) Extra class (k + 1)+ work is served Extra class (k + 1)+ work is served Initial workload, x

Fig. 2 The class (k + 1)+ workload process after the transformation, V ((k+1)+) (t)
By conditioning on the amount of class k+ workload at a class k fluid drop arrival, we have that with G ((k+1)+) (y, x) being the distribution of the busy period of the class (k + 1)+ workload starting at level x, defined by where τ ((k+1)+) = inf(t > 0 : V ((k+1)+) (t) = 0) is the first time point when the workload process reaches level 0.
To get rid of the integral, we apply a transformation that transforms the analysis problem of the class k waiting time to the analysis problem of the accumulated reward of an appropriately constructed fluid flow model over a busy period. The main idea (introduced in [18]) is to start this fluid model at level 0 and add a phase whose role is to set the fluid level according to the initial workload at class k fluid arrival, having densityv (k) (x) (Fig. 2). The main benefit of this transformation is that it makes it possible to use the results of Sect. 3.2, where the busy period with 0 initial level is characterized.
The generator of the background process, the fluid rates and the reward rates of this special fluid flow model are The first state group of the background process is responsible for the accumulation of the workload a class k fluid drop finds in the system upon its arrival. Observe that the probability density of the time spent here equals e K (k+) t A (k+) R (k) /λ (k) , which, when started from the vector κ (k+) , equalsv (k) (t), according to (23). (Remember that κ (k+) is the initial vector of the stationary solution of the workload process; see (23).) While staying in this state group the workload increases with rate 1 (hence the I matrix in the first block of C), and during this period no reward is accumulated, since the time spent here is not part of the waiting time. (There is a 0 in the corresponding block in the matrix D.) When the appropriate workload level is reached, a transition occurs to the second group of states. In this part, the fluid level in this special model represents the workload of the system ahead of the tagged class k fluid drop. If class k has the highest priority in the system, then the workload ahead of the tagged fluid drop decreases linearly by rate 1. If there is higher priority work, class (k + 1)+ fluid flows in the system, then the workload brought by them increases the workload ahead of the tagged class k fluid drop. The rate at which the workload changes during this period is given by the diagonal elements of 1 d R ((k+1)+) − I. The reward rate is 1 in each state of the modulating CTMC; thus, the reward measures the time spent in this part of the state space. The time spent until the workload ahead of the tagged fluid drop reaches level 0 is identical to the waiting time.
Consequently, the waiting time is equal to the amount of reward accumulated during the busy period, starting from the first state group; thus, the initial phase distribution vector isκ Based on the above reasoning, we can state the following corollary.

Corollary 2 The LST of the distribution function of T , f * T (s) and its mth moment, E(T m ), are expressed by
and the order-n approximation of the distribution function, F (n) where the matrices * (s), (m) and n (t) are given by the reward analysis of the fluid model defined by the matrices (25), as detailed in Sect. 3.2.
In (29), the first termα (k) 1 reflects that the waiting time is zero if the workload is zero when the fluid drop arrives.

Queue length analysis of the fluid priority queue
To analyze the stationary queue length for each fluid class of the system, we first characterize the queue length at fluid drop departure instants. Then, we provide the relation between the queue length distribution at departures and the queue length distribution at random point in time.

Queue length distribution at fluid drop departures
The approach for characterizing the amount of class k fluid in the system at class k fluid drop departures is similar to the one presented in Sect. 4.2. A special Markovmodulated fluid model is constructed such that the reward accumulated during its busy period is equal to the queue length at departures. The matrices defining this special fluid model are Observe that these matrices are similar to (25), and the interpretation is similar as well. The first state group sets the initial workload seen by a class k fluid drop arrival. The second state group follows the workload ahead of the fluid drop until it leaves the system. With the given reward rate matrix D, however, the reward accumulated over the busy period measures the amount of class k fluid that arrives until the tagged fluid drop leaves; hence, it provides the class k queue length at the departure of the tagged class k fluid drop.

Corollary 3 Denoting the class k queue length at departures by X , the LST of the distribution function by f * X (s), its mth moment by E(X m ), and the order-n approximation of its distribution function by F
The definition and the computation of the matrices * (s), (m) and n (x) of the fluid model characterized by (30) are detailed in Sect. 3.2.
The first term in (33),α (k) , is the probability that the class k+ workload is zero at a class k fluid drop arrival, implying that the queue length is zero, too.

Queue length distribution at random point in time
In the case of non-fluid queues with discrete customers, the relation between the queue length distribution at departures and at random point in time is known [19]. For fluid queues, however, such a result is not available in the literature.
Let Y denote the class k queue length at random point in time. The density function, the distribution function and the probability mass at level 0 are denoted by f Y (x), F Y (x) and p Y , respectively. The following theorem establishes a relationship between the density and distribution functions of X and Y. Before stating the theorem, we introduce the vector forms of these quantities that include the state of the background process as well, characterizing lim t→∞ {Y(t), J (t)}. These row vector quantities are defined by

Theorem 3 The density and distribution functions of X and Y are related by
and for the masses at level 0 we have Proof The proof is inspired by [20] and is based on simple balance equations. Assume that the size of the fluid drops is δ. Expressing the rate at which the state {Y ∈ [0, x), J (t) = j} is left (for x > 0) we get where the first term is the rate of moving from {Y ∈ [0, x), J (t) = j} to {Y ≥ x, J (t) = j} due to a fluid drop arrival, and the second term is the rate of leaving {Y ∈ [0, x), J (t) = j} due to the state change of the background process. For x > 0, the rate at which the system enters state {Y ∈ [0, x), J = j} is where P(X ∈ [x − δ, x), J = j), in the first term, is the fraction of fluid drops that leave the system behind with fluid level < x.
] j , then switching to matrix notation, gives (34). The relation between the probability masses of X and Y at level 0 can be obtained using similar arguments. The rate of leaving state which, letting δ → 0, establishes (35).
It is worth noting that (35) is the continuous counterpart of [21,Theorem 5], while the derivative of (34) is the one of [21,Theorem 6]. The relation between X and Y in the Laplace transform domain, given by the following corollary, is similar to the result obtained for discrete-state (non-fluid) queues in [21,Equation (30)].

Corollary 4 The LSTs of the distribution functions of X and Y, denoted by f
Next, the class k queue length moments are derived based on the ones embedded at departure epochs. The moments are obtained by taking the derivatives of the LST, i.e., The following theorem (which is the continuous counterpart of [21, Section 4.1]) provides a recursive algorithm to compute the vectors y (m) .

Theorem 4 The vectors y (m) can be obtained recursively by
for m > 0, where y (m) 1 is given by The vectors x (m) are defined by x (m) = d m ds m f * X (s)| s→0 ; hence, For m = 0, we have y (0) = π , where π is the stationary probability vector of Q.
Proof For m = 0, the elements of the vector y (0) are the stationary probabilities of the background process, from which y (0) = π follows. Taking the derivative of (39) m times (m > 0) and setting s = 0 leads to Since the matrix Q is not invertible, y (m) cannot be expressed from (44) directly. Let us subtract y (m) 1π from both sides of (44). We get Observing that 1π − Q is invertible and that 1π(1π − Q) −1 = 1π provides (41). Multiplying both sides of (41) by R (k) 1 from the right yields Exploiting that multiplying (44) by 1 from the right gives y (m) R (k) 1 = λ (k) x (m) 1 and that λ (k) = π R (k) 1, (42) is proven.
With the theorems and corollaries presented above, it is possible to compute the LST and the moments of the queue length distribution. Finally, it remains to obtain the order-n approximation of the class k queue length distribution in the time domain, F (n) Y (x), based on Erlangization. Assume that the order-n approximation of the queue length distribution at fluid drop departures is available in the form of (33), F (n) (16)). The following theorem demonstrates how to obtain the approximation of F (n) Y (x) from the sameˆ (n/x) matrices.

Theorem 5 The order-n approximation of the distribution function F
for n > 1, and for n = 1, where ν = n/x.

Proof
The idea behind Erlangization is to express the distribution function as the limit of F (n) Multiplying both sides of (34) by an Erlang density f E(n,ν) (u) and taking the integral from 0 to ∞ gives which, when integrating both sides by parts, transforms to . The terms in the square brackets are zero, since F X (0) = p X and F Y (0) = p Y hold, and f E(n,ν) (u)| u=∞ = 0 holds, too. The remaining terms can be expressed as which equals (46) since F (n) (33) and (16)). Following similar steps for n = 1 gives

Numerical examples
The proposed algorithms have been implemented in the MATLAB environment and are available to download 1 . This implementation uses the ADDA algorithm [22] to solve the Riccati equations, and the lyap function of MATLAB to solve the Sylvester equations (the lyap function is based on the Hessenberg-Schur algorithm [23]). All numerical experiments have been executed on an average PC with a CPU clocked at 3.4 GHz and 4 GB of RAM.

Example 1
In the first example, there are three fluid types and the background process has four states. The generator of the CTMC and the fluid input rates in various states are given by the matrices and the service rate is d = 4. With these parameters, the utilization of the system is λ/d = 0.875. The algorithms returned the waiting time and queue length moments without a measurable delay. For the numeric results, see Table 1.
Obtaining the distributions is, however, more demanding computationally. The accuracy and the execution speed depend on how many terms are taken into account during the Erlangization. Our experience is in line with the earlier reports (see [15]),  namely that Erlangization is able to provide relatively accurate results with a small number of terms (the parameter n). Figure 3 depicts the distribution functions of class 1 waiting times for different values of n. For n = 10, the approximation has some visible error, but the approximate distribution functions for n ≥ 25 are almost indistinguishable from each other. Profiling this example revealed that the bottleneck of the Erlangization algorithm is the computation of the twofold summation at line 5 in (17). This term, however, is zero when the workload process of class 2 and 3 jobs does not have zero rates. The workload process does have a zero rate in this example, since in state 3 we have r 3 to a value different from 4 eliminated the zero state, leading to an order of magnitude faster execution times 2 ( Table 2).

Example 2
The aim of the second example is to demonstrate how the speed of the algorithms scales with the size of the input. Let us define the model as the superposition of N identical on-off sources with exponentially distributed on and off periods, with transition rates denoted by ν and γ , respectively. Each source generates fluid with rate ρ/N while being in the on state. The matrix parameters of the superposed fluid model are given by where N determines the size of the model. The matrix Q N is the generator of the background process, and the matrix R N is the diagonal matrix of the overall fluid rates. If we assume that 10% of the fluid generated by the sources is class 1 fluid, 30% is class 2 fluid and 60% is class 3 fluid, the per class fluid input rates are R In this example, the parameters are ν = 1/N and γ = 1/(3N ) (the reason for making them dependent on N is to keep the values of the generator in the same range irrespective of the size), and the fluid rate parameter is ρ = 50. With these settings, the mean fluid arrival rate is λ = 12.5.
We have computed the first five moments of the queue length of all job classes with different model sizes (controlled by the parameter N ) and with different utilizations u (controlled by the service rate d = λ/u). The execution times are depicted in Fig. 4, which reports the average of multiple repetitions for every point, together with the bootstrap-t confidence intervals (that are very tight).
As visible in the figure, the presented algorithms are able to solve huge models with 1000 states in reasonable time, without any numerical problems. The execution times depend on the utilization as well, since this affects the number of positive and negative states. At low utilization, obtaining five moments of the queue length for all job classes requires 10 s for N = 1000, while at high utilization about 40 s are required for the largest model. Taking the derivatives of F(v) is trivial if all entries of D are nonzero. If this is not the case, the matrices involved must be re-partitioned according to the sign of the reward rates, and hence, Note that this partitioning is not the same as the one assumed in Sect. 3. In Sect. 3, the partitioning is done according to the fluid rates; here, it is done according to the reward rates. According to [24], block matrices can be inverted using Introducing X(v) = (F ++ − vD + + F +0 (−F 00 ) −1 F 0+ ) −1 , this means that Since the only term in (50) which depends on v is X(v), the derivatives of (vD − F) −1 can now be easily obtained since Examining (15), it can be observed that the entries of the matrix n (t) are the probabilities that the accumulated reward is less than an order-n Erlang-distributed random variable with parameter ν. Hence, the matrix n (t) is obtained by creating a process where the reward accumulation process and the Erlang distribution are running in parallel and are competing. To this end, the state space of the background process (given by the generator matrix F) has to be extended such that it does not only keep track of the state of J (t), but also counts the number of events generated by a Poisson process having parameter ν. At the end of the busy period, the probabilities of those states, where the counter is less that n, contribute to n (t). To take the state-dependent reward rates into account, the transition rate incrementing the counter must be ν · d i if J (t) = i. Thus, the time is accelerated for the counting process if the reward rate is high in the current state, and it is slowed down if the reward rate is low. The generator matrix and the fluid rates of the modified process, which also includes the counter, have the following block structure: where being in block i at time t means that the value of the counter is i at time t. By partitioning the state space according to the sign of the fluid rates, the matrices F(ν) andC becomê where each of the 3 × 3 partitions of the matricesF(ν) andĈ have size infinity.