Fundamental diagram estimation by using trajectories of probe vehicles

The fundamental diagram (FD), also known as the flow--density relation, is one of the most fundamental concepts in the traffic flow theory. It describes the relation between equilibrated flow, density, and speed in traffic flow. Conventionally, FDs are estimated by using fixed roadside sensors such as detectors. On the other hand, there is a paucity of studies on FD estimation by using mobile sensors, namely probe vehicles. Probe vehicles enable the collection of traffic data from the entire road network, and thus they are useful sensors for large-scale traffic management. In this study, a framework of FD estimation by using probe vehicle data is developed. It determines FD parameters based on trajectories of randomly sampled vehicles and a given jam density that is easily inferred by other data sources. Specifically, an algorithm for estimating a triangular FD based on actual, potentially noisy traffic data is developed. The algorithm is empirically validated by using real-world probe vehicle data on a highway. The results suggest that the algorithm accurately and robustly estimates the FD parameters.


Introduction
The fundamental diagram (FD), also known as the flow-density relation, is literally one of the most fundamental concepts in the traffic flow theory.An FD describes the relation between flow and density in steady traffic (that is sometimes referred as equilibrium or stationary traffic), i.e., traffic in which all the vehicles exhibit the same constant speed and spacing (Daganzo, 1997;Treiber and Kesting, 2013;Yan et al., 2018).In theory, an FD itself contains useful information on traffic features, such as the value of free-flow speed and flow capacity, and distinction between free-flow and congested regimes.Observational studies also indicated a clear relation between flow and density in actual traffic at near-steady states (Cassidy, 1998).Additionally, macroscopic traffic flow dynamics are modeled by combining an FD and other principles-the most widely known example is the Lighthill-Whitham-Richards (LWR) model (Lighthill and Whitham, 1955;Richards, 1956).Furthermore, FDs explain microscopic vehicle behavior to a certain extent (Newell, 2002;Duret et al., 2008;Jabari et al., 2014).Thus, FDs are utilized by a wide variety of academic and practical purposes in traffic science and engineering fields such as traffic flow modeling/analysis/simulation (Daganzo, 1994(Daganzo, , 1997(Daganzo, , 2006;;Treiber and Kesting, 2013), traffic control (Papageorgiou et al., 2003), traffic state estimation (Seo et al., 2017a), and dynamic traffic assignment in network (Szeto and Lo, 2006).
In order to estimate a parametric FD in actual traffic, it is necessary to collect data from traffic, assume the functional form of its FD, and estimate the FD parameters by fitting the curve to the data in general.Typically, FDs are estimated by using roadside sensors (e.g., cameras and detectors) from the era of Greenshields (1935).This is a straightforward method because usual roadside sensors measure traffic count and occupancy that are closely related to flow and density, respectively, at their location.
The limitation of the roadside sensor-based FD estimation is evident: it is unable to estimate FDs where sensors are not installed.Thus, it is difficult to determine FDs on roads without sensors such as arterial roads.It is also difficult to identify the exact locations and characteristics of bottlenecks even on freeways with some sensors.The issue poses substantial difficulties for practical implementations of the aforementioned traffic science and engineering works because the installation and operational costs of roadside sensors are generally high.
Probe vehicles received increasing attention currently (Sanwal and Walrand, 1995;Zito et al., 1995;Herrera et al., 2010;Shladover, 2017). 1 A well-known advantage of probe vehicles as a traffic monitoring tool is their significantly broad data collection range in the spatiotemporal domain when compared with roadside sensors.For example, one of the most typical applications is traffic state estimation using partial measurement data (Nanthawichit et al., 2003;Seo et al., 2017a); probe vehicles enabled the estimation of a traffic state in a wide spatiotemporal domain.However, existing studies on these types of applications of probe vehicles often assume their FDs exogenously to relate flow and density to speed measured by probe vehicles (e.g., Nanthawichit et al., 2003;Yuan et al., 2012;Wada et al., 2015;Canepa and Claudel, 2017).This is a critical limitation for application of probe vehicles because FD estimation itself is a challenging problem as discussed above.Especially, there is a paucity of research on FD estimation by using probe vehicle data.Thus, systematic and computational approaches for FD estimation are desirable given the high availability of probe vehicle data.
The aim of the study is to propose methodology to estimate FDs by using probe vehicle data.In order to enable the estimation, we allow the methodology to rely on minimum exogenous assumptions such as FD's functional form and value of a parameter of the FD (i.e., the jam density that is inferred from external knowledge).Moreover, a computable algorithm for FD estimation is developed and validated by using real-world data to illustrate the empirical validity of the proposed methodology.
The rest of the paper is organized as follows.The literature on FD estimation is reviewed and originality of this study is highlighted in Section 2. Identifiability of an FD by using probe vehicle data under idealized conditions is rigorously established in Section 3. Specifically, sufficient conditions with respect to input data, FD's functional form, and exogenous assumptions to identify an FD are discussed for general cases and explained in detail for a triangular FD case.An algorithm for estimating a triangular FD by using actual, noisy traffic data is developed in Section 4. The proposed algorithm is validated by applying it to real-world probe vehicle data in Section 5. Section 6 concludes this paper.The frequently used abbreviations and notations are listed in Appendix A.

Literature review on FD estimation
FD estimation by using detectors is extensively studied since 1935.Many approaches use traffic data aggregated for a short time period (e.g., traffic count for 30 sec-5 min periods) as inputs of estimation and extract statistical relation from them (see Treiber and Kesting, 2013;Qu et al., 2017;Knoop and Daamen, 2017, and references therein).Other approaches use disaggregated data that distinguish each individual vehicle in measurement area (e.g., Duret et al., 2008;Chiabaut et al., 2009;Coifman, 2014Coifman, , 2015;;Jin et al., 2015).The latter approaches enable the precise identification of steady traffic states.These methodologies rely on the fact that flow and density or trajectories of all the vehicles are easily obtained by detectors.However, this is not the case for conventional probe vehicles: only trajectories of randomly sampled vehicles are collected by probe vehicles.Thus, these methodologies are not applicable for FD estimation by using probe vehicle data.
There is a paucity of research on FD estimation by using probe vehicle data.A few probe vehicle-based traffic state estimation methods simultaneously estimate FDs and traffic state (e.g., Tamiya and Seo, 2002;Seo et al., 2015;Sun et al., 2017).However, they rely on additional data sources, such as detector measurement or vehicle spacing in addition to probe vehicle's positioning.Appendix in Herrera et al. (2010) mentioned an idea for manual inference of an FD by using probe vehicle data under special conditions.However, it was not formulated and not computable.Their idea involves manually detecting a shockwave between a saturated free-flowing traffic and congested traffic, deriving the backward wave speed of a triangular FD by manually measuring the speed of the shockwave, and subsequently calculating the other FD parameters by using a given jam density.However, detection of a shockwave, detection of saturated traffic, and identification of the speed of shockwave are not trivial tasks in actual traffic.

Preliminaries 3.1.1. Assumptions
We assume that traffic dynamics are described by the LWR model with a homogeneous FD and first-in first-out (FIFO) condition without a source or sink.It is expressed as follows: where q(t, x) and k(t, x) denote flow and density, respectively, at time-space point (t, x), Q denotes a flow-density FD, and θ denotes a vector representing the FD parameters (Lighthill and Whitham, 1955;Richards, 1956).
The functional form of an FD Q(k; θ) is assumed as a single-valued, continuous, and piecewise differential.The functional form and parameter values of each of the piecewise differential parts of an FD are assumed as not identical to each other.Additionally, Q(0; θ) = 0 is assumed to hold, which is obviously reasonable.Note that many of commonly used FD functional forms, such as the triangular FD (Newell, 1993), satisfy these assumptions.
A probe vehicle dataset is assumed to contain continuous trajectories of randomly sampled vehicles without any errors.It is represented as {X(t, n) | ∀n ∈ P}, where X(t, n) denotes the position that vehicle n exists at time t, and P denotes a set of all the probe vehicles.The penetration rate of probe vehicles is unknown.

Definitions
Let traffic state in a time-space region be defined by Edie's generalized definition (Edie, 1963), namely where A denotes a time-space region, and q(A), k(A), and v(A) represent flow, density, and speed, respectively, in the time-space region A, N(A) denotes a set of all the vehicle in A, d n (A) and t n (A) denote distance traveled and time spent, respectively, by vehicle n in A, and |A| denotes the area of A. If traffic in A is steady, is satisfied by the definition of the LWR model (2).Without loss of generality, the FD is defined as follows.The total number of parameters of the FD is denoted by g.The number of the piecewise differential parts of the FD is h, and its j-th part (labeled by the increasing order of k) includes g j parameters and the function is denoted by Q j (hereafter, a piecewise differential part is simply referred to as a part).Thus, h j=1 g j = g = |θ| holds.The h-th part includes at least one parameter, that is, the jam density κ.An example of the FD is illustrated in Fig. 1.
A traffic state (q, k) is said to belong to j-th part of an FD if k ∈ (k min j , k max j ) holds, where k min j denotes the minimum or lower-bound density of j-th part and k max j denotes the maximum or upper-bound density of j-th part.We assume that an observed traffic state in j-th part does not belong to the extrapolation of (j − 1)-th part for the purpose of simplicity.2Density Flow jam density κ part: j = 1 Let us consider a probe pair, that is, a pair of arbitrary two probe vehicles.A pair is denoted by m, and the preceding probe vehicle in the pair m is denoted by m − whereas the following probe vehicle is denoted by m + .Note that arbitrary and unknown number of vehicles may exist between a probe pair.
Let (q m (A), k m (A)) be probe traffic state (PTS) of probe pair m defined as following Edie's definition.Similarly, probe-speed is also defined as The PTS denotes flow and density of probe vehicle m + only in time-space region A. If (i) traffic in A is steady, (ii) distance traveled by each vehicle in A is equal to the others, and (iii) probe vehicle m + exists in A, then the following relations are satisfied: where |N(A)| denotes the size of N(A) (i.e., number of the vehicles in A).
Let c m − 1 denote the number of vehicles between probe pair m, and a m denote a time-space region whose boundary is on trajectories of m + and m − (the region includes m + but not m − ).Since source/sink and FIFO-violations are absent, c m is constant and always holds irrespective of shape, time, and location of a m .If (i) traffic in a m is steady and (ii) distance traveled by each vehicle in a m is equal to the others, then is satisfied from Eqs. ( 4), (7a), and (7b).Hereafter, a PTS of pair m under the aforementioned conditions (i) and (ii) is simply referred to as steady PTS.Equation ( 9) means that steady PTSs follow a relation that is similar (in the geometric sense) to the FD with a scale of 1/c m on a flow-density plane.Hereafter, this relation is referred to as probe fundamental diagram (PFD) of m.Note that the PFD of probe pair m is also expressed as follows: where θ m denotes a pair-specific parameter vector, because the FD and PFD are similar.
In order to ensure that a m easily satisfies the aforementioned conditions (i) and (ii), we specify a m as follows: where i denotes an index number, τ i denotes a reference time for a mi , and φ and ∆t are the given parameters.An arbitrary value is acceptable for τ i and φ, and an arbitrary non-zero value is acceptable for ∆t if the corresponding trajectory X(t, n) exists.The shape of a m is illustrated in Fig. 2. The specification (11) guarantees that if traffic in a m is steady, then distance traveled by each vehicle in a m is equal to the others, which is d m + (a mi ).This feature ensures that the following discussion is concise without sacrificing significant generality.
An FD is said to be identifiable by using a given dataset (a set of traffic states or a set of PTSs) if the value of its parameter vector θ is determined based on the dataset.Similarly, a PFD of certain probe pair m is said to be identifiable by using the given dataset (a set of PTSs of m) if the value of its parameter vector θ m is determined based on the dataset.
For the readers convenience, important abbreviations and notations used throughout this paper are summarized in Appendix A.
Figure 2: Shape of the time-space region a m .

Identifiability of FD based on probe vehicle data
Given the aforementioned assumptions, a sufficient condition for the identification of an FD is stated as follows.
Proposition 1.An FD with a known functional form is identifiable if g j or more different steady traffic state data (q i , k i ) belonging to j-th part of the FD are available for all j.
Proof.First, let us consider determination of belongingness of traffic state data to j-th part, when belongingness of data to l-th part for all l < j is determined if j > 1.The shape of j-th part can be uniquely determined using first g j data, as Q j is assumed as differential.Once the j-th part is determined, the belongingness of other data belonging to j-th part (if any) are also resolved, as Q j is assumed as differential, Q l (∀l = j) is assumed as not identical to Q j , and we assumed that a traffic state datum in (j + 1)-th part does not belong to the extrapolation of j-th part.As a consequence, if j < h, we know that the data whose belongingness is not determined are belonging to (j + 1)-th or later parts; if j = h, the belongingness of all the data is determined.
By executing the above procedure with j := 1, and subsequently iteratively execute it with j := j + 1 until j = h, the belongingness of all the data is determined.Now, the parameters of the function Q j can be determined based on the known functional form and the given traffic state data for j-th part for all j, since g j or more different data are available for each j-th part.This is directly related to conventional FD estimation by using detector data.Note that, in reality, it is difficult to rigorously assess whether the condition-g j or more different steady traffic state data (q i , k i ) belonging to j-th part of the FD are available for all j-is actually satisfied based on an available dataset.However, the condition is likely to be satisfied when abundant data were collected, and it is easy to collect these types of data in today's data age.In this proof, we used the assumption that a traffic state datum in (j + 1)-th part does not belong to the extrapolation of j-th part.This is actually a fairly reasonable assumption, because Q j is not identical to Q j+1 and it is very rare to observe a data belonging an extrapolation of the other parts.
By substituting "FD" with "PFD of a probe pair m" and "traffic state (q i , k i )" with "PTS (q mi , k mi )" in Proposition 1, the sufficient condition for the identification of a PFD by using probe vehicle data is stated as follows.
Corollary 1.An PFD of a probe pair m with a known functional form is identifiable if g j or more different steady PTS data (q mi , k mi ) belonging to j-th part of PFD are available for all j.
Now, we consider the identifiability of an FD when a corresponding PFD is identified.
Proposition 2. An FD with a known functional form is identifiable if (i) a corresponding PFD of a probe pair is identified and (ii) the value of the jam density is known.
Proof.Since an PFD of a pair, say m, is identified, the jam density of the PFD, κ/c m , is known.Thus, the value of c m is determined by comparing κ/c m to the known jam density κ.Thus, the other FD parameters are determined based on Eq. ( 9).
In fact, this FD parameter identification from known PFD parameters is very simple-a PFD is similar (in the geometric sense) to an FD with a scale of 1/c m as shown in Eq. ( 9); thus, the FD parameters are simply obtained once c m and the PFD parameters are known.
From Corollary 1 and Proposition 2, we establish the following corollary on the identifiability of FD based on probe vehicle data: Corollary 2. An FD with known functional form is identifiable if (i) g j or more different steady PTS data belonging to j-th part of the FD are available for all j for a probe pair and (ii) the value of the jam density is known.
Besides, even if the value of the jam density is not known, some of the FD parameters are determined by following a similar approach.Proposition 3. It is possible to determine parameters of an FD without the number-of-vehicle dimensions if g j or more different steady PTS data belonging to j-th part of the FD are available for all j for a probe pair.
Proof.Given this condition, the PFD of a probe pair, say m, is identified.Since a PFD is similar (in the geometric sense) to an FD in a flow-density plane, only parameters with number-of-vehicle dimensions differ between the PFD and FD.Thus, parameters of the FD without number-of-vehicle dimensions are the same as those of the PFD.
A typical example of such determinable parameters is the free-flow speed.
In summary, an FD with known functional form is identifiable if probe vehicles traveled through traffic with sufficiently wide variety of regimes (density) and the jam density is known.The first condition is realized when probe vehicle data were collected for a long time.The second condition is realized by several methods as discussed later in Section 3.4.The probe vehicle data determine the shape of an FD while the jam density determines the size of an FD.

Geometric interpretation with triangular FD case
A detailed explanation is given with a case with the triangular FD (Newell, 1993).A triangular FD is defined as where u denotes the free-flow speed and w denotes the backward wave speed.Equation (12a) represents a free-flowing regime, whereas Equation (12b) represents a congested regime.In this case, θ consists of u, w, and κ; and g = 3, h = 2, g 1 = 1, and g 2 = 2 hold.Note that it is not necessary for an analyst to know whether a specific PTS belongs to free-flowing or congested regime.The belongingness is determined as a result of the identification of the FD.Geometric interpretation of the triangular FD identification is explained using Fig. 3 as follows.As explained previously, traffic state S m (a mi ) follows the FD if region a mi is steady.Similarly, PTS S m (a mi ) is equal to S m (a mi )/c m and follows the PFD if region a mi is steady.In Fig. 3, regions a m1 and a m4 are steady; and thus their PTS are on the PFD.Conversely, regions a m2 and a m3 are non-steady; and thus their PTS are not on the PFD.Steadiness of region a mi can be inferred by using trajectories of probe vehicles.The necessary conditions for the steadiness are that the two probe vehicles' speeds are time-invariant, and they are equal to each other.
Suppose that the probe vehicles m − and m + traveled to another congested region, say a m5 (the gray dot in Fig. 3).In this case, the PFD is uniquely determined by selecting values of u, w, and c m such that all the points S m (a m1 ), S(a m4 ), and S(a m5 ) are on the PFD q = Q(k; u, w, κ/c m ).In other words, a triangle that satisfy following conditions can be determined uniquely: its vertex is at point (0, 0), one of its edges is on line q = 0, and its another two edges pass the blue point, red point, and gray point.Thus, the FD is identified from the probe vehicle data and the jam density.

Discussion
The proposed method identifies all the parameters of an FD with the exception of its jam density by using probe vehicle trajectory data.The value of the jam density is required as known prior to the FD identification.There are several methods to obtain the value: common knowledge and remote sensing would be especially useful at this moment.
The jam density is considerably time-and space-independent variable when compared to the other FD parameters such as capacity.Thus, it is reasonable to set the jam density based on a value that is observed at nearby locations.Alternatively, remote sensing (e.g., image recognition) from satellites is used to measure the jam density directly.Although remote sensing does not provide us time-continuous traffic state variables, such as flow and speed (because the time interval between two measurements is extremely long, few hours to few days at least), it measures space-continuous density accurately (McCord et al., 2003).Thus, the jam density could be inferred as an upper limit of such measured density.
Besides, for the triangular FD case, the proposed method identifies the values of the free-flow speed and the backward wave speed even if the value of the jam density is unknown as stated in Corollary 4. Since these two variables are important to characterize traffic flow, this property would be useful for various purposes such as traffic speed reconstruction and travel time estimation (Treiber et al., 2011).
A limitation of the proposed method is that a homogeneous FD is assumed.It implies that the method in itself is unable to consider a bottleneck.This issue could be resolved by extending the method: for example, it is possible to estimate section-dependent FDs by applying the proposed method to each section iteratively.Another limitation is that the proposed method does not consider multi-valued and/or non-continuous FDs such as those involving a capacity drop.

Estimation algorithm for actual traffic
Actual traffic data are not exactly consistent with the theoretical LWR model, and thus it is not possible to directly apply the methodology described in Section 3 to actual data-all the data will be considered as non-steady.The inconsistency is mainly due to two sources.First, various traffic phenomena, such as the existence of non-equilibrium flow, bounded acceleration, heterogeneity among vehicles, FIFO violation in multi-lane sections, and the resulting stochasticity of macroscopic features, are not considered by the LWR model (Jabari et al., 2014;Jin et al., 2015;Coifman, 2015;Qu et al., 2017;Jin and Laval, 2018).Second, data always include measurement noise.
In order to account for this challenge, this section presents an estimation algorithm for a triangular FD (Newell, 1993) from actual traffic data using maximum likelihood estimation.Specifically, a method of extracting near-steady traffic state data from probe vehicle data is presented in Section 4.1.Then, an estimation method for free-flow speed and backward wave speed based on near-steady traffic state data is presented in Sections 4.2 and 4.3.The triangular FD is employed owing to its simplicity, good physical interpretation of parameters, and popularity in various applications (e.g., Daganzo, 2006).
In order to simplify the mathematical representation, the triangular FD is characterized by free-flow speed, u, backward wave speed, w, and, y-intercept of the congestion part of FD, α as illustrated in Fig. 4. Note that α ≡ κw holds, and thus, the characterization does not contradict the assumptions in Section 3.1.

Extraction of near-steady PTSs
Suppose that multiple PTSs S m (a mi ) for multiple m and i were calculated from probe vehicle data based on the definition ( 5) and ( 11).This set of PTSs is referred to as raw PTS set.
It is necessary to extract PTSs that are near-steady, meaning that they nearly satisfy the conditions in Corollary 3 from the raw PTS set.To do this, two filtering methods are employed.The first method is based on the linearity of the trajectories of two probe vehicles in pair m.Specifically, S m (a mi ) is considered as near-steady if trajectories of pair m in region a mi satisfy a condition based on the coefficient of variation: where σ m (a mi ) and vm (a mi ) denote the standard deviation and mean, respectively, of the instantaneous speed of probe vehicles in pair m in region a mi , and θ steady denotes a given threshold.
In the second method, {S m (a mi ) | ∀i} of each m is verified as to whether it is likely to satisfy the conditions in Corollary 3. Data from the probe pair m are discarded if {S m (a mi ) | ∀i} included no S m (a mi ) with speed lower than u min or no S m (a mi ) with speed faster than u min , where u min denotes a given lower bound for free-flow speed.3This corresponds to conditions (i) and (ii) in Corollary 3. Additionally, data from probe pair m are discarded if correlation between probe-flow and probe-density in {S m (a mi ) | ∀i, v(a mi ) ≤ u min } exceeds θ corr , or standard deviation in {v m (a mi ) | ∀i, v(a mi ) ≤ u min } is smaller than θ std , where θ corr and θ std denotes given thresholds.This corresponds to condition (ii) in Corollary 3.
Hereafter, a set of PTS data obtained by the above procedure is referred to as near-steady PTS set.A set of probe pairs that collected elements in the near-steady PTS set is denoted as M , and a near-steady PTS set obtained by pair m is denoted as I m .

Backward wave speed estimation
Backward wave speed is first estimated separately from the free-flow speed.This is because it is easy to determine whether a PTS datum definitely belongs to congested regimes based on the lower bound of the free-flow speed u min . 4pecifically, the backward wave speed w is estimated by performing a linear regression on (k(a mi ), q(a mi )) for each i ∈ I m , m ∈ M , and v m (a mi ) ≤ u min ; and the final estimate is given as the mean of all the regression results.Evidently, condition (ii) in Corollary 3 is essential to estimate w appropriately by the method.

Free-flow speed estimation by using EM algorithm
The free-flow speed u is subsequently estimated.To do this, an iterative algorithm categorized as expectation-maximization (EM) algorithm (Dempster et al., 1977) is constructed.In this algorithm, free-flow u is estimated in conjunction with the y-intercept value α m and other variables, such as standard deviation of free-flowing and congested part of an FD, based on the near-steady PTS set.

Likelihood function
Assume that a probe-flow of a steady PTS belonging to regime s (free-flowing, s = F, or congested, s = C) follows a normal distribution in which the mean is the true PFD and standard deviation is r m σ s , where σ s denotes a standard deviation of traffic state from an FD, r m denotes a pair-specific normalizing parameter defined as qm /q, qm denotes the mean of probe-flow of pair m, and q denotes the mean of probe-flow of all of the probe pairs.Thus, likelihood of regime s of PTS datum i of pair m is expressed as where N denotes the probability distribution function of a normal distribution.Note that variables u, α m , s, and σ s are unknown, while the rest are known.The likelihood of the mixture distribution for all the near-steady PTS set is given by where π s denotes the contribution of state s in the mixed distribution (Bishop, 2006).The intuitive meaning of L is a likelihood of observing PTS (q mi , k mi ) set when parameters5 π, u, w, {α m }, {σ s } are given.
Our problem is to estimate the latent parameters π, u, {α m }, {σ s } from the observed near-steady PTS (q mi , k mi ) set.This can be solved by applying an EM algorithm.If the complete dataset (Bishop, 2006) is available (i.e., latent variable z mis that is 1 if S mi belongs to regime s, 0 otherwise is given), then Eq. ( 15) is transformed By taking its expectation of logarithm, one obtains6 where p mi (s) denotes an expectation of z mis .The EM algorithm determines the value of latent parameters π, u, w, α, σ, and p by maximizing Eq. ( 17) with an iterative procedure termed as E-steps and M-steps.

E-step
In the E-step, p mi (s) is updated under given π s , u, w, α m , σ.It is expressed as

M-Step
In the M-step, π s , u, {α m }, {σ s } is updated such that the likelihood ( 17) is maximized under given p mi (s).Specifically, it is expressed as max π,u,{αm},{σs} where u min and u max are non-negative constants.In order to solve problem ( 19), the method of Lagrange multiplier is adopted.The Lagrangian for the problem is defined as where each λ with subscript denotes a Lagrangian multiplier corresponding to each constraint in problem (19).By adopting the standard procedure that exploits the Karush-Kuhn-Tucker condition (Karush, 1939;Kuhn and Tucker, 1951), this problem is analytically solved.The solution is summarized as follows: with and with and Notice that the right-hand side of Eq. ( 22) is indefinite if p mi (F) = 0 holds for all m, i; in other words, it is not possible to determine the free-flow speed if data are not observed from the free-flowing regime.This corresponds to condition (i) in Corollary 3.

Summary of the algorithm
Input for the proposed algorithm is as follows: • Probe vehicle trajectory data The entire procedure of the proposed algorithm is summarized as follows: Step 1 Near-steady PTS set construction step Step 1.1 Construct raw PTS set by calculating PTSs for all m + ∈ P, m − = m − − ∆m, ∆m ∈ ∆M , t ∈ T , ∆t ∈ ∆T , and φ ∈ Φ based on the given probe vehicle data and definition ( 5) and ( 11), where ∆M , ∆T , and Φ denote sets of given parameter values.
Step 1.2 Construct near-steady PTS set by removing specific PTSs from the raw PTS set.Specifically, remove non-steady PTSs based on Eq. ( 13) with given θ steady and nonappropriate PTSs that do not satisfy Corollary 3 based on given θ corr and θ std .
Step 2 Backward wave speed estimation step Step 2.1 Estimate the backward wave speed w by performing a linear regression on PTSs with v ≤ u min in the near-steady PTS set.
Step 3 Free-flow speed estimation step Step 3.1 Define the likelihood function (17).Give an initial state on p mi (s).
Step 3.3 Update p mi (s) by executing the E-step.
Step 3.4 Check the convergence: if the change rate of all the parameters is lower than a given threshold ε, then it converges.If it converges, output the current parameter values and halt the algorithm.If not, go to Step 3.2.

Discussion
The proposed method estimates mean free-flow speed and backward wave speed of a triangular FD along with other auxiliary variables.Two of the auxiliary variables correspond to the standard deviation of flow from the mean FD; therefore, the method also captures the stochasticity of an FD to a certain extent.It is possible to develop estimation algorithms for other FD functional forms by using similar methods.
Although the method requires a few input parameters, their values are determined easily (i.e., no need of fine tuning) and the estimation result would be insensitive to most parameters with the exception of θ steady and κ.The parameter θ steady , threshold for near-steady state identification, must be important for the proposed method, since it is relevant to the original definition of the FD.This issue is investigated by performing a sensitivity analysis in Section 5.3.The jam density could be determined reasonably as discussed in Section 3.4; moreover, even if the jam density is not known, the free-flow speed and the backward wave speed are estimated as discussed in Section 3. The other parameters are not essential.Any combinations of ∆m, ∆t, φ, θ corr , and θ std are acceptable as long as it produces sufficient number of near-steady PTSs from a given probe vehicle dataset such that the conditions (i) and (ii) of Corollary 3 are satisfied.It is not necessary for the upper and lower bounds for free-flow speed, u max and u min , to be close to the actual free-flow speed.The upper bound u max is only required to eliminate non-realistic local optima with unrealistically high u (e.g., over 1000 km/h) as discussed subsequently in this section.The lower bound u min is required to determine congested traffic speed that exhibits significantly slower speed than the free-flow speed; therefore low values (e.g., 60 km/h for a highway) are acceptable.
The objective function of the problem ( 16) is not necessarily convex nor unimodal.Thus, local optima that does not represent an FD may exist.For example, if the upper limit for free-flow speed (19b) does not exist, obvious local optima are obtained at u with unrealistically high speeds, such as 1000 km/h, which classifies all the PTSs belonging to congested regime.The constrain (19b) is useful to eliminate such local optima.
The proposed method estimates u and w separately although it is possible to formulate a joint estimation problem of u, w, and other auxiliary variables as a single optimization problem similar to (19) (Seo et al., 2017b;Kawasaki et al., 2017).However, the separate estimation approach is adopted for the sake of empirical performance, because it was found that a joint estimation approach is difficult to solve stably due to numerous number of local optima, originated from fluctuations in actual traffic (Kawasaki et al., 2017).
The relation between the proposed method and the manual inference method proposed by Appendix in Herrera et al. (2010) is as follows.They share a few of the ideas: that is, they both estimate triangular FD parameters, namely u and w, based on probe vehicle data and a given jam density.Meanwhile, essential differences exist.First, the proposed method is based on the result of more general identifiability of an FD as discussed in Section 3. Thus, the proposed method does not rely on the detection of a shockwave nor the detection of saturated free-flowing traffic (which are used by Herrera et al. (2010)), which are special phenomena and difficult to automatically detect by using probe vehicles.Furthermore, the proposed method is a computable algorithm where discrepancies between the LWR theory and actual traffic are considered based on the steadiness of traffic.

Empirical validation
In this section, validation results of the estimation algorithm based on real-world data are presented.

Datasets
We used traffic data collected at approximately 4 km length section in an inter-city highway near Tokyo, Japan (specifically, the section between Kawasaki Junction and Yokohama-Aoba Interchange in Tomei Expressway).In this section, on-ramps, off-ramps, and merge/diverging sections are absent.The number of lanes is three, and the speed regulation was 100 km/h.According to the statistics, congestion occurs frequently during weekends and holidays mainly due to tourists in this section.The data collection duration was weekends and holidays from April 2016 to September 2016 and corresponded to a total of 61 days.
As probe vehicle data, we used probe vehicle data provided by Fujitsu Traffic & Road Data Service Limited.The data were collected by connected vehicles consisting of logistic trucks and vans who were connected to the Internet for fleet management purposes.The data included map-matched location information of each vehicles for every 1 s.Probe vehicle data collected on a day (April 9, 2016) are visualized in Fig. 5 as time-space diagrams; and it indicate that some probe vehicles slowed down due to congestion approximately between 9 and 11 a.m.The number of trips performed by probe vehicles during the aforementioned period was 23 103 trips.This corresponds to mean headway time between two consecutive probe vehicles 3.8 min.However, many of probe vehicles avoided congested time periods (this is confirmed by Fig. 5); therefore, headway time in congested traffic should exceed the value.In order to validate the estimation results, loop detector data provided by Central Nippon Expressway Company Limited were used as a reference.The detector is located in the middle of the section.It measures vehicle type-specific flow and occupancy for every 5 min.Density was obtained from the data by using the method proposed by Cassidy and Coifman (1997).The flow-density plot based on the detector data is shown in Fig. 6.A clear triangular shape is observed.According to the plot, the flow capacity was approximately 1500 veh/h and the jam density was approximately 100 veh/km.These values are slightly lower than those of usual highway.This might be due to the property of drivers during the period in which most of them were potentially unexperienced drivers (i.e., sightseeing travelers in weekends) and/or possible miss detection by the loop detector.Nevertheless, the values of the free-flow speed and the backward wave speed seems reasonable; therefore, this data can be considered as useful to validate the proposed method.A comparison between speed measured by probe vehicles and detector for each time period is shown in Fig. 7. Essentially, the probe vehicle data exhibits a tendency similar to that of the detector data.However, in the high-speed regime, probe vehicle data tend to be slightly slower than the detector data.The bias is potentially because the probe vehicles consisted of logistic trucks and vans that are typically slower than average cars.From these results, we conclude that the probe vehicle data and detector data were sufficiently consistent in terms of speed measurement.

Parameter setting
The time-space regions a were constructed based on all the combinations (∆m, ∆t, φ) of the following values: ∆m ∈ {1, 2, 3, 4}, ∆t ∈ {1, 2, 3, 4} (s), φ ∈ {5, 10, 15, 20, 30} (km/h).The threshold for near-steady state identification, θ steady , was selected from 0.005 to 0.3.We performed a sensitivity analysis on this parameter by selecting each of them.The cleaning parameters for PTSs were set as θ corr = −0.8 and θ std = 1 (km/h).The upper and lower bound for free-flow speed were set as u min = 60 (km/h) and u max = 120 (km/h), respectively, and they are considered as reasonably mild constraints for a highway with 100 km/h speed regulation.The jam density κ was set to 100 veh/km based on Fig. 6.The initial state of the EM algorithm is p mi (F) = 1 and p mi (C) = 0 for m and i with v m (a mi ) ≤ u min , and p mi (F) = 0 and p mi (C) = 1 for m and i with v m (a mi ) > u min .The convergence check parameter of the EM algorithm, ε, was set to 0.01.

Estimation results
First, an estimation result under a moderate setting, θ steady = 0.15, is presented.Given the parameter setting, 1762 probe pairs produced multiple near-steady PTSs that appeared to satisfy the conditions in Corollary 3, and the total number of PTSs was 437 337.The estimated values were free-flow speed u = 81.2(km/h) and backward wave speed w = 17.44 (km/h).The values of estimates are summarized in Tab. 1 along with other parameter settings.
The estimation result with θ steady = 0.15 is illustrated in Fig. 8.The solid red line denotes mean estimates for u and w, the dashed lines represent the mean plus/minus the standard deviations σ F and σ C of the estimates, and the blue dots represent the detector data for reference.The mean estimate indicates good agreement with detector data.Additionally, the ratio of detector measurements that fall into the range between mean plus/minus standard deviation of estimates was 57.2%.The estimated free-flow speed appears slightly slower than that of detector data; and this is potentially due to the biased speed of probe vehicles (Fig. 7).
All the estimation results with different θ steady are summarized in Tab. 1.Based on these results, sensitivity on θ steady is investigated.With respect to the estimated u and w, similar results were obtained in all cases with the exception of θ steady = 0.3 and θ steady = 0.05.The case with θ steady = 0.3 resulted in low u and excessively high w.This is due to the low steadiness in the PTSs filtered by θ steady = 0.3; as a result, many of PTSs covered both the free-flowing and congested regime that are not appropriate for the FD estimation.The case with θ steady = 0.05 resulted in slightly low w; and this potentially suggests that the estimates are slightly unstable if the number of PTSs are small.With respect to the estimated σ s , a clear tendency was obtained: namely, they tended to decrease when the θ steady decreased.This is an expected result since low θ steady tends to eliminate scattered PTSs in the flow-density plane that typically exhibits relatively low steadiness.

Discussion
The results indicated that the proposed method estimates reasonable values of the free-flow speed and the backward wave speed, depending on the value of threshold for the steadiness θ steady .According to Tab. 1, the method derives fairly good estimates when θ steady is less than 0.2, which is a wide range for this type of variable (i.e., coefficient of variation).Conversely, as qualitatively expected, a case with excessively high θ steady results in improper estimates because it did not properly remove non-steady traffic data.
In addition to the mean FD, the proposed method also estimates standard deviation in the FD, implying that it can capture stochasticity of the FD.The standard deviation in free-flow regime is almost always lower than that of the congested regime; this is a reasonable result considering the nature of traffic flow.When compared to the detector data, the cases with θ steady in [0.1, 0.2] derived standard deviations similar to those obtained in detector data.
With respect to the required amount of probe vehicle data, the results suggest that probe vehicle data with a few minutes of mean headway collected for few months are sufficient to estimate the FD.If the penetration rate is lower than this case, the FD estimation is possible by collecting the data for a longer period (as long as probe vehicles were randomly sampled).With respect to the section length required to estimate an FD, the results suggest that the length of a few kilometers is sufficient for the employed data.It is difficult to precisely derive the minimum requirement because it depends on traffic conditions and the amount of data.Generally, if traffic congestion occurs frequently or the amount of probe vehicle data is large, then some of the probe vehicles are likely to satisfy the conditions (i) and (ii) of Corollary 3 (experience several traffic states) in a short road section.For this dataset, the minimum would be less than a kilometer, because some of the probe vehicles experience several traffic states at one kilometer length segment from the downstream end of this section as illustrated in Fig. 5.
The optimization problem may involve unrealistic local optima as discussed in Section 4.5.However, in the empirical analysis, the proposed method always yielded solution with physically reasonable results with the exception of the excessively high θ steady case.This is potentially due to the separated estimation approach and the constraint on free-flow speed that we employed.

Conclusion
A theory and an algorithm for fundamental diagram (FD) estimation by using trajectories of probe vehicles is presented.They estimate FD parameters (e.g., free-flow speed and backward wave speed in the case with a triangular FD) based on the jam density and time-continuous trajectories of randomly sampled vehicles.The algorithm is developed to estimate a reasonable FD from actual noisy traffic data.The algorithm was empirically and quantitatively validated based on real-world probe vehicle data on a highway.The results suggested that the algorithm accurately and robustly estimates triangular FD parameters and captures the stochasticity of the FD.
Two future research directions are especially considerable.The first involves the development of a method that considers spatially heterogeneous FDs such that a bottleneck is endogenously detected.This is possible by extending the method to estimate location-dependent u and w and global c m .The second involves the incorporation of a jam density inference method based on remote sensing by using satellites (c.f., Section 3.4), in order to make the method independent of exogenous knowledge.For this purpose, the "small satellites" that attracted significant attention recently (Sandau, 2010) will be useful owing to their flexible and low-cost sensing capability.The latter direction is now being investigated by the authors (Seo and Kusakabe, 2018, submitted).Threshold for near-steady traffic detection i Index for datum a mi i-th time-space region between probe vehicles in pair m q m (a mi ) Probe-flow of pair m in region a mi k m (a mi ) Probe-density of pair m in region a mi v m (a mi ) Probe-speed of pair m in region a mi S m (a mi ) PTS of pair m in region a mi .S m (a mi ) ≡ (q m (a mi ), k m (a mi )) q mi Probe-flow of datum i of pair m. q mi ≡ q m (a mi )

Figure 1 :
Figure 1: Example of a considered FD with g = 6 and h = 3.
Fig. 3 shows traffic flow consisting of free-flow traffic and congested traffic as a time-space diagram (top) and a flow-density plane (bottom).In the time-space diagram, probe vehicles consisting of pair m are denoted by solid lines, non-probe vehicles are denoted by dashed gray lines, and a shockwave is denoted by an arrow.In the flow-density plane, dots indicate traffic states or PTSs, solid lines indicate transition of traffic states; the dashed line denotes the FD, and the thick dashed line denotes the PFD.Notice that the FD and PFD are similar (in the geometric sense) with scale of 1/c m as in Eq. (9).The blue area (top) and dot (bottom) indicate a steady region with free-flow traffic, the green ones indicate non-steady regions, and the red ones indicate steady region with congested traffic.

Figure 3 :
Figure 3: Illustration of triangular FD identification.Top: Probe pair and PTS in time-space diagram.Bottom: FD, PFD, and PTS on flow-density plane.

Figure 4 :
Figure 4: Triangular FD with its three parameters u, w, and α.Black curve denotes FD, and red curve denotes the corresponding PFD of pair m.

Figure 5 :
Figure 5: Probe vehicle trajectories on a day.

Figure 6 :
Figure 6: Reference flow-density relation based on the detector data.

Figure 7 :
Figure 7: Comparison between speed measured by probe vehicles and detector.
Distance traveled by vehicle n in time-space region A t n (A) Time spent by vehicle n in region A S Traffic state.S ≡ (q, k) Q Density-to-flow FD function Q j Function of j-th piecewise differentiable part of an FD g Number of parameters of an FD g j Number of parameters of j-th piecewise differentiable part of an FD h Number of piecewise differentiable parts of an FD κ Jam density u Free-flow speed of a triangular FD w Backward wave speed of a triangular FD α y-intercept of the congestion part of flow-density triangular FD. α ≡ wκ m Probe pair (i.e., pair of two probe vehicles) m − Preceding probe vehicle in pair m m + Succeeding probe vehicle in pair m c m Number of vehicles between probe vehicles in pair m plus one φ A parameter (angle) used to define a shape of region a ∆t A parameter (width) used to define a shape of region a α m y-intercept of the congestion part of flow-density triangular PFD of pair m. α m ≡ wκ/c m θ steady

•
Jam density: κ • Lower and upper bounds for free-flow speed: u min and u max • Threshold for near-steady state θ steady • Cleaning parameter for PTSs based on Corollary 3: θ corr , θ std .• Parameter sets determining shape of time-space region a: ∆m, ∆t, and φ • Initial state of the EM algorithm: p mi (s) • Threshold for the convergence of the EM algorithm: ε F • Standard deviation of congested flow from mean FD: σ C • Number of vehicles between the probe pairs: c m

Table 1 :
Summary of estimation results.
mi Probe-density of datum i of pair m. k mi ≡ k m (a mi )S miPTS of datum i of pair m. S mi ≡ (q mi , k mi ) = S m (a mi ) s Contribution of state s in the mixed distribution z mis 1 if S mi belongs to regime s, 0 otherwise p mi (s) Expectation of z mis