An exact root-free method for the expected queue length for a class of discrete-time queueing systems

For a class of discrete-time queueing systems, we present a new exact method of computing both the expectation and the distribution of the queue length. This class of systemsincludesthebulk-servicequeueandtheﬁxed-cycletrafﬁc-light(FCTL)queue, whichisabasicmodelintrafﬁc-controlresearchandcanbeseenasanon-exhaustive time-limitedpollingsystem.Ourmethodavoidsﬁndingtherootsofthecharacteristic equation,whichenhancesboththereliabilityandthespeedofthecomputationscom-paredtotheclassicalroot-ﬁndingapproach.Werepresentthequeue-lengthexpectation inanexactclosed-formexpressionusingacontourintegral.Wealsointroduceseveral realisticmodiﬁcationsoftheFCTLmodel.FortheFCTLmodelforaturningﬂow, weproveadecompositionresult.Thisallowsustoderiveaboundonthedifference betweenthebulk-serviceandFCTLexpectedqueuelengths,whichturnsouttobe smallinmostoftherealisticcases.


Introduction
In the analysis of many queues, the probability-generating function (pgf) of the queue length is represented as a rational function with several unknown variables.The classical way of finding the unknowns is to use the analyticity of the pgf inside the unit disk; see [3,6].The analyticity implies that each zero of the denominator, i.e., each root of the corresponding characteristic equation, inside the closed unit disk is also a zero of the numerator, which yields a system of linear equations for the unknowns.In the case of a stable system and distinct zeros, it is possible to find all the required variables.The main drawback of this approach is finding the zeros, since, in general, no closed-form expression exists.Moreover, the precision of the found zeros also influences the result.
To avoid these problems, some authors suggest approximations of the queue length; see, for example, [11,16,20].Another approach is to use the matrix-analytic method; see [13].However, this also requires solving a (matrix) equation.
For a certain class of discrete-time queueing systems, we present an exact root-free approach to find the expected queue length and the queue-length distribution.The key point of our method is that the unknown probabilities depend on the zeros of the denominator in a symmetric way and that, as we show, the actual values of the zeros are not important.We apply the residue theorem and find these unknown probabilities in terms of contour integrals.The mean queue length is, in fact, a function of only one such contour integral.This makes our method computationally efficient in comparison with the standard root-finding method.Moreover, our method gives a closed-form solution without implicitly defined variables such as the roots of the characteristic equation in the root-finding method or the rate matrix in the matrix-analytic approach.Thus, our solution can be used, for example, to find asymptotic results for heavy load.
The considered class of queueing systems includes, among others, the bulk-service queue, see [3], and the fixed-cycle traffic-light (FCTL) queue, see [6].The server in the bulk-service queue serves a fixed-size batch of customers each time unit.If there are fewer customers in the queue, all of them are served.In each time unit, a random number of customers arrives, where the arrivals during different time units are independent of each other.The bulk-service queue also appears as an embedded Markov chain for the continuous-time multi-server M/D/s queue, see [10], or M X /G B /1 queue with batch arrivals and batch service.The FCTL model is the basic and key model in trafficcontrol research and can be seen as a polling system with non-exhaustive time-limited service discipline.The service, visiting and switching times of the server are fixed and deterministic.Unlike most polling systems, the server continues to serve even if the queue empties, as it has no knowledge of the queue length (there are no detectors), cf.[2,7].However, if the queue empties during the green (service) time, then, during the rest of the green time, all arriving customers are served instantly, i.e., arriving vehicles proceed without delay.We can also apply our method for queueing systems, where the pgf of the analyzed random variable depends on the minimal polynomial of the unknown zeros, i.e., the polynomial that has the same zeros (with the same multiplicity) and has the smallest power.An example of such a system is the discrete-time G/G/1 queue; see [18].
In this paper, we also consider a modification of the FCTL model for a lane with turning flow, where, unlike the FCTL model, the customers arriving during green time and finding an empty queue require the same service time as the queued vehicles.For this new model, we prove a decomposition result similar to [5].Moreover, we prove that given the same arrival process and server capacity, the expected bulk-service queue is bounded from below by the expected overflow queue for the FCTL model and from above by the expected overflow queue for the FCTL model for turning flow, where the overflow queue is the queue just after service.These results give us a bound on the difference between these models.We propose several other extensions of the FCTL queue, which include the effects of traffic interruption, for example, due to actuated control for pedestrians or uncertain departure times.
In the numerical results, we compare our method with the root-finding and matrixanalytic approaches.We observe that the root-finding approach can be numerically quite unreliable, while the matrix-analytic approach is much slower than our method.We also apply our method to the FCTL queue and its extensions.In particular, we investigate the impact of the variability of the arrival process on the expected queue length and consider the difference between straight-going and turning flows.In addition, we study effects of traffic disruption.
The paper is structured as follows: The method in its general form is presented in Sect. 2. In Sect.3, we apply our method to the bulk-service queue and compare it with the root-finding and matrix-analytic approaches in terms of speed and reliability.In Sect.4, we present several realistic extensions of the FCTL model.Finally, in Sect.5, we present the conclusions and propose future research directions.

Method
In this section, we propose a generic method to analyze various discrete-time queueing systems, such as the bulk-service queue and the FCTL queue; see Sects. 3 and 4. For these systems, the pgf of the queue length is a rational function with unknown probabilities in the numerator.The classical way of determining the unknowns requires finding the roots of the characteristic equation.In our method, we use the special structure of the numerator for these systems to find an alternative root-free solution.Let us consider the function X (z) of the following form: where x k are unknown coefficients that we want to determine; D(z), B(z) and f (z) are known functions such that B(1) = 1, f (1) = 0, and, for each zero z l ∈ {z : |z| 1, z = 1} of the denominator, f (z l ) = 0; and X (1) = 1.For the analysis in this section, we do not need X (z) to be a pgf.In what follows, we suppose that the following assumptions on the functions B(z), D(z), f (z) and X (z) hold: Assumption 1 (Analyticity assumption) For some α > 0, the functions B(z) and D(z) are analytic and the function f (z) is two times differentiable in the disk D 1+α = {z : |z| < 1 + α}.The function X (z) is continuous in the closed unit disk.
For any z = w, z, w ∈ D1 , the function B(z) satisfies Note that from Assumption 2, the following statement is deduced immediately.

Corollary 1
The equation z = B(z) has only one root in D1 , namely 1.
Remark 1 Without loss of generality, we can assume that D(0) = 0. Otherwise, both denominator and numerator can be divided by z n for some n > 0 such that the new D(z) satisfies the condition D(0) = 0.

Coefficients expressed in terms of roots
In this subsection, we discuss how x k , k = 0, . . ., g − 1, depend on the roots of (2). Let Since we know that D(0) = 0, then z k = 0 for each k = 0, . . ., g − 1. Rewrite X (z) in the following way: Since z k = 0 and f (z k ) = 0, k = 1, . . ., g − 1, the numerator is equal to 0 for z = z k and Consider the following polynomial: The function h(y) is a polynomial of degree g − 1, and y k , for k = 1, . . ., g − 1, are zeros of this polynomial.Note that whenever z k is a multiple root of the equation D(z) = 0, y k is a multiple root of the polynomial h(y) with the same multiplicity.Moreover, according to Assumption 2, y k = y l if z k = z l .Thus, h(y) can be written as By applying Vieta's formulas (see [19]) to the polynomial h(y), we get, for k = 1, . . ., g − 1, that where are elementary symmetric polynomials.Here, for notational convenience, we assume that σ 0 (y 1 , . . ., y g−1 ) = σ 0 = 1.
Eq. ( 4) gives us x k up to a normalization constant, which we derive from the equation X (1) = 1 by applying L'Hôpital's rule:

Coefficients expressed in terms of contour integrals
Using (4), we can find the coefficients x k if we know the values of the elementary symmetric polynomials σ k .We will represent these symmetric polynomials as functions of the following symmetric polynomials: for k = 1, . . ., g − 1. Below, in Theorem 1, we provide a way to find η k without computing the roots of (2).The proof of the theorem requires the following auxiliary result, which can be obtained using L'Hôpital's rule.
Lemma 1 Consider a function F(z) that is analytic in a neighborhood of z k , where k = 0, . . ., g − 1.Then, the residue of the function D (z) D(z) F(z) at z k is equal to where m z k is the multiplicity of the zero z k .
Theorem 1 Consider ε > 0 and δ > 0 such that there are no roots of Eq. (2) in D1+ε \ D1 and Dδ .Then, where S r = {z : |z| = r }.  1.Thus, by the residue theorem and Lemma 1, the difference of the integrals on the right-hand side of ( 7) is equal to

Remark 2
Strictly speaking, we also need ε < α; see Assumption 1.However, in many applications, α is very big or infinite, and, therefore, we omit this requirement hereafter.In Appendix B, we elaborate on how to choose ε and δ.
The following lemma provides a way of finding σ k , k = 1, . . ., g − 1, using η k , k = 1, . . ., g − 1.We omit the proof of this lemma since it requires only careful computation of the monomials' coefficients on the right side of Eq. (8).

Expectation expressed in terms of roots
Now, consider X (1).For queueing systems, this represents the expected queue length.

Expectation expressed in terms of contour integrals
We can use a contour integral to find Let ε be as above.Then where is the residue of the function D (z) By Corollary 1, we know that there are no roots of the equation z = B(z) inside D1 except 1.Thus, we do not need to subtract any other residues of the function inside the integral.
Remark 3 When X (z) is a pgf of the queue length, it is also possible to find the variance of the queue length, i.e., X (1) + X (1) − (X (1)) 2 , in the same way, but since it is a lengthy expression, we omit it here.The application of contour integrals goes further.Under certain conditions, the function X (z) can be written as an exponent of a contour integral, yielding a Pollaczek integral; see [4], where authors use our factorization result (3).
Note that the expectation X (1) represents the average queue length and is a real number.Thus, for B (1) ∈ R, only the real part of the integral (12) should be computed.

Algorithms
In this subsection, we summarize our contour-integral algorithms to compute X (1) and the coefficients x k using the results of Sects.2.1, 2.2, 2.3 and 2.4.For computational questions such as how to choose the parameters ε and δ or how to compute a complex integral, we refer to Appendix B.

Algorithm 1 (Computation of X (1))
1. Find ε > 0 such that there are only g roots of Eq. ( 2) in the disk D1+ε , using one of two ways described in Sect.B.3. 2. Compute (the real part of) the integral where Remark 4 Note that we parametrized S 1+ε by z(ϕ) = (1+ε)e iϕ .Therefore, dz = i zdϕ and where F(z) represents the integration function.
In Appendix B.1, we consider a simplified algorithm for the case when B(z) has no roots inside D1 , which is equivalent to B(−1) > 0 if B(z) is a pgf.Fortunately, in many applications, the function B(z) has this property.Therefore, it is often better to use the simplified algorithm as described in Appendix B.1.

Method application in queueing theory
In this subsection, we give some remarks concerning the application of our method.First, we study a special subclass of functions X (z) of the form (1) for which it is easy to check Assumption 2. This subclass generalizes the form of the queue-length pgf for the bulk-service and FCTL-type queues considered in Sects.3 and 4.Then, we describe the application to queueing systems that do not have a pgf of the form (1).
Suppose that D(z) = z g − A(z) (16) and the functions A(z), B(z) are pgfs.Consider the following assumption.
The following two lemmas show that Assumptions 1 and 3 imply Assumption 2. Therefore, if the denominator of X (z) has the form ( 16), and Assumptions 1 and 3 hold, we can apply our method.

Lemma 3 Consider a pgf A(z) that is analytic in the disk
For the proof of Lemma 3, we refer to [1].The proof of Lemma 4 is given in Appendix A.

Remark 5
In our examples, the second part, B (1) < 1, of the stability assumption will follow from the first part, A (1) < g.However, for the general case, we have to assume that B (1) < 1 to prove the second part of Assumption 2; see the proof of Lemma 4.
Note that for many queueing systems, it was proven that the unknown part of the pgf can be expressed in terms of the function where z k , k = 1, . . ., g − 1 are the only roots of a certain equation in D1 \{1}.For such systems, our approach can be used to find the pgf, see Eqs. ( 4), ( 7) and ( 8), or the expectation, see (10) and (12).One of the examples is a discrete-time G/G/1 queue; see [18].For the G/G/1 queue, the generating function of the limiting waiting-time distribution depends either on h(z) or on 1/(z g h(1/z)), depending on whether the pgf of the inter-arrival or service time is a rational function, i.e., is equal to where A 1 (z) and A 2 (z) are polynomials of (at most) power g.Similar results hold for the queue-length pgf in the case of a bulk-service queue with a random bulk size.The reason for this is that these two systems are governed by the same equation: For the bulk-service queue, X n is the queue length at time n, and A n and B n are the numbers of arrivals and departures at time n, respectively.For the G/G/1 queue, X n and A n are the waiting time and the service time of the nth customer, respectively, and B n is the inter-arrival time between the nth and (n + 1)st customers.

Method comparison
In this section, we consider a discrete-time bulk-service queue.For this example, we compare our contour-integral approach in terms of speed and reliability with the root-finding and matrix-analytic methods.In Sect.3.1, we describe the model.Then, we apply the considered methods in Sects.3.2, 3.3 and 3.4.Finally, we present the numerical results in Sect.3.5 and discuss the differences between the approaches in Sect.3.6.

Discrete bulk-service queue
Consider a discrete-time queue with bulk service.This means that in each time unit the server serves g customers in the queue, and during this time A b new customers arrive with a pgf A b (z).The server serves only those customers that are present in the queue at the beginning of the time unit.If there are fewer than g customers in the queue, all of them are served.Denote by X b (z) the pgf of the queue length at the beginning of a time unit in the stationary state.Let q k be the stationary probability of having k customers in the queue at the beginning of a time unit.Then (see, for example, [3]), we obtain The stability condition in this model is A b (1) < g.If we consider the pgf X − b (z) of the queue length after service but before arrivals, we get In the following subsections, we apply our method, the root-finding and matrixanalytic approaches to find the expectation (X − b ) (1).

Our method
At first glance, the pgf of the queue length, see (17), is not in the form ( 1), but we can rearrange it as follows: Hence, we take and q 1 = x 1 .In the case A b (0) = 0, we can reduce the complexity as before, since q l , as the probability of having l customers in the queue, is zero for l = 0, . . ., n − 1, where n is the least possible number of customers that arrive in one time slot.Thus, we suppose that A b (0) = 0. Hence, z k = 0 for every k = 1, . . ., g − 1 and the stability assumption holds.Therefore, we can use our method if A b (z) is analytic in the disk D 1+α for some α > 0. For the case of X − b (z), the only change would be that From ( 14), we get that the expected queue length is given by the following formula: where I is given by and ε is chosen as in Theorem 1.Note that ( 18) and ( 19) are equivalent to the following formula: Thus, for the case of X − b (z), we get

Root-finding approach
The root-finding method requires a numerical algorithm to find the roots of (2).Then, one needs to pick the roots that are inside the unit circle, namely the roots z 1 , . . ., z g−1 .
The coefficients x k in (1) then follow from a linear system of equations.For the case of the bulk-service queue, we get Note that the matrix is singular if two roots coincide.If the roots are close, the resulting values of x l depend heavily on the precision of the roots.Given values of x l , we compute the probabilities q k = x k − x k−1 , where x −1 = 0, and find the expected queue length as follows: Alternatively to solving (21), we can directly use the roots; see Eq. (11).For the bulk-service, we get The possible drawbacks of both approaches are that the root-finding algorithm can fail to find exactly g − 1 roots inside the unit circle and that the precision of the roots can be unsatisfactory.

Matrix-analytic approach
It is possible to represent the bulk-service queue as an M/G/1-type queue; see [13].
Consider the Markov chain of the queue length at service completion, i.e., before the arrivals.Then, the transition matrix is a block matrix where and a k is the probability of having k arrivals, with a k = 0 for k < 0. Using the notation for M/G/1-type queues, instead of the queue length l we consider a two-dimensional process (n, m), where the level state n is equal to the integer part of the quotient l/g, i.e., l/g , and the background state m is the remainder of such a division, i.e., l − l/g g.The steady-state distribution is found using the matrix G defined as the minimal nonnegative solution of the following equation: For numerical computation, this infinite summation should be truncated.The stationary probabilities π n = (π ng , . . ., π ng+g−1 ) for level n are computed recursively using Ramaswami's formula, see [8]: where and π 0 follows from the following system of linear equations: with B0 = B 0 + Ā2 G and e = (1, . . ., 1).The average queue length is approximated by where the value for N can be chosen so that the contribution of the vector π N is small, i.e., In [15], an alternative way was proposed to find the average queue length for M/G/1-type queues.For the bulk-service queue, we can use a simplified version of this approach.The average queue length is computed using π 0 and π * = ∞ i=1 π i , which can be found from the following system of linear equations: Note that the average queue length is the sum of the average level and the average background state, i.e., where b = (0, 1, . . ., g − 1).Thus, it is only left to find r = ∞ k=1 kπ k using In the following subsection, we numerically compare the considered approaches in terms of speed and reliability.

Numerical results
In this subsection, we compare our approach with the root-finding and matrix-analytic methods for the bulk-service queue.We focus on speed and reliability, i.e., the absence of failures such as a wrong number of roots, or unrealistic result (negative, imaginary or not accurate queue length); see Table 1.We generate 10,000 random parameters g = 2, . . ., 30, c = g+1, . . ., 70 and load ρ ∈ [0, 0.99].We set A b (z) = (ρgz/c+1− ρg/c) c , i.e., we consider binomial arrivals.We have used standard python methods for computing integrals, finding polynomial roots and solving linear systems of equations.
For the root-finding technique, the use of the linear system (21) leads to poor reliability.In almost 12% of the cases, the approach failed to give a realistic result, where in 3% of the cases, the number of roots in the unit circle was wrong, and there was no result at all.The root-finding algorithm together with (23) is two times faster and fails in only 4% of the cases.Note that the failures occur more often for large batch sizes; see Fig. 2. For a batch size of 30, the root-finding method with (21) fails in more than half of the cases.Our approach has no fails for all values of batch sizes and produces realistic results.
The matrix-analytic approach always produces a nonnegative queue length.However, the use of (24) leads to a serious underestimation of the queue length in 10% of the cases.This happens due to slow convergence of the distribution tail for high loads  There are 10,000 random combinations of g, c and λ.For the speed reference of our method, only the real part of the integral was calculated.The real and imaginary parts of the expected queue length are denoted by Re(EL) and Im(EL), respectively.The absolute difference is calculated between our and other approaches Fig. 2 The percentage of failed runs for root-finding algorithm depending on the batch size g (ρ > 0.9).The calculation of the aggregated distribution (given by π 0 , π * ) leads to results that (almost) coincide with the results of our approach, suggesting that both approaches produce accurate results.The main drawback is the slow computational speed, which is mainly due to the computation of the matrix G, which we find by recursively applying

Fig. 3
The average runtime as a function of g for the matrix-analytic approach using (24), T D , or (25), T A , and for our integral approach (20), T I with G 0 = 0, until absolute values of entries in the matrix G n+1 − G n are less than ε M = 0.1 10 .The runtime is longer for bigger sizes of the matrices, i.e., for bigger g; see Fig. 3.We conclude that our contour-integral method outperforms the considered approaches both in speed and reliability.In the following subsection, we summarize the advantages of our method compared to the existing approaches.

Discussion
The main contribution of the paper is our explicit method of finding the queue-length distribution.Both the root-finding approach and the matrix-analytic approach require the computation of implicitly defined variables, i.e., the roots z 1 , . . ., z g−1 or matrix G.For the root-finding approach, this step is error-prone, resulting in poor reliability of the method, while for the matrix-analytic approach, the main issue is the speed.Note that explicit closed-form results are much easier to manipulate to show a number of structural results such as convexity and sensitivity to input parameters or asymptotic behavior.
In Sect.3.5, we used a simple example of a bulk-service queue with binomial arrivals.For such arrivals, the pgf A b (z) is a polynomial, and A n = 0 for n c/g + 1.Thus, we use standard python methods to find the roots of the polynomials and do not need to truncate the arrival distribution to find the matrix G, which might result in an underestimation of the queue length.For our approach, the required pgf A b (z) is usually easy to compute even if it is not a polynomial; see Sect.4.5 below.Therefore, we do not have to choose a truncation bound as in the existing methods.
We note that, for certain arrival distributions, there are formulas to compute the roots; see, for example, [10].These formulas contain infinite summations, which should be truncated for numerical results.As we saw in Sect.3.5, a poor precision of the roots may lead to unrealistic results, which can also happen in the case of a badly chosen truncation bound.For the case of the bulk-service queue, root-free results were obtained in [9].The average queue length is represented as a (double) infinite summation, which also leads to a truncation problem.
Our results are applicable to a quite general class of queueing systems.This includes several M/G/1-type queues, such as bulk-service queues and FCTL-type queues described in Sect.4, but also discrete G/G/1 queues, see Sect.2.6, for which the matrix-analytic approach is not applicable.

FCTL-type queues
In this section, we consider several traffic-light models and apply our method to them.First, we explain the basic FCTL model in detail and use our approach to find the average queue length; see Sect.4.1.Then, we propose several modifications of this model to include the following realistic situations: 1. the lane serves a turning flow, so vehicles slow down before departing; 2. pedestrian and/or bicycle traffic lights have an actuated control; 3. just after the intersection there is a bridge or a railway and a part of the green time may be lost; 4. another lane at the intersection does not have a fixed length of the green time due to an actuated control; 5. times between departures differ due to driver distraction and/or vehicle condition and length.
In Sect.4.2, we consider the turning flow and prove a decomposition result.Even though situations 2-4 differ a lot, they can be modeled in the same way by assuming that the green and red times consist of a random number of time intervals; see Sect.4.3.In Sect.4.4, we suggest a way to model the distraction of the drivers.We conclude this section by presenting numerical examples in Sect.4.5.

FCTL model
Consider an isolated intersection with fixed-cycle traffic-light control, i.e., the green and red times of each lane are fixed.This system is modeled using a fixed-cycle trafficlight (FCTL) model, which is a basic model in traffic theory and has been extensively studied; see, for example, [6,17].Due to the fact that the green and red times are fixed, the queues at different lanes are independent and can be considered separately.Thus, we focus on one approaching lane and consider the queue length for this lane.The analysis for the other lanes is the same.Suppose that each delayed vehicle needs the same time τ to depart from the lane.Thus, if there is a long queue at the beginning of the green time, we see a departure every τ seconds.In what follows, we suppose that time is split in time intervals, each of τ seconds.The effective green time consists of g ∈ N time intervals and the effective red time of r ∈ N time intervals.Thus, not more than g delayed vehicles can depart during the green time.We suppose that each cycle starts with g green time intervals and then switches to r red time intervals, where by green and red time intervals we mean time intervals during green and red time, respectively.Together, this gives c = g + r time intervals in a cycle.
Let Y n,m be the number of arrivals during the nth time interval of the mth cycle with pgf Y n,m (z).We make the following two assumptions: Assumption 4 (Independent arrivals assumption) The numbers of arrivals Y n,m , for m = 1, 2, . .., n = 0, . . ., c − 1, are identically distributed and independent of each other.
Thus, we can denote Y n,m simply as Y and Y n,m (z) as Y (z).This assumption is realistic for an intersection that lies far enough from other signal-controlled intersections; for example, if the intersection is isolated.However, if the distance between intersections is short, then the vehicles arrive in platoons and this assumption does not hold.Following [17], we add the so-called FCTL assumption: Assumption 5 (FCTL assumption) If a vehicle arrives during the green time and finds an empty queue, then it proceeds without delay.
Assumption 5 means that if the queue is cleared before the end of the green time, the vehicles that arrive during the remaining green time will immediately depart.Therefore, the queue length will be 0 until the end of the green time.This assumption is quite realistic for a straight-going flow since vehicles that find no queue will proceed without stopping and, thus, are able to depart at the free-flow speed.However, for turning flows, this assumption is not realistic, especially for right turn (or left turn in the case of left-hand traffic), because the turning vehicles need to slow down.Thus, in the next subsection, we also consider the FCTL model with the one-vehicle assumption (see Assumption 6 below) instead of the FCTL assumption.
Denote by X n (z) the pgf of the queue length in the stationary state at the beginning of the nth time interval during an arbitrary cycle.The pgf exists in the case of a stable system, i.e., when the possible number of departures is greater than the expected number of arrivals: cY (1) < g.
Let q n be the stationary probability of having an empty queue at the beginning of the nth time interval during a cycle, i.e., q n = X n (0).Then, under Assumptions 4 and 5, we get The first equation in ( 26) is based on the fact that if there is at least one vehicle in the queue at the beginning of the green time interval, then one vehicle departs and Y vehicles arrive.In terms of pgfs, this means multiplying by Y (z)/z.However, if the queue is empty, it remains empty throughout the time interval.Thus, q n should not be multiplied by any additional function.During the red time, Y vehicles arrive each time interval and none depart.Therefore, we only multiply by Y (z) in the final two equations of (26).
This gives us, after some rearrangement, the pgf of the overflow queue, defined as the queue length at the beginning of the red time: The expression has the form (1), where ) and x k = q k .Note that by Corollary 1, f (z) = 0 inside and on the unit circle only for z = 1.Thus, f (z k ) = 0 for each root z k of (2), k = 1, . . ., g −1.We need to assume that Y (z) is an analytic function in the disk D 1+α for some α > 0.Then, Assumption 1 holds.Since we consider only stable systems, we assume that A (1) = cY (1) < g, which means that B (1) = Y (1) < 1, and, therefore, Assumption 3 holds.
One can prove, see, for example, [6], that the expected overflow queue EX g = X g (1) and the expected queue length at an arbitrary point during the cycle are connected in the following way: Now consider X g (1).From ( 14), we get the following theorem.
Theorem 2 For the FCTL model, the overflow queue is given by where and ε is chosen as described in Theorem 1.

FCTL model for the turning flow
For some cases, it may turn out that the FCTL assumption does not hold.For example, in the case of a right turn at an intersection (or a left turn in the case of left-hand traffic), vehicles that find the queue empty need to decelerate almost to the speed of a delayed vehicle.Thus, rather than the FCTL assumption, it is logical to make the following assumption: Assumption 6 (One-vehicle assumption) If a number of vehicles arrives during the green time interval and finds the queue empty, then only one of the vehicles proceeds without delay, and the others form a queue.
This assumption means that even if the queue was cleared during the previous green time intervals, the overflow queue may be nonempty.The stability condition in this model coincides with the stability condition of the FCTL model under the FCTL assumption.Equations (26) will change to Indeed, if there are no vehicles, then max{Y − 1, 0} vehicles will queue, with pgf (Y (z) − Y (0))/z + Y (0).Therefore, we obtain for n < g that . This gives us the above result after a small rearrangement.Hence, the pgf of the overflow queue is expressed as As in the case of the FCTL assumption, we can use our method for the FCTL model with the one-vehicle assumption if the system is stable, i.e., g < cY (1), and the function Y (z) is analytic in some disk D 1+α .Note that in the case of Bernoulli arrivals, Assumptions 5 and 6 coincide.In this case, Y (z) = (1−Y (0))z+Y (0), and, thus, z−Y (z) = Y (0)(z−1).Note also that even though the equation for the roots is the same, generally the coefficients q k may differ from those under the FCTL assumption.This happens because one of the roots is 1 and plugging it in the numerator gives 0 without resulting in an additional equation to find q k .The required gth equation comes from the fact that X g (1) = 1.This normalization equation has different forms for Assumptions 5 and 6.Under Assumption 5, we get that However, under Assumption 6, the normalization equation is of the following form: We see that these models coincide only if 1 − Y (0) = Y (1), i.e., this happens only in the case of Bernoulli arrivals.
From Eqs. (34) and ( 35) and the fact that the roots are the same in both models, it follows that where q k,fctl (resp., q k,1v ) is the stationary probability of having an empty queue at the beginning of the kth time interval during a cycle in the case of Assumption 5 (resp., 6).Moreover, by comparing Eqs. ( 27), ( 33) and ( 26), (32), we get that in the case of a stable system, for each k = 0, . . ., c − 1, where X k,fctl (z) and X k,1v (z) are pgfs at the beginning of the kth time interval for the FCTL model and the one-vehicle model, respectively.Note that X dif (z) is independent of g and c with X dif (1) = 1 and ) .Consider the case g = c = 1.In the case of the FCTL assumption, the queue is always empty once it becomes empty since there is no red time during which vehicles can form a queue.Therefore, in the case of stability, i.e., Y (1) < 1, we get X 0,fctl (z) = 1.Thus, we get X 0,1v = X dif (z), i.e., X dif (z) is the pgf of the queue with g = c = 1 and arrivals Y (z).Hence, we get the following decomposition result (compare with the decomposition results in polling systems, [5]): Theorem 3 For arbitrary g c, in the case of a stable system, the queue with the one-vehicle assumption can be considered as a sum of two independent queues: the FCTL queue with the same g, c and the one-vehicle queue with g = c = 1.
Note that X dif (z) can be also viewed as a pgf of the queue length of a bottleneck, for example, when part of the road has low speed limit and all the vehicles decelerate before this place.
We can as well change the model to allow more than one vehicle, but not all, to pass the junction if the queue is empty.Then, the changes will be only in the function f (z), and we can easily use our method if Assumptions 1 and 3 hold.For the same reason, there would be a similar decomposition result.The expected queue length in this new case will be less than the expected queue length in the case of the onevehicle assumption but more than the expected queue length in the case of the FCTL assumption.

Remark 6
The bulk-service model can be seen as an intermediate step between the considered traffic-light models.More precisely, for A b (z) = Y (z) c , it can be shown using a sample path argument that the average queue length satisfies the following inequalities: In traffic terms, the bulk-service model can be seen as a queue with service at the end of the green time that serves g (or fewer) vehicles and then waits one cycle for new arrivals.
From (36), the difference between the two FCTL models, described by Theorem 3, gives an upper bound on the difference between the FCTL and bulk-service models.In Sect.4.5.2,we show that most of the time this difference is small.

Random green and red times
Suppose we know that with probability θ r ,g the red time consists of r time intervals and the green time after this red time consists of g time intervals, r , g ∈ N ∪ {0}, r + g = 0. Thus, cycles start with a red time, followed by a green time.The lengths of the red time and the green time are independent of the lengths of red and green times before this cycle.We assume that the length of the green time is at most N , i.e., ∞ k=0 ∞ l=N +1 θ k,l = 0.For simplicity, we consider the smallest such N , so ∞ k=0 θ k,N > 0. As before, in each time interval, Y vehicles arrive with pgf Y (z).We also use the FCTL assumption.One can easily alter the formulas in the case of the one-vehicle assumption.
Let X (z) be the pgf of the overflow queue and X r ,g k (z) be the (conditional) pgf of the queue length at the beginning of the kth green time unit given that there are r red time units and g green time units in this cycle.Note that when the choices of r and g are known, the pgf of the queue length changes as in (26).Thus, where p r ,g k is the probability of having an empty queue at the beginning of the kth green time interval in the case of a cycle with r red and g green time intervals.Hence, with probability θ r ,g , the pgf of the queue length at the beginning of the next red time will be Therefore, the pgf of the overflow queue satisfies the following equation: Isolating the function X (z) gives us After multiplying both numerator and denominator by z N and making some rearrangements, we get where q k = r ,g θ r ,g p r ,g g−1−k and A(z) = r ,g θ r ,g (Y (z)) g+r z N −g .As we see, we can apply our method if A(z) and Y (z) are analytic in some disk D 1+α , A (1) < N and Y (1) < 1.
Note that if the summation 1 (z) = r ,g θ r ,g (g where G and R are random variables that denote the length of the green and red times during a cycle.Thus, if A (1) < N , we get E(G + R)Y (1) < EG.It immediately follows that Y (1) < 1.We see once again that the stability assumption holds in the case of a stable system.Now, we give three examples of possible applications of this model.

Interruption by pedestrians and/or cyclists
Let us focus on one lane and, for simplicity, on one source causing interruptions.Suppose that with probability p there is an arrival of a cyclist during a cycle.Suppose that we need t c green time intervals for switching on and off the cyclists' green time.
The green time for pedestrians and cyclists can be provided using the green time g of this lane or can be added as extra time to the cycle time c = g + r .In the first case, we have θ r ,g = 1 − p, θ r +t c ,g−t c = p, and in the second case, we have θ r ,g = 1 − p, θ r +t c ,g = p.

Interruption by trains and boats
This case is more-or-less the same as the previous one.The only difference is in the length of interruption.We assume that interruption happens during a cycle with probability p.We also assume that if it happens, then all green time is effectively red.Therefore, θ r ,g = 1 − p, θ c,0 = p.

Vehicle-actuated control on another lane
Suppose that another lane has a vehicle-actuated control, and, thus, the length of its green time is not deterministic.In this case, to find an approximation for the queue length on the considered lane, we assume that the lengths of red times during different cycles are independent of each other.Then, if we know the distribution of the green time for the other lane, we get a distribution of the red time for the considered lane.Let θ r ,g = p r , where p r is the probability of the red time length r , and g is fixed.In the case of a fixed cycle, we take θ r ,c−r = p r .

Remark 7
The analysis of the actuated-controlled lane is complicated due to the fact that the arrival and service processes are connected.We leave this problem for future research.

Uncertain departure time
In this subsection, we consider the following extension of the FCTL model.Suppose that, during a green time interval, the driver of the departing vehicle may be distracted with a fixed probability p ∈ [0, 1].Therefore, during the green time, even in the case of a nonempty queue, vehicles do not depart each time interval.Then, each driver has a geometrically distributed number of tries to depart from the queue.This can be modeled using the FCTL model in the following way: Suppose that during the green time there is an extra arrival with probability p.In this way, we compensate for the uncertainty in the departures.The distribution of the queue length in this model will be the same as the distribution in our new model.So, we consider the case that the number of the arrivals during the green time interval has pgf Y (z)( pz + (1 − p)).However, during the red time, we have still Y (z) arrivals.Hence, the pgf of the overflow queue is given by the following formula: Using the same arguments as in Sect.4.1, one can see that the pgf is of the general form described in Sect. 2. Thus, it is possible to use the same techniques in the case when Y (z) is analytic in some D 1+α , and the system is stable, i.e., cY (1) + gp < g.

Numerical results
In this subsection, we present the numerical results for FCTL-type models.First, we consider the impact of the variability of the arrival process for the FCTL queue.Then, we discuss the difference between the FCTL assumption and the one-vehicle assumption.Finally, we study the impact of the traffic interruptions caused by cyclists and uncertain departure times.Further numerical results can be found in the technical report [14].

Variability of the arrival process
In this subsection, we show the impact of the variability of the arrival process on the queue length for the FCTL model.In what follows, λ denotes the arrival rate per time interval, i.e., λ = Y (1).In Table 2, we summarize pgfs and variances of the considered arrival processes.
Note that Bernoulli arrivals have the least possible variance for a fixed λ.Indeed, the variance of an arbitrary arrival process with rate λ is equal to Y (1)+Y ( 1)−(Y (1)) 2 λ−λ 2 .In what follows, the parameter n of both the binomial and the negative binomial distribution is set to be equal to 2.
First, we focus on the expected delay per vehicle for various types of arrivals.The expected delay is computed using Little's law ED = EL/λ, where EL = EL fctl is computed by (31).In Fig. 4, the expected delay is plotted as a function of the load ρ = cλ/g for c = 60 and g = 5 or 40.The expected delay is given in seconds, rather than time intervals, and one time interval is set to be equal to 2 s.
Our first observation is that the higher the variability of the arrival process, the longer the expected delay.However, for each green time g, the difference is significant only for a load higher than 0.8.As we can see, for high load, the relative difference between expected delays for various types of arrivals is increasing as the green time increases.However, the absolute difference is decreasing.For load ρ = 0.9833 (the highest load in the figures) and g = 5, 15, 30, 40, the absolute differences are given in Table 3.In Fig. 5 The probabilities q k of a queue being empty, k = 0, . . ., g − 1, for g = 10, c = 20, n = 2 and a λ = 0.2, b λ = 0.45 the table, we use ED bern , ED binom , ED pois , ED negbin to express the expected delay in the cases of Bernoulli, binomial, Poisson and negative binomial arrivals, respectively.
From the table, we see that the absolute difference is decreasing as the green time increases.In many approximations (see [12,16]), the dependence on the variance is set to be linear for the fixed arrival rate, green times and cycle times.However, we see that this difference increases for higher variance.
Next, we consider the probabilities q k of a queue being empty for k = 0, . . ., g − 1.In Fig. 5, for g = 10, c = 20 and λ = 0.2, 0.45, these probabilities are plotted as functions of k.For fixed g, c and λ, the sum of the probabilities is the same for various types of arrival, but the distributions of this sum are different.As we see, for all arrival rates, at the beginning of the green time the probability of a queue being empty is lower for the arrival processes that have lower variability, whereas at the end of green time the situation is reversed.The graphs look very different depending on the arrival rates.For low arrival rates, the graphs are concave and the difference between graphs for different arrival types is small.For high arrival rates, the graphs are convex and the relative difference between arrival types is larger.

Remark 8
The speed and reliability of our method, which we showed in Sect.3.5, allows us to consider the green time allocation problem for an isolated traffic-light intersection.This is the problem of choosing the green time length (value of g) for different lanes under constraints on the total green time; see [14].
Fig. 6 The expected extra delay due to the one-vehicle assumption as compared to the FCTL assumption

Comparing the FCTL and one-vehicle assumptions
Let us first consider X dif (1)/λ, i.e., the expected difference in the delay between FCTL and one-vehicle models.We plot this as a function of the arrival rate λ = Y (1) in Fig. 6.
Note that if at an intersection there are at least two conflicting "main" phases, then both of them use less than half of the cycle time.Thus, in the case of a stable system, the arrival rate is lower than 0.5.If we consider, for example, Poisson arrivals, the expected extra delay is 1 time interval, i.e., not more than 2 s.Therefore, for most of the cases the extra delay will be very short.Similarly, in these cases, the extra queue length is less than one vehicle.
As we have shown, the absolute difference between FCTL and one-vehicle models is quite small.However, it is interesting to consider the relative difference.For the same settings as in Fig. 4, we plot the relative difference in the expected delay in Fig. 7.
For a shorter green time, the delay observed in the FCTL model is already very long and the expected difference X dif (1)/λ is small.Thus, the relative difference is very small.For a longer green time, the delay is shorter and the arrival rate is higher, so the relative difference is higher.As we see, for g = 40, this relative difference reaches 10% for negative binomial arrivals with n = 2.

Disruption of the traffic
In this subsection, we consider disruption of traffic due to cyclists (or pedestrians).Suppose that cyclists need five time intervals, i.e., 10 s, to cross the road.There are two ways to provide the required green time.We can either shorten the green time of one lane (thus making room for cyclists to cross), or, alternatively, extend the total cycle time and, hence, add extra red time to all lanes.Let p be the probability of cyclists arrival during the cycle.In Fig. 8, we plot the overflow queue as a function of the rate for different p and g and fixed c = 60.In each graph we plot only up to load ρ = 0.975.The arrival process is assumed to be Poisson.As we see, the first way to deal with cyclists is highly disadvantageous for the lane.It significantly decreases the capacity of the lane and increases the overflow queue.The effect is stronger for a shorter green time.The conclusion is that it is better, if possible, to increase the cycle time.

Uncertain departure times
Consider the FCTL model with uncertain departure times that was presented in Sect.4.4.In Fig. 9, we plot the expected overflow queue as a function of the arrival rate for different probabilities of departure.We suppose that, on average, a vehicle departs in 2 s.The length of the time interval is set to be equal to τ = 2 p seconds, where p is the probability of departure.The arrival process is assumed to be Poisson.We fixed the cycle time (in seconds) and consider two different green times.Depending on p = 1, 5/6, 5/7, there are on average five departures in 5, 6 or 7 time intervals.The uncertainty in departure times does not influence the capacity of the system but only increases the overflow queue and, consequently, the delay.The overflow queue is a bit longer for a green time of 10 s (about 2 vehicles for the highest considered load 0.9833).The effect of the uncertainty becomes noticeable for loads over 0.8.Note that the average overflow queue is almost 0 for loads up to 0.7-0.8(depending on the green time), which means that most of the vehicles are served during one cycle time after the arrival.

Conclusions
We have presented a new method to calculate the expectation and distribution of the queue length for several discrete-time queueing systems.The solution is given in terms of contour integrals without any implicitly defined variables.We successfully applied this method to the bulk-service queue and several traffic-light models, providing insights in traffic behavior.For the FCTL model with the one-vehicle assumption, we proved a decomposition result, which allowed us to give a bound on the difference between the FCTL and bulk-service queue-length expectations, which turns out to be small for many relevant applications.
We have shown that our method outperforms the classic root-finding and matrixanalytic methods in terms of speed and reliability.Using our method, the queue-length expectation can be found more than 7 and 17 times faster compared to the root-finding and matrix-analytic methods, respectively.Moreover, we found that the root-finding approach can fail to give a realistic result.This happens more often for the case of many roots, where up to 10% of the considered cases result in no answer and 40% in the negative or nonreal (complex) expected queue length.The speed of the matrixanalytic approach depends on the size of the block matrices, which coincides with the number of roots, and can be up to 54 times slower than our approach.Future work will be on extending our integral method to a wider range of systems.This can be rewritten as Here, we used that B (1) < B(1) + b 0 .Thus, we proved for each z ∈ D1 inequality (39), which is equivalent to (38).
We proved that d|z z| dt > 0 for each z = z(t) ∈ D 1+δ \{0}.Since b 0 + B(0) = 0, we get that z(t) = 0. Therefore, any root of Eq. (37) for t > 0 goes outside D1+δ .Hence, the number of roots inside D 1 decreases when t increases.But, for small t > 0, we can use Rouché's theorem to show that there is only one root of Eq. (37).Thus, we have proved Lemma 4. (41) Hence, we do not need to compute two integrals as in (7).As before, we need to subtract the residue at 1, which is equal to 1.This leads to the following simplified algorithm to find unknown variables x k : Algorithm 3 (Computation of x k , k = 0, . . ., g − 1, in the case of zero-free B(z)) 1. Find ε > 0 such that there are only g roots of Eq. ( 2) in the disk D1+ε , using one of two ways described in Sect.B. The necessary and sufficient condition for a pgf B(z) to be zero-free in D 1 is given by the following lemma.1− < 1.Therefore, according to Lemma 4, C(z) has not more than one zero in the unit disk.However, C(t) < 0 < C(−1), C (1).Thus, by the intermediate value theorem, it has at least two zeroes.We get a contradiction.

Fig. 4
Fig. 4 The expected delay as a function of load for Bernoulli, binomial, Poisson and negative binomial arrivals with c = 60, n = 2 and a g = 5, b g = 40

Fig. 7 Fig. 8 Fig. 9
Fig. 7 The relative difference (in %) in the expected delay between the FCTL and one-vehicle models as a function of load for binomial, Poisson and negative binomial arrivals for c = 60, n = 2 and a g = 5, b g = 40

Lemma 7 1 −
The pgf B(z) is zero-free in D 1 if and only if B(−1) > 0.Proof Suppose that B(−1) > 0 and there is a zero t of B(z) in D 1 .Since B(z) has at most one zero and B(t) = B(t), this zero is on the real line.Consider some 0 < < min(B(−1), 1 − B (1)) and the function C(z) = B(z)− B (1) > , the function C(z) is a pgf and C (1) = B(1)

Table 1
Comparison of speed and reliability for the expected queue length

Table 2
The

Table 3
The differences (in seconds) in the expected delay for different types of arrival process, c = 60 and load ρ = 0.9833