A variational formulation for higher order macroscopic traffic flow models: Numerical investigation

Transportation Research Part B 70 (2014) 112–133

Contents lists available at ScienceDirect

Transportation Research Part B journal homepage: www.elsevier.com/locate/trb

A variational formulation for higher order macroscopic traffic flow models: Numerical investigation Guillaume Costeseque a,b,⇑, Jean-Patrick Lebacque b a b

Université Paris-Est, Ecole des Ponts ParisTech, CERMICS, 6 & 8 avenue Blaise Pascal, Cité Descartes, Champs sur Marne, 77455 Marne la Vallée Cedex 2, France Université Paris-Est, IFSTTAR, GRETTIA, 14-20 Boulevard Newton, Cité Descartes, Champs sur Marne, 77447 Marne la Vallée Cedex 2, France

a r t i c l e

i n f o

Article history: Received 15 March 2014 Received in revised form 8 July 2014 Accepted 25 August 2014

Keywords: Traffic flow Hamilton–Jacobi equation Lax–Hopf algorithm Lagrangian data

a b s t r a c t This paper deals with numerical methods providing semi-analytic solutions to a wide class of macroscopic traffic flow models for piecewise affine initial and boundary conditions. In a very recent paper, a variational principle has been proved for models of the Generic Second Order Modeling (GSOM) family, yielding an adequate framework for effective numerical methods. Any model of the GSOM family can be recast into its Lagrangian form as a Hamilton–Jacobi equation (HJ) for which the solution is interpreted as the position of vehicles. This solution can be computed thanks to Lax–Hopf like formulas and a generalization of the inf-morphism property. The efficiency of this computational method is illustrated through a numerical example and finally a discussion about future developments is provided. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction 1.1. General background In order to get a realistic estimation of the real-time traffic states on networks, traffic operators and managers need macroscopic traffic flow models. These models must be simple, robust, allowing to get solutions at a low computational cost. The main macroscopic models are based on conservation laws or hyperbolic systems (see Hoogendoorn and Bovy (2001); Kuhne and Michalopoulos (1997) for traffic aspects and Garavello and Piccoli (2006) for mathematical aspects). The seminal LWR model (for Lighthill–Whitham and Richards) was proposed in Lighthill and Whitham (1955), Richards (1956) as a single conservation law with unknown the vehicle density. This model based on a first order Partial Differential Equation (PDE) is very simple and robust but it fails to recapture some empirical features of traffic. In particular, it does not allow to take into account non-equilibrium traffic states mainly in congested situation. More sophisticated models referred to as higher order models were developed to encompass kinematic constraints of real vehicles or also the wide variety of driver behaviors, even at the macroscopic level. In this paper we deal with models of the Generic Second Order Modeling (GSOM) family. Even if these models are more complicated to deal with, they permit to reproduce traffic instabilities (such as the so-called stopand-go waves, the hysteresis phenomenon or capacity drop) which move at the traffic speed and differ from kinematic waves (Zhang, 2002) (see also Lebacque et al. (2007) and references therein). Before the wide propagation of internet handsets, traffic monitoring has mainly been built on dedicated infrastructure which imply quite important installation and maintenance costs. Traffic flow monitoring and management has been deeply ⇑ Corresponding author at: Université Paris-Est, Ecole des Ponts ParisTech, CERMICS, 6 & 8 avenue Blaise Pascal, Cité Descartes, Champs sur Marne, 77455 Marne la Vallée Cedex 2, France. E-mail address: [email protected] (G. Costeseque). http://dx.doi.org/10.1016/j.trb.2014.08.012 0191-2615/Ó 2014 Elsevier Ltd. All rights reserved.

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

113

modified with the development of new technologies in mobile sensing aiming to provide a quite important quantity of floating car data. Traffic flow models are needed to be well suited such that managers could use both Eulerian and Lagrangian data for improving traffic state estimation. The term Eulerian refers to ‘‘classical’’ fixed equipment giving records of occupancy or flow of vehicles on a freeway section. This kind of measurements come from e.g. fixed inductive loop detectors, Radio Frequency Identification (RFID) transponders, radars or video cameras. By opposite, the term Lagrangian is used to characterize data coming from sensors which move within the measured field of interest. Lagrangian data are provided by on board mobile sensors such as Global Positioning Systems (GPS) or GPS-enabled smartphones, without being exhaustive. While the idea of monitoring traffic using mobile sensors appeared less than ten years ago with the popularization of the mobile internet and cellular devices, there exists a fast growing literature about how to integrate Lagrangian data into classical macroscopic traffic flow models. The process of incorporating Eulerian and Lagrangian data into a mathematical model to improve the modeling is called data estimation or equivalently inverse modeling. According to the major UC Berkeley field experiment named Mobile Century and then Mobile Millennium investigating Lagrangian sensing with GPS-enabled smartphones, it has been shown that even a 2% to 5% penetration rate of probe vehicles into the driver population provides sufficient and accurate data for estimating traffic velocity or density on highways (Herrera et al., 2010; Herrera and Bayen, 2010; Work et al., 2009). Nevertheless, it has been demonstrated in Piccoli et al. (2012) that the quality of estimation for higher-order traffic quantities including vehicle acceleration/deceleration, emission and fuel consumption rates is dramatically affected when the penetration rate of probe vehicles or the sampling frequency of the current mobile sensors decrease. However on board devices propose a real breakthrough in traffic monitoring by providing a very cheap and efficient way to collect traffic data. 1.2. Motivation In order to improve traffic states estimation from Lagrangian data, we propose to deal with macroscopic traffic flow models of the GSOM family. As these models combine the simplicity of the LWR model with the dynamics of driver specific attributes, we are able to recapture more specific phenomenon with a higher accuracy. While methods of data assimilation have been only developed for first order models up to now, this work presents a new algorithm to reconstruct traffic states from both Eulerian and Lagrangian data. We take advantage of a very recent article (Lebacque and Khoshyaran, 2013) in which a variational principle has been proved for models of the GSOM family. 1.3. Organization of the paper The rest of this paper is structured as follows. Section 2 presents more in detail the GSOM models and sheds a specific light on the LWR model which is widely used in traffic engineering. The variational principles for the GSOM models including LWR model are briefly recalled in Section 3. Section 4 is devoted to the presentation of the main elements of our computational method. Finally, Section 4.4 proposes a numerical example of our method. 2. GSOM traffic flow models 2.1. Formulation of GSOM models In Lebacque et al. (2005, 2007), the authors introduce a general class of macroscopic traffic flow models called the Generic Second Order Models (GSOM) family. Any model of the GSOM family can be stated in conservation form as follows

8 Conservation of vehicles; > < @ t q þ @ x ðqv Þ ¼ 0 @ t ðqIÞ þ @ x ðqv IÞ ¼ quðIÞ Dynamics of the driver attribute I; > : v ¼ Iðq; IÞ Fundamental diagram;

ð2:1Þ

where q stands for the density of vehicles, v for the flow speed (equal to the mean spatial velocity of vehicles), x and t for position and time. The variable I is a specific driver attribute which can represent for example the driver aggressiveness, the driver destination or the vehicle class. The flow-density fundamental diagram (FD) is defined by

F : ðq; IÞ # qIðq; IÞ: Notice moreover that it was shown in Lebacque et al. (2007) that the notions of Supply and Demand functions defined in Lebacque (1996) for the classical LWR model could be extended to the GSOM family. The GSOM models admit two kinds of waves:  Kinematic waves or 1-waves as in the seminal LWR model: a wave propagates density variations at speed m1 ¼ @ q Fðq; IÞ while the driver attribute I is continuous across such a wave.  Contact discontinuities or 2-waves: a wave propagates variations of driver attribute I at speed m2 ¼ Iðq; IÞ while the flow speed is constant across such a wave.

114

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

2.2. Examples The GSOM family recovers a wide range of existing models:  The LWR model (Lighthill and Whitham, 1955; Richards, 1956) itself is simply a GSOM model with no specific driver attribute, expressed as follows



  

 



@ t q þ @ x ðqv Þ ¼ 0 Conservation of vehicles; v ¼ Iðq; xÞ Fundamental diagram:

ð2:2Þ

The fundamental diagram (FD) for the LWR model states that traffic flow is always at an equilibrium state. It is commonly assumed that the flow is an increasing function of density between zero (corresponding to an empty section) and a critical density and then the flow decreases until the maximal density (corresponding to a bumper-to-bumper situation). However the FD shape is always a subject of debates (see for instance Fan and Seibold (2012)) and there exists a wide variety in the literature encompassing concave and triangular flow functions (see Fig. 1 and also Chapter 3 of Garavello and Piccoli (2006) for additional examples). The LWR model with bounded acceleration proposed in Lebacque (2002, 2003, 2007) is also a GSOM model in which the propagated driver attribute is simply the speed of vehicles. The ARZ model (standing for Aw and Rascle (2000) and Zhang (2002)) for which the driver attribute is taken as I ¼ v  V e ðqÞ that is Iðq; IÞ ¼ I þ V e ðqÞ. The Generalized ARZ model proposed in Fan et al. (2013) that can be also seen as a particular case of the model described in Zhang et al. (2009). These models introduce an interaction mechanism between two different FDs for distinguish equilibrium and non-equilibrium states. Multi-commodity models (multi-class, multi-lanes) of Jin and Zhang (2004), Bagnerini and Rascle (2003) or Herty et al. (2008). It encompasses also the model of Klar et al. (2003). The Colombo 1-phase model deduced in Lebacque et al. (2007) from the 2-phase model of Colombo (2003). In this case, the driver attribute I is a scalar which is non-trivial in congested situation. In fluid area, the model follows the classical LWR model. The stochastic GSOM model of Khoshyaran and Lebacque (2009). The driver attribute I is a random variable depending on the vehicle index N and on the random event x such that I ¼ IðN; t; xÞ. The random perturbations do not affect the vehicle dynamics but affect the driver perception and its behavior. The interested reader is referred to Lebacque and Khoshyaran (2013) and references therein for more details on examples.

3. Variational principles in traffic flow modeling In traffic flow literature, the variational formulation was first conjectured by Newell in Newell (1993) for the LWR firstorder traffic flow model. It was then properly established and generalized by Daganzo in Daganzo (2005a), Daganzo (2005b), Daganzo and Menendez (2005), Daganzo (2006). In this section, we first present the variational formulation of the GSOM family lead by the ideas developed for the LWR model. Then the numerical methods for solving these variational formulations are introduced. 3.1. Variational formulation of the GSOM family Guided by the ideas developed for the LWR model, variational formulations have been recently developed and proved for models of the GSOM family in both Eulerian (Li and Zhang, 2013) and Lagrangian frameworks (Lebacque and Khoshyaran, 2013).

Fig. 1. Illustrations of some flow functions F for the LWR model: Greenshields (left), triangular (center) and exponential (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

115

3.1.1. Lagrangian setting In Lebacque and Khoshyaran (2013), the authors prove the existence of a variational principle for the GSOM models family expressed under its Lagrangian form, introducing r :¼ 1=q the spacing and N the vehicle label

8 Conservation of vehicles; > < @t r þ @N v ¼ 0 @ t I ¼ uðN; I; tÞ Dynamics of the driver attribute I; > : v ¼ Vðr; IÞ ¼ Ið1=r; IÞ Fundamental diagram:

ð3:3Þ

Remark 3.1 (Vehicle labeling). Note that vehicles are labeled in the classical traffic sense, according to their passing time through a reference spatial point. Thus, the label axis is in the opposite sense that the spatial axis. Rt Considering the position X ðN; tÞ ¼ 1 v ðN; sÞds, we obviously have that

v ¼ @t X

and r ¼ @ N X:

ð3:4Þ

Hence the system (3.3) could be written as a Hamilton–Jacobi equation satisfied by X

@ t X  WðN; @ N X ; tÞ ¼ 0;

ð3:5Þ

where W denotes the speed-spacing fundamental diagram of vehicle N at time t. It is defined such that

WðN; r; tÞ :¼ Vðr; IðN; tÞÞ ¼ Ið1=r; IðN; tÞÞ; where the driver attribute IðN; tÞ solves the following nonlinear first order ODE

  @ t IðN; tÞ ¼ uðN; I; tÞ;   IðN; 0Þ ¼ i ðNÞ; for any N: 0 By classical results on optimal control problems i.e. dynamic programming on Hamilton–Jacobi–Bellman equations, we can check that the solution X of (3.5) satisfies, for any point ðN T ; TÞ 2 R  ½0; þ1Þ, the following relation

RT XðNT ; TÞ ¼ min t0 MðN; u; tÞdt þ nðN0 ; t0 Þ; u;ðN 0 ;t 0 Þ   N_ ¼ u    Nðt0 Þ ¼ N0 ; NðTÞ ¼ NT   ðN ; t Þ 2 J 0

ð3:6Þ

0

where M is the Legendre–Fenchel transform of W according to its second variable (see Fig. 2 for an illustration of the Legendre–Fenchel transform), defined such that

MðN; u; tÞ ¼ sup fWðN; r; tÞ  urg:

ð3:7Þ

r2R

Optimality conditions. Notice that the variational formulation of GSOM models (3.6) is only available if and only if W is Lipschitz continuous w.r.t. to its first and third variables, and if r # WðN; r; tÞ is strictly concave, such that the Legendre–Fenchel transform of M gives back the function W. Moreover, we need the initial-boundary condition n to be Lipschitz continuous. In the previous formulation (3.6), J is the locus of initial and boundary conditions and nðN 0 ; t0 Þ is the initial data of the position of vehicle N 0 at time t 0 for any ðN 0 ; t0 Þ 2 J . The optimal trajectories are given by the characteristics (Eulerian and Lagrangian characteristics are equal) described by the following coupled ODEs

(

N_ ¼ @ r WðN; r; tÞ; r_ ¼ @ N WðN; r; tÞ;

ð3:8Þ

Fig. 2. Legendre transform M of function W.

116

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

where u ¼ N_ is the command of the optimal control problem (3.6). The system (3.8) is a simple system of ODEs in the ðN; rÞ plane. 3.1.2. Eulerian setting In Li and Zhang (2013), the authors prove the existence of a variational formulation of a class of models from the GSOM family. Unlike the paper of Lebacque and Khoshyaran (2013), these models are expressed from the Eulerian viewpoint as a system of two conservation laws



@ t q þ @ x f ðq; sÞ ¼ 0 Conservation of vehicles; @ t s þ @ x gðq; sÞ ¼ 0

Dynamics around the equilibrium:

ð3:9Þ

The first conservation law is obviously satisfied by the vehicle density q and the corresponding flux f is the product of density times speed f ðq; v Þ ¼ qv . The second conservation law is satisfied by a variable s (a non-equilibrium measure) which can be replaced by s ¼ qI yielding the GSOM family without source terms uðIÞ ¼ 0. The corresponding flux g needs to satisfy some specific conditions to ensure that the wave speed is less or equal to traffic speed. The interested reader is referred to Li and Zhang (2013) for more details. The paper (Li and Zhang, 2013) highlights that taking uðIÞ ¼ 0 is equivalent to say that the driver attribute I is invariant along vehicle trajectories which seems to be corroborated for instance by the work of Duret et al. (2008) on the NGSIM I-80 trajectories data-set (in congested situation). For the models (3.9) where uðIÞ ¼ 0, the authors show in Li and Zhang (2013) a variational formulation based on Lax–Hopf like formulas for both N q and N s defined as the cumulative quantity of respectively q and s. Optimality conditions. The Lax–Hopf like formulas are ensured through the following assumptions: we need flux f to be strictly convex in its first variable and flux g to be convex in its second variable. We also need the initial-boundary data to be Lipschitz continuous. When taking into account a non trivial relaxation term (uðIÞ – 0), the authors assume that the problem reduces to solving a LWR model in large time because the relaxation term induces an exponential decay in time of the difference between (3.9) and a regular LWR model (2.2). Remark 3.2 (Variational formulations of the LWR model). The variational theory for the LWR relies on the three dimensional representation of traffic flow or the so-called Moskowitz function (Moskowitz and Newan, 1963; Mazaré et al., 2011) yielding Hamilton–Jacobi formulations in Eulerian and Lagrangian coordinates (see also Laval and Leclercq (2013) and references therein). Under suitable assumptions, the solutions to the HJ equations could be determined thanks to a variational formulation, known as Lax–Hopf formula in the mathematics literature (Evans, 2009). There exist different methods to obtain the variational formulation of these problems, including calculus of variations, dynamic programming (Daganzo, 2005a,b) and viability theory (Claudel and Bayen, 2010a; Han et al., 2012). Notice that from a mathematical point of view it was rigorously established in Han et al. (2012) that the (viscosity) solutions of HJ equations in Eulerian and Lagrangian coordinates are strictly equivalent. 3.2. Computational methods for GSOM models 3.2.1. Lax–Hopf algorithms in LWR case Up to our best knowledge, existing numerical procedures for data assimilation use the seminal paper (Claudel and Bayen, 2010a) which proposes a semi-explicit form of the solution to Hamilton–Jacobi equations with concave flow-density fundamental diagrams and with any piecewise affine (PWA) initial and boundary conditions. They introduce a generalized Lax– Hopf formula and the inf-morphism property (Aubin et al., 2008; Aubin et al., 2011) to compute the solution by taking the infimum of all solutions associated with simpler partial problems. It follows computational methods that are commonly referenced as ‘‘grid free’’ schemes (Claudel and Bayen, 2010b; Mazaré et al., 2011). Unlike dynamic programming methods (Daganzo, 2005a), Lax–Hopf algorithms have been proved to be exact in general cases of concave flow-density fundamental diagrams (i.e. for convex Hamiltonians) (Claudel and Bayen, 2010a). Some elements of comparison are given in Mazaré et al. (2011) between dynamic programming, Lax–Hopf algorithm and more classical methods such as Godunov scheme (Daganzo, 1995; Lebacque, 1996) and wave tracking algorithm (Henn, 2005; Holden and Risebro, 2011; Wong and Wong, 2002). The Lax–Hopf algorithm has been extended to the Lagrangian version of the LWR problem in Han et al. (2012) for a triangular flow-density FD and for piecewise affine initial and boundary conditions. 3.2.2. Review of computational methods for GSOM models There already exists some works on computational methods for models of the GSOM family (Mammar et al., 2009; Li and Zhang, 2013; Qiu et al., 2013). However these existing methods are developed in the Eulerian framework which does not seem to be the best one to deal with Lagrangian data. Moreover the paper (Qiu et al., 2013) deals with the LWR model with bounded acceleration which is a very specific GSOM model. The algorithm is very similar to the ones developed in Claudel and Bayen (2010b), Mazaré et al. (2011) but it is not applicable to general models of the GSOM family.

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

117

In Li and Zhang (2013) the authors present a numerical scheme that is based on a discretization of Eulerian time and space domain. The idea of the scheme (close to dynamic programming) is then to update the cumulative quantities

Nq :¼

Z

þ1

qðy; tÞdy and Ns :¼

x

Z

þ1

sðy; tÞdy;

x

on the vertexes of a mesh (whatever its shape if this mesh is dense enough) following the optimal paths in an iterative way because of the coupling of the equations to solve. The values of q (resp. s) are deduced from a first order approximation of the derivative of N q (resp. N s ). Notice that this computational method is not exact in general. 3.2.3. A finite difference scheme A numerical scheme based on a uniform grid of the Lagrangian domain with steps Dt and Dn, was presented in Lebacque and Khoshyaran (2013). We define for any n; t P 0

X tn :¼ X ðnDn; t DtÞ: As usual with finite difference scheme, we need to introduce a Courant–Friedrichs–Lewy (CFL) condition to ensure the scheme to be monotone and convergent

Dn P supj@ r WðN; r; tÞj: Dt N;r;t Then we consider the classical first order finite difference scheme as follows

  X tn1  X tn t : X tþ1 ¼ X þ D t W n D n; ; t D t n n Dn

ð3:10Þ

By construction the above scheme can be interpreted as the seminal Godunov (finite volume) scheme (3.11) for the Lagrangian formulation of the GSOM model (see Costeseque et al. (2014))

rtþ1 ¼ r tn þ n

 Dt  t W n1  W tn ; Dn

ð3:11Þ

where we have defined the discrete spacing and respectively the numerical speed as follows

(

rtn :¼ W tn

X tn1 X tn Dn



;

:¼ W nDn; rtn ; t Dt :

The upstream and downstream boundary conditions for the finite difference scheme (3.10) are fully described in Lebacque and Khoshyaran (2013). They match the Bardos–LeRoux–Nédélec or Dubois–LeFloch boundary conditions and it was already shown that such conditions are equivalent to prescribe supply and demand conditions (see Lebacque and Khoshyaran (2013) and references therein). There exists many numerical methods for dealing with hyperbolic systems of conservation laws. Our scope is not to make a review of these techniques because we are mainly focused on the Lagrangian reformulation of such systems with application to traffic. Among the very few works dealing with such an approach, we can quote (Greenberg, 2002 and Zhang et al., 2012) which are focused on the numerical resolution of the Aw-Rascle model (Aw and Rascle, 2000). The former deals with a downwind monotone scheme while the latter introduces a semi-discrete scheme. Let us also mention the semi-implicit finite volume scheme presented in Fan et al. (2013) for a generalized Aw and Rascle (2000), Zhang (2002) model. Up to our best knowledge, there is no other existing work on computational methods for the generic expression of GSOM models recast into the Lagrangian framework. The method of finite differences does not enjoy the semi-analytical expression of the solution of the HJ Eq. (3.5) thanks to Lax–Hopf formula (4.12). Hence, we think that monotone or semi-discrete schemes cannot recover shocks waves and contact discontinuities as well as we can hope with Lax–Hopf formula. In the next section, we provide a semi-analytical algorithm to compute the solution based on that representation formula. 4. A ‘‘grid free’’ scheme We are now interested in numerical methods to efficiently recover the solution of the Hamilton–Jacobi problem (3.5). The computation of numerical solutions of the HJ equation has already attracted an important interest in the mathematical community. The majority of numerical schemes which were proposed to solve HJ equations are based on finite difference methods, Semi-Lagrangian schemes and discontinuous Galerkin methods (the interested reader is referred to Claudel and Bayen (2010b), Costeseque et al. (2014) and references therein).

118

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

4.1. Computational strategy The idea to compute the solution of the Eulerian GSOM model (3.9) is to recast it under its Lagrangian form (3.3). As it was shown in the previous section, the position X of vehicle N at time t solves the Hamilton–Jacobi Eq. (3.5). As explained in Lebacque and Khoshyaran (2013), the Initial and Boundary Value Problem (3.5) admits a quite simple representation formula (3.6), very similar to the Hopf–Lax formula presented for the LWR model. The Hopf–Lax formula (3.6) can be simplified because it is well-known from optimal control theory, that the optimal trajectories for which the minimum is attained are the characteristics (3.8). Hence we have to solve the following system of coupled ODEs (3.8). Then the generalized Hopf–Lax formulation (3.6) can be recast as follows

RT XðNT ; TÞ ¼ min t0 MðN; @ r WðN; r; tÞ; tÞdt þ nðN0 ; t 0 Þ; ðN 0 ;r0 ;t 0 Þ   NðtÞ _ ¼ @ r WðN; r; tÞ  _  rðtÞ ¼ @ N WðN; r; tÞ   Nðt 0 Þ ¼ N0 ; rðt0 Þ ¼ r0 ; NðTÞ ¼ NT   ðN ; r ; t Þ 2 K 0 0 0

ð4:12Þ

where K is the set of initial/boundary values obtained by combining the initial values of J with the initial/boundary value r 0 deduced from the initial values nðN 0 ; t 0 Þ. According to the principle of inf-morphism property (Aubin et al., 2008; Aubin et al., 2011; Claudel and Bayen, 2010a; Claudel and Bayen, 2010b), if the initial/boundary condition data n are located on a union (non necessarily disjoint) of sets

K¼

[

Kl ;

l

it suffices to solve partial problems on each set Kl and to compute the minimum of the solutions to these sub-problems

XðNT ; TÞ ¼ minX l ðN T ; TÞ; l

ð4:13Þ

with

RT X l ðNT ; TÞ :¼ min t0 MðNl ; @ r WðNl ; rl ; tÞ; tÞdt þ nðN0 ; t0 Þ: ðN0 ;r 0 ;t0 Þ _  Nl ðtÞ ¼ @ r WðNl ; r l ; tÞ  _  rl ðtÞ ¼ @ N WðN l ; rl ; tÞ   Nl ðt 0 Þ ¼ N0 ; r l ðt 0 Þ ¼ r 0 ; Nl ðTÞ ¼ NT   ðN0 ; r 0 ; t0 Þ 2 Kl

ð4:14Þ

In the remaining of this article, we will apply the inf-morphism property by considering initial and boundary conditions which are piecewise affine (PWA). Thus we will calculate the solution generated by each piece using (4.14) and then apply (4.13) in order to obtain the solution of (3.6). 4.2. Algorithm for piecewise affine value conditions Hereafter, we study separately the different elements which contribute to the value of the solution of (3.5) in its Lagrangian setting. We distinguish  the initial condition at time t ¼ t 0 describing the initial position of vehicles nðN; t0 Þ for any considered N,  the ‘‘upstream’’ boundary condition that is the trajectory of the first vehicle N ¼ N 0 traveling on the section nðN 0 ; tÞ for any considered t,  and internal boundary conditions given for instance by cumulative vehicle counts at fixed location X ¼ x0 . Note that Lagrangian data which are individual vehicle trajectories are considered as a particular case of ‘‘upstream’’ boundary condition for a given N P N 0 (see Remark 4.7). While the Lax–Hopf algorithm can handle infinite horizon problems either in the Eulerian or in the Lagrangian framework, we restrict ourselves to finite values for a convenient numerical implementation. Without loss of generality, we assume that N 2 ½N 0 ; N max  with N max < þ1 and similarly t 2 ½t0 ; tmax  with tmax < þ1. We recall that according to Remark 3.1, in all the following figures, the N-axis and x-axis are increasing in opposite directions. Hence, the highest labels actually match vehicles that are located further upstream. 4.2.1. Initial conditions In this case, at t ¼ t 0 , the positions nðn; t0 Þ of vehicles n are given. We have that

r0 ðNÞ ¼ @ N nðN; t0 Þ;

for any N:

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

119

The initial conditions for the characteristics are the couples ðN; r0 ðNÞÞ. First, we need to discretize the set of vehicle labels into intervals ½np ; npþ1  of length Dn, for p ¼ 1; . . . ; P, with the additional convention

n1 :¼ N0

and nPþ1 :¼ Nmax ;

in such a way that the dynamics u of the driver attribute I can be approximated by up in the interval ½np ; npþ1 

uðN; I; tÞ ¼ up ðI; tÞ; for any N 2 ½np ; npþ1 : If the discrete step Dn is small enough, we could also assume that the initial data are piecewise constant, say



IðN; t0 Þ ¼ I0;p ; rðN; t0 Þ ¼ r0;p

for any N 2 ½np ; npþ1 ;

;

or equivalently we define

8 R < I0;p :¼ D1n nnpþ1 IðN; t 0 ÞdN; p ; R : r 0;p :¼ 1 npþ1 rðN; t 0 ÞdN; Dn np

for any p ¼ 1; . . . ; P:

We will now first define properly the initial condition and then we calculate the solution generated by the pth component (p ¼ 1; . . . ; P). Definition 4.1 (PWA initial condition). Let t0 P 0 be fixed. Then the pth component of the initial condition (p ¼ 1; . . . ; P) is given by

X ini ðN; t 0 Þ ¼ r 0;p ðN  np Þ þ ap ;

for any N 2 ½np ; npþ1 :

Without loss of generality, we state that

a1 ¼ X ini ðN0 ; t0 Þ ¼ 0: To ensure continuity of the initial data X ini on ½N 0 ; N max , we require that

ap ¼

p1 X r0;l Dn;

for any p ¼ 2; . . . ; P:

l¼1

In a first step, we want to compute the characteristics generated by the initial conditions (IðN; t0 Þ ¼ i0;p and rðN; t0 Þ ¼ r 0;p given on the interval ½np ; npþ1  at time t0 ). In the general case, we would have

IðN; tÞ ¼ Ip ðtÞ for any N 2 ½np ; npþ1 : In the following, we consider a characteristic denoted by NðtÞ for t 2 ½t 0 ; t max  and we distinguish two cases according to the location of this characteristic. On the interval N 2np ; npþ1 ½. Consider Xp :¼ ftjNðtÞ 2np ; npþ1 ½g for any p ¼ 1; . . . ; P. We then have to solve the following nonlinear first order ODE

(

I_p ðtÞ ¼ up ðIp ðtÞÞ;

for t 2 Xp ;

Ip ðt 0 Þ ¼ I0;p : On the interval np ; npþ1 ½, we have that IðN; tÞ ¼ Ip ðtÞ is independent of N i.e. @ N I ¼ 0. Then from (3.8), we deduce

r_ ¼ @ N WðN; r; tÞ ¼ @ I VðN; IðN; tÞÞ @ N IðN; tÞ ¼ 0: Thus rðN; tÞ on any characteristic in Xp . Hence it simply suffices to solve the second ODE in (3.8) that is

N_ ¼ @ r WðN; r; tÞ to get the equation of the characteristic on the interval np ; npþ1 ½. We finally get

NðtÞ ¼

Z

t

@ r WðN; r; sÞds þ NðuÞ;

for any t; u 2 Xp :

ð4:15Þ

u

To considerably simplify the presentation, we will restrict ourselves to the case of systems with no relaxation that is the dynamics of driver attributes up ðI; tÞ ¼ 0 at any time t P t0 and for any interval ½np ; npþ1  with p ¼ 1; . . . ; P. If we assume that the dynamics is null, then it obviously follows that I_p ¼ 0 which leads to

Ip ðtÞ ¼ I0;p

for any t 2 Xp :

ð4:16Þ

120

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

Remark 4.2 (Particular choice of u ¼ 0). In traffic modeling, this choice is particularly relevant when the driver attribute I represents vehicle (or driver) properties that do not evolve in time, for instance vehicle type (car, truck, bus) or car mode (taxi, electric vehicle, etc.), vehicle emission or consumption rates, driver’s origin–destination pattern, driver’s turning preference for junction modeling and so on. The driver attribute is conserved at any time and along any characteristic inside each ‘‘strip’’ ½np ; npþ1   ½t0 ; tmax  of width Dn. By the way, the simplification (4.16) allows to consider characteristics which are straight lines while there are not in the general case. At the edge N ¼ np . The edge N ¼ np separates two different traffic states that exhibit two different speed-spacing fundamental diagrams. We need to take care of the kinematic waves that can appear at this locus. As illustrated on Fig. 3, we set the characteristic speed related to the spacing r0;p as follows

mp :¼ @ r Vðr0;p ; I0;p Þ; and we also define the speed of the ‘‘refracted’’ characteristic wave through the discontinuity of I such that





mp :¼ @ r V r0;p ; I0;p :

ð4:17Þ

This characteristic speed is related to a new spacing value r0;p defined as follows

r 0;p

¼

8  < V 1 Vðr 0;p1 ; I0;p1 Þ; I0;p if : þ1

Vðr0;p1 ; I0;p1 Þ < supVðr; I0;p Þ; r2R

ð4:18Þ

else:

Remark 4.3 (Transported spacing value along a characteristic). We extend the idea of what happens at each discontinuity of I by introducing a function

s : ðr0;p ; qÞ # rq0;p

ð4:19Þ

that maps a given (initial) spacing r 0;p to its projection rq0;p (which is the spacing value transported along the corresponding characteristic) in a strip ½nq ; nqþ1   ½t0 ; tmax , with any q P p and p 2 f1; . . . ; Pg. The transported value can be computed thanks to a recursive composition of function s on consecutive strips

rp0;q :¼ sðr0;p ; qÞ ¼ sðsð. . . sðr 0;p ; p þ 1Þ; . . . ; q  1Þ; qÞ;   where rpþ1 0;p :¼ s r 0;p ; p þ 1 ¼ r 0;pþ1 can be easily computed according to (4.18), for any p 2 f1; . . . ; Pg. Indeed, the speed is conserved through a discontinuity of I i.e.



 qþ1 V r q0;p ; I0;q ¼ V r 0;p ; I0;qþ1 ;

for any q P p:

Obviously we also have that rp0;p :¼ sðr0;p ; pÞ ¼ r 0;p . Regarding the values of the initial spacing r0;p and r0;p1 , we can distinguish two cases that could occur starting to an edge N ¼ np :  either mp > mp and in this case there will be a shock wave. The characteristics carrying respectively the initial states r 0;p1 and r0;p will cross each other, offering a partial superposition of both characteristics domains (see (a) on Fig. 4). The solution is then simply computed thanks to the inf-morphism property. Hence we only consider the two characteristics (indicated as r and s on Figs. 6, 7 and 9) that encompass the domain of influence generated by the initial spacing r 0;p .  or mp < mp and in this case there will be a rarefaction fan. The characteristics waves carrying the initial states r 0;p1 and r 0;p diverge (see (b) on Fig. 4). The solution on the area between both extreme characteristics may not be computed and even the inf-morphism property will not be able to recover it. That is why we need to consider the first two characteristics (as above) and a third wave corresponding to the first characteristic matching with the initial state r0;p1 (labeled { on Figs. 6, 7 and 9).

Fig. 3. Schematic view of what happens for the characteristic wave generated from the edge N ¼ np and passing through a discontinuity of I.

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

121

Fig. 4. Characteristics through a discontinuity of I of r.

Notice that by convention, for a given p ¼ 1; . . . ; P, we add the rarefaction wave { separating the two states at N ¼ np (if any) to the solution generated by the initial condition on ½np ; npþ1 . Remark 4.4 (Degenerate case for characteristics through a discontinuity of I). When passing through a discontinuity of I (assume ðaÞ and ðbÞ the states on both sides of the discontinuity as illustrated on Fig. 5), the characteristics speed may be changed. Nevertheless it is well known (Rankine–Hugoniot jump condition) that through such a discontinuity, the traffic speed should be unmodified

v ðaÞ :¼ VðrðaÞ ; IðaÞ Þ ¼ VðrðbÞ ; IðbÞ Þ ¼: v ðbÞ : ðbÞ

In case of v ðaÞ > v max :¼ supr2R Vðr; IðbÞ Þ, it is obvious that the equality between these two speeds cannot be complied. In an Eulerian setting, it means that the group of upstream vehicles ðbÞ cannot accelerate enough (even increasing their speeds ðbÞ to v max ) to catch up with the downstream vehicles ðaÞ. This induces the apparition of a vacuum area between both groups of vehicles. There is a rarefaction wave as vehicles ðbÞ accelerate to attain their top speed and a shock wave following the last vehicle of type ðaÞ. Notice that the vacuum states ða1 Þ and ðb1 Þ defined on Fig. 5 coincide for r ¼ þ1 (or equivalently at q ¼ 0). It is noteworthy that the vacuum area is not visible in Lagrangian viewpoint (see Fig. 5 (b)) while it matches the rarefaction fan ðbÞ and characteristics of type (a) in Eulerian viewpoint (see the area delimited by characteristic with speed v max and dashðaÞ dotted characteristic with speed v on Fig. 5 (d)). The contribution X p of the initial condition defined on ½np ; npþ1   ft0 g for any p 2 f1; . . . ; Pg is then computed as follows: (i) Initialize the partial solution X p at þ1 on the whole computational domain ½N 0 ; N max   ½t0 ; tmax . (ii) Determine the number of characteristics to compute (two or three) according to what could occur at the edge N ¼ np . (iii) Compute the equation NðtÞ of each characteristic while t 6 tmax and N 6 N max .

Fig. 5. Critical ‘‘vacuum’’ case appearing from special condition values.

122

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

(iv) Calculate the (exact) solution X p all along each computed characteristic generated by the interval ½np ; npþ1   ft0 g, namely characteristics r and s, and { whenever it appears (see Fig. 6). To this aim, we use the generalized Lax–Hopf formula (4.14) which gives that

_ tÞ: X_ p ðtÞ ¼ MðNðtÞ; NðtÞ; The interested reader is referred to Brenier (1989), Corrias (1996), Lucet (1997) for additional information about the Legendre–Fenchel transform and fast algorithms for its numerical computation. (v) Compute the exact value at any point within the characteristics strip (delimited by characteristics r and s on Fig. 6) using the fact that for any point ðN; tÞ belonging to the characteristic strip, the position at this point can be deduced by a simple translation h iof the position on characteristic r (see Fig. 6). Indeed we have for any N 2 ½nq ; nqþ1  with q P p ð2Þ ð1Þ and any t 2 t N ; t N

 Z ð1Þ X p ðN; tÞ ¼ X p N; t N þ

t

t ð1Þ

VðrðN; sÞ; IðN; sÞÞds ð4:20Þ

N





 ð1Þ ð1Þ ¼ X p N; tN þ V rq0;p ; I0;q t  t N ;

ðiÞ

where we recall that up ðI; tÞ ¼ 0 for any couple ðI; tÞ and r q0;p is computed according to (4.19). The time t N ; i ¼ 1; 2 corresponds to the time when the characteristic  i crosses the line fNg  ½t0 ; tmax . (vi) In the case of a rarefaction fan, evaluate the value of X p at each point within the influence domain of the considered initial condition (illustrated on Fig. 6) by an interpolation technique based on triangular meshes (the value at each triangle vertex is exact).

Algorithm 1. Pseudo-code for the computation of X p on the computational domain ½N 0 ; N max   ½t0 ; tmax , under the initial condition X ini . Input: N 0 ; N max ; t0 ; t max ; Dn, N p and r0;p 1: X p þ1 2: Compute lp and r 0;p according to (4.17) and (4.18) 3: if lp > lp then 4: ncharac ¼ 2 5: else 6: ncharac ¼ 3 7: end if 8: for j ¼ 1 to ncharac do 9: for t ¼ t 0 to tmax do 10: Compute N ðjÞ ðtÞ using (4.15) 11: Compute X p ðN ðjÞ ðtÞ; tÞ using (4.14) 12: end for 13: end for 14: for N ¼ N p to N max do ð2Þ ð1Þ 15: for t ¼ t N to tN do 16: Compute X ðN; tÞ using (4.20) 17: end for 18: end for 19: if ncharac ¼ 3 then 20: for N ¼ N p to N max do ð1Þ ð0Þ 21: for t ¼ t N to tN do 22: Compute X ðN; tÞ using interpolation on the three closest exact values 23: end for 24: end for 25: end if Output: X p ð:; :Þ over ½N 0 ; N max   ½t 0 ; t max 

{Input label domain, time domain and initial condition} {Initialization of the position function to infinity} {Determination of the number of characteristic curves to compute}

{Iteration on the characteristics} {Iteration on time} {Computation of the characteristic equation} {Computation of the position along the characteristic}

{Iteration on greater labels} {Iteration on time} {Computation of the position within the characteristics strip}

{Case of the rarefaction fan} {Iteration on greater labels} {Iteration on time} {Computation of the position within the rarefaction fan}

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

123

Fig. 6. Domain of influence of the initial data ½np ; npþ1   ft 0 g.

Fig. 7. Domain of influence of the upstream boundary data fN ¼ N 0 g  ½tq ; tqþ1 .

4.2.2. ‘‘Upstream’’ boundary conditions In the Lagrangian setting, boundary conditions describe floating vehicle conditions in Eulerian setting. Indeed such a condition is equivalent to consider the trajectory t # nðN 0 ; tÞ of a given vehicle N 0 . Then we have that

v 0 ðtÞ :¼ @ t nðN0 ; tÞ ¼ Vðr0 ðN0 ; tÞ; IðN0 ; tÞÞ;

for any t P t 0 ;

which can be solved and yields a unique solution (because V is strictly increasing)

r0 ðtÞ ¼ V 1 ðv 0 ðtÞ; IðN0 ; tÞÞ;

for any t P t 0 :

The initial conditions for the characteristics are the couples ðN 0 ; r 0 ðtÞÞ. As for initial condition, we introduce a discrete time step Dt and we consider a full discretization of the time domain ½t 0 ; t max  into segments ½tq ; tqþ1  of length Dt with q ¼ 1; . . . ; Q . We have by convention

tq¼1 :¼ t 0

and t Qþ1 :¼ t max :

We still consider that IðN; tÞ is piecewise constant w.r.t. N and constant w.r.t. time t (u ¼ 0) such that IðN 0 ; tÞ ¼ I0;p¼1 . If Dt is small enough, we can consider that for any q 2 f1; . . . ; Q g, we have



v 0 ðtÞ ¼: v 0;q ; r 0 ðtÞ ¼ V 1 ðv 0;q ; I0;p¼1 Þ ¼: r 0;q ;

for any t 2 ½t q ; tqþ1 :

Definition 4.5 (PWA upstream boundary condition). Let N 0 2 Z be fixed. Then the qth component of the upstream boundary condition (q ¼ 1; . . . ; Q ) is given by

X up ðN0 ; tÞ ¼ v 0;q ðt  t q Þ þ bq ;

for any t 2 ½tq ; t qþ1 :

Without loss of generality, we set b1 such that

b1 :¼ X up ðN0 ; t 0 Þ ¼ 0: To ensure continuity of the upstream boundary data X ini on ½t 0 ; t max , we require that

bq ¼

q1 X l¼1

v 0;l Dt;

for any q ¼ 2; . . . ; Q :

124

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

The contributions of upstream boundary conditions are computed in a similar way than those of initial condition described in the previous subsection. The domain of influence of upstream boundary conditions is simply separated into what happens on each interval tq ; tqþ1 ½ and at the edge t q for any q 2 f1; . . . ; Q g (see Fig. 7). It is worth noting that we have a rarefaction wave at t ¼ t q for any q 2 f1; . . . ; Q g if and only if v 0;q1 < v 0;q (or equivalently r0;q1 < r0;q ). Note that in this case, the rarefaction wave denoted { separating the two states at t ¼ tq (see Fig. 7) is added by convention, to the solution computed from the boundary condition on ½tq ; tqþ1 . Algorithm 2. Pseudo-code for the computation of X q on the computational domain ½N 0 ; N max   ½t0 ; t max , under the upstream boundary condition X up . Input: N 0 ; N max ; t0 ; t max ; Dt, tq and r 0;q 1: X q þ1 2: if r 0;q1 > r0;q then 3: ncharac ¼ 2 4: else 5: ncharac ¼ 3 6: end if 7: for j ¼ 1 to ncharac do 8: for t ¼ t 0 to tmax do 9: Compute N ðjÞ ðtÞ using (4.15) 10: Compute X p ðN ðjÞ ðtÞ; tÞ using (4.14) 11: end for 12: end for 13: for N ¼ N p to N max do ð1Þ ð2Þ 14: for t ¼ t N to tN do 15: Compute X ðN; tÞ using (4.20) 16: end for 17: end for 18: if ncharac ¼ 3 then 19: for N ¼ N p to N max do ð0Þ ð1Þ 20: for t ¼ tN to t N do 21: Compute X ðN; tÞ using interpolation on the three closest exact values 22: end for 23: end for 24: end if Output: X q ð:; :Þ over ½N 0 ; N max   ½t 0 ; t max 

{Input label domain, time domain and initial condition} {Initialization of the position function to infinity} {Determination of the number of characteristic curves to compute}

{Iteration on the characteristics} {Iteration on time} {Computation of the characteristic equation} {Computation of the position along the characteristic}

{Iteration on greater labels} {Iteration on time} {Computation of the position within the characteristics strip}

{Case of the rarefaction fan} {Iteration on greater labels} {Iteration on time} {Computation of the position within the rarefaction fan}

4.2.3. Internal boundary conditions Note that internal boundary condition is here understood in the Lagrangian framework. It does not correspond to internal vehicle trajectories which can also be incorporated in the algorithm (see the previous subsection and Remark 4.7). Notice that internal boundary condition is called mixed condition in Lebacque and Khoshyaran (2013) to avoid any misunderstanding. We assume that the data comes from vehicles NðtÞ located at a point nðNðtÞ; tÞ at time t. Thus we get

d _ nðNðtÞ; tÞ ¼ NðtÞr 0 ðtÞ þ WðNðtÞ; r 0 ðtÞ; tÞ: dt

ð4:21Þ

As not all data are compatible (in the sense that they need to comply to the traffic dynamics laws), we need to introduce compatibility conditions. From Legendre–Fenchel transform, we have

n o _ _ þ WðNðtÞ; r; tÞ ¼ MðNðtÞ; NðtÞ; tÞ; sup NðtÞr r2R

hence the equation yielding r 0 ðtÞ (4.21) admits a solution only if

d _ tÞ: nðNðtÞ; tÞ 6 MðNðtÞ; NðtÞ; dt

ð4:22Þ

_ _ P 0 for any time t. This case occurs when the data originates at a ¼ 0 and NðtÞ Let us describe the specific situation of nðtÞ fixed measurement point x0 like a fixed detector data that is nðNðtÞ; tÞ ¼ x0 for all time t. Then NðtÞ represents the cumulative

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

125

flow at point x0 . In this case the above compatibility condition (4.22) is satisfied if moreover we assume that N_ 6 Fmax ðIÞ. It is _ quite natural to impose such a condition because NðtÞ is equal to the instantaneous traffic flow which is evidently bounded by Fmax ðIÞ (see Fig. 8). Moreover such a condition ensures that the Lagrangian stays non-negative which is coherent with its interpretation as a instantaneous speed cost. As illustrated on Fig. 8, there are two solutions (except in the very particular case of N_ ¼ Fmax ), one under-critical and the other over-critical according to classical traffic definitions. Recall that the speed-spacing fundamental diagram W depends on the vehicle N. If we propagate both pairs of characteristics, then the inf-morphism property will automatically select a single solution which matches the over-critical situation (because the speed is lower than in the under-critical situation). It may happen that the congested solution is not the good one. Then to avoid any mistake, we assume that the fixed detector gives the cumulative flow and also the traffic flow speed. In this way, we can select one traffic state only. In our case, we are interested in including some Eulerian data coming from classical fixed sensors like inductive loop detectors. To achieve this aim, Eulerian data become internal boundary condition into the Lagrangian framework. We assume that a Eulerian sensor located at a fixed position x0 gives us the incremental cumulative vehicle count which is then _ interpolated in a piecewise affine function NðtÞ for t 2 ½t 0 ; t max . We then define by f0;p the value of NðtÞ for any t 2 ½tp ; tpþ1  and we set ^t p the time such that Nð^t p Þ ¼ np for any p 2 f1; . . . ; Pg. Definition 4.6 (PWA internal boundary condition). Let x0 P 0 be fixed. Then the pth component of the internal boundary condition (p ¼ 1; . . . ; P) is given by

( X int ðn; tÞ ¼ x0 ;

with

 t ¼ f0;p n  np þ cp ; n 2 ½np ; npþ1 :

Without loss of generality, we consider that c1 ¼ ^t1 P t 0 is given. To ensure continuity of the trajectory t # nðtÞ on which the internal boundary data X int is prescribed, we require that

cp ¼ c1 þ

p1 X f0;l ðnlþ1  nl Þ;

for any p ¼ 2; . . . ; P:

l¼1

We assume that N is piecewise affine on each discrete segments ½np ; npþ1  that is equivalent to say that N_ is piecewise constant (i.e. constant in each strip ½np ; npþ1   ½t 0 ; t max  for p ¼ 1; . . . ; P). It is easy to deal with this case in the algorithm because the computational steps are similar to both previous cases for initial and upstream boundary conditions. The only difference resides in the fact that we have a characteristic strip that matches either the under-critical or the over-critical traffic state. It is simple to verify that in the first (under-critical) case, the pth component generates characteristics which are emitted with speeds munder > N_ while it is the inverse in the other case (see Fig. 9). From a traffic point of view, it is relevant with observations stating that a deceleration wave spreads over vehicles located further upstream at a slower speed than an acceleration wave. Hence, from a computational viewpoint, if we are in the under-critical (resp. over-critical) case, then we will use the same strategy than Algorithm 1 (resp. Algorithm 2). Remark 4.7 (Extensions). We consider some extensions to the previous data assimilation:  Lagrangian boundary conditions: The Algorithm 2 for upstream boundary conditions can be extended to any individual trajectory condition. Any Lagrangian data giving the trajectory of a vehicle N I with N I 2 ½N 0 ; N max  will strictly generate the same kind of calculations than the ones described for the upstream boundary condition (see SubSection 4.2.2). It is then easy to deal with data coming from mobile sensors moving within the traffic stream.  Eulerian ‘‘downstream’’ boundary conditions: It is noteworthy that downstream boundary condition in the Eulerian framework, say data coming from a fixed location xdown , can be accommodated as an internal boundary condition in the Lagrangian framework (see SubSection 4.2.3).

_ _ Fig. 8. Existence of two solutions corresponding to a condition nðtÞ ¼ 0 and NðtÞ P 0.

126

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

Fig. 9. Domain of influence of an internal boundary data (in red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The different elements of the computations for piecewise affine value conditions are recapitulated in the following pseudo-code (see Algorithm 3). Algorithm 3. Pseudo-code implementation for the Lax–Hopf based computation of the position function X on the computational domain ½N 0 ; N max   ½t 0 ; tmax  prescribed by the user.

Input: N 0 ; N max ; t0 and t max 1: X þ1 2: for p ¼ 1 to P do 3: for T ¼ t 0 to tmax do 4: for N T ¼ N p to N max do 5: compute X ini p ðN T ; TÞ using Algorithm 1 6: if X ini p ðN T ; TÞ < X then 7: X X ini p ðN T ; TÞ 8: end if 9: end for 10: end for 11: end for 12: for q ¼ 1 to Q do 13: for T ¼ t q to tmax do 14: for N T ¼ N 0 to N max do 15: compute X up q ðN T ; TÞ using Algorithm 2 16: if X up q ðN T ; TÞ < X then 17: X X up q ðN T ; TÞ 18: end if 19: end for 20: end for 21: end for 22: for p ¼ 1 to P do 23: for T ¼ t p to tmax do 24: for N T ¼ N p to N max do 25: compute X int p ðN T ; TÞ 26: if X int p ðN T ; TÞ < X then 27: X X int p ðN T ; TÞ 28: end if 29: end for 30: end for 31: end for Output: X ð:; :Þ over ½N 0 ; N max   ½t0 ; tmax 

{Input label domain, time domain} {Initialization of the position function to infinity} {Iteration on the initial conditions} {Iteration on time} {Iteration on greater labels} {Component induced by the initial condition r0;p } {Update the position function}

{Iteration on upstream boundary conditions} {Iteration on greater times} {Iteration on labels} {Component induced by the upstream boundary condition X N0 ; tq } {Update the position function}

{Iteration on internal boundary conditions} {Iteration on greater times} {Iteration on greater labels} {Component induced by the internal boundary condition X Np ; tp } {Update the position function}

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

127

Remark 4.8. This pseudo-code underlines the property of the Lax–Hopf algorithm to treat in a parallel way the different value conditions (initial, upstream boundary and internal boundary). 4.3. Instantiation In order to simplify the computations and to ease the presentation of the following example, we consider that the driver attribute I is piecewise constant with respect to N at initial time t ¼ t 0 . It means that there exist some platoons of vehicles which share the same driver attribute. It could be for example some vehicles of the same kind (cars and trucks), or vehicles that go to the same destination, or vehicles that have the same desired maximal speeds. Actually, in this numerical example, we have considered Dn ¼ 1 so the attribute I is constant for each vehicle. This choice is particularly convenient for our example but it is not restrictive. As u ¼ 0, the individual attribute is conserved at any time t P t0 (see Fig. 11 (left)). We consider for this numerical example a Colombo 1-phase model (Lebacque et al., 2007) given by

( Wðr; N; tÞ :¼ Vðr; IÞ ¼

V  br

if r P r crit ðIÞ;  1 ðI þ q rÞ 1  rR else;

where b; q ; V and R are given parameters of the model. More precisely, q stands for the theoretical maximal flow, while V and R describe respectively the maximal speed and density. The critical spacing rcrit ðIÞ separating the free and congested cases is computed as the inverse of the critical density

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi "  ffi#

2 1 q q I  V þ I : qcrit ðIÞ :¼ V þ  I  4q b  R R R 2 b  RI The numerical values considered in this example are given below

8 > V ¼ 25 m s1 > > > R ¼ 1=5 veh m1 > > > :b R ¼ 1=30 veh m1

8 < b ¼ VbV ¼ 90 m2 s1 bR and : q ¼ 1 veh s1

The flow-density and speed-spacing fundamental diagrams are given on Fig. 10. Notice that these numerical values are satisfactory for a two-lanes road section. The driver attribute I varies between 0 and Imax ¼ 5 which is designed such that the flow-density fundamental diagram Fðq; IÞ stays concave with respect to q. The Lagrangian function M is given by

( Mðu; IðN; tÞÞ :¼ sup

) sup g 1 ðr; u; IÞ;

rPr crit ðIÞ

sup r min 6r6r crit ðIÞ

g 2 ðr; u; IÞ

with

Fig. 10. Speed-spacing fundamental diagram Vðr; IÞ (left) and flow-density fundamental diagram Fðq; IÞ (right).

128

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

(

g 1 ðr; u; IÞ :¼ V  br  ru;  1 g 2 ðr; u; IÞ :¼ ðI þ q r Þ 1  rR  ru:

Some simple algebra leads to

8 qffiffi b > < g 1 ðr crit ðIÞÞ if < rcrit ðIÞ; p sup g 1 :¼

qffiffi > b rPr crit ðIÞ : g1 else p and

8 > > g 2 ðr min Þ if > > > < sup g 2 :¼ g 2 ðr crit ðIÞÞ if > r min 6r6r crit ðIÞ >

qffiffi > > b > : g2 if p

q < p and

qffiffiffiffiffiffiffiffiffiffiffiffi

I Rðpq Þ

< r min ;

qffiffiffiffiffiffiffiffiffiffiffiffi I q < p and > rcrit ðIÞ; Rðpq Þ qffiffiffiffiffiffiffiffiffiffiffiffi I q < p and r min 6 Rðpq  Þ 6 r crit ðIÞ: 

q P p or if



The Lagrangian M is then computed for our numerical instantiation, as illustrated on Fig. 11 (right). We consider piecewise affine initial conditions that is we prescribe each vehicle position Xðt; nÞ at time t ¼ t 0 for n 2 ½N 0 ; N max  according to a normal distribution of spacing between rmin ¼ 5 m and a 50 m ‘‘maximal’’ spacing. The initial condition for this numerical test are plotted on Fig. 12. Moreover we prescribe piecewise affine upstream boundary condition that is the trajectory Xðt; nÞ of the first vehicle N ¼ N 0 for t 2 ½t 0 ; t max . This upstream boundary condition is illustrated on Fig. 13. We also have included two (fictive) vehicle trajectories for N ¼ 25 and N ¼ 46 as Lagrangian internal boundary conditions. They are plotted on Fig. 14. For these trajectories, we have considered randomly distributed speed values that are piecewise constant. Note that the speed distribution is chosen such that the mean value is equal to 18 m s1 in a first time and then to 10 m s1. In that way, one can verify the trajectories are piecewise affine and averaging on time, vehicles are globally decelerating. 4.4. Numerical result and interpretation The solution is thus given by applying the computation algorithm described in Algorithm 3. The simulations deal with 50 vehicles on a 2 min run. We also compare the solution with initial and upstream boundary conditions and the solution with initial, upstream and internal boundary conditions. The numerical solutions X are respectively plotted on Fig. 15 (a) and (b). In Fig. 15 (c) and (d), we have also plotted the speeds of vehicles which are computed as the discrete gradient of the solution X ðN; tÞ according to time t and for any vehicle N. One can observe the shock and rarefaction waves generated from initial and boundary conditions. Note that the internal boundary conditions are not necessarily defined for every t 2 ½t 0 ; t max . It is the reason why there is an area up to N ¼ 46 and t ¼ 110 that matches a high speed rarefaction fan since the solution is no more constrained by the internal boundary condition. One can note that the integration of two slower vehicles (N ¼ 25 and N ¼ 46) in the traffic flow implies the modification of the solution. The values of X for N P 25 and then N P 46 are notably decreased. Indeed slower vehicles impose lower

Fig. 11. Initial values of driver attribute IðN; t0 Þ (left) and Lagrangian function MðN; u; tÞ (right).

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

129

Fig. 12. Initial conditions for the GSOM PDE at t ¼ t0 .

Fig. 13. Boundary conditions for the GSOM PDE at N ¼ N 0 .

speeds to the following vehicles as it could be noticed on vehicle trajectories plotted on Fig. 16. Thus, data assimilation allows to modify the previous over-estimation we made by considering only initial and upstream boundary conditions. We also include a fictitious Eulerian data, i.e. Cumulative Vehicle Counts NðtÞ, coming from a fixed location x0 ¼ 0. The numerical solution is plotted on Fig. 17 (a) and the corresponding speed locus are given on Fig. 17 (b). One can see that according to the traffic speed, the domain of influence is either above (under-critical case) or below (over-critical case) the Eulerian curve. As it is clearly visible for the speeds (see Fig. 17 (b)) in comparison to the ‘‘basic’’ case (see Fig. 15 (d)), the propagation of under-critical values (resp. over-critical values) matches with high speeds (resp. low speeds). Vehicle trajectories are also corrected according to the position values deduced from the Eulerian data. In the numerical example above, we do not consider real data assimilation. The next step for numerical results should be the assimilation of real data in a model of the GSOM family. As it was described in the previous section, vehicle trajectories and fixed detector measurements can be easily incorporated in the algorithm. However we need to identify a driver attribute I and to know how it evolves on the computational domain. That is the tough point because it is hard to measure or to evaluate. Notice that an example of speed data assimilation with the ARZ model which is a model of the GSOM family can be found in Mammar et al. (2009). Concerning the identification of I, a method could be deduced from Fan et al. (2013), exploiting data-fitted fundamental diagrams. In the case of a scalar attribute I, we can use for instance speed and spacing measurements to evaluate the value of I.

130

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

Fig. 14. Lagrangian internal boundary conditions for the GSOM PDE at N ¼ 25 and N ¼ 46.

5. Discussion and conclusion In this paper, we are interested in the assimilation of different kinds of data into a generic class of macroscopic traffic flow models for improving traffic state estimation. The data come from mixed sources including Lagrangian vehicle trajectories and Eulerian cumulative vehicle counts obtained from fixed detectors. We then consider macroscopic traffic flow models of the Generic Second Order Modeling (GSOM) family in the Lagrangian system of coordinates which seems to be the most adapted framework for dealing with such data. We describe a computational method providing solutions to these models under piecewise affine initial and boundary conditions. The numerical method is based on the variational theory which has been extended to the GSOM family very recently (Lebacque and Khoshyaran, 2013). Solution can be computed thanks to Lax–Hopf like formula and a generalization of the inf-morphism property (Claudel and Bayen, 2010a). It is possible to considerably reduce the number of integral curves by following only the characteristics which are the optimal trajectories. Another computational trick is to solve separately many partial problems (by discretizing initial and boundary conditions) instead of the general problem. We are aware that there exists some natural relationships between our method based on the Hamilton–Jacobi formulation of the models and shock-fitting (Wong and Wong, 2002) or front-tracking (Henn, 2005) methods for hyperbolic systems of conservation laws. This could be an interesting avenue of investigation. However we have not deepened these relationships in the present work. One reason is that there is a multiplication of the number of waves in the framework of GSOM

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

Fig. 15. Numerical solution for GSOM PDE obtained on the computation domain ½N 0 ; N max   ½t 0 ; tmax .

Fig. 16. Vehicle trajectories obtained for integers of the GSOM numerical solution Xðt; nÞ.

131

132

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

Fig. 17. Numerical solution for GSOM PDE obtained on the computation domain ½N 0 ; N max   ½t 0 ; tmax  with Lagrangian and Eulerian data assimilation.

models, since the fundamental diagram is driver-dependent. Another reason is that for such methods, we need to restrict ourselves to piecewise linear fundamental diagrams while it is not the case with our methodology. Extensions of this work could include more general assumptions on the dynamics u – 0 of driver attributes, always with the assumption of piecewise affine value conditions. We are convinced that the adaptation of our methodology is quite straightforward for Lagrangian source term uðN; tÞ – 0. It is more challenging for Eulerian source term, say /ðx; tÞ – 0, that can model junctions for instance. It is the subject of active ongoing researches. Additionally, the computational benefits of the variational method over the finite difference method which is also presented in the paper need to be illustrated on several comparisons. We expect that the variational method has a lower computational cost (due to parallel computing e.g.) and a globally higher accuracy compared with the other method. A still open but interesting question is the application of the algorithm to data assimilation on road networks. The challenging point is to deal with fundamental diagrams which will depend on the drivers attribute but also on the spatial position. To our best knowledge Lax–Hopf formula seems not to be so tractable in that framework because the formula does not simplify so much for space dependent Hamiltonians. One considered possibility is to use the extended concepts of supply and demand functions to models of the GSOM family (Lebacque et al., 2008). This subject is being currently investigated by the authors. Aknowledgments The authors are very grateful to both anonymous reviewers that give valuable comments and insights on the first version of this paper. The first author would like to thank Jia Li, PhD candidate at UC Davis, California, USA, for deeper insight of his work. This work was partially supported by the ANR (Agence Nationale de la Recherche) through HJnet project ANR-12-BS01– 0008-01. References Aubin, J.-P., Bayen, A.M., Saint-Pierre, P., 2008. Dirichlet problems for some Hamilton–Jacobi equations with inequality constraints. SIAM Journal on Control and Optimization 47 (5), 2348–2380. Aubin, J.-P., Bayen, A.M., Saint-Pierre, P., 2011. Viability Theory: New Directions. Springer. Aw, A., Rascle, M., 2000. Resurrection of ‘‘second order’’ models of traffic flow. SIAM Journal on Applied Mathematics 60 (3), 916–938. Bagnerini, P., Rascle, M., 2003. A multiclass homogenized hyperbolic model of traffic flow. SIAM Journal on Mathematical Analysis 35 (4), 949–973. Brenier, Y., 1989. Un algorithme rapide pour le calcul de transformées de Legendre–Fenchel discretes. Comptes rendus de l’Académie des sciences. Série 1, Mathématique 308 (20), 587–589. Claudel, C.G., Bayen, A.M., 2010a. Lax–Hopf based incorporation of internal boundary conditions into Hamilton–Jacobi equation. Part I: Theory. IEEE Transactions on Automatic Control 55 (5), 1142–1157. Claudel, C.G., Bayen, A.M., 2010b. Lax–Hopf based incorporation of internal boundary conditions into Hamilton–Jacobi equation. Part II: computational methods. IEEE Transactions on Automatic Control 55 (5), 1158–1174. Colombo, R.M., 2003. Hyperbolic phase transitions in traffic flow. SIAM Journal on Applied Mathematics 63 (2), 708–721. Corrias, L., 1996. Fast Legendre–Fenchel transform and applications to Hamilton–Jacobi equations and conservation laws. SIAM Journal on Numerical Analysis 33 (4), 1534–1558. Costeseque, G., Lebacque, J.-P., Monneau, R., 2014. A convergent scheme for Hamilton–Jacobi equations on a junction: application to traffic. Numerische Mathematik, 1–43. Daganzo, C.F., 1995. The cell transmission model, part II: network traffic. Transportation Research Part B 29 (2), 79–93.

G. Costeseque, J.-P. Lebacque / Transportation Research Part B 70 (2014) 112–133

133

Daganzo, C.F., 2005a. A variational formulation of kinematic waves: basic theory and complex boundary conditions. Transportation Research Part B 39 (2), 187–196. Daganzo, C.F., 2005b. A variational formulation of kinematic waves: solution methods. Transportation Research Part B 39 (10), 934–950. Daganzo, C.F., 2006. On the variational theory of traffic flow: well-posedness, duality and applications. Networks and Heterogeneous Media 1 (4), 601–619. Daganzo, C.F., Menendez, M., 2005. A variational formulation of kinematic waves: Bottleneck properties and examples. In: Transportation and Traffic Theory. Flow, Dynamics and Human Interaction. 16th International Symposium on Transportation and Traffic Theory. Duret, A., Buisson, C., Chiabaut, N., 2008. Estimating individual speed-spacing relationship and assessing ability of Newell’s car-following model to reproduce trajectories. Transportation Research Record: Journal of the Transportation Research Board 2088 (1), 188–197. Evans, L.C., 2009. Partial Differential Equations (Graduate Studies in Mathematics, Vol. 19), American Mathematics Society. Fan, S., Herty, M., Seibold, B., 2013. Comparative Model Accuracy of a Data-fitted Generalized Aw-Rascle-Zhang model. arXiv Preprint arXiv:1310.8219. Fan, S., Seibold, B., 2012. A Comparison of Data-fitted First Order Traffic Models and Their Second Order Generalizations via Trajectory and Sensor Data. arXiv Preprint arXiv:1208.0382. Garavello, M., Piccoli, B., 2006. Traffic Flow on Networks. American institute of mathematical sciences, Springfield, MO, USA. Greenberg, J.M., 2002. Extensions and amplifications of a traffic model of Aw and Rascle. SIAM Journal on Applied Mathematics 62 (3), 729–745. Han, K., Yao, T., Friesz, T.L., 2012. Lagrangian-based Hydrodynamic Model: Freeway Traffic Estimation. arXiv preprint arXiv:1211.4619. Henn, V., 2005. A wave-based resolution scheme for the hydrodynamic LWR traffic flow model. In: Traffic and Granular Flow03. Springer, pp. 105–124. Herrera, J.C., Bayen, A.M., 2010. Incorporation of lagrangian measurements in freeway traffic state estimation. Transportation Research Part B 44 (4), 460– 481. Herrera, J.C., Work, D.B., Herring, R., Ban, X.J., Jacobson, Q., Bayen, A.M., 2010. Evaluation of traffic data obtained via GPS-enabled mobile phones: the Mobile Century field experiment. Transportation Research Part C 18 (4), 568–583. Herty, M., Kirchner, C., Moutari, S., Rascle, M., 2008. Multicommodity flows on road networks. Communications in Mathematical Sciences 6 (1), 171–187. Holden, H., Risebro, N.H., 2011. Front Tracking for Hyperbolic Conservation Laws, vol. 152. Springer. Hoogendoorn, S.P., Bovy, P.H., 2001. State-of-the-art of vehicular traffic flow modelling. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering 215 (4), 283–303. Jin, W., Zhang, H.M., 2004. Multicommodity kinematic wave simulation model for network traffic flow. Transportation Research Record: Journal of the Transportation Research Board 1883 (1), 59–67. Khoshyaran, M.M., Lebacque, J.-P., 2009. A stochastic macroscopic traffic model devoid of diffusion. In: Traffic and Granular Flow07. Springer, pp. 139–150. Klar, A., Greenberg, J., Rascle, M., 2003. Congestion on multilane highways. SIAM Journal on Applied Mathematics 63 (3), 818–833. Kuhne, R., Michalopoulos, P., 1997. Continuum Flow Models. Traffic Flow Theory: A State of the Art Report. Revised Monograph on Traffic Flow Theory. Laval, J.A., Leclercq, L., 2013. The Hamilton–Jacobi partial differential equation and the three representations of traffic flow. Transportation Research Part B 52, 17–30. Lebacque, J., Mammar, S., Haj-Salem, H., 2008. An intersection model based on the GSOM model. In: Proceedings of the 17th World Congress, The International Federation of Automatic Control, Seoul, Korea. pp. 7148–7153. Lebacque, J.-P., 1996. The Godunov scheme and what it means for first order traffic flow models. In: International Symposium on Transportation and Traffic Theory. pp. 647–677. Lebacque, J.-P., 2002. A two phase extension of the LWR model based on the boundedness of traffic acceleration. In: Transportation and Traffic Theory in the 21st Century. Proceedings of the 15th International Symposium on Transportation and Traffic Theory. Lebacque, J.-P., 2003. Two-phase bounded-acceleration traffic flow model: analytical solutions and applications. Transportation Research Record: Journal of the Transportation Research Board 1852 (1), 220–230. Lebacque, J.-P., Haj-Salem, H., Mammar, S., 2005. Second order traffic flow modeling: supply-demand analysis of the inhomogeneous Riemann problem and of boundary conditions. In: Proceedings of the 10th Euro Working Group on Transportation (EWGT) 3. Lebacque, J.-P., Khoshyaran, M.M., 2013. A variational formulation for higher order macroscopic traffic flow models of the GSOM family. Transportation Research Part B 57, 245–265. Lebacque, J.-P., Mammar, S., Salem, H.H., 2007. Generic second order traffic flow modelling. In: Transportation and Traffic Theory 2007. Papers Selected for Presentation at ISTTT17. Leclercq, L., 2007. Bounded acceleration close to fixed and moving bottlenecks. Transportation Research Part B 41 (3), 309–319. Li, J., Zhang, H., 2013. The variational formulation of a non-equilibrium traffic flow model: theory and implications. Procedia-Social and Behavioral Sciences 80, 327–340. Lighthill, M.J., Whitham, G.B., 1955. On kinematic waves II. A theory of traffic flow on long crowded roads. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 229 (1178), 317–345. Lucet, Y., 1997. Faster than the fast Legendre transform, the linear-time Legendre transform. Numerical Algorithms 16 (2), 171–185. Mammar, S., Lebacque, J.-P., Salem, H.H., 2009. Riemann problem resolution and Godunov scheme for the Aw-Rascle-Zhang model. Transportation Science 43 (4), 531–545. Mazaré, P.-E., Dehwah, A.H., Claudel, C.G., Bayen, A.M., 2011. Analytical and grid-free solutions to the Lighthill–Whitham–Richards traffic flow model. Transportation Research Part B 45 (10), 1727–1748. Moskowitz, K., Newan, L., 1963. Notes on freeway capacity. Highway Research Record (27). Newell, G.F., 1993. A simplified theory of kinematic waves in highway traffic, Part I: general theory, Part II: Queueing at freeway bottlenecks, Part III: Multidestination flows. Transportation Research Part B 27 (4), 281–313. Piccoli, B., Han, K., Friesz, T.L., Yao, T., 2012. Second Order Models and Traffic Data from Mobile Sensors. arXiv Preprint arXiv:1211.0319. Qiu, S., Abdelaziz, M., Abdellatif, F., Claudel, C.G., 2013. Exact and grid-free solutions to the Lighthill–Whitham–Richards traffic flow model with bounded acceleration for a class of fundamental diagrams. Transportation Research Part B 55, 282–306. Richards, P.I., 1956. Shock waves on the highway. Operations Research 4 (1), 42–51. Wong, S., Wong, G., 2002. An analytical shock-fitting algorithm for LWR kinematic wave model embedded with linear speed–density relationship. Transportation Research Part B 36 (8), 683–706. Work, D.B., Tossavainen, O.-P., Jacobson, Q., Bayen, A.M., 2009. Lagrangian sensing: traffic estimation with mobile devices. In: American Control Conference, 2009. ACC’09. IEEE, pp. 1536–1543. Zhang, H.M., 2002. A non-equilibrium traffic model devoid of gas-like behavior. Transportation Research Part B 36 (3), 275–290. Zhang, P., Wong, S., Dai, S., 2009. A conserved higher-order anisotropic traffic flow model: description of equilibrium and non-equilibrium flows. Transportation Research Part B 43 (5), 562–574. Zhang, P., Wu, C.-X., Wong, S., 2012. A semi-discrete model and its approach to a solution for a wide moving jam in traffic flow. Physica A: Statistical Mechanics and its Applications 391 (3), 456–463.