VRP in urban areas to optimize costs while mitigating environmental impact

Nowadays, the need to think about sustainable mobility, both goods and people, is widely recognized. For this reason, many recent papers have moved in this direction. In this context, particular attention is now devoted to urban mobility, mainly from a smart city perspective. The present work focuses on sustainable urban freight distribution and proposes a variant of the VRP, which presents some innovative aspects. The goal is to minimize the routes’ cost components, including traveling and external costs due to environmental issues, depending on the chosen vehicles and the different urban streets to cross. In addition, restrictions on the maximum duration of each route to ensure frequent sanitation of vehicles used for deliveries, as required from the beginning of the COVID-19 pandemic, are imposed. The distribution network is modeled by a weighted digraph for which some properties are proved. To face the problem, we present a mixed-integer linear programming model, a math-heuristic associated with it, and a memetic algorithm approach. The results of the reported computational experimentation with random instances specifically tailored for the problem show the efficiency of the proposed methods. Further, test cases based on data of the distribution network of two B2C companies operating in the city of Genoa, Italy, proved the effective application of the proposed methods in the direction of sustainable urban distribution plans.


Introduction
From the last few years, almost daily, media reports of damage caused by weather and atmospheric phenomena. Public opinion agrees on the need to change their lifestyles towards more environmentally sustainable behavior, while both scientists and politicians call for a reduction of pollutants to protect the environment. Above all, logistics and transport enterprises cannot ignore the need to establish shipping and routing plans aimed at reducing, among others, emissions, noise, and accidents along their paths. In this direction, recent literature has presented many works, having as one of the main goals the minimization of the mentioned above external Carmine Cerrone and Anna Sciomachen have contributed equally to this work. costs and their evaluation. More precisely, some papers in the last few years proposed either different ways to calculate and evaluate externalities (see, e.g., Bektaş et al. 2019;Bigazzi and Figliozzi 2013;Dekker et al. 2012;Petro and Konečnỳ 2019) or pricing policies to reduce them (see, among others, Kellner et al. 2017;Chang et al. 2018). Other research works presented strategies to incentivize vehicles or transport modalities with lower environmental impact (Okushima 2018;Xu et al. 2019). In particular, taking into account also the social need and the impact that distribution of goods has on citizens, in the current decade, many studies have been proposed aimed at sustainable urban mobility, especially from a smart city perspective (Carrabs et al. 2017;Kauf 2019;Szymczyk and Kadłubek 2019;Strale 2019).
It is true that smart mobility is a crucial topic able to impact the well-being of citizens and the environmental policies of the cities (Taniguchi et al. 2012). As part of the policies implemented in the context of smart mobility to reduce the environmental impact, the distribution of goods, and more precisely the home delivery due to e-commerce, has focused a lot of attention. The goal is twofold: to analyze the likely impact of distribution operations in urban areas upon pollu-tion and to see how sustainable freight distribution can be assumed (Figliozzi 2010;Cerulli et al. 2018;Heshmati et al. 2019).
From an operations research and management science point of view, the Vehicle Routing Problem (VRP) is widely used in decision-making processes about the distribution of goods and services. Today, models and solution methods of VRP are encountered very frequently in the optimization of distribution plans where sustainable factors have to be taken into proper account (Bektaş and Laporte 2011;Behnke and Kirschstein 2017;Demir et al. 2014;Ehmke et al. 2018;Erdogan and Miller-Hooks 2012;Koyuncu and Yavuz 2019;Kramer et al. 2015;Yu et al. 2019). In the very recent literature, VRPs have been proposed in which, in addition to the need to plan the distribution of goods in urban areas in a sustainable manner, they also consider the need to contain the spread of the COVID-19 virus (Chen et al. 2020;Cerrone et al. 2021). In accordance with this topic, the present work focuses on sustainable urban freight distribution within a smart city framework and proposes a variant of the Green VRP (GVRP), in which mobility models are suggested by the needs arising from the COVID-19 pandemic are also taken into accounts. In particular, in the proposed GVRP, we have to choose from an available fleet of non-homogeneous vehicles, those vehicles to deliver goods of a large-scale retailer company in an urban area. Vehicles differ for their size, capacity, and pollution class. In addition, government traffic restrictions and street limitations are given. The goal is to minimize the cost components, which include traveling, loading, and external costs, due to environmental issues, depending on the type of chosen vehicles and urban streets to cross. The present work originates from and extends the results presented in Carrabs et al. (2017), where the authors proposed a Mixed Integer Linear Programming (MILP) model for defining delivery plans in a Business to Consumer (B2C) context. More precisely, we present a refinement of the model given in Carrabs et al. (2017), by including valid inequalities and examining additional constraints, depending on the vehicle type and EURO-Class. The minimization of emissions and pollution reduction is the main objective or part of the generalized cost function in many VRPs. In literature, these problems have often been called Pollution Routing Problem (PRP): Erdogan andMiller-Hooks (2012), Xiao et al. (2012). In this paper, we propose two different approaches derived from the recent literature to solve our GVRP separating two of the most crucial aspect of any optimization technique: the improving phase and the feasibility restoring.
The organization of the paper is as follows. In Sect. 2 we introduce in detail the underlying optimization problem together with the required notation. Some properties concerning the derived network model are also given. In Sect. 3, we give the formulation of the related MILP model. In Sect. 4 the proposed heuristics are given. Section 5 presents the com-putational experimentations, based on both random instances and real size instances derived from a retail company located in Genoa, Italy. Finally, we give some conclusions and outline future work.

Problem definition
As already mentioned, in the present work, we are involved with the definition of sustainable distribution plans of goods in urban areas. This problem can be considered as a variant of the so-called GVRP. More precisely, given a number of customers located in a metropolitan area requiring a known amount of goods and a fleet of not homogeneous vehicles, we have to determine i) a set of routes originating and terminating at a depot and servicing all the customers and ii) which vehicles to select for traveling along the routes. Vehicles differ for their size, capacity, cost, and pollution class. The goal is to minimize the total cost given by traveling, loading, and pollutant cost components. Note that this goal aims at favoring the use of the most eco-friendly vehicles, especially in the city center and close to schools and parks (Carrabs et al. 2014).
To model the problem we are focused on, a digraph G = (N , A) defines the underlying urban transportation network, where N = {d} ∪ U is the set of n + 1 nodes, representing the set U of customers and the depot node d; A is the set of m paths of the city connecting pair of nodes.
For each customer i ∈ U , we indicate with q i the required amount of units of goods. A value l i j is associated with each arc (i, j) ∈ A, representing the length of the corresponding road connections.
Set V contains the fleet of vehicles. To each vehicle, v ∈ V has associated a quintuple of parameters (w, e, k, c, s) specifying, respectively, its class of weight, pollution emission type, capacity, cost/km, and setup cost. All vehicles are grouped into classes t, such that given two vehicles Note that the class of weight of the vehicles determines the subset of urban paths that are inhibited to them, while the type refers to their category in terms of polluting emissions. Finally, note that the cost of each vehicle depends on both its class of weight and the pollution emission type, which refers to the standard EURO-classes from III to VI (see Sect. 2.1).
The solution to the problem is the set of minimum-cost vehicle routes for delivering the required goods to all customers. As in the classical VRP, some conditions have to be satisfied: (i) all nodes, except node d, must be served only once by a vehicle; (ii) all vehicles return to node d after completing their service; (iii) all vehicles have a maximum capacity k v .
It should be noted that in this work, the time windows for the delivery of goods are given implicitly. In fact, the delivery plans are defined considering only the time slots of a day in which in the city of reference is allowed the transit of commercial vehicles. In particular, the time slots available for the delivery of goods vary from 2 to 4 h depending on the parameters e and w associated with each vehicle v ∈ V . Thus, let L v be the maximum length of a tour traveled by vehicle v, ∀v ∈ V , depending on the dimension of the time slot and the parameters w v and e v . In order to obtain a feasible solution, it is important to take into account the service time of each customer. Considering the average speed of the vehicles, in this paper, the service time is transformed into a distance value L c , so that it can be included in the constraint on the maximum allowed length L v of each tour. Furthermore, during this period of the COVID-19 pandemic, the presence of these constraints allows the sanitation of the vehicle (and drivers) several times during the day.

Cost of arcs
In this work, depending on the pollution emission type e of each vehicle v, the air pollution (a v ) is quantified economically, and its impact is added to the cost of each arc in the graph, also considering the marginal climate change costs (m v ). We calculate the costs of the arcs, according to the results shown in Cerulli et al. (2018). In particular, the air pollution and marginal climate change costs are considered with respect to three different cases, corresponding to EURO classes III, IV, and V-VI of the vehicles, respectively. The costs due to air pollution and marginal climate change of the fleet of vehicles expressed in Euro/km are reported in Table 1, and refer to the two main types of commercial vehicles, namely light commercial vehicles (LCV) and heavy good vehicle (HGV). The costs reported in Table 1 have been derived from (Ricardo 2014).
In this scenario the cost of crossing arc (i, j) ∈ A by vehicle v is calculated as: The parameter δ i j ≥ 1 is associated with each arc (i, j) ∈ A and is used to increase the traversing cost of a street that passes too close to schools, parks, or other places of interest.
In many real-world scenarios, crossing a road by vehicles with a specific class of weight (w) could be prohibited (see Fig. 1). Considering that a cost c i jv is associated with each vehicle v ∈ V , in case the road is prohibited for v, logically we set c i jv = +∞ (see Fig. 1). Obviously in all the procedures implemented in this work, arcs with cost c i jv = +∞ are not created.

Graph reduction
Graph G = (N , A) represents the road network within which our vehicles move. In Sect. 2 we have specified that the set N of nodes is equal to {d} ∪ U . This implies that excluding the depot d, all the nodes of the graph represent customers, and there are no transit nodes.
To transform a generic street networkĜ = (N ,Â) in graph G, we use the following procedure (see Algorithm 1 -Graph Reduction).
In lines 4 and 5 of the algorithm, according to the customers and the depot ofĜ, the corresponding nodes are created in G. At line 6, for each pair of nodes i, j ∈ N the arcs (i, j) and ( j, i) are created in A. On lines 8 and 9 the shortest paths between nodesî,ĵ ∈N related to nodes i, j ∈ N are computed. On lines 10 and 12, the cost of each arc (i, j) ∈ A is created in relation to the specific vehicle v. Finally, lines 11 and 13 show the procedure for calculating the route length of the paths between each pair of nodes in G on the basis of the average speed of the vehicles and the service time at the nodes, as described above. Figure 2 reports an example of how the graph reduction algorithm works. The nodes of the street networkĜ ( Fig. 2a) are the depot d, the customersn 2 andn 3 , and the transit nodeŝ n 1 andn 4 . The dashed arc (n 1 ,n 3 ) represents a street prohibited to the transit of HGV vehicles. The weights on the arcs of Fig. 1 Example of calculation of the arcs' cost for two vehicles v 1 and v 2 . Dashed lines represent streets that cannot be traveled by HGV vehicles. Therefore, arc (1, 2) is not compatible with the weight class w v2 . δ 2d is set to 2 since arc (2, d) is adjacent to a school Fig. 2 An example of the Graph Reduction Algorithm. In G the cost c dn2v1 of crossing arc (d, n 2 ) ∈ A with vehicle v 1 , represents inĜ the cost of the shortest path (d,n 1 )(n 1 ,n 2 ) using the same vehicle v 1 .
The cost c dn2v2 of crossing arc (d, n 2 ) ∈ A with vehicle v 2 , represents inĜ the cost of the shortest path (d,n 1 )(n 1 ,n 3 )(n 3 ,n 4 )(n 4 ,n 2 ) using vehicle v 2 G are those derived in Fig. 1. Figure 2b depicts the resulting graph G, which is our referring network model. Arcs of G are the solution of the shortest path problem between customers and customers and the depot, according to lines 8 and 9 of the Graph Reduction algorithm. Note that the cost of the path from the depot to customer n 2 traveled by vehicle v 2 is much higher than the same path traveled by vehicle v 1 , due to the traffic limitation of graphĜ.

Graph reduction property
The graph reduction function allows us to process smaller inputs speeding up the solution approaches. This section will analyze the impact of the transformation on the optimal solution produced by an exact approach.
Proposition 1 LetĜ be the graph of the road network and G = Graph Reduction(Ĝ) the reduced graph obtained by applying the algorithm (Graph Reduction). The value of the optimal solution calculated on graph G is greater than or equal to the value of the optimal solution calculated on grapĥ G. Fig. 3 a Street network graph; nodesn 1 andn 2 are customers. Maximum length imposed for each route equal to 6.Ŝ = (d,n 1 ,n 3 ,n 2 ) optimal solution inĜ, cost = 4, length = 6. b Graph Reduction Algorithm example from graphĜ to graph G. S = (d, n 2 , n 1 ) optimal solution in G, cost = 6, length = 6. c Graph Reduction L1 Algorithm example from graphĜ to graph G . S = (d, n 1 , n 2 ) optimal solution in G, cost = 3, length = 6 Proof Starting from a solution S in G it is always possible to create a solutionŜ inĜ with the same objective function value. The assumption is trivial since each route R ∈ S can be transformed into a routeR ∈Ŝ by replacing the arcs of R with the corresponding shortest paths inĜ, thus keeping the cost of the route, its travel time and its total demand unchanged. This assures that it is not possible to create a solution in G better than the optimal solution presented inĜ. To prove this, we will use the example reported in Fig. 3b.

Proposition 2 LetĜ be the graph of the road network and G = Graph Reduction(Ĝ) the reduced graph obtained by applying the algorithm (Graph Reduction). If the cost and the length of each arc are such that:
then the cost of the optimal solution calculated on graph G is equal to the cost of the optimal solution calculated on grapĥ G.
Proof GivenŜ a feasible solution of our problem on grapĥ G, letR be a sequence of nodes representing a generic route belonging toŜ associated with a generic vehicle v ∈ V . Removing fromR all nodes with null demand, we will get a new sequence of nodes R. Considering that by construction, between each pair of successive nodes in R there is a path inĜ, so R is a route also in G. Given two consecutive nodes (n i , n j ) in R, inR there is the path composed of nodes (n i , . . . , n j ). Further, since each arc of G represents the minimum cost path between the same pair of nodes in G, then the cost c i jv of the arc (n i , n j ) traveled with vehicle v ∈ V is less than or equal to the cost of the route from node n i to node n j inĜ with the same vehicle. This means that the cost of the route R will always be less than or equal to the cost of the routeR. Based on inequality (2) we can say that the travel time of the arc (n i , n j ) is less than or equal to the travel time of the path from node n i to node n j inĜ. This means that, according to the given travel time constraint, the length (travel time) of route R will be less than or equal to that of routeR. Since routesR and R have the same customer demand, the maximum capacity of the vehicle v ∈ V is never exceeded. With this simple conversion we have shown that it is possible to transform each solutionŜ in G into a solution S in G having a value less or equal to the objective function. Now using proposition (2.3) our proof is complete.
for all v ∈ V do 5:Ŝ P i, j,v = Shortest path fromî toĵ inĜ with respect to the length of the arcs. 6:Ŝ P j,i,v = Shortest path fromĵ toî inĜ with respect to the length of the arcs. 7: l jiv = (n1,n2)∈Ŝ P j,i ln 1n2 9: end for 10: end for 11: Return G.
Starting from the Graph Reduction algorithm, we created the Graph ReductionL1 algorithm by changing the calculation of the length of the arcs in G. In particular ∀(i, j) ∈ A, ∀v ∈ V , l ilv = minimum path between i and j inĜ with respect to the length of the arcs inÂ.
Proposition 3 LetĜ be the street network graph and G = Graph ReductionL1(Ĝ) the reduced graph obtained by applying the algorithm (Graph Reduction L1). Then the cost of the optimal solution calculated on graph G is less or equal to the cost of the optimal solution calculated on graphĜ.
Proof GivenŜ a feasible solution of our problem on grapĥ G, letR be a sequence of nodes representing a generic route belonging toŜ associate with a generic vehicle v ∈ V . Removing fromR all the nodes with null demand, we will get a new sequence of nodes in R. Considering that, by construction, between each pair of successive nodes in R there is a path inĜ, R is a route also in G. Given two consecutive nodes (n i , n j ) in R, inR there is the path composed of nodes (n i , . . . , n j ). Since in G the cost of each arc is the minimum cost path between the same pair of nodes inĜ, the cost c i jv of the arc (n i , n j ) traveled with vehicle v ∈ V is less than or equal to the cost of the route from node n i to node n j inĜ, using vehicle v ∈ V ; this means that the cost of route R will always be less than or equal to the cost of routeR.
Moreover, the length l i jv of the arc (n i , n j ) traveled with vehicle v ∈ V is less than or equal to the cost of the route from node n i to node n j inĜ; this means that the length of route R will always be less than or equal to the length of routê R, respecting the maximum duration constraint. Considering also that between routesR and R the sum of the customer demand does not change, then the maximum capacity of vehicle v ∈ V is never exceeded. With this simple transformation we have shown that it is possible to transform each solution S inĜ into a solution S in G with objective function value less than or equal to that ofŜ. To prove that it is possible to create a better solution in G than the optimal solution inĜ, we will use the example in Fig. 3c.

Mixed integer linear programming model
In order to model the distribution problem in urban areas defined in the previous session, we define the following decisional variables.
, representing the amount of goods carried by vehicle vfrom ito j.
Let us now give the proposed MILP model of the problem that is derived from Cerulli et al. (2018). Note that the given MILP model refers to a decision planning problem of a single time slot within a working day.
v∈V i∈N In the above model, the objective function (3) minimizes the total distribution cost to serve all customers, given by traveling, loading, and environmental cost components, as defined in Sect. 2.1, plus the setup cost associated with the activation of each vehicle. Constraints (4) impose that a customer is served exactly once by a vehicle, while constraints (5) impose that if a vehicle delivers goods to a customer, it must depart from it.
Constraints (6) are the so-called flow-commodity constraints, specifying the difference between the quantity of goods a vehicle carries before and after visiting a node. Constraints (7) guarantee that the capacity of the selected vehicle is never exceeded. Constraints (8) impose that the value y i jv must be at least equal to demand q j of node i if a vehicle v crosses the arc (i, j) ∈ A. Constraints (9) ensure that the maximum length of the route traveled by each vehicle v is no longer than a given value L v . Finally, (10) and (11) define the decision variables.
To reduce the symmetry of the above model (3)-(11), for each class of vehicles t we added the following family of constraints (12) to the model. Suppose that in a given class there are h vehicles {v 1 , v 2 , . . . , v h }. For this class of vehicles we added (h − 1) anti-symmetry constraints. The aim is to avoid the generation of identical solutions which differ in the allocation of routes to different vehicles of the same class. In particular, to create these constraints, we associate a unique integer number id (d, j) to each outgoing arc from the depot.
Constraint (12) allows us to create a route/vehicle assignment order. The ordering ensures that in the case of h equivalent vehicles, the route vehicle assignment will always be unique. This is true because thanks to our reduction algorithm (see Sect. 2.2), the arc leaving the depot is associated with the first customer served on the route, making the integer id(d, j) unique for each route.

Heuristic approach
In our study, we present two different approaches to solve our GVRP problem: a math-heuristic approach and a memetic approach.

Math-heuristic approach
The math-heuristic methods proposed in the literature to solve a different variant of the VRP could be classified in the following three classes (see the survey by Archetti and Speranza 2014).
• Decomposition approaches: the problem is split into smaller and simpler subproblems, and a specific solution method is applied to each subproblem. • Improvement heuristics: use a mathematical programming model to improve a solution found by a different heuristic approach. • Column generation-based approaches: Branch-and-price algorithms. Such algorithms use a set partitioning formulation by associating a binary variable to each possible route (column). Due to the exponential number of variables, the solution of the linear relaxation of the formulation is performed through column generation.
Our approach can be considered belonging to class 2 with a substantial modification: we designed the heuristic in order to give more importance to the optimization of the objective function at the expense of its feasibility. The restoration of the feasibility is entrusted to the mathematical model. In detail, our math-heuristic (MAT) approach uses the mathematical model to solve a "reduced" version of the instance to be solved.
For each instance to be solved, we perform the iterated local search (ILS) technique several times, using a multi-start approach, and analyze all the different solutions obtained. As a result of this process, we try to identify how often each arc is present or not in the generated solutions. In this way, we reduce the size of the instance by eliminating arcs that have a low probability of being part of the optimal solution and leaving (promising) arcs that have a good probability of being part of the optimal solution. To implement this reduction, for each arc, we count the number of times it is present in the final ILS solution: those arcs that exceed a predetermined threshold are left in the graph, while those below this threshold are eliminated. Obviously, this reduction process is carried out by preserving the connection of the graph. In the implementation of the proposed technique, It is chosen a dynamic threshold for each node, such that at least five incoming and five outgoing arcs are present for each node of the graph after the reduction.

Memetic approach
An extensive review and analysis of the main meta-heuristics applied to the various VRP problems were presented in the paper by Vidal et al. (2013Vidal et al. ( , 2019. The meta-heuristics used to solve VRP and which have given the best results are Tabu search, simulated annealing, deterministic annealing (Li et al. 2012), genetic algorithms (Minocha et al. 2011), ant colony (Donati et al. 2008), and variable neighborhood search implemented by Hemmelmayr et al. (2009).
The memetic approach we designed is obtained by embedding our simple local search (LS) in a genetic framework. In particular, we apply the LS to each chromosome that our GA produces. We use this memetic approach instead of a classic GA to guarantee the feasibility of the solutions generated. In fact, constraint (9) in the model (3)-(12), which limits the maximum length of each route, makes it very hard to get feasible solutions. For this reason, the crossover and mutation operators we designed are highly destructive and, in general, do not guarantee the feasibility of the solutions produced. Then, we combine the power of the GA that tries to improve the solutions without caring too much about the Each solution S is composed of a set of routes R. A vehicle and a sequence of nodes representing the path are associated with each route. The first element of each route is always the depot d. An example of a solution to our problem is reported in Fig. 4.
In the remainder of this paper, we denote | S | the number of routes contained in a solution S, and we use | R | to indicate the number of customers in the route R.

Constructive algorithm
The purpose of this algorithm is to create an initial solution (not necessarily feasible) to be used as input for the successive approaches. Until all customers are not served, this algorithm: (i) creates an empty route; (ii) associates a vehicle with the route; (iii) adds other customers to the route until the vehicle has residual capacity. The route always starts from the depot. Then, iteratively, the customers are added, choosing among the not served customers one, close enough to the last customer reached. Note that this scheme does not take into account constraints (9) of the MILP model, while constraints (7) are considered in the construction phase of the solution. A detailed pseudo-code of the algorithm is shown in Algorithm 3.
This algorithm creates a random solution S. As it is possible to notice for the schema reported in Algorithm 3, this solution depends on the parameter px. More precisely, with px = 1, we have a solution in which, for every route, the Nearest Neighbor policy is used. Increasing the value of px, we can get very different solutions from each other. Since this procedure takes into account only the capacity of the vehicle and not the maximum travel time of a route, it is possible to find unfeasible solutions to the problem. It would be easy to remedy the problem, but since the solutions produced are only useful as initial solutions for the algorithms proposed in the next sections, a degree of unfeasibility has proved to be experimentally useful to guarantee a wider exploration of the feasible region.

Fitness function
In order to compare two solutions and determine the best, it is usually sufficient to use the objective function associated Remove a random vehicle v from V . 5: Create an empty route R.

6:
Add R in the solution S. 7: Set v in R. 8: Add the depot d to R. 9: while U = ∅ do 10: let n the last node added to R. 11: n next = Get N extCustomer(U , v, n, px) 12: if k v ≥ q n + i∈R q i then 13: Add n next to R. 14: Remove n next from U , 15:  Example of neighborhood Swap, the customer n 2 ∈ R 1 is exchanged with the customer n 8 ∈ R 2 with the problem. Unfortunately, when algorithms have to compare solutions that are not always feasible, the use of the objective function alone is not enough. In fact, the objective function value may not be calculable in the case of an unfeasible solution, or this value may not be representative of the goodness of the solution.
In this paper, we use the objective function as a fitness function in the case of a feasible solution. In the case of an unfeasible solution, we add a penalty to the fitness function useful to represent the degree of the infeasibility of the solution.
Given the solution S, S z represents the objective function value of the solution. In the case of an unfeasible solution, S 6 represents the total violation in terms of a vehicle's maximum capacity (violation of constraint 6), and S 8 represents the total violation in terms of routes length (violation of constraint 8).
The fitness function f f computed for solution S is equal to:

Local search
A local search algorithm was developed to make the solution proposed by the constructive approach feasible and improve its quality. The local search is a technique that makes use of the neighborhood of a solution. Starting from an initial solution, this is iteratively improved by exploring its neighborhood until the technique reaches a local minimum. In order to describe the implemented algorithm, it is necessary to introduce the neighborhood used. In our implementation, we used four different neighborhoods combined in sequence. Before describing the algorithm, we introduce the used neighborhood.
Neighborhood 1: One-Change. For each route R ∈ S, for each customer n ∈ R, n is removed from R and reinserted in R in every possible position. This neighborhood will therefore contain a number of solutions equal to | S | * | R | * (| R | −1).

Neighborhood 2: 2-OPT.
For each route R ∈ S, for each couple of arcs (n i , n j ), (n k , n h ) ∈ R, the arcs (n i , n j ) and (n k , n h ) are removed from R, the path from n j to n k in R is reversed and the arcs (n i , n k ), (n j , n h ) are added to R to close the route. This neighborhood will therefore contain a number of solutions equal to | S | * | R | * (| R | −1). Figure 5 shows an example of this operator.
Neighborhood 3: Swap. For each pair of routes R 1 , R 2 ∈ S, for each pair of customers n 1 ∈ R 1 , n 2 ∈ R 2 , in the route R 1 the customer n 1 is replaced by n 2 and in the route R 2 the customer n 2 is replaced by n 1 . This neighborhood will therefore contain a number of solutions equal to | S | * (| S | −1) | R 1 | * | R 2 |. Figure 6 shows an example of this operator.
Neighborhood 4: Migration. For each pair of routes R 1 , R 2 ∈ S, for each customer n 1 ∈ R 1 , n 1 is removed from R 1 and added to R 2 in every possible position. This neighborhood will therefore contain a number of solutions The local search algorithm is summarized in the flowchart shown in Fig. 7.

Genetic algorithm
Genetic algorithms are powerful metaheuristic algorithms, initially introduced by Holland (1975) and largely used to solve complex combinatorial optimization problems. Holland was inspired by Darwin's theory of evolution, based on the idea of natural selection. Genetic algorithms somehow reproduce the mechanisms which characterize the evolution of life forms. Basically, they simulate the evolution of a population of individuals, each of which represents a feasible solution to the problem. The algorithm, starting with an initial population, performs a certain number of iterations, producing at each step a new generation of individuals created on the basis of the previous one. In this section, we describe a genetic algorithm for our problem. First of all, we need to specify what is the form of a solution, or rather how we encode the chromosome which represents a solution. Then, we can define the fitness function and the stopping criteria and describe how the selection, crossover, and mutation operations are performed.

The chromosome
The purpose of the chromosome is to represent a generic solution to the problem. The genetic algorithm uses the same  Table 2 Random choice of the assignment of vehicles to each node representation of the solution introduced in Sect. 4 and shown in Fig. 4.

The fitness function
The fitness function is used to associate a numerical value to each chromosome. This value is essential to be able to compare the chromosomes. Usually, this function is a direct consequence of the objective function of the problem, but it can manage the unfeasibility of some solutions through an appropriate penalty. The genetic algorithm described in this work uses the fitness function presented in Sect. 4.4.

The crossover
In a genetic algorithm, the crossover is usually used to create individuals of generation (i + 1) starting from some individuals of generation i. Our implementation of the crossover requires 2 chromosomes called parents and produces as output a new chromosome called son. Our crossover procedure is very simple; it starts from the parent chromosomes C 1 , C 2 and creates the child chromosome C 3 by iterating over all the nodes of the graph excluding the depot. For each node n i the vehicle that will serve it in C 3 is chosen randomly between the vehicle that serves n i in C 1 or C 2 . Figure 8 shows an example of input chromosomes for the crossover. Table 2 shows for each node, the vehicle that serves it in C 1 and C 2 . Figure 9 shows the possible crossover resulting by choosing randomly from Table 2 a vehicle for each node.

The mutation
The mutation procedure used is very invasive but, at the same time, elementary to describe. For √ | N | + 1 nodes, a new vehicle is randomly set among all those available. All routes associated with nodes for which the vehicle has been  Table 2 modified are shuffled. The reasons for using such a strong mutation will be explained later.

The initial population
The initial population is created by using the constructive algorithm (Algorithm 2) introduced in Sect. 4.3, setting the parameter px = 5. The population is composed of 100 chromosomes. The size of the population, as well as all the other parameters of the genetic algorithm, has been chosen following an intense manual tuning phase.

The Genetic Algorithm scheme
The GA starts from an initial population P created as described above. Then, two chromosomes C 1 , C 2 are chosen using a tournament selection. The tournament selection "TournamentSelectionGood" requires as input the population P and an integer denoted el. The procedure randomly selects a number of chromosomes from P equal to el, and returns the one with the best fitness function. The next step involves creating the C 3 chromosome by applying the crossover operator. The chromosome C 3 is improved by applying the previously described local search algorithm. If C 3 is already present in the population, this chromosome is reconstructed using the same technique already used for the creation of the initial population. If C 3 is not present in the population, its fitness is compared to the median fitness of the population. If the fitness of C 3 is greater than the median fitness of the population, the mutation operator is applied to chromosome C 3 . Otherwise, we proceed by selecting a chromosome to be removed from P to empty a slot for C 3 . C 4 is the chromosome to be removed from the population. It is selected through the "TournamentSelectionBad" procedure which behaves extensively like the previous tournament selection except for the fact that this procedure chooses the chromosome with worst fitness value. The procedure restarts again until the stopping criterion is reached; otherwise, the chromosome in the current population whose fitness value is the best (the minimum) is returned. This genetic algorithm scheme has the characteristic of not making the population evolve for generations but of allowing a constant evolution of the current population. To forcefully introduce the concept of generation into the population, we can say that the population of generation i differs from that of generation i + 1 by exactly two chromosomes. The operating scheme of G A is summarized in Fig. 10. The extreme randomness of the operators used within this genetic algorithm is due to the

Computational experiments
In this section, we will report the experimental results of the tests performed. For these experiments, we used input instances of three types: (i) instances from paper (Cerulli et al. 2018), (ii) instances randomly generated, (iii) instances related to the city of Genoa. All the code has been implemented in Java, using CPLEX version 12.8 as a solver (with the default setting on all parameters and without thread limit). All tests were performed on a laptop with the following hardware features: MacBook Pro, processor Intel i9 2.9GHz, 32GB RAM. For all the experiments, we set the following operating scenario: 2 h time window, vehicle speed 40K m/h, 7 min service time. This scenario implied that L v = 80K m, L c = 5K m. The local search algorithm was iterated 100 times (the same value as the GA population), and the best objective function value and the sum of the running times of the 100 iterations is shown in the tables reported in the next subsections.
The first set of instances used for our experiments is described in the paper of Cerulli et al. (2018) and is representative of a real problem of distribution of goods within the city of Genoa. Thanks to our mathematical model, which makes use of the reduction phase of the input graph, we were able to solve optimally all the input instances in a computation time less than 1 s, compared to the minutes or hours needed for the model proposed in the original work. For this reason, we realized that in order to test the developed approach, larger instances had to be used. In the next sections, we will introduce these instances describing their characteristics.

Random instances
This new set of instances was generated randomly; the motivation is that we wanted to carry out test cases useful to validate different characteristics of our algorithms. The approach used to create the instances can be schematized as follows: (i) place all customers within a 1000 × 1000 size grid randomly; (ii) the distance between customers is set equal to the respective Euclidean distance; (iii) the cost of each arc is set equal to its length multiplied by a value in the range [0.216: 0.432] randomly chosen for each type of vehicle (see Fig. 1); (iv) the demand of each customer is set to a random value in the range [1, 10]; (v) for each type of vehicle, the capacity is set randomly so that the sum of all the capacities of all the vehicles is double compared to the total customer demand. This leads to the generation of instances for which at least half of the vehicles must be selected. Two sets of instances were generated using this algorithm, refer-   Tables 3  and 4, respectively. For each scenario, five different instances have been created; the values reported are the average of the generated instances. With reference to the medium size instance (see Table  3), the proposed model (3)-(12) is able to obtain the optimal solution in less than 1 h of CPU time in scenarios with up to 4 vehicles and a maximum of 20 customers. The LS algorithm is very fast and produces good solutions but, unfortunately, sometimes fails to identify a feasible solution (see the cases highlighted in bold in Table 3). Although the GA algorithm is computationally more expensive, it always finds a feasible solution, which is almost always the optimal solution; when the model fails to find the optimal solution, GA improves in a few seconds the results obtained by the model in 1 h. The math-heuristic in a computational time lower than the mathematical model produces good quality results. In two scenarios, MAT produces better solutions than GA.
In the case of large instances (see Table 4), the genetic algorithm produces very good quality solutions in a maximum CPU time of about 2 min. LS is very fast but, in cases where it is very difficult to find a feasible solution, often it remains trapped in an unfeasible local minimum without being able to provide a solution (see the three cases having the value "No" in Table 4). Since we noticed that the proposed LS algorithm is very fast and produces good solutions in only 100 iterations, as a further test, we wanted to verify what would happen by increasing this number to a higher running time than GA. The answer is that at 100,000 iterations, the LS becomes slower than GA and produces many unfeasible solutions (12), never reaching the quality of GA's solution. These results are summarized in Table 5.

Genoa instances
This set of instances is composed of graphs extracted from the road network of the city of Genoa. To create the graphs,  Table 3 for the meaning of the columns Using the instances shown in Table 4 we used the data of OpenStreetMap contributors (2017). We generated directed graphs, taking into account the direction of the roads. For each vertex, we have the corresponding (x, y) coordinate in a 2D plane. Each arc of the graph is associated with its length on the actual street network. All the graphs are strongly connected. The generation of the instances uses the algorithm described in Capobianco et al. (2017) and used in Cerrone et al. (2019). The computational results on this set of instances are shown in Table 6. The complete graph of the instance of the city of Genoa consists of 8660 vertices and 16158 arcs. Note that this graph has been obtained by applying the graph reduction algorithm described in Sect. 2.2. From Table 6, readers can note that for scenarios up to about 20 customers, the proposed MILP model is able to find the optimal solution in an acceptable CPU time. The local search algorithm cannot get feasible solutions for the largest instances (6 vehicles, 63 and 66 customers), while the results reported in Table 6 confirm the good behavior of the GA algorithm. It is worth remembering that in all instances, in the graph representing the city of Genoa, the real distances of the arcs have been weighed not only with external costs due to environmental factors (see Sect. 2.1) but also with a penalty due to the transit of commercial vehicles near schools or other places of interest (see formula 2). For this reason, we wondered if the lower bound computed by using Proposition 3 (see Sect. 2.3) would produce a solution having the optimal value of the objective function lower than that found by the model (3)-(12) in the first three instances, for which we computed the optimal solution. Fortunately, the lower bound given by applying Proposition 3 gave the same solution for the three instances as the MILP model. This means that in these cases, the reduction algorithm did not remove the optimal solution by reducing the size of the graph. Figure 11 reports the graph of the city of Genoa. The green nodes represent the customers' positions. The single black node represents the depot. The red points indicate the locations of the schools within the city. Before applying our algorithm to these instances, all the red nodes have been removed, and all the arcs within a distance less than 100Mt from these nodes have been marked as "not recommended" by setting the corresponding multiplier δ = 2 instead of δ = 1 for all other arcs.

Conclusion
In this paper, we focused on a variant of the Green Vehicle Routing Problem (GV R P) aimed at defining grocery distribution plans in urban areas, according to a vision of sustainable mobility, which is increasingly necessary today. In addition to the typical constraints of the problem, we included in the objective function both travel costs and externality costs due to environmental factors, which depend both on the type of vehicle used and the type of road crossed. Further, the arcs of the graph near schools, parks, or other places of interest have been associated with a penalty that increases their cost, thus encouraging the search for alternative routes. Finally, the delivery plans are defined considering only the time slots of a day in which in the city of reference is allowed the transit of commercial vehicles; in particular, from 2 to 4 h time slots are considered, including the service time. We presented a MILP model and two heuristic approaches, a math-heuristic and a memetic algorithm, both using a local search procedure. In particular, we imposed a constraint that limits the maximum length of a tour traveled by any vehicle to face two different issues on the urban freight distribution: the transit of commercial vehicles in the city and, in this period of the COVID-19 pandemic, the possibility of sanitizing the vehicle (and drivers) several times during the day. The reported computational experiments based on random instances specifically tailored for the problem show the efficiency of the proposed solution method. The validation of the proposed algorithms performed with the graph of the  | V | indicates the number of available vehicles. For this experiment, we have two vehicle classes, LCV and HGV. | U | indicates the number of customers. For the MILP model, for the local search (LS), for the genetic algorithm (GA), and the metaheuristic (MAT), both the value of the objective function and the running time are shown. Times are shown in seconds; for the MILP and for MAT, the maximum running time is limited to 3600 and 600 s, respectively. L B represents the value produced by the linear relaxation of the model city of Genoa allows us to say that the proposed algorithms are a valuable aid in the definition of sustainable distribution plans in urban areas, also considering that the proposed model is easily adaptable to different urban realities and mobility rules.
In the future, we would like to investigate the proposed math-heuristic approach further. In particular, we want to make the mathematical model "tighter" and efficient and want to deepen the ways of reducing the size of the input graph based on the solutions proposed by the local search procedure.
Funding Open access funding provided by Università degli Studi di Genova within the CRUI-CARE Agreement. The authors did not receive support from any organization for the submitted work Data Availability The datasets analyzed during the current work are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Ethics approval This article does not contain any studies with human participants or animals performed by any of the authors.

Code availability
Code developed for the model will be made available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.