Properties of the LWR model with time delay
Simone Göttlich Elisa Iacomini Thomas Jung
math.NA cs.NA 35L65, 90B20, 65M06

 




[1]University of Mannheim, Department of Mathematics, 68131 Mannheim, Germany (goettlich@uni-mannheim.de, eiacomin@mail.uni-mannheim.de)
[4]Fraunhofer Institute ITWM, 67663 Kaiserslautern, Germany (thomas.jung@itwm.fraunhofer.de)


Properties of the LWR model with time delay
    Simone Göttlich[1],   Elisa Iacomini[1],   Thomas Jung[4]
     December 30, 2023
=============================================================




In this article, we investigate theoretical and numerical properties of the first-order Lighthill-Whitham-Richards (LWR) traffic flow model with time delay. Since standard results from the literature are not directly applicable to the delayed model, we mainly focus on the numerical analysis of the proposed finite difference discretization. The simulation results also show that the delay model is able to capture Stop & Go waves. 


AMS Classification. 35L65, 90B20, 65M06   	 


Keywords. Macroscopic traffic flow models, hyperbolic delay partial differential equation, numerical simulations 





§ INTRODUCTION



Nowadays traffic models have become an indispensable tool in the urban and extraurban management of vehicular traffic.
Understanding and developing an optimal transport network, with efficient
movement of traffic and minimal traffic congestions, will have a great socio- economical impact on the society. This is why  in the last decades an intensive research activity in the field of traffic flow modelling flourished.

Literature about traffic flow is quite large and many methods have been developed resorting to different approaches.
Starting from the natural idea of tracking every single vehicle, several microscopic models, based on the idea of Follow-the-Leader, grew-up for computing positions, velocities and accelerations of each car by means of systems of ordinary differential equations (ODEs) <cit.>. Other ways go from kinetic <cit.> to macroscopic fluid-dynamic and measures approaches <cit.>, focusing on averaged quantities, such as 
the traffic density and the speed of the traffic flow, by means of systems of hyperbolic partial differential equations (PDEs), in particular conservation laws. In this way we loose the detailed level of vehicles' description, indeed they become indistinguishable from each other. 
The choice of the scale of observation mainly depends on the number of the involved vehicles, the size of the network and so on.

In this paper we deal with the macroscopic scale, in particular we will focus on first order macroscopic models. The most relevant model in this framework is the LWR model, introduced by Lighthill, Whitham <cit.> and Richards <cit.> in the '50. The main idea underlying this approach is that the total mass has to be preserved, since cars can not disappear. Moreover, in this model the mean velocity is supposed to be dependent on the density, thus is closing the equation. 
On the other hand the lacks of the LWR model are well-known. For example, it fails to generate capacity drop, hysteresis, relaxation, platoon diffusion, or spontaneous congestions like Stop & Go waves, that are typical features of traffic dynamics.
These drawbacks are due to the fact that the LWR model represents a simplification of the reality, assuming that accelerations are instantaneous and traffic is described only at the equilibrium.

In order to overcome these issues, second order models have been proposed, see <cit.>. They take into account
the non-equilibria states, assuming that accelerations are not instantaneous. To do this, the equation that describes the variation of the velocity in time has to be added to the system, replacing the typical given law of the fist order models.
Other ways are also possible to improve first order models, just considering phase transition models <cit.> or multi-scale approaches <cit.>.

Instead of switching to second order models, we propose a first order macroscopic model with a time delay term in the flux function, for taking into account that the velocity can not change instantaneously. In this framework the delay represents the reaction time of both drivers and vehicles.

At a microscopic level, a model with time delay appears for the first time in the work done by Newell <cit.>, then similar models are presented in <cit.>. The mathematical tools needed in this framework are not systems of ODEs anymore, but systems of delay differential equations (DDEs), particular differential equations in which the derivative of the unknown function at a certain time is given in terms of the values of the function at previous times.
Macroscopic models can be derived from microscopic description following a well-known procedure described in <cit.>.
Depending on how to treat the delay term, one can recover different macroscopic models, as in <cit.>, in which a Taylor's approximation is applied to the delay term and the obtained model is a diffusive LWR type model. 
On the other hand, we want to keep the delay in the explicit form, and therefore avoid the diffusion approximation. The model derived in <cit.> will be studied in details in the following, investigating carefully its theoretical and numerical features. 

Several delayed-systems are presented in literature, since many phenomena need some transient to become visible or effective: the study of the evolution of the HIV in medicine <cit.>, cell population dynamics in biology <cit.>, the feedback control loops in control engineering <cit.>, and many applications in mechanics and economics <cit.>, but to the authors best knowledge, they are closer to delayed parabolic partial differential equations or to delayed ordinary differential equations, i.e. they are studied only at a microscopic level.

In this work instead we deal with a delayed hyperbolic partial differential equation. We will point out similarities and differences with the undelayed model in order to catch the effect of the delay on traffic dynamics, both from theoretical and numerical points of view.
Moreover, since we are interested in reproducing real traffic phenomena, the numerical tests are mainly focused on traffic instabilities.
In particular we investigate the phenomenon of Stop & Go (S&G) waves, which are a typical feature of congested traffic and represent a real danger for drivers. 
They lead not only to safety hazard, but they also have a negative impact on fuel consumption and pollution.
Indeed a S&G wave is detected when vehicles stop and restart without any apparent reason, generating a wave that travels backward with respect to the cars' trajectories.


Since modeling properly this phenomenon is crucial for developing techniques
aimed at reducing it, a considerable literature is growing up on this topic. This means that a lot of models have been developed in the last years, i.e. <cit.>, 
and also several real experiments took place, just see <cit.>.

In this framework, our aim is to investigate if our delayed model is able to capture the S&G phenomena and, therefore, to present an easy to use algorithm able to reproduce S&G waves at a macroscopic level. 
Indeed from the numerical point of view, just an altered Lax Friedrichs method will be employed to compute the evolution of the density.
In order to validate our model, several numerical tests will be provided for
comparing our delayed model with the existing ones.




  §.§.§ Paper organization

In Section 2, we introduce the delayed model and investigate its theoretical properties, as the conservation of mass, the positivity and the boundedness of the solution. After that, we focus on the numerical aspects, presented in Section 3, proposing a suitable numerical scheme and checking the theoretical features still hold. Section 4 is completely devoted to the numerical tests.





§ THE DELAYED TRAFFIC FLOW MODEL





In macroscopic models <cit.>, traffic is described in terms of macroscopic variables such as density ρ=ρ(x,t), that is the number of vehicles in a kilometre, and the mean velocity V=V(x,t) at the point x∈ at time t>0. 

The LWR model, introduced by Lighthill, Whitham <cit.> and Richards <cit.>, is one of
the oldest and still most relevant first order macroscopic models for traffic flow.
The natural assumption that the total mass is conserved along the road is closed by the assumption that the velocity V=V(ρ) is given as function of the density ρ:

    ∂_tρ(x,t)+∂_x(ρ(x,t)  V(ρ(x,t)))=0
    ρ(x,0)=ρ^0(x).


A lot of possible choices for the function V(ρ) are available in the literature, i.e. the Greenshields function <cit.> which proposes a linear velocity function:

    V(ρ)=ρ_max( 1- ρ/ρ_max).


In order to simplify the notation, we will consider the normalized quantities  ρ_max=V_max=1.









Aiming to overcome the drawbacks of LWR model presented in the introduction, we propose a first order macroscopic model with a time delay term in the flux function, for taking into account that the velocity can not change instantaneously.

In this framework the delay represents the reaction time of both drivers and vehicles.
Such a model has been recovered from a delayed microscopic model, as shown in <cit.> keeping the delay in the explicit form. Assuming T≥0 as the time delay, we consider:

    ∂_t ρ(x,t) + ∂_x (ρ(x,t)  V(ρ(x,t-T)))=0.

We will call this model delayed LWR model. Note that in the limit case of T=0 the classical LWR model is recovered and therefore, it can be interpreted as a generalization of the LWR model. On the other hand,
if the delay is too large and there are suitable initial conditions, cars can overtake or crash each other, for example when a vehicle suddenly
brakes and the following car is not reacting in time to slow down.

Note that in order to guarantee the well-posedness of the problem, we have to provide an initial history function as initial data defined on [-T,0], thus we need ρ^0(x,t) defined on t∈ [-T,0] when starting at t=0.



 §.§ Properties of the model

After introducing the delayed model, we want to investigate its properties. Since this model can be seen as a generalization of the classical LWR model, i.e. when T=0, it is natural to investigate how its properties differ from
the undelayed model.



  §.§.§ Conservation of mass

In  the  framework  of  conservation  laws  and  traffic flow models the conservation of the total mass is a crucial property which has to be guaranteed. 
For the LWR model (<ref>), we have one equation and one
conserved quantity, i.e. ρ. Indeed cars do not appear or disappear, they can only enter and leave at the boundaries.
We note that the structure of the equation in the delayed and undelayed model stays
the same and that we have a flux function that is now dependent on two variables.


The delayed LWR model (<ref>) conserves the quantity ρ(x,t).


	We integrate the equation (<ref>) over an arbitrary space interval [a,b] and get
	
    d/dt∫_a^bρ(x,t)dx   =-∫_a^b∂_x(ρ(x,t)V(ρ(x,t-T)))dx
       =ρ(a,t)V(ρ(a,t-T))-ρ(b,t)V(ρ(b,t-T)).

	Since ∫_a^bρ(x,t)dx is the amount of density in the interval [a,b], d/dt∫_a^bρ(x,t)dx denotes the change over time for the density. Therefore, the density only changes due to the flux at the boundaries a and b for every space interval.


We see that the density is still conserved in the delayed model, which is very important for its reliability. The introduction of an explicit time delay therefore does not destroy this property.



  §.§.§ Positivity

Another property one would ensure is the positivity of the solution. Indeed we want the density to stay positive, as negative densities have no physical meaning. 
 
	Assume we have initial data with non-negative density ρ. For the delayed LWR model, then the density stays non-negative.


	We rewrite (<ref>) as
	
    ∂_tρ(x,t)=-(ρ(x,t)∂_xV(ρ(x,t-T))+V(ρ(x,t-T))∂_xρ(x,t)).

	For the density to become negative, we need to have ρ=0 and ∂_tρ<0. Hence, assume ρ(x^*,t^*)=0 at an arbitrary point (x^*,t^*). This means, due to our assumption, that ρ(x^*,t^*) is a minimum, since the density is positive for t≤ t^* and therefore, if the derivative exists, ∂_xρ(x^*,t^*)=0. Plugging in then gives us
	
    ∂_tρ(x^*,t^*)=-(0 ∂_xV(ρ(x^*,t^*-T))+V(ρ(x^*,t^*-T)) 0)=0.

	We have therefore shown that ρ can not become negative.


	The velocity V in this model is a function of ρ and can be chosen and altered depending on the needs. The properties regarding the velocity in the first order model can therefore be acquired by choosing a suitable function V.  For example, we can have lower and upper bounds for the velocity by defining V to be cut at the bounds.




  §.§.§ Upper bound


The last property we want to investigate is the boundedness of the solution. In particular, we want to know if there is a maximal density. For the undelayed model, this is guaranteed. For the delayed model, we need to check if this still true. 

	Assume V is monotone decreasing and V(ρ_max)=0 for ρ_max, the maximal density in the classical LWR model. The delayed first order model (<ref>) has no maximal density ρ_max.


	Assume we have a maximal density ρ_max. The velocity function V is chosen in such a way that V(ρ_max)=0 and monotone decreasing. Then, for an arbitrary point (x^*,t^*) where ρ(x^*,t^*)=ρ_max we have
	
    ∂_t ρ(x^*,t^*)=-ρ(x^*,t^*)∂_xV(ρ(x^*,t^*-T))-V(ρ(x^*,t^*-T))∂_xρ(x^*,t^*).

	Since ρ(x^*,t^*) is the maximal density, ∂_xρ(x^*,t^*)=0 if the derivative exists and we have
	
    ∂_t ρ(x^*,t^*)=-ρ(x^*,t^*)∂_xV(ρ(x^*,t^*-T))

	left. We know ρ(x^*,t^*)≥0 and this means that the sign of ∂_tρ(x^*,t^*) is only dependent on ∂_xV(ρ(x^*,t^*-T)).
	
	In the undelayed case, we know that ∂_xV(ρ(x^*,t^*))>0, since V is monotone decreasing and ρ(x^*,t^*) is the maximal ρ. 
	
	In the delayed case, we do not have knowledge if ρ(x^*,t^*-T) is maximal, so we can in general say nothing about ∂_xV(ρ(x^*,t^*-T)). This means, in general, ρ>ρ_max is possible.


	Regarding the positivity, we claimed that the choice of V in the first order model is a key to guarantee a positive velocity. We here see, due to the fact that ρ overshoots any ρ_max, that the classical choices for V must be altered to avoid negative velocities, i.e. we need to cut the function.


If the density ρ>ρ_max, the model is not reliable any more. On the other hand this situation could not be avoided since rear-end collisions are actually possible in real situations.






§ NUMERICAL DISCRETIZATION


After the investigations on the analytical properties of the delayed LWR model, let us focus on its numerical counterpart.

Since (<ref>) is a hyperbolic partial differential equation, we can employ the Lax-Friedrichs method for the numerical approximation.
To do that, we first introduce space and time steps Δ x, Δ t >0  and a grid in space {x_i= i  , i ∈ℤ} and time {t^n= n  , n ∈ℕ}.
Discretized variables are expressed by ρ_i^n, where i is the space and n the time index.  Also, we have Δ t ≤ T to be able to treat the delay.

The Lax-Friedrichs method for (<ref>) is stated by:

    ρ_i^n+1=1/2(ρ_i+1^n+ρ_i-1^n)-Δ t/2Δ x(f(ρ_i+1^n)-f(ρ_i-1^n)).

Now using the structure of (<ref>), we can identify a flux function f(ρ(x,t-T),ρ(x,t))=V(ρ(x,t-T))ρ(x,t). Plugging this into the Lax-Friedrichs method, we end up with an altered Lax-Friedrichs method

    ρ_i^n+1=1/2(ρ_i+1^n+ρ_i-1^n)-Δ t/2Δ x(f(ρ_i+1^n-T_Δ,ρ_i+1^n)-f(ρ_i-1^n-T_Δ,ρ_i-1^n)),

where T_Δ is the number of steps that make up the time delay T.
In order to guarantee the well-posedness of the discrete problem, we have to provide an initial history function as initial data defined on [-T,0], as we said above for the continuous problem. The simplest choice one can do is to consider ρ^0(x,t) as a constant function on t∈ [-T,0] when starting at t=0. In the following we will assume that ρ^0(x,0) is constant in t for t∈ [-T,0].
The Lax-Friedrichs method has a CFL condition in the classical case, which is given as Δ t≤Δ x/max_k(λ_k), where λ_k are the eigenvalues of the jacobian matrix of f. We also expect to find a CFL condition in the delayed case, but a priori it is not clear how this condition may look like. In the following, we want to investigate some properties of this method, and in this process we will find an appropriate CFL condition. For the sake of the calculations, we assume the velocity function V to be the Greenshields function, or a cut variation of it, where we have |V(ρ)|≤|ρ_max| with V_max=1.



 §.§ Properties of the discretization




  §.§.§ Conservation of mass


First, we check if the conservation property is preserved from the numerical scheme. Here, we assume the density to be on a compact support, so we do not have infinite density initially.
We get

    Δ x∑_j ρ^n+1_j=Δ x∑_j 1/2(ρ_j+1^n+ρ_j-1^n)-Δ t/2Δ x(V(ρ_j+1^n-T_Δ)ρ_j+1^n-V(ρ_j-1^n-T_Δ)ρ_j-1^n),

where the part 
    ∑_j 1/2(ρ_j+1^n+ρ_j-1^n)=∑_j ρ_j^n
 and the part 
    ∑_j -Δ t/2Δ x(V(ρ_j+1^n-T_Δ)ρ_j+1^n-V(ρ_j-1^n-T_Δ)ρ_j-1^n)
 is a telescope sum and equals zero due to the compact support. This gives us 
    Δ x∑_j ρ^n+1_j=Δ x∑_j ρ^n_j
 and therefore conservation.



  §.§.§ Positivity


We show that Lemma <ref> holds also at the discrete level under a certain CFL condition. Starting with (<ref>), we see
that 1/2(ρ_j+1^n+ρ_j-1^n) is always positive, since we assume ρ^n to be positive. W.l.o.g. we can even say it is bigger than 0, since if it is 0, ρ_j^n+1 will also be zero. So to guarantee positivity, we need to guarantee

    Δ t/2Δ x(V(ρ_j+1^n-T_Δ)ρ_j+1^n-V(ρ_j-1^n-T_Δ)ρ_j-1^n)≤1/2(ρ_j+1^n+ρ_j-1^n).

We introduce a CFL-condition, namely Δ t_n≤Δ x/max{|ρ^n|,|ρ^n-T_Δ|}. Therefore, we get for the left-hand-side

    1/2max{|ρ^n|,|ρ^n-T_Δ|}(V(ρ_j+1^n-T_Δ)ρ_j+1^n-V(ρ_j-1^n-T_Δ)ρ_j-1^n)
    ≤   V(ρ_j+1^n-T_Δ)/2max{|ρ^n|,|ρ^n-T_Δ|}ρ_j+1^n-V(ρ_j-1^n-T_Δ)/2max{|ρ^n|,|ρ^n-T_Δ|}ρ_j-1^n
    ≤   max{|ρ^n|,|ρ^n-T_Δ|}/2max{|ρ^n|,|ρ^n-T_Δ|}ρ_j+1^n+max{|ρ^n|,|ρ^n-T_Δ|}/2max{|ρ^n|,|ρ^n-T_Δ|}ρ_j-1^n
    
    =   1/2(ρ_j+1^n+ρ_j-1^n),

which shows the positivity for this CFL-condition. Here we use that |V(ρ)|≤max{|ρ^n|,|ρ^n-T_Δ|} if ρ is positive.



  §.§.§ L^∞-Bound


Focusing on the boundedness of the discrete solution, we look for an estimate in the norm |·|_L^∞. Since we have positivity, only an upper bound for ρ in (<ref>) is required. We assume that the data at time t_n has an upper bound which we denote with ρ^n_max. We further denote V(ρ_j^n-T_Δ)ρ_j^n=f(ρ_j^n-T_Δ,ρ_j^n) and get for a special ξ=(ξ_1,ξ_2) using the mean value theorem:

    ρ_j^n+1   =1/2(ρ_j+1^n+ρ_j-1^n)-Δ t/2Δ x(f(ρ_j+1^n-T_Δ,ρ_j+1^n)-f(ρ_j-1^n-T_Δ,ρ_j-1)^n)
       = 1/2(ρ_j+1^n+ρ_j-1^n)+Δ t/2Δ x(∇ f(ξ_1,ξ_2)[         ρ_j-1^n-ρ_j+1^n; ρ_j-1^n-T_Δ-ρ_j+1^n-T_Δ ])
       =1/2(ρ_j+1^n+ρ_j-1^n)+Δ t/2Δ x(V(ξ_2)(ρ_j-1^n-ρ_j+1^n)+ξ_1(ρ_j+1^n-T_Δ-ρ_j-1^n-T_Δ))
       ≤1/2(ρ_j+1^n+ρ_j-1^n)+1/2(ρ_j-1^n-ρ_j-1^n-T_Δ+ρ_j+1^n-T_Δ-ρ_j+1^n)
       =1/2(2ρ_j-1^n-ρ_j-1^n-T_Δ+ρ_j+1^n-T_Δ)
       ≤ρ_max^n+1/2(ρ_j+1^n-T_Δ-ρ_j-1^n-T_Δ)
       ≤3/2max{|ρ^n|,|ρ^n-T_Δ|}


We use again the positivity and the CFL-condition.
We can find a different estimate for ρ^n+1_j, depending on T. Note that the estimate in (<ref>) is more accurate. 
We start in the same way as above, and have 

    ρ_j^n+1≤1/2(2ρ_j-1^n-ρ_j-1^n-T_Δ+ρ_j+1^n-T_Δ).
 Rearranging gives us

    1/2(ρ_j-1^n+ρ_j-1^n-ρ_j-1^n-T_Δ+ρ_j+1^n-T_Δ)=1/2(ρ_j-1^n-ρ_j-1^n-T_Δ+ρ_j+1^n-T_Δ+ρ_j-1^n),

and by again using the mean value theorem for ρ_j-1(t) we get

    ρ_j^n+1≤1/2(T∂_tρ_j-1+ρ_j-1^n+ρ_j+1^n-T_Δ)≤max{|ρ^n|,|ρ^n-T_Δ|}+1/2T||∂_tρ||_L^∞







  §.§.§ TV Bound


Let us look for an estimate on the Total Variation for the method (<ref>). The first thing we check is the difference between the velocity function in two cells.

    V(ρ_j+1^n-T_Δ)-V(ρ_j^n-T_Δ)≤(1-ρ_j+1^n-T_Δ)-(1-ρ_j^n-T_Δ)=ρ_j^n-T_Δ-ρ_j+1^n-T_Δ=-Δ_j+1/2^n-T_Δ

Then, we look at the difference between two neighboring cells, where we use (<ref>), (<ref>) and the notation Δ_j+1/2^n=ρ_j+1^n-ρ_j^n:

    Δ_j+1/2^n+1   =ρ_j+1^n+1-ρ_j^n+1
       =1/2(Δ_j+3/2^n+Δ_j-1/2^n)
       -Δ t/2Δ x(V(ρ_j+2^n-T_Δ)Δ_j+3/2^n-V(ρ_j^n-T_Δ)Δ_j-1/2^n-ρ_j+1^nΔ_j+3/2^n-T_Δ+ρ_j-1^nΔ_j-1/2^n-T_Δ)

The Total Variation at time t_n+1 (denoted by TV(ρ^n+1_Δ)) is given as ∑_j|Δ_j+1/2^n|. For the next step and with (<ref>) as well as shifting the indices, we get 

    ∑_j|Δ_j+1/2^n+1|   =∑_j|ρ_j+1^n+1-ρ_j^n+1|
    
    =∑_j(1   +2|Δ t/2Δ xρ_j+1^n-T_Δ|+2|Δ t/2Δ xρ_j+1^n|)max{|Δ_j+1/2^n|,|Δ_j+1/2^n-T_Δ|}.

If we introduce a CFL-condition of the type Δ t≤2Δ x/max{|ρ_j^n-T_Δ|,|ρ_j^n|} (i.e. not as strong as above), the dependency on j disappears and we get

    (5+1/max{|ρ_j^n-T_Δ|,|ρ_j^n|})∑_jmax{|Δ_j+1/2^n|,|Δ_j+1/2^n-T_Δ|}.

By assumption, TV(ρ_Δ^n) and TV(ρ_Δ^n-T_Δ) are finite, so we can estimate further

    TV(ρ_Δ^n+1)≤2(5+1/max{|ρ_j^n-T_Δ|,|ρ_j^n|})max{TV(ρ_Δ^n),TV(ρ_Δ^n-T_Δ)}.
 
For the estimate in time, we look at 

    ∑_j|ρ_j^n+1-ρ_j^n|.

By simply plugging in (<ref>) and using the same CFL-condition, we can write

    ∑_j|ρ_j^n+1-ρ_j^n|   ≤∑_j1/2(|ρ_j+1^n|+|ρ_j-1^n|)+Δ t/2Δ x(|V(ρ_j+1^n-T_Δ)||ρ_j+1^n|+|V(ρ_j-1^n-T_Δ)||ρ_j-1^n|)+|ρ_j^n|
       ≤∑_j2max{|ρ_j^n-T_Δ|,|ρ_j^n|}+2+2max{|ρ_j^n-T_Δ|,|ρ_j^n|}
       =∑_j4max{|ρ_j^n-T_Δ|,|ρ_j^n|}+2.

We now have a BV-Bound in space as well as in time, which gives us all the desired  BV estimates. 



  §.§.§ Time span


With the CFL-condition we introduced, the time step can become smaller every step, since we have no maximum principle. Here, we want to see if we can actually reach every time horizon. Therefore, we plug our L^∞-bound into the CFL and get

    Δ t≤Δ x/max{|ρ^n|,|ρ^n-T_Δ|}≤Δ x/(3/2)^n||ρ^0||_L^∞.

Now, the time horizon we reach with n time steps is given by

    t_n=∑_i=1^nΔ t=Δ x/||ρ^0||_L^∞∑_i=1^n(2/3)^i.

So for infinite time steps n→∞, 

we end up with a geometric series which converges to 3Δ x/||ρ^0||_L^∞. That is the time horizon we can guarantee with this estimate. With the alternative estimate that depends on T, the time horizon we can reach is dependent also on T. With n time steps, we can now reach the time horizon

    t_n=Δ x∑_i=1^n1/ρ^0_max+i1/2||∂_tρ||T.

This is basically a harmonic series shifted and with a factor, but it is divergent to ∞ if the factor is not zero. This also means, that for small T, the steps can be larger. Furthermore, with this second estimate, we can guarantee to reach every time horizon.



§ NUMERICAL RESULTS

This section is devoted to the numerical simulation results for the model presented above, focusing in particular on the S&G waves phenomenon, a typical feature of congested traffic, detected when vehicles stop and restart without any apparent reason, generating a wave that travels backward with respect to the cars' trajectories.



Starting from empirical observations and the work done in <cit.>, let us assume the velocity function as follow: 

    V(ρ)= 
    V_max       ρ≤ρ_f 
    α ( 1/ρ - 1/ρ_c )         ρ_f<ρ<ρ_c  
    
    0        ρ≥ρ_c

where α>0 is a parameter, ρ_c∈ (0,ρ_max] and ρ_f ∈ [0,ρ_max) are two density thresholds. In particular ρ_c represents the so-called safe distance at the macroscopic level: if the density is higher than ρ_c, vehicles do not respect the safe distance so the desired velocity has to be 0, indeed they should stop. 

On the other hand, if the density is very low, which means that vehicles are far enough from each others, the desired velocity is the maximum one.

Note that (<ref>) respects the hypothesis |V(ρ)|≤ |ρ_max| choosing V_max=1. Moreover, depending on the choice made for α, the velocity function can be discontinuous.

For the discretization, let us assume the space interval [a,b]=[0,1], =0.02 and periodic boundary conditions. Moreover, the time step  is chosen in such a way the CFL condition is satisfied.
The density thresholds are ρ_c=0.75 and ρ_f=0.2 as real data suggests.
 The delay term depends on the CFL condition and the initial data, for this reason each numerical test has its time delay interval which ensures the reliability of the model, i.e. ρ≤ρ_max. However, in general, one has to assume T_Δ one order of magnitude greater than  to see the effect, i.e. T_Δ≈ (10,20).



 §.§ Backward propagation of Stop & Go waves

In order to be more comprehensive as possible in reproducing S&G waves, let us describe first the backward propagation of the perturbation. After that we will focus also on the triggering of this phenomenon. 




  §.§.§ Test 0


In order to point out the crucial role played by the delay in this framework, let us compare the evolution of the density obtained with the delayed model and the classical LWR model, or, in other words, when T_Δ=0.

Assume as initial data ρ^0(x)=5/8+1/8sin(2π x). Moreover the time step is =0.01 and the delay T_Δ=15 in the delayed case.



It is evident how the LWR model smears out the perturbations in the initial data and after a certain time the density becomes constant on the whole road, see Fig. <ref>(right). On the other hand, the delayed model preserves the perturbations and also makes them increase as usually happens in traffic evolution, Fig. <ref>(left).


 We have to be very careful in choosing the delay term. Indeed if the delay is too small, we recover a situation very similar to the LWR model, but, on the other hand, if the delay is too high, i.e. T_Δ=18 , the hypothesis on the model are no longer satisfied and the density grows more than 1, so the model has no sense anymore, see Fig. <ref>. 


In the following, we will consider the delay as the maximum allowed by the model feasibility. 





  §.§.§ Test 1


In this numerical simulation we want to reproduce with our model the tests presented in <cit.>. Starting from the same initial data our aim is to recover a similar behaviour for the density. In <cit.>, a nonlinear 2-equations discrete velocity model is implemented to compute the density evolution:

    ∂_t ρ + ∂_x q    =0
    ∂_t q + Hq/1-ρ∂_x ρ + (1-Hq/1-ρ∂_x q  )    =-1/ε(q-F(ρ)),

where q is the flux, F is the fundamental diagram and H is a measure for the look ahead and the nonlocality of the equations. Such a model converges to the LWR type equations in the relaxation limit but shows also similarities to the Aw-Rascle model. The main difference to the delayed model is that at the macroscopic level they are composed of a system of two PDEs while the delayed model is described only by a conservation law with a time delay in the velocity term.

Let us assume ρ^0(x,k)=5/8+1/8sin(2kπ x) for k=1,2 and =0.01. We consider T_Δ=16 and T_Δ=22 , respectively if k=1,2.



In Fig. <ref>, the density values on the (x,t)-plane are represented by the colors as shown in the colorbar. We obtain persistent waves where the number of waves is directly related to the perturbations in the initial condition, as in <cit.>.
Numerically, we recover this result if T_Δ∈ (11, 16) for k=1 and T_Δ∈ (18, 22) for k=2.



  §.§.§ Test 2

 We consider the discrete delayed model (<ref>) with the velocity function (<ref>) and initial data:
 
    ρ^0(x)=
         0.6         x< 0.5
    
         0.1         x≥ 0.5.

 Let us assume the time delay as T_Δ=10 and the time step =0.01 such that the CFL condition is satisfied.
 Moreover, the coefficient α is chosen in such a way as to make the velocity function continuous.
 







Looking at Fig. <ref>(left), one can immediately note that the initial slowdown in the first half of the road is increasing in time until vehicles completely stop since ρ≥ρ_c, as the colorbar suggests. This is a typical S&G wave  behaviour, as we can see also in <cit.>, where they compute the density evolution considering another extension of the LWR model. Indeed they assume that the velocity term depends on ρ, the reaction time τ and the derivative in space of the velocity itself ∂_x V(ρ):


    ∂_tρ + ∂_x ( ρ V (ρ/1-τ∂_x V(ρ)))=0.

 
Investigating the density profile at the end of the simulation, we are able to recognize a well-defined S&G wave profile, see Fig. <ref>(right), as described in <cit.>.



Moreover, let us note that if the delay is smaller, i.e. T_Δ=4, the model is no longer able to preserve the perturbation and the density profile becomes smoother, see Fig. <ref>.



Choosing the delay term 7  <T_Δ< 11, we are able to recover the S&G behaviour.



 §.§ Triggering of Stop & Go waves

Let us now focus on the triggering of S&G waves. 
Starting from a small perturbation, i.e. a slowdown, in the initial data, our aim is to find out if it is possible to recover a S&G wave. In this direction, let us consider the example presented in <cit.>, in which the initial data is given by:

    ρ^0(x)=
        0.35         1.34≤ x≤ 1.342
    
        0.2         elsewhere,

and Dirichlet boundary conditions.
The initial data is modelling a small region, a cell, where vehicles are moving slower than elsewhere and therefore the density is higher in that cell. This slowdown can be caused by the presence of sags, road sections in which gradient changes significantly from downwards to upwards <cit.>, or the presence of a school zone in which the velocity has to be reduced, <cit.>.

In <cit.>, the density evolution is computed by coupling the LWR model with a second order microscopic model, specifically conceived to reproduce S&G waves.
Instead of switching to multiscale models in which we have to manage the microscopic model too, let us see if we can recover the density evolution with the delayed model.

Assuming =0.009 and T_Δ=21, we note that the initial perturbation increases and moves backward as the time increases, as it happens in the multiscale case. 


Therefore, the delayed model is able to reproduce the triggering of S&G waves too, not only their backward propagation.



§ CONCLUSION

In this paper, we have introduced theoretical and numerical properties of the delayed LWR traffic model. While the derivation of the model has been done in <cit.>, the numerical behavior of the delayed model has not been studied intensively before. 

Starting from the undelayed scenario, we investigated the theoretical features of the delayed model to point out the numerical properties of the scheme and we proposed an altered Lax-Friedrichs method to compute the evolution of the density.
The key observation therefore is that the delayed model is really able to reproduce Stop & Go waves for the right choice of parameters. Comparisons to already existing results from the literature also underline this characteristic. 

Future works will include the extension to networks as well as parameter estimation techniques to determine the time delay.

  
abbrv

