Steady-state analysis of single exponential vacation in a $PH/MSP/1/\infty$ queue using roots

TL;DR

This paper analyzes a single-server queue with phase-type arrivals, Markovian service, and exponential vacations, deriving steady-state distributions and performance measures using root-based methods, with approximations and numerical insights.

Contribution

It introduces a root-based analytical approach for steady-state analysis of PH/MSP/1 queues with vacations, including new approximations and numerical evaluations.

Findings

Derived explicit steady-state distributions and performance metrics.

Developed heavy- and light-traffic approximations.

Provided numerical results illustrating parameter effects.

Abstract

We consider an infinite-buffer single-server queue where inter-arrival times are phase-type ($PH$), the service is provided according to Markovian service process $(MSP)$, and the server may take single, exponentially distributed vacations when the queue is empty. The proposed analysis is based on roots of the associated characteristic equation of the vector-generating function (VGF) of system-length distribution at a pre-arrival epoch. Also, we obtain the steady-state system-length distribution at an arbitrary epoch along with some important performance measures such as the mean number of customers in the system and the mean system sojourn time of a customer. Later, we have established heavy- and light-traffic approximations as well as an approximation for the tail probabilities at pre-arrival epoch based on one root of the characteristic equation. At the end, we present numerical results in the form of tables to show the effect of model parameters on the performance measures.

Equations

$$F_A(x)$$

$$\text{mbox}and{f}_A(x)$$

$$\frac{1}{\lambda }=\mathbf{\alpha }{\int }_0^{\infty }x{e}^{\mathbf{T}x}dx(-\mathbf{T}){\mathbf{e}}_{\eta }$$

$$P_{i,j}(n,t)=Pr\{N(t)=n,J(t)=j|N(0)=0,J(0)=i\},1\le i,j\le m.$$

$${\mathbf{P}}^{\ast }(z,t)=e^{\mathbf{L}(z)t},|z|\le 1,t\ge 0,$$

$$f_V^{\ast }(s)=\frac{\gamma }{\gamma +s},f_V(x)=\gamma e^{-\gamma x},F_V(x)=1-e^{-\gamma x}.$$

$${\mathbf{S}}_n={\int }_0^{\infty }\mathbf{P}(n,t)d{F}_A(t),n\ge 0.$$

$$\mathbf{S}(z)$$

$$F_{\hat{A}}(x)=F_{\tilde{A}}(x)={\int }_0^x\lambda (1-F_A(y))dy,$$

$$F_{\hat{V}}(x)$$

$$f_{\hat{V}}(x)$$

$${\mathbf{\Omega }}_n$$

$${\mathbf{S}}_n$$

$${\tilde{P}}_{ij}(n,t+\Delta t)$$

$$\frac{d}{dt}{\tilde{P}}_{ij}(n,t)$$

$$\frac{d}{dt}\tilde{\mathbf{P}}(n,t)$$

$$\omega$$

$$\tau$$

$$\omega$$

$$\tau$$

$$F_{A^{+}}(x)$$

$$F_{{\hat{A}}^{+}}(x)=F_{{\tilde{A}}^{+}}(x)={\int }_0^x{\lambda }_1(1-F_{A^{+}}(y))dy,$$

$${\mathbf{\alpha }}_1$$

$${\mathbf{T}}_1$$

$${\mathbf{V}}_n$$

$${\mathbf{V}}_n^{\ast }$$

$${\mathbf{\Delta }}_n$$

$${\mathbf{\Delta }}_0$$

$${\mathbf{\Delta }}_n$$

$${\mathbf{V}}_n^{\ast }$$

Full text

Keywords: Markovian service process ($MSP$), renewal input, infinite-buffer, exponential single vacation, log-normal inter-arrival, phase-type approximation

1. 1 Introduction

2. 2 Description of the model

3. 3 The vector generating function of the number of customers served during an inter-arrival and other related probability matrices

4. 4 Analysis of $GI/MSP/1/\infty$ queue with single vacation

5. 4.1 Stationary system-length distribution at pre-arrival epoch

6. 4.2 Stationary system-length distribution at arbitrary epoch

7. 4.3 Queue-length distribution at post-departure epoch and their relation with pre-service epoch

8. 5 Performance measures

9. 5.1 **Waiting-time analysis **

10. 5.2 **Expected length of busy and idle periods **

11. 6 Approximation for system-length distributions based on one root

12. 7 Numerical results and discussion

13. 8 Conclusions and future scope

1 Introduction

In recent times, queueing models with non-renewal arrival and service processes have been used to model networks of complex computer and communication systems. Traditional queueing analysis using Poisson processes is not powerful enough to capture the correlated nature of arrival (service) processes. The performance analysis of correlated type of arrival processes may be done through some analytically tractable arrival process viz., Markovian arrival process ( $MAP),$ see Lucantoni et al. ${\hbox{\cite[cite]{[\@@bibref{}{nu4}{}{}]}}}.$ The $MAP$ has the property of both time varying arrival rates and correlation between inter-arrival times. To consider batch arrivals of variable capacity, Lucantoni [20] introduced batch Markovian arrival process ( $BMAP)$. The processes $MAP$ and $BMAP$ are convenient representations of a versatile Markovian point process, see Neuts [23] and Ramaswami [27]. Like the $MAP$, Markovian service process $(MSP)$ is a versatile service process which can capture the correlation among the successive service times. Several other service processes, e.g., Poisson process, Markov modulated Poisson process ($MMPP$) and phase-type ($PH$) renewal process can be considered as special cases of $MSP$. For details of $MSP$, the readers are referred to Bocharov [3] and Albores and Tajonar [1]. The analysis of finite-buffer $G/MSP/1/r~{}(r\leq\infty)$ queue has been performed by Bocharov et al. [4]. The same queueing system with multiple servers such as $GI/MSP/c/r$ has been analyzed by Albores and Tajonar [1]. Gupta and Banik [17] analyzed $GI/MSP/1$ queue with finite- as well as infinite-buffer capacity using a combination of embedded Markov chain and supplementary variable method.

During the last two decades, queueing systems with vacations have been studied extensively. For more details on this topic, the readers are referred to a recent book by Tian and Zhang [30] and references therein. An extensive amount of literature is available on infinite- and finite-buffer $M/G/1$- and $GI/M/1$-type queueing models with multiple vacations, see first few chapters of [30], Karaesmen and Gupta [18] and Tian et al. [29]. However, limited studies have been done on $GI/M/1$ queue with single vacation, see Chapter 4 of [30]. In the past few years, there is a growing trend to analyze queueing models with renewal or non-renewal arrival and service processes with server vacation, see, e.g., Lucantoni et al. [22] and Shin and Pearce [26]. The analysis of phase-type server vacation for the case of $GI/M/1$ queue has been carried out by Chen et al. [7]. Baba [2] analyzes $M/PH/1$ queue where the server is allowed to take working vacations as well as vacation interruptions. Samanta [28] discussed a discrete-time $GI/Geo/1$ queue with single geometric vacation time. Recenly, Chaudhry et al. [9, 12] discussed $GI/MSP/1/\infty$ queues with single and batch arrivals using the roots method, respectively.

In this paper, we carry out the analytic analysis of the $PH/MSP/1/\infty$ queue with exponential single vacation through the calculation of roots of the denominator of the underlying vector generating function of the steady-state probabilities at pre-arrival epoch. In this connection, the readers are referred to Chaudhry et al. [9, 10, 12], Tijms [31] and Chaudhry et al. [11] who have used the roots method. The roots can be easily found using one of the several commercially available packages such as Maple and Mathematica. The algorithm for finding such roots is available in some papers, e.g., see Chaudhry et al. [11]. The purpose of studying this queueing model using roots is that we obtain computationally simple and analytically closed form solution to the infinite-buffer $PH/MSP/1$ queue with the vacation time following exponential distribution. It may be remarked here that the matrix-geometric method (MGM) uses iterative procedure to get steady-state probabilities at the pre-arrival epochs. Further, it is well known that for the case of the MGM it is required to solve the non-linear matrix equation with the dimension of each matrix in this equation being the number of service-phases involved in a $PH/MSP/1$ queue. In the case of the roots method, we do not have to investigate the structure of the transition probability matrices (TPM) at the embedded pre-arrival epochs. It may be mentioned here that the basic idea of correlated service was first introduced by Chaudhry [13]. Further, it may be remarked here that the analysis of the infinite-buffer queues with renewal input and exponential service time under exponential server vacation(s) has been carried out by Tian and Zhang [30], see Chapter 4. The queueing model that we are going to consider has non-renewal service ($MSP$) and exponential single vacation time. In addition, we discuss several other quantitative measures such as system-length distribution at a post-departure epoch and expected busy and idle periods. Later, we have established heavy- and light-traffic approximations as well as an approximation for the tail probabilities at pre-arrival epoch based on one root of the characteristic equation. Finally, some numerical results have been presented which may help researchers/practitioners to tally their results with those of ours.

2 Description of the model

Let us consider a single-server infinite-buffer queueing system with the server’s single vacation. The inter-arrival time of customers, the service time of a customer and the vacation time of the server are represented by the generic random variables (r.v.’s) $A$, $S$ and $V$, respectively. Let $F_{X}(x)$ denote the distribution function (D. F.) of the random variable $X$ with $f_{X}(x)$ and $f_{X}^{*}(s)$ the corresponding probability density function (p.d.f.) and Laplace-Stieltjes transform (LST), respectively. The inter-arrival time $A$ is assumed to have a general distribution with p.d.f. $f_{A}(x)$, D. F. $F_{A}(x)$ and LST $f_{A}^{*}(s)$.

Arrivals. The inter-arrival times are assumed to be independent and identically-distributed (i.i.d.) random variables and they are independent of the service process as well as vacation time. The inter-arrival time distribution $PH$ is an important special case of general distribution as the distribution possesses nice vector and matrix form representation. Several probability distributions such as Earlang, hyper-exponential, generalized Earlang, Coxian etc. can be treated as special cases of $PH$-distribution. It may be noted here that $PH$-distribution is a special case of general distribution. If the inter-arrival times follow $PH$-type distribution with irreducible representation $({\bm{\alpha}},{\bm{T}})$, where ${\bm{\alpha}}$ & ${\bm{T}}$ are a vector and a matrix of dimension $1\times\eta$ and $\eta\times\eta$, respectively, the p.d.f. and D.F. of inter-arrival times are given by

[TABLE]

where ${{\bm{T}}}^{0}$ is a non-negative vector and satisfies ${\bm{T}}{\bm{e}}_{\eta}+{{\bm{T}}}^{0}={\bm{0}}$ and ${\bm{e}}_{\eta}$ is an $\eta\times 1$ vector with all its elements equal to 1. Throughout the paper we write a subscript as the dimension of the column vector ${\bm{e}}$ and sometimes we write ${\bm{e}}$ by dropping its subscript. The mean inter-arrival time during a normal busy period is given by

[TABLE]

Services. The customers are served singly according to the continuous-time Markovian service process ($MSP$) with matrix representation $({\bm{L}}_{0},{\bm{L}}_{1})$. The $MSP$ is a generalization of the Poisson process where the services are governed by an underlying $m$-state Markov chain. For more details on $MSP$, the readers are referred to recent papers by Chaudhry et al. [9, 12]. Let $N(t)$ denote the number of customers served in $t$ units of time and $J(t)$ the state of the underlying Markov chain at time $t$ with its state space $\{i:1\leq i\leq m\}$. Then $\{N(t),J(t)\}$ is a two-dimensional Markov process with state space $\{(\ell,i):\ell\geq 0,1\leq i\leq m\}$. Average service rate of customers $\mu^{\star}$ (the so called fundamental service rate) of the stationary $MSP$ is given by \mu^{\star}=\overline{\mbox{\boldmath\pi}}{\bm{L}}_{1}{\bm{e}}, where \overline{\mbox{\boldmath\pi}}=[\overline{\pi}_{1},\overline{\pi}_{2},\ldots,\overline{\pi}_{m}] with $\overline{\pi}_{j}$ denoting the steady-state probability of servicing a customer in phase $j~{}(1\leq j\leq m)$. The stationary probability row-vector \overline{\mbox{\boldmath\pi}} can be calculated from \overline{\mbox{\boldmath\pi}}{\bm{L}}={\bm{0}} with \overline{\mbox{\boldmath\pi}}{\bm{e}}=1, where ${\bm{L}}={\bm{L}}_{0}+{\bm{L}}_{1}$. The customers are served singly according to a $MSP$ with steady-state mean service time $1/\mu^{*}$.

Now, let us define $\{{\bm{P}}(n,t):n\geq 0,t\geq 0\}$ as the $m\times m$ matrix whose $(i,j)$th element is the conditional probability defined as

[TABLE]

Let ${\bm{P}}(n,t),n\geq 0,t\geq 0$ be the $m\times m$ matrices whose elements are $P_{i,j}(n,t)$. Then using Chaudhry et al. [9, 12], it may be derived that

[TABLE]

where ${\bm{L}}(z)={\bm{L}}_{0}+{\bm{L}}_{1}z$ and ${\bm{P}}^{\ast}(z,t)=\sum\limits_{n=0}^{\infty}{\bm{P}}(n,t)z^{n},\quad|z|\leq 1.$

Vacations. The server is allowed to take a single vacation whenever the system becomes empty. On return from a vacation if the server finds the system nonempty he will serve the customers present in the queue, otherwise the server waits for a customer to arrive and the system continues in this manner. For an exponential single vacation time represented by the r.v. $V$, the LST, p.d.f. and D.F. are given as follows:

[TABLE]

where $1/\gamma~{}(>0)$ is assumed as the mean vacation time. The Vacation times are independent of the arrival as well as of the service processes. The traffic intensity is given by $\rho=\lambda E(S)=\lambda/\mu^{*}$ which is also independent of the vacation process.

3 The vector generating function of the number of customers served during an inter-arrival and other related probability matrices

Let ${\bm{S}}_{n}~{}(n\geq 0)$ denote the matrix of order $m\times m$ whose $(i,j)$th element represents the conditional probability that during an inter-arrival period $n$ customers are served and the service process passes to phase $j$, provided at the initial instant of the previous arrival epoch there were at least $n$ customers in the system and the service process was in phase $i$. Then

[TABLE]

If ${\bm{S}}(z)$ is the matrix-generating function of ${\bm{S}}_{n}$, where $S_{i,j}(z)~{}(1\leq i,j\leq m)$ are the elements of ${\bm{S}}(z)$, then, using (6) and (4), we get

[TABLE]

The evaluation of the matrices ${\bm{S}}_{n}$ can be carried out along the lines proposed by Lucantoni [20]. For the sake of completeness, we have given the procedure of obtaining ${\bm{S}}_{n}$, see Lucantoni [20]. One may note that the computation of ${\bm{S}}(z)$ using Equation (7) may be cumbersome. However, the following scheme may be efficient and is given by

[TABLE]

where ${\bm{S}}_{n}$ may be obtained as proposed in Chaudhry et al. [10].

We further introduce a few more notations which are required for the rest of the analysis of the queueing model under consideration. Now from renewal theory of semi-Markov process, if we let $\widehat{A}$ and $\widetilde{A}$ denote the remaining and elapsed times of an inter-arrival time, respectively, then

[TABLE]

which we use while deriving the expression for ${\bm{\Omega}}_{n}$ in Equation (12). Similar to the case of inter-arrival time, if we let $\widehat{V}$ denote the remaining vacation time, then

[TABLE]

and

[TABLE]

As above, we introduce the matrices ${\bm{\Omega}}_{n}~{}(n\geq 0)$ of order $m\times m$ whose $(i,j)$th element represents the limiting probability that $n$ customers are served during an elapsed inter-arrival time of the arrival process with the service process being in phase $j$, given that there were at least $(n+1)$ customers in the system with the service process being in phase $i$ at the beginning of the inter-arrival period. Then, from Markov renewal theory as given in Chaudhry and Templeton [8, p. 74-77], we have

[TABLE]

The matrices ${\bm{\Omega}}_{n}$ can be expressed in terms of the matrices ${\bm{S}}_{n}$ and their relationship discussed in [12] is as follows:

[TABLE]

where ${\it 1_{\{n\geq 1\}}}$ is an indicator function and takes value 1 if the the condition $n\geq 1$ is satisfied, otherwise it takes value 0.

Further, let $\widetilde{P}_{ij}(n,t)$ be the conditional probability that at least $n$ customers are served in $(0,t]$ and the service process is in phase $j$ at the end of the $n$th service completion, given that there were $n$ customers in the system and the service process was in phase $i$ at time $t=0$. The probabilities $\widetilde{P}_{ij}(n,t),~{}n\geq 1,~{}t\geq 0,$ then satisfy the equations

[TABLE]

with the initial condition $\widetilde{P}_{ij}(n,0)=0,~{}n\geq 1.$ Rearranging the terms and taking the limit as $\Delta t\to 0$, it reduces to

[TABLE]

for $t\geq 0,~{}1\leq i,j\leq m$, with the initial conditions $\widetilde{P}_{ij}(n,0)=0$. This system may be written in matrix notation as

[TABLE]

with $\widetilde{{\bm{P}}}(n,0)={\bm{0}},~{}~{}n\geq 1$.

Let $\omega$ denote the probability that $\widehat{V}$ exceeds an inter-arrival time $A$, then

[TABLE]

Similarly, if we let $\tau$ denote the probability that $\widehat{V}$ exceeds $\widehat{A}$, then

[TABLE]

Remark 2.1: $\omega$ and $\tau$ may be derived in slightly different way. In the following we present a slightly different derivation for $\omega$ and $\tau$ may be done similarly.

[TABLE]

One of the frequently used inter-arrival time is phase-type renewal process which also serves as a special case of several other inter-arrival time distributions and is well-known in the literature. Therefore, we state the above formulae (17) and (16) for the case of phase-type inter-arrival time by the following theorem.

Theorem 3.1

If inter-arrival time follows a $PH$-distribution with irreducible representation $({\bm{\alpha}},{\bm{T}})$, where ${\bm{\alpha}}$ and ${\bm{T}}$ are of dimension $\eta$, then the expressions for $\omega$ and $\tau$ are as follows.

[TABLE]

Proof:* Using the definition of $\omega$ and $\tau$, after little algebraic manipulation the the results (35) and (19) may be obtained.*

In the following, we further define a few notations which are required to analyze the queueing model under consideration. If we let $A^{+}=A-\widehat{V}$ with $A-\widehat{V}>0$, i.e., inter-arrival time is greater than remaining vacation time, then

[TABLE]

where $A^{+}$ may be called excess inter-arrival time. Further, let us denote $\widehat{A}^{+}$ and $\widetilde{A}^{+}$ as the remaining and elapsed times of the excess inter-arrival time random variable $A^{+}$, respectively, then

[TABLE]

where $\lambda_{1}=\int_{0}^{\infty}x\ dF_{A^{+}}(x)$ is the mean of the random variable $A^{+}.$ For an important special case of phase-type inter-arrival time the distribution of excess inter-arrival time $A^{+}$ is also phase-type which is proved in the following theorem.

Theorem 3.2

If inter-arrival times follow a $PH$-distribution with irreducible representation $({\bm{\alpha}},{\bm{T}})$, where ${\bm{\alpha}}$ and ${\bm{T}}$ are of dimension $\eta$, then the distribution of $A^{+}$ is also phase-type with representation $({\bm{\alpha}}_{1},{\bm{T}}_{1})$, where ${\bm{\alpha}}_{1}$ and ${\bm{T}}_{1}$ are of dimension $\eta$ and are given by

[TABLE]

Proof:* Simple algebraic calculation will give the proof.*

Further, we introduce a few more matrices which are required for the analysis of this queueing model. With vacation ending, let ${\bm{V}}_{n}~{}(n\geq 0)$ denote the matrix of order $m\times m$ whose $(i,j)$th element represents the conditional probability that $n$ customers are served during an excess inter-arrival time $A^{+}$ and the service process passes to phase $j$, provided at the initial instant of the previous arrival epoch there were at least $(n+1)$ customers in the system and the server was on a vacation with the service phase $i$. Then

[TABLE]

Now, with vacation ending, let ${\bm{V}}_{n}^{*}$ denote the matrix of order $m\times m$ whose $(i,j)$th element represents the probability that at least $n$ customers are served during an excess inter-arrival period $A^{+}$ and the service process is in phase $j$ with the server going on vacation at the end of the $n$-th service completion, provided at the initial instant of previous arrival epoch there were exactly $n$ customers in the system and the server was on a vacation with the service phase $i$. Then similar to the results derived above, we obtain

[TABLE]

Further, with vacation ending, let ${\bm{\Delta}}_{n}~{}(n\geq 0)$ denote the matrix of order $m\times m$ whose $(i,j)$th element represents the limiting probability that $n$ customers are served during an elapsed excess inter-arrival time $\widetilde{A}^{+}$ with the service process being in phase $j$, given that there were at least $(n+1)$ customers in the system with the server being on a vacation with the service phase $i$ at the beginning of the inter-arrival period. Then we have the following expression for ${\bm{\Delta}}_{n}$

[TABLE]

The relationships among the matrices ${\bm{V}}_{n}$, ${\bm{V}}_{n}^{*}$ and ${\bm{\Delta}}_{n}$ can be derived as follows.

[TABLE]

and

[TABLE]

Further, using [12], it can be shown that

[TABLE]

Similarly, let ${\bm{\Delta}}_{n}^{*}~{}(n\geq 1)$ denote the $m\times m$ matrix whose $(i,j)$th element represents the limiting probability that $n$ or more customers have been served during an elapsed excess inter-arrival time $\widetilde{A}^{+}$ and the service process is in phase $j$ with the server going on vacation at the end of $n$th service completion, provided at the previous arrival epoch the server was on a vacation with service phase $i$ and the arrival lead the system to state $n$ or more customers. Then we can write

[TABLE]

Now using the procedure discussed in [12], one may derive the following relation:

[TABLE]

Finally, we define a few notations which are required to analyze the queueing model under consideration. Let $A^{++}=A-(\widehat{V}+V)$ given that $A-(\widehat{V}+V)>0$. It is needless to mention that since $V$ is exponentially distributed, the distribution of $\widehat{V}+V$ will be Erlang of order two, which is a phase-type distribution with two states. Let the phase type representation of $\widehat{V}+V$ be denoted as \mbox{\boldmath\beta}=\left[\matrix{1.0&0.0}\right] with {\bm{U}}=\left[\matrix{-\gamma&\gamma\cr 0.0&-\gamma}\right]. To calculate the distribution function of $A^{++}$, we need to define $\omega_{2}$ which denotes the probability that $\widehat{V}+V$ exceeds an inter-arrival time $A$. Then,

[TABLE]

Similarly, if we let $\tau_{2}$ denote the probability that $\widehat{V}+V$ exceeds $\widehat{A}$, then

[TABLE]

Further, if an inter-arrival time is following $PH$ distribution with the above representation, then following the derivation as presented in Theorem 3.1, we have

[TABLE]

Now the distribution function of $A^{++}$ may be derived as

[TABLE]

Similarly, let us denote $\widehat{A}^{++}$ and $\widetilde{A}^{++}$ as the remaining and elapsed times of the random variable $A^{++}$, then

[TABLE]

where $\lambda_{2}=\int_{0}^{\infty}x\ dF_{A^{++}}(x)$ is the mean of the random variable $A^{++}.$ For an important special case of phase-type inter-arrival time, the distribution of $A^{++}$ is also phase-type whose representation can be obtained following Theorem 3.2 as $({\bm{\bm{\alpha}}}_{2},{\bm{T}}_{2})$, where

[TABLE]

Further, with another vacation ending, let ${\bm{C}}_{n}~{}(n\geq 0)$ denote $m\times m$ matrix whose $(i,j)$th element represents the conditional probability that $n$ customers are served during an excess inter-arrival time $A^{++}$ and the service process passes to phase $j$ with the server becoming idle (after completing service of the $n$-th customer), provided at the initial instant of the previous arrival epoch there were at least $(n+1)$ customers in the system and the server was on a vacation with the service phase $i$. Then

[TABLE]

Similarly, with the vacation ending, let ${\bm{C}}_{n}^{*}$ denote the matrix of order $m\times m$ whose $(i,j)$th element represents the probability that at least $n$ customers are served during an excess inter-arrival period $A^{++}$ and the service process is in phase $j$ with the server becoming idle (after completing service of the $n$-th customer), provided at the initial instant of previous arrival epoch there were exactly $n$ customers in the system and the server was on a vacation with the service phase $i$. Then, we can define

[TABLE]

Also, with the vacation ending, let ${\bm{\Phi}}_{n}~{}(n\geq 0)$ denote $m\times m$ matrix whose $(i,j)$th element represents the conditional probability that $n$ customers are served during an elapsed excess inter-arrival time $\widetilde{A}^{++}$ and the service process passes to phase $j$ with the server becoming idle (after completing service of the $n$-th customer), provided at the initial instant of the previous arrival epoch there were at least $(n+1)$ customers in the system and the server was on a vacation with the service phase $i$. Then,

[TABLE]

Similarly, with the vacation ending, let ${\bm{\Phi}}_{n}^{*}~{}(n\geq 1)$ denote the $m\times m$ matrix whose $(i,j)$th element represents the limiting probability that at least $n$ customers have been served during an elapsed excess inter-arrival time $\widetilde{A}^{++}$ and the service process is in phase $j$ with the server becoming idle (after completing service of the $n$-th customer), provided at the previous arrival epoch the server was on a vacation with service phase $i$ with the arrival leading the system to state $n$ customers. Then we can write

[TABLE]

One may note here that the matrices ${\bm{V}}_{n}~{}(n\geq 0)$ are required to obtain other matrices by using the above relations (27)-(28). The relationships among the matrices ${\bm{C}}_{n},~{}{\bm{C}}_{n}^{*},~{}{\bm{\Phi}}_{n}$ and ${\bm{\Phi}}_{n}^{*}$ can be similarly obtained as for the matrices ${\bm{V}}_{n},~{}{\bm{V}}_{n}^{*},~{}{\bm{\Delta}}_{n}$ and ${\bm{\Delta}}_{n}^{*}$. These relationships are given below.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The matrices ${\bm{C}}_{n}$ are calculated exactly the same way as we derive the matrices ${\bm{S}}_{n}$, see Chaudhry et al. [csa15]. For more information on $MSP$, readers are referred to Bocharov [3], Albores and Tajonar [1] and Gupta and Banik [17].

4 Analysis of $GI/MSP/1/\infty$ queue with single vacation

We consider a $GI/MSP/1/\infty$ queueing system with single vacation as described above. In the following subsections we obtain steady-state distributions for this queueing system at different epochs considering $\rho<1$.

4.1 Stationary system-length distribution at pre-arrival epoch

Consider the system just before arrival epochs which are taken as embedded points. Let $t_{0},t_{1},t_{2},...$ be the time epochs at which arrivals occur and $t_{k}^{-}$ the time instant before $t_{k}$. The inter-arrival times $T_{k+1}=t_{k+1}-t_{k},~{}k=0,1,2,\ldots$ are i.i.d.r.v.’s with common distribution function $F_{A}(x)$. The state of the system at $t_{k}^{-}$ is defined as $\zeta_{k}=\{N_{t_{k}^{-}},J_{t_{k}^{-}},\xi_{t_{k}^{-}}\}$ where $N_{t_{k}^{-}}$ is the number of customers $n~{}(\geq 0)$ present in the system including the one currently in service. Whereas $J_{t_{k}^{-}}=\{j\},~{}1\leq j\leq m,$ denotes phase of the service process and $\xi_{t_{k}^{-}}=l=0$ or 1 indicates that the server is on vacation $(l=0)$ or busy $(l=1)$. In the limiting case, we define the following probabilities:

[TABLE]

where $\pi_{j,0}^{-}(n)$ represents the probability that there are $n~{}(\geq 0)$ customers in the system just prior to an arrival epoch of a customer when the server is on vacation with phase of the service process $j$. Similarly, $\pi_{j,1}^{-}(n)$ denotes the probability that there are $n~{}(\geq 0)$ customers in the system just prior to an arrival epoch of a customer when the server is in a busy (when $n\geq 1$) or dormant (when $n=0$) state with phase of the service process $j$. Let \mbox{\boldmath\pi}_{0}^{-}(n) and \mbox{\boldmath\pi}_{1}^{-}(n) be the row vectors of order $1\times m$ whose $j$-th components are $\pi_{j,0}^{-}(n)$ and $\pi_{j,1}^{-}(n),$ respectively.

Observing the state of the system at two consecutive embedded points, we have an embedded Markov chain whose state space is equivalent to $\Omega=\{(k,j,0),~{}k\geq 0,1\leq j\leq m,\}\cup\{(n,j,1),~{}n\geq 1,~{}1\leq j\leq m\}.$ Observing the system at two consecutive embedded Markov point, we have the following system of vector difference equations

[TABLE]

Multiplying (53) by $z^{n}$, summing from $n=1$ to $\infty$, after adding (52) and using the vector-generating function \mbox{\boldmath\pi}_{1}^{-*}(z)=\sum_{n=0}^{\infty}\mbox{\boldmath\pi}_{1}^{-}(n)z^{n}, we obtain

[TABLE]

leading to

[TABLE]

where {\bm{Y}}=\sum_{k=0}^{\infty}\mbox{\boldmath\pi}_{0}^{-}(k){\bm{C}}_{k+1}^{*}+\sum_{n=0}^{\infty}\mbox{\boldmath\pi}_{1}^{-}(n)\Big{(}\omega({\bm{V}}_{n+1}^{*}+\sum_{j=0}^{n}{\bm{V}}_{j})+(1-\omega)\sum_{j=0}^{n}{\bm{V}}_{j}-\sum_{i=0}^{n}{\bm{S}}_{i}\Big{)}. For further analysis, we first determine an analytic expression for each component of \mbox{\boldmath\pi}_{1}^{-*}(z). Each component $\pi_{j,1}^{-*}(z)$ defined as $\pi_{j,1}^{-*}(z)=\sum_{n=0}^{\infty}\pi_{j,1}^{-}(n)z^{n}$ of the VGF \mbox{\boldmath\pi}_{1}^{-*}(z) given in (55) being convergent in $|z|\leq 1$ implies that \mbox{\boldmath\pi}_{1}^{-*}(z) is convergent in $|z|\leq 1$. As each element of ${\bm{S}}(z^{-1})$ is a rational function, see Chaudhry et al. [10]. Therefore, each element of $det[{\bm{I}}_{m}-z{\bm{S}}(z^{-1})]$ is also a rational function and we assume that

[TABLE]

Equation (55) can be rewritten element-wise as

[TABLE]

where $\xi_{j}(z)$ is the $j$-th component of \Big{(}\sum_{j=0}^{\infty}\sum_{i=j+1}^{\infty}\mbox{\boldmath\pi}_{0}^{-}(j)(1-\omega){\bm{V}}_{i}z^{j-i+1}+\sum_{k=0}^{\infty}\mbox{\boldmath\pi}_{0}^{-}(k){\bm{C}}_{k+1}^{*}\\ +\sum_{n=0}^{\infty}\mbox{\boldmath\pi}_{1}^{-}(n)\Big{(}\omega({\bm{V}}_{n+1}^{*}+\sum_{j=0}^{n}{\bm{V}}_{j})+(1-\omega)\sum_{j=0}^{n}{\bm{V}}_{j}-\sum_{i=0}^{n}{\bm{S}}_{i}\Big{)}\Big{)}Adj[{\bm{I}}_{m}-z{\bm{S}}(z^{-1})]\varphi(z). To evaluate the vector in the numerator of equation (55), we show that the equation $det[{\bm{I}}_{m}z-{\bm{S}}(z)]=0$ has exactly $m$ roots inside the unit circle $|z|=1$, see Chaudhry et al. [9] Let these roots be $\gamma_{i}~{}(1\leq i\leq m)$. Now, consider the zeros of the function $d(z)$. Since the equation $det[{\bm{I}}_{m}z-{\bm{S}}(z)]=0$ has $m$ roots $\gamma_{i}$ inside the unit circle, the function $det[{\bm{I}}_{m}-z{\bm{S}}(z^{-1})]$ has $m$ zeros $1/\gamma_{i}$ outside the unit circle $|z|=1$. As $\pi_{j,1}^{-*}(z)$ is an analytic function of $z$ for $|z|\leq 1$, applying the partial-fraction method, we obtain

[TABLE]

where $k_{ij}$ are constants to be determined. Now, collecting the coefficient of $z^{n}$ from both sides of (57), we have

[TABLE]

Now we assume \mbox{\boldmath\pi}_{0}^{-}(0) as

[TABLE]

where $b_{1},b_{2},\ldots,b_{m}$ are $m$ arbitrary positive constants to be computed as described below. Hereafter, we substitute \mbox{\boldmath\pi}_{0}^{-}(0) from Equation (59) into the Equation (51) and obtain

[TABLE]

From (59) and (60) we are able to express \mbox{\boldmath\pi}_{0}^{-}(n)~{}(n\geq 0) in terms of the $m$ constants $(b_{1},b_{2},\ldots,b_{m})$ and $\omega$ as defined above.

Next using (58) in (52) and (53) for $n=0,2,\ldots,m-1$, we have

[TABLE]

where \mbox{\boldmath\pi}_{0}^{-}(n)~{}(n\geq 0) should be used in terms of $\omega$ and the constants as given in Equations (59) and (60). Now Equations (61) and (62) give $m^{2}$ simultaneous equations in $m(m+1)$ unknowns, $k_{ij}$’s $(1\leq i\leq m,~{}1\leq j\leq m)$ and $b_{i}~{}(1\leq i\leq m)$. The other $m$ equations can be obtained through equating the corresponding components of both sides of the vector Equation (50) as follows.

[TABLE]

where \mbox{\boldmath\pi}_{0}^{-}(k)~{}(k\geq 0) should be used in terms of $\omega$ and the constants as given in Equations (59) and (60). Finally, we have a total of $m^{2}+m=m(m+1)$ equations with $m(m+1)$ unknowns $k_{ij}$’s $(1\leq i\leq m,~{}1\leq j\leq m)$ and $b_{i}~{}(1\leq i\leq m)$. One may note here that we ignore any one component Equation of (62) for $n=m-1$ which is a redundant equation and instead we use the normalization condition given by

[TABLE]

Above normalization condition can be simplified by putting $z=1$ in Equation (57) leading to

[TABLE]

Thus solving these $m(m+1)$ equations, we get $m(m+1)$ unknowns.

4.2 Stationary system-length distribution at arbitrary epoch

We now derive explicit expressions for the steady-state queue-length distribution. Define $\pi_{i,l}(n)=$ steady-state probability that $n~{}(\geq 0)$ customers are in the system at an arbitrary epoch with server busy $(l=1)$ or on vacation $(l=0)$ and the phase of the service process is $i~{}(1\leq i\leq m).$ In other words, \mbox{\boldmath\pi}_{l}(n)=[\pi_{1,l}(n),\pi_{2,l}(n),\ldots,\pi_{m,l}(n)],~{}~{}n\geq 0;l=1 or $n\geq 0;l=0,$ at an arbitrary epoch. Here \mbox{\boldmath\pi}_{l}(n) is an arbitrary epoch stationary probability vector whose $j$-th component $\pi_{j,l}(n)~{}(1\leq j\leq m)$ is the steady-state probability that $n$ customers are in the system with server busy $(l=1)$ or vacation $(l=0).$ The classical argument based on renewal theory relates the steady-state system-length distribution at an arbitrary epoch to that at the corresponding pre-arrival epoch. Using similar results of Markov renewal theory and semi-Markov processes, see, e.g., Çinlar [14] or Lucantoni and Neuts [21], we obtain

[TABLE]

Note that since the service process is interrupted during periods in which the server is on a vacation or in an idle state, it follows that

[TABLE]

which in turn implies that

[TABLE]

for some positive constant $C$. Thus, by post multiplying the members of the previous equation by ${\bm{e}}$, we conclude that

[TABLE]

Therefore, using (72) in (71) we obtain

[TABLE]

where \rho^{{}^{\prime}}=\sum_{n=1}^{\infty}\mbox{\boldmath\pi}_{1}(n){\bm{e}} represents the probability that the server is busy. The above result (73) is useful while performing numerical calculations.

4.3 Queue-length distribution at post-departure epoch and their relation with pre-service epoch

In this subsection, we derive the probabilities for the states of the system immediately after a service completion takes place. Let \mbox{\boldmath\pi}^{+}(n)=[\pi^{+}_{1}(n),\pi^{+}_{2}(n),\ldots,\pi^{+}_{m}(n)],~{}n\geq 0, be the $1\times m$ vector whose $i$-th component $\pi^{+}_{i}(n)$ represents the post-departure epoch probability that there are $n$ customers in the queue immediately after a service completion of a customer and the server is in phase $i$. The post-departure epoch thus occurs immediately after the server has either reduced the queue or become idle. Hence, using level-crossing arguments given in Chaudhry and Templeton [8, p. 299], we have

[TABLE]

It may be noted that \sum_{n=1}^{\infty}\mbox{\boldmath\pi}_{1}(n){\bm{L}}_{1}{\bm{e}}=\mu^{*}\rho^{{}^{\prime}}, which represents the departure rate when the server is busy.

Let \mbox{\boldmath\pi}^{s-}(n)=[\pi^{s-}_{1}(n),\pi^{s-}_{2}(n),\ldots,\pi^{s-}_{m}(n)],~{}n\geq 1, be the $1\times m$ vector whose $i$-th component $\pi^{s-}_{i}(n)$ represents the pre-service epoch probability that there are $n$ customers in the queue immediately before a service of a customer takes place and the server is in phase $i$. The argument used to find post-departure epoch probabilities may be based on the distribution of the probabilities for the system at pre-service epoch of a customer, the instant in time immediately before a real service of a customer starts. Using the above arguments, we obtain the following result.

[TABLE]

5 Performance measures

As state probabilities at various epochs are known, performance measures can be easily obtained. The average number of customers in the system (queue) at an arbitrary epoch are given by

[TABLE]

5.1 **Waiting-time analysis **

In this section, we obtain the LST of waiting-time distribution of a customer who is accepted in the system. Let \mbox{\boldmath\phi}_{k}(\theta) be the LST of the probability that $k$ customers will be served within a time $x$ and the service process upon completion of service passes to phase $j$, provided $k$ customers were in the system and the service process was in phase $i$ at the beginning of service. Since the probability that the service of a customer is completed in the interval $(x,x+dx]$ is given by the matrix $e^{{\bm{L}}_{0}x}{\bm{L}}_{1}dx$ and the total service time of $k$ customers is the sum of their service times, {\mbox{\boldmath\phi}}_{1}(s), the LST of service time with corresponding phase change, is given by

[TABLE]

Further, we also need the LST of remaining vacation-time. It is given by

[TABLE]

Let $W_{s}^{*}(s)$ denote the LST of the actual waiting time distribution of an arbitrary customer in system. The LST of the waiting-time distribution in system can be derived as follows:

[TABLE]

where bpi_{0}^{-*}(z)=\sum_{j=0}^{\infty}\mbox{\boldmath\pi}_{0}^{-}(j)z^{j},~{}|z|\leq 1, defined similarly as \mbox{\boldmath\pi}_{1}^{-*}. One can find mean waiting time in the system from Equation (80) by differentiating it and putting $s=0$. It is given by

[TABLE]

From the Little’s law, we can also get mean sojourn time as $W_{s}(LL)=\frac{L_{s}}{\lambda}$ which may serve one of the verifications while performing numerical computation.

5.2 **Expected length of busy and idle periods **

Since for this system, in the limiting case, the proportions of times the server is busy and idle are $\rho^{{}^{\prime}}$ and $1-\rho^{{}^{\prime}}$, respectively, we have

[TABLE]

where $B$ and $I$ are random variables denoting the lengths of busy and idle periods, respectively. We first discuss the mean busy period $E(B)$, which is comparatively easy to evaluate. Let $N_{q}(t)$ denote the number of customers in system at time $t$ and $\xi_{q}(t)$ be the state of the server, i.e., busy $(=1)$ or idle $(=0)$. $\{N_{q}(t),\xi_{q}(t)\}$ enters the set of busy states, $\Upsilon\equiv\{(0,1),(1,1),(2,1),\ldots\}$ at the termination of an idle period. The conditional probability that $\{N_{q}(t),\xi_{q}(t)\}$ enters $(0,0),$ given that $\{N_{q}(t),\xi_{q}(t)\}$ enters $\Upsilon$, is therefore C\mbox{\boldmath\pi}^{+}(i){\bm{e}},~{}i\geq 0, where C=\frac{1}{\mbox{\boldmath\pi}^{+}(0){\bm{e}}}. Now $\{N_{q}(t),\xi_{q}(t)\}$ enters $(i,1),~{}i\geq 0,$ irrespective of customers’ arrival during a service time, which may happen in expected time $E(S)$. Thus

[TABLE]

Substituting $E(B)$ from Equation (83) in Equation (82), we obtain

[TABLE]

6 Approximation for system-length distributions based on one root

In this context one may note that the heavy-traffic approximation was first investigated by J.F. C. Kingman (see [19]) who showed that when the utilisation parameter ($\rho$) of an $M/M/1$ queue is near 1 a scaled version of the queue length process can be accurately approximated by a reflected Brownian motion. In this direction one can get several approximate results such as the tail probabilities of the queue-length distribution at a pre-arrival epoch, heavy- or light-traffic behaviour of the queue-length distributions based on the real root of $det[{\bm{I}}_{m}-z{\bm{S}}(z^{-1})]=0$ which is closest to 1 and outside $|z|\leq 1$. The existence of such a root has been discovered since long time in the literature, e.g., Feller ([16], pg. 276-277) calculated the tail probabilities using a single root of the denominator which is smallest root in absolute value. Also Chaudhry et al. [11] have given a formal proof of the existence of such a root. One may note that it is not difficult to calculate this root numerically. In this context it is worth mentioning that sometimes an approximate value of this root may be used to get desired queue-length distributions and this approximate value may be obtained in the following way. We investigate the approximate root inside $|z|\leq 1$ by expanding the matrix of the left-hand side of the characteristic equation $det[{\bm{I}}_{m}z-{\bm{S}}(z)]=0$ in powers of $\rho$ as

[TABLE]

where $\overline{{\bm{H}}}$ is some unknown matrix and ${\bm{S}}^{{}^{\prime}}(.),~{}{\bm{S}}^{{}^{\prime\prime}}(.)$ are the successive differentiation of ${\bm{S}}(.)$ of order one and two, respectively. Also one may note that in Equation (85), as usual, $o(x)$ represents a function of $x$ with the property that $\frac{o(x)}{x}\rightarrow 0$ as $x\rightarrow 0.$ Multiplying the right-hand side of the above Equation (85) by the vector \overline{\mbox{\boldmath\pi}} from left and the vector ${\bm{e}}$ from right, we may write the characteristic equation as follows

[TABLE]

where \overline{c}=\overline{\mbox{\boldmath\pi}}\overline{{\bm{H}}}{\bm{e}} and is a finite constant. Now an approximate value of the root for the above described three cases is obtained as follows.

• •

Light-traffic case: Applying $\rho\rightarrow 0+$ in Equation (86) gives

[TABLE]

which gives an approximate value of this root.

• •

Heavy-traffic case: Replacing $\rho$ by $(1-\rho)$ in Equation (86) and applying the condition $\rho\rightarrow 1-$, we obtain

[TABLE]

Solving (88) we get the desired value.

• •

Tail probabilities at a pre-arrival epoch: If $\rho_{1}$ denotes any arbitrary offered load numerically very close to $\rho$ then we replace $\rho$ by $(\rho-\rho_{1})$ in Equation (86) and applying the condition $\rho\rightarrow\rho_{1}$ to get

[TABLE]

Hence, we can obtain the desired root, say $z_{1}$, by solving the Equation (89) for $z$.

Finally, it may be noted that once we obtain an approximate value of a root, we can obtain the exact root through various numerical methods. The equation $det[{\bm{I}}_{m}z-{\bm{S}}(z)]=0$ may be used to find the original root which is closest to 1 in all the above described cases. For the sake of completeness we present below the procedure to calculate tail probabilities at a pre-arrival epoch based on this one root. To get the tail probabilities, assume

[TABLE]

where $n_{\epsilon}$ is chosen as the smallest integer such that $|(\pi_{j,1}^{-}(n)-p_{j,1}^{a1}(n))/\pi_{j,1}^{-}(n)|<\epsilon,$ i.e., $|1-\frac{p_{j,1}^{a1}(n)}{\pi_{j,1}^{-}(n)}|<\epsilon.$ But, since the probability $p_{j,1}^{a1}(n)$ follows a geometric distribution with common ratio $z_{1}$, it is better to choose $n_{\epsilon}$ such that $|\frac{\pi_{j,1}^{-}(n)}{z_{1}\pi_{j,1}^{-}(n-1)}-1|<\epsilon.$ The approximation gets better if more than one root, in ascending order of magnitude, is used. It should however be mentioned that those roots that occur in complex-conjugate pairs should be used in pairs. Thus, the tail probabilities using three roots can be approximated by

[TABLE]

where $z_{i}~{}(i=1,2,3)$ are the roots in ascending order of magnitude and $n_{\epsilon}^{1}$ may be chosen by $|1-\frac{p_{j,1}^{a3}(n)}{\pi_{j,1}^{-}(n)}|<\epsilon,~{}n>n_{\epsilon}^{1}.$ Similar procedure may be adopted to calculate queue-length distributions for the cases of light- and heavy-traffics. It may be remarked here that this root can also be obtained accurately by simply using high precision of the software packages mentioned earlier. Similar way we can compute the waiting time distribution based on a few number of roots in case of light- and heavy-traffic.

7 Numerical results and discussion

To demonstrate the applicability of the results obtained in the previous sections, some numerical results have been presented in two self explanatory tables. At the bottom of the tables, several performance measures are given. Since various distributions can be either represented or approximated by $PH$-distribution, we take inter-arrival time distribution to be of $PH$-type having the representation $({\bm{\alpha}},{\bm{T}})$, where ${\bm{\alpha}}$ and ${\bm{T}}$ are of dimension $\eta$. Then ${\bm{S}}(z)$ can be derived as follows using the procedure adopted in [10].

[TABLE]

with ${\bm{L}}(z)\oplus{\bm{T}}=({\bm{L}}(z)\otimes{\bm{I}}_{\nu})+({\bm{I}}_{m}\otimes{\bm{T}})$,  where $\oplus$ and $\otimes$ are used for Kronecker product and sum, respectively. For the derivation of ${\bm{S}}(z)$, see Chaudhry et al [10]. Knowing that each element of ${\bm{L}}(z)$ is a polynomial in $z$, each element of ${\bm{L}}(z)\oplus{\bm{T}}$ is also a polynomial in $z$ and hence the determinant of $({\bm{L}}(z)\oplus{\bm{T}})$ is a rational function in $z$. Thus, from the above expression for ${\bm{S}}(z)$, we can immediately say that each element of ${\bm{S}}(z)$ is a rational function in $z$ with the same denominator. One may note that in case the degree of the polynomials in each element of ${\bm{S}}(z)$ is very high, it may be difficult or time consuming to calculate the roots of the characteristic equation

[TABLE]

This difficulty may be minimized by calculating ${\bm{S}}(z)$ in the following way.

[TABLE]

where ${\bm{S}}_{n}$ may be obtained as proposed in Lucantoni [20].

We have carried out extensive numerical work based on the procedure discussed in this paper by considering different service matrices $MSP~{}({\bm{L}}_{0},{\bm{L}}_{1})$ and phase-type inter-arrival time distribution $PH({\bm{\alpha}},{\bm{T}})$. All the calculations were performed on a PC having Intel(R) Core 2 Duo processor @1.65 GHz with 8 GB DDR2 RAM using MSPLE 18. Further, though all the numerical results were carried out in high precision, they are reported here in 6 decimal places due to lack of space.

In Table 1, we have presented various epoch probabilities for a $PH/MSP/1/\infty$ queue with exponential single vacation using our method described in this paper. Vacation time is following exponential distribution with average number of vacations per unit of time is $\gamma=1.8$. Inter-arrival time is $PH$-type and its representation is given by

[TABLE]

,

[TABLE]

with $\lambda=0.259558$. The $MSP$ matrices as

[TABLE]

[TABLE]

with stationary mean service rate $\mu^{*}=1.121972$, lag-1 correlation coefficient $0.618173$ between successive service times and \overline{\mbox{\boldmath\pi}}=\left[\matrix{0.264645&0.253046&0.254961&0.227348}\right] so that $\rho=\lambda/(\mu^{*})=0.231341.$ To calculate system-length distribution we need to calculate the roots of

[TABLE]

where $m=4$ as given above. Here ${\bm{S}}(z)$ may be obtained by Equation (LABEL:numa_1) or (94) for $N=70$, see Equation (92). The $m=4$ roots of (95) inside $|z|<1$ are evaluated. The corresponding $k_{ij}~{}(1\leq i\leq 4,~{}1\leq j\leq 4)$ and $b_{i}~{}(1\leq i\leq 4)$ values are calculated using the procedure described in Section 4.1, see Appendix A. Now using Equation (58), (59) and (60), one can obtain system-length distribution at pre-arrival epoch and after that using relations (66)-(69) the arbitrary epoch probabilities may be derived, see Table 1.

Table 1: System-length distributions at pre-arrival and arbitrary epoch.

[TABLE]

[TABLE]

It may be noted that in the above numerical experiment, we can find the conditional probability that the server is busy in phase $i,~{}i=1,2.$ It is given by

[TABLE]

which matches with \overline{\mbox{\boldmath\pi}} up to almost 2 decimal places. As shown in the above table, Little’s law is satisfied up to two digits. These, to some extent, support the validity of our analytical as well as numerical results.

Next one may note here that this queueing model deals with generally distributed inter-arrival time distribution. Theoretically, it is possible to approximate any non-negative distribution arbitrarily closely by a $PH$-type distribution, see Bobbio and Telek [5], Bobbio et al. [6] and references therein. Using these methods available in the literature, one may estimate a general inter-arrival time distribution by a $PH$-type distribution and calculate stationary distribution at various epochs as stated above. As a demonstration, $PH$ approximation of inter-arrival time distributions which are not phase-type, we consider the following example.

Similar to the above tables, in Table 2, we have presented various epoch probabilities for $LN/MAP/1/\infty$ queue with exponential single vacation, where $LN$ stands for log-normal inter-arrival time distribution. Vacation time is exponential with the stationary mean vacation rate $\gamma=1.7$. Inter-arrival time is $LN$-type and its probability density function and distribution functions are given by $f_{A}(x)=\frac{1}{x\alpha\sqrt{2\pi}}e^{-\frac{(\ln(x)-\beta)^{2}}{2\alpha^{2}}},~{}\alpha=1.04,~{}\beta=0.215,~{}x>0,$ and $F_{A}(x)=\int_{0}^{x}f_{A}(u)du,~{}x>0,$ with $\lambda=0.469635$, respectively. The $MAP$ matrices have 4 phases in this case and their representation is given by

[TABLE]

[TABLE]

with stationary mean service rate $\mu^{*}=1.112288$, lag-1 correlation coefficient 0.179796 between successive service times and \overline{\mbox{\boldmath\pi}}=\left[\matrix{0.264645&0.253046&0.254961&0.227348}\right] so that $\rho=\lambda/(\mu^{*})=0.422224.$ . In the following we state a procedure to discuss a $PH$ approximation of the $LN$-type inter-arrival time. As discussed by Bobbio and Telek [5], Bobbio et al. [6] and Pulungan [25], we assume a minimal acyclic phase-type canonical representation of the log-normal inter-arrival time distribution as follows:

[TABLE]

and

[TABLE]

where $\overline{t}_{1},~{}\overline{t}_{2},~{}\overline{t}_{3}>0,$ $\overline{t}_{1}\leq\overline{t}_{2}\leq\overline{t}_{3},$ $\alpha_{1},~{}\alpha_{2},~{}\alpha_{3}\geq 0,$ and ${\bm{\alpha}}{\bm{e}}_{\eta}=1.0,~{}\eta=3.$ The probability density and distribution function of this phase-type representation is given by

[TABLE]

where ${{\bm{T}}}^{0}$ is non-negative vector and satisfies ${\bm{T}}{\bm{e}}_{\eta}+{{\bm{T}}}^{0}={\bm{0}}$ and $\eta=3$. The moments of this acyclic phase-type representation may obtained by the following formulae:

[TABLE]

Using the log-normal inter-arrival time probability density function we can calculate moments by the following formulae:

[TABLE]

The first four moments of the original log-normal distribution have been targeted to match with corresponding moment of the acyclic phase-type representation. This multi-objective problem can be reduced to a single objective problem by assigning weights to objectives. The non-linear programming (NLP) problem can be formulated as:

[TABLE]

where $w_{i}$’s ($i=1,~{}2,~{}3,~{}4$) are the weights. We always set $w_{i}=1,~{}i=1,~{}2,~{}3,~{}4$ except in cases when numerical values of moments are very high. For example when numerical value of fourth moment is very high we set $w_{4}=0$ and NLP problem gives better result by assigning $w_{4}=0$ in lieu of $w_{4}=1.$ One may note that after solving the above NLP we can completely specify an approximate acyclic phase-type representation of the given weibull distribution. For better approximation one may increase the order of acyclic phase-type distribution and in that case obtaining a solution to the above NLP may pose some problem due to increase in the number of variables. To overcome this problem we supply eigenvalues of the matrix ${\bm{T}}$, i.e., $\overline{t}_{1}$ $\overline{t}_{2}$ and $\overline{t}_{3}$ using the following procedure. We calculate 20 moments from $f_{A}(x)$ to construct an approximate $f_{A}^{*}(s)$ through the Padé approximation [2/3] as:

[TABLE]

where $[2/3]$ stands for a rational function with degree of numerator polynomial 2 and degree of denominator polynomial 3 in the variable $s$. Now after equating denominator of (104) equal to zero we obtain three roots as $-0.017943,~{}-0.112322$ and $-0.548702.$ Hereafter, we assign eigenvalues of the matrix ${\bm{T}}$ as $\overline{t}_{1}=0.017943,~{}\overline{t}_{2}=0.112322$ and $\overline{t}_{3}=0.548702.$ After this, our job is to solve the above NLP and obtain the vector ${\bm{\alpha}}.$ Solution of above NLP with assignment of weights as $w_{i}=1,~{}i=1,~{}2,~{}3,~{}4$ gives $\alpha_{1}=0.000024,~{}\alpha_{2}=0.034291$ and $\alpha_{3}=0.965685.$ Therefore, we are able to specify the vector ${\bm{\alpha}}$ as well as the matrix ${\bm{T}}$ corresponding to an approximate phase-type distribution of the log-normal inter-arrival time distribution. We calculate $\lambda$ using these ${\bm{\alpha}}$ and ${\bm{T}}$ values and found $\lambda=0.469635$ which is exactly the same as the one calculated using the original log-normal density function given above. In the following we present graphs of the density and distribution functions (see Figure 1 and 2, respectively) for the original log-normal inter-arrival time and the corresponding phase-type approximation.

Now we consider these approximate phase-type representations for the corresponding inter-arrival and vacation time distributions and using similar procedure as described in Table 1, we have computed stationary system length distribution at pre-arrival and arbitrary epochs, see Table 2.

Table 2: System-length distributions at pre-arrival and arbitrary epoch.

[TABLE]

[TABLE]

It may be noted that in the above numerical experiment, we can find the conditional probability that the server is busy in phase $i,~{}i=1,2.$ It is given by

[TABLE]

which is \overline{\mbox{\boldmath\pi}} up to almost 2 decimal places as was anticipated. Also, Little’s law is satisfied up to two digits after rounding off. These facts confirm our analytical as well as numerical results.

8 Conclusions and future scope

In this paper, we have successfully analyzed the $PH/MSP/1/\infty$ queue with single exponential vacation. We have suggested a procedure to obtain the steady-state distributions of the number of customers in the system at pre-arrival, arbitrary and post-departure epochs. Similar kind of analysis that is described in this paper may work for the corresponding queueing systems under multiple exponential vacations of the server or for batch arrival or batch service queues, i.e., $PH^{[X]}/MSP/1/\infty$ queue or $PH/MSP^{(a,b)}/1/\infty$ queue with exponential single or multiple vacations. One may be interested in analyzing the same queueing model with different type of vacation policies, e.g., multiple adaptive vacation(s) and working vacation(s). Another area of interest may be to find the approximations for the tail of the waiting-time distribution a swell as an approximation for the waiting-time distribution in cases of heavy- and light-traffic. These problems are left for future investigations.

Appendix A

The roots used in Table 1 are given as follows: $\gamma_{1}=-0.010443,~{}\gamma_{2}=0.853818,~{}\gamma_{3}=0.160075-0.064040\mathrm{i},~{}\gamma_{4}=0.160075+0.064040\mathrm{i}.$ The corresponding $k_{ij}~{}(1\leq i\leq 4,~{}1\leq j\leq 4)$ and $b_{i}~{}(1\leq i\leq 4)$ values are as follows: $k_{1,1}=0.250531,~{}k_{1,2}=0.009814,~{}k_{1,3}=0.017357,~{}k_{1,4}=0.093024,~{}k_{2,1}=0.001666,~{}k_{2,2}=0.009563,~{}k_{2,3}=0.007588,~{}k_{2,4}=0.002873,~{}k_{3,1}=0.083053+0.092411\mathrm{i},~{}k_{3,2}=-0.008081-0.013702\mathrm{i},~{}k_{3,3}=-0.009387-0.050870\mathrm{i},~{}k_{3,4}=-0.047920-0.161839\mathrm{i},~{}k_{4,1}=0.083053-0.092411\mathrm{i},~{}k_{4,2}=-0.008081+0.013702\mathrm{i},~{}k_{4,3}=-0.009387+0.050870\mathrm{i},~{}k_{4,4}=-0.047920+0.161839\mathrm{i}$ and $b_{1}=0.394155,~{}b_{2}=0.001935,~{}b_{3}=0.007826,~{}b4=0.001078$.

The roots used in Table 2 are given below. $\gamma_{1}=0.194296,~{}\gamma_{2}=0.831411,~{}\gamma_{3}=0.328896-0.025835\mathrm{i}$ and $\gamma_{4}=0.328896+0.025835\mathrm{i}.$ The corresponding $k_{ij}~{}(1\leq i\leq 4,~{}1\leq j\leq 4)$ and $b_{i}~{}(1\leq i\leq 4)$ values are as follows: $k_{1,1}=-0.205961,~{}k_{1,2}=0.001301,~{}k_{1,3}=-0.036359,~{}k_{1,4}=0.147868,~{}k_{2,1}=0.005287,~{}k_{2,2}=0.024221,~{}k_{2,3}=0.015824,~{}k_{2,4}=0.010081,~{}k_{3,1}=0.234467-0.549961\mathrm{i},~{}k_{3,2}=-0.012014-0.060543\mathrm{i},~{}k_{3,3}=0.025634+0.033050\mathrm{i},~{}k_{3,4}=-0.078568+0.512859\mathrm{i},~{}k_{4,1}=0.234467+0.549961\mathrm{i},~{}k_{4,2}=-0.012014+0.060543\mathrm{i},~{}k_{4,3}=0.025634-0.033050\mathrm{i},~{}k_{4,4}=-0.078568-0.512859\mathrm{i}$ and $b_{1}=0.190353,~{}b_{2}=0.001284,~{}b_{3}=0.028939,~{}b_{4}=0.000851$.

Acknowledgement The second author was supported partially by NSERC under research grant number RGPIN-2014-06604.