1 Introduction

In this era of rapidly emerging new technologies for urban traffic control, the vast majority of urban traffic intersections are still priority-based unsignalized intersections, in which the major roads have priority over the minor roads. Due to the increasing level of traffic congestion, the need for state-of-the-art quantitative analysis methods is greater than ever. Performance analysis of unsignalized traffic intersections is traditionally based on so-called gap acceptance models, with roots that can be traced back to classical queueing models. The basis of these models is the assumption that the crossing decision of a driver on the minor road is based on the gap between two successive vehicles on the major road (Drew 1968). This model can be applied in other contexts as well, e.g. when analyzing freeways (Drew et al. 1967a, b; Munjal and Pipes 1971) and pedestrian crossings (Tanner 1951; Wei et al. 2018).

In this paper, we contribute to the modeling and analysis of unsignalized intersections. In classical gap acceptance models vehicles on the minor road accept any gap greater than the critical gap (or critical headway), and reject gaps below this threshold, where the gap is defined as the time between two subsequent vehicles on the major road. However, in reality, drivers typically do not need the full critical gap; the remaining part can be used by the next vehicle on the minor road. In the literature, the part of the critical gap that is really used by the vehicle is referred to as the merging time.

Although gap acceptance models have improved greatly since their introduction more than fifty years ago, there are still some fundamental limitations to their practical applicability. As pointed out recently by Liu et al. (2014), two of its main limitations are the inability to incorporate (1) different driver behaviors and (2) heterogeneous traffic; see also Prasetijo (2007) and Prasetijo and Ahmad (2012). The main reason why these features could not be included is that the currently used analysis techniques are based on a single-server queueing model with exceptional first service (see also Welch 1964; Yeo 1962; Yeo and Weesakul 1964). In the gap acceptance literature, this model is commonly referred to as the M/G2/1 queue, a term seemingly introduced by Daganzo (1977). Importantly, in this queue, the ‘service times’ (corresponding to the time required to search for a sufficiently large gap and crossing the intersection or, depending on the application, merging with the high-priority traffic flow) are assumed to be independent. Incorporating more realistic features, such as driver impatience or different types of driver behavior, creates dependencies that make the model significantly more difficult to analyze.

An important development, that helps us overcome this obstacle, is the recently developed framework for the analysis of queueing models that allows a highly general correlation structure between successive service times (the M\(^X\)/SM2/1 queue discussed in Abhishek et al. 2018). Although set up with general applications in mind, this framework turns out to be particularly useful in road traffic models, where the dependencies can be used to model clustering of vehicles on the major road, to differentiate between multiple types of driver behavior, or to account for heterogeneous traffic. The present paper focuses on the latter two features.

In the existing literature, various quantitative methods have been used to study unsignalized intersections. The two most common procedures are empirical regression techniques and gap-acceptance models (cf. Brilon et al. 1997), but other methods have been employed as well (for example the additive conflict flows technique by Brilon and Wu 2002). A topic of particular interest concerns the capacity of the minor road (Brilon and Miltner 2005), which is defined as the maximum possible number of vehicles per time unit that can pass through an intersection from the minor road. Other relevant performance measures are the queue length and the delay on the minor road.

One of the first studies using queueing models for unsignalized intersections was by Tanner (1962), who assumed a constant critical gap and move-up time. Tanner first determines the mean delay of the low priority vehicles and characterizes the capacity of the minor road as the arrival rate of the low-priority vehicles, at which the mean delay grows beyond any bound. Initially one worked with a constant critical gap (for all users), resulting in a nice, compact, closed-form expression for the capacity (cf. Welch 1964; Yeo 1962; Yeo and Weesakul 1964). Then one empirically observed that the statistical variations of the drivers’ critical gaps have an impact; Siegloch (1973) and Harders (1968, 1976) developed some theoretically based variations and corrections. One of the adaptations was applied to the capacity estimate in which the old number is multiplied with \(f:=1-10^{-7}\cdot q^2\), where q is the arrival rate on the major road, to account for various sources of randomness. Later studies (see Brilon 1988) mention that the approach with a universal factor f does not always work. In addition, it is empirically observed that drivers tend to finally accept gaps of a value they have rejected before (Abou-Henaidy et al. 1994; Findeisen 1971; Retzko 1961; Tonke 1983; Tupper 2011). This resulted in models that also included driver impatience (cf. Abhishek et al. 2016; Drew et al. 1967a, b; Weiss and Maradudin 1962), which complicates the analysis drastically (cf. Heidemann and Wegmann 1997), in particular when trying to maintain a realistic merging behavior.

Tanner’s model has been generalized in various ways (Catchpole and Plank 1986; Cheng and Allam 1992; Hawkes 1965; Heidemann 1994; Wegmann 1991; Weiss and Maradudin 1962) by allowing random critical gaps and move-up times, also analyzing performance measures such as the queue length and the waiting time on the minor road. Heidemann and Wegmann (1997) give an excellent overview of the earlier existing literature. Moreover, they add a stochastic dependence between the critical gap and the merging (move-up) time, and study the minor road as an M/G2/1 queue, to determine the queue length, the delay and the capacity on the minor road. They distinguish between three types of driver behavior. The first is the basic model in which all low-priority drivers are assumed to use the same fixed critical gap. In the second model, critical gaps are random, with drivers sampling a new critical gap at each new attempt. This is typically referred to as inconsistent driver behavior. The third model is also known as the consistent model, in which a random critical gap is sampled by each driver for his first attempt only, and the driver then uses that same value at his subsequent attempts. The present paper generalizes all three gap acceptance models into a more realistic model that covers various realistic driving behavior features.

Although the literature on gap acceptance models is relatively mature and the existing gap acceptance models have proven their value, they are sometimes criticized for looking like rigorous mathematics, but in reality being based on pragmatic simplifications. As a consequence, the produced results might be of a correct magnitude, but are, however, only of approximative nature (Brilon and Wu 2002). The main goal of this paper is to increase the significance of gap acceptance models by taking away some of these concerns by including essential new features that are typically encountered in practice. In more detail, the main novelties of this paper are the following.

  • First, we incorporate driver impatience behavior while allowing for a realistic merging behavior. Put more precisely, we distinguish between the critical gap and the merging time, thus allowing multiple vehicles to use a sufficiently large gap. Incorporating this feature is particularly challenging in models with driver impatience.

  • Secondly, we allow for multiple classes of gap acceptance behavior, enabling us to distinguish between different driver types and/or different vehicle types.

  • Thirdly, we use a queueing model in which vehicles arriving in an empty queue on the minor road have different service-time distributions than the queueing vehicles (where ‘service time’ is meant in the sense introduced above). This setup facilitates the analysis of the queue length on the minor road as well as the service-time distribution of an arbitrary low-priority vehicle. The capacity for the minor road vehicles is then derived from the expectation of the service time for queuers.

The remainder of this paper is organized as follows. In Sect. 2 we introduce our model. Then the full queue-length analysis is presented in Sect. 3, whereas the capacity of the minor road is determined in Sect. 4. In Sect. 5 numerical examples demonstrate the impact of driver behavior on the capacity of the minor road.

2 Model description

In this section, we provide a detailed mathematical description of our new model. We study an unsignalized, priority controlled intersection where drivers on the major road have priority over the drivers on the minor road. For notational convenience, we will focus on the situation with one traffic stream on each road, but several extensions also fall in our framework (see, for example, Fig. 1). Vehicles on the major road arrive according to a Poisson process with intensity q. On the minor road, we have a batch Poisson arrival process with \(\lambda \) denoting the arrival intensity of the batches (platoons) and B denoting the (random) platoon size. We assume that traffic on the major road is not hindered by traffic on the minor road, which is a reasonable assumption in most countries (but not all, see Liu et al. 2014).

Fig. 1
figure 1

The situation analyzed in this paper

Full size image

The drivers on the minor road cross the intersection as soon as they come across a sufficiently large gap between two subsequent vehicles on the major road. Any gap that is too small for a driver on the minor road is considered a failed attempt to cross the intersection and adds to the impatience of the driver.

In the existing literature, three model variants have been introduced: (1) the standard model with constant critical gaps, (2) inconsistent gap acceptance behavior where a new random critical gap is sampled at each attempt, and (3) consistent behavior where drivers sample a random critical gap which they use for all attempts (cf. Abhishek et al. 2016; Heidemann and Wegmann 1997). In this paper, we apply a framework that allows us to create a single model that generalizes all three aforementioned model variations, with a combination of consistent and inconsistent behavior, allowing for heterogeneous traffic and driver impatience. This is a major enhancement of the existing models, as such a general model has not been successfully analyzed so far. Unfortunately, this requires a significant adaptation of the underlying M/G2/1 queueing model that has been the basis of all gap acceptance models so far.

We will now describe the full model dynamics in greater detail. We distinguish between \(R\in \{1,2,\ldots \}\) driver/vehicle profiles. The profile of any arriving vehicle on the minor road is modeled as a random variable that equals r with probability \(p_r\) (with \(\sum _{r=1}^R p_r=1\)). The profile represents the consistent part of the gap selecting behavior. Each driver profile has its own critical gap acceptance behavior, allowing us to distinguish between drivers that are willing to accept small gaps and (more cautious) drivers that require longer gaps, but also to allow for heterogeneous traffic with, for example, trucks requiring larger gaps than cars or motorcycles.

We introduce \(T_{(i,r)}\) to denote the critical gap of a driver with profile r during his ith attempt. The impatient behavior is incorporated in our framework by letting the critical gap \(T_{(i,r)}\) depend on the attempt number \(i=1,2,3,\dots \). This generalizes classical gap acceptance models where the critical gap is the same for all driver types and throughout all attempts.

Finally, we introduce the inconsistent component of our model by allowing the critical gap of a driver of type r in attempt i (\(T_{(i,r)}\), that is) to be a discrete random variable which can take any of the following values:

$$\begin{aligned} T_{(i,r)}= {\left\{ \begin{array}{ll} u_{(i,1,r)} &{}\quad \text { with probability }p_{(i,1,r)},\\ u_{(i,2,r)} &{}\quad \text { with probability }p_{(i,2,r)},\\ \vdots &{}\\ u_{(i,M,r)} &{}\quad \text { with probability }p_{(i,M,r)}, \end{array}\right. } \end{aligned}$$

with \(\sum _{k=1}^M p_{(i,k,r)}=1\) for all (i, r). Inconsistent behavior is used in the existing literature when the driver samples a new random critical gap after each failed attempt. It is generally being criticized for not being very realistic because it is unlikely that a driver’s gap acceptance behavior fluctuates significantly throughout multiple attempts. In our model, however, it is an excellent way to include some randomness in the gap selection process (which is due to the fact that not every driver from the same profile will have exactly the same critical gap at his ith attempt), because we can limit the variability by cleverly choosing the values \(u_{(i,k,r)}\).

Realistic values for the critical gaps \(u_{(i,k,r)}\) can be empirically obtained by measuring the lengths of all rejected and accepted gaps for each vehicle type. One needs to distinguish between systematic fluctuations due to different driver/vehicle types and (smaller) random fluctuations within the same driver profile. Since the index i represents impatience, it makes sense that \(u_{(i,k,r)}\) is decreasing in i. Moreover, since r denotes the driver profile, \(u_{(i,k,r)}\) is relatively large for driver profiles r that represent slow vehicles/drivers and vice versa. The variability in \(u_{(i,1,r)}, \dots , u_{(i,M,r)}\) can be limited to avoid unrealistic fluctuations. For example: if driver profile \(r=1\) corresponds to trucks and \(r=2\) corresponds to motorcycles, then it may be quite realistic that

$$\begin{aligned} \min _{1\le k\le M} u_{(i,k,1)} > \max _{1\le k\le M} u_{(i,k,2)}. \end{aligned}$$

Paraphrasing, the slowest motorcycle will always be faster than the fastest truck.

To make the model even more realistic, we assume that the actual vehicle merging time, denoted by \(\Delta _r\) for profile r, is less than the critical gap \(T_{(i,r)}\). As a consequence, the remaining part of the critical headway can be used by following drivers.

At first sight, our model has a few seeming limitations. In the first place, we assume the critical gap has a discrete distribution with M possible values. In the second place, as will be discussed in great detail in the next section, the analysis technique requires that there is a maximum number of possible attempts: \(i\in \{1,2,\dots ,N\}\). It is important to note, however, that both issues do not have big practical consequences since one can choose M and N quite large. This does come at the expense of increased computation times, obviously.

3 Queue length analysis

The main objective of this section is to determine the stationary queue length on the minor road at departure epochs of low-priority vehicles. In the next section, we will use these results to determine the service-time distribution of the low-priority vehicles, which enables us to evaluate the capacity of the minor road.

3.1 Preliminaries

We start by describing the queueing process on the minor road. All arriving vehicles experience three stages before leaving the system: queueing, scanning and merging. The queueing phase is defined as the time between joining the queue and reaching the front of the queue (for queuers this is the moment that the preceding driver accepts a gap and departs). The scanning phase commences when the driver starts scanning for gaps and ends when a sufficiently large gap is found. The scanning phase is followed by the merging phase, which ends at the moment that the vehicle merges into the major stream. Note that phases 1 and 2 may have zero length. The length of phase 1 is called the delay. The length of phase 3 for a vehicle of profile r is denoted by the profile-specific constant \(\Delta _r\). We define the service time as the total time spent in phases 2 and 3 (scanning and merging).

Remark 1

In some papers, an additional move-up phase is distinguished but, depending on the application (freeway or unsignalized intersection), this phase can be incorporated in one of the other phases. See Heidemann and Wegmann (1997, Section 3) for more details.

The queue-length analysis is based on observing the system at departure epochs. In contrast to the classical M/G2/1 queueing models, in our current model we now have dependent service times. This dependence is caused by the fact that the length of the lag (i.e., the remaining gap) left by the previous driver, due to the impatience, now depends on the attempt number and the profile of the previous driver. Using the notation introduced in the previous section: we need to know the critical gap \(T_{(i, r)}\) of the previous driver of profile r that led to a successful merge/crossing of the intersection. A queueing model with this type of dependencies is called a queueing model with semi-Markovian service times, sometimes referred to as the M/SM/1 queue. Although several papers have studied such a queueing model (cf. Abhishek et al. 2017; Adan and Kulkarni 2003; inlar 1967; de Smit 1986; Gaver 1963; Neuts 1966, 1977a, b; Purdue 1975), we still need to make several adjustments in order to make it applicable to our situation. The analysis below is based on the framework that was recently introduced in Abhishek et al. (2017).

Denote by \(X_n\) the queue length at the departure epoch of the nth vehicle, i.e. the number of vehicles behind the nth driver at the moment that he merges into the major stream. We use \(G^{(n)}\) to denote the service time of the nth vehicle and \(A_n\) to denote the number of arrivals (on the minor road) during this service time. Since vehicles arrive in batches (platoons), we denote by \(B_n\) the size of the batch in which the nth vehicle arrived. Let \(X(\cdot )\), \(A(\cdot )\), and \(B(\cdot )\) denote, respectively, the probability generating functions (PGFs) of the limiting distributions of these random variables,

$$\begin{aligned} X(z)=\lim _{n\rightarrow \infty }{\mathbb {E}}[z^{X_n}], \quad A(z)=\lim _{n\rightarrow \infty }{\mathbb {E}}[z^{A_n}], \quad B(z)=\lim _{n\rightarrow \infty }{\mathbb {E}}[z^{B_n}], \qquad |z| \le 1, \end{aligned}$$

whereas \({\tilde{G}}(\cdot )\) is the limiting Laplace-Stieltjes transform (LST) of \(G^{(n)}\):

$$\begin{aligned} {\tilde{G}}(s) =\lim _{n\rightarrow \infty }{\mathbb {E}}[e^{-s G^{(n)}}], \qquad s\ge 0. \end{aligned}$$

Note that the arrivals constitute a batch Poisson process. Therefore, A(z) can be expressed in terms of B(z) and \({\tilde{G}}(s)\),

$$\begin{aligned} A(z)={\tilde{G}}(\lambda (1-B(z))), \end{aligned}$$
(3.1)

but finding an expression for \({\tilde{G}}(s)\) is quite tedious. For this reason, we split the analysis into two parts. This section discusses how to find X(z) in terms of A(z), while Sect. 4 shows how to find G(z) for the gap acceptance model considered in this paper.

The starting-point of our analysis is the following recurrence relation:

$$\begin{aligned} X_n = \left\{ \begin{array}{l l} X_{n-1}-1+A_n &{} \quad \text {if }X_{n-1} \ge 1 \\ A_n +B_n-1&{} \quad \text {if }X_{n-1} =0 \end{array} \right. , ~~~ n=1,2,3,\dots . \end{aligned}$$
(3.2)

As argued before, \(\{X_n\}_{n\ge 1}\) is not a Markov chain. In order to obtain a Markovian model, we need to keep track of all the characteristics of the \((n-1)\)th driver. Let \(J_n= (i,k, r)\) contain these characteristics of the \((n-1)\)th driver, where i is the ‘succeeded attempt number’, k indicates which random gap the driver sampled, and r denotes the profile of the driver.

In Abhishek et al. (2018) it is shown that the ergodicity condition of this system is given by

$$\begin{aligned} \lim _{n\rightarrow \infty } {\mathbb {E}}[A_n|X_{n-1}\ge 1]<1. \end{aligned}$$
(3.3)

Throughout this paper, we assume that this condition holds.

3.2 A queueing model with semi-Markovian service times and exceptional first service

To find the probability generating function (PGF) of the queue-length distribution at departure epochs, we use the framework introduced in Abhishek et al. (2018). There an analysis is presented of the M\(^X\)/SM2/1 queue, which is a very general type of single-server queueing model with batch arrivals and semi-Markovian service times with exceptional service when the queue is empty, extending earlier results (cf. also Abhishek et al. 2017; inlar 1967; Gaver 1963; Neuts 1966) to make it applicable to situations like the one discussed in the present paper. To improve the readability of this paper, we briefly summarize the most important results from Abhishek et al. (2018) that are also valid for our model.

In Abhishek et al. (2018) the type of the nth customer is denoted by \(J_{n}\in \{1,2, \dots , N\}\). In order to apply these results correctly to our model, we need to make one small adjustment. In order to determine the service time of the nth customer, we need to know exactly how much of the critical gap of the \((n-1)\)th customer remains for this nth customer. As it turns out in the next section, this requires the following knowledge about the \((n-1)\)th customer:

  1. (A)

    we need to know the sampled critical gap \(u_{(i,k,r)}\) of the \((n-1)\)th customer;

  2. (B)

    we need to know whether the \((n-1)\)th customer emptied the queue at the minor road and the nth customer was the first in a new batch arriving some time after the departure of his predecessor.

To start with the latter: a consequence of (B) is that we need a queueing model with so-called exceptional first service. As a consequence of (A), we use the triplet (i, k, r), extensively discussed in the previous subsection, to denote the ‘customer type’, which should be interpreted as a vehicle of profile r that accepted the ith gap, where it sampled \(u_{(i,k,r)}\) as its critical headway. It means that if the nth customer is of type (i, k, r), then his predecessor (which is the \((n-1)\)th customer) succeeded in his ith attempt, while being of profile r and having critical gap \(u_{(i,k,r)}\).

A translation from our model to the model in Abhishek et al. (2018) can easily be made by mapping our \({\bar{N}}:=NMR\) customer types onto their N customer types (for instance, if \(j=(i-1)MR+(k-1)R+r\), then \(j\in \{1,2,\dots ,{\bar{N}}\}\) can serve as the customer type in Abhishek et al. 2018).

Mimicking the steps in Abhishek et al. (2018), it is immediate that

$$\begin{aligned}&{\mathbb {E}}[z^{X_{n}}1_{\{J_{n+1}=(j,l,r_1)\}}] \nonumber \\&\quad =\,\frac{1}{z}\sum _{r_0=1}^R\sum _{k=1}^{M}\sum _{i=1}^{N} {\mathbb {E}}[z^{X_{n-1}}1_{\{J_{n}=(i,k,r_0)\}}]{\mathbb {E}}[z^{A_{n}}1_{\{J_{n+1}=(j,l,r_1)\}}|J_{n}=(i,k,r_0),X_{n-1} \ge 1]\nonumber \\&\qquad +\frac{1}{z}\sum _{r_0=1}^R\sum _{k=1}^{M}\sum _{i=1}^{N}\Big (B(z){\mathbb {E}}[z^{A_{n}}1_{\{J_{n+1}=(j,l,r_1)\}}|J_{n}=(i,k,r_0),X_{n-1}=0]\nonumber \\&\qquad -{\mathbb {E}}[z^{A_{n}}1_{\{J_{n+1}=(j,l,r_1)\}}|J_{n}=(i,k,r_0),X_{n-1}\ge 1]\Big ){\mathbb {P}}(X_{n-1}=0,J_{n}=(i,k,r_0)),\nonumber \\ \end{aligned}$$
(3.4)

for \(n\in {{\mathbb {N}}}\), \( j\in \{1,2,\dots ,N\}\), \(l\in \{1,2,\dots ,M\}\), and \(r_1\in \{1,2,\dots ,R\}\).

To further evaluate these expressions, we need to introduce some additional notation. Define, for \(i,j\in \{1,2,\dots ,N\}\)\(k,l\in \{1,2,\dots ,M\}\), \(r_0,r_1\in \{1,2,\dots ,R\}\), and \(|z| \le 1\),

$$\begin{aligned} A^{(j,l,r_1)}_{(i,k,r_0)}(z)&=\lim _{n\rightarrow \infty }{\mathbb {E}}[z^{A_{n}}1_{\{J_{n+1}=(j,l,r_1)\}}|J_n=(i,k,r_0),X_{n-1}\ge 1],\\ A^{*(j,l,r_1)}_{\,\,(i,k,r_0)}(z)&=\lim _{n\rightarrow \infty }{\mathbb {E}}[z^{A_{n}}1_{\{J_{n+1}=(j,l,r_1)\}}|J_n=(i,k,r_0),X_{n-1}=0], \end{aligned}$$

where \(A^*\) corresponds to the number of arrivals during the (exceptional) service of a vehicle arriving when there is no queue in front of him. In addition,

$$\begin{aligned} f_{(j,l,r_1)}(z)&= \lim _{n \rightarrow \infty }{\mathbb {E}}[z^{X_n}1_{\{J_{n+1}=(j,l,r_1)\}}], \end{aligned}$$
(3.5)

such that

$$\begin{aligned} f_{(j,l,r_1)}(0)&= \lim _{n \rightarrow \infty } {\mathbb {P}}(X_n=0,J_{n+1}=(j,l,r_1)). \end{aligned}$$
(3.6)

The above definitions entail that

$$\begin{aligned} X(z)=\sum _{r_1=1}^R\sum _{l=1}^{M}\sum _{j=1}^{N}f_{(j,l,r_1)}(z). \end{aligned}$$
(3.7)

In steady state, Eq. (3.4) thus leads to the following system of \({\bar{N}}\) linear equations: for any j,  l, and \(r_1\),

$$\begin{aligned}&\left( z-A^{(j,l,r_1)}_{(j,l,r_1)}(z)\right) f_{(j,l,r_1)}(z)-\sum _{\begin{array}{c} r_0=1,\\ r_0\ne r_1 \end{array}}^RA^{(j,l,r_1)}_{(j,l,r_0)}(z)f_{(j,l,r_0)}(z) \nonumber \\&\qquad -\sum _{r_0=1}^R\sum _{\begin{array}{c} k=1,\\ k\ne l \end{array}}^{M}A^{(j,l,r_1)}_{(j,k,r_0)}(z)f_{(j,k,r_0)}(z) -\sum _{r_0=1}^R\sum _{k=1}^{M}\sum _{\begin{array}{c} i=1,\\ i\ne j \end{array}}^{N}A^{(j,l,r_1)}_{(i,k,r_0)}(z)f_{(i,k,r_0)}(z)\nonumber \\&\quad =\sum _{r_0=1}^R\sum _{k=1}^{M}\sum _{i=1}^{N}(B(z)A^{*(j,l,r_1)}_{\,\,(i,k,r_0)}(z)-A^{(j,l,r_1)}_{(i,k,r_0)}(z)) f_{(i,k,r_0)}(0). \end{aligned}$$
(3.8)

From the above we conclude that when we know \(A^{(j,l,r_1)}_{(i,k,r_0)}(z)\) and \(A^{*(j,l,r_1)}_{\,\,(i,k,r_0)}(z)\), the PGFs of our interest can be evaluated. In the next section, we show how to find the Laplace-Stieltjes transforms of the conditional service-time distributions, which lead to \(A^{(j,l,r_1)}_{(i,k,r_0)}(z)\) and \(A^{*(j,l,r_1)}_{\,\,(i,k,r_0)}(z)\).

We refer to Abhishek et al. (2018, Section 3) for more details on how to write the linear system of equations (3.8) in a convenient matrix form and solve it numerically. In Abhishek et al. (2018), and in more detail in Abhishek et al. (2017), it is also discussed how the PGF of the queue-length at arbitrary epochs (denoted by \(X^{\mathrm{arb}}\)) can be found; in particular the clean relation between \(X^{\mathrm{arb}}\) and X:

$$\begin{aligned} {\mathbb {E}}[z^{X^{\mathrm{arb}}}]={\mathbb {E}}[z^{X}] \frac{{\mathbb {E}}[B](1-z)}{1-B(z)} \end{aligned}$$
(3.9)

is useful in this context. The interpretation is that the number of customers at an arbitrary epoch is equal to the number of customers left behind by an arbitrary departing customer, excluding those that arrived in the same batch as this departing customer. We refer to Abhishek et al. (2017) for more details and a proof of (3.9).

4 Capacity

In this section, we derive the Laplace–Stieltjes transform of the service-time distribution of our gap acceptance model. This service-time analysis serves two purposes. Firstly, it is used to complete the queue-length analysis of the previous section. Secondly, it is an important result by itself, because it is directly linked to the capacity of the intersection, which is defined as the maximum possible number of vehicles per time unit that can pass from the minor road. An alternative but equivalent definition of capacity, also employed by Heidemann and Wegmann (1997), is the maximum arrival rate for which the corresponding queue remains stable. Combining (3.1) and (3.3), and denoting

$$\begin{aligned} g:=\lim _{n\rightarrow \infty }{\mathbb {E}}[G^{(n)}\,|\,X_{n-1}\ge 1], \end{aligned}$$

we can rewrite the stability condition as

$$\begin{aligned} \lambda {\mathbb {E}}[B]\,g<1, \end{aligned}$$
(4.1)

where \(\lambda \) is the arrival rate of batches, implying that the arrival rate of individual vehicles is \(\lambda {\mathbb {E}}[B]\). As a consequence, the capacity of the minor road is given by

$$\begin{aligned} C=\frac{1}{g}. \end{aligned}$$
(4.2)

Every driver (of profile r) samples a random \(T_{(i,r)}\) at his ith attempt, irrespective of his previous attempts. However, he uses only a constant \(\Delta _r\) of the gap to cross the main road. As a consequence, when accepting the gap, the remaining part \((T_{(i,r)}-\Delta _r)\) can be used by subsequent drivers. In order to keep the analysis tractable, we need the following assumption.

Assumption 1

(Limited gap reusability assumption) We assume that at most one driver can reuse the remaining part of a critical gap accepted by his predecessor. The means that \((T_{(i,r)}-\Delta _r)\) should not be large enough for more than one succeeding vehicle to use it for his own critical headway.

To avoid confusion, we stress that this assumption does not mean that a large gap between successive vehicles on the major road cannot be used by more than one vehicle from the minor road. It does mean, however, that the difference between critical gaps of two successive vehicles cannot be so large that the second vehicle’s critical gap is completely contained in its predecessor’s critical gap (excluding the merging time).

Remark 2

This assumption seems to be realistic in most, but not all, cases, because one can envision for instance situations in which (parts of) a critical gap accepted by, say, a slow truck might be reused by two fast accelerating vehicles following this truck. It is noted, however, that if Assumption 1 is violated, despite our method no longer being exact, the computed capacity is still very close to the real, simulated capacity; numerical evidence of this property is provided in Sect. 5. It implies that our method can still be used as a very accurate approximation in those (rare) cases where Assumption 1 is violated.

We recall that there is a complicated dependence structure between service times of two successive vehicles, due to our assumption that a vehicle can use part of the gap left behind by its predecessor combined with the assumption that drivers have different profiles and tend to become impatient. This obstacle can be overcome by analyzing the service time of a vehicle of type \((j, l, r_1)\) in the situation that its predecessor was of type \((i, k, r_0)\).

We proceed by introducing some notation. Let \(E_y\) be the event that there is a predecessor gap available of length y. Denote by \(\pi _{(i,k,r_0)}\) the probability that an arbitrary driver with profile \(r_0\) succeeds in his ith attempt with critical headway \(u_{(i,k,r_0)}\). In addition, let

$$\begin{aligned} {\tilde{G}}_{(j,l,r_1)}(s,y)&={\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r_1)\}}\,|\,E_y\right] , \end{aligned}$$
(4.3)

for \(0\le y \le u_{(i,k,r_0)}-\Delta _{r_0}\), \(j=1,2,\dots ,N\), and \(l=1,2,\dots ,M\). The LST of an arbitrary service time can be found by conditioning on the type of the current vehicle and its predecessor:

$$\begin{aligned} {\tilde{G}}(s)={\mathbb {E}}[e^{-sG^{(n)}}]=\sum _{r_0=1}^R\sum _{r_1=1}^R\sum _{j=1}^N\sum _{l=1}^M\sum _{i=1}^N\sum _{k=1}^M{\tilde{G}}^{(j,l,r_1)}_{(i,k,r_0)}(s){\mathbb {P}}(J_n=(i,k,r_0)), \end{aligned}$$
(4.4)

where

$$\begin{aligned} {\tilde{G}}^{(j,l,r_1)}_{(i,k,r_0)}(s)={\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r_1)\}}|J_n=(i,k,r_0)\right] , \end{aligned}$$
(4.5)

and

$$\begin{aligned} {\mathbb {P}}(J_n=(i,k,r_0))=\pi _{(i,k,r_0)}p_{(i,k,r_0)}p_{r_0}. \end{aligned}$$
(4.6)

We are now ready to state the main result of this paper. Define

$$\begin{aligned} {\bar{f}}_{(i,k,r)}:= \frac{f_{(i,k,r)}(0)}{\pi _{(i,k,r)}p_{(i,k,r)}p_{r}},\,\,\, {\bar{u}}_{(i,k,r)}:= u_{(i,k,r)}-\Delta _{r}. \end{aligned}$$

Theorem 1

The partial service-time LST of a customer, jointly with the event that he is of type \((j,l,r_1)\), given that his predecessor was of type \((i,k,r_0)\), is given by

$$\begin{aligned}&{\tilde{G}}^{(j,l,r_1)}_{(i,k,r_0)}(s)\,=\,{\tilde{G}}_{(j,l,r_1)}(s,{\bar{u}}_{(i,k,r_0)})\left( 1-{\bar{f}}_{(i,k,r_0)}\right) \nonumber \\&\quad +\Bigg (\int _{x=0}^{{\bar{u}}_{(i,k,r_0)}}\lambda e^{-\lambda x}{\tilde{G}}_{(j,l,r_1)}(s,{\bar{u}}_{(i,k,r_0)} -x)\mathrm{d}x + {\tilde{G}}_{(j,l,r_1)}(s,0) e^{-\lambda {\bar{u}}_{(i,k,r_0)}}\Bigg ){\bar{f}}_{(i,k,r_0)}, \end{aligned}$$
(4.7)

where

\(\circ \) :

\(f_{(i,k,r_0)}(0)\) is obtained by solving the system of equations (3.8), and

\(\circ \) :

expressions for \({\tilde{G}}_{(j,l,r_1)}(s,y)\) and \(\pi _{(i,k,r_0)}\) are provided in Lemmas 1 and 2 below.

Proof

This partial LST \({\tilde{G}}^{(j,l,r_1)}_{(i,k,r_0)}(s)\) can be expressed in terms of \({\tilde{G}}_{(j,l,r_1)}(s,y)\) by first distinguishing between the case where there was no queue upon arrival of the nth customer (\(X_{n-1}=0\)), and the case where there was a queue (\(X_{n-1}\ge 1\)). If the customer arrives in an empty system at time t and the previous arrival took place (merged on the major road) at time \(t-x\), then we need to check whether the remaining gap \(u_{(i,k,r_0)}-\Delta _{r_0}-x={\bar{u}}_{(i,k,r_0)}-x\) is greater than zero, or not. If the system was not empty at arrival time, the lag is simply \({\bar{u}}_{(i,k,r_0)}\) due to Assumption 1 (we provide more details on this in Remark 3). This, combined with (3.6) and (4.6), yields

$$\begin{aligned}&{\tilde{G}}^{(j,l,r_1)}_{(i,k,r_0)}(s) \\&\quad = \lim _{n\rightarrow \infty }{\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r_1)\}}|X_{n-1}\ge 1,J_n=(i,k,r_0)\right] {\mathbb {P}}(X_{n-1}\ge 1|J_{n}=(i,k,r_0)) \\&\qquad +\lim _{n\rightarrow \infty }{\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r_1)\}}|X_{n-1}=0,J_n=(i,k,r_0)\right] {\mathbb {P}}(X_{n-1}=0|J_{n}=(i,k,r_0))\\&\quad = {\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r_1)\}}| E_{(u_{(i,k,r_0)}-\Delta _{r_0})}\right] \left( 1-{\bar{f}}_{(i,k,r_0)}\right) \\&\qquad + \Bigg (\int _{x=0}^{\infty }\lambda e^{-\lambda x}{\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r_1)\}}| E_{\max (u_{(i,k,r_0)}-\Delta _{r_0}-x,0)}\right] \mathrm{d}x \Bigg ){\bar{f}}_{(i,k,r_0)}, \end{aligned}$$

which can easily be rewritten to (4.7). \(\square \)

It thus remains to find expressions for \({\tilde{G}}_{(j,l,r_1)}(s,y)\) and \(\pi _{(i,k,r_0)}\), which we do in the following two lemmas. We denote by \(E^{(n)}_{r}\) the event that the nth driver does not succeed in his first attempt, given that he has profile r.

Lemma 1

The conditional service-time LST of a vehicle, given the event \(E_y\), is given by

$$\begin{aligned} {\tilde{G}}_{(1,l,r)}(s,y)&= p_{r}p_{(1,l,r)}e^{-q\max \{u_{(1,l,r)}-y,0\}}e^{-s\Delta _{r}},\\ {\tilde{G}}_{(j,l,r)}(s,y)&= p_{r}p_{(j,l,r)}\left( 1-{\mathbb {E}}\left[ e^{-(s+q)\max \{T_{(1,r)}-y,0\}}\right] \right) \nonumber \\&\quad \times e^{-(s(y+\Delta _{r})+qu_{(j,l,r)})} \left( \frac{q}{s+q}\right) ^{j-1}\prod \limits _{m=2}^{j-1}\left( 1-{\mathbb {E}}\left[ e^{-(s+q)T_{(m,r)}}\right] \right) , \end{aligned}$$

for \(j=2,3,\dots ,N\).

Proof

For \(j=1\), the expression \({\tilde{G}}_{(j,l,r)}(s,y)\) simply follows from computing the probability that a customer is of type rand merges successfully during its first attempt with critical gap \(u_{(1,l,r)}\), for \(l=1,2,\dots ,M\), and faces a gap that is greater than \(u_{(1,l,r)}-y\). In this case, the customer’s service time is \(\Delta _{r}\). For the case \(j\ge 2\) we find

$$\begin{aligned}&{\tilde{G}}_{(j,l,r)}(s,y) \\&\quad =p_{r}\sum _{l_1=1}^{M}p_{(1,l_1,r)}1_{\{u_{(1,l_1,r)}> y\}}\int _0^{u_{(1,l_1,r)}-y}qe^{-qt} e^{-s(y+t)} {\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r)\}}|E^{(n)}_{r}\right] \mathrm{d}t \\&\quad =p_{r}\sum _{l_1=1}^{M}p_{(1,l_1,r)}1_{\{u_{(1,l_1,r)}> y\}} \frac{qe^{-sy}}{s+q}(1-e^{-(s+q)(u_{(1,l_1,r)}-y)}){\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r)\}}|E^{(n)}_{r}\right] \\&\quad =p_{r}\sum _{l_1=1}^{M}p_{(1,l_1,r)} \frac{qe^{-sy}}{s+q}(1-e^{-(s+q)\text {max}\{u_{(1,l_1,r)}-y,0\}}){\mathbb {E}}\left[ e^{-sG^{(n)}}1_{\{J_{n+1}=(j,l,r)\}}|E^{(n)}_{r}\right] \\&\quad =p_{r}\frac{qe^{-sy}}{s+q}\left( 1-{\mathbb {E}}\left[ e^{-(s+q)\text {max}\{T_{(1,r)}-y,0\}}\right] \right) {\mathbb {E}}\left[ e^{-sG^{(n)}} 1_{\{J_{n+1}=(j,l,r)\}}|E^{(n)}_{r}\right] ; \end{aligned}$$

here, it should be kept in mind that

$$\begin{aligned} {\mathbb {E}}\left[ e^{-sG^{(n)}} 1_{\{J_{n+1}=(j,l,r)\}}|E^{(n)}_{r}\right] \end{aligned}$$

does not depend on n and can be rewritten as

$$\begin{aligned}&\left( \prod \limits _{m=2}^{j-1}\sum _{k=1}^{M}p_{(m,k,r)} \int _0^{u_{(m,k,r)}}qe^{-qt}e^{-st}\mathrm{d}t \right) p_{(j,l,r)} \int _{u_{(j,l,r)}}^\infty qe^{-qt}e^{-s\Delta _{r}}\mathrm{d}t \\&\quad =p_{(j,l,r)}e^{-(s\Delta _{r}+qu_{(j,l,r)})} \left( \frac{q}{s+q}\right) ^{j-2}\prod \limits _{m=2}^{j-1}\sum _{k=1}^{M}p_{(m,k,r)}(1-e^{-(s+q)u_{(m,k,r)}}) \\&\quad =p_{(j,l,r)}e^{-(s\Delta _{r}+qu_{(j,l,r)})} \left( \frac{q}{s+q}\right) ^{j-2}\prod \limits _{m=2}^{j-1}\left( 1-\sum _{k=1}^{M}p_{(m,k,r)}e^{-(s+q)u_{(m,k,r)}}\right) \\&\quad =p_{(j,l,r)}e^{-(s\Delta _{r}+qu_{(j,l,r)})} \left( \frac{q}{s+q}\right) ^{j-2}\prod \limits _{m=2}^{j-1}\left( 1-{\mathbb {E}}\left[ e^{-(s+q)T_{(m,r)}}\right] \right) . \end{aligned}$$

Slightly rewriting this last expression completes the proof of Lemma 1. \(\square \)

The next lemma shows how to compute the probabilities \(\pi _{(i,k,r)}\). We denote by \(\tau _q\) the generic interarrival time between two vehicles on the major road, which is exponentially distributed with parameter q. Let \({{\hat{v}}}\) be defined as \(\max \{v,0\}\), and \({\check{v}}\) as \(\max \{-v,0\}\). In addition, \(v^{(1,k,r)}_{(i,l,r_0)}=u_{(1,k,r)}-u_{(i,l,r_0)}+\Delta _{r_0}\).

Lemma 2

The probability that a driver of profile r is served in his ith attempt with critical headway \(u_{(i,k,r)}\) is given by

$$\begin{aligned} \pi _{(i,k,r)}=P_{(i,k,r)}\prod \limits _{m=1}^{i-1}\left( 1-\sum _{k_m=1}^Mp_{(m,k_m,r)}P_{(m,k_m,r)}\right) , \end{aligned}$$
(4.8)

where, for \(i\in \{2,3,\dots ,N\}\) and \(k\in \{1,2,\dots ,M\}\),

$$\begin{aligned} P_{(i,k,r)}={\mathbb {P}}(\tau _q\ge u_{(i,k,r)})=e^{-qu_{(i,k,r)}}. \end{aligned}$$
(4.9)

The remaining \(P_{(1,k,r)}\) can be found by solving the system of linear equations

$$\begin{aligned}&\Bigg (1 +p_{(1,k,r)}\left( -p_{r}e^{-q{{\hat{v}}}^{(1,k,r)}_{(1,k,r)}}+c_{(k,r,r)} \right) \Bigg )P_{(1,k,r)} \nonumber \\&\qquad +\sum _{\begin{array}{c} r_0=1,\\ r_0\ne r \end{array}}^R\Bigg (-p_{r_0}e^{-q{\hat{v}}^{(1,k,r)}_{(1,k,r_0)}}+c_{(k,r_0,r)} \Bigg )p_{(1,k,r_0)}P_{(1,k,r_0)}\nonumber \\&\qquad +\sum _{r_0=1}^{R}\sum _{\begin{array}{c} l_1=1,\\ l_1\ne k \end{array}}^M\Bigg (-p_{r_0}e^{-q {\hat{v}}^{(1,k,r)}_{(1,l_1,r_0)}}+c_{(k,r_0,r)} \Bigg )p_{(1,l_1,r_0)}P_{(1,l_1,r_0)} \nonumber \\&\quad =\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N}\Bigg (1-e^{-\lambda {\check{v}}^{(1,k,r)}_{(i,l,r_0)}} -e^{-q{\hat{v}}^{(1,k,r)}_{(i,l,r_0)}}+\frac{\lambda }{\lambda +q}e^{-(\lambda +q){\check{v}}^{(1,k,r)}_{(i,l,r_0)}-qv^{(1,k,r)}_{(i,l,r_0)}} \nonumber \\&\qquad +\frac{q}{\lambda +q}e^{-\lambda {\bar{u}}_{(i,l,r_0)}-qu_{(1,k,r)}} \Bigg )f_{(i,l,r_0)}(0) +\sum _{r_0=1}^Rc_{(k,r_0,r)}, \end{aligned}$$
(4.10)

for \(k\in \{1,2,\dots ,M\}\) and \( r\in \{1,2,\dots ,R\}\), where

$$\begin{aligned} c_{(k,r_0,r)}=\sum _{l=1}^M\sum _{i=2}^{N}p_{r_0}p_{(i,l,r_0)}e^{-q(u_{(i,l,r_0)}+{\hat{v}}^{(1,k,r)}_{(i,l,r_0)})} \prod \limits _{m=2}^{i-1}\left( 1-{\mathbb {E}}[e^{-qT_{(m,r_0)}}]\right) . \end{aligned}$$

Proof

Expression (4.8) simply follows from the definition of \(\pi _{(i,k,r)}\), by multiplying the probabilities that the driver does not succeed in attempts \(1, 2, \dots , i-1\) (each time distinguishing between all M possible random critical gap values) and finally multiplying with the probability that he succeeds in the ith attempt with critical headway \(u_{(i,k,r)}\).

For \(i\in \{2,3,\ldots ,N\}\) it is easy to determine \(P_{(i,k,r)}\) because we do not have to take into account any gap left by the predecessor; we simply compute the probability that the critical gap \(u_{(i,k,r)}\) is smaller than the remaining interarrival time \(\tau _q\) (which, due to the memoryless property, is again exponentially distributed with rate q).

The case \(i=1\) is considerably more complicated. Using (4.6) and noting that \(\pi _{(1,k,r)}=P_{(1,k,r)}\), we can find a system of equations for \(P_{(1,k,r)}\) by conditioning on the type of the predecessor of the vehicle under consideration. Define

$$\begin{aligned} \chi _{(i,l,r_0)}^{(j,k,r)}:= {\mathbb {P}}(J_n=(j,k,r)|J_{n-1}=(i,l,r_0)), \end{aligned}$$

and

$$\begin{aligned} \phi _{(i,l,r_0)}^{(j,k,r)}&:= {\mathbb {P}}(J_n=(j,k,r)|J_{n-1}=(i,l,r_0),X_{n-2}=0), \end{aligned}$$
(4.11)
$$\begin{aligned} \psi _{(i,l,r_0)}^{(j,k,r)}&:= {\mathbb {P}}(J_n=(j,k,r)|J_{n-1}=(i,l,r_0),X_{n-2}\ge 1). \end{aligned}$$
(4.12)

Evidently,

$$\begin{aligned} P_{(1,k,r)}&=\frac{1}{p_{(1,k,r)}p_r}{\mathbb {P}}(J_n=(1,k,r)) \nonumber \\&=\frac{1}{p_{(1,k,r)}p_r}\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N} \chi _{(i,l,r_0)}^{(1,k,r)}\,{\mathbb {P}}(J_{n-1}=(i,l,r_0)) \nonumber \\&=\frac{1}{p_{(1,k,r)}p_r}\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N}\chi _{(i,l,r_0)}^{(1,k,r)}\,\pi _{(i,l,r_0)}p_{(i,l,r_0)}p_{r_0} \nonumber \\&=\frac{1}{p_{(1,k,r)}p_r}\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N}\Big (\phi _{(i,l,r_0)}^{(1,k,r)}\,{\mathbb {P}}(X_{n-2}=0|J_{n-1}=(i,l,r_0))\nonumber \\&\quad + \psi _{(i,l,r_0)}^{(1,k,r)}\,{\mathbb {P}}(X_{n-2}\ge 1|J_{n-1}=(i,l,r_0))\Big )\pi _{(i,l,r_0)}p_{(i,l,r_0)}p_{r_0}, \end{aligned}$$
(4.13)

which can be further evaluated, using notation introduced before, as

$$\begin{aligned}&\frac{1}{p_{(1,k,r)}p_r}\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N} \Bigg (\phi _{(i,l,r_0)}^{(1,k,r)}\, {{\bar{f}}_{(i,l,r_0)} } +\psi _{(i,l,r_0)}^{(1,k,r)}\left( 1-{\bar{f}}_{(i,l,r_0)}\right) \Bigg ){\bar{\pi }}_{(i,l,r_0)} , \end{aligned}$$

with \({\bar{\pi }}_{(i,l,r_0)}:= \pi _{(i,l,r_0)}p_{(i,l,r_0)}p_{r_0}\). We now consider the individual terms separately. Observe that, by conditioning on the value of \(\tau _q\),

$$\begin{aligned}&\frac{1}{p_{(1,k,r)}p_r}\phi _{(i,l,r_0)}^{(1,k,r)}\, {{\bar{f}}_{(i,l,r_0)} }{\bar{\pi }}_{(i,l,r_0)}\\&\quad =\Bigg (\int _{0}^{\max \{{\bar{u}}_{(i,l,r_0)}-u_{(1,k,r)},0\}}\lambda e^{-\lambda x}\mathrm{d}x \\&\qquad +\int _{\max \{{\bar{u}}_{(i,l,r_0)}-u_{(1,k,r)},0\}}^{{\bar{u}}_{(i,l,r_0)}}{\mathbb {P}}(\tau _q \ge u_{(1,k,r)}-{\bar{u}}_{(i,l,r_0)}+x)\lambda e^{-\lambda x}\mathrm{d}x \\&\qquad + \int _{{\bar{u}}_{(i,l,r_0)}}^{\infty }{\mathbb {P}}(\tau _q \ge u_{(1,k,r)})\lambda e^{-\lambda x}\mathrm{d}x \Bigg )f_{(i,l,r_0)}(0). \end{aligned}$$

Also,

$$\begin{aligned}&\frac{1}{p_{(1,k,r)}p_r}\psi _{(i,l,r_0)}^{(1,k,r)}\, (1-{{\bar{f}}_{(i,l,r_0)} }){\bar{\pi }}_{(i,l,r_0)} \\&\quad =\,{\mathbb {P}}\left( \tau _q \ge \max \{u_{(1,k,r)}-{\bar{u}}_{(i,l,r_0)},0\}\right) \left( \pi _{(i,l,r_0)}p_{(i,l,r_0)}p_{r_0}-f_{(i,l,r_0)}(0)\right) . \end{aligned}$$

Combining the above, it follows that \(P_{(1,k,r)}\) equals, with \(v^{(1,k,r)}_{(i,l,r_0)}\) as defined above,

$$\begin{aligned}&\sum _{r_0=1}^R \sum _{l=1}^M\sum _{i=1}^{N}\Bigg (1-e^{-\lambda {\check{v}}^{(1,k,r)}_{(i,l,r_0)}} -e^{-q{\hat{v}}^{(1,k,r)}_{(i,l,r_0)}}+\frac{\lambda }{\lambda +q}e^{-(\lambda +q){\check{v}}^{(1,k,r)}_{(i,l,r_0)}-qv^{(1,k,r)}_{(i,l,r_0)}} \\&\quad +\frac{q}{\lambda +q}e^{-\lambda (u_{(i,l,r_0)}-\Delta _{r_0})-qu_{(1,k,r)}} \Bigg )f_{(i,l,r_0)}(0)+\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N}e^{-q{\hat{v}}^{(1,k,r)}_{(i,l,r_0)}}\pi _{(i,l,r_0)} p_{(i,l,r_0)}p_{r_0}. \end{aligned}$$

These expressions can be written in a more convenient form. To this end, it is noted that

$$\begin{aligned}&\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N}e^{-q {\hat{v}}^{(1,k,r)}_{(i,l,r_0)}}\pi _{(i,l,r_0)}p_{(i,l,r_0)}p_{r_0} \\&\quad =\sum _{r_0=1}^R\sum _{l=1}^Me^{-q{\hat{v}}^{(1,k,r)}_{(1,l,r_0)}} P_{(1,l,r_0)}p_{(1,l,r_0)}p_{r_0} +\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=2}^{N}e^{-q{\hat{v}}^{(1,k,r)}_{(i,l,r_0)}}\pi _{(i,l,r_0)}p_{(i,l,r_0)}p_{r_0} \\&\quad =\sum _{r_0=1}^R\sum _{l=1}^Me^{-q{\hat{v}}^{(1,k,r)}_{(1,l,r_0)}}P_{(1,l,r_0)}p_{(1,l,r_0)}p_{r_0} \\&\qquad +\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=2}^{N}e^{-q{\hat{v}}^{(1,k,r)}_{(i,l,r_0)}}p_{(i,l,r_0)}p_{r_0}P_{(i,l,r_0)} \prod \limits _{m=1}^{i-1}\left( 1-\sum _{l_m=1}^Mp_{(m,l_m,r_0)}P_{(m,l_m,r_0)}\right) . \end{aligned}$$

After some basic algebraic manipulations, and using (4.9) we obtain

$$\begin{aligned}&\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=1}^{N}e^{-q {\hat{v}}^{(1,k,r)}_{(i,l,r_0)}}\pi _{(i,l,r_0)}p_{(i,l,r_0)}p_{r_0}=\sum _{r_0=1}^R\sum _{l_1=1}^M\Bigg ( e^{-q{\hat{v}}^{(1,k,r)}_{(1,l_1,r_0)}}\nonumber \\&\quad -\sum _{l=1}^M\sum _{i=2}^{N}p_{(i,l,r_0)}e^{-q(u_{(i,l,r_0)}+{\hat{v}}^{(1,k,r)}_{(i,l,r_0)})}\prod \limits _{m=2}^{i-1} \left( 1-\sum _{l_m=1}^Mp_{(m,l_m,r_0)}e^{-qu_{(m,l_m,r_0)}}\right) \Bigg )P_{(1,l_1,r_0)}p_{(1,l_1,r_0)}p_{r_0} \nonumber \\&\quad +\sum _{r_0=1}^R\sum _{l=1}^M\sum _{i=2}^{N}p_{r_0}p_{(i,l,r_0)}e^{-q(u_{(i,l,r_0)}+{\hat{v}}^{(1,k,r)}_{(i,l,r_0)})}\prod \limits _{m=2}^{i-1}\left( 1-\sum _{l_m=1}^Mp_{(m,l_m,r_0)}e^{-qu_{(m,l_m,r_0)}}\right) . \end{aligned}$$
(4.14)

Using Eq. (4.14) in (4.13), and after simplification, we obtain the linear system of equations (4.10), which proves Lemma 2. \(\square \)

The expressions for \(A^{(j,l,r_1)}_{(i,k,r_0)}(z)\) and \(A^{*(j,l,r_1)}_{\,\,(i,k,r_0)}(z)\), required for the queue-length analysis in Sect. 3, follow from substituting \(s=\lambda (1-B(z))\) in the two main terms in Eq. (4.7):

$$\begin{aligned} A^{(j,l,r_1)}_{(i,k,r_0)}(z)&= \, {\tilde{G}}_{(j,l,r_1)}(\lambda (1-B(z)),{\bar{u}}_{(i,k,r_0)}),\\ A^{*(j,l,r_1)}_{\,\,(i,k,r_0)}(z)&= \int _{x=0}^{{\bar{u}}_{(i,k,r_0)}}\lambda e^{-\lambda x}{\tilde{G}}_{(j,l,r_1)}(\lambda (1-B(z)),{\bar{u}}_{(i,k,r_0)} -x)\mathrm{d}x \\&\quad + {\tilde{G}}_{(j,l,r_1)}(\lambda (1-B(z)),0) e^{-\lambda {\bar{u}}_{(i,k,r_0)}}. \end{aligned}$$

Remark 3

We now discuss Assumption 1 in more detail, and how it plays a role in our result. The assumption entails that the gap \(T_{(i,r_0)}-\Delta _{r_0}\) left by a driver of profile \(r_0\), while being successful in the ith attempt, can be used by the subsequent driver only. If we would allow more than one following vehicle use a gap left behind by a merging vehicle, the analysis becomes much more complicated due to the fact that we need to distinguish between all possible cases where multiple drivers fit into this lag.

For example, in the proof of Theorem 1 we use explicitly that, if the system was not empty at arrival time, the gap left behind by a predecessor of type \((j, l, r_1)\) is simply \({\bar{u}}_{(j,l,r_1)} = u_{(j,l,r_1)}-\Delta _{r_1}\). However, when the predecessor succeeded at his first attempt (\(j=1\)) and required a critical gap of \(u_{(1,l,r_1)}\) that is smaller than the remaining gap of the predecessor’s predecessor, the gap left behind by the predecessor is equal to

$$\begin{aligned} u_{(i,k,r_0)} - \Delta _{r_0} - \Delta _{r_1} > u_{(1,l,r_1)} - \Delta _{r_1}, \end{aligned}$$

meaning that the current vehicle and its predecessor both used the same gap \(u_{(i,k,r_0)}\). In this sense, Assumption 1 ensures the tractability of the model.

Assumption 1 has not been stated in terms of the model input parameters, and is therefore not straightforward to check. To remedy this, we give an equivalent definition in terms of the critical headway and the merging time. Assumption 1 is satisfied if and only if

$$\begin{aligned} u_{(1,l,r_1)}\ge {\bar{u}}_{(i,k,r_0)}= u_{(i,k,r_0)}-\Delta _{r_0}, \end{aligned}$$
(4.15)

for all \(i\in \{1,2,\dots ,N\}\), \(k,l\in \{1,2,\dots ,M\}\), and \(r_0,r_1\in \{1,2,\dots ,R\}\). In the next section, we will numerically show that in cases that Assumption 1 is not fulfilled, the estimated capacity of the minor road is slightly smaller, but still highly accurate. The reason why the estimated capacity is a lower bound for the true capacity, is the fact that we waste some capacity by not allowing more than one vehicle to use the remaining part of a critical gap. The astute reader might have noticed that the maximum operators in Lemma 1 are not really needed due to condition (4.15). It turns out, however, that by preventing the corresponding terms from becoming negative, the capacity is approximated much better even when the condition is violated.

5 Numerical results

The analysis from the previous sections facilitates the evaluation of the performance of the system, including the assessment of the sensitivity of the capacity when varying model parameters. In this section, we present numerical examples to demonstrate the impact of these model parameters and of different driver behavior.

5.1 Example 1

In this illustrative example, we distinguish between two driver profiles. Profile 1 represents ‘standard’ traffic, whereas profile 2 can be considered as ‘slower’ traffic (for example large, heavily loaded vehicles). We assume that the ratio between profile 1 and profile 2 vehicles is \(90\%/10\%\). The fact that profile 1 vehicles are faster than profile 2 vehicles, is captured in their merging times (respectively \(\Delta _1\) and \(\Delta _2\)) and in the length of their critical gaps. From the profile 1 drivers, we assume that 40% need a gap of at least 5 s, upon arrival at the intersection, while the remaining 60% need a gap of at least 6 s. If, however, the drivers do not find an acceptable gap right away, they will grow impatient and be more and more prepared to accept slightly smaller gaps. We introduce an impatience rate \(\alpha \in (0,1)\) that determines how fast critical gaps decrease. Profile 2 vehicles are assumed to be slower, meaning that they have a need critical gaps of respectively 8 s (50%) or 9 s (50%) at their first attempt. Summarizing, we have the following model parameters in this example:

$$\begin{aligned} \begin{array}{|r|r|} \hline \text {Profile 1} &{} \text {Profile 2} \\ \hline p_1 = 0.9 &{} p_2=0.1\\ \hline p_{(i,1,1)}=0.4 &{} p_{(i,1,2)}=0.5 \\ p_{(i,2,1)}=0.6 &{} p_{(i,2,2)}=0.5 \\ \hline u_{(1,1,1)}= 5.0 &{} u_{(1,1,2)}= 8.0 \\ u_{(1,2,1)}= 6.0 &{} u_{(1,2,2)}= 9.0 \\ \hline \end{array} \end{aligned}$$

and, due to the impatience, we have for both profiles:

$$\begin{aligned} u_{(i+1,k,r)}=\alpha (u_{(i,k,r)}-\Delta _r)+\Delta _r, \quad i=1,2,\dots ,N-1,~ k=1,2,~ r=1,2. \end{aligned}$$
(5.1)

We vary \(\alpha \), \(\Delta _1\) and \(\Delta _2\) to gain insight into the impact of these parameters on the capacity of the minor road, while also varying q, the arrival rate on the major road. We take \(\alpha \in \{0.6,0.8,0.9,1.0\}\), \(\Delta _1\in \{4,5\}\) and \(\Delta _2\in \{5,6,7\}\) and vary q between 0 and 1000 vehicles per hour. Note that \(\alpha =1\) corresponds to a model without impatience. The results are depicted in Fig. 2a, b. In Fig. 2a, we can observe how the capacity of the minor road increases when merging time corresponding to at least one of the driver profiles decreases. From Fig. 2b we conclude by how much the capacity of the minor road increases when drivers become more impatient. The biggest capacity gain is caused by a decrease in \(\Delta _1\) from 5 to 4, obviously, because profile 1 constitutes the vast majority of all vehicles.

Fig. 2
figure 2

Capacity of the minor road (veh/h) as a function of the flow rate on the major road (veh/h) in Example 1

Full size image

So far we have only studied the impact of the model parameters on the capacity of the minor road. In practice, queue length and delay distributions are other important performance measures. In this numerical example, we compute the distribution of the minor road queue length and show how it depends on the batch size distribution. In our current example, we fix \(q=200\), \(\alpha =0.7\), \(N=10\), and we take an arrival rate of 300 vehicles per hour, arriving in batches of two vehicles on average. We compare the case where every batch consists of exactly two vehicles with random batch sizes where one, two, or three vehicles arrive simultaneously, each with probability 1 / 3. The results are shown in Table 1. We observe that the mean queue length depends on the distribution of the batch sizes, whereas the capacity only depends on the mean batch size. In this experiment, the differences are relatively small, but still the probability of having a queue length of more than five vehicles is almost twice as large (0.029) for random batch sizes than for fixed batch sizes (0.017). For other instances, the difference between the performance metrics in the two columns can be substantially more pronounced.

Table 1 Numerical results for Example 1
Full size table

Remark 4

Obtaining the numerical results for the queue lengths turned out to be significantly harder than computing only the service-time related performance measures (such as the capacity). We had to apply several numerical enhancements to the standard procedures. The most efficient way (that we found) to solve the system of equations (3.8) was taking Taylor series approximations of the \(A^{(j,l,r_1)}_{(i,k,r_0)}(z)\) and \(A^{*(j,l,r_1)}_{\,\,(i,k,r_0)}(z)\), and using a mix of numerical procedures and symbolic computations.

5.2 Example 2

The model parameters in Example 1 are such that condition (4.15) is met. In practice, however, situations might arise where this is not the case. The purpose of this example is to quantify errors made when computing the capacity using the methods proposed in Sect. 4, since these methods rely on condition (4.15). We take the same settings as in Example 1, but we adjust the critical gaps for profile 2 drivers from 8 and 9 to, respectively, 10 and 12 s. We keep the original merging times for both profile drivers of \(\Delta _1=4\) s and \(\Delta _2=5\) s. In order to quantify the errors, we use both the analysis from Sect. 4 and a computer program that simulates the model and gives accurate results. Table 2 compares the capacities obtained with both methods, for various values of q and \(\alpha \). It immediately becomes apparent that the relative error is below 0.5% in all cases. We conducted similar experiments for different settings and in all cases, we obtained relative errors below 2%, which clearly justifies using the analysis of this paper in all practical situations, including those violating (4.15). As expected (see Remark 3), the approximated capacities are (slightly) underestimating the actual capacities.

Table 2 Capacity of the minor road (veh/h) for various values of the flow rate, q, on the major road (veh/h) in Example 2
Full size table

5.3 Example 3

We demonstrate our methodology in an example with realistic parameters based on the empirical data study by Tupper (2011), who recorded data on more than ten thousand critical gaps in Massachusetts and Oregon. In this study, it is argued that the standard correction factors suggested in the Highway Capacity Manual (HCM 2010) do not always adequately describe the true impact of various factors on the critical gap. In Tupper (2011) the factors that influence this critical gap are studied; in this example, we focus on vehicle types (car, van, SUV, or Truck) and the age of the driver (teen, adult, or elderly). We base the critical gaps in this example on Tupper (2011, Chapter 6); for ease, we consider the case of no interaction between the driver type and vehicle type (but evidently any type of interaction could be included). In our experiment, we consider a situation in which 10% of the drivers are teens, 70% are adults, and 20% are elderly people; in addition, 70% of the vehicles are cars, while vans, SUVs, and trucks each constitute 10% of the total traffic. As a result, we have \(R=3\times 4=12\) driver profiles. We matched the averaged values as much as possible with data from Tupper (2011).

Table 3 Input data for Example 3: mean critical headways (in seconds) for the first attempt per combination of vehicle type and gender and their (weighted) averages
Full size table

We now explain how we incorporate impatience and random variations in the model. Denote the mean critical headways for the first attempt (see Table 3) by \(u_r\), \(r=1,2,\dots ,R\). Due to impatience, the mean critical headway should decrease from 6.5 s (first attempt) to 5.5 s (second attempt), 5.25 s (third attempt), and 5 s (all subsequent attempts). For this purpose, we define an impatience effect \(\beta _1=0\), \(\beta _2=6.5-5.5=1\), \(\beta _2=6.5-5.25=1.25\), and \(\beta _n=6.5-5=1.5\) for \(n=3,4,\dots ,N\). We take the maximum number of attempts \(N=100\), which is more than sufficient to obtain accurate values for the capacities. Finally, we incorporate randomness by adding \(-\,1\), 0, or 1, each with probability 1 / 3, unless it would result in a critical headway less than 2.5 s, in which case we take 2.5. These effects combined result in the following random critical gaps:

$$\begin{aligned} T_{(i,r)} = \max (2.5, u_r-\beta _i+\epsilon ), \end{aligned}$$

with

$$\begin{aligned} \epsilon = {\left\{ \begin{array}{ll} -\,1 &{} \quad \text {w.p. }1/3,\\ 0 &{} \quad \text {w.p. }1/3,\\ 1 &{} \quad \text {w.p. }1/3,\\ \end{array}\right. } \end{aligned}$$

for \(i=1,\dots ,N\) and \(r=1,\dots ,R\). The values for \(\Delta _r\) are chosen to be equal to the lowest possible critical headway that can be realized per driver profile:

$$\begin{aligned} \Delta _r=\min _{i,k} u_{(i,k,r)}=\max (2.5, u_r-\beta _N-1)=\max (2.5, u_r-2.5), \end{aligned}$$

which varies from 2.5 to 5.625 s. The capacities for this example are given in Table 4. We have used two methods to compute these capacities:

  1. (1)

    The full model discussed in this paper, with profiles, impatience, and randomness;

  2. (2)

    An aggregated model where critical headways are averaged over all profiles.

The reason to compare these two models is to check whether one could also determine the capacity by just averaging over all driver profiles and use simpler models; a procedure that was shown in the HCM to perform well in certain circumstances. For a fair comparison, we have taken \(\Delta =\Delta _1=4.0175\) s in model 2, which is equal to the value found in the full model. The results in Table 4 indicate that the aggregated model overestimates the capacity. For small values of q, the approximation error decreases and even disappears in the limit \(q\downarrow 0\), where the capacity approaches \(3600/\Delta =896.1\) vehicles per hour, as we have witnessed in Example 1. We conclude that for these parameter values, indeed, the simpler model is a reasonable approximation for the full model. This may be true for these parameters, but not in general, and therefore (if the increased computational complexity is no issue) we always recommend using the full model.

Table 4 Capacities (veh/h) for Example 3
Full size table

6 Discussion and concluding remarks

In this paper, we have presented a gap acceptance model for unsignalized intersections that considerably generalizes the existing models. The model proposed incorporates various realistic aspects that were not taken care of in previous studies: driver impatience, heterogeneous driving behavior, and the service time being dependent on the vehicle arriving at an empty queue or not. Despite the rather intricate system dynamics, we succeed in providing explicit expressions for the stationary queue-length distribution of vehicles on the minor road, which facilitate the evaluation of the corresponding capacity (of the minor road, that is).

We have concluded our paper by presenting a series of numerical results, which are representative of the extensive experiments that we performed. Our techniques facilitate the quantitative evaluation of the system at hand, including the assessment of the sensitivity of the capacity when varying the model parameters (which is typically considerably harder when relying on simulation).

An interesting direction in which to pursue this line of research, would be to allow a more general arrival process on the major road. In the present paper, we assume that the arrival process on the major road is a Poisson process. In practice, however, there might be platoon forming on this road and several papers have shown that this clustering of vehicles will influence the capacity of the minor road (cf. Abhishek et al. 2019; Tanner 1962; Wegmann 1991; Wu 2001). It would certainly be interesting to combine the framework in this paper with, for example, the gap-block arrival process in Cowan (1975) or Markov platooning introduced in Abhishek et al. (2019).