Comparative Accuracy of a Data-Fitted Generalized Aw-Rascle-Zhang Model

Comparative Model Accuracy of a Data-Fitted Generalized Aw-Rascle-Zhang Model

Shimao Fan Department of Civil and Environmental Engineering

University of Illinois at Urbana Champaign
205 N. Mathews Ave
Urbana, IL 61801
Michael Herty Department of Mathematics
RWTH Aachen University

Templergraben 55
D-52056 Aachen
 and  Benjamin Seibold Department of Mathematics
Temple University

1805 North Broad Street
Philadelphia, PA 19122

The Aw-Rascle-Zhang (ARZ) model can be interpreted as a generalization of the Lighthill-Whitham-Richards (LWR) model, possessing a family of fundamental diagram curves, each of which represents a class of drivers with a different empty road velocity. A weakness of this approach is that different drivers possess vastly different densities at which traffic flow stagnates. This drawback can be overcome by modifying the pressure relation in the ARZ model, leading to the generalized Aw-Rascle-Zhang (GARZ) model. We present an approach to determine the parameter functions of the GARZ model from fundamental diagram measurement data. The predictive accuracy of the resulting data-fitted GARZ model is compared to other traffic models by means of a three-detector test setup, employing two types of data: vehicle trajectory data, and sensor data. This work also considers the extension of the ARZ and the GARZ models to models with a relaxation term, and conducts an investigation of the optimal relaxation time.

Key words and phrases:
traffic model, Lighthill-Whitham-Richards, Aw-Rascle-Zhang, generalized, second order, fundamental diagram, trajectory, sensor, data, validation
2000 Mathematics Subject Classification:
35L65; 35Q91; 91B74

1. Introduction

The mathematical modeling of vehicular traffic flow knows a variety of types of descriptions (see [38, 7] for review papers): microscopic (e.g., [61, 54, 4, 37]), which model the individual vehicles and their interactions by ODE; cellular (e.g., [52, 29, 22, 65, 1]), which divide the road into cells and prescribe stochastic rules how vehicles advance through cells; and continuum. This latter class divides into mesoscopic/gas-kinetic (e.g., [39, 60, 46, 42, 41, 40]) and macroscopic (e.g., [49, 63, 69, 58, 59, 47, 44, 45, 20, 3, 22, 30, 31, 8, 6]), i.e., fluid-dynamical, models. Among the (inviscid) macroscopic models one distinguishes between first-order models based on scalar hyperbolic equations and second-order models comprised of systems of hyperbolic equations. Specific examples for the latter are the Payne-Whitham model [58, 72], two-phase models [16], and the Aw-Rascle-Zhang model [3, 33, 74, 8].

All of these types of traffic models are of practical importance, however—due to their different mathematical structure—for different purposes. For instance, microscopic models are well-suited for traffic simulation, i.e., the “in silico” study of a specific scenario; cellular models reproduce jamming behavior while being simple to implement and easy to parallelize; mesoscopic models provide a high modeling flexibility; and macroscopic models provide a suitable framework for the incorporation of on-line data. Moreover, there are mathematical relations between these types of models. For example, microscopic models as well as cellular models converge to mesoscopic or macroscopic models in the limit of vanishing cell size or vehicle spacing, respectively [19, 22, 1, 12]. Similarly, macroscopic models arise as suitable limits of mesoscopic models [53, 46].

Microscopic and cellular models are nowadays widely used in traffic engineering, and their combination with data is ubiquitous. In contrast, certain types of continuum models have been studied mathematically, but very little work has been conducted on their validation with traffic data. Examples of macroscopic first-order models used in traffic engineering practice are the Mobile Century project [2] and the Mobile Millennium project [51], including approaches based on the reformulation of the Lighthill-Whitham-Richards model (2) in terms of Hamilton-Jacobi equations [5, 6], or in terms of the velocity variable [73]. Further examples are [9, 10] and the references therein. Those projects focus on the assimilation of data, the reconstruction of traffic states (e.g., from cell phone data), and the combination of macroscopic models with filtering techniques.

In contrast, here different types of traffic models are generated via historic data, and then their predictive accuracy is investigated using time-dependent data. Moreover, the main focus lies on second-order models; and first-order models are considered mainly for comparison purposes. The contributions of this paper are: (i) the design of a generalized Aw-Rascle-Zhang model and the analysis of its mathematical properties; (ii) a systematic methodology to construct data-fitted first-order and second-order traffic models, using historic fundamental diagram data; (iii) the validation of first- and second-order macroscopic traffic models via time-dependent trajectory and sensor data, and the comparison of the predictive accuracy of different models; and (iv) the investigation of the optimal relaxation time in second-order models with a relaxation term.

This paper is organized as follows. In §2 an overview over existing macroscopic traffic models is provided, including a discussion of some of the modeling shortcomings of the Aw-Rascle-Zhang model. We then introduce the generalized Aw-Rascle-Zhang model as an approach that addresses the shortcomings, and discuss its mathematical properties. The fitting of the model parameters and functions is then described in §3. Given historic fundamental diagram data in the flow rate vs. density plane, we systematically construct data-fitted first- and second-order macroscopic models. In §4 the numerical methods used to conduct the model validation and comparison are presented. Unlike studies of cell-transmission models [19], in this paper all studies are carried out in a macroscopic sense; in particular the governing PDE are numerically solved with high enough accuracy such that the numerical approximation errors are negligibly small relative to the model errors. The comparison of the models on a three-detector test setup [21] is then carried out in §5. In addition to the macroscopic traffic models, we also consider a predictor that simply interpolates the traffic state from the boundaries. For vehicle trajectory data, and for sensor data, the predictive accuracies of the models are compared, and their reproduction of features in the traffic states are studied. In §6 we then extend the studies to data-fitted second-order models with relaxation terms. In particular, we study the dependence of the model accuracy on the relaxation time at which drivers adjust their driving behavior. Finally, in §7 we present the conclusions from our studies.

2. Existing and New Macroscopic Traffic Models

Common to all macroscopic traffic models is the continuity equation


which gives the conservation of vehicles. In (1), the vehicle density is , and the vehicle velocity field is , where is the position along the road, and is time. If the road has multiple lanes (in a given direction), we consider these aggregated into the scalar field quantities and .

2.1. The Lighthill-Whitham-Richards Model and Flow Rate Functions

The simplest macroscopic traffic model, the Lighthill-Whitham-Richards (LWR) model [49, 63], is obtained by assuming a functional relationship between and , i.e., . This turns equation (1) into a scalar hyperbolic conservation law


where the flux is given by the flow rate function . Because the LWR model (2) is a closed model consisting of a single equation, it is denoted a first order model. The velocity function is commonly assumed to be decreasing in with for some maximal vehicle density .

Popular examples of flow rate functions are the Greenshields flux [35], in which is a quadratic function, and the Newell-Daganzo flux [55, 19], in which is a piecewise linear function. While these different choices of functions lead to well-posed first-order models, the second-order models derived below call for further properties that the function must satisfy. In particular, the velocity function must nowhere be constant, because otherwise hyperbolicity would be lost (see the analysis in §2.5.2). This rules out the Newell-Daganzo flux; and as a consequence, in this paper we consider flow rate functions (12) that resemble the shape of the Newell-Daganzo flux, but that are strictly concave.

2.2. The Aw-Rascle-Zhang Model

The strict functional relationship between and is loosened in second order models, which augment (1) by an evolution equation for the velocity field. Payne and Whitham proposed a model [58, 72] in which the vehicle velocity relaxes towards a velocity function, while also being affected by a “traffic pressure” that is analogous to a pressure in fluid dynamics. Because this pressure can lead to unrealistic solutions (vehicles going backwards on the road, shocks that overtake vehicles, etc., see [20]), Aw and Rascle [3], and independently Zhang [74], proposed a different form of “pressure” that remedies the shortcomings of the Payne-Whitham model. The homogeneous Aw-Rascle-Zhang (ARZ) model reads as


where we call the hesitation function.111The function is sometimes called “pressure”, and denoted , even though it does not play the role of a pressure in the equations. We assume that and use the gauge . The addition of a relaxation term (analogous to the Payne-Whitham model) yields the inhomogeneous ARZ model [33, 62]


We call the desired velocity function or the equilibrium velocity function, and the relaxation time scale.

As one can easily verify, the homogeneous ARZ model (3) possesses no mechanism to make drivers move when starting with all vehicles at rest, i.e., . In turn, the inhomogeneous ARZ model (4) does. We therefore expect the homogeneous ARZ model to yields reasonable results only when the traffic flow is close to its equilibrium state, i.e., . Yet, in general the inhomogeneous ARZ model has the potential to yield more realistic predictions.

Remark 1.

The conservative form of (4) is given by


where the two conserved variables are and , and is called the equilibrium curve that the momentum density relaxes to.

Remark 2.

Various authors (e.g., [33, 67]) have proposed to choose the functions and dependent on each other, namely . In this case, the solution relaxes towards an equilibrium state (see §2.5.3). In this paper, we follow this philosophy, since it generates both functions from the same data-fitting procedure. However, it should be noted that it is in principle perfectly reasonable to choose and independently of each other. As shown in [34, 28, 66, KasimovRosalesSeiboldFlynn2013], such a choice can generate (whenever ) instabilities and self-sustaining traveling wave solutions that model phantom traffic jams and traffic waves, respectively.

2.3. Interpretation of ARZ as Generalization of LWR

As pointed out in [47, 3, 8, 26], the homogeneous ARZ model (3) can be interpreted as a generalization of the LWR model, by introducing the variable . Thus, system (3) takes the form


The interpretation of (6) is that is advected with the flow , i.e., it moves with the vehicles. One can therefore interpret as a quantity that is associated with each vehicle, and that is influencing the velocity. Since for , we have and thus , we call the empty road velocity. The actual vehicle velocity is given by its empty road velocity, reduced by the hesitation .

Using this interpretation, the homogeneous ARZ model generalizes the LWR model (2), as follows. Given a (decreasing) LWR velocity function , we define (clearly, and ). Then, model (6) possesses a one-parameter family of velocity curves, namely , and the LWR velocity curve is one of them, namely the one corresponding to . The same behavior translates to the flow rate curves that live in the fundamental diagram ( vs. ). The ARZ model (6) possesses a one-parameter family of flow rate curves, namely , and the LWR flow rate curve is one of them, namely the one corresponding to . This has been observed by Lebacque [47].

The aforementioned relationship between the LWR and the ARZ model is shown in Fig. 2. The single velocity curve (left panel) and flow rate curve (right panel), shown in red, is one representative of the family of curves, shown in black, that the ARZ model possesses. By construction, each velocity curve in the ARZ model is merely a vertical translation of the LWR velocity curve. Hence, the homogeneous ARZ model possesses the same number of parameter functions (namely: a single one) as the LWR model.

In line with Remark 2, i.e., by choosing , we can extend model (6) to the inhomogeneous case, yielding


which adds a temporal relaxation of each vehicle’s empty road velocity towards a uniform value . In other words: the dynamics, that can be on any velocity curve , are driven towards the LWR velocity function , i.e., the red curves in Fig. 2—unless the system is driven away from equilibrium by another effect, such as boundary conditions (see §2.5.3).

2.4. Generalized ARZ Model

The interpretation of the ARZ model (3) as possessing a family of velocity curves (see (6)) reveals a fundamental shortcoming of the model: due to the additive relationship between velocity, empty road velocity, and hesitation, there is not a unique maximum density , at which the flow stagnates. On the contrary, as the plots in Fig. 2 indicate, variations in can lead to significant variations in the density at which . However, since in reality the maximum density is largely a property of the road, it should not depend (at least not strongly) on the velocity that drivers assume when alone on the road. In order to remedy this shortcoming, the relationship between , , and must be generalized.

To that end, we consider a generalized Aw-Rascle-Zhang (GARZ) model, which is a representative of the class of generic second order models (GSOM), proposed by Lebacque, Mammar, and Haj-Salem [48]. Specifically, the homogeneous ARZ model (6) generalizes to


where we impose the following requirements on the velocity function , and the associated generalized flow rate function :

  • , i.e., vehicles never go backwards on the road.

  • , i.e., we gauge the convected quantity to play the role of the empty road velocity, as in the ARZ model (6).

  • for , i.e., each flow-rate curve is strictly concave. This condition implies (see Lemma 1) in particular that , i.e., each velocity curve is strictly decreasing w.r.t. the density.

  • , i.e., a faster empty road velocity results in a faster velocity for all possible densities.

  • , i.e., for , the concavity of and the slope of hold with an equality sign.

Lemma 1.

Consider a function , and let . If everywhere, then everywhere.


The function satisfies: (i) , and (ii) everywhere. Hence, everywhere. ∎

In order to define an inhomogeneous GARZ model, an equilibrium velocity curve must be specified. We assume that it is a member of the family of velocity curves defined by , i.e.

for some equilibrium empty road velocity . We choose to generalize (4) to the GARZ case as follows:


Note that it would alternatively be conceivable to propose a relaxation in (9) of the form


While for the ARZ model, the forms (4) and (7) are equivalent, for the GARZ model the relaxations (9) and (10) are not. Specifically, if , then both forms relax to the same limit, but at different rates. A simple Taylor expansion yields that for nearby , the relaxation in (9) happens times as rapidly as the relaxation in (10).

Finally, in line with (5), the GARZ model (9) is meant to be interpreted in the conservative form


where the two conserved variables are and . Moreover, and .

2.5. Properties of the GARZ Model

Most properties of the classical ARZ model transfer over to its generalization, the GARZ model. Here we only collect relevant results, many of which have been presented in [48], or that are relatively straightforward generalizations of the results given in [3, 33, 62]. The theoretical results presented below are in particular important for the data-fitting methodologies conducted in §3, and for the interpretation of the results obtained in §6.

2.5.1. Regions of GARZ Variables and Inverse Velocity Functions

Because of their relations and because of their physical meaning, the quantities , , and cannot assume any arbitrary values. In this paper, we assume that there is a unique stagnation density , at which vehicles come to a stop, independent of their empty road velocity, i.e., for all . We therefore have and , where the latter inequality follows from . Moreover, we assume that there is a minimum and maximum empty road velocity, i.e., . Because and , the function can be “inverted” to define the functions (and their domains):

Here the functions and are defined as the unique solutions to the problems:

  1. given and , find , s.t. ;

  2. given and , find , s.t. .

From the fact that the quantity is transported with the flow (while possibly relaxing to some ), and from the solution of the Riemann problems of the GARZ model (see below), it follows that the dynamics of the model never generate values or . Hence, analogous to the ARZ model (cf. [3]), the domain is an invariant region.

2.5.2. Characteristics and Associated Fields

The homogeneous part of the GARZ model (11) is a conservation law of the form


The Jacobian of the flux function is

and its eigenvalues and associated eigenvectors are


Hence, like the ARZ model (5), the GARZ model (11) is strictly hyperbolic for . One of its characteristic velocities, , is slower than the vehicles (i.e., ) and its associated field is genuinely nonlinear (i.e., it corresponds to shocks and rarefaction waves, see below). Its other characteristic velocity equals the vehicle velocity and its associated characteristic field is linearly degenerate (i.e., its associated waves are contact discontinuities that are transported with the flow).

Lemma 2.

Both characteristic velocities are strictly decreasing w.r.t. , i.e., and .


We have that . Since is assumed concave w.r.t. , it follows that . Moreover, by Lemma 1. ∎

We continue with the discussion of the characteristic fields. The scalar function satisfies , and it is a Riemann invariant to . Hence, across waves of the first family, the empty road velocity is constant. The field associated with the second eigenvalue is linearly degenerate. It is given by and across waves of the second family, the velocity is constant. The solution to a Riemann problem, i.e., the Cauchy problem to system (11) on the real line with discontinuous piecewise constant initial data , where is the Heaviside function, generalizes naturally from the ARZ model as well. In general, the solution is obtained by superposition of simple waves connecting different constant states: from a given left state to a given right state via an intermediate state that is connected to by a 1-wave (i.e., a Lax-shock or rarefaction associated to the first characteristic field), and to by a 2-wave (i.e., a contact discontinuity). Since the GARZ system (11) is—as the ARZ model—of Temple class [68], shocks and rarefaction wave curves in phase space coincide. Moreover, due to Lemma 2, a simple wave of the first family is either a shock or a rarefaction wave.

In the phase space -plane we discuss the shape of the characteristic fields. The second fields are parallel to the -axis, and the first fields are the contours . Since by assumption and , the contours of always have a finite and truly positive slope in the -plane. Thus, for any two states and that satisfy , where , there is a unique intermediate state , defined via the inverse function given in §2.5.1. Moreover, because is decreasing with (see Lemma 2), the Lax entropy conditions [23] imply that for the 1-wave is a shock wave (moving with speed , given by the Rankine-Hugoniot conditions [23]), while for it is a rarefaction wave. The condition means that drivers on the left wish to drive at least as fast as the vehicles on the right are driving. If this is not the case, i.e., if , then there is no non–negative density at which the 1-wave and the 2-wave intersect. Here, a vacuum state will be generated, analogously to the construction for the ARZ model [3, 62]. The left state is connected by a rarefaction wave to a left vacuum state ; this state is connected to a right vacuum state via another rarefaction (which is feasible because ); and this then connects to the right state via a 2-contact discontinuity.

2.5.3. Relaxation of GARZ to LWR

Smooth solutions to the Cauchy problem of the inhomogeneous GARZ model (9) relax in time towards solutions of the LWR model (2), because tends to along characteristic curves. Note that for general relaxation systems, convergence to a first-order equation is only warranted if a sub-characteristic condition is satisfied, cf. [50, 14, 66]. Here, we are in the characteristic case. Moreover, shocks of the hyperbolic system (11) are also shocks of the LWR model (2). Therefore, for the Cauchy problem, GARZ solutions converge to LWR solutions as , and/or as .

Figure 1. Riemann problem at right domain boundary. Consider a constant initial state in the domain, and prescribed boundary data with . The LWR model preserves the constant state (the boundary data is projected onto the equilibrium curve vertically along the dashed line). In contrast, the GARZ model generates a new state at the boundary (projected onto the equilibrium curve along a ray through the origin), from which a shock moves into the domain, thus changing the initial state to the new boundary state.

In contrast, for initial-boundary-value problems on a bounded domain , this last property is in general not true. To highlight this fact we consider a simplified setting depicted in Fig. 1. Let constant boundary data , , , be given. With this data, we can solve the LWR model (2) and the inhomogeneous GARZ model (11). Then, in general solutions to the latter problem do not converge to solutions of the former, even in the limit . Consider a constant state inside the domain. At the outflow boundary , let a state be given where . The LWR model only uses the density information, and thus the constant state is preserved. In contrast, the GARZ model uses the full state in the Riemann problem. Its solution yields an intermediate state , whose density is determined via the relation . This intermediate state connects to the boundary state via a contact discontinuity (a 2-wave) moving with speed , and to the interior state via a shock (a 1-wave) that moves with speed

which in the situation depicted in Fig. 1 moves into the domain, since . Thus, after some time in the GARZ model the intermediate state is observed within the domain. Note that this argument holds independent of the value of the relaxation time in the model (9).

3. Data-Fitted Traffic Models

In this section we describe how the parameter functions of the traffic models presented in §2 can be fitted to historic fundamental diagram data. We assume that flow rate vs. density pairs are given from long-term measurements (commonly obtained via stationary sensors). As visible in the right panel of Fig. 2, these data (gray dots) tend to exhibit a relatively clear functional relationship between and for low densities. In turn, for medium densities, a significant spread is visible, i.e., a single -value corresponds to many different flow rates . Finally, for large densities, very few data points are available at all.

3.1. Data-Fitting for the LWR and ARZ Models

The first-order LWR model (2) must represent these data via a single function . As the spread of the data cannot be captured, it is reasonable to find a function that lies “in the middle” of the cloud of data points. Specifically, we employ the approach presented in [26]. First, since the stagnation density is not represented well via data, we prescribe it as a fixed constant, given by a typical vehicle length of 5 meters, plus 50% of additional safety distance,

Second, a three-parameter family of smooth and strictly concave flow rate curves is selected as



Each flow rate function in this family vanishes for and . The three free parameters allow for controlling three important features of : the value of maximum flow rate (mainly determined by ), the critical density (mainly controlled by ), and the “roundness” of the curve, i.e., how rapidly the slope transitions from positive to negative near (dominated by ).

Third, from this three-parameter family of flow rate curves, the one is selected that is the closest to the data points in a least-squares sense, i.e. we solve


In the right panel of Fig. 2, the resulting least-squares fit to the given gray data points, called , is depicted by the red curve. The red curve in the left panel represents the resulting velocity function .

As described in §2.3, the ARZ model (6) generalizes the LWR model to a one-parameter family of velocity curves (black curves in the left panel of Fig. 2) and flow rate curves (right panel of Fig. 2). Due to this property, the ARZ model possesses the same amount data-fitted parameters as the LWR model.

An interpretation of the ARZ family of velocity curves is that different -values represent different types of drivers; the larger , the faster the corresponding drivers tend to drive. As motivated in §2.4, this captures the spread in the fundamental diagram (which is desirable), but it also results in vastly varying stagnation densities for different types of drivers. This last property is unrealistic, as the maximum density is a property of the road, rather than of the behavior of drivers [26]. We therefore need to construct a family of curves that are not simple shifts of each other, as done below.

Figure 2. Velocity vs. density (left panel) and flow rate vs. density (right panel) curves of the smooth three-parameter model (12), fitted with historic fundamental diagram data (gray dots), for the ARZ model.

3.2. Data-Fitting for the GARZ Model

The GARZ model (8) is based on a generalized velocity function , where—as for the ARZ model— represents different types of drivers, and thus parameterizes a families of velocity and flow rate curves, respectively. We construct these families of curves by generalizing the least-squares fit (13) to a weighted least-squares fit, as follows.

Given a weight parameter , we consider the minimization problem



For , problem (14) reduces to (13), i.e., the LWR equilibrium curve is recovered. For , data below the curve is penalized more, and consequently the resulting curve moves downwards. In turn, if , curves above the equilibrium curve are obtained.

The weighted least-squares problem (14) generates a one-parameter family of curves , parameterized by ; and consequently it also generates a family of velocity functions

In this paper, we restrict to the case that the velocity curves in the family are non-intersecting, i.e.,


Note that in general (i.e., for a general family of flow rate functions, and for general data points), property (15) is not necessarily satisfied; and furthermore, problem (14) need not have a unique solution. However, in all cases studied in this paper, the minimization problem (14) does have a unique solution; and property (15) is in fact satisfied.

Figure 3. Family of velocity vs. density (left panel) and flow rate vs. density (right panel) curves generated from the weighted least square (WLSQ) algorithm in constructing the velocity function in the GARZ model.

While problem (14) is defined for all , values of that are extremely close to 0 or 1 tend to lead to unreasonable curves. The reason is that in the limit , the resulting curve is the lowest curve that has no data points above it, and as a result it adjusts to outliers in the data (similar arguments hold for ). We therefore define a lower/upper flow rate curve, such that 99.9% of all data points lie above/below it. Consequently, together with the equilibrium curve, we have the following three special flow rate curves

where here we use and . In the right panel of Fig. 3 these three curves are depicted by the lowest black curve, the red curve, and the uppermost black curve, respectively.

Even though the parameter defines a family of flow rate curves as desired, it has the shortcoming that it does not have an immediate interpretation as a property of traffic flow. We therefore re-parameterize the family in terms of the empty road velocity , as follows. For any , we define as the resulting slope of the curve at , i.e.,

Due to property (15), the relationship is strictly increasing, and thus can be inverted into , defined on the interval , where and . Using this re-parameterization, we obtain a generalized flow rate function

and a generalized velocity function

as used in the GARZ model (8). The properties of assumed in §2.4 are satisfied by construction.

3.3. Domain Extension for the GARZ Model

The systematical construction of a generalized velocity function , presented in §3.2, is in line with the regions of the GARZ variables defined in §2.5.1. In particular, the function is defined only for . However, when applying the GARZ model in a forward computation, velocity data may be provided through initial and boundary conditions that lie outside of the domain of . In order to make sense of the model for such data, we effectively extend the domain of the function via a projection of such data, as follows.

Given a density–velocity pair , where , we define a projected velocity as

and thus obtain the extended function

that is defined for arbitrary velocity values. This simple projection (densities are left unchanged, and velocities are moved onto the lowest or highest curve, respectively) provides a constant extension of the function beyond its domain . Note that the range of remains unaffected as being , and consequently the function is not invertible outside of .

4. Numerical Methods

All models are approximated numerically using a finite volume method on a regular grid of cell size and time step , chosen so that the CFL condition [18]

is satisfied, where is the largest wave speed (see §2.5.2 for the characteristic velocities of the models). In all examples throughout this paper, the grid resolution is chosen small enough () so that the numerical approximation errors are much smaller than the model errors. Hence, the studies are conducted truly on the continuum level.

The first order model (2) is solved using Godunov’s method [32]. For the second order models, we have to account for the fact that the inhomogeneous GARZ model (11) becomes stiff if is small. Hence, we employ a semi-implicit finite volume scheme that treats the nonlinear hyperbolic terms explicitly and the relaxation terms implicitly (to prevent a time step restriction ). The update rule of a state in cell from time to the state at time reads as


Here denotes the numerical flux for the quantity through the boundary between cells and ; the other fluxes are defined accordingly. Moreover, is the model’s two-parameter flow rate function, and is the equilibrium flow rate function.

As in the Godunov method [32] one could use the exact solution to the Riemann problem (cf. [3, 24]) to define the numerical fluxes. However, this would require the inversion of the velocity function , which is costly for the GARZ model. A less expensive approach, employed here, is to define the numerical fluxes via the HLL approximate Riemann solver [36], which approximates the true Riemann problem by a single constant intermediate region. Note that due to the fine grid resolution, and due to the fact that initial and boundary conditions are continuous, spurious overshoots that may occur in the velocity (cf. [13]) are negligibly small.

Since the inhomogeneous model (11) possesses a relaxation only in the momentum equation, the time update of the density variable, given by (16), is fully explicit. Therefore, in the time update of the generalized momentum, given by (17), the quantity is known. Specifically, the numerical scheme is implemented in three steps:

  1. Based on the data at time , the fluxes are computed.

  2. The new density states are computed via (16).

  3. The new generalized momenta are computed according to (17). On the cell , the new state is obtained as the solution of the scalar nonlinear equation , where


    The root of (18) is found up to machine accuracy within a few Newton steps, using as the starting guess.

In the special case of the ARZ model (5), the update (17) is given by the explicit formula


It should further be remarked that the semi-implicit scheme (16) and (17) is equivalent to the fractional step approach that first approximates the homogeneous part of (11) via a forward Euler step, and then approximates the relaxation part via a backward Euler step. Hence, in the limit (while fixed), the scheme amounts to simply projecting onto the equilibrium curve in the relaxation step, i.e., equation (17) turns into . In turn, in the homogeneous case, i.e., , the relaxation terms are simply omitted.

The boundary data are provided by introducing a ghost cell adjacent to the outermost grid cell (on either side of the domain), in which the boundary state is assumed. The numerical fluxes in (16) and (17) then by construction pick up the information corresponding to waves that enter the computational domain.

5. Validation and Comparison of Models via Real Data

In the following, we validate the presented models by studying how well they reproduce the evolution of real traffic data, and in that process we compare the predictive accuracy of the models. A particular focus lies on the investigation of the extent to which the GARZ model, that addresses various shortcoming of traditional models, improves the actual model agreement with real traffic data. We conduct the validations using the NGSIM trajectory data set [71] and the RTMC sensor data set [56].

The test framework considered here is based on the methodology presented in [26] and further developed in [27]: on a segment of highway, a three-detector test problem [21] is formulated. At each end of the segment, the traffic state is (at all times) provided to the traffic model, which is advanced forward in time (using the numerical methods described in §4) inside the segment. The predictions that the traffic model produces in time are then compared to real data inside the segment, and the deviation between predicted and real traffic states is used to quantify the model error.

5.1. Treatment of Data

As described in [26], continuous field quantities and are constructed from the NGSIM vehicle trajectory data [70], using kernel density estimation [57, 64]. In this approach, given vehicle locations (including “ghost vehicle” positions, obtained via reflection at the boundaries, see [43]), density and flow rate functions are obtained as superpositions of Gaussian profiles,

and the velocity field is then given by . The kernel width is chosen  meters. These field quantities then define initial conditions () and boundary conditions (when evaluated at the segment boundary positions) for the traffic models, and reference states for the validation (inside the segment for ). Before the boundary data can be provided to the traffic model, one additional processing step must be applied to address spurious fast oscillations in the reconstructed boundary data (due to variations in the starting and end position of each vehicle trajectory in the data set): the time domain is divided into intervals of length 15 seconds, and on each interval the boundary data is replaced by a cubic polynomial that is a least-squares fit to the data, under the constraint that the resulting evolution is globally .

For the RTMC sensor data, vehicles densities and flow rates are given at three sensor positions, aggregated in intervals of length 30 seconds. Temporally continuous quantities and at a sensor position are generated via cubic spline interpolation (in time) of the aggregated information. One shortcoming of the RTMC data is the absence of a reliable initial state (because information is given only at the sensor positions). This problem is circumvented by running the models forward through an initialization phase (5 minutes), before the actual model comparison is started. During this phase, the boundary data has time to move into the domain and create a reasonable initial state for the actual validation.

5.2. Quantification of Model Errors

The quantification of the deviation of the model predictions from the actual data requires two aspects to be specified: first, which field quantities to consider and how to combine them into a single quantity; and second, if data is available at multiple positions and/or times, how to combine these multiple pieces of information into a single quantity?

Regarding the choice of field quantities, in this study we are interested in models that predict traffic densities (as required for instance for ramp metering) and velocities (as required for instance for travel time estimates) accurately. Since densities and velocities have different physical units, suitable normalization constants must be found, so that the deviations in each quantity contribute with equal influence to the total model error.

Given model predictions and , and real data and , we define a space-and-time-dependent error measure as


where the normalization constants and represent the ranges of the fundamental diagram data, as defined below. Note that various choices of normalization constants have been proposed in the literature. For instance, in [11] the absolute errors in density and velocity are scaled with and , respectively. A shortcoming of this choice is that for traffic flow at low densities, errors in density get divided by a very small number and thus significantly amplified. An alternative choice is employed in [26] by using and as normalization constants. However, these tend to give too much influence to velocity errors, because even in moving congested traffic flow, tends to be significantly smaller than .

Figure 4. Construction of ranges in density and velocity in historic fundamental diagram data, using an upper density (red line) and lower and upper velocity boundaries (blue lines). Data points with densities below a threshold (black line), as well as outliers, are systematically excluded.

In line with [27], we argue that balanced weights are given when the error in each field quantity is related to the maximum variation that the respective quantity exhibits in the historic fundamental diagram. In order to exclude the influence of outliers in the data, we conduct the following four-step approach. First, all data points with are neglected. The rationale is that these data do not contribute any meaningful information about the spread in the traffic states, and moreover such low density values are not meaningful in the context of a macroscopic description of traffic flow. In Fig. 4, this boundary is depicted by the black line. Second, similar to the method presented in [9], the upper density boundary is defined such that 99.9% of the remaining data points lie below it (red line in Fig. 4). Third, the lower (upper) velocity boundary () is defined such that 99.9% of the remaining data points lie above (below) it (blue lines in Fig. 4). Fourth, we define the data ranges

Regarding the norms and averages to measure the model errors, we use the following expressions. On a segment and time interval , spatial and spatio-temporal averages are considered


Moreover, for the RTMC data set, the temporal error at a sensor position inside the road segment on a given day is considered, as well averages over multiple days

Figure 5. Flow rate vs. density curves for the models LWR (left column), ARZ (middle column), and GARZ (right column), together with the traffic states observed in the test cases (gray dots). The results show the NGSIM data (top row) and the RTMC data (bottom row).

5.3. List of Models

We compare the following four models in terms of their predictive accuracy of the real data:

  1. Interpolation: A predictor that, at any instance in time, constructs the traffic density and velocity via direct linear interpolation of the boundary conditions, i.e., on the road segment , the predicted state is and . Of course, this predictor is not an actual traffic model. However, due to its simplistic nature, it represents an important means of comparison.

  2. LWR: The first-order model (2), in which only the density state is evolved, based on the data-fitted equilibrium velocity curve , resulting from the least-squares fit (13) of the family (12) to the fundamental diagram data.

  3. ARZ: The second-order ARZ model (6) that generalizes the least-squares fitted flow rate curve of the LWR model to a family of curves .

  4. GARZ: The second-order generalized ARZ model (8), whose generalized velocity function is obtained via a weighted least-squares fit (14) of the family (12) to the fundamental diagram data.

The fundamental diagram curves of the three traffic models are shown in Fig. 5, overlayed on top of the traffic state data that are actually observed in the test cases (gray dots). The top row of figures corresponds to the NGSIM data, and the bottom row represents the RTMC data. The LWR model is shown on the left, the ARZ model in the middle, and the GARZ model on the right.

Regarding the reproduction of real traffic data, it should be stressed that the models/predictors differ in the way they use data. In the model generation step, the Interpolation predictor uses no historic data; the LWR and ARZ model are based on the same least-squares fit; and the GARZ model employs more information from the fundamental diagram data due to the weighted least-squares fit. In turn, during the advance forward in time, the Interpolation predictor uses two pieces of information ( and ) at each boundary. In contrast, the LWR model uses only one piece of information () at one of the two boundaries (assuming the traffic states at the two boundaries are either both in free flow or both congested). Finally, the ARZ and GARZ model use a total of two pieces of information through the boundary conditions.

Figure 6. Model comparison on the NGSIM 5:15pm–5:30pm data set. In each panel, the time-evolutions of the model predictions (colored dashed curve) and measured data (solid gray curve) at the middle position of the study area are shown (top box: , bottom box: ). The four panels correspond to: Interpolation (top left, green), LWR (top right, blue), ARZ (bottom left, red), GARZ (bottom right, black).
Figure 7. Comparison of traffic models for the data sets NGSIM 4:00pm–4:15pm (at the top), NGSIM 5:00pm–5:15pm (in the middle), and NGSIM 5:15pm–5:30pm (at the bottom). In each panel, the time-evolution of spatially averaged errors (left top box), measurement traffic density data (left bottom box), and space-time average errors (right box) are shown in log-scale for Interpolation, LWR, ARZ, and GARZ.

5.4. Model Validation with the NGSIM Trajectory Data

As in [26], we consider a segment of 450 meters in length in the domain of the NSGIM vehicle trajectories. Data for three time intervals is available: 4:00pm–4:15pm, 5:00pm–5:15pm, and 5:15pm–5:30pm. However, because our model validation requires a traffic state be defined on the complete study segment, slightly shorter study time intervals must be chosen that guarantee that recorded vehicles are present everywhere on the road. We choose: 4:00:30pm–4:14:00pm, 5:00:30pm–5:13:30pm, and 5:15:30pm–5:28:00pm. The parameters of the traffic models, obtained by fitting to the fundamental diagram data provided with NSGIM, are given in the NGSIM row of Table 1. In line with the temporal resolution of the data, we generate density and velocity functions (real data and model predictions) in intervals of 0.1 seconds.

Figure 6 displays the time-evolution (5:15pm–5:30pm) of the traffic density and velocity that are predicted by the selected models, in comparison with the real evolution of these quantities, at the mid-point of the study area. Both the LWR and the ARZ model do reproduce the general trends present in the true density evolution, albeit with a systematic delay of 30–60 seconds (see [27] for a detailed discussion on this delay). One difference between LWR and ARZ is that the latter reproduces, modulo the delay, the velocity evolution better than the former. This is particularly visible in the predicted velocity in the plateau between 5:19:00pm and 5:22:30pm. In turn, the GARZ model captures the evolution of density and velocity significantly better than the other two traffic models. There is still a systematic delay in the model predictions, but this delay is very small. Finally, the Interpolation predictor captures the large–scale features of the real data as well. However, it also exhibits an oscillatory behavior of a larger frequency. This is due to the fact that the linear interpolations transmit temporal oscillations from both boundaries immediately to the observation site. In contrast, the traffic models pick up fewer boundary data and furthermore oscillations can turn into N-waves and thus reduce in magnitude as they move into the domain.

Data set (veh/h/lane) (km/h) (km/h)
NGSIM 247.38 23.41 0.16 39.49 82.01
RTMC 316.46 23.91 0.16 75.81 102.82
Table 1. Model parameters for the data-fitted models for the two data sets. Here, , , and are the free parameters of the equilibrium flow rate curve in the family (12). Moreover, and denote the boundaries of the empty road velocity for the GARZ model, as described in §3.2.

To quantify the predictive accuracy of the different models, we turn to the error metric (20). Figure 7 shows the model errors for the NGSIM data sets: 4:00pm–4:15pm (top), 5:00pm–5:15pm (middle), and 5:15pm–5:30pm (bottom). In each panel, the time-evolution of the spatially averaged error (21) is shown (top-left box), as well as the spatio-temporal average error (22) (top-right box). The models are: Interpolation (thick solid yellow), LWR (solid red), ARZ (dashed blue), and GARZ (thin solid black). In each bottom-left box, the time-evolution of the average traffic density is shown. The numerical values of the model errors (22) are given in rows 2–4 of Table 2.

Data set Interp. LWR ARZ GARZ
NGSIM 4:00–4:15 0.151 (+10%) 0.181 (+31%) 0.153 (+11%) 0.138
NGSIM 5:00–5:15 0.160 (+25%) 0.161 (+26%) 0.174 (+35%) 0.129
NGSIM 5:15–5:30 0.168 (+30%) 0.162 (+25%) 0.228 (+76%) 0.129
RTMC congested 0.203 (+14%) 0.228 (+24%) 0.208 (+13%) 0.184
RTMC non-cong. 0.108 (+63%) 0.081 (+26%) 0.067 (+4%) 0.064
Table 2. Spatio-temporal average errors of the traffic models and the Interpolation predictor for the NGSIM data sets (rows 2–4) and the RTMC data (rows 5–6), separated into congested and non-congested days. In each row, the parentheses denote the excess error relative to the best model (GARZ in all cases).

The time-evolution of the errors (21) confirm several of the above observations, such as the highly oscillatory nature of the Interpolation predictor, and the good accuracy of the GARZ model. Another particularly visible effect is the bad performance of the ARZ model for large densities: for 4:00–4:15, the ARZ model yields smaller errors than LWR; in contrast, for 5:15–5:30, the ARZ leads to larger errors than LWR. Moreover, the peaks in the ARZ errors coincide with local maxima in the traffic density. These observations confirm that: a) the ARZ model has the potential to improve upon the LWR model, and b) the non-uniform stagnation density of the ARZ model significantly affects its predictive accuracy. In contrast, the GARZ model inherits the advantages of the ARZ model for low densities, and furthermore remedies its shortcoming for high densities.

Regarding the performance of the Interpolation predictor; its accuracy is, except for the large oscillations, surprisingly good. Specifically, its average errors (22) are similar to those of the LWR and ARZ models. The following explanations can be provided to address this observation:

  1. The considered road segment is very short. Hence, the coherence between the boundaries and the inside of the domain is high, and the actual delays due to finite-speed information propagation are small. One can expect that on longer road segments, the finite-speed transport of information becomes more relevant, and therefore interpolation significantly loses in accuracy.

  2. As described in §5.3, the Interpolation predictor picks up twice (four times) as much data from the boundaries as the second (first) order models. One could therefore argue that the traffic models are as accurate as the interpolation, while utilizing less input data.

  3. The Interpolation predictor is on par with an “incomplete” model (LWR does not evolve velocities) and a “defective” model (ARZ is flawed for large densities). In turn, it does not achieve the accuracy of the GARZ model.

Figure 8. Model comparison on March 26, 2003 in the RTMC data set. In each panel, the time-evolutions (4:05pm–5:00pm) of the model predictions (colored dashed curve) and measured data (solid gray curve) at the middle sensor are shown (top box: , bottom box: ). The four panels correspond to: Interpolation (top left, green), LWR (top right, blue), ARZ (bottom left, red), GARZ (bottom right, black).
Figure 9. Comparison of models for the RTMC sensor data on 45 week days with congested traffic (top) and 29 days with non-congested traffic (bottom). The left boxes show the time-averaged error for each day (23), and the right boxes show the fully averaged errors (24) of congested and non-congested days. The errors for Interpolation, LWR, ARZ, and GARZ are depicted in log-scale.

5.5. Model Validation with the RTMC Sensor Data

In line with the test described in [26], we consider a 1.214 km long segment of highway on the I-35W, Minneapolis. Two traffic sensors are at the ends of the study segment, and one sensor is inside of the segment (roughly in the middle). The parameters of the data-fitted models, obtained from one-year data at the middle sensor, are given in the RTMC row of Table 1. On each weekday (Monday–Friday) between 01/01/2003 and 04/14/2003, we consider the onset of afternoon rush hour between 4pm and 5pm, and we divide the 74 days into 45 days with congested traffic (the space-time averaged traffic density exceeds 20 veh/km/lane) and 29 days with non-congested traffic (the remaining ones). As described in §5.1, the traffic models are run through an initialization phase 4:00pm–4:05pm, in which boundary data create a realistic state inside the segment. The actual validation is then conducted in 4:05pm–5:00pm.

Analogously as for the NGSIM data, we first consider the temporal evolution of the traffic states that the models predict and study the qualitative behavior of the models. We look at a day (03/26/2003) on which congestion builds up rapidly (between 4:35pm and 4:40pm) at the sensor position, so that any delays in the model behavior become visible. The results, shown in Fig. 8, confirm the observations made for the NGSIM data: a) LWR and ARZ predict similar densities, but ARZ does a better job at also capturing velocities correctly; b) both LWR and ARZ propagate information too slowly, resulting in the onset of congestion to be predicted 5 minutes too late; c) the GARZ model is not perfect (it still predicts the onset of congestion 2 minutes too late), but it captures the general trends in both variables quite nicely. This last point is particularly visible in the velocity profile in the congested state in 4:42pm–5:00pm.

One aspect that is different from the NGSIM test is that the Interpolation predictor performs visibly worse than the traffic models during the low density and high velocity period 4:05pm–4:35pm. This demonstrates the assertion made in §5.4, that simple interpolation performs much worse on longer road segments.

The average model errors obtained with the RTMC data are shown in Fig. 9. The top panel collects the 45 congested days, and the bottom panel contains the 29 non-congested days. In the left boxes the time-averaged errors (23) in 4:05pm–5:00pm are shown for each day. The right boxes show the resulting total errors (24), resulting from averaging over all study days. The numerical values, as well as the excess errors of the models relative to GARZ, are shown in rows 5–6 of Table 2.

The results of the low density days demonstrate quite clearly that a) traffic models yield noticeably better predictions than simple interpolation; and b) second-order models reproduce the real traffic behavior better than the first-order LWR model. The fact that the GARZ model does not differ much from the ARZ model results from the fact that for low densities, the two families of flow rate curves are very similar (see Fig. 5). In turn, the results of the high density days confirm that a) as expected, the quality of the ARZ model worsens relative to the other models; while b) the GARZ model does not suffer from the same amount of accuracy deterioration than ARZ. A somewhat unexpected observation is that for the high-density days, the Interpolation predictor does not perform worse than LWR and ARZ. The reason for this lies is the fact that LWR and ARZ propagate information too slowly, and thus capture features with a delay (see Fig. 8). In contrast, interpolation propagates information instantaneously—which is obviously unrealistic, but here happens to lead to less severe errors than the spurious delays in LWR and ARZ.

6. Inhomogeneous ARZ and GARZ Models

Thus far we have restricted our attention to homogeneous second-order models, in which each vehicle remains for all times attached to its specific velocity curve . However, it is plausible that in real traffic, drivers vary their empty road velocity over time, and that the overall traffic dynamics tend towards an equilibrium velocity curve . We therefore extend our model validations to the inhomogeneous ARZ model (7), denoted ARZ-, and the inhomogeneous GARZ model (9), denoted GARZ-. In both cases the equilibrium velocity curve is the LWR velocity curve.

As argued in §2.5.3, for a Cauchy problem the solutions of the inhomogeneous second-order models converge (as for fixed; or as for fixed) to solutions of the first-order LWR model. In the presence of boundary data, the same statement holds; however, the limits of the second-order model solutions are LWR solutions with different boundary conditions than the LWR solutions that we consider here. Consequently, the models ARZ- and GARZ- are not merely perturbations of the LWR model, but instead could reproduce the dynamics of real traffic better than LWR and better than homogeneous second-order models.

Figure 10. Model errors of the inhomogeneous ARZ and GARZ models, ARZ- (dashed blue) and GARZ- (solid black), as functions of the relaxation time . Shown are the time-averaged errors (23) for one day in the RTMC data set (top left), and spatio-temporally averaged errors (22) for the three NSGIM data sets: 4:00pm–4:15pm (top right), 5:00pm–5:15pm (bottom left), and 5:15pm–5:30pm (bottom right). Also shown are the errors with the homogeneous (i.e., ) ARZ model (blue square) and GARZ model (black circle), the error obtained with the LWR model (red square), and the choices of that yield the smallest model error (red star and red plus).

The inhomogeneous second-order models are based on the same data-fitted functions as their homogeneous counterparts. The only new parameter is the relaxation time . Since historic fundamental diagram data is commonly devoid of temporal information relevant for traffic dynamics, there is no canonical way to obtain the value of from historic data. Regarding realistic choices for , the only information found in the literature is that it cannot be much smaller than 3 seconds due to physical restrictions of the vehicle engines. Due to this absence of a good modeling principle for the value of , we conduct our model validation procedure for multiple instances of the inhomogeneous ARZ and GARZ models, where we let range from milliseconds to days. For each test, whichever choice of yields the smallest model error is then the optimal ARZ- (respectively optimal GARZ-) model.

The described -study is conducted for one day (January 8, 2003) in the RTMC data set, and for the three NGSIM data sets. For the RTMC data we evaluate the temporally averaged error (23), and for the NSGIM data the spatio-temporally averaged errors (22). The results of the study are shown in Fig. 10. In each test case, the ARZ- and GARZ- models are computed for many values of (blue and black functions, respectively). In addition, the (homogeneous) ARZ and GARZ models are computed (blue square and black circle). As expected, they agree with ARZ- and GARZ-, respectively, for . Moreover, the LWR model is computed (red square). As argued above, due to the presence of boundary data with , its results are different from ARZ- and GARZ- for . Finally, in each test case the particular is marked (red star and red plus) for which ARZ- and GARZ-, respectively, yield the smallest model errors. It is apparent that in all four cases shown in Fig. 10, the inhomogeneous models possess an optimal relaxation time .

NGSIM 4:00–4:15 0.153 (+13%) 0.141 (+4%) 0.138 (+2%) 0.135
NGSIM 5:00–5:15 0.174 (+43%) 0.137 (+13%) 0.129 (+6%) 0.122
NGSIM 5:15–5:30 0.228 (+87%) 0.165 (+35%) 0.129 (+6%) 0.122
RTMC congested 0.208 (+18%) 0.192 (+8%) 0.184 (+4%) 0.177
RTMC non-cong. 0.067 (+6%) 0.066 (+5%) 0.064 (+2%) 0.063
Table 3. Spatio-temporal average errors of the homogeneous ARZ and GARZ models, together with the inhomogeneous versions, ARZ- and GARZ-, with optimal choices of the relaxation time . Shown are the results for the three NGSIM data sets (rows 2–4), as well as the RTMC data, separated into congested and non-congested days (rows 5–6). In each row, the parentheses denote the excess error relative to the best model (GARZ- in all cases). For the RTMC data, the optimal