QUEUE LENGTH ANALYSIS OF A MARKOV-MODULATED VACATION QUEUE WITH DEPENDENT ARRIVAL AND SERVICE PROCESSES AND EXHAUSTIVE SERVICE POLICY

. The paper introduces a class of vacation queues where the arrival and service processes are modulated by the same Markov process, hence they can be dependent. The main result of the paper is the probability generating function for the number of jobs in the system. The analysis follows a matrix- analytic approach. A step of the analysis requires the evaluation of the busy period of a quasi birth death process with arbitrary initial level. This element can be useful in the analysis of other queueing models as well. We also discuss several special cases of the general model. We show that these special settings lead to simpliﬁcation of the solution.


1.
Introduction. The importance of vacation queues comes from their diverse application fields: modeling various computer systems, telecommunication protocols, manufacturing, logistics, etc. A large number of vacation queue variants have been analyzed in the literature, with different combinations of arrival processes, service time distributions, vacation time distributions and service policies, see [8] and [17] for recent surveys. The M/G/1 vacation model has already been extensively studied in the 1980's. For details on analysis works on classical vacation models the reader is referred to the survey of Doshi [5] or the book of Takagi [16].
Later on another vacation queue variants have been analyzed in the literature, like the GI/M/1 vacation queue in the work of Baba [3] or the discrete-time GI/Geom/1 vacation queue in the paper of J. H. Li et al. [10].
For more suitable modeling of manufacturing systems, real-life scenarios and specific modern telecommunication networks the basic vacation models have been generalized in several directions, like multi-server vacation models, working vacation models and application specific extended vacation models. For results on such models the reader is referred to [18], [14] and to surveys [8] and [17].
Due to the versatility of the Markovian Arrival Process (MAP) [11], vacation queues with MAP input and general service times have also been investigated in several past papers [4,12,13]. As a special variant of the MAP/G/1 vacation model, the MAP/PH/1 vacation model have been analyzed by A. S. Alpha for exhaustive time-limited and gated time-limited service disciplines in [1,2] and by C. Goswami and N. Selvaraju in [6]. All of these models assume independent service times. Up to the best knowledge of the authors the MAP/MAP/1 vacation queue, the vacation model where both the arrival and service processes are MAPs (hence, raising the opportunity of considering correlated inter-arrival and service times), has not been considered yet in the literature.
In this paper we consider a more general class of exhaustive vacation queues with dependent Markov modulated arrival and service processes (which we shortly refer to as Markov modulated vacation queue). The analysis of this model class requires the introduction of a new analysis element, the evaluation of the busy period of quasi birth death (QBD) processes with arbitrary initial level. As the main result of the paper we provide a computationally efficient expression of the probability-generating function (PGF) of the number of jobs in the general system. Additionally, we investigate the simplifications applicable for special model variants (e.g., for MAP/MAP/1 vacation model). In the last part of the paper we provide numerical examples and investigate the effects of different vacation distributions on the mean number of jobs.
The rest of the paper is organized as follows. In section 2 we introduce the general exhaustive vacation model with dependent Markov modulated arrival and service processes and discuss several special model variants. The analysis part of the work is presented in section 3. After a step-by-step establishment of the relations for characterizing the transient behavior of the model we derive the expression of the PGF of the number of jobs in the system. In section 4 the analysis of some special cases of the general model are provided. Numerical examples close the paper in section 5.
2. Model description. We consider the Markov modulated vacation queue. This model falls in the class of single server FCFS queues with multiple vacations, exhaustive discipline [17] and dependent Markov modulated arrival and service processes. According to the rule of the exhaustive service discipline the server serves the jobs in the queue until it gets idle, then the server leaves for vacation for an independent and identically distributed random amount of time. If the queue is idle at the end of the vacation the server leaves for a new vacation, otherwise it starts serving the jobs in the queue. The random vacation time, its probability density function (pdf) and its Laplace transform (LT) are denoted byσ, σ(t) and σ * (s) = E(e −sσ ), respectively.
The arrivals and services are characterized by seven matrices: L v , F v , B s , L s , F s , Π vs and Π sv . Their interpretations are as follows.
• During the vacations the arrivals are given by a MAP, where the entries of L v are the rates of transitions without a job arrival, and the entries of F v are the rates of transitions that are accompanied by a job arrival. Matrix L v + F v is therefore the generator of the continuous time Markov chain (CTMC) with N v states which modulates the arrivals during the vacation. • When the server is back from the vacation and serves jobs, the queue behaves like a quasi birth-death (QBD) [9] process, where the matrices B s , L s and F s contain the transition rates associated with a service completion, without service completion and job arrival, with a job arrival, respectively. In this case the generator of the modulating CTMC is B s + L s + F s , and it has N s states.
• The transition probabilities between the vacation and service periods are given by N v × N s stochastic matrix Π vs , whose entries are the probabilities of the state transitions occurring at the end of the vacation period. The probabilities of (extra) transitions (additional to those of B s ) between the service and vacation periods are given by the stochastic matrix Π sv of size N s × N v .
We impose the following assumptions on the model: The phase process of the QBD characterized by matrix B s + L s + F s and the phase process of the MAP characterized by matrix L v + F v are irreducible.
A.2 The rates of transitions in the QBD and in the MAP during the vacation are finite as well as the mean vacation time is positive and finite.
A. 3 The arrival process during the vacation and the sequence of the vacation periods are mutually independent.
A. 4 The model is stable. The QBD describing the service period is stable if and only if the stationary drift of its level process is negative [9]. Hence the necessary and sufficient condition of the stability of this vacation model is where 1 is the column vector of ones and α s is the stationary distribution of the phase process of the QBD and can be obtained as the unique solution of the system of linear equations α s (B s + L s + F s ) = 0 and α s 1 = 1. During the vacation period the queue length is monotone increasing, but due to the fact that the mean vacation time is finite the mean number of customers arrived during the vacation period is finite as well.
If the QBD describing the service period is positive recurrent the customers which arrived during the vacation period are served during the consecutive service period in finite time and this way the overall vacation model remains positive recurrent if the service period is positive recurrent.
2.1. Special model variants. The introduced Markov modulated vacation queue is a general model which covers a number of special cases. Some of those special cases are already discussed in the literature, e.g., MAP/PH/1 vacation queue and some have not been considered before, e.g., the QBD and the MAP/MAP/1 vacation queues. Below we list some special cases and define the corresponding model parameters.
• Markov modulated vacation queue without special phase change: In this case no additional phase transitions occur at the start and at the end of the vacation periods, i.e., Π vs = Π sv = I. As a consequence the matrices F v , L v , F s , L s and B s have the same size, N s = N v . • QBD vacation queue: A special case of the Markov modulated vacation queue is the QBD queue with vacations, which is similar to an ordinary QBD queue, but the transitions representing a service event in the service period are not accompanied by a service event in the vacation periods. The matrices characterizing the system are F v = F s = F, L s = L, B s = B, L v = B + L and Π vs = Π sv = I. Observe that the service times and the arrival times can be (cross) dependent in this case. • QBD vacation queue with independent arrival and service processes: A special case of the QBD vacation queue is when two independent Markov chains  Figure 1. Subset relations of the considered special vacation queue models modulate the arrival and the service processes. In this case the individual processes, i.e. the arrival process and the service process can be correlated, but the arrival and the service processes are independent. Let the arrivals be generated by a MAP with matrices (D 0 , D 1 ) and the service times be generated by a MAP with matrices (S 0 , S 1 ). In this case the characterizing matrices are • MAP/MAP/1 vacation queue: This model is very similar to the previous one.
A MAP with matrices (D 0 , D 1 ) generates the arrivals and a MAP with matrices (S 0 , S 1 ) generates the service times. The only difference compared to the previous model is the behavior during the vacation time: in the MAP/MAP/1 vacation queue the phase of the service process is frozen during the vacation. Consequently, the matrices characterizing the queue are Similar to the previous case, both the inter-arrival and service times can be correlated, but the arrival and service processes do not dependent on each other. In contrast with the QBD vacation queue with independent arrival and service processes the MAP/MAP/1 vacation queue is not a special case of the QBD vacation queue, but it is a special case of the Markov modulated vacation queue without special phase change. • MAP/PH/1 vacation queue: When the arrivals are generated by a MAP with matrices (D 0 , D 1 ) and the service time is phase type distributed with with initial vector α and generator matrix A, the characterizing matrices are In this case the service times are independent and identically distributed.
The subset relation of the listed special model sets is depicted in Figure 1. Apart of these standard cases the flexibility of the general model is exemplified by the following special models.
Vacation period Vacation period Service period Vacation period Figure 2. Cycles in the evolution of the queue • M/PH/1 queue in a random environment: Assume that the environment is given by a two-state CTMC with transition rates ν 1 , ν 2 . In state 1 the arrival rate is λ 1 , and the service times are phase-type (PH) distributed with initial vector α and generator matrix c 1 A. In state 2 the arrival rate is λ 2 and the service times are given by (α, c 2 A). The matrices characterizing the vacation queue are where 1 is the column vector of ones of compatible size. • MAP/MAP/1 vacation queue with discouraged customers: In this vacation queue the arrival process is different in the vacation and the service periods, due to the fact that discouraged customers do not attend the system during the vacation period. If the ratio of discouraged customers is then the characterizing matrices are The following section provides the analysis of the general vacation model. Some of the special cases are considered in section 4.
3. The number of jobs in the system. To characterize the number of jobs in the system, let us introduce the two dimensional process where N (t) denotes the number of jobs (also referred to as levels) and J (t) denotes the state of the modulating CTMC (also referred to as phase) at time t. For the analysis of X (t) the evolution of the queue is divided to cycles, as shown in Figure 2.
Each cycle starts with a vacation period, which is followed by a service period, and the cycle ends when the last job leaves the system. Note that a cycle can also be degenerate, if no jobs arrive during the vacation period, there is no service period (see cycle i − 1 in Figure 2).
The stationary probability that there are jobs in the system is closely related to M , the mean time spent at level in a stationary cycle.
where row vector β of size N s is the stationary phase distribution at the end of a service period (hence, βΠ sv is the one at the beginning of a vacation period), matrix P (v) (t) characterizes the number of arrivals up to time t during the vacation period, defined as and [H m, ] i,j is the mean times spent in level and phase j in the service period starting from level m and phase i. The first term of (2), M (v) , corresponds to the vacation period, and the second term, M (s) , to the service period. From M the Over the next subsections closed form formulas are provided for all ingredients of the solution, thus • for β in Section 3.2, • for H m, in Section 3.3, • for the PGFs of M (v) and M (s) in terms of in Section 3.4, finally, the pieces are put together in Section 3.5.
3.1. The evolution of the number of jobs during the vacation period. The evolution of the number of jobs during the vacation period resembles to the counting process of a MAP given by matrices with initial condition P (v) (0) = δ 0, I, where δ denotes the Kronecker delta. Similar to [9,Sec. 3], multiplying the th equation by z , summing up and solving the differential equation gives the generating function 3.2. The phase distribution at the end of the service period. The phase at the end of the cycles, J (t), forms an embedded discrete time Markov chain (DTMC).
The probability matrices characterizing the number of arriving jobs and the phase transitions during the vacation period are If m jobs are in the queue when the system enters the service period, the phase transitions over the busy period are given by G m , where matrix G is the minimal non-negative solution to the matrix-quadratic equation 0 = B s + L s G + F s G 2 ( [9]). Thus, the transition probability matrix of the DTMC, denoted by Q, is expressed by The stationary distribution of Q, denoted by β, is determined by the linear system βQ = β, β1 = 1.
Proof. In the first step the matrix-exponential representation for β Applying the column stacking vec operator and the identity vec Transposing (3), Kronecker multiplying both sides by G m , summing from m = 1 to ∞ and adding (4) gives Inserting this solution back to (8), multiplying it by σ(x) and taking the integral provides the theorem.
3.3. The mean time spent in different levels during the service period.
As a new contributions of the paper we derive matrix H m, , thus the mean time spent in various phases of level starting from level m in a QBD characterized by matrices B s , L s and F s . It is known that QBDs have a matrix-geometric stationary distribution with matrix R, and that the mean time spent at different phases of level starting from level 0 before returning to level 0 is given by R [9]. However, in our vacation queue the starting level after a vacation is not 0, but the number of arrivals accumulated during the vacation, which is denoted by m. According to our best knowledge, this measure has not been investigated yet.
For m > 0, we define matrix P where u marks the beginning and Θ marks the end of the service period, thus Θ = min{t : and for = 1 we have with initial values P where matrix Ψ is given as and matrices R and S are determined as the minimal non-negative solutions to the quadratic equations Proof. Integrating the differential equation (9) and (10) The solution of (14)-(17) is given by a matrix-geometric combination where matrices R and S are obtained such that the regular equations (15) and (

VACATION QUEUE WITH DEPENDENT ARRIVAL AND SERVICE PROCESSES 1373
The solution of Ψ is given by (12) and the solution of Φ is where we exploited various identities of the fundamental matrices of QBDs. Finally using the expressions of H m, from (19) and (18) as well as (12) and (20) establish the theorem. Theorem 3.2 reflects a quantitative behavior which is present also in the piecewise homogeneous QBD precess with a single point of inhomogeneity. Above the point of inhomogeneity (m < in our case) a matrix geometric behavior characterizes H m, with initial matrix Ψ(I − S m R m ) and matrix geometric term R −m . While below the point of inhomogeneity ( < m in our case) a combination of two matrix geometric series characterizes H m, with matrix geometric terms R and S m− .

Corollary 1. The generating function of H m, with respect to is expressed by
Proof. The generator function is derived by routine calculations as The following lemma provides some relations among the matrices playing an important role in the system, which will be referred several times in the sequel.
Proof. We note that for a stable vacation queue matrix I − SR and I − RS are non-singular, due to the fact that the spectral radius of S and R are one and less than one, respectively. Multiplying both sides of (12) by S(−SF s − L s − RB s ) from the left and exploiting that −S 2 F s − SL s = B s due to (13) gives which is equivalent to (22). (23) can be proven similarly, by multiplying (12) by R(−SF s − L s − RB s ) from the left.
Proving the third equality is a bit more involved. Making use of relation U = L s + RB s (valid for all QBDs, see [9]) and H 1,1 = (−U) −1 (by stochastic interpretation, i.e., both represent the mean time spent at level 1 before returning to level 0, starting from level 1) leads to since H 1,1 = Ψ(I − SR) according to (11).
Additional to matrices R and S, we introduce matrix G as the minimal nonnegative solutions to the quadratic matrix equation [9] The next two propositions are necessary to arrive to short queue length formulas at the end of the section.
Proof. The proof is based on simple algebraic manipulations. Based on Lemma 3.3, where σ * (M) with square matrix M is defined by Proof. By applying (5) in the first term of (2) the generating function M (v) * (z) can be expressed as The expression (27) comes by applying the shorthand notation σ * (M) in (29). M (s) * (z) is obtained by substituting H * m (z) provided by Corollary 1 into the definition (2). We get where we made use of Proposition 1 to write G m Ψ instead of ΨS m . After this swap, however, it can be recognized that the integral is equal to matrix Q (introduced by (6)), for which βQ = β holds. This observation simplifies M (s) * (z) to The final expression (28) is obtained by applying Proposition 2 in (30).
We note that the relative simplicity of (28) is due to the matrix swapping relation of Proposition 1.

3.5.
The generating function of the number of jobs in the system. The previous results lead to the following main theorem. Let q i be the stationary probability that there are i jobs in the system for i ≥ 0.
Theorem 3.5. The generating function of the number of jobs in the system, q * (z) = ∞ i=0 z i q i , is given by where β is determined by (7) and constant c satisfies lim z→1 q * (z) = 1.
Proof. The PGF of the mean time spent at different levels in the stationary cycle, The probability q i is proportional to M i , therefore q * (z) can be expressed as where c is an appropriate normalization constant. Inserting (27) and (28) into (32) provides the theorem.
Taking the derivatives of q * (z) at z → 1 gives the factorial moments of the number of jobs in the queue. 4. Analysis of special model variants. In this section we present simpler analysis results for special sub-classes of the general Markov modulated vacation queue.

4.1.
Stationary distribution at service completion. For the general Markov modulated vacation queue the phase distribution at the end of service, β, can be computed according to (7). Unfortunately, the numerical analysis based on (7) is a computationally demanding task, due to the Kronecker product in the evaluation of σ * (L T v ⊗ I + F T v ⊗ G). The following theorem presents a much simpler computation of β for a subset of Markov modulated vacation queues.
Theorem 4.1. In the special case of Markov modulated vacation queues, where G commutes with Π vs and G commutes with (L v + F v G), the stationary distribution at service completion, β, is the unique solution of Proof. The stationary distribution at service completion, β, satisfies the equation βQ = β. Starting from this and using (6) this equation for β can be rearranged as

VACATION QUEUE WITH DEPENDENT ARRIVAL AND SERVICE PROCESSES 1377
where in the last step we made use of the commutativity of G and Π vs . Next we derive an explicit form for m (x)G m for the case when the commutativity of G and (L v + F v G) holds. Multiplying differential equations (3) and (4) by G m from the right and summing them with regards to m leads to where the last step makes use of the commutativity of G and The theorem is proven by applying (36) in (34) and using the shorthand notation for the matrix LT.
The result of Theorem 4.1 makes a more efficient computation of β possible, as the matrices involved are much smaller due to the lack of Kronecker operations. For the special case when even the supplementary condition Π sv = Π vs = I holds, β can be obtained even simpler. Proof. Applying Π sv = Π vs = I to (33) leads to a simplified equation for β as βσ * (L v + F v G) = β. It can be seen from the Taylor expansion of the exponential term in the matrix LT that β fulfills the equation β(L v + F v G) = 0. Making use of the commutativity of G and (L v + F v G) we can write The matrix (L v +F v G) is a proper generator matrix, since it can contain negative elements only in the diagonal and (L v + F v G)1 = L v 1 + F v 1 = 0, since G is stochastic. Therefore the homogenous equation β(L v + F v G) = 0 determines β (up to a normalization constant) uniquely, which together with (38) implies that βG = β.
It is an interesting qualitative property in Theorem 4.2 that the stationary phase distribution at the end of service is independent of the vacation time distribution. This way Theorem 4.2 provides an easy to check condition (the commutativity of 2 matrices) for the stationary distribution to be independent of the vacation time and to be easy to compute.
The practical importance of Theorem 4.2 comes from the fact that the required commutativity holds in a number of practically important cases. E.g., in case of a MAP/MAP/1 vacation queue matrix (L v + F v G) = (D 0 ⊗ I) + (D 1 ⊗ I)G and matrix G do commute according to [7,Theorem 9], and the same relation applies for the MAP/PH/1 vacation queue as well. On the contrary, the commutativity relation of Theorem 4.2 does not hold in general for the QBD vacation queue and for the QBD vacation queue with independent arrival and service processes.

4.2.
Special formulas for the number of jobs. In the special case of Markov modulated vacation queue without special phase change (where Π sv = Π vs = I) Theorem 3.5 simplifies to the following corollary.
Corollary 2. When Π sv = Π vs = I the generating function of the number of jobs in the system is Further more, in case of the QBD vacation queue Theorem 3.5 further simplifies.
Both of the above expressions are products of transform domain expressions, which can be seen as factorization property. [4]. The second term of the product has a nice stochastic interpretation, because the vector generating function of the stationary queue length distribution of a QBD (without vacations) with regular matrices B, L, F and matrix L 0 = L + B at level zero is π * (z) = (1 − z)B(B + zL + z 2 F) −1 1.
Considering the subset relations of the set of special models depicted in Figure 1, we can see that the 2 branches have got different benefits. The MAP/MAP/1 vacation queue and the MAP/PH/1 vacation queue fulfills the commutativity condition of Theorem 4.2 and consequently the computation of vector β is computationally easy, while the QBD vacation queue and the QBD vacation queue with independent arrival and service processes do not fulfills the commutativity condition of Theorem 4.2 in general, but exhibit the factorization property.
4.3. The number of jobs in case of phase-type distributed vacations. If the vacations are phase-type (PH) distributed, the analysis of the number of jobs in the system is much simpler. Assuming that the density function of the vacations is σ(x) = αe Ax a, the queue length process X (t) is a CTMC with a QBD structured generator Q, which is

VACATION QUEUE WITH DEPENDENT ARRIVAL AND SERVICE PROCESSES 1379
The blocks of the QBD in the regular part (non-zero levels) are defined by The states belonging to level i are partitioned into two groups. The first state group corresponds to the vacation, the second to the service periods. The end of the vacation triggers a transition to the second state group and the phases are changed according to matrix Π vs .
Moving from the second group to the first one and initiating a vacation period can occur only when the last departure leaves the system (the QBD reaches level 0). The blocks for the irregular level 0 are The stationary distribution of X (t) is matrix-geometric where R is the minimal non-negative solution to the matrix-quadratic equation 0 = F + RL + R 2 B, while p 0 and p 1 are obtained by solving the linear system representing the equilibrium equations for levels 0 and 1 and the normalization condition: 5. Numerical example. The presented procedure has been implemented in Mathematica with high precision arithmetic 1 . σ * (X) for any matrix X is obtained by numerical integration, and the moments of the queue length are calculated by numerical derivation of q * (z). The first numerical experiment investigates the effect of the mean and the distribution of the vacation time on the mean number of jobs in the system. We consider the following matrices defining a QBD vacation queue.
In this example the arrival and service processes are dependent as much as possible: a service can occur only in state 1, and an arrival only in state 2. The computation has been performed for the following type of vacation distributions: Uniform distribution, Exponential distribution, and Weibull distribution with shape parameter k = 1/2. The mean number of jobs is depicted in Figure 3. As expected, the number of jobs in the system is the highest when the vacation times are Weibull distributed which has the heaviest tail.
The next experiment demonstrates the effect of Theorem 4.2, i.e., that β, the stationary phase distribution at the end of the service period does not depend on the vacation time distribution in the MAP/MAP/1 case. We have set up simple examples for all the special cases and computed β with various vacation time distributions. The results are shown in Table 1, according to which β does depend on  Table 1. Vector β as a function of the vacation distribution the vacation time distribution in all special cases considered in the paper, except the MAP/MAP/1 one. 6. Conclusions. In this paper we considered a general class of vacation queues with dependent arrival and service processes. Several previously studied vacation queues can be obtained as special cases of the general model as it is depicted on Figure 1. Based on the principles of the matrix-analytic approach we presented a computationally efficient description of the stationary queue length distribution in transform domain, which is based on an essential matrix swapping relation provided by Proposition 1. Finally, a numerical example demonstrates the stochastic behavior of the considered class of vacation queues, e.g. the effect of vacation time distribution. Unfortunately, we do not have nice stochastic interpretations for the majority of compact matrix expressions which were obtained in an essentially algebraic way. It is among our future research plans to find such interpretations which make these results easier to use.