Genetic Algorithm for Grocery Delivery Systems

(1)

Genetic Algorithm for

Grocery Delivery Systems

Vincent Velthuis 6381448 / 10120750

Bachelor thesis Credits: 18 EC

Bachelor Opleiding Kunstmatige Intelligentie University of Amsterdam Faculty of Science Science Park 904 1098 XH Amsterdam Supervisor Dr. S. van Splunter Informatics Institute Faculty of Science University of Amsterdam Science Park 904 1098 XH Amsterdam June 29th, 2017

(2)

Abstract

Grocery delivery is a service provided by retail companies. The optimisation of delivery routes is an important part of the delivery process. An optimal route minimises travel time and maximises efficiency of resource usage. The aims of this study were to assess what improvement in travel time could be made by Genetic Algorithms (GA) compared with random delivery route solutions, and to assess how many generations were required to obtain optimal results over computation time. Genetic algorithms apply evolutionary principles of selection, genetic recombination and mutation for general problem solving purposes. In this study, random delivery routes were constructed for a commercial dataset and evolved for 1000 to 10000 generations, using three different GA designs: uniform crossover, two-point crossover, and 50% truncation selection. For all GA, the travel time improved compared with the random delivery routes. Best results were obtained using the two-point crossover, with an average reduction of travel time by 23.7 minutes (7.5%). Results were not improved by increasing the number of generations beyond 1000, and show that GA can make a feasible delivery route solution in a fraction of the computation time of brute-force calculations.

Acknowledgements

The author would like to thank J. Prins for providing the dataset, drs. H. Kaufmann for planning support, Y.M.A.A. Galama for programming advice and coffee breaks and A. Kok for proofreading, insightful discussions and mo-tivational support. Special thanks goes out to dr. S. van Splunter for regular feedback, quick responses and overall guidance.

(3)

3 Method and Approach 5 3.1 VRP parameters . . . 5 3.2 Initial population . . . 5 3.3 Reproductive process . . . 6 3.4 Dataset . . . 8 3.4.1 Cleanup . . . 8 3.5 Implementation . . . 9 4 Results 9 4.1 Experiment setup . . . 9 4.2 Experiment results . . . 9 5 Conclusion 11 6 Discussion 11 6.1 Crossover . . . 12 6.2 Benchmark . . . 13 6.3 Processing time . . . 13

(4)

1 Introduction

Grocery delivery is an additional service provided by retail companies, where a customer can make an online selection of products which are to be delivered to them within a specified time slot. The retailer collects the products at a depot and delivers them to the customer at a selected location and time. Each time slot can contain multiple customers and therefore a route connecting all customers and the depot must be constructed. Because the retailer has a finite amount of time and resources (e.g., vehicles, staff), the optimisation of these routes is an important part of the delivery process.

For this project, a route optimisation method is proposed for use in grocery delivery. A dataset was provided by a grocery delivery service, allowing the investigation of optimisation of their own routing approach. The provided data consists of historical data points of previous deliveries. These provided routes are used in this thesis to evaluate delivery routes or solutions generated by genetic algorithms. The provided routes are further referenced as benchmark routes. The dataset is described in more detail in subsection 3.4.

An optimal route minimises the travel time and maximises efficiency of resource usage, i.e. the number of used vehicles and how many orders they contain. An Artificial Intelligence approach is perfectly suited for these min/max optimi-sation criteria. In this thesis the heuristic frameworks of the Vehicle Routing Problem (VRP) and Genetic Algorithm (GA) are used, further described in section 2.

The research questions for this project are: What improvement in travel time can be made by Genetic Algorithms compared with random route solutions with efficiency constraints? And how many generations are required to obtain optimal results over computation time?

2 Theoretical Framework

2.1 Travelling Salesman Problem

To understand the Vehicle Routing Problem, the Travelling Salesman Problem (TSP) [Rego et al., 2011] needs to be defined first. In a TSP, the goal is to find a route solution (or tour) that minimises the travel distance it takes an agent to visit all vertices (customers) in a graph and return to its starting point (depot). The solution to the TSP is NP-hard, with n! solutions for a problem with n vertices. Heuristic approximation of the optimal solution is therefore better suited to tackle these problems than brute-force calculations. This thesis uses the upper bound of the Nearest Neighbour Algorithm (NNA) as its heuristic approximation of the tour time. The NNA creates a tour by travelling to the nearest unvisited vertex from the current vertex, until all vertices have been visited. From the last vertex, the agent returns to the start vertex.

However, the NNA does not necessarily result in the optimal or shortest route. For example, a drawback of the NNA is that there is no control over the choice of the final edge; it may be proportionally longer than expected. Normally, this

(5)

drawback is circumvented by exploring the graph from different starting vertices, but since the starting vertex is set at the depot, this method is not applicable. Still, if the solution from the NNA is assumed to be an approximation of the optimal route, the decrease in complexity is a big improvement compared with a randomly generated route.

2.2 Vehicle routing problem

A Vehicle Routing problem (VRP) [Laporte, 1992] is a generalised version of the TSP, capable of handling solutions with multiple vehicles, under a set of side constraints. The problem is to find a set of least-cost vehicle routes such that:

• each vertex in the graph is served exactly once by exactly one vehicle • each vehicle route starts and ends at the depot

• all side constraints are satisfied

The aforementioned side constraints for the VRP can include for instance a maximum number of orders per vehicle, or a maximum route duration. These side constraints are specific to each unique instance of the VRP. For this project, the side constraints are explained in subsection 3.1.

2.3 Genetic Algorithms

A VRP can have multiple solutions, and to investigate the fitness potential of the benchmark routing solution, this project uses the Genetic Algorithms heuristics [Baker and Ayechew, 2003] to explore possible solutions to the VRP. The Genetic Algorithm (GA) is a general method of problem solving based on the evolutionary principles of selection, genetic recombination and mutation. It starts with an initial population of chromosomes, each representing a possi-ble distribution of the orders over the different vehicles for the given propossi-blem. For each chromosome, a fitness function quantifies it’s efficiency (or fitness) at handling the problem. Over the course of several generations, the population evolves by generating offspring based on recombination of the genetic code of a selection of the fittest members. Offspring are inserted into the population only if their fitness is higher than that of any population member, creating a higher mean fitness over each iteration. Further elaboration on the specific parameters for the GA will be discussed in section 3.

2.4 Heuristics for the VRP

Other possible heuristics that were considered for this project are Ant Colony Optimization[Dorigo and St¨utzle, 2010][Bell and McMullen, 2004], Greedy Ran-domized Search Procedure (GRASP) [Baker and Carreto, 2003], Simulated An-nealing [Zeng et al., 2005], Tabu Search [Gendreau et al., 1994] and Variable Neighborhood Search [Bula et al., 2017]. These heuristics were ommited for the reasons as discussed below.

(6)

Ant Colony Optimization was proven to be less effective when measured against other state-of-the-art algorithms for the TSP [Dorigo and St¨utzle, 2010, p. 3], and therefore for the VRP. GRASP needs a functioning user interface and re-lies on user interaction which is too time consuming to set up for this project. Simulated Annealing and Tabu Search are described as improvement heuris-tics, constructing new route solutions from a given starting route. These were omitted because of their strong reliance on a near-optimal route, which is not available when vertices are highly variable, as is the case in the grocery deliv-ery service. Variable Neighborhood Search was discarded because the literature mostly described computation on VRPs that were much larger than ours (20000 vertices against 18) and it was unknown if the complexity of that heuristic was needed for this project.

Genetic Algorithms are suitable for solving small VRPs [Baker and Ayechew, 2003, p. 797]. They are not reliant on a near-optimal reference route, because they construct their own random population of routes. Moreover, GAs can be de-veloped for this problem within the time constraints of the project. These rea-sons, combined with the extensive documentation of the example application in [Baker and Ayechew, 2003], make the GA an suitable choice for this project.

3 Method and Approach

3.1 VRP parameters

To describe the specific problem adressed in this thesis formally, we define it as a VRP. Let G = (V, A) be a graph with A the arc set and the vertex set V = {0, ..., n}, where 0 is the depot and the other vertices are the customers. The non-negative distance matrix D = (dij) can be interpreted as the travel

time (in seconds) over arc (i, j), i 6= j. The demand qi is 1 for all vertices, and

the sum of demands Q on any vehicle should not exceed the vehicle capacity Qmax = 6. For each delivery, a processing time p = 840(sec) is added to the

travel time, based on a prior research by the delivery company. Each vehicle must deliver their load within a time limit of 2 hours (7200 seconds) to guarantee a freshness of goods. The fleet of vehicles is homogeneous, since the capacity and fixed cost f are the same for all vehicles.

The graph was constructed using data from Google MapsAPI, constructing a distance matrix between all customers and the depot. This matrix contains the travel time in seconds between all locations. We represent the distance in time, because our specific problem uses a time constraint for each vehicle and calculation is simplified by using this quantity. If a vehicle should exceed their capacity limit, a penalty is given in the form of an increased unfitness (see subsection 3.3). This allows for solutions with higher vehicle load to exist, while at the same time discouraging overloading the vehicles.

3.2 Initial population

The GA maintains a population of possible solutions, which will be further referred to as chromosomes. The structure of the chromosomes for this project

(7)

should be capable of representing the mapping of each delivery location to a particular vehicle. If the vehicles are numbered as v ∈ [1, 3], a chromosome of length G looks as in table Table 1. Each gene value represents the vehicle to which the customer at that index is mapped, e.g. customer 1 is mapped to vehicle 1, customer 2 to vehicle 3. Index 0 is not represented in the chromosome, because it refers to the depot which doesn’t have to be mapped to any vehicle.

Table 1: Chromosome representation. index 1 2 3 4 ... G-1 G gene 1 3 2 1 ... 2 1

The initial population is created by seeding N chromosomes of length G. (N = 100, G ∈ [1, 18]) Each chromosome is seeded with random integers, without a predefined distribution. Double occurring chromosomes are reseeded, ensuring unique chromosomes in the population. In all VRPs, the maximum number of available vehicles at each depot is represented in the chromosomes, although the ideal number of vehicles is dependent of the size of G and Qmax. For instance,

a VRP with 11 orders and a maximum of 3 vehicles will in most cases suffice with just 2 vehicles. Using more vehicles than necessary is penalised by adding a loading time l = 900(sec) to the total travel time for each vehicle present in the chromosome. This will help the GA to converge to solutions using minimal vehicles.

3.3 Reproductive process

Each chromosome in the population is evaluated by the fitness function, re-turning a fitness and unfitness value. The fitness value represents the objective function value, in the case of the VRP the total travel time. The total travel time is an addition of the vehicle travel time for all used vehicles. The individual vehicle travel times are calculated as an addition of the individual travel times, drop allowances and loading times for each order in the vehicle V , as in:

Vt=

X

i,j

(dij) + G ∗ p + l

The unfitness value represents the level of constraint violations. If a vehicle exceeds the maximum vehicle capacity or the vehicle travel time limit, it receives a penalty relative to the severity of the violation. The unfitness is expressed as a proportion of any excess load or of travel time over the maximum value. In an ideal solution, the unfitness value should be 0; there is no theoretic maximum to the unfitness.

Parent chromosomes are selected from the population by binary tournament. Two random chromosomes are selected, and the individual with the highest fitness is selected as the first parent. The process is repeated to determine the second parent.

The two parents selected for reproduction are recombined into offspring by means of either uniform or 2-point crossover. Uniform crossover creates off-spring by randomly selecting genes from each parent for each index position.

(8)

With 2-point crossover, two cuts are made at random points in the chromosomes of the parents and the intermediate genes are switched. Both crossover methods are shown in Figure 1. After applying either both crossover method, offspring are mutated by switching the values of two random genes. These mutations are performed to increase the genetic diversity in the overall population, making it easier for the GA to avoid local optima.

Figure 1: Crossover Methods, adapted from [Mehboob et al., 2016]. Each generation yields two possible offspring genomes, which are disregarded if their complete genetic code is already present in the population. If the offspring has a novel gene combination, it is evaluated for insertion into the population by a replacement scheme. For each offspring, the population is split into four subsets with respect to the child’s fitness and unfitness value. Representing the offspring by c and the fitness and unfitness value of each population member (p ∈ P ) as f (p) and u(p) respectively, the subsets are defined as

S1= {p ∈ P : f (p) ≥ f (c), u(p) ≥ u(c)},

S2= {p ∈ P : f (p) < f (c), u(p) ≥ u(c)},

S3= {p ∈ P : f (p) ≥ f (c), u(p) < u(c)},

S4= {p ∈ P : f (p) < f (c), u(p) < u(c)}

The population member to be replaced is chosen from the first non-empty set in the order S1, S2, S3. If all these sets are empty, the offspring does not enter

the population. In the chosen subset, the member with the worst unfitness is selected for replacement, with ties resolved by worst fitness. This gives the GA

Table 2: GA design choices. Design choice Method used (exp) Node sorting No sorting

Population seeding Random Integers

Parent selection Binary Tournament (1,3), 50% Truncation (2) Crossover Uniform Crossover(1,2), 2-point Crossover(3) Mutation 2-node switch

(9)

a preference to explore solutions that are within the constraints, even though an improvement in unfitness could be coupled with a reduced fitness. Allowing the fitness to decrease or unfitness to increase makes it easier for the GA to escape from local optima.

In summary, the design for the GA as discussed above is given in Table 2.

3.4 Dataset

A set of historical data was provided for this project by a grocery delivery service, wanting to evaluate their own routing software. The company remains undisclosed, and the dataset ensures the anonimity of their customers. Each customer is represented as a postal code, which gives an indication of their location but no further details regarding their personal information. There are 9 different depots from which orders are delivered over 4 time slots on each day. Each depot is represented as a code instead of a location. At each depot, 2 or 3 vehicles are accessible, corresponding with a maximum time slot capacity of 12 or 18 orders respectively. Each time block at each depot is regarded as a single VRP, making the computational effort relatively small for each problem.

3.4.1 Cleanup

he original .xls dataset consists of 6878 data points, each containing a date, depot number, time slot, order number and postal code. A Dutch postal code consists of four digits and two letters, e.g., 2017AM. The first two digits indicate a region, the last two digits indicate a district and the two letters indicate a neighbourhood and street. Each four-digit combination is serviced by exactly one depot. The dataset had to be processed to filter out several errors. Of the 6878 data points, 9 were discarded due to absence of letters in the postal code. 376 data points which were indexed at the wrong depot were allocated to the correct depot. An additional dataset was created where those data points were instead deleted. The difference in Figure 2 between the complete dataset (blue)

(10)

and correct dataset (red) bars shows the amount of orders that were wrongfully indexed at each depot. The green bar represents the dataset after allocation, which was used for the experiments of this project.

3.5 Implementation

The dataset, initially in.xls, was converted to .csv for easier processing in Python. A .csv (Comma Seperated Value) can be loaded into python as a numpy array with numpy’s genfromtxt() function. To allow seperate testing and training, the data was divided into 9 smaller sets, each containing the data of one depot. To evalutate fitness of routes in each time block, each depot-specific numpy array was processed by a conversion program convert.py. This programs takes a .npy file containing all orders for each depot and exports it as a .p (pickle dump) file containing order locations and the corresponding distance matrix between all locations. The name of the file contains information about the depot and timeblock, e.g. ’44422 2-27 09-11.p’. The output progress is slow due to query limitations of 2500 requests per developer key per day to the Google Maps Directions API, which is used to construct the distance matrix. The progress was sped up by using three seperate API developer keys instead of one, increasing the daily query limit from 2500 to 7500. Of the provided postal codes, 83 were unavailable in Google Maps. The associated postal codes have been documented and changed to the nearest existing postal code.

4 Results

4.1 Experiment setup

The impact of using the GA was assessed by comparing fitness of routes cal-culated by the GA routes with the best routes from the (random) initial pop-ulation. This evaluation was performed by running 5 iterations of the GA per time block and denoting the best known member of the initial population and the population after a given number of generations (1000, 2500, 5000 and 10000 generations). Time blocks with 0-6 orders were omitted, because they can be calculated without a GA by putting all orders in one vehicle and performing a NNA over the data matrix for that vehicle. For each run, the absolute and rela-tive improvement rates were recorded from the average fitness of the 5 random routes to those of the GA routes.

4.2 Experiment results

Three experiments were performed, each experiment differing in one design choice from the original design as described in section 3. The results for all experiments can be assessed in Table 3.

For experiment 1, the original GA design as described in section 3 was used, with a binary tournament selection and uniform crossover method. On average,

(11)

(a) Relative improvement (b) Absolute improvement

Figure 3: Improvement rates for Exp1 (5000 generations).

fitness improved by 3.48 % or 10.51 minutes per time block. The best improve-ment was docuimprove-mented with 1000 generations, showing an average improveimprove-ment of 3.57 % or 10.79 minutes. For these settings, larger generation sizes do not nec-essarily yield positive results. This is reflected in the decrease in improvement for 10000 generations.

The relative improvement of the GA route over the best random solution varies with the number of orders, with larger order sizes yielding lower rates of im-provement (Figure 3a).

The mean improvement rate remains positive for most time blocks, although there are some negative values for all runs. The absolute fitness improvement seems constant across order size (Figure 3b).

Table 3: Relative and absolute improvement rates for the experiments.

Number of Exp 1 Exp 2 Exp 3

generations % min SE % min SE % min SE 1000 3.57 10.79 0.24 3.18 9.80 0.23 7.35 23.20 0.43 2500 3.42 10.25 0.25 2.99 9.16 0.23 7.55 23.86 0.44 5000 3.50 10.61 0.26 3.10 9.43 0.22 7.56 23.95 0.46 10000 3.43 10.37 0.25 3.13 9.56 0.23 7.55 23.90 0.45 Average 3.48 10.51 3.10 9.49 7.50 23.73

For the second experiment, the original design was used, with a change in the selection method. The binary tournament was replaced with a 50% truncation method. For this selection method, the population is split between the 50% fittest population members and the other 50% unfittest members. A random parent is selected from the subset of fittest population members.

For every number of generations, the relative and absolute improvement was below that those of experiment 1. The mean improvement rate of all time

(12)

blocks was 3.10 % or 9.49 minutes.

The third experiment was a continuation of the design of experiment 1, with a 2-point crossover instead of uniform crossover. This experiment yielded better average results than the first two experiments, with an average improvement rate of 7.50% or 23.73 minutes.

5 Conclusion

The results in section 4 show that experiment 3 yields the best improvement rate overall. From this, it can be concluded that for this dataset, the 2-point crossover method is preferred over uniform crossover. The average improvement rate of the 2-point crossover is almost twice as large as the uniform crossover, in both relative as absolute measurements.

Between experiment 1 and experiment 2, the binary tournament selection method can be assessed as a better alternative than a 50% truncation.

From the spread around the average improvement (X ± SE) subsequent gener-ations, it can be concluded that a GA with 2,500 or more generations yields no significant improvement over 1,000 generations. In every experiment, 1,000 gen-erations are sufficient in obtaining a decent approximation of a well-constructed route.

6 Discussion

To evaluate the applicability of the GA, the stakeholders of this system need to be examined first. On the one hand, there are the academic stakeholders, interested not so much as to if a system works, but how it works and why it is better than any given benchmark method. On the other hand, there are the business stakeholders. For them, the most important output is a working system with a distinct added value over the original routing methods. That added value can manifest as any of the side parameters of a VRP or as the amount of calculation time needed.

From the academic point of view, it can be argued that the complexity of the delivery problems presented is too low to support the use a GA. The small-est problem addressed by Baker & Ayechew had 5 vehicles and 50 customers [Baker and Ayechew, 2003, p .797]. In comparison, this project uses only 3 ve-hicles and a maximum of 18 customers. Is it justified to use a GA for these relatively small problems?

Answers to these questions about applicability can be found by taking the inter-est of the business stakeholders in mind. These stakeholders do not necessarily want the best route, they can make good use of an approximation that is cal-culated much more quickly. Time is a valuable resource, even more so because of the needed flexibility in which a route needs to be constructed when orders are added to a time block.

(13)

good solution short, independent of the number of orders; it converges to a proper solution in the same number of generations regardless of the problem size. Within this project, the largest GA problem has 18 orders to be distributed over 3 busses. This distribution can be done in 18₆ ∗ 12

6 = 1.72 ∗ 10

7_{ways, and}

that is only if the solution uses exactly 6 orders in the first 2 busses. The routes within the individual busses are not taken into account in this calculation, the actual problem size is even bigger. For this same problem, a GA can make an approximation of a good solution in a fraction of the computation time. A feasible solution can be found with 100 random routes in the initial population, and a run time of 1,000 generations. This structuring of the data makes it so that in an 10,000 times smaller computation time, a feasible solution estimate can be found.

Should a company choose to expand to delivery blocks with a larger amount of customers, the same principles of the smaller GA can be integrated into the larger system. When looking at a possible larger problem, for instance 24 customers over 4 vehicles, a brute-force method has to investigate 2.3 ∗ 109data points, while the GA computation time remains of a similar size.

6.1 Crossover

Looking at the results, the GA with a 2-point crossover method yielded better results than a uniform crossover. It was originally hypothesised that a uniform crossover should converge to an ideal solution faster than a 2-point crossover.

Figure 4: Order locations and polar an-gles around depot P.

Since the GA does not sort the or-der locations in any way, the distance between two orders in the chromo-some contains no information about the distance between order locations in reality. Therefore, the sequence of genes should in theory be incon-sequential to the fitness of the chro-mosome. The results of experiment 3 show however that this is not the case, yielding a higher improvement rate for the 2-point crossover. From this it can be concluded that the structure of the chromosome se-quences is more relevant than initially thought. For future research, it could be investigated if the sorting of or-der locations according to a common principle increases the improvement rate. When using postal codes as a location index this sorting can be done with a standard character sort, because of the inherent clustering of locations built in to the postal code

(14)

be-cause 1521AL is by definition located in the same neighbourhood and street. A drawback for this method is that neighbouring postal areas don’t necessar-ily have neighbouring numbers, e.g., 1521 is bordered by 1561 and 1531, but not 1522. An alternative method would be to calculate the polar angle of each order relative to the depot and sort the orders according to those angles [Baker and Ayechew, 2003]. This clusters orders together based on their direc-tions from the depot, but does not take their relative distance into account. This method is preferably used for dataset with a central depot with a circular spread of orders, but may not be useful when this is not the case.

For instance, in Figure 4 the polar angle of orders C and K are relative to P are close together, although the orders themselves are not.

6.2 Benchmark

Originally, the results of the GA were scheduled to be compared with the original routes used by the company which provided the dataset. These routes are further referenced as benchmark routes.

Because the original routing method is a ’black box’, it can not be answered to its full extent why the results of GA are better of worse than a benchmark method. The results of this thesis show that the GA is not an improvement over the benchmark method, but it gives no information about the assumptions made for the benchmark method and how similar assumptions would influence the results for the GA. From personal correspondence with some of the drivers, the suggestion was made that the benchmark system was often structured in such a way as to avoid traffic jams which occurred on regular intervals at specific locations. Since our data matrices are constructed at one specific time point, it is unknown how applicable it remains when used on different times of day. Further research into these questions needs to be done, before any definitive statement about the GA’s added value. A potential solution for this would be to construct the data matrices at their corresponding time slots, which would take the average road usage into account for that particular time. Another solution would be to construct the data matrices at different times of day, regardless of the original time block, and compare the improvement rates between those systems.

6.3 Processing time

After the experiments for this project were concluded, a visualisation of a so-lution was prepared for a presentation. This image (as shown in Figure 5) brought to light a flaw in the implementation. On the right a GA solution to one particular problem is shown, on the left the benchmark solution for the same problem. The GA solution had a better fitness than that of the benchmark route, although the benchmark was clearly better structured.

After debugging, it was found that the processing time for each order was not set at 840 seconds as documented in this thesis, but at 480 seconds. This meant that for each order, the fitness function failed to incorporate 6 minutes of loading time, adding op to a maximum of 36 per vehicle. This error had

(15)

the effect that the evolutionary pressure of the unfitness was not sufficient. A chromosome could be assessed to be within the vehicle time constraints, while actually it was unfit for the time constraints. This occurs most with time blocks around the maximum capacity, e.g., 17 or 18 orders. With the error in the GA, there was not enough selective pressure from the unfitness value to evolve to a well-structured solution.

(a) Benchmark (b) GA solution (pi= 480)

(c) GA solution (pi= 840)

Figure 5: Example of order distribution among vehicles (G = 18). After the error was fixed, the resulting route for this problem were better struc-tured than in the GA with the incorrect processing time, as shown in 5c. However, within the time constraints of the project, there was no time to rerun the entirety of the experiments with this parameter change.

References

[Baker and Ayechew, 2003] Baker, B. M. and Ayechew, M. (2003). A genetic algorithm for the vehicle routing problem. Computers & Operations Research, 30(5):787 – 800.

[Baker and Carreto, 2003] Baker, B. M. and Carreto, C. A. (2003). A visual interactive approach to vehicle routing. Computers & Operations Research, 30(3):321 – 337.

[Bell and McMullen, 2004] Bell, J. E. and McMullen, P. R. (2004). Ant colony optimization techniques for the vehicle routing problem. Advanced Engineer-ing Informatics, 18(1):41 – 48.

[Bula et al., 2017] Bula, G. A., Prodhon, C., Gonzalez, F. A., Afsar, H. M., and Velasco, N. (2017). Variable neighborhood search to solve the vehicle rout-ing problem for hazardous materials transportation. Journal of Hazardous Materials, 324, Part B:472 – 480.

(16)

[Dorigo and St¨utzle, 2010] Dorigo, M. and St¨utzle, T. (2010). Ant Colony Op-timization: Overview and Recent Advances, pages 227–263. Springer US, Boston, MA.

[Gendreau et al., 1994] Gendreau, M., Hertz, A., and Laporte, G. (1994). A tabu search heuristic for the vehicle routing problem. Management Science, 40(10):1276–1290.

[Laporte, 1992] Laporte, G. (1992). The vehicle routing problem: An overview of exact and approximate algorithms. European Journal of Operational Re-search, 59(3):345 – 358.

[Mehboob et al., 2016] Mehboob, U., Qadir, J., Ali, S., and Vasilakos, A. (2016). Genetic algorithms in wireless networking: techniques, applications, and issues. Soft Computing, 20(6):2467–2501.

[Rego et al., 2011] Rego, C., Gamboa, D., Glover, F., and Osterman, C. (2011). Traveling salesman problem heuristics: Leading methods, implementations and latest advances. European Journal of Operational Research, 211(3):427 – 441.

[Zeng et al., 2005] Zeng, L., Ong, H. L., and Ng, K. M. (2005). An assignment-based local search method for solving vehicle routing problems. Asia-Pacific Journal of Operational Research, 22(01):85–104.

Genetic Algorithm for Grocery Delivery Systems