Robust meal delivery service for the elderly: a case study in Hong Kong

Service providers in a community center in Hong Kong deliver meals to community-dwelling elderly, first from a central kitchen to intermediate depots by a van and then to the homes of the elderly via walking. We propose a modified two-echelon vehicle routing model with concerns of both delivery efficiency and workload fairness among workers, incorporating important practical aspects, such as continuity of care and unique features of buildings and served elderly. Notably, we employ robust optimization to address service time uncertainties that differentiate between frail and ordinary elderly. The robust model can be transformed into a mixed integer program, for which we provide two decomposition-based approaches to accelerate computation. Through a real-data case study, we verify the effectiveness of the proposed models. We show that robust solutions can protect against service time variations and achieve better performance while incurring a small additional cost over deterministic ones. We provide insights into choosing the level of conservatism and human resource planning for practitioners.

are high-rise, and several clients may live in the same building. In practice, home care services delivered by an IHCST from a center to clients are made by two types of transportation, i.e., first by van and then on foot. In this work, we model such kind of service deliveries as a modified two-echelon vehicle routing problem (2E-VRP). In classic 2E-VRPs in the literature, the freight delivery from the depot to the customers is managed by shipping the freight through intermediate depots 7 . In terms of the home care services in HK, the transportation network is decomposed into two levels, where the 1st level connects a center to the intermediate depots (say van stops or other types of drop-off locations) and the 2nd connects the intermediate depots to the buildings of clients. Without loss of generality, we will use meal delivery as a case to illustrate home care services delivered by IHCSTs in the remaining context, as shown in Fig. 1.
To assign delivery tasks and plan delivery routing for workers on a complex two-echelon network is challenging (as shown in Fig. 1). The current practice of IHCSTs still relies on a manual approach to generate task assignments and routing solutions, which is experience-based, with resulting solutions always far from optimal. It is highly demanding to propose a scientific approach to it. Based on consultations with a center manager, we note that working times spent on deliveries variate among workers mainly due to uncertain service times, which has been one of the main reasons for high turnover rates because of uneven and over workload, causing the feeling of unfairness among workers. In addition, as the workers may share the same van to return to the delivery center after finishing the assigned service tasks, ensuring a fair workload is also for efficiency. Thus, both center managers and workers find fairness to be very important. In this study, we use a min-max measure over working times spent by different workers (i.e., to minimize the "makespan"). According to historical observations in practice, service times largely variate case by case, and the fluctuations differentiate between ordinary cases and frail cases. A frail elderly, as opposed to an ordinary elderly, typically requires more time when receiving the same service, and service time for them is more likely to show larger variations. For the case of meal delivery, uncertain service time may lead to scenarios in which the realized workloads of some workers are extremely long, and the time of some clients to receive the meal has been out of the proper lunch or dinner time window. The abovementioned challenges have been a huge burden on already cash-strapped IHCSTs. In this study, our objective is to provide a scientific approach for robust route design and task assignments against uncertain service times.

Related work
The existing studies on home-care operations focus on vehicle routing and scheduling to perform various services at clients' homes/sites using operations-research modeling (see reviews 3,4 ). Utilizing mathematical modeling and computing methods has shown an increased operational efficiency by 10-15% 8 , and savings ranging from 5 to 20% of the total costs 9 , compared with manual approaches. The home care service problems in the literature are formulated as a classic vehicle routing problem (VRP) 3,10 . Differently, due to the mentioned unique features of buildings in HK, we model the meal delivery service problem as a modified 2E-VRP. The two-echelon version of the VRP (i.e., 2E-VRP) was considered in freight distribution systems for congested urban areas 11,12 and was first formally studied in 2011 7 , where the delivery from the depot (i.e., delivery center) to the customers is managed by rerouting and consolidating the freight through different intermediate satellites. In our meal-delivery problem where van stops correspond to the intermediate satellites (as shown in Fig. 1), the focal setting has some salient differences compared with classic 2E-VRPs. First, the demand at the 2nd-level served points in classic 2E-VRPs cannot be split among different service providers 7 . On the contrary, the same building of clients can be served by more than one worker in our setting, e.g., buildings 1, 7, and 8 in Fig. 1. Thus, the aggregated demand in a building should be split into client-based points in our model. Second, in our meal-delivery case, the delivery  www.nature.com/scientificreports/ is allowed to go directly from the delivery center to the clients, e.g., the routing of 0-5-4 in Fig. 1. For classic 2E-VRPs, instead, the freight must be consolidated from the depot to a satellite and then delivered from the satellite to the end customers. Third, the existing studies on classic 2E-VRPs consider deterministic settings or depend on a two-stage stochastic setting for tactical planning 13 . We focus on 2E-VRPs under uncertainty, which has been largely ignored.
To address the service time uncertainty and protect against the extremely long working time for some workers, we utilize robust optimization techniques in our modified 2E-VRP incorporating uncertainty. Robust optimization has been widely applied to cases with partial knowledge of the underlying distribution 14 . Various studies have developed effective reformulations and algorithms for generic optimization problems with uncertain parameter sets [15][16][17][18] , with applications in many different areas including transportation 4 , engineering 19 , and healthcare setting 20 . In our meal-delivery case, we specify two uncertain-but-bounded perturbation sets separately for service times of the ordinary and frail elderly, where parameters of the uncertainty set are obtained based on historical observations and consultation. For the frail elderly, due to the different magnitudes of mobility limitations (e.g., meal-on-wheel 21 ), the service time is relatively long with a greater perturbation range compared with the ordinary ones. We employ a budget-constrained robust model 22 , which is featured with a budget on uncertainty to control the level of conservatism (i.e., to control the percentage of served clients' service time that can take their extreme duration values). The resulting solution, subject to the budget constraint, will become safe and perform well even under highly uncertain scenarios. In this study, we separately specify the budgets on the service time uncertainty for the ordinary and frail elderly.
Our proposed robust optimization model can be applied to routine community-based home care services in densely populated and aging cities (e.g., HK). The contributions of this work are summarized below.
• Due to the unique features of residential estates in HK, we first model the home care services as a modified two-echelon vehicle routing problem (2E-VRP), incorporating concerns of both efficiency and fairness under a meal-delivery case. As to be introduced in the next section, our model also incorporates other important practical aspects, such as "continuity of care" and features of buildings (e.g., setup time, worker capacity, with or without a lift). • We propose a robust 2E-VRP model under service time uncertainty that is differentiated between frail and ordinary elderly. We show that the robust model can be transformed into a mixed integer program (MIP), for which we provide two decomposition-based approaches to ensure computation efficiency. • Through a real-data meal-delivery case study, we verify the effectiveness of our modified 2E-VRP model and its robust version, with improving on average 4.7-7.0% delivery efficiency compared with the operations implemented by IHCSTs. The robust solutions can protect against service time variations and achieve better performance while incurring a small additional cost over deterministic optimal solutions. We also provide managerial insights regarding risk management and long-term human resource planning.
The remainder of the paper is organized as follows. First, the problem is described, and formulations (deterministic and robust counterparts) are presented. Second, two decomposition-based algorithms are developed for the computation of the resulting MIP. Third, the numerical results of the case study and discussions are presented. Last, conclusions and directions for future work are discussed.

Model formulations
In this section, the model development for meal-delivery planning is presented. First, we define our nominal (deterministic) model, followed by its robust counterpart and its linear equivalent.
Problem description. We use a network G = (N, A) to describe our problem where N is a node set and A is an arc set. We define K as the set of meal-delivery workers and B as the set of buildings that need to be served. In our meal-delivery network, N = {0, e * } ∪ B ∪ S , i.e., the node set is comprised of the delivery center (node 0), the dummy ending point e * , a set of served buildings B, and a set of van stops S that can be traversed. Van stops are intermediate parking areas or stations pre-specified for delivery vans, similar to bus stops. Let L be the full client set and L i be the set of clients to be served at building i ∈ B . The number of clients to be served in building i ∈ B is |L i | . Before delivery, the information about clients from all different buildings is known, and the demanding services must be fulfilled. All workers (i.e., ∀k ∈ K ) start at the same time from the same delivery center (node 0). As mentioned, clients' buildings are reached through a 2-echelon delivery: first via van and then on foot. We have the arc set A = {(i, j) : i, j ∈ N, i � = j} , where the related distance for each arc is expressed in travel time τ ij . We refer A S ( A S ⊂ A ) to the set of arcs linking van stops and the delivery center. Specifically, for (i, j) ∈ A S , we refer τ ij to the travel time by van, and for arcs not belonging to A S , we refer τ ij to the travel time on foot. The route via van is comprised of some linkages among |S| van stops in sequence. For each delivery (i.e., lunch or dinner), each worker is assumed to get off at no more than one van stop. In other words, once a worker gets off the van at a certain stop, he/she cannot go back onto the van before he/she finishes the assigned delivery jobs, i.e., visiting and serving his/her assigned clients on foot sequentially. Note that workers can also directly go from the delivery center to reach clients on foot. This permission largely expands the solution space of the 2-echelon setting, bringing challenges in computation. After finishing all of the assigned jobs, the workers will go back to the delivery center. The time spent on the back journey (i.e., τ ie * ) is calculated by the shortest time spent from node i to the delivery center (on foot and by van). As customers may share the same van and thus need to wait for the ones who are relatively late to arrive at van stops, without further justification, we use a dummy point e * instead for simplicity. www.nature.com/scientificreports/ The objective of task assignments and route design is to maximize delivery efficiency while ensuring the fairness of working time among workers. The total working time for each worker includes the time spent on the road, entering a building, and serving frail and ordinary elderly. The delivery center promises the clients to obtain their meal at a regular 1-hour time interval (e.g., 11:30 -12:30 for lunch and 17:30-18:30 for dinner). Although the objective of the current study is worker-orient, improving the delivery efficiency is also to ensure that clients can receive the meals on time. Workload fairness among workers is important in many planning, scheduling, and resource allocation problems 23 . When route design and task assignments are done through mathematical optimization rather than by hand, perceived political issues can be avoided, since workers no longer doubt the fairness of the method used to allocate workload 24 . A min-max fairness approach can balance the workload among workers in home-care service problems 25 . In our model, we adopt such a min-max measure to balance the fairness of workload among workers (to be introduced in formula (1)). Besides, we consider continuity of care, i.e., "having the same care provider who knows and follows the patient ′′ , which is important for home-care service [26][27][28][29] because it fosters improved communication, trust, and a sustained sense of responsibility 30 . Continuity of care can also be reflected in a caregiver-client relationship where there is a preference for each caregiver serving the same, or similar, set of clients 27 . In our model, we assume that each worker has a prespecified set of his/her familiar clients, and additional penalty time is incurred when a worker serves a client he/she is unfamiliar with. In the context of meal delivery, no penalty reflects the worker's ability to quickly and easily find his/her way to a client's living place. With such a consideration, the worker becomes "heterogeneous" in terms of the delivery time for familiar or unfamiliar clients, which brings new challenges in modeling and computation.
The set notations, parameters, and decision variables used throughout the paper are summarized in Tables 1, 2 and 3. We first have a deterministic model (DM) where the service time is assumed to be certain.    www.nature.com/scientificreports/  (1) is to minimize the maximum working time spent among workers. The total working time for each worker k includes the time duration on the road (i.e., time spent on van and on foot), setup time to enter buildings, and service time for the frail elderly and ordinary elderly. Constraints (2) and (4) ensure that every worker starts at the delivery center (i.e., node 0) and finishes at the dummy node. Constraints (3) and (5) take care of flow conservation for each worker k in nodes and arcs, respectively. Constraint (6) ensures that one worker should take on one van to traverse between van stops. Constraints (7)-(9) ensure the network flow balance of vans. Constraint (10) is the capacity for each van. Constraint (11) ensures that all demand is fulfilled, and constraint (12) ensures that visiting a building without filling any demand is not allowed. Constraint (13)-(15) are the capacity constraints on the total demand filled by each worker, demand filled by workers serving a building without a lift, and the number of workers to enter into the same building, respectively. Constraint (15) is called demand splitting where α is the parameter to control the magnitude of splitting, i.e., at most ⌈|L i |/(αC)⌉ workers can serve building i. This splitting constraint is just following the practice. Constraint (16) is for sub-tour elimination. As described above, the DM model balances the delivery efficiency and fairness of workloads among workers, incorporating some important practical aspects, including building-specific features (e.g., setup time to enter into a building, workload capacity without a lift, demand splitting) and client-specific features (e.g., "continuity of care", frail and ordinary cases). The DM model can also be reformulated as the following MIP formulation.
Model formulation under uncertainty. As introduced, the service times for clients are highly unpredictable with great variations, leading to uneven realized workloads among workers and further exacerbating delivery efficiency and staff shortage due to the high turnover rate. We model the uncertainty in service times for the ordinary and frail elderly (i.e., µ o l and µ f l ) separately as serving a frail elderly takes a relatively long time with a greater fluctuation range observed from historical practice. We employ box perturbation sets to model the service time uncertainty with µ k l ∈ [µ k l −μ k l , µ k l +μ k l ] , where µ k l is the nominal value and μ k l is the positive radius of perturbation set 18  To control the level of conservatism, we further employ a budget-constrained approach that retains the advantages of computationally demanding linear programming 18 . We use parameters Ŵ k o and Ŵ k f to control the total deviation of the uncertain parameters (i.e., the budget of uncertainty) for the ordinary and the frail elderly, respectively. The goal of the budget-constrained approach is to protect against all cases in which up to ⌊Ŵ k o ⌋ coefficients (i.e., µ k l ) of ordinary service jobs are allowed to change by μ o and one coefficient changes by at most and Ŵ k f is likewise for the frail service jobs. We then have the following robust optimization model (RO): www.nature.com/scientificreports/ The budget constraints (28) and (29) in the RO are to model the magnitude of conservatism against uncertainty, with β k o for the ordinary case and β k f for the frail case, respectively. Each of them involves an inner optimization problem. Take the inner optimization problem in (28) as an example, S k is a decision variable set ( S k ⊂ L o ) corresponding to selecting up to ⌊Ŵ k o ⌋ coefficients and t k is the decision variable to select one coefficient that changes by at most (Ŵ k o − ⌊Ŵ k o ⌋)μ o . Next, we show that the RO formulation has an equivalent MIP. Given y , problem (28) can be transformed 22 as which is a linear program and the associated dual is: where w k l and r k o are the dual variables. By strong dual theorem, objective values of the original problem and dual problem coincide. In the same way, we have β k f = min r k ∀l ∈ L f } and its dual format. Therefore, we have the following equivalent MIP formulation of the RO model: The MIP formulation (35)-(41) can be directly solved by some commercial optimizer (e.g., CPLEX or Gurobi). For large-scale instances, it is even time-consuming to find a feasible solution.
Two decomposition-based approaches. The existing studies on algorithms for the deterministic 2E-VRP either consider metaheuristics such as adaptive large neighborhood search algorithm 31 or column generation-based algorithms 32 . These approaches cannot be directly applied to our setting as we consider many important and practical features in our modified 2E-VRP model under uncertainty. We propose two decomposition-based approaches to accelerate the computation of the MIP.
First, we consider an easy-to-implement decomposition-based heuristic (DH) to reduce the complexity of the problem. We first decompose the problem into two segments: buildings directly served by workers on foot and buildings served by vans then on foot, according to the physical distances from the delivery center to the buildings. Such a segmentation strategy has been shown to be effective in practice 33 . After decomposing the service region into two segments (say on-foot segment and by-van segment), one needs to further assign workers to two segments. In our heuristic, we first determine the number of assigned workers according to the total worker-client percentage and then determine the on-foot workers based on the familiar set of workers. The worker with more familiar clients living in the on-foot segment is more likely to be selected as an on-foot worker. This strategy can largely reduce the solution space, making the resulting two subproblems solvable within adequate time. such a segmentation arises naturally since there is already a distinction between buildings served on foot and by van in practice. This practice-driven decomposition heuristic can also be extended to solving highly large-scale instances by decomposing the problem into more segments. (30) (2)− (21). www.nature.com/scientificreports/ Second, through utilizing the structure of the MIP formulation, we consider a Benders decomposition (BP) approach 34,35 as the BP technique can well separate discrete variables of routing and assignment decisions and continuous variables from the transformed RO model. We briefly describe the BP framework for our mealdelivery context. In the MIP, given routing variables x and assignment variables ȳ and z , the resulting subproblem becomes the following linear programming: Let ξ and ζ be the dual variable associated with constraints (43)-(45). The dual of the SP formulation is given by: By strong duality, the optimal objective value of the primal problem (i.e., SP) is equal to that of its dual problem (i.e., DSP). Combining the resulting solution (w * ,r * ,T * ) with (x,ȳ,z) , we obtain a feasible solution, providing an upper bound of the optimal objective value for the MIP. Note that the feasible region defined by (48)-(53) does not depend on the value of (x, y, z) . It is also straightforward to check that the optimal values of the SP and DSP are finite subject to non-empty region. Let ξ j (j ∈ J ) denote all extreme points of the feasible region in the DSP, with [ξ j ] k and [ζ j ] k l be the associated values. The DSP then can be reformulated as: Therefore, the MIP can be reformulated as a master problem (MP) of the BP approach: where constraint (56) is also called Benders cuts. It may be impractical to enumerate all extreme points. One can start from the MP with some bender-cut constraints associated with a subset of extreme points, and further solve a dual subproblem to add new Benders cuts to the MP. Various alternative BP strategies have been introduced (as reviewed in the paper 35 ). We adopt the classic one to accelerate the computation as shown in Algorithm 1. As the www.nature.com/scientificreports/ DH can solve the problem within minutes, we utilize it to find the initial feasible solution when implementing the BP algorithm in our numerical test.

Case study
We carry out a real-world case study of meal delivery by the IHCST in HK. First, we utilize real data to test the effectiveness of our modified 2E-VRP model, compared with the current manual approach. Second, we test how our robust model protects against uncertainty, i.e., how the delivery efficiency changes with the chosen level of robustness. Third, we illustrate and provide insight into how large is the price of robustness in our meal-delivery case. Finally, we discuss some perspectives and insights in terms of future demand growth. All the computational experiments are coded in Python programming language 36 and the Gurobi 7.5.1 solver is used on a Windows 7 PC with an Intel Core i7 processor, 3.60 GHz CPU, and 16GB RAM. An IHCST serves elderly clients living in neighboring communities in HK, with daily meal deliveries, including lunch and dinner. Some clients ask for both lunch and dinner and some ask for either. We collected related data and information (in 2017-2018) from the following sources: (i) consultations with delivery center managers, (ii) field observations from regular delivery sessions that provide time spent on routes and meal deliveries, (iii) maps that cover the delivery routes, van/bus stops, and buildings of clients, and (iv) a database of the clients, their living units and regular delivery staff of the IHCST. Travel time and service time duration for the frail and ordinary elderly were measured and recorded by a team who followed the delivery workers and marked the routes. The task assignment and routing by the IHCST were generated manually largely based on workers' experience. We observed that the nearby buildings were served by workers directly on foot, while buildings farther away were reached by workers who first take via a delivery van and then on foot from the stops to their destinations. Compared with the manual approach in practice, we utilize the available data and measurements to test the effectiveness of our modified 2E-VRP models (i.e., the DM and RO). To ensure a fair comparison, we simulate with 5,000 replications that ordinary and frail elderly randomly realize to take their extreme values in terms of service time, and compare the average realized working time spent by workers under different approaches. Except for the uncertain service time, the daily elements of planning (e.g., the traveling time for each link and aggregate demand) are relatively stable in a short time. Hence, we focus on a typical day's planning in December 2017 without losing generality. In total, the IHCST serves around 70 elderly clients spreading across 17 different buildings. Following the practitioners, the demand splitting parameter α is set as 0.5 and the additional penalty time is 1 min. In terms of the uncertain service time, we obtain the average service time for the ordinary and frail elderly (i.e., μ o l and μ f l ) based on the historical observations. As the service time is with highly unstable perturbations, we consider two sets of radius parameters for the test: a mild perturbation with μ o = 80 and μ f = 160 and a more severe perturbation with μ o = 120 and μ f = 200 (expressed in seconds). These perturbations were discussed and estimated together with the center's staff and were taken as base cases. Take the first set of parameters as an example, we use a perturbation set of Comparison of results. The results associated with two perturbation sets are shown in Table 4 (with μ o = 80 and μ f = 160 ) and Table 5 (with μ o = 120 and μ f = 200 ). The first column distinguishes between two delivery sessions (lunch and dinner). The second column shows the percentage of ordinary and frail elderly (i.e., % o and % f ) taking their extreme values in the simulation. We have . The next three columns show the realized objective values averaged over 5000 simulation runs of three approaches (i.e., manual approach, our DM, and RO), followed by three columns showing percentage change (i.e., Gap M,DM , Gap M,RO and Gap DM,RO ) between the three approaches with "M ′′ denoting the manual approach. The manual approach is adopted by IHCSTs to generate task assignments and routing solutions based on experience. The bottom row shows grand averages across all shown realizations. As shown in Table 4 and 5, the average realized time spent by workers is about 3.0-4.5 min less for the DM and RO compared to the manual solution, with 4.7-7.0% improvement. Specifically, the improvement can be as large as 14.1% (for the case of % o , % f = 1.0 ). When % o , % f = 0.0 (i.e., no one takes the extreme value), the RO is equivalent to the DM, with a 4.6% improvement compared with the manual approach. We first conclude that our modified 2E-VRP model is effective in modeling meal delivery for the elderly in HK. As expected, the benefit of the RO becomes more evident when a higher percentage of clients take their extreme values, showing larger improvements over the DM. The DM outperforms the RO when the proportion of extreme values is small. This illustrates the price of robustness (to be discussed below). When the service time is subject to more severe perturbations (as shown in Table 5), the improvement of the RO is more significant, with an average of 7% improvement compared with the current manual approach. The improved efficiency is mainly attributable to our practical model capturing important real-world attributes. In practice, the written working time limit on each delivery session www.nature.com/scientificreports/ for each worker is an hour, corresponding to the length of the proper time window for customers to have meals. As shown in Tables 4 and 5, the realized time of the manual approach is very possible to be larger than 60 mins, especially when the service time is subject to a more severe perturbation. This means that some workers are overloaded and may feel an unfair workload, further leading to a high turnover rate and exacerbating the staff shortage issue. On the contrary, the realized time of our RO model is within 60 mins for most of the scenarios. Regarding workload fairness, we also calculate the variance of realized time spent among the delivery workers, as shown in Fig. 2. We observe that by fixing different parameters and taking an average of the variance of the realized time among workers, the variance value through the manual approach can be reduced significantly by introducing our models, especially the RO model. This means that the min-max objective format in our proposed models can well address the fairness issue regarding the workload of the delivery workers. Note that we also tested more instances on other various days and observed similar results as in Tables 4 and 5 and Fig. 2. Therefore, we further confirm the effectiveness of the RO model in protecting against uncertain service time, with a balance between delivery efficiency and workload fairness.

Price of robustness.
Given the effectiveness of introduced robustness in our modified 2E-VRP model, we further analyze the price of robustness in terms of objective value, i.e., the reduction in the objective value from its nominal optimal value to its robust optimal value. As the levels of Ŵ k o and Ŵ k f increase, the solution becomes more conservative and more robust against uncertain service times. Table 6 shows the price of robustness at different levels of chosen conservatism (as shown in the column "Robustness ′′ ), where the "Obj ′′ columns refer to the objective value of the RO. Note that this objective value is the worst-case value, different from the average "realized value ′′ in Tables 4 and 5. As shown in Table 6, the price of robustness ranges from 8.4 to 16.9 min. The benefit of using the RO as opposed to the DM increases with the perturbation scale and with more extreme realizations. This is particularly visible in the last column of Tables 4 and 5. Table 6 demonstrates that the price   www.nature.com/scientificreports/ of robustness does not increase further after a certain level of robustness. Hence, whether we choose to be robust against 40% or against 80% of clients taking their extreme values, the worst-case objective value will be the same, although their realized values may differ. This provides us an insight that the center managers can be much more conservative without paying a high price of robustness when implementing the RO model. The results inform practitioners to be risk-averse toward uncertain service times and uneven workloads. The robustness parameters Ŵ o and Ŵ f control the degree of conservatism that the user wants to exert. In practice, as we don't know the exact magnitude of the permutation set in advance, we need to consider some practical methods to determine the parameters, although the price of robustness is not high. Some standard procedures could be followed to choose a proper level of conservatism against uncertainty 37 . One can consider a validation method by introducing a randomly generated validation set of historical observations to help determine these parameters. For example, following the standard procedure, we can divide the historical data into a training set and a validation set with a ratio of 0.8:0.2.  www.nature.com/scientificreports/ Computation efficiency. In terms of computation efficiency, we report and compare the objective values obtained by our two decomposition-based approaches (i.e., DH and BD) with a CPU time limit set as 15 min and directly by the solver (i.e., Gurobi) with a CPU time limit set as 30 and 60 min in Table 7. We observe that our approaches (DH15 and BD15 in Table 7) can effectively accelerate the computation, obtaining better results within the less CPU time limit compared with being directly solved by Gurobi (i.e., Solver30). The practical planning of routing and assignment is operated daily, and thus the computation time before daily delivery should be controlled within a short time. Notably, the decomposition-based heuristic (i.e., DH) performs so efficiently, obtaining high-quality solutions within just a few minutes.
Sensitivity analysis. We also conduct sensitivity analysis concerning having ± 1 and ± 2 workers. Figure 3 shows the realized objective value for 0,  Demand growth. Demand for home care services is expected to grow, and we assess the effect on our model's objective value and computation efficiency. We have a test session with a demand of 98 clients spread across 20 buildings. In our test, we observe that the RO directly solved by Gurobi cannot achieve the optimal (or even obtain very bad solutions) within the 1-h time limit for all of the instances. With our DH approach, we can find better solutions for the RO within several minutes. Table 8 shows average realized values after 5000 simulations with different levels of conservatism under the perturbation set [80,160]. We observe that the RO (solved by the DH) can achieve more benefit than the DM with on average 5.6% improvement in realized value when the demand becomes higher in comparison with small-scale problems (as shown in Tables 4 and 5). We also observe that the average realized makespan is much shorter than 60 mins, meaning that the human resources can be better utilized under scientific approaches and the positive effect of the modified 2E-VRP modeling is amplified with growing demand. This provides insights into long-term human resource planning for decision-makers. That is, the parameter α in staffing formula i∈B ⌈|L i |/(α × C)⌉ can be modified to a larger value (with 0.5 ≤ α ≤ 1 ) when the demand is growing. The resource aggregation effect under scientific approaches can partially offset challenges faced in demand growth for meal delivery.

Concluding remarks
The current context provides a case study of robust optimization applied to a practical problem of home care services, which is uniquely modeled as a modified 2-echelon vehicle routing problem (2E-VRP). We provide a scientific approach for decision-makers to determine the appropriate routing and workload allocation to improve delivery efficiency, with different levels of conservatism. Specifically, our model and analysis have incorporated important practical aspects, such as fairness among workers, "continuity of care", and features of buildings. Notably, we differentiate between the frail and ordinary elderly in terms of uncertain service time, which is mainly addressed in our focal robust 2E-VRP model. We also provide practical decomposition-based approaches to accelerate the computation of the robust 2E-VRP model, which is transformed into a mixed integer program. Through numerical analysis, we verify the value and effectiveness of our proposed 2E-VRP models in the mealdelivery context, compared with the current practice. The structures of our models can also be applied to general community-based home care services (such as housekeeping and medical examination) in any densely populated and aging city like Hong Kong. Our case study results show that robust solutions can protect against service time variations and realize better performance while incurring a small additional cost over deterministic optimal solutions. As demonstrated, the benefit of using robust optimization also increases as larger perturbations apply. With our scientific approach against uncertainty, the improvement in delivery efficiency can further reduce the huge burden on the already cash-strapped delivery team and better coordinate with resources in integrated care services. The results and analysis also shed light on risk management (e.g., determining the level of conservatism) and long-term human resource planning. One of the future directions is to integrate multiple home care services into a single optimization effort. It would be advantageous if different services, such as housekeeping, medical examination, and meal delivery, can be coordinated and optimized simultaneously 3 . In this case, a van assisting workers in improving on-road transportation between different buildings (instead of walking) can be considered to improve the efficiency of the delivery system. The second direction is to investigate delivery economies where vehicles and staff are tracked in real-time and planning decision epochs are much more dynamic and granular. Considering a dynamically updated delivery system may better protect against real-time variable service times. As real-time decision-making largely depends on the internet of things (e.g., sensors) to obtain real-time information, designing a system that is remotely programmable to implement optimization algorithms within the sensor-based network can also be an interesting direction 38 . Another practical direction is to incorporate the elements of ride-hailing service 39 into the meal-delivery system. A discussion on efficiency and fairness under such a new context would be promising. www.nature.com/scientificreports/

Data availibility
The onsite delivery data and simulation data that support the findings of this study are available at reasonable request from the corresponding author. The data are not publicly available because they contain information that could compromise the privacy of research participants. Due to the nature of privacy issues, supporting data about the detailed information of buildings and clients in this study are not available. www.nature.com/scientificreports/