Classical Regeneration Based Approaches for Output Analysis of Multiserver Systems Simulation Irina Peshkova Petrozavodsk State University, Lenina Pr. 33, Petrozavodsk, Russia iaminova@petrsu.ru http://www.petrsu.ru Abstract. We present and discuss three types of regeneration that are based on classical regeneration for steady-state simulation and output analysis of multiserver systems. First method is artificial regeneration based on exponential splitting property of some distributions. Regen- erative envelopes allow to construct upper and lower bound for QoS characteristic of the system. Resampling of New Better than Used inter- arrivals or New Worse than Used service times accelerates simulation. All this techniques are reduced to classical regeneration and allow to receive shorter regeneration cycles. Keywords: Multiserver Systems · Regenerative Simulation· Artificial Regeneration · Splitting · Regenerative Envelopes · NBU distributions. 1 Introduction Modern computer systems, as a rule, have a multiserver or multicore architec- ture, which is caused both by the physical features of the modern technological process of microchip production and by the need to create high-performance computers suitable for solving computationally complex problems. Data trans- fer to the users in modern networks also requires coordinated performance of several devices, as well as work with several servers, for solving problems of re- ducing delays and optimizing throughput. By this reason, reliable estimation of Quality of Service (QoS) parameters of multiserver systems is an important problem to study. Regeneration of stochastic process is one of the powerful and efficient meth- ods of steady-state simulation and output analysis of complex queueing sys- tems [3,5,8,7,6,11]. Standard classical regeneration method uses visits of the process (describ- ing the systems behavior) to a fixed return state as regeneration points. The trajectory between such points can be divided into independent (for classical regeneration) and identically distributed groups (in general, they can be weakly dependent for wide-sense regeneration). It allows to apply well-developed classi- cal statistical methods based on the special form of the Central Limit Theorem (CLT) for confidence estimation of QoS parameters. However, in multiserver systems standard regeneration is not guaranteed by stability condition of the system and regeneration points may be too rare to be useful in practice or even they might not exist. We study an alternative constructive method of obtaining the regeneration points. The so-called artificial regeneration uses the splitting property [2] and can be applied for estimation of multiserver systems performance in two ways. First method is based on constructing a new, equivalent to the original, system by using the splitting property, which allows to change the internal structure of continuous valued random variables in such a way to obtain a sequence of times at which the distribution has some specific properties, and these times occur with positive probability (possibly by enriching the probability space). Second technique uses coupling method and the epochs when the process has a given splitting distribution as regenerative times [1]. If the splitting is based on expo- nential distribution, then the system at this particular points has memoryless property, and thus, regenerates in classical sense in both cases. The regenerative envelopes method is reduced to replacing the original system based on the coupling method [9] and the monotonicity properties of the studied processes with a pair of minorant and majorant systems that regenerate in the standard classical sense [10]. Then we can use upper (lower) bounds of confidence intervals for QoS characteristic of majorant (minorant) system for the confidence interval construction of the original system. We discuss the possibility of accelerating the simulation for multiserver sys- tems in which service times and iterarrival times have New Worse than Used distribution or New Better than Used distribution. If service times have New Worse than Used distribution we resample service time each time when the re- maining service time exceeds a fixed constant a. If interarrival times has New Better than Used distribution we resample input each time when the remaining interarrival time is less than a fixed constant b. In both cases we obtain more frequent classical regeneration that allows to accelerate simulation. 2 Artificial regeneration based on splitting property in multiserver systems The density f of continuous random variable (r.v.) T can be split if there exists some 0 < p < 1 and density f0 , such that: f > pf0 . (1) We define a new Bernoulli random variable I called splitting indicator such that P(I = 1) = p. We construct a new splitted r.v. T 0 (stochastically equivalent to the r.v. T ) as follows: T 0 = IT0 + (1 − I)T1 , (2) where T0 has density f0 , and T1 has density f − pf0 f1 = . (3) 1−p Now we can replace the original r.v. T with a triple [T0 , T1 , I]. The particular case of splitting, the exponential splitting, could be efficiently applied in simulation due to the memoryless property of exponential distribution. We say that a positive valued r.v. T with density f is exponentially splitted if there exist constants λ > 0, τ0 > 0 and 0 < p < 1 such that f (x) > pλe−λ(x−τ0 ) , x > τ0 . (4) That is, splitting representation T 0 of exponentially splitted r.v. T is a triple (T0 , T1 , I) where T0 is left truncated exponentially distributed r.v. with trunca- tion point τ0 , rate λ, density  0, x 6 τ0 , f0 (x) = λe−λ(x−τ0 ) , x > τ0 , and T1 has density ( f (x) 1−p , x 6 τ0 , f1 (x) = f (x)−pλe−λ(x−τ0 ) 1−p , x > τ0 . Now we describe the application of exponential splitting technique to discrete event simulation of multiserver systems. Consider a stochastic process Θ = {X(t), T(t)}t>0 , with discrete (vector) component X = (X1 , . . . , Xn ) and continuous component T = (T1 , . . . , Tm ) > 0 (hereafter we denote vector-valued variables in bold and assume the dimension is clear from context, if not sub-indexed explicitly). Moreover, we assume that T decreases linearly with time, that is, T(t + δ) = T(t) − δ1, for δ > 0, and define an event as a time epoch ti such that Tj ∗ (ti ) = 0 for some j ∗ ∈ {1, . . . , m} (in this case we say that the ith event occurring at ti is of j ∗ th type). We assume that the discrete component X is changed only at event epochs, e.g. by some recurrent relation X(ti +) = Gj ∗ (X(ti −)) , (5) where Gj ∗ is the recursion related to type j ∗ event, or, in general, assume the change is governed by some Markovian probability Pj ∗ (x, x0 ) = Pj ∗ {X(ti +) = x0 |X(ti −) = x}. The component X(t) does not change if T(t) > 0, that formally means X(t) = X(ti +), i = max{k : tk < t}. We also assume that in the vector T the zeroed at ti − component Tj ∗ is initial- ized at ti + from some distribution with density fj ∗ , that is, fj ∗ (u, x, x0 ) = P(Tj ∗ (ti +) ∈ du|X(ti −) = x, X(ti +) = x0 ). (6) Now we assume, that if X(t) = x∗ for some specific x∗ , then the process Θ allows to perform exponential splitting of the continuous component T, that is for each continuous component Tj there exist constants τ0 (j), λ(j), p(j) such that (4) holds for the density fj (u, x, x∗ ), with corresponding density f0,j (u) = λj e−λj (u−τ0 (j)) and f1,j defined appropriately. At that, we introduce a split pro- cess Θ0 = {X(t), T0 (t), Z(t)}t>0 enhanced with auxiliary phase variable Z ∈ {0, 1, 2}m denoting the phases of components of T defined componentwise in a recursive manner as follows if X(ti +) = x∗ and Ij ∗ (ti +) = 1,  1, Zj ∗ (ti +) = (7) 2, otherwise, Zk (ti +) = Zk (ti −), k 6= j ∗ , (8) where ti is the epoch of type j ∗ event, and Ij ∗ is the splitting indicator corre- sponding to the component Tj0∗ with success probability p(j ∗ ). However, starting from Zj ∗ (ti +) = 1 at time ti , the component Zj ∗ makes a step between the events: 1, ti < t < ti + τ0 (j ∗ )  Zj ∗ (t) = (9) 0, t > ti + τ0 (j ∗ ), that is, the phase Zj ∗ is changed from one to zero when not less than τ0 (j ∗ ) time passes since the epoch ti corresponding to the last event of j ∗ th type (we call this time epoch a pseudo-event). In case Zj 6= 1, the component Zj does not change during inter-event time. We define the artificial regeneration epochs as the (increasing) sequence of pseudo-event epochs {βk }k>1 βk+1 = min {t > βk : X(t+) = x∗ , Z(t+) = 0} . (10) This method of artificial regeneration could be applied to performance es- timation of m-server system with energy efficiency control and single waiting room. The discrete and continuous components of the simulated process Θ we define in the following way. The discrete components τ0 (t) is the queue size at time t, and Xi (t) is the mode of ith server, i = 1, . . . , m, X1 (t) = 1 if the server is active, and X1 (t) = 0 if the server is in the sleep mode. Continuous components: T0 (t) is the time before the next arrival epoch, and Ti (t) is the remaining times of current activity (service/sleep) of ith server, i = 1, . . . , m. For the regeneration epoch, we set x∗ = {k, 1, . . . , 1} for some k > 0 and recall that artificial regeneration occurs at such pseudo-event epoch t that X(t+) = x∗ and all the exponentially split components T0 , . . . , Tm are in exponential phases, that is, Z(t+) = 0. The second approach of defining the regeneration times is to construct stochas- tically equivalent system with the same input, but sampling from either f0 or f1 at every random variable initialization time epoch. Then the original r.v. T is re- placed by a stochastically equivalent splitting representation T 0 = IT0 +(1−I)T1 . Coupling method allows to construct new artificial system, which is equivalent to the original system, and by this reason has equivalent QoS parameters. To construct the artificial regeneration points we suggest the following method. Sample the sequence of splitting indicators Ii , i > 1 as a Bernoulli 0-1 r.v. with success probability p. Sample service times T̃j either from density f0 , or from f1 , using relation (2). Fix integer ν0 (fixed number of customers in system) and define n βn+1 = inf k > βn : νk = ν0 6 m, o Ii = 1, i ∈ Mk , Mβn ∩ Mk = ∅ , n > 0, (11) where Mn = {i : ti 6 tn < tdi } is the set of the customers presented in new system at arrival instant tn and tdi is the departure instant of the ith customer. Note that at the constructed regeneration epoch, the remaining work possesses memoryless property, since all the tasks being served have exponentially dis- tributed service times. It is easy to see that Θn := {νn , Ij , j ∈ Mβn }, n > 1 is a Markov process, and that distribution of Θβn is independent of n > 1 and pre-history {Θk , k < βn }. 3 Regenerative envelopes for multiserver systems simulation Now we consider an infinite buffer FCFS m-server GI/G/m queueing system Σ, with the renewal input with instants tn , the iid interarrival times Ln = tn+1 − tn and the iid service times Sn , n > 1. Denote νn the number of customers at instant t− − n , Qn the number of customers waiting in the queue at instant tn , so Qn = max(0, νn − m). The main idea of the regenerative envelopes method is the construction of two new systems: a majorant system Σ and a minorant system Σ with the same input as in the original system Σ, but (stochastically) smaller (in minorant system) or larger (in majorant system) service times. We describe the construction of m-server majorant system Σ̃. The corre- sponding variables in the system Σ̃ we endow with tildes. Assume that the service time distributions in both systems are ordered as FS (x) 6 FS̃ (x), x > 0, that is S >st S̃ (stochastically). Moreover, by construction, L =st L̃. Then the coupling method allows to take t̃n = tn and Sn > S̃n w.p.1. Define the set Mn = {i : ti 6 tn < z i } of the customers, which are being served in the system Σ at instant tn , S i (n) – the remaining service time of customer i at instant t−n in Σ. Fix arbitrary integer ν0 > 0 (number of customers in the system), constants 0 6 a 6 b < ∞ and define n \ o β n+1 = inf k > β n : ν k = ν0 , S i (k) ∈ (a, b), Mβ n Mk = ∅, i ∈ Mk . (12) Then, at the (discrete) instant k = β n+1 , we replace all the non-zero remain- ing times S i (k), i ∈ Mk , by the upper bound b. It is easy to see that Θn := {ν n , S i (n), i ∈ Mn }, n > 1, is a Markov pro- cess and the distribution of Θβ n is independent of n and pre-history {Θk , k < β n }, n > 1. The process {Θn , n > 1} classically regenerates at the instants {β k } and has the i.i.d regeneration cycles {Θk , β n 6 k < β n+1 } with i.i.d cycle lengths αn = β n+1 − β n , n > 1. In a minorant system Σ we define the regenerative moments β n , n > 0 in a similar way (with possibly other values a, b, ν0 ). At each such instant β n we replace the remaining times S i (k), i ∈ Mk , by the lower bound a. The systems Σ, Σ, Σ have zero initial state, and moreover (in an evident notation) Ln = Ln = Ln , S n > Sn > Sn , n > 1. Then it follows from coupling property that ν n > νn > ν n , Qn > Qn > Qn , n > 1, (13) and it allows us to use regenerative simulation of Σ and Σ to obtain an upper and lower bounds of the mean stationary queue size (number of customers) in original system. EQ 6 EQ 6 EQ. 4 Accelerated simulation of multiserver networks with New Better than Used distributions In this section we will discuss how the properties of New Better (Worse) than Used distributions of input or service times could be applied to accelerate regen- eration simulation. We say that the distribution F of r.v. T has increasing (decreasing) failure rate IFR (DFR) if F (t + x) − F (t) Ft (x) = P (T 6 t + x|T > t) = P (T 6 x) = (14) F (t) is increasing (decreasing) in t. Remark, that Ft (x) f (t) lim = = r(t). (15) x→0 x F (t) Rt − r(x)dx Function r(t) is called failure rate and it is easy to see that F (t) = e o . Weibull and Gamma distributions with shape parameter α > 1 are IFR, with α < 1 are DFR. It is known [4], that if ET = µ1 and F is IFR (DFR), then −t F (t) > (6) e µ1 , t < µ1 . (16) We say that the distribution F of r.v. T has increasing (decreasing) failure rate in average IFRA (DFRA) if Zt 1 r(u)du increases in t. (17) t 0 Distribution F of r.v. T is called New Better (Worse) than Used NBU (NWU) if F t 6 F or F (t + x) 6 F (t) · F (x). (18) The following relation between IFR, IFRA and NBU holds. IFR ⇒ IFRA ⇒ NBU. Let S be the (generic) service time (with distribution B), L – (generic) in- terarrival time (with distribution A) in m-server system Σ. Define r.v. St - remaining service time at t− and B t (x) = P (S > t + x|S > t) = P (St > x), r.v. Lt - remaining interarival time at t− and At = P (L > t + x|L > t) = P (Lt > x). If service times S have NBU distribution or interarrival times L have NWU distribution in the original system Σ then we can apply coupling method to construct a new system Σ̃ using property (14) of such distributions. If S has NWU (NBU) distribution, then B t > (6) B or P (St > x) > (6) P (S > x). Then we construct new system Σ̃ with the same input, but each time when St > a = const we resample S with d.f. B. Using classical regeneration approach for new system Σ̃ we obtain lower (upper) bound for QoS parameter of the original system Σ (for example, queue size). If L has NBU (NWU) distribution, then At 6 (>) A or P (Lt > x) 6 (>) P (L > x). Then we construct new system Σ̃ with the same service times, but each time when Lt 6 a = const we resample L with d.f. A. In this case we can use classical regeneration approach for new system Σ̃ to obtain lower (upper) bound for QoS parameter of the original system Σ. In the case of NWU service times or NBU interarrival times we obtain more frequent classical regenerations that allows to accelerate simulation. 5 Conclusion We present three techniques which allow to construct and speed-up classic regen- erations in multiserver systems: artificial regeneration based on splitting prop- erty, regenerative envelopes and accelerated simulation for systems with NWU service times of NBU input Our future work will focus on the comparison of the simulation results efficiency. ACKNOWLEDGEMENTS The publication was carried out under state order to the Karelian Research Centre of the Russian Academy of Sciences (Institute of Applied Mathematical Research KRC RAS). This research is partially supported by RFBR project 18- 07-00147 and by the Ministry of Education and Science of Russian Federation (project no. 14.580.21.0009, unique identifier RFMEFI58017X0009). References 1. Andradóttir, S., Calvin, J.M., Glynn, P.W.: Accelerated Regeneration for Markov Chain Simulations. Probability in the Engineering and Informational Sciences 9(04), 497–523 (1995). https://doi.org/10.1017/S0269964800004022 2. Andronov, A.: Artificial regeneration points for stochastic simulation of complex systems. In: Simulation Technology: Science and Art. 10th European Simulation Symposium ESS’98. pp. 34–40. SCS, Delft, The Netherlands (1998) 3. Asmussen, S.: Applied probability and queues. Springer, New York (2003) 4. Barlow, R.E., Proschan, F.: Mathematical theory of reliability. With contributions by L. C. Hunter. The SIAM Series in Applied Mathematics, John Wiley and Sons (1965) 5. Crane, M.A., Lemoine, A.J.: An Introduction to the Regenerative Method for Simulation Analys, Lecture Notes in Control and Information Sciences, vol. 4. Springer, Berlin, Heidelberg (1977). https://doi.org/10.1007/BFb0007339 6. Foss, S.G., Kalashnikov, V.V.: Regeneration and renovation in queues. Queue- ing Systems 8(1), 211–223 (1991), http://link.springer.com/article/10.1007/ BF02412251 7. Glynn, P.: Wide-sense regeneration for Harris recurrent Markov pro- cesses: an open problem. Queueing Systems 68(3-4), 305–311 (2011). https://doi.org/10.1007/s11134-011-9238-x 8. Glynn, P.W., Iglehart, D.L.: Simulation methods for queues: An overview. Queue- ing Systems 3(3), 221–255 (1988), http://link.springer.com/article/10.1007/ BF01161216 9. Morozov, E.: Coupling and stochastic monotonicity of queueing processes. PetrSU, Petrozavodsk (2013) 10. Morozov, E., Rumyantsev, A., Peshkova, I.: Monotonicity and stochastic bounds for simultaneous service multiserver systems. In: Ultra Modern Telecommuni- cations and Control Systems and Workshops (ICUMT), 2016 8th International Congress on. pp. 294–297. IEEE (2016), http://ieeexplore.ieee.org/abstract/ document/7765374/ 11. Sigman, K., Wolff, R.W.: A review of regenerative processes. SIAM Review 35(2), 269–288 (1993)