Comparative Model Accuracy of a DataFitted Generalized AwRascleZhang Model
Abstract.
The AwRascleZhang (ARZ) model can be interpreted as a generalization of the LighthillWhithamRichards (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 AwRascleZhang (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 datafitted GARZ model is compared to other traffic models by means of a threedetector 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, LighthillWhithamRichards, AwRascleZhang, generalized, second order, fundamental diagram, trajectory, sensor, data, validation2000 Mathematics Subject Classification:
35L65; 35Q91; 91B741. 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/gaskinetic (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., fluiddynamical, models. Among the (inviscid) macroscopic models one distinguishes between firstorder models based on scalar hyperbolic equations and secondorder models comprised of systems of hyperbolic equations. Specific examples for the latter are the PayneWhitham model [58, 72], twophase models [16], and the AwRascleZhang 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 wellsuited 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 online 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 firstorder 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 LighthillWhithamRichards model (2) in terms of HamiltonJacobi 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 timedependent data. Moreover, the main focus lies on secondorder models; and firstorder models are considered mainly for comparison purposes. The contributions of this paper are: (i) the design of a generalized AwRascleZhang model and the analysis of its mathematical properties; (ii) a systematic methodology to construct datafitted firstorder and secondorder traffic models, using historic fundamental diagram data; (iii) the validation of first and secondorder macroscopic traffic models via timedependent trajectory and sensor data, and the comparison of the predictive accuracy of different models; and (iv) the investigation of the optimal relaxation time in secondorder 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 AwRascleZhang model. We then introduce the generalized AwRascleZhang 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 datafitted first and secondorder macroscopic models. In §4 the numerical methods used to conduct the model validation and comparison are presented. Unlike studies of celltransmission 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 threedetector 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 datafitted secondorder 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
(1) 
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 LighthillWhithamRichards Model and Flow Rate Functions
The simplest macroscopic traffic model, the LighthillWhithamRichards (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
(2) 
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 NewellDaganzo flux [55, 19], in which is a piecewise linear function. While these different choices of functions lead to wellposed firstorder models, the secondorder 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 NewellDaganzo flux; and as a consequence, in this paper we consider flow rate functions (12) that resemble the shape of the NewellDaganzo flux, but that are strictly concave.
2.2. The AwRascleZhang 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 PayneWhitham model. The homogeneous AwRascleZhang (ARZ) model reads as
(3) 
where we call the hesitation function.^{1}^{1}1The 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 PayneWhitham model) yields the inhomogeneous ARZ model [33, 62]
(4) 
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
(5) 
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 datafitting 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 selfsustaining 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
(6) 
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 oneparameter 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 oneparameter 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
(7) 
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 AwRascleZhang (GARZ) model, which is a representative of the class of generic second order models (GSOM), proposed by Lebacque, Mammar, and HajSalem [48]. Specifically, the homogeneous ARZ model (6) generalizes to
(8) 
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 flowrate 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.
Proof.
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:
(9) 
Note that it would alternatively be conceivable to propose a relaxation in (9) of the form
(10) 
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).
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 datafitting 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:

given and , find , s.t. ;

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
where
The Jacobian of the flux function is
and its eigenvalues and associated eigenvectors are
and  
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 .
Proof.
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 1wave (i.e., a Laxshock or rarefaction associated to the first characteristic field), and to by a 2wave (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 1wave is a shock wave (moving with speed , given by the RankineHugoniot 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 1wave and the 2wave 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 2contact 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 firstorder equation is only warranted if a subcharacteristic 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 .
In contrast, for initialboundaryvalue 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 2wave) moving with speed , and to the interior state via a shock (a 1wave) 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. DataFitted 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 longterm 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. DataFitting for the LWR and ARZ Models
The firstorder 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 threeparameter family of smooth and strictly concave flow rate curves is selected as
(12) 
where
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 threeparameter family of flow rate curves, the one is selected that is the closest to the data points in a leastsquares sense, i.e. we solve
(13) 
In the right panel of Fig. 2, the resulting leastsquares 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 oneparameter 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 datafitted 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.
3.2. DataFitting 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 leastsquares fit (13) to a weighted leastsquares fit, as follows.
Given a weight parameter , we consider the minimization problem
(14) 
where
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 leastsquares problem (14) generates a oneparameter 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 nonintersecting, i.e.,
(15) 
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.
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 reparameterize 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 reparameterization, 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 semiimplicit 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
(16)  
(17) 
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 twoparameter 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:

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

The new density states are computed via (16).
In the special case of the ARZ model (5), the update (17) is given by the explicit formula
(19) 
It should further be remarked that the semiimplicit 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 threedetector 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 leastsquares 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 spaceandtimedependent error measure as
(20) 
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 .
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 fourstep 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 spatiotemporal averages are considered
(21)  
(22) 
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
(23)  
(24) 
5.3. List of Models
We compare the following four models in terms of their predictive accuracy of the real data:

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.

ARZ: The secondorder ARZ model (6) that generalizes the leastsquares fitted flow rate curve of the LWR model to a family of curves .
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 leastsquares fit; and the GARZ model employs more information from the fundamental diagram data due to the weighted leastsquares 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.
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 timeevolution (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 midpoint 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 Nwaves 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 
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 timeevolution of the spatially averaged error (21) is shown (topleft box), as well as the spatiotemporal average error (22) (topright box). The models are: Interpolation (thick solid yellow), LWR (solid red), ARZ (dashed blue), and GARZ (thin solid black). In each bottomleft box, the timeevolution 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  noncong.  0.108  (+63%)  0.081  (+26%)  0.067  (+4%)  0.064 
The timeevolution 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 nonuniform 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:

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 finitespeed information propagation are small. One can expect that on longer road segments, the finitespeed transport of information becomes more relevant, and therefore interpolation significantly loses in accuracy.

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.

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.
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 I35W, 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 datafitted models, obtained from oneyear 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 spacetime averaged traffic density exceeds 20 veh/km/lane) and 29 days with noncongested 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 noncongested days. In the left boxes the timeaveraged 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) secondorder models reproduce the real traffic behavior better than the firstorder 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 highdensity 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 secondorder 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 secondorder models converge (as for fixed; or as for fixed) to solutions of the firstorder LWR model. In the presence of boundary data, the same statement holds; however, the limits of the secondorder 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 secondorder models.
The inhomogeneous secondorder models are based on the same datafitted 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 spatiotemporally 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 .
Data set  ARZ  ARZ  GARZ  GARZ  

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  noncong.  0.067  (+6%)  0.066  (+5%)  0.064  (+2%)  0.063 