Performance Analysis of Multi-Server Queueing System Operating under Control of a Random Environment

1. Introduction

Since the early 1900th the Erlang multi-server queueing systems with losses ( $B$ -model or $$M/M/N/0$$ system) and with an infinite size buffer ( $C$ -model or $$M/M/N$$ system) provided good mathematical tools for capacity planning and performance evaluation in the classic telephone networks for many years. Good quality of the loss probability forecasting in real world networks based on the formulas obtained for the $$M/M/N/0$$ system and the delay prediction based on the formula obtained for the $$M/M/N$$ system was a rather surprising because the requirement that inter-arrival and service times have an exponential distribution, which is imposed in the $$M/M/N/0$$ and $$M/M/N$$ models, seems to be too strict. The interest of mathematicians to the fact of good matching of the calculated under debatable assumptions characteristics to their measured value in real world systems have lead to the following two results.

By efforts of many mathematicians (A. Ya. Khinchin and B. I. Grigelionis first of all), it was proved that the superposition of a large number of independent flows having uniformly small intensity approaches to the stationary Poisson input when the number of the superposed inputs tends to infinity. It explains the fact that the flows in classic telephone networks (where flows are composed by small individual flows from independent subscribers) have the exponentially distributed inter-arrival times.

Concerning the service time distribution, situation was more complicated. The real-life measurements have shown that the service (conversation) time can not be well approximated by means of the exponentially distributed random variable. So, due to the good matching of results obtained for the $$M/M/N/0$$ queue performance characteristics to characteristics of real systems modelled by such a queue, the hypothesis has arisen that the stationary state distribution in the $$M/M/N/0$$ queue is the same as the one in the $$M/G/N/0$$ queue conditional that the average service times in both models coincide. This property was called the invariant (or insensitivity) property of the model with respect to the service time distribution. The work [10] by B. A. Sevastjanov is the first one where this property was proven strictly.

So, the question why the Erlang's models give very good results for a practice was highlighted. The special books containing the tables for a loss probability under the given values of the number $N$ of channels and intensities of the input and service exist. Different design problems (for a fixed value of permissible loss probability, to find the maximal intensity of the flow, which can be served by the line consisting of a fixed number of channels under the fixed average service time, or to find the necessary number of channels sufficient for transmission of the flow with a fixed intensity, etc) are solved by means of these tables.

However, the flows in the modern telecommunication networks have lost the nice properties of their predecessors in the old classic networks. In opposite to the stationary Poisson input (stationary ordinary input with no aftereffect), the modern real life flows are non-stationary, group and correlated. The $$BMAP$$ (Batch Markovian Arrival Process) arrival process was introduced as a versatile Markovian point process ( $$VMPP$$ ) by M.F. Neuts in the 70th. The original development of $$VMPP$$ contained extensive notations; however these notations were simplified greatly in [7] and ever since this process bears the name $$BMAP$$ . The class of $$BMAP$$ s includes many input flows considered previously, such as stationary Poisson ( $M$ ), Erlangian ( $E_k$ ), Hyper-Markovian ( $HM$ ), Phase-Type ( $PH$ ), Interrupted Poisson Process ( $IPP$ ), Markov Modulated Poisson Process ( $MMPP$ ). Generally speaking, the $$BMAP$$ is correlated, so it is ideal to model correlated and (or) bursty traffic in modern telecommunication networks.

As it was mentioned above, the question why the inter-arrival times in the classical networks have the exponential distribution was answered in literature. However, Erlang's assumption that the service time distribution has the exponential distribution is not supported by the real networks measurements. In the case of the $$M/M/N/0$$ system, good fitting of performance measures of this system with the respective measures of real world systems is easy explained by Sevastjanov's result. But in the case of the $$M/M/N$$ system, it was necessary to generalize results by Erlang to the cases of another, than exponential, service time distributions. This work was started by Erlang who offered so called Erlang's distribution. He introduced Erlangian of order $k$ distribution as a distribution of a sum of $k$ independent identically exponentially distributed random variables (phases). Further, so called phase type ( $PH$ ) distribution was introduced into consideration as the straightforward generalization of Erlangian distribution, see, e.g., [8]. $PH$ distribution includes as the special cases the exponential, Erlangian, Hyper-exponential, Coxian distributions. In our chapter we assume that service times at the fixed operation mode of the system have $PH$ distribution.

It follows from discussion above that it is interesting to extend investigation of Erlang's models to the case of the $$BMAP$$ input and $PH$ type service process. This work was started by M. Combe in [1] and V. Klimenok in [5] where the $$BMAP/M/N$$ and $$BMAP/M/N/0$$ models, respectively, were investigated.

In paper [2], we investigated the $$BMAP/PH/N/0$$ model having no buffer. It was shown there that the stationary distribution of the system states essentially depends on the shape of the service time distribution and so Sevastjanov's invariant property does not hold true in the case of the general $$BMAP$$ arrival process. Here we analyze the $$BMAP/PH/N/L$$ system with a finite buffer and the $$BMAP/PH/N$$ system with infinite buffer. Simultaneously, we make one more essential generalization of the model under study. Motivation of this generalization is as follows.

Even if one will use such general models of the arrival and service process as the $$BMAP$$ and $PH$ , he may fail in application to practical systems. The reason is the following. Assumption that the input flow is described by the $$BMAP$$ allows to take into consideration a burstiness, an effect of correlation in the arrival process and variation of inter-arrival times. Assumption that the service process is described by the $PH$ distribution allows to take into consideration variation of service times. But the $$BMAP$$ arrival process and $PH$ service process are assumed to be stationary and independent of each other within the borders of the models of $$BMAP/PH/N/L$$ type, $$0\le L\le \infty \mathrm{.}$$ While in many real world systems the input and service processes are not absolutely stable and may be mutually dependent. They may be influenced by some external factors, e.g., the different level of the noise in the transmission channel, hardware degradation and recovering, change of the distance by a mobile user from the base station, parallel transmission of high priority information, etc. Information transmission channel modeled by means of the $$BMAP/PH/N/L$$ queueing system can be a part of complex communication network. The rest of the network may essentially vary characteristics of the arrival and service process in this system by means of: (i) changing the bandwidth of the channel (due to reliability factors or the needs to provide good quality of service in another parts of the network when congestion occurs); (ii) changing the mean arrival rate due breakdowns, overflow or underflow of alternative information transmission channels. Thus, to get the mathematical tool for adequate modeling such information transmission channels, more complicated queues than the $$BMAP/PH/N/L$$ queueing system should be analyzed. These queues, in addition to the account of complicated internal structure of the arrival and service processes by means of considering the $$BMAP$$ and $PH$ must take into account the influence of random external factors. In some extent, it can be done by means of analyzing the models of queues operating in a random environment. Such an analysis is the topic of this chapter.

Importance of investigation of the queues operating in a random environment ( $RE$ ) drastically increased in the last years due to the following reason. The flows of information in the modern communication networks are essentially heterogeneous. Some types of information are very sensitive with respect to a delay and an jitter but tolerant with respect to losses. Another ones are tolerant with respect to the delay but very sensitive with respect to the loss of the packets. So, different schemes of the dynamic bandwidth sharing among these types exist and are developing. They assume that, in the case of congestion, transmission of the delay tolerant flows is temporarily postponed to provide better conditions for transmission of the delay sensitive flows. Analysis of such schemes requires the probabilistic analysis of the multi-dimensional processes describing transmission process of the different flows. This analysis is often impossible due to the mathematical complexity. In such a case, it is reasonable to decompose a simultaneous consideration of all flows to separate analysis of the processes of transmission of the delay sensitive and the delay tolerant flows. To this end, we model transmission of the delay sensitive flows in terms of the queues with the controlled service or (and) arrival rate where the service or the arrival rate can be changed depending on the queue length or the waiting time. Redistribution of a bandwidth to avoid congestion for the delay sensitive flows causes a variation, at random moments, of an available bandwidth for the delay tolerant flows. Correspondingly, the queues operating in a random environment naturally arise as the mathematical model for the delay tolerant flows transmission.

Mention that the $$BMAP/PH/N/0$$ model operating in the $RE$ was recently investigated in [4]. Short overview of the recent research of queues operating in the RE can be found there. In this chapter, we consider the models $$BMAP/PH/N/L$$ and $$BMAP/PH/N$$ operating in the $RE$

2. The Mathematical Model

We consider the queueing system having $N$ identical servers. The system behavior depends on the state of the stochastic process (random environment) $r_{t}t\ge 0$ which is assumed to be an irreducible continuous time Markov chain with the state space $$\{1\dots R\}$$ $$R\ge 2$$ and the infinitesimal generator $$Q$$

The input flow into the system is the following modification of the $$BMAP$$ . In this input flow, the arrival of batches is directed by the process ${\nu }_{t}t\ge 0$ (the underlying process) with the state space $$\{\mathrm{0,1,}\dots W\}$$ Under the fixed state $$r$$ of the $RE$ this process behaves as an irreducible continuous time Markov chain. Intensities of transitions of the chain ${\nu }_{t}t\ge 0$ which are accompanied by arrival of $$k$$ -size batch, are described by the matrices $$D_k^{(r)}k\ge 0$$ $r=\overline{1R}$ with the generating function $$D^{(r)}(z)={\displaystyle \sum _{k=0}^{\infty }{D}_k^{(r)}}{z}^k\left|z\right|\le 1$$ The matrix $$D^{(r)}(1)$$ is an irreducible generator for all $r=\overline{1R}$ Under the fixed state $$r$$ of the random environment, the average intensity $${\lambda }^{(r)}$$ (fundamental rate) of the $$BMAP$$ is defined as $${\lambda }^{(r)}={\theta }^{(r)}{\left(D^{(r)}(z)\right)}^{\prime }|{}_{z=1}e$$ and the intensity $${\lambda }_b^{(r)}$$ of batch arrivals is defined as $${\lambda }_b^{(r)}={\theta }^{(r)}\left(-D_0^{(r)}\right)e$$ Here the row vector $${\theta }^{(r)}$$ is the solution to the equations $${\theta }^{(r)}{D}^{(r)}(1)=0{\theta }^{(r)}e=1$$ $$e$$ is a column vector of appropriate size consisting of 1's. The variation coefficient $$c_{var}^{(r)}$$ of intervals between batch arrivals is given by

$${\left(c_{var}^{(r)}\right)}^2=2{\lambda }_b^{(r)}{\theta }^{(r)}{(-D_0^{(r)})}^{-1}e-1$$ E1

while the correlation coefficient $$c_{cor}^{(r)}$$ of intervals between successive batch arrivals is calculated as

$$c_{cor}^{(r)}=\left({\lambda }_b^{(r)}{\theta }^{(r)}{(-D_0^{(r)})}^{-1}(D^{(r)}(1)-D_0^{(r)}){(-D_0^{(r)})}^{-1}e-1\right)/{\left(c_{var}^{(r)}\right)}^2$$ E2

At the epochs of the process $r_{t}t\ge 0$ transitions, the state of the process ${\nu }_{t}t\ge 0$ is not changed, but the intensities of its transitions are immediately changed.

The service process is defined by the modification of the $PH$ -type service time distribution. Service time is interpreted as the time until the irreducible continuous time Markov chain $m_{t}t\ge 0$ with the state space $$\{\mathrm{0,1,}\dots M+1\}$$ reaches the absorbing state $M+1$ Under the fixed value $$r$$ of the random environment, transitions of the chain $m_{t}t\ge 0$ within the state space $$\{1\dots M\}$$ are defined by an irreducible sub-generator $S^{(r)}$ while the intensities of transition into the absorbing state are defined by the vector $S_0^{(r)}=-S^{(r)}e$ At the service beginning epoch, the state of the process $m_{t}t\ge 0$ is chosen according to the probabilistic row vector $${\beta }^{(r)}r=\overline{1R}$$ It is assumed that the state of the process $m_{t}t\ge 0$ is not changed at the epoch of the process $r_{t}t\ge 0$ transitions. Just the exponentially distributed sojourn time of the process $m_{t}t\ge 0$ in the current state is re-started with a new intensity defined by the sub-generator corresponding to the new state of the random environment $r_{t}t\ge 0$

The system under consideration has $$L0\le L\le \infty$$ waiting positions. In the case of an infinite buffer ( $$L=\infty$$ ) all customers are always admitted to the system. In the case of a finite $$L$$ the system behaves as follows. If the system has all servers being busy at a batch arrival epoch, the batch looks for the available waiting position, and occupies it in case of success. If the system has all servers and all waiting positions being busy, the batch leaves the system forever and is considered to be lost. Due to a possibility of the batch arrivals, it can occur that there are free servers or waiting positions in the system at an arrival epoch, however the number of these positions is less than the number of the customers in an arriving batch. In such situation the acceptance of the customers to the system is realized according to the partial admission ( $PA$ ) discipline (only a part of the batch corresponding to the number of free servers is allowed to enter the system while the rest of the batch is lost), the complete rejection ( $CR$ ) discipline (a whole batch leaves the system if the number of free servers is less than the number of customers in the batch), complete admission ( $CA$ ) discipline (a part of the batch corresponding to a number of free servers starts the service immediately while the rest of the batch waits for a service in the system in some special waiting space). All these disciplines are popular in the real life systems and got a lot of attention in the literature. Here, we consider all these disciplines.

Our aim is to calculate the stationary state distribution and main performance measures of the described queueing model.

For the use in the sequel, let us introduce the following notation:

$e_n(0_n)$ is a column (row) vector of size $$n$$ consisting of 1's (0's). Suffix may be omitted if the dimension of the vector is clear from context; $$I\left({\text{O}}_{}\right)$$ is an identity (zero) matrix of appropriate dimension (when needed the dimension of this matrix is identified with a suffix); $diag\left\{a_{k}k=\overline{1K}\right\}$ is a diagonal matrix with diagonal entries or blocks $a_k;$ $$\otimes$$ and $$\oplus$$ are symbols of the Kronecker product and sum of matrices; ${\Omega }^{\otimes l}=\underset{l}{\underbrace{\Omega \otimes \dots \otimes \Omega }}l\ge 1{\Omega }^{\otimes 0}=1$ ${\Omega }^{\oplus l}={\displaystyle \sum _{m=0}^{l-1}{I}_{n^m}\otimes \Omega \otimes I_{n^{l-m-1}}}l\ge 1$ where $n$ is the dimension of square matrix $\Omega$ ;

$$D(z)={\displaystyle \sum _{k=0}^{\infty }diag\left\{D_k^{(r)}r=\overline{1R}\right\}}{z}^k$$ E3

$${\mathcal{D}}_k^{(n)}=diag\left\{D_k^{(r)}\otimes I_{M^n}r=\overline{1R}\right\}n=\overline{0N}k\ge 0$$ E4

$${\mathcal{D}}^{(n)}(z)={\displaystyle \sum _{k=0}^{\infty }{\mathcal{D}}_k^{(n)}{z}^k}n=\overline{0N};$$ E5

$${ℬ}_l^{(n)}=diag\left\{I_{\overline{W}}\otimes I_{M^n}\otimes {\left({\beta }^{(r)}\right)}^{\otimes l}r=\overline{1R}\right\}n=\overline{1N}\overline{W}=W+1$$ E6

$${\mathcal{S}}^{(n)}=diag\left\{I_{\overline{W}}\otimes {\left(S^{(r)}\right)}^{\oplus n}r=\overline{1R}\right\}n=\overline{1N};$$ E7

$$\mathcal{S}=diag\left\{S^{(r)}r=\overline{1R}\right\}$$ E8

$${\mathcal{S}}_0^{(n)}=diag\left\{I_{\overline{W}}\otimes {\left(S_0^{(r)}\right)}^{\oplus n}r=\overline{1R}\right\}n=\overline{1N}$$ E9

$${\overline{\mathcal{S}}}_0^{(N)}=diag\left\{I_{\overline{W}}\otimes {\left(S_0^{(r)}{\beta }^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}$$ E10

$${\mathcal{C}}^{(n)}=Q\otimes I_{\overline{W}}\otimes I_{M^n}+{\mathcal{D}}_0^{(n)}+{\mathcal{S}}^{(n)}n=\overline{0N}$$ E11

$${\overline{\mathcal{C}}}^{(n)}=Q\otimes I_{\overline{W}}\otimes I_{M^{min\{nN\}}}+{\mathcal{D}}_0^{(min\{nN\})}+{\displaystyle \sum _{k=N+L-n+1}^{\infty }{\mathcal{D}}_k^{(min\{nN\})}}+{\mathcal{S}}^{(min\{nN\})}n=\overline{0N}$$ E12

$${\overline{\mathcal{C}}}^{(N+L)}=Q\otimes I_{\overline{W}}\otimes I_{M^N}+{\displaystyle \sum _{k=0}^{\infty }{\mathcal{D}}_k^{(N)}}+{\mathcal{S}}^{(N)}$$ E13

3. Process of the System States

It is easy to see that operation of the considered queueing model is described in terms of the regular irreducible continuous-time Markov chain

$${\xi }_t=\{n_t{r}_t{\nu }_t{m}_t^{(1)}\dots m_t^{(\min \{n_{t}N\})}\}t\ge 0$$ E14

where

$n_t$ is the number of customers in the system, where $n_t=\overline{0N+L}$ in case of $PA$ and $CR$ disciplines and $n_t\ge 0$ in case of $CA$ discipline; $r_t$ is the state of the random environment, $r_t=\overline{1R}$ ${\nu }_t$ is the state of the $$BMAP$$ underlying process, ${\nu }_t=\overline{0W}$ $m_t^{(n)}$ is the phase of $PH$ service process in the $n$ th busy server, $m_t^{(n)}=\overline{1M}$ $n_t=\overline{1N}$ (we assume here that the busy servers are numerated in order of their occupying, i.e. the server, which begins the service, is appointed the maximal number among all busy servers; when some server finishes the service, the servers are correspondingly enumerated) at epoch $tt\ge 0$

Let us enumerate the states of the chain ${\xi }_{t}t\ge 0$ in the lexicographic order and form the row vectors $p_n$ of probabilities corresponding to the state $n$ of the first component of the process ${\xi }_{t}t\ge 0$ Denote also $p=\left(p_0{p}_1{p}_2\dots \right)$

It is well known that the vector $p$ satisfies the system of the linear algebraic equations (so called equilibrium equations or Chapman-Kolmogorov equations) of the form:

$$pA=0pe=1$$ E15

where $A$ is the infinitesimal generator of the Markov chain ${\xi }_{t}t\ge 0$

Structure of this generator and methods of system (1) solution vary depending on the admission discipline.

3.1. The Case of Partial Admission Discipline

Lemma 1. Infinitesimal generator $A$ of the Markov chain ${\xi }_{t}t\ge 0$ in the case of partial admission discipline has the following block structure:

$$A={\left(A_{n{n}^{\prime }}\right)}_{n{n}^{\prime }=\overline{0N+L}}=$$ E16

$$=\left(\begin{array}{ccccccccc}{\mathcal{C}}^{(0)}& {\Psi }_{\mathrm{1,1}}^{(0)}& \dots & {\Psi }_{N-1N-1}^{(0)}& {\Psi }_{NN}^{(0)}& {\Psi }_{N+1N}^{(0)}& \dots & {\Psi }_{N+L-1N}^{(0)}& {\widehat{\Psi }}_{N+LN}^{(0)}\\ {\mathcal{S}}_0^{(1)}& {\mathcal{C}}^{(1)}& \dots & {\Psi }_{N-\mathrm{2,}N-2}^{(1)}& {\Psi }_{N-1N-1}^{(1)}& {\Psi }_{NN-1}^{(1)}& \dots & {\Psi }_{N+L-\mathrm{2,}N-1}^{(1)}& {\widehat{\Psi }}_{N+L-1N-1}^{(1)}\\ O& {\mathcal{S}}_0^{(2)}& \dots & {\Psi }_{N-\mathrm{3,}N-3}^{(2)}& {\Psi }_{N-\mathrm{2,}N-2}^{(2)}& {\Psi }_{N-1N-2}^{(2)}& \dots & {\Psi }_{N+L-\mathrm{3,}N-2}^{(2)}& {\widehat{\Psi }}_{N+L-\mathrm{2,}N-2}^{(2)}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ O& O& \dots & {\mathcal{C}}^{(N-1)}& {\Psi }_{\mathrm{1,1}}^{(N-1)}& {\Psi }_{\mathrm{2,1}}^{(N-1)}& \dots & {\Psi }_{L\mathrm{,1}}^{(N-1)}& {\widehat{\Psi }}_{L+\mathrm{1,1}}^{(N-1)}\\ O& O& \dots & {\mathcal{S}}_0^{(N)}& {\mathcal{C}}^{(N)}& {\mathcal{D}}_1^{(N)}& \dots & {\mathcal{D}}_{L-1}^{(N)}& {\widehat{\mathcal{D}}}_L^{(N)}\\ O& O& \dots & O& {\overline{\mathcal{S}}}_0^{(N)}& {\mathcal{C}}^{(N)}& \dots & {\mathcal{D}}_{L-2}^{(N)}& {\widehat{\mathcal{D}}}_{L-1}^{(N)}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ O& O& \dots & O& O& O& \dots & {\overline{\mathcal{S}}}_0^{(N)}& {\overline{\mathcal{C}}}^{(N+L)}\end{array}\right)$$ E17

where

$${\Psi }_{m{m}^{\prime }}^{(k)}={\mathcal{D}}_m^{(k)}{ℬ}_{m^{\prime }}^{(k)}k=\overline{0N-1}m\ge 1m^{\prime }=\overline{1N}$$ E18

$${\widehat{\Psi }}_{m_1{m}^{\prime }}^{(k)}={\displaystyle \sum _{m=m_1}^{\infty }{\Psi }_{m{m}^{\prime }}^{(k)}}{\widehat{\mathcal{D}}}_{m_1}^{(k)}={\displaystyle \sum _{m=m_1}^{\infty }{\mathcal{D}}_m^{(k)}}$$ E19

Proof of the Lemma follows from analysis of Markov chain ${\xi }_{t}t\ge 0$ transitions during an infinitesimal interval. Block entries of the generator have the following meaning. The non-diagonal entries of the matrix ${\mathcal{C}}^{(k)}$ define intensity of transition of the components $\{r_t{\nu }_t{m}_t^{(1)}\dots m_t^{(n_t)}\}$ of the Markov chain ${\xi }_{t}t\ge 0$ which do not lead to the change of the number $k$ of busy servers. The diagonal entries of the matrix ${\mathcal{C}}^{(k)}$ are negative and define, up to the sign, intensity of leaving the corresponding states of the Markov chain ${\xi }_{t}t\ge 0$ The entries of the matrix $${\Psi }_{m{m}^{\prime }}^{(k)}={\mathcal{D}}_m^{(k)}{ℬ}_{m^{\prime }}^{(k)}$$ define intensity of transitions of the components $\{r_t{\nu }_t{m}_t^{(1)}\dots m_t^{(n_t)}\}$ of the Markov chain ${\xi }_{t}t\ge 0$ which are accompanied by arrival of $m$ customers and occupying $m^{\prime }$ servers conditional the number of busy servers is $k$ The entries of the matrix ${\mathcal{S}}_0^{(k)}$ define intensity of transitions, which are accompanied by a departure of a customer, conditional the number of busy servers is $k$

To solve system (1) with the matrix $A$ defined by Lemma 1, we use the effective numerically stable procedure developed in [2] that exploits the special structure of the matrix $A$ (it is upper block Hessenberg) and probabilistic meaning of the unknown vector $p$ This procedure is given by the following statement.

Theorem 1. In case of partial admission, the stationary probability vectors

$p_{i}i=\overline{0N+L}$ are computed as follows:

$$p_l=p_0{F}_{l}l=\overline{1N+L}$$ E20

where the matrices $F_l$ are calculated recurrently:

$$F_l=\left({\overline{A}}_{0l}+{\displaystyle \sum _{i=1}^{l-1}{F}_i{\overline{A}}_{il}}\right){\left(-{\overline{A}}_{ll}\right)}^{-1}l=\overline{1N+L-1}$$ E21

$$F_{N+L}=\left(A_{0N+L}+{\displaystyle \sum _{i=1}^{N+L-1}{F}_i{A}_{iN+L}}\right){\left(-A_{N+LN+L}\right)}^{-1}$$ E22

the matrices ${\overline{A}}_{IN+L}$ are calculated from the backward recursions:

$${\overline{A}}_{IN+L}=A_{iN+L}i=\overline{0N+L}$$ E23

$${\overline{A}}_{il}=A_{il}+{\overline{A}}_{il+1}{G}_{l}i=\overline{0l}l=N+L-1N+L-\mathrm{2,}\dots \mathrm{,0}$$ E24

the matrices $G_{i}i=\overline{0N+L-1}$ are calculated from the backward recursion:

$$G_i={\left(-A_{i+1i+1}-{\displaystyle \sum _{l=1}^{N+L-i-1}{A}_{i+1i+1+l}{G}_{i+l}{G}_{i+l-1}\dots G_{i+1}}\right)}^{-1}{A}_{i+1i}$$ E25

$$i=N+L-1N+L-\mathrm{2,}\dots \mathrm{,0}$$ E26

the vector $p_0$ is calculated as the unique solution to the following system of linear algebraic equations:

$$p_0{\overline{A}}_{\mathrm{0,0}}=0p_0\left({\displaystyle \sum _{l=1}^{N+L}{F}_{l}e+e}\right)=1$$ E27

3.2. The Case of Complete Rejection Discipline

Lemma 2. Infinitesimal generator

$A$ of the Markov chain ${\xi }_{t}t\ge 0$ in the case of complete rejection discipline has the following block structure:

$$A={\left(A_{n{n}^{\prime }}\right)}_{n{n}^{\prime }=\overline{0N+L}}=$$ E28

$$=\left(\begin{array}{ccccccccc}{\overline{\mathcal{C}}}^{(0)}& {\Psi }_{\mathrm{1,1}}^{(0)}& \dots & {\Psi }_{N-1N-1}^{(0)}& {\Psi }_{NN}^{(0)}& {\Psi }_{N+1N}^{(0)}& \dots & {\Psi }_{N+L-1N}^{(0)}& {\Psi }_{N+LN}^{(0)}\\ {\mathcal{S}}_0^{(1)}& {\overline{\mathcal{C}}}^{(1)}& \dots & {\Psi }_{N-\mathrm{2,}N-2}^{(1)}& {\Psi }_{N-1N-1}^{(1)}& {\Psi }_{NN-1}^{(1)}& \dots & {\Psi }_{N+L-\mathrm{2,}N-1}^{(1)}& {\Psi }_{N+L-1N-1}^{(1)}\\ O& {\mathcal{S}}_0^{(2)}& \dots & {\Psi }_{N-\mathrm{3,}N-3}^{(2)}& {\Psi }_{N-\mathrm{2,}N-2}^{(2)}& {\Psi }_{N-1N-2}^{(2)}& \dots & {\Psi }_{N+L-\mathrm{3,}N-2}^{(2)}& {\Psi }_{N+L-\mathrm{2,}N-2}^{(2)}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ O& O& \dots & {\overline{\mathcal{C}}}^{(N-1)}& {\Psi }_{\mathrm{1,1}}^{(N-1)}& {\Psi }_{\mathrm{2,1}}^{(N-1)}& \dots & {\Psi }_{L\mathrm{,1}}^{(N-1)}& {\Psi }_{L+\mathrm{1,1}}^{(N-1)}\\ O& O& \dots & {\mathcal{S}}_0^{(N)}& {\overline{\mathcal{C}}}^{(N)}& {\mathcal{D}}_1^{(N)}& \dots & {\mathcal{D}}_{L-1}^{(N)}& {\mathcal{D}}_L^{(N)}\\ O& O& \dots & O& {\overline{\mathcal{S}}}_0^{(N)}& {\overline{\mathcal{C}}}^{(N)}& \dots & {\mathcal{D}}_{L-2}^{(N)}& {\mathcal{D}}_{L-1}^{(N)}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ O& O& \dots & O& O& O& \dots & {\overline{\mathcal{S}}}_0^{(N)}& {\overline{\mathcal{C}}}^{(N+L)}\end{array}\right)$$ E29

The proof of the Lemma is analogous to the proof of the previous Lemma and takes into account the fact that the number of customers in the system does not change when the number of customers in an arriving batch exceeds the number of free servers.

To solve system (1) with the matrix $A$ defined by Lemma 2, we also use the procedure described by Theorem 1.

3.3. The Case of Complete Admission Discipline

Lemma 3. Infinitesimal generator $A$ of the Markov chain ${\xi }_{t}t\ge 0$ in the case of complete admission discipline has the following block structure:

$$A={\left(A_{n{n}^{\prime }}\right)}_{n{n}^{\prime }\ge 0}=$$ E30

$$=\left(\begin{array}{cccccccccc}{\mathcal{C}}^{(0)}& {\Psi }_{\mathrm{1,1}}^{(0)}& \dots & {\Psi }_{NN}^{(0)}& {\Psi }_{N+1N}^{(0)}& \dots & {\Psi }_{N+L-1N}^{(0)}& {\Psi }_{N+LN}^{(0)}& {\Psi }_{N+L+1N}^{(0)}& \dots \\ {\mathcal{S}}_0^{(1)}& {\mathcal{C}}^{(1)}& \dots & {\Psi }_{N-1N-1}^{(1)}& {\Psi }_{NN-1}^{(1)}& \dots & {\Psi }_{N+L-\mathrm{2,}N-1}^{(1)}& {\Psi }_{N+L-1N-1}^{(1)}& {\Psi }_{N+LN-1}^{(1)}& \dots \\ O& {\mathcal{S}}_0^{(2)}& \dots & {\Psi }_{N-\mathrm{2,}N-2}^{(2)}& {\Psi }_{N-1N-2}^{(2)}& \dots & {\Psi }_{N+L-\mathrm{3,}N-2}^{(2)}& {\Psi }_{N+L-\mathrm{2,}N-2}^{(2)}& {\Psi }_{N+L-1N-2}^{(2)}& \dots \\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮\\ O& O& \dots & {\Psi }_{\mathrm{1,1}}^{(N-1)}& {\Psi }_{\mathrm{2,1}}^{(N-1)}& \dots & {\Psi }_{L\mathrm{,1}}^{(N-1)}& {\Psi }_{L+\mathrm{1,1}}^{(N-1)}& {\Psi }_{L+\mathrm{2,1}}^{(N-1)}& \dots \\ O& O& \dots & {\mathcal{C}}^{(N)}& {\mathcal{D}}_1^{(N)}& \dots & {\mathcal{D}}_{L-1}^{(N)}& {\mathcal{D}}_L^{(N)}& {\mathcal{D}}_{L+1}^{(N)}& \dots \\ O& O& \dots & {\overline{\mathcal{S}}}_0^{(N)}& {\mathcal{C}}^{(N)}& \dots & {\mathcal{D}}_{L-2}^{(N)}& {\mathcal{D}}_{L-1}^{(N)}& {\mathcal{D}}_L^{(N)}& \dots \\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \dots \\ O& O& \dots & O& O& \dots & {\mathcal{C}}^{(N)}& {\mathcal{D}}_1^{(N)}& {\mathcal{D}}_2^{(N)}& \dots \\ O& O& \dots & O& O& \dots & {\overline{\mathcal{S}}}_0^{(N)}& {\overline{\mathcal{C}}}^{(N+L)}& O& \dots \\ O& O& \dots & O& O& \dots & O& {\overline{\mathcal{S}}}_0^{(N)}& {\overline{\mathcal{C}}}^{(N+L)}& \dots \\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots \end{array}\right)$$ E31

Essential difference of complete admission discipline is that the state space of the Markov chain ${\xi }_{t}t\ge 0$ is infinite and this makes its analysis more complicated. However, the block rows, except the first $N+L$ boundary block rows, have only two non-zero blocks and this Markov chain behaves as Quasi-Death process when the state of the first component $n_t$ of the Markov chain ${\xi }_{t}t\ge 0$ if greater than $N+L$ It allows to construct effective stable algorithm for calculation of the stationary distribution of this Markov chain. Note, that although the state space of the Markov chain ${\xi }_{t}t\ge 0$ is infinite, this Markov chain is ergodic under the standard assumptions about the parameters of the $$BMAP$$ input, the $PH$ type service and the random environment. The algorithm for calculation of the stationary distribution is given in the following statement.

Theorem 2. In case of complete admission discipline, the stationary probability vectors

$p_{l}l\ge 0$ are calculated as follows:

$$p_l=p_0{F}_{l}l\ge 1$$ E32

where the matrices $F_l$ are calculated recurrently

$$F_l=\left({\overline{A}}_{0l}+{\displaystyle \sum _{i=1}^{l-1}{F}_i{\overline{A}}_{il}}\right){\left(-{\overline{A}}_{ll}\right)}^{-1}l\ge 1$$ E33

the matrices ${\overline{A}}_{iI}$ are calculated as:

$${\overline{A}}_{il}=A_{il}+{\displaystyle \sum _{k=1}^{\infty }{A}_{il+k}}{G}^{\max \{0l+k-N-L\}}{G}_{\min \{N+Ll+k\}-1}{G}_{\min \{N+Ll+k\}-2}\dots G_{l}i=0\dots ll\ge 1$$ E34

the matrix $G$ has a form

$$G=-{\left({\overline{\mathcal{C}}}^{(N+L)}\right)}^{-1}{\overline{\mathcal{S}}}_0^{(N)}$$ E35

the matrices $G_{i}i=\overline{0N+L-1}$ are calculated from the backward recursion

$$G_i=-{\left(A_{i+1i+1}+{\displaystyle \sum _{l=i+2}^{\infty }{A}_{i+1l}{G}^{\max \{0l-N-L\}}}{G}_{\min \{N+Ll\}-1}{G}_{\min \{N+Ll\}-2}{G}_{i+1}\right)}^{-1}{A}_{i+1i}$$ E36

the vector $p_0$ is the unique solution of the system:

$$p_0{\overline{A}}_{\mathrm{0,0}}=0p_0\left({\displaystyle \sum _{l=1}^{\infty }{F}_{l}e+e}\right)=1$$ E37

The proof of the Theorem follows from the theory of multi-dimensional Markov chains with continuous time, see [6]. It is worth to note that Neuts' matrix $G$ which is usually found numerically as solution to matrix equation, see [9], here is obtained in the explicit form.

3.4. The Case of an Infinite Size of a Buffer $$\left(L=\infty \right)$$

The system under consideration in this section has an infinite waiting space. If an arriving batch of customers sees idle servers, a part of the batch corresponding to the number of free servers occupy these servers while the rest of the batch joins the queue. If the system has all servers being busy at a batch arrival epoch, all customer of the batch go to the queue.

Lemma 3. Infinitesimal generator

$A$ of the Markov chain ${\xi }_{t}t\ge 0$ has the following block structure:

$$A={\left(A_{n{n}^{\prime }}\right)}_{n{n}^{\prime }\ge 0}=$$ E38

$$=\left(\begin{array}{ccccccccc}{\mathcal{C}}^{(0)}& {\Psi }_{\mathrm{1,1}}^{(0)}& {\Psi }_{\mathrm{2,2}}^{(0)}& \dots & {\Psi }_{N-1N-1}^{(0)}& {\Psi }_{NN}^{(0)}& {\Psi }_{N+1N}^{(0)}& {\Psi }_{N+\mathrm{2,}N}^{(0)}& \dots \\ {\mathcal{S}}_0^{(1)}& {\mathcal{C}}^{(1)}& {\Psi }_{\mathrm{1,1}}^{(1)}& \dots & {\Psi }_{N-\mathrm{2,}N-2}^{(1)}& {\Psi }_{N-1N-1}^{(1)}& {\Psi }_{NN-1}^{(1)}& {\Psi }_{N+1N-1}^{(1)}& \dots \\ O& {\mathcal{S}}_0^{(2)}& {\mathcal{C}}^{(2)}& \dots & {\Psi }_{N-\mathrm{3,}N-3}^{(2)}& {\Psi }_{N-\mathrm{2,}N-2}^{(2)}& {\Psi }_{N-1N-2}^{(2)}& {\Psi }_{NN-2}^{(2)}& \dots \\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& \ddots \\ O& O& O& \dots & {\mathcal{C}}^{(N-1)}& {\Psi }_{\mathrm{1,1}}^{(N-1)}& {\Psi }_{\mathrm{2,1}}^{(N-1)}& {\Psi }_{\mathrm{3,1}}^{(N-1)}& \dots \\ O& O& O& \dots & {\mathcal{S}}_0^{(N)}& {\mathcal{C}}^{(N)}& {\mathcal{D}}_1^{(N)}& {\mathcal{D}}_2^{(N)}& \dots \\ O& O& O& \dots & O& {\overline{\mathcal{S}}}_0^{(N)}& {\mathcal{C}}^{(N)}& {\mathcal{D}}_1^{(N)}& \dots \\ O& O& O& \dots & O& O& {\overline{\mathcal{S}}}_0^{(N)}& {\mathcal{C}}^{(N)}& \dots \\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right)$$ E39

In what follows we perform the steady state analysis of the Markov chain having generator of form (2). To this end, we use the results for continuous time multi-dimensional Markov chain ( $$QTMC$$ ) presented in [6].

Theorem 3. The necessary and sufficient condition for existence of the Markov chain

${\xi }_{t}t\ge 0$ stationary distribution is the fulfillment of the inequality

$$\rho =\lambda /\overline{\mu }1$$ E40

where

$$\lambda =x_1{D}^{\text{'}}{(z)}^{}|{}_{z=1}^{}e$$ E41

$$\overline{\mu }=x_{2}diag\left\{{\left(S_0^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}e$$ E42

the vectors $x_{n}n=\mathrm{1,2}$ are the unique solutions to the following systems of linear algebraic equations:

$$x_1\left(Q\otimes I_{\overline{W}}+D(1)\right)=0x_{1}e=1$$ E43

$$x_2\left(Q\otimes I_{M^N}+diag\left\{{\left(S^{(r)}+S_0^{(r)}{\beta }^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}\right)=0x_{2}e=1$$ E44

Proof. Using the results of [6], we directly obtain the desired condition in the form of inequality

$$x{\left[{\mathcal{C}}^{(N)}z+z{\mathcal{D}}^{(N)}(z)\right]}^{\prime }{}_{z=1}e0$$ E45

where $x$ is the unique solution to the system

$$x\left(Q\otimes I_{\overline{W}{M}^N}+{\mathcal{D}}^{(N)}(1)+{\mathcal{S}}^{(N)}+{\overline{\mathcal{S}}}_0^{(N)}\right)=0xe=1$$ E46

It is easy to show that inequality (7) is reduced to the following inequality:

$$x_1{D}^{\text{'}}(z)|{}_{z=1}^{}e{x}_{2}diag\left\{{\left(S_0^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}e$$ E47

where $x_1=x\left(I_{R\overline{W}}\otimes e_{M^N}\right)$

$$x_2=x\left(I_R\otimes e_{\overline{W}}\otimes I_{M^N}\right)$$ E48

To get the equations for the row vectors $x_1$ and $x_2$ we multiply equation (8) by the matrices $I_{R\overline{W}}\otimes e_{M^N}$ and $I_R\otimes e_{\overline{W}}\otimes I_{M^N}$ respectively. After multiplication and some algebra we obtain equations (5), (6) for the vectors $x_1$ and $x_2$ So, inequality (9) is equivalent to inequality (3) and the theorem is proved.

The value $\rho$ has a meaning of the system load. In what follows we assume inequality (3) be fulfilled.

To solve system (1) with the matrix $A$ defined by (2), we use the effective numerically stable procedure [6] based on the account special structure of the matrix $A$ notion of the censored Markov chain and probabilistic meaning of the unknown vector $p$ For more detail see [6]. This procedure is given by the following statement.

Theorem 4. The stationary probability vectors

$p_{l}l\ge 0$ are calculated as follows:

$$p_l=p_0{F}_{l}l\ge 1$$ E49

where the matrices $F_l$ are calculated recurrently:

$$F_l=\left({\overline{A}}_{0l}+{\displaystyle \sum _{i=1}^{l-1}{F}_i{\overline{A}}_{il}}\right){\left(-{\overline{A}}_{ll}\right)}^{-1}l\ge 1$$ E50

the matrices ${\overline{A}}_{iI}$ are calculated as:

$${\overline{A}}_{il}=A_{il}+{\overline{A}}_{il+1}{G}_{l}0\le i\le lN$$ E51

$${\overline{A}}_{il}=A_{il}+{\overline{A}}_{il+1}Gl\ge \max \left\{iN\right\}i\ge 0$$ E52

the matrix $G$ is calculated from the equation

$$G={\left(-{\displaystyle \sum _{l=1}^{\infty }{A}_{N+1N+l}}{G}^{l-1}\right)}^{-1}{A}_{N+1N}$$ E53

the matrices $G_{i}i=\overline{0N-1}$ are calculated from the backward recursion:

$$G_i=-{\left(A_{i+1i+1}+{\displaystyle \sum _{l=i+2}^{\infty }{A}_{i+1l}{G}^{\max \{0l-N\}}}{G}_{\min \{Nl\}-1}{G}_{\min \{Nl\}-2}\dots G_{i+1}\right)}^{-1}{A}_{i+1i}$$ E54

the vector $p_0$ is calculated as the unique solution to the following system of linear algebraic equations:

$$p_0{\overline{A}}_{\mathrm{0,0}}=0p_0\left({\displaystyle \sum _{l=1}^{\infty }{F}_{l}e+e}\right)=1$$ E55

4. Performance Measures

Having the probability vector $p$ been computed, we are able to calculate performance measures of the considered model. The main performance measure in the case of a finite buffer is the probability $P_{loss}$ that an arbitrary customer will be lost (the loss probability).

Theorem 5. The loss probability

$P_{loss}$ is calculated as follows

(i) in the case of $PA$ discipline

$$P_{loss}=1-\frac{1}{\lambda }{\displaystyle \sum _{i=0}^{N+L-1}{p}_i}{\displaystyle \sum _{k=0}^{N+L-i}(k+i-N-L){\mathcal{D}}_k^{(i)}e}$$ E56

(ii) in the case of $CR$ discipline

$$P_{loss}=1-\frac{1}{\lambda }{\displaystyle \sum _{i=0}^{N+L-1}{p}_i}{\displaystyle \sum _{k=0}^{N+L-i}k{\mathcal{D}}_k^{(i)}e}$$ E57

(iii) in the case of $CA$ discipline

$$P_{loss}=1-\frac{1}{\lambda }{\displaystyle \sum _{i=0}^{N+L-1}{p}_i}{\displaystyle \sum _{k=0}^{\infty }k{\mathcal{D}}_k^{(i)}e}$$ E58

Proofs of formulae (10) - (12) are analogous. So, we will prove only formula (10). According to a formula of the total probability, the probability $P_{loss}$ is calculated as

$$P_{loss}=1-{\displaystyle \sum _{i=0}^{N+L-1}{\displaystyle \sum _{k=1}^{\infty }{P}_k{P}_i^{(k)}{R}^{(ik)}}}$$ E59

where $P_k$ is a probability that an arbitrary customer arrives in a batch consisting of $k$ customers; $P_i^{(k)}$ is a probability to see $i$ servers being busy at the epoch of the $k-$ size batch arrival; $R^{(ik)}$ is a probability that an arbitrary customer will not be lost conditional it arrives in a batch consisting of $k$ customers and $i$ servers are busy at the arrival epoch.

It can be shown that

$$P_i^{(k)}=\frac{p_i{\mathcal{D}}_k^{(i)}e}{x_1{\mathcal{D}}_k^{(0)}e}i=\overline{0N+L-1}k\ge 1$$ E60

$$P_k=\frac{k{x}_1{\mathcal{D}}_k^{(0)}e}{x_{1}D(z)|{}_{z=1}^{\text{'}}e}=k\frac{x_1{\mathcal{D}}_k^{(0)}e}{\lambda }k\ge 1$$ E61

$$R^{(ik)}=\{\begin{array}{l}1k\le N+L-i\\ \frac{N+L-i}{k}kN+L-ii=\overline{0N+L-1}\end{array}$$ E62

By substituting (14)-(16) into (13) after some algebra we get (10). $$$$

Some performance measures for the case $L=\infty$ are presented below.

The probability to see $i$ customers in the system

$$p_i=p_{i}ei\ge 0$$ E63

The mean number of customers in the system

$$L_{queue}={\displaystyle \sum _{i=0}^{\infty }i}{p}_{i}e$$ E64

The probability to see $n$ busy servers

$$p_n=p_{n}en=\overline{0N-1}{p}_N={\displaystyle \sum _{n=N}^{\infty }{p}_{n}e}$$ E65

The mean number of busy servers

$$N_{busy}={\displaystyle \sum _{n=1}^{N}n{p}_{n}e+N{\displaystyle \sum _{n=N}^{\infty }{p}_{n}e}}$$ E66

The mean number $N_{idle}$ of idle servers

$$N_{idle}=N-N_{busy}$$ E67

The vector ${\widehat{p}}_n$ whose $\left(\overline{W}(r-1)+\nu +1\right)$ -th entry is the joint probability to see $n$ busy servers, the random environment in the state $r$ and the process ${\nu }_t$ in the state $\nu$

$${\widehat{p}}_n=p_n\left(I_{R\overline{W}}\otimes e_{M^n}\right)n=\overline{0N-1}{\widehat{p}}_N={\displaystyle \sum _{n=N}^{\infty }{p}_n\left(I_{R\overline{W}}\otimes e_{M^n}\right)}$$ E68

The vector of conditional means of the number of busy servers under the fixed states of the random environment

$$n={\displaystyle \sum _{n=1}^{\infty }\min \left\{nN\right\}}{\widehat{p}}_n{}_n\left(I_R\otimes e_{\overline{W}}\right)diag\left\{q_r^{-1}r=\overline{1R}\right\}$$ E69

The vector $p^{(a)}(n)$ whose $\left(\overline{W}(r-1)+\nu +1\right)$ -th entry is the joint probability that an arbitrary arriving call sees $n$ busy servers and the random environment in the state $r$ and the state of the process ${\nu }_t$ becomes $\nu$ after the arrival epoch

$$p^{(a)}(n)={\lambda }^{-1}{\widehat{p}}_n{D}^{\text{'}}(z)|{}_{z=1}^{}n=\overline{0N}$$ E70

The probability $p^{(a)}(n)$ that an arbitrary arriving call sees $n$ busy servers

$$p^{(a)}(n)=p^{(a)}(n)en=\overline{0N}$$ E71

The vector $p_b^{(a)}(n)$ whose $\left(\overline{W}(r-1)+\nu +1\right)$ -th entry is the joint probability that an arbitrary arriving batch of size $k$ sees $n$ busy servers and the random environment in the state $r$ and the state of the process ${\nu }_t$ becomes $\nu$ after the arrival epoch

$$p_b^{(a)}(kn)={\lambda }_b^{-1}{\widehat{p}}_{n}diag\left\{D_k^{(r)}r=\overline{1R}\right\}n=\overline{0N}k\ge 1$$ E72

where

$${\lambda }_b=x\left(D(1)-D(0)\right)e$$ E73

The probability $p_b^{(a)}(kn)$ that an arbitrary arriving batch of size $k$ sees $n$ busy servers

$$p_b^{(a)}(kn)=p_b^{(a)}(kn)en=\overline{0N}k\ge 1$$ E74

The vector $p_b^{(a)}(n)$ whose $\left(\overline{W}(r-1)+\nu +1\right)$ -th entry is the joint probability that an arbitrary arriving batch sees $n$ busy servers and the random environment in the state $r$ and the state of the process ${\nu }_t$ becomes $\nu$ after the arrival epoch

$$p_b^{(a)}(n)={\lambda }_b^{-1}{\widehat{p}}_n\left(D(1)-D(0)\right)n=\overline{0N}$$ E75

The probability $p_b^{(a)}(n)$ that an arbitrary arriving batch sees $n$ busy servers

$$p_b^{(a)}(n)=p_b^{(a)}(n)en=\overline{0N}$$ E76

The probability $P_{imm}$ that an arbitrary customer will enter the service immediately upon arrival (without visiting a buffer)

$$P_{imm}={\lambda }^{-1}{\displaystyle \sum _{i=0}^{N-1}{p}_i}{\displaystyle \sum _{k=0}^{N-i}(k+i-N){\mathcal{D}}_k^{(i)}e}$$ E77

5. Actual Sojourn Time

Let ${\nu }_a(s)Res0$ be the Laplace-Stieltjes transform ( $LST$ ) of the sojourn time distribution and ${\overline{\nu }}_a$ be the mean sojourn time of the arbitrary customer in the system.

Theorem 6. The Laplace-Stieltjes transform

${\nu }_a(s)$ is calculated as follows

$${\nu }_a(s)=\frac{1}{\lambda }\{{\displaystyle \sum _{i=0}^{N-1}{p}_i}{\displaystyle \sum _{k=1}^{\infty }\min \left\{kN-i\right\}{\mathcal{D}}_k^{(i)}\left(I_R\otimes e_{\overline{W}{m}^i}\right)}+$$ E78

$$+{\displaystyle \sum _{i=0}^{\infty }{p}_i}{\displaystyle \sum _{k=\max \{1N-i+1\}}^{\infty }{\mathcal{D}}_k^{(\min \{iN\})}{ℬ}^{(\max \{0N-i\})}{\displaystyle \sum _{l=\max \{1N-i+1\}}^k{\left(ℱ(s)\right)}^{i-N+l}\left(I_R\otimes e_{M^N}\right)}}\}\mathscr{H}(s)e_R$$ E79

where

$$\mathscr{H}(s)=-diag\left\{{\beta }^{(r)}r=\overline{1R}\right\}{\left(sI-(Q\otimes I_M+\mathcal{S})\right)}^{-1}diag\left\{S_0^{(r)}r=\overline{1R}\right\}$$ E80

$$ℱ(s)=-{\left(sI-\left(Q\otimes I_{M^N}+diag\left\{{\left(S^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}\right)\right)}^{-1}diag\left\{{\left(S_0^{(r)}{\beta }^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}$$ E81

$${ℬ}^{(n)}=diag\left\{e_{\overline{W}}\otimes I_{M^{N-n}}\otimes {\left({\beta }^{(r)}\right)}^{\oplus n}r=\overline{1R}\right\}n=\overline{0N}$$ E82

$$\mathcal{S}=diag\left\{S^{(r)}r=\overline{1R}\right\}$$ E83

Proof. We derive the expression for the

$LST$ ${\nu }_a(s)$ by means of the method of collective marks (method of additional event, method of catastrophes) for references see, e.g. [3], [11]. To this end, we interpret the variable $s$ as the intensity of some virtual stationary Poisson flow of catastrophes. So, ${\nu }_a(s)$ has the meaning of probability that no one catastrophe arrives during the sojourn time of an arbitrary customer. Then, the proof of the theorem follows from the formula of total probability if we analyze the states of the system at an arbitrary customer arrival epoch and take into account the probabilistic meaning of the involved matrices. The matrix $\mathscr{H}(s)$ is the matrix $LST$ of an arbitrary customer service time distribution. It is the $$R$$ -size square matrix whose $$(r{r}^{\prime })$$ entry is a probability that during the service time of a customer a catastrophe does not arrive and the process $r_{t}t\ge 0$ transits from the state $r$ to the state $r^{\prime }r{r}^{\prime }=\overline{1R}$ It is defined by the formula:

$$\mathscr{H}(s)=diag\left\{{\beta }^{(r)}r=\overline{1R}\right\}{\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-st}}{e}^{(Q\otimes I_M+diag\{S^{(r)}r=\overline{1R}\})t}dtdiag\left\{S_0^{(r)}r=\overline{1R}\right\}$$ E84

Analogously, the entries of the matrix $LST$ $ℱ(s)$ are the probabilities of no catastrophe arrival and corresponding transitions of the process $\{r_t{m}_t^{(1)}\dots m_t^{(N)}\}t\ge 0$ during the time interval from an arbitrary moment when all $N$ servers are busy till the first epoch when one of these servers finishes the service of a customer. This matrix is defined by the formula:

$$ℱ(s)={\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-st}}{e}^{(Q\otimes I_{M^N}+diag\{{(S^{(r)})}^{\oplus N}r=\overline{1R}\})t}dtdiag\left\{{\left(S_0^{(r)}{\beta }^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}$$ E85

Theorem 7. The mean sojourn time

${\overline{\nu }}_a$ of an arbitrary customer in the system is calculated by

$${\overline{\nu }}_a=-\frac{1}{\lambda }\{{\displaystyle \sum _{i=0}^{\infty }{p}_i}{\displaystyle \sum _{k=\max \{1N-i+1\}}^{\infty }{\mathcal{D}}_k^{(\min \{iN\})}{ℬ}^{(\max \{0N-i\})}}{\displaystyle \sum _{l=\max \{1N-i+1\}}^k{\displaystyle \sum _{m=0}^{i+l-N-1}{\left(ℱ(0)\right)}^m{ℱ}^{\prime }(0)\left(I_R\otimes e_{M^N}\right)+}}$$ E86

$$+{\displaystyle \sum _{i=0}^{N-1}{p}_i}{\displaystyle \sum _{k=1}^{\infty }\min \{kN-i\}{\mathcal{D}}_k^{(i)}\left(I_R\otimes e_{\overline{W}{M}^i}\right){\mathscr{H}}^{\prime }(0)}+$$ E87

$$+{\displaystyle \sum _{i=0}^{\infty }{p}_i}{\displaystyle \sum _{k=\max \{1N-i+1\}}^{\infty }{\mathcal{D}}_k^{(\min \{iN\})}{ℬ}^{(\max \{0N-i\})}}{\displaystyle \sum _{l=\max \{1N-i+1\}}^k{\left(ℱ(0)\right)}^{i-N+l}\left(I_R\otimes e_{M^{\mathcal{N}}}\right){\mathscr{H}}^{\prime }(0)}\}{e}_R$$ E88

where

$${\mathscr{H}}^{\prime }(0)=-diag\left\{{\beta }^{(r)}r=\overline{1R}\right\}{\left[Q\otimes I_M+\mathcal{S}\right]}^{-2}diag\left\{S_0^{(r)}r=\overline{1R}\right\}$$ E89

$${ℱ}^{\prime }(0)=-{\left[Q\otimes I_{M^N}+diag\left\{{\left(S^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}\right]}^{-2}diag\left\{{\left(S_0^{(r)}{\beta }^{(r)}\right)}^{\oplus N}r=\overline{1R}\right\}$$ E90

Proof. To get expression (18) for ${\overline{\nu }}_a$ we differentiate (17) at the point $s=0$ and use the formula ${\overline{\nu }}_a=-{\overline{\nu }}_a^{\text{'}}(0)$

6. Numerical Examples

The goal of the numerical experiments is to demonstrate the feasibility of the proposed algorithms for computing the stationary distributions of the number of customers and the sojourn time in the system and to give some insight into behavior of the considered queueing systems. In particular, the following issues are addressed:

Comparison of the mean sojourn time of an arbitrary customer and the probability of immediate access to the servers in the systems with varying traffic intensities and different coefficients of correlation in the $BMAP$ s (experiment $\#1$ );

Comparison of the mean sojourn time of an arbitrary customer and the probability of immediate access to the servers in the original system in a $RE$ and in more simple queueing systems for different system loads (experiments $\#2$ );

Demonstration of possible positive effect of redistribution of traffic between the peak traffic periods and normal traffic periods (experiment $\#3$ );

Comparison of the exact value of performance measures of the system in a $RE$ and their simple engineering approximations in cases of slowly and quickly varying $RE$ (experiment $\#4$ );

Investigation of the rate of convergence of the mean sojourn time and the probability of immediate access in the system with the finite buffer to corresponding performance measures of the system with an infinite buffer when the buffer size increases (experiment $\#5$ );

Demonstration of the possibility to apply the presented results for optimization of the number of servers in the system (experiment $\#6$ ).

In numerical examples, we consider the systems operating in the $RE$ which has two states ( $R=2$ ). The generator of the random environment is $Q=\left(\begin{array}{cc}-5& 5\\ 15& -15\end{array}\right)$ The stationary distribution of the $RE$ states is defined by the vector $q=\left(\mathrm{0.75,}0.25\right)$ The number of servers is $N=3$

In the presented examples, we will use several different $MAP$ s and $BMAP$ s for description of the arrival process and two $PH$ type distributions for description of the service processes under the fixed value of the $RE$ For the use in the sequel, let us define these processes.

We consider four arrival processes $MA{P}_{r}r=\overline{\mathrm{1,4}}$ $MA{P}_r$ is defined by the matrices $D_0^{(r)}$

$$D_1^{(r)}r=\overline{\mathrm{1,4}}$$ E91

where

$$D_0^{(1)}=\left(\begin{array}{ccc}-3.9& 0.15& 0.15\\ 0.13& -0.6& 0.1\\ 0.15& 0.14& -0.5\end{array}\right)D_1^{(1)}=\left(\begin{array}{ccc}3.5& 0.08& 0.02\\ 0.03& 0.3& 0.04\\ 0.02& 0.06& 0.13\end{array}\right)$$ E92

$$D_0^{(2)}=\left(\begin{array}{ccc}-6.4& 0.1& 0.1\\ 0.04& -0.6& 0.1\\ 0.07& 0.1& -0.44\end{array}\right)D_1^{(2)}=\left(\begin{array}{ccc}6.06& 0.12& 0.02\\ 0.01& 0.4& 0.05\\ 0.01& 0.06& 0.2\end{array}\right)$$ E93

$$D_0^{(3)}=\left(\begin{array}{ccc}-2.9& 0.73& 0.77\\ 0.87& -3.06& 0.53\\ 0.85& 0.5& -2\end{array}\right)D_1^{(3)}=\left(\begin{array}{ccc}0.68& 0.45& 0.27\\ 0.48& 1.08& 0.1\\ 0.35& 0.05& 0.25\end{array}\right)$$ E94

$$D_0^{(4)}=\left(\begin{array}{ccc}-1.3& 0.21& 0.17\\ 0.16& -2.04& 0.21\\ 0.01& 0.16& -1.3\end{array}\right)D_1^{(4)}=\left(\begin{array}{ccc}0.46& 0.32& 0.14\\ 0.13& 1.34& 0.2\\ 0.02& 0.01& 1.1\end{array}\right)$$ E95

All these $MAP$ s have fundamental rate ${\lambda }^{(r)}=1.25$ The $MA{P}_1$ has the squared variation coefficient ${\left(c_{var}^{(1)}\right)}^2=4$ and the coefficient of correlation of the lengths of successive inter-arrival times $c_{cor}^{(1)}=0.2$ For the rest of the $MAP$ s, the corresponding parameters are: ${\left(c_{var}^{(2)}\right)}^2=4$ $c_{cor}^{(2)}=0.3$ ${\left(c_{var}^{(3)}\right)}^2=1.097$ $c_{cor}^{(3)}=0.0052$ ${\left(c_{var}^{(4)}\right)}^2=1.037$

$$c_{cor}^{(4)}=0.0065$$ E96

Based on these $MAP$ s, we construct batch flows $BMAP$ s as follows. If the $MAP$ is defined by the matrices $D_0^{(r)}$ and $D_1^{(r)}r=\overline{\mathrm{1,4}}$ then the $BMAP$ having the maximal size of a batch equal to $K$ is defined by the matrices $D_0^{(r)}$ $D_k^{(r)}=D_1^{(r)}{q}^{k-1}(1-q)/(1-q^K)k=\overline{1K}r=\overline{1R}$ where

$$q=0.9$$ E97

Following this way, we construct the $BMA{P}_1$ $BMA{P}_2$ $BMA{P}_3$ $BMA{P}_4$ flows based on the $MA{P}_1$ $MA{P}_2$ $MA{P}_3$ $MA{P}_4$ correspondingly, with $K=5$ Note that the coefficients of variation and correlation of all $BMAP$ s are the same as these coefficients for the corresponding $MAP$ s. Fundamental rate ${\lambda }^{(r)}$ and the mean batch size ${\overline{k}}^{(r)}$ of the $BMAP$ s are the following: ${\lambda }^{(1)}={\lambda }^{(2)}={\lambda }^{(3)}={\lambda }^{(4)}=3.488$

$${\overline{k}}^{(1)}={\overline{k}}^{(2)}={\overline{k}}^{(3)}={\overline{k}}^{(4)}=1.989$$ E98

The $P{H}_{r}r=\overline{\mathrm{1,2}}$ service processes are defined by the vectors ${\beta }^{(1)}=\left(10\right)$ ${\beta }^{(2)}=\left(\mathrm{0.2,}0.8\right)$ and the matrices

$$S^{(1)}=\left(\begin{array}{cc}-4& 4\\ 0& -4\end{array}\right)S^{(2)}=\left(\begin{array}{cc}-10& 2\\ 2& -20\end{array}\right)$$ E99

The mean rates of service are ${\mu }^{(1)}=\mathrm{2,}{\mu }^{(2)}=14$ The coefficients of variation of the service time distribution are defined by ${\left(c_{var}^{(1)}\right)}^2=0.5$

$${\left(c_{var}^{(2)}\right)}^2=1.24$$ E100

In the first experiment, we compare the dependence of ${\overline{\nu }}_a$ and $P_{imm}$ on the system load $\rho$ for the $BMAP$ s with different correlations.

In the experiment we use service processes defined by $P{H}_1$ and $P{H}_2$ and four different input flows which are described by $BMA{P}_1$ $BMA{P}_2$ $BMA{P}_3$ and $BMA{P}_4$ having the same mean fundamental rate equal to 3.488 but different correlation coefficients.

We consider three queueing systems which have different combinations of the $BMAP$ s under the first and second states of the $RE$

The input flow in the first system is defined by $BMA{P}_1$ and $BMA{P}_2$ These $BMAP$ s have large coefficients of correlation $c_{cor}^{(1)}=0.2$ and

$$c_{cor}^{(2)}=0.3$$ E101

The input flow in the second system is defined by $BMA{P}_3$ and $BMA{P}_4$ These $BMAP$ s have small coefficients of correlation $c_{cor}^{(3)}=0.0052$ and

$$c_{cor}^{(4)}=0.0065$$ E102

In the third system the input is defined by $BMA{P}_1$ and $BMA{P}_4$ The correlation coefficients of these $BMAP$ s differ significantly.

Figures 1 and 2 show the dependence of the mean sojourn time ${\overline{\nu }}_a$ and the probability $P_{imm}$ on the system load $\rho$ for all these systems. Variation of the value of $\rho$ in all experiments is performed by means of multiplying the entries of the matrices, which define the corresponding $BMAP$ by some varying factor $\gamma$ This implies the increase of the fundamental rate of all the $BMAP$ by a factor $\gamma$ Service time distributions are not modified. It is clear from Figure 1 that correlation in $BMAP$ has a great impact on the sojourn time in the system. An increase of correlation at least in one of the $BMAP$ s describing input in the system implies an increase of the sojourn time in the system in all range of the system load.

Figure 1.

Mean sojourn time in the system as a function of the system load for different correlations in the $BMAP$ s

Figure 2.

Probability of immediate access to the servers in the system as a function of the system load for different correlations in the $BMAP$ s

In the second experiment we compare the values ${\overline{\nu }}_a$ and $P_{imm}$ in the $BMAP/PH/N$ system operating in the $RE$ and in more simple queueing systems which can be considered as its simplified analogs. The first type analog is the $M^X/PH/N$ system in the $RE$ where, under the fixed value of the $RE$ the input flow is a group stationary Poisson with the same batch size distribution and intensity equal to fundamental rate of the corresponding $BMAP$ in the original system. The second type analog is the system $M^X/M/N$ with parameters of arrival and service processes which are obtained by means of averaging, according to stationary distribution of the $RE$ parameters of the original system.

Input flow is described by $BMA{P}_1$ and $BMA{P}_2$ Service processes are $P{H}_1$ and

$$P{H}_2$$ E103

Figures 3 and 4 show the dependence of the the mean sojourn time ${\overline{\nu }}_a$ and the probability $P_{imm}$ on the value of $\rho$

Figure 3.

Mean sojourn time in original system and more simple queueing systems

Figure 4.

Probability of immediate access to the servers in original system and more simple queueing systems

It can be seen from Figures 3 and 4 that an approximation of the mean sojourn time and the probability that an arbitrary call reaches the server immediately by means of their values in some specially constructed more simple queueing system can be rather bad.

The idea of the third experiment is the following. Let us assume that the $RE$ has two states. One state corresponds to the peak traffic periods, the second one corresponds to the normal traffic periods. Service times during these periods are defined by $P{H}_1$ and $P{H}_2$ distributions. Arrivals during these periods are defined by the stationary Poisson flow with the rates ${\lambda }_1$ and ${\lambda }_2$ correspondingly and initially we assume that ${\lambda }_1\gg {\lambda }_2$ It is intuitively clear that if it is possible to redistribute the arrival processes (i.e., to reduce the arrival rate during the peak periods and to increase it correspondingly during the normal traffic periods) without changing the total average arrival rate, the mean sojourn time in the system can be reduced. In real life system such a redistribution is sometimes possible, e.g., by means of controlling tariffs during the peak traffic periods. The goal of this experiment is to show that this intuitive consideration is correct and to illustrate the effect of the redistribution.

We assume that the averaged arrival rate $\lambda$ should be 12.5 and consider four different situations: a huge difference of arrival rates ${\lambda }_1=50{\lambda }_2$ a very big difference ${\lambda }_1=10{\lambda }_2$ a big difference ${\lambda }_1=3{\lambda }_2$ and equal arrival rates ${\lambda }_1={\lambda }_2$ The generator of the random environment is $Q=\left(\begin{array}{cc}-15& 15\\ 5& -5\end{array}\right)$

It can be seen from Figures 5 and 6 that the smoothing of the peak rates can cause essential decrease of the mean sojourn time and the increase of the probability that an arbitrary call reaches the server immediately upon arrival in the system.

Figure 5.

Mean sojourn time as a function of system load for different relations of arrival rates

Figure 6.

Probability of immediate access to the servers as a function of system load for different relations of arrival rates

In the second experiment, we have seen that an approximation of the system performance measures by means of their values in more simple queueing system can be bad. However, it is intuitively clear the following. If the random environment is "very slow" (the rate of the $RE$ is much less then the rates of the input flow and the service processes), an approximation called below as "mixed system" can be applied successfully. This approximation consists of calculation of the system characteristics under the fixed states of the RE and their averaging by the RE distribution. If the random environment is "very fast", approximation called below as "mixed parameters" can be successfully applied. This approximation consists of averaging parameters of the arrival and service processes by the distribution of the $RE$ and calculation of performance measures in $BMAP/PH/N$ system with the averaged arrival and service rates.

In the fourth experiment, we show numerically that sometimes the described approximations make sense. However, in situations when environment is neither "very slow" nor "very fast", these approximations can be very poor. We consider the $RE$ s with different rate which are characterized by the generators of the form $Q^{(k)}=\left(\begin{array}{cc}-5& 5\\ 15& -15\end{array}\right){10}^k$

We vary the parameter $k$ from -7 to 4 what corresponds to the variation of the $RE$ rate from "very slow" to "very fast". In this and further experiments, the input flow is described by the $BMA{P}_1$ and $BMA{P}_2$ and the service process is defined by the $P{H}_1$ and $P{H}_2$ The results are presented in Figures 7, 8 and 9. In application of "mixed system" approximation, the averaged arrival rates under both states of the RE are equal to 3.488. The averaged service rate is equal to 2 at the first state of the $RE$ and is equal to 14 at the second state. The mean sojourn times of an arbitrary customer at these states are equal to 3.5297 and 0.0998, respectively; the probabilities of immediate access to the servers are equal to 0.2021 and 0.81399; the mean numbers of customers in the system are equal to 12.311 and 0.3482. The averaged, according to the stationary distribution of the $RE$ mean sojourn time of an arbitrary customer is equal to 2.6722 and the probability that an arbitrary arriving customer sees an idle server is equal to 0.355. In application of "mixed parameters" approximation, the averaged, according to the stationary distribution of the $RE$ arrival rate is equal to 3.778 while the averaged service rate is equal to 4.0625. The value of the mean sojourn time of an arbitrary customer in the system with averaged arrival and service rates is equal to 0.7541, the probability of immediate access to the servers is equal to 0.4168, the mean number of customers in the system is equal to 2.849.

Figures 7, 8 and 9 confirm the hypothesis that the first type approximation ("mixed system") is good in case of "very slow" $RE$ and the second one ("mixed parameters") can be applied to case of "very fast" $RE$ But sometimes the second type approximation is not very good (see Figures 8 ) because it is not quite clear how to make averaging of service intensity. Simple averaging of service rates under the different states of the $RE$ may be not correct when the load of the system is not high because there are time intervals when the system is empty and no service is provided. It is worth to note also that there is an interval for $RE$ rate (interval $k\in \left[-\mathrm{3,}0\right]$ ) where one should not use the values of the system performance measures calculated based on the considered approximating models. The use of these values can lead to the large relative error. Thus, Figures 7, 8 and 9 confirm the importance of investigation implemented in this chapter. Simple engineering approximations can lead to unsatisfactory performance evaluation and capacity planning in real world systems.

Figure 7.

Mean sojourn time of an arbitrary customer as a function of the $RE$ rate

Figure 8.

Probability of immediate access to the servers as a function of the $RE$ rate

Figure 9.

The mean number of customers in the system $L_{queue}$ as a function of the $RE$ rate

In the fifth experiment we compare the mean number of customers, probability of immediate access to the servers and loss probability in the $BMAP/PH/N$ and $BMAP/PH/N/L$ systems operating in the $RE$ for different values $L$ of the buffer capacity and different customers admission discipline.

Figure 10.

Mean number of customers in the system as a function of the buffer capacity $L$

Figure 11.

Probability of immediate access to the servers as a function of the buffer capacity $L$

Figure 12.

Loss probability as a function of the buffer capacity $L$

Looking at Figures 10-12, it should be noted that the rate of convergence of the curves corresponding to the disciplines $PA$ and $CR$ to their limits defined by the system with an infinite buffer is not very high. When we further increase the value $L$ , we discover that even for the buffer capacity $L$ about 5000, the difference is not negligible. So, estimation of performance measures of the system with an infinite buffer by the respective measures of the system with a finite buffer can be not very good. This explains why we made the separate analysis of the system with an infinite buffer.

Finally, in the sixth experiment we consider the next optimization problem:

$$J\left(N\right)=c_1{\overline{\nu }}_a+c_{2}N\to \underset{N}{max}$$ E104

where

${\overline{\nu }}_a$ is the mean sojourn time in the system, $N$ is the number of servers, $c_1$ is the charge for an unit of customer sojourn time in the system, $c_2$ is the cost of a server maintenance per unit of time.

It is clear that this problem is not trivial. When the number of servers is small, the cost of servers maintenance is also small, but the mean sojourn time is large. If we increase the number of servers, the mean sojourn time decreases while the cost of servers maintenance increases.

Let us assume that the cost coefficients be fixed as $c_1=5$ and $c_2=3$ Service time distribution at both states of the $RE$ is exponential with intensities ${\mu }^{(1)}=1$

$${\mu }^{(2)}=7$$ E105

On Figure 13, dependence of the cost criterion $J\left(N\right)$ on $N$ is presented along with the dependences of the summands $c_1{\overline{\nu }}_a$ and

$$c_{2}N$$ E106

Figure 13.

Criterion $J\left(N\right)$ as a function of number of servers in the system

Based on Figure 13 , one can conclude that our analysis allows effectively solve the problems of the system design and that the optimal value of the cost criterion (in this example it is provided by $N=4$ ) can be significantly smaller than the values of the cost criterion for other values on $N$

7. Conclusion

The $BMAP/PH/N/L$ system operating in a finite state space Markovian random environment is investigated for the finite and infinite buffer capacity. The joint stationary distribution of the number of the customers in the system, the state of the random environment, and the states of the underlying processes of arrival and service processes is calculated. The analytic formulas for performance measures of the system are derived. The Laplace-Stieltjes transform of sojourn time distribution is derived and the mean sojourn time is calculated. Selected results of numerical study are presented. They show an impact of the correlation in arrival process, illustrate the poor quality of the system characteristics approximation by means of more simple models, confirm the positive effect of the traffic redistribution between the peak and normal operation periods. The results can be used for the optimal design, capacity planning, and performance evaluation of real world systems in which operation of the system can be changed depending on some external factors.

Acknowledgments

This work was supported by the Korea Research Foundation Grant Funded by the Korean Government (MOEHRD)(KRF-2008-313-D01211).