Logistic Optimization for Multi Depots Loading Capacitated Electric Vehicle Routing Problem From Low Carbon Perspective

In this paper, a multi depots capacitated electric vehicle routing problem where client demand is composed of two-dimensional weighted items (2L-MDEVRP) is addressed. This problem calls for the minimization of the transportation distance required for the delivery of the items which are demanded by the clients, carried out by a fleet of electric vehicles in several depots. Since the 2L-MDEVRP is an NP-hard problem, a heuristic algorithm combined variable neighborhood search algorithm (VNS) and space saving heuristic algorithm (SSH) is proposed. The VNS algorithm is used to solve the vehicle routing problem (VRP) sub-problem, and the SSH algorithm is used to solve the bin packing problem (BPP) sub-problem. We propose the space saving heuristic to find the best matching solution between the next loading item and the feasible loading position. The SSH-VNS algorithm is tested by using benchmark instances available from the literature. The results show that the SSH-VNS algorithm has better performance compared with other published results for solving capacity vehicle routing problem (CVRP) and two-dimensional capacity vehicle routing problem (2L-CVRP). Some new best-known solutions of the benchmark problem are also found by SSH-VNS. Moreover, the effectiveness of the proposed algorithm on 2L-MDEVRP is demonstrated through numerical experiments and a practical logistic distribution case. In the last section, the managerial implications and suggestions for future research are also discussed.


I. INTRODUCTION
Vehicle routing problem (VRP) is an important and typical distribution optimization problem in which a fleet of vehicles is required to deliver items demanded by clients at a minimized total cost [1]. Since the complexity in an actual environment, the objectives and constraints encountered are highly variable. Dantzig and Ramser [2] pioneered the study on ''truck dispatching'' problems. Clarke and Wright [3] were the first to incorporate more than one vehicle in the problem formulation. This was considered one of the first studies in VRP literature. Other versions of the VRP emerged in the early 1970s for different types of problems and The associate editor coordinating the review of this manuscript and approving it for publication was Sabah Mohammed . subjects, among all the above variants, researchers have been particularly interested in the CVRP which is the combi-nation of VRP and bin packing problem (BPP). The research of CVRP has far-reaching practical significance, because only by taking both loading and routing into consideration can we make sure the delivery route is the most economic and the items are completely and reasonably loaded into the vehicles. In recent years, some attention has been paid to 2L-CVRP, in which client demands are defined as sets of non-stackable rectangular weighted items, such as furniture, home appliances, or breakables. The 2L-CVRP was first presented by Iori and Vigo [4]. There are mainly four variants of the 2L-CVRP, namely, 2L-sequential oriented loading (2|SO|L), 2L-sequential non-oriented (rotated) loading (2|SR|L), 2L-unrestricted oriented loading (2|UO|L), and 2L-unrestricted non-oriented (rotated) loading (2|UR|L). Leung et al. [5] studied the two-dimensional loading heterogeneous fleet vehicle routing problem (2L-HFVRP), wherein clients were served by a heterogeneous fleet of vehicles and proposed a simulated annealing with heuristic local search (SA-HLS). Zhu et al. [6] addressed a multi depot CVRP wherein the client demands comprised of two-dimensional items (2L-MDCVRP). Dominguez et al. [7] considered a realistic extension of the classical vehicle routing problem wherein both the delivery and pickup demands comprised non-stackable items. They presented a hybrid algorithm that integrated biased-randomized versions of vehicle routing and packing heuristics into a large neighborhood search metaheuristic framework. The 2L-CVRP in an uncertain environment was also studied [8]. With another meta-heuristic for solving 2L-CVRP, Zachariadis et al. [9] proposed the promise routing-memory packing (PRMP).
Recently, more and more researches have been dedicated to the electric vehicle routing problem (EVRP) since concerns have been raised about the related technology and development of electric vehicles for the international automobile market with the global financial crisis, environmental degradation, and energy depletion problems [10]- [13]. Electric vehicles have been used extensively for city distribution logistics. The driving distance of electric vehicles has improved significantly since the development of battery technology, yet not rapidly to meet the demand of large city distribution. Therefore, at current state, changing battery or recharging in public fast-recharging stations are the best feasible solutions to extend the driving distance of electric vehicles. Artmeier et al. [14] studied the most economical one in terms of the energy consumption. This was considered the first attempt to introduce electric vehicles. Later, Conrad and Figliozzi [15] proposed a rechargeable vehicle routing problem with time windows, which became the basis for the studied EVRP with time windows. Schneider et al. [16] studied a vehicle routing problem with intermediate stops considering necessary visits at intermediate locations as the EVRP required.
In recent years, studies have focused on various EVRPs. Sassi et al. [17] addressed a vehicle routing problem with a mixed fleet of conventional and heterogeneous electric vehicles, denoted as VRP-MFHEV, and proposed a multi start iterated Tabu Search (ITS) based on the large neighborhood search (LNS). Desaulniers et al. [18] considered four variants of the EVRP with time windows. They found that allowing multiple as well as partial recharges helped in reducing routing costs and the number of employed vehicles, com-pared to variants with single and with full recharges.
More recently, Sweda et al. [19] studied the problem of finding an optimal adaptive routing and recharging policy for an electric vehicle in a network =. They developed algorithms for finding an optimal a priori routing and recharging policy and then presented solution approaches to an adaptive problem that was built on the priori policy. Subsequently, they presented two heuristic methods for finding adaptive policies, one with adaptive recharging decisions only and another with both adaptive routing and recharging decisions. Barco et al. [20] proposed a model based on the longitudinal dynamics equation of motion and estimated the energy consumption of each Battery Electric Vehicle (BEV) with a case study of an airport shuttle service scenario used to demonstrate the feasibility of the proposed methodology. Montoya et al. [21] extended current EVRP models to consider nonlinear recharging functions. They proposed a hybrid metaheuristic that combined simple components from literature and components specifically designed for this problem. They found that neglecting nonlinear recharging could lead to infeasible or overly expensive solutions.
As discussed above, studies on EVRP only focused on the routing optimization. Few researches are conducted in the combination of EVRP and BPP. Nevertheless, an overall consideration given to loading and routing will contribute to a more comprehensive and reasonable logistics optimization. Besides, there are two main differences between electric vehicle distribution and conventional gasoline vehicle distribution. First, the effect of items' weight on the battery consumption should be considered in the electric vehicle distribution. As a client is served, the items' weight is reduced, and the battery consumption decreased accordingly. Therefore, it is not reasonable to set the battery consumption between two clients as a fixed value. Second, electric vehicles need to go to recharging stations for recharging or for battery replacement during the distribution considering mileage limitations. Thus, it is very important to decide when and where to charge or replace the batteries in the distribution. However, according to our observation, relatively few studies have been conducted on these issues, which is the primary motivation of this research.
The remainder of this paper is organized as follows. The model development in the next section introduces the formulation notations, constraints and complete mathematical model of 2L-MDEVRP. Section 3 introduces a heuristic algorithm combining variable neighborhood search algorithm (VNS) and space saving heuristic algorithm (SSH) for solving the 2L-MDEVRP model. Section 4 introduces three groups of experiments based on benchmark instances including multi deports CVRP, EVRP and 2L-CVRP. Section 5 presents a case study of SH company. In the last section, the key contributions of this research are discussed along with the related managerial implications.

II. MODEL DEVELOPMENT
Considering the characteristics of electric vehicle for distribution, we propose a model for multi depots electric vehicle routing problem with two-dimensional loading constraints. In this problem, the distribution tasks are executed by the same type of electric vehicles. The electric vehicles belong to several depots in different locations. Each electric vehicle is subject to capacity constraint and electricity constraint.
An example of the routing solution and loading solution of 2L-MDEVRP is shown in Fig. 1. Fig. 1(a) illustrates the 2L-MDEVRP that includes two depots D1 and D2, two recharging stations C1 and C2, nine clients and twenty-three items. It can be seen from Fig.1(a) that there are two routes in this example: (D1-1-2-C1-3-D1) and (D2-6-7-8-9-C2-P4-5-D2). After serving clients 1 and 2, the vehicle from depot D1 has to charge in recharging station C1 to continue serving client 3 and subsequently return to depot D1. Similarly, the vehicle from depot D2 has to enter the recharging station C2 for recharging. Fig. 1(b) shows two feasible loading solutions for the above two routes, respectively.

A. FORMULATION NOTATIONS
The following notations are used in the model presented in this paper: The 2L-MDEVRP model considered in this paper is based on some constraints which can be classified into two categories: routing constraints and loading constraints.

1) ROUTING CONSTRAINTS
The routing constraints of 2L-MDEVRP are given as follows: •Closed-tour constraints: one client must be served by a vehicle only once; every vehicle can visit all depots or recharging stations at most once; every vehicle must leave from the specific node where it enters.
•Sub-tour elimination constraint: there is no sub-tour in each tour.
•Demand non-split constraint: all items of one client should be loaded in a single vehicle.
Routing constraints that constitute the feasible domain of the routing solution can be widely found in many VRP related studies [22]. According to the characteristics of the 2L-MDEVRP, we give the mathematical formulas of routing constraints as follows: Equation (1) and (2) denote the conversion formulas of variable a th ii and variable b th i , and the purpose of these two equations is to simplify modeling process by increasing the number of variables.
Equation (3) expresses that a client can be served only once. Equation (4) expresses that every vehicle can visit all depots or recharging stations at most once. Equation (5) expresses that every vehicle must leave from the specific node where it enters. Equation (3)-(5) satisfy the closed-tour constraints. Equation (6) guarantees that no sub-tour will be generated. The solution that satisfying the (3)-(6) constitutes a Hamiltonian cycle. Equation (7) expresses that the weight of items loaded on each vehicle cannot exceed the capacity of the vehicle. Equation (8) expresses that the area of items loaded on each vehicle cannot exceed the area of loading surface of the vehicle. Equation (7) and (8) guarantee that the loaded items cannot exceed the vehicle capacity.

2) LOADING CONSTRAINTS
To ensure that all items can be loaded into the vehicles as required, the loading constraints are introduced as follows: •Rectangle constraint: the loading surface of each vehicle and item can be considered as rectangle.
•Parallel constraint: all items must be loaded with their edges parallel to the edges of the vehicles.
•Directional constraint: all items must be loaded and unloaded straight from the rear door.
•Boundary constraint: all items are not allowed to surpass loading surface area.
•Last-in-first-out (LIFO) constraint: items are not allowed to be rearranged at client sites.
In order to describe the loading constraints in equations, we need to introduce a Cartesian coordinate system which is adopted with its origin in the container's front-left corner, and let (x, y) be the possible coordinates where the front-left corner of an item can be placed. These positions along axes L and W of the container belong to the sets: X = {0, 1, . . . ,L − min(l nk ) and Y = {0, 1, . . . ,W − min(w nk ), respectively. So, the container can be divided into many grids, as shown in Fig. 2. If front-left corner of item k belonging to client n is placed on point C (x nk , y nk ), then the rectangle of ABCD will be occupied.
It's easy to find that the rectangle constraint, parallel constraint and directional constraint are natural satisfaction. Then we need to address other loading constraints as following.

a: BOUNDARY CONSTRAINT
This constraint can be expressed as the coordinate occupied by any item that cannot exceed the loading surface. If item is not rotated, the x-coordinate of front-left corner plus the length of item must be less than carriage length, and the y-coordinate of front-left corner plus the width of item must be less than carriage width. If item is rotated, the x-coordinate of front-left corner plus the width of item must be less than carriage width, and the y-coordinate of front-left corner plus the length of item must be less than carriage length.
There are four possible position relationships between two items as shown in Fig. 3 Considering the simplicity of the mathematical model, we integrate the above four position relationships into one formula. Set item A is the item k belonging to client n, and item B is the item n belonging to client n , n − n + k − k > 0. VOLUME 8, 2020 Firstly, we define ll kk nn and ww kk nn as follows: ll kk nn denotes the value of the length projection on x-axis of item A and item B minus the sum length of the two items. ww kk ii denotes the value of the width projection on y-axis of item A and item B minus the sum width of the two items. If the two items don't overlap, then the maximal value between ll kk ii or ww kk ii is not less than zero. Then we get the constraint as following: And b th n = b th n = 1 express that client n and client n are served by the same electric vehicle.

c: LIFO CONSTRAINT
Usual and practical request in transportation is that, the items belonging to current visiting client can be unloaded from the rear door of vehicle, without having to move items belonging to successive clients along the route. This implies that the portion of loading surface between each item of the client being served and the opening of the vehicle must be empty. The rear door of vehicle is exactly x-axis, as shown in Fig. 3. According to the requirements of LIFO constraint, when vehicle h belonging to depot t travels from client n to client n , we can get the constraint as following: max (ww kk nn , (x nk , x n k )) ≥ 0, k = 1, 2, . . . m n , k = 1, 2, . . . , m n , a th nn = 1, (12) when

III. SOLUTION ALGORITHM
In this section, we introduce the hybrid heuristic algorithm for solving the 2L-MDEVRP model.

A. VARIABLE NEIGHBORHOOD SEARCH
VNS is a metaheuristic for solving the combinatorial and global optimization problems based on the principle of systematic changes of neighborhoods within the search. Many extensions of VNS have been studied for solving large scale instances [23], [24].
In the classical implementation of VNS algorithm, four key components should be specified: (i) method to construct an initial solution; (ii) neighborhood structures NS; (iii) shaking process; (iv) local search. The algorithm framework of VNS is shown by Wei et al. (2015). With an initial solution, a so-called shaking step in which randomly selects a solution from the first neighborhood is performed followed by applying an iterative improvement algorithm. This procedure will repeat until a new incumbent solution is found. Otherwise, one switches to the next larger neighborhood and performs a shaking step followed by the iterative improvement. Once a new incumbent solution is identified, one starts with the first neighborhood. Otherwise one proceeds with the next neighborhood, and so forth.

1) INITIAL SOLUTION
There are many methods that can be used to construct an initial solution, such as minimum spanning tree, random procedure, saving algorithm and so on. However, all these methods are not suitable for dealing with the 2L-MDEVRP proposed in this paper since they cannot select the best time for recharging stations to join the route. To solve this problem, a method of generating initial solution based on scanning algorithm is proposed. The basic idea is to set a rule to scan nodes one by one and classify them into the current circuit as shown in Fig. 4. The concrete implementation procedure of scanning algorithm is as following: (1) The polar coordinate is used to represent depots, recharging stations and clients. We randomly select a depot as the pole of polar coordinate, and draw a line connecting the pole and any client as the polar axis.
(2) The polar axis is rotated clockwise or counterclockwise, and connected to the nodes scanned by the polar axis (neglect the other depots) to generate a route that named current route.
(3) Depending on the vehicle capacity, loading area and battery power, it determines whether the node swept by the polar axis joins the current route. If the swept node is at recharging station, then it is added to the current route directly, and the battery power of electric vehicle is set to full.
(4) If there is no more nodes that can be added to the current route, then we generate a feasible route, and eliminate the clients on the current route from the polar coordinate.
Repeat step (1) to (4) until there is no client node on the polar coordinate.

2) NEIGHBORHOOD STRUCTURES
Six neighborhoods include the 1-1 interchange (swap), two types of the 2-0 shift, the 2-1 interchange, and two types of the perturbation are used in this paper (i.e. k max=6). The order of the neighborhoods is as follows: the 1-1 interchange is set as N 1 , the 2-0 shift of type 1 is set as N 2 , the 2-1 interchange is set as N 3 , the perturbation of type 1 is set as N 4 , the perturbation of type 2 is set as N 5 , and the 2-0 shift of type 2 is set as N 6 . The six neighborhoods are briefly described as follows: The 1-1 interchange (the swap procedure): aims to identify a feasible solution by swapping a pair of clients from two routes. This procedure starts from taking a random client from a random route and tries to swap it systematically with other clients of all other routes. This procedure will not stop until a feasible move is identified. The 2-0 shift (type 1 and type 2): at the beginning of type 1, two consecutive random clients from a random route are identified, and then are checked for possible insertion in other routes. This procedure will repeat until a feasible move is identified. Type 2 is similar to type 1 except that the two clients are considered for insertion into two different routes. These insertion moves are performed in a systematic manner. The 2-1 interchange: This type of insertion attempts to shift two consecutive random clients from a randomly chosen route to another route selected systematically while getting one client from the receiver route until a feasible move is obtained.
Perturbation mechanisms (type 1 and type 2): This scheme was first proposed by Salhi and Rand in VRP by considering three routes simultaneously [25]. The process is to systematically take a client from a route and relocate it into another route without considering capacity and time constraints in the receiver route. A client from this receiver route is then shifted to the third route when both capacity and time constraints for the second and the third routes are not violated. This is the perturbation of type 1. An extension of this perturbation is to shift two consecutive clients from a route instead of removing one client only. We call this as the perturbation of type 2, and these moves are evaluated systematically.

3) SHAKING PROCESS
Shaking is a key process in the VNS algorithm design which aims to extend the current solution search space and reduce the possibility falling into the local optimal solution.
We propose two parameters in the shaking procedure: solution S, and the number of iterations k. In each iteration, a route is chosen randomly, and then it's part of random length is chosen so that a swap-location is not included in the chosen part. Then each client from the selected part is moved to another route, keeping the feasibility and deteriorating the current objective value as low as possible. The obtained solution is identified as the new incumbent solution S, and the whole process will repeat until the maximum number of iterations (i.e., k) is reached.

4) LOCAL SEARCH OPERATOR
In a VNS algorithm, local search procedures can search the neighborhood of a new solution space obtained through shaking in order to achieve a locally optimal solution. Local search is the most time-consuming part in the entire VNS algorithm framework and can decide the final solution quality, so computational efficiency must be considered in the design process of local search algorithm.
In this paper, insertion, exchange and 2-opt are used for local search operator in order to obtain the good quality local optimal solution in a short period.

B. SPACE SAVING HEURISTIC
The optimized routing solution can be obtained by VNS algorithm given in section 4. In this section, we will discuss how to get the feasible loading solution. The quantity and shape loaded on each vehicle are fixed, so the packing problem in 2L-MDEVRP is just a sub problem of traditional packing problem. There are two loading steps of 2L-MDEVRP: determine the items order and determine the feasible loading position list. In addition, we process the space saving heuristic to find the best matching solution between the next loading item and the feasible loading position.

1) LOADING ORDER
For a certain tour, the visiting order of all clients on this tour is fixed. For sequential model described in section 2, other items in the vehicle cannot be moved when items of one client are being unloaded. So the items in a vehicle are sorted in descending order of visit, called OV. Then the order of items belonging to different clients is fixed, and the order of items of same client is not fixed.
For items of one client in sequential model or items of all clients in unrestricted model, we use two orders O1, O2, and O3 to determine their final order. For O1, O2, and O3, items are sorted in descending order of area lw, length l and width w, respectively. O1 is prior to O2 unless the bottom areas of the two items are the same, and O2 is prior to O3 unless the lengths of the two items are the same. The orders of the items loaded for the example of section 2 can be identified according to the above rules:

2) FEASIBLE LOADING POSITION
The feasible loading position list means all positions that can place items, which will be changed after loading an item. Generally speaking, there will be more than one feasible If an item is closer to the rear door than all adjacent items, make a straight line along the internal edge close to the rear door, as Line1 and Line2 shown in Fig. 5(b).

3) SPACE SAVING HEURISTIC
If the loading orders of items are kept the same while the feasible loading position list is constantly updated, it is likely to occur space wasting, which eventually results in the failure to load all the items into the vehicle. As shown in Fig. 6(a), if we load the next item R1 according to the loading order, then the loading surface N will be wasted and no item can be loaded on it. However, if we load item R2 as Fig. 6(b) or load item R3 as Fig. 6(c) before loading item R1, then the waste of loading surface will be avoided. So, we propose a heuristic matching method of loading items and feasible loading position based on space saving.
The matching fitness values (MFV) are used to indicate the matching degree of loading items and feasible loading positions. There are ten situations on the relationship with loading items and feasible loading positions as shown in Fig. 7(a) to Fig. 7(j), where the grey area denotes already placed items and the white area denotes space. There are some independent spaces below the dotted line enclosed by the edges of carriage and the loaded items as shown in Fig. 7. It's better to make full use of these independent spaces, thereby we can avoid the waste of space to a great extent as shown in Fig. 6.

IV. NUMERICAL EXPERIMENTS
This section presents the computational results based on the widely used benchmark instances and some new generated instances. The SSH-VNS algorithm was coded in MATLAB, an all experiments were executed on Intel Core i7-8700K 4.8GHz CPU and 16GB RAM running the Windows10 operating system.

A. RESULTS FOR CVRP AND 2L-CVRP INSTANCES
To test the proposed SSH-VNS algorithm's performance on 2L-CVRP problem we applied it to 180 2L-CVRP benchmark instances introduced by Iori and Vigo [4]. These instances were derived from 36 CVRP in-stances, described by Toth and Vigo [26], where the customer demand is expressed as a set of two-dimensional, weighted and rectangular items. To generate the aforementioned item sets, five classes of the item demand characteristics are introduced [4]. Class 1 defines a single 1 1 item for each client corresponding to basic CVRP instances, so actually there is no difference between class 1 and CVRP instances, and the experimental results of class 1 can be considered as experimental results of CVRP Classes 2-5 contain in-stances with non-unit item sizes. The item size has been randomly generated to define ''vertical'' instances (width greater than length), ''homogenous'' instances (square) and ''horizontal'' instances (length greater than width). These instances can be seen at URL: http://www.or.deis.unibo.it/research.html. Since there is one depot and no recharging station in the 180 instances, we set the number of depots as 1 and set the number of recharging stations as 0.
Average results for classes 1-5 are presented in Table 1. Row ''Num Best'' gives the number of times a method provides the best solution and ''Num Record'' gives the number of times the method gives a solution strictly better than all other methods. As can be observed, the SSH-VNS algorithm refreshes the record of 11% instances (20 out of 180) and achieved an average 1.2% improvement of the best solutions obtained for instances of Classes 1-5. Thus, the proposed SSH-VNS has a good performance for 2L-CVRP, as well as reducing the overall distances between the depot and the clients.
The routing solution provided in this research utilizes electric vehicles to accomplish delivery assignments. Because of the various structure of power generation and different levels of clean energy generation, at the moment, there are substantial variations in the testing and calculating of carbon dioxide emissions of electric vehicles in numerous nations in the whole world. Based on the data from Bloomberg New Energy Finance, the average carbon dioxide emission of gasoline vehicles in the world is 248.5 g / mile. Meantime, the carbon dioxide emissions of electric vehicles in China are 189 g / mile, and the carbon dioxide emissions of electric vehicles in the United States are 147 g / mile. Besides, the carbon dioxide emissions of Japanese electric vehicles are 142 g / mile, German electric vehicles are 140 g / mile, and the British electric vehicles are 76 g / mile. The case analyzed in this paper is in Beijing, China, therefore we picked 189 g / mile as the electric vehicles' coefficient of carbon dioxide emission factor. According to the coefficient, we contrast the energy consumption of complete 2L-CVRP instances adopting electric vehicles and conventional gasoline vehicles.  Hence, we discover that the average carbon dioxide emissions of every instance of the electric vehicle are more than 28%. Among them, 23.9% of the contribution arises from the electric vehicles' lower coefficient of carbon dioxide emission, and 4.1% of the contribution originates from shorter delivery distances.
The parameters of the electric vehicle are set as follows: Q=25, M=1000, η = 0.8, α = 0.08, β = 120. The maximum travelling distance of the electric vehicle without any item is about 200, and the maximum travelling distance of the electric vehicle with full items is about 70.
The optimization results of instance 0703E are shown in Fig. 8. There are five routes in Fig. 8, three of them are processed through depot D1 and two of them are through depot D2. There are two routes through recharging stations: route 1 through recharging station twice and route 4 through recharging station once. It is worth noting that although the distances of route 4 and route 5 are very close, due to the lower weight of loading items, so it is not necessary for route 5 to go through the recharging station. The details of the five routes are shown in Table 3.
For the distribution lines of electric vehicles, which require to charge the extra driving distance in and out of the charging station is increased. Consequently, the driving distance for the electric vehicle is usually higher than the driving distance for the gasoline vehicle. The best solution of instance 0703 was given in Dominguez et al. [7]. The entire distance driven by gasoline vehicles on the additional route of instance 0703 is 709.72. Furthermore, the total distance traveled by the electric vehicle, of instance 0703E extra route is 763.76. If a mile is the distance unit, the carbon emissions between gasoline vehicles and electric vehicles are 176.37kg and 144.35kg. Moreover, the general carbon emissions of electric vehicles and gasoline vehicles are diminished by 18.1%. It means that even if the driving distance needs to increase due to midway charging, the overall view is that the performance of the electric vehicle can still notably diminish carbon emissions.

V. CASE STUDY
SH is a medium-sized logistics enterprise in Beijing, the main business of which is to distribute power distribution cabinets to construction sites in the city. This logistics enterprise has three depots and fifty-two electric vehicles for the distribution. Considering the height of the power distribution cabinets, stacking is not allowed, hence making this a two-dimensional routing problem. When the battery power is not enough to complete the distribution task, it is necessary to go to the recharging station for recharging. At the end of December 2017, there were 1771 public recharging stations in Beijing. Considering the recharging cost, recharging matching, traffic restrictions and convenience, seven public recharging stations are selected as the recharging nodes for all electric vehicles of company SH. In this case, company SH needs to deliver several power distribution cabinets to 17 construction sites. Fig. 9 shows the locations of the depots, recharging stations, and clients. We need to design the shortest distribution route, select the most appropriate recharging stations, and set a loading plan for each vehicle to facilitate the delivery services for clients in need.  Deliveries of company SH in Beijing are handled by a fleet of electric vehicles, which use lithium battery pack as the power source. As shown in Fig. 10, the red circles represent the depots, the black squares represent the construction sites, and the blue triangles represent the recharging stations.
The demands are variety types of power distribution cabinets which can be seen as different sizes of two-dimensional rectangular items since one cannot be stacked on the top of the other one. The vehicles are identical, while they have a weight capacity and a rectangular two-dimensional loading surface. To reduce the real problem to 2L-MDEVRP, we do not consider time windows and pickup of damaged cabinets at the construction sites. The distance between any two nodes (depots, recharging stations and construction sites) is set to the recommended route length of Google Maps, while elevation and slope are not considered in this case.  The data of this instance such as distance, demand and item size are shown in Table 4 -Table 7. The parameters of vehicles: empty weight 4.6t, maximum load weight 2.8t, length 4150mm, width 3300mm, battery capacity 75kWh.  The values of the electric vehicle related to parameters mentioned in (2) are set as: g = 9.8, f = 0.01, C D = 0.3, A = 1.8, ρ = 1.2, η te = 0.85, η m = 0.85, P accessory = 1000. There are five types of power distribution cabinets: XL-21, XL-51, JXF1000, JXF2000 and JXF3000. The size and weight of each type of cabinet are shown in Table 4. The stock of each depot is shown in Table 5, and the demand of each construction site is shown VOLUME 8, 2020 in Table 6. The real distance between any two nodes is shown in Table 7. As the distance unit applied, in this case, is kilometers, the coefficient of carbon dioxide emission demands to convert. Eventually, the converted coefficient of it is 117.44g / km.
By running the optimization algorithm, we get the optimal solution of the practical instance which has four tours as shown in Fig.10. As presented in Table 8, the four tours go through four different recharging stations once in order to complete distribution tasks. The more details of all tours are shown in Table 8.

VI. CONCLUSION
This study investigates the 2L-MDEVRP which is a variant of VRP. The main differences of the 2L-MDEVRP compared to the VRP are: there is more than one depot; the delivery vehicles are pure electric vehicles; there are several recharging stations that can be used for recharging the electric vehicles; the demand of clients consists of weighted, two-dimensional, rectangular items. The aims of 2L-MDEVRP are generating the minimum distance routes, feasibly packing items onto the loading surfaces of the vehicles, and the order of recharging stations on the tours.
The 2L-MDEVRP is of particular theoretical interest as it combines two frequently studied combinatorial optimization problems, namely the VRP and the two-dimensional BPP which are both NP-hard problems. Spontaneously, the 2L-MDEVRP is also a NP-hard problem. To solve this NP-hard optimization problem, this study proposes a performance meta-heuristic algorithmic framework with VNS. On one hand, this study improves a new method for the generation of initial solution to get a better initial solution. On the other hand, this study designs a new heuristic loading algorithm named SSSH which successfully enhances the probability of satisfying the feasible constraints. The numerical experiment reveals that the proposed algorithm can effectively solve the 2L-CVRP and 2L-MDEVRP. Compared with gasoline vehicles, the electric vehicles can remarkably decrease in carbon emissions. Moreover, this study demonstrates the optimization process and optimized solution of a practical distribution case.
Based on the results of this research, the authors plan to construct the model of multi depots electric vehicle routing problem with three-dimensional loading constraints and apply the proposed new method in the future follow-up research. In addition, the authors also consider a further study of a new variant with one additional constraint: time windows.
XIAONING ZHU received the Ph.D. degree in management science and engineering from the University of Science and Technology Beijing, in 2012. She is currently an Associate Professor of technical economy with the University of Science and Technology Beijing. Her current research interests are logistics management and low-carbon management.
RUI YAN received the Ph.D. degree in management science and engineering from the University of Science and Technology Beijing, in 2014. He currently holds a Postdoctoral position with the Beijing Institute of Technology and an Associate Professor of management science and engineering with the University of Science and Technology Beijing. His current research interests are logistics management and low-carbon management.
ZHAOCI HUANG is currently pursuing the degree in mathematics and applied mathematics with the University of Science and Technology Beijing.
WENCHAO WEI received the Ph.D. degree in applied economics from the Catholic University of Leuven, in 2015. He is currently an Associate Professor of logistics management with Beijing Jiaotong University. His current research interests are logistics management and production optimization.
JIAQIN YANG received the Ph.D. degree in business management from Georgia State University. He is currently a Tenured Professor of management school with Georgia State University. His current research interests are operation management and strategic management.
SHAMSIYA KUDRATOVA received the Ph.D. degree in management science and engineering from the University of Science and Technology Beijing, in 2019. She is currently an Associate Professor of business administration with the Jiangxi University of Science and Technology. Her current research interests are investment project management and sustainable finance. VOLUME 8, 2020