THE MULTI-PERIOD ELECTRIC HOME HEALTHCARE ROUTEING AND SCHEDULING PROBLEM WITH NONLINEAR CHARGING

Electric vehicles (EVs) offer ideal opportunities to reduce the negative effects of transport. Demographic changes have led to an increase in the demand for home healthcare services and the transporting of healthcare staff. This paper is an attempt to introduce and address the multi-period electric home healthcare routeing and scheduling problem with nonlinear charging (E-HHCRSP-NL) to handle decisions such as allocating the EVs used by healthcare staff for clients based on their characteristics, determining the EV routes, tracking their charge state, and choosing the most convenient charging locations. The problem is solved using a newly developed hybrid evolutionary variable neighbourhood search algorithm (HE-VNS). The computational results of the newly created instances indicate the effectiveness of the proposed approach.


INTRODUCTION
Home healthcare (HHC) can be defined as the provision of a wide range of healthcare services to patients at their homes rather than in hospitals. An aging population, increasing incidences of diseases and disabilities, and patients' preferences for independent living are the key drivers that have led to the need for HHC services. This market was worth $220.67 billion in 2016, and is expected to reach $517.23 billion by 2025 [1]. Thus decision-makers look for alternative ways to decrease healthcare-related costs and to improve the quality of care. The cost of HHC services is more affordable than those associated with inpatient residential facilities and long-term care hospitals. Moreover, these services help patients to recover more easily, regain their independence, and improve the quality of their lives as a whole. However, it is estimated that the 23 organisation for economic co-operation and development (OECD) countries will have a shortage of 2.5 million nurses by 2030 [2]. In 2018, HHC workers travelled about eight billion miles, with an average travelled daily distance of almost 11 miles [3] to visit six to eight clients, and a driving time that took up nearly 25 per cent of the working time every day [4].
The transportation sector is one of the crucial parts of any economy, with an infrastructure that is directly related to competitiveness and growth [5]. Efficient transportation can decrease costs and unintended consequences in different sectors [6]. Its energy consumption leads to greater environmental degradation; and so traffic congestion, air and noise pollution, and greenhouse gas (GHG) emissions are mainly attributed to transportation activities. The use of electric vehicles (EVs) can help to reduce these side effects [7] with 93 numerous benefits over conventional vehicles, such as lower running costs, cheaper maintenance, and safety improvements. For a sustainable system, economic, environmental, and social pillars should be considered in this framework. Transportation using EVs is one of the possible solutions [8].
The present article addresses this task by introducing a new challenge -a multi-period electric home healthcare routeing and scheduling problem with nonlinear charging (E-HHCRSP-NL). This problem considers decisions such as the allocation of the EVs (nurses) to clients by considering their features and characteristics within a predetermined time slot; the determination of the route each EV has to cover each day within the planning horizon; and the determination of the charge state and charging spots for each EV. The developed model differs from the HHC literature in the following aspects: first, the sustainability concept is not ignored; in other words, healthcare staff drives EVs to perform all the assigned jobs. Thus using EVs for HHC helps to reduce the environmental impact of transportation. Existing models covered the linear charging function and/or full charging strategy. Second, as opposed to the current literature, our model considers a nonlinear charging function and partial charging strategy.

Related work
The HHCRSP is an extension of two well-known NP-hard problems; the vehicle routeing problem with time windows (VRPTW), and the nurse rostering problem (NRP). The HHCRSP covers a set of assignments as well as temporal and geographic constraints. The presence of nurses' characteristics and of clients' requirements and preferences increases the complexity of the problem. Also, the different range of requests should be carried out by well-qualified nurses, whose allocation increases the satisfaction of clients. Furthermore, gender, language, pet ownership, etc. are defined as the characteristics or preferences to be taken into account for an appropriate assignment; otherwise incompatibilities may occur. Clients tend to seek services from the same care worker; this is known as loyalty. A minimum number of different workers should therefore be allocated to the same client during the planning horizon. As in the classical VRPTW, each client and nurse has their own time frame. Nurses can serve clients within predetermined time windows, depending on the types of the contracts. Throughout the planning horizon, the requests of clients are defined as the frequency of visits -e.g., once-a-week visits or every-other-day visits. In other words, a service that is requested every day corresponds to seven visits during one week. Nurses are not allowed to work every day, and so the maximum weekly working times have to be determined and implemented. Some jobs require the presence of more than one worker to be performed. In such cases, healthcare workers should start their jobs simultaneously. All existing publications in the literature have so far considered mainly the travel costs but ignored the ecological and social criteria. Therefore new transportation technologies, such as EVs or zero-or low-emission vehicles, should be incorporated into the planning to achieve a more sustainable transportation model. For a detailed review related to HHC, see [9].
To tackle the multi-period problem, a series of exact and metaheuristic solution procedures is applied. A branch-price-and-cut algorithm [10,11] is implemented for the exact solution procedure. Metaheuristics can be classified into two groups: single solution-based and hybrid-based. Tabu search [12], adaptive large neighbourhood search [13], and variable neighbourhood search (VNS) [10] are also employed.
In order to deal with the negative effects of transportation, zero-or low-emission vehicles have received growing interest in recent years. Erdoğan and Miller-Hooks [14] developed the green vehicle routeing problem (GVRP), in which vehicles run with alternative fuels such as biodiesel, ethanol, and natural gas. Schneider et al. [15] proposed the electric VRPTW (E-VRPTW), which is an extension of the GVRP. For interested readers, related studies for the GVRP and E-VRPTW can be found in [16][17][18][19]. Although the use of EVs has good potential, two major difficulties should be solved to eliminate their technical limitations: their limited driving ranges and their long charging times. In the literature, the proposed E-VRP (or E-VRPTW) models focus on energy consumption, charging technologies, capacity of the EVs, and location of charging stations. Most of these papers consider energy consumption in relation to the distance travelled [14,15]. On the other hand, some of the papers assume that energy consumption is not constant and depends on many factors [20,21]. The battery charging strategy is a crucial challenge to determine where, when, and how EVs will recharge. In the literature, authors assume either full or partial charging strategies. In the first strategy, the battery is fully recharged with a linear charging function when an EV visits a charging station [20,22]. Battery swapping is also taken into account, in which case the charging time is a constant [14,18,23,24]. For the partial charging strategy, the duration of the charging time is based on the battery levels when an EV arrives at a charging station and then leaves it. In other words, the battery level is defined as a decision variable. While most models propose a partial charging strategy with a linear charging function [19,[25][26][27], a partial charging strategy with a nonlinear function approximation was also recently considered [28]. 94

Contribution of this article
We combine two streams of research, one on EVRPTW and the other on HHCRSP, and introduce a new problem: the multi-period electric home healthcare routeing and scheduling problem with nonlinear charging (E-HHCRSP-NL). The aim of our study is to develop a mathematical model and a solution approach based on a hybrid evolutionary variable neighbourhood search algorithm (HE-VNS) for the E-HHCRSP-NL. The contributions of this article can be summarised as follows:  It develops a mathematical model of the multi-period problem. Nonlinear charging and a partial charging strategy are also integrated into the mixed integer programming (MIP) model we present.  We introduce a new set of instances for this new problem.  An effective HE-VNS is proposed to solve the problem. The proposed new algorithm involves the evolutionary process as opposed to the periodic HHC problems in the literature. Our algorithm is a novel approach, since there has not been any study that employs such an evolutionary VNS for solving any variant of the HHCRSP. The computational experiments indicate that the proposed approach can solve small instances optimally. A sensitivity analysis is also carried out to assess the impact of the number of charging stations and the effects of fleet mix.
The remainder of this article is organised as follows. In Section 2, the MIP formulation is introduced in detail. Section 3 describes the HE-VNS algorithm to solve the problem. Section 4 provides computational experiments for the new benchmark instances. Finally, conclusions and a short summary of our work are given in Section 5.

PROBLEM DESCRIPTION AND MODEL
In this section we introduce the multi-period E-HHCRSP-NL, and then present a mathematical model formulation.

Problem definition
Let J be a set of nodes representing the jobs, N a set of EVs (i.e., nurses), F a set of charging stations, and (hk, h'k) nodes denoting the home location (depot) for EV k. While hk refers to the starting depot node for EV k, h'k means the ending depot node for the same vehicle. Let V be the union of J and F', which is a set of dummy charging stations denoting the multiple visits of F.
 be a set of arcs connecting the nodes. In addition, let C be a set of clients, each of whom can be visited multiple times during a day. Thus let T be a set of days, and the requested jobs are scheduled considering the different types of client schemes. In other words, the service frequency is determined by every-day visits, twice-a-week visits, every-other-day visits, and once-a-week visits. For example, a client demands a specific job -e.g., checking blood pressure, insulin shot -and it is necessary to visit that client every day. Each daily visit is defined as a single job. Moreover, a daily request from a client generates seven visits during a week. Thus the parameter θig represents service requests. θig=1 if a job i is required on day g; otherwise, zero.
The travel time between any two nodes i and j is represented by sij. Time windows are represented by  [ , ] ii , and the service time of job i on day g is defined by dig. In order to perform various types of jobs, a certain order or coordination is needed. Multiple EVs should thus visit the same client at the same time. Therefore a set of synchronised jobs, P, is defined. In order to perform a job, an EV has an equal or higher qualification level. The hierarchical qualification levels are represented by EVs and jobs -q'k and qirespectively. Similarly, the characteristics of EVs and jobs are represented by m'k and mi respectively. Gender, pet ownership and allergy, and smoking habits are defined as the characteristics that are considered while allocating the appropriate assignment. Loyalty, or the number of different EVs treating/visiting client l, is determined by the parameter δl. For each EV, a maximum of two consecutive work days is allowed, and the integer constant τk is used for this purpose.
The charging function is concave [29][30][31], which can be approximated employing a piecewise linear function [28,32]. Similar to Montoya, Guéret, Mendoza, and Villegas [28], we employ a similar function for the EVs, as represented in Figure 1. In this plot, parameters Ωiek and Ψiek refer to the charge level and the charging time for the break point e of charging station i for EV k, where B is the set of breakpoints of the piecewise linear approximation. The maximum energy capacity and the consumption rate of EV k are denoted as Yk and rk respectively. The energy consumption of vehicles is related to the distance between any two nodes. It is assumed that vehicles start routeing with a fully charged battery. To formulate the E-HHCRSP-NL, we introduce the following decision variables. Let xijng be equal to 1 if EV k moves directly from node i to node j on day g, and 0 otherwise. Let λlkg be equal to 1 if EV k visits client l on day g, and 0 otherwise. Let µkg be equal to 1 if EV k works on day g, and 0 otherwise. Let ting be the scheduling variable indicating the starting time at node i on route k on day g. The non-negative variables yikg and zikg track the battery level of EV k at node i on day g, and the battery level when EV k departs from charging station i on day g respectively. Let Δikg be the time spent at charging station i of EV k on day g. The charge duration Δikg depends on the remaining battery level and the length of staying time at the charging station. Thus the non-negative variables Фikg and Гikg are defined to calculate Δikg. While variable Фikg is the battery charging time when EV k departs from charging station i on day g, Гikg is the battery charging time when EV k arrives at charging station i on day g. Variables σiekg and ξiekg are the coefficient of the break point e in the piecewise linear approximation, when EV k arrives at and departs from the station i on day g. Finally, let ςiekg and wiekg be 1 if the charge level of EV k is between Ωie,k-1 and Ωiek, when EV k arrives at and departs from station i on day g. The notations used in the following mathematical model are summarised in Table 1.

Example
A numerical example of the medium-term E-HHCRSP-NL is shown in Figure 2, which represents a solution to an instance with 10 jobs (a synchronous job), three EVs (two different types), and two stations. While the squares represent the homes, the circles are jobs and stations. Routeing and charging decisions are illustrated in this example. Each line shows different days. The starting and ending node of each route is the same depot (home location). In order to perform synchronous Job 2, EVs 1 and 2 start simultaneously. Four types of client schemes are defined as the visiting patterns or healthcare service frequency. These are as follows: every-day visits (i.e., Job 1), every-other-day visits (i.e., Job 6), twice-a-week visits (i.e., Job 2), and once-a-week visits (i.e., Job 3). As stated before, an EV is not allowed to work more than two consecutive days. For instance, EV 1 works on Days 1-2, and after two work days, EV 1 has to rest. Considering that Job 1 is requested by Client 1, this job is performed by three different EVs, and the loyalty is calculated as 2 for Client 1 throughout the planning horizon. The battery states of the vehicles are represented as the percentage values. On Day 1, Jobs 1, 2, and 4 are assigned to EV 1 (type 1). EV 1 arrives at the charging station with a 25 per cent charge level and leaves from the station with an 85 per cent partially charged battery. The charging function is presented in Figure 3; it maps the battery levels z111 and y111 to the charging times Г111 and Ф111 to compute the charge duration Δ111 (=Г111-Ф111) from the piecewise linear charging function. On Day 1, EV 2 does not need to visit any station. On Day 3, EV 2 departs from Station 2 with a fully charged battery. The same vehicle visits Station 2 twice during its route. In other words, an EV can visit a station multiple times. The MIP model follows: , , , ,    , , , , ( , ) and ( , ) , , , , ( 1) , , , , , , , , , , , , , Objective function (1) seeks to minimise the total travel time. Constraint (2) ensures that each job is taken care of on the predetermined days. Constraint (3) means that a charging station does not need to be a part of a solution. The flow conservation is guaranteed by constraint (4). The hierarchical qualification levels are defined for the EVs (nurses), and jobs 1 to 5 are taken into account by employing constraint (5). In order to perform a job, each EV has an equal or higher qualification level. Constraint (6) considers the number of different EVs visiting client l during the planning horizon. The loyalty is bounded by constraint (7). The characteristics of the EVs and the job specifications are considered by constraints (8). The working days of EVs are computed by constraint (9). Each EV is allowed to work no more than a maximum of two consecutive working days (constraint 10), in which case the next day is defined as off-duty. All nodes should be visited within the predetermined time window (constraint 11). The starting times at nodes are considered via constraints (12) and (13). Constraint (14) guarantees that two different EVs start synchronous jobs simultaneously. The state of charge is tracked by constraints (15) and (16). Constraint (17) states that the battery charge level of EVs cannot exceed the maximum battery capacity when departing from the charging station. It is obvious that the charging time is based on the energy (charge) level, which is defined by a piecewise linear approximation. When an EV arrives at a charging station, constraints (18) to (24) consider the charge level and charging time; and when it departs from the station, constraints (25) to (31) determine the charging level and time. The duration of charging is defined by constraint (32). Constraints (33) to (36) set the domains of the decision variables. Finally, (37) and (38) are non-negativity constraints.

DESCRIPTION OF THE HYBRID EVOLUTIONARY VARIABLE NEIGHBOURHOOD SEARCH ALGORITHM
This section presents a detailed description of our hybrid evolutionary variable neighbourhood search algorithm, called HE-VNS, for the periodic HHCRSP-NL. The VNS algorithm was first proposed by Mladenović and Hansen [33], and is widely used as a local search methodology due to its simplicity, efficiency, robustness, and user-friendliness. Other successful applications have also been found in scheduling [34] and vehicle routeing [35][36][37][38][39][40]. The principle of a systematic change of neighbourhood is the backbone of the VNS algorithm, which is used not only to find a local optimum, but also to escape from being stuck in the valleys that involve local minima [41]. While the local search corresponds to the deterministic phase, escaping from the local minima refers to the stochastic phase [33]. The solution quality of the VNS algorithm is based on a single solution, which is a drawback because it restricts the exploration of a larger search space. On the other hand, the HE-VNS algorithm presented in this paper searches in a given population of potential solutions -as in the population-based meta-heuristics. In other words, our algorithm enhances the VNS by means of an evolutionary strategy [42,43]. In this way, the HE-VNS explores a population of candidate solutions. Here, VNS is initialised by different solutions at each iteration, leading to a higher diversification power, in contrast with the classical VNS. Furthermore, the variable neighbourhood descent (VND) is integrated into the proposed algorithm. The VND is the deterministic version of VNS, and the neighbourhood structure is ordered to obtain a nested strategy.
The structure of the HE-VNS is represented in the algorithm. The initial population of potential solutions is generated in line 1. The fitness of each solution in the population is calculated, and the solutions are sorted in an ascending order according to their fitness (line 2). Once solutions have had their fitness evaluated, they are selected through a roulette-wheel selection (RWS) method to form the next generation in the evolution cycle (line 4). A partially mapped crossover (PMX) operator creates new offspring from the selected two parents (line 5). Each new created offspring is sent to intensification (line 7). After this procedure, each of the offspring is evaluated for its charging condition (line 8). This procedure is integrated to eliminate infeasibility resulting from either insufficient or excessive charging conditions. The fitness of each offspring is calculated in line 9. If the new solutions provide improvements, then they are accepted and inserted into the population. This replacement is carried out in line 10. During this process, infeasible solutions are also accepted as in line 1. When the number of iterations without improvement is met, the process is terminated. The population is sorted, and the current best solution is selected (line 12) and sent to the VND algorithm (line 13). While the VNS algorithm runs fewer iterations to eliminate premature convergence, the VND focuses not only on intensification, but also on infeasibility.
Algorithm : The HE-VNS algorithm 1: Generate a population using a construction heuristic 2: Calculate the fitness of each solution and sort them in ascending order 3: while (until the stopping criteria is met) do 4: Select two solutions from a population using the RWS mechanism 5: Apply crossover and form two new solutions using PMX 6: For each solution 7: Perform the VNS 8: Evaluate the charging condition 9: Calculate the fitness 10: Replacement, if the new solution provides improvement, then accept and insert in the population 11: end while 12: Sort population in an ascending order and select the current best solution xbest 13: Apply the VND algorithm 14: Return best solution xbest In order to generate an initial population, we consider a three-stage initialisation procedure. While the first stage considers the routeing decisions, the scheduling is considered in the second stage. The infeasibility is penalised in the last stage. Inıtially, the requested jobs are classified into days -that is, the number of jobs to be performed daily is determined. The algorithm assigns the requested jobs to the EVs randomly. Scheduling variables are then determined. Time windows and work day limitations (i.e., a maximum of two consecutive days) are also considered in this step. Finally, the violation of constraints is penalised in the last stage. A solution does not necessarily satisfy constraints such as qualifications, time windows, and battery capacity; these violations are converted into a penalty cost [18]. 100 Crossover is carried out between two selected solutions by exchanging the subparts of their chromosomes to yield a new generation. This procedure is a crucial component of the evolutionary algorithm, because it aids in transferring the information from the parent to create a new generation. We employ a PMX [42] that operates as follows. First, two random crossover points are selected. Then the subparts between these points are exchanged, and the remaining parts are replaced with partial mapping.
The shaking phase is the stochastic component of the VNS algorithm. When there is no improvement in the objective, a new solution is randomly generated in a given neighbourhood. A swap operator is used in this phase. A pair of nodes within a single route is exchanged without considering the objective. In order to improve the quality of the solutions, a local search technique is applied. The search operation is achieved by two operators whose order of neighbourhood is crucial in the VNS and VND, and they are ordered according to the size of the neighbourhood [41]. By means of the swap operator, the two randomly chosen nodes are first selected and then exchanged. A relocation operator selects one position and one node randomly, then the chosen node is moved to position. The randomly selected node can be a job, a depot, or a station; it is thus not necessary to define different types of operators for single-route or multi-route improvement.
In order to satisfy the charging constraint, we employ an integrated operator that involves a station insertion and removal neighbourhood structure. When a negative battery state is observed, the nearest station is inserted into the route to eliminate the violation of this constraint. Similarly, a station is removed from the route if there is no need to visit a station with a high or full battery level. This operator also determines the necessary duration of recharge at the station for each EV by calculating the battery level at the arrival depot. The charging condition is evaluated after applying the VNS (line 8), and within the VND algorithm (line 13).
This algorithm can tackle infeasible solutions during the search. As mentioned before, a penalty cost is calculated for each solution, thus guiding the search through both feasible and infeasible regions. Once any constraint is satisfied or violated, a penalty factor based on the number of iterations is updated.

COMPUTATIONAL EXPERIMENTS
In this section we present the computational results of the algorithm for solving a set of newly generated benchmarking instances. The algorithm is run on a PC with Intel Core i7-3770K, 3.5 GHz, 16GB RAM. The mathematical model presented in Section 2 is solved using CPLEX 12.8.

Data and experimental setting
In order to test the algorithm, a set of new test instances is generated, since there are no existing benchmarking instances for the E-HHCRSP-NL. We extend the instances from the single-period HHCRSP work of Hiermann, Prandtstetter, Rendl, Puchinger and Raidl [44] and Erdem and Bulkan [45] to generate a series of small-and large-scale instances from them to maintain the same ratio of the given features (time windows, the characteristics of clients and nurses, the qualification levels of EVs, job specifications). Table  A.1 in the Appendix represents the small-and large-scale instances generated. Five per cent of all jobs is considered synchronous, and is integrated into each of the generated instances. The depots, clients, and stations are distributed in an area of 400 km 2 around Vienna. The distance between any two locations is defined as a time unit, which equals five minutes. The ranges of EVs are regarded as 130 km and 180 km respectively [46]. In the same way as in the literature, here the energy consumption rate is multiplied by the distance between any two locations. As mentioned earlier, the piecewise linear function is fitted as the charging function similar to the work of Montoya, Guéret, Mendoza and Villegas [28] and Uhrig, Weiß, Suriyah and Leibfried [47]. It takes four and five hours to fully charge type 1 and type 2 EVs respectively. The service times are uniformly distributed in Grand View Research [1] and Bruglieri, Colorni and Luè [30], which are equal to the numbers between the interval of 5 and 150 minutes'. The frequency of visits is also generated randomly from the four types of client schemes -every-day visits, twice-a-week visits, everyother-day visits, and once-a-week visits -and are generated randomly while keeping the ratios of 40 per cent, 20 per cent, 20 per cent, and 20 per cent respectively. Finally, the loyalty parameter is determined as 2 in the 1-8 instances, 3 in the 9-16 instances, and 4 in the remaining instances. In order to tune the parameters, a series of computations is carried out; the selected parameters are shown in Table 2.

Computational results
The computational experiments are first conducted to assess the performance of the HE-VNS algorithm. The small-scale instances with up to 10 EVs and 111 visits are solved with both the HE-VNS and CPLEX. Table 3 presents the comparison results given by the HE-VNS and CPLEX. Here, for CPLEX, the name, the solution (Obj), and the run-time in seconds are reported. For the HE-VNS, the best solution (Best), the percentage deviation from CPLEX solution Dev (%), the average results of ten runs (Avg), and the run-time are shown. Finally, the averages of the solutions and of the run-times are given at the end of the table.
The first eight instances and instance 14 could be solved optimally. It is obvious that the first eight instances involve fewer EVs and jobs. The optimal solutions are thus found in less than five minutes. Instances 21, 23, and 24 could not be solved, and no feasible solutions were generated. The rest of the instances are shown in the third column, with an asterisk denoting that the instance was not solved optimally within one hour. Concerning the solution quality, the proposed algorithm was able to solve the small-scale instances optimally. The average of the best results of the HE-VNS is calculated as 210.1 seconds, with an average computational time of 267.5 seconds. Other experiments were carried out on the large-scale instances, with up to 50 EVs and 546 visits to evaluate the performance of the algorithm. However, the large-scale instances could not be solved by CPLEX. For each instance, Table 4 contains the best solution, the average solution, and the run-time for ten runs. The last row indicates the average of the best results and the average solutions. The average of the best results in this approach is calculated as 283. 5

Sensitivity analysis
In this section, the set of EV types is decreased to a single type in order to analyse the effect of fleet mix with different EVs. Table 5 indicates the best and average results over ten runs for each EV type. While the second column represents the best results for type 1 EVs, the seventh column refers to the best results for type 2 EVs. The third and sixth present the average objective values. When Table 5 is investigated further, it is observed that employing type 2 EVs provides improvements on average. In other words, the solutions with type 2 EVs are not only good as an objective value, but also result in a lower run-time. Since the range and the recharge duration of the EVs are considered in the proposed model, the EV type with a greater range or more battery capacity contributes more to the results.
The impact of the number of charging stations on the solutions is also investigated on randomly selected instances. For each instance, the number of stations is predetermined and fixed. Here the number of stations is changed (both increased and decreased). The best solutions and the average results are presented in Table 6. When it is investigated further, the results of instance 5 show greater improvement than the rest. As the number of stations increases, the quality of solutions initially gets slightly worse, but then improves on average. In general, when the number of stations increases, the solution improves, as expected. However, in terms of best solutions, on average the vehicles usually stop at the charging station less than once. Hiermann et al. [22] concluded that this is due to the use of a fleet mix. In other words, in good solutions, the EVs with more battery capacity or a higher range do not need to visit the charging station.

CONCLUSION
This paper introduces a new home healthcare routeing and scheduling problem using EVs. The periodic E-HHCRSP-NL considers routeing, scheduling, and charging decisions with nonlinear charging. The developed model aims to reduce the total travel time while considering a series of constraints, such as time windows, synchronisation, the characteristics of EVs (nurses), jobs specifications, competencies, loyalty, and working days. It is shown that transportation using EVs can help us to reduce the negative impacts of road transportation, thus taking into account the sustainability concept in this work. To solve the problem, we proposed a hybrid evolutionary variable neighbourhood search algorithm (HE-VNS). While diversification is achieved using an evolutionary strategy, intensification is carried out with the VNS/VND. Our proposed algorithm is able to search through both feasible and infeasible regions. Therefore, in small instances, the HE-VNS is able to yield optimal solutions. Our algorithm can also find good solutions in a reasonable amount of time. A sensitivity analysis is performed on the results to investigate both the effects on the number of charging stations and the effect of a fleet mix of different EVs. The results of our computational experiments indicate that the EVs with more battery capacity or a higher range do not need to visit the charging stations in good solutions.