Analysis of Instantaneous Feedback Queue with Heterogeneous Servers

Abstract

A system with heterogeneous servers, Markov Modulated Poisson flow and instantaneous feedback is studied. The primary call is serviced on a high-speed server, and after it is serviced, each call, according to the Bernoulli scheme, either leaves the system or requires re-servicing. After the completion of servicing of a call in a slow server, according to the Bernoulli scheme, it also either leaves the system or requires re-servicing. If upon arrival of a primary call the queue length of such calls exceeds a certain threshold value and the slow server is free, then the incoming primary call, according to the Bernoulli scheme, is either sent to the slow server or joins its own queue. A mathematical model of the studied system is constructed in the form of a three-dimensional Markov chain. Approximate algorithms for calculating the steady-state probabilities of the models with finite and infinite queues are proposed and their high accuracy is shown. The results of numerical experiments are presented.

1. Introduction

When building mathematical models of communication processes, as well as processes related to the processing of parts in production systems, inventory management processes, etc., it becomes necessary to take into account the repeated processing of some calls. For example, in many information transmission systems, erroneously transmitted data (packets, frames, etc.) are re-transmitted, since the performance of such systems is often evaluated by the reliability of data transmission. Similarly, in manufacturing systems, if the parts produced have certain defects, then in some situations they need to be re-processed.

Taking into account that the effect of repeated call servicing leads to the need to use models of systems with feedback (Feed Back Queue, FBQ), such models were first introduced in [1,2]. It is necessary to distinguish between Instantaneous Feed Back Queue (IFBQ) systems, where a retry occurs immediately after the service is complete, and Delayed Feed Back Queue (DFBQ) systems, in which a replay of a service request occurs after a certain positive time.

The current state of the problem of studying FBQ models is described in detail in [3], so here, we will not dwell on the presentation of the results known in this direction. In this work, it is noted that the overwhelming majority of works study FBQ models without a buffer for waiting for calls. At the same time, in many systems, buffers are organized for waiting of calls in queue. Therefore, to adequately describe the operation of such systems, it becomes necessary to study FBQ models with buffers.

Based on the above facts, in this paper we study the IFBQ model with buffers. Note that there are few works in the available literature devoted to the study of such models. Among them, we note the works [4,5,6,7,8,9,10,11,12,13]. Simple one-dimensional IFBQ models with one server and impatient calls using various mechanisms for keeping them in the queue at the moments of the end of the admissible waiting time are studied in [4,5,6,7,8,9,10]. A model with two heterogeneous servers and a bounded queue, in which incoming calls with known probabilities are assigned to the servers, was studied in [11]. Note that in [4,5,6,7,8,9,10,11] the primary calls and the calls that require repeated servicing are not distinguished. Therefore, using the approach proposed in them, it is impossible to find the distribution of the number of calls that require repeated servicing. A more complex IFBQ model with one server and bounded queue was studied in [12], where the incoming stream is a Markov Arrival Process (MAP) stream and the service time of calls has a phase-type distribution function. It is shown that the mathematical model of the system is a certain four-dimensional Markov Chain (MC) and the performance measures of the system are calculated. In the recent work [13], the IFBQ model with one server and two Poisson streams was studied. It is considered that for each type of call there are separate buffers of infinite size, while only high priority calls can repeat requests for servicing, according to the Bernoulli scheme. It is assumed that calls that require re-servicing are sent to the queue of low priority calls, while models with preemptive and non-preemptive priorities are studied. In [12,13], the matrix-geometric method [14] is used to study the proposed models.

An analysis of the available works showed that they are based on a number of assumptions. The main ones are the following: (1) initial and repeated applications are identical in all respects; (2) calls can only repeat a request once for service; (3) calls of both types (primary and secondary or feedback calls) are transmitted by a single server.

However, in many real systems, these assumptions are not met, and therefore, in order to improve the adequacy of the models, this work studies an IFBQ in which the above assumptions are not met. In this case, it is assumed that the servers are specialized in the type of requests, and they are heterogeneous, i.e., have different service rates (earlier systems with heterogeneous servers in the absence of feedback were studied in many works, see, for example, [15] and its bibliography). In addition, in contrast to the standard assumption that the incoming stream is Poisson, here we study a model with a Markov Modulated Poisson Process (MMPP) stream [16]. Note that the proposed model is studied using the method of hierarchical space merging of a multi-dimensional MC [17]. Application of this method makes it possible to develop effective numerical procedures for calculating the performance measures of the system under study.

The paper is structured as follows. Section 2 contains a description of the system under study and a statement of the problem. In Section 3, a mathematical model of the system in the form of a 3D MC is developed and its generating matrix is determined. Here, the ergodicity condition of the model is obtained and an approximate method for calculating the steady-state probabilities and performance measures of the system is developed. Section 4 presents the results of numerical experiments. Concluding remarks are given in Section 5.

2. Description of the System with Instantaneous Feedback and Formulation of the Problem

As in [17], it is assumed that input flow is MMPP with parameters $\left(\mathsf{\Sigma },\hspace{0.17em}\mathsf{\Lambda }\right)$. This means that the MC, which controls the intensities of the incoming flow, has a generating matrix (GM) $\mathsf{\Sigma }=\Vert {\sigma }_{ij}\Vert$ of dimension $N\times N$, where ${\sigma }_{ij}$ determines the transition intensity from state $i$ to state $j$, $i, j=1,\hspace{0.17em}\hspace{0.17em}2,\hspace{0.17em}\hspace{0.17em}\dots ,\hspace{0.17em}\hspace{0.17em}N$, and ${\sigma }_{ii}=-{\displaystyle \sum _{j=1,\hspace{0.17em}j\ne i}^N{\sigma }_{ij}}$. It is assumed that when MC is in a state $n$, the intensity of the incoming stream is equal to ${\lambda }_n,\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N$, and when the state of the MC changes, the intensity of the incoming flow also changes instantly. The vector of intensities of the incoming flow is denoted by $\mathsf{\Lambda }=\left({\lambda }_1,\hspace{0.17em}{\lambda }_2,\hspace{0.17em}\dots \hspace{0.17em},\hspace{0.17em}{\lambda }_N\right)$.

The system contains two servers, a high-speed (F-server) and a low-speed (S-server) server, while primary calls (p-calls), as a rule, are served in the F-server; after the completion of servicing, the p-calls, according to the Bernoulli scheme or with probability $\alpha$, leave the system, or, with probability $1-\alpha$, are sent to be served in the S-server. Calls that require re-servicing (secondary calls, s-calls) can queue up in front of the S-server. It is assumed that s-calls can repeatedly require re-servicing, i.e., after the completion of service, each s-call is independent of other calls according to the Bernoulli scheme or with the probability $\beta$ leaves the system permanently or with a complementary probability $1-\beta$ instantly requires re-servicing in the S-server. The service times of calls in the F-server and in the S-server are random variables that have an exponential distribution with parameters ${\mu }_f$ and ${\mu }_s$, respectively, where ${\mu }_s<{\mu }_f$.

If at the moment of receipt of a p-call the number of calls in the queue in front of the F-server is higher than a certain threshold value, $M>0$, and at the same time the S-server is not busy, then the received p-call either with probability $\sigma$ is sent for service to the S-server or with a complementary probability $1-\sigma$ it joins its queue.

Serving different types of calls and changing the states of the MC, which controls the intensity of the incoming flow, are random processes independent of each other.

The problem is to find the joint distribution of the states of the MMPP flow and the number of calls of each type in the system. The solution to this problem will allow us to find the performance measures of the system, i.e., the average number of p-calls $\left(L_p\right)$ and s-calls $\left(L_s\right)$ in the system and the intensity of p-calls served in the S-server $\left(R_{ps}\right)$.

3. Calculation of the Steady-State Probabilities and Performance Measures

The state of the system is determined by the 3D vector $(n,\hspace{0.17em}k,\hspace{0.17em}r)$, where $n$ is the state of the MC that controls the intensity of the incoming flow, $k$ is the total number of calls in front of the F-server and calls in it, and $r$ is the total number of calls in front of the S-server and calls in it. Then, the state space of this 3D MC is defined as the Cartesian product of three sets:

$$E=\left\{1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N\right\}\times \left\{0,\hspace{0.17em}1,\hspace{0.17em}\dots ,\hspace{0.17em}\right\}\times \left\{0,\hspace{0.17em}1,\hspace{0.17em}\dots ,\hspace{0.17em}\right\}.$$

(1)

The intensity of transition from state $(n,\hspace{0.17em}k,\hspace{0.17em}r)$ to state $(n^{\prime },\hspace{0.17em}{k}^{\prime },\hspace{0.17em}{r}^{\prime })$ will be denoted by $q((n,k,r),(n^{\prime },k^{\prime },r^{\prime }))$. These values are calculated as follows:

• transitions $(n,\hspace{0.17em}k,\hspace{0.17em}r)\to \left(n^{\prime },k,r\right),\hspace{0.17em}{n}^{\prime }\ne n$ are carried out with intensity ${\sigma }_{n,n^{\prime }}$ when the state of the MC that controls the intensity of the incoming MMPP flow changes;
• transitions $(n,\hspace{0.17em}k,\hspace{0.17em}r)\to \left(n,k+1,r\right),\hspace{0.17em}r>0,$ and $(n,\hspace{0.17em}k,\hspace{0.17em}0)\to \left(n,k+1,0\right),\hspace{0.17em}k<M,$ are carried out with intensity ${\lambda }_n$ when a p-call arrives;
• transitions $\left(n,\hspace{0.17em}k,\hspace{0.17em}0\right)\to \left(n,k+1,\hspace{0.17em}0\right),\hspace{0.17em}k\ge M,$ are carried out with intensity ${\lambda }_n\left(1-\sigma \right)$ when a p-call arrives;
• transitions $\left(n,\hspace{0.17em}k,\hspace{0.17em}0\right)\to \left(n,k+1,\hspace{0.17em}1\right),\hspace{0.17em}k\ge M,$ are carried out with intensity ${\lambda }_n\sigma$ when a p-call arrives;
• transitions $(n,\hspace{0.17em}k,\hspace{0.17em}r)\to \left(n,k-1,r\right),\hspace{0.17em}k>0,\hspace{0.17em}$ are carried out with intensity ${\mu }_f\alpha$ upon completion of servicing the p-call in the F-server and leaving the system;
• transitions $(n,\hspace{0.17em}k,\hspace{0.17em}r)\to \left(n,k-1,r+1\right),\hspace{0.17em}k>0,$ are carried out with intensity ${\mu }_f\left(1-\alpha \right)$ upon completion of servicing the p-call in the F-server and returning to the S-server for re-servicing;
• transitions $(n,\hspace{0.17em}0,r)\to \left(n,0,r-1\right),\hspace{0.17em}r>0,\hspace{0.17em}$ are carried out with intensity ${\mu }_s\beta$ upon completion of servicing the p-call in the S-server.

Therefore, the positive elements of the GM of the given 3D MC are determined from the following relations:

$$q\left(\left(n,k,r\right),\left(n^{\prime },k^{\prime },r^{\prime }\right)\right)=\{\begin{array}{ll}{\sigma }_{n{n}^{\prime }},\hfill & \mathrm{if} n^{\prime }\ne n,\hspace{0.17em}{k}^{\prime }=k,\hspace{0.17em}{r}^{\prime }=r,\hfill \\ {\lambda }_n,\hfill & \begin{array}{l}\mathrm{if} r>0,\hspace{0.17em}{n}^{\prime }=n,\hspace{0.17em}{k}^{\prime }=k+1,\hspace{0.17em}{r}^{\prime }=r\\ \mathrm{or} k<M,\hspace{0.17em}r=0,n^{\prime }=n,\hspace{0.17em}{k}^{\prime }=k+1,\hspace{0.17em}{r}^{\prime }=r,\end{array}\hfill \\ {\lambda }_n\sigma ,\hfill & \mathrm{if} r=0,\hspace{0.17em}{n}^{\prime }=n,\hspace{0.17em}\hspace{0.17em}{k}^{\prime }=k,\hspace{0.17em}\hspace{0.17em}{r}^{\prime }=1,\hfill \\ {\lambda }_n\left(1-\sigma \right),\hfill & \mathrm{if} k\ge M,\hspace{0.17em}r=0,n^{\prime }=n,\hspace{0.17em}{k}^{\prime }=k+1,\hspace{0.17em}{r}^{\prime }=r,\hfill \\ {\mu }_f\alpha ,\hfill & \mathrm{if} k>0,\hspace{0.17em}{n}^{\prime }=n,\hspace{0.17em}{k}^{\prime }=k-1,\hspace{0.17em}{r}^{\prime }=r,\hfill \\ {\mu }_f\left(1-\alpha \right),\hfill & \mathrm{if} k>0,\hspace{0.17em}{n}^{\prime }=n,\hspace{0.17em}{k}^{\prime }=k-1,\hspace{0.17em}{r}^{\prime }=r+1,\hfill \\ {\mu }_s\beta ,\hfill & \mathrm{if} k=0,\hspace{0.17em}r>0,\hspace{0.17em}{n}^{\prime }=n,\hspace{0.17em}{k}^{\prime }=0,\hspace{0.17em}{r}^{\prime }=r-1.\hfill \end{array}$$

(2)

Let us denote by $p(n,k,r)$ the probability of state $(n,k,r)\in E$. The condition for the existence of a stationary regime is obtained below.

Finding the steady-state probabilities is sufficient to calculate the performance measures of the system under study. Thus, the average number of p-calls $\left(L_p\right)$ and s-calls $\left(L_s\right)$ in the system are defined as the mathematical expectation of the corresponding random variables:

$$L_p={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^{\infty }k{\displaystyle \sum _{r=0}^{\infty }p\left(n,k,r\right)}}}\hspace{0.17em};$$

(3)

$$L_s={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{r=1}^{\infty }r{\displaystyle \sum _{k=0}^{\infty }p\left(n,k,r\right)}}}\hspace{0.17em}.$$

(4)

The intensity of p-calls, which are served in the S-server $\left(R_{ps}\right)$, is defined as:

$$R_{ps}=\sigma {\displaystyle \sum _{n=1}^N{\lambda }_n}{\displaystyle \sum _{k=M}^{\infty }p\left(n,\hspace{0.17em}k,\hspace{0.17em}0\right)}\hspace{0.17em}.$$

(5)

Using the method of multidimensional generating functions to find the steady-state probabilities faces a number of methodological and computational difficulties. On this basis, below is described an alternative method for solving this problem, based on the hierarchical space merging method of a multidimensional MC [17].

By taking into account that changing the states of the MC, which controls the intensity of the incoming flow, are independent on other parameters of the system, consider the following splitting of the state space (1):

$$E={\displaystyle \underset{n=1}{\overset{N}{\cup }}{E}_n}\hspace{0.17em},\hspace{0.17em}{E}_n{\displaystyle \cap E_{n^{\prime }}=\varnothing \hspace{0.17em}, \mathrm{if} n\ne n^{\prime },}$$

(6)

where $E_n=\left\{\left(n,k,r\right)\in E:\hspace{0.17em}\hspace{0.17em}k=0,\hspace{0.17em}1,\hspace{0.17em}\dots ;\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots \right\},\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N.$

All states from the class $E_n$ are combined into one merged state, $<n>$, and on the basis of splitting (6) in the state space (1), the merge function $U_1\left(n,k,r\right)=<n>\hspace{0.17em},\hspace{0.17em}\left(n,k,r\right)\in E_n\hspace{0.17em}$ is determined. Let us denote the set of merged states $<n>$ by ${\mathsf{\Omega }}_1=\left\{<n>\hspace{0.17em}:\hspace{0.17em}\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N\right\}\hspace{0.17em}.$

Then, the approximate values of the steady-state probabilities of the original model, denoted by $\tilde{p}\left(n,k,r\right),$ are defined as (see [17])

$$\tilde{p}\left(n,k,r\right)={\rho }_n\left(k,r\right){\pi }_1\left(<n>\right)\hspace{0.17em},$$

(7)

where ${\rho }_n\left(k,r\right)$ is the probability of the state $(k,r)$ inside a split model with state space $E_n$, and ${\pi }_1\left(<n>\right)$ is the probability of the merged state $<n>\in \mathsf{\Omega }{\hspace{0.17em}}_1.$

As noted above, the transitions between the states of the MC that control the intensity of the incoming flow do not depend on the statuses of the F-server and S-server, and therefore, the probabilities of the states ${\pi }_1\left(<n>\right),\hspace{0.17em}<n>\in {\mathsf{\Omega }}_1\hspace{0.17em},$ are defined by its GM $\mathsf{\Sigma }$. Consequently, to find the steady-state probabilities of the original model, we only need to determine the stationary distributions of a 2D MC with state spaces $E_n,\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N$ (see Formula (7)). To solve this problem, a merge procedure (second level of the hierarchy) is applied to this 2D MC. Since all split models with state spaces $E_n$ have identical GM, then we fix the value of the parameter $n,\hspace{0.17em}\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N$.

Below we assume that $\alpha >>1-\alpha$ and $\sigma <<1-\sigma$, i.e., $\alpha >>0.5$ and $\sigma <<0.5$, and consider the following partition in the class $E_n$:

$$E_n={\displaystyle \underset{r=0}{\overset{\infty }{\cup }}{E}_n^r}\hspace{0.17em},\hspace{0.17em}{E}_n^r{\displaystyle \cap E_n^{r^{\prime }}=\varnothing \hspace{0.17em},\mathrm{if} r\ne r^{\prime },}$$

(8)

where $E_n^r=\left\{\left(k,r\right)\in E_n:\hspace{0.17em}\hspace{0.17em}k=0,\hspace{0.17em}1,\hspace{0.17em}\dots \right\},\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots$.

In accordance with the above accepted assumptions, transition intensities between states within each split class $E_n^r$ are essentially larger than transition intensities between states from different splitting classes. On the basis of splitting (8), all states from class $E_n^r$ are combined into one merged state $<r>$, and in the state space $E_n$ the merge function $U_2\left(k,r\right)=<r>\hspace{0.17em},\hspace{0.17em}\left(k,r\right)\in E_n^r\hspace{0.17em}$ is determined. Let us denote the set of merged states $<r>$ by ${\mathsf{\Omega }}_2=\left\{<r>\hspace{0.17em}:\hspace{0.17em}\hspace{0.17em}n=0,\hspace{0.17em}1,\hspace{0.17em}\hspace{0.17em}2,\hspace{0.17em}\dots \right\}\hspace{0.17em}.$

According to [17], we have:

$${\rho }_n\left(k,r\right)={\rho }_n^r\left(k\right){\pi }_2^n\left(<r>\right)\hspace{0.17em},$$

(9)

where ${\rho }_n^r\left(k\right)$ is the probability of the $(k,r)$ inside a split model with state space $E_n^r$, and ${\pi }_2^n\left(<r>\right)$ is the probability of the merged state $<r> \in {\mathsf{\Omega }}_2\hspace{0.17em}.$

In class $E_n^r,\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots ,$ in the all state vectors, the second component is to $r$. Therefore, when studying the split models with space $E_n^r$, every state $(k,r)$ can only be specified by the first component, i.e., in further the states $(k,r)\in E_n^r$ are denoted as $k,k=0,\hspace{0.17em}1,\hspace{0.17em}\dots$. Then, from relations (2) we conclude that the intensities of transitions between the states of the split model with the state space $E_n^r,\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots$ depend on the parameter $r,\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots$ and are defined as follows:

Case $r=0\hspace{0.17em}:$

$$q_n^0\left(k,k^{\prime }\right)=\{\begin{array}{ll}{\lambda }_n,\hfill & \mathrm{if} k<M,\hspace{0.17em}{k}^{\prime }=k+1,\hfill \\ {\lambda }_n\left(1-\sigma \right),\hfill & \mathrm{if} k\ge M,\hspace{0.17em}{k}^{\prime }=k+1,\hfill \\ {\mu }_f\alpha ,\hfill & \mathrm{if} k^{\prime }=k-1.\hfill \end{array}$$

(10)

Cases $r>0\hspace{0.17em}:$

$$q_n^r\left(k,k^{\prime }\right)=\{\begin{array}{ll}{\lambda }_n,\hfill & \mathrm{if} k^{\prime }=k+1,\hfill \\ {\mu }_f\alpha ,\hfill & \mathrm{if} k^{\prime }=k-1.\hfill \end{array}$$

(11)

From (10) we obtain that under the fulfilling of the condition ${\nu }_n<{\left(1-\sigma \right)}^{-1},\hspace{0.17em}{\nu }_n={\lambda }_n/{\mu }_f\alpha \hspace{0.17em},$ state probabilities of a split model with state space $E_n^0$ are defined as follows:

$${\rho }_n^0\left(k\right)=\{\begin{array}{ll}{\nu }_n^k{\rho }_n^0\left(0\right),\hfill & \mathrm{if} 0\le k\le M,\hfill \\ \frac{1}{{\left(1-\sigma \right)}^M}{\left({\nu }_n\left(1-\sigma \right)\right)}^k{\rho }_n^0\left(0\right),\hfill & \mathrm{if} k>M\hspace{0.17em},\hfill \end{array}$$

(12)

where ${\rho }_n^0\left(0\right)$ is determined from the normalization condition, i.e., ${\displaystyle \sum _{k=0}^{\infty }{\rho }_n^0\left(k\right)}=1\hspace{0.17em}.$ After standard transformations we get

$${\rho }_n^0\left(0\right)={\left(\frac{1-{\nu }_n^{M+1}}{1-{\nu }_n}+{\nu }_n^{M+1}\frac{1-\sigma }{1-{\nu }_n\left(1-\sigma \right)}\right)}^{-1}.$$

(13)

Remark 1.

Since the condition${\nu }_n<{\left(1-\sigma \right)}^{-1}$is satisfied for every$n,\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N$, then we get that the following condition must be satisfied:

$$\underset{n=\overline{1,N}}{\max }\left\{{\nu }_n\right\}<{\left(1-\sigma \right)}^{-1}\hspace{0.17em}.$$

(14)

From (11) we obtain that under the fulfilling of the condition${\nu }_n<1$, state probabilities of split models with state spaces$E_n^r$are not dependent on index$r,\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots ,$and are determined as follows:

$${\rho }_n^r\left(k\right)=\left(1-{\nu }_n\right){\nu }_n^k\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=0,\hspace{0.17em}1,\hspace{0.17em}\hspace{0.17em}\dots ,\hspace{0.17em}$$

(15)

Remark 2.

Since the condition${\nu }_n<1$should be satisfied for each$n,\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N$, then we get that the following condition must be satisfied:

$$\underset{n=\overline{1,N}}{\max }\left\{{\nu }_n\right\}<1.$$

(16)

Then, combining conditions (14) and (16), we obtain the first ergodicity condition for the model, i.e., condition (16) must be satisfied.

Taking into account (2), (12), (13), and (15), we obtain that the intensities of transitions between states${\mathsf{\Omega }}_2$are defined as

$$q_n\left(<r>,<r^{\prime }>\right)=\{\begin{array}{ll}{\eta }_n\left(0,\hspace{0.17em}1\right),\hfill & \mathrm{if} r=0,\hspace{0.17em}{r}^{\prime }=1,\hfill \\ {\mu }_f\left(1-\alpha \right){\nu }_n,\hfill & \mathrm{if} r>0{,}^{\prime }{r}^{\prime }=r+1,\hfill \\ {\mu }_s\beta ,\hfill & \mathrm{if} r>0,\hspace{0.17em}{r}^{\prime }=r-1.\hfill \end{array}$$

(17)

Taking into account (12) and (13), the value${\eta }_n\left(0,\hspace{0.17em}1\right)$in Formula (17) is defined as follows:

$$\begin{array}{l}{\eta }_n\left(0,\hspace{0.17em}1\right)={\mu }_f\left(1-\alpha \right){\displaystyle \sum _{k=1}^{\infty }{\rho }_0^0\left(k\right)}+{\lambda }_n\sigma {\displaystyle \sum _{k=M}^{\infty }{\rho }_n^0\left(k\right)}={\mu }_f\left(1-\alpha \right){\nu }_n+\\ +{\lambda }_n\sigma {\rho }_n^0\left(0\right)\frac{1}{{\left(1-\sigma \right)}^M}{\displaystyle \sum _{k=M}^{\infty }{\left({\nu }_n\left(1-\sigma \right)\right)}^k}={\mu }_p\left(1-\alpha \right){\nu }_n+{\lambda }_n\sigma {\rho }_n^0\left(0\right)\hspace{0.17em}\frac{{\nu }_n^M}{1-{\nu }_n\left(1-\sigma \right)}\hspace{0.17em}.\\ \end{array}$$

From (17) we obtain that under the condition${\psi }_n<1,\hspace{0.17em}{\psi }_n={\mu }_f\left(1-\alpha \right){\nu }_n/{\mu }_s\beta ,$state probabilities of a merged model with a state space${\mathsf{\Omega }}_2$are defined as follows:

$${\pi }_2^n\left(<r>\right)=\{\begin{array}{ll}{\theta }_n{\pi }_2^n\left(<0>\right),\hfill & \mathrm{if} r=1,\hfill \\ {\theta }_n{\psi }_n^{r-1}{\pi }_2^n\left(<0>\right),\hfill & \mathrm{if} r>1\hspace{0.17em},\hfill \end{array}$$

(18)

where${\theta }_n=\frac{{\eta }_n\left(0,\hspace{0.17em}1\right)}{{\mu }_s\beta }$and${\pi }_2^n\left(<0>\right)$is determined from the normalization condition, i.e.,${\displaystyle \sum _{r=0}^{\infty }{\pi }_2^n\left(<r>\right)}=1\hspace{0.17em}.$After standard transformations we get

$${\pi }_2^n\left(<0>\right)=\frac{1-{\psi }_n}{1-{\psi }_n+{\theta }_n}.$$

(19)

Remark 3.

Since the condition${\psi }_n<1$should be satisfied for each$n,\hspace{0.17em}n=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N$, then we obtain the second ergodicity condition of the model:

$$\underset{n=\overline{1,N}}{\max }\left\{{\nu }_n\right\}<\frac{{\mu }_s\beta }{{\mu }_f\left(1-\alpha \right)}\hspace{0.17em}.$$

(20)

Combining relations (16) and (20), we obtain the following condition for the ergodicity of the model:

$$\underset{n=\overline{1,N}}{\max }\left\{{\nu }_n\right\}<\min \left\{\frac{{\mu }_s\beta }{{\mu }_f\left(1-\alpha \right)}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}1\right\}.$$

(21)

Remark 4.

From (21) we conclude that the ergodic condition does not depend on the parameters of the GM of the MC that controls the intensity of the MMPP flow.

Thus, when the ergodicity condition of model (21) is satisfied, taking into account relations (7), (9), (12), (13), (18), and (19), the approximate values of the steady-state probabilities of the original 3D MC are found.

Further, using steady-state probabilities, approximate values of performance measures (3)–(5) can be calculated. Indeed, from (7) and (9) we conclude that these performance measures are calculated as follows:

$$L_p\approx {\displaystyle \sum _{n=1}^N{\pi }_1\left(<n>\right)}{\displaystyle \sum _{k=1}^{\infty }k{\displaystyle \sum _{r=0}^{\infty }{\rho }_n^r\left(k\right){\pi }_2^n\left(<r>\right)}}\hspace{0.17em};$$

(22)

$$L_s\approx {\displaystyle \sum _{n=1}^N{\pi }_1\left(<n>\right){\displaystyle \sum _{r=1}^{\infty }r{\pi }_2^n\left(<r>\right)}}\hspace{0.17em};$$

(23)

$$R_{ps}\approx \sigma {\displaystyle \sum _{n=1}^N{\pi }_1\left(<n>\right){\displaystyle \sum _{k=M}^{\infty }{\rho }_n^0\left(k\right){\pi }_2^n\left(<0>\right)}}\hspace{0.17em}.$$

(24)

Then, after certain mathematical calculations from (22)–(24) we get

$$L_p\approx {\displaystyle \sum _{n=1}^N{\pi }_1\left(<n>\right)}\left({\mathsf{\Psi }}_n^1+{\mathsf{\Psi }}_n^2+{\mathsf{\Psi }}_n^3\right)\hspace{0.17em},$$

(25)

where

$$\begin{array}{l}{\mathsf{\Psi }}_n^1={\displaystyle \sum _{k=1}^{M}k{\nu }_n^k\left({\rho }_n^0\left(0\right)+\left(1-{\nu }_n\right)\left(1-{\pi }_2^n\left(<0>\right)\right)\right)}\hspace{0.17em},\\ \hspace{0.17em}{\mathsf{\Psi }}_n^2={\rho }_0^0\left(0\right){\pi }_2^n\left(<0>\right)\frac{{\nu }_n^2}{{\left(1-\sigma \right)}^{M-2}}\frac{M\left(1-{\nu }_n\left(1-\sigma \right)\right)+1}{{\left(1-{\nu }_n\left(1-\sigma \right)\right)}^2}\hspace{0.17em}\hspace{0.17em},\\ {\mathsf{\Psi }}_n^3=\left(1-{\pi }_2^n\left(<0>\right)\right)\frac{{\nu }_n^2}{1-{\nu }_n}\frac{M\left(1-{\nu }_n\right)+1}{{\left(1-\sigma \right)}^M}\hspace{0.17em};\end{array}$$

$$L_s\approx {\displaystyle \sum _{n=1}^N{\pi }_1\left(<n>\right){\pi }_2^n\left(<0>\right)\frac{{\theta }_n}{{\left(1-{\psi }_n\right)}^2}}\hspace{0.17em};$$

(26)

$$R_{ps}\approx \sigma {\displaystyle \sum _{n=1}^N\frac{{\nu }_n^M}{1-{\nu }_n\left(1-\sigma \right)}{\pi }_1\left(<n>\right){\pi }_2^n\left(<0>\right){\rho }_n^0\left(0\right)}\hspace{0.17em}.$$

(27)

It is important to note that the proposed approach can be used to analyze a similar model with limited buffers in front of servers of different types. Indeed, let the sizes of buffers in front of flow and slow servers be equal to $R_f<\infty$ and $R_s<\infty$, respectively.

The state space of this model is finite-dimensional and is defined as $E_1=\left\{1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}N\right\}\times \left\{0,\hspace{0.17em}1,\hspace{0.17em}\dots ,\hspace{0.17em}{R}_f\right\}\times \left\{0,\hspace{0.17em}1,\hspace{0.17em}\dots ,\hspace{0.17em}{R}_s\right\}\hspace{0.17em}.$ Elements of the GM of the corresponding 3D MC are defined similarly to relation (2) by taking into account that $1\le M\le R_f$.

Note that for any positive values of the initial parameters of the system, a stationary mode exists in this Markov chain, since it is irreducible and has a finite state space. In this model, along with performance measures (3)–(5), a new measure appears, i.e., the probability of losing p-calls, $P_p$. It is defined as follows:

$$P_p={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{r=1}^{R_s}p\left(n,R_f,r\right)}}\hspace{0.17em}.$$

(28)

To find the steady-state probabilities of the appropriate 3D MC, an exact method based on the use of a system of equilibrium equations can be used here. This approach is effective only for models of moderate dimensions. For large-scale models, the approximate approach described above can be used.

Let’s briefly consider the applications of the approximate approach. Here, splitting (6) is also considered at the first level of the hierarchy, and splitting (8) is used at the second level of the hierarchy, but it is taken into account that the infinity sign is replaced by the parameter $R_s\hspace{0.17em}.$ Probabilities of states inside split classes $E_n^r,\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots ,\hspace{0.17em}{R}_s$ are calculated as follows:

Case $r=0\hspace{0.17em}:$

$${\rho }_n^0\left(k\right)=\{\begin{array}{ll}{\nu }_n^k{\rho }_n^0\left(0\right),& \mathrm{if}\hspace{0.17em}0\le k\le M,\\ \frac{1}{{(1-\sigma )}^M}{\left({\nu }_n(1-\sigma )\right)}^k{\rho }_n^0\left(0\right),& \mathrm{if}\hspace{0.17em}M<k\le R_f\end{array}\hspace{0.17em}$$

(29)

where ${\rho }_n^0\left(0\right)={\left(\frac{1-{\nu }_n^{M+1}}{1-{\nu }_n}+{\nu }_n^{M+1}\frac{1-a^{R_f-M}}{1-a}\right)}^{-1},\hspace{0.17em}a={\nu }_n\left(1-\sigma \right)$;

Cases $r>0\hspace{0.17em}:$

$${\rho }_n^r\left(k\right)={\nu }_n^k\frac{1-a^{R_f+1}}{1-a}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=0,\hspace{0.17em}1,\hspace{0.17em}\hspace{0.17em}\dots ,\hspace{0.17em}{R}_f\hspace{0.17em}.$$

(30)

The probabilities of merged states in this model are determined similarly to (18). The normalizing constant, ${\pi }_2^n\left(<0>\right)$, in this model is determined from the normalization condition, i.e., ${\displaystyle \sum _{r=0}^{R_s}{\pi }_2^n\left(<r>\right)}=1\hspace{0.17em}.$ Further, as above, using Formula (29)–(31), the steady-state probabilities are determined and the approximate values of the performance measures of the system (3)–(5) and (28) are calculated.

4. Numerical Results

For the sake of brevity, we consider here the results of numerical experiments only for the model with infinite buffers. Numerical experiments carried out here have two purposes: (1) to show the high accuracy of the developed algorithm for the approximate calculation of steady-state probabilities of the system under study; (2) to study the behavior of the performance measures of the system with respect to changes in the threshold parameter $M$.

The accuracy of the approximate algorithm is estimated using simulation. In this case, the closeness of the results obtained using different approaches is estimated using the cosine similarity norm, i.e.,

$$\Vert {\rm N}\Vert =\frac{\left(\mathit{p},\hspace{0.17em}\tilde{\mathit{p}}\right)}{\left|\mathit{p}\right|\left|\tilde{\mathit{p}}\right|}\hspace{0.17em},$$

(31)

where $\mathit{p}=\left(p\left(n,k,r\right)\hspace{0.17em}:\hspace{0.17em}\left(n,k,r\right)\in E\right)$ and $\tilde{\mathit{p}}=\left(\tilde{p}\left(n,k,r\right)\hspace{0.17em}:\hspace{0.17em}\left(n,k,r\right)\in E\right)$ are vectors of exact and approximate values of steady-state probabilities, respectively, $(\mathit{p},\tilde{\mathit{p}})$ denotes the dot product of vectors $\mathit{p}$ and $\tilde{\mathit{p}}$; $\left|\mathit{p}\right|$ and $\left|\tilde{\mathit{p}}\right|$ are the magnitudes of vectors $\mathit{p}$ and $\tilde{\mathit{p}}$ respectively.

Note that, as a rule, the cosine similarity norm is used to determine the orientation of two vectors, and not to compare their values. However, in our case, this measure adequately estimates the proximity of the end points of the vectors $\mathit{p}$ and $\tilde{\mathit{p}}$, since according to the normalizing condition we have ${\displaystyle \sum _{\left(n,k,r\right)\in E}p\left(n,k,r\right)}={\displaystyle \sum _{\left(n,k,r\right)\in E}\tilde{p}\left(n,k,r\right)}=1$; in other words, the endpoints of these vectors are in the same hyper plane.

First, consider the case when the values of the rates of the incoming flow are fixed, but the rates of service for different types of servers change. In this case, in the experiments carried out, the initial data of the hypothetical model are determined as follows. The GM of the MC with three states that controls the intensity of the incoming MMPP flow is defined as:

$${\mathsf{\Sigma }}_1=\Vert \begin{array}{rrr}\hfill -34& \hfill 20& \hfill 14\\ \hfill 18& \hfill -32& \hfill 14\\ \hfill 4& \hfill 16& \hfill -20\end{array}\Vert$$

The values of incoming flow rates are $\mathsf{\Lambda }=\left(15,\hspace{0.17em}10,\hspace{0.17em}5\right)$. The parameters of Bernoulli schemes were chosen as follows: $\alpha =0,\hspace{0.17em}8\hspace{0.17em};\hspace{0.17em}\beta =0,\hspace{0.17em}4\hspace{0.17em};\hspace{0.17em}\sigma =0,\hspace{0.17em}1\hspace{0.17em}.$

Comparative analysis of the results of the approximate algorithm and the simulation method for this series of computational experiments are shown in

Table 1. Accuracy of the developed algorithm for calculation of steady-state probabilities. Case of fixed values of intensities of Markov Modulated Poisson Process (MMPP) flow.

$\left({\mathit{\mu }}_{\mathit{f}},\hspace{0.17em}{\mathit{\mu }}_{\mathit{s}}\right)$	$\mathit{M}$	Values of Norm (31)
(60, 45)	1	0.9945
	2	0.9940
	3	0.9914
(60, 50)	1	0.9931
	2	0.9948
	3	0.9927
(60, 55)	1	0.9943
	2	0.9951
	3	0.9954
(65, 45)	1	0.9936
	2	0.9944
	3	0.9936
(65, 50)	1	0.9941
	2	0.9958
	3	0.9953
(65, 55)	1	0.9952
	2	0.9949
	3	0.9958
(70, 45)	1	0.9946
	2	0.9955
	3	0.9950
(70, 50)	1	0.9944
	2	0.9968
	3	0.9960
(70, 55)	1	0.9964
	2	0.9970
	3	0.9960

. We also performed computational experiments for the cases when the values of the rates of the incoming flow and the rate of service of different types of servers change. In three series of experiments, the parameters of the Bernoulli scheme remained unchanged. In the first series of experiments, the GM of MC with three states that controls the intensity of the incoming MMPP flow was chosen in the same way as before (i.e., it is equal to ${\mathsf{\Sigma }}_1$); in this case, the values of the intensities of the incoming flow were chosen as ${\mathsf{\Lambda }}_1=\left(10,\hspace{0.17em}5,\hspace{0.17em}3\right)$. In the second and third series of experiments, the GM of MC with three states that control the intensities of the incoming MMPP flow and the corresponding values of the intensities of the incoming flow were chosen as follows:

$$\begin{gathered} {\mathsf{\Sigma }}_2=\Vert \begin{array}{rrr}\hfill -25& \hfill 15& \hfill 10\\ \hfill 13& \hfill -23& \hfill 10\\ \hfill 15& \hfill 11& \hfill -26\end{array}\Vert \hspace{0.17em}\hspace{0.17em},\hspace{0.17em}{\mathsf{\Lambda }}_2=\left(15,\hspace{0.17em}12,\hspace{0.17em}7\right)\hspace{0.17em}; \\ {\mathsf{\Sigma }}_3=\Vert \begin{array}{rrr}\hfill -27& \hfill 12& \hfill 15\\ \hfill 9& \hfill -29& \hfill 20\\ \hfill 8& \hfill 15& \hfill -23\end{array}\Vert \hspace{0.17em}\hspace{0.17em},\hspace{0.17em}{\mathsf{\Lambda }}_3=\left(20,\hspace{0.17em}10,\hspace{0.17em}5\right)\hspace{0.17em}. \end{gathered}$$

The results of these experiments are shown in

Table 2. Accuracy of the developed algorithm for calculation of steady-state probabilities. Case of varied values of intensities of MMPP flow.

$\mathsf{\Lambda }$	$\left({\mathit{\mu }}_{\mathit{f}},\hspace{0.17em}{\mathit{\mu }}_{\mathit{s}}\right)$	$\mathit{M}$	Values of Norm (31)
(10, 5, 3)	(60, 45)	1	0.9979
		2	0.9979
		3	0.9973
	(65, 50)	1	0.9983
		2	0.9984
		3	0.9971
	(70, 55)	1	0.9986
		2	0.9986
		3	0.9986
(15, 12, 17)	(60, 45)	1	0.9924
		2	0.9934
		3	0.9914
	(65, 50)	1	0.9953
		2	0.9949
		3	0.9945
	(70, 55)	1	0.9955
		2	0.9962
		3	0.9964
(20, 10, 5)	(60, 45)	1	0.9892
		2	0.9868
		3	0.9867
	(65, 50)	1	0.9903
		2	0.9904
		3	0.9917
	(70, 55)	1	0.9933
		2	0.9938
		3	0.9936

. Table 1 and Table 2 show that the developed algorithm has high accuracy, since the value of the norm (31) is practically equal to one.

Remark 5.

Since the system is infinite-dimensional (with respect to the second and third components of the state vector), when calculating the norm (31), the maximum number of calls of each type in the system from above is limited by sufficiently large finite values. Such a replacement is justified, since when these values exceed certain (sufficiently large) values, the corresponding state probabilities become infinitely small values, i.e., practically equal to zero.

The developed approximate formulas also make it possible to study the behavior of the similarity norm (31) with respect to changes in any (structural and load) parameters of the system. Due to the limited scope of work and for the sake of concreteness, only the results are presented here that show the behavior of this quantity with respect to the change in the parameter $\sigma$. The corresponding results are shown in Figure 1. It can be seen from the graph that with an increase in the value of this parameter, the accuracy of the formulas for calculating the approximate values of the steady-state probabilities of the original 3D MC deteriorates. In other words, the smaller the value of the specified parameter, the higher the accuracy of the proposed formulas. However, even in the worst cases, when the values of this parameter are close to unity, the values of norm (31) turn out to be greater than 0.9, i.e., the developed approximate formulas have a sufficiently high accuracy.

This behavior of norm (31) with respect to a parameter has a completely logical explanation. Indeed, it is known from the theory of space merging of MCs [3] that the lower the intensity of transitions between split classes of states, the higher the accuracy of the space merge algorithms. For the model studied here, at small values of the parameter $\sigma$, the intensities of transitions between split classes $E_n^r\hspace{0.17em},\hspace{0.17em}r=0,\hspace{0.17em}1,\hspace{0.17em}\dots$ (see Formula (8)) turn out to be small, and with its growth the intensity of transitions between split classes increases.

Note that a similar behavior of norm (31) is observed with respect to changes in the parameters $\alpha$ and $\beta$, i.e., with increasing of parameter $\alpha$ the value of the norm is systematically increasing, while with increasing of parameter $\beta$, on the contrary, it systematically decreases (note that, in contrast to parameter $\sigma$, ranges of change in parameters $\alpha$ and $\beta$ and are determined from the ergodicity condition for the system, see (21)).

Behavior of performance measures relative to changes in the threshold parameter $M$, with the above parameters of the MMPP flow and with $\left({\mu }_f,\hspace{0.17em}{\mu }_s\right)=\left(70,\hspace{0.17em}55\right),\hspace{0.17em}\hspace{0.17em}\sigma =0.4$, are shown in Figure 2. The function $L_p$ is non-decreasing (Figure 2a), which should be expected, since with an increase in the parameter $M$ the chances of initial calls joining their queue grows. For the same reason, the functions $L_s$ (Figure 2b) and $R_{ps}$ (Figure 2c) are non-increasing.

Note that with increasing parameter $\beta$ the load of the S-server increases, and thus the chances of p-calls being served in this server decreases, and therefore, the number of such calls in front of the F-server increases (see Figure 2a); with increasing parameter $\alpha$ the load of the S-server decreases, and thus, the number of such calls in the system decreases (see Figure 2b). The values of function R_ps at different ratios of parameters $\alpha$ and $\beta$ almost do not differ from each other (see Figure 2c).

In the numerical experiments discussed above, the autocorrelation of MMPP flow is close to zero. However, the great interest represents the estimation of accuracy of the developed approximate algorithm for the model with MMPP flow for which arrivals are strongly autocorrelated. For this reason, we consider a model where the infinitesimal matrix of a Markov chain with two states that control intensities of MMPP flow is defined as follows:

$${\mathsf{\Sigma }}_4=\Vert \begin{array}{ll}-0.002\hfill & 0.002\hfill \\ 0.001\hfill & -0.001\hfill \end{array}\Vert$$

Let us say that the vector of intensities is equal to $\mathsf{\Lambda }=\left(1,\hspace{0.17em}0.1\right)$. Note that this MMPP flow has an autocorrelation of about 0.4 for the lag−1.

The parameters of the Bernoulli schemes were chosen as above, i.e., $\alpha =0.8\hspace{0.17em};\hspace{0.17em}\beta =0.4\hspace{0.17em};\hspace{0.17em}\sigma =0.\hspace{0.17em}1\hspace{0.17em}.$ Results of numerical experiments for this case are given in

Table 3. Accuracy of the developed algorithm for calculation of steady-state probabilities. Case of autocorrelated arrivals of MMPP flow.

${\mathit{\mu }}_{\mathit{f}}$	${\mathit{\mu }}_{\mathit{s}}$	$\mathit{M}$	Values of the Norm (31)
60	40	1	0.9908
	40	3	0.7369
	45	2	0.9999
	50	1	0.9973
	55	2	0.9692
	55	3	0.7635
65	40	3	0.8021
	45	1	0.9749
	45	2	0.9917
	45	3	0.9839
	50	2	0.7357
	55	2	0.9699
	55	3	0.7803
70	40	2	0.8742
	45	2	0.9995
	45	3	0.9955
	50	1	0.8454
	50	2	0.9950
	50	3	0.9981
	55	2	0.9930
75	40	2	0.9990
	45	2	0.9543
	55	1	0.7905
	55	2	0.8041
	55	3	0.9967

. Comparative analysis of Table 3 with Table 1 and 2 shows that the accuracy of the developed approximate algorithm for such MMPP flow is slightly lower than for an MMPP flow with small autocorrelation coefficient. At the same time, it has a sufficiently high accuracy for practical applications.

5. Conclusions

The paper studies a model of a system with heterogeneous servers, MMPP flow, and instant feedback. After the completion of servicing in the high-speed server, primary calls according to the Bernoulli scheme either leave the system or immediately require re-servicing. Repeated (feedback) calls are served in a low-speed server, and after the service is completed, feedback calls can be repeated many times. If, at the time of the arrival of the primary call, the low-speed server is empty and the number of calls in the queue in front of the high-speed server exceeds a certain threshold value, then the incoming call, according to the Bernoulli scheme, is either sent for service to the low-speed server or joins its own queue. It is believed that primary and feedback calls can form queues of infinite or finite length.

For the case of infinite buffers, it is shown that the mathematical model of the system under study is a certain three-dimensional Markov chain with an infinite-dimensional state space. The ergodicity condition for the model is found and an approximate algorithm for calculating the steady-state probabilities of the corresponding Markov chain is developed. Using simulation experiments, the high accuracy of the developed algorithm has been shown.

A direction of further research might be indicated via investigation of similar model with MAP flow and PH(phase type) distribution of service time as well as the solving of the problem of finding the optimal value (in some sense) of the introduced threshold parameter.