\section{Introduction}

The home service industry has been rapidly growing worldwide due to emergent changes in family structures, work obligations and extended work hours, aging population, and outspread of chronic diseases, among others.  The United States alone annually spends about $\$82$ billion on home healthcare services and $\$300$ billion on home repairs and maintenance \citep{Lunden, zhan2020home}. In 2018, the global home service market was valued at around $\$282$ billion and is expected to grow by $\$18.9$\% annually from 2019-2026, reaching $\$1133.4$ billion by 2026 (Verified Market Research, 2019). Therefore, the development of efficient home service scheduling and planning tools is essential to support decision-making in all areas of the home service industry.

Home service planning process often includes the following steps. First, customers request their home services on their preferred dates through a booking system. Second, given each date’s requested services, a home service provider assigns the available providers to customers. Third, the home service provider sets the appointment time for each customer and determine each operator’s route. Fourth, the service provider communicates the appointment times to the customers. On the day of service, the operator executes the schedule and visits customers one by one. In this paper, we consider the perspective of a home service provider who needs to make the following planning decisions for each operator and assigned set of customers: (1) how to route the operator, and (2) how to assign appointment times to customers. Henceforth, for simplicity, we use the term ``\textit{routing}'' to refer to routes that specify the operator movement between customers or the sequence of visits on the operator’s schedule, and ``\textit{appointment time}'' to refer to the scheduled start time of the service of a customer. The set of appointment times of customers represent a schedule.

We formally refer to this problem as the home service operator routing and appointment scheduling (HRAS) problem. Specifically, given a set of customers within a service region that an operator needs to serve, the HRAS problem seeks to determine the operator's route and customers' appointment times. We consider both random customers' service time and travel time between customers. If the operator arrives at the customer's location before the scheduled service start time, then the operator needs to wait until the scheduled start time (i.e., remains idle). On the other hand, the customer needs to wait if the operator arrives after the scheduled start time. The operator has fixed service hours beyond which s/he experiences overtime. The quality of routing and scheduling decisions is a function of the total operational cost, consisting of customers' waiting time, and the operator's travel time, idle time, and overtime.

HRAS is a challenging optimization problem for several reasons. First, it is a complex combinatorial optimization problem that requires deciding the operator's route and assigning appointment times to each customer in the route simultaneously. Second, suppose we fix the appointment times. In this case, HRAS becomes similar to the traveling salesman problem (TSP) and the single-vehicle routing problem with time constraints, which are known to be challenging optimization problems (see \citealp{cook2011pursuit}). Third, for a fixed route, HRAS reduces to the well-known challenging  single-server stochastic appointment scheduling (AS) problem with random service time, which seeks a sequence of appointment times that minimizes customers' waiting times and the provider's idling time and overtime \citep{ahmadi2017outpatient, robinson2003scheduling}.  The rationale for considering travel time or travel distance as an additional optimization criterion in HRAS is to minimize service delivery or operational cost by controlling operator's amount of compensation for the distance or time traveled \citep{grieco2020operational,kandakoglu2020decision}.

One can use stochastic programming (SP) to model uncertain travel and service time. However, the applicability of SP is limited to the case in which we know the distribution of uncertainty or we have a large data sample. In most real-world applications, it is challenging to estimate the actual distribution of uncertain parameters accurately, especially with limited data during the planning process. Moreover, future uncertainty realizations may differ from all training samples. If we optimize according to a biased data sample from a misspecified distribution, the resulting optimal routing and scheduling decisions may have a disappointing out-of-sample performance (e.g., excessive overtime) under the true distribution.

Alternatively, one can construct an \textit{ambiguity set} (i.e., a family) of all possible distributions that share some partial information about the random parameter. Then, using the ambiguity set, one can formulate a distributionally robust optimization (DRO) problem to minimize the worst-case (i.e., maximum) expectation of second-stage random cost over all distributions residing within the ambiguity set.  Note that in DRO, the optimization is based on distributions within an ambiguity set, i.e., the distribution is a decision variable in DRO. DRO has recently received substantial attention in various application domains due to the following benefits. First, DRO models acknowledge and hedge against the presence of distributional uncertainty \citep{kuhn2019wasserstein}. Therefore, DRO solutions often faithfully anticipate the possibility of disappointing out-of-sample consequences and depending on the ambiguity set, they often guarantee an out-of-sample cost that falls below the worst-case optimal cost \citep{kuhn2019wasserstein,rahimi2020comprehensive}.

Second, by allowing random parameters to follow a distribution defined in the ambiguity set, DRO alleviates the unrealistic assumption of the decision-maker's complete distributional knowledge. Third, intuitive and easy-to-approximate statistics could be used to construct ambiguity sets  that a decision-maker may approximate and change in the model. For example, decision-makers could estimate the average service and travel time based on their experience, a prediction model, or a  small  data set. Then, one could construct a mean-support ambiguity set on service and travel times where the support could represent the error margin and variability in the estimates that we seek protection against. Fourth, some DRO models of real-world problems are more tractable than their SP counterparts (see, e.g, \citealp{jian2017integer, saif2020data,  shehadeh2020distributionally,  shehadeh2020distributionallyTucker, wang2019distributionally, wang2020distributionally,  zhang2017distributionally}).

The ambiguity set is a key ingredient of DRO models that should capture the true distribution with a high degree of certainty and be computationally manageable (i.e., allow for a tractable DRO model). There are several approaches to construct the ambiguity set. Most applied DRO literature employs moment-based ambiguity,  which incorporates all distributions sharing certain moments (e.g., mean). This is because many tools have been developed to derive tractable and solvable reformulations of moment-based DRO models. In most cases, moment-based DRO models do not scale with the number of scenarios and this leads to a more efficient computational approach. However, moment-based DRO models often do not have asymptotic properties  because the ambiguity sets only incorporates partial distributional information (i.e., moment information).

Recent DRO studies have shifted toward data-driven, distance-based DRO approaches such as $\phi$--divergence \citep{jiang2016data} and Wasserstein distance \citep{esfahani2018data, gao2016distributionally}, that construct ambiguity sets in the vicinity of a reference (e.g., empirical) distribution. One advantage of the Wasserstein ambiguity set is that it could incorporate small-sized data in the ambiguity set and  enjoys asymptotic properties \citep{kuhn2019wasserstein, esfahani2018data}. Recent results indicate that Wasserstein ambiguity set with carefully chosen radius may contain the unknown true distribution with a high probability and is richer than other divergence-based ambiguity sets.  Despite the potential advantages, to the best of our knowledge, there are no moment-based, Wasserstein-based, or other DRO approaches for the HRAS problem we study in this paper. This inspires us to consider this paper's main question: \textit{what are the computational and operational values of employing DRO models to address both service and travel time uncertainty, compared with the classical SP approach for HRAS}?. To answer this question, we design and rigorously analyze an SP and two DRO models based on the mean-support ambiguity and Wasserstein ambiguity sets and demonstrate where significant performance improvements can be gained with each proposed model and derive insights into the problem.

\subsection{Contributions}
In this paper, we address the uncertainty and distributional ambiguity of random travel and service times in HRAS. We summarize our main contributions as follows.

(1)  We propose the first two \textit{\underline{d}istributionally robust \underline{h}ome service \underline{r}outing \underline{a}nd \underline{s}cheduling} (DHRAS) models that seek optimal routing and scheduling decisions to minimize the worst-case sum of the operater's travel time, idle time, and overtime costs, and customers' waiting time costs. We consider two popular types of ambiguity sets for the DRO models, namely, the mean-support (M-DHRAS model) and 1-Wasserstein (W-DHRAS model) ambiguity set.  We also propose SP model for HRAS.

(2) We derive equivalent mixed-integer linear programming (MILP) reformulations of the min-max M-DHRAS and W-DHRAS models. The reformulations can be implemented and solved by off-the-shelf optimization software, thereby enabling practitioners to use these models.

(3) To investigate the value of the DRO approach for HRAS, we conduct an extensive numerical experiment comparing the proposed models operational and computational performance. Our results demonstrate that: (a) W-DHRAS yields robust decisions with both small and large data size, and enjoys both asymptotic consistency and finite-data guarantees; (b) M-DHRAS produces the most conservative schedules and it performs well when the distribution of travel time changes dramatically (e.g., actual travel times are longer); (c) W-DHRAS solutions have better operational performance than the SP solutions both under perfect and misspecified distributions, even when only a small data set is available; (d) DRO models produce more reliable solutions than the SP model; (e) the proposed SP and DRO models are computationally efficient  under realistic settings.

To the best of our knowledge, according to the recent review of operational research in home health care applications \citep{grieco2020operational} and our literature review in Section~\ref{sec2:Literature}, our paper is the first to attempt to investigate and compare theoretically and numerically DRO models for this specific HRAS problem. The remainder of this paper is organized as follows. In Section~\ref{sec2:Literature}, we review the relevant literature. In Section~\ref{sec3:Formulations}, we formally define the DHRAS models, present and analyze the M-DHRAS and W-DHRAS models. In Section~\ref{sec4:Exp}, we present computational results and managerial implications. Finally, we draw conclusions and discuss future directions in Section~\ref{sec5:conclusion}.

\section{Relevant Literature}\label{sec2:Literature}
In this section, we focus primarily on the literature mostly relevant to our problem, i.e., papers that apply stochastic optimization to address HRAS problems similar to ours.  For comprehensive recent surveys of operations research methods applied to decisions in home health care, we refer readers to \cite{fikar2017home,grieco2020operational, gutierrez2013home}. If we fix the appointment time, our problem reduces to the single-vehicle routing problem. That is, we need to design the optimal route of a single vehicle to visit all the customers by minimizing the costs of travel. Recent survey of the literature on these problems includes \cite{costa2019exact, oyola2018stochastic}. On the other hand, if we fix the route, our problem is similar to the single-server appointment scheduling problem. That is, we need to decide the optimal appointment time for the customers to minimize the total cost of customer waiting and server idling and overtime costs (see, e.g., \citealp{berg2014optimal, shehadeh2019analysis}). Appointment times can be fixed time points (see, e.g., \citealp{braekers2016bi, demirbilek2019dynamically, milburn2012operations}) or time windows (see, e.g., \citealp{lee2013continuous, mankowska2014home, yuan2015branch} and references therein). Herein, we focus on the former case under random service and travel times.  We refer to \cite{pinedo2016scheduling} for a comprehensive survey of scheduling theory, applications and methods.

SP approaches for HRAS include \cite{liu2019branch,yuan2015branch,zhan2018vehicle}. The recent work of  \cite{zhan2020home} is closely related to ours. Specifically, \cite{zhan2020home} consider a single operator HRAS problem in which they assume a \textit{deterministic travel time }between customers and that \textit{each customer's service time follows a fully known probability distribution}.  Accordingly,  \cite{zhan2020home} propose a two-stage stochastic mixed-integer linear program,  which seeks first-stage routing and scheduling decisions that minimize the operator traveling cost and an expectation of the second-stage cost of the operator idling and customer waiting. However, they do not incorporate the operator overtime cost in the objective, which is one important source of operational expenses. Note that ignoring random travel time may lead to sub-optimal solutions with, for example, excessive customer waiting time or provider overtime.

The SP approach assumes that the decision-maker is risk-neutral and knows the random parameters' distributions with certainty or can fully estimate them. In practice, it is implausible that decision-makers have sufficient high-quality data to infer the true distribution of random parameters, especially in health care applications. Given distributional ambiguity, if we calibrate a model to a misspecified (i.e., biased) distribution, the resulting optimal SP decisions may have a disappointing out-of-sample performance under the true distribution.  Various robust approaches have been proposed to model the risk-averse nature of decision-makers and uncertain parameters based on partial information of their distributions. The classical robust optimization (RO) approach assumes that uncertain parameters resided on an uncertainty set of possible outcomes with some structure (e.g., polyhedron, see, \citealp{ben2015deriving,bertsimas2004price,soyster1973convex}). Accordingly, optimization in RO is based on the worst-case scenario within the uncertainty set.

Recently, \cite{shi2019robust} propose an RO approach for home health care service.  They consider the problem with multiple caregivers and the goal is to find the optimal route for each caregiver and appointment time for each patient. They propose a robust model on both unknown service and travel time and discuss the heuristic solution approaches, with the objective function being the travel cost and delay cost with respect to the schedules. They find that, as compared with the SP model in \cite{shi2018modeling}, solutions from RO model has a better out-of-sample performance. For instance, the probability of visiting the customers timely is higher.

By focusing the optimization on the worst-case scenario, classical RO approaches often yield over-conservative solutions and suboptimal decisions for other more-likely scenarios \citep{chen2020robust, delage2018value}. Distributionally robust optimization (DRO) is another approach for modeling uncertainty that bridges SP and RO by overcoming their shortcoming and has recently gained significant attention. In DRO, we model the distribution of uncertainty as a decision variable that belongs to an ambiguity set (i.e., a family of all possible distributions of uncertainty). We then optimize based on the distribution within this set \citep{rahimian2019distributionally}.

 There are two common types of ambiguity sets: moment-based and distance-based ambiguity sets (see, e.g., \citealp{rahimian2019distributionally}). In this paper,  we address the distributional ambiguity of service time and travel time in the single-operator HRAS using these two types of ambiguity set, namely the M-DHRAS and W-DHRAS models, in the absence of abundant historical data. We reformulate them into MILPs, which could be solved by off-the-shelf optimization software. To the best of our knowledge and according to the recent surveys \cite{fikar2017home,grieco2020operational}, there is no DRO approach for the specific HRAS problem that we address.  Moment-based DRO models for single-server appointment scheduling includes \citep{jiang2017integer, kong2013scheduling, mak2015appointment}.   Recently, \cite{jiang2019data} propose the distance-based DRO model based on the Wasserstein metric for single-server appointment scheduling for a fixed sequence of customers with random service time and no-show. They derive tractable reformulations of their model under the 1--Wasserstein and 2--Wasserstein ball.

We compare our work with \cite{jiang2019data} and \cite{zhan2020home}, which are recent single-server studies relevant to our work (see Table \ref{table:literature_comparison} in \ref{Appx:LR_Compare} for a quick comparison). First, the three papers address stochastic service time. Second, \cite{jiang2019data} study a classical single-server appointment scheduling problem in which they need to determine the appointment time for a fixed sequence of customers. The server does not need to travel to different customers. Hence, no routing decisions or travel time consideration is necessary. In other words, \cite{jiang2019data}'s model cannot be used for routing and appointment scheduling. Third, \cite{zhan2020home} assume that the travel time is deterministic, which as mentioned earlier, may lead to sub-optimal and unrealistic solutions. In addition,  \cite{zhan2020home} do not address the potential distributional ambiguity of service time while \cite{jiang2019data} address this.  In contrast to \cite{zhan2020home} we consider both stochastic service and travel times. Moreover, we address both uncertainty and distributional ambiguity of these parameters via SP and DRO. Note that if we fix the travel time as zero, our model reduces to a scheduling model. Thus, in this case, our models generalize that of \cite{jiang2019data} by incorporating sequencing decisions and hence, can be used for sequencing and scheduling. Various appointment scheduling studies \citep{ahmadi2017outpatient, berg2014optimal, cayirli2006designing, cayirli2008assessment, creemers2012optimal,  shehadeh2019analysis} demonstrate the benefit of sequencing customers appointments based on their characteristics for improving performance and reducing costs as compared with fixed sequence approach that \cite{jiang2019data} adopt.

\section{Formulation and Analysis}\label{sec3:Formulations}
In this section, we formulate the DHRAS problem using moment  (M-DHRAS) and Wasserstein  (W-DHRAS) ambiguity sets. We define DHRAS and provide a two-stage stochastic programming formulation in Section~\ref{sec3:Defintions}. We analyze and derive equivalent mixed-integer linear programming reformulations of M-DHRAS in Section~\ref{MDHRAS} and W-DHRAS in Section~\ref{WDHRAS}.

\vspace{1mm}
\noindent \textbf{Notation:} For $a,b \in \mathbb{Z}$, we define $[a]:=\{1,2,\dots, a \}$ and $[a,b]_\mathbb{Z}:=\{ c \in \mathbb{Z}: a \leq c \leq b \}$, i.e., $[a,b]_\mathbb{Z}$ represent the set of running integer indices $\{a, a+1, a+2, \dots, b \}$. For a real number $a$, we define $(a)^+=\max\{a,0\}$. We use boldface notation to denote vectors, e.g., $\pmb{d}:=[d_1, d_2, \ldots, d_N]^\top$.

\begin{table}[t]
\small
\center
\renewcommand{\arraystretch}{0.5}
\caption{Notation.}
\begin{tabular}{ll}
\hline
 \multicolumn{2}{l}{\textbf{Indices}} \\
$i$ & index of customer, $i=1,\dots,N$ \\
$j$ & index of service position, $j=1,\dots,N$ \\
\multicolumn{2}{l}{\textbf{Parameters and sets}} \\
$N$ & number of customers \\
$L$ & standard work time \\
$c^{\mbox{\tiny u}}_j$ & idling cost associated to the early arrival at the $j$th customer's location\\
$c^{\mbox{\tiny w}}_j$ & waiting cost of the $j$th customer \\
$c^{\mbox{\tiny o}}$ & overtime cost\\
$\lambda$ & travel time cost \\
$d_i$ & service time of customer $i$ \\
$t_{i,i'}$ & travel time from customers $i$ to $i'$ (with $0$ denoting the depot)\\
$\underline{d}_i$/$\overline{d}_i$ & lower/upper bound of the service time of customer $i$ \\
$\underline{t}_{i,i'}$/$\overline{t}_{i,i'}$ & lower/upper bound of the travel time from customers $i$ to $i'$ \\
\multicolumn{2}{l}{\textbf{First-stage decision variables}} \\
$x_{i,j}$ & binary variable equal to $1$ if customer $i$ is the $j$th customer served and $0$ otherwise \\
$a_j$ & appointment time of the $j$th customer\\
\multicolumn{2}{l}{\textbf{Second-stage decision variables}} \\
$u_j$ & idling time due to an early arrival at the $j$th customer's location \\
$w_j$ & waiting time of the $j$th customer ($w_{N+1}$ as the overtime)\\
\hline
\end{tabular}\label{table:notation}
\end{table}

\subsection{\textbf{Definitions and assumptions}}\label{sec3:Defintions}

We consider a set of $N$ customers that need to be served within a given day by a single home service provider (operator). Traveling time, $t_{i,i'}$, between each pair of customers $i$ and $i'$ ($i,i'=0, \dots, N$) is random, with $i=0$ representing the service provider's office (depot).  Service time $d_i$ of each customer $i \in [N]$ is also random. Note that we may only have partial information on the distributions or possibly small data available on these random parameters. We seek two sets of decisions:  (1) service provider's visiting sequence (routes), (2) customers appointment times (schedule).  The objective is to minimize the sum of customers' waiting time, and the provider's traveling time, idling time, and overtime.

For all $i \in [N]$ and $j \in [N]$, we let the binary decision variable $x_{i,j}$ equal 1 if customer $i$ is the $j$th customer in the operator's route/schedule, and zero otherwise.  For all $j \in [N]$, we let the continuous variable $a_j$ represent the appointment starting time of the $j$th customer.  The feasible region of variable $\pmb{x}$ is defined in \eqref{eq:RegionX} such that each customer is assigned to one position in the operator's route/schedule, and each position is assigned to one customer. The feasible region of  $\pmb{a}$ is defined in \eqref{eq:RegionA} such that all appointments are scheduled within the provider's service hours [$0, L$].
\allowdisplaybreaks
\begin{align}
\mathcal{X}&=\left\{ \pmb{x}:  \begin{array}{l} \sum_{i=1}^N x_{i,j} =1, \  \forall j \in [N] \\
 \sum_{j=1}^N x_{i,j} =1,  \ \forall  i \in [N] \\
x_{i,j} \in \{0, 1 \}, \  \forall i \in [N], \ \forall j \in [N] \end{array} \right\}  \label{eq:RegionX}
\end{align}
\begin{align}
 \mathcal{A}&= \left\{ \pmb{a}: \begin{array}{l}
 \  0 \leq a_j \leq L, \  \ \forall j \in [N] \\
  \ a_{j} \geq a_{j-1}, \ \forall j \in [2, N]_\mathbb{Z} \end{array} \right\} \label{eq:RegionA}
\end{align}

Due to random travel and service times, one or multiple of the following scenarios may happen: (1) operator arrives at the customer's location before the scheduled service start time ($a_j$), and thus s/he remains idle until the scheduled start time; (2)  the operator arrives after the scheduled start time of the customer, and thus the customer incurs waiting cost; (3) provider works overtime beyond his/her scheduled $L$ to finish all appointments. Let the continuous decision variable $w_j$ represent the waiting time of the $j$th customer, for all $j \in [N]$ and $w_{N+1}$ represent the operator's overtime. For all $j \in [N]$, let the continuous decision variable $u_j$ represent the provider's idle time before the start time of the $j$th appointment. For all $i \in [N]$, let the random parameter $d_i$ represent the service duration of customer $i$. For all $i \in [N]$ and $i' \in [N]$, let the random parameter $t_{i,i'}$ represent the travel time between $i$ and $i'$. Given a fixed $\pmb{x} \in \mathcal{X}$,  $\pmb{a} \in \mathcal{A}$ and a joint realization $\pmb{\xi}:=[\pmb{t}, \pmb{d} ]^\top$, we can compute the operational costs (as a function of travel time, idle time, and overtime, and customers waiting time), using the following linear program:
\allowdisplaybreaks
\begin{subequations}\label{LP:Sec}
\begin{align}
  f(\pmb{x}, \pmb{a}, \pmb{\xi} ):=  & \min_{\pmb{u},\,\pmb{w}} \ \Bigg \{  \sum \limits_{j=1}^N (c^{\mbox{\tiny w}}_j w_j+c^{\mbox{\tiny u}}_j u_j \big) +c^{\mbox{\tiny o}} w_{N+1} + \lambda A \Bigg \}  \label{2nd:Obj}\\
& \ \text{s.t.}  \ \ w_1-u_1= \sum \limits_{i=1}^N t_{0,i}x_{i,1}-a_1, \label{2nd:C1}\\
&  \ \ w_j-w_{j-1}-u_j=a_{j-1}-a_j+\sum_{i=1}^N d_ix_{i,j-1}+ \sum_{i=1}^N \sum_{i'\neq i} t_{i,i'} x_{i,j-1}x_{i,j}, \ \forall j\in [2, N]_\mathbb{Z}, \label{2nd:C2}\\
& \ \   w_{N+1}-w_N-u_{N+1}= a_N-a_{N+1}+\sum_{i=1}^N d_ix_{i,N}, \label{2nd:C3}\\
& \ \ A= \sum \limits_{j=2}^N   \sum_{i=1}^N \sum_{i'\neq i} t_{i,i'} x_{i,j-1}x_{i',j}+\sum_{i=1}^N(t_{0,i}x_{i,1}+t_{i,0}x_{i,N}), \label{2nd:C5}\\
& \ \ (w_j,u_j)\geq 0, \ \ \forall j \in [N+1], \label{2nd:C4}
\end{align} \label{2ndStage
\end{subequations}
where $c^{\mbox{\tiny w}}_j$, $c^{\mbox{\tiny u}}_j$, $c^{\mbox{\tiny o}}$, and $\lambda$ are respectively the non-negative unit penalty costs of waiting, idling,  overtime, and travel time for all $j\in [N]$.  Also, we let $a_{N+1}=L$ and $x_{i,N+1}=0$ for all $i\in[N]$. The objective function \eqref{2nd:Obj} minimizes a linear cost function of waiting, idling, overtime, and travel time. Constraint \eqref{2nd:C1}  yields either provider's idle time before the scheduled time of the first appointment or the waiting time of the first customer. Constraint \eqref{2nd:C2} yields either the waiting time of the $j$th customer or the provider's idle time if s/he arrives at the $j$th customer, respectively, after or before the scheduled start time $a_j$ of the $j$th customer.  Constraints \eqref{2nd:C3} yields either the overtime or the schedule earliness. Constraint \eqref{2nd:C5} computes the operator's total travel time.

The SP model in \eqref{SP} seeks to find $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$ that minimizes the expectation of the cost $f(\pmb{x},\pmb{a}, \pmb{\xi}$) subject to random $\pmb{\xi}$ with a known joint probability distribution $\mathbb{P}_{\pmb{\xi}}$.
\begin{align}
& Z^*:= \min \limits_{\pmb{x} \in \mathcal{X}, \pmb{a} \in \mathcal{A}} \mathbb{E}_{\mathbb{P}_{\pmb{\xi}}} [f(\pmb{x},\pmb{a}, \pmb{\xi}) ]  \label{SP}
\end{align}

\subsection{\textbf{DHRAS over moment ambiguity (M-DHRAS)}}\label{MDHRAS}

In this section, we consider the case when $\mathbb{P}_{\pmb{\xi}}$ is not fully known or hard to estimate. We assume that we know or can approximate the mean values $\mathbb{E}_{\mathbb{P}_{\pmb{\xi}}}(\pmb{\xi}):=\pmb{\mu}=[\mub^{\mbox{\tiny d}}, \mub^{\mbox{\tiny t}}]^\top$, lower bounds $[\underline{\pmb{d}},\underline{\pmb{t}}]^\top $ and upper bounds $[\overline{\pmb{d}},\overline{\pmb{t}}]^\top$ of $[\pmb{d},\pmb{t}]^\top$. Mathematically, we consider the support $\mathcal{S}= \mathcal{S}^{\mbox{\tiny d}} \times \mathcal{S}^{\mbox{\tiny t}}$ of $\pmb{\xi}$, where $\mathcal{S}^{\mbox{\tiny d}}$ and $\mathcal{S}^{\mbox{\tiny t}}$ are respectively the supports of random parameters $\pmb{d}$ and $\pmb{t}$ defined in \eqref{Support}.
\begin{subequations} \label{Support}
\begin{align}
& \mathcal{S}^{\mbox{\tiny d}}:=\left\{\, \pmb{d}  \geq 0: \begin{array}{l} \underline{d}_{i} \leq d_{i} \leq \overline{d}_{i},  \ \forall i \in [N]\end{array} \right\}\,,\\
& \mathcal{S}^{\mbox{\tiny t}}:=\left\{\, \pmb{t}  \geq 0: \begin{array}{l} \underline{t}_{i, i'} \leq t_{i, i'} \leq \overline{t}_{i,i'},  \ \forall i \in [0,N],\, i'\in[0,N],\, i' \neq i\end{array} \right\}\,.
\end{align}
\end{subequations}
The assumption of known mean and support is motivated by the fact that experts may be able to approximate the average service or travel time based on their experience (or from prediction models). The range can represent the error margin of the mean values estimates and variability of the parameters which we seek protection against. Experts may also be able to provide an upper and lower bound on service time. Therefore, we can use such information to define a mean-support ambiguity set of uncertain parameters distributions. Using $\pmb{\mu}$ and $\mathcal{S}$, we consider the following mean-support ambiguity set  $\mathcal{F}(\mathcal{S}, \pmb{\mu})$:
\begin{align}\label{eq:Mambiguity}
\mathcal{F}(\mathcal{S}, \pmb{\mu}) := \left\{ \mathbb{P} \in \mathcal{P}(\mathcal{S}): \begin{array}{l} \int_\mathcal{S} d\mathbb{P} = 1\\ \mathbb{E_P}(\pmb{\xi}) = \pmb{\mu} \end{array} \right\},
\end{align}
where $\mathcal{P}(\mathcal{S})$ represents the set of all probability distributions supported on $\mathcal{S}$. Using ambiguity set $\mathcal{F}(\mathcal{S},\pmb{\mu})$, we formulate M-DHRAS as follows:
\begin{align}
(\text{M-DHRAS}) &\ \ \ \min \limits_{\pmb{x} \in \mathcal{X}, \pmb{a} \in \mathcal{A}} \bigg \{ \sup_{\mathbb{P} \in \mathcal{F}(\mathcal{S},\pmb{\mu}) }\mathbb{E}_{\mathbb{P}} [f(\pmb{x},\pmb{a}, \pmb{\xi}) ] \bigg\}. \label{M-DHRAS}
\end{align}
Formulation \eqref{M-DHRAS} seeks first stage decisions ($\pmb{x}, \pmb{a}$) that minimizes the worst-case (maximum) expectations of the second-stage operational cost over $\mathcal{F}(\mathcal{S}, \pmb{\mu})$.

\begin{rem}
Model \eqref{M-DHRAS} could hedge against distributional uncertainty, particularly when limited (or no) information on random parameters is available. The use of mean-support ambiguity set is prevalent due to the following reasons. First, only intuitive statistics, namely the mean and the range, are required. They can be interpreted easily as mean represents the distribution centrality while the range represents the dispersion of the distribution. Decision-makers could approximate these values in various ways (e.g., estimation from the limited available data or from expert knowledge). Second, mathematically, it allows us to derive a tractable MILP reformulation. Various studies have demonstrated that incorporating higher moments or random parameters often undermines the computational tractability of the DRO model due to the non-linearity in higher moments.

\end{rem}

\subsubsection{\textbf{MILP reformulation of M-DHRAS}}\label{sec:ReforM_DHRAS}

Note that model \eqref{M-DHRAS} is not directly solvable in the presented form due to the minimax structure. In this section, we use the recourse problem properties to derive an equivalent MILP reformulation of model \eqref{M-DHRAS}. First, for a fixed $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$, we rewrite the inner maximization problem $\sup_{\mathbb{P} \in \mathcal{F}(\mathcal{S},\pmb{\mu}) }\mathbb{E}_{\mathbb{P}} [f(\pmb{x},\pmb{a}, \pmb{\xi}) ]$ as follows:
\begin{subequations}\label{MDHRAS:InnerMax}
\begin{align}
& \max   \ \mathbb{E_P}[f(\pmb{x},\pmb{a}, \pmb{\xi})]  \\
& \ \text{s.t.}  \ \ \mathbb{E_P}(\pmb{\xi}) = \pmb{\mu}, \label{MDHRAS::InnerMax-C1}  \\
& \ \ \ \ \  \ \ \mathbb{E_P}[\mathds{1}_\mathcal{S}(\pmb{\xi})] = 1, \label{MDHRAS::InnerMax-C2}
\end{align}
\end{subequations}
where $\mathds{1}_\mathcal{S}(\pmb{\xi})= 1$ if $\pmb{\xi} \in \mathcal{S}$ and $\mathds{1}_\mathcal{S}(\pmb{\xi}) = 0$ if $\pmb{\xi} \notin \mathcal{S}$.  As we show in the proof of Proposition~\ref{Prop1M:DualMinMax} in \ref{Appx:ProofOfProb1M}, problem \eqref{MDHRAS:InnerMax} is equivalent to the deterministic problem \eqref{eq:FinalDualInnerMax-1}.

\begin{prop}\label{Prop1M:DualMinMax}
For any $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$, problem \eqref{MDHRAS:InnerMax} is equivalent to
\begin{align} \label{eq:FinalDualInnerMax-1}
&  \min_{\pmb{\alpha},\,\pmb{\rho}}  \ \Bigg \{ \sum \limits_{i=1}^N \mu_i^{\mbox{\tiny d}} \rho_i+ \sum_{i=1}^N \sum_{i'=1}^N \mu_{i,i'}^{\mbox{\tiny t}}  \alpha_{i,i'} +\max \limits_{\pmb{\xi} \in \mathcal{S}  }  \Bigg \{  f(\pmb{x},\pmb{a}, \pmb{\xi})- \sum \limits_{i=1}^N d_i \rho_i- \sum_{i=1}^N \sum_{i'=1}^N t_{i,i'} \alpha_{i,i'}  \Bigg \}  \Bigg\} \\
& \ \ \textup{s.t.} \ \ \pmb{\alpha}\in\mathbb{R}^{(N+1)\times(N+1)},\, \pmb{\rho}\in\mathbb{R}^N. \nonumber
\end{align}
\end{prop}

 Note that $ f(\pmb{x},\pmb{a}, \pmb{\xi})$  is a minimization problem, and thus in \eqref{eq:FinalDualInnerMax-1} we have an inner max-min problem, which is not suitable to solve using standard solution methods. For a given solution $(\pmb{x},\pmb{a})$ and realized value of $\pmb{\xi}$, $ f(\pmb{x},\pmb{a}, \pmb{\xi})$ is a linear program (LP). The dual of $  f(\pmb{x},\pmb{a}, \pmb{\xi})$ is as follow

\begin{subequations}\label{Dual:Sec}
\begin{align}
  f(\pmb{x}, \pmb{a},\pmb{\xi} ):=  & \max_{\pmb{y}} \Bigg \{ \Bigg(\sum_{i=1}^N t_{0,i}x_{i,1}-a_1\Bigg)y_1+  \sum_{j=2}^{N+1} \Bigg ( a_{j-1}-a_j+\sum_{i=1}^N d_ix_{i,j-1}+ \sum_{i=1}^N \sum_{i'\neq i} t_{i,i'} x_{i,j-1}x_{i,j}\Bigg)y_j  \nonumber\\
  & \qquad \ \ \  \ + \lambda \Bigg[\sum \limits_{j=2}^N   \sum_{i=1}^N \sum_{i'\neq i} t_{i,i'} x_{i,j-1}x_{i',j}+\sum_{i=1}^N\Big(t_{o,i}x_{0,i}+t_{i,0}x_{i,N}\Big) \Bigg] \Bigg \}  \label{Dual2nd:Obj}\\
  & \text{s.t.}  \ \  \pmb{y}\in\mathcal{Y}:= \Big\{\pmb{y} \mid  0\leq y_{N+1} \leq c^{\mbox{\tiny o}}, \ -c^{\mbox{\tiny u}}_j\leq y_j \leq c^{\mbox{\tiny w}}_j+ y_{j+1}, \ \forall j \in [N]\Big\}, \label{Dual2nd:C1}
\end{align} \label{Dual2ndStage}
\end{subequations}
\hspace{-1.2mm}where $y_1, \ldots, y_{N+1}$ are the dual variables associated with constraints \eqref{2nd:C1}--\eqref{2nd:C3}.  Given $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$ and $\pmb{\xi}$, the objective function in \eqref{Dual2nd:Obj} is linear (convex) in $\pmb{y}$. Hence, the inner maximization problem in \eqref{Dual:Sec} is a convex maximization problem. It follows from the fundamental convex analysis that it suffices to consider the extreme points of polytope $\mathcal{Y}$. This motivates us to leverage the properties of the extreme points of $\mathcal{Y}$ in deriving an equivalent LP reformulation of \eqref{Dual:Sec}. We formally prove this in Proposition~\ref{Prop2M} (see \ref{Appx:ProofProp2M} for a detailed proof).

\begin{prop}\label{Prop2M}
For any $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$, the inner maximization problem in \eqref{eq:FinalDualInnerMax-1} is equivalent to .
\begin{subequations}\label{DHRASMINLP}
\begin{align}
 \min_{\pmb{\beta},\,\pmb{\gamma},\,\pmb{\delta}}  &\  \sum \limits_{j=1}^{N+2} \beta_j  +\sum_{i=1}^N\left[(\lambda-\alpha_{i,0})\underline{t}_{i,0}+\Delta t_{i,0}(\lambda-\alpha_{i,0})^+\right]x_{i,N}  \label{DHRAS_MINLP_Obj}\\
\textup{s.t.} & \ \beta_1 \geq -a_1  \pi_{1,1} +\sum_{i=1}^N \big [\underline{t}_{0,i} (\pi_{1,1}+\lambda-\alpha_{0,i})+\Delta t_{0,i} \gamma_{0,i,1,1}\big] x_{i,1} , \label{DHRASMINLP:Con1}\\
& \ \sum_{j=1}^v \beta_j \geq -a_1  \pi_{1,v} + \sum_{j=2}^{\min(v,N+1)} \big ( a_{j-1}-a_j \big) \pi_{j,v}   + \sum_{i=1}^N \big [\underline{t}_{0,i} (\pi_{1,v}+\lambda-\alpha_{0,i})+\Delta t_{0,i} \gamma_{0,i,1,v}\big] x_{i,1}  \nonumber \\
& \ \qquad \   + \sum_{j=2}^{\min(v,N+1)} \sum_{i=1}^N \sum_{i' \neq i} \big[ \underline{t}_{i,i'} (\pi_{j,v}+\lambda-\alpha_{i,i'})+ \Delta t_{i,i'} \gamma_{i,i',j,v}  \big]x_{i,j-1}x_{i', j} \nonumber\\
& \ \qquad \   + \sum_{j=2}^{\min(v,N+1)} \sum_{i=1}^N \big[ \underline{d}_i (\pi_{j,v}-\rho_i)+ \Delta d_i \delta_{i,j,v} \big ] x_{i,j-1}, \ \
 \forall v \in [2, N+2]_{\mathbb{Z}} ,  \label{DHRASMINLP:Con2} \\
& \ \sum_{j=k}^v\beta_j \geq \sum_{j=k}^{\min(v,N+1)} \big ( a_{j-1}-a_j\big)  \pi_{j,v}+\sum_{j=k}^{\min(v,N+1)} \sum_{i=1}^N \big[ \underline{d}_i (\pi_{j,v}-\rho_i)+ \Delta d_i \delta_{i,j,v} \big ] x_{i,j-1} \nonumber\\
& \ \qquad \    \ + \sum_{j=k}^{\min(v,N+1)} \sum_{i=1}^N \sum_{i' \neq i}  \bigg\{ (\pi_{j,v}+\lambda-\alpha_{i,i'})\underline{t}_{i,i'}+\Delta t_{i,i'}\gamma_{i,i',j,v}  \bigg\}x_{i,j-1}x_{i',j},   \nonumber \\
& \ \qquad \    \  \forall k\in[2,N+1]_\mathbb{Z},\,\forall v\in[k,N+2]_\mathbb{Z} , \label{DHRASMINLP:Con3} \\
& \ \beta_{N+2} \geq 0 , \label{DHRASMINLP:Con4} \\
& \ \gamma_{0,i,1,v}\geq 0,\quad \gamma_{0,i,1,v}\geq \pi_{1,v}+\lambda-\alpha_{0,i},\quad \forall i\in[N],\,\forall v\in[N+2] ,\label{DHRASMINLP:Con5} \\
& \ \gamma_{i,i',j,v}\geq 0,\quad \gamma_{i,i',j,v}\geq \pi_{j,v}+\lambda-\alpha_{i,i'},\quad \forall i\in[N],\,\forall i'\in[N]\setminus\{i\}, \nonumber \\
& \hspace{70mm} \forall j\in[2,N+1]_\mathbb{Z},\,\forall v\in[j,N+2]_\mathbb{Z} , \label{DHRASMINLP:Con6}  \\
& \ \delta_{i,j,v}\geq 0,\quad \delta_{i,j,v}\geq \pi_{j,v} -\rho_i,\quad\forall i\in[N],\,\, \forall j\in[2,N+1]_\mathbb{Z},\,\forall v\in[j,N+2]_\mathbb{Z}  ,\label{DHRASMINLP:Con7}
\end{align}
\end{subequations}
where $\pi_{j,v}=-c_v^u+\sum_{l=j}^{v-1}c_l^w$ for $1\leq j\leq v\leq N+2$, Let $\Delta d_i=\overline{d}_i-\underline{d}_i$, and $\Delta t_{i,i'}=\overline{t}_{i,i'}-\underline{t}_{i,i'}$.
\end{prop}

With the use of Proposition \ref{Prop2M}, we can combine the minimization over $\pmb{\alpha}$ and $\pmb{\rho}$ in problem \eqref{eq:FinalDualInnerMax-1}. Corollary $\ref{Coro3}$ provides the reformulation of problem \eqref{eq:FinalDualInnerMax-1}  (see \ref{apdx: Coro_proof} for a proof).
\begin{cor} \label{Coro3}
Problem \eqref{eq:FinalDualInnerMax-1} is equivalent to
\begin{subequations}\label{InnerMINLPnew}
\begin{align}
 \min_{\pmb{\alpha},\,\pmb{\rho},\,\pmb{\beta},\,\pmb{\gamma},\,\pmb{\delta}} & \  \sum \limits_{i=1}^N \mu_i^{\mbox{\tiny d}} \rho_i+ \sum_{j=2}^N\sum_{i=1}^N \sum_{i'\ne i} \mu^{\mbox{\tiny t}}_{i,i'}  \alpha_{i,i'}x_{i,j-1}x_{i',j}  + \lambda \sum \limits_{i=1}^N \mu^{\mbox{\tiny t}}_{i,0}  x_{i,N} + \sum \limits_{i=1}^N \mu^{\mbox{\tiny t}}_{0,i} \psi_{0,i}  + \sum \limits_{j=1}^{N+2} \beta_j  \label{InnerMINLPnew:Obj}\\
\textup{s.t. \hspace{3mm}} & \ \pmb{\alpha}\in\mathbb{R}^{(N+1)\times(N+1)},\, \pmb{\rho}\in\mathbb{R}^N,\, \pmb{\beta}\in\mathbb{R}^{N+2},\, \eqref{DHRASMINLP:Con1} - \eqref{DHRASMINLP:Con7}
\end{align}
\end{subequations}
\end{cor}

Combining  $\sup_{\mathbb{P} \in \mathcal{F}(\mathcal{S},\pmb{\mu}) }\mathbb{E}_{\mathbb{P}} [f(\pmb{x},\pmb{a}, \pmb{d}) ]$  in the form of \eqref{InnerMINLPnew} with the outer minimization problem in \eqref{M-DHRAS}, we obtain a mixed-integer nonlinear program (MINLP).  Due to its large size and page limitation, we present this MINLP in \ref{Appx:Mac_MDHRAS}, where we also derive the following exact and tight MILP reformulation of this this MINLP (equivalently, M-DHRAS model).
\allowdisplaybreaks
\begin{subequations}\label{M_DHRAS_L}
\begin{align}
 \min & \  \sum \limits_{i=1}^N \mu^{\mbox{\tiny d}}_i \rho_i+ \sum_{j=2}^N\sum_{i=1}^N \sum_{i'\ne i} \mu^{\mbox{\tiny t}}_{i,i'} \eta_{i,i',j-1,j}  + \lambda \sum \limits_{i=1}^N \mu^{\mbox{\tiny t}}_{i,0} x_{i,N} + \sum \limits_{i=1}^N \mu^{\mbox{\tiny t}}_{0,i} \alpha_{0,i} x_{i,1}  + \sum \limits_{j=1}^{N+2} \beta_j  \label{M_DHRAS_L:obj}\\
 \text{s.t.}  & \ \pmb{x} \in \mathcal{X}, \ \pmb{a} \in \mathcal{A},\, \pmb{\alpha}\in\mathbb{R}^{(N+1)\times(N+1)},\, \pmb{\rho}\in\mathbb{R}^N,\, \pmb{\beta}\in\mathbb{R}^{N+2}, \\
& \ \text{constraints }\eqref{M-DHRAS_MILP_Con2}-\eqref{M-DHRAS_MILP_Con5},\,\eqref{DHRASMINLP:Con5}-\eqref{DHRASMINLP:Con7},\,\eqref{MILP-MAC1M}-\eqref{MILP-MAC11M}.
\end{align}
\end{subequations}

\subsection{\textbf{DHRAS over 1-Wasserstein ambiguity (W-DHRAS)}}\label{WDHRAS}

In this section, we consider the case that $\mathbb{P}_{\pmb{\xi}}$ may be observed via a  small finite set $\big\{\hat{\pmb{\xi}}^1, \ldots, \hat{\pmb{\xi}}^R\big\}$ of $R$ i.i.d. samples, which may come from the limited historical realizations or a reference empirical distribution. Note that the empirical distribution defined on this data (or any other reference distribution) could serve as an estimator of the true distribution. To hedge against the estimation error, we robustify the optimization problem against all distributions $\mathbb{P}_{\pmb{\xi}}$ close to the empirical distribution $\mathbb{P}_{\pmb{\xi}}^R$.  Specifically, we construct a set of all distributions $\mathcal{P}(\mathcal{S})$ supported on the support $\mathcal{S}$ of the unknown true distribution for which all distributions $\mathbb{P}_{\pmb{\xi}}$ have 1-Wasserstein distance that is less than or equal to $\epsilon$ from the reference distribution. We use 1-Wasserstein distance (i.e., we use $\ell_1$--norm in the definition of Wasserstein metric) because it often admits tractable reformulation in most real-world applications (see, e.g., \citealp{Daniel2020, hanasusanto2018conic, jiang2019data, saif2020data}). Next, we present our W-DHRAS model.

  Suppose that two probability distributions $\mathbb{Q}_1$ and $\mathbb{Q}_2$ are defined on a common support $\mathcal{S} \in \mathbb{R}^N$, and let $|| \cdot ||_p$ represent the $p$-norm on $\mathbb{R}^N$ with $p  \geq 1$. Suppose that random vectors $\pmb{\xi}_1$ and $\pmb{\xi}_2$ follow $\mathbb{Q}_1$ and $\mathbb{Q}_2$ respectively. Then, the Wasserstein distance between $\mathbb{Q}_1$ and $\mathbb{Q}_2$, denoted as $W_p(\mathbb{Q}_1, \mathbb{Q}_2)$, represents the minimum transportation cost of moving from $\mathbb{Q}_1$ to $\mathbb{Q}_2$, where the cost of moving from $\pmb{\xi}_1$ to $\pmb{\xi}_2$ is measured by the norm $\norms{\pmb{\xi}_1-\pmb{\xi}_2}_p$. Mathematically,
\begin{equation}
W_p(\mathbb{Q}_1, \mathbb{Q}_2):= \bigg ( \inf_{\Pi \in \mathcal{P}(\mathbb{Q}_1, \mathbb{Q}_2)} \mathbb{E}_{\Pi} \big[ \norm{\pmb{\xi}_1 -\pmb{\xi}_2}_p^p\big]  \bigg)^\frac{1}{p},
\end{equation}
where $\mathcal{P}(\mathbb{Q}_1, \mathbb{Q}_2) $ is the set of all joint distributions of ($\pmb{\xi}_1$, $\pmb{\xi}_2$) with marginals $\mathbb{Q}_1$ and $\mathbb{Q}_2$. Since we only observe a  set $\{\hat{\pmb{\xi}}^1, \ldots, \hat{\pmb{\xi}}^R\}$ of $R$ i.i.d. samples, we consider the following $p$-Wasserstein ambiguity set
\begin{align}\label{c}
\mathcal{F}_p (\hat{\mathbb{P}}_{\pmb{\xi}}^R, \epsilon)= \left\{ \mathbb{Q}_{\pmb{\xi}} \in \mathcal{P}(\mathcal{S}): W_p(\mathbb{Q}_{\pmb{\xi}}, \hat{\mathbb{P}}_{\pmb{\xi}}^R) \leq \epsilon\right\},
\end{align}
where $ \mathcal{P}(\mathcal{S})$ is the set of all probability distributions on $\mathcal{S}$, $\hat{\mathbb{P}}_{\pmb{\xi}}^R=\frac{1}{R} \sum_{r=1}^R \delta_{\hat{\pmb{\xi}}^r}$ is the empirical distribution of $\pmb{\xi}$ based on the $R$ i.i.d samples with $\delta$ being the Dirac measure, and $\epsilon >0$ is the radius of the ambiguity set. The set $\mathcal{F}_p (\hat{\mathbb{P}}_{\pmb{\xi}}^R, \epsilon)$ can be viewed as the $p$-Wasserstein ball of radius $\epsilon$ centered at the empirical distribution $\hat{\mathbb{P}}_{\pmb{\xi}}^R$. Therefore, in some sense, one can think Wasserstein ball as the set of all distributions under which our estimation error is below $\epsilon$, where $\epsilon$ is the estimation error we seek protection against. A larger radius $\epsilon$ indicates that we seek more robust solutions. Using the ambiguity set $\mathcal{F}_p(\hat{\mathbb{P}}_{\pmb{\xi}}^R, \epsilon)$, we formulate W-DHRAS as follows:
\begin{align}
(\text{W-DHRAS}) \ \ \ \hat{Z}(R, \epsilon)&= \min \limits_{\pmb{x} \in \mathcal{X}, \pmb{a} \in \mathcal{A}} \Bigg \{ \sup_{\mathbb{Q}_{\pmb{\xi}} \in \mathcal{F}_p (\hat{\mathbb{P}}_{\pmb{\xi}}^R, \epsilon) }\mathbb{E}_{\mathbb{P}_{\pmb{\xi}}} \big[f(\pmb{x},\pmb{a}, {\pmb{\xi}}) \big] \Bigg \}. \label{W-DHRAS}
\end{align}

In data-driven approaches such as W-DHRAS, we often seek asymptotic consistency. Specifically, we expect that as the sample size $R$ increases to infinity, the optimal value of the problem \eqref{W-DHRAS} $\hat{Z}(R, \epsilon)$ converges to $Z^*$ (the optimal value of the  SP model in \eqref{SP} with perfect knowledge of $\mathbb{P}_{\pmb{\xi}}$), and an optimal solution $(\pmb{x},\pmb{a})$ of  (W-DHRAS) converges to an optimal solution of problem \eqref{SP}. Additionally, if $\hat{Z}(R, \epsilon)>Z^*$ almost surely, then W-DHRAS provides a safe upper bound guarantee on the expected total cost with any finite data size $R$. Recall that our support set $\mathcal{S}$ is non-empty, convex and compact. As such, we can use the existing theory in establishing the asymptotic consistency and finite sample guarantee of our W-DHRAS.  We refer the readers to  \ref{Proof_Lemma1}--\ref{Proof_Thrm2} for adapted proofs of these results.

\subsubsection{\textbf{MILP reformulation of W-DHRAS}}
In this section, we derive an equivalent MILP reformulation of our W-DHRAS model. First, we consider the inner maximization problem of W-DHRAS for a fixed $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$.
\begin{align}
\sup_{\mathbb{Q}_{\pmb{\xi}}\in \mathcal{F}_1 (\hat{\mathbb{P}}_{\pmb{\xi}}^R, \epsilon) }\mathbb{E}_{\mathbb{Q}_{\pmb{\xi}}} [f(\pmb{x},\pmb{a}, \pmb{\xi}) ] \label{InnerMaxW}
\end{align}
 In Proposition~\ref{Prop1}, we present an equivalent dual formulation of  \eqref{InnerMaxW}  (see \ref{Proof_Prop1} for a  proof).
\begin{prop}\label{Prop1}
The optimal value of formulation  \eqref{InnerMaxW} equals that of the following formulation:
\begin{equation}\label{DualOfInner}
\inf_{\rho\geq 0} \Bigg\{ \epsilon \rho + \frac{1}{R} \sum_{r=1}^R \sup_{\pmb{\xi} \in \mathcal{S}}  \bigg\{ f(\pmb{x}, \pmb{a}, \pmb{\xi})-\rho \norms{\pmb{\xi} -\hat{\xib}^r }_1  \bigg\}  \Bigg\}
\end{equation}
\end{prop}

Note that $f(\pmb{x}, \pmb{a}, \pmb{\xi})-\rho \norms{\pmb{\xi} -\hat{\xib}^r}_1$ is neither convex nor concave in $\pmb{\xi}$. Thus, formulation \eqref{DualOfInner} is potentially challenging to solve because it requires solving $R$ non-convex  optimization problems. Given that the support of service duration and travel times are rectangular and finite, we next show that we can recast these problems as linear programs for fixed $\rho$ and $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$. In what follows, we provide a high-level road map of the reformulation, and we relegate the notationally heavy proofs and details to Appendices. First, for fixed $(\pmb{x},\pmb{a})\in\mathcal{X}\times\mathcal{A}$ and $\rho \geq 0$, we denote $g_r(\rho, \pmb{x}, \pmb{a})= \sup_{\xi \in \mathcal{S}}  \{ f(\pmb{x}, \pmb{a}, \pmb{\xi})-\rho \norms{ \pmb{\xi} -\hat{\pmb{\xi}}^r }_1  \}$. Given the dual formulation of $f(\pmb{x}, \pmb{a}, \pmb{\xi})$ in \eqref{Dual2ndStage}, and using the dual of $f(\pmb{x},\pmb{a},\pmb{\xi})$ and using the same proof techniques in Proposition \ref{Prop2M}, we can reformulate $g_r(\rho, \pmb{x}, \pmb{a})$ (see \ref{Proof_Prop4M} for a detailed proof).

\begin{prop}\label{Prop4M}
The maximization problem $g_r(\rho, \pmb{x}, \pmb{a})$ is equivalent to
\allowdisplaybreaks
\begin{subequations}\label{W_inner_r_LP}
\begin{align}
 \min_{\pmb{\beta}^r,\,\pmb{u}^r,\,\pmb{\nu}^r}  &\  \sum \limits_{j=1}^{N+2} \beta^r_j + \sum_{i=1}^N u_{i,0}^r x_{i,N}  \\
\textup{s.t.    } & \ \beta^r_1 \geq -a_1  \pi_{1,1} +\sum_{i=1}^N u_{0,i,1,1}^r x_{i,1} , \label{W_inner_r_LP:Con1} \\
& \ \sum_{j=1}^v \beta^r_j \geq -a_1  \pi_{1,v} + \sum_{j=2}^{\min(v,N+1)} \big ( a_{j-1}-a_j \big) \pi_{j,v} + \sum_{i=1}^N u_{0,i,1,v}^r x_{i,1} \nonumber \\
& \ \qquad \   + \sum_{j=2}^{\min(v,N+1)} \sum_{i=1}^N \sum_{i' \neq i} u_{i,i',j,v}^r x_{i,j-1}x_{i', j}  + \sum_{j=2}^{\min(v,N+1)} \sum_{i=1}^N \nu_{i,j,v}^r x_{i,j-1}, \ \
 \forall v \in [2, N+2]_{\mathbb{Z}},   \label{W_inner_r_LP:Con2} \\
& \ \sum_{j=k}^v\beta^r_j \geq \sum_{j=k}^{\min(v,N+1)} \big ( a_{j-1}-a_j\big)  \pi_{j,v}+\sum_{j=k}^{\min(v,N+1)} \sum_{i=1}^N \nu_{i,j,v}^r x_{i,j-1}  \nonumber \\
& \ \qquad \   + \sum_{j=k}^{\min(v,N+1)} \sum_{i=1}^N \sum_{i' \neq i} u_{i,i',j,v}x_{i,j-1}x_{i',j},\quad  \forall k\in[2,N+1]_\mathbb{Z},\,\forall v\in[k,N+2]_\mathbb{Z}  ,\label{W_inner_r_LP:Con3} \\
& \ \beta^r_{N+2} \geq 0 , \label{W_inner_r_LP:Con4} \\
& \ u^r_{i,0} \geq \lambda\hat{t}^r_{i,0},\quad  u^r_{i,0} \geq \lambda\overline{t}_{i,0}-\rho(\overline{t}_{i,0}-\hat{t}^r_{i,0}),\quad\forall i\in[N], \label{W_inner_r_LP:Con5}  \\
& \ u_{0,i,1,v}^r  \geq (\pi_{1,v}+\lambda)\underline{t}_{0,i}-\rho(\hat{t}_{0,i}^r-\underline{t}_{0,i}),\quad  u_{0,i,1,v}^r  \geq (\pi_{1,v}+\lambda)\hat{t}_{0,i}^r,\nonumber\\
& \ u_{0,i,1,v}^r  \geq (\pi_{1,v}+\lambda)\overline{t}_{0,i}-\rho(\overline{t}_{0,i}-\hat{t}^r_{0,i}),\quad \forall i\in[N],\, v\in[N+2], \label{W_inner_r_LP:Con7} \\
& \ u_{i,i',j,v}^r   \geq (\pi_{j,v}+\lambda)\underline{t}_{i,i'}-\rho(\hat{t}_{i,i'}^r-\underline{t}_{i,i'}),\quad u_{i,i',j,v}^r   \geq (\pi_{j,v}+\lambda)\overline{t}_{i,i'}-\rho(\overline{t}_{i,i'}-\hat{t}^r_{i,i'}), \nonumber\\
& \ u_{i,i',j,v}^r   \geq (\pi_{j,v}+\lambda)\hat{t}_{i,i'}^r,\quad \forall i\in[N],\, i'\in[N]\setminus\{i\},\, j\in[2,N+1]_\mathbb{Z},\, v\in[j,N+2]_\mathbb{Z} ,\label{W_inner_r_LP:Con9} \\
& \ \nu_{i,j,v}^r    \geq \pi_{j,v}\underline{d}_{i}-\rho(\hat{d}_{i}^r-\underline{d}_{i}),\quad
 \nu_{i,j,v}^r    \geq \pi_{j,v}\overline{d}_{i}-\rho(\overline{d}_{i}-\hat{d}^r_{i}), \nonumber \\
& \ \nu_{i,j,v}^r    \geq \pi_{j,v}\hat{d}_{i}^r,\quad \forall i\in[N],\, j\in[2,N+1]_\mathbb{Z},\, v\in[j,N+2]_\mathbb{Z}. \label{W_inner_r_LP:Con11}
\end{align}
\end{subequations}
\end{prop}

Summing $g_r(\rho,\pmb{x},\pmb{a})$ in the form of Proposition \ref{Prop4M} over $r$ and combining it with the outer minimization in \eqref{DualOfInner} and \eqref{W-DHRAS}, we derive MINLP reformulation of the W-DHRAS. Again due to its large size, we present this MINLP in \ref{Appx:Mac_WDHRAS}, where we also derive the following tight and exact MILP reformulation of the W-DHRAS model.
\allowdisplaybreaks
\begin{subequations}\label{W_DHRAS_L}
\begin{align}
 \min &\ \ \epsilon \rho + \frac{1}{R} \sum_{r=1}^R \left(\sum_{i=1}^N \psi^r_{i,0} + \sum_{j=1}^{N+2}\beta_j^r \right)  \\
\text{s.t.} &\ \  \pmb{x}\in\mathcal{X},\,\pmb{a}\in\mathcal{A},\, \rho\geq 0,\, \pmb{\beta}^r\in\mathbb{R}^{N+2},\quad\forall r\in[R], \\
& \ \ \text{constraints }\eqref{WMILP-MAC1M} -\eqref{WMILP-MAC9M}, \,\eqref{W_inner_r_LP:Con5}-\eqref{W_inner_r_LP:Con11},\, \eqref{W-DHRAS_MILP_Con2}-\eqref{W-DHRAS_MILP_Con5},\, \forall r\in[R]. \label{W_DHRAS_L:Con2}
\end{align}
\end{subequations}

\begin{rem}
The use of Wasserstein ambiguity is useful in modeling uncertainty with limited data. In this work, we only consider the $1$-Wasserstein ambiguity set (i.e., $p=1$) since for general $p>2$, it usually results in a non-linear and complex reformulation. In the HRAS problem, routing decision is binary and hence, it could result in mixed-integer non-linear formulation that is more computationally expensive than MILP in general. Also, it is not reasonable to consider $p=\infty$ since the term $\norm{\pmb{\xi}_1-\pmb{\xi}_2}_\infty$ only measure the maximum difference of one particular entry. However, in our problem, every random parameter plays a role in the optimization model and differences in each service and travel time should be considered.
\end{rem}

\section{Numerical Experiments} \label{sec4:Exp}
Our computational study's primary objective is to compare the performance of the proposed DRO models (M-DHRAS and W-DHRAS) and a sample average approximation (SAA) model, and derive insights into HRAS. The SAA model solves model \eqref{SP} with $\mathbb{P}_{\pmb{\xi}}$ replaced by an empirical distribution based on $N$ samples of the random parameters (see \ref{apdx: SAA} for the formulation). For simplicity, we call the SAA model as the SP model. We focus on HRAS instances where the sample size is possibly small or we do not have enough data to model distributions accurately, which is often seen in healthcare applications. However, we also test the computational performance of SP under a large sample size. In Section \ref{subsection:Expt_Setup}, we describe the set of HRAS instances that we constructed and discuss other experimental setup. In Section \ref{subsection:Wass_Eps}, we examine the choice of $\epsilon$ in W-DHRAS model and the corresponding effect on the out-of-sample simulation performance. In Section~\ref{subsection:AppointmentTime}, we analyze the appointment time structure. In Section \ref{subsection:Out-of-sample}, we analyze the optimal solutions of the models and then compare their out-of-sample simulation performance. We also discuss the reliability of the models in Section \ref{subsection:Reliability}. In Section \ref{subsection:CPU_Time}, we compare the computational performance of the three models.

\subsection{\textbf{Description of the experiments}} \label{subsection:Expt_Setup}
We construct HRAS instances based on benchmarks, parameters settings and assumptions made in recent related literature (see, e.g., \citealp{jiang2019data}, \citealp{zhan2020home}, \citealp{zhan2018vehicle}). The average number of customers that an operator may visit per day is often less than six in home health care and banking, and less than 10 in repair service \citep{NAH2010, yuan2015branch,zhan2020home}. Accordingly, we consider problem instances with 6, 8 and 10 customers. For example, it is often not feasible for a caregiver to visit more than 10 patients (often 6)  a day considering care duration at each location and travel time between patients in home health care applications. However, we also test the computational performance of the proposed models under unrealistic larger instances of HRAS. We consider two different cost structures for waiting, idling and overtime in the objective function:  (a) $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ (\citealp{jiang2017integer}, \citealp{jiang2019data}), and (b) $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$ (\citealp{shehadeh2020distributionally}). For the transportation cost, we consider $\lambda\in\{0.5,1,2\}$ (\citealp{zhan2020home}).  We set  $L$ to 8 hours as in \cite{zhan2018vehicle}.

We use the lognormal distribution (\citealp{jiang2019data}) for the service time $d_i$ truncated on the interval $[10,50]$ with mean $\mu$ and $\sigma=0.5\mu$, where $\mu$ is generated from $U[25,35]$ ($U[a,b]$ refers to uniform distribution over the interval $[a,b]$). Our model works with any choices of the range $[10,50]$. Previous studies such as \cite{zhan2020home} assumes deterministic travel time. In our study, we generate the random travel time $t_{i,i'}$ from $U[15,25]$. That is, we assume that customers are fairly separated and traveling from one place to the other takes $20$ minutes on average (though our models can solve instances with any ranges and distribution of travel time). This is also consistent with prior and recent literature. \cite{nikzad2021matheuristic} particularly point out that customers within a service region form a basic unit or a cluster that share the same distribution of travel time, which is seen in urban areas. We round each generated parameter to the nearest integer.

We use the same upper and lower bounds of service time $[\underline{d}_i,\overline{d}_i]=[10,50]$ and travel time $[\underline{t}_{i,i'},\overline{t}_{i,i'}]=[15,25]$ in M-DHRAS and W-DHRAS. In the M-DHRAS model, the mean parameters $\mub^{\mbox{\tiny d}}$ and $\mub^{\mbox{\tiny t}}$ in the ambiguity set are set as the sample mean of the data. We also introduce new symmetry-breaking constraints (see \ref{apdx: SymBreak}) to break the symmetry in the routing decision. We implemented the three proposed models in AMPL programming language and use CPLEX (version 12.7.0.0) solver with the default setting. The relative MIP gap tolerance is set to $0.02$ while most of the instances have a terminal relative MIP gap tolerance very close to $0$. All the experiments were conducted on a computer with AMD Opteron 2.0 GHz CPU and 16 Gb memory. The time limit for solving each instance is set to 2 hours.

\subsection{\textbf{Effect of $\epsilon$ in W-DHRAS model}} \label{subsection:Wass_Eps}

In the W-DHRAS model, there is one parameter in the uncertainty set that serves as an input: the Wasserstein ball's radius $\epsilon$. In this section, we demonstrate the effect of $\epsilon$ on the out-of-sample performance of the W-DHRAS's optimal solution, $(\hat{\pmb{x}} (\epsilon, R),\hat{\pmb{a}} (\epsilon, R))$, with respect to the radius $\epsilon$. For illustrative purposes, we focus on one instance of $6$ customers with cost structure $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ and $\lambda=0.5$.

For each $\epsilon \in \{0.01, 0.02, \dots, 0.09, 0.1, \dots, 0.9,1,\dots,10\}$ (i.e., log-scaled interval as in \citealp{esfahani2018data} and \citealp{jiang2019data}), we evaluate the out-of-sample performance  as follows. First, for each customer, we randomly generate $30$ data sets of service duration and travel time scenarios, each consisting of $R\in \{ 5,20,50  \}$ scenarios. We generate these data sets under the same settings described in Section \ref{subsection:Expt_Setup}. Second, we solve the W-DHRAS model in \eqref{W_DHRAS_L} using the generated data sets under each of the candidate Wasserstein radius $\epsilon$. Finally, we fix the first-stage variables to the optimal solution of each instance, and then re-optimize the second-stage of the SP using 10,000  out-of-sample (unseen) data.  This is to compute the corresponding out-of-sample overtime, idle time, travel time, and waiting time and hence, the second-stage cost.
\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
       \includegraphics[scale=0.48]{CV_1.png}
       \caption{$R=5$}
       \label{fig:Wass_eps_5}
    \end{subfigure}
    %
    \begin{subfigure}[b]{0.3\textwidth}
       \includegraphics[scale=0.48]{CV_2.png}
       \caption{$R=20$}
       \label{fig:Wass_eps_20}
    \end{subfigure}
    %
    \begin{subfigure}[b]{0.3\textwidth}
       \includegraphics[scale=0.48]{CV_3.png}
       \caption{$R=50$}
       \label{fig:Wass_eps_50}
    \end{subfigure}
    %
    \caption{Performance of W-DHRAS with different choices of $\epsilon$}
    \label{fig:Wass_eps}
\end{figure}

Figure \ref{fig:Wass_eps} illustrates the out-of-sample performance with $R\in\{5,20,50\}$ under different choices of $\epsilon$. The blue line represents the mean of the 30 out-of-sample costs, while the shaded region is the area between the 20th and 80th percentiles of the 30 out-of-sample costs. It is quite evident that the out-of-sample performance first decreases with $\epsilon$ and then increases after some value of $\epsilon$. This pattern is often observed in the literature (see, e.g., \citealp{esfahani2018data} and \citealp{jiang2019data}). The values of $\epsilon$ which gives the smallest out-of-sample cost with $R=5$, $20$ and $50$ are respectively $5$ (Figure \ref{fig:Wass_eps_5}), $0.6$ (Figure \ref{fig:Wass_eps_20}) and $0.5$ (Figure \ref{fig:Wass_eps_50}) respectively. This decrease in $\epsilon$ with the increase in $R$ is not surprising. Intuitively, a small sample provides little information on the true distribution, and thus a larger $\epsilon$ produces robust solutions that better hedge against ambiguity. In contrast, with a larger sample, we have more information from the data, and so we can make less conservative decisions using a smaller $\epsilon$ value. As such, one should choose a larger (smaller) $\epsilon$ with a small (large) sample. Indeed, we can see that the out-of-sample cost for $R=20$ is smaller than that of $R=5$ in most of the cases, which indicates that a larger sample size could give a better out-of-sample performance. Moreover, this shows that when sufficient data is available, one could employ the SP model directly (with $\epsilon=0$). Thus, as mentioned earlier, W-DHRAS is particularly useful when there is a small data on random parameters.

In practice, decision-makers do not often have optimization expertise (or time) to run the above procedure (or other iterative or cross validation procedures). Additionally, we do not have full distributional information most of the time but only a small set of data samples. In the following experiments, we pick three different values of $\epsilon$, namely $0.5$, $5$, and $50$, which captures different extents of robustness of the W-DHRAS model. For brevity, we label them as W-DHRAS(0.5), W-DHRAS(5) and W-DHRAS(50) respectively.

\subsection{\textbf{Appointment time structure}}\label{subsection:AppointmentTime}
In this section, we analyze the optimal appointment structure produced by the SP, M-DHRAS, and W-DHRAS models.  For illustrative purposes and brevity, we present results under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ and $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$ with $\lambda=2$. We observe similar results with other choices of $\lambda$ (see \ref{apdx:AT_Structure}). Figures~\ref{fig:AppointmentTime1_lam2} and \ref{fig:AppointmentTime2_lam2} present the optimal schedules of the operator produced by the SP and DRO models under  $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ and $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$, respectively. The point $(x,y)=(i$,  inter-arrival time$)$ of every schedule in each subfigure corresponds to the mean optimal inter-arrival time (i.e., differences between the scheduled service start times of two consecutive customers, $I_j=a_j-a_{j-1}$ for $i\in[N]$ with $a_0=0$).
\begin{figure}[t]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{AT_Cost1_R50_lambda2.png}
    \caption{Mean inter-arrival times $a_j-a_{j-1}$ with $R=50$ under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ and $\lambda=2$}
    \label{fig:AppointmentTime1_lam2}
\end{figure}
\begin{figure}[t!]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{AT_Cost2_R50_lambda2.png}
    \caption{Mean inter-arrival times $a_j-a_{j-1}$ with $R=50$ under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$ and $\lambda=2$}
    \label{fig:AppointmentTime2_lam2}
\end{figure}

We first analyze the results under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ in Figure~\ref{fig:AppointmentTime1_lam2}. First, all models assign less time for the first customer than the subsequent customers and more time between the first and second customers. The models distribute (roughly) equal time between the middle customers (e.g., customers 2-5 when $N=6$ and customers 3-7 when $N=8$), and schedule less time between the last 2 to 3 customers, especially when $N=10$ for M-DHRAS, though SP allocates (roughly) equal time for both middle and last customers. Second, it is evident that M-DHRAS is the most conservative model and tends to schedule more time between customers when $N=6$ than the SP and W-DHRAS models. Also, M-DHRAS model yields a dome-shaped inter-arrival time structure (i.e., the time between customers first increases then decreases) when $N=8$ and $N=10$. Third, the SP model always schedules less time between the customers than the DRO models (except possibly the last few customers). This is not surprising since the SP model seeks risk-neutral decisions. Fourth, the W-DHRAS model yields a similar pattern as the SP model; however, it schedules more time between customers than the SP as $\epsilon$ increases. This makes sense because as $\epsilon$ increases, the W-DHRAS model becomes more conservative and schedules more time between appointments. This demonstrates that the W-DHRAS model is more conservative in hedging ambiguity than the SP model but less conservative than the M-DHRAS model. In the next section, we show that the W-DHRAS appointment pattern yields the best out-of-sample performance under various distributions.

Next, we analyze the results under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$ in Figure \ref{fig:AppointmentTime2_lam2}. We observe different patterns under this cost structure. All models assign more time for the first few customers but gradually less time for the middle and last few customers. The pattern is different from the cost structure, mainly due to a significantly higher idling cost relative to the waiting cost. However, we can still conclude that M-DHRAS is the most conservative model, but instead of allocating more time, it allocates less time between customers to hedge against scenarios that yields long idling time. Third, we observe that the dome-shaped inter-arrival time structure for all models, which is particularly prominent in $N=6$. Forth, we note that W-DHRAS models with $\epsilon=0.5$ and $\epsilon=5$ roughly allocates the same amount of time between customers as the SP model. However, with $\epsilon=50$, W-DHRAS model tends to allocate more time for the first few customers but fewer time for the last few customers.

\subsection{\textbf{Out-of-sample performance}}\label{subsection:Out-of-sample}
In this section, we compare the optimal solutions of the three models we analyzed in the previous section under limited historical data.  We test the out-of-sample performance (i.e.,  the objective value obtained by simulating the optimal solution of a model under a larger unseen data) as follows. For each $N\in\{6,8,10\}$ and $R \in \{5, 10, 20, 50\}$, we generate $30$ SP, M-DHRAS and W-DHRAS instances using the same parameters settings described in Section \ref{subsection:Expt_Setup}. Note that the sample sizes are not large to mimic the situation where the decision-making possesses limited data. We solve each instance and obtain the optimal first-stage $(\pmb{x},\pmb{a})$ decision. Fixing the first-stage decision to $(\pmb{x},\pmb{a})$, we then re-optimize the second-stage of the SP using the following five sets of  $10,000$ samples.  In the first set, we assume that the data we rely on in the optimization comes from the true distribution. We use data from Sets 2-5 to evaluate the performance of the models under the case when the data we rely on in the optimization may follow a biased or misspecified distribution different from the true distribution.

\begin{enumerate}\itemsep0em
    \item [Set 1.] We assume perfect information for the distributions. That is, we generate the $10,000$ samples from the same distribution we use in the optimization as discussed in Section \ref{subsection:Expt_Setup}. This simulation assumes that the data comes from the true unknown distribution.

    \item [Set 2.] In this set, we assume that we have misspecified the distribution of the travel time in the optimization. We generate $t_{i,i'}$ from $U[25,35]$ instead of $U[15,25]$. That is, the average travel time takes 10 minutes longer than usual. This situation (shift in the travel time range) might be seen in practice due to unexpected traffic congestion (e.g., caused by traffic accidents, weather conditions etc.).

    \item [Set 3.]In this set, we assume that we have misspecified both the service and travel times distributions in the optimization. Specifically, we follow a similar out-of-sample simulation testing procedure described in \cite{wang2020distributionally} and perturb the support of the random travel and service times by a parameter $\delta$ as $[(1-\delta)$ lower bound, $(1+\delta)$ upper bound], where $\delta\in\{0.1,0.25,0.5\}$. A higher value of $\delta$ corresponds to a higher variation level.
\end{enumerate}

For brevity, we focus on HRAS instances of $N = 6$ customers (i.e., the average number of customers that a home service operator often visits per day; see Section \ref{subsection:Expt_Setup}) and present results under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ and $\lambda=2$. We observe similar out-of-sample performance with $N=8$ and $N=10$, and under other values of $\lambda$ as well as two additional sets of misspecified distribution (see these results in  \ref{apdx:OutSample}). Figure \ref{fig:OC_1} shows the results using Set 1, where the shaded region represents the 20\% and 80\% percentiles of the out-of-sample cost. Obviously, the optimal solutions of the M-DHRAS model have poor out-of-sample performance. Recall from Section \ref{subsection:AppointmentTime} that M-DHRAS schedules a long time between customers, which, as shown in Figure \ref{fig:IWO_1}, yields a significant amount of operator idle time (indicating poor utilization of the operator's time). The out-of-sample performance of the W-DHRAS depends on $\epsilon$, which intuitively follows from the observation that the optimal solutions depend on $\epsilon$ (Section \ref{WDHRAS}). When $\epsilon=0.5$, the performance of W-DHRAS is similar to SP, while a more conservative choice of $\epsilon$ such as $50$ yields a poorer performance as compared to SP. When $\epsilon=5$,  W-DHRAS  outperforms SP, especially when the sample size is small. This demonstrates that the W-DHRAS model has a superior performance when the data size is small (i.e., when there is limited data on random parameters).
\begin{figure}[t]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{OC_DataI_Cost1_lambda2_N6.png}
    \caption{Out-of-sample cost with cost structure (a) and $\lambda=2$ under Set 1}
    \label{fig:OC_1}
\end{figure}
\begin{figure}[t]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{IW_DataI_Cost1_lambda2_N6.png}
    \caption{Mean idle time, wait time and overtime with cost structure (a) and $\lambda=2$ under Set 1}
    \label{fig:IWO_1}
\end{figure}

Next, we discuss results for the case when we have misspecified the distributions. Figure \ref{fig:OC_2} shows the performance under Set 2, where a longer travel time is encountered. In contrast to the results in Figure \ref{fig:OC_1}, M-DHRAS and W-DHARS(50) outperforms SP, W-DHRAS(0.5) and W-DHRAS(5). This demonstrates that M-DHRAS and W-DHRAS(50), which are more conservative than the other models, have superior performance in environments where travel time increases by a lot.  This is not suprising since M-DHRAS and W-DHRAS(50) schedule longer time between customers that could hedge against the increase in travel time. The W-DHRAS model outperforms the SP model and the M-DHRAS model when $\epsilon=50 $ (a conservative choice of $\epsilon$).  We can observe from Figure \ref{fig:IWO_2} the trade-off between idle and wait time. Under Set 2, SP model results in a significant amount of waiting time but a small idling time compared with the two DRO models.
\begin{figure}[t]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{OC_DataII_Cost1_lambda2_N6.png}
    \caption{Out-of-sample cost with cost structure (a) and $\lambda=2$ under Set 2}
    \label{fig:OC_2}
\end{figure}
\begin{figure}
    \hspace{-10mm}
    \includegraphics[scale=0.60]{IW_DataII_Cost1_lambda2_N6.png}
    \caption{Mean idle time, wait time and overtime with cost structure (a) and $\lambda=2$ under Set 2}
    \label{fig:IWO_2}
\end{figure}

We also demonstrate one of the cases that the supports are perturbed. Figures \ref{fig:OC_3a} to \ref{fig:OC_3c} show the results for Set 3 with different choices of $\delta$. We observe that M-DHRAS performs better (with respect to SP) when $\delta$ increases, which corresponds to a larger deviation. The W-DHRAS(5) model consistently outperforms the SP model. We also observe that the W-DHRAS model has a larger edge over SP when the sample size is small. This is consistent with what we have discussed in previous cases. These results demonstrate that the DRO models perform better than the SP model when there is a large deviation from the sample distribution.

\begin{figure}[t]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{OC_DataIVa_Cost1_lambda2_N6.png}
    \caption{Out-of-sample cost with cost structure (a) and $\lambda=2$ under Set 3 ($\delta=0.10$)}
    \label{fig:OC_3a}
\end{figure}
\begin{figure}[t]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{OC_DataIVb_Cost1_lambda2_N6.png}
    \caption{Out-of-sample cost with cost structure (a) and $\lambda=2$ under Set 3 ($\delta=0.25$)}
    \label{fig:OC_3b}
\end{figure}
\begin{figure}[t!]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{OC_DataIVc_Cost1_lambda2_N6.png}
    \caption{Out-of-sample cost with cost structure (a) and $\lambda=2$ under Set 3 ($\delta=0.50$)}
    \label{fig:OC_3c}
\end{figure}

\subsection{\textbf{Reliability of the models}}\label{subsection:Reliability}
In this section, we analyze the reliability of the three models. Home service operators often need to estimate the budgets on the operational costs ahead of the actual service provision. The optimal value of an HRAS optimization model serves as a relevant estimate of the associated cost if we implement the corresponding optimal decision $(\pmb{x},\pmb{a})$. Intuitively,  when the estimated cost is smaller than the actual cost, the operator may run into financial problems related to a budget deficit and poor planning decisions. A risk-averse operator would then prefer implementing decisions that provide an upper bound on the estimated cost, i.e., seek solutions with higher reliability. Hence, we make use of the reliability measure. For any given decision $(\pmb{x},\pmb{a})$, reliability is the probability that the optimal value from the model is greater than or equal to the actual cost $\mathbb{E}[f(\pmb{x},\pmb{a},\pmb{\xi})]$. Mathematically, if $v(R)$ is the optimal value and $(\pmb{x}(R),\pmb{a}(R))$ is the corresponding optimal solution with $R$ samples, $$\text{reliability}=\mathbb{P}_{\pmb{\xi}}^R\Big(v(R) \geq \mathbb{E}_{\mathbb{P}_{\pmb{\xi}}}[f(\pmb{x}(R),\pmb{a}(R),\pmb{\xi})]\Big),$$ where $\mathbb{P}_{\pmb{\xi}}^R$ is the product measure of $R$ copies of $\mathbb{P}_{\pmb{\xi}}$. The actual cost $\mathbb{E}_{\mathbb{P}_{\pmb{\xi}}}[f(\pmb{x}(R),\pmb{a}(R),\pmb{\xi})]$ is estimated from $10,000$ out-of-sample scenarios and the reliability is computed from $30$ instances.

Figures \ref{fig:Reliability1} and \ref{fig:Reliability2} show the results for the two cost structures $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ and $(1,5,7.5)$ respectively with $\lambda=2$. We observe that the reliability is, in general, increasing with $R$. This is reasonable since we have a larger amount of historical data. The reliability of SP is the lowest among the three models, which demonstrates that SP with a small number of historical scenarios may not be able to provide a robust cost estimate. On the other hand, the reliability of W-DHRAS(50) consistently gives the best reliability result. This makes sense since $\epsilon=50$ corresponds to a relatively robust model. It is interesting that although M-DHRAS appears to be the most conservative model from the out-of-sample cost perspective, the reliability result is not the best. Overall, this demonstrates that DRO models have a better reliability compared to SP.
\begin{figure}[t]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{Re_Cost1_lambda2.png}
    \caption{Reliability under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$ and $\lambda=2$}
    \label{fig:Reliability1}
\end{figure}
\begin{figure}[t!]
    \hspace{-10mm}
    \includegraphics[scale=0.60]{Re_Cost2_lambda2.png}
    \caption{Reliability under $(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$ and $\lambda=2$}
    \label{fig:Reliability2}
\end{figure}

\subsection{\textbf{CPU time}}\label{subsection:CPU_Time}

In this section, we analyze the computational time for solving the constructed practical instances using the three models. Specifically, we examine the effect of $N$ (number of customers) and $R$ (number of scenarios) on the solution time for the three models as follows. Table $\ref{table:CPU_Time}$ shows the mean CPU time for solving $30$ instances using two different cost structures with $\lambda=2$.

\begin{table}[t!]\centering\small
\renewcommand{\arraystretch}{0.4}
\caption{CPU time in seconds for solving the three models are shown and the number in parentheses followed by W-DHRAS model is the choice of $\epsilon$. The reported times are for $\lambda=2$ with two different cost structures. }
\begin{tabular}{@{}lrrrr|rrrr@{}} \toprule
$N = 6$ & \multicolumn{4}{c}{$(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$} & \multicolumn{4}{c}{$(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$}  \\
\cmidrule{2-5}  \cmidrule{6-9}
CPU Time (in s) & $R=5$ & $R=10$ & $R=20$ & $R=50$ & $R=5$ & $R=10$ & $R=20$ & $R=50$\\
\cmidrule{1-5} \cmidrule{6-9}
SP           & 0.25 & 0.30 & 0.37  & 0.42   & 0.23 & 0.21 & 0.33  & 0.44   \\
M-DHRAS      & 0.46 & 0.48 & 0.47  & 0.48   & 0.46 & 0.46 & 0.46  & 0.46   \\
W-DHRAS(0.5) & 2.21 & 5.24 & 18.57 & 102.52 & 1.95 & 5.04 & 18.94 & 108.99 \\
W-DHRAS(5)   & 2.33 & 6.02 & 20.93 & 125.19 & 1.99 & 5.22 & 20.26 & 121.32 \\
W-DHRAS(50)  & 2.34 & 5.34 & 22.93 & 134.77 & 2.48 & 7.00 & 25.11 & 170.84 \\
\toprule
$N = 8$ & \multicolumn{4}{c}{$(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$} & \multicolumn{4}{c}{$(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$}  \\
\cmidrule{2-5}  \cmidrule{6-9}
CPU Time (in s) & $R=5$ & $R=10$ & $R=20$ & $R=50$ & $R=5$ & $R=10$ & $R=20$ & $R=50$\\
\cmidrule{1-5} \cmidrule{6-9}
SP           & 0.28  & 0.33  & 0.39   & 0.76    & 0.28 & 0.32  & 0.39   & 0.79   \\
M-DHRAS      & 1.11  & 1.05  & 1.07   & 1.10    & 0.90 & 0.92  & 0.91   & 0.93   \\
W-DHRAS(0.5) & 6.14  & 16.70 & 90.07  & 583.90  & 5.29 & 15.29 & 81.90  & 514.50 \\
W-DHRAS(5)   & 7.42  & 20.42 & 101.00 & 681.52  & 5.59 & 18.05 & 85.59  & 530.59 \\
W-DHRAS(50)  & 10.05 & 27.10 & 145.41 & 1295.76 & 8.98 & 24.90 & 114.88 & 718.17 \\
\toprule
$N = 10$ & \multicolumn{4}{c}{$(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(2,1,20)$} & \multicolumn{4}{c}{$(c^{\mbox{\tiny w}}_j,c^{\mbox{\tiny u}}_j,c^{\mbox{\tiny o}})=(1,5,7.5)$}  \\
\cmidrule{2-5}  \cmidrule{6-9}
CPU Time (in s) & $R=5$ & $R=10$ & $R=20$ & $R=50$ & $R=5$ & $R=10$ & $R=20$ & $R=50$\\
\cmidrule{1-5} \cmidrule{6-9}
SP           & 0.376 & 0.43  & 0.54   & 1.04    & 0.35  & 0.42  & 0.54   & 1.07    \\
M-DHRAS      & 2.05  & 2.02  & 2.11   & 1.97    & 1.79  & 1.83  & 1.86   & 1.89    \\
W-DHRAS(0.5) & 17.84 & 56.71 & 249.84 & 1801.04 & 15.07 & 51.71 & 220.39 & 1415.10 \\
W-DHRAS(5)   & 18.22 & 61.54 & 260.82 & 1922.29 & 15.54 & 53.34 & 228.30 & 1523.63 \\
W-DHRAS(50)  & 19.98 & 70.72 & 297.87 & 2490.39 & 16.41 & 59.70 & 270.15 & 1919.69 \\
\bottomrule
\end{tabular}
\label{table:CPU_Time}
\end{table}

We first observe that the SP model takes the shortest time to solve for all instances and most of them can be solved within one second. Moreover, SP requires a longer solution time with a larger number of scenarios. The M-DHRAS model has a slightly longer solution time than the SP model but, since it only depends on the mean and range of the sample, it has a consistent performance under all values of $R$. Solution times of the W-DHRAS model are longer than those of the M-DHRAS and SP models. Solution times of the W-DHRAS model increase with $R$ and varies across three choices of $\epsilon$. We attribute the difference in solution times between the SP and W-DHRAS models to their respective sizes (i.e., the number of variables and constraints).  For fixed $R$, the W-DHRAS model has more variables and constraints than the SP model. As pointed out by previous papers (e.g. \citealp{artigues2015mixed,klotz2013practical,shehadeh2020distributionally}), an increase in model size often suggests an increase in solution time for the LP relaxation and, thus, the MILP model’s overall solution time.

Recall that W-DHRAS is theoretically more useful when there is only a small data on random parameters. When there is a large data set, W-DHRAS may converge to SP with $\epsilon$ close to zero and our SP is efficient with large $R$ (see \ref{apdx:SolTime}). Thus, although W-DHRAS takes the longest time among the three, it can solve realistic instances within a reasonable time, especially when $R$ is small. Specifically, when $N=6$ (common in home healthcare applications), W-DHRAS solution times range from 2 seconds ($R=5$) to 3 minutes ($R=50$). When $N=10$ (the maximum number of customers in HRAS applications; see discussion and references in Section~\ref{subsection:Expt_Setup}), W-DHRAS solution times range from 15 seconds ($R=5$) to around 30 minutes ($R=50$ and $\epsilon=0.5, 5$).

In \ref{apdx:SolTime}, we demonstrate the SP and M-DHRAS abilities to solve (unrealistic) large instances with a large number of customers and scenarios. For example, the solutions times of the SP with $N=15$ customers and $R=500$ scenarios are less than 1.1 minutes. Solution times of the M-DHRAS for the same instance are negligible.

\section{Conclusion} \label{sec5:conclusion}

In this paper, we address the uncertainty and distributional ambiguity of service and travel times in an HRAS problem. We propose and analyze the first SP model and two DRO models for HRAS with random service and travel times. We use two popular ambiguity sets in the DRO models: the mean-support ambiguity set (in the M-DHRAS model) and the $1$-Wasserstein ambiguity set (in the W-DHRAS model).  We derive equivalent mixed-integer linear programming (MILP) reformulations of both DRO models that can be implemented and efficiently solved using off-the-shelf optimization software.  We conduct extensive experiments comparing the proposed models. Our results demonstrate that: (1) W-DHRAS yields robust decisions with both small and large data size, and enjoys both asymptotic consistency and finite-data guarantees; (2) M-DHRAS produces the most conservative schedules and it performs well when the distribution of travel time changes dramatically (e.g., actual travel times are longer); (3) W-DHRAS solutions have better operational performance than the SP solutions both under perfect information and misspecified distributions, even when only a small data set is available; (4) DRO models produce more reliable solutions than the SP model; (5) the proposed SP and DRO models are computationally efficient under realistic HRAS settings. Thus, we conclude that the proposed DRO models are useful if the distributions of the random service and travel time are unknown and only a limited amount of data is available. The resulting optimal solutions are robust to misspecifications of the underlying distribution with high reliability. This addresses the primary goal of our paper to investigate the value of the DRO approach for HRAS and compare it with the SP approach.

Under the circumstance that only limited data is available, decision making process under uncertainty (in both service and travel times) is typically challenging. This is particularly true in the home health care context where the service time data is often not available. Except in the case that the underlying distribution coincides with the sample data distribution, the standard SP approach appears to be vulnerable (e.g., poor performance and reliability). The use of DRO models could remedy the distribution misspecification problem and provide a more robust solution, leading to a reduction in out-of-sample costs and improvement in reliability. This is typically true when we experience changes in the underlying distribution. While M-DHRAS yields the most conservative solution (e.g., large inter-arrival times), the model yields superior performance when there is a large perturbation in distribution (e.g., drastic increase in travel time). For instance, when the operator is planning for a day with possibly severe traffic congestion or weather condition (and hence, travel time is highly random), adopting M-DHRAS solutions will yield better operational performance. On the other hand, under normal operations, adopting W-DHRAS solutions may yield better operational performance  (e.g., with an empirical choice of $\epsilon=5$) and reduce out-of-sample costs due to distribution misspecification with a better reliability than SP.

 Our model can serve as a building block for the following future extensions and areas of research in terms of home service aspects, constraints, and various sources of uncertainties. First, we want to generalize our distributionally robust approach to the case when we have multiple service providers and include decisions such as (1) determining the number of service teams to hire,  (2) assigning the hired service teams to customers, (3) constructing routes for the service teams, and (4) determining the customers' appointment times. Second, we aim to model other sources of uncertainty, such as last-minute customer cancelations and operator's cancelation on the day of service. Third, we also aim to incorporate customers' preferences on appointment times.

\bibliographystyle{elsarticle-harv}
\linespread{1}