The Pheromone-Based Harmony Search Algorithm for the Asymmetric Traveling Salesman Problem

: This paper presents a modiﬁcation of the Harmony Search algorithm (HS) adjusted to an effective solving of instances of the Asymmetric Traveling Salesman Problem. The improvement of the technique spans the application of a pheromone, which, by serving the role of long-term memory, enables the improvement of the quality of determined results, especially for tasks characterized by a signiﬁcant number of vertices. The publication includes the results of tests that suggest the achievement of effectiveness improvement through the modiﬁcation of the HS and recommendations concerning the proper conﬁguration of the algorithm.


Introduction
Due to the considerable complexity of many problems occurring in business practice, the use of classical deterministic algorithms that enable determination of exact solutions becomes impossible in an acceptable time. For this reason, approximate methods are increasingly used, among which one can distinguish metaheuristics. These algorithms have found application in the optimization of logistic processes, e.g., Nondominated Sorting Genetic Algorithm II and Multiobjective Particle Swarm Optimization were used for sustainable ship routing and scheduling with draft restrictions [1], Non-Dominated Sorting Genetic Algorithm II was applied for the design of the Shuttle-Based Storage and Retrieval System [2] and Particle Swarm Optimization-Composite Particle was used for sustainable integrated dynamic ship routing and scheduling optimization [3], as well as sustainable maritime inventory routing problem with time window constraints [4].
The Harmony Search method (HS) is an interesting metaheuristic, which despite causing a great many controversies (it is falsely [5] accused of a lack of innovativeness and similarity to the Evolution Strategies (ES) [6]), found application in the process of solving many utilitarian problems (it was applied, among others, in the process of flood control management [7], in the fourth-party logistics routing problems [8] and in the optimization of allocation of containers around a harbor [9]). The character of the technique, implying its preferential usage to the continuous space, imposes the application of sophisticated approaches aimed at adjusting the algorithm to combinatorial optimization problems, at the same time influencing the appearance of difficulties connected with an effective designing of a method of the HS adjusted to solving a significant number of practical problems.
One of the most tested concepts of combinatorial optimization is the Traveling Salesman Problem (TSP). It is characterized by unquestionable utilitarian meaning-its asymmetric version (ATSP) models the linear infrastructure located in urbanized areas, thus becoming groundwork for the planning of routes of vehicles collecting municipal waste [10] and waste electrical and electronic equipment [11]. The use of ATSP goes beyond the area of reverse logistics (that can be defined as 'the

A Classical Harmony Search Algorithm
According to [19] we can assume that the aesthetic quality of a musical instrument is determined by its pitch (that in Musical Instrument Digital Interface is often represented as a numerical scale), timbre and amplitude. Timbre is determined by the harmonic content that depends on the pitch or frequency range of the particular instrument. When we adjust the pitch, we are trying to change the frequency.
Zong Woo Geem in [20] suggested a controversial metaheuristic that was based on the similarity of the process of jazz improvisation to looking for global optimum through the means of algorithmic methods. The technique assumes the existence of a structure HM (referred to as the Harmony Memory), which stores HMS harmonies that consist of a set number of pitches (representing the values of decision variables of a given result). The particular HM elements are interpreted as a complete solution of the problem, the value of the objective function of which is defined on the basis of its component parts.
The initial content of the Harmony Memory is generated randomly, after which it is subjected to sorting on the basis of appropriate objective function values (in such a way as to make it possible that the result in the first position can be characterized by the best result). The realization of the described steps brings about the launch of an iterative creation of subsequent solutions.
The procedure of creating a new solution applies the knowledge stored in HM and is based on an analogy to the process of harmony improvisation in music. The construction of such a solution consists in an iterative selection of the next pitch, according to two parameters-HMCR (Harmony Memory Consideration Rate) and PAR (Pitch Adjustment Rate). On the basis of probability HMCR the pitch i is selected, using the values located in position i in harmonies that belong to HM (otherwise random generation of an acceptable value takes place). Creating a solution based on a component HM a modification of the pitch with the set probability PAR may occur (the change of value takes place on the basis of parameter bw, the value of which depends on the representation of the problem). To sum up the random selection of the pitch value is done with probability 1 − HMCR, using an unmodified value from HM may occur with probability HMCR · PAR, whereas modified value from HM is chosen with probability HMCR · (1 − PAR).
After generating the subsequent solution, comparison of its objective function value with the appropriate parameter describing the component located in the last position in HM is held. In the situation of determining a better result, it replaces the worst result stored in the Harmony Memory and another sorting of the HM elements is conducted.
The procedure of generating a new solution is conducted by IT iterations, after which the return of the best result occurs (located in the first position in the Harmony Memory). Pseudocode of the HS was presented in Algorithm 1.

Algorithm 1
The Harmony Search pseudocode (based on [15] for i = 0; i < n; i + + do 8: Choose random r ∈ (0, 1) 9: if r < HMCR then 10: H[i]=choose randomly available pitch on position i in HM 11: Choose random k ∈ (0, 1) 12: if k < PAR then 13: α = bw· random ∈ (−1, 1)  24: iterations + + 25: end while 26: return HM[0] The existence of a visible similarity between the HS and the ES caused a great deal of controversy (e.g., [6]), however in the reference article by [5] a crucial difference in the means of functioning of exploration and exploitation, present in the quoted techniques, is emphasized. The HS has three operators, the occurrence frequency of which is controlled by the values of particular parameters describing the algorithm, whereas the strategy (µ + 1) bases its functioning only upon two obligatory operations -crossover and mutation.

Formulating Asymmetric Traveling Salesman Problem
Based on article [21] the following definition of the TSP was adapted: for a directed graph G = (V, Ar), with arc weights c ij (i, j ∈ {1, 2, . . . , n}), a route is being searched (directed cycle including all n cities) characterized by a minimum length. The asymmetric version assumes the possibility of occurrence of inequality c ij = c ji .
The decision variable x ij , represented by the presence of an edge between the vertex i and j in the selected solution, adapts the following values: The objective function, which assumes the minimization of the values of edges included in the solution, was formulated as: Limiting conditions-ensuring that the salesman visits every city on his way exactly one time-were presented in the following way: In order to avoid the occurrence of solutions representing separate cycles, instead of a connected one, it is crucial to introduce additional limitations called Miller, Tucker and Zemlin (MTZ): MTZ assumes the use of u i variables to determine the order in which each vertex i is visited by salesman during the trip. According to it, the numbers from 1 to n are being assigned to the vertices in such a way that the numbering corresponds to the order of the vertices in the route.

Harmony Search Hybridization
The adaptation of a given metaheuristic approach to solve a given optimization problem is a kind of a metaphor, which brings about a suitable reaction of the system to the conditions of the changing environment, which, in turn, is the surrounding of the analyzed problem. In order to prepare a suitable ability of fast reaction in approaches solving difficult optimization problems, a great many scientists aim at creating new algorithms, which cope with the difficulties of limitations, or the multitude of the analyzed dependencies between data. All the factors contribute to a snowballing increase in the number of various metaheuristic approaches. However, it is not always that new solutions skillfully cope with finding the best solution. Usually, there is a tendency to create a general approach that would correspond with a broad range of applications. However, knowing the No Free Lunch theorem (according to it, there is no universal, best optimization algorithm for all tasks, the theorem is described in the paper [22]), we should firmly state that it is impossible. As a result, approaches spanning several metaheuristics into one, referred to as hybridization, are used. Such a connection can have a homogeneous or non-homogeneous character and as a result, it can be realized in a sequential or parallel way. Most frequently, however, a hybrid of optimization approaches is created as a consequence of a lack of some mechanisms, essential to solve a given problem and simultaneously absent in the homogeneous approach.
Looking at hybridization in the HS one might be sensitive to the following issues: lack of a visible mechanism of machine learning (or lack of experience) during solving a given optimization problem, lack of a clear criterion of accepting worse solutions, or finally lack of official local search/greedy algorithm, which would enable fast convergence of the hybrid algorithm. Owing to that, in this part of our work we will concentrate on the possibility of connecting the HS with other techniques copying nature. The following will be briefly discussed: Ant Colony Optimization (ACO), Neural Networks (NN), Particle Swarm Optimization (PSO), Genetic Algorithm (GA) and others.
One of the approaches to the HS and the ACO is the solution of the problem connected with an optimal installation location of structural dampers [23]. Two ideas stemming from the ACO were added to the HS: pheromone trail and some value of the heuristic function. The information was automatically ascribed to every component of the solution. Owing to such an approach to the process of obtaining a better solution over a short period of time was accelerated.
A hybrid connection of the HS with the NN took place in the problem of selecting a proper proportion while producing industrial concrete mix [24]. The optimization of the suggested mix is obtained owing to the refining of this mix through steering the proportions using Neural Networks.
The connection of the HS with an immune system can be analyzed on the basis of the problem of scheduling the work of nurses, which belongs to the NP-complete class. Owing to such an approach a proper balance between local and global search was enabled and simultaneously premature convergence was prevented.
In paper [25] a new approach was described -a connection of the HS with the Cuckoo Search (CS) in the 0-1 knapsack problem. Similarly, as in the previously analyzed article, the Global-best HS (GHS) was responsible for exploration and the CS for the exploitation of the solution space. The hybrid CSGHS was definitely better than the binary version of the CS.
In the Optimal Power Flow Problem, the BCGAs-HSA methodology was suggested in [26], which connected the binary coding of the GA with the HS algorithm. In both cases for the GA and the HS, the hybrid yielded better results in every instance of the analyzed problem.
The most complicated hybrid was applied to the truss structures problem, where the HS was connected with the PSO and the ACO [27]. In this approach, the PSO with the passive congregation was applied. The PSOPC was used for global search and the ACO was applied to updating the position of particles searching the solution space. The Fly-back mechanism was also used here. The above hybrid was compared with various variants of the PSO.
It can be noticed that more and more often [28] the HS is treated, with some modifications taken from the Swarm Intelligence, as a technique of global search, hence there are numerous modifications connected with the local search. The missing learning mechanism is located/observed by adding an approach connected with swarm algorithms (bees)-HHSABC [29]. In the global optimization, the hybrid turned out to be one of the best hybrids solving this problem. What is worth mentioning is the dynamic form of adapting the values of parameters in the HS, which contributes to the increase of adaptability in generality.

Pheromone-Based Harmony Search for ATSP
This section comprises three parts. The first one presents modifications that enable adjustment of the HS to solving instances of the ATSP, the second presents the application of a pheromone in Ant Systems adjusted to searching the route of the salesman, and finally the third one-the structure of the PBHS.

Harmony Search Adjusted to ATSP
This work was based on a modification of the algorithm suggested in [15]. It assumes representing each pitch by integers, corresponding with the numbers of particular cities, which the salesman is to visit. The order of their occurrence-in a harmony representing the complete route-points towards the sequence of the travel.
While creating a new harmony the sequence of vertices is being considered, which is realized through choosing the next value of pitch based on the generated list of available nodes, occurring in the memorized solution immediately after the last city, which belongs to the result under construction. On the basis of the created structure a city is chosen according to the roulette wheel selection (probability of accepting a given element is dependent on the value of solution objective function, represented by the length of the route, analogically to the approach presented in [30]) or a random non-visited node is drawn (in case when the created list of vertices is empty). As a modification of pitch-connected with the parameter PAR-a selection of city (conducted within the area of available nodes), located in the closest proximity of the most recently visited city in the created solution, was adapted.
What was introduced in order to avoid premature convergence is a possibility to reset elements HM, when the particular number of iterations R is executed, from the last replacing of the result in the harmony memory. The mechanism assumes storing the best results and drawing the remaining solutions.

The Role of Pheromone in as Adjusted to Solving TSP
One of the key elements describing the AS is the existence of a pheromone, owing to which the particular agents (the number of which is m) communicate with one another. The route chosen more often by the ant is characterized by a higher intensification of a pheromone, encouraging, as a result, other insects to follow the popular route. The probability of an agent going from the city i to a not visited city j, in time t, is defined on the basis of the following formula: where: τ ij (t) is the size of pheromone trail on the edge connecting the node i and j in time t, η ij is the inverse of the weight of edge between the vertex i and j (short length routes are promoted), n is the number of cities, and parameter β is used to steer the significance of heuristic and pheromone. Among the various approaches to updating the pheromone, what deserves particular attention is an effective [17] Ant-cycle System, in which the agent k leaves the pheromone after constructing the entire route. ∆τ k ij (t, t + n) adapts the following values: where: Q 3 is a constant (usually taking on value 1 [31]), and L k -length of the route constructed by an ant k. After executing n steps the updating of the pheromone trail occurs, according to the following formula: where: ∆τ ij (t, t + n) = ∑ m k=1 ∆τ k ij (t, t + n) and ρ is parameter representing evaporation of pheromone (value (1 − ρ) ∈ [0, 1] [31]; used to reduce the impact of historical results and promote the latest one).
The initial size of the pheromone trail τ 0 can be established regardless of the characteristic of a given instance of the problem (e.g., assuming that τ 0 = 1, according to paper [32]) or more sophisticated approaches can be applied, such as e.g., the way suggested in paper [33] based on the Nearest Neighbor Algorithm (NNA). The second of the quoted methods assumes the usage of the following formula: where L nn is the length of the route constructed by the NNA.

Pheromone-Based Harmony Search for ATSP
In order to increase the effectiveness of the HS, suggested in [15], the process of constructing a new harmony was complemented with the knowledge accumulated by the means of a pheromone (performing the role of long-term memory). To preserve the idea describing the HS, activation of the described mechanism on the basis of HMCR probability was assumed -when the list including successors of the recently visited cities is empty, a selection of the next city, on the basis of pheromone trail, takes place. The updating of the HS global memory is conducted according to Formula (7), after every iteration has been run. Based on the reference books it was assumed that Q 3 = 1 [31] and τ 0 = 1 [32]. The probability of choosing a pitch in position s in harmony H is defined on the basis of a modified Formula (5) (based on the Simple-ACO approach, suggested in [32]), i.e., where: value i corresponds with the number of a city that is in harmony H in position s − 1, and A is a set of cities available for visiting. The introduction of pheromone memory is connected with the creation of a structure, the size of which is n 2 , and which can influence the problematic implementation of the suggested solution designed for large ATSP instances. In order to eliminate the indicated imperfection it is possible to apply selective pheromone memory, according to the approach suggested in [34]. Pseudocode of the PBHS for ATSP was presented in Algorithm 2.

Research Methodology
The effectiveness of the suggested approach, taking into consideration the application of a pheromone, was tested on the basis of nineteen instances of the ATSP the characteristic of which was presented in Table 1 (on the basis of [35], where the optimal (best obtainable) route lengths for particular tasks are presented; they are available in the TSPLIB library, described in [36]). In addition, the effectiveness of the methods was verified on the basis of a case study, involving the visit of 40 points located in the Silesian Voivodeship in Poland (it can model e.g., problems related to the distribution of many products; the distance matrix has been prepared with an accuracy to the one meter). A case study map with the locations marked is shown in Figure 1. All tasks were solved 30 times by the HS algorithm and three variants of the PBHS (differing in the value of the parameter ρ). To solve the case study we also used the NNA that is often used when companies do not use the right optimization software.
On the basis of work [15] the following parameter values were set: R = 1000, HMS = 5, HMCR = 0.98 and PAR = 0.25. Additionally, when analyzing the convergence of the technique (presented in [37]) it was assumed that IT = 1,000,000, and as the value of the parameter ρ: 0.25, 0.5 and 0.75 were selected (tests were conducted for each of them in order to select the recommended approach to the designed method).
The algorithms analysed in this work were implemented in language C#, and tests were conducted on a Lenovo Legion Y520 laptop (with: Intel Core i7-7700HQ, 32   Both the average error (defined on the basis of Formula (10)) and the actual time of running algorithms were subjected to evaluation. The second measure makes it possible to take into consideration time overhead connected with the introduced modifications and as a result enables defining of the usefulness of the technique in practical adaptations. In addition, we analyzed the number of iterations, after which the algorithm achieved convergence for tasks from TSPLIB.
The effectiveness of the proposed method was assessed by comparing the obtained average error with the results achieved by the following algorithms appearing in the literature on the subject: We chose them, because the first two methods are population-based (like the HS), while the rest of them are classic algorithms commonly used to solve the ATSP instances.

Results
This section comprises two parts. The first one presents results for tasks from TSPLIB library, while the second presents the results for the analyzed case study.

Tasks from TSPLIB
The breakdown of the average error achieved by the particular methods was presented in Table 2 (the best results between the listed methods are in bold). Based on it, it was established that the algorithm devoid of a pheromone defined more favorable results than the PBHS only for two benchmark tests described by the means of a relatively small number of vertices (36 and 45). The effectiveness of the suggested modification was especially emphasized for instances described by minimum 171 nodes (for task rbg358 limitation of average error by almost 25% was noticed). Regardless of the applied value of parameter ρ, usage of the pheromone substantially increased the effectiveness for the tested method (average error was limited from 13.42 to 10.15, 10.92 and 11.68% for ρ = 0.25, ρ = 0.5 and ρ = 0.75 respectively).
Taking into consideration the total average error and number of tasks, for which the lowest result was observed -among the results defined by the analyzed methods -the most favorable PBHS variant is the technique the parameter value of which is ρ = 0.25. The smallest effectiveness is characteristic of the PBHS equipped in a fast pheromone evaporation mechanism (ρ = 0.75), which constructed the best route only for the p43 test (difference of the average error between the created solution and results achieved by the remaining methods equals merely 0.01%). The breakdown of the average error achieved by PBHS and selected methods from the literature on the subject was presented in Table 3 (the best results between the listed methods are in bold) and Figure 2. On their basis, the effectiveness of the proposed approach was found to be significant in relation to the selected techniques described in the literature on the subject. Only for the largest tasks the method obtained a less favorable average error, indicating the need to increase the exploitation process that dominates in GLS and HC.  The breakdown of the running time of particular methods was presented in Table 4. According to the achieved results, the total time of running the method increased by 50% in relation to the technique devoid of the long-term memory. What deserves particular attention is the increase in the time of realizing IT iterations by the PBHS for ρ = 0.25, in relation to the variants with different values of the indicated parameter. The shortening of the time of running the tasks was due to optimization conducted by the compiler, as a result of which the multiplication operation (of the value (1 − ρ)) was substituted by shifting (for ρ = 0.5 and ρ = 0.75 was 2 −1 and 2 −2 achieved). The insignificant differences which appeared when the PBHS was running with ρ = 0.5 and ρ = 0.75 were probably caused by an unequal processor overload.
The number of iterations after which different PBHS variants have reached convergence is shown in Figure 3 (using the popular box plot, described e.g., in the paper [38]). On their basis, there were no significant differences between the individual PBHS settings to achieve convergence. In addition, the scatter of iterations indicates the choice of the not too low value of IT for the tested set (only for the smallest task br17 it could be reduced without any loss of efficiency).
The average error-defined by the particular methods-was subjected to the Wilcoxon Signed-Rank Test, with the application of the R package and the wilcox.test function (applying the following values of parameters: paired = TRUE, alternative = "less", exact = F, correct = F). The following pairs of methods were used in the process (marked as O1 and O2), and the assumed test significance level value was 0.05 (achieved p-values, which were characterized by a lower result, were distinguished by bold and they cause to accept the alternative hypothesis, according to which O1 achieved lower results than O2). The results of the test were presented in Table 5. On their basis, it was established that taking into consideration only the value of the objective function of the solution, it is recommended to use the PBHS with ρ = 0.25. Additionally, no grounds were found to negate the hypothesis which claims that the usage of a pheromone makes it possible to achieve better results than the ones defined by the HS without global memory.   Figure 4 presents a comparison of results (total distance) obtained for a case study by PBHS, HS and NNA (starting to create a solution from node number 1). Based on them, the effectiveness of the developed PBHS was found -it allowed to shorten the average length of the route from 278.12 (HS results) to 276.67 (ρ = 0.25), 276.05 (ρ = 0.5) and 275.88 (ρ = 0.75) km. It is also worth emphasizing the higher predictability of PBHS results (smaller spread) than HS results.

Case Study
Limiting the average length of the route by 21.58 km, in relation to the result constructed by NNA, indicates significant benefits that may arise as a result of the use of PBHS in business practice. The relatively small reduction in the length of the route by PBHS compared to the HS result may allow for a significant reduction of costs when frequent route optimization is required.
In contrast to previous results, the recommended ρ value for the case study is 0.75, however, the slight difference between the length of the obtained routes indicates the possibility of the influence of nondeterminism on the results.

Conclusions and Future Work
As a result of the conducted research the effectiveness of the application of pheromone in HS adapted to solving the ATSP was confirmed. The combination of the advantages of the algorithm proposed in [15] with the version of the Simple-ACO (suggested in [32]), allowed to reduce the average error from 13.42 to 10.15%.
A comparison of the results obtained with the results of selected methods from the subject literature indicates significant effectiveness of the proposed approach, allowing it to be used in business practice. In addition, the presented case study showed significant advantages of the method and verified its suitability both in terms of execution time and the possibility of optimizing the route.
For objectives requiring the relatively fast achievement of a result setting the value of the ρ parameter in PBHS at 0.5 is recommended. In the remaining cases, the application of ρ = 0.25 is recommended. Further work concerning PBHS can concentrate on a more detailed adjustment of the value of parameter ρ and conducting the hybridization of the technique with different methods, aimed at increasing the effectiveness of the suggested approach. Additionally, what is recommended is to prepare a parallel version of the PBHS (e.g., based on the approach presented in [39]) and adjusting the algorithm to the remaining combinatorial optimization problems. It is also worth checking the efficiency of the algorithm for another τ 0 value (e.g., using the approach proposed in [33]). Comprehensive studies should also be performed to compare the effectiveness of the proposed PBHS and other HS variants, including Improved HS [40], GHS [28] and Self-adaptive GHS [41].