The adaptation of the harmony search algorithm to the ATSP with the evaluation of the influence of the pitch adjustment place on the quality of results

ABSTRACT The paper is an extended version of the conference article, which presents a modification of the Harmony Search algorithm, adapted to the effective resolution of the asymmetric case of the Traveling Salesman Problem. The efficacy of the proposed approach was measured with benchmarking tests and in a comparative study based on the results obtained with the Nearest Neighbor Algorithm, Greedy Local Search and Hill Climbing. The discussion also embraced the study of the convergence of the proposed algorithm and the analysis of the impact of the pitch adjustment place on the quality of the solutions.


Introduction
The Harmony Search (HS) algorithm is a promising metaheuristic used to solve a variety of optimization problems (it has been successfully used in the design of steel frames by Degertekin, 2008, the routing optimization in 4PL with time windows by Bo, Huang, Ip, & Wang, 2009, the Internet routing by Forsati, Haghighat, & Mahdavi, 2008, the optimization of container storage in a harbor area by Ayachi, Kammarti, Ksouri, & Borne, 2010 and in the flexible job shop scheduling problem with multiple objectives by Gao et al., 2016). Its results are usually characterized with the favorable value of the objective function, while at the same time they are achieved in a relatively short time. The nature of the algorithm, which determines the preferential use of the method for continuous optimization problems, requires the application of sophisticated approaches that allow for its adaptation to other ways of representing a number of important issues in business practice.
The Traveling Salesman Problem (TSP) is a classical combinatorial optimization problem that involves finding the shortest Hamiltonian cycle in a complete weighted graph. Its popularity stems from the fact that it belongs to the class of NP-hard problems and that it can model a variety of utilitarian issuesits asymmetric variant (characterized by the possibility of varying weights of edges connecting the same nodes), representing line infrastructure located in urban areas, has become the basis for many logistical problems (it models, for example, the process of planning the mobile collection of waste electrical and electronic equipment, which has been described by Mrówczyńska & Nowakowski, 2015;Nowakowski, Szwarc, & Boryczka, 2018, and the process of transport activities related to the acquisition of municipal waste, described by Płaczek & Szołtysek, 2008;Syberfeldt, Rogstrom, & Geertsen, 2015).
Taking into account the relatively good research results concerning the use of the Harmony Search algorithm to solve many practical problems, a number of controversies over HS (according to Weyland, 2010 it is assumed that the method is non-innovative and it is only a special case of Evolution Strategies), and the utilitarian importance of the asymmetric variant of the Traveling Salesman Problem, we chose to conduct a study aimed at adapting the cited metaheuristics to the combinatorial optimization problem. The proposed topic was also discussed by Antosiewicz, Koloch, and Kamiński (2013), but the results presented by them show the ineffectiveness of the method adapted to TSP, implying the need for an innovative approach to the design of the algorithm facilitating the process of planning the traveling salesman's route.
The paper consists of the following parts: the first part introducing the topic, the second part presenting the Harmony Search algorithm, the third part formulating the asymmetric problem of the traveling salesman, the fourth part proposing the approach to adapt the metaheuristic, the fifth part describing the methodology of the study, the sixth part discussing the results, and the seventh part concentrating on the conclusions and planned work.
The paper is an extended version of an article authored by Boryczka and Szwarc (2018), additionally enriched with the study of the convergence of the proposed algorithm and the analysis of the impact of the pitch adjustment place on the quality of the solutions in our approach to the HS design. In Section 1 we have added information about new applications of the HS and additional references for ATSP applications. We have changed the formulation of the ATSP in Section 3 and added the last paragraph and the Figure 3 in Section 5. Additionally, in the Results we have addded the last four paragraphs and Tables 6-10. Article was also extended by the last paragraph in Section 7.

Harmony search algorithm
The Harmony Search technique, proposed in the paper of Geem, Kim, and Loganathan (2001), is based on the similarity of the jazz improvization process to the search for a global optimum by algorithmic methods. It assumes the existence of a HM structure (referred to as harmony memory), which stores HMS harmonies (Panchal, 2009 said that usually from 4 to 10), consisting of a specified number of pitches (representing the values of the decision variables of a given result). Each HM element is interpreted as a complete solution to a problem whose objective function value is determined based on its components.
The initial harmony memory content is generated randomly and later aligned according to the appropriate objective function values (adapted to the analyzed problem), which describe individual harmonies (so that the first result in HM is the most favorable). These steps initiate the iterative creation of successive solutions.
The procedure for creating a new solution uses the knowledge accumulated in HM and is based on the analogy to the improvization of harmony in music. The development of a solution involves an iterative selection of the next pitch, according to two parameters -HMCR (a harmony memory considering rate; as Ayachi et al., 2010 noticed its value is usually within the range of 0.7 to 0.99) and PAR (a pitch adjustment rate; as Ayachi et al., 2010 said it is often from 0.1 to 0.5). Based on the probability of HMCR, the pitch i is selected, using the values in the i position in harmonies belonging to HM (otherwise the value is generated randomly). Creating a solution based on the HM component, the pitch can be modified with a defined probability of PAR (the change in value is based on the bw parameter, the value of which depends on the representation of a problem).
When the next solution is generated, a comparison of its objective function value with the relevant parameter describing the HM component in the last position is made. When a more favorable result is identified, it replaces the worst result in harmony memory and HM elements are rearranged. The procedure for generating a new solution is performed by IT iterations, and then the best result (in the first position in harmony memory) is returned. The pseudocode of the method is shown in Figure 1.

The formulation of the asymmetric traveling salesman problem
Based on the paper of Öncan, Altınel, and Laporte (2009), the following TSP formulation was adapted: for a directed graph G = (V, A), with weighted arcs represented as c ij (where i, j [ {1, 2, . . . , n}), a route (a directed cycle comprising all n cities) of minimal length is sought. The asymmetric variant of the Traveling Salesman Problem (ATSP) is characterized with the possibility of the occurrence of inequality c ij = c ji .
A decisive variable x ij representing the edge between vertices i and j in the solution found adapts the following values: x ij = 1 when the edge (i, j) is part of the route constructed, 0 otherwise. (1) The objective function was formulated in the following way: The constraints ensuring exactly one visit of the traveling salesman to every city were presented in the following way: x ij = 1, j = 1, . . . , n, n j=1 x ij = 1, i = 1, . . . , n. (3) Such a formulation of the Traveling Salesman Problem could result in the occurrence of solutions representing separate cycles instead of 1 cycle, so it is necessary to introduce additional constraints (MTZ):

The proposal of the approach to HS design
According to the proposed approach to the design of the HS method (tailored to solve the Asymmetric Traveling Salesman Problem), each pitch is represented by integers corresponding to the number of the individual cities visited by the sales agent. The order of their occurrencein harmony representing the complete routeindicates the sequence of a journey. Taking into account the nature of the optimization problem under study, the position occupied by a given city in harmony is deemed irrelevant. It is necessary, however, to consider the sequence of vertices, by selecting the next pitch value based on the generated list of available nodes, occurring in the solutions recorded directly after the last city that belongs to the constructed result. Based on the structure created, the city is selected according to the roulette wheel method (the probability of acceptance of a given item is dependent on the value of the objective function of the solution represented by the length of the route, similarly to the approach discussed in the paper of Komaki, Sheikh, & Teymourian, 2014), or any unvisited node is drawn (when the list of vertices is empty). As a modification of the pitchrelated to the PAR parameterwe opted for the choice (made within the nodes available) of the city nearest to the last visited site in the solution being created (the results of empirical studies assuming the alignment of the representation problem to continuous space and the use of the bw parameter showed the ineffectiveness of this approach in the case of ATSP).
In order to avoid premature convergence (the situation in which the algorithm gets stuck in the local optimum), we introduced the option of resetting the HM elements at the time of execution of a specified number of R iterations from the last replacement of the result in HM. This mechanism assumes that the best result is kept and the remaining solutions are drawn, forcing exploration of the space of solutions. Pseudocode of the proposed approach to the HS design is presented in Figure 2.

Empirical research methodology
The nineteen tasks representing the Asymmetric Traveling Salesman Problem were selected as the 'test bed'. Their characteristics are presented in Table 1. It is assumed that statistically significant resultsfor non-deterministic algorithmscan be achieved by repeating calculations at least ten times for each instance of the problem (Talbi, 2009), so each test was solved 30 times.
The values of the HS parameters were as follows: R=1000, HMS=5, HMCR=0.98, PAR=0.25 (empirical studies were carried out for their designation, the fragmented results of which are shown in Table 2). It was assumed that the algorithm would be executed for 10 minutes, during which the value of the objective function of the best solution found after 2, 6, and  10 minutes would be determined (the same time interval between measurements allows for the observation of changes in the dynamics of the improvement of the solution). Average error was calculated as follows: Algorithms were implemented in C# and the study was conducted on a Lenovo Y50-70 laptop, the parameters of which are presented in Table 3.
For comparative purposesto provide background for the results generated by HSthe solutions obtained by the Nearest Neighbor Algorithm (NNA), Greedy Local Search (GLS) and Hill Climbing (HC) were selected. The first method always started constructing a route from city number one.
The Greedy Local Search variant is characterized with the acceptance of the first designated neighborhood solution, which is described with the more favorable value of the objective function, while Hill Climbing grants acceptance when the entire neighborhood is reviewed. The local search methods used in the study were based on the result found by the Nearest Neighbor Algorithm, and the neighborhood of the current solution was defined as a set of routes differing from it by two cities only, while the initial vertex remains unchanged.
The obtained results are also compared with the solutions presented in literature (Osaba et al., 2014). Due to the different formulation of the algorithm stop condition, they should not be treated as the basis for a direct comparison of the efficiency of metaheuristics, but should only be used to estimate the quality of the proposed solutions.
In order to determine the impact of the pitch adjustment operation on the effectiveness of the proposed HS design approach, an average error was determined for two additional HS variants: HS ′ (in which the operation which depends on the probability PAR can only be executed when the criterion r ≥ HMCR is met) and HS ′′ (in which the pitch adjustment can occur regardless of the HMCR probability). To sum up, HS ′ assumes that the greedy choice can be performed with PAR · (1 − HMCR) probability, while in HS ′′ this operation is  executed with PAR probability (in the classic HS, the probability is equal to PAR · HMCR). The pseudocode of the analyzed techniques is presented in Figure 3.

Results
The average error determined by the methods studied is shown in Table 4. The table also includes the reference results obtained with the Genetic Algorithm (GA) and the Adaptive Multi-Crossover Population Algorithm (AMCPA; results were processed based on Osaba et al., 2014). Based on the results obtained, it was found that the proposed approach to the design of the HS algorithm is characterized with high efficacy. For most instances of the problemthe number of cities below 100the average percentage surplus of the objective function value, in relation to the optimum, was decidedly less than 5% (regardless of the time  needed to apply the method). HS also rated the best results (among the analyzed algorithms; in terms of the objective function value) for 74% of the benchmarking tests. For each of the analyzed time intervals in which the measurement was performed, HS obtained the most favorable results in terms of the average percentage surplus of the objective function values (within the methods under study), ranging from 15.73% to 13.54% (respectively for 2 and 10 minutes).
It is of particular interest that the efficiency of GLS and HC for tasks described by a minimum of 323 vertices is relatively high. Due to the multitude of permissible solutions, the proposed method, together with AMCPA and GA, yielded results with significantly higher values of the objective function. Table 5 presents the values of the objective function determined by the tested solution methods. Accordingly, it was found that NNA, GLS and HC did not provide the optimal solution in any of the analyzed benchmarking tests. The proposed approach to the HS design made it possible to obtain the most favorable result for seven tasks, every time finding a given result in the first two minutes of running the algorithm. In time, one can observe a decrease in the dynamics of the improvement of the solution by HS (the visualization of the dependence is shown in Figure 4), but the process still produces the expected effects while avoiding stagnationit would probably be possible to determine solutions optimal for each test by performing the method longer. In addition, the best results for the ftv38, ry48p, ftv64, and ftv70 tasks are characterized with the objective function value diverging from the optimum only to a small extent, which might demonstrate the need to improve the mechanism of exploiting the method (for example, by hybridizing it with other techniques).
Based on the above compilation, we argue that the algorithm is characterized by a significant standard deviation from the objective function value of the obtained solutions, thus   implying the observable non-determinism of the method leading to significantly different results after the same period of time (in consequence, preventing their prediction). Table 6 shows the compilation of iterations after which the HS algorithm has converged. On its basis, it was found that in order to achieve stagnation (even for relatively small instances of the problem) it is necessary to perform a significant number of iterations. Table 7 presents a compilation of the average error by HS, HS ′ and HS ′′ . Based on the results presented in Table 7, it was found that the lowest total percentage value surplus of the objective function in relation to the optimum is characterized by the variant HS ′′ , whose significant efficiency was recorded for problem instances, which were described by above 43 vertices (which points to the effectiveness of the application heuristics of the Nearest Neighbor, indicating the direction of searching the space for solutions). Regardless of the time of performing the techniques, the reduction in the likelihood of making a pitch adjustment operation (the HS ′ method) negatively affected the quality of the obtained results. Table 8 and 9 show analyzed values of the objective function determined by the HS ′ and the HS ′′ .
To determine the recommended place of pitch adjustment operation, the Wilcoxon Signed-Rank Test was used for the percentage surplus of the objective function value, set by the analyzed HS variants after ten minutes. The value of 0.05 was adapted as a significance level (a lower result indicates that there are no grounds to undermine the alternative hypothesis according to which H1 yielded lower results than H2). The results are presented in Table 10. Based on their analysis, it is recommended that HS ′′ should be used, while the HS ′ variant should be avoided.

Conclusions and future work
The obtained results indicate the relatively good effectiveness of the proposed approach (within the methods under study, it obtained the bestin terms of the total percentage surplus of the objective function valueresults, characterized with the deviation from the optimum ranging from 13.54% to 15.73%, depending on the running time of the algorithm), encouraging further work on its refinement and the adaption of the method to solving other combinatorial problems.
Accounting for a variety of the obtained results close to the optimum and the observations discussed in the paper of Wang, Gao, and Zenger (2015), it is assumed that HS has a weak mechanism of exploitation (its inefficiency was revealed in particular for the tasks where the number of vertexes was above 300, for which the method had far worse results than the results obtained by simple local search methods that were based on the NNA), implying the need to conduct research on the hybridization of this method with other techniques to eliminate the imperfection. The high effectivness of the HS ′′ variant achieved in the process of solving the instances of the problem that are described by a significant number of vertices and its lower effectiveness (compared to HS) for relatively small problems points to the potential appplicability of a dynamic value of the PAR parameter, which would combine the strengths of the two approaches. Regardless of the scale of the problem under analysis, we have noticed that the reduction in the likelihood of making a pitch adjustment operation (the HS ′ method) negatively affected the quality of the obtained results.

Disclosure statement
No potential conflict of interest was reported by the authors.

Notes on contributors
Urszula Boryczka currently works at the Institute of Computer Science, University of Silesia in Katowice. She completed her Master's degree at University of Silesia in Sosnowiec in 1984 and she received Ph.D. degree at University in Wroclaw in 1993. She holds Doctor's of Science (Habilitation) degree in Computer Science from the Polish Academy of Science in Warsaw in 2009. Urszula Boryczka is a head of Division of Algorithmic and Computational Intelligence. She does research in algorithms and artificial intelligence, especially in computational swarm intelligence. Swarm intelligence is the discipline that deals with natural and artificial systems composed of many individuals that coordinate using decentralized control and self-organization. In particular, the discipline focuses on the collective behaviors that result from the local interactions of the individuals with each other and with their environment. Examples of systems studied by swarm intelligence are colonies of ants and bee colonies, flocks of birds and many others. ACO, BCO, PSO, CS and BA are examples of such systems examined in global or combinatorial optimization systems. Her current project is 'Ant colony algorithms in clustering problems'. Another issue is a 'Harmony search' which is a new metaheuristics in her study. This technique is now examined as a hybrid approach with other optimization techniques in aTSP problems.
Krzysztof Szwarc is currently a research assistant at the University of Silesia in Katowice. He received a B.Eng. degree in Logistics (2015) and an M.Eng. degree in Transport (2016)