Traffic-flow models: analysis, estimation and
control
by
Chao Xia
B.S., University of Science and Technology of China, Hefei, China, 2012
M.Sc., Brown University, Providence, RI, 2013
A dissertation submitted in partial fulfillment of the
requirements for the degree of Doctor of Philosophy
in the Division of Applied Mathematics at Brown University
PROVIDENCE, RHODE ISLAND
May 2017
c Copyright 2017 by Chao Xia
Abstract of “ Traffic-flow models: analysis, estimation and control ” by Chao Xia,
Ph.D., Brown University, May 2017
My research focuses on the study of traffic-flow models and their applications. Macro-
scopic and microscopic models are the two main approaches: macroscopic models
describe the spatial quantities of traffic, such as density, velocity and flux; while mi-
croscopic models simulate the behavior of individual cars based on their interaction.
For the macroscopic model, we study the Lighthill-Whitham equation, and ac-
count for multiple traffic scenarios by modifying the original Lighthill-Whitham equa-
tion. We also study several microscopic car-following models: the optimal velocity
model, the full velocity difference model, the modified GHR model and the intelligent
driver model. The main research work include:
• Investigate the collision behavior of the microscopic car-following model. We
theoretically prove the collision-free property of several car-following models
through fast-slow system technique, and also carry out numerical simulations
to provide a valid reference to the dynamics of traffic collisions.
• Apply data assimilation technique (ensemble Kalman filter and particle filter)
to estimating the traffic states and uncertain parameters. An augmented ap-
proach is proposed to simultaneously assimilate the Eulerian sensor data and
Lagrangian GPS data.
• Study the phenomenon of capacity discharge in the lane-drop scenario. Macro-
scopically, we model the lane-drop scenario with inhomogeneous Lighthill-
Whitham equation, and then proposed two controlling strategies to guide ve-
hicles smoothly through the bottleneck: (1) change driving habit through fun-
damental diagram; (2) merge vehicles in advance through virtual lane usage.
This dissertation by Chao Xia is accepted in its present form
by the Division of Applied Mathematics as satisfying the
dissertation requirement for the degree of Doctor of Philosophy.
Date
Bj¨orn Sandstede, Ph.D., Advisor
Recommended to the Graduate Council
Date
Constantine Dafermos, Ph.D., Reader
Date
Matthew Harrison, Ph.D., Reader
Approved by the Graduate Council
Date
Andrew G. Campbell, Dean of the Graduate School
iii
Vitae
Education
Brown University, Providence, RI, 2012-2017
Ph.D. in Applied Mathematics, expected May 2017
M.Sc. in Applied Mathematics, 2013
University of Science and Technology of China, Hefei, China, 2008-2012
B.S. in Statistics, School for Gifted Young, 2012
Publications
C Xia and B Sandstede. Traffic control in lane-drop scenarios. In preparation, 2017.
C Xia and B Sandstede. A study of collision behavior of car-following models.
Preprint, 2017.
C Xia, C Cochrane, J DeGuire, G Fan, E Holmes, M McGuirl, P Murphy, J Palmer,
P Carter, L Slivinski, and B Sandstede. Assimilating Eulerian and Lagrangian data
iv
in traffic-flow models. Physica D: Nonlinear Phenonmena. Under review, 2016.
P Wu, C Xia, F Liu, Y Wu and Y He. An integrated water strategy based on the
current circumstances in China. Applied Mathematical Modeling 40 (2016) 8108-
8124.
K Judd, T Stemler and C Xia. Paleoclimate variability and frequency selection by
stochastic resonance in Pitchfork Bifurcations. Physical Review E. Submitted.
Selected Presentations
”Traffic flow: estimation and control.” Graduate Student Seminar of the Division of
Applied Mathematics, Brown University, Providence, RI, November 2016.
”Introduction to data Assimilation for traffic estimation.” Oberwolfach Seminar:
Data Assimilation: The Mathematics of connecting Dynamical System to Data,
Oberwolfach, Germany, May 2016.
”Lagrangian data assimilation in traffic-flow models.” (Poster) IPAM Workshop of
Traffic Estimation in UCLA, Los Angeles, CA, Octorber 2015.
”Data assimilation for traffic states and parameters estimation.” (Poster) SIAM
Conference on Applications of Dynamical Systems, Snowbird, UT, May 2015.
v
Honors and Awards
2016 NSF Supplemental Funding (summer internship)
2012 Outstanding Undergraduate Research Project, University of Science and Tech-
nology of China
2011 Third Prize for China Undergraduate Mathematical Contest Modeling
2011 Certificate of Excellence from University of Western Australia, University of
Western Australia
Professional Experience
Intern, JP. Morgan Chase & Co., 2016 Summer
Intern, Liberty Mutual Insurance, 2015 Summer
Visitor, University of Western Australia, 2011 Summer
Teaching Experience
Teaching Assistant, Brown University
APMA1650 (Statistical Inference), Fall 2013
vi
APMA0650 (Essential Statistics), Spring 2014
Professional Associations
American Mathematical Society (AMS)
Association for Women in Mathematics (AWM)
Society for Industrial and Applied Mathematics (SIAM)
vii
Dedication
To my loving father, who nurtured my curiosity in mathematics when I was in
primary school, and supported my pursuit in mathematics in higher education.
viii
Acknowledgements
I would first like to thank my advisor Bj¨orn Sandstede, for his mentorship in my thesis
topic and his support in my professional development. He has his own philosophies
of supervising students, and his responsibility, dedication and enthusiasm in research
work has a positive impact on my experience as a graduate student. Aside from being
an advisor in my research work, he is willing to provide resources and opportunities
to develop my expertise in industrial job application. It is truly my honor to have
Bj¨orn as my advisor in my graduate study.
I would also like to thank Professor C. Dafermos and Professor M. Harrison for
donating their time to serve on my thesis committee. Thanks to the valuable com-
ments and feedback for my thesis. I would like to thank the collaboration in the
Summer RTG Program on ”Integrating Dynamics and Stochastics”, which enlight-
ened my interest in the traffic-flow models.
I would also like to thank the team members of the research group for the infor-
mative group meetings, and helpful discussions. In particular, I would like to thank
Al, Paul and Veronica for sharing a lovely office, varied opinions and study resources.
Finally, I want to express my thanks to my family: my husband Wenzhe, my
parents Linyan and Qiuqin, for their continued care and love. Special thanks to
the companion and encouragement from my husband in my graduate study, and his
sharing of my happiness and sadness.
ix
Contents
Vitae iv
Dedication viii
Acknowledgments ix
1 Introduction 1
1.1 Motivation and preliminaries . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Thesis objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Review: macroscopic traffic models 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Lighthill-Whitham-Richard model . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Fundamental diagram . . . . . . . . . . . . . . . . . . . . . . 11
2.2.2 Characteristic curve . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Traffic scenarios based on LWR model . . . . . . . . . . . . . . . . . 17
2.3.1 Normal traffic . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3.2 On-/off-ramps . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.3 Traffic lights . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3.4 Construction zone . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.5 Traveling bottleneck . . . . . . . . . . . . . . . . . . . . . . . 22
3 Review: microscopic traffic models 24
3.1 Introduction to microscopic traffic models . . . . . . . . . . . . . . . 25
3.2 Car-following models . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2.1 The optimal velocity model . . . . . . . . . . . . . . . . . . . 28
3.2.2 The velocity difference model . . . . . . . . . . . . . . . . . . 29
3.2.3 The modified GHR model . . . . . . . . . . . . . . . . . . . . 30
3.2.4 The intelligent driver model . . . . . . . . . . . . . . . . . . . 31
4 Dynamics of car-following models 33
x
4.1 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.1.1 Stability analysis of car-following model . . . . . . . . . . . . 34
4.2 Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2.1 Flux-density relation . . . . . . . . . . . . . . . . . . . . . . . 39
4.2.2 Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.3 Collision behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3.1 Simulation of two-vehicle system . . . . . . . . . . . . . . . . 45
4.3.2 Proof of collision-free behavior . . . . . . . . . . . . . . . . . . 49
4.3.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5 Data assimilation for traffic flow estimation 69
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2 Underlying macroscopic traffic model . . . . . . . . . . . . . . . . . . 72
5.3 Eulerian and Lagrangian observations . . . . . . . . . . . . . . . . . . 73
5.3.1 Continuous observation models . . . . . . . . . . . . . . . . . 74
5.3.2 Discretized observation models . . . . . . . . . . . . . . . . . 77
5.4 Ensemble methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.4.1 Particle filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.4.2 Ensemble Kalman filter . . . . . . . . . . . . . . . . . . . . . . 82
5.4.3 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . 86
5.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.5.1 Performance criteria . . . . . . . . . . . . . . . . . . . . . . . 87
5.5.2 Traffic scenarios with true underlying model . . . . . . . . . . 88
5.5.3 Traffic scenarios with unknown underlying model . . . . . . . 99
5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6 Traffic control in lane-drop scenarios 110
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.2 Inhomogeneous LWR model . . . . . . . . . . . . . . . . . . . . . . . 112
6.2.1 Model description . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.2.2 Solutions to Riemann problems . . . . . . . . . . . . . . . . . 117
6.2.3 Traffic pattern on a ring road . . . . . . . . . . . . . . . . . . 124
6.3 Controlling strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
6.3.1 Virtual lane usage . . . . . . . . . . . . . . . . . . . . . . . . . 140
6.3.2 Reshaped fundamental diagram . . . . . . . . . . . . . . . . . 143
6.3.3 Optimization in controlling strategies . . . . . . . . . . . . . . 145
6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
7 Conclusion 153
Bibliography 158
xi
List of Tables
4.1 Parameters from auto-calibration of the fundamental diagram . . . . 42
4.2 Parameter values for car-following models. Type (a) parameters de-
termine the flux-density relation, while type (b) parameters describe
the driving habit and sensitivity terms. . . . . . . . . . . . . . . . . . 44
4.3 Setting of time steps in Monte Carlo simulation . . . . . . . . . . . . 64
6.1 Parameters from auto-calibration of the fundamental diagram . . . . 116
xii
List of Figures
2.1 Fundamental diagram from [64]. Flow-density data obtained from
sensor measurement data, aggregated over time intervals of ∆t =
30s (left). The free flow curve and a synchronized flow region of
fundamental diagram (right). . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Velocity functions and respective fundamental diagrams for Green-
shields (left) and Daganzo-Newell (right). . . . . . . . . . . . . . . . . 13
2.3 Characteristic curves of the Riemann problem when wave speeds (a)
λL = λR , (b) λL > λR and (c) λL < λR . . . . . . . . . . . . . . . . . 15
2.4 On-/off-ramps allow vehicles to enter or exit a controlled-access highway. 19
2.5 Deceleration factor a(x, t; x`1 ) in traffic light scenario: (a) yellow light
period; (b) Red light period; (c) Green light period. . . . . . . . . . . 21
3.1 Notations for car-following models: vehicles are numbered from front to
back such that vehicle n follows vehicle n − 1. For vehicle n, xn denotes
position, vn denotes velocity, `n denotes length and dn denotes gap. . . . . 27
3.2 Optimal velocity function in respect to parameter V0 , ∆s and β . . . 29
4.1 Notations for car-following models: vehicles are numbered from front to
back such that vehicle n follows vehicle n − 1. For vehicle n, xn denotes
position, vn denotes velocity, `n denotes length and dn denotes gap. . . . . 34
4.2 Automatic calibration of the fundamental diagram from [24]. (a) The
scatter of flux-density data collected by PeMS. (b) The estimation of
main parameters by regression. . . . . . . . . . . . . . . . . . . . . . 42
4.3 Calibration for flux-density relation. (a) The flux-density relation
for (a1) in equation (4.9). (b) The flux-density relation for (a2) in
equation (4.9). The blue solid line represents the flux-density relation
for the car-following models, while the red dashed line represents the
relation from PeMS traffic data. . . . . . . . . . . . . . . . . . . . . . 43
4.4 Phase portraits of car-following models. The x-axis represents the gap
d2 between vehicles, and the y-axis represents the speed difference w2 .
The red dashed line d2 = 0 represents the critical situation of collisions. 47
4.5 Analysis of collision-prone behavior in optimal velocity model and full ve-
locity difference model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.6 Slope comparison. System 1 (dashed green line) represents the fast system
(4.15) and system 2 (solid blue) represents the original car-following system
(4.14). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
xiii
4.7 Trajectories of car-following models. The x-axis represents the bins of
collisions observed, and the y-axis represents the counts of collision in
percentage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.8 The relationship of collision counts over time step h. The red dots
represent the simulation results, and the blue line is a fitted line from
linear regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.9 Region partition of initial condition (d, w): the red checks and dia-
monds area is collision-prone; the blue lines area is collision-free. . . . 67
4.10 (d, w) collision region for different time step h (h1 : h2 : · · · : h8 : h9 =
1 : 2 : · · · : 8 : 9). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1 Evolution into a shock in the viscous Lighthill–Whitham equation (5.1).
The shock is moving backward with gradually decreasing amplitude. . 72
5.2 Comparison of the scaled Kalman gain matrices with localization applied
to the forecast covariance matrix (left) and directly to the Kalman gain
matrix (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.3 Traffic state estimation in the traffic light scenario. Shown is the rela-
tive RMSE for assimilating no data (yellow dashed dot), Eulerian data
(blue dashed), Lagrangian data (magenta dotted), and both Eulerian
and Lagrangian data (green solid) using the (a) EnKF and (b) PF. . 91
5.4 Traffic state estimation comparison for Lagrangian data assimilation
based on different observation data from normal traffic flow modeled
by (2.10). Shown is the relative RMSE for the (a) EnKF and (b)
PF, where we assimilate position data only (solid), velocity data only
(dotted), and combined position and velocity data (dashed). . . . . . 93
5.5 Evaluation of various sensors configurations using the EnKF. Results
for different sensor locations near (a) on-ramps, (b) off-ramps, and (c)
bottlenecks are shown, where the solid curves correspond to sensor lo-
cations upstream to the target location, the dotted curves represent
sensors at the target location, and the dashed curves represent loca-
tions downstream to the target location. . . . . . . . . . . . . . . . . 95
5.6 Estimation of traffic states and parameters in the ramps scenario:
shown are the relative RMSE for assimilating (a) Eulerian and (b)
Lagrangian observations as well as the results of parameter estimation
for (c) Eulerian and (d) Lagrangian observations. Estimated are the
maximal density ρmax , the maximal velocity vmax , the flux of the on-
ramp ρon off
1 , and the flux of the off-ramp ρ1 . PEOFF represents that
the parameters are known, while PEON represents that the parameter
estimation is considered. . . . . . . . . . . . . . . . . . . . . . . . . 97
5.7 Estimation of traffic states and parameters in the traveling bottleneck
traffic scenario: shown are the relative RMSE for assimilating (a)
Eulerian and (b) Lagrangian observations and the results of parameter
estimation for (c) Eulerian and (d) Lagrangian observations, where we
estimate the maximal density ρmax , the maximal velocity vmax , and the
time-dependent bottleneck location xb1 . . . . . . . . . . . . . . . . . . 98
5.8 Estimating microscopic data in scenario 1 using the EnKF. True and
estimated density values are represented by blue pluses and green
circles, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.9 Estimating microscopic data in scenario 2 using the EnKF. True and
estimated density values are represented by blue pluses and green
circles, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
xiv
5.10 A strip of highway I-35E in Minnesota. Cars are moving from left to right.
Sensors, on-ramps, and off-ramps are labelled using blue stars, red pluses,
and green circles, respectively. . . . . . . . . . . . . . . . . . . . . . . . 105
5.11 Relative RSME (a) and parameter estimate (b) are shown for data
assimilation of real traffic data taken during late night. Observations
are taken from all sensors (solid), half the sensors (dotted), or no
sensors (dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.12 Relative RSME (a) and parameter estimate (b) are shown for data
assimilation of real traffic data collected during rush hour. Observa-
tions are taken from all sensors (solid), half the sensors (dotted), or
no sensors (dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.1 Lane-drop scenario from two lanes to one lane. . . . . . . . . . . . . . . 115
6.2 Fundamental diagram ϕ(ρ). The slope of the free flow region is the max-
imal velocity vmax , and the slope of the congested flow is the congestion
parameter ω. The maximal density is ρmax , and the maximal flux is ϕmax . 116
6.3 The flux-density relation of ϕtotal (ρtotal , I(x)) in lane-drop scenario.
The red solid line represents the fundamental diagram for road of one
lane, while the green dashed dotted line is for road of two lanes. . . . 117
6.4 The wave solution to Riemann problem of inhomogeneous LWR. ρL
and ρR are the initial conditions. ρA and ρB are the transitional traffic
states in the upstream and downstream of the standing wave. . . . . 122
6.5 The weak solution to Riemann problem (6.15) of inhomogeneous LWR
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.6 A ring road with two branches: one-lane branch with length L1 and
two-lane branch with length L2 . The ratio is L2 : L1 = 1 : 1. . . . . . 124
6.7 Initial conditions of the Riemann problem on a ring road with two
branches. The first branch is of two lanes from −L2 to 0, and the
second branch is of one lane from 0 to L1 . The boundary is periodic. 125
6.8 The fundamental diagram of the ring road. There are two points C1
and C2 which can provide the road capacity ϕmax . . . . . . . . . . . . 126
6.9 The weak solution to Riemann problem (6.20) on a ring road. The
plots are vehicle density, velocity and flux from top to bottom. . . . . 127
6.10 The weak solution to Riemann problem (6.21) on a ring road. The
plots are vehicle density, velocity and flux curves from top to bottom. 128
6.11 The weak solution to Riemann problem (6.22) on a ring road. The
plots are vehicle density, velocity and flux curves from top to bottom. 130
6.12 The weak solution to Riemann problem ( ρ1 +ρ 2
2
= 4.615ρC1 ) on a ring
road. The plots are vehicle density, velocity and flux curves from top
to bottom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
6.13 The fundamental diagram of the ring road. The density ρ˜2 has equal
vehicle flux to the given initial density ρ1 . The density ρ˜1 has equal
vehicle flux to the given initial density ρ2 . . . . . . . . . . . . . . . . 132
6.14 The weak solution to Riemann problem (6.22) on a ring road. The
plots are vehicle density, velocity and flux curves from top to bottom. 133
6.15 The weak solution to Riemann problem (6.24) on a ring road. The
plots are vehicle density, velocity and flux curves from top to bottom. 135
6.16 A ring road with two branches: k1 -lane branch with length Lk1 and
k2 -lane branch with length Lk2 . The ratio is Lk2 : Lk1 = λ2 : λ1 . . . . 135
xv
6.17 The evolution of the solution to Riemann problem (6.29). In the
top figures, a traveling wave moving forward with speed vmax . In the
bottom figures, the solution consists of a shock wave moving backward,
a standing wave at the origin and a traveling wave moving forward. . 138
6.18 The wave solution to (6.29). The plots are density curve, velocity
curve and flux curve from top to bottom. . . . . . . . . . . . . . . . . 139
6.19 Comparison of fundamental diagram: (a) lane number I(x) decreases
from 2 to 1; (b) lane number I(x) decreases from 2 to 1.5, then to 1. . 141
6.20 The wave solution to (6.29) after introduce a virtual lane I(x) = 1.5.
The plots are density curve, velocity curve and flux curve from top to
bottom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
6.21 Comparison of fundamental diagram: (a) original shape with transi-
tional state A; (b) alternative shape with transitional state A0 . . . . . 144
6.22 The wave solution to (6.29) after reshaping the fundamental diagram.
The plots are density curve, velocity curve and flux curve from top to
bottom. The original wave solution is in blue dashed line. . . . . . . . 144
6.23 Lane number I(x) (top) and average velocity vave (x) (bottom) using
virtual lane usage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
6.24 The relation between lane number I(x) and velocity vave (x). (a)&(b)
The slopes sLA and sAL corresponding to lane number I(x); (c) fun-
damental diagram with virtual lane. . . . . . . . . . . . . . . . . . . 148
6.25 The velocity gap and shock wave speed in reshaped fundamental di-
agram. A0 is the intersection of the alternative fundamental diagram
with line BA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
xvi
Chapter One
Introduction
2
1.1 Motivation and preliminaries
Traffic flow theory and modeling started in the 1930. This field has gained consider-
able attention as overall traffic demand has increased and more data as well as easy
access to computing power has become available. Traffic is everywhere in our daily
life. There is a growing need for accurate traffic information so that the public can
efficiently schedule their trips, and the government can better provide traffic control
strategies. Thus, modeling, analyzing, estimating and controlling the dynamics of
traffic flow are of great importance. First of all, I would like to have a overview of
traffic flow: (a) traffic-flow models (b) traffic-flow observation.
Mathematical models for traffic flow come in many different flavors. Common
models range from microscopic models for individual cars to macroscopic models
for car densities and possibly other quantities. As the name indicates, macroscopic
models formulate the relationship among the spatial quantities such as density, ve-
locity and flow. In contrast, microscopic models simulate individual vehicle-driver
unites, so the dynamic variables of the models represent properties like the position
and velocity of individual vehicles.
The macroscopic traffic models are usually described in a partial-differential equa-
tion (PDE) to build a relationship between spatial quantities [39, 54, 63]. The
most well-known macroscopic model is the Lighthill-Whitham-Richard (LWR) model
[54, 63], which describes the the evolution of traffic by the vehicle density accord-
ing the conservation law. The LWR model can be modified to account for multiple
traffic scenarios through velocity-density relation or flux-density relation. The com-
mon traffic scenarios include normal traffic, traffic lights, on-/off-ramps, construction
zones, traveling bottlenecks and lane drops.
3
The microscopic traffic models include the acceleration and lane-changing mod-
els [61]. The acceleration models are formulated by coupled ordinary differential
equations, which describe the variables of individual vehicles, such as positions and
velocities. The lane-changing models involve the vehicles interaction across different
lanes, and specify the rules for merging. The microscopic models are able to provide
traffic simulation in a realistic road profile, and the well-known traffic simulators in-
clude AIMSUN (advanced interactive microscopic simulator for urban and nonurban
networks), PARAMICS (parallel microscopic simulation) and VISSIM (Verkehr in
Stadten simulation) [59]. The dynamics of the microscopic models is an interesting
area so that the models can better simulate the realistic traffic.
Different aspects of traffic dynamics are captured by different measurement meth-
ods. Traffic-flow observations come from a variety of sources and are available as
functions of time. Stationary sensors, such as induction loops or cameras, provide
the flux, average velocity, and local density of cars that move pass the fixed sensor
location. GPS data from cell phone or navigation devices, on the other hand, provide
information about the positions and velocities of individual cars that move with the
traffic flow. We refer to observations that come from a fixed observation location as
Eulerian observations and to observations that come from parcels that move with the
traffic flow as Lagrangian observations. Both Eulerian and Lagrangian observations
have their own advantages and disadvantages with respect to collection cost, data
accuracy, traffic coverage or even driver privacy. A lot of comparative studies are
done to compare the properties of both types of observation. They are used in dif-
ferent circumstances for calibration, estimation and prediction based on the impacts
of their advantages and disadvantages.
4
1.2 Thesis objectives
There are a lot of interesting traffic topics, ranging from the theoretic study of
traffic-flow models, to their application in realistic traffic scenarios. There are three
objectives related to traffic-flow modeling in my thesis:
• Studying the dynamics of car-following models
Car-following models are widely used in traffic simulation because of its vari-
ability in the formula of ordinary differential equations. The traffic dynamics
of these models are of great importance, and thus have received plenty of study
from scholars. However, the collision behavior of car-following models remains
a challenging problem. One difficulty is that the traffic system of multiple
vehicles is complicated when all interactions among vehicles are taken into
consideration. Previous study for making inferences about the collision behav-
ior depended on vast numerical simulations. On one side, the study cannot
provide a theoretical analysis of the car-following models from its essence. On
the other side, the inferences from the simulations could be misleading in the
existence of numerical error. In my thesis, we carry out a theoretical analysis of
the collision behavior of four well-known car-following models. We show that
the modified GHR model and the intelligent driver model are collision-free un-
der all traffic circumstances. We also show that the numerical errors introduce
collisions that the model doesn’t support.
• Assimilating Eulerian and Lagrangian data in traffic-flow models
Data assimilation of traffic flow remains a challenging problem. One diffi-
culty is that data come from different sources ranging from stationary sensors
and camera data to GPS and cell phone data from moving cars. Sensors and
5
cameras give information about traffic density, while GPS data provide infor-
mation about the positions and velocities of individual cars. Previous methods
for assimilating Lagrangian data collected from individual cars relied on spe-
cific properties of the underlying computational model or its reformulation in
Lagrangian coordinates. These approaches make it hard to assimilate both Eu-
lerian density and Lagrangian positional data simultaneously. In this thesis, we
propose an alternative approach that allows us to assimilate both Eulerian and
Lagrangian data. We show that the proposed algorithm is accurate and works
well in different traffic scenarios and regardless of whether ensemble Kalman
or particle filters are used. We also show that the algorithm is capable of
estimating parameters and assimilating real traffic observations and synthetic
observations obtained from microscopic models.
• Providing traffic control strategies to maximize capacity in lane-drop bottle-
neck
Lane drop is a location where number of lanes provided decreases. Empiri-
cal observations at lane-drop bottleneck revealed that there is a drop in the
bottleneck discharge rate when queues form in the upstream of the lane drop.
Even though the research results show that lane changes are the main cause of
the drop in discharge rate, the relevant traffic control to maximize the capacity
is still a challenging problem. One difficulty is that the interaction between
lanes are complicated to quantify, and that there are a lot of parameters to
specify. Previous methods to maximize capacity either attempted to develop
cooperative lane-changing models to reduce the queues in the upstream of lane
drop or formulated a conservation equation for individual lane to optimize the
density transfer between adjacent lanes. These approaches focus on the inter-
action rules between lanes, which are complicated to quantify, and not target
to maximize capacity directly. In my thesis, we introduce the inhomogeneous
6
LWR model to account for the change of lane number, and use one continuity
equation to model the total density. When the maximal capacity is guaranteed
by the inhomogeneous LWR model, we attempt to smooth the velocity curve
so that the embedded vehicles are comfortable to follow it. We show the traf-
fic pattern by the inhomogeneous LWR model and its property to maximize
capacity. We also show two controlling strategies to smooth and narrow the
velocity gap in the solution to the inhomogeneous LWR model.
Outlines
The work is organized as follows. In Chapter 2, we provide an overview of macro-
scopic traffic models (LWR model) as well as the fundamental diagram, and a discus-
sion of multiple traffic scenarios based on LWR model, which will be used to test the
efficacy of our proposed assimilation method in Chapter 5. In Chapter 3, we provide
an overview of microscopic traffic models, and give a detailed introduction to the
well-known car-following models, along with notations and setup for the following
chapter. Chapter 4 provides the analysis of the stability and collision behavior of the
car-following models. The theoretical proof as well as the numerical implementation
is shown to be in agreement. In Chapter 5, we have a detailed discussion of the
traffic observations and data assimilation approaches. Then we proposed an alterna-
tive approach that allows us to assimilate both Eulerian and Lagrangian data, and
show the accuracy and efficacy of this approach in multiple traffic scenarios. The
lane-drop bottleneck is investigated in Chapter 6. The inhomogeneous LWR model
is introduced to account for the multi-lane traffic scenarios. We provide two control-
ling strategies to maximize the capacity while providing a comfortable velocity curve
for vehicles to follow. Finally, Chapter 7 includes an overall discussion, conclusion
and directions for future work.
Chapter Two
Review: macroscopic traffic
models
8
2.1 Introduction
A macroscopic traffic model is a mathematical traffic model that formulates the
relationships among traffic flow characteristics like density, flow, mean speed of a
traffic stream. The characteristics are locally aggregated, which vary across space
and time, i.e., they correspond to dynamic fields. Thus, macroscopic models are
able to describe collective phenomena such as the evolution of congested regions or
the propagation of traffic wave ([73]). Macroscopic models [19, 42, 54, 63] describe
traffic flow analogously to liquids or gases in motion. Hence they are sometimes
called hydrodynamics models.
Since the legendary paper [54, 63] by Lighthill, Whitham and Richards, dynam-
ical macroscopic traffic flow modeling became a central focus for both theoretical-
and application-oriented research. This is a first order model based on the scalar
conservation law. The first-order model was complemented by second-order model
by Payne [60], in an attempt to avoid some known deficiencies of the first-order
model. However, in Daganzo’s note [20], it described the logical flows of the higher
order continuum models. In our research, we will focus the first-order model [54, 63]
and its variants.
Conventionally, the macroscopic models are derived from integrating the micro-
scopic traffic flow models and converting the single-entity level characteristics to
system level characteristics ([5, 25]). In [5], it establishes a connection between a
microscopic follow-the-leader model based on ordinary differential equation and a
semidiscretization of a macroscopic continuum model based on a conservation law.
They also show rigorously that, at least in the homogeneous case, the macroscopic
model can be viewed as the limit of the time discretization of the microscopic model
9
as the number of vehicles increase, with a scaling in space and time.
The foundations of macroscopic model are the hydrodynamic relation “flow equals
density times speed” and the continuity equation, which describes the temporal
evolution of the traffic flow characteristics. The vehicle flux-density relation is called
the fundamental diagram of traffic, which is a modeling choice. Traffic flux refers
to the number of passing vehicles per unit time, and traffic density refers to the
number of vehicles per unit length. Then the vehicle velocity is defined such as the
flow-density relation is satisfied. More details about the fundamental diagram will
be provided in section §2.2.1.
The macroscopic models are particularly useful when one is interested in macro-
scopic quantities and the microscopic effects (lane changes, driver-vehicle type, ac-
celeration) need not be considered. Some main applications of macroscopic model
include:
• Macroscopic models can be used to study the spatiotemporal evolution of con-
gested traffic pattern, and simulate the effects of traffic flow breakdown caused
by high traffic load, bottleneck and disturbance of drivers.
• Macroscopic models can be modified to account for realistic road profiles, such
as traffic lights, on-/off-ramps, change in the number of lanes, intersections,
construction zone and so on.
• Flow and aggregated speed data from stationary detectors and trajectories of
moving cars allow us to estimate and predict the traffic states. Since most
of this information is not well known, it is typically more common to use
macroscopic models for traffic estimation.
• From the analysis of the spatiotemporal dynamics on highways, we can im-
plement model-based traffic flow optimization to increase the efficiency and
10
stability of traffic flow.
The traffic pattern and wave propagation are well studied in [49, 54, 63, 73]. In this
chapter, we touch on the variations of macroscopic model to account for road pro-
files mentioned above. In chapter §5, we will apply the Lighthill-Whitham-Richards
model to provide traffic evolution at later times, and focus on the estimation of traffic
states and uncertain parameters. In chapter §6, we study the lane-drop bottleneck
and discuss the potential optimization strategies based on the Lighthill-Whitham-
Richards model.
2.2 Lighthill-Whitham-Richard model
This section reviews the theory of scalar first-order conservation law, known as the
Lighthill-Whitham-Richard (LWR) partial differential model [54, 63]. The LWR
model describes the evolution of traffic by the vehicle density ρ(x, t) at location x
and time t. For simplicity, we consider a ring road of length L but emphasize that
roads with other boundary conditions can also be treated by our approach. In the
absence of sinks and sources, the conservation law for the vehicle density ρ(x, t) with
periodic boundary conditions is given by
∂t ρ(x, t) + ∂x ϕ(ρ(x, t)) = 0, x ∈ T, (2.1)
ρ(x, 0) = ρ0 (x),
where ρ0 (x) denotes the initial data, and T denotes the circle of circumference L
that reflects the periodic boundary conditions we employ. The equation describes
the rate of change of density in terms of gradients of the flow.
11
The flux function ϕ(ρ) expresses the dependence of the vehicle flux ϕ on the
density ρ; this relationship is usually referred to as the fundamental diagram. We
assume that ϕ is defined on the interval [0, ρmax ], where ρmax is the maximal density
achievable on the road.
For traffic flow, we can write ϕ(ρ) = ρV (ρ), where the function V (ρ) relates
velocity v and density ρ for densities ρ ∈ [0, ρmax ]. The function V (ρ) is a mod-
eling choice: examples include the Daganzo–Newell velocity function [76] and the
Greenshields affine velocity function [36].
2.2.1 Fundamental diagram
Fundamental diagram of traffic flow is the plot of the traffic flux ϕ versus the traffic
density ρ. According to [64], the two-phase traffic theory divides traffic flow into
free flow for low densities, and congested flow for large densities. Figure 2.1 is
the fundamental diagram from Seibold’s paper [64]. In the free flow regime, flux
increases as density increases, while in the congested flow regime, flux decrease as
density increases.
The flux-density relation can be mathematically written as ϕ(ρ) = ρV (ρ), where
velocity function V (ρ) relates velocity v and density ρ. The algebraic expression of
the velocity function is a modeling choice, and it is typically constructed to fit ex-
perimental data. We will briefly introduce the the Daganzo–Newell velocity function
[76] and the Greenshields affine velocity function [36], which are used in §5 and §6.
12
Figure 2.1: Fundamental diagram from [64]. Flow-density data obtained from sensor measurement
data, aggregated over time intervals of ∆t = 30s (left). The free flow curve and a synchronized
flow region of fundamental diagram (right).
Greenshields velcoity function: It is one of the earliest introduced velocity
functions by Greenshields. It expresses a linear relationship between the speed and
density: the velocity linearly decreases as the density increases.
ρ
v = VG (ρ) = vmax 1 − (2.2)
ρmax
where vmax is the maximal free-flow velocity, and ρmax is the maximal density achiev-
able on the road. This model remains useful because of its simplicity. The velocity
function and respective fundamental diagram are shown in Figure 2.2.
Daganzo-Newell velocity function: It is a widely used velocity function by
assuming a constant velocity in free-flow regime, and a hyperbolic velocity in the
congestion regime:
v
max , if ρ ≤ ρc
v = VDN = (2.3)
ρmax
−ω 1 −
, if ρ > ρc
ρ
13
where ρc is called the critical density that separates the free-flow regime and conges-
tion regime. The Daganzo-Newell velocity gives a triangular fundamental diagram.
The velocity function and respective fundamental diagram are shown in Figure 2.2.
V ( ⇢)
V ( ⇢)
⇢ ⇢
( ⇢)
( ⇢)
⇢ ⇢
Figure 2.2: Velocity functions and respective fundamental diagrams for Greenshields (left) and
Daganzo-Newell (right).
2.2.2 Characteristic curve
The LWR equation is a first-order PDE, we could use the method of characteristics to
discover curves along which the PDE becomes an ODE. The LWR equation without
viscous term can be written in a quasilinear form:
∂ρ ∂ρ
+ ϕ0 (ρ) =0 (2.4)
∂t ∂x
14
The characteristic equations for the original system with initial conditions are:
∂t
∂s
= 1, t(0) = 0
∂x
∂s
= ϕ0 (ρ), x(0) = x0
∂ρ = 0, ρ(0) = ρ(x0 , 0)
∂s
Because of t = s from the first equation, x and ρ can be expressed as functions of t:
x(t) = ϕ0 (ρ)t + x0
ρ(t) = ρ(x0 , 0) = ρ0 (x0 )
Therefore, we get this relation ρ(x(t), t) = ρ0 (x0 ) = ρ0 (x − ϕ0 (ρ)t). In this case, the
characteristic lines are straight lines with slope ϕ0 (ρ), and the value of ρ remains
constant along the characteristic line. Therefore, the solution of Lighthill–Whitham
model comes from following characteristic curve back from (x, t) to a point when
t = 0:
ρ(x, t) = ρ0 (x − ϕ0 (ρ(x, t))t) (2.5)
Riemann problem: The initial problem (2.4) with discontinuous initial condition
of the form
ρ L ,
x≤0
ρ(x, 0) = (2.6)
ρ R ,
x>0
is called a Riemann problem. By using the characteristic method presented above,
we can obtain two characteristic curves:
x = ϕ0 (ρL )t + x0
x = ϕ0 (ρR )t + x0 .
15
The slopes of the characteristic curve can be understood as wave speed. We use new
notations λL and λR to represent the wave speeds in Riemann problem.
Solutions Based on the relation of λL and λR , the solution to the Riemann problem
can be traveling wave, shock wave and rarefaction. The relations include λL = λR ,
λL > λR and λL < λR , which are plotted in Figure 2.3.
t t t
0 x 0 x 0 x
(a) λL = λR (b) λL > λR (c) λL < λR
Figure 2.3: Characteristic curves of the Riemann problem when wave speeds (a) λL = λR , (b)
λL > λR and (c) λL < λR
We substitute the Greenshields velocity function into the flux function, and
obtain ϕ(ρ) = ρvmax (1 − ρ/ρmax ). The slope of the characteristic curve becomes
ϕ0 (ρ) = vmax (1 − 2ρ/ρmax ), which is a linear decreasing function respect to density
ρ. The Riemann solutions include:
• Traveling wave (λL = λR = λ)
The wave doesn’t change the speed from upstream to downstream, so it is a
traveling wave with speed λ. The solution can be expressed as
ρ(x, t) = ρ0 (x − λt) (2.7)
• Shock wave (λL > λR )
When the wave speed in upstream is greater than downstream, a discontinuity
16
will be generated and propagated. This is shock wave.
The Rankine-Hugoniot jump condition determines the position of a shock at a
given time. This condition emerges when one consider the equation in integral
form by integrating the LWR respect to x:
Z x2
d
ρ(x, t)dx + ϕ(ρ)|xx21 = 0
dt x1
Then we can obtain the wave speed of the shock
ϕ(ρL ) − ϕ(ρR ) λL + λR
λ= = (2.8)
ρL − ρR 2
The discontinuity propagating with speed λ satisfies the entropy condition
because λL > λ > λR holds.
• Rarefaction (λL < λR )
The whole situation changes. Even though there is a shock wave solution
satisfying the Rankine-Hugoniot condition, the shock wave doesn’t satisfy the
entropy condition because of λL < λ < λR .
A weak solution satisfying the entropy condition is rarefaction wave:
ρL x < λL t
ρ(x, t) = w(x/t) λL t ≤ x ≤ λR t , (2.9)
ρ R
x > λR t
where w(.) is a smooth function with w(λL ) = ρL and w(λR ) = ρR .
So far, we are discussing the homogeneous LWR model and the solutions to the
Riemann problem. In Chapter 6, we will introduce the inhomogeneous LWR model
17
to account for the change of lane number.
2.3 Traffic scenarios based on LWR model
As we know, the original LWR equation (2.1) models the conservation of traffic flow
in absence of sources or sinks. In addition, the velocity function V (ρ) is assumed
not dynamic, which is unable to describes a dynamic relation between velocity and
density. Therefore, the original LWR equation is limited to model more complicated
road profiles.
The main road profiles include normal traffic (§ 2.3.1), on-/off-ramp (§ 2.3.2),
traffic lights (§ 2.3.3), construction zone (§ 2.3.4), traveling wave (§ 2.3.5), change in
the number of lanes (§ 6.2), and so on. In this section, we modify the existing LWR
model to account for effects of the first five traffic scenarios. The modifications make
the macroscopic model more accurate to describe the realistic traffic. The change in
the number of lanes is a particular traffic scenario we will discuss in more details in
§ 6.2.
2.3.1 Normal traffic
In the LWR model (2.1), the adaption of speed to traffic density is instantaneous,
because the relation v(x, t) = ϕ(ρ(x, t))/ρ(x, t) holds if function ϕ and ρ are appro-
priately defined. This means that it takes zero time for a driver to change driving
speed according to the density. It is normal that some fine structure of traffic is
missing in the LWR model.
18
In the normal traffic, the drivers need response time to adjust the velocity, so
the relation v(x, t) = ϕ(ρ(x, t))/ρ(x, t) doesn’t hold when density experiences big
changes. In [81], the researcher exploit the relation between driver memory and
viscosity to develop a viscous continuum model, which takes the drivers’ response
time into consideration through a diffusion term on the right-hand-side of the LWR
equation. Therefore, we introduce the viscous LWR model to describe the normal
traffic:
ρt + (ρVG (ρ))x = ερxx , (2.10)
where ε is the diffusion coefficient, which acts as a smoothing factor in the model.
The diffusion term with diffusion coefficient ε has two other advantages: this
term accounts for low-level noise when updating the traffic states at later times;
from the perspective of numerical solution, the introduction of small diffusion term
can produce a weak solution that satisfies the entropy condition.
The viscous LWR model (2.10) will act as the base model for the following traffic
scenarios: on-/off-ramps, traffic lights, construction zone and traveling bottleneck.
2.3.2 On-/off-ramps
Ramp is a road junction of the highway system: an on-ramp provides access to the
specific part of a road system, while an off-ramp is one-way lane for departing a
main highway. The sketch of on-/off-ramps is shown in Figure 2.4. Mathematically,
on-/off-ramps imply additional in-/outflows, which have to be added to the section
boundaries where the on-/off-ramps are.
19
xo↵
j xon
j
Figure 2.4: On-/off-ramps allow vehicles to enter or exit a controlled-access highway.
Let {xon off
i }i∈I and {xj }j∈J be the positions of on-ramps and off-ramps, and
{ϕon off
i (t)}i∈I and {ϕj (t)}j∈J are time-dependent flows of on-ramps and off-ramps
respectively. The ramp flow is positive for on-ramps, and negative for off-ramps.
Instead of using a delta function δ(x − xon off
i ) and δ(x − xi ) to indicate the loca-
tions of sources sinks, we introduce sech functions f (x − xon on
i ) = sech(x − xi ) and
f (x−xoff off
i ) = sech(x−xi ) to represent the spread effect of ramps. This consideration
comes from the phenomenon that ramps would affect the areas nearby rather the
specific points.
By addition source and sink terms to the viscous LWR equation (2.10), we get
the model for on-/off-ramps scenario:
X X
ρt + (ρVG (ρ))x = ερxx + ϕon on
i (t)f (x − xi ) + ϕoff off
j (t)f (x − xj ), (2.11)
i∈I j∈J
In this scenario, spatial shocks are generated around the on-ramp and off-ramp.
Specifically, a density bump appears around the on-ramp, while a density valley
appears around the off-ramp. More details about the ramp positions and ramp flows
will be mentioned in § 5 when estimating traffic states in on-/off-ramps scenario.
20
2.3.3 Traffic lights
The traffic lights on the road force vehicles to decelerate and accelerate periodi-
cally, resulting a periodic oscillations of the vehicle density along the road. In this
circumstance the velocity function is dynamic corresponding to the color of traffic
lights.
In the traffic light scenario, we place a traffic light at position x`1 on the road,
and set three sub-periods for yellow light T1y , red light Tir and green light T1g . We
introduce a deceleration factor to account for the effects from traffic light, and the
velocity function becomes:
v(x, t) = VG (ρ(x, t))a(x, t; x`1 ) (2.12)
where a(x, t; x`1 ) is the deceleration factor from a traffic light at position x`1 . The
function a(x, t; x`1 ) is periodic in time t with period T , which is the total length of
three sub-periods T = Tiy + Tir + Tig . In Figure 2.5, The traffic light cycle has a
period of T = 6 minutes, which consists of yellow light time T1y , red light time T1r
and green light time T1g as follows:
• During yellow lights, drivers within Ly1 miles (see Figure 2.5(a)) notice the
yellow lights, and they are decelerating. The drivers out of the range will drive
normally.
• During red lights, cars within Lr1 miles (see Figure 2.5(b)) are forced to stop,
and those following within 2Lr1 are decelerating when approaching the traffic
light.
• During green lights T1g (see Figure 2.5(c)), traffic flows normally and the ve-
21
locity is exactly equal to the predetermined velocity function.
(a) a(x, t; x`1 ), t 2 T1y (b) a(x, t; x`1 ), t 2 T1r (c) a(x, t; x`1 ), t 2 T1g
1 1 1
0 0 0
0 x`1 Ly1 x`1 0 x`1 2Lr1 x`1 Lr1 x`1 0 x`1
(a) Yellow light period T1y (b) Red light period T1r (c) Green light period T1g
Figure 2.5: Deceleration factor a(x, t; x`1 ) in traffic light scenario: (a) yellow light period; (b) Red
light period; (c) Green light period.
The modified LWR model for traffic lights is mathematically given by:
ρt + ρVG (ρ(x, t)) a(x, t; x`1 ) x = ερxx ,
(2.13)
(x, t) ∈ (x`1 − Ly1 , x`1 ) × T1y
0.5
(x, t) ∈ (x`1 − Lr1 , x`1 ) × T1r
0
`
a(x, t; x1 ) = (2.14)
(x`1 − Lr1 − x)/Lr1 (x, t) ∈ (x`1 − 2Lr1 , x`1 − Lr1 ) × T1r
1
otherwise,
More details about the length of lights (T1y , T1r , T1g ) and the effect region (Ly1 , Lr1 ) will
be mentioned in § 5 when estimating traffic states in traffic light scenario.
2.3.4 Construction zone
When approaching a construction zone, we can see signs or signals like “SLOW”,
“ROAD WORK AHEAD”, “SINGLE LANE AHEAD” or even “STOP”. Around the
construction zone, typically the drivers are experiences slow speed and lane merge.
In this circumstance, the velocity function is dynamic corresponding to effects of
construction zone.
Let xb1 be the location of a construction zone, which is a constant. We introduce
22
a bottleneck factor a(x − xb1 ) to account for the slow-down effects from construction
zone. Then the velocity function becomes:
ρt + ρVG (ρ)a(x − xb1 ) x = ερxx ,
(2.15)
where xb1 is the location of the bottleneck. The quantity a(x − xb1 ) is the bottleneck
factor, which can be written as the product of severity coefficient and spread effect
a(x − xb1 ) = cf (x − xb1 ). The severity coefficient c reflects the degree of influence of
the construction zone or traffic accident, and f (x − xb1 ) represents the spread effect
of the bottleneck. In this scenario, a stationary density bottleneck appears around
xb1 .
More details about the bottleneck position xb1 , severity coefficient c and the spread
effect f (x−xb1 ) will be mentioned in § 5 when estimating traffic states in construction
zone.
2.3.5 Traveling bottleneck
We also construct a traveling bottleneck scenario, for which we use a model similar
to that of the stationary bottleneck in §2.3.4. To account for a traveling bottleneck,
we allow the position of the bottleneck to be time-dependent, and the macroscopic
model can be described as:
ρt + ρVG (ρ)a(x − xb1 (t)) x = ερxx ,
(2.16)
23
where xb1 (t) is the location of traveling bottleneck, which is assumed to be time-
dependent.
The traveling bottleneck scenario can be used to model slow truck, moving con-
struction zone, moving snowplows and so on. In this traffic scenario, we would like
test the assimilation efficacy of our algorithm in time-dependent parameter xb1 (t).
More details about the traveling bottleneck position xb1 (t) will be mentioned in
§ 5 when estimating traffic states and parameters.
2. Parameter estimation for real traffic data In addition to estimating sim-
ulated traffic states and parameters, we are interested in applying the developed
approach to real traffic data. In this section, we use data from the Minnesota De-
partment of Transportation [70].
Chapter Three
Review: microscopic traffic models
25
3.1 Introduction to microscopic traffic models
A microscopic traffic model is an import class of model in traffic simulation. In con-
trast to macroscopic traffic model, microscopic model simulates individual vehicle-
driver entities, so the dynamic variables of the models represent microscopic proper-
ties like the position and velocity of individual vehicles. Conventionally, the micro-
scopic models can be used to derive the macroscopic models through micro-macro
transition.
The mathematical formulations include car-following models and cellular automa-
ton models.
• The car-following models are time-continuous models, which are defined by
ordinary differential equations describing the vehicles’ positions and velocities.
It is assumed the accelerations of the drivers depend only on their own ve-
locities, the velocity of the leading vehicle, and the distance to the leading
vehicle. Usually the driving behavior can be affected by multiple leading ve-
hicles. Then the acceleration function can be generalized to account for more
variables. Well-known car-following models include the optimal velocity model
[8], the velocity difference model [46], the GHR model [32] and the intelligent
driver model [72].
• The cellular automaton models describe traffic dynamics in a completely dis-
crete way: space is subdivided into cells, time into time steps, and speed or
acceleration are integer multiples of the corresponding units. Each road section
can either be occupied by a vehicle or empty. Examples of cellular automaton
models include Nagel-Schreckenberg model[57], and other refined models.
26
Compared to the car-following models, the advantages and disadvantages of the
cellular automaton models are due to the discrete scaling of time, space and state
variables. Even though the cellular automaton models are easy and fast to simulate,
they lack robustness because of the discrete nature. Therefore, we will focus on
car-following models in our study.
Traffic simulation is an interesting and important application of car-following
models. Car-following models, together with lane-changing models [61] are devel-
oped and implemented in a microscopic simulation framework. In [59], a compara-
tive evaluation of car-following behavior in a number of well-known traffic simulators,
including AIMSUN (advanced interactive microscopic simulator for urban and nonur-
ban networks), PARAMICS (parallel microscopic simulation) and VISSIM (Verkehr
in Stadten simulation).
The dynamics of car-following models is another important area for investigation.
In Chapter 4, we will discuss the stability and collision behavior of car-following
models.
3.2 Car-following models
Microscopic models describe traffic flow dynamics in terms of single vehicles. As
the most important representatives of microscopic traffic flow models, car-following
models describe traffic dynamics from the perspective of individual driver-vehicle
units with interaction to the leading ones.
In continuous-time models, these car-following models are defined by coupled
ordinary differential equations, and the drivers’ responses are informed by their own
27
velocities, headways and the velocities of the leading vehicles. Many car-following
models are of the form:
vn+1 vn vn 1
vehicle n + 1 vehicle n vehicle n 1
xn+1 xn xn 1
`n dn
Figure 3.1: Notations for car-following models: vehicles are numbered from front to back
such that vehicle n follows vehicle n − 1. For vehicle n, xn denotes position, vn denotes
velocity, `n denotes length and dn denotes gap.
dxn (t)
x˙ n (t) = = vn (t),
dt (3.1)
dvn (t)
v˙ n (t) = = faccel (vn , vn−1 , dn ) , dn = xn−1 − xn − `n
dt
where
xn : the position of vehicle n,
vn : the speed of vehicle n,
vn−1 : the speed of the leading vehicle n − 1,
`n : the length of vehicle n,
dn : the gap between vehicle n and the leading vehicle n − 1,
faccel : the acceleration function.
In recent years, car-following models have received increasing attentions both in
mathematics and engineering. Many car-following models have been developed and
studied in the past half century. In [14, 17, 39], the dynamical behavior and proper-
ties of car-following models have been reviewed, including the formula, parameters,
steady state, stability, trajectories, etc.
28
In this section, we will focus on four well-known continuous car-following models,
which are frequently used in traffic simulations. These models are (i) the optimal ve-
locity model, (ii) the full velocity difference model, (iii) the modified Gazis-Herman-
Rothery model, and (iv) the intelligent driver model.
3.2.1 The optimal velocity model
The optimal velocity model was originally proposed by Bando et al. in [8], which
assumed that each vehicle has a distance-dependent optimal velocity V(·). In this
model, the driver controls the acceleration to maintain an optimal velocity according
to the motion of the leading vehicle. In addition, a constant adaption time τ is
introduced to reflect the driver’s sensitivity to adapt to the optimal velocity. The
optimal velocity model can be written as:
x˙ n (t) = vn (t)
(3.2)
V(dn ) − vn
v˙ n (t) = ,
τ
where dn is the gap to the leading vehicle, V(·) is the optimal velocity function
that depends on the distance to the leading vehicle, and τ is the adaption time.
The acceleration equation describes the adaption of actual speed vn to the optimal
velocity V(dn ) on a time scale given by the adaption time τ . The smaller the adaption
time τ is, the faster the vehicle adapts to the optimal velocity. The optimal velocity
function V(·) is a modeling choice: common choices are hyperbolic tangent and
piecewise function. In this study, we will use the hyperbolic tangent function by
Bando et al.:
tanh(dn /∆s − β) + tanh β
V(dn ) = V0
1 + tanh β
29
where V0 is the desired speed, ∆s is the transition width and β is the form factor.
Figure 3.2 shows how these parameters change the shape of optimal velocity function
V(·). The bigger the desired speed V0 is or the smaller the transition width ∆s is,
the greater the optimal velocity will be. For the form factor β, larger values shift
the optimal velocity curve to right. Overall the optimal velocity model is easy to
understand and implement, but it has a strong dependency on the parameters and
simulation results can be unrealistic.
80 80 80
Optimal Velocity (mph)
Optimal Velocity (mph)
Optimal Velocity (mph)
60 60 60
40 40 40
Desired speed V0 = 75mph Transition width s = 8m Form factor = 1.2
Desired speed V0 = 80mph Transition width s = 10m Form factor = 1.5
20 20 20 Form factor = 1.8
Desired speed V0 = 85mph Transition width s = 12m
0! 0! 0!
0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50
Headway (m) Headway (m) Headway (m)
(a) Desired speed V0 (b) Transition width ∆s (c) Form factor β
Figure 3.2: Optimal velocity function in respect to parameter V0 , ∆s and β
3.2.2 The velocity difference model
The main deficiency of optimal velocity model is that it does not take the speed of
the leading vehicle into consideration. Thus the full velocity mode was proposed in
[46] by extending the optimal velocity model with an additional linear stimulus for
the speed difference. On one hand, if the speed of vehicle n is larger than that of
the leading vehicle, additional deceleration will be generated to slow it down. On
the other hand, if the following vehicle is slower than the leading vehicle, additional
acceleration will be generated to speed it up. Thus, the full velocity difference model
has better stability than the optimal velocity model. The mathematical formula is
30
as follows:
x˙ n (t) = vn (t),
(3.3)
V(dn ) − vn
v˙ n (t) = − γ(vn − vn−1 )
τ
where (vn − vn−1 ) is the velocity difference to the leading vehicle n − 1, and γ is the
sensitivity factor for speed difference. Parameter γ measures the driver’s sensitivity
to speed difference: the greater γ is, the more sensitivity the driver has. When γ = 0,
the full velocity difference model (3.3) reduces to the optimal velocity model (3.2).
According to the simulations in [46], the full velocity difference model predicts the
correct delay time of car motion and kinematic wave speed at jam density. However,
the sensitivity term γ does not depend on the gap dn , which can generate a significant
and unrealistic deceleration even when the leading vehicle is far away.
3.2.3 The modified GHR model
The Gazis-Herman-Rothery (GHR) model is the most well-known model from the
late fifties by Gazis et al. ([32]). This model is based on the intuition that the
acceleration is proportional to the velocity difference; it also includes an explicit
time delay to account for driver’s reaction time. The general formula of the GHR
model is a delay ODE:
v˙ n (t + T ) = λ(vn−1 − vn )
where T is the reaction time, and λ is a relaxation rate. Here we modified the
original GHR model based on the following considerations: firstly we replace the
delay formula by the adaption to the optimal velocity with reaction time τ ; then we
31
η
use the reciprocal spacing for the relaxation rate, i.e., λ = ; lastly the sensitivity
dn
vn − vn−1
−η to velocity difference is bounded from above in order to avoid unrealistic
dn
acceleration generated by negative speed difference vn − vn−1 and small gap dn .
With these sensible modifications, the stimulus is able to characterize the distance-
dependent sensitivity to velocity difference as well as realistic acceleration. The
modified formula can be written as:
x˙ n (t) = vn (t),
V(dn ) − vn
vn − vn−1
(3.4)
v˙ n (t) = + min −η ,A ,
τ dn
where η is a constant sensitivity parameter, and A is the upper bound of acceleration.
The modified GHR model is an improvement over the optimal velocity model and
the full velocity difference model, since it generates distance-dependence stimulus to
velocity difference.
3.2.4 The intelligent driver model
The intelligent driver model is a novel continuous model by considering the real
driving behavior, which was proposed by Treiber et al. in the work [72]. In this
model, the driving behaviors such as keeping a ”safe distance”, driving at a desired
speed and preferring comfortable accelerations are taken into consideration. The
main idea behind this model is to superpose (I) the acceleration tendency on a
free road by comparing the current velocity to the desired velocity, and (II) the
deceleration tendency in the presence of interactions by comparing the current gap
32
to the desired safe distance. The mathematical formula is as follows:
x˙ n (t) = vn (t),
" δ # " ∗ 2 #
vn d (vn , vn − vn−1 ) (3.5)
v˙ n (t) = a 1 − + −a ,
V0 dn
| {z } | {z }
I II
where a is the comfortable acceleration, V0 is the desired velocity, δ is the accel-
eration exponent and d∗ is the desired distance. In model (3.5), term I represents
the acceleration tendency with exponent δ controlling the reduction of acceleration
when velocity approaches the desired velocity V0 , and term II represents deceleration
tendency away from the leading vehicle by comparing the current gap dn and the
desired distance d∗ . The form of desired distance d∗ is also a modeling choice, and
we use the expression from book [73] in this paper:
vn (vn − vn−1 )
∗
d (vn , vn − vn−1 ) = d0 + max 0, vn T + √ ,
2 ab
where d0 is the minimum gap to keep, b is the comfortable deceleration, T is the
reaction time of drivers. Larger minimum gap d0 , longer time gap T , or smaller
deceleration b will induce larger desired distance d∗ , and thus larger gap dn for a
driver to keep. The intelligent driver model is simple, and only has a few intuitive
parameters with realistic values and reproduces collective dynamics.
Chapter Four
Dynamics of car-following models
34
4.1 Stability analysis
Instability of traffic flow results in traffic waves, which is also called stop-and-go
waves. This is mainly caused by the delays in adapting the speed to the actual traffic
condition, which are consequences of finite acceleration and deceleration capability.
In this section, we would like to make a stability analysis of the car-following
models without using a Taylor approximation.
4.1.1 Stability analysis of car-following model
The system of N vehicles on a ring road of length L can be describe as:
vn+1 vn vn 1
vehicle n + 1 vehicle n vehicle n 1
xn+1 xn xn 1
`n dn
Figure 4.1: Notations for car-following models: vehicles are numbered from front to back
such that vehicle n follows vehicle n − 1. For vehicle n, xn denotes position, vn denotes
velocity, `n denotes length and dn denotes gap.
When studying the interactions between two adjacent vehicles, the headway
rather than the vehicle length matters. Without loss of generosity, we can set the
`n = 0. By introducing the headway dn between vehicle n − 1 and n, the coupled
ordinary differential equations become:
d˙n (t) = vn−1 (t) − vn (t)
(4.1)
v˙ n (t) = faccel (vn , vn−1 , dn ) ,
35
The steady-state gap is de = L/N and the steady-state velocity ve satisfies faccel (ve , ve , de ) =
0, from which an explicit expression of ve can be obtained as ve = Ve (de ). As in the
analysis for stability, we assume small deviations yn and un from steady state de and
ve .
dn (t) = de + yn (t)
vn (t) = ve + un (t)
(i) By substituting these in, we get odes for the deviations:
dyn ddn
= = un−1 − un
dt dt (4.2)
dun dvn
= = av un + al un−1 + ad yn
dt dt
∂faccel ∂faccel ∂faccel
where av = ∂vn e
|, al = ∂vn−1 e
| and ad = ∂dn e
|. In addition, moving along
the space of steady-state solutions by simultaneously changing dn and vn =
vn−1 must not change the acceleration faccel , so we can obtain a relationship
between partial derivatives in steady-state:
0 = ad dde + (av + al )dve = ad dde + (av + al )Ve0 (de )dde
resulting in
ad = −Ve0 (de )(av + al ) (4.3)
It is assumed that ad > 0, av < 0 and al > 0. Because of equation (4.3), it is
easy to know that av + al < 0.
36
(ii) We can solve equation (4.2) by using Fourier-Ansatz:
yn (t) = yˆeλt+ink , un (t) = uˆeλt+ink
Inserting the ansatzs into odes (4.2) results in
−ik
λ 1−e yˆ
· =0
−ik
−ad λ − (av + al e ) uˆ
.
The resulting solvability condition assumes the form of a quadratic equation
1 p
λ2 + p(k)λ + q(k) = 0 or λ1/2 = −p(k) ± p2 (k) − 4q(k)
2
where p(k) = −av − al e−ik > 0 and q(k) = ad (1 − e−ik ).
p
(iii) Stability condition is equivalent to Re(λ1/2 ) < 0, i.e. Re( p2 (k) − 4q(k)) <
Re(p(k))
Rewrite p2 (k) − 4q(k) as Reiϕ = X + iY , then it becomes
p √
Re( p2 (k) − 4q(k)) = R cos ϕ/2 < Re(p(k)) (4.4)
with
Re(p(k) = −av − al cos(−k)
X = R cos(ϕ) = (av + al cos(−k))2 − a2l sin2 (−k) − 4ad (1 − cos(−k))
Y = R sin(ϕ) = 2al sin(−k)(av + al cos(−k)) + 4ad sin(−k) (4.5)
37
Take squares on inequality (4.4):
R cos2 (ϕ/2) = R(1 + cos(ϕ))/2 < Re(p(k))2
Use cos ϕ = X/R and simplify the inequality as:
R < 2Re(p(k))2 − X
Take squares on both sides and substitute R2 = X 2 + Y 2 :
X 2 + Y 2 < (2Re(p(k))2 − X)2
Remove the items X 2 on both sides:
Y 2 < 4Re(p(k))2 − 4XRe(p(k)) (4.6)
Next we will plug in the formula of Y and Re(p(k)), and simplify the expression
Y 2 − (4Re(p(k))2 − 4XRe(p(k)))
(4.7)
= 16ad [ad sin2 (−k) + (av + al cos(−k))(avd − av )(1 − cos(−k))] < 0
It is assumed that ad > 0, av < 0, al > 0 and av + al < 0. Therefore, the
inequality (4.7) becomes:
av + al cos(−k) 1 − cos(−k)
Ve0 (de ) < (al − av ) 2
av + al
| {z } | sin {z
(−k)
}
I II
Because of av + al cos(−k) ≤ av + al < 0, the part I ≥ 1 with equality holding
if cos(−k) = 1. The part II is ∞ when cos(k) = ±1, and it is 1/(1 + cos(−k))
38
when cos(−k) 6= ±1. Then we can get that II > 1/2 for sure.
In sum, the stability condition is
1
Ve0 (de ) ≤ (al − av ) (4.8)
2
There are three main factors determining the stability of traffic flow: the suf-
ficient sensitivity Ve0 (de ), speed sensitivity av and lead speed sensitivity al .
System of greater anticipation to the lead car (bigger al ) will favor stability,
which means anticipative driver can better adjust deviations and thereby the
system is stable.
(iv) the stability condition for specific models
• Optimal Velocity Model(OVM): Ve = V, av = −1/τ, al = 0
1
V 0 (de ) <
2τ
• Full Velocity Difference(FVDM): Ve = V, av = −1/τ − γ, al = γ
1
V 0 (de ) < +γ
2τ
• Modified Gazis-Herman-Potts Model(GHP): Ve = V, av = −1/τ −1/de , al =
1/de
1 1
V 0 (de ) < +
2τ de
• Intelligent Driver Model(IDM):
(s0 + T ve )( ab ve + aT ) 2a ve 3
p
Ve0 (de ) ≤ +
d2e v0 v0
39
with
s0 + V T
Ve−1 (v) = p
1 − (v/v0 )δ
3
(s0 + T ve )( ab ve + 2aT )
p
4a ve
av = − −
v0 v0 d2e
pa
(s0 + T ve ) b ve
al =
d2e
4.2 Model calibration
In this section, the parameters of the car-following models are calibrated to the
empirical traffic data. This procedure of calibration makes the car-following mod-
els not only realistic to reflect driving behaviors, but also comparable to simulate
microscopic traffic.
The model parameters may fall into two types: (a) parameters that determine
the macroscopic flux-density relation (fundamental diagram); (b) parameters that
do not contribute the flux-density relation, but instead describe the driving habit
and sensitivity: such as adaption time and comfortable acceleration.
4.2.1 Flux-density relation
Flux-density formula
Fundamental diagram of traffic flow is the plot of the traffic flux Q (the number of
passing vehicles per unit time) versus the traffic density ρ (the number of vehicles
40
per unit length). According to [64], the two-phase traffic theory divides traffic flow
into free flow for low densities, and congested flow for large densities. In the free
flow regime, flux increases as density increases, while in the congested flow regime,
flux decrease as density increases. Please see Figure 4.2.
Next we will induce the macroscopic flux-density relation from the microscopic
car-following models through steady-state equilibrium. Firstly, the steady-state
speed ve and steady-state gap de satisfy the following relation:
faccel (ve , ve , de ) = 0 ⇒ ve = Ve (de ) or de = De (ve ),
where De is the inverse function of Ve . The equilibrium formula of the car-following
models are as follows: (a1) the optimal velocity model, the full velocity difference
model, and the modified GHR model; (a2) the intelligent driver model
tanh(dn /∆s − β) + tanh β
(a1) Ve (de ) = V(de ) = V0
1 + tanh β
δ !− 12
ve
(a2) De (ve ) = (s0 + ve T ) 1 −
V0
Secondly, a micro-macro relation between gap de and the macroscopic density ρ can
be built as:
1
ρ = ρ(de ) = ,
de + `
where ` is the average length of vehicles. Thus by combining the relations above,
we obtain the flux-density relation for the car-following models: (a1) the optimal
velocity model, the full velocity difference model, and the modified GHR model;
41
(a2) the intelligent driver model
(a1) Q(ρ) = ρ(de ) × Ve (de )
(4.9)
−1
(a2) Q(ρ) = ρ(de ) × De (de )
Flux-density data
We now apply the formula above to calibrating the car-following models using data
that is provided by California Department of Transportation. The traffic data is
collected by the Caltrans Performance Measurement System (PeMS) [69], which
provides flow, speed and occupancy data across the vehicle detector stations (vds)
in the form of time series over days of operation. These sensors span the the freeway
system across all major metropolitan areas of the State of California.
We recall the results from [24], which used data from a 28 mile stretch of Freeway
I-880S (between 29th Avenue on-ramp in Oakland and the Auto Mall Parkway on-
ramp in Fremont) [24]. The data is aggregated over time intervals of 5 minutes over
98 days that were identified as functioning over 80% between February 2007 and
March 2008. The scatter plot of flux-density (density data is computed based on
flow and occupancy) data is shown in Figure 4.2(a).
The paper [24] presented a method for automated and empirical calibration of
freeway traffic flow characters. The parameters for flux-density relation are shown
in Figure 4.2(b) and Table 4.1, including free-flow speed vfree , maximal flux Qmax ,
congestion speed ω and jam density ρmax .
42
Flow (vehicles per hour per lane)
2500
Flow (vehicles per hour per lane)
2500 slope vfree
2000 bound Qmax
2000
1500
1500 slope w
1000
1000
500 intercept ⇢max
500
0
0 0 50 100 150 200 250
0 20 40 60 80 100 120
Density (vehicles per mile per lane)
Density (vehicles per mile per lane)
(a) Data for fundamental diagram (b) Estimated parameters using fundamen-
tal diagram
Figure 4.2: Automatic calibration of the fundamental diagram from [24]. (a) The scatter of
flux-density data collected by PeMS. (b) The estimation of main parameters by regression.
Description Parameter Value Unit
Free-flow speed vfree 63.3 miles/hour
Maximal flux Qmax 2031 vehicles/hour
Congestion speed param- ω 10.1 miles/hour
eter
Jam density ρmax 232 vehicles/mile
Table 4.1: Parameters from auto-calibration of the fundamental diagram
4.2.2 Model calibration
All parameters for the car-following models in § 3.2 consist of: (a) the parameters
that determine the flux-density relation; (b) the parameters that describe the driving
habit and sensitivity. For type (a), the parameters are estimated by tuning the flux-
density curve to mimic the desired triangle shape in Figure 4.2. The flux-density
relations produced by the car-following models with the calibrated parameters are
shown in Figure 4.3.
For type (b), the parameters take the typical values that were used or estimated
in the research work [37, 50]. In [50], the car-following behavior is studied on the
basis of (publicly available) trajectory datasets recorded by a vehicle equipped with
an radar sensor. The car-following models (except the modified GHR model) were
calibrated by minimizing the deviation between observed dynamics and the simulated
43
2500 2500
2000 2000
Flux (vehicles/hour)
Flux (vehicles/hour)
1500 1500
1000 1000
500 500
0 0
0 50 100 150 200 250 0 50 100 150 200 250
Density (vehicles/mile) Density (vehicles/mile)
(a) Calibration for optimal velocity model (b) Calibration for intelligent driver model
Figure 4.3: Calibration for flux-density relation. (a) The flux-density relation for (a1) in equation
(4.9). (b) The flux-density relation for (a2) in equation (4.9). The blue solid line represents the
flux-density relation for the car-following models, while the red dashed line represents the relation
from PeMS traffic data.
trajectory. The remainder parameter for the modified GHR model is from [37]. The
data used for calibration are the trajectories of 4, 733 vehicles provided through the
FHWA Next Generation Simulation (NGSIM) project.
4.3 Collision behavior
As the main component of traffic systems, car-following models have become of in-
creased importance in traffic simulation and safety research. In [14, 39, 59], the
dynamical behaviors and properties of car-following models have been well studied
and reviewed. Among these traffic dynamics, we are most interested in the colli-
sion behavior (collision-prone behavior and collision-free behavior) of car-following
models, which focuses on the critical situation of potential collisions between two
adjacent vehicles.
The collision behavior of traffic models are closely related to the traffic safety
studies [11, 13, 75], and the collision analysis could act as a guide for microscopic
simulation, especially for lane-changing models [61]. Scholars are interested in the
44
Description Para Value Type/Ref
Optimal velocity model
Adaption time τ 1.5 s (b)[50]
Desired speed v0 68 mph (a)
Transition width ∆s 30 m (a)
Form factor β 0.5 (a)
Vehicle length ` 7m (a)
Others
Sensitivity factor γ 0.65 s−1 (b)[50]
Sensitivity parameter η 12 m/s (b)[37]
Acceleration bound A 3 m/s2 (b)[50]
Intelligent driver model
Desired velocity v0 69 mph (a)
Acceleration a 1.5m/s2 (b)[50]
Acceleration exponent δ 8 (a)
Minimum gap s0 2m (a)
Safe time gap T 1.5 s (a)
Deceleration b 3 m/s2 (b)[50]
Vehicle length ` 5m (a)
Table 4.2: Parameter values for car-following models. Type (a) parameters determine the flux-
density relation, while type (b) parameters describe the driving habit and sensitivity terms.
drivers for collision-free behavior, which could be used to develop the sensitivity
terms in new traffic models [18]. In addition, artificial intelligence, as the ultimate
goal of traffic studies, would make a good use of these collision-free behavior to
enhance the cooperative collision avoidance in highways [12]. However, evaluating
the collision behavior of different car-following models is a complex and crucial task,
which has to determine the effects of specific components and parameters. Thus,
most assertions of the collision behaviors are based on the microscopic simulations
rather than the theoretical analysis.
In this section, we would like to carry out a theoretical analysis of the collision
behavior of four well-known continuous car-following models: the optimal velocity
model [8], the full velocity difference model [46], the modified GHR model [32] and
the intelligent driver model [72]. There are two main goals of the work: (1) identify
whether the car-following models are collision-prone or collision-free with theoretical
45
analysis; (2) simulate the collision behavior of the car-following models for a com-
parison between theoretical and computational results. Specifically, we will show
that the optimal velocity and the full velocity difference models are collision-prone
regardless of the parameters, and that the modified GHR and intelligent driver mod-
els are unconditionally collision-free regardless. In addition, our study discovers that
the simulation is not valid to assert the collision behavior, since the numerical errors
could introduce collisions even though the theoretical results do not support them.
4.3.1 Simulation of two-vehicle system
Collision behavior can be described as collision-prone or collision-free: collision-prone
means that the vehicles generate collisions in critical scenarios, while collision-free
means that the vehicles can avoid collisions appropriately in any critical scenario.
In the context of the car-following models in (3.1): a model is called collision-free
if collision (i.e. xn (t) = xn+1 (t) for two adjacent vehicles n and n + 1 at some
time t) cannot occur; a model is called collision-prone if collisions could occur under
some circumstances. In this section, we provide an intuition of the collision behavior
through the phase portrait analysis.
Collision simulations
Consider a simple scenario: there are two vehicles on a road. The leading vehicle is
moving with a constant speed c, and the other vehicle is following it according to
a car-following model. This traffic system can be simplified by introducing the gap
46
d2 = x1 − x2 and the speed difference w2 = v1 − v2 :
x˙ 1 (t) = c,
d˙2 (t) = w2 (t),
v˙ (t) = 0,
1
=⇒ (4.10)
w˙ 2 (t) = −faccel (c − w2 , c, d2 ),
x˙ 2 (t) = v2 (t),
v˙ 2 (t) = faccel (v2 , v1 , d2 ),
We are interested in the behavior when the following vehicle approaches the leading
one with a faster speed (i.e. d2 (0) > 0, w2 (0) < 0). If the gap d2 (t) between vehicles
goes negative some time, we claim that the collision-prone behavior is observed.
In our simulation, we take the parameters for the car-following models from Ta-
ble 4.2, and assume that the leading vehicle moves with a constant speed c = 10mph.
The phase portraits (d2 , w2 ) of the car-following models are shown in Figure 4.4, with
the collision boundary d2 = 0 highlighted by the red dashed line. We can conclude:
• It is obvious to observe that the optimal velocity model and the full veloc-
ity difference model are collision-prone in this simply two-vehicle system. In
the following subsection, we will justify that these two models are generally
collision-prone independent of parameter settings.
• Collisions are not observed in the modified GHR model and the intelligent
driver model in this simple two-vehicle system. In §4.3.2, a rigorous proof will
be provided to show that these two models are collision-free in multiple-vehicle
system.
47
30 30
20 20
10 10
w2 (mph)
w2 (mph)
0 0
-10 -10
-20 -20
-30 -30
-10 0 10 20 30 40 50 -10 0 10 20 30 40 50
d2 (m) d2 (m)
(a) Optimal velocity model (b) Full velocity difference model
30 30
20 20
10 10
w2 (mph)
w2 (mph)
0 0
-10 -10
-20 -20
-30 -30
-10 0 10 20 30 40 50 -10 0 10 20 30 40 50
d2 (m) d2 (m)
(c) Modified velocity model (d) Intelligent driver model
Figure 4.4: Phase portraits of car-following models. The x-axis represents the gap d2 between
vehicles, and the y-axis represents the speed difference w2 . The red dashed line d2 = 0 represents
the critical situation of collisions.
48
Collisions in the optimal velocity and full velocity difference models
The two-vehicle system (4.10) is specified for the optimal velocity model (with γ = 0)
and the full velocity difference model:
d˙2 = w2 ,
(4.11)
w˙ 2 = − 1 (V(d2 ) + w2 − c) − γw2 .
τ
The phase portrait of the system is sketched in Figure 4.5. The trajectory through
the origin is critically collision-free (d2 ≥ 0, but d2 (t) = 0 for some time t). For the
trajectories to its left, the gap d2 will definitely go negative, and thus generating
collisions. The region (d2 , w2 ) that could lead to negative gap is called the collision
region, which is the area to the left of the critical trajectory in Quadrant IV.
30
20
10
w2 (mph)
0
-10
-20
-30
-10 0 10 20 30 40 50
d2 (m)
Figure 4.5: Analysis of collision-prone behavior in optimal velocity model and full velocity
difference model.
Let us investigate the field direction along the y axis, i.e. the collision boundary
d2 = 0:
49
• Origin (d2 = 0, w2 = 0): the filed direction is up because of (d˙2 , w˙ 2 ) = (0, τc >
0). Along the x axis with small gap d2 < V −1 (c), the field direction is also
upright. Please see the black arrows in Figure 4.5.
• Negative y axis (d2 = 0, w2 < 0): the filed direction is up-left because of
(d˙2 , w˙ 2 ) = w2 < 0, τc − ( τ1 + γ)w2 > 0 . Here the constants c, τ, γ are all pos-
itive. Thus, the gap d2 will go negative soon and collisions occur. Please see
the magenta arrows in Figure 4.5.
• Positive y axis near the origin (d2 = 0, w2 > 0): the filed direction is always
to the right because of d˙t = w2 > 0. Thus the two-vehicle system goes back to
the normal situation from collisions. Please see the green arrows in Figure 4.5.
Based on the analysis above, the optimal velocity model and full velocity differ-
ence model are collision-prone, independently of parameter settings. However, the
collision region can be reduced through appropriate parameter settings, such as the
sensitivity parameter γ. By adding the speed sensitivity −γw2 , the full velocity dif-
ference model has a much smaller collision region compared to the optimal velocity
model (see Figure 4.4). This is why the full velocity model is often claimed to be
collision-free in simulation [73] even though it is collision-prone in theory.
4.3.2 Proof of collision-free behavior
In this section, we aim to prove that the modified GHR model and the intelligent
driver model are collision-free. The analysis approach developed in this paper can
be applied to other car-following models to study their collision behavior. There
are two traffic scenarios studied here: the collision-free proof for traffic on a straight
road is presented first, while the last part focuses on traffic on a ring road.
50
The collision behavior of the car-following models can be concluded as:
Theorem 4.1. The optimal velocity model (3.2) and the full velocity difference model
(3.3) are collision-prone, which is independent of parameter settings.
Theorem 4.2. The modified GHR model (3.4) and the intelligent driver model (3.5)
are collision-free under all traffic scenarios no matter it is a straight road or a ring
road.
Traffic on a straight road
There are n vehicles on a straight road: vehicles are labeled 1, 2, . . . , n from front to
back, and the leading vehicle 1 is moving at a constant speed c. As below, we define
xi : the position of vehicle i, i = 1, . . . , n
vi : the speed of vehicle i, i = 1, . . . , n
di = xi−1 − xi : the headway of vehicle i, i = 2, . . . , n
wi = vi−1 − vi : the speed difference of vehicles i and i − 1, i = 2, . . . , n
The full system of n vehicles can be expressed in terms of xi and vi , which can be
also written in terms of di and wi :
d˙2 = w2
x˙ 1 = c
v˙ 1 = 0
w˙ 2 = −faccel (c − w2 , c, d2 )
x˙ 2 = v2 =⇒
d˙3 = w3 (4.12)
v˙ 2 = faccel (v2 , c, d2 )
w˙ 3 = faccel (v2 , v1 , d2 ) − faccel (v3 , v2 , d3 )
. . .
· · ·
51
The idea behind the proof is mathematical induction. Firstly, we show that the
system of 2 vehicles is collision-free. We then show that system of j vehicles is also
collision-free given the system of j − 1 is known to be collision-free, which is the
inductive step. In the proof, we only consider the critical scenario that the following
vehicle is approaching the leading vehicle with a faster speed but a small gap ahead.
Otherwise collisions never happen.
Modified GHR model As mentioned earlier, we focus on the critical scenario
when (1a) wj < 0 and the gaps are small. Then the full system for the modified
GHR model can be simplified as:
d˙2 = w2
w˙ 2 = −η wd22 − τ1 (V(d2 ) + w2 − c)
d˙3 = w3 (4.13)
w˙ 3 = −η wd33 + η wd22 + τ1 (V(d2 ) − V(d3 ) − w3 )
· · ·
where the stimulus −η wdii can provide strong repulsions when the gaps di are very
small, and thus these terms will dominate in a critical scenario. We use the technique
of fast system to analyze the modified GHR model by substituting dj = εqj and
52
t = εT :
q20 = w2
w20 = −η wq22 − τε (V(εq2 ) + w2 − c)
0
q3 = w 3 (4.14)
w3 = −η wq33 + η wq22 + τε (V(εq2 ) − V(εq3 ) − w3 )
0
· · ·
where the terms contain ε can be treated as small perturbations. When ε = 0, the
system (4.14) becomes the fast system:
q20 = w2
w20 = −η wq22
q30 = w3 (4.15)
w30 = −η wq33 + η wq22
· · ·
Next we will divide the proof into two parts: firstly we would like to show that
the fast system (4.15) is collision-free; then we go back to the original system (4.14),
and show that the perturbations do not impact the collision-free property.
(i) The fast system could be generally written as
Q0j = Wj ,
(#)
Wj0 = αj Wj−1 − βj Wj ,
Qj−1 Qj
αj ≥ 0, βj > 0
53
We start with a general form of fast system with non-negative parameters αj
and positive parameter βj . These parameters can be specified based on the
car-following models. The solution of system (#) satisfies:
Q0j (#) Wj Wj−1 Wj
βj + Wj0 = βj + αj − βj
Qj Qj Qj−1 Qj
Wj−1
= αj
Qj−1
Q0j−1
= αj
Qj−1
Taking the integral on both sides gives:
d d
(βj ln Qj + Wj ) = (αj ln Qj−1 ), (4.16)
dt dt
Thus during the period of (1a) wj < 0, we conclude the following:
αβj
Wj (0) − Wj (T )
Qj−1 (T ) j
Qj (T ) = Qj (0) exp
Qj−1 (0) βj
αβj (4.17)
Qj−1 (T ) j Wj (0)
≥ Qj (0) exp
Qj−1 (0) βj
Next we will show that the gaps {Qj (T )}j=2,3,··· ,n never go negative during the
critical traffic scenario given the positive initial conditions {Qj (0) > 0}j=2,3,··· ,n :
• For the second vehicle with α2 = 0, β2 = η, the gap Q2 satisfies:
W2 (0)
Q2 (T ) ≥ Q2 (0) exp = B2
η
which guarantees that the second vehicle never collides with the leading
vehicle, and the positive lower bond is B2 .
54
• For the third vehicle with α3 = β3 = η, the gap Q3 satisfies:
Q2 (T ) W3 (0)
Q3 (T ) ≥ Q3 (0) exp = B3
Q2 (0) η
where the gap for the second vehicle Q2 (t) has been proved positive for
any time T no-negative. Thus, the third vehicle never collides with the
second vehicle in a critical scenario.
• ···
• For the last vehicle n, we can similarly prove that it will not collide with
its leading vehicle n − 1. Thus, by using this induction method, the fast
system (4.15) is surely collision-free.
(ii) Next we analyze the perturbation terms in equation (4.14) and evaluate their
impact in the collision-free property. As mentioned earlier, collisions are only
possible when the gaps are very small 0 < {qj }j=2,3,··· ,n 1 . Thus, we study
the perturbation terms only when the gaps are less than the predetermined
thresholds.
• For the second vehicle, the perturbation term is in bold as follows:
q20 = w2
w20 = −η w2 − ε (V(εq2 ) + w2 − c),
q2 τ
c≥0
The threshold for q2 is (1b) 0 < εq2 ≤ V −1 (w2 ) , which defines a corner
area in Quadrant IV enclosed by the curve q2 = 0, w2 ≤ 0 and the curve
V(εq2 ) + w2 = 0. Then the perturbation term satisfies:
ε (1b) ε
− (V(εq2 ) + w2 − c) ≥ − (0 − c) ≥ 0
τ τ
55
Take the differential of trajectory w2 (q2 ):
dw2 w0 η 1 ε (1a)(1b) η
= 02 = − + − (V(εq2 ) + w2 − c) ≤ − (4.18)
dq2 q2 q2 w 2 τ q2
The slope along the trajectory w2 (q2 ) is negative, which is steeper than
that of the fast system (4.15). The comparison of the slopes is demon-
strated in Figure 4.6. Once a trajectory w2 (q2 ) enters the threshold area
(defined by 0 < εq2 ≤ V −1 (w2 )), it will cross the trajectories of the fast
system nearby to the upright direction. As proved earlier, the fast system
is collision-free, and its trajectories never go across the collision boundary
dj = 0. Therefore, it is easy to observe that the original system is more
conservative than the fast system under the critical circumstance.
wi
system 1
system 2
slopes
0 di
Figure 4.6: Slope comparison. System 1 (dashed green line) represents the fast system
(4.15) and system 2 (solid blue) represents the original car-following system (4.14).
• For the third vehicle, the perturbation term is in bold as follows:
q30 = w3
w30 = η w2 − η w3 + ε (V(εq2 ) − V(εq3 ) + w3 )
q2 q3 τ
56
It is known that the second vehicle never collide, i.e. q2 (T ) > 0, so we set
the threshold for gap q3 as (1c) 0 < q3 ≤ q2 . Otherwise, collisions can not
occur. Then the perturbation term satisfies:
ε (1c) ε (1a)
(V(εq2 ) − V(εq3 ) − w3 ) ≥ w3 ≥ 0
τ τ
Take the differential of trajectory w3 (q3 ):
dw3 w0 ηw2 η 1 ε
= 03 = − + (V(εq2 ) − V(εq3 ) − w3 )
dq3 q3 w 3 q 2 q3 w 3 τ
(1a)(1b) ηw
(4.19)
2 η
≤ −
w3 q2 q3
The slope along the trajectory w3 (q3 ) is negative and steeper than that of
the fast system. Thus, we could infer that the third vehicle never collide
based on the similar analysis and comparison of slopes in Figure 4.6.
• ···
• For the last vehicle n, we can prove the same results. Thus, the original
system (4.14) of the modified GHR model is collision-free.
In summary, the modified GHR model has been proved collision-free on a straight
road.
Intelligent driver model The analysis for the intelligent driver model is similar
to that in the previous section. We reform the formula of the intelligent driver
model under the critical traffic scenario (1a) wj < 0, and compare it with the proved
collision-free fast system.
57
The full system for the intelligent driver model can be simplified as:
d˙2 = w2
δ ∗ 2
s (v ,−w )
v
w˙ 2 = a −1 + V0 +2 2 2
d2
d˙3 = w3 (4.20)
δ δ ∗ 2 ∗ 2
v v s (v ,−w ) s (v ,−w )
w˙ 3 = a V0 − V0 + −
3 2 3 3 2 2
d3 d2
· · ·
As mentioned earlier, collision are only possible when the gaps are very small 0 <
dj 1, so we study the system only when the gaps are less than the predetermined
thresholds.
• For the second vehicle, we only consider the circumstance of small gap with
(2b) d2 < min(1, d0 /2). The acceleration of the speed difference can be bounded
from below:
δ √ !2
d0 + v2 T − v2 w2 /2 ab
(4.20) v2
w˙ 2 = a −1 + +
V0 d2
√ !2
(2b) s0 s0 v2 T v2 w2 /2 ab
≥ a −1 + + + +
2d2 2d2 d2 d2
2 √ !2
(2b) s0 v2 w2 /2 ab
≥ a + (4.21)
2d2 d2
2 2
(2b) s v 1 + w22
≥ a min 0 , 2
4 4ab d2
2 2
s v −2w2
≥ a min 0 , 2
4 4ab d2
w2
= −β2 ,
d2
where the right-hand-side is a formula of fast system as in (#) by introducing
58
s20 v22
α2 = 0 and β2 = 2a min ,
4 4ab
. Take the differential of trajectory w2 (d2 ):
dw2 w˙ 2 (4.21)(2a) β2
= ≤ − (4.22)
dd2 d˙2 d2
The slope along the trajectory w2 (d2 ) is negative and steeper than that of the
fast system. Thus the second vehicle never collides based on the analysis and
comparison of slopes in Figure 4.6.
• For the third vehicle, we only consider the circumstance of small gap with
(2c) d3 < min(d2 , 1, d0 /2). The acceleration of the speed difference can be
bounded from below:
δ δ √ !2 √ !2
(4.20) v3 v2 d0 + v3 − v3 w3 /2 ab d0 + v2 − v2 w2 /2 ab
w˙ 3 = a − + −
V0 V0 d3 d2
√ !2 √ !2
(2a) d0 + v3 − v3 w3 /2 ab d0 + v2 − v2 w2 /2 ab
≥ a −
d3 d2
" √ ! √ !#
(2c) d0 + v3 − v3 w3 /2 ab d0 + v2 − v2 w2 /2 ab
≥ 4a −
d3 d2
√ √ !
(2c) d0 + v3 d0 + v2 −v3 w3 /2 ab −v2 w2 /2 ab
≥ 4a − + −
d3 d3 d3 d2
√ √ !
(2a) −v3 w3 /2 ab −v2 w2 /2 ab
≥ 4a −
d3 d2
w2 w3
= α3 − β3
d2 d3
(4.23)
where the right-hand-side is also a formula of the fast system as in (4.15) by
√ √
2v2 a 2v3 a
introducing α = √
b
and β3 = √ .
b
Take the differential of trajectory
59
w3 (d3 ):
dw3 w˙ 3 (4.23)(2a) w2 1
= ≤ α3 − β3 (4.24)
dd3 d˙3 w 3 d2 d3
The slope along the trajectory w3 (d3 ) is negative and steeper than that of the
fast system. Thus the third vehicle never collide based on the analysis and
comparison of slopes in Figure 4.6.
• ···
• For the last vehicle n, we can prove the same results. Thus, the original system
(4.20) of the intelligent driver model is collision-free.
Traffic on a ring road
There are n vehicles on a ring road with the same notations as on a straight road.
Unlike the straight road scenario, there is no leading vehicle moving at a constant
c, and each vehicle can be influenced by the rest of vehicles rather than only the
leading vehicles. Despite the increased complexity of the ring road scenario, similar
mathematical analysis can be applied to prove the collision-free property.
The full system of n vehicles can be expressed in terms of xi and vi , which can
60
be also written in terms of di and wi :
d˙1 = w1
x˙ 1 = v1
v˙ 1 = faccel (v1 , vn , d1 )
w˙ 1 = faccel (vn , vn−1 , dn ) − faccel (v1 , vn , d1 )
x˙ 2 = v2 =⇒
d˙2 = w2
v˙ 2 = faccel (v2 , v1 , d2 )
w˙ 2 = faccel (v1 , vn , d1 ) − faccel (v2 , v1 , d2 )
. . .
· · ·
(4.25)
The idea behind the proof is mathematical induction. Specifically, the proof can
be divided into three steps:
(i) The vehicles are on a ring road, so the sum of the gaps is equal to the length
of road L all the time:
d1 (t) + d2 (t) + · · · + dn (t) = L
The averaged gap is L/n. At time t0 , there is at least one vehicle k0 , whose
headway is above the averaged gap L/n (Pigeonhole Principle). Since the
vehicle speed is bounded by the maximal speed V0 , there must exist a time
period ∆ = L/nV0 such that the vehicle k0 will not collide with its leader.
Thus, vehicle k0 is collision-free during time [t0 , t1 ], where t1 = t0 + ∆.
(ii) Given a collision-free vehicle k0 is collision-free, we intend to show that its
following vehicles are all collision-free during the period [t0 , t1 ]. Firstly, we
consider its following vehicle k0 + 1. For any time t ∈ [t0 , t1 ], the following
vehicle k0 + 1 must stay in one of the following circumstances:
61
(a) wk0 +1 ≥ 0 or dk0 +1 ≥ min(dk0 , d0 /2, 1): there is no need to worry about
collisions between vehicle k0 and k0 + 1. There is because positive speed
difference wk0 +1 is increasing the gap dk0 +1 , and that collisions cannot
occur with the gap dk0 +1 larger than the predetermined threshold.
(b) wk0 +1 < 0 and 0 < dk0 +1 < min(dk0 , d0 /2, 1): according the proof in
§4.3.2, the slope of trajectory wk0 +1 (dk0 +1 ) is negative and steeper than
that of the fast system (#), which has already shown collision-free. Thus,
the vehicle k0 + 1 never collides based on the comparison of slopes in
Figure 4.6.
By mathematical induction, we can repeat this process and show that all ve-
hicles are collision-free during the period [t0 , t1 ].
(iii) By repeating step 1 and 2, we can show that the system is collision-free during
[t1 , t2 ], [t2 , t3 ], · · · . The length of each period is ∆ = L/nV0
(iv) It is obvious ∆ + ∆ + ∆ + · · · = ∞, so the traffic system is theoretically
collision-free forever.
4.3.3 Numerical results
The main goal is to study the collision-prone and collision-free behaviors of car-
following models: the optimal velocity model and the full velocity difference model
are shown collision-prone independent of parameter settings through the phase por-
trait analysis in §4.3.1; while the modified GHR model and the intelligent driver
model are shown collision-free through mathematical analysis in §4.3.2.
In this section, we carry out the traffic simulations to test the consistency of
62
collision behaviors to the theoretical analysis. The parameter settings are based on
the calibration results in §4.2. In addition, we investigate how the numerical errors
potentially influence the collision behaviors, and discover that the numerical errors
can introduce collisions even though the model does not support them.
Monte Carlo simulation: counts of collisions
We put 28 vehicles on a ring road of 350m length: the initial positions and velocities
are randomly picked according to the uniform distribution; with the same initial
conditions, the vehicles will move 20 seconds based on the four car-following models.
In order to make the results representative of the collision behaviors, we implement
the Monte Carlo simulation with 500 iterations. We count the number of collisions
by analyzing the vehicle trajectories, and the resulting distribution is in Figure 4.7.
• The averaged number of collisions in the optimal velocity model reaches as
many as 25.85, while it is only 3.87 for the full velocity difference model.
This is consistent to our phase portrait analysis in §4.3.1. For the full velocity
difference model, it is uncommon to observe collisions in simulation if the initial
conditions are not extremely critical. This is why this model is often claimed
to be collision-free in simulation [73].
• There are a few collisions in the modified GHR model and the intelligent driver
model, even though they are collision-free theoretically. This is because the
numerical errors bring collisions that are not supposed to exist. Thus, the col-
lision behaviors in simulation cannot used to conclude whether a car-following
model is collision-prone or collision-free.
63
14 30
12
25
10
20
Frequency (%)
Frequency (%)
8
15
6
10
4
2 5
0 0
5 10 15 20 25 30 35 40 45 0 1 2 3 4 5 6 7 8 9
Number of Collisions Number of Collisions
(a) Optimal velocity model (b) Full velocity difference model
80 80
70 70
60 60
Frequency (%)
Frequency (%)
50 50
40 40
30 30
20 20
10 10
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5 3
Number of Collisions Number of Collisions
(c) Modified GHR model (d) Intelligent driver model
Figure 4.7: Trajectories of car-following models. The x-axis represents the bins of collisions
observed, and the y-axis represents the counts of collision in percentage.
64
Why can collisions occur in simulations of collision-free models?
It is an interesting phenomena that collisions still numerically exist in the theoreti-
cally collision-free car-following models. An intuition for this comes from the phase
portraits in Figure 4.4: when the trajectory gets close to the y-axis (di = 0), it is
easy to be pushed across the y-axis by numerical errors. Therefore, we are interested
in studying whether the collisions in the simulations of collision-free car-following
models could diminish or even disappear with finer meshes of time.
Monte Carlo simulation The basic setting is the same: there are28 vehicles
on a ring road of 350m length with randomly picked initial conditions. We run
the numerical scheme for different time step h, and repeat this process with 4000
iterations for different initial conditions. The meshes of time are getting finer and
finer according to the following table:
Total time T 20s 20s 20s ··· 20s
Mesh points n 400 800 1200 ··· 6000
Time step h = T /n 0.0500s 0.0250s 0.0167s ··· 0.0033s
Table 4.3: Setting of time steps in Monte Carlo simulation
There are two ways to record the trajectories of vehicles, and thus there are two
ways to count collisions respectively:
(i) Update the positions of vehicles every h time, but only record the trajectories
at the fixed time points [0.05s, 0.10s, 0.15s, · · · , 20s], i.e. the recording time
interval is ∆T = 0.5s. Then the number of collisions are counted based on the
trajectories that are linearly interpolated among these 400 fixed time points.
(ii) Update the positions of vehicles every h time, and record the trajectories every
65
time step until 20s [h, 2h, 3h, · · · , nh = 20s], i.e. the recording time interval is
∆T = h. Then the number of collisions are counted based on the trajectories
that are linearly interpolated among these n = 20/h time points.
The main difference between the two methods above is the time interval ∆T for
recording trajectories and measuring collisions. The intuition tells that the more
frequent we check vehicles’ positions, the more collisions we observe. The Monte
Carlo simulations are implemented in the modified GHR model with Euler explicit
scheme, and the results are shown in Figure 4.8.
0 -0.5
Monte Carlo Monte Carlo
Linear Regression (slop=1.932) Linear Regression (slop=0.9981)
-1 -1
log(number of collisions)
log(number of collisions)
-2 -1.5
-3 -2
-4 -2.5
-5 -3
-6 ! ! -3.5 !
-6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 -6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5
log(time step h) log(time step h)
(a) log-log plot with ∆T = 0.05 (b) log-log plot with ∆T = h
Figure 4.8: The relationship of collision counts over time step h. The red dots represent the
simulation results, and the blue line is a fitted line from linear regression.
The simulation results are valid to show that the collisions caused by the numer-
ical errors will diminish with smaller time step h. Thus, by choosing a sufficiently
small time step h in simulations will effectively avoid these unrealistic collisions from
the collision-free car-following models. In Figure 4.8, we also observe that the num-
ber of collisions is approximately proportional to the exponentiation of time step h:
for the first counting method with ∆T = 0.05s, Counts ∝ O(h2 ), and for the second
counting method with ∆T = h, Counts ∝ O(h1 ).
Even thought it is hard to figure out how the collisions are related to numerical
66
errors, we guess that they are proportional to the truncation error of the Euler explicit
scheme. The local truncation error of Euler explicit scheme is h2 , which coincides
with the first counting method, and the global truncation error is h1 , which also
coincides with the second counting method. An descriptive explanation might be:
for a numerical scheme with accuracy O(hp ), the local truncation error is proportional
to O(hp+1 )), which directly affects the collision behaviors in one time step h. Then by
checking M points in the vehicle trajectories, we could observe M O(hp+1 ) collisions
in the simulations. For the first counting method with time interval ∆T = 0.05, we
check M = 400 equally distributed points along the trajectories, and thus observe
around 400O(hp+1 ) = O(hp+1 ) in simulations. For the second counting method with
time interval ∆T = h, we check n = 20/h equally distributed points along the
trajectories, and thus observe around 20/hO(hp+1 ) = O(hp ).
Investigation for Euler explicit scheme In the previous part, we come up
with a guess that the number of observed collisions is related to the accuracy of
the numerical scheme. The general proof should be complicated and beyond the
scope of this paper. In this part, we would like to provide a brief discussion about
Euler explicit scheme in a simplified traffic situation. Here we only consider the
second counting method with a time interval ∆T = h, and show that the number of
collisions is proportional to O(h1 ).
We only put two vehicles on a straight road with the leading vehicle moves at
a constant speed c, and study what (d, w) values can generate collisions. In order
to have a collision within a time step h, the following relation d(t) + w(t)h < 0
must hold. This relation sketches a triangle region, and a collision will definitely
occur within one time step h if current (d, w)) falls into it (red checks in Figure 4.9).
In addition, any initial conditions (d(0), w(0)) whose trajectory interacts with this
67
triangle area will definitely generate a collision, too. Thus, the (d, w) values that
generates collisions are in the triangle region plus the supplementary region (red
diamonds in Figure 4.9).
0 d
c
w
V0
(V0 c)h
Figure 4.9: Region partition of initial condition (d, w): the red checks and diamonds area is
collision-prone; the blue lines area is collision-free.
It is observed that the area of this (d, w) region depends on the picked time step
h, so we could study how the collisions diminishes by studying how the region area
shrink when time step h decreases. The regions are sketched in red for different time
step h in Figure 4.10. By linear regression, the region area is proportional to O(h1 ).
4.3.4 Discussion
We aim to study the collision behavior of the car-following models. Specifically, we
studied four well-known continuous car-following models: the optimal velocity model,
the full velocity difference model, the modified GHR model, and the intelligent driver
model. The parameter setting is based on the model calibration using historical
traffic data. The optimal velocity model and the full velocity difference model has
been shown collision-prone independent of parameter settings through the phase
portrait analysis in §4.3.1, while the modified GHR and the intelligent driver model
68
d
h1 h2 h3
w
h4 h5 h6
h7 h8 h9
Figure 4.10: (d, w) collision region for different time step h (h1 : h2 : · · · : h8 : h9 = 1 : 2 : · · · : 8 :
9).
are proved to be collision-free through mathematical analysis in §4.3.2.
In addition to the theoretical analysis of the collision property, we carry out ex-
periments to simulate the trajectories of the car-following models. The simulation
results are consistent with our theoretical analysis except that numerical errors intro-
duce collisions that the model does not support. More simulations are implemented
to show that the number of the observed collisions is proportional to the truncation
error of the numerical scheme.
The study of collision behavior provides a theoretical reference for the car-
following models. The collision behavior will be taken into consideration when de-
signing new models or studying artificial intelligence. In the future work, we will
combine these collision-free models with lane-changing model to construct a complete
traffic system.
Chapter Five
Data assimilation for traffic flow
estimation
70
5.1 Introduction
Estimating the traffic state on a network of highway is a challenging problem. With
the availability of stationary sensor data, traffic cameras, and GPS cell phone data,
more and more data can now be used to predict car density, velocity and travel
time in real time. Predicting traffic involves (i) a model that provides the traffic at
later times based on given initial condition and (ii) techniques to incorporate and
assimilate actual traffic-state observations into the initial conditions to improve the
model forecast. To formulate a framework for estimating traffic flow, we therefore
need a mathematical model of traffic flow, a sense of what types of observations are
available, and a scheme for assimilating these observations into the model. We now
discuss these three ingredients in turn.
Mathematical models for traffic flow come in many different flavors. Common
models range from cellular automata or microscopic car-following models for indi-
vidual cars to macroscopic discrete or partial-differential equation (PDE) models for
car densities, velocities, and possibly other quantities [39, 54, 63]. Knowledge of the
positions and velocities of all cars on a road or road network is required to success-
fully employ a microscopic model. Since most of this information is not well known,
it is typically not feasible to use such models for traffic flow prediction. Thus,
macroscopic models are used more frequently for traffic estimation. A commonly
used minimal macroscopic PDE model is the Lighthill–Whitham equation [54]. This
model is discussed in more detail in §2.
Traffic-flow observations come from a variety of sources and are available as
functions of time. Stationary sensors, such as induction loops or cameras, provide
the flux, average velocity, and local density of cars that move past the fixed sensor
71
location. GPS data from cell phones or navigation devices, on the other hand,
provide information about the positions and velocities of individual cars that move
with the traffic flow. We refer to observations that come from a fixed observation
location as Eulerian observations and to observations that come from parcels (cars)
that move with the traffic flow as Lagrangian observations.
There are a number of different data-assimilation techniques available that help
incorporate traffic observations at discrete times into an underlying model. Most of
these are based on a Bayesian statistics analysis that treats the forecast from the
model as the prior distribution and then calculates a posterior distribution based on
the available observations. Commonly, traffic estimators employ ensemble Kalman
filters for this calculation. The main limitation of the existing data-assimilation
techniques stems from the fact that sensor and GPS data require very different
assimilation approaches: it is therefore technically involved to use existing techniques
to assimilate Eulerian and Lagrangian observations simultaneously.
The main contribution of our research is to propose a method that provides
an efficient alternative approach to the assimilation of Lagrangian GPS data, whilst
keeping the ability to simultaneously assimilate Eulerian sensor data. In this part, we
will discuss the efficacy and accuracy of the proposed method for assimilating sensor
and GPS data in traffic flow. Specifically, we will demonstrate that the proposed
algorithm is robust with respect to different traffic scenarios and its implementation
using ensemble Kalman filter or particle filter. Thus, our method provides a viable
alternative to existing data assimilation techniques in their ability to assimilate both
Eulerian and Lagrangian data.
72
5.2 Underlying macroscopic traffic model
We use the viscous Lighthill-Whitham equation as the underlying core macroscopic
mathematical model of traffic flow. The velocity v and density ρ is linearly related
by the Greenshields function.
ρt + (ρVG (ρ))x = ερxx , (5.1)
ρ
VG (ρ) = vmax 1 − , (5.2)
ρmax
where vmax is the maximal free-flow velocity. Traffic jams manifest themselves as
shock waves in (5.1), see Figure 5.1, which travel with characteristic speed c, where
positive speeds correspond to propagation of traffic jams against the traffic flow.
40 (a) t = 1 min 36 (b) t = 60 mins (c) t = 120 mins
33
36 34 Density(cars/mile)
Density(cars/mile)
Density(cars/mile)
32
32
32 31
30
30
28
28
29
24 26
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Position(miles) Position(miles) Position(miles)
Figure 5.1: Evolution into a shock in the viscous Lighthill–Whitham equation (5.1). The shock
is moving backward with gradually decreasing amplitude.
Given a solution ρ(x, t) of (5.1), we can recover the movement of individual cars.
Indeed, if (pi (t))i=1,2,···Nc denote the positions of Nc individual cars as functions of
time, then these functions satisfy the ordinary differential equations
dpi
= V (ρ(pi , t)), i = 1, 2, · · · , Nc , (5.3)
dt
which can be solved to provide car positions for all times given their initial locations
at t = 0.
73
The macroscopic model (5.1) can be used to describe normal traffic flow in the
absence of ramps, traffic lights, and construction zones. In §2, we have added ad-
ditional terms to the model (5.1) to capture the effects of on- and off-ramps, traffic
lights, road construction sites, and traveling bottlenecks. The modifications allow
us to simulate common traffic scenarios and therefore enable us to better test and
validate the data assimilation strategies proposed here. In the next section, we will
focus on (5.1) but stress that the techniques outlined in this section will also apply
to the modifications of Lighthill-Whitham model.
5.3 Eulerian and Lagrangian observations
There are two main categories of observation models in traffic state estimation: one
is Eulerian observation in which data is collected from stationary sensors on the
road such as cameras, induction loops, and radar; the other type is Lagrangian
observation, which includes the position and velocity data from moving vehicles
equipped with GPS devices or cell phones.
As explained earlier, it is technically involved to use existing techniques to as-
similate Eulerian and Lagrangian observation simultaneously. If we use a partial
differential equation for the car density ρ(x, t) at position x and time t as our un-
derlying model, then sensor data produce estimates for the car density ρ(x0 , t) at
the sensor location x0 , which can be compared directly to the prediction made by
the model. GPS data, on the other hand, are Lagrangian in nature, that is, they
provide the location x(t) and velocity v(t) of an individual car as it moves with the
traffic but cannot produce an estimate of the car density at that moving location.
In particular, such data cannot be assimilated directly into a density-based model.
74
Two approaches have been considered to accomplish the assimilation of La-
grangian data. First, the underlying model can be written in terms of Lagrangian
coordinates [79, 80] that move with the speed of the underlying cars, which makes
it possible to assimilate Lagrangian observations. However, it is then much more
difficult to assimilate Eulerian observations (such as sensor data); moreover, the
approach outlined in [79, 80] was applied only to certain types of Eulerian observa-
tions (namely the velocities of cars measured at fixed sensor locations, but not the
densities or fluxes recorded at these fixed locations). Second, one can convert the
partial differential equation for the density to an equation for the velocity by in-
verting the density-velocity function V (ρ) and then assimilate velocity observations
[76, 77]; however, this inversion may not be possible for other models.
We now introduce our proposed approach that can efficiently assimilate Eulerian
and Lagrangian data simultaneously.
5.3.1 Continuous observation models
Eulerian observation data: Sensors deployed at fixed positions along roads are
capable of counting the number of vehicles passing by as well as reporting the flux
at their location as a function of time. Thus, if there are Ns sensors positioned at
the fixed locations (qi )i=1,2,··· ,Ns on a ring road T, then the complete system of traffic
model (M) and Eulerian observation operators (D) is given by
(M) ρt (x, t) + ∂x ϕ(ρ(x, t)) = ερxx (x, t), x∈T
(5.4)
ϕ
(D) h (ρ) := (ϕ(ρ(qi , t)))i=1,2,··· ,Ns .
75
An advantage of Eulerian sensor data is their high accuracy: modern traffic sensors
are so sensitive that the probability of missing or over-counting a vehicle is below
0.1% ([1]). However, the Eulerian data is expensive to collect, and fixed locations
of sensors limit the observation coverage. Note also that the Eulerian observation
function hϕ (·) maps the vehicle density ρ to the vehicle flux ϕ(ρ) and is therefore
typically nonlinear: the nonlinearity of the observation function will influence the
efficacy of data assimilation, especially for the ensemble Kalman filter, and requires
special treatment which will be discussed in §5.4.2.
Lagrangian observation data: Lagrangian observations consist typically of the
positions and velocities of moving vehicles collected by GPS devices or cell phones.
If there are Nc cars equipped with these devices, and their positions and velocities
are denoted by (p, v) := (pj , vj )j=1,2,··· ,Nc , then the complete system of the traffic
model (M), capable of tracking the macroscopic traffic density and the positions of
the Nc individual cars, and the Lagrangian observation operators (D) is given by
ρt (x, t) + ∂x ϕ(ρ(x, t)) = ερxx (x, t),
x∈T
(M)
p˙j (t) = vj (t) = V (ρ(pj (t), t)),
j = 1, 2, · · · , Nc
(5.5)
hp (p, v) = p
(D)
hv (p, v) = v.
Thus, in contrast to previous methods for assimilating Lagrangian data that
relied on its reformulation in Lagrangian coordinates or reshaping the macroscopic
traffic model in terms of velocity or spacing [79, 80], our approach adds the governing
equations for the positions of individual cars to the macroscopic model by exploiting
that our constitutive law relates density and velocity.
76
The main advantage of Lagrangian data is the broad coverage of traffic states as
the tracked vehicles travel along the road. The disadvantage is that density informa-
tion is not collected directly. Furthermore, the observation noise of Lagrangian data
is relatively large (see below for details) and depends, for instance, on the carrier’s
signal and weather conditions.
Eulerian and Lagrangian observation data: We now outline how Eulerian and
Lagrangian data can be assimilated simultaneously by concatenating the systems
(5.4) and (5.5) that we introduced above. The resulting system for the traffic flow
model (M) and the observation operator (D) is given by
ρt (x, t) ερxx (x, t) − ∂x ϕ(ρ(x, t))
(M) p˙ (t) =
j V (ρ(pj (t), t))
j = 1, 2, · · · , Nc
vj (t) V (ρ(pj (t), t)) j = 1, 2, · · · , Nc
(5.6)
ϕ
h (ρ) (ϕ(ρ(qi , t)))i=1,2,··· ,Ns
(D) hp (p, v) =
p .
v
h (p, v) v
The framework proposed in (5.6) is capable of assimilating both Eulerian sensor and
Lagrangian GPS data as it is a closed system for the density ρ and the locations and
velocities of the Nc vehicles from which GPS data are gathered. We note that we
will consider Nc to be fixed but emphasize that our framework allows this quantity
to depend on function of time as long as it changes slowly compared to the time
scale over which GPS data are collected.
Observation errors: We assume that the observation errors for sensor, GPS po-
sition, and GPS velocity data are independent and normally distributed with zero
77
mean. As mentioned above, the observation noise ηiϕ for sensor data has a variance
of 0.1% ϕ(qi , t). Based on the current accuracy of GPS data, we assume that position
and velocity measurements have errors with standard deviation of 5.12m ([35]) and
0.0707m/s ([58, 62]), respectively. Using these assumptions, the covariance matrix
of the observation errors can be written as
0.1% (ϕ(qi , t))i=1,...,Ns 0 0
R(t) =
0 2
(5.12m) INc 0 ,
(5.7)
0 0 (0.0707m/s)2 INc
where IN denotes the N × N identity matrix.
5.3.2 Discretized observation models
We now discuss the numerical implementation of the continuous scheme (5.6), which
is obtained by discretization in space and time. Recall that we consider a ring road
T of length L.
First, we consider the discretization of the traffic model (M) given by
ρt (x, t) ερxx (x, t) − ∂x ϕ(ρ(x, t))
p˙ (t) = V (ρ(p (t), t)) j = 1, 2, · · · , Nc (5.8)
j j
vj (t) V (ρ(pj (t), t)) j = 1, 2, · · · , Nc .
For the spatial discretization, we pick a large integer M , define a spatial step size
∆x := L/M , and choose M consecutive mesh points xm ∈ T that are equally spaced
along the ring road with distance ∆x. Time is discretized by choosing a small positive
time step ∆T and evaluating solutions at times t = n∆T for integers n ∈ N.
78
The traffic state at time t = n∆T consists of the densities at each mesh point
xm ∈ T as well as the GPS positions and velocities of tracked vehicles and is therefore
given by
X(t) := [ρ(x1 , t), . . . , ρ(xM , t), p1 (t), . . . , pNc (t), v1 (t), . . . , vNc (t)]T . (5.9)
| {z } | {z } | {z }
densities GPS positions GPS velocities
The discretized system obtained from solving (5.8) numerically can then be written
as
X(t) = f (X(t − ∆T ); θ(t − ∆T )), (5.10)
where θ(t) denotes system parameters, which may be constant or vary in time. The
function f is obtained using the following numerical schemes. First, we discuss
how we discretize the viscous Lighthill–Whitham equation for the density ρ(x, t)
that appears in (5.8): for periodic boundary conditions, we approximate the spatial
derivatives in the viscous Lighthill–Whitham equation for the density ρ(x, t) using
spectral differentiation matrices in physical space as outlined in [71, (3.10) in §3]; for
non-periodic boundary conditions, we discretize the spatial derivatives using finite
differences. Next, we describe how we evaluate the equations for the positions and
velocities (pj , vj ) that appear in (5.8): to determine the values ρ(pj , t) of the density
at arbitrary positions pj , we interpolate ρ from its known values ρ(xm , t) at the
mesh points xm using trigonometric interpolation ([40]). Finally, we integrate the
resulting system of differential equations from time t − ∆T to time t using a third-
order Runge–Kutta method.
We now move to the observations. Using the observation errors introduced above,
79
Eulerian and Lagrangian observations can be written as
Y(t) = h(X(t)) + η(t), η(t) ∼ N (0, R(t)), (5.11)
where h(X(t)) with X(t) given in (5.9) is defined via
h(X(t)) := (ϕ(ρ(q1 , t)), . . . , ϕ(ρ(qNs , t)), p1 (t), . . . , pNc (t), v1 (t), . . . , vNc (t))T . (5.12)
In order to evaluate ρ(qi , t) at the location qi of a sensor, we can either arrange
that qi coincides with a mesh point xm or again use trigonometric interpolation
([40]). We remark that the expression in (5.12) can be modified easily to account for
the situation where we have observations from sensors only or from GPS data only
available for data assimilation.
5.4 Ensemble methods
The goal of data assimilation is to accurately predict the states and parameters of a
numerical model by incorporating available information obtained from observations
into the model. In our work, this technique is carried out to estimate the traffic
states and uncertain parameters in the traffic flow models described above. A general
introduction to data assimilation and filters is provided in this section.
We first review two common sequential data assimilation methods, the particle
filter (PF) and the ensemble Kalman filter (EnKF). The goal of each of these meth-
ods is to estimate the Bayesian prior and posterior probability distributions using an
ensemble of possible states, as they change in time. First, define the state of interest
to be X, and the observation Y to be a noisy function of the state: Y = h(X) + η,
80
where the observation error is normally distributed η ∼ N (0, R) with N (·, ·) denot-
ing the Gaussian distribution. Then, we are interested in estimating the posterior
probability distribution p(X|Y) given a prior distribution p(X) based on our current
knowledge of the state, and a likelihood distribution p(Y|X) of the observations Y
given the state X, based on our knowledge of the noise in the observations. Bayes’
rule gives the true posterior distribution in this case:
Z
p(Xt |Y1:t ) ∝ p(Yt |Xt ) p(Xt |Xt−1 )p(Xt−1 |Y1:t−1 ) (5.13)
where p(Y|X) is Gaussian likelihood with Y|X ∼ N (h(X), R). In both the PF and
EnKF methods, these distributions are approximated by a weighted ensemble of the
state X given by {Xi , wi }N
i=1 which implies the distribution
e
Ne
X Ne
X
wi δ(X − Xi ) with wi = 1 , (5.14)
i=1 i=1
where δ(X − Xi ) is the Dirac delta centered at Xi . Generally for the EnKF, wi =
1/Ne . Between the observation times, the weights are kept fixed and the state
variables are evolved according to the dynamics of the system. The main difference
between the two methods comes at the time when observations are available. This
is described in the next two subsections. We will use the notation that {Xfi , wif }N e
i=1
is the ensemble from the prior distribution p(X) whereas {Xai , wia }N
i=1 is from the
e
posterior distribution p(X|Y).
5.4.1 Particle filter
The posterior is approximated by updating the weights but leaving the particle
positions fixed, i.e., Xai = Xfi := Xi . The updated weights are obtained by applying
81
Bayes’ rule to the weights as follows:
p(Y|Xi )wif
wia = PNe f
, (5.15)
j=1 p(Y|Xj )wj
that is, the updated weights are found by multiplying the likelihood of that particle
by the previous weight and normalizing to sum to 1. This is the simplest imple-
mentation of the PF, also known as sequential importance sampling ([27, 34]). For
a more complete derivation of the particle filter and the proposal distributions, see
[26, 66, 68].
Due to the finite nature of the approximation and the recursive updating of the
weights, sequentially applying this algorithm eventually leads to one particle with
very high weight, while the rest of the particles have almost zero weight (so-called
“filter divergence” or “weight collapse”.) To prevent the weights of the ensemble con-
centrating on a single state, various resampling methods for the ensemble states may
be used ([74]). The basic idea behind each of these methods is to monitor when a pre-
determined threshold is hit (for example, when the “effective sample size” becomes
small ([51])), and to then resample the particles from the discrete approximation
of the posterior distribution and reset all weights to 1/Ne . We approximate the
effective sample size to be
1 N
e for wj = 1/Ne
Neff ≈ PNe = (5.16)
i=1 wi2
1
for wj = δij
thresh
as in [51]. We will apply resampling when Neff < Neff for a predetermined thresh-
thresh
old Neff that will typically be a small fraction of the total number of particles
Ne .
82
One major drawback of the PF is that it has been shown to fail in high dimen-
sions ([67]). Thus, in the case where the state of interest is a spatially-discretized
function over some domain, the PF is intractable. However, when the state dimen-
sion is small enough and the number of particles is large enough, the PF can provide
an accurate approximation to the exact Bayesian posterior distribution. This is es-
pecially useful if this distribution is skewed or multimodal, which is often the case
when the dynamical system governing the state is nonlinear.
5.4.2 Ensemble Kalman filter
Like the PF, the EnKF [28, 29] employs an ensemble of state vectors {Xi }i=1...Ne to
represent the posterior distribution; however, unlike the PF, the ensemble members
are equally weighted for the entire assimilation window. Instead of updating the
weights at analysis times, the members themselves are updated according to an
ensemble approximation of the traditional Kalman filter update step, given here by
the so-called perturbed observation EnKF [15, 29, 43]
Xai = Xfi + K(Y − HXfi + ηi ), (5.17)
−1
K = Pf HT HPf HT + R , (5.18)
where Xfi is the forecast of the ith ensemble member, Xai is the ith updated (analysis)
ensemble member, H is the observation operator, K is the Kalman gain matrix, R
is the covariance of the observation error, and ηi ∼ N (0, R) are the observation
perturbations. Pf is the forecast ensemble covariance given by
eN Ne
1 X 1 X
f
P = (Xf − X
¯ f )(Xfi − X
¯ f )T , ¯ f
X = Xf . (5.19)
Ne − 1 i=1 i Ne i=1 i
83
Nonlinear observation function: We note that the update step requires that the
observation operator H is linear. However, the observation function given in (§5.12
equation 5.12) is nonlinear, and we will therefore outline how we adjust the update
step to accommodate the nonlinear observation function. Instead of calculating the
operator H and the covariance matrix Pf separately when estimating the state Xai
in (5.17) and (5.18), we follow [56] and treat Pf HT and HPf HT as follows as single
entities:
e N Ne
1 X 1 X
HP H = f T
(Zfi − Z
¯ f )(Zfi − Z
¯ f )T , ¯ f
Z := Zfi , (5.20)
Ne − 1 i=1 Ne i=1
e N Ne
1 X 1 X
P H =f T
(Xf − X
¯ f )(Zfi − Z
¯ f )T , ¯ f
X := Xf , (5.21)
Ne − 1 i=1 i Ne i=1 i
where Zfi := h(Xfi ) denotes the ensembles projected into observation space. Simi-
larly, the term HXfi in (5.17) is replaced by Zfi .
Localization: There has been a lot of work towards improving the EnKF to make
it feasible for very high dimensional systems; in particular, covariance inflation and
localization have provided significant improvement of the performance of the EnKF
(see [4, 38, 43, 44]). Since the ensemble size is often much smaller than the dimension
of the state, spurious correlation will arise in the sample covariance matrices. Local-
ization is a method in which these spatially long-range correlations are diminished:
in particular, this is often implemented via a Schur product of the sample covariance
matrix and some sort of cutoff matrix.
One common localization method applies an exponential decay function, such
as the Gaspari–Cohn function in [31], to the forecast covariance matrix Pf . In our
formulation of the Kalman gain matrix, the forecast covariance matrix Pf is not
84
directly accessible as we provide only HPf HT and Pf HT through (5.20) and (5.21),
respectively. As in [44], we could now apply elementwise localization separately to
each of these two matrices: localization of Pf HT would account for the assumption
that observations at a given location would only affect the density and car positions
near that location, while localization of HPf HT would correspond to localizing the
effect observations at different positions have on each other.
Alternatively, if we interpret the operator K = (Kij ) loosely as a weight matrix
where the (i, j)th entry Kij determines how much the jth entry of the observation
vector affects the ith entry of the state vector, then we can apply localization directly
to the matrix K to ensure that observations at a certain location affect only nearby
states. Specifically, we assume that only mesh points within 0.5 miles of the location
of a sensor are affected by Eulerian observations. We then define an exponential
decay localization function for the jth column of the Kalman gain matrix K via
e−d|xm −(qj +s)| with |xm − qj | < 0.5 miles, (5.22)
where xm is the coordinate of the mth mesh point, qj is the coordinate of the jth
sensor, and d is a decay coefficient that describes the decay rate of the localization
function. Since cars detected by a sensor will have travelled a certain distance s
during time step ∆T , we center the localization at qj + s. Throughout, we use the
parameters d = 0.5 and s = 0.35 miles. We use the same localization function also
for Lagrangian observations: we set d = 1.2 and choose s = 0 as there will be no
time delay for positions and velocities collected from GPS data.
Figure 5.2 shows that applying localization directly to the Kalman gain matrix
turns out to be very similar to applying localization to the forecast covariance matrix.
Indeed, our numerical results show that the observation covariance HPf HT + R is
85
1 1
0.8 0.8
50 0.6 50 0.6
0.4 0.4
100 0.2 100 0.2
0 0
150 −0.2 150 −0.2
−0.4 −0.4
200 200
−0.6 −0.6
−0.8 −0.8
250 250
−1 −1
5 10 15 5 10 15
Figure 5.2: Comparison of the scaled Kalman gain matrices with localization applied to
the forecast covariance matrix (left) and directly to the Kalman gain matrix (right).
dominated by its diagonal elements: using the notation I for the identity matrix,
we then see that the Kalman gain matrix, with localization applied to the forecast
covariance matrix, can be approximated by the Kalman gain matrix with direct
localization:
−1 −1
loc ◦ (Pf HT ) I ◦ (HPf HT ) + R = loc ◦ Pf HT I ◦ (HPf HT ) + R
−1
≈ loc ◦ Pf HT HPf HT + R .
Inflation: The EnKF is not immune to filter divergence, where the prior ensemble
spans a smaller and smaller space until the observations no longer have any effect
in the analysis step ([4]). The EnKF is prone to artificial collapse of the covariance
matrix. To address this issue, the covariance may be “inflated” before the Kalman
update step is performed. We use two different inflation algorithms for the EnKF.
First, we employ a constant scalar inflation for the simulations of multiple traffic
86
scenarios as the localization algorithm described above does not significantly tighten
the covariance matrix. Second, for the computations with realistic traffic data, we
use adaptive inflation as described in [2, 3].
5.4.3 Parameter estimation
The macroscopic model we utilize may contain parameters whose values are not
known: for instance, the expression (5.2) for the Greenshields velocity function con-
tains the parameters vmax and ρmax , which will generally be unknown.
We briefly describe three approaches for how data assimilation may be used
to estimate parameters in addition to the state of interest. One approach is the
augmented-vector method, in which the parameters to be estimated are appended
to the state vector ([6, 55]) and are updated concurrently. Another approach is the
“separated” method [9], in which the state is updated first, followed by updating the
parameters, using localization on the state and non-global parameters. The third
approach applies a two-stage filter [65] that consists of applying the EnKF to the
state estimate and the PF to the parameters.
There are then several methods that may be used to model the evolution of the
parameters between assimilation steps. The parameters may be kept fixed between
time steps (persistence model, see [23]), but this can lead to degeneracy (sample
attrition) when using PF methods. In [55], the authors suggest an artificial evolution
of parameters between observations, in which the parameter evolution model is given
87
by
θt+1 = θt + ξt+1 , (5.23)
ξt+1 ∼ N (0, Wt+1 ), (5.24)
where Wt+1 is a specified covariance matrix. In [78], multiplicative parameters were
estimated by evolving the parameters under a temporally smoothed version of the
persistence model. It was shown that this results in smooth temporal variation of
the parameters, which in turn prevents blow-up of the model.
The augmented-vector approach is used to estimate parameters, and artificial
evolution is introduced to time-varying parameters.
5.5 Numerical results
5.5.1 Performance criteria
In this paper, we use both absolute and relative root-mean-square-errors (RMSE)
as performance indicators to assess different scenarios. The absolute performance
index at time tn is defined as
s
PM
m=1 (ρm (tn ) − ρˆm (tn ))2
RMSEA = , (5.25)
M
while the corresponding relative performance index is given by
RMSEA
RMSER = , (5.26)
ρmax
88
where ρmax is the maximal density in the macroscopic model. Here, M is the number
of mesh points used to discretize a given highway stretch, ρm (tn ) denotes the true
density at the mth mesh point xm at time tn , and ρˆm (tn ) denotes the corresponding
estimation from data assimilation. For those traffic scenarios with unknown maximal
density ρmax , the denominator in the relative RMSE is replaced by the average
PM
m=1 ρm (tn )/M of the true density.
5.5.2 Traffic scenarios with true underlying model
The main goal of this part is to apply data assimilation to estimate traffic states
and parameters using the theory described in § 5.3 and § 5.4. Specifically, we will
investigate
(i) assimilation of Eulerian and Lagrangian data in dynamic traffic flow;
(ii) efficacy of Lagrangian data assimilation;
(iii) impact of sensor location on data assimilation;
(iv) parameter estimation of traffic flow models.
We note that the first three topics are investigated under the assumption that the
estimator has the exact values of model parameters, and the emphasis is placed on
the capability of traffic state estimators. However, in the last topic, the parameter
estimation is activated with the intent of investigating the capability of parameter
estimators.
89
Basic setup: In the simulations outlined in § 5.5.2, the underlying model for the
data assimilation is the same as the model generating the truth; in each case, the
viscous Lighthill-Whitham equation (5.1) is used with additional terms to simulate
various traffic scenarios, to be outlined individually below. We will show that data
assimilation can be used to reduce the error in simulating these scenarios when
imperfect initial data is used. The initial data for density is set as the truth plus
10% noise through Fourier coefficients of the density curve.
For each traffic scenario and associated macroscopic system, we consider traffic
flow on a ring road of length L = 50 miles, which is equally divided into 256 mesh
points for the numerical computation. The truth is generated by initializing the
system with the profile
ρ0 (x) = 0.5ρmax + 0.4ρmax sech(x − L/2), (5.27)
and the parameters
ρmax = 45 cars/mile, vmax = 75 miles/hour, ε = 0.1. (5.28)
For the data assimilation algorithm, the underlying model is taken to be the
same as that generating the synthetic truth; however, the initial conditions for the
ensemble are the noisy density curves after adding noise directly through the Fourier
coefficient. Aside from the fourth investigation, where parameter estimation is em-
ployed, we use the same systems parameters in both the generation of the truth and
the data assimilation.
The Eulerian and/or Lagrangian observations are collected and assimilated every
one minute, and the entire simulation has 180 total updates. For the Eulerian
90
observation system, we place Ns = 8 stationary sensors equidistantly along the
ring road. For the Lagrangian observation system, we pick Nc = 15 cars with GPS
devices that are initially equally spaced along the ring road. The ensemble size is 30
for the EnKF and 300 for the PF because of the high dimension of the state vector.
In each of the following investigation, we will outline the traffic scenarios to
be studied and the associated macroscopic models, followed by the results of data
assimilation in estimating the traffic states which emerge.
1. Data assimilation of Eulerian and Lagrangian data First, we demon-
strate that the proposed algorithm can be used to assimilate Eulerian, Lagrangian,
and combined Eulerian-Lagrangian observations into a macroscopic model to accu-
rately estimate traffic states. To test the algorithm, we consider a ring road with
a traffic light for our simulations and assimilate sensor and GPS data. In order to
better evaluate the performance of the algorithm, we will compare the results also
to those obtained by simulating the macroscopic model with the same initial data
and parameter values but without assimilating observations.
The traffic light scenario is described in § 2.3.3 by the following equation
ρt + ρVG (ρ(x, t)) a(x, t; x`1 ) x = ερxx ,
(x, t) ∈ (x`1 − Ly1 , x`1 ) × T1y
0.5
(x, t) ∈ (x`1 − Lr1 , x`1 ) × T1r
0
`
a(x, t; x1 ) =
(x`1 − Lr1 − x)/Lr1 (x, t) ∈ (x`1 − 2Lr1 , x`1 − Lr1 ) × T1r
1
otherwise,
The traffic light is placed in the middle of the road with x`1 = 25. We set the lights
91
length (|T1y |, |T1r |, |T1y |) = (10, 190, 400) seconds and response distances (Ly1 , Lr1 ) =
(1, 0.8) miles in order to generate significantly oscillating traffic flow.
We use the above model and the initial conditions and parameters from the basic
setup to generate the truth. The same system and parameters are used for the
underlying data assimilation model, with perturbed initial conditions.
Results: We carried out simulations using the algorithm described in § 5.3 with
the EnKF and PF for (i) Eulerian observations only, (ii) Lagrangian observations
only, (iii) both Eulerian and Lagrangian observations, and (iv) no observations. The
simulation results based on the traffic lights scenario are shown in Figure 5.3: panels
(a) and (b) contain the results for the EnKF and PF, respectively. The relative
RMSE for the cases (i)-(iii) ranges between 1% and 2% after 3 hours, while the
relative RMSE for case (iv) (no observations) is above 7%. The results suggest that
our algorithm is capable of assimilating any combination of Eulerian and Lagrangian
data and performs well independently of whether we use the EnKF or PF.
(a) (b)
0.08 0.08
0.07 0.07
Relative RMSE
Relative RMSE
0.06 0.06
No DA No DA
0.05 EnKF Sensors 0.05 PF Sensors
EnKF GPS PF GPS
0.04 EnKF Sensors&GPS 0.04 PF Sensors&GPS
0.03 0.03
0.02 0.02
0.01 0.01
0 ! 0 !
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time(hours) Time(hours)
Figure 5.3: Traffic state estimation in the traffic light scenario. Shown is the relative RMSE for
assimilating no data (yellow dashed dot), Eulerian data (blue dashed), Lagrangian data (magenta
dotted), and both Eulerian and Lagrangian data (green solid) using the (a) EnKF and (b) PF.
92
2. Efficacy of Lagrangian data assimilation In practice, either position or
velocity data can be collected through GPS devices. Thus, we are interested in
the efficacy and accuracy of different Lagrangian observations. According to the
statistical data of GPS devices on the internet, the error of position data has a
standard deviation around 5.12 meters [35], while the error of velocity data has a
standard deviation around 0.0707 m/s [58, 62]. In this section, both the EnKF
and PF are used to assimilate three groups of different Lagrangian observations in
a normal traffic scenarios: (1) positions of vehicles, (2) velocities of vehicles, (3)
combined positions and velocities of vehicles.
we simulate a scenario of unimpeded traffic flow, which we model via the viscous
Lighthill-Whitham system (2.10)
ρt + (ρVG (ρ))x = ερxx
with no additional terms added. The initial conditions and parameters are as out-
lined in basic setup.
Results: The simulation results are shown in Figure 5.4. Figure 5.4(a) shows the
relative RMSE of Lagrangian data assimilation for different observation data using
the EnKF. The error for assimilating position data (solid blue) approaches 3.5%, the
error for assimilating velocity data (dashed blue) approaches 0.4% and the error for
assimilating both data (dotted blue) approaches 0.1%. In contrast, Figure 5.4(b)
shows the comparison based on the PF: the error for assimilating velocity data
(dashed magenta) approaches 0.25%, and the errors for assimilating position data
and both data (solid and dotted magenta lines) approach 0.4%. The phenomenon
that assimilating both is less accurate initially is related to the “curse of dimension-
93
ality” of the PF: the dimension of observations is doubled when assimilating both
data, but we are unable to correspondingly increase the ensemble size to compensate.
0.045 0.035
(a) (b)
0.04
0.03
0.035
0.025
Relative RMSE
0.03
Relative RMSE
0.025 EnKF Tracking Position 0.02 PF Tracking Position
EnKF Tracking Velocity PF Tracking Velocity
EnKF Tracking Both PF Tracking Both
0.02
0.015
0.015
0.01
0.01
0.005 0.005
0 ! 0 !
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time(hours) Time(hours)
Figure 5.4: Traffic state estimation comparison for Lagrangian data assimilation based on different
observation data from normal traffic flow modeled by (2.10). Shown is the relative RMSE for the
(a) EnKF and (b) PF, where we assimilate position data only (solid), velocity data only (dotted),
and combined position and velocity data (dashed).
For both the EnKF and PF, we can see that assimilating velocity data is more
accurate compared to assimilating position data, which is consistent with the fact
that velocity data has smaller observation noise than position data. However, the PF
has a stronger tolerance of observation noise since the difference of error between as-
similating velocity data and position data is much smaller compared to those by the
EnKF. In addition, the performance of the EnKF can be improved with combined
observation data, while the performance of the PF is limited by the number of par-
ticles with combined observation data. Therefore, we will use combined Lagrangian
data for EnKF, while only use velocity data for PF in other simulations.
3. Impact of sensor location on data assimilation When the Department of
Transportation installs sensors on the road, an important consideration is where to
install them in order to optimize the amount of traffic information collected by them.
In this section, we are interested in studying the impact of the configuration of sensors
94
in the presence of on-/off-ramps or bottlenecks in order to efficiently place sensors
near ramps and bottlenecks. We will use two traffic scenarios for investigation: an
on-/off-ramps scenario and a stationary bottleneck scenario.
The on-/off-ramps scenario is described in § 2.3.2 with the following formula:
X X
ρt + (ρVG (ρ))x = ερxx + ϕon on
i (t)f (x − xi ) + ϕoff off
j (t)f (x − xj ),.
i∈I j∈J
we place one on-ramp at xon off
1 = 37.5 miles and one off-ramp at x1 = 12.5 miles on
the road, and set constant flow ϕon off
i (t) = 2ρmax for the on-ramp and ϕi (t) = 4ρmax
for the off-ramp. Three comparative experiments are designed where a sensor is
installed (1) at the ramp, (2) 0.5 miles upstream to the ramp, and (3) 0.5 miles
downstream to the ramp.
The stationary bottleneck scenario is describe in § 2.3.4 with the following for-
mula:
ρt + ρVG (ρ)a(x − xb1 ) x = ερxx .
We set the bottleneck at xb1 = 25 miles, severity coefficient c = 1 and spread effect
f (x) = 1 − 0.5sech(x). In this scenario, a stationary density bottleneck appears
around xb1 . Three comparative experiments are designed where a sensor is installed
(1) at the bottleneck, (2) one mile upstream to the bottleneck, and (3) one mile
downstream to the bottleneck.
Results: The simulation results for ramps are shown separately in Figure 5.5(a)
and 5.5(b), and those for the bottleneck are shown in Figure 5.5(c). In Figure 5.5(a),
it is observed that the performance of data assimilation is improved by moving a
95
0.06
(a) On-ramp Upstream (b) O↵-ramp Upstream (c) Bottleneck Upstream
0.045 At on-ramp At o↵-ramp At bottleneck
Downstream 0.05 Downstream 0.05 Downstream
0.04
0.035
Relative RMSE
0.04
Relative RMSE
Relative RMSE
0.04
0.03
0.025 0.03 0.03
0.02
0.02 0.02
0.015
0.01
0.01 0.01
0.005
0 ! 0 ! 0 !
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time(hours) Time(hours) Time(hours)
Figure 5.5: Evaluation of various sensors configurations using the EnKF. Results for different
sensor locations near (a) on-ramps, (b) off-ramps, and (c) bottlenecks are shown, where the solid
curves correspond to sensor locations upstream to the target location, the dotted curves represent
sensors at the target location, and the dashed curves represent locations downstream to the target
location.
sensor from downstream to upstream of the on-ramp. While in Figure 5.5(b), unlike
the on-ramp case, the performance is improved by moving a sensor from upstream
to downstream of the off-ramp. In Figure 5.5(c), it shows that the performance
is improved by moving a sensor from upstream to downstream of the bottleneck.
By studying the shape of density on the ring road, we find that a density bump is
generated gradually upstream to the on-ramp. Therefore, installing a sensor slightly
upstream can provide more traffic information and thereby improve the performance
of the EnKF. In contrast, a density valley is generated gradually downstream to the
off-ramp or the bottleneck, so installing a sensor slightly downstream can provide
more traffic information.
In summary, the locations of sensors do impact the performance of data assimila-
tion when using the EnKF: the simulation results suggest installing sensors upstream
of on-ramps and downstream of off-ramps and bottlenecks. In contrast, the PF per-
forms well regardless of the locations of sensors (simulation results not shown).
4. Parameter estimation of traffic-flow models In previous parts, we showed
that our algorithm successfully assimilates traffic states when the parameters in the
96
underlying traffic model are known. However, in practical applications, these pa-
rameters are not known and need to be estimated when traffic data are assimilated.
Hence, we are interested in the efficacy of data assimilation when traffic parameter
estimation is activated. In this part, we consider two traffic scenarios for investiga-
tion: an on-/off-ramps scenario and traveling bottleneck scenario. For each scenario,
the basic setup for generating the truth and the underlying model for data assimi-
lation is the same as in the setup; however, the parameters in the underlying traffic
model for data assimilation are now taken as unknown.
The model for the on-/off-ramps scenario is as described in the 3rd case for on/off-
ramps, and the unknown parameters include maximal density ρmax , maximal velocity
vmax , the fluxes of on-ramp ρon off
1 and the fluxes of off-ramp ρ1 . These parameters are
taken to be constant throughout the simulation.
The traveling bottleneck scenario is described in § 2.3.5 with the formula:
ρt + ρVG (ρ)a(x − xb1 (t)) x = ερxx .
The position of the bottleneck is assumed to move back and forth periodically ac-
cording to the formula xb1 (t) = L/2 + L/4 cos(2πt/3). The unknown parameters
include maximal density ρmax , maximal velocity vmax , taken to be constant, and the
bottleneck location xb1 (t), which is time-dependent.
Results: As shown in Figure 5.6(a), the relative RMSE by using the EnKF (blue
dashed line) and the PF (magenta dashed line) hovers around 3% and 2%, respec-
tively, after 3 hours of assimilating Eulerian observations. In Figure 5.6(b), the
relative RMSE by using the EnKF (blue dash line) and the PF (magenta dash line)
is around 10% and 20%, respectively, after 3 hours of assimilating Lagrangian ob-
97
servations. Observations of the flux, which is the product of density and velocity,
can balance the estimation of ρmax and vmax and thereby provide accurate estimation
with an error below 5%. However, Lagrangian observations, especially of velocity
data, cannot provide as good a balance between density and velocity as Eulerian ob-
servations, and we indeed find that Lagrangian observations are less efficient when
the parameters are unknown. We also see in Figure 5.6(c) and 5.6(d) that the PF
estimator (magenta dashed) is more sensitive when considering the estimated fluxes
as these oscillate around the true flux value near the ramps.
0.08 0.25
(a) EnKF Sensor PEOFF (b) EnKF Tracking PEOFF
EnKF Sensor PEON EnKF Tracking PEON
PF Sensor PEOFF PF Tracking PEOFF
PF Sensor PEON 0.2 PF Tracking PEON
0.06
Relative RMSE
Relative RMSE
0.15
0.04
0.1
0.02
0.05
0 ! 0 !
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Times(hours) Times(hours)
50 85 60 85
(c1) Est of ⇢max (c2) Est of vmax (d1) Est of ⇢max (d2) Est of vmax
Velocity(miles/hour)
True
Velocity(miles/hour)
Density(cars/mile)
Density(cars/mile)
True EnKF
48 EnKF
80
55 PF
80
PF
46 75
75 50
44 70
!
42 ! 70 ! 45 65
2.8 5 2.4 4.5
(c3) Est of ⇢on (c4) Est of ⇢o↵ (d3) Est of ⇢on (d4) Est of ⇢o↵
Flux(cars/hour·mile)
Flux(cars/hour·mile)
Flux(cars/hour·mile)
1 1 1 1
Flux(cars/hour·mile)
2.6 2.2
4.5
2.4 2 4
4
2.2 1.8
2! 3.5 ! 1.6 ! 3.5 !
0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3
Figure 5.6: Estimation of traffic states and parameters in the ramps scenario: shown are the
relative RMSE for assimilating (a) Eulerian and (b) Lagrangian observations as well as the results
of parameter estimation for (c) Eulerian and (d) Lagrangian observations. Estimated are the
maximal density ρmax , the maximal velocity vmax , the flux of the on-ramp ρon
1 , and the flux of the
off-ramp ρoff
1 . PEOFF represents that the parameters are known, while PEON represents that the
parameter estimation is considered.
In the bottleneck scenario, the conclusions are similar as those in the ramps
scenario. The relative RMSEs are below 5% in Figure 5.7(a) and below 9% in
98
Figure 5.7(b), which means our approach has provided a very good state estimator.
In addition, Figure 5.7(c) and 5.7(d) demonstrate the dynamical tracking capability
of our approach for the time-varying parameter xb1 (t).
0.08 0.1
(a) (b)
0.07 0.09
0.08
0.06
0.07
Relative RMSE
Relative RMSE
0.05
0.06
EnKF Sensor PEOFF EnKF Tracking PEOFF
0.04 EnKF Sensor PEON 0.05 EnKF Tracking PEON
PF Sensor PEOFF PF Tracking PEOFF
PF Sensor PEON 0.04
0.03 PF Tracking PEON
0.03
0.02
0.02
0.01
0.01
0 ! 0 !
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Times(hours) Times(hours)
50 85 49 85
(c1) Est of ⇢max True (c2) Est of vmax (d1) Est of ⇢max True (d2) Est of vmax
Velocity(miles/hour)
Velocity(miles/hour)
Density(cars/mile)
Density(cars/mile)
EnKF EnKF
48 80 48
PF PF
80
46 75 47
75
44 70 46
42 65 45 70
40 40
Location(miles)
Location(miles)
30 30
20 20
(c3) Est of xb1 (d3) Est of xb1
10 ! 10 !
0 1 2 3 0 1 2 3
Figure 5.7: Estimation of traffic states and parameters in the traveling bottleneck traffic scenario:
shown are the relative RMSE for assimilating (a) Eulerian and (b) Lagrangian observations and
the results of parameter estimation for (c) Eulerian and (d) Lagrangian observations, where we
estimate the maximal density ρmax , the maximal velocity vmax , and the time-dependent bottleneck
location xb1 .
In summary, our approach is capable of providing accurate estimations for both
constant and time-varying parameters. Generally Lagrangian observations are more
efficient when estimating traffic states, but Eulerian observations appear to be slightly
more efficient when parameter estimation is activated.
99
5.5.3 Traffic scenarios with unknown underlying model
The simulations in §5.5.2 used the same traffic model for (i) generating the truth
and (ii) assimilating data in our algorithm. The main goal of this section is to test
the implementation of data assimilation to estimate traffic states when the truth is
not generated by the underlying model for the data assimilation scheme.
More specifically, in case 1, we will use data assimilation to estimate traffic states
using observations generated from a microscopic traffic flow model. Then, in case 2,
we will apply parameter estimation to real traffic data obtained from the Minnesota
Department of Transportation.
5. Truth generated from microscopic models We are interested in testing
the sensitivity with respect to the model used in the data assimilation algorithm in a
controlled environment. In this section, we generate synthetic observations using an
extended optimal-velocity type microscopic model and assimilate the data using the
macroscopic LWR equation under two different scenarios of highway traffic. These
scenarios, to be described in more detail below, are obtained by different choices of
the parameters and initial conditions in the microscopic model.
We begin with an overview of the microscopic model, followed by the results
when using data assimilation to estimate traffic states in each of the scenarios. In
our simulations below, we will study each scenario separately using the EnKF to
assimilate Lagrangian observations (vehicle positions) from the microscopic model
to estimate the macroscopic traffic density. To compare the predicted traffic densities
with the truth, we transform the microscopic traffic data into continuous densities
on the ring road.
100
The microscopic optimal-velocity model attempt to model traffic through differ-
ential equations governing the evolution of the individual cars’ positions and veloci-
ties. In general, these models allow for each car to adjust its velocity and accelera-
tion based on the position and/or velocity of its neighbors. We focus on an extended
optimal-velocity type model introduced in [16] in which the drivers adjust their ve-
locity according to not only the headway but also the relative velocity to the car in
front. In this model, the position pn (t) and the individual target headway sn (t) of
the nth vehicle evolve according to
pn+1 − pn − sn
τ p¨n = Vg tanh + V0 − p˙n (5.29)
`0
Vg
αs˙ n = s − sn ) − β (p˙n+1 − p˙n ) ,
(¯
`0
where the dot symbol denotes d/dt. Here τ is the reaction time, and the right hand
side of the equation is the so-called optimal velocity function. The adjustment V
depends on the difference of the headway pn+1 − pn from the target headway sn
that encompasses the car length plus a safety distance. A common choice for this
function is V(u) = tanh(u), which we will use in this paper. Since V(0) = 0, the
quantity V0 represents an optimal velocity at which the car drives when the headway
is equal to the optimal headway sn . The quantity Vg represents a velocity gain, and
by considering an infinite headway, we have V(∞) = 1, so that the quantity Vg + V0
gives an effective speed limit for the drivers. The parameter `0 is a characteristic
length scale which describes the pace at which the driver attempts to achieve the
optimal velocity.
The second equation describes the evolution of the individual target headway
sn (t): the term s¯ − sn describes the relaxation of sn (t) to an optimal headway s¯,
while the term β(p˙n+1 − p˙n ) takes into account that drivers may increase speed
101
if the car in front does so. The parameters α and β are dimensionless: α is a
measure of the overall adjustment time of the individual headways, while β measures
how proactively drivers react to changes of the relative velocity. Note that setting
(α, β) = 0 recovers the standard optimal-velocity model [8].
We consider the situation of N cars driving on a circular road of fixed length L
under the generalized OV model described above which is achieved by adding the
periodicity conditions
pn+N = pn + L, sn+N = sn . (5.30)
It was shown in [16] that the model (5.29) admits both free-flow and traffic
jam (traveling wave) solutions. In this paper, we take the following choice of the
parameters Vg , V0 , `0 , s¯, τ which were obtained from Japanese highway data ([7, 30]):
Vg = 37.57 mi/hr, V0 = 34.30 mi/hr, `0 = 38.15 ft, τ = 0.5 s, s¯ = 82.00 ft ,
(5.31)
and we consider two different scenarios of highway traffic, to be described below,
obtained by different choices of the parameters α, β, N, L and initial conditions. In
the first scenario, a square wave evolves into a shock with gradually decreasing
amplitude (a solution that can also qualitatively be found in the macroscopic LWR
model), while the second scenario describes the emergence of a traveling wave with
fixed amplitude (which cannot be generated by the LWR model).
Results (square wave scenario): The first scenario is that of a long circular
road consisting of a large number of cars engaged in free flow interacting with a
102
small traffic jam region. We set (α, β) = (1, 4) and the length L = 22 miles. To
initialize the positions of the cars, we place cars equally spaced along the road with
a fixed density of 40 cars per mile except for a two mile stretch in which the density
is increased to 130 cars per mile for a total of N = 1060 cars. This gives that the
initial headways are sf0 = 131.95 ft in the free flow region and sj0 = 40.61 ft in the
traffic jam region. The initial velocities of the cars are taken to be
! !
sf0 − s¯ sj0 − s¯
v0f = Vg tanh + V0 , v0j = Vg tanh + V0 , (5.32)
`0 `0
for the free flow and jam regions, respectively. We initialize the system with the
above conditions and compute the local density at each time step using a circular
kernel density estimator in Matlab. In terms of density, the initial condition resem-
bles a square wave which evolves into a shock with decreasing amplitude over time
(Figure 5.8).
The simulation results for the first scenario are shown in Figure 5.8. We observe
that the estimated density (green circled line) gradually converges to the true value
(blue starred line). This is plausible as the viscous Lighthill–Whitham model can
generate a traveling shock with decreasing amplitude as demonstrated in §5.2. With
the aid of data assimilation, the macroscopic model provides an accurate estimation
for the traveling shock wave transferred from microscopic traffic data.
6. Truth generated from microscopic models
103
130 75
(a) t = 1 min True density (b) t = 20 mins True density
Est density Est density
110
65
Density(cars/mile)
Density(cars/mile)
90
55
70
45
50
30 ! 35 !
0 2 4 6 8 10 12 14 16 18 20 22 0 2 4 6 8 10 12 14 16 18 20 22
Position(miles) Position(miles)
(c) t = 40 mins True density (d) t = 60 mins True density
52 Est density 50 Est density
Density(cars/mile)
50 49
Density(cars/mile)
48 48
46 47
44 ! 46 !
0 2 4 6 8 10 12 14 16 18 20 22 0 2 4 6 8 10 12 14 16 18 20 22
Position(miles) Position(miles)
Figure 5.8: Estimating microscopic data in scenario 1 using the EnKF. True and estimated density
values are represented by blue pluses and green circles, respectively.
Results (traveling wave scenario): The second scenario describes the propa-
gation of a traffic jam or traveling wave solution. We set (α, β) = (4, 0.01) and the
length L = 2.5 miles. We place N = 150 cars offset from an equal spacing of length
(n − 1)L/N on the road with initial positions and velocities
(n − 1)L 2π(n − 1) s0 − s¯
pn (0) = + 10 sin , vn (0) = Vg tanh + V0 , (5.33)
N N `0
where s0 = 87.97 ft is the initial average headway. Figure 5.9 shows the evolution of
the local density of this system. After a short time, a traffic jam, or traveling wave,
emerges that travels backwards against the flow of traffic.
104
Also shown in Figure 5.9 are the results of the simulation using data assimilation.
We observe that the estimated density (green circles) is smoothed out rapidly and
does not converge to the true density value (blue pluses). Unlike the previous case,
the macroscopic model cannot provide an accurate estimation even though data
assimilation is used. This is not unexpected as the viscous Lighthill–Whitham model
cannot generate a traveling square wave shape with stable amplitude and is therefore
not capable of reproducing the correct density profile.
100 110
(a) t = 1 min True density (b) t = 5 mins True density
Est density Est density
90
Density(cars/mile)
90
Density(cars/mile)
80
70
70
60
50 50
40 ! !
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Position(miles) Position(miles)
130 130
(c) t = 10 mins True density (d) t = 15 mins True density
Est density Est density
110 110
Density(cars/mile)
Density(cars/mile)
90 90
70 70
50 50
! !
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Position(miles) Position(miles)
Figure 5.9: Estimating microscopic data in scenario 2 using the EnKF. True and estimated density
values are represented by blue pluses and green circles, respectively.
In summary, a macroscopic model may produce inaccurate results if the model
is not capable of generating dynamical features present in the traffic flow data that
need to be assimilated.
2. Parameter estimation for real traffic data In addition to estimating sim-
ulated traffic states and parameters, we are interested in applying the developed
105
approach to real traffic data. In this section, we use data from the Minnesota De-
partment of Transportation [70].
A freeway strip of I-35E with length 7.175 miles between Main Street and Hwy
96E is considered. There are 14 consecutive detectors S1535-S1548 along this stretch,
as well as three on-ramps and one off-ramp (see Figure 5.10). We picked two repre-
sentative time periods from a day for simulation: one is late night (00:01am-00:30am)
when free flow is expected, and the other is rush hour (05:31pm-06:00pm) when traf-
fic jams are expected. The Minnesota Department of Transportation provides traffic
density, velocity and flow data collected from the sensors, and the density data of
ramps. All traffic data are collected each minute.
A Strip of Highway I-35E
!
0 1 2 3 4 5 6 7
Sensors On-ramps O↵-ramps
Figure 5.10: A strip of highway I-35E in Minnesota. Cars are moving from left to right.
Sensors, on-ramps, and off-ramps are labelled using blue stars, red pluses, and green circles,
respectively.
We use a macroscopic model that takes the on- and off-ramps into account,
and assimilate the flux data collected from the sensors using the EnKF. We also
estimate the unknown maximal density ρmax and the maximal velocity vmax . The
density profiles of the ramps are artificially taken as unknown and are replaced by
the average density value during a period. In contrast to the simulation presented
in previous sections, we use an adaptive inflation algorithm [2] here.
Results: The simulation results for late night traffic are shown in Figure 5.11.
The relative RMSE by assimilating all observation approaches 15%, while the rela-
106
tive RMSE by not assimilating observations approaches 40%. The parameter esti-
mation in 5.11(b) shows that estimated maximal velocity vˆmax is approximately 85
miles/hour, and the maximal density ρˆmax decreases to below 30 cars/mile. This is
a reasonable description of late night traffic.
The simulation results for rush hour traffic are shown in Figure 5.12. The ap-
proach results in a decrease of error from 15% to 8%. Since the flow of ramps is quite
stable during rush hour, the average density value for ramps are close to the true
density value which explains why our approach results in a less drastic improvement
when compared to the previous case. In Figure 5.12(b), the estimated maximal
velocity vˆmax is still approximately 85 miles/hour, but the maximal density ρˆmax
increases to 150 cars/mile, which is characteristic of traffic flow during rush hour.
In summary, the efficacy of our approach has been validated by the real traffic
data from Minnesota Department of Transportation. Estimation for traffic states is
significantly improved in both late night and rush hour conditions, and the estimation
for parameter is consistent with characteristics of traffic scenarios.
Velocity(miles/hour)
0.7 100
(a) Full Observations
(b1)Estimated maximal velocity
Half Observations
0.6 No Observations 90
0.5 80
Relative RMSE
0.4 70
0 5 10 15 20 25 30
Time(minutes)
0.3 100
Density(cars/mile)
(b2)Estimated maximal density
0.2
50
0.1
0 ! 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Time(minutes) Time(minutes)
Figure 5.11: Relative RSME (a) and parameter estimate (b) are shown for data assimilation of
real traffic data taken during late night. Observations are taken from all sensors (solid), half the
sensors (dotted), or no sensors (dashed).
107
Density(cars/mile) Velocity(miles/hour)
0.35
100
(a) Full Observations (b1)Estimated maximal velocity
Half Observations 95
0.3 No Observations
90
0.25
85
Relative RMSE
0.2 80
0 5 10 15 20 25 30
Time(minutes)
0.15
200
(b2)Estimated maximal density
0.1 150
0.05 100
0 ! 50
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Time(minutes) Time(minutes)
Figure 5.12: Relative RSME (a) and parameter estimate (b) are shown for data assimilation of
real traffic data collected during rush hour. Observations are taken from all sensors (solid), half
the sensors (dotted), or no sensors (dashed).
5.6 Discussion
We studied several aspects of traffic flow prediction: mathematical traffic models,
observation models, and data assimilation techniques. Specifically, we reviewed the
Lighthill-Whitham macroscopic model that is used as the basic underlying traffic
model in data assimilation (§5.2). We then discussed two types of observations:
Eulerian sensor data and Lagrangian GPS data, and developed a formulation that
is capable of assimilating the Eulerian and Lagrangian observations simultaneously
(§5.3). We also investigated how our approach could be used to estimate traffic
states as well as parameters using the EnKF and PF (§5.4).
The initial motivation of this paper was to propose an efficient approach that
could assimilate both Eulerian sensor data and Lagrangian GPS data simultane-
ously. The idea behind our approach is to append the differential equations for the
positions and velocities of the vehicles to the macroscopic traffic model in order to
solve them simultaneously. Unlike the previous methods, our approach is capable of
handling flux data, position data and velocity data efficiently without reformulation
108
in Lagrangian coordinates or reshaping the macroscopic traffic model. This approach
has been shown to accurately estimate traffic states in different traffic scenarios with
true underlying traffic model (Case 1, Case 4). We also studied how the choice of
filters and observations affect data assimilation: compared to the EnKF, the PF is
less sensitive to observation noise and sensor locations, but its computation cost is
higher (Case 2); the positions of sensor locations impact the accuracy of the EnKF,
and relevant suggestions for installing sensors are provided (Case 3). However, this
approach is less accurate when the true traffic model is unknown and an estimated
underlying traffic model is used (Case 5, Case 6).
One limitation we reported on in §5.5.3 is that the accuracy of data assimilation
is reduced significantly if the underlying traffic model cannot reproduce the actual
traffic flow from which observations are collected: in such a situation, the model
error is very large and will therefore dominate the overall error. Hence, to increase
accuracy in such cases requires the development of better, more accurate traffic flow
models that are capable of reproducing a wider range of traffic flows.
Another limitation is that the proposed approach for assimilating Lagrangian
observations requires knowledge of the velocity of a vehicle at a given position: in our
macroscopic model for the density, the velocity was given explicitly as a function of
the density, which allowed us to access the velocity in our algorithm. If the underlying
traffic-flow model does not provide information about the velocity of vehicles at any
specified position, our method for assimilating Lagrangian observations will not work.
There are various extensions that could be implemented to improve the frame-
work presented here. First, the proposed data assimilation approach was only applied
to simplified traffic scenarios such as on and off ramps, traffic lights, and bottlenecks
that were implemented in the LWR equation. In practice, the road network is typ-
109
ically more complicated and may include, for instance, intersections and multiple
lanes. An exciting extension would be to expand the above framework to more com-
plex and realistic traffic scenarios. Second, it would be interesting to see whether the
algorithm could be extended to use parameter estimation to relay traffic information
to drivers. For instance, the period of traffic lights could be estimated by data as-
similation, and thus the waiting time could be provided to drivers. Data assimilation
could also be used to predict the location and duration of congestion caused by slow
trucks or road constructions.
Chapter Six
Traffic control in lane-drop
scenarios
111
6.1 Introduction
A lane drop is defined as a location on a highway where the number of lanes provided
decreases. The reasons for a lane drop can be varied: the road construction may
close the usage of some lanes during road work; an unexpected traffic accident also
restricts the usage of the shoulder lane; the highway may reduce the number of lanes
when space is limited, and so on.
The merge bottleneck from lane drop is an interesting traffic phenomena. An
increase in lane-changing maneuvers will form queues in the upstream of the lane
drop, and the queues are propagating backward. Traffic was studied upstream and
downstream of a bottleneck in the lane drops. Recent empirical observations [10] at
freeway merge bottlenecks revealed: a drop in the bottleneck discharge rate when
queues form upstream, which is also called capacity drop. Further, it is shown that
the bottleneck’s discharge flow is about 10% lower than the prevailing flow observed
prior to queue formation. Upon bottleneck activation, flow reduction occurs sequen-
tially in time and space. In [53], numerical simulation results were in agreement with
the findings of the empirical studies at merge bottlenecks. In addition, it strongly
suggested that lane changes are the main cause of the drop in discharge rate.
Two types of approaches have been considered to reduce the capacity drop of
merge bottleneck. One type of approach attempts to develop cooperative lane-
changing models so that the lane-changing behavior is more efficient. Like in [41, 48],
cooperative lane-changing systems were developed microscopically in order to smooth
the merges. However, this microscopic approach is complicated to specify and there
are lot of parameters to quantify. The other type of approach [52] embeds the ve-
hicles as discrete particles in a multi-lane stream using the kinematic wave model.
112
It attempts to maximize capacity through efficiently controlling the lane-changing
ratio between adjacent lanes. However, this macroscopic approach formulates a con-
servation equation for individual lane.
In this chapter, we would like to introduce a macroscopic approach to model and
analyze the multi-lane traffic, by using only one continuity equation. Then we will
develop controlling strategies based on this continuity equation to maximize the road
capacity in lane drop. The chapter is organized as follows. In § 6.2, we describe the
inhomogeneous LWR equation that used to model multi-lane traffic, and analyze the
solutions to Riemann problems. In § 6.3, we develop two controlling strategies in
order to maximize road capacity in lane drop. Then we conclude with a discussion
of our results and open problems in § 6.4.
6.2 Inhomogeneous LWR model
In this section we will first briefly review the homogeneous LWR model, and then
introduce the inhomogeneous LWR model to describe the multi-lane traffic by taking
the number of lanes into consideration(§ 6.2.1). In addition, we analyze the solutions
to the Riemann problems of inhomogeneous LWR equation by applying the supply–
demand theory (§ 6.2.2). In § 6.2.3, we explore the traffic pattern on a ring road and
the efficiency of LWR equation to organize the road capacity.
The LWR model studies the traffic dynamics based on the conservation law. We
have described the LWR model in § 2.2. In this model, the traffic is described by
113
the vehicle density ρ and vehicle flux ϕ:
∂ρ ∂ϕ
+ = 0, (6.1)
∂t ∂x
where flux is a function of density via ϕ(ρ). The solutions to the Riemann problems
of homogeneous LWR model are already discussed in § 2.2 using the method of char-
acteristics, which include traveling wave, shock wave and rarefaction. In addition,
we generated multiple traffic scenarios by adding additional terms to capture the
respective effects.
The vehicle flow is defined as a function of traffic density ρ: ϕ = ϕ(ρ), which
is independent of the location x. Therefore, the LWR model is only for modeling a
homogeneous road, where flow function is uniform with respect to location x. That
is, the LWR model has limited modeling capacity when the traffic scenario becomes
complicated.
In this section, we will focus on the traffic scenarios with multiple lanes, where
the number of lanes are location dependent, and introduce the inhomogeneous LWR
model to account for the changes of the number of lanes.
6.2.1 Model description
We are interested in modeling the traffic scenarios with multiple lanes macroscopi-
cally. One well-known method to model multi-lanes was developed by Daganzo in
[21, 22, 53]. It formulated a conservation equation for each lane, and the right hand
side of each conservation equation is the net lane-changing rate onto the current lane
from adjacent lanes. The lane-changing rate consists of the sources from neighbor
114
lanes and sink to neighbor lanes, which depends on the nearby traffic at a specific
location and time.
However, the interaction between lanes can be very complicated to quantify in
lane-drop scenarios if we use a conservation equation for each lane. In addition, we
are only interested in the dynamics of the aggregated traffic characteristics instead of
the inner interaction between lanes. Therefore, we use another approach mentioned
in [73] (chapter 7.2), which formulates the averaged density and flux with only one
continuity equation:
∂(Iρ) ∂(Iϕ)
+ =0 (6.2)
∂t ∂x
where ρ is the averaged density per lane and ϕ is the averaged flux per lane. The
function I(x) is the number of lanes at location x, which models the the lane changes
before a lane closure or after a new lane open. The number of lanes is integer in
reality, but I(x) can be non-integer to account for partial lane usage. For example,
a value of I = 1.5 of a realistic 2-lane road indicates the second/shoulder lane is less
frequently used, and the flow on this lane is only half of the averaged flow on the
other lane.
By introducing the notations for total density and total flux
ρtotal (x, t) = I(x)ρ(x, t),
ϕtotal (x, t) = I(x)ϕ(x, t),
the conservation law (6.2) for multi-lane traffic can be rewritten in a form that is
similar to the homogeneous LWR model:
∂ρtotal ∂ϕtotal
+ = 0. (6.3)
∂t ∂x
115
Unlike the homogeneous LWR model, the total vehicle density ρtotal and total vehicle
flux ϕtotal are not directly related with respected to a fixed flux-density relation.
Instead, they are related corresponding to a lane-dependent flux-density relation:
ϕtotal (x, t) = ϕtotal (ρtotal , I(x))
ρtotal
= I(x)ϕ . (6.4)
I(x)
A lane-drop traffic scenario
In the whole chapter, we will focus on a simplified lane-drop traffic scenario where it
changes from two lanes to one lane. The analysis of solutions and traffic controlling
strategies in the following sections are also based on this specific lane-drop scenario.
Extension to generalized multi-lane traffic scenarios will be discussed in final section
§ 6.4. The lane-drop scenario is shown in the figure below:
I(x) = 2 I(x) = 1
0 x
Figure 6.1: Lane-drop scenario from two lanes to one lane.
The lane drop is at the origin with x = 0, and there are two lanes in the upstream
x < 0 and one lane in the downstream x > 0.
The fundamental diagram is a modeling choice, which is the plot of traffic flux
ϕ versus traffic density ρ. According to [64], the two-phase traffic theory divides
traffic flow into free flow for low densities, and congested flow for large densities,
116
which together form a reversed λ shape. The fundamental diagram we use to model
multi-lane traffic is shown in Figure 6.2. The shape of the fundamental diagram is
determined by the parameters from [24]. (see model calibration in § 4.2)
max
co
ng
vm flow
est
ed
ax
!
ee
flo
w
fr
⇢
⇢max
Figure 6.2: Fundamental diagram ϕ(ρ). The slope of the free flow region is the maximal
velocity vmax , and the slope of the congested flow is the congestion parameter ω. The
maximal density is ρmax , and the maximal flux is ϕmax .
Description Parameter Value Unit
Free-flow speed vmax 63.3 miles/hour
Maximal flux ϕmax 2031 vehicles/hour
Congestion speed parameter ω −10.1 miles/hour
Jam density ρmax 232 vehicles/mile
Table 6.1: Parameters from auto-calibration of the fundamental diagram
The critical density ρcrit is the divider of free flow and congested flow in the
fundamental diagram. It is easy to calculate the critical density via
ρcrit = Qmax /vmax ≈ 32.09vehicles/mile (6.5)
The critical density is a very important bifurcation point in multi-lane traffics. We
will talk about it with more details when analyzing the solution to multi-lane traffic
in the following subsection.
117
By using the Equation (6.4), we obtain the fundamental diagram for multiple-
lane traffic. It is of the same reversed λ shape, but scaled in both x-axis and y-axis
according to the lane number I(x). The Figure 6.3 shows the fundamental diagram
for the simplified lane-drop scenario in Figure 6.1.
total
2 I(x) = 1
max
I(x) = 2
ax
!
vm
max
⇢total
⇢max 2⇢max
Figure 6.3: The flux-density relation of ϕtotal (ρtotal , I(x)) in lane-drop scenario. The red solid
line represents the fundamental diagram for road of one lane, while the green dashed dotted line is
for road of two lanes.
6.2.2 Solutions to Riemann problems
Next, we will investigate the solutions to the inhomogeneous LWR equation. The
problem with discontinuous initial traffic states and lane number of the following
form:
ρ L ,
x≤0 2,
x≤0
ρtotal (x, 0) = , I(x) = . (6.6)
ρ R ,
x>0 1,
x>0
is called a Riemann problem to the inhomogeneous LWR model. In § 2.2, we have
reviewed the characteristic method, and used it to solve the homogeneous LWR
model. The solutions include traveling wave, shock wave and rarefaction. For the
inhomogeneous LWR model (6.2), there are some analytical and numerical studies
[42, 45, 47] to investigate the solutions to Riemann problems. In our study, we adopt
118
the method called wave separation, which was introduced in article [42]. Specifically,
the separation method includes:
(i) first, we introduce a stationary wave between the road branches of different
number of lanes. The stationary wave is subject to entropy condition. More
details about the entropy condition between road branches are explained as
follows.
(ii) then, we use the characteristic method to solve the homogeneous LWR model
for each road branch with constant number of lanes. The solution for each
road branch is also subject to the entropy condition. The entropy condition
for homogeneous LWR is listed earlier in § 2.2.
According to the article [42], the solutions to Riemann problems of the inhomoge-
neous LWR model exist and are unique if the following two entropy conditions are
satisfied:
Entropy 1: the waves from left (upstream) to right (downstream) should increase
their wave speeds so that they don’t cross each other.
Entropy 2: the standing wave can not cross the transition curve ρtotal /I(x) = ρcrit .
The left side of the transition curve is ρtotal /I(x) < ρcrit , and the right side of the
transition curve is ρtotal /I(x) ≥ ρcrit .
According to [47], the entropy condition can be transferred into a supply-demand
framework. Supply–demand framework is relatively easy to understand the traf-
fic flow. Another advantage about supply–demand framework is that it can guide
you to find standing wave between branches through analyzing the boundary flux.
119
Therefore, we will use it to construct the standing wave that satisfies the entropy
conditions, and then get the unique solution to the Riemann problem. Next let us
briefly review the supply–demand theory in traffic.
Supply and demand function
We start with showing the supply and demand function in multi-lane traffic scenarios.
The traffic supply function computes the traffic flow to supply for the upstream cell:
the supply of the free flow is the maximal flow, while the supply of the congested
flow is only itself. The mathematical expression is as follows:
S(ρtotal ) = I(x)ϕ(max(ρcrit , ρtotal /I(x)))
(6.7)
I(x)ϕ(ρcrit ),
ρtotal /I(x) < ρcrit
=
I(x)ϕ(ρtotal /I(x)),
ρtotal /I(x) ≥ ρcrit
The traffic demand function computes the traffic flow to demand from the down-
stream cell: the demand of the free flow is itself, while the demand of the congested
flow is the maximal flow. The mathematical expression is as follows:
D(ρtotal ) = I(x)ϕ(min(ρcrit , ρtotal /I(x)))
(6.8)
I(x)ϕ(ρtotal /I(x)),
ρtotal /I(x) ≤ ρcrit
=
I(x)ϕ(ρcrit ),
ρtotal /I(x) ≥ ρcrit
Let two adjacent cells be i and i + 1 with total vehicle density ρtotal,i and ρtotal,i+1 .
According to the supply–demand theory, the boundary flux is determined by the
supply from downstream and the demand from upstream. Specifically, the boundary
flux from cell i to cell i + 1 is the smaller of the supply of the cell i + 1 and the
120
demand of the cell i. This is in agreement to our intuition.
ϕ∗i,i+1 = min(S(ρtotal,i+1 ), D(ρtotal,i )) (6.9)
The Godunov method is used to solve the LWR model, which is a conservative
numerical scheme suggested by S. K. Godunov [33]. It is a conservative finite-volume
method which solves exact, or approximate Riemann problems at each inter-cell
boundary. The Godunov method is equivalent to the supply–demand theory when
calculating the boundary flux. The formula of Godunov method can be rewritten
as:
∆t ∗
ρj+1 j
ϕi,i+1 − ϕ∗i−1,i
total,i = ρtotal,i − (6.10)
∆x
Analytic solution to Riemann problem
It is convenient to find the unique solution to the Riemann problem by using the
idea of boundary flux. Given two initial vehicle density ρL and ρR in the Riemann
problem, we can first calculate the boundary flux at the origin x = 0. That is
ϕ∗L,R = min(S(ρR ), D(ρL )) by the supply–demand theory. As the standing wave has
speed 0 at the origin, the flux close to the lane drop must equals to the boundary
flux.
Next we have to figure out the two transitional vehicle densities in the upstream
ρA and downstream ρB of the standing wave. At least one of the transitional densities
121
can be determined by the calculated boundary flux as follows:
if ϕ∗L,R = D(ρL ) =⇒ ρA = min(ρcrit , ρL /I(x)) (6.11)
if ϕ∗L,R = S(ρR ) =⇒ ρB = max(ρcrit , ρR /I(x)) (6.12)
For instance, if ϕ∗L,R = D(ρL ) < S(ρR ), then we can figure out ρA by using the
formula above. Next, we have to specify the other transitional density ρB . First,
the vehicle flow corresponding to ρB should be the same as the boundary flow ϕ∗L,R
because of the standing wave at origin. As the fundamental diagram has a reversed
λ shape, there are at most two admissible densities for ρB . Then we can apply the
two entropy conditions Entropy 1 and Entropy 2 to pick the unique transitional
density for ρB .
The last step is to solve two Riemann problem for homogeneous LWR equations
in the upstream and downstream branches. The characteristic method can be applied
here. The two Riemann problems are:
ρL ,
x≤0
ρ(x, 0) = (6.13)
ρA ,
x = 0−
ρ ,
B x = 0+
ρ(x, 0) = (6.14)
ρR ,
x>0
The typical density solution to the Riemann problem is sketched in Figure 6.4.
There is a standing wave in the origin that is subject to entropy condition. Two
wave solutions are in upstream and downstream based on the characteristic method.
122
⇢A ⇢total
⇢L
⇢B
⇢R
I(x) = 2 0 I(x) = 1 x
Figure 6.4: The wave solution to Riemann problem of inhomogeneous LWR. ρL and ρR are the
initial conditions. ρA and ρB are the transitional traffic states in the upstream and downstream of
the standing wave.
Examples
We will show how to apply the wave separation method to solve the following Rie-
mann problem:
2ρ
crit , x<0 2,
x<0
ρtotal (x, 0) = , I(x) = . (6.15)
0,
x>0 1,
x>0
First, we have to compute the boundary flux by using the supply–demand theory:
D(ρL ) = 2ϕmax
S(ρR ) = ϕmax
ϕ∗L,R = min(S(ρR ), D(ρL )) = ϕmax
The computations show that the standing wave has boundary flux of ϕmax . The
downstream transitional vehicle density is ρB = max(ρR , ρcrit ) = ρcrit because of
ϕ∗L,R = S(ρR ). The upstream transitional vehicle density can be ρcrit or ρcrit + ρmax
because their respective vehicle fluxes equal the boundary flux ϕ∗L,R . By checking
the two entropy conditions, we find that only ρA = ρcrit + ρmax satisfies the entropy
conditions.
123
The initial traffic conditions as well as the transitional traffic states are indicated
in the fundamental diagram in Figure 6.5(a). The wave solution by using the sep-
aration method is illustrated in Figure 6.5(b). It is a shock wave between state L
ϕtotal (ρL ,2)−ϕtotal (ρA ,2)
and A, which is moving backward with speed ω = ρL −ρA
. The standing
wave is stationary at the origin. It is a traveling wave between state B and R, which
is moving forward with speed vmax . We will justify that this weak solution satisfies
both entropy conditions.
total
L I(x) = 1
I(x) = 2 ⇢total
A
B A
L
B
R
R ⇢total
I(x) = 2 0 I(x) = 1 x
(a) Traffic states in fundamental diagram. (b) Wave solution in density ρtotal .
Figure 6.5: The weak solution to Riemann problem (6.15) of inhomogeneous LWR model.
Entropy 1: the waves from left (upstream) to right (downstream) should increase
their wave speeds so that they don’t cross each other. The waves’ speeds from
upstream to downstream are ω, 0 and vmax . It is obvious that the relation
ω < 0 < vmax holds.
Entropy 2: the standing wave can not cross the transition curve ρtotal /I(x) = ρcrit .
The two transitional states ρA and ρB around the standing wave both satisfy
ρtotal /I(x) ≥ ρcrit . That means they do not across the transition curve.
124
6.2.3 Traffic pattern on a ring road
With the ability to solve the Riemann problems of inhomogeneous LWR model, we
can explore the traffic in more complicated road structure. We are most interested
in the traffic pattern by the inhomogeneous LWR on a ring road with two branches:
one branch is of one lane, and the other branch is of two lanes. On one side, we could
observe how the inhomogeneous LWR model generates congestion pattern. On the
other side, we would like to understand how LWR model efficiently organizes and
distributes traffic.
The ring road profile is shown in Figure 6.16. There are two road branches: the
one-lane branch is of length L1 and the two-lane branch is of length L2 . Here, we
just set the length ratio L2 : L1 = 1 : 1 for investigating the traffic pattern, and then
make a conclusion for a general case in the end of this subsection.
L2 : L1 = 1 : 1
Figure 6.6: A ring road with two branches: one-lane branch with length L1 and two-lane branch
with length L2 . The ratio is L2 : L1 = 1 : 1.
The initial densities on the ring road are described in equation (6.16) and Fig-
ure 6.7. This Riemann problem has periodic boundary condition at x = L1 / − L2
125
(we will use x = L1 to represent the boundary in the following discussion).
ρ 2 ,
−L2 < x < 0 2,
−L2 < x < 0
ρtotal (x, 0) = , I(x) = . (6.16)
ρ 1 ,
0 < x < L1 1,
0 < x < L1
Road capacity refers the maximal traffic flow (vehicles per hour) obtainable on a
⇢total
⇢2
⇢1
x
L2 I(x) = 2 0 I(x) = 1 L1
Figure 6.7: Initial conditions of the Riemann problem on a ring road with two branches. The
first branch is of two lanes from −L2 to 0, and the second branch is of one lane from 0 to L1 . The
boundary is periodic.
given road using all available lanes. Therefore, the two-lane branch has double road
capacity of the one-lane branch. The road capacity of this ring road is restricted by
the one-lane branch, so the overall road capacity is only ϕmax . In the fundamental
diagram of the ring road, there are two points that can provide the exact road
capacity ϕmax . Please see Figure 6.8 for positions of these two critical points C1 and
C2 . The analysis shows that the averaged density (ρ1 + ρ2 )/2 determines the traffic
pattern of the ring road, and there are three regimes controlled by the points C1 and
C2 . The three regimes are:
ρ1 + ρ2
Regime 1: ≤ ρC1 (6.17)
2
ρ1 + ρ2 ρC + ρC2
Regime 2: ρC1 < ≤ 1 (6.18)
2 2
ρ1 + ρ2 ρC1 + ρC2
Regime 3: > (6.19)
2 2
with ρC1 =ρcrit , ρC2 = ρcrit + ρmax
126
total
2 I(x) = 1
max
I(x) = 2
C1 C2
max
⇢total
⇢crit ⇢crit + ⇢max
Figure 6.8: The fundamental diagram of the ring road. There are two points C1 and C2 which
can provide the road capacity ϕmax .
Next we will show the traffic pattern in each regime as well as provide the analytic
solution to the Riemann problems:
ρ1 +ρ2
• Regime 1: 2
≤ ρ C1
ρC1 is the critical density in the fundamental diagram, where maximal traffic
flux is achieved. In a one-lane road, when the vehicle density is less than
ρC1 , it is free flow; when the vehicle density is above ρC1 , it is congested flow.
Therefore, the intuition tells us that it should be free flow everywhere if the
initial average density is not more than the critical density ρC1 .
ρ1 +ρ2
For the Regime 2
≤ ρC1 , there is at least one density (ρ1 or ρ2 ) should be
less than or equal to the critical density ρC1 . There are two sub cases for this
Regime 1: (a) both ρ1 and ρ2 are less than or equal to the critical density; (b)
only one density of ρ1 and ρ2 is less than or equal to the critical density.
– Case 1: ρ1 ≤ ρC1 and ρ2 ≤ ρC1
This is a trivial case. When the initial densities on both road branch are
less than or equal to the critical density ρC1 , it is a traveling wave with
127
maximal speed vmax . The traveling wave consists of two traffic states:
ρ1 and ρ2 . The distance between the two discontinuities are exactly the
length L1 of road branch.
In Figure 6.9, we show the traveling wave solution to the following Rie-
mann problem (6.20) in which both ρ1 and ρ2 are less than or equal to
the critical density ρC1 . The traveling wave in Figure 6.9 has two density
states: 0.5ρC1 and 0.8ρC1 . Even though we switch the initial density for
two branches, the stationary wave solution is the same.
0.5ρ
C1 0.8ρC1 , −L2 < x < 0
ρtotal (x, 0) = or . (6.20)
0.8ρC1
0.5ρC1 , 0 < x < L1
⇢total
⇢C 1
vmax vmax
x
vave
vmax
x
total
max
vmax vmax
x
L2 I(x) = 2 0 I(x) = 1 L1
Figure 6.9: The weak solution to Riemann problem (6.20) on a ring road. The plots are vehicle
density, velocity and flux from top to bottom.
– Case 2: ρ1 > ρC1 or ρ2 > ρC1
In this case, one of the branches has an initial density above the critical
density ρC1 . Even though the analytic solution is more complicated in the
128
beginning, the stationary solution is still a traveling wave with maximal
speed vmax . The traveling wave consists of two traffic states: min(ρ1 , ρ2 )
and ρC1 . The distance between the two discontinuities is not necessarily
the length of roach branch, but is uniquely determined by the conservation
of vehicles. Let D be the length with density ρC1 , then it must satisfy the
relation ρ1 L1 + ρ2 L2 = min(ρ1 , ρ2 )(L1 + L2 − D) + ρC1 D.
In Figure 6.10, we show that the traveling wave solution to the following
Riemann problem (6.21) in which only one of ρ1 and ρ2 is less than or
equal to the critical density ρC1 . The traveling wave in Figure 6.9 has two
density states: 0.5ρC1 and 1ρC1 . The length with density ρC1 is 1.4L1 .
Even though we switch the initial density for two branches, the stationary
wave solution is the same.
1.2ρ
C1 0.5ρC1 , −L2 < x < 0
ρtotal (x, 0) = or . (6.21)
0.5ρC1
1.2ρC1 , 0 < x < L1
⇢total
⇢C 1
vmax vmax
x
vave
vmax
x
total
max
vmax vmax
x
L2 I(x) = 2 0 I(x) = 1 L1
Figure 6.10: The weak solution to Riemann problem (6.21) on a ring road. The plots are vehicle
density, velocity and flux curves from top to bottom.
129
ρ1 +ρ2 ρC1 +ρC2
• Regime 2: ρC1 < 2
≤ 2
ρC1 is the critical density in the fundamental diagram of one-lane road where
maximal traffic flux is achieved. In the fundamental diagram of two-lane road,
there are two densities ρC1 and ρC2 whose respective fluxes are exactly the
maximal traffic flux ϕmax of the one-lane road. When the average road density
is above the critical density ρC1 , it is impossible to see free flow everywhere.
Congestion exists somewhere.
ρ1 +ρ2 ρC1 +ρC2
For the condition 2
≤ 2
, there is at least one inequality of ρ1 ≤ ρC1
and ρ2 ≤ ρC2 hold. When ρ1 ≤ ρC1 , the capacity supply of the one-lane
road branch is ϕmax . When ρ2 ≤ ρC2 , the capacity supply of the two-lane
road branch is above ϕmax . Therefore, there is at least one road branch that
can provide capacity ϕmax , and generate free flow with density ρC1 . Then the
free flow stays on the road branch of one lane after traveling. In contrast,
congestion occurs on the road branch of two lanes in the upstream of the
lane drop location (origin x = 0). Because the average density is bounded by
ρC1 +ρC2
2
, the congestion will only occur on the road branch of two lanes, and
not propagate backward to the branch boundary x = L1 . The free flow on the
road branch of one lane is not influences by the congestion.
The stationary weak solution is a standing wave, which consists of two traffic
states: ρC1 and ρC2 . The density ρC2 only occurs on the road branch of two
lanes in the upstream of the lane drop location. The two discontinuities are
both standing waves. The length of the density ρC2 is uniquely determined
by the conservation of vehicles. Let D be the length with density ρC2 , then it
must satisfy the relation ρ1 L1 + ρ2 L2 = ρC1 (L1 + L2 − D) + ρC2 D. Even though
it is not free flow everywhere, we can observe the maximal vehicle flux ρmax
everywhere.
130
The parameter setting of the fundamental diagram indicates ρC2 ≈ 8.23ρC1 .
ρ1 +ρ2
Then the Regime 2 becomes ρC1 < 2
≤ 4.615ρC1 . In Figure 6.11, we show
the stationary wave solution to the following Riemann problem (6.22) in which
ρ1 + ρ2 = 7ρC1 . The stationary wave solution in Figure 6.11 has two density
states: ρC1 and ρC2 . Even though we switch the initial density for two branches,
the stationary wave solution is the same.
3ρ
C1 4ρC1 , −L2 < x < 0
ρtotal (x, 0) = or . (6.22)
4ρC1
3ρC1 , 0 < x < L1
⇢total
⇢C 2
⇢C 1
x
vave
vmax
x
total
max
x
L2 I(x) = 2 0 I(x) = 1 L1
Figure 6.11: The weak solution to Riemann problem (6.22) on a ring road. The plots are vehicle
density, velocity and flux curves from top to bottom.
ρ1 +ρ2
The extreme situation is shown in Figure 6.12 with initial condition 2
=
4.615ρC1 . The stationary wave has two density states: ρC1 and ρC2 . The
discontinuities are standing with zero speed. The congestion makes up the
whole road branch of two lanes.
ρ1 +ρ2 ρC1 +ρC2
• Regime 3: 2
> 2
As mentioned earlier, in the fundamental diagram of two-lane road, there are
131
⇢total
⇢C 2
⇢C 1
x
vave
vmax
x
total
max
x
L2 I(x) = 2 0 I(x) = 1 L1
Figure 6.12: The weak solution to Riemann problem ( ρ1 +ρ 2
2
= 4.615ρC1 ) on a ring road. The
plots are vehicle density, velocity and flux curves from top to bottom.
two densities ρC1 and ρC2 whose respective fluxes are exactly the maximal
traffic flux ϕmax of the one-lane road. When the average road density is above
ρC1 +ρC2
2
, it is impossible to see a sustainable free flow on the road branch of
one lane. That is, congestion could exist everywhere on the ring road.
The stationary solution consists of two shock waves, which are moving back-
ward with speed ω on each road branch. The shock wave evolves into another
shock wave when it hits the origin (x = 0) or the boundary (x = L1 ). In total,
there are four types of shock waves even though we can observe only two shock
waves at a time. For each road branch, the shock discontinuity consist of two
fixed density states, which are uniquely determined by the initial conditions.
Next, we will introduce the shock waves case by case.
ρ1 +ρ2 ρC1 +ρC2
For the condition 2
> 2
, there is at least one inequality of ρ1 ≥ ρC1
and ρ2 ≥ ρC2 hold. The two cases for Regime 3 are: (a) both inequalities hold
for ρ1 and ρ2 ; (b) only one inequality holds for ρ1 or ρ2 .
– Case 1: ρ1 ≥ ρC1 and ρ ≥ ρC2
132
The road capacity of both branches are restricted by the supply ϕtotal (ρ1 , 1)
and ϕtotal (ρ2 , 2). Then the stationary solution has a capacity switch be-
tween ϕtotal (ρ1 , 1) and ϕtotal (ρ2 , 2) on each road branch. There are two
corresponding densities called ρ˜1 and ρ˜2 which satisfy
ϕtotal (ρ1 , 1) = ϕtotal (ρ˜2 , 2),
ϕtotal (˜
ρ1 , 1) = ϕtotal (ρ2 , 2).
Please see Figure 6.13 for the positions of these densities in fundamental
diagram. The density ρ˜2 has equal flux to density ρ1 , and the density
ρ˜1 has equal flux to density ρ2 . The shock wave on the road branch of
one lane is made up with ρ1 and ρ˜1 , and the shock wave on the road
branch of two lanes is made up with ρ2 and ρ˜2 . In total, there are four
possible types of shock waves on the ring road, and they will disappear
and regenerate when hitting the origin or boundary. The evolution rules
are: shock(ρ2 , ρ˜2 )↔ shock(ρ1 , ρ˜1 ) and shock(˜
ρ2 , ρ2 )↔ shock(˜
ρ1 , ρ1 ).
total
2 I(x) = 1
max
I(x) = 2
C1 C2
max
⇢total
⇢1 ⇢˜1 ⇢˜2 ⇢2
Figure 6.13: The fundamental diagram of the ring road. The density ρ˜2 has equal vehicle flux to
the given initial density ρ1 . The density ρ˜1 has equal vehicle flux to the given initial density ρ2 .
In Figure 6.14, we show the shock wave solution to the following Riemann
problem (6.23) in which ρ1 > ρC1 and ρ2 > ρC2 . Both shock waves are
133
moving with speed ω. When the shock wave hits the origin (x = 0) or
the boundary (x = L1 ), it disappears on the current road branch, but
evolves into another shock on the other road branch. In this case, the
distance between the two observed shocks are exactly the length L1 of
roach branch because the two shocks waves are generated simultaneously.
ρ
C2 + 2ρC1 , −L2 < x < 0
ρtotal (x, 0) = . (6.23)
2ρC1 ,
0 < x < L1
⇢total
⇢2
! ⇢˜2 ⇢˜1
⇢1 !
⇢C 1
x
vave
vmax
!
!
x
total
max
! !
x
L2 I(x) = 2 0 I(x) = 1 L1
Figure 6.14: The weak solution to Riemann problem (6.22) on a ring road. The plots are vehicle
density, velocity and flux curves from top to bottom.
– Case 2: ρ1 < ρC1 or ρ < ρC2
In this case, one of the inequality ρ1 < ρC1 or ρ < ρC2 holds. Even though
the analytic solution is more complicated at the beginning, the stationary
solution still consists of two shocks wave. Because of the inequality ρ1 <
ρC1 or ρ < ρC2 , the capacity supply could somehow reach the maximal
flux ϕmax . Without loss of generosity, we just assume ρ1 < ρC1 fo the
following analysis. Then the stationary solution has a capacity switch
134
between ϕmax and ϕtotal (ρ2 , 2). The four densities states in the stationary
solution are ρC1 , ρC2 , ρ2 and ρ˜1 :
ϕtotal (ρC1 , 1) = ϕtotal (ρC2 , 2),
ϕtotal (˜
ρ1 , 1) = ϕtotal (ρ2 , 2).
In Figure 6.15, we show the piecewise shock wave solution to the following
Riemann problem (6.24) in which ρ1 ≤ ρC1 and ρ2 > ρC2 . Both shock
waves are moving with speed ω. When the shock wave hits the origin
(x = 0) or the boundary ( x = L1 ), it disappears on the current road
branch, but evolves into another shock on the other road branch. The
evolution rules are: shock(ρ2 , ρC2 )→ shock(ρC1 , ρ˜1 ) and shock(ρC2 , ρ2 )→
shock(˜
ρ1 , ρC1 ). In this case, the distance between the two observed shocks
are not necessary the length L1 of roach branch, but is uniquely deter-
mined by the conservation of vehicles. Let D be the length with capacity
ρC1 +ρC2
ϕmax , then it must satisfy the relation ρ1 L1 +ρ2 L2 = 2
D+ ρ2 +˜
2
ρ1
(L1 +
L2 − D).
ρ
C2 + 2ρC1 , −L2 < x < 0
ρtotal (x, 0) = . (6.24)
0.5ρC1 ,
0 < x < L1
Generalization of traffic pattern on a ring road
The discovery of traffic pattern on a ring road can be generalized to any ring road pro-
file. The general road profile is shown in Figure 6.16. There are two road branches:
the k1 -lane branch is of length Lk1 , and the k2 -lane branch is of length Lk2 with the
135
⇢total
⇢2
! ⇢C 2 ⇢˜1
⇢C 1
⇢C 1 !
x
vave
vmax
!
!
x
total
max
! !
x
L2 I(x) = 2 0 I(x) = 1 L1
Figure 6.15: The weak solution to Riemann problem (6.24) on a ring road. The plots are vehicle
density, velocity and flux curves from top to bottom.
length ratio Lk2 : Lk1 = λ2 : λ1 . The Riemann problem is described in (6.25) with
periodic boundary condition at x = Lk1 or −Lk2 .
L k2 : L k1 = 2 : 1
Figure 6.16: A ring road with two branches: k1 -lane branch with length Lk1 and k2 -lane branch
with length Lk2 . The ratio is Lk2 : Lk1 = λ2 : λ1 .
ρ 2 ,
−Lk2 < x < 0 k2 ,
−Lk2 < x < 0
ρtotal (x, 0) = , I(x) = . (6.25)
ρ 1 ,
0 < x < Lk1 k1 ,
0 < x < L k1
136
Without loss of generosity, we just assume k2 > k1 . The road capacity of the
k1 -lane branch is k1 ρmax . In the fundamental diagram of the ring road, there are two
points that can provide exact k1 ρmax capacity. These two points are called C1 and
λ1 ρ1 +λ2 ρ2
C2 . The averaged density λ1 +λ2
determines the traffic pattern of the ring road.
There are three regimes controlled by the points C1 and C2 .
λ1 ρ1 + λ2 ρ2
Regime 1: ≤ ρ C1 (6.26)
λ1 + λ2
λ1 ρ1 + λ2 ρ2 λ1 ρC1 + λ2 ρC2
Regime 2: ρC1 < ≤ (6.27)
λ1 + λ2 λ1 + λ2
λ1 ρ1 + λ2 ρ2 λ1 ρC1 + λ2 ρC2
Regime 3: > (6.28)
λ1 + λ2 λ1 + λ2
with ρC1 =k1 ρcrit , ρC2 = k1 ρcrit + (k2 − k1 )ρmax
6.3 Controlling strategies
In this section, we will develop two controlling strategies based on the multi-lane
LWR model such that the wave solution can provide smooth velocity solution without
losing the road capacity. The vehicles are embedded in this macroscopic stream and
their movements obey the velocity solution.
The study of traffic pattern on a ring road indicates that the multi-lane LWR
model is efficient in organizing traffic density, and thus capable to maximize the
capacity. In order to better show our proposed controlling strategies, we use a traffic
scenario to explain: there is a straight road with a lane drop at the origin (x = 0).
A wave of vehicles (density 2ρcrit ) is approaching the lane drop. The scenario can be
137
mathematically formulated as:
2ρ
crit , x << 0 2,
x<0
ρtotal (x, 0) = ρ0 (x) = , I(x) = . (6.29)
0,
otherwise 1,
x>0
Then complete LWR model as well as the initial conditions and boundary conditions
is as follows:
∂ρtotal ∂ϕtotal
+ =0 (6.30)
∂t ∂x
IC: ρtotal (x, 0) = ρ0 (x) (6.31)
BC: ρtotal (−∞, t) = 2ρcrit , ρtotal (∞, t) = 0 (6.32)
Let us analyze the wave solution to the equation above. First, it is equivalent to a
Riemann problem of homogeneous LWR equation in the upstream branch (x < 0).
The solution is a traveling wave with speed vmax . When the traffic density 2ρcrit
reaches the lane drop (x = 0), it becomes a Riemann problem of inhomogeneous
LWR equation with initial conditions: ρL = 2ρcrit and ρR = 0. We have figured out
the unique solution that satisfy the entropy conditions. (Please see § 6.2.2 for the
solution analysis.) The evolution of the solution is sketched in the Figure 6.17.
The stationary solution consists of three waves: the shock wave ρL − ρA moving
backward with speed ω, the standing wave ρA −ρB at the lane drop, and the traveling
wave ρB − ρR moving forward with speed vmax .
With the density solution to the multi-lane LWR equation, we can compute the
respective velocity/flux solution by using the velocity/flux function. By taking the
138
⇢total ⇢total
⇢L ⇢L
vmax
⇢R
I(x) = 2 0 I(x) = 1 x I(x) = 2 0 I(x) = 1 x
⇢A ⇢total ⇢A ⇢total
! !
⇢L ⇢L
⇢B ⇢B
vmax vmax
⇢R ⇢R
I(x) = 2 0 I(x) = 1 x I(x) = 2 0 I(x) = 1 x
Figure 6.17: The evolution of the solution to Riemann problem (6.29). In the top figures, a
traveling wave moving forward with speed vmax . In the bottom figures, the solution consists of a
shock wave moving backward, a standing wave at the origin and a traveling wave moving forward.
lane number I(x) into consideration, the velocity and flux functions become:
ρtotal
ϕtotal (ρtotal , I(x)) = I(x)ϕ
I(x)
ϕtotal ,
ρtotal
ρtotal 6= 0
vave (ρtotal , I(x)) =
vmax ,
ρtotal = 0
We plot the velocity and flux curves together with the density solution to the multi-
lane LWR equation in Figure 6.18. From the flux curve in the bottom, we observe
that the total flux ϕtotal equals to the road capacity ρmax in the downstream even
though there is a queue in the upstream. In other words, the phenomena of capacity
drop is not observed in the multi-lane LWR model.
This is an exciting discovery. If the vehicles on the road can adjust the speed
based on the velocity solution from multi-lane LWR model and use zipper merging,
then there is no so-called capacity drop. This maximizes the capacity usage of
highways, and efficiently avoids delays. It is not hard to realize if all the vehicles are
self-driving cars or they can communicate with a central control system. The vehicle
at position x and time t moves with speed vave (x, t) of the velocity solution. Because
139
⇢total
⇢A
!
⇢L
⇢B
vmax ⇢R
x
vave
vL vmax vB vR
vA
x
total
max
L
A B
R
x
I(x) = 2 0 I(x)R= 1
Figure 6.18: The wave solution to (6.29). The plots are density curve, velocity curve and flux
curve from top to bottom.
the density–velocity relation follows the fundamental diagram, the vehicles always
stay at a safe velocity based on current headways. Therefore, we are interested
in guiding vehicles to go through the lane-drop bottleneck by obeying the velocity
solution from the multi-lane LWR.
However, the velocity curve generated by the multi-lane LWR system is not
applicable for a vehicle to follow because of the huge velocity jumps. Even though,
the self-driving cars can adjust speed fast enough to follow the velocity curve, it
definitely sacrifices the comfort of the passengers by using extreme acceleration or
deceleration.
Therefore, our goal is to build a multi-lane LWR system, which can generate
efficient and comfortable average velocity vave (x, t) for vehicles to follow while max-
imizing the capacity. The desired average velocity vave (x, t) is a smoothly changing
curve that avoids extreme acceleration and deceleration. Here, we still use the fun-
damental diagram to restrict the velocity–density relation in order to guarantee the
140
safety. There are two potential controlling strategies to generate a comfortable aver-
age velocity: (a) lower the velocity jumps vL /vA and vB /vA ; (b) smooth the velocity
jumps vL /vA and vB /vA .
In the next two subsections, we will show the controlling strategies through pro-
viding the traffic explanation and mathematical formula.
6.3.1 Virtual lane usage
When the vehicles approach the lane drop, the vehicles on the dropping lane will
merge into other existing lanes gradually. The common senses tell us that drivers
should reduce the usage of the dropping lanes and merge into other lanes earlier.
Otherwise, the vehicles have to cut in by forcing the vehicles on the target lane to
stop. The rushed cut-in can disrupt the ongoing traffic flow and bring stop-and-go
wave.
The idea of early merge when approaching the lane drop can be mathematically
understood as virtual lane usage. In the original multi-lane LWR model, the number
of lanes I(x) jumps from 2 to 1 immediately at the origin (x = 0). It means drivers
are still taking full use of the shoulder lane until the lane drop. As explained in
§ 6.2, I(x) can be non-integer to account for partial lane usage. For example, a
value of I(x) = 1.5 represents the shoulder lane is less frequently used. In order to
demonstrate the idea of virtual lane usage, the lane number I(x) should gradually
decreases from 2 to 1.
We pick a road strip of length L in the upstream of the lane drop, and set the
lane number I(x) as 1.5 instead of 2. The idea is shown in Figure 6.19. As shown in
141
Figure 6.19(a), the density gap between state L and A is determined by the difference
of lane number I(x). When we pick a road strip of 1.5 lanes in Figure 6.19(b), there
are more transitional states between the gap L/A. The fundamental diagram of the
virtual 1.5 lanes is in the orange dashed line.
total total
L I(x) = 1 L I(x) = 1
I(x) = 2 I(x) = 2
E C
B A B D A
R ⇢total R ⇢total
(a) Original fundamental diagram (green (b) Introduced fundamental diagram (or-
dashed dotted line) for one-lane and two- ange dashed line) for 1.5-lane traffic in the
lanes traffic in multi-lane LWR model. upstream of the lane drop.
Figure 6.19: Comparison of fundamental diagram: (a) lane number I(x) decreases from 2 to 1;
(b) lane number I(x) decreases from 2 to 1.5, then to 1.
Next let us analyze the wave solution to the multi-lane LWR system with virtual
lane usage. There are two lane drops in the road profile now. First, the solution
is a traveling wave with speed vmax in the upstream branch (x < −L). Then, it
becomes a Riemann problem of inhomogeneous LWR model with initial conditions:
2ρcrit and 0 in the branch (−∞, 0). By using the supply–demand theory, it is not
hard to compute the boundary flux and figure out the two transitional states C and
E. The density E becomes a new traveling wave with speed vmax in the branch
(−L, 0). When the state E arrives the second lane drop x = 0, another Riemann
problem of inhomogeneous LWR model is generated with initial conditions: ρE and
0. By using the supply–demand theory, we figure out the two transitional states D
and B. The shock wave E/D move backwards with speed ω, and evolves into shock
wave C/A across the first lane drop (x = −L). The wave solution, including density,
velocity and flux, is plotted in Figure 6.20.
142
⇢total
⇢A
⇢C !
⇢L ! ⇢D
⇢B
vmax
⇢R
x
vave
vL vmax vB vR
vD
vC
vA
x
total
C max
L
A D B
R
x
I(x) = 2 I(x) = 1.5 I(x) = 1
Figure 6.20: The wave solution to (6.29) after introduce a virtual lane I(x) = 1.5. The plots are
density curve, velocity curve and flux curve from top to bottom.
There are two shock waves L/C and C/A moving backward with the same speed
ω. There are two standing wave at the lane drop location x = 0 and x = −L0 .
It is a traveling wave between state B and R in the downstream, which is moving
forward with speed vmax . We will justify that this weak solution satisfies both entropy
conditions.
Entropy 1: the waves from left (upstream) to right (downstream) should increase
their wave speeds so that they don’t cross each other. The waves’ speeds from
upstream to downstream are still ω, 0 and vmax . It is obvious that the relation
ω < 0 < vmax holds.
Entropy 2: the standing wave can not cross the transition curve ρtotal /I(x) = ρcrit .
The two transitional states ρA and ρD at lane drop x = 0 satisfy ρtotal /I(x) ≥ ρcrit .
The two transitional states ρD and ρB at lane drop x = 0 satisfy ρtotal /I(x) ≥ ρcrit .
That means that they do not across the transition curve.
143
6.3.2 Reshaped fundamental diagram
When the vehicles approach the lane drop, congestions occur because of non-uniform
lane-changing behaviors. The common senses tell us that drivers should slow down
when approaching the lane drop in attempt to leave enough headway for the merging
purpose. For example, drivers usually drive 50mph when the headway is 20m. Then
the drivers should drive only 40mph with the same headway when there is a lane
drop ahead. This action can avoid over-congestion so that vehicles are relatively
further away from each other.
The change of driving habit when approaching the lane drop can be mathemat-
ically understood as the change of the velocity–density relation in the regime of
congested flow in the fundamental diagram. The change of velocity–density relation
is equivalent to the change of the flux–density relation, so we can reshape the flux-
density relation around the lane drop area. Therefore, an alternative fundamental
diagram is introduced as the controlling strategy in order to narrow the density gap
ρL /ρA , and thus narrow the velocity gap vL /vA .
The comparison is shown in Figure 6.21. As shown in Figure 6.21(a), the density
gap between state L and A is determined by the shape of the flux-density relation
in the congested flow. When we propose an alternative flux-density relation in Fig-
ure 6.21(b) under the original relation, the density gap is narrowed. The alternative
flux-density relation is in the orange dashed line.
Next let us analyze the wave solution to the multi-lane LWR system with the
modified flux-density relation. First, the solution is a traveling wave with speed
vmax in the upstream branch (x < 0). Then, it becomes a Riemann problem of
inhomogeneous LWR model with initial conditions: ρL = 2ρcrit and ρR = 0. By
144
total total
L I(x) = 1 L I(x) = 1
I(x) = 2 I(x) = 2
B A B A0 A
R ⇢total O R ⇢total
(a) Original fundamental diagram (green (b) Alternative fundamental diagram (or-
dashed dotted line) for two-lane traffic in ange dashed line) for two-lane traffic near
multi-lane LWR model. lane drop.
Figure 6.21: Comparison of fundamental diagram: (a) original shape with transitional state A;
(b) alternative shape with transitional state A0 .
using the supply–demand theory, it is not hard to compute the boundary flux and
figure out the two transitional states A0 and B. The wave solution, including density,
velocity and flux, is plotted in Figure 6.22.
⇢total
⇢A 0
⇢L !0
⇢B
vmax ⇢R
x
vave
vL vmax vB vR
v A0
x
total
max
L
A0 B
R
x
I(x) = 2 0 I(x) = 1
Figure 6.22: The wave solution to (6.29) after reshaping the fundamental diagram. The plots are
density curve, velocity curve and flux curve from top to bottom. The original wave solution is in
blue dashed line.
It is a shock wave between state L and A0 , which is moving backward with speed
ϕL −ϕA0
ω0 = ρL −ρA0
. The standing wave is stationary at the origin. It is a traveling wave
between state B and R, which is moving forward with speed vmax . We will justify
145
that this weak solution satisfies both entropy conditions.
Entropy 1: the waves from left (upstream) to right (downstream) should increase
their wave speeds so that they don’t cross each other. The waves’ speeds from
upstream to downstream are ω 0 , 0 and vmax . It is obvious that the relation
ω 0 < 0 < vmax holds.
Entropy 2: the standing wave can not cross the transition curve ρtotal /I(x) = ρcrit .
The two transitional states ρA0 and ρB around the standing wave both satisfy
ρtotal /I(x) ≥ ρcrit . That means that they do not across the transition curve.
6.3.3 Optimization in controlling strategies
In the last two parts, we introduced idea of traffic control via reshaped fundamental
diagram and virtual lane usage. In this part, we would like to provide more insights
for the parameter setting in the two controlling strategies.
virtual lane usage
The goal of virtual lane usage is to smooth the velocity jump so that the vehicles are
comfortable to follow. Consider the velocity jump that consists the higher velocity
vL and the lower velocity vA . The jump vL /vA is moving backward with speed ω,
and the jump vA /vL is standing with speed zero. Please see Figure 6.23 for lane
number I(x) and average velocity vave (x).
First, we investigate the comfort requirement on the velocity solution vave , then
146
we attempt to obtain the comfort requirement on the virtual lane I(x) through the
relation between vave (x) and I(x).
I
2
1
s
x
vave
vL vL
!
vA
sLA sAL
x
Figure 6.23: Lane number I(x) (top) and average velocity vave (x) (bottom) using virtual lane
usage.
The comfortable acceleration and deceleration are denoted by a and b (see § 4.2).
The comfort requirements on velocity vave include: the deceleration from vL to vA
should not be greater than |b|; the acceleration from vA to vL should not be greater
than a. The deceleration can be computed as:
vave (x) − vave (x + (vave (x) + ω)∆t)
lim
∆t→0 ∆t
0
vave (x) − vave (x) − (vave (x) + ω)vave (x)∆t (6.33)
= lim
∆t→0 ∆t
0
= − (vave (x) + ω)vave (x)
The acceleration can be computed as:
vave (x + vave (x)∆t) − vave (x)
lim
∆t→0 ∆t
0
vave (x) + vave (x)vave (x)∆t − vave (x) (6.34)
= lim
∆t→0 ∆t
0
=vave (x)vave (x)
147
Therefore, the requirements on the velocity vave can be written as:
0
(vave (x) + ω)vave (x) ≤ |b| (6.35)
0
vave (x)vave (x) ≤ a (6.36)
In Figure 6.23, let s represent the span function I(x) takes from 2 to 1, sLA represent
the span from velocity vL to vA , and sAL represent the span from velocity vA to vL .
The length of the spans satisfy the following relation:
sAL = 1s (6.37)
vmax + |w|
sLA = s (6.38)
vmax
The shapes of sAL and sLA in velocity vave (x) are determined by the lane number
I(x) via the fundamental diagram. Then we would like to derive the relation between
lane number I(x) and the shapes of sAL and sLA by using the fundamental diagram.
The relation is shown in Figure 6.24.
The virtual lane I(x) = k will bring two transitional states: state C in slope
sLA and state D in slope sAL . The density and flux of these states can be inferred
form the fundamental diagram, and then we can compute the velocity by using
vave = ϕtotal /ρtotal :
kϕmax
ϕC = kϕmax , ρC = kρcrit + (2 − k)ρmax =⇒ vC =
kρcrit + (2 − k)ρmax
ϕmax
ϕD = ϕmax , ρD = ρcrit + (k − 1)ρmax =⇒ vD =
ρcrit + (k − 1)ρmax
By incorporating the scaling factor in Equation (6.37) and (6.38) into consideration,
the relation between velocity lane number I(x) and velocity vave (x) can be explicitly
expressed as follows. We are only interested in the shapes of sAL and sLA , so neither
148
I(x) I(x)
k k
vL vD vL
vC
vave (x) vave (x)
vA vA
(a) slope sLA (b) slope sAL
total
L I(x) = 1
I(x) = 2
C
B D A
O ⇢total
(c) fundamental diagram
Figure 6.24: The relation between lane number I(x) and velocity vave (x). (a)&(b) The slopes
sLA and sAL corresponding to lane number I(x); (c) fundamental diagram with virtual lane.
the time term nor the position shift is included in the expression.
I(x)ϕmax
vave (x) = (6.39)
I(x)ρcrit + (2 − I(x))ρmax
vmax + |w|
ϕmax
vave x = (6.40)
vmax ρcrit + (I(x) − 1)ρmax
By substituting the relation into the requirements on the velocity vave , we can get
the requirements on the lane number I(x). From the perspective of comfort, there is
a variety of functions I(x) that could provide smooth velocity changes that satisfy
the comfort requirements in Equations (6.35) and (6.36). From the perspective of
efficiency, we would like to find a function I(x) that has the shortest span s. Thus,
149
the road design becomes an optimization question:
0
(v
ave (x) + ω)vave (x) ≤ |b|
min s(I(x)) given (6.41)
0
vave (x)vave
(x) ≤ a
This method can be generalized to other velocity gaps. The key consideration of
comfort is enforced through the comfortable acceleration a and deceleration b.
reshaped fundamental diagram
The goal of reshaped fundamental diagram is to narrow the velocity gap vL /vA .
However, the speed of the shock wave will increase when the gap is narrowed. In
Figure 6.25, we show the contrary relation between vA0 and ω 0 :
total
L I(x) = 1
I(x) = 2
↵1
B A0 A
O ↵2 ⇢total
Figure 6.25: The velocity gap and shock wave speed in reshaped fundamental diagram. A0 is the
intersection of the alternative fundamental diagram with line BA.
The velocity at state A0 is equivalent to the slope of line OA, and the speed
of shock wave is equivalent to the slope of line LA0 . It can be observed that: the
velocity gap vL − vA0 is narrowed when A0 moves aways from state A; in contrast, the
150
shock wave speed ω 0 increases when A0 moves aways from state A. How to choose
an alternative fundamental diagram in this controlling strategy depends on the op-
timization question we are interested. Here I only provide two potential objective
functions:
• Usually drivers prefer not to vary velocity dramatically, so the smaller gap
vL − vA0 is preferred.
• As mentioned in the virtual lane usage, another objective could be the com-
fortable experience is acceleration and deceleration.
From the practical perspective, the objective function could be a combination of
both velocity gap and the span. The weights w1 and w2 are used to balance the two
sub-objectives. Like the virtual lane usage, we prefer to have smaller span s. Unlike
the analysis earlier, we should use the reshaped fundamental diagram to derive the
relation between vave and I(x). The optimization question becomes:
(vave (x) + ω 0 )vave
0
(x) ≤ |b|
min w1 (vL − vA0 ) + w2 s(I(x)) given (6.42)
0
vave (x)vave
(x) ≤ a
6.4 Discussion
In this chapter, we aim to study the lane-drop bottleneck, and provide relevant
controlling strategies to maximize capacity. In the introduction, we described the
lane-drop bottleneck, and the capacity loss phenomena. Observations show that
the queues formed by the lane-changing behavior will bring capacity loss in the
downstream of lane drop.
151
In § 6.2, we introduced a macroscopic approach to model the multi-lane traffic,
which uses only one continuity equation (multi-lane LWR model). Then we analyzed
the solutions to the multi-lane LWR model by applying the supply–demand theory.
The study of traffic pattern demonstrated that the multi-lane LWR model is effi-
cient in organizing the traffic density and maximizing the road capacity. Therefore,
we attempted to modify the multi-lane LWR equation, such that the vehicles can
comfortably follow the velocity solution to the modified equation.
We proposed two controlling strategies in § 6.3 to overcome the velocity gap in
the original velocity solution. The strategy of virtual lane usage represents the idea
of early merge when approaching the lane drop. It aims to smooth the velocity gap
in order to generate comfortable acceleration and deceleration for the vehicles. The
lane usage I should gradually decrease in the upstream of the lane drop under the
comfort requirements. The strategy of reshaped fundamental diagram represents the
idea of slow-down driving when approaching the lane drop. It attempts to narrow
the velocity gap by sacrificing the speed of the backward shock wave. The slope of
the congested flow in the fundamental diagram should be appropriately reshaped to
narrow the velocity gap.
The biggest advantage of the proposed traffic control is that it efficiently max-
imizes the capacity through controlling the macroscopic flow using the multi-lane
LWR model. The two control strategies represent the common senses in lane drop,
which is easy to understand and quantify. The design of the target multi-lane LWR
system can be concluded as an optimization question, and the objective function is
a modeling choice.
The proposed controlling strategies can applied to artificial intelligence of traffic.
In the future, the vehicles will be self-driving or they could communicate with a
152
central control system. Under this circumstance, their movements can obey the
velocity solution to the multi-lane LWR system. That is, the control in traffic flow
is reflected on the movement of individual vehicle.
Chapter Seven
Conclusion
154
In the project of traffic-flow models, I have a systematic study of both macroscopic
and microscopic traffic models as well as their dynamics and applications. In this
thesis, we mainly investigate three interesting traffic-flow problems: the collision be-
havior of car-following models, traffic estimation using data assimilation techniques,
and the traffic control in lane-drop bottleneck. The discoveries and results are sum-
marized as follows:
Collision behavior of car-following models
First, we studied the collision behavior of four well-known car-following models:
the optimal velocity model, the full velocity difference model, the modified GHR
model and the intelligent driver model. The parameter setting is based on the model
calibration using historical traffic data.
• We applied the phase portrait techniques to showing that the optimal velocity
model and the full velocity difference model are collision-prone. This collision-
prone behavior is independent of parameter setting.
• We also applied the fast-slow system techniques and mathematical induction
to proving that the modified GHR model and the intelligent driver model are
collision-free in any traffic scenarios.
• The simulations results are in agreement with the theoretical analysis. How-
ever, collisions can be occasionally observed from the trajectories generated by
the collision-free models: the modified GHR model and the intelligent driver
model.
• More simulations are implemented to show that the numerical errors introduce
155
collisions that the model does not support. In addition, the number of observed
collisions is proportional to the truncation error of the numerical scheme.
In the study of the collision behavior of car-following model, we identified the driver
to prevent collisions in critical situation, which can be taken into consideration when
designing new acceleration models in the future. In addition, the collision-free model
can be used in artificial intelligence, in which the self-driving vehicles can move
according to the collision-free acceleration models.
In the future work, we would like to study the collision behavior of lane-changing
models so that we can combine these safe car-following model and lane-changing to
construct a complete traffic system.
Traffic estimation using data assimilation
Second, we investigated the data assimilation techniques which can be used to esti-
mate the traffic states and uncertain parameters, including ensemble Kalman filter
and particle filter. The initial motivation was to propose an efficient approach that
could assimilate both Eulerian and Lagrangian GPS data simultaneously.
• We proposed an alternative approach that allows us to assimilate both Eulerian
and Lagrangian GPS data simultaneously. The idea behind it is to append the
differential equations for the positions and velocities of the vehicles to the
macroscopic traffic model in order to solve them simultaneously.
• For ensemble Kalman filter and particle filter, we used custom localization,
inflation, weight collision diagnose and resampling methods, which are suitable
for traffic flow estimation.
156
• We added additional terms to the LWR model to capture the effects of on-
/off-ramps, traffic lights, road construction and traveling bottlenecks. The
proposed approach is accurate and works well in different traffic scenarios. In
addition, the proposed approach can improve the accuracy when applied to
realistic traffic data from Minnesota Department of Transportation.
• We compared the performance of ensemble Kalman filter and particle filter.
The proposed approach is accurate and works well regardless of which ensemble
Kalman filter or particle filter is used. Compared to ensemble Kalman filter,
the particle filter is less sensitive to observation noise and sensor locations, but
its computation cost is higher.
One limitation we reported is that the accuracy of data assimilation can be reduced
significantly if the underlying traffic model cannot reproduce the actual traffic flow
from which observation are collected. Hence, more accurate traffic models and cali-
bration are needed to improve the estimation accuracy.
There are various extension that could be implemented in the future work. An
exciting extension would be to apply the proposed data assimilation approach to
more complicated and realistic road network. Second extension would be to use
parameter estimation to relay traffic information to drivers. For instance, data as-
similation can be used to predict the location and duration of congestion caused by
slow trucks or road constructions.
Traffic control in lane-drop bottleneck
We also studied the lane-drop bottleneck, in which there is capacity drop when
queues form in the upstream of the lane drop. Capacity drop brings unnecessary
157
congestion, and reduce the usages of roads. Therefore, we would like to propose
some controlling strategies in order to maximize the capacity.
• We introduced only one continuity equation (inhomogeneous LWR model) to
model the multi-lane traffic. This macroscopic approach focuses on the aggre-
gated traffic characteristics instead of the interactions between lanes.
• The multi-lane LWR model was shown efficient in organizing the traffic density,
and thus maximize the traffic flow in the lane drop. Therefore, we would like
to have vehicles adjust their speeds based on the velocity solution to the multi-
lane LWR model.
• We applied two controlling strategies to adjust the velocity curve such that it
is comfortable for vehicles to follow. The first strategy targeted to smooth the
velocity gaps in the velocity solution through virtual lane usage. The second
strategy attempted to narrow the velocity gaps in the velocity solution through
reshaped fundamental diagram.
In the future work, we would like to combine the proposed controlling strategies
with zipper merging to build an intelligent control system. The acceleration and
lane-changing behavior of individual vehicles are communicated through this system
in order to reduce traffic delay in lane-drop bottleneck. In addition, we look forward
to applying this approach to other traffic bottlenecks, for instance ramp bottleneck
and accident bottleneck.
Bibliography
[1] Administration, F.H.. Over-roadway sensor technologies.
http://www.fhwa.dot.gov/policyinformation/pubs/vdstits2007/05pt2.cfm;
2017.
[2] Anderson, J.L.. An adaptive covariance inflation error correction algorithm for
ensemble filters. Tellus 2007;59A:210–224.
[3] Anderson, J.L.. Spatially and temporally varying adaptive covariance inflation
for ensemble filters. Tellus A 2009;61(1):72–83.
[4] Anderson, J.L., Anderson, S.L.. A Monte Carlo implementation of the non-
linear filtering problem to produce ensemble assimilations and forecasts. Mon
Wea Rev 1999;127(12):2741–2758.
[5] Aw, A., Klar, A., Rascle, M., Materne, T.. Derivation of continuum traffic flow
models from microscopic follow-the-leader models. SIAM Journal on Applied
Mathematics 2002;63(1):259–278.
[6] Baek, S.J., Hunt, B.R., Kalnay, E., Ott, E.. Local ensemble Kalman filtering
in the presence of model bias. Tellus 2006;58A(3):293–306.
[7] Bando, M., Hasebe, K., Nakanishi, K., Nakayama, A.. Analysis of optimal
velocity model with explicit delay. Phys Rev E 1998;58(5):5429.
[8] Bando, M., Hasebe, K., Nakayama, A., Shibata, A., Sugiyama, Y.. Dynam-
ical model of traffic congestion and numerical simulation. Physical Review E
1995;51(2):1035.
[9] Bellsky, T., Berwald, J., Mitchell, L.. Nonglobal parameter estimation using
local ensemble Kalman filtering. Mon Wea Rev 2014;142(6):2150–2164.
[10] Bertini, R.L., Leal, M.T.. Empirical study of traffic features at a freeway lane
drop. Journal of Transportation Engineering 2005;131(6):397–407.
[11] Bevrani, K., Chung, E., Miska, M.. Evaluation of the ghr car following model
for traffic safety studies. In: Proceedings of the 25th ARRB Conference. ARRB
Group Ltd; 2012. p. 1–11.
158
159
[12] Biswas, S., Tatchikou, R., Dion, F.. Vehicle-to-vehicle wireless communication
protocols for enhancing highway traffic safety. IEEE communications magazine
2006;44(1):74–82.
[13] Bonsall, P., Liu, R., Young, W.. Modelling safety-related driving behaviourim-
pact of parameter values. Transportation Research Part A: Policy and Practice
2005;39(5):425–444.
[14] Brackstone, M., McDonald, M.. Car-following: a historical review. Transporta-
tion Research Part F: Traffic Psychology and Behaviour 1999;2(4):181–196.
[15] Burgers, G., van Leeuwen, P.J., Evensen, G.. Analysis scheme in the ensemble
Kalman filter. Mon Wea Rev 1998;126:1719–1724.
[16] Carter, P., Christiansen, P., Gaididei, Y., Gorria, C., Sandstede, B.,
Sorensen, M., Starke, J.. Multi-jam solutions in traffic models with velocity-
dependent driver strategies. SIAM J Appl Math 2014;74(6):1895–1918.
[17] Chandler, R.E., Herman, R., Montroll, E.W.. Traffic dynamics: studies in
car following. Operations research 1958;6(2):165–184.
[18] CHUNG, S.B., SONG, K.H., HONG, S.Y., KHO, S.Y.. Development of
sensitivity term in car-following model considering practical driving behavior of
preventing rear end collision. Journal of the Eastern Asia Society for Trans-
portation Studies 2005;6:1354–1367.
[19] Daganzo, C.F.. The cell transmission model: A dynamic representation of high-
way traffic consistent with the hydrodynamic theory. Transportation Research
Part B: Methodological 1994;28(4):269–287.
[20] Daganzo, C.F.. Requiem for second-order fluid approximations of traffic flow.
Transportation Research Part B: Methodological 1995;29(4):277–286.
[21] Daganzo, C.F.. A behavioral theory of multi-lane traffic flow. part i: Long
homogeneous freeway sections. Transportation Research Part B: Methodological
2002;36(2):131–158.
[22] Daganzo, C.F.. A behavioral theory of multi-lane traffic flow. part ii: Merges
and the onset of congestion. Transportation Research Part B: Methodological
2002;36(2):159–169.
[23] Dee, D.P., Da Silva, A.M.. Data assimilation in the presence of forecast bias.
Q J Roy Meteor Soc 1998;124(545):269–295.
[24] Dervisoglu, G., Gomes, G., Kwon, J., Horowitz, R., Varaiya, P.. Automatic
calibration of the fundamental diagram and empirical observations on capacity.
In: Transportation Research Board 88th Annual Meeting. volume 15; 2009. .
[25] Di Francesco, M., Rosini, M.D.. Rigorous derivation of nonlinear scalar conser-
vation laws from follow-the-leader type models via many particle limit. Archive
for Rational Mechanics and Analysis 2015;217(3):831–871.
160
[26] Doucet, A., Freitas, N.D., Gordon, N.. An introduction to sequential Monte
Carlo methods. In: Sequential Monte Carlo methods in practice. Springer; 2001.
p. 3–14.
[27] Doucet, A., Godsill, S., Andrieu, C.. On sequential Monte Carlo sampling
methods for Bayesian filtering. Stat Comput 2000;10:197–208.
[28] Evensen, G.. Sequential data assimilation with a nonlinear quasi-geostrophic
model using Monte Carlo methods to forecast error statistics. J Geophys Res:
Oceans 1994;99(C5):10143–10162.
[29] Evensen, G.. The ensemble Kalman filter: theoretical formulation and practical
implementation. Ocean dynamics 2003;53(4):343–367.
[30] Gaididei, Y.B., Gorria, C., Berkemer, R., Kawamoto, A., Shiga, T., Chris-
tiansen, P.L., Sørensen, M.P., Starke, J.. Controlling traffic jams by time-
modulating the safety distance. Phys Rev E 2013;.
[31] Gaspari, G., Cohn, S.E.. Construction of correlation functions in two and
three dimensions. Q J Roy Meteor Soc 1999;125(554):723–757.
[32] Gazis, D.C., Herman, R., Rothery, R.W.. Nonlinear follow-the-leader models
of traffic flow. Operations research 1961;9(4):545–567.
[33] Godunov, S.K.. A difference method for numerical calculation of discontin-
uous solutions of the equations of hydrodynamics. Matematicheskii Sbornik
1959;89(3):271–306.
[34] Gordon, N.J., Salmond, D.J., Smith, A.F.M.. Novel approach to nonlinear-
non-Gaussian Bayesian state estimation. IEE Proc-F 1993;140(2):107–113.
[35] Gov, G.. GPS accuracy. http://www.gps.gov/systems/gps/performance/accuracy/.
[36] Greenshields, B.D., Channing, W., Miller, H., et al. A study of traffic
capacity. In: Highway research board proceedings. National Research Council
(USA), Highway Research Board; volume 1935; 1935. .
[37] Hamdar, S., Mahmassani, H.. From existing accident-free car-following mod-
els to colliding vehicles: exploration and assessment. Transportation Research
Record: Journal of the Transportation Research Board 2008;(2088):45–56.
[38] Hamill, T.M., Whitaker, J.S., Snyder, C.. Distance-dependent filtering of
background error covariance estimates in an ensemble Kalman filter. Mon Wea
Rev 2001;129(11):2776–2790.
[39] Helbing, D.. Traffic and related self-driven many-particle systems. Rev of Mod
Phys 2001;73(4):1067.
[40] Henrici, P.. Applied and computational complex analysis, discrete Fourier
analysis, Cauchy integrals, construction of conformal maps, univalent functions.
volume 3. John Wiley & Sons, 1993.
161
[41] Hidas, P.. Modelling lane changing and merging in microscopic traffic simula-
tion. Transportation Research Part C: Emerging Technologies 2002;10(5):351–
371.
[42] Holden, H., Risebro, N.H.. A mathematical model of traffic flow on a network of
unidirectional roads. SIAM Journal on Mathematical Analysis 1995;26(4):999–
1017.
[43] Houtekamer, P.L., Mitchell, H.L.. Data assimilation using an ensemble Kalman
filter technique. Mon Wea Rev 1998;126:796–811.
[44] Houtekamer, P.L., Mitchell, H.L.. A sequential ensemble Kalman filter for
atmospheric data assimilation. Mon Wea Rev 2001;129(1):123–137.
[45] Isaacson, E., Temple, B.. Nonlinear resonance in systems of conservation laws.
SIAM Journal on Applied Mathematics 1992;52(5):1260–1278.
[46] Jiang, R., Wu, Q., Zhu, Z.. Full velocity difference model for a car-following
theory. Physical Review E 2001;64(1):017101.
[47] Jin, W.L., Chen, L., Puckett, E.G.. Supply-demand diagrams and a new
framework for analyzing the inhomogeneous lighthill-whitham-richards model.
In: Transportation and Traffic Theory 2009: Golden Jubilee. Springer; 2009. p.
603–635.
[48] Kato, S., Tsugawa, S., Tokuda, K., Matsui, T., Fujii, H.. Vehicle con-
trol algorithms for cooperative driving with automated vehicles and intervehi-
cle communications. IEEE Transactions on Intelligent Transportation Systems
2002;3(3):155–161.
[49] Kerner, B.S.. The physics of traffic: empirical freeway pattern features, engi-
neering applications, and theory. Springer, 2012.
[50] Kesting, A., Treiber, M.. Calibrating car-following models by using trajectory
data: Methodological study. Transportation Research Record: Journal of the
Transportation Research Board 2008;(2088):148–156.
[51] Kong, A., Liu, J., Wong, W.. Sequential imputations and Bayesian missing
data problems. J Am Stat Assoc 1994;89:278–288.
[52] Laval, J., Cassidy, M., Daganzo, C.. Impacts of lane changes at merge
bottlenecks: a theory and strategies to maximize capacity. In: Traffic and
Granular Flow05. Springer; 2007. p. 577–586.
[53] Laval, J.A., Daganzo, C.F.. Lane-changing in traffic streams. Transportation
Research Part B: Methodological 2006;40(3):251–264.
[54] Lighthill, M.J., Whitham, G.B.. 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 1955;229(1178):317–345.
[55] Liu, J., West, M.. Combined parameter and state estimation in simulation-
based filtering. In: Sequential Monte Carlo methods in practice. Springer; 2001.
p. 197–223.
162
[56] Mandel., J.. Efficient implementation of the ensemble Kalman filter. University
of Colorado at Denver and Health Sciences Center, Center for Computational
Mathematics, 2006.
[57] Nagel, K., Schreckenberg, M.. A cellular automaton model for freeway traffic.
Journal de physique I 1992;2(12):2221–2229.
[58] Net, G.I.. GPS speed. http://www.gpsinformation.net/main/gpsspeed.htm.
[59] Panwai, S., Dia, H.. Comparative evaluation of microscopic car-
following behavior. Intelligent Transportation Systems, IEEE Transactions on
2005;6(3):314–325.
[60] Payne, H.J.. Models of freeway traffic and control. Mathematical models of
public systems 1971;.
[61] Rahman, M., Chowdhury, M., Xie, Y., He, Y.. Review of microscopic
lane-changing models and future research opportunities. IEEE transactions on
intelligent transportation systems 2013;14(4):1942–1956.
[62] Replay, G.A.. GPS and speed measurement.
http://gpsactionreplay.free.fr/index.php?menu=6&choice=5.
[63] Richards, P.I.. Shock waves on the highway. Oper Res 1956;4(1):42–51.
[64] Rosales, R.R., Kasimov, A.R., Flynn, M.R., Seibold, B.. Constructing
set-valued fundamental diagrams from jamiton solutions in second order traffic
models 2013;.
[65] Santitissadeekorn, N., Jones, C.K.R.T.. Two-state filtering for joint state-
parameter estimation. arXiv preprint arXiv:14035989 2014;.
[66] Snyder, C.. In: Proceedings of the ECMWF Seminar on Data Assimilation for
Atmosphere and Ocean. 2011. .
[67] Snyder, C., Bengtsson, T., Bickel, P., Anderson, J.. Obstacles to high-
dimensional particle filtering. Mon Wea Rev 2008;136(12):4629–4640.
[68] Snyder, C., Bengtsson, T., Morzfeld, M.. Performance bounds for particle
filters using the optimal proposal. Mon Wea Rev 2015;143(11):4750.
[69] of Transportation, C.D.. The freeway performance measurement (pems).
http://pems.dot.ca.gov/; .
[70] of Transportation, M.D.. Traffic forecasting & analysis.
http://www.dot.state.mn.us/traffic/data/; .
[71] Trefethen, L.N.. Spectral methods in MATLAB. SIAM, 2000.
[72] Treiber, M., Hennecke, A., Helbing, D.. Congested traffic states in empirical
observations and microscopic simulations. Physical Review E 2000;62(2):1805.
[73] Treiber, M., Kesting, A.. Traffic flow dynamics. Traffic Flow Dynamics: Data,
Models and Simulation, Springer-Verlag Berlin Heidelberg 2013;.
163
[74] van Leeuwen, P.J.. Particle filtering in geophysical systems. Mon Wea Rev
2009;137:4089–4114.
[75] Vogel, K.. A comparison of headway and time to collision as safety indicators.
Accident analysis & prevention 2003;35(3):427–433.
[76] Work, D.B., Blandin, S., Tossavainen, O.P., Piccoli, B., Bayen, A.M.. A traf-
fic model for velocity data assimilation. Appl Math Res Express 2010;2010(1):1–
35.
[77] Work, D.B., Tossavainen, O.P., Blandin, S., Bayen, A.M., Iwuchukwu,
T., Tracton, K.. An ensemble Kalman filtering approach to highway traffic
estimation using GPS enabled mobile devices. In: 2008 47th IEEE Decis. Contr.
P. IEEE; 2008. p. 5062–5068.
[78] Yang, X., Delsole, T.. Using the ensemble Kalman filter to estimate multi-
plicative model parameters. Tellus 2009;61A(5):601–609.
[79] Yuan, Y., Van Lint, H., Van Wageningen-Kessels, F., Hoogendoorn, S..
Network-wide traffic state estimation using loop detector and floating car data.
J Intell Transport S: Technology, Planning, and Operations 2014;18(1).
[80] Yuan, Y., Van Lint, J.W.C., Wilson, R.E., Van Wageningen-Kessels, F.,
Hoogendoorn, S.P.. Real-time Lagrangian traffic state estimator for freeways.
IEEE Trans Intell Transp Syst 2012;13(1):59–70.
[81] Zhang, H.M.. Driver memory, traffic viscosity and a viscous vehicular traffic
flow model. Transport Res B-Meth 2003;37(1):27–41.