1 Introduction
Transportation is an important factor for today’s smart society. Public transportation, for example, is used by almost the whole population, and it directly affects the quality of life. However, modelling and planning such complex transport system is a very challenging task. Here, this paper focuses on another sort of transport: the transportation in the business world. Due to the rapid advancement of technologies, logistics systems have become very important for media companies. The fact that anyone in the world can be well connected has led to complex transport networks that are very demanding and are becoming increasingly important. Therefore, an efficient logistic network can make a huge difference for companies and relevant business operations.
To cite one fact that highlights the importance of logistics in this sector, in some businesses, such as groceries delivery, the distribution costs can lead to an increase in the product price up to 70% [1]. Thanks to cases like this, it is obvious to show the importance of this sector.
Therefore, this paper focuses on a realworld logistic problem, and its effective treatment. The realworld situation dealt in this work is related to the daily delivery of newspapers. Specifically, the object of this study is a mediumsized newspaper distribution company, which offers its services at the regional level. This company has certain mandatory principles, which must be taken into account when performing the daily planning tasks of deliveries. One of these principles is a strict recycling system. Another one is the treatment of the different towns, or cities, as separate units. Besides that, the company considers certain factors for the scheduling process, as variable travel times (depending on the hour they are carried out), or the transit prohibition in certain streets. Although this paper is focused on a specific company located in Bizkaia (Spain), it is noteworthy that the goal of this work is to develop a model that is applicable to every company of similar characteristics.
Therefore, the aim of this paper is to address efficiently this Newspaper Distribution System with Recycling Policy (NDSRP). For this purpose, the NDSRP has been modeled as a MultiAttribute, or RichVehicle Routing Problem (RVRP). Nowadays, as indicated in [2] or [3], RVRPs form a hot topic in the scientific community. These kinds of problems are special cases of the wellknown conventional vehicle routing problem (VRP) [4], with the distinction of having multiple constraints and complex formulations. This sort of problems can have a great scientific interest, since such NPHard problems present a tough challenge to solve. Furthermore, their social interest is also high, as their applicability to realworld situations is greater than the conventional versions of routing problems.
To be more specific in the present paper, an Asymmetric and Clustered Vehicle Routing Problem with Simultaneous Pickup and Deliveries, Variable Costs and Forbidden Paths (ACVRPSPDVCFP) was proposed to address the presented NDSRP. Some examples of recently developed RVRPs can be [5] or [6]. In the former work, the authors propose an RVRP with hard and soft time windows, heterogeneous fleet, customers priorities and vehiclecustomer constraints. The solution proposed by the authors of this paper has already been integrated into the optimization tool of a fleet management system used in the Canary Islands. In the latter work, an RVRP was proposed to deal with the perishable food management. In this work, a heterogeneous fleet sitedependent VRP with multiple time windows was presented. Another example of recently developed RVRP can be the proposed in [7]. In their paper, an RVRP was developed for the olive oil collection in Tunisia. The RVRP designed in this work was a multiproduct, multiperiod and multicompartment VRP. These were some of the examples that justified the increasing interest of RVRP in the scientific community. For further information on RVRP, readers can refer to the surveys [2] and [3].
Though some appropriate methods can be found in the literature to address such complex problems, arguably the most successful techniques to solve these RVRPs are the heuristics and metaheuristics [8]. In this study, the attention has been focused in the last ones. In fact, to solve the ACVRPSPDVCFP proposed in this paper one metaheuristic has been developed. Some classical examples of metaheuristics are the Tabu search [9] and simulated annealing (SA) [10]
as local search methods, and genetic algorithm (GA)
[11, 12], ant colony optimization [13]and particle swarm optimization
[14] as population ones. Though such methods were proposed some decades ago, they still attract the attention in the scientific community [15, 16, 17] and metaheuristics, especially new metaheuristic algorithms and their proper implementations still form a hot topic in the field. In fact, many different metaheuristics have been proposed in the last two decades, which have been successfully applied to various problems and fields. Some examples of these techniques are the imperialist competitive algorithm, presented in 2007 by Gargari and Lucas [18], the bat algorithm, proposed by Yang in 2010 [19], or the harmony search, presented in 2001 by Geem et al. [20].Another type of technique that shows a good performance applied to routing problems are the memetic algorithms [21, 22]. In [23], for example, a hybrid genetic algorithm is presented for solving a large class of vehicle routing problems with timewindows. On the other hand, in [24] a hybrid genetic algorithm is presented to solve an RVRP composed by multidepot and periodicity features. In [25] a decomposition based memetic algorithm is proposed for a multiobjective VRP with time windows. Another kind of RVRP is solved by a hybrid techniques in [26]. In this case, the characteristics of the problem are clustered backhauls and 3D loading constraints. Finally, a multiperior VRP with profit is addressed by a memetic algorithm in [27]. Works cited in this paragraph are some recent and interesting examples of the whole literature, some other recent examples can be found in, for example, [28, 29, 30].
Therefore, in this paper one metaheuristic proposed a few years ago is used to solve the presented problem. This technique is the firefly algorithm (FA), proposed by Yang in 2008 [31]. This natureinspired algorithm is based on the flashing behaviour of fireflies, which acts as a signal system to attract other fireflies. As can be seen in several surveys [32, 33], the FA has been successfully applied to many different optimization fields and problems since its proposal. In addition, it still attracts a lot of interests in the current scientific community [34, 35, 36]. Nevertheless, the FA has never been applied to any RVRP. This lack of works, along with the growing scientific interest in bioinspired algorithms, and the good performance shown by the FA since its proposal in 2008, has motivated its use in this study.
In this way, the main novelties and contributions of this paper can be listed in the following way. It is noteworthy that these originalities and contributions have motivated the realization of this work:

In this paper, a DPRP is dealt with using an RVRP. As will be mentioned in the following section, the problem of newspaper distribution has been previously addressed in the literature, but never using an RVRP as complex as the one presented in this paper. Furthermore, it is the first time that an RVRP of this characteristics is proposed in the literature. For this reason, a benchmark composed by 15 instances has been developed, based on real geographic locations.

In this work, a discrete version of the FA (DFA) is presented for solving the proposed RVRP. Until the time of writing, the FA has never been applied to any RVRP. Anyway, the main novelty of the presented DFA is not only its application field. The technique developed in this work uses the Hamming Distance function to measure the distance between two fireflies of the swarm. This approach has been rarely used previously, as well as the move strategy used by the fireflies, which is based on evolution strategies. All these characteristics are described in following sections. Additionally, in order to prove that the DFA is a promising technique to solve the designed ACVRPSPDVCFP, the results obtained by the DFA are compared with the ones obtained by the evolutionary algorithm (EA) based on mutations and the evolutionary simulated annealing (ESA) [37].
The rest of the paper is organized as follows. In Section 2, the realworld problem addressed in this work is described. In Section 3, the RVRP proposed to deal with the real problem is described in detail. Furthermore, in Section 4, the developed DFA is described in detail. After that, the experimentation carried out is detailed in Section 5, as well as the benchmark designed for the presented RVRP. This paper then concludes with some brief conclusions and further research topics in Section 6.
2 Newspapers Distribution with Recycling Policy
As has been pointed in the previous section, the realworld situation faced in this work is related to the newspaper distribution. More specifically, the object of study is a mediumsized newspaper distribution company. The area of the coverage of this company is at a provincial level. This means that the company has to serve a set of customers geographically distributed in separate towns and cities. The company in question has some principles, which are the base of their logistics planning. The first principle is to treat towns, or cities, as separate units. In this way, if one vehicle enters a city, or a town, it is forced to serve each and every customer located therein. Therefore, one vehicle is banned for entering one city or town whether it does not have sufficient capacity to meet the demand of all customers deployed there.
On the other hand, due to the current environmental requirements, the company has a simple but robust paper recycling policy. In this case, the objects to recycle are the newspapers that were not sold the previous day. Thus, as can be deduced, vehicles not only have to meet the delivery demands of the customers. Besides that, they have to collect at each point those newspapers that have not been sold the day before.
In addition, the company has to take into account certain factors in the routes planning process. The first one is related to the hours at which the deliveries and collections are done. The service is performed daily during morning from 6:00am to 15:00. Within this time window exists one range, established between 8:00am and 10:00am, considered as “peak hours”. In this way, traveling costs from one point to another are greater if they are performed at “peak hours”. Besides this, with the aim of respecting all the traffic rules, vehicles cannot go through prohibited roads.
Throughout the past few decades, the problem of the newspaper delivery has been addressed many times in the literature. In 1996, Ree and Yoon presented a twostage heuristic for a newspaper distribution problem, taking as a case study of a major press company of Korea [38]. The problem used in this study is a multiVRP with time windows, which is faced in two different stages by the proposed heuristic. On the other hand, in [39] a free newspaper delivery problem is addressed. In their study, the problems was decomposed in two phases. In the first stage, the delivery plan was created, based on concept taken from the Inventory Routing Problem [40]. Once the delivery plan was fixed, the resulting problem was solved in the second stage as a variant of the conventional VRP with time windows (VRPTW) [41]. They used a heuristic for solving this second stage problem. The same VRPTW was used in [42] to face a case study of the newspaper delivery problem focused in the city of Bangkok and a modified sweep algorithm was presented to solve the proposed VRPTW. These papers are some of the examples that can be found in the literature. Many others can be found, considering, for example, the problem of delivering newspapers to private subscribers [43, 44].
Regarding the problem addressed in this paper, some originalities and novelties are proposed: the features of the variable travel times and forbidden paths have never been used for the newspaper distribution system, as well as the recycling policy applied in this work. In addition, as we will see in the rest of the paper, the RVRP proposed in this work has a great number of constraints, making easier its application to the real world. Finally, the technique developed in this work solves the proposed ACVRPSPDVCFP in only one phase, in contrast with the approaches presented in [38, 39]. This is an advantage because solving the problem in only one phase, the runtimes and the computational effort can be considerably reduced.
3 Asymmetric and Clustered Vehicle Routing Problem with Simultaneous Pickups and Deliveries, Variable Costs and Forbidden Paths
As has been stated earlier in this paper, the realworld situation of the newspaper distribution has been modeled as an RVRP. In this section, a detailed description of the presented problem is provided. First of all, in Section 3.1, the basic features of the problem are detailed, followed by the mathematical formulation of the proposed RVRP is depicted in Section 3.2.
3.1 General description of the problem
The proposed RVRP has the following general characteristics. It is noteworthy that all these features proposed here are to take into account the conditions stated in Section 2. Besides those, some additional restrictions have been considered with the intention of developing a model closer to real conditions.

Asymmetry: The traveling costs in the proposed RVRP are asymmetric. This means that the traveling cost from any node to another node is different from the reverse trip cost. It is important to highlight that, in the problem proposed in this paper, this asymmetry appears in every nodetonode travel. This feature is not common in most routing problems that can be found in the literature, and it brings realism to the problem. Anyway, asymmetric costs have been applied previously in the literature in a different context [45, 46, 47]. It is noteworthy that this feature is very valuable in realworld applications.

Clusterized: With this attribute, the different clients that make up the system are grouped in different sets or clusters. In this case, each cluster represents a city or town. The condition that vehicles must meet is the following: if one vehicle meets the demand of any customer belonging to any cluster, this vehicle must serve each and every customer of this cluster. In this way, one vehicle cannot enter a cluster if it does not have enough capacity to serve all customers in this city or town. This feature has been used previously in many studies in a different context. [48, 49, 50].

Simultaneous Pickup and Delivery: This property is an adaptation of the often used pickup and delivery system of some routing problems [51, 52]. Basically, this system consists in the existence of two types of nodes: the delivery nodes and the collect nodes. The former ones are those points in where newspapers are delivered. On the other hand, in collect nodes, newspapers are collected, with the intention of bringing them back to the center for their subsequent recycling.
In addition, it is important to highlight that, due to the simultaneous nature of the realworld situation, one customer can ask for both delivery and collection of newspapers. In this way, both delivery nodes and deliverycollect nodes can be found in the system. On the other hand, it is assumed that all customers request the delivery of newspapers. For these reason, customers which demand only the collection of newspapers have not been taken into account.

Variable Travel Times: In real transportation situations, the travel between two points does not always take the same cost, either temporarily or economically. In many cases, this cost is subject to some external variables, such as the hour, the traffic or the weather. With the intention of adding more realism to the problem, this situation has been represented in the RVRP proposed in this work. To this end, it has established a schedule between 6:00am and 15:00. Within this schedule, two different timeperiods have been established, called “peak hours” and “offpeak hours”. Peak hours are composed by two hours, between 8:00am and 10:00am. All trips performed in this time window imply a higher cost, comparing with the costs of conducting the same trip in “offpeak” period. A similar characteristic has been previously used a few times in the literature [53].

Forbidden Paths: In the real world, it is common to find oneway roads, where the traffic in a particular direction is prohibited. There are also pedestrian streets, where vehicles cannot go through. With the aim of recreating this common situation, the problem has certain arcs which cannot be used in the final solution. A similar philosophy has been used previously in some studies in the literature [54].
With the above assumptions and simplifications, the proposed ACVRPSPDVCFP is an RVRP, whose objective is to find a set of routes, trying to minimize the total cost of the complete solution, and taking into account the two different kinds of nodes, respecting the restrictions of the clusters and vehicles capacities () and not traveling through any forbidden path. As can be seen, being an RVRP problem, the ACVRPSPDVCFP has multiple constraints. It is important to point out that all these restrictions reduce the size of the search space which encompasses the feasible solutions. Anyway, all these constraints make the process of generating feasible solutions and successors to be a very complex task.
3.2 Mathematical description of the problem
The presented ACVRPSPDVCFP can be defined on a complete graph where is the set of vertices which represent the nodes of the system. On the other hand, is the set of arcs which represent the interconnections between nodes. Each arc has an associated cost . Due to the asymmetry feature, . It is noteworthy that, in this case, the cost of traveling from to is always different from the cost of traveling from to . The cost of traveling via a forbidden arc is infinite. In this way, it is ensured that they will not appear in the final solution. Furthermore, the vertex represents the depot, and the rest are the visiting points. Additionally, is divided into mutually exclusive nonempty subsets, , each one for each cluster. These subsets hold the following conditions:
(1) 
(2) 
It is noteworthy that only contains , which depicts the depot. The remaining nodes are divided into clusters. Besides that, each node has assigned two kinds of demands: one associated with the delivery 0 and the other with the pickup 0.
The proposed ACVRPSPDVCFP can be mathematically formulated in the following way. It is important to highlight that depicts the demand pickedup in clients routed up to node (including node ) and transported in arc . On the other hand, represents the demand to be delivered to clients routed after node and transported in arc [55]. Furthermore,
, is a binary variable, which takes 1 value whether vehicle
enters the cluster , and 0 in other case. Finally, the binary variable is 1 if the vehicle uses the arc , and 0 otherwise.The main problem is now to minimize:
(3) 
where
(4) 
(5) 
(6) 
(7) 
This is subject to the following constraints:
(8) 
(9) 
(10) 
(11) 
(12) 
(13) 
(14) 
(15) 
(16) 
(17) 
(18) 
(19) 
The first clause represents the objective function, which is the sum of the costs of all routes of the solution, and it must be minimized. The formulas (4), (5), (6) and (7) depict the nature of the variables , , and . Equations (8) and (9) assure that all the nodes are visited exactly once. On the other hand, constraints (10) and (11) ensure that the total amount of vehicles leaving the depot, and the number of vehicles that return to it is the same. Besides, the correct flow of each route is ensured by functions (12) and (13).
Additionally, restrictions (14) and (15) guarantee that the flows for the collection and delivery demands, respectively, are properly conducted. These formulas ensure that both demands are satisfied for every customer. Furthermore, constraint (16) assures that the total capacity of any vehicle is never exceeded; it also establishes that pickup and delivery demands will only be transported using arcs included in the solution [55].
On the other hand, functions (17) and (18) ensure that every trip between one node and another has not an infinite cost. Thus, it is guaranteed that forbidden paths will not form part of the final solution. Finally, restriction (19) assures that only one vehicle enters every cluster. This function, together with the above mentioned (8) and (9), ensures that all the customers belonging the same cluster are visited by the same vehicle.
4 Firefly algorithm
As has been stated in the introduction, a discrete firefly algorithm (DFA) is proposed in this work to address the designed ACVRPSPDVCFP. In this section, the description of the classic FA is shown first (Section 4.1). Then, the proposed DFA is described in detail in Section 4.2.
4.1 Classic Firefly Algorithm
The basic FA was first developed by XinShe Yang in 2008 [31, 56], and it was based on the idealized behaviour of the flashing characteristics of fireflies. To properly understand the algorithm, it is important to highlight the following three idealized rules [31]:

All the fireflies of the swarm are unisex, and one firefly will be attracted to other ones regardless of their sex.

Attractiveness is proportional to the brightness, which means that, for any two fireflies, the brighter one will attract the less bright one. The attractiveness decreases as the distance between the fireflies increases. Furthermore, if one firefly is the brightest one of the swarm, it moves randomly.

The brightness of a firefly is directly determined by the objective function of the problem under consideration. In this manner, for a maximization problem, the brightness can be proportional to the objective function value. On the other hand, for a minimization problem, it can be the reciprocal of the objective function value.
In Algorithm 1, the pseudocode of the basic FA is shown, which was proposed by Yang in [31]. In line with this, there are three important factors to consider in the FA: the attractiveness, the distance and the movement. In the basic version of the FA these factors are tackled in the following way. First of all, the attractiveness of a firefly is determined by its light intensity, and it can be calculated using this formula:
(20) 
On the other hand, in the basic FA the distance between any two fireflies and is calculated using the Cartesian distance, and it is computed by the following equation:
(21) 
where is the th component of the spatial coordinate of the th firefly. Finally, the movement of a firefly toward any other brighter firefly is determined by this formula:
(22) 
where is the randomization parameter and
is a random number uniformly distributed in [0,1]. On the other hand, the second term of the equation stems from the attraction assumption.
4.2 The proposed Discrete Firefly Algorithm
It is noteworthy that the original FA was developed primarily for solving continuous optimization problems. For this reason, the classic FA cannot be applied directly to solve the proposed ACVRPSPDVCFP. Therefore, some modifications of the original FA are needed in order to prepare it for addressing such ACVRPSPDVCFP problem.
First of all, in the proposed DFA, each firefly in the swarm represents a possible and feasible solution for the ACVRPSPDVCFP. All the fireflies are initialized randomly. Additionally, as has been detailed in Section 3, the sum of the costs of all routes of the solution has been used as the objective function. Therefore, the ACVRPSPDVCFP is a minimization problem, in which the fireflies with a lower objective function value are the most attractive ones. In addition, the concept of light absorption is also represented in this version of the FA. In this case, = 0.95, and this parameter is used in the way as can be seen in Equ. (22). This parameter has been set following the guidelines proposed in several studies of the literature [56, 31].
Besides that, the distance between two different fireflies is represented by the Hamming Distance. The Hamming distance between two fireflies is the number of noncorresponding elements in the sequence. In the proposed ACVRPSPDVCFP, the comparison is made cluster by cluster. For example, taking into account two random fireflies, and one random cluster composed by 8 nodes:
the Hamming Distance between and for the cluster would be 4. This same comparison is made for every cluster. In this way, the total distance between fireflies and is the sum of all the distances for every cluster.
Finally, the movement of a firefly attracted to another brighter firefly is determined by
(23) 
(24) 
where is the Hamming Distance between firefly and firefly , and is the iteration number. In this case, the length of the movement of a firefly will be a random number between 2 and . As for the movement function, the Insertion Function has been used. This function selects and extracts one random node from a random route. After that, this node is reinserted in a random position inside the cluster of the selected node. This function takes into account the capacity constraint, in order not to create infeasible solutions.
Following the same philosophy as other previously developed and published FA for the Traveling Salesman Problem [57, 58], fireflies in the proposed DFA do not have directions to move. Instead, fireflies move using evolution strategies. In this way, each firefly moves using times the Insertion Function, generating potential successors. After these movements, the best one is performed, generating the new firefly.
The pseudocode of the presented DFA is shown in Algorithm 5. In lines 13 the initialization phase of the algorithm is carried out, in which fireflies are initialized and evaluated. Besides, the parameter is initialized as previously described. In addition, in lines 1012 the movement process is performed. In line 10 the distance between the selected and is calculated via Hamming Distance. Once the distance is obtained, the parameter is calculated, which is a random number between 2 and . Finally, the movement is performed in line 12 using the Insertion Function as explained before. After this movement process, fireflies are evaluated in line 14 and ranked in line 17. This iterative procedure is repeated until termination criterion is reached.
5 Experimentation
The computational experiments carried out in this study are described in this section. First of all, the details of the proposed benchmark for the ACVRPSPDVCFP are detailed (Section 5.1). After that, the results obtained by the developed DFA for this benchmark are presented (Section 5.2). In order to prove that the DFA is a promising metaheuristic for solving routing problems, the results obtained by the DFA have been compared with the ones obtained by the EA and the ESA. In addition, two different statistical tests have been conducted, with the aim of obtaining rigorous and fair conclusions (Section 5.3).
5.1 The proposed benchmark for the VRP
The type of RVRP proposed in this paper has never been treated before in the literature. It is for this reason that there is no benchmark available in the literature for the ACVRPSPDVCFP. In this work, a benchmark composed by 15 instances is proposed. At the same time, these instances are composed of 50 to 100 nodes. Every node represents a customer, and all nodes are placed in real geographical locations, which are located in the province of Bizkaia, Spain. In addition, the maximum number of clusters has been established as 10, existing also instances with 5, and 8 of them. In the Figure 2, a map with the geographical locations of the depot, the customers and the clusters are shown. This map has been made using Open Streep Maps technology, via uMap tool^{1}^{1}1http://umap.openstreetmap.fr.
The clusters have been organized in order of appearance. That is, the nodes 110 compose cluster 1, nodes 1120 for cluster 2, and so on. It is important to highlight that every cluster has the same number of customers in all instances. Besides that, as has been pointed out in Section 3.1, each customer has assigned two kinds of demands: one associated with the delivery and the other with the pickup . The assignments of these demands have been performed following this procedure:
(25) 
(26) 
(27) 
(28) 
where =0 and =0, taking into account that is considered the depot. In addition, the costs of traveling from any customer to other customer have been established following the procedure depicted in Algorithm 3. By this method, the asymmetry characteristic is met. It is important to highlight that these costs are assigned for the “offpeak” period. These costs are incremented when they are conducted on the “peak” period, following the procedure shown in Algorithm 4. In is noteworthy that, with the intention of simplifying the complexity of the problem, the traveling time between any node and is the same as its traveling cost (in seconds).
Finally, depending on the instance, some paths are chosen in each cluster to be forbidden. In Table 1, the characteristics of all the instances developed for the benchmark are summarized. In order to understand the content of this table correctly, the following clarifications should be made: Osaba_50_1_1 and Osaba_50_1_2 are comprised by 5 clusters, which are the clusters {1, 3, 5, 7, 9}. On the other hand, Osaba_50_2_1 and Osaba_50_2_2 are made up by clusters {2, 4, 6, 8, 10}. Moreover, clusters in Osaba_50_1_3 and Osaba_50_1_4 are composed by 5 nodes. In this case, this customers are the first 5 of every cluster. This is in contrast with what happens in Osaba_50_2_3 and Osaba_50_2_4, where the 10 clusters are comprised of the last 5 clients of every of them. Lastly, to complete all Osaba_80_X instances, the first 8 clusters, or nodes (depending on the case) have been chosen.With the aim of allowing the replication of this experimentation, it is noteworthy that the benchmark developed is available under request to the corresponding author of this paper, or via Web^{2}^{2}2http://research.mobility.deustotech.eu/media/publication _resources/Instances_Osaba_ACVRPSPDVCFP.rar.
Instance  Nodes  Clusters  Vehic. capacity  Forbidden paths 

Osaba_50_1_1  50  5  240  5 
Osaba_50_1_2  50  5  160  10 
Osaba_50_1_3  50  10  240  5 
Osaba_50_1_4  50  10  160  10 
Osaba_50_2_1  50  5  240  5 
Osaba_50_2_2  50  5  160  10 
Osaba_50_2_3  50  10  240  5 
Osaba_50_2_4  50  10  160  10 
Osaba_80_1  80  8  240  5 
Osaba_80_2  80  8  160  10 
Osaba_80_3  80  10  240  5 
Osaba_80_4  80  10  160  10 
Osaba_100_1  100  10  140  5 
Osaba_100_2  100  10  260  10 
Osaba_100_3  100  10  320  10 
5.2 Results
All the tests conducted in this work have been performed on an Intel Core i5 – 2410 laptop, with 2.30 GHz and a RAM of 4 GB. Java has been used as the programming language. All the 15 instances described in the previous section have been used in this experiment. Every instance has been run 20 times. As has been said before, the results obtained by the DFA are compared with the ones obtained by the EA and the ESA. The reason why these two techniques have been used for this experimentation can be summarized as follows: First of all, both metaheuristics are wellknown, and they have been frequently used to solve routing problems. Proving that the DFA outperforms these two techniques can be concluded that it is a promising technique to solve RVRPs. On the other hand, all these three techniques have two similarities: all of them base the movement of their individuals on a shortstep operator, and they are easy to implement and can be adapted to solve new problems.
EA  ESA  

Parameter  Value  Parameter  Value 
Population size  100  Population size  100 
Mutation functions  Insertion Function  Successor Function  Insertion Function 
Mutation prob.  1.0  Temperature  
Survivor func.  70% Elitist  30% Random  Cooling constant  0.95 
It is worth highlighting that, as far as possible, the same operators and similar parameters have been used for all the algorithms implemented for the experimentation. In this way, the intention is to conclude which algorithm obtains better results using similar operators and a similar number of times. Furthermore, with the intention of facilitating the replicability of this study, the parameters used for the EA and ESA are shown in Table 2. It is worth pointing out that all the individuals are randomly generated. Additionally, as for the termination criterion, every algorithm finishes its execution when there are function evaluations without improvements in the best solution, where is the size of the problem. In addition, the parameters used for DFA are the described in Section 4.2. In case of DFA, a population of 100 fireflies has also been used.
Finally, in this study, the permutation codification has been used for the representations of the solutions. Thus, each solution is encoded by an unique permutation of numbers, which represents the different routes that compose that solution. Besides that, with the aim of distinguishing the routes in one solution, they are separated by zeros.
Before starting the tests, a small study about the parametrization of the DFA is shown. It should be taken into account that performing a comprehensive study on the parameterization of the algorithm would be very extensive. That study has been planned as future work, since the objective of this paper is to present the NDSRP model, and to demonstrate that the DFA shows an adequate performance applied to an RVRP. For this reason, in this paper a small portion of that study is shown, in order to justify the population size used. For that, we compare 4 different versions of the DFA, each one with a different population size. Each version is called , where is the popualtion size. For this small test, 4 different instances have been used. In Table 3 the results of this experimentation are depicted. In this table, the results average (avg.), and average runtime (Time, in seconds) are shown. In addition, the last row of the table represents the average ranking for each alternative.
Instance  

Name  Avg.  Time  Avg.  Time  Avg.  Time  Avg.  Time 
Osaba_50_1_1  51945.7  14.8  51561.3  25.8  50989.5  37.9  50934.3  68.9 
Osaba_50_1_2  57398.7  15.8  56721.8  26.9  56203.8  35.1  56213.7  71.8 
Osaba_80_3  92990.39  33.7  91663.8  47.3  89512.0  75.3  89531.0  112.3 
Osaba_100_1  110206.94  51.4  108241.6  81.8  107799.5  112.8  107745.7  176.4 
Ranking  4  3  1.5  1.5 
Some conclusions can be drawn if the data presented in Table 3 are analyzed. It can be seen a slight trend of improvement in the results when the population size increases. Even so, this fact involves a significant increase of runtime, which is not directly proportional to the improvement in results. In this paper, to achieve the proposed, the option which best balances the runtime and results quality has been selected. This option is . This version has acceptable execution times, and is the alternative that more improvement offers regarding its previous version. , on the other hand, needs very high execution times without offering significant improvements in the results quality.
Once carried out this small study, the results obtained by the three techniques for the proposed benchmark are summarized in Table 4
. In this table, the results average (avg.), standard deviation (S. dev), average runtime (Time, in seconds) and average convergence time (C. T, in seconds) are shown. Additionally, the best results averages have been represented bolded. Besides that, the best results found for each instance are represented in Table
5. This table also shows the number of vehicles needed to perform every solution, and the technique which reached each of these results. Since this is the first appearance of the ACVRPSPDVCFP in the literature, these solutions are considered the best ones found until the publication of this paper.Instance  DFA  ESA  EA  

Name  Avg.  S. dev.  Time  C. T.  Avg.  S. dev.  Time  C. T.  Avg.  S. dev.  Time  C. T. 
Osaba_50_1_1  50989.5  234.9  37.9  26.2  51632.6  479.0  31.0  22.7  51569.3  649.7  38.4  26.8 
Osaba_50_1_2  56203.8  203.4  35.1  25.4  56929.4  459.1  31.7  21.0  57008.8  540.9  34.0  26.4 
Osaba_50_1_3  71730.0  1486.2  37.6  22.3  72298.3  1484.4  32.1  17.4  72490.7  1413.4  36.1  21.8 
Osaba_50_1_4  78883.8  1193.4  35.7  23.9  79168.0  1788.8  29.4  20.0  79207.9  1612.4  34.8  24.5 
Osaba_50_2_1  49276.0  392.4  38.6  22.3  49751.3  538.3  33.0  19.8  49416.8  760.6  36.9  20.4 
Osaba_50_2_2  54589.6  628.1  37.8  28.5  54871.6  621.0  31.8  23.7  55007.3  837.0  38.1  27.9 
Osaba_50_2_3  69631.4  2400.1  38.4  25.9  71352.6  2883.4  33.6  21.6  71286.3  3037.6  37.5  25.1 
Osaba_50_2_4  80543.7  1512.3  36.1  22.6  81122.0  1675.3  32.4  19.1  81297.3  1894.3  36.0  23.1 
Osaba_80_1  82307.8  1043.4  72.4  57.6  80834.2  1801.4  67.0  53.4  81779.6  2018.4  71.8  58.0 
Osaba_80_2  89324.9  698.0  74.0  60.7  89989.6  1134.4  67.6  54.3  90090.6  1032.5  73.4  59.4 
Osaba_80_3  89512.0  1414.0  75.3  60.4  89431.2  2680.1  68.5  55.2  89883.1  2949.6  75.6  59.8 
Osaba_80_4  104601.9  1299.2  75.0  64.8  105141.7  1963.3  69.4  61.3  106689.3  1832.9  74.8  65.3 
Osaba_100_1  107799.5  1501.8  136.1  112.8  109183.6  1540.4  129.3  108.1  109614.4  1700.7  140.7  116.7 
Osaba_100_2  100522.9  1683.0  138.0  109.2  101708.0  1644.1  130.4  102.0  101908.1  1699.0  141.0  111.9 
Osaba_100 3  95553.7  1470.6  139.7  105.3  95641.3  2687.3  128.0  98.4  95893.7  2472.0  140.3  108.4 
Name  Best Result  Vehicles  Technique 

Osaba_50_1_1  50641.46  2  DFA 
Osaba_50_1_2  55923.48  3  DFA 
Osaba_50_1_3  68535.19  2  DFA 
Osaba_50_1_4  76276.83  3  DFA 
Osaba_50_2_1  48819.71  2  DFA 
Osaba_50_2_2  53876.86  3  DFA 
Osaba_50_2_3  65059.03  2  DFA 
Osaba_50_2_4  77497.06  3  DFA 
Osaba_80_1  77660.84  3  ESA 
Osaba_80_2  87835.50  4  DFA 
Osaba_80_3  83713.72  3  ESA 
Osaba_80_4  101497.11  5  DFA 
Osaba_100_1  105165.35  5  DFA 
Osaba_100_2  97924.51  4  DFA 
Osaba_100 3  88966.48  3  ESA 
5.3 Analysis and Discussion
Analyzing the results in Table 4, the first conclusion that can be drawn is the following one: DFA outperforms clearly the other algorithms in terms of results. Specifically, the DFA performs better than ESA in 86.66% of the instance (13 out of 15), and in 93.33% of the cases comparing with the EA (14 out of 15). Additionally, the supremacy of DFA can also be seen in Table 5, having obtained the best solution in 80% of the instances (12 out of 15). In relation with the data shown in Table 5, in Figure 3 and Figure 4, the partial representation of the best solutions found by the DFA for the instances Osaba_80_2 and Osaba_100_1 are shown. It is noteworthy that, due to the complex nature of the problem, the complete representation of these solutions would lead to maps with some overlapping lines, complicating the visibility of the whole solution. Therefore, these solutions are shown at clusterlevel, showing also several clusters in detail.
Another important factor that is worth mentioning is the robustness of the DFA in relation to the other techniques. It should be clarified that the robustness is the capacity of a technique to obtain similar solutions in every run. As can be seen in Table 4, the standard deviation of the results obtained by the DFA is lower than the ones presented by the other metaheuristics in most instances (12 out of 15). This means that the quality of the solutions provided by the DFA move in a narrow range. This characteristic gives robustness and reliability to the algorithm, something crucial if the technique is applied a real environment.
Regarding the runtimes, in Table 4 it can be seen how DFA and EA have similar execution times, while the ESA takes less time than its competitors. These same conclusions can be drawn looking at the convergence times, where ESA shows that it needs less time to reach the final solution. This fact can be analyzed in the following way: The DFA needs more time to reach its final solution, showing a better exploration capacity than the ESA. On the other hand, it shows a better exploitation capacity than the EA, since it obtains better results needing similar execution and convergence times.
Besides that, two different statistical tests have been conducted in order to obtain rigorous and fair conclusions. It is important to clarify that the result averages obtained by each technique have been use to perform these tests. The guidelines given by Derrac et al. in [59]
have been followed to perform this statistical analysis. First of all, the Friedman’s nonparametric test for multiple comparisons has been used to check if there are any significant differences among all the techniques. The resulting Friedman statistic has been 17.73. Taking into account that the confidence interval has been stated at the 99% confidence level, the critical point in a
distribution with 2 degrees of freedom is 9.21. Since 17.73
9.21, it can be concluded that there are significant differences among the results reported by the three compared algorithms, being DFA the one with the lowest rank. Finally, regarding this Friedman’s test, the computed pvalue has been 0.000141.To evaluate the statistical significance of the better performance of DFA, the Holm’s posthoc test has been conducted using DFA as control algorithm. The unadjusted and adjusted pvalues obtained through the application of Holm’s posthoc procedure can be seen in Table 7. Analyzing this data, and taking into account that all the pvalues are lower than 0.05, it can be concluded that DFA is significantly better than ESA and EA at a 95% confidence level.
Algorithm  Average Ranking 

DFA  1.2 
ESA  2.0667 
EA  2.7333 
Algorithm  Unadjusted  Adjusted 

ESA  0.017622  0.017622 
EA  0.000027  0.000054 
6 Conclusions and Further Work
In this work, a new version of the newspaper delivery problem with recycling policy has been tackled. This problem has been modelled as a rich vehicle routing problem, specifically, as an asymmetric and clustered vehicle routing problem with simultaneous pickup and deliveries, variable costs and forbidden paths. It is the first time that a problem like this is presented in the literature, for this reason, a benchmark composed by 15 instances has been also developed, using realworld geographical locations. To deal with such a complex problem, a discrete firefly algorithm has been developed. This application can be considered as the first application of a firefly algorithm to any rich vehicle routing problem. To prove that the proposed DFA is a promising technique, its performance has been compared with those obtained by two other wellknown techniques: the evolutionary algorithm, and the evolutionary simulated annealing. The DFA has outperformed these two classic metaheuristics.
As for future work, it is intended to extend the application of the DFA to other complex realworld situations, related to transportation and logistics. In addition, comparison with more algorithms and techniques will be carried out. Various improvements will be investigate so as to see if the results shown in this work for the ACVRPSPDVCFP can be improved. Besides this, it is intended to conduct a thorough study on the parameterization of the DFA (analyzing, for example, the time complexity). This study has not been done in this work because the main objective is to present the NDSRP model, and to demonstrate that the DFA shows an adequate performance applied to an RVRP. Finally, it would be useful to perform ACVRPSPDVCFP instances using random changes in the asymmetric travel costs.
Acknowledgement
This project was supported by the European Union’s Horizon 2020 research and innovation programme through the TIMON: Enhanced real time services for optimized multimodal mobility relying on cooperative networks and open data project (636220). As well as by the projects TEC201345585C22R from the Spanish Ministry of Economy and Competitiveness, and PC201371A from the Basque Government.
References
 [1] Golden, B.L., Wasil, E.A.: Or practice—computerized vehicle routing in the soft drink industry. Operations research 35(1) (1987) 6–17
 [2] Vidal, T., Crainic, T.G., Gendreau, M., Prins, C.: Heuristics for multiattribute vehicle routing problems: A survey and synthesis. European Journal of Operational Research 231(1) (2013) 1–21
 [3] Lahyani, R., Khemakhem, M., Semet, F.: Rich vehicle routing problems: From a taxonomy to a definition. European Journal of Operational Research 241(1) (2015) 1–14
 [4] Toth, P., Vigo, D.: The vehicle routing problem. Society for Industrial and Applied Mathematics (2001)

[5]
de Armas, J., MeliánBatista, B., MorenoPérez, J.A., Brito, J.:
Gvns for a realworld rich vehicle routing problem with time windows.
Engineering Applications of Artificial Intelligence
42 (2015) 45–56  [6] Amorim, P., Parragh, S.N., Sperandio, F., AlmadaLobo, B.: A rich vehicle routing problem dealing with perishable food: a case study. Top 22(2) (2014) 489–508
 [7] Lahyani, R., Coelho, L.C., Khemakhem, M., Laporte, G., Semet, F.: A multicompartment vehicle routing problem arising in the collection of olive oil in tunisia. Omega 51 (2015) 1–10
 [8] CaceresCruz, J., Arias, P., Guimarans, D., Riera, D., Juan, A.A.: Rich vehicle routing problem: Survey. ACM Computing Surveys (CSUR) 47(2) (2014) 32
 [9] Glover, F.: Tabu search, part i. ORSA Journal on computing 1(3) (1989) 190–206
 [10] Kirkpatrick, S., Gellat, C., Vecchi, M.: Optimization by simmulated annealing. science 220(4598) (1983) 671–680

[11]
Goldberg, D.:
Genetic algorithms in search, optimization, and machine learning.
AddisonWesley Professional (1989)  [12] De Jong, K.: Analysis of the behavior of a class of genetic adaptive systems. PhD thesis, University of Michigan, Michigan, USA (1975)
 [13] Dorigo, M., Blum, C.: Ant colony optimization theory: A survey. Theoretical computer science 344(2) (2005) 243–278

[14]
Kennedy, J., Eberhart, R., et al.:
Particle swarm optimization.
In: Proceedings of IEEE international conference on neural networks. Volume 4., Perth, Australia (1995) 1942–1948
 [15] Rodriguez, A., Gutierrez, A., Rivera, L., Ramirez, L.: Rwa: Comparison of genetic algorithms and simulated annealing in dynamic traffic. In: Advanced Computer and Communication Engineering Technology. Springer (2015) 3–14
 [16] Cao, B., Glover, F., Rego, C.: A tabu search algorithm for cohesive clustering problems. Journal of Heuristics (2015) 1–21
 [17] İnkaya, T., Kayalıgil, S., Özdemirel, N.E.: Ant colony optimization based clustering methodology. Applied Soft Computing 28 (2015) 301–311

[18]
AtashpazGargari, E., Lucas, C.:
Imperialist competitive algorithm: an algorithm for optimization
inspired by imperialistic competition.
In: IEEE Congress on Evolutionary Computation. (2007) 4661–4667
 [19] Yang, X.S.: A new metaheuristic batinspired algorithm. In: Nature inspired cooperative strategies for optimization. Springer (2010) 65–74
 [20] Geem, Z.W., Kim, J.H., Loganathan, G.: A new heuristic optimization algorithm: harmony search. Simulation 76(2) (2001) 60–68
 [21] Moscato, P., Cotta, C.: A gentle introduction to memetic algorithms. In: Handbook of metaheuristics. Springer (2003) 105–144
 [22] Nalepa, J., Blocho, M.: Cooperation in the parallel memetic algorithm. International Journal of Parallel Programming (2014) 1–28
 [23] Vidal, T., Crainic, T.G., Gendreau, M., Prins, C.: A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with timewindows. Computers & Operations Research 40(1) (2013) 475–489
 [24] Vidal, T., Crainic, T.G., Gendreau, M., Lahrichi, N., Rei, W.: A hybrid genetic algorithm for multidepot and periodic vehicle routing problems. Operations Research 60(3) (2012) 611–624
 [25] Qi, Y., Hou, Z., Li, H., Huang, J., Li, X.: A decomposition based memetic algorithm for multiobjective vehicle routing problem with time windows. Computers & Operations Research 62 (2015) 61–77
 [26] Bortfeldt, A., Hahn, T., Männel, D., Mönch, L.: Hybrid algorithms for the vehicle routing problem with clustered backhauls and 3d loading constraints. European Journal of Operational Research 243(1) (2015) 82–96
 [27] Zhang, Z., Che, O., Cheang, B., Lim, A., Qin, H.: A memetic algorithm for the multiperiod vehicle routing problem with profit. European Journal of Operational Research 229(3) (2013) 573–584
 [28] Nagata, Y., Bräysy, O., Dullaert, W.: A penaltybased edge assembly memetic algorithm for the vehicle routing problem with time windows. Computers & Operations Research 37(4) (2010) 724–737
 [29] Nalepa, J., Blocho, M.: Adaptive memetic algorithm for minimizing distance in the vehicle routing problem with time windows. Soft Computing (2015) 1–19
 [30] Marinakis, Y., Marinaki, M., Spanou, P.: A memetic differential evolution algorithm for the vehicle routing problem with stochastic demands. In: Adaptation and Hybridization in Computational Intelligence. Springer (2015) 185–204
 [31] Yang, X.S.: Natureinspired metaheuristic algorithms. Luniver press, UK (2008)
 [32] Fister, I., Yang, X.S., Fister, D., Fister Jr, I.: Firefly algorithm: A brief review of the expanding literature. In: Cuckoo Search and Firefly Algorithm. Springer (2014) 347–360
 [33] Fister, I., Fister Jr, I., Yang, X.S., Brest, J.: A comprehensive review of firefly algorithms. Swarm and Evolutionary Computation (2013)
 [34] Ma, Y., Zhao, Y., Wu, L., He, Y., Yang, X.S.: Navigability analysis of magnetic map with projecting pursuitbased selection method by using firefly algorithm. Neurocomputing (2015)
 [35] Liang, R.H., Wang, J.C., Chen, Y.T., Tseng, W.T.: An enhanced firefly algorithm to multiobjective optimal active/reactive power dispatch with uncertainties consideration. International Journal of Electrical Power & Energy Systems 64 (2015) 1088–1097
 [36] Zouache, D., Nouioua, F., Moussaoui, A.: Quantuminspired firefly algorithm with particle swarm optimization for discrete optimization problems. Soft Computing (2015) 1–19
 [37] Yip, P.P., Pao, Y.H.: Combinatorial optimization with use of guided evolutionary simulated annealing. IEEE Transactions on Neural Networks 6(2) (1995) 290–295
 [38] Ree, S., Yoon, B.S.: A twostage heuristic approach for the newspaper delivery problem. Computers & industrial engineering 30(3) (1996) 501–509
 [39] Archetti, C., Doerner, K.F., Tricoire, F.: A heuristic algorithm for the free newspaper delivery problem. European Journal of Operational Research 230(2) (2013) 245–257
 [40] Campbell, A., Clarke, L., Kleywegt, A., Savelsbergh, M.: The inventory routing problem. In: Fleet management and logistics. Springer (1998) 95–113
 [41] Kallehauge, B., Larsen, J., Madsen, O.B., Solomon, M.M.: Vehicle routing problem with time windows. Springer (2005)
 [42] Boonkleaw, A., Suthikarnnarunai, N., Srinon, R.: Strategic planning and vehicle routing algorithm for newspaper delivery problem: Case study of morning newspaper, bangkok, thailand. In: Proceedings of the World Congress on Engineering and Computer Science. Volume 2. (2009) 1067–1071
 [43] Hurter, A.P., Van Buer, M.G.: The newspaper production/distribution problem. Journal of Business Logistics 17 (1996) 85–107
 [44] Van Buer, M.G., Woodruff, D.L., Olson, R.T.: Solving the medium newspaper production/distribution problem. European Journal of Operational Research 115(2) (1999) 237–253
 [45] Laporte, G., Mercure, H., Nobert, Y.: An exact algorithm for the asymmetrical capacitated vehicle routing problem. Networks 16(1) (1986) 33–46
 [46] Toth, P., Vigo, D.: A heuristic algorithm for the symmetric and asymmetric vehicle routing problems with backhauls. European Journal of Operational Research 113(3) (1999) 528–543
 [47] Herrero, R., Rodríguez, A., CáceresCruz, J., Juan, A.A.: Solving vehicle routing problems with asymmetric costs and heterogeneous fleets. International Journal of Advanced Operations Management 6(1) (2014) 58–80
 [48] Chisman, J.A.: The clustered traveling salesman problem. Computers & Operations Research 2(2) (1975) 115–119
 [49] Battarra, M., Erdogan, G., Vigo, D.: Exact algorithms for the clustered vehicle routing problem. Operations Research 62(1) (2014) 58–71
 [50] Vidal, T., Battarra, M., Subramanian, A., Erdogan, G.: Hybrid metaheuristics for the clustered vehicle routing problem. Computers & Operations Research 58(1) (2014) 87–99
 [51] Wang, C., Mu, D., Zhao, F., Sutherland, J.W.: A parallel simulated annealing method for the vehicle routing problem with simultaneous pickup–delivery and time windows. Computers & Industrial Engineering 83 (2015) 111–122
 [52] Li, J., Pardalos, P.M., Sun, H., Pei, J., Zhang, Y.: Iterated local search embedded adaptive neighborhood selection approach for the multidepot vehicle routing problem with simultaneous deliveries and pickups. Expert Systems with Applications 42(7) (2015) 3551–3561
 [53] Haghani, A., Jung, S.: A dynamic vehicle routing problem with timedependent travel times. Computers & operations research 32(11) (2005) 2959–2986
 [54] Villeneuve, D., Desaulniers, G.: The shortest path problem with forbidden paths. European Journal of Operational Research 165(1) (2005) 97–107
 [55] Montané, F.A.T., Galvao, R.D.: A tabu search algorithm for the vehicle routing problem with simultaneous pickup and delivery service. Computers & Operations Research 33(3) (2006) 595–619
 [56] Yang, X.S.: Firefly algorithms for multimodal optimization. In: Stochastic algorithms: foundations and applications. Springer (2009) 169–178
 [57] Jati, G.K., et al.: Evolutionary discrete firefly algorithm for travelling salesman problem. Volume 6943. Springer (2011)
 [58] Zhou, L., Ding, L., Qiang, X.: A multipopulation discrete firefly algorithm to solve tsp. In: BioInspired ComputingTheories and Applications. Springer (2014) 648–653
 [59] Derrac, J., García, S., Molina, D., Herrera, F.: A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm and Evolutionary Computation 1(1) (2011) 3–18
Comments
There are no comments yet.