1. INTRODUCTION
Cost of transportation has a significant effect on enduser price of the goods, and optimizing the distribution system is a core objective of distribution companies. Today, thousands of companies involved in collecting and distributing goods are in need of an efficient serviceproviding plan in order to reduce their operational costs and satisfy their customers. Therefore, optimizing the distribution system is a major goal of goods and services distribution companies. Evidences show that optimizing service providing and destitution systems plays an essential role in reducing the transportation costs (Toth and Vigo, 2002;Abdollahbeigi, 2021). Regarding the customer satisfaction, it can be said that avoiding losing only 5% of the customers shall result in increasing annual income of service providers and commodity suppliers by 2585% depending on their field of activity. Generally, customer satisfaction in the main factor customer retention. Customer satisfaction shall also result in increased customers’ lifespan (SaorínIborra and Cubillo, 2019;Gaikwad, 2021). Due to promotion of urbanization, many efforts should be made in transportation sector to respond to challenges such as climatic changes, limited access to fossil energy sources, and airquality concerns. European Union (EU) intends to reduce greenhouse gas production by 20% until 2020, and by 40% until 2030, compared to 1990. This is a radical challenge to the transportation sector that produces 20% of the total greenhouse gasses (Schiffer and Walther, 2017;European Commission, 2014;Sokil et al., 2020). Therefore, considering what mentioned above, development of environmentfriendly logistic systems, addressing the significant effects of transportation activities on air pollution and global warming has neem promoted vastly in the recent decades. There has been also a growing social awareness of environmental issues, new regulations, and social incentives that aim at progress of research and development (R&D) activities and investments on transportation and distribution models (Sinha and Labi, 2011). There are plenty of studies on the socalled “Green Logistics”, majority of which aim to promote sustainability of transportation in consideration of environmental and social concerns. Replacing fossil fuel vehicles with the electrical vehicles is an important domain of Green Logistics studies. For instance, United States Environment Protection Agency (USEPA) and National Highway Traffic Safety Administration (NHTSA) are cooperating in development of a new generation clean vehicles, from smallest commercial cars to heaviest trucks, to reduce emission of greenhouse gasses and optimize fuel consumption (European Commission, 2014). These instances show the necessity of studies on environmentfriendly transportation system. Accordingly, the present study attempts to provide a twoobjective model for routing capacitated electric vehicles by a rechargeable time window in order to solve some of the problems in this field. In the proposed model, the vehicles are capacitated and it is difficult to access them in the time window. A hard time window is also considered for providing services to the customers in order to increase customer satisfaction. Additionally, in the second objective function, a deadline is set for the vehicle to arrive at the destination depot; or in other words, a soft time windows is considered to reduce airpollution and achieve customer satisfaction by reducing the time required to reach each customer. Vehicle routing is a key factor in goods distribution management, and its objective function and limitations in goods and services distribution systems vary depending on structural and operational specifications of the distribution system. These differences result in a wide range of circumstances with different operational and structural specifications. Moreover, regardless of applicability of this feature for distribution practitioners, its complex theoretical structure has attracted the attention of researchers in various fields of science. Therefore, this subject has an extensive literature, and tracing its trend of development and evolution is a complicated process. A research on electric vehicle routing problem time window (EVRPTW) and capacitated vehicle routing problem (CVRP) will be explained in this chapter.
Vehicle routing problem (VRP) has been investigated repeatedly by scholars of different fields of science, and many people have attempted to solve this problem. Puchkova et al. (2020) conducted the first study on vehicle routing problem, aiming to solve the Travelling Salesman Problem (TSP) (Puchkova et al., 2020). Clarke and Wright (1964) proposed an innovative method of solving the vehicle routing problem for the first time. VRP was investigated in the early 1970s under titles such as Fleet Routing, Distribution Management, Transportation Network Design, and Routing Public Vehicles. Golden et al. (1977) were the first to use the conventional title “Vehicle Routing Problem” to name their study. After that, other studies aimed at optimizing the distribution systems used the same title. As the studies continued in the early 1980s, the concept of time windows was introduced to the literature of this problem (Solomon, 1984;Ahmadi Daryakenari and Nasiri, 2021;Allakhverdieva, 2020).
Erdoğan and MillerHooks (2012) introduced Green Vehicle Routing Problem (GVRP), where the vehicles can get fuel in the fuel stations proposed along the route. They proposed a mixed integer linear programing (MILP) for minimizing the traveled distance of a noncapacitated green vehicle with limited travel time and a given timespan. Fueling time is fixed and the vehicles are recharged fully. Felipe et al. (2014) developed the green vehicle routing problem using multiple technologies and partial chargers, where charging stations use different charging technologies and vehicle can be recharged partially. They proposed a MILP formulation whose objective function is minimizing the total fueling cost as well as the total traveled distance. They also proposed an innovative local search method as well as a simulated annealing sequence to solve this problem and showed that SA algorithm works more effectively in solving largescale problems. Schneider et al. (2014) presented an electric vehicle routing problem with time window and recharging feature. They solved various small, medium, and largescale problems, including 5 to 175 customer conflicts, using different algorithms and concluded that VNS/TS algorithm has a better performance compared to the other algorithms. So far, the first attempt to solve this problem optimally in a reasonable calculation period in mediumscale instances was made by Koç and Karaoglan (2016). They presented a mathematical formulation that uses a reduced number of variable and limitations and produces a series of valid inequalities. A branch and cut algorithm is used in combination with a simulated annealing sequence that optimizes the first response and upper bounds. Xia and Konak (2016) discusses capacitated vehicle routing in their article. Considering time window for vehicles and the time required to provide service to the customers, capacitated vehicles and traffic conditions along the route were the distinguished parts of their study. They used a combined algorithm consisting of MIP and repeated local searches to solve the proposed problem. Leggieri and Haourai (2017) in their article presented a compressed formulation for the capacitated and noncapacitated vehicles routing problem. In order to optimize the proposed model, they limited their formulation using consequent linearization. Keskin and Çatay (2018) proposed an adaptive large neighborhood search algorithm for solving fastcharging green vehicles routing problem, and took various vehicle recharging methods into consideration. It is included both conventional and electric vehicles in their research and introduce the green vehicle routing problem with mixed vehicle fleet, partial battery recharging, and time window (Macrina et al., 2019).
Zhang et al. (2020) developed a fuzzy mathematical model for vehicle routing problem. In their study, the time window of visiting each customer is taken into consideration and battery charging and replacement stations are defined separately. The vehicle routing problem with random waiting time was developed. In this study, simulation approach was used to confront the uncertainty and an innovative method was developed for solving the mathematical model (Keskin and Çatay, 2018).
The literature review shows that there are a number of gaps in the field of vehicle routing problem (VRP). The present study attempts to fill some of these gaps. Generally, innovations that distinguish this study from the others are as follows:

•Designing a model of capacitated electric vehicles routing problem with hard and soft time windows and battery recharging feature

•Considering a soft time window for the vehicle to reach the destination depot and regarding minimization of this time as an objective function
2. MATERIALS AND METHODS
Today, issues such as reducing costs, improving customers’ satisfaction by ontime delivery of goods/services, and reducing pollutions are the important concerns of transportation industry. The model proposed in this dissertation assumes using electric vehicles to provide services to customers, which shall contribute to reducing pollutions. Different vehicle sizes are used to deliver goods, so that the delivery vehicle size suites quantity of the dispatched goods. Using the same size vehicles will result is using multiple vehicles to deliver large products, or using oversized vehicles to deliver smaller products, which leads to increased costs and environmental risks. In addition, in order to reduce the costs, one of the objective functions consists of reducing travel costs, fuel capacity, and fuel consumption. In order to improve customers’ satisfaction by ontime delivery of goods a limited hard time window is assumed for the vehicles to start the service. In the second objective function, the time to reach the last depot has been minimized.
There are a number of assumptions to explain the mathematical model, where Points “C” represent “Customer” and Points “R” show the “Recharging Stations”. This routing problem includes only one depot. The proposed routing problems consists of G = (V, A) Graph, where $V=\left\{0\right\}{{\displaystyle \cup}}^{\text{}}C{{\displaystyle \cup}}^{\text{}}R$ , which indeed shows the existing knots consisting of depot, customers, and researching stations; and A shows the series of curves. Capacitated vehicles with different capacities start to travel from the depot. Each vehicle has a level of energy consumption that is determined by the distance it travels.

1. The vehicle tour starts from the departure depot and the vehicle is charged at the beginning.

2. The tour ends in the depot and the arrival time has time window that is soft and can exceed the allowance; but this value is included in the second objective function.

3. Each vehicle is selected for only one route.

4. The time to reach the customer to provide services has a hard window and cannot exceed a given timespan.

5. Availability of the vehicles has a given timespan.

6. There are several vehicle types with different capacities.
Indexes

i and j : Knot Counters

K : Vehicle Counter

F' : Series of Recharging Stations

F'_{0} : Series of Recharging Stations plus Zero Depot

V : Series of Customers V= {1, 2,…, N}

V_{0} : Series of Customers plus Zero Depot

V' : Series of Customers plus Recharging Stations $\text{V}\text{'}=\text{V}{{\displaystyle \cup}}^{\text{}}\text{F}\text{'}$

V'_{0} : Series of Customers plus Recharging Stations plus Zero Depot $\text{V}{\text{'}}_{0}={\text{V}}^{\text{'}}{{\displaystyle \cup}}^{\text{}}\left\{0\right\}$

V'_{N+1} : Series of Customers plus Recharging Stations plus N +1 Depot $\text{V}{\text{'}}_{0}={\text{V}}^{\text{'}}{{\displaystyle \cup}}^{\text{}}\{\text{N}+1\}$

V’_{0N+1} : Series of Customers plus Recharging Stations Plus Zero and N +1 Depot ${\text{V}}^{\text{'}}{}_{0,N+1}={\text{V}}^{\text{'}}{{\displaystyle \cup}}^{\text{}}\{0,\text{N}+1\}$
Parameters

T_{max} : Maximum Time for the Vehicles to Reach the Depot

h_{k} : Energy Consumption Rate of Vehicke k

d_{ij} : The Distance between Knots i and j

t_{ijk} : The Time Needed to Travel the Distance Between Knots i and j by Vehicle K

S_{i} : Time of Providing Services to the i^{th} Customer

E_{i} : Earliest Time of Providing Services to the i^{th} Customer

l_{i} : Latest Time of Providing Services to the i^{th} Customer

q_{i} : i^{th} Customer’s Query

g : Recharging Time Rate

Q_{k} : Fuel Capacity of Vehicle k

C_{k} : Capacity of Vehicle k

B_{k} : Beginning of Availability of Vehicle k

E_{k} : End of Availability of Vehicle k

M_{k} : Overloading Penalty PerVehicle PerKnot

O_{k} : Remaining Charge Penalty PerVehicle Per Knot
Variables

τ_{ik} : The time for Vehicle k to Reach Knot i

u_{ik} : Remaining Load of Vehicle k When Arriving Knot i

y_{ik} : Remaining Charge of Vehicle k When Arriving Knot i

Cijk : Equals 1 if the Distance between Knots i and j is Travelled by Vehicle k; Otherwise Equals 0
Phrase (1) is the first objective function, which is reducing the electric energy consumption by reducing the travelled distance; because consuming electric energy requires use of fossil fuels for power generation. Therefore, this model attempts to minimize the travel distance to reduce the fuel consumption. The costs are also minimized by minimizing the electric power consumption and remaining load of the vehicles at the meeting knots, as the objective is keeping to extra loads, which result is increased energy consumption, in the vehicles. Assuming that battery level and load of the vehicle are intended to be zero in the last depot, cost of these two items is taken into consideration in the objective function. Phrase (2) is the second objective function that is minimizing delayed arrival of vehicles in (N+1)^{th} depot. This objective function is rarely taken into consideration in other studies, and arrival time of vehicles to destination depot has not been used as an independent objective function.
Phrase (3) shows the meeting knots to assure that each costumer is met exactly once by a vehicle to respond the customers’ queries. Phrase (4) explains that the vehicle has maximum one output at the recharging station. Phrase (5) shows that input and output of a knot should be the same; in other words, a vehicle that enters a knot must exit that know. Phrase (6) shows number of routes or number of vehicles; i.e. there is a route per each employed vehicle. Phrase (7) shows the arrival time in consecutive knots, when we are in the depot or customer knot, the time of arrival to the next knot must be greater that of the previous knot. Phrase (8) show the time of arrival to consecutive knots when we are in the recharging station knot. Phrase (9) shows the arrival time in customer knot with hard time window; because each customer needs the selected good or product in a given timespan. Phrase (10) shows the vehicle load in consecutive knots; i.e. when the load (less than or equal to the vehicle capacity) is embarked in the origin depot, the load is reduced after reaching each customer and disembarking the customer’s order. Phrase (11) shows the load at the origin depot, which equals maximum to the vehicle capacity. Phrase (12) shows the vehicle battery level in consecutive knots, when the origin is the customer’s location. Phrase (13) shows the vehicle battery level in consecutive knots, when the origin is the recharging station or the 0^{th} depot. This shows that when the vehicle departs the recharging station or 0^{th} depot, the battery has the maximum level of energy. Phrase (14) shows availability of the vehicle with time window, and this time window is a hard time window.
3. RESULTS AND DISCUSSION
The proposed model is solved using Epsilon Constraint Method (EC) and NSGAII Metaheuristic Algorithm, which will be explained below:
Epsilon Constraint Method (EC)
General form of a MODM problem is as follows:
Assume that the first objective is selected as the main objective and the other objectives are applied to the upper bounds of the epsilon constraints and problem constraints. In this case, EC method is used to work out the following oneobjective model:
where the first objective is regarded as the main objective, and second to n^{th} objectives are limited to maximum value of e_{i}. It seems necessary to explain that if either of objectives (for example f_{k}(x)) is maximization, the general form is still applicable if ${{f}^{\prime}}_{i}\left(x\right)={f}_{i}\left(x\right)$ is defined and regarded as the k^{th} objective. Another solution is to define the constraint related to this objective as:
A decimal string is used to present the answer, which is regarded as the answer presentation code, and implies one answer to the problem when decoded as follows:
Length of this string is the answer to number of cities plus recharging stations plus first vehicle. Putting this string into an ascending order shall result in a permutation.
Larger sums of cities and recharging points are selected as dividers. In the permutation, the numbers between two dividers show the route of a vehicle. In the example shown in the tables, the first vehicle is not used; cities 2, 3, 5, and 7 are in the route of the second vehicle, and cities 1, 4, and 6 are in the route of the third vehicle. Other variables are determined considering the travel route, and regulations related to charging/discharging the vehicle, loading capacity, and other variables are calculated in a cycle; i.e. they are determined at beginning of the route and the new battery levels and loads are calculated after leaving each city. The next part provides an overall explanation of NSGAII Metaheuristic Method and its application in solving the twoobjective optimization problem of this study.
The major steps of running NSGAII to solve a MODM problem is as follows (Jie et al., 2019):

•Developing initial population

•Calculating fitness criteria

•Nondominated sorting of population and calculating crowing distance

•Recombination and mutation to produce new children

•Recombination and mutation operator
In this operator, first two parents are selected for reproductions. The parents are named X1 and X2 respectively. Each parent is a string with a length of N. Then a string with length of N is produced randomly by numbers between zero and one, named α.
The two parents’ corporate two children named Y1 and Y2, each one of whom has a share of each parent depending on α; i.e. each gene of each child is made of the parents’ genes in the following formula. Therefore, two children are inherited from the two parents.
Mutation operator is applied to a part of the population, depending on the mutation rate, and each mutation adds a new member to the mutated population. These could be fitter or worse than the existing members could. Mutation gives the algorithm the chance of moving from convergence to local optimization. In this operator, after selecting a member for mutation, a number of the member’s genes are selected considering the mutation step Mu and change randomly.

•Combination of initial population and the population formed by recombination and mutation

•Replacing the parent population with the best members of the combined population from the earlier stages

•All stages will be repeated until the desired generation (or optimum conditions) are achieved.
In the first stage, the lower ranked members replace the earlier parents, and then will be sorted based on crowding distance. The initial population and population formed by recombination and mutation are first classified based on their rank, and those ranked lower are removed from the population. In the next stage, the remaining population are sorted based on crowing distance. Here, sorting is done within the same field. The process is shown in Figures 1 and 2.
Different criteria can be taken into consideration for termination of the algorithm, and usually one criterion is used for termination. One criteria can be the conditions where the algorithm cannot optimize the best answer after a given number of repetitions.
Another criterion can be the conditions where average fitness of existing answers in the current population is the same as or very close to the best answer. Alternatively, we can decide from the beginning how many times the algorithm can run (this is the most popular criterion). In the proposed NSGAII method, the termination condition is running the algorithm for a given number of times. This number is regarded a NSGAII parameter like mutation rate, recombination rate, population, etc.
In this part, a smallscale problem is designed to validate the model. The problem is solved in Gams commercial software using epsilon constraint method and NSGAII algorithm. Pareto borders formed by the two methods and the routes travelled by the vehicle for each Pareto point are shown. This problem consists of five customer knots, 2 charging station knots, and 3 capacitated vehicles. Information of the problem are shown in Table 4. Before evaluating performance of the proposed algorithm, their parameters are adjusted using Taguchi method. Results of adjusting the proposed algorithm parameters are shown in the following tables.
The Pareto front produced using EC method and NSGAII algorithm are shown in Figure 3. Routes travelled by the vehicles are also shown in Figure 3.
In this part, 20 problems in different scales with increased number of knots and vehicles are solve using NSGAII algorithm and EC method in order to evaluate performance of the proposed algorithm. All results are extracted from running the program in a notebook system with Intel ® Core^{™} i5 processor, 4 Gb RAM, and Windows 10 Ultimate Microsoft operating system. Considering the fact that timing and routing models are NPHard, increasing scape of the problem results in increased calculation time. In addition, as scale of the problem is increased, the software cannot find the optimum answer and GAMS answers will be no more optimal. Version 2014a of MATLAB software is used to present the proposed algorithm. Now, the question is how can we use the algorithm outputs shown is the said tables to evaluate the algorithm and compare it to the precise method. If the problem is a oneobjective optimization problem (optimizing z_{1}), it is clear that the best solution is the one that offers a solvable value for which z1 is smaller. However, in MODM problem, the evaluation method is different and the problem cannot be evaluated as a oneobjective problem. Generally, in evaluation of MODM solutions, it is important to note that any solution that offers maximum spread or diversity (MD), mean ideal distance (MID), mean RAS and spacing (S) values, and mean algorithm running time (T) (which is used to compare the calculation pace) has a better performance. Accordingly, the following criteria are defined to evaluate two MODM methods and comparing them (Madera et al., 2017).
Algorithms have different running times, and the researcher uses metaheuristic algorithm to work out the answer in a reasonable time. Therefore, any algorithm that finds the answer in a shorter time span is better for the researcher.
In these criteria, first, an ideal answer is assumed for the problem, and then the average deviation of Pareto answers from the ideal answer is calculated (Deb et al., 2002). The ideal answer that is shown by I_{sol} is the state where both answers are optimal at the same times (I_{sol}=(min(z_{1}), min(z_{2})). If origin of the coordinates is regarded as the ideal answer, I_{sol}=(0,0). Assuming that F(A) stands for praetor front of solution A, then the mean ideal distance criterion is calculated using Equation (19):
where $\left\rightpa{I}_{sol}{}_{2}$ is the Euclidean distance and $pa\u03f5F(A)$ shows the distance between the algorithm result and the ideal answer? It is clear that the smaller the MID criterion is the better the algorithm would be (Golden et al., 1977).
Maximum spread/diversity criterion measures diameter of the spatial cube used by the final objective values for series of nondominated answers. The criterion is calculated using Equation (20). The larger this value is the better the solution would be (Deb et al., 2002).
Another evaluation criterion proposed in the present study is RAS. RAS equation is shown below (Deb et al., 2002)
where ${F}_{i}=\text{min}\{{f}_{1i},{f}_{2i}\}.$
Smaller RAS values are preferred. All indexes consider positive objective values; therefore, the negative objectives must be converted to the positive zone of the axis (Deb et al., 2002).
After defining the evaluation indexes to evaluate performance of the algorithm, 20 problems were designed as shown in Table 5 that contains the problems specifications. The algorithm was run 10 times for each problem, and the average value was reported as the final answer of each solution. Values of all indexes for NSGAII algorithm are shown in Table 5. It seems necessary to mention that GAMS software has been able to solve only the first four problems. Diagrams of the four metrics selected for the fi rst four problems is shown in Figure 4.
4. CONCLUSION
Capacitated electric vehicles routing problem is an important and frequently used problem in supply chain and logistics management. In the present study, time window constraints are used to bring the problem closer to the real world and make it more functional. Hard time windows were used for providing services to the customers and soft time windows were used for arrival of vehicles in the origin depot, and this time is regarded as the second objective function. The proposed model is solved using the precise epsilon constraint (EC) method and considering the fact that the mode is NPHard in medium and large scales, it is impossible to find optimal answers that are positioned on the overall front. Therefore, metaheuristic methods were used to solve this problem. NSGAII is the metaheuristic method used in this study. Afterward, answers and algorithm operators used in the study were explained, and parameters were adjusted using Taguchi method. After comparing results of the algorithm with those of the precise method, efficiency of the proposed method in solving largescale problems was validated and finally sensitivity of the model was analyzed.