Introduction

Modern computer systems, as a rule, have a multiserver or multicore architecture, 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 transfer to the users in modern networks also requires coordinated performance of several devices, as well as work with several servers, for solving problems of reducing 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 methods of steady-state simulation and output analysis of complex queueing systems [3,5,8,7,6,11].

Standard classical regeneration method uses visits of the process (describing 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 classical 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 exponential 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 systems 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 remaining 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.

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 f 0 , such that:

f pf 0 .(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 (stochastically equivalent to the r.v. T ) as follows:

T = IT 0 + (1 − I)T 1 ,(2)

where T 0 has density f 0 , and T 1 has density

f 1 = f − pf 0 1 − p .(3)

Now we can replace the original r.v. T with a triple [T 0 , T 1 , 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 of exponentially splitted r.v. T is a triple (T 0 , T 1 , I) where T 0 is left truncated exponentially distributed r.v. with truncation point τ 0 , rate λ, density

f 0 (x) = 0, x τ 0 , λe −λ(x−τ0) , x > τ 0 ,

and T 1 has density

f 1 (x) = f (x) 1−p , x τ 0 , 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 = (X 1 , . . . , X n ) and continuous component T = (T 1 , . . . , T m ) 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 t i such that T j * (t i ) = 0 for some j * ∈ {1, . . . , m} (in this case we say that the ith event occurring at t i 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(t i +) = G j * (X(t i −)) ,(5)

where G j * is the recursion related to type j * event, or, in general, assume the change is governed by some Markovian probability

P j * (x, x ) = P j * {X(t i +) = x |X(t i −) = x}.

The component X(t) does not change if T(t) > 0, that formally means

X(t) = X(t i +), i = max{k : t k < t}.

We also assume that in the vector T the zeroed at t i − component T j * is initialized at t i + from some distribution with density f j * , that is,

f j * (u, x, x ) = P(T j * (t i +) ∈ du|X(t i −) = x, X(t i +) = x ). (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 T j there exist constants τ 0 (j), λ(j), p(j) such that (4) holds for the density f j (u, x, x * ), with corresponding density f 0,j (u) = λ j e −λj (u−τ0(j)) and f 1,j defined appropriately. At that, we introduce a split process Θ = {X(t), T (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

Z j * (t i +) = 1, if X(t i +) = x * and I j * (t i +) = 1, 2, otherwise,(7)Z k (t i +) = Z k (t i −), k = j * ,(8)

where t i is the epoch of type j * event, and I j * is the splitting indicator corresponding to the component T j * with success probability p(j * ). However, starting from Z j * (t i +) = 1 at time t i , the component Z j * makes a step between the events:

Z j * (t) = 1, t i < t < t i + τ 0 (j * ) 0, t t i + τ 0 (j * ),(9)

that is, the phase Z j * is changed from one to zero when not less than τ 0 (j * ) time passes since the epoch t i corresponding to the last event of j * th type (we call this time epoch a pseudo-event). In case Z j = 1, the component Z j 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 estimation 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 X i (t) is the mode of ith server, i = 1, . . . , m, X 1 (t) = 1 if the server is active, and X 1 (t) = 0 if the server is in the sleep mode. Continuous components: T 0 (t) is the time before the next arrival epoch, and T i (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 T 0 , . . . , T m are in exponential phases, that is, Z(t+) = 0.

The second approach of defining the regeneration times is to construct stochastically equivalent system with the same input, but sampling from either f 0 or f 1 at every random variable initialization time epoch. Then the original r.v. T is replaced by a stochastically equivalent splitting representation T = IT 0 +(1−I)T 1 . 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 I i , i 1 as a Bernoulli 0-1 r.v. with success probability p. Sample service times Tj either from density f 0 , or from f 1 , using relation (2). Fix integer ν 0 (fixed number of customers in system) and define

β n+1 = inf k > β n : ν k = ν 0 m, I i = 1, i ∈ M k , M βn ∩ M k = ∅ , n 0, (11)

where

M n = {i : t i t n < t d i }

is the set of the customers presented in new system at arrival instant t n and t d i 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 distributed service times. It is easy to see that Θ n := {ν n , I j , 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 }.

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 t n , the iid interarrival times L n = t n+1 − t n and the iid service times S n , n 1. Denote ν n the number of customers at instant t − n , Q n the number of customers waiting in the queue at instant t − n , so

Q n = 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 corresponding variables in the system Σ we endow with tildes. Assume that the service time distributions in both systems are ordered as F S (x) F S (x), x 0, that is S st S (stochastically). Moreover, by construction, L = st L. Then the coupling method allows to take tn = t n and S n Sn w.p.1. Define the set M n = {i : t i t n < z i } of the customers, which are being served in the system Σ at instant t n , 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 a b < ∞ and define

β n+1 = inf k > β n : ν k = ν 0 , S i (k) ∈ (a, b), M β n M k = ∅, i ∈ M k .

(12) Then, at the (discrete) instant k = β n+1 , we replace all the non-zero remaining times S i (k), i ∈ M k , by the upper bound b.

It is easy to see that Θ n := {ν n , S i (n), i ∈ M n }, n 1, is a Markov process 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 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 ∈ M k , by the lower bound a.

The systems Σ, Σ, Σ have zero initial state, and moreover (in an evident notation) L n = L n = L n , S n S n S n , n 1. Then it follows from coupling property that

ν n ν n ν n , Q n Q n Q n , 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 EQ 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 regeneration simulation. We say that the distribution F of r.v. T has increasing (decreasing) failure rate IFR (DFR) if

F t (x) = P (T t + x|T > t) = P (T x) = F (t + x) − F (t) F (t) (14) is increasing (decreasing) in t. Remark, that lim x→0 F t (x) x = f (t) F (t) = r(t).(15)

Function r(t) is called failure rate and it is easy to see that

F (t) = e − t o r(x)dx

. 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

F (t) ( ) e −t µ 1 , t < µ 1 .(16)

We say that the distribution F of r.v. T has increasing (decreasing) failure rate in average IFRA (DFRA) if

1 t t 0 r(u)du increases in t. (17) Distribution F of r.v. T is called New Better (Worse) than Used NBU (NWU) if F t F or F (t + x) 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) interarrival time (with distribution A) in m-server system Σ. Define r.v. S tremaining service time at t − and B t (x) = P (S t + x|S t) = P (S t x), r.v. L t -remaining interarival time at t − and A t = P (L t + x|L t) = P (L t 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 ( ) B or P (S t x) ( ) P (S x).

Then we construct new system Σ with the same input, but each time when S t 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 A t ( ) A or P (L t x) ( ) P (L x).

Then we construct new system Σ with the same service times, but each time when L t 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.

Conclusion

We present three techniques which allow to construct and speed-up classic regenerations in multiserver systems: artificial regeneration based on splitting property, 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.