MARKOVIAN RETRIAL QUEUES WITH TWO WAY COMMUNICATION

. In this paper, we ﬁrst consider single server retrial queues with two way communication. Ingoing calls arrive at the server according to a Poisson process. Service times of these calls follow an exponential distribution. If the server is idle, it starts making an outgoing call in an exponentially distributed time. The duration of outgoing calls follows another exponential distribution. An ingoing arriving call that ﬁnds the server being busy joins an orbit and retries to enter the server after some exponentially distributed time. For this model, we present an extensive study in which we derive explicit expressions for the joint stationary distribution of the number of ingoing calls in the orbit and the state of the server, the partial factorial moments as well as their generating functions. Furthermore, we obtain asymptotic formulae for the joint stationary distribution and the factorial moments. We then extend the study to multiserver retrial queues with two way communication for which a necessary and suﬃcient condition for the stability, an explicit formula for average number of ingoing calls in the servers and a level-dependent quasi-birth-and-death process are derived.

1. Introduction. Recently, retrial queues are paid much attention because they have applications in performance analysis of various systems such as call centers, computer networks and telecommunication systems [3,12,17,28]. Retrial queues are characterized by the fact that customers (i.e., calls) that cannot receive service upon arrival enter a virtual orbit and retry for service again after some random time. The arrival flow from the orbit makes the underlying Markov chain of retrial queues to be nonhomogeneous. As a result, analysis of retrial queues is much more difficult than that of the corresponding queueing models without retrials and explicit results are obtained only in a few special cases [3,12,23,24].
Hypergeometric functions and their special versions play an important role in the derivation of analytical solutions for retrial queues. In fact, the stationary characteristics of the system state of the conventional M/M/1/1 retrial queue are expressed in terms of special hypergeometric functions [3,12,24]. A review of the The rest of the paper is organized as follows. Section 2 describes the model in detail. Section 3 is devoted to the main results of this paper in which an extensive study of the M/M/1/1 retrial queue with two way communication is presented. In Section 4, we discuss an extension to a multiserver retrial queue with two way communication and obtain some explicit results. Finally, we conclude our paper and present some future research topics in Section 5.
2. Model description and preliminaries. In this section, we present the mathematical description of the single server retrial queue with two way communication in detail and provide some preliminaries which will be used in the main results presented in Section 3. We separate the multiserver model to Section 4 because the methodology for the multiserver model is different from that for the single server model.
2.1. Queueing model. We consider a single server retrial queue with two way communication. Primary ingoing calls arrive at the server according to a Poisson process with rate λ. An ingoing call that sees the server being busy enters an orbit and retries to occupy the server after an exponentially distributed time with mean 1/µ. In addition, we assume that if the server is idle then it makes an outgoing call after an exponentially distributed time with mean 1/α. The service times of the ingoing and the outgoing calls are exponentially distributed with mean 1/ν 1 and 1/ν 2 , respectively. See Figure 1 for transitions among states. In what follows, we consider the system under the stability condition which will be derived later. Furthermore, let π i,j = lim t→∞ Pr(S(t) = i, N (t) = j), i = 0, 1, 2, j ∈ Z + , denote the joint stationary distribution of the system state.
Summing up equations (4)- (6) and rearranging the result, we obtain Dividing both sides of the above formula by (z − 1) yields Note that equation (7) represents a balance between the flows coming into and out the orbit.

Hypergeometric functions.
In this section, we give a brief summary on hypergeometric functions, which will be used to obtain our explicit results in the sequel. For a complex number x, let denote the Pochhammer symbol, where N = {1, 2, . . . }. Then, for complex numbers a, b, c and z, the hypergeometric function F (a, b; c; z) is defined as It should be noted that where the last equality follows from the generalized Newton binomial formula. This formula will be frequently used throughout the paper.
Proposition 2.1 (Theorem VI.12, p. 434 in [13]). Let a(z) = ∞ n=0 a n z n and b(z) = ∞ n=0 b n z n denote two power series with radii of convergence r a > r b ≥ 0, respectively. Assume that b(z) satisfies the test Then the coefficients of the product g(z) = a(z)b(z) satisfy provided that a(r b ) = 0, where [z n ]g(z) denotes the coefficient of z n in the power series expansion of g(z) and x n ∼ y n is defined by lim n→∞ x n /y n = 1.
Proposition 2.2 (See p. 377 in [13]). For any complex number a whose real part is positive, we have where Γ(a) is the Euler Gamma function defined as Corollary 2.1. For any complex number a whose real part is positive and positive number γ, we have Proof. Let x = γz. We have where the ∼ follows from Proposition 2.2.
Proposition 2.3. Let {a n , b n , a n , b n ; n ∈ Z + } denote sequences of real numbers such that a n ∼ a n , b n ∼ b n , n → ∞, lim n→∞ a n b n = 0.
Then, we have a n + b n ∼ b n as n → ∞.
Proof. The result is straightforward from the definition. Indeed, we have We apply these three propositions to derive asymptotic results in Section 3.6.
3.1. Generating functions. First, we derive explicit expressions for the partial generating functions.
Theorem 3.1. Explicit expressions for the partial generating functions are as follows: where , Remark 3.1. We observe that Π 0 (z) is expressed in terms of a product of two special hypergeometric functions. Note also that ρ and σ denote the traffic intensity of ingoing calls and the traffic intensity of outgoing calls, respectively.

Stationary distribution.
Our goal in this section is to derive explicit expressions for {π i,j ; i = 0, 1, 2, j ∈ Z + }.
Theorem 3.2. Explicit expressions for the stationary distribution are given by for j ∈ Z + .

Factorial moments. We now deal with the partial factorial moments
The moments of order k = 0 are trivially given by . It is easy to see that and therefore M i k (k ∈ Z + ) can be obtained from the coefficient of z k in the series Π i (1+z). Using the same techniques as used in the proof of Theorem 3.2, we obtain the following results.
Theorem 3.3. The partial factorial moments are given by for k ∈ Z + , where

JESUS R. ARTALEJO AND TUAN PHUNG-DUC
3.4. Recursive formulae. In the previous sections 3.2 and 3.3, we have derived explicit expressions for the stationary distribution and the partial factorial moments. However, since the obtained formulae involve the auxiliary quantities C 1 , C 2 , D 1 and D 2 , which may be either positive or negative, a computational implementation based on these formulae may be numerically unstable. In order to resolve this problem, we next present simple recursive schemes to compute the stationary probabilities and their partial factorial moments. It should be noted that these recursive formulae can be used for both symbolic or numerical implementations.
Theorem 3.4. The stationary probabilities can be computed from the following recursive formulae: We remember that π 0,0 was given in (14).
Theorem 3.5. We have the following recursive formulae for the partial factorial moments: Expressions for the moments of order k = 0 were given in (19).

3.5.
First moments and cost model.
In this section, for the sake of completeness, we summarize some simple formulae for the first moments of the number of customers in the orbit. By combining (19) and (32), we have It follows from (33) and (38) that Finally, from (34), (38) and (39), we obtain The above formulae (38)-(40) are consistent with those derived in the existing literature [5,11].
3.5.2. Cost model. Let U denote the utilization of the server, i.e., where (19) is used in the second equality. From formulae (38)-(40), we also obtain the mean number of customers in orbit which is given by From a management point of view, we need to minimize the idle ratio of the server, 1 − U , but at the same time, from a quality of service (QoS) point of view, we also need to minimize E[N ]. Thus, our objective is to find an optimal σ which satisfies both these needs. For this purpose, we consider the following minimization problem: where A and B are positive costs and ρ, µ, λ and ν 2 are kept constant. Since ν 2 remains constant, to minimize with respect to σ in fact amounts to minimize with respect to the outgoing call rate α. Remark 3.3. Our motivation for the above choice of the cost function is that the outgoing call rate α is directly under the control of the server. The chance for controlling other system parameters (e.g., the ingoing call rate, λ, or the retrial rate, µ) is much more reduced in practice. However, we may consider other static optimization problems by introducing in the cost function other performance measures (e.g., blocking probabilities, waiting time indicators). A more sophisticated dynamic approach based on Markov decision processes may also be considered but it is not the objective of this paper.
We now express f (σ) as .
Thus, we have 3.6. Asymptotic analysis. The limitations of the explicit formulae, as long as they are expressed as finite sums involving positive and negative quantities, were mentioned in Section 3.4. The alternative recursive scheme is helpful for computing the stationary probabilities but it does not provide explicit expressions. In what follows, we supplement the exact and recursive formulae by deriving simple asymptotic formulae both for the stationary distribution and the factorial moments. In addition to their inherent value as limiting results, the asymptotic formulae are simpler than their explicit counterparts (see Theorems 3.6 and 3.7). As a result, a sensitivity analysis of the main performance measures can be carried out easier when it is based on the asymptotic formulae rather than on the explicit ones.
3.6.1. Asymptotic formulae for the stationary distribution.
Theorem 3.6. We have the following asymptotic results for the stationary distribution. If ρ > θ, then we have as n → ∞. For the case θ > ρ, we have as n → ∞.
Remark 3.4. We observe that π 0,n and π 2,n have the same order in the case of ρ > θ, while π 1,n and π 2,n have the same order in the case of ρ < θ.
Remark 3.5. We confirm that our asymptotic results for the case ρ > θ are consistent with those derived by Kim et al. [15] for the conventional M/M/1/1 retrial queue without outgoing calls (i.e., α = 0).
Finally, we derive asymptotic formulae for {π 1,n ; n ∈ Z + }. It follows from (11) and (12) that (50) Applying the same techniques as used in the derivation of (41) for the two components in the right hand side of (50), we find that as n → ∞. Therefore, from (51), (52) and Proposition 2.3, we obtain (43). Asymptotic formulae (44), (45) and (46) for the case θ > ρ can be derived following similar arguments to that used above for the case ρ > θ. Thus, we omit the proof.
We can now derive the following asymptotic results for the factorial moments using the same techniques as used in Section 3.6.1.
Remark 3.6. Similar to the stationary distribution, we also observe that M 0 n and M 2 n have the same order in the case of ρ > θ, while the order of M 1 n and of M 2 n is the same in the case of ρ < θ.
Proof. From (11), (12) and (13), we observe that the structure of Π i (1 + z) (i = 0, 1, 2) is similar to that of Π i (z) where the poles are replaced by ρ −1 and θ −1 instead of ρ −1 and θ −1 . Therefore, using the same arguments as used in the proof of Theorem 3.6 and the fact that yields the first ∼ in (53) to (58). The second ∼ of (53) to (58) follows from the first ∼ and the Stirling formula: n! ∼ √ 2π n n+ 1 2 e n , n → ∞.
3.6.3. Numerical validation of asymptotic formulae. In this section, we present some numerical examples to show the tail asymptotic behavior of the join stationary distribution. We set µ = 1, ν 1 = 1 and α = 1. Figure 2 represents {π i,n ; i = 0, 1, 2, n ∈ Z + } against n for the case λ = 0.9 and ν 2 = 2.5 for which ρ > θ. Figure 3 shows {π i,n ; i = 0, 1, 2, n ∈ Z + } against n for the case λ = 0.1 and ν 2 = 0.01 for which θ > ρ. In both figures, {π i,n ; i = 0, 1, 2, n ∈ Z + } computed by recursive formulae presented in Theorem 3.4 and those calculated by asymptotic formulae in Theorem 3.6 are plotted. We observe in both figures that the curves by asymptotic formulae are well fitted to those by the recursive formulae when the number of customers in the orbit is large. The observation suggests that these asymptotic formulae can be used to estimate exact values with high accuracy for the case of a large number of retrial customers. We also observe from Figure 2 that π 1,n dominates π 0,n and π 2,n and that the curves of π 0,n and π 2,n are asymptotically parallel when n is large. These observations agree with the asymptotic results presented in (41), (42) and (43), where ρ −1 is the dominant pole. On the other hand, Figure 3 shows that the curves of π 1,n and π 2,n are asymptotically parallel and that π 0,n is dominated by π 1,n and π 2,n when n is large. These results are consistent with asymptotic formulae (44), (45) and (46), where θ −1 is the dominant pole. π 0,n exact π 1,n exact π 2,n exact π 0,n asymptotic π 1,n asymptotic π 2,n asymptotic  π 0,n exact π 1,n exact π 2,n exact π 0,n asymptotic π 1,n asymptotic π 2,n asymptotic

4.
A multiserver retrial queue with two way communication.

4.1.
Queueing model. In this section, we present a generalization of our single server model. In particular, we consider an M/M/c/c retrial queue with two way communication, where definitions of λ, α, ν 1 , ν 2 and µ are the same as those of the single server model. An arriving customer that finds all the servers being busy joins the orbit. The behavior of the servers in the M/M/c/c retrial queue is the same as that of the single server model, i.e., each idle server makes an outgoing call in an exponentially distributed time with mean 1/α. For this model, first we obtain an explicit formula for the average number of ingoing calls in the servers. Second, we establish a necessary and sufficient condition for the ergodicity. Third, we express the underlying Markov chain as a level-dependent quasi-birth-and-death (QBD) process.

4.2.
Markov chain and ergodic condition. Let S 1 (t), S 2 (t) and N (t) denote the numbers of ingoing calls and outgoing calls in the servers and the number of ingoing calls in the orbit at time t ≥ 0, respectively. It is easy to see that the process {χ(t); t ≥ 0} = {(S 1 (t), S 2 (t), N (t)); t ≥ 0} forms a Markov chain on the state space S defined by S = {(i, j, k); i = 0, 1, . . . , c, j = 0, 1, . . . , c − i, k ∈ Z + }.   Proof. A proof of Theorem 4.2 is presented in Appendix C.

4.3.
Level-dependent QBD process. It is easy to see that {χ(t); t ≥ 0} forms a level-dependent QBD process, where N (t) and {S 1 (t), S 2 (t)} are referred to as the level and the phase, respectively. The infinitesimal generator of the process is given by where O denotes a matrix with an appropriate size whose entries are all zero, while Q k,k−1 (k ∈ N), Q k,k and Q k,k+1 (k ∈ Z + ) are explicitly given in Appendix D. Since the stability condition and explicit block matrices of the infinitesimal generator are explicitly given, we could apply several approximation methods [22,25] in order to obtain numerical results. However, because the objective of this paper is to derive explicit results, we omit here a detailed numerical analysis.

5.
Conclusion and future work. We have analyzed the M/M/1/1 retrial queue with two way communication in detail. In particular, we have derived explicit expressions for the stationary distribution and the partial factorial moments. We have also derived recursive formulae based on which both numerical and symbolic algorithms can be implemented. In addition, a cost model has been presented in order to find the optimal rate of outgoing calls. Furthermore, some simple asymptotic formulae for the stationary distribution and partial factorial moments have also been obtained.
As for the M/M/c/c retrial queue, we have established a necessary and sufficient condition for the stability and have derived an explicit formula for the average number of ingoing calls in the servers. Furthermore, we have presented a leveldependent QBD process of the model for which numerical analysis could be carried out by several methods presented in literature [22,25].
For the future work, we pay attention to the consideration of impatient customers in multiserver retrial queues with two way communication. Another extension is the consideration of single server retrial queues with two way communication with MAP arrivals and more general service time distributions of ingoing and outgoing calls. comments and suggestions which improved the paper. J.R. Artalejo was supported by the Government of Spain (Ministry of Science and Innovation) and the European Commission through project MTM2008-01121. T. Phung-Duc was supported in part by Japan Society for the Promotion of Science, Grant-in-Aid for JSPS Fellows (No. 22.470).
Appendix A. Special cases. This section presents results for the special cases ν 1 = ν 2 and ν 1 = λ + ν 2 . While the former corresponds just with the particular case where both the ingoing and the outgoing calls receive identical service times, the latter implies some minor mathematical differences.
A.1. The case ν 1 = ν 2 . In the special case ν 1 = ν 2 , the stationary probabilities {π 0,j ; j ∈ Z + } reduce to and their generating function is given by leading to the factorial moments On the other hand, the probabilities π 1,j and π 2,j (j ∈ Z + ) can be expressed in terms of π 0,j by the same formulae as (21) and (22). It is easy to confirm that these results are consistent with those derived in Section 3 where D 1 = 0 due to ν 1 = ν 2 , and also with those presented by Falin [11].
Remark A.1. We notice that the case ν 1 = ν 2 , as it appears in the literature [11], does not distinguish if the service in progress corresponds either to an ingoing call or to an outgoing call. In this short section, we have showed that the results in this paper are helpful to keep knowledge of the identity of the call that is receiving service.
A.2. The case ν 1 = λ + ν 2 . Finally, we consider the case ν 1 = λ + ν 2 . It follows from (7), (13) and the first equality in (15) that the differential equation for Π 0 (z) is given by A comparison with equation (17), for the case ν 1 = λ + ν 2 , shows that the point ρ −1 is now a pole of order 2 in the right hand side of (59) instead of order 1 as in the right hand side of (17). The solution of (59) is given by where we have used (19) is used in the second equality and Remark A.2. We observe that the generating function is expressed in terms of a product of a special hypergeometric function and an exponential function.
The first partial moments for the number of customers in the orbit are given by Remark A.3. It should be noted that Algorithm 1 and Algorithm 2 can be applied for the above special cases ν 1 = ν 2 and ν 1 = λ + ν 2 , where π 0,0 is computed by (14) and (60), respectively.
The transition probability matrix P = (p (i,j,k),(i ,j ,k ) ) of {ζ n ; n ∈ Z + } is given either by p (i,j,k),(i ,j ,k ) = q (i,j,k),(i ,j ,k ) q i,j,k , if (i , j , k ) = (i, j, k) or by 0 if (i , j , k ) = (i, j, k). We observe that sup (i,j,k)∈S q i,j,k = +∞, inf (i,j,k)∈S q i,j,k > λ > 0, which guarantees that {χ(t); t ≥ 0} is ergodic if the embedded chain is ergodic. To apply the Foster's criterion to {ζ n ; n ∈ Z + }, we consider a test function f (i, j, k) = ai + bj + k, where a ≥ 0 and b ≥ 0 will be appropriately determined later. Let φ(i, j, k) denote the mean drift associated with f (i, j, k), so we have follows: where N