Simulation of multiclass retrial system with coupled orbits Evsey Morozov1,2 , Taisia Morozova2 , and Ioannis Dimitriou3 1 Institute of Applied Mathematical Research, Karelian Research Centre of the RAS, Petrozavodsk, Russia 2 Petrozavodsk State University, Petrozavodsk, Russia, tiamorozova@mail.ru 3 University of Patras, Greece Abstract. In this work, we verify by simulation some recent theoretical results describing the dynamics of the the retrial system with coupled orbits. In such a system, retransmission rate of customers blocked in a virtual orbit depends in general on the binary state, busy or idle, of other orbits. We consider a system with N classes of customers, where an ar- riving customer which meets server busy, joins the corresponding orbit depending on the class of customer. The top (oldest) blocked customer makes an attempt to enter server, with class-dependent exponential time between attempts. At that the retrial rate is defined by the current states (busy or idle) of other orbits. To verify theoretical results, we simulate single-server retrial system with 3 classes of customers following indepen- dent Poisson inputs, while service times are class-dependent and have general distributions. In particular, we verify necessary and sufficient stability conditions and focus on the analysis of symmetric model. Nu- merical experiments confirm theoretical analysis. Keywords: Multiclass retrial queues · Stability · Constant retrial rates · Coupled orbit queues · Cognitive network. 1 Introduction This research is devoted to verification by simulation a few performance re- sults obtained in previous work [17,16]. Moreover, we focus on verification of the obtained stability conditions of the system with coupled orbits. In general, sufficient and necessary stability conditions are different, and the study of the ”gap” between these conditions is important in practice. Indeed, as we show, suf- ficient condition is redundant, and possible extension of stability region seems to be useful to increase the efficiency (throughput) of the system. In this work we pay the main attention to the so-called symmetric model, in which all classes of customers have the same parameters. This scenario allows to simplify anal- ysis and detect some properties which turn out to be useful in a more general setting. In simulation, we focus on the symmetric system with three classes of customers which follow independent Poisson inputs. It is worth mentioning that this research complements and develops previous works [17,16]. To study classical retrial queues, we mention the books [12,1], and the survey papers [2,13]. Also the stability analysis of a multi-class retrial queue with con- stant retrial rates, which do not depend on the states of other orbits, has been developed in [4,5]. As to application of the system with coupled orbits, we mention the mod- elling of wireless multiple access systems, in particular, relay-assisted cognitive cooperative wireless systems [18]. Moreover, in the modern cognitive radio [15] there exists a possibility to dynamically adjusts retransmission rates to improve spectrum utilization [6,8,11]. Furthermore, as it has been mentioned in [16], this model is suitable to describe dynamics of cellular networks, in which the trans- mission rate in a particular cell decreases as the number of users in the neighbor- ing cells increase [6]. A similar effect is observed in the processor sharing models [7,14]. This paper is organized as follows. In Section 2 we describe the basic model. In Section 3 we summarize the main theoretical results obtained in [17] which we verify by simulation in Section 4. In particular, we introduce the symmetric model with coupled orbits. In section 4 we verify theoretical results simulating 3- class symmetric model with exponential and Pareto service times. In particular, we estimate the stationary probability that a fixed orbit is busy and server is idle, and demonstrate the correctness of the lower and upper bounds of this probability. Also the accuracy of stability conditions is studied, using the ”gap” between the necessary and sufficient stability conditions. Actually simulation shows that the necessary stability condition is in fact stability criterion of the system. 2 Description of the model We study a single-server with no buffer for the waiting customers, and with three classes of customers. Nevertheless, it is worth mentioning that the theoretical results, which we will verify, hold for an arbitrary number N of customer classes. By this reason we will formulate below theoretical results for this general setting. The class-i customers form Poisson input with rate λi , i = 1, . . . , N . Because all inputs are assumedPto be independent, then the summary input is Poisson as well, with rate λ := i λi , in which an arbitrary arrival belongs to class i with the probability pi =: λi /λ, i = 1, . . . , N. Then interarrival times of the input are exponential with generic element τ and expectation Eτ = 1/λ ∈ (0, ∞). We (i) also assume service times of class-i customers, {Sn , n ≥ 1}, to be independent identically distributed (iid) with service rate 1 γi =: ∈ (0, ∞), i = 1, . . . , N. ES (i) A customer, meeting server busy joins a class-dependent virtual orbit, where joins the end of the orbit queue. At that the head customer waiting in orbit i makes retrial attempts until he finds server idle to occupy it. The distance between attempts are exponentially distributed with rate µJ(i) , where J(i) is the current configuration of the orbits: busy or empty. Thus each orbit acts as a FIFO queueing system with state-dependent ”service” rate µJ(i) . This dependence is a key new property of the model. To be more precise, for each i, we define (N − 1)-dimensional vectors J(i) = {j1 , . . . , ji−1 , ji+1 , . . . , jN } with binary components jk ∈ {0, 1}, where the ith component is omitted. If the k-th orbit is currently busy, we put jk = 1, otherwise, jk = 0. Each vector J(i) is called configuration. For each i, we introduce the set G(i) = {J(i)} of possible configurations. It is assumed that there is given constant µJ(i) , retransmission rate from orbit i, if current configuration is J(i). We denote Mi the set of rates for all configurations belonging to G(i). This construction, proposed in [17], considerably generalizes the setting stud- ied in previous works [16,11,10,9]. In these works, it is assumed that orbit i has rate µi if at least one (other) orbit is busy, otherwise, the rate is µ∗i , i = 1, . . . , N . Thus, in setting in [16], each set Mi = {µ∗i , µi }. In this work we continue to study th new general setting from [17] with focus on verification some bounds and sta- bility conditions proved in [17], for the system with three classes of customers. Before to give main stationary performance measures to be verified by simu- lation, we mention that the main stochastic processes describing the dynamics of the system, such like accumulated work (workload) orbit size, ect., are regenera- tive, with regeneration instants Tn . A regeneration occurs when a new customer meets an idle system [3]. The distances Tn+1 − Tn are iid regeneration periods, and we denote T the generic period. A queuing process is called positive recurrent if the mean generic period is finite, that is ET < ∞ [3]. Under positive recurrence, there exists the stationary regime of the system [3]. 3 Preliminary results Now we define the main stationary performance metrics and give necessary and sufficient stability conditions proved in [17] by regenerative method. Let I(t) be the summary idle time of the server in interval [0, t], then busy time is defined as B(t) = t − I(t). If the system is positive recurrent, then there exist the limits, with probability (w.p.) 1, I(t) lim = P0 = 1 − Pb , t→∞ t where P0 is the stationary idle probability of the server, and B(t) Pb = lim t→∞ t is the stationary busy probability of the server. Denote Bi (t) the summary time, in interval [0, t], when the server is occupied by class-i customers. It is proved in [17] that the stationary probability that the server is occupied by class-i customer is defined as Bi (t) (i) lim = Pb = ρi , i = 1, . . . , N. t→∞ t Denote the traffic intensity for each class, ρi = λi /γi , i = 1, . . . , N, and summary traffic intensity X ρ= ρi . P (i) Because Pb = i Pb then it follows that Pb = ρ. Now we introduce the maximal possible rate from orbit i: µ̂i = max µJ(i) . J(i)∈G(i) The following statement, which has been proved in [17], contains the necessary stability (positive recurrence) condition of our system. Theorem 1. If the N -class retrial system with coupled orbits is positive recur- rent, then N X h µ̂i i Pb = ρi = ρ ≤ min < 1. (1) i=1 1≤i≤N λi + µ̂i To formulate the next statement and sufficient stability conditions, we de- note, for each class i, the minimal retrial rate µ0i = min µJ(i) . J(i)∈G(i) (i) Also let P0 be the stationary probability that server is idle and orbit i is busy. The following statement is proved in [17]. Theorem 2. The following inequalities hold λi (i) λi ρ ≤ P0 ≤ 0 ρ, i = 1, . . . , N. (2) µ̂i µi Now we formulate the sufficient stability condition [17]. Theorem 3. The sufficient stability condition of N -class retrial system with coupled orbits is N X λ ρi + max 0 < 1. 1≤i≤N µi + λ i=1 This condition implies a negative drift of the workload process, and as a result, positive recurrence of the system, and can be written as µ0i X   ρ= ρi < min . (3) i i λ + µ0i To compare (3) with the necessary stability condition (1), we define the gap between the necessary and sufficient conditions, µ̂i µ0i ∆ =: min − min > 0. (4) i λi + µ̂i i λ + µ0i An important special class constitute the symmetric systems. To explain sym- metry for the system with coupled orbits, we first note that this system becomes standard retrial system with constant rate provided µJ(i) ≡ µi , implying equal- ity µ̂i = µ0i = µi . The latter standard retrial system with constant rate, becomes symmetrical, if the corresponding parameters are equal, that is λi ≡ λ, γi ≡ γ, µi ≡ µ. However, the notion symmetrical system with coupled orbits is more flexible. To explain it in more detail, we define, for each i, the set of vectors describing retrial rates for each configuration J(i) = {j1 , . . . , jN }. Because in our case N = 3, then the capacity |J(i)| = 4 for each i = 1, 2, 3. To compose J(i), we use lexicographical order, that is, for each orbit i, and two remaining orbits j < k, with k, j 6= i, the following four configurations J(i) are possible: Mi =: {(ij = 0, ik = 0), (ij = 1, ik = 0), (ij = 0, ik = 1), (ij = 1, ik = 1)}. (5) (Recall that ik = 1 means that orbit k is busy, while ik = 0 means that it is empty.) We denote µi00 , µi10 , µi01 , µi11 , the retrial rates corresponding to each configuration in (5). Then we obtain that, in the symmetrical system, all sets G(i) = {µi00 , µi10 , µi01 , µi11 }, i = 1, 2, 3, are identical, although rates within G(i) may differ. One can expect that this structure leads to similar behavior of the orbits, and it is confirmed below by simulation. Note that for the symmetric coupled orbits, the difference (4) becomes µ̂ µ0 ∆ := − , λ/N + µ̂ λ + µ0 where µ̂ = µ̂i , µ0 = µ0i . Finally, for symmetric classical (non-coupled) orbits, µ̂ = µ0 and (4) becomes µ µ ∆= − . λ/N + µ λ + µ In the next section, containing simulation results, we analyze the symmetric model and leave studying a more general model for a further work. Remark. For non-coupled orbits, µJ(i) = µi for all configurations J(i), and each busy orbit i has a fixed retrial rate µi . Then relation (2) becomes (i) λi Pb = µi P0 , and we obtain the following explicit formula for the stationary probability that orbit queue i is busy and server is idle, see [17]: N (i) λi X P0 = ρk , i = 1, . . . , N. µi k=1 4 Simulation results In this section we verify by simulation some obtained above theoretical results for three classes of customers considering symmetric model. To this end, we define the following variables, µ̂i µ0i Γ1 := min − ρ, Γ2 := min − ρ. i λi + µ̂i i λ + µ0i which delimit the boundary of stability region. More exactly, if Γi > 0, i = 1, 2, then both stability conditions (1) and (3) are satisfied. If Γi < 0, i = 1, 2, then instability of the orbits is expected. However, if the intermediate case Γ1 > 0, Γ2 < 0 holds, then, as simulation shows, the orbits are stable in all experiments. It indicates that the necessary stability condition is in fact stability criterion, while sufficient stability condition is redundant. (However we can not prove it strictly for the system with general service times.) These observations are illustrated by simulation below. We emphasize again that in this work we pay attention simulation symmetric model with coupled orbits. Now we present numerical results obtained by simulation, to verify theoreti- cal analysis of stationary regime and necessary and sufficient conditions (1), (3). Everywhere we use the black, grey and dotted curve to demonstrate the dynam- ics of the 1st, 2nd and 3rd orbit, respectively. (The axis t counts the number of discrete events: arrivals, departures, attempts, in the applied discrete-event simulation algorithm.) First we perform some experiments for the completely symmetric model and demonstrate stability/instability of all orbits. This analysis also shows the re- dundancy of sufficient condition (3) for stability because simulation stays all orbits stable even for Γ1 < 0. In the 1st experiment, we use the following input and service rates, λ1 = λ2 = λ3 = 3, γ1 = γ2 = γ3 = 15, and the following retrial rates: M1 = {µ100 = 20, µ110 = 30, µ101 = 15, µ111 = 25}, M2 = {µ200 = 20, µ210 = 30, µ201 = 15, µ211 = 25}, M3 = {µ300 = 20, µ310 = 30, µ301 = 15, µ311 = 25}. (6) Thus, with those parameters we receive ρ = 0.6 and Γ1 = 0.3, Γ2 = 0. Therefore, both stability conditions are satisfied and all orbits are stable, as expected, and we can see that at Fig. 1 7 6 5 4 N(t) 3 2 1 0 0 500 1000 1500 2000 2500 3000 t Fig. 1. The symmetric system, exponential service time. Condition (1) and (3) hold: Γ1 > 0, Γ2 > 0; all orbits are stable. In the following experiments we show the results for Pareto service times (Pareto model), and not equal service rates, ceteris paribus. That is we still keep a symmetry in the input rates. Fig. 4 shows the dynamics of the orbits in Pareto model with service time distribution xi Fi (x) = 1 − ( 0 )α , x ≥ xi0 (Fi (x) = 0, x ≤ xi0 ), x and expectation α xi0 ES (i) = , α > 1, xi0 > 0, i = 1, 2, 3. α−1 We select α = 2 and the following values of the shape parameter xi0 for orbit i = 1, 2, 3, respectively: 1 xi0 = , i = 1, 2, 3. 24 This choice gives the following service rates γ1 = γ2 = γ3 = 12, ceteris paribus. Here ρ = 0.76, implying Γ1 = 0.14 and Γ2 = −0.16. Thus condition (1) holds while condition (3) is violated. As we see, all three orbits remain stable, however the stability is reached at a higher level. In the 3rd experiment, shown on the Fig.3, we further increase service rate 20 15 N(t) 10 5 0 0 5000 10000 15000 t Fig. 2. The symmetric system, Pareto service time. Condition (1) holds, condition (3) is violated: Γ1 > 0, Γ2 < 0; all orbits are stable at a higher level. (ceteris paribus) γ1 = γ2 = γ3 = 10. Here we obtain ρ = 0.9 and it gives Γ1 = 0, Γ2 = −0.3. Thus both conditions (1) and (3) are violated. (Γ1 = 0 is called boundary case.) As we see on Fig. 3, all orbits become now unstable. Thus in the these experiments, a gradual decreasing service rates (which im- plies reduction Γ1 and Γ2 ) makes orbits unstable only if the necessary condition (1) is violated. So we suggest that condition (3) is redundant, and moreover that the necessary condition (1) is indeed stability criterion. In the following experiment we demonstrate the estimation the stationary (1) probability P0 for the symmetric model with exponential service time. In that experiment we use the following input and service rates, λ1 = λ2 = λ3 = 3γ1 = γ2 = γ3 = 15, while the retrial rates remain (6). Thus in this case µ0i and µ̂i are different, and exact value of the target probability is unknown. However, by the positive recurrence, the sample mean estimate still converges to a limit. Fig. 4 shows that (1) the sample mean estimator of P0 satisfies the corresponding inequality in (2). (2) (3) (The dynamics of the estimators of P0 and P0 is similar and is omitted.) 40 30 N(t) 20 10 0 0 1000 2000 3000 4000 5000 6000 7000 t Fig. 3. The symmetric system, Pareto service time. Conditions (1) and (3) are violated: Γ1 < 0, Γ2 < 0; all orbits are unstable. 5 Conclusion In this work, we simulate a 3 -class symmetric retrial system with independent Poisson inputs and the coupled orbits to verify some theoretical results found earlier. In this system, a new customer meeting server busy joins the correspond- ing infinite capacity orbit. The retrial rate from orbit i depends on the current configuration of other orbits: busy or idle. We verify by simulation some sta- tionary performance measures and the accuracy of the found earlier stability conditions of this model. ACKNOWLEDGEMENTS The study was carried out under state order to the Karelian Research Centre of the Russian Academy of Sciences (Institute of Applied Mathematical Research KRC RAS). The research of EM is partly supported by Russian Foundation for Basic Research, projects 18-07-00147, 18-07-00156. The research of TM is supported by Petrozavodsk State University and Russian Foundation for Basic Research, project 18-07-00147. References 1. Artalejo, J.R., Gómez-Corral, A.: Retrial Queueing Systems: A Computational Approach. Springer-Verlag Berlin Heidelberg (2008), https://doi.org/10.1007/ 978-3-540-78725-9 0.30 0.25 probability of busy orbit 1 and idle server 0.20 0.15 0.10 0.05 0.00 500 1000 1500 2000 2500 3000 3500 t (1) Fig. 4. The symmetric system, Pareto service time. Estimation the probability P0 = P(busy orbit 1, idle server). 2. Artalejo, J.: Accessible bibliography on retrial queues: Progress in 2000-2009. Mathematical and Computer Modelling pp. 9–10 (2010) 3. Asmussen, S.: Applied probability and queues. Springer, New York (2003) 4. Avrachenkov, K., Morozov, E., Nekrasova, R., Steyaert, B.: Stability analy- sis and simulation of N-class retrial system with constant retrial rates and poisson inputs. Asia-Pacific Journal of Operational Research 31(2) (2014). https://doi.org/10.1142/S0217595914400028 5. Avrachenkov, K., Morozov, E., Steyaert, B.: Sufficient stability conditions for multi-class constant retrial rate systems. Queueing Systems 82(1-2), 149–171 (Feb 2016). https://doi.org/10.1007/s11134-015-9463-9, http://link.springer. com/10.1007/s11134-015-9463-9 6. Bonald, T., Borst, S., Hegde, N., Proutiere, A.: Wireless data performance in multi- cell scenarios. Proc. ACM Sigmetrics/Performance ’04 pp. 378–388 (2004) 7. Bonald, T., Massoulié, L., Proutiére, A., Virtamo, J.: A queueing analysis of max- min fairness, proportional fairness and balanced fairness. Queueing Syst. (2006) 8. Borst, S., Jonckheere, M., Leskela, L.: Stability of parallel queueing systems with coupled service rates. Discrete Event Dyn. S. pp. 447–472 (2008) 9. Dimitriou, I.: Modeling and analysis of a relay-assisted cooperative cognitive net- work. Springer (2017) 10. Dimitriou, I.: A queueing system for modeling cooperative wireless networks with coupled relay nodes and synchronized packet arrivals. Perform. Eval. (2017). https://doi.org/10.1016/j.peva.2017.04.002 11. Dimitriou, I.: A two class retrial system with coupled orbit queues. Prob. Engin. Infor. Sc. pp. 139–179 (2017) 12. Falin, J., Templeton, J.G.C.: Retrial Queues. Chapman and Hall/CRC (1997) 13. Kim, J., Kim, B.: A survey of retrial queueing systems. Annals of Operations Research pp. 3–36 (2016) 14. Liu, X., Chong, E., Shroff, N.: A framework for opportunistic scheduling in wireless networks. Comp. Netw. pp. 451–474 (2003) 15. Mitola, J., Maguire, G.: Cognitive radio: making software radios more personal. IEEE Pers. Commun. 6(4) pp. 13–18 (1999) 16. Morozov, E., Dimitriou, I.: Stability analysis of a multiclass retrial system with cou- pled orbit queues. Proceedings of 14th European Workshop, EPEW 2017, Berlin, Germany, September 7-8, 2017 (2017). https://doi.org/10.1007/978-3-319-66583- 2-6 17. Morozov, E., Morozova, T.: Analysis of a generalized system with coupled orbits. Proceedings of Fruct23, Bologna (2018) 18. Sadek, A., Liu, K., Ephremides, A.: Cognitive multiple access via cooperation: Protocol design and performance analysis. IEEE Trans. Infor. Th. 53(10) pp. 3677– 3696 (2007)