{'fraction_non_alphanumeric': 0.048105525675310234, 'fraction_numerical': 0.013716534959887023, 'mean_word_length': 4.170343327637098, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Optimization Models for a Home Service Routing and Appointment Scheduling Problem with Random Travel and Service Times

We study a routing and appointment scheduling problem with uncertain service and travel times arising from home service practice. Specifically, given a set of customers within a service region that an operator needs to serve, we seek to find the operator's route and time schedule. 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. We propose and rigorously analyze a stochastic programming model and two distributionally robust optimization (DRO) models to solve the problem, assuming a known and unknown service and travel time distribution, respectively. We consider two popular types of ambiguity sets for the DRO models, namely, the mean-support and 1-Wasserstein ambiguity set. 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, thereby enabling practitioners to use these models. In an extensive numerical experiment, we investigate the proposed models computational and operational performance, demonstrating where significant performance improvements can be gained with each proposed model and derive insights into the problem.


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 (Lunden, 2019;Zhan et al., 2021). 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 "routing" to refer to routes that specify the operator movement between customers or the sequence of visits on the operator's schedule, and "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 Cook, 2011). 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 (Ahmadi-Javid et al., 2017;Robinson and Chen, 2003). 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 (Grieco et al., 2020;Kandakoglu et al., 2020).
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 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 (Kuhn et al., 2019). Therefore, DRO solutions often faithfully anticipate the possibility of disappointing out-ofsample consequences and depending on the ambiguity set, they often guarantee an out-of-sample cost that falls below the worst-case optimal cost (Kuhn et al., 2019;Rahimi and Gandomi, 2020).
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, Jiang et al., 2017a;Saif and Delage, 2021;Shehadeh et al., 2020;Shehadeh and Tucker, 2020;Wang et al., 2019Wang et al., , 2020Zhang et al., 2017).
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 φ-divergence (Jiang and Guan, 2016) and Wasserstein distance (Mohajerin Esfahani and Kuhn, 2018;Gao and Kleywegt, 2016), 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 (Kuhn et al., 2019;Mohajerin Esfahani and Kuhn, 2018). 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: 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 meansupport ambiguity and Wasserstein ambiguity sets and demonstrate where significant performance improvements can be gained with each proposed model and derive insights into the problem.

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 distributionally robust home service routing and scheduling (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 minmax 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 (Grieco et al., 2020) and our literature review in Section 2, 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 2, we review the relevant literature. In Section 3, we formally define the DHRAS models, present and analyze the M-DHRAS and W-DHRAS models. In Section 4, we present computational results and managerial implications. Finally, we draw conclusions and discuss future directions in Section 5.

Relevant 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 Fikar and Hirsch (2017) 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., Berg et al., 2014;Shehadeh et al., 2019). Appointment times can be fixed time points (see, e.g., Braekers et al., 2016;Demirbilek et al., 2019;Milburn, 2012) or time windows (see, e.g., Lee et al., 2013;Mankowska et al., 2014;Yuan et al., 2015 and references therein). Herein, we focus on the former case under random service and travel times. We refer to Pinedo (2016) for a comprehensive survey of scheduling theory, applications and methods. consider a single operator HRAS problem in which they assume a deterministic travel time between customers and that each customer's service time follows a fully known probability distribution.
Accordingly, Zhan et al. (2021) 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) ap-proach assumes that uncertain parameters resided on an uncertainty set of possible outcomes with some structure (e.g., polyhedron, see, Ben-Tal et al., 2015;Bertsimas and Sim, 2004;Soyster, 1973). Accordingly, optimization in RO is based on the worst-case scenario within the uncertainty set.
Recently, Shi et al. (2019) 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 Shi et al. (2018), 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 (Chen et al., 2020;Delage and Saif, 2018). 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 (Rahimian and Mehrotra, 2019).
There are two common types of ambiguity sets: moment-based and distance-based ambiguity sets (see, e.g., Rahimian and Mehrotra, 2019). 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 Fikar and Hirsch (2017); Grieco et al.
(2020), there is no DRO approach for the specific HRAS problem that we address. Moment-based DRO models for single-server appointment scheduling includes (Jiang et al., 2017b;Kong et al., 2013;Mak et al., 2015). Recently, Jiang et al. (2019) 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 Jiang et al. (2019) and Zhan et al. (2021), which are recent singleserver studies relevant to our work (see Table A.3 in Appendix A for a quick comparison). First, the three papers address stochastic service time. Second, Jiang et al. (2019) 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, Jiang et al.
(2019)'s model cannot be used for routing and appointment scheduling. Third, Zhan et al. (2021) assume that the travel time is deterministic, which as mentioned earlier, may lead to sub-optimal

Formulation and Analysis
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 3.1. We analyze and derive equivalent mixed-integer linear programming reformulations of M-DHRAS in Section 3.2 and W-DHRAS in Section 3.3.

Definitions and assumptions
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, . . . , N ) is random, with i = 0 representing the service provider's office (depot). Service time d i of each customer i ∈ [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 ∈ [N ] and j ∈ [N ], we let the binary decision variable x i,j equal 1 if customer i is the jth customer in the operator's route/schedule, and zero otherwise. For all j ∈ [N ], we let the continuous variable a j represent the appointment starting time of the jth customer. The feasible region of variable x x x is defined in (1) 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 a a a is defined in (2) such that all appointments are scheduled within the provider's service hours [0, L].
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 f (x x x, a a a, ξ ξ ξ) := min where c w j , c u j , c o , and λ are respectively the non-negative unit penalty costs of waiting, idling, overtime, and travel time for all j ∈ [N ]. Also, we let a N +1 = L and x i,N +1 = 0 for all i ∈ [N ]. The objective function (3a) minimizes a linear cost function of waiting, idling, overtime, and travel time.
Constraint (3b) yields either provider's idle time before the scheduled time of the first appointment or the waiting time of the first customer. Constraint (3c) yields either the waiting time of the jth customer or the provider's idle time if s/he arrives at the jth customer, respectively, after or before the scheduled start time a j of the jth customer. Constraints (3d) yields either the overtime or the schedule earliness. Constraint (3e) computes the operator's total travel time.
The SP model in (4) seeks to find (x x x, a a a) ∈ X × A that minimizes the expectation of the cost f (x x x, a a a, ξ ξ ξ) subject to random ξ ξ ξ with a known joint probability distribution P ξ ξ ξ .

DHRAS over moment ambiguity (M-DHRAS)
In this section, we consider the case when P ξ ξ ξ is not fully known or hard to estimate. We assume that we know or can approximate the mean values E P ξ ξ ξ (ξ ξ ξ) 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 µ µ µ and S, we consider the following mean-support ambiguity set F(S, µ µ µ): where P(S) represents the set of all probability distributions supported on S. Using ambiguity set E P [f (x x x, a a a, ξ ξ ξ)] .
Formulation (7) seeks first stage decisions (x x x, a a a) that minimizes the worst-case (maximum) expectations of the second-stage operational cost over F(S, µ µ µ).
Remark 1. Model (7) 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.

MILP reformulation of M-DHRAS
Note that model (7) 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 (7). First, for a fixed (x x x, a a a) ∈ X × A, we rewrite the inner maximization problem sup P∈F (S,µ µ µ) E P [f (x x x, a a a, ξ ξ ξ)] as follows: where 1 S (ξ ξ ξ) = 1 if ξ ξ ξ ∈ S and 1 S (ξ ξ ξ) = 0 if ξ ξ ξ / ∈ S. As we show in the proof of Proposition 1 in Appendix B, problem (8) is equivalent to the deterministic problem (9).
Proposition 1. For any (x x x, a a a) ∈ X × A, problem (8) is equivalent to Note that f (x x x, a a a, ξ ξ ξ) is a minimization problem, and thus in (9) we have an inner max-min problem, which is not suitable to solve using standard solution methods. For a given solution (x x x, a a a) and realized value of ξ ξ ξ, f (x x x, a a a, ξ ξ ξ) is a linear program (LP). The dual of f (x x x, a a a, ξ ξ ξ) is as follow where y 1 , . . . , y N +1 are the dual variables associated with constraints (3b)-(3d). Given (x x x, a a a) ∈ X × A and ξ ξ ξ, the objective function in (10a) is linear (convex) in y y y. Hence, the inner maximization problem in (10) is a convex maximization problem. It follows from the fundamental convex analysis that it suffices to consider the extreme points of polytope Y. This motivates us to leverage the properties of the extreme points of Y in deriving an equivalent LP reformulation of (10). We formally prove this in Proposition 2 (see Appendix C for a detailed proof).
Proposition 2. For any (x x x, a a a) ∈ X × A, the inner maximization problem in (9) is equivalent to . where With the use of Proposition 2, we can combine the minimization over α α α and ρ ρ ρ in problem (9).
Corollary 3 provides the reformulation of problem (9) (see Appendix D for a proof).
Corollary 3. Problem (9) is equivalent to min α α α, ρ ρ ρ, β β β, γ γ γ, δ δ δ in the form of (12) with the outer minimization problem in (7), we obtain a mixed-integer nonlinear program (MINLP). Due to its large size and page limitation, we present this MINLP in Appendix E, where we also derive the following exact and tight MILP reformulation of this this MINLP (equivalently, M-DHRAS model).

DHRAS over 1-Wasserstein ambiguity (W-DHRAS)
In this section, we consider the case that P ξ ξ ξ may be observed via a small finite set ξ ξ ξ 1 , . . . ,ξ ξ ξ R 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 P ξ ξ ξ close to the empirical distribution P R ξ ξ ξ . Specifically, we construct a set of all distributions P(S) supported on the support S of the unknown true distribution for which all distributions P ξ ξ ξ have 1-Wasserstein distance that is less than or equal to from the reference distribution. We use 1-Wasserstein distance (i.e., we use 1 -norm in the definition of Wasserstein metric) because it often admits tractable reformulation in most real-world applications (see, e.g., Duque and Morton, 2020;Hanasusanto and Kuhn, 2018;Jiang et al., 2019;Saif and Delage, 2021). Next, we present our W-DHRAS model.
Suppose that two probability distributions Q 1 and Q 2 are defined on a common support S ∈ R N , and let || · || p represent the p-norm on R N with p ≥ 1. Suppose that random vectors ξ ξ ξ 1 and ξ ξ ξ 2 follow Q 1 and Q 2 respectively. Then, the Wasserstein distance between Q 1 and Q 2 , denoted as W p (Q 1 , Q 2 ), represents the minimum transportation cost of moving from Q 1 to Q 2 , where the cost of moving from ξ ξ ξ 1 to ξ ξ ξ 2 is measured by the norm ξ ξ ξ 1 − ξ ξ ξ 2 p . Mathematically, where P(Q 1 , Q 2 ) is the set of all joint distributions of (ξ ξ ξ 1 , ξ ξ ξ 2 ) with marginals Q 1 and Q 2 . Since we only observe a set {ξ ξ ξ 1 , . . . ,ξ ξ ξ R } of R i.i.d. samples, we consider the following p-Wasserstein ambiguity set where P(S) is the set of all probability distributions on S,P R ξ ξ ξ = 1 R R r=1 δξ ξ ξ r is the empirical distribution of ξ ξ ξ based on the R i.i.d samples with δ being the Dirac measure, and > 0 is the radius of the ambiguity set. The set F p (P R ξ ξ ξ , ) can be viewed as the p-Wasserstein ball of radius centered at the empirical distributionP R ξ ξ ξ . Therefore, in some sense, one can think Wasserstein ball as the set of all distributions under which our estimation error is below , where is the estimation error we seek protection against. A larger radius indicates that we seek more robust solutions.
Using the ambiguity set F p (P R ξ ξ ξ , ), we formulate W-DHRAS as follows: 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 (16)Ẑ(R, ) converges to Z * (the optimal value of the SP model in (4) with perfect knowledge of P ξ ξ ξ ), and an optimal solution (x x x, a a a) of (W-DHRAS) converges to an optimal solution of problem (4). Additionally, ifẐ(R, ) > 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 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 Appendix F-Appendix H for adapted proofs of these results.

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 (x x x, a a a) ∈ X × A.
In Proposition 4, we present an equivalent dual formulation of (17) (see Appendix I for a proof).
Proposition 4. The optimal value of formulation (17) equals that of the following formulation: Note that f (x x x, a a a, ξ ξ ξ) − ρ ξ ξ ξ −ξ ξ ξ r 1 is neither convex nor concave in ξ ξ ξ. Thus, formulation (18) 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 ρ and (x x x, a a a) ∈ X × 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 (x x x, a a a) ∈ X × A and ρ ≥ 0, we denote g r (ρ, x x x, a a a) = sup ξ∈S {f (x x x, a a a, ξ ξ ξ) − ρ ξ ξ ξ −ξ ξ ξ r 1 }. Given the dual formulation of f (x x x, a a a, ξ ξ ξ) in (10), and using the dual of f (x x x, a a a, ξ ξ ξ) and using the same proof techniques in Proposition 2, we can reformulate g r (ρ, x x x, a a a) (see Appendix J for a detailed proof).
Proposition 5. The maximization problem g r (ρ, x x x, a a a) is equivalent to Summing g r (ρ, x x x, a a a) in the form of Proposition 5 over r and combining it with the outer minimization in (18) and (16), we derive MINLP reformulation of the W-DHRAS. Again due to its large size, we present this MINLP in Appendix K, where we also derive the following tight and exact MILP reformulation of the W-DHRAS model.
Remark 2. 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 = ∞ since the term ξ ξ ξ 1 − ξ ξ ξ 2 ∞ 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.

Numerical Experiments
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 (4) with P ξ ξ ξ replaced by an empirical distribution based on N samples of the random parameters (see Appendix L 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 4.1, we describe the set of HRAS instances that we constructed and discuss other experimental setup. In Section 4.2, we examine the choice of in W-DHRAS model and the corresponding effect on the out-of-sample simulation performance. In Section 4.3, we analyze the appointment time structure. In Section 4.4, 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 4.5. In Section 4.6, we compare the computational performance of the three models.

Description of the experiments
We construct HRAS instances based on benchmarks, parameters settings and assumptions made 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.

Effect of in W-DHRAS model
In the W-DHRAS model, there is one parameter in the uncertainty set that serves as an input: the Wasserstein ball's radius . In this section, we demonstrate the effect of on the out-ofsample performance of the W-DHRAS's optimal solution, (x x x( , R),â a a( , R)), with respect to the radius . For illustrative purposes, we focus on one instance of 6 customers with cost structure , 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 ∈ {5, 20, 50} scenarios. We generate these data sets under the same settings described in Section 4.1. Second, we solve the W-DHRAS model in (20) using the generated data sets under each of the candidate Wasserstein radius . Finally, we fix the first-stage variables to the optimal solution of each instance, and then re-optimize the secondstage 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. in with the increase in R is not surprising. Intuitively, a small sample provides little information on the true distribution, and thus a larger 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 value. As such, one should choose a larger (smaller) 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 = 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 , 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.

Appointment time structure
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 w j , c u j , c o ) = (2, 1, 20) and (c w j , c u j , c o ) = (1, 5, 7.5) with λ = 2. We observe similar results with other choices of λ (see Appendix N). Figures 2 and 3 present the optimal schedules of the operator produced by .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 ∈ [N ] with a 0 = 0).
We first analyze the results under (c w j , c u j , c o ) = (2, 1, 20) in Figure 2. 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 increases. This makes sense because as 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 w j , c u j , c o ) = (1, 5, 7.5) in Figure 3. 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 = 0.5 and = 5 roughly allocates the same amount of time between customers as the SP model. However, with = 50, W-DHRAS model tends to allocate more time for the first few customers but fewer time for the last few customers.

Out-of-sample performance
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 ∈ {6, 8, 10} and R ∈ {5, 10, 20, 50}, we generate 30 SP, M-DHRAS and W-DHRAS instances using the same parameters settings described in Section 4.1. 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 (x x x, a a a) decision. Fixing the first-stage decision to (x x x, a a 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. 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 4.1. This simulation assumes that the data comes from the true unknown distribution.
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 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).
Next, we discuss results for the case when we have misspecified the distributions. Figure 6 shows the performance under Set 2, where a longer travel time is encountered. In contrast to the results We also demonstrate one of the cases that the supports are perturbed. Figures 8 to 10 show the results for Set 3 with different choices of δ. We observe that M-DHRAS performs better (with respect to SP) when δ 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.

Reliability of the models
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 (x x x, a a 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 (x x x, a a a), reliability is the probability that the optimal value from the model is greater than or equal to the actual cost E[f (x x x, a a a, ξ ξ ξ)].
Mathematically, if v(R) is the optimal value and (x x x(R), a a a(R)) is the corresponding optimal solution with R samples, where P R ξ ξ ξ is the product measure of R copies of P ξ ξ ξ . The actual cost E P ξ ξ ξ [f (x x x(R), a a a(R), ξ ξ ξ)] is estimated from 10, 000 out-of-sample scenarios and the reliability is computed from 30 instances. (1, 5, 7.5) respectively with λ = 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 = 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.

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 2 shows the mean CPU time for solving 30 instances using two different cost structures with λ = 2.
We first observe that the SP model takes the shortest time to solve for all instances and most 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 close to zero and our SP is efficient with large R (see Appendix P). Thus, although W-DHRAS takes the longest time among the three, it can solve realistic instances within a reasonable time, especially when R

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 mod- 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 = 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. work in multiple aspects, which are two recent single-server papers relevant to our work. f (x x x, a a a, ξ ξ ξ) dP (B.1a) Letting ρ i , α i,i , and θ be the dual variables associated with constraints (B.1b), (B.1c) and (B.1d) respectively, we present the dual of problem (B.1): where ρ ρ ρ ∈ R N , α α α ∈ R (N +1)×(N +1) and θ ∈ R are unrestricted in sign, and constraint (B.2b) is associated with the primal variable P. Under the standard assumptions that (1) µ d i lies in the interior of the set { S d i dQ : Q is a probability distribution over S}, and (2) µ t i,i lies in the interior of the set { S t i,i dQ : Q is a probability distribution over S}, strong duality holds such that (B.1) and (B.2) equal in optimal objective value (Bertsimas and Popescu, 2005;Shehadeh et al., 2020). Note that for fixed (ρ ρ ρ, α α α, θ), constraint (B.2b) is equivalent to Since we are minimizing θ in (B.2), the dual formulation of (B.1) is equivalent to: This completes the proof.
Appendix C. Proof of Proposition 2 Proof. In view of the objective function (10a), we consider all the terms involving y and define the function H as We are maximizing H(y y y) over a polyhedral set y y y ∈ Y. Note that H is a linear (convex) function in Y and we are maximizing H over a convex compact set Y. From basic convex analysis, we know that there exists an optimal solution at some extreme point of Y. Recall the definition of Y. We Note that an extreme point of Y satisfies (i) y N +2 = −c u N +2 = 0 and (ii) for j ∈ [N + 1], the variable constraint on y j is binding either at the lower bound or the upper bound. If y j is binding at the upper bound, it does depend on y j+1 but if it is binding at the lower bound, it just takes the value of −c u j . From this, we can construct a one-to-one correspondence between the extreme points in Y and partitions of the set {1, . . . , N + 2}. That is, given an interval [k, v] Z in [N + 2], y v is binding at the lower bound, i.e. y v = −c u v , and for j ∈ [k, v − 1], y j is binding at the upper bound, i.e. y j = −c u j + v−1 l=j c w l . Therefore, for notational simplicity, we define We can reformulate our optimization problem over Y to optimizing over partitions of {1, . . . , N + 2}. To do so, let b k,v be the binary variable with 1 indicating that the interval [k, v] belongs to an element of a partition for 1 ≤ k ≤ v ≤ N + 2. The condition With the use of the equality for any extreme point y, we have and the problem maximizing H(y y y) over y y y ∈ Y is equivalent to the following integer program, We remark that the last summation in (C.1a) is summed from j = 2 up to N + 1 instead of N since x i,N +1 = 0 for all i. This makes the last summation term can behave similarly to the previous terms, which simplifies our discussion. Now, we consider the problem With the use of (C.1), we can reformulate it as By the construction of S, we can take maximum over each d i and t i,i . To achieve this, we first gather all the terms involving d i and t i,i separately and take supremum accordingly. We can reformulate the objective function (C.2a) as where we focus on the third to the last term. We first define ∆d i = d i − d i and ∆t i, For the forth term, for any given x x x ∈ X and feasible b, since x i,1 = 1 and b 1,v = 1 for exactly one i and v, say i 1 and v 1 , we have For the third term, notice that for any given x x x ∈ X , i and j have a one-to-one correspondence. That is, for any i, we can identify exactly one j = j i such that x i,j−1 = 1. Then, since b k,v = 1 only for one pair of (k, v), say (k i , v i ). Hence, we have The last term involves the quadratic term x i,j−1 x i ,j . For any given pair (i, i ), there exists at most one j, say j i , such that x i,j−1 x i ,j = 1. If there does not exist such an j i , the triple summation term is just zero. Otherwise, there exists such an j i and the last term reads With the use of these two observations, we have Eventually, we have reformulated the objective function using (C.4), (C.5), (C.6) and (C.7). Note that from (C.2b), the coefficient matrix is totally unimodular and hence, we can relax the integer constraint (C.2c) to b k,v ≥ 0. This reformulates (C.2) into an LP in b b b. One can observe that we put b k,v to the outermost position in all the terms so that we are able to collect the terms with the same b k,v by interchanging the summations. Indeed, for some arbitrary terms z jkv , the triple summation can be written as Then, we are able to take the dual of the resulting LP. Let β j be the dual variable associated to the constraint (C.2b) for j = 1, . . . , N + 2. Thus, we arrive at the following program.
This completes the reformulation.

Appendix D. Proof of Corollary 3
Proof. In view of Propositions 1 and 2, problem (9) (i.e., the inner problem sup in (7)) is equivalent to Noting that the variables α i,0 for i ∈ [N ] appear in the objective function only and it is easy to see that the minimizer is given by λ. Indeed, consider the terms involving α i,0 , i.e., which is non-increasing in α i,0 . Plugging in the optimal α i,0 = λ to the objective function, we arrive at the desired reformulation.

Appendix E. Details of the final MINLP and MILP reformulation of M-DHRAS
From Corollary 3, we can reformulate M-DHRAS as the following MINLP.
s.t. x x x ∈ X , a a a ∈ A, α α α ∈ R (N +1)×(N +1) , ρ ρ ρ ∈ R N , β β β ∈ R N +2 , (E.1b) Note that this is non-linear due to the interaction terms (such as x i,j−1 x i ,j in the objective, ,v x i,1 and ρ i x i,j−1 in the constraints). To linearize this MINLP formulation, we We also introduce McCormick inequalities (E.2a) to (E.2k) for these variables.
The notations · and · are the upper and lower bounds for the variable respectively. Therefore, we obtain an MILP reformulation of the M-DHRAS model.
Note that the McCormick inequalities often rely on big-M coefficients (i.e., variable lower and upper bounds) which take large values and can undermine the computational efficiency. Next, we derive tight bounds of these big-M coefficients to strengthen formulation.
Proof. We only present the case for ρ i and the remaining two cases are similar. From the model in Proposition 1 and its reformulation (C.2) with (C.4) to (C.7), the terms involving ρ i can be summarized as If ρ i ≥ P u 1 , we can simplify the terms as Similarly, if ρ i ≤ P l 1 , then the terms read This shows that there exists a minimizer over ρ i lying in [P l 1 , P u 1 ] and we can set ρ i = P l 1 and ρ i = P u 1 .
Corollary 7. We have lower bounds δ i,j,v = 0, γ i,i ,j,v = 0 and γ 0,i,1,v = 0 with upper bounds Proof. We only present the case for δ i,j,v and the remaining two cases are similar. Recall δ i,j,v = (π j,v − ρ i ) + and from Proposition 6, we have where the lower bound is no greater than 0 and the upper bound is no less than 0.  2019)). Suppose that S is non-empty, convex and compact. Then there exist nonegative constants c 1 and c 2 such that, for all R ≥ 1 and β ∈ (0, min{1, c 1 }), where P R ξ ξ ξ represents the product measure of R copies of P ξ ξ ξ and R (β) = log(c 1 β −1 ) For completeness, we provide the proof of Lemma 1 of Jiang et al. (2019) as detailed in their paper, which is adapted from Theorem 2 in Fournier and Guillin (2015).
To sum up, we arrive at  , 2018). Suppose that the support S is non-empty, convex and compact.
Proof. Recall the dual of f (x x x, a a a, ξ ξ ξ).
f (x x x, a a a, ξ ξ ξ) = max For a fixed pair of (x x x, a a a) ∈ X × A, this is an LP in y y y. Note that Y is bounded, which implies that f (a a a, x x x, ξ ξ ξ) is finite. Also, since X × A and S are bounded sets, the function f (x x x, a a a, ξ ξ ξ) is bounded on X × A × S. Immediately, we have |f (x x x, a a a, ξ ξ ξ)| ≤ M (1 + ξ ξ ξ ), where M is just the upper bound on f . Next, we claim that f (x x x, a a a, ξ ξ ξ) is continuous on X × A × S. For simplicity, write the objective function of the dual problem as b(x x x, a a a, ξ ξ ξ) y, where the i-th entry of b(x x x, a a a, ξ ξ ξ) is the coefficient associated to y i . By fundamental theorem of LP, instead of maximizing over the entire Y, we can maximize over the set of finite extreme points of Y. Hence, f (x x x, a a a, ξ ξ ξ) is the maximum of finitely many linear functions and the continuity follows. Finally, the result follows directly from Theorem 3.6 of Mohajerin Esfahani and Kuhn (2018) jerin Esfahani and Kuhn, 2018). For any β ∈ (0, 1), let (x x x(R, R (β R )),â a a(R, R (β R ))) represent an optimal solution of W-DHRAS with ambiguity set F p (P R ξ ξ ξ , R (β R )). Then, Proof. (adapted from Jiang et al., 2019). By the assumption on support S and Lemma 1, all conditions of Theorem 3.5 in Mohajerin Esfahani and Kuhn (2018) are satisfied. Therefore, the conclusions of Theorem 2 hold valid.
Appendix J. Proof of Proposition 5 Proof. With reference to (C.1), we can reformulate max ξ ξ ξ∈S {f (x x x, a a a, ξ ξ ξ) − ρ ξ ξ ξ −ξ ξ ξ r 1 } as As in the proof of Proposition 2, we group terms with the same t i,i and d i and find an explicit formula for their supremum. First, we define the following quantities for r ∈ [R].
Note that for the terms with t i,0 , Since t i,0 − ρ|t i,0 −t r i,0 | is a piecewise linear function on the intervals [t i,0 ,t r i,0 ] and [t r i,0 , t i,0 ], and it is strictly increasing on [t i,0 ,t r i,0 ], the maximum is attained either att r i,0 or t i,0 . That is, For the terms with t 0,i , we have where the last equality follows from the fact that (π 1,v + λ)t 0,i − ρ|t 0,i −t r 0,i | is a piecewise linear function on the interval [t 0,i ,t r 0,i ] and [t r 0,i , t 0,i ] and the maximum is attained when t 0,i takes any one of the values from t 0,i ,t r 0,i and t 0,i . Using the same argument, we obtain the expressions for the terms with t i,i and d i as follows.
With the use of (J.3) to (J.7), we arrive at the following integer program.
Since the coefficient matrix associated to the constraints (J.8b) is totally unimodular, we can relax the integral constraints on b k,b by b k,v ≥ 0. Taking the dual of the resulting LP in (J.8) with β r j being the dual variable associated to the constraint (J.8b), we arrive at the following LP.
Finally, we can treat u u u r and ν ν ν r as variables and introduce the following constraints.
This completes the proof.

Appendix K. Details of the final MINLP and MILP reformulation of W-DHRAS
Using Proposition 5, we can reformulate W-DHRAS as the following MINLP.
Therefore, we obtain an MILP reformulation of the W-DHRAS model.
s.t. x x x ∈ X , a a a ∈ A, ρ ≥ 0, β β β r ∈ R N +2 , ∀r ∈ [R], (K.3b) Finally, we give tight big-M coefficients in the McCormick inequalities. From the proof of Proposition 5 in Appendix J, we have the explicit formula for the quantities u u u r and ν ν ν r in (J.2).
Immediately, the lower bounds are u r and ν r i,j,v = π j,vd r i . Since the expressions in (J.2) are non-increasing in ρ and we have ρ ≥ 0, the upper bounds are simply u r

Appendix L. Sample averaging approximation approach
In the numerical experiments, we compare DRO models with the sample averaging approximation (SAA) approach. Suppose we have a set of R scenarios {ξ ξ ξ 1 , . . . ,ξ ξ ξ R }. The SAA approach is to solve the original problem by replacing the true distribution with the empirical distribution P R ξ ξ ξ . That is, min x x x,a a a,u u u,w w w

Appendix M. Symmetry-breaking constraints
To enhance the tractability of the models, we introduce the symmetry-breaking constraints.
The idea behind these constraints follows from prior appointment scheduling observations that customers of the same type (i.e., requesting the same service) have the same service time distribution and the observation from Nikzad et al. (2021) that customers within a service region form a basic unit or cluster that shares the same travel time distribution. In this case, the route within the same group does not matter since the service and travel time distributions are the same. This means, with the presence of homogeneous groups, we could be able to eliminate some of the equialent routes and hopefully, improve the model's tractability. In particular, we focus on the presence of one homogeneous group, which is common when services are provided within a service region.
Suppose that except the depot (node 0), the N customers form a homogeneous group. That is, d i are distributionally the same over i ∈ [N ] and t i,i are distributionally the same over i ∈ [N ] and i ∈ [N ] \ {i}. In this case, we argue that the service provider only needs to make decision on the nodes departing from and entering to the depot. We summarize this result formally in the following lemma.
Lemma 8. Assume that there is only one homogeneous group. For a fixed route (i 1 , . . . , i N ) denoted as x x x, i.e. x i j ,j = 1 for j ∈ [N ], define x x x Π by another route with x Π(i j ),j , where (Π(i j )) j∈[N ] is a permutation of [N ] with Π(i 1 ) = i 1 and Π(i N ) = i N . Then, f (x x x, a a a, ξ ξ ξ) is distributionally the same as f (x x x Π , a a a, ξ ξ ξ) for any a a a ∈ A.
Proof. Recall that f (x x x, a a a, ξ ξ ξ) is the second-stage cost, which is a function of the idling time u i for i ∈ [N ], the waiting time and overtime w i for i ∈ [N + 1]. Define the same quantities u Π i and w Π i for the decision x x x Π . We will show, by induction, that w j and u j are distributionally the same as w Π j and u Π j respectively. Note that w 1 and u 1 are distributionally the same as w Π 1 and u Π 1 respectively since by assumption, the first visiting customer is the same. Next, assume that w j−1 and u j−1 are distributionally the same as w Π j−1 and u Π j−1 respectively. Note that By induction assumption and the assumption that the distributions of travel and service times are the same, w j equals w Π j in distribution. It is easy to see that, with the same argument, u j equals u Π j in distribution. This completes the proof.
From Lemma 8, we have E[f (x x x, a a a, ξ ξ ξ)] = E[f (x x x Π , a a a, ξ ξ ξ)] for any permutation with fixed entering and departing node. As an illustration of the symmetry, consider that we have a homogeneous group of 6 customers. Then, the routes (2, 1, 3, 4, 6, 5) and (2, 4, 3, 6, 1, 5) give the same expected cost. We can impose the lexicographic order of the route for the middle customers (not the first and the last customer), i.e., the service provider will first visit customers with a smaller index i. In this example, we only consider the route (2, 1, 3, 4, 6, 5) and eliminate any other possible permutations of the middle customers. Hence, we can introduce the following symmetry-breaking constraints.
x 1,j ≥ x 1,j+1 , ∀j ∈ [2, N − 2] Z (M.1a) x i,j ≤ i−1 l=1 x l,j−1 , ∀i ∈ [2, N ] Z , j ∈ [3, N − 1] Z (M.1b) If customer 1 is neither the first nor the last customer, constraints (M.1a) enforces customer 1 to be the second customer (see example matrix X 1 , where each entry represents x i,j ). Otherwise, the i.e. x i,j = 1 for some j ∈ [3, N − 1] Z , we require that a customer with index less than i must be served at position j − 1 (see X 1 , X 2 and X 3 ). By imposing these constraints, the solver solely determines the entering and departing node without optimizing any permutations for the middle customers. from the same distribution we use in the optimization as discussed in Section 4.1. This simulation assumes that the data comes from the true unknown distribution.
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.).
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 Wang et al. (2020) and perturb the support of the random travel and service times by a parameter δ as [(1 − δ) lower bound, (1 + δ) upper bound], where δ ∈ {0.1, 0.25, 0.5}. A higher value of δ corresponds to a higher variation level.
Set 4. In this set, we assume that we have misspecified the service time distribution in the optimization (i.e., lognormal is not the true distribution). We simulate d i from the U-shaped beta distribution B(0.5, 0.5) on [10, 50] (Jiang et al., 2019), which has the same mean and range as the in-sample distribution.
Set 5. We generate data similar to Set 3 but only changing the service time.
We show the results of the 9 out-of-sample testing sets in each figure in the following order: Similar patterns could be observed as in the case with N = 6. First, M-DHRAS model is the most conservative model, which yields the highest out-of-sample cost for almost all the cases.
Second, we observe that W-DHRAS shares a similar out-of-sample performance as SP (mainly due to a small choice of ) for N = 8 while for N = 10, the W-DHRAS model consistently produces a lower out-of-sample cost than the SP model. Indeed, for N = 10, the operator has a very tight schedule. A small deviation from the in-sample distribution (from the empirical data) results in a significant amount of overtime. Therefore, W-DHRAS could perform better with its ability to hedge against unfavorable scenarios.  the results for the SP model. We observe that all the instances could be solved efficiently (typically less than one minute). We observe that the SP can solve all instances quickly in less than a minute.
This demonstrates that our SP model can solve large instances with an unrealistic number of customers even with many scenarios. Table P.5 shows the computational time for M-DHRAS for different numbers of customers. Recall that this model is not a sample-based model, i.e., it does not depend on the scenarios directly, and hence, we only report the time with R = 5. We observe that the M-DHRAS model can solve all instances efficiently in a few seconds. In Table P.6, we also present the computational time for W-DHRAS with N = 15 under limited data setting R ∈ {5, 10} (where W-DHRAS could yield a better performance). We observe that the model could be solved in a reasonable time (within 4 minutes). In Table P.6, we also present the computational time for W-DHRAS with N = 15 under limited data setting R ∈ {5, 10} (where W-DHRAS could yield a better performance). We observe that the model can solve these instances in a reasonable time within 4 minutes.