Retrial Queues with Unreliable Servers and Delayed Feedback

In this paper, models of unreliable multi-server retrial queues with delayed feedback are examined. The Bernoulli retrial is allowed upon the arrival of both primary (from outside) and feedback customers (from orbit), as well as the Bernoulli feedback that may occur after each service in this system. Servers can break down both during the service of customers and when they are idle. If a server breaks down during the service of a customer, then the interrupted customer, in accordance with the Bernoulli scheme, decides either to leave the system or join a common orbit of retrial and feedback customers. An approximate method, based on the space merging approach of three-dimensional Markov chains, is proposed for the calculation of the steady-state probabilities, as well as performance measures of the system. The results of the numerical experiments are demonstrated.


Introduction
The theory of retrial queue (RQ) is widely implemented for modeling the functioning of real-world systems. State of the art theory of RQ might be found in the recent review [1]. The general property of RQ is the following: if upon arrival of a primary customer, all servers of the system and/or waiting rooms are full, then a primary customer can join the group of customers that are waiting for the service outside of the system. This virtual waiting room is called orbit for retrial customers. As a rule, in the classical theory of RQ, it is assumed that the system has no information about the number of retrial customers in the orbit, and retrial customers generate their requests independently of each other.
In classical RQ, both primary and retrial customers leave the system after getting the required service. However, in many cases, some customers require an arbitrary number of optional services aside from the essential services. In such systems, some customers provide feedback to the system for optional service. These customers are referred to as feedback customers. Note that a customer can provide feedback either instantaneously after leaving the server or after some delay. Retrial queues with delayed feedback (RQwDFB) are very useful tools for modeling stochastic processes arising in communication networks, call centers, queueing-inventory systems, etc. This paper analyzes the steady-state behavior of an unreliable multi-server retrial queueing system, first with essential services and then with an arbitrary number of optional services.
In the above-mentioned papers, the models of RQwDFB with reliable servers are investigated. However, in practice, the servers of the systems are subject to random

Model Description
A pictorial representation of the investigated RQwDFB with unreliable servers is shown in Figure 1.
First, consider the model with a finite capacity of orbit for retrial customers, i.e., assume that the maximal size of orbit is equal to R, R < ∞. The system contains K > 1 identical and unreliable servers. It is assumed that servers can independently break down in both busy and idle modes. The failure rate of servers is different in their various statuses, i.e., failures of servers in both busy and idle statuses occur via two independent Poisson processes with rates of θ 1 and θ 2 , respectively. Each server also has its own repairman and the repair times for servers are exponentially distributed with the common rate ξ. If the failure of a server occurs during the busy status (active breakdown), then a customer that is preempted due to a server failure either leaves the system (i.e., is lost) with a probability (w.p.) γ r or enters the retrial orbit w.p. 1 − γ r , provided that there are r customers in the orbit, r = 0, 1, . . . , R − 1 ; γ R = 1. First, consider the model with a finite capacity of orbit for retrial customers, i.e., assum that the maximal size of orbit is equal to . The system contains 1 > K identi and unreliable servers. It is assumed that servers can independently break down in bo busy and idle modes. The failure rate of servers is different in their various statuses, i failures of servers in both busy and idle statuses occur via two independent Poiss processes with rates of 1 θ and 2 θ , respectively. Each server also has its own repairm and the repair times for servers are exponentially distributed with the common rate ξ the failure of a server occurs during the busy status (active breakdown), then a custom that is preempted due to a server failure either leaves the system (i.e., is lost) with probability (w.p.) r γ or enters the retrial orbit w.p.
Primary customers (p-customers) form Poisson flow with rate λ . An incoming customer is accepted for service if at this moment the total number of busy and fail servers is less than K ; otherwise, they either leave the system w.p. r α or go to or w.p.
provided that there are r customers in the orbit, ; 1 ..., , Primary customers (p-customers) form Poisson flow with rate λ. An incoming pcustomer is accepted for service if at this moment the total number of busy and failed servers is less than K; otherwise, they either leave the system w.p. α r or go to orbit w.p. 1 − α r provided that there are r customers in the orbit, r = 0, 1, . . . , R − 1 ; α R = 1.
Sojourn times of retrial customers (r-customers) in the orbit are exponentially distributed with a mean η −1 and here a classical retrial scheme is considered, i.e., if there are r customers in the orbit, then the arrival intensity of r-customers is rη (linear retrial rate). If upon the arrival of r-customer there is no free and serviceable server, then they either leave the orbit (i.e., are lost) w.p. σ r or return to the orbit w.p. 1 − σ r provided that there are r customers in the orbit, r = 1, . . . , R.
After the completion of the service process, the customer either leaves the system forever w.p. β r or goes to the orbit for repeated service w.p. 1 − β r provided that there are r customers in the orbit, r = 0, 1, . . . , R − 1 ; β R = 1. customers in the orbit, r = 0, 1, . . . , R − 1; β = 1. Customers who require re-service are called feedback customers (f -customers). For simplicity, in the orbit, both kinds of customers (i.e., r-customers and f -customers) are treated as the same and they are referred to as r-customers. Here, multiple re-servicing processes (feedbacks) are allowed.
For the model with an infinite orbit capacity for retrial and feedback customers, i.e., when R = ∞, we assume that α r < 1 , β r < 1 , γ r < 1 for any value of r.

of 23
The service times for all kinds of customers (p-, r-, and f -customers) are assumed to be exponentially distributed with a common mean µ −1 .
All times involved in the model, as well as the arrival of p-customers from the outside, the arrival of r-customers from the orbit, and the service of customers, are mutually independent of each other.
Our goal is to determine the stationary distribution of the given system, as well as to calculate its performance measures: loss probabilities of both p-customers and r-customers, the average number of both failed and busy servers, and the average number of r-customers in the orbit.

Balance Equations Method for the Model with a Finite Orbit Size
The state of the system is defined by a three-dimensional vector (b, f , r), where b is the number of busy servers, f is the number of failure servers, and r is the number of r-customers in the orbit. In other words, the mathematical model of the given system is 3D MC with the following state space (SS): where It is clear that in the state (b, f , r) ∈ S, the number of idle and serviceable servers is equal to K − b − f . From (1), we conclude that the SS of the investigated 3D MC geometrically represents a prism with height R. Consider the problem of constructing the generator of the indicated 3D MC.
Let us fix the value of r and first consider determining transition intensities between the states within a subset (level), S r . In order to be short, a transition from . Let the initial state be (b, f , r) ∈ S. As a result of analyzing possible transitions, we have the following relations: is carried out upon departure of the customer from the system or their loss due to a server's active breakdown, where µ r = µβ r + θ 1 γ r ; is carried out when the failure of a server occurs in its idle status; is carried out when the repair of a failed server is completed; Now consider determining the transition intensities between the states of the various subsets (levels). Note that only transitions between states of neighboring levels are allowed. For such kinds of transitions, we have the following relations: is carried out upon the joining of the arrived p-customer to orbit when b + f = K and r < R; is carried out when the serviced customer feedbacks to the orbit when r < R; → (b, f + 1, r + 1) is carried out when the customer goes to orbit due to a server's active breakdown when r < R; is carried out when an arrived r-customer is lost, if b + f = K, when r > 0.
The transition intensity from state (b, f , r) to state (b , f , r ) is denoted by q((b, f , r), (b , f , r )). From the above relations, we conclude that the positive elements of the generator of the investigated 3D MCs are determined as follows: Hereinafter, the equality of vectors means that their corresponding components are equal to one another.
From the relations (2), we conclude that states of constructed finite 3D MCs are communicated with each other, i.e., a stationary mode exists in this system. Let p(b, f , r) denote the stationary probability of the state (b, f , r) ∈ S. These probabilities can be found using the following balance equations: where S + (b, f ,r) is the set of those states of S, which can be reached from the state (b, f , r) in one step, and S − (b, f ,r) is the set of those states of S, from which one can get to the state (b, f , r) in one step. Note that the sets S + (b, f ,r) and S − (b, f ,r) are determined from the relations (2). The dimension of the balance Equations (3) and (4) is equal to (K + 1)(K + 2)(R + 1)/2. Modern software packages allow the linear equations of arbitrary (finite) dimensions to be solved.
The desired performance measures are calculated via steady-state probabilities. The average number of busy (NBS av ) and failed servers (NFS av ), as well as the average number of r-customers in the orbit (NCO av ), are calculated as the mathematical expectations of the corresponding random variables, i.e., The loss probabilities of p-customers P p and r-customers (P r ) are calculated as Hereinafter δ(i, j) are Kronecker's delta.

The Space Merging Method for the Model with an Infinite Orbit Size
Methodological difficulties occur in solving the investigated problem for the model with an infinite orbit size, i.e., when R = ∞. Now consider this problem for the case where the indicated above probabilities, α r , β r , γ r , σ r , are state-independent, i.e., assume that The mathematical model of the given system is 3D MC with the following state space (SS): where (10), we conclude that the SS of the investigated 3D MC geometrically represents an infinite prism. The generator of the presently examined infinite 3D MC is determined similarly to (2), wherein the right side condition, r < R, should be omitted.
For simplicity, let p(b, f , r) again denote the stationary probability of the state (b, f , r) ∈ E. Below, we show that for the case of the linear retrial rate, the stationary mode exists for any positive values of input parameters; however, in the case of the constant retrial rate, the fulfillment of some stability conditions is required.
From the form of the SS (8), we conclude that the constructed 3D MC does not belong to the class of quasi-birth and death (QBD) processes (see [26]). Therefore, unfortunately, to calculate the steady-state probabilities of the investigated 3D MC, we cannot apply well-worked matrix-geometric methods. For this reason, we propose an alternate approach below to solve the indicated problem. The proposed approach is based on the space merging principles of a multi-dimensional MC and it allows for the calculation of steadystate probabilities and performance metrics via explicit formulas. Moreover, the obtained formulas contain quantities that are tabulated, i.e., Erlang's B-formula as well as the formula for calculating the average number of busy servers in Erlang's classical model.
It is known that satisfying the following condition is required for the correct application of the SMM: the state space of the system should be decomposed into classes in such a way that the rates of transitions between states within classes are much larger than rates of transitions between states from different classes. The fulfillment of this condition can be ensured at certain ratios between the initial parameters of the system under study. So, suppose that the arrival rate of p-customers is much larger than the arrival rate of r-customers, i.e., η << λ.
When this assumption is fulfilled, the transitions between the states within classes E r (see (10)) are much larger than transition classes between the states of different classes. By using this fact, the following merge function in the SS (10) is defined: where < r > is the merged state, which includes all micro-states from the class E r , r = 0, 1, . . . The set of merged states < r > is denoted by According to hierarchical SMM (see [7,27]), the approximate values of the steady-state probabilities of the initial 3D MC, denoted by p(b, f , r), are calculated as follows: where ρ r (b, f ) is the probability of the state, (b, f , r), within a splitting model with SS E r , and π 1 (< r >) is the probability of the merged state < r >∈ Ω 1 . From (12), conclude that to calculate the steady-state probabilities of the initial 3D MC, it is necessary to find the steady-state probabilities, ρ r (b, f ), of the infinite number of 2D MCs with a finite SS as well as the steady-state probabilities, π 1 (< r >), of one 1D MC with the SS Ω 1 .
First, consider the calculation of the steady-state probabilities, ρ r (b, f ), of the 2D MCs with the SS E r . The transitions diagram for the splitting model with the SS E r is shown in Figure 2. , do not depend on r because, below, w fix the value of r and omit the subscript r in the notation of these quantities.
The calculation of the steady-state probabilities, ( ) f b, ρ , for a moderately sized S r S , might be performed via solving balance equations. However, it is possible to develo explicit formulas for their approximate calculations. For this purpose, the mergin procedure (the second level of the hierarchy) is applied to these 2D MCs.
For the correct application of the SMM in the second level, it is assumed th 2 1 θ θ >> . Note that this assumption is a realistic one since, in real situations, the failur intensity of servers in a busy mode is much higher than for those in an idle mode. Unde this assumption, consider the following splitting of r S : ..., , 1 , 0 : , The transition intensity from state (b, f , r) ∈ E r to state (b , f , r) ∈ E r is denoted by q r ((b, f , r), (b , f , r)). From (2), we conclude that these quantities are determined as follows: where µ = µβ + θ 1 γ. From (13), we conclude that q r ((b, f , r), (b , f , r)) are identical for any value of r, i.e., the steady-state probabilities, ρ r (b, f ), do not depend on r because, below, we fix the value of r and omit the subscript r in the notation of these quantities.
The calculation of the steady-state probabilities, ρ(b, f ), for a moderately sized SS, S r , might be performed via solving balance equations. However, it is possible to develop explicit formulas for their approximate calculations. For this purpose, the merging procedure (the second level of the hierarchy) is applied to these 2D MCs.
For the correct application of the SMM in the second level, it is assumed that θ 1 >> θ 2 . Note that this assumption is a realistic one since, in real situations, the failure intensity of servers in a busy mode is much higher than for those in an idle mode. Under this assumption, consider the following splitting of S r : where Similar to (11), based on splitting (14) in the SS, E r , the following merge function is defined: where << f >> is a merge state which includes all micro-states from the subclass E f r . The set of merged states << f >> is denoted by Ω 2 = {<< f >> : f = 0, 1, . . . , K}.
In accordance with [7,27], we have: where ρ f (b) is the probability of the state (b, f ) within the splitting model with space E f r , and π 2 (<< f >>) is the probability of the merged state << f >>∈ Ω 2 .
In a subclass, E  (2) we conclude that these quantities are independent of the parameter r, and are determined as follows: From (16), we conclude that the state probabilities, ρ f (b), of the splitting models with the SS E f r , f = 0, 1, . . . , K − 1, coincides with the state probabilities of Erlang's model, M/M/K-f/K-f, with load ν = λ/(µβ + θ 1 γ), i.e., the desired probabilities are calculated by well-known Erlang's formulas: The splitting model with the SS E K r contains only one state (0, K, r), i.e., for this model, we set ρ K (0) = 1.
By using (12), (15), (17), and (21), the approximate values of the steady-state of the initial 3D MC are found. So, for any positive values of the loading parameters in this system, there exists stationary mode, and after some standard mathematical transformations, we obtain the following approximate formulas for calculating the performance measures: Note. The average number of busy (NBS av ) and failed servers (NFS av ) (see Formulas (24) and (25)) does not depend on the parameters η and σ. These facts are explained by the accepted assumption that η << λ, i.e., the impact of the indicated parameters on two performance measures is not essential. Other performance measures depend on all system parameters.

Numerical Results
In this section, the qualitative behavior of the performance measures is explored through a few numerical experiments. Due to the size limitations for the present paper, only the results for the model with an infinite orbit size have been considered here.
It can be noticed that the proposed algorithms allowed for a simple investigation of the behavior of the performance measures over any parameter of the model without any computational difficulties. For the sake of brevity, we have omitted the results of numerical experiments that demonstrate the high accuracy of the algorithm developed for the approximate calculation of steady-state probabilities of the system under study (appropriate results might be found in [9]). Moreover, for brevity, we only demonstrate experiments that present the behavior of performance measures with respect to the number of servers over various values of the probabilities α, β, γ and σ.
For the numerical calculations, the values of the model parameters are selected as follows: λ = 5, µ = 3, θ 1 = 5, θ 2 = 1, ξ = 1, η = 1 . Note that since the performance measures depend on many parameters, the analysis of the numerical experiments indicated below concerns only the selected values of the initial data. However, in some cases, the conclusions are general.  Figure 3a,b that both PB p and PB r decrease with respect to K, where an increase of α is favorable for PB r , whereas PB p increases with respect to α; for large values of K, i.e., when K ≥ 7, the impact of the values of α is negligible. It is seen from Figure 3c,d that both the NBS av and the NFS av increase, and their values are almost independent of α. Figure 3e exhibits the impact of K and α over the NCO av . It is seen from the graph that the values of the decrease with respect to both K and α. As in Figure 3a,b, here, when K ≥ 7, the impact of the values of α on the values of the NCO av is negligible. decrease with respect to both K and α . As in Figure 3a,b, here, when K ≥ 7, the impact of the values of α on the values of the av NCO is negligible.     PB is decreasing by one with respect to β . As was expected, the av NBS increases with respect to K and decreases when β increases (see Figure 4c). Figure 4d shows that the av NFS increases with respect to K, and its values are almost independent of β . However, it is seen from Figure 4e Figure 4a,b that both PB p and PB r decrease, and their values are almost independent on β; some dependence on β is observed for PB r under small values of K, i.e., under small values of K, the performance measure PB r is decreasing by one with respect to β. As was expected, the NBS av increases with respect to K and decreases when β increases (see Figure 4c). Figure 4d shows that the NFS av increases with respect to K, and its values are almost independent of β. However, it is seen from Figure 4e that the values of the NCO av increase with respect to K and decrease when β increases, and that the impact of the values of β on the values of the NCO av is essential for a large value of K.        Figure 5a,b, we see that both PB p and PB r decrease with respect to K, and their values are almost independent on γ. Surprisingly, the NBS av increases when γ increases (see Figure 5c). Such behavior of the NBS av can be explained as follows: for selected values of the load parameters of the system, an increase in the probability of a customer leaving the system from a broken server leads to an increase in the chance of p-customers being accepted, i.e., the average number of busy servers increases. Though, for other values of the load parameters, this performance measure might decrease when γ increases. Figure 5d demonstrates that the NFS av increases with respect to K and its values are almost independent of γ. The impact of γ on the values of the NCO av is negligible, and it decreases when K increases (see Figure 5e).    = 0.7. Plot 6(a) exhibits that PB p is almost not affected by an increment in σ, and it decreases with the increasing number of servers. Conversely, plot 6(b) exhibits that, for a small value of K, PB r is affected by an increment in σ, but when K ≥ 7, the impact of the values of σ is negligible and with increasing σ, PB r increases. Note that both the NBS av and the NFS av increase with respect to K, and the values are almost independent of σ (see Figure 5c,d). It is seen from Figure 6e that for a small value of K, the NCO av is affected by an increment in σ, but when K ≥ 4, the impact of the values of σ is negligible and with increasing σ, the NCO av decreases.

Conclusions
The present paper proposes mathematical models of unreliable multi-server RQwDFB. We studied both models with common finite and infinite orbit capacity for retrial and feedback customers. In a model with a finite orbit capacity, it is assumed that the probabilities of retrial, feedback, and balking from the orbit depend on the current number of customers in the orbit. For such a kind of model, we propose an exact method for calculating its steady-state probabilities as well as it performance measures. For a model with an infinite orbit size, an approximate space merging method is developed for calculating its steady-state probabilities as well as its performance measures.
As for directions for further research, it is proposed that models of RQwDFB with MAP flow and phase-type distributions of service time be studied. Practical interests are represented by the problems of optimizing RQwDFB regarding the selected criteria for the quality of their functioning. These problems are the subjects of future works.

Conclusions
The present paper proposes mathematical models of unreliable multi-server RQwDFB. We studied both models with common finite and infinite orbit capacity for retrial and feedback customers. In a model with a finite orbit capacity, it is assumed that the probabilities of retrial, feedback, and balking from the orbit depend on the current number of customers in the orbit. For such a kind of model, we propose an exact method for calculating its steady-state probabilities as well as it performance measures. For a model with an infinite orbit size, an approximate space merging method is developed for calculating its steady-state probabilities as well as its performance measures.
As for directions for further research, it is proposed that models of RQwDFB with MAP flow and phase-type distributions of service time be studied. Practical interests are represented by the problems of optimizing RQwDFB regarding the selected criteria for the quality of their functioning. These problems are the subjects of future works.