Estimation of metro network passenger delay from individual trajectories

Smart card data enables the estimation of passenger delays throughout the public transit network. However, this delay is measured per passenger trajectory and not per network component. The implication is that it is currently not possible to identify the contribution of individual system components – stations and track segments – to overall passenger delay and thus prioritize investments and disruption management measures accordingly. To this end, we propose a novel method for attributing passenger delays to individual transit network elements from individual passenger trajectories. We decompose the delay along a passenger trajectory into its corresponding track segment delay, initial waiting time and transfer delay. Using these delay components, we construct a solvable system of equations, using which the delays on each network component can be computed. The estimation method is demonstrated on one year of data from the Washington DC metro network. Our approach produces promising results by compressing millions of individual trajectories into 3D networks, leading to a dimensionality reduction of 94%. Moreover, the mean slack variable value (that can be interpreted as proxies for estimation errors) is smaller than five seconds per passenger, and has the desired positive sign for almost 90% of all travelers. Applications using the estimation results include revealing network-wide recurrent delay patterns, modeling delay propagation and detecting disruptions.


Introduction
Service reliability is known to be one of the most important determinants of transit performance, ridership and user satisfaction. Traditionally, service reliability has been focused on vehicles rather than customers using measures such as headway adherence or schedule punctuality (Trépanier et al., 2009). However, there is a strong shift towards using passenger-oriented measures for quantifying the transit performance (Zhao et al., 2013;Hendren et al., 2015). The data with which many of these quantities can be directly or indirectly estimated are already collected by many transport authorities through different sensors and information sources. Some of the well-known and increasingly used data sources for public transport network are infrastructure (stations and track segments) and service network (line) information, timetable data, automatic vehicle location (AVL) data which contains real-time locations of transit vehicles (Moreira-Matias et al., 2015) and by implication the realization of the schedule and automatic fare collection (AFC) or smart card data with origin-destination specific information of passengers (Pelletier et al., 2011).
An increasing body of literature is available on the use of smart card data in public transit; reviews of which can be found in e.g. Pelletier et al. (2011) and Koutsopoulos et al. (2017). Most of these studies focus on origin-destination matrix estimation (Nassir et al., 2011;Gordon et al., 2013), extracting passenger patterns (Ma et al., 2013;Bhaskar and Chung, 2014;, network performance analysis (Ma et al., 2015; and passenger route choice determinants (Yap et al., 2018). AFC systems can be classified based on: (i) whether they include tap-in only or tap-in and tap-out records; (ii) whether fare validation is performed upon boarding (and possibly alighting) the transit vehicle or upon entering (and possibly leaving) a transit station. Depending on the two former aspects, destination (for boarding only) and vehicle (for station validation) inferences might be applied. Methods have been developed to infer the alighting station of a given tap-in record (Trépanier et al., 2007;Sánchez-Martínez, 2017). In addition, new methods have been proposed for the train inference of each passenger, referred as Passenger-to-Train-Assignment (Zhao et al., 2016;Zhu et al., 2017).
By comparing the timetable and AVL data, the vehicle delay of the individual transit vehicles can be directly determined. However, in a shift towards passenger-oriented measures, it is more interesting and relevant to investigate the passenger delay for a given line, network segment or transfer station as it encompasses the delay incurred by a passenger within the entire public transit system and not just the on-board delay. Clearly, with the increasing availability of AFC data and progress in Passenger-to-Train-Assignment research, it is possible to estimate and quantify the delay experienced by each passenger. However, it remains unknown how much of the passenger delay can be attributed to individual network elements. In particular, in the context of transit, delays can be associated with different travel time segments -initial waiting time, on-board time and transfer times.
To this end, we study passenger delays on individual network elements by fusing AVL and AFC data. The key contribution of this paper thus is a new data-driven method to derive PT network delays from individual trajectories. The methodology has similarities with the data-driven OD estimation method we propose in Krishnakumari et al., in that we construct a solvable system of equations utilizing all the information at hand without making more assumptions than strictly needed. The usage of Passenger-to-Train-Assignment data makes this study unique because this type of information is relatively new. To the authors' knowledge, this is the first study to explore the potential application of this unique data set. In this work, we distinguish two different types of passenger delay in relation to the public transit network: average passenger delay and total passenger delay. We define the average passenger delay as the delay incurred by a passenger while traversing a track segment (link), station (node) or trajectory. The total passenger delay is the total delay experienced by all the passengers that traverse that link, node or trajectory. Thus, the total passenger delay is a function of the number of passengers that traverse that network element during a given time period.
In the remainder of this paper, we show that realized passenger-trajectories (resulting from AFC data and Passenger-to-Train-Assignment) and schedule information are sufficient for estimating average and total passenger delays for all the different network elements. The estimation improves by adding a constraint derived from AVL data, because these are-in our case-readily available. With these constraints for each passenger and each vehicle, a solvable system of equations can be formulated. There are various applications for the resulting network indicators, such as identifying key bottlenecks and critical network elements, prioritizing investments and maintenance of assets such as switches, modeling delay propagation through the network and an automated disruption detection. We demonstrate the estimation framework for the metro network of Washington DC using one year of data.
The paper is organized as follows: Section 2 describes the overall estimation framework; in Section 3 we apply the framework on the Washington DC metro network. We outline the network and data used in this section and present the estimation results. We offer conclusions and a discussion on further research avenues in Section 4.

Methodology
For convenient reference, the notation used for recurrent variables in the methodology is first presented as follows:

Problem formulation
where N is the total number of passengers, is defined based on two sets, one of stops and one of lines, each of which are a subset of S and L, respectively. The combination of which allows one to define a third set: the set of track segments (subset of E) that is traversed by the passenger. As a result, we obtain the initial stop (first element in the stop set), the intermediate transfer stops (second to one before last element in the stop set) and the link set. The stop set is defined as In the following, we assume that the origin and destination station, and the trajectory r s s n , , o d of each passenger are known. The passenger trajectory contains all the trip-legs for a given passenger journey with the corresponding time stamps such as the tap-in, tap-out and transfer information. Depending on the fare validation scheme, these variables can either be directly observed or need to be inferred. If the destination of the passenger is not known, we can apply the by now well-established destination inference methods (Trépanier et al., 2007;Sánchez-Martínez, 2017). If the trajectory of the passenger is not known, one can employ trajectory inference methods such as Passenger-to-Train-Assignment  and ODX method (Sánchez-Martínez, 2017). In case of vehiclebased tap-validation, there is need to also determine whether two subsequent transactions constitute parts of a single journey by performing a transfer inference (Yap et al., 2017). Depending on the fare validation scheme, inferences (destination, trajectory and/ or transfer) are performed so that the trajectory r s s n , , o d of each passenger is given as input to the method proposed in this study, so that the passenger delay d s s k n , , , o d is defined as: P. Krishnakumari, et al. Transportation Research Part C 117 (2020)  , o d can be defined in several ways. In this study, the scheduled travel time include walking times (at both ends as well as at the transfer location if any), maximum waiting times assumed to equal the service headway (of the first leg and any other leg in case of transfers) and the scheduled on-board travel time between boarding and alighting stations (for each leg). These values were made available by the local public transport authority. We then formally define the scheduled travel time as follows: at the origin station and transfer stations denoted by s i . The headway between transit services is included in calculating the schedule travel time due to the definition of on-time journey considered in this study. Note that a passenger is considered to be on-time with regards to the maximum scheduled journey time even he/she just missed a transit vehicle but was able to catch the next service of that line and the headway between the services account for that fallback time.
Given the individual trajectories of the passengers, the aim is to decompose the delay experienced on a given passenger journey into the corresponding network elements they traversed along their trajectory. The estimated delays for each network element are assumed common (or average) for all individuals. We assume that the   P. Krishnakumari, et al. Transportation Research Part C 117 (2020) Based on the type of the fare validation scheme, the passenger delay definition may differ. For example, the initial waiting time delay for a surface bus system cannot be observed from the AFC data and hence, the waiting delay component will then be removed from (4). The delay definition in (4) can be further simplified or expanded based on the specifics of the transit system subject to investigation.
The inequality in (4) is due to potential non-observable personal travel components which may add additional delays such as performing an activity at one of the stations within the gated area. While on-board, waiting and transfer delays relate to those delays that can be attributed to service elements, any other delay -including possibly walking time related -is assumed unsystematic and is therefore attributed to the individual. Formulating the relationship between these personal delays and the three network related delay components by means of an inequality constraint allows us to perform the delay inference without ignoring unobserved personal delay components. We define this unobserved error term as a slack variable for each passenger; thus (4) can be reformulated as: An additional constraint can also be formulated, that pertains to the on-board delay component, which must be equal to the transit vehicle delay between the corresponding trip segments. This vehicle delay can be directly inferred from the AVL data and the schedule information as follows: This results in an additional set of equations to solve for the unknowns in (5) which reads where

Formulating a solvable system of equations
We decompose each of the passenger delay into delay at three network elements -origin stations, track segments and transfer stations and each transit vehicle delay into delay per track segments. With (5), now each passenger trajectory can be written as a linear combination of these three passenger delay component types with an additional slack variable for the error. Our main hypothesis is that formulating (5) for all passenger trajectories leads to a potentially solvable system of equations, since each passenger's trip between k and + k t s s k n , , , o d serves as a constraint for all other trips that traverse one or more common network elements during this trip.
We demonstrate this point with an example. Consider the toy network shown in Fig. 1. If a passenger traveling from s 1 to s 5 experiences no delay while a passenger traveling from s 1 to s 6 experiences a delay of say 5 min, then this delay is probably due to delay occurring between stations s 4 and s 6 . Thus, given a sufficient number of passengers, there exists a bounded solution for (5) Krishnakumari, et al. Transportation Research Part C 117 (2020) Note that the stations and transfer stations are directed for public transport network. The number of station directions depends on the number of its outgoing neighbors. Thus, for station node s 4 in Fig. 1, there are three directions possible, either towards s s , 3 5 or s 6 . This implies that the number of directed stations in a public transport network is the same as the number of directed edges. We thus refine (8) and (12) based on these directed travel nodes as follows: Next, we can reformulate the passenger delay given in (5) Similarly, we can formulate the transit vehicle delay given in (7) as: Based on the example network from Fig. 1 We can generalize (18) and (19) to N passengers and V vehicles respectively to build a system of equations and thus a matrix equation as: where P. Krishnakumari, et al. Transportation Research Part C 117 (2020) The x is shorthand for the delay attributed to each network element; C contains ones or zeros (LHS of Eqs. (18) and (19)); is the slack variable for each passenger and B corresponds to the passenger and vehicle delay (RHS of Eqs. (18) and (19)), respectively. A zero in the C matrix implies that the corresponding network element does not contribute to the respective passenger or vehicle delay in B and a value of one imply that the corresponding network element is part of that passenger's or vehicle's journey and hence contributes to the delay. Thus, Eqs. (18) and (19) } respectively. The matrix equality given in Eq. (23) can be solved using a constrained linear least squares solution (Altman and Gondzio, 1999) with the lower bound set to 0 to ensure a non-negative solution as follows: In case a non-negative solution does not exist, an ordinary least square solution is computed and the negative values of the delay matrix are ignored when computing the estimation error. The minimum requirement to solve the system of equations for a given time period is that at least one passenger traverses each of the network elements.

Evaluation metrics
An optimal solution is reached when is minimized. Thus, we use the slack variable in Eq. (27) to evaluate the estimation results. The delay estimates x of the individual network elements are used to reconstruct the individual passenger delay by multiplying C and x. This can be used to find the estimation error or the slack variable as follows: This estimation error evaluates the discrepancy between the reconstructed individual total passenger delay and the ground-truth passenger delay at the journey level. Based on our problem formulation, a positive slack variable is favored as this implies that our estimates are less than or equal to the observed delay of the passenger. This is reasonable as some of the delay may be attributed to individual-related factors rather than a delay associated with the respective network element. Conversely, a negative slack variable value implies that our estimation attributes to network elements a total delay that exceeds the delay experienced by a passenger. Consequently, the expected slack distribution is skewed towards positive values whereas the error distribution is assumed to be Gaussian in the least square approach (Hayashi, 2000). However, we choose to employ the least square approach as a first estimation and evaluate if it can reasonably solve our system of equations for our formulation.

Application
In this section, we demonstrate our estimation method on a real-world application. For this, we first explain the data and network used (Section 3.1). Thereafter, we provide some useful descriptive statistics of the data (Section 3.2) and finally, we present the results of the estimation (Section 3.3) and analyze its performance (Section 3.4).

Data
We demonstrate the estimation approach on a data set of smart card data from the Washington DC Area Metro Transit Authority P. Krishnakumari, et al. Transportation Research Part C 117 (2020) 102704 (WMATA) in the United States. The data is composed of one year of smart card data from 19 August 2017 to 28 August 2018 for the entire metro network of Washington DC which contains the Passenger-to-Train-Assignment outputs derived from an application of the so-called ODX method described in Sánchez-Martínez (2017). In general, the ODX method infers the origin, destination and transfer of a given passenger. For Washington metro, the tap-in and tap-out location and time stamp of a passenger are directly available. Therefore, the method only infers the paths and transfers of the (possibly multi-leg) journeys by minimizing a generalized disutility function. The method searches for the path that minimizes the time between the tap-out and the latest possible train arrival time. This is under the assumption that passengers do not engage in non-trip related activities at the destination. Thus, for each passenger, the method provides the following trip legs which includes both the location and the corresponding timestamps access (directly observed), egress (directly observed), ungated transfer and train ride itself. In addition, the dataset also includes the rail movement data, schedule information and disruption log file. The metro network is comprised of 6 lines, 91 stations, 186 links and 9 transfer stations as shown in Fig. 2.

Descriptive statistics
The network has an average ridership of 438 000 rides per day and a total of 157 million rides during the entire study period with an average journey time of 28 min per passenger. Of these, 14% of the passengers experience a delay. For those experiencing a delay, the mean delay is 6 min. Moreover, 39% of passenger trips include transfers with an average of 1.14 trip legs per passenger journey. A breakdown of the number of passengers for different time periods is shown in Fig. 3. There are two distinct peaks in the passenger distribution -morning and afternoon peak.

Estimation results
In this section, we present the results of the delay estimation for the Washington metro network. We chose a temporal aggregation of 30 min for k in (18), since the maximum headway between transit vehicles is 20 min and choosing a temporal aggregation smaller than that would imply that there would be no vehicles between some OD pairs, hence no passengers and consequently no system of  P. Krishnakumari, et al. Transportation Research Part C 117 (2020) 102704 equations. Having an aggregation of 30 min ensures that at least one transit vehicle per time period is included in the system of equations as represented in (18). Fig. 4 shows the estimation results of the average passenger delay for the three network elements for a selected weekday (Thursday), i.e. March 1, 2018, with a temporal aggregation of 30 min. There is no significant track segment delay for this particular day. The waiting time and transfer delay are mapped on a link rather than a node as we incorporate directionality of the node to distinguish journeys between different lines. This allows us to build these three compact 3D graphs for visualizing the delay propagation of each day for each of the delay components. This can be used for evaluating the performance of the metro network on a given day or over a long period of time or estimate the passenger delay incurred between any given origin-destination (OD) pair. Moreover, there were 596 000 rides recorded on this particular weekday and we were able to represent the dynamics of the network using these three 3D networks with the dimensions 3 × 48 × 186 (3 × temporal aggregation × number of links), thus leading to a dimensionality reduction of about 95% (26784/596000).
A 3D graph of a weekend day (Saturday), i.e. March 3, 2018, in the same month is also shown in Fig. 4. This exhibits a significant delay in the red line compared to a normal weekday. From the disruption log file of WMATA, we are able to attribute this delay to late track clearing in the morning which had a cascading effect on the rest of the line. One of the main findings from our delay estimation is that not all the delays have a corresponding explanation in the disruption log file of WMATA, which is maintained manually, and vice versa. Thus, our estimation can be used to enrich the log file with additional incidents that caused passenger delay as well as quantify the consequences of these incidents in terms of the number of passengers affected, the average passenger delay and the spatial extent of the impact of the incident.
We applied the estimation framework on the 359 days of smart card data with around 157 million rides and compressed each day into the three 3D delay networks, leading to an overall dimensionality reduction of 94%. We also used the estimates to explore the distribution of the average and total passenger delays for the entire analysis period decomposed into different network elements and different time periods as shown in Fig. 5. The average passenger delay, shown in Fig. 5(a), is stable compared to the morning and evening peak in the total passenger delay distribution in Fig. 5(b). Thus, delays occur at all time periods, whereas more passengers are affected during the peak periods as can be expected. Figs. 5(c) and (d) show the contribution of each network element delay to the overall delay. In the case of average passenger delay, 59% of the delay is associated with the initial waiting time. However, when the number of passengers affected is considered, track segment delay contributes the most with 41% of the delay.

Validation
We evaluate the validity of the estimation by reconstructing the individual passenger delay and calculating the slack variable. A slack value of 0 implies that the delay experienced by a passenger based on the AFC data is the same as obtained by summing over our estimates for the respective journey travel components. We do except a slight variation in the slack value of each passenger due to heterogeneity in user behavior. This is confirmed by inspecting the distribution of the slack variable using boxplot for the above mentioned selected weekday and weekend day shown in Figs. 6(a) and (b).
In Figs. 6 compact boxplot with small variation is desired as this implies that the delay estimated for all passengers at the time . This is expected as more passengers mean more equations and thus more confidence in the estimates and thus a lesser dispersion. The effect of the number of passengers is also evident from the distribution of the slack for the weekday. The weekday has a smaller number of passengers relative to the weekday and thus the boxplot is also relatively compact. Moreover, the distinctively different passenger distribution than the weekday connect to the slack variable result. Another important finding from Fig. 6(a) and (b) is that some passengers have a negative slack. This is presumably primarily due to the time aggregation. Our estimation is aggregated for a certain k time interval. Assume that a disruption occurred at + k i time whereas a passenger departed at time < + k i, thus avoiding the disruption. However, due to the level of aggregation, that trip will be  P. Krishnakumari, et al. Transportation Research Part C 117 (2020) 102704 categorized as being affected by the disruption, since most of the passengers in that time period experienced that delay. The slack distribution for all the days across different time periods is shown in Fig. 7. As can be seen from the figure, most of the values are distributed around zero and our estimation framework achieves a mean slack of 0.072 min across all the passengers which is promising. Additionally, only 10.76% of all passengers had a negative slack value, with a mean of -5 min. Moreover, when complemented with the passenger distribution shown in Fig. 3, it is evident that the boxplot at a time period with a large number of passengers are more compact with less dispersion compared to the time period with less number of passengers. This is, also, in line with the finding from the individual weekday and weekend.
In order to assess the validity of using the least square method to solve the system of equations, we conducted model diagnostics for the error terms, i.e.the slack variables. There are three key assumptions about the error term that are necessary for the least square approach to provide unbiased, efficient and linear estimators: a) errors follow identical distributions, b) errors are independent, c) errors are normally distributed. We obtain a p-value greater than 0.05 for the Breush-Pagan and Durbin-Watson tests, thus confirming the null hypotheses of homoscedasticity and independence, respectively. Moreover, the Shapiro-Wilk normality test value confirms what we observed in the slack variable distribution that the errors are not normally distributed. It is unclear whether non-normalcy would substantially affect our results although there is evidence to suggest that when the normality assumption is not valid but the other assumptions are, the estimates are still consistent (Box, 1976). Different approaches can be taken to remedy this such as reformulating the problem, warping or employing distribution independent methods.
Since there is no ground truth data for the passenger delay associated with each network element and the estimation framework are based on inferred variables, the estimates might be biased. The biases can stem from either the inferred trajectory resulting from the Passenger-to-Train assignment or from the least square approach itself. The in-vehicle time and transfer time would not be affected by these biases as the system of equations are also constrained by the AVL data (essentially serve as the ground truth data for these estimates). As for the bias from least squares, it can mainly occur when the solution is not unique. This can happen when data is insufficient, but this bias can be avoided (e.g. time slices definition) and relates to the robustness of the approach.

Conclusion
In this study, we propose a new estimation method to map passenger delay into network elements. The outputs of our method can aid in measuring network performance for any given origin-destination pair of the public transportation network and in prioritizing measures for improving service robustness. We decompose the delay along a passenger trajectory into its corresponding track segment delay, initial waiting time and transfer delay. We demonstrate the method using one-year data from the Washington metro network.
Our method estimates how passenger delays are distributed across network elements and can thus assess the contribution of each network element to system performance. We were able to achieve a dimensionality reduction of 94% by representing these individual trajectories as 3D networks. The estimation results show that the confidence of the estimation, measured based on the compactness of the boxplot, increases with an increase in the number of passengers. Overall, our estimation framework achieved a mean slack of 0.072 min or less than 5 s per passenger, which is very promising.
These network element estimates are used to analyze the average and total passenger delay of the metro network. The average passenger delay is more or less stable throughout the day whereas total passenger delay, which accounts for the number of passengers affected, contains two distinct peaks. The two peaks correspond to the morning peak and evening peak with the evening peak associated with greater passenger delay than the morning peak. The initial waiting time contributes the most to the average passenger delay with 59% whereas the track segment delay contributes the most to the total passenger delay with 41% of the total delay being attributed to it.
The estimates generated by the method proposed in this paper open new avenues for future research. In particular, a more sophisticated distribution independent method can be used to solve the system of equations instead of the least square approach. Furthermore, better error assessment for each individual network element estimates can be done using additional personal data such as disaggregated personal travel diaries and video data. The estimation framework can also be adapted to allow non-uniform temporal aggregation to reduce the share of passenger delay records for which the slack variable value is negative. The delay estimates at the individual network element level can be used to reveal hourly, daily, weekly or seasonal delay patterns across the metro network. Moreover, the estimation approach is easily transferable to other locations. For this, it is necessary to ensure that the system of equations is solvable and that a different temporal aggregation might be needed depending on the headway between the transit vehicles. Furthermore, it can help in understanding the delay propagation through a network and potentially contribute to the prediction of such delays and their spill-over impacts. Such advancements can ultimately help in supporting decision making in relation to improved service robustness at both the tactical level (e.g. locating switches to allow for short-turning) and operational level (e.g. disruption mitigation strategies including information provision and resource allocation).