The multi-period home healthcare routing and scheduling problem with electric vehicles

The home healthcare (HHC) sector covers many services including the medical, therapeutic, and non-medical services. Due to providing services to patients in their own environment, it is considered to be convenient, cost efficient, and effective (Fikar and Hirsch 2017). The global HHC sector is one of the fastest growing and highest valued sectors. This is not surprising as the world’s population ages and their demand for medical services increases. According to the World Health Organization (WHO), there were in total 703 million individuals aged 65 years or over in the world in 2019. This number is expected to double by the year 2050. Furthermore, the total global HHC market was valued at USD 281.8 billion in 2019 and it is projected to raise at a compound annual growth rate of 7.9% from 2020 to 2027 (Grand View Research 2020). A number of key factors are driving this growth. These factors include ever growing government investments, increase in chronic diseases across the regions as well as growing awareness of the benefits of the HHC (Market Data Forecast 2020).

Recent global pandemic also showed the importance of the HHC. Thousands of individuals who have contracted COVID-19 but do not require hospitalization may need home healthcare, as well as those individuals who were hospitalized with the virus may need home healthcare during their post-discharge period. Furthermore, some older adults, who are less willing to move into nursing homes following the news regarding tens of thousands of deaths in nursing homes, may also choose the HHC instead (The Conversation 2020).

The HHC scheduling is inherently a complicated problem. Healthcare professionals, nurses in this case, have to be routed to many patients at different locations during the planning horizon. This scheduling therefore involves the temporal (time windows on when to see the patients), assignment (right nurses to right patients), and geographic types (patients at different locations) of constraints (Cisse et al. 2017). Furthermore, similar to other sectors, in the HHC sector, the efficient use of resources (mainly, healthcare professionals) is a need. A study on time consumption in the HHC sector in Norway showed that the average driving time accounted for 20% of the total shift time (Holm et al. 2014). Together with the high costs associated with the HHC services, this shows that there is still room for improvement in the HHC sector. This is what the research area known as the home healthcare routing and scheduling problem (HHCRSP) aims to do; improving efficiency and decreasing costs in the HHC sector.

The HHCRSP is an extension of the vehicle routing problem (VRP) as it involves routing of multiple nurses. VRP is a topic that has been studied more than six decades by many academics. After 2000 s, the interest in the topic has greatly increased, mainly due to its increasing economic importance. The multi-period HHCRSP is a further extension of the well-known and studied aforementioned VRP, where the planning horizon of the problem is composed of a number of periods (e.g., days) instead of a single period. In real-world circumstances, depending on their health conditions, patients may require multiple visits even during the same working day. As an example, blood pressure checks and blood glucose monitoring require multiple visits on the same working day. Therefore, HHC services are typically planned weekly and one week ahead to inform the nurses (Trautsamwieser et al. 2014) and the patients to be served, resulting in the multi-period HHCRSP.

With ever-increasing environmental concerns in recent years, more research has been published on how to integrate environmental aspects to VRP and its variants (Bektaş et al. 2019). Moreover, electric vehicles (EVs) are becoming more popular as people get more environmentally conscious. As an example, the global electric car fleet exceeded 7.2 million in 2019, it was 40% higher than the previous year’s figure. China was the world’s largest market with nearly 3.4 million total electric cars in 2019 followed by Europe and USA with 1.7 million and 1.5 million, respectively (International Energy Agency 2020). Furthermore, governments across the globe are implementing policies that limit the sale of fossil fuel-based vehicles in order to further incentivize the use of more environmental friendly vehicles. Such is the case with the British Government policy that would end fossil fuel car sales by 2035 (Transportation and Environment 2020). This causes increase in adaption and usage of EVs in the HHC sector. Accordingly, this study considers the use of EVs as means of transportation for the nurses.

This paper investigates the multi-period HHCRSP with electric vehicles (MP-HHCRSP-E). Even though the problem addressed in this paper is observed in other on premise service sectors, MP-HHCRSP-E has healthcare services specific objectives and constraints—that will be further explained in detail—such as; the possibility of multiple visit requirement for a single patient during the planning horizon and even on the same day, patient and nurse skill/competency requirements, time windows for both nurses and patients, and patient specific service times that are particular to the healthcare sector.

To the best of our knowledge, despite the fact there exists extensive amount of research into each part of the problem separately, there is no research that investigates the impact of all these aforementioned factors integrated together. Next, we provide a brief literature review of the studies on the multi-period HHCRSP and on the electric VRP (E-VRP) with a summary on the scientific contributions of this study to the literature.

1.1 A brief literature review on the multi-period HHCRSP

In the HHCRSP literature, less research has been conducted on the multi-period version compared to the single-period version. Next, we review them based on their objectives, constraints, and solution approaches compared to this study.

In multi-period HHCRSP literature, several studies considered different objectives such as travel time minimization (Begur et al. 1997; Bowers et al. 2015), traveling cost minimization (Bard et al. 2014a, b), traveling distance minimization (Maya Duque et al. 2015; Nickel et al. 2012), and minimization of the number of assigned healthcare workers (Barrera et al. 2012; Rodriguez et al. 2015). Compared to the single-period problems, even though in general the travel time was still a major focus of the multi-period papers, a greater focus was on staffing related issues. Considering all costs associated with the HHC services, the objective is to minimize the total cost that comprises the fixed cost of utilizing healthcare nurses, the energy charging costs, the costs associated with depot-to-nurse home transfer services, and the costs of a patient left unserved. We see that studies have considered constraints such as time windows (Trautsamwieser et al. 2014) and task skill requirements (Bard et al. 2014a, b; Cappanera et al. 2015).

Skill level requirements have been used in various ways. In a study by Bard et al. (2014a), therapists could only serve the patients if they had the required certification for the therapy type the patients requested. In a study by Cappanera et al. (2015), patients and nurses had skill level requirements, where nurses could work on patient tasks with equal and/or lower level of skill requirements. This is also how our problem employs the nurse and patient task skill level requirements, as nurses could only serve the patients with lower or equal skill level requirements.

Majority of the past studies have assumed that all the tasks had to be served. Nevertheless, similar to our study, a number of studies (e.g., Manerba et al. (2016)) considered the cases where the workload of nurses was more than they could handle, as a result the patients to be served had to be selected. Similarly, Cinar et al. (2021) considered a case where the number of tasks that are to be served by each nurse is very high, as a result some patients could not be served. The authors used a time-dependent prize matrix to effectively assign patients to the right days in addition to finding optimal routes for the patients for each day in the planning horizon.

A number of past studies used instances provided by a real home healthcare provider (Cinar et al. 2021). Similarly, Maya Duque et al. (2015) used the real-world data to provide a solution to a nonprofit Belgian HHC organization that provides service at different provinces. Quintanilla et al. (2020) also used the real-world data of a Spanish hospital that utilized a particular taxi company to serve HHC patients. Nevertheless, there are also studies that used both real-world data and randomly generated data together (Nickel et al. 2012; Begur et al. 1997) similar to our study.

Depending on the level (tactical or operational) of decisions, some studies used planning horizons longer than a month. As an example, considering a planning horizon of two or three months, Salman et al. (2021) researched into the possibility of mobile health service delivery to Syrian migrant agricultural workers, who might otherwise not have an easy access to such services, using mobile clinics that would provide service during the day time and return to a depot after work hours. Considering a short-to-medium term problem in this study, we consider a planning horizon of multiple consecutive days.

There have also been past studies which research into the continuity of care during MP-HHC problems (Savaşer et al. 2022). In their study, patients in villages at rural area served by the same practitioner at all visits. Although this study does not enforce such a constraint, it can easily be incorporated as described in Sect. 2.3.

Exact solution algorithms were commonly used in past articles (Cappanera et al. 2015; Rodriguez et al. 2015), in conjunction with heuristic approaches (Trautsamwieser et al. 2011; Bard et al. 2014a; Shao et al. 2012). Table 1 summarizes the literature on the multi-period HHCRSP. Readers are referred to Fikar and Hirsch (2017) for a detailed review of the solution approaches used in past research for both multi-period and single-period HHCRSP.

1.2 A brief literature review on the E-VRP

The E-VRP is a natural extension of the classical VRP that encompasses the specific constraints and optimization issues that are associated with EVs such as the limited battery capacity, charging time, and locations related constraints (Schneider et al. 2014). Unlike fossil fuel-based vehicles, EVs have more specific constraints such as the range limits due to battery capacity, variation and relative slow speed of recharge energy, and limited charging station (CS) availability (Pelletier et al. 2016). In recent years, there has been plenty of research into this topic (Koç et al. 2019; Lee 2020).

Several past studies considered EVs with homogeneous fleet compositions (Keskin et al. 2018). A number of studies also considered different technology and capacity options for charging EVs such as full charging policy using both identical technologies (Schneider et al. 2014) and different technologies (Basso at al. 2019) as well as policies that allow for partial charging with both identical technologies (Erdem and Koç 2019) and different technologies (Koç et al. 2019). As lack of sufficient infrastructure for CS is one of the problems that especially affect long distance travel of EVs, there has been research into the distribution of CS for EVs (Gacias and Meunier 2015). Time windows (Keskin et al. 2021) and total shift duration (Ceselli et al. 2021) are some of the other operational constraints that have been considered. Due to the complex nature of the problem, many authors used a variety of heuristics.

Number of past studies which utilized exact approaches as a solution method in E-VRP is relatively small (Tahami et al. 2020; Desaulniers et al. 2016). On the contrary, many past studies utilized heuristics to solve the E-VRP as we do in this study. Such heuristic approaches included adaptive large neighborhood search (ALNS) (Goeke and Schneider 2015; Keskin et al. 2018; Koç et al. 2019), adaptive variable neighborhood search (AVNS) (Bruglieri et al. 2017), and iterated local search (ILS) (Montoya et al. 2017).

In our study, we combine the aspects of both E-VRP with MP-HHCRSP. Our study investigates the impact of MP-HHCRSP with use of EVs that have their own specific constraints such as the ones related to battery capacity and recharge energy variations. Our study, to the best of our knowledge, is the only study that combines these two research topics in VRP in one.

1.3 Scientific contributions and structure of the paper

This paper makes three main scientific contributions. Firstly, we introduce the multi-period electric home healthcare routing problem with time windows, which includes different possible starting and ending locations, several charger technologies, patient competency requirements, and nurse competency levels. Secondly, we propose an ALNS algorithm to effectively solve the problem. It starts with forming the initial solution using a construction heuristic algorithm, and then improves the solutions through sequential neighborhood change intensification mechanism as well as problem specific destroy and repair operators. Thirdly, we deeply analyze the problem parameters to provide managerial insights for HHC providers using EVs.

The remainder of this paper is structured as follows. Section 2 provides a formal description as well as the mathematical formulation of the problem. Section 3 contains a detailed description of the ALNS algorithm. Section 4 presents the computational experiments, followed by our conclusions in Sect. 5.

This section presents the problem definition, notation, and mathematical formulation.

2.1 Problem definition and notation

The problem aims to construct the weekly routes of healthcare nurses that provide service to patients in a geographic region. The patients that need to be served during the planning horizon are indicated by the set \({\mathcal {I}}\). The planning horizon that typically corresponds to a working week comprises multiple consecutive days and is indicated by the set \({\mathcal {D}}\). Each nurse is assigned to an EV from a homogeneous EV fleet and uses it throughout the planning horizon. The term EV will be used instead of the term nurse in the reminder of the paper. The set of EVs is denoted by \({\mathcal {K}}\). Although the EV fleet is homogeneous, as the nurses differ by their competency levels, the set \({\mathcal {K}}\) is a heterogeneous set. Each EV \(k \in {\mathcal {K}}\) has to be collected and deposited back to the depot at the beginning and end of the planning horizon by each nurse. At the beginning of the planning horizon, all EVs are fully charged and ready to be used by the nurses. The battery capacity and the amount of energy consumed per unit distance of an EV are denoted by Y and c, respectively.

During the working week, each EV \(k \in {\mathcal {K}}\) could be charged at a CS by using one of the three different technologies, normal, fast, and super-fast, as in Keskin et al. (2018). The set of charger technologies is indicated by \({\mathcal {M}} = \{1, 2, 3\}\) where \(m = 1\), \(m = 2\), and \(m = 3\) specify normal, fast, and super-fast chargers, respectively. The set of CSs is denoted by \({\mathcal {S}}\). In order to track the multiple visits of a CS by the same EV, we introduce multiple copies of each CS. The set of dummy nodes indicating the multiple visits of a CS in \({\mathcal {S}}\) is denoted by \({\mathcal {S}}^\prime\). The charging rate of charger type \(m \in {\mathcal {M}}\) is defined as \(r_{m}\). The unit energy cost of using charger technology \(m \in {\mathcal {M}}\) is defined by the parameter \(\alpha ^C_m\). A nurse may charge her/his EV during the working day between consecutive patient visits using any of the charging technologies. Furthermore, after completing daily visits, a nurse may charge her/his EV at the depot using normal charging technology. When an EV is deposited to the depot at the end of a working day, the nurse should take it from the depot at the beginning of the upcoming work day. A fixed transportation fare support, denoted by \(\alpha ^{S}_k\), is paid to nurse \(k \in {\mathcal {K}}\) to compensate her/his travel from/to depot to/from home.

Since each patient \(i \in {\mathcal {I}}\) may request service more than once during the planning horizon, a copy of each patient \(i \in {\mathcal {I}}\) is introduced for each of her/his service request and the corresponding task set is indicated by \({\mathcal {I}}^\prime\). The binary parameter \(e_{ij}\) indicates whether tasks \(i, j\in {\mathcal {I}}^\prime\) correspond to the same patient. Whether task \(i \in {\mathcal {I}}^\prime\) should be served on day \(d\in {\mathcal {D}}\) is indicated by the binary parameter \(d_{id}\). The parameter \(\lambda _{k}\) and \(h_{k}\) indicate the starting and ending home location of the nurse using EV \(k \in {\mathcal {K}}\), respectively. The starting and ending depots of EV \(k \in {\mathcal {K}}\) are represented by \({0_k}\) and \({n_k}\), respectively. Therefore, the set \({\mathcal {N}}_{k} = {\mathcal {I}}^\prime \cup {\mathcal {S}}^\prime \cup \{0_k\, n_k\, \lambda _{k}, h_{k}\}\) comprises all nodes that can be visited by EV \(k \in {\mathcal {K}}\) during the planning horizon.

The problem is defined on a complete directed graph \({\mathcal {G}} = ({\mathcal {N}}, {\mathcal {A}})\), where \({\mathcal {N}}= \bigcup _{k \in {\mathcal {K}}} {\mathcal {N}}_{k}\) is the set of nodes and \({\mathcal {A}} =\{(i, j): i, j\in {\mathcal {N}}\}\) is the set of arcs. Here, \(s^d_{ij}\) and \(s_{ij}\) are indicated as the distance and travel time of an arc \((i, j) \in {\mathcal {A}}\), respectively. On each day of the planning horizon, each EV \(k \in {\mathcal {K}}\) is available to work at a specific time interval indicated by \([a^\prime _{k},b^\prime _{k}]\) and overtime is not allowed. Each patient task \(i \in {\mathcal {I}}^\prime\) must be visited within a specific time interval indicated by \([a_{i},b_{i}]\). If an EV that is assigned to task \(i\in {\mathcal {I}}^\prime\) arrives the location of the corresponding patient before the start time of the time window of the task, i.e., \(a_{i}\), then it should wait until \(a_{i}\). The service duration of task \(i \in {\mathcal {I}}^\prime\) is denoted by \(\tau _i\).

Each task has a competency level requirement defined through the parameter \(q_{i}\). The competency level of EV \(k \in {\mathcal {K}}\) is defined by the parameter \(q^\prime _k\). A task \(i \in {\mathcal {I}}^\prime\) can be served by an EV \(k \in {\mathcal {K}}\) having a competency level of at least \(q_{i}\), i.e., \(q^\prime _k \ge q_i\). Due to shortage in the workforce, some patients might not be served within the planning horizon. Unable to serve a patient of task \(i \in {\mathcal {I}}\) results in a monetary loss of \(\alpha ^{R}_i\). Utilizing nurse \(k \in {\mathcal {K}}\) incurs a fixed cost of \(\alpha ^{N}_k\).

In order to make the problem definition clear, Fig. 1 visualizes a feasible solution to a three-day long planning horizon HHCRSP having three EVs. In the figure, the home location that corresponds to the starting node \(\lambda _{k}\) and ending node \(h_{k}\) for each EV \(k \in {\mathcal {K}}\) is shown through home symbol. There are four patients indicated by P1,..., P4 and three tasks for each patient indicated by T1,..., T3. Hence, T1P1 corresponds to task 1 of patient 1. The dashed arcs in the figure correspond to the transfers of nurses from/to their homes to/from the hospital/depot the nodes, i.e., \(0_k\) or \(n_k\), by taxi and the related transfer fee is given next to each dashed arc. The competency level, \(q^\prime _k\), of each nurse \(k \in {\mathcal {K}}\) is given next to the associated EV symbol. The time interval \([a_{i}, b_{i}]\) and the competency requirement level \(q_{i}\) of each task \(i \in {\mathcal {I}}^\prime\) are provided next to the node corresponding to the task. The percentages written at the beginning and ending points of the arcs show the charging status of EVs at the corresponding points.

A visualization of a three-day long planning horizon problem

Full size image

For the solution provided in Fig. 1, on the first day, nurses/EVs 2 and 3 are transferred from their homes to the hospital/depot for a cost of 50 units each, which is considered as the fixed one way transfer cost. Nurse 2 then performs the first task of patient 3 (indicated by T1P3), returns to the hospital with 60% charge left, and finally is transferred back home. Nurse 3 performs the first task of patient 2 and returns back home with 50% of charge capacity left in her EV. On the second day, nurse 3 departs from her own home to perform the third task of patient 2 and afterward charges the EV and returns back home with 51% charge. Nurses 1 and 2 are transferred from their home to the hospital by transfer services and after performing their tasks and returning the EVs back to the hospital, they are again transferred back to their homes. On the third and last day, only nurse 3 works and after serving three tasks and parking the EV at the hospital, she is transferred back to her home. Due to not having the required competency level of 3, no nurse is able to serve any of the tasks of patient 4.

2.2 Mathematical formulation

The sets, parameters, and decision variables of the mathematical model are presented below, followed by the objective function and the constraints.

Sets

\({\mathcal {K}}\) :

Set of EVs

\({\mathcal {D}}\) :

Set of days

\({\mathcal {S}}\) :

Set of CSs

\({\mathcal {S}}^\prime\) :

Set of multiple copies of all CSs in \({\mathcal {S}}\)

\({\mathcal {M}}\) :

Set of charger technologies

\({\mathcal {I}}\) :

Set of patients

\({\mathcal {I}} ^\prime\) :

Set of multiple copies of all patients in \({\mathcal {I}}\)

\({\mathcal {N}}_{k}\) :

Set of all nodes including the starting/ending depot of EV k and the starting/ending home of the nurse assigned to EV k, i.e., \({\mathcal {N}}_{k}={\mathcal {I}}^\prime \cup {\mathcal {S}}^\prime \ \cup \{0_k, n_k, \lambda _{k}, h_{k}\}, \; \forall k\in {\mathcal {K}}\)

\({\mathcal {N}}\) :

Set of all nodes, i.e., \({\mathcal {N}}= \cup _{k \in {\mathcal {K}}} {\mathcal {N}}_{k}\)

\({\mathcal {A}}\) :

Set of arcs between nodes

Parameters

\(s^d_{ij}\) :

Distance of arc \((i,j) \in {\mathcal {A}}\)

\(s_{ij}\) :

Traveling time of arc \((i,j) \in {\mathcal {A}}\)

\(e_{ij}\) :

Binary parameter indicating whether tasks \(i, j\in {\mathcal {I}}^\prime\) correspond to the same patient

\({[}{{a}_{i}},{{b}_{i}}{]}\) :

Time window of task \(i \in {\mathcal {I}}^\prime\)

\({[}{{a}^{\prime }_{k}},{{b}^{\prime }_{k}}{]}\) :

Time window of nurse using EV \(k \in {\mathcal {K}}\)

\(\lambda _{k}\) :

Starting home of EV \(k \in {\mathcal {K}}\)

\(h_{k}\) :

Ending home of EV \(k \in {\mathcal {K}}\)

\(0_k\) :

Starting depot of EV \(k \in {\mathcal {K}}\)

\(n_k\) :

Ending depot of EV \(k \in {\mathcal {K}}\)

\(\tau _i\) :

Service duration of task \(i \in {\mathcal {I}}^\prime\)

\(d_{id}\) :

Binary parameter indicating whether task \(i \in {\mathcal {I}}^\prime\) should be served on day \(d\in {\mathcal {D}}\)

\(q_{i}\) :

Competency level of task \(i \in {\mathcal {I}}^\prime\)

\(q^\prime _k\) :

Competency level of EV \(k \in {\mathcal {K}}\)

Y :

Battery capacity of an EV

c :

Amount of energy consumed per unit distance taken by any EV

\(r_{m}\) :

Recharging rate using charger type \(m \in {\mathcal {M}}\)

\(\alpha ^{N}_k\) :

Fixed cost of utilizing a healthcare nurse using EV \(k \in {\mathcal {K}}\)

\(\alpha ^C_m\) :

Unit energy cost of using charger technology \(m \in {\mathcal {M}}\)

\(\alpha ^{S}_k\) :

Fixed cost of utilizing a transfer service for \(k \in {\mathcal {K}}\)

\(\alpha ^{R}_i\) :

Cost of task left unserved \(i \in {\mathcal {I}}\)

M :

a sufficiently large number

Decision variables

\(x_{ijkd}\):

1 if EV \(k \in {\mathcal {K}}\) travels on arc \((i,j) \in {\mathcal {A}}\) on day \(d\in {\mathcal {D}}\), 0 otherwise

\(\delta _{ik}\):

1 if EV \(k \in {\mathcal {K}}\) serves task \(i \in {\mathcal {I}}^\prime\), 0 otherwise

\(\beta _{ikd}\):

1 if EV \(k \in {\mathcal {K}}\) is recharged with normal charger at CS \(i\in {\mathcal {S}}^\prime\) on day \(d\in {\mathcal {D}}\), 0 otherwise

\(\beta ^\prime _{ikd}\):

1 if EV \(k \in {\mathcal {K}}\) is recharged with fast charger at CS \(i\in {\mathcal {S}}^\prime\) on day \(d\in {\mathcal {D}}\), 0 otherwise

\(\psi _{i}\):

1 if task \(i \in {\mathcal {I}}\) is unserved, 0 otherwise

\(z_{k}\):

1 if EV \(k \in {\mathcal {K}}\) is utilized, 0 otherwise

\(\mu _{kd}\):

1 if EV \(k \in {\mathcal {K}}\) needs to be transferred to her/his home on day \(d\in {\mathcal {D}}\), 0 otherwise

\(\rho _{kd}\):

1 if EV \(k \in {\mathcal {K}}\) needs to work on day \(d\in {\mathcal {D}}\), 0 otherwise

\(y_{ikd}\):

State of charge (SoC) of EV \(k \in {\mathcal {K}}\) at the arrival time of node \(i \in {\mathcal {N}}\) on day \(d\in {\mathcal {D}}\)

\(g_{ikd}\):

SoC of EV \(k \in {\mathcal {K}}\) at the departure time of CS \(i \in {\mathcal {S}}^\prime\) on day \(d\in {\mathcal {D}}\)

\(t_{ikd}\):

Starting time of EV \(k \in {\mathcal {K}}\) at node \(i \in {\mathcal {N}}\) on day \(d\in {\mathcal {D}}\)

\(w_{ikd}\):

Charging duration of EV \(k \in {\mathcal {K}}\) at CS \(i\in {\mathcal {S}}^\prime\) on day \(d\in {\mathcal {D}}\)

\(\sigma _{ikmd}\):

Amount of recharged energy at CS \(i\in {\mathcal {S}}^\prime\) of EV \(k \in {\mathcal {K}}\) using charger \(m \in {\mathcal {M}}\) on day \(d\in {\mathcal {D}}\)

subject to

The objective function (1) aims to minimize the total cost, which comprises four terms. The first term corresponds to the total fixed cost of healthcare nurses providing service to patients. The second term refers to the total cost of recharged energy. The third term corresponds to the total monetary loss corresponding to unserved patients. The last term corresponds to the total cost of utilizing transfer services for the nurses.

Constraints (2) guarantee that each task of an accepted patient is served by an EV and any task of an unserved patient should not be served by any EV. Constraints (3) ensure that a served task should be in the route of an EV in the day that the task should be served. Constraints (4) and (5) determine whether a nurse is utilized. In constraints (4), M corresponds to a sufficiently large number. Constraints (6) guarantee the conservation of flow for the EVs in the model. Constraints (7) and (8) ensure that an EV tour must be completed within its daily working time and that each accepted patient must be visited in his/her time window, respectively, where M corresponds to a sufficiently large number. Constraints (9) and (10) calculate the node visit times for each day of the planning horizon. Constraints (11)–(16) determine the earliest and the latest times that EVs can be in the depot or home considering the transfer time to their home if transfer services are used. Constraints (17) and (18) compute the SoC of an EV at the arrival time of a node. Constraints (19) compute and ensure the SoC limits of EVs at the time of departure from a CS or a task. Constraints (20) and (21) calculate the amount of energy charged and the charging duration for an EV that visits a CS on a day in the planning horizon, respectively. Constraints (22)–(24) select the charging technology as normal, fast, or super-fast, respectively, for an EV charged at a CS node. Constraints (25) and (26) are used to compute and limit the departure SoC from depot or home based on the previous day’s transfer information of nurses. Constraints (27)–(29) calculate the recharged energy in the depot and guarantee that EVs are fully charged on the next day if they are left in the depot at the end of a day. Constraints (30) are used to transfer the remaining charge to the next day for an EV that has returned home. Constraints (31) guarantee that each accepted task is assigned to a nurse having an appropriate competency level for that task. Constraints (32)–(35) guarantee that only when the last node of the tour is the depot for an EV, the transfer fee is charged and likewise only when the EV starts its tour from the depot, the transfer fee is charged, respectively. Constraints (36)–(38) are used to determine whether an EV works on a given day. Finally, constraints (39)–(47) define the domains of the decision variables.

2.3 Considering the continuity of care constraints

The proposed model can be easily adjusted to handle the cases where the patients prefer the same nurse at all her/his visits, through including the following constraints:

Being in the class of VRPs, the addressed problem is an NP-hard problem. Therefore, finding a feasible solution with readily available commercial solvers within acceptable time, especially for real-size problem instances, is not achievable. To obtain good solutions within reasonable time for real-size instances, we propose a heuristic algorithm that finds an initial solution through a construction heuristic (CH) (Sect. 3.1) and improves the initial solution using an ALNS heuristic (Sect. 3.2) supported by a sequential neighborhood change procedure (SNC) (Sect. 3.3).

3.1 Construction heuristic

The proposed CH has two phases, where the tasks are assigned to EVs in the first phase (see Algorithm 1) and the schedule of each EV is determined in the second phase (see Algorithm 2).

The first phase of the CH starts with calculating the scores for each task \(i \in {\mathcal {I}}^\prime\) and each EV \(k \in {\mathcal {K}}\) (lines 1–2 in Algorithm 1). The score \(P_i^{(j)}\) of task \(i \in {\mathcal {I}}^\prime\) is calculated based on Eqs. (49)–(54). The tasks that have (50) earlier start time (\(a_{i}\)) and (51) earlier end time (\(b_{i}\)), (52) longer service time (\(\tau _i\)), (53) higher level of competency requirement (\(q_{i}\)), and (54) smaller number of tasks for the same patient (\(e_{ij}\)) have higher scores. Where each of these properties is associated with a term (\(\varphi\)) that has a weight (\(\omega\)) in the score calculation equation provided in (49). We note that in the proposed CH, the task scores (\(P_i^{(j)}\)’s) are used for not prioritizing but selecting the next task to be scheduled in the second phase of the CH for the tasks that are assigned to the same nurse in the first phase and have the same time window start time.

The score \(P_{k}^{(n)}\) of EV \(k \in {\mathcal {K}}\) is calculated using Eqs. (55–57). EVs that have (56) longer time windows \([a^\prime _{k^\prime },b^\prime _{k^\prime }]\) and (57) higher competency levels (\(q^\prime _{k^\prime }\)) get higher scores. Where each of these is associated with a term (\(\varphi\)) that has a weight (\(\omega\)) in the score calculation equation provided in (57).

Subsequently, during the first phase, first, for each task-EV pair, the competency and time window compatibility is determined (lines 3–11 in Algorithm 1) and next, starting from the first day of the planning horizon, task-EV assignment is determined based on the competency and time window compatibility (lines 12–19 in Algorithm 1).

Second phase of the CH (see Algorithm 2) determines the schedule of EVs given the set of tasks assigned to each EV for each day. The EV charge requirements (lines 19–27), starting locations for each day (lines 6–17), ending location for each day (lines 46–52), and time windows of nurses and tasks (lines 29–45) are considered when forming the schedule and maintaining its feasibility. During this phase if any task for any patient remains uncovered due to time window incompatibility with the initially assigned nurse (line 31), the task is tried to be assigned to another competency compatible EV whose schedule is not determined yet (lines 32–38). If no such EV exists, the patient of that task is unserved (line 35).

3.2 Adaptive large neighborhood search (ALNS)

ALNS is an effective heuristic that is used to solve the VRP variants. It was first developed by Ropke et al. (2006) and has been successfully used in multitude of routing variations in many different studies (Goeke and Schneider 2015; Keskin et al. 2021; Koç et al. 2019).

Starting with the solution of the CH, namely \(s^{initial}\); at each iteration, ALNS employs a destroy and a repair operator selected from the destroy operators set D and the repair operators set R, respectively to generate a new solution.

Five destroy operators in the set D are described below.

1. (i)

   Random route removal: This operator randomly selects a route of a randomly selected EV and removes all of its tasks.

2. (ii)

   Worst waiting time removal: The task with the longest waiting time among the tasks of a randomly selected EV is removed.

3. (iii)

   Worst distance removal: The task with the furthest distance from the previous and the next nodes in a route of a randomly selected EV is removed.

4. (iv)

   At least \(n \%\) random route destroy: At least \(n \%\) of the total task visits of a randomly selected EV are removed.

5. (v)

   Random CS removal: A randomly selected CS visit of a randomly selected EV is removed.

The selected repair heuristic destroys the current solution until a specific percentage of task or CS visits, referred to as the degree of destruction and denoted by \(\mu\), is removed. The removed tasks are inserted into the unserved task list \(\Upsilon\).

Four repair operators in the set R are described below.

1. (i)

   Random task insertion: A randomly selected task from the unserved task list \(\Upsilon\) is inserted into a randomly selected position in the route of a compatible EV, if possible. This process is repeated until all unserved tasks in \(\Upsilon\) are traversed.

2. (ii)

   Waiting time-based insertion: A randomly selected task from the unserved task list \(\Upsilon\) is inserted into a position in the route of an EV route so that the smallest possible waiting time for the EV occurs.

3. (iii)

   Distance-based insertion: A randomly selected task from the unserved task list \(\Upsilon\) is inserted into a position in the route of an EV route so that the total distance between the previous and next nodes is the shortest possible.

4. (iv)

   CS insertion: Calculating the SoC for each task visit of each EV, a CS visit is inserted just before each task visit having negative SoC for a CS having the shortest distance to the node of the corresponding task visit, using the slowest charging technology possible.

The roulette-wheel selection approach is used to select a destroy and a repair operator among the corresponding sets based on their weights, namely \(W_h\)’s. After a destroy and a repair operator are used to obtain \(s^{candidate}\) from \(s^{current}\), the scores of the corresponding heuristics, namely \(S_h\)’s, are incremented by \(\Delta S_h\), which is determined by using the following Eq. (58).

where \(c(s^{best})\), \(c(s^{candidate})\), \(c(s^{current})\) correspond to the best, candidate, and current objective function values and \(\Lambda _1\), \(\Lambda _2\) and \(\Lambda _3\) are set as \(\Lambda _1>\Lambda _2>\Lambda _3\). The scores of the operators \(S_h\) are initially zero and reset in every \({\Omega }\) iterations. The weights of destroy and repair operators, i.e., \(W_h\)’s, are initially equal and adjusted in every \({\Theta }\) iterations based on Eq. (59)

where the parameter \(\zeta\) corresponds to the weight of the historical performance of the heuristic and the parameter n(h) denotes the number of times that h is used in the last \({\Theta }\) iterations.

The acceptance criteria of the simulated annealing algorithm are employed for the acceptance of the new solution. T and \(C_r\) denote the current temperature and the cooling rate, respectively. The initial temperature is set to a value so that a candidate solution, which is \(50\%\) worse than \(s^{initial}\), could be accepted with \(50\%\) probability.

As the stopping condition, maximum number of iterations, namely \(iter^{max}\), is used and iter denotes the iteration number. In every \(\gamma\) iterations, to further improve the current solution, the SNC procedure (Sect. 3.3) is used.

Algorithm 3 describes the general framework of the proposed ALNS.

3.3 Sequential neighborhood change procedure (SNC)

In order to improve the current solution s, the proposed ALNS algorithm makes use of a SNC, which is an effective and widely used approach to improve the performance of heuristics (Hansen et al. 2019). The procedure systemically changes the neighborhood structure on the current solution s, where at each iteration a neighborhood structure is used from the neighborhood set \(N^{SNC}\). Starting from the first neighborhood in \(N^{SNC}\) as the current neighborhood, at any iteration if a neighboring solution \(s^n\) in the current neighborhood yields a better solution then the current solution s, the current solution s is updated and the next iteration continues with the first neighborhood. If not, the process continues with the next neighborhood in \(N^{SNC}\) until all neighborhoods in \(N^{SNC}\) are traversed.

In the proposed SNC procedure, the following ten neighborhood structures are used (i.e., \(|N^{SNC}|=10\)).

1. (i)

   Unserved patient insertion: A patient that has been left unserved in the current solution s is randomly inserted into the route of a compatible EV.

2. (ii)

   Horizontal relocation: The position of a randomly selected task served by an EV is changed within its current tour.

3. (iii)

   Horizontal swap: The position of two tasks served by the same EV on the same day are swapped.

4. (iv)

   Vertical relocation: A randomly selected task served by an EV is removed from its current tour and inserted into another randomly selected EV’s route on the same day.

5. (v)

   Vertical swap: Two tasks served by two different randomly selected EVs on the same day are swapped.

6. (vi)

   Home-depot swap: A randomly selected EV’s home-return for a day and home-depart for the next day are replaced with depot-return for the corresponding day and depot-depart for the next day.

7. (vii)

   EV route destruction: All patient visits of a randomly selected EV are removed and inserted into the route of another compatible randomly selected EV, if possible.

8. (viii)

   CS technology swap: The charging technology of a CS visit of a randomly selected EV is changed.

9. (ix)

   Depot-home swap: A randomly selected EV’s depot-return for a day (except the last day) and depot-depart for the next day are replaced with home-return for the corresponding day and home-depart for the next day.

10. (x)

    Random CS insertion: For the task visits having a state of charge (SoC) of below \(50\%\), the closest CS by distance is selected and inserted just before the task visit for a randomly selected EV.

Algorithm 4 presents the pseudo-code of the SNC procedure.

This section presents the results of the computational experiments to assess the performance of the metaheuristic compared to the MIP formulation and to obtain insights on problem parameters from realistic size problem instances. The mathematical model was solved using IBM CPLEX 20.10. A PC with 32 GB RAM and i7 2.3 GHz processor was used to carry out all the experiments. The ALNS algorithm was implemented in Python. Ten replications of the ALNS are performed for each instance. The best result of 10 consecutive runs, together with the total run time of all runs of the ALNS, is reported for each instance.

The rest of this section is structured as follows. In Sect. 4.1, we describe the data set generation and experimental setting. The comparative results for the MIP and the ALNS on small- and medium-size instances are presented in Sect. 4.2. Results on large-size instances are reported in Sect. 4.3. Finally, effects of the changes in problem parameters are analyzed in Sect. 4.4.

4.1 Data sets and experimental settings

We considered both literature and real-world data to generate our comprehensive benchmark instance set. Nurse weekly and daily working hours were generated based on the maximum possible weekly nurse working hours without overtime (Medisözlük 2020), and are set to 8 h daily and 40 h weekly. Task service times were randomly selected between 15 and 60 min. The planning horizon was taken as 2 or 3 days for small-instances and as 5 days, i.e., one working week, for medium- and large-size instances. All costs have been taken in the local currency of Turkish Lira (TRY). Nurse competency/skill levels and task competency/skill requirement levels were considered as 1, 2, and 3, where the average weekly wage for each level was considered to be 2000, 2250, and 2500, respectively (Vergi Net 2020; Türkiye Belediyeler Birliği 2020).

Among the set of existing EV models, an EV model with a high availability and relatively low price point was chosen (e.g., Opel Corsa-e) for the homogeneous EV fleet of our problem. All EVs had a range of 275 km when fully charged (Electric Vehicle Database 2020). The vehicle speed was fixed to 60 km/hour, consuming one unit of energy per unit distance traveled. The charging rates for normal, fast, and super-fast chargers were set to 3.47, 0.62, and 0.28 units of energy per minute, respectively. The costs of normal, fast, and super-fast chargers were taken as 1, 1.1, and 1.2 units, respectively (Keskin et al. 2018).

Locations of tasks, depot, hospital, and nurse homes were generated using real-world neighborhood location data of the city of Ankara, Turkey (TUIK 2020). The CS locations were determined by considering the actual CS locations in Ankara (Zorlu Energy Solutions 2020).

The task and nurse home locations for instances were determined using random generation method or population-based weighted random generation method. In random generation method, the locations were randomly selected from all available neighborhood coordinates. In population-based weighted random generation method, the locations from neighborhoods with higher population were given higher probability to be selected as home or task locations.

The unserved patient penalty was rounded to be 5000 (Istanbul Medical Center 2018) as a life time loss of not serving the patient in question. Lastly, the one way transfer cost for the nurses was taken to be 50.

In order to test the performance of the proposed heuristic and to determine the problem sizes that the MIP can solve to optimality in a reasonable time, we generated a total of 24 small-size instances with 3, 4, and 5 tasks. Half of these small-size instances have a planning horizon of 2 days where the remaining ones have a planning horizon of 3 days. A total of 42 medium-size instances having 9 to 19 tasks were generated. A total of 24 large-size instances with 38, 50, 63, and 74 tasks were generated. For each data set, in half of the instances population-based weighted random generation method is used, and the random generation method is used for the remaining.

The naming convention of the instances is as follows; ‘p’ stands for the # of patients, ‘n’ stands for the # of nurses, ‘j’ stands for the # of tasks, ‘d’ stands for the # of days, ‘cs’ stands for the # of CS, and ‘R’ or ‘W’ stands for the location generation method (random or weighted random, respectively). For example, an instance named as 3p4n9j5d3cs-R-1 corresponds to the 1st random generation method used problem instance with 3 patients, 4 nurses, 9 tasks, 5 days, and 3 CSs. Instances can be accessed at Yazır et al. (2023).

To further investigate the effect of changes in various parameters (such as the patient locations, the patient service durations, the patient service time windows, and the number of available nurses) of the problem, we generated 30 additional instances classified in the following five groups.

• Group-random-loc (G1): Patient, nurse home, and depot (i.e., hospital) locations are randomly selected.

• Group-control (G2): Patient, nurse home, and depot locations are selected using population-based weighted random generation method, where locations with higher population density have higher probability to be selected as patient and nurse home locations.

• Group-high-var (G3): Patient, nurse home, and depot locations were selected based on population-based weighted random generation method and the task service durations (\(\tau _i\)’s) were randomly set between 10 and 90 min. Therefore, compared to Group-control instances, the instances in this set have larger task service duration variance.

• Group-less-EV (G4): Patient, nurse home, and depot locations were selected based on population-based weighted random generation method and the number of nurses is lower than the other groups with an average number of tasks per nurse is 75 % larger than that of Group-control.

• Group-long-tw (G5): Patient, nurse home, and depot locations were selected based on population-based weighted random generation method and the ending times of the time windows of tasks are set to the largest possible value, i.e., \(b_i=\max _{k\in K}b'_k\), \(\forall i\in I'\). Therefore, compared to Group-control instances, the instances in this set have larger task time windows.

Furthermore, each group is divided into two subgroups based on their competence level requirements. In one subgroup, there is 1 competency level, i.e., each task can be served by any nurse, in the other there are three competency levels as in the original problem. These differences within groups are indicated by the ‘s’ in the problem. Where, for example, instance G1-7p7n16j1s5d3cs corresponds to a Group-random-loc, 7 patients, 7 nurses, 16 tasks, 1 competency level, 5 days, randomly generated problem with 3-CSs.

The values of the parameters of the CH and ALNS are provided in Table 2.

4.2 Results on small- and medium-size instances

The small- and medium-size instances were solved by the MIP with a run time limit of three hours. Table 3 presents the summary of the MIP and ALNS results in terms of average, minimum, and maximum objective function value, i.e., total cost in Turkish Lira (in TRY), and run time (in seconds) for small- and medium-size instances classified according to location generated method as random and weighted random (w-random, in short). Table 3 also reports the average, minimum, and maximum percentage deviation of objective function values of the ALNS from that of the MIP, where the percentage deviation is calculated as Dev. (%)= {((ALNS  Obj. /MIP  Obj.)*100)-100}. Tables 11 and 12 in Appendix show the detailed comparative results on small- and medium-size instances.

Table 3 shows that both the MIP and ALNS find the optimal result for small-size instances. However, the MIP finds the optimal solution in more than 100 s for four small-size instances (see Table 11) with a run time average of 448 s for all small-size instances. On the contrary, the ALNS finds the optimal solution within 30 s for any small-size instance with a run time average of 10 s.

For medium-size instances, the MIP could not find an optimal solution within the 3-hour run time limit. On the contrary, the ALNS obtained the same or better results within 437 s on average. The objective function value of the ALNS solutions is 11% on average (48% max.) better than that of the MIP solutions. As expected, as the problem size increases, the run time of the MIP increases and the solution quality for a limited run time decreases. For all eight instances (7p7n16j5d3cs-R-1, 7p7n16j5d3cs-W-1, 8p8n17j5d3cs-R-3, 8p8n17j5d3cs-W-1, 8p8n18j5d3cs-R-3, 8p 8n18j5d3cs-W-1, 9p9n19j5d3cs-R-2, and 9p9n19j5d3cs-W-3 see Appendix Table 13), ALNS found better solutions than the MIP. The ALNS solutions utilize less EVs and/or leave less patients unserved (e.g., 8p8n18j5d3cs-R-3).

In both small- and medium-size instances, the average objective function value of w-random instances is lower than that of random instances. This is due to the fact that in random instances patient and nurse home locations are further away from each other compared to w-random instances.

4.3 Results on large-size instances

Table 4 presents the ALNS results for large-size instances in terms of the total cost in Turkish Lira (TRY), run time (in seconds), and the percentage distribution of the total cost terms, i.e., recharged energy costs (in column Charge Cost (%)), cost of utilizing transfer service for the nurses (in column Transfer Cost (%)), fixed cost of utilizing healthcare nurses (in column EV Cost (%) as the term EV is used instead of nurse), the cost of unserved patients (in column Unserved Cost (%)), and the number of unserved patients. It is also divided into two parts, the upper part provides the results for the instances where there are unserved patients, the lower part corresponds to the results of the instances where there are no unserved patients.

According to Table 4, for large-size instances, the average run time of the ALNS is 1032 s and the maximum run time is slightly longer than 30 min, which indicates that the proposed heuristic is highly suitable for practical application. In only 6 of 24 large-size instances (see the upper of the table), some patients are unserved leading to a positive unserved cost having a share of 46.8% (on average) in total cost. For the remaining large-size instances where there are no unserved patients (see the lower part of the table), the fixed cost of healthcare nurses has the largest share (77.2% on average, 86.0% max., and 68.7% min.) in total cost, followed by the recharged energy cost (16.3% on average, 22.2% max., and 10.4% min.). For these instances, the charging cost and transfer cost, together, has a total share of 22.8% on average (33.1% max., and nearly 13.7% min.), demonstrating the need for optimization in routing and charging decisions.

When the total costs of random and w-random instances are compared, we see that the total cost of the w-random instances is 8038 TRY on average, which is much less than that of random instances, which is 13826 TRY. This is due to the fact that when the patient locations are totally random, the task and nurse home locations could be further away from each other compared to those generated based on population density. This leads to larger distances in between locations and higher recharging costs for random instances.

Table 5 provides detailed information on the workload of the nurses during the planning horizon for large-size instances in terms of the percentage of nurses utilized among all nurses (in column Nurses utilized (%)), the number of tasks assigned per nurse (as min., max., and avg. over all nurses), the number of days a nurse is required to work (as min., max., and avg. over all nurses), as well as the average utilized nurse time (in column Total) and its distribution among service time and travel time (in columns Travel and Service, respectively).

According to Table 5, some nurses are not utilized in some of the instances and the number of nurses utilized is considerably higher for w-random instances compared to random instances. Utilized nurses work for 4.4 days and serve 19.7 on average throughout the 5-day planning horizon. For random instances, the average number of assigned tasks per nurse is 16.6, which is considerably less than that of w-random instances, which is 22.7. On average, nearly 24% of working time of a utilized nurse is spent on traveling and nearly 30% working time of a utilized nurse is spent on serving tasks. For random instances, the average travel time of nurses (28.9%) has a considerably larger share in total working time compared to that of w-random instances (nearly 18.3%). The utilized nurse time is nearly 53.6% on average, 21% minimum and 72.4% maximum, demonstrating a considerable amount of idle time, which might be due to nurse-patient competency and time window availability mismatch.

To further analyze the nurse-patient competency mismatch, Table 6 presents the comparison of nurse and patient competency needs in large-size instances. It reports the number of unserved patients, the number of unserved tasks, the competency level of assigned nurses (as min., max., and avg.), the average competency level of all nurses, the average competency requirement of unserved tasks, and the average competency level of all tasks.

According to Table 6, for the six large-size instances for which there are unserved patients, the average competency level of unserved tasks is larger than the average competency level of all available nurses. According to Table 5, in none of these six large-size instances, all nurses are utilized.

For the instances that there are no unserved patients (except two instances, 18p4n50j5d3cs-W-3 and 33p6n74j5d3cs-R-3), the average competency level of all nurses is larger than the average competency level of tasks. For these two instances (18p4n50j5d3cs-W-3 and 33p6n74j5d3cs-R-3), the average competency level of assigned nurses is larger than the average competency level of tasks. Hence, more patient tasks could be served and unserved costs could be decreased by increasing the competency levels of nurses. This would be beneficial when the cost of providing training a nurse is relatively lower than the cost incurred when a patient is unserved.

4.4 Analysis on varying problem parameters

This section analyzes the effects of changes in the parameter values such as patient locations, patient service durations, patient service time windows, number of available nurses, and existence of multiple competency levels. The instances in Group-random-loc (G1) to Group-long-tw (G5) are used in this analysis.

Table 7 presents the total cost, the distribution of cost terms, the average number of nurses utilized, the average number of unserved tasks, and the average number of working days of nurses for group instances averaged for each group. Table 14 in Appendix reports these results individually for each instance of each group. In addition, Table 8 presents detailed information on average workload of the nurses during the planning horizon for each group in terms of the average number of nurses utilized, the number of tasks assigned per nurse (as min., max., and avg. over all nurses), the number of days a nurse is required to work (as min., max., and avg. over all nurses), as well as the average utilized nurse time (in column Total) and its distribution among service time and travel time (in columns Travel and Service, respectively) averaged for each group.

Table 15 in Appendix reports these results individually for each instance of each group. In Tables 7 and 8, the averages for the instances having only 1 competency level and the averages for the instances having three competency levels are also provided.

According to Table 7, the average total cost is the highest for Group-high-var instances, where the task service duration variance is larger. The larger task service duration variance in Groups-3 instances results in higher number of utilized nurses (the second largest number of utilized nurses among all groups), higher number of unserved tasks, and lower average workload level per nurse (see Table 8). Except Group-control and Group-high-var, all tasks are served in all instances. Although there is no unserved task cost incurred in Group-random-loc instances, the larger average total cost in this group of instances is due to larger recharging costs that result from random selection of locations of patients. Compared to Group-random-loc instances, the average total cost of Group-control instances is lower due to the closer locations of patients and nurses (see the share of travel in total workload of utilized nurses in Table 8). In Group-random-loc instances where the patient and nurse locations are randomly generated, the time spent on traveling is largest (see Table 8). The results for Group-long-tw instances have the smallest average total cost demonstrating the negative effect of strict time windows of patients on the service costs. In Group-long-tw instances, the tasks have more flexible time windows and there is larger room for schedule and route optimization for nurses allowing for utilizing less number of nurses with larger workload efficiency (see the average workload of utilized nurses in Table 8).

According to Tables 7 and 8, less number of nurses are utilized and the working times of utilized nurses are more efficiently planned for the instances having only 1 competency level. It is also clear that average total costs of instances with only 1 competency level is much lower owing to both higher nurse occupancy and lower unserved costs compared to instances with 3 competency levels.

4.5 Analysis on the impact of ALNS components

Additional tests were performed on 24 randomly selected instances (6 small-size instances, 12 medium-size instances, and 6 large-size instances) to assess the performance of the ALNS components, such as the SNC procedure, and the destroy and repair heuristics.

Table 9 reports the average objective function value, i.e., total cost in TRY, of the ALNS with and without the SNC procedure and the average percentage deviation for each instance size group, where the average percentage deviation is calculated as Dev. (%)= {((Without SNC Obj.) /(With SNC Obj.)*100)-100}.

Table 9 shows that the SNC procedure has a positive impact on the quality of the ALNS solutions and the impact gets larger as the problem size increases, reaching nearly 32% improvement on average in total cost in large-size instances.

Table 10 reports the percentage of iterations each destroy or repair operator is utilized among the total ALNS iterations averaged for each instance size group.

According to Table 10, the random CS removal and CS insertion operators are the least selected destroy and repair operators, respectively. This result is expected as the CS visits are rare relative to the patient visits in any nurse route and there are limited number of CSs available in Ankara decreasing the number of alternative CS locations for a required CS visit. The rest of the repair operators are relatively evenly disturbed in large-size instances. Among the rest of destroy operators, the random route removal has been most utilized operator for all instance size groups, demonstrating the improvement potential of staffing level decisions on the total cost.

We have introduced the multi-period home healthcare routing problem with electric vehicles, time windows, durations, fast chargers, and competency levels. The aim of the problem is to minimize the total cost that arises during the planning horizon, including cost of utilizing healthcare nurses, recharge energy, monetary loss due to not providing service to patients, and finally utilization of transfer service for the nurses. We have developed a powerful ALNS algorithm. Our metaheuristic has specially tailored 10 step SNC intensification procedure as well as five destroy and four repair operators. To validate the ALNS, we have conducted rigors testing on small-, medium-, and large-size instances as well as sensitivity analysis in order to examine the effects of the changes in the problem variables. Comparing the results on small- and medium-size instances has shown that the ALNS yielded lower objective functions within much shorter time frames, compared to a commercial solver used to solve the mathematical model. Our ALNS is suitable for practical use.

Our analysis has shown that the cost of HHC providers could be reduced if patient visiting hours can be easily adjusted by the providers. However, in practice, patient visit times may have to be determined more precisely and in advance according to patient specific needs and could be inflexible.

Closer examination of the results shows the importance of competency levels of nurses. Lack of nurse competency has been identified as a major driving factor for large costs. One can come to conclusion that when the cost of increasing the competency level of nurses, e.g., through training, is lower than the cost of not providing a service to patients, HHC providers should consider increasing the competency level of nurses.

Nevertheless, the fleet used in this problem is a homogeneous one, when another EV with different battery capacity and consumption rate is used the results of the problem could change. Using a heterogeneous fleet composition with various EVs will undoubtedly influence the problem. Future work could make use of a heterogeneous fleet composition while assigning different model EVs to different nurses. Another possible future study might be to investigate how to provide an alternative form of service to unserved patients, such as phone or online checks for the patients wherever possible or divert task of unserved patients to another service provider.