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.