Mo et al.

Evaluation of Public Transit Systems

Evaluation of Public Transit Systems under Short Random Service Suspensions: A
Bulk-Service Queuing Approach

Baichuan Mo Department of Civil and Environmental Engineering, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA Li Jin* UMich Joint
Institute, Shanghai Jiao Tong University, Shanghai 200240, China;
Tandon School of Engineering, New York University, Brooklyn, NY 11201, USA
Haris N. Koutsopoulos Department of Civil and Environmental Engineering,
Northeastern University, Boston, MA 02115, USA Zuo-Jun Max Shen Department of
Industrial Engineering and Operations Research, University of California,
Berkeley, Berkeley, CA 94720, USA Jinhua Zhao Department of Urban Studies and
Planning, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

This paper proposes a stochastic framework to evaluate the performance of
public transit systems under short random service suspensions. We aim to
derive closed-form formulations of the mean and variance of the queue length
and waiting time. A bulk-service queue model is adopted to formulate the
queuing behavior in the system. The random service suspension is modeled as a
two-state (disruption and normal) Markov process. We prove that headway is
distributed as the difference between two compound Poisson exponential random
variables. The distribution is used to specify the mean and variance of queue
length and waiting time at each station with analytical formulations. The
closed-form stability condition of the system is also derived, implying that
the system is more likely to be unstable with high incident rates and long
incident duration. The proposed model is implemented on a bus network. Results
show that higher incident rates and higher average incident duration will
increase both the mean and variance of queue length and waiting time, which
are consistent with the theoretical analysis. Crowding stations are more
vulnerable to random service suspensions. The theoretical results are
validated with a simulation model, showing consistency between the two
outcomes.

Bulk service queuing; Random disturbances; Public transit

##  1 Introduction

Public transit systems (PTSs) play a crucial role in cities worldwide,
transporting people to jobs, homes, outings, and other activities. However,
PTSs are usually susceptible to unplanned delays and service disruptions,
which happen frequently in PTSs. According to Mo et al. (2022), there are on
average 75 incidents happening in the Chicago urban rail system per day and
more than 75% of them are less than 5 minutes. Causes for these short-term
suspensions can be signal system failures, passenger behavior, and
infrastructure problems. For this reason, it is important to recognize how a
PTS is affected by these short-term service suspensions.

In this study, we consider two key performance metrics of PTSs: queue length
and waiting time. Specifically, we model a PTS as a bulk-service queue and aim
to derive the closed-form formulations for the mean and variance of
passengers’ queue length and waiting time at a station under random service
suspensions. To this end, we derive a stability criterion for the passenger
queues under the influence of suspensions, which quantifies the throughput
loss due to suspensions. We also characterize the steady-state distribution of
passenger queues, which naturally leads to the quantification of the
aforementioned metrics.

Queuing behavior at a public transit (PT) station is usually modeled as a
bulk-service queue model (Powell 1981, Islam et al. 2014, Wang et al. 2014).
Bulk service means that customers are served in groups rather than
individually. At a PT station, with the arrival of vehicles (e.g., buses or
trains), a group of passengers will board (i.e., being served in groups). If
the vehicle capacity is less than the number of customers waiting, some
customers are left behind (Kahraman and Gosavi 2011). Most of the previous
studies on a bulk service model for PTSs focus on stations (Selvi and
Rosenshine 1983, Powell 1985, Wang et al. 2014). Islam et al. (2014) used a
Markov model to extend the station-level analysis to the route level. However,
these studies all considered PTSs under normal operating conditions. The
studies of PTSs under service suspensions using queuing analysis are limited.
Regarding the treatment of service disruptions in bulk-service queue models,
Madan (1989) first considered a single channel bulk service queue subject to
interruptions. They assumed there are two states (work and repair) in the
system and derived the probability generating function (PGF) of queue length
using steady-state equations. Many researchers extended Madan (1989)’s
framework by considering more channels (Singh and Ram 1991), more
heterogeneous states (Madan 1992, Ayyappan and Karpagam 2020), different
service interruption assumptions (Jayaraman et al. 1994), and different repair
policies (Tadj and Choudhury 2009, Tadj et al. 2012). However, all these
studies assumed that the service is offered with a fixed batch size (i.e.,
fixed capacity), which is not valid for PTSs where the available vehicle
capacity for boarding is a random variable depending on the current vehicle
load. Besides, all these studies used steady-state equations to derive
multiple PGFs of queue length under different system states (e.g., work and
repair). Results are usually mathematically tedious and the queue length and
waiting time can only be analyzed with a very small service batch size (e.g.,
Madan (1989) only analyzed the problem with service batch size equal to 1 and
2, for batch size more than 3, the closed-form formulas are hard to derive).
Finally, previous studies usually consider the breakdown of servers. However,
there is no straightforward way to map the “breakdown of servers” to a PTS
with valid real-world assumptions because, in a PTS, the assumption of an
independent server is invalid due to inter-station passenger flows.

To fill the research gaps, we propose a bulk queue service-based framework to
describe the passenger and vehicle dynamics for a PTS and analyze the system
performance under short random service suspensions. The objective of this
study is to derive the stability condition of a PTS and the mean and variance
of passengers’ queue length and waiting time for each station under random
suspensions. This analysis provides important insights into how short-term
service disruptions impact PTSs’ performance. The results are helpful for the
future design of PT’s control and planning strategies.

Two building blocks of this study are the work by Powell (1985) and Islam et
al. (2014). Powell (1985) proposed a bulk service queue model for
transportation terminals (i.e., station-level) with analytical queue length
and waiting time formulations under normal conditions using transform methods
(as opposed to steady-state equations methods) and Islam et al. (2014)
extended the analysis from station-level to route-level. In this study, we
explicitly model the random service suspension in a single-route PTS (in
reality, it represents a bus route or one-directional rail line, which is a
basic element of more complex PT networks). Different from typical service
interruption studies where servers may break down, we assume a vehicle in the
PTS may suffer from random suspensions. A detailed discussion of this
assumption is provided in Section 3.2, where we show how it corresponds to
many real-world situations and can be seen as the first step toward a general
incident representation in PTSs. Under this assumption, we extend Powell
(1985) and Islam et al. (2014)’s work to obtain the mean and variance of
passengers’ queue length and waiting time at each station in the single route
PTS by analyzing the headway distribution under random service suspensions.

The major contribution of this paper is fourfold:

  1. 1.

This is the first study to explore analytically the bulk-service queuing
problem involving short random service suspensions applied to PTSs. We model
the service suspension in PTSs by analyzing vehicles’ speed profiles, which is
a novel and practical way to consider “server breakdown” in PTSs.

  2. 2.

We prove that the headway under random service suspensions can be represented
as the difference between two compound Poisson exponential variables. We
assume there is no vehicle overtaking and approximate the headway distribution
as a zero-inflated truncated normal distribution to obtain a closed-form
moment-generating function. Based on this we derive the PGF and corresponding
moments of the number of arrival passengers within a headway (these are
critical components for the bulk-service queue model). This is a new
analytical contribution to the bulk-service queuing theory.

  3. 3.

We introduce a Markov chain model to capture the inter-station passenger flow
dynamics based on Islam et al. (2014)’s work, which extends the typical bulk-
service queuing analysis from the station level to the route level.

  4. 4.

We propose an interpolation-based roots-solving method to find all complex
roots for this study’s model specification. Roots-solving is an essential step
to obtain the queue length and waiting time for the bulk-service queuing
model.

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

##  2 Literature review

###  2.1 Bulk queue models

In the bulk service queuing literature, customers are served in a batch of
fixed or variable lengths. The service rate may depend on the number of
customers waiting for service. The motivation for this model rises from
addressing problems in manufacturing systems, elevators, transport systems,
etc.

Bailey (1954) originated the study of bulk queues by considering a system with
simple Poisson arrivals at a server that serves, at particular points in time,
all waiting customers up to a fixed capacity $c$. If no customers are waiting,
a zero number of customers are served, implying that the server is never idle.
The queue, denoted by $M/G^{c}/1$, is described using an embedded Markov chain
defined at points of service completions. Immediately following Bailey (1954),
Downton (1955) obtained the waiting time distribution of bulk service queues
by considering random arrivals and random service time distribution. Jaiswal
(1960) confirmed the results in Downton (1955). He derived the waiting time
distribution using the embedded Markov-chain approach.

The general bulk service rule was first introduced by Neuts (1967), where a
server, upon finishing a batch, may remain idle if there are fewer than $m$
customers waiting for service. Thus all departing batches from the queue have
at least $m$ customers, although no more than the service capacity.

Along with and after those milestone studies, papers have appeared which can
be differentiated on the basis of the queuing types (arrival process, service
process, number of servers), objectives (queues, waiting times, busy periods,
etc.), the time domain of the solution (i.e., steady-state or transient), and
the method of solution (transforms or direct numerical methods). Chaudhry and
Templeton (1983) and Sasikala and Indhira (2016) provide a more complete
review of the developments in bulk service queue models.

###  2.2 Random service disruptions

The subject of queuing systems wherein the service channel is subject to
breakdowns is a popular subject that has received a lot of attention in the
past fifty years. For a recent survey of the related literature, readers can
refer to Krishnamoorthy et al. (2014).

However, most of the research on this topic deals with models where the server
serves the customers one at a time. The related literature on bulk service is
limited. Madan (1989) studied a single-channel queueing system with Poisson
arrivals and exponential service in batches of fixed size. The system is
subject to random interruptions with an operating state and a repairing state.
Both the operating times and the repair times of the service channel are
assumed to be exponential. Madan (1992) generalized the model in Madan (1989)
to the case where the repairs are performed in two phases. Singh and Ram
(1991) extended the model in Madan (1989) by considering a system with three
identical channels, with operating and repair times for all three service
channels distributed exponentially. Jayaraman et al. (1994) considered a
single-server queueing system with general bulk service. Arrivals are Poisson
but alternate between two modes according to whether the server is operational
or in the failed state. The duration of the operating and repair periods are
exponential and phase-type distributions, respectively. Tadj and Choudhury
(2009) analyzed a bulk service queueing system with an unreliable server,
Poisson input, and general service and repair times. Tadj et al. (2012)
considered a bulk service queuing system where service is provided to groups
of customers of fixed size. Service consists of two consecutive phases and may
take a vacation following the second phase of service. While providing
service, the server may break down and a delay period precedes the repair
period.

###  2.3 Queuing models in public transit systems

Queuing theory in PTSs is usually conducted at the station level, aiming at
obtaining the mean queue length and waiting time. In the case of regular
services where headways are equal, assuming that a) passengers arrive at stops
according to a Poisson process and b) passengers can be served by the first
arriving vehicle, the mean waiting time of passengers ($\mathbb{E}[W]$) is
given by:

$\displaystyle\mathbb{E}[W]=H/2,$ (1)

where $H$ is the service headway and $W$ is the passenger waiting time. This
is the most widely used queuing assumption in transit studies (Dial 1967,
Clerq 1972, Wirasinghe 1980). However, in the case where service is not
reliable, the assumption of regular service can be problematic. Numerous
models have been proposed to account for the stochastic nature of headways
(Welding 1957, Osuna and Newell 1972). A well-known model proposed by Osuna
and Newell (1972) with Poisson arrival passengers and stochastic headways is

$\displaystyle\mathbb{E}[W]=\frac{1}{2}\cdot\left[\mathbb{E}[H]+\frac{\text{Var}[H]}{\mathbb{E}[H]}\right],$
(2)

where $\mathbb{E}[H]$ and $\text{Var}[H]$ are the expectation and variance of
headways, respectively. In the case of regular services, the variance is zero
and the model reverts to Equation 1.

However, the results in Eq. 1 and 2 do not consider the vehicle capacity
(i.e., they assume all passengers can board the first vehicle). In a congested
PTS, passengers may be left behind due to limited vehicle capacity, leading to
an increase in waiting times (Mo et al. 2020). Bulk service queue models have
been applied in PTSs to capture the effects of capacity constraints. Powell
(1981, 1983, 1985) used a bulk service queue model to calculate the passenger
queue length and waiting times at public transportation terminals. The closed-
form mean and variance for these two quantities are derived using a transform
method. Rapoport et al. (2010) studied bulk service queues with constant or
variable capacity and exogenously determined arrival times (e.g., passenger
arrivals based on smart card data). Wang et al. (2014) proposed a bulk service
and batch arrival queuing model with reneging behavior to estimate passengers’
waiting for public transport services.

All the aforementioned studies consider the queuing analysis at the station
level. The extension of queuing analysis from a station level to a route level
is not a trivial problem. First, the boarding and alighting behavior at
upstream stations affect the available capacity distribution at downstream
stations. Second, headways may be correlated across stations, leading to
different headway distributions for different stations (Marguier 1985, Hickman
2001). To address this problem, Islam et al. (2014, 2015) proposed a Markov
model to combine the Powell (1981) and Hickman (2001)’s approaches and used a
bulk service model to analyze system performance at the route level. However,
a limitation of their research is that the calculation of headway correlation
does not consider the vehicle capacity (though the capacity constraint is
considered in the queuing behavior), resulting in the inconsistency of model
assumptions. Also, they assume that headways follow the Erlang distribution,
which leads to model tractability but is not consistent with empirical
observations (Bellei and Gkoumas 2010).

Our paper can be seen as an extension of Powell (1985)’s and Islam et al.
(2014)’s work to incorporate random service suspensions in a PTS with more
consistent and reasonable assumptions. And we also characterize the headway
distribution under service suspensions.

###  2.4 Service interruptions in public transit systems

Studies on service interruption in PTSs can be categorized into two groups:
impact analysis and operations control. Impact analysis studies have used a
variety of methods to analyze the impact of service disruptions on performance
and level of service. Of these methods, the most common is based on graph
theory, surveys, simulation, and empirical data. Graph theory-based methods
usually derive resilience or vulnerability indicators based on the network
topology (Yin et al. 2016, Zhang et al. 2018, Xu et al. 2015, Berdica 2002).
These methods are effective for understanding high-level network properties
related to incidents. Survey-based methods investigate passenger behavior and
opinions during incidents (Currie and Muir 2017, Murray-Tuite et al. 2014,
Fukasawa et al. 2012, Teng and Liu 2015, Lin et al. 2018). Passengers’
individual-level behavior is analyzed and understood using econometric models.
Simulation-based methods simulate passenger flows on the transit network under
incident scenarios (Balakrishna et al. 2008, Suarez et al. 2005, Hong et al.
2018). The empirical data-based methods use smart card and vehicle location
data to analyze real-world incident impacts (Sun et al. 2016, Tian and Zheng
2018, Mo et al. 2022). These studies can output many metrics of interest such
as vehicle load, travel delays caused by incidents, distribution of the
impact, etc. Studies focusing on operations control under service disruptions
address aspects including shuttle bus design (Jin et al. 2016, Luo et al.
2019), vehicle holding (O’Dell and Wilson 1999), integrating local services
(Jin et al. 2014), and timetable adjustment (Kroon and Huisman 2011).

The analysis presented in this paper belongs to the “impact analysis”
category, which aims to obtain stability conditions for PTSs and the mean and
variance of passengers’ queue length and waiting time of each station under
short random suspensions. None of the previous studies has used the bulk
service model for this type of analysis.

##  3 Model

###  3.1 Single-route public transit system and vehicle movements

Consider a single-route PTS with $N$ stations as shown in Figure 1. Vehicles
are dispatched from a transportation hub (also referred to as station 0) and
travel from station 1 to station $N$. At a specific station $n$, we assume
that passenger arrivals follow a Poisson process with a fixed rate
$\lambda^{(n)}$ during the time period of interest. When a vehicle arrives at
station $n$, each passenger in the vehicle has a probability of $\alpha^{(n)}$
to alight. Thus, the number of alighting passengers at station $n$ follows a
binomial distribution. Poisson arrivals and binomial alighting are two common
assumptions in much of the PT-related literature (Hickman 2001). In this
study, we do not consider reneging behavior of passengers (i.e., passengers
may leave the system if they have waited for too long) since the focus of the
paper is on “short” service suspensions and we assume passengers choose to
wait. Empirical studies (Sun et al. 2016, Rahimi et al. 2019) show that
passengers start to leave the system only when delays are large (e.g., 30
minutes or more). Incorporating balking and reneging is outside the scope of
this paper and can be a future extension of this work.

(a)  Figure 1: Schematic presentation of a single-route public transit system

Let $l=1,2,...$ be a superscript denoting the vehicle run number (or vehicle
ID). Smaller $l$ means vehicles are dispatched at an earlier time. Figure 2
summarizes the vehicle and passenger interactions at station $n$ over time.
Let $t_{A}^{(n,l)}$ be the time that vehicle $l$ arrives at station $n$, and
$t_{D}^{(n,l)}$ the time that vehicle $l$ departs station $n$. $H^{(n,l)}$ is
the headway between the preceding vehicle $l-1$ and vehicle $l$, as they
depart from stop $n$ (i.e., $H^{(n,l)}=t_{D}^{(n,l)}-t_{D}^{(n,l-1)}$). When a
vehicle arrives at station $n$, some of the onboard passengers alight first,
then the queuing passengers start to board. Let $Q^{(n,l)}$ be the number of
queuing passengers when vehicle $l$ arrives at station $n$, $R^{(n,l)}$ the
number of left-behind passengers when vehicle $l$ departs station $n$, and
$Y^{(n,l)}$ the number of passengers arriving between $t_{D}^{(n,l)}$ and
$t_{A}^{(n,l+1)}$. By definition,

$\displaystyle Q^{(n,l+1)}=R^{(n,l)}+Y^{(n,l)}.$ (3)

(a)  Figure 2: Diagram of vehicles and passengers interaction at station $n$
in the time dimension

In this study, we assume that the dwell time (i.e.,
$t_{D}^{(n,l)}-t_{A}^{(n,l)}$) is negligible compared to the vehicle travel
time ($t_{A}^{(n+1,l)}-t_{D}^{(n,l)}$) such that the number of passengers
arriving during the dwell time is zero (same assumption as in Powell (1981)).
Then, given the headway $H^{(n,l+1)}$, $Y^{(n,l)}|H^{(n,l+1)}$ follows a
Poisson distribution with parameter $\lambda^{(n)}H^{(n,l+1)}$:

$\displaystyle
Y^{(n,l)}\;|\;H^{(n,l+1)}\sim\text{{Poi}}(\lambda^{(n)}H^{(n,l+1)}).$ (4)

In other words, $Y^{(n,l)}$ can be seen as the number of arriving passengers
within a headway (i.e., $H^{(n,l+1)}$).

From the vehicle’s perspective, let $S^{(n,l)}$ be the number of available
spaces after passengers alighting from vehicle $l$ at station $n$, $G^{(n,l)}$
the number of remaining passengers on vehicle $l$ after passengers alighting
at station $n$. By definition,

$\displaystyle G^{(n,l)}=C-S^{(n,l)},$ (5)

where $C$ is the capacity of vehicles. Denote $V^{(n,l)}$ as the vehicle load
(i.e., number of onboard passengers) when vehicle $l$ departs station $n$
(i.e., the vehicle load when it arrives at station $n+1$). Then, the number of
alighting passengers from vehicle $l$ at station $n$ given $V^{(n-1,l)}$
follows a binomial distribution:

$\displaystyle\left(V^{(n-1,l)}-G^{(n,l)}\right)\;|\;V^{(n-1,l)}\sim\text{{Bin}}(V^{(n-1,l)},\alpha^{n}).$
(6)

###  3.2 Random service suspensions and vehicle speed profile

Let us assume that there are random service suspensions when a vehicle travels
in the system. Given these disturbances, the speed curve of vehicle $l$ from
station $n$ to $n+1$ can be described by the red line in Figure 3. Every
random incident causes a speed reduction or stop of the vehicle. In reality,
these incidents can be caused by many reasons. For example, in a bus system,
they may be caused by traffic congestion or accidents, drivers’ or passengers’
behavior, vehicle engine issues, etc. In a rail system, the reasons may be
signal failures, infrastructure problems, and drivers’ or passengers’
behavior. The speed curve is a general representation of different incidents,
interruptions, suspensions, or disruptions that impede the vehicle’s movement.

(a)  Figure 3: Schematic speed curve of vehicle $l$ traveling from station $n$
to $n+1$

The actual vehicle speed profile under interruptions can be complicated. To
facilitate mathematical modeling, we assume that the speed of a vehicle under
random interruptions can be approximated by an impulse function (blue line in
Figure 3). The impulse function separates the vehicle trajectory into
traveling and stopping phases, denoted as the normal and disruption states,
respectively. In the normal state, a vehicle travels at a constant speed. Once
an incident happens, the vehicle stops immediately and enters the disruption
state. We assume that, in a sufficiently small time interval, $\Delta$, the
probability of incident occurrence is $\gamma\Delta$. Furthermore, the
duration of an incident follows an exponential distribution with rate $\theta$
(i.e., mean of $\frac{1}{\theta}$). Then, the state of a vehicle is a two-
state Markov process (Figure 4) with the state space of {Normal, Disruption}.

(a)  Figure 4: Transition diagram of vehicle states

For the two-state Markov process, the duration of disruption and normal states
follows the exponential distribution with rates $\theta$ and $\gamma$,
respectively.

Approximating the actual speed curve as an impulse function can be seen as the
first step toward a general incident representation in PTSs. Actually, any
type of incident can be represented as a mixture of different types of normal
and disruption states. The normal and disruption states can be defined with
heterogeneous occurrence probabilities and duration for different categories
of incidents, which results in a more sophisticated speed curve
representation.

###  3.3 Headway under random service suspensions

Under the assumption of an impulse-function speed profile, all vehicles have
the same fixed travel speed under the normal state. Therefore, if there is no
incident in the system, all stations have the same deterministic nominal
headway (denoted as $\bar{H}$). The relationship among $\bar{H}$, route cycle
time $\bar{E}$ (i.e., the time that a vehicle travels from the transportation
hub to the last station and returns to the hub), and fleet size (denoted as
$\bar{F}$) for the route is

$\displaystyle\bar{H}=\frac{\bar{E}}{\bar{F}}$ (7)

With random service suspensions, the route cycle time would increase. There
are two possible responses for the transit agency: 1) To maintain the same
planned headway $\bar{H}$, the agency needs to increase the fleet size (i.e.,
number of vehicles) for the route. 2) With the same fleet size (i.e., limited
resources), the agency would have to increase the planned headway $\bar{H}$.
In this paper, we consider the second scenario because it reflects incidents’
impact on headway and service performance, which is more relevant to this
paper’s topic.

Therefore, we assume that at the route planning stage, transit agencies have
an estimate of the average delay in the cycle time, $\bar{D}$. Let $I^{(n,l)}$
be the total duration of all incidents happening during the vehicle $l$’s
travel time from the transportation hub to station $n$ (a random variable).
Then $\mathbb{E}[I^{(N,l)}]$ is the expected incident duration for a vehicle
traveling from the transportation hub to the last station $N$. Assuming the
road conditions for two directions of the route are the same, then the total
estimated delay for the cycle trip is

$\displaystyle\bar{D}$ $\displaystyle=2\cdot\mathbb{E}[I^{(N,l)}]$ (8)

It is worth noting that some transit agencies may plan the headway by assuming
a larger delay (e.g., not the mean, but the 85% percentile). Hence, we may
also formulate $\bar{D}$ as a general function of $\mathbb{E}[I^{(N,l)}]$. In
this study, we adopt Eq. 8 for simplicity. Then, the incident-adjusted planned
headway (denoted as $\bar{H}^{\text{Adj}}$) is

$\displaystyle\bar{H}^{\text{Adj}}=\frac{\bar{E}+\bar{D}}{\bar{F}}=\bar{H}+\frac{2\cdot\mathbb{E}[I^{(N,l)}]}{\bar{F}}$
(9)

where $\frac{2\cdot\mathbb{E}[I^{(N,l)}]}{\bar{F}}$ is the planned headway
adjustment term due to incidents. Note that we assume $I^{(N,l)}$ are
identically distributed for all $l$. So the incident-adjusted planned headway
is not affected by vehicle ID. In this study, we assume that the single-route
PTS will dispatch vehicles based on the incident-adjusted planned headway
$\bar{H}^{\text{Adj}}$ and that all dispatches are on time.

Let $T^{(n)}$ be the travel time for vehicles from the transportation hub to
station $n$ when there is no incident ( a fixed constant in this study due to
the fixed speed assumption). Without loss of generality, let us assume vehicle
$(l-1)$ departs from the transportation hub at time $0$. Considering random
service suspensions, vehicle $(l-1)$’s departure time from station $n$ is:

$\displaystyle t_{D}^{(n,l-1)}=T^{(n)}+I^{(n,l-1)}$ (10)

Note that the dwell time is ignored as we assumed before.

Given that the incident-adjusted planned headway is
$\bar{H}+\frac{2\cdot\mathbb{E}[I^{(N,l)}]}{\bar{F}}$, the planned departure
time of vehicle $l$ from station $n$ is

$\displaystyle
t_{D}^{(n,l)}=\bar{H}+\frac{2\cdot\mathbb{E}[I^{(N,l)}]}{\bar{F}}+T^{(n)}+I^{(n,l)}$
(11)

Therefore, with random incidents, the actual headway of vehicle $l$ at station
$n$ is

$\displaystyle
H^{(n,l)}=t_{D}^{(n,l)}-t_{D}^{(n,l-1)}=\bar{H}+\frac{2\cdot\mathbb{E}[I^{(N,l)}]}{\bar{F}}+I^{(n,l)}-I^{(n,l-1)}.$
(12)

In this study, we assume $I^{(n,l)}$ and $I^{(n,l-1)}$ are independent. This
assumption facilitates closed-form derivations. In reality, if the incidents
are caused by road congestion or infrastructure issues, it is possible that
the incident durations for two consecutive vehicles passing through the same
route segment are correlated. However, addressing the correlation is not a
trivial problem in the bulk service queue model (Powell 1981) and is beyond
the scope of this paper.

##  4 Analysis

The objective of this study is to derive the stability conditions of a PTS
(Proposition 4.15) and the mean and variance of passengers’ queue length and
waiting (Propositions 4.5 and 4.6) time at each station under random service
suspensions. Figure 5 shows how the distributions of different random
variables (particularly, $S^{(n,l)},V^{(n,l)},Q^{(n,l)}$) are calculated. The
major calculation consists of three parts:

  * •

Given the distribution of $V^{(n-1,l)}$ (i.e., vehicle load when vehicle $l$
departs station $n-1$), calculate the distribution of $S^{(n,l)}$ (i.e., the
number of available space after passengers alighting from vehicle $l$ at
station $n$). The details are shown in Section 4.1

  * •

Given the distribution of $S^{(n,l)}$, calculate the distribution of
$Q^{(n,l)}$ (i.e., the number of queuing passengers at station $n$ when
vehicle $l$ arrives) and the mean and variance of queue length and waiting
time at station $n$. This is discussed in Section 4.3.

  * •

Given the distribution of $S^{(n,l)}$ and $Q^{(n,l)}$, calculate the
distribution of $V^{(n,l)}$, which is discussed in Section 4.2

(a)  Figure 5: Analysis framework

With the three components, we can derive the distribution of
$S^{(n,l)},Q^{(n,l)},V^{(n,l)}$ for all $n=1,...,N$ given the distribution of
$V^{(0,l)}$ (i.e., vehicle load when vehicle $l$ arrives at the first station,
it is always zero by definition). Note that, in this section, we focus on the
steady-state distribution of these variables (i.e., $l\to\infty$).

After obtaining the corresponding distributions, we discuss the stability
conditions in Section 4.4 and summarize the approach in Section 4.5.

###  4.1 Available vehicle space steady-state distribution

In this section, we aim to derive the steady-state distribution of $S^{(n,l)}$
given the steady-state distribution 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 all $k=0,1,...,C$. Assuming that
the steady state probabilities for all variables exist (the stability
condition will be discussed in Section 4.4), we have
$v^{(n)}_{k}:=\lim_{l\to\infty}v^{(n,l)}_{k}=\mathbb{P}(V^{(n)}=k)$,
$s^{(n)}_{k}:=\lim_{l\to\infty}s^{(n,l)}_{k}=\mathbb{P}(S^{(n)}=k)$, and
$g^{(n)}_{k}:=\lim_{l\to\infty}g^{(n,l)}_{k}=\mathbb{P}(G^{(n)}=k)$, where
$V^{(n)}=\lim_{l\to\infty}V^{(n,l)}$, $S^{(n)}=\lim_{l\to\infty}S^{(n,l)}$,
and $G^{(n)}=\lim_{l\to\infty}G^{(n,l)}$.

###### Proposition 4.1

[Distribution of $S^{(n,l)}$ given $\mathbb{P}[V^{(n-1,l)}]$]: $\forall$ n =
1,..,N, given the distribution of $V^{(n-1)}$ (i.e.,
$v^{(n)}:=[v^{(n-1)}_{0},...,v^{(n-1)}_{C}]\in\mathbb{R}^{C+1}$), the
distribution of $S^{(n)}$ (i.e.,
$s^{(n)}:=[s^{(n-1)}_{0},...,s^{(n-1)}_{C}]\in\mathbb{R}^{C+1}$) is given as:

$\displaystyle s_{k}^{(n)}$ $\displaystyle=g_{C-k}^{(n)}\quad\quad\forall
k=0,1,...,C,$ (13)

where $g^{(n)}:=[g^{(n)}_{0},...,g^{(n)}_{C}]\in\mathbb{R}^{C+1}$ and

$\displaystyle g^{(n)}$ $\displaystyle=v^{(n-1)}A^{(n)},$ (14)

$A^{(n)}$ is a $(C+1)\times(C+1)$ matrix with the element in row $i$ and
column $j$ equal to $a_{ij}^{(n)}$, and $a_{ij}^{(n)}$ is defined as

$a_{ij}^{(n)}=\left\\{\begin{aligned} &1,&\quad\text{ if $i=0$ and $j=0$}\\\
&\binom{i}{i-j}(\alpha^{(n)})^{i-j}(1-\alpha^{(n)})^{j},&\quad\text{ if $i\geq
j$ and $i,j\neq 0$}\\\ &0,&\quad\text{
otherwise}\end{aligned}\right.\quad\quad\forall\;i,j=0,1,...,C$ (15)

###  4.2 Vehicle load steady-state distribution

In this section, we derive the steady-state distribution of $V^{(n,l)}$ given
the steady-state distribution of $G^{(n,l)}$ and $Q^{(n,l)}$. Define
$q^{(n)}_{k}:=\lim_{l\to\infty}q^{(n,l)}_{k}=\mathbb{P}(Q^{(n)}=k)$, where
$Q^{(n)}=\lim_{l\to\infty}Q^{(n,l)}$ and
$q^{(n,l)}_{k}=\mathbb{P}(Q^{(n,l)}=k)$. Denote the first $C$ elements of the
steady-steady queue length distribution as $q_{0:C-1}^{(n)}$, where
$q_{0:C-1}^{(n)}=[q_{0}^{(n)},...,q_{C-1}^{(n)}]\in\mathbb{R}^{C}$.

###### Proposition 4.2

[Distribution of $V^{(n,l)}$ given $\mathbb{P}[G^{(n,l)}]$ and
$\mathbb{P}[Q^{(n,l)}]$]: $\forall$ n = 1,..,N, given the distribution of
$G^{(n)}$ (i.e., $g^{(n)}$) and $q_{0:C-1}^{(n)}$, the distribution of
$V^{(n)}$ can be expressed as:

$\displaystyle v^{(n)}=g^{(n)}B^{(n)}$ (16)

where $B^{(n)}$ is a matrix with the element in row $i$ and column $j$ equal
to $b_{ij}^{(n)}$:

$b_{ij}^{(n)}=\left\\{\begin{aligned} &q^{(n)}_{j-i},&\quad\text{ if $0\leq
i\leq j<C$}\\\ &1-\sum_{k=0}^{C-i-1}q^{(n)}_{k},&\quad\text{ if $j=C$ and
$0\leq i<C$}\\\ &1,&\quad\text{ if $i=j=C$}\\\ &0,&\quad\text{
otherwise}\end{aligned}\right.\quad\quad\forall\;i,j=0,1,...,C$ (17)

###  4.3 Queuing analysis at a station

In this section, assuming that we know the distribution of $S^{(n)}$ (i.e.,
$s^{(n)}=[s^{(n)}_{0},...,s^{(n)}_{C}]\in\mathbb{R}^{C+1}$), our goal is to
derive $q_{0:C-1}^{(n,l)}$ and the mean and variance of passenger queue length
and waiting time.

####  4.3.1 Probability generating function of queue length

We start with deriving the probability generating function (PGF) for
$Q^{(n)}$, where $Q^{(n)}=\lim_{l\to\infty}Q^{(n,l)}$.

###### Proposition 4.3

[PGF of $Q^{(n)}$]: $\forall$ n = 1,..,N, given the distribution of $S^{(n)}$
(i.e., $s^{(n)}$), the PGF of $Q^{(n)}$ can be expressed as:

$\displaystyle
Q(z)=\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}},$
(18)

where $Y{(z)}$ is the PGF of $Y^{(n)}$ and
$Y^{(n)}=\lim_{l\to\infty}Y^{(n,l)}$ is the number of arrival passengers at
station $n$ within a headway at the steady state.

In Eq. 18, there are $C$ unknown variables, $q_{0}^{(n)},...,q_{C-1}^{(n)}$.
Note that $q_{C}^{(n)}$ does not appear in $Q(z)$ because when $u=C$ and
$i=C$, we have $q_{C}^{(n)}(z^{C}-z^{C-u+i})\equiv 0$. To quantify $Q(z)$,
Rouche’s theorem is used (Beardon 2019). Let $\textsc{Num}(z)$ and
$\textsc{Den}(z)$ be the numerator and denominator of $Q(z)$ (i.e.,
$Q(z)=\frac{\textsc{Num}(z)}{\textsc{Den}(z)}$). As shown in Powell (1981),
one can prove that $\textsc{Den}(z)$ (i.e.,
$\frac{z^{C}}{Y(z)}-\sum_{u=0}^{C}s_{u}^{(n)}z^{C-u}$) has exactly $C$ complex
roots within (or on) the unit circle on a complex plane using Rouche’s
theorem. Notice that for any $z\in\mathbb{C}$ that satisfies $|z|\leq 1$,
where $\mathbb{C}$ is the set of complex numbers, the generating function
$Q(z)$ must be analytic. Therefore, if $z^{*}$ is the root of
$\textsc{Den}(z)$ (i.e., $\textsc{Den}(z^{*})=0$), it should also be the root
of $\textsc{Num}(z)$ (i.e., $\textsc{Num}(z^{*})=0$) such that $Q(z)$ is
analytic (Rudin 2006). Hence, one can solve for
$q_{0}^{(n)},...,q_{C-1}^{(n)}$ using the following two steps:

  * •

Step 1: Solve $\textsc{Den}(z)=0$ for $C$ different roots
$z^{*}_{0},...,z^{*}_{C-1}\in\mathbb{C}$ that satisfy $|z_{i}^{*}|\leq 1$,
$\forall\;0\leq i\leq C-1$. Note that $z=1$ is always a root of
$\textsc{Den}(z)$. But it does not give information about
$q_{0}^{(n)},...,q_{C-1}^{(n)}$ as $\textsc{Num}(1)=0$ is naturally satisfied.
Hence, we adopt the convention that $z^{*}_{0}=1$.

  * •

Step 2: Combining $\textsc{Num}(z_{i}^{*})=0$ ($\forall\;1\leq i\leq C-1$) and
$Q(1)=1$, solve for $q_{0}^{(n)},...,q_{C-1}^{(n)}$ (there are $C$ system
equations and $C$ unknown variables). Note that when $z\to 1$, both
$\textsc{Num}(z)$ and $\textsc{Den}(z)$ approach 0. Therefore, using
L’Hopital’s rule,

$\displaystyle\lim_{z\to 1}Q(z)=\lim_{z\to
1}\frac{\textsc{Num}^{\prime}(z)}{\textsc{Den}^{\prime}(z)}=\frac{\sum_{u=0}^{C}s_{u}^{(n)}\left[\sum_{i=0}^{u}q_{i}^{(n)}(u-i)\right]}{\bar{S}^{(n)}-\bar{Y}^{(n)}}=1$
(19)

where $\bar{S}^{(n)}=\sum_{u=0}^{C}us^{(n)}_{u}=\mathbb{E}[{S}^{(n)}]$,
$\bar{Y}^{(n)}=Y^{\prime}(1)=\mathbb{E}[{Y}^{(n)}]$. Eq. 19 is the equation
used to solve for $q_{0:C-1}^{(n)}$ (instead of directly using $Q(1)=1$).

####  4.3.2 Queue length distribution

Though $q_{0}^{(n)},...,q_{C-1}^{(n)}$ can be obtained by solving $C$ system
equations as mentioned in Section 4.3.1, we provide a simpler way to calculate
$q_{0:C-1}^{(n)}$, which is known as matching the polynomial coefficients.

###### Proposition 4.4

[Distribution of $Q^{(n)}$ given $\mathbb{P}[S^{(n)}]$]: $\forall$ n = 1,..,N,
given the distribution of $S^{(n)}$ (i.e., $s^{(n)}$), all complex roots of
$\textsc{Den}(z)$ (i.e., $z^{*}_{0},...,z^{*}_{C-1}$), and $\bar{Y}^{(n)}$, if
$s_{C}^{(n)}>0$, then $q_{0:C-1}^{(n)}$ can be solved as:

$\displaystyle
q_{0}^{(n)}=\frac{1}{s_{C}^{(n)}}(\bar{S}^{(n)}-\bar{Y}^{(n)})\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1},$
(20)

and

$\displaystyle q_{0:C-1}^{(n)}=\tilde{\eta}^{(n)}(\Lambda^{(n)})^{-1},$ (21)

where
$\tilde{\eta}^{(n)}=[s_{C}^{(n)}q_{0}^{(n)}\eta_{0}^{(n)},s_{C}^{(n)}q_{0}^{(n)}\eta_{1}^{(n)},...,s_{C}^{(n)}q_{0}^{(n)}\eta_{C-1}^{(n)}]\in\mathbb{R}^{C}$
and

$\Lambda^{(n)}=\left[\begin{array}[]{ccccccc}s_{C}^{(n)}&s_{C-1}^{(n)}&s_{C-2}^{(n)}&&...&&s_{1}^{(n)}\\\
0&s_{C}^{(n)}&s_{C-1}^{(n)}&&...&&s_{2}^{(n)}\\\
...&0&s_{C}^{(n)}&&...&&s_{3}^{(n)}\\\ 0&...&0&&...&&s_{4}^{(n)}\\\
0&0&...&&...&&...\\\ 0&0&0&&...&&s_{C}^{(n)}\\\
\end{array}\right]\in\mathbb{R}^{C\times C}.$ (22)

$\eta_{j}^{(n)}$ is the polynomial coefficient of $z^{j}$ in
$\prod_{i=0}^{C-1}\left(1-\frac{z}{z_{i}^{*}}\right)$ (i.e.,
$\sum_{j=0}^{C}\eta_{j}^{(n)}z^{j}:=\prod_{i=0}^{C-1}\left(1-\frac{z}{z_{i}^{*}}\right)$).
As $z_{i}^{*}$ is specified for station $n$, a superscript $n$ is added to the
coefficients.

Note that assuming $s_{C}^{(n)}>0$ in Proposition 4.4 is not restrictive
because otherwise we can reduce $C$ such that $s_{C}^{(n)}>0$ always holds.

####  4.3.3 Analytical formulation of mean and variance of queue length and
waiting time

After solving for $q_{0}^{(n)},...,q_{C-1}^{(n)}$, $Q(z)$ is determined. The
expectation and variance of the queue length at station $n$ can be written by
definition as:

$\displaystyle\mathbb{E}[Q^{(n)}]$
$\displaystyle=\sum_{k=0}^{\infty}kq_{k}^{(n)}=\left.\frac{dQ(z)}{d{z}}\right|_{z=1}$
(23) $\displaystyle\text{Var}[Q^{(n)}]$
$\displaystyle=\mathbb{E}[(Q^{(n)})^{2}]-\mathbb{E}[Q^{(n)}]^{2}=\left.\frac{d^{2}Q(z)}{d{z^{2}}}\right|_{z=1}+\mathbb{E}[Q^{(n)}]-\mathbb{E}[Q^{(n)}]^{2}.$
(24)

###### Proposition 4.5

[Mean and variance of queue length]: $\forall$ n = 1,..,N, given the
distribution of $S^{(n)}$ and the expression of $Y(z)$, $\mathbb{E}[Q^{(n)}]$
and $\text{Var}[Q^{(n)}]$ can be calculated as:

$\displaystyle\mathbb{E}[Q^{(n)}]$
$\displaystyle=\frac{\bar{\bar{S}}^{(n)}+\bar{\bar{Y}}^{(n)}+(\bar{S}^{(n)}-\bar{Y}^{(n)})[1+2(\bar{S}^{(n)}-C)]-(\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}}{2(\bar{S}^{(n)}-\bar{Y}^{(n)})}+\sum_{i=1}^{C-1}\frac{1}{1-z_{i}^{*}}$
(25) $\displaystyle\text{Var}[Q^{(n)}]$
$\displaystyle=\frac{1}{12(\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}}\bigg{[}-4(\bar{\bar{\bar{S}}}^{(n)}-\bar{\bar{\bar{Y}}}^{(n)})(\bar{S}^{(n)}-\bar{Y}^{(n)})+3(\bar{\bar{S}}^{(n)}+\bar{\bar{Y}}^{(n)})^{2}$
$\displaystyle-[6({\bar{\bar{S}}}^{(n)}-{\bar{\bar{Y}}}^{(n)})-1](\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}-(\bar{S}^{(n)}-\bar{Y}^{(n)})^{4}\bigg{]}-\sum_{i=1}^{C-1}\frac{z_{i}^{*}}{(1-z_{i}^{*})^{2}}$
(26)

where $\bar{\bar{S}}^{(n)}$ and $\bar{\bar{\bar{S}}}^{(n)}$ (resp.
$\bar{\bar{Y}}^{(n)}$ and $\bar{\bar{\bar{Y}}}^{(n)}$) are the second and
third central moments of $S^{(n)}$ (resp. $Y^{(n)}$).

###### Proposition 4.6

[Mean and variance of waiting time]: $\forall$ n = 1,..,N, given the
distribution of $S^{(n)}$ and the expression of $Y(z)$, the mean and variance
of waiting time at station $n$ (denoted as $W^{(n)}$) is given as:

$\displaystyle\mathbb{E}[W^{(n)}]$
$\displaystyle=\frac{\bar{Q}_{t}^{(n)}}{\lambda^{(n)}}$ (27)
$\displaystyle\text{Var}[W^{(n)}]$
$\displaystyle=\frac{\bar{\bar{Q}}_{t}^{(n)}-\bar{Q}_{t}^{(n)}}{(\lambda^{(n)})^{2}}$
(28)

where ${Q}_{t}^{(n)}$ is the queue length at an arbitrary time point (as
opposed to $Q^{(n)}$ which is the queue length at the time of vehicle
arrival). $\bar{Q}_{t}^{(n)}$ and $\bar{\bar{Q}}_{t}^{(n)}$ are defined as

$\displaystyle\bar{Q}_{t}^{(n)}$
$\displaystyle=\mathbb{E}[Q^{(n)}]-\bar{Y}^{(n)}+\frac{1}{2}\left({\bar{\bar{Y}}^{(n)}}/{\bar{Y}^{(n)}}+\bar{Y}^{(n)}-1\right)$
(29) $\displaystyle\bar{\bar{Q}}_{t}^{(n)}$
$\displaystyle=\text{Var}[Q^{(n)}]-\bar{\bar{Y}}^{(n)}+\frac{1}{12(\bar{Y}^{(n)})^{2}}\left[4\bar{Y}^{(n)}\bar{\bar{\bar{Y}}}^{(n)}+6(\bar{Y}^{(n)})^{2}\bar{\bar{Y}}^{(n)}-(\bar{Y}^{(n)})^{2}+(\bar{Y}^{(n)})^{4}-3(\bar{\bar{Y}}^{(n)})^{2}\right]$
(30)

Eq. 27 is the application of Little’s law. Proposition 4.6 is directly
obtained from Powell (1985).

###### Remark 4.7

The formulation of $\mathbb{E}[Q^{(n)}]$, $\text{Var}[Q^{(n)}]$,
$\mathbb{E}[W^{(n)}]$, and $\text{Var}[W^{(n)}]$ in this study are equivalent
to Powell (1985) because in his paper the $M/G^{[S]}/1$ bulk queue model was
considered, where $G^{[S]}$ represents a general (i.e., arbitrary) bulk-
service distribution, which includes the service distribution incorporating
random service suspension considered in this study. However, this does not
lower the contribution of this paper because to implement these equations, the
formulation of $Y(z)$ needs to be specified. In the next section, we show how
random service suspension introduces a new distribution for $Y^{(n)}$, which
has not been considered in the literature.

####  4.3.4 Headway distribution

According to Propositions 4.4 to 4.6, to calculate $q^{(n)}_{0:C-1}$ and the
mean and variance of queue length and waiting time, it is essential to specify
$Y(z)$ (i.e., the PGF of the number of passengers arriving within a headway).
According to Eq. 4, taking $l\to\infty$ gives that $Y^{(n)}|H^{(n)}$ is a
Poisson random variable with parameter $\lambda^{(n)}H^{(n)}$ . Therefore, we
first consider the distribution of $H^{(n)}$ under the random service
suspension.

According to the discussion in Section 3.3, the actual headway for vehicle $l$
at station $n$ is
$H^{(n,l)}=\bar{H}+\frac{2\cdot\mathbb{E}[I^{(N,l)}]}{\bar{F}}+I^{(n,l)}-I^{(n,l-1)}$,
where $I^{(n,l)}$ is the total duration of incidents for vehicle $l$ during
its travel from the transportation hub to station $n$. Since $\bar{H}$ and
$\mathbb{E}[I^{(n,l)}]$ are constants, obtaining the headway distribution is
equivalent to quantifying the distribution of $I^{(n,l)}-I^{(n,l-1)}$.

Notice that $I^{(n,l)}$ and $I^{(n,l-1)}$ are i.i.d for all $l$ by our
assumption. It is useful to first consider the distribution of $I^{(n,l)}$.

###### Proposition 4.8

[Distribution of incident duration]: The total incident duration for vehicle
$l$ during its travel from the transportation hub to station $n$ (i.e.,
$I^{(n,l)}$) follows a compound Poisson-Exponential distribution with Poisson
rate $\gamma T^{(n)}$ and exponential rate $\theta$. Mathematically,

$\displaystyle I^{(n,l)}=\sum_{i=1}^{K}X_{i}\quad\quad\text{ where $X_{i}\sim$
{Exp}}(\theta)\;\;\forall i=1,...,K,\text{ and
}K\sim\text{{Poi}}(\gamma{T}^{(n)})$ (31)

The moment generating function (MGF) of a compound Poisson-Exponential
variable can be written as (Hogg et al. 2010)

$\displaystyle M_{I^{(n,l)}}(t)=\mathbb{E}[e^{tI^{(n,l)}}]=e^{\gamma
T^{(n)}(\frac{\theta}{\theta-t}-1)}\quad\quad\forall\;t<\theta$ (32)

Similarly, the MGF of $-I^{(n,l-1)}$ is

$\displaystyle M_{-I^{(n,l-1)}}(t)=\mathbb{E}[e^{-tI^{(n,l-1)}}]=e^{\gamma
T^{(n)}(\frac{\theta}{\theta+t}-1)}\quad\quad\forall\;t>-\theta$ (33)

From the MGF of $I^{(n,l)}$, we obtain $\mathbb{E}[I^{(N,l)}]=\frac{\gamma
T^{(N)}}{\theta}$. Then the headway equation (Eq. 12) becomes

$\displaystyle H^{(n,l)}=\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}}+I^{(n,l)}-I^{(n,l-1)}$ (34)

The following proposition provides the headway distribution:

###### Proposition 4.9

[MGF of headway]: Under the setting of this study, $\forall$ n = 1,..,N, the
MGF of $H^{(n)}$ can be expressed as

$\displaystyle M_{H^{(n)}}(t)=e^{t(\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}})}e^{\gamma T^{(n)}(\frac{2t^{2}}{\theta^{2}-t^{2}})}$
(35)

From the MGF of $H^{(n)}$, we can obtain the corresponding mean and variance
of headway as:

$\displaystyle\mathbb{E}[H^{(n)}]$ $\displaystyle=\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}}$ (36) $\displaystyle\text{Var}[H^{(n)}]$
$\displaystyle=\frac{4T^{(n)}\gamma}{\theta^{2}}$ (37)

###### Remark 4.10

The results show that random suspensions can increase the mean and variance of
headway. The impact on mean headway is through the increase in cycle time at
the route planning stage. The headway variance will increase with a higher
incident rate ($\gamma$) and higher average incident duration
($\frac{1}{\theta}$). Meanwhile, our model also captures the headway variance
propagation along stations as observed in many previous studies (Andersson and
Scalia-Tomba 1981, Hickman 2001): $\text{Var}[H^{(n)}]$ increase with the
station index $n$ (due to the increase in $T^{(n)}$).

However, the support of the derived headway distribution is $\mathbb{R}$,
meaning that $H^{(n)}$ can be negative due to the overtaking of vehicles. The
negative value of $H^{(n)}$ will cause problems in the definition of $Y^{(n)}$
(i.e., the number of arrival passengers within a headway). To address this
problem, we assume that drivers are not allowed to overtake the preceding
vehicles. This is true for the subway systems. Many transit agencies also use
this policy for bus operations. Given this assumption, the support of
$H^{(n)}$ becomes $[0,+\infty]$. Whenever $H^{(n)}<0$, the actual headway will
be 0 since the successor vehicle will not pass through the predecessor and
they will arrive at the station simultaneously (i.e., bus bunching). Hence,
the new truncated headway, denoted as $\hat{H}^{(n)}$, has a zero-inflation
mixture distribution:

$\hat{H}^{(n)}=\left\\{\begin{aligned} 0&\quad\text{ if $H^{(n)}\leq 0$}\\\
H^{(n)}&\quad\text{ otherwise}\end{aligned}\right.$ (38)

The zero-inflation truncated headway distribution is also observed in the
previous empirical study assuming no overtaking (Bellei and Gkoumas 2010).

However, to the best of the author’s knowledge, there is no closed-form MGF
for $\hat{H}^{(n)}$ because the difference between two compound Poisson-
exponential random variables has no closed-form probability density function.
Therefore, to have a tractable headway distribution, we have to approximate
$H^{(n)}$ with other distributions for which the corresponding zero-inflation
truncated distribution has analytical MGF.

In this study, we approximate the distribution of $H^{(n)}$ with normal
distribution for two reasons: 1) $I^{(n,l)}$ can be seen as the summation of a
large number of i.i.d random variables (Eq. 31) when the incident frequency is
high (i.e., $K$ is large, which is true for this study because we are
considering high-frequency short random disturbance). Hence, from the Central
Limit Theorem (CLT), we may approximate $I^{(n,l)}$ as normally distributed,
which leads to $H^{(n)}$ being normally distributed as well. Approximating the
headway disturbance as a normal random variable with the CLT was also used in
Daganzo (2009). 2) After approximating $H^{(n)}$ as a normal random variable
(denoted as $H^{(n)}_{\text{Normal}}$) with the same mean and variance, the
first three moments of $H^{(n)}$ and $H^{(n)}_{\text{Normal}}$ are the same.
This shows that the distribution of $H^{(n)}$ is similar to normal. Appendix
16 shows detailed derivations and numerical experiments to validate the
approximation.

Now let us consider a zero-inflation truncated distribution of
$H^{(n)}_{\text{Normal}}$ with support $[0,+\infty]$ and a probability mass
concentrated at zero. Denote the truncated random variable as
$\hat{H}^{(n)}_{\text{Normal}}$.

###### Proposition 4.11

[MGF of approximated headway]: Under the setting of this study, $\forall$ n =
1,..,N, the MGF of $\hat{H}^{(n)}_{\text{Normal}}$ can be expressed as

$\displaystyle
M_{\hat{H}^{(n)}_{\text{Normal}}}(t)=\Phi\left(\frac{-(\bar{H}\theta+\frac{2\gamma
T^{(N)}}{\bar{F}})}{2\sqrt{T^{(n)}\gamma}}\right)+e^{t(\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}})}e^{\gamma
T^{(n)}(\frac{2t^{2}}{\theta^{2}})}\left[1-\Phi\left(\frac{-(\bar{H}\theta+\frac{2\gamma
T^{(N)}}{\bar{F}})}{2\sqrt{T^{(n)}\gamma}}-\frac{2t\sqrt{T^{(n)}\gamma}}{\theta}\right)\right]$
(39)

where $\Phi(\cdot)$ is the cumulative density function (CDF) of a standard
normal distribution.

Based on the MGF of $\hat{H}^{(n)}_{\text{Normal}}$, notice that
$\left[1-\Phi\left(\frac{-\mu}{\sigma}\right)\right]=\Phi\left(\frac{\mu}{\sigma}\right)$,
we can get the corresponding mean and variance as follows.

$\displaystyle\mathbb{E}[\hat{H}^{(n)}_{\text{Normal}}]$
$\displaystyle=\mu\cdot\Phi\left(\frac{\mu}{\sigma}\right)+\sigma\cdot\phi\left(\frac{-\mu}{\sigma}\right)$
(40) $\displaystyle\text{Var}[\hat{H}^{(n)}_{\text{Normal}}]$
$\displaystyle=\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}$
(41)

where $\phi(\cdot)$ is the probability density function (PDF) of a standard
normal distribution. It is not clear how incidents will affect the mean
headway from Eq. 40 directly. However, the following proposition shows that
the mean headway increases as incident frequency ($\gamma$) and average
incident duration ($\frac{1}{\theta}$) increase.

###### Proposition 4.12

[Impact of incidents on headway]: The mean of the zero-inflation truncated
headway (i.e, either $\mathbb{E}[\hat{H}^{(n)}]$ or
$\mathbb{E}[\hat{H}^{(n)}_{\text{Normal}}]$) increases with the increase in
incident intensity (i.e., increase in $\gamma$ or $\frac{1}{\theta}$, or
both).

Proposition 4.12 is useful for the analysis of system stability with respect
to incidents, which is shown in Section 4.4.

####  4.3.5 Distribution of $Y^{(n)}$

The distribution of $Y^{(n)}$ is derived by assuming the headway is
$\hat{H}^{(n)}_{\text{Normal}}$ (instead of $H^{(n)}$, which may be negative).
To derive the PGF of $Y^{(n)}$, the following lemma is introduced.

###### Lemma 4.13

For two arbitrary random variable $U$ and $V$, assume that

  * •

there is a $\delta>0$ such that for $t$ in $(-\delta,\delta)$, the MGF of
$U|V$ is $M_{U|V}(t)=C_{1}(t)e^{C_{2}(t)V}$, where $C_{1}(t)$ and $C_{2}(t)$
are finite functions of $t$ that do not depend on $V$,

  * •

and the MGF of $V$, $M_{V}(\cdot)$, exists and $M_{V}[C_{2}(t)]$ is finite for
$t$ in $(-\delta,\delta)$.

Then the MGF of $U$ is given by

$\displaystyle M_{U}(t)=C_{1}(t)M_{V}[C_{2}(t)],\quad\quad-\delta<t<\delta.$
(42)

The proof of Lemma 4.13 can be found in Villa and Escobar (2006) Result 1.

###### Proposition 4.14

[PGF of $Y^{(n)}$]: Under the setting of this study, $\forall$ n = 1,..,N, the
PGF of $Y^{(n)}$, $Y(z)$, can be expressed as

$\displaystyle
Y(z)=\Phi\left(\frac{-\mu}{\sigma}\right)+e^{\mu\lambda^{(n)}(z-1)+\frac{\sigma^{2}(\lambda^{(n)}z-\lambda^{(n)})^{2}}{2}}\left[1-\Phi\left(\frac{-\mu}{\sigma}-\sigma\lambda^{(n)}(z-1)\right)\right]$
(43)

where $\mu=\bar{H}+\frac{2\gamma T^{(N)}}{\theta\bar{F}}$ and
$\sigma=\frac{2\sqrt{T^{(n)}\gamma}}{\theta}$ are the mean and standard
deviation of ${H}^{(n)}_{\text{Normal}}$, respectively.

From Eq. 105, we can obtain $\bar{Y}^{(n)}$, $\bar{\bar{Y}}^{(n)}$, and
$\bar{\bar{\bar{Y}}}^{(n)}$ by taking corresponding derivatives. The
expression of $\bar{Y}^{(n)}$ is shown below. The expressions for
$\bar{\bar{Y}}^{(n)}$ and $\bar{\bar{\bar{Y}}}^{(n)}$ are complicated and thus
omitted.

$\displaystyle\bar{Y}^{(n)}$
$\displaystyle=\left(\mu\cdot\Phi\left(\frac{\mu}{\sigma}\right)+\sigma\cdot\phi\left(\frac{-\mu}{\sigma}\right)\right)\cdot\lambda^{(n)}$
(44)

####  4.3.6 Solving for the roots

With the expression of $Y(z)$, the only unknown parts for the queue length
calculation (Eq. 25) are $z_{0}^{*},...,z_{C-1}^{*}$, which can be obtained by
solving the nonlinear equation $\textsc{Den}(z)=0$ (see Section 4.3 for
details). It is well known in the queuing literature that solving for the
roots of $\textsc{Den}(z)$ is practically difficult because typical
optimization algorithms usually only find only one root, while we need to find
all $C$ roots within the unit circle. This is especially changeling for $Y(z)$
with complex expressions because the objective function can be highly
nonlinear (such as $Y(z)$ in this study).

In this study, we propose an interpolation-based searching algorithm to
efficiently find all roots of $\textsc{Den}(z)$ within the unit circle. The
key idea is to intelligently set up different initial values for a general
root-solving algorithm (such as trust-region and Levenberg-Marquardt
algorithms) and let the algorithm converge to different solutions (i.e.,
different roots). We elaborate on the algorithm details in Appendix 18. The
numerical testing shows that our algorithm is able to find desired roots for
all testing scenarios in Section 5.1. It outperforms the methods in Powell
(1985) and Wilson (2014), where both of them have cases of not being able to
find all roots.

###  4.4 Stability condition

For all the derivations above, we assume that the steady-state distributions
of all variables exist. This triggers the discussion about the stability
condition. At the station level, the stability condition is described in
Proposition 4.15.

###### Proposition 4.15

[Stability condition]: Under the setting of this study, the bulk-service
queuing system at station $n$ is stable if and only if

$\displaystyle\rho^{(n)}=\frac{\bar{Y}^{(n)}}{\bar{S}^{(n)}}=\frac{\left(\mu\cdot\Phi\left(\frac{\mu}{\sigma}\right)+\sigma\cdot\phi\left(\frac{-\mu}{\sigma}\right)\right)\cdot\lambda^{(n)}}{\sum_{u=0}^{C}s_{u}^{(n)}u}=\frac{\lambda^{(n)}\cdot\mathbb{E}[\hat{H}^{(n)}_{\text{Normal}}]}{\sum_{u=0}^{C}s_{u}^{(n)}u}<1$
(45)

where $\rho^{(n)}$ is the utilization ratio for station $n$.

Proposition 4.15 is intuitive as it indicates that station $n$ is stable if
the average number of passengers arrived within a headway is smaller than the
average available capacity for each arrival vehicle (after alighting). From
Proposition 4.8, we know that a higher rate of incidents (i.e., larger
$\gamma$) and higher duration of incidents (i.e., higher $\frac{1}{\theta}$)
increase $\mathbb{E}[\hat{H}^{(n)}_{\text{Normal}}]$, which makes the system
more likely to be unstable. Hence, the above result quantifies the throughput
loss due to incidents. There are some remarks for Proposition 4.15.

###### Remark 4.16

As $\rho^{(n)}$ depends on $s^{(n)}$ and $s^{(n)}$ depends on the roots (i.e.,
$z_{0}^{*},...,z_{C-1}^{*}$) at station $n$, there is no direct way to judge
the stability at station $n$ without iterating the previous $n-1$ stations.
But for the first station ($n=1$), we have $s^{(1)}_{C}=1$ and $s^{(1)}_{u}=0$
for all $u=0,...,C-1$. Then Eq. 45 reduces to
$\rho^{(1)}=\frac{\lambda^{(n)}\cdot\mathbb{E}[\hat{H}^{(n)}_{\text{Normal}}]}{C}$,
which can be used to assess the stability directly.

###### Remark 4.17

Proposition 4.15 only discusses the stability at the station level. At the
route level, a route is considered stable if “all stations in the route are
stable”. Mathematically, a route is stable if and only if
$\rho^{(n)}<1,\forall\;n=1,2,...,N$.

###### Remark 4.18

It is worth discussing the relationship of stability of stations $n$ and
$n-1$. If station $n-1$ is stable, then $s^{(n)}$ can be calculated as
described in Section 4.1, and the stability of station $n$ can be evaluated
accordingly. However, if station $n-1$ is not stable, station $n$ may be
stable because there may be passengers alighting at station $n$. For this
situation, we have $v^{(n-1)}_{C}=1$ and $v^{(n-1)}_{k}=0$ for all
$k=0,1,...,C-1$. Then $s^{(n)}$ is determined by the alighting rate at station
$n$. It is easy to verify that in this situation
$\bar{S}^{(n)}=\alpha^{(n)}C$. And the stability condition is
$\rho^{(n)}=\frac{\lambda^{(n)}\cdot\mathbb{E}[\hat{H}^{(n)}_{\text{Normal}}]}{\alpha^{(n)}C}<1$.

###  4.5 Summary of calculation procedure

So far, we have derived the calculation process for key variables of interest.
Algorithm 1 summarizes the calculation procedure, which iterates through the
$N$ stations of the route. This is more efficient and provides more analytical
insights than a simulation model.

Algorithm 1 Performance indicators calculation procedure

1:Initialize $v^{(0)}_{0}=1$ and $v^{(0)}_{k}=0\quad\quad\forall k=1,...C$.

2:for $n=1:N$ do

3: $g^{(n)}=v^{(n-1)}A^{(n)}$ $\triangleright$ Eq. 14

4: $s_{k}^{(n)}=1-g_{C-k}^{(n)}\quad\quad\forall k=0,1,...,C$ $\triangleright$
Eq. 13

5: Calculate $\bar{S}^{(n)},\bar{\bar{S}}^{(n)},$ and
$\bar{\bar{\bar{S}}}^{(n)}$ based on $s^{(n)}$.

6: Calculate $\bar{Y}^{(n)},\bar{\bar{Y}}^{(n)},$ and
$\bar{\bar{\bar{Y}}}^{(n)}$ $\triangleright$ Section 4.3.4

7: if $\bar{Y}^{(n)}<\bar{S}^{(n)}$ then $\triangleright$ Station $n$ is
stable

8: Solve the roots $z_{0}^{*},...,z^{*}_{C-1}$ for the denominator of $Q(z)$
in Eq. 18$\triangleright$ Section 4.3.1

9: Calculate $q_{0}^{(n)},...,q_{C-1}^{(n)}$ based on
$z_{0}^{*},...,z^{*}_{C-1}$ $\triangleright$ Section 4.3.2

10: Calculate $\mathbb{E}[Q^{(n)}],\text{Var}[Q^{(n)}],\mathbb{E}[W^{(n)}]$,
and $\text{Var}[W^{(n)}]$ $\triangleright$ Eq. 25 \- 28

11: $v^{(n)}=g^{(n)}B^{(n)}$ $\triangleright$ Eq. 16. $B^{(n)}$ is a function
of $q_{0:C-1}^{(n)}$

12: else$\triangleright$ Station $n$ is not stable

13: $q^{(n)}_{k}=0\quad\quad\forall k=0,1,...,C-1$

14: Set $\mathbb{E}[Q^{(n)}],\text{Var}[Q^{(n)}],\mathbb{E}[W^{(n)}]$, and
$\text{Var}[W^{(n)}]$ to infinity

15: $v^{(n)}_{C}=1$ and $v^{(n)}_{k}=0\quad\quad\forall k=0,1,...,C-1$

##  5 Numerical example

###  5.1 Experimental design

To test the proposed framework, we use an example bus route adapted from Islam
et al. (2015) and Hickman (2001). There are 10 stations and the attributes for
each station are shown in Table 1. The layout of the bus route is shown in
Figure 6, where we assume the no-incident travel time between two consecutive
stations is 5 minutes, the total cycle time without incident is $\bar{E}=100$
min, and travel time from the transportation hub to the last station is
$T^{(N)}=50$ minutes.

Table 1: Example bus system parameters

Station ID |  $\lambda^{(n)}$ (passengers/min) | $\alpha^{(n)}$ | Station ID |  $\lambda^{(n)}$ (passengers/min) | $\alpha^{(n)}$
---|---|---|---|---|---
1 | 0.75 | 0 | 6 | 1 | 0.8
2 | 1.5 | 0 | 7 | 0.75 | 0.5
3 | 0.75 | 0.1 | 8 | 0.5 | 0.1
4 | 3 | 0.25 | 9 | 0.2 | 0.75
5 | 1.5 | 0.25 | 10 | 0 | 1

(a)  Figure 6: Case study route layout

To test the sensitivity of performance indicators to different parameters, we
consider different values of $C$, $\theta$, $\gamma$, $\bar{H}$, and demand
(Table 2). The demand is adjusted by a scaling factor that is applied to the
arrival rates $\lambda^{(n)}$ in Table 1. The fleet size $\bar{F}$ is
determined as $\frac{\bar{E}}{\bar{H}}$. When the sensitivity testing is
conducted for one parameter (e.g., $C$), other parameters (e.g., $\theta$,
$\gamma$, $\bar{H}$, and the demand factor) are set to their reference values
for comparison.

Table 2: Scenario design

Parameters | Value space | Reference value
---|---|---
$C$ | {30, 34, 38} | 34
$\gamma$ (/min) | {0, 1/10, 1/5, 1/3} | 1/5
$\theta$ (/min) | {2, 1 ,1/2} | 1
$\bar{H}$ (min) | {2, 4, 7} | 6
Demand factor | {0.2, 0.4, 0.6, 0.8, 1} | 0.8

###  5.2 Performance indicators

The mean and standard deviation of queue length for each station under
different testing scenarios are shown in Figure 7. Generally, for all
scenarios, the queue length patterns are consistent with the congestion
patterns we expect given the passenger arrival and alighting rates. That is,
the expected queue length is relatively higher at stations 2 and 8. The
expected queue length at the last station is always zero as its passenger
arrival rate is 0.

Figure 7a shows the queue length patterns with respect to bus capacity. The
system is not very sensitive to bus capacity. The reason is that under the
reference scenario, the system is not congested and capacity is not fully
utilized. Thus, increasing capacity does not affect the queuing distribution.
Figure 7b shows the impact of incident occurrence rate $\gamma$ on queue
length. When there is no random suspension in the system ($\gamma=0$), the
expected queue length at station 8 is 4.5. As the frequency of incidents
increases, the system becomes more congested with longer expected queue
lengths and higher variances. When the incident frequency increases to $1/3$
per minute on average ($\gamma=1/3$), the expected queue length at station 8
is increased to 8.3 units. Similar results can be observed for the duration of
incidents (Figure 7c). When the average incident duration is 30 seconds
($\theta=2$), $\mathbb{E}(Q^{(8)})=5.0$. When the average incident duration is
2 minutes ($\theta=1/2$), $\mathbb{E}(Q^{(8)})$ increases to 12.6. The impacts
of $\theta$ and $\gamma$ on queue length are both more significant at crowded
stations. The impact of $\bar{H}$ is shown in Figure 7d. As expected, higher
headway means a lower service rate and thus a higher expected queue length. As
$\bar{H}$ increases from 2 minutes to 7 minutes, the queue length at station 8
increases from 4.1 to 9.7. The impact of the demand factor (Figure 7e) shows
similar patterns. As the demand factor increases from 0.5 to 1.0, the queue
length at station 8 increases from 4.1 to 8.3. The impact of $\bar{H}$ and the
demand factor are relatively similar for crowded and uncrowded stations.

(a) Sensitivity on $C$

(b) Sensitivity on $\gamma$

(c) Sensitivity on $\theta$

(d) Sensitivity on $\bar{H}$

(e) Sensitivity on demand factor

Figure 7: Mean and standard deviation of queue length (the shaded part is
$0.2\times$standard deviation)

Figure 8 shows the mean and standard deviation of passenger waiting time for
the different scenarios. We observe that the downstream stations generally
have higher waiting time expectations and variances due to the headway
variance propagation. For congested stations, such as stations 3 and 8, extra
waiting times are observed due to passengers left behind with capacity
constraints.

Figure 8a shows the impact of capacity on waiting time. Similar to the results
on queue length, the impact is not very significant. The impacts of $\gamma$
and $\theta$ on waiting times are shown in Figure 8b and 8c, respectively. As
increases in $\gamma$ and $1/\theta$ result in an increase in expected
headway, the mean waiting times at all stations are increased. The impacts on
crowding stations are more significant. When $\gamma=0$, there is no incident
in the system. In this case, there are no left behind or headway
irregularities at any stations and their expected waiting times are all equal
to 2 minutes (i.e., $\frac{1}{2}\bar{H}$, as no incidents mean all stations
have the same fixed headway). When $\gamma$ increases to $1/5$, station 3 has
left behind passengers and the waiting time is increased to 4.6 minutes. When
$\theta$ decreases (i.e., mean incident duration increases) from $2$ to $1/2$,
the expected waiting time at station 8 increases from 3.0 to 11.8 minutes.
Changes in $\bar{H}$ have the most direct impact on the expected waiting time.
The increase in planned headway causes an increase in waiting time for all
stations. There are a few left-behind passengers observed at stations 3 and 8
when $\bar{H}=7$ min. Finally, as demand increases, the waiting time increases
only if there are left behind (e.g., when demand factor $=1$) because it does
not change the headway distribution. At station 3, the increase in the demand
factor from 0.5 to 1.0 results in an increase in the expected waiting time
from 3.5 to 4.2 minutes.

(a) Sensitivity on $C$

(b) Sensitivity on $\gamma$

(c) Sensitivity on $\theta$

(d) Sensitivity on $\bar{H}$

(e) Sensitivity on demand factor

Figure 8: Mean and standard deviation of waiting time (the shaded part is
$0.2\times$standard deviation)

###  5.3 Comparison between simulated and theoretical results

To validate the theoretical results, we develop a simulation model to
calculate the expectation and variance of queue length and waiting time. The
simulation procedure is shown in Appendix 20.

We compare the simulation and theoretical results for the reference parameter
setting (Table 2). A total of $50,000$ vehicle runs are simulated. The
comparisons of mean and standard deviation for queue length and waiting time
are shown in Figure 9. We observe that the simulation and theoretical results
match well, validating the theoretical model’s correctness. However, the
theoretical results slightly overestimate the mean and variance of the queue
length and waiting time. The main reason may be the approximation of headway
distribution as normal. As shown in Figure 11, the actual headway has more
probability density concentrated at the mean (i.e., more peakedness), implying
that the actual headway has less probability of deviating from the planned
one, thus the simulation scenario may have a smaller queue length and waiting
time.

(a) $\mathbb{E}[Q^{(n)}]$

(b) $\text{Std dev}[Q^{(n)}]$

(c) $\mathbb{E}[W^{(n)}]$

(d) $\text{Std dev}[W^{(n)}]$

Figure 9: Comparison between simulation and theoretical results (reference
scenario)

##  6 Conclusion and discussion

This paper proposes a stochastic framework to evaluate the performance of PTSs
under short random service suspensions. Specifically, we analyze the system
stability conditions and derive closed-form formulations for the mean and
variance of queue length and waiting time at each station. The derived
stability conditions are intuitive and imply that the system is more likely to
be unstable with high incident rates, high incident duration, high demand, low
service frequency, and low vehicle capacity. The proposed model is implemented
using an example bus network adapted from the literature. A sensitivity
analysis of different parameters (such as incident rate, incident duration,
vehicle capacity, etc.) was conducted. The results show that congested
stations (i.e., stations with high demand rates) are more vulnerable to random
service suspensions. The results are validated with a simulation model,
showing consistency between theoretical and simulation outcomes.

The proposed model has several potential applications. 1) It can facilitate
the design and planning of PTS with the consideration of random system
interruptions, such as the design of headways and the determination of vehicle
capacity. Moreover, the estimated queue length can be used to evaluate the
layout and capacity of congested stations. 2) The model can be used to monitor
system performance and identify critical stations by inputting the historical
demand and incident information. 3) The model can support efficient cost-
benefit analysis of approaches to improve services using estimates of waiting
time and queue length. For example, the model can answer that, to control the
waiting time within a threshold, what is the most cost-effective way (e.g.,
increase vehicle size, decrease headway, or increase maintenance frequency to
reduce the random suspension rate). In summary, the efficient calculation of
the system’s performance indicators can be used in public transit planning,
operations, and management applications.

Future studies can address a number of aspects. First, the model can be
extended from the route-level to the network-level. The main difference
between route level and network level is the consideration of transfer
passengers. A straightforward way is to incorporate the transfer demand as
part of the arrival demand. But additional transfer-related parameters (which
should be connected with the alighting rate) need to be specified in the
model. Second, like many previous random service disruption papers (Section
2.2), the model can be extended to consider partial interruptions (as opposed
to fully stopped as assumed in this paper). With partial interruptions,
vehicles can still have positive speed at the disruption state. The headway
distribution assumption needs to be revised. Third, as mentioned before, this
paper assumes no balking and reneging behavior of passengers. Future studies
may extend the model by considering a more complicated passenger-side
behavior.

This work was in part supported by SJTU UM Joint Institute, J. Wu & J. Sun
Endowment Fund, C2SMART University Transportation Center, and NYU Tandon
School of Engineering.

##  7 Proof of Proposition 4.1

When vehicle $l$ arrives at station $n$, by definition, there are
$V^{(n-1,l)}$ number of passengers in the vehicle. Given that there are $i$
passengers on board when vehicle $l$ arrives at station $n$, let the
probability that there are $j$ passengers remaining on the vehicle be
$a_{ij}^{(n)}$. $a_{ij}^{(n)}$ also represents the probability of $i-j$
passengers alighting, which follows a binomial distribution with parameters
$i$ and $\alpha^{(n)}$ (if $i\geq j$ and $i,j\neq 0$). Hence, $a_{ij}^{(n)}$
can be expressed as Eq. 15. Then we have

$\displaystyle g^{(n,l)}=v^{(n-1,l)}A^{(n)}$ (46)

where $v^{(n-1,l)}=[v^{(n-1,l)}_{0},...,v^{(n-1,l)}_{C}]\in\mathbb{R}^{C+1}$,
$g^{(n,l)}=[g^{(n,l)}_{0},...,g^{(n,l)}_{C}]\in\mathbb{R}^{C+1}$. According to
the relationship between $S^{(n,l)}$ and $G^{(n,l)}$ as shown in Eq. 5, the
distribution of the number of available spaces after alighting is simply

$\displaystyle s_{k}^{(n,l)}=g_{C-k}^{(n,l)}\quad\quad\forall k=0,1,...,C$
(47)

Note that Eq. 46 and 47 hold for all $l$. Since we assume the steady state
distributions exist, letting $l\to\infty$ on both sides of Eq. 46 and 47
completes the proof.

##  8 Proof of Proposition 4.2

Let $b_{ij}^{(n,l)}$ be the probability that the load of vehicle $l$ is $j$
after passenger boarding given that there are $i$ passengers onboard after
alighting (i.e., $G^{(n,l)}=i$) at station $n$. Hence, if $0\leq i\leq j<C$,
$b_{ij}^{(n,l)}$ is simply the probability that there are $j-i$ passengers in
the queue (such that after boarding there are $j$ passengers on the vehicle):

$\displaystyle b_{ij}^{(n,l)}=q^{(n,l)}_{j-i},$ $\displaystyle\quad\text{ if
$0\leq i\leq j<C$}.$ (48)

If $j=C$ and $0\leq i<C$, the vehicle reaches capacity after boarding. Then
$b_{ij}^{(n,l)}$ should be the probability that the number of passengers in
the queue is greater than or equal to $C-i$ (i.e., one minus the probability
that there are less than or equal to $C-i-1$ passengers in the queue). This
leads to:

$\displaystyle b_{ij}^{(n,l)}=1-\sum_{k=0}^{C-i-1}q^{(n,l)}_{k},$
$\displaystyle\quad\text{ if $j=C$ and $0\leq i<C$}.$ (49)

When $i=j=C$, we simply have $b_{ij}^{(n,l)}=1$ because regardless of the
number of waiting passengers in the queue, nobody can board as the vehicle is
full. Given Eq. 17, the vehicle load distribution can be calculated as

$\displaystyle v^{(n,l)}_{j}=\sum_{i=1}^{C}g_{i}\cdot
b_{ij}^{(n,l)}\quad\quad\forall j=0,1,...,C$ (50)

Notice that Eq. 50 holds for all $l$. As we assume that the steady-state
distributions exist, taking $l\to\infty$ for both sides of Eq. 50 and
rewriting it in a matrix form completes the proof.

##  9 Proof of Proposition 4.3

The proof follows a similar idea in Powell (1981). The difference from Powell
(1981) is that we consider an arbitrary vehicle capacity distribution
$s_{0}^{(n)},...,s_{C}^{(n)}$, while in Powell (1981) the capacity is fixed.
Note that Powell (1985) provided an equivalent formulation as Eq. 18 with
variable vehicle capacities using the transform of $S^{(n)}$.

From the relationship in Eq. 3, we have

$\displaystyle q^{(n,l+1)}_{k}=\sum_{i=0}^{k}r^{(n,l)}_{i}y_{k-i}^{(n,l)}$
(51)

where $r_{k}^{(n,l)}:=\mathbb{P}(R^{(n,l)}=k)$ and
$y_{k}^{(n,l)}:=\mathbb{P}(Y^{(n,l)}=k)$ for all non-negative integers $k$.
Let $S^{(n,l)}$ be the number of available spaces for train $l$ when it
arrives at station $n$. Given that $S^{(n,l)}=u$, if $u$ is greater than or
equal to $Q^{(n,l)}$, all passengers can board, and there is no left behind.
Then we have

$\displaystyle\left.r^{(n,l)}_{0}\right|_{S^{(n,l)}=u}=\mathbb{P}(u\geq
Q^{(n,l)})=\sum_{k=0}^{u}q^{(n,l)}_{k}$ (52)

where
$\left.r^{(n,l)}_{k}\right|_{S^{(n,l)}=u}=\mathbb{P}(R^{(n,l)}=k\;|\;S^{(n,l)}=u)$.
If $u$ is less than $Q^{(n,l)}$, only $u$ passengers can board and there are
$Q^{(n,l)}-u$ number of left-behind passengers. So

$\displaystyle\left.r^{(n,l)}_{i}\right|_{S^{(n,l)}=u}=\mathbb{P}(Q^{(n,l)}-u=i)=q^{(n,l)}_{u+i}\quad\quad\;i=1,2,...$
(53)

Based on Eq. 52 and 53, Eq 51 can be reformulated as

$\displaystyle
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)$
(54)

Assume the steady-state probabilities for all variables exist, we have
$\lim_{l\to\infty}q^{(n,l)}_{k}=q^{(n)}_{k}$,
$\lim_{l\to\infty}s^{(n,l)}_{k}=s^{(n)}_{k}$,
$\lim_{l\to\infty}y^{(n,l)}_{k}=y^{(n)}_{k}$. Taking the limit of $l$ for both
sides of Eq. 54 leads to

$\displaystyle
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)}$
(55)

Assume the probability generating function (PGF) for $Q^{(n)}$, $R^{(n)}$ and
$Y^{(n)}$ are $Q(z)$, $R(z)$, and $Y(z)$, respectively, where $Q^{(n)}$,
$R^{(n)}$ and $Y^{(n)}$ and the steady state random variables of $Q^{(n,l)}$,
$R^{(n,l)}$ and $Y^{(n,l)}$. Hence,

$\displaystyle Q(z)$ $\displaystyle=\sum_{k=0}^{\infty}q^{(n)}_{k}z^{k}$ (56)
$\displaystyle R(z)$
$\displaystyle=\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|_{S^{(n)}=u}$
(57) $\displaystyle Y(z)$ $\displaystyle=\sum_{k=0}^{\infty}y^{(n)}_{k}z^{k}$
(58)

Substituting Eq. 52 and 53 into 57 results in

$\displaystyle R(z)$
$\displaystyle=\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)}$
(59)
$\displaystyle=\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]$
(60)

Notice that $Q^{(n)}=R^{(n)}+Y^{(n)}$ (this is obtained by taking the limit of
$l$ for Eq. 3). Since $R^{(n)}$ and $Y^{(n)}$ are independent, we have

$\displaystyle Q(z)=R(z)Y(z)$ (61)

Combining Eq. 60 and 61 obtains

$\displaystyle Q(z)$
$\displaystyle=\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}}}$
$\displaystyle=\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}}$
(62)

##  10 Proof of Proposition 4.4

Though $q_{0}^{(n)},...,q_{C-1}^{(n)}$ can be obtained by solving $C$ system
equations as mentioned in Section 4.3, we attempt to provide a more direct way
to calculate $q_{0:C-1}^{(n)}$ in this section.

Using the fact that the numerator of $Q(z)$ is in the polynomial order of $C$,
$Q(z)$ can be reformulated in terms of $z_{1}^{*},...,z_{C-1}^{*}$ as:

$\displaystyle Q(z)$
$\displaystyle=\frac{(z-1)\prod_{i=1}^{C-1}(z-z_{i}^{*})\sum_{u=0}^{C}s^{(n)}_{u}\sum_{i=0}^{u}q_{i}^{(n)}}{\frac{z^{C}}{Y(z)}-\sum_{u=0}^{C}s_{u}^{(n)}z^{C-u}}$
(63)

When $z\to 1$, both $\textsc{Num}(z)$ and $\textsc{Den}(z)$ approach 0. We
also have the fact that $\lim_{z\to 1}Q(z)=1$. Therefore, using L’Hopital’s
rule we have

$\displaystyle\lim_{z\to 1}Q(z)=1=\lim_{z\to
1}\frac{\textsc{Num}^{\prime}(z)}{\textsc{Den}^{\prime}(z)}=\frac{\prod_{i=1}^{C-1}(1-z_{i}^{*})\sum_{u=0}^{C}s^{(n)}_{u}\sum_{i=0}^{u}q_{i}^{(n)}}{\sum_{u=0}^{C}s_{u}^{(n)}u-Y^{\prime}(1)}$
(64) $\displaystyle\Rightarrow$
$\displaystyle\sum_{u=0}^{C}s^{(n)}_{u}\sum_{i=0}^{u}q_{i}^{(n)}=\frac{\sum_{u=0}^{C}s_{u}^{(n)}u-Y^{\prime}(1)}{\prod_{i=1}^{C-1}(1-z_{i}^{*})}$
(65)

Define $\bar{Y}^{(n)}:=Y^{\prime}(1)=\mathbb{E}[Y^{(n)}]$ as the mean number
of arrival passengers within a headway at station $n$,
$\bar{S}^{(n)}:=\sum_{u=0}^{C}s^{(n)}_{u}u=$ as the mean number of available
spaces in an arriving bus at station $n$. Substituting Eq. 65 into 63, $Q(z)$
can be rewritten as

$\displaystyle Q(z)$
$\displaystyle=\frac{(\bar{S}^{(n)}-\bar{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}}$
(66)

Comparing Eq. 66 and 18, let the numerators of two equations be equal, we have

$\displaystyle(\bar{S}^{(n)}-\bar{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]$
(67)

As the LHS and RHS of Eq. 67 are both polynomials about $z$, the coefficients
of each polynomial in $z$ must be equal. By matching the coefficients of
$z^{0}$, we have

$\displaystyle(-1)(\bar{S}^{(n)}-\bar{Y}^{(n)})\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1}=-q_{0}^{(n)}s_{C}^{(n)}$
(68)

which leads to

$\displaystyle
q_{0}^{(n)}=\frac{1}{s_{C}^{(n)}}(\bar{S}^{(n)}-\bar{Y}^{(n)})\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1}$
(69)

Note that $\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1}$ is always greater
than 0 because 1) when $C$ is odd, as the complex roots appear as conjugates,
$\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1}>0$. 2) when $C$ is even,
besides $z_{0}^{*}=1$, there exists another real root on the negative real
axis (denoted as $z_{\frac{C}{2}}^{*}$, where $-1\leq z_{\frac{C}{2}}^{*}<0$).
So, we have $\frac{z_{\frac{C}{2}}^{*}}{z_{\frac{C}{2}}^{*}-1}>0$, which leads
to $\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1}>0$.

To validate Eq. 67, consider the fixed capacity situation where
$s_{C}^{(n)}=1$ and $\bar{S}^{(n)}=C$. Then Eq. 67 reduces to

$\displaystyle\left.q_{0}^{(n)}\right|_{s_{C}^{(n)}=1}=(C-\bar{Y}^{(n)})\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1}$
(70)

This is the same as Chaudhry et al. (1987). Now we will derive
$q_{1:C-1}^{(n)}$. Observing that the numerator of Eq. 66 can be rewritten as

$\displaystyle(\bar{S}^{(n)}-\bar{Y}^{(n)})(z-1)\prod_{i=1}^{C-1}\frac{z-z_{i}^{*}}{1-z_{i}^{*}}$
$\displaystyle=\frac{1}{s_{C}^{(n)}}(\bar{S}^{(n)}-\bar{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)}$
$\displaystyle=s_{C}^{(n)}q_{0}^{(n)}(z-1)\prod_{i=1}^{C-1}\left(1-\frac{z}{z_{i}^{*}}\right)$
$\displaystyle=-s_{C}^{(n)}q_{0}^{(n)}\prod_{i=0}^{C-1}\left(1-\frac{z}{z_{i}^{*}}\right)$
(71)

Define
$\prod_{i=0}^{C-1}\left(1-\frac{z}{z_{i}^{*}}\right):=\sum_{j=0}^{C}\eta_{j}z^{j}$,
where $\eta_{j}$ is the polynomial coefficient of $z^{j}$. For the RHS of Eq.
67, the polynomial coefficient of $z^{C-k}$ is
$-\sum_{u=k}^{C}s_{u}^{(n)}q_{u-k}^{(n)}$. And from Eq. 71, the polynomial
coefficient of $z^{C-k}$ is $-s_{C}^{(n)}q_{0}^{(n)}\eta_{C-k}$. Matching the
coefficient of the same order of $z$ leads to

$\displaystyle
s_{C}^{(n)}q_{0}^{(n)}\eta_{C-k}=\sum_{u=k}^{C}s_{u}^{(n)}q_{u-k}^{(n)}\quad\quad\;k=1,2,...,C-1$
(72)

To validate Eq. 72, consider the fixed capacity situation where
$s_{C}^{(n)}=1$ and $s_{k}^{(n)}=0,\forall\;0\leq k<C$. then Eq. 72 reduces to

$\displaystyle
q_{C-k}^{(n)}=q_{0}^{(n)}\eta_{C-k}\quad\quad\;k=1,2,...,C-1\quad\text{ if
}{s_{C}^{(n)}=1}$ (73)

which is the same as Chaudhry et al. (1987).

Eq. 72 can be expressed in a matrix form by adding
$s_{C}^{(n)}q_{0}^{(n)}\eta_{0}^{(n)}=s_{C}^{(n)}q_{0}^{(n)}$ (note that
$\eta_{0}^{(n)}=1$ by definition):

$\displaystyle\tilde{\eta}^{(n)}=q_{0:C-1}^{(n)}\Lambda^{(n)}$ (74)

where
$\tilde{\eta}^{(n)}=[s_{C}^{(n)}q_{0}^{(n)}\eta_{0}^{(n)},s_{C}^{(n)}q_{0}^{(n)}\eta_{1}^{(n)},...,s_{C}^{(n)}q_{0}^{(n)}\eta_{C-1}^{(n)}]\in\mathbb{R}^{C}$
and

$\Lambda^{(n)}=\left[\begin{array}[]{ccccccc}s_{C}^{(n)}&s_{C-1}^{(n)}&s_{C-2}^{(n)}&&...&&s_{1}^{(n)}\\\
0&s_{C}^{(n)}&s_{C-1}^{(n)}&&...&&s_{2}^{(n)}\\\
...&0&s_{C}^{(n)}&&...&&s_{3}^{(n)}\\\ 0&...&0&&...&&s_{4}^{(n)}\\\
0&0&...&&...&&...\\\ 0&0&0&&...&&s_{C}^{(n)}\\\
\end{array}\right]\in\mathbb{R}^{C\times C}$ (75)

As $s_{C}^{(n)}>0$ is a known condition, the triangular matrix $\Lambda^{(n)}$
is invertible. Thus, we have

$\displaystyle q_{0:C-1}^{(n)}=\tilde{\eta}^{(n)}(\Lambda^{(n)})^{-1}$ (76)

##  11 Proof of Proposition 4.5

The derivation follows the same idea in Powell (1981). These results are
equivalent to Powell (1985) which considered the general bulk-service queue
model (but Powell (1985) did not provide detailed proof in the paper).

Here we try to provide analytical formulations of $\mathbb{E}[Q^{(n)}]$ and
$\text{Var}[Q^{(n)}]$. The key is to find $Q^{\prime}(1)$ and
$Q^{\prime\prime}(1)$. The derivation follows a similar idea in Powell (1981).

Let
$A(z)=\frac{(\bar{S}^{(n)}-\bar{Y}^{(n)})(z-1)}{\frac{z^{C}}{Y(z)}-\sum_{u=0}^{C}s_{u}^{(n)}z^{C-u}}$
and $B_{i}(z)=\frac{z-z_{i}^{*}}{1-z_{i}^{*}}$, then
$Q(z)=A(z)\prod_{i=1}^{C-1}B_{i}(z)$. Based on the fact that $B_{i}(1)=1$ and
$Q(z)=1$, we must have $A(1)=1$. Hence,

$\displaystyle Q^{\prime}(1)$
$\displaystyle=A^{\prime}(1)B_{1}(1)...B_{C-1}(1)+A(1)B_{1}^{\prime}(1)...B_{C-1}(1)+...+A(1)B_{1}(1)...B_{C-1}^{\prime}(1)$
$\displaystyle=A^{\prime}(1)+\sum_{i=1}^{C-1}B_{i}^{\prime}(1)$ (77)

Since $B_{i}^{\prime}(1)=\frac{1}{1-z_{i}^{*}}$, the problem now becomes
finding $A^{\prime}(1)$. Again, let $A(z)=\frac{A_{1}(z)}{A_{2}(z)}$. Then,

$\displaystyle
A^{\prime}(z)=\frac{A_{1}^{\prime}(z)A_{2}(z)-A_{1}(z)A_{2}^{\prime}(z)}{(A_{2}(z))^{2}}$
(78)

Notice that when $z\to 1$, the numerator and denominator of $A^{\prime}(z)$
approach 0 (because $A_{1}(1)=0$ and $A_{2}(1)=0$). Therefore, applying
L’Hopital’s rule yields:

$\displaystyle
A^{\prime}(z)=\frac{A_{1}^{\prime\prime}(z)A_{2}(z)-A_{1}(z)A_{2}^{\prime\prime}(z)}{2A_{2}(z)A_{2}^{\prime}(z)}$
(79)

Again we have $0/0$ when $z\to 1$ because $A_{1}^{\prime\prime}(z)=0$ and
$A_{2}(1)=0$. Applying L’Hopital’s rule once more gives:

$\displaystyle
A^{\prime}(z)=\frac{-A_{1}^{\prime}(z)A_{2}^{\prime\prime}(z)-A_{1}(z)A_{2}^{\prime\prime\prime}(z)}{2A_{2}^{\prime}(z)A_{2}^{\prime}(z)+2A_{2}(z)A_{2}^{\prime\prime}(z)}$
(80)

Substituting $z=1$ leads to

$\displaystyle
A^{\prime}(1)=\frac{-A_{1}^{\prime}(1)A_{2}^{\prime\prime}(1)}{2(A_{2}^{\prime}(1))^{2}}$
(81)

Based on the fact that $Y(1)=1$, $Y^{\prime}(1)=\bar{Y}^{(n)}$,
$Y^{\prime\prime}(1)=\mathbb{E}[({Y}^{(n)})^{2}]-\bar{Y}^{(n)}$, we have

$\displaystyle A_{1}^{\prime}(1)$ $\displaystyle=\bar{S}^{(n)}-\bar{Y}^{(n)}$
(82) $\displaystyle A_{2}^{\prime}(1)$
$\displaystyle=\left.\frac{Cz^{C-1}Y(z)-Y^{\prime}(z)z^{C}}{Y(z)^{2}}-\sum_{u=0}^{C}(C-u)s_{u}^{(n)}z^{C-u-1}\right|_{z=1}=\bar{S}^{(n)}-\bar{Y}^{(n)}$
(83) $\displaystyle A_{2}^{\prime\prime}(1)$
$\displaystyle=\left.\frac{C{(C-1)}z^{C-2}}{Y(z)}-\frac{2Y^{\prime}(z)Cz^{C-1}}{Y(z)^{2}}-\frac{Y^{\prime\prime}(z)z^{C}}{Y(z)^{2}}+\frac{2Y^{\prime}(z)^{2}z^{C}}{Y(z)^{3}}-\sum_{u=0}^{C}(C-u)(C-u-1)s_{u}^{(n)}z^{C-u-2}\right|_{z=1}$
$\displaystyle=C(C-1)-2\bar{Y}^{(n)}C-\mathbb{E}[({Y}^{(n)})^{2}]+2(\bar{Y}^{(n)})^{2}+\bar{Y}^{(n)}-C^{2}+C+2C\bar{S}^{(n)}-\bar{S}^{(n)}-\mathbb{E}[(S^{(n)})^{2}]$
$\displaystyle=-2\bar{Y}^{(n)}C-\mathbb{E}[({Y}^{(n)})^{2}]+2(\bar{Y}^{(n)})^{2}+\bar{Y}^{(n)}+2C\bar{S}^{(n)}-\bar{S}^{(n)}-\mathbb{E}[(S^{(n)})^{2}]$
(84)

Substituting Eq. 82, 83, and 84 into Eq. 81 results in

$\displaystyle
A^{\prime}(1)=\frac{2\bar{Y}^{(n)}C+\mathbb{E}[({Y}^{(n)})^{2}]-2(\bar{Y}^{(n)})^{2}-\bar{Y}^{(n)}-2C\bar{S}^{(n)}+\bar{S}^{(n)}+\mathbb{E}[(S^{(n)})^{2}]}{2(\bar{S}^{(n)}-\bar{Y}^{(n)})}$
(85)

Therefore, we have

$\displaystyle\mathbb{E}[Q^{(n)}]=\frac{2\bar{Y}^{(n)}C+\mathbb{E}[({Y}^{(n)})^{2}]-2(\bar{Y}^{(n)})^{2}-\bar{Y}^{(n)}-2C\bar{S}^{(n)}+\bar{S}^{(n)}+\mathbb{E}[(S^{(n)})^{2}]}{2(\bar{S}^{(n)}-\bar{Y}^{(n)})}+\sum_{i=1}^{C-1}\frac{1}{1-z_{i}^{*}}$
(86)

To validate this formulation, let us consider a fixed capacity situation with
$s_{C}^{(n)}=1$. Then $\bar{S}^{(n)}=C$, $\mathbb{E}[(S^{(n)})^{2}]=C^{2}$.
Then Eq. 86 reduces to

$\displaystyle\left.\mathbb{E}[Q^{(n)}]\right|_{s_{C}^{(n)}=1}=\frac{C-C^{2}+2\bar{Y}^{(n)}C+\mathbb{E}[({Y}^{(n)})^{2}]-2(\bar{Y}^{(n)})^{2}-\bar{Y}^{(n)}+}{2(C-\bar{Y}^{(n)})}+\sum_{i=1}^{C-1}\frac{1}{1-z_{i}^{*}}$
(87)

which is equivalent to Powell (1981)’s.

According to Eq. 24, the key to obtain $\text{Var}[Q^{(n)}]$ is to calculate
$Q^{\prime\prime}(1)$. Taking the logarithm of
$Q(z)=A(z)\prod_{i=1}^{C-1}B_{i}(z)$ gives

$\displaystyle\log Q(z)=\log A(z)+\sum_{i=1}^{C-1}\log B_{i}(z)$ (88)

Taking derivatives of both sides leads to

$\displaystyle\frac{Q^{\prime}(z)}{Q(z)}=\frac{A^{\prime}(z)}{A(z)}+\sum_{i=1}^{C-1}\frac{B^{\prime}_{i}(z)}{B_{i}(z)}$
(89)

Taking derivatives again:

$\displaystyle\frac{Q^{\prime\prime}(z)}{Q(z)}-\frac{Q^{\prime}(z)^{2}}{Q(z)^{2}}=\frac{A^{\prime\prime}(z)}{A(z)}-\frac{A^{\prime}(z)^{2}}{A(z)^{2}}+\sum_{i=1}^{C-1}\left(\frac{B^{\prime\prime}_{i}(z)}{B_{i}(z)}-\frac{B^{\prime}_{i}(z)^{2}}{B_{i}(z)^{2}}\right)$
(90)

Solving for $Q^{\prime\prime}(z)$ and letting $z=1$ gives:

$\displaystyle
Q^{\prime\prime}(1)=\mathbb{E}[Q^{(n)}]^{2}+A^{\prime\prime}(1)-A^{\prime}(1)^{2}+\sum_{i=1}^{C-1}\left(B^{\prime\prime}_{i}(1)-B^{\prime}_{i}(1)^{2}\right)$
(91)

Notice that $B^{\prime\prime}_{i}(1)=0$ ($\forall i=1,...,C-1$) and
$\mathbb{E}[Q^{(n)}]=Q^{\prime}(1)$. Substituting Eq. 77 and 91 into Eq. 24
gives

$\displaystyle\text{Var}[Q^{(n)}]=A^{\prime\prime}(1)-A^{\prime}(1)^{2}+A^{\prime}(1)+\sum_{i=1}^{C-1}\left(B_{i}^{\prime}(1)-B^{\prime}_{i}(1)^{2}\right)$
(92)

Now we only need to solve for $A^{\prime\prime}(1)$. The process is similar to
finding $A^{\prime}(1)$. Applying L’Hopital’s rule five times to Eq. 80 and
substituting $z=1$ leads to

$\displaystyle
A^{\prime\prime}(1)=\frac{-2A_{2}^{\prime}(1)A_{2}^{\prime\prime\prime}(1)+3A_{2}^{\prime\prime}(1)^{2}}{6A_{2}^{\prime}(1)}$
(93)

Notice that the derivation process uses $A_{1}^{\prime\prime}(z)=0$,
$A_{1}(1)=0$, $A_{2}(1)=0$, and $A_{1}^{\prime}(1)=A_{2}^{\prime}(1)$. Details
are omitted due to the tedious mathematical manipulation. To obtain
$A_{2}^{\prime\prime\prime}(1)$, taking derivative of Eq. 84 gives:

$\displaystyle A_{2}^{\prime\prime\prime}(1)$
$\displaystyle=\left[\frac{C{(C-1)(C-2)}z^{C-3}}{Y(z)}\right.-\frac{3Y^{\prime}(z)C(C-1)z^{C-2}}{Y(z)^{2}}-\frac{3Y^{\prime\prime}(z)Cz^{C-1}}{Y(z)^{2}}+\frac{4Y^{\prime}(z)^{2}Cz^{C-1}Y(z)}{Y(z)^{4}}$
$\displaystyle-\frac{Y^{\prime\prime\prime}(z)z^{C}}{Y(z)^{2}}+\frac{2Y(z)Y^{\prime}(z)Y^{\prime\prime}(z)z^{C}}{Y(z)^{4}}+\frac{4Y^{\prime\prime}(z)Y^{\prime}(z)z^{C}+2Cz^{C-1}Y^{\prime}(z)^{2}}{Y(z)^{3}}-\frac{6Y(z)Y^{\prime}(z)^{3}z^{C}}{Y(z)^{6}}$
$\displaystyle-\left.\sum_{u=0}^{C}(C-u)(C-u-1)(C-u-2)s_{u}^{(n)}z^{C-u-3}\right]_{z=1}$
$\displaystyle=C(C-1)(C-2)-3\bar{Y}^{(n)}C(C-1)-3{Y}^{\prime\prime}(1)C+6(\bar{Y}^{(n)})^{2}C-Y^{\prime\prime\prime}(1)+6\bar{Y}^{(n)}Y^{\prime\prime}(1)$
$\displaystyle-6(\bar{Y}^{(n)})^{3}-(C^{3}-3C^{2}+2C)+(2+3C^{2}-6C)\bar{S}^{(n)}+(3-3C)\mathbb{E}[(S^{(n)})^{2}]+\mathbb{E}[(S^{(n)})^{3}]$
(94)

Notice that
$Y^{\prime\prime\prime}(1)=\mathbb{E}[(Y^{(n)})^{3}]-3\mathbb{E}[(Y^{(n)})^{2}]+2\bar{Y}^{(n)}$.
Hence,

$\displaystyle A_{2}^{\prime\prime\prime}(1)$
$\displaystyle=3C^{2}\bar{S}^{(n)}-3C^{2}\bar{Y}^{(n)}-6C\bar{S}^{(n)}-3C\mathbb{E}[(S^{(n)})^{2}]-3C\mathbb{E}[(Y^{(n)})^{2}]+6C(\bar{Y}^{(n)})^{2}+6C\bar{Y}^{(n)}+2\bar{S}^{(n)}$
$\displaystyle+3\mathbb{E}[(S^{(n)})^{2}]+\mathbb{E}[(S^{(n)})^{3}]+6\mathbb{E}[(Y^{(n)})^{2}]\bar{Y}^{(n)}+3\mathbb{E}[(Y^{(n)})^{2}]-\mathbb{E}[(Y^{(n)})^{3}]-6(\bar{Y}^{(n)})^{3}-6(\bar{Y}^{(n)})^{2}-2\bar{Y}^{(n)}$
(95)

Substituting Eq. 83, 84, and 95 into Eq. 93 results in

$\displaystyle A^{\prime\prime}(1)$
$\displaystyle=\left[6C^{2}(\bar{S}^{(n)})^{2}-12C^{2}\right.(\bar{S}^{(n)})(\bar{Y}^{(n)})+6C^{2}(\bar{Y}^{(n)})^{2}-6C(\bar{S}^{(n)})\mathbb{E}[(S^{(n)})^{2}]-6C(\bar{S}^{(n)})\mathbb{E}[(Y^{(n)})^{2}](\bar{Y}^{(n)})^{2}$
$\displaystyle+12C(\bar{S}^{(n)})+6C\mathbb{E}[(S^{(n)})^{2}](\bar{Y}^{(n)})+6C\mathbb{E}[(Y^{(n)})^{2}](\bar{Y}^{(n)})-12C(\bar{Y}^{(n)})^{3}-(\bar{S}^{(n)})^{2}-2(\bar{S}^{(n)})\mathbb{E}[(S^{(n)})^{3}]$
$\displaystyle-12(\bar{S}^{(n)})\mathbb{E}[(Y^{(n)})^{2}](\bar{Y}^{(n)})+2(\bar{S}^{(n)})\mathbb{E}[(Y^{(n)})^{3}]+12(\bar{S}^{(n)})(\bar{Y}^{(n)})^{3}+2(\bar{S}^{(n)})(\bar{Y}^{(n)})$
$\displaystyle+3\mathbb{E}[(S^{(n)})^{2}]^{2}+6\mathbb{E}[(S^{(n)})^{2}]\mathbb{E}[(Y^{(n)})^{2}]-12\mathbb{E}[(S^{(n)})^{2}](\bar{Y}^{(n)})^{2}+2\mathbb{E}[(S^{(n)})^{3}](\bar{Y}^{(n)})+3\mathbb{E}[(Y^{(n)})^{2}]^{2}$
$\displaystyle-2\mathbb{E}[(Y^{(n)})^{3}](\bar{Y}^{(n)})-\left.(\bar{Y}^{(n)})^{2}\right]\bigg{/}6(\bar{S}^{(n)}-\bar{Y}^{(n)})$
(96)

Now with Eq. 96 and 85 we have

$\displaystyle
A^{\prime\prime}(1)-A^{\prime}(1)^{2}+A^{\prime}(1)=\left[(\bar{S}^{(n)})^{2}\right.-4\bar{S}^{(n)}\mathbb{E}[(S^{(n)})^{3}]-24\bar{S}^{(n)}\mathbb{E}[(Y^{(n)})^{2}]\bar{Y}^{(n)}+4\bar{S}^{(n)}\mathbb{E}[(Y^{(n)})^{3}]$
$\displaystyle+24\bar{S}^{(n)}(\bar{Y}^{(n)})^{3}-2\bar{S}^{(n)}\bar{Y}^{(n)}+3\mathbb{E}[(S^{(n)})^{2}]^{2}+6\mathbb{E}[(S^{(n)})^{2}]\mathbb{E}[(Y^{(n)})^{2}]-12\mathbb{E}[(S^{(n)})^{2}](\bar{Y}^{(n)})^{2}+4\mathbb{E}[(S^{(n)})^{3}]\bar{Y}^{(n)}$
$\displaystyle+3\mathbb{E}[(Y^{(n)})^{2}]^{2}+12\mathbb{E}[(Y^{(n)})^{2}](\bar{Y}^{(n)})^{2}-4\mathbb{E}[(Y^{(n)})^{3}]\bar{Y}^{(n)}-12(\bar{Y}^{(n)})^{4}+\left.(\bar{Y}^{(n)})^{2}\right]\bigg{/}12(\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}$
(97)

Slight manipulation of Eq. 97 leads to

$\displaystyle A^{\prime\prime}(1)-A^{\prime}(1)^{2}+A^{\prime}(1)$
$\displaystyle=\frac{-4(\bar{\bar{\bar{S}}}^{(n)}-\bar{\bar{\bar{Y}}}^{(n)})(\bar{S}^{(n)}-\bar{Y}^{(n)})+3(\bar{\bar{S}}^{(n)}+\bar{\bar{Y}}^{(n)})^{2}-[6({\bar{\bar{S}}}^{(n)}-{\bar{\bar{Y}}}^{(n)})-1](\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}-(\bar{S}^{(n)}-\bar{Y}^{(n)})^{4}}{12(\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}}$
(98)

Observe that
$B_{i}^{\prime}(1)-B_{i}^{\prime}(1)^{2}=\frac{-z_{i}^{*}}{(1-z_{i}^{*})^{2}}$.
Therefore, substituting Eq. 97 into 92 gives the final results:

$\displaystyle\text{Var}[Q^{(n)}]$
$\displaystyle=\frac{1}{12(\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}}\left\\{-4(\bar{\bar{\bar{S}}}^{(n)}-\bar{\bar{\bar{Y}}}^{(n)})(\bar{S}^{(n)}-\bar{Y}^{(n)})+3(\bar{\bar{S}}^{(n)}+\bar{\bar{Y}}^{(n)})^{2}\right.-$
$\displaystyle\left.[6({\bar{\bar{S}}}^{(n)}-{\bar{\bar{Y}}}^{(n)})-1](\bar{S}^{(n)}-\bar{Y}^{(n)})^{2}-(\bar{S}^{(n)}-\bar{Y}^{(n)})^{4}]\right\\}-\sum_{i=1}^{C-1}\frac{z_{i}^{*}}{(1-z_{i}^{*})^{2}}$
(99)

##  12 Proof of Proposition 4.8

When there are no incidents in the system, vehicle $l$ reaches station $n$
after $T^{(n)}$ time units. Since the system can only switch to the incident
state from the normal state, the number of incident occurrences, $K$, follows
a Poisson distribution with rate $\gamma T^{(n)}$. The vehicle stopping time
for the $i$-th incident, $X_{i}$, follows an exponential distribution with
rate $\theta$ (i.e., mean $\frac{1}{\theta}$). Therefore, the duration of all
incidents is $I^{(n,l)}=\sum_{i=1}^{K}X_{i}$, where
$X_{i}\sim\text{{Exp}}(\theta)\;\;\forall i=1,...,K$, and
$K\sim\text{{Poi}}(\gamma T^{(n)})$

##  13 Proof of Proposition 4.9

$\displaystyle M_{H^{(n,l)}}(t)$
$\displaystyle=\mathbb{E}[e^{tH^{(n,l)}}]=\mathbb{E}[e^{t(\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}}}e^{tI^{(n,l)}}e^{-tI^{(n,l-1)}}]$
$\displaystyle=e^{t(\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}}}\mathbb{E}[e^{tI^{(n,l)}}]\mathbb{E}[e^{-tI^{(n,l-1)}}]$
$\displaystyle=e^{t(\bar{H}+\frac{2\gamma T^{(N)}}{\theta\bar{F}}}e^{\gamma
T^{(n)}(\frac{\theta}{\theta-t}-1)}e^{\gamma
T^{(n)}(\frac{\theta}{\theta+t}-1)}$ $\displaystyle=e^{t(\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}})}e^{\gamma T^{(n)}(\frac{2t^{2}}{\theta^{2}-t^{2}})}$
(100)

where Eq. 100 is because of the independence between $I^{(n,l)}$ and
$I^{(n,l-1)}$. As this equation holds for all vehicles $l$, the MGF of
$H^{(n)}$ (i.e., $l\to\infty$) is $M_{H^{(n)}}(t)=M_{H^{(n,l)}}(t)$.

##  14 Proof of Proposition 4.11

Let $\mu$ and $\sigma^{2}$ be the mean and variance of
${H}^{(n)}_{\text{Normal}}$, respectively, where $\mu=\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}}$ and $\sigma=\frac{2\sqrt{T^{(n)}\gamma}}{\theta}$.
The MGF of $\hat{H}^{(n)}_{\text{Normal}}$ can be derived as

$\displaystyle M_{\hat{H}^{(n)}_{\text{Normal}}}(t)$
$\displaystyle=\mathbb{E}[e^{t\hat{H}^{(n)}_{\text{Normal}}}]=\mathbb{P}[H^{(n)}_{\text{Normal}}\leq
0]\cdot
e^{0}+\int_{0}^{+\infty}e^{tz}\cdot\phi_{{H}^{(n)}_{\text{Normal}}}{(z)}\cdot
dz$
$\displaystyle=\Phi(\frac{-\mu}{\sigma})+\frac{1}{\sigma\sqrt{2\pi}}\int_{0}^{+\infty}e^{tz+\frac{(z-\mu)^{2}}{-2\sigma^{2}}}dz$
$\displaystyle=\Phi(\frac{-\mu}{\sigma})+e^{\mu
t+\frac{\sigma^{2}t^{2}}{2}}\left[1-\Phi(\frac{-\mu}{\sigma}-\sigma t)\right]$
(101)

where Eq. 101 follows the same derivation of a truncated normal distribution
(Burkardt 2014). Subsisting the value of $\mu$ and $\sigma$ completes the
proof.

##  15 Proof of Proposition 4.12

The strict mathematical proof can be done by taking the derivative of
$\mathbb{E}[\hat{H}^{(n)}_{\text{Normal}}]$ in terms of $\gamma$ or
$\frac{1}{\theta}$ and show that it is always positive. However, in this
study, we adopt a more intuitive graphical proof, which is easier for
understanding.

As shown in Figure 10, consider an arbitrary truncated headway distribution
(shown in the red line, denoted the headway as $\hat{H}_{\text{Red}}$). When
the incident intensity increases, according to Eqs. 36 and 37, both $\mu$ and
$\sigma$ increase. Let us first consider the increase in $\sigma$ and assume
$\mu$ does not change (which corresponds to the scenario where
$\bar{F}\to\infty$). Then the distribution will become the blue curve (denote
the corresponding headway as $\hat{H}_{\text{Blue}}$). Note that
$\hat{H}_{\text{Red}}$ and $\hat{H}_{\text{Blue}}$ have the same peak value,
but since $\hat{H}_{\text{Blue}}$ has longer positive tail, we have
$\mathbb{E}[\hat{H}_{\text{Blue}}]>\mathbb{E}[\hat{H}_{\text{Red}}]$. Next,
let us also consider the incident’s impact on the increase in $\mu$ as well.
The distribution is shown by the green curve (denoted the headway as
$\hat{H}_{\text{Green}}$). Since $\hat{H}_{\text{Blue}}$ and
$\hat{H}_{\text{Green}}$ has the same $\sigma$, but $\hat{H}_{\text{Green}}$
has higher $\mu$ (shifted right), we have
$\mathbb{E}[\hat{H}_{\text{Green}}]>\mathbb{E}[\hat{H}_{\text{Blue}}]$. Hence,
$\mathbb{E}[\hat{H}_{\text{Green}}]>\mathbb{E}[\hat{H}_{\text{Red}}]$, showing
that the increase in incident intensity will increase $\mu$ and $\sigma$, thus
increase the mean of the truncated headway.

(a)  Figure 10: Illustration for the impact of incidents on expected headway.
As the probability mass at zero does not contribute to the expectation
calculation, it is not shown in the figure.

##  16 Approximating headway as a normal distributed random variable

We observe that the third central moment of $H^{(n)}$, which is a measure of
skewness, is $\text{Skewness}[H^{(n)}]=0$, implying that $H^{(n)}$ is
symmetric. Moreover, the MGF of the normal distribution of
$H^{(n)}_{\text{Normal}}$ with the same mean and variance is

$\displaystyle M_{H^{(n)}_{\text{Normal}}}(t)=e^{t(\bar{H}+\frac{2\gamma
T^{(N)}}{\theta\bar{F}})}e^{\gamma T^{(n)}(\frac{2t^{2}}{\theta^{2}})},$ (102)

which is very similar to Eq. 35 (the MGF of $H^{(n)}$). Therefore, it is
reasonable to approximate the distribution of $H^{(n)}$ as a normal
distribution with the same mean and variance. Note that the first three
moments of $H^{(n)}$ and $H^{(n)}_{\text{Normal}}$ are the same. And the
corresponding forth moments (i.e., Kurtosis) are:

$\displaystyle\text{Kurtosis}[H^{(n)}]$
$\displaystyle=\frac{48(T^{(n)}\gamma)^{2}+48T^{(n)}\gamma}{\theta^{4}}$ (103)
$\displaystyle\text{Kurtosis}[H^{(n)}_{\text{Normal}}]$
$\displaystyle=\frac{48(T^{(n)}\gamma)^{2}}{\theta^{4}}$ (104)

which means that the distribution of $H^{(n)}$ may have heavier tails and
peakedness compared to $H^{(n)}_{\text{Normal}}$.

Figure 11 empirically compares the distribution of $H^{(n)}$ and
$H^{(n)}_{\text{Normal}}$ with various values of $T^{(n)},\theta$, and
$\gamma$. The histogram of $H^{(n)}$ is generated by sampling variables from
the associated exponential and Poisson distributions to get the compound
distribution. Results show that the normal distribution approximates the
original distribution well. As expected, $H^{(n)}$ shows more peakedness than
$H^{(n)}_{\text{Normal}}$.

(a) Example 1

(b) Example 2

(c) Example 3

Figure 11: Empirical validation for approximating the headway distribution as
normal

##  17 Proof of Proposition 4.14

Recall that $Y^{(n)}|\hat{H}^{(n)}_{\text{Normal}}$ is a Poisson random
variable with parameter $\lambda^{(n)}\hat{H}^{(n)}_{\text{Normal}}$. So, the
MGF of $Y^{(n)}|\hat{H}^{(n)}_{\text{Normal}}$ is
$M_{Y^{(n)}|\hat{H}^{(n)}_{\text{Normal}}}(t)=\exp[{\lambda^{(n)}\hat{H}^{(n)}_{\text{Normal}}(e^{t}-1)}]$.
Based on Lemma 4.13, setting $C_{1}(t)=1$ and
$C_{2}(t)=\lambda^{(n)}(e^{t}-1)$, we conclude that the MGF of $Y^{(n)}$ is

$\displaystyle
M_{Y^{(n)}}(t)=\Phi\left(\frac{-\mu}{\sigma}\right)+e^{\mu\lambda^{(n)}(e^{t}-1)+\frac{\sigma^{2}(\lambda^{(n)}e^{t}-\lambda^{(n)})^{2}}{2}}\left[1-\Phi\left(\frac{-\mu}{\sigma}-\sigma\lambda^{(n)}(e^{t}-1)\right)\right]$
(105)

Substituting $t=\log z$ in Eq. 105 completes the proof.

##  18 Interpolation-based searching algorithm for root solving

Notice that $\textsc{Den}(z)=0$ is equivalent to find
$z_{0}^{*},...,z_{C-1}^{*}$, such that

$\displaystyle\frac{1}{Y(z_{k}^{*})}-S(1/z_{k}^{*})=0\Leftrightarrow
J(z_{k}^{*})=1\quad\quad\forall\;k=0,...,C-1$ (106)

where $J(z):=Y(z)S(1/z)$. Taking the logarithm of both sides of Eq. 106 and
matching the real and imaginary parts gives:

$\left\\{\begin{aligned} &\text{Re}[\log(J(z))]=0\\\
&\text{Im}[\log(J(z))]=0\end{aligned}\right.$ (107)

where $\text{Re}[\cdot]$ and $\text{Im}[\cdot]$ represent the real and
imaginary part of a complex number. Eq. 107 can be solved efficiently with
many optimization algorithms (such as trust-region and Levenberg-Marquardt
algorithms). However, as there are $C$ optimal solutions for this problem with
$|z^{*}|\leq 1$, the challenge is how to select different initial values so as
to find all solutions.

It can be empirically observed that the distribution of the $C$ solutions has
an oval-like shape. Figure 12 shows some examples of the solution distribution
with different values of $\rho^{(n)}$ (where
$\rho^{(n)}={\bar{Y}^{(n)}}/{\bar{S}^{(n)}}$ is the utilization ratio of a
bulk service queuing system) and $s^{(n)}$. It is found that the closer
$\rho^{(n)}$ is to 1 (resp. 0), the closer the shape of the root distribution
is to an ellipse (resp. circle). The value of $s^{(n)}$ (i.e., available
capacity distribution) can also slightly affect the root distribution.

(a) $\rho^{(n)}=0$, $C=40$

(b) $\rho^{(n)}=0.15$, $C=40$

(c) $\rho^{(n)}=0.64$, $C=40$

(d) $\rho^{(n)}=0.74$, $C=40$

Figure 12: Examples of root distribution

We first express the complex number in polar coordinate system with
$z=r\exp[\varphi i]$, where $i=\sqrt{-1}$, $r$ is the length from $z$ to the
origin, and $\varphi$ is the angle. Eq. 107 now has $C$ optimal solutions
$(r_{k}^{*},\varphi_{k}^{*})$ for $k=0,1,...,C-1$, where $0\leq r_{k}^{*}\leq
1$ and $0\leq\varphi_{k}^{*}<2\pi$. Note that $z_{0}^{*}=1$ corresponds to
$r_{0}^{*}=1$ and $\varphi_{0}^{*}=0$. Another property is that the roots must
appear as conjugate pairs. Hence, if $(r^{*},\varphi^{*})$ is a root and
$0<\varphi^{*}<\pi$, then $(r^{*},2\pi-\varphi^{*})$ is also a root.

The proposed search algorithm has two steps. The first step is referred to as
“clockwise searching”, which is adapted from the numerical method in Powell
(1985). The empirical observation (Figure 12) shows a rough relationship that
$r_{k+1}^{*}-r_{k}^{*}\approx r_{k}^{*}-r_{k-1}^{*}$, especially for small
$\rho^{(n)}$. This is equivalent to

$\displaystyle r_{k+1}^{*}\approx 2r_{k}^{*}-r_{k-1}^{*}$ (108)

Eq. 108 provides a way to determine the initial value for solving for the
$k+1$-th root when the $k$-th and $k-1$-th roots are available. As we already
know $r_{0}^{*}=1$ and $\varphi_{0}^{*}=0$, we first set the initial value for
solving for the second root as $r_{1}^{\text{Ini}}=1-0.5\rho^{(n)}$ and
$\varphi_{1}^{\text{Ini}}=3\pi/C$. This is motivated by the shape of the root
distribution with respect to $\rho^{(n)}$. Then $r_{1}^{\text{Ini}}$ and
$\varphi_{1}^{\text{Ini}}$ are used as the initial value to solve for
$r_{1}^{*}$ and $\varphi_{1}^{*}$ based on Eq 107. For $k\geq 2$, the initial
values for solving the for $k$-th root are set to
$r_{k}^{\text{Ini}}=r_{k-1}^{*}+(r_{k-1}^{*}-r_{k-2}^{*})$,
$\varphi_{k}^{\text{Ini}}=\varphi_{k-1}^{*}+(\varphi_{k-1}^{*}-\varphi_{k-2}^{*})$
according to Eq. 108.

However, only performing step 1 (i.e., Powell (1985)’s method) may not find
all $C$ distinct roots. Figure 13 shows some examples of the comparison
between roots found in step 1 and all roots. We observe that when $\rho^{(n)}$
is relatively large (i.e., the system is relatively congested), the clockwise
search does not perform well because the approximate relationship in Eq. 108
does not hold. Even when $\rho^{(n)}$ is relatively small, it is also possible
that some roots do not perfectly fit the oval-like shape (such as Figure 13a),
resulting in the failure of step 1 to find all roots.

(a) $\rho^{(n)}=0.17$, $C=40$

(b) $\rho^{(n)}=0.48$, $C=36$

(c) $\rho^{(n)}=0.64$, $C=40$

(d) $\rho^{(n)}=0.85$, $C=40$

Figure 13: Comparison between roots found with clockwise search and all roots

Therefore, we propose a second step called “interpolation search”. Let the set
of found roots from step 1 be
$\mathcal{Z}^{(0)}=\\{(r_{0}^{(0)},\varphi_{0}^{(0)}),(r_{1}^{(0)},\varphi_{1}^{(0)}),...,(r_{M_{0}}^{(0)},\varphi_{M_{0}}^{(0)})\\}$,
where $M_{0}=|\mathcal{Z}^{(0)}|$ is the number of roots from step 1. Without
loss of generality, assume that the elements in $\mathcal{Z}^{(0)}$ are
clockwise ranked (i.e.,
$\varphi_{0}^{(0)}<\varphi_{1}^{(0)}<...<\varphi_{M_{0}}^{(0)}$). The
interpolation search is described in Algorithm 2. The main idea is to perform
interpolation between any two adjacent roots that are already found. The
interpolated points are set as initial values and fed into Eq. 107 to solve
for new distinct roots. Then we update the set of roots with the new distinct
roots or perform a finer (i.e., larger $L$) interpolation if no distinct roots
are found. This process is repeated until there are $C$ distinct roots found.
In Algorithm 2, $L$ is a parameter controlling how many points to interpolate
between two known roots, and $\epsilon$ is a predetermined probability
threshold to add randomness in the search process.

Algorithm 2 Interpolation searching

1:Initialize $\mathcal{Z}^{(0)}$, $M_{0}$, $\epsilon$. Initialize $L=2$,
$k=0$.

2:while $M_{k}<C$ do

3: Initialize $\mathcal{Z}^{\text{Ini}}$ as an empty set.

4: for $i=1:M_{k}$ do

5: for $d=1:L-1$ do

6: $r^{\text{Ini}}=r_{i}^{(k)}+d\cdot\frac{r_{i+1}^{(k)}-r_{i}^{(k)}}{L}$;
$\varphi^{\text{Ini}}=\varphi_{i}^{(k)}+d\cdot\frac{\varphi_{i+1}^{(k)}-\varphi_{i}^{(k)}}{L}$

7: Draw a random value $w$ uniformly from $[0,1)$

8: if $w<\epsilon$ then

9: Draw a random value $\delta_{1}$ uniformly from
$[-\frac{|r_{i+1}^{(k)}-r_{i}^{(k)}|}{2L},\frac{|r_{i+1}^{(k)}-r_{i}^{(k)}|}{2L}]$

10: $r^{\text{Ini}}=r^{\text{Ini}}+\delta_{1}$

11: Draw a random value $\delta_{2}$ uniformly from
$[-\frac{|\varphi_{i+1}^{(k)}-\varphi_{i}^{(k)}|}{2L},\frac{|\varphi_{i+1}^{(k)}-\varphi_{i}^{(k)}|}{2L}]$

12: $\varphi^{\text{Ini}}=\varphi^{\text{Ini}}+\delta_{2}$

13: Add $(r^{\text{Ini}},\varphi^{\text{Ini}})$ into
$\mathcal{Z}^{\text{Ini}}$.

14: Initialize $\mathcal{Z}^{\text{temp}}$ as an empty set.

15: for all $z^{\text{Ini}}$ in $\mathcal{Z}^{\text{Ini}}$ do

16: Solve Eq. 107 using $z^{\text{Ini}}$ as the initial value, obtaining
$z^{*}_{\text{temp}}$. Let its conjugate be $\bar{z}^{*}_{\text{temp}}$.

17: If $z^{*}_{\text{temp}}$ ($\bar{z}^{*}_{\text{temp}}$) not in
$\mathcal{Z}^{(k)}$, add it to $\mathcal{Z}^{\text{temp}}$, otherwise do
nothing.

18: $\mathcal{Z}^{(k+1)}=\mathcal{Z}^{(k)}\cup\mathcal{Z}^{\text{temp}}$ and
rank all elements in $\mathcal{Z}^{(k+1)}$ clockwise

19: Denote $\mathcal{Z}^{(k+1)}$ as
$\\{(r_{0}^{({k+1})},\varphi_{0}^{({k+1})}),...,(r_{M_{k+1}}^{({k+1})},\varphi_{M_{k+1}}^{({k+1})})\\}$

20: $k=k+1$

21: if $M_{k+1}$ = $M_{k}$ then

22: $L=L+1$

##  19 Proof of Proposition 4.15

The stability condition is equivalent to
$\mathbb{P}(Q^{(n)}=0)=q_{0}^{(n)}>0$. In Eq. 20, we notice that
$\prod_{i=1}^{C-1}\frac{z_{i}^{*}}{z_{i}^{*}-1}$ is always greater than 0 (see
Appendix 10 for details), and $s_{C}^{(n)}>0$ is a known condition. Therefore,
$q_{0}^{(n)}>0$ if and only if $\bar{Y}^{(n)}<\bar{S}^{(n)}$ (i.e.,
$\rho^{(n)}<1$), which completes the proof.

##  20 Simulation procedure for comparison

For each vehicle $l$ at each station $n$, we generate the total duration of
incidents $I^{(n,l)}$ as a compound Poisson exponential variable to get the
arrival time. Since no overtaking is allowed, the arrival time at station $n$
cannot be earlier than vehicle $l-1$. When a vehicle arrives at a station,
passengers board based on the first-come-first-serve (FCFS) principle up to
the vehicle’s capacity $C$. Queue lengths at vehicle arrival and passenger
waiting times are recorded during the simulation. To ensure the system reaches
steady-state conditions, the first 10% records are dropped.

Algorithm 3 Simulation procedure

1:Initialize model parameters: $C$, $\gamma$, $\theta$, $\bar{H}$, Demand
factor. Set the total number of vehicles $L$.

2:for $l$ = 1:$L$ do

3: Get vehicle dispatch time as $DT^{(l)}$

4: for $n$ = $1:N$ do

5: Sample total incident duration $I^{(n,l)}$ from a compound Poisson
exponential distribution

6: $t^{(n,l)}_{D}=\min{\\{DT^{(l)}+T^{(n)}+I^{(n,l)},t^{(n,l-1)}_{D}\\}}$

7: Headway for vehicle $l$ at station $n$ is $t^{(n,l)}_{D}-t^{(n,l-1)}_{D}$

8: Sample the arrival passengers within the headway as a Poisson process based
on $\lambda^{(n)}$.

9: Record queue length (including left behind passengers from the last run)

10: Alight passengers based on the binomial distribution with parameter
$\alpha^{(n)}$

11: Board passengers based on FCFS principle up to the vehicle capacity

12: Record left behind passengers and passengers’ waiting time

13:Drop the first 10% records. Calculate
$\mathbb{E}[Q^{(n)}],\text{Var}[Q^{(n)}],\mathbb{E}[W^{(n)}]$, and
$\text{Var}[W^{(n)}]$ based on the recorded samples for $n=1,...,N$

## References

  * Andersson and Scalia-Tomba (1981) Andersson PÅ, Scalia-Tomba GP (1981) A mathematical model of an urban bus route. _Transportation Research Part B: Methodological_ 15(4):249–266.
  * Ayyappan and Karpagam (2020) Ayyappan G, Karpagam S (2020) Analysis of a bulk service queue with unreliable server, multiple vacation, overloading and stand-by server. _International Journal of Mathematics in Operational Research_ 16(3):291–315.
  * Bailey (1954) Bailey NT (1954) On queueing processes with bulk service. _Journal of the Royal Statistical Society: Series B (Methodological)_ 16(1):80–87.
  * Balakrishna et al. (2008) Balakrishna R, Wen Y, Ben-Akiva M, Antoniou C (2008) Simulation-based framework for transportation network management in emergencies. _Transportation Research Record_ 2041(1):80–88.
  * Beardon (2019) Beardon AF (2019) _Complex analysis: The argument principle in analysis and topology_ (Courier Dover Publications).
  * Bellei and Gkoumas (2010) Bellei G, Gkoumas K (2010) Transit vehicles’ headway distribution and service irregularity. _Public transport_ 2(4):269–289.
  * Berdica (2002) Berdica K (2002) An introduction to road vulnerability: what has been done, is done and should be done. _Transport policy_ 9(2):117–127.
  * Burkardt (2014) Burkardt J (2014) The truncated normal distribution. _Department of Scientific Computing Website, Florida State University_ 1–35.
  * Chaudhry and Templeton (1983) Chaudhry M, Templeton JG (1983) First course in bulk queues. .
  * Chaudhry et al. (1987) Chaudhry ML, Madill B, Briere G (1987) Computational analysis of steady-state probabilities of m/g a, b/1 and related nonbulk queues. _Queueing systems_ 2(2):93–114.
  * Clerq (1972) Clerq F (1972) A public transport assignment method. _Verkeerstechniek, Netherlands_ 23(6).
  * Currie and Muir (2017) Currie G, Muir C (2017) Understanding passenger perceptions and behaviors during unplanned rail disruptions. _Transportation Research Procedia_ 25:4396–4406.
  * Daganzo (2009) Daganzo CF (2009) A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons. _Transportation Research Part B: Methodological_ 43(10):913–921.
  * Dial (1967) Dial RB (1967) Transit path finder algorithm. _Highway Research Record_ (205).
  * Downton (1955) Downton F (1955) Waiting time in bulk service queues. _Journal of the Royal Statistical Society: Series B (Methodological)_ 17(2):256–261.
  * Fukasawa et al. (2012) Fukasawa N, Yamauchi K, Murakoshi A, Fujinami K, Tatsui D (2012) Provision of forecast train information and consequential impact on decision making for train-choice. _Quarterly Report of RTRI_ 53(3):141–147.
  * Hickman (2001) Hickman MD (2001) An analytic stochastic model for the transit vehicle holding problem. _Transportation Science_ 35(3):215–237.
  * Hogg et al. (2010) Hogg RV, Tanis EA, Zimmerman DL (2010) _Probability and statistical inference_ (Pearson/Prentice Hall Upper Saddle River, NJ, USA:).
  * Hong et al. (2018) Hong L, Gao J, Zhu W (2018) Self-evacuation modelling and simulation of passengers in metro stations. _Safety science_ 110:127–133.
  * Islam et al. (2014) Islam MK, Vandebona U, Dixit VV, Sharma A (2014) A bulk queue model for the evaluation of impact of headway variations and passenger waiting behavior on public transit performance. _IEEE Transactions on Intelligent Transportation Systems_ 15(6):2432–2442.
  * Islam et al. (2015) Islam MK, Vandebona U, Dixit VV, Sharma A (2015) A model to evaluate the impact of headway variation and vehicle size on the reliability of public transit. _IEEE Transactions on Intelligent Transportation Systems_ 16(4):1840–1850.
  * Jaiswal (1960) Jaiswal N (1960) Time-dependent solution of the bulk-service queuing problem. _Operations Research_ 8(6):773–781.
  * Jayaraman et al. (1994) Jayaraman D, Nadarajan R, Sitrarasu M (1994) A general bulk service queue with arrival rate dependent on server breakdowns. _Applied mathematical modelling_ 18(3):156–160.
  * Jin et al. (2014) Jin JG, Tang LC, Sun L, Lee DH (2014) Enhancing metro network resilience via localized integration with bus services. _Transportation Research Part E: Logistics and Transportation Review_ 63:17–30.
  * Jin et al. (2016) Jin JG, Teo KM, Odoni AR (2016) Optimizing bus bridging services in response to disruptions of urban transit rail networks. _Transportation Science_ 50(3):790–804.
  * Kahraman and Gosavi (2011) Kahraman A, Gosavi A (2011) On the distribution of the number stranded in bulk-arrival, bulk-service queues of the m/g/1 form. _European journal of operational research_ 212(2):352–360.
  * Krishnamoorthy et al. (2014) Krishnamoorthy A, Pramod PK, Chakravarthy SR (2014) Queues with interruptions: a survey. _Top_ 22(1):290–320.
  * Kroon and Huisman (2011) Kroon L, Huisman D (2011) Algorithmic support for railway disruption management. _Transitions Towards Sustainable Mobility_ , 193–210 (Springer).
  * Lin et al. (2018) Lin T, Srikukenthiran S, Miller E, Shalaby A (2018) Subway user behaviour when affected by incidents in toronto (subwait) survey—a joint revealed preference and stated preference survey with a trip planner tool. _Canadian Journal of Civil Engineering_ 45(8):623–633.
  * Luo et al. (2019) Luo C, Li X, Zhou Y, Caunhye AM, Alibrandi U, Aydin NY, Ratti C, Eckhoff D, Bojic I (2019) Data-driven disruption response planning for a mass rapid transit system. _Smart Transportation Systems 2019_ , 205–213 (Springer).
  * Madan (1989) Madan K (1989) A single channel queue with bulk service subject to interruptions. _Microelectronics Reliability_ 29(5):813–818.
  * Madan (1992) Madan K (1992) A bulk queueing system with random failures and two phase repairs. _Microelectronics Reliability_ 32(5):669–677.
  * Marguier (1985) Marguier PHJ (1985) _Bus route performance evaluation under stochastic considerations_. Ph.D. thesis, Massachusetts Institute of Technology.
  * Mo et al. (2020) Mo B, Ma Z, Koutsopoulos HN, Zhao J (2020) Capacity-constrained network performance model for urban rail systems. _Transportation Research Record_ 2674(5):59–69.
  * Mo et al. (2022) Mo B, Von Franque MY, Koutsopoulos HN, Attanucci JP, Zhao J (2022) Impact of unplanned long-term service disruptions on urban public transit systems. _IEEE Open Journal of Intelligent Transportation Systems_ 3:551–569.
  * Murray-Tuite et al. (2014) Murray-Tuite P, Wernstedt K, Yin W (2014) Behavioral shifts after a fatal rapid transit accident: A multinomial logit model. _Transportation research part F: traffic psychology and behaviour_ 24:218–230.
  * Neuts (1967) Neuts MF (1967) A general class of bulk queues with poisson input. _The Annals of Mathematical Statistics_ 38(3):759–770.
  * Osuna and Newell (1972) Osuna E, Newell GF (1972) Control strategies for an idealized public transportation system. _Transportation Science_ 6(1):52–72.
  * O’Dell and Wilson (1999) O’Dell SW, Wilson NH (1999) Optimal real-time control strategies for rail transit operations during disruptions. _Computer-aided transit scheduling_ , 299–323 (Springer).
  * Powell (1981) Powell WB (1981) _Stochastic delays in transportation terminals: New results in the theory and application of bulk queues_. Ph.D. thesis, Massachusetts Institute of Technology.
  * Powell (1983) Powell WB (1983) Bulk service queues with deviations from departure schedules: The problem of correlated headways. _Transportation Research Part B: Methodological_ 17(3):221–232.
  * Powell (1985) Powell WB (1985) Analysis of vehicle holding and cancellation strategies in bulk arrival, bulk service queues. _Transportation Science_ 19(4):352–377.
  * Rahimi et al. (2019) Rahimi E, Shamshiripour A, Shabanpour R, Mohammadian A, Auld J (2019) Analysis of transit users’ waiting tolerance in response to unplanned service disruptions. _Transportation Research Part D: Transport and Environment_ 77:639–653.
  * Rapoport et al. (2010) Rapoport A, Stein WE, Mak V, Zwick R, Seale DA (2010) Endogenous arrivals in batch queues with constant or variable capacity. _Transportation Research Part B: Methodological_ 44(10):1166–1185.
  * Rudin (2006) Rudin W (2006) _Real and complex analysis_ (Tata McGraw-hill education).
  * Sasikala and Indhira (2016) Sasikala S, Indhira K (2016) Bulk service queueing models–a survey. _International Journal of Pure and Applied Mathematics_ 106(6):43–56.
  * Selvi and Rosenshine (1983) Selvi A, Rosenshine M (1983) A queueing system for airport buses. _Transportation Research Part B: Methodological_ 17(6):427–434.
  * Singh and Ram (1991) Singh I, Ram C (1991) Three-server bulk service queue with service interruptions and exponential repairs. _Microelectronics Reliability_ 31(2-3):257–259.
  * Suarez et al. (2005) Suarez P, Anderson W, Mahal V, Lakshmanan T (2005) Impacts of flooding and climate change on urban transportation: A systemwide performance assessment of the boston metro area. _Transportation Research Part D: transport and environment_ 10(3):231–244.
  * Sun et al. (2016) Sun H, Wu J, Wu L, Yan X, Gao Z (2016) Estimating the influence of common disruptions on urban rail transit networks. _Transportation Research Part A: Policy and Practice_ 94:62–75.
  * Tadj and Choudhury (2009) Tadj L, Choudhury G (2009) A quorum queueing system with an unreliable server. _Applied mathematics letters_ 22(11):1710–1714.
  * Tadj et al. (2012) Tadj L, Choudhury G, Rekab K (2012) A two-phase quorum queueing system with bernoulli vacation schedule, setup, and n-policy for an unreliable server with delaying repair. _International Journal of Services and Operations Management_ 12(2):139–164.
  * Teng and Liu (2015) Teng J, Liu WR (2015) Development of a behavior-based passenger flow assignment model for urban rail transit in section interruption circumstance. _Urban Rail Transit_ 1(1):35–46.
  * Tian and Zheng (2018) Tian X, Zheng B (2018) Using smart card data to model commuters’ responses upon unexpected train delays. _2018 IEEE International Conference on Big Data (Big Data)_ , 831–840 (IEEE).
  * Villa and Escobar (2006) Villa ER, Escobar LA (2006) Using moment generating functions to derive mixture distributions. _The American Statistician_ 60(1):75–80.
  * Wang et al. (2014) Wang Y, Guo J, Ceder AA, Currie G, Dong W, Yuan H (2014) Waiting for public transport services: Queueing analysis with balking and reneging behaviors of impatient passengers. _Transportation Research Part B: Methodological_ 63:53–76.
  * Welding (1957) Welding P (1957) The instability of a close-interval service. _Journal of the operational research society_ 8(3):133–142.
  * Wilson (2014) Wilson H (2014) Complex zeros of functions. https://www.mathworks.com/matlabcentral/fileexchange/25178-complex-zeros-of-functions, mATLAB Central File Exchange. Retrieved December 6, 2020.
  * Wirasinghe (1980) Wirasinghe S (1980) Nearly optimal parameters for a rail/feeder-bus system on a rectangular grid. _Transportation Research Part A: General_ 14(1):33–40.
  * Xu et al. (2015) Xu X, Chen A, Jansuwan S, Heaslip K, Yang C (2015) Modeling transportation network redundancy. _Transportation research procedia_ 9:283–302.
  * Yin et al. (2016) Yin H, Han B, Li D, Wang Y (2016) Evaluating disruption in rail transit network: a case study of beijing subway. _Procedia Engineering_ 137:49–58.
  * Zhang et al. (2018) Zhang Dm, Du F, Huang H, Zhang F, Ayyub BM, Beer M (2018) Resiliency assessment of urban rail transit networks: Shanghai metro as an example. _Safety Science_ 106:230–243.