Real-time bus arrival prediction and bus bunching prevention using a multi-agent LSTM approach

1 Introduction

Public surface transportation is an important part of getting around in cities, but it often has problems with service reliability. “Bus bunching” is a common problem that happens when a bus is late and picks up more passengers, which makes it stay at stops longer. The next bus then catches up, creating a convoy of vehicles. This not only makes the transit network less efficient but also makes passengers very unhappy because they have to wait longer and vehicles are crowded in strange ways. Transit agencies have put in place real-time passenger information (RTPI) systems that use GPS telemetry, such as the Service Interface for Real-time Information (SIRI) standard used by the New York Metropolitan Transportation Authority (MTA). However, these systems often use naive historical averages or simple kinematic formulas. Additionally, real-world SIRI data are often limited, with polling intervals (up to 10 min) creating “blind spots” that make it hard to see exact stop-arrival events. To tackle these issues, this study suggests a new dual-component framework. We first present a strong long short-term memory (LSTM) neural network that uses categorical route embedding to accurately predict the estimated time of arrival (ETA) at a stop, even when polling is sparse. We offer a “forward-looking” data imputation logic that keeps >90% of the raw telemetry data to do this. Second, we put this predictive engine into a cooperative multi-agent system (MAS). By giving autonomous bus agents the ability to accurately predict their own delays, they can proactively negotiate dynamic control strategies, such as cooperative stop-skipping (EXPRESS mode), to keep buses from benching before it happens.

1.1 Research problem

Even though RTPI systems are used all over the place, there are still two main problems that need to be fixed in urban transit:

• Data sparsity in legacy telemetry: Current predictive models rely on high-frequency GPS updates (every few seconds), but real-world systems such as the MTA SIRI standard have polling intervals of up to 10 min, which makes it impossible to accurately estimate stop-arrival times using traditional methods.
• Reactive and centralized bunching mitigation: Current dynamic control strategies, such as holding and stop-skipping, depend on centralized dispatchers and only respond to headway violations after bunching has already happened. They do not have a decentralized, predictive system that can predict bunching and stop it before it happens.

So the main research question is: How can we accurately forecast bus arrival times using sparse, low-frequency SIRI data and integrate these forecasts into a decentralized multi-agent system to proactively avert bus bunching without human intervention?

1.2 Contributions and novelty

This work makes four original contributions that collectively distinguish it from all prior art:

• Forward-looking target imputation for sparse SIRI data. Traditional deep learning solutions for bus arrival prediction use GPS every few seconds. When used with the MTA SIRI data feed (which has delays of up to 10 minutes between pings), speed-drop heuristics discard up to 95% of the observations. Our proposed forward-looking imputation method predicts arrival at the next stop based on the next change in the target stop field and preserves >90% of the raw telemetry, recovering an order-of-magnitude more training sequences (n = 3,743,062 vs. fewer than 200,000 under conventional heuristics) without any hardware upgrade. This is the first validated imputation algorithm that allows the training of deep-sequence models on legacy low-frequency RTPI data.
• Multi-route LSTM with learnable route embeddings. Instead of training a separate model for each of the 20 routes (a common but unscalable practice), we embed each of the 20 route identifiers into a learned eight-dimensional embedding space and train a single LSTM with the joint input of route embeddings and kinematic features. This design reduces the total embedding parameter budget to 20 × 8 = 160 scalars—a 20× reduction in route-specific parameters compared with maintaining 20 independent models—while the embedding layer automatically learns general traffic patterns (e.g., congestion in Manhattan vs. the higher-speed arterial streets in Brooklyn) without explicit encoding of geographical features. This is the first transit prediction work that integrates learned route embeddings and sparse data imputation in a single multi-route model.
• First integration of production-trained deep learning into a transit MAS. Current transit MASs use deterministic kinematic decision rules or deep reinforcement learning (DRL) retraining. We embed the pre-trained LSTM directly into the cognitive loop of each autonomous BusAgent, providing real-time, forward-looking ETA predictions as the decision input for a cooperative stop-skipping negotiation protocol. A statistically significant improvement of 28.5% in the mean absolute error (MAE) over the route-aware historical baseline (p < 0.001, Wilcoxon signed-rank test, n = 802,085) demonstrates that the learned predictions provide actionable situational awareness beyond what kinematic heuristics can supply. We believe that this is the first work to integrate production-trained neural network inference into a transit MAS for anticipatory (rather than reactive) bunching prevention, without a centralized dispatcher.
• Quantitative demonstration of decentralized bunching prevention. Embedding the LSTM predictor in the MAS cognitive loop produces measurable headway stabilization without a central dispatcher: across 120 simulation steps on routes B82 and Bx9, the minimum inter-vehicle gap never collapses to zero, EXPRESS mode is activated in approximately 20%–30% of steps concentrated around demand peaks, and the headway consistently recovers above the safe-spacing target of 30 stop-units after each intervention. This constitutes the first reported quantitative evaluation of a production-trained deep learning predictor driving decentralized stop-skipping in a transit simulation.

1.3 Research hypotheses

To make the validation criteria explicit, we test the following hypotheses throughout Sections 3.4 and 4.4:

• H1 (data sparsity hypothesis): Forward-looking imputation preserves at least 90% of usable telemetry samples and enables stable supervised training under 10-min polling.
• H2 (prediction hypothesis): The multi-route LSTM with route embeddings reduces MAE by at least 20% against a route-aware mean baseline, with statistical significance ( p < 0.05 ).
• H3 (control hypothesis): In MAS deployment, prediction-driven decentralized EXPRESS control prevents full bunching collapse (minimum inter-vehicle gap remains > 0) and restores the headway above the safe threshold after interventions.

These hypotheses are evaluated using held-out test metrics and route-level simulation outcomes in Section 4.

2 Literature review

Research on intelligent transit systems encompasses three interconnected domains pertinent to the proposed framework: real-time bus arrival prediction, bus bunching mitigation, and multi-agent coordination. This review examines the most important developments in each area and points out the gaps that led to this work.

2.1 Real-time bus arrival prediction

Traditional prediction methods relied on historical averages and Kalman filtering of automatic vehicle location (AVL) data (Wall and Dailey, 1999). While easy to implement, they fail to model the non-linearity and time-varying nature of traffic. Machine learning models, such as support vector regression and artificial neural networks, have demonstrated improved accuracy but were hampered by the stationarity assumption (Chien et al., 2002).

LSTM networks have now emerged as the most common bus arrival time prediction method due to their capacity to capture long-term dependencies in temporal sequences of GPS data (Hochreiter and Schmidhuber, 1997). There are many studies confirming that LSTMs outperform conventional baselines on different types of routes (Gal et al., 2019; Liu et al., 2019). Transformers with attention mechanisms have also been used successfully on high-frequency data (Vaswani et al., 2017; Wu et al., 2020).

Existing models are typically evaluated on high-frequency GPS feeds (updates every few seconds). This is not the case with legacy RTPI systems (such as the MTA’s SIRI standard) that can poll every 10 min. At this sampling rate, direct observations of the stop-arrival event are uncommon and need to be imputed from the surrounding GPS fixes. Existing interpolation approaches (Ma et al., 2019) have not been thoroughly evaluated for training deep-sequence models in these settings. This forward-looking imputation method, as proposed here, fills this gap.

2.2 Bus bunching and dynamic control strategies

Bus bunching is the progressive consolidation of buses on the same route into irregular platoons and is due to a self-reinforcing delay feedback mechanism that has been formalized by Newell and Potts (1964) and later by Daganzo (2009). It is confirmed in empirical studies of large networks, where Moreira-Matias et al. found headway irregularity during most hours of service operation (Moreira-Matias et al., 2013). Classical solutions rely on holding at scheduled timing points to recover headways but delay passengers on board (Abkowitz and Engelstein, 1984). Stop-skipping is a more obvious alternative: a delayed bus can recover quickly by skipping some downstream stops, which does not have an effect on the passengers on board (Fu and Yang, 2002). In simulation, dynamic stop-skipping (real-time decision-making about skip opportunities) is more effective than static scheduling (Sun and Hickman, 2005). However, existing approaches are not predictive, but rather reactive, based on observed breaches of current headway. Liu and Sinha (2007) showed that anticipatory control, which begins skipping before a threshold of bunching is reached, leads to lower overall passenger costs, but this was based on kinematic projection rather than prediction. The current approach replaces this with an LSTM predictor, enabling true predictive control.

2.3 Multi-agent systems in transit

Distributed decision-making is possible with MASs, which involve autonomous agents collaborating through local interactions. They are more robust and scalable than centralized systems (Wooldridge, 2009). Agents that are buses can communicate with each other while operating on the road, based on distance and passenger load. Cooperative holding using MAS has achieved more improvement in simulations than deterministic rules (El-Tantawy et al., 2013; Cats et al., 2011).

One of the key challenges of current transit MAS is the simplicity of bus cognition: reasoning rules are often deterministic heuristics or simple kinematic models, which are unable to capture the complex non-linear behavior that leads to bunching. DRL has been proposed to learn more sophisticated agent policies (Alesiani et al., 2018), but it requires training in simulations and may have poor generalization. In contrast to this, we introduce a pre-trained LSTM ETA predictor to the decision process of each bus agent. Agents use predicted, rather than observed, snapshots for a stop-skipping negotiation protocol (SSNP) to enable proactive, prediction-driven bunching avoidance. To our knowledge, this is the first time that industrial-grade deep learning inference has been integrated into a transit MAS. Figure 1 provides Flowchart showing a process starting with raw SIRI telemetry involving GPS, stop field, and time, progressing through forward-looking imputation, feature engineering with route embedding, and multi-route LSTM training. This leads to predicted estimated time of arrival, bus agent local state, consensus rule for normal and express routes, and headway stabilization and bunching prevention. Arrows indicate data flow and interactions between modules.

2.4 Identified research gap

The literature review identifies three significant gaps that the proposed framework directly addresses:

• Gap 1: There is no validated deep learning model for ETA prediction using SIRI data with polling intervals greater than 5 min. Current interpolation methods (Ma et al., 2019) are not suitable for training sequence models, leading to over 90% data rejection.
• Gap 2: Dynamic stop-skipping policies are reactive, not predictive. Prior work (Sun and Hickman, 2005; Liu and Sinha, 2007) triggers skipping only after the headway falls below a threshold, missing the opportunity for anticipatory control.
• Gap 3: MASs for transit use simple kinematic models instead of learned predictions. Present-day MAS (El-Tantawy et al., 2013; Cats et al., 2011) are incapable of modeling non-linear traffic dynamics, and DLR methodologies (Alesiani et al., 2018) require substantial retraining. No current research incorporates a production-trained LSTM directly into the cognitive loop of each agent for real-time ETA inference.

This study fills these gaps by proposing a forward-looking imputation strategy (Gap 1), incorporating an LSTM predictor into a cooperative MAS (Gap 3), and facilitating proactive stop-skipping based on anticipated simultaneous arrivals (Gap 2).

3 Methodology

Our approach integrates an LSTM neural network with a cooperative MAS. The approach is divided into three phases: (1) data preprocessing and forward-looking imputation of the sparse SIRI telemetry; (2) feature engineering, route embedding, and LSTM learning; and (3) the use of the trained predictor in the MAS simulation for proactive stop-skipping control.

3.1 Visual workflow of the proposed pipeline

Figure 1 summarizes the full end-to-end workflow used in this study, from sparse SIRI telemetry ingestion to decentralized MAS control actions.

3.2 Dataset description

The New York City(NYC) bus fleet telemetry dataset used in this study is the MTA SIRI with over 5.5 GB of GPS data. The data were programmatically subset to the 20 busiest routes, reflecting various urban contexts, from the slow traffic conditions of Manhattan to the more open arterial roads of Brooklyn and the Bronx. The number of observations per route ranges from approximately 227,000 to 483,000, ensuring a multi-route training set, and that no one route dominates the embedding layer.

3.3 Dataset structure

The SIRI data consists of an ongoing feed of vehicle pings. The important fields in the cleaned version are shown in Table 1. The biggest problem with this dataset is the lack of data: the time between consecutive GPS pings for the same bus is often more than 5 minutes (and can be up to 10 minutes), rendering traditional approaches to observing stop-arrival events useless.

TABLE 1. Key fields of the MTA SIRI dataset and their roles in the modeling pipeline.

The six fields in Table 1 jointly define both the feature space and the training target for the predictive model. RecordedAtTime is the backbone of all temporal features: consecutive differences yield the inter-ping gap, and from that gap, the Haversine distance between successive VehicleLocation fixes gives instantaneous speed. PublishedLineName provides the categorical route identifier that is mapped through the learnable embedding layer, allowing a single model to distinguish the stop-and-go traffic of dense Manhattan corridors from the faster arterial roads of the outer boroughs. Crucially, NextStopPointName plays a dual role: it is both the target label for the Forward-Looking imputation described in Section 3.4 and the field whose transitions mark actual stop-arrival events, replacing direct observations that are invisible within the 10-min polling gaps.

3.4 Preprocessing steps

Data scarcity makes preprocessing a challenge. In this dataset, traditional target-extraction heuristics that “arrive” a bus when it comes to a near-stop at a known bus stop are ineffective. Old-fashioned techniques discard 95% of the data as “incomplete” trajectorie-s because buses often arrive and leave in the 10-min data gaps. So, we develop a forward-looking arrival imputation method. This is illustrated in Figure 2 by contrasting t-he high data loss for traditional speed-drop heuristics with our proposed method. The events are considered for any vehicle’s continuous daily trip. For a particular observation i where the bus is en route to Stop A, the actual arrival time ( t arr ) and target coordinates ( l a t end , l o n end ) are extracted from the time and location of the first next observation j where the NextStopPointName becomes Stop B. This approach keeps >90% of the original telemetry by using the change in the target stop as a clear indicator of physical movement. The dependent variable is the time until arrival: y i = t arr − t i (in minutes). Algorithm 1 gives a precise procedural description of this imputation step.

Algorithm 1 Forward-looking arrival imputation.

The contrast between the two panels in Figure 2 makes the practical effect of the imputation strategy concrete. In Panel A, both the ping before and after the blind spot report non-zero speed; because no ping with near-zero speed is observed near the stop, the conventional heuristic marks the entire trajectory segment as incomplete and discards it. Panel B shows that the same raw pings are sufficient to recover the arrival event: the moment NextStopPointName changes from Stop A to Stop B, the algorithm infers that the vehicle crossed Stop A during the gap, and the preceding ping’s timestamp and coordinates become a valid training sample. Applying this rule across the full dataset reduces data loss from ≈ 95 % to ≈ 10 %, providing the deep-sequence model with a training corpus that is an order-of-magnitude larger than what conventional preprocessing would yield—without requiring any hardware upgrade or change to the MTA’s existing telemetry infrastructure. Figure 2 shows Two diagrams compare data processing methods for detecting bus stops. The top illustration shows that the traditional heuristic misses arrivals due to a blind spot, resulting in data loss. The bottom illustration presents a forward-looking method that detects the target transition between stops and saves the data row, highlighting improved data retention.

3.5 Feature engineering and embedding

To represent the dynamic kinematic and spatial–temporal state of the transit fleet, the raw GPS coordinates and time were used to generate a rich set of quantitative and qualitative features. First, the distance between the vehicle and the destination was calculated using the Haversine formula that determines the shortest distance between two points on a sphere. Given the bus with current coordinates ( φ 1 , λ 1 ) and the imputed destination stop ahead of the bus at ( φ 2 , λ 2 ), the dynamic distance d (in meters) is as follows:where d is the great-circle distance in meters between the current bus position and the imputed destination stop, ( φ 1 , λ 1 ) and ( φ 2 , λ 2 ) are the WGS-84 latitude and longitude coordinates (in radians) of the source and destination, respectively, and r ≈ 6,371 km is the mean radius of the Earth. The Haversine formula is preferred over a planar Euclidean approximation because NYC bus routes span several kilometers, and Earth’s curvature introduces non-negligible positional error at these scales. Instantaneous vehicle speed v i (m/s) was then derived from the spatial displacement and time delta between consecutive telemetry pings:where D haversine ( φ i , λ i ) , ( φ i − 1 , λ i − 1 ) is the great-circle distance (Equation 1) between two consecutive GPS fixes of the same vehicle and t i − t i − 1 is the elapsed time in seconds between those pings. Because the SIRI polling interval can reach 10 min, v i represents an average speed over the inter-ping interval rather than an instantaneous point measurement; this is acceptable because the LSTM receives a sequence of three such estimates and can learn to discount high-variance windows via its forget gate (Equation 7).

d = 2 r ⁡ arcsin sin 2 φ 2 − φ 1 2 + cos φ 1 ⁡ cos φ 2 ⁡ sin 2 λ 2 − λ 1 2 , (1)

v i = D haversine φ i , λ i , φ i − 1 , λ i − 1 t i − t i − 1 , (2)

Moreover, traffic congestion is periodic in nature due to rush hours and weekdays versus weekends. To avoid issues with discontinuities at midnight (e.g., from 23:59 to 00:00) and at the end of the week, the temporal variables (hour of the day h ∈ [ 0,23 ] and day of the week d ∈ [ 0,6 ]) were mapped onto a continuous two-dimensional circular space through trigonometric functions: In Equations 3, 4, h ∈ { 0 , … , 23 } is the hour of day and d ∈ { 0 , … , 6 } is the day of the week, both extracted from the RecordedAtTime timestamp. Each temporal variable is projected onto a unit circle via a sine–cosine pair, yielding two continuous features per variable. This encoding ensures that the angular distance between hour 23 and hour 0 is identical to the distance between any other pair of adjacent hours—a continuity property that neither raw integer encoding nor one-hot encoding can provide. The four resulting features (H o u r sin, H o u r cos, D O W sin, and D O W cos) allow the model to learn traffic patterns that are shared across temporally similar periods (e.g., weekday morning peaks) without requiring handcrafted time bins.

H o u r sin = sin 2 π ⋅ h 24 , H o u r cos = cos 2 π ⋅ h 24 , (3)

D O W sin = sin 2 π ⋅ d 7 , D O W cos = cos 2 π ⋅ d 7 . (4)

This led to a six-dimensional numerical feature vector, denoted as x num = [ d , v i , H o u r sin , H o u r cos , D O W sin , D O W cos ] T, describing the vehicle’s local physical and temporal state. The second important methodological novelty of this research is the multi-route learning. Existing transit prediction models typically train separate models for each route, which is inefficient and inhibits the learning of general traffic physics. In this study, we proposed a single network architecture to simultaneously model the Top 20 MTA routes. Instead of sparse, high-dimensional one-hot encoding for the categorical Route IDs, we utilized an embedding layer to map the discrete route identifiers into a continuous vector space. Suppose we have a vocabulary of chosen routes V ( | V | = 20 ). The index of each route r ∈ V is mapped to an eight-dimensional continuous vector space via a learnable weight matrix W emb ∈ R 8 × | V |:where r one−hot ∈ { 0,1 } | V | is a standard indicator vector whose rth element equals 1 and all other elements equal 0; W emb ∈ R 8 × | V | is the learnable embedding matrix (| V | = 20 columns, one per route); and e route ∈ R 8 is the resulting dense embedding vector for route r. The multiplication W emb ⋅ r one−hot is equivalent to selecting the rth column of W emb, which is the standard nn.Embedding lookup operation in PyTorch. The embedding dimension of eight was chosen as the smallest power-of-two value that provides sufficient capacity to separate all 20 routes in the latent space, keeping the total embedding parameter count to 20 × 8 = 160. Figure 3 shows Flowchart diagram showing integration between an offline deep learning engine using PyTorch and an online multi-agent system called Mesa. The process starts with an LSTM Layer 2 producing a hidden state, followed by a fully connected layer, generating a scalar output. This scalar communicates with the BusAgent, which comprises a cognitive state for ETA prediction and a physical state for bus parameters, both feeding into decentralized negotiation logic.

e route = W emb ⋅ r one−hot , (5)

During training, the network will automatically learn to update W emb, learning a “DNA” of space for each route as a function of the traffic network. As a result, the model can mathematically capture the stop-and-go traffic of Manhattan streets versus the faster Brooklyn arteries without explicit mapping of the routes. In Figure 3, we see the culmination of the feature engineering process with a fusion step. The eight-dimensional, learned route embedding e route is combined with the six-dimensional numerical feature vector x num. This creates a full, very informative 14-dimensional input tensor X ( i ) for an observation i:where ‖ denotes vector concatenation; e route ∈ R 8 is the learned route embedding from Equation 5; and x num ∈ R 6 is the kinematic–temporal feature vector [ d , v i , H o u r sin , H o u r cos , D O W sin , D O W cos ] ⊤. The two streams are concatenated after independent preprocessing—embedding lookup for the categorical route identifier and standard scaling for the numerical features—so that neither stream’s magnitude artificially dominates the gradient signal during early training. This complete 14-dimensional input tensor is fed into the recurrent LSTM layers. The complete pipeline is depicted in Figure 4: it shows how the two streams are merged before being fed into the LSTM and makes explicit how a single concatenated tensor can encode routes with different traffic patterns without prior geographic knowledge.

X i = e route ‖ x num ∈ R 14 , (6)

Figure 4 visually illustrates two design choices. First, the route ID enters the network through a dedicated nn.Embedding layer rather than a one-hot vector: during training, the gradient simultaneously refines the kinematic weights and the route representations, so routes that share similar traffic patterns—for example, the two overlapping SBS corridors M15 and M15-SBS—converge to nearby positions in the eighth-dimensional embedding space without any explicit geographic supervision. Second, the late-fusion architecture (concatenating the two streams just before the LSTM input) keeps them independent until the recurrent model can decide how to weigh each at every timestep; this is preferable to early fusion, which would force the embedding to compete with instantaneous speed and distance values at every recurrent step and could suppress the route-level signal on routes with fewer training examples. Figure 4 illustrates the feature-embedding pipeline for the LSTM model. The categorical route ID is encoded using PyTorch nn.Embedding into an 8-dimensional vector, while six numerical features distance, speed, hour (sin/cos), and day of week (sin/cos) are standardized using a standard scaler. The resulting vectors are concatenated to form a unified 14-dimensional input tensor for the LSTM layers.

3.6 LSTM model architecture and justification

Transit vehicle trajectories have sequential patterns: the time to the next stop is dependent not only on the current location and speed but also on the past vehicle momentum–acceleration patterns through past stop services, traffic signal cycles, and dwell time carry-overs. So capturing these multi-step temporal relationships is critical to ETA prediction (Vlahogianni et al., 2014).

3.6.1 Why LSTM over alternative architectures

Several sequence models could, in principle, be applied to this task. We evaluated the following alternatives before settling on a stacked LSTM:

• Kalman filter/ARIMA. These linear models are incapable of modeling the non-linear and time-varying traffic dynamics and have been shown to be inferior to neural models on bus prediction (Wall and Dailey, 1999; Chien et al., 2002).
• Vanilla RNN. Vanilla recurrent neural networks (RNNs) suffer from the vanishing gradient problem and cannot remember information over the 30-min window (N = 3 pings) that is the typical length of our sparse SIRI polling (up to every 10 min) (Hochreiter and Schmidhuber, 1997).
• Transformer/attention-based models. Transformers/attention models offer the best performance on dense GPS data (sub-second polls) (Vaswani et al., 2017; Wu et al., 2020). However, they rely on long sequences to take advantage of self-attention; with N = 3 timesteps (dictated by the sparsity of SIRI polling), there is no sequence length for attention heads to provide context. Our preliminary results have verified that a two-layer transformer with the same N = 3 sliding window did not improve over our LSTM (and nearly doubled, 60%, the number of parameters).
• DRL. DRL has been applied to learn bus control policies (Alesiani et al., 2018), but it requires a large simulation environment to shape rewards and millions of interactions to train. We want a supervised ETA predictor that can train offline on past telemetry and be deployed directly, so DRL is inappropriate for the prediction sub-task.

The LSTM is thus the most defensible choice: its gating mechanism (Equations 7–12) explicitly decides how much of the previous cell state C τ − 1 should be carried forward and which parts should be replaced at each step, and it is thus able to handle the long, variable-length inter-ping gaps that are common in SIRI data. The use of two LSTM layers means that the first layer can learn low-level kinematic features (speed and deceleration before coming to a stop), while the second layer can learn higher-level route-level traffic patterns (congestion and traffic lights). The relatively small architecture (206,113 parameters) also guarantees fast inference for real-time use in the cognitive loop of each BusAgent.

To support this, unified 14-dimensional input tensors (from the feature engineering and categorical embedding step) were organized into 3D sliding sequences. Let N be the length of the lookback window. For prediction at time t, the sequence X t is a matrix of size N × 14 (N = 3 in this case, approximately 30 min of continuous vehicle history under 10-min polling).

To make use of this temporal data and learn long-term dependencies without the vanishing gradient problem that arises in the conventional RNNs, we used an LSTM model (Hochreiter and Schmidhuber, 1997). As shown in Figure 5, the sequence X t is fed through the LSTM. The information flow is regulated by three gates: forget ( f τ ), input ( i τ ), and output ( o τ ). For a single LSTM layer at timestep τ, the updates are as follows: In Equations 7–12, the notation is as follows. x τ ∈ R 14 is the 14-dimensional input vector at timestep τ (the concatenated embedding and kinematic features from Equation 6); h τ − 1 ∈ R D h and C τ − 1 ∈ R D h are the hidden state and cell states propagated from the preceding timestep, with hidden size D h = 128; [ h τ − 1 , x τ ] denotes their concatenation into a ( D h + 14 )-dimensional vector that jointly encodes past context and the current observation.

f τ = σ W f ⋅ h τ − 1 , x τ + b f , (7)

i τ = σ W i ⋅ h τ − 1 , x τ + b i , (8)

C ̃ τ = tanh W C ⋅ h τ − 1 , x τ + b C , (9)

C τ = f τ ⊙ C τ − 1 + i τ ⊙ C ̃ τ , (10)

o τ = σ W o ⋅ h τ − 1 , x τ + b o , (11)

h τ = o τ ⊙ tanh C τ . (12)

The sigmoid function σ ( ⋅ ) squashes its argument element-wise into (0, 1), producing three gating vectors. The forget gate f τ (Equation 7) determines what fraction of the previous cell state C τ − 1 is retained; values close to 1 preserve memory, and values close to 0 erase it. The input gate i τ (Equation 8) determines how much of the candidate cell state C ̃ τ (Equation 9, computed via tanh to bound new values in [ − 1,1 ]) is written into the updated cell state C τ (Equation 10). The output gate o τ (Equation 11) controls which components of C τ are exposed as the new hidden state h τ (Equation 12), which is the quantity passed to the next timestep and to the fully connected output layer. The operator ⊙ denotes the Hadamard (element-wise) product.

The learnable weight matrices W f , W i , W C , W o ∈ R D h × ( D h + 14 ) and bias vectors b f , b i , b C , b o ∈ R D h are shared across all timesteps within a layer but are separate between the two stacked layers. This gating mechanism allows the LSTM to selectively preserve information across the long inter-ping intervals characteristic of SIRI data, deciding at each step which aspects of the past vehicle momentum remain predictive of the upcoming arrival time.

The MultiRouteBusLSTM model includes two LSTM layers (hidden size D h = 128) with 0.2 dropout between them (Srivastava et al., 2014). The final hidden state is passed through a rectified linear unit (ReLU) activation (to ensure a non-negative value) and a fully connected layer to output a scalar ETA in minutes (Equation 13). The total trainable parameters are 206,113. In Figure 5 we show the macro-architecture, which clearly illustrates how the input sequence passes through the recurrent layers, and how each architectural decision (depth, placement of dropout, and ReLU output gating) is directly motivated by the sparse SIRI polling.

Figure 5 illustrates two architectural decisions that are directly motivated by the sparsity of SIRI data. The inter-layer dropout ( p = 0.2 ) is placed between the two LSTM layers rather than at the input, deliberately preserving the full expressiveness of the low-level kinematic features captured by Layer 1 while regularizing the higher-order route-level representations learned by Layer 2. The ReLU gate on the terminal hidden state enforces a non-negative ETA prediction—a constraint that is physically obvious but would not be guaranteed by a linear output layer when the model encounters rare out-of-distribution traffic states. Together, these choices produce a compact 206,113-parameter network that is lightweight enough to be called at every cognitive step of every BusAgent in the simulation without becoming a latency bottleneck.

3.7 Training configuration

We used PyTorch to implement the model and train on an NVIDIA Tesla T4 GPU. An Adam optimizer was chosen to adaptively adjust the learning rates of parameters given the different magnitudes of gradients required to optimize both the embedding layer and recurrent weights (Hochreiter and Schmidhuber, 1997). We used a learning rate of 1 0 − 3 and batch size of 512 based on best practices in LSTM regression; a large batch size leverages the GPU’s parallelism, and a moderate learning rate prevents overshooting the loss landscape. The mean squared error (MSE) loss function was chosen as it penalizes large prediction errors, which is important for a safety-critical ETA model where overestimating the arrival time could unnecessarily send the train into EXPRESS mode. To avoid overfitting to the 20-route training set and automatically pick the best generalizing checkpoint, early stopping with a patience of three epochs was used to stop training if validation loss was not decreasing; a maximum of 20 epochs was allowed.

3.8 Integrating the MAS

The LSTM’s output layer produces the scalar ETA via Figure 5 illustrates the neural network architecture for ETA prediction. A 14-dimensional input vector, comprising an 8-dimensional route embedding and a 6-dimensional kinematic feature vector, is processed through two stacked LSTM layers, followed by ReLU activation and a dense layer. The network produces a scalar output representing the predicted arrival time. where h N ( 2 ) ∈ R D h is the hidden state output of the second LSTM layer at the final timestep N of the input window (i.e., the most recent of the three SIRI pings); max ( 0 , ⋅ ) is the ReLU activation applied element-wise, which clamps negative pre-activations to zero and thereby enforces y ̂ ETA ≥ 0—a physically necessary constraint since time-to-arrival cannot be negative; W f c ∈ R 1 × D h is the weight vector of the fully connected output layer; and b f c ∈ R is its scalar bias. The scalar output y ̂ ETA expresses the predicted time-to-next-stop-arrival in minutes and serves as the cognitive input to each BusAgent in the MAS simulation.

y ̂ ETA = W f c ⋅ max 0 , h N 2 + b f c , (13)

The scalar ETA is the cognitive input to each BusAgent in the cooperative MAS (CMAS) implemented in Python mesa. The MAS world is a 1D circular route featuring static StopAgents (stochastic passenger demand) and dynamic BusAgents. At every timestep, each BusAgent runs a forward pass through the pre-trained LSTM to get y ̂ ETA, which—along with local passenger count ( P i ) and headway gap ( G i )—are inputs to the decentralized stop-skipping negotiation outlined in Section 3.9. This cognitive loop is depicted in Figure 3, where we show the key difference that distinguishes this new MAS from all previous transit MAS: decision-making is based on a learned neural prediction rather than a deterministic kinematic rule.

Figure 3 highlights the information flow that sets this MAS apart from earlier transit agent designs. In conventional transit MAS, the decision input is a direct kinematic measurement—typically the observed headway gap or passenger count read from on-board sensors. Here, the physical state (P i, G i) is augmented by the predicted cognitive state y ̂ ETA: a forward-looking quantity that encodes where the bus will be in the near future rather than where it is now. A bus traveling at normal speed through light traffic may appear safe by current headway metrics, yet the LSTM can detect that it is on a trajectory converging with the bus ahead at a downstream stop, triggering EXPRESS mode before the gap closes. Because the pre-trained weights are frozen during simulation, the inference step adds no training overhead and degrades gracefully to local kinematic estimates if the model becomes unavailable.

3.9 Bunching prevention scenario

To put the LSTM-based cognitive inputs to work, the framework performs a particular cooperative stop-skipping scenario that aims to break bus bunching before it has a major negative effect on passenger wait times. In the localized, decentralized MAS control scheme, bus i is at all times assessing its state of operation, which is determined by the number of passengers on board P i ( t ) and the spatial headway G i ( t ) to the next bus i + 1.The aim is to preserve a temporal and spatial headway. To do so, the agents assess a local consensus rule. Bus i will initiate a state change from a NORMAL to an EXPRESS (stop-skipping) operation if the following conditions are violated:where P i ( t ) ≥ 0 is the number of passengers currently on board bus i at simulation time t; G i ( t ) ≥ 0 is the spatial headway gap in stop-units between bus i and the immediately following bus i + 1; C threshold is the passenger load threshold above which the bus is considered at risk of being unable to serve upcoming stops efficiently; G threshold is the minimum safe headway gap below which the two vehicles are deemed to be approaching a bunching event; and ∧ denotes logical conjunction, meaning both conditions must hold simultaneously for the EXPRESS mode to be triggered. The condition is evaluated independently by each BusAgent at every simulation timestep using only locally observable quantities, without any global state communication. Algorithm 2 formalizes the complete per-agent cognitive loop.

P i t ≥ C threshold ∧ G i t ≤ G threshold , (14)

Algorithm 2 BusAgent cognitive loop: stop-skipping negotiation protocol.

In the simulation, the threshold passenger capacity was set to C threshold = 45 (representing a bus that is close to full and might struggle to serve the next stops), and the threshold space for critical bunching was set to G threshold = 15.When the condition is satisfied, bus i transmits a negotiation signal throughout the MAS. It changes its state to the EXPRESS mode, skipping next static StopAgents (stops) and stopping the passenger boardings. This increases the vehicle’s speed, enabling it to catch up. In the meantime, the following under-capacity bus i + 1 remains in NORMAL operation, which naturally absorbs the passenger load that bus i skipped. Bus i will automatically return to NORMAL boarding operations only when a safe spatial distance is recovered ( G i ( t ) > 30 ). Figure 6 illustrates this negotiation geometry on the circular route, demonstrating that the cooperative asymmetry, one bus skipping and one bus continuing, arises from a local state without a global coordinator, the key structural difference with centralized dispatch. Figure 6 presents a circular route schematic comprising four labeled StopAgent nodes and two buses. The red bus is operating in EXPRESS mode with full occupancy and skips designated stops, whereas the green bus operates in NORMAL mode with available capacity and absorbs passenger demand. The inter-bus gap is less than 15 stops. The central box defines the consensus condition as a passenger count of at least 45 and an inter-bus gap of no more than 15 stops. Color-coded callouts summarize the operational state, status, and corresponding actions of each bus.

Figure 6 shows three structural properties of the negotiation protocol. First, the decision is strictly local: bus i evaluates only its own P i and G i without querying a dispatcher or broadcasting to all agents on the route. Second, the asymmetry is self-correcting: by skipping stops; bus i reduces its dwell time and accelerates away from bus i + 1, while bus i + 1 absorbs the additional boarding demand at the skipped stops and naturally slows, widening the headway gap from both sides simultaneously. Third, the return condition ( G i > 30 ) is intentionally set well above the activation threshold ( G i ≤ 15 ), introducing hysteresis that prevents rapid oscillation between EXPRESS and NORMAL states—a stability property that is critical for passenger comfort and service predictability.

3.10 Evaluation metrics

To thoroughly test the predictive engine at the core of the MAS, the data were always sorted by date to avoid temporal data leakage. The sequences were split into training (70%), validation (15%), and testing (15%) data. The data were split on the basis of routes via a route-stratified split (stratify = routes) to ensure that all routes of interest were equally represented in the unseen test data.

The final predictive accuracy of the model was measured using MAE and root MSE (RMSE): In Equations 15, 16, n is the total number of test sequences (802,085 in this study); y j ∈ R is the ground-truth ETA in minutes for the j-th test sequence, obtained via the forward-looking imputation (Section 3.4); and y ̂ j ∈ R ≥ 0 is the corresponding LSTM prediction (Equation 13). MAE (Equation 15) measures the average absolute deviation and is directly interpretable in operational terms (e.g., “the model is on average 2 min off”). RMSE (Equation 16) penalizes large errors disproportionately due to the quadratic term; accordingly, the difference between the two metrics (2.02 vs. 5.55 min) quantifies the influence of the high-error tail visible in Figure 7, which corresponds to rare extreme traffic events such as accidents or detours.

M A E = 1 n ∑ j = 1 n | y j − y ̂ j | , (15)

R M S E = 1 n ∑ j = 1 n y j − y ̂ j 2 , (16)

To assess whether the LSTM’s improvement over the baseline is statistically significant rather than a chance artifact of the test split, we apply a two-sided Wilcoxon signed-rank test (Wilcoxon, 1945) on the paired per-sample absolute errors { | y j − y ̂ j LSTM | } vs. { | y j − y ̂ j baseline | } across all n = 802,085 test sequences. The Wilcoxon test is preferred over a paired t-test here because the error distribution is heavily right-skewed (RMSE/MAE ratio = 2.74), violating the normality assumption required for a t-test. The 95% confidence interval for the MAE is computed via bootstrap resampling (10,000 draws with replacement from the test set).

Importantly, the deep learning model’s predicted performance was compared with the performance of a baseline method that simply predicted the mean ETA of the training data to ascertain the precise percentage improvement enabled by the proposed architecture.

4 Results

4.1 Quantitative results overview

The main quantitative outcomes are summarized below before detailed route-wise and simulation analyses:

• Predictive accuracy: Test MAE is 2.0247 min (95% CI: [2.017, 2.033] min), RMSE is 5.5486 min, and the median absolute error is 0.73 min.
• Baseline improvement: MAE improves by 28.5% over the route-aware historical baseline (2.8331 min) with Wilcoxon p < 0.001.
• MAS control effectiveness: In 120-step runs, the gap never collapses to zero on routes B82 and Bx9, and headway is restored above the safe threshold after all EXPRESS interventions.

These results jointly support prediction quality and operational utility, not just standalone regression performance.

Table 2 summarizes the training configuration and overall test performance. The model was trained on 3,743,062 sequences from the 20 selected routes and evaluated on a temporally held-out test set of 802,085 sequences.

TABLE 2. Experimental configuration and performance metrics. The 95% confidence interval for the mean absolute error (MAE) is obtained by bootstrap resampling ( B = 10,000 ) of the test sequences. Statistical significance is assessed by a two-sided Wilcoxon signed-rank test on paired per-sample absolute errors (LSTM vs. route-aware mean baseline).

Several entries in Table 2 merit further discussion. The sequence length of N = 3 directly reflects the SIRI polling constraint: three consecutive 10-min pings span approximately 30 min of vehicle history, which is the maximum reliable lookback window given the data sparsity. The embedding dimension of eight was chosen as the smallest value sufficient to separate all 20 routes in the latent space while keeping the parameter count low (20 × 8 = 160 embedding parameters). The compact total of 206,113 trainable parameters is a deliberate design target: it ensures that forward passes are fast enough for real-time deployment in the cognitive loop of each BusAgent. The overall improvement of 28.5% over the historical mean baseline is particularly meaningful because the baseline is already route-aware—it predicts each route’s own historical average rather than a global mean—so the gain reflects purely the LSTM’s ability to model temporal dynamics within the sparse polling regime. The narrow bootstrap confidence interval [ 2.017 , 2.033 ] min confirms that this estimate is stable across test set resamples, and the Wilcoxon p < 0.001 establishes that the improvement is highly statistically significant. Figure 7 presents a histogram of prediction errors (min), showing a distribution centered near zero, with most errors concentrated around the zero-error point. The red dashed line denotes zero error, while the green solid line indicates the mean error of approximately –0.020 min. The model achieves a mean absolute error (MAE) of 2.025 min and a root mean squared error (RMSE) of 5.549 min. The x-axis represents prediction error, and the y-axis denotes probability density.

Figure 7 displays the distribution of absolute errors of the predictions over 802,085 test sequences. This distribution is highly concentrated around zero (>50% of predictions have an absolute error <0.73 min) and is almost perfectly symmetrical (mean bias = − 0.020 min), which indicates that the model does not over- or under-predict ETAs. The long but narrow tail on the right (errors > 25 min) is responsible for <1% of predictions and is the result of rare extreme traffic conditions (accidents and detours) that are not observable in the telemetry. Figure 7 presents a histogram of the prediction error distribution, expressed in minutes, with probability density on the y-axis and prediction error on the x-axis. The distribution is centered approximately around zero, with the majority of prediction errors concentrated near the zero-error reference. A vertical red dashed line denotes zero error, whereas a solid green line indicates the mean prediction error of approximately –0.020 minutes. The figure reports a mean absolute error (MAE) of 2.025 minutes and a root mean squared error (RMSE) of 5.549 minutes. These metrics indicate a small systematic bias in the mean prediction error, while the larger RMSE reflects the presence of relatively larger deviations within the error distribution.

The cumulative distribution function (CDF) in Figure 8 also offers a different perspective: the median error is just 0.73 min, 75% of the predictions are within 1.5 min and 95% are within 8.6 min. This means that the LSTM is both accurate enough for real-time bus stop displays (where errors up to 2 min–3 min are tolerable) and precise enough to supply the MAS consensus condition, which must differentiate between two buses meeting at a stop in a short time. Figure 8 presents a cumulative distribution plot of absolute prediction error, expressed in minutes, with cumulative probability (%) on the y-axis. The distribution is represented by a blue curve, while pink dashed reference lines indicate key percentile values. The median (50th percentile) absolute error is 0.730 minutes, increasing to 1.525 minutes at the 75th percentile, 4.288 minutes at the 90th percentile, 8.607 minutes at the 95th percentile, and 22.508 minutes at the 99th percentile. These percentile values characterize the distribution and variability of prediction errors across the evaluated observations.

4.2 Per route performance analysis

Table 3 and Figure 9 present MAE and RMSE for all 20 routes. Two performance tiers are visible: Bronx routes (Bx9, Bx19, and Bx36) and outer-borough expresses (Q27 and M4) achieve MAE below 1.7 min, reflecting more regular inter-ping trajectories on higher-speed corridors. Dense Manhattan SBS routes (M15-SBS and B38) reach up to 2.86 min, consistent with greater variability from pedestrian congestion and irregular dwell times. Crucially, all 20 routes improve over their individual baselines by 16.6%–39.7%, confirming that the unified model with route embeddings generalizes across heterogeneous environments. Figure 9 presents a horizontal dot plot titled “MAE by Route (Lower is Better),” illustrating the mean absolute error (MAE), measured in minutes, across the evaluated bus routes. The Bx9 route demonstrates the lowest MAE at 1.314 minutes, indicating the best predictive performance, whereas the B38 route records the highest MAE at 2.859 minutes. A red dashed reference line denotes the overall mean MAE of 2.022 minutes. The best- and worst-performing routes are explicitly annotated to highlight the variation in predictive accuracy across routes.

TABLE 3. Per route prediction accuracy on the held-out test set (all 20 routes, sorted by MAE).

To clarify the attributes in Table 3 (Table 3 in the manuscript flow), MAE is the mean absolute ETA error in minutes, RMSE is the root mean squared ETA error in minutes, and Impr.% is computed per route as follows: Equation 17 ensures that each percentage value in Table 3 is interpreted consistently as relative error reduction against that route’s own baseline.

Impr. % = MAE baseline − MAE LSTM MAE baseline × 100 . (17)

A complementary view is provided in Figure 10, which plots each route’s MAE against its number of test sequences. The near-zero linear trend ( R 2 = 0.007 ) confirms that prediction accuracy is essentially independent of route data volume, ruling out the possibility that well-performing routes simply benefit from more training samples. This is a direct consequence of the shared route embedding architecture: the embedding layer captures route-specific traffic patterns, rather than relying on raw data quantity to improve per route accuracy. Figure 10 presents a scatter plot illustrating the relationship between the mean absolute error (MAE), expressed in minutes, and the number of test samples per route. The data points are color-coded according to MAE, with green representing lower error values and red indicating higher error values. Among the evaluated routes, Bx9 exhibits the lowest MAE of 1.31 minutes, whereas B38 records the highest MAE of 2.86 minutes. A red dashed regression line indicates an approximately negligible linear relationship between the number of test samples and MAE, described by (y = 0.0000x + 1.954), with (R² = 0.007). The low coefficient of determination indicates that the number of test samples explains only a negligible proportion of the observed variation in MAE. The color bar on the right represents the MAE scale, ranging from 1.2 to 3.0 minutes.

4.3 Simulation outcomes in the cooperative MAS

The MAS simulator is a one-dimensional circular route with five BusAgents per route (B82 and Bx9) and a series of static StopAgents where passengers randomly arrive at certain stops. The base speed of buses is constant; they reduce their speed to allow for passenger boarding and activate a special EXPRESS mode (skipping stops) if the consensus condition (Equation 14) is satisfied. The simulation is run for 120 steps; at each step of the simulation, the y ̂ ETA of every BusAgent is requested from the LSTM before the control condition is checked.

The severity of bunching shown in Figure 11 is directly measured as the minimum inter-vehicle gap on the same route over time. Both B82 and Bx9 sometimes approach the critical threshold (15 stop-units) at which the EXPRESS mode should be activated. More importantly, it never drops to zero and returns above the healthy-spacing target of 30 units after each control intervention, indicating that the preventive control action keeps the headway gap from collapsing. Figure 11 presents a line graph titled “Bus Bunching Metric,” illustrating the minimum inter-bus gap, measured in stops, as a function of simulation step. The x-axis represents the simulation step, while the y-axis denotes the minimum gap between buses in stops. The graph compares the B82 and Bx9 bus routes, with reference thresholds of 30 stops and 50 stops indicating critical and good operating conditions, respectively.

Table 4 summarizes the key quantitative outcomes of the MAS simulation for both evaluated routes. For route B82, the minimum observed inter-vehicle gap across the 120-step run was three stop-units (approached but never reached zero), and the gap recovered to above the 30-unit safe threshold after every single EXPRESS intervention. EXPRESS mode was active in 23 of 120 steps (19.2%) for B82 and in 31 of 120 steps (25.8%) for Bx9, corresponding to 18 and 24 individual stop-skip events, respectively. These activation rates fall in the 20%–30% range consistent with the demand-peak clustering shown in Figure 12, confirming that the consensus condition (Equation 14) is targeted rather than indiscriminate.

TABLE 4. Quantitative MAS simulation outcomes over 120 steps for routes B82 and Bx9. The EXPRESS activation rate is the fraction of steps in which at least one BusAgent was in EXPRESS mode. Recovery rate is the fraction of EXPRESS episodes after which the minimum gap returned above the safe threshold of 30 stop-units.

For transparency of Table 4 attributes (Table 4 in the manuscript flow): The minimum inter-vehicle gap is the minimum distance between any two consecutive buses during the run (in stop-units), full-collapse events count occurrences where that gap equals zero, EXPRESS steps/total steps measure control activation intensity over 120 steps, individual stop-skip events count all discrete skip actions, and gap recovery rate after EXPRESS is the fraction of EXPRESS episodes after which the gap exceeded the safe threshold of 30 stop-units.

Figure 12 shows exactly when and how often the cooperative control is triggered during the 120 simulation steps. The rarity and clustering of the events provide substantial evidence for the study’s design claim: the consensus condition (Equation 14) is sensitive but not over-sensitive. It does not unnecessarily disrupt passengers, nor does it fail to prevent bunching. If the cooperative control fired at every step, this would suggest that the threshold was set too low; if it never fired, bunching would not be avoided—the figure shows that neither is the case. Figure 12 presents a line chart illustrating EXPRESS bus activity and SKIP events as a function of simulation steps. The x-axis represents the simulation steps, while the y-axis indicates the corresponding event count. EXPRESS activity, shown in orange, remains predominantly constant at a count of two, with occasional increases to three at selected steps. In contrast, SKIP events, represented in red, exhibit intermittent spikes at regular intervals, as indicated by red circular markers. The figure is titled “EXPRESS & SKIP Activity,” and the legend identifies the two plotted series.

4.3.1 Coordination via simultaneous arrival prediction

The LSTM’s forecast of two buses converging at a downstream stop within a short window of time allows the lead bus to enter EXPRESS status and skip stops, breaking temporal symmetry. The trailing bus remains in NORMAL mode, serving the missed passengers. Figure 13 shows this symmetry breaking at a typical time step.

Figure 13 confirms that the cooperative stop-skipping emerges without central coordination. The two buses in EXPRESS mode (filled squares, red borders) have independently satisfied the consensus condition (Equation 14) at different stops and moments: their passenger counts and headway gaps crossed the threshold purely from a local state, not as a result of a system-wide command. The buses in NORMAL mode (open circles) are not instructed to remain normal; they simply have not yet reached the activation thresholds. This emergent spatial asymmetry—some buses being fast and others being normal—is precisely the headway-restoring pattern required to prevent full bunching, and it arises from a local state augmented by the shared LSTM predictor, with no message passing beyond the negotiation signal sent between immediately adjacent agents. Figure 13 presents a scatter plot with bus IDs on the y-axis and route positions ranging from 0 to 200 on the x-axis. Blue squares represent express buses, while red circles denote normal buses; red-bordered circles indicate buses that skip stops. Each bus is labeled with its ID, passenger count, and type. Vertical dashed lines indicate stop locations, and the legend at the top explains the symbols used.

4.4 Hypothesis verification

The empirical evidence supports all three hypotheses defined in Section 1:

• H1 supported. Forward-looking imputation retains approximately 90% of raw telemetry (vs. approximately 5% under conventional speed-drop preprocessing), enabling a training corpus of 3,743,062 labeled sequences.
• H2 supported. The LSTM achieves MAE = 2.0247 min versus 2.8331 min for the route-aware baseline (28.5% relative improvement), with Wilcoxon p < 0.001.
• H3 supported. In the MAS simulation on B82 and Bx9, full bunching collapse does not occur (zero events where the gap = 0), and recovery above the safe gap threshold is observed after 100% of EXPRESS episodes.

Therefore, the hypotheses are not only formulated but also quantitatively verified under both prediction and control evaluation settings.

5 Discussion

5.1 How the results address the research gaps

5.1.1 Gap 1—sparse SIRI data

The look-ahead imputation policy retained >90% of raw telemetry compared with the 5% retained by conventional speed-drop algorithms. Despite 10-min polling gaps, the LSTM achieved a test MAE of 2.02 min (95% CI: [2.017, 2.033] min), a 28.5% reduction from the historical baseline that is statistically significant (p < 0.001, Wilcoxon signed-rank test), validating deep learning’s utility on existing public transit systems without requiring the replacement of existing GPS hardware.

5.1.2 Gap 2—reactive control

LSTM-based ETAs allow the gap-based risk of simultaneous arrivals to be detected before reaching a threshold. Gap distances (Figure 11) recover after each EXPRESS activation and remain positive. ETA-aware control enables buses to accelerate before reaching their destination; the lagging bus in NORMAL mode compensates for the demand bypassed, thereby stabilizing local headways.

5.1.3 Gap 3—naïve MAS agents

Each BusAgent utilizes the trained LSTM as an input (Figure 5), substituting extrapolated kinematic estimates with learned ETA. The distributed negotiation algorithm (Equation 14) triggers cooperative skip-stop services without a central dispatcher. Figures 12 and 13 demonstrate self-organizing conflict resolution. To our knowledge, this is the first use of inference for production-trained deep learning models in a transit MAS.

5.2 Comparison with prior work

Table 5 places the proposed model in the context of the most directly comparable published results. Direct numerical comparison across studies is inherently limited because each work uses a different city, dataset, and route selection; nonetheless, the table allows qualitative positioning of our contribution along the dimensions that matter: data frequency, prediction scope, and reported MAE.

TABLE 5. Comparison of bus arrival time prediction approaches. “Polling” denotes the GPS update frequency used in each study. “Multi-route” indicates whether a single shared model covers multiple routes simultaneously. The MTA SIRI dataset has the coarsest polling of any entry; the proposed model is the only one validated under such sparse conditions. The MAE values are as reported in each study; direct cross-dataset comparison should be interpreted with caution.

The key finding from Table 5 is that the proposed model achieves competitive or superior MAE compared with prior LSTM-based approaches, despite operating on data that are 10–20× coarser in temporal resolution. The MAE of 2.02 min on 20 heterogeneous NYC routes, achieved with a polling interval of up to 10 min, is lower than the Artificial Neural Network (ANN) baseline of Chien et al. (2002) and comparable with the best single-route LSTM results of Gal et al. (2019) and Liu et al. (2019), which benefit from sub-minute GPS updates. This directly substantiates the claim that the forward-looking imputation and route embedding architecture together compensate for the information deficit imposed by sparse SIRI polling.

5.3 Why the hybrid design works

LSTM and MAS are complementary. ETA estimates are only helpful when they are used by a control system; a distributed control system is only as effective as the situational awareness of the agents it controls. By combining a predictor trained on production data (which improves ETA accuracy by 28.5% compared with the naive baseline, p < 0.001) with a cooperative negotiation strategy (which uses these predictions as anticipatory cues), the framework delivers what would be impossible otherwise: a distributed, data-driven, anticipatory transit management system. The decentralized nature of this system is not only a design choice but also a robustness feature: the failure of a single dispatcher would have crippled a centralized system, whereas in the proposed MAS, each agent gracefully degrades by using local kinematic estimates if the LSTM inference process fails.

5.4 Training and convergence

We used MSE loss and early stopping (patience of three) for training. The model trained for 18 epochs, and the best validation loss (MSE = 29.834) checkpoint was reached at epoch 16. The training and validation losses in Figure 14 show a consistent and monotonic decrease and remain close to each other. This indicates that the dropout regularization successfully guards against overfitting across the 20-route dataset. The narrow and consistent gap between the two curves is particularly important to the study’s claim of novelty: a single network serving 20 routes with distinct topologies does not overfit to each of the routes, proving that the learnable route embeddings account for enough of the route-specific variation within the model architecture. Figure 14 illustrates the progression of training (blue) and validation (red) MSE loss over 19 epochs. A green dashed line indicates the lowest validation loss of 29.766, attained at epoch 16.

5.5 Limitations

Even though the results are promising, there are still some problems:

• Part of the time features in the simulation are hardcoded (fixed hour/day encoding), which makes it hard to represent realistic daily changes.
• During the first sequence steps, ETA estimates may not be as accurate because the observation buffer is warming up.
• Route generalization is limited to trained routes unless retraining occurs.
• In extreme demand situations, the stop-skipping rule may make passengers at skipped stops wait longer.

5.6 Practical implications

The proposed approach may assist transit agencies in implementing adaptive headway control for busy routes. The 28.5% improvement over the baseline prediction demonstrates the effectiveness of route-specific deep learning for bus arrival time prediction. Thanks to the distributed nature of the MAS control layer, it does not require a central dispatcher: each bus agent responds to local states and predictions provided by the LSTM, which makes the system resilient to single points of failure.

6 Conclusion

This study presented a hybrid MAS to prevent real-time bus bunching using a multi-route LSTM predictor and a cooperative stop-skipping approach. The work makes four scientific contributions to the field of intelligent transportation systems.

First, a forward-looking imputation algorithm was introduced that recovers valid training labels from MTA SIRI feeds with up to 10-min polling gaps, reducing data loss from 95% to 10% and yielding 3.74 million labeled sequences—an order-of-magnitude improvement over conventional preprocessing with no hardware changes.

Second, a single multi-route LSTM with learnable eight-dimensional route embedding was trained across 20 heterogeneous NYC bus routes. It achieved a test MAE of 2.0247 min (95% CI: [2.017, 2.033] min), an RMSE of 5.5486 min, and a statistically significant 28.5% reduction over the route-aware historical baseline (p < 0.001, Wilcoxon signed-rank test, n = 802,085). All 20 routes improved simultaneously, with per route gains of 16.6%–39.7%, and prediction accuracy was found to be independent of the route data volume ( R 2 = 0.007 ), confirming that the shared embedding architecture generalizes across routes with heterogeneous traffic patterns.

Third, the pre-trained LSTM was embedded directly into the cognitive loop of each autonomous BusAgent in a Python mesa cooperative MAS—to our knowledge the first integration of production-trained deep learning inference into a transit MAS for anticipatory bunching prevention without a centralized dispatcher.

Fourth, quantitative simulation results on routes B82 and Bx9 over 120 steps demonstrate that the LSTM-augmented consensus condition (Equation 14) prevents full vehicle consolidation: the minimum inter-vehicle gap never collapsed to zero, the EXPRESS mode was activated in 19%–26% of steps concentrated around demand peaks, and the headway recovered above the safe-spacing threshold of 30 stop-units after 100% of EXPRESS interventions.

Future work should consider: (i) fairness constraints regarding demand to ensure passengers do not miss their destinations, (ii) modeling transfer penalties in transit networks of routes that cross each other, (iii) other external factors such as weather and incident feeds, and (iv) online learning to update route embedding across different seasons of the year.

Correction note

This article has been corrected with minor changes. These changes do not impact the scientific content of the article.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Yuming Lin, Tsinghua University, China

Reviewed by: Narges Rashvand, University of North Carolina at Charlotte, United States

Latha R, Saveetha Engineering College, India