OF A DISCRETE-TIME QUEUE WITH GENERAL SERVICE DEMANDS AND PHASE-TYPE SERVICE CAPACITIES

In this paper, we analyze a non-classical discrete-time queueing model where customers demand variable amounts of work from a server that is able to perform this work at a varying rate. The service demands of the customers are integer numbers of work units. They are assumed to be independent and identically distributed (i.i.d.) random variables. The service capacities, i.e., the numbers of work units that the server can process in the consecutive slots, are also assumed to be i.i.d. and their common probability generating function (pgf) is assumed to be rational. New customers arrive in the queueing system according to a general independent arrival process. For this queueing model we present an analysis method, which is based on complex contour integration. Expressions are obtained for the pgfs, the mean values and the tail probabilities of the customer delay and the system content in steady state. The analysis is illustrated by means of some numerical examples.

1. Introduction.In many queueing applications, where some kind of customers require some kind of service from a given service facility (or "server"), the amounts of service that the customers demand are different, and the rate at which the server is able to provide that service also varies over time.In some applications, the service rate varies between only 2 possible values, corresponding to "on" and "off" states of the server.This is typically modeled using the concept of vacations or service interruptions (see e.g.[13,14,26,27]).In other applications, where the server is not always either on or off but may also be working at a reduced (or accelerated) speed, the service rate can take on a range of possible values.Example applications where both the service demands and the service rates vary over a range of possible values include packet-switched routers where the packet sizes are variable and the available bandwidth fluctuates over time due to network congestion or the time-varying or fading nature of wireless channels (see e.g.[5,10,19,20]), manufacturing plants where the production capacity varies due to maintenance actions, gradual deterioration or failures of machines (see e.g.[4,15]), web services where the available processing power fluctuates due to background processes or shared hosting, etc.These variable service requirements and variable service rates are traditionally modeled using the single notion of "service time", which is the amount of time that the server needs to fully serve one customer (see e.g.[2,17,30]).This single notion, however, is not always sufficient to model both effects.This is because if the amount of service that each customer requires (which we refer to as the "service demand" of the customer) varies from customer to customer, and the amount of work that the server is able to perform per time unit (which we refer to as the "service capacity" of the server) varies over time, then there may be a non-trivial correlation between the consecutive service times.Indeed, the service time of a customer depends on the service capacity of the server during the service of that customer, which in turn may be correlated with the service capacity during the service of the next customer, which in turn influences the service time of that next customer.The degree of correlation between consecutive service times may depend on many system parameters, including the load and the arrival process in general.Indeed, if the load is very low, then there will most often be long idle times between the service of one customer and the next, so there will likely be little or no correlation between the service times.On the contrary, if the load is high, the amount of correlation may be very large.These kinds of effects cannot be modeled by the many classical queueing models that take the arrival process and the service times to be independent of each other.
In the scientific literature some papers do exist on continuous-time queuing models where the variable service demands of customers and the variable service rates of the server are modeled explicitly (see e.g.[6,18,21,24,28]).However, the discretetime equivalents with variable service demands and variable service capacities have received only little attention so far.
In this paper, we therefore study a discrete-time queueing model that explicitly models both the service demands of the customers and the service capacity of the server.Specifically, in our queueing model, time is divided into fixed-length intervals, referred to as (time) slots, and both the service demands of the customers and the service capacities in each time slot are assumed to be integer numbers of "work units".The service demands of the customers are assumed to be independent and identically distributed (i.i.d.) from customer to customer, and the service capacities are i.i.d.from slot to slot.Finally, also the numbers of customer arrivals during the consecutive slots are assumed to be i.i.d.random variables.
With a focus on discrete-time queueing models, existing related work considered various restrictions for the distribution of the service capacities.In [8] and [31], it was assumed that the service capacities follow a geometric distribution.In [9], the model was analyzed under the restriction that the service capacities are deterministically equal to a given constant.In [32], Yao et al. obtained a relationship between the customer-delay and the system-content distributions for the specific case of constant service demands and constant service capacities.In [11] we analyzed this model with the restriction that the distribution of the service capacities has finite support.Finally, in [12], we considered the case where the probability generating function (pgf) of the service capacities is a rational function; the analysis in [12] was, however, restricted to the steady-state customer delay.
In our present paper, we now again consider service capacities with a rational pgf and we extend the analysis of [12] to include the system content as well.Note that all phase-type distributions have a rational pgf (see e.g.[23]), and note in particular that all the restrictions on the service capacities in previous work imply that the service capacities follow a phase-type distribution.Therefore, our present paper can be seen as a generalization of the papers [8,31,9,11,12].
The paper is organized as follows.In Section 2, we describe the considered queueing model in more detail.Next, in Section 3, we outline the analysis of the queueing model and present expressions for the pgfs of the unfinished work in the system, the delay of an arbitrary customer, and the system content in steady state.The calculation of moments of the system content and the delay is discussed in Section 4, and approximations for the tail probabilities of the system content and the delay are derived in Section 5. Section 6 is devoted to a remarkable property, referred to as the "invariance property".We discuss some numerical examples in Section 7 and conclusions are given in Section 8.

2.
Queueing model description.In this paper, we study a discrete-time queueing model, where time is divided into contiguous fixed-length intervals, referred to as (time) slots.New customers arrive in the queueing system according to a general independent arrival process.This means that the numbers of arriving customers during the consecutive slots are i.i.d.from slot to slot.We denote pgf of the number of customers arriving in an arbitrary slot by A(z).The mean number of customers arriving per slot, the so-called mean arrival rate, is denoted by λ A (1).
Each customer has a service demand, expressed as a positive integer number of work units.This is exactly the amount of work that the server will have to perform, possibly over the course of multiple time slots, to completely serve the customer.The service demands of the customers are assumed to be i.i.d.from customer to customer.The common pgf of the service demands of the customers is denoted by S(z).The mean service demand is denoted by τ S (1).
The number of work units that the server can execute in a time slot is referred to as the service capacity of the server during that time slot.These service capacities are assumed to be non-negative integers that are i.i.d.from slot to slot.We denote their common pgf by R(z), which is assumed to be a rational function.The mean service capacity (per slot) is denoted by µ R (1).We also introduce the mutually prime polynomials P R (z) and Q R (z) such that and we let m denote the degree of Q R (z).We assume that the numbers of arrivals in each slot, the service demands of the customers, and the service capacities during each slot, are mutually independent random variables.
The server cannot initiate the service of a customer during the arrival slot of that customer.Stated otherwise, the service of a customer can start at the earliest during the slot following his arrival slot, even if the customer arrives in an empty system.Customers from the queue are served sequentially by the server in first-come-firstserved (FCFS) order.In each slot, no more work units are executed than the available service capacity for that slot.If the available service capacity during a slot is less than the (remaining) service demand of the customer currently in service, then that customer's service simply continues in the next slot, with a reduced remaining service demand.Conversely, if the service capacity is larger than the remaining service demand of the customer currently in service, the server will completely serve that customer and use its remaining service capacity to immediately (i.e., still during the same slot) start also the service of the next customer in the queue (if any).This is repeated until either the whole service capacity of the server during that slot has been used or there are no customers left in the queue that still require service.
3. Queueing analysis.In this section, we present the analysis of the above queueing model.Specifically, we derive expressions for the steady-state pgfs of the unfinished work in the system, the customer delay and the system content.Additional details on the analysis method are given in the appendices.
3.1.Pgf of the unfinished work.As a first step in the analysis, we derive an expression for the pgf U (z) of the unfinished work, i.e., the sum of the (remaining) service demands of all the customers in the system, at the beginning of an arbitrary slot in steady state.This pgf is a useful intermediate result for the delay analysis of Section 3.2.
We begin our derivation by introducing a notation for some of the random variables pertaining to the state of the system in slot k.We let U k denote the unfinished work at the beginning of slot k, we define A k as the number of customer arrivals during slot k, the variable S k,i represents the service demand of the ith (i = 1, 2, ..., A k ) customer entering the system during slot k and R k is the service capacity during slot k.In every slot k, the following system equation then holds between these random variables: where (...) + = max(..., 0).Let us now assume that the system is stable, i.e, that the equilibrium condition λτ < µ is satisfied.Taking the z-transform of both sides of equation (2), taking the limit for k → ∞ and using the fact that all the above random variables pertaining to the same slot k (i.e., U k , A k , R k and the S k,i 's) are all independent of each other, we then easily obtain In Appendix A, we show that under the assumption that the service-capacity pgf R(z) is a rational function (see (1)), equation (3) can be further transformed into the following expression for the pgf U (z): where S −1 R denotes the set of poles of R(1/z), µ ξ denotes the multiplicity of a pole ξ ∈ S −1 R , N − T denotes the set of zeros of T (z) Q R (z) − A(S(z))P R (z) inside or on the unit circle, excluding the zero at z = 1, and n ξ denotes the multiplicity of a zero ξ ∈ N − T .
3.2.Pgf of the customer delay.We now turn to the analysis of the customer delay.More specifically, we derive in this section an expression for the pgf D(z) of the delay D C that an arbitrary customer C experiences in the system in steady state, under a FCFS scheduling discipline.This delay is measured as the number of slots between the end of the arrival slot of the customer and the end of the slot during which the customer leaves.Note that the customer delay cannot be 0, since a customer that arrives during a slot cannot receive any service during that same slot.Before proceeding with the actual delay analysis, we first derive the pgf of a related quantity, V C , the unfinished work observed by the customer C upon arrival.It is defined as the total number of work units present in the system just after the arrival slot of customer C, but to be executed before or during the service of customer C. Work units belonging to any customer arriving after customer C are not counted, even if some of those work units are executed while customer C is technically still in the system, i.e., during the last slot of the service of customer C. Mathematically, V C is therefore defined as where J denotes the arrival slot of customer C, F C is the number of customers that arrive in slot J but are to be served before C, and S J,i is the service demand of the ith customer in slot J.It is well-known (see e.g.[25]) that for any queue with independent, ordered arrivals, the pgf of F C is given by Using this property and equation ( 3), the pgf V (z) of V C then follows immediately as • (A(S(z)) − 1) The delay D C of customer C is related to the quantity V C as follows.Customer C will still be in the system at the start of a slot if and only if fewer than V C work units have been executed since the end of the arrival slot J of C. In other words, In Appendix B, the relationship (9) between random variables D C and V C is transformed into the following relationship between their pgfs: where L is a contour around the origin such that ∀ξ ∈ L : Next, under the assumption of a rational pgf R(1/z), it is shown in Appendix B that the relation (10) further simplifies to where the functions α k (z) are the m zeros for ξ of The relationship (11) is valid for all z for which the zeros α k (z) are distinct, which is the case for all but at most 2m − 1 values of z.
Substituting the expression (8) that we previously found for the pgf V (z), and using R(1/α k (z)) = 1/z (see (12)) to simplify the result, we finally obtain This expression still contains the functions α k (z).These are the zeros of an mdegree polynomial with coefficients that depend on z, and for these zeros a closedform solution is only available for specific classes of distributions (see Appendix D).If no closed-form solution is available for these zeros, then inverting the pgf D(z) analytically is very difficult.However, inverting this pgf numerically, using methods such as those described in [1], is always very straightforward.Additionally, expression ( 13) can be used to calculate the expected value and other moments of the delay, as detailed in Section 4, and may also be used to derive a dominant-pole approximation for the tail probabilities of the delay, as explained in Section 5.

Pgf of the system content.
The last quantity that we study in this paper is the system content.Specifically, in this section, we obtain an expression for the pgf B(z) of the system content, i.e., the total number of customers present in the queueing system, at the beginning of an arbitrary slot in steady state.
To do so, we first relate the distribution of the system content at the beginning of an arbitrary slot to the system-content distributions at arrival and departure instants.It is a well-known fact (see e.g.[22]) that in any system with ordered arrivals and departures, the pgf B a (z) of the system content just before the arrival of an arbitrary customer is equal to the pgf B d (z) of the system content just after the departure of an arbitrary customer, i.e., B a (z) = B d (z).To relate B(z) and B a (z), we simply use (6) and find Secondly, we relate the system content just after the departure of an arbitrary customer C to the delay D C of that customer.Since customers are served in FCFS order, the system content just after the departure of C is given by where G C is the number of customers arriving in the same slot as customer C, but to be served after C, and A C,k is the number of customers arriving in the kth slot after the arrival slot of customer C. Note that G C and D C are conditionally independent when given F C .Also note that similar to (6), the joint distribution of G C and F C can be related to the distribution of the number of customers A arriving in an arbitrary slot in steady state, as follows: Therefore, we can derive the pgf B d (z) of the system content at departure instants as where D f (z) denotes the conditional pgf of the delay of an arbitrary customer C, given that F C takes the integer value f .Next, using similar methods as before, the conditional pgf D f (z) and, in view of the intermediate results ( 14) and ( 17), also the pgf B(z) can be related to the pgf U (z) of the unfinished work at the beginning of an arbitrary slot.In particular, in Appendix C, the following relationship is obtained: where L is a contour around the origin such that ∀ξ ∈ L : Moreover, it is shown in Appendix C that for a rational pgf R(z) of the service capacities, the contour integral in this expression can be further simplified to where the functions β k (z) are defined as α k (A(z)), k = 0, 1, ..., m − 1, i.e., the zeros for ξ of The relationship (19) is valid for all z for which these zeros β k (z) are distinct, which is the case for all z but an isolated set; see Appendix D for further discussion on these functions β k (z).Finally, combining equation (19) with the result (4) for U (z), we obtain the following expression for B(z): 4. Moments of the system content and the customer delay.In this section, we derive explicit expressions for the first-order moments, i.e., the expected values, of the system content at the beginning of an arbitrary slot and the delay of an arbitrary customer in steady state, and we explain how to obtain expressions for higher-order moments of these quantities.
An expression for the mean system content can be obtained by evaluating the first derivative of B(z) to z at z = 1.This is however not entirely trivial, as it requires the determination of the derivatives of the implicitly defined functions β k (z) at z = 1.These latter derivatives can be calculated by differentiating both sides of equation (20), which yields Higher-order derivatives of β k (z) may be obtained by differentiating both sides of equation ( 22) repeatedly.
Evaluation of the derivative of B(z) at z = 1 also requires the use of l'Hôpital's rule, since it can be seen from equation (20) that one of the zeros β k (z) equals 1 when z = 1.This zero must have multiplicity 1, since R (1) > 0, so we may unambiguously denote this zero as β 0 (z), so that β 0 (1) = 1.Under a few assumptions which we will discuss shortly, we may now obtain an expression for the mean system content by differentiating (21) to z, substituting z = 1, and using l'Hôpital's rule on the term for k = 0.This yields The assumptions that must hold are the following: • All zeros β k (1) must be distinct.If they are not, then equation ( 21) is not applicable for z = 1.The mean value of the system content may then be obtained by evaluating the derivative with respect to z at z = 1 of both sides of equation (72) instead.However, since this is relatively difficult in general, it may be easier to approximate B (1) numerically by using a finite difference method in that case.• The term for k = 0 must be the only term of equation ( 21) that requires l'Hôpitals rule for the evaluation of the derivative at z = 1.It can be shown that the assumption does not hold if and only if the greatest common divisor g of the period of the service demands and the period of the service capacities is greater than 1.In that case, the simple division of all service demands and capacities by g yields a queueing system that behaves exactly the same way, and in particular it has the same moments of the system content and the delay.For this equivalent queueing system, g will be equal to 1, so the expression (23) (applied to the equivalent queueing system) may be used to calculate the mean system content.The mean value D (1) of the customer delay may be obtained similarly by evaluating the derivative of equation ( 13) at z = 1.The obtained expression is the same as (23) multiplied by a factor 1/λ, in agreement with Little's law.
By using the moment-generating property of probability generating functions, also higher-order moments of the system content and delay may be obtained.This requires the evaluation of higher-order derivatives of B(z) or D(z) at z = 1.The calculations, e.g. of the variances of the system content and the delay, are rather tedious and the resulting expressions are too long to be included here.We do, however, show some numerical results in Section 7.

5.
Tail probabilities of the system content and the customer delay.In this section, we describe dominant-pole approximations for the tail probabilities of the system content and the customer delay in steady state.As explained in [7], if z X is the smallest positive real-valued pole of the pgf X(z) of a random variable X, then the tail probabilities of X can be approximated as and where C X is the residue of X(z) at z = z X .These approximations thus require the determination of 2 constants: z X and C X .We will describe how to determine these constants, first for the system content and then for the customer delay.
5.1.System content.In general, the dominant pole z B of the system content has to be found numerically.One way of doing this is by finding the smallest positive real-valued zero of B(1/z).Since this is a real-valued zero, many methods can be used for this, such as the bisection method or the Illinois algorithm.
If it is known a priori that the zeros β k (z B ) are distinct, e.g. if the service-capacity distribution belongs to one of the classes of distributions described in Appendix D, then it is not necessary to use the lengthy expressions for B(z) obtained in Section 3.3.Instead, the following theorem is available for a much faster determination of z B : Proof.From expression (21) for B(z) it is clear that z B must be one of the following: Therefore, A(z B ) would have to be 0.But by ( 14), (17) and the fact that D f (0) = 0, B(z B ) would then have to be 0 as well, but z B is supposed to be a pole of B(z).Once z B is known, C B can be calculated.We again assume that the zeros β k (z B ) are distinct.It is easy to check numerically which terms in the summation over k in (21) contribute to the pole of B(z) at z = z B , i.e., for which values of k z B = S(β k (z B )).We denote the set of these k as K.The residue C B of B(z) at z = z B can now be calculated by summing the individual residues of these terms in (21).We find where we also used the fact that if k ∈ K, then S(β k (z B )) = z B to simplify the result.

Customer delay.
For the delay, the dominant pole z D and the residue C D can be found as follows.Determining z D could in principle be done similarly to determining z B by finding the smallest positive real-valued zero of 1/D(z).However, in Appendix B we showed that the radius of convergence R D of D(z) is given by R D = 1/R(1/R V ).This implies that z D = 1/R(1/z V ), where z V is the dominant pole of V (z).From ( 7) and (4) it is easy to see that the dominant pole z V of V (z) is given by the smallest positive real-valued zero of 1−R(1/z)A(S(z)).Since this does not contain any implicitly defined functions, finding z V and using it to calculate z D is typically easier than finding z D directly.
To determine the residue C D , we will use the relation (11).Note that the factor α k (z) − 1 in the denominator of (11) does not yield a pole of D(z) at z = z D since that would imply z D = 1.Then, using arguments similar to those used in Theorem 5.1, we find that if the zeros α k (z D ) are distinct, then the factor R (1/α k (z)) does not contribute to the dominant pole either, so that α k (z D ) must be a pole of V (z) for some k.Denote the set of all z that are both a pole of V (z) and equal to α k (z D ) for some k ∈ [0, m − 1] as V.Then, we can calculate the residue C D as where we have used the facts that if 6. Invariance property.In [8] it was first observed that when the service demands follow a shifted geometric distribution (with minimum value 1) and the service capacities follow a geometric distribution, i.e., when and then the distributions of the system content and the delay depend on τ and µ only through their ratio τ /µ.This remarkable property was termed the "geometric invariance property" in [8].From numerical results obtained for the more general queueing model considered in this paper (see also further in Section 7), we discovered that this may also be true for some other service-demand and service-capacity distributions.In this section, we therefore first explore other conditions under which the "invariance property" holds, i.e., under which the distributions of the system content and the delay only depend on τ and µ through their ratio τ /µ.If the service demands follow a shifted geometric distribution, then due to the memoryless property of the geometric distribution, each work unit of a customer's service demand has the same probability 1/τ of being the last work unit of that customer's demand, regardless of how many other work units of service that customer has already received.Therefore, if the system has a certain service capacity of n work units in a given slot, then the number of customers that can leave in that slot is binomially distributed with parameters n and 1/τ .Hence, the pgf R(z) of the number of customers that can leave in an arbitrary slot k can be calculated as To determine when the invariance property holds, let c denote the ratio τ /µ.Then (30) can be rewritten as Since the pgf R(z) describes the number of customers that can be served each slot, which is independent from slot to slot, this pgf fully characterizes the service process.Therefore, if the service-capacity distribution is such that (31) only depends on c, and not on the individual value of µ, then the invariance property holds.This is at least the case for the following distributions: 1.The negative binomial distribution, with pgf with as special case the geometric distribution (m = 1).

The binomial distribution, with pgf
with as special case the Bernoulli distribution (m = 1).3. The Poisson distribution, with pgf R(z) = exp(µ(z − 1)).Next, if the service capacities follow a geometric distribution, then again due to the memoryless property of the geometric distribution, regardless of how many work units of service have already been executed during a slot, there will be service capacity remaining in that slot with probability µ/(1 + µ) and no remaining capacity with probability 1/(1 + µ).Therefore, the number of slots it takes to serve each work unit of a customer's service demand simply follows a geometric distribution with parameter 1/(1 + µ) and mean 1/µ.Hence, if a customer has a certain service demand n, then his "service time" is negative binomially distributed with parameters 1/(1 + µ) and n and mean n/µ.Note that this service time is independent from one customer to the next, and can be equal to 0 slots.This allows us to calculate the pgf S(z) of the service time of an arbitrary customer C as Again denoting the fraction τ /µ by c, we rewrite this as Similarly as before, the pgf S(z) fully characterizes the service process, so that if the service-demand distribution is such that expression (35) does not depend on τ but only on c, then the behavior of the queue only depends on the value of τ (or µ) through the ratio τ /µ, i.e., the invariance property holds.This is for instance the case for the shifted negative binomial distribution, with pgf with as special case the shifted geometric distribution (m = 1).
For the above combinations of service-demand and service-capacity distributions we know for sure that the invariance property holds.There may, however, be other combinations for which it holds, and in particular, there may be combinations for which the invariance property holds where neither the service demands nor the service capacities are geometrically distributed.The invariance property, however, does not hold in general.This will be illustrated in the next section, where numerical results (see Figure 3) indicate that the property does not hold anymore for the combination of shifted geometric demands and shifted geometric or deterministic service capacities.A full classification of all the conditions under which the invariance property holds thus remains an open question.

7.
Numerical examples and discussion.In this section, we present a few numerical examples that illustrate the behavior of the queueing system, and in particular the impact of the service-capacity distribution on several performance measures of the system.Throughout this section, we consider Poisson arrivals with mean arrival rate λ, i.e., A(z) = e λ(z−1) .
In the first example, shown in Figures 1 and 2, we study the impact of the service-capacity distribution on the mean and the variance of the customer delay under varying loads ρ = λτ /µ.Here, the load ρ is varied by varying λ, the service demands are constant, with τ = 11, and 4 different service-capacity distributions are considered, all with mean µ = 10: deterministic capacities with pgf geometric capacities with pgf and a weighted mixture of 2 geometric distributions with means 5 and 30 such that the overall mean µ = 10, with corresponding pgf .
From Figure 1 it can be seen that, generally, a higher variance of the service capacity leads to a higher mean delay.This is also what one would expect intuitively, since more variability on the service capacities is in general expected to lead to a burstier service process, which should in turn cause longer queues.However, we notice in Figure 1 also one case where the opposite is true: under low load, the system with negative binomially distributed service capacities turns out to have a lower mean delay than the system with deterministic capacities, despite the fact that the deterministic distribution has the lowest variance.Under high load it is the other way around again.
This observation can be explained as follows.When the load ρ is very close to 0, almost every customer arrives in an empty system.Due to the fact that the service demands are deterministically equal to 11 work units, any customer arriving in an empty system will be served in one slot if and only if the service capacity in the next slot is at least 11 work units.Since this is never the case for deterministic service capacities (with µ = 10), the minimum delay in this case is 2 slots.Indeed, in Figure 1 it can be seen that for deterministic service capacities, the mean delay goes to 2 as ρ goes to 0. For negative binomial service capacities however, the probability that the service capacity is at least 11 work units in a slot is approximately 40.4%, so that around 40.4% of the customers arriving in an empty system will be served in one slot, and if the load is low enough then the mean delay will be lower than 2 slots, as can be observed in Figure 1.On the contrary, when the load is high, relatively few customers arrive in an empty system, so the benefit that negative binomial service capacities give to these customers becomes insignificant.In this case, the system with deterministic service capacities will outperform the system with negative binomial service capacities, due to the lack of randomness in the service process, which allows customers to be served more regularly.
From the above discussion, we conclude that the queueing model studied in this paper cannot be reduced to a classical model with independent service times, because the "mean service time" depends on the arrival process and on the load ρ in particular.For instance, for deterministic service capacities, this mean service time approaches 2 as ρ approaches 0, whereas it approaches 11/10 as ρ approaches 1.
In Figure 2, we show the variance of the customer delay for the same system parameters as in Figure 1.We observe that the relative ordering of the variances of the customer delay for the various service-capacity distributions remains the same for all values of the load ρ, i.e., the lines in Figure 2 do not cross, unlike those in Figure 1.In Figure 2, a higher variance of the service capacity always corresponds to a higher variance of the customer delay.In particular, the variance of the customer delay is lower for deterministic service capacities than for negative binomial service capacities for all values of ρ, including those near 0. In fact, as ρ approaches 0, then for deterministic service capacities, the variance of the customer delay approaches 0 as well, which is not the case for negative binomial service capacities.This is because if the load is very low, then with deterministic service capacities almost all customers have a delay of 2 slots, whereas with negative binomial service capacities, around 40.4% of the customers have a delay of 1 slot, the rest has a delay of 2 or more slots, as discussed earlier.
In a second example, shown in Figure 3, we keep the load ρ and the mean arrival rate λ fixed at ρ = λ = 0.9, so we keep the ratio τ /µ = 1 and we scale τ and µ together to observe the impact of their actual values on the mean system content.The service demands have a shifted geometric distribution (with minimum value 1) and 5 different service-capacity distributions are considered, all with mean µ = τ : deterministic, negative binomial (with m = 3), geometric, shifted geometric (with minimum value 1), and binomial (with m = 10) service capacities.
We clearly see that for geometric service capacities (and shifted geometric demands), the mean system content does not depend on the actual values of τ and µ but merely on their ratio, as predicted by the invariance property discussed in Section 6.This means that higher service demands of the customers are in this case exactly compensated by proportionally equally higher service capacities of the system.However, note that while the invariance property holds for geometric service capacities, it does not hold for shifted geometric service capacities, as the mean system content in the latter case clearly depends on µ.The same is true for deterministic service capacities.Finally, from Figure 3 it can also be seen that the invariance property also holds for binomial and negative binomial service capacities, in agreement with the discussion in Section 6.In our final numerical example, shown in Figure 4, we consider Poisson arrivals with λ = 3 and service demands that are uniformly distributed between 1 and 10 work units.The service capacities follow a negative binomial distribution with µ = 10 and parameter m, and we examine the influence of the parameter m on the tail probabilities of the system content.We observe from Figure 4 that an increase in the parameter m reduces the probability that the system content becomes large.This could be expected intuitively, since an increase of m corresponds to a decrease of the variance of the service capacities, which makes the system serve the customers at a more regular rate, which in turn decreases the probability of a large system content.
We moreover see that the impact of increasing m is the largest when m is very small.For instance, the increase from m = 1 to m = 2 corresponds to a relatively large performance gain, whereas the step from m = 10 to m = 20 results in a comparatively small performance gain.This is because there is an upper limit to the amount of performance that can be gained by increasing m further and further.In Figure 4, we also show the tail probabilities for the case of deterministic service capacities, which can be considered as negative binomial service capacities with parameter m = ∞.The full line in Figure 4 therefore represents a lower bound on the tail probabilities that can be achieved by changing the parameter m.

Conclusion.
In this paper, we analyzed a discrete-time queueing model with general service demands and service capacities with rational pgf.Our main results are the obtained analytical expressions for the pgfs of the steady-state customer delay and the system content, the expressions for the mean values of the system content and delay, and the dominant-pole approximations for the tail probabilities of the system content and the delay.We also found new sufficient conditions under which the invariance property holds.
While the considered model is very general, allowing arbitrary distributions of the number of arrivals per slot and the service demands of the customers, and arbitrary phase-type distributions for the service capacity per slot, a restriction of the model is that the service capacities are assumed to be independent from slot to slot.Studying the effects of correlation between service capacities is a compelling direction for future research.of R(z) that of R(z) k are equal.By means of (50), we can rewrite equation (49) as The above infinite summation of contour integrals is equal to the contour integral of the infinite series (i.e., we may "swap" the summation and integration symbols) if the contour L is chosen such that the resulting infinite series is uniformly convergent.
It is important to question when such a contour can be constructed.
If |z| < R D , we may therefore construct L as described above and bring the summations in (51) behind the integral.We obtain Substituting z = 0 in (52), in view of D(0) = 0, we get Using this result in (52) again, we find In order to further simplify the above expression, we now split the integrand into two terms, as follows: The latter term has no poles inside L, since L was chosen such that ∀ζ ∈ L : |zR(ζ)| < 1, which implies (due to Rouché's theorem) that 1 − zR(ζ) has no zeros inside L, and since the simple zero of the denominator at ζ = 1 (if that would be inside L) is canceled by the zero of the numerator at ζ = 1.We conclude that the contribution of the latter term to the value of the contour integral in (53) is zero.Therefore we can rewrite (53) as Moreover, we change the integration variable in (55) to ξ = 1/ζ (which yields a factor −1/ξ 2 in the integrand), and we invert the integration path L into L but still integrate in counterclockwise sense (which yields an extra factor of -1, since the inversion of L is a clockwise path).This then leads to the desired relationship between the pgfs D(z) and V (z).
If the pgf R(z) of the service capacities is a rational function, the general relationship (47) can be further transformed into where the functions α k (z) are the m zeros for ξ of 1 − zR(1/ξ).The relationship (56) is valid for all z for which the zeros α k (z) are distinct, which is the case for all but at most 2m − 1 values of z.
Proof.If the pgf of the service capacities is a rational function, i.e., if R(1/z) is given by (1), equation (47) can be rewritten as We now focus on the poles of the integrand in (57) inside the contour L .Since the service demand of each customer is at least 1 work unit, S(0) must equal 0, and by ( 7) V (0) must equal 0 as well.The zero of the factor V (ξ) in the numerator of the integrand at ξ = 0 then ensures that the factor ξ in the denominator does not cause a pole of the integrand at ξ = 0. Furthermore, the factor Q R (ξ) − P R (ξ) in the numerator ensures that the factor (ξ − 1) in the denominator does not cause a pole of the integrand at ξ = 1.Finally, since the contour L was chosen such that ∀ξ ∈ L : |ξ| < R V , V (ξ) has no poles inside L either.Therefore, the only poles of the integrand in (57) inside L are the zeros for ξ of or equivalently, of 1 − zR(1/ξ).
(59) Since (58) is a polynomial in ξ of degree m, (58) has exactly m zeros for ξ.We denote these zeros for a given value of z by α k (z), k = 0, 1, ..., m−1.It can easily be seen that all these zeros lie inside L .Indeed, the contour L was chosen such that ∀ξ ∈ L : |zR(1/ξ)| < 1.This implies that |zP R (ξ)| < |Q R (ξ)|, so using Rouché's theorem we can say that (58) has as many zeros inside L as Q R (ξ).But all m zeros (counting with multiplicities) of Q R (ξ) must lie inside L , because L was chosen such that ∀ξ ∈ L : |ξ| > 1/R R .This means that (58) must have exactly m zeros for ξ inside L .
To calculate the value of the contour integral in (57), we can therefore apply Cauchy's residue theorem (see e.g.[16]).Note that the zeros α k (z) are not necessarily distinct.For a given value of z, let α(z) denote the set of zeros for ξ of (58).The pgf of D(z) is then obtained as where the residue at a pole ξ * with multiplicity m * is given by Since all quantities in expression (60) are known or can be calculated numerically (when z is known), this expression may be used to evaluate D(z) for any z.However, due to the (m * − 1)st derivative with respect to ξ in (61), the evaluation of D(z) may be difficult in practice if the zeros α k (z) are not distinct.

Theorem 5 . 1 .
If the zeros β k (z B ) are distinct, then the dominant pole z B of B(z) is the smallest positive real-valued zero from all the zeros of z − S(β k (z)), k = 0, 1, ..., m − 1.
From this contradiction we conclude that β k (z B ) = ζ for all k and ζ.(b) A pole of β k (z) for some k.However, this would not lead to a pole of B(z) because the factors β k (z) and β k (z) − ξ in the numerator have combined multiplicity m, while the factors β k (z) − ζ in the denominator also have combined multiplicity m.(c) A zero of R (1/β k (z)) for some k.However, this is not possible due to our assumption that the zeros β k (z B ) are distinct (see (20) and Appendix D).(d) A pole of S(β k (z))/(z − S(β k (z)) for some k.This can only occur if z is a zero of z − S(β k (z)).Since (d) is the only possibility, we conclude that z B is a zero of z − S(β k (z)) for some k ∈ [0, m − 1].Furthermore, since z B is the dominant pole, it must be the smallest positive real-valued z for which z = S(β k (z)) for some k ∈ [0, m − 1].

Figure 1 .
Figure 1.Mean customer delay versus the load ρ, for Poisson arrivals with varying λ, deterministic service demands with τ = 11 and various service-capacity distributions (as indicated), all with mean µ = 10.

Figure 2 .
Figure 2. Variance of the customer delay versus the load ρ, for Poisson arrivals with varying λ, deterministic service demands with τ = 11 and various service-capacity distributions (as indicated), all with mean µ = 10.

Figure 3 .
Figure 3. Mean system content versus the mean service demand τ for Poisson arrivals with λ = 0.9, shifted geometric service demands and various service-capacity distributions (as indicated), with mean µ = τ .

Figure 4 .
Figure 4. Dominant-pole approximation of the tail probabilities of the system content, for Poisson arrivals with λ = 3, uniformly distributed service demands from 1 to 10 work units, and negative binomial service capacities with µ = 10 and various values of the parameter m, as well as deterministic service capacities.
This is the case if ∀ζ ∈ L : |1/ζ| < R V and |zR(ζ)| < 1.The former condition imposes a lower bound on |ζ|, whereas the latter imposes an upper bound on |R(ζ)| that depends on z.Since this upper bound is most severe when ζ is real and positive, and since R(ζ) is an increasing function of ζ on the part of the real axis where 0 ≤ ζ < R R , the bounds can be rewritten as R(1/R V ) < R(|ζ|) < |1/z|.We conclude that a contour can be constructed if and only