key: cord-0232578-0diptjz1
authors: Wang, Ridong; Shehadeh, Karmel S.; Xie, Xiaolei; Li, Lefei
title: Integrated Home Care Staffing and Capacity Planning: Stochastic Optimization Approaches
date: 2022-03-28
journal: nan
DOI: nan
sha: a1a44ba25daf1eba0821edcd850368b46491090e
doc_id: 232578
cord_uid: 0diptjz1

We propose stochastic optimization methodologies for a staffing and capacity planning problem arising from home care practice. Specifically, we consider the perspective of a home care agency that must decide the number of caregivers to hire (staffing) and the allocation of hired caregivers to different types of services (capacity planning) in each day within a specified planning horizon. The objective is to minimize the total cost associated with staffing (i.e., employment), capacity allocation, over-staffing, and under-staffing. We propose two-stage stochastic programming (SP) and distributionally robust optimization (DRO) approaches to model and solve this problem considering two types of decision-makers, namely an everything in advance decision-maker (EA) and a flexible adjustment decision-maker (FA). In the EA models, we determine the staffing and capacity allocation decisions in the first stage before observing the demand. In the FA models, we decide the staffing decisions in the first stage. Then, we determine the capacity allocation decisions based on demand realizations in the second stage. We derive equivalent mixed-integer linear programming (MILP) reformulations of the proposed nonlinear DRO model for the EA decision-maker that can be implemented and efficiently solved using off-the-shelf optimization software. We propose a computationally efficient column-and-constraint generation algorithm with valid inequalities to solve the proposed DRO model for the FA decision-maker. Finally, we conduct extensive numerical experiments comparing the operational and computational performance of the proposed approaches and discuss insights and implications for home care staffing and capacity planning.

capacity planning (HSCP ) problem.

Unfortunately, the HSCP problem is a challenging stochastic optimization problem for several reasons. First, it is a complex combinatorial optimization problem that requires deciding the number of different (cross-trained) caregivers to be hired and allocating (i.e., assigning) the capacity of hired caregivers to different types of services simultaneously. Second, there is significant variability in demand for home care services and service duration. In addition, our analysis of real-life data from a home care service provider suggest that there is a wide range of possible probability distributions for modeling the variability in the demand, indicating distributional ambiguity, i.e., uncertainty in probability distribution (see this analysis in Section 2).

Such uncertainty and distributional ambiguity increase the complexity of modeling and solving the HSCP problem. However, ignoring uncertainty may lead to sub-optimal decisions and the inability to meet customer demand and high costs, among other disappointing consequences of adopting deterministic solutions. Failure to meet customer demand, in particular, may lead to adverse outcomes, especially in healthcare, as it impacts population health. It also impacts customers' satisfaction and thus the reputation of the home service providers.

Motivated by these important issues and motivated by our collaboration with a home service provider in Beijing, in this paper, we propose new stochastic optimization approaches for modeling and solving the HSCP problem under uncertainty. Specifically, given sets of service types, days in the planning horizon, and caregivers and their types (i.e., training), our models solves the following decisions problems simultaneously (a) a staffing problem that determines the number of caregivers to hire within the planning horizon, and (b) an allocation problem that determines the allocation of caregivers capacities (i.e., daily service hours) to service types during the planning horizon. We consider two types of decision-makers: Everything in Advance (EA) and Flexible Adjustment (FA) decision-makers. The (EA) decision-maker decides the number of caregivers to hire and their daily capacity allocation to service types at the start of the planning horizon before the demand and service time are realized. In contrast, the (FA) decision-maker decides the number of caregivers to hire at the beginning of the planning horizon. Then, s/he allocates the hired caregivers to service types once the demand is realized.

To address uncertainty, we propose and analyze two-stage stochastic programming (SP) and distributionally robust optimization (DRO) models for both decision-makers, assuming known and ambiguous distribution, respectively. We obtain near-optimal solutions of our SP models using sample average approximation (SAA). We derive a mixed-integer linear programming (MILP) reformulation of our proposed DRO model for the (EA) decision-maker that can be implemented and efficiently solved using off-the-shelf optimization software. We propose a computationally efficient column-and-constraint generation (C&CG) method to solve our proposed DRO model for the (FA) decision-maker. In addition, we derive valid inequalities to accelerate the C&CG convergence.

We conduct extensive computational experiments comparing the proposed methodologies' computational and operational performances demonstrating where significant performance improvements can be gained. Additionally, we derive several managerial insights relevant to practice. To the best of our knowledge, and according to our literature review in Section 3, our paper is the first to propose stochastic optimization models and methods to model and solve the HSCP problem under uncertainty considering both (EA) and (FA) decision-makers.

The remainder of this paper is organized as follows. In Section 2, we present a motivating example for modeling uncertainty and distributional ambiguity. In Section 3, we review the relevant literature. In Section 4 we formally define the HSCP problem. In Sections 5-6, we present our proposed SP and DRO models. In Section 7, we present our C&CG algorithm and strategies to speed convergence. In section 8, we present our computational results and managerial implications. Finally, we draw conclusions and discuss future directions in Section 9.

This paper is based on our research with Pinetree Care, one of China's leading home care agencies. Pinetree Care provides more than 5 million care sessions to over 600,000 families through a network of more than 10,000 caregivers. One of the main challenges that Pinetree Care and other home care companies face when making staffing and capacity allocation decisions is the uncertainty in service demand and duration. To analyze this uncertainty, we use data from Pinetree Care collected electronically from June 2017 through December 2019, representing ∼205 K care sessions in Beijing.

In the top panel of Figure 1 , we present the empirical and fitted probability distributions for actual demands for rehabilitation nursing and health status assessment, and in Appendix A we present the associated the goodness-of-fit metrics (e.g., Negative of the Log-Likelihood Indeed, future uncertainty is often not distributed in the past. These results support prior numerical studies that show that different distributions can typically explain raw data of uncertain parameters, indicating distributional ambiguity. Finally, these results motivate us to adopt a DRO approach to obtain robust decisions that could hedge against distributional ambiguity.

Operations research in home care service can be divide into three decision levels: strategical level, tactical level and operational level (Restrepo et al., 2020) . In this section, we focus primarily on the literature in tactical level which our problem lies in, i.e., papers that apply stochastic optimization to address HSCP problems similar to ours. For comprehensive recent surveys of operations research methods applied to decisions in home health care, we refer readers to Gutiérrez & Vidal (2013) , Fikar & Hirsch (2017) , Grieco et al. (2021) and Di Mascolo et al. (2021) . We also refer to Michael (2018) for a comprehensive survey of scheduling theory, applications and methods.

Existing stochastic home care service staffing and capacity planning models often consider the demand uncertainty only. SP approaches for HSCP include Rodriguez et al. (2015) , Restrepo et al. (2020) and Zheng et al. (2021) . Rodriguez et al. (2015) proposed a two-stage stochastic programming (SP) approach to solve the staff dimensioning problem in home care services considering random demand. In their two-stage approach, the first stage calculates (near-)optimal levels of resources for possible demand scenarios, while the second stage computes the optimal number of caregiver for each profession. Notably, Restrepo et al. (2020) is the first to study integrated staffing and scheduling problem in the context of home healthcare. They considered staffing and scheduling decisions in the first stage and the temporary reallocation of caregivers to neighboring districts in the second stage. Their study show that the total cost presents an increase when there is less flexibility associated with the allocation of schedules.

The recent work of Zheng et al. (2021) is closely related to ours. Specifically, Zheng et al.

(2021) considered a multi-resource and multi-demand network matching problem in which they assume deterministic service duration and that the demands of services follow a fully known probability distribution. Accordingly, Zheng et al. (2021) proposed a two-stage stochastic MILP which seeks the optimal service authorization and capacity planning decision in first stage, then the service resource allocation in the second stage. In their model, the goal is to maximize the total expected revenue netted by the authorization cost, payment and the penalty cost of service shortage. However, they do not incorporate the penalty cost of overstaffing in the objective, which is an important source of operational expenses (Restrepo et al., 2020) . In addition, ignoring random service duration may lead to sub-optimal solutions.

These SP models assume that the decision-makers are risk-neutral and have complete knowledge about the underlying uncertainty through a known probability distribution (Shapiro et al., 2014) . Hence, the applicability of the SP approach is limited to the case in which the decision-makers have sufficient data to model the distribution of random parameters which is implausible in practice. Given the distributional ambiguity (see our analysis in Section 2), if we calibrate a model to a misspecified distribution, the resulting optimal SP decisions may have a disappointing out-of-sample performance under the true distribution which is termed the optimizer's curse (Smith & Winkler, 2006) .

Robust optimization (RO) is alternative technique to model, analyze and optimize decisions under uncertainty and ambiguity (where the distributions are unknown). It assumes that the decision maker is risk-averse and has no distributional knowledge about the underlying uncertainty, except for its support (partial information), and the model minimizes the worstcase cost over an uncertainty set (Ben-Tal et al., 2009) . In home care service area, several studies propose robust approaches to address the uncertainty (Carello & Lanzarone, 2014; Cappanera et al., 2018; Shi et al., 2019; Cappanera & Scutellà, 2021) . Notably, Cappanera et al. (2018) addressed uncertainty of patient demand over a multiple-day time horizon and jointly studied the assignment, scheduling and routing decisions via a non-standard cardinalityconstrained robust approach. Shi et al. (2019) proposed a robust model for both uncertain service and travel time and discuss the heuristic solution approaches to find the optimal route for each caregiver and appointment time for each patient.

Although RO is a powerful technique in dealing with uncertainty in optimization, its solutions can be too conservative (Thiele, 2010; Roos & den Hertog, 2020) . Additionally, RO does not fully utilize any distributional information of uncertainty. Hence, RO may yield poor expected performance and sub-optimal decisions for other more-likely scenarios. For comprehensive recent surveys of robust optimization, we refer to Bertsimas et al. (2011 ), Gabrel et al. (2014 and Gorissen et al. (2015) .

Distributionally robust optimization (DRO) is another approach for modeling uncertainty that bridges the gap between the conservatism of robust optimization and the specificity of stochastic programming (Goh & Sim, 2010) . In DRO, we assume only partial distributional information is known. Then, we model the distribution of uncertainty as a decision variable that belongs to an ambiguity set (i.e., a family of distributions characterized through certain known properties of the unknown data-generating distribution, Esfahani & Kuhn, 2018) . Optimization is based on the worst-case distribution within the ambiguity set. And we can use easy to compute information such as the mean and range of random parameters to construct the ambiguity sets and build DRO models. We refer to Rahimian & Mehrotra (2019) for a comprehensive survey on DRO techniques.

DRO has the following striking benefits. First, it adopts a worst-case approach regularizing the optimization problem and thereby mitigating the optimizer's curse characteristic for SP.

Second, DRO avoids the well-known over-conservatism and the poor expected performance of RO. Finally, DRO models are often more computationally tractable than SP and RO (Esfahani & Kuhn, 2018) . We observe that despite the potential advantages, only one study, Tsang & Shehadeh (2021) , using DRO approach to hedge against uncertainty in home care service area.

However, their study focus on the routing and appointment scheduling problem for a singleprovider. To the best of our knowledge, there is no DRO approach for the specific HSCP problem that we address in this paper.

In this work, we develop two DRO models to consider the realistic lack of distributional information and compare these models to SP models. We compare our work with Zheng et al.

(2021) and Restrepo et al. (2020) , which are recent studies relevant to our work. First, the three papers both address stochastic demand. However, Zheng et al. (2021) and Restrepo et al. (2020) assumed that the service time is deterministic and do not address the potential distributional ambiguity of demand which is captured by real world data (see our analysis in Section 2). As mentioned earlier, this may lead to sub-optimal and unrealistic solutions.

In contrast to them, we extend the considered problem by incorporating uncertainty and distributional ambiguity for demand and service duration. Second, we propose different SP and DRO models for two types of decision-makers (EA and FA). Third, Zheng et al. (2021) and Restrepo et al. (2020) assumed that the caregivers have single skill. In contrast to these papers, we consider cross-trained caregivers representing flexible capacity. To the best of our knowledge, our work is the first to consider cross-training in the HSCP problem.

Let us now start to introduce our HSCP problem. Suppose that there is a home care company that provides L different types of home care services (e.g., rehabilitation nursing, health status assessment, etc. caregivers is h k . The fixed cost of hiring type k caregiver is c k .

The demand, d l,t , for service type l ∈ [L] in each day t ∈ [T ] is random. Service time, s l,t , for each service type l ∈ [L] in each day t ∈ [T ] is also random. Due to the uncertainty of demand and service time, one or multiple of the following scenario can happen: (1) overstaffing, or and, (2) under-staffing (i.e., failure to meet customers demand). The over-staffing and under-staffing unitary penalty costs are respectively c o l,t and c u l,t , for all l ∈ [L] and t ∈ [T ]. A complete listing of the problem parameters can be found in Table 1. Given T , K, and L, our models determine the number of caregivers to hire and the capacity of type k caregivers allocated to type l ∈ [L] in each day t ∈ [T ]. In particular, we consider two types of decision-makers. The first decides the number of caregivers to hire of each type and their daily capacity allocation to service types at the beginning of the planning horizon before the demand and service time are realized. We call this decision-maker Everything in Advance (EA) decision-maker. The second decision-maker decides the number of caregivers to hire before realizing the demand (i.e., at the beginning of the planning horizon). S/he then allocates these caregivers to service types once the demand is realized. We call this decision maker Flexible Adjustment (FA) decision-maker.

In the next sections, we propose two-stage SP and DRO models for (EA) and (FA), considering the case of known and unknown distributions of demand and service time, respectively. the skills set of type k caregivers c k the fixed cost for hire type k full-time caregivers c k,l,t the unitary service cost of type k full-time caregivers for serve type l customers in day t c u l,t the unitary penalty cost due to shortage of service capacity of type l service c o l,t the unitary penalty cost due to over-staffing of service capacity of type l service w the upper bound number of full-time caregivers w the lower bound number of full-time caregivers h k the service capacity (service hours) of one type k full-time caregiver in one day Scenario-dependent Parameters d l,t the demand quantity for service type l in day t s l,t the service time for service type l of caregiver in day t

In the first stage of the (EA) model, we decide the number of caregivers to hire from each type k ∈ [K] and their capacity allocation to service type l ∈ [L] in each day t ∈ [T ]. And in the second stage, after the demand and service time are realized, we compute over-staffing and under-staffing. The objective is to minimize the total costs: fixed staffing (or hiring) and allocation (assignment) costs and a weighted sum of random costs related to over-staffing and under-staffing. In the first stage of the (FA) model, we decide the number of caregivers to hire, and in the second stage we decide capacity allocation and compute over-staffing and understaffing. The objective of this model is to minimize the total costs: fixed hiring cost and a weighted sum of the random cost related to capacity allocation (assignment) cost, over-staffing cost, and under staffing cost. 

where for each feasible decision x, y and a joint realization of uncertain parameters ξ := [d, s] , our second-stage (recourse) formulation is as follows:

Formulation (1) seeks to find staffing and allocation decisions (x, y x, y x, y) that minimize the sum of (1) fixed employment and assignment costs (first two terms), and (2) expectation of the second stage (recourse) random cost Q(x , y , ξ) (third term) subject to uncertainty ξ ∼ P.

The second stage recourse function minimize the random over-staffing (first-term) and understaffing (second-term) costs. Constraints (1b) ensure that the total number of staff is at least w and at most w (these can be adjusted by the decision-maker depending on relevant practical constraints). Constraints (1c) ensure that capacity allocations to service types does not exceed the available capacity. Constraints (1d) specifies feasible ranges of the decision variables.

Second stage constraint (2b) computes the over-staffing and under-staffing.

In this section, we present our proposed two-stage DRO formulation of the HSCP with (EA) decision-maker, which does not assume that the probability distributions of demand and service durations are known. We, however, assume that we know the mean and support of these random parameters. (Recall that these parameters can be estimated based on, e.g., expert knowledge, where support could represent the dispersion or range of variability that we seek protection against). First, let us define some additional parameters, sets, and notations defining our ambiguity set and DRO model. We let µ d l,t and µ s l,t respectively represent the mean value of d l,t and s l,t , for all l ∈ [L] and t ∈ [T ]. The random vector ξ is a measurable function ξ : Ω → U with a measurable space (Ω, F), where U = U d ×U s is the bounded support of ξ. U d and U s are respectively the support set of d d d and s s s defined as follows:

We define P(U) as the set of probability distributions supported on U, and E P as the expectation under P. Finally, we donate µ := E P [ξ] = [µ d , µ s ] . Using this notation, we construct the following mean-support ambiguity set:

Using the ambiguity set F(U, µ), we formulate our DRO model of the HSCP with (EA)

decision-maker (denoted as E-DHSCP) as follow:

where X := {(x, y): (1a)-(1d)} is the feasible region of first stage decision variables. Formulation (6) seeks first stage decisions x and y that minimizes the summation of fixed hiring cost (first term), the capacity allocation cost (second term) and the worst-case (maximum) expectations of the second-stage operational cost (third term), where the expectation is taken over all distributions residing in F(U, µ).

Recall that Q(·) is defined by a minimization problem; hence, in (6), we have an inner max-min problem. As such, it is not straightforward to solve (6) in its presented form. In this section, we derive an equivalent formulation of the min-max model (6) that is solvable. First,

in Proposition 1, we present an equivalent reformulation of the inner maximization problem sup P∈F (U ,µ) E P [Q(x , y , ξ)] in (6) (see Appendix B for a detailed proof).

Proposition 1. For any (x, y) ∈ X , problem sup

Again, the problem in (7) involves an inner max-min problem that is not straightforward to solve in its presented form. Next, we use the structural properties of the recourse Q(·) to derive an equivalent linear programming reformulation of the inner maximization problem in (7). Note that for fixed x, y, ξ, Q(x , y , ξ) is a linear program (LP). The dual formulation of Q(x , y , ξ) is as follows:

where ρ l,t is the dual variables associated with constraints (2b). In view of the dual formulation in (8), the inner problem max (d,s)∈U Q(x, y, ξ) − L l=1 T t=1 (d l,t α l,t + s l,t β l,t ) in (7) is equivalent to problem (9) in Proposition 2 (see Appendix C for a detailed proof).

Proposition 2. For fixed x, y, α, and β, problem max

Replacing the inner problem max (7) with its equivalent reformulation in (9), then combining with the outer minimization problems in (7) and (6), we derive the following equivalent MILP reformulation of the E-DHSCP model.

(10c)

6. Models for the (FA) Decision-maker 6.1. Two-stage SP for (FA)

In this section, we present our proposed two-stage SP formulation of the HSCP problem considering the (FA) decision-maker. Specifically, we keep x k , ∀k ∈ [K] as first stage decisions, and move y k,l,t to the second-stage, for all k ∈

. In addition, we define nonnegative continuous variables o k,t representing the capacity surplus of type k ∈ [K] caregiver

. This is to account scenarios where the total capacity of hired caregivers is larger than the total realized demand. Accordingly, our F-SP model can now be stated as follows:

where for each feasible decision x and a joint realization of uncertain parameters ξ := [d, s] , our second-stage (recourse) formulation is as follows:

Formulation (11) seeks to find staffing decisions (x x x) that minimize the sum of (1) fixed employment cost (first term), and (2) expectation of the second stage (recourse) random cost Q(x, ξ) (second term) subject to uncertainty ξ ∼ P. The second stage recourse function minimize the random assignment cost (first-term), over-staffing (second-term) and under-staffing (third term) costs. Constraints (12b) compute the over-staffing. Constraints (12c) compute the under-staffing.

In this section, we present our proposed two-stage DRO formulation of the HSCP with (FA) decision-maker. Using the ambiguity set F(U, µ) defined in (5), we formulate our DRO model of the HSCP with (FA) decision-maker (denoted as F-DHSCP) as follow:

where X 2 := {(11b) − (11c)} is the feasible region of first stage decision variables. Formulation (13) seeks the first stage decisions x that minimize the sum of fixed employment cost (first term) and the worst-case (maximum) expectations of the second-stage operational cost (second term), where the expectation is taken over all distributions residing in F(U, µ).

Using the same techniques in Proposition 1, we derive the following equivalent reformulation (13).

Note that it is not straightforward to solve formulation (14) in its presented form due to the inner max-min problem. Therefore, our goal is to derive an equivalent formulation of (14) that is solvable. As a first step, we re-write the recourse Q A (x, ξ) using its equivalent dual formulation as (note that for fixed x and ξ, Q A (x, ξ) is a bounded LP):

where ρ l,t and λ k,t is the dual variables associated with constraints (12b) and (12c). Replacing (14) with its dual formulation in (15), we derive the following equivalent reformulation of (14).

where h(x, d, s, α, β) that is solvable (see Appendix D for a detailed proof).

Proposition 3. Let ∆d l,t = d l,t − d l,t and ∆s l,t = s l,t − s l,t , for all l ∈ [L] and t ∈ [T ]. Then, fixed x, α, and β, solving max

h(x, d, s, α, β) is equivalent to solving the following MILP:

Note that McCormick inequalities (18f)-(18k) in formulation (18) involve big-M coefficients ρ l,t , for all l ∈ [L] and t ∈ [T ], that can undermine computational efficiency if they are set too large. Therefore, in Proposition 4, we derive tight lower bounds (i.e., ρ l,t ) on variables ρ l,t to strengthen the MILP reformulation (see Appendix D for a detailed proof).

Proposition 4. ρ l,t = min (16) with (18), then combining the inner problem in the form of (17) with the outer minimization problem in the F-DHSCP model in (13), we derive the following equivalent reformulation of the F-DHSCP model:

d l,t α l,t + ∆d l,t g l,t α l,t + s l,t β l,t + ∆s l,t z l,t β l,t + d l,t s l,t ρ l,t + ∆s l,t d l,t r l,t + ∆d l,t s l,t q l,t + ∆d l,t ∆s l,t e l,t : (18b) − (18l) ,t ∈ [T ].

(19b)

In this section, we present a Monte Carlo Optimization algorithm to obtain near-optimal solutions to the SP models (Section 7.1), and a column-and-constraint generation (C&CG) algorithm to solve the F-DHSCP model (Section 7.2). Note that the E-DHSCP model can be implemented and solved directly using off-the-shelf optimization software.

Note that it is difficult to obtain an exact optimal solution to the two-stage SP models in (1) and (11). Hence, we first present a Monte Carlo approach base on Homem-de Mello &

Bayraksan (2014), Kleywegt et al. (2002) and Shehadeh et al. (2021) to obtain near-optimal solutions. We replace the distribution of random parameters with a empirical distribution N independent and identically distributed (i.i.d) samples of services demands and durations at first and then we solve the sample average approximation (SAA) formulations of (1) and (11) in Appendix G.

Algorithm 2 in Appendix G summarizes the Monte Carlo Optimization (MCO) algorithm.

We start with a small initial sample size. In first step, we will solve the SAA formulation with designated sample size N and evaluate the solution x k N by Monte Carlo simulation with the generated N random samples and record the objective value v k N of SAA and v k N of Monte Carlo simulation. We will run step 1 with K replicates. In second step, we calculate the average of v k N and v k N as v N and v N . According to Mak et al. (1999) , v N and v N are statistical lower and upper bound of the optimal value of original formulation. In step 3, we calculate the approximate optimality index as

AOI N approximately estimate the optimality gap of lower bound and upper bound of original formulation. Therefore, if AOI N is smaller than the very small predetermined termination tolerance (note that can be zero), the algorithm will terminate and output the termination sample size N , objective value v N , v N and approximate optimality index AOI N . Otherwise, we will increase the sample size (N ← 2N ) then go to step 1.

The F-DHSCP formulation (19) involves a maximization problem in constraints (19b) and

thus we cannot solve it directly using standard techniques. In this section, we develop a C&CG algorithm to solve the F-DHSCP model (19). The motivation of this algorithm is as follows.

Considering the inner maximization problem max

Similarly, the optimal s * ∈ U s would only take value s s s or s s s. Thus, the inner maximization over (d, s) ∈ U could be taken over all the combination of d d d and s s s.

However, only a few critical scenarios of d and s play a role in deriving an optimal solution to the F-DHSCP model. As such, instead of solving a problem with exponentially many constraints, we use C&CG to identify these scenarios and obtain an optimal solution.

Algorithm 1 presents our C&CG algorithm which is based on the original C&CG of Zeng & Zhao (2013) and is implemented in a master-sub problem framework. The master problem (20) employs scenario-based constraints (20d)-(20f) and the sub-problem generates demand and service time scenarios and optimality cuts from the worst-case distribution. The algorithm proceeds as follows. At each iteration, we first solve the master problem, which includes only a subset of demand and service time scenarios, and obtain the optimal solution (i.e., hiring decisions). Thus, by considering a subset of demand and service time scenarios, the master problem is a relaxation of the original problem and thus provides a lower bound on the optimal value of F-DHSCP. Then, given the optimal solution from the master problem, we identify demand d * d * d * and service time s * s * s * scenarios by solving the sub-problem. Note that solutions of the master problem are feasible to the original problem, then by solving the sub-problem using these solutions, we obtain an upper bound. Third, we pass the identified demand and service time scenarios to the master problem on an as-needed basis, along with introducing relevant second-stage variables and constraints. Then, we solve the master problem again with the new information (a larger set of scenarios O) from the sub-problems. This process continues until the gap between the lower and upper bound on the optimal value of the problem obtained in each iteration satisfies a predetermined termination tolerance (note that we set this tolerance to 0.02). Recall that the F-DHSCP formulation's recourse (second-stage) problem is feasible for every feasible first-stage decision and ((d, s) ∈ U). Thus, we do not need to add any feasibility cuts (Zeng & Zhao, 2013) . Moreover, we only have a finite number of scenarios.

Thus, the algorithm terminates in a finite number of iterations (see Zeng & Zhao, 2013 for detailed proofs).

Although C&CG theoretically provides a good way to solve our problem, it sometimes does not perform well from a computational viewpoint. In this section, we aim to incorporate more second-stage information into the master problem by exploiting the properties of the recourse problem.

(1) Valid lower bound (LB) inequalities

Algorithm 1 Column-and-constraint generation (C&CG).

Step 1.

Step 2. Master Problem. Solve the following master problem:

Record an optimal solution (x * , α * , β * , δ * ) and set LB = Z * .

Step 3: Solve sub-problem.

3.1 With (x, α, β) fixed to be (x * , α * , β * ), solve the following problem (21) for all t ∈ [T ].

3.2 Record the optimal solution (ρ * t , λ * t , g * t , z * t ) and optimal value W * t .

≤ ε stop and return x * as the optimal solution of (13); else go to step 4.

Step 4. Column-and-Constraint Generation Routine

, and the following constraints to the master problem:

Observe that the Dirac measure on µ µ µ lies in F(U, µ µ µ). Therefore, we have

where the lower bound is a deterministic problem with a single scenario ξ ξ ξ = µ. This means we can impose the following constraint

which serves as a global lower bound on the objective of the DRO model. Second, from (22), for any given first-stage decision, the recourse value with ξ ξ ξ = µ provides a lower bound.

Therefore, in the initialization step of the C&CG method, we could include the scenario ξ ξ ξ = µ.

(2) Valid upper and lower bound of variables α and β

According to the structure of (16), we propose valid upper and lower bounds on variables α and β in Proposition 5 (see Appendix F for a detailed proof). In Section 8.4, we demonstrate the computational advantages gained by incorporating these valid inequalities.

Proposition 5. α l,t and β l,t are respectively valid upper bound on variables α l,t and β l,t . α l,t and β l,t , defined in (24), are respectively valid lower bound on variables α l,t and β l,t .

In this section, we compare the computational and operational performance of the proposed models for each decision maker and derive insights into the HSCP problem. In Section 8.1, we describe the set of test instances that we constructed and discuss other experimental setups.

In Section 8.2, we obtain near-optimal solutions to the E-SP and F-SP via their SAA formulations. In Section 8.3, we compare solution times of the proposed models. In Section 8.4, we demonstrate efficiency of the proposed valid inequalities for the F-DHSCP model. In Section 8.5, we compare optimal solutions of the E-SP and E-DHSCP models and their out-of-sample performance. Finally, we compare the optimal solutions of the F-DHSCP and F-SP models and their out-of-sample performance in Section 8.6.

To the best of our knowledge, there are no standard benchmark instances in the literature for the HSCP problem. Therefore, we construct 15 HSCP instances based on the same param- service types, we have c u l,t + c o k,t > c k,l,t for all k, l and t. This is because it is always optimal to set y k,l,t = 0 when c u l,t + c o k,t < c k,l,t . We followed the same procedures in the stochastic scheduling, home care service, and Similarly, for each instance, we generated N scenarios of s n 1,1 , . . . , s n L,T , for n = 1, . . . , N , by following LogN distributions with the generated µ µ µ s and σ σ σ s . We truncated these scenarios on [s, s]=U [20, 80] . We round each generated parameter to the nearest integer (Tsang & Shehadeh, 2021) . We chose the LogN distribution to generate the in-sample data because it is typically used to model service time and demand in the home care, scheduling, and other literature (Cayirli et al., 2006; Jiang et al., 2019; Tsang & Shehadeh, 2021) .

According to staffing standards for various types of elderly care institutions in China, the minimum staffs number for serving less than 100 customers is three or five. Hence, we set the minimum number of caregivers, w, to 3. We implement all models and solution methods using the Python 3.8.8 programming language calling Gurobi V9.1.2 as a solver with default settings. We ran all experiments on a computer with Intel(R) Xeon(R) E3-1231-v3 processor with two 3.40GHz CPUs and 32 GB shared RAM memory.

For each instance in Table H .1, we implemented the MCO algorithm in Section 7.1 with N o = 10, N = 1000 and K = 10. 

) and the 95% confidence intervals (95%CI) of the statistical lower and upper bounds (v N and v N , receptively) on the objective values of the E-SP and F-SP models for each instance under (c u l,t , c o l,t (c o k,t )) = (20, 2) and truncated LogN distribution. We observe similar results under different choices of these cost parameters. In the next section, we present and analyze the ranges of solution time of the SAA formulations.

Herein we summarize the key findings of the MCO algorithm. Clearly, N = 320 and N = 100 scenarios of the demand and service duration are sufficient to obtain a near-optimal solution to the E-SP and F-SP models via their SAA formulations, respectively. First, AOI N at N = 320 and N = 100 ranges from 0 to 0.0098 and from 0 to 0.0055 for the E-SP and F-SP models, (100) are very tight (i.e., have a small variance). These results qualify v E N and v F N as tight statistical estimates for the lower and upper bounds on the optimal value of each instance, respectively. In addition, we also use truncated normal distribution and uniform distribution to generate the sample data. We obtain similar conclusion for these distributions (see Appendix I and Appendix J).

It is worth mentioning that solving these models with larger N resulted in longer solution times without consistent and significant improvements in the relative approximation gaps.

However, as shown in Appendix I-Appendix J, solution times are still reasonable under N = 500 scenarios. And the AOI N is close to the value at termination sample size.

In this section, we analyze the solution times of the proposed SP and DRO models. In addition to the based range of the demand described in Section 8.1, we consider another range [40, 100] , which corresponding to a higher amound of demand and variation level. For each of the 15 instances in Table H .1 and demand range, we generated and solved 10 instances as described in Section 8.1 for a total of 300 random instances. Then, we solved these instances using the proposed models for each decision-maker.

Let us first compare the computational performance of the E-SP and E-DHSCP models. Table 2 presents the minimum (Min), average (Avg), and maximum (Max) CPU time (in seconds) for these models. We first observe that the computational effort tend to increase with the scale (L × K × T ). Second, both models can quickly solve all of the instances to optimality. However, the E-DHSCP model takes a shorter time than the E-SP model to solve each instance. These results are consistent under both ranges of the demand.

Next, we compare the solution times of the F-SP and F-DHSCP models presented in Table 3. We make the following observations from this table. First, it is clear that solution times of these models are much longer than the E-SP and E-DHSCP models. For example, the average solution time of instance 15 (largest instance) using the E-SP and F-SP models Overall, these results demonstrate the computational efficiency of the proposed approaches.

To confirm our conclusion for the scenario-based SAA models, we solve the SAA formulations of E-SP and F-SP models with 500 samples. We solve 6 instances with L = 6 under these 

In this section, we investigate the efficiency of the proposed valid inequalities (23)- (24) for the E-DHSCP and C&CG in Section 7.2.1. Given that we observe similar results for all instances, for brevity and illustrative purposes, we use Instance 15 (L = 6, K = 8, T = 180) in this experiment. We separately solve the E-DHSCP model via C&CG with (w/) and without It is obvious that both the lower bound and gap values converge faster when we introduce inequalities (23)-(24) into the master problem. Moreover, because of the better bonding effect, the algorithm converges to the optimal solution in fewer iterations and shorter solution times.

For example, the algorithm takes 531 seconds and 19 iterations to solve this instance w/ inequalities (23)-(24) and terminates at the specified time limit (i.e., 7200 seconds) with gap 99.95% w/o these inequalities.

In this section, we analyze the optimal staffing patterns of the E-SP and E-DHSCP models (Section 8.5.1) and their out-of-sample performances (Section 8.5.2). For illustrative purposes and brevity, we present results for Instance 7 (L = 4, K = 8, T = 30) and Instance 10 (L = 6, K = 6, T = 30), which corresponds to middle-sized instances. In addition, we consider the following types of caregivers. In Instance 7 (which consists of K = 8 types of caregivers), we let caregivers types k ∈ {1, . . . , 4} represent specialized caregivers, meaning they can only provide one type of service. And we let type k ∈ {5, . . . , 8} represents cross-trained caregivers who are trained or qualified for two types of services. In Instance 10 (which consist of K = 6 types of caregivers), all caregivers are cross-trained.

In this section, we compare the optimal staffing pattern (i.e., the optimal number of hired We observe the following from Figure 3a . First, both models hire more caregivers under range 2 than range 1. This makes sense because, in range 2, both the volume and variability of the demand are higher. Second, it is clear that both models hire more caregivers as the understaffing cost increases from 1 to 25. However, the number of hired caregivers stabilize after c u l,t = 10. Third, when c u l,t ≤ 5 (i.e., low under-staffing cost), both models hire approximately the same number of caregivers. In contrast, when c u l,t > 5, the E-DHSCP models hire more caregivers than the E-SP model to mitigate excessive under-staffing. These results indicate that c u l,t = 10 may be a suitable penalty cost of under-staffing for (EA) decision-makers. We observe a different pattern when we fix c u l,t = 10 and increase the penalty for over-staffing c o l,t from 1 to {5, 10, 15, 20, 25}. Specifically, as shown in Figure 3b both models hire less caregivers as increasing c o l,t . And the number of hired caregivers stabilize after c o l,t = 10. Next, we analyze the detailed staffing patterns under cost structure (c u l,t , c o l,t ) = (10, 1) presented in Figure K .2 (see Appendix K). First, if there are only cross-trained caregivers (Instance 10), all models tend to hire all types of caregivers. Second, if there are both specialized and cross-trained caregivers (Instance 7), the E-SP model tend to hire specialized caregivers under both demand ranges while the E-DHSCP model hire a mixture of specialized and crosstrained caregivers. In addition, the E-DHSCP model hire more cross-trained caregivers under demand range 2 than under demand range 1. This could indicate that cross-trained caregivers are more favorable under a higher variation level.

In this section, we compare the operational performance of the optimal solutions of the E-DHSCP and E-SP via out-of-sample simulation. 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. First, we solve each instance and obtain the optimal first-stage decisions. Second, we fix the obtained first-stage decisions, then re-optimize the second stage of the SP using the following sets of N = 10, 000 out-of-sample data to compute the second stage cost (i.e., under-staffing plus over-staffing costs). We compute the out-of-sample total cost as first-stage cost plus the out-of-sample second-stage cost.

Set 1. We assume that we have perfect distributional information. That is, we use the same distributions (LogN) and parameter settings that we use for generating the N data in the optimization to generate N data. This is to simulate the performance when the true distribution is the same as the one used in the optimization.

Set 2. In this set, we assume that we have used misspecified distributions of the demand and service duration in the optimization (motivated by our analysis in Section 2). We follow a similar out-of-sample simulation testing procedure described in Shehadeh (2020) and Wang et al. (2020) and perturb the support of the random demand and service durations by a parameter ∆ and used parameterized uniform distribution U [(1 − ∆) lower bound, (1+∆) upper bound] to generate the N = 10, 000 data. We apply ∆ ∈ {0, 0.1, 0.25, 0.5}.

A higher ∆ corresponds to a higher variation level.

For brevity, we focus on Instance 10 and present results under (c u l,t , c o l,t ) = (10, 1). We observe similar results for Instance 7 (see Appendix L.2). First, let us analyze simulation results under Set 1 (i.e., perfect distributional information case). Figure L.1 in Appendix L.1 presents histograms of the out-of-sample total cost, second-stage costs, and under-staffing. It is clear from this figure that the E-DHSCP model yields a higher total cost and lower secondstage cost than the E-SP model on average and at all quantiles. This makes sense because, as we discussed in Section 8.5.1, the E-DHSCP tends to hire more caregivers to hedge against under-staffing (reflected by significantly lower under-staffing on average and all quantile than the E-SP model in Figure L.1c) . However, this results in a higher fixed cost associated with hiring more caregivers. Next, we investigate the value of distributional robustness from the perspective of out-ofsample disappointment, which measures the extent to which the out-of-sample cost exceeds the model's optimal value (Van Parys et al., 2021; Wang et al., 2020) . Following Shehadeh (2020) and Wang et al. (2020) , we define the out-of-sample disappointment as max 0, Vout−Vopt Vopt × 100%, where V opt and V out are as the model's optimal value and the out-of-sample objective value, respectively. A disappointment of zero indicates that the model's optimal value is equal to or larger than the out-of-sample (actual) cost, suggesting that the model is more conservative and avoids underestimating costs. In contrast, a larger disappointment implies a higher level of over-optimism because, in this case, the actual cost (V out ) of implementing the optimal solution of a model is larger than the estimated cost (V opt ). These results confirm that the DRO model provides a more robust estimate of the actual cost that we will incur in practice. In a decision-making context where the goal is to minimize expenses, disappointments are more harmful than positive surprises (overestimated costs).

In this section, we analyze the optimal staffing patterns of the F-SP and F-DHSCP models (Section 8.6.1) and their out-of-sample performance (Section 8.6.2). We present the results for

Instance 7 and Instance 10 as defined in Section 8.5.

In this section, we compare the optimal staffing patterns of the F-SP and F-DHSCP under different values of the under-staffing c u l,t and over-staffing c o k,t penalty costs. In Figure 6a , we present the optimal number of hired caregivers under demand range 1 and range 2 with c o k,t = 1 and c u l,t ∈ {2, 3, 4, 5, 6, 7} for Instance 10. We observe the same staffing patterns for Instance 7 (see Figure M .1 in Appendix M).

We observe the following from Figure 6 . First, both models hire more caregivers under demand range 2 than under range 1 and the number of hired caregivers stabilize after c u l,t = 5. Second, when c u l,t ≤ 3, both models hire approximately the same number of caregivers. In contrast, when c u l,t > 3 the F-DHSCP models hire more caregivers than the E-SP models. These results indicate that c u l,t = 5 may be a suitable penalty cost of under-staffing for the (FA) decision-makers. We observe a different pattern when we fix c u l,t = 5 and increase the penalty for over-staffing c o k,t from 2 to {3, 4, 5, 6, 7}. Specifically, as shown in Figure 6b both models hire less caregivers as c o k,t increases. And the number of hired caregivers stabilize after c o k,t = 3. Next, we analyze the detailed staffing patterns under cost structure (c u l,t , c o k,t ) = (5, 2) presented in Figure M .2 (see Appendix M). Compared with EA models, both F-SP and F-DHSCP models prefer cross-trained caregivers under two demand ranges for Instance 7. Figure 7a , the F-DHSCP model yields a lower total cost than the F-SP model on average and at all quantiles. This makes sense because F-DHSCP make more flexible capacity allocation, which reduce the number of hired caregivers then the fixed cost associated with hiring more caregivers.

Finally, it is clear from Figure 8 that the F-DHSCP solutions yield smaller out-of-sample disappointments on average and at all quantiles than the F-SP solutions under all levels of variation (∆). Notably, the disappointment of the F-DHSCP solutions is smaller than 50% under a higher variation ∆ = 0.5. These results confirm that the F-DHSCP model provides a more robust estimate of the actual cost that we will incur in practice.

In this paper, we proposed stochastic optimization methodologies for the HSCP problem, which seeks to find the number of caregivers to hire and their allocation decisions to different types of services in each day within a specified planning horizon. The objective is to minimize the total cost associated with staffing (i.e., employment), capacity allocation, over-staffing, and under-staffing. To address uncertainty in service duration and the demand for home service, we proposed and analyzed SP and DRO models considering two types of decision-makers, namely everything in advance decision-maker (for which we proposed E-SP and E-DHSCP models) and a flexible adjustment decision-maker (FA) (for which we proposed (F-SP and F-DHSCP). We derive equivalent MILP reformulations of the proposed nonlinear E-DHSCP model for the EA decision-maker that can be implemented and efficiently solved using off-the-shelf optimization software. We propose a computationally efficient column-and-constraint generation algorithm with valid inequalities to solve the proposed F-DHSCP model for the FA decision-maker.

Finally, we conduct extensive numerical experiments comparing the operational and computational performance of the proposed approaches and discuss insights and implications for home care staffing and capacity planning. The results demonstrate (1) the DRO models tend to hedge against uncertainty and ambiguity by hiring more caregivers and thus have a higher fixed cost but significantly lower operational (i.e., under-staffing) costs than the SP models, (2) the DRO solutions are more robust under various distributions with substantially lower disappointments and thus higher reliability than the SP models, (3) the computational efficiency of the proposed approaches, indicating their applicability in practice. We also derive insights into the HSCP under different settings. In particular, we show how one can use the proposed models to obtain staffing and capacity allocation decisions under different cost structures that the decision-maker can easily change in the models. These results address the primary objective of this paper. That is, to investigate the value of SP and DRO approaches for this specific HSCP problem. More broadly, our results motivate the need for considering different approaches when modeling uncertainty and draw attention to the need to model the distributional ambiguity of uncertain problem data in strategic real-world stochastic optimization problems such as the HSCP problem.

We suggest the following areas for future research. First, we want to extend our approach by incorporating other multi-modal sources of uncertainty (e.g., travel time). Second, our models can be considered as the first step toward building data-driven and robust home care service planning including the creating routing plans for caregivers and assigning appointment times for customers, among other important considerations. Finally, we aim to incorporate other aspects, such as the possibility of hiring part-time caregivers and the unexpected availability of the hired full-time caregivers on the day of service.

Zheng, C., Wang, S., Li, N., & Wu, Y. (2021) . Stochastic joint homecare service and capacity planning with nested decomposition approaches. European Journal of Operational Research, . Proof. For fixed x and y, we can formulate problem sup P∈F (U ,µ) E P [Q(x , y, ξ)] as following linear optimization problem.

Note that U is compact, and the recourse is a continuous function in d and s. Under standard assumptions that µ d l,t belongs to the interior of set { U d l,t dQ : Q is a probability distribution over U } for all l ∈ [L] and t ∈ [T ] and that s l,t belongs to the interior of set { U s l,t dQ : Q is a probability distribution over U } for all l ∈ [L] and t ∈ [T ] , slater condition will hold (Boyd et al., 2004) . Then, strong duality holds and the worst-case expectation equals:

where α l,t , β l,t and θ are the dual variables associated with constraints (B.1b), (B.1c), and (B.1d). Note that we can rewrite constraint (B.2b) as follows:

Given that we are minimizing θ in (B.2a), we can rewrite formulation (B.2) as follows:

(d l,t α l,t + s l,t β l,t )

Proof. Not that problem Q(x , y , ξ) in (8) is separable by each t and l. Therefore, w.l.o., we can rewrite Q(x , y , ξ) = l t Q l,t (x , y , d l,t , s l,t ), where for fixed x x x, y y y, l, t, d l,t ∈ [d l,t , d l,t ], and s l,t ∈ [s l,t , s l,t ].

h l,t (d l,t , s l,t , y l,t , α l,t , β l,t ),

where h l,t (d l,t , s l,t , y l,t , α l,t , β l,t ) := max

It is straightforward to verify that for fixed x x x, y y y, d l,t ∈ [d l,t , d l,t ], and s l,t ∈ [s l,t , s l,t ] function (C.3) is convex in variable ρ l,t . Hence, problem (C.3) is a convex maximization problem. It follows from the fundamental convex analysis (see Boyd et al. (2004) ) that here exists an optimal solutionρ l,t at one of the extreme points of polyhedron P l,t := {−c o l,t ≤ ρ l,t ≤ c u l,t }. In any extreme point, constraint −c o l,t ≤ ρ l,t ≤ c u l,t is binding at either the lower bound or upper bound. Thus, given that d l,t ∈ [d l,t , d l,t ], and s l,t ∈ [s l,t , s l,t ], it is easy to verify that the optimal value η * l,t of problem max

h l,t (d l,t , s l,t , y l,t , α l,t , β l,t ) in (C.2) equals to:

h l,t (d l,t , s l,t , y, α l,t , β l,t ) is equivalent to the following minimization problem.

η l,t ≥ d l,t s l,t c u l,t − k:l∈R k y k,l,t c u l,t − d l,t α l,t − s l,t β l,t , (C.5e) η l,t ≥ d l,t s l,t c u l,t − k:l∈R k y k,l,t c u l,t − d l,t α l,t − s l,t β l,t , (C.5f) η l,t ≥ d l,t s l,t c u l,t − k:l∈R k y k,l,t c u l,t − d l,t α l,t − s l,t β l,t , (C.5g)

Replacing max h(·) in C.2 with its equivalent reformulation in (C.5) and summing over all l ∈ [L] and t ∈ [T ], we derive the following equivalent reformulation.

Proof. Not that problem Q A (x , ξ) in (15) is separable by each t. Therefore, w.l.o., we can

3) is not solvable directly due to the trilinear terms between the dual variables ρ t , d t , and s t . Note that we could first perform maximization over d t and s t , i.e.,

The above problem is separable in l. Thus, it is easy to verify that the optimal d l,t is attained at either d l,t = d l,t or d l,t = d l,t , and the optimal s l,t is attained at either s l,t = s l,t or s l,t = s l,t . Accordingly, we define binary variables g l,t = 1 if d l,t = d l,t and g l,

And for all ∈ [L], we define binary variables z l,t = 1 if s l,t = s l,t and z l,t = 0 if s l,t = s l,t . Then, we have:

With (D.5) and (D.6), we drive the following equivalent reformulation of h t (x, d t , s t , α t , β t ) as (D.7). max ρt,λt, gt,zt

∆s l,t z l,t β l,t + d l,t α l,t + ∆d l,t g l,t α l,t + s l,t β l,t + L l=1 d l,t s l,t ρ l,t + ∆s l,t d l,t z l,t ρ l,t + ∆d l,t s l,t g l,t ρ l,t + ∆d l,t ∆s l,t g l,t z l,t ρ l,t (D.7a)

where ∆d l,t = d l,t − d l,t and ∆s l,t = s l,t − s l,t .

Then, we define variables v l,t = g l,t z l . And we introduce the following McCormick inequal-

Accordingly, we can derive the following mixed integer quadratic programming (MIQP) reformulation of (D. x k h k λ k,t − L l=1 d l,t α l,t + ∆d l,t g l,t α l,t + s l,t β l,t + ∆s l,t z l,t β l,t + L l=1 d l,t s l,t ρ l,t + ∆s l,t d l,t z l,t ρ l,t + ∆d l,t s l,t g l,t ρ l,t + ∆d l,t ∆s l,t v l,t ρ l,t In formulation D.9, there are three bilinear terms in objective function (D.9a). To linearize the problem, we introduce three decision variables e l,t = z l,t ρ l,t , q l,t = g l,t ρ l,t , r l,t = v l,t ρ l,t and the following McCormick inequalities for them. 

{ d l,t α l,t + ∆d l,t g l,t α l,t + s l,t β l,t + ∆s l,t z l,t β l,t } + L l=1 d l,t s l,t ρ l,t + ∆s l,t d l,t r l,t + ∆d l,t s l,t q l,t + ∆d l,t ∆s l,t e l,t (D.11a) x k h k λ k,t − L l=1 d l,t α l,t + ∆d l,t g l,t α l,t + s l,t β l,t + ∆s l,t z l,t β l,t + L l=1 d l,t s l,t ρ l,t + ∆s l,t d l,t r l,t + ∆d l,t s l,t q l,t + ∆d l,t ∆s l,t e l,t (D.12a) c

In formulation (G.1) and (G.2), we replace all scenario-dependent parameters, variables and constraints with scenario index n ∈ [N ]. Parameters d l,t and s l,t are replaced by d n l,t and s n l,t to represent the realized services demands and duration in scenario n ∈ [N ]. And we use variables y n k,l,t , o n l,t , and u n l,t to represent the second stage decision variables under scenarios n ∈ [N ].

According to Law of Large Number, the optimal objective values of formulation (G.1) will converge to the optimal values of formulation (1) with probability one as N → ∞ (Homem- de Mello & Bayraksan, 2014; Shehadeh et al., 2021) . However, the computational effort and solution time of the SAA formulations will increase when we increase the sample size. Hence, we use Algorithm 2 to determine an appropriate sample size N and obtain near-optimal solutions of our model via SAA formulation.

In Algorithm 2, we start with a small initial sample size. In first step, we will solve the SAA formulation with designated sample size. In step 1.1, we generate N random samples of demands and service durations. Then, in step 1.2, we will solve the SAA formulation with samples generated in step 1.1 and record its optimal objective value v k N and solution x k N . In step 1.3, we generate N random samples of demands and service duration. We evaluate the solution x k N by Monte Carlo simulation with the generated N random samples and record the objective value v k N . We will run step 1 with K replicates. All results are recorded. In second step, we calculate the average of v k N and v k N as v N and v N . According to Shehadeh et al. (2021) , v N and v N are statistical lower and upper bound of the optimal value of original formulation. In step 3, we calculate the approximate optimality index as

AOI N approximately estimate the optimality gap of lower bound and upper bound of original formulation. Therefore, if AOI N is smaller than the predetermined termination tolerance , we believe the optimality gap is small enough. The algorithm will terminate and output the termination sample size N , objective value v N , v N and approximate optimality index AOI N .

Otherwise, we will increase the sample size then go to step 1.

In Appendix I and Appendix J, we present the result of Approximate Optimality Index and 95% Confidence Interval of the statistical lower bound and upper bound on the objective value, termination sample size and CPU time for 15 random generated instances with lognormal, normal and uniform distributions.

Algorithm 2 Monte Carlo optimization (MCO) Method 1: Initial sample size N 0 (10), the number of replicates K(10), the number of scenarios in the monte carlo simulation step, N (1000), and a termination tolerance (0.01). Step 1: MCO Procedure;

for k = 1, ..., K do 9:

Step 1.1 Scenario Generation 10:

Step 1.2 Solving the SAA formulation 12:

Solve the SAA formulation with the scenarios generated in Step 1.1 and record the corresponding optimal objective value v k N and optimal staffing number (x) k N .

Step Step 2:

Step 3: Compute the Approximate Optimality In this section, we show the convergence results of SAA formulation (G.2). Without spacial mention, we will use the termination MIP Gap 0.1 % (We allowed all instances to run for one day, however, the MIP Gaps of some instances remained at around 0.1 %). And if there are more than 50000 iterations with no MIP Gap improvement, the Gurobi will also terminate. Appendix M. Results of optimal staffing patterns for FA models 

Robust optimization

Theory and applications of robust optimization

Convex optimization

Addressing consistency and demand uncertainty in the home care planning problem

Demand uncertainty in robust home care optimization

A cardinality-constrained robust model for the assignment problem in home care services

Designing appointment scheduling systems for ambulatory care services. Health care management science

Canada's elder care crisis: Addressing the doubling demand

Routing and scheduling in home health care: A literature survey and bibliometric analysis

Data-driven distributionally robust optimization using the wasserstein metric: Performance guarantees and tractable reformulations

Home health care routing and scheduling: A review

Recent advances in robust optimization: An overview

Cost of care survey

Distributionally robust optimization and its tractable approximations

A practical guide to robust optimization

Operational research applied to decisions in home health care: A systematic literature review

Home health care logistics management problems: A critical review of models and methods

Data-driven distributionally robust appointment scheduling over wasserstein balls

The sample average approximation method for stochastic discrete optimization

Robust nurse-to-patient assignment in home care services to minimize overtimes under continuity of care

Monte carlo bounding techniques for determining solution quality in stochastic programs

The home health care routing and scheduling problem with interdependent services. Health care management science

Monte carlo sampling-based methods for stochastic optimization

Scheduling: theory, algorithms, and systems

Care is essential

Distributionally robust optimization: A review

Home healthcare integrated staffing and scheduling

Staff dimensioning in homecare services with uncertain demands

Home health-care network design: Location and configuration of home health-care centers

Reducing conservatism in robust optimization

Lectures on Stochastic Programming: Modeling and Theory

Distributionally robust optimization approaches for a stochastic mobile facility routing and scheduling problem

Using stochastic programming to solve an outpatient appointment scheduling problem with random service and arrival times

A robust optimization for a home health care routing and scheduling problem with consideration of uncertain travel and service times

The optimizer's curse: Skepticism and postdecision surprise in decision analysis

A note on issues of over-conservatism in robust optimization with cost uncertainty. Optimization

Stochastic optimization models for a home service routing and appointment scheduling problem with random travel and service times

From data to decisions: Distributionally robust optimization is optimal

Distributionally robust hub location

Washington state plan on aging

Solving two-stage robust optimization problems using a column-andconstraint generation method

Vehicle routing and appointment scheduling with team assignment for home services

Home service routing and appointment scheduling with stochastic service times

We want to thank all of our colleagues who have contributed significantly to the related literature. Dr. Karmel S. Shehadeh dedicates her effort in this paper to every little dreamer in the whole world who has a dream so big and so exciting. Believe in your dreams and do whatever it takes to achieve them-the best is yet to come for you.

Proof. According to (15b), we have ρ l,t ≤ c k,l,t − λ k,t , ∀k ∈ [K], l ∈ R k . Hence, we have (E.1). ρ l,t ≤ min k:l∈R k (c k,l,t − λ k,t ) (E.1)As we consider a maximization problem in (15) and d l,t s l,t > 0 always holds, it is always optimal to increase ρ l,t in the feasible region. Hence, we have ρ * l,t = min k:l∈R k (c k,l,t − λ k,t ) or ρ * l,t = c u l,t . Then, according to constraints (15d), we have c k,l,t − λ k,t ≥ c k,l,t − c o k,t . This leads to ρ * l,t = min k:l∈R k (c k,l,t − λ k,t ) ≥ min k:l∈R k (c k,l,t − c o k,t ). As we assume c u l,t + c o k,t > c k,l,t , we have min k:l∈R k (c k,l,t − c o k,t ) ≤ c u l,t . Finally. we have that ρ l,t ≥ min k:l∈R k (c k,l,t − c o k,t ) in the optimal solution.Appendix F. Proof of Proposition 5Proof. Observe from the objective of problem (16) that variable α l,t are multiplied by parameters µ d l,t and variables d l,t , for all l ∈ [L], t ∈ [T ]. And so, for fixed β and x, the joint contribution α and d to the objective of problem (16) equals:where K is a non-negative term that equal to:First, we suppose that α l,t > 0. If α l,t > s l,t c u l,t where c u l,t is the upper bound of ρ l,t . The optimal solution of d * l,t will be d l,t . In this case, α l,t contributes to the objective value of problem (16) i.e., α l,t improves the objective value of problem (16). It follows that without loss of optimality α l,t = s l,t c u l,t is valid upper bound on variables α l,t . Second, we suppose that α l,t < 0. If α l,t < −s l,t |ρ l,t | which leads to s l,t ρ l,t − α l,t ≥ 0, the optimal solution of d * l,t always is d l,t . In this case, α l,t contributes to the objective value of problem (16) i.e., α l,t improves the objective value of problem (16). It follows that without loss of optimality α l,t = −s l,t |ρ l,t | is valid lower bound on variables α l,t .Follow the same proof process, we can prove the β l,t and β l,t , defined in (24), are respectively valid upper and lower bounds on variable β l,t . In this section, we show the convergence results of SAA formulation (G.1). without spacial mention, we will use the termination MIP Gap 1 % (We allowed all instances to run for one day, however, the MIP Gaps of some instances remained at around 1 %). And if there are more than 50000 iterations with no MIP Gap improvement, the Gurobi will also terminate. Appendix K. Results of optimal staffing patterns for EA models