Home healthcare routing and scheduling of multiple nurses in a dynamic environment

Human resource planning in home healthcare is gaining importance day by day since companies in developed and developing countries face serious nurse and caregiver shortages. In the problem considered in this paper, the decision of patient assignment must be made immediately when the patient request arrives. Once patients have been accepted, they are serviced at the same days, times and by same nurse during their episode of care. The objective is to maximise the number of patient visits for a set of nurses during the planning horizon. We propose a new heuristic based on generating several scenarios which include current schedules of nurses, the new request under consideration, as well as randomly generated future requests to solve three decision problems: first, do we accept the patient? If so, which nurse services the patient? Finally, which days and times are weekly visits of the patient assigned to? We compare our approach with a greedy heuristic from the literature by considering some real-life aspects such as clustered service areas and skill requirements, and empirically demonstrate that it achieves significantly higher average daily visits and shorter travel times compared to the greedy method.


Introduction
Home Healthcare (HHC), also referred to as in-home care, social care, or domiciliary care, is becoming one of the most important components of healthcare. HHC helps hospitals and retirement homes to free capacity and decrease care delivering cost (Hall 2012). The most crucial objective of HHC is to ensure people who need medical attention and daily care receive high-standard home services. According to patients' needs, nurses, physicians, doctors, and operators visit patients' homes periodically and provide services. Many elderly, people who are chronically ill, and individuals with disabilities receive HHC services (CMS 2008).
In 2008, the United States (US) saved $25 billion in hospital cost thanks to HHC services according to the National Association for Home Care and Hospice. HHC companies employed 1.8 million caregivers in 2014 and it is estimated that 500,000 more jobs are potentially created by 2024 (http://www.hcaoa .org/asset s/1/27/Value _of_Home_Care___SECUR ED.pdf). 40% of adults aged 65+ already receive HHC. The majority of HHC users are people with an average age of 69. 70% of all Americans aged 65 and older will need HHC service at some point in their lives (http:// www.hcaoa .org/asset s/1/27/Value _of_Home_Care___SECUR ED.pdf). Due to some factors such as aging population, chronic diseases, insufficient capacity of hospitals, etc., it is projected that the demand for HHC doubles by 2040 (Landers 2010). The following information shows why HHC is steadily gaining more importance in the US: • The number of people aged 65 and over will be 56 million by 2020, numbers reach 84 million by 2040 (http://www.censu s.gov/popul ation /www/proje ction s/ usint erimp roj). • Care of a patient in the home costs only $45,000 per year for an average of 44 h of care per week while $91,250 are spent for a patient receiving care in a nursing home (http://www.hcaoa .org/asset s/1/27/Value _of_Home_Care___SECUR ED.pdf). • Home-based health technologies cost $3 billion in 2007 versus $7.7 billion in 2012 (Hall 2012). • The percentage of American adults who are chronically ill is more than 50% (AHRQ 2007).
In this study, we focus on the acceptance decision of a request as well as the scheduling. In the literature, an acceptance policy is occasionally discussed in different areas such as public transportation (Ritzinger et al. 2016) and vehicle dispatching (Ichoua et al. 2000). The HHC provider has to decide whether or not to accept the patient and, if accepted, assign suitable appointment days, times and a nurse. In each shift, all nurses start from their homes, visit scheduled patients at the agreed appointment times, serve them for the prescribed time, and finally return to their homes. The objective of this problem is to maximise the average number of daily visits performed by nurses. We require that patients are serviced at the same days and times during their service horizon, which is called consistent, periodic Vehicle Routing Problem (VRP). Furthermore, the problem is dynamic, and acceptance and assignment time decisions have to be made as soon as patient requests arrive. Although some papers provide solutions to this problem by proposing greedy algorithms, these algorithms do not consider or only partially consider future demand. We propose a Scenario Based Approach for Multiple Nurses (SBAM) which simulates several scenarios, scheduling randomly generated and actual requests with a simple but fast heuristic, and analysing results to decide whether to accept the new patient and at which appointment day, time and nurse. Demirbilek et al. (2018) already studied the Scenario Based Approach (SBA) for a single nurse and showed that this approach is superior to greedy approaches under different scenarios. While it is possible to turn a multi-nurse problem into several single nurse problems by dividing the service region into several districts, such districting is difficult in practice, and would not be applicable if nurses have different skills. Our main contributions are the following: • A heuristic for the dynamic HHC nurse routing and scheduling problem that anticipates future demand implications and that can handle multiple nurses and different skill levels, • Empirical insight how much better it is to plan nurses' routing and scheduling without restricting nurses to districts. • Insights to algorithmic performance under different conditions such as clustered service areas, different service times and service horizons. • Empirically demonstrated improvement over a benchmark heuristic proposed by Bennett and Erera (2011).
In Sect. 2, we present a literature review on HHC nurse routing and scheduling problems, and Dynamic VRP. In Sect. 3, we define the problem. Our solution approach for the single-and multi-nurse cases are presented in Sect. 4 and 5, respectively. The approaches are empirically evaluated and compared with other approaches in Sect. 6. Finally, we conclude our study and outline future opportunities in Sect. 7.

Literature review
In this section, we go over the most relevant studies in HHC and on the Dynamic Vehicle Routing Problem (DVRP) which is closely related to our problem and solution methodology. Note that some parts of this section are based on Demirbilek et al. (2018).

HHC studies
HHC-related models were first investigated by Begur et al. (1997). They constructed a decision support system for a home healthcare company to optimise their routing and rostering operations without considering time windows. We refer readers to Fikar and Hirsch (2017), Mutingi and Mbohwa (2013), and Cisse et al. (2017) for state-of-the-art reviews of the models and algorithms that have been reported in the HHC nurse routing and scheduling literature, and concentrate on some key papers.
Tables 1, 2, 3 and 4 represent a classification of publications in terms of planning horizon, uncertainty, objectives/performance measures, constraints, and solution methodologies. Redjem and Marcon (2016) present a HHC service problem  Home healthcare routing and scheduling of multiple nurses…  that covers multiple visits to the same patient in a day and temporal dependencies based on ordering some tasks in terms of their relations with each other. They aim to minimise waiting and travel time of caregivers under hard time window constraint. They present a two-stage caregiver routing heuristic: At the first stage, the shortest travel time is found for each caregiver without coordination of patients and sequencing restrictions. At the second stage, all assumptions and constraints are integrated into the final solution. According to their results, their algorithm is very efficient in terms of computational time. However, it does not take temporal dependencies sufficiently into account. Rest and Hirsch (2016) introduce a HHC model which considers routing and scheduling on a daily basis. Their objective is to minimise travel and waiting times. They consider maximum daily working time with shifts, breaks after a certain time work, clients and workers' satisfaction factor, and caregivers' qualification levels as constraints. Public transportation such as bus, train, and metro is taken into account to model travelling of caregivers between their visits and real time tables derived from public transport service are used when calculating time-dependent travel times. They propose three Tabu Search based solution methods for the scheduling problem. Guericke and Suhl (2017) develop a HHC model also considering work regulations and legal requirements. They take into account break-rest times, weekly work durations, and shift rotations according to laws and regulations in Germany to be able to investigate their influence on results. They propose a mixed integer linear solution for a small-size problem setting. Moreover, an adaptive large neighborhood search based heuristic is provided to cope with a real-size complex problem in a reasonable computational time. The heuristic method performs well on run time compared to the mixed integer program. Wirnitzer et al. (2016) develop a monthly nurse rostering model for a HHC company to optimise scheduling activities done manually before. They propose five mixed integer programming formulations. Each one has a different objective function targeting continuity of care from the perspective of patients with same hard constraints such as breaks, maximum daily and weekly working times, patient/nurse preferences, and shift rotations. According to results based on randomly generated data derived from real-world input and the company's data, all models outperform the manual planning in reasonable times. Moreover, they compare the results of models in terms of computational times, the number of assigned and switched nurse. Cappanera and Scutellà (2014) develop a model that matches operators' skill levels to patients' needs. They propose an Integer Linear Programming formulation to solve this assignment problem including scheduling and routing factors. They use real data derived from HHC providers in Italy to evaluate their model and observe that it works successfully on large HHC instances. Based on their work, Yalçındağ et al. (2016) develop a two-phase solution methodology. The main difference between the two studies is that Yalcindag et al. decompose the joint approach of Cappanera and Scutella that includes assignment, scheduling, and routing solutions in order to solve large instances in a computationally efficient way. They test several two-phase combinations in which each phase can cover one or two of assignment, scheduling, and routing solutions. Computational results based on large real world instances suggest that the two-phase method (assignment and scheduling in phase 1;

3
Home healthcare routing and scheduling of multiple nurses… routing in phase 2) is more computationally efficient than the single-phase method by Cappanera and Scutella. Issabakhsh et al. (2018) present a robust mathematical model for patients who need peritoneal dialysis in their home. According to their model, patients have different requirements such as collection of their urine or blood samples, visits by nurses and technicians, and deliveries of some necessary medicines. Because of such necessities, they have to take into account not only depots and patients' locations, but also dialysis centres and laboratories. Furthermore, they consider some constraints such as labs to be visited after collecting blood or urine tests, and nurses and technicians who have to be picked up from a dialysis center before visiting patients. Since just-on-time visits are a crucial factor for peritoneal dialysis, they develop a robust optimisation model to handle uncertainty in travel times. Even under the scenario with the highest uncertainty, the difference between the robust and deterministic result is less than 1.2%. Frifita et al. (2017) develop a model for a HHC problem with time windows and synchronization which means multiple caregivers visit a patient at the same time. They propose a general variable neighbourhood search method to be able to minimise travel times of caregivers. The proposed methodology is compared to a mixed integer programming model and a heuristic method for a variety of real-life instances. Their method is fast compared to the mixed integer programming model and provides results close to the optimal solution. Du et al. (2017) present a HHC scheduling optimization problem with known demands and service capabilities. Their aim is to minimise the total service cost that includes travel, service, and penalty costs by considering qualification matches, time windows, and service priority based on seriousness of the patients' conditions. They develop an integer programming model and propose a genetic algorithm with a local search method in order to solve the problem. They compare results of the proposed solution method with a commercial software using a case study in China. They report that the genetic algorithm with local search method provides fast and reasonable results for their real-life data.
Existing papers in the literature generally focus on static problem settings for which all patient requests are fully known, but requests arrive dynamically during the service horizon in the real world. Additionally, they do not consider any acceptance policy. We have found only the study of Bennett and Erera (2011) and Demirbilek et al. (2018) which consider dynamic patient arrivals. Bennett and Erera (2011) present a rolling horizon myopic planning approach for the single nurse HHC problem. They propose a capacity based insertion heuristic when integrating a new patient request to the existing schedule by considering the nurse's distribution of idle time explicitly. Furthermore, they model the problem as dynamic and periodic with fixed appointment times, which means that patients arrive dynamically and they are assigned to days over a predetermined number of weeks to visit. Their objective is to maximise the number of patients being served by a nurse. However, the proposed distance and capacity heuristics are greedy algorithms which try to minimise insertion costs whenever a new request arrives without considering future requests or only partially. Moreover, these heuristics accept all requests and ignore that to reject a request now can allow to accept more requests in the future. If a request of a patient located far from the tour is rejected, more closer requests in the future can be assigned to the tour. In other words, we spend time serving patients instead of travelling between far locations. Demirbilek et al. (2018) used Scenario Based Approach (SBA) for a dynamic patient set. To be able to mimic future demand, SBA generates scenarios including random requests and the current one, and constructs tours modelling a nurse's daily visits. After that, SBA analyses when/how many times the current request is scheduled to decide acceptance and schedule days and times. They show that their algorithm outperforms the distance and capacity heuristics under similar simulation settings. Both studies focused on only a single nurse case. However, companies employ many nurses and simply dividing the whole service region into some subregions according to the number of nurses does not work properly if we do not consider all nurses and future expected demand simultaneously when deciding the acceptance of a patient and scheduling. Therefore, in this paper, we extend SBA from Demirbilek et al. (2018) to take all nurses and future demand into account at the same time during the decision process. Furthermore, we test our algorithm under several scenarios such as different skill levels, clustered service areas, and different service times and horizons.

DVRP studies
In contrast to the classical VRP, real-world applications often force decision makers to design routing plans online where the visit of the next customer is decided as soon as it becomes available. DVRP studies are rooted in Wilson and Colvin (1977). They employ a greedy insertion heuristic to put dynamically arriving requests into a tour for a single vehicle. Interested readers can find detailed literature reviews on the DVRP in Ritzinger et al. (2016), Thomas (2010), Pillac et al. (2013). Because the DVRP literature is vast, we only discuss some papers whose solution methods are related to our solution methodology. Ichoua et al. (2006) suggest a Tabu Search method to exploit probabilistic knowledge about future request arrivals when assigning new customers now. They propose a waiting strategy that vehicles wait at their current locations based on knowledge about future requests if there is a time gap until the next customer service. Hvattum et al. (2006) propose a multi-stage stochastic programming model and a heuristic method. The heuristic generates scenarios including scheduled visits and random customers sampled from known distributions. Each sample scenario is solved as a deterministic VRP and common features in the sample scenario solutions are employed to construct routes. Bent and Van Hentenryck (2004) model DVRP with time windows and aim to maximise the number of daily visits. They propose a multiple scenario approach based on generating routing plans including both known and future customers. A distinguished plan selected by a consensus function in terms of the smallest travel cost is employed for decision making. The multiple scenario approach is tested against greedy approaches under dynamism varying between 30 and 80%. According to results, the new methodology shows superior performance for DVRP models. Although some similarities exist between the solution methodology of Bent et al. and our approach, there are also some differences. First, we do not have any waiting strategy since accepted requests in our study start to be serviced the week after acceptance. We have to consider multiple weekly visits associated with the same patient while each request in their study is scheduled only one time in the service horizon. Moreover, all weekly visits belonging to a request are taken into consideration at the same time in the scenario generation phase when finding the most suitable request. Next, they use time windows while we use exact times for appointments. Finally, the degree of dynamism in their study varies between 30 and 80%. In our problem, all patient requests are dynamic.

Problem definition
The problem we consider is the multiple nurses HHC routing and scheduling problem in a dynamic environment. We follow the definition in Demirbilek et al. (2018) but extend it to allow for multiple nurses.
Nurses Each patient is visited by one out of a set of nurses in a defined geographic service area. Each working day of a nurse is divided into k equally-spaced time slots to schedule patient visits. A set of possible appointment times for nurse n, K n , can be defined: where b is the earliest time for an appointment and is the time between appointment times. Travel time between patient i and j is denoted by m(g i , g j ) in minutes where g i represents the location of patient i. All travel times are always rounded up to the nearest multiple of a time slot.
Patients Inter-arrival times between patients' requests are exponentially distributed with known parameters over the planning horizon. A patient's request i from location g i contains weekly service frequency f i , episode of care ec i that represents how many weeks patient i needs care, and service duration per visit sd i .
Dynamics The problem is dynamic in that there are many acceptance/rejection decisions during the planning horizon. Thus, the solution depends on our scenarios. At each stage (a request arrives), decisions are whether or not the request is accepted, and if so, which day combination, time slot and nurse it should be assigned to. Patients that cannot be scheduled are rejected. Constraints • Let i and j be two consecutive appointments in the schedule of nurse n, and let g i and g j represent locations of the patients assigned to those appointments. Every route for that day is feasible, if and only if for each two consecutive appointments, i and j. • A visit, representing a duty at a patient's home, has to be carried out as often as determined by its frequency and episode. • Weekly visit days can be selected for each patient.
• Patients, if accepted, must be serviced at same days, times and by same nurses every week during their service horizon. • Patients must be serviced by nurses whose qualification levels match the requirements of patients. • Nurses start tours from their homes and end tours at their homes again within the shift's time window. • Nurses have to visit patients in their scheduled time periods.
Objective The objective is maximisation of patient visits during the planning horizon. This is different from maximisation of the number of patients served since patients need different numbers of visits. If T represents a set of patients accepted over the planning horizon, our objective is:

Solution methods for the single nurse case
In this section, we are going to explain and compare solution approaches proposed by Bennett and Erera (2011) and SBA (Demirbilek et al. 2018) for the single nurse dynamic HHC problem since there are many similarities between the proposed solution methodology and SBA as well as the distance heuristic for a single nurse and multiple nurses. Bennett and Erera (2011) proposed two different approaches, the distance and capacity heuristics, for routing and scheduling of a single nurse. The distance heuristic (DH) assigns visits by minimising insertion costs while the capacity heuristic (CH) additionally tries to allocate enough slack time for longer distances that an additional visit may be included later to break up the journey. Patients are assigned to exact visit times and days in terms of their weekly visit needs and have to be serviced at same time and days in each week during their service horizon. According to an empirical comparison reported in Bennett and Erera (2011), CH is superior to DH in small areas while DH provides higher average daily visits in large areas and under high demands.

Scenario based approach (SBA)
DH only assigns visits of patients between assigned pairs of visits without considering future patients' requests. CH takes future requests into account, but only by making sure that long driving distances have sufficient slack to allow insertion of additional customers. Patients are always inserted if this is feasible. Demirbilek et al. (2018) developed a new solution methodology in order to decide acceptance or rejection of patients as well as days and times their visits take place by explicitly considering future demand. Whenever a new patient arrives, SBA generates a number of scenarios that include all existing patients' visits, the new request under consideration, and randomly generated visits mimicking future demand (see Sects. 5.3.1 and 6.1 for more details). A nurse tour is constructed with those visits for each scenario. Finally, the algorithm checks whether or not the current patient has been scheduled in at least one scenario. If not, the patient is rejected. Otherwise, visit days and times are determined based on the most frequently scheduled days and times over all scenarios. They developed two different versions of SBA, Daily SBA (DSBA) and Weekly SBA (WSBA). Tours constructed based on scenarios for each day do not consider weekly visits of each randomly generated request in DSBA. In other words, different requests are generated for each day even though these requests, mimicking future patients, need one, two, or three weekly visits. WSBA is more realistic in the sense that it considers all weekly visits of randomly generated requests and the current patient at the same time in scenarios. Although WSBA is intuitively more accurate for this problem, empirical results of WSBA and DSBA are not statistically significantly different while the computational cost of WSBA is relatively high compared to DSBA. Therefore, Demirbilek et al. (2018) employ DSBA for further experiments. Test results show that DSBA increases average daily visits by up to 10% and reduces travel times per patient by up to 13% compared to DH and CH.

Solution methods for multiple nurses case
One way to adapt single nurse approaches to the case of multiple nurses is to split a service region into several districts. We explain the drawbacks of this intuitive approach, and present alternative heuristics that operate without territorial assignments.

Nurse districting problem
Districting problems, also called territory design, territory alignment, zone, or sector design, are concerned with defining areas in a geographical region in order to distribute scarce resources into those areas effectively. Effectiveness depends on some criteria such as balance, contiguity, and compactness. Balance can be described in terms of workloads of workers and the number of customers. Contiguity and compactness are related to the geographical shape and boundaries of territories, and has effect on travel times (Kalcsics 2015). The districting problem has a broad range of application areas such as political, school, waste collection, police patrolling area districtings. The districting problem for HHC is to define subregions that minimise travel times, balance workloads of nurses, and maximise acceptance of patients. It is not straightforward since demand and population fluctuations in subregions cause workload inequities between nurses. In addition, if more qualified nurses are scarce resources, the decision of which subregions they are assigned to is a challenging problem. Since the subregions are fixed, they cannot react to the dynamically changing requests.
We investigate how average daily visits and travel times are affected if we use DH and SBA for multiple nurses servicing their assigned territories (as in the nurse districting problem) compared to using SBAM where nurses work without any territorial restrictions. In this study, we use idealistic settings which include a square service area, equally-sized subregions, and equal expected demand for each subregion. In real-life, districting problems are complicated and there are many studies related to determine the optimal area size for each subregion by considering different constraints (Kalcsics 2015). We refer readers to Blais et al. (2003) and Nasir and Dang (2017) for more information on the districting problem and resource dimensioning in HHC.

Distance and capacity heuristics for multiple nurses
The distance heuristic for multiple nurses (DHM) is a greedy method which assigns a new request between the pair of already scheduled patients with the smallest insertion cost/additional travel time. If the Euclidean distance between request j and its predecessor i and successor i + 1 are represented as m(g i , g j ) and m(g j , g i+1 ) and the Euclidean distance between its predecessor and successor is m(g i , g i+1 ) , the insertion cost C is calculated as: Whenever a new patient arrives to the system, the algorithm calculates the cost of insertion of that patient between all requests assigned already in each day of the week and nurse. After that, the method finds intervals with the cheapest insertion costs according to visit frequency of the patient for each nurse and sums them up to be able to select the nurse with the smallest insertion cost. Finally, all visits are scheduled to those cheapest days and time slots of the nurse during the service horizon of the patient. If there are several nurses who have the same insertion costs, as a tie-breaker, we assign the visit to the nurse with the fewest scheduled patient visits to balance the workload of nurses.
Both Bennett and Erera (2011) and Demirbilek et al. (2018) show that the distance heuristic outperforms the capacity heuristic under high demand, large uniform and large uniform-clustered areas for the single nurse case. In this study, we test our algorithm only in the large area and high demand case since we have many nurses. Therefore, we only compare our algorithm with the distance heuristic.

Scenario based approach for multiple nurses (SBAM)
We explain SBAM in two sections. First for the case where all nurses have the same skills and are interchangeable, then for the case where nurses have different skill levels. We propose some modifications to SBAM to test a new assignment method and propose a simple strategy in order to choose the right scheduling method under different demand settings.

3
Home healthcare routing and scheduling of multiple nurses…

SBAM without qualification constraint
Let us start by elaborating on the main differences between SBAM and SBA. First, all nurses are taken into consideration at the same time when generating scenarios. Second, DSBA has an important drawback when multiple nurses are taken into consideration. It can happen that the new patient is allocated a visit by one nurse on one day, and another nurse on another day, which violates the consistency constraint. According to our experiments, the possibility to encounter the situation above highly increases under low demand. Therefore, we develop SBAM based on WSBA by considering all nurses' schedules, the randomly generated requests and the new patient weekly visits at the same time during the scenario phase.
The algorithm starts as soon as we receive a patient request. The first question is whether or not we accept that patient. If we accept, we need to decide what time/ days and which nurse to assign. SBAM starts generating a given number of scenarios. Each scenario includes random requests sampled from known distributions of next week's demand as well as the current patient request under consideration. To find how many random requests we need to generate, we calculate the expected weekly demand. For example, if the inter-arrival time is 255 min, we expect 2 patients for each day and 10 patients for a week. It means that 10 requests are generated for each scenario and the current request is added to them.
Next the algorithm tries to assign each request in the scenario into the existing tours of the nurses. To find which nurse's tour is the most suitable for weekly visits of a request, the algorithm looks at the smallest cost over all possible insertions in each nurse route and calculates the total insertion cost for all weekly visits and the average insertion cost per visit. The request that has the cheapest average insertion cost is assigned to the most suitable nurse and removed from the scenario pool. These iterations last until all requests in the scenario are assigned or no free capacity remains for any nurse.
After we repeat this procedure for each scenario, we record how many times the patient request has been assigned. If the request has not been assigned in any scenario, we reject it. Otherwise, we find the nurse to whom the request was most frequently assigned and assignment days/times of that nurse. Algorithm 1 shows the pseudo code for SBAM. "nCombinations" represents the days on which visits of a patient can be scheduled and "nReqInTour" represents how many times the request has been assigned over all scenarios.
Let's give a concrete example to make the method more understandable. We have three nurses, A, B, and C, and a request R (ReqR) that needs three visits per week, Monday, Wednesday, and Friday. We generate scenarios to decide whether or not we accept request R. The algorithm generates 5 random requests, R1 to R5, from the distribution of next week's demand. Table 5 shows the assignment costs of random request 2 (R2) calculated with the cheapest insertion heuristic for each visit day and nurse. After finding the total assignment cost by summing up daily costs for each nurse, we select the lowest one. In this example, Nurse A provides the lowest insertion cost.
In Table 6, the "Visits" column shows how many weekly visits patients need. "Nurse" and "Total Cost" show the selected nurse at the previous step and the total insertion cost of all visits. "Average Cost" is calculated by dividing the total cost by the number of weekly visits. We need average cost to compare visit cost of different requests to select the cheapest one. In the example, RandomR2 has the cheapest insertion cost. Thus the algorithm chooses it to assign it to the weekly schedule with all its visits. It is removed from the scenario before iteration 2. Same calculation procedure is repeated in iteration 2 and results in Table 6 are observed. As can be seen, our actual request has the cheapest insertion cost this time and it is assigned to the weekly schedule with its three visits and removed from the scenario before iteration 3. Iterations last until all requests are scheduled or no more request can be inserted. After repeating this procedure for a predefined number of scenarios, we determine how many times the request has been assigned to each nurse and which time slots it has been scheduled in for each visit day. Finally, the algorithm selects the nurse/time slot combination.

SBAM with qualification constraint
In Sect. 5.3.1, we assumed that all nurses are homogeneous in terms of their skill levels. The assumption may be justifiable in real-life since a company can be specialized for only one type of nursing service. However, companies often provide a range of services to patients. If all workers in a company are just assigned to visits that exactly match their skill level, we can consider this problem as homogeneous nurses that we have dealt with before. All we need to do is to find demand for all different visits, distinguish nurses in terms of their qualifications, and construct schedules and routes for them separately. However, caregivers qualified for a particular level can also perform visits of lower qualification levels. For example, if a company employs two different nurses, senior and junior, senior nurses can be asked to perform some visits that junior nurses can perform due to lack of junior nurses at that time or lack of demand for visits that only senior nurses are qualified for. On the other hand, junior nurses are not allowed to perform some visits that require higher qualifications. Of course, hourly cost of more skilled workers' qualifications is higher than that of less skilled workers. Therefore, decision makers have to take this into account when maximising the number of visits.
In this case, we examine different assignment strategies, "Mixed", "Homo", and "Prior". The "Mixed" strategy represents mixed assignments where a patient can be assigned to a nurse if the nurse has a sufficient skill level. "Homo" is a homogeneous assignment where patients can only be assigned to nurses who are exactly matched in terms of skill levels. At last, "Prior", gives priority to patients who demand higher skilled service when assigning them in the scenario generation phase and we can assign patients who need lower skill to higher skilled nurses if and only if there is no further patient who needs higher skilled nurses in the scenario. Finally, we compare results to DHM by considering only mixed assignments.

Empirical evaluation
In this section, we present numerical results on the performance of the proposed heuristic. We begin by describing the general simulation settings used in most of the remaining sections.

Simulation set-up for homogeneous nurses case
Experiments are carried out according to different demand patterns and service area types for SBAM and the distance heuristic. We run 30 simulations for each experiment. Each simulation horizon is 360 working days. A day is composed of 34 time slots ( k = 34 ) each lasting 15 min. The nurse works between 08.00 and 16.30 each day during the planning horizon. A 20-day warm-up period is set at the beginning of the simulation. Inter-arrival times between requests are exponentially distributed with mean 340 (low demand), 255 (medium demand), and 150 (high demand) min. Each patient has to be serviced over 4 weeks with a frequency 1, 2, or 3 visits per week with probabilities 0.05, 0.35, and 0.60, respectively. The first visit starts the following week after the request has been accepted. Visit durations are deterministic and take 30 min. Each arriving customer request is equally likely to arise from a square region X ∈ [0, 60] and Y ∈ [0, 60] . However, in real-life, the number of patient requests arriving from one region can be higher than from another region. Maybe some regions are slightly or not populated. To be able to test our algorithm under those conditions, we cluster patient requests in three rectangular subregions with given coordinates X ∈ [8, 20] and Y ∈ [3, 15] , X ∈ [23, 35] and Y ∈ [18,30] , and X ∈ [43, 57] and Y ∈ [33, 47] as shown in Fig. 1. Two different cluster types are proposed for the other two scenarios. Patient requests arrive only from these subregions equally likely in Cluster 1. In Cluster 2, 70% of patient requests arise from those subregions while 30% arrive from the remaining area of the whole service region. Simulation parameters are shown in Table 7. 75 scenarios are generated for each decision and we accept a patient if he/she can be scheduled at least once over 75 scenarios (Demirbilek et al. 2018). Note that always accepting the patient would be similar to DH and leads to inferior results. We have two different set-ups for visit days each patient can be assigned to according to his weekly visit frequency. For day set 1, each patient can be scheduled any combination of days in the week. Day combinations for a patient can be n f when f represents the visit frequency and n shows the number of days (Monday, Tuesday, Wednesday, Thursday, Friday). We also use another day set-up, day set 2, which does not allow to schedule sequential days when the visit frequency of a patient is two or three since it is not realistic to perform some visits such as cooking, bathing, etc. the first two or three days at the beginning of a week and to do nothing at the remaining days. Therefore, only the following visit day combinations can be assigned to a patient who needs two visits in a week, [(Monday, Friday), (Monday, Thursday), (Monday, Wednesday), (Tuesday, Friday), (Tuesday, Thursday)] and a patient who needs three visits in a week, (Monday, Wednesday, Friday). There are three nurses in the tests and their homes are located at (10, 10), (30, 30), and (40, 50) in the region. To understand differences between results of our algorithm and the distance heuristic, we conduct independent samples t test and calculate p values for each pair. Results written in bold in tables mean that they are statistically better than those from other methods.
The parameter settings such as the simulation horizon, length of time slot, weekly visit frequencies and probabilities are taken from the study of Bennett and Erera (2011) since they derived these parameters from real data. However, we employ different inter-arrival times to model patient arrivals since we have multiple nurses and want to see how SBAM works under low and high demands where acceptance rates change between 60 and 100%. The HHC routing and scheduling problem is a rich problem in terms of objectives and constraints. It is also hard to consider each constraint and performance measure in the same problem. For example, nurses need breaks during working days and patients can be allowed to choose a nurse during their service horizons. Both could be adapted in our study easily. In our problem, nurses work between 8.00 and 16.30 in each day. We can add 30-min breaks around noon and change afternoon schedules accordingly. This could be either a hard constraint (there has to be a break from 12.00-12.30), or a flexible constraint (there has to be a 30 min break sometime between 12.00 and 14.00). This is straightforward to include in the construction heuristic for the tours. The point here is that our algorithm provides flexibility and we can change or add many things into the problem with small modifications. Similarly, our patients can be allowed to choose a nurse. All we do is to find the best visit times in the schedule of selected nurses if feasible. Moreover, our approach does not depend on many of the parameter settings, and should also be applicable for other time intervals, different appointment durations, and service horizons.

Results for homogeneous nurses case
We evaluate our heuristic using two measures: average daily visits (which we seek to maximize) and travel time per visit. Table 8 reports the result of this evaluation for the case of 3 nurses and day set 1. When the demand intensity is low (interarrival time is 340 min) for each regional demand case, there are no meaningful differences between average daily visits of the two methods, SBAM and DHM, since acceptance rates for both approaches are 100%. In other words, nurses have sufficient capacity to accept all patients. However, when we look at travel time per visit, SBAM reduces it by 24%, 51%, and 34% depending on the regional demand scenario. When the demand intensity is higher (inter-arrival times 255 or 150 min), SBAM accommodates more average daily visits than DHM. Particularly if patient demand arrives only from subregions (Cluster 1), increases in average daily visits are higher. SBAM improves average daily visits by 11% in the high demand scenario. It is important to emphasise that acceptance rate is around 85% in that case. When demand increases, the gap between average daily visits of the two methods increases as well. All differences between results except those for average daily visits under the low demand are statistically significant. Table 9 reports the result of this evaluation for the case of 3 nurses and day set 2. Because we can assign patients to only some specific day combinations, average daily visits decreases compared to day set 1. However, we can observe the same overall picture as before. In low demand scenarios, SBAM reduces travel time per visit significantly for uniform and clustered demand scenarios even though average daily visits of both methods are similar. Under the high demand, SBAM improves average daily visits by 11%, 19%, and 16% for each regional demand scenario respectively. Again, all differences between results except those for average daily visits for the low demand are statistically significant.
Another important issue for HHC providers is to ensure fair workloads for their caregivers. It is not acceptable for a member of staff to have to work substantially more than others over a prolonged period of time without any compensation. In the literature, some studies are devoted to balance workloads of workers (Hertz and Lahrichi 2009). Although our objective is to maximise total daily visits, we should examine workloads of nurses explicitly. We find a visit range for each trial, which is the difference between maximum and minimum daily visits of nurses. For example, if Nurse 1, 2, and 3 visit on average 4, 5, and 6 patients in a trial, the range for this trial is 2 visits. After that, we calculate average range for 30 trials and a confidence interval for the mean. Means closer to zero are more desirable since it points out that average daily visits of nurses are equal. "Range" in Tables 8 and 9 show average ranges and confidence intervals with 95% confidence level. The only case when the difference between daily visits of nurses clearly exceeds one visit is day set 1/ low demand/uniform scenario. In other scenarios, particularly high demands, differences between daily visits of nurses are very small. On the other hand, DHM causes unbalanced workloads in many scenarios. Differences between daily visits of nurses exceed one visit and sometimes two visits. Because SBAM considers remaining workloads of nurses when assigning randomly generated and actual requests, we can observe more balanced workloads overall compared to DHM that assigns requests only by considering insertion costs. We tested the algorithms for 3 nurses so far. For higher number of nurses, we increase demand proportionally in order to keep acceptance rates under 100%. Otherwise, we cannot understand whether or not our algorithm is superior since there are sufficient nurses to accept almost all requests for both methods. However, higher demand causes significant computational times. Thus, we also increase visit times and service horizon to 50 min and 8 weeks respectively. So caregivers spend 50 min for each visit and patients are served during 8 weeks instead of 4 weeks. Locations of nurses and patients are uniformly distributed in the service region and scenarios are tested for day set 2. Other simulation parameters are same as in Table 7. We test SBAM for 12 and 24 nurses. In this trial, all nurses are homogeneous in terms of their qualifications. Location of nurses are assigned uniformly across the service area. Tables 10 and 11 show total daily visits and travel times per visit for 12 and 24 nurses, respectively. For 12 nurses, there is no statistical difference between total daily visits of SBAM and DHM since the acceptance rate is 100% for both when inter-arrival time is 255 min. For other cases, SBAM is able to schedule significantly more visits than DHM, and the improvement is even larger for total travel times. All differences are statistically significant.
We have similar results for 24 nurses in Table 11. Although the demand is more than 5 requests every day (inter-arrival time is 100 min), acceptance rates are 100%. When the acceptance rate is around 98%, SBAM provides 4 more visits than DHM every day. It is obvious based on previous experiments that the gap between total visits of SBAM and DHM is getting larger when the demand increases. In addition, travel times per visit shorten more than 50% with SBAM.
A key requirement for our algorithm is that it should be possible to respond to patients quickly. To show how computational times increase proportionally with increasing demand and the number of nurses, we provide Fig. 2. This graphic shows average execution times of acceptance and assignment decisions based on 100 patients, day set 2, uniform service area, 4-week episode of care, and different interarrival times and number of nurses. Our algorithm is clearly more sensitive to increasing demand than number of nurses. However, the average decision time in the highest case (Interarrival time is 25 min) for a patient is still around only 1 min, which should be acceptable in practice.

Simulation set-up and results for nonhomogeneous nurses case
We test three assignment strategies for two different types of patients (Type 1 and Type 2) according to their need of nurses, junior and senior. Moreover, we have two different demand estimations for patients. First, we assume that 60% of patients need nurses qualified at least junior level and others need senior nurses. Second, 80% of patients need nurses qualified at least junior level. There are 9 nurses, 6 junior and 3 senior. Locations of nurses and patients are uniformly distributed in the service region and strategies are tested for day set 2. Nurses spend 50 min for each visit and patients are served over 8 weeks. Interarrival times are exponentially distributed with 150 min. Other simulation parameters are identical to Table 7.
In Tables 12 and 13, we show some results where "TotalAcc" refers to the proportion of patients accepted. "Type1Acc and Type2Acc" represent acceptance rates of Type 1 and Type 2 patients, respectively. According to Table 12, homogeneous assignment, "Homo", provides the highest acceptance rate for Type 2 patients while keeping total daily visits and total acceptance at the same level (even slightly better  but not statistically significant) with other strategies. This is expected since we know that although 40% of patients need senior nurses, only 33% of nurses have this skill level. Thus, dedicating senior nurses to Type 2 patients causes the highest acceptance rate for those patients. The second best, "Prior", shows that giving absolute priority to Type 2 patients and assign Type 1 patients to senior nurses only after there is no Type 2 patient in the scenario seems a good strategy since we do not know exact demands, but we can guess more or less which type of service has a higher demand. Therefore, "Prior" looks applicable under different demands and when acceptance of some kind of patients is more valuable. Table 13 shows total daily visits, travel times, and acceptance rates for 3 scenarios and DHM when 80% of patients need junior nurses and 20% need senior nurses. In this setting, total daily visits in "Homo" is lower than other strategies since dedicating all senior nurses to only Type 2 patients causes ineffective utilization. Although the difference between total daily visits of "Prior" and "Mixed" is not statistically significant, the acceptance rate of "Prior" is significantly higher than "Mixed". However, when comparing the acceptance rates of Type 2 patients in "Prior" with "Homo", the gap is massive since the possibility of accepting Type 1 patients during scenario generation phase is high even if the algorithm gives priority to Type 2 patients. Total travel time in "Prior" is higher than other strategies except "Homo". This shows that we have to accept longer travel times to accept less suitable requests for tours.
Although "Prior", where the algorithm gives priority to patients who demand higher skilled nurses in scenario generation phase provides robust results in two different demand structures, mixed and homogeneous assignment can be reasonable according to targets of companies. For example, according to acceptance rates in Table 13, if the company charges Type un1 patient service hour 100$, the most profitable strategy is mixed assignment if hourly price for Type 2 patient is up to 108$. If hourly price for Type 2 patient is considered between 108$ and 176$, "Prior" gives the highest profit. After 176$, the best strategy is to use homogeneous assignment. This is valid if the company only considers profit maximisation. If the aim is to visit as many patient as possible, the mixed assignment can be chosen.

Simulation set-up and results for nurse districting problem
We tested two different nurse sets with 2 and 4 nurses. In the first scenario, nurses are located in the centre of their service regions. The whole service region is divided  and (45, 45). In the second scenario, nurses are not located in the centre of their service regions. For two nurses, they are located at (5, 10) and (60, 25), and four nurses, their locations are (5, 10), (55, 30), (0, 50), and (60, 25). Inter-arrival times between patient arrivals are 255 and 510 min for two nurse and 100 and 200 min for four nurse cases. The expected demand is divided equally among subregions. We compare this case with the case where nurses can service the whole service area without any restriction. Of course, the overall area and locations of nurses are identical. Moreover, we also test DH performance for the case that nurses only service in their subregions and the case that they can visit patients in the whole service area (DHM). All results are statistically significantly different according to the t test .  Tables 14 and 15 show average daily visits and travel times per visit for SBA, DH, DHM and SBAM. SBAM clearly increases average daily visits compared to SBA. As expected, districting the service area reduces travel times per visit. For four nurses, the fact that nurses travel in the whole area increases travel times per visit by up to 23% while average daily visits rise by around 5%. It is considerable to sacrifice the small amount of visits in order to reduce long travel times. Note that nurses are located at the centre of their subregions and the demand in each subregion is equal in Scenario 1. On the other hand, DH increases average daily visits when nurses service in the whole area rather than only their own subregions as seen in Table 14. However, under high demand and 4 nurses, districting areas gives better results. The most important reason is that DHM fails to balance workload of nurses. Therefore, once workload of a nurse is quickly filled, patients are assigned to nurses whose tours are not suitable. This causes higher travel times and less number of visits. Table 15 clearly shows quite long travel times when four nurses service in the whole region under high demand compared to their assignment to small regions.
Although it is desirable that nurses are located close to the centre of their service regions, it is not always possible. Scenario 2 represents the situation where nurses are not located close to the centre of their service areas. Tables 16 and 17 represent average daily visits and travel times per visit for SBA and SBAM in Scenario 2. Although SBAM increases average daily visits in each case as in Scenario 1, increments are higher compared to results in Table 14. Moreover, changing the location of nurses affects results in SBA more than in SBAM. For example, average daily visits of nurses by SBA decrease around 10% while visits by SBAM decrease only 2%. On the other hand, DH shows the same pattern as in Scenario 2. When there are four nurses and the demand is high, assigning nurses to districts increases average daily visits and lessens travel times per visit.
Overall, results show that considering the whole service area, demand over the whole service area and all nurses at the same time when routing and scheduling provides better results in terms of average daily visits in each scenario. Moreover, SBAM seems more robust against changes of nurses' locations and expected demands. It is notable that Scenario 1, equal size territories, equal demand and central nurse locations, is an idealistic condition for a districting problem. Therefore, average daily visits and travel times per visit are near optimal. When some conditions change as in Scenario 2, results for districting approach also deteriorate.

Conclusion and future work
Importance of home healthcare (HHC) is growing day by day since populations of developed and even developing countries are getting older quickly and the number of hospitals, retirement homes, and medical staff do not increase at the same rate.
In this study, we consider that patients have to be serviced on the same days and times by the same nurse during their entire service horizon. Furthermore, acceptance and assignment time decisions have to be made as soon as patient arrives. We propose a Scenario Based Approach for Multiple Nurses (SBAM), which is based on generating several scenarios that include current schedules of nurses, randomly generated requests mimicking future demand, and the patient under decision, constructing tours by using a simple but fast heuristic, and analysing results to decide acceptance, nurse, and appointment days and times. The basic idea behind the algorithm is to run a number of simulations (scenarios) to see how many times the patient's visit is assigned among all requests and which time slot and nurse it is allocated to most frequently. Based on this information, we decide to accept or reject the patient and time slots and nurse.
Compared to our previous approach for a single nurse (Demirbilek et al. 2018), we make important modifications to the algorithm to be able to work with multiple nurses. First, we consider all nurses' current schedules when assigning current and random requests at the same time in the scenario generation phase. Second, we calculate assignment costs by considering all weekly visits of a request and compare average costs to choose a request in a scenario. On the other hand, we test how SBAM works based on all nurses service in the whole area versus SBA by assigning each nurse to a district. SBAM provides significantly higher average daily visits and acceptance rates albeit longer travel times. After that, we test SBAM against DHM for three, twelve, and twenty-four nurses under different demand structures, inter-arrival times and day sets. Results show that SBAM significantly increases average daily visits and decreases daily travel times per visit for each scenario. In high demands, average daily visits increase around 20% and daily travel times per visit are reduced by up to 50%. Finally, we test our algorithm if nurses are not homogeneous in terms of their skill level. Three different strategies in terms of assignment structure are reviewed. The purpose is to propose different strategies to decision makers according to their needs or targets. We show how to choose a strategy in terms of hourly service prices of different patients.
In future research, we plan to extend this study to include aspects of revenue management. In that case, patient preferences in terms of visit days and times are evaluated and suitable prices are determined according to optimal schedules. Furthermore, some patient visits need more than one nurse due to their complexities and optimising routing and scheduling of nurses under this constraint seems quite interesting as well as a challenge for future research.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.