1 Introduction↩︎

Public transit systems (PTSs) connect people with jobs, homes, and other daily activities and therefore play an important role in cities worldwide. Their service, however, is susceptible to frequent unplanned delays and disruptions. For example, [1] report an average of 75 incidents per day in the Chicago urban rail system, of which more than 75% last less than 5 minutes. These short suspensions may result from signal failures, passenger behavior, or infrastructure problems. Quantifying these effects is therefore important for assessing transit resilience.

In this study, we use two operational resilience metrics for PTSs, queue length and waiting time. Specifically, we model a PTS as a bulk service queue and aim to derive closed form expressions for the mean and variance of passenger queue length and waiting time at a station under random service suspensions. To this end, we derive a stability criterion for passenger queues under a renewal representation of suspension induced headways. The representation captures throughput loss when a nonpositive raw headway is converted to bunching and the lost separation is not recovered in a later headway. We also characterize the steady state distribution of passenger queues, which naturally leads to the quantification of these metrics.

Queueing behavior at a public transit station is commonly represented by a bulk service queue, in which a vehicle serves waiting passengers up to its available capacity [2]. Most transit applications analyze one station, while [3] extended the analysis to a route through a Markov representation of downstream passenger flows. These studies focus on normal operations. Models of bulk service under disruptions instead represent a server with operating and repair states and usually assume a fixed batch size [4], [5]. That assumption does not represent a transit route where the available boarding capacity varies with the onboard load and passenger alighting. Mapping a server breakdown to route operations is also difficult because passenger flows link the service capacities of consecutive stations along the route.

To fill these research gaps, we develop a bulk service queueing framework that links suspension induced headways, stochastic passenger arrivals, station queues, and the vehicle capacity propagated from upstream stations. The route calculation allows available boarding capacity to vary with onboard load and passenger alighting. Tracking these interactions station by station shows how short disruptions erode transit resilience through queue growth, waiting time variability, and station instability. These measures can inform future control and planning strategies for public transit.

Within this framework, we explicitly model random service suspensions in a single route PTS. Such a system represents a bus route or one direction of a rail line and is a basic element of more complex public transit networks. Each vehicle may experience random suspensions while traveling along the route. Section 3.2 explains the operational meaning of this assumption, shows how it corresponds to many real world situations, and presents the assumption as a first step toward a general incident representation for PTSs.

Our contributions are summarized as follows.

1. We develop an analytical bulk service queueing framework for quantifying line level public transit resilience under short random service suspensions. In particular, we represent service suspensions through a vehicle speed process with normal and disruption states. The speed process gives a direct operational meaning to the abstract notion of a server breakdown in public transit operations under short suspensions.

2. We develop an analytically tractable renewal approximation for headways under random service suspensions. The common marginal distribution is constructed as the positive part of a raw headway formed from the incident adjusted scheduled headway and the difference between two independent compound Poisson exponential variables. Under the approximation, headways are independent across vehicle runs, and departure epochs accumulate from these nonnegative headways. As a result, the analytical service sequence prevents overtaking without imposing a pathwise first in, first out (FIFO) recursion. We then use a censored normal approximation with a point mass at zero to obtain a closed form moment generating function and derive the probability generating function (PGF) and moments of passenger arrivals within a headway. This is a new analytical contribution to bulk service queueing under short service suspensions.

3. We develop a route level probability propagation model for inter station passenger flows. The model recursively links passenger alighting, available vehicle capacity, station queues, and departing vehicle loads. These links allow the resilience indicators to be computed sequentially at every station along the route.

4. We develop an interpolation based root solving algorithm that makes the analytical queueing indicators computationally implementable for realistic vehicle capacities. The algorithm finds all complex roots required by the transform based bulk service queue model, addressing a numerical bottleneck that standard single start nonlinear solvers and existing root search heuristics may miss.

The rest of this paper is organized as follows. Section 2 reviews the literature on the bulk service queue problem, random service disruptions, and queueing models for PTSs. Section 3 presents the model settings for a single route system with random service suspensions. Section 4 presents the major analytical results and their derivations. Section 5 provides numerical examples to illustrate the theoretical results and assesses the proposed approach using simulation. Section 6 concludes the paper and discusses future research directions.

2 Literature review↩︎

2.1 Bulk service queues↩︎

Bulk service queues serve customers in groups rather than individually. [6] formulated an early fixed capacity model with Poisson arrivals and an embedded Markov chain at service completion epochs. For transportation terminals, [2] developed transform methods for random passenger arrivals and bulk service. These models provide the queueing foundation for capacity constrained passenger boarding.

2.2 Random service disruptions↩︎

Models of queues with service breakdowns typically alternate between operating and repair states. [4] introduced this structure to a fixed batch bulk service queue, while [5] review the broader breakdown literature. These models are not directly applicable to a transit route because the service capacity at a station depends on the onboard load propagated from upstream stations.

2.3 Queueing models in public transit systems↩︎

Transit queueing analyses are usually conducted at the station level. Under regular service, suppose that passengers arrive according to a Poisson process, can board the first arriving vehicle, and face a deterministic headway \(H\). Their mean waiting time is then \[\begin{align} \mathbb{E}[W] = H/2, \label{eq95mean95wt951} \end{align}\tag{1}\] where \(H\) is the service headway and \(W\) is the passenger waiting time. With stochastic headways, [7] obtained \[\begin{align} \mathbb{E}[W] = \frac{1}{2} \cdot \left[\mathbb{E}[H] +\frac{\operatorname{Var}[H]}{\mathbb{E}[H]} \right], \label{eq95mean95wt952} \end{align}\tag{2}\] where \(\mathbb{E}[H]\) and \(\operatorname{Var}[H]\) are the expectation and variance of headways, respectively. When headways are deterministic, \(\operatorname{Var}[H]=0\), and the model reduces to Eq. 1 .

Both Eqs. 1 and 2 assume that every passenger can board the first vehicle. Once capacity becomes binding, however, some passengers are left behind and must wait for a later vehicle [8]. Transform based bulk service models account for this additional delay and provide closed form moments of queue length and waiting time under capacity constraints [9].

Extending these models from one station to a route requires the downstream propagation of passenger loads and available capacity. [10] characterized route level headway and passenger flow dynamics, and [3] combined this structure with a bulk service model. Their formulation assumes Erlang headways and does not jointly propagate headway correlation and vehicle capacity. Empirical headway distributions can differ substantially from the Erlang distribution [11]. Our model instead derives a suspension induced headway distribution and links it to downstream capacity propagation.

2.4 Service interruptions in public transit systems↩︎

Research on public transit interruptions includes impact analysis and operations control. Impact studies measure observed operational effects [1]. Control studies examine responses such as substitute bus services [12] and passenger path recommendations under demand or behavior uncertainty [13], [14]. The present study complements these approaches by deriving analytical resilience indicators, including stability conditions and station level moments of queue length and waiting time, within a route level capacity constrained model.

3 Model↩︎

This section defines the transit line, passenger flow primitives, and stochastic suspension process used in the analysis. We begin with the single route bulk service setting and describe how vehicles and passengers interact at each station. We then represent short service suspensions with a two state vehicle speed process and derive the stochastic headway used in the station level queueing analysis.

3.1 Single route public transit system and vehicle movements↩︎

Consider the single route PTS with \(N\) stations shown in Figure 1. Vehicles depart from a transportation hub, denoted by station 0, and visit stations 1 through \(N\) in sequence. At station \(n\), passenger arrivals follow a Poisson process with a fixed rate \(\lambda^{(n)}\) during the period of interest. When a vehicle arrives, each onboard passenger independently alights with probability \(\alpha^{(n)}\). Consequently, the number of alighting passengers at station \(n\) follows a binomial distribution. Poisson arrivals and binomial alighting are common assumptions in the public transit literature [10]. Because the paper focuses on short service suspensions, we assume that passengers continue to wait and do not renege. The assumption is consistent with evidence that passenger abandonment becomes more important during substantially longer delays [15]. Models of longer disruptions could relax the assumption by incorporating balking and reneging. Section 7 of the Electronic Companion (EC) appendix summarizes the modeling assumptions and their operational interpretations for the current setting.

Let \(l=1,2,\ldots\) denote the vehicle run index, with a smaller value indicating an earlier dispatch. Figure 2 summarizes the vehicle and passenger interactions at station \(n\) over time. Let \(t_A^{(n,l)}\) be the arrival time of vehicle \(l\) at station \(n\) and \(t_D^{(n,l)}\) its departure time. For \(l=2,3,\ldots\), we use \(H^{(n,l)}\) provisionally for the service interval between vehicle runs \(l-1\) and \(l\). Section 3.3 distinguishes its stationwise FIFO realization from the analytical renewal approximation. When a vehicle arrives, onboard passengers alight before queueing passengers board in their order of arrival. Let \(Q^{(n,l)}\) be the number of queueing passengers when vehicle \(l\) arrives at station \(n\), \(R^{(n,l)}\) the number left behind when it departs, and \(Y^{(n,l)}\) the number arriving between \(t_D^{(n,l)}\) and \(t_A^{(n,l+1)}\). By definition, \[\begin{align} Q^{(n,l+1)} = R^{(n,l)} + Y^{(n,l)}. \label{eq95relation95QRY} \end{align}\tag{3}\]

Following [2], we assume that dwell time is negligible relative to travel time and that no passenger arrives during dwell. Under this assumption, the number of passenger arrivals depends on the time between consecutive vehicles. For any service interval of length \(h\geq0\), the conditional arrival count satisfies \[\begin{align} Y^{(n,l)} \; |\; (\text{service interval}=h) \sim \text{{\fontfamily{qcs}\selectfont Poi}}(\lambda^{(n)} h). \label{eq95Y95dist} \end{align}\tag{4}\] Thus, \(Y^{(n,l)}\) counts passenger arrivals during the realized service interval.

To describe the same interaction from the vehicle side, let \(S^{(n,l)}\) be the number of available spaces after passengers alight from vehicle \(l\) at station \(n\), and let \(G^{(n,l)}\) be the number of passengers who remain onboard. By definition, \[\begin{align} G^{(n,l)} := C - S^{(n,l)}, \label{eq95Sn95Gn} \end{align}\tag{5}\] where \(C\) is the vehicle capacity. Let \(V^{(n,l)}\) denote the vehicle load when vehicle \(l\) departs station \(n\), which is also its load when it arrives at station \(n+1\). Conditional on \(V^{(n-1,l)}\), the number of passengers who alight from vehicle \(l\) at station \(n\) follows a binomial distribution. \[\begin{align} \left(V^{(n-1,l)} - G^{(n,l)}\right)\;|\;V^{(n-1,l)} \sim \text{{\fontfamily{qcs}\selectfont Bin}}(V^{(n-1,l)}, \alpha^{(n)}). \end{align}\]

3.2 Random service suspensions and vehicle speed profile↩︎

We assume that random service suspensions may occur as a vehicle travels through the system. The red curve in Figure 3 (a) shows a possible speed profile for vehicle \(l\) between stations \(n\) and \(n+1\) under such disturbances. Each incident slows or stops the vehicle. In bus systems, incidents may result from traffic congestion, accidents, vehicle engine problems, or driver and passenger behavior. In rail systems, they may result from signal or infrastructure failures and driver or passenger behavior. The speed profile provides a general representation of incidents, interruptions, suspensions, or disruptions that impede vehicle movement. We refer to the triggering event as an incident and to the resulting interruption of vehicle movement as a service suspension.

An actual vehicle speed profile under service suspensions can be complicated. For analytical tractability, we approximate it by the impulse function shown by the blue line in Figure 3 (a). The impulse approximation separates the trajectory into traveling and stopping phases, which we call the normal and disruption states. In the normal state, a vehicle travels at a constant speed. When an incident occurs, the vehicle stops immediately and enters the disruption state. Over a sufficiently small interval \(\Delta\), the probability of an incident is \(\gamma\Delta\). Once the vehicle enters the disruption state, the incident duration follows an exponential distribution with rate \(\theta\) and mean \(1/\theta\). The vehicle state therefore follows the two state Markov process in Figure 3 with state space \(\{\text{Normal},\text{Disruption}\}\). Equivalently, the durations of the normal and disruption states are exponentially distributed with transition rates \(\gamma\) for incident occurrence and \(\theta\) for service recovery.

The impulse function provides a first step toward a more general incident representation. Different incident categories could be represented by normal and disruption states with heterogeneous occurrence rates and durations, producing a richer speed profile.

3.3 Headway under random service suspensions↩︎

Under the assumed impulse function speed profile, all vehicles have the same fixed travel speed in the normal state. Consequently, when no incident occurs, all stations have the same deterministic headway \(\overline{H}\), which we call the nominal scheduled headway. To relate this headway to route planning, let \(\overline{E}\) denote the route cycle time required for a vehicle to travel from the transportation hub to the last station and return to the hub. Given route fleet size \(\overline{F}\), the scheduled headway is \(\overline{H} := \overline{E}/\overline{F}\). Random service suspensions increase the expected route cycle time and therefore require a planning response. A transit agency can maintain the scheduled headway \(\overline{H}\) by increasing the route fleet size, or it can maintain the same fleet size by increasing the scheduled headway. We consider the second response because it makes the effects of incidents on headway and service performance explicit, which is the focus of this paper.

Following this planning response, we assume that the transit agency has an estimate \(\overline{D}\) of the average incident delay in the cycle time. Let \(I^{(n,l)}\) be the total duration of all incidents while vehicle \(l\) travels from the transportation hub to station \(n\). Then \(\mathbb{E}[I^{(N,l)}]\) is the expected incident duration for travel from the hub to the last station \(N\). Assuming the same road conditions in both directions, the total estimated delay for the cycle trip is \[\begin{align} \overline{D} &:= 2 \cdot \mathbb{E}[I^{(N,l)}]. \label{eq95e95delay} \end{align}\tag{6}\] Some transit agencies may plan for a delay above the mean, such as the 85th percentile. In that case, the cycle delay could be formulated as a more general function of the incident duration distribution. We use the expected delay in Eq. 6 to keep the subsequent analysis transparent. After incorporating this delay, the incident adjusted scheduled headway is \(\overline{H}^{\text{Adj}} := (\overline{E}+\overline{D})/\overline{F} = \overline{H}+2\mathbb{E}[I^{(N,l)}]/\overline{F}\). The term \(2\mathbb{E}[I^{(N,l)}]/\overline{F}\) is therefore the scheduled headway adjustment associated with incidents. Because \(I^{(N,l)}\) is identically distributed across vehicle runs, the adjusted scheduled headway does not depend on the vehicle run index. The single route PTS then dispatches vehicles on time at intervals of \(\overline{H}^{\text{Adj}}\).

Let \(T^{(n)}\) be the travel time from the transportation hub to station \(n\) when there is no incident. This travel time is fixed under the constant speed assumption. Since dwell time is negligible, the arrival and departure times at a station are identical in the model. Taking the dispatch time of vehicle 1 as time zero, the free departure time of vehicle \(l\) at station \(n\) is \[\begin{align} \hat{t}_D^{(n,l)} := (l-1)\overline{H}^{\text{Adj}}+T^{(n)}+I^{(n,l)}. \label{eq95free95departure} \end{align}\tag{7}\] The free departure time in Eq. 7 is the time that would result if vehicle order were unrestricted. Under this construction, a sufficiently long incident affecting vehicle \(l-1\) can yield \(\hat{t}_D^{(n,l)}<\hat{t}_D^{(n,l-1)}\). Vehicle \(l\) would then appear to overtake vehicle \(l-1\), although vehicle \(l-1\) was dispatched first.

To rule out such reversals, we assume that vehicles serve each station in dispatch order and that overtaking is not allowed. The FIFO rule represents rail lines where passing is infeasible and bus services where passing is prohibited or uncommon at the route scale. It also keeps the vehicle index aligned with the service order observed by waiting passengers. If the free departure time of vehicle \(l\) precedes the realized departure time of vehicle \(l-1\), vehicle \(l\) is held until vehicle \(l-1\) departs. The realized departure times at station \(n\) therefore satisfy \[\begin{align} t_D^{(n,1)} := \hat{t}_D^{(n,1)}, \quad t_D^{(n,l)} := \max\left\{\hat{t}_D^{(n,l)},t_D^{(n,l-1)}\right\}, \qquad \forall l\geq 2. \label{eq95fifo95departure} \end{align}\tag{8}\] Under Eq. 8 , the stationwise FIFO headway is \(H_{\mathrm{FIFO}}^{(n,l)}:=t_D^{(n,l)}-t_D^{(n,l-1)}=\max\{0,\hat{t}_D^{(n,l)}-t_D^{(n,l-1)}\}\). This headway depends on the realized departure time of the preceding vehicle and can therefore inherit delays from several earlier vehicles at the same station. Direct use of this recursion would therefore require the departure times of multiple vehicle runs to be analyzed jointly.

To avoid this recursive dependence, we impose the following approximation before constructing the analytical headway distribution.

Although the incident draws are independent of the exogenous Poisson passenger arrival process, the number of passengers arriving within an interval depends on its length through Eq. 4 . Separately, on time dispatch supports the assumed independence of the analytical headway from prior queue and vehicle load states. Independence across vehicle runs is an analytical approximation that is most defensible when the total incident delay of one vehicle rarely extends across several scheduled headways. The condition is consistent with service suspensions that are short relative to the scheduled headway. The approximation is also suitable for headway based operations in which holding or dispatch control resets the next service interval from the current departure [16]. It becomes less suitable when a long delay propagates through many subsequent vehicles. Section 5 evaluates this limitation using a simulation that retains the recursive FIFO dynamics in Eq. 8 .

Under Assumption 1, the common marginal distribution is obtained from the difference between the free departure times of a generic pair of consecutive vehicles. Accordingly, subtracting the two free departure times in Eq. 7 gives the raw headway \[\begin{align} \widehat{H}^{(n,l)} := \hat{t}_D^{(n,l)}-\hat{t}_D^{(n,l-1)} = \overline{H}^{\text{Adj}} + I^{(n,l)} - I^{(n,l-1)}, \qquad l=2,3,\ldots. \label{eq95headway95raw} \end{align}\tag{9}\] The deterministic term in Eq. 9 is \(\overline{H}^{\text{Adj}}=\overline{H}+2\mathbb{E}[I^{(N,l)}]/\overline{F}\). For one marginal headway calculation, the two incident durations are independent. In a physical vehicle sequence, adjacent raw headways share one incident duration and are therefore dependent. The renewal model draws a new pair of incident durations for each analytical headway. If the difference between the two durations is sufficiently negative, \(\widehat{H}^{(n,l)}\) becomes negative and cannot represent elapsed time between service events. We therefore define the analytical headway as \[\begin{align} H^{(n,l)} := \max\{0,\widehat{H}^{(n,l)}\}, \qquad l=2,3,\ldots. \label{eq95renewal95headway} \end{align}\tag{10}\] Equation 10 is not algebraically equivalent to the recursive FIFO rule in Eq. 8 . A nonpositive raw headway becomes zero and represents bunching, but the next headway is redrawn instead of inheriting the accumulated delay. The approximation preserves vehicle order through nonnegative analytical service intervals, but it omits serial dependence and the propagation of accumulated delay to later vehicles.

4 Analysis↩︎

The analysis derives the stability condition in Proposition [prop95stability] and the mean and variance of passenger queue length and waiting time in Propositions [prop95q95length] and [prop95waiting]. Figure 4 places these station level results within the route calculation by summarizing the propagation of the distributions of \(S^{(n,l)}\), \(V^{(n,l)}\), and \(Q^{(n,l)}\). The calculation consists of three linked updates.

• Given the distribution of \(V^{(n-1,l)}\), calculate the distribution of the available spaces \(S^{(n,l)}\) after alighting at station \(n\), as shown in Section 4.1.

• Given the distribution of \(S^{(n,l)}\), calculate the distribution of the queue \(Q^{(n,l)}\) and the moments of queue length and waiting time, as shown in Section 4.3.

• Given the distributions of \(S^{(n,l)}\) and \(Q^{(n,l)}\), calculate the distribution of the departing vehicle load \(V^{(n,l)}\), as shown in Section 4.2.

Together, these components yield the distributions of \(S^{(n,l)}\), \(Q^{(n,l)}\), and \(V^{(n,l)}\) for all \(n=1,\ldots,N\). The recursion starts from \(V^{(0,l)}=0\), which represents an empty vehicle departing the hub for the first station. The analysis focuses on steady state distributions as \(l\to\infty\).

4.1 Available vehicle space distribution↩︎

This subsection derives the steady state distribution of \(S^{(n,l)}\) from that of \(V^{(n-1,l)}\). Define \(v_k^{(n,l)} := \mathbb{P}(V^{(n,l)}=k)\), \(s_k^{(n,l)} := \mathbb{P}(S^{(n,l)}=k)\), and \(g_k^{(n,l)} := \mathbb{P}(G^{(n,l)}=k)\) for \(k=0,1,\ldots,C\). When the steady state probabilities exist, define \(v_k^{(n)} := \lim_{l\to\infty}v_k^{(n,l)}\), \(s_k^{(n)} := \lim_{l\to\infty}s_k^{(n,l)}\), and \(g_k^{(n)} := \lim_{l\to\infty}g_k^{(n,l)}\). Let \(V^{(n)}\), \(S^{(n)}\), and \(G^{(n)}\) denote random variables with these steady state distributions. Section 4.4 establishes the stability condition under which the distributions exist.

4.2 Vehicle load distribution↩︎

This subsection states the steady state update for \(V^{(n,l)}\) in terms of the distributions of \(G^{(n,l)}\) and \(Q^{(n,l)}\); Section 4.3 subsequently derives the required queue distribution. Combining these two marginal distributions requires the following approximation for route propagation.

Assumption 2 is exact at the first station because vehicles arrive from the hub without passengers and \(G^{(1)}=0\). At a downstream station, the remaining onboard load results from boarding at upstream stations followed by binomial alighting, whereas the local queue results from the station specific Poisson arrival process and passengers previously left behind at that station. These quantities therefore arise from different passenger flow components. Independent renewal headways, alighting at intermediate stations, and repeated queue clearance under stable operation weaken the dependence induced by common headways and upstream capacity constraints. The approximation is most appropriate for short suspensions and stable stations, and it becomes weaker when long disruptions or persistent oversaturation jointly produce high onboard loads and long local queues. It permits the two marginal distributions to be propagated separately along the route. Section 5 assesses the approximation through recursive simulation.

Define \(q_k^{(n,l)} := \mathbb{P}(Q^{(n,l)}=k)\) and, when the steady state distribution exists, \(q_k^{(n)} := \lim_{l\to\infty}q_k^{(n,l)}\). Let \(Q^{(n)}\) denote a random variable with probabilities \(\mathbb{P}(Q^{(n)}=k)=q_k^{(n)}\). The vector \(q_{0:C-1}^{(n)} := [q_0^{(n)},\ldots,q_{C-1}^{(n)}]\in\mathbb{R}^{C}\) contains its first \(C\) probabilities.

4.3 Queueing analysis at a station↩︎

Given the distribution \(s^{(n)} := [s_0^{(n)},\ldots,s_C^{(n)}]\in\mathbb{R}^{C+1}\) of \(S^{(n)}\), this subsection identifies the characteristic roots and determines \(q_{0:C-1}^{(n)}\) and the mean and variance of passenger queue length and waiting time. The transform calculations in this subsection are conditional on station stability, so the steady state queue and its PGF are well defined. Section 4.4 establishes the stability condition independently from the transient queue recursion. For compact notation, let \(Y^{(n)}\) denote the steady state number of passenger arrivals within a headway at station \(n\), and write \(\overline{S}^{(n)}=\mathbb{E}[S^{(n)}]\). The quantities \(\overline{\overline{S}}^{(n)}\) and \(\overline{\overline{\overline{S}}}^{(n)}\) denote the second and third central moments of \(S^{(n)}\), respectively, and the same convention applies to \(Y^{(n)}\).

4.3.1 Probability generating function and characteristic roots.↩︎

We first derive the PGF of the steady state queue \(Q^{(n)}\) and then identify the characteristic roots required by the queue probability calculation in Section 4.3.2. Throughout this station level derivation, \(Q(z)\), \(Y(z)\), \(\mathrm{\small Num}(z)\), and \(\mathrm{\small Den}(z)\) refer to station \(n\), whose superscript is suppressed.

Let \(\mathrm{\small Num}(z)\) and \(\mathrm{\small Den}(z)\) denote the numerator and denominator of \(Q(z)\) in Eq. 15 . Because \(Q(z)\) is a PGF, it is analytic for \(|z|<1\) and bounded for \(|z|\leq1\). Every nonunit root of \(\mathrm{\small Den}(z)\) in the closed unit disk must therefore be canceled by a root of \(\mathrm{\small Num}(z)\) with the same multiplicity. The unit root is also a root of \(\mathrm{\small Num}(z)\) because direct substitution gives \(\mathrm{\small Num}(1)=0\). Since \(\mathrm{\small Num}(z)\) is a polynomial of degree at most \(C\), the \(C\) characteristic roots counted according to multiplicity determine its factorization up to a scale factor. Section 4.3.2 uses this factorization to solve for \(q_0^{(n)},\ldots,q_{C-1}^{(n)}\).

To locate the required roots, we multiply \(\mathrm{\small Den}(z)\) by \(Y(z)\) to obtain the characteristic function \[\begin{align} \mathcal{F}^{(n)}(z) := z^C-Y(z)\sum_{u=0}^{C}s_u^{(n)}z^{C-u}. \label{eq95characteristic95function} \end{align}\tag{16}\] The result of [17] implies that \(\mathcal{F}^{(n)}(z)\) has exactly \(C\) roots in the closed unit disk when roots are counted according to multiplicity. Because \(s_C^{(n)}>0\) and \(Y(0)>0\), zero is not a root of \(\mathcal{F}^{(n)}(z)\). At every nonzero root, \(Y(z)\) is nonzero, so the same root and multiplicity apply to \(\mathrm{\small Den}(z)\). In particular, \(z=1\) is a simple root because \([\mathcal{F}^{(n)}]'(1)=\overline{S}^{(n)}-\overline{Y}^{(n)}>0\). Interior roots need not be distinct [18]. We therefore solve \(\mathcal{F}^{(n)}(z)=0\) for the distinct root locations \(\zeta_1,\ldots,\zeta_J\) in the closed unit disk, where \(J\) denotes the number of distinct locations. Let \(m_1,\ldots,m_J\) denote their corresponding multiplicities, which satisfy \(\sum_{j=1}^{J}m_j=C\). For use in the queue calculation, we list the roots as \(z_0^*,\ldots,z_{C-1}^*\) according to multiplicity and adopt the convention \(z_0^*=1\). Section 4.3.2 uses this list to determine \(q_0^{(n)},\ldots,q_{C-1}^{(n)}\), while Section 4.3.6 presents the numerical procedure for constructing it.

4.3.2 Queue length distribution.↩︎

Section 4.3.1 identifies the \(C\) characteristic roots required to factor the numerator of \(Q(z)\). We now translate that factorization into the queue probabilities \(q_0^{(n)},\ldots,q_{C-1}^{(n)}\). First, the normalization condition \(Q(1)=1\) determines the scale factor. Matching the resulting polynomial coefficients then yields the required probabilities.

The condition \(s_C^{(n)}>0\) holds at the first station because every vehicle arrives empty. At a downstream station, it holds whenever there is a positive probability that all onboard passengers alight or that the vehicle arrives empty. In a degenerate case with \(s_C^{(n)}=0\), the station level transform can instead use the largest \(u\) for which \(s_u^{(n)}>0\) as the effective maximum service capacity without changing the queue process.

4.3.3 Analytical formulation of mean and variance of queue length and waiting time.↩︎

Once \(q_0^{(n)},\ldots,q_{C-1}^{(n)}\) have been obtained, the PGF \(Q(z)\) is fully determined. We can therefore recover the expectation and variance of queue length at station \(n\) from its first two derivatives as follows. \[\begin{align} \mathbb{E}[Q^{(n)}] &= \sum_{k = 0} ^ \infty k q_k^{(n)} =\left.\frac{dQ(z)}{d{z}}\right\vert_{z = 1},\\ \operatorname{Var}[Q^{(n)}] &= \mathbb{E}[(Q^{(n)})^2] - \mathbb{E}[Q^{(n)}]^2 =\left.\frac{d^2Q(z)}{d{z^2}}\right\vert_{z = 1} + \mathbb{E}[Q^{(n)}] - \mathbb{E}[Q^{(n)}]^2. \end{align}\]

We next express waiting time in terms of the queue length observed at an arbitrary time.

Under the Poisson arrival and first come, first served assumptions, Eq. 19 applies Little’s law, and Proposition [prop95waiting] follows from [9].

4.3.4 Headway distribution.↩︎

As defined in Eq. 10 , the analytical headway is \(H^{(n,l)} := \max\{0,\widehat{H}^{(n,l)}\}\). For the steady state analysis, we suppress the vehicle run index and write the common marginal variables as \(H^{(n)}\) and \(\widehat{H}^{(n)}\). We derive the distribution of \(H^{(n)}\) in two steps. We first characterize the distribution of the raw headway \(\widehat{H}^{(n)}\) and then take its positive part to obtain the distribution of \(H^{(n)}\). This distribution is required in Propositions [prop95solve95q] to [prop95waiting] because Eq. 4 implies that \(Y^{(n)}\mid H^{(n)}\) is Poisson with parameter \(\lambda^{(n)}H^{(n)}\).

Equation 9 shows that the raw headway depends on the difference between two independent copies of the total incident duration for travel from the hub to station \(n\). We therefore begin the derivation by characterizing the distribution of \(I^{(n,l)}\).

The moment generating function (MGF) of a compound Poisson exponential variable is \[\begin{align} M_{I^{(n,l)}}(t) = \mathbb{E}[e^{tI^{(n,l)}}] = e^{\gamma T^{(n)} (\frac{\theta}{\theta - t} - 1)}, \qquad t < \theta. \end{align}\] Similarly, the MGF of \(-I^{(n,l-1)}\) is \[\begin{align} M_{-I^{(n,l-1)}}(t) = \mathbb{E}[e^{-tI^{(n,l-1)}}] = e^{\gamma T^{(n)} (\frac{\theta}{\theta + t} - 1)}, \qquad t > -\theta. \end{align}\] The MGF of \(I^{(n,l)}\) gives \(\mathbb{E}[I^{(N,l)}] = \gamma T^{(N)}/\theta\). Using this expectation in Eq. 9 and the independence of the two incident duration copies gives the raw headway MGF in Proposition [prop95headway95distribution1].

From the MGF of \(\widehat{H}^{(n)}\), we obtain the mean and variance of the raw headway as \[\begin{align} \mathbb{E}[\widehat{H}^{(n)}] &= \overline{H} + \frac{2 \gamma T^{(N)}}{\theta\overline{F}},\\ \operatorname{Var}[\widehat{H}^{(n)}] &= \frac{4 T^{(n)} \gamma}{\theta^{2}}. \end{align}\]

Having characterized the raw headway, the second step applies the positive part transformation in Eq. 10 . A direct MGF representation of the resulting headway involves nested integrals over the two compound Poisson distributions. Although this representation is exact, it does not reduce to a tractable finite expression for general parameter values. Moreover, the root calculations would require repeated numerical integration or infinite series evaluation at every complex argument. To avoid this computational burden, we first approximate the raw headway \(\widehat{H}^{(n)}\) by a normal distribution with the same mean and variance and then take its positive part.

We use the normal approximation for two reasons. First, \(I^{(n,l)}\) is a random sum of independent exponential durations. A normal approximation becomes more accurate as the expected incident count increases, consistent with the setting of frequent short disturbances. A related central limit approximation for headway disturbances appears in [16]. Second, let \(\widehat{H}^{(n)}_{\text{Normal}}\) denote a normal variable with the same mean and variance as \(\widehat{H}^{(n)}\). Both variables are symmetric and have a third central moment of zero. Moreover, the MGF of \(\widehat{H}^{(n)}_{\text{Normal}}\) is \[\begin{align} M_{\widehat{H}^{(n)}_{\text{Normal}}}(t) = e^{t(\overline{H} + \frac{2 \gamma T^{(N)}}{\theta\overline{F}})} e^{\gamma T^{(n)}(\frac{2t^2}{\theta^2})}, \end{align}\] which has the same second order expansion around \(t=0\) as the exact MGF in Eq. 22 . The fourth central moment is \(48(T^{(n)}\gamma)^2/\theta^4\) for \(\widehat{H}^{(n)}_{\text{Normal}}\) and \([48(T^{(n)}\gamma)^2+48T^{(n)}\gamma]/\theta^4\) for \(\widehat{H}^{(n)}\). The exact raw headway may therefore have heavier tails. Figure 5 numerically compares the two distributions for several values of \(T^{(n)}\), \(\theta\), and \(\gamma\). The histograms are generated by sampling from the associated exponential and Poisson distributions. The normal distribution provides a close approximation in these examples while preserving the expected difference in tail behavior.

We next take the positive part of \(\widehat{H}^{(n)}_{\text{Normal}}\). Denote the resulting headway approximation by \(H^{(n)}_{\text{Normal}}\). It has support \([0,+\infty)\) and a probability mass at zero.

For compact notation, let \(\mu:=\overline{H}+2\gamma T^{(N)}/(\theta\overline{F})\) and \(\sigma:=2\sqrt{T^{(n)}\gamma}/\theta\) denote the mean and standard deviation of \(\widehat{H}^{(n)}_{\text{Normal}}\). Differentiating the MGF and using \(1-\Phi(-\mu/\sigma)=\Phi(\mu/\sigma)\) then gives the corresponding mean and variance of the approximated headway. \[\begin{align} \mathbb{E}[H^{(n)}_{\text{Normal}}] &= \mu \cdot \Phi\left(\frac{\mu}{\sigma} \right) +\sigma\cdot \phi\left(\frac{-\mu}{\sigma}\right), \label{eq95mean95normal95tr95hdw}\\ \operatorname{Var}[H^{(n)}_{\text{Normal}}]&= \mu \sigma \phi\left(\frac{-\mu}{\sigma}\right) +\Phi\left({\frac{\mu}{\sigma}}\right) \left(\mu^{2} + \sigma^{2}\right) \;- \left( \mu \Phi\left(\frac{\mu}{\sigma}\right) + \phi\left(\frac{-\mu}{\sigma}\right) \sigma\right)^{2}, \end{align}\tag{24}\] where \(\phi(\cdot)\) is the probability density function of a standard normal distribution. For \(\gamma=0\), these expressions are interpreted by continuity, giving mean \(\overline{H}\) and variance zero. The following proposition characterizes how incident frequency \(\gamma\) and mean incident duration \(1/\theta\) affect the expected headway under the normal approximation.

Proposition [prop95headway95increase] supports the comparative stability analysis in Section 4.4.

4.3.5 Distribution of passenger arrivals.↩︎

We can now return to the passenger arrival distribution that enters the queueing model. Recall that \(H^{(n)}_{\text{Normal}}\) is obtained by taking the positive part of the normal approximation to the raw headway \(\widehat{H}^{(n)}\). Averaging the conditional Poisson PGF over this headway distribution gives the following result for \(Y^{(n)}\).

Differentiating the PGF in Proposition [prop95arrival95pax] at \(z=1\) gives the first three factorial moments. We then convert these quantities to \(\overline{Y}^{(n)}\), \(\overline{\overline{Y}}^{(n)}\), and \(\overline{\overline{\overline{Y}}}^{(n)}\). The first moment is shown below. The second and third central moments follow from the same differentiation and conversion, but their longer expressions are omitted for readability. \[\begin{align} \overline{Y}^{(n)} &= \left(\mu \cdot \Phi\left(\frac{\mu}{\sigma} \right) +\sigma\cdot \phi\left(\frac{-\mu}{\sigma}\right)\right) \cdot \lambda^{(n)}. \end{align}\] When \(\gamma=0\), this expression is interpreted by continuity and gives \(\overline{Y}^{(n)}=\lambda^{(n)}\overline{H}\).

4.3.6 Solving for the roots.↩︎

With \(Y(z)\) now specified, we turn to the \(C\) roots of \(\mathcal{F}^{(n)}(z)\) in the closed unit disk that are required by the queue length and waiting time formulas. A nonlinear solver with a single initialization generally returns only one root, and repeated initializations may converge to the same location. The expression for \(Y(z)\) makes this problem more difficult for the vehicle capacities considered here. We therefore develop a certified interpolation search that combines the structured initialization of [9] with deterministic interpolation and multiplicity certification. The algorithm supplies the required roots and reports a diagnostic failure when certification is unsuccessful.

Step 1. Initial root search. Following the structured angular initialization of [9], we place \(C\) initial values at equally spaced angles \(2\pi(k+1/2)/C\) for \(k=0,\ldots,C-1\). Let \(\rho^{(n)}:=\overline{Y}^{(n)}/\overline{S}^{(n)}\) denote utilization at station \(n\). Their common radius is \(1-0.5\rho^{(n)}\), bounded between 0.05 and 0.98. From each initial value, a local solver evaluates the real and imaginary parts of \(\mathcal{F}^{(n)}(z)\) in Eq. 16 . Candidates that pass the convergence, residual, and unit disk checks are retained together with their conjugates and grouped into distinct root clusters. Polar extrapolation from consecutive clusters provides additional initial values for the initial search.

The root locations often form an oval pattern. Figure 6 illustrates this pattern for different utilization ratios \(\rho^{(n)}\) and available capacity distributions. As \(\rho^{(n)}\) increases, the root locations generally move from a nearly circular pattern toward a more elongated pattern.

Step 2. Interpolation search. The initial search may return fewer than \(C\) roots counted according to multiplicity. We order the distinct root clusters by complex argument and interpret adjacent pairs circularly. At interpolation level \(L\), the algorithm interpolates both the modulus and the unwrapped argument of every adjacent pair at \(L-1\) locations. It places initial values at the interpolated modulus and at radial factors 0.85 and 1.15, with the modulus bounded by one. These radial values prevent a straight chord between two roots from moving all initial values toward the origin. The algorithm applies the local solver, refines each successful candidate, and updates the root clusters. If an iteration finds no new root location, it increases \(L\). Because the coefficients of \(\mathcal{F}^{(n)}(z)\) are real, the conjugate of every nonreal candidate is also evaluated. Figure 7 compares the distinct locations from the initial search with the root set after interpolation and certification.

Step 3. Multiplicity certification and termination. Finding \(C\) distinct locations is not a valid stopping rule because a root may have multiplicity greater than one. Let \(\mathcal{Z}^{(k)} := \{\zeta_0^{(k)},\ldots,\zeta_{M_k-1}^{(k)}\}\) contain the \(M_k\) distinct root clusters after iteration \(k\). For cluster \(j\), define the counterclockwise circular contour \(\Gamma_j^{(k)} := \{\zeta_j^{(k)}+r_j^{(k)}\exp(\varphi\sqrt{-1})\mid0\leq\varphi\leq2\pi\}\). Its radius is smaller than half the distance to the nearest other cluster and no greater than 0.05. We use 0.05 when only one cluster is available. Thus, \(\Gamma_j^{(k)}\) is tied to candidate \(\zeta_j^{(k)}\) and separates it from other candidates. We use the contour only when \(\mathcal{F}^{(n)}(z)\neq0\) for every \(z\) on \(\Gamma_j^{(k)}\). \[\begin{align} m_j^{(k)} := \frac{1}{2\pi\sqrt{-1}}\int_{\Gamma_j^{(k)}}\frac{[\mathcal{F}^{(n)}]'(z)}{\mathcal{F}^{(n)}(z)}\,dz. \label{eq95root95multiplicity} \end{align}\tag{25}\] Because \(\mathcal{F}^{(n)}(z)\) has no poles, \(m_j^{(k)}\) equals the number of roots enclosed by \(\Gamma_j^{(k)}\), counted according to multiplicity. It is also the winding number of \(\mathcal{F}^{(n)}(z)\) around the origin along this contour. A value of one certifies a simple root in cluster \(j\), while a value of two certifies a root of multiplicity two. We compute this value from the unwrapped change in the argument of \(\mathcal{F}^{(n)}(z)\) and repeat the calculation with finer discretizations and smaller contours. A cluster is certified when these calculations return the same positive integer.

The algorithm terminates when every cluster is certified and the sum of the certified multiplicities equals \(C\). Otherwise, it returns to Step 2 until a search limit is reached and then reports a diagnostic failure. This certification follows numerical treatments of transform equations with possible multiple roots [18].

The complete procedure is provided in Algorithm 13 of the EC appendix.

4.4 Stability condition↩︎

The preceding transform calculations characterize the queue conditional on the existence of a steady state. We now establish when this steady state exists directly from the transient queue recursion. We call a station stable when its queue process admits a proper stationary probability distribution.

Assumption 3 separates the propagation of available capacity along the route from the evolution of the queue at an individual station. Upstream passenger flows determine \(s^{(n)}\), so available capacity remains endogenous across stations. After \(s^{(n)}\) is propagated to station \(n\), the model uses its complete distribution to generate repeated service opportunities that are exogenous to the local queue, thereby carrying the effect of upstream crowding into the stationary single station analysis.

Assumption 3 is exact at the first station because every arriving vehicle has capacity \(C\). Under persistent upstream saturation, the approximation is also reasonable because consecutive vehicles depart at capacity and independent binomial alighting determines the available capacity at station \(n\). Outside these cases, independent renewal headways generate new passenger arrivals for successive vehicles, while alighting at intermediate stations weakens the persistence of upstream loads. These mechanisms limit dependence in available capacity when suspensions are short and rarely affect several vehicle runs. When long disruptions or common operating conditions sustain correlated vehicle loads, the approximation becomes weaker. Section 5 assesses it through recursive simulation.

The proof in the EC appendix writes the residual queue as a Lindley recursion [19] and applies the stability result of [20]. It therefore establishes the existence of the stationary distribution without using \(q_0^{(n)}\) or any other quantity derived under stationarity.

Proposition [prop95stability] states that station \(n\) is stable when the mean number of passenger arrivals within a headway is smaller than the mean available capacity of an arriving vehicle after alighting. Proposition [prop95headway95increase] shows that a higher incident rate or a longer mean incident duration increases \(\mathbb{E}[H^{(n)}_{\text{Normal}}]\) under the conditions stated there, making instability more likely. The result therefore quantifies the throughput loss induced by incidents.

4.5 Summary of calculation procedure↩︎

Algorithm 8 combines the preceding analytical results into a sequential calculation over the \(N\) stations of the route.

5 Numerical experiments↩︎

This section illustrates the analytical framework on a bus route with ten stations and assesses the resulting resilience indicators against recursive FIFO simulation. The experiments isolate the effects of vehicle capacity, incident frequency, incident duration, scheduled headway, and demand, with particular attention to stations with limited capacity slack. All analytical results use the certified interpolation root search described in Section 4.3.6. Section 8 of the EC appendix reports the numerical settings used for the search.

5.1 Experimental design↩︎

We use an example bus route adapted from [21]. Table 1 reports the attributes of its ten stations. Figure 9 shows the route layout. Travel time without incidents is 5 minutes between consecutive stations, total cycle time is \(\overline{E}=100\) minutes, and travel time from the hub to the last station is \(T^{(N)}=50\) minutes.

To test the sensitivity of resilience indicators to different parameters, we consider different values of \(C\), \(\theta\), \(\gamma\), \(\overline{H}\), and demand in Table 2. The demand factor multiplies the base arrival rates \(\lambda^{(n)}\) in Table 1. Keeping the route cycle time at \(\overline{E}=100\) minutes, we select the integer fleet sizes \(\overline{F}\in\{50,25,14\}\), which yield \(\overline{H}=\overline{E}/\overline{F}\in\{2,4,50/7\}\) minutes. The reference scenario sets \(\overline{F}=25\) and \(\overline{H}=4\) minutes; each test varies one parameter and retains all other reference values.

5.2 Performance indicators↩︎

Figure 10 shows the mean and standard deviation of queue length under the tested scenarios. The patterns are consistent with the passenger arrival and alighting rates. Expected queues are relatively high at stations 2 and 8, whereas the expected queue at station 10 is zero because \(\lambda^{(10)}=0\).

Figure 10 (a) shows limited sensitivity to vehicle capacity within the tested range because capacity is not fully utilized under the reference setting. Figure 10 (b) shows the effect of incident occurrence rate \(\gamma\). When \(\gamma=0\), the expected queue length at station 8 is 4.5 passengers. As \(\gamma\) increases, expected queue lengths and their standard deviations rise. At \(\gamma=1/3\) per minute, the expected queue at station 8 reaches 7.2 passengers. Incident duration produces a similar pattern in Figure 10 (c). When mean incident duration is 30 seconds and \(\theta=2\), \(\mathbb{E}[Q^{(8)}]=5.0\). When mean incident duration is 2 minutes and \(\theta=1/2\), this expectation increases to 9.3. Both effects are larger at crowded stations.

Figure 10 (d) shows that a longer scheduled headway lowers the service rate and increases the expected queue. As \(\overline{H}\) increases from 2 minutes to \(50/7\) minutes, the expected queue at station 8 increases from 4.1 to 9.9 passengers. Figure 10 (e) shows a similar pattern for demand. Increasing the demand factor from 0.5 to 1.0 raises the expected queue at station 8 from 4.1 to 8.2 passengers. The effects of \(\overline{H}\) and demand are more uniform across stations than the effects of incident frequency and duration.

Figure 11 shows the mean and standard deviation of passenger waiting time. Station 10 is omitted because its zero passenger arrival rate produces no waiting time observations. Downstream stations generally have higher waiting time means and standard deviations because headway variability accumulates along the route. At congested stations such as stations 3 and 8, capacity constraints leave some passengers behind and further increase waiting time.

Figure 11 (a) shows that capacity has a limited effect on waiting time within the tested range. Figures 11 (b) and 11 (c) show the effects of \(\gamma\) and \(\theta\). Higher values of \(\gamma\) and \(1/\theta\) increase the expected headway and, in turn, mean waiting time at stations with passenger demand, particularly at crowded stations. When \(\gamma=0\), there are no incidents, left behind passengers, or headway irregularities. Expected waiting time at stations 1 through 9 is then \(\overline{H}/2=2\) minutes under the reference headway of 4 minutes. When \(\gamma\) increases to \(1/5\), station 3 has left behind passengers and its expected waiting time increases to 3.5 minutes. When \(\theta\) decreases from 2 to \(1/2\), expected waiting time at station 8 increases from 3.0 to 8.8 minutes.

A larger scheduled headway directly increases waiting time at stations with passenger demand, and passengers are left behind at stations 3 and 8 when \(\overline{H}=50/7\) minutes. Demand does not alter the headway distribution, so its effect on waiting time remains limited until capacity constraints leave passengers behind. At station 8, increasing the demand factor from 0.5 to 1.0 raises expected waiting time only from 4.77 to 4.83 minutes within the tested range.

5.3 Comparison between simulation and analytical results↩︎

To assess the renewal headway, marginal route decoupling, and normal distribution approximations, we develop a recursive FIFO simulation model for the expectation and variance of queue length and waiting time. Unlike the analytical renewal model, the simulation carries the effect of a delayed vehicle to its successors through realized departure times and jointly evolves passenger queues and vehicle loads; see Section 10 of the EC appendix.

The assessment covers seven settings from Table 2, including the reference setting and variations in capacity, headway, demand, incident frequency, and incident duration. For each setting, we conduct 20 independent replications of 10,000 vehicle runs using distinct random seeds and exclude the first 10% of vehicle runs. Figure 12 compares the analytical results with the means across replications for the reference setting. The whiskers give 95% confidence intervals (CIs) based on variation across replications. Station 10 has \(\lambda^{(10)}=0\), so no passenger arrives there and no waiting time observation exists. We therefore report waiting time through station 9. All error calculations also use stations 1 through 9 so that the four performance measures have common support.

Table ¿tbl:tab95sim95validation? reports the relative \(\ell_1\) error in percentage terms and the average half width (HW) of the station CIs. For a given performance metric, \(x_{\mathrm{Sim}}^{(n)}\) denotes the mean estimate across the 20 independent simulation replications at station \(n\), and \(x_{\mathrm{Ana}}^{(n)}\) denotes the corresponding result obtained from the analytical model. The relative \(\ell_1\) error is calculated as \(100\sum_{n=1}^{9}|x_{\mathrm{Sim}}^{(n)}-x_{\mathrm{Ana}}^{(n)}|/\sum_{n=1}^{9}|x_{\mathrm{Ana}}^{(n)}|\). Under the reference setting, the queue length errors are 5.64% for the mean and 3.27% for the standard deviation. The corresponding waiting time errors are 3.87% and 1.26%.

Across all tested settings, the waiting time errors remain below 8.21%. The largest discrepancy is the 16.58% error in mean queue length when the mean incident duration is 2 minutes. This setting is also where the short suspension premise and the independent renewal approximation are most strained. The average confidence interval half width over the seven settings is at most 0.013 passengers for mean queue length, 0.025 passengers for its standard deviation, 0.039 minutes for mean waiting time, and 0.052 minutes for its standard deviation. The sampling uncertainty is therefore small relative to the main analytical differences.

The results provide numerical evidence for the accuracy of the approximations within the tested parameter range, but they do not establish exact equivalence between the analytical renewal process and the recursive FIFO process. One source of discrepancy is the normal approximation of the raw headway distribution. As shown in Figure 5, the compound Poisson based raw headway concentrates more probability density near the mean than the normal approximation and therefore assigns a lower probability to large deviations from the scheduled headway. A second source is the independent redraw of each analytical headway, whereas the recursive FIFO simulation carries delay forward.

6 Conclusion and discussion↩︎

This paper proposes a stochastic framework for quantifying the resilience of PTSs under short random service suspensions. Specifically, we introduce an independent renewal approximation for nonnegative headways and derive stability conditions and closed form expressions for the mean and variance of queue length and waiting time at each station. The common headway marginal preserves the incident effects implied by a two state vehicle speed process, while the renewal and route decoupling assumptions make the bulk service analysis tractable. The stability conditions imply that the system is more likely to be unstable under high incident rates, long incident durations, high demand, low service frequency, and low vehicle capacity. We implement the model on an example bus route adapted from the literature, and the sensitivity analysis shows that congested stations are more vulnerable to random service suspensions. Comparison with a recursive FIFO simulation across seven parameter settings quantifies the approximation error and shows that the discrepancy increases when incident durations move beyond the intended short suspension setting.

The proposed model can support public transit planning, operations, and management in several ways. For the design and planning of resilient systems under random system interruptions, the model can inform headway design and vehicle capacity determination. Its queue length estimates can be used to evaluate the layout and capacity of congested stations. Historical demand and incident information can be used with the model to monitor system resilience and identify critical stations. The estimates of waiting time and queue length can also support cost benefit analyses of service improvements. For example, the model can identify whether a larger vehicle, shorter headway, or higher maintenance frequency is the most cost effective way to keep waiting time below a threshold.

Several extensions are possible. A network level model would differ primarily through the inclusion of transfer passengers. Transfer demand could enter the arrival process, but the model would require additional parameters related to transfers and alighting rates. The renewal approximation could also be replaced by a dependence aware model that carries recursive FIFO delay and the joint passenger queue and vehicle load state. Another extension would consider the partial interruptions studied in Section 2.2. Vehicles retain a positive speed during a partial interruption, which requires a revised headway distribution. A final extension would relax the assumption of no passenger balking or reneging and incorporate more complicated passenger behavior.

7 Modeling assumptions↩︎

Table 3 summarizes route scope and stochastic primitives, and Table 4 summarizes dependence and service sequence approximations.

8 Certified interpolation root search↩︎

Algorithm 13 gives the complete procedure for the three search steps described in Section 4.3.6 and reports the parameter settings used for every station and numerical scenario.

9 Proofs↩︎

9.1 Proof of Proposition [prop95space95dist]↩︎

Suppose that vehicle \(l\) carries \(i\) passengers when it reaches station \(n\). During alighting, each passenger independently remains onboard with probability \(1-\alpha^{(n)}\). The remaining load therefore follows a binomial distribution with parameters \(i\) and \(1-\alpha^{(n)}\). For every \(j=0,\ldots,i\), the corresponding conditional probability is \(a_{ij}^{(n)}=\binom{i}{j}(1-\alpha^{(n)})^j(\alpha^{(n)})^{i-j}\), as stated in Eq. 12 .

The law of total probability gives \(g_j^{(n,l)}=\sum_{i=0}^{C}v_i^{(n-1,l)}a_{ij}^{(n)}\). In vector form, this relation is \[\begin{align} g^{(n,l)} = v^{(n-1,l)}A^{(n)}, \label{eq95g1} \end{align}\tag{27}\] where \(v^{(n-1,l)} = [v^{(n-1,l)}_0,\ldots,v^{(n-1,l)}_C]\) and \(g^{(n,l)} = [g^{(n,l)}_0,\ldots,g^{(n,l)}_C]\). Equation 5 shows that \(S^{(n,l)}=C-G^{(n,l)}\). Hence, the distribution of the available spaces after alighting satisfies \[\begin{align} s_k^{(n,l)} = g_{C-k}^{(n,l)} \quad\quad \forall k = 0,1,\ldots,C. \label{eq95s1} \end{align}\tag{28}\]

Equations 27 and 28 hold for every vehicle run. When the steady state limits defined in Section 4.1 exist, taking \(l\to\infty\) converts these run specific relations into their steady state counterparts and thereby establishes Proposition [prop95space95dist].

9.2 Proof of Proposition [prop95load95dist]↩︎

After alighting and boarding at station \(n\), the steady state vehicle load satisfies \(V^{(n)}=\min\{C,G^{(n)}+Q^{(n)}\}\). Define \(b_{ij}^{(n)}:=\mathbb{P}(V^{(n)}=j\mid G^{(n)}=i)\). If \(0\leq i\leq j<C\), then \(b_{ij}^{(n)}=\mathbb{P}(Q^{(n)}=j-i\mid G^{(n)}=i)\). By Assumption 2, this conditional probability equals the marginal probability \(q_{j-i}^{(n)}\), which gives the first case in Eq. 14 .

If \(j=C\) and \(0\leq i<C\), the vehicle reaches capacity when \(Q^{(n)}\geq C-i\). Assumption 2 therefore gives \(b_{iC}^{(n)}=\mathbb{P}(Q^{(n)}\geq C-i)=1-\sum_{k=0}^{C-i-1}q_k^{(n)}\). If \(i=j=C\), the vehicle is already full after alighting, so \(b_{CC}^{(n)}=1\). All other transitions are infeasible and have probability zero.

The law of total probability now gives \(v_j^{(n)}=\sum_{i=0}^{C}g_i^{(n)}b_{ij}^{(n)}\) for every \(j=0,\ldots,C\). Writing these equations in matrix form yields \(v^{(n)}=g^{(n)}B^{(n)}\), which completes the proof.

9.3 Proof of Proposition [prop95pgf95of95q]↩︎

The proof adapts the transform argument of [2] to the arbitrary available capacity distribution \(s_0^{(n)},\ldots,s_C^{(n)}\), consistent with the variable capacity transform of [9].

We begin with the one step transition of the queue. From the relationship in Eq. 3 , \[\begin{align} q^{(n,l+1)}_k = \sum_{i=0}^{k} r^{(n,l)}_i y_{k-i}^{(n,l)}, \label{eq95queueing95prob1} \end{align}\tag{29}\] where \(r_k^{(n,l)} := \mathbb{P}(R^{(n,l)}=k)\) and \(y_k^{(n,l)} := \mathbb{P}(Y^{(n,l)}=k)\) for every nonnegative integer \(k\). For \(u=0,\ldots,C\), also write \(\left.r^{(n,l)}_k\right\vert_{S^{(n,l)}=u} := \mathbb{P}(R^{(n,l)}=k\mid S^{(n,l)}=u)\). Assumption 3 makes \(S^{(n,l)}\) independent of \(Q^{(n,l)}\). Conditional on \(S^{(n,l)}=u\), all passengers can board when \(u\geq Q^{(n,l)}\), so no passenger is left behind. Therefore, \[\begin{align} \left.r^{(n,l)}_0\right\vert_{S^{(n,l)}=u} = \mathbb{P}(u \geq Q^{(n,l)}\mid S^{(n,l)}=u) = \mathbb{P}(u \geq Q^{(n,l)}) = \sum_{k = 0}^u q^{(n,l)}_k. \label{eq95relation95rq1} \end{align}\tag{30}\] If \(u<Q^{(n,l)}\), only \(u\) passengers can board and \(Q^{(n,l)}-u\) passengers remain. Hence,

\[\begin{align} \left.r^{(n,l)}_i\right\vert_{S^{(n,l)}=u} = \mathbb{P}(Q^{(n,l)} - u = i\mid S^{(n,l)}=u) = \mathbb{P}(Q^{(n,l)} - u = i) = q^{(n,l)}_{u+i}, \qquad i=1,2,\ldots. \label{eq95relation95rq2} \end{align}\tag{31}\]

Using Eqs. 30 and 31 , Eq. 29 becomes \[\begin{align} q^{(n,l+1)}_k = \sum_{u=0}^C s^{(n,l)}_u \left(\sum_{i=0}^{u} q^{(n,l)}_i y_{k}^{(n,l)} + \sum_{i=u+1}^{u+k} q^{(n,l)}_i y_{k-i+u}^{(n,l)}\right). \label{eq95queueing95prob2} \end{align}\tag{32}\]

Assuming that the steady state probabilities exist, we have \(\lim_{l\to\infty}q_k^{(n,l)}=q_k^{(n)}\), \(\lim_{l\to\infty}s_k^{(n,l)}=s_k^{(n)}\), and \(\lim_{l\to\infty}y_k^{(n,l)}=y_k^{(n)}\). Taking limits on both sides of Eq. 32 gives \[\begin{align} q^{(n)}_k = \sum_{u=0}^C s^{(n)}_u\sum_{i=0}^{u} q^{(n)}_i y_{k}^{(n)} + \sum_{u=0}^C s^{(n)}_u\sum_{i=u+1}^{u+k} q^{(n)}_i y_{k-i+u}^{(n)}. \end{align}\]

Let \(Q(z)\), \(R(z)\), and \(Y(z)\) denote the probability generating functions (PGFs) of the steady state random variables \(Q^{(n)}\), \(R^{(n)}\), and \(Y^{(n)}\), respectively. Suppressing the station superscript, we write \[\begin{align} Q(z) &= \sum_{k=0}^{\infty} q^{(n)}_k z^k,\\ R(z) &= \sum_{k=0}^{\infty} r^{(n)}_k z^k = \sum_{k=0}^{\infty} z^k \sum_{u=0}^C s^{(n)}_u\cdot \left.r^{(n)}_k\right\vert_{S^{(n)}=u}, \label{eq95PGF95R}\\ Y(z) &= \sum_{k=0}^{\infty} y^{(n)}_k z^k. \end{align}\tag{33}\] Substituting Eqs. 30 and 31 into Eq. 33 gives \[\begin{align} R(z) &= \sum_{u=0}^C s^{(n)}_u \sum_{i=0}^{u}q_i^{(n)} + \sum_{k=1}^{\infty} z^k \sum_{u=0}^C s^{(n)}_u q_{k+u}^{(n)}\\ & = \sum_{u=0}^C s^{(n)}_u\left[\sum_{i=0}^{u}q_i^{(n)} + \frac{1}{z^u}Q(z) - \frac{1}{z^u}\sum_{i=0}^u q_i^{(n)} z^i \right]. \label{eq95R95intermsof95Q} \end{align}\tag{34}\]

Taking the steady state limit in Eq. 3 gives \(Q^{(n)}=R^{(n)}+Y^{(n)}\). For a fixed vehicle run, the residual queue \(R^{(n)}\) is determined by the queue and service outcomes through the current departure. Under Assumption 1, \(Y^{(n)}\) is generated by the exogenous Poisson process during the next independently drawn headway and is therefore independent of this history. Consequently, \(R^{(n)}\) and \(Y^{(n)}\) are independent in the analytical model, and \[\begin{align} Q(z) = R(z)Y(z). \label{eq95Q95eq95RY} \end{align}\tag{35}\] Combining Eqs. 34 and 35 yields \[\begin{align} Q(z) &= \frac{Y(z)\sum_{u=0}^{C}s^{(n)}_u\left[\sum_{i=0}^u q_i^{(n)} (1-\frac{z^i}{z^u}) \right]}{1 - \sum_{u=0}^{C}s_u^{(n)} \frac{Y(z)}{z^u}} \nonumber\\ &= \frac{\sum_{u=0}^{C}s^{(n)}_u\left[\sum_{i=0}^u q_i^{(n)} (z^C-z^{C-u+i}) \right]}{\frac{z^C}{Y(z)} - \sum_{u=0}^{C}s_u^{(n)} z^{C-u}}. \end{align}\]

9.4 Proof of Proposition [prop95solve95q]↩︎

Section 4.3.1 identifies the \(C\) characteristic roots of the denominator but does not determine the queue probabilities. We now use these roots to factor the numerator of \(Q(z)\) and then match its polynomial coefficients to obtain \(q_0^{(n)},\ldots,q_{C-1}^{(n)}\).

The numerator of \(Q(z)\) is a polynomial of degree at most \(C\). In Eq. 15 , the term with \(i=u\) vanishes, so the coefficient of \(z^C\) in this numerator is \(\kappa^{(n)}:=\sum_{u=1}^{C}s_u^{(n)}\sum_{i=0}^{u-1}q_i^{(n)}\). Analyticity requires the numerator to cancel all \(C\) characteristic roots according to multiplicity. Therefore, the root \(z_i^*\) is repeated in the following product whenever its multiplicity exceeds one, and \(Q(z)\) can be written as \[\begin{align} Q(z) &= \frac{\kappa^{(n)}(z-1)\prod_{i=1}^{C-1}(z-z_i^*)}{\frac{z^C}{Y(z)} - \sum_{u=0}^{C}s_u^{(n)} z^{C-u}}. \label{eq95Qz2} \end{align}\tag{36}\]

As \(z\to1\), both \(\mathrm{\small Num}(z)\) and \(\mathrm{\small Den}(z)\) approach zero, while \(Q(z)\to1\). Applying l’Hôpital’s rule to this limit gives the following normalization identity. \[\begin{align} &\lim_{z\to1} Q(z) = 1 = \lim_{z\to1} \frac{\mathrm{\small Num}'(z)}{\mathrm{\small Den}'(z)} = \frac{\kappa^{(n)}\prod_{i=1}^{C-1}(1-z_i^*)}{\sum_{u=0}^{C}s_u^{(n)}u - Y'(1) } \\ \Rightarrow & \kappa^{(n)} = \frac{\sum_{u=0}^{C}s_u^{(n)}u - Y'(1) }{\prod_{i=1}^{C-1}(1-z_i^*)}. \label{eq95sq} \end{align}\tag{37}\]

Recall that \(\overline{Y}^{(n)} := Y'(1)=\mathbb{E}[Y^{(n)}]\) is the mean number of passenger arrivals within a headway and \(\overline{S}^{(n)} := \sum_{u=0}^{C}s_u^{(n)}u=\mathbb{E}[S^{(n)}]\) is the mean number of available spaces in an arriving vehicle. Substituting Eq. 37 into Eq. 36 gives \[\begin{align} Q(z) &= \frac{(\overline{S}^{(n)} - \overline{Y}^{(n)}) (z-1)\prod_{i=1}^{C-1}\frac{z-z_i^*}{1-z_i^*}}{\frac{z^C}{Y(z)} - \sum_{u=0}^{C}s_u^{(n)} z^{C-u}}. \label{eq95Qz3} \end{align}\tag{38}\]

Equating the numerators in Eqs. 38 and 15 gives \[\begin{align} (\overline{S}^{(n)} - \overline{Y}^{(n)}) (z-1)\prod_{i=1}^{C-1}\frac{z-z_i^*}{1-z_i^*} = \sum_{u=0}^{C}s^{(n)}_u\left[\sum_{i=0}^u q_i^{(n)} (z^C-z^{C-u+i}) \right]. \label{eq95Qz95num} \end{align}\tag{39}\] Both sides of Eq. 39 are polynomials in \(z\), so their coefficients must be equal. Matching the constant terms on both sides gives the first coefficient relation. \[\begin{align} (-1)(\overline{S}^{(n)} - \overline{Y}^{(n)})\prod_{i=1}^{C-1}\frac{z_i^*}{z_i^*-1} = -q_0^{(n)}s_C^{(n)}. \end{align}\] Solving this coefficient relation for the empty queue probability gives \[\begin{align} q_0^{(n)} = \frac{1}{s_C^{(n)}} (\overline{S}^{(n)} - \overline{Y}^{(n)})\prod_{i=1}^{C-1}\frac{z_i^*}{z_i^*-1}. \end{align}\]

The coefficients of the characteristic equation are real, so every nonreal root and its multiplicity are paired with the corresponding conjugate root. The product in Eq. 17 is therefore real. Under stability and irreducibility, the stationary empty queue probability is positive, which establishes the strict positivity of the product.

As a consistency check, consider fixed capacity with \(s_C^{(n)}=1\) and \(\overline{S}^{(n)}=C\). Under this specialization, Eq. 39 reduces to the following fixed capacity expression. \[\begin{align} \left.q_0^{(n)}\right\vert_{s_C^{(n)} = 1} = (C - \overline{Y}^{(n)})\prod_{i=1}^{C-1}\frac{z_i^*}{z_i^*-1}. \end{align}\] This expression agrees with [22]. We next derive \(q_{1:C-1}^{(n)}\). The numerator of Eq. 38 can be rewritten as \[\begin{align} (\overline{S}^{(n)} - \overline{Y}^{(n)}) (z-1)\prod_{i=1}^{C-1}\frac{z-z_i^*}{1-z_i^*} &= \frac{1}{s_C^{(n)}} (\overline{S}^{(n)} - \overline{Y}^{(n)}) \prod_{i=1}^{C-1}\frac{z_i^*}{z_i^*-1} \prod_{i=1}^{C-1}\frac{z_i^* - z}{z_i^*} (z-1)s_C^{(n)} \nonumber \\ &= s_C^{(n)} q_0^{(n)} (z-1) \prod_{i=1}^{C-1} \left(1 - \frac{z}{z_i^*}\right)\nonumber \\ &= - s_C^{(n)} q_0^{(n)} \prod_{i=0}^{C-1} \left(1 - \frac{z}{z_i^*}\right). \label{eq95Qz95num952} \end{align}\tag{40}\] Define \(\prod_{i=0}^{C-1}(1-z/z_i^*) := \sum_{j=0}^{C}\eta_j^{(n)}z^j\), where \(\eta_j^{(n)}\) is the coefficient of \(z^j\) for station \(n\). On the right hand side of Eq. 39 , the coefficient of \(z^{C-k}\) is \(-\sum_{u=k}^Cs_u^{(n)}q_{u-k}^{(n)}\). In Eq. 40 , the corresponding coefficient is \(-s_C^{(n)}q_0^{(n)}\eta_{C-k}^{(n)}\). Matching these coefficients gives \[\begin{align} s_C^{(n)} q_0^{(n)} \eta_{C-k}^{(n)} = \sum_{u=k}^C s_u^{(n)} q_{u-k}^{(n)}, \qquad k=1,2,\ldots,C-1. \label{eq95q95sq} \end{align}\tag{41}\] For fixed capacity, \(s_C^{(n)}=1\) and \(s_k^{(n)}=0\) for \(0\leq k<C\). Equation 41 then reduces to

\[\begin{align} q_{C-k}^{(n)} = q_0^{(n)} \eta_{C-k}^{(n)}, \qquad k=1,2,\ldots,C-1, \quad \text{if }s_C^{(n)}=1, \end{align}\] which agrees with the result of [22].

Adding the identity \(s_C^{(n)}q_0^{(n)}\eta_0^{(n)}=s_C^{(n)}q_0^{(n)}\), where \(\eta_0^{(n)}=1\), converts Eq. 41 into the following matrix representation used to solve the queue probabilities. \[\begin{align} \widetilde{\eta}^{(n)} = q_{0:C-1}^{(n)}\Lambda^{(n)}, \end{align}\] where \(\widetilde{\eta}^{(n)} := [s_C^{(n)}q_0^{(n)}\eta_0^{(n)},s_C^{(n)}q_0^{(n)}\eta_1^{(n)},\ldots,s_C^{(n)}q_0^{(n)}\eta_{C-1}^{(n)}]\in\mathbb{R}^{C}\) and

\[\Lambda^{(n)} := \left[ \begin{array}{ccccc} s_C^{(n)} & s_{C-1}^{(n)} & \cdots & s_2^{(n)} & s_1^{(n)} \\ 0 & s_C^{(n)} & \cdots & s_3^{(n)} & s_2^{(n)} \\ \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & \cdots & 0 & s_C^{(n)} & s_{C-1}^{(n)} \\ 0 & \cdots & \cdots & 0 & s_C^{(n)} \end{array}\right] \in \mathbb{R}^{C\times C}.\] Because \(s_{C}^{(n)} > 0\), the triangular matrix \(\Lambda^{(n)}\) is invertible. Therefore, \[\begin{align} q_{0:C-1}^{(n)} = \widetilde{\eta}^{(n)} (\Lambda^{(n)})^{-1}. \end{align}\]

9.5 Proof of Proposition [prop95q95length]↩︎

The derivation adapts the transform argument of [2] and is consistent with the general bulk service result of [9]. For \(X\in\{S,Y\}\), the symbols \(\overline{X}^{(n)}\), \(\overline{\overline{X}}^{(n)}\), and \(\overline{\overline{\overline{X}}}^{(n)}\) denote its mean, second central moment, and third central moment, respectively. Define \[\begin{align} \mathcal{A}(z) &:= \frac{(\overline{S}^{(n)}-\overline{Y}^{(n)})(z-1)}{\frac{z^C}{Y(z)}-\sum_{u=0}^{C}s_u^{(n)}z^{C-u}}, & \mathcal{B}_i(z) &:= \frac{z-z_i^*}{1-z_i^*}. \label{eq95AB95definition} \end{align}\tag{42}\] Equation 38 gives \(Q(z)=\mathcal{A}(z)\prod_{i=1}^{C-1}\mathcal{B}_i(z)\). At \(z=1\), \(Q(1)=\mathcal{A}(1)=\mathcal{B}_i(1)=1\) and \(\mathcal{B}_i'(1)=1/(1-z_i^*)\). Differentiating the product then gives \(Q'(1)=\mathcal{A}'(1)+\sum_{i=1}^{C-1}1/(1-z_i^*)\).

Let \(\mathcal{A}_1(z):=(\overline{S}^{(n)}-\overline{Y}^{(n)})(z-1)\) and let \(\mathcal{A}_2(z)\) denote the denominator of \(\mathcal{A}(z)\) in Eq. 42 . We have \(\mathcal{A}_1(1)=\mathcal{A}_2(1)=0\) and \(\mathcal{A}_1'(1)=\mathcal{A}_2'(1)=\overline{S}^{(n)}-\overline{Y}^{(n)}>0\). A Taylor expansion of \(\mathcal{A}_2(z)\) around \(z=1\) gives the following two derivative identities. \[\begin{align} \mathcal{A}'(1) &= -\frac{\mathcal{A}_2''(1)}{2\mathcal{A}_2'(1)}, & \mathcal{A}''(1) &= \frac{3[\mathcal{A}_2''(1)]^2-2\mathcal{A}_2'(1)\mathcal{A}_2'''(1)}{6[\mathcal{A}_2'(1)]^2}. \label{eq95A95derivatives} \end{align}\tag{43}\] Differentiating \(\mathcal{A}_2(z)\) and expressing the resulting raw moments in terms of central moments gives \[\begin{align} \mathcal{A}'(1)=\frac{\overline{\overline{S}}^{(n)}+\overline{\overline{Y}}^{(n)}+(\overline{S}^{(n)}-\overline{Y}^{(n)})[1+2(\overline{S}^{(n)}-C)]-(\overline{S}^{(n)}-\overline{Y}^{(n)})^2}{2(\overline{S}^{(n)}-\overline{Y}^{(n)})}. \end{align}\] Substitution into \(Q'(1)\) yields the mean in Eq. 18 .

For the variance, differentiating \(\log Q(z)=\log \mathcal{A}(z)+\sum_{i=1}^{C-1}\log \mathcal{B}_i(z)\) twice and then applying \(\mathcal{A}(1)=\mathcal{B}_i(1)=1\) and \(\mathcal{B}_i''(1)=0\) yields the following identity. \[\begin{align} \operatorname{Var}[Q^{(n)}] = \mathcal{A}''(1) - \mathcal{A}'(1)^2 + \mathcal{A}'(1) + \sum_{i=1}^{C-1} \left(\mathcal{B}_i'(1) - [\mathcal{B}_i'(1)]^2 \right). \label{eq95var95Q2} \end{align}\tag{44}\]

Using Eq. 43 , differentiating \(\mathcal{A}_2(z)\) through third order, and then converting the resulting raw moments to central moments yields the following expression. \[\begin{align} &\mathcal{A}''(1)-\mathcal{A}'(1)^2+\mathcal{A}'(1) \nonumber\\ &=\frac{1}{12(\overline{S}^{(n)}-\overline{Y}^{(n)})^2}\bigg[ -4(\overline{\overline{\overline{S}}}^{(n)}-\overline{\overline{\overline{Y}}}^{(n)})(\overline{S}^{(n)}-\overline{Y}^{(n)}) +3(\overline{\overline{S}}^{(n)}+\overline{\overline{Y}}^{(n)})^2 \nonumber\\ &\qquad-[6(\overline{\overline{S}}^{(n)}-\overline{\overline{Y}}^{(n)})-1](\overline{S}^{(n)}-\overline{Y}^{(n)})^2 -(\overline{S}^{(n)}-\overline{Y}^{(n)})^4\bigg]. \label{eq95Q95var95first95part} \end{align}\tag{45}\] Finally, \(\mathcal{B}_i'(1)-[\mathcal{B}_i'(1)]^2=-z_i^*/(1-z_i^*)^2\). Substituting this identity and Eq. 45 into Eq. 44 yields the variance in Proposition [prop95q95length].

9.6 Proof of Proposition [prop95hdw95distribution]↩︎

Without incidents, vehicle \(l\) requires \(T^{(n)}\) time units to reach station \(n\). Incident onsets occur only during normal travel and follow a Poisson process with rate \(\gamma\), so the number of incidents during this travel follows a Poisson distribution with mean \(\gamma T^{(n)}\). Independent exponential incident durations with rate \(\theta\) then yield the compound Poisson representation in Eq. 21 .

9.7 Proof of Proposition [prop95headway95distribution1]↩︎

The raw headway moment generating function (MGF) is \[\begin{align} M_{\widehat{H}^{(n,l)}}(t) &= \mathbb{E}[e^{t \widehat{H}^{(n,l)}}] = \mathbb{E}\left[e^{t(\overline{H} + \frac{2 \gamma T^{(N)}}{\theta\overline{F}} + I^{(n,l)} - I^{(n,l-1)})}\right] \notag\\ & = e^{t(\overline{H} + \frac{2 \gamma T^{(N)}}{\theta\overline{F}})} \mathbb{E}[e^{tI^{(n,l)}}] \mathbb{E}[e^{-tI^{(n,l-1)}}] \notag\\ & = e^{t(\overline{H} + \frac{2 \gamma T^{(N)}}{\theta\overline{F}})} e^{\gamma T^{(n)} (\frac{\theta}{\theta - t} - 1)} e^{\gamma T^{(n)} (\frac{\theta}{\theta + t} - 1)} \notag\\ & = e^{t(\overline{H} + \frac{2 \gamma T^{(N)}}{\theta\overline{F}})} e^{\gamma T^{(n)} (\frac{2t^2}{\theta^2 - t^2})}, \qquad |t|<\theta. \label{eq95H95mgf95ind} \end{align}\tag{46}\] The factorization in Eq. 46 follows from the independence of the two copies \(I^{(n,l)}\) and \(I^{(n,l-1)}\) used to define one generic raw headway. Because this is the common marginal law used for every renewal draw, \(M_{\widehat{H}^{(n)}}(t)=M_{\widehat{H}^{(n,l)}}(t)\). Assumption 1 separately imposes independence across successive headway draws in the analytical renewal sequence.

9.8 Proof of Proposition [prop95headway95distribution]↩︎

Let \(\mu = \overline{H} + \frac{2 \gamma T^{(N)}}{\theta\overline{F}}\) and \(\sigma = \frac{2 \sqrt{T^{(n)} \gamma}}{\theta}\) denote the mean and standard deviation of \(\widehat{H}^{(n)}_{\text{Normal}}\), respectively. Because \(H^{(n)}_{\text{Normal}}\) is the positive part of this variable, its MGF is \[\begin{align} M_{H^{(n)}_{\text{Normal}}}(t) &= \mathbb{E}[e^{t H^{(n)}_{\text{Normal}}}] = \mathbb{P}[\widehat{H}^{(n)}_{\text{Normal}}\leq0]\cdot e^{0} + \int_{0}^{+\infty} e^{t z} \cdot \phi_{\widehat{H}^{(n)}_{\text{Normal}}}{(z)} \cdot dz \notag\\ &= \Phi(\frac{-\mu}{\sigma}) + \frac{1}{\sigma\sqrt{2\pi}} \int_{0}^{+\infty} e^{t z + \frac{(z - \mu)^2}{-2\sigma^2}} dz \notag\\ & = \Phi(\frac{-\mu}{\sigma}) + e^{\mu t + \frac{\sigma^2t^2}{2}}\left[1 - \Phi(\frac{-\mu}{\sigma} - \sigma t)\right]. \label{eq95positive95normal95mgf} \end{align}\tag{47}\] Completing the square in the exponent gives the final equality in Eq. 47 . Substituting the definitions of \(\mu\) and \(\sigma\) into this expression establishes Proposition [prop95headway95distribution].

9.9 Proof of Proposition [prop95headway95increase]↩︎

For a normal variable with mean \(\mu\) and standard deviation \(\sigma>0\), Eq. 24 gives the expectation of its positive part. Differentiating this expectation with respect to its two parameters gives \[\begin{align} \frac{\partial\mathbb{E}[H^{(n)}_{\text{Normal}}]}{\partial\mu} &=\Phi\left(\frac{\mu}{\sigma}\right)>0, & \frac{\partial\mathbb{E}[H^{(n)}_{\text{Normal}}]}{\partial\sigma} &=\phi\left(\frac{\mu}{\sigma}\right)>0. \label{eq95positive95normal95derivatives} \end{align}\tag{48}\] The mean \(\mu=\overline{H}+2\gamma T^{(N)}/(\theta\overline{F})\) and standard deviation \(\sigma=2\sqrt{T^{(n)}\gamma}/\theta\) are nondecreasing in both \(\gamma\) and \(1/\theta\). The mean is strictly increasing in \(\gamma\), and both parameters are strictly increasing in \(1/\theta\) when \(\gamma>0\). Equation 48 establishes the stated monotonicity. The case \(\gamma=0\) follows by continuity. Figure 14 illustrates these two effects.

9.10 Proof of Proposition [prop95arrival95pax]↩︎

Conditional on \(H^{(n)}_{\text{Normal}}\), the variable \(Y^{(n)}\) is Poisson with parameter \(\lambda^{(n)}H^{(n)}_{\text{Normal}}\). The law of iterated expectations therefore gives \(Y(z)=\mathbb{E}\{\exp[\lambda^{(n)}H^{(n)}_{\text{Normal}}(z-1)]\}=M_{H^{(n)}_{\text{Normal}}}[\lambda^{(n)}(z-1)]\). For \(\gamma>0\), substituting \(t=\lambda^{(n)}(z-1)\) into the MGF in Proposition [prop95headway95distribution] gives the stated PGF. When \(\gamma=0\), \(H^{(n)}_{\text{Normal}}=\overline{H}\) almost surely, so the conditional Poisson PGF directly gives \(Y(z)=\exp[\lambda^{(n)}\overline{H}(z-1)]\).

9.11 Proof of Proposition [prop95stability]↩︎

We establish stability directly from the transient queue dynamics. By the definitions of the queue and available capacity, \(R^{(n,l)}=\max\{0,Q^{(n,l)}-S^{(n,l)}\}\). Combining this relation with Eq. 3 gives \[\begin{align} R^{(n,l+1)}=\max\{0,R^{(n,l)}+Y^{(n,l)}-S^{(n,l+1)}\}. \label{eq95lindley95residual} \end{align}\tag{49}\] Following [19], Eq. 49 represents the residual passenger queue after successive vehicle service events as the standard Lindley recursion.

Assumption 1, the independent increments of the Poisson passenger arrival process, and Assumption 3 imply that the net inputs \(Y^{(n,l)}-S^{(n,l+1)}\) are identically and independently distributed across vehicle runs. They are exogenous to the current residual queue and have finite mean \(\overline{Y}^{(n)}-\overline{S}^{(n)}\). The stability theorem of [20] therefore applies directly to Eq. 49 . If \(\overline{Y}^{(n)}<\overline{S}^{(n)}\), the recursion has a proper stationary solution. Because \(Q^{(n,l+1)}=R^{(n,l)}+Y^{(n,l)}\), the queue observed at vehicle arrival also has a proper stationary distribution.

If \(\overline{Y}^{(n)}>\overline{S}^{(n)}\), the strong law of large numbers gives a positive long run drift for the partial sums of the net inputs, and the residual queue grows without bound. At the boundary \(\overline{Y}^{(n)}=\overline{S}^{(n)}\), the nonconstant independent net inputs in Assumption 3 define a critical Lindley recursion. The critical case has no proper stationary probability distribution under the same result of [20]. Hence, the station is stable if and only if \(\overline{Y}^{(n)}<\overline{S}^{(n)}\). Dividing both sides by \(\overline{S}^{(n)}>0\) and substituting the expressions for the two means gives Eq. 26 .

Upstream marginal propagation supplies \(s^{(n)}\), and the proof requires only that the resulting station service sequence satisfy Assumption 3. The assumption makes this sequence exogenous to the local queue history. The station by station argument therefore does not presuppose the existence of \(q_0^{(n)}\) or another stationary queue probability at the station under consideration.

10 Simulation procedure for comparison↩︎

The simulation deliberately relaxes the analytical renewal, marginal route decoupling, and exogenous station service sequence assumptions. For each vehicle, we generate independent compound Poisson exponential incident durations on successive route segments and cumulatively sum them to obtain \(I^{(n,l)}\) at each station. Realized departure times follow the recursive FIFO rule in Eq. 8 . Vehicle \(l\) cannot depart station \(n\) before vehicle \(l-1\), so an extreme delay can propagate through several successors. Passenger queues and vehicle loads are updated jointly rather than through independent marginal distributions. When a vehicle arrives, passengers board according to the first come, first served principle up to capacity \(C\). Queue lengths at vehicle arrival and waiting times of boarded passengers are recorded during the simulation.

We evaluate the reference setting and six variations that cover low capacity, long headway, high demand, no incidents, frequent incidents, and long incident duration. Each setting uses 20 independent random seeds. A replication contains 10,000 vehicle runs, and the first 1,000 runs are excluded to reduce initialization effects. We calculate each performance measure separately within a replication and then use variation across replications to construct a 95% confidence interval. The relative \(\ell_1\) errors compare the mean across replications with the analytical result over stations 1 through 9. Station 10 is excluded from the waiting time evaluation because its zero arrival rate produces no waiting time observations.

References↩︎