## 1. INTRODUCTION

Recently, the concern about environmental issues, especially global warming and climate change has increased which appeared through the studies focusing on the green concept. It is a concept used to describe any alternatives aimed to eliminate or reduce waste whatever it is emissions, chemical, or solid waste during its implementation, therefore reduce their effects on the environment. One of the most important concepts used in the industry is the green supply chain management considering the minimization of waste along the product’s life cycle including procurement, manufacturing, distribution and reverse logistics (Ninlawan *et al*., 2010). In this work, we concentrate on a distribution stage that is because of the emissions generated by the transportation sector, in united states the significant percentage of greenhouse gas emission was sheared by transportation compared with the other economic sectors in 2017 (U.S Environmental Protection Agency, 2019). Also, according to a European Union report, the percentage of emissions generated by main sources such as energy industries and manufacturing was reduced from 1990 to 2016 except that was shared by transport, which road transport represent almost threequarters of it, increased by 30% (Climate change, 2016).

There are many studies focused on improving the logistics and distribution system to achieve their goals in line with reducing their impact on the environment, win- win strategy, the Green Vehicle Routing Problem (GVRP) was introduced to minimize the emission and fuel consumption (Park and Chae, 2014). In parallel with that, the EVs have appeared as an alternative to conventional vehicles that utilized fossil fuel.

According to International Energy Agency (IEA, 2017) the number of registered EV increased in the past few years, approximately reached 750 thousand sales in 2016, that shows the growth in environmental concerns and the major role of EV in reducing the pollution. Furthermore, as result of development EV's technology the usage of EV not limited to private transportation, they were used in public transportation in different cities such as Germany, USA, China (Hu *et al*., 2018), Recently Jordan started using it and they can be used in different distribution operation (Keskin and Çatay, 2016). Also, some major distrusted company such as FedEx have started adding electric vans to their fleet.

Even the EVs do not produce pollutants; they have limitations such as limited driving range which means vehicle needs to recharge along the route, and they are recharged by electrical energy that is sometimes generated by non-clean methodologies, as in Jordan for example. Some of these non-clean methodologies are by using nuclear plant, natural gas, and coal which represents around 37% of the whole generated electrical energy worldwide (World Coal Association, 2018). Accordingly, utilizing electric vehicles in an optimized manner will reduce the number of recharge times and consequently reduce the consumption of non-clean electrical energy.

The remainder of the paper is organized as follow. Section two represents a review of relevant literature. The problem description is presented in section three including the discussion of the mixed integer linear programming formulation for the model. In section four, a numerical example is presented. The structure of a solution approach is discussed in section five. Section six describes the numerical experiments to evaluate the performance of the proposed heuristic. Section seven provides the conclusion and future work. Finally, a real case study is discussed in the last section.

## 2. RELATED LITERATURE

In this section, we review literature that studied the Electric Vehicle Routing Problem (EVRP). Firstly, the vehicle routing problem is an optimization problem to determine the routes for vehicle fleet serving certain customers with given demand and optimizing one or more objectives (Laporte *et al*., 1992). The EVRP is one of the recent types of VRP in which a fleet of EVs is used.

Erdoğan and Miller-Hooks (2012) introduced a model of Green Vehicle Routing Problem (GVRP) which is, to the best of our knowledge, the first model involving limitation in a fueling capacity of the vehicle and in the availability of fuel station. The objective function of the model is a minimization of the total distance traveled by a vehicle with a given consumption rate. Two heuristics proposed to solve a real-world problem instance the Modified Clarke and Wright Saving heuristic (MCWS) and the Density-Based Clustering Algorithm (DBCA). These algorithms provided results with difference percentages ranging form 0% to 11.5% compared to the CPLEX results in small instances, where the percentages rise up to 22% when clustered customers are used in larger instances. The similarity of complication between this problem and EVRP inspired researchers to extend Erdoğan and Miller-Hooks (2012) works in different ways, such as Schneider *et al*. (2014) and Felipe *et al*. (2014). Schneider *et al*. (2014) introduced the Electric Vehicle Routing Problem with Time Window and Recharge Station (EVRPTW) with the possibility to recharge the vehicle to fully charge within a time depended to a vehicle charge level at the station, also vehicle capacity constraint and customer time window were interpolated into the model. Schneider *et al*. (2014) proposed a hybrid heuristic, a Variable Neighborhood Search (VNS) algorithm combined with Tabu Search heuristic, which provides results with gap between 0% to slightly less than 3%.

Similarly to previous studies, Felipe *et al*. (2014) approached a model of Green Vehicle Routing Problem with Multiple Technologies and Partial Recharges (GVRPMTPR). In this model, the possibility of partial recharge by multiple technologies was interpolated which made the problem more realistic and saving in potential cost. Moreover, a vehicle could recharge in a depot. The model objective function was a minimization of total cost. Felipe *et al*. (2014) used a constructive and local search heuristic to solve the problem. Hiermann *et al*. (2016) developed a model for the Electric Fleet Size and Mix Vehicle with Time Windows and Recharging Stations (E-FSMFTW) which aim to minimize acquisition cost and the total distance traveled. According to this model, the mix vehicles which are different in capacity and cost are used. Hiermann *et al*. (2016) proposed an Adaptive Large Neighbourhood Search (ALNS) with local search and labelling procedures. This hybrid heuristic was able to solve the instances with deviation less than 3% compared with branch and price.

In the same year, Keskin and Çatay (2016) adapted the E-VRPTW in (Schneider *et al*., 2014) and allowed partial recharging at stations to introduce the EVRPTWPR. ALNS approach was proposed to solve the problem. Results showed improvement in the solution quality up to 3% compared to the best known solutions. Bruglieri *et al*. (2015) solved the EVRPTW-PR by using a Variable Neighbourhood Search Branching (VNSB). The objective was to minimize the total travel, waiting, and recharging times. The suggested VNSB was able to solve the pro posed models with gap less than 1% and less computational efforts compared with CPLEX results. As an extension to previous research, Bruglieri *et al*. (2017) proposed a Three-Phases Matheuristic (TPM) to solve the problem. According to the result on some benchmark instances, TPM was able to provide better solutions better than VNSB, which was proposed in Bruglieri *et al*. (2015). More recently, Zhenfeng *et al*. (2017) considered driving cost and punishment cost for any delay on the customer time window to develop the EVRPTW model. Furthermore, a genetic algorithm was proposed to determine the solutions. The proposed algorithm was used to solve some small instances. Zhang *et al*. (2018) developed the EVP model to minimize energy consumption as an objective function instead of the classical function which is minimizing total travelled distance. Energy consumption was represented in the model as the function of vehicle weight, cargo weight, speed, and road friction. Also, they determined the effectiveness of that by solving the problem of both the distance minimizing and the energy minimizing function for the same instants then compared the total energy consumption of two objective functions. The result showed that the classical function had a high energy consumption by 16.44% on average. Moreover, they tested the performance of an ant colony (AC) algorithm as a solution method by comparing it with an ALNS on small and large size problem instances. In small size problem instances, the AC solution was the closest to the optimal solution, also it was better in quality and computing time in large size problem instances. Keskin and Çatay (2016) used the same strategy to illustrate the effeteness of combining an ALNS with an exact method and considering fast recharge at recharge stations to find the optimal solution which has a minim energy cost for the Electric Vehicle Routing Problem with Time Windows and Fast Chargers (EVRPTW-FC). The result of small size problem instances showed that the proposed method outperformed CPLEX in both quality and computing time. Whereas the result of large size problem instances showed the benefit of using fast recharge in reducing the number of needed vehicle and reducing the total energy cost. Ferro *et al*. (2018) introduced a new EVRP model considering different consumption and charge mode, which are selected depending on vehicle weight and speed and the battery charge at recharging node. The objective function is to minimize the cost associated with the total travelling distance and cost of recharging energy depending on the selected recharge mode.

The proposed model in this study is a mixed integer linear programming model that builds on the classical NP-hard vehicle routing problem with incorporation of other distribution decisions, avoiding traveling electric vehicles while they are fully loaded, and recharging them at proper location and time. This, up to our knowledge, has not been considered in the previous studies. Additionally, a new heuristic, different than the approaches that were used in literature, will be constructed to solve the resulting model in reasonable computational effort based on solving the model iteratively and optimally in the last phase according to some specific variable values which are determined greedily in previous phases. Table 1 shows the studies with respect to different EV criteria with a comparison with this study.

## 3. PROBLEM DESCRIPTION AND MATHEMATICAL FORMULATION

The distribution system in this problem consists of a single manufacturing plant, one single depot, and a heterogeneous fleet of EVs represented by the different loading and battery capacity. The vehicle departs from a depot with battery charged at any quantity then the battery level starts decreasing proportionally with travel time and the vehicle load, also it may visit the recharging station to be able to continue its tour. The first decision that should make is what vehicles leave the depot and the quantity to fill these vehicles. After that, vehicles move to customer's nodes to unload their loads, and if needed, they come back to the depot to reload.

The type of delivery used in the distribution is a split delivery which differs than other VRP that assuming each customer should be visited by exactly one vehicle; thus, a customer demand should not exceed the vehicle capacity. In the split delivery, a customer could be serviced by more than one vehicle, which may reduce the number of the vehicle utilized, the total travelled distance and improve the overall efficiency of a distribution system by reducing the vehicle’s idle time.

The objective function of the proposed model is to minimize the unsatisfied demand and the total traveling time and to eliminate unneeded movements for each vehicle as much as possible. Where, the unmet demand can be caused by some factors in the logistic system, such as time horizon limitation or lack of manufacturing productivity. Moreover, the model could be used to minimize the total cost of the consumption energy indirectly by multiplying the unit cost of using the vehicle by the overall distribution time for each vehicle.

### 3.1 Notation – Parameters

Let set *N* represents all vertices in the network; where *F* represents the charging stations node, *d *represents the depot node, M represents the set of customers node and *i*′ represents the dummy node that is used to facilitate the mathematical modelling, as it will be discussed later in the article. Since the same vehicle may visit the recharge station multiple times or different vehicles may visit it at the same time, *L _{f}* is defined to denote the capacity of each recharging station

*f*∈

*F*.

The time horizon T consists of discrete integer periods with equal length, *T* = {1, 2, …, |*T*|}. The parameters related to customer *i*∈*M* are *q _{i}* and

*P*, representing customer demand and customer priority, respectively. Each vehicle in set

_{i}*V*is associated with capacity

*C*and battery capacity

_{v}*Q*also it leaves the depot with battery charged to

_{v}*y*level.

_{v}Furthermore, the energy is consumed at a rate of *R* per unit of time, per unit of demand. *G* denotes the recharging rate. The travel time between node *i* ∈*N* to node *j* ∈ *N* is represented by *τ _{ijv}*. Finally,

*S*denotes the available amount in the depot. Table 2 summarize the sets and parameters.

_{d}### 3.2 Notation – Decision Variables

In this section, decision variables are presented, which are required to give more details about the solution obtained by the proposed model. There are two types of variables, integer and binary variable. Where, the integer variables are ${Q}_{vt}^{P},{Q}_{ivt}^{D},{D}_{i},M{R}_{v},{Z}_{vt},{A}_{vt}$ representing pickup, delivery, deviation, maximum return time, battery charge, actual charging rate variables, respectively, and the binary variables are the routing variables *x _{ijvt}*. Table 3 summarizes these variables.

### 3.3 Mathematical Model

The MILP in this problem including constraints and objective function is described in this section. The constraints are divided to three categories, routing, pickup and delivery, and charging constraint. The objective function will be described after introduced the constraints.

#### 3.3.1 Routing Constraints

The objective function subject to the routing constraints to manage and control the vehicle's movements from leaving the depot until return to the depot at the end of the route. These constraints are presented below:

The first constraint (1) is required to make sure that each vehicle *v* should be in one place at each time *t *, which means it cannot do more than one job at the same time. Constraints (2) - (5) manage the movements from and to the depot. Specifically, constraint (2) ensures that each vehicle *v* starts from the dummy node *i*′ and departures from it just once in all time periods, as well as, constraint (3) imposes that the vehicle returns to *i*′ once in all time periods. Constraints (4) and (5) are used to determine the first and the last movement for each vehicle v assigned to routes. Particularly, constraint (4) imposes that each vehicle v should start its route by traveling from node *i*′ to the depot at the beginning of its time horizon, and constraint ((5)) imposes that it should end its route by traveling from the depot to node *i*′ at the end of its time horizon. Depot is not directly used to start the vehicles’ tour, but the dummy node is used for this purpose to satisfy the constraints (6) – (7), as will be explained below.

Constraints (6) and (7) are continuity constraints which are required to keep the continuity of routes for each vehicle *v*. Specifically, constraint (6) states that each vehicle *v* travels from node *j* to node *k* only if it comes from another node *i*, where the time needed to reach the node *j* less or equal the current time (t) minus the travel time between *j* and *k*, for ease to understand suppose the time, when the vehicle reaches node *k*, equal (t = 9) and *τ _{jkv}* = 4. Then, it should reach node

*j*coming from

*i*at time $r\le (t-{\tau}_{jkv})$ ,

*r*=5 . Constraint (7) states that each vehicle

*υ*reaches node

*j*, should leave it to another node.

There are two types of vehicle movements related to the depot, the first type of movement is when the vehicle leaves the depot at the beginning of time horizon, and the second type is when the vehicle returns to the depot to be re-supplied and leaves again to perform other distribution plan. Using the depot without the dummy node is not allowed, as Constraints (6-7) do not allow the first type of movement. This is because when the vehicle leaves the depot at the beginning, as it does not come from other node. Thus, the dummy node is added to allow these constraints to be satisfied.

To eliminate the effect of this added dummy node, Constraints (2-5) enforce all vehicles to leave the dummy to one direction only; which is toward the real depot, and when the travel time between the dummy and depot is zero. Dummy node can be removed only if the vehicles are not allowed to re-visit the depot and i’ is replaced by d in constraints (6-7).

Constraint (8), presented below, is used to ensure that each vehicle *v* returns to the depot before the time horizon *T *ends, with no need to consume all *T*. Note that each vehicle can visits the depot more than one time during *T* to reload.

Constraint (9), presented below, is required to define the maximum time that vehicle *v* need to return to the depot at the end of its tour.

As mentioned before, the recharging station can be visited by the same vehicle more than one time, or by different vehicles within the time horizon *T*. Therefore, constraint (10) enforces that the number of visits to each recharging station *f* ∈ *F*, dose not exceed its capacity *L _{f}*.

#### 3.3.2 Pickup and Delivery Constraints

In this subsection, the pickup and delivery constraints are presented, which are used to manage the amounts picked up or delivered by the vehicles within the time horizon *T*. These constraints are shown below:

Constraint (11) enforces that the total amount delivered by all vehicles to all customers at the end of *T* should not exceed the available amount at the depot. Constraint (12) is required to define the deviation in demand for each *j *∈*M *, which equals the difference between the customer demand and the total amount delivered to customer *j *by all vehicles at the end of time *T *. As mention before, the reasons that may lead to the deviation in demand are the limitation in time horizon or lack of manufacturing productivity. It is worth to mention that these variables, *D _{i}*, can be eliminated from the model, as they strict equal to the right hand side of Constraints (12) and they appear only in the objective function. But if they are eliminated, the objective function should be updated by including the right hand side of Constraint (12), and this term should be explained as deviations again. Thus, using Di variables is to improve the understandability of the model.

Constraints (13) and (14) are required to balance the flow of amount picked up and delivered between nodes, within the time horizon *T *. Specifically, constraint (13) is used to ensure that the total picked up amount loaded into each vehicle *v* up to time *s*≤*t*−1 should be greater than the total amount delivered to all customers by the time *s *≤*t *. Constraint (14) ensures that the total amount picked up by each vehicle *v *should be equal the total amount delivered to all customers at the end of time horizon *T *.

Constraint (15), presented below, enforces that the difference between the total picked up amount and the delivered ones to all customers by each vehicle *v *up to time *t *does not exceed the vehicle capacity, where s≤ t .

Finally, constraint (16) and (17) are called big M constraints, which are required to ensure that each vehicle *v* dose not deliver or pick up any amount from any node unless it visits that node:

Constraint (16) ensures that each vehicle *v* doesdeliver any amount to any customer at time *t* unless it reaches that customer at the same time from any node. *M*_{1} in the constraint is a big number, and it may define as the demand of each customer $({M}_{1}={q}_{k}\hspace{0.17em}\forall k\in M)$, which is a maximum amount that the vehicle may deliver to the customer. Similarly, constraint (17) is required to enforces that vehicle *v* can pick up amount from the depot only if it reaches the depot. The *M*_{2} value in the constraint equal the maximum capacity of vehicle *v*$\left({M}_{2}={C}_{v}\hspace{0.17em}\forall v\in V\right)$.

#### 3.3.3 Charging Constraints

The final set of constraints are used to determine the battery charge level for each vehicle at each *t* ∈ *T*. These constraints are presented below: Constraint (18) implies that each vehicle departs from depot at the beginning of time horizon with level power equal ** y_{v}**, which means that it is recharged before departing the depot to any level of power could be less than or equal to battery capacity

**.**

*E*_{v}

Constraint (19) determines the actual recharging rate of each vehicle at recharging station at time *t*, where the partial recharge is considered at the recharging station. Constraint (20) is used to determine the battery charge of vehicle *v *at time *t*, which equals the previous battery charge plus the consumption rate per unit multiplied by the difference between the picked up and delivered amount plus the charging rate up to this time minus the consumption rate of vehicle weight. Without using the variables *A _{vt}*, vehicle can go to a charging station only if the remaining energy in the battery plus the full charging rate equal to the battery capacity, and the partial charging rate will not be allowed. Let’s explain this by an example, if the remaining energy in a vehicle battery is 20%, and the charging rate at a station is enough to charge 90% of this battery in a time period, then, without using

*A*, this vehicle can not go to the station until it consumes 10% more.

_{vt}With using *A _{vt}*, the vehicle can go to the charging station and

*A*can take a value of 80%. The other benefit of using such variable, is to increase the efficiency of vehicles, as any vehicle, if it has free or ideal time, it can get some energy charged to be ready for any new distribution plans without need to wait until it consumes most of its energy load.

_{vt}Constraint (21) ensures that the power level of vehicle ** v** at each time should not exceed the battery capacity

**.**

*E*_{v}#### 3.3.4 Objective Function

The objective function of this EVRP almost follows the same objective used in (Al Theeb *et al*., 2019), which is to determine the optimal route for each vehicle which ensures the customer’s satisfaction by delivering demands based on the customer's priority and has a minimum total time and cost. The formulation of this objective is:

There are three terms in this formulation, the first term is the main one, which represents the customer’s satisfaction by minimizing the deviation of demand *D _{i}* ordered by customer

*i*based on his priority as much as possible. The definition of

*D*is detailed in constraint (11). The second term is to minimize the time that vehicle

_{i}*v*needed to return to the depot at the end of its route.

*MR*is detailed in constraint (9), which ensure vehicle

_{v}*υ*not do any unnecessary movement during the route. The third term is to prevent any extra tours of vehicle

*v*. The second and third terms are supporter terms. They have less importance than the first one and to ensure that

*α*and

*β*must be small number depending on the problem size.

## 4. HEURISTIC SOLUTION APPROACH

As the proposed model is NP-hard, it can not be solved to optimality in reasonable computational efforts for large scale problems. In this section, the Tabu Search with Adaptive Random Search (TSARS) is proposed to solve the different size of data sets. It is required to obtain solutions with high quality in reasonable computational time compared to the time needed by the exact method, which is clarified in the next chapter. The basic concept of Tabu Search (TS) is described in (Glover, 1986). Gendreau *et al*. (1994) are the first who use TS to solve VRP. Moreover, the development of TS has been widely used on optimization problem in many applications involving integer programming (Wang *et al*., 2018), scheduling (Li and Gao, 2016), and traveling sales-man (Lin *et al*., 2016). Regard the using of Adaptive Random Search (ARS), many authors have used this technique to solve NP-hard models, both (Al Theeb and Murray, 2017) and (Al Theeb *et al*., 2020) have used this approach which is the closest to the ARS used in this research. The main advantage of the solution approach, that is used in this research, is to combine both ARS and Tabu search to create TSARS which shows better results compared to the using ARS alone.

The proposed heuristic passes through three main phases. In the first phase, routes are constructed based on greedy random search and feasible solutions are generated. Then, in the second phase, neighbours are generated by using swapping techniques. Finally, the third phase utilizes the tabu search.

### 4.1 Generate Solutions and Neighbours

The heuristic approach begins by applying the solution generation function, as in Algorithm 1. This is achieved by defining the updated demand (*UD _{j}*), which equals to the actual demand at the beginning, and the importance of each customer (

*Imp*)

_{i}*i*∈

*M*as in equations 23 and 24. Then, these nodes are arranged in descending order based on their importance value

*Imp*. Parameter

_{i}*H*in equation 24 is a uniform random number generated within a range of [0.8, 1.2]. The benefit of using

*H*is to add randomness to the calculations and to get slightly different results in each iteration. Next, the vehicle loop starts considering the maximum number of visits for the vehicle

*v*before it returns to the

*i*′ at the end of its tour or to the depot to reload again, this number is represented by

*dec*depending on the vehicle capacity and average customers demand as in equation 24. The parameter

*dec*is just an approximate estimate for the number of visits and does not mean that demand is evenly distributed to the visited customers. Note that, the first movement of the vehicle should be from the depot. Therefore, it is cited as the

*current*location for

*v*as in line 7.

Inside the *v* loop, the time loop begins for such vehicle considering two conditions. The first one is to test if the c*ounter* is less than or equal the *dec* value, it indicates that such vehicle can transfer from the current node to the selected customer node *i*∈*arranged _{i}*. The second condition is to test the feasibility of visiting this node, such as testing if this node feasible to be visited before the time horizon ends. If these conditions are satisfied, then the vehicle travels from current node to the highest importance customer node. Therefore, ${x}_{ijvt}^{new}$ is updated to 1, not that new new ${x}_{ijvt}^{new}={x}_{divt}^{new}$ for the first route of vehicle

*v*. If these conditions are not satisfied, this means that the tloop will do nothing at this value of

*t*, as the vehicle needs longer time to reach the node

*i*at

*t*. So, we might find that, for example, there are a visit at to node

*i*at

*t*= 2, and then at other visit to node

*j*at

*t*= 4, but not

*t*= 3, as the vehicle needs 2 time units to reach from

*i*to

*j*, and at

*t*= 3 it will be in the road between them.

When the vehicle reaches the customer node, this implies that such customer received all or some of his demand. Therefore, the updated demand for such customer should be dropped down, as in equation 26, also, the importance of such customer node is decreased or dropped down, so *Imp _{i}* should be recalculated to avoid being revisited by other vehicles in the same period unless its importance excuses that.

After that, customer nodes are rearranged based on the *Imp _{i}* value in descending order, and the number of visits

*counter*is incremented by 1 after each visit. The whole process is repeated for the same vehicle until the

*counter*value becomes greater or equal than the

*dec*value, which means that such vehicle cannot travel to any customer node. This is because such vehicle battery does not have enough charge to visit the customer node or the such vehicle should return to the depot to reload again in case of the determined visits have not been finished yet.

The vehicle *v* is assigned to visit recharging station *i *∈ *F* to recharge its battery, as in line 19, since the *counter* value is greater than or equal to the *dec *value. Then such visit is tested by feasibility condition. This condition is to test if the vehicle is feasible to visit the recharging station before the time horizon ends by checking the inequality now ${t}^{now}+{\tau}_{currentiv}\le \left|T\right|$. If this condition is satisfied, the visit for the recharging station is accomplished. Then, ${x}_{ijvt}^{new}$ value updated to 1, and the current location is updated, where *current* =*i*; *i* ∈ *F* . Further, the importance value for each customer (*Imp _{i}*) should be recalculated again, as in equation (24), because of the travel time between the last selected node and the recharging node is changed. After that, the customer nodes are rearranged again in descending order based on the

*Imp*value. The algorithm can record the stopping times while the times during the travel are not recorded and appear as doing nothing inside the for-loop. Hence, if both ifstatements in lines 10 and l9 are not satisfied and the counter is not updated, the for loop will do nothing until the if-statement in line 10 is satisfied. Doing nothing means that the vehicle travelling in the road during.

_{i}Once the recharging station is visited, if the battery capacity (*E _{v}*) of vehicle

*v*at each time is greater than the recharging rate

*G*by 2 or equal it, as line 27, this indicates that such vehicle should spend extra time period at the recharging station. If it is less than 2, as in line 30, then such vehicle departs the recharging station without spending any extra time. However, this concept can be easily changed and vehicles can spend more time periods if the time needed to charge is greater than one extra time period by changing the if statements. If the vehicle fails to make the visit to charging station at all the remaining time periods, this means that the remaining time is limited and it returns to the depot.

The whole procedure is repeated until all vehicles are selected. Binary variables are known and set to zero or ones for each vehicle at the end of the time horizon, if the value of the binary variable ${x}_{ijvt}^{new}$ equal zero, then the value of original binary variable *x _{ijvt}* is added as constraint to the model, such that (

*x*= 0 ∀

_{ijvt}*i*,

*j*∈

*N*,

*i*&

*j*≠

*i*′ ,

*v*∈

*V*,

*t*∈

*T*). while, if the value of the binary variable new ${x}_{ijvt}^{new}$ equal one, then the value of original binary variable

*x*is added as constraint to the model, such that (

_{ijvt}*x*= 1 ∀

_{ijvt}*i*,

*j*∈

*N*,

*i*&

*j*≠

*i*′#x2032;,

*v*∈

*V*,

*t*∈

*T*).

After all binary variables are set to zeroes or ones, the model is solved to determine the total delivered and picked-up amounts for each vehicle *v*. It is important to mention that the quantity delivered to each customer is determined by the exact solution provided by CPLEX after the binary variables are fixed to zeros and ones. The updated demand in Equation 26 helps in determining the length of tours while the routes are constructed and not necessarily that the CPLEX will produce the same approximate delivered demand when solve the model.

After the model is exactly solved by CPLEX at fixed binary variables, the function, named “generate neighbours” is applied to generate other solutions, as in Algorithm 2. The delivered amount by each vehicle can be calculated, as per equation 27 in line 1. Then, select two different binary variables randomly, representing the successful visit node for example ${x}_{ijvt}\hspace{0.17em}\hspace{0.17em}{x}_{{i}^{\prime}{j}^{\prime}{v}^{\prime}{t}^{\prime}}$, if it feasible to swap between them, such as (${x}_{i{j}^{\prime}vt},\hspace{0.17em}\hspace{0.17em}{x}_{{i}^{\prime}j{v}^{\prime}{t}^{\prime}}$).

### 4.2 Tabu Search Phase

To perform the tabu search, the function of neighbour generation is repeated *n* times to find a new route for vehicle *υ*, and Algorithm 3 is applied. In this proposed algorithm 3, the best candidate is initially considered as the solution produced by the end of phase 1, neighbourhood is defined as any solution can be reached by swapping two nodes in the route as described in phase 2, Tabu list is initially empty and is filled by the solutions that are worse more than the current best candidate as in Line 8, size of tabu list is 20 routes for each vehicle, aspiration criterion is an optional phase in TS and it is not defined in this algorithm. This is because any solution in Tabu list will not be attractive to be accepted as the objective value is used in this algorithm to judge the acceptance of solutions, and as the solutions in Tabu list have objective values worse than other solutions; the tabu listed solutions will not be considered anymore.

After the neighbour route is found, this route is tested according to tabu search by a condition which is used to check if the route in the tabu list. If this condition is not satisfied, then the model should be resolved to find another complete solution with the whole delivered amounts and to find objective function value. Then, the Tabu list might be updated as it will contain all routes for each vehicle that their new delivered amount less than the current one (worse), which means that such vehicle delivers less amount if it takes the new route. Finally, the whole solving procedure is repeated until the termination criterion, this implies that until the best solution found, which achieves the objective function of the model is selected and saved as the final solution. Termination criteria are selected as follows; the main algorithm number 1 is repeated 50 iterations, Algorithm 2 is repeated 3 times the number of nodes in each vehicle’s route, and Algorithm 3 is applied after finishing Algorithm 2. Both Algorithms 2 and 3 are repeated 20 iterations (*n* in Algorithm 3) inside each iteration of the main Algorithm 1.

## 5. NUMERICAL ANALYSIS

In this section, different size data sets are solved by CPLEX within a limited time of three hours to test the performance of the proposed model. Then they are solved by the proposed heuristic explained in Section four to test its performance by comparing the results produced from it with the results obtained by CPLEX in term of solution quality and computation time.

The size of data sets is small, medium, and large classified based on previous studies. Each one of them contains 20 generated data sets to test the performance of the proposed model and the heuristic approach. The timelimited based on an actual situation is three hours. The experiments are implemented on an HP-Tower workstation Z240 with a Xeon-Intel process running Ubuntu Linux 16.04.1 LTS in 64-bit mode.

The input parameters are randomly generated using a C++ programming language with CPLEX Concert Technology version 12.6.3 and follow a uniform distribution. The value of parameters that have a significant effect on the computation time follow a different uniform distribution depending on the size of the data set. Table 4 shows these input parameters.

Other parameters are considered as fixed parameters, which have less effects on computation time and are independent of the problem size. These parameters are the capacity of each vehicle, which is generated in the range of unif~(1000, 2000), the demand requested by each cus tomer, which is selected in the range of unif~(600, 1400), the customer priority which has the rang of unif~(1, 4), and the travel times between nodes, which selected in the range of unif~(1, 3).

In the next sub section, the different data size sets are solved by the proposed heuristic approach. Then the generated results are compared with that obtained by CPLEX in term of solution quality and computation time to test the performance of the solution approach.

### 5.1 Small Size Problems

The small size data sets are solved by using the proposed solution approach. These sets are used to test the optimal results produced by CPLEX versus heuristic. The gap can be defined as a difference between the objective function value (OFV) generated from CPLEX solution, and the OFV resulted from the solution approach; it is calculated by the following formulation:

Figure 1 shows the comparison between the performance of CPLEX and the proposed heuristic in terms of the gap and computation time. As shown in the figure, the proposed heuristic solved 16 sets of 20 with zero gaps, which means they have the same OFV of what the CPLEX found, while the other 3 data sets are solved with average gap equals 4%. In term of computation time, it can be noticed that CPLEX solves the small data sets within the average time less than what heuristic need but saved time could be negligible.

### 5.2 Medium Size Problems

The medium size data sets worked by CPLEX are solved by the proposed heuristic method. Equation 27 is used to calculate the gap. Figure 2 shows the comparison between the results of both methods CPLEX and proposed heuristic. From results, the proposed heuristic obtained the same or better results compared to CPLEX for 15 of 20 data sets, which means that both CPLEX and heuristic approach provide competitive results in term of solution quality, while the heuristic approach outperforms CPLEX in term of computation time, it consumes on average 102 seconds, which represents 1.5% of the time consumed by CPLEX to solve the same data sets.

### 5.3 Large Size Problems

For large sets, CPLEX consuming all allowable time in average to find a feasible solution for five out of twenty large instances. In this sub-section, these data instances are solved by the proposed Heuristic. Equation 27 is used to calculate the gap. Figure 3 shows the comparison between the performance of CPLEX and the proposed heuristic in terms of the gap and computation time. As seen, the proposed heuristic is able to produce solutions for all large data sets. The gap is not applicable for 15 data sets, although the heuristic found a feasible solution for all data sets, this is due to the inability of CPLEX to find a feasible solution for these data sets, as previously explained. Thus, they are excluded from Figure 3. The proposed heuristic consumes 0.81 hour on average, which represents 20% of what CPLEX consumed. Also, there are three instances have a negative gap, which means the heuristic was able to find a better solution for them. From this section, it can conclude that the proposed heuristic approach outperforms CPLEX for solving many instances in term of quality and computational time. Moreover, both CPLEX and the heuristic solution become unable to find the optimal solution when the size of data sets increases.

At the end of results section, the convergence curve is shown below in Figure 4 for a random selected large scale data set. The aim of this figure is to show the improvement of objective function values versus the iterations to figure out that the improvement is achieved by a high percentage of iterations.

## 6. CASE STUDY

Farm dairy company was established in 1994 in the capital of Jordan, Amman, specialized in the production of dairy products such as cheese and yogurt. It distributes the products to hotels, restaurants, malls, and markets throughout the kingdom. According to the priority policy of the company, the specific vehicles are assigned to serve the customers based on their priority, where hotels have a high priority, then restaurants and mall, and markets in the last.

The distribution system of the company inside Am-man is selected to define it as the problem because it covers the largest number of customers who consume the largest amount of daily production. Figure 5 shows the distribution of nodes throughout Amman, where a big red sign, small red signs, green signs represent the location of the depot, the recharging stations, and the customers, respectively. Noted that the data collected for a specific day. As CPLEX fails to find an optimal solution for large size problems in a reasonable time limit, the clustering method is used to divide the nodes to clusters based on the demand and how close the nodes from each other's, where each cluster cover approximately 25% of the total demand to maintain the realistic of the problem. That could help the CPLEX to find the optimal solution or close to optimal in a reasonable time for each cluster, where one vehicle is assigned to serve the customers inside each cluster.

First, the problem is directly solved by using the model as EVRP, where the electric vehicles are used. Then it is solved as VRP by the same model, where the conventional vehicles are used, but the starting charge y_v for each vehicle should redefine to value ensuring that the vehicle does not need to visit the recharging station during its tour which means it is operated as old fuelled vehicles. The comparison between the two problems is presented in the next subsection.

### 6.1 Results of Case Study

CPLEX is used to solve the problem for each cluster, assuming that the electric and conventional vehicle delivers the same amount for each cluster and travels with constant speed equals 40km/hr, to ensure the effectiveness of the comparison. The travel cost for each vehicle is calculated by equation 29, which equals the total distance travelled by each vehicle multiplying by the consumption rate and the unit cost. Note that the consumption rate for the conventional vehicle is Litre per km, while is kWh per km for the electrical one. According to the Ministry of Energy and Mineral Resources in Jordan, the unit cost is 0.59 JD/ L, 0.135 JD/ kWh for the conventional and electrical vehicle, respectively (IMDb, 2019;Memr.gov.jo, 2019).

The company uses the Hyundai HD truck to distributes the products, where the consumption rate for this type of vehicle is 0.175 L/ km, while the consumption rate of the electrical vehicle is 0.48 kWh/ km. The total travel cost for each vehicle *υ* is presented in Figure 6. As seen, it can be concluded that using the electric vehicle instead of the conventional one provides saving in total travelled cost by 11%, 9%, 14%, 12% for each vehicle, respectively. Where the company can save 12.39 JD for that day, approximately 322.14 JD in month throughout Amman only, noted that the workdays of the company in week is 6 days. Moreover, the CPLEX can find the optimal or close to optimal solution for each cluster consuming 2.3 days.

## 7. CONCLUSIONS

In this research, an optimization model for EVRP is presented. The concept of a heterogeneous electric vehicle with different load capacity and battery charge capacity used to transfer the commodities from the depot of a company to multiple customers using the split delivery system which allows vehicles to return to the depot to reload if needed is concerned. The proposed MILP determines the optimal route for each vehicle that minimizes the total amount of unmet demand and minimize the overall distribution time.

CPLEX is used as an exact method to solve the different size of data sets to test the performance of the proposed model. According to the result as the size of data set increased the CPLEX become unable to find the optimal solution within the limited time. There for a new heuristic approach developed to solve, which obtains high quality solution within less time than CPLEX needed.

Numerical analysis is conducted on different size of data sets. Each size consists of 20 data sets randomly generated. For small size data sets, CPLEX finds an optimal solution within an average time of 5 seconds which is less than 2 second of an average time needed by the heuristic. For medium-size data sets, the proposed heuristic consumes on average 102 second to produces the same or better solution compared to CPLEX, which represents 1.5% of an average time needed by CPLEX. For large size data sets, CPLEX produces sub-optimal solution for only 5 of 20 data sets in an average time of 3 hours, while the proposed heuristic finds solutions for all large data sets within average time of .8 hours. Therefore, the proposed heuristic outperforms CPLEX for many instances.

A distributed system of Farm dairy company throughout Amman is represented as real case problem. It can be classified as a large size problem; therefore, the cluster method, based on demand and distance, is used first to divide it into four clusters, that helps CPLEX to find an optimal or close to optimal solutions within competitive time, where one vehicle is assigned to serve each cluster. Then a comparison between VRP and EVRP in term of the total traveling time is presented. Based on the result using the electric vehicle instead of the conventional one provides saving in total traveled cost by 9%- 14%.

Finally, the problem can be extended by developing a mathematical model which allows multiple depots and time window for each customer node. Also, different objective functions can be used be side the minimizing the deviation, such as minimizing the number of vehicles used. Finally, in this problem, the recharge stations are assumed available in all time, which is not true in real life. The vehicles spend more time because of the queue in the station, or they may drive to another station. Therefore, the availability of recharging stations can be investigated in future work.