Next Article in Journal
A Transfer Learning Approach: Early Prediction of Alzheimer’s Disease on US Healthy Aging Dataset
Previous Article in Journal
Effort-Aware Fault-Proneness Prediction Using Non-API-Based Package-Modularization Metrics
Previous Article in Special Issue
Inter-Departure Time Correlations in PH/G/1 Queues
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Regenerative Analysis and Approximation of Queueing Systems with Superposed Input Processes

by
Irina Peshkova
1,†,
Evsey Morozov
1,2,3,† and
Michele Pagano
4,*,†
1
Department of Applied Mathematics and Cybernetics, Petrozavodsk State University, Lenin str. 33, 185910 Petrozavodsk, Russia
2
Institute of Applied Mathematical Research, Karelian Research Centre of Russian Academy of Sciences, Pushkinskaja str. 11, 185910 Petrozavodsk, Russia
3
Moscow Center for Fundamental and Applied Mathematics, Lomonosov Moscow State University, GSP-1, Leninskie Gory, 119991 Moscow, Russia
4
Department of Information Engineering, University of Pisa, Via G. Caruso 16, 56122 Pisa, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2024, 12(14), 2202; https://doi.org/10.3390/math12142202
Submission received: 23 May 2024 / Revised: 23 June 2024 / Accepted: 25 June 2024 / Published: 13 July 2024
(This article belongs to the Special Issue Advances in Queueing Theory, 2nd Edition)

Abstract

:
A single-server queueing system with n classes of customers, stationary superposed input processes, and general class-dependent service times is considered. An exponential splitting is proposed to construct classical regeneration in this (originally non-regenerative) system, provided that the component processes have heavy-tailed interarrival times. In particular, we focus on input processes with Pareto interarrival times. Moreover, an approximating G I / G / 1 -type system is considered, in which the independent identically distributed interarrival times follow the stationary Palm distribution corresponding to the stationary superposed input process. Finally, Monte Carlo and regenerative simulation techniques are applied to estimate and compare the stationary waiting time of a customer in the original and in the approximating systems, as well as to derive additional information on the regeneration cycles’ structure.

1. Introduction

In this research, we consider a single-server FIFO (First In First Out) queueing system i = 1 n G I i / G / 1 , with a superposition of n independent renewal stationary inputs (components) corresponding to customers of different classes. The superposition may include component processes with general (component-dependent) interarrival times. The service time distributions are also assumed to be general and class-dependent. (In what follows, we use the terms ‘distribution’ and ‘distribution function’ as synonymous terms.)
The superposed input process system has been intensively studied since the 1970s (see, for example, [1,2,3,4,5,6,7], among others). It is well known that, unless the components are Poissonian, the superposed input process system is not (classically) regenerative. On the other hand, the regeneration property is quite important for both the stability and performance analysis of the queueing system. Indeed, it admits accurate estimation of stationary performance metrics on the basis of regenerative simulation [8,9,10,11], which allows one to use the independence of regeneration cycles (of a basic process) in the framework of the regenerative variant of the Central Limit Theorem (CLT) in order to estimate the required metrics with a given precision. (In Section 5, we discuss this topic in more detail.) What is especially important is that the values of a queueing process, within a regeneration cycle, are typically highly correlated, and a direct application of the classical variant of the CLT is not justified.
In general, the superposed process generated by independent renewal processes can be transformed, under mild assumptions, into a regenerative process with the so-called one-dependent regeneration cycles (see [6,12]). This construction, which is based on splitting and coupling procedures, is rather complicated and hardly applied in practice to evaluate the required performance measures by simulation. (We will discuss it in Section 3.) Moreover, to apply this approach for accurate estimation, we must overcome an even bigger problem; we need to construct regeneration events of the queueing process (not only of the input process), which, in turn, requires the synchronization of a regeneration point of the input process with an empty state of the system. For this reason, we instead will base our analysis on the construction of the artificial regeneration obtained by exponential splitting. (For a deeper understanding of the splitting method, we refer readers to the fundamental monograph [13].) A key property opening the possibility of this construction in the initially non-regenerative process is that the interarrival distributions of the component processes are assumed to be heavy-tailed. Thus, in this research, we focused on the case where the superposition consists of n 1 Poisson inputs and n 2 inputs with Pareto interarrival times. The splitting approach, in keeping the distributional properties of the original process, allows one to use regenerative simulation for accurate estimations of the required stationary performance metrics. To the best of our knowledge, this approach is realized for the first time in systems with a superposed input process, and this finding is the main contribution of this research.
Our second main purpose is to study the basic system by means of another classically regenerative system with a single renewal input process. It is not a new approach as, for instance, it has been realized in a series of papers, such as [3,5,7]. In these works, the so-called stationary-interval method (based on the Palm approach) was used to approximate the distribution of the interarrival times in the superposition process by a stationary distribution with independent identically distributed (iid) intervals. We emphasize that such a replacement is indeed an approximation because the successive intervals in the superposed process are dependent in general. The latter issue has been analyzed in [1], in which the correlation structure of the successive interevent times was considered in detail. Another motivation of the approximation, as we mentioned, is that the target processes in the basic system (for example, workload or queue size) are not regenerative; however, a regenerative structure is highly desirable both to establish the stability of the system and to accurately estimate its stationary performance metrics. An issue in the study of queueing systems with superposed input (with independent renewal components) concerns the approximation of such a process by a Poisson process. It is a well-known problem going back to A. Khintchine [14]. More exactly, when the number of summands (components) increases, while the (properly defined) traffic intensity ρ remains bounded (and not close to 1), then the superposed process approaches a Poisson process, with the rate being the sum of the rates of the summands. On the other hand, it has been shown in [2] that, as the number n of the summands increases, this approximation works well only if the product n ( 1 ρ ) 2 1 . The latter result indicates that the Poissonian approximation is not suitable when the number of components is small and/or the traffic intensity ρ is close to 1, that is, the system is close to ‘saturation’. The saturated system (in the so-called heavy-traffic regime) with the superposed input process is studied in the paper [15]. Notice also that the basic results for ordinary renewal process, obtained originally by W. Smith [16] and D. Cox [17], have been extended to the superposition of independent renewal processes in [4].
Thus, our second important contribution is the detailed analysis of the approximating G I / G / 1 -type system when the interarrival time distribution is a mixture containing Pareto distributions, and when the service time distribution is also a mixture of distributions that describe the service times of the customers in the component processes. Indeed, we show that, in the approximating system, the stationary workload distribution is analytically available when the service times of the different classes of customers have a common exponential distribution (Section 4.2); this is also a new contribution of this research. Moreover, we compare the accuracy of the approximation in terms of the Kolmogorov distance between the (empirical) distributions of the stationary waiting times in the original and approximated systems (Section 6.1). In particular, in spite of the presence of heavy-tailed components in interarrival distribution, the simulation results show quite acceptable distance values (less than 5 % ) for most of the considered cases (when the traffic intensity ρ < 0.6 ), and it is consistent with the results of the abovementioned works [3,5,7]. Hence, the proposed approximation may be acceptable if the splitting construction turns out to be too complicated, e.g., if the number of component processes is very large. Another possible application of the approximating system is as a base for appropriate lower/upper bounds for the stationary workload process, exploiting the known monotonicity properties of the classical queueing system under some mild assumptions. (In this regard, see [18] and recent papers of one of the authors [19,20] on the monotonicity of systems with mixed service times.)
As mentioned above, the approximation of the superposed input process by a renewal process is indeed a mixed distribution, consisting of an appropriately weighted sum of the stationary distributions of the remaining renewal times (integrated-tail distributions) in the component stationary renewal processes. This distribution, based on the Palm probability concept, is well known (see, for instance, [7,17]) and is the exact one-dimensional distribution of the interarrival times, which, however, does not capture possible dependencies between the neighboring intervals. For this reason, reducing the superposed input process to an ordinary renewal process is indeed an approximation. We note that, in general, a dependence between the adjacent intervals cannot be ignored in the careful study of the superposed input, as shown in [5].
The rest of the paper is organized as follows. In Section 2, we give some important preliminary results from the theory of renewal and point processes with elements of the Palm theory. In Section 3, we outline the construction of one-dependent regeneration using the approach from [12] and highlight the difficulties arising in its practical realization in simulation. Nevertheless, we use some basic ideas of this approach to construct regeneration events by the exponential splitting of heavy-tailed interarrival times.
In Section 4, we examine an approximating G I / G / 1 system with iid interarrival times, generated by the Palm distribution of the stationary intervals in the original superposed input process. We then focus on special cases where we find explicit expressions for the Laplace–Stieltjes transform (LST) of the stationary workload. In Section 5, we introduce exponential splitting and explain how it is used to obtain classical regenerations in the superposed process generated by Pareto components. Next, in Section 6, we analyze our theoretical findings through discrete-event simulations. Specifically, we first evaluate the accuracy of the input process approximation in terms of Kolmogorov distance for different system parameter settings. Then, we apply exponential splitting to construct a confidence estimate of the empirical mean of the stationary waiting time in the original system. Finally, Section 7 concludes the paper by providing an overview of its main contributions.

2. Preliminary Results from the Theory of Renewal and Point Processes

In this section, we present theoretical results from the theory of renewal and point processes [21,22]. These results are then utilized in the analysis of queueing systems with the superposed input process generated by independent renewal processes. Among the most important sources on renewal theory, besides the pioneering work [23], we also mention the recognized monographs [24,25].
To start, let us construct a stationary renewal process. To do this, consider iid non-negative random variables (r.v.s) { X k , k 2 } and an independent r.v. X 1 0 . We denote by F the distribution function of an r.v. X, which is a stochastic copy of any X k , k 2 . It is assumed that the renewal rate λ = 1 / E X < . Let
G ( x ) = λ 0 x ( 1 F ( u ) ) d u , x 0 ,
be the distribution of the first r.v. X 1 . Also, denote by T k = 1 i k X i the renewal instants, where k 1 . Hence, N s ( t ) : = max { k > 0 : T k t } represents the number of renewal events that occur in the time interval [ 0 , t ] , with N s ( 0 ) = 0 . The process { N s ( t ) , t 0 } is a time-stationary renewal process [22,26]. Specifically, the stationary renewal process satisfies E N s ( t ) = λ t , t 0 , and the distribution of the increments N s ( t + u ) N s ( u ) does not depend on u 0 for each fixed t . There are other equivalent statements to identify a stationary renewal process (see, for example, Chapter 2 in [25]). It is also known that (1) represents the distribution of the stationary remaining (and attained) renewal time, and it is referred to as the integrated-tail distribution [12].
In the following analysis, we will consider n independent stationary renewal processes that model the arrivals of different classes of customers. Then, the renewal points { T k ( i ) , k 1 } correspond to the arrival instants of class-i customers from process i, and let X ( i ) denote the generic interarrival time (with a rate of λ i = 1 / E X ( i ) ) in process i, where i = 1 , , n . These processes will compose a merged superposed input process. To provide a more detailed description, let T = { T k ( i ) , k 1 ; i = 1 , , n } be the overall set of renewal points for all n processes, and let { T k } be the ordered sequence of these renewal instants. In other words, T 1 is the first smallest element of T , T 2 is the second smallest element, and so on.
Assuming, for simplicity, that there are no batch arrivals, it follows that the sequence { T k } satisfies
0 < T 1 < T 2 <
It is clear that, unlike a renewal process, the distances between events (arrivals) in the superposed process are, in general, dependent.
Since the component renewal processes are assumed to be time-stationary, then the superposed process, defined (in evident notation) as
N s ( t ) : = N s ( 1 ) ( t ) + + N s ( n ) ( t ) , t 0 ,
is also time-stationary (counting), but not renewal in general. (A detailed proof can be found, for example, in [21], Chapter 1.)
To explain this process in more detail, we recall some findings from the theory of stationary point processes, while maintaining the previously established notation. Following [21,25], we consider a set of points (occurrence times) { T k } on the positive real line that satisfy (2), with the condition that T 0 = 0 is an occurrence time. This process is assumed to be interval-stationary and is referred to as the stationary Palm process. We also consider a counting process { N ( t ) , t 0 } , which records the number of points in the interval [ 0 , t ] , formally defined as
N ( t ) = max { k 0 : T k t } , t 0 ,
with the assumption that N ( 0 ) = 1 . Let X k : = T k T k 1 denote the interval between the kth and ( k 1 ) st points, where k 1 . These intervals are identically distributed with distribution function F. This new process is also referred to as the stationary Palm process [7,21,22].
In the following discussion, we will distinguish between the Palm process (4) and the related time-stationary (also counting) process { N s ( t ) , t 0 } . Note that an interval sequence { X i } (forming the Palm process) is stationary if the joint distribution of the k-tuple ( X i 1 + h , , X i k + h ) is independent of the positive integer h for all integers k > 0 and all non-negative ( i 1 , , i k ) . It is worth pointing out that the Palm process { N ( t ) } and the associated sequence { X k } are usually not both stationary at the same time [7]. However, there is a one-to-one correspondence between the stationary counting process { N ( t ) } and the stationary sequence { X k } based on the Palm theory [21]. To establish this correspondence, let us introduce the Palm function  φ k ( t ) = P ( N ( t ) = k ) , which represents the probability of having k points ( k 0 ) in a time interval of length t provided that the beginning of the interval is an occurrence point. Hence, Palm functions describe the dynamics of the Palm process (4). The distribution of the process { N s ( t ) } associated with the Palm process (4) can be expressed via Palm functions by the following relations [7,14]:
P ( N s ( t ) = k ) = λ 0 t ( φ k 1 ( u ) φ k ( u ) ) d u , k 1 ,
P ( N s ( t ) = 0 ) = 1 λ 0 t φ 0 ( u ) d u = : 1 G ( t ) ,
G ( t ) = λ 0 t ( 1 F ( x ) ) d x , t 0 ,
implying, in particular, the following key connection:
φ 0 ( t ) = 1 F ( t ) = : F ¯ ( t ) , t 0 ,
where, as already stated, F is the distribution function of the intervals in the stationary Palm process. (In what follows, for any distribution function F, we denote by F ¯ its tail.) Note that, by construction, the stationary process { N s ( t ) } satisfies both required conditions (see [7]):
E N s ( t ) < , lim t N s ( t ) / t = : λ < .
Recall that { N s ( i ) ( t ) } represents the ith component (stationary renewal process) and we will supply with index i the corresponding quantities in process i. Note that, by (1), the first interval in the process i has distribution
G i ( t ) = λ i 0 t ( 1 F i ( u ) ) d u ,
and rate λ i = 1 / E X ( i ) , i = 1 , , n . Let λ = λ 1 + + λ n . It is known [7,14,21] that the Palm functions φ k ( t ) of the superposed stationary Palm process { N ( t ) } are expressed via the original time-stationary processes { N s ( i ) ( t ) } and their stationary Palm versions { N ( i ) ( t ) } as follows:
φ k ( t ) = i = 1 n λ i λ P j = 1 i 1 N s ( j ) ( t ) + N ( i ) ( t ) + j = i + 1 n N s ( j ) ( t ) = k = i = 1 n λ i λ P ( N ( i ) ( t ) = k i ) j , j i n P ( N s ( j ) ( t ) = k j ) ,
where the summation is taken over all different partitions of k 0 such that k = k 1 + + k n . In particular, by (8), the tail of the distribution F satisfies [7]:
F ¯ ( t ) = φ 0 ( t ) = P ( N ( t ) = 0 ) = i = 1 n λ i λ ( 1 F i ( t ) ) j , j i G j ( t ) = i = 1 n λ i λ ( 1 F i ( t ) ) j , j i λ j t ( 1 F j ( u ) ) d u , t 0 .
It is easy to verify that F ¯ ( 0 ) = φ ( 0 ) = 1 , as required. We emphasize that the Palm process { N ( t ) } is interval-stationary, while the process { N s ( t ) } is time-stationary. Furthermore, since, for a single component process with interevent distribution F and rate λ , we have
λ 0 F ¯ ( t ) d t = 1 ,
Equation (5) implies that P ( N s ( t ) = 0 ) = G ( t ) , satisfying (6). Additionally, expression (10) has a straightforward intuitive interpretation.
Since a renewal process is a particular case of the point process, it is worth mentioning that F represents the exact marginal (one-dimensional) distribution of the interarrival time in the superposed Palm process { N ( t ) } . However, this distribution does not capture the dependence between adjacent intervals. As a result, expression (10) only provides an approximation of the original superimposed process, embedded at the instances of the points’ occurrences.
Therefore, it is crucial to examine the accuracy of this approximation using various metrics. Specifically, as detailed in Section 6, the approximation (10) of the original i = 1 n G i / G i / 1 system yields an acceptable Kolmogorov distance between the empirical distributions of the stationary waiting times when ρ < 0.6 .
Remark 1. 
Representations (10), (9) hold for any superposition in which the components { N s ( i ) ( t ) , t 0 } are time-stationary processes, not necessary renewal. In this general case, the intensities are defined as λ i : = E N s ( i ) ( 1 ) , i = 1 , , n [7].

3. Splitting of a Superposed Process

In general, for superposed processes, it is only possible to construct the so-called one-dependent regenerations. In this section, we will briefly describe this construction following [12]. Roughly speaking, in a one-dependent regenerative process, two consecutive regeneration cycles are dependent, but the cycles separated by at least one cycle are independent. However, as we will discuss below, this approach is challenging to apply in simulations. Therefore, in Section 5, we will introduce exponential splitting, which can be effectively applied when the renewal intervals follow heavy-tailed distributions [22].
Let us describe the construction of one-dependent regenerations for the superposition of two renewal processes, with inter-renewal time distributions F i , i = 1 , 2 . Firstly, we assume that both these processes are positive recurrent, i.e., the mean inter-renewal time is finite for both processes. Furthermore, let us assume that F 1 satisfies the following condition:
F 1 ( d x ) c d x for x ( a , a + 2 b )
for some positive constants a , b , and c. This means that F 1 is a particular case of a spread-out distribution F. (In the general case, some convolution of F with itself must satisfy (11); see [22].)
Denote by R i ( t ) the (right-continuous) remaining renewal time at instant t in the ith renewal process, i = 1 , 2 ; t 0 . Then, under assumption (11), the distribution of R 1 ( t ) has the following smoothness property:
P ( R 1 ( t ) B ) δ ν ( B ) , t C ,
for each Borel set B on ( 0 , b ) , where ν is the uniform measure on ( 0 , b ) , δ > 0 , and C is some constant. (Condition (12) is sometimes called minorization.) That is, for large enough t, the distribution of R 1 ( t ) has a uniform component. Let R t = ( R 1 ( t ) , R 2 ( t ) ) and denote by ( t n ( i ) , n 1 ) the renewal points of the process i = 1 , 2 . Let
t 0 = 0 , L 0 = R 1 ( t 0 ) , t 1 = min ( t n ( 2 ) : t n ( 2 ) t 0 + L 0 + C ) ,
where C satisfies (12). Also, denote by { U i , i 1 } the iid r.v.s with uniform distribution ν (in what follows, we will denote it as U ν .) Also, let { V i , i 1 } be iid 0 1 Bernoulli r.v.s such that P ( V = 1 ) = δ . Take
R 1 ( t 1 ) = U 1 V 1 + ( 1 V 1 ) Z 1 ,
where the r.v. Z 1 has the rest distribution (see (12))
P ( R 1 ( t ) · ) δ ν ( · ) 1 δ .
Continuing as in (13), define recursively
t k = min ( t i ( 2 ) : t i ( 2 ) t k 1 + L k 1 + C ) , L k = R 1 ( t k ) = U k V k + ( 1 V k ) R k , k 1 ,
where U k , V k , R k are independent of all U l , V l , R l , l < k . Eventually, let
T ^ 1 = min ( t i : V i = 1 ) , T ^ k + 1 = min ( t i > T ^ k : V i = 1 ) , k 1 .
Then, the sequence { T ^ k } { t k ( 2 ) } and moreover, at each instant T ^ k , the remaining times vector has distribution
R T ^ k = ( R 1 ( T ^ k ) , R 2 ( T ^ k ) ) ν F 2 , k 1 ,
which is independent of the pre-history prior to instant T ^ k 1 . (Here, ν F 2 denotes the distribution generated by the marginal distributions ν and F 2 .) Thus, { T ^ k } are one-dependent regeneration points for the superposed remaining time process { R t } , and, according to (14)–(16), they constitute a subsequence of the renewal points of input 2. The extension of this construction to arbitrary n renewal processes is straightforward [12].
Remark 2. 
It is worth pointing out that, under classical regeneration, the distribution of R T ^ k in (17) would be independent of the pre-history prior to instant T ^ k .
The primary challenge in practically implementing regenerations in the superposed process is that the constants in Equation (12) are not explicitly defined. An even more difficult problem is the need to synchronize these regeneration points with empty states of the queueing system to obtain regenerations of the entire system. However, in some cases, particularly when the inter-renewal intervals follow heavy-tailed distributions, this difficulty can be overcome. In this case, as we will show in Section 5, by utilizing the aforementioned construction of one-dependent regenerations, the splitting procedure can be effectively applied to construct the classical regenerations of the superposed process.

4. GI / G / 1 Approximation

In this section, we assume that in the basic system i = 1 n G i / G i / 1 (further denoted by Σ ) with n independent renewal processes, the ith component has iid interarrival times { τ k ( i ) , k 2 } with distribution A i , and the distribution of the first interarrival time (up to the first arrival) satisfies (cf. (1), (8))
P ( τ 1 ( i ) x ) = 1 E τ 2 ( i ) 0 x ( 1 A i ( u ) ) d u , i = 1 , , n .
Customers from input i (called class-i customers) have iid (class-dependent) service times { S k ( i ) , k 1 } with a general distribution B i , i = 1 , , n . In the following analysis, we will focus on exponential and Pareto distributions for A i . Then, for specific cases, we will construct the Palm distribution (10) to use it as the distribution of the iid interarrival times in an approximating G I / G / 1 queueing system (denoted by Σ ^ ).

4.1. Exponential–Pareto Input Process

Assume that out of the total n inputs, the first n 1 components have exponential interarrival times, τ k ( i ) E x p ( λ i ) , k 1 ; i = 1 , , n 1 ,
A i ( x ) = 1 e λ i x , x 0 , λ i > 0 , i = 1 , , n 1 ,
while the remaining inputs have interarrival times following a type II Pareto distribution, P a r e t o ( x 0 , α j ) , j = 1 , , n 2 , i.e.,
P ( τ k ( j ) x ) = : A j ( x ) = 1 x 0 x 0 + x α j , x 0 , k 1 , x 0 > 0 , α j > 1 , j = 1 , , n 2 .
Let μ i : = 1 / E S k ( i ) denote the mean service rate of the i-th class customers, where i = 1 , , n . In the following, we use index i to represent the quantities associated with the Poisson inputs and index j to represent the quantities related to Pareto input processes. By substituting the distributions (18) and (19) into Formula (10), we can obtain the distribution of the stationary renewal input process in the following form:
A ( x ) = 1 p e λ x x 0 x 0 + x α n 2 ( 1 p ) e λ x x 0 x 0 + x α n 2 + 1 ,
where
λ = λ 1 + + λ n 1 , α = α 1 + + α n 2 ,
and
p = λ x 0 λ x 0 + α n 2 .
Note that the tail distribution of (20) is the following mixture:
A ¯ ( x ) = p F ¯ 1 ( x ) + ( 1 p ) F ¯ 2 ( x ) ,
of the two-tail distributions
F ¯ 1 ( x ) = e λ x x 0 x 0 + x α n 2 , F ¯ 2 ( x ) = e λ x x 0 x 0 + x α n 2 + 1 ,
with the mixing proportion p defined by (21). Denote by τ the (generic) interarrival time following distribution (20) in the renewal process defined by distribution (22). Let distribution F 1 define an r.v. Y, and distribution F 2 define an r.v. Z. Then, the r.v. τ can be expressed as a two-component mixture:
τ = I Y + ( 1 I ) Z ,
where I is the indicator function with P ( I = 1 ) = p and
Y = min ( Y 1 , Y 2 ) with Y 1 E x p ( λ ) , Y 2 P a r e t o ( x 0 , α n 2 ) , Z = min ( Y 1 , Z 1 ) with Z 1 P a r e t o ( x 0 , α n 2 + 1 ) .
Moreover, the tail distribution of the service time is also determined by the following mixture of the predefined class-dependent distributions { B k } , namely:
B ¯ ( x ) = i = 1 n 1 p i B ¯ i ( x ) + j = 1 n 2 p j B ¯ j ( x ) ,
where
p i = λ i x 0 λ x 0 + α n 2 , i = 1 , , n 1 ; p j = α j 1 λ x 0 + α n 2 , j = 1 , , n 2 .
Thus, instead of the original system Σ , we consider a classical system Σ ^ , in which the renewal input process has interarrival time distribution (20) and the service times have a tail distribution defined by the mixture (25). Now, we give the stability condition of the new system Σ ^ . This system is regenerative, and regenerations (of all processes describing its dynamics) occur when the arrivals meet the system idle. Denote by T the generic regeneration period, i.e., the distance between two neighbor regeneration points. It is well known that the stability of a regenerative system is equivalent to positive recurrence, in which case E T < [22]. Denote by S the service time with tail distribution (25); recalling (24), the stability criterion of the system Σ ^ is (see Section 2.2, Chapter 2 in [26]):
ρ : = E S E τ < 1 ,
where
E S = i = 1 n 1 p i E S ( i ) + j = 1 n 2 p j E S ( j ) = i = 1 n 1 p i μ i + j = 1 n 2 p j μ j .
The first step in calculating E τ is to determine E Y :
E Y = E min ( Y 1 , Y 2 ) = 0 x f Y 1 ( x ) x f Y 2 ( y ) d y d x + 0 y f Y 2 ( y ) y f Y 1 ( y ) d x d y = 0 x λ e λ x x 0 α n 2 ( x 0 + x ) α + n 2 d x + 0 x e λ x ( α n 2 ) x 0 α n 2 ( x 0 + x ) α + n 2 1 d x = 1 λ 1 ( α n 2 ) e λ x 0 ( λ x 0 ) α n 2 Γ ( α + n 2 , λ x 0 ) .
Similarly, we find
E Z = x 0 e λ x 0 ( λ x 0 ) α n 2 Γ ( α + n 2 , λ x 0 ) ,
resulting in
E τ = p E Y + ( 1 p ) E Z = p λ = x 0 λ x 0 + α n 2 .
Let
ρ i : = E S ( i ) E τ ( i ) = λ i μ i , i = 1 , , n 1 , ρ j : = E S ( j ) E τ ( j ) = α j 1 x 0 μ j , j = 1 , , n 2 ;
then, it is easy to check, using (27) and (28), that
ρ = i = 1 n 1 ρ i + j = 1 n 2 ρ j .
Remark 3. 
Condition (26) is also the stability criterion of the original system with the stationary input process, not necessary renewal. In this case, the proof is based on the so-called Loynes’ construction; see, e.g., Chapter 2 in [21].

4.2. Case 1: Common Exponential Service Time

An important special case of the system is the model i = 1 n G i / M / 1 in which all arrivals have the same (iid) exponential service times, i.e., μ 1 = = μ n = μ . Such a setting has been motivated in the paper [5] where it is compared with an M / M / 1 queue for n > 10 . In this case, the approximating system is a stable classical G I / M / 1 queue ( ρ < 1 ), for which the analytical form of the stationary workload distribution is known [27]. Additionally, the LST obtained for the input stream of the form (20) allows us to compare simulation results for the waiting time distribution of the queueing system fed by the actual superposed input and its renewal approximation with the corresponding analytical solution (see Section 6.1).
Denote by ψ X ( z ) the LST of an r.v. X with distribution F:
ψ F ( z ) = E e z X = 0 e z x d F ( x ) , z 0 .
Then, the distribution of the stationary waiting time W can be expressed as
F W ( x ) = 1 σ e μ ( 1 σ ) x , x > 0 ,
where σ is the root inside the unit circle of the equation
z ψ A ( μ ( 1 z ) ) = 0 .
After some simple algebra, we can obtain the LST ψ A ( z ) of the distribution (20) as follows:
ψ A ( z ) = 0 e z x d A ¯ ( x ) = 1 z 0 e z x A ¯ ( x ) d x = 1 p z 0 e z x e λ x x 0 x 0 + x α n 2 d x ( 1 p ) z 0 e z x e λ x x 0 x 0 + x α n 2 + 1 d x = 1 p z e x 0 ( z + λ ) x 0 α n 2 ( z + λ ) α n 2 1 Γ ( α + n 2 + 1 , x 0 ( z + λ ) ) ( 1 p ) z e x 0 ( z + λ ) x 0 α n 2 + 1 ( z + λ ) α n 2 Γ ( α + n 2 , x 0 ( z + λ ) ) ,
where
Γ ( ξ , x ) = x e t t ξ 1 d t ,
is the upper incomplete Gamma function. Finally, applying the known property of the Gamma function
Γ ( ξ + 1 , x ) = ξ Γ ( ξ , x ) + x ξ e x ,
we obtain
ψ A ( z ) = 1 p z z + λ + z z + λ e x 0 ( z + λ ) x 0 α n 2 ( z + λ ) α n 2 × Γ ( α + n 2 , x 0 ( z + λ ) ) p ( α n 2 ) ( 1 p ) x 0 ( z + λ ) .
The explicit expression (31) of the LST allows us to find the moments of the r.v. τ . For example, the mean is given by
E τ = ψ A ( 0 ) = p / λ = x 0 λ x 0 + α n 2 ,
and this matches the expression in (28). Then, by utilizing the LST (31), we can numerically solve (30) and use the solution σ to obtain the target distribution function in (29).

4.3. Case 2: Class-Independent Inputs

Assume that
λ 1 = = λ n 1 = : λ and α 1 = = α n 2 = : α .
Then, the tail distribution of the interarrival time becomes (cf. (23)):
A ¯ ( x ) = p e λ n 1 x x 0 x 0 + x n 2 ( α 1 ) + ( 1 p ) e λ n 1 x x 0 x 0 + x n 2 ( α 1 ) + 1 ,
where
p = n 1 x 0 λ n 1 x 0 λ + n 2 ( α 1 ) .
Since there are only two customer classes, the service time distribution becomes a two-component mixture:
B ¯ ( x ) = p B ¯ 1 ( x ) + ( 1 p ) B ¯ 2 ( x ) , x 0 .

4.4. Case 3: Two Pareto Inputs

Consider the superposition of two (stationary) renewal processes with interarrival times τ ( i ) P a r e t o ( x 0 , α i ) , α i > 1 , i = 1 , 2 . By using (10), it is easy to show that the stationary sequence of interarrival times in the new system Σ ^ follows a Pareto distribution, i.e.,
A ¯ ( x ) = x 0 x 0 + x α 1 + α 2 1 , x 0 .
In this two-class scenario, the service time distribution can be expressed as a two-component mixture with the mixing proportion
p = α 1 1 α 1 + α 2 2 ,
and the stability condition (26) becomes
ρ = α 1 1 x 0 μ 1 + α 2 1 x 0 μ 2 < 1 .

4.5. Case 4: n 1 = 1 , n 2 = 1

Finally, we consider the superposition of two inputs with exponential and Pareto interarrival times, i.e.,
A 1 ( x ) = 1 e λ x and A 2 ( x ) = 1 x 0 x 0 + x α .
The interarrival time Palm distribution (10) becomes
A ¯ ( x ) = p e λ x x 0 x 0 + x α 1 + ( 1 p ) e λ x x 0 x 0 + x α ,
where
p = λ x 0 α 1 + λ x 0 .
We define F ¯ 1 ( x ) and F ¯ 2 ( x ) as:
F ¯ 1 ( x ) = e λ x x 0 x 0 + x α 1 and F ¯ 2 ( x ) = e λ x x 0 x 0 + x α .
Then, the r.v. τ with distribution satisfying (34) can be expressed as a two-component mixture (24) of an r.v. Y with distribution F 1 and an r.v. Z with distribution F 2 (with the mixing proportion p satisfying (35)), where
Y = min ( Y 1 , Y 2 ) with Y 1 E x p ( λ ) , Y 2 P a r e t o ( x 0 , α 1 )
and
Z = min ( Y 1 , Z 1 ) with Z 1 P a r e t o ( x 0 , α ) .
In this case, the stability condition (26) becomes
ρ = λ μ 1 + α 1 x 0 μ 2 < 1 .

5. Exponential Splitting and Regenerative Estimation

In this section, we first present the construction of artificial regeneration points based on the exponential splitting of the interarrival time density. Then, we outline the regenerative approach to estimate stationary performance indicators of the queueing system with the superposed input process. (For details on regenerative estimation, see, for instance, Chapter IV in [11]).

5.1. Exponential Splitting

The splitting method, in its general form, is described in [13], while the specific variant used in this paper relies on an idea presented in [28]. The method involves replacing the original process with a stochastically equivalent process (defined on an enlarged probability space) with some desired property. In our case, it is the memoryless property that leads to a regenerative structure. This allows us to apply regenerative simulation for constructing confidence estimates and accurately analyzing the stationary performance metrics of the original process. Next, we will explain the splitting procedure, which is the core of the method.
Assume that the density g of an absolutely continuous r.v. T satisfies the inequality
g δ g 0 ,
where δ ( 0 , 1 ) is a constant and g 0 is some density. Now, let us define a new Bernoulli r.v. I T (referred to as the splitting indicator) such that P ( I T = 1 ) = δ . We can construct a new r.v. T as follows:
T = I T T 0 + ( 1 I T ) T 1 ,
where the r.v. T 0 has density g 0 , and the r.v. T 1 has density (c.f. (14))
g 1 = f δ g 0 1 δ .
As we mentioned earlier, the purpose of such a transformation is to let the r.v. T have a desired property of the distribution g 0 with a probability of at least δ . In our case, it is the memoryless property that allows us to construct regeneration points. To achieve this, we focus on exponential splitting; a positive r.v. T is exponentially split if there exist constants η > 0 and δ ( 0 , 1 ) such that
g ( x ) δ η e η x , x 0 .
That is, the splitting representation T of the r.v. T consists of an exponential r.v. T 0 with rate η (and density g 0 ( x ) = η e η x , x 0 ), and an r.v. T 1 which, according to (38), has density
g 1 ( x ) = g ( x ) δ η e η x 1 δ , x 0 .
The applicability of exponential splitting to the superposition of two Pareto inputs (see Section 6.2) is based on the following lemma.
Lemma 1. 
Assume that an r.v. T follows a P a r e t o ( x 0 , α ) distribution, with the density function given by
g ( x ) = α x 0 α ( x 0 + x ) ( α + 1 ) , α > 0 , x 0 > 0 , x 0 .
Then, under the condition
α + 1 η x 0 α / δ ,
where η and δ are defined in (39), the r.v. T has distribution (37), where T 0 is an exponential r.v. with parameter η and inequality (39) holds. Additionally, T 1 has the following distribution:
G 1 ( x ) = 1 1 1 δ x 0 x 0 + x α + δ 1 δ e η x , x 0 ,
with density g 1 satisfying (40).
Proof. 
Recall that the failure rate of an r.v. X 0 with distribution F and density f is defined as:
r X ( x ) = f ( x ) F ¯ ( x ) , x 0 .
Consider now two r.v.s X E x p ( η ) and Y P a r e t o ( x 0 , α + 1 ) . For any x 0 , the failure rate order, denoted by X r Y , implies the following inequality between the failure rates:
r X ( x ) = η α + 1 x 0 α + 1 x 0 + x = r Y ( x ) , x 0 .
This implies the stochastic order, X Y [18]. In turn, it means that the corresponding tail distributions satisfy the following inequality:
e η x x 0 x 0 + x α + 1 , x 0 .
If x 0 δ η α , then the inequality
x 0 δ η α e η x x 0 x 0 + x α + 1
holds as well. Therefore,
δ η e η x α x 0 x 0 x 0 + x α + 1 = g ( x ) ,
and inequality (39) is satisfied. □
Figure 1 illustrates the role of conditions (42) showing the graphs of the Pareto density (41) (with x 0 = 1 and α = 3 ) and exponential functions δ η e η x for different values of η and δ .
The plots highlight that the exponential function is upper-bounded by the given Pareto distribution only when both inequalities (42) are met (blue curve with δ = 0.75 , η = 4 ). Instead, the green and orange curves, violating the l.h.s. and r.h.s of (42), respectively, do not satisfy the splitting condition (39).

5.2. Regenerative Simulation and Estimation

From the simulation viewpoint, splitting means that, instead of generating an r.v. T with density g, a triple ( T 0 , T 1 , I ) is generated and then the r.v. T is constructed as in (37). The property (36), together with the specially constructed measure (38), ensures that the simulated value T is stochastically equivalent to T.
Denote by { t k } the arrival instances in the superposed input process. To construct regeneration points of the input process, denoted by { γ k } , we first select and fix an arbitrary non-Poisson component process i 0 (with the arrival instances { t k ( i 0 ) t k ( 0 ) , k 1 } ). Now consider the sequence of events:
E k ( 0 ) = { the   j th   component   process   has   exponential   phase   at   instant t k ( 0 ) , j i 0 } , k 1 .
This event means that at the kth arrival instant of input i 0 , the interarrival times in all other components have an exponential ‘phase’, obtained by splitting (i.e., the corresponding splitting indicators are equal to 1 simultaneously).
For any integer k 1 , we introduce the index n 0 ( k ) as t k = t n 0 ( k ) ( 0 ) . Now, we define the regeneration instances { γ k } of the superposed input as follows: γ 0 = 0 ,
γ k + 1 = min { i : 1 ( t n 0 ( i ) ( 0 ) > γ k ) · 1 ( E n 0 ( i ) ( 0 ) ) = 1 } , k 0 .
It is easy to verify that the difference Δ k : = γ k + 1 γ k is the number of arrivals in the superposed process between the kth and the ( k + 1 ) st occurrence of the events E ( 0 ) . These events are generated by the arrivals of class- i 0 customers, which meet the interarrival times of all other components in the exponential phase. It is easy to understand that { γ k } are indeed the regeneration points of the superposed input process, Δ k is the (discrete) length of the kth regeneration cycle of the input process, and the increments { Δ k } are iid.
If there are n 2 2 Pareto input components, then, according to the described construction, the density of each component j i 0 must satisfy the inequalities
η j x 0 α j + 1 , δ j η j x 0 α j ,
for some constants δ j ( 0 , 1 ) , η j > 0 , j = 1 , , n 2 .
For each j i 0 , the interarrival times τ k ( j ) are drawn with probability δ j from the distribution E x p ( η j ) , and with probability 1 δ j from the distribution (43) with parameters α j , x 0 , δ j , η j . Then, we construct the instances { γ k } following (45).
Now, we return to the system Σ described in Section 4. In this system, the classical regeneration points are arrival points satisfying the following two conditions: (i) they are regeneration instants of the input process (i.e., { γ k } ); (ii) at these arrival instants, the server is idle, i.e.,
β 0 = 0 , β k + 1 = min { γ i > β k : W γ i = 0 } , k 0 ,
where W i is the waiting time of the ith customer in the queue. Let Δ ^ k : = β k + 1 β k be the length of the kth regeneration cycle of type (47); then, { Δ ^ k } are iid.
To estimate the workload process by the regenerative method [11], we define the iid sequence { Y j } of the accumulated workload within regeneration cycles,
Y j = i = β j β j + 1 1 W i ,
where W i is the waiting time of the ith customer within the regeneration cycle j , j 1 . It is worth mentioning that we use a regenerative version of the CLT (called RCLT) because correlations between the values of W i belonging to the same regeneration cycle are not captured by the classical CLT. To obtain a regenerative estimate of the mean stationary waiting time, the iid observations ( Y j , Δ ^ j ) , j 1 are used to estimate the ratio [11]
r = E Y 1 E Δ ^ 1 .
Then, applying the RCLT to the iid zero mean variables { Y i r Δ ^ i , i 1 } and assuming that E ( Y 1 r Δ ^ 1 ) 2 < , we obtain the ( 1 2 ϑ ) % confidence interval (CI) for the parameter r in the following form [11]:
r ˜ N ± h ϑ S ˜ N 2 Δ ˜ N N ,
where r ˜ N : = Y ˜ N / Δ ˜ N , Δ ˜ N and Y ˜ N are the sample means of { Δ ^ i } and { Y i } , respectively, N is the number of regeneration cycles, S ˜ N 2 is the estimate of the second moment
E ( Y 1 r Δ ^ 1 ) 2 = Var ( Y 1 ) 2 r cov ( Y 1 , Δ ^ 1 ) + r 2 Var ( Δ ^ 1 ) ,
and
h ϑ = Φ 1 ( ( 1 ϑ ) / 2 ) ,
where Φ ( x ) is the Laplace function.
The finiteness of E ( Y 1 r Δ ^ 1 ) 2 , containing the ‘cycle’ variables Y 1 and Δ ^ 1 , can be easily expressed in terms of the moments of the queueing parameters (service and interarrival times). Indeed, it follows from [29] that (50) is finite if both the service time and the interarrival time have finite second moments. In the simulation settings reported in Section 6.2, this condition holds since we consider exponential service times and the superposition of two Pareto inputs, which is characterized by a Pareto distribution of the interarrival time with shape parameter α eq = α 1 + α 2 1 3 in all our experiments.
Finally, the finiteness of the term cov ( Y 1 , Δ ^ 1 ) follows by the finiteness of E ( Y 1 2 ) , E ( Δ ^ 1 2 ) and Hölder’s inequality.
Remark 4. 
It is useful to note that the width of the constructed CI, and hence the accuracy of the estimation, depends critically on the number of observed regeneration cycles N rather than on the number of observations, as in classical iid estimation based on the CLT. In general, this makes estimation through regenerative simulation more time-consuming, which is expected when dealing with correlated data. However, in all our experiments, a large number of regeneration cycles was obtained in a very limited time.
Remark 5. 
In general, the ratio estimate r, while strongly consistent, is biased. This bias arises from estimating a nonlinear function (the ratio r) based on the ‘scale’ of regeneration cycles. The detailed proof of the interval estimate (49) can be found in [11] (Propositions 4.1 and 4.2).

6. Simulation Results

The aim of this section is two-fold: to highlight the accuracy of the approximations presented in Section 4, focusing on the waiting time distribution for various types of arrival processes, and to verify the applicability of exponential splitting to superposed processes, as described in Section 5.

6.1. Waiting Time Analysis in Presence of Superposed Processes

Initially, we consider queueing systems fed by two independent input processes with Pareto interarrival times. In all the simulations, 10 10 arrivals are generated and exponential service times are assumed.
Using the notation from Section 4.4, we compare the actual arrival process with a renewal process having the same interarrival time distribution, specifically P a r e t o ( x 0 , α 1 + α 2 1 ) . The first set of simulations, carried out with α 1 = 2 and x 0 = 1 , evaluates the accuracy of the approximation by comparing the Kolmogorov distance between the corresponding waiting time distributions for different values of ρ and α 2 . As highlighted in Figure 2, the approximation works better for low utilization values and high α 2 values, i.e., when more moments of the underlying Pareto distribution are finite.
The assumption of common exponential service times allows us to compare the simulation estimates with the exact expression for the distribution function of the stationary waiting time, given by Formula (29). The results are summarized in Table 1, where W, W A , and W T denote the waiting time distributions when the queue is fed by the actual superposed process, its renewal approximation, and the corresponding theoretical approximation (using distribution (29) with the appropriate value of the parameter σ ), respectively. For completeness, the values of μ and σ are also reported. The results in the table confirm the validity of the simulation approach (implemented using ad hoc software written in C) and the trends (in terms of α 2 and ρ ) highlighted in Figure 2.
The distribution functions for the worst and best cases are shown in Figure 3 and Figure 4, respectively. In both cases, W T completely overlaps with W A (refer to the values of d ( W T , W A ) in Table 1). The figures highlight the effect of correlation in the input process on the probability of an empty system, π 0 = 1 σ = F W ( 0 ) , and, more generally, on the shape of the distribution of W.
Apart from the superposition of two Pareto distributions, it is also worth considering the superposition of two processes with Pareto and exponential interarrival time distributions. In these cases, the situation is even better than for the superposition of two Pareto inputs, as highlighted by Table 2, which reports the simulation parameters and the Kolmogorov distance values (as before, 10 10 arrivals have been generated). For the sake of brevity, the waiting time distribution is shown only for the worst case, corresponding to α = 2 , x 0 = 1 , λ = 1 , and μ = 3 (see Figure 5).
To complete the overview of the approximations described in Section 4.4, we also consider the case of class-dependent service times. Table 3 reports the Kolmogorov distance d ( W , W A ) for different values of the service rates μ 1 and μ 2 (and hence of the utilization ρ ). Since the focus is on the impact of the service rate (which is stochastically equivalent to a hyperexponential distribution with parameters μ 1 , μ 2 , and p), we consider two Pareto distributions with the same parameters α 1 = 2 , α 2 = 3 , and x 0 = 1 in all cases. As expected, the values of d ( W , W A ) are in good agreement with those in Figure 2 for α 2 = 3 and the corresponding values of ρ .
Going back to the example at the beginning of this section (see Table 1 for the simulation settings), Table 4 compares the mean values of the actual waiting time W and the waiting time W A of the approximating system. While the Kolmogorov distance, reported in the first column of the table for clarity, was generally less than 5% in most cases, the relative error (RE) can still be significant, consistent with findings in the literature [3,5,7]. This highlights the importance of using such an approximation carefully, as a rough/preliminary estimate when accurate estimation is challenging.

6.2. Exponential Splitting

We consider the original system Σ fed by two Pareto input processes (with α 1 = 2 , α 2 = 3 , and x 0 = 1 ) and a common exponential service rate μ .
We carry out several sets of simulation to evaluate the effect of both the splitting parameters (namely, the choice of the i 0 process, δ and η ) and the traffic intensity ρ .
First, we select the first component of the superposed input (with α 1 = 2 ) as i 0 -process, so the exponential splitting is applied to the other flow (with α 2 = 3 ). In the following tables, we denote by RegIn the number of regenerations in the superposed input (i.e., the number of events (45)), and by RegW the number of regenerations (47) of the queueing system.
Table 5 and Table 6 refer to μ = 4 (i.e., ρ = 0.75 ) and μ = 6 (i.e., ρ = 0.5 ), respectively. In both tables, the estimations of the mean waiting time are the same for all the choices of the splitting parameters and coincide with the results obtained through traditional simulation ( 1.358 and 0.2597 for μ = 4 and μ = 6 , respectively), so they are not explicitly reported. To evaluate the goodness of the estimates, the last column reports the half-widths of the 99 % CI for the mean waiting time, estimated over RegW intervals (corresponding to a fixed number of arrivals, namely 10 10 as already stated).
By comparing Table 5 and Table 6, we see that for fixed values of δ and η , the number of input regeneration points RegIn is approximately the same. This fact is not surprising since the input process does not depend on the service rate μ . On the other hand, RegW increases with μ since a lower traffic intensity implies a higher probability of an idle server. This increment, together with the decrease in the mean waiting time, leads to a difference in the width of the CIs of around one order of magnitude.
Finally, the simulation results in the case of splitting of the first component (with α 1 = 2 ) are reported in Table 7 and Table 8, which confirm the previous observations about RegIn , RegW , and the width of the CIs.
Focusing on the choice of the splitting parameters, the following remarks can be drawn from the simulation results. On the one hand, the larger the δ parameter (for a fixed value of η ), the more frequent the exponential phase. As a result, RegIn also increases. On the other hand, the higher the η parameter (for a fixed δ ), the lower the number of regenerations RegIn . Moreover, simulation results show that the highest number of regenerations RegW of the queueing system is achieved when equality is met in both inequalities (42) for all inputs i i 0 . It is mentioned in [11] that more frequent regenerations normally reduce the required simulation time, but not always. These considerations are intuitive and may be valuable when, for various reasons, the simulation time is limited and the number of regenerations becomes critical for constructing the CI.

7. Conclusions

This paper deals with queueing systems fed by superposed input processes generated by independent stationary renewal sources. It is well known that, in general, such systems are not classically regenerative; however, several analysis methods have been proposed in the literature and we discuss the two most promising ones. First, we describe a general approach that allows for the construction of the so-called one-dependent regeneration. Then, we focus on the superposed process containing heavy-tailed Pareto components and apply exponential splitting to construct classical regenerations of the system. These regenerations are used for confidence estimation of the mean waiting time by the regenerative version of the CLT. This is the main theoretical contribution of the research. Moreover, we identify the conditions under which exponential splitting results in the most frequent regenerations.
We also study an approximating G I / G / 1 system with a single renewal process. This construction is based on Palm theory, and the corresponding stationary interarrival distribution takes the form of a mixture of the interarrival time distributions of the component processes. As a side contribution, in the special case of exponential, class-independent service times, the stationary waiting time distribution in the approximating system is explicitly derived.
Finally, an important element of this research is that the theoretical findings are supported by comprehensive numerical experiments, including both traditional Monte Carlo and regenerative simulation techniques.

Author Contributions

Conceptualization, writing—original draft preparation, E.M. and I.P.; writing—review and supervision, project administration, E.M.; editing, simulation and visualization M.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by the Italian Ministry of Education and Research (MIUR) in the framework of the FoReLab project (Departments of Excellence) and by the University of Pisa in the framework of the PRA_2022_64 “hOlistic Sustainable Management of distributed softWARE systems (OSMWARE)” project.

Data Availability Statement

Data (simulation results) are contained within the article.

Acknowledgments

We would like to thank anonymous referees for their helpful comments. We acknowledge the use of Grammarly (https://app.grammarly.com/) and ChatGPT (https://chat.openai.com/) for English correction and rephrasing in some sections of the manuscript. These tools were not used to generate new content or ideas; they were employed solely to enhance the smoothness and correctness of the existing text.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
FIFOFirst In First Out
iidindependent and identically distributed
r.v.random variable
LSTLaplace–Stiltjes Transform
CLTCentral Limit Theorem
RErelative error
CIconfidence interval

References

  1. Lawrence, A.J. Dependency of Intervals Between Events in Superposition Processes. J. R. Stat. Soc. Ser. B (Methodol.) 1973, 35, 306–315. [Google Scholar] [CrossRef]
  2. Newell, G.F. Approximations for Superposition Arrival Processes in Queues. Manag. Sci. 1984, 30, 623–632. [Google Scholar] [CrossRef]
  3. Whitt, W. Characterizing superposition arrival processes in packet multiplexers for voice and data. J. Sel. Areas Commun. 1986, 4, 833–846. [Google Scholar]
  4. Lam, T.; Lehoczky, J.P. Superposition of Renewal Processes. Adv. Appl. Probab. 1991, 23, 64–85. [Google Scholar]
  5. Albin, S.L. On Poisson Approximations for Superposition Arrival Processes in Queues. Manag. Sci. 1982, 28, 126–137. [Google Scholar] [CrossRef]
  6. Asmussen, S.; Schmidli, H.; Schmidt, V. Tail probabilities for non-standard risk and queueing processes with subexponential jumps. Adv. Appl. Probab. 1999, 31, 422–447. [Google Scholar] [CrossRef]
  7. Whitt, W. Approximating a point process by a renewal process, I: Two basic methods. Operat. Res. 1982, 30, 125–147. [Google Scholar] [CrossRef]
  8. Shedler, G.S. Regenerative Stochastic Simulation; Academic Press Inc.: Cambridge, MA, USA, 1993. [Google Scholar]
  9. Glynn, P.W. Some topics in regenerative steady-state simulation. Acta Appl. Math. 1994, 34, 225–236. [Google Scholar] [CrossRef]
  10. Glynn, P.W.; Iglehart, D.L. Conditions for the applicability of the regenerative method. Manag. Sci. 1993, 39, 1108–1111. [Google Scholar] [CrossRef]
  11. Asmussen, S.; Glynn, P.W. Stochastic Simulation: Algorithms and Analysis; Springer Nature: Cham, Switzerland, 2007. [Google Scholar]
  12. Sigman, K. One-dependent regenerative processes and queues in continuous time. Math. Oper. Res. 1990, 15, 175–189. [Google Scholar] [CrossRef]
  13. Thorisson, H. Coupling, Stationarity, and Regeneration, Probability and Its Applications; Springer: New York, NY, USA, 2000. [Google Scholar]
  14. Khintchine, A.Y. Mathematical Methods in Queueing Theory; Griffin: London, UK, 1960. [Google Scholar]
  15. Whitt, W. Queues with superposition arrival processes in heavy traffic. Stoch. Process. Their Appl. 1985, 21, 81–91. [Google Scholar] [CrossRef]
  16. Smith, W.L. Renewal theory and its ramifications. J. R. Stat. Soc. (Ser. B) 1958, 20, 243–302. [Google Scholar] [CrossRef]
  17. Cox, D.R. Renewal Theory; Methuen: London, UK, 1962. [Google Scholar]
  18. Whitt, W. Comparing counting processes and queues. Adv. Appl. Prob. 1981, 13, 207–220. [Google Scholar] [CrossRef]
  19. Peshkova, I.; Morozov, E.; Maltseva, M. On Comparison of Multiserver Systems with Exponential-Pareto Mixture Distribution. In Computer Networks. CN 2020. Communications in Computer and Information Science; Gaj, P., Gumiński, W., Kwiecień, A., Eds.; Springer: Cham, Switzerland, 2020; Volume 1231. [Google Scholar]
  20. Peshkova, I.V.; Morozov, E.V. On Comparison of Multiserver Systems with Multicomponent Mixture Distributions. J. Math. Sci. 2022, 267, 260–272. [Google Scholar] [CrossRef]
  21. Baccelli, F.; Bremaud, P. Elements of Queueing Theory. Palm Martingale Calculus and Stochastic Recurrences, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2003. [Google Scholar]
  22. Asmussen, S. Applied Probability and Queues, 2nd ed.; Springer: New York, NY, USA, 2003. [Google Scholar]
  23. Smith, W.L. Regenerative stochastic processes. Proc. R. Soc. (Ser. A) 1955, 232, 6–31. [Google Scholar]
  24. Gut, A. Stopped Random Walks Limit Theorems and Applications, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
  25. Serfozo, R.F. Basics of Applied Stochastic Processes; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
  26. Morozov, E.; Steyaert, B. Stability Analysis of Regenerative Queueing Models; Springer: Cham, Switzerland, 2021. [Google Scholar]
  27. Cohen, J.W. The Single-Server Queue, 2nd ed.; North-Holland: Amsterdam, The Netherlands, 1982. [Google Scholar]
  28. Andronov, A. Artificial regeneration points for stochastic simulation of complex systems. In Simulation Technology: Science and Art, 10th European Simulation Symposium; Society for Computer Simulation International: Delft, The Netherlands, 1998; pp. 34–40. [Google Scholar]
  29. Thorisson, H. The queue GI/GI/1: Finite moments of the cycle variables and uniform rates of convergence. Stoch. Process. Their Appl. 1985, 19, 85–99. [Google Scholar] [CrossRef]
Figure 1. Graphical interpretation of inequalities (39).
Figure 1. Graphical interpretation of inequalities (39).
Mathematics 12 02202 g001
Figure 2. Kolmogorov distance between the waiting time distributions of the superposed process and its approximation as a function of ρ and α 2 .
Figure 2. Kolmogorov distance between the waiting time distributions of the superposed process and its approximation as a function of ρ and α 2 .
Mathematics 12 02202 g002
Figure 3. Waiting time distributions of the superposed process W and its approximation W A for α 1 = 2 , α 2 = 2 , x 0 = 1 , and μ = 3 .
Figure 3. Waiting time distributions of the superposed process W and its approximation W A for α 1 = 2 , α 2 = 2 , x 0 = 1 , and μ = 3 .
Mathematics 12 02202 g003
Figure 4. Waiting time distributions of the superposed process W and its approximation W A for α 1 = 2 , α 2 = 10 , x 0 = 1 , and μ = 20 .
Figure 4. Waiting time distributions of the superposed process W and its approximation W A for α 1 = 2 , α 2 = 10 , x 0 = 1 , and μ = 20 .
Mathematics 12 02202 g004
Figure 5. Waiting time distributions of the superposed process W and its approximation W A for α = 2 , x 0 = 1 , λ = 1 , and μ = 3 .
Figure 5. Waiting time distributions of the superposed process W and its approximation W A for α = 2 , x 0 = 1 , λ = 1 , and μ = 3 .
Mathematics 12 02202 g005
Table 1. Kolmogorov distance for different simulation settings (in all cases, α 1 = 2 and x 0 = 1 ).
Table 1. Kolmogorov distance for different simulation settings (in all cases, α 1 = 2 and x 0 = 1 ).
α 2 ρ μ σ d ( W , W A ) d ( W , W T ) d ( W T , W A )
80.8100.81960.023690.023681.307 × 10 5
30.7540.80780.071270.071261.306 ×  10 5
20.6730.76410.082890.082881.030 ×  10 5
30.560.57660.026330.026326.125 ×  10 6
40.580.55870.017570.017565.919 ×  10 6
50.5100.54750.012530.012526.097 ×  10 6
100.5200.52440.004010.004005.044 ×  10 6
30.3100.36550.009490.009492.680 ×  10 6
Table 2. Kolmogorov distance for different simulation settings.
Table 2. Kolmogorov distance for different simulation settings.
α x 0 λ μ σ d ( W , W A ) d ( W , W T ) d ( W T , W A )
21130.70670.0640.0649.84 ×  10 6
21140.54800.0320.0325.97 ×  10 6
310.540.68790.0230.0237.56 ×  10 6
24130.42520.0010.0062.73 ×  10 6
Table 3. Kolmogorov distance for different simulation settings with class-dependent service time.
Table 3. Kolmogorov distance for different simulation settings with class-dependent service time.
μ 1 μ 2 ρ d ( W , W A )
340.8330.1178
450.6500.0540
460.5830.0454
470.5360.0403
480.5000.0370
490.4720.0350
670.4520.0226
680.4170.0201
690.3890.0184
Table 4. Comparison of the mean waiting times for the simulation settings in Table 1.
Table 4. Comparison of the mean waiting times for the simulation settings in Table 1.
d ( W , W A ) mean ( W ) mean ( W A ) RE (%)
0.023690.496600.454208.54
0.071271.358301.0506922.65
0.082891.484571.0795427.28
0.026330.259720.2269312.63
0.017570.173730.158228.93
0.012530.129580.121016.62
0.004010.056420.055122.29
0.009490.061930.057596.99
Table 5. Regenerative simulation for i 0 = 1 and μ = 4 ( ρ = 0.75 ).
Table 5. Regenerative simulation for i 0 = 1 and μ = 4 ( ρ = 0.75 ).
δ η RegInRegW99%CI
0.104166.65  × 10 6 13.5 ×  10 6 4.59 ×  10 4
0.334555.61 ×  10 6 45.1 ×  10 6 4.55 ×  10 4
0.504833.40 ×  10 6 67.5 ×  10 6 4.58 ×  10 4
0.7541250.05 ×  10 6 101.3 ×  10 6 4.58 ×  10 4
0.105133.29 ×  10 6 9.1 ×  10 6 4.34 ×  10 4
0.335444.41 ×  10 6 30.4 ×  10 6 4.56 ×  10 4
0.505666.67 ×  10 6 45.6 ×  10 6 4.55 ×  10 4
0.106111.11 ×  10 6 6.5 ×  10 6 4.46 ×  10 4
Table 6. Regenerative simulation for i 0 = 1 and μ = 6 ( ρ = 0.5 ).
Table 6. Regenerative simulation for i 0 = 1 and μ = 6 ( ρ = 0.5 ).
δ η RegInRegW99%CI
0.104166.67 ×  10 6 42.7 ×  10 6 4.29 ×  10 5
0.334555.61 ×  10 6 142. 4 ×  10 6 4.59 ×  10 5
0.504833.32 ×  10 6 213. 6 ×  10 6 4.20 ×  10 5
0.7541250.06 ×  10 6 320.5 ×  10 6 4.57 ×  10 5
0.105133.34 ×  10 6 29.8 ×  10 6 4.47 ×  10 5
0.335444.41 ×  10 6 99.5 ×  10 6 4.27 ×  10 5
0.505666.60 ×  10 6 149.3 ×  10 6 4.55 ×  10 5
0.106111.10 ×  10 6 22.1 ×  10 6 4.58 ×  10 5
Table 7. Regenerative simulation for i 0 = 2 and μ = 4 ( ρ = 0.75 ).
Table 7. Regenerative simulation for i 0 = 2 and μ = 4 ( ρ = 0.75 ).
δ η RegInRegW99%CI
0.5031111.11 ×  10 6 88.7 ×  10 6 4.56 ×  10 4
0.104166.65 ×  10 6 11.1 ×  10 6 4.51 ×  10 4
0.334555.57 ×  10 6 37.1 ×  10 6 4.61 ×  10 4
0.504833.35 ×  10 6 55.6 ×  10 6 4.55 ×  10 4
0.106111.11 ×  10 6 5.6 ×  10 6 4.54 ×  10 4
0.206222.18 ×  10 6 11.3 ×  10 6 4.60 ×  10 4
0.336370.41 ×  10 6 18.8 ×  10 6 4.56 ×  10 4
Table 8. Regenerative simulation for i 0 = 2 and μ = 6 ( ρ = 0.5 ).
Table 8. Regenerative simulation for i 0 = 2 and μ = 6 ( ρ = 0.5 ).
δ η RegInRegW99%CI
0.5031111.14 ×  10 6 282.5 ×  10 6 4.52 ×  10 5
0.104166.65 ×  10 6 36.7 ×  10 6 4.41 ×  10 5
0.334555.53 ×  10 6 122.4 ×  10 6 4.23 ×  10 5
0.504833.36 ×  10 6 183.6 ×  10 6 4.18 ×  10 5
0.106111.09 ×  10 6 19.5 ×  10 6 3.84 ×  10 5
0.206222.21 ×  10 6 39.0 ×  10 6 4.41 ×  10 5
0.336370.35 ×  10 6 65.0 ×  10 6 4.43 ×  10 5
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Peshkova, I.; Morozov, E.; Pagano, M. Regenerative Analysis and Approximation of Queueing Systems with Superposed Input Processes. Mathematics 2024, 12, 2202. https://doi.org/10.3390/math12142202

AMA Style

Peshkova I, Morozov E, Pagano M. Regenerative Analysis and Approximation of Queueing Systems with Superposed Input Processes. Mathematics. 2024; 12(14):2202. https://doi.org/10.3390/math12142202

Chicago/Turabian Style

Peshkova, Irina, Evsey Morozov, and Michele Pagano. 2024. "Regenerative Analysis and Approximation of Queueing Systems with Superposed Input Processes" Mathematics 12, no. 14: 2202. https://doi.org/10.3390/math12142202

APA Style

Peshkova, I., Morozov, E., & Pagano, M. (2024). Regenerative Analysis and Approximation of Queueing Systems with Superposed Input Processes. Mathematics, 12(14), 2202. https://doi.org/10.3390/math12142202

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop