A Hybrid Genetic Algorithm for a Multi-Depot Vehicle Routing
Problem Variant with Split Delivery and an Inventory Restriction
Nika de Vries
1721712
Master’s thesis Econometrics, Operations Research and Actuarial Studies
Supervisors:
Prof. Dr. K.J. Roodbergen (RUG)
V. Ponsioen (Districon)
Y. Blom (Districon)
Co-assessor
A Hybrid Genetic Algorithm for a Multi-Depot Vehicle Routing
Problem Variant with Split Delivery and an Inventory Restriction
Nika de Vries
April 7, 2014
Abstract
Preface
With this thesis I will conclude my masters in Econometrics, Operations Research and Actuarial Studies. My research period started with an internship at Districon, however due to personal circumstances I finished the last part consisting of computational experiments at home. Within this thesis the reader will find a result of six months workload on the subject of home delivery distribution.
Although due to circumstances I could not utilize the internship opportunity as much as I had expected and would have liked, it was a very valuable experience. Im am very glad with the practical knowledge I was able to gain. Without the support and involvement during my recovery period I wouldn’t have been able to finish the project successfully. I would like to thank all my former Districon colleagues for this. My special thanks goes to Yvonne Blom and Victor Ponsioen, for the understanding and helping me to make the right decisions with regards to my thesis as well as on a personal level.
Of course, my thanks also goes out to Professor Kees Jan Roodbergen. In his supervision I appreciate the space and flexibility he gave me to work independently and tuned to my level of health. Moreover, I am very grateful for his mathematical expertise. This definitely improved the quality of my research. I also would like to thank Professor Ruud Teunter for his feedback, which inspired me to make some final improvements.
Also, I would like to thank my boyfriend, friends an family for their unconditional and great support.
Contents
1 Introduction 7 1.1 Situation . . . 7 1.2 Motivation . . . 7 1.3 Research goal . . . 9 1.4 Thesis outline . . . 10 2 Literature Review 11 3 Problem Formulation 14 3.1 Assumptions . . . 14 3.2 Methodology . . . 15 3.3 Notation . . . 17 3.4 Mathematical model . . . 18 4 Genetic Algorithm 20 4.1 Search Space and Solution Representation . . . 224.2 Algorithm Notation . . . 24
4.3 Generation of Initial Population . . . 25
4.3.1 Assignment . . . 25
4.3.2 Routing . . . 26
4.4 Evaluation . . . 26
4.4.1 Total Actual Costs . . . 27
4.4.2 Total Penalty Costs . . . 27
4.4.3 Diversity Contribution . . . 28
4.4.4 Evaluation Function . . . 29
4.5 Education . . . 29
4.6 Generation of Offspring . . . 30
4.6.1 Stopping Criterion Generation . . . 30
4.6.2 Parent Selection . . . 31
4.6.3 Crossover . . . 31
4.6.4 Routing . . . 35
4.6.5 Evaluation . . . 35
4.6.6 Education . . . 35
4.6.7 Condition Triggering Diversification . . . 35
4.6.8 Stopping Criteria Generation . . . 35
4.7 Stopping Criteria Genetic Algorithm . . . 36
5 Model Extensions 37 5.1 Implementing a Trade-off Cost for Allocating Inventory . . . 37
5.2 Vehicles Available Regarded as Input for the Model . . . 37
5.3 Adapting Penalty Structure to User Preferences . . . 37
5.4 Multiple Central Distribution Centers . . . 38
6 Parameter Specification 39
6.1 Random Generated Data . . . 39
6.2 Fixed Parameters . . . 39
6.3 Final Algorithm Structure . . . 40
6.4 Algorithmic Parameter Tuning . . . 43
7 Results 45 7.1 Results Benchmark against ILP-formulation . . . 45
7.2 Results Algorithmic Benchmark . . . 49
7.3 Effects of a Dynamic Combined Direct and Indirect Distribution System . . . 50
7.4 Effects of Inventory Allocation . . . 51
8 Conclusion 52 9 Limitations and Further Research 53 10 Appendix 58 10.1 Algorithmic Parameter Tuning: Inventory Spreading . . . 58
10.2 Results Algorithmic Benchmark: Unequal Inventory Spreading . . . 58
1
Introduction
1.1
Situation
E-commerce has grown significantly over the past couple of years. Research shows that online shopping turnover in the Netherlands septupled, with an increase from 1.5 billion in 2000 to 9.8 billion in 2012. Causes are both the increase in the number of online customers and moreover the growth of the number of orders per customer. These ongoing developments are causing a gradual shift in the partition of market shares. In the 2005 retail industry, off-line shopping had a market share of 96.3% and online shopping 3.7%. Every year this division changes by on average more than 1%, resulting in a 89.2%- 10.8% balance for 2011 [41].
Many retailers try to benefit from this development by becoming a clicks and mortar company, i.e. not only sell products in physical stores, but also on the Internet. Traditional physical stores introduced webshops and additionally pure e-tailers are opening physical stores. These trends have caused an increase in the number of multi-channel retailers. These retailers commonly use a multi-channel sales strategy, where each channel controls its own operations. Channel preferences between customers vary and moreover individual customers increasingly become multi-channel shoppers, preferring different channels at different moments and different stages of the shopping process. Therefore by using both an online and off-line sales channel a broader range of cus-tomers can be reached [2].
This relatively new sales channel requires a different approach than traditional logistics. One essential factor is the distribution of products to customers [2]. It is increasingly common that (online) purchased products are delivered at home. Research shows 90% of the customers want to receive their orders at home [36]. This contrasts the traditional form of distribution where a product is delivered to a physical store and customers buy it there to take care of the ‘last mile delivery’ themselves. Another important aspect concerning multi-channelling is inventory management. Decisions with respect to inventory pooling, combining inventory for in-store sales and home deliveries have to be made. Moreover, in order to be able to deliver a product to a customer, a source with sufficient inventory must be selected.
1.2
Motivation
Consider the scenario with one central distribution center and multiple intermediate points. The most common method for (online) home deliveries is direct distribution. This means direct delivery from a central distribution center (DC) to customers with vans or trucks. This is a static assignment strategy, since home deliveries are always delivered from the same prespecified source. An alternative manner for home deliveries is indirect distribution. In this strategy products are first distributed with larger vehicles to an intermediate point (for example a store, cross dock center, hub, depot). Then, from this intermediate point home deliveries are routed making use of smaller vehicles [21].
(a) Direct distribution (b) Indirect distribution
Figure 1: Comparison distribution types for home deliveries
Implementing an indirect distribution system might induce lower costs than executing a direct distribution strategy. For example, when considering a retailer with multiple physical stores located throughout a certain area, these stores could be used as intermediate points for indirect distribution. These local stores are replenished via large trucks from a central distribution center on a regular basis, to supply their usual in-store demand. When using excess truck capacity to distribute products for home delivery to the stores (i.e. the intermediate points), the products are closer to the customers, without (significant) additional cost. Consequently, shorter delivery routes can be achieved and delivery speed can be improved, leading to cost savings and service improvement. Moreover, an indirect distribution system might result in cost savings because of economies of scale [2].
Combining the two approaches (direct and indirect distribution) is also possible. However, even if retailers already use the more advanced combined strategy, this is generally done using a static assignment strategy. For each customer location it is prespecified from which source possible demand will be supplied. This in contrast to dynamic assignment strategies, for which home delivery responsibilities are determined in real time for each incoming demand. Research shows that dynamic assignment strategies can significantly reduce total costs [29]. Therefore different structures for home delivery will be analysed in this research; a dynamic strategy combining direct and indirect distribution versus i) a static assignment strategy considering all possible sources and ii) the practically common static assignment strategy ’direct distribution from the DC’.
in-ventory levels are an important decision factor for the distribution process. A static distribution strategy does not consider the current inventory positions of the sources. This increases the need to backorder or the number of lost sales. A dynamic distribution strategy assigns orders in a smarter way, satisfying the inventory restriction and thereby reducing the number of backorders or lost sales. To be able to implement a real time dynamic distribution system incorporating an inventory restriction is essential.
Acknowledging the importance of retailers having a competitive e-fulfilment strategy, both in terms of costs savings and service levels, this research analyses distribution strategies combined with an inventory aspect. To link the research fields of distribution and inventory management consistent to our problem definition, we focus on the distribution problem approach with an inventory extension aspect. In the literature, the problem of determining how to distribute prod-ucts to a set of customers when there are multiple possible sources, is known as the Multi-Depot Vehicle Routing Problem (MDVRP). The solution to the problem consists of two decisions:
i) the assignment of customers to the sources, ii) the route scheme for each source.
The objective of the MDVRP is to minimize total cost or travel time/ distance, while satisfying certain constraints. For example, all customers should be served, the vehicle capacity cannot be exceeded or the route cannot take longer than a certain duration. The regular MDVRP is a well known problem in scientific literature and as many of its variants is still frequently discussed. The literature related to this problem and literature focussing on a problem combining routing and inventory management over time, the Inventory Routing Problem, will be presented in the Literature Review (Section 2).
The key difference with regards to current literature is that this thesis extends the regular MDVRP with an inventory restriction (at a specific point in time). The current inventory levels of sources are input for the model and cannot be violated to fulfil demand. When inventory is limited allowing split delivery solutions is essential. Therefore, this research takes a broader ap-proach than the regular MDVRP by allowing split delivery solutions between sources. The split delivery option in a multi-depot vehicle routing structure has not received much attention yet, as will be discussed in Section 2. Moreover, in the assumed situation the sources, the distribution center and the intermediate points, are related over time (similar as in [21]). The intermediate points are replenished from the central distribution center, which also can be a direct source for home deliveries. In contrast, the general MDVRP assumes unrelatedness of the sources.
1.3
Research goal
The research goal can be summarized as follows:
1.4
Thesis outline
2
Literature Review
The introduction described the literature concerning the background of the problem setting. Relevant literature on the invasion of e-commerce and the need for sophisticated distribution and inventory models (by retailers) were presented. This section provides a literature review for the MDVRP, a widely studied combinatorial optimization problem. The purpose of this review section is to i) describe and give an overview of currently available techniques for solving the (MD)VRP and ii) specifically focus on genetic algorithms methods to solve the MDVRP. The MDVRP is an extension of the VRP; it considers a system with multiple sources instead of one. The first version of the VRP was presented as a generalization of the likewise exten-sively studied TSP, the Travelling Salesman Problem [12]. For a given set of cities the TSP is concerned with finding the shortest route that visits each city exactly once and ends the route at the originating location. The TSP was generalized to the VRP by adding customer demands and imposing a vehicle load restriction on the fleet. Later on the original VRP was extended with a route duration restriction [22]. Ever since the VRP and many of its variants have been researched extensively because they are found to be widely applicable in practice. The VRP played a significant role in the planning of distribution, transport and logistics in many sectors such as maintenance planning, task sequencing, waste collection and mail delivery [18], [28]. Both the TSP and VRP are NP-hard problems and therefore difficult to solve to optimality [33] Quite a few exact algorithms were proposed for the VRP, such as dynamic programming [8], integer programming, branch and bound [7], branch and cut methods [24], [32] and network flow formulation [5]. Moreover, for the MDVRP, although to a lesser extent, a number of exact methods were proposed [6]. However, since the MDVRP is a complex problem, solving it by an exact algorithm is very time and memory consuming. The instance sizes exact methods can handle are therefore limited. Consequently, MDVRP variants are in practice mostly solved using algorithms that can efficiently and effectively provide near optimal solutions [19].
Several other extensions of the VRP have also gained significant research attention. These other VRP variants include: i) the Capacitated VRP (CVRP), which is the most common and basic variant (vehicle load is restricted), ii) the Heterogeneous VRP (HVRP), in which there is a mixed fleet of vehicles (different types of vehicles each with its own capacity), iii) the Vehicle Routing Problem with Pickups and Deliveries (VRPPD), in which besides customers demanding products they can also return products, similarly as in iv) the VRP with backhauls (VRPB), however in this variant all deliveries must be completed before any pickups can be done and v) the VRPTW, which is the VRP with time windows [45], [22].
Another variant, closely resembling our problem definition is the MDSDVRP, the Multi-Depot Split Delivery Vehicle Routing Problem, which was recently introduced [17]. This VRP variant combines the MDVRP and SDVRP, the Split Delivery Vehicle Routing Problem. In this last variant the delivery of a customer order can be split between multiple vehicles. The MDSDVRP allows split deliveries among vehicles based at the same depot and also between vehicles based at different depots. In contrast, this research only allows splitting deliveries between vehicles from different depots.
with routing. It is shown that a one-stage approach can provide better results. The performance for the two stage approach is highly dependent on the assignment stage. The performance de-creases when the instance size gets larger, as more wrong assignments are made [26].
Besides the exact methods, several heuristics for solving the MDVRP have been proposed. We differentiate between two classes: i) classical heuristics (for example, the savings method [9] and sweep algorithm [44], [13]) and ii) modern meta-heuristics, that provide a deep exploration of promising solution space regions. The latter category usually provides better results at the cost of an increasing computation time [23]. Several meta-heuristic used to solve large-scale practical applications of the (MD)VRP are: i) tabu search, ii) ant colony optimization, iii) particle swam optimization, iv) genetic algorithm and v) simulated annealing [27]. Tabu search enhances the performance of local search techniques by implementing memory structures. This meta-heuristic will not consider a solution repeatedly, as considered solutions are marked as tabu for a certain specified period [14]. The tabu search method has been used for the MDVRP first by [37] and [11]. This thesis will implement an idea based on tabu search. After a certain threshold of generations it will be checked if newly generated offspring resembles one of the chromosomes in the current population. The threshold will be set according to the idea that the probability of creating similar solutions increases with the generation number. For similar chromosomes several steps of the algorithm can be skipped, which saves computation time.
This thesis focusses on the genetic algorithm (GA) which has been successfully applied to a broad range of complex optimization problems [45]. The success of the GA is because of charac-teristics as simplicity, ease of operation and great flexibility [19]. Moreover, it is demonstrated that a genetic algorithm variant, namely the ‘Hybrid Genetic Search Algorithm with Adaptive Diversity Control’ (HGSADC) by [43], is the most successful approach to date for solving the MDVRP. Before that a tabu search method penalizing recurring assignments and thereby grad-ually diversifying the search by [11] was the best method for solving the MDVRP for a long time. Together with other neighbourhood-based methods such as the tabu search method by [37] and the simulated annealing method by [27] it was the most efficient method. These methods were outperformed by the adaptive large neighbourhood search method (ALNS) from [34], which implements the ruin-and-recreate idea (i.e. removal and insertion techniques) with an adaptive selection of operators.
To the best of our knowledge there are only a few other evolutionary approaches for the MDVRP to date. Two of them are based on taking advantage of geometric aspects within the problem [40], [33]. The first time a hybrid genetic algorithm was used for solving the MDVRP was by [19] in 2008. Two hybrid genetic algorithms differing in the generation of initial solutions were constructed. In one of the variants initial solutions are randomly generated in the other variant the savings algorithm by [9] and the nearest neighbour algorithm were applied to generate ini-tial solutions. Their analysis showed that using the latter iniini-tial solution generation method is superior to randomly generating initial solutions. However, these genetic algorithms are both of less quality than the ‘HGSADC’ by [43]
Split algorithm [35]) to find the optimal segmentation of the giant-tour into routes [43]. An interesting aspect of the ’HGSADC’ is the diversity management mechanism, consisting of i) the search space including both feasible and infeasible solutions, ii) a fitness measure based on both cost and a diversity contribution and iii) refreshing the population according to a prespecified rule on the number of iterations without improvement in the current generation. Because of the good results of [43], [15], this thesis will also implement a diversity management mechanism with the three discussed criteria.
To the best of our knowledge there is no literature on the MDVRP variant with an inven-tory restriction at a specific point in time. A field that does combine distribution with inveninven-tory management (over time) is the Inventory Routing Problem (IRP); here three decisions need to be made: i) when to serve customers (inventory points)?, ii) how much to deliver? and iii) which delivery routes [10]? The customers in the IRP are inventory points having local inventories and facing stochastic demand. The literature on this problem is very broad, as there is a large variety in possible assumptions when combining inventory management and routing. Assumptions are made on the supply chain type, the length of the planning horizon, the demand type (deter-ministic or stochastic), inventory management decisions on lost sales, emergency transshipments and backorders, the vehicle fleet (homogeneous versus heterogeneous and single, multiple versus unconstrained) and the solution approach. This makes it hard to exactly classify the problem Because long term dynamic and stochastic problems are very difficult to solve, the approaches found in current literature for solving the IRP have simplified the problem by considering shorter time horizons, assuming deterministic demand or looking at a simplified network with a single distribution center [4].
Because the MDVRP with the possibility of split delivery between sources in itself is already a complex routing problem this research assumes inventory is known (at a specific point in time). This research looks at the distribution of products at a specific point in time, namely the daily planning for home deliveries using current inventory levels. This contrasts the above specification of the IRP which takes in to account a specified time horizon and decides on how much and when to deliver. Moreover, whereas the IRP considers the routing from a central distribution center to inventory points having local inventories and facing stochastic demand, this research analyses the distribution from the sources to home delivery customers with known demand.
3
Problem Formulation
This section will carefully explain the assumptions and the methodology used to determine a model that satisfies the research goal. Thereafter the notation and the mathematical formal problem statement will be presented.
3.1
Assumptions
In this section, the assumptions made to define the scope of the model are presented. To provide a clear overview of the problem we first describe the most important set of assumptions defining the structure of the model. Several extensions of the model, made to provide more practical relevance, will be described in Section 5. As previously mentioned, this research analyses an MDVRP variant with an inventory restriction. This extension makes the model more realistic compared to other models, because in order to fulfil an order the source needs to have sufficient inventory. Similar as in the regular MDVRP a maximum vehicle capacity and maximum route duration is taken into account. To provide a broad approach split delivery solutions from multiple sources are allowed and a system with multiple types of sources is considered. The sources are related and home deliveries can be fulfilled via direct or indirect distribution. The implications of relaxing some of the restrictive assumptions are described in Section 9.
• A system with one central distribution center (DC), several intermediate points and mul-tiple (in-store and home delivery) customers is assumed. An extension of the model will allow the option of multiple central distribution centers (Section 5.4).
• The model considers the distribution of a single product type.
• The vehicle capacity of trucks supplying intermediate point is assumed to be unrestricted. This means there is always enough remaining capacity in the trucks to distribute home delivery orders to intermediate points.
• The objective of the model is to minimize total costs. Different cost components are included, namely transport cost and handling cost. An extension of this objective will include a (trade-off) cost associated with allocating a product to a home delivery instead of reserving it for in-store demand (Section 5.1). The transport cost is an all-in cost per kilometer to represent fuel cost/travelling time. Handling cost only arise at intermediate points and are presented as a fixed cost per unit of volume. The cost structure is set up to be able to take into account the distribution to the intermediate points, which is part of a time aspect not considered in this analysis. This way we try to provide a good building block for further analysis. No cost concerning the transport of products from the distribution center to intermediate points are included. This is because it is assumed that existing transport lanes (for normal store replenishment) are used and the truck capacity corresponding to these existing transport lanes is unrestricted.
• Home deliveries must be satisfied and are handled on a daily basis. This implicitly implies that the model can only allow scenarios where the total inventory exceeds total demand. Relaxing this assumption will be discussed in Section 5.5.
• The number of vehicles at the sources is unrestricted. An extension of the model where the number of vehicles is an input for the model is discussed in Section 5.2.
3.2
Methodology
This section describes the method we will use to search for near optimal solutions for the MDVRP with related sources and restricted inventory. As the constructed MDVRP variant is a general-ization of the VRP, which in turn is a generalgeneral-ization of the TSP [25], the problem is concluded to be NP-hard. Therefore a hybrid genetic algorithm will be presented, which has been successfully applied to a broad range of complex optimization problems [45], such as the MDVRP [43]. This type of algorithm belongs to the larger class of evolutionary algorithms, which find solutions to optimization problems using techniques inspired by natural evolution. Some techniques, or so called genetic operators, are for instance selection, crossover and mutation. This algorithm type (hybrid) combines components of several heuristics, to utilize their different properties. It combines the exploration breadth of population based evolutionary search and the improvement capabilities of local search (neighbourhood-based) meta-heuristics.
Artificial evolution as an optimization method was introduced in the sixties and became widely known through the work of Holland [20]. A genetic algorithm is a search heuristic based on imitating the evolutionary system. The heuristic starts with an initial population of chromo-somes, where each chromosome in the population is a candidate solution for the problem. Each chromosome has a set of properties, according to which an evaluation (fitness) of the solution can be made. This fitness measure represents the quality of a candidate solution. The procedure is an iterative process, in which the population evolves over time. The population in each iter-ation is called a generiter-ation. For each generiter-ation, the fitness of every individual (chromosome) is evaluated. According to their fitness, chromosomes in a population are selected for genetic operations as selection, crossover and mutation. This way the higher the quality of an indi-vidual the higher the probability of being selected and/or mutated to form a new population. This enhances evolution towards better solutions. The algorithm terminates when for instance a maximum number of generations has been formed or a satisfactory fitness level has been attained for the population. The individual with the highest fitness in the final generation is selected as the best solution for the problem.
com-ponent assures that the best total cost solution always advances to the next generation. The other chromosomes are removed. After a prespecified stopping criterion is reached the algorithm returns the best chromosome of the final generation. The stopping criterion can either be a limit on the calculation time or on the number of iterations. The best selected chromosome should represent a near optimal solution. For this MDVRP variant it should give a close to optimal fea-sible routing scheme, specifying the assignment of customer demands to sources and the routes from each source.
The algorithm described in this research is based on the hybrid genetic algorithms proposed for the MDVRP with homogeneous and related sources by respectively [21] and [43]. The hy-brid genetic search technique introduced by [43] ’Hyhy-brid Genetic Search with Adaptive Diversity Control meta-heuristic’ (HGSADC) is able to efficiently address several classes of vehicle rout-ing problems, namely the multi-depot vehicle routrout-ing problem (MDVRP), the periodic vehicle routing problem (PVRP) and the multi-depot periodic vehicle routing problem (MDPVRP). For these problems few efficient algorithms are currently available [43]. [43] shows that their method performs very well, ”for all currently available benchmark instances for the three problem classes, HGSADC either identifies the best known solutions, including the optimal ones, or new best so-lutions in terms of computational efficiency and solution quality”. An important aspect for our genetic algorithm, is the search space that includes both feasible and infeasible solutions with respect to the route duration constraint. Meta-heuristic literature indicates that allowing con-trolled exploration of infeasible solutions ensures a more easily transition between structurally different solutions and thereby enhances the search performance [15].
The key contribution of this research is that, different from current literature, the MDVRP is extended with an inventory aspect. In this thesis it is assumed that all sources have restricted inventory. Decisions on how to assign (and thereby also distribute) products are influenced by this restriction. Therefore the algorithm proposed in this thesis is able to take into account the inventory levels at a specific point in time as input for the model. This implies that only sources with sufficient inventory can be selected in the assignment scheme of a solution. This has a great influences on the solution structure and several steps of the algorithm, as the assignment of the initial population and the crossover phase.
3.3
Notation
Table 1 gives an overview of the notation used. All parameters, variables and decision variables will be declared. We aimed at selecting the notation in accordance with existing literature.
Indices
i, j Indices for the locations (customers and sources), i, j ∈ V ` Index representing the originating depots for routes, ` ∈ VD
k Index representing routes, k ∈ K Parameters
n Number of customer locations δ Number of depots
ct All-in cost for transport per kilometer
ch Fixed handling cost per product
Q Vehicle capacity in units of volume (products) T Maximum route duration in minutes
U Size of the (homogeneous) vehicle fleet Variables
qi Product demand of customer i in units of volume, i ∈ VC
τi Service time needed to serve customer i in minutes, i ∈ VC
Ij Inventory level of source j, j ∈ VD
dij Travel distance in kilometers between location i and j, i, j ∈ V
tij Travel time in minutes between location i and j, i, j ∈ V
ψ`k Total time for route k originating at depot `, ` ∈ VD, k ∈ K
Decision variables
x`kij Binary decision variable indicating whether connection (i, j) is included in route k starting at
depot `, ` ∈ VD, k ∈ K, i, j ∈ V
u`ki Positive integer decision variable denoting the position of customer location j in route k,
starting at depot l, l ∈ VD, k ∈ K, j ∈ V
yli Nonnegative integer decision variable indicating how much of the demand qi of customer i
is supplied by depot l, l ∈ VD, i ∈ VC
Other notation G Complete graph
V Set containing all locations
VC Set containing all customer locations
VD Set containing all depot locations
VR Set containing all central distribution center source locations
VI Set containing all intermediate point source locations
K Set consisting of all routes
S Search space consisting of all solutions
σ Solution to the MDVRP variant with inventory restriction, σ ∈ S |A| The cardinality of set A is denoted with |A|.
3.4
Mathematical model
The situation we consider can be described as an extension of the MDVRP, which in turn is one of the numerous researched extensions of the VRP. We will show that allowing split delivery so-lutions can be described by an integer linear programming model (ILP), by changing the format of the order input as compared to the model only allowing complete solutions.
Let G = (V, A) be a complete graph with V the total set of vertices and A the set of arcs. The set of vertices can be split up into two sets, V = VD ∪ VC, where VD is the set of depots
(all sources) and VC is the set of customer locations. In contrast to the regular MDVRP, this
research assumes two types of sources and therefore VD = VR∪ VI. By specifying two types
of sources we can model the relatedness of these sources. VR corresponds to the set of central
distribution centers, while VI represents the set of intermediary point sources. We let δ denote
the total number of depots and n the number of customer locations. Hence, the set of vertices V consists of locations V = {1, ..., δ, δ + 1, ..., δ + n}. Each customer i, i ∈ VC has a positive
demand and a service time, defined by qi and τi respectively. In this extension of the MDVRP
each source j, j ∈ VD has an inventory level denoted by Ij. The demand that can be fulfilled by
a source j is limited by this inventory level Ij
The arcs denoted with (i, j), i, j ∈ V indicate a direct travel possibility from location i to location j, with corresponding travel time tij and travel distance dij. The duration of a single
route consists of both the total travel time and the total service time to serve all customers included in the route. This duration is restricted by T , the maximum route duration. Moreover, the vehicle load is restricted by Q, the vehicle capacity.
The MDVRP can have several objectives, of which minimizing total cost and minimizing to-tal travel time are the most common. We focus on a cost minimizing objective: minimize transport and holding cost. The goal of the MDVRP is to determine a minimal cost routing scheme that satisfies the route duration, vehicle capacity and inventory constraints and serves all customers. As discussed in the assumptions we have two types of cost: i) transport costs (an all-in cost per kilometer) denoted by ct and handling costs at intermediate point sources (a cost
per product) denoted by ch (note that chonly applies to depots j ∈ VI). Let U denote the size
of the homogeneous vehicle fleet at the depots. As the number of vehicles available is assumed to be unlimited we have U = n. Let K denote the set of routes of a solution (index k). Hence, the maximum number of routes available at source l and thus the cardinality of set K is equal to n. Mathematical model: only complete order solutions allowed
• Let x`kij, ` ∈ VD, k ∈ K, i, j ∈ V denote the binary decision variable that indicates whether
arc (i, j) is included in route k starting at depot l of the routing scheme considered. x`kij,
` ∈ VD, k ∈ K, i, j ∈ V equals 1 if location j is visited directly after location i in route k
starting at depot ` and 0 otherwise.
• Let u`ki, ` ∈ VD, k ∈ K, i ∈ VC denote the position of customer location i in route k
starting at depot `.
In this formulation, for all customers i ∈ VC, the following holds:
X k∈K X j∈V x`kji=
0 ⇒ customer i is not served by source `
Now the mathematical formulation of the MDVRP with related sources and an inventory con-straint that only allows complete order solutions can be stated as follows (based on [43], [38], [6], [39]): min X `∈VD X k∈K X i∈V X j∈V ctdijx`kij+ X `∈VI X k∈K X i∈VC X j∈V chqix`kij (2) subject to X `∈VD X k∈K X i∈V x`kij = 1 j ∈ VC (3) X j∈V x`klj ≤ 1 ` ∈ VD, k ∈ K (4) X j∈V x`kij = 0 ` ∈ VD, i ∈ VD\ `, k ∈ K (5) X j∈V x`kji− X j∈V x`kij = 0 ` ∈ VD, i ∈ V, k ∈ K (6) X i∈V X j∈V (tij+ τi)x`kij ≤ T ` ∈ VD, k ∈ K (7) X i∈VC X j∈V qix`kij ≤ Q ` ∈ VD, k ∈ K (8) X k∈K X i∈VC X j∈V qix`kij ≤ I` ` ∈ VD (9) 2 ≤ u`ki≤ n ` ∈ VD, k ∈ K, i ∈ VC (10) u`ki− u`kj≤ n(1 − x`kij) ` ∈ VD, k ∈ K, i ∈ VC (11) x`kij ∈ {0, 1} ` ∈ VD, k ∈ K, i ∈ V, j ∈ V (12)
Objective (2) minimizes total transport and handling costs. Constraint (3) ensures that all cus-tomers are visited exactly once. Constraint (4) guarantees that routes (vehicles) are used at most once (single use or remain unused). Constraint (5) is necessary to make sure that all routes start and end at the same depot. Flow conservation (also ensuring routes are finished) is modelled by constraint (6). Constraints (7), (8) and (9) respectively enforce limits on the route duration, the vehicle load and the inventory usage. Constraints (10) and (11) make sure that sub-tours are eliminated and constraints (12) restricts the decision variables to be binary.
Mathematical model: complete and split delivery solutions allowed
4
Genetic Algorithm
This section will present a detailed description of the genetic algorithm proposed for the MDVRP with related sources and restricted inventory. Moreover, the solution representation, the search space and the algorithmic notation will be discussed.
In the first phase of the algorithm an initial population of feasible and infeasible individuals with respect to route duration is generated. These chromosomes are kept in two separate sub-populations. The chromosomes present a solution to the MDVRP variant, i.e. a distribution scheme that is feasible with respect to the inventory restriction. Each candidate solution is then evaluated according to a fitness measure which represents the quality of a solution. This measure is based on total cost (actual and penalty) and the diversity of a solution. In the next phase offspring is created by combining two parent chromosomes.
Parents are selected according to a tournament selection procedure. This is a genetic opera-tor that enhances fitness proportionate selection, i.e. individuals with a higher fitness value have a higher probability of selection and vice versa. In this selection mechanism a tournament is played between a specified number of randomly selected individuals. The individual with the highest fitness value ‘wins’ and is selected for recombination. After selecting the first parent, this chromosome is deleted from the set of possible parents used for selecting the second parent. By incorporating the tournament size as a parameter in the genetic algorithm, we are able to control the selection pressure, which heavily influences the GA’s convergence rate [31]. If the selection pressure is too low, it takes unnecessarily longer for the GA to find the near optimal solution. In contrast, if the selection pressure is too high, there is an increased probability that the algorithm converges prematurely to a suboptimal solution.
General Scheme Algorithm: split deliveries allowed 1 Get order and parameter input
2 Stopping Criteria Genetic Algorithm
While both the maximum computation time and the prespecified number of iterations without improvement have not been reached do
3 Generation of Initial Population
While the initial population has not reached its maximum size do
3.1 Assignment: assign customers using a distance-based assignment formula
3.2 Routing: determine a routing scheme by applying the Nearest Neighbour algorithm and perform a feasibility check
End while
4 Evaluation
5 Education (local search procedure) 6 Generation of Offspring
6.1 While both subpopulations have not reached their maximum size do
6.2 Parent Selection: select parent chromosomes using tournament selection operator 6.3 Crossover: recombine parents using the crossover with insertions mechanism (specified
in Table 7) to generate offspring
6.4 Routing: determine a routing scheme corresponding to the offspring’s assignment scheme obtained in 6.3 by applying the Nearest Neighbour algorithm and
perform feasibility check
6.5 Evaluation: evaluate offspring using the constructed fitness measure 6.6 Education: educate offspring (local search procedure)
6.7 If the current best solutions has not improved for the last prespecified number of iterations then diversify the population by selecting survivors and return to 2.1 6.8 If maximum subpopulation size reached then select survivors from current total
population (parents and offspring), then return to 5.1 6.9 End while (Stopping Criteria Generation)
7 End while (Stopping Criteria Genetic Algorithm)
Figure 2: Flow Chart Algorithm
4.1
Search Space and Solution Representation
This section discuss the search space and the solution representation. The search space S is defined as the set of feasible and infeasible solutions σ ∈ S, with respect to the route duration. With respect to the vehicle load and the inventory restriction only feasible solutions are allowed. Let K(σ) represent the set of routes composing solution σ. As previously defined, V = VD∪ VC,
where V is the set of all locations, VDis the set of all source locations and VCthe set of all customer
locations. Each route k ∈ K starts at a source vk
` (v`∈ VD), visits a sequence of nk customer
lo-cations vk
1, ..., vnkk∈ VC, and returns to the originating source v
k
`. A route is characterized by the
load q(k) =Pnk
i=1qvk
i and total duration ψ(k) = t(k) + τ (k), the sum of total travelling time and
total service time. The total travelling time t(k) is given by t(k) = tvk `v1k+
Pnk
i=1tvk
ivki+1+ tvnkv`.
The total service time τ (k) for a route equals τ (k) =Pnk
i=1τvk
i, the sum of the service times for
each customer in the route.
The algorithm allows a structure with split delivery from multiple sources. These are solutions where one or more customer locations are served partially by multiple sources. This structure, although more memory exhaustive and introducing a greater model complexity, is preferred over a structure allowing only complete delivery solutions because i) split delivery might be the only possible solution (because of the inventory restriction) and ii) scenarios may exists for which costs savings using a split delivery solution instead of a single complete delivery solution outweigh the disadvantage of the service level decrease.
the stated problem and therefore consist of two parts: i) source assignments and ii) the routes for each source. The combination of these two parts presents the final solution.
• Part i: assignments. To be able to incorporate a structure which includes the possibility of split delivery solutions, the (partial) demand sizes fulfilled for each source have to be saved instead of the binary decision answering if the assignment is true or false. This part of the chromosome consists of segments of variable length, each segment provides the customer assignments and corresponding fulfilled demand sizes for one source. By making the length segments dependent of the number of customer assignments, the coding remains flexible which saves computation time.
• Part ii: routes. Contains for each source ` a sequence of customers v with trip delimiters, obtained by concatenating all routes from source ` and not removing visits to sources. To provide a clear idea of the solution representation we will now illustrate it with a simple example. Consider a system with three sources (these can either be DC’s or intermediate points), each having six products in stock and ten customers with varying demand sizes reflected by the vector [1 1 3 2 1 4 1 1 1 1]. The system (left figure) and a possible assignment scheme (right figure) satisfying the inventory restriction for the given system are presented in Figure 3. Green dots (located at sources) represent products in inventory, while red dots (located at customers) represent product demands.
Figure 3: Example: System Overview (left) and possible assignment scheme (right)
The solution representation for chromosome part i, corresponding to the assignment solution given in Figure 3 (right), will now look as presented in Table 3. It can be seen that this part of the chromosome consists of segments of variable length, which for each source specify the customers and corresponding (partial) order sizes linked to the source. Note that this assignment is an example of a complete order solution assignment; all customer demands are served completely from a single source.
Source ` 1 2 3
Customers c {1 (1) ,3 (3) ,7 (1) ,8 (1)} {2 (1) ,4 (2) ,10 (1) } {5 (1) ,6 (4) ,9 (1)} Table 3: Chromosome: assignment representation (complete order solution)
Figure 4: Possible Assignment Scheme (left) and Corresponding Routing Scheme (right)
Source ` 1 2 3
Customers c `1 1 3 7 8 `1 `2 4 10 2 `2 `3 9 6 5 `3
Table 4: Chromosome: routes representation
4.2
Algorithm Notation
Table 5 gives an overview of the algorithmic notation used. We aimed at selecting the notation in accordance with existing literature.
Parameters
µ Minimal subpopulation size
λ Number of offspring generated per generation
nclose Number of close individuals considered for diversity evaluation
nc Proportion of close individuals considered for diversity evaluation, such that nclose= nc · µ,
nc ∈ [0, 1]
h Granularity threshold; defining the size of the neighbourhood set containing closest locations considered in the improvement phase, h ∈ [0, 1]
Tmax Maximal duration in minutes defining a stopping criteria of the algorithm
ITactual Current number of iterations (chromosomes created)
ITtot Total number of iterations allowed defining a stopping criteria of the algorithm
ITactdiv Actual number of iterations without improvement in a row
ITdiv Total number of iterations without improvement allowed in a row, defining the condition to
diversify the population
φ Diversification parameter, such that ITdiv= φ · ITtot, φ ∈ [0, 1]
T S Tournament size in the parent selection phase
Pf eas Proportion of feasible solutions generated in the assignment phase
Pe Constant probability operator for education
Binh Binary decision specifying if the inheritance option providing more convergence is on or off,
Binh∈ 0, 1
4.3
Generation of Initial Population
The first phase of the algorithm entails the generation of an initial population of chromosomes. The assignment schemes of the initial chromosomes are generated based on a distance formula, thereafter routes are made by implementing the nearest neighbour algorithm, the initial chro-mosomes are evaluated and are considered for education. Similar as in [43], the size of the initial population is 4µ chromosomes, where µ is the minimal size of the subpopulations. In the initial generations of chromosomes it might be possible that one of the subpopulations contains less than µ chromosomes. This can be regulated by the parameter Pf eas, specifying the proportion
of feasible solutions generated in the routing phase.
4.3.1 Assignment
The assignment procedure of the initial population of the heuristic is presented in Table 6. The assignment procedure only generates feasible solutions with respect to the inventory constraint. Therefore this assignment procedure is different from the assignment procedures in current liter-ature as they do not take the inventory levels into account. Assignment schemes are created by randomly assigning customers to sources using a distance based formula. For a selected customer i, in the first phase only sources with sufficient inventory are considered for possible assignments. This makes the assignments biased towards complete solutions. Assigning customers based on distance will decrease the appearance of infeasible solutions, ceteris paribus the parameter set-tings and the input.
The distance based formula is such that the source with the highest distance to the selected customer is selected with the lowest probability and vice versa. First, the distances from the selected customer to possible sources are found, di`, for all sources ` ∈ VD. Then the probability
of assigning customer i to source ` is found by an adapted version of the roulette wheel approach by [16] defined by: P (`) = X `∈VD di` di` X `∈VD X `∈VD di` di`
Assignment Scheme
STEP 1: Fill Sources With Complete Customer Orders 1.1 For each customer c selected in random order do
1.2 If there exists sources ` for which 0 <Demand(c)<Inventory(`)
1.3 Assign c to one of the sources with sufficient inventory using a distance based formula 1.4 Decrease inventory source ` and demand selected customer c
1.5 End if
1.6 End for
STEP 2: Assign Customers With Unfulfilled Demand (⇒ Partial Orders)
2.1 For each customer c, with unfulfilled demand after step 1, selected in random order do 2.2 While Demand(c)> 0 and exist sources ` with Inventory(`)> 0
2.3 Assign customer c randomly to one of these sources or using a distance based formula 2.4 Decrease inventory source ` and demand selected customer c
2.5 End while
2.6 End for
Table 6: Assignment Procedure Initial Population (split deliveries allowed)
4.3.2 Routing
Now that the assignment schemes of the initial population are generated, routes are created using the nearest neighbour algorithm. This approach is similar as in [21], we do however extent the nearest neighbour algorithm with a route duration constraint to control the feasibility of routes created. Moreover, a threshold parameter, specifying a generation number, for checking resemblance of the newly created offspring is implemented. When the threshold is reached, the assignment scheme of the offspring will be compared with the assignment schemes available in the current parent population. If a resembling assignment scheme is found, the route solution is directly copied, saving computational time. A threshold is implemented, since the probability of creating a similar chromosome increases as the algorithm proceeds and converges.
In the routing phase the parameter Pf eas specifies the proportion of feasible solutions with
respect to the maximal route duration generated in this phase. This parameter provides con-trol over the number of feasible and infeasible solutions allowed. As mentioned before literature indicates that allowing a controlled exploration of infeasible solutions enhances the search per-formance.
At the end of the routing procedure chromosomes will be checked on their feasibility with respect to the route duration, similarly as in [43]. If al routes of the routing scheme satisfy the maximum route duration the chromosome is marked feasible. If one or more routes of the routing scheme does not satisfy the maximum route duration specified, the chromosome is infeasible. According to their feasibility they will be placed in the corresponding subpopulation.
4.4
Evaluation
are given by:
Total Costs(σ) = Total Actual Costs(σ) + Total Penalty Costs(σ). 4.4.1 Total Actual Costs
The actual costs are build up by the transportation costs and the handling costs (only for fulfil-ment by intermediate points), which is similar as in [21]. The transportation cost of a solution are found by multiplying the total distance of all routes in the solution by an all-in transportation cost per kilometer ct. The transport cost of solution σ are given by:
Transport Costs(σ) = X `∈VD X k∈K X i∈V X j∈V ctdijx`kij(σ),
where dij is the travel distance between location i and location j, and x`kijis the binary decision
variable indicating whether arc (i, j) is included in route k starting at depot `. The handling costs, with cost parameter chper product, are given by
Handling Costs(σ) = X
`∈VI
X
i∈VC
chy`i(σ),
where yi` is the nonnegative integer decision variable indicating how much of the demand qi of
customer i is fulfilled by source `.
Hence, the total actual costs are calculated by: Total Actual Costs(σ) = X
`∈VD X k∈K X i∈V X j∈V ctdijx`kij(σ) + X `∈VI X i∈VC chy`i(σ).
4.4.2 Total Penalty Costs
The total penalty costs are the cost concerned with violating the route duration constraint. In-feasible solutions are penalized more heavily as the algorithm advances, as we do not consider them as suitable final solutions. An extension of the model concerned with additional penalty costs for violating the available number of vehicles is presented in Section 5.2.
The penalty costs for infeasibility are given by the total amount of time in minutes (of all routes) that exceeds the maximum route duration, multiplied by a dynamically changing penalty param-eter ω. By incorporating a dynamic penalty paramparam-eter the preference of feasible chromosomes increases during the execution of the algorithm. This set-up is a small adaptation of the one pre-sented in [21], as we base it on the number of iterations (i.e. creating one chromosome, similar as in [43]) instead of the number of iterations without improvement. Initially the penalty parameter is set to 1 and thereafter is adjusted according to the number of chromosomes generated. The structure of the dynamic parameter, ω ∈ [1, 2], is as follows:
ω = 1 + ITactual ITtot
,
where ITactual reflects the current iteration number and ITtot is one of the stopping criteria of
the algorithm, namely the total number of allowed iterations. Now the total penalty costs are given by:
Total Penalty Costs(σ) = ω · X
`∈VD
X
k∈K
where ψ`k is the total duration (driving and service time) of route k originating from depot ` and
T is the specified maximum route duration. We would like to note that the penalty structure can easily be adapted to user preferences as explained in Section 5.3.
4.4.3 Diversity Contribution
The diversity contribution of chromosome reflects its contribution in variation to the current gene pool. It is calculated as the average difference in assignment schemes between the considered chromosome and its nclose closest neighbours, similar as in [43]. Closest neighbours are defined
as the chromosomes with the most resembling assignment scheme as the chromosome considered. The number of closest neighbours considered is set by nclose= nc · µ, where nc is the proportion
of close individuals considered for diversity evaluation, nc ∈ [0, 1] and nclose is rounded. Note
that since the diversity contribution of a chromosome depends on all chromosomes present in the current population, this procedure can only be called once the assignment schemes of all chromosomes in the initial population are known. The average difference θ(σ1, σ2) between two
complete delivery solutions σ1 and σ2 is calculated by ([43]):
θ(σ1, σ2) = 1 n X i∈VC [Iχi(σ1)6=χi(σ2)],
where n is the number of customers with demand, and I is an indicator function returning one if the source assignment χi(σ1) of customer i in chromosome σ1, is not similar to the assignment
χi(σ2) in chromosome σ2. The diversity contribution ∆(σ1) of chromosome σ1 to the current
population can now be stated as: ∆(σ1) = 1 nclose X σ∈Nclose θ(σ1, σ),
where Nclose is the set with the nclose closest assignment scheme neighbours of the chromosome
considered, σ1.
When split delivery solutions are involved the average difference needs to be calculated dif-ferently, as an assignment for a single customer can contain multiple sources. The new approach comes down to finding the sum of the cardinalities of the difference sets of source assignments and dividing the result by two to avoid double counting. These difference sets are sets of variable length containing all customers assigned to the source in only one of the solutions considered. Let U`(σ
1, σ2), I`(σ1, σ2), D`(σ1, σ2) respectively denote the union set, the intersection set and the
difference set of source ` corresponding to solutions σ1and σ2. Now the new average difference,
θ∗ can be stated as:
θ∗(σ1, σ2) = 1 n " X `∈VD |U`(σ 1, σ2) − I`(σ1, σ2) | # /2 = 1 n " X `∈VD |D`(σ 1, σ2) | # /2
Consequently, the new diversity contribution ∆∗of solution σ1to the current population can be
When calculating the average difference between two complete solutions we stick to the previously explained method of finding ∆, which showed to be faster in AIMMS. The new method involves extensively checking cardinalities of indexed sets (sets dependent on multiple indices), which was found to be time consuming in AIMMS. We note that this new method is also valid for calculating the average difference between two complete solutions; it gives the same result as the standard method.
4.4.4 Evaluation Function
The evaluation function is based on the total costs and the diversity contribution. This fitness measure is used to perform two selections, i.e. parents for crossovers and chromosomes to ad-vance to the next generation (survivors). The function is based on [21], although implemented differently:
EV (σ) = ξtc· [|P op| + 1 − RAN K (Total Costs(σ))] + ξdc· [|P op| + 1 − RAN K (−∆(σ))] ,
where ξtcand ξdcare dynamic weight parameters, RAN K (Total Costs(σ)) and RAN K (−∆(σ))
present the rank of chromosome σ in the current population with respect to the total costs and diversity contribution respectively and |P op| is the size of the current population. By basing the function on the number of ranks to distribute in total minus the appropriate rank of the chromosome considered, we i) avoid favouring newly created offspring whilst ii) maintaining the goal to give the highest fitness value to the best chromosomes. To give the highest rank to the fittest individual the ranks are based on the total cost and the negative diversity contribution. The dynamic weight parameters are set such that the weight on the cost increases and the weight on the diversity contribution decreases during the execution of the algorithm. These parameters are defined similarly as in [21]:
ξtc= 0.5 + 0.5 · ITactual ITtot , ξdc= 0.5 − 0.5 · ITactual ITtot
4.5
Education
After evaluation chromosomes are educated using several improvement procedures (local search procedures), as specified in [43]. This way the routes constructed by the nearest neighbour algo-rithm can be improved. Important to note is that the education operator is modelled such that it can be called on a source level. The improvement steps considered only apply between routes originating at the same source. Thereby assignment schemes remain intact and feasibility with respect to the inventory constraint is guaranteed.
Let u be a customer location, i.e. u ∈ VC and let k(u) denote the route containing customer
location u. Let (u1, u2) denote the partial route from u1 to u2. Note that this means u1 and
u2 are served in the same route. Define the neighbourhood of customer u as the hn closest
locations as well as the considered source.
The improvement phase iterates for each source in random order over all assigned customers u and each of its neighbours v. The improvement steps require that u 6= y, v 6= x and that only when considered in the called improvement step, the successors are customer locations, i.e. x, y ∈ VC. The improvement steps are randomly investigated and for each selected customer u
and neighbour v the first improvement in total costs is implemented. Whether or not an improve-ment is made is only based on the total costs, as assignimprove-ment schemes and thereby the diversity contribution does not change by the procedures. The process of selecting u and v continues until all customers assigned to a source are considered. Then the operator moves to the next source to continue the improvement phase. The considered route improvement moves are the following: M1 Remove u and place it after neighbour v
M2 Remove (u, x) and place (u, x) after v M3 Remove (u, x) and place (x, u) after v M4 Swap u and v
M5 Swap (u, x) and v M6 Swap (u, x) and (v, y)
M7 Replace (u, x) and (v, y) by (u, v) and (x, y) M8 Replace (u, x) and (v, y) by (u, y) and (x, v)
The first three moves are insertions and the next three moves are swaps. The last two moves are 2-opt intra-route moves if k(u) = k(v) and 2-opt inter-route moves if k(u) 6= k(v). Stating it this way, all moves can be applied indifferently on whether locations are in the same or in a different route originating at the same source.
To identify whether a step indeed results in an improvement, after each iteration the new total cost need to be calculated. An improvement is only implemented if the new distribution scheme remains feasible with respect to the vehicle capacity. For every improvement the feasibility check needs to be applied and the subpopulations will be updated accordingly. The education operator is called with probability Pe.
4.6
Generation of Offspring
4.6.1 Stopping Criterion Generation
The stopping criteria specified are similar as in [43]. The population of parent chromosomes is ex-tended by creating offspring while the two subpopulations have not reached their maximum size. When the stopping criteria is triggered it will first be checked if the number of iterations in a row without improvement, ITactdiv, exceeds the allowed number of iterations in a row before
diver-sification takes place, ITdiv. If the check is true a small number of survivors is selected based on
chromosomes for the new generation without needing to verify whether they are improvements. Consequently, the time consuming diversity contribution only needs to be called once for the complete new population. This instead of calling the procedure per newly generated chromo-some (to verify if it is an improvement) and again once all chromochromo-somes are created, since the diversity contribution needs to be up to date before selecting survivor chromosomes, as in [21].
4.6.2 Parent Selection
To select parents for recombination we implement the tournament selection operator, as discussed in [31]. In contrast to [43] that implements binary tournaments, in this research the tournament size will be a tuning parameter. This way the selection pressure can be controlled. By increasing the tournament size, the degree to which fitter individuals are favoured is increased. Vice versa, decreasing the tournament size will lead to a lower selection pressure. The convergence rate of a GA is proven to be largely determined by the selection pressure. Incorporating the tournament selection operator with the tournament size as parameter, we can control the convergence rate of the genetic algorithm. Therefore, we select this approach instead of the also well known roulette wheel approach from [16], implemented by [21]. In this last approach the probability for an individual to participate in a crossover is directly proportional to his relative fitness measure and therefore the convergence rate cannot be controlled easily.
In the tournament selection procedure two tournaments of size T S are played. All chromosomes in the current population are candidates in the first tournament. T S number of individuals are randomly selected from the candidate pool. The candidate with the highest fitness value wins the tournament and as a prize gets to be the first parent. Then the second tournament of size T S is played similarly. However, the selected first parent chromosome is excluded from the tournament pool. The outcome of the second tournament completes the couple selected for recombination.
4.6.3 Crossover
After selecting two parents their genetic material will be recombined to produce offspring. This paper proposes a new crossover with insertions mechanism designed to transmit good groups of customer assignments for a source, while enabling assignment recombinations and satisfying the inventory constraint. Compared to the research of [21] we need a more complex crossover mechanism, since feasibility already has to be accomplished during crossover (in the assignment phase) instead of being able to create a feasible solution after crossover in the routing phase (by for example starting a new route when the route duration of vehicle capacity is exceeded). The method is based on the ’periodic crossover with insertions’ (PIX) crossover procedure from [43]. However, the structure is completely different since i) this research does not consider a periodic problem, ii) we have assignment crossovers instead of route crossovers (sequences of vis-its) and most importantly iii) this paper takes into account the additional restriction of limited inventory for each source. This last characteristic complicates the crossover significantly, as only feasible solutions with respect to this restriction are allowed. Moreover, thought has to be given on how to deal with crossovers regarding split delivery parent solutions.
children chromosomes to inherit genetic material from both parents in nearly equal proportions, while the second capability is guaranteed by inheriting most material from one parent” [43].
Crossover with Insertions: split deliveries allowed
STEP 0: INHERITANCE RULE
0.1 Draw two random numbers between 0 and # sources according to a uniform distribution. Let n1 and n2 be respectively the smallest and the largest of these numbers
0.2 Randomly select n1 sources to form the set Λ1
0.3 Randomly select n2− n1 remaining sources to form the set Λ2
0.4 The remaining # sources−n2sources form the set Λmix
STEP 1: INHERIT GENETIC MATERIAL FROM P1
1.1 For each source, selected in random order, contained in the set:
1.2 Λ1: For each customer assignment in Parent 1 selected in random order
1.3 If customer demand satisfies inventory left in source s after inherited assignments from P1then copy complete customer assignment from P1 to C
1.4 End for
1.5 Λmix: Draw a random number, n3, between 0 and the # of sources contained in the
set Λmixaccording to a uniform distribution
1.6 Randomly select n3 customer assignments to copy from P1to C
1.7 Assign each of the n3customer assignments in random order to C, if the customer
demand satisfies the inventory left in source s after previous assignments from P1
1.8 End for
STEP 2: INHERIT GENETIC MATERIAL FROM P2
2.1 For each source s selected in random order 2.2 For each customer assignment in P2
2.3 If customer not yet served and demand satisfies inventory left in source s after inherited assignments from P1and possibly P2then copy customer assignment from P2to C
2.4 End for
2.5 End for
STEP 3: COMPLETE CUSTOMER SERVICES
3.1 Assign complete orders
3.2 For each customer with unfulfilled demand after inheritance from P1 and P2, selected in
random order
3.3 Construct a set containing all sources with sufficient inventory for the complete order. If this set is not empty then:
3.4 Option i) Randomly select one of the possible sources
Option ii) Randomly select one of the parent 1 sources if possible, else apply option i)
3.5 End for
3.6 Assign partial orders if necessary (results in split deliveries)
3.7 If there exist customers with unfulfilled demand then assign customers with unfulfilled demand, selected in random order, to sources with inventory left (now: random). This step (if necessary) assigns partial demands to multiple sources.
In the first step, STEP 0, inheritance rules for the considered crossover are specified. Random numbers are drawn to construct three mutually exclusive sets of sources. A set with sources for which assignments which will completely (if possible) be taken over from parent one, Λ1. A set
with sources for which material from both parents will be used, Λmix. And a set with sources
for which parent one material will not be considered, Λ2.
In STEP 1 the first parent is targeted. For each randomly selected source either the complete material (if the inventory allows) will be copied if the source is in set Λ1 or a random number
of customers if the source belongs to Λmix. For each selected source filled with customers the
copying takes place by going randomly over the customers. Since we want to limit the number of split delivery solutions, we only allow split delivery parts in the genetic material to advance one generation. Therefore, in STEP 1 only complete customer orders are copied. If a customer is selected for which a partial demand is assigned to the source, the complete customer demand (if possible) will be assigned in the offspring chromosome. Consequently, for each iteration its necessary to check if copying is still possible, as the inventory is limited.
In STEP 2 the offspring inherits material from parent two. To make the offspring take over as much material from its parents as possible, this step considers all sources. Again iteratively all sources and corresponding customers are considered. For each iteration we need to check if the considered customer still has demand (it might already been copied from parent one) and if the considered source has enough inventory left. Sources in Λ1 and Λmix might already been
filled to a great extent.
Figure 5: Crossover with Insertions Mechanism (complete order solution example) In the example of the crossover mechanism presented in Figure 5 the inheritance rules of STEP 0 have specified the following sets: Λ1 = {`2}, Λmix = {`1, `3} and Λmix = ∅. In STEP 1
inheritance of parent 1 results in completely taking over customers 2 and 3 which are assigned in source `2∈ Λ1. Moreover, customers assigned to sources in Λmixare partially copied in this step.
This leads to copying customers 5, 6 and 8 from source `1 ∈ Λmix and none of the customers
of source `3 ∈ Λmix. Each iteration is done step by step to check if the available inventory
is sufficient to serve the considered customer. Then in STEP 2 inheritance from parent two takes place. The sources are considered in random order, in this case the order is `1, `3, `2. At
each iteration its necessary to check if the randomly selected customer still has demand and if this demand is less then the inventory of the source considered. The order of customers is as follows: `1 : 9, 5, 10, 7, 8 followed by `3 : 6, 1 and then `2 : 3, 2. In the figure we find that only
4.6.4 Routing
The nearest neighbour algorithm is applied to obtain a routing scheme corresponding to the assignment scheme of the newly generated offspring. The feasibility of the newly chromosome is is investigated. The size of the corresponding subpopulation is updated.
4.6.5 Evaluation
Once all chromosomes of the new generation are created evaluation of the complete new gener-ation takes place. First the costs and diversity contribution are calculated. Next, the fitness of the chromosomes can be calculated.
4.6.6 Education
Educate chromosomes using route improvement procedures (local search).
4.6.7 Condition Triggering Diversification
When the maximum number of iterations without improvement in a row is reached the current population will be diversified. After a new generation of chromosomes is completed, it will be checked if the current number of iterations without improvement in a row ITactdiv exceeds the
allowed number of iterations in a row ITdiv [43]. This in contrast to [21] where the condition
is based on the iterations without improvement in a row within one generation. We specify ITdiv = φ · ITtot, where φ ∈ [0, 1] is the diversify parameter. This specification avoids drastic
behaviour in calling the diversification procedure. If using the count within a generation, the parameter ITdiv has to be set very sharply. In this specification it can easily happen that the
diversification procedure is called intensively once good solutions are found. During the algo-rithm, more good solutions are found and therefore the probability of a new solution being an improvement gets smaller. This leads to a random search in the solution space and does not provide justice to the genetic algorithm setting. Hence, we call the diversification procedure in a controlled manner. The counter ITactdiv is set to zero every time diversification takes place.
When the diversification operator is called, similar as in [43], the fittest µ/3 chromosomes from each population are selected as survivors. Additionally, elite preservation is incorporated to also select the best total cost chromosome as survivor. This way its prevented that excellent solutions found in first phases of the algorithm are lost because of the significant weight of the diversity contribution on the fitness. With the group of survivors the process re-enters the stage of cre-ating an ‘initial’ population of chromosomes. The population is diversified until it has size 4µ again, whereafter the process continues with this new population.
4.6.8 Stopping Criteria Generation
4.7
Stopping Criteria Genetic Algorithm
The algorithm finishes when one of the following two stopping criteria is triggered: • The specified maximum execution duration in minutes Tmaxis reached
