Vehicle routing decision support for a local retailer

One of the most challenging decisions that has to be made routinely by dispatch managers at distribution centres of warehousing and distribution businesses in the retail sector involves the assignment of delivery vehicles to service customers exhibiting demand for retail goods and the subsequent routing of these delivery vehicles to the various customers and back again. Perhaps surprisingly, these dispatch managers do not always use vehicle routing software to schedule goods deliveries to customers, instead often relying on teams of human schedulers who perform this task manually. The reason for not using such software is usually a perception that it may be diﬃcult to integrate the software with existing enterprise resource planning systems already in use. In such cases, estimates of potential costs savings that may be brought about by appropriate software is often required before the dispatch department will risk the signiﬁcant step towards investing in vehicle routing planning software. Dispatch managers may then employ these cost savings estimates in cost-beneﬁt trade-oﬀ analyses. This paper contains a practical case study in which the potential cost savings of a vehicle routing optimisation approach are quantiﬁed for a large retail distribution centre in the South African Western Cape in a bid to support its decision as to whether or not to invest in vehicle routing planning software.


Introduction
It is well known that distribution centres (DCs) form a crucial part of a retail supply chain, facilitating the shipment of retail goods from suppliers to end customers. Products are usually received from different suppliers at a DC which acts as a central hub at which all of these products are consolidated, sorted, picked and dispatched in appropriate quantities to customers [22]. This type of distribution hub makes sense from both the perspectives of reducing the lead time associated with customer demand fulfilment and leveraging economies of scale by the retailer in question.
The management of operations at a retail DC is typically complex, requiring challenging decisions in respect of DC layout, inventory replenishment, and customer order picking, to name but a few. But perhaps one of the most challenging decisions that routinely has to be made by the dispatch departments of warehousing and distribution businesses in the retail sector involves the assignment of delivery vehicles to service customers exhibiting demand for retail goods and the subsequent routing of these delivery vehicles from the DC to the various customers and back to the DC again. The SPAR Group Ltd is one such a warehousing and distribution business. It is listed on the Johannesburg Securities Exchange (JSE) [31]. Headquartered in Durban, South Africa, SPAR Group Ltd operates in nine different countries, has ten DCs and serves 3 768 SPAR retail stores in total.
One of its DCs is located in Phillipi, in the Western Cape. Servicing 284 South African SPAR stores in January 2020 alone, the SPAR Western Cape DC scheduled between 174 and 303 individual deliveries on any given day of the month, save for Saturdays [1]. Deliveries are usually scheduled the day before, and no deliveries are made on Sundays. The stores serviced by the SPAR Western Cape DC are situated in the Western Cape and portions of the Northern Cape. The DC has access to an active fleet of more than 40 vehicles of varying sizes, as may be seen in Figure 1, which services stores within a travel radius of approximately 150km from the DC [10]. Additional third-party vehicles may, however, also be hired and are typically used to service more distant stores. Delivery vehicle dispatch decisions at retail DCs typically involve solving instances of the celebrated vehicle routing problem (VRP). The VRP is a combinatorial optimisation problem in which a set routes is sought for a fleet of delivery vehicles that collectively services a number of customers 1 and which minimises the total distance travelled by all the vehicles. The archetypal VRP prototype is the capacitated VRP with time windows (CVRPTW) in which each vehicle additionally has a specified capacity (upper bound) in terms of the volume of goods that it can transport to customers and in which visits by vehicles assigned to service customers have to be timed so that these customers are serviced within pre-specified time windows. Large instances of the CVRPTW are notoriously difficult to solve exactly and the highly combinatorial nature of these problem instances makes it extremely difficult to construct high-quality approximate solutions manually. A variety of specialised commercial software options therefore exist for supporting delivery vehicle dispatch decisions at DCs [9,23,24,28,29].
The SPAR Western Cape DC does not, however, use vehicle routing software to schedule its deliveries. Instead, a team of human schedulers performs this task manually. While this may come as a surprise to the reader, it is our experience that this seemingly precarious situation is often the case in the South African retail sector. Reasons for not using specialised software to streamline vehicle routing dispatch decisions may be related to the often prohibitive cost of such software or may simply be the result of such software being difficult to integrate with existing enterprise resource planning (ERP) systems already in use. It is often difficult for a retailer to take the significant step towards investing in VRP planning software for the aforementioned reasons, among others.
We were therefore tasked by the SPAR Western Cape DC to assist its management in determining whether such an investment step would be appropriate in terms of the expected cost savings that may be brought about by VRP planning software. We did this by quantifying the range of potential savings in the cost associated with its delivery vehicle routing decisions that may be realised if a VRP optimisation approach were to be adopted when performing these decisions as opposed to performing these decisions manually based on (considerable) historical route scheduling experience. We anticipate that the experience that we gained during this process may be useful for operations research consultants and so we decided to write up this cost saving quantification process as a practical case study.
The remainder of this paper is organised as follows. Upon having conducted a brief review in §2 of the literature with the aim of identifying powerful CVRPTW solution approaches capable of solving realistically sized CVRPTW instances approximately, we derive a mixed-binary programming CVRPTW model in §3 that we subsequently set about solving in the context of the SPAR Western Cape DC dispatch decisions with a view to quantify the aforementioned potential cost savings. Thereafter, we show that the computational complexity of this model is such that its exact solution seems to be out of reach of a state-of-the-art mixed-integer programming solver and standard personal computing technology. Based on our findings during an exploratory literature review, we describe an acclaimed hybrid metaheuristic in §4 selected for the purpose of solving approximately three representative CVRPTW model instances in a special case study in §5, involving dispatch data obtained from the SPAR Western Cape DC. In doing so, we demonstrate that the hybrid metaheuristic is capable of achieving considerable dispatch cost savings. The paper finally closes with a number of concluding remarks in §6.

Literature review
In their most basic form, exact methods for solving VRP instances are based on the branch-and-bound method, the branch-and-cut method and the branch-and-price method. The last half-century has, however, seen a monumental volume of research on the design and improvement of methods for solving VRP instances exactly. While efforts continue to develop more powerful exact mathematical programming techniques, usually involving the implementation of ever-more specialised and sophisticated types of cuts, Toth and Vigo [35] have claimed that these exact methods can only reasonably be applied to VRP instances of approximately one hundred customers or fewer. For larger problem instances, the use of (meta)heuristics becomes necessary in order to obtain high-quality solutions within a reasonable time.
Heuristics have been a part of VRP research for as long as the VRP itself has been researched explicitly. In their landmark 1959 paper, Dantzig and Ramser [8] described the VRP for the first time, as well as a simple heuristic for solving VRP instances that seeks to eliminate partial solutions through trial and error. Just as newer and better exact VRP algorithms have been proposed over the years, so VRP heuristics have also made significant advances. One of the most popular VRP heuristics, however, remains the savings method proposed in 1964 by Clark and Wright [4].
Simple heuristics are, however, constrained in how useful they are on their own. In 1986, Glover [11] coined the term metaheuristic to describe a category of heuristic search methods that pursue alternative solutions even after having encountered local optima. These methods are often considered to be heuristics guiding other heuristics. Since the dawn of the new millennium, metaheuristics that have appeared in the literature generally reside in one of two major classes: Trajectory-based metaheuristics and population-based metaheuristics. The former iteratively seek to improve a single solution by exploring solutions locally, while the latter maintain a population of simultaneous solutions that are recombined iteratively to create new populations of solutions [40].
Recently, hybrid metaheuristics have risen to prominence in the VRP literature. These methods seek to combine the strengths of both trajectory-based and population-based approaches in order to provide more powerful and more flexible algorithms. Hybrid metaheuristics can be applied to a large variety of VRP variations without the need for implementing major algorithmic changes in order to accommodate particular problem attributes. Their strengths include allowing for VRP instances to be solved in consistently shorter times and returning approximate solutions that are less than a percent away from the best solutions documented for well-known benchmark problems [35].
Vidal et al. [40] identified the most common metaheuristics applied to VRP instances exhibiting a variety of attributes in an important survey paper of 2013. They documented the frequencies of use of these metaheuristics, as summarised in Table 1. Hybridisations of these metaheuristics have also been applied frequently (thirty nine of the sixty five algorithms surveyed, in fact, featured some form of hybridisation).  Both Toth and Vigo [35] and Vidal et al. [40] presented computational results obtained by a variety of metaheuristics in the context of the two most commonly adopted VRP benchmark data sets. The algorithmic performance data reviewed by Toth and Vigo [21] are recounted here, as much of these data appear to be the same as those considered by Vidal et al. [40]. In particular, the performance data presented were obtained by the metaheuristics listed in Table 2 in the context of a data set proposed earlier by Golden et al. [12] and a data set proposed by Christofides et al. [2]. The former data set contains the larger instances of the two, with twenty instances each containing between 240 and 483 customers, while the latter data set consists of fourteen instances, each containing between 50 and 200 customers. Since the number of customers encountered later in the case study of this paper ranges from 200 to 300 customer stops, considering the larger set of instances would seem more appropriate in the context of this paper.
Acronym Reference Algorithmic approach CLM01 Cordeau et al. [5] Tabu search TV03 Toth and Vigo [34] Granular tabu search RDH04 Reimann et al. [27] Ant colony optimisation T05 Tarantilis [33] Adaptive memory + tabu search MB07 Mester and Bräysy [19] Evolutionary algorithm + local search PR07 Pisinger and Ropke [25] Adaptive large neighbourhood search NB09 Nagata and Bräysy [21] Hybrid genetic algorithm P09 Prins [26] Greedy randomised adaptive search procedure + evolutionary local search GGW10 Groër et al. [13] Route-to-route + evolutionary computation ZK10 Zachariadis and Kiranoudis [41] Guided local search + tabu search GGW11 Groër et al. [14] Parallel route-to-route CM12 Cordeau and Maischberger [6] Parallel iterated tabu search JCL12 Jin et al. [15] Parallel cooperative tabu search VCGLR12 Vidal et al. [38] Hybrid genetic algorithm SUO13 Subramanian et al. [32] Set partitioning + iterated local search Table 2: Top-performing VRP metaheuristics compared by Toth and Vigo [35]. Adapted from Toth and Vigo [35] and from Vidal et al. [40]. Table 3 contains a summary of the average percentage gap achieved between each algorithm's best solution and the overall best solution known for the benchmark instances, along with the actual run time and normalised 2 run time in each case. These average gap percentages and normalised run times are presented graphically in Figure 2.
According to the results in Figure 2, the algorithms achieving the most impressive trade-offs between computation time expended and solution quality achieved (P09, VCGLR12-A, VCGLR12-B and MB07) are all hybrid metaheuristics. A plausible explanation for the success of hybrid metaheuristics is that they combine the strength of population-based algorithms in terms of a thorough exploration of the search space with the strength of trajectory-based algorithms in terms of exploitation of local optima.
In 2013, Vidal et al. [39] proposed the Hybrid Genetic Search with Advanced/Adaptive Diversity Control (HGSADC) algorithm for solving CVRPTW instances. This algorithm is similar to the one proposed the year before by Vidal et al. [38], but features the addition of decomposition phases to the HGSADC algorithm which are beneficial when solving large problem instances. Results obtained by the original version of the algorithm (the square data points VCGLR12-A and VCGLR12-B) appear on the Pareto frontier of the   Table 2 when solving the VRP benchmark instances of Golden et al. [12]. The percentage gap from the best-known solution, the normalised run time, and the clock time are shown. Adapted from Toth and Vigo [35].
benchmark results in Figure 2 (the data points in red). The HGSADC algorithm achieved the highest-quality solutions within a reasonable time (of respectively no more than 75 minutes and 17 minutes, on average, per instance). These data made the 2013 algorithmic variation our algorithm of choice when solving the (relatively large) CVRPTW instances considered in the case study of this paper.

Mathematical model
In the spirit of reproducibility and for the sake of future comparisons of results, this section is devoted to a detailed discussion on the CVRPTW model considered in this paper. The model is derived in the form of a mixed-binary programming problem in §3.1 and its implementation in a state-of-the-art mixed-integer solver and subsequent verification are described in §3.2. The complexity of the implementation is finally measured in §3.3. Average gap (%)

GGW11-
Normalised running �me (s) Figure 2: Computational results reported for the metaheuristics listed in Table 3 when solving the benchmark instances of Golden et al. [12], displayed in a log-linear scatter graph. The percentage gaps from best-known solutions and normalised run times (measured in seconds) are displayed. Adapted from Toth and Vigo [35].

Model derivation
The CVRPTW model considered in this paper is derived in this section. The various input parameters required for an instance of the model are described in §3.

Model parameters
Let the set L = {0, 1, . . . , m, m + 1} index m customers who have to be serviced from a central depot. As is customary in the VRP literature, the depot itself is represented by the indices 0 and m + 1 for vehicle departure and return, respectively. "Customers" 0 and m + 1 therefore share the same location. Also, denote the number of pallets of demand experienced by customer j ∈ L \ {0, m + 1} by Q j , with the demand of the depot being zero.
Denote the start and end times of the time window during which a vehicle may service customer j ∈ L by W S j and W E j , respectively, both measured in hours after midnight. If any customer j does not specify such a time window, the parameter values W S j = 0 and W E j ≥ 24 are assumed (this is the case for the depot as well, as it is always open). Denote the average unload time per pallet at customer j ∈ L \ {m + 1} by U j , measured in hours, and the expected time required to load a pallet at the depot by U 0 .
Moreover, let the set V = {1, 2, . . . , n} index the vehicles available to service customers, and denote the cost per kilometre and the cost per hour expended by vehicle k ∈ V by C D k and C T k , respectively. Suppose the pallet loading capacity of vehicle k ∈ V is C P k , and that the expected travel distance and travel time from customer i ∈ L to customer j ∈ L are denoted by d ij and t ij , respectively.

Model variables
Define the binary decision variable and capture the number of pallets delivered to customer j ∈ L \ {m + 1} (or being loaded at the depot in the case of j = 0) by vehicle k ∈ V in an integer decision variable p jk . This variable allows for split deliveries, as the demand experienced by some customers may not be serviceable by a single vehicle.
While it can, in fact, be derived from the decision variables and model parameters, an auxiliary variable s jk is defined in order to simplify the model constraints, capturing the expected time (measured in hours) at which vehicle k ∈ V will begin its service at customer j ∈ L. A vehicle may arrive at a customer before the service window of that customer starts, but can only start to unload (commence service) within the window. Note that s 0k represents the time at which vehicle k ∈ V starts loading at the depot and that s m+1,k denotes the time at which the vehicle returns to the depot. A second auxiliary variable z ijk is required to enable the objective function to be written in a linear form. This non-negative variable either represents the difference between the times when vehicle k ∈ V begins service at customers i ∈ L and j ∈ L if vehicle k indeed travels from customer i to customer j, or is considered meaningless otherwise.

Model constraints
A number of constraints are imposed in the model so as to ensure the practical feasibility of solutions. In order to ensure that the demand of customer j ∈ L \ {0, m + 1} is met, the constraint set is enforced. Moreover, to allow later constraints to accommodate the time it takes to load each vehicle at the depot, the number of pallets loaded into vehicle k at the depot is required to match the sum total of all pallets being delivered by the vehicle. The constraint enforces this. A vehicle, of course, cannot service a customer if it does not visit this customer. In the constraint set the right-hand side equals zero if vehicle k is not assigned to service customer j. The constraint set therefore prevents the vehicle from servicing that customer. Each vehicle should also depart from the depot exactly once if it is assigned to service customers at all, and it may not depart with a pallet load greater than the vehicle's capacity. The constraint set ensures this. The depot itself should be excluded from the set of customers in this case. Similarly, each vehicle is required to return to the depot exactly once after having departed from the depot. Applying the constraint set j∈L\{0,m+1} x 0jk − i∈L\{0,m+1} enforces this, while the constraint set j∈L\{0,m+1} allows vehicles to depart from the depot only once each. The constraint sets and j∈L x m+1,jk = 0, k ∈ V are furthermore included in a bid to reduce the model solution time. These constraints are not necessary to find feasible solutions, as x i0k and x m+1,jk are meaningless in the final output. Informal empirical experimentation, however, showed that the inclusion of constraints (7)-(8) leads to a significant decrease in the time required to solve the model.
By ensuring that the difference between the sum total of all arrivals at customer j and the sum total of all departures from customer j remains zero for each vehicle k, the constraint set limits each vehicle to depart from a given customer the same number of times that it arrives there (either once or not at all). Constraint set (9) also prevents vehicles from "departing" and subsequently "arriving" at the same customer (thus avoiding loops), as the exclusion of customer j from the first sum means that the constraint set is not satisfied if x jjk = 1.
With constraints in place ensuring that demand is met and that vehicles are assigned routes that begin and end at the depot, the following constraints ensure solution feasibility with respect to the time windows. The auxiliary variable s jk , representing the time at which vehicle k starts servicing customer j, has to be constrained to within the given customer's time window. The constraint set ensures that if a vehicle arrives at a customer before the start of that customer's time window, it can only start unloading at (i.e., servicing) the customer once its service time window begins. Each side of the constraint set (10) could have been multiplied by i∈L x ijk to render s jk = 0 if vehicle k does not visit customer j. This would, however, have introduced non-linearity into the model. The value of the variable s jk can instead simply be considered meaningless if vehicle k does not visit customer j.
It should further be ensured that every vehicle arrives at each customer to which it has been assigned service before the end of the customer's time window. While constraint set (10) ensures that the service only starts within the customer's time window, the constraint set ensures that the expected time at which service starts at customer i, plus the expected unloading time at customer i, plus the expected travel time from customer i to customer j is no more than the time at which service at customer j starts. The right-hand side of the above constraint set ensures that the constraint is only enforced if vehicle k travels from customer i to customer j. This requires the constant The constraint set is required to allow the objective function to be represented linearly. In this constraint set, the value z ijk must be at least the difference between between s jk and s ik (the expected times at which vehicle k begins service at customers j and i, respectively) if the vehicle travels from customer i to customer j (i.e., if x ijk = 1). Recall that z ijk is non-negative.
where the maximum and minimum are taken over all customers j ∈ L. The fact that z ijk will be included in the objective function, which has to be minimised, means that the variable itself will be minimised to zero if the value of x ijk is zero.
It is implicitly assumed that multiple vehicles can unload simultaneously at the same customer. For all time-based parameters and variables, it is possible to consider a period longer than a day by extending W S j and W E j for all appropriate customers to the desired number of hours after the reference time 0. There may, however, only be one time window per customer for the entire period under consideration.

Objective function
The objective pursued in the model is to This cost function takes into account the distance over which each vehicle travels and the time that it spends waiting and unloading at customers, as well as the time it spends travelling between customers.
If vehicle k does not travel from customer i to customer j (i.e., if x ijk = 0), there is no contribution the objective function associated with that vehicle. If, however, vehicle k does indeed travel from customer i to customer j (i.e., if x ijk = 1 and z ijk > 0), the expected cost of both distance and time are taken into consideration. The cost of travelling is accounted for by multiplying C D k (the cost per kilometre of vehicle k) by d ij (the expected travel distance from customer i to customer j). The cost of time expenditure is accounted for by multiplying C T k (the cost per hour of vehicle k) by z ijk (the expected time difference between commencing service at customer j and commencing service at customer i).

Model implementation and verification
The model derived in §3.1 was implemented in IBM ILOG's CPLEX Optimisation Studio 12.10 (CPLEX), which is a state-of-the-art software suite for solving mixed-integer and linear programming problems. After constructing an initial solution by root node processing, CPLEX applies the branch-and-cut method to find an optimal solution to a mixed integer programming problem instance.
We subsequently verified our model implementation by applying the model to eighty small, randomly generated hypothetical test problem instances. These test instances were generated by sampling methodically randomised input parameters. Any set of input data for which CPLEX found no feasible solution was discarded upon which the problem instance was randomised again.
The following method was utilised during this randomisation process: 1. Three vehicle types were conceived and assigned the pallet capacities 16, 22, and 30, respectively. 2. Each vehicle type was assigned a random distance cost (a real number between 0.7 and 1.5) a random time cost (a real number between 7 and 15). 3. A predetermined number of vehicles were created, each assigned properties (cost-per-km, cost-per-hour, and pallet load capacity) matching a random vehicle type. 4. "Customers" 0 and m + 1 were both assigned coordinates (0, 0) in the Cartesian plane. 5. Each customer was assigned random coordinates (real numbers between −100 and 100) in the Cartesian plane.
6. The distances between locations were calculated as the Euclidean distances between the coordinates of customers in the Cartesian plane. 7. The travel times between customers were calculated as the distance divided by 80. 8. The sum of all vehicle capacities was divided by the number of customers. Each customer was assigned a random integer demand between 1 and this number. 9. With a probability N vehicles /(2N customers ), customers were assigned 1-hour delivery windows, the starts of which were integers between 5 and 23. Here and in the remainder of the paper N x denotes "the number of x." If, however, customer j was not assigned a delivery time window, its window start and end times, W S j and W E j , were set to 0 and 24, respectively. 10. Every customer was assigned a random per-pallet unload time between 0.01667 and 0.1667, corresponding to 1 and 10 minutes, respectively.
For each test problem instance generated as described above, the model solution output data were inspected and compared with the input parameters so as to verify solution feasibility. In particular, the output data aspects that were inspected for validity included the following: • Each vehicle had to be assigned a circular route, departing first from the depot ("customer" 0), visiting every customer that it had to service exactly once along its route, and finally had to return to the depot ("customer" m + 1) again. • The sum of all delivery quantities assigned to a vehicle had to be no more than its capacity. • The sum of all delivery quantities received by a customer had to be equal to the demand of that customer. • The service start time at any given customer had to be within the delivery window of that customer. • The service start time of any customer had to be at least as large as the service start times of all preceding customers visited along the vehicle's route.
Upon inspection, it was found that every solution returned by the model implemented in CPLEX satisfied all of the above criteria, and so the model was considered verified in the program testing with test data sense, as recommended by Kendall and Kendall [17]. Many of the parameter values selected for the above verification process were arbitrary, serving merely as a means to verify that the solutions returned by the model implementation indeed satisfy all the model constraints. For example, the probability of customers having delivery time windows associated with them could have been raised, but this would increase the portion of random problem instances that have no feasible solutions.
Some of the numbers mentioned above were, however, not chosen completely arbitrarily. For example, the impact of time and distance on the objective function was balanced, with roughly the same order of magnitude for the final values. Moreover, the range of possible vehicle capacity values was based on SPAR Western Cape DC's actual largest and smallest vehicle capacities. The DC's customers who actually have negotiated time windows associated with them also all have one-hour time windows. Finally, stores' per-pallet unload times often vary between one and ten minutes.

Model time complexity
In order to quantify the correlation between the problem instance dimensions and the required corresponding solution time, batches of test problem instances were generated randomly and solved. Each batch consisted of ten randomised problem instances, and the solution times were recorded. The numbers of customers and vehicles varied between batches, but were kept fixed within each batch. In particular, the number of customers varied from six to fifteen, while the number of vehicles was taken as three or five. The mean solution times associated with these batches are illustrated graphically in Figure 3.  Note that the vertical axis of the graph in Figure 3 is equipped with a logarithmic scale. The solution time durations therefore appear to follow an approximately linear upward trend, which very clearly demonstrates that the computational complexity of the model grows exponentially as the number of customers increases, for a fixed number of vehicles, while the computational complexity of the model also increases as the number of vehicles increases. Hence an exact solution of the model would not seem to be scalable to real-world problem instances which easily involve 300 customers and as many as 40 vehicles. This demonstrates the need for a metaheuristic solution approach that can provide reasonable approximate model solutions in polynomial time.

Approximate model solution procedure
As motivated in §2, we elected to solve the realistically sized instances of the CVRPTW model of the previous section by means of the state-of-the-art HGSADC algorithm. The algorithm can reportedly handle very large and/or very constrained problems exhibiting a variety of attributes without requiring significant modifications. Vidal et al. [39], in fact, used the algorithm to solve CVRPTW instances exhibiting both multi-depot and multi-period attributes. This section is devoted to a description of the working of the metaheuristic in §4.1, an overview of our implementation of the algorithm in §4.2 and a discussion in §4.3 on how this implementation was verified.

Working of the metaheuristic
The high-level working of the HGSADC algorithm is summarised in pseudo-code form in Algorithm 1. The algorithm employs a number of parameters, including an allowable number ItNI of non-improving iterations, the maximum run time Tmax allocated, the probability Prep that an infeasible solution is repaired, an allowable number Itdiv of non-improving iterations that may elapse before diversification is initiated, and an allowable number Itdec of iterations before decomposition is performed. Create an offspring C via crossover with P 1 and P 2 ;

5
Educate C by applying a local search procedure; Insert C into infeasible sub-population; 8 Repair with probability Prep; 9 if C feasible then 10 Insert C into feasible sub-population; 11 if maximum sub-population size reached then 12 Select survivors and cull sub-population; 13 if best solution not improved for Itdiv iterations then Use HGSADC to solve each sub-problem; 20 Reconstitute three solutions, and insert them into the population; 21 Return best feasible solution; The algorithm allows for the consideration of "solutions" that are infeasible with respect to load, duration and time window constraints. Such infeasible solutions are, however, penalised in the objective function. Load and duration penalties are increasing functions of the cargo a vehicle carries or time it drives beyond a predetermined maximum duration. Time windows are penalised based on the amount of "time warp" a vehicle has to perform. Time warp is the opposite of waiting time. If a vehicle arrives at a customer after the end of the customer's service window, it has to "warp" back to the end of the service window. The time between when the vehicle arrives at the customer and the end of the customer's service window determines the penalty in this case. For the remainder of that route, the vehicle is assumed to have arrived at the customer at the end of the service window.
Each candidate solution is encoded as a set of chromosomes. The number of chromosomes depends on the attributes of the CVRPTW instance under consideration. Vidal et al. [38,39] included three chromosomes, namely period, depot and giant tour chromosomes. The period chromosome contains a list of integers for each customer, representing time periods during which the customers are serviced. The depot chromosome contains an integer for each customer representing the depot from which that customer is serviced. Further secondary chromosomes may be added to these in order to represent more attributes, or these chromosomes can be removed if they are not necessary; the algorithm is flexible in this respect.
The giant tour chromosome is the most consequential. All the other chromosomes may be reconstructed from this one. For each combination of period, depot and other problem attributes, the giant tour chromosome contains a sequence of values representing the order of customer visitation. Individual tours contain customers serviced by multiple vehicles, with no indication of which parts of these tours are serviced by which individual vehicles.
In order to evaluate the cost of a candidate solution, the tour has to be partitioned into routes. A split algorithm was proposed by Vidal et al. [40] for this purpose (i.e., to find optimal routes by splitting a giant tour into sub-tours for the available vehicles).
The An m × n matrix is maintained whose entry in row k and column t is denoted by potential[k, t] and represents the expected cost of servicing all customers from the start of a tour until customer t has been serviced, if k vehicles are utilised. The function is used to add the cost cost(i, j) of travelling from customer i to customer j to the potential cost array only if this does not violate the capacity constraint of the vehicle in question. Otherwise, the cost is considered infinite, which leads the algorithm to assign the next stop in the tour to another vehicle. The quantity sumLoad [x] in the function f (i, j) above denotes the vehicle cargo after having visited customer x.
Another m × n matrix is maintained whose entry in row k and column i is denoted by pred[k, i] and represents the customer at the end of the sub-tour immediately preceding the one in which customer i is included, when k vehicles are used. A domination function is employed to compare the cost of moving to the next customer from two different customers and returns True only if customer i is better than customer j as a predecessor.
Once the split algorithm has been executed, the optimal number k of vehicles can be found by comparing the cost of service values found in potential[k, N customers − 1]. Once In order to evaluate any individual solution P , the HGSADC algorithm employs the biased fitness measure biasedFitness(P ) = costRank(P ) + diversityRank(P ) × 1 − N elite N individuals which accounts for both the solution cost and the diversity contribution of P to the population of solutions. After splitting and evaluating all solutions in the population, the individuals within this population are ranked according to their cost. Solutions are also assigned diversity rankings based on the well-known Hamming distance between the secondary chromosomes (the period and depot chromosomes). While the secondary chromosomes are not binary, the distance between two solutions is considered to increase by one for each corresponding element in the chromosomes that is not identical (e.g., each customer that is not served during the same set of periods). Each solution is assigned a distance score, after which the diversity ranks are determined, where the solution incurring the largest distance receives the best rank. The impact of the diversity rank on the biased fitness is based on both the number N elite of elite individuals to keep when culling a population and the number N individuals of individuals currently in the population. This impact is always less than the impact of the cost rank, causing the lowest-cost solutions always to be retained from one generation to the next.
The algorithm initialises its solution population by randomly generating a number of individuals equal to four times a pre-specified minimum sub-population size. The algorithm Algorithm 3: Pseudo-code description of the PIX operator of Vidal et al. [38]. 1 Step 0: Inheritance Rule; 2 td ← number of depot and period combinations; 3 Pick two random numbers between 0 and td according to a uniform distribution. Let n 1 and n 2 be respectively the smallest and largest of these numbers; 4 Randomly select n 1 (depot, period) couples to form the set Λ 1 ; 5 Randomly select n 2 − n 1 remaining couples to form the set Λ 2 ; 6 The remaining td − n 2 couples make up the set Λ mix ; 7 Step 1: Inherit data from P 1 ; 8 for each (depot, period) (d, p), belonging to sets Λ 1 and Λ mix do 9 Λ 1 : Copy the sequence of customer visits from V d,p (P 1 ) to V d,p (C); 10 Λ mix : Randomly select two chromosome-cutting points α and β (according to a uniform distribution), then copy the α to β sub-sequence of V d,p (P 1 ) to V d,p (C); 11 Step 2: Inherit data from P 2 ; 12 for each (depot, period) (d, p) ∈ Λ 2 ∪ Λ mix selected in random order do 13 Consider each customer visit i in V d,p (P 2 ) and copy it to the end of V d,p (C) if the following two conditions are met; 14 1) The depot choice δ i (C) for the child C is equal to d or undefined (meaning no visit to i has been copied to C yet); 15 2) One visit pattern of customer i, at least, contains the set π i (C) ∪ p of visit periods; 16 Step 3: Complete customer services; 17 Perform the Split algorithm and extract the routes for each (depot, period) pair; 18 if the service-frequency requirements are satisfied for all customers then 19 Stop; 20 else 21 while customers with unsatisfied service-frequency requirements exist do 22 Randomly select a customer i for which service-frequency requirements are not satisfied; 23 Let F be the set of admissible (depot, period) combinations (d, p) with respect to its pattern list L i and the visits already included in C. Let ψ(i, d, p) be the minimum penalised cost for the insertion of customer i from depot d in period p; 24 Insert i into (o * , l * ) = arg min (d,p)∈F {ψ(i, d, p)}; performs several steps during each iteration. The first step involves selecting two parent solutions. Two pairs of solutions are chosen randomly from the pool of solutions currently within both the feasible set and the infeasible set. The solutions in each pair are compared, and the one with the better biased fitness value is selected as a parent.
A child solution is created by applying a crossover operator to the giant tour chromosome with the aid of the information contained in the secondary chromosomes. The 2012 paper by Vidal et al. [38] contains a new periodic crossover with insertions (PIX) operator, for which a pseudo-code description is provided in Algorithm 3. This operator pursues a balance between search space exploration and attaining improvements on existing solutions. An offspring solution inheriting similar amounts of genetic material from both parents will be substantially different from both of them, while one that inherits mostly from one parent will differ only slightly from that parent. The PIX procedure employs randomness to determine how much material to inherit from each parent instead of following predefined rules by which offspring inherit a fixed amount of material from each parent.
After a child solution has been generated, there is a chance to educate it, based on an education probability parameter. The education procedure is what differentiates the HGSADC algorithm from standard genetic algorithms, as the procedure applies a limited local search to the solutions generated in order to find quick and easy improvements. Education makes use of two local search procedures: Route improvement (RI) and pattern improvement (PI). These procedures are applied in the order RI, PI, RI.
The RI procedure seeks to optimise each tour separately. It considers various insertion, swap and 2-opt moves [7]. For each vertex in a route, a granularity threshold is employed to restrict the search neighbourhood to the closest vertices [34]. The procedure randomly iterates over each vertex and each of its neighbours, attempting to apply the moves in a random order. The procedure terminates when all possible moves have been attempted in respect of all neighbours of the current vertex without uncovering an improving move.
The PI procedure attempts to move customers between depots and periods. Iterating over each customer in a random order, it computes the minimum cost of satisfying the service requirements of the customer from the current depot and within the current period, and again with each other possible combination of depot and period. If there is another combination that yields a lower minimum cost of serving the customer, the customer is inserted in the position that achieves the minimum cost.
After having implemented the RI, PI, RI education procedure, the child solution may or may not be feasible. It is inserted into the appropriate sub-population. Based on a repair probability, infeasible solutions may be repaired. In order to repair an individual, the penalty parameters are multiplied by 10 and the education procedure is performed again. If the individual is still infeasible, the penalty parameters are multiplied by 100 and the procedure is repeated. If a feasible individual is produced by the repair procedure, this individual is added to the feasible population without modifying or removing the original infeasible individual from the infeasible sub-population.
When either sub-population reaches the maximum allowed size, it is culled down to the minimum population size. All solutions are ranked according to their biased fitness values based on both cost and diversity contributions. A list of clones is compiled, where individuals are considered clones if either their depot and period chromosomes or their solution costs match. The clone with the worst biased fitness value is removed until there are either no more clones, or until the minimum population size is reached. If there are no clones, or all clones have been removed, the remaining individuals with the worst biased fitness values are repeatedly removed until the minimum population size is reached.
Every 100 iterations, the capacity and duration penalty parameters are adjusted, based on the proportions of naturally-feasible individuals with respect to vehicle capacity and route duration over the last 100 iterations. The parameters are adjusted slightly in order to achieve a target proportion of naturally-feasible individuals.
If more than Itdiv iterations have been performed without encountering any improving solution, both sub-populations are culled to one third of the normal minimum size. A number of new solutions equal to four times the minimum population size are then generated randomly and added to the appropriate sub-populations. This introduces a large number of novel solutions to the populations.
The main contribution of the 2013 paper by Vidal et al. [39] is the algorithm's decomposition procedure. It is performed once every Itdec iterations, and only on problem instances containing more than 120 customers. A random individual is selected from the 25% best feasible solutions and split into smaller parts, thus inducing sub-problems. Each tour in this individual's giant tour chromosome is iterated across sequentially, generating a new set of tours. Every time the set reaches more than 120 customers, or all tours have been visited, a new sub-problem is created from this set of tours and solved by performing another instance of the HGSADC algorithm, with the maximum non-improving iterations parameter set to Itdec/2. Once every sub-problem has been solved, the top three solutions from each sub-problem are combined into three new solutions which are then introduced into the appropriate main populations.
Finally, once the algorithm has performed ItNI iterations during which no new best solution update has occurred or if the specified maximum run time is reached, it returns the best solution encountered during the entire search.

Metaheuristic implementation
We implemented the HGSADC algorithm described in the previous section in Python 3.8, using the CPython interpreter. A handful of publicly available libraries were installed on the interpreter, including OpenPyXl and NumPy which provide useful tools for scientific computing and provided us with the capability of reading input from and writing output to Microsoft Excel files. The sequence of steps followed within each algorithmic iteration are almost identical to those described in Algorithm 1, with the addition of programming logic facilitating a comparison between generated solutions and the currently best-known solution. We adopted an object-oriented approach in order to achieve modularity and clarity of implementation.
We omitted the chromosomes accounting for the multi-period and multi-depot attributes in our implementation, replacing them with one capable of indicating which type of vehicles from the DC's heterogeneous delivery fleet should service each customer. The tours within the giant tour chromosome were therefore appropriately indexed by vehicle types, rather than by depots and periods. The capacity, cost per kilometre and cost per hour incurred by any vehicle that services customers is determined by its type. The implementation also accounts for the number of vehicles of each type that are available for service according to the input data as well as for the possibility of split deliveries, which often occur at the SPAR Western Cape DC. This occurs when a customer exhibits a larger demand than can be satisfied by a single vehicle, in which case multiple vehicles are required to service the customer.
The customer completion procedure of Algorithm 3 was found to be insufficient for our specific context, as it was designed only to check that customers are visited during as many periods as required by the customer. A wholly new customer completion procedure, described in pseudo-code form in Algorithm 4, was therefore designed and implemented to ensure that all customers' demands are met exactly without disrupting the relative positions of split deliveries within their tours.
The split algorithm of Vidal [36] proved the most significant challenge to implement. Fortunately, Vidal made his C++ implementation of the algorithm available on GitHub [37]. This, together with a slightly different implementation by Juppe [16], provided enough guidance for us to implement the algorithm anew. In our implementation, the algorithm takes as input single tours and vehicle type information, and then returns as output lists of vehicle routes, where each route is a sequence of customers to be visited by a single vehicle.
Another aspect of note about our implementation is the way in which departure times are calculated. Although Vidal et al. [39] described how time-warps may be used to penalise solutions that involve vehicles arriving at customers after the end of their time windows, they did not describe how to determine the corresponding vehicle departure times. Although algorithms are described in the literature for determining optimal departure times (see, for example, [18]), these algorithms were deemed computationally too expensive to integrate with the HGSADC algorithm. A greedy heuristic was therefore designed and implemented in which vehicles are assumed to service customers at the latest possible time. Each vehicle's departure time from the depot is set such that it is expected to arrive at the first customer with a time window along its route at the end time of this customer's time window. This may result in a number of slightly sub-optimal departure times, and would require SPAR schedulers who use the model to determine the departure times for vehicles along routes without time windows themselves.

Implementation verification
A careful verification of the metaheuristic implementation described in the previous section is important, as it provides an indication of whether the solutions produced by the metaheuristic are of a sufficient quality to be recommended to the SPAR Western Cape DC. The exact CPLEX model solution implementation of §3.2 was employed as reference point in order to verify the aforementioned implementation of the metaheuristic.
Recall that during the verification of the exact model solution procedure in §3.2, randomised test instances involving small numbers of customers and vehicles were generated. The input and output data of these test instances were recorded. A script was implemented that would read these recorded data and sequentially execute the metaheuristic implementation of §4.2 for each instance. The metaheuristic was subjected to a maximum allowable number of 5 000 non-improving iterations (no maximum permissible run time was imposed).
The metaheuristic was verified by comparing the cost of its best solution for each test instance with the cost of an exact solution to the same instance. Occasionally the two solution approaches resulted in slightly different evaluations of the same solution. This occurred in only one of sixty of the random test instances involving 6-12 customers, but in all twenty test instances involving 15 customers, and is due to the fact that the mathematical model of §3.1 optimises each vehicle's departure time (when solved exactly). The metaheuristic of §4.1, however, adopts a greedy heuristic to guesstimate a good departure time of the vehicles, as described in the previous section. This may possibly degrade the (accidental) optimality of some solutions returned by the metaheuristic, but such degradations were observed to have a negligible impact in the case of the large problem instances for which it was intended.
In order to facilitate a fairer comparison between the solution costs for verification purposes, costs associated with the exact solutions to the mathematical model were subsequently also evaluated by means of the heuristic employed in the metaheuristic implementation of §4.2. Table 4 contains the mean differences in cost between the solutions returned by the metaheuristic and the corresponding exact solutions for the randomised test instances of §3.2. The mean difference between the solution costs and the number of test instances for which the metaheuristic produced identical solutions to those returned by CPLEX are also shown in the table. Solutions were considered identical for practical purposes when the difference in cost, as evaluated by the metaheuristic, was less than 0.001%. The difference in cost between solutions was calculated as (Meta Cost -Exact Cost)/(Exact Cost).  The metaheuristic obtained 59 solutions that are identical to those returned by CPLEX out of the 80 small random test instances considered, with an overall average optimality gap of 0.56%. The metaheuristic is therefore considered to have been verified successfully. Figure 4 contains a plot of the mean solution times incurred by the metaheuristic for these test instances on the same processor as that used in §3.2. The mean solution times incurred by CPLEX are also plotted by means of dashed lines in the same figure for comparison purposes. When comparing the mean solution times of the two solution implementations, it is worth noting that CPLEX utilises multiprocessing, while the metaheuristic implementation only utilises a single thread. CPLEX can therefore take advantage of all the available processing power, while the metaheuristic can only use a maximum of a quarter of the processing power available on the quad-core CPU upon which the test instances were solved. Additionally, the exact solution times of the instances fluctuated to a far greater degree than did the approximate solution times. CPLEX ran so long for some of the larger test instances that the branch-and-cut procedure either exhausted all the computer memory or was interrupted by Windows updates after many hours. Such instances were re-randomised and re-run in an attempt at ensuring a fairer comparison. Although the computation time required by the metaheuristic (to find approximate model solutions) would seem large in comparison with that required by CPLEX (to find exact model solutions) in Figure 4, the slope of increase in the former times (the solid line graphs in the figure) are not as steep as that of the latter times (the broken line graphs in the figure). Moreover, it is not visible in Figure 4, but the solution times required by CPLEX rise substantially as the number of customers increases beyond fifteen. In fact, for eighteen customers, CPLEX could not obtain exact solutions within 100 000 seconds (27.78 hours). In contrast, the solution times incurred by the metaheuristic are considered to be acceptable, even for relatively large problem instances.

Case study
The efficacy of the mathematical model of §3 and the approximate solution methodology of §4 are evaluated in this section in the form of a case study involving real data obtained from the SPAR Western Cape DC for three work days during 2019. This case study allows for a comparative evaluation of the vehicle routes recommended by the model with the status quo. The section opens in §5.1 with a background discussion on the data pertaining to the case study and a discussion in §5.2 on the assumptions made in order to implement the model of §3 and the metaheuristic of §4 in the context of real data. The results returned by the metaheuristic for the three days in question are presented and discussed in some detail in §5.3.

Background information and input data
As mentioned in the introduction, the SPAR Western Cape DC serves more than 300 customers in the Western Cape and Northern Cape provinces of South Africa. Human schedulers construct daily vehicle routing schedules by hand, aiming to satisfy all customer demand on any given day. These schedulers utilise a predetermined delivery schedule to guide their routing decisions. This schedule is based on one-hour delivery windows starting at different times of the day, which are distinct for different days of the week, because some customers only specify delivery time windows on certain days of the week. Past experience of the schedulers in conjunction with this schedule allows for the construction of high-quality delivery vehicle routes that enable the DC to operate effectively in a relatively cost-efficient manner.
The DC has more than 40 vehicles in its active delivery fleet, most of which consist of interchangeable horses and trailers of varying capacities, but the DC also temporarily hires third-party vehicles and drivers, as required, in order to facilitate a flexible service capacity. Vehicles that return to the depot upon completing delivery routes may be assigned further routes to service on the same day.
All deliveries are scheduled on a central Google Sheets document visible to multiple parties within the DC. Once a vehicle has departed from the DC, information about the vehicle's route is removed from the live document and stored in an archive. The 2019 archive shows that the number of deliveries scheduled per day varied between 2 and 369, with a mean of 169 (190 when excluding Sundays). These numbers allowed us to identify three days -a relatively quiet day, an average day, and a busy day -on which to base the case study of this section. Information pertaining to the numbers of deliveries, delivery locations, vehicle routes, in-house and third-party vehicles, as well as demand volumes, experienced on these three days are shown in Table 5.  The archived data for the three days listed in Table 5 were formatted into the appropriate format required by our metaheuristic implementation in order to facilitate an evaluation of the performance of the model of §3 and the metaheuristic of §4 in the context of real data.
For the sake of transparency and reproducibility, these data are available electronically [30].
The vehicle type information summarised in Table 6

Simplifying assumptions
The model of §3 does not, of course, represent all the complexities of the real-world situation exactly, and so our most important model simplifications are addressed in this section in order to facilitate a meaningful evaluation of its efficacy.
Each of the three problem instances (days) considered in the case study covers a period of time large enough for vehicles to service multiple routes upon having returned to the DC in order to load new pallets and having been assigned new drivers. These nuances were excluded from the scope of the small test problem instances considered in §3.2 and §4.3, but are accounted for during the evaluation process documented here.
As mentioned, the DC has the option of hiring additional vehicles in order to complement its own fleet, as necessary. The linear-time unlimited-fleet split algorithm [36] was therefore used when solving the three problem instances considered in this section, as opposed to the polynomial-time limited-fleet split algorithm used during the verification performed in §4.3. This achieved a significant decrease in the required solution time for the larger two problem instances.
The evaluation of the cost associated with solutions returned by the metaheuristic was also configured to estimate the number of vehicles required by summing the expected durations of routes together with additional "wasted" time between routes, and dividing this sum by a standard number of hours of service per vehicle per day. Face validation of the cost estimation thus produced was performed by February [10], the transport controller at the SPAR Western Cape DC at the time of writing. He is responsible for research into the costing of the DC's deliveries and routing decisions, and thoroughly understands the method of cost estimation utilised at the DC. He confirmed that the cost estimations returned by the method described above were close enough to the estimates produced by the DC's ERP systems to facilitate meaningful comparisons.
In reality, hiring vehicles is often cheaper for the DC than using its own vehicles, largely due to how well the DC pays its drivers. It is, however, still preferable for the DC to use its own vehicles and drivers before hiring third-party resources. A cost multiplier was therefore applied to any vehicles used beyond the DC's own active delivery fleet. It is assumed that the first routes appearing in any recommended tour will be assigned to the DC's own vehicles, until all of these vehicles have been exhausted, upon which the cost multiplier is applied to later routes that exceed the total service time estimated to be available to the DC's own vehicle fleet.
The total loads of routes in the DC archive often appear to be greater than the capacities of the vehicles assigned to routes. This is feasible in reality if some of the pallets destined for customers along the vehicle's delivery route are not completely full. The employees loading the vehicles are able to consolidate these partially full pallets in a bid to reduce the total pallet count. In contrast, the model has no way of determining how full pallets are. It would not, however, be sensible, when evaluating the performance of the model, to assume that vehicles have soft load capacities. The demands of customers along such a route were instead reduced slightly until the vehicle's capacity was no longer exceeded. Care was taken not to reduce any customer's demand to zero, thus ensuring that any feasible solution would still visit all customers.
Sometimes distinct customers -such as a TOPS and a SUPERSPAR -would occupy the same physical location, but are nevertheless treated as separate customers by the DC's ERP systems. In such cases, the schedulers will, of course, schedule stops at these customers to occur immediately adjacent to one another along a delivery vehicle route. In order not to increase the problem instance dimensions unnecessarily, demand and consecutive stops for customers at the same location were merged when performing the case study of this section. The values in the Locations column in Table 5 represent these merged customer locations.
The DC's geo-tracking system was able to provide average unload times per pallet for most locations, but not all. For locations for which unload times were not available, unload times were assumed to be 5.5 minutes per pallet, on average.
We retrieved the expected travel distances and travel times between customers using the Bing Maps API [20]. Since no vehicle departure or arrival times were available prior to solving the model, the results requested from the API conform to the expected average travel times between the locations. In addition, 45 minutes of rest time were added to the durations of routes for every five hours of continuous driving. Eight hours of rest time were added to the duration of routes for every ten hours of total driving time required.
Since the case study of this section is aimed at ascertaining the potential cost savings of the DC when using the model proposed in this paper to guide its delivery vehicle routing decisions, the usual customer service time windows were ignored. As mentioned, these scheduled delivery windows merely serve to guide the human schedulers' routing decisions. We believe that greater cost improvements can be realised by relaxing these windows.

Numerical results
Despite employing the linear-time split algorithm, our metaheuristic implementation was found to execute slower than expected when solving the three case study problem instances. It created, educated and evaluated approximately 100 solutions every two minutes on a 3.41 GHz Intel i5-7500 processor with 8 GB of RAM running in a Windows 10 operating system. The model was initialised with parameter values specifying that at most 2 500 non-improving iterations and at most 7 200 seconds of run time should be allowed. The thresholds for performing a problem decomposition phase were reduced to 1 000 nonimproving iterations and 60 customers.
In order to take advantage of the DC schedulers' experience and knowledge, the algorithm was seeded with a solution derived from the DC's archived data in addition to the population of random solutions it generates upon initialisation. This seeded solution was used during parent selection and was iteratively improved alongside the novel randomly generated solutions.
The vehicle routes actually implemented at the DC during the high-demand day (26 November 2019) are tabulated in  in the appendix at the end of the paper. In total, there were 211 stops along 128 routes in this set of vehicle routes. For each route, the anonymised store names are provided in the order visited by the vehicle upon leaving the DC and before returning to the DC, together with the number of pallets delivered at each location along the route in brackets, as well as the total travel distance, travel time, and service cost associated with the route. The route costs tabulated are as they appear before having applied the multiplier for hired vehicles. The total estimated cost of implementing these routes was R 1 564 771.87.
The status quo solution in Tables 11-13 was included in the initial population of solutions with which the algorithm of §4 was initialised. The vehicle routes appearing in the best solution returned upon termination of the algorithm for the high-demand day are summarised in Tables 7-9 in the same format as that used in the appendix. In total, there are 170 stops along 107 routes in the vehicle route set returned by the algorithm, and the total anticipated cost of implementing these routes is R 1 410 468.59. This represents a R 169 522.45 cost saving, which equates to 10.83%, over the route set actually implemented.
The structure of the solution in Tables 7-9 suggests that the DC should have relied more heavily on its larger vehicles to perform deliveries on the day under consideration, according to the metaheuristic. The consolidation of multiple smaller loads into larger loads allows for a large reduction in cost, as it reduces the frequency with which the depot must be visited. The solutions returned by the metaheuristic, in fact, exhibit an increasing reliance on the smaller vehicle types as demand increases and the cost of hiring additional vehicles becomes more significant.
The numbers of stops and routes, as well as the estimated total costs, of the status quo and the metaheuristic-recommended route sets are provided in Table 10 for all three days considered in this case study. As is clear from the difference columns of the table, the estimated costs associated with solutions returned by the metaheuristic of §4 are significantly better than those associated with the actual archived route sets implemented.    The relative percentage improvements for the different days of the case study make sense when compared with one another. The second of these days, 30 October 2019, was an average day, with demand levels similar to the vast majority of days in the year. The DC's schedulers are accustomed to dealing with this level of demand and their guiding time-window schedule is focused towards handling this type of day. The first day, 7 October 2019, however experienced a lower demand than usual. The schedulers' reliance on pre-determined delivery time windows may have impacted negatively on their overall routing efficiency on this day. Finally, 26 November 2019 was part of a higher-than-average demand period just before Black Friday. When faced with such a peak in demand, mistakes and inefficiencies are amplified. For example, one customer in that day's archive appears to have been served four times, with vehicles departing along routes consisting only of visiting this customer (at 9:55, 12:10, 20 :30, and 22:20). None of these vehicles carried full loads. It is not hard for the metaheuristic (with perfect information on that day's demand) to find better routing solutions.  It is acknowledged that the three problem instances of this case study were solved outside of the busy and chaotic context of managing operations at the DC. In reality, schedulers will not have perfect information on a day's demand levels. Moreover, the expected distances and durations of travel and service are only estimates, and there will invariably be other delays and complications beyond the scope of the model of §3 to cope with. It is therefore expected that a portion of the savings indicated in Table 10 will be lost to these real-world inefficiencies when applying the metaheuristic-recommended routing decisions in practice. Nevertheless, there is clearly a potential for the SPAR Western Cape DC to realise considerable cost savings if their schedulers were to utilise the model and accompanying metaheuristic of this paper to augment daily vehicle routing decisions.
Finally, the modelling approach is rather flexible and can be applied in a great variety of situations. It may, for example, be used to plan the servicing of known demand for a full day (or more), over the next twelve hours (the minimum time required at the DC before picking of customer orders may take place), or over a period of a few hours in advance (including only the demand serviced by one or two picking waves in the DC). Time windows can, furthermore, be customised independently for each customer, allowing for experimentation and further optimisation in cases where customers are not insistent on receiving service during specified periods of time. The maximum allowable number of non-improving iterations and the maximum run time can also be customised as required, allowing the metaheuristic to be left running for longer periods of time with a view to find better solutions, especially when exploring the search spaces of large problem instances. The algorithm can also be allowed to use an unlimited fleet of vehicles, or be forced to utilise only a limited fleet (which is only recommended for smaller problem instances, as it otherwise significantly slows down the algorithmic implementation). The costs of using different vehicle sizes, as well as the cost multipliers applied for hired vehicles, can be adjusted to influence the proportions of different vehicle sizes utilised.

Conclusion
In this paper, we formulated a CVRPTW model appropriate for use at the SPAR Western Cape DC, implemented this model in CPLEX, and verified the implementation in respect of a sample of small randomised test instances. As expected, however, the branch-and-cut method employed by CPLEX required solution times that increase exponentially as the problem instance dimensions increase linearly. This indicated the need for a metaheuristic approach towards solving realistically sized problem instances (approximately).
A state-of-the-art hybrid metaheuristic, called the HGSADC algorithm, was therefore implemented in Python with minor alterations (such as those applied to the customer completion procedure and the secondary chromosomes utilised). This metaheuristic implementation was again verified in the context of the same small test instances solved by CPLEX. While the metaheuristic also exhibited solution times that increase exponentially with linear increases in the model instance dimensions, the rate of increase was far smaller than that incurred by CPLEX.
A practical case study was finally performed by applying the HGSADC algorithm to historical routing data obtained from the SPAR Western Cape DC, pertaining to three workdays of representative dispatch workloads. After making a number of necessary simplifications, the costs of the historically implemented routes for these workdays were evaluated and compared with the costs associated with routing solutions recommended by the metaheuristic. Significant cost savings of the order of 5-10% were thus realised. While admittedly a simplification of reality, the cost estimation procedure employed during this comparison was face-validated by the transport controller at the DC, lending credence to the cost savings reported.
The potential benefits of utilising the metaheuristic to assist the dispatch department of the DC in routing decisions are compelling, but caution is advised in accepting the reported savings at face value. Many simplifications were necessary during the derivation of the CVRPTW model, including the implicit assumption that all input information pertaining to a problem instance is deterministic. Execution of the recommended routing decisions in the real world may therefore well result in slight degradations of the expected savings. The reported daily cost savings are nevertheless significant and warrant consideration of an integration of commercial vehicle routing decision support software into the ERP systems already in use at the DC instead of relying entirely on manual vehicle routing.