On the BMAP_1, BMAP_2/PH/g, c retrial queueing system

In this paper, we analyze a retrial queueing system with Batch Markovian Arrival Processes and two types of customers. The rate of individual repeated attempts from the orbit is modulated according to a Markov Modulated Poisson Process. Using the theory of multi-dimensional asymptotically quasi-Toeplitz Markov chain, we obtain the stability condition and the algorithm for calculating the stationary state distribution of the system. Main performance measures are presented. Furthermore, we investigate some optimization problems. The algorithm for determining the optimal number of guard servers and total servers is elaborated. Finally, this queueing system is applied to the cellular wireless network. Numerical results to illustrate the optimization problems and the impact of retrial on performance measures are provided. We find that the performance measures are mainly affected by the two types of customers' arrivals and service patterns, but the retrial rate plays a less crucial role.


1.
Introduction. Most queueing systems assume that the arrival process is a stationary Poisson process. But such a process does not characterize the typical features of traffic in modern telecommunication networks such as correlation and burstiness. Bursty traffic does appear in telecommunication networks, such as orthogonal frequency division multiplexing networks (see [27] and [24]). Klemm et al. [14] pointed out that the batch Markovian arrival process (BMAP) is a useful and appropriate model for describing such features in Internet-related systems. For example, BMAPs are a better representative of the future traffic characteristics than the traditional Poisson traffic (see [26]). The BMAP is a generalization of many well-known processes including the Markovian arrival process (MAP), the Markov-modulated Poisson process (MMPP), the PH-renewal process and the Poisson process. Furthermore, the BMAP preserves the tractable Markovian structure.
The origins of the BMAP can be traced to the development of the versatile Markovian point process by Neuts [21]. The currently used denotations were introduced in [22]. It is widely employed in research for queueing systems, inventory systems, reliability engineering, manufacturing systems, computer communication systems, insurance problems and so on. Chakravarthy [4] gave a detailed survey of the results about the queues which the input flow is a BMAP.
Retrial queues are very good mathematical models for cellular wireless networks, telephone switching systems, local area networks under the protocols of random multiple access, computer systems for competing to gain service from a central processor unit, etc. Over recent years it has been a rapid growth in the literature on retrial queues. The reader can refer to the survey papers by Falin [7] and Artalejo [1,2] for a review of main results and methods. Most of the papers and books are focused on single server retrial queues. Multi-server retrial queues are more complicated for research and even the model M/M/N is studied until now only algorithmically via different truncation schemes. But for the finite source MAP/PH/N retrial queue, Wu et al. [28] gave a detailed analysis of this queueing system operating in a random environment. The multi-server retrial queues with the BMAP are still not well investigated in literature. In [3], the classical BMAP/PH/N retrial queue was analyzed by means of discrete-time multi-dimensional asymptotically quasi-Toeplitz Markov chain. Recently, Klimenok et al. [15] analyzed the BMAP/PH/N retrial queue with impatient customers. Later, Kim et al. [11] studied the BMAP/PH/N retrial queue with Markovian flow of breakdowns. Subsequently, Kim et al. [12] considered the BMAP/PH/N retrial queueing system operating in a Markovian random environment. Dudin and Klimenok [5] investigated the BMAP/PH/N retrial queue with Markov modulated retrials. Particularly, for the retrial queues with two types of customers, the single server retrial queues have received considerable attention in the literature. For example, using the matrix analytic method, Martin and Artalejo [19] studied an M 1 , M 2 /G 1 , G 2 /1/1 retrial queue with impatient customers. Recently, Wu and Lian [29] introduced negative customers into the retrial queues with two types of customers. Morozov et al. [20] considered a single-server retrial queue with multiple classes of customers. Our model is very close to the model of Kim et al. [13]. They analyzed a MAP/M/N retrial queue with two types of customers. However, our model has essential differences in comparison to that paper: batch arrivals, MMPP retrials, PH-type service. In fact, in the multiple access cellular system, many users make calls independently at the same time. Hence, we use the BMAP to model the arrival process. MMPP is a Poisson process modulated by a random environment. That is to say, the external environment can change the retrial rate which is more realistic. In addition, PH distribution is a generalization of exponential distribution and can approach any general distribution.
We successfully studied this generalized model by means of asymptotically quasi-Toeplitz Markov chains. Here, we briefly introduce the continuous-time multidimensional asymptotically quasi-Toeplitz Markov chains. For a comprehensive and excellent overview of multi-dimensional asymptotically quasi-Toeplitz Markov chains with discrete and continuous time, readers may refer to the paper [16]. The next definition and propositions follow [16]. Let ξ t = {i t , a t }, t ≥ 0, be a regular irreducible continuous time Markov chain with the state space Ω = {(i, a), a ∈ L i , i = 0, 1, · · · , i 0 ; (i, a), a ∈ L, i > i 0 }. Enumerate the states of the chain {ξ t , t ≥ 0} as follows: the states (i, a) are numerated in ascending order of the component i and for fixed i, the states (i, a) are arranged in lexicographic order. We denote the block matrix Q = (Q i,j ) as the generator of the chain, where the block Q i,j is formed by intensities q (i,a);(j,b) of the chain, transition from the state (i, a) to the state (j, b).
Y k z k be the matrix generating function characterizing the limiting Markov chain. The following two propositions present the ergodicity conditions for the Markov chain {ξ t , t ≥ 0} in terms of the generating function Y (z) based on whether the matrix Y (1) is an irreducible matrix or not.
is an irreducible matrix and assume that (i) the series ∞ n=1 nQ i,i+n−1 e converges for i = 0, i 0 , (ii) the series ∞ n=1 nQ i,i+n−1 e converges for i > i 0 , and there exists a positive integer i * > i 0 such that this series converges uniformly in the region i ≥ i * , (iii) the sequence Λ −1 i , i ≥ 0, has an upper bound, then the sufficient condition for the Markov chain {ξ t , t ≥ 0} ergodicity is the fulfillment of the inequality If Y (1) is a reducible matrix, then it can be presented in a normal form which has the following structure: is a reducible matrix of the above form and the conditions of (i)-(iii) of Proposition 1 are satisfied, then the sufficient condition for the Markov chain {ξ t , t ≥ 0} ergodicity is the fulfillment of the inequality The rest of the paper is organized as follows. The model description is given in Section 2. The Markov chain associated with the system is analyzed in Section 3 where we formulate the model as a multi-dimensional asymptotically quasi-Toeplitz Markov chain. The algorithms for calculating the stationary state probabilities are elaborated in Section 4. A set of performance measures is demonstrated in Section 5. Two optimization problems and an algorithm for calculating the optimal values are provided in Section 6. An application for the model under discussion and some numerical examples are shown in Section 7. Finally, a summary of the results is presented in Conclusion.
2. Model description. We consider a multi-server retrial queueing system. Two types of customers (ordinary customers and priority customers) arrive to the system according to two independent BMAPs (denoted by BM AP 1 and BM AP 2 ). The notion of the BMAP and its detailed description were given by Lucantoni [18]. We denote the directing process of the BM AP 1 by {ν t , t ≥ 0} with the state space {1, · · · , W } and the directing process of the BM AP 2 by {γ t , t ≥ 0} with the state space {1, · · · , V }. The two directing processes behave as two independent irreducible continuous-time Markov chains. Denote the matrix sequences (D n : n ∈ N 0 ) and (E n : n ∈ N 0 ) as the sequences of characterizing matrices for the BM AP 1 and BM AP 2 , respectively. Introduce the matrix generating functions D(z) = ∞ n=0 D n z n and E(z) = ∞ n=0 E n z n , |z| < 1. We assume the matrix functions D(z) and E(z) satisfy all assumptions of [18]. Therefore, the matrix D(1) is the generator of the process {ν t , t ≥ 0} and the matrix E(1) is the generator of the process {γ t , t ≥ 0}. Suppose that the two arrival processes BM AP 1 and BM AP 2 start with the initial phase probability distributions α 1 and α 2 , respectively. Let θ 1 and θ 2 be the stationary probability vectors of D(1) and E(1), respectively. Then λ 1 = θ 1 D (1)e is the stationary arrival rate of the ordinary customers and λ 2 = θ 2 E (1)e is the stationary arrival rate of the priority customers. Moreover, we define λ b1 = θ 1 (−D 0 )e and λ b2 = θ 2 (−E 0 )e as the intensities of group ordinary arrivals and group priority arrivals, respectively. Here e is a column vector of appropriate size with all elements equal to 1.
The service facility consists of c (c ≥ 2) identical servers which operate independently of each other. All customers have the same independent, identical PH distributed service times S which is governed by the irreducible continuous time Markov chain {m t , t ≥ 0} with the state space {1, · · · , M }. The transitions of the Markov chain {m t , t ≥ 0}, which do not lead to service completion, are defined by the irreducible matrix S of size M × M . The transitions of the Markov chain {m t , t ≥ 0}, which lead to service completion, are defined by the vector S 0 = −Se. At the service beginning epoch, the state of the process {m t , t ≥ 0} is determined by the probabilistic row vector ς of size 1 × M . The service time distribution function has the form B(x) = 1 − ςe Sx e. The mean rate of service is µ = [−ςS −1 e] −1 . A more detailed information about the PH type distribution can be seen in the books [23] and [17].
We assume there are c − g (1 ≤ g ≤ c − 1) servers (called guard servers) opening only for the priority customers. When a batch of priority customers arrives at the system and there are several servers being idle, the priority customers occupy the corresponding number of the servers. If the number of the idle servers is insufficient (or all the servers are busy) the rest of the batch (or all the batch) enters the retrial group (called orbit). When a batch of ordinary customers arrives at the system and meets no more than g − 1 busy servers, the ordinary customers occupy the corresponding number of the servers (If the number of the idle servers is insufficient, the rest of the batch goes to the orbit); otherwise, the whole batch of ordinary customers is blocked and goes to the orbit. These customers in orbit are said to be repeated customers and the orbit capacity is assumed to be unlimited. For all customers in the orbit, the individual repeated attempts are governed by a common MMPP. The directing process of the MMPP is denoted by {r t , t ≥ 0} which is an irreducible continuous time Markov chain with the state space {1, · · · , R}. Transition intensities of the Markov chain {r t , t ≥ 0}, which are not related with customers retrials, are defined by the matrix T 0 , and transition intensities, which are accompanied by retrial, are described by the diagonal matrix T 1 = diag{σ 1 , · · · , σ R }. The matrices T = T 0 + T 1 is the infinitesimal generators of the processes {r t , t ≥ 0}. That is, when the process r t stays in the state r, r = 1, · · · , R, each customer from the orbit makes attempts to seek the service at exponentially distributed time intervals with mean 1/σ r . Hence, the total rate of retrials is equal to nσ r , σ r > 0, when the process r t is in the state r and the orbit size (the number of calls on the orbit) is equal to n, n > 0. Let θ 0 be the stationary probability vector of T . Then σ = θ 0 T 1 e is the stationary retrial rate of the repeated calls in the orbit. On retrial, an orbiting customer obtains service immediately from one of the idle servers if the number of busy servers is no more than g − 1; otherwise it will return to the orbit for a later retrial.
For the use in the sequel, let us introduce the following notations: e is a column vector of ones of suitable size. When needed we will identify the dimension of this vector with a subindex; I is an identity matrix of appropriate size. When needed we will identify the dimension of this matrix with a subindex; O L is a zero matrix of size L × L; diag{a r , r = 1, L} is a diagonal matrix with diagonal entries a r , r = 1, L; ⊗ and ⊕ are the symbols of Kronecker product and sum of matrices, see, e.g. [8]; I n m ⊗ Ω ⊗ I n l−m−1 , l ≥ 1, for the matrix Ω having n rows.
Our work objectives are to calculate the stationary state distribution, to derive the main performance measures of the system and optimize the values of g and c.
3. Analysis. In this section, we adopt a layered approach for the description of the state space to construct the infinitesimal generator Q of the Markov chain which can represent the model. Let * o t be the number of repeated customers in the orbit, at the epoch t, t ≥ 0, o t ≥ 0; * b t be the number of busy servers, at the epoch t, t ≥ 0, b t = 0, c; * r t be the state of the directing process of the MMPP of repeated customers, at the epoch t, t ≥ 0, r t = 1, R; * ν t be the state of the directing process of the BMAP of ordinary customers, at the epoch t, t ≥ 0, ν t = 1, W ; * γ t be the state of the directing process of the BMAP of priority customers, at the epoch be the state of the directing process of the service on the jth busy server, at Here, we assume that the busy servers are numerated in order of their occupying, i.e., the server, which begins the service, is appointed the maximal number among all busy servers; when some server finishes the service and is released, the servers are correspondingly enumerated.
Then, the process of the model can be described by the regular irreducible continuous time multi-dimensional Markov chain: and the dimension of the state space of the Markov chain {ξ t , t ≥ 0} is equal to We note that the dimension K can be very large. E.g., if the array (c, R, W, V, M ) is equal to (10, 2, 2, 2, 2), then K = 16376.
We suppose that the states of the Markov chain {ξ t , t ≥ 0} are enumerated in lexicographic order. By means of considering the probabilities of the Markov chain transitions during an infinitesimal time interval, the infinitesimal generator matrix Q of the Markov chain {ξ t , t ≥ 0} has the following structure: where the blocks Q ij of size K × K are defined by: Remark 1. When E n = 0, n = 0, 1, · · · and c = g, the present system reduces to a retrial BMAP/PH/c queueing system with Markov modulated retrials. It may be noted that the results above after putting E n = 0, n = 0, 1, · · · and c = g agree with the results presented in [5].
In the following sections, we analyze the queueing system by using the theory of multi-dimensional asymptotically quasi-Toeplitz Markov chains.
The diagonal blocks of the generator Q are given by Denote the diagonal entries of the matrices D 0 , E 0 and S as {−λ Then the diagonal matrix Λ i in Definition 1.1 is given by: Further, introduce the matricesĪ, Υ and Γ k , k ≥ 0, of size K × K: and ∞ k=0 Y k e = e.
It can be easily verified that the Markov chain {ξ t , t ≥ 0} satisfies Definition 1.1. Hence, it belongs to the class of asymptotically quasi-Toeplitz Markov chain (see [3] and [16]). From (1)-(3), we get where Θ l (z), l = l, l = g, c, I RW ⊗ E r ⊗ I M l ⊗ ς ⊗r , l = l + r, r = g + 1 − l, c − 1 − l, l = 0, g − 1, Substituting the expressions of the matrices Υ,Ī, ∆ and Γ(z) in to (4), we see that the matrix function Y (z) is of the form It is seen from (5) that the matrix Y (z) is reducible and the matrix Y (1) contains only one irreducible stochastic diagonal block Y 22 (1). According to the Proposition 2, a sufficient condition for ergodicity of the Markov chain {ξ t , t ≥ 0} is the fulfillment of the following inequality: Next, we assume that the system is stable. Denote the steady state probabilities of the Markov chain as Let the row-vector P i of dimension K denote the stationary-state probabilities P i = P i, b, ν, γ, m (1) , · · · , m (b) , i ≥ 0, and define the infinite-dimensional probability vector P = (P 1 , P 2 , · · · ). When the system is stable, P is the unique solution to PQ = 0 and Pe = 1. However, it is still an open problem to express P in the closed form (e.g., generating function), and even though Q is highly structured, P cannot be expressed in a tractable analytical form. In the following, we will adopt the algorithm of asymptotically quasi-Toeplitz Markov chains to solve the equilibrium equation PQ = 0.

Algorithm.
In this section, we use the numerically stable algorithm which has been elaborated in [16] to calculate the stationary distribution of the system. The algorithm is based on censoring technique of the asymptotically quasi-Toeplitz Markov chain and consists of some essential steps which are given in the following algorithm. Step 1: Calculate the matrices G k from the following recursive equation Step 2: Calculate the matricesH i,j , j ≥ i, i ≥ 0,using the formulaē Step 3: Calculate the matrices F j , j ≥ 0, according to the recursive relation Step 4: Calculate the vector P 0 as the unique solution to the system of the linear equations Step 5: Calculate the vectors P j , j ≥ 1, as

Remark 2. We can eliminate the infinite sums in
Step 4 as the following method: pick N > k 0 , and replace the second equation in Step 4 with P 0 N j=0 F j e = 1, such that the inequality P N e < 0 is satisfied, where 0 is some preassigned small number. Note that the inverse matrices that appear in the algorithm exist and are nonnegative. So the computation of the algorithm is stable. The algorithm can be realized in framework of extension of "SIRIUS++" software (see [6]) that was developed using the results of the asymptotically quasi-Toeplitz Markov chains. As such the determination of the stationary distribution is straightforward. 5. Performance measures. Having the stationary distribution P i , i ≥ 0 been calculated, we are able to calculate some stationary performance measures of the system under consideration in this section.
Let Denote P (i, b) the joint probability to have i repeated customers in the orbit and b customers on the servers at arbitrary time. From the results of the above section, we get  1. The probability that there are i repeated customers in the orbit at arbitrary time 2. The probability that there are b busy servers at arbitrary time

The mean number of customers (ordinary customers and priority customers) in the system
6. The blocking probability for an arbitrary ordinary customer 7. The blocking probability for an arbitrary batch of ordinary customers

The blocking probability for an arbitrary priority customer
9. The blocking probability for an arbitrary batch of priority customer Moreover, using the approach based on Markov renewal processes, we can obtain the following statement.
Proposition 3. The generating function Π(z) of the queue state stationary distribution at the departure epochs is characterized by the following functional-differential equation: Proof. The proof can be implemented similarly as the proof of [3] and is omitted here.
6. Optimization. In this section, we borrow some ideas and methods from [25] to investigate the optimization problem that how to choose the optimal values of g and c so as to minimize the priority customer blocking probability as well as the ordinary customer blocking probability. Therefore, we obtain a multi-objective optimization problem in which the number of total servers c and the parameter g which is equal to the number of non-guard servers are the decision variables. However, due to the complexity, we may fix c and only consider g as the decision variable to simply the optimization problem. In the case of the two objectives, we can choose either P b1 or P b2 as the objective function to be minimized and impose a constraint on the other one. If we regard the two blocking probabilities as two binary functions of g and c two variables, then we write two representative optimization problems as follows: (I) Given T k , D k , E k , k = 0, 1, · · · , S and c determine the optimal integer values of g such that Mininize P b1 (g) subject to P b2 (g) ≤ p 0 .
(II) Given T k , D k , E k , k = 0, 1, · · · , S determine the optimal integer values of g and c such that Mininize c subject to P b1 (g, c) ≤ p 1 , and P b2 (g, c) ≤ p 2 .
Here, the constants p 0 , p 1 and p 2 are three preassigned numbers. We first consider the optimization problem (I). It is clear that P b2 (g) > P b2 (g −1) (see Fig. 3). Thus, we can get the largest value of g such that P b2 (g) ≤ p 0 . Then based on the monotonicity property P b1 (g) > P b1 (g + 1), we know that such a value of g will minimize P b1 (g). Therefore, we can obtain the optimal value of g using a simple one dimensional search over the range {1, · · · , c − 1} for g such that g * = max{g|P b2 (g) ≤ p 0 }.
Then we study the optimization problem (II). Since we cannot derive the analytical solution of this system, in the next section, we will resort to the numerical examples to illustrate the optimization problem (II). Here, we provide an algorithm to obtain the minimum value of c. Algorithm 6.1. Solving the optimization problem (II) Step 1: Set the values p 1 and p 2 and let c=2; Step 2: For g = 1 to c − 1, find the feasible region for g: Step 3: If Θ = ∅, then stop and obtain c * = c. Otherwise, put c = c + 1 and goto step 2; Step 4: Obtain c * = c.

Application to cellular wireless networks and numerical examples.
In this section, we present an application to cellular wireless networks for the model under discussion and some numerical results that show the effect of the system parameters on the performance measures and channel allocation schemes in the cellular mobile networks. In the cellular mobile wireless network, the regional service area is divided into multiple adjacent cells, each of them served by a base station (BS) with a limited number of channels. Here, we assume the number of channels is equal to c, c ≥ 2. A mobile subscriber who has wireless terminal such as mobile phone for voice service, wireless data terminal for data service, dual terminal for voice and data, and so on, communicates via radio links to a BS, one for each cell. Generally, there are two types of calls in a specific cell in the cellular mobile networks: originating calls, i.e., calls originating in that cell, and handoff calls, i.e., current ongoing calls caused by mobile users moving from one cell to that cell. We assume that the originating calls arrive at the system according to a BMAP with the characterizing matrix sequence (D n : n ∈ N 0 ) and the handoff calls arrive at the system according to a BMAP with the characterizing matrix sequence (E n : n ∈ N 0 ). When a batch of new calls originates in that cell, some idle channels are required to the BS for the new calls setup. If there are idle channels in the BS of that cell, then the BS assigns some idle channels to the originating calls; otherwise some new calls (maybe all the batch) in that batch are blocked and the blocked mobile subscribers may try their luck later. Similarly, if some idle channels are available in the BS of that cell, the batch of handoff calls is successfully handed over without an interruption; otherwise some handoff calls (or all the batch) are blocked and retry for service after a short period. Any blocked call (originating or handoff) retries for service where the rates of the retrials depend on a changing environment. Thus, we use the MMPP with the characterizing matrices (T 0 , T 1 ) to model the retrial process. The quality of service is disturbed seriously if a handoff call becomes blocked when it crosses the cell boundary. Generally, the base station gives priority to a handoff call over an originating call. The usual way to do this is the scheme with guard channels, see [9] and [10]. Here, we assume the number of non-guard channels is equal to g. Thus, the number of guard channels is equal to c−g. Furthermore, we suppose the conversation time of an ongoing call (originating or handoff) by a channel is of PH type with a representation (ς, S). Then the cellular mobile wireless network can be modeled as a BM AP 1 , BM AP 2 /P H/g, c retrial queueing system where the rate of individual repeated attempts from the orbit is modulated according to a MMPP. The originating calls can be regarded as ordinary customers and the handoff calls can be regarded as priority customers.
In We suppose that handoff calls and originating calls arrive to the system according to two independent special BMAPs which are characterized by the following matrices: This means the stationary arrival rates of the originating calls and the handoff calls are λ 1 = 10.6364λ o and λ 2 = 1.6667λ h , respectively. Moreover, we assume the repeated calls in the orbit repeat their attempts to reach a server according to an exponential distribution and the mean retrial rate for the mobile engaged in the cell is σ = 13.2857λ r . The service time is also exponentially distributed with rate µ = 8.1288.
The joint distribution of the number of busy channels and the number of repeated calls in the orbit is presented in Tab. 1 with (c, g, λ o , λ h ) = (8, 6, 2, 2). From Tab. 1, we observe that the joint probability that there are three busy channels and zero repeated call at arbitrary time is equal to 0.2148 which is the maximum value and the probability that there is no repeated call in the orbit is equal to 0.9084 which is the maximum value of the marginal probabilities (lateral). The probability that there are three busy channels at arbitrary time is equal to 0.2206 which is the maximum value of the marginal probabilities (longitudinal).
The effects of the retrial rate and the arrival rate of the originating calls on the optimal value of g for the optimization problem (I) with the set of parameters (c, λ h , θ 0 ) = (20, 10, 0.0001) are reported in Tab. 2. We observe that the optimal value g * decreases monotonously as the arrival rate of the handoff calls increases. From this, we can predict that the optimal value g * increases monotonously as the arrival rate of the originating calls increases. However, the retrial rate of the         repeated calls in the orbit has no significant impact on the optimal value g * . The appearance of this phenomenon may be based on the fact that the mean number of repeated calls is very small compared to the mean number of arriving calls. Hence it may be concluded that the arrival rates of the two types of calls are the major factors in deciding the optimal numbers of guard channels. Tab. 3 shows the effects of the arrival rates of the originating calls and the handoff calls on the optimal value c * for the optimization problem (II), where we set (λ r , p 1 , p 2 ) = (10, 0.001, 0.0001). As is to be expected, when the arrival rates of the two types of calls increase, the optimal value c * increases. Moreover, we observe that the optimal value c * increases faster with the arrival rates of the originating calls than that with the arrival rates of the handoff calls. Therefore, we understand that the arrival of originating calls plays a most significant role in deciding the optimal value c * . In other word, this type of calls whose arrival rate is largest will most strongly affect the optimal design of the system. Fig. 1 illustrates the dependence of the mean number of repeated calls in the orbit on the value of parameter g where we set (c, λ o , λ h ) = (10, 2, 2). As it can be seen from Fig. 1 that L orb decreases monotonously as the value g increases, especially when g ≤ 6. But when g ≥ 7, the curves of the values L orb are almost overlapping with x-axis. That is to say, there are very few repeated customers in the orbit when the number of guard channels is small. Figs. 2 and 3 show the dependence of the blocking probability for an arbitrary originating call and handoff call on the value of parameter g respectively when (c, λ o , λ h ) = (10, 2, 2). From Figs. 2 and 3, we find that the blocking probability for an arbitrary originating call P b1 decreases monotonously, but the blocking probability for an arbitrary handoff call P b2 increases monotonously as the value g increases, which agrees with the intuitive expectations. In addition, when g ≤ 6, the curves of the values P b2 are almost overlapping with x-axis. It means that the blocking probability for an arbitrary handoff call P b2 is roughly equal to zero when the number of guard channels is large.
We must note that, in each of the pictures, the three curves are almost overlapping. This phenomenon tells us that the retrial rate σ has little impact on the system performance measures. From what has been described above, we can draw the conclusion that the performance measures are mainly affected by the two types of calls' arrivals and service patterns, but the retrial rate plays a less crucial role. 8. Conclusions. In this paper, we have investigated the BMAP/PH/N type retrial queue with two types of customers which is a generalization of the MAP/M/N retrial queue with two types of customers studied in [13]. The behavior of this model is described by a multi-dimensional continuous-time asymptotically quasi-Toeplitz Markov chain. An algorithm for computing the stationary distribution is presented. Expressions for calculation of main performance measures of the model are derived. The optimization problem how to choose the optimal values of g and c is discussed and an algorithm is also provided. An application of the model to the cellular wireless network is implemented. The dependences of the main performance characteristics of the model on its parameters are graphically demonstrated. We obtain the results that the retrial pattern has very little impact on the optimal values of guard channels and total channels, the number of busy channels and the blocking probabilities for the two types of calls. Future work will consider the BMAP/PH/N retrial queueing systems with impatient calls and their applications to computer communication networks.