Mathematical Modeling and Simulations

of Traffic Flow

by
Jiah Song

A dissertation submitted in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy
(Applied and Interdisciplinary Mathematics)
in The University of Michigan
2019

Doctoral Committee:
Professor Smadar Karni, Co-Chair
Professor Romesh Saigal, Co-Chair
Professor Robert Krasny
Professor Philip Roe
 Jiah Song

[email protected]

ORCID iD: 0000-0003-4317-2286

c Jiah Song 2019

To Leah and Seonghoon with all my love.

ii
 ACKNOWLEDGEMENTS

First, I want to express my sincere gratitude to my advisors Professors Smadar

Karni and Romesh Saigal. I am especially thankful for Professor Karni’s patient

guidance, constant support, and kind encouragement. She is a role model of me who

can be incredibly patient, brilliant and successful. Professor Saigal is an incredible

mentor who provided inspired and valuable advice whenever I seek for a help. I

especially thank for his help in finding a job. I was very lucky to have them as my

advisors who are immensely approachable. We had so many productive conversation

about the research.

I would like to thank my committee members Professors Philip Roe and Robert

Krasny for taking their time to allow me to present my research, and offering useful

comments to this dissertation.

I also thank all the friends that I have made here in Michigan, and all the support

from the Math department. Thanks to my family and friends in South Korea who

always believe in me and support me with your unconditional love.

Lastly, Leah and Seonghoon. Leah, you gave me a new life. Thank you for making

me a good person, my forever sunshine. All the work in this dissertation went along

with you for the past three years. Seonghoon, this dissertation would have not been

able to be complete without your support. I look forward to spending our time after

graduation. You are my only hero.

iii
 TABLE OF CONTENTS

DEDICATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii

ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

CHAPTER

I. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Hyperbolic Conservation Laws . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Mathematical Theory . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.2 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2 Traffic Field Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

II. Microscopic Car-Following Model . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2 Numerical Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

III. Macroscopic Traffic Flow Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.1 Scalar Conservation Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.1.1 Derivation of Equation . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.1.2 Fundamental Relation and a New Fundamental Diagram . . . . . . 26

iv
 3.1.3 Riemann Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.1.4 Generalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.2 System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.2.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.2.2 A New Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.3 New Multilane Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.3.1 Lane Change Conditions . . . . . . . . . . . . . . . . . . . . . . . . 46

3.3.2 Multilane Scalar and System Models . . . . . . . . . . . . . . . . . 47

IV. Numerical Method and Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.1 Numerical Method Revisit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.2 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

V. Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

v
 LIST OF FIGURES

Figure

1.1 Shock jump condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2 Example I.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.3 Example I.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.4 Schematic for numerical computation . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.5 Study area schematic and camera coverage . . . . . . . . . . . . . . . . . . . . . . . 11

1.6 Two velocity visual representations of US 101 (lane 1) . . . . . . . . . . . . . . . . 11

1.7 Two velocity visual representations of I-80 (lane 2) . . . . . . . . . . . . . . . . . . 12

2.1 Car-following schemetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.2 Uniform flow (a), disturbed flow with one entering car (b) and leaving car (c). . . . 20

2.3 Top: car trajectories with #1 car trajectory (bold line). Bottom: speed of #1 car. 20

3.1 The Greenshields’ model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.2 Example III.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.3 Example III.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.4 Some visual representations of Table 3.1. . . . . . . . . . . . . . . . . . . . . . . . . 26

3.5 US 101 data represented in v − ρ plane, original in black and smoothed data in

green. Fitted curve is in blue (Left). Same data in f − ρ plane (Right). . . . . . . 28

3.6 Reconstruction and errors of US 101 . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.7 I-80 data represented in v − ρ plane, original in black and smoothed data in green.

Fitted curve is in blue (Left). Same data in f − ρ plane (Right). . . . . . . . . . . 29

3.8 Reconstruction and errors of I-80 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.9 A practical shape of the new FD (3.7) . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.10 Road example of the variable speed limit . . . . . . . . . . . . . . . . . . . . . . . . 33

3.11 Riemann solution structure in terms of wave propagation . . . . . . . . . . . . . . 33

3.12 Graphs of admissible fL and fR in solid curves . . . . . . . . . . . . . . . . . . . . 33

vi
3.13 Riemann solution locus and evolution of density. . . . . . . . . . . . . . . . . . . . 35

3.14 Convex hull, evolution of density, and characteristics of Example III.5 . . . . . . . 36

3.15 Convex hull, evolution of density, and characteristics of Example III.6 . . . . . . . 37

3.16 Headway in microscopic (a) and macroscopic (b) description . . . . . . . . . . . . . 43

3.17 Numerical study on diffusive effect . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.18 Diagram illustrating lane changes in three lanes . . . . . . . . . . . . . . . . . . . . 48

4.1 Instability of PDE system of (3.26) . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4.2 Example IV.2 results by (3.26) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.3 Example IV.2 results by (3.26) with three-phase V (ρ) . . . . . . . . . . . . . . . . 55

4.4 Example IV.3 results by (3.31) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.5 Example IV.3 results by (3.31) with (3.18) . . . . . . . . . . . . . . . . . . . . . . . 59

4.6 Example IV.4 results by (3.33). In this simulation, R = 0. . . . . . . . . . . . . . . 61

4.7 Example IV.4 results by (3.33). In this simulation, R = 0.3. . . . . . . . . . . . . . 62

4.8 Example IV.4 results in lane 1 superimposed with and without effect of lane changes. 63

4.9 Trajectories of vehicles in lane 1 according to the computation of Example IV.4 . 63

A.1 Optimized m and α, and its relation with ρ . . . . . . . . . . . . . . . . . . . . . . 67

vii
 ABSTRACT

This dissertation concerns modeling and computation of macroscopic fluid-like traffic

flow on a single lane road as well as on multilane, coupled with lane changes. We

start with single lane models and consider the scalar model, expressing car mass

conservation, proposed in the 1950s (Lighthill-Whitham 1955, Richards 1956). This

requires a velocity-density closure relation, called Fundamental Diagram. We pro-

pose a more realistic closure relation, based on field data fitting. Using a relationship

between car headway and flow density, we develop and extend to a new 2x2 system

that includes an intra-lane acceleration equation reflecting microscopic flow charac-

teristics. The models are then generalized to multilane traffic flow and incorporate

mass exchange terms, as well as inter-lane acceleration terms. We present two types

of multilane models that are nonlinear hyperbolic equations: (i) scalar models and

(ii) 2x2 systems. Both types are formulated by establishing lane changing conditions

from drivers’ perspectives and incorporating them into source terms. Consideration

for lane changing includes the potential for velocity gain and available space. Using

Roe-type upwind scheme, we reproduce realistic flow patterns such as stop-and-go

flow and study the influence of lane changing and lane reduction on flow capacity.

viii
 CHAPTER I

Introduction

The scope of this dissertation encompasses mathematical models, numerical com-

putations of traffic flow, and understanding its essential fundamentals based on field

data. Chapter I presents the mathematical structure and its numerical aspects in

section 1.1 and the description of data being utilized in section 1.2.

Two classes of models are taken into account in describing traffic flow: micro-

scopic and macroscopic descriptions. The former is based on individual vehicles’

motion, and the latter is concerned with average behavior. Our interest is mainly in

macroscopic understanding that naturally leads to the hyperbolic partial differential

equations (PDEs). The following, therefore, sheds light on theory and numerical

methods of nonlinear hyperbolic equations.

1.1 Hyperbolic Conservation Laws

1.1.1 Mathematical Theory

We start by stating a system of hyperbolic conservation laws with an initial con-

dition

 ∂t W + ∂x F (W ) =

0
(1.1)
W (x, 0) = W0 (x)



1
 2

where t is the time, x is the space, W = W (x, t) is a vector of n conserved variables,

and F (W ) is a vector of n fluxes. For brevity, x and t are omitted. Here

W = (w1 , w2 , · · · , wn )T ,

and

F (W ) = (f1 (W ), f2 (W ), · · · , fn (W ))T .

We deal with one-dimensional case in x but note that this may be extended to multi-

dimensions as needed. Many practical problems, including gas dynamics, a spiral

galaxy, shallow water waves, traffic flow etc, involve conserved quantities. Therefore,

the field of hyperbolic conservation laws provides the theoretical underpinnings for

modeling them.

In a quasilinear form, (1.1) can be put into the following frame:


 ∂t W + A(W )∂x W =

0
(1.2)
W (x, 0) = W0 (x)



where  
 A11 ... A1n 
∂F (W )  
A(W ) = F 0 (W ) = =  . . . Aij = ∂fi

... 
∂W  ∂wj 
 
An1 ... Ann

is the Jacobian matrix of F (W ). (1.1) is linear if A(W ) returns a constant matrix, or

nonlinear otherwise. (1.1) is hyperbolic if A(W ) has n real eigenvalues and linearly

independent n eigenvectors. Derivation of hyperbolic systems can be found in many

books [32, 35, 51, 53]. We will later focus on the derivation of scalar equation in

terms of traffic flow in section 3.1.1.

Hyperbolic PDE features characteristic curves along which information propagates

 3

at a finite speed. Let us take the Burgers’ equation for illustrating characteristics

 ∂t w + ∂x f (w) =

0
(1.3)
w(x, 0) = w0 (x)



w2
where f (w) = 2
. Characteristic curve, x(t), is chosen so that

d dx
dt
w(x(t), t) = ∂t w + ∂ w
dt x

= 0.

Note that w is constant along the curve satisfying x0 (t) = f 0 (w) = w, and the curve

is a straight line. For smooth initial condition, we can construct a solution implicitly

with characteristics:

w(x, t) = w0 (x − f 0 (w)t) = w0 (x − wt).

This seems like we completely solved a Burgers’ equation, but indeed it is only

valid up until the solution is smooth. When faster information along its character-

istics overtakes slower one, the characteristics intersect. So what will happen is the

smoothness breaks, and the PDE becomes invalid because of a multi-valued solution.

Allowing discontinuity, so-called shock wave, is a way to extend the solution beyond

the time when the solution becomes multi-valued. The jump condition across the

discontinuity relies on an integral form of the conservation laws.

The integral form of (1.1) on domain [x1 , x2 ] and time interval [t1 , t2 ] reads
Z x2 Z x2 Z t2  
(1.4) W (x, t2 ) dx = W (x, t1 ) dx + F (W (x1 , t)) − F (W (x2 , t)) dt.
x1 x1 t1

This provides an insight of shock jump condition. As Figure 1.1 indicates, in an

infinitesimal area including a shock path, (1.4) can be approximated by

WL ∆x ≈ WR ∆x + ∆t(FL − FR )
 4

t
Shock

WL
∆t
WR

∆x

Figure 1.1: Shock jump condition

which leads to the Rankine-Hugoniot jump condition

∆x
[F ] = FR − FL = (WR − WL ) = s[W ]
∆t

∆x
where [·] = (·)R − (·)L and shock speed s ≈ ∆t
. This shows that the shock jump

condition ensures the conservation laws.

The solution that is piecewise smooth separated by discontinuity curves and sat-

isfies PDE in smooth parts and shock jump condition across discontinuities is called

a weak solution. Weak solutions satisfy

Z +∞ Z ∞ Z +∞
(1.5) W φt dx + F (W )φx dt = W (x, 0)φ(x, 0) dx
−∞ 0 −∞

for all test functions φ that are smooth and of compact support. See [51, 35] for

details. Note that smoothness in W is not necessary anymore by transferring its

requirement to test functions φ. It turns out that weak solutions are not unique. The

physically relevant solution can be obtained by the principle of stability or checking

the Entropy condition. One can either regularize the initial data by a small amount

or the equation itself by adding a small viscous term, and then apply limiting process

to find the vanishing viscosity solution. This is based on the principle that solution

should be stable to small perturbations: if problem is changed by a small amount,

then the solution should change by a small amount. The Entropy condition is the

condition that picks a solution of physical relevance among many weak solutions. A
 5

discontinuity propagating with speed s satisfies the Lax Entropy condition if

f 0 (wL ) > s > f 0 (wR )

for convex f [32]. Geometrically it says that shock wave is admissible if characteristics

converge into the shock. A more general entropy condition due to Oleinik, valid for

non-convex f , says about the entropy solution with speed s: all discontinuities have

the property that

f (w) − f (wL ) f (w) − f (wR )
≥s≥
w − wL w − wR
for all w between wL and wR [51]. Another approach to formulate the Entropy

condition is to use entropy functions and entropy fluxes. See [35] for details.

The Riemann Problem associated with a hyperbolic conservation laws is to find

a solution, called the Riemann solution, with an initial condition:


 WL , x < 0

(1.6) W (x, 0) =
 WR , x > 0.


Understanding the Riemann solution provides both insights into the physical system

itself and critical numerical perspectives due to the emergence of shock phenomenon.

Example I.1. Consider Burgers’ equation (1.3) with


 wL = 1, x < 0

w(x, 0) =
 wR = 0, x > 0.


Characteristics are straight lines starting from the initial profile, see Figure 1.2a.

They are converging as fast information with speed f 0 (1) catches up with the slow

information with speed f 0 (0). Therefore, the Riemann solution is a shock wave

connecting two states wL and wR , traveling with speed

1 2
[f ] w − 12 wL2
2 R 1
s= = = ,
[w] wR − wL 2

see Figure 1.2b. 

 6

f’(wL) f’(wL)
t ? f’(wR)
t
f’(wR)

x x
x=0 x=0
(a) Characteristics creating a shock wave

1.2

0.8

0.6

0.4

0.2

0
t=0
t=150
-0.2
-100 -80 -60 -40 -20 0 20 40 60 80 100

(b) Riemann solution

Figure 1.2: Example I.1

Example I.2. Now consider (1.3) with


 wL = 0,

x<0
w(x, 0) =
 wR = 1,

x > 0.

Characteristics are nondecreasing as x increases, see Figure 1.3a. One may con-
f (wL )−f (wR )
sider a weak solution that contains a discontinuity with speed wL −wR
. However,

this solution violates the Entropy condition. Therefore, the Riemann solution is a

rarefaction wave as in Figure 1.3b 

 7

f’(wL)

t f’(wR)

x
x=0
f’(wL)
? (x)
t f’(wR)
(o)
x
x=0
f’(wL)

t f’(wR)

x
x=0
(a) Characteristics creating a rarefaction wave

1.2
t=0
t=80
1

0.8

0.6

0.4

0.2

-0.2
-100 -80 -60 -40 -20 0 20 40 60 80 100

(b) Riemann solution

Figure 1.3: Example I.2

1.1.2 Numerical Methods

Foundational numerical methods in particular finite volume method (FVM) is

focused here and will be used throughout this dissertation for PDE-based computa-

tions. Let us consider

(1.7) ∂t W + ∂x F (W ) = 0.

FVM deals with a cell average by defining it at time tn over [xj− 1 , xj+ 1 ] as
2 2

Z xj+ 1
1 2
Wjn ≈ W (x, tn ) dx.
∆x xj− 1
2
 8

Wjn+1
n+1
t

n
Wj-1n Wjn Wj+1n
t
. . .
xj-1 xj-1/2 xj xj+1/2 xj+1
x

Figure 1.4: Schematic for numerical computation

Here ∆x is a mesh width and ∆t is a time step. FVM approximates the integral

form (1.4) of (1.7) on the shaded region in Figure 1.4. One numerical feature in

conservation laws is to maintain its conservation. A conservative method requires

∆t
(1.8) Wjn+1 = Wjn − [F 1 − Fj− 1 ]
∆x j+ 2 2

where Fj+ 1 is a numerical flux defined at the cell interface as some function of,
2

for example, Wjn and Wj+1

n
. Different choices of numerical flux provide different

methods. There is a necessary stability condition, called CFL condition, that the

domain of dependence of numerical method should contain the domain of dependence

of the PDE. This condition was originally derived for finite difference methods for

PDEs by Courant, Friedrichs, and Lewy [13]. Nonlinear hyperbolic systems yields:

∆t
(1.9) λp (Wjn ) ≤1
∆x

for all p at each Wjn , where λp is the p th eigenvalue of the Jacobian matrix F 0 (W ).

Lax-Wedroff theorem says that if numerical solutions by a conservative method con-

verge to some function as the grid is refined, then this function is a weak solution of

the conservation law [31]. Godunov’s scheme is one conservative method proposed in

1959 [18]. It numerically solves exact Riemann solution on every cell interface, which

is expensive. We introduce the Upwind scheme by linearizing the Jacobian matrix

 9

at each cell interface developed by Roe [48, 49, 50]. This is one kind of approximate

Riemann solvers. In a quasilinear form, (1.7) becomes

(1.10) ∂t W + A(W )∂x W = 0.

Roe’s method approximates A(W ) at xj+ 1 by Ā(Wj , Wj+1 ) satisfying

2

1. Ā(Wj , Wj+1 )(Wj − Wj+1 ) = F (Wj ) − F (Wj+1 ),

2. Ā(Wj , Wj+1 ) is diagonalizable with real eigenvalues, and

3. Ā(Wj , Wj+1 ) → F 0 (W̄ ) smoothly as Wj , Wj+1 → W̄ .

Unless we deal with a scalar model, the choice of Ā is not unique. The scheme can

be expressed by

∆t n + o
(1.11) Wjn+1 = Wjn − Āj− 1 (Wjn − Wj−1
n
) + Ā−
j+ 21
(W n
j+1 − Wj
n
)
∆x 2

and

X X
(1.12) Ā+ ∆W = αk λ̄k r̄k , Ā− ∆W = αk λ̄k r̄k
λ̄k >0 λ̄k <0

where ∆(·) = (·)i+1 − (·)i . Here r̄k /λ̄k are the eigenvectors/eigenvalues of local

linearization Ā at xj+ 1 . Because of the second requirement, wave strengths αk can

2

be uniquely determined:
X
∆W = αk r̄k .
k

The first requirement implies that this method guarantees conservation.

When there is a source term, the hyperbolic balance laws in a quasilinear form

read

(1.13) ∂t W + A(W )∂x W = S(W ).

 10

Source term, S(W ), can be treated by means of Upwind as well. The scheme includes

the source term by being projected onto the eigenstructure of Ā [50]:

∆t n + − −

o
Wjn+1 = Wjn − +
Āj− 1 ∆W n − ∆xS̄j− n


(1.14) 1 + Āj+ 12
∆W − ∆x S̄ j+ 12
∆x 2 2

where

(1.15)
X X
Ā+ ∆W − ∆xS̄ + = (αk λ̄k − ∆xβk )r̄k , Ā− ∆W − ∆xS̄ − = (αk λ̄k − ∆xβk )r̄k .
λ̄k >0 λ̄k <0

Wave strengths αk and βk can be determined:

X X
∆W = αk r̄k , S̄ = βk r̄k
k k

where S̄ is some averaging source term at the cell interface.

Closing Remarks: More details on theory and different numerical methods includ-

ing higher order schemes can be found in [32, 51, 53, 35].

1.2 Traffic Field Data

This section describes traffic field data used in this dissertation. Dataset was

collected under Next Generation SIMulation (NGSIM) project1 . Data from the in-

nermost lane of US Highway 101 (US 101) and lane 22 of Interstate 80 (I-80) were

used because it is less likely to have trucks or large vehicles that may have differ-

ent driving maneuver than regular cars, and less vehicles changing lanes. US 101

data records vehicle trajectories and instantaneous velocity on the southbound of US

101 in Los Angeles, CA, between 7:50 a.m. and 8:35 a.m. on June 15, 2005. The

study area is approximately 2, 100 feet in length. Vehicle trajectories of I-80 data are

recorded on northbound direction of Interstate 80 in San Francisco, CA, between 5:00

1 NGSIM dataset is open to public for the purpose of research: https://data.transportation.gov.
2 Lane 1 of I-80 is designated as high-occupancy vehicle lane.
 11

(a) US 101 (b) I-80

Figure 1.5: Study area schematic and camera coverage

p.m. and 5:30 p.m. on April 13, 2005. The site is approximately 1, 650 feet in length.

More specific study area location is in Figure 1.5. Both US 101 and I-80 data provide

the precise location of every vehicle within the study area every one-tenth of a second

along with detailed lane positions. See Figure 1.6a for the first 3 minutes of vehicle

trajectories. Our objective is to extract a manageable macroscopic information

from this microscopic-collected dataset without losing much information. Therefore

2000
45
1800
40
1600
35
1400
30
1200
25
1000
20
800
15
600

400 10

200 5

500 1000 1500 2000 2500

(a) Vehicle trajectories color-coded with instanta- (b) Mean velocity map (mph) from 7:50 am - 8:35 am
neous velocity (ft/sec) from 7:50 am - 7:53 am

Figure 1.6: Two velocity visual representations of US 101 (lane 1)

 12

1600

1400 25

1200
20

1000

15
800

600 10

400
5

200

200 400 600 800 1000 1200 1400 1600 1800

(a) Vehicle trajectories color-coded with instanta- (b) Mean velocity map (mph) from 5:00 pm - 5:30 pm
neous velocity (ft/sec) from 5:00 pm - 5:02 pm

Figure 1.7: Two velocity visual representations of I-80 (lane 2)

the overall emphasis in this section is on the method how we transform given micro-

scopic format into macroscopic format. The variables describing macroscopic flow

are: density ρ(x, t)= number of vehicles per unit length, average velocity v(x, t), and

flux or the rate of flow f (x, t)= number of vehicles per unit time. Moreover, they are

related in a fundamental way that f (x, t) = ρ(x, t)v(x, t). This implies that two vari-

ables are enough to describe traffic flow. Hereafter, we drop (x, t) for easy reading.

Later we dive into the continuity equation for vehicle conservation that relates these

variables. To compute macroscopic variables ρ and v given the trajectory dataset,

the t − x plane is discretized by choosing rectangles of ∆x = 50 (feet) and ∆t = 30

(sec)3 unless we specify otherwise, and the discrete mesh points (xi , tj ) by xi = i∆x

and tj = j∆t are defined. Within each rectangle [xi , xi + ∆x) × [tj , tj + ∆t), we

average the measurement of the speeds of individual vehicles using the space mean

speed. Space mean speed is defined as the harmonic mean of speeds traversing a
∆x Pn n
point xi + 2
during time ∆t: v = 1/vi
. f is computed by counting the number
i=1

3 ∆x and ∆t need to be carefully chosen. If this is too small, information is not aggregated enough and if this is

too large, important information may be overlooked.

 13

∆x
of vehicles, denoted by N , that pass a point xi + 2
during time ∆t:
Z tj +∆t
(1.16) f dt = N
tj

at x = xi + ∆x
2
for all i, j. Then ρ by the fundamental relation ρ = fv . These values

are assigned to each rectangle, see Figure 1.6b for an example of average velocity v

over entire space and time for US101. We found that the on-ramp is controlled by

a ramp meter4 , and this probably accounts for the quasi-periodic pattern in Fig 1.6.

The mathematical models of traffic flow that are studied in this dissertation generally

neglect lane changing, so we only studied the data from lane 1, which has only one

neighboring lane to change into. In fact , it can be deduced from the data in [8] that

fewer than half of the vehicles observed did change lanes during the observation, so

the data are quite suitable for our analysis. Corresponding data for I-80 is in Figure

1.7. All testing that follows in section 3.1.2 is done on these aggregated ρ and v.

Closing Remarks:

1. Theory, computation, and applications of traffic flow are increasingly recog-

nized as autonomous and connected vehicles have been a rising issue. Full under-

standing of mixed flow of human and autonomous vehicles requires good knowledge

of human-driving models, which are main theme in this work, and from which trans-

portation engineers can benefit.

2. The rest of the dissertation is organized as follows. Chapter II reviews micro-

scopic models and presents its interesting simulation. Chapter III concerns macro-

scopic models comprehensively: it opens up from existing macroscopic models, and

closes by proposing new multilane models. Chapter IV lists the proposed model-

based numerical results along with numerical methods. Finally Chapter V concludes
4 This can be deduced from Google Earth.
 14

the dissertation by summarizing the contribution of this work. Future research topics

and potential applications are also highlighted in Chapter V.

 CHAPTER II

Microscopic Car-Following Model

As mentioned in Chapter I, mathematically traffic flow can be described by sys-

tems of ordinary differential equations (ODEs), known as microscopic car-following

models, or PDEs that are macroscopic fluid-like models. The particle-based model

illustrates the motion of individual vehicles by taking its acceleration into account.

Although the main interest in this dissertation pertains to macroscopic PDE mod-

els, we start with understanding how individual vehicles act and communicate each

other. This will be used to relate particle-level description to fluid-level description.

For this purpose, we review some microscopic models and present some numerical

results.

2.1 Literature Review

Car-following model, sometimes referred to as the follow-the-leader model, has

been extensively developed since 1950s. Its underlying idea is that (n)th vehicle’s

response is based on its sensitivity and stimulus from (n+1)th vehicle that is assumed

to be placed ahead, see Figure 2.1.

(2.1) Responsen (t) = Sensitivityn (t − τn ) · Stimulusn (t − τn )

15
 16

hn

xn-1 xn xn+1 xn+2

Figure 2.1: Car-following schemetics

where τn is reaction time of (n)th driver. Given assumptions that vehicles are iden-

tical and drivers react in the same way, same sensitivity for instance, movement of

vehicles may be governed by fundamental relation of (2.1). In 1958, Chandler et al.

studied

(2.2) ẍn (t) = α[ẋn+1 (t − τn ) − ẋn (t − τn )]

where xn (t) is the position of (n)th car at time t, and α is a constant in unit of

time−1 [10]. ˙ indicates a derivative taken with respect to t. Responsen in (2.1)

is replaced with (n)th drivers’ acceleration, Sensitivityn with a constant parameter,

and Stimulusn with relative speed between (n + 1)th and (n)th vehicles. This simply

illustrates that (n)th vehicle accelerates (or decelerates) when it travels slower (or

faster) than (n + 1)th vehicle. Parameter α in the model is measured by field traffic

data. This model is simple but has a drawback that it has a constant sensitivity

that is too limited. Since this model, there have been many improvements. In 1959,

Gazis et al. proposed

α
(2.3) ẍn (t) = (ẋn+1 (t − τn ) − ẋn (t − τn ))
hn (t − τn )

where hn (t)=xn+1 (t) − xn (t) is the headway at time t, see again Figure 2.1, and α

is a constant [16]. This incorporates headway into sensitivity: (n)the vehicle is less

affected by the car ahead when headway is large enough, and vice versa. General

Motors Nonlinear Model was followed by Gazis et al. in 1961:

ẋn (t)β
(2.4) ẍn (t) = α (ẋn+1 (t − τn ) − ẋn (t − τn ))
hn (t − τn )γ
 17

where α, β, γ are constants to be determined [17]. Note that (2.2) and (2.3) are

special forms of (2.4). More realistic ideas were followed. In 1968, Bexelius developed
N
X
(2.5) ẍn (t) = αi (ẋn+i (t − τn ) − ẋn (t − τn )).
i=1

Here (n)th car’s response is not just due to (n + 1)th car, but also to (n + i)th cars

for i = 2, 3, ..., N weighed by αi [7]. In 1993, Ozaki suggested to treat acceleration

and deceleration differently, and empirically estimated parameters using (2.4) [45].

Recent work proposed by Zhang in 2002 is

ẋn+1 (t) − ẋn (t)

(2.6) ẍn (t) =
τ (hn (t))

where τ is the response time [59]. Rather than having a constant factor in sensitiv-

ity and reaction time, Zhang uses an idea of response time, denoted by τ , varying

with headway. (2.6) was used to derive macroscopic acceleration equation or the

continuum version:

h
(2.7) ∂t v + v∂x v ≈ ∂x v
τ (h)

where v represents the mean velocity and h is the macroscopic headway. Zhang

formulated
h
∂x v = −ρV 0 (ρ)∂x v
τ (h)
where V is some equilibrium velocity that depends on the local density. ρV 0 (ρ) can

be interpreted as the traffic sound speed at which small traffic disturbances propagate

relative to a moving traffic stream.

In 1995, Bando et al. suggested a model that is independent of ẋn+1 :

(2.8) ẍn (t) = α[VB (hn (t)) − ẋn (t)]

where α is a constant. Here they designated VB as some desired velocity function.

This model says that what a driver really follows is some optimal velocity he/she has
 18

in mind that depends on the front headway. For this, this model is called Optimal

Velocity (OV) model. In spite of its simple form, this was a big leap towards a

realistic traffic flow because it supported the formation of stop-and-go flow pattern

[2, 4]. Note that OV model does not take the reaction time into consideration, so

the consequence of this model implies instantaneous reaction. In 2000, Berg et al.

also made an effort to translate a set of OV ODEs to PDE using a series expansion

of the headway in terms of ρ [6]. Berg formulated

h2 ∂x ρ h3 ∂xx ρ
(2.9) hρ + + + ··· = 1
2! 3!

by expending
Z h(x,t)
(2.10) ρ(x + y, t) dy = 1
0

in powers of y, and integrating. Keeping up to cubic terms in (2.9), expansion in a

series in terms of ρ takes the following approximation:

1 ∂x ρ ∂xx ρ (∂x ρ)2

(2.11) h≈ − − + + ··· .
ρ 2ρ3 6ρ4 2ρ5

And then a continuum version was derived

h ∂ ρ ∂ ρ (∂ ρ)2 i
x xx x
(2.12) ∂t v + v∂x v = α[V̄ (ρ) − v] + αV̄ 0 (ρ) + −
2ρ 6ρ2 2ρ3

where V̄ (ρ) is related to VB in a way that V̄ (ρ) = VB (1/ρ).

In 2001, Jiang et al. combined the OV model with the classical model by :

(2.13) ẍn (t) = κ[V (hn ) − ẋn (t)] + λ[ẋn+1 (t) − ẋn (t)]

where V (hn ) is some equilibrium velocity pertaining to headway, and κ and λ are

sensitivity constants. This combined model supported not only the formation of

stop-and-go congestion but also hysteresis loop1 [26]. We refer to [3, 43, 44] and the

references therein for the study with addition of delay time.

1 This is a loop typically shown in speed-density plane that accounts for acceleration and deceleration branches

separately.
 19

2.2 Numerical Result

We carry out numerical experiments of (2.13) first to obtain an insight of the

vehicle motions that cause inevitable instability by perturbing initial condition. With

the assumption of N vehicles on a ring road, initial condition is the uniform flow with

a small disturbance by one vehicle coming in (Test 1) as in Figure 2.2. To be specific,

Test 1 has an initial disturbance:

xn (0) = n NL for n < 50

xn (0) = (n − 12 ) NL for n = 50

xn (0) = (n − 1) NL for n > 50

where L is the total length of the ring road, and n = 1, · · · , N + 1. This initial small

perturbation can be viewed as random disturbance or consequence of lane change of

n = 50th car. In this simulation, L = 1500 (m), N = 100,

V (h) = 6.75 + 7.91 tanh[0.13(h − 5) − 1.57] (m/s),

κ = 0.41 s−1 , and λ = 0.5 s−1 were taken. With this choice of parameters, numerical

experiments show that V is nonnegative. It, however, could go negative by varying

the value of λ. For details, see hysteresis loops in [26].

Test 1 exhibits that the cars start off with decelerating because of the one car

merged, and all cars experience break-acceleration-break travel at the end, see Fig-

ure 2.3. Interestingly, this also happens when one car is removed, see Test 2 in

Figure 2.3. A real experiment was similarly set up and two supplementary videos

can be found in [52].

In [2, 26], bifurcation study was done under linear stability analysis:

κ
(2.14) V 0 (h) < + λ.
2
 20

ow ow ow
n of fl n of fl n of fl
tio tio tio

ec
ec

ec
#1 #N #1 #N+ #1 #N-

Dir
Dir

Dir
1 1

#N

#2

#N
#2

#2
#N . . .
-1 . . .

-2 . . .
.
..
..

..
. .

leaving
coming in coming in
(a) Uniform flow (b) Test 1 (c) Test 2

Figure 2.2: Uniform flow (a), disturbed flow with one entering car (b) and leaving car (c).

1000 1000

900 900

800 800

700 700

600 600

500 500

400 400

300 300

200 200

100 100

0 0
0 500 1000 1500 0 500 1000 1500

14 14

12 12

10 10

8 8

6 6

4 4

2 2

0 0
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000

(a) Test 1: one car moving into lane, total N + 1 vehicles (b) Test 2: one car leaving out, total N − 1 vehicles

Figure 2.3: Top: car trajectories with #1 car trajectory (bold line). Bottom: speed of #1 car.

If (2.14) is satisfied, then small disturbance always dies out. If not, then small

disturbance may result in stop-and-go formation.

From macroscopic perspective, this test indeed indicates a richness of adding the

acceleration equation when modeling is concerned. Later we will take (2.13) to

convert to PDE and then reproduce this stop-and-go pattern.

 21

Closing Remarks: Microscopic approaches are, however, computationally expen-

sive, as each vehicle has an ODE to be solved at each time step. Therefore as the

number of vehicles increases, so does the size of system to be solved. For this rea-

son, it is desirable to use the macroscopic models if a good model can be found

satisfactorily.
 CHAPTER III

Macroscopic Traffic Flow Model

Chapter III is concerned with macroscopic models. We begin with the scalar

conservation laws of vehicles. In order for this core equation to be solvable, it requires

a closure that relates mean velocity and density. A new fundamental relation is

proposed. A discussion of discontinuous flux and non-convex flux is included. Then

we extend to newly-developed 2x2 system model. Built on all of these, Chapter III

closes with presenting two types of multilane models: scalar model (Type 1) and

system model (Type 2).

3.1 Scalar Conservation Laws

In macroscopic models the interest is the collective average of traffic flow and

its evolution. The underlying fundamental principle is the car mass conservation,

resulting in a conservation law for the car density. In this section, we derive the

mass conservation equation and discuss simple traffic examples.

3.1.1 Derivation of Equation

Consider a finite road segment [x1 , x2 ] during a time period [t1 , t2 ]. Then
Z x2
total number of vehicles in [x1 , x2 ] at time t = ρ(x, t) dx
x1

22
 23

where 0 ≤ ρ ≤ ρmax . ρmax is the state at which vehicles are touching bumper to
ρ
bumper. ρ = 0 implies an empty road. By dividing ρ by ρmax , 0 ≤ ρ̃ = ρmax
≤ 1. We

drop ˜ throughout the dissertation. Assuming that there is no on/off-ramp, the total

number of vehicles in this section can change only because of traffic flowing across

the boundaries x1 or x2 :
Z x2 Z x2 Z t2  
(3.1) ρ(x, t2 ) dx = ρ(x, t1 ) dx + ρ(x1 , t)v(x1 , t) − ρ(x2 , t)v(x2 , t) dt.
x1 x1 t1

This is the integral form of conservation laws. For differentiable ρ and v, the funda-

mental theorem of calculus leads

Z t2 Z x2
(3.2) ∂t ρ + ∂x (ρv) dx dt = 0.
t1 x1

Because x1 , x2 , t1 and t2 were chosen arbitrarily, we reach a conclusion that the

integrand in (3.2) vanishes, yielding the well-known LWR model (Lighthill-Whitham

1955 and Richards 1956) [37, 47]:

(3.3) ∂t ρ + ∂x (ρv) = 0.

This classical LWR model, a nonlinear first order PDE, treats traffic flow as a one-

directional compressible fluid and studies properties induced by the interactions of

vehicles as a whole.

We assume that vehicles drive at around speed limit of the road in light traffic,

and slow down in increasing traffic. Therefore the average velocity of vehicles can

be an a-priori relationship that describes v as a function of local ρ, v = V (ρ).

Depending on a variety of ways of formulating the equilibrium state velocity V (ρ),

this simple but insightful scalar model allows various traffic waves to arise: shocks,

rarefaction waves, or compound waves. This relationship between v and ρ is, in

general, called ’the speed-density relationship’ in traffic flow or Fundamental diagram

 24

60 5000

4500
50
4000

3500
40
3000

30 2500

2000
20
1500

1000
10
500

0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 3.1: The Greenshields’ model

(FD)1 equivalently. This relation follows the general rules that v is non-increasing

i.e. V 0 (ρ) ≤ 0. Moreover, v tends to flatten out to zero speed when ρ approaches 1

while v stays near speed limit when ρ is very low and close to 0. The LWR model

with a proper speed-density relation has been an indispensable building block since

its appearance. The LWR model (3.3) with the closure relation becomes

(3.4) ∂t ρ + ∂x ρV (ρ) = 0.

A simple relationship was proposed by Greenshields et al. [20]:

(3.5) V (ρ) = vmax (1 − ρ).

Together with (3.5), (3.4) is a nonlinear scalar equation. See Figure 3.1 for the

speed-density relationship as well as FD of (3.5).

Closing Remarks: We close this section by revisiting simple Riemann problem

associated with shock and rarefaction waves in traffic flow.

1 In fact, FD refers to a diagram describing flux and ρ relation. Because of f = ρv, speed-density relationship and

FD are interchangeably used.

 25

f’(ρR)
t f’(ρL) t

x x
x=0 x=0
(a) Characteristics creating a shock wave (b) Visualization of vehicle trajectories

Figure 3.2: Example III.1

Example III.1. Consider (3.4) with (3.5) with Riemann data: ρL = 0.5 and ρR = 1.

This models the situation in which vehicles moving at speed VL > 0 unexpectedly

encounter a bumper-to-bumper traffic jam and slam on their brakes, instantaneously

coming to a stop while the density jumps from ρL to ρR . We notice that character-

istics in each of the regions where ρ is constant go into the shock as time evolves as

in Figure 3.2a, and vehicle trajectories are shown in Figure 3.2b. Therefore, shock

wave forms and propagates backwards with speed vmax (1 − ρL − ρR ). 

Example III.2. Consider instead with ρL = 0.8 and ρR = 0.3.

The vehicles to the left move initially slow but can begin to accelerate once the

vehicles in front of them begin to speed up. Under the assumption of V (ρ) here, each

driver can speed up only by allowing the distance between her/him and the leading

car to increase, and so we see a gradual acceleration and spreading out of vehicles.

Vehicles are rarefied as this wave passes through, see Figure 3.3. 

Lorem ipsum

f’(ρR)
f’(ρL)
t
t

x x
x=0 x=0
(a) Characteristics creating a rarefaction wave (b) Visualization of vehicle trajectories

Figure 3.3: Example III.2

 26

3.1.2 Fundamental Relation and a New Fundamental Diagram

Since the earliest work measured and quantified empirically in 1930s by Green-

shields et al. [20], there has been many studies on the relationship between flow

speed and flow density due to easier access to field data, summarized in Table 3.1. A

fit to the Lincoln tunnel data was found by Greenberg in 1959 [19], Underwood and

Drake et al. proposed exponential forms separately [54, 15], and Kerner et al. also

deployed an exponential form in their simulation [27]. Del Castillo et al. proposed a

double exponential form [9]. Traffic flow is typically divided into two states, free flow

and congested flow. Based on this, some relationships consist of two phases [57, 38].

See Figure 3.4 for a comparison. Here vf refers to a free flow speed.

Relationship

Greenshields et al. (1935) vf (1 − ρ)

Greenberg (1959) vf ln(1/ρ)
Underwood (1961) vf exp(−ρ)
Drake et al. (1966) vf exp(−(ρ/3)2 /2)
Del Castillo et al. (1995) vf (1 − exp(1 − exp(7(1/ρ − 1)/vf )))
Yang et al. (2011) min{vf , αρm }

Table 3.1: A list of speed-density relationships.

60
Greenshields 1935
Drake 1967
50 Del Castillo 1995
Yang 2011

40

30

20

10

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 3.4: Some visual representations of Table 3.1.

 27

We focus on the two-phase model by Yang et al. . The model is

(3.6) V = min{vf , αρm }

where m < 0 and α > 0. This model was fitted to 2009 traffic data from Interstate

95 in Virginia and α, m were optimized in [57, 11]. In congested flow, m turned out

to be less than −1, which makes f turn from linear to convex as it switches from free

to congested flow. One main drawback is that f in this model immediately decreases

as soon as v drops down from vf . But an observation from traffic data shows that

f actually increases, albeit rather slowly even when v decreases from vf , and then

f starts to decrease as v further decreases. To capture this characteristic from the

data observation, we propose a three-phase model

(3.7) V = min{vf , α1 exp (α2 ρ)ρm1 ρ+m2 }

where α2 ≤ 0 and m1 ≤ 0. Note that this has an analogous structure of (3.6), but is

different in that parameters are replaced by some function of ρ, α(ρ) = α1 exp (α2 ρ),

and m by m(ρ) = m1 ρ + m2 . Choice of such functions are from the data observation

and this is presented in Appendix A. See Figure 3.9 for a schematic of three-phase

FD. It is not easy to catch an improvement in v − ρ plane, but it is noticed that in

f − ρ plane it fixes the drawback with this simple modification.

We take logarithm of (3.7) and reformulate the optimization problem as

min k ln v − ln α1 − α2 ρ − m1 ρ ln ρ − m2 ln ρk22
m1 ,m2 ,α1 ,α2

with α2 ≤ 0, m1 ≤ 0. In Figure 3.5, the black dots associated with v are basically

identical to dataset in Figure 1.6b but displayed for the FD purpose. A pair of (v, ρ)

at a space and time is represented as a dot in v − ρ plane, and it is aggregated over

all domain in Figure 3.5. Linear regression was performed on smoothened data that
 28

50 3000

40
original 2500

smoothed 2000
30
fit 1500
20
1000

10
500

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Figure 3.5: US 101 data represented in v − ρ plane, original in black and smoothed data in green.
Fitted curve is in blue (Left). Same data in f − ρ plane (Right).

obtained by moving window with size of 4 data points and returning its mean. Recall

that the US101 data was collected during peak time so there is a lack in free flow

data. For this reason, we fit the congested regime only. The fit result to congested

part is

(3.8) v = 22.5 exp (−3.3ρ) ρ0·ρ−0.4 .

The fit (3.8) gives expected velocities and therefore the model predicts expected

traffic patterns. Figure 3.5b shows the corresponding f − ρ diagram. It switches

from concave to decreasing convex as ρ increases. Convexity change is one of key

features differing from the two-phase FD. This would have one additional increasing

linear when the free flow part is completed, and this will complete a three phase FD.

Having a non-convex flux in a scalar conservation laws is a completely new regime

to look at, which will be addressed later. With the new three-phase results, we

are capable of reconstructing the predicted mean velocity and analyzing the error2 .

Figure 3.6 reports them. In reconstruction, there are some missing data due to the

lack of data availability for free flow. As two measures of performance, histogram of

errors and its qq-plot describe left light long tail and right heavy tail.

2 By error, it means predicted velocity - observed velocity.

 29

2000 2000
45 45
1800 1800
40 40
1600 1600
35 35
1400 1400
30 30
1200 1200
25 25
1000 1000
20 20
800 800
15 15
600 600

400 10 400 10

200 5 200 5

500 1000 1500 2000 2500 500 1000 1500 2000 2500

(a) Mean velocity map (b) Predicted mean velocity by (3.8)

QQ Plot of Sample Data versus Standard Normal

15

10

Quantiles of Input Sample

0

-5

-10

-15

-20

-25

-30

-35
-4 -3 -2 -1 0 1 2 3 4
Standard Normal Quantiles

(c) Histogram of error (d) Corresponding qq-plot

Figure 3.6: Reconstruction and errors of US 101

The same methodology was also performed on I-80 data. Result is below:

(3.9) v = 20.1 exp (−3.5ρ) ρ0·ρ−0.34 ,

and see Figure 3.7 and 3.8.

Here we present a practical shape of the new FD, so that this can be easily

tractable and applicable in PDE-based problems. For this, we constraint for FD to

30 2500

25 original 2000

20 smoothed
1500
15 fit
1000
10

500
5

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Figure 3.7: I-80 data represented in v − ρ plane, original in black and smoothed data in green.
Fitted curve is in blue (Left). Same data in f − ρ plane (Right).
 30

1600 1600

1400 25 1400 25

1200 1200
20 20
1000 1000

15 15
800 800

600 10 600 10

400 400
5 5
200 200

200 400 600 800 1000 1200 1400 1600 1800 200 400 600 800 1000 1200 1400 1600 1800

(a) Mean velocity map (b) Predicted mean velocity by (3.9)

QQ Plot of Sample Data versus Standard Normal

15

10

Quantiles of Input Sample

0

-5

-10

-15

-20

-25
-4 -3 -2 -1 0 1 2 3 4
Standard Normal Quantiles

(c) Histogram of error (d) Corresponding qq-plot

Figure 3.8: Reconstruction and errors of I-80

be continuous and differentiable. Assume that vf and ρc 3 are provided. A posted

speed limit is a good candidate, or average of max and min speed limit is another

good candidate for vf , and ρc usually assumed to be 10-14% of maximal density. In

US101 and I-80 data, it turned out that m1 ≈ 0. But if we had enough data in free

flow part and the transition part between free flow and congested regime, then m1

may have played a role. Therefore in practical form, we keep m1 . Once m1 and m2

are determined, assigning

m1 ρc + m2
α1 = vf exp(−α2 ρc )ρ−m
c
1 ρc −m2
and α2 = −m1 log(ρc ) −
ρc

will create a smooth flux function as in Figure 3.94 . This is an exact analogous to

put constraints that

V (ρc ) = vf and V 0 (ρc ) = 0.

3A critical density, denoted by ρc , is a break density point that divides a free flow and congested flow in FD.
4 Itis shifting down so that f (ρmax ) = 0 and this makes sense if we consider a static homogeneous flow in
equilibrium state.
 31

concave
50 10
convex
45 9

40 8

35 7

30 6

25 5

20 4

15 3

10 2

5 1

0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 3.9: A practical shape of the new FD (3.7)

In one practical FD, vf = 50, ρc = 0.1, m1 = −0.4, and m2 = 0.4. The maximal flux
√
−α2 − α22 −4m1 (m2 +1)
is achieved at ρ ≈ 2m1
≈ 0.30. Observe that there are three states

of traffic condition, free flow and two types of congested flow that are divided by

the change in convexity nature in flux. Flux concavity changes to flux convexity at

ρ ≈ 0.525 . The practical form may be changed depending on different choices of

parameters such as vf , but the general shape with three phases- free, two congested

in concave to convex- should remain as in Figure 3.9.

Closing Remarks: One can illustrate traffic flow under certain circumstances by

constructing an appropriate form of V (ρ). For night time driving example, see [34].

5 This requires a simple algebraic calculation that needs to solve a 4th degree polynomial.
 32

3.1.3 Riemann Problem

In this section, we discuss Riemann problems with spatially varying flux and

non-convex flux, which lead to more complicated wave structure [35].

Spatially Varying Flux

It is not unusual that we encounter variable speed limit area either reduced or

increased speed limit zone when driving. For instance this may be for the purpose

of real-time control on the highway traffic volume, or for the specific purpose such

as school zone or highways with changing surface conditions [40]. Sometimes, even

the road that is unpaved or uneven forces drivers to reduce their speed. This results

in average velocity of traffic that depends not only on the local density but also on

the space.

Imagine that a road has a speed limit 55 mph (x < 0) and it drops down to 35

mph (x > 0), see Figure 3.10. Consider

(3.10) ∂t ρ + ∂x ρV (ρ, x) = 0

where

(3.11) V (ρ, x) = Vmax (x)(1 − ρ)

and

 55,

x<0
(3.12) Vmax (x) =
 35,

x > 0.

Solution of the Riemann problem consists of

1. left-going and right-going waves (either shock or rarefaction wave), and

2. a stationary discontinuity at x = 0.
 33

REDUCE
SPEED
AHEAD

speed limit = 55mph 35mph

x=0

Figure 3.10: Road example of the variable speed limit

S/R SD S/R

ρL,M ρR,M

t ρL ρR

x
x=0

Figure 3.11: Riemann solution structure in terms of wave propagation

.

See Figure 3.11 for the solution structure. In Figure 3.11, S implies shock, R rar-

efaction, SD stationary discontinuity, and ρL,M /ρR,M denote the intermediate states

to the left/right of x = 0, respectively. Let us denote the flux for x < 0 by fL ,

and for x > 0 by fR . Since fL governs the propagation of left-moving waves and fR

governs the propagation of right-moving waves, some parts of the flux curves are not

admissible. The discontinuous flux is illustrated in Figure 3.12 where the locus of

admissible curves are in solid curves.

14

12

10
fL

6
fR
4

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 3.12: Graphs of admissible fL and fR in solid curves

.
 34

Rankine-Hugoniot jump condition

fL (ρL,M ) − fR (ρR,M ) = s(ρL,M − ρR,M )

across a stationary discontinuity with s = 0 implies that f is continuous at x = 0, or

(3.13) fL (ρL,M ) = fR (ρR,M ).

Depending on the choice of ρL and ρR , different wave patterns may arise. We

illustrate two cases below. The Riemann solution is illustrated geometrically in

Figure 3.13.

Example III.3. Solve (3.10) with (3.11), given


 0.4, x < 0

ρ(x, 0) =
 0.3, x > 0.


The Riemann solution contains a left-going shock from a state ρL = 0.4 to a

fL (ρL )−fL (ρL,M )
state ρL,M ≈ 0.86 , with speed ρL −ρL,M
, a stationary jump from ρL,M = 0.8 to

ρR,M = 0.5 at x = 0 along with a rarefaction fan from ρR,M = 0.5 to ρR = 0.3. The

traffic slows down since the speed limit decreases at x = 0, and this creates a shock

wave going backwards, and then gradually speeds up through a rarefaction wave. 

Example III.4. Consider


 0.9,

x<0
ρ(x, 0) =
 0.3,

x > 0.

Then the Riemann solution contains a left-going rarefaction wave, a stationary

discontinuity, and a right-going rarefaction wave. 

6 This is one critical value satisfying fL (ρ) = fR ( 12 ).
 35

14

12
(ρL,fL(ρL)) . 1

0.9
t=0
t=2

0.8

10
. . 0.7

6
.
(ρR,fR(ρR))
0.6

0.5

0.4
fL(ρL,M)=fR(ρR,M)
4 0.3

0.2
2
0.1

0 0
0 0.2 0.4 0.6 0.8 1 -100 -50 0 50 100

(a) Example III.3

14 1
t=0
0.9 t=2
12
0.8
fL(ρL,M)=fR(ρR,M)

. . .
10 0.7

0.6
8
0.5

.
(ρR,fR(ρR))
6
0.4
(ρL,fL(ρL))
4 0.3

0.2
2
0.1

0 0
0 0.2 0.4 0.6 0.8 1 -100 -50 0 50 100

(b) Example III.4

Figure 3.13: Riemann solution locus and evolution of density.

Non-convex Flux

Recall that the data-driven three-phase flux in section 3.1.2 turned out to be

non-convex with two inflection points. The solution to the Riemann problem can be

determined from the graph of f (ρ) by the convex-hull7 method [35]. If ρR < ρL , we

construct the convex hull of the set {(ρ, y) : ρR ≤ ρ ≤ ρL and y ≤ f (ρ)} and look at

the upper boundary of this set to determine the structure of the Riemann solution. If

ρL < ρR , then we follow instead the convex hull of the set above the graph y = f (ρ),

{(ρ, y) : ρL ≤ ρ ≤ ρR and y ≥ f (ρ)}, and read off the lower boundary of this set.

Note that if f (ρ) is convex (or concave), then the convex hull construction gives

either a single shock or a single rarefaction as we expected. Let us take a look at

7 The convex hull is defined as the smallest convex set containing the original set.
 36

10 1
t=0
9 0.9 t=10

8 0.8

7 0.7
Convex hull
6

5
. 0.6

0.5
(ρL,f(ρL))
f’(ρL) f’(ρR)
4
.
(ρR,f(ρR))
0.4

3 0.3
t
2 0.2

1 0.1
x
0 0 x=0
0 0.2 0.4 0.6 0.8 1 -100 -50 0 50 100

Lorem ipsum
(a) Example III.5 (b) Characteristics

Figure 3.14: Convex hull, evolution of density, and characteristics of Example III.5

two examples to seek for Riemann solution by constructing convex hull.

Example III.5. Solve

(3.14) ∂t ρ + ∂x f (ρ) = 0

with

 0.15,

x<0
ρ(x, 0) =
 0.7,

x>0

where f (ρ) is the three-phase FD in Figure 3.9.

By constructing the convex hull, it is composed of a straight line segment from

(ρL , f (ρL )) to (ρR , f (ρR )), see Figure 3.14a. The straight line represents a shock

jumping from ρL to ρR . Therefore, the Riemann solution contains a shock wave

going backwards, see Figure 3.14a and its corresponding characteristics in b. 

Example III.6. Consider


0.7 x<0


ρ(x, 0) =
 0.15,

x > 0.

By constructing the convex hull, it follows from (ρR , f (ρR ) to (ρ∗ , f (ρ∗ )) and then

has a straight line segment down to (ρL , f (ρL )), see Figure 3.15a. The segment where
 37

10 1
t=0
9 0.9 t=2

8 0.8
f’(ρ*)
7 .
(ρR,f(ρR))
. (ρ*,f(ρ*))
0.7

6 0.6

5 0.5

4
Convex hull
. 0.4 f’(ρR)
3 0.3
f’(ρL)
2 0.2
t
1 0.1

0
0
.
ρR
0.2 0.4 0.6
.
ρL
0.8 1
0
-100 -50 0 50 100
x=0
x

Lorem ipsum
(a) Example 6 (b) Characteristics

Figure 3.15: Convex hull, evolution of density, and characteristics of Example III.6

boundary follows is the rarefaction wave and the straight line is a shock jumping

from ρ = ρ∗ to ρ = ρL , where ρ∗ is the inflection point. Hence, the Riemann

solution contains a shock wave propagating to the left followed by a rarefaction, see

Figure 3.15a, as a compound wave. Figure 3.15b is its characteristics. 

Example III.7. Consider


 0.8,

x<0
ρ(x, 0) =
 0.7,

x>0

Both states are in a very congested flow and ρL = 0.8 > ρR = 0.7. In this case,

we end up having a Riemann solution with a shock traveling to the left. This is

the entropy-satisfying solution that makes sense mathematically, however we would

expect a rarefaction fan in reality. This will always happen if flux has the convex

portion. 

Closing Remarks: Variations to scalar model do exist. Appropriate source terms

can be put into model to investigate, for example, on/off-ramps effects [35]. Adding

a source term to capture the stochastic nature of traffic is another way of approach

[57]. Some generalizations of scalar models will briefly follow.

 38

3.1.4 Generalization

One generalization based on LWR model is to consider different types of vehi-

cles known as multi-class or n-population. Heterogeneous vehicles’ characters are

reflected in that automobiles allow a higher speed limit than trucks, and that the

length of automobiles is smaller than that of trucks, etc [56, 58, 5]. Let us take a

look one model of 2-population in [5]

(3.15) ∂t ρi + ∂x (ρi vi ) = 0 for i = 1, 2

 
r
where ρi is the density of vehicles belonging to the ith class. vi = vi (r) = Vi 1− rmax

where r = l1 ρ1 +l2 ρ2 , li is the average length of vehicles in ith class, and Vi is maximal

speed of class i. Assuming li = 1 for i = 1, 2 and rmax = 1, (3.15) becomes

 2
X 
(3.16) ∂t ρi + ∂x ρi Vi (1 − ρj ) = 0 for i = 1, 2.
j=1

This can also viewed as using ’weighted maximal speed’ as follows

(3.17) ∂t ρ + ∂x (ρVm (1 − ρ)) = 0

ρ1 ρ2
where Vm = ρ 1
V + ρ 2
V and ρ=ρ1 + ρ2 . It has an analogous form as (3.4) but

has weighed elements in Vm . In [56], they developed another 2-population model

in which numerical simulations showed two-capacity regimes, hysteresis and platoon

dispersion. The idea in [58] is that the performance of vehicles belonging to different

classes is different in free flow8 but indistinguishable in heavy traffic 9 .

Another generalization is to consider flow of traffic over more than a single lane.

This can be done by introducing a road configuration n(x) to indicate the number

of lanes and consider the lane-aggregate model

(3.18) ∂t n(x)ρ + ∂x n(x)ρV (ρ) = 0
8 This is mainly due to the fact that Vm for automobiles is greater than that for trucks.
9 Both type of vehicles travel at the same group speed
 39

with ρ the average lane car density across n lanes. Varying n(x), one may simulate

lane closure or open-up [33].

Closing Remarks: Despite of variations in LWR model, scalar conservation laws

alone fails to explain complicated realistic phenomena, such as stop-and-go oscilla-

tions [36].

3.2 System

In the LWR model, v = V (ρ) implies that velocity adjusts instantaneously to

density. This is clearly an over simplification. Vehicles react to flow conditions,

they accelerate or decelerate in response to traffic flow around them, possibly in

consideration of vehicle mechanical constraints. The LWR model cannot describe

these processes. In recent years, higher-order models10 have been proposed that

treat velocity v = v(x, t) as a separate dependent variable and describes acceleration,

response time etc.

3.2.1 Literature Review

There are numerous 2x2 systems: some are inspired by gas dynamics and others

are derived from the microscopic acceleration equation [46, 55, 27, 1, 6, 12, 25, 59].

Payne in 1971 and Whitham in 1974 (PW) [46, 55] first proposed a model inspired

by gas dynamics

∂t ρ + ∂x (ρv) = 0
(3.19)
V (ρ)−v P 0 (ρ)∂x ρ
∂t v + v∂x v = T
− ρ

where the acceleration is a summation of relaxation term and pressure term. Here

T is the relaxation time that takes to relax to some equilibrium velocity V (ρ). By
10 From a mathematical point of view, it means ’a system of two-(or three-) equation models’ rather than ’second-(or

high-) order models’. By convention these models are sometimes referred to as high-order models.
 40

convention, V (ρ) can be replaced by one choice among many in section 3.1.2. P (ρ) is

a pressure11 satisfying P (ρ) = ργ , γ > 0. The PW model couples velocity dynamics

as a second equation in order to link ρ and v, viewing them as independent variables.

It is worthwhile to point out that the speed-density relation, V (ρ), plays a part even

in 2x2 system, and this re-emphasize the importance of its study. The properties of

the system are controlled by the eigenvalues of the system. These eigenvalues, also

known as characteristic speeds, determine how traffic information is propagated in

the stream. It turns out that

p p
λ1 = v − P 0 (ρ) < λ2 = v + P 0 (ρ).

We note that the second eigenvalue is greater than the traffic velocity v. It means

that some part of the information travels faster than average flow speed. As a result,

waves associated with the second characteristic always reach vehicles from ’behind’.

This happens in fluid dynamics but is an unrealistic property to traffic. A driver is

an anisotropic particle that responds to frontal stimuli as opposed to the fact that a

fluid particle may respond to stimuli both from the front and behind. It also turned

out that PW model may produce negative travel velocity, i.e. wrong-way travel [14].

In 2000, Aw and Rascle (AR) proposed a new system to address these drawbacks

by simply replacing the space derivative with the convective derivative [1]. Replacing

∂x (P (ρ)) by ρ(∂t +v∂x )P (ρ) in (3.19) leads to AR model ignoring the relaxation term:

∂t ρ + ∂x (ρv) = 0
(3.20)
∂t (v + P (ρ)) + v∂x (v + P (ρ)) = 0

where P (ρ) = ργ , γ > 0, is an increasing function. They proposed five traffic model

principles including so-called Anisotropy condition that prevents rare disturbances

11 In physics, pressure depends on Newton’s second law of motion because a difference in pressure across a surface

implies a difference in force, which can result in an acceleration, if there is no additional force to balance it.
 41

from propagating forward. This is a condition that eigenvalues of the system are at

most the traffic velocity v

λ ≤ v.

With a choice of P (ρ) above, AR model ensures the anisotropic property:

λ1 = v − ρP 0 (ρ) < λ2 = v.

Therefore, every driver reacts to what happens in front and is not affected by dis-

turbance from behind. Qualitative properties that braking/accelerating produces

shock/rarefaction waves is also one of principles, which will be discussed later.

Around the same period, independently, Zhang proposed a model in 2002 [59]:

∂t ρ + ∂x (ρv) = 0
(3.21)
∂t v + (v + ρV 0 (ρ))∂x v = 0.

A micro-to-macro link is established for this model, meaning that the second equation

is derived from the microscopic car-following model (2.6), which naturally carries the

anisotropic property: v + ρV 0 (ρ) = λ1 < λ2 = v.

Closing Remarks: To derive the macroscopic acceleration equation, vehicle-level

models should guide us a reasonably correct way. First of all, traffic flow is not really

the physics as in gas dynamics. It is difficult to interpret pressure (or pressure-

like term) in traffic flow. Second of all, microscopic acceleration model (2.13) can

reproduce realistic flow patterns such as the appearance of stop-and-go pattern.

Existing 2x2 systems have the general form

∂t ρ + ∂x (ρv) = 0
(3.22) V (ρ)−v
∂t v + v∂x v = T
+ ...
|{z} .
| {z }
reaction terms
relaxation term
 42

Different models can be formulated for the acceleration equation depending on de-

scriptions towards reaction term from different routes. Recall that Example III.7

produced an opposite formation of wave to what was expected due to the convex

part. In 2x2 system, there would be no issue using the convex FD.

3.2.2 A New Model

Microscopic car-following model (2.13) governs the dynamic of each vehicles’ po-

sitions and velocities in a way that it depends on the headway and the relative speed.

First our goal is to embed this idea into the macroscopic level. We assume that all

drivers react in a standard way. To transform vehicle acceleration equation (2.13) to

macroscopic level, we introduce macroscopic variable v(x, t) that is smooth. Then

define v(xn (t), t) = ẋn (t) and v(xn (t)+h(t), t) = ẋn+1 (t) where h(t) = xn+1 (t)−xn (t).

We then expand

h2 h3
ẋn+1 (t) − ẋn (t) = v(xn (t) + h(t), t) − v(xn (t), t) = h∂x v + ∂xx v + ∂xxx v + ...
2 6

and truncate the expansion ignoring cubic terms to yield the so-called non-equilibrium

model

V (ρ) − v h h2
(3.23) ∂t v + v∂x v ≈ + ∂x v + ∂xx v.
T τ 2τ

Here V (ρ) is an optimal velocity function, T is a relaxation time and τ a reaction

time.

We denote by s the Space between vehicles, that is the distance between the front

bumper of one vehicle to the rear bumper of the vehicle ahead of it, and by h the

Headway, that is the distance from front bumper to front bumper. We therefore

have, h = s + lc where lc is the car length. See Fig. 3.16a. We further assume that

there is a minimal safety space between vehicles at which the vehicles will come to a
 43

10
s min = lc

8 2s min = lc

h( )
h
s 4

lc 2

0
0 0.2 0.4 0.6 0.8 1

xn xn+1
(a) Headway and its relations (b) Headway function of density

Figure 3.16: Headway in microscopic (a) and macroscopic (b) description

standstill, and denote this distance by smin (for example, a certain fraction of a car

length). We relate microscopic headway to macroscopic density as follows

smin
(3.24) h(ρ) = + lc .
ρ

Putting (3.23) together with the conservation of car density leads to the viscous 2x2

system

∂t ρ + ∂x (ρv) = 0
(3.25)
h(ρ) V (ρ) − v h2 (ρ)∂xx v
∂t v + (v − )∂x v = + .
τ T 2τ
(3.25) satisfies the anisotropy condition:

h(ρ)
λ1 = v − < λ2 = v.
τ

No information travels faster than the average traffic velocity. Moreover,

∇λ1 · r1 < 0 and ∇λ2 · r2 = 0

where ri are eigenvectors

   
1  1 
r1 =  , r2 =  
 
− h(ρ)
τρ
0

corresponding to λi for i = 1, 2. The first condition, genuinely nonlinear field, satisfies

the property that flow deceleration leads to shock waves (converging characteristics)
 44

and flow acceleration gives rise to rarefaction waves (diverging characteristics). The

second condition, linearly degenerate field, is associated with contact discontinuity,

in the context of traffic flow it describes density change carried with average car

speed.

Non-conservative system

The second equation in (3.25) describes the rules governing car acceleration and

does not reflect a conservation law. Acceleration is the natural quantity to use to

describe the traffic flow dynamics, yet as a consequence, the resulting 2x2 system is

not in conservation form. Clearly, the left-hand side of (3.25) by itself, being in non-

conservation form, is not valid to predict flow patterns involving shock waves. We

also note that the right-hand side of (3.25) includes two terms: a relaxation term and

a diffusion term. We have conducted numerical simulations of (3.25) including both

relaxation and diffusion terms, and compared them to simulations with the relaxation

term alone. In Figure 3.17, we present numerical Riemann solution for three cases:

hyperbolic system itself, denoted by 0 RHS, hyperbolic system with relaxation term,

Relaxed RHS, and hyperbolic system with relaxation and diffusion terms, Relaxed +

Viscous RHS. We concluded that the ’hyperbolic’ solution is substantially modified

by the inclusion of the relaxation term, but once relaxation is included, the diffusion

has only a negligible effect on the solution, and may often be dropped. The model

we propose is therefore

∂t ρ + ∂x (ρv) = 0
(3.26)
h(ρ) V (ρ) − v
∂t v + (v − )∂x v =
τ T
h(ρ)
where τ
is the speed of propagation of small disturbance mostly due to the response

explained by the relative speed between cars.

 45

0.8 0.45
0 RHS 0 RHS
Relaxed RHS 0.4 Relaxed RHS
Relaxed+Viscous RHS Relaxed+Viscous RHS
0.6
0.35

0.3
0.4
0.25

0.2 0.2
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000

x x
1.6
1.5 0 RHS
1.5 Relaxed RHS
Relaxed+Viscous RHS
1.4
v

v
1
1.3
0 RHS
Relaxed RHS 1.2
0.5 Relaxed+Viscous RHS
1.1
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000

x x
(a) Initial ρL > ρR (b) Initial ρL < ρR

Figure 3.17: Numerical study on diffusive effect

3.3 New Multilane Models

In this section we generalize single lane models to multilane, we formulate lane

changing conditions and incorporate inter lane exchange terms. We propose two

types of a multilane model: the first consists of only mass conservations referred to

as Type 1, and the other is in a form of systems, based on the new 2x2 system (3.26),

referred to as Type 2. Other multilane models were studied in [39, 22, 23, 29, 24,

30, 21]. We assume that

1. We consider only one vehicle class whose length is lc . smin can be assumed to

be constant provided one class of vehicle is considered. Note, however, that it

may vary when multi-class plays a role. For example, trucks require larger smin

than automobiles.

2. ρi is the density in lane i for i = 1, · · · , n.

3. V (ρi ) is the optimal velocity in lane i that is a given function of ρi while vi is

the average velocity that is independent of the surrounding ρi .

 46

3.3.1 Lane Change Conditions

We assume that from a driver’s standpoint, lane switch from lane i to a target

lane j (either i − 1 or i + 1) is oftentimes triggered if an extra speed is gained. This

first condition can be formulated by

(3.27) vj > (1 + α)vi

where α ∈ [0, 1] that tends to be relatively small. Note that vi should be regarded as

V (ρi ) in scalar models. (3.27) says that one gains at least α% of current speed vi by

switching to lane j. Once this necessary incentive is met, drivers’ decision to switch

lanes is prompted by checking a sufficient factor: whether or not there is enough

headway (or empty space) in lane j for oneself to merge in. With (3.24), the second

lane change condition now may be expressed as

(3.28) h(ρj ) > hmin + vi τlc .

τlc is a time that takes for changing lanes. The right-hand side in (3.28) expresses

the physical distance to merge for a driver traveling at vi . It consists of the minimal

headway that is required even with vi = 0 plus vi τlc that is a distance to complete

the change. (3.28) says that the current headway in target lane j, h(ρj ), has to be

strictly larger than the required distance. τlc is assumed to be constant in this work,

but this can be generalized more. For example, defensive drivers with speed vi , who

would require larger τlc , are likely to need more distance than aggressive drivers.

This way, τlc could depend on drivers’ characteristics, or could depend on velocity.

This linear assumption in vi can also be found in the reaction thresholds in [29], [28].

When considering three lanes or more, the priority order of lane change in the

middle lane is designated as follows: when a driver needs to speed up, the person
 47

attempts toward the left lane first. If a driver finds a spot, then the driver will

switch the lane. However, if this cannot be conducted because either of lane change

condition is not satisfied, the person then attempts to the right. One would stay its

lane if changing to the right lane is not even allowed.

Closing Remarks: With the two lane change conditions in mind, we now first

generalize LWR model to multilane model. For this, we extend the idea of single

mass conservation into n mass conservations in which each equation represents the

flow in each lane, coupled through source terms representing gain or loss of inter-lanes

due to lane switches.

3.3.2 Multilane Scalar and System Models

Multilane Scalar Model (Type 1)

In Type 1, we essentially have n mass conservations regarding each lane for each

conserved quantity, ρi for i = 1, · · · , n. ρi is then coupled through source terms

representing mass gain or loss by lane changes. It takes the form

(3.29) ∂t ρi + ∂x ρi V (ρi ) = Lossi + Gaini

where

si,i−1 + si,i+1
(3.30a) Lossi = − Rρi
τlc
si−1,i si+1,i
(3.30b) Gaini = Rρi−1 + Rρi+1
τlc τlc

for i = 1, · · · , n. This is a hyperbolic set of balance laws. Here ρ0 = ρn+1 = 0. Loss

and gain involve mass exchanges only when lane change happens. Even when lane

change conditions are met, not all drivers may choose to switch lanes. For this, we

use R to denote the fraction of drivers that may choose to switch. τlc is the time for
 48

lane 1 ρ1 , v1

lane 2 S1,2 ρ2 , v2
S3,2

lane 3 ρ3 , v3

Figure 3.18: Diagram illustrating lane changes in three lanes

changing lanes and si,j is

 V (ρj ) > (1 + α)V (ρi )

n 1 if
si,j = h(ρj ) > hmin + V (ρi )τlc

0 otherwise

si,j is an indicator function of lane-changing. si,j may be expressed concisely as

n  o
max 0, min sign(V (ρj ) − (1 + α)V (ρi ), sign(h(ρj ) − hmin − V (ρi )τlc ) .

See Figure 3.18 for a three lanes schematic. In the simulations in this work, R is

constant for simplicity. Note that s0,1 = s1,0 = sn+1,n = sn,n+1 = 0. Switching on

and off the source terms in (3.29), we are capable of differentiating the effect of lane

changes. Also the fact that summing up over all lanes results in zero source tells

that there is a conservation over all lanes but not in-between lanes because of lane

switches. We also deal with effects of parameters (lc , smin , τlc , etc) when we visit

section 4.2. Three-lane example is following

s1,2 s2,1
∂t ρ1 + ∂x ρ1 V (ρ1 ) = − Rρ1 + Rρ2
τlc τlc

s1,2 s2,1 + s2,3 s3,2
(3.31) ∂t ρ2 + ∂x ρ2 V (ρ2 ) = Rρ1 − Rρ2 + Rρ3
τlc τlc τlc

s3,2 s2,3
∂t ρ3 + ∂x ρ3 V (ρ3 ) = − Rρ3 + Rρ2 .
τlc τlc
Recall that V (ρi ) in Type 1 locally depends on nearby ρi . Oftentimes, the average

speed in real traffic does not reflect ρ explicitly. We now build up the Type 2

multilane model based on the 2x2 system that considers this element.
 49

Multilane System Model (Type 2)

Including the lane changing wisdom in the same way, we extend the single-lane

2x2 system to multilane:

(3.32a) ∂t ρi + ∂x ρi vi = Lossi + Gaini
h(ρi ) V (ρi ) − vi
∂t vi + (vi − )∂x vi = +
τ T
si−1,i ρi−1 si+1,i ρi+1
(3.32b) R (vi−1 − vi ) + R (vi+1 − vi ).
τlc ρi τlc ρi

v0 = vn+1 = 0. In (3.32b), there are reaction terms associated with intra-lane

acceleration: one is acceleration towards an equilibrium speed and the other is an

approximate acceleration towards the speed of car ahead. And the exchange terms

involving (vi −vj ) are inter-lane acceleration towards the speed in the target lane due

to the lane changing. If a car in the other lane with a slower speed moves into the

target lane, it causes some deceleration. This way v’s in different lanes are allowed

to communicate each other through this term as far as lane change is concerned. For

a three-lane example, this reads

(3.33)
s1,2 s2,1
∂t ρ1 + ∂x (ρ1 v1 ) = − Rρ1 + Rρ2
τlc τlc
h(ρ1 ) V (ρ1 ) − v1 s2,1 ρ2
∂t v1 + (v1 − )∂x v1 = + R (v2 − v1 )
τ T τlc ρ1
s1,2 s2,1 + s2,3 s3,2
∂t ρ2 + ∂x (ρ2 v2 ) = Rρ1 − Rρ2 + Rρ3
τlc τlc τlc
h(ρ2 ) V (ρ2 ) − v2 s1,2 ρ1 s3,2 ρ3
∂t v2 + (v2 − )∂x v2 = + R (v1 − v2 ) + R (v3 − v2 )
τ T τlc ρ2 τlc ρ2
s2,3 s3,2
∂t ρ3 + ∂x (ρ3 v3 ) = Rρ2 − Rρ3
τlc τlc
h(ρ3 ) V (ρ3 ) − v3 s2,3 ρ2
∂t v3 + (v3 − )∂x v3 = + R (v2 − v3 ).
τ T τlc ρ3
Multilane model with more than three lanes can be formulated as desired.
 50

Closing Remarks: Multilane model especially benefit from understanding traffic

bottleneck. Recall that chronic traffic bottleneck occurs most of the time in the

cases of lane-closing, accidents, ramps etc, where it seems that the multilane model

should play a role. Studying lane-by-lane traffic flow capturing lane switch effects

also provides us a more accurate understanding of overall traffic behavior.

 CHAPTER IV

Numerical Method and Result

Chapter IV includes numerical methods and PDE-based computations built on

Chapter III.

4.1 Numerical Method Revisit

Upwind Scheme Revisit

Consider 2x2 model (3.26) in matrix form

∂t W + A(W )∂x W = S(W )

where
     
 ρ   v ρ 0
W =   , A(W ) =   , and S(W ) =  .
  
v 0 v − h(ρ)
τ
V (ρ)−v
T

The eigenstructure is given by the matrices

   
h(ρ)
 1 1   v− τ 0 
R=  and Λ =  
h(ρ)
− τρ 0 0 v

so that A(W )R = RΛ. We use a Roe-type upwind scheme with the source term by

projecting S(W ) onto the eigenvectors of A(W ) as illustrated in section 1.1.2. A(W )

is linearized about averages at the cell interface

ρi+1 + ρi vi+1 + vi
ρ̄ = , v̄ = .
2 2

51
 52

The source term is approximated by

Si+1 + Si
S̄ = .
2

Wave strengths are then

τ ρ̄∆v
α1 = −h(ρ̄)
, α2 = ∆ρ − α1
τ ρ̄S2
β1 = −h(ρ̄)
, β2 = −β1

where ∆( ) = ( )i+1 − ( )i . This scheme can be also applicable to Type 1 and 2

models. For example, the three lane model of Type 2 (3.33) has eigenstructure
 
 1 1 0 0 0 0 
 
 h(ρ1 ) 
 − 0 0 0 0 0 
 τ ρ1 
 
 0 0 1 1 0 0 
 
R=



 0
 0 − h(ρ 2)
τ ρ2
0 0 0 

 
 
 0 0 0 0 1 1 
 
 
0 0 0 0 − h(ρ 3)
τ ρ3
0

and
 h(ρ1 ) h(ρ2 ) h(ρ3 ) 
Λ = diag v1 − , v1 , v2 − , v2 , v3 − , v3 .
τ τ τ

Wave strengths are

τ ρ̄1 ∆v1 τ ρ̄2 ∆v2

α1 = −h(ρ̄1 )
, α2 = ∆ρ1 − α1 , α3 = −h(ρ̄2 )

τ ρ̄3 ∆v3
α4 = ∆ρ2 − α3 , α5 = −h(ρ̄3 )
, α6 = ∆ρ3 − α5

τ ρ̄1 S̄2 τ ρ̄2 S̄4

β1 = −h(ρ̄1 )
, β2 = −β1 , β2 = −h(ρ̄2 )

τ ρ̄3 S̄6
β4 = −β3 , β5 = −h(ρ̄3 )
, β6 = −β5 .

Example results displayed in the following are by this first-order scheme.

 53

4.2 Numerical Examples

We first present the stop-and-go flow with (3.26), followed by a school zone ex-

ample. In examples, we assume that lc = 1 and τ = 1. Note that different FDs are

chosen in different examples.

Example IV.1.

Stop-and-go flow by System (3.26): we carried out the stop-and-go phenomenon

by a microscopic model in section 2.2. In this example, we show that the 2x2 model

(3.26) regenerates the unpredictable appearance of a traffic congestion comparable

to stop-and-go flow. For this, we locally perturb uniform flow at three different

places with the Gaussian pulse. A little later it amplifies the perturbations as in

the ODE system, and then in the long run, solution remains periodic, capturing

break-acceleration-break flow, see Figure 4.1. In this simulation,

V (ρ) = Vmax 1/(1 + e(ρ−0.25)/0.06 ) − 3.72 · 10−6

was used as in [27] by Kerner et al. . smin = 1/3, and T = 10 were used. We

found that the number of initial perturbations corresponds to the number of break-

acceleration flows, and the magnitude of the perturbation decides how fast the break-

acceleration flow reaches. Larger perturbation, faster the flow stabilizes. Note that

the initial condition needs to be chosen so that

h(ρ(x, 0))
(4.1) v(x, 0) − ≤ v(x, 0) + ρ(x, 0)V 0 (ρ(x, 0)) ≤ v(x, 0)
lc

is satisfied. See [27, 25] for details. 

 54

1
t= 0
0.8 t= 101.1284
t= 850
0.6

0.4

0.2

0
0 100 200 300 400 500 600 700 800 900 1000

3
t= 0
t= 101.1284
t= 850
2

0
0 100 200 300 400 500 600 700 800 900 1000

(a) Profiles of ρ and v at t = 0, 100, 850 (b) Time evolution of ρ and v

Figure 4.1: Instability of PDE system of (3.26)

Example IV.2.

School Zone Area by System (3.26): we simulate a school zone area that is the

road that has a portion of variable speed limit, see the geometric structure in Fig-

ure 4.2a. We model this by spatially varying the speed limit in V (ρ) in (3.26). In

this simulation, Greenshields’s model

V (ρ) = Vmax (1 − ρ)

was used and Vmax is reduced by 1/3 in the school zone. smin = 1 and T = 1/2

were used. What cars experience is slowing down as entering the school zone, which

creates car builds-ups going backwards, and then speeding up once they leave the

school zone, see Figure 4.2b and 4.2c. The newly-developed V (ρ) that has three

phases was used and results are computed in Figure 4.3. This three phase V (ρ)

added more diffusive effect on the results.

This can be applicable to cases: the area where locally bad weather is affected,

or the road that has a poor condition such as unpaved road so speed limit has to be

reduced [40, 35]. 

 55

SCHOOL
REDUCE Vmax =Vmax(x)
SPEED
AHEAD

Vmax >0
Vmax(x)
x
(a) Schematic for single-lane traffic flow through school zone

1
t= 0
0.8 t= 200

0.6

0.4

0.2

0
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000

x x
3

2
v
v

1
t= 0
t= 200
0
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000

x x
(b) Profiles of ρ and v at t = 0, 200 (c) Time evolution of ρ and v

Figure 4.2: Example IV.2 results by (3.26)

0.8
t= 0
t= 200
0.6

0.4

0.2

0
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000

x x
3

2
v
v

1
t= 0
t= 200
0
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000

x x
(a) Profiles of ρ and v at t = 0, 200 (b) Time evolution of ρ and v

Figure 4.3: Example IV.2 results by (3.26) with three-phase V (ρ)

 56

We now present examples of both Type 1 and 2 models. We assume that lc = 1

and τlc = 1 for Type 1 and τ = 1 for Type 2.

Example IV.3.

Lane Closure by (3.31): in this example, we consider the three-lane traffic flows

where a lane drop is included. This circumstance may be attributed to a non-uniform

geometry of roadways or an accident which necessitates a lane block and therefore

lane switches. It is a trouble spot where there is chronic congestion or bottleneck, so

it has been attracted attention to many researchers, see [41, 28, 42]. See Figure 4.4a

for the road designed for this simulation. This lane-drop is then modeled by reducing

the speed limit in lane 1 gradually to zero over the taper from three lanes to two. In

this simulation,

V (ρ) = Vmax (1 − ρ),

and Vmax in lane 2 and 3 are assigned to be 10% less than lane 1. Different speed

limits are imposed as to the very left lane is reserved for faster-moving or passing

vehicles under traffic laws. smin = 2, R = 0.5, and α = 0.01. Putting initial density

and speed uniformly in lane 2 and 3, and linearly decreasing density and speed to

zero over the road taper in lane 1, reports the result in Figure 4.4b and its time

evolution in Figure 4.4c. Because lane 1 ends, lane switch occurs to the right lanes,

and this leads to traffic bottlenecks. We found that different initial data makes the

results ended differently. Lighter initial density less than about ρ = 0.3 does not

result in traffic bottleneck at all, as opposed to traffic bottleneck is observed over all

lanes with denser initial density. Our results show that shock waves in ρ in all lanes

propagates backwards with similar speed but little different strength. Larger smin ,

closer the shock strength in ρi for all i. This is because it leaves more minimal space
 57

to allow cars to merge in, and this reduces the difference in shock strengths in each

lanes. This also can be observed in (3.28) by rearranging,

1 
smin − 1 > V (ρi ).
ρj

As smin gets larger, more chances of lane change condition to be satisfied and it

affects the results in speed and strength of shock.

As a reference computation, we add a comparison with the scalar lane-aggregate

model (3.18) by constructing n(x) according to the road set up, see Figure 4.5a.

It illustrates the accuracy and robustness of the model (3.31) especially the source

terms, see Figure 4.5b. Propagation of shock speed and strength are not exactly

same because the scalar lane-aggregate model (3.18) and Type 1 multilane model

(3.31) in its summation form do not fully agree. 

 58

Vmax(x)
LANE ENDS
x
MERGE
RIGHT

Vmax =0

Vmax =Vmax(x) in lane 1

(a) Schematic for multilane traffic flow through lane closure

1
lane1 lane1
1.2
0.8 lane2 lane2
lane3 1 lane3

0.6 0.8

0.4 0.6

0.4
0.2
0.2

0 0
0 20 40 60 80 100 0 20 40 60 80 100

(b) Profiles of ρ and ρV (ρ) in each lane at t = 70

(c) Time evolution of ρ and ρv

Figure 4.4: Example IV.3 results by (3.31)

 59

n(x) =3
n(x) =2
n(x)
x
(a) n(x) in (3.18) for the lane closure example

2.5 2.8

2.6
2
2.4
1.5
2.2

1 2

1.8
0.5
0 20 40 60 80 100 0 20 40 60 80 100

(b) Comparison of cumulative ρ and ρV (ρ) at t = 70 in lane drop simulation solved

by (3.18), shown in dashed line, and by (3.31), shown in solid line. In this simulation
speed limit in all lanes is assigned to be identical.

Figure 4.5: Example IV.3 results by (3.31) with (3.18)

Example IV.4.

Variable Speed Limit in Multilane by (3.33): in this example, we compute similar

scenario of Example IV.2 but on a multilane road. Imagine a case where a road work

is going on the shoulder of one side, which implicitly have drivers slow down to

work-zone speed limit when workers are present. See Figure 4.6a for an example of

three-lane roadways that has a construction area on a shoulder of lane 1. Traffic flow

in this environment can be computed by changing speed limit in V (ρ) in (3.33). In

this computation, we present a comparison with and without lane changing effect.

When the lane changing is not allowed, only intra-lane acceleration plays a role with

no vehicle exchange and consequently no inter-lane acceleration is involved. See

Figure 4.6b. In this simulation,

V (ρ) = Vmax (1 − ρ)

was used. smin = 1, T = 1/3, τlc = 2, α = 0.1, and R = 0. Although vehicles

gradually move downstream in lane 1 causing rarefaction wave, vehicles are backed
 60

up backwards forming relaxed shock-like wave due to the reduction in speed limit

ahead. Meanwhile, with homogenous initial conditions, ρ and v keep uniform flow

in the other lanes due to no interruption. See Figure 4.6c.

When the lane change is allowed, vehicles move away from lane 1 to lane 2 and

lane 3, and eventually stabilizes. See Figure 4.7. We put results with and without

the lane change effect together. The traffic queue in lane 1 when lane changing was

not allowed is found to be alleviated because of lane changes, see Figure 4.8. We plot

a trajectory of a vehicle in lane 1 calculated by v1 computation, when lane changing

is allowed. This is then compared by the trajectory when the lane changing is not

allowed, see Figure 4.9. This provides it takes 29.60 unit time faster to get through

the congestion when lane changing is allowed. 

 61

Vmax(x)
SHOULDER CONSTRUCTION x

Vmax >0

Vmax =Vmax(x) in lane 1

(a) Schematic for Example IV.4

0.8 1.5
t=0 t=0
0.7 t=500 t=500

0.6

0.5 1

0.4

0.3

0.2 0.5
0 50 100 150 200 250 300 0 50 100 150 200 250 300

(b) Profiles of ρ and v at t = 0, 500 in lane 1

0 50 100 150 200 250 300 0 50 100 150 200 250 300

0 50 100 150 200 250 300 0 50 100 150 200 250 300

(c) Time evolution of ρ and v

Figure 4.6: Example IV.4 results by (3.33). In this simulation, R = 0.

 62

0.9 1.6
lane1 lane1
0.8 lane2 1.4 lane2
lane3 lane3
0.7
1.2
0.6
1
0.5

0.4 0.8

0.3 0.6
0 50 100 150 200 250 300 0 50 100 150 200 250 300

(a) Profiles of ρ and v at t = 500

0 50 100 150 200 250 300 0 50 100 150 200 250 300

0 50 100 150 200 250 300 0 50 100 150 200 250 300

0 50 100 150 200 250 300 0 50 100 150 200 250 300

(b) Time evolution of ρ and v

Figure 4.7: Example IV.4 results by (3.33). In this simulation, R = 0.3.

 63

0.9 1.5
Lane change ON Lane change ON
0.8 Lane change OFF Lane change OFF

0.7

0.6 1

0.5

0.4

0.3 0.5
0 50 100 150 200 250 300 0 50 100 150 200 250 300

Figure 4.8: Example IV.4 results in lane 1 superimposed with and without effect of lane changes.

300

250

200

150

100

50 Lane change ON
Lane change OFF
0
0 50 100 150 200 250 300

Figure 4.9: Trajectories of vehicles in lane 1 according to the computation of Example IV.4

Closing Remarks: Without traffic field data, it is hard to quantitatively conclude

that results reflect the real traffic. However, it demonstrates that models developed

in this dissertation are practicable to simulate interesting traffic examples. Moreover,

all results show good qualitative agreement with daily life traffic experience.
 CHAPTER V

Conclusion and Future Work

The field of traffic flow theory, computation, and its application is increasingly rec-

ognized as the era of autonomous vehicles approaches. Full understanding of mixed

flow of human and autonomous vehicles requires good understanding of human-

driving models, which are main theme of this dissertation. Furthermore, realistic

and more accurate models make predictions better, and reduce the amount of money

spent on gas, driving times and drivers’ frustration levels because it improves the

overall flow mobility. All above has been a motivation of this work.

We summarize the main contributions and new results of this dissertation, and

discuss potential future applications. Data-driven fundamental relation of density

and mean velocity was proposed, which is essential for macroscopic traffic flow models

to be effective. We developed an acceleration equation that reflects microscopic

characteristics, by making a link between microscopic and macroscopic variables,

headway and density. The models were then generalized to multilane traffic models

incorporating several features for lane changing: comparison of speeds and headways

between lanes, mass exchanges, and inter-lane acceleration effects due to the car

exchanges. We presented various numerical simulations, implemented by the Roe-

type upwind scheme, including a lane reduction example from three to two lanes and

64
 65

variable speed limits area in multilane.

This work opens up possibilities for other traffic challenges such as a network

of roads, a road with signals, and roundabout problem. This work can be useful to

control the flow. This is because prediction and estimation of throughput of the road

can be used to enhance to maximize the flux without delaying vehicles on a road. This

study may also deduce the lane-wise stop-and-go flow in multilane. In addition to

this, multilane models can combine with multi-class to better understand the complex

traffic phenomenon. Connected and autonomous vehicles will be on highways in

a very near future. For this, modeling multilane with multi-class of platoons of

connected and autonomous vehicles and human driving, which may require different

lane-changing maneuver, is inevitable. Lastly, one ultimate goal may be to compare

the multilane model to field data, which remains as future work.

APPENDIX

66
 67

APPENDIX A

Data-driven three-phase FD

Using US 101 data, we optimize m and α in (3.6) by moving a window over the entire

time and space. We obtain m and α by solving the optimization problem:

min k ln v − m ln ρ − ln αk22
m,α

with a constraint m ≤ 0. The result of the optimized m and α of US 101 data are

shown in Figure A.1a and A.1c. Patterns in m and α do align with v. m and α

are inversely related to ρ. Visual representation of the two parameters against ρ in

Figure A.1b and A.1d supports the generalization form in (3.7).

Figure A.1: Optimized m and α, and its relation with ρ

BIBLIOGRAPHY

68
 69

BIBLIOGRAPHY

[1] A. Aw and M. Rascle. Resurrection of ”Second Order” Models of Traffic Flow. SIAM Journal
on Applied Mathematics, 60(3):916–938, January 2000.
[2] M. Bando, K. Hasebe, A. Nakayama, A. Shibata, and Y. Sugiyama. Dynamical model of traffic
congestion and numerical simulation. Physical Review E, 51(2):1035–1042, February 1995.
[3] Masako Bando, Katsuya Hasebe, Ken Nakanishi, and Akihiro Nakayama. Analysis of optimal
velocity model with explicit delay. Physical Review E, 58(5):5429, 1998.
[4] Masako Bando, Katsuya Hasebe, Ken Nakanishi, Akihiro Nakayama, Akihiro Shibata, and
Yūki Sugiyama. Phenomenological Study of Dynamical Model of Traffic Flow. Journal de
Physique I, 5(11):1389–1399, November 1995.
[5] Sylvie Benzoni-Gavage and Rinaldo M. Colombo. An $n$-populations model for traffic flow.
European Journal of Applied Mathematics, 14(5):587–612, October 2003.
[6] Peter Berg, Anthony Mason, and Andrew Woods. Continuum approach to car-following mod-
els. Physical Review E, 61(2):1056–1066, February 2000.
[7] Sten Bexelius. An extended model for car-following. Transportation Research, 2(1):13–21,
March 1968.
[8] Inc. Cambridge Systematics. NGSIM U.S. 101 Data Analysis.
https://data.transportation.gov/Automobiles/Next-Generation-Simulation-NGSIM-Vehicle-
Trajector/8ect-6jqj, December 2005.
[9] J. M. Del Castillo and F. G. Benı́tez. On the functional form of the speed-density rela-
tionship—I: General theory. Transportation Research Part B: Methodological, 29(5):373–389,
October 1995.
[10] Robert E. Chandler, Robert Herman, and Elliott W. Montroll. Traffic Dynamics: Studies in
Car Following. Operations Research, 6(2):165–184, April 1958.
[11] Kang-Ching Chu, Li Yang, R. Saigal, and K. Saitou. Validation of stochastic traffic flow
model with microscopic traffic simulation. In 2011 IEEE Conference on Automation Science
and Engineering (CASE), pages 672–677, August 2011.
[12] R. M. Colombo. A 2x2 hyperbolic traffic flow model. Mathematical and Computer Modelling,
35(5):683–688, March 2002.
[13] R. Courant, K. Friedrichs, and H. Lewy. On the Partial Difference Equations of Mathematical
Physics. IBM Journal of Research and Development, 11(2):215–234, March 1967.
[14] Carlos F. Daganzo. Requiem for second-order fluid approximations of traffic flow. Transporta-
tion Research Part B: Methodological, 29(4):277–286, August 1995.
[15] J. L. Drake and Joseph L. Schofer. A Statistical Analysis of Speed-Density Hypotheses. High-
way Research Record 154, pages 53–87, 1966.
 70

[16] Denos C. Gazis, Robert Herman, and Renfrey B. Potts. Car-Following Theory of Steady-State
Traffic Flow. Operations Research, 7(4):499–505, August 1959.
[17] Denos C. Gazis, Robert Herman, and Richard W. Rothery. Nonlinear Follow-the-Leader Mod-
els of Traffic Flow. Operations Research, 9(4):545–567, August 1961.

[18] S. K. Godunov. A difference method for numerical calculation of discontinuous solutions of

the equations of hydrodynamics. page 37, 1959.
[19] Harold Greenberg. An Analysis of Traffic Flow. Operations Research, 7(1):79–85, February
1959.

[20] B.d. Greenshields, J.r. Bibbins, W.s. Channing, and H.h. Miller. A study of traffic capacity.
Highway Research Board proceedings, 1935, 1935.
[21] Arvind Kumar Gupta and Sapna Sharma. Analysis of the wave properties of a new two-lane
continuum model with the coupling effect. Chinese Physics B, 21(1):015201, 2012.
[22] Dirk Helbing and Andreas Greiner. Modeling and simulation of multilane traffic flow. Physical
Review E, 55(5):5498–5508, May 1997.
[23] Edward N. Holland and Andrew W. Woods. A continuum model for the dispersion of traffic
on two-lane roads. Transportation Research Part B: Methodological, 31(6):473–485, November
1997.

[24] Haijun Huang, Tieqiao Tang, and Ziyou Gao. Continuum modeling for two-lane traffic flow.
Acta Mechanica Sinica, 22(2):131–137, April 2006.
[25] Rui Jiang, Qing-Song Wu, and Zuo-Jin Zhu. A new continuum model for traffic flow and
numerical tests. Transportation Research Part B: Methodological, 36(5):405–419, June 2002.
[26] Rui Jiang, Qingsong Wu, and Zuojin Zhu. Full velocity difference model for a car-following
theory. Physical Review E, 64(1):017101, June 2001.
[27] Boris S. Kerner and Peter Konhäuser. Cluster effect in initially homogeneous traffic flow.
Physical Review E, 48(4):R2335, 1993.
[28] A. Klar and R. Wegener. A Hierarchy of Models for Multilane Vehicular Traffic I: Modeling.
SIAM Journal on Applied Mathematics, 59(3):983–1001, January 1998.
[29] Axel Klar and Raimund Wegener. A hierarchy of models for multilane vehicular traffic II:
Numerical investigations. SIAM Journal on Applied Mathematics, 59(3):1002–1011, 1998.
[30] Jorge A. Laval and Carlos F. Daganzo. Lane-changing in traffic streams. Transportation
Research Part B: Methodological, 40(3):251–264, March 2006.

[31] Peter Lax and Burton Wendroff. Systems of conservation laws. Communications on Pure and
Applied Mathematics, 13(2):217–237, May 1960.
[32] Peter D. Lax. Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock
Waves. SIAM, January 1973.

[33] Chin Jian Leo and Robert L. Pretty. Numerical simulation of macroscopic continuum traffic
models. Transportation Research Part B: Methodological, 26(3):207–220, June 1992.
[34] Randall J LeVeque. Some Traffic Flow Models Illustrating Interesting Hyperbolic Behavior.
page 11, 2001.

[35] Randall J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University
Press, August 2002. Google-Books-ID: O ZjpMSZiwoC.
 71

[36] Jia Li, Qian-Yong Chen, Haizhong Wang, and Daiheng Ni. Analysis of LWR model with fun-
damental diagram subject to uncertainties. Transportmetrica, 8(6):387–405, November 2012.
[37] M. J. Lighthill and G. B. Whitham. 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, 1955.
[38] Xiao-Yun Lu, P Varaiya, and Roberto Horowitz. Fundamental Diagram modeling and analysis
based NGSIM data. In 12th IFAC Symposium on Control in Transportation System, pages
367–374, 2009.
[39] Panos G. Michalopoulos, Dimitrios E. Beskos, and Yasuji Yamauchi. Multilane traffic flow
dynamics: Some macroscopic considerations. Transportation Research Part B: Methodological,
18(4):377–395, August 1984.
[40] S. Mochon. An analysis of the traffic on highways with changing surface conditions. Mathe-
matical Modelling, 9(1):1–11, January 1987.

[41] P. K. Munjal, Yuan-Shih Hsu, and R. L. Lawrence. Analysis and validation of lane-drop effects
on multi-lane freeways. Transportation Research, 5(4):257–266, December 1971.
[42] K. Nassab, M. Schreckenberg, A. Boulmakoul, and S. Ouaskit. Effect of the lane reduction in
the cellular automata models applied to the two-lane traffic. Physica A: Statistical Mechanics
and its Applications, 369(2):841–852, September 2006.

[43] Gábor Orosz, R. Eddie Wilson, and Bernd Krauskopf. Global bifurcation investigation of an
optimal velocity traffic model with driver reaction time. Physical Review E, 70(2), August
2004.
[44] Gábor Orosz, R. Eddie Wilson, Róbert Szalai, and Gábor Stépán. Exciting traffic jams:
Nonlinear phenomena behind traffic jam formation on highways. Physical Review E, 80(4),
October 2009.
[45] H. OZAKI. Reaction and Anticipation in the Car-Following Behavior. Proc.of 12th Interna-
tional Symposium on Theory of Traffic Flow and Transportation, pages 349–366, 1993.
[46] Harold J. Payne. MODELS OF FREEWAY TRAFFIC AND CONTROL. MATHEMATICAL
MODELS OF PUBLIC SYSTEMS, 1971.
[47] Paul I. Richards. Shock Waves on the Highway. Operations Research, 4(1):42–51, February
1956.
[48] P. L Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal
of Computational Physics, 43(2):357–372, October 1981.

[49] P. L Roe. Fluctuation and signals - a framework for numerical evolution problems. Numerical
Methods for Fluid Dynamics, (Academic Press):219–257, 1982.
[50] P. L. Roe. Upwind differencing schemes for hyperbolic conservation laws with source terms.
In Claude Carasso, Denis Serre, and Pierre-Arnaud Raviart, editors, Nonlinear Hyperbolic
Problems, Lecture Notes in Mathematics, pages 41–51. Springer Berlin Heidelberg, 1987.
[51] Joel Smoller. Shock Waves and Reaction—Diffusion Equations. Springer Science & Business
Media, December 2012. Google-Books-ID: lOIlBQAAQBAJ.
[52] Yuki Sugiyama, Minoru Fukui, Macoto Kikuchi, Katsuya Hasebe, Akihiro Nakayama, Kat-
suhiro Nishinari, Shin-ichi Tadaki, and Satoshi Yukawa. Traffic jams without bottle-
necks—experimental evidence for the physical mechanism of the formation of a jam. New
Journal of Physics, 10(3):033001, 2008.
 72

[53] Eleuterio F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Prac-
tical Introduction. Springer Science & Business Media, April 2013. Google-Books-ID: zkLt-
CAAAQBAJ.
[54] R.T. Underwood. Speed, volume, and density relationships: Quality and theory of traffic flow.
Yale Bureau of Highway Traffic, New Haven, CT, pages 141–188, 1961.
[55] G. B. Whitham. Linear and nonlinear waves. John Wiley &amp; Sons, 1974.
[56] G. C. K Wong and S. C Wong. A multi-class traffic flow model – an extension of LWR model
with heterogeneous drivers. Transportation Research Part A: Policy and Practice, 36(9):827–
841, November 2002.

[57] Li Yang, Romesh Saigal, Chih-Peng Chu, and Yat-Wah Wan. Stochastic model for traffic flow
prediction and its validation. In Transportation Research Board 90th Annual Meeting, 2011.
[58] H. Zhang and W. Jin. Kinematic Wave Traffic Flow Model for Mixed Traffic. Transportation
Research Record: Journal of the Transportation Research Board, 1802:197–204, January 2002.

[59] H. M. Zhang. A non-equilibrium traffic model devoid of gas-like behavior. Transportation

Research Part B: Methodological, 36(3):275–290, March 2002.