Abstract

We consider a single-server discrete-time queueing system with N sources, where each source is modelled as a correlated Markovian customer arrival process, and the customer service times are generally distributed. We focus on the analysis of the number of customers in the queue, the amount of work in the queue, and the customer delay. For each of these quantities, we will derive an expression for their steady-state probability generating function, and from these results, we derive closed-form expressions for key performance measures such as their mean value, variance, and tail distribution. A lot of emphasis is put on finding closed-form expressions for these quantities that reduce all numerical calculations to an absolute minimum.

1. Introduction

In this paper, we analyze a queueing model with generally distributed customer service times and a correlated arrival process. This kind of model is, for instance, useful in assessing the performance of a packet-switched communication infrastructure, where messages carried by the network can have a variable transmission time, such as the current Internet, or a dedicated packet-based network carrying video-on-demand (VoD). Multiplexer queues and/or switches and routers in such a network can in general be modelled by means of a discrete-time queueing system, where new packets (i.e., “customers”) are generated by a superposition of individual traffic sources. The size of the packets is proportional to their transmission times, whose distribution depends on the specific application under consideration. Analyzing such a system is mandatory in the design and evaluation of these networks. However, this can be a difficult task because of the time-correlated behaviour that each of the individual sources may exhibit.

To facilitate the queueing analysis, a source is often modelled as a discrete-time Markov modulated arrival process, such as the discrete-time batch Markovian arrival process (D-BMAP; [1, 2]) with the Markov modulated Bernoulli process (MMBP; [35]) and the switched batch bernoulli process (SBBP; [6, 7]) as frequently used special cases. As a consequence, the analysis of a queueing model with a specific SBBP [8] or MMBP [4, 911] source model or with a general D-BMAP (e.g., [1, 1221], and the vast amount of references therein) has become the focus of many research papers. In addition to the characteristics of the arrival process, such queueing models can also be classified by means of the service process. For customer service times equal to one slot, single-server queueing models are considered in [4, 9, 12, 15], while [18, 21] focus on the multiserver case. For the single-server case, geometrically [1, 10] and phase-type [13] distributed customer service times have been considered as well, while these quantities are generally distributed in [8, 11, 14, 17, 19, 20]. In this paper, we investigate the scenario where arrivals are generated by an aggregation of sources, where each source is characterized by means of a (two-state) D-BMAP, and the sources are assumed to be homogeneous, that is, they are all described by the same stochastic process. The service times on the other hand are assumed to be generally distributed.

Although the queueing model that is studied is not entirely new, the novelty of this contribution lies in the way by which the problem is tackled and solved. A widely used method to deal with correlated arrival processes is the matrix-analytic solution technique [22, 23], which provides a broad and generic framework for the analysis of a wide set of queuing systems, both for finite- [4, 912, 14, 15, 19] and infinite-capacity [8, 24] queues. This is based on an intelligent matrix representation of the quantities that are needed in these analyses and making full use of the properties that these matrices exhibit. One impediment it may nevertheless sometimes suffer from is the size of the matrices that must be manipulated, which is (roughly speaking) determined by the state space of the system. An alternative solution technique, adopted in [1, 1618], is the so-called spectral decomposition solution method; however, the main results that are generated still contain operations that can be computationally demanding, such as Krönecker sums and products, or inversions, of matrices that can have large dimensions. The approach that we adopt—based on the work presented in [21] for single-slot service times—is related to the latter one, albeit that we follow a different path by invoking as much as possible a representation by means of (joint) probability generating functions (pgfs) in our derivations presented in the remainder of this paper. The main advantage lies in the observation that the formulae for the performance indices that we finally obtain, such as the moments of the queue contents, unfinished work, and customer delay, are—to a large extent—expressed directly in terms of the system parameters and therefore often lend themselves to an interpretation that is helpful for understanding the mechanisms that determine the queue behaviour. In addition, in order to establish (semi)analytic expressions for the relevant performance measures that pertain to the queue behaviour, a (possibly large) set of boundary probabilities has to be calculated numerically. To reduce the numerical work as much as possible, we introduce approximations for these quantities that allow us to compute the aforementioned performance measures without great difficulty. By means of some numerical examples we will demonstrate the efficiency and accurateness of these approximations, as well as the influence of the system parameters on the queue behaviour.

2. System Description

We consider a discrete-time queueing model, that is, a system where the time axis is divided into fixed-length contiguous intervals, referred to as slots. The system consists of one single server and an infinite waiting room for customers awaiting service. Customers arrive in the queue according to a correlated Markovian arrival process, described further on. It is assumed that the service of a customer requires a positive integer number of slots, described by a general service time distribution, and can start (and end) at slot marks only, that is, at the time points between consecutive slots, meaning that the customer service is synchronized with respect to the slot marks. This is not necessarily the case as far as customer arrivals are concerned: customers may enter the queue at any continuous time instant. Nevertheless, due to the synchronous nature of a customer’s service, it can commence at the earliest beginning of the slot following its arrival (if it arrives in an empty queue). Therefore, the precise details of the position of the customer arrival instants within a slot are irrelevant, and it suffices to characterize the arrival process by a random variable describing the total number of customer arrivals during a slot. We will now discuss the customer arrival and service processes in more detail.

Consider a buffer model fed by identical and independent sources generating customers with variable service times. Each source is modelled as a 2-state Markov modulated arrival process, where the state of a source during a slot is represented by , . The customer arrival process during a slot is completely characterized by the matrix State transitions occur at slot boundaries, and let us denote by , , the one-step transition probability at the end of slot , where of course . The elements of are then given by where , , is the probability generating function (pgf) describing the number of customers generated during a slot by a source, given that the source is in state during the tagged slot and was in state during the preceding slot, that is, where , , represents the number of customer arrivals generated by source during slot and denotes the expected value of the tagged quantity. For the Markov modulated Bernoulli process (MMBP), the number of customers generated by a source during any slot is either zero or equal to one, meaning that each of the generating functions is a linear function of , that is, for some parameters satisfying . Although attention is focused on this specific arrival model in the numerical examples, the analysis is general and can also be applied when the ’s have a more complex form. Moreover, the notational conventions have been chosen such that the derivations and results that are presented in the subsequent sections can be readily extended to an -state D-BMAP as well (see also [20]), . We confine ourselves to the case of in this paper in order to enhance the readability of the derivations that follow hereafter, while their level of complexity is that of arbitrary .

The aggregate customer arrival process is fully determined, once the probability generating matrix has been specified. Let us define , as the number of sources in states during slot . Note that, due to the fact that the total number of sources equals , these random variables satisfy In the rest of the paper, we will denote by the column vector (where represents the matrix transposition operation), and similarly we write . Let us also denote by the column vector the matrix product and, for a set of random variables , define the joint generating function as Then, with the previous definitions, it is not difficult to show that the joint probability generating function of the random variables denoting the number of customer arrivals and the state of the arrival process during a slot can be written in terms of the state of the arrival process during the previous slot, leading to where represents the total number of customer arrivals during slot , that is, Equation (2.8) fully describes the number of customers generated during consecutive slots by the customer sources.

Starting from some initial state, the customer arrival process will evolve to a stochastic equilibrium after a sufficiently long period of time, and we define as the joint probability generating function of the number of sources in states and during an arbitrary slot. If we define and as the steady-state probability that the Markov state of a source during an arbitrary slot is and , respectively, and as the vector , which is the solution of the matrix equation leading to where is the column vector with all elements equal to 1, then equals It is further assumed that the service times of consecutive customers that arrive in the queue form a set of independent and identically distributed random variables, denoted by , with common probability mass function and corresponding probability generating function Obviously, we assume that the service time of a customer is at least one slot.

Due to the infinite-queue capacity, the system will reach a stochastic equilibrium only if the equilibrium condition is satisfied. Defining as the load of the system, this implies that must hold, where denotes the mean number of customer arrivals per slot and per source, and represents the mean customer service time. Taking into account the customer arrival process described in this section, it follows that can be calculated from where primes denote derivatives with respect to the argument.

3. The Joint pgf of the State Vector

3.1. System Equations

We first establish the system equations that control the evolution of the number of customers in the queue during consecutive slots. This is inspired by the work presented in [25], where the case of an i.i.d. (i.e., uncorrelated) customer arrival process was considered and analysed by means of the so-called supplementary variable method. Let us therefore define the random variable as the system contents at the beginning of slot , which is the number of customers in the queue at the beginning of slot , including the one being served (in case of a nonempty queue). In addition to , we also need information about the amount of time the customer in service still requires at the beginning of slot before being completely served. We therefore define the random variable as the residual service time, which is the number of slots required at the beginning of slot to complete the service of the customer being served in case of a nonempty queue; when , we automatically have .

Together with the definitions of the previous section, we can now establish the system equations. We must distinguish between the three following cases.(1). When the queue is empty at the beginning of a slot, the number of customers in the queue at the beginning of the next slot equals the number of new arrivals, and we find that If no customers have entered the queue during slot , the residual service time remains zero, otherwise it is equal to the service time of a new customer, which leads to (2). This implies that the customer in service is completely served at the end of slot : Similar to the previous case, if no customers have entered the queue during slot and the customer in service was the only one in the queue, then the queue becomes empty and the residual service time equals zero at the beginning of the next slot, otherwise the service of a new customer commences, that is, (3). In this case, the customer being served receives an extra slot of service without the service being completed, and we obtain

3.2. Derivation of a Functional Equation

Clearly, in view of (3.1)–(3.5), we need to keep track of the random variables and in our state description. Consequently, due to the correlated first-order Markov nature of the process controlling the number of customer arrivals during a slot, it becomes clear that the set of random variables constitutes a four-dimensional Markov state description of the system at the beginning of consecutive slots; this set of four random variables will be referred to as the state vector. Note that, in view of (2.5), one of the random variables can be omitted from the state description. However, due to reasons of symmetry, we prefer to maintain both random variables in the state description.

Let us now define the joint probability generating function of as Let us, for the sake of notation, define a partial generating function of a random variable if an event occurs as From system equations (3.1) and (3.2), taking into account that , we then first of all derive that In a similar way, we obtain from (3.3) and (3.4) Finally, (3.5) can also be translated into a relation between -transforms, leading to Summation of (3.8)–(3.10) yields Note that, given that is known, the value of (and the number of customer arrivals ) is independent of . Therefore, applying (2.8), the previous relation can be transformed into where First of all, due to the foregoing definitions, we have that . If we insert into (3.12), we derive that Let us now assume that the equilibrium condition (2.15) is satisfied, implying that all the functions that occur in (3.12) evolve to a steady-state limit, which will be reflected in the remainder by suppressing the subscript (which expressed the time dependence in the previous derivations). Together with the previous relation, we finally obtain the functional equation Equation (3.15), together with definition (3.13), defines a functional equation for , the steady-state joint generating function of the random variables , that is, the system contents and residual service time at the beginning of an arbitrary slot, and the state of the Markovian customer arrival process during the preceding slot.

3.3. Solution of the Functional Equation

In order to solve the functional equation, we follow a similar approach to that in [21]. First, if we substitute consecutive arguments , into (3.15), we obtain We now let approach infinity. Let us therefore express in terms of its eigenvalues and eigenvectors as . For further details concerning the eigenvalues and eigenvectors of , we refer to the appendix. Define the functions ,  ,  , as In the above expression, the functions are the components of the column eigenvectors of , given by (A.5). One can easily derive that this relation implies that these functions can be written as where and . In view of definition (3.17) and defining for , the first sum in the right-hand side of (3.16) can be written as where and constitute the th eigenvalue and corresponding row eigenvector of . The Perron-Fröbenius (P-F) eigenvalue is denoted by (see the appendix for further details). Due to definition (3.17), this can be transformed into Let us now consider values of and for which . Such values of and exist; this is for instance the case if , and and (e.g., see [26]) since unless some special cases for the generating functions defined in (2.3) are considered—for example, the case that the number of customer arrivals per slot is always a multiple of —a situation that falls beyond the scope of this paper. The sum for in the above expression then converges, and we obtain In a similar way we can derive that where In addition, note that which, for , becomes equal to 0. Summarizing, in view of the previous results, (3.16) becomes where has been defined as Expression (3.27) for the joint probability generating function still contains the unknown constants , as well as the unknown functions . Let us therefore, for each value of , consider values of for which , which will be represented by . Note that if , due to (3.22) this solution for also satisfies , and, for these values of and , is finite. This means that, if we multiply both hand sides of (3.27) by and consider , then the left-hand side becomes equal to zero, implying that the same must hold for the right-hand side. For each value of , this leads to the following equation: If we insert this equation into (3.27), we finally obtain This final result for expresses the joint generating function of the state of the arrival process during an arbitrary slot, and the queue contents and residual service time at the beginning of the next slot, in terms of the unknown constants . In the next section, it will become clear how these can be calculated or approximated. Also, starting from this result, we will establish expressions for the generating functions of quantities such as the queue contents, unfinished work, and message delay, and their corresponding moments and tail distribution.

3.4. Calculation of the Unknowns

We can calculate these probabilities, defined in (3.19), by expressing that is analytic when the complex variables and both lie inside the unit circle, that is, and . We have already expressed in the previous section that, for and is finite. In view of expression (3.30), the only remaining potential singularities inside the complex unit disk are those values of (and ) for which . Indeed, taking into account the definition for , then from (3.22) and Rouché’s theorem it follows that this equation has exactly one solution inside the complex unit circle. Let us, for each value of , denote this solution by , that is, Then, by expressing that must also be a zero of the corresponding numerator in expression (3.30) for , we obtain where Note that (since , see (A.6), implying ), which leads to no additional equation for the unknowns. Together with the normalization condition, which yields (as we will see further on), (3.32) forms a set of linear equations for the unknowns that can easily be solved.

However, in deriving numerical results for the performance measures considered further in this paper, we want to avoid all numerical calculations as much as possible. This makes sense in the particular case that is large, that is, when a high number of sources are multiplexed, since the number of equations to be solved equals . Therefore we also present an approximation for the boundary probabilities . As before, denote by the random variables and the number of customer arrivals during a tagged arbitrary slot and the queue contents at the beginning of the following slot; also, represents the state of the Markovian arrival process during the tagged slot. We define the steady-state joint probability Obviously, implies that there have been no customer arrivals during the tagged slot, that is, . It is therefore clear that the inequality holds. We will show by some numerical examples that approximating the boundary probabilities by yields excellent approximations for the performance measures of interest. Here, we have taken into account that so that substituting by this formula in (3.30) still yields a normalized generating function.

An additional advantage of the approximation proposed here is that it avoids the explicit calculation of the ’s, as well as the functions from (3.18), that occur in expression (3.30) for . From definition (3.17) and expression (2.8) it is not difficult to show that where represents the th column eigenvector of . In view of the previous, this implies that the following approximation may be applied: To conclude, note that it immediately follows from (2.8) that satisfies

4. The Queue Contents

4.1. Derivation of the Probability Generating Function

First of all, due to property (A.2), we have that . Therefore, inserting and into (3.30), we obtain the following expression for the probability generating function describing the number of customers in the queue at the beginning of an arbitrary slot: From (A.7) (see the appendix) and (2.15), it follows that =, implying that the normalization condition for indeed requires that , that is, . From (4.1), one can easily derive the performance measures of interest related to the queue contents, as shown next.

4.2. Moments of the Queue Contents

The mean and variance of the queue contents can be calculated by taking the appropriate derivatives of expression (4.1) for . First, note that (A.9) implies that (where represents the discrete Krönecker-delta function, which equals 1 if is 0 and 0 otherwise). Since this, in turn, implies that , and defining , we obtain the following expression for mean and variance of the queue contents: where, with (A.7) and (2.15), In these expressions, primes denote derivatives with respect to the argument, as before. From definition (3.33) for and the observation that , we also find that If we choose to use approximation (3.37) for in order to avoid the numerical calculation of the ’s, this becomes In addition, the term in the expression for is given by which can be considered as the first derivative with respect to of for , where, using (3.22) and the definition for and , satisfies where we have also taken into account that . Taking the first-order derivative, we thus obtain, keeping in mind that and (see (A.9)-(A.10)) By using expressions (4.1)–(4.8), we have now expressed the moments of the queue contents exclusively in terms of the system parameters and the derivatives of the P-F eigenvalue and eigenvector given in (A.7)–(A.11).

4.3. Tail Distribution of the Queue Contents

It has been observed in many cases (see, e.g., [16, 20, 2629]) that approximating the tail distribution of the queue contents by a geometric form may be quite accurate if the dominant singularity (i.e., the singularity with the smallest modulus) of is a pole; as shown in [20], this will be the case if where is the minimum of the radii of convergence of the pgf’s . If each of the pgf’s and either is analytic in the entire complex plane or has a dominant singularity that is a pole, then corresponds to a pole as well and the above inequality is always satisfied, since will then tend to for . Only if corresponds to a branch point, it may be so that this inequality is not satisfied; this calls for a case-by-case investigation of the corresponding tail behaviour of the queue contents and falls beyond the scope of the current paper.

Approximating the distribution of the queue contents by a geometric form corresponds to approximating by where we are particularly interested in sufficiently large values of . In this expression, is the pole with the smallest modulus of , which is a zero of the denominator for in expression (4.1) for , which yields and is a real and positive number larger than 1. We assume that this is a pole with multiplicity 1; note that the above formula can be readily extended to the case where has multiplicity larger than 1 should such a scenario occur. Furthermore, the constant can be calculated from the residue theorem, leading to Combining these relations, we find the following approximate expression for the probability that the queue contents exceed a certain threshold : Again, we can either use the exact value for calculated from (3.33) or the approximation calculated from (3.37).

5. The Unfinished Work

5.1. Derivation of the Probability Generating Function

The unfinished work in the queue at the beginning of a slot is defined as the amount of work in the queue at the beginning of the slot, which is the number of slots required to empty the queue if no more customers were to arrive during the subsequent slots. The amount of work required to process a customer whose service has not started yet is given by its service time described by . On the other hand, at the beginning of a slot, the customer currently being served (if any) still requires slots before it is entirely sent, where equals the residual service time. The unfinished work at the beginning of slot is therefore given by where each of the ’s represents a customer service time, and is therefore described by the pgf . Consequently, the corresponding steady-state pgf is given by The observation that can be written as together with the previous relations, leads to the following expression for the steady-state probability generating function describing the unfinished work at the beginning of an arbitrary slot: Again, due to and , note that , as expected. We would like to point out that each of the denominators in the above expression also has a zero inside the unit disk , which leads to an additional linear equation for the boundary probabilities (since is bounded inside the complex unit circle ), as explained in Section 3.4. However, note that if is a solution of with , we then have that , implying that where satisfies (3.35), and we obtain exactly the same set of linear equations as in (3.37) for the boundary probabilities.

5.2. Moments of the Unfinished Work

In a similar way to that in Section 4.2, if we define and take the appropriate derivatives with respect to of expression (5.4) for , we find an expression for the mean and variance of the unfinished work: where The derivatives of and in the previous expressions are given by (4.4) (or (4.5) in case the approximations for the boundary probabilities are applied) and (4.8), respectively.

5.3. Tail Distribution of the Unfinished Work

Under similar assumptions to those outlined in Section 4.3, we can adopt a geometric tail approximation for the distribution of the unfinished work. We then obtain an approximation for the probability that the unfinished work exceeds a threshold : where the smallest pole of is the solution of , which is a real and positive quantity larger than 1, and, similar to (5.5), can also be calculated from . Also note that , implying that the above formula can be expressed in terms of the smallest pole of as well.

6. The Customer Delay

6.1. Derivation of the Probability Generating Function

We define the customer delay as the number of slots an arbitrary customer remains in the queue, including its service time, that is, the number of slots between the end of its arrival and departure slot, as depicted in Figure 1. Consider an arbitrary tagged customer that enters the queue during a slot, say slot . Note that this slot is not an arbitrary slot but a slot where an arbitrary customer arrives (which for instance implies that at least one customer arrives). Due to the FCFS service discipline, the delay of a customer is determined by the amount of work in the queue upon the tagged customer’s arrival. In particular, if we denote by the unfinished work at the beginning of slot and by the number of customers that have arrived during the same slot as the tagged one and will be served before it, then we find that the delay of the tagged customer is given by where the random variables represent the service times of the tagged customer and the customers that have arrived during its arrival slot and will be served before it, which are all described by the pgf . Let us denote by the steady-state joint generating function of the random variables describing the unfinished work at the beginning and the number of customer arrivals during an arbitrary slot , that is, and by the steady state joint pgf of . Then, similar to [30], we can establish the following relation: System equation (6.1) can be transformed into a relation between -transforms, which together with (6.3), leads to an expression for the steady-state pgf describing the delay of an arbitrary customer On the other hand, since the unfinished work satisfies the system equation and assuming that a stochastic equilibrium is reached, this implies that the pgf (given by (5.4)) can also be written as (again, all ’s are described by ) In combination with the previous expression for (and making use of ), this finally leads to where in view of This formula for , which expresses the delay of an arbitrary customer in terms of the unfinished work at the beginning of an arbitrary slot and which apparently is independent of the details of the customer arrival process (correlation, etc.) was also derived in [31] via an alternative analysis. From this result, it is not difficult to express the performance measures concerning the customer delay in terms of the corresponding performance measures for the unfinished work, as shown next.

6.2. Moments of the Customer Delay

By taking the first two order derivatives of (6.7)-(6.8) with respect to for , we can establish the following relations between the mean and variance of the customer delay and the mean and variance of the unfinished work: Through some simple calculations, one can verify that Little’s theorem indeed holds, that is, .

6.3. Tail Distribution of the Customer Delay

Again, applying an analogous technique as the one used for determining the tail distribution of the queue contents and unfinished work, the asymptotic behavior of the customer delay tail distribution can be calculated from where is the smallest pole of which, in view of (6.7)-(6.8), leads to the observation that , and the probability that the customer delay exceeds a threshold is approximated by

7. Some Numerical Examples

We will illustrate the results derived throughout the previous sections by a few numerical examples. Let us therefore consider a two-state ON/OFF MMBP arrival process for each source, which means that that is, customers are generated at a rate during an ON period (corresponding with state in the modulating Markov chain) and no customers are generated during an OFF period (corresponding with state ). Let us also denote by and the mean length of an ON and OFF period, respectively, that is, and . Furthermore, we will assume that customers have a constant service time equal to slots, implying that . The load of the system is then given by and any subset containing five parameters of the set fully describes the customer arrival process in the queue. Such a model has for instance been used to effectively model a voice-over-IP (VoIP) packet arrival stream in a network node [3].

In the remainder we consider a multiplexer queue fed by identical sources and focus on the customer delay behaviour, which is the primary performance measure for VoIP applications. Similar results can be readily obtained for the queue contents and unfinished work. For the arrival model described above, we have plotted the mean and variance of the the moments of the customer waiting time (which is defined as the number of slots a customer waits in the queue before its service commences and is therefore equal to the customer delay minus the customer service time) versus in case of and for increasing customer service time in Figures 2 and 4, respectively. Varying the value of the load in this case means that the arrival rate is adjusted. From these figures we can conclude that the moments of the customer waiting time strongly depend on the customer service time and increase for increasing values of . In these figures, we have also included the curves when approximation (3.35) is used for the boundary probabilities (dashed lines), which leads to approximation (4.5) for the derivatives of ; the absolute error between the exact and approximate results is so small that the latter are hardly visible in Figures 2 and 4. Therefore, in Figures 3 and 5, we have also plotted the relative error (RE) that is made when using the approximate formula for the mean and variance of the waiting time. Apparently, except for , the RE of mean and variance appears to be a decreasing function of that approaches zero for . In addition, the larger the customer service times are, the smaller the relative error generally tends to be. We also observed that, while for the mean value the approximation always forms an upper bound compared to the exact value, this is not the case as far as the variance is concerned.

In Figure 6, we consider constant values for and and increasing values of the mean length of an active period. Hence, the probability that a source is in the ON state during an arbitrary slot tends to 0 for . Clearly, the mean value of the customer waiting time, which is plotted versus in these figures, strongly depends on the mean length of an active period. In addition, we observe from these curves that the error that is made by using the approximation for the boundary probabilities (dashed line) can be significant when the load is low and the arrival rate and are high. Another feature is that the mean waiting time not necessarily tends to 0 for . Indeed, even under these conditions, we sporadically encounter the situation where a source becomes active and generates customers at a rate during a geometrically distributed ON period with mean length , and the intersection with the ordinate designates the mean waiting time for customers generated by such a single source.

In Figures 7, 8 and 9, we have plotted the tail distribution of the customer waiting time for a wide range of system parameters. Each curve for the exact value of the parameters of the geometric tail behaviour (full line) is accompanied by the corresponding approximation (dashed line) obtained by using approximation (3.37) for . We may conclude that, for the wide range of system parameters considered in these figures, the approximate results are extremely accurate. Furthermore, it is observed that the approximation forms an upper bound, and, consequently, is on the safe side. The only numerical calculation required for this upper bound is that of the dominant pole of the pgf describing the steady-state waiting time distribution. Insertion of the system parameters in the corresponding formulas then immediately yields these results.

In Figure 7, we have plotted the probability that the customer waiting time exceeds a certain threshold versus , for a multiplexer queue with sources, a system load , an arrival rate during an ON period, and a mean length of an ON period , for varying values of the customer service time . Again, as was also observed when discussing the moments, for these parameter values the results for the quantile of the waiting time (defined as the smallest value of for which ) significantly depend on the customer service time. From this figure, it also becomes apparent that the quantile of the waiting time is roughly proportional to , and a wide range of case studies has revealed that this observation is not limited to the case of constant customer service times.

In Figure 8, we show the results for the customer waiting time in case of , , and and increasing values of the load (which means that the arrival rate during the ON periods increases), while the load is fixed to and in Figure 9. As can be expected, the slope of the customer waiting time tail distribution increases as both and increase. We also should point out that, whereas the accuracy of the approximate results for the moments of the customer waiting time—calculated by invoking the approximation for the boundary probabilities—deteriorates for increasing values of (and ) (e.g., Figure 6), this is not so as far as the tail distribution of the customer waiting time is concerned, and the approximate data that we thus compute remain highly accurate.

8. Conclusion

In this paper we have investigated and analysed the discrete-time (2-state) D-BMAP/G/1 queueing model. By combining a solution method that was previously developed for the D-BMAP/D/c queueing system, with the supplementary variable technique that allows us to cope with generally distributed customer service times, we were able to derive an expression for the steady-state pgfs of the queue contents, unfinished work, and customer delay. This has led to (semi)closed-form expressions for a number of interesting performance measures such as the mean value, variance, and tail distribution of these quantities. In addition to the analysis of the system, a lot of emphasis was put on finding closed-form expressions for these quantities that reduce all numerical calculations to an absolute minimum, at the expense of introducing approximations for a number of boundary probabilities that must normally be calculated by solving a—potentially large—set of linear equations. The only numerical computation that needs to be carried out is that of a dominant pole, situated on the positive real axis, when one calculates tail probabilities, and the complexity does not depend on the size of the state space of the model under consideration. To conclude, the accuracy of this approach, as well as the impact of the system parameters on the customer delay performance, was evaluated through a considerable set of numerical examples.

Appendix

The Eigenvalues and Eigenvectors of

As it became clear from the analysis throughout this report, an important role is played by the eigenvalues and eigenvectors of . Let us denote by the diagonal matrix containing the eigenvalues and by and the matrices containing the corresponding right column and left row eigenvectors, respectively, that is, with the diagonal unity matrix. The normalization constants of the left row and right column eigenvectors are chosen such that We will denote by () the right column (left row) eigenvector corresponding to .

Taking into account expression (2.1) for , the characteristic equation of this matrix can be written as Solution of the characteristic equation of leads to In addition, solving (2.16) and (3.1) yields the following expression for : In a similar way, we can also calculate the elements of ; however, the explicit expression of this matrix is not required in our derivations. Note that is the Perron-Fröbenius (P-F) eigenvalue of and is an important quantity for the analysis in this paper. Similarly, will be referred to as the P-F eigenvector of .

In our expressions for the mean and variance of the queue contents, unfinished work and customer delay, we also require the first-, second-, and third-order derivatives with respect to for of the P-F eigenvalue. Taking the appropriate derivatives of the characteristic equation (A.3), we obtain with (A.6) where was defined in (2.16). The th derivative of with respect to for is readily calculated in terms of the th derivative of from the relation

Similarly, we also require the first- and second-order derivatives with respect to for of the P-F eigenvector. First of all, we point out that, due to (A.5), satisfies Then, (Strictly speaking, the second column of is no longer an eigenvector of for . However, it is not difficult to verify that all the results derived throughout this paper still hold when taking the limit .) taking the appropriate derivatives of the equations generated by (A.2) and the relation , we obtain where

Disclosure

J. Walraevens and D. Fiems are Postdoctoral Fellows of the Research Foundation-Flanders (F.W.O.).