Fast Calibration of Car Following models to Trajectory data using the Adjoint Method

TL;DR

This paper introduces an efficient adjoint method for calibrating car-following models to trajectory data, significantly reducing computational cost and improving accuracy compared to traditional gradient-free optimization techniques.

Contribution

It develops an analytical gradient calculation using the adjoint method for car-following model calibration, enhancing efficiency and accuracy over existing methods.

Findings

Adjoint method reduces calibration computational cost.

Quasi-Newton algorithms outperform genetic algorithms.

Calibration accuracy is improved with the adjoint approach.

Abstract

Before a car-following model can be applied in practice, it must first be validated against real data in a process known as calibration. This paper discusses the formulation of calibration as an optimization problem, and compares different algorithms for its solution. The optimization consists of an arbitrary car following model, posed as either an ordinary or delay differential equation, being calibrated to an arbitrary source of trajectory data which may include lane changes. Typically, the calibration problem is solved using gradient free optimization. In this work, the gradient of the optimization problem is derived analytically using the adjoint method. The computational cost of the adjoint method does not scale with the number of model parameters, which makes it more efficient than evaluating the gradient numerically using finite differences. Numerical results are presented which show that quasi-newton algorithms using the adjoint method are significantly faster than a genetic algorithm, and also achieve slightly better accuracy of the calibrated model.

Tables

Table 1. Table 1: Table of speeds for adjoint method compared to finite differences for computing the gradient. Relative time refers to the equivalent number of objective function evaluations.

Parameters	Obj. time (s)	Adjoint time (s)	Adjoint relative time	Finite time (s)	Finite relative time
5	.0062	.0250	4.03	.0390	6.3
10	.0109	.0438	4.02	.1250	11.5
15	.0172	.0686	3.99	.2811	16.3

Table 2. Table 2: Algorithms on the Pareto front. All algorithms tested shown in appendix.

Algorithm	% found global opt	Avg. RMSE (ft)	Avg. Time (s)	Avg. RMSE / % over global opt	Initial Guesses	# Obj. Evals	# Grad. Evals
Adj BFGS-0	94.0	6.46	7.42	.10 / 2.0%	3	307.1	307.1
Adj BFGS-7.5	85.2	6.56	4.16	.20/6.8%	1.67	165.0	165.0
GA	95.0	6.47	33.5	.10/2.2%	-	5391	0
Fin BFGS-0	94.3	6.46	11.4	.10/2.2%	3	311.5	311.5
Fin BFGS-7.5	85.7	6.55	6.32	.19/6.3%	1.66	164.8	164.8
Adj TNC-0	90.7	6.43	6.87	.05/2.4%	3	285.7	285.7
Adj TNC-7.5	82.7	6.48	3.82	.11/5.0%	1.63	152.3	152.3
Adj TNC-$\infty$	77.3	6.62	2.26	.25/7.7%	1	94.1	94.1

Table 3. Table 3: Algorithms on the Pareto front depend on how “finding the global optimum” is defined (without that metric, TNC with the adjoint method would be the only algorithm on the pareto front). Set notation indicates the strategies used for initial guesses.

Tolerance for global opt (ft)	Algorithms on the pareto front
0	Adj TNC-{0, 7.5, Inf}, Adj BFGS-{0, 7.5, Inf}
1/24	Adj TNC-{0, 7.5, Inf}, Adj BFGS-{0, 7.5, Inf}, Fin BFGS-{0, 7.5}
1/12	Adj TNC-{0, 7.5, Inf}, Adj BFGS-{0, 7.5}, Fin BFGS-{0, 7.5}, GA
1/4	Adj TNC-{0, 7.5, Inf}, GA
1/2	Adj TNC-{0, 7.5, Inf}

Table 4. Table 4: Compares the number of equivalent objective evaluations for the calibration problem for an increasing number of parameters, relative to the Adj TNC-0 algorithm.

# parameters	Adj TNC	Fin TNC	Fin TNC relative evals.	GA	GA relative evals.
5	1098	1790	1.63	5256	4.79
10	1200	3300	2.75	37715	31.4
15	1780	7139	4.01	135485	76.1

Table 5. Table 5: The full table showing all variants of the 7 algorithms considered.

Algorithm	% found global opt	Avg. RMSE (ft)	Avg. Time (s)	Avg. RMSE / % over global opt	Initial Guesses	# Obj.	# Grad.	# Hess.
Adj BFGS-0	94.0	6.46	7.42	.10 / 2.0%	3	307.1	307.1	0
Adj BFGS-7.5	85.2	6.56	4.16	.20/6.8%	1.67	165.0	165.0	0
Adj BFGS-$\infty$	75.6	6.96	2.55	.59 / 15.4%	1	107.3	107.3	0
GA	95.0	6.47	33.5	.10/2.2%	-	5391	0	0
NM	39.4	7.25	6.85	.88/21.0%	1	1037	0	0
Fin BFGS-0	94.3	6.46	11.4	.10/2.2%	3	311.5	311.5	0
Fin BFGS-7.5	85.7	6.55	6.32	.19/6.3%	1.66	164.8	164.8	0
Fin BFGS-$\infty$	77.0	6.88	3.91	.52/13.8%	1	107.3	107.3	0
Adj GD-0	20.2	7.58	25.5	1.21/28.2%	3	1374	542.5	0
Adj GD-7.5	13.4	7.83	18.1	1.46/40.8%	1.93	1094	436.6	0
Adj GD-Inf	8.8	8.90	8.75	2.53/67.2%	1	719.3	282.6	0
Adj SQP-0	29.1	8.37	21.1	2.01/61.7%	3	287.6	83.8	80.8
Adj SQP-7.5	21.5	8.59	15.4	2.22/71.4%	2.17	244.8	71.3	69.1
Adj SQP-$\infty$	10.0	11.1	5.87	4.70/140.7%	1	141.6	39.4	38.4
SPSA	4.2	53.2	38.5	46.8/786%	1	6001	0	0
Fin GD-$\infty$	9.1	8.86	14.7	2.50/66.2%	1	657.3	258.5	0
Fin SQP-$\infty$	2.6	9.61	15.8	3.25/81.8%	1	179.4	67.2	66.2
Adj TNC-0	90.7	6.43	6.87	.05/2.4%	3	285.7	285.7	0
Adj TNC-7.5	82.7	6.48	3.82	.11/5.0%	1.63	152.3	152.3	0
Adj TNC-$\infty$	77.3	6.62	2.26	.25/7.7%	1	94.1	94.1	0
Fin TNC-0	36.3	6.84	13.2	.47/11.0%	3	302.4	302.4	0
Fin TNC-7.5	24.9	7.05	8.35	.69/20.8%	1.82	183.2	183.2	0
Fin TNC-$\infty$	18.0	7.63	4.49	1.27/33.2%	1	100.8	100.8	0

Equations

$${\ddot{x}}_n(t)=h({\dot{x}}_{n-1}(t-\tau ),x_{n-1}(t-\tau ),{\dot{x}}_n(t-\tau ),x_n(t-\tau ),p)$$

$${\dot{x}}_n(t)=\min \{v_f,\frac{x_{n-1}(t)-x_n(t)-s_{\text{jam}}}{\alpha }\}.$$

$$x_n(t)=x_{n-1}(t-\alpha )-s_{\text{jam}}$$

$${\dot{x}}_n(t)=\left(\begin{array}{c}{\dot{x}}_n^{\ast }(t)\\ {\dot{v}}_n(t)\end{array}\right)=\left(\begin{array}{c}{v}_n(t)\\ h^{\ast }(x_{n-1}(t-{\tau }_n),x_n(t-{\tau }_n),p_n)\end{array}\right)=h_n(x_n(t),x_n(t-{\tau }_n),x_{L(n)}(t-{\tau }_n),p_n)$$

$$\underset{p}{\min }$$

$${\dot{x}}_i(t)-h_i(x_i(t),x_{L(i)}(t),p_i)=0,t\in [t_i,T_{i-1}],i=1,\dots ,n$$

$$x_i(t_i)-{\hat{x}}_i(t_i)=0i=1,\dots ,n$$

$$b_{\mathrm{low}}\le p\le b_{\mathrm{high}}$$

$$\underset{p}{\min }$$

$${\dot{x}}_i(t)-h_i(x_i(t),x_i(t-{\tau }_i),x_{L(i)}(t-{\tau }_i),p_i)=0,t\in [t_i+{\tau }_i,T_{i-1}^{\ast }],i=1,\dots ,n$$

$$x_i(t)-{\hat{x}}_i(t)=0t\in [t_i,t_i+{\tau }_i],i=1,\dots ,n$$

$$b_{\mathrm{low}}\le p\le b_{\mathrm{high}}$$

$${\dot{x}}_i(t)={\dot{\hat{x}}}_i(t),t\in [T_{i-1},T_i]$$

$$h_i(x_i,x_{L(i)},p_i)=\{\begin{array}{ll}{h}_{i,2}& g_{i,2}(x_i,x_{L(i)},p_i)\ge 0\\ h_{i,1}& \text{ otherwise}\end{array}$$

$$\frac{\partial x_i}{\partial p}=\left(\begin{array}{cccc}\partial x_i^{\ast }/\partial p^1& \partial x_i^{\ast }/\partial p^2& \dots & \partial x_i^{\ast }/\partial p^m\\ \partial v_i/\partial p^1& \partial v_i/\partial p^2& \dots & \partial v_i/\partial p^m\end{array}\right)$$

$$L=\int_{t_{n}}^{T_{n-1}}dt\bigg{[}f(x_{n},\hat{x}_{n},t)+\lambda^{T}(t)\big{(}\dot{x}_{n}(t)-h_{n}(x_{n}(t),x_{L(n)}(t),p_{n})\big{)}\bigg{]}$$

$$\dfrac{d}{dp}F=\dfrac{d}{dp}L=\int_{t_{n}}^{T_{n-1}}dt\bigg{[}\dfrac{\partial f}{\partial x_{n}}\dfrac{\partial x_{n}}{\partial p}+\lambda^{T}(t)\bigg{(}\dfrac{\partial\dot{x}_{n}}{\partial p}-\dfrac{\partial h_{n}}{\partial x_{n}}\dfrac{\partial x_{n}}{\partial p}-\dfrac{\partial h_{n}}{\partial p}\bigg{)}\bigg{]}$$

$$\int_{t_{n}}^{T_{n-1}}\lambda^{T}(t)\dfrac{\partial\dot{x}_{n}}{\partial p}dt=\lambda^{T}(t)\dfrac{\partial x_{n}}{\partial p}\bigg{|}^{T_{n-1}}_{t_{n}}-\int_{t_{n}}^{T_{n-1}}\dot{\lambda}^{T}(t)\dfrac{\partial x_{n}}{\partial p}dt$$

$$\dfrac{d}{dp}F=\int_{t_{n}}^{T_{n-1}}dt\bigg{[}\left(\dfrac{\partial f}{\partial x_{n}}-\dot{\lambda}^{T}(t)-\lambda^{T}(t)\dfrac{\partial h_{n}}{\partial x_{n}}\right)\dfrac{\partial x_{n}}{\partial p}-\lambda^{T}(t)\dfrac{\partial h_{n}}{\partial p}\bigg{]}$$

$$\frac{\partial f}{\partial x_n}-{\dot{\lambda }}^T(t)-{\lambda }^T(t)\frac{\partial h_n}{\partial x_n}=0,t\in [T_{n-1},t_n]$$

$$\frac{d}{dp}F=-{\int }_{t_n}^{T_{n-1}}{\lambda }^T(t)\frac{\partial h_n}{\partial p}dt$$

$$L=\sum_{i=1}^{n}\left(\int_{t_{i}}^{T_{i-1}}f(x_{i},\hat{x}_{i})+\lambda_{i}^{T}\big{(}\dot{x}_{i}-h_{i}(x_{i},x_{L(i)},p_{i})\big{)}dt+\int_{T_{i-1}}^{T_{i}}\lambda_{i}^{T}(\dot{x}_{i}-\dot{\hat{x}}_{i})dt\right)$$

$$\frac{d}{dp}F=\frac{d}{dp}L=\sum _{i=1}^n\left({\int }_{t_i}^{T_{i-1}}\frac{\partial f}{\partial x_i}\frac{\partial x_i}{\partial p}-{\lambda }_i^T(t)\frac{\partial h_i}{\partial x_i}\frac{\partial x_i}{\partial p}-{\lambda }_i^T(t)\frac{\partial h_i}{\partial p}dt\right)$$

$$+\sum _{i=1}^n\left({\int }_{t_i}^{T_i}{\lambda }_i^T(t)\frac{\partial {\dot{x}}_i}{\partial p}-\text{mathbbm}1(G(i,t)\ne 0){\lambda }_{G(i)}^T(t)\frac{\partial h_{G(i)}}{\partial x_i}\frac{\partial x_i}{\partial p}dt\right)$$

$$-{\dot{\lambda }}_i^T(t)+\text{mathbbm}1(t\le T_{i-1})\left(\frac{\partial f}{\partial x_i}-{\lambda }_i^T(t)\frac{\partial h_i}{\partial x_i}\right)-\text{mathbbm}1(G(i,t)\ne 0){\lambda }_{G(i)}^T(t)\frac{\partial h_{G(i)}}{\partial x_i}=0,t\in [T_i,t_i],\forall i$$

$$\frac{d}{dp}F=\sum _{i=1}^n\left(-{\int }_{t_i}^{T_{i-1}}{\lambda }_i^T(t)\frac{\partial h_i}{\partial p}dt\right)$$

$${\ddot{x}}_n(t)=c_4(V(s)-{\dot{x}}_n(t))$$

$$V(s)=c_1[\tanh (c_{2}s-c_3-c_5)-\tanh (-c_3)]$$

$$f(x_i,{\hat{x}}_i,t)=(x_i^{\ast }(t)-{\hat{x}}_i^{\ast }(t){)}^2$$

$$F=\sum _{i=1}^n\sum _{t=t_i}^{T_{i-1}}f(x_i,{\hat{x}}_i,t)$$

$$\text{RMSE}=\sqrt{\frac{F}{{\sum }_{i=1}^n{T}_{i-1}-t_i}}$$

$$\text{relative error}=\frac{||{\nabla }_{fd}F-\nabla F|{|}_2}{||{\nabla }_{fd}F|{|}_2}$$

$$\dot{x}(t)=h(t,x(t)),x(t_0)=x_0$$

$$\left|\frac{\partial h}{\partial x}\right|\le C$$

$$\underset{ϵ\to 0}{\lim }\{{\theta }_j(p+ϵ)\}=\{{\theta }_j(p)\}$$

Full text

Fast Calibration of Car Following models to Trajectory data using the Adjoint Method

Ronan Keane

Systems Engineering, Cornell University, Email: (rlk268@cornell.edu)

H. Oliver Gao

Civil and Environmental Engineering, Cornell University, Email: (hg55@cornell.edu)

Abstract

Before a car-following model can be applied in practice, it must first be validated against real data in a process known as calibration. This paper discusses the formulation of calibration as an optimization problem, and compares different algorithms for its solution. The optimization consists of an arbitrary car following model, posed as either an ordinary or delay differential equation, being calibrated to an arbitrary source of trajectory data which may include lane changes.

Typically, the calibration problem is solved using gradient free optimization. In this work, the gradient of the optimization problem is derived analytically using the adjoint method. The computational cost of the adjoint method does not scale with the number of model parameters, which makes it more efficient than evaluating the gradient numerically using finite differences. Numerical results are presented which show that quasi-newton algorithms using the adjoint method are significantly faster than a genetic algorithm, and also achieve slightly better accuracy of the calibrated model.

1 Introduction

1.1 The Case for Car Following and Trajectory Data

A complete traffic simulator consists of many different modules, including algorithms to determine the demand for different road sections, complex rules to govern lane changes, mergers, and diverges, a model to describe traffic flow and driving behaviors, and more. Out of all these modules, the description of traffic flow is the most fundamental. The decision of what model is used to describe traffic flow has impacts on the simulator’s strengths, weaknesses, and overall performance. In particular, one important distinction is whether to simulate at the macroscopic or microscopic level.

At the microscopic level, car following models are often used. Car following models describe how a vehicle behaves as a function of its surroundings; typically they are expressed as differential equations where any vehicle is coupled to the vehicle directly ahead (i.e. its lead vehicle). Many car following models are based on common sense rules of driving. Some simple examples are adopting a speed based on the distance between you and your leader, or accelerating based on the difference in speeds between you and your leader.

Whereas microscopic traffic models describe the motions of individual vehicles, their macroscopic counterparts treat traffic as a fluid. At the macroscopic level, there are no vehicles, only traffic densities. Congested road sections with high densities slowly trickle forward, while uncongested traffic at low densities flow freely. This connection to fluid dynamics was originally made by [1], and the influence of the resulting traffic flow theory can hardly be overstated. In the more than six decades since, numerous works have attempted to improve upon the so-called Lighthill-Whitham-Richards (LWR) model in order to better describe various complex elements of traffic flow. Among those, [2], [3], [4], [5], [6] are perhaps some of the most notable. Nonetheless, when it comes to macrosimulation, the LWR model, and its discrete numerical formulation ([7]), is still widely used today.

For all its merits, the LWR model is not without its limitations. A key relationship between flow and density known as the “fundamental diagram”, must be supplied. More importantly, it is assumed vehicles are homogeneous and exhibit static behavior. In real traffic, drivers behave differently; vehicles have varying performances; traffic demand fluctuates; controlled intersections and lane changes have complex effects. Traffic flow is incredibly dynamic and inhomogeneous. The shortcomings of the LWR model become apparent in practice (see [4, 8, 9]). In the fundamental diagram, instead of equilibrium curves, one finds a wide scattering of points. The time evolution of real traffic flow shows a much richer picture than what the LWR model predicts. Examining the evolution of a vehicle trajectory in the speed-spacing plane, one sees complex hysteresis loops which are inconsistent with current macroscopic theory ([10], [11] or Fig. 2 are examples of such plots). Given these limitations, one would hope for some car following model that could explain all these phenomena.

However, existing works are not prepared to offer guidelines as to which car following models are best. It is also not clear as to how much more accurate car following models are compared to macro-models (like LWR). Answering either of these pressing questions is difficult: the traffic microsimulation literature consists of literally hundreds of different models, many of which have never actually been thoroughly validated ([12]). Works such as [13], [14], [15], and [16] have attempted to address this long-standing problem of car following models but have had varying results as to which models are best.

Drawing conclusions on the efficacy of different car following models requires both a large vehicle trajectory dataset (e.g. [17] or [18]) to facilitate comparison, as well as a reliable methodology for model calibration and validation. The aim of this work is to develop new methodology to improve the calibration of car following models. However, one should keep in mind that this only addresses one half of the challenge; equally important is the creation of new datasets and new approaches for collecting accurate and meaningful traffic data.

There is a need for an increased use of trajectory data to study traffic flow. This need goes far beyond the validation of existing car following models. As new technologies like autonomous vehicles, connected cities and on-demand mobility become fully realized, the LWR model assumption of homogeneous drivers exhibiting static behavior becomes increasingly unrealistic. The limitations of macrosimulation and the uncertainty of new technology are both good reasons to suspect a growing demand for microsimulation. At the same time, vehicle trajectory data (which is needed for microsimulation/car following) is easier to get than ever before because of new technologies like computer vision, lidar/radar, and GPS. Compared to the more commonly used loop detector data, trajectory data describes traffic flow at much higher temporal and spatial resolutions. These new sources of data, still largely unexploited, can help advance both micro- and macro- theory.

1.2 Calibration Using Trajectory Data

It is relatively simple to write down some mathematical equations and call them a traffic model. Far more difficult is applying this model in practice. In order to validate a model against real data, one must find the model parameters that best explain the data: this is known as calibration. Solving the calibration problem requires both a dataset describing the quantity of interest, as well as a mathematical model that explains the data. This paper deals with car-following models that are posed as differential equations, and may include elements such as delays, bounds on acceleration/speed, and multiple driving regimes. The car-following model is to be calibrated against trajectory data, time series data that describes the motion of individual vehicles.

Existing literature has covered various aspects of the trajectory data calibration problem: the use of meta-models, sensitivity analysis, optimization algorithms, measures of performance (e.g. speed or spacing data), and goodness-of-fit function (i.e. loss function). [19], [20], and [21] are good works focused primarily on methodology. [22], [23], [15] and [16] are good works focused on both calibration and the benchmarking of various models.

In this paper calibration is posed as an optimization problem. For a wide class of models, being calibrated to trajectory data gathered from any source, the gradient of the optimization problem has been derived in an efficient way using the adjoint method. This enables the calibration problem to be solved by more efficient optimization routines which use gradient information. Several state of the art algorithms are benchmarked in order to determine both the accuracy of the fit they can achieve, as well as the computational effort required to apply the algorithm. The paper is organized as follows. Section 2 introduces car following models and discusses current research directions for improving car following models. Section 3 describes the optimization problem and discusses some issues arising in its solution. Section 4 applies the adjoint method to the optimization problem and compares it to finite differences. Section 5 compares different algorithms for solving the calibration problem.

2 Car-Following Models

We will consider any car-following model that can be expressed as either an ordinary or delay differential equation. This describes the majority of models in the literature, as well as the majority of models used in commercial microsimulation packages ([24]). We deal with models of the following form

[TABLE]

where $x_{n}(t)$ is the position of vehicle $n$ at time $t$. $h$ is the car-following model, which depends on parameters $p$, the “following” vehicle $x_{n}$ and “lead” vehicle $x_{n-1}$. $\tau$ is a time-delay parameter which represents human reaction time. Eq. (2) states that a driver’s behavior (how they choose to accelerate) is determined from the speed and position of their own car, the speed and position of the leading car, and parameters $p$. Given a lead vehicle trajectory $x_{n-1}(t)$, parameters $p$, and initial condition $x_{n}(t_{n})$, the car following model $\eqref{1}$ will produce the follower trajectory $x_{n}(t)$. The presence of nonzero $\tau$ leads to a delay differential equation (DDE). Many car-following models (including most of those used in commercial software) treat reaction time explicitly by using a DDE.

Besides ODEs and DDEs, some car-following models are formulated as stochastic differential equations or derived using nonparametric techniques. We briefly discuss these types of models in section 2.3 but the current methodology does not address their calibration.

2.1 Ordinary Differential Equations

Eq. $\eqref{1}$ gives a general form for car-following models. The simplest type of model is a first-order ODE where the speed of the following vehicle $n$ is completely determined by the space headway $(x_{n-1}-x_{n})$. In this type of model, each driver adopts a speed based on the distance between their car and the leading car.

This simple kind of model is important because of its connection to macroscopic theory and the LWR model. Specifically, if one assumes all vehicles are identical, and interprets the reciprocal of density as headway, then a functional relationship between traffic flow and density can be recast as a relationship between vehicle speed and headway ([25]). Under these assumptions, the fundamental diagram can be thought of as describing a first-order car-following model. Consider the Newell car following model an example

[TABLE]

$v_{f},\alpha$ and $s_{\text{jam}}$ are parameters describing a triangular fundamental diagram ([22]). Namely $v_{f}$ is the speed of a vehicle during free flow, $s_{\text{jam}}$ is the distance between two completely stopped vehicles, and $-s_{\text{jam}}/\alpha$ is the shockwave speed.

The Newell model can also be written as a trajectory translation model, where the trajectory of a following vehicle is equivalent to the lead vehicle trajectory, translated over both space and time.

[TABLE]

We refer to Eq. (2.1) as Newell (TT) and Eq. (2.1) as Newell (DE), where the acronyms stand for trajectory translation and differential equation, respectively. In the special case where the TT and DE models both use a numerical time step equal to $\alpha$, and vehicles are always in the congested regime (i.e. not moving at their free flow speeds), the two models become equivalent ([25]). In this special case, either model gives trajectories that exactly obey the LWR model with a triangular fundamental diagram. In general however, the two models are different (see for example Fig. 2). The advantages of the Newell (TT) or (DE) models are that they contain a small number of parameters, and those parameters are easily interpretable under macroscopic theory.

More complicated car-following models (compared to those of Newell) are sometimes criticized for two main reasons: first, for not having enough empirical support for their validity, and second, for having too many parameters. Despite these criticisms, more complicated models have behaviors that attract researchers. For example, the optimal velocity model (OVM) proposed by [26] exhibits a spontaneous phase transition between free flow and congested states, which may describe “phantom” traffic jams ([27]). Another example is the intelligent driver model (IDM) proposed by [28], which can qualitatively reproduce a number of complex (empirical) traffic states. Besides more recent models like OVM or IDM, there are also a plethora of older models to complicate the mix ([12]). The default models in most modern traffic simulator packages (for example, VISSIM or AIMSUN) are still older car-following models ([24]), presumably because it is unclear if the newer car-following models offer any benefits.

Fig. 2 compares each of the four models just mentioned by calibrating them to a single vehicle trajectory from the reconstructed NGSim data ([29]). The calibration minimized the root mean square error (RMSE) between the simulation and measurements. The vehicle (ID 1013) was recorded for only 55 seconds, but even in that time many hysteresis loops were recorded in the speed-headway plane. For the NGSim data, which was recorded on the freeway in congested conditions, trajectories that qualitatively resemble this one are the rule rather than the exception. The different models all give a different caricature of the oscillatory behavior observed. In particular, the Newell (TT) model (with a simulation time step of .1 s) appears to better capture the seemingly chaotic nature of the oscillations, whereas the more complicated (and computationally expensive to calibrate) IDM or OVM models result in a lower RMSE and better capture the width of the hysteresis (the differences in speeds at a given headway). The point of this figure is not to make a case for one model over the other, but rather to point out the interesting behavior shown even in a single short trajectory. It also raises an important question which existing research seems unprepared to answer.

With so many different models, each of varying complexities, which ones best describe traffic flow? Our work on calibration is motivated by this question.

2.2 Delay-Differential Equations

Compared to ODEs, DDEs are more complicated to analyze and solve. Car-following models typically only consider the simplest case, where time delay (reaction time) is a constant. It seems the distinction between ODEs and DDEs has been given limited importance in the calibration literature, even though the addition of time delay as a parameter causes extra complications. As for the literature as a whole, the effect of nonzero time delay on car-following models is still not fully understood.

When studying the behavior of car-following models, one common approach is to consider a line of cars with identical model parameters on a closed loop. In this special setting, one can find an “equilibrium” solution where all vehicles have the same speed and same headway, with zero acceleration. The equilibrium solution is contrasted with an oscillatory/periodic solution, where the speed of a vehicle can oscillate in a manner reminiscent of stop-and-go or congested traffic. Since most analytic (i.e. not numerical) results come from studying the behavior around the equilibrium solution, the stability of the equilibrium solution is particularly important. Adjusting parameters in a model can be said to have a stabilizing (destabilizing) effect if the equilibrium solution will persist (degenerate) under a wider variety of conditions. Intuitively, one expects that drivers having slower reaction times will have a destabilizing effect, and this seems to be generally correct.

It has been shown for both the optimal velocity model (OVM) and intelligent driver model (IDM) that the inclusion of time delays has a destabilizing effect (see [30], [31], [32]). But time delay is not the only model parameter that can have such effects. For example, [33] showed that increasing (decreasing) the OVM’s “sensitivity parameter” has stabilizing (destabilizing) effects (the sensitivity parameter is proportional to how strongly the driver accelerates). Time delay can also have more subtle effects: [31] found the inclusion of time delays for OVM caused bistability between equilibrium and oscillatory solutions.

In [34] the OVM is considered with different reaction times $\leq 0.3$ seconds. In this range, no qualitative changes to the model are observed. The authors state that changing the sensitivity parameter can account for the effects of reaction time. [23] and [30] have done experiments on OVM with larger reaction times ($\approx 1$ seconds) and found that while it is true that reaction time is not especially important at low values, once increased past some threshold it will start to cause unstable behavior leading to vehicle collisions. A similar behavior was found in [35] for the Gipps car-following model. These findings show the inclusion of reaction time will have a destabilizing effect, but large reaction times seem to have especially unstable effects, perhaps owing to some sort of phase transition.

In practice, one can observe that calibrated values for reaction time greatly vary across the literature, even for the same model. Consider that [36],[15], and [16] all calibrate the Gipps model, but all find very different means and variances for drivers’ reaction times. One found the mean to be 0.57 s with a variance of 0.02 s. Another found the mean to be 1.73 s with variance 1.35 s.

Besides reaction time, another interesting extension to a model is the inclusion of “anticipation” effects. Examples of anticipation are a driver estimating the future speed/position of a leading vehicle, or a following vehicle reacting to multiple vehicles instead of just its leader. [32] showed that anticipation has a stabilizing effect for IDM, and [37] showed the same for OVM. Of course, due to nonlinearities, stabilizing and destabilizing effects will not simply cancel out.

2.3 Stochastic and Data-Driven Models

In this paper we consider the calibration of car following models formulated as ODEs or DDEs. This section discusses some recent research directions in which car following models are formulated as stochastic differential equations (SDEs), or derived from data using nonparametric models.

Papers using SDEs to model traffic flow typically argue that the stochastic models are better able to describe the formation and growth of congestion. Besides being able to achieve a closer fit to data, authors have concluded that deterministic models cannot reproduce some qualitative features of congestion. In particular, some recent field experiments involving 25 and 51 human driven vehicles on a straight road ([38]) showed that the standard deviation of traffic speeds initially grow in a concave way, whereas deterministic models seem to only be able to show convex growth.

SDEs in the literature are usually based on modifying an existing deterministic ODE model rather than created as an entirely novel model. Two common stochastic mechanisms are acceleration noise, and action points. For acceleration noise, white noise is added to the output of the deterministic model. The white noise may be sampled from either a uniform or normal distribution. Closely related to acceleration noise, is the case where the noise is added to a model parameter instead of directly to the output of the model. In [39] and [40] gaussian acceleration noise was added to the Newell model and IDM respectively. In 2D-IDM from [41] uniform noise is added to a model parameter (in that case, the desired time headway parameter).

Action points are implemented by creating a model regime where the driver is insensitive to changes in the headway and/or lead vehicle velocity. Thus when inside the regime, a driver will maintain their current behavior until a sufficiently large difference in the state is achieved, at which point the driver will react to the now significantly different driving situation. The reader is referred to [42] or [43] for specific details on implementing action points.

The motivation for including white noise is to model random human behavior and reflect the fact that humans cannot perfectly control their vehicles, whereas the motivation for action points is that drivers are not able to perceive minute changes to the traffic situation and thus a sufficiently large difference is required before they can start to react. [42] shows that action points and white noise have similar effects on traffic simulations, and suggests that both mechanisms should be included.

Data-driven car following models are based on machine learning models and algorithms, and as a result their formulation is significantly different than the ODE/DDE/SDE discussed so far. However, despite their different formulation, the inputs (current headway, follower and leader velocities) and output (current follower acceleration) are typically very similar. Papers in the literature have reported that data-driven models can achieve better fits to empirical data compared to parametric models such as OVM or IDM. The discussion of these machine learning techniques is outside the scope of this paper but we mention the different approaches which have been successfully applied. In [44] $k$-nearest neighbor (knn) is used, and [45] uses locally weighted regression (loess), which can be thought of as a generalization of knn. [46] discuss using feedforward neural networks with a single hidden layer and ReLU nonlinearity. That paper also uses recurrent neural networks (RNNs) to overcome some difficulties of the feedforward neural networks. [47] uses long short-term memory neural networks, a special kind of RNN. [48] use deep reinforcement learning (deep RL) and policy gradients, and do a comparison between deep RL, IDM, RNN, loess, and feedforward neural networks.

3 Calibration

It is useful to rewrite the second order differential equation Eq. $\eqref{1}$ as a system of first order equations:

[TABLE]

where $x_{n}$ is now a vector consisting of the position and velocity of vehicle $n$ at time $t$. $x_{n}^{*}$ and $v_{n}$ refer to the position and speed for vehicle $n$. $p_{n}$ denotes the parameters for vehicle $n$. $L(n,t)$ refers to the leader of vehicle $n$ at time $t$, and we denote the lead trajectory as $x_{L(n)}$ (time dependence suppressed) to reflect the fact that the index corresponding to the lead vehicle may change due to lane changing. Eq. (3) is what we refer to when we say “car-following model” and its solution $x_{i}(t)$ is a “simulated trajectory”.

The calibration is done for some arbitrarily sized platoon of $n$ vehicles, indexed as $1,\ldots,n$. Each vehicle has some individual parameters $p_{i}$ and a corresponding car following model $h_{i}$, so all vehicles have their own set of parameters, and $p_{i}$ only describes vehicle $i$. Let $p$ be a concatenated vector of all the parameters, so $p=[p_{1},p_{2},\ldots,p_{n}]$. Calibration finds the parameters $p$ which minimize the total loss between the simulated trajectories $x_{i}(t)$ and the corresponding measured vehicle trajectories, which are denoted $\hat{x}_{i}(t)$. The loss function is denoted as $f(x_{i},\hat{x}_{i},t)$ and represents the goodness of fit of the simulated trajectory $x_{i}$ at time $t$. The measurements $\hat{x}_{i}(t)$ are known in the interval $[t_{i},T_{i}]$ for each $i$.

In order to compute the simulated trajectories for each vehicle $i$, we require a lead vehicle trajectory $x_{L(i)}(t)$ as well as an initial condition $x_{i}(t_{i})$. The time interval that the lead vehicle trajectory is known is denoted $[t_{i},T_{i-1}]$. $x_{i}(t_{i})$ is the initial condition, which is specified as $x_{i}(t_{i})=\hat{x}_{i}(t_{i})$ for an ODE model.

For an ODE, the optimization problem is

[TABLE]

The objective function $F$ is defined as the integral of the loss function $f$ over time, summed over each simulated vehicle. The constraints are that the simulated trajectories $x_{i}$ obey the car-following model $h_{i}$ from time $t_{i}$ to $T_{i-1}$, with the initial conditions $x_{i}(t_{i})=\hat{x}_{i}(t_{i})$. There may also be some box constraints $b_{\rm low},b_{\rm high}$ which prevent parameters from reaching unrealistic values.

The DDE case is similarly defined:

[TABLE]

where $T_{i-1}^{*}=\min\{T_{i-1}+\tau_{i},T_{i}\}$, and the reaction time for vehicle $i$, $\tau_{i}$, is assumed to be a model parameter. $x_{i}(t)-\hat{x}_{i}(t)=0$ defines the history function for $x_{i}$ (DDEs have history functions as their initial conditions).

3.1 Downstream Boundary Conditions for Car Following Models

In Eqs. (5) and (6) we have specified the initial conditions, which give the times vehicles enter the simulation and their position/speed upon entering. These conditions might also be referred to as upstream boundary conditions since they state how vehicles enter the simulation. For video data such as NGSim, which is recorded on a fixed section of road, downstream boundary conditions also need to be specified. Video data has staggered observation times ($T_{i-1}<T_{i}$) meaning that a vehicle’s leader leaves the study area before the vehicle itself. Figure 3 gives a visual depiction of this.

Car following models require a lead vehicle trajectory as input, therefore once the lead vehicle leaves the simulation, the car following model can no longer be used. Downstream boundary conditions then give a rule for the dynamics a vehicle follows once its lead vehicle leaves the simulation (i.e. they specify the dynamics for $x_{i}(t),\ \ t\in[T_{i-1},T_{i}]$). The the boundary conditions we implement in this paper are

[TABLE]

meaning that when a vehicle’s leader is no longer available, the vehicle will move with whatever speed was recorded in the data. Under these boundary conditions, the velocity $v_{i}(t)$ and acceleration $\dot{v}_{i}(t)$ will be piecewise continuous (they will have a jump discontinuity at $T_{i-1}$). The objective $F$ and gradient $dF/dp$ are still continuous under this boundary condition (see sections 3.3 and 4.1.2 for the relevant discussion on discontinuities in car following models and how they impact the calibration problem). Note that for a DDE, the time $T_{i-1}^{*}$ is used instead of $T_{i-1}$.

It is necessary to specify the downstream boundary conditions for two reasons. First, it allows vehicles to exit the simulation area. Second, it allows congestion to propogate into the simulation area from downstream. This happens when the boundary conditions specify a reduction in speed.

3.2 Reaction Time as a Parameter

Simulated trajectories are computed numerically with a specified time discretization. Since we want to compare the simulated trajectories with the actual measurements, it makes sense for the simulation to use the same time discretization as the data. For example, if the data is recorded every 0.1 seconds, the simulated trajectories should use a time step of 0.1 seconds. If the simulation used a different time step (e.g. 0.09), then the loss can only be calculated at a few points (e.g. 0.9, 1.8) unless we define some interpolation procedure.

For a DDE, the time discretization and reaction time need to be multiples of each other. If the simulation uses a time step of 0.1 seconds, then the reaction time will have to be an integer multiple of 0.1 (e.g. 0.1, 0.2). We can either:

Enforce that every vehicle’s reaction time is a integer multiple of 0.1 2. 2.

Allow reaction time to be a continuous variable and define the loss function through interpolation.

If we choose 1, it leads to a mixed integer nonlinear programming problem (MINLP). This is a much harder problem than the original nonlinear program. One approach used in [16] would be to split up all vehicles into small platoons (either pairs or triplets) and use a brute force search for reaction time. However, even when the platoons are small, using brute force vastly increases computation because the calibration has to be repeated for every possible combination of reaction times.

The other option is to use an approach similar to [23] where reaction time is a continuous variable and interpolation is used on the simulated trajectories. In general, any sort of interpolation could be used, but linear interpolation is the most common. In this approach, interpolation needs to be used on the lead trajectory and initial history function in order to define the model. Assuming the numerical time step is not equal to the delay, interpolation will need to be used on the simulated trajectory as well.

One issue with interpolation is how to compute the loss function. Specifically, assuming the simulation is defined on a different set of times than the measurements, should the simulation be interpolated onto the measurement’s times or vice versa? We suggest doing both. This avoids the possibility of error cancelation occurring during the interpolation and causing the loss function to be erroneously small. In this scenario it would be preferable to reweight the loss function as well since the loss will be computed over more points.

3.3 Multiple Regimes, Lane Changing, and Discontinuities in Car Following Models

We say a model has multiple regimes when the model is expressed either as a piecewise function, or by using $\min$ and $\max$ functions. These types of modifications are frequently applied in practice. Perhaps the most common are examples like Eq. (2.1), where the maximum speed is explicitly constrained by some constant $v_{f}$. For second order models, often the acceleration will be constrained between minimum and maximum accelerations $a_{\rm min},a_{\rm max}$. Besides these simpler regimes which implement bounds on the acceleration and/or velocity, there are also models which use different regimes to describe different car following behavior. The Gipps model ([49]) is one such example, which has one regime to describe normal car following behavior, and a second regime which implements an emergency braking behavior. Other examples of multi-regime models are the constant acceleration heuristic proposed in [50], the follower stopper controller proposed in [51], or any model implementing action points. We consider

[TABLE]

as a general form of a two regime car following model, with the two regimes $h_{i,1}$ and $h_{i,2}$. The condition for the associated regime to become activated is $g_{i,2}$ (which is possibly a vector valued function). For example with Eq. (2.1): $h_{i,1}=(x_{L(n)}-x_{n}-s_{\rm jam})/\alpha$, $h_{i,2}=v_{f}$ and $g_{i,2}=h_{i,1}-h_{i,2}$. There must always be $m-1$ such conditions for a model with $m$ regimes, so a model with more regimes will simply have $h_{i,2},\ldots,h_{i,m}$ and $g_{i,2},\ldots,g_{i,m}$.

In the case of a first order model’s velocity being bounded, or a second order model’s acceleration being bounded, the resulting model $h_{i}$ is everywhere continuous (because both regimes will be equal when the switching condition $g_{i,2}$ is met). In general a multi-regime model need not be everywhere continuous. Gipps and follower stopper are both examples of a multi-regime model being everywhere continuous; the constant acceleration heuristic or a model implementing action points are examples not guarenteed to be continuous when switching between regimes. Of course, even for a multi-regime model which is everywhere continuous, its derivatives will have discontinuities when switching between regimes. Note that a multi-regime model written using $\min$ and $\max$ will always be continuous. The reason why we prefer using the form Eq. (8) is because 1) it is more general, 2) we can explicitly refer to the switching condition $g_{i,2}$ and 3) the regime of the model is explicitly stated, which avoids the problem of the regime being ambiguous when the quantities inside a $\min/\max$ are equal.

Lane changing is another issue which causes discontinuities in car following models. Specifically, whenever a lead vehicle changes, i.e. $L(i,t)$ changes, then the lead vehicle trajectory $x_{L(i)}$ will have a jump discontinuity, and therefore $h_{i}$ will have a jump discontinuity as well.

Given that lane changing always causes discontinuities and multi-regime models may cause discontinuities, one may worry whether or not the objective $F$ will have discontinuities, and whether or not the gradient $dF/dp$ is even guarenteed to exist. We give sufficient conditions for the continuity of $F$ and $dF/dp$ in section 4.1.2.

3.4 Algorithms for Solving the Calibration Problem

Eqs. $\eqref{5}$ and $\eqref{6}$ are nonlinear and non-convex optimization problems that must be solved numerically. Current approaches usually rely on gradient-free optimization, with genetic algorithm (GA) and Nelder-Mead simplex algorithm (NM) being the two most commonly applied algorithms ([12, 52, 13, 14, 53, 23, 15, 16]). These methods have the disadvantage of being slow to converge, and require many objective function evaluations. The objective $F$ is relatively expensive to compute in this context because we first must solve a system of differential equations (i.e. compute the simulated trajectories). As an alternative to directly solving the optimization problems $\eqref{5}$ and $\eqref{6}$, parameters can be estimated directly from data using statistics techniques ([54, 55, 56, 57]) but this approach is not guaranteed to find local minima. Also, models may have parameters which do not have clear physical interpretations, and therefore cannot be directly estimated.

Gradient-based optimization algorithms can potentially solve the calibration problem faster than gradient-free approaches. However, even though most optimization algorithms use gradient information, not many gradient-based algorithms have been examined in the calibration literature. To the authors’ knowledge, the only gradient-based algorithms that have been considered are simultaneous perturbation stochastic approximation (SPSA) (see [58]), Lindo’s Multistart (a commercial solver’s global optimization routine, see [59]), and sequential quadratic programming using filtering (SQP-filter, see [60]). In [61] SQP-filter was found to give better solutions than NM or GA, as well as solve the problem up to 10 times faster; even better results were found when it was additionally combined with gradient-free algorithm DIRECT ([62]). [35] did not consider computational speed, but found the Multistart algorithm had a much higher chance of finding the true global minimum compared to GA or NM.

Even if one only knows how to explicitly evaluate the objective function $F$, it is not necessary to use gradient-free optimization because the gradient (the vector $dF/dp$) can be calculated numerically. To do this, one varies the parameters $p$ one at a time, and recomputes the objective each time a parameter is varied. Then the gradient can be calculated using finite differences. We refer to this as the finite difference approach for computing the gradient, and it requires us to recompute the simulated trajectories of every single vehicle $m$ additional times, where $m$ is the number of parameters. When using a forward euler scheme with stepsize $\Delta t$ for an ODE model, the computation cost for computing the simulated trajectories once is $\mathcal{O}(T(n))$, where $T(n)=\sum_{i=1}^{n}(T_{i-1}-t_{i})/\Delta t$. So computing the gradient with finite differences has a complexity of $\mathcal{O}(mT(n))$. The adjoint method has complexity $\mathcal{O}(T(n))$, and so it is the prefered method for computing the gradient when $m$ is large ([63, 64]).

4 Adjoint Method and Gradient-Based Optimization

In this paper we are primarily concerned with applying the adjoint method to Eq. $\eqref{5}$ in order to efficiently evaluate the gradient $dF/dp$ for an ODE model. The gradient can then be used with any gradient-based optimization routine. In this section we derive the gradient and compare it to using finite differences; we consider different optimization algorithms in the next section.

The adjoint method can also be applied to the calibration problem $\eqref{6}$ where the model includes time delay. It can also be used to calculate the hessian. Both of these extensions follow the same techniques elaborated in this section, and are included in the appendix.

We first consider the calibration of a single vehicle ($n=1$) before dealing with the general case of an arbitrary number of vehicles.

4.1 Adjoint Method for ODEs

In our notation, $x_{i}$, $h_{i}$, and $p$ are vectors, and $f$ and $F$ are scalars. A derivative of a scalar with respect to a vector is a vector, so for example $\partial f/\partial x_{i}=[\partial f/\partial x^{*}_{i},\partial f/\partial v_{i}]$. A derivative of a vector with respect to a vector is a matrix, where the rows correspond to the vector being differentiated, and the columns correspond to the vector being differentiated with respect to. For example,

[TABLE]

where $p^{j}$ represents the $j{\rm th}$ parameter. Additionally, the superscript T indicates transpose.

Define the augmented objective function, denoted $L$

[TABLE]

in this equation, $\lambda(t)$ is referred to as the adjoint variable. $\lambda(t)$ has the same dimension as the state variable $x_{n}(t)$, so it is a length two vector for a second order model. Note that since $\dot{x}_{n}-h_{n}(\cdot)=0$, we have that $L=F$. These equalities are always true regardless of the choice of $\lambda(t)$.

Following the observation that $L=F$, it follows that $dL/dp=dF/dp$, so differentiating Eq. (4.1) with respect to the parameters $p$ gives

[TABLE]

where quantities are evaluated at time $t$ if not explicitly stated. Note in this case there are no terms depending on the lead trajectory $x_{L(n)}$ since the platoon consists only of $x_{n}$ and the lead trajectory is regarded as fixed.

Eq. $\eqref{11}$ cannot be computed in the current form because we have no way to calculate ${\partial x_{n}}/{\partial p}$. The adjoint method works by choosing the the adjoint variable $\lambda(t)$ in a way which allows us to avoid the computation of the unknown ${\partial x_{n}}/{\partial p}$. When doing this, one will find that the correct choice of $\lambda(t)$ results in solving a similar system to the original constraint, but “backwards” (e.g. when the model is an ODE, $\lambda(t)$ will end up being defined as an ODE which is solved backwards in time). The equations which define the adjoint variable are known as the adjoint system.

To derive the adjoint system, first apply integration by parts:

[TABLE]

The quantity ${\partial x_{n}(t_{n})}/{\partial p}$ is known because it simply depends on the initial conditions; based on (5) it is simply zero. We will choose $\lambda(T_{n-1})=0$ to eliminate the ${\partial x_{n}(T_{n-1})}/{\partial p}$ term. This choice becomes the initial conditions for the adjoint system; note that $L=F$ regardless of how the initial conditions for $\lambda$ are chosen. Now combining Eqs. $\eqref{11}$ and $\eqref{12}$:

[TABLE]

Then to avoid calculating $\partial x_{n}/\partial p$ we will enforce:

[TABLE]

Eq. $\eqref{13}$ defines a new differential equation in $\lambda(t)$. Our choice $\lambda(T_{n-1})=0$ is the corresponding initial condition. Then after solving for $\lambda(t)$ we can calculate the desired gradient:

[TABLE]

4.1.1 N Car Platoon

Now we will apply the adjoint method to a platoon of arbitrary size, and sketch a general algorithm for computing the gradient. Define the augmented objective function and proceed as before. Again, we suppress $t$ dependence for clarity.

[TABLE]

[TABLE]

Here, $\lambda_{n}$ corresponds to the adjoint variables for vehicle $n$. Similar to how $L(n,t)$ was defined, $G(n,t)$ is defined as the index of vehicle $n$’s follower at time $t$, so that $x_{G(n)}$ is the following trajectory of vehicle $n$. If vehicle $n$ has no follower at time $t$, or if the follower is not part of the platoon being calibrated, let $G(n,t)$ return zero. $\mathbbm{1}$ is the indicator function. When $n>1$, the downstream boundary conditions must be put into the augmented objective function, because a vehicle $i$ may affect its follower in the time $[T_{i-1},T_{i}]$, even though it only contributes to the loss function during $[t_{i},T_{i-1}]$.

Apply integration by parts, take the initial condition as $\lambda_{i}(T_{i})=0$, and gather terms with $\partial x_{i}/\partial p$ to form the adjoint system:

[TABLE]

Comparing this adjoint system to the one for a single vehicle (4.1), the difference here is that there is an extra contribution which occurs from the coupling between vehicles. This coupling term occurs only when a vehicle acts as a leader for another vehicle in the same platoon. After solving the adjoint system,

[TABLE]

All the partial derivatives in Eqs. $\eqref{16},\eqref{18}$ are computed exactly for a given loss function and car-following model. Solving for $x_{i}(t),\lambda_{i}(t),\text{ and }dF/dp$ is done numerically with an appropriate time discretization. The same discretization should be used for all quantities, so if one uses a forward euler scheme for $x_{i}$, the same forward euler scheme should be used for $\lambda_{i}$ and the integral in $F$ should use a left Riemann sum. Note as well that for a multi-regime model, one should keep the regime for each timestep in memory, because otherwise the switching condition $g_{i,n}$ will have to be evaluated at each timestep when computing the adjoint variables.

We define a general algorithm for computing $F$, $\dfrac{d}{dp}F$:

**Inputs: **measured trajectory data $\hat{x}_{1}(t),\ldots,\hat{x}_{n}(t)$, any necessary lead vehicle trajectories $x_{L(i)}(t)$ for $i\in[1,\ldots,n],\ L(i,t)\notin[1,\ldots,n]$. Models $h_{1},\ldots,h_{n}$ with parameters $p_{1},\ldots,p_{n}$.

For each vehicle $i\in[1,\ldots,n]$:

1. (a)

Compute simulated trajectory $x_{i}(t)$, $t\in[t_{i},T_{i-1}]$ by solving $\dot{x}_{i}(t)=h_{i}(x_{i},x_{L(i)},p_{i})$ with intial condition $x_{i}(t_{i})=\hat{x}_{i}(t_{i})$. 2. (b)

Compute $x_{i}(t)$, $t\in[T_{i-1},T_{i}]$ using downstream boundary conditions Eq. (7) 3. (c)

Keep in memory: $x_{i}(t)$, $G(i,t)$, and the regime of $h_{i}$ at each time 2. 2.

Compute objective $F=\sum_{i}^{n}\int_{t_{i}}^{T_{i-1}}f(x_{i},\hat{x}_{i},t)dt$ 3. 3.

For each adjoint variable $i\in[n,\ldots,1]:$ compute $\lambda_{i}(t)$, $t\in[T_{i},t_{i}]$ by solving adjoint system Eq. $\eqref{16}$ with initial condition $\lambda_{i}(T_{i})=0$. 4. 4.

Compute gradient $dF/dp=-\sum_{i}^{n}\int_{t_{i}}^{T_{i-1}}\lambda_{i}^{T}\partial h_{i}/\partial p\ dt$

**Outputs: ** $F$, $dF/dp$

4.1.2 Sufficient Conditions for a Continuous Objective and Gradient

Theorem 1**.**

The following are sufficient conditions for $F$ and $dF/dp$ to be continuous:

the loss function $f(x_{i},\hat{x}_{i})$ and its partial derivative $\partial f/\partial x_{i}$ are piecewise continuous on $[t_{i},T_{i-1}]$ for all $i$ 2. 2.

the car following model $h_{i}(x_{i},x_{L(i)},p_{i})$ and its partial derivatives $\partial h_{i}/\partial x_{i}$, $\ \partial h_{i}/\partial p$, and $\partial h_{i}/\partial x_{L(i)}$ are piecewise continuous on $[t_{i},T_{i-1}]$ for all $i$ 3. 3.

$\partial h_{i}/\partial x_{i}$* must be able to be bounded by a constant* 4. 4.

for a multi-regime model, the Eq. 17 must be satisfied

where it is additionally assumed that $x_{L(i)}(t),L(i,t)\notin[1,\ldots,n]$ (any lead vehicle trajectories which are required but not simulated) are piecewise continuous and $\hat{x}_{i}(t),i\in[1,\ldots,n]$ are continuous.

Corollary 1**.**

If in addition to $F$ and $dF/dp$ being continuous, $\partial f/\partial x_{i},\ \partial h_{i}/\partial x_{L(i)},\text{and}\ \partial h_{i}/\partial p$ are bounded, then $F$ is Lipschitz continuous.

We conclude that a multi-regime model need not be continuous when switching between regimes, and lane changing, which will always cause jump discontinuities in $h_{i}(x_{i},x_{L(i)},p_{i})$, is allowed. Note that by piecewise continuity we mean a model (or loss function) can switch regimes only a finite number of times in a finite time interval (and each regime is continuous); additionally, the time interval of each regime must be continuous with respect to the model parameters. Refer to the the full proof in the appendix for further details. We also show in the appendix that the objective is Lipschitz continuous for the OVM, which was the model used for the numerical experiments in this paper.

4.2 Implementation of Adjoint Method and Calibration Problem

To test the adjoint method and optimization algorithms, the calibration problem Eq. (5) was implemented in python for the optimal velocity model. The equations below show the functional form of the model as originally formulated in [57] and [26].

[TABLE]

$c_{1}$ through $c_{5}$ are the 5 parameters of the model. $s$ is the space headway, defined as $s=x_{L(i)}-x_{i}-l_{L(i)}$, where $l_{i}$ is the length of vehicle $i$. All vehicles have their own unique parameters, so the total number of parameters in the optimization problem is 5 times the number of vehicles.

The optimal velocity model and its adjoint system were both discretized with a forward euler scheme. The data used was the reconstructed NGSim data ([29]), so a time step of .1 seconds was used to be consistent with this data. The loss function $f$ used was square error, which was discretized with a left Riemann sum.

[TABLE]

which is equivalent to minimizing the root mean square error (RMSE)

[TABLE]

Since only a car following model is being calibrated, any lane changes in the measurements will be the same lane changes in the simulation. This also means that the leader/follower relationships in the simulation are the same as they are in the data. This allows trajectories to be calibrated at the level of individual vehicles, since it does not make sense to directly compare the simulated and measured trajectories of a vehicle if the simulation and measurements have different lead vehicles.

4.3 Adjoint method compared to other methods for calculating the gradient

Many optimization algorithms will automatically use finite differences to compute the gradient if an explicit gradient function is not given. Simultaneous perturbation is another option when using an optimization algorithm like the one outlined in [58]. Then, assuming that one has chosen to use a gradient-based algorithm for calibration, it is of practical interest to investigate how the adjoint method compares to both of these methods. Specifically, how does the speed and accuracy of these three methods compare in practice?

4.3.1 Comparison of Speed

The computational cost of an objective function evaluation is $\mathcal{O}(T(n))$, where $T(n)$ is the number of simulated timesteps, added over all $n$ vehicles. The theoretical computational complexities for the adjoint method is $\mathcal{O}(T(n))$, while forward differences has $\mathcal{O}(mT(n))$, where $m$ is the number of parameters. The simultaneous perturbation method has $\mathcal{O}(kT(n))$ complexity, where the final gradient is the average of $k$ estimates of the gradient. Note that these computational complexities are only for the gradient calculation and do not consider the complexity of whatever optimization algorithm the gradient is used with.

To test these theoretical results in practice, the gradient of the calibration problem was calculated and the platoon size was varied. Since the OVM was used, the total number of parameters is 5 times the number of vehicles. We considered a platoon having between 1 and 15 vehicles. The results are shown in the figure 4.

Figure 4’s left panel shows the total CPU time, measured in seconds, needed to calculate the gradient of the calibration problem using the three different methods. The simultaneous perturbation was done with both 1 and 3 trials, denoted SPSA and SPSA-3. The right panel shows the time needed to calculate the gradient divided by the time needed to calculate the objective. To minimize the randomness in the CPU times, each data point was repeated 10 times and the results were averaged. Some variation is still present because some vehicles are observed for shorter times. The same 15 vehicles were used in the experiment and always added in the same order. It should also be stated that for either a finite differencing scheme or for the adjoint method, one needs to first compute the objective before computing the gradient. Therefore, it should be implicitly understood that “time to calculate the gradient”, really refers to “time to calculate both the objective and gradient” for those two methods. This shows that using a finite difference method quickly becomes intractable when the number of parameters increases. The time increases quadratically in the left panel because both $m$ and $T(n)$ are increasing. In the right panel, when the time needed to calculate the gradient is divided by the time needed to calculate the objective, one sees the complexity growing as a linear function of $m$. In that figure, we also see that the complexities for the adjoint method and simultaneous perturbation do not depend on $m$.

The table 1 shows the cpu times for 5-15 parameters. One can observe that this implementation of the adjoint method for OVM requires the equivalent of 4 objective function solves, whereas forward differences requires $m+1$ (where there are $m$ parameters).

4.3.2 Comparison of Accuracy

For the calibration problem, there is no way to calculate the gradient exactly because there is no closed form solution for the car following model. To evaluate the accuracy of the adjoint method and simultaneous perturbation, we will therefore treat the result from forward differences as the true solution. In the experiment, points were evenly sampled between an initial starting guess and the global minimum of the calibration problem for a single vehicle. At each of these points in the parameter space, the relative error was calculated

[TABLE]

where $\nabla_{fd}F$ is the gradient calculated with forward differences, $\nabla F$ is the gradient calculated with the adjoint method or simultaneous perturbation, and $||\cdot||_{2}$ denotes the 2 norm. The results of the experiment are shown below.

For the adjoint method, the relative error is usually small, on the order of $10^{-6}$. When the parameters start to become close to a local minimum, the relative error starts to rapidly increase to the order of $10^{-4}-10^{-3}$. When at the optimum, the relative error becomes as high as .3. The reason why the relative error increases is because as the parameters approach the minimum, the norm of the gradient decreases. The right panel shows the relationship between the log of the norm of the gradient calculated with finite differences, and the log of the relative error. One can see that they are inversely proportional to each other. This makes theoretical sense: we expect the residual $\nabla_{fd}F-\nabla F$ to always be on approximately the same order of magnitude since it depends mainly on truncation error which does not change with the parameters. Then, it follows that the relative error is dominated by the $||\nabla_{fd}F||$ term in the denominator, which is what we observe.

The relative error for simultaneous perturbation follows the same shape as the adjoint method, but is always about 6 orders of magnitude larger. Of course, you would not use simultaneous perturbation as a substitute for the finite difference gradient, since it is used with its own optimization algorithm. Nonetheless, this shows that although the two methods (adjoint and simultaneous perturbation) both have a flat computational cost with the number of parameters, they have different uses. The reason why the the error for simultaneous perturbation is quite high in this context is because the gradient is scaled very poorly for the OVM model, meaning the parameters have very different sensitivities. This poor scaling is a common feature for car following models.

5 Benchmarking different Algorithms and the Adjoint method

The ideal optimization algorithm would be one that can not only solve the problem quickly, but also one that is able to avoid bad solutions corresponding to local extrema. Because existing literature has focused on gradient-free optimization algorithms, not much is known about using gradient/Hessian-based optimization to solve the calibration problem. The only works the authors know of in the car following calibration literature that consider gradient based optimization are [61] and [35]. Both those works found the gradient based algorithms to work best overall. This result is consistent with other engineering applications; for example [65] and [66] deal with aerodynamic shape optimization and material science respectively, and report that gradient-based optimization can offer significant speed increases.

To establish the benefits of gradient based optimization and the adjoint method, we consider 2 gradient-free and 5 gradient-based optimization algorithms. The different algorithms, their specific implementations, and abbreviations used throughout the rest of the paper, are detailed below.

• •

Genetic Algorithm (GA) - an initial population of trial solutions is randomly selected throughout the parameter space. In each iteration, the population is altered through stochastic processes referred to as “mutation” and “crossover”, and then updated in a selection step. An implementation following [67] is part of the scipy package in python.

• •

Nelder-Mead Simplex (NM) - a given initial guess is used to define a collection of points (a simplex). In each iteration, the corners of the simplex are updated according to deterministic rules. Implemented in scipy following [68]. NM works on unconstrained problems, so a large penalty term was added to the objective if the parameters went outside the bounds.

• •

Limited memory BFGS for bound constraints (BFGS) - at every iteration, gradient information from the previous iterates is used in the BFGS formula to form an approximation to the hessian. A quadratic approximation to the objective is formed, and the generalized cauchy point of this approximation is then found, which defines a set of active bound constraints. A minimization is then performed over the variables with inactive constraints, while the variables with active constraints are held fixed. The solution to this last minimization finally defines the search direction for the algorithm, and a line search is performed to determine the step length. The algorithm, referred to as l-bfgs-b, is due to [69] and [70] and wrapped as part of the scipy package.

• •

Truncated Newton Conjugate (TNC) - at every iteration, instead of explicitly computing the Newton direction $-\nabla^{2}F^{-1}\nabla F$, the conjugate gradient method is applied to the newton system $\nabla^{2}Fd=-\nabla F$ (where $d$ is the search direction to be determined). Conjugate gradient is itself an iterative algorithm, and it is terminated early (truncated) to give an approximate solution to the newton system which is also guaranteed to define a direction of descent. A line search is then performed to define the step length. The conjugate gradient algorithm only needs hessian-vector products which can be computed efficiently using finite differences, so the full hessian never needs to be calculated. An implementation due to [71] is wrapped as part of scipy. That implementation uses a conjugate gradient preconditioner based on the BFGS update, and the search direction is projected to enforce the bound constraints.

• •

Simultaneous Perturbation Stochastic Approximation (SPSA) - At each iteration, the gradient is approximated using simultaneous perturbation. Then the algorithm moves in the direction of the negative gradient with a fixed step length, so no line search is performed. The size of the step length is gradually reduced in subsequent iterations. Implemented in python following [58]; the step sizes and number of iterations were manually tuned by hand.

• •

Gradient Descent (GD) - At every iteration, the negative gradient is computed and scaled according to the spectral step length (also known as the Barzilei Borwein method). The scaled gradient is then projected onto the bound constraints to define the search direction, and a nonmonotone backtracking linesearch is used to define the step length. Implemented according to [72] in python. The algorithm was tested with different line searches (backtracking, weak/strong wolfe, watchdog, see [73]) and the nonmonotone backtracking gave the best performance.

• •

Sequential quadratic programming (SQP) - At every iteration the hessian is computed and the search direction is defined by the newton direction projected onto the bound constraints. This algorithm explicitly computes the hessian instead of approximating it. Because the hessian may fail to be positive definite, the search direction is safeguarded so that if the newton direction fails to define a direction of descent, the algorithm instead searches in the direction of the negative gradient. A small regularization is also added to the hessian so that it can always be inverted. The algorithm is a simple python implementation of sequential quadratic programming based on [73]. We again found the nonmonotone backtracking linesearch to be most effective.

In summary, out of all the algorithms tested: NM and GA are gradient-free; SPSA and GA are stochastic; BFGS and TNC are quasi-newton; GD and SPSA do not use hessian information; SQP explicitly computes the hessian; GA is the only global algorithm.

For the TNC, BFGS, GD, and SQP algorithms, we considered two separate modifications for each algorithm. First, either the finite difference method or the adjoint method can be used to compute the gradient, and we denote this by either “Fin” or ”Adj” respectively. Second we considered starting from multiple initial guesses. Multiple initial guesses are specified, and after each calibration, if the RMSE was above a specified threshold, the next initial guess would be used. The same 3 initial guesses were used in the same order for all algorithms requiring an initial guess. So with a threshold of 0, 3 guesses are always used, and with an arbitrarily large threshold, only a single guess will be used. The thresholds of 0, 7.5 and $\infty$ were used, and denoted in reference to the algorithm. The initial guesses used are based off of the calibrated parameters found from [57].

5.1 Evaluating Algorithm Performance

As detailed in section 4.2, we tested the calibration problem on the OVM car following model using squared error in distance as the loss function (equivalent to minimizing the RMSE). In this experiment, for each algorithm, every vehicle in the reconstructed NGSim dataset was calibrated using the true measurements of its lead trajectory. In total, there were 7 unique algorithms and 23 total algorithms tested when including the variants described above.

The calibration results are evaluated primarily based on three different measures: the average calibration time for a single vehicle, the % of the time the algorithm found the global minimum, and the average RMSE of a calibrated vehicle. Since there is no way to know the actual global minimum for the calibration problem, we regard the global minimum for a vehicle as being the best result out of all the algorithms, and we regard an algorithm as having found the global minimum if its RMSE is within 1/12 ft (1 inch) of this best result.

Based on these three metrics, all algorithms on the pareto front were identified and are shown in the table 2. In figure 6 every algorithm is plotted in the RMSE - time and % found global optimum - time plane. The algorithms on the pareto front are circled in black.

The only three algorithms on the pareto front were the GA, TNC, and BFGS. GA gave the best % of finding the global optimum. TNC with the adjoint method and three guesses gave the best overall RMSE, and was the fastest algorithm overall when used with the adjoint method and a single guess. BFGS was slightly slower than TNC and has a slightly higher RMSE, but it has a higher chance of finding the global optimum.

The conclusion of what algorithms are on the pareto front depends strongly on how “finding the global optimum” is defined. Here, we treat an algorithm as having found the global optimum if it gives a result within an inch (1/12) of the best RMSE. Let us refer to this value of 1/12 as the “tolerance” for finding the global optimum. If a smaller tolerance is used, then GA will no longer be on the pareto front. If a larger tolerance is used, GA and BFGS will no longer be on the pareto front. For example, with a tolerance of 0, the pareto front consists of only the adj TNC and adj BFGS variants, and adj BFGS-0 gives the highest chance of finding the global optimum. With a tolerance of 1/2, the adj TNC variants are the only algorithms on the pareto front. The effects of changing the tolerance are shown in the table 3. Adj TNC will always be on the pareto front regardless of the tolerance because it gives the best average RMSE for any given speed. No algorithms other than adj TNC, BFGS, and GA are ever on the pareto front.

The behavior shown in table 3 is explained by figure 7. That figure shows how the % found global optimum metric changes with tolerance. For example, a tolerance of 1/12 corresponds to -1.08 in a log scale, and the top right panel shows that at that point, the GA has the highest chance of finding the global optimum. The same panel shows that if the tolerance increases, adj TNC-0 will have the highest chance, and that adj BFGS-0 will have the highest chance if the tolerance decreases.

BFGS typically found the best solution out of all the algorithms tested, but its distribution of RMSE has a heavier tail in the sense that there were also some vehicles for which BFGS couldn’t find a good solution. TNC is most consistent because its distribution shows less of a heavy tail (it is the first to reach 100% found global optimum as the tolerance increases). The GA has a tail similar to BFGS and performs well when only a moderate tolerance is needed.

The appendix includes a plot that shows the distributions of RMSE (on a non-log scale), and the differences between the three algorithms are essentially indistinguishable at that resolution.

Overall, these differences between the three algorithms (GA, adj TNC-0, adj BFGS-0) are relatively minor, with all three giving essentially the same average RMSE (GA had a slightly worse RMSE than the other two). It’s not clear to the authors whether the subtle differences between them would cause any differences in practice. However, the BFGS and TNC algorithms both give a significant speed increase; the adjoint method with BFGS or TNC is about 5 times faster than the GA on average. Even when using finite differences instead of the adjoint method, BFGS is still about 3 times faster than the GA (finite differences with TNC did not work well). A much larger speed increase can be realized with only a small trade off in accuracy. Adj TNC-Inf had an average RMSE of 6.62 compared to 6.47 for the GA, but is 15 times faster. Adj TNC-7.5 had an average RMSE of 6.48 while being 9 times faster than the GA.

The comparison between the different strategies for initial guesses (0, 7.5, Inf) shows the trade off between speed and accuracy when starting a local search based algorithm multiple times. For car following models, it seems that even a suboptimal local minimum still gives a reasonable accuracy. The difference in accuracy made by starting a algorithm multiple times is smaller than the difference from using another algorithm, so regardless of the initial guess, the local search seems to converge to similar parameters. We also saw that a good local search algorithm (in this case, the quasi-newton algorithms) can give results equivalent to a global search algorithm (the GA) when a small number of initial guesses are used. Of course, each initial guess used greatly increases the time needed for the algorithm, so a hybrid strategy similar to the one proposed in [61] can be considered if a large number of initial guesses are needed to achieve good performance.

Comparing the results of using the adjoint method to finite differences, the differences between the two methods depends heavily on the algorithm used. For BFGS, using finite differences instead of the adjoint method gives essentially the same result, in the same number of objective/gradient evaluations, and the difference is speeds is consistent with the speed increase found in section 4. For TNC, the algorithm does not perform well when finite differences is used.

All the attention so far has been concentrated on the TNC, BFGS, and GA algorithms since those gave good results. None of the other algorithms performed well. NM is designed for unconstrained problems, so it is unsurprising that it did not perform as well as an algorithm designed specifically to deal with box constraints (recall that here the box bounds were enforced by adding a penalty term). GD suffered from the problem widely reported in the literature where it takes a very large number of iterations to fully converge. The majority of the time, the algorithm terminates due to exceeding the maximum number of evaluations allowed. For the SQP (recall in this case we explicitly computed the Hessian), we found that the newton direction often fails to define a direction of descent, and the algorithm simply searches in the direction of the gradient instead. Then, the algorithm spends a lot of computational time computing the Hessian but ends up not using it for anything, so explicitly computing the hessian seems to not work poorly for car following calibration; BFGS updating to estimate the Hessian, or a truncated newton method for estimating the newton direction should be used instead. SPSA performed worst out of all algorithms tested because of the vastly different sensitivities for car following models, as well as the rapidly changing magnitude of the gradient. A variant of SPSA such as c-SPSA ([74]) or W-SPSA ([75]) may have overcome these issues.

The analysis in this section has shown the value of gradient based algorithms and the adjoint method; compared to gradient free methods, they can give the same overall performance while offering a significant speed increase. We found that for the calibration of a single vehicle, the truncated newton (TNC) algorithm with the adjoint method gives a slightly more accurate calibration than a genetic algorithm, while also being 5 times faster. TNC is 15 times faster than the genetic algorithm if $\approx$ 2% decrease in accuracy is acceptable (when only a single initial guess is used). Actually, as we will see in the next subsection, when the calibration problem becomes larger, the benefit of gradient based algorithms and the adjoint method becomes even larger.

5.2 Calibration of larger platoons

Define platoon size as the number of vehicles calibrated at the same time ($n$ in Eq. (5)). The simplest case is when $n=1$: a single vehicle is calibrated at a time (this was the strategy used in the above section). Taking $n=1$ is how most of the literature treats car following calibration. Moreover, trajectories are typically calibrated to the measured, as opposed to simulated, lead trajectory. To the authors’ knowledge, [16] is the only paper which has taken a platoon size of larger than 1 ($n=2$ in that paper) and [46] is the only paper which considers the issue of predicting on the simulated vs. measured lead trajectory. In actual traffic simulations, all vehicles are simulated at the same time, and measured trajectories aren’t assumed to be available. In that case, errors from one vehicle can accumulate in the following vehicles and produce unrealistic effects. This motivates the question of calibrating platoons of vehicles, where each vehicle uses the simulated lead trajectory.

In this section, we are interested in two questions regarding the calibration for a platoon size of $n>1$. First, for the different methods considered so far (gradient free, gradient based with finite differences, adjoint method) how does the time needed to solve the calibration problem scale with the number of parameters. Second, what is the benefit (in terms of RMSE) from considering a larger platoon of vehicles.

A platoon of 100 vehicles was formed from the NGSim data. The platoon was formed by randomly selecting the first vehicle to be calibrated, and then repeatedly adding vehicles which have their leaders in the platoon. We considered calibrating these 100 vehicles with a platoon size $n=1,2,\ldots 10$. When the platoon size doesn’t divide 100, the last platoon used will be the remainder (e.g. for $n=3$, the last platoon is of size 1). Here the calibration is done sequentially, so each vehicle is calibrated to the simulated trajectory of its leader.

To compare how different algorithms scale, we performed calibration using the GA, Adj TNC, and Fin TNC algorithms, as shown in figure 8. For each platoon of a specific size, the number of equivalent objective evaluations is recorded. Gradient evaluations are converted into objective evaluations following the results of 4.3.1 (adjoint method can compute the gradient and objective in the equivalent of 4 objective evaluations, finite differences computes the gradient and objective in the equivalent of $m+1$ objective evaluations, where $m$ is the number of parameters). The results are shown in figure 8.

Assuming that a gradient based algorithm converges in approximately the same number of gradient/objective evaluations for either using finite differences or the adjoint method, the equivalent number of objective evaluations can be easily converted between the two methods. Namely, if using the adjoint method requires $g(m)$ total computational effort for computing the gradient, for some arbitrary function $g$, then finite differences will require $\frac{m+1}{4}g(m)$. So using finite differences will always scale an order of $m$ worse than the adjoint method. For the GA, it is clear that its cost is growing extremely fast with the number of parameters, and we only considered a platoon up to size 3 because of this reason. It was about 75 times as expensive as using TNC with the adjoint method for a problem size of only 15 parameters. We conclude that when the number of parameters becomes large, using the adjoint method with gradient based optimization can offer an arbitrarily large speed increase compared to either gradient free optimization or finite differences.

As for the effect of platoon size on the overall RMSE of all 100 vehicles, the results are shown in figure 9. Note overall RMSE is different than average RMSE because vehicles are weighted according to how many observations they have (see section 4.2). Also, the RMSE is expected to be higher here compared to the previous section since the calibration is being done sequentially using the leader’s simulated trajectory instead of the leader’s measured trajectory. The left panel shows the overall RMSE versus platoon size. The right panel shows the relative improvement in RMSE. To calculate the relative improvement, the percent change in RMSE compared to a platoon size of 1 is calculated for a specific platoon. Then the average and standard deviation of the percent change is calculated over all platoons of the specific size (e.g. 33 size 3 platoons or 10 size 10 platoons). The standard deviations are quite large: for any specific platoon, the improvement could be as large as 50%, or the RMSE could sometimes be significantly worse. For a platoon size $n=2$ or $3$, we saw a neglible overall improvement. In both those cases, we saw that the majority of platoons had an improvement of around $8\%-9\%$, but there were some platoons which achieved a significantly worse fit (20%-40%+ worse). These significantly worse fits suggest that in those cases, the algorithm was converging to bad local minima which apparently were avoided when calibrating vehicles one at a time. As the platoon size increased, the improvement became more consistent, and for platoon size $n=9$ and $10$ we saw an overall $13.8\%$ and $17.8\%$ improvement respectively. The authors stress that this improvement was achieved without changing the algorithm or the model at all. The only thing that was changed was the platoon size ($n$ in Eq. (5)).

Overall, these results suggest that considering a larger platoon size can improve calibration results, and further research on larger platoon sizes should be considered. In figure 9 we saw that the GA starts to perform better than TNC when the platoon size increases. We also saw that for larger platoon sizes, there were some platoons which ended up with a significantly worse fit compared to the case of $n=1$. This suggests that the local search algorithm (TNC) was having problems converging to bad local minima for $n>1$. We therefore conclude that while using multiple guesses was adequate for the calibration of a single vehicle, a more sophisticated strategy is needed when multiple vehicles are being calibrated.

6 Conclusion

We considered an optimization based formulation of the calibration problem for car following models being calibrated to trajectory data. It was shown how to apply the adjoint method to derive the gradient or hessian for an arbitrary car following model formulated as either an ordinary or delay differential equation. Several algorithms for solving the calibration problem were compared, and it was found that the best overall algorithms were the genetic algorithm, l-bfgs-b, and truncated newton conjugate.

For the calibration of a single vehicle at a time, it was found that using the adjoint method and a quasi-newton method gives slightly better performance than a genetic algorithm, and is 5 times faster. As the number of parameters increases, the speed increase offered by the adjoint method and gradient based optimization can become arbitrarily large (for 15 parameters, a truncated newton algorithm using the adjoint method was 75 times faster than a genetic algorithm). The adjoint method will always perform better than finite differences as the number of parameters becomes large, because computing the gradient with the adjoint method has a flat cost with respect to the number of parameters, whereas using finite differences has a cost which increases linearally with the number of parameters. Numerical experiments used the reconstructed NGsim data and the optimal velocity model.

There are two main directions for future research. In this paper we consider only the calibration of a car following model, while a full microsimulator has several other important modules, such as route choice or lane changing decision models. The question of how to apply the adjoint method to those other components can be considered in the future so that the methodology could be applied to a full microsimulation model. There are also more questions regarding the calibration of car following models. Future research can consider the effect of a larger platoon size on the calibration problem, so that multiple vehicles (which share leaders) are calibrated at a time. Another interesting question is using the new highly accurate vehicle trajectory data ([18]) to compare and validate different car following models. Since it was found during this study that lane changing vehicle trajectories are not described well by car following models, another question is how to incorporate lane changing dynamics into an arbitrary car following model, and the authors have recently proposed a method for doing so in [76].

Appendix A - Proof of Theorem 1 and Corollary 1

Definition 1**.**

We say that the function $h(y(t),z(t))$ is piecewise continuous on an interval $\mathcal{I}$ if there are a finite number of points $t\in\mathcal{I}$ where $h(y(t),z(t))$ is not continuous.

Definition 2**.**

(Existence and Uniqueness Theorem.) Define the initial value problem

[TABLE]

where $h(t,y(t))$ is a continuous function on the region $\mathcal{R}=\{(t,x):|t-t_{0}|<s,|x-x_{0}|<y\}$, and $h$ is uniformly Lipschitz in $x$, meaning

[TABLE]

for some constant $C$ in the region $\mathcal{R}$. Then the solution $x(t)$ exists and is unique on the interval $|t-t_{0}|<s*$ for some $0<s^{*}<s$.

Proof 1** (Proof of Theorem 1.).**

*We say that there is a discontinuity for vehicle $i$, either in its car following model or its adjoint variable, if at time $t$: a) its lead vehicle changes, b) its model regime changes, c) the vehicle changes from the car following model to the downstream boundary condition, or d) the loss function $f$ is not continuous at time $t$. We will refer to these different causes of discontinuities as types a) - d) respectively. Define the times $\theta_{j}$ for $j=0,\ldots,m$ as the set of all times when when any vehicle $x_{i},i=1,\ldots,n$ experiences a discontinuity. The set $\{\theta_{j}\}$ is of measure zero because each vehicle can contribute only a finite number of discontinuities. The assumption that $f$ and $h_{i}$ are piecewise continuous and switch regimes only a finite number of times means a vehicle can contribute only a finite number of type b) or d) discontinuities. A vehicle switches to the downstream boundary only once, contributing a single type c) discontinuity. And based on the physical nature of the problem it can be assumed that a vehicle can change its lead vehicle only a finite number of times in a finite interval.

Some of the times $\theta_{j}$ may depend on the switching condition $g$ in Eq. (8). These times may therefore depend on the parameters $p$, and we use the notation $\theta_{j}(p)$ to refer to the switching times for parameters $p$. Besides the set $\{\theta_{j}\}$ being finite, we also require that*

[TABLE]

*meaning that the times of discontinuities must be continuous with respect to the parameters.

Let the set $\{\theta_{j}\}$ be ordered in increasing time,and let it include the first and last times of the simulation $\theta_{0}=\min\{t_{1},\ldots,t_{n}\}$ and $\theta_{m}=\max\{T_{1},\ldots,T_{n}\}$. Since all vehicles have the same lead vehicle on each time interval $[\theta_{j},\theta_{j+1}]$, on that interval there exists some ordered set of vehicles $\chi_{j}$ which contains all vehicles which are simulated in that time interval, ordered such that the first vehicle has no simulated vehicles as leaders (i.e. the vehicle farthest downstream) and the last vehicle has no simulated vehicles as followers (i.e. the vehicle farthest upstream). Following the definition it is true that for $i\in\chi_{j}$ and $t\leq T_{i-1}$ all the partial derivatives $\partial f/\partial x_{i}$, $\partial h_{i}/\partial x_{i}$, $\partial h_{i}/\partial p$, $\partial h_{i}/\partial x_{L(i)}$ as well as $f$ and $h_{i}$ are continuous with respect to their arguments, except possibly at the endpoints.

First we show that $x_{i}(t),i\in\chi_{0}$ is continuous on the interval $[\theta_{0},\theta_{1}]$. This is guaranteed by applying the existence and uniqueness theorem in order to each vehicle. Since by construction $h_{i}(x_{i},x_{L(i)},p_{i})$ is continuous in this interval excluding the right endpoint, and $\partial h_{i}/\partial x_{i}$ is bounded, then the initial value problem $\dot{x}_{i}=h_{i}$ with initial condition $x_{i}(t_{i})=\hat{x}_{i}(t_{i})$ has a solution, namely

[TABLE]

*which is valid for $t_{i}\leq t<\theta_{1}$. Any discontinuity which may exist at $\theta_{1}$ is not felt by the integral since the point has measure 0. Since we solve for each vehicle following the order $\chi_{0}$, the continuity of $x_{L(i)}$ for $L(i)\in\chi_{0}$ is already obtained by the time it is needed to compute its follower.

For some arbitrary interval $[\theta_{j},\theta_{j+1}]$, where the previous interval is already solved, take the left sided limit of the previous interval as the initial condition, i.e. $x_{i}(\theta_{j})=\underset{\epsilon\rightarrow 0}{\lim}\ x_{i}(\theta_{j}-\epsilon)$. This point will always be defined by the previous interval. Then the solution on the new interval is defined as*

[TABLE]

*Note that any possible discontinuities in $h_{i}$ at times $\theta_{j}$ and $\theta_{j+1}$ do not affect the integral because they are of measure 0. Applying the existence and uniqueness theorem again, we have that $x_{i}(t)$ is continuous on $[\theta_{j},\theta_{j+1}]$. Our choice of initial conditions guarentees continuity on the larger interval $[\theta_{j-1},\theta_{j+1}]$. Thus it is clear from an induction argument that $x_{i}(t)$ will be continuous on the interval $[t_{i},\theta_{j+1}]$, despite any possible discontinuities in $h_{i}$.

The exception to the continuity of $x_{i}$ is at time $T_{i-1}$, where the downstream boundary conditions will cause a jump discontinuity. Upon reaching the $\theta_{j}=T_{i-1}$, the remaining trajectory of $x_{i}$ is to be taken as*

[TABLE]

*which guarantees continuity only in position (i.e. $x_{i}^{*}$, the first component of vector $x_{i}$) and the speed $v_{i}$ will have a jump discontinuity.

Thus we proven that the simulated trajectories $x_{i}$ are continuous on the interval $[t_{i},T_{i-1}]$ (excluding $T_{i-1}$) by repeatedly applying the existence and uniqueness theorem on the intervals $[\theta_{j},\theta_{j+1}]$ for an arbitrary $j$.

Then for the objective*

[TABLE]

*since $f$ is assumed to have a finite number of discontinuities in $t$, $\hat{x}_{i}$ is assumed to be continuous, and we have shown $x_{i}(t)$ is continuous on the interval $t_{i}\leq t<T_{i-1}$, it follows that $F$ is continuous.

To show that $\lambda_{i}(t)$ is continuous, apply the existence and uniqueness theorem starting with the last interval $[\theta_{m},\theta_{m-1}]$, going backwards over vehicles on $\chi_{m-1}$. Define*

[TABLE]

and for an arbitrary interval $[\theta_{j},\theta_{j-1}]$ define

[TABLE]

*where additionally for the first interval for $\lambda_{i}$ (the interval where $\theta_{j}=T_{i}$), the initial condition should be $\lambda_{i}(\theta_{j})=0$. The justification for why $\lambda_{i}(t)$ is continuous follows the same argument as for $x_{i}$. The Lipschitz condition for $\lambda_{i}$ necessary to apply the existence and uniqueness theorem for $\lambda_{i}(t)$ is the same as the Lipschitz condition for $x_{i}$ ($\partial h_{i}/\partial x_{i}$ is bounded). Additionally, by assumption we have all the partial derivatives of $h_{i}$ are continuous on each interval excluding the endpoints. Since the adjoint system contains the term $\mathbbm{1}(G(i,t)\neq 0)\lambda_{G(i)}^{T}(t)\dfrac{\partial h_{G(i)}}{\partial x_{i}}$ the $\lambda_{i}$ need to be solved backwards in the order $\chi_{j}$ so that the follower’s adjoint variable is already proven as continuous by the time it is needed. Also note that $\dfrac{\partial h_{G(i)}}{\partial x_{i}}$ is simply $\partial h_{k}/\partial x_{L(k)}$ where $k=G(i)$.

Because $\lambda_{i}$ is continuous on $[t_{i},T_{i}]$, $x_{i}$ is continuous on $[t_{i},T_{i-1})$, and $\partial h_{i}/\partial p$ contains a finite number of discontinuities, it follows that $\lambda_{i}^{T}(t)\partial h_{i}/\partial p$ is a piecewise continuous function with a finite number of discontinuities. Then the gradient*

[TABLE]

*is continuous since the number of discontinuities over the integral is of measure 0. ∎

Proof 2** (Proof of Corrolary 1.).**

Showing $F$ is lipschitz continuous is equivalent to showing that $dF/dp$ is bounded. To do this define

[TABLE]

Then by definition we have

[TABLE]

and by Gronwall’s inequality in integral form

[TABLE]

*To show that all $\lambda_{i}(t)$ are bounded, prove that each $\lambda_{i},i\in\chi_{j}$ are bounded in the interval $[\theta_{j},\theta_{j-1}]$ going backwards over the order specified by $\chi_{j-1}$, starting with $j=m$. Since we assumed that $\partial f/\partial x_{i}$ and $\partial h_{i}/\partial x_{i}$ are both bounded, and for the last vehicle $k$ in $\chi_{j}$, $G(i,t)=0$, then it follows from Eq. (18) that $\lambda_{k}(t)$ is bounded in the interval $[\theta_{j},\theta_{j-1}]$. Then after showing this for $k$, it can be shown for the other vehicles in $\chi_{j}$. Repeating this process for each interval shows that each $\lambda_{i}(t)$ is bounded on $[T_{i},t_{i}]$.

Then for the gradient Eq. (15) since we showed $\lambda_{i}(t)$ is bounded, and $\partial h_{i}/\partial p_{i}$ is bounded by assumption, then it follows that $dF/dp$ is also bounded, and hence $F$ is Lipschitz continuous. ∎*

A.1 - Example with OVM

Written as a first order model OVM is

[TABLE]

Where $p_{i}=(c_{1},c_{2},c_{3},c_{4},c_{5})>\textbf{0}$ and $s_{i}(t)=x_{i}^{*}(t)-x_{L(i)}(t)-l$ (where $l$ is the length of the lead vehicle). Since it is a single regime model, it will be everwhere continuous. The partial derivative

[TABLE]

so it can be bounded by a constant for an arbitrary vector space. The other partial derivatives can be bounded as

[TABLE]

And so the objective will be Lipschitz continuous assuming that the headway $s_{i}(t)$ and velocity $v_{i}(t)$ are finite.

Appendix B - Supplemental figures/tables for comparison of algorithms

Appendix C - Adjoint method for gradient with DDE model

We will now deal with calculating the gradient for Eq. (6). As before, if not explicitly stated, quantities can be assumed to be evaluated at time $t$. Define the augmented objective function:

[TABLE]

and since $F=L$

[TABLE]

where $S(i)$ is defined as the set of all followers for vehicle $i$ and we have chosen $\lambda_{i}^{T}(T_{i-1}^{*})=0$. For a DDE, use a change of coordinates on terms containing $\partial x_{i}(t-\tau)/\partial p$.

[TABLE]

Where we make use of the fact that $\partial x_{i}(t)/\partial p=0$ for $t\in[t_{i},t_{i}+\tau_{i}]$ since the initial history function doesn’t depend on the parameters. Now substituting Eq. $\eqref{29}$ into Eq. $\eqref{28}$ and grouping terms, one obtains the adjoint system:

[TABLE]

with initial conditions $\lambda_{i}(T_{i})=0$. For a DDE model, the adjoint system also becomes a DDE. After solving for the adjoint variables, one can compute the gradient:

[TABLE]

Appendix D - Extension of adjoint method for Hessian

In this section, for a scalar being differentiated with respect to two vectors, the result is a matrix where the rows correspond to the first vector, and the columns correspond to the second vector. For example,

[TABLE]

For a matrix being differentiated to two vectors, the result is a 3-dimensional array which should be interpreted as a vector of matrices. For example

[TABLE]

where $h_{i}^{1}$ and $h_{i}^{2}$ refer to the two components of $h_{i}$. So the normal rules of linear algebra can be applied, e.g.

[TABLE]

where $\lambda_{i}^{1}$ and $\lambda_{i}^{2}$ are the two components of the adjoint variable.

The adjoint method extends in a straightforward way for computing the Hessian. Using finite differences would have complexity $\mathcal{O}(m^{2}T(n))$, whereas the adjoint method can be shown to have complexity $\mathcal{O}(mT(n))$ ([77]).

When calculating the gradient, the benefit of the adjoint method is that you avoid calculating $\partial x_{i}/\partial p$ terms. When calculating the Hessian, the adjoint method will be used to avoid calculating $\partial^{2}x_{i}/\partial p^{2}$ terms. We will still be left with the unknowns $\partial x_{i}/\partial p$, which need to be calculated using the “variational” approach, also known as the “forward sensitivity” approach (see [64, 63]).

To obtain the variational system, differentiate the model Eq. (3) with respect to the parameters

[TABLE]

and then switch the order of differentiation on the left hand side to obtain the variational system

[TABLE]

Eq. $\eqref{24}$ gives ${\partial x_{i}(t)}/{\partial p}$ starting from the initial conditions ${\partial x_{i}(t_{i})}/{\partial p}$. Note that the variational system for $x_{i}$ is only defined up to time $T_{i-1}$ because $\partial x_{i}(t)/\partial p=\partial x_{i}(T_{i-1})/\partial p$ for $t\in[T_{i-1},T_{i}]$. Now starting from Eq. (4.1), differentiate twice instead of once

[TABLE]

The adjoint variables are defined the same way whether the gradient or Hessian is being computed. Eq. (15) gives the adjoint system for the unmodified objective $F$. After solving for all the adjoint variables, in addition to solving the variational system Eq. (24), one can calculate the Hessian:

[TABLE]