Enhancing heuristic bubble algorithm with simulated annealing

: In this study, a new way to improve the Heuristic Bubble Algorithm (HBA) is presented. HBA is a nature-inspired algorithm, which is a new approach to and initially implemented for, vehicle routing problems of pickup and delivery (VRPPD). Later, it was reinforced to solve other routing problems, such as vehicle routing problem with time windows (VRPTW), and vehicle routing problem with stochastic demands (VRPSD). HBA is a greedy algorithm. It will mostly find local optimal solutions. The proposed method is an improvement over HBA enabling it to reach the global minimum. It uses specialized simulated annealing methods in its operators. A well-known data-set is used to benchmark the proposed method. Better results over HBA and some best results in literature are recorded.


Introduction
Vehicle routing problem (VRP), or "Truck Dispatching Problem", is a deterministic polynomial-time hard (NP-hard) problem and can be formulated as a generalization of the "Travelling Salesman Problem" (TSP). TSP, although its beginning is unclear (Hoffman et al., 1986), was first publicized by Flood (1955).

ABOUT THE AUTHOR
Ali Gunes is a professor of Computer Engineering Department of Istanbul Aydin University. He is also the Head of Distance Education Research and Application Centre at Istanbul Aydin University. He currently teaches, consults and conducts research on computer software models and software development, application of information technologies on different sectors and new learning technologies in computer education. His minor is distance education, e-learning and transportation of educational media on different platforms. He has been studying computer software modelling and development, integrated education and technological developments on education for more than 30 years in the academic field. He has written more than 40 articles and more than 10 books and other academic works, on computer software, integrated software development and e-learning published in several journals and publishing companies.

PUBLIC INTEREST STATEMENT
The importance of transportation services for businesses is very great. Nowadays, the issue of transportation of goods and services contains many problems that enterprises need to resolve. In this regard two most important problems of transportation that can be done timely and at optimal cost. One of the most important issues in determining the cost of transportation is to allow goods to be made via the shortest possible route. Many studies have been done and several models and algorithms have been proposed to solve this problem. In this study, a new way to improve the Heuristic Bubble Algorithm (HBA) is presented. HBA is a new approach to and initially implemented for, vehicle routing problems of pickup and delivery. The proposed method is an improvement over HBA. A well-known data-set is used to benchmark the proposed method. Better results over HBA and some best results in literature are recorded.
The problem is formulated as such: • G: a Hamiltonian path or a complete graph, • V: set of vertices, • E: set of edges, • c ij : cost linked with each element inside E, is given. c ij is the cost to traverse from k ∈ V to m ∈ V (Zambito, 2006). VRP, on the other hand, is first presented by Dantzig and Ramser (1959) and is formulated as: • R: a set of routes, • V: that is served by a fleet of vehicles, • P: that is visiting a set of customers, • C: with some constraints. (Kumar & Panneerselvam, 2012). Thus, each of the routes inside a VRP solution becomes a TSP instance.
Simulated Annealing (SA) has been proved to solve VRPPD problems (Czech & Czarnas, 2002). In this study we extended E-HBA with SA and applied it to the entire schedule resulting from a single run of the HBA algorithm.
Section 2 of this paper discusses SA as it is applied to other problems. Section 3 describes the E-HBA. Section 4 presents the way SA is applied to E-HBA, as well as the results from Gehring and Homberger's (1999) data-sets. Section 5 is a computational study and Sections 6 and 7 end with sample results, conclusions and future directions.

Simulated annealing
SA was first formulated by Černý (1985) and Kirkpatrick, Gelatt, and Vecchi (1983), as a re-engineering of a Monte-Carlo Method (Metropolis-Hastings algorithm) (Metropolis, Rosenbluth, Rosenbluth, & Teller, 1953). In Kirkpatrick et al. (1983), interestingly, it was first implemented to solve a TSP problem using a probability provided by the Boltzmann-Gibbs distribution.
SA is a local search method capable of escaping the local minimum by taking a worse solution according to a probability (Lin, Yu, & Chou, 2009). It has been applied successfully to complicated combinatorial optimization problems (Abramson, 1991;Jayaraman & Ross, 2003;McKendall, Shang, & Kuppusamy, 2006).
The terminology is taken from metallurgical engineering. Annealing is the process of slow-cooling metals to produce better aligned crystallization. SA uses a similar method for getting close to the global minima, by slowly decreasing the heat. In each run, SA accepts a new best solution from the neighbouring solutions. If the solution is better, it is accepted as the new best solution. If it is worse, then according to a probability it may be accepted to continue with. By doing so, SA tries to escape from the local minima. This new solution becomes the starting point of the next run.
The probability calculations come from Boltzmann-Gibbs distribution (Černý, 1985;Kirkpatrick et al., 1983) which is expressed as where E is the energy of the state, k b is the Boltzmann constant and T is the heat. The ratio of two states in the distribution is called Boltzmann Factor and is given as which in turn simply can be expressed by where Δ is the change in energy.

Enhanced heuristic bubble algorithm
HBA is an iterative heuristic inspired by the division and union of water bubbles (Sakalli et al., 2013). Bubbles are the goods that are carried between distribution centres. During transportation, goods are combined and separated from each other to better serve scheduling needs. In each run, HBA tries to find the bubble that, according to the objective function, is the best route in the entire route set which can be constructed using the remaining orders. This structure is best suited for pickup and delivery problems.
HBA was enhanced with Split, Result Elimination (RE) and Swap Operators (Savran et al., 2014), hence the new name "Enhanced HBA" (E-HBA) Each addition provides a better enhancement over initial bubbles.
Next we will briefly describe the Merge Operator, as it is fundamental in understanding E-HBA. Because we have integrated SA as a RE operator, we will give a short explanation on this operator.

Merge Operator
The Merge Operator is the core of the HBA. Together with the Join Operator, it is the initial implementation of the algorithm. It searches the entire solution space for a single route and, according to the criteria presented to it, finds the best one. The length of the route is the main input provided for the operator. Owing to the combinatorial nature of the problem, small lengths are usually preferred. When a new route is found, unscheduled orders are processed to find the next route. During this process, the first route found is the overall "best route". Later routes are called "worse routes"; because they are constructed after better routes preceding them, and they have not been proved to be the best in the solution space. If the first route were to be another route, later routes could be constructed differently. In this, the merge is blind; it just tries to find the best route from remaining orders, without thinking about later routes.

Result elimination
As we have indicated, the routes that are generated later in the merge process are called the "worse routes". The last route is called the "worst route". What RE does is to eliminate worse routes by combining them with better routes; if possible, combining them with the best route. When doing this, parts of the routes can be replaced. RE tries all the orders in the worse route one by one, and adds them to the better route. If any of the orders added to the better route makes the schedule better, it accepts the change and removes the order from the worse route and adds it to the better route. Thus, the worse route either diminishes completely or shrinks. During these operations, a couple of options manage how the RE is executed. For example, index selection is the process of finding the index on the best route to insert the order from the worse route. FIRST_FOUND method inserts the order into the first index found. This is a search process that tries to find indexes that do not violate constraints and makes the schedule better. FIRST_FOUND is very fast but cannot guarantee the best solution. On the other hand, the BEST_FOUND method considers all the indexes and inserts the order into the best index that provides maximum gain.

Applying simulated annealing to E-HBA
SA is applied to E-HBA as a RE Operator. As indicated, RE always tries to eliminate the worst routes, which are the later routes on the schedule. This is because when the E-HBA tries to construct the routes, it tries to find the best available from all the iterations, and so the remaining routes become worse according to the objective function, which is a simple distance function in our case.
Previous RE Operators were blind even if a worse objective would provide a better objective later in the iteration. In this situation, the new RE Operator, which is designed according to SA principles, can accept a worse objective according to the SA probability. This method, after applying it to the resulting schedule, can provide a better solution to the entire schedule, which we have shown below using Gehring and Homberger's (1999) data-set (referred to as GH in the coming sections). Appendix A gives the entire schedule for two of the best results we found (C1_10_2 and C1_10_6, respectively). However, we did not put all of the result schedules for the sake of brevity.
SA has a random change principle which takes a new solution from the result neighbourhood. Here our next neighbour is taken from the output of RE. In previous operators, a worse route could not be accepted, but this brand new operator can accept a worse route that could lead to a better solution when the heat goes to the minimum.
Basic algorithm steps are as follows (see Appendix B for the terminology); (1) Initially run a simple clustering algorithm (K-Means) to generate stop clusters.
(2) Set Merge Join distance (MJD) to a minimum value such that only neighbours in the same clusters are selected.
(5) If enough or no changes continue with 6, else go to 2.
(6) Now we have a reasonable schedule at hand.
(a) Select two routes according to the RE selection rules, (b) Select some stops from the first route according to RE line selection rules, remove and try to insert them to the other route, (c) If constraints are not violated, continue else go to 8-a.
(d) If new distance is better then accept and go to 8-a.
(e) If the new distance is not better, calculate SA principles and decrease the heat if probability allows the new state.
(f) If the heat reaches a minimum, go to 9.
(g) Else go to 8-a.
(9) If enough or no changes in the routes, change RE settings, continue with 8, if no more settings could be changed exit.

Computational study
After initial experimentations, we concluded that when enabling each of the features of simulated annealing, it was best just trying with the minimum setting. Thus, our SA engine only displays two settings to be arranged.
The two controls both arrange the same setting for the number of runs needed. If the "estimated number of runs" control is changed, then "Initial Temperature" also changes. In simulated annealing, it is crucial to set the Initial Temperature and the minimum temperature appropriately. They must be arranged according to the problem at hand. "Initial Temperature" is the heat to be set for the system. Other settings for SA are taken as shown below: is the heat changer in the SA system. If a new result is found, then the heat is multiplied by this value and the system continues to cool. is the minimum temperature we allow. "Default Heat" is the initial heat if no settings are given. These values are constant in our system. The only thing we change is the Initial Temperature and with it, the number of runs. The number of runs and the heat are calculated as follows in pseudo code, respectively; Since we have previous results, we can deduce appropriate values for SA to continue. What we did in each step of the algorithms is as follows: • If the sixth step is executed for the first time, initial temperature is set to 6. In turn sets the number of run to around 18,000. With these values we are trying to give SA enough room and length so that large gains can be obtained. With the first run, we know that the system is still hot and SA must know this.
• After each subsequent run, this value is changed until we arrive at a good result or no other changes take place.
If we give initial heat larger values than 5, then SA accepts too many poor results so that the overall distance increases and little or no gain is observed.
Our algorithm is coded in C# programming language on an Intel Core i7 2.4-GHz computer. We performed around 30 runs for each instance.

Sample results
Here we will give short descriptions on how we arranged our parameters to obtain the results. Also, the proposed method has many options that can be tweaked to make the compromises among time, distance or objective function at hand.
We take the routes as Hamiltonian paths; that is, the vehicles return to their bases. When returning, the distance and the time taken are added to the resulting value. Distances and the time they take are set to the same value, as is done in other papers. For example, if a single navigation takes ten units of length, then the time it takes that distance is again taken as ten units of time.
It can be seen from the previous best results that most of them are not directly comparable. Because the parameters they have taken are so different, one needs to fix the options to compare the two results chosen. To make the comparison between the proposed methodology and the others, we tried to find better results using the values they had as parameters to their algorithms.

Distance
The GH data-set contains many individual instances. Each configuration contains end points (or customers/orders) ranging from 200 to 1,000. Each instance may have a different configuration in terms of distance, location and time window. The routes constructed must have Hamiltonian path features. A Hamiltonian path is a graph where each vertex is visited only once and the route must construct a cycle. That is to say, the route must return to the first point or stop. There is only one central depot and many individual end points where orders must be delivered within the time window of each end point. Each customer may have orders in different volumes. The data-set contains the location in discrete X and Y coordinate values. The distance is calculated using Euclidian 2D arithmetic. The duration of a path is equal to the distance value. Each instance has a maximum number of vehicles with a maximum capacity. Only one vehicle is allowed to carry an order. An order cannot be carried with different vehicles in the same schedule (that is to say, no split may occur).
C1_10 instances are clustered instances that allow around 90 routes. They can be called lightly clustered (as opposed to C2_10 instances) instances. C2_10 instances are massively clustered instances and only allow around 30 routes. There are other instance types, such as R1_10, R2_10, RC1_10 and RC2_10, but in this study, we concentrated on C1_10 and C2_10 instances. Table 1 shows some benchmark outputs that we obtained from the data-set. With these results in overall distance, our algorithm is 1.03% better. In overall number of routes our algorithm is 0.33% worse.
ii. Join Operator The Merge Operator cannot construct longer routes in an acceptable computational time. The longer the route length, the longer it will take to generate a single route, because it tries every combination of stops that could be constructed. Usually, when the data-set is good enough, merge will quickly generate routes with four stops, in around one hundred milliseconds. Getting longer routes generally makes the processing time longer and more expensive. To overcome that, a local search method called "join" is created as an operator. This operator searches the remaining solution space for a better candidate with a longer route length.
There are two options that can be set for the Join Operator: the first is the "delta distance". Delta distance is computed between three stops. When the Join Operator tries to insert a new stop into the route, it selects a "navigation" edge. A "navigation" edge can be thought of as a directed graph edge, going from one stop (vertex) to another in the same route. The candidate stops that can be inserted between these two stops (creating two edges), is taken from the delta distance set. A "delta distance set" for an edge is the set of all the stops that, when inserted in between, will not increase the route more than the "delta distance".
In Equation (3), ΔDS is the delta distance set available for the selected edge. D is the two-dimensional distance matrix, v f and v t are from and to vertices (stops), respectively, v s is the newly generated vertex to be created, Δ is the delta distance setting provided to the HBA and S is the set of all stops available.
The other option provided is the "Alpha Distance" setting, which is used for performance reasons. The S set above indicates all the stops available to be joined to the selected edge of the route. Alpha distance filters those stops, yielding a smaller set called "Alpha distance set". Then when calculating delta distance, the stops provided will be taken from this set.
In Equation (4), αDS is the alpha distance set available for v f ; the "from stop". From stop or from vertex is the left-hand side of the edge, when the route is thought out as a directed graph from left to right. D is the two-dimensional distance matrix, s is a stop that is near to v f . is the alpha distance setting provided to the HBA, and again, S is the set of all stops available.
The time it takes to make a single join to a resulting route is considerably less than generating the same route using merge, as the maximum number of operations that will be executed is only proportional to the number of stops available, not to any combinatorial combination of them.
iii. Result elimination settings

Cluster distance
Cluster distance is the same as delta distance in the Join Operator.

Route selection
Decides which routes should be selected for the worse route's order to be distributed upon. Centre of Gravity, Centre of Mass, ALL etc. are the basic options.

Order selection
Decides on the orders of the worse route to be removed and added to the better route. Multiple All, Multiple Random, Single etc. are basic options.

iv. Split Operator
Split Operator is another operator that executes after Merge and Join Operators, together with RE and Swap. This operator also tries to make the schedule better by utilizing cross sections of routes. If two routes visit the same stop, and if the time window and other constraints are preserved, then one route can leave some of its orders on the stop and the other route can pick those orders and take them to their customer locations, but only if the resulting routes creates a better schedule.

v. Swap Operator
Swap Operator is similar to RE Operator. The only difference is that it can exchange the order between routes. That is to say, if an order is selected in one route, then another order can be selected on the other routes. These two orders can be exchanged if the resulting schedule will become better. This order selection is stochastic. Swap Operator may not select two orders but one from a single route and insert it into the other route. If this is the case, it is technically what the Result Elimination does.

vi. Merge Join Distance
Merge Join Distance (MJD) is the maximum distance a stop can increase the route length. If a stop should be added to a route, it cannot increase the route distance by more than the MJD.
For example, Let us assume that • a route has two stops, a and b, • and the route has a length of 10.
• we assume that previously we have set the MJD as 4.
Later we needed to add a new stop, c, to the end after b. The distance between b and c, d [b, c], is 5. Because the MJD is 4, c cannot be added to the route.
Also, if we needed to add a new stop, d, to the end after b and the distance between b and d, d [b, c], is 3. Because the MJD is 4, d can be added to the route.

vii. Finalize Operator
It is always not possible to route all the orders in a set. Usually some orders will remain as unscheduled. It is because if they are routed, they will violate some constraint. Finalize Operator will schedule those orders even if they break the validity of the schedule. This is done partly because usually after running other operators, such as result elimination, these orders can be inserted into a better valid route to remove the violation.