Randomized Threshold Strategy for Providing Flexible Priority in Multi-Server Queueing System with a Marked Markov Arrival Process and Phase-Type Distribution of Service Time

Abstract

In this paper, we analyze a multi-server queueing system with a marked Markov arrival process of two types of customers and a phase-type distribution of service time depending on the type of customer. Customers of both types are assumed to be impatient and renege from the buffers after an exponentially distributed number of times. The strategy of flexible provisioning of priorities is analyzed. It assumes a randomized choice of the customers from the buffers, with probabilities dependent on the relation between the number of customers in a priority finite buffer and the fixed threshold value. To simplify the construction of the underlying Markov chain and the derivation of the explicit form of its generator, we use the so-called generalized phase-type distribution. It is shown that the created Markov chain fits the category of asymptotically quasi-Toeplitz Markov chains. Using this fact, we show that the considered Markov chain is ergodic for any value of the system parameters and compute its stationary distribution. Expressions for key performance measures are presented. Numerical results that show how the parameters of the control strategy affect the system’s performance measurements are given. It is shown that the results can be used for managerial purposes and that it is crucial to take correlation in the arrival process into account.

1. Introduction

Queueing theory is a rapidly developing branch of applied probability that helps solve problems related to the sharing and scheduling of finite resources that have to be used by competing users. Systems with heterogeneous customers have been investigated less in the existing queueing literature. The varying requirements of several types of customers for the required system resource, indicators of the quality of service, and their different importance and value all lead to their needing different treatments. In particular, a powerful tool for the enhancement of the operation of such systems is the assignment of different priorities in accessing the system resources for different types of customers. The theory of priority queues is pretty well developed in the case of single-server queues (see, e.g., [1,2,3,4,5]). However, the operation of many real-world systems can be adequately described only by multi-server priority queueing systems. The problem of optimal management in such systems is much more complicated than in single-server systems due to the larger dimensionality of the stochastic process that describes the dynamics of the system. However, the importance of consideration of such priority queueing systems made them a popular subject of research; see, e.g., the papers [6,7,8,9,10,11,12,13,14,15,16] published within the last two years.

There are a lot of various priority schemes (preemptive and non-preemptive, alternating, static and dynamic, changing and accumulating, etc.) known in the literature; some references can be found, e.g., in [17]. All these schemes have their own advantages and disadvantages. For example, the static priorities assuming strict adherence to the established rules of picking up from the queue are more easily realized in real systems. The preemptive priorities suggesting service interruption in the case of the arrival of a customer who has a higher priority than the one receiving service are fine for high-priority customers but are discriminatory with respect to low-priority customers. The dynamic priorities (see, e.g., [18,19],) give more flexibility but require efforts to monitor the system states. Therefore, to better manage the operation of different real-world systems, different new priority schemes are generated in the literature. For example, the scheme proposed in [17] suggests the maintenance of auxiliary finite input buffers with different rates of customer transfer to the main buffer. A proper choice of the capacity of the buffers and the corresponding rates allows one to reach any desired degree of non-preemptive priority for one class of customers over another. Another flexible scheme for providing non-preemptive priority was offered in [20]. That scheme assumes the randomized choice of the next customer for service, with probabilities of choice depending on the length of the queue of the priority customers. The scheme was applied and analyzed for the single-server queue. In this paper, we extend the analysis to the case of a multi-server queue. It is worth noting that there exist works (see, e.g., [21]) where a priority to certain flows, e.g., flows of handover users in cellular mobile networks or primary users in cognitive radio systems, is provided via different kinds of schemes of reservation of some servers of the multi-server queueing system for service of priority users. Here, we do not touch on such a possibility. Priority is provided via the proper strategy of customers picking up from the queues.

It is already well known that flows of customers receiving service in many telecommunication systems, contact centers, hospitals, etc., exhibit a so-called bursty nature. The instantaneous arrival rate may significantly fluctuate during system operation. This makes the very popular-in-the-literature, stationary Poisson model of the arrival process a poor descriptor of real flows. The use of this model is not suitable for reliable prediction of the main performance measures of the systems. It may lead to huge errors in the estimation of the required system resources sufficient for providing the desired quality of service. Therefore, in this paper, we assume that the arrival flows of two types of customers receiving service in the system are defined by the much more general Marked Markov Arrival Process ($MMAP$), see, e.g., [22,23,24]. The $MMAP$ is an extension of the Markov Arrival Process ($MAP$) for accounting for the possibility of heterogeneous customer arrivals. In contrast to the stationary Poisson arrival process, the $MAP$ is suitable for modeling flows with a non-zero coefficient of correlation and high variation in inter-arrival times. For more information about the basic properties, particular cases, and possible applications, as well as about the analysis of the systems with the $MAP$, see, e.g., [24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44]. The multi-server priority queues with $MMAP$ arrival process have been analyzed, e.g., in [9,13,16,17,20,45,46,47].

The model considered in this paper assumes a so-called phase-type ($PH$) distribution of service times of both types (see, e.g., [24,48,49,50,51,52]). This distribution essentially generalizes the well-known exponential distribution but still allows us to obtain nice analytical results. It can be applied to approximate any non-negative random variable distribution. In particular, it allows taking into account the variance of a random variable while the exponential distribution has a coefficient of variation equal to 1. However, the use of this distribution instead of the exponential distribution for the description of, e.g., service time, leads to the necessity of extending the dimension of the state space of the Markovian process describing the dynamics of a queueing system. It is necessary to keep track of the state of the underlying process (phase) of the $PH$ distribution in each busy server. If the considered queueing system has many servers, it essentially complicates the analysis of the system. If the number of identical servers is equal to N while the number of phases is equal to $M,$ then the direct way of accounting for the phase of the underlying processes of service on all servers leads to the necessity of operating with matrices of size $M^N$ describing the dynamics of the service process. When at least one of the numbers N and M is relatively large (while the number N of servers in real systems can be pretty large), operating with matrices of such a size can be difficult. This is why in this paper, we use another way to monitor the service process. This way assumes an account of the number of servers providing service at each of M possible phases. In the case when M is relatively small, the state space of the service processes becomes much smaller, namely, $\left(\genfrac{}{}{0pt}{}{N+M-1}{M-1}\right)=\frac{(N+M-1)!}{N!(M-1)!}.$ For example, in many cases, the value of M equal to 2 allows us to have a good catch on the variance of service time. For $M=2,$ the size of the matrices that have to be treated is equal to $N+1$, which may be much less than $M^N.$ This way of describing the service process in multi-server queues was offered in [53,54].

In those papers, the recursive formulas for the computation of the matrices defining transition rates of the multi-dimensional Markov process, which describes the service process in many servers without a service completion and at the moments of new service beginning or ending, are presented. Formulas for computation of these matrices for the value n of the number of currently busy servers (among N available servers), $n=\overline{1,N},$ depend on both n and $N.$ This is not convenient, especially when it is required to make computations for several values of N, e.g., in the process of computing the minimum required number N of servers to guarantee the fixed values of the service quality indicators. In [55], the formulas derived in [53,54] are modified to eliminate the use of the value of N in recursive computations.

In the considered model, it is assumed that the parameters of $PH$ distribution of the service time of customers of different types may be different. To avoid the necessity of separately monitoring the number of each type of customer in service and make the analysis easier, we use the notation of the generalized phase-type distribution introduced in [56,57].

We take into consideration the possibility of customers departing from the buffers due to impatience. Many real-world systems are inherently characterized by the impatience of their customers. For example, customers (information units) in telecommunication networks can depart from the waiting area due to information obsolescence, users’ departure from the service zone, etc. In contact centers, impatience is explained by the psychology of users. In retail networks and health systems (such as blood banks or organ transplantation hospitals), impatience is explained by the restricted term of the suitability of the items required for service. Recent lists of related research on queues with impatient customers can be found, e.g., in [58,59,60,61].

Account of customers’ impatience makes the considered model potentially useful; beyond the evident applications in telecommunications, it is also useful in the following systems. When two competing flows of passengers or vehicles, which arrive from different directions, merge into one flow, a popular mechanism of merging implementation consists of admitting a certain percentage of users from one direction and another. The considered model can help answer the question of how to optimally redistribute the percentage when one of the queues becomes large. In call centers, agents simultaneously handle the voice requests, which are admitted via the buffer of a finite capacity, as well as the requests arriving via the Internet that are placed in the buffer of a potentially infinite capacity. Over time, agents select for service waiting voice requests and Internet requests in some proportion. However, when the buffer for voice requests becomes close to its capacity (exceeds some threshold) and there is a high risk of their loss, the agents can decrease the percentage of Internet users admitted for service and, correspondingly, increase the percentage of voice users.

The structure of this paper is as follows. In Section 2, the model of the considered queueing system is completely defined. In Section 3, a multi-dimensional Markov chain describes how the system operates. The generator of this chain is derived there, and the problems of the existence of the stationary distribution of this Markov chain and the computation of this distribution are briefly discussed. In Section 4, formulas for the computation of some performance measures of the system are presented. Some numerical illustrations are presented in Section 5. The paper is concluded in Section 6.

2. Mathematical Model

We consider an N-server queuing system with two buffers, the structure of which is displayed in Figure 1.

The system processes customers of two types. The customer’s arrival in the system is described by the $MMAP$ which is defined by the underlying continuous-time Markov chain ($MC$) ${\nu }_t,t\ge 0,$ with the state space $\{1,2,\cdots ,W\}$, and the square matrices $D_0,D_1,$ and $D_2$ of size $W.$ The negative diagonal entries of the matrix $D_0$ define, up to the sign, the transition rates of the $MC$${\nu }_t,t\ge 0,$ from the corresponding states. The non-diagonal entries of the matrix $D_0$ define the transition rates that do not lead to the arrival of any customer. The entries of the matrices $D_l,l=1,2,$ are the transition rates of the $MC$${\nu }_t,t\ge 0,$ that cause the arrival of a type-l customer. More information about the $MMAP$ and its characteristics can be found, e.g., in [24].

Let us denote the total customers’ arrival rate as $\lambda$ and type-l customers’ arrival rate as ${\lambda }_l,l=1,2.$ These intensities are determined by the formulas:

$$\lambda =\mathit{\theta }(D_1+D_2)\mathbf{e},{\lambda }_1=\mathit{\theta }{D}_1\mathbf{e},{\lambda }_2=\mathit{\theta }{D}_2\mathbf{e}$$

where the stationary probability vector $\mathit{\theta }$ of process ${\nu }_t$ is determined to be the one solution to the system

$$\mathit{\theta }(D_0+D_1+D_2)=\mathbf{0},\mathit{\theta }\mathbf{e}=1$$

where the row vector $\mathbf{0}$ consists of 0 s and the column vector $\mathbf{e}$ consists of 1 s.

The service time of a type-$l,l=1,2,$ customer has the phase-type distribution $P{H}_l$ that is defined by the underlying continuous-time $MC$$m_t^{\left(l\right)},t\ge 0,$ with the state space $\{1,2,\cdots ,M_l\}$, the vector ${\mathbf{b}}_l$ and the matrix $S_l$. The stochastic row vector ${\mathbf{b}}_l$ determines the initial state of the process $m_t^{\left(l\right)}$ at the service beginning epoch, and the subgenerator $S_l$ defines rates of transition of the process $m_t^{\left(l\right)}$ inside the state space $\{1,2,\cdots ,M_l\}.$ The service time finishes when the process $m_t^{\left(l\right)}$ transits to the absorbing state. The mean service time of type-l customers is calculated as ${\mathbf{b}}_l{(-S_l)}^{-1}\mathbf{e}.$ For further details on the $PH$ distribution, its moments, and specific cases, see, for example, [24].

If a customer of any type arrives at the system when there is a free server, this customer immediately starts service. A type-2 customer joins the infinite buffer if there are no available servers during its arrival period. Let us call this buffer Buffer-2. If there are no free servers during the arrival epoch of a type-1 customer, this customer tries to join the finite buffer of capacity R (Buffer-1). This trial will be successful in the case of a not full buffer. The arriving customer is lost forever if Buffer-1 is full during the type-1 customer arrival moment.

If during the epoch when some server finishes the service only customers of one type are in the buffers, the customer of such a type is chosen for service. If there are customers of two types in the buffers at the service completion moment, the type of customer that will receive service next is defined according to the following rule. If the number of type-1 customers in the buffer is less than the parameter K, a type-1 customer is chosen for service with the probability $q_1$, and with the complimentary probability $1-q_1$, a type-2 customer starts service. If the number of type-1 customers in the buffer is greater or equal to K, a type-1 customer is chosen for service with the probability $q_2,q_2\ge q_1,$ and with the complimentary probability $1-q_2$ a type-2 customer starts service. We assume that the parameters $K,$ $q_1$, and $q_2$ are control parameters, and the task of modeling is to develop a mathematical tool that allows us to find the optimal values of these parameters according to a fixed cost criterion.

Each type of customer is assumed to be impatient and depart from the system without service if the service does not start within an exponentially distributed time since the arrival moment. The total rate of customers’ reneging from the corresponding buffer depends on the number of customers staying in the buffer. We assume that if there are r customers in Buffer-1, the total rate of reneging by type-1 customers is equal to ${\alpha }_r,{\alpha }_r\ge 0,r=\overline{1,R}.$ If there are i customers in the infinite Buffer-2, the total rate of type-2 customers reneging is equal to ${ϵ}_i,{ϵ}_i\ge 0,$ and ${ϵ}_i\to \infty$ if $i\to \infty .$

Let us analyze the stationary behavior of the model under consideration.

3. The Process of the System States and Its Analysis

The operation of the system under study is described by the following regular irreducible $(M_1+M_2+4)$-dimensional continuous-time $MC$

$${\zeta }_t=\{i_t,r_t,n_t,{\nu }_t,{\sigma }_t^{\left(1\right)},\cdots ,{\sigma }_t^{\left(M_1\right)},{\gamma }_t^{\left(1\right)},\cdots ,{\gamma }_t^{\left(M_2\right)}\},t\ge 0,$$

where at time $t,t\ge 0,$

•

$i_t$ is the number of customers in the system, $i_t\ge 0;$

•

$r_t$ is the number of type-1 customers in Buffer-1, $r_t=\overline{0,\min \left\{\max \right\{0,i-N\},R\}}$;

•

$n_t$ is the number of type-1 customers in service, $n_t=\overline{0,\min \{N,i_t\}}$;

•

${\nu }_t$ is the state of the underlying process of the $MMAP$ of customers, ${\nu }_t=\overline{1,W};$

•

${\sigma }_t^{\left(m\right)}$ is the number of type-1 customers at the m-th phase of the $P{H}_1$ service process, ${\sigma }_t^{\left(m\right)}=\overline{0,n_t},m=\overline{1,M_1},{\displaystyle \sum _{m=1}^{M_1}}{\sigma }_t^{\left(m\right)}=n_t;$

•

${\gamma }_t^{\left(l\right)}$ is the number of type-2 customers at the l-th phase of the $P{H}_2$ service process, ${\gamma }_t^{\left(l\right)}=\overline{0,\min \{N,i_t-n_t\}},l=\overline{1,M_2},{\displaystyle \sum _{l=1}^{M_2}}{\gamma }_t^{\left(l\right)}=\min \{N,i_t-n_t\}.$

One can see that the process ${\zeta }_t$ is complicated, and its analysis may be cumbersome. To make the investigation easier, we use the notion of the generalized phase-type distribution, see, [56,57]. Instead of the separate consideration of service processes of each type of customer, we suggest that the customer’s service time at a server has the generalized $PH$ ($GPH$) distribution having the parameters $({\mathit{\beta }}_1,{\mathit{\beta }}_2,S)$ where the vectors ${\mathit{\beta }}_1$ and ${\mathit{\beta }}_2$ and subgenerator S are defined by the formulas:

$${\mathit{\beta }}_1=({\mathbf{b}}_1,\underset{M_2}{\underbrace{0,0,\dots ,0}}),{\mathit{\beta }}_2=(\underset{M_1}{\underbrace{0,0,\dots ,0}},{\mathbf{b}}_2),S=\left(\begin{array}{cc}{S}_1& O\\ O& S_2\end{array}\right).$$

The duration of service is governed by the underlying process $m_t$ with the state space $\{1,2,\cdots ,M\},$ where $M=M_1+M_2.$ The state of the process $m_t$ during the service beginning epoch is defined according to the probabilistic vector ${\mathit{\beta }}_l,l=1,2,$ if a type-l customer starts service. The transition rates into the absorbing state (which implies service completion on one of the busy servers) are given by the entries of the column vector ${\mathbf{S}}_0=\left(\begin{array}{c}-S_1\mathbf{e}\\ -S_2\mathbf{e}\end{array}\right).$

The use of $GPH$ allows us to significantly simplify the process of defining the dynamics of the system and, correspondingly, its analysis.

The following regular irreducible continuous-time $MC$ can be used to describe the behavior of the system under consideration:

$${\xi }_t=\{i_t,r_t,{\nu }_t,{\eta }_t^{\left(1\right)},\cdots ,{\eta }_t^{\left(M\right)}\},t\ge 0,$$

The meaning of the first three components $\{i_t,r_t,{\nu }_t\}$ is the same as for the $MC$${\zeta }_t$. The component ${\eta }_t^{\left(m\right)}$ of the vector-valued process ${\mathit{\eta }}_t=\{{\eta }_t^{\left(1\right)},\cdots ,{\eta }_t^{\left(M\right)}\},t\ge 0,$ denotes the number of customers at the m-th phase of the $GPH$ service process at the moment t, ${\eta }_t^{\left(m\right)}=\overline{0,\min \{i_t,N\}},m=\overline{1,M},{\displaystyle \sum _{m=1}^M}{\eta }_t^{\left(m\right)}=\min \{i_t,N\}.$

To describe the transition rates of the $MC$${\xi }_t,$ let us first introduce the following auxiliary matrices:

$A_i=A_i\left(S\right),i=\overline{1,N},$ define the rates of transitions of the process ${\mathit{\eta }}_t,t\ge 0,$ which do not lead to the end of service provided that i servers are busy;

$P_i\left({\mathit{\beta }}_1\right),P_i\left({\mathit{\beta }}_2\right),i=\overline{0,N-1},$ define the transition probabilities of the process ${\mathit{\eta }}_t,t\ge 0,$ when i servers are busy and a type-1 or type-2 customer is chosen for service, correspondingly;

$L_i=L_i\left({\mathbf{S}}_0\right),i=\overline{1,N},$ define the transitions rates of the process ${\mathit{\eta }}_t,t\ge 0,$ which cause the finish of service if i servers are busy;

${\mathsf{\Delta }}_i,i=\overline{1,N},$ define the rates of exit of the process ${\mathit{\eta }}_t,t\ge 0,$ from its states when i customers are in service.

Note that the recursive algorithms for the computation of the matrices $A_i,L_i,$${\mathsf{\Delta }}_i,i=\overline{1,N},$ and $P_i\left({\mathit{\beta }}_1\right),P_i\left({\mathit{\beta }}_2\right),i=\overline{0,N-1},$ are available in paper [55] and represent the advanced modification of the algorithms earlier presented in [53,54] and used, e.g., in [62,63,64].

We enumerate the components $\{i_t,r_t,{\nu }_t\}$ of the $MC$${\xi }_t,t\ge 0,$ in the direct lexicographic order and the components $\{{\eta }_t^{\left(1\right)},\cdots ,{\eta }_t^{\left(M\right)}\}$ in reverse lexicographic order. All states with the value i of the first component $i_t$ are called the level i of the $MC$${\xi }_t.$

We refer to the chain’s infinitesimal generator as G. All possible transition rates of the chain under study during an infinitesimally short time interval are defined as the entries of the generator G.

The notation used up until this point, and the new notation are compiled in

Table 1. Notation.

N	the number of servers
$D_0,D_1,D_2$	the $W\times W$ matrices that define the $MMAP$
${\lambda }_l,l=1,2$	the average arrival rate of type-l customers
$\mathbf{0}$	row vector which consists of 0s
$\mathbf{e}$	column vector which consists of 1s
$({\mathit{\beta }}_l,S_l),l=1,2$	an irreducible representation of the type-l customer $PH$ service
	time distribution
R	Buffer-1 capacity
K, $q_1$ and $q_2$	control parameters that define the choice of customers for service
${\alpha }_r,r=\overline{1,R}$	the total rate of reneging by type-1 customers if there are
	r customers in Buffer-1
${ϵ}_i,i\ge 1$	the total rate of reneging by type-2 customers if there are
	i customers in Buffer-2
$({\mathit{\beta }}_1,{\mathit{\beta }}_2,S)$	the parameters that define $GPH$ service
	time distribution with M states
S0	$\left(\begin{array}{c}-S_1\mathbf{e}\\ -S_2\mathbf{e}\end{array}\right)$
$A_i=A_i\left(S\right),i=\overline{1,N}$	the matrix that defines the rates of transitions that do not lead
	to the end of service provided that i servers are busy
$P_i\left({\mathit{\beta }}_1\right),P_i\left({\mathit{\beta }}_2\right),$	the matrix that defines the transition probabilities of the process
$i=\overline{0,N-1}$	${\eta }_t,t\ge 0,$ when i servers are busy and a type-1 or type-2
	is chosen for service, correspondingly
$L_i=L_i\left({\mathbf{S}}_0\right),$	the matrix that defines the transition rates that cause the finish
$i=\overline{1,N}$	of service if i servers are busy
${\mathsf{\Delta }}_i,i=\overline{1,N}$	the matrix that defines the rates of exit of the service
	process from its states when i customers are in service
$T_N$	$\left(\genfrac{}{}{0pt}{}{N+M-1}{M-1}\right)=\frac{(N+M-1)!}{N!(M-1)!}$
I	the identity matrix
O	a zero matrix
⊗ and ⊕	the symbols of Kronecker product and sum of matrices,
	see, e.g., [65,66,67]
$\mathrm{diag}\{a_1,\cdots ,a_l\}$	a block-diagonal matrix with the diagonal blocks $a_1,\cdots ,a_l$
$C_l^{\left(1\right)}$	$\mathrm{diag}\{0,{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_l\},l=\overline{1,R-1}$
$C_l^{\left(2\right)}$	$\mathrm{diag}\{{ϵ}_l,{ϵ}_{l-1},\cdots ,{ϵ}_1,0\},l=\overline{1,R-1}$
C	$\mathrm{diag}\{0,{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_R\}$
${\overline{C}}_l$	$\mathrm{diag}\{{ϵ}_l,{ϵ}_{l-1},\cdots ,{ϵ}_{l-R-1},{ϵ}_{l-R}\},l\ge R$
$E_l^{+},l=\overline{0,R-1}$	an $(l+1)\times (l+2)$ matrix having all zero
	elements except the elements ${\left(E_l^{+}\right)}_{k,k+1},k=\overline{0,l},$
	which are equal to 1
${\overline{E}}_l^{+},l=\overline{0,R-1}$	an $(l+1)\times (l+2)$ matrix having all zero elements except the
	elements ${\left({\overline{E}}_l^{+}\right)}_{k,k},k=\overline{0,l},$ which are equal to 1
$E^{+}$	a square matrix of size $R+1$ having all zero elements except the
	elements ${\left(E^{+}\right)}_{k,k+1},k=\overline{0,R-1},$ which are equal to 1
$\widehat{I}$	a square matrix of size $R+1$ having all zero elements except
	the one ${\left(\widehat{I}\right)}_{R,R}=1$
$E_l^{-},l=\overline{1,R}$	an $(l+1)\times l$ matrix having all zero elements except the
	element ${\left(E_l^{-}\right)}_{k,k-1},k=\overline{1,l},$ which are defined as follows:
	${\left(E_l^{-}\right)}_{k,k-1}=\left\{\begin{array}{cc}{q}_1,& ifk\in \{1,2,\cdots ,K-1\},k\ne l,\\ q_2,& ifk\in \{K,K+1,\cdots ,l-1\}\\ 1,& ifk=l.\end{array}\right.$
${\overline{E}}_l^{-},l=\overline{1,R}$	an $(l+1)\times l$ matrix having all zero elements except the
	elements ${\left({\overline{E}}_l^{-}\right)}_{k,k},k=\overline{0,l-1},$ which are defined as follows:
	${\left({\overline{E}}_l^{-}\right)}_{k,k}=\left\{\begin{array}{cc}1,& ifk=0,\\ 1-q_1,& ifk\in \{1,2,\cdots ,K-1\},\\ 1-q_2,& ifk\in \{K,K+1,\cdots ,l-1\},\end{array}\right.$
$I_l^{-},l=\overline{1,R}$	an $(l+1)\times l$ matrix having all zero elements except the elements
	${\left(I_l^{-}\right)}_{k,k-1},k=\overline{1,l},$ which are equal to 1;
${\overline{I}}_l^{-},l=\overline{1,R}$	an $(l+1)\times l$ matrix having all zero elements except the
	elements ${\left({\overline{I}}_l^{-}\right)}_{k,k},k=\overline{0,l-1},$ which are equal to 1;
$E^{-}$	a square matrix of size R + 1 having all zero elements except the elements ${\left(E^{-}\right)}_{k,k-1},k=\overline{1,R},$ which are defined as follows:
	${\left(E^{-}\right)}_{k,k-1}=\left\{\begin{array}{cc}{q}_1,& ifk\in \{1,2,\cdots ,K-1\},\\ q_2,& ifk\in \{K,K+1,\cdots ,R\}.\end{array}\right.$
${\overline{E}}^{-}$	a square matrix of size $R+1$ having all zero elements except the
	elements ${\left({\overline{E}}^{-}\right)}_{k,k},k=\overline{0,R-1},$ which are defined as follows:
	${\left({\overline{E}}^{-}\right)}_{k,k}=\left\{\begin{array}{cc}1,& ifk=0,\\ 1-q_1,& ifk\in \{1,2,\cdots ,K-1\},\\ 1-q_2,& ifk\in \{K,K+1,\cdots ,R\},\end{array}\right.$
$I^{-}$	a square matrix of size $R+1$ having all zero elements except the
	elements ${\left(I^{-}\right)}_{k,k-1},k=\overline{1,R},$ which are equal to 1

for clarity and to make the description of the shape of the blocks of the generator G simpler.

Theorem 1.

The generator G has the following block-three-diagonal form:

$$G=\left(\begin{array}{cccccc}{G}_{0,0}& G_{0,1}& O& O& O& \dots \\ G_{1,0}& G_{1,1}& G_{1,2}& O& O& \dots \\ O& G_{2,1}& G_{2,2}& G_{2,3}& O& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right).$$

The non-zero blocks $G_{i,j},i,j\ge 0,$ determining the transition rates from level i to level j are defined as follows:

$$G_{0,0}=D_0,$$

$$G_{i,i}=D_0\oplus (A_i+{\mathsf{\Delta }}_i),i=\overline{1,N},$$

$$G_{i,i}=I_{i-N+1}\otimes (D_0\oplus (A_N+{\mathsf{\Delta }}_N))-C_{i-N}^{\left(1\right)}\otimes I_{W{T}_N}-C_{i-N}^{\left(2\right)}\otimes I_{W{T}_N},i=\overline{N+1,N+R-1},$$

$$G_{i,i}=I_{R+1}\otimes (D_0\oplus (A_N+{\mathsf{\Delta }}_N))-C\otimes I_{W{T}_N}-{\overline{C}}_{i-N}\otimes I_{W{T}_N}+\widehat{I}\otimes D_1\otimes I_{T_N},i\ge N+R,$$

$$G_{i,i+1}=D_1\otimes P_i\left({\mathit{\beta }}_1\right)+D_2\otimes P_i\left({\mathit{\beta }}_2\right),i=\overline{0,N-1},$$

$$G_{i,i+1}=E_{i-N}^{+}\otimes D_1\otimes I_{T_N}+{\overline{E}}_{i-N}^{+}\otimes D_2\otimes I_{T_N},i=\overline{N,N+R-1},$$

$$G_{i,i+1}=E^{+}\otimes D_1\otimes I_{T_N}+I_{R+1}\otimes D_2\otimes I_{T_N},i\ge N+R,$$

$$G_{i,i-1}=I_W\otimes L_i,i=\overline{1,N},$$

$$G_{i,i-1}=E_{i-N}^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_1\right)+{\overline{E}}_{i-N}^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_2\right)+$$

$$+C_{i-N}^{\left(1\right)}{I}_{i-N}^{-}\otimes I_{W{T}_N}+C_{i-N}^{\left(2\right)}{\overline{I}}_{i-N}^{-}\otimes I_{W{T}_N},i=\overline{N+1,N+R},$$

$$G_{i,i-1}=E^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_1\right)+{\overline{E}}^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_2\right)+$$

$$+C{I}^{-}\otimes I_{W{T}_N}+{\overline{C}}_{i-N}\otimes I_{W{T}_N},i>N+R.$$

Proof. 

By analyzing every potential transition of the $MC$$x{i}_t$ during an interval of infinitesimal length and reconstructing the rates of such transitions in block matrix form, the proof of the theorem is implemented.

Customers of both types arrive and are served one at a time. The possibility of more than one consumer entering the system or being served in a moment that lasts an infinite amount of time is negligibly small. Hence, the generator G has a block-tridiagonal structure, and all the matrices $G_{i,j}$ are zero ones for i and j, such as $|i-j|>1.$

The negative diagonal entries of the matrices $G_{i,i},i\ge 0,$ define, up to the sign, the total rate of leaving the corresponding state. The $MC$${\xi }_t$ can exit from its state in the following cases:

(1) The underlying arrival process changes its state. The rates of the arrival process transitions are defined by the modulus of the corresponding diagonal entries of the matrix $D_0$ if $i=0$ (the system is empty), $D_0\otimes I_{T_i}$ for $i=\overline{1,N},$ and $I_{\min \left\{\max \right\{0,i-N\},R\}+1}\otimes D_0\otimes I_{T_N}$ for $i\ge N+1.$ Note that if the underlying arrival process makes a transition from some state to the same state with a type-1 customer arrival when Buffer-1 is full, the customer is lost, and the exit from the corresponding state does not occur. The rates of these transitions are given by the corresponding diagonal elements of the matrix $\widehat{I}\otimes D_1\otimes I_{T_N},i\ge N+R.$

(2) The transition in the service process on one of the servers if there are busy servers in the system. The rates of these transitions are defined by the modulus of the corresponding diagonal elements of the matrix $I_W\otimes {\mathsf{\Delta }}_i,i=\overline{1,N}$ (if i customers are served in the system and there is no customer in the buffers) and by the matrix $I_{\left(\min \right\{\max \{0,i-N\},R\}+1)W}\otimes {\mathsf{\Delta }}_N,i\ge N+1$ (if there are customers in service and in the buffers).

(3) One of the customers leaves the system due to impatience. If a type-1 customer leaves Buffer-1, the rates of these transitions are defined by the diagonal elements of the matrices $C_{i-N}^{\left(1\right)}\otimes I_{W{T}_N}$ for $i=\overline{N+1,N+R-1}$ and of the matrix $C\otimes I_{W{T}_N}$ for $i\ge N+R.$ The diagonal elements of the matrices $C_{i-N}^{\left(2\right)}\otimes I_{W{T}_N},i=\overline{N+1,N+R-1},$ and of the matrix ${\overline{C}}_{i-N}\otimes I_{W{T}_N},i\ge N+R$, are the transition rates in the case of a type-2 customer leaving Buffer-2 due to impatience.

The non-diagonal elements of the blocks $G_{i,i}$ determine the transition rates of the $MC$${\xi }_t$ that do not lead to the change in the number of customers in the system $i\ge 0.$ These transitions are as follows:

(1) The transition in the underlying arrival process without the arrival of a new customer (the rates of these transitions are the non-diagonal elements of the matrix $D_0$ if $i=0$, $D_0\otimes I_{T_i}$ for $i=\overline{1,N},$ and $I_{\min \left\{\max \right\{0,i-N\},R\}+1}\otimes D_0\otimes I_{T_N}$ for $i\ge N+1$).

(2) The arrival of a type-1 customer when N customers are being served in the system and Buffer-1 is full. Note that in such a case, the arriving customer is not accepted into the system. The corresponding rates are given by the elements of the matrix $\widehat{I}\otimes D_1\otimes I_{T_N},i\ge N+R.$

(3) The change of the state of the service underlying process without service completion. The rates of these transitions are defined by the non-diagonal elements of the matrices $I_W\otimes A_i,i=\overline{1,N}$ (if i servers are busy but the buffers are empty) and of the matrices $I_{\left(\min \right\{\max \{0,i-N\},R\}+1)W}\otimes A_N,i\ge N+1$ (if all servers are busy and there are customers in the buffers).

Using all these reasonings, we prove the form of the blocks $G_{i,i},i\ge 0.$

The elements of the matrices $G_{i,i+1},i\ge 0,$ define the transition rates of the components of the $MC$${\xi }_t$ that lead to the increase in the number of customers i in the system by one. This can happen only if a new customer is admitted to the system. When there is an idle server at an arrival moment, the transition rates are given by the elements of the matrices $D_l\otimes P_i\left({\mathit{\beta }}_l\right),i=\overline{0,N-1},$ if type-$l,$ $l=1,2,$ customer arrives. The Kronecker multiplier $P_i\left({\mathit{\beta }}_l\right)$ reflects the fact that the arriving customer immediately starts service and, therefore, installation of the initial state of its service underlying process is required. If all servers are busy and an arriving type-1 customer joins Buffer-1, the elements of the matrices $E_{i-N}^{+}\otimes D_1\otimes I_{T_N}$ define the corresponding rates. The Kronecker multiplier $P_i\left({\mathit{\beta }}_l\right)$ is replaced here by the Kronecker multiplier $I_{T_N}$ because the new service does not start, and no transition of the service underlying process occurs. If a type-2 customer arrives at the system and is appended to Buffer-2, the corresponding rates are defined by the matrices ${\overline{E}}_{i-N}^{+}\otimes D_2\otimes I_{T_N},i=\overline{N,N+R-1},$ and $I_{R+1}\otimes D_2\otimes I_{T_N},i\ge N+R.$ The matrices $E_{i-N}^{+}$ and ${\overline{E}}_{i-N}^{+}$ reflect the fact of increasing the number of customers of the corresponding type in the buffer.

Next, we explain the formation of the blocks $G_{i,i-1},i\ge 1.$ Their elements specify the transition rates of the components of the $MC$${\xi }_t$ that lead to the decrease in the number of customers i in the system by one. These transitions are the following:

(1) A customer leaves the system due to the successful completion of the service. In this case, the corresponding rates are given as the elements of the matrices: (i) $I_W\otimes L_i,i=\overline{1,N},$ if there are i customers in service and the buffers are empty; (ii) $E_{i-N}^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_1\right),i=\overline{N+1,N+R},$ and $E^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_1\right),i>N+R,$ if all servers are busy and there are customers in the buffers, and a type-1 customer from Buffer-1 is chosen for service. Arguing in a similar way, we obtain matrices ${\overline{E}}_{i-N}^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_2\right),i=\overline{N+1,N+R},$ and ${\overline{E}}^{-}\otimes I_W\otimes L_N{P}_{N-1}\left({\mathit{\beta }}_2\right),i>N+R,$ that are defined the corresponding rates if a type-2 customer from Buffer-2 goes to service.

(2) A customer leaves the buffer due to impatience. If a type-1 customer leaves the system, the corresponding rates are given by the matrices $C_{i-N}^{\left(1\right)}{I}_{i-N}^{-}\otimes I_{W{T}_N}$ for $i=\overline{N+1,N+R}$ and by the matrix $C{I}^{-}\otimes I_{W{T}_N}$ for $i>N+R.$ If a type-2 customer is lost due to impatience, the related rates are given by the matrices $C_{i-N}^{\left(2\right)}{\overline{I}}_{i-N}^{-}\otimes I_{W{T}_N},i=\overline{N+1,N+R},$ and ${\overline{C}}_{i-N}\otimes I_{W{T}_N},i>N+R.$ □

Having obtained the generator G of the $MC$${\xi }_t,$ it is possible to analyze this $MC.$

Theorem 2.

The $MC$${\xi }_t$ is ergodic for any choice of the system parameters.

Proof. 

It is easy to prove that, due to the assumption that ${ϵ}_i\to \infty$ when $i\to \infty ,$ the following limits exist:

$$Y_0=\underset{i\to \infty }{lim}{R}_i^{-1}{G}_{i,i-1}=I,$$

$$Y_1=\underset{i\to \infty }{lim}{R}_i^{-1}{G}_{i,i}+I=O,$$

$$Y_2=\underset{i\to \infty }{lim}{R}_i^{-1}{G}_{i,i+1}=O$$

where the matrix $R_i=-I⨀G_{i,i}.$ Here, ⨀ is the symbol of the Hadamard product of matrices, see [67]. Therefore, the $MC$${\xi }_t,t\ge 0,$ fits the category of asymptotically quasi-Toeplitz Markov chains ($AQTMC$), see the definition of $AQTMC$ in [68].

As follows from [68], a sufficient ergodicity condition of the $MC$${\xi }_t$ can be presented in the form

$$\mathbf{y}{Y}_0\mathbf{e}>\mathbf{y}{Y}_2\mathbf{e}$$

where the row-vector $\mathbf{y}$ can be obtained as the only solution

$$\mathbf{y}(Y_0+Y_1+Y_2)=\mathbf{y},\mathbf{y}\mathbf{e}=1.$$

Since $Y_0=I,Y_1=O,Y_2=O,$ one can see that the ergodicity condition turns to $1>0$ and is fulfilled for any system parameters. □

It follows from Theorem 2, that the following limits (stationary probabilities of the states of the $MC$${\xi }_t$) exist:

$$\pi (i,r,\nu ,m)=\underset{t\to \infty }{lim}P\{i_t=i,r_t=r,{\nu }_t=\nu ,{\eta }_t^{\left(1\right)}={\eta }^{\left(1\right)},\cdots ,{\eta }_t^{\left(M\right)}={\eta }^{\left(M\right)}\},$$

$$i\ge 0,r=\overline{0,\min \left\{\max \right\{0,i-N\},R\}},\nu =\overline{1,W},$$

$${\eta }^{\left(m\right)}=\overline{0,\min \{i,N\}},m=\overline{1,M},\sum _{m=1}^M{\eta }^{\left(m\right)}=\min \{i,N\}.$$

Let us form the row vectors ${\mathit{\pi }}_i,i\ge 0,$ of these probabilities enumerated in the same order as the components of the $MC$${\xi }_t$.

To obtain the stationary probability vectors ${\mathit{\pi }}_i,i\ge 0,$ we have to find the only solution to the following system:

$$({\mathit{\pi }}_0,{\mathit{\pi }}_1,\cdots )G=\mathbf{0},({\mathit{\pi }}_0,{\mathit{\pi }}_1,\cdots )\mathbf{e}=1$$

called the equilibrium or Chapman–Kolmogorov equations.

This system has an infinite number of unknowns and equations, and the generator G does not have a level-independent structure. Thus, the problem of solving it is not easy. To solve this system, we recommend using the numerically stable algorithm presented in [69]. In contrast to the general algorithm for $AQTMC$ proposed in [68] and designed for the upper-Hessenberg structure of the generator, the algorithm presented in [69] is oriented, namely, to solving Chapman–Kolmogorov equations with the block tridiagonal structure of the generator, which has the $MC$${\xi }_t$ under study.

4. System Performance Characteristics

The average number of customers in the system is computed by

$$N_{sys}=\sum _{i=1}^{\infty }i{\mathit{\pi }}_i\mathbf{e}.$$

The mean number of type-1 customers in Buffer-1 can be found as

$$N_{buf}^{\left(1\right)}=\sum _{i=N+1}^{\infty }\sum _{r=1}^{\min \{i-N,R\}}r\mathit{\pi }(i,r)\mathbf{e}.$$

The mean number of type-2 customers in Buffer-2 can be found as

$$N_{buf}^{\left(2\right)}=\sum _{i=N+1}^{\infty }\sum _{r=0}^{\min \{i-N-1,R\}}(i-N-r)\mathit{\pi }(i,r)\mathbf{e}.$$

The mean number of customers in the buffers is

$$N_{buf}=\sum _{i=N+1}^{\infty }(i-N){\mathit{\pi }}_i\mathbf{e}=N_{buf}^{\left(1\right)}+N_{buf}^{\left(2\right)}.$$

The loss probability of an arbitrary type-1 customer upon arrival due to Buffer-1 overflow can be found as

$$P_1^{ent-loss}=\frac{1}{{\lambda }_1}\sum _{i=N+R}^{\infty }\mathit{\pi }(i,R)(D_1\otimes I_{T_N})\mathbf{e}.$$

The loss probability of an arbitrary type-1 customer due to impatience in Buffer-1 can be found as

$$P_1^{imp-loss}=\frac{1}{{\lambda }_1}\sum _{i=N+1}^{\infty }\sum _{r=1}^{\min \{i-N,R\}}{\alpha }_r\mathit{\pi }(i,r)\mathbf{e}.$$

The loss probability of an arbitrary type-2 customer due to impatience in Buffer-2 can be found as

$$P_2^{imp-loss}=\frac{1}{{\lambda }_2}\sum _{i=N+1}^{\infty }\sum _{r=0}^{\min \{i-N-1,R\}}{ϵ}_{i-N-r}\mathit{\pi }(i,r)\mathbf{e}.$$

The loss probability of an arbitrary customer due to impatience can be found as

$$P^{imp-loss}=\frac{{\lambda }_1{P}_1^{imp-loss}+{\lambda }_2{P}_2^{imp-loss}}{\lambda }.$$

The rate of the output flow of successfully served type-1 customers can be found as

$${\lambda }_1^{out}=\sum _{i=1}^{\infty }{\mathit{\pi }}_i(I_{\left(\min \right\{\max \{0,i-N\},R\}+1)W}\otimes L_{\min \{i,N\}}^{\left(1\right)})\mathbf{e}.$$

The rate of the output flow of successfully served type-2 customers can be found as

$${\lambda }_2^{out}=\sum _{i=1}^{\infty }{\mathit{\pi }}_i(I_{\left(\min \right\{\max \{0,i-N\},R\}+1)W}\otimes L_{\min \{i,N\}}^{\left(2\right)})\mathbf{e}.$$

Remark 1.

The matrices $L_i^{\left(1\right)}=L_i^{\left(1\right)}\left({\mathbf{S}}_0^{\left(1\right)}\right)$ and $L_i^{\left(2\right)}=L_i^{\left(1\right)}\left({\mathbf{S}}_0^{\left(2\right)}\right),i=\overline{1,N},$ can be calculated using the algorithm from [55]. Their elements define the rates of transitions of the process ${\mathit{\eta }}_t,t\ge 0,$ which lead to the service completion of a type-1 and type-2 customer correspondingly on the condition that i servers are busy. The matrices ${\mathbf{S}}_0^{\left(l\right)},l=1,2,$ are given as 

$${\mathbf{S}}_0^{\left(1\right)}=\left(\begin{array}{c}-S_1\mathbf{e}\\ {\mathbf{0}}^T\end{array}\right),{\mathbf{S}}_0^{\left(2\right)}=\left(\begin{array}{c}{\mathbf{0}}^T\\ -S_2\mathbf{e}\end{array}\right).$$

The rate of the output flow of successfully served customers can be found as

$${\lambda }^{out}=\sum _{i=1}^{\infty }{\mathit{\pi }}_i(I_{\left(\min \right\{\max \{0,i-N\},R\}+1)W}\otimes L_{\min \{i,N\}})\mathbf{e}={\lambda }_1^{out}+{\lambda }_2^{out}.$$

The probability of an arbitrary type-1 customer loss can be found as

$$P_1^{loss}=1-\frac{{\lambda }_1^{out}}{{\lambda }_1}=P_1^{ent-loss}+P_1^{imp-loss}.$$

The probability of an arbitrary type-2 customer loss can be found as

$$P_2^{loss}=1-\frac{{\lambda }_2^{out}}{{\lambda }_2}=P_2^{imp-loss}.$$

The probability of an arbitrary customer loss can be found as

$$P^{loss}=1-\frac{{\lambda }^{out}}{\lambda }=\frac{{\lambda }_1{P}_1^{loss}+{\lambda }_2{P}_2^{loss}}{\lambda }.$$

The probability that the system is empty at an arbitrary moment can be found as

$$P_{idle}={\mathit{\pi }}_0\mathbf{e}.$$

The probability that an arbitrary type-1 customer starts service immediately upon arrival can be found as

$$P_1^{imm}=\frac{1}{{\lambda }_1}\sum _{i=0}^{N-1}{\mathit{\pi }}_i(D_1\otimes I_{T_i})\mathbf{e}.$$

The probability that an arbitrary type-2 customer starts service immediately upon arrival can be found as

$$P_2^{imm}=\frac{1}{{\lambda }_2}\sum _{i=0}^{N-1}{\mathit{\pi }}_i(D_2\otimes I_{T_i})\mathbf{e}.$$

The probability that an arbitrary customer starts service immediately upon arrival can be found as

$$P^{imm}=\frac{1}{\lambda }\sum _{i=0}^{N-1}{\mathit{\pi }}_i[(D_1+D_2)\otimes I_{T_i}]\mathbf{e}=\frac{{\lambda }_1{P}_1^{imm}+{\lambda }_2{P}_2^{imm}}{\lambda }.$$

5. Numerical Example

In the numerical examples, we investigate the influence of the control parameters $q_1$, $q_2$, and K on the system operation. Since we can only present three-dimensional graphs, first, we examine the impact of the probabilities $q_1$ and $q_2$ on the main performance characteristics of the system under the fixed value of K. Then, we examine the influence of the probability $q_1$ and threshold K under the fixed value of probability $q_2.$ Finally, we illustrate the importance of considering the system with the $MMAP.$ Namely, we illustrate that the approximation of the $MMAP$ by the superposition of two stationary Poisson arrival processes may lead to significant errors in the evaluation of the key performance measures of the system.

Experiment 1. Let us assume that the customers arrive at the system according to the $MMAP$ arrival process that is defined by the following matrices:

$$D_0=\left(\begin{array}{cc}-4.20885& 0.0884212\\ 0.0442106& -0.4863169\end{array}\right),$$

$$D_1=\left(\begin{array}{cc}1.90105& 0.0707369\\ 0.00265263& 0.0893054\end{array}\right),D_2=\left(\begin{array}{cc}2.0160099& 0.132632\\ 0.00530527& 0.344843\end{array}\right).$$

The total rate of type-1 and type-2 customers’ arrivals to the system is $\lambda =1.$ The coefficient of correlation of successive inter-arrival times in this arrival process is $0.256958,$ the squared coefficient of variation is $2.61753$. The average rate of type-1 customers’ arrivals ${\lambda }_1=0.377074,$ the average rate of type-2 customers’ arrivals is ${\lambda }_2=0.622926$.

The service process for type-1 customers $P{H}_1$ is defined by the vector ${\mathbf{b}}_1=(0.2,0.8)$ and the matrix $S_1=\left(\begin{array}{cc}-6.16045& 0.616045\\ 0.0123209& -1.23209\end{array}\right).$ The mean service time of type-1 customers is $0.7$, and the squared coefficient of variation of the service times is 1.24723.

The service process for type-2 customers $P{H}_2$ is defined by the vector ${\mathbf{b}}_2=(0.4,0.6)$ and the matrix $S_2=\left(\begin{array}{cc}-3.50388& 0.175194\\ 0.00350388& -0.700775\end{array}\right).$ The mean service time of type-2 customers is 1, and the squared coefficient of variation of the service times is 1.61091.

Let us fix the control parameter $K=7$ and suggest that the number of servers N is equal to 4 and the capacity R of Buffer-1 is 10. Let the rates of impatience of type-1 customers and type-2 customers be defined as ${\alpha }_r=0.01r,r=\overline{1,R},$ and ${ϵ}_i=0.03i,i\ge 0,$ respectively.

Let us vary the probability $q_1$ over the interval [0, 1] with step 0.05 and the probability $q_2$ over the interval $[q_1,1]$ also with step 0.05.

Figure 2 presents the dependence of the probability $P_1^{ent-loss}$ of an arbitrary type-1 customer loss upon arrival due to Buffer-1 overflow on the parameters $q_1$ and $q_2.$

One can see from Figure 2 that in the considered example, the loss probability $P_1^{ent-loss}$ decreases with growth in the probabilities $q_1$ and $q_2.$ It is evident that when the probabilities $q_1$ and $q_2$ grow, type-1 customers are chosen for service more often, and the probability of meeting Buffer-1 full upon type-1 customer arrival decreases. The minimal value of the loss probability $P_1^{ent-loss}$ is achieved for $q_1=q_2=1$ (what corresponds to the standard non-preemptive priority of type-1 customers) and is equal to 0.000018797.

Figure 3 illustrates the dependence of the probability $P_1^{imp-loss}$ of an arbitrary type-1 customer loss due to impatience on the parameters $q_1$ and $q_2.$

One can see from Figure 3 that in the considered example, the loss probability $P_1^{imp-loss}$ sharply decreases with growth in the probability $q_1$ and slightly decreases with growth in the probability $q_2.$ If the probability $q_1$ is small and the number of customers in Buffer-1 is less than $K=7$, type-1 customers are rarely chosen for service. Thus, the number of type-1 customers staying in Buffer-1 is bigger for small values of $q_1$, and more customers leave the buffer due to impatience.

The dependence of the total probability $P_1^{loss}$ of an arbitrary type-1 customer loss on the parameters $q_1$ and $q_2$ is presented in Figure 4.

The minimal value of the loss probability $P_1^{loss}$ is also reached for $q_1=q_2=1$ and is equal to 0.0011.

Figure 5 illustrates the dependence of the probability $P_2^{imp-loss}$ of an arbitrary type-2 customer loss due to impatience on the parameters $q_1$ and $q_2.$

This probability evidently increases with growth in the probabilities $q_1$ and $q_2$ since the growth in these probabilities leads to worse conditions for type-2 customer service. More customers of this type stay in the infinite buffer, and consequently, more customers leave it due to impatience.

The dependence of the total probability $P^{loss}$ of an arbitrary customer loss on the parameters $q_1$ and $q_2$ is presented in Figure 6.

In the considered case, the probability $P^{loss}$ decreases with growth in the probability $q_2$ and increases with growth in the probability $q_1.$ Note that for various input parameters, the dependence of the probability $P^{loss}$ can be absolutely different. This is because the increase in probabilities $q_l,l=1,2,$ implies lower probabilities of type-1 customers lost upon arrival and due to impatience but a higher probability of loss of type-2 customers. In the considered case, the minimal value of the loss probability $P^{loss}$ is reached for $q_1=0$ and $q_2=1$ and is equal to 0.003326. So, to have minimal customer loss, we should always choose type-2 customers, if any, for service if the number of type-1 customers in the buffer is less than $K=7$, and always choose type-1 customers otherwise. However, in real-world systems, the customers may have different importance to the system; that is, the charges for losing customers of different types may essentially differ.

Let us assume that the following economic criterion describes the system’s operating quality:

$$E=E(q_1,q_2)=c_1{\lambda }_1{P}_1^{loss}+c_2{\lambda }_2{P}_2^{imp-loss}$$

where $c_1$ is a fee paid by the system in the case when a type-1 customer is lost, and $c_2$ is a fee in the case when a type-2 customer from Buffer 2 is lost due to impatience.

The average losses of the system per unit of time are described by the economic criterion E. In order to optimize the system’s performance, we have to determine the probabilities $q_1$ and $q_2$ for which the economic criterion E assumes the lowest value.

Let us assume the following cost coefficients: $c_1=18,c_2=10.$

Figure 7 presents the dependence of the cost criterion E on the parameters $q_1$ and $q_2$.

The optimal value of the economic criterion is $E(q_1,q_2)=0.0446312$ and is achieved when $q_1=0.35$ and $q_2=1.$ Thus, if there are type-1 and type-2 customers in the buffers, it is reasonable to choose the type-1 customer for the next service with a probability of 0.35 if the number of type-1 customers in Buffer-1 is less than $K=7$, and always choose the type-1 customer otherwise.

Experiment 2. In the second experiment, we analyze the dependence of the system performance measures on the control parameters $q_1$ and $K.$ To this end, let us fix the probability $q_2=1$, i.e., when the number of customers in the Buffer-1 exceeds the threshold $K,$ then type-1 customers are always picked up for service. Additionally, let us increase the system load: we decrease the number of servers N to 3 and increase the capacity of Buffer-1 R from 10 to 20. The rest of the parameters are assumed to be the same as in the first experiment, including the arrival and service processes. We vary the control parameter K over the interval [2, 20] with step 1 and the probability $q_1$ over the interval [0, 1] with step 0.1.

Figure 8 illustrates the dependence of the probability $P_1^{ent-loss}$ of an arbitrary type-1 customer loss upon arrival due to Buffer-1 overflow on the parameters K and $q_1.$

One can see from Figure 8 that the probability $P_1^{ent-loss}$ decreases with an increase in the probability $q_1$ and increases with a growth in the parameter K. The maximal value of the probability $P_1^{ent-loss}$ is equal to $0.0018176$ and is reached for $q_1=0$ and $K=20$, i.e., when there are type-2 customers in the buffer, type-1 customers are chosen for service only if Buffer-1 is full. The minimal value of this probability is equal to $8.43\times {10}^{-6}$ and is reached for $q_1=q_2=1.$ Note that in the case of $q_1=q_2,$ the parameter K does not impact the system’s operation, and the system’s performance measures do not depend on K.

The dependence of the probability $P_1^{imp-loss}$ of an arbitrary type-1 customer loss due to impatience in Buffer-1 on the parameters K and $q_1$ is presented in Figure 9.

The probability $P_1^{imp-loss}$ also increases with the decrease of the probability $q_1$ and increases with growth in the parameter K. The same behavior we can also observe in Figure 10 that illustrates the dependence of the probability $P_1^{loss}$ of an arbitrary type-1 customer loss on the parameters K and $q_1.$

The loss probability $P_1^{loss}$ takes its minimal value 0.002738 for $q_1=1.$ This loss mainly occurs due to type-1 customers’ impatience.

The probability $P_2^{imp-loss}$ of an arbitrary type-2 customer loss due to impatience in Buffer-2 behaves oppositely. This probability grows with an increase in the probability $q_1$ and a decrease in the parameter K. The shape of the dependence of the probability $P_2^{imp-loss}$ on the parameters K and $q_1$ is presented in Figure 11.

The dependence of the probability $P_{loss}$ of an arbitrary customer loss on the values K and $q_1$ is presented in Figure 12.

The minimal value of the probability $P_{loss}$ is equal to 0.009641 for $K=16$ and $q_1=0.$

The quality of system operation can be described by the economic criterion $E=E(q_1,K)$, which is given in the same way as the criterion from the first experiment. The coefficients $c_1$ and $c_2$ also coincide with the cost coefficients from the previous experiment. The dependence of the cost criterion $E=E(q_1,K)$ on the parameters $q_1$ and K is presented in Figure 13.

The economic criterion takes the optimal value $0.134653$ for $K=11$ and $q_1=0.$

Experiment 3. In the third experiment, instead of a correlated $MMAP$ arrival process, we consider two stationary Poisson arrival processes and will show that accounting for only the mean arrival rates and ignoring the correlation and variation of the inter-arrival times can lead to significant errors in assessing the performance of the system under consideration.

Consider the $MMAP$ arrival process characterized by the value $W=1$ with the following matrices: $D_0=(-1),D_1=\left(0.377074\right),D_2=\left(0.622926\right).$ It defines a heterogeneous stationary Poisson arrival process with two types of customers having arrival rates ${\lambda }_1=0.377074$ and ${\lambda }_2=0.622926.$ The other parameters coincide with those presented in the second experiment. The control parameters $q_1$ and K are varied in the same way as in the previous experiment.

Figure 14, Figure 15 and Figure 16 show the dependence of the probabilities $P_1^{loss}$, $P_2^{imp-loss}$ and the cost criterion E on the parameters K and $q_1$ in the case of the heterogeneous stationary Poisson arrival flow.

One can draw the conclusion that the correlation and variance of the arrival process have a significant impact on the systems’ performance indicators by comparing Figure 10 and Figure 14, Figure 11 and Figure 15, and Figure 13 and Figure 16, respectively. If someone tries to evaluate the system with correlated inter-arrival times using the model with the heterogeneous stationary Poisson arrival flow, he/she can obtain huge errors in the estimation. In the case of the heterogeneous stationary Poisson arrival flow, the minimal value of the probability $P_1^{loss}$ is equal to 0.000262, while in the case of the $MMAP$ the minimal value of this probability is 0.002738. The minimal value of the probability $P_2^{imp-loss}$ is equal to 0.00087 for heterogeneous stationary Poisson arrival flow, while in the case of the $MMAP$ the minimal value of this probability is 0.00717. That means that the real values of the main loss probabilities can be 10 times higher than the values expected based on the model with heterogeneous stationary Poisson arrival flow.

The minimal value of the cost criterion $E^{*}$ in the case of the heterogeneous stationary Poisson arrival flow is equal to 0.00833988, which is 16 times less than this value for the case of the $MMAP$ arrival process.

Therefore, ignorance of positive correlation leads to an overly optimistic prediction of the performance characteristics of the system and an underestimation of the requirements for providing the desired quality of service, such as service rates or the number of servers. This, along with the observation that flows in many real-world systems exhibit correlation and high variability in inter-arrival times, clearly motivates the necessity of analyzing systems with the $MMAP$.

6. Conclusions

We analyzed the flexible randomized threshold strategy of priority access in a multi-server queueing system with the $MMAP$ arrival process and the phase-type distribution of the service times dependent on the type of customer. The impatience of customers of both types during their stay in the buffers is taken into account. Analysis of the multi-dimensional Markov chain describing the dynamics of the system is implemented via the use of the generalized $PH$ distribution and the results known for asymptotically quasi-Toeplitz Markov chains. Numerical results are presented. They give some insight into the influence of the parameters of the control strategy. The necessity of accounting for correlation in the arrival process is shown.

Analysis can be extended to the cases of several thresholds and to the systems in which service requires the presence of some additional inventory, e.g., energy harvested online.