Generalization Second Order Macroscopic Traffic Models via Relative Velocity of the Congestion Propagation

Abstract

This paper presents a generalized second-order hydrodynamic traffic model. Its central piece is the expression for the relative velocity of the congestion (compression wave) propagation. We show that the well-known second-order models of Payne–Whitham, Aw–Rascal and Zhang are all special cases of the featured generalized model, and their properties are fully defined by how the relative velocity of the congestion is expressed. The proposed model is verified with traffic data from a segment of the Interstate 580 freeway in California, USA, collected by the California DOT’s Performance Measurement System (PeMS).

1. Introduction

In this paper, we introduce a new second-order hydrodynamic traffic model that generalizes the existing second-order models. The key feature of the new model is the notion of the relative velocity of the congestion (compression wave) propagation. This quantity can be obtained empirically from traffic detector measurements. Moreover, knowledge of the congestion relative velocity obviates the need for the fundamental diagram. We test the proposed model through numerical experiments using typical traffic detector data from the Performance Measurement System (PeMS) ([1]), on the I-580 freeway in California, USA.

An extensive development of gas dynamics (generalized conservation laws, stable differencing schemes) occurred in the 1950s, when the first macroscopic (hydrodynamic) models were developed. In those models, traffic flow was likened to a flow of “motivated” compressible liquid. For instance, the Lighthill–Whitham–Richards (LWR) model ([2,3]) describes the traffic flow through the law of conservation of vehicles. This model postulates the unique dependence of traffic flow on traffic density, which is called a fundamental diagram, and assumes a unique relationship between traffic speed $v(t,x)$ and density $\rho (t,x)$: $v(t,x)=V\left(\rho \right(t,x\left)\right)$. Here, $\rho (t,x)$ denotes vehicle density—the number of vehicles per unit length of the road at time t around the position x; and $v(t,x)$ denotes mean traffic speed at time t around the position x. The function $V\left(\rho \right)$ is non-increasing: $\frac{\partial V}{\partial \rho }\le 0$. We denote the flow (i.e., the number of vehicles passing a reference point in unit time) as $Q\left(\rho \right)=\rho V\left(\rho \right)$. The function is referred to as the fundamental diagram (although sometimes this name is reserved for $V\left(\rho \right)$). The vehicle conservation law in the LWR model is expressed in the differential form of a continuity equation with zero right-hand side:

(1)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v\right)}{\partial x}=\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho V\left(\rho \right)\right)}{\partial x}=\frac{\partial \rho }{\partial t}+\frac{Q\left(\rho \right)}{\partial x}=0\hfill \\ v(t,x)=V\left(\rho \right(t,x\left)\right)\hfill \end{array}\right.$

Later, the class of macroscopic models has expanded to systems of non-linear hyperbolic second-order PDEs, where the unique dependence of traffic density and speed is no longer assumed ([4,5,6,7,8,9,10,11,12]). These models differ in the way that they describe the dependency of traffic flow (or velocity) on density. Here, we will discuss some popular second-order models and show how they can be generalized using the relative velocity of the congestion.

We start with the Payne–Whitham model ([4,13]) with zero right-hand side to simplify the subsequent calculations:

(2)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=0\hfill \\ \frac{\partial v}{\partial t}+V\frac{\partial v}{\partial x}+\frac{c_0^2}{\partial \rho }\frac{\partial \rho }{\partial x}=0\hfill \end{array}\right.$

where $c_0<0$ is the so-called traffic ‘sound speed’, and generalize it for the arbitrary form of the relative velocity of the congestion, $c\left(\rho \right)$:

(3)$c{\left(\rho \right)}^2=\frac{\partial P\left(\rho \right)}{\partial \rho }\Rightarrow P\left(\rho \right)={\int }_0^{\rho }c{\left(\tilde{\rho }\right)}^{2}d\tilde{\rho }$

With $P\left(\rho \right)$ defined in (3), model (2) can be rewritten ([8]):

(4)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+\left(\frac{\partial Q}{\partial \rho }\right)\frac{\partial \rho }{\partial x}=\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}+\rho \frac{\partial v}{\partial x}=0\hfill \\ \frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\frac{c{\left(\rho \right)}^2}{\rho }\frac{\partial \rho }{\partial x}=\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\frac{1}{\rho }\frac{\partial P}{\partial \rho }\frac{\partial \rho }{\partial x}=\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\frac{1}{\rho }\frac{\partial P}{\partial x}=0\hfill \end{array}\right.$

which generalizes (2). The proposed model (4) can also be written in the conservative form:

(5)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v\right)}{\partial x}=0\hfill \\ \frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2+P\left(\rho \right))\right)}{\partial x}=0\hfill \end{array}\right.$

The work ([8]) also proposed the model (4) in anisotropic form, defined by system:

(6)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}+\rho \frac{\partial v}{\partial x}=0\hfill \\ \frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+c\left(\rho \right)\frac{\partial v}{\partial x}=\frac{\partial v}{\partial t}+\left(v+\rho \frac{\partial V\left(\rho \right)}{\rho }\right)\frac{\partial v}{\partial x}=0\hfill \end{array}\right.$

where the relative velocity of the congestion propagation is expressed through the derivative of the equilibrium speed: $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$

Note that system (6) is almost identical to the system proposed in ([7]) when its second equation is in a non-conservative form. The difference is in the form of $c\left(\rho \right)$. In ([7]), the smooth function $P\left(\rho \right)$ is used instead of the equilibrium speed $V\left(\rho \right)$:

(7)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}+\rho \frac{\partial v}{\partial x}=0\hfill \\ \frac{\partial v}{\partial t}+\left(v-\rho \frac{\partial P\left(\rho \right)}{\rho }\right)\frac{\partial v}{\partial x}=\frac{\partial v}{\partial t}+\left(v+\rho \frac{\partial V\left(\rho \right)}{\rho }\right)\frac{\partial v}{\partial x}=0\hfill \end{array}\right.$

From the physical standpoint, the use of the equilibrium speed $V\left(\rho \right)$ in (7) instead of an arbitrary smooth increasing function $P\left(\rho \right)\sim {\rho }^{\gamma },\gamma >0$ is more justified (see ([8])), as the function $V\left(\rho \right)$ is generally not smooth in points of transition between different traffic phases.

A system similar to (6), but with an additional relaxation term $\beta (\rho ,v)\left(V\right(\rho )-v)$ in the right-hand side, was introduced in ([10,11]):

(8)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}+\rho \frac{\partial v}{\partial x}=0\hfill \\ \frac{\partial v}{\partial t}+\left(v+\rho \frac{\partial V\left(\rho \right)}{\partial \rho }\right)\frac{\partial v}{\partial x}=\beta (\rho ,v)(V\left(\rho \right)-v)\hfill \end{array}\right.$

where the equilibrium speed is expressed as $V\left(\rho \right)=V_{max}{\left(1-{\left(\frac{\rho }{{\rho }_{max}}\right)}^{n_1}\right)}^{n_2}$ in ([10,11]); and as $V\left(\rho \right)=V_{max}\left(1-exp\left(-\frac{\lambda }{V_{max}}\left(\frac{1}{\rho }-\frac{1}{{\rho }_{max}}\right)\right)\right)$ in ([11]).

The rest of the paper is organized as follows. Section 2 has two parts: in Section 2.1 the relative velocity of congestion propagation is expressed using the model variables, traffic velocity v and density $\rho$; Section 2.2 proposes the generalized second-order model and shows, using the relative velocity of the congestion (compression wave) propagation, that the existing second-order models are special cases of the proposed one. Section 3 presents the discretized numerical method for model simulation. Section 4 discusses the simulation results, which support the proposition that different second-order models with the same expression for the relative velocity of the congestion behave similarly. Section 5 concludes the paper.

2. Generalized Model

2.1. Relative Velocity of the Congestion Propagation

First, we need to find a formal way to express the relative velocity of the congestion propagation, $c\left(\rho \right)$, or traffic pressure $P\left(\rho \right)={\int }_0^{\rho }c{\left(\tilde{\rho }\right)}^{2}d\tilde{\rho }$ using the model variables, traffic speed v and density $\rho$.

From the continuity equation, we know that:

$\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial \rho }\frac{\partial \rho }{\partial x}=0\Rightarrow \frac{1}{{\rho }_x}=-\frac{Q_{\rho }}{{\rho }_t}$

Using this, we express partial derivative $P_{\rho }=\frac{\partial P}{\partial \rho }$ via $P_x=\frac{\partial P}{\partial x}$:

(9)$P_{\rho }=\frac{\partial P}{\partial \rho }=\frac{\partial P}{\partial x}\frac{\partial x}{\partial \rho }=\frac{P_x}{{\rho }_x}=-\frac{Q_{\rho }}{{\rho }_t}{P}_x$

Then, we express $P_x=\frac{\partial P}{\partial x}$ from the second equation of (5):

$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2+P\left(\rho \right)\right)}{\partial x}=\frac{\partial Q}{\partial t}+\frac{\partial \left(Qv\right)}{\partial x}+\frac{\partial P}{\partial x}=0\Rightarrow$

$-P_x=-\frac{\partial P}{\partial x}=\frac{\partial Q}{\partial t}+\frac{\partial \left(Qv\right)}{\partial x}=\frac{\partial Q}{\partial t}+\frac{\partial \left(Qv\right)}{\partial Q}\frac{\partial Q}{\partial x}=\left(\frac{\partial \rho }{\partial t}=-\frac{\partial Q}{\partial x}\right)=\frac{\partial Q}{\partial t}-\frac{\partial \left(Qv\right)}{\partial Q}\frac{\partial \rho }{\partial t}$

and substitute $P_x$ into the Equation (9):

$\frac{\partial P}{\partial \rho }=-\frac{Q_p}{{\rho }_t}{P}_x=\frac{Q_p}{{\rho }_t}\left(\frac{\partial Q}{\partial t}-\frac{\partial \left(Qv\right)}{\partial Q}\frac{\partial \rho }{\partial t}\right)=\frac{\partial Q}{\partial \rho }\left(\frac{\partial Q}{\partial \rho }-\frac{\partial \left(Qv\right)}{\partial Q}\right)={\left(\frac{\partial Q}{\partial \rho }\right)}^2-\frac{\partial \left(Qv\right)}{\partial \rho }$

Taking ${\left(\frac{\partial Q}{\partial \rho }\right)}^2$ and $\frac{\partial \left(Qv\right)}{\partial \rho }$ from the known relations:

${\left(\frac{\partial Q}{\partial \rho }\right)}^2={\left(\frac{\partial \left(\rho v\right)}{\partial \rho }\right)}^2={\left(v+\rho \frac{\partial v}{\partial \rho }\right)}^2=v^2+2\rho v\frac{\partial v}{\partial \rho }+{\rho }^2{\left(\frac{\partial v}{\partial \rho }\right)}^2$

$\frac{\partial \left(Qv\right)}{\partial \rho }=\frac{\partial \left(\rho v^2\right)}{\partial \rho }=v^2+\rho \frac{\partial \left(v^2\right)}{\partial \rho }=v^2+2\rho v\frac{\partial v}{\partial \rho }$

and using them in the expression of $P_{\rho }=\frac{\partial P}{\partial \rho }$, we obtain the equation for $c{\left(\rho \right)}^2$:

(10)$c{\left(\rho \right)}^2=\frac{\partial P\left(\rho \right)}{\partial \rho }={\left(\frac{\partial Q}{\partial \rho }\right)}^2-\frac{\partial \left(Qv\right)}{\partial \rho }={\rho }^2{\left(\frac{\partial v}{\partial \rho }\right)}^2$

The sign: $c\left(\rho \right)=\pm \sqrt{\frac{\partial P\left(\rho \right)}{\partial \rho }}$ remains to be determined. For that, we substitute expression for $c{\left(\rho \right)}^2$ in the second equation of (4):

(11)$\begin{array}{c}\hfill \frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\frac{c{\left(\rho \right)}^2}{\rho }\frac{\partial \rho }{\partial x}=\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\rho \frac{\partial v}{\partial \rho }\left(\frac{\partial v}{\partial \rho }\frac{\partial \rho }{\partial x}\right)=\\ \hfill \frac{\partial v}{\partial t}+\left(v+\rho \frac{\partial v}{\partial \rho }\right)=\frac{\partial v}{\partial t}+\left(v+c\left(\rho \right)\right)\frac{\partial v}{\partial x}=0\end{array}$

and finally obtain the formula for the relative velocity of the congestion propagation:

(12)$c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$

2.2. Model Equations

The vehicle conservation law in the LWR model (1) is expressed by the continuity equation in the differential form with zero right-hand side. Accounting for variations in vehicle numbers, for example, due to lane changes or freeway entrances and exits, Equation (1) takes the form:

(13)$\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+\rho \frac{\partial v}{\partial x}+v\frac{\partial \rho }{\partial x}=f_0$

where $f_0$ is the number of vehicles entering (positive) or leaving (negative) in unit time. Equation (13) by itself is not sufficient for an adequate description of all the traffic flow phases ([5]). To correct this, we use the differential transformation of the conservation law ([14]). We multiply (13) by $\frac{\partial v}{\partial \rho }$:

$\frac{\partial v}{\partial \rho }\left(\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}\right)=\frac{\partial v}{\partial \rho }\frac{\partial \rho }{\partial t}+\frac{\partial v}{\partial \rho }\left(\rho \frac{\partial v}{\partial x}+v\frac{\partial \rho }{\partial x}\right)=\frac{\partial v}{\partial \rho }{f}_0$

and get:

(14)$\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\rho \frac{\partial v}{\partial \rho }\frac{\partial v}{\partial x}=\frac{\partial v}{\partial t}+\left(v+\rho \frac{\partial v}{\partial \rho }\right)\frac{\partial v}{\partial x}=\frac{\partial v}{\partial \rho }{f}_0$

Recall that $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$ is the relative velocity of the congestion propagation. Now, we can rewrite Equations (13) and (14) in their final form:

(15)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+\rho \frac{\partial v}{\partial x}+v\frac{\partial \rho }{\partial x}=f_0\hfill \\ \frac{\partial v}{\partial t}+\left(v+\rho \frac{\partial v}{\partial \rho }\right)\frac{\partial v}{\partial x}=\frac{\partial v}{\partial t}+\left(v+c\left(\rho \right)\right)\frac{\partial v}{\partial x}=\frac{\partial v}{\partial \rho }{f}_0\hfill \end{array}\right.$

Thus, we constructed a second-order macroscopic model that does not assume a unique dependency between traffic speed and density, as in the LWR. The right-hand side $f_0$ of the first equation of (15) accounts for the number of vehicles incoming (positive) or leaving (negative) in unit time. Looking at the right-hand side of the velocity equation in (15), we see the impact of these vehicles on traffic speed.

The number of left and right boundary conditions in (15) depends on the signs of the eigenvalues of the system: ${\lambda }_1=v,{\lambda }_2=v+\rho \frac{\partial v}{\partial \rho }$. Since ${\lambda }_1=v$ is always $\ge 0$, we need to know the sign of ${\lambda }_2=v+\rho \frac{\partial v}{\partial \rho }$. At a freeway entrance, there are two boundary conditions, if ${\lambda }_1>0$, ${\lambda }_2>0$; one, if ${\lambda }_1>0$, ${\lambda }_2\le 0$; or zero, if ${\lambda }_1=0,{\lambda }_2\le 0$. The opposite is true for a freeway exit: zero boundary conditions, if both ${\lambda }_1$, ${\lambda }_2\ge 0$; one, if ${\lambda }_1\ge 0,{\lambda }_2<0$. Therefore, we can use time-dependent traffic flow $Q\left(t\right)$ and speed $v\left(t\right)$ as boundary conditions. Besides the boundary conditions, we set the following initial conditions:

$\rho (x,0)={\rho }_0\left(x\right),v(x,0)=v_0\left(x\right)$

Several shortcomings of the Payne–Whitham model (2) were pinpointed in ([5]). In particular, it was mentioned that, with strong spatial disturbances in initial conditions, negative values of velocity and/or density might appear in the solution. According to the model, the traffic is influenced by the vehicles traveling behind it, which was not possible in the case of a one-lane road. These inconsistencies arise because, initially, the system (2) (or (4)) was used to model barotropic gas that is isotropic, and all directions of motion for barotropic gas are equally probable.

Mathematically, this relates to the fact that the relative velocity of the congestion propagation, $c\left(\rho \right)=\sqrt{\frac{\partial P\left(\rho \right)}{\partial \rho }}=\rho \frac{\partial v}{\partial \rho }$, is squared in the second equations of systems (4) and (5). Thus, these equations are not sensitive to the sign of the relative velocity and are agnostic of the congestion direction.

To resolve this issue, in system (11), we turned from $c{\left(\rho \right)}^2$ to $c\left(\rho \right)$ and obtained the resulting second-order macroscopic model (15) with zero right-hand side. It should also be noted that once we found the expression for the relative velocity of the congestion, $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$, the eigenvalues of the isotropic system (4), ${\lambda }_{1,2}=v\pm c\left(\rho \right)$, are transformed to ${\lambda }_1=v,{\lambda }_2=v+c\left(\rho \right)\le v$, in the system (15).

These results are stated in the following theorem.

Theorem 1.

Any macroscopic second-order hyperbolic model for the description of traffic flow dynamics can be formulated as

$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}+\rho \frac{\partial v}{\partial x}=f_0\hfill \\ \frac{\partial v}{\partial t}+\left(v+c\left(\rho \right)\right)\frac{\partial v}{\partial x}=\frac{\partial v}{\partial \rho }{f}_0+f_1\hfill \end{array}\right.$

by choosing a specific form of the relative velocity of the congestion propagation, $c\left(\rho \right)$. The term $f_0$ in the right-hand side of the model accounts for the number of vehicles incoming or leaving in unit time. The term $f_1$ plays the role of a relaxation term when necessary.

Remark 1.

Note that ([8]) was the first to propose an expression for the relative velocity of the congestion propagation in the form: $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$. This differs from (12) in that, instead of the observed traffic flow velocity, it uses the equilibrium velocity $V\left(\rho \right)$. Additionally, the expression for the relative velocity of the congestion propagation, $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$, in ([8]), was derived from the car-following model ([15]), while we obtained our expression $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$ from the system (5), with no additional restrictions on values of traffic speed.

The approach of ([8]) is justified, since we need to calculate the relative velocity of the congestion for each road segment and we have only two ways of doing this:

• First, we can define the relative velocity of the congestion propagation by using equilibrium velocity, as was carried out in ([8]), $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$. The function for equilibrium velocity $V\left(\rho \right)$ (fundamental diagram) can be empirically set using traffic detector data for a long period of time for each segment of the road.

• Second, we can do this without using equilibrium velocity or any form of a fundamental diagram. We approximate the value of the relative velocity of the congestion propagation using traffic density and velocity measured at the current time instant:

  (16)$c\left(\rho \right)=\frac{{\rho }_{in}+{\rho }_{out}}{2}\left(\frac{v_{out}-v_{in}}{{\rho }_{out}-{\rho }_{in}}\right)$

In (16), ${\rho }_{in},v_{in}$ and ${\rho }_{out},v_{out}$ denote the measured values of traffic densities and velocities at the current time instant at the entrance and exit of the chosen road segment. If ${\rho }_{in}={\rho }_{out}$, we can use the value of the relative velocity of congestion propagation at the previous time instant to avoid the division by zero.

As an example of our theorem application, we show how the model equations in the conservative form, proposed in ([7]),

(17)$\left\{\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v\right)}{\partial x}=0\hfill \\ \frac{\partial \left(\rho (v+P(\rho \left)\right)\right)}{\partial t}+\frac{\partial \left(\rho v(v+P(\rho \left)\right)\right)}{\partial x}=0\hfill \end{array}\right.$

can be derived from the system (5), using the expression for the relative velocity of the congestion in the form: $c\left(\rho \right)=-\rho \frac{\partial P\left(\rho \right)}{\partial \rho }$. The density conservation equations in both systems are identical, so we will focus on the momentum equations:

$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2+P\left(\rho \right))\right)}{\partial x}=\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2\right)}{\partial x}+\frac{\partial P\left(\rho \right)}{\partial \rho }\frac{\partial \rho }{\partial x}=0\iff$

$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2\right)}{\partial x}+c^2\left(\rho \right)\frac{\partial \rho }{\partial x}=0\iff$

$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2\right)}{\partial x}+{\rho }^2{\left(\frac{\partial P\left(\rho \right)}{\partial \rho }\right)}^2\frac{\partial \rho }{\partial x}=0\iff$

(18)$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2\right)}{\partial x}+{\rho }^2\frac{\partial P\left(\rho \right)}{\partial \rho }\frac{\partial P\left(\rho \right)}{\partial x}=0$

From the first equation of (5), using $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }=-\rho \frac{\partial P\left(\rho \right)}{\partial \rho }$, we obtain:

$\frac{\partial \rho }{\partial t}+\frac{\partial Q}{\partial x}=\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}+\rho \frac{\partial v}{\partial x}=\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}+\rho \frac{\partial v}{\partial \rho }\frac{\partial \rho }{\partial x}=0\Rightarrow$

$\frac{\partial \rho }{\partial t}+\left(v-\rho \frac{\partial P\left(\rho \right)}{\partial \rho }\right)\frac{\partial \rho }{\partial x}=0\Rightarrow$

$\frac{\partial P\left(\rho \right)}{\partial \rho }\left(\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial x}-\rho \frac{\partial P\left(\rho \right)}{\partial x}\right)=0\iff$

$\frac{\partial P\left(\rho \right)}{\partial t}+v\frac{\partial P\left(\rho \right)}{\partial x}-\rho \frac{\partial P\left(\rho \right)}{\partial \rho }\frac{\partial P\left(\rho \right)}{\partial x}=0\iff$

(19)${\rho }^2\frac{\partial P\left(\rho \right)}{\partial \rho }\frac{\partial P\left(\rho \right)}{\partial x}=\rho \left(\frac{\partial P\left(\rho \right)}{\partial t}+v\frac{\partial P\left(\rho \right)}{\partial x}\right)$

By substituting ${\rho }^2\frac{\partial P\left(\rho \right)}{\partial \rho }\frac{\partial P\left(\rho \right)}{\partial x}=\rho \left(\frac{\partial P\left(\rho \right)}{\partial t}+v\frac{\partial P\left(\rho \right)}{\partial x}\right)$ into the last of (18), we get:

$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2\right)}{\partial x}+\rho \left(\frac{\partial P\left(\rho \right)}{\partial t}+v\frac{\partial P\left(\rho \right)}{\partial x}\right)=0\iff$

$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho v^2\right)}{\partial x}+\rho \left(\frac{\partial P\left(\rho \right)}{\partial t}+v\frac{\partial P\left(\rho \right)}{\partial x}\right)+P\left(\rho \right)\left(\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v\right)}{\partial x}\right)=0\iff$

$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho P\left(\rho \right)\right)}{\partial t}+\frac{\partial \left(\rho v^2\right)}{\partial x}+\frac{\partial \left(\rho vP\left(\rho \right)\right)}{\partial x}=0\iff$

(20)$\frac{\partial \left(\rho \left(v+P\left(\rho \right)\right)\right)}{\partial t}+\frac{\partial \left(\rho v\left(v+P\left(\rho \right)\right)\right)}{\partial x}=0$

Therefore, we arrive at the momentum equation of (17), QED.

Reversing this logic, we can easily obtain system (5) from the system (17) using the same expression for the relative velocity of the congestion, $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }=-\rho \frac{\partial P\left(\rho \right)}{\partial \rho }$.

Thus, we see that the existing second-order hydrodynamic traffic models can be generalized using the relative velocity of the congestion propagation.

Additionally, it is important to know why the transition from macroscopic first-order models to second-order models is useful: adding the momentum equation into the system obviates the need for a fundamental diagram. Instead, the relative velocity of the congestion propagation can be obtained from traffic detector measurements. This enables the model to instantly adjust to real-time road conditions.

3. Computational Method

The second-order model (15) is of a hyperbolic type. A variety of finite-difference methods for solving such systems exist. By introducing vectors $W={\left[\rho ,v\right]}^T$ and $f={\left[f_0,\frac{\partial v}{\partial \rho }{f}_0+f_1\right]}^T$, the system (15) can be written in the vector form:

(21)$\frac{\partial W}{\partial t}+A\frac{\partial W}{\partial x}=\frac{\partial W}{\partial t}+\left({\mathsf{\Omega }}^{-1}\mathsf{\Lambda }\mathsf{\Omega }\right)\frac{\partial W}{\partial x}=f$

with the Jacobian:

(22)$A=\left[\begin{array}{cc}v& \rho \\ 0& v+\rho \frac{\partial v}{\partial \rho }\end{array}\right]=\left[\begin{array}{cc}v& \rho \\ 0& v+c\left(\rho \right)\end{array}\right]={\mathsf{\Omega }}^{-1}\mathsf{\Lambda }\mathsf{\Omega }$

In (22), $\mathsf{\Lambda }=\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]$ is a diagonal matrix of eigenvalues of the Jacobian A; $\mathsf{\Omega }=\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\end{array}\right]$ is the matrix of left eigenvectors of A; and ${\mathsf{\Omega }}^{-1}$ is the inverse of $\mathsf{\Omega }$.

Next, we write a finite-difference approximation on the grid shown in Figure 1 by making the values of the equation variables constant at the edges of the grid:

$\frac{W_m^{n+1}-W_m^n}{\Delta t}+A\left(\frac{W_{m+\frac{1}{2}}^{n+\frac{1}{2}}-W_{m-\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta x}\right)=\frac{W_m^{n+1}-W_m^n}{\Delta t}+\frac{{\left(AW\right)}_{m+\frac{1}{2}}^{n+\frac{1}{2}}-{\left(AW\right)}_{m-\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta x}=f_m^n$

Here, $m=1,\dots ,M$ is the grid node index along the x-axis, $n=1,\dots ,N$, is the grid node index along the t-axis, $\Delta t$ and $\Delta x$ are numerical integration steps bounded via the Courant condition: $\Delta t·\underset{m}{max}\left|{\left({\lambda }_{1,2}\right)}_m^n\right|<\Delta x$.

For system (15), the whole family of difference schemes can be expressed as follows (Kholodov and Tsybulin, 2018):

(23)$W_m^{n+1}=W_m^n-\frac{\Delta t}{\Delta x}\left(G_{m+\frac{1}{2}}^{n+\frac{1}{2}}-G_{m-\frac{1}{2}}^{n+\frac{1}{2}}\right)+f_m^n.$

The choice of the approximation method for $G_{m\pm \frac{1}{2}}^{n+\frac{1}{2}}={\left(AW\right)}_{m\pm \frac{1}{2}}^{n+\frac{1}{2}}$ in (23) is crucial for obtaining the scheme with the desired properties. When choosing a difference scheme, we presume that the solution for the model equations at every road segment is defined by the change in computed values in its boundary points. In this case, we can choose the conservative monotonic characteristic method of the first-order approximation ([16]) with the following expressions for $G_{m\pm \frac{1}{2}}^{n+\frac{1}{2}}$:

(24)$G_{m\pm \frac{1}{2}}^{n+\frac{1}{2}}=\frac{1}{2}{A}_{m\pm \frac{1}{2}}^n\left(W_m^n+W_{m\pm 1}^n\right)\pm \frac{1}{2}{\left({\mathsf{\Omega }}^{-1}\left|\mathsf{\Lambda }\right|\mathsf{\Omega }\right)}_{m\pm \frac{1}{2}}^n\left(W_m^n-W_{m\pm 1}^n\right)$

These could also be more complex expressions, which allow for the building of higher-order (hence, higher-accuracy) difference schemes for a given stencil. Variable values in the neighboring nodes $m\pm \frac{1}{2}$ in (24) can be computed using a simple linear interpolation without loss of accuracy.

$A_{m\pm \frac{1}{2}}^n=\left[\begin{array}{cc}\frac{v_m^n+v_{m\pm 1}^n}{2}& \frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2}\\ 0& \frac{v_m^n+v_{m\pm 1}^n}{2}+c\left(\frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2}\right)\end{array}\right],$

${\mathsf{\Lambda }}_{m\pm \frac{1}{2}}^n=\left[\begin{array}{cc}{\lambda }_{1,m\pm \frac{1}{2}}^n& 0\\ 0& {\lambda }_{2,m\pm \frac{1}{2}}^n\end{array}\right]=\left[\begin{array}{cc}\frac{v_m^n+v_{m\pm 1}^n}{2}& 0\\ 0& \frac{v_m^n+v_{m\pm 1}^n}{2}+c\left(\frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2}\right)\end{array}\right],$

${\mathsf{\Omega }}_{m\pm \frac{1}{2}}^n=\left[\begin{array}{c}{\omega }_{1,m\pm 1/2}^n\\ {\omega }_{2,m\pm 1/2}^n\end{array}\right]=\left[\begin{array}{cc}\frac{2}{{\rho }_m^n+{\rho }_{m\pm 1}^n}c\left(\frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2}\right)& -1\\ 0& 1\end{array}\right],$

(25)${\left({\mathsf{\Omega }}^{-1}\right)}_{m\pm \frac{1}{2}}^n=\left[\begin{array}{cc}\frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2c\left(\frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2}\right)}& \frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2c\left(\frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2}\right)}\\ 0& 1\end{array}\right].$

To obtain the characteristic form of compatibility equations, equivalent to the system (21) along the characteristics ${\lambda }_{1,2}=\frac{dx}{dt}$, we multiply the system (21) from the left by the eigenvector matrix $\mathsf{\Omega }$ (i.e., we compose 2 linearly independent combinations of the original equations) and bring it to the characteristic form:

$\mathsf{\Omega }\left(W_t+\left({\mathsf{\Omega }}^{-1}\mathsf{\Lambda }\mathsf{\Omega }\right)W_x\right)=\mathsf{\Omega }{W}_t+\mathsf{\Lambda }\mathsf{\Omega }{W}_x=\mathsf{\Omega }f.$

Then, in the scalar form, this yields:

(26)$\left({\omega }_i·W_t\right)+{\lambda }_i\left({\omega }_i·W_x\right)={\omega }_i·\left(\frac{\partial W}{\partial t}+{\lambda }_i\frac{\partial W}{\partial x}\right)=\left({\omega }_i·f\right),i=1,2$

Each of the two equations (26) is an ordinary differential equation along the characteristic ${\lambda }_{1,2}=\frac{dx}{dt}$.

We usually need to numerically integrate the compatibility conditions (26) at boundary points along the characteristics entering the integration area. For example, at the beginning of the freeway segment at boundary point $\left(t^{n+1},x_1\right)$ with negative eigenvalue ${\left({\lambda }_2\right)}_{\frac{3}{2}}^n=\frac{v_1^n+v_2^n}{2}+c\left(\frac{{\rho }_1^n+{\rho }_2^n}{2}\right)<0$, we use:

(27)$\left\{\begin{array}{c}{\left({\omega }_2\right)}_{\frac{3}{2}}^n·\left(\frac{W_1^{n+1}-W_1^n}{\Delta t}+{\left({\lambda }_2\right)}_{3/2}^n\frac{W_2^n-W_1^n}{\Delta x}\right)={\left({\omega }_2\right)}_{\frac{3}{2}}^n·f_1^n\hfill \\ v_1^{n+1}=v\left(t^{n+1},x_1\right)\hfill \end{array}\right.$

together with the boundary condition for traffic velocity $v\left(t,x_1\right)$.

At the end of the freeway segment, we do the same for the positive ${\left({\lambda }_2\right)}_{M-\frac{1}{2}}^n=\frac{v_M^n+v_{M-1}^n}{2}+c\left(\frac{{\rho }_M^n+{\rho }_{M-1}^n}{2}\right)>0$:

(28)${\left({\omega }_{1,2}\right)}_{M-\frac{1}{2}}^n·\left(\frac{W_M^{n+1}-W_M^n}{\Delta t}+{\left({\lambda }_{1,2}\right)}_{M-\frac{1}{2}}^n\frac{W_M^n-W_{M-1}^n}{x}\right)={\left({\omega }_{1,2}\right)}_{M-\frac{1}{2}}^n·f_1^n,$

for both compatibility equations, because ${\left({\lambda }_1\right)}_{m\pm \frac{1}{2}}^n=\frac{v_m^n+v_{m\pm 1}^n}{2}>0$ is always positive.

We should reiterate that the relative velocity of the congestion, $c\left(\rho \right)$ in (22)–(28) can be defined using equilibrium velocity $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$ as was done in ([8]). The function for equilibrium velocity $V\left(\rho \right)$ (fundamental diagram) can be empirically set using historic traffic detector data for each road segment. We propose the new approach without using equilibrium velocity or any form of a fundamental diagram $Q\left(\rho \right)$. We approximate the value of the relative velocity of the congestion propagation, $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$, using traffic density and velocity observed by traffic detection at the current time instant:

(29)$c{\left(\rho \right)}_{m\pm \frac{1}{2}}^n={\left(\rho \frac{\partial v}{\partial \rho }\right)}_{m\pm \frac{1}{2}}^n=\frac{{\rho }_m^n+{\rho }_{m\pm 1}^n}{2}\left(\frac{v_m^n-v_{m\pm 1}^n}{{\rho }_m^n-{\rho }_{m\pm 1}^n}\right).$

This approach enables our model (15) to instantly adjust to the changing road situation at runtime.

4. Numerical Results

To verify the proposed model, we carried out numerical experiments using traffic detector data for a segment of the I-580 freeway in California, USA, obtained from PeMS ([1]). We used two loop detectors, denoted $\#1$ and $\#2$, separated by $\sim 1$ km segment of the 4-lane freeway with no entrances or exits, as shown in Figure 2.

The data of the stationary traffic detectors $\#1$ and $\#2$ were used to construct the corresponding fundamental diagrams for equilibrium velocity $V\left(\rho \right)$. According to the three-phase traffic theory ([17]), the following traffic phases are distinguished:

• Free flow: $0\le \rho <{\rho }_1$

• Synchronized flow: ${\rho }_1\le \rho <{\rho }_2$

• Wide moving jam: ${\rho }_2\le \rho \le {\rho }_{max}$

For each phase, we define separate functions $V\left(\rho \right)$ and $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$, stitching them together at the transition points:

• Free flow: $\left\{\begin{array}{c}V\left(\rho \right)={\alpha }_2\rho +{\alpha }_1\hfill \\ c\left(\rho \right)={\alpha }_2\rho \hfill \end{array}\right.0\le \rho <{\rho }_1$

• Synchronized flow: $\left\{\begin{array}{cc}& V\left(\rho \right)={\beta }_2\rho +{\beta }_1+\frac{{\beta }_0}{\rho }\hfill \\ & c\left(\rho \right)={\beta }_2\rho -\frac{{\beta }_0}{\rho }\hfill \end{array}\right.,{\rho }_1\le \rho <{\rho }_2$

• Wide moving jam: $\left\{\begin{array}{cc}& V\left(\rho \right)=c_{*}\left(\frac{{\rho }_{max}}{\rho }-1\right)\hfill \\ & c\left(\rho \right)=-\frac{c_{*}{\rho }_{max}}{\rho }\hfill \end{array}\right.,{\rho }_2\le \rho \le {\rho }_{max}$

For the traffic detectors showcased in Figure 2, we have the following coefficients and parameters, aggregated over four lanes of I-580:

detectors $\#1$: ${\rho }_1=0.084$ veh./m, ${\alpha }_1=49.6$, ${\alpha }_2=-293.2$, ${\rho }_2=0.141$ veh./m, ${\beta }_0=2.49$, ${\beta }_1=-4.9$, ${\beta }_2=1.6$, $c_{*}=4.20$ m/s, ${\rho }_{max}=0.58$ veh./m;

detectors # 2: ${\rho }_1=0.076$ veh./m, ${\alpha }_1=50.2$, ${\alpha }_2=-295.7$, ${\rho }_2=0.165$ veh./m, ${\beta }_0=2.30$, ${\beta }_1=-0.4$, ${\beta }_2=-27.6$, $c_{*}=3.58$ m/s, ${\rho }_{max}=0.58$ veh./m;

For the separate lane of the I-580 freeway, the following parameters are different. For example, in the first lane:

detectors $\#1$: ${\rho }_1=0.023$ veh./m, ${\alpha }_1=52.5$, ${\alpha }_2=-1026.9$, ${\rho }_2=0.036$ veh./m, ${\beta }_0=0.46$, ${\beta }_1=22.3$, ${\beta }_2=-585.5$, $c_{*}=4.57$ m/s, ${\rho }_{max}=0.145$ veh./m;

detectors $\#2$: ${\rho }_1=0.021$ veh./m, ${\alpha }_1=53.2$, ${\alpha }_2=-1075.0$, ${\rho }_2=0.040$ veh./m, ${\beta }_0=0.58$, ${\beta }_1=10.3$, ${\beta }_2=-361.6$, $c_{*}=3.96$ m/s, ${\rho }_{max}=0.145$ veh./m;

The equilibrium velocity $V\left(\rho \right)$ — solid green line and the relative velocity of the congestion propagation $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$ — dashed blue line are shown as the functions of density for both detectors, together with the historic data for the one-year period in Figure 3 for the first lane and in Figure 4 aggregated over all four lanes.

Data (traffic flow and velocity, with 30-second temporal resolution) from detector $\#1$ were used as the initial and left boundary conditions, and data from the downstream detector $\#2$ were used for verification of the model output. We simulated a 24-h interval for a single weekday for all lanes of I-580. The results are shown in Figure 5 for the first lane and in Figure 6 aggregated over four lanes in different subplots: top left, traffic flows; top right, flow relative error in logarithmic scale; bottom left—traffic velocities; bottom right — velocity relative error in the logarithmic scale.

The simulation results obtained with a different expression for the relative velocity of the congestion propagation are shown in different colors: proposed macroscopic second-order model, where we approximate the value of relative velocity of the congestion, $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$, using instantaneous measurements of traffic density and speed (29), is shown in green; the same model (15), with equilibrium velocity $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$ as an empirical function of density obtained from the historic detector data, is shown in blue. The reference results (data of detector $\#2$) are depicted by a grey dashed line.

The simulation shows a good but not exact match with the results obtained using the same second-order macroscopic model (15) with a different form of the relative velocity of the congestion propagation, $c\left(\rho \right)$. We also see that both approaches have a good match with the reference results.

The root mean square error (RMSE) of velocity for the results of compared models and the detector $\#2$ measurements in Figure 5 is $0.25406$ m/s for the model (15) with $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$; and $0.31976$ m/s, if $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$. For the flow, RMSE is $0.00387$ veh./s, if $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$; and $0.00383$ veh./s, if $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$. For the results in Figure 6, velocity RMSE is $0.19646$ m/s, if $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$; and $0.2331$ m/s, if $c\left(\rho \right)=\rho \frac{\partial v}{\partial \rho }$. The flow RMSE is $0.01237$ veh./s if and $0.01301$ veh./s, respectively.

The difference in the results between these two approaches come from the time structure. In the first case, we have the continuous empirical function $c\left(\rho \right)=\rho \frac{\partial V\left(\rho \right)}{\partial \rho }$, with values dependent only on traffic density at each time instant. In the second case, we have a discrete function $c\left(\rho \right)=\frac{{\rho }_{in}+{\rho }_{out}}{2}\left(\frac{v_{out}-v_{in}}{{\rho }_{out}-{\rho }_{in}}\right)$, with values dependent both on traffic density and velocity. An instant error generally occurs in detector measurements, which is smoothed out in the first case, when we use the fundamental diagram.

5. Conclusions

In this paper, we generalized the second-order hydrodynamic macroscopic traffic models in various formulations [4,5,6,7,8,9,10,11]. Existing second-order macroscopic models describe traffic as non-linear systems of hyperbolic equations (for density and speed), which differ in the way they account for dependency between traffic flow (or velocity) and density. We have shown that all these second-order macroscopic models can be unified using the expression for the relative velocity of the congestion propagation.

To verify the proposed methodology, we conducted numerical experiments by simulating the segment of the I-580 freeway in California, USA, using data from PeMS ([1]). The simulations show that the results obtained for the same macroscopic model when using different forms of the relative velocity of the congestion propagation are slightly different. This suggests that the properties of every phenomenological model are defined by the form of the relative velocity of the congestion propagation—namely, the form of its dependency on density.

It is important to understand why the transition from macroscopic first-order models to second-order models is useful: adding the momentum equation into the system obviates the need for a fundamental diagram. Instead, the relative velocity of the congestion propagation can be obtained from traffic detector measurements. This enables the model to instantly adjust to real-time road conditions.

Author Contributions

Conceptualization, Y.K. and A.K.; methodology, Y.K. and A.A.; software, A.A.; validation, A.A. and V.K.; formal analysis, Y.K. and A.A.; investigation, Y.K. and A.A.; resources, Y.K.; data curation, A.A. and V.K.; writing—original draft preparation, Y.K.; writing—review and editing, Y.K. and A.K.; visualization, A.A. and V.K.; supervision, Y.K.; project administration, Y.K.; funding acquisition, Y.K. All authors have read and agreed to the published version of the manuscript.

Funding

The work was supported by grant RSCF 14-11-00877.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Contour bounding an integration cell.

Figure 2. Traffic detectors #1 and #2 chosen in the segment of the I-580 freeway in California, USA. Data from detector #1 were used as a left boundary condition for modeling, and data from the downstream detector #2 were used for the verification of modeling results.

Figure 3. Equilibrium speed V(ρ) (solid green line) and relative velocity of the congestion propagation cρ=ρ∂V(ρ)∂ρ (dashed blue line) as the functions of density for the first lane of I-580 together with the historic data for the one-year period. Left: detector #1, right: detector #2.

Figure 4. Equilibrium speed V(ρ) (solid green line) and relative velocity of the congestion propagation cρ=ρ∂V(ρ)∂ρ (dashed blue line) as functions of density aggregated over four lanes of I-580, together with the historic data for the one-year period. Left: detector #1, right: detector #2.

Figure 5. The comparison of calculated flows (top) and velocities (bottom) from the first lane of I-580 and observed values from detector #2 (dashed grey lines). On the right, flow (top) and velocity (bottom) relative errors in logarithmic scale. The simulation results obtained with a different expression for the relative velocity of the congestion propagation, cρ, are shown in different colors: green, model (with c(ρ)=ρ∂vρ∂ρ from (29); blue, the same model (15) with c(ρ)=ρ∂Vρ∂ρ as an empirical function of density.

Figure 6. The comparison of calculated flows (top) and velocities (bottom) aggregated over four lanes of I-580 and observed values from detector #2 (dashed grey lines). On the right, flow (top) and velocity (bottom) relative errors in logarithmic scale. The simulation results obtained with a different expression for the relative velocity of the congestion propagation, c(ρ), are shown in different colors: green, model (with c(ρ)=ρ∂vρ∂ρ from (29); blue, the same model (15) with c(ρ)=ρ∂Vρ∂ρ as an empirical function of density.