# An exact upper bound for the minimum size of a path system that weakly
separates a clique

George Kontogeorgiou and Maya Stein Center for Mathematical Modeling (CNRS
IRL2807), University of Chile. Supported by ANID Basal Grant CMM FB210005.
Department of Mathematical Engineering and Center for Mathematical Modeling
(CNRS IRL2807), University of Chile. Supported by FONDECYT Regular Grant
1221905, by ANID Basal Grant CMM FB210005, and by MSCA-RISE-2020-101007705
project RandNET.

(March 2024)

###### Abstract

We show that for each natural number $n$ the clique $K_{n}$ has a weakly
separating path system of size $n+2$, improving the previous best known bound
of $(1+o(1))n$. Since $n-1$ is a lower bound, this is almost tight.

##  1 Introduction

Let $G$ be a graph. A _path system_ of $G$ is a set of paths in $G$. A path
system $\mathcal{P}$ of $G$ is _weakly separating_ , or simply, _separating_ ,
if for every pair of edges of $G$ there exists a path in $\mathcal{P}$ that
contains exactly one of them. Let $f:G\rightarrow\mathbb{N}$ be the minimum
size of a path system that separates $G$. We are interested in the problem of
determining $f$ for complete graphs.

The study of general separating set systems was initiated by Rényi in the
1960s [8]. Various versions concerning the separation of edges by subgraphs
were introduced around the 2000s by computer scientists who were seeking
efficient ways to detect faulty links in networks [4, 5, 9, 11]. The problem
of determining a tight upper bound for the size of separating path systems of
graphs on $n$ vertices (i.e. determining $\min\\{f(G):G\text{ is an $n$-vertex
graph}\\}$) has become very popular in the combinatorics community during the
last decade.

Among the first to consider this problem were Falgas-Ravry, Kittiparson,
Korándi, Letzter and Narayanan [2], who proved a bound of $O(n\log n)$ and
conjectured one of $O(n)$. Letzter [6] came very close by proving that
$O(n\log^{*}n)$ is sufficient. The conjecture was eventually settled by
Bonamy, Botler, Dross, Naia and Skokan [1], who proved an upper bound of
$19n$. They further conjectured that $(1+o(1))n$ is optimal.

Restricting the question to the class of complete graphs is an appropriate
direction to follow, not only because of the reduced complexity but also
because of the fact that complete graphs have the highest known lower bound
for $f$ among all graph classes. It is therefore natural to suspect that they
are extremal for this problem. Specifically, the lower bound $f(K_{n})\geq
n-1$ is tight and easy to obtain [2]. As for the upper bound, owing to work by
Fernandes, Oliveira Mota and Sanhueza-Matamala [3], it is known that

$f(K_{n})=(1+o(1))n.$

The authors of [3] employ hypergraph matching theory in their proof which
allows them to prove their bound for all dense regular graphs. Our approach to
the problem is very different, and more alike to that of Wickes [10]. Her work
relies on _generating paths_ (for a definition see below). Namely, Wickes
showed that every complete graph on a prime number of vertices has a
generating path, and that the set of all the rotations of a generating path is
a separating path system. Hence, $f(K_{n})\leq n$ when $n$ is prime. She also
proved that $f(K_{n})\leq n$ when $n-1$ is prime or when $n\leq 20$.

We give a similar upper bound for $f(K_{n})$ for all $n\in\mathbb{N}$.

######  Theorem 1.

For every natural number $n$, $f(K_{n})\leq n+2$.

In combination with the lower bound of $n-1$, our result nearly resolves the
problem of the size of minimum separating path systems for complete graphs.
Moreover, our proof is both short and entirely constructive, providing
concrete and relatively easy to describe examples of separating path systems
for all complete graphs.

In brief, we decompose $n$ into appropriate prime numbers $p_{k}$ and a small
rest, and then apply Wicke’s construction to each of the $K_{p_{k}}$’s. Then
begins the main work of our proof, which is to carefully modify the resulting
path system so that it becomes separating and has the required number of
paths.

##  2 Definitions and Preliminaries

A _labelled clique_ is a complete graph $K_{n}$ together with a labelling of
$V(K_{n})$ with numbers from an interval of the naturals of length $n$. The
latter naturally induces on $V(K_{n})$ a cyclic ordering and a rotation
action. Given a labelled clique $K_{n}$, we partition $E(K_{n})$ in _types_ as
follows: the type of an edge $\\{x_{1},x_{2}\\}$ is defined to be
$\min\\{z_{1},z_{2}\\}$, where $z_{i}:=(x_{i}-x_{3-i})(\text{mod
$n$})\in[\frac{n-1}{2}]$ for $i=1,2$. Edges of the same type will be called
_homotypical_. It is also useful to define the _clockwise distance_ of two
homotypical edges contained in a labelled clique as the minimum order of a
rotation that maps one to the other.

For our immediate purposes, it suffices to define a generating path $P$ of
$K_{p}$, where $p$ is a prime number, as one that exhibits the following
properties:

  1. (GP1)

$P$ contains exactly one edge of type $1$;

  2. (GP2)

$P$ contains exactly two edges of each other type;

  3. (GP3)

any two pairs of homotypical edges in $P$ have distinct clockwise distances.

From these properties it is also clear that in the path system generated by
the rotations of $P$ each edge of type $1$ is contained in exactly one path,
whereas every other edge is contained in exactly two paths.

In the following proof, we work under the convention that every path has a
first and a last vertex (although which is which will often not be relevant).
We denote $e(P,P^{\prime})$ the edge connecting the last vertex of $P$ to the
first vertex of $P^{\prime}$.

##  3 Proof

Let $3\leq n=\sum_{k=1}^{m}p_{k}+b$ with $p_{k}$ primes. For brevity, let
$s_{k}$ denote the partial sums of $\sum_{k=1}^{m}p_{k}$. It is a corollary of
[7] that we can choose each $p_{k}$ so that $n-s_{k}\leq\frac{p_{k}-3}{2}$ and
$b\in\\{3,...,19\\}$.

We define pairwise disjoint induced cliques of $K_{n}$ indexed by, and with
sizes corresponding to, the terms of the above sum. We label the vertices of
each clique $K_{p_{k}}$ using the interval $[s_{k-1}+1,s_{k}]$. We call these
labelled cliques _levels_ , and we say that an edge _touches_ a level if its
smaller end is contained in that level. Note that exactly the edges that are
contained in a level have a type in it. Therefore, whenever we refer to the
type of an edge, it will always be an edge contained in a level and we will
always mean its type in that level.

Let

$\mathscr{B}:=\\{\mathbf{B_{m+1}^{s_{m}+1}},...,\mathbf{B_{m+1}^{n}}\\}$

be a separating path system for $K_{b}$.

Let $P_{k}^{s_{k-1}+1}$ be a generating path for the clique $K_{p_{k}}$. Let

$\mathcal{P}_{k}:=\\{P_{k}^{s_{k-1}+1},...,P_{k}^{s_{k}}\\}$

be the separating path system of $K_{p_{k}}$ generated by $P_{k}^{s_{k-1}+1}$.
Let $s_{k-1}+i$ be the first vertex of $P_{k}^{s_{k-1}+i}$.

We construct a separating path system $\mathscr{Q}$ for $K_{n}$ as follows.

To begin, we note that each generating path $P_{k}^{s_{k-1}+1}$ defines a
permutation $\sigma_{k}\in S_{\frac{p_{k}-1}{2}}$ that maps each number in
$\\{2,...,\frac{p_{k}-1}{2}\\}$ to the clockwise distance in $K_{p_{k}}$ of
the homotypical edges of $P_{k}^{s_{k-1}+1}$ that have it as their type, and
maps $1$ to the unique number in $[\frac{p_{k}-1}{2}]$ that is not the
clockwise distance in $K_{p_{k}}$ of any homotypical pair of edges in
$P_{k}^{s_{k-1}+1}$. Based on this observation, for each $p_{k}$ we define an
injection
$\phi_{k}:\\{s_{k}+1,...,n\\}\rightarrow\\{2,...,\frac{p_{k}-1}{2}\\}$ in the
following way. While there does not exist a cycle of $\sigma_{k}$ with at
least as many elements distinct from $1$ as there are remaining undefined
instances of $\phi_{k}$, we arbitrarily choose a cycle of $\sigma_{k}$ that
does not contain $1$ and arbitrarily decide the preimages of all of its
elements, avoiding $s_{k}+1$. Then, we choose a cycle of $\sigma_{k}$ with at
least as many elements distinct from $1$ as there are remaining vertices in
$\\{s_{k}+1,...,n\\}$ with undefined images under $\phi_{k}$, and we map all
of these remaining vertices to consecutive elements of the chosen cycle
distinct from $1$, starting with $s_{k}+1$. This concludes the definition of
$\phi_{k}$.

For each path $P_{k}^{s_{k-1}+i}\in\mathcal{P}_{k}$ and each type
$\phi_{k}(s_{k}+j)\in im(\phi_{k})$, we remove exactly one edge of type
$\phi_{k}(s_{k}+j)$ in $P_{k}^{s_{k-1}+i}$, different for each $i$, and join
the ends of that edge to $s_{k}+j$. For example, given an edge of type
$\phi_{k}(s_{k}+j)$, say $\\{s_{k-1}+x,s_{k-1}+(x+\phi_{k}(s_{k}+j)\text{(mod
$p_{k}$))}\\}$, that is contained in $P_{k}^{s_{k-1}+1}$, then in each
$P_{k}^{s_{k-1}+i}$ we can apply the treatment described above to the edge
$\\{s_{k-1}+(x+i-1\text{(mod
$p_{k}$)}),s_{k-1}+(x+i-1+\phi_{k}(s_{k}+j)\text{(mod $p_{k}$)})\\}$. We name
the resulting path $\mathbf{P_{k}^{s_{k-1}+i}}$, and define

$\mathscr{P}_{k}:=\\{\mathbf{P_{k}^{s_{k-1}+1}},...,\mathbf{P_{k}^{s_{k}}}\\}.$

It is worth noting that every path
$\mathbf{P_{k}^{s_{k-1}+i}}\in\mathscr{P}_{k}$:

  * (P1)

has exactly one edge in level $k$ of each type in $im(\phi_{k})\cup\\{1\\}$,
which appears in no other path of $\mathscr{P}_{k}$ \- in particular, this
implies that each such edge is contained in exactly one path in
$\mathscr{P}_{k}$;

  * (P2)

has exactly two edges in level $k$ of each other type, which it inherits from
$P_{k}^{s_{k-1}+i}$ \- in particular, this implies that every such edge is
contained in exactly two paths of $\mathscr{P}_{k}$;

  * (P3)

has exactly two edges incident to each vertex $s_{k}+j$, one of which it
shares with $\mathbf{P_{k}^{s_{k-1}+(i+\phi_{k}(s_{k}+j)\text{(mod $p$)})}}$
and the other with $\mathbf{P_{k}^{s_{k-1}+(i-\phi_{k}(s_{k}+j)\text{(mod
$p$)})}}$.

We also observe that, for each vertex $s_{k}+j$, there do not exist
homotypical edges in any path of $\mathscr{P}_{k}$ that have clockwise
distance $\phi_{k}(s_{k}+j)$ in $K_{p_{k}}$, since in every such path one edge
of each such type has been substituted by a path of length $2$ passing through
another vertex in $\\{s_{k}+1,...,n\\}$ as described above - we denote this
property (P4). The only possible exception to (P4) is for the vertex
$s_{k}+1$, if the last cycle of $\sigma_{k}$ that we picked in the process of
defining $\phi_{k}$ does not lie entirely in $im(\phi_{k})$ and does not map
$1$ to $\phi_{k}(s_{k}+1)$. In that case, every path in $\mathscr{P}_{k}$
contains two homotypical edges, of type say $c_{k}\notin
im(\phi_{k})\cup\\{1\\}$, that have clockwise distance $\phi_{k}(s_{k}+1)$ in
$K_{p_{k}}$.

We continue by constructing sets of paths

$\mathcal{R}_{k}:=\\{R_{k}^{s_{k}+1},...,R_{k}^{n}\\},$

where $R_{k}^{s_{k}+j}$ contains all the edges of $K_{p_{k}}$ of type
$\phi_{k}(s_{k}+j)$ but one. We also define a path $S_{k}$ that contains all
the edges of type $1$ but one, and a path $T_{k}$ that contains all the edges
of type $c_{k}$ but one. Specifically, we pick the one edge excluded by each
path in $\mathcal{R}_{k}\cup\\{S_{k},T_{k}\\}$ as follows:

  * (P5)

for $\mathcal{R}_{k}\cup\\{S_{k}\\}$, we pick them from pairwise distinct
paths of $\mathscr{P}_{k}$;

  * (P6)

for $T_{k}$, we pick the one that is at the same two paths of
$\mathscr{P}_{k}$ as $e(R_{k}^{s_{k}+1},P_{k+1}^{s_{k}+1})$ (or
$e(R_{m}^{s_{m}+1},\mathbf{B_{m+1}^{s_{m}+1}})$, if $k=m$), if such exists -
otherwise we choose arbitrarily.

We may impose (P5) because each path of $\mathscr{P}_{k}$ has exactly one edge
of each type in $im(\phi)\cup\\{1\\}$ and it is different from those of all
other paths in $\mathscr{P}_{k}$, due to (P1). We may impose (P6) because each
edge incident to $s_{k}+1$ in a path of $\mathscr{P}_{k}$, including
$e(R_{k}^{s_{k}+1},P_{k+1}^{s_{k}+1})$ (or
$e(R_{m}^{s_{m}+1},\mathbf{B_{m+1}^{s_{m}+1}})$), is contained in the same two
paths of $\mathscr{P}_{k}$ as a different unique edge of type $c_{k}$, if
$s_{k}+1$ is an exception to (P4), or as none of the edges of type $c_{k}$, if
(P4) applies to $s_{k}+1$.

We then extend each path $R_{k}^{s_{k}+j}$ by the edge
$e(R_{k}^{s_{k}+j},P_{k+1}^{s_{k}+j})$ if $k\leq m-1$ and $j\leq p_{k+1}$, or
by the edge $e(R_{k}^{s_{k}+j},R_{k+1}^{s_{k}+j})$ if $k\leq m-1$ and $j\geq
p_{k+1}+1$, or by the edge $e(R_{m}^{s_{m}+j},\mathbf{B_{m+1}^{s_{m}+j})}$ if
$k=m$, to obtain a path $\mathbf{R_{k}^{s_{k}+j}}$. We define

$\mathscr{R}_{k}:=\\{\mathbf{R_{k}^{s_{k}+1}},...,\mathbf{R_{k}^{n}}\\}.$

For $k\leq m-1$ we also extend $S_{k}$ by $e(S_{k},S_{k+1})$ and $T_{k}$ by
$e(T_{k},T_{k+1})$ to obtain paths $\mathbf{S_{k}}$ and $\mathbf{T_{k}}$,
respectively. Exploiting the fact that each path has two ends, we impose that
the starting vertex of $T_{k+1}$ be not $s_{k}+1$, i.e. $e(T_{k},T_{k+1})$ is
not incident to $s_{k}+1$ (P7).

Finally, we define

$\mathscr{Q}:=\\{\mathbf{Q^{1}},...,\mathbf{Q^{n}}\\}\cup\\{\mathbf{S},\mathbf{T}\\},$

where:

  * •

$\mathbf{S}$ is the union of all the paths $\mathbf{S_{k}}$;

  * •

$\mathbf{T}$ is the union of all the paths $\mathbf{T_{k}}$;

  * •

$\mathbf{Q^{\ell}}$ is the union of all the boldface paths of exponent $\ell$,
i.e. $\mathbf{Q^{\ell}}$ is the union of all paths $\mathbf{R^{\ell}_{i}}$
that were defined, and of $\mathbf{P^{\ell}_{k}}$ or $\mathbf{B^{\ell}_{m+1}}$
if these were defined, where $k$ is the appropriate index.

The reader can check that $\mathscr{Q}$ is indeed a path system, as its
elements are concatenations of pairwise internally disjoint paths. It remains
to show that it is separating.

We begin by noting that $\mathbf{Q^{1}}\cup...\cup\mathbf{Q^{\ell}}$, where
$\ell=s_{k}$, contains exactly the edges that touch one of the first $k$
levels. Therefore, edges that touch distinct levels are separated. Moreover,
edges in level $m+1$ are separated by $\mathscr{B}$. On the other hand, all
the edges that touch a level $k\leq m$ are contained in the paths of
$\mathscr{Q}_{k}:=\mathscr{P}_{k}\cup\mathscr{R}_{k}\cup\\{\mathbf{S_{k}},\mathbf{T_{k}}\\}$.
Moreover, any path $P\in\mathscr{Q}_{k}$ extends to a path $P\in\mathscr{Q}$
such that $P^{\prime}\cap\mathscr{Q}_{k}=P$. It is therefore sufficient to
show that $\mathscr{Q}_{k}$ separates the edges that touch level $k$.

Let $e,e^{\prime}$ be two distinct edges that touch level $k$ but are not
contained in it. Then the only paths that may contain them both are those of
$\mathscr{P}_{k}$, since every other path of $\mathscr{Q}_{k}$ contains at
most one such edge. So, suppose that $\mathbf{P_{k}^{s_{k-1}+i}}$ is a path
containing both $e$ and $e^{\prime}$.

If $e$ and $e^{\prime}$ have a common greater end, say $s_{k}+j$, then by (P3)
the only other path of $\mathscr{P}_{k}$ that contains, say, $e$ is
$\mathbf{P_{k}^{s_{k-1}+(i-\phi_{k}(s_{k}+j)\text{(mod $p_{k}$)})}}$, and the
only other path of $\mathscr{P}_{k}$ that contains $e^{\prime}$ is
$\mathbf{P_{k}^{s_{k-1}+(i+\phi_{k}(s_{k}+j)\text{(mod $p_{k}$)})}}$. But
these are distinct paths, since $2\phi_{k}(s_{k}+j)\neq 0\text{(mod
$p_{k}$)}$. So $e$ and $e^{\prime}$ are separated.

If $e$ and $e^{\prime}$ have distinct greater ends, say $s_{k}+j_{1}$ and
$s_{k}+j_{2}$, then again by (P3) the only other path of $\mathscr{P}_{k}$
that contains $e$ is
$\mathbf{P_{k}^{s_{k-1}+(i\pm\phi_{k}(s_{k}+j_{1})\text{(mod $p_{k}$)})}}$,
and the only other path of $\mathscr{P}_{k}$ that contains $e^{\prime}$ is
$\mathbf{P_{k}^{s_{k-1}+(i\pm\phi_{k}(s_{k}+j_{2})\text{(mod $p_{k}$)})}}$.
Again, since $\phi_{k}(s_{k}+j_{1})$ and $\phi_{k}(s_{k}+j_{2})$ are distinct
elements of $\\{2,...,\frac{p_{k}-1}{2}\\}$, these two paths are distinct,
separating $e$ and $e^{\prime}$.

Hence, every two edges that touch level $k$ but are not contained in it are
separated.

Subsequently, let $e$ and $e^{\prime}$ both be contained in level $k$. If one
of them is contained in one of the paths in
$\mathscr{R}_{k}\cup\\{\mathbf{S_{k}}\\}$ and the other is not contained in
the same path, then of course they are separated. For this to not happen, one
of the following must occur.

It may be that $e$ and $e^{\prime}$ are both contained in the same path of
$\mathscr{R}_{k}\cup\\{\mathbf{S_{k}}\\}$. But then they are homotypical, and
their type is in $im(\phi_{k})\cup\\{1\\}$. By (P1), that means that each of
$e$ and $e^{\prime}$ belongs to a different path of $\mathscr{P}_{k}$, so they
are separated.

It may be that $e$ and $e^{\prime}$ are not contained in any of the paths in
$\mathscr{R}_{k}\cup\\{\mathbf{S_{k}}\\}$ because their types are not in
$im(\phi_{k})\cup\\{1\\}$. But then by (P2) each of them appears in exactly
two paths of $\mathscr{P}_{k}$, corresponding to the two paths in which it
appears in $\mathcal{P}_{k}$, therefore they are separated.

It may be that $e$ is the unique edge of some type in
$im(\phi_{k})\cup\\{1\\}$ that is not contained in any paths in
$\mathscr{R}_{k}\cup\\{\mathbf{S_{k}}\\}$. If the same holds for $e^{\prime}$,
then, by (P5), $e$ and $e^{\prime}$ are still separated by a path in
$\mathscr{P}_{k}$. If instead $e^{\prime}$ is of a type not in
$im(\phi_{k})\cup\\{1\\}$, then by (P2) $e^{\prime}$ appears in two paths in
$\mathscr{P}_{k}$ and by (P1) $e$ appears only in one, so again they are
separated.

Therefore, every two edges contained in level $k$ are separated.

Finally, let $e$ be contained in level $k$ and let $e^{\prime}$ touch level
$k$ but not be contained in it. Let $s_{k}+j$ be the greater end of
$e^{\prime}$. Let $\mathbf{P_{k}^{s_{k-1}+i}}$ be a path of $\mathscr{P}_{k}$
containing $e^{\prime}$, and suppose that it also contains $e$. Then by (P3)
and without loss of generality $e^{\prime}$ is also contained in
$\mathbf{P_{k}^{s_{k-1}+(i+\phi_{k}(s_{k}+j)\text{(mod $p_{k}$)})}}$. If
$j\neq 1$, then by (P4) there do not exist homotypical edges in any path of
$\mathscr{P}_{k}$ that have clockwise distance $\phi_{k}(s_{k}+j)$ in
$K_{p_{k}}$, so $e\notin\mathbf{P_{k}^{s_{k-1}+(i+\phi_{k}(s_{k}+j)\text{(mod
$p_{k}$)})}}$ and the two edges are separated. The same holds if $j=1$ and
either $s_{k}+1$ is not an exception to $\mathbf{(P4)}$ or $e$ is not of type
$c_{k}$. If $j=1$, $e$ is of type $c_{k}$ and $e\in\mathbf{T_{k}}$, then again
$e$ and $e^{\prime}$ are separated by $T_{k}$, since
$e^{\prime}\notin\mathbf{T_{k}}$ due to (P7). Finally, if $j=1$, $s_{k}+1$ is
an exception to $\mathbf{(P4)}$ and $e$ is the unique edge of type $c_{k}$ not
in $\mathbf{T_{k}}$, then by (P6) either $e$ and $e^{\prime}$ are separated by
a path in $\mathscr{P}_{k}$, or else $e^{\prime}\in\mathbf{R_{k}^{s_{k}+1}}$,
so $e$ and $e^{\prime}$ are again separated.

This concludes the case analysis, and the proof.

## References

  * [1] Marthe Bonamy, Fábio Botler, François Dross, Tássio Naia, and Jozef Skokan.  Separating the edges of a graph by a linear number of paths.  Advances in Combinatorics.  https://doi.org/10.19086/aic.2023.6.
  * [2] Victor Falgas-Ravry, Teeradej Kittipassorn, Daniel Korándi, Shoham Letzter, and Bhargav P. Narayanan.  Separating path systems.  J. Comb., 5(3):335–354, 2014.
  * [3] Cristina G. Fernandes, Guilherme Oliveira Mota, and Nicolas Sanhueza-Matamala.  Separating path systems in complete graphs.  https://arxiv.org/abs/2312.14879, 2023.
  * [4] Nicholas J. A. Harvey, Mihai Patrascu, Yonggang Wen, Sergey Yekhanin, and Vincent W. S. Chan.  Non-adaptive fault diagnosis for all-optical networks via combinatorial group testing on graphs.  IEEE INFOCOM, pages 697–705, 2007.
  * [5] Iiro Honkala, Mark Karpovsky, and Simon Litsyn.  Cycles identifying vertices and edges in binary hypercubes and 2-dimensional tori.  Discrete Appl. Math., 129(2-3):409–419, 2003.
  * [6] Shoham Letzter.  Separating paths systems of almost linear size.  https://arxiv.org/abs/2211.07732, 2022.
  * [7] Jitsuro Nagura.  On the interval containing at least one prime number.  Proceedings of the Japan Academy, Series A, 28(4):177–181, 1952.
  * [8] Alfred Renyi.  On random generating elements of a finite boolean algebra.  Acta Sci. Math. (Szeged), (22):75–81, 1961.
  * [9] Janos Tapolcai, Lajos Ronyai, and Pin-Han Ho.  Link fault localization using bi-directional m-trails in all-optical mesh networks.  IEEE Transactions on Communications, 61(1):291–300, 2013.
  * [10] Belinda Wickes.  Separating path systems for the complete graph.  Discrete Math., 347(3):113784, 2023.
  * [11] Lev Zakrevski and Mark Karpovsky.  Fault-tolerant message routing for multiprocessors.  Lecture Notes in Computer Science, International Parallel Processing Symposium, 1388:714–730, 1998.


A branch-cut-and-price algorithm for the dial-a-ride problem with minimum
disease-transmission risk

Shuocheng Guo Department of Civil, Construction and Environmental Engineering,
The University of Alabama, Tuscaloosa, AL 35487<EMAIL_ADDRESS>Iman Dayarian
Culverhouse College of Business, The University of Alabama, Tuscaloosa, AL
35487<EMAIL_ADDRESS>Jian Li School of Transportation Engineering,
Tongji University, Shanghai, China<EMAIL_ADDRESS>Xinwu Qian Department
of Civil, Construction and Environmental Engineering, The University of
Alabama, Tuscaloosa, AL 35487<EMAIL_ADDRESS>
This paper investigates a novel variant of the dial-a-ride problem (DARP),
namely Risk-aware DARP (RDARP), which focuses on the bi-objective DARP by
minimizing the total routing cost and the maximum exposed risk for onboard
passengers. The RDARP requires that the exposed risk of each passenger is
minimized along the trip while satisfying existing constraints of the DARP. We
model the RDARP using a three-index arc-based formulation and reformulate it
into an equivalent min-max trip-based formulation, which is solved optimally
using a tailored Branch-Cut-and-Price (BCP) algorithm. The BCP algorithm
adopts the Column Generation method by decomposing the problem into a master
problem and a subproblem. The subproblem takes the form of an Elementary
Shortest Path Problem with Resource Constraints and Minimized-Maximum Risk
(ESPPRCMMR). To solve the ESPPRCMMR efficiently, we develop a new linear-time
risk calibration algorithm and establish families of resource extension
functions in compliance with the risk-related resources. Additionally, we
extend the applicability of our setting and solution approach to address a
generic bi-objective DARP, namely Equitable DARP (EDARP), which aims to
distribute ride times among passengers equitably. To demonstrate the
effectiveness of our approach, we adopt a real-world paratransit trip dataset
to generate the RDARP instances ranging from 17 to 55 heterogeneous passengers
with 3 to 13 vehicles during the morning and afternoon periods. Computational
results show that our BCP algorithm can optimally solve 50 out of 54 real-
world instances (up to 55 passengers) within a time limit of one hour. dial-a-
ride problem; branch-cut-and-price; infectious disease; demand-responsive
transit; equity

##  1 Introduction

The _dial-a-ride problem_ (DARP) generalizes the _pick-up and delivery
problems with time windows_ (PDPTWs) with the maximum ride time constraints
(cordeau2007dial), and has seen extensive real-world applications such as
dial-a-ride services, on-demand mobility services, and paratransit services.
For the DARP, each passenger must be picked up and dropped off within a
specific time interval by the same vehicle. The focus of the DARP is primarily
serving demand at a low cost (ho2018survey) by transporting more passengers at
a time with limited resources.

Our study is motivated by the emerging real-world challenges posed by the
COVID-19 pandemic on demand-responsive transit (DRT) operations. Specifically,
the paratransit services, which provide essential mobility services to
individuals eligible under the Americans with Disabilities Act, serve a
vulnerable population such as people with disability, people with pre-existing
health conditions, and senior citizens (nadtc2021). The conventional DARP
setting, which involves transporting more passengers at a time, facilitates
the spread of infectious diseases and is unsuitable during pandemics,
especially when considering the compound exposed risk towards the vulnerable
population who require access to healthcare facilities
(chen2020transportation, nie2022impact). Consequently, there is an urgent need
for DRT service providers such as paratransit, ridesharing, and flexible-route
transit to develop risk-aware operation strategies that provide a safe travel
environment for the vulnerable population in a cost-effective manner.

In response to the emerging disease transmission concerns during the COVID-19
pandemic, we propose the risk-aware routing and scheduling framework for DRT
services. Specifically, we introduce the Risk-aware DARP (RDARP), aiming to
minimize both the total travel cost and the maximum individual exposed risk
among onboard passengers (the reason for considering the min-max exposed risk
will be explained in Remark 3.1). This framework has significant implications
for mitigating the risk of ongoing waves of COVID-19 outbreak, seasonal
communicable diseases (e.g., flu and respiratory infections), and future
unknown diseases that may arise. Additionally, our model can be broadly
extended to other practical applications that require minimal interactions
among people or goods, such as food delivery service for fresh produce and
frozen/hot food, the DRT service with co-rider preferences, patient transport
service, and DARP with an equitable level of service (with the min-max
objective of excessive ride time or detour rate).

This study proposes the mathematical model of the RDARP and a tailored
algorithm to optimally solve it. The RDARP includes provisions for exposed
risk among passengers with different levels of risk, subject to the
propagation of disease risks along the route and the maximum cumulative
exposed risk constraints, in addition to the conventional constraints in the
DARP.

The main contributions of our study are summarized below:

  * •

We propose the RDARP with a particular focus on the en-route exposed risk in
addition to the total travel cost, which is formulated as a bi-objective
optimization problem. Sets of risk-aware routing strategies on the Pareto
front are provided to aid decision-making for DRT service.

  * •

We solve the RDARP by means of the column generation method under a Branch-
Cut-and-Price (BCP) framework. Our BCP model efficiently handles the non-
linearity of the risk-related constraints using tailored valid inequalities
and a labeling algorithm.

  * •

We design a tailored labeling algorithm to provide an exact solution by
leveraging the properties of the dynamic time window in the typical DARP and
disease transmission among the onboard passengers. We consider the exposed
risk among passengers as a function of onboard duration consisting of travel
time, service time, and waiting time, and we introduce a family of risk-
related resources with associated resource extension functions and dominance
rules to accurately capture the minimum possible exposed risk.

  * •

We propose a model extension of our RDARP called Equitable DARP, which aims to
distribute ride times equitably among the passengers with minor modifications
to the RDARP.

  * •

We conduct extensive numerical experiments using the real-world paratransit
trip log data in the central Alabama area. In particular, the trade-off
between travel cost and maximum individual exposed risk is broadly discussed
upon the results of sensitivity analyses and the Pareto fronts.

The remainder of this paper is organized as follows. We review the related
studies in Section 2. Section 3 formally introduces the arc-based formulations
of the RDARP and the equivalent trip-based formulation, which is decomposed
based on the Dantzig-Wolfe (DW) decomposition method. In Section 4, we propose
the labeling algorithm and a risk calibration algorithm to solve the RDARP.
Section LABEL:sec:branch_and_price discusses the branch-and-price algorithm.
The model extension is discussed in Section LABEL:sec:edarp. Section
LABEL:sec:results introduces the data instances and presents the computational
results and discussion. Finally, Section LABEL:sec:conclusion concludes our
study and proposes several future directions.

##  2 Literature Review

Surveys on the DARP are provided by cordeau2007dial and most recently,
ho2018survey. In general, the DARP can be divided into static and dynamic
modes based on whether the demand is pre-assigned (cordeau2007dial). The
static DARP is associated with pre-scheduled demand, e.g., dial-a-ride service
and paratransit service (fu1999improving, gupta2010improving). On the other
hand, the dynamic DARP aims at providing high-quality service to immediate or
short-term demand. Practical applications can be seen in the shared taxi
service (hosni2014shared), on-demand mobility system (sayarshad2018scalable),
and accessible taxi services (dikas2018scheduled). The objectives of DARP
mainly include minimizing the operating costs and/or user inconvenience (for a
comprehensive survey of the objectives in DARP, see paquette2009quality,
molenbruch2017typology), which can be handled by weighting or
$\epsilon$-constraint methods (demir2014bi). In addition, the min-max
objectives have been investigated in the Vehicle Routing Problem with load
balancing or workload equity (matl2018workload), i.e., minimizing the maximum
service completion time (perrier2008vehicle), maximum routing cost
(carlsson2009solving).

The exact algorithms for the DARP mainly focus on the Branch-and-Cut (BC) and
Branch-and-Price (BP) algorithms that are associated with the arc-based and
trip-based formulations, respectively. In the arc-based formulation, each
location is associated with the arrival time and load to ensure the
feasibility regarding the time window, capacity, and ride time constraints. To
tighten the linear programming (LP) relaxation of arc-based formulation,
families of valid inequalities with separation heuristics were developed
(cordeau2006branch, ropke2007models, ropke2009branch), including subtour
elimination constraints, rounded capacity constraints, and infeasible path
elimination constraints (IPEC). Later, BC algorithms were applied to more
practical cases in terms of different user types (parragh2011introducing),
multi-depot setting (braekers2014exact), lunch-break for drivers
(liu2015branch), and driver consistency (paquette2013combining). On the other
hand, the trip-based formulation relies on the BP algorithms. The BP algorithm
for PDPTW was first proposed by dumas1991pickup and extended by
savelsbergh1998drive. ropke2009branch incorporated a set of valid inequalities
in the BC algorithm into the BP procedure, yielding a BCP algorithm. In
addition, gschwind2015effective proposed an exact BCP algorithm that handles
the maximum ride time constraints. The feasibility check of the ride time
constraint is conducted by inspecting the arrival time and the latest delivery
time for each onboard passenger. luo2019two applied the BCP algorithm to solve
the DARP under sets of practical constraints using IPEC and Benders cuts.

Finally, it is worth noting that our concept of ”risk while traveling” differs
from existing studies. The term can be found in the cash transportation
problem (allahyari2021novel), where en-route risk relies only on simple
propagation of the start-of-service time without contagion among onboard items
and varying waiting times. On the contrary, our study considers potential
waiting times and realistic risk-related constraints that are applicable to
real-world operations. Additionally, accurate calculation of the minimal
exposed risk depends on the exact start-of-service time, which has not been
studied previously. Furthermore, although min-max DARP has been investigated
in the equity issues (matl2018workload, lehuede2014multi), no existing exact
BCP algorithms are available to find the optimal solution nor the Pareto
front.

##  3 Problem Setting and Formulations

We develop a risk-aware routing and scheduling framework, motivated by the
emerging challenges posed by the COVID-19 pandemic on DRT operations. This
section starts with the problem definition and the unique features of modeling
the exposed risk among onboard passengers during travel. We then propose an
arc-based formulation for the problem, which aids in understanding its
structure and serves as a stepping stone for the trip-based formulation
detailed next. Additionally, we underscore in Remark 3.1 the significance of
the min-max exposed risk objective in the epidemic process over a physical
contact network.

###  3.1 Problem Definition

Our RDARP is extended based on the static multi-vehicle DARP with additional
consideration of exposed risk among onboard passengers. In our RDARP,
passenger trips are made via reservations with known pick-up and drop-off
locations, as well as the desired pick-up and/or drop-off time windows
depending on outbound or inbound trips. Let $n$ be the number of request
locations. The RDARP is defined on a complete directed graph denoted by
$\mathcal{G}=(\mathcal{N},\mathcal{A)}$. The node set $\mathcal{N}$ includes
four subsets: $\\{0\\}$ and $\\{2n+1\\}$ to represent the origin and final
depots which may coincide; $\mathcal{P}=\\{1,2,\ldots,n\\}$ and
$\mathcal{D}=\\{n+1,n+2,\ldots,2n\\}$ to denote the requests’ pick-up and
drop-off nodes, respectively. Hence, each request is associated with the pick-
up node $i$ and drop-off node $n+i$. Each request may include one or multiple
passengers, where $w_{i}$ indicates the number of passengers associated with
location $i$, and therefore, $w_{i}=-w_{n+i},~{}\forall i\in\mathcal{P}$. For
each node $i\in\mathcal{P}\cup\mathcal{D}$, the time window is represented by
$\left[a_{i},b_{i}\right]$, denoting the earliest and latest start-of-service
time at location $i$. A non-negative service duration $s_{i}$ is considered
for request locations $i\in\mathcal{P}\bigcup\mathcal{D}$, and
$s_{0}=s_{2n+1}=0$ at the depots. Let $\mathcal{K}:=\\{1,2,\ldots,K\\}$ denote
the set of identical vehicles with the capacity of $W_{\max}$. Each arc
$(i,j)\in\mathcal{A}$ is associated with travel time $t_{ij}$. The maximum
allowed ride time for passenger $i$ is $L^{i}_{\max}$. In the context of
disease exposed risk, we define a static risk score $r_{i}$ for every pick-up
node $i\in\mathcal{P}$, which is a parameter estimated by the passenger’s
demographics and trip purpose. For passenger $i\in\mathcal{P}$, analogous to
the passenger count ($w_{i}=-w_{n+i}$), we denote $r_{i}=-r_{n+i}$ associated
with the pick-up and drop-off nodes, which helps to track the onboard
passengers’ risk scores.

For a specific vehicle $k\in\mathcal{K}$, let the binary variable
$x_{ij}^{k}=1$ if vehicle $k$ traverses arc $(i,j)$. For each node
$i\in\mathcal{N}$, let $A_{i}^{k}$ denote the start-of-service time and
$L_{i}^{k}$ represents the ride time for request $i$. Note that, at node
$i\in\mathcal{N}$, we allow the early arrival, where one will wait until the
start of the time window $a_{i}$. Besides, one can wait and delay the start of
the service at the latest to $b_{i}$ to align with the feasibility of maximum
ride times. Let $W_{i}^{k}$ and $R_{i}^{k}$ denote the number of onboard
passengers and total risk scores of onboard passengers at node
$i\in\mathcal{N}$, respectively. At each node $i\in\mathcal{N}$, $W_{i}^{k}$
and $R_{i}^{k}$ are propagated in a similar fashion by addition and
subtraction of the onboard passengers’ demand $w_{i}$ and risk score $r_{i}$
based on the pick-up and drop-off sequences.

To formulate the risk components in the RDARP, we define two additional
variables $Q_{i}^{k}$ and $H_{i}^{k}$. Specifically, $Q_{i}^{k}$ is the
cumulative exposed risk along the vehicle $k$’s route after visiting node $i$.
Note that pathogen can remain in the vehicle environment for hours, the
cumulative risk at destination depot $Q_{2n+1}^{k}$ can then be used to
restrain excessive long trips (e.g., not exceeding $Q_{\max}$) and help guide
the regular disinfection of vehicles. For each traversed arc $(i,j)$, the
exposed risk is associated with onboard passengers’ risk score $R_{i}^{k}$ and
onboard duration on arc $(i,j)$, including travel time, service time, and
potential waiting time or delay. Moreover, $H_{i}^{k}$ is the individual
exposed risk for request $i$’s trip from $i$ to $n+i$ on vehicle $k$, which is
equivalent to the increment between $Q_{i}^{k}$ and $Q_{n+i}^{k}$ after
subtracting the self-transmission regarding the risk score $r_{i}$ and onboard
duration from $i$ to $n+i$. Finally, let $\overline{H}$ be a nonnegative
variable representing an upper bound of the exposed risk among all passengers.

###  3.2 Arc-based Formulation

We first present the RDARP as a bi-objective mixed-integer non-linear
programming problem, established based on the three-index arc-based
formulation (cordeau2006branch).

$\displaystyle\min_{x}~{}$
$\displaystyle{\Biggl{\\{}\sum_{k\in\mathcal{K}}\sum_{i\in\mathcal{N}}\sum_{j\in\mathcal{N}}t_{ij}x_{ij}^{k},\
\overline{H}\Biggl{\\}}}$ (1a) s.t.
$\displaystyle\sum_{k\in\mathcal{K}}H_{i}^{k}\leq\overline{H}$
$\displaystyle\forall i\in\mathcal{P}~{}{}$ (1b)
$\displaystyle\sum_{k\in\mathcal{K}}\sum_{j\in\mathcal{N}}x_{ij}^{k}=1$
$\displaystyle\forall i\in\mathcal{P}~{}$ (1c)
$\displaystyle\sum_{j\in\mathcal{N}\setminus\\{2n+1\\}}x_{ji}^{k}-\sum_{j\in\mathcal{N}\setminus\\{2n+1\\}}x_{j,n+i}^{k}=0$
$\displaystyle\forall i\in\mathcal{P},k\in\mathcal{K}~{}$ (1d)
$\displaystyle\sum_{j\in\mathcal{N}\setminus\\{2n+1\\}}x_{ji}^{k}-\sum_{j\in\mathcal{N}\setminus\\{0\\}}x_{ij}^{k}=0$
$\displaystyle\forall i\in\mathcal{P}\cup\mathcal{D},k\in\mathcal{K}~{}$ (1e)
$\displaystyle\sum_{j\in\mathcal{N}}x_{0j}^{k}=1$ $\displaystyle\forall
k\in\mathcal{K}~{}$ (1f) $\displaystyle\sum_{j\in\mathcal{N}}x_{j,2n+1}^{k}=1$
$\displaystyle\forall k\in\mathcal{K}~{}$ (1g) $\displaystyle A_{j}^{k}\geq
A_{i}^{k}+s_{i}+t_{ij}+M_{1}\left(x_{ij}^{k}-1\right)$ $\displaystyle\forall
i\in\mathcal{N},j\in\mathcal{N},k\in\mathcal{K}~{}$ (1h) $\displaystyle
a_{i}\leq A_{i}^{k}\leq b_{i}$ $\displaystyle\forall
i\in\mathcal{N},k\in\mathcal{K}~{}$ (1i) $\displaystyle
L_{i}^{k}=A_{n+i}^{k}-\left(A_{i}^{k}+s_{i}\right)$ $\displaystyle\forall
i\in\mathcal{P},k\in\mathcal{K}~{}$ (1j) $\displaystyle t_{i,n+i}\leq
L_{i}^{k}\leq L^{i}_{\max}$ $\displaystyle\forall
i\in\mathcal{P},k\in\mathcal{K}~{}$ (1k) $\displaystyle W_{j}^{k}\geq
W_{i}^{k}+w_{j}+M_{2}\left(x_{ij}^{k}-1\right)$ $\displaystyle\forall
i\in\mathcal{N},j\in\mathcal{N},k\in\mathcal{K}~{}$ (1l) $\displaystyle 0\leq
W_{i}^{k}\leq W_{\max}$ $\displaystyle\forall
i\in\mathcal{N},k\in\mathcal{K}~{}$ (1m) $\displaystyle R_{j}^{k}\geq
R_{i}^{k}+r_{j}+M_{3}\left(x_{ij}^{k}-1\right)$ $\displaystyle\forall
i\in\mathcal{N},j\in\mathcal{N},k\in\mathcal{K}~{}$ (1n) $\displaystyle
Q_{j}^{k}\geq
Q_{i}^{k}+\left(A_{j}^{k}-A_{i}^{k}\right)R_{i}^{k}+M_{4}\left(x_{ij}^{k}-1\right)$
$\displaystyle\forall i\in\mathcal{N},j\in\mathcal{N},k\in\mathcal{K}~{}$ (1o)
$\displaystyle 0\leq Q_{2n+1}^{k}\leq Q_{\max}$ $\displaystyle\forall
k\in\mathcal{K}~{}$ (1p) $\displaystyle
H_{i}^{k}=Q_{n+i}^{k}-Q_{i}^{k}-\left(A_{n+i}^{k}-A_{i}^{k}\right)r_{i}$
$\displaystyle\forall i\in\mathcal{P},k\in\mathcal{K}~{}$ (1q) $\displaystyle
x_{ij}^{k}\in\\{0,1\\}$ $\displaystyle\forall
i\in\mathcal{N},j\in\mathcal{N},k\in\mathcal{K}~{}$ (1r)

The objective function (1a) aims at minimizing the total travel time and the
maximum exposed risk among all passengers. Constraints (1b) serves as a min-
max linearization, where the individual exposed risk for passenger
$i\in\mathcal{P}$ is bounded by $\overline{H}$. Constraints (1c) and (1d)
ensure that each passenger is visited exactly once and one passenger’s pick-up
and drop-off nodes are served by the same vehicle. Constraints (1e) are the
flow conservation constraints. Constraints (1f) and (1g) enforce that each
vehicle departs and arrives at the depot exactly once. Constraints (1h) ensure
the arrival time consistency by considering the travel time and service time,
which is bounded by Constraints (1i). Constraints (1j) impose the ride time
for each passenger, bounded by Constraints (1k). Constraints (1l) define the
conservation of the onboard number of passengers at each intermediate node,
and Constraints (1m) limit the total number of onboard passengers by the
vehicle capacity. The remaining constraints are associated with contagion risk
calculation and propagation, and readers may refer to Appendix
LABEL:sec:risk_measure_calculation for an illustrative example. Specifically,
Constraints (1n) define the cumulative per-unit-time risk score at the
intermediate node, which is further used in Constraints (1o) to calculate the
propagation of cumulative risk by multiplying the lapse of in-vehicle time
between two consecutive nodes. To avoid excessive viruses suspended in the
vehicle, we enforce a maximum cumulative risk limit $Q_{\max}$ by Constraints
(1p). With cumulative risk in Constraints (1o), the individual risk exposure
can therefore be calculated by subtracting the amount of self-exposure from
the total risk as in Constraints (1q).

Next, we remark on the rationale behind the inclusion of minimizing the
_maximum_ individual exposed risk as one of the objectives and the challenges
in solving the arc-based formulation above.

###### Remark 3.1

(The role of minimized maximum individual risk) To understand the use of
$\overline{H}:=\max\limits_{i\in\mathcal{P}}(\sum\limits_{k\in\mathcal{K}}H_{i}^{k})$
in the objective function, one may consider the disease-spreading process in
the dial-a-ride system as an epidemic process over physical contacts among all
passengers. In such a contact network, all passengers are the nodes, and the
contact duration between a node pair serves as the weight on edges. Note that
while some of the passengers $i\in\mathcal{P}$ may not be physically adjacent
at the same time, they may still be connected due to the existence of the
vehicles $k\in\mathcal{K}$, e.g., the interaction of drivers and cleaning
technicians at the depot. The spread of an infectious disease can therefore be
modeled as a percolation process (meyers2007contact) and the level of risk is
determined by a network invasion threshold similar to the idea of the basic
reproduction number $\overline{r}_{0}$, which is used to estimate the speed of
disease spreading (for a comprehensive survey, see guerra2017basic). In
particular, the disease transmission between two nodes (passengers $i$ and
$j$) takes place with a transmission rate $\zeta_{ij}$ that is positively
correlated to the contact duration between two passengers and an infected
individual can fully recover at rate $\overline{r}_{c}$. Specifically,
qian2021scaling showed that the disease can be eliminated if the degree
distribution of the network satisfies:

$\max_{i\in\mathcal{N}}\left(\sum_{j\in\mathcal{N}}\zeta_{ij}\right)<\overline{r}_{c}$
(2)

where the left-hand side denotes the maximum individual cumulative exposed
risk among all passengers (see Proposition 1 in qian2021scaling), which
coincides $\overline{H}$ in our case. Minimizing $\overline{H}$ can be
considered as minimizing the upper bound of the disease threshold, where the
optimization of the latter is challenging (overton1988minimizing). Instead,
the use of $\overline{H}$ provides a linear function that can be embedded
nicely into the RDARP framework and is shown to reach an exponential reduction
in the number of infections (qian2021connecting).

Finally, the arc-based formulation is known to be inefficient due to the weak
LP relaxation (barnhart1998branch) in the BC process. Our RDARP problem
further exacerbates the weakness due to the introduction of bilinear
constraints (Constraints (1o)) as the product of two variables, arrival time
and cumulative risk score. This will introduce the spatial branch and bound
algorithm to handle the bilinear constraints inside the BC framework, which
requires exponential computation time for single-objective problems, letting
alone the further challenge arising from the bi-objective function and the
min-max term in our setting. We hence turn to the BCP framework and show that
the problem can be efficiently solved with a new labeling algorithm by
resolving the bilinear constraints within the pricing subproblems.

###  3.3 Trip-based Formulation

In this section, we reformulate the RDARP into a set-partitioning problem
based on the DW Decomposition (dantzig1960decomposition), which takes the form
of a trip-based formulation.

Let $\Omega$ be the set of elementary trips satisfying the RDARP constraints
(see Eqs. (1d)-(1q)). A feasible route
$r:=\left(0,\cdots,2n+1\right)\in\Omega$ must satisfy the following conditions
for node $i\in r$: (1) if $i\in\mathcal{P}$ then $n+i\in r$; if
$i\in\mathcal{D}$ then $i-n\in r$, (2) $A_{i}\in[a_{i},b_{i}]$, (3) $W_{i}\leq
W_{\max}$, (4) $0<L_{i}\leq L_{\max}^{i},\forall i\in\mathcal{P}$, and (5)
$Q_{i}\leq Q_{\max}$. The travel cost of route $r\in\Omega$ is denoted by
$c_{r}$; parameter $H_{ir}$ is the individual exposed risk of request
$i\in\mathcal{P}$ on route $r$. If request $i$ is not part of route $r$, we
set $H_{ir}=0$. The maximum exposed risk on route $r$ is bounded by variable
$\overline{H}$. Parameter $\alpha_{ir}$ is a binary parameter indicating
whether request $i$ is visited in route $r$. Binary variable $\lambda_{r}$
equals 1 if route $r$ is selected to be part of the solution and $0$
otherwise. The trip-based formulation takes the following form.

$\displaystyle\min$
$\displaystyle\Biggl{\\{}\sum_{r\in\Omega}c_{r}\lambda_{r},\
{\overline{H}}\Biggl{\\}}~{}$ (3a) s.t.
$\displaystyle{\sum_{r\in\Omega}H_{ir}\lambda_{r}\leq\overline{H}}$
$\displaystyle\forall i\in\mathcal{P}$ (3b)
$\displaystyle\sum_{r\in\Omega}\alpha_{ir}\lambda_{r}=1$ $\displaystyle\forall
i\in\mathcal{P}$ (3c) $\displaystyle\sum_{r\in\Omega}\lambda_{r}\leq K$ (3d)
$\displaystyle\lambda_{r}\in\\{0,1\\}$ $\displaystyle\forall r\in\Omega~{}$
(3e)

where the exposed risk of each request is bounded by $\overline{H}$ in
Constraints (3b). Constraints (3c) guarantee that each request is served
exactly once while Constraint (3d) poses a cap on the number of vehicles used,
which cannot exceed the fleet size, $K$. The trip-based formulation decomposes
the problem into trips $r\in\Omega$, each operated by one vehicle, with route-
dependent risk information $H_{ir}$ for passenger $i$ when served on trip $r$.
Note that in the arc-based formulation, the arrival time to each passenger is
regulated using a decision variable, and therefore, the model can arrange the
arrival times such that the maximum individual exposed risk is minimized. The
trip-based formulation does not possess a mechanism to arrange the arrival
times at passenger locations along the route to guarantee that the maximum
individual exposed risk is minimized. Therefore, the time schedule along a
trip $\Gamma_{r}$ shall be optimized to meet the Minimized-Maximum Risk (MMR)
condition, defined as follows.

######  Definition 3.2 (MMR condition)

A route $r$, together with a given arrival time schedule $\Gamma_{r}$, are
said to satisfy the MMR condition if the maximum total exposed risk of any
individual on the route, $\max\limits_{i\in\mathcal{P}\cap r}H_{ir}$, can no
further be reduced by changing $\Gamma_{r}$.

###### Proposition 3.3

The arc-based and trip-based formulations are equivalent if every route
$r\in\Omega$ satisfies the MMR condition.

###### Proof 3.4

Proof. See Appendix LABEL:proof:equivalence.

Based on the above statements, the key to the trip-based formulation is the
accurate calibration of MMR, which can be further decomposed into the
identification of a feasible $\Gamma_{r}$ that minimizes the maximum
individual risk along the route. Note that a passenger’s exposed risk is
accumulated from the contact times with other passengers, while the contact
time can be broken down into arc-level travel time, node-level service time,
and additional waiting time due to the early arrival at a pick-up node. Given
a route $r$, we observe that both travel and service times are constant.
Therefore, the only component of a route that can potentially be altered to
minimize the maximum risk is the waiting times before the start of service at
a location. This can be made possible by delaying the service time at the
predecessor nodes but at the risk of violating the latest arrival time and/or
maximum trip length constraints for some of the passengers. Therefore, the
calibration of MMR requires strategically adjusting the start-of-service time
at the pick-up nodes while ensuring the trip remains feasible for all
passengers.

Nevertheless, ensuring the MMR condition for all serviced passengers given a
route $r$ is a nontrivial task for the following two reasons. First, there is
an inherent trade-off between earlier and later start-of-service times.
Delaying a pick-up can reduce onboard passengers’ waiting time and decrease
exposed risk. However, it may also lead to delayed arrival at subsequent nodes
hence the violation of the latest arrival time constraints. Second, due to the
interlacing of pick-up and drop-off sequences, the optimal time delay becomes
a dynamic programming problem that requires new algorithm development for the
pricing subproblem. Existing DARP studies based on trip-based formulation
(ropke2009branch, gschwind2015effective) only concern the feasibility of
$\Gamma_{r}$ for generated trips, without considering optimization around
$\Gamma_{r}$ which is essential to meet MMR condition at the trip level. In
the following section, we will first introduce the basic BCP components for
RDARP, and we will then return to our specific algorithm for optimizing
schedule $\Gamma_{r}$ in Section 4.1.3.

##  4 Solution Approach

Since the cardinality of $\Omega$ is extremely large, the approach used is
based on a Column Generation (CG) algorithm to solve a restricted linear
relaxation of the trip-based formulation, namely the _restricted linear master
problem_ (RLMP). Specifically, the RLMP is limited to a subset of feasible
trips, which is of moderate size and can be directly solved by the simplex
algorithm. The CG method is an iterative procedure embedded in the branch-and-
bound framework. At each branching node, the RLMP and the pricing subproblem
(SP) are solved iteratively. The RLMP is first solved with a set of initial
columns, and the resulting dual variables are used to generate feasible
columns to be included in the basis in the SP. The RLMP is then resolved with
the updated set of columns. This process is repeated until no valid columns
with negative reduced costs are identified. If an integer solution is obtained
at the root node, it is considered optimal. Otherwise, the branching procedure
is initiated, where we compare and update the best integer solution with the
incumbent solutions identified at branching nodes. During the branching
procedure, invalid columns are eliminated and new columns may be generated in
compliance with the branching rules. Finally, the best IP solution is reported
to be the optimal solution of the original problem (dayarian2015column).

#####  Restricted Linear Master Problem (RLMP).

The RLMP is constructed by relaxing the integrality constraints (3e) in the
trip-based formulation. Instead of handling an exponential number of columns,
the RLMP is restricted to a subset of feasible trips
$\Omega^{\prime}\subset\Omega$. In addition, we adopt the
$\epsilon$-constraint method to handle the bi-objective nature of the
objective function (3a). The $\epsilon$-constraint method proceeds by
optimizing one objective function while considering the other objective
function as a constraint with a predefined bound. Let $\epsilon^{\text{risk}}$
and $\epsilon^{\text{cost}}$ be the predetermined upper bounds of the maximum
exposed risk and total travel cost. We will consider two RLMPs,
$\BFP^{\text{cost}}(\epsilon^{\text{risk}})$ and
$\BFP^{\text{risk}}(\epsilon^{\text{cost}})$, aiming to minimize the total
travel cost and maximum exposed risk while imposing the other objectives with
the upper bounds of $\epsilon^{\text{risk}}$ and $\epsilon^{\text{cost}}$,
respectively. The RLMPs take the form of

$\displaystyle\BFP^{\text{cost}}(\epsilon^{\text{risk}}):$
$\displaystyle\min~{}\sum_{r\in\Omega^{\prime}}c_{r}\lambda_{r}$
$\displaystyle\BFP^{\text{risk}}(\epsilon^{\text{cost}}):\min~{}\overline{H}$
(4a) s.t. $\displaystyle~{}\sum_{r\in\Omega^{\prime}}\alpha_{ir}\lambda_{r}=1$
$\displaystyle\forall i\in\mathcal{P}$
$\displaystyle\text{s.t.}~{}\sum_{r\in\Omega^{\prime}}\alpha_{ir}\lambda_{r}=1$
$\displaystyle\forall i\in\mathcal{P}$
$\displaystyle{(\pi_{i}\in\mathbb{R})}~{}$ (4b)
$\displaystyle\sum\limits_{r\in\Omega^{\prime}}\lambda_{r}\leq K$
$\displaystyle\sum\limits_{r\in\Omega^{\prime}}\lambda_{r}\leq K$
$\displaystyle{(\mu\leq 0)}~{}$ (4c)
$\displaystyle\sum_{r\in\Omega^{\prime}}H_{ir}\lambda_{r}\leq\epsilon^{\text{risk}}$
$\displaystyle\forall i\in\mathcal{P}$
$\displaystyle\sum_{r\in\Omega^{\prime}}H_{ir}\lambda_{r}\leq\overline{H}$
$\displaystyle\forall i\in\mathcal{P}$ $\displaystyle{(\rho_{i}\leq 0)}~{}$
(4d)
$\displaystyle\sum_{r\in\Omega^{\prime}}c_{r}\lambda_{r}\leq\epsilon^{\text{cost}}$
$\displaystyle{(\xi\leq 0)}~{}$ (4e) $\displaystyle\lambda_{r}\geq 0$
$\displaystyle\forall r\in\Omega^{\prime}$ $\displaystyle\lambda_{r}\geq 0$
$\displaystyle\forall r\in\Omega^{\prime}$ (4f)

where $\pi:=\\{\pi_{i}\in\mathbb{R}:\forall i\in\mathcal{P}\\}$, $\mu\leq 0$,
$\rho:=\\{\rho_{i}\leq 0:\forall i\in\mathcal{P}\\}$, and $\xi\leq 0$ are the
dual variables associated with Constraints (4b)-(4e), respectively. In
$\BFP^{\text{cost}}(\epsilon^{\text{risk}})$, one can convert the set
partitioning constraints (4b) to set covering constraints. However, in
$\BFP^{\text{risk}}(\epsilon^{\text{cost}})$, using the set covering
constraint may result in multiple visits in some extreme cases. For instance,
given the riskiest passenger with strictly positive risk $\overline{H}>0$, the
passengers with no co-riders ($H_{ir}=0$) can be visited multiple times while
satisfying Constraints (4c)-(4e). Therefore, we keep the set partitioning
form.

We use the $\epsilon$-constraint method to obtain the exact Pareto optimal
solutions. We adopt an exact method (berube2009exact) by iteratively solving
$\BFP^{\text{cost}}(\epsilon^{\text{risk}})$ and
$\BFP^{\text{risk}}(\epsilon^{\text{cost}})$. At each iteration, the objective
function value of the former is considered as the upper bound of
$\epsilon^{\text{cost}}$ for the latter, and the latter is solved to avoid the
dominated solutions (berube2009exact). For brevity, we present the detailed
algorithm with an illustration in Appendix
LABEL:sec:Pareto_front_approximation.

##### Pricing Subproblem (SP).

The columns in the subset $\Omega^{\prime}$ are generated iteratively by
solving the SP. The SP aims to identify the columns with negative reduced cost
or prove the optimality of RLMP if none is found. The reduced cost functions
$\tilde{c}_{r}^{\text{cost}}$ and $\tilde{c}_{r}^{\text{risk}}$, associated
with route $r\in\Omega$, are shown in Eqs. (5a) and (5b). Further, the reduced
cost functions can be rewritten back to arc-level to incorporate the RDARP
Constraints (1d)-(1q). Given an independent trip $r$, we omit the dimension
$k$ in $x_{ij}^{k}$ and Constraints (1d)-(1q). The reduced costs take the form
of

$\displaystyle{\mathbf{P}_{\text{cost}}(\epsilon^{\text{risk}}):\tilde{c}_{r}^{\text{cost}}=c_{r}-\sum_{i\in\mathcal{P}}\alpha_{ir}\pi_{i}-\sum_{i\in\mathcal{P}}H_{ir}\rho_{i}-\mu=\sum_{(i,j)\in
r}\tilde{c}_{ij}^{\text{cost}}x_{ij}-\sum_{i\in\mathcal{P}}H_{ir}\rho_{i}-\mu}~{}$
(5a)
$\displaystyle{\mathbf{P}_{\text{risk}}(\epsilon^{\text{cost}}):\tilde{c}_{r}^{\text{risk}}=-\sum_{i\in\mathcal{P}}\alpha_{ir}\pi_{i}-\sum_{i\in\mathcal{P}}H_{ir}\rho_{i}-\mu-\xi
c_{r}=\sum_{(i,j)\in
r}\tilde{c}_{ij}^{\text{risk}}x_{ij}-\sum_{i\in\mathcal{P}}H_{ir}\rho_{i}-\mu}~{}$
(5b)

where $\tilde{c}_{ij}^{\text{cost}}=t_{ij}-\pi_{i},\forall i\in\mathcal{P}$
and $\tilde{c}_{ij}^{\text{cost}}=t_{ij}$ otherwise;
$\tilde{c}_{ij}^{\text{risk}}=-\xi t_{ij}-\pi,\forall i\in\mathcal{P}$ and
$\tilde{c}_{ij}^{\text{risk}}=-\xi t_{ij}$ otherwise. Based on the expressions
above, the pricing subproblem can be formulated as follows.

$\displaystyle\min$ $\displaystyle{\sum\limits_{(i,j)\in
r}\tilde{c}_{ij}x_{ij}-\sum_{i\in\mathcal{P}}H_{ir}\rho_{i}-\mu}~{}$ (6a) s.t.
$\displaystyle\eqref{eq:served_by_one}-\eqref{eq:direct_infection}$ (6b)

where the arc-level reduced cost $\tilde{c}_{ij}=\tilde{c}_{ij}^{\text{cost}}$
in $\mathbf{P}_{\text{cost}}(\epsilon^{\text{risk}})$ and
$\tilde{c}_{ij}=\tilde{c}_{ij}^{\text{risk}}$ in
$\mathbf{P}_{\text{risk}}(\epsilon^{\text{cost}})$.

The SP is solved in the same graph $\mathcal{G}=(\mathcal{N},\mathcal{A})$ by
replacing the arc cost $t_{ij}$ with $\tilde{c}_{ij}$, for arc
$(i,j)\in\mathcal{A}$. The SP with Constraints (1d)-(1m) is a typical
Elementary Shortest Path Problem with Resource Constraints (ESPPRC), which is
strongly NP-hard (dror1994note). To guarantee the MMR condition (see
Definition 3.2 and Proposition 3.3), we impose additional constraints on our
SP, resulting in the ESPPRC with the MMR condition (ESPPRCMMR). In contrast to
previous studies that consider a linear cost function of travel time $t_{ij}$
(fukasawa2016branch, luo2017branch), our study considers the exposed risk,
which depends on the sequence of pick up and drop off of passengers and the
varying onboard times including the potential waiting time due to early
arrivals. To address this problem, we propose a dynamic programming-based
algorithm, specifically the labeling algorithm, which is customized to satisfy
the MMR condition in our RDARP.

###  4.1 Labeling Algorithm

We solve the SP (formulated as ESPPRCMMR) by means of the labeling algorithm
(dumas1991pickup, feillet2004exact). In the labeling algorithm, each label
represents a partial trip from origin $0$ to node
$i\in\mathcal{N}\setminus\\{0\\}$ that is extended along the arcs in graph
$\mathcal{G}=(\mathcal{N},\mathcal{A})$. The labeling algorithm aims to
generate a complete trip that satisfies the RDARP constraints and contributes
to the LP relaxation of the RLMP. Specifically, for each potential feasible
label extension along arc $(i,j)\in\mathcal{A}$, the stored information
(namely resources) in the label are updated following resource extension
functions (REFs) (desaulniers1998unified) (see Sections 4.1.1-4.1.3). Based on
the updated resources, the feasibility of the label is checked regarding the
RDARP constraints. However, the enumeration of label extensions can still
result in an exponential size of explorations. In this regard, label dominance
rules are required to terminate unpromising labels during the extension
procedure (see Section LABEL:subsec:dominance_darp).

In our labeling algorithm, we define a label as follows.

###### Definition 4.1

Let $l$ denote a label, representing a partial trip that starts from origin
depot $0$ and resides at the current node $\eta_{l}$. The label can be encoded
as

$l:=\left(\underbrace{\eta_{l},\tilde{c}_{l},A_{l},B_{l},W_{l},\mathcal{V}_{l},\mathcal{O}_{l}}_{\text{PDPTW-
related}},\underbrace{\mathcal{O}_{l}^{a},R_{l},Q_{l},\\{h_{l}^{o},d_{l}^{o}:o\in\mathcal{O}_{l}\cup\mathcal{V}_{l}\\}}_{\text{Risk-
related}},\underbrace{\\{B_{l}^{o},D_{l}^{o}(A_{l}),D_{l}^{o}(B_{l}^{o}):o\in\mathcal{O}_{l}\\}}_{\text{Ride
time-related}}\right)$

  * •

$\eta_{l}$: the node where label $l$ resides

  * •

$\tilde{c}_{l}$: reduced cost after serving node $\eta_{l}$

  * •

$A_{l}$: the earliest possible start-of-service time at node $\eta_{l}$ given
a partial path $l$

  * •

$\Gamma_{l}$: a sequence of time schedules that record the feasible start-of-
service times, denoted as
$\Gamma_{l}^{0}<\Gamma_{l}^{i}<\Gamma_{l}^{j}<\cdots$ for route
$(0,i,j,\cdots)$

  * •

$B_{l}$: the latest possible start-of-service time at node $\eta_{l}$ given a
partial path $l$

  * •

$W_{l}$: number of requests on board at node $\eta_{l}$

  * •

$\mathcal{V}_{l}$: the set of serviced requests along the partial trip until
node $\eta_{l}$

  * •

$\mathcal{O}_{l}$: the set of open requests, where request
$o\in\mathcal{O}_{l}$ has been picked up but not dropped off

  * •

$\mathcal{O}_{l}^{a}$: the set of associated requests at $\eta_{l}$, which
includes the requests that have been dropped off but _were exposed to the open
requests $o\in\mathcal{O}_{l}$_

  * •

$R_{l}$: sum of risk scores for all onboard requests $o\in\mathcal{O}_{l}$ at
node $\eta_{l}$

  * •

$Q_{l}$: cumulative risk at node $\eta_{l}$

  * •

$h_{l}^{o}$: individual exposed risk of request
$o\in\mathcal{O}_{l}\cup\mathcal{V}_{l}$ at node $\eta_{l}$

  * •

$d_{l}^{o}$: available delay buffer for request
$o\in\mathcal{O}_{l}\cup\mathcal{V}_{l}$ when starting from node $\eta_{l}$
(see Section LABEL:sec:risk_min_without_associated)

  * •

$D_{l}^{o}(\tau)$: latest possible drop-off time for request
$o\in\mathcal{O}_{l}$ as a function of start-of-service time $\tau$

  * •

$B_{l}^{o}$: a breaking point representing a feasible start-of-service time at
node $\eta_{l}$ for an open request $o\in\mathcal{O}_{l}$; If $A_{l}\geq
B_{l}^{o}$, it implies that $D_{l}^{o}(A_{l})$ is bounded by the latest time
window for pick-up ($b_{o}$) or drop-off ($b_{o+n}$).

In particular, $\eta_{l}$ is recorded for label $l$ to identify the labels
with the same ending node, and $\tilde{c}_{l}$ is used to prioritize the
labels with lower reduced cost. The resources $A_{l}$ and $W_{l}$ are
considered to keep track of capacity and time window constraints.
$\mathcal{V}_{l}$ and $\mathcal{O}_{l}$ are designed for pairing and
precedence constraints. Set $\mathcal{O}_{l}^{a}$ is introduced to ensure that
the adjustment for the time schedule is in compliance with the MMR condition
for both open requests as well as the associated requests, where an
illustrative example is shown in Figure 1. For example, at node $j+n$,
passenger $k$ is the only onboard passenger hence $\mathcal{O}_{l}=\\{k\\}$.
Due to exposure between $j$ and $k$ prior to the arrival, $j$ is an associated
request to $k$ and $\mathcal{O}_{l}^{a}=\\{j\\}$.

Figure 1: Illustration for the set of open and associated requests on an
interlaced
trip$0$$i$$\mathcal{O}_{l}=\\{i\\}$$\mathcal{O}_{l}^{a}=\varnothing$$i+n$$\mathcal{O}_{l}=\varnothing$$\mathcal{O}_{l}^{a}=\varnothing$$j$$\mathcal{O}_{l}=\\{j\\}$$\mathcal{O}_{l}^{a}=\varnothing$$k$$\mathcal{O}_{l}=\\{j,k\\}$$\mathcal{O}_{l}^{a}=\varnothing$$j+n$$\mathcal{O}_{l}=\\{k\\}$$\mathcal{O}_{l}^{a}=\\{j\\}$$m$$\mathcal{O}_{l}=\\{k,m\\}$$\mathcal{O}_{l}^{a}=\\{j\\}$

In addition, resources $B_{l}$, $D_{l}^{o}(\tau)$, and $B_{l}^{o}$ are
introduced to handle the maximum ride time constraint. The main idea is to
delay the start-of-service time such that potential waiting time can be
reduced and total ride time is more likely to be feasible for the unknown
subsequent nodes.

To ensure feasibility with respect to the risk-related components, we propose
four risk-related resources, $R_{l}$,
$\\{h_{l}^{o},d_{l}^{o}:o\in\mathcal{O}_{l}\cup\mathcal{V}_{l}\\}$, and
$Q_{l}$. The resource $R_{l}$ keeps track of total risk scores of visited
requests $i\in\mathcal{O}_{l}\cup\mathcal{V}_{l}$. The individual exposed risk
$h_{l}^{o}$ is propagated at the arc level along the partial path of label $l$
by considering the risk scores of co-riders and the duration of onboard
travel, which can be calibrated by utilizing the delay buffer $d_{l}^{o}$.
Lastly, $Q_{l}$ is used to verify the feasibility of the maximum cumulative
risk constraints ($Q_{l}\leq Q_{\max}$).

Algorithm 1 Labeling algorithm

$\mathcal{G}=(\mathcal{N},\mathcal{A})$,
{$a_{i},b_{i},n_{i},s_{i},r_{i}:\forall i\in\mathcal{N}\\},\\{t_{ij}:\forall
i,j\in\mathcal{A}\\}$ $\triangleright$ Input
Label $l^{\star}$ with most negative reduced cost $\triangleright$ Output

1:Initiate the label,
$l_{0}\leftarrow\left(0,0,a_{0},b_{0},0,\\{0,b_{0}-a_{0}\\},\varnothing,\varnothing,0,0,\varnothing,\varnothing\right)$

2:Initialize a priority queue:
$\mathcal{Q}\leftarrow(\tilde{c}_{l_{0}},l_{0})$. Initialize a collection of
sets of associated labels with the same ending node,
$\mathcal{L}=\\{\mathcal{L}\left(i\right):i\in\mathcal{N}\\}$. Initialize
$\tilde{c}^{\star}=0$.

3:while $\mathcal{Q}$ is not empty do

4: $l\leftarrow\mathcal{Q}.\text{pop}()$

5: for all $j\in\mathcal{N}$ do

6: Update
$\eta_{l^{\prime}},A_{l^{\prime}},W_{l^{\prime}},\mathcal{V}_{l^{\prime}},\mathcal{O}_{l^{\prime}}$
following the REF rules in Eqs. (7a)-(7e).$\triangleright$ PDPTW-related
Updates

7: if $(\eta_{l},\eta_{l^{\prime}})$ is a PDPTW-feasible extension
then$\triangleright$ DARP-related Updates

8: for all $o\in\mathcal{O}_{l^{\prime}}$ do

9: Update
$B_{l^{\prime}}^{o},D_{l^{\prime}}^{o}(A_{l^{\prime}}),D_{l^{\prime}}^{o}(B_{l^{\prime}}^{o})$
following REFs in Eqs. (8a)-(8d)

10: end for

11: if $(\eta_{l},\eta_{l^{\prime}})$ is a DARP-feasible extension
then$\triangleright$ RDARP-related Updates

12: Update
$\mathcal{O}_{l^{\prime}}^{a},R_{l^{\prime}},Q_{l^{\prime}},\\{h_{l^{\prime}}^{o},d_{l^{\prime}}^{o}:o\in\mathcal{O}_{l}\cup\mathcal{O}_{l}^{a}\\},\tilde{c}_{l^{\prime}}$
following REFs in Eqs.
(LABEL:ref:associated_request_set)-(LABEL:REF:reduced_cost) and Algorithm
LABEL:alg:contact_duration_calibriation

13: if $(\eta_{l},\eta_{l^{\prime}})$ is a RDARP-feasible extension then

14:
$l^{\prime}\leftarrow\left(\eta_{l^{\prime}},\tilde{c}_{l^{\prime}},A_{l^{\prime}},W_{l^{\prime}},\mathcal{V}_{l^{\prime}},\mathcal{O}_{l^{\prime}},\mathcal{O}_{l^{\prime}}^{a},\\{h_{l^{\prime}}^{o},d_{l^{\prime}}^{o}:o\in\mathcal{O}_{l}\cup\mathcal{O}_{l}^{a}\\},R_{l^{\prime}},Q_{l^{\prime}},\\{B_{l^{\prime}}^{o},D_{l^{\prime}}^{o}(A_{l^{\prime}}),D_{l^{\prime}}^{o}(B_{l^{\prime}}^{o}):o\in\mathcal{O}_{l^{\prime}}\\}\right)$

15: Check if $l^{\prime}$ is dominated by the existing labels in
$\mathcal{L}\left(\eta_{l^{\prime}}\right)$ following Eqs.
(LABEL:dom:dominance_rule_darp0)-(LABEL:dom:dominance_rule_darp3)

16: if $l^{\prime}$ is not dominated then

17:
$\mathcal{Q}\leftarrow\mathcal{Q}\cup\\{(\tilde{c}_{l^{\prime}},l^{\prime})\\}$,
$\mathcal{L}\left(\eta_{l^{\prime}}\right)\leftarrow\mathcal{L}\left(\eta_{l^{\prime}}\right)\cup\\{l^{\prime}\\}$

18: end if

19: for all $l_{k}\in\mathcal{L}\left(\eta_{l^{\prime}}\right)$ do

20: if $l^{\prime}$ dominates $l_{k}$ then

21: $\mathcal{Q}\leftarrow\mathcal{Q}\setminus\\{(\tilde{c}_{k},l_{k})\\}$,
$\mathcal{L}\left({\eta_{l^{\prime}}}\right)\leftarrow\mathcal{L}\left(\eta_{l^{\prime}}\right)\setminus\\{l_{k}\\}$

22: end if

23: end for

24: if $\eta_{l^{\prime}}=2n+1$ and
$\tilde{c}_{\eta_{l^{\prime}}}<\tilde{c}^{\star}$ then

25: Update the current best reduced cost and label:
$\tilde{c}^{\star}=\tilde{c}_{\eta_{l^{\prime}}}$, $l^{\star}=l^{\prime}$

26: end if

27: end if

28: end if

29: end if

30: end for

31:end while

Our labeling algorithm proceeds as shown in Algorithm 1. The labeling
algorithm describes a label extension procedure from the origin depot $0$ to
final depot $2n+1$. For the label extension from $l$ to $l^{\prime}$, the
feasibility checks and the updates of related resources consist of three
levels, which are detailed as follows in compliance with line numbers in
Algorithm 1. PDPTW (line 7): The PDPTW-related resources are first updated to
the current node $\eta_{l^{\prime}}$. An extension is said to be PDPTW-
feasible if it satisfies the time window, capacity, as well as pairing and
precedence constraints, corresponding to
$A_{l}+s_{\eta_{l}}+t_{\eta_{l},\eta_{l^{\prime}}}\leq b_{\eta_{l^{\prime}}}$,
$W_{l}+w_{\eta_{l^{\prime}}}\leq W_{\max}$,
$\eta_{l^{\prime}}-n\in\mathcal{O}~{}\text{if}~{}\eta_{l^{\prime}}\in\mathcal{D}$
$\eta_{l^{\prime}}\not\in\mathcal{V}_{l}\text{~{}if~{}}\eta_{l^{\prime}}\in\mathcal{P}$,
and $\mathcal{O}=\varnothing$ if $\eta_{l^{\prime}}=2n+1$, respectively. DARP
(line 11): For the ride time constraint, we verify the existence of a feasible
time schedule by checking if there is a valid drop-off time $A_{l^{\prime}}$
within the feasible range of the time window and ride time constraints, such
that $A_{l^{\prime}}\leq D_{l^{\prime}}^{o}(A_{l^{\prime}})$. RDARP (line 13):
The risk-related resources are updated to obtain the individual risks and
cumulative exposed risk along the route. The extension is feasible with
respect to RDARP if $Q_{l^{\prime}}\leq Q_{\max}$. Finally, the label
dominance rules (lines 16-26) are conducted to eliminate unpromising labels,
thus accelerating the labeling procedure. The labeling algorithm terminates
when all non-dominated labels are extended and reached node $2n+1$. The final
output is a set of RDARP-feasible trips. If no trip with negative reduced cost
$\tilde{c}_{l}<0$ is found, then the RLMP is proven to be optimal.

We note that the proposed labeling algorithm is tailored to handle the risk
calibration process in the SP with three major differences as compared to
existing ones for DARP (ropke2009branch, qu2015branch, gschwind2015effective).
First, we add a family of risk-related resources, $R_{l}$, $Q_{l}$,
$h_{l}^{i}$, and $d_{l}^{i}$, to capture the risk propagation dynamics and
update individual exposed risk $h_{l}^{i}$ with label extensions. Second, we
introduce new resources, dominance rule, and risk calibration algorithm
collectively reduce the search space and guarantee the simultaneous
satisfaction of risk-related constraints and the standard DARP constraints
(ensuring feasibility). Finally, we propose a risk calibration algorithm
(proceed in line 14 in Algorithm 1, see details in Algorithm
LABEL:alg:contact_duration_calibriation) that tracks and updates the delaying
potential for each request, which is used to determine the optimal
$\Gamma_{r}$ during label extension and guarantee the satisfaction of the MMR
condition (ensuring optimality).

##### Resource Extension Functions.

The REFs either initiate the resources or update the resources following a
label extension. By updating the resources in the RDARP, the feasibility of
the partial route (the PDPTW constraints, maximum ride time constraints, and
maximum cumulative risk constraints) is sequentially checked. The REFs for the
RDARP are described in three parts: PDPTW-related resources, DARP-related
resources, and RDARP-related resources in Sections 4.1.1-4.1.3.

####  4.1.1 REF for PDPTW.

While extending labels $l$ to $l^{\prime}$, we update the first set of
resources following the REFs (7a)-(7e), which have been applied to the PDPTW.

$\displaystyle\eta_{l^{\prime}}$ $\displaystyle\leftarrow\eta_{l}$ (7a)
$\displaystyle A_{l^{\prime}}$
$\displaystyle\leftarrow\max\\{A_{l}+s_{\eta_{l}}+t_{\eta_{l},\eta_{l^{\prime}}},a_{\eta_{l^{\prime}}}\\}$
(7b) $\displaystyle W_{l^{\prime}}$ $\displaystyle\leftarrow
W_{l}+w_{\eta_{l^{\prime}}}$ (7c) $\displaystyle\mathcal{V}_{l^{\prime}}$
$\displaystyle\leftarrow\mathcal{V}_{l}\cup\\{\eta_{l^{\prime}}-n\\}\text{~{}if~{}}\eta_{l^{\prime}}\in\mathcal{D}$
(7d) $\displaystyle\mathcal{O}_{l^{\prime}}$
$\displaystyle\leftarrow\begin{cases}\mathcal{O}_{l}\cup\\{\eta_{l^{\prime}}\\}&\text{if}~{}\eta_{l^{\prime}}\in\mathcal{P}\\\
\mathcal{O}_{l}\setminus\\{\eta_{l^{\prime}}-n\\}&\text{if}~{}\eta_{l^{\prime}}\in\mathcal{D}\\\
\end{cases}$ (7e)

Eq. (7a) updates the current node from $\eta_{l}$ to $\eta_{l^{\prime}}$. Eqs.
(7b)-(7c) update the earliest arrival time and number of requests onboard. Eq.
(7d) updates the set of serviced nodes by adding the corresponding request
$\eta_{l^{\prime}}-n$ after dropping off at node
$\eta_{l^{\prime}}\in\mathcal{D}$. Finally, Eqs. (7e) record the set of open
requests after serving node $\eta_{l^{\prime}}$.

####  4.1.2 REF for DARP.

If the arc extension $(\eta_{l},\eta_{l^{\prime}})$ is feasible according to
the PDPTW constraints, we update the set of resources for all open requests
$o\in\mathcal{O}_{l^{\prime}}$ to account for constraints related to maximum
ride time and risk. For each $o\in\mathcal{O}_{l^{\prime}}$, we update the
resources related to ride time as the REFs below.

$\displaystyle B_{l^{\prime}}^{o}$
$\displaystyle\leftarrow\begin{cases}\min\\{B_{l^{\prime}},b_{\eta_{l^{\prime}}+n}-s_{\eta_{l^{\prime}}}-L_{\max}^{\eta_{l^{\prime}}}\\}&\text{if}~{}\eta_{l^{\prime}}=o\\\
\max\\{A_{l^{\prime}},\min\\{B_{l}^{o}+s_{\eta_{l}}+t_{\eta_{l},\eta_{l^{\prime}}},B_{l^{\prime}}\\}\\}&\text{otherwise}\\\
\end{cases}$ (8a) $\displaystyle D_{l^{\prime}}^{o}(A_{l^{\prime}})$
$\displaystyle\leftarrow\begin{cases}\min\\{A_{l^{\prime}}+s_{\eta_{l^{\prime}}}+L_{\max}^{\eta_{l^{\prime}}},b_{\eta_{l^{\prime}}+n}\\}&\text{if}~{}\eta_{l^{\prime}}=o\\\
D_{l}^{o}(A_{l})+\min\\{A_{l^{\prime}}-s_{\eta_{l}}-t_{\eta_{l},\eta_{l^{\prime}}}-A_{l},B_{l}^{o}-A_{l}\\}&\text{otherwise}\\\
\end{cases}$ (8b) $\displaystyle D_{l^{\prime}}^{o}(B_{l^{\prime}}^{o})$
$\displaystyle\leftarrow\begin{cases}B_{l^{\prime}}^{o}+s_{\eta_{l^{\prime}}}+L_{\max}^{\eta_{l^{\prime}}}&\text{if}~{}\eta_{l^{\prime}}=o\\\
D_{l}^{o}(B_{l}^{o})-\max\\{0,B_{l}^{o}+s_{\eta_{l}}+t_{\eta_{l},\eta_{l^{\prime}}}-B_{l^{\prime}}\\}&\text{otherwise}\\\
\end{cases}$ (8c) $\displaystyle B_{l^{\prime}}$
$\displaystyle\leftarrow\begin{cases}b_{\eta_{l^{\prime}}}&\text{if}~{}\eta_{l^{\prime}}\in\mathcal{P}\\\
\min\\{b_{\eta_{l^{\prime}}},D^{\eta_{l^{\prime}}-n}_{l}(B_{l}^{\eta_{l^{\prime}}-n})\\}&\text{if}~{}\eta_{l^{\prime}}\in\mathcal{D}\\\
\end{cases}$ (8d)

The REFs in Eqs.(8a)–(8d) were originally proposed by gschwind2015effective to
handle the maximum ride time constraints. The main idea is to update two
distinct time stamps for each open request on the partial route, for every
label extension, to specify the range of the latest possible drop-off time of
the particular open request. These time stamps are referred to as dynamic time
windows (DTW) proposed by gschwind2015effective. Specifically, in Appendix
LABEL:appendix:description_ride_time_resources, we present two particular
cases to better understand the properties of DTWs: (1) initiating the DTW for
a new open request, and (2) updating the DTW for an existing open request.

####  4.1.3 REF for RDARP

In the labeling algorithm, the risk calibration process is invoked when a
feasible DARP extension is identified along a partial path. The purposes of
the risk calibration, upon the extension, are to identify possible waiting due
to early arrival, compute the maximum available resources for adjusting start-
of-service time, and update the up-to-date value of MMR. To accomplish the
tasks, we will make use of three auxiliary variables as follows:

  * •

$\Delta d_{\eta_{l},\eta_{l^{\prime}}}^{i}$: the start-of-service delay buffer
for request $i$ on label extension $(\eta_{l},\eta_{l^{\prime}})$, which
corresponds to the maximum amount of time that an open or associated request
$i\in\mathcal{O}_{l}\cup\mathcal{O}_{l}^{a}$ can delay its start-of-service at
node $\eta_{l}$, while extending to $\eta_{l^{\prime}}$ without violating the
arrival time and maximum ride time constraints. We use $\Delta
d_{\eta_{l},\eta_{l^{\prime}}}$ without $i$ in the superscript to represent
the set $\\{\Delta
d_{\eta_{l},\eta_{l^{\prime}}}^{i}:i\in\mathcal{O}_{l}\cup\mathcal{O}_{l}^{a}\\}$

  * •

$\Delta\omega_{\eta_{l},\eta_{l^{\prime}}}$: waiting time due to the early
arrival at node $\eta_{l^{\prime}}$ that can be offset by delaying the start-
of-service at node $\eta_{l}$

  * •

$\delta_{\eta_{l},\eta_{l^{\prime}}}^{i}$: the actual amount of start-of-
service delay for each request $i\in\mathcal{O}_{l}\cup\mathcal{O}_{l}^{a}$ to
achieve the MMR condition along the partial route associated with label
$l^{\prime}$. We use $\delta_{\eta_{l},\eta_{l^{\prime}}}$ without $i$ in the
superscript to represent the set
$\\{\delta_{\eta_{l},\eta_{l^{\prime}}}^{i}:i\in\mathcal{O}_{l}\cup\mathcal{O}_{l}^{a}\\}$

With the above notation, for a label $l$ extending from $\eta_{l}$ to
$\eta_{l^{\prime}}$, the anticipated maximum risk
$\BFH_{\eta_{l},\eta_{l^{\prime}}}(\cdot)$ can be minimized following the
Bellman Equation as below:

$\BFH_{\eta_{l},\eta_{l^{\prime}}}(l,\Delta
d_{\eta_{l},\eta_{l^{\prime}}})=\min_{\delta_{\eta_{l},\eta_{l^{\prime}}}}\big{\\{}\Delta\BFH_{\eta_{l},\eta_{l^{\prime}}}(\delta_{\eta_{l},\eta_{l^{\prime}}})+\BFH_{\eta_{l^{\prime}},*}(l^{\prime},\Delta
d_{\eta_{l^{\prime},*}})\big{\\}}$ (9)

where we use $\Delta\BFH_{\eta_{l},\eta_{l^{\prime}}}$ to denote the realized
risk increment to the maximum risk by extending from $\eta_{l}$ to
$\eta^{\prime}_{l}$, and $\BFH_{\eta_{l^{\prime}},*}$ as the future
realization of maximum risk by further extension from $\eta^{\prime}_{l}$ with
remaining delaying resources $\Delta d_{\eta_{l^{\prime},\star}}$. In general,
solving Eq. (9) needs to consider the trade-off between spending resources
(e.g., the allowable start-of-service delay $\Delta
d_{\eta_{l},\eta_{l^{\prime}}}$) for immediate risk reduction and the
potential gain from future extensions. Nevertheless, we can leverage problem-
specific structures to decompose the dependency between current and future
risk realizations based on the following series of propositions.

###### Proposition 4.2

Both $\Delta\BFH_{\eta_{l},\eta_{l^{\prime}}}$ and
$\BFH_{\eta_{l^{\prime}},\star}$ are non-increasing function with respect to
the increase in $\delta_{\eta_{l},\eta_{l^{\prime}}}$.