Study on Transient Queue-Size Distribution in the Finite-Buffer Model with Batch Arrivals and Multiple Vacation Policy

The transient behavior of the finite-buffer queueing model with batch arrivals and generally distributed repeated vacations is analyzed. Such a system has potential applications in modeling the functioning of production systems, computer and telecommunication networks with energy saving mechanism based on cyclic monitoring the queue state (Internet of Things, wireless sensors networks, etc.). Identifying renewal moments in the evolution of the system and applying continuous total probability law, a system of Volterra-type integral equations for the time-dependent queue-size distribution, conditioned by the initial buffer state, is derived. A compact-form solution for the corresponding system written for Laplace transforms is obtained using an algebraic approach based on Korolyuk’s potential method. An illustrative numerical example presenting the impact of the service rate, arrival rate, initial buffer state and single vacation duration on the queue-size distribution is attached as well.


Introduction
Mathematical models of queueing systems with limited capacity of buffers accumulating customers waiting for service are widely used in practice. Among many different areas of their application, it is worth mentioning issues related to modeling the functioning of production systems, problems arising in transport issues or logistics. However, particular applications of systems with finite waiting rooms have been found in modeling the functioning of telecommunication and computer network nodes. Indeed, phenomena, such as the merging of many independent streams of incoming packets (packet queueing), packet processing (service) but also possible overflows of accumulating buffers and packet losses because of too much disproportion between the intensity of the input traffic and the speed of processing, occur there. Dynamically developing the branch of wireless communication (Wi-Fi, LTE, sensor networks, etc.) forces the search for solutions that, in addition to modeling the packet processing itself, ensure the implementation of tasks resulting from the need for energy saving. In the field of modeling the energy saving mode in wireless network nodes, queueing systems with various types of mechanisms of the periodic suspension of customer processing are used. Among many solutions, it is worth mentioning the N-policy threshold-type discipline, originally proposed in [1], in which the processing restarts after an idle period simultaneously with the Nth arrival, and the T-policy, introduced in [2], in which the server is being activated exactly T time units after the last busy period. Moreover, in a fundamental paper [3], single and multiple vacation policies were proposed in which a single or repeated independent vacations are taken until at least one customer present in the buffer is detected. The results obtained in [1,3] for the M/G/1-type queue were next generalized in [4,5] for the batch Poisson arrival stream and for the M/M/c queue, respectively. An overview of results obtained by using the time and single vacation duration are generally distributed. In this article, we study the time-dependent queue-size distribution in the M X /G/1/N-type queueing model with a multiple vacation policy and an exhaustive service. Applying a theoretical approach based on identifying Markovian moments in the evolution of the system, we construct a Volterratype system of integral equations for the transient queue-size distribution conditioned by the buffer filling level at starting time t = 0. For the corresponding system written for Laplace transforms (LTs for short) we use the idea of Korolyuk's potential (see [39]) to obtain the solution in the compact form. So, the main contribution of the paper is the closed-form representation for the LT of the conditional queue-size distribution in the non-stationary state of the system. The formula is written by means of a certain functional sequence (in [39] called a potential) with terms, which can be obtained recursively or with the help of power expansion. Numerical examples illustrate analytical results, namely the impact of the vacation duration, arrival intensity, service speed, buffer size and the initial system state on the behavior of the queue-size distribution.
The remaining part of the article is organized as follows. In Section 2, we provide a detailed description of the considered queueing model and state the necessary nomenclature. In Section 3, we obtain the system of integral equations for the conditional queue-size distribution in the transient state. In Section 4, we sketch the idea of Korolyuk's potential approach and use it to obtain the compact-form solution of the corresponding system written for LTs. Section 5 is devoted to numerical computations illustrating theoretical results. To visualize the appropriate probabilities we use one of the algorithms of numerical LT inversion. Moreover, comparing numerical results to the simulation is conducted. Finally, Section 6 contains a summary and short conclusions.

Mathematical Description of the Model
We deal with the M X /G/1/N-type queueing model in which groups (batches) of packets arrive according to a compound Poisson process with intensity λ. Group sizes are generally-distributed random variables and p k stands for the probability that the arriving group consists of k packets exactly, where ∑ ∞ k=1 p k = 1. The processing is organized according to the natural FIFO discipline and the service time of an individual packet is a random variable with a general-type cumulative distribution function (CDF for short) F(·). The system capacity equals N,, that is, there is an accumulating buffer with N − 1 places for waiting packets and one place reserved for the packet being processed. The considered queueing system is governed by the multiple vacation policy. Namely, every time the server becomes idle (at the completion epoch of each busy period) the server takes a vacation of a random duration with a general CDF G(·). After the vacation, the buffer state is examined. If there are no packets waiting for service, the next independent vacation is started and so on. The last vacation is the one at the completion epoch of which at least one packet waiting for processing is detected. Successive server vacations started as long as the buffer is empty are called a multiple vacation period (MVP for short). So, the operation of the system can be observed in successive busy periods, during which the processing of packets is uninterrupted, followed by MVPs during which the processing is suspended.
We end this section by introducing necessary notations. We denote by f (·) and g(·) Laplace-Stieltjes transforms (LST for short) of CDFs F(·) and G(·), respectively. Similarly, tails of CDFs F(·) and G(·) are denoted by F(·) and G(·), respectively. We define the i-fold Stieltjes transform of G(·) with itself as follows: where t > 0 and i ≥ 1. Similarly, p i * j stands for the jth term of the i-fold convolution of the sequence (p k ) with itself, namely: where i ≥ 1, and I{·} is the indicator defined as : Finally, let us denote by X(t) the number of packets present in the system at time t (including the one being processed at this time, if any).

Equations for Conditional Queue-Size Distribution
In this section, identifying Markovian moments in the evolution of the system, we build a system of integral equations for the transient queue-size distribution conditioned by the buffer state at the starting time t = 0. We introduce the following notation: where t > 0 and 0 ≤ m, n ≤ N. Let us consider firstly the case of the system being empty at the starting epoch t = 0. In such a case we assume that an MVP begins at this moment, so we identify the starting time t = 0 with the completion of a busy period. Note that the continuous total probability law leads to the following integral equation (having in mind that λ > 0): Let us comment briefly on the above formula. The moment of the first packet arrival is denoted by y, while the last single vacation before time t (distributed with a CDF G(·)) completes at time u. Two components in the first square brackets on the right side of (5) correspond to the situation where the first MVP completes before time t (the completion epoch of the MVP is denoted by u + v). The first summand relates to the case in which there is at least one place free in the accumulating buffer at time u + v (this moment also coincides with the starting moment of the busy period). So, if k ≤ N − 1 packets enter at time y, then next, at most N − k − 1 packets arrive during a time period of length u + v − y. The second summand describes the case of the buffer being saturated at the starting time of the busy period. The second square brackets on the right side of (5) represent the situation in which the MVP ends after time t, however, at least one packet arrives before t. If m ≤ N − 1 then the number of arrivals up to the moment t equals m exactly (the first summand). In the case m = N, the number of arrivals must be greater than or equal to N (the second summand).
Finally, the expression I{m = 0}e −λt corresponds to the case of no arrivals before time t. In this case, X(t) = m if and only if m = 0. Consider now the case where the buffer contains at least one packet at the starting moment t = 0. As it is well known, (X(t n +)), where t 1 < t 2 < . . . , are successive service completion (departure) epochs, is an embedded Markov chain (EMC) in the M X /G/1-type queue. The continuous version of the total probability law applied with respect to the first departure epoch after the starting moment t = 0 leads to the following system of integral equations: where λ > 0 stands for the group arrival intensity and n ∈ {1, . . . , N}. Indeed, the integral on the right side of (6) relates to the situation in which the first packet leaves the system at time y < t : the first summand in brackets corresponds to the case where there is at least one place before that moment in the accumulating buffer; in the second summand, the buffer is saturated just before the first departure, so after the first service completion epoch the state of the system equals N − 1. The other two components after the integral on the right side of (6) represent the situation in which the first service ends after time t.
Let us introduce the LT of P n (t, m) as: and define for s > 0 the following functional sequences: where Note that the system (5)-(6) leads now to the following corresponding one written for LTs: and where n ∈ {1, . . . , N}.
We will now rewrite the last system of Equations (16) and (17) in a specific form that will be relevant in the context of the method that we will use to obtain the solution. By introducing the following notation: Equation (17) will take the form: where n ∈ {0, . . . , N − 1} and Similarly, Equation (16) can be rewritten as:

Main Analytical Result
In this section, we use Korolyuk's potential approach to obtain the solution of the system (19) and (21) in a compact form. In [39], the following linear system with an infinite number of equations is considered: where (α n ) and (ψ n ) are two known number sequences, and (q n ) is the sequence of unknowns. Such a system has infinitely many solutions; however, as was proved in [39], each solution of (22) can be written as follows: where L ∈ R and, in general, it can be different for different solutions, and the sequence (R k ), called a potential in [39], is defined by the sequence of coefficients (α n ) in two equivalent ways. Firstly, its successive terms can be obtained recursively, one by one, in the following way: where k ≥ 2. Alternatively, successive terms of the sequence (R k ) can be obtained by means of power expansion. Indeed, defining: one can easily show that In consequence, expanding the right side of (26) in the Maclaurin series and next comparing coefficients at successive powers of z on the left and right sides of (26), we obtain: Note that the systems of Equations (19) and (22) have the same form. However, there are two main differences. Firstly, functions in the system (22) depend on some arguments: we have α k = α k (s), q k = q k (s, m) and ψ n = ψ n (s, m). Secondly, the number of equations in (19) compared to (22) is finite. To obtain the representation for the solution of (19), we can use the result (23) assuming now that L = L(s, m) and R k = R k (s). As it turns out, Equation (21) written for n = N will, however, serve as the boundary condition allowing for the explicit determination of L(s, m). We prove the following main theorem: In the M X /G/1/N-type queueing model with multiple vacation policy the representation for the Laplace transform P n (s, m) of the conditional transient queue-size distribution is the following: (29) and the formulae for α k (s), β k (s), γ(s), δ(s, m), k (s, m) are given in (8), (9), (10), (11) and (12), respectively, and the sequence R k (s) is defined recursively as follows (compare (24)): where k ≥ 1.
To obtain the explicit representation for Q n (s, m) we need the formula for L(s, m). Taking n = 0 in (32), we have: Similarly, substituting n = 0 into (19), we get: where (see (20)) Due to the fact that ∑ ∞ k=0 α k (s) = f (s) we have from (34) Using now (33) and (35) in (32), we obtain where functions Λ n (s) and χ n (s, m) are defined in (29) and (30), respectively. Now let us use Equation (21) as a specific boundary condition in order to find the formulae for Q 0 (s, m). Applying (36) in (21) gives: Now, referring to (18) and (36), completes the proof.

Numerical Study
In this section, we investigate numerically the impact of the main "input" system parameters (like arrival intensity, service rate or mean vacation duration) on the queuesize distribution. To obtain the transient queue-size distribution from the formula (28) for particular system parameters we use the Gaver-Stehfest algorithm of the numerical Laplace transform inversion described in detail in [40,41] (see also [42,43]). In computations, we deal with a model described as follows: • Poisson arrivals of packets of sizes 100 B; • CDF of the processing time is a mixture of two exponential distributions (second-order hyper-exponential distribution) and is defined as: where α, µ 1 and µ 2 are given; • CDF of a single server vacation has 2-Erlang distribution with fixed parameter ξ, namely: • buffer size equals N = 10 packets (so 1 kB). A hyper-exponential probability distribution of single packet's processing time can be used in modeling the situation in which a single packet (customer, task, message, etc.) can receive two types of service: normal (typical) and special, for which the average duration times differ significantly. The value of parameter α indicates the fraction of incoming packets, which receive typical service that takes on average µ −1 1 time units. The other packets are processed according to a special mode.
Since during the server idle mode (multiple vacation period) the service station can perform other tasks (such as e.g., background tasks of the CPU processor), the 2-Erlang distribution can be used to approximate a single vacation period duration: we can assume that it consists of two consecutive subperiods of random length. During the first of them, the task is prepared to run in the background, and its direct execution takes place in the second period.
The result obtained in Theorem 1 is of a general nature (i.e., CDFs F(·), G(·) and the sequence (p k ) are general). However, in numerical computations we use the fact that the right side of (28) is a rational function of variable s. If P n (s, m) is not a rational function of s we can use, for example, the classical technique of Padé approximation of a given Laplace transform by an appropriate rational function (see, e.g., [44]).
Hence, µ is equal to 800, 600 and 480 packets/s that correspond to processing rates 640, 480 and 384 kb/s, respectively. Keeping constant the arrival intensity λ = 600 packets/s, we have the traffic load of the system = 0.75, 1 and 1.25, respectively. Assume that the distribution of the arriving group size is given by the sequence p k = {0.5, 0.5, 0, . . .} so the mean group size equals 1.5. Moreover, let ξ = 1000 so the mean single vacation duration equals 0.002 s. In Figure 1 Observe that the probability P{X(t) = 1 | X(0) = 0} is the highest just after the opening of the system. This is due to the fact that, after starting, the server immediately goes into the vacation mode due to the lack of packets to handle.

Impact of Arrival Rate
Analyze now the impact of the arrival rate λ on the transient queue-size distribution. Consider three scenarios in which λ = 450, 600 and 750 packets/s, which corresponds to data transfer rates of 360, 480 and 600 kb/s. The parameters of the CDF of the processing time are α = 0.6, µ 1 = 680 and µ 2 = 510, which gives the mean service rate µ = 600. The values of the traffic load for successive scenarios are 0.75, 1 and 1.25, respectively. Moreover, as previously, ξ = 1000 so the average single vacation duration equals 0.002 s. Taking (p k ) = {0.5, 0.25, 0.25, 0, . . .}, in Figure 2 probabilities P{X(t) = 1 | X(0) = 0} are presented.
As time passes, the probability decreases until the system reaches equilibrium. The nature of the decreasing probability depends on the arrival intensity: the higher the λ, the faster the system reaches equilibrium.

Impact of Initial Buffer State
Investigate now the impact of the initial buffer state n on the transient probability P{X(t) = 1 | X(0) = n} for n = 0, 5 and 10. Assume λ = 450 packets/s, which corresponds to the arrival rate 600 kb/s. Moreover, take the service rate µ = 600 packets/s so the offered traffic load equals 0.75. Taking the same parameter ξ = 1000 of the single vacation duration CDF and the sequence (p k ) = {0.5, 0.25, 0.25, 0, . . .} describing sizes of arriving groups, the transient behavior of the considered probability is visualized in Figure 3. As can be seen, for small values of time t, the influence of the initial buffer state is clearly visible and the differences between corresponding probabilities are relatively large.

Impact of Single Vacation Duration
Finally, let us check the response of the transient queue-size distribution to changes in the average length of a single vacation. Consider three different scenarios in which the arrival intensity λ takes, as previously, values 450, 600 and 750 packets/s and the service rate µ = 600 packets/s. In consequence, = 0.75, 1 and 1.25, respectively. In each scenario we analyze the behavior of the transient probability P{X(t) = 1 | X(0) = 0} for three different values of parameter ξ = λv of the single vacation duration CDF, namely 500, 1000 and 1500 (corresponding to mean single vacation durations 0.004, 0.002 and 0.0013 s). The size of an arriving group of packets is distributed according to sequence Let us comment briefly on the final results. As one can observe, the longer the mean single vacation duration (the smaller the value of parameter ξ = λv), the smaller the probability P{X(t) = 1 | X(0) = 0}. Indeed, longer vacations are conducive to the accumulation of more packets in the buffer before the start of the processing. This phenomenon is particularly visible at the lowest traffic load value = 0.75.

Numerical Computations vs. Simulation Results
In this subsection, we compare numerical results with the simulation model created in Objective Modular Network Testbed in C++ (OMNeT++) [45]. We modify the standard FIFO model imbedded in OMNeT++ discrete event simulator in a way to transmit and queue batches of packets according to the following scenarios: • The transient probability P{X(t) = 1 | X(0) = 0} is calculated for λ = µ = 600 packets/s (resulting in = 1). We take ξ = 1000 as a parameter of the single vacation duration CDF (corresponding to mean single vacation duration 0.002 s) and the sequence (p k ) = {0.5, 0.25, 0.25, 0, . . .} describing sizes of arriving groups (so the mean group size equals 1.75); • As in the previous scenario except for the sequence (p k ) = {1, 0, 0, . . .} (single arrivals).
The simulation model is extended by a time step of t c = 0.0005 s, describing the time moments k · t c = t ∈ [0, 0.06] for k = 0, 1, 2, . . . , 60, at which the number of packets present in the system is measured. The simulation run is divided into 10,000 individual trials. Each attempt is a separate time period t ∈ [0, 0.06] in which the random processing is simulated in accordance with the multiple vacation policy assumed in our queueing model. A unique seed is randomly selected for each sample and used to initialize the pseudorandom number generator representing random arrival and service as well as the lengths of individual vacation periods. Figures 7 and 8 confirm the compliance of the numerical calculations with the simulation results. The maximum difference in the probability values between the numerical calculations and the simulation results measured at times k · t c equals 0.00599106. (2) as a statistical result of a 10,000th random sample-the red squares in the diagram.

Conclusions
In this article, transient queue-size distribution in the batch-arrival system with multiple vacation policy is studied. For the purposes of the analysis, the approach based on the idea of the embedded Markov chain, continuous formula of total probability, integral equations and linear algebra is proposed. As the main result, the closed-form formula for the Laplace transform of the transient queue-size distribution conditioned by the initial buffer state is obtained. In the numerical study, the impact of the service rate, arrival intensity, initial buffer state and mean single vacation duration on the studied characteristic is investigated for examples of the considered queueing model. The considered queueing model has many potential applications in modeling the operation of production systems, computer and telecommunication networks with an energy saving mechanism based on cyclic monitoring the state of the queue of packets (messages, tasks, etc.) waiting for service. During the multiple vacation period, the service station may perform other types of tasks, such as, for example, the time the CPU is running in the background, installing software updates or maintaining the server. In applications in modeling the operation of an automated production line, during downtime, a single machine (service station) can be used in secondary production; during this time, the machine can also be calibrated or the settings can be changed. Due to the fact that, in some practical applications, interarrival times are not exponentially distributed, it is planned in future to extend the analysis for the case of a non-Poisson arrival stream (e.g., hyper-exponential or a general renewal one).