Mo et al.

Individual Path Recommendation

Individual Path Recommendation Under Public Transit Service Disruptions
Considering Behavior Uncertainty

Baichuan Mo Department of Civil and Environmental Engineering, Massachusetts
Institute of Technology, Cambridge, MA 02139 Haris N. Koutsopoulos Department
of Civil and Environmental Engineering, Northeastern University, Boston, MA
02115 Zuo-Jun Max Shen Department of Industrial Engineering and Operations
Research, University of California, Berkeley, Berkeley, CA 94720 Jinhua Zhao
Department of Urban Studies and Planning, Massachusetts Institute of
Technology, Cambridge, MA 20139

This study proposes a mixed-integer programming formulation to model the
individual-based path (IPR) recommendation problem during public transit
service disruptions with the objective of minimizing system travel time and
respecting passengers’ path choice preferences. Passengers’ behavior
uncertainty in path choices given recommendations is also considered. We model
the behavior uncertainty based on the passenger’s prior preferences and
posterior path choice probability distribution with two new concepts:
$\epsilon$-feasibility and $\Gamma$-concentration, which control the mean and
variance of path flows in the optimization problem. We show that these two
concepts can be seen as a way of approximating the recourse function (expected
system travel time) in a two-stage stochastic optimization. It is proved that
these two concepts help to bound the difference between the approximated
recourse function and the exact one. Additional theoretical analysis shows
that $\epsilon$-feasibility and $\Gamma$-concentration can be seen as an
approximation of expectation and chance constraints in a typical stochastic
optimization formulation, respectively. The proposed IPR problem with behavior
uncertainty is solved efficiently with Benders decomposition. The model is
implemented in the Chicago Transit Authority (CTA) system with a real-world
urban rail disruption as the case study. Results show that the proposed IPR
model significantly reduces the average travel times compared to the status
quo and outperforms the capacity-based benchmark path recommendation strategy.
We also show that incorporating behavior uncertainty with respect to responses
to information achieves lower system travel times than assuming that all
passengers would follow the recommendations. In terms of respecting people’s
preferences, results show that it is possible to make recommendations so that
most of the passengers (e.g., more than 70%) use their preferred paths while
only increasing the system travel time by 0.51%.

Individualized recommendation; Mixed-integer programming; Behavior
uncertainty.

##  1 Introduction

###  1.1 Background and challenges

With aging systems and near-capacity operations, service disruptions often
occur in urban public transit (PT) systems. These incidents may result in
passenger delays, cancellation of trips, and economic losses (Cox et al.
2011).

During a significant disruption where the service is interrupted for a
relatively long period of time (e.g., 1 hour), affected passengers usually
need to find an alternative path or use other travel modes (such as transfer
to another bus route). However, due to a lack of knowledge of the system
(especially during incidents), the routes chosen by passengers may not be
optimal or even cause more congestion (Mo et al. 2022c). For example, during a
rail disruption, most of the passengers may choose bus routes that are
parallel to the interrupted rail line as an alternative. However, given
limited bus capacity, parallel bus lines may become oversaturated and
passengers have to wait for a long time to board due to being denied boarding
(or left behind).

One of the strategies to better guide passengers is to provide path
recommendations so that passenger flows are re-distributed in a better way and
the system travel times are minimized. This can be seen as solving an optimal
passenger flow distribution (or assignment) problem over a public transit
network. However, different from the typical flow redistribution problem,
there are several unique characteristics and challenges for the path
recommendation problem under PT service disruptions.

  * •

Passengers may have different preferences on different alternative paths. This
heterogeneity suggests that we cannot treat a group of passengers simply as
flows. Individualization is needed in the path recommendation design.

  * •

Passengers may not follow the recommendation. When providing a specific path
recommendation to a passenger, their actual path choice is uncertain (though
the recommendation may change their preferences). This behavior uncertainty
brings challenges to the recommendation system design and has not been
considered in the path recommendation literature. In the context of
individualization, the behavior uncertainty is also individual-specific, which
requires a more granular modeling approach.

###  1.2 Organization and contributions

To tackle these challenges, this study proposes an individual-based path
recommendation model to reduce system congestion during public transit
disruptions considering behavior uncertainty. We first formulate an optimal
flow problem as a linear program based on the model of Bertsimas et al.
(2020), which solves the optimal path flows for each OD pair and time interval
that minimize the system travel time. Then, we add the recommendation decision
variables, $x_{p,r}$ (binary variable indicating whether path $r$ is
recommended passenger $p$), and associated constraints to capture the behavior
uncertainty. The behavior uncertainty is modeled with a conditional path
choice probability distribution for each passenger given their received path
recommendation. We introduce two new concepts: $\epsilon$-feasible flows and
$\Gamma$-concentrated flows, to connect the optimal flow problem with the
conditional path choice probabilities. We show that these two concepts can be
seen as a way of approximating the recourse function (expected system travel
time) in a two-stage stochastic optimization. It is proved that these two
concepts help to bound the difference between the approximated recourse
function and the exact one. Additional theoretical analysis shows that
$\epsilon$-feasibility and $\Gamma$-concentration are approximations of
expectation and chance constraints in a typical stochastic optimization
formulation, respectively. The individual path recommendation problem with
behavior uncertainty is a mixed-integer program. We solve it efficiently with
Benders decomposition. The proposed approach is implemented in the Chicago
Transit Authority (CTA) system with a real-world urban rail disruption as the
case study.

The main contributions of this paper are as follows:

  * •

To the best of the authors’ knowledge, this is the first article dealing with
individual path recommendations under public transit service disruptions
considering behavior uncertainty. Previous studies only considered uncertainty
in demand (Mo et al. 2022a) or incident duration (Tan et al. 2020).

  * •

To model behavior uncertainty, this paper proposes a framework with prior path
utility and posterior path choice distribution given recommendations. We use
two new concepts: $\epsilon$-feasibility and $\Gamma$-concentration, to
control the mean and variance of path flows due to behavior uncertainty and
transform these two requirements to linear constraints in the optimization
model using Chebyshev’s inequality.

  * •

The proposed concepts can be seen as a way of approximating the recourse
function (expected system travel time) in a two-stage stochastic optimization.
It is proved that these two concepts help to bound the difference between the
approximated recourse function and the exact one. Additional theoretical
analysis shows that $\epsilon$-feasibility and $\Gamma$-concentration are
approximations of expectation and chance constraints in a typical stochastic
optimization formulation, respectively.

  * •

Benders decomposition (BD) is used to solve the mixed-integer individual path
recommendation problem efficiently. Under BD, the master problem becomes a
small-scale integer program and the sub-problem reduces to a linear program. A
series of theoretical analyses are provided to show the connections between
the proposed concepts and typical stochastic and robust optimization.

The remainder of the paper is organized as follows. The literature review is
discussed in Section 2. In Section 3, we describe the problem conceptually and
analytically. Section 4 develops the solution methods, including the optimal
flow problem formulation, modeling of the behavior uncertainty, theoretical
analysis of the proposed concepts, and Benders decomposition. In Section 5, we
apply the proposed model on the CTA system as a case study. The model results
are analyzed in Section 6. Finally, we conclude the paper and summarize the
main findings in Section 7. All mathematical proofs can be found in
appendices.

##  2 Literature review

###  2.1 Individualized recommendations system

Individualized recommendations design is a popular topic in the field of
computer science and operations research, with many real-world implications
such as Ads ranking (Richardson et al. 2007, Khalid et al. 2013), mobile news
recommendations (Yeung and Yang 2010), travel recommendations (Majid et al.
2013), etc. Most of these recommendation systems focus on individual
preference maximization, which, in return, can increase indicators of interest
such as click-through rate (CTR) and conversion rate. However, in the context
of path recommendations under disruptions, though respecting passengers’
preferences is important, the ultimate goal is to minimize the system travel
time and mitigate the impact of disruptions, which is different from typical
recommendation design literature. Another difference is that, the typical
recommendation systems are usually designed with machine learning algorithms
trained with the real-world user and system interaction data because they have
to learn users’ preferences based on their interaction histories. However, in
this study, the system travel time can be evaluated using a network loading
model. This implies that, instead of using machine learning models, we can use
an optimization formulation to determine the individualized path
recommendations that minimize system travel time.

In summary, different from the typical individualized recommendation system
literature, this study focuses on system-level objectives instead of
individual-level preferences. It leverages an optimization model to design the
recommendation, rather than machine learning models.

###  2.2 Path recommendations during disruptions

Most previous studies on path recommendation under incidents are like
designing a “trip planner”. That is, the main objective is to find available
routes or the shortest path given an OD pair when the network is interrupted
by incidents. For example, Bruglieri et al. (2015) designed a trip planner to
find the fastest path in the public transit network during service disruptions
based on real-time mobility information. Böhmová et al. (2013) developed a
routing algorithm in urban public transportation to find reliable journeys
that are robust for system delays. Roelofsen et al. (2018) provided a
framework for generating and assessing alternative routes in case of
disruptions in urban public transport systems. To the best of the authors’
knowledge, none of the previous studies have considered path recommendations
aiming to minimize the system-wide travel time.

Providing path recommendations during disruptions is similar to the topic of
passenger evacuation under emergencies. The objective of evacuation is usually
to minimize the total evacuation time. For example, Abdelgawad and Abdulhai
(2012) developed an evacuation model with routing and scheduling of subway and
bus transit to alleviate congestion during the evacuation of busy urban areas.
Wang et al. (2019) proposed an optimal bus bridging design method under
operational disruptions on a single metro line. Tan et al. (2020) propose an
evacuation model with urban bus networks as alternatives in the case of common
metro service disruptions by jointly designing the bus lines and frequencies.

However, although these passenger evacuation studies focus on minimizing the
system travel time, there are several differences from this paper. First, in
our paper, the service disruption is not as severe as the emergency situation.
We assume the service will recover after a period of time and passengers are
allowed to wait. They do not necessarily need to cancel trips or follow the
evacuation plan as assumed in previous evacuation studies. Second, in this
article, we do not adjust the operations on the supply side. Instead, we focus
on providing information to the passengers to better utilize the existing
resources/capacities of the system. Third, as mentioned before, this paper
considers passenger heterogeneity and focuses on individual-level path
recommendations, while previous evacuation papers simply model passengers as
flows. Besides, we also assume that passengers may not follow the
recommendation (i.e., behavior uncertainty), which has not been considered in
any evacuation paper before.

###  2.3 Behavior uncertainty

Behavior uncertainty is a well-known challenge in transportation modeling
(Mahmassani 1984). Typically, passenger’s behavior is modeled using various
econometrics approaches (Ben-Akiva et al. 1985, Train 2009, Mo et al. 2021) or
machine learning models (Mirchevska 2013, Wang et al. 2020). These models
output the probability distribution for the passenger’s possible behavior. At
the aggregate level, there are numerous studies using the predicted demand for
different transportation applications taking demand uncertainty into
consideration, such as ride-sharing (Guo et al. 2021), transit route planning
(Yoon and Chow 2020), and supply chain management (Jung et al. 2004).

However, at the individual level, the number of studies is limited. The main
reason is that individual-level decision-making is usually discrete. So it is
challenging to use typical robust optimization to address discrete uncertain
variables (Subramanyam et al. 2021). In terms of stochastic optimization, the
number of possible scenarios increases exponentially with the number of
individuals in the system. Some studies use simulation to incorporate
individual-level behavior uncertainty. For example, Horne et al. (2005) use a
discrete choice model to simulate how different hybrid energy-economy policies
can motivate users’ responses. However, to incorporate behavior uncertainty in
an optimization model (such as the individual path recommendation model in
this paper), new modeling techniques are needed.

Another difference in this study compared to previous literature is that the
behavior uncertainty (i.e., passenger’s response to the recommendation) makes
the decision variables (i.e., passenger flow) become random variables. Typical
robust optimization or stochastic optimization usually assumes the parameters
of constraints are random variables, but not the decision variable.

##  3 Problem description

###  3.1 Conceptual description

Consider a service disruption in an urban rail system. During the disruption,
some stations in the incident line (or the whole line) are blocked. Passengers
in the blocked trains are usually offloaded to the nearest platforms. To
respond to the incident, some operating changes are made, such as dispatching
shuttle buses, rerouting existing services, short-turning in the incident
line, headway adjustment, etc. Assume that all information about the operating
changes is available. These changes define a new PT service network and
available path sets. Our objective is to develop an individual-based path
recommendation model that, when an incident happens, provides a recommended
path to every passenger who uses their phones, websites, or electronic boards
at stations to enter their origin, destination, and departure time. The
recommendation considers the individual’s preferences and behavioral
histories. Hence, passengers with the same origin, destination, and departure
time may get different recommended paths. The overall system aims to minimize
the total travel time for all passengers, including passengers in nearby lines
or bus routes without incidents (note that these passengers may experience
additional crowding due to transfer passengers from the incident line).

Figure 1 shows a simple example of the path recommendation problem. In this
example, Rail Line 1 has an incident and cannot provide service for a period
of time. Both of the two passengers at station A want to go to station C.
Assuming that they request path recommendations. The alternative paths include
using the bus route (blue dashed line), using Rail Line 2 (green dashed line),
or waiting for the system to recover (i.e., still using Rail Line 1). Note
that using either the bus route or Rail Line 2 will take away capacity from
passengers who originally use these two services (i.e., the orange passengers
in the figure). Hence, the model should consider the total travel time of all
four passengers in the system to design recommendation strategies.

Figure 1: Example of the individual path recommendation problem

Moreover, as mentioned in the introduction, behavior uncertainty needs to be
considered. In this example, if we recommend a passenger use a bus route,
he/she may not follow the recommendation and choose Rail Line 2 instead.

###  3.2 Analytical description

Let us divide the analysis period into several time intervals with equal
length $\tau$ (e.g., $\tau=5$ min). Let $t$ be the integer time index. $t=1$
is the start of the incident and $t\leq 0$ indicates the time before the
incident. Let $\mathcal{P}$ be the set of passengers that will receive path
recommendations. We assume $\mathcal{P}$ is known as we can obtain passengers’
requests before running the model. Given the revised operation during the
incident, let $\mathcal{R}_{p}$ be the feasible path set for each passenger
$p\in\mathcal{P}$. Note that $\mathcal{R}_{p}$ includes all feasible services
that are provided by the PT operator. A path $r\in\mathcal{R}_{p}$ may be
waiting for the system to recover (i.e., using the incident line), or transfer
to nearby bus lines, using shuttle services, etc. We do not consider non-PT
modes such as TNC or driving for the following reasons: 1) This study aims to
design a path recommendation system used by PT operators. The major audience
should be all PT users. Considering non-PT modes needs the supply information
of all other travel modes and even consider non-PT users (such as the impact
of traffic congestion on drivers), which is beyond the scope of this study.
Future research may consider a multi-modal path recommendation system. 2)
Passengers using non-PT modes can be simply treated as demand reduction for
the PT system. So their impact on the PT system can still be captured.

Given a passenger $p\in\mathcal{P}$, we aim to determine $x_{p,r}$ for each
$p$, where $x_{p,r}$ indicates whether path $r\in\mathcal{R}_{p}$ is
recommended to passenger $p$ or not. Assume only one path is recommended to
each passenger, we have

$\displaystyle\sum_{r\in\mathcal{R}_{p}}x_{p,r}=1\quad\forall p\in\mathcal{P}$
(1)

Note that we can relax this assumption by designing the recommendation to a
passenger as a “composition” including multiple paths or travel times. This
generalization is discussed in Section 10.1.

$\mathcal{P}$ includes passengers with different origins, destinations, and
departure times. If an incident ends at $t^{\text{end}}$, the recommendation
should consider a time horizon after $t^{\text{end}}$ because there is
remaining congestion in the system. Hence, we provide recommendations until
time $T^{D}>t^{\text{end}}$ (e.g., $T^{D}$ can be one hour after
$t^{\text{end}}$). Therefore, the departure times for passenger
$p\in\mathcal{P}$ range from $[1,T^{D}]$ ($T^{D}$ and $t^{\text{end}}$ are
both time indices).

The recommendation model will be solved at $t=1$ and will generate the
recommendation strategies
$\boldsymbol{x}=(x_{p,r})_{p\in\mathcal{P},r\in\mathcal{R}_{p}}$ for
passengers who depart at time $t\in[1,T^{D}]$. In reality, the model can be
implemented in a rolling horizon manner. Specifically, at each time interval
$t\geq 1$, we first update the demand and supply information in the system,
including new demand estimates, new to-be-recommended passenger set
$\mathcal{P}$, newly available path sets $\mathcal{R}_{p}$, new service routes
and frequencies, new incident duration estimates, new onboard passenger
estimates, etc. Based on this information, we solve the model to obtain
recommendations for passengers with departure time in $[t,T^{D}]$. But we only
implement the recommendation strategies for passengers who depart at the
current time $t$.

Therefore, in the following formulation, we only focus on solving the model at
$t=1$, which is the start of the incident. The whole analysis period includes
warm-up and cool-down periods to better estimate the system states (e.g.,
vehicle loads, passenger travel times, etc.). Therefore, the analysis period
is defined as $[t^{\text{min}},T]$, where $t^{\text{min}}<1$ (time before the
incident) and $T>T^{D}$. For example, $t^{\text{min}}$ and $T$ can be one hour
before and after $t=1$ and $T^{D}$, respectively. And we define all time
intervals in the analysis period as
$\mathcal{T}=\\{t^{\text{min}},t^{\text{min}}+1,...,T\\}$. The overall path
recommendation framework can be summarized in Figure 2.

Figure 2: Problem description and model framework

##  4 Formulation

In this section, we elaborate on the detailed formulation of the individual
path recommendation model. Section 4.1 develops an optimization model to solve
the optimal flow distribution over a public transit network with disruptions.
Section 4.2 describes how passengers’ behavior uncertainties (i.e., non-
compliance to recommendation) are modeled based on a random utility
maximization framework. Section 4.3 provides the overall formulation of the
individual path recommendation model by combining the optimal flow model in
Section 4.1 and the behavior uncertainty component in Section 4.2. Section 4.5
shows how the individual path recommendation model can be solved efficiently
using Benders decomposition.

The notations used in the paper are summarized in Appendix 9.

###  4.1 Optimal flow during disruptions

In this section, we formulate a linear programming (LP) model adapted from
Bertsimas et al. (2020) to solve the optimal flow distribution in a public
transit system with service disruptions. Consider an OD pair $(u,v)$ and
departure time $t$. Let $\mathcal{R}^{u,v}$ be the set of feasible paths for
OD pair $(u,v)$. Define $q_{t}^{u,v,r}$ (resp. $f_{t}^{u,v,r}$) as the number
of passengers in (resp. not in) $\mathcal{P}$ with OD pair $(u,v)$ and
departure time $t$, who use path $r\in\mathcal{R}^{u,v}$. Specifically,
$q_{t}^{u,v,r}$ represents the passenger flows that receive recommendations
while $f_{t}^{u,v,r}$ those do not. Hence, the total path flow in
$r\in\mathcal{R}^{u,v}$ is $q_{t}^{u,v,r}+f_{t}^{u,v,r}$. Let $d_{t}^{u,v}$ be
the total demand of OD pair $(u,v)$ at time $t$, we have

$\displaystyle
q_{t}^{u,v,r}+f_{t}^{u,v,r}=d_{t}^{u,v}\quad\forall(u,v)\in\mathcal{W},t\in\mathcal{T}$
(2)

where $\mathcal{W}$ is the set of all OD pairs. As we focus on path
recommendations for $\mathcal{P}$, in this study, $q_{t}^{u,v,r}$ is the
decision variable while $f_{t}^{u,v,r}$ is a known constant (i.e., the
estimated demand information). For mathematical convenience, we define
$\mathcal{F}$ as the set of all triplets $(u,v,r)$ in the system. And the
objective in this section is to find the optimal flows $q_{t}^{u,v,r}$
($\forall(u,v,r)\in\mathcal{F},t\in\mathcal{T}$) that minimize the total
system travel time.

Consider a path $r$ for OD pair $(u,v)$. A path may include multiple legs,
where each leg is associated with the service in a rail or a bus line. For
example, the path in Figure 3 (indicated by green arrows) has two legs: the
first one in the rail line and the second in the bus line. Every leg has a
boarding and an alighting station. For example, Leg 1 (resp. 2) in this
example has boarding station A (resp. C) and alighting station B (resp. D).
Let $\mathcal{I}^{u,v,r}=\\{1,...,|\mathcal{I}^{u,v,r}|\\}$ be the set of legs
for path $r$. We use a four-element tuple $(u,v,r,i)$ to represent a leg $i$
of path $r$ for OD pair $(u,v)$, where $i\in\mathcal{I}^{u,v,r}$.

Figure 3: Definition of paths and legs

Let $\Delta^{u,v,r,i}_{t}$ (resp. $\delta^{u,v,r,i}_{t}$) be the travel time
between the terminal and the boarding (resp. alighting) station of leg
$(u,v,r,i)$ for a vehicle departing from the terminal at time $t$. Hence, the
vehicle’s arrival time at the boarding (resp. alighting) station of leg
$(u,v,r,i)$ is $t+\Delta^{u,v,r,i}_{t}$ (resp. $t+\delta^{u,v,r,i}_{t}$).
$\delta^{u,v,r,i}_{t}-\Delta^{u,v,r,i}_{t}$ represents the total in-vehicle
time of leg $(u,v,r,i)$ for the vehicle. Define $z^{u,v,r,i}_{t}$ (decision
variable) as the total number of onboard passengers in leg $(u,v,r,i)$ who
board a vehicle that had departed from the terminal at time $t$.

There are three types of constraints for the network flow description: 1)
existing flow constraints, 2) vehicle capacity constraints, and 3) flow
conservation constraints.

Existing flows constraints: Although the path recommendations start at time
$t=1$, there are passengers that already boarded the vehicles. Ignoring these
existing flows may lead to an overestimation of the system’s available
capacity. To capture the existing onboard flows at $t=1$, we define the set of
onboard flow indices as

$\displaystyle\Omega_{1}=\\{(u,v,r,i,t):t+\Delta^{u,v,r,i}_{t}\leq 1\leq
t+\delta^{u,v,r,i}_{t}\\}$ (3)

And the existing flow constraints can be expressed as

$\displaystyle
z_{t}^{u,v,r,i}=\hat{z}_{t}^{u,v,r,i}\quad\forall(u,v,r,i,t)\in\Omega_{1}$ (4)

where $\hat{z}_{t}^{u,v,r,i}$ are constants that capture the existing onboard
flows when the incident happens. These flows can be directly obtained from a
simulation model or real-time passenger counting data.

Capacity constraints: Transit vehicles have limited capacity. Consider a
vehicle departing at time $t$ on line $l$ (referred to as vehicle $(l,t)$). We
denote its total number of onboard passengers at time $t^{\prime}$ as
$O_{l,t,t^{\prime}}$. Specifically, $O_{l,t,t^{\prime}}$ can be expressed as

$\displaystyle
O_{l,t,t^{\prime}}(\boldsymbol{z})=\sum_{\\{(u,v,r,i,t)\in\texttt{Onboard}(l,t^{\prime})\\}}z_{t}^{u,v,r,i}\quad\forall
l\in\mathcal{L},\forall t\in\mathcal{T},t^{\prime}=t,t+1,...,T_{l,t}$ (5)

where $T_{l,t}$ is the time index that vehicle $(l,t)$ arrives at the last
station of line $l$,
$\boldsymbol{z}=(z_{t}^{u,v,r,i})_{t\in\mathcal{T},(u,v,r)\in\mathcal{F},i\in\mathcal{I}^{u,v,r}}$.
$\texttt{Onboard}(l,t^{\prime})$ is the set of onboard flow indices for
vehicle $(l,t)$, defined as

$\displaystyle\texttt{Onboard}(l,t^{\prime})=\\{(u,v,r,i,t):\text{Leg
}(u,v,r,i)\text{ on line $l$, and }t+\Delta_{t}^{u,v,r,i}\leq t^{\prime}\leq
t+\delta_{t}^{u,v,r,i}\\}$ (6)

Then the capacity constraint is:

$\displaystyle O_{l,t,t^{\prime}}(\boldsymbol{z})\leq K_{l,t}\quad\forall
l\in\mathcal{L},t\in\mathcal{T},t^{\prime}=t,t+1,...,T_{l,t}$ (7)

where $K_{l,t}$ is the capacity of the vehicle $(l,t)$. $\mathcal{L}$ is the
set of all lines.

Flow conservation constraint: There are two different flow conservation
constraints: 1) flow conservation at origin stations and 2) at transfer
stations. To ensure the origin flow conservation, the cumulative number of
arrival passengers should be larger than the cumulative number of boarding
passengers at an origin at any time. This indicates that not all arrival
passengers can board due to potentially being left behind because of capacity
constraints.

The number of arriving passengers (i.e., demand) for $(u,v,r)$ at time $t$ is
$q_{t}^{u,v,r}+f_{t}^{u,v,r}$. And the number of boarding passengers at the
origin station (i.e., $u$) at time $t$ is $z_{t^{\prime}}^{u,v,r,1}$ (i.e.,
the first leg) with $t^{\prime}+\Delta^{u,v,r,1}_{t^{\prime}}=t$. $t^{\prime}$
is the vehicle departure time from the terminal and
$t^{\prime}+\Delta^{u,v,r,1}_{t^{\prime}}$ is the time when the vehicle
arrives at the boarding station. Therefore, the origin flow conservation
constraint can be written as:

$\displaystyle\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\Delta^{u,v,r,1}_{t^{\prime}}\leq
t\\}}z_{t^{\prime}}^{u,v,r,1}\leq\sum_{t^{\prime}=t^{\text{min}}}^{t}\left(f_{t^{\prime}}^{u,v,r}+q_{t^{\prime}}^{u,v,r}\right)\quad\forall(u,v,r)\in\mathcal{F},t\in\mathcal{T}$
(8)

Now consider the flow conservation at a transfer station. All arrival
passengers at a transfer station of a path are the onboard passengers from the
last leg. Therefore, we use a similar way to define the transfer flow
conservation: the cumulative number of onboard passengers from the last leg
should be larger than the cumulative number of boarding passengers at the
transfer station. And the number of boarding passengers at the transfer
station is simply $z_{t^{\prime}}^{u,v,r,i}$ with $i\geq 2$. Hence, flow
conservation constraints at a transfer station are:

$\displaystyle\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\Delta_{t^{\prime}}^{u,v,r,i}\leq
t\\}}z_{t^{\prime}}^{u,v,r,i}\leq\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\delta_{t^{\prime}}^{u,v,r,i-1}\leq
t\\}}z_{t^{\prime}}^{u,v,r,i-1}\quad\forall(u,v,r)\in\mathcal{F},i\in\mathcal{I}^{(u,v,r)}\setminus\\{1\\},t\in\mathcal{T}$
(9)

Note that $z_{t^{\prime}}^{u,v,r,i}$ is defined as the onboard passengers for
vehicles departing at time $t^{\prime}$. Therefore,
$t^{\prime}+\delta_{t^{\prime}}^{u,v,r,i-1}$ is the alighting time for
passengers at leg $i-1$ (which is also the transfer demand arrival time at leg
$i$ as we assume transfer walk time is within a time interval $\tau$ and is
negligible). $t^{\prime}+\Delta_{t^{\prime}}^{u,v,r,i}$ is the boarding time
for passengers at leg $i$.

The objective is to minimize the total travel time for all passengers in the
system. Total travel time can be decomposed into waiting time and in-vehicle
time.

In-vehicle time: Total in-vehicle time is simply the onboard flow multiplied
by the travel time on each leg:

$\displaystyle
IVT(\boldsymbol{z})=\sum_{(u,v,r)\in\mathcal{F}}\sum_{i\in\mathcal{I}^{u,v,r}}\sum_{t\in\mathcal{T}}z^{u,v,r,i}_{t}\cdot
T^{\text{IVT}}_{u,v,r,i,t}$ (10)

where $T^{\text{IVT}}_{u,v,r,i,t}$ is the in-vehicle time of leg $(u,v,r,i)$
of the vehicle departing at time $t$.

Waiting time: There are two causes of waiting time: 1) waiting time because of
vehicle headways, and 2) waiting time resulting from being left behind. During
a specific time interval $t$, all left behind passengers would have a waiting
time of $\tau$. All boarding passengers, assuming uniform arrival, have an
average waiting time that is half of the time interval (i.e.,
$\frac{\tau}{2}$). Therefore, the total waiting time for passengers at station
$s$ and time $t$ can be formulated as

$\displaystyle
WT_{s,t}=\tau(AD_{s,t}+XD_{s,t}-BD_{s,t})+\frac{\tau}{2}(BD_{s,t+1}-BD_{s,t})$
(11)

where $AD_{s,t}$ represents the cumulative arriving demand at station $s$ up
to time $t$, $XD_{s,t}$ represents the cumulative transferring demand at
station $s$ up to time $t$, and $BD_{s,t}$ represents the cumulative boarded
demand at station $s$ up to time $t$. Hence, $(BD_{s,t+1}-BD_{s,t})$
represents the total number of boarding passengers at time $t$ and station
$s$, and $(AD_{s,t}+XD_{s,t}-BD_{s,t})$ represents the total number of left
behind passengers at station $s$ and time $t$. Finally, the total system
waiting time is

$\displaystyle
WT(\boldsymbol{q},\boldsymbol{z})=\sum_{s\in\mathcal{S}}\sum_{t=1}^{T}WT_{s,t}$
(12)

where
$\boldsymbol{q}=(q_{t}^{u,v,r})_{t\in\mathcal{T},(u,v,r)\in\mathcal{F}}$.

The cumulative arriving demand $AD_{s,t}$ is simply all arriving passengers
with origin $s$ up to time $t$:

$\displaystyle
AD_{s,t}=\sum_{\\{(u,v,r):u=s\\}}\sum_{t^{\prime}=t^{\text{min}}}^{t}\left(f_{t^{\prime}}^{u,v,r}+q_{t^{\prime}}^{u,v,r}\right)\quad\forall
s\in\mathcal{S},t\in\mathcal{T}$ (13)

where $\mathcal{S}$ is the set of all stations.

The cumulative transferring demand is all passengers alighting at station $s$
from their previous leg $i-1$ for their next leg $i$:

$\displaystyle
XD_{s,t}=\sum_{\\{(u,v,r,i)\in\texttt{Xth}(s)\\}}\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\delta^{u,v,r,i-1}_{t^{\prime}}\leq
t\\}}z_{t^{\prime}}^{u,v,r,i-1}\quad\forall t=t^{\text{min}},...,T$ (14)

where $\texttt{Xth}(s)$ is the set of legs that transfer at station $s$.

The cumulative boarded demand is all passengers that successfully board a
vehicle at station $s$ at time $t$. Define $\texttt{Bdat}(s)$ as the set of
all legs with boarding station $s$, we have

$\displaystyle
BD_{s,t}=\sum_{\\{(u,v,r,i)\in\texttt{Bdat}(s)\\}}\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\Delta^{u,v,r,i}_{t^{\prime}}\leq
t\\}}z_{t^{\prime}}^{u,v,r,i}\quad\forall t=t^{\text{min}},...,T$ (15)

Taking everything into consideration, the total travel time in the system is
$WT(\boldsymbol{x},\boldsymbol{z})+IVT(\boldsymbol{z})$. The optimal flow
problem is:

$\displaystyle(OF)\quad\min_{\boldsymbol{q},\boldsymbol{z}}\quad$
$\displaystyle WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})$ (16a)
s.t. $\displaystyle O_{l,t,t^{\prime}}(\boldsymbol{z})\leq K_{l,t}\quad\forall
l\in\mathcal{L},t\in\mathcal{T},t^{\prime}=t,t+1,...,T_{l,t}$ (16b)
$\displaystyle\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\Delta^{u,v,r,1}_{t^{\prime}}\leq
t\\}}z_{t^{\prime}}^{u,v,r,1}\leq\sum_{t^{\prime}=t^{\text{min}}}^{t}\left(f_{t^{\prime}}^{u,v,r}+q_{t^{\prime}}^{u,v,r}\right)\quad\forall(u,v,r)\in\mathcal{F},t\in\mathcal{T}$
(16c) $\displaystyle\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\Delta_{t^{\prime}}^{u,v,r,i}\leq
t\\}}z_{t^{\prime}}^{u,v,r,i}\leq\sum_{\\{t^{\prime}:t^{\text{min}}\leq
t^{\prime}+\delta_{t^{\prime}}^{u,v,r,i-1}\leq
t\\}}z_{t^{\prime}}^{u,v,r,i-1}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\forall(u,v,r)\in\mathcal{F},i\in\mathcal{I}^{(u,v,r)}\setminus\\{1\\},t\in\mathcal{T}$
(16d)
$\displaystyle\sum_{r\in\mathcal{R}^{u,v}}q^{u,v,r}_{t}+f^{u,v,r}_{t}=d^{u,v}_{t}\quad\forall(u,v)\in\mathcal{W},t\in\mathcal{T}$
(16e) $\displaystyle
z_{t}^{u,v,r,i}=\hat{z}_{t}^{u,v,r,i}\quad\forall(u,v,r,i,t)\in\Omega_{1}$
(16f) $\displaystyle z^{u,v,r,i}_{t}\geq 0\quad\forall
t\in\mathcal{T},(u,v,r)\in\mathcal{F},i\in\mathcal{I}^{u,v,r}$ (16g)
$\displaystyle q^{u,v,r}_{t}\geq 0\quad\forall
t\in\mathcal{T},(u,v,r)\in\mathcal{F},$ (16h)

As the objective function is minimizing the system travel time, this
formulation will automatically load passengers to a train as long as there is
available capacity (Bertsimas et al. 2020).

Path travel time calculation: It is worth noting that Eq. 16 does not
explicitly output the travel time of different paths. The travel time of a
path $(u,v,r)$ for trips departing at time $t$ (denoted as $TT^{u,v,r}_{t}$)
has to be obtained from the network flow patterns after solving Eq. 16.
Specifically, consider the group of passengers using path $(u,v,r)$ and
departing at time $t$. Their arrival time at the destination (denoted as
$AT^{u,v,r}_{t}$) can be calculated as

$\displaystyle
AT^{u,v,r}_{t}=\min\Bigg{\\{}\tilde{t}\in\mathcal{T}^{u,v,r}_{t}:\sum_{t^{\prime}=t^{\text{min}}}^{t}\left(f_{t^{\prime}}^{u,v,r}+q_{t^{\prime}}^{u,v,r}\right)\leq\sum_{t^{\text{min}}\leq
t^{\prime}+\delta^{u,v,r,|\mathcal{I}^{u,v,r}|}_{t^{\prime}}\leq{\tilde{t}}}z^{u,v,r,|\mathcal{I}^{u,v,r}|}_{t^{\prime}}\Bigg{\\}}\quad$
$\displaystyle\qquad\qquad\qquad\forall t\in\mathcal{T},(u,v,r)\in\mathcal{F}$
(17)

where $\mathcal{T}^{u,v,r}_{t}$ is the set of possible arrival time indices,
defined as $\mathcal{T}^{u,v,r}_{t}=\\{t^{\prime}:t\leq t^{\prime}\leq T\\}$.
Eq. 17 represents the travel time calculation with cumulative demand curves at
origins and destinations.
$\sum_{t^{\prime}=t^{\text{min}}}^{t}\left(f_{t^{\prime}}^{u,v,r}+q_{t^{\prime}}^{u,v,r}\right)$
is the cumulative demand up to time $t$ at the origin.
$\sum_{t^{\text{min}}\leq
t^{\prime}+\delta^{u,v,r,|\mathcal{I}^{u,v,r}|}_{t^{\prime}}\leq{\tilde{t}}}z^{u,v,r,|\mathcal{I}^{u,v,r}|}_{t^{\prime}}$
is the cumulative passengers arriving at the destination up to time
$t^{\prime}$. When the cumulative arrivals at the destination are greater or
equal to the cumulative demand at the origin (up to time $t$), all passengers
finish the trip. So taking the minimum over $t^{\prime}$ gives the arrival
time for passengers departing at $t$. The path travel time is then simply:

$\displaystyle TT^{u,v,r}_{t}=AT^{u,v,r}_{t}-t\quad\forall
t\in\mathcal{T},(u,v,r)\in\mathcal{F}$ (18)

Figure 4 illustrates the travel time calculation.

Figure 4: Travel time calculation

Incident specification: Eq. 16 is a general formulation of the optimal flow
problem. Now we will introduce how the incident-specific information is
incorporated into this problem. We assume the incident causes a service
disruption in a specific line (if only several stations are interrupted, we
can separate the line into multiple lines so that the assumption always
holds). The service disruption in a line can be seen as stops of vehicles for
a period of time. The vehicle stopping can be captured by the parameters
$\Delta^{u,v,r,i}_{t}$, $\delta^{u,v,r,i}_{t}$, and $K_{l,t}$. Specifically, a
long stop due to an incident can be seen as an increase in travel time from
the terminal to downstream stations (i.e., increase in $\Delta^{u,v,r,i}_{t}$
and $\delta^{u,v,r,i}_{t}$). Moreover, since there is no vehicle dispatching
during the incident, we set $K_{l,t}=0$ for the corresponding time and line.
In this way, we can model the incident without changing the formulation.

###  4.2 Behavior uncertainty

Consider a passenger $p$ with a path set $\mathcal{R}_{p}$. Their inherent
preference (utility) of using path $r$ is denoted as $V_{p}^{r}$. If path
$r^{\prime}$ was recommended, the impact of the recommendation on the utility
of path $r$ is denoted as $I_{p,r^{\prime}}^{r}$. Hence, his/her overall
utility of using path $r$ can be represented as

$\displaystyle
U_{p}^{r}=V_{p}^{r}+\sum_{r^{\prime}\in\mathcal{R}_{p}}x_{p,r^{\prime}}\cdot
I_{p,r^{\prime}}^{r}+\xi_{p}^{r}\quad\forall
r\in\mathcal{R}_{p},\;p\in\mathcal{P}.$ (19)

where $\xi_{p}^{r}$ is the random error. $x_{p,r^{\prime}}=1$ if passenger $p$
is recommended path $r^{\prime}$, otherwise $x_{p,r^{\prime}}=0$. Let
$\pi_{p,r^{\prime}}^{r}$ be the conditional probability that passenger $p$
chooses path $r$ given that the recommended path is $r^{\prime}$. Assuming a
utility-maximizing behavior, we have

$\displaystyle\pi_{p,r^{\prime}}^{r}=\mathbb{P}(V_{p}^{r}+I_{p,r^{\prime}}^{r}+\xi_{p}^{r}\geq
V_{p}^{r^{\prime\prime}}+I_{p,r^{\prime}}^{r^{\prime\prime}}+\xi_{p}^{r^{\prime\prime}},\;\forall
r^{\prime\prime}\in\mathcal{R}_{p})$ (20)

Different assumptions for the distribution of $\xi_{p}^{r}$ can lead to
different expressions. For example, if $\xi_{p}^{r}$ are i.i.d. Gumbel
distributed, the choice probability reduces to multinomial logit model (Train
2009) and we have

$\displaystyle\pi_{p,r^{\prime}}^{r}=\frac{\exp(V_{p}^{r}+I_{p,r^{\prime}}^{r})}{\sum_{r^{\prime\prime}\in\mathcal{R}_{p}}\exp(V_{p}^{r^{\prime\prime}}+I_{p,r^{\prime}}^{r^{\prime\prime}})}$
(21)

The value of $V_{p}^{r}$ and $I_{p,r^{\prime}}^{r}$ can be calibrated using
data from individual-level surveys or smart card data, which deserves separate
research. When developing the individual path recommendation model, we assume
$\pi_{p,r^{\prime}}^{r}$ is known. Figure 5 shows an example of the
conditional probability matrix. The specific values assume that paths with
recommendations are more likely to be chosen.

Figure 5: Example of conditional path choice probability

The conditional probability $\pi_{p,r^{\prime}}^{r}$ captures the individual’s
inherent preference for different paths as well as the response to the
recommendation system. It varies across individuals and reflects their
behavioral uncertainties. This study focuses on design path recommendation
systems based on the value of $\pi_{p,r^{\prime}}^{r}$.

###  4.3 Individual path recommendation

Let $\mathbbm{1}_{p,r^{\prime}}^{r}$ be the indicator random variable
representing whether passenger $p$ actually chooses path $r$ or not given that
he/she is recommended path $r^{\prime}$. By definition,
$\mathbbm{1}_{p,r^{\prime}}^{r}$ is a Bernoulli random variable with
$\mathbb{E}[\mathbbm{1}_{p,r^{\prime}}^{r}]=\pi_{p,r^{\prime}}^{r}$ and
$\text{Var}[\mathbbm{1}_{p,r^{\prime}}^{r}]=\pi_{p,r^{\prime}}^{r}\cdot(1-\pi_{p,r^{\prime}}^{r})$

Therefore, the actual flow for path $(u,v,r)$ at time $t$ is

$\displaystyle
Q_{t}^{u,v,r}=\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\mathbbm{1}_{p,r^{\prime}}^{r}$
(22)

$Q_{t}^{u,v,r}$ is also a random variable.
$\mathcal{P}^{u,v}_{t}\subseteq\mathcal{P}$ is the set of passengers with OD
pair $(u,v)$ arriving at the system at time interval $t$ that receive path
recommendations. $\mathcal{R}^{u,v}$ is the set of paths of OD pair $(u,v)$.
The mean and variance of the actual flow is

$\displaystyle{\mu}_{t}^{u,v,r}(\boldsymbol{x}):=\mathbb{E}\left[Q_{t}^{u,v,r}\right]=\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}$
(23)
$\displaystyle({\sigma}_{t}^{u,v,r}(\boldsymbol{x}))^{2}:=\text{Var}\left[Q_{t}^{u,v,r}\right]=\sum_{p\in\mathcal{P}^{u,v}_{t^{\prime}}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}\cdot(1-\pi_{p,r^{\prime}}^{r})$
(24)

Note that Eqs. 24 is based on the fact that
$x_{p,r^{\prime}}^{2}=x_{p,r^{\prime}}$ and
$\text{Cov}[\mathbbm{1}_{p,r^{\prime}}^{r},\mathbbm{1}_{p,r^{\prime\prime}}^{r}]=0$
if $r^{\prime}\neq r^{\prime\prime}$.

In an optimization model, we cannot use a random variable (e.g., actual flow)
as the decision variable. Therefore, let us treat $\boldsymbol{q}$ in Eq. 16
as a realization of the random variable flow
$\boldsymbol{Q}=(Q_{t}^{u,v,r})_{t\in\mathcal{T},(u,v,r)\in\mathcal{F}}$. To
make $\boldsymbol{q}$ a reasonable realization, some constraints need to be
considered between the value of $\boldsymbol{q}$ and the distribution of
$\boldsymbol{Q}$. We define two new concepts: “$\epsilon$-feasibility” and
“$\Gamma$-concentration”.

######  Definition 4.1 ($\epsilon$-feasible flows)

A flow $q_{t}^{u,v,r}$ is $\epsilon$-feasible if and only if

$\displaystyle\left|q_{t}^{u,v,r}-{\mu}_{t}^{u,v,r}(\boldsymbol{x})\right|\leq\epsilon_{t}^{u,v,r},\quad\forall(u,v,r)\in\mathcal{F},t=t^{\text{min}},...,T$
(25)

where $\epsilon_{t}^{u,v,r}$ is a small positive constant. This means that
$\boldsymbol{q}$ is close to the expectation of the actual flow under
recommendation strategy $\boldsymbol{x}$.

######  Definition 4.2 ($\Gamma$-concentrated flows)

A flow $q_{t}^{u,v,r}$ is $\Gamma$-concentrated if and only if it is
$\epsilon$-feasible and for any constant $a>\epsilon_{t}^{u,v,r}$, we have

$\displaystyle\mathbb{P}\left[\left|Q_{t}^{u,v,r}-{q}_{t}^{u,v,r}\right|\geq
a\right]\leq\left(\frac{\Gamma_{t}^{u,v,r}}{a-\epsilon_{t}^{u,v,r}}\right)^{2}\quad\forall(u,v,r)\in\mathcal{F},t=t^{\text{min}},...,T$
(26)

where $\Gamma_{t}^{u,v,r}$ is a small positive constant. This means that the
probability that $Q_{t}^{u,v,r}$ and ${q}_{t}^{u,v,r}$ are very different
(i.e., with difference greater than $a$) is bounded from above, suggesting
that $Q_{t}^{u,v,r}$ is concentrated around $q_{t}^{u,v,r}$.

###### Remark 4.3

The logic of using $\boldsymbol{q}$ as the decision variable and defining the
above two concepts is as follows. The objective of this study is to find the
best recommendation strategy $\boldsymbol{x}$ that minimizes the system travel
time. The system travel time is a function of network flows. Given a
recommendation strategy $\boldsymbol{x}$, the actual flow $\boldsymbol{Q}$ is
a random variable, which cannot be directly used in the optimization model (as
decision variables) to evaluate the system travel time. Hence, we assume that
$\boldsymbol{q}$ in Eq. 16 is a realization of the actual flow (deterministic
variable). We also add two constraints to $\boldsymbol{q}$ so that
$\boldsymbol{q}$ is close to the mean of the actual flow, and the distribution
of the actual flow is concentrated around $\boldsymbol{q}$. Then, using
$\boldsymbol{q}$ to evaluate the system travel time is similar to that of
using the actual flows (or taking the expectation). In section 4.4.2, we show
that these two concepts help to bound the difference between the expected
system travel time and the system travel time evaluated using
$\boldsymbol{q}$.

Note that one may argue that we can directly use
${\mu}_{t}^{u,v,r}(\boldsymbol{x})$ as decision variables to represent network
flows and eliminate $\boldsymbol{q}$. This idea is essentially equivalent to
setting $\epsilon_{t}^{u,v,r}=0$ and does not consider the concentration
property (i.e., $\Gamma_{t}^{u,v,r}=+\infty$), which is a special case of our
framework. Our framework has more advantages in controlling the variance.
Specifically, ignoring the $\Gamma$-concentration may make the model
recommendation strategies meaningless. Consider an extreme scenario that there
is a recommendation strategy $\boldsymbol{x}$, under which the actual flow is
uniformly distributed in $[0,1]$. Further, assume that the system travel time
is simply a linear function of the actual flow, say the factor is 1 (i.e., the
system travel time is also uniformly distributed in $[0,1]$). Suppose that the
recommendation strategy $\boldsymbol{x}$ minimizes the expected system travel
time (now the system travel time is
$1\times{\mu}_{t}^{u,v,r}(\boldsymbol{x})=0.5$). However, as the actual flow
can be any value between 0 and 1 with equal probability, we know that the
actual system travel time can also be any value between 0 and 1. Hence, the
recommendation strategy $\boldsymbol{x}$, though minimizing the expected
system travel time, is meaningless because there are too many variations for
the actual system travel time under this recommendation.
$\Gamma$-concentration is an important property to ensure that the
distribution of actual flows is not too dispersed111In reality, as the actual
flow is the summation of many Bernoulli random variables, the coefficient of
variation will shrink with the increase in passenger size. So in the case of a
large number of passengers, the $\Gamma$-concentration should be naturally
satisfied so that the recommendation strategy $\boldsymbol{x}$ is solved based
on a reliable estimate of the system travel time. In Section 4.4.2, we show
that controlling the variance is important if we wish to minimize the expected
system travel time. In Sections 4.4.1 and 4.4.3, we also elaborate on how
these two concepts relate to typical stochastic optimization methods.

We will therefore incorporate $\epsilon$-feasibility and
$\Gamma$-concentration as constraints into the optimization formulation (Eq.
16). It turns out that both of them can be modeled as linear constraints.
$\epsilon$-feasibility (Eq. 25) can be easily transformed into a linear
constraint by eliminating the absolute value. To incorporate
$\Gamma$-concentration (Eq. 26), the following Proposition is used:

###### Proposition 4.4

The $\Gamma$-concentration inequality (Eq. 26) holds if the variance of
$Q_{t}^{u,v,r}$ is bounded from above by $(\Gamma_{t}^{u,v,r})^{2}$.
Mathematically:

$\displaystyle\sum_{p\in\mathcal{P}^{u,v}_{t^{\prime}}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}\cdot(1-\pi_{p,r^{\prime}}^{r})\leq(\Gamma_{t}^{u,v,r})^{2}$
(27)

For modeling convenience, we set
$\epsilon_{t}^{u,v,r}=\epsilon\cdot{\mu}_{t}^{u,v,r}(\boldsymbol{x})$ and
$\Gamma_{t}^{u,v,r}=\Gamma\cdot{d}_{t}^{u,v}$, where $\epsilon$ and $\Gamma$
are hyper-parameters determining how close and concentrated the actual flow
should be. Then the final constraint becomes:

$\displaystyle(1-\epsilon)\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}\leq
q_{t}^{u,v,r}\leq(1+\epsilon)\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}$
(28)

and

$\displaystyle\sum_{p\in\mathcal{P}^{u,v}_{t^{\prime}}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}\cdot(1-\pi_{p,r^{\prime}}^{r})\leq(\Gamma\cdot{d}_{t}^{u,v})^{2}$
(29)

Both constraints are linear and can be added into Eq. 16.

Besides the total system travel time, many recommendation systems also aim to
respect passengers’ preferences. That is, if possible, a path with high
inherent utility $V_{p}^{r}$ should be recommended. Hence the following term
is added to the objective function.

$\displaystyle\max\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}x_{p,r}\cdot
V_{p}^{r}\Longleftrightarrow\min\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}-x_{p,r}\cdot
V_{p}^{r}$ (30)

The final individual path recommendation (IPR) model can be formulated as:

$\displaystyle(IPR)\quad\min_{\boldsymbol{x},\boldsymbol{q},\boldsymbol{z}}\quad$
$\displaystyle
WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})+\Psi\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}-x_{p,r}\cdot
V_{p,r}$ (31a) s.t. $\displaystyle\text{Constraints
}(\ref{const_OF1})-(\ref{const_OF2})$ (31b)
$\displaystyle(1-\epsilon)\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}\leq
q_{t}^{u,v,r}\leq(1+\epsilon)\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\forall
t\in\mathcal{T},(u,v,r)\in\mathcal{F}$ (31c)
$\displaystyle\sum_{p\in\mathcal{P}^{u,v}_{t^{\prime}}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}\cdot(1-\pi_{p,r^{\prime}}^{r})\leq(\Gamma\cdot{d}_{t}^{u,v})^{2}\quad\forall
t\in\mathcal{T},(u,v,r)\in\mathcal{F}$ (31d) $\displaystyle\sum_{r\in
R_{p}}x_{p,r}=1\quad\forall p\in\mathcal{P}$ (31e) $\displaystyle
x_{p,r}\in\\{0,1\\}\quad\forall p\in\mathcal{P},r\in\mathcal{R}_{p}$ (31f)

where $\Psi$ is a hyper-parameter to adjust the scale and balance the trade-
off between system efficiency and passenger preferences.

###  4.4 Discussions on $\epsilon$-feasibility and $\Gamma$-concentration

####  4.4.1 Connections to two-stage stochastic optimization

As the problem needs to solve both recommendation strategy $\boldsymbol{x}$
and path flow $\boldsymbol{q}$, a possible alternative formulation is a two-
stage stochastic optimization, where the first stage is to determine the
recommendation and the second stage is to determine the path flow:

$\displaystyle\quad\min_{\boldsymbol{x}}\quad$
$\displaystyle\Psi\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}-x_{p,r}\cdot
V_{p,r}+\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}})]$
(32a) s.t. $\displaystyle\text{Constraints
}(\ref{const_sum_x1})-(\ref{const_IRR2})$ (32b)

where $\boldsymbol{Q}|_{\boldsymbol{x}}$ is the path flow (random variable)
given recommendation strategy $\boldsymbol{x}$ and
$STT(\boldsymbol{Q}|_{\boldsymbol{x}})$ is the actual system travel time (STT)
defined as

$\displaystyle
STT(\boldsymbol{Q}|_{\boldsymbol{x}})=\min_{\boldsymbol{q},\boldsymbol{z}\in\mathcal{X}^{\text{OF}}(\boldsymbol{Q}|_{\boldsymbol{x}})}WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})$
(33)

where
$\mathcal{X}^{\text{OF}}(\boldsymbol{Q}|_{\boldsymbol{x}})=\\{(\boldsymbol{q,z}):\text{Constraints
}(\ref{const_OF1})-(\ref{const_OF2}),\boldsymbol{q}=\boldsymbol{\boldsymbol{Q}|_{\boldsymbol{x}}}\\}$
is the feasible region of the optimal flow problem (Eq. 16) with an additional
constraint of fixing $\boldsymbol{q}$ to the value of
$\boldsymbol{Q}|_{\boldsymbol{x}}$.

Our formulation differs from the typical two-stage stochastic optimization
problem because the distribution of the uncertain parameter
$\boldsymbol{Q}|_{\boldsymbol{x}}$ also depends on the first-stage decision
$\boldsymbol{x}$. While the typical two-stage stochastic optimization assumes
the uncertain parameters in the second stage have some predetermined exogenous
distribution (Ahmed 2010). The difference comes from that
$\boldsymbol{q}=\boldsymbol{\boldsymbol{Q}|_{\boldsymbol{x}}}$ essentially
makes the decision variable become a random variable, as we mentioned in
Section 4.3. Section 4.4.3 discusses how to deal with random decision
variables in stochastic optimization and how that connects with
$\epsilon$-feasibility and $\Gamma$-concentration.

A typical way to solve the two-stage stochastic optimization is to construct
an approximation $\hat{STT}(\boldsymbol{x})$ for
$\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}})]$.
Then we solve

$\displaystyle\quad\min_{\boldsymbol{x}}\quad$
$\displaystyle\Psi\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}-x_{p,r}\cdot
V_{p,r}+\hat{STT}(\boldsymbol{x})$ (34a) s.t. $\displaystyle\text{Constraints
}(\ref{const_sum_x1})-(\ref{const_IRR2})$ (34b)

for the first stage instead222The approximation is then updated based on the
second-stage solutions. From this perspective, we can treat
$\epsilon$-feasibility and $\Gamma$-concentration as a way of constructing
$\hat{STT}(\boldsymbol{x})$, that is,

$\displaystyle\quad\hat{STT}(\boldsymbol{x})=\min_{\boldsymbol{q,z}}\quad$
$\displaystyle WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})$ (35a)
s.t. $\displaystyle\text{Constraints }(\ref{const_OF1})-(\ref{const_OF2})$
(35b) $\displaystyle\text{Constraints }(\ref{const_IRR1})-(\ref{const_Var})$
(35c)

Therefore, combining Eqs. 34 and 35 as a one-stage optimization problem yields
our IPR formulation (Eq. 31).

However, a natural question would be how good the approximation is. Section
4.4.2 discusses the system travel time difference bounds that partially answer
this question.

####  4.4.2 Difference between the expected system travel time and the IPR
objective function

Let $(\boldsymbol{{q}^{*},{z}^{*},{x}^{*}})$ be the optimal solution of Eq.
31. Given the optimal recommendation strategy $\boldsymbol{x}^{*}$, we denote
the random variable of path flow as $\boldsymbol{Q}|_{\boldsymbol{x}^{*}}$.
Let $\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\cdot)$ be the density
function of $\boldsymbol{Q}|_{\boldsymbol{x}^{*}}$. The expectation of
$STT(\boldsymbol{Q}|_{\boldsymbol{x}^{*}})$ is

$\displaystyle\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}^{*}})]=\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})}\left[\min_{\boldsymbol{q},\boldsymbol{z}\in\mathcal{X}^{\text{OF}}(\hat{\boldsymbol{q}})}WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})\right]\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
(36)

where

$\displaystyle\mathcal{Q}({\boldsymbol{x}^{*}})=\\{\hat{\boldsymbol{q}}\geq
0:\hat{q}_{t}^{u,v,r}=\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}^{*}\cdot\hat{\mathbbm{1}}_{p,r^{\prime}}^{r},\;\forall\;\hat{\mathbbm{1}}_{p,r^{\prime}}^{r}\in\\{0,1\\},\sum_{r\in\mathcal{R}_{p}}\hat{\mathbbm{1}}_{p,r^{\prime}}^{r}=1,\;\forall(u,v,r)\in\mathcal{F},t\in\mathcal{T}\\}$
(37)

is the set of all possible values of network flows given recommendation
$\boldsymbol{x}^{*}$.

$\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}^{*}})]$
is usually the indicator of the system performance. However, our model is
optimized with Eq. 31, where the minimization is conducted over the model-
evaluated system travel time with optimal realized flow $\boldsymbol{q}^{*}$:

$\displaystyle
STT({\boldsymbol{q}^{*}})=\min_{\boldsymbol{q},\boldsymbol{z}\in\mathcal{X}^{\text{OF}}(\boldsymbol{q}^{*})}WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})$
(38)

It is worth analyzing the relationship between the model-evaluated system
travel time ($STT({\boldsymbol{q}^{*}})$) and the expected system travel time
($\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}^{*}})]$).
This analysis tells us how well our proposed approach can approximate the real
system performance indicator.

We first introduce a Lemma based on Berge’s Maximum Theorem (Berge 1957).

###### Lemma 4.5

$STT(\hat{{\boldsymbol{q}}})=\min_{\boldsymbol{q},\boldsymbol{z}\in\mathcal{X}^{\text{OF}}(\hat{\boldsymbol{q}})}WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})$
is continuous in terms of $\hat{{\boldsymbol{q}}}$ if the set of optimal flows
is bounded (i.e., there are a limited number of flow patterns that permits the
optimal system travel time).

Lemma 4.5 implies that a small change in the path flows only results in small
changes in system travel time. Since the system travel time is usually bounded
from above given a finite-scale transit network, unit flow changes should not
yield infinite changes in the system travel time. Hence,
$STT(\hat{{\boldsymbol{q}}})$ should have bounded gradient. Combining the
continuity property in Lemma 4.5, we conclude that
$STT(\hat{{\boldsymbol{q}}})$ is Lipschitz continuous. That is, there exists a
constant $L$ such that, for any network flows $\boldsymbol{q}_{1}$ and
$\boldsymbol{q}_{2}$, we have

$\displaystyle|STT({\boldsymbol{q}}_{1})-STT({\boldsymbol{q}}_{2})|\leq
L\cdot\left\lVert\boldsymbol{q_{1}}-\boldsymbol{q_{2}}\right\rVert_{1}$ (39)

###### Proposition 4.6

Let $(\boldsymbol{{q}^{*},{z}^{*},{x}^{*}})$ be the optimal solution of Eq.
31. $\boldsymbol{Q}|_{\boldsymbol{x}^{*}}$ is the random variable of path
flows. The difference between the model-evaluated system travel time and the
expected system travel time is bounded from above if

  * •

the set of optimal flows is bounded (to implement Lemma 4.5),

  * •

the network flows are bounded from above (i.e., there exist
$\boldsymbol{q}^{\text{Max}}<\infty$ such that
$\hat{\boldsymbol{q}}\leq\boldsymbol{q}^{\text{Max}}\;\forall\;\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})$),
and

  * •

$|\boldsymbol{q}^{*}-\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]|\leq\boldsymbol{\epsilon}$
($\epsilon$-feasibility) and
$\text{Var}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\leq\boldsymbol{\Gamma}^{2}$
($\Gamma$-concentration)

where
$\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]=(\mathbb{E}[{Q}_{i}|_{\boldsymbol{x}^{*}}])_{i\in\mathcal{F}\times\mathcal{T}}$
is the element-wise expectation vector,
$\boldsymbol{\epsilon}=(\epsilon_{i})_{i\in\mathcal{F}\times\mathcal{T}}$,
$\boldsymbol{\Gamma}^{2}=(\Gamma_{i}^{2})_{i\in\mathcal{F}\times\mathcal{T}}$.
The bound of the difference is determined by both $\boldsymbol{\epsilon}$ and
$\boldsymbol{\Gamma}$. Mathematically,

$\displaystyle\big{|}\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}^{*}})]-SST(\boldsymbol{q}^{*})\big{|}\leq
2L\cdot\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}+L\cdot\big{(}\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}+\left\lVert\boldsymbol{q}^{\text{Max}}\right\rVert_{1}+2\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}\big{)}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(40)

###### Remark 4.7

Proposition 4.6 shows that even if the model is optimized on a realization of
the system travel time (not the expectation), as long as we impose the
$\epsilon$-feasibility and $\Gamma$-concentration, the model-evaluated system
travel time and the expected system travel will be similar if
$\boldsymbol{\epsilon}$ and $\boldsymbol{\Gamma}$ are small.

####  4.4.3 Random decision variables in stochastic optimization

The proposed $\epsilon$-feasibility and $\Gamma$-concentration are used to
solve optimization problems with random decision variables. Though typical
stochastic optimization methods generally deal with random parameters (either
in the objective functions or in constraints), the methods may also apply to
cases where decision variables are random. In this section, we analyze a
general optimization problem with random decision variables from the
stochastic optimization point of view, trying to construct the connections
between the proposed two concepts with typical stochastic optimization
methods.

Consider a general optimization problem where the decision variable
$\boldsymbol{Y}=(Y_{i})_{i=1,...,n}\in\mathbb{R}^{n}$ is a random variable
with density function ${f}(\cdot\mid\boldsymbol{\theta})$ (Eq. 41).

$\displaystyle\min_{\boldsymbol{Y}\sim{f}(\cdot\mid\boldsymbol{\theta})}$
$\displaystyle\quad g(\boldsymbol{Y})$ (41a) s.t. $\displaystyle\quad
h_{j}(\boldsymbol{Y})\leq{b}_{j}\quad\forall j\in\mathcal{J}$ (41b)

where $\boldsymbol{\theta}$ is the parameter for the density function.
$g(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the objective function,
$h_{j}(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the constraint function,
$b_{j}\in\mathbb{R}$ is the constraint parameter, $\mathcal{J}$ is the set of
constraints indices. The typical way to transform this problem into a
deterministic problem is to take the expectation of the objective function and
constraints (or consider the probability guarantee of the constraints with a
pre-defined parameter $\eta$, such as $\eta=0.95$), as shown in Eq. 42. That
is, instead of solving for random variable $\boldsymbol{Y}$, we treat the
distribution parameters $\boldsymbol{\theta}$ as the decision variables
(Hernandez 2018), which is deterministic.

$\displaystyle\min_{\boldsymbol{\theta}}$
$\displaystyle\quad\mathbb{E}_{\boldsymbol{Y}\sim{f}(\cdot\mid\boldsymbol{\theta})}[g(\boldsymbol{Y})\mid\boldsymbol{\theta}]$
(42a) s.t.
$\displaystyle\quad\mathbb{E}_{\boldsymbol{Y}\sim{f}(\cdot\mid\boldsymbol{\theta})}[h_{i}(\boldsymbol{Y})]\leq{b}_{j}\quad\forall
j\in\mathcal{J}\qquad\text{(Expectation constraints)}$ (42b) or/and
$\displaystyle\mathbb{P}_{\boldsymbol{Y}\sim{f}(\cdot\mid\boldsymbol{\theta})}[h_{j}(\boldsymbol{Y})\leq{b}_{j}]\geq{\eta}\quad\forall
j\in\mathcal{J}\qquad\text{(Chance constraints)}$ (42c)

However, the formulations in Eq. 42 are in general hard to solve except that
we have the closed-form expressions for
$\mathbb{E}_{\boldsymbol{Y}\sim{f}(\cdot\mid\boldsymbol{\theta})}[\cdot]$ and
$\mathbb{P}_{\boldsymbol{Y}\sim{f}(\cdot\mid\boldsymbol{\theta})}[\cdot]$ (or
using some approximation techniques for the constraints).

In this study, we propose two concepts, $\epsilon$-feasibility and
$\Gamma$-concentration, to model random decision variables. The following
propositions discuss how these two concepts are related to the stochastic
optimization formulation (Eq. 42).

###### Proposition 4.8

The $\epsilon$-feasibility constraint is an approximation to the expectation
constraints (Eq. 42b). Define:

$\displaystyle
G_{\text{SO}}=\min_{\boldsymbol{\theta}}\\{\mathbb{E}[g(\boldsymbol{Y})]:\;\mathbb{E}[h_{j}(\boldsymbol{Y})]\leq{b}_{j},\;\forall
j\in\mathcal{J}\\}$ (43) $\displaystyle
G_{\text{EP}}(\boldsymbol{\epsilon})=\min_{\boldsymbol{y},\boldsymbol{\theta}}\\{g(\boldsymbol{y}):\;h_{j}(\boldsymbol{y})\leq{b}_{j},\;\forall
j\in\mathcal{J},|{\boldsymbol{y}}-\mathbb{E}[\boldsymbol{Y}]|\leq\boldsymbol{\epsilon}\\}$
(44)

where ${\boldsymbol{y}}=(y_{i})_{i=1,...,n}\in\mathbb{R}^{n}$ is a realization
of $\boldsymbol{Y}$.
$\boldsymbol{\epsilon}=(\epsilon_{i})_{i=1,...,n}\in\mathbb{R}^{n}$ is a
vector of small constants (i.e., $\epsilon$-feasibility). All expectations are
taken over $\boldsymbol{Y}\sim{f}(\cdot\mid\boldsymbol{\theta})$.
$\mathbb{E}[\boldsymbol{Y}]=(\mathbb{E}[{Y}_{i}])_{i=1,...,n}$ is the element-
wise expectation vector of $\boldsymbol{Y}$. $G_{\text{SO}}$ is the optimal
solution of the stochastic optimization problem with expectation constraints.
$G_{\text{EP}}$ is the optimal solution of the proposed approach with the
$\epsilon$-feasibility constraint. If $\boldsymbol{\epsilon}=0$, and
$g(\cdot)$ and $h_{j}(\cdot)$ are both convex functions (corresponding to the
convex optimization), we have

$\displaystyle G_{\text{EP}}(\boldsymbol{\epsilon}=0)\leq G_{\text{SO}}$ (45)

That is, the proposed approach is a lower bound of the stochastic optimization
problem.

###### Remark 4.9

The proposition is related to the certainty-equivalent (or mean-field) variant
of a stochastic optimization problem. Consider a special case where both
$g(\cdot)$ and $h_{j}(\cdot)$ are linear, we would have
$G_{\text{EP}}(\boldsymbol{\epsilon}=0)=G_{\text{SO}}$. Note that Eq. 31 is an
integer linear programming. Hence, setting $\boldsymbol{\epsilon}=0$ makes Eq.
31 a special version of the stochastic optimization formulation with
expectation constraints.

###### Proposition 4.10

The $\Gamma$-concentration constraint is an approximation for the chance
constraints (Eq. 42c). Define:

$\displaystyle\tilde{G}_{\text{SO}}=\min_{\boldsymbol{\theta}}\\{\mathbb{E}[g(\boldsymbol{Y})]:\;\mathbb{E}[h_{j}(\boldsymbol{Y})]\leq{b}_{j},\;\mathbb{P}[h_{j}(\boldsymbol{Y})\leq{b}_{j}]\geq{\eta},\;\forall
j\in\mathcal{J}\\}$ (46)
$\displaystyle\tilde{G}_{\text{EP}}=\min_{\boldsymbol{\theta}}\\{\mathbb{E}[g(\boldsymbol{Y})]:\;\mathbb{E}[h_{j}(\boldsymbol{Y})]\leq{b}_{j},\;\forall
j\in\mathcal{J},\text{Var}[\boldsymbol{Y}]\leq\boldsymbol{\Gamma}^{2}\\}$ (47)

where
$\text{Var}[\boldsymbol{Y}]=(\text{Var}[Y_{i}])_{i=1,..,n}\in\mathbb{R}^{n}$
is the element-wise variance vector of $\boldsymbol{Y}$.
$\boldsymbol{\Gamma}^{2}=(\Gamma_{i}^{2})_{i=1,..,n}\in\mathbb{R}^{n}$ is a
vector of squared small constants (i.e., $\Gamma$-concentration). If
$h_{j}(\cdot)$ is Lipschitz continuous (i.e., there exists a positive constant
$C$ such that, for all $\boldsymbol{y_{1}}$ and $\boldsymbol{y_{2}}$,
$|h_{j}(\boldsymbol{y_{1}})-h_{j}(\boldsymbol{y_{2}})|\leq
C\left\lVert\boldsymbol{y_{1}}-\boldsymbol{y_{2}}\right\rVert_{2}$), and
$\boldsymbol{\Gamma}>0$ is sufficiently small, we have

$\displaystyle\text{Var}[\boldsymbol{Y}]\leq\boldsymbol{\Gamma}^{2}\Rightarrow\mathbb{P}[h_{j}(\boldsymbol{Y})\leq{b}_{j}]\geq{\eta}$
(48)

That is, the $\Gamma$-concentration constraints can lead to chance
constraints.

###### Remark 4.11

Proposition 4.10 does not require $h_{j}(\cdot)$ to be a convex function.
Instead, we only need it to be Lipschitz continuous. If $\boldsymbol{\Gamma}$
is not small enough, $\text{Var}[\boldsymbol{Y}]\leq\boldsymbol{\Gamma}^{2}$
cannot result in the chance constraint. We can only obtain a weaker condition
(i.e., looser than the chance constraint). Details can be found in Appendix
15.

####  4.4.4 Connections to robust optimization

Essentially, $\epsilon$-feasibility and $\Gamma$-concentration define a set
for the uncertain path flow $\boldsymbol{Q}$. The uncertainty set is

$\displaystyle\Lambda(\boldsymbol{x})=\\{\boldsymbol{q}\geq
0:|\boldsymbol{q}-\boldsymbol{\mu}(\boldsymbol{x})|\leq\boldsymbol{\epsilon},\text{Var}[\boldsymbol{Q}|_{\boldsymbol{x}}]\leq\boldsymbol{\Gamma}\\}.$
(49)

This set is determined by the endogenous decision variable $\boldsymbol{x}$.
And the optimization is conducted over the uncertainty set
$\Lambda(\boldsymbol{x})$. This concept is related to two-stage robust
optimization with decision-dependent uncertainties (Zeng and Wang 2022), where
the uncertainty set in the second stage is determined by the first-stage
decision variables. However, in our model, the randomness is in the decision
variables rather than in the parameters. Hence, instead of optimizing under
the worst case in $\Lambda(\boldsymbol{x})$, we directly optimize over
$\Lambda(\boldsymbol{x})$ (equivalent to the “best-case”), which avoids the
complexity of min-max formulation in robust optimization.

###  4.5 Solving the problem by Benders decomposition

Eq. 31 is a mixed-integer linear programming (MILP). The structure of Eq. 31
allows us to efficiently solve it by Benders decomposition (BD) (Benders
1962). The basic idea of BD is to decompose the problem into a master problem
and a subproblem and solve these problems iteratively. The decision variables
are divided into difficult variables, which in our case are the binary
variables $\boldsymbol{x}$, and a set of easier variables, the continuous
$\boldsymbol{q}$ and $\boldsymbol{z}$. At each iteration, the master problem
determines one possible leader decision $\boldsymbol{x}$. This solution is
used in the subproblem to generate optimality-cuts or feasibility-cuts, which
are added to the master problem.

Interestingly, in this study, the master problem decides the recommendation
strategies, which is a MILP of a smaller scale and can be solved efficiently
using existing solvers. The subproblem reduces to the optimal flow problem
(Eq. 16) with one more linear constraint (still linear programming). This
format makes the BD an appropriate algorithm for the original problem.

####  4.5.1 Subproblem

The subproblem is derived by fixing the decision variables $\boldsymbol{x}$,
and only considering the components including $\boldsymbol{q}$ and
$\boldsymbol{z}$.

$\displaystyle[SP(\boldsymbol{x})]\quad\min_{\boldsymbol{q},\boldsymbol{z}}\quad$
$\displaystyle WT(\boldsymbol{q},\boldsymbol{z})+IVT(\boldsymbol{z})$ (50a)
s.t. $\displaystyle\text{Constraints }(\ref{const_OF1})-(\ref{const_OF2})$
(50b) $\displaystyle\text{Constraint }(\ref{const_IRR1})$ (50c)

The objective of the dual problem of Eq. 50 is

$\displaystyle
D(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\iota},\boldsymbol{\kappa},\boldsymbol{\rho};\boldsymbol{x})=$
$\displaystyle\sum_{l\in\mathcal{L}}\sum_{t\in\mathcal{T}}\sum_{t^{\prime}=t}^{T_{l,t}}K_{l,t}\alpha_{l,t,t^{\prime}}+\sum_{(u,v,r)\in\mathcal{F}}\sum_{t\in\mathcal{T}}\sum_{t^{\prime}=t^{\text{min}}}^{t}f_{t^{\prime}}^{u,v,r}\beta^{u,v,r}_{t}$
$\displaystyle+\sum_{(u,v,r,i,t)\in\Omega_{1}}\hat{z}_{t}^{u,v,r,i}\gamma_{t}^{u,v,r,i}+\sum_{(u,v)\in\mathcal{W}}\sum_{t\in\mathcal{T}}d^{u,v}_{t}\iota^{u,v}_{t}$
$\displaystyle+\sum_{(u,v,r)\in\mathcal{F}}\sum_{t\in\mathcal{T}}\kappa_{t}^{u,v,r}\cdot(1-\epsilon)\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}$
$\displaystyle+\sum_{(u,v,r)\in\mathcal{F}}\sum_{t\in\mathcal{T}}\rho_{t}^{u,v,r}\cdot(1+\epsilon)\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\pi_{p,r^{\prime}}^{r}$
(51)

where
$\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\iota}$
are the dual variables associated with constraints 16b, 16c, 16f, 16e,
respectively. $\boldsymbol{\kappa}$, $\boldsymbol{\rho}$ are the dual
variables associated with constraint 31c. Let
$\boldsymbol{\Theta}:=(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\iota},\boldsymbol{\kappa}$,
$\boldsymbol{\rho})$. If the dual problem of Eq. 50 is feasible and bounded
with a solution $\boldsymbol{\Theta}^{*}$, the following optimality cut is
added to the master problem:

$\displaystyle Z\geq D(\boldsymbol{\Theta}^{*};\boldsymbol{x})$ (52)

where $Z$ is a decision variable in the master problem. If the dual problem of
Eq. 50 is unbounded, and $\boldsymbol{\Theta}^{*}$ is an optimal extreme ray
of the dual, the following feasibility cut is added to the master problem:

$\displaystyle D(\boldsymbol{\Theta}^{*};\boldsymbol{x})\leq 0$ (53)

####  4.5.2 Master problem

Let $\mathcal{A}^{\text{O}}$ be the set of solutions $\boldsymbol{\Theta}^{*}$
of optimality cuts and $\mathcal{A}^{\text{F}}$ be the set of solutions
$\boldsymbol{\Theta}^{*}$ of feasibility cuts. At each iteration of the BD, a
cut based on the solution of the subproblem is added to the respective set,
and the corresponding master problem is defined as follows:

$\displaystyle[MP(\mathcal{A}^{\text{O}},\mathcal{A}^{\text{F}})]\quad\min_{\boldsymbol{x},Z}$
$\displaystyle\quad\Psi\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}-x_{p,r}\cdot
V_{p,r}+Z$ (54a) s.t. $\displaystyle Z\geq
D(\boldsymbol{\Theta}^{*};\boldsymbol{x})\quad\forall\boldsymbol{\Theta}^{*}\in\mathcal{A}^{\text{O}}$
(54b) $\displaystyle D(\boldsymbol{\Theta}^{*};\boldsymbol{x})\leq
0\quad\forall\boldsymbol{\Theta}^{*}\in\mathcal{A}^{\text{F}}$ (54c)
$\displaystyle\text{Constraints }(\ref{const_Var})-(\ref{const_IRR2})$ (54d)

Note that the master problem has a smaller scale compared to the original
problem (because there are no $\boldsymbol{z}$ and $\boldsymbol{q}$), which
can be solved efficiently.

####  4.5.3 Convergence

Let $(\boldsymbol{x}^{(k)},Z^{(k)})$ and
$(\boldsymbol{q}^{(k)},\boldsymbol{z}^{(k)})$ be the solutions of the master
problem and subproblem, respectively, in the $k$-th iteration. Then, the upper
($UB^{(k)}$) and lower ($LB^{(k)}$) bounds at the $k$-th iteration are given
by:

$\displaystyle UB^{(k)}$
$\displaystyle=\Psi\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}-x_{p,r}^{(k)}\cdot
V_{p,r}+WT(\boldsymbol{q}^{(k)},\boldsymbol{z}^{(k)})+IVT(\boldsymbol{z}^{(k)})$
(55) $\displaystyle LB^{(k)}$
$\displaystyle=\Psi\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}-x_{p,r}^{(k)}\cdot
V_{p,r}+Z^{(k)}$ (56)

$LB^{(k)}$ will keep increasing as $k$ increases because more cuts are added
to the master problem. $UB^{(k)}$ does not necessarily decrease at every
iteration. The convergence criterion is

$\displaystyle\text{Gap}^{(k)}=\frac{UB^{(k)}-LB^{(k)}}{LB^{(k)}}\leq\text{Predetermined
threshold}$ (57)

##  5 Case study

###  5.1 Case study design

####  5.1.1 CTA Blue Line disruption

For the case study, we consider an actual incident in the Blue Line of the
Chicago Transit Authority (CTA) urban rail system (Figure 6). The incident
starts at 8:14 AM and ends at 9:13 AM on Feb 1st, 2019 due to infrastructure
issues between Harlem and Jefferson Park stations (the red X in the figure)
that led to a whole Blue Line suspension. During the disruption (morning
hours), the destination for most of the passengers is the “Loop” in the CBD
area in Chicago. There are four alternative paths to the Loop: 1) using the
Blue Line (i.e., waiting for the system to recover), 2) using the parallel bus
lines, 3) using the North-South (NS) bus lines to transfer to the Green Line,
and 4) using the West-East (WE) bus lines to transfer to the Brown Line. Based
on the service structure, the route sets $\mathcal{R}^{(u,v)}$ for each OD
pair $(u,v)$ can be constructed.

Figure 6: Case study network

In the case study, we divide the time into $\tau=5$ mins equal-length
intervals, and focus on solving the problem at $t=1$ (i.e., beginning of the
incident). We assume that the set of passengers to receive recommendations
($\mathcal{P}$) consists of all passengers with their intended origins at the
Blue Line and destinations in the Loop. A simulation model (Mo et al. 2020) is
used to get the system state up to time $t=1$ (i.e., the incident time 8:14
AM) and generate $\hat{z}_{t}^{u,v,r,i}$ and $\Omega_{1}$. The recommendation
strategy covers passengers departing between $t=1$ and $T^{D}=23$,
approximately one hour after the end of the incident (9:13 AM). The analysis
period is set as $t^{\text{min}}=-13$ and $T=34$, approximately one hour
before $t=1$ and after $T^{D}$, providing enough buffer (warm-up and cool-down
time) for passengers in $\mathcal{P}$ to finish their trips. As demand and
incident duration predictions are out of the scope of this paper, we simply
use the actual demand and incident duration for all experiments. Our other
work (Mo et al. 2022a) proposes to use robust and stochastic optimization to
address demand and incident duration uncertainty, respectively.

####  5.1.2 Conditional probability matrix $\boldsymbol{\pi}$

In this section, we describe how to generate the synthetic conditional
probability matrix $\boldsymbol{\pi}$ used for the case study. During the
incident, CTA does not provide specific path recommendation information. For
every individual, we assume that their actual path choices (referred to as the
“status quo” choices) reflect their inherent preferences. 16 presents the
method of inferring passengers’ status quo choices during the disruption using
smart card data (Mo et al. 2022b). The basic idea is to track their tap-in
records when entering the Blue Line and nearby bus routes, and compare them
with their historical travel histories to get the transfer information.

Given the status quo choices, we assume that the “true” passenger $p$’s
inherent preference for path $r$ is given by

$V_{p}^{r}=\left\\{\begin{aligned} &1+v_{p}^{r}&\quad\text{ if $r$ is $p$'s
actual path choice}\\\ &v_{p}^{r}&\quad\text{
otherwise},\end{aligned}\right.\quad\forall\;p\in\mathcal{P},r\in\mathcal{R}_{p}$
(58)

where $v_{p}^{r}$ is drawn uniformly from $\mathcal{U}[0,1]$. Eq. 58 indicates
every path has a random utility $v_{p}^{r}$ normalized to 0 $\sim$ 1\. And the
chosen path has an additional utility value of 1. We assume that the impact of
the recommendation of $r^{\prime}$ on the utility of path $r$ is

$I_{p,r^{\prime}}^{r}=\left\\{\begin{aligned} &\text{Drawn from
}\mathcal{U}[0,5]&\quad\text{ if $r=r^{\prime}$}\\\ &0&\quad\text{
otherwise},\end{aligned}\right.\quad\forall\;p\in\mathcal{P},r,r^{\prime}\in\mathcal{R}_{p}$
(59)

Eq. 59 means that the utility of the path recommended (i.e., $r=r^{\prime}$)
has an additional positive impact drawn uniformly from $\mathcal{U}[0,5]$. The
utilities of paths not being recommended ($r\neq r^{\prime}$) do not change.
Given Eqs. 58 and 59, we can generate the conditional probability
$\boldsymbol{\pi}$ using Eq. 21. It is worth mentioning that the above
assumptions for generating synthetic passenger prior preferences are based on
two reasonable principles: 1) Passenger’s actual chosen path has a higher
inherent utility. 2) Recommendations of a path can increase its probability of
being chosen.

###  5.2 Parameter settings

The $\epsilon$-feasibility and $\Gamma$-concentration parameters are set as
$\epsilon=0.05$ and $\Gamma=0.3$, indicating 5% and 30% variation constraints
in mean and variance. The convergence gap threshold for Benders decomposition
is set as $1\times 10^{-8}$. The post-adjustment updating step is set as
$\lambda_{k}=\frac{1}{4}$ based on numerical tests.

###  5.3 Benchmark models

There are two benchmark path choice scenarios we use for comparison purposes:

Status-quo path choices. This scenario provides the status quo situation which
does not include any recommendations. It represents the worst case. In this
scenario, no behavior uncertainty is considered because this is based on the
actual path choices realized by passengers.

Capacity-based path recommendations. The capacity-based path recommendations
aim to recommend passengers to different paths according to the available
capacity of paths. Specifically, for a path in OD pair $(u,v)$ and time $t$,
its capacity is the total available capacity of all vehicles passing through
the first boarding station of the path during the time period. For example,
for a path consisting of an NS bus route and the Green Line, the path capacity
is the total available capacity of all buses at the boarding station of the NS
bus route during time interval $t$. The available capacity can be obtained
from a simulation model using historical demand as the input or using
historical passenger counting data. The available capacity for the Blue line
(the incident line) depends on modified operations during the incident (i.e.,
the service suspension is considered). When no vehicles operate in the Blue
line during time interval $t$, the path capacity is zero.

###  5.4 System travel time evaluation

Given a recommendation strategy $\boldsymbol{x}$, as mentioned above, the
actual system travel time is a random variable because of the passenger
behavior uncertainty. To obtain the mean and standard deviation of the system
travel time, we generate multiple passenger choice realizations based on
$\boldsymbol{\pi}$ and $\boldsymbol{x}$. For each generated passenger choice
($\hat{\mathbbm{1}}_{p,r^{\prime}}^{r}$), the realized path flows are

$\displaystyle\hat{q}_{t}^{u,v,r}=\sum_{p\in\mathcal{P}^{u,v}_{t}}\sum_{r^{\prime}\in\mathcal{R}^{u,v}}x_{p,r^{\prime}}\cdot\mathbbm{1}_{p,r^{\prime}}^{r}\quad\forall(u,v,r)\in\mathcal{F},t\in\mathcal{T}.$
(60)

The system travel time for the above passenger choice realization is
calculated by solving the optimal flow problem (Eq. 16) with the constraints
${q}_{t}^{u,v,r}=\hat{q}_{t}^{u,v,r}$ for all $(u,v,r)\in\mathcal{F}$ and
$t\in\mathcal{T}$. This process is repeated with multiple realizations,
providing the sample mean and standard deviation of the system travel time
under recommendation strategy $\boldsymbol{x}$.

###  5.5 Experimental design

As this paper considers various components (such as optimal flow optimization,
passengers’ path preferences, behavior uncertainty, etc.), it is useful to
test different components separately to identify the impact of each one.
Hence, we design the following test cases, each one with specific parameter
settings to systematically evaluate the impacts of each component.

Model performance compared to benchmark models. The most straightforward model
validation is to evaluate the effect of reducing system travel time. In this
test case, we set $\Psi=0$, meaning that we ignore the passengers’ preferences
and focus only on minimizing system travel time. The results of this test case
are discussed in Section 6.2

The benefit of considering behavior uncertainty. In this test case, we
evaluate the importance of incorporating behavior uncertainty in the model.
The model without behavior uncertainty assumes that passengers take the
recommended path. The recommendation strategy is obtained by solving Eq. 31
with $\pi_{p,r^{\prime}}^{r}=1$ if $r=r^{\prime}$. Similarly, we set $\Psi=0$.
Note that, when we evaluate the recommendation strategy, the behavior
uncertainty is still considered in generating the system travel time (see
Section 5.4). The results of this test case are shown in Section 6.3

Impact of considering passenger preference. In all the above tests, $\Psi=0$
is used, focusing on the system travel time. In this test case, we evaluate
the model performance under different values of $\Psi$ in order to assess the
impact of considering passenger preferences. The results of this test case are
discussed in Section 6.4.

##  6 Results

###  6.1 Model convergence

Figure 7 shows the convergence of the BD algorithm. As expected, the lower
bound of the model keeps increasing, while the upper bound, after dropping
significantly in early iterations, exhibits some fluctuations. The model
converges after 28 iterations with a relative gap of less than 1$\times
10^{-8}$. The number of optimality cuts was 28 and no feasibility cut was
generated.

Figure 7: Convergence of the Benders decomposition

Table 1 compares the computational time of the Benders decomposition and off-
the-shelf solvers. The BD algorithm was implemented using Julia 1.6 with the
Gurobi 9.1 solver (Gurobi Optimization, LLC 2021) on a personal computer with
the I9-9900K CPU. The total computational time is 17.8 seconds (master problem
8.2 seconds + subproblem 9.6 seconds), which is more efficient than directly
using the Mixed integer programming (MIP) solvers, including Gurobi (Gurobi
Optimization, LLC 2021), CPLEX (Cplex 2009), GLPK (GNU Linear Programming Kit)
(Makhorin 2008), and CBC (Coin-or branch and cut) (Forrest and Lougee-Heimer
2005).

Table 1: Computational time comparison

Solver | CPU time (sec) | Gap | Solver | CPU time (sec) | Gap
---|---|---|---|---|---
BD | 17.8 | 0.000% | Gurobi | 55.1 | 0.000%
CPLEX | 65.7 | 0.000% | CBC | 425.4 | 0.000%
GLPK | 562.6 | 0.000% |  |  |

###  6.2 Model performance compared to benchmark models

In this section, we compare the system travel time under the proposed
individual path recommendations (without post-adjustment) and two benchmark
models. All travel times (except for the status quo that is deterministic) are
calculated based on 10 replications using the randomly sampled actual path
choices based on the given recommendation (see Section 5.4).

Table 2 shows that the proposed model (IPR) significantly reduces the average
travel time in the system compared to the status quo. Specifically, there is a
6.6% reduction in travel times of all passengers in the system. And for
passengers in the incident line (i.e., passengers who received the
recommendation, $\mathcal{P}$), the average travel time reduction is 19.0%.
Our model also outperforms the capacity-based benchmark path recommendation
strategy, which reduces the travel time of all passengers by 2.5% and incident
line passengers by 15.9%. It is also worth noting that the standard deviation
is small, meaning that variations due to behavior uncertainty are not
significant.

Table 2: Average travel time comparison for different models

Models | Average travel time (all passengers) | Average travel time (incident line passengers)
---|---|---
Mean (min) | Std. (min) | Mean (min) | Std. (min)
Status quo | 28.318 | N.A. | 40.255 | N.A.
Capacity-based | 27.609 (-2.5%) | 0.033 | 33.848 (-15.9%) | 0.165
IPR model | 26.457 (-6.6%) | 0.018 | 32.626 (-19.0%) | 0.187
| Numbers in parentheses represent percentage travel time reduction compared
to the status quo
---

###  6.3 Benefits of considering behavior uncertainty

In this section, we aim to compare the model with and without considering the
behavior uncertainty. The model without behavior uncertainty assumes that all
passengers follow the recommended path when designing the recommendation (but
they may not in reality).

Table 3 shows the comparison of average travel time for the two models. As
expected, considering behavior uncertainty in the path recommendation design
achieves smaller travel time for all passengers and incident line passengers.
Note that, though the 0.93% reduction (around 15 seconds saving per passenger)
is relatively small, considering a large number of passengers in the system,
the total travel time savings are still significant.

Table 3: Average travel time comparison with and without behavior uncertainty
(BU)

Models | Average travel time (all passengers) | Average travel time (incident line passengers)
---|---|---
Mean (min) | Std. (min) | Mean (min) | Std. (min)
IPR model (w.o. BU) | 26.706 | 0.026 | 32.852 | 0.122
IPR model (w. BU) | 26.457 (-0.93%) | 0.018 | 32.626 (-0.69%) | 0.187
| Numbers in parentheses represent percentage travel time reduction compared
to the IPR model w.o. BU
---

###  6.4 Impact of respecting passenger’s prior preferences

In this section, we evaluate the impact of different values of $\Psi$ in terms
of respecting passengers’ prior preferences. Besides the system travel time,
we also evaluate the total utility, defined as the sum of the prior utilities
of the recommended path:

$\displaystyle
TU(\boldsymbol{x})=\sum_{p\in\mathcal{P}}\sum_{r\in\mathcal{R}_{p}}x_{p,r}\cdot
V_{p,r}.$ (61)

Note that the maximum value of $TU(\boldsymbol{x})$ is achieved when every
passenger is recommended with their preferred path (i.e., the path with the
highest prior utility, $V_{p,r}$). Denote this maximum value as
$TU^{\text{max}}$. The relative ratio of total utility,
$\frac{TU(\boldsymbol{x})}{TU^{\text{max}}}$, represents the fraction of the
total (prior) utility that the recommendation has achieved.

Another indicator is the number of passengers recommended with their preferred
path (denoted as $NP(\boldsymbol{x})$). Similarly, we also define the
proportion of passengers recommended with their preferred path (i.e.,
$\frac{NP(\boldsymbol{x})}{|\mathcal{P}|}$, where $|\mathcal{P}|=5,827$ in the
case study).

Figure 8 shows the results for different values of $\Psi$. The x-axis is
plotted in a log-scale. In Figure 8a, the average travel time for all
passengers and incident-line passengers increases with the increase of $\Psi$,
which is as expected because the larger value of $\Psi$ means that the
recommendation generation focuses more on satisfying passenger’s inherent
preferences rather than minimizing the system travel time. Similarly, in
Figure 8b, as expected, both $TU(\boldsymbol{x})$ and $NP(\boldsymbol{x})$
increase with the increase in $\Psi$. When $\Psi=10^{5}$, the average travel
time of the incident line passengers increased by 21.3%, which is close to the
status quo scenario. This is because we generate passengers’ prior utilities
based on the status quo choices. Figure 8b shows that nearly all passengers in
$\mathcal{P}$ are recommended with their preferred path when $\Psi=10^{5}$.

Figure 8 illustrates the trade-off between respecting passengers’ preferences
and reducing system congestion. When the value of $\Psi$ is relatively small
(e.g., less than $10^{3}$), increasing $\Psi$ can effectively increase the
total utility and number of passengers recommended with their preferred path.
Meanwhile, the system travel time only slightly increases. But when $\Psi$ is
large (e.g., greater than $10^{4}$), increasing $\Psi$ significantly increases
the system travel time, but the impact on increasing the passenger’s utility
is limited. The reason may be that, in the system, there are some passengers
whose preferred paths are not at the capacity bottlenecks. Hence, when $\Psi$
is small, the optimal solution recommends those passengers use their preferred
paths without significantly impacting the system travel time. When $\Psi$ is
large, passengers are recommended to use their preferred paths even if these
paths are highly congested, causing a significant increase in the system
travel time. The results imply that a reasonable value of $\Psi$ should be
relatively small. With small $\Psi$, most of the passengers (e.g., more than
70%) are recommended to use their preferred paths without significantly
reducing the system efficiency.

(a) Average travel time

(b) Total utility and number of passengers being recommended with preferred
path

Figure 8: Impact of different values of $\Psi$ on results. The percentage
change in Figure (a) is compared with the scenario of $\Psi=0$. The percentage
in parentheses in Figure (b) represents the relative ratio of total utility
and proportion of passengers recommended with their preferred path,
respectively.

##  7 Conclusion and discussion

This study proposes a mixed-integer programming formulation to model the
individual-based path recommendation problem during PT service disruptions
with the objective of minimizing total system travel time and respecting
passengers’ path choice preferences. Passengers’ behavior uncertainty in path
choices given recommendations is also considered in the formulation. We first
formulate the optimal flow distribution problem in PT systems as linear
programming, which outputs the optimal path flows for each OD pair and time
interval that minimize the total system travel time. Then, we model the
behavior uncertainty based on the passenger’s prior preferences and posterior
path choice probability distribution with two new concepts:
$\epsilon$-feasible flows and $\Gamma$-concentrated flows, which control the
mean and variance of path flows in the optimization problem. We show that
these two concepts can be transformed into linear constraints using
Chebyshev’s inequality. Besides, we show that these two concepts can be seen
as a way of approximating the recourse function (expected system travel time)
in a two-stage stochastic optimization. It is proved that these two concepts
help to bound the difference between the approximated recourse function and
the exact one. Additional theoretical analysis shows that
$\epsilon$-feasibility and $\Gamma$-concentration are approximations of
expectation and chance constraints in a typical stochastic optimization
formulation, respectively. The individual path recommendation problem with
behavior uncertainty is solved using Benders decomposition (BD) efficiently.
The master problem of BD is small-scale integer programming and the subproblem
of BD reduces the optimal flow problem that is a linear program. The BD is
more efficient than many off-the-shelf MIP solvers.

The proposed approach is demonstrated in a case study using data from a real-
world urban rail disruption in the CTA system. The results show that the
proposed IPR model significantly reduces the average travel times in the
system compared to the status quo. Specifically, there is a 6.6% reduction in
travel times for all passengers in the system. Passengers in the incident line
(i.e., passengers who received the recommendation), experience a 19.0% average
travel time reduction. Our model also outperforms the capacity-based benchmark
path recommendation strategy. Compared to the model that assumes all
passengers would follow the recommendations, considering behavior uncertainty
in the path recommendation design can achieve smaller system travel time. In
terms of respecting passengers’ preferences, we show that it is possible that
most of the passengers (e.g., more than 70%) are recommended with their
preferred paths while only increasing the system travel time by 0.51%.

Following the discussion in Section 10, future studies can be pursued in the
following directions. First, as shown in Section 10.1, it is possible to
extend the current framework with more complex recommendation compositions.
The challenges in implementing the more general framework stem from the
quantification of the posterior path choice probabilities. Future studies may
conduct corresponding surveys to calibrate passengers’ responses to the
recommendations. Besides, future studies may consider different sources of
uncertainty (including incident duration, in-vehicle time, demand, etc.) for a
more realistic modeling framework.

##  8 Acknowledgement

The authors would like to thank the Chicago Transit Authority (CTA) for their
support and data availability for this research.

##  9 Notation

Table 4: Notation summary Notation |  Description
---|---
_Model Parameters_
$(u,v,r,i)$ |  The $i$-th leg of path $r$ for OD pair $(u,v)$
$t$ |  Integer time index, $t=1$ represents the start of the incident. Non-positive time indices indicate the time before the incident.
$T^{D}$ |  Time index at which the recommendation system stops working
$t^{\text{min}}$ |  Start time index of the whole analysis period (negative by definition)
$t^{\text{end}}$ |  End time index of the incident
$T$ |  End time index of the whole analysis period (greater than $T^{D}$)
$\mathcal{T}$ |  The set of time indices of analysis and $\mathcal{T}=\\{t^{\text{min}},t^{\text{min}}+1,...,T\\}$
$\tau$ |  The time duration that each time index represents
$\mathcal{P}$ |  The set of passengers that will receive the path recommendation
$\mathcal{R}_{p}$ |  The set of feasible paths for passenger $p\in\mathcal{P}$
$\mathcal{R}^{u,v}$ |  The set of feasible paths for OD pair $(u,v)$
$\Delta^{u,v,r,i}_{t}$ |  Travel time between the terminal station and the boarding station of leg ($(u,v,r,i)$) for vehicle departing at time $t$
$\delta^{u,v,r,i}_{t}$ |  Travel time between the terminal station and the alighting station of leg ($(u,v,r,i)$) for vehicle departing at time $t$
$f_{t}^{u,v,r}$ |  Number of passengers with OD pair ($u,v$) and departure time $t$ and using path $r$ who are not provided path recommendations
$d_{t}^{u,v}$ |  Total number of passengers with OD pair ($u,v$) and departure time $t$
$\Omega_{1}$ |  The set of onboard flow indices at time $t=1$
$\hat{z}^{u,v,r,i}_{t}$ |  Number of onboard passengers in the vehicle departing at time $t$ in leg $(u,v,r,i)$
$\mathcal{P}^{u,v}_{t}$ |  The set of passengers with OD pair $(u,v)$ and departure time $t$ that will receive the path recommendation (a subset of $\mathcal{P}$)
$\mathcal{F}$ |  The set of all $(u,v,r)$ indices
$\mathcal{S}$ |  The set of all stations (stops)
$\mathcal{W}$ |  The set of all OD pairs
$\mathcal{L}$ |  The set of all transit lines (routes) in the system
$\mathcal{I}^{u,v,r}$ |  The set of legs for path $r$ of OD pair $(u,v)$
Vehicle $(l,t)$ |  The vehicle departing at time $t$ on line $l$
$T_{l,t}$ |  The time that vehicle $(l,t)$ arrives the last station of line $l$
$O_{l,t,t^{\prime}}$ |  Total number of onboard passengers at time $t^{\prime}$ for vehicle $(l,t)$
$K_{l,t}$ |  The capacity of vehicles $(l,t)$
$T^{\text{IVT}}_{u,v,r,i,t}$ |  In-vehicle time of leg $(u,v,r,i)$ of vehicle departing at time $t$
$AD_{s,t}$ |  Cumulative arriving demand at station $s$ up to time $t$
$XD_{s,t}$ |  Cumulative transferring demand at station $s$ up to time $t$
$BD_{s,t}$ |  Cumulative boarded demand at station $s$ up to time $t$
$TT^{u,v,r}_{t}$ |  Travel time of a path $(u,v,r)$ at time $t$
$\mathcal{T}^{u,v,r}_{t}$ |  The set of possible arrival times for path $(u,v,r)$ at time $t$
$AT^{u,v,r}_{t}$ |  Arrival time at the destination for the group of passengers using path $(u,v,r)$ and departing at time $t$
$V_{p}^{r}$ |  Passenger $p$’s inherent preference (utility) of using path $r$
$I_{p,r^{\prime}}^{r}$ |  The impact of the recommendation of path $r^{\prime}$ for passenger $p$ on his/her utility of path $r$
$\pi_{p,r^{\prime}}^{r}$ |  Conditional probability for passenger $p$ to choose path $r$ given that he/she is recommended with path $r^{\prime}$
${\mu}_{t}^{u,v,r}$ |  Expectation of $Q^{u,v,r}_{t}$
${\sigma}_{t}^{u,v,r}$ |  Standard deviation of $Q^{u,v,r}_{t}$
$\epsilon_{t}^{u,v,r}$ |  Threshold parameter for $\epsilon$-feasibility
$\Gamma_{t}^{u,v,r}$ |  Threshold parameter for $\Gamma$-concentration
$\Psi$ |  The hyper-parameter to adjust the scale and balance the trade-off between system efficiency and passenger preference
$\mathcal{X}^{\text{OF}}(\hat{\boldsymbol{q}})$ |  The feasible region of the optimal flow problem with an additional constraint ${\boldsymbol{q}}=\hat{\boldsymbol{q}}$
$\mathcal{Q}({\boldsymbol{x}^{*}})$ |  The set of all possible flow patterns given recommendation strategy $\boldsymbol{x}^{*}$
$\eta$ |  Probability guarantee in the chance constraints
$L,C$ |  Parameters for Lipschitz continuity
_Random Variables_
$\mathbbm{1}_{p,r}^{r^{\prime}}$ |  Binary variable indicating whether passenger $p$ has chosen route $r$ or not
$Q^{u,v,r}_{t}$ |  Actual flow for path $(u,v,r)$ at time $t$
_Decision Variables for Optimization Models_
$x_{p,r}$ |  Binary variable indicating whether recommending passenger $p$ to use route $r$ or not
$q^{u,v,r}_{t}$ |  Number of passengers in $\mathcal{P}$ with OD pair ($u,v$) and departure time $t$ and using path $r$
$z^{u,v,r,i}_{t}$ |  Number of passengers boarding a vehicle that had started at time $t$ on leg $(u,v,r,i)$
$y^{u,v,r,t}_{\tilde{t}}$ |  Binary decision variable indicating whether the arrival time of passenger departing at time $t$ using $(u,v,r)$ is $\tilde{t}$ or not
$Z$ |  Decision variable in the master problem of the BD, representing the tentative objective function of the subproblem
$(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\iota},\boldsymbol{\kappa}$, $\boldsymbol{\rho})$ |  Dual variables in the subproblem of the BD

##  10 Model extension

In this section, we discuss several extensions of the model to accommodate
more realistic/general scenarios.

###  10.1 Generalization of recommendations

In this study, we assume the information given to passengers is a recommended
path. In reality, the recommendation system may provide a bundle of
recommended paths with information like estimated in-vehicle time, waiting
time, travel cost, etc. The proposed framework can be extended to handle
different recommendation typologies. Figure 9 shows an example where the
recommendation system will provide a composition of path and travel time
information, where each composition can include different paths, different
estimated waiting/in-vehicle times, etc. Then, we can change $x_{p,r}$ to
$x_{p,c}$, where $x_{p,c}$ indicates whether we will present composition $c$
to passenger $p$. Similarly, each $c$ is associated with a conditional
probability $\pi_{p,c}^{r}$ as shown in Figure 9 (the probability for
passenger $p$ to choose path $r$ given that he/she is recommended composition
$c$).

Figure 9: Illustration of the generalized recommendation typology.
$\mathcal{C}_{p}$ is the predetermined recommendation composition sets for
passenger $p$

In this way, we only need to calibrate $\pi_{p,c}^{r}$ and predetermine the
composition set $\mathcal{C}_{p}$ for each passenger $p$. The overall
framework proposed above can be easily adapted to the new recommendation
typology by replacing $x_{p,r}$ and $\pi_{p,r^{\prime}}^{r}$ with $x_{p,c}$
and $\pi_{p,c}^{r}$, respectively.

###  10.2 Feedback and rolling-horizon

As mentioned in Section 3.2, the whole path recommendation problem should be
solved in a rolling-horizon manner. At each time interval $t\geq 1$, we update
the demand, supply, and system state information, and solve the proposed
framework above to get a recommendation strategy $\boldsymbol{x}$. But we only
implement the $x_{p,r}$ for $p\in\mathcal{P}^{u,v}_{t},\;\forall(u,v)$ (i.e.,
passengers departing at current time $t$).

The rolling horizon requires updating the estimated demand and system state
information. The recommendation system can ask for passenger feedback to
facilitate the estimation. For example, after providing a recommendation, we
can ask the passenger to respond whether he/she will actually use it or not.
This feedback can be used to update the demand predictions.

##  11 Proof of Proposition 4.4

From the triangular inequality, we have:

$\displaystyle\underbrace{\left|Q_{t}^{u,v,r}-{q}_{t}^{u,v,r}\right|}_{\text{LHS}}$
$\displaystyle\leq\left|Q_{t}^{u,v,r}-{\mu}_{t}^{u,v,r}(\boldsymbol{x})\right|+\left|{\mu}_{t}^{u,v,r}(\boldsymbol{x})-{q}_{t}^{u,v,r}\right|$
$\displaystyle\leq\underbrace{\left|Q_{t}^{u,v,r}-{\mu}_{t}^{u,v,r}(\boldsymbol{x})\right|+\epsilon_{t}^{u,v,r}}_{\text{RHS}}$
(62)

As $\text{LHS}\leq\text{RHS}$, the probability measure satisfies (for all
$a>\epsilon_{t}^{u,v,r}$):

$\displaystyle\mathbb{P}[\text{LHS}\geq a]\leq\mathbb{P}[\text{RHS}\geq a]$
(63)

Notice that

$\displaystyle\mathbb{P}[\text{RHS}\geq
a]=\mathbb{P}\left[\left|Q_{t}^{u,v,r}-{\mu}_{t}^{u,v,r}(\boldsymbol{x})\right|\geq
a-\epsilon_{t}^{u,v,r}\right]\leq\frac{({\sigma}_{t}^{u,v,r}(\boldsymbol{x}))^{2}}{(a-\epsilon_{t}^{u,v,r})^{2}}$
(64)

Eq. 64 is based on Chebyshev’s inequality. Therefore,

$\displaystyle\mathbb{P}[\text{LHS}\geq
a]=\mathbb{P}[\left|Q_{t}^{u,v,r}-{q}_{t}^{u,v,r}\right|\geq
a]\leq\frac{({\sigma}_{t}^{u,v,r}(\boldsymbol{x}))^{2}}{(a-\epsilon_{t}^{u,v,r})^{2}}$
(65)

Comparing Eqs. 65 and 26, we know that to satisfy Eq. 26, we only need
${\sigma}_{t}^{u,v,r}(\boldsymbol{x})\leq\Gamma_{t}^{u,v,r}$, which completes
the proof.

##  12 Proof of Lemma 4.5

$STT(\hat{{\boldsymbol{q}}})$ is obtained by solving a linear programming.
$\hat{{\boldsymbol{q}}}$ is a parameter in the constraints. For every
$\hat{{\boldsymbol{q}}}\geq 0$, the problem is feasible because the physical
meaning of the optimal flow problem (i.e., assigning flows to the network) as
long as the system has enough capacity (i.e., dispatching enough vehicles).
Hence, the lemma directly follows Theorem 1 in Martin (1975), which implements
Berge’s Maximum Theorem in parametric linear programming.

##  13 Proof of Proposition 4.6

$\displaystyle\left|\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}^{*}})]-SST(\boldsymbol{q}^{*})\right|=\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})}|SST(\hat{\boldsymbol{q}})-SST(\boldsymbol{q}^{*})|\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle\leq\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})}L\cdot\left\lVert\hat{\boldsymbol{q}}-\boldsymbol{q}^{*}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
(66)

Let us divide the support of the random variable
$\boldsymbol{Q}|_{\boldsymbol{x}^{*}}$ as three mutually exclusive subsets:

$\displaystyle\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}=\mathcal{Q}({\boldsymbol{x}^{*}})\cap\\{\hat{\boldsymbol{q}}:\;0\leq\hat{\boldsymbol{q}}<\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{\epsilon}\\}$
(67)
$\displaystyle\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Mid}}=\mathcal{Q}({\boldsymbol{x}^{*}})\cap\\{\hat{\boldsymbol{q}}:\;\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{\epsilon}\leq\hat{\boldsymbol{q}}\leq\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]+\boldsymbol{\epsilon}\\}$
(68)
$\displaystyle\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Geq}}=\mathcal{Q}({\boldsymbol{x}^{*}})\cap\\{\hat{\boldsymbol{q}}:\;\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]+\boldsymbol{\epsilon}<\hat{\boldsymbol{q}}\leq\boldsymbol{q}^{\text{Max}}\\}$
(69)

where
$\mathcal{Q}({\boldsymbol{x}^{*}})=\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}\cup\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Mid}}\cup\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Geq}}$

We can calculate the summation over these three subsets separately:

(1) Bounds on the summation over
$\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}$:

$\displaystyle\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}L\cdot\left\lVert\hat{\boldsymbol{q}}-\boldsymbol{q}^{*}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle\leq\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}L\cdot\left(\left\lVert\hat{\boldsymbol{q}}-\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}+\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{q}^{*}\right\rVert_{1}\right)\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
(70)

which is followed by the triangle inequality. Notice that

$\displaystyle\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}L\cdot\left\lVert\hat{\boldsymbol{q}}-\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})\leq
L\cdot\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle=L\cdot\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}\cdot\mathbb{P}\big{[}\boldsymbol{Q}|_{\boldsymbol{x}^{*}}\leq\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{\epsilon}\big{]}\leq
L\cdot\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(71)

where the last inequality is the result of the following:

$\displaystyle\mathbb{P}\big{[}\boldsymbol{Q}|_{\boldsymbol{x}^{*}}\leq\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{\epsilon}\big{]}\leq\mathbb{P}\big{[}|\boldsymbol{Q}|_{\boldsymbol{x}^{*}}-\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]|\geq\boldsymbol{\epsilon}\big{]}$
$\displaystyle\leq\sum_{i\in\mathcal{F}\times\mathcal{T}}\mathbb{P}[|{Q_{i}}|_{\boldsymbol{x}^{*}}-\mathbb{E}[{Q_{i}}|_{\boldsymbol{x}^{*}}]|\geq\epsilon_{i}]\leq\sum_{i\in\mathcal{F}\times\mathcal{T}}(\Gamma_{i})^{2}=\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(72)

where the inequality follows by the union bound and the $\Gamma$-concentration
property. Similarly, we have

$\displaystyle\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}L\cdot\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-{\boldsymbol{q}}^{*}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})\leq\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}L\cdot\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle=L\cdot\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}\cdot\mathbb{P}\big{[}\boldsymbol{Q}|_{\boldsymbol{x}^{*}}\leq\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{\epsilon}\big{]}\leq
L\cdot\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(73)

where the first inequality is due to the $\epsilon$-feasibility.

Therefore, combining Eqs. 71 and 73 leads to

$\displaystyle\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}L\cdot\left\lVert\hat{\boldsymbol{q}}-\boldsymbol{q}^{*}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})\leq
L\cdot\big{(}\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}+\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}\big{)}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(74)

(2) Bounds on the summation over
$\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Mid}}$:

$\displaystyle\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Mid}}}L\cdot\left\lVert\hat{\boldsymbol{q}}-\boldsymbol{q}^{*}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle\leq\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}}L\cdot\left(\left\lVert\hat{\boldsymbol{q}}-\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}+\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{q}^{*}\right\rVert_{1}\right)\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle\leq\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Mid}}}L\cdot(\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}+\left\lVert\boldsymbol{\epsilon}\right\rVert_{1})\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})\leq
2L\cdot\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}$ (75)

(3) Bounds on the summation over
$\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Geq}}$:

Similar to the proof of $\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Leq}}$,
notice that

$\displaystyle\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Geq}}}L\cdot\left\lVert\hat{\boldsymbol{q}}-\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})\leq
L\cdot\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Geq}}}\left\lVert\boldsymbol{q}^{\text{Max}}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle=L\cdot\left\lVert\boldsymbol{q}^{\text{Max}}\right\rVert_{1}\cdot\mathbb{P}\big{[}\boldsymbol{Q}|_{\boldsymbol{x}^{*}}\leq\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{\epsilon}\big{]}\leq
L\cdot\left\lVert\boldsymbol{q}^{\text{Max}}\right\rVert_{1}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(76)

Combining Eqs. 76 and 73, we have

$\displaystyle\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Geq}}}L\cdot\left\lVert\hat{\boldsymbol{q}}-\boldsymbol{q}^{*}\right\rVert_{1}\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle\leq\sum_{\hat{\boldsymbol{q}}\in\mathcal{Q}({\boldsymbol{x}^{*}})^{\text{Geq}}}L\cdot\left(\left\lVert\hat{\boldsymbol{q}}-\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}+\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]-\boldsymbol{q}^{*}\right\rVert_{1}\right)\cdot\mathbb{P}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}(\hat{\boldsymbol{q}})$
$\displaystyle\leq
L\cdot\big{(}\left\lVert\boldsymbol{q}^{\text{Max}}\right\rVert_{1}+\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}\big{)}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(77)

In summary, combining the summation over three mutually exclusive sets, we
have:

$\displaystyle\left|\mathbb{E}_{\boldsymbol{Q}|_{\boldsymbol{x}^{*}}}[STT(\boldsymbol{Q}|_{\boldsymbol{x}^{*}})]-SST(\boldsymbol{q}^{*})\right|\leq
2L\cdot\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}+L\cdot\big{(}\left\lVert\mathbb{E}[\boldsymbol{Q}|_{\boldsymbol{x}^{*}}]\right\rVert_{1}+\left\lVert\boldsymbol{q}^{\text{Max}}\right\rVert_{1}+2\left\lVert\boldsymbol{\epsilon}\right\rVert_{1}\big{)}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(78)

##  14 Proof of Proposition 4.8

When $\boldsymbol{\epsilon}=0$, we have
${\boldsymbol{y}}=\mathbb{E}[\boldsymbol{Y}]$. Then:

$\displaystyle
G_{\text{EP}}(\boldsymbol{\epsilon}=0)=\min_{\boldsymbol{y},\boldsymbol{\theta}}\\{g(\mathbb{E}[\boldsymbol{Y}]):\;h_{j}(\mathbb{E}[\boldsymbol{Y}])\leq{b}_{j}\\}$
(79)

According to Jensen’s inequality, we have:

$\displaystyle
g(\mathbb{E}[\boldsymbol{Y}])\leq\mathbb{E}[g(\boldsymbol{Y})],\quad
h_{j}(\mathbb{E}[\boldsymbol{Y}])\leq\mathbb{E}[h_{j}(\boldsymbol{Y})],\;\forall
j\in\mathcal{J}$ (80)

Therefore, the proposed approach has a smaller objective function and a larger
feasible space, which makes it a lower bound of the stochastic optimization
problem (Eq. 43).

##  15 Proof of Proposition 4.10

Step 1: We first show that if $\text{Var}[\boldsymbol{Y}]$ is bounded, then
$\text{Var}[h_{j}(\boldsymbol{Y})]$ is also bounded.

Notice that for any random variable ${X}$, we have
$\text{Var}[{X}]=\mathbb{E}[{X}^{2}]-(\mathbb{E}[{X}])^{2}\leq\mathbb{E}[{X}^{2}]$.
Hence, if we take
${X}=h_{j}(\boldsymbol{Y})-h_{j}(\mathbb{E}[\boldsymbol{Y}])$, we get

$\displaystyle\text{Var}[h_{j}(\boldsymbol{Y})-h_{j}(\mathbb{E}[\boldsymbol{Y}])]=\text{Var}[h_{j}(\boldsymbol{Y})]\leq\mathbb{E}[(h_{j}(\boldsymbol{Y})-h_{j}(\mathbb{E}[\boldsymbol{Y}]))^{2}]$
(81)

From the Lipschitz continuity of $h_{j}(\cdot)$, we have

$\displaystyle|h_{j}(\boldsymbol{Y})-h_{j}(\mathbb{E}[\boldsymbol{Y}])|\leq
C\left\lVert\boldsymbol{Y}-\mathbb{E}[\boldsymbol{Y}]\right\rVert_{2}$ (82)

which further yields:

$\displaystyle\mathbb{E}[(h_{j}(\boldsymbol{Y})-h_{j}(\mathbb{E}[\boldsymbol{Y}]))^{2}]\leq
C^{2}\mathbb{E}[\left\lVert\boldsymbol{Y}-\mathbb{E}[\boldsymbol{Y}]\right\rVert^{2}_{2}]=C^{2}\cdot\sum_{i=1}^{n}\mathbb{E}\left[(Y_{i}-\mathbb{E}[Y_{i}])^{2}\right]=C^{2}\sum_{i=1}^{n}\text{Var}[Y_{i}]$
(83)

Combining with Eq. 81, we have

$\displaystyle\text{Var}[h_{j}(\boldsymbol{Y})]\leq
C^{2}\sum_{i=1}^{n}\text{Var}[Y_{i}]\leq
C^{2}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$ (84)

Step 2: We then show that $\text{Var}[h_{j}(\boldsymbol{Y})]\leq
C^{2}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$ can lead to an
approximation for the chance constraint if $\boldsymbol{\Gamma}$ is
sufficiently small.

Consider Chebyshev’s inequality, for a given positive number $a>0$:

$\displaystyle\mathbb{P}[|h_{j}(\boldsymbol{Y})-\mathbb{E}[h_{j}(\boldsymbol{Y})]|>a]\leq\frac{\text{Var}[h_{j}(\boldsymbol{Y})]}{a^{2}}$
(85)

Eq. 85 implies

$\displaystyle\mathbb{P}[|h_{j}(\boldsymbol{Y})-\mathbb{E}[h_{j}(\boldsymbol{Y})]|\leq
a]\geq
1-\frac{\text{Var}[h_{j}(\boldsymbol{Y})]}{a^{2}}\;\Rightarrow\;\mathbb{P}[h_{j}(\boldsymbol{Y})\leq
a+\mathbb{E}[h_{j}(\boldsymbol{Y})]]\geq
1-\frac{\text{Var}[h_{j}(\boldsymbol{Y})]}{a^{2}}$ (86)

Since we know that $\mathbb{E}[h_{j}(\boldsymbol{Y})]\leq{b}_{j}$, Eq. 86
yields:

$\displaystyle\mathbb{P}[h_{j}(\boldsymbol{Y})\leq a+b_{j}]\geq
1-\frac{\text{Var}[h_{j}(\boldsymbol{Y})]}{a^{2}}\geq
1-\left(\frac{C}{a}\right)^{2}\cdot\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}^{2}$
(87)

Let us pick
$a=\frac{C\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}}{\sqrt{1-\eta}}$, we
have

$\displaystyle\mathbb{P}\left[h_{j}(\boldsymbol{Y})\leq\frac{C\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}}{\sqrt{1-\eta}}+b_{j}\right]\geq\eta$
(88)

Therefore, when $\boldsymbol{\Gamma}$ is sufficiently small, we would have
$\mathbb{P}\left[h_{j}(\boldsymbol{Y})\leq\frac{C\left\lVert\boldsymbol{\Gamma}\right\rVert_{2}}{\sqrt{1-\eta}}+b_{j}\right]\approx\mathbb{P}[h_{j}(\boldsymbol{Y})\leq
b_{j}]$. In this case, we derive the chance constraints from the
$\Gamma$-concentration constraints.

##  16 Inference of status quo choices

The status quo path choice inference method is based on our previous study (Mo
et al. 2022b), which is also similar to the trip-train method used for
destination inference in open public transit systems (i.e., no tap-out).

[In the system when the incident happens]: Consider a passenger
$p\in\mathcal{P}$ with an incident line tap-in record before the end of the
incident, meaning that he/she were in the transit system when the incident
happens. We then track his/her next tap-in record. If he/she next tap-in is a
transfer at a nearby bus or rail station, we can identify his/her chosen path
based on the transfer station. We can also identify the waiting passenger if
he/she continues to use the incident line to his/her intended destination
inferred by his/her next tap-in records.

[Out of the system when the incident happens]: For a passenger
$p\in\mathcal{P}$ with only a tap-in record in nearby bus or rail stations.
He/she may be affected by the incident to change the tap-in station, or just
use the service as a normal commute. To identify whether he/she was affected,
we extract his/her travel histories on previous days without incidents to get
the normal commute trajectories. If his/her tap-in time and location on the
incident day have never appeared in the historical records before, we treat
him/her as a passenger affected by the incident and identify his/her chosen
path based on the tap-in station.

For passengers in $\mathcal{P}$ without next tap-in records or travel
histories, we randomly assign him/her a status quo path based on the
proportion of inferred passengers.

## References

  * Abdelgawad and Abdulhai (2012) Abdelgawad H, Abdulhai B (2012) Large-scale evacuation using subway and bus transit: approach and application in city of toronto. _Journal of Transportation Engineering_ 138(10):1215–1232.
  * Ahmed (2010) Ahmed S (2010) Two-stage stochastic integer programming: A brief introduction. _Wiley encyclopedia of operations research and management science_ 1–10.
  * Ben-Akiva et al. (1985) Ben-Akiva ME, Lerman SR, Lerman SR, et al. (1985) _Discrete choice analysis: theory and application to travel demand_ , volume 9 (MIT press).
  * Benders (1962) Benders JF (1962) Partitioning procedures for solving mixed-variables programming problems. _Numerische mathematik_ 4(1):238–252.
  * Berge (1957) Berge C (1957) Two theorems in graph theory. _Proceedings of the National Academy of Sciences_ 43(9):842–844.
  * Bertsimas et al. (2020) Bertsimas D, Sian Ng Y, Yan J (2020) Joint frequency-setting and pricing optimization on multimodal transit networks at scale. _Transportation Science_ 54(3):839–853.
  * Böhmová et al. (2013) Böhmová K, Mihalák M, Pröger T, Srámek R, Widmayer P (2013) Robust routing in urban public transportation: How to find reliable journeys based on past observations. _ATMOS-13th Workshop on Algorithmic Approaches for Transportation Modelling, Optimization, and Systems-2013_ , volume 33, 27–41 (Schloss Dagstuhl―Leibniz-Zentrum fuer Informatik).
  * Bruglieri et al. (2015) Bruglieri M, Bruschi F, Colorni A, Luè A, Nocerino R, Rana V (2015) A real-time information system for public transport in case of delays and service disruptions. _Transportation Research Procedia_ 10:493–502.
  * Cox et al. (2011) Cox A, Prager F, Rose A (2011) Transportation security and the role of resilience: A foundation for operational metrics. _Transport policy_ 18(2):307–317.
  * Cplex (2009) Cplex II (2009) V12. 1: User’s manual for cplex. _International Business Machines Corporation_ 46(53):157.
  * Forrest and Lougee-Heimer (2005) Forrest J, Lougee-Heimer R (2005) Cbc user guide. _Emerging theory, methods, and applications_ , 257–277 (INFORMS).
  * Guo et al. (2021) Guo X, Caros NS, Zhao J (2021) Robust matching-integrated vehicle rebalancing in ride-hailing system with uncertain demand. _Transportation Research Part B: Methodological_ 150:161–189.
  * Gurobi Optimization, LLC (2021) Gurobi Optimization, LLC (2021) Gurobi Optimizer Reference Manual. URL https://www.gurobi.com.
  * Hernandez (2018) Hernandez H (2018) Introduction to randomistic optimization. _ForsChem Research Reports_ 11\.
  * Horne et al. (2005) Horne M, Jaccard M, Tiedemann K (2005) Improving behavioral realism in hybrid energy-economy models using discrete choice studies of personal transportation decisions. _Energy Economics_ 27(1):59–77.
  * Jung et al. (2004) Jung JY, Blau G, Pekny JF, Reklaitis GV, Eversdyk D (2004) A simulation based optimization approach to supply chain management under demand uncertainty. _Computers & chemical engineering_ 28(10):2087–2106.
  * Khalid et al. (2013) Khalid O, Khan MUS, Khan SU, Zomaya AY (2013) Omnisuggest: A ubiquitous cloud-based context-aware recommendation system for mobile social networks. _IEEE Transactions on Services Computing_ 7(3):401–414.
  * Mahmassani (1984) Mahmassani HS (1984) Uncertainty in transportation systems evaluation: issues and approaches. _Transportation planning and technology_ 9(1):1–12.
  * Majid et al. (2013) Majid A, Chen L, Chen G, Mirza HT, Hussain I, Woodward J (2013) A context-aware personalized travel recommendation system based on geotagged social media data mining. _International Journal of Geographical Information Science_ 27(4):662–684.
  * Makhorin (2008) Makhorin A (2008) Glpk (gnu linear programming kit). _http://www. gnu. org/s/glpk/glpk. html_ .
  * Martin (1975) Martin D (1975) On the continuity of the maximum in parametric linear programming. _Journal of Optimization Theory and Applications_ 17(3):205–210.
  * Mirchevska (2013) Mirchevska V (2013) _Behavior Modeling by Combining Machine Learning and Domain Knowledge_. Ph.D. thesis, PhD Thesis, IPS Jožef Stefan, Ljubljana, Slovenia.
  * Mo et al. (2022a) Mo B, , Koutsopoulos H, Zhao J (2022a) Robust path recommendations during public transit disruptions under demand uncertainty. _Working paper_.
  * Mo et al. (2022b) Mo B, Koutsopoulos H, Zhao J (2022b) Inferring passenger responses to urban rail disruptions using smart card data: A probabilistic framework. _Transportation Research Part E: Logistics and Transportation Review_ 159:102628.
  * 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_ 0361198120914309\.
  * Mo et al. (2022c) Mo B, Von Franque MY, Koutsopoulos HN, Attanucci JP, Zhao J (2022c) Impact of unplanned long-term service disruptions on urban public transit systems. _IEEE Open Journal of Intelligent Transportation Systems_ 3:551–569, URL http://dx.doi.org/10.1109/OJITS.2022.3199108.
  * Mo et al. (2021) Mo B, Wang QY, Moody J, Shen Y, Zhao J (2021) Impacts of subjective evaluations and inertia from existing travel modes on adoption of autonomous mobility-on-demand. _Transportation Research Part C: Emerging Technologies_ 130:103281.
  * Richardson et al. (2007) Richardson M, Dominowska E, Ragno R (2007) Predicting clicks: estimating the click-through rate for new ads. _Proceedings of the 16th international conference on World Wide Web_ , 521–530.
  * Roelofsen et al. (2018) Roelofsen D, Cats O, van Oort N, Hoogendoorn S (2018) Assessing disruption management strategies in rail-bound urban public transport systems from a passenger perspective. _Proceedings of the 14th Conference on Advanced Systems in Public Transport (CASPT), Brisbane, Australia_.
  * Subramanyam et al. (2021) Subramanyam A, Mufalli F, Laínez-Aguirre JM, Pinto JM, Gounaris CE (2021) Robust multiperiod vehicle routing under customer order uncertainty. _Operations Research_ 69(1):30–60.
  * Tan et al. (2020) Tan Z, Xu M, Meng Q, Li ZC (2020) Evacuating metro passengers via the urban bus system under uncertain disruption recovery time and heterogeneous risk-taking behaviour. _Transportation research part C: emerging technologies_ 119:102761.
  * Train (2009) Train KE (2009) _Discrete choice methods with simulation_ (Cambridge university press).
  * Wang et al. (2019) Wang J, Yuan Z, Yin Y (2019) Optimization of bus bridging service under unexpected metro disruptions with dynamic passenger flows. _Journal of Advanced Transportation_ 2019\.
  * Wang et al. (2020) Wang S, Mo B, Zhao J (2020) Deep neural networks for choice analysis: Architecture design with alternative-specific utility functions. _Transportation Research Part C: Emerging Technologies_ 112:234–251.
  * Yeung and Yang (2010) Yeung KF, Yang Y (2010) A proactive personalized mobile news recommendation system. _2010 Developments in E-systems Engineering_ , 207–212 (IEEE).
  * Yoon and Chow (2020) Yoon G, Chow JY (2020) Contextual bandit-based sequential transit route design under demand uncertainty. _Transportation Research Record_ 2674(5):613–625.
  * Zeng and Wang (2022) Zeng B, Wang W (2022) Two-stage robust optimization with decision dependent uncertainty. _arXiv preprint arXiv:2203.16484_ .