A novel design of differential evolution for solving discrete traveling salesman problems

Journal Pre-proof A novel design of differential evolution for solving discrete traveling salesman problems Ismail M. Ali, Daryl Essam, Kathryn Kasmarik PII:

S2210-6502(19)30446-8

DOI:

https://doi.org/10.1016/j.swevo.2019.100607

Reference:

SWEVO 100607

To appear in:

Swarm and Evolutionary Computation BASE DATA

Received Date: 10 June 2019 Revised Date:

19 September 2019

Accepted Date: 25 October 2019

Please cite this article as: I.M. Ali, D. Essam, K. Kasmarik, A novel design of differential evolution for solving discrete traveling salesman problems, Swarm and Evolutionary Computation BASE DATA (2019), doi: https://doi.org/10.1016/j.swevo.2019.100607. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

A Novel Design of Differential Evolution for Solving Discrete Traveling Salesman Problems Ismail M. Alia,∗, Daryl Essama , Kathryn Kasmarika a School

of Engineering and Information Technology (SEIT), University of New South Wales, Canberra, Australia

Abstract Differential evolution is one of the most powerful and popular evolutionary algorithms, which is primarily used to solve continuous-based optimization problems. Although differential evolution has been considered unsuitable for solving many permutation-based real-world combinatorial problems, several efforts for designing an efficient discrete version of the differential evolution algorithm have occurred in recent years. This paper introduces a novel discrete differential evolution algorithm for improving the performance of the standard differential evolution algorithm when it is used to solve complex discrete traveling salesman problems. In this approach, we propose to use a combination of, 1) an improved mapping mechanism to map continuous variables to discrete ones and vice versa, 2) a k-means clustering-based repairing method for enhancing the solutions in the initial population, 3) an ensemble of mutation strategies for increasing the exploration capability of the algorithm. Finally, for improving the local capability of the proposed algorithm when solving discrete problems, two well-known local searches have also been adapted. To judge its performance, the proposed algorithm is compared with those of 27 state-of-the-art algorithms, for solving 45 instances of traveling salesman problems, with different numbers of cities. The experimental results demonstrated that our technique significantly outperforms most of the comparative methods, in terms of the average errors from the best-known solutions, and achieved very competitive results with better computational time than others. Keywords: Discrete differential evolution, Combinatorial optimization problems, Traveling salesman problems, Local search, Mapping method

1. Introduction Differential evolution (DE) is a well-known stochastic population-based search technique, which belongs to the family of evolutionary algorithms. It was ini∗ Corresponding

author Email addresses: [email protected] (Ismail M. Ali), [email protected] (Daryl Essam), [email protected] (Kathryn Kasmarik)

Preprint submitted to Elsevier

October 29, 2019

tially developed by Storn and Price [55]. DE was inspired by the Darwinian theory of evolution, which mimics the biological model of evolution and natural selection. DE is a fast and simple technique that performs well on a wide variety of continuous optimization problems. It has been successfully applied in various domains and fields, such as communication [60, 11], power systems [44], signal processing [13, 9], pattern recognition [31], and mechanical engineering [41]. In its process, DE iteratively evolves a population of candidate solutions throughout several generations using evolutionary operations, such as mutation, crossover and the biological concept of natural selection, to obtain better individuals. This allows individuals with higher quality to survive from generation to generation, where the quality of each individual can be calculated by a predefined fitness/objective function. DE has been successfully applied for solving a wide range of optimization problems, which are characterized by continuous variables/domains. This is because of the learning concept of its mutation operator, which depends on the Euclidean distance measurement between more than one individual. This learning concept has the capability to efficiently explore a search space and allows DE to find new candidate solutions. Moreover, many modifications have been applied to its crossover and mutation operators to further improve the performance of the DE while solving problems in continuous domains [69, 49]. On the other hand, DE has not generally been considered as applicable for solving optimization problems characterized by permutation-based combinatorial variables. This DE deficiency has made it inapplicable, for some time, for solving several real-world combinatorial optimization problems (COPs). Due to their practical utility in various disciplines [46], such as scheduling [4, 30, 45], control systems engineering [25, 22], manufacturing [43], and power systems and electrical engineering [8, 51], COPs have been widely studied using several algorithms. In the literature, many evolutionary algorithms have been introduced for solving COPs. In this paper, these algorithms are classified into two types: discrete-based and continuous-based algorithms. In the algorithms of the first type, the evolving operators are primarily suited for discrete problems. The solutions, in this type, are perturbed without change to their discrete format, by simply swapping the positions of two genes or changing the position of a gene by inserting it in a new position, and so on. While in the second type, the operators generate the new solutions by predefined mathematical formulas, which are commonly working on real numbers. This makes these operators more suitable for solving optimization problems with continuous variables. However, modifications can be applied to algorithms in type one (and two) to be able to solve continuous (and discrete) problems, respectively. Genetic algorithms (GAs) and ant colony optimization (ACO) algorithms are examples of type one algorithms, whereas DE, particle swarm optimization (PSO) and bat algorithms are examples of type two. In this study and among the COPs, we will focus on the traveling salesman problem (TSP). TSP is a well-known combinatorial optimization problem with discrete decision variables that has been proven to be an NP-complete problem [48]. In recent years, the PSO algorithm has been widely used for solving dif2

ferent COPs [61, 63, 45] and particularly TSPs [38, 14, 33]. The most recent work was proposed in 2019 [29], where PSO was combined with GA for solving standard instances of TSPs. In 2018, a PSO with a new learning concept that applies the acceptance criterion of the simulated annealing algorithm, was also proposed for solving TSP [66]. Moreover, many other metaheuristics have been proposed with the same aim, such as the discrete bat algorithm [47], ACO algorithms [15, 28], the artificial bee colony algorithm [16], GAs [12, 1, 37], simulated annealing [26, 39], neural networks [23] and many other heuristics algorithms [52]. However, improvements for increasing their performances in many aspects, such as convergence speed, computational time and their ability to solve large-sized discrete problems are still needed. Additionally, a limited number of DE works for TSPs can be found in the literature. In 2010, an improved DE with a new operator for regulating integer sequences to a mutation process, and that applied the Liuhai crossover instead of the original crossover operator, was introduced for solving TSPs. In 2011, a hybrid differential evolution (HDE) algorithm for solving TSPs of middle-size was proposed. The authors applied position-ordering mapping, by involving a hill-climbing algorithm. However, those two papers have only been tested on a small number of benchmark instances, where the largest problem was “KroA150” with 150 cities. Recently, in 2016, a DE algorithm, called real adjacency matrix-coded differential evolution algorithm (RAMDE), was proposed for solving TSPs. In this work, the individuals are represented by a swarm of real adjacency matrices, and the arithmetical/evolving operators of DE are executed in the form of real matrices. Although the algorithm could achieve solutions better than some other algorithms in the literature with lower fitness evaluation numbers, the algorithm could not reach the optimal solution for any solved instances. On the other hand, recent discrete versions of DE have been successfully applied for solving motor train-sets [64], multi-skill resource-constrained project [42] and flow shop [68] scheduling problems with very promising results. Also, a binary version of DE has been proposed for solving different instances of binary knapsack problems with outstanding performances compared with other algorithms [5]. The limitations of applying discrete versions of DE for solving TSPs and the promising results achieved by discrete DE in other problems, have formed the main motivation of this work. In this paper, a novel discrete version of DE (NDDE) for solving permutationbased TSPs is introduced. The proposed DE incorporates several efficient components, that each has its own contribution to improving DE’s performance. Firstly, an enhanced mapping method is applied to the randomly generated initial population with continuous variables, for mapping those continuous variables to discrete. Secondly, an improved k-means clustering is adapted as a repairing method for enhancing the discrete solutions in the initial phase. Then, a copy of the modified discrete population is generated in continuous format by the proposed backward mapping method, in order to enable DE to continue its evolutionary processes (such as mutation and crossover) in the continuous landscape. After that, an ensemble of DE continuous mutation operators is proposed for increasing the exploration capability of DE. On the other hand 3

and for improving the exploitation of DE, the well-known 3-opt [35] and an enhanced double bridge [54] local search are applied. Finally, the performance of our proposed NDDE is compared with those of 27 state-of-the-art algorithms for solving instances of TSPs with different dimensions. 11 of the comparative algorithms are based on GA, as it is mainly designed for problems in a discrete domain. The experimental results show that NDDE performs well on all the test instances and that it can hence obtain high-quality solutions. The rest of this paper is organized as follows. Section 2 describes the TSP definition and its mathematical formulation, and shows the original DE structure in the continuous domain. Section 3 represents the proposed framework of the discrete DE algorithm. Section 4 discusses the experimental analyses and results, and conducts several comparisons with other algorithms. Finally, conclusions, which are drawn from this study, and future work, are provided in Section 5. 2. Preliminaries 2.1. Traveling salesman problem The traveling salesman problem (TSP) is a well-known combinatorial optimization problem with discrete decision variables. TSP has been proven to be an NP-complete problem [48]. In this paper, the symmetric TSP is used. It considers a salesman who must travel between N cities with the primary target of visiting each city during the trip, and returning back to the original city, with minimum distance. Cities must be visited only once and can be visited in any order (i.e. the distance from node 1 to node 2 equals distance from node 2 to node 1). Each city is connected to other cities/nodes by connections, such as roads, airplanes, or railway. Each of those connections has its own cost, which may be one or more. Those costs may take several forms. For instance, the cost of an airplane or train ticket, or the length of the road, or the time required to finish the trip. In this paper, the base version of TSP, which represents the cost as a distance between cities, is used. So that the salesman needs to keep the distance he/she travels as low as possible. Since TSP is an NP-hard combinatorial optimization problem [32], and thus the running time of an algorithm used for solving a TSP increases exponentially with the number of cities, it has attracted several mathematicians and computer scientists. Moreover, it has several real-world applications in science and engineering. For example, vehicle routing, scheduling, transportation and logistics, the design of hardware devices and radio electronic systems, and computer networks [34]. Mathematically, the objective function of TSP can be described as the minimization of the total distances between the given visited cities, in addition to the distance to return to the start point/city, as follows: M inimize : D, where

4

D = dN,1 +

N −1 X

di,i+1

(1)

i=1

where D is the total distance of the trip, N is the number of cities to be visited, i = 1, 2, ..., N , dN,1 is the return distance from the last city (N ) to the first one, and di,i+1 is the distance between two consecutive cities i and i + 1. 2.2. Differential evolution (DE) in continuous space DE optimizes a problem by iteratively trying to improve its candidate solutions by using DE’s operators. In DE, the quality of each solution can be measured by a given fitness function. It uses three evolving operators, mutation, crossover and selection, to guide the search to find (near-) optimal solutions. DE is considered a powerful tool for solving optimization problems in continuous spaces. The structure of DE, shown in Algorithm 1, is as follows: firstly, an initial population with a predetermined size (P S) is generated and − then each individual (→ xi ), which consists of N variables, is evolved in the iteration phase using three evolutionary operators. These operators are discussed in the following subsections: Algorithm 1 Structure of original DE 1: Initial random population (P S) 2: while termination condition not met do 3: for each individual i (i = 1, 2, . . . , P S) do 4: Mutation; 5: Crossover; 6: Selection; 7: end for 8: end while

2.2.1. Mutation operation − − For each target vector (→ x i ), a mutant vector (→ v i ) is generated, which in its simplest form [55], is calculated as: → − − − − v i,g+1 = → x r1 ,g + F × (→ x r2 ,g − → x r3 ,g ), r1 6= r2 6= r3

(2)

where F is the mutation parameter, or weighting factor, that controls the am− − plification of the differential variation (→ x r2 ,g − → x r3 ,g ) and generally lies within → − → − → − the range of [0, 2], and x r1 ,g , x r2 ,g , x r3 ,g are three randomly chosen vectors − which are not equal to either each other or the target vector (→ x i,g ). Moreover, the mutant vector of any generation can be produced by incorporating the best target vector in that generation. This DE mutation scheme can be applied through the DE/best/1 strategy [19] as: → − − − − v i,g+1 = → x best,g + F × (→ x r1 ,g − → x r2 ,g ) 5

(3)

− where → x best,g is the best individual vector in generation g. 2.2.2. Crossover operation − In a crossover operation, trial vectors (→ u ) are generated by combining the target and mutant vectors (offspring) according to a predefined possibility. The two most well-known types of crossover in the literature are binomial and exponential. In a binomial crossover operation, the trial vectors are generated according to equation 4:  −  → v ji,g+1 , rand(j) ≤ CR or j = aj  → − u ji,g+1 = (4)   j − → x i,g , otherwise where j=1,2,. . . ,n and i = 1, 2, . . . , P S, with CR the crossover possibility in the range of [0,1], rand(j) is the j th evaluation of a uniform random number generated within [0,1] and aj is a randomly selected dimension, to ensure that − at least one element of → u i,g+1 is chosen from the mutant vectors [55, 50]. The exponential crossover first selects a starting cut point, which is randomly chosen from {1, 2, ..., N }, and e consecutive components (counted in a circular manner) are taken from the mutant vector. The probability of replacing the nth element in the sequence {1, 2, ..., e}, where e ≤ N , decreases exponentially with the increase of n. The pseudo-code for the exponential crossover is shown in Algorithm 2, where j = 1, 2, . . . , N ; i = 1, 2, . . . , P S; P S is the population size; CR is the crossover possibility in the range of [0,1], rand is a uniform random number generator within [0,1] [56]. Algorithm 2 Exponential crossover pseudo-code 1: i ← 1 2: for i from 1 to P S do 3: j ←randomly chosen from [1, N ] → − − 4: u ij,g ← → x ij,g 5: e ←1 6: while rand[0, 1] ≤ CR and e ≤ N do → − − 7: u ij,g ← → v ij,g 8: j ← (j + 1) mod N 9: e←e+1 10: end while 11: end for

2.2.3. Selection operation In a selection operation, a comparison of each trial vector and its corresponding target vector is conducted, based on the fitness values (f (...)). This process

6

Figure 1: Example of 3-Opt modifications [38]

is used to determine whether either one can survive to the next generation (g + 1), by using the greedy selection strategy shown in equation 5:  → − → − → −    u i,g+1 , f ( u i,g+1 ) ≤ f ( x i,g ) → − (5) x i,g+1 =   − → x i,g , otherwise As the processes of the evolutionary operations (i.e., mutation, crossover and selection) continue as long as the averaged cost function is more than a predefined cost value (objective function), the point of termination of the DE algorithm can be determined by the maximum allowable value of the averaged cost function. 2.3. Local searches 2.3.1. 3-Opt technique In TSPs, each solution represents a tour (T ) and each T has a fitness value, which is the total distance to be traveled to complete this tour. In a local opt search method, a notion of a neighborhood is firstly defined on the set of all tours. For instance, the neighborhood of a tour (T ) can be identified as all those tours which can be obtained by changing at most K edges of T . For a specific tour if no tour in its neighborhood is shorter than this tour, this means that no improvements can be achieved and this tour is considered to be locally optimal. The 3-Opt algorithm is a special case of the K-opt algorithm proposed in [24]. In 3-Opt, three connections/edges in a network/tour are deleted and the tour is reconnected in all other possible ways, and then each potential re-connection is evaluated in order to find the optimum connection that achieves the minimum total distance. A different set of three connections are then selected and the same process is repeated. An example of 3-Opt is given in Figure 1. 2.3.2. Double bridge technique The double bridge technique is another case of the K-opt algorithm, where K= 4. It is made by perturbing the current local optimum with the same 7

process of 3-opt technique, but with a slightly larger number of edges to be disconnected and reconnected than 3-opt. Expanding the number of edges to be disconnected and reconnected aims to reduce the chances of becoming stuck in a local optimum. An example that illustrates the process of the standard double bridge technique [54] is provided in Figure 2.

Figure 2: Double bridge move [54]

2.4. Mapping methods 2.4.1. Ranked-order value (ROV) As the objective function of TSP requires discrete numbers to represent the cities for measuring the total distance between them (fitness value), the ranked-order value (ROV) [36] mapping method that is based on the random key representation [10] is applied. In ROV, the rank of each continuous value of an individual represents an index of a city. Then, the continuous values are sorted in ascending order, so that the smallest continuous value of the individual is assigned the smallest rank value of one, and the second smallest position value has a rank value of two, and so on. An example of using ROV to map the continuous values of an individual to discrete variables is given in Table 1. Table 1: Ranked-order value (ROV) mapping method example

Order/Index j

1

2

3

4

5

6

Continuous values

2.54

1.52

2.13

3.48

2.37

2.87

Ranked-order values

4

1

2

6

3

5

2.4.2. Best-matched value (BMV) The produced individuals from the DE mutation are used to be vectors with continuous variables. Recently, a new mapping method, called Best-matched value (BMV) mapping method has been proposed to be used with DE algorithm when solving discrete optimization problems [6]. BMV maps continuous values of each individual to discrete variables, and at the same time, directs the produced discrete individual towards optimality by using the best solution in each generation as a guide. Therefore, the exploitation process of the DE can

8

be improved by further exploring the current solution. The main steps of BMV is shown in Algorithm 3. Algorithm 3 Main steps of BMV mapping method 1: BMV ranks the continuous variables, of the generated individuals from the DE mutation, from the smallest to largest. 2: An integer vector is produced by replacing each continuous value by its ranked discrete index. 3: According to a pre-defined probability, some individuals are randomly chosen. 4: According to another probability, the values of some genes from the selected individuals are set to be same as the corresponding values of the best individual gene.

3. Proposed method In this section, a detailed description of our novel discrete DE (NDDE) is shown. The design of the NDDE’s components, namely: mapping mechanism, repairing method, ensemble of mutations, local searches and other components are given and explained in the following subsections. Moreover, the general framework of the proposed NDDE is provided in Figure 3 and the pseudo-code of it’s implementation is given in Algorithm 4.

M axF it cf e PS N DP op CP op BMV ROV f (i) ClusP rob cp g → − x i,g F CR rand() AE%

Table 2: Nomenclature for proposed NDDE Allowed maximum number of calls to the fitness evaluation function Current number of calls to the fitness evaluation function Population size Problem’s dimension Discrete population of solutions Continuous population of solutions Best-matched value mapping method Ranked-order value mapping method Fitness value of vector i Probability of applying proposed clustering technique to solutions in the population Number of clusters/groups used in k-means repairing method Index of current generation Current vector i in generation g Mutation parameter Crossover possibility Uniform random number [0,1] Average error of a solution from the optimal one

3.1. Individuals representation and repairing method Initially, every individual in the population is represented by a vector of real numbers (continuous) in the range [0, 1], where the length of each vector equals the number of cities+1 in the processed TSP. After that, the mapping technique 9

Start Ensemble of mutations Generate initial population with continuous variables No

DE/rand/1

DE/best/2

Ranked-order value mapping to obtain discrete population

Trigonometric mutation

Exponential crossover Apply k-means clustering for improving initial population

BMV mapping to obtain discrete population

Fitness evaluation Fitness re-evaluation Remapping modified discrete population to continuous one

Pairwise selection Calc. SuccRate of each mutation using Equ. 10

Satisfy stop criteria

Same best fitness repeated in 10 generation

Yes

No

Apply double bridge LS to all individuals with same fitness

Apply 3-Opt LS to the best solution found

Yes No

Mod(cycle, 2)=0

Yes Calc. the average SuccRate of each mutation and select the one with the highest value

Stop

DE’s basic step

NDDE’s additional step

Figure 3: Framework of the novel discrete differential evolution (NDDE)

10

Algorithm 4 Pseudo-code of the framework of NDDE 1: Read data file with the coordinates of the cities for each TSP 2: Set all NDDE’s parameters (subsection 4.1) 3: Generate initial random population with P S individuals with continuous variables (CP op) (subsection 3.1) 4: Map continuous variables of CP op to discrete ones using ROV to generate a new discrete population (DP op) (subsection 3.1.1) 5: Repair selected solutions from DP op using the k-means clustering repairing method (subsection 3.1.2) 6: Evaluate the fitness values f (...) of solutions in DP op and rank them according to their f (...) (subsection 3.1.3) 7: cf e ← cf e + P S 8: Re-map DP op from discrete to continuous variables and update the CP op in order to maintain the original continuous flow of DE (subsection 3.1.4) 9: while cf e ≤ M axF it or no improvement occurs after 100 iterations (subsection 3.7) do 10: for each individual in CP op do 11: apply ensemble of mutations (subsection 3.2) 12: apply crossover (subsection 3.3) 13: Map the new CP op (N ewCP op) using BMV mapping method to update DP op (N ewDP op) (subsection 3.4) 14: Re-evaluate the fitness of the N ewDP op’s solutions (subsection 3.4) 15: cf e ← cf e + P S 16: apply a pair-wise selection to compare the f (N ewDP op) and f (DP op) (subsection 3.5) 17: DP op ←solutions with lower distances (fitness values) 18: Rank DP op according to the fitness values of the solution with first solution is the best one 19: Apply bi-local searches (subsection 3.6) 20: end for 21: end while

11

is applied to map the continuous population of the solutions (CP op) to a discrete population (DP op), by converting continuous variables to discrete vectors. Then, the k-means clustering technique is applied, as a repairing method, for enhancing the quality of the individuals in the initial population. The detailed processes of these components are shown as follows: 3.1.1. Mapping method In order to deal with the continuous variables of the generated individuals, the ROV mapping method is applied. The explanation of ROV is illustrated in sub-section 2.4.1. The process of ROV is repeated until all the continuous values are converted to a series/path of discrete cities. In order to maintain the computational time, the mapping process of each individual is executed immediately after the initialization of it. So, by the end of the initialization process of the initial population, a continuous and its discrete populations are generated. 3.1.2. k-means clustering technique After the initial population is generated and mapped to discrete vectors, k-means clustering is applied to improve the quality of some solutions of the generated population. It is used according to a predefined percentage of the population size (ClusP rob). Number of points (cp), which can be calculated by equation 6, are generated over the geometric problem space. These points are considered the centroid points of number of groups (equal to cp). The main idea here is to first divide the problem into cp sub-problems (groups), and then solve each one to find the optimal sub-tour in each group. After that, recombining all sub-tours to form a (near-)optimal completed tour. In order to do that, cp locations, which can be calculated by equation 6, are generated over the geometric problem space. √ cp = round( N + 0.5)

(6)

where N is the number of the cities and round is a function that rounds the resulted number to nearest integer. A number of groups that equal to cp, are formed with those location points as the centroids. Each group attracts the points (cities), which are closest to its centroid point. After all the points are attached to the groups, a greedy path is formed between those points to find the shortest path in each group. The main process of the greedy algorithm that is used in each group is as follows: at each step of the journey, visit the nearest unvisited city. After that, the distances between the cp centroid points are calculated, and the two groups with the smallest distance between their centroid points are then merged to form one larger group. The merge between groups is done by joining two points/cities of a group with their closest two cities of the other group after breaking their old connections with their other local cities. The merge process between two groups is illustrated in Figure 4. Then, a centroid point for the new merged group is calculated, the distances between it and the other centroids are measured and

12

other two groups are selected to be merged. This process is repeated until only one group, containing all the cities (the completed tour) with a smaller total distance, is formed.

Group A

Group B

Group A

Group B

Group C

Figure 4: Merger process between two sub-paths

To summarise the process of the proposed repairing method, Figure 5 shows the steps of the k-means clustering and how it can be applied in the problem’s search space. In this figure, sub-figure (a) provides the distribution of the generated cp centroid points over the search space, which contains 50 candidate solutions. Sub-figure (b) shows the cp(=7) groups that are formed, the centroid of each group, and the shortest path between the cities in each group. Lastly, sub-figure (c) explains the first merger process between groups 1 and 3 that have the smallest distance between their centroids. In order to judge its effect on improving the initial solution, Figure 6 compares the total distances provided by a random individual from the initial population before and after applying the k-means clustering repairing method. From this figure, we can notice that the required distance value (fitness value) of the “before solution” to visit all the cities is 1,668, whereas the “after solution” is only 613. This suggests that the proposed repairing method can improve the quality of the solutions in the initial population by 63.25% (calculated as ( 1668−613 1668 ) × 100). The pseudo-code of the proposed k-means clustering repairing method is shown in Algorithm 5.

13

(a) Distribution of the centroid points

(b) The shortest paths in each group

3

1

(c) Merging two subgroups Figure 5: Steps of the k-means clustering repairing method

Figure 6: Solution before and after applying the repairing method

14

Algorithm 5 Pseudo-code of the k-means clustering repairing method 1: for i from 1 to P S do 2: if rand ≤ClusP rob then 3: Calculate the cp centroid points by equation 6 4: for j from 1 to cp (for each group) do 5: Add all cities that the Euclidean distance between their coordinates are very close to the coordinate of j th cp (apply k-means function) 6: Find the shortest path between the added cities in j th group 7: end for 8: Calculate the Euclidean distance between the coordinates of all the centroid points (cp) 9: while cp > 1 do 10: merge the two groups that have the smallest distance value between them by joining the closest two cities of each and breaking their current links with their local cities 11: Calculate the centroid point of the new group 12: cp ← cp − 1 13: Calculate the Euclidean distance between the coordinates of all the centroid points (cp) 14: end while 15: end if 16: end for 3.1.3. Fitness evaluation and sorting In this step, the quality of each individual in the discrete population is assessed by measuring the fitness values, which can be calculated using equation 1. As the objective function of TSP is to minimize the total distance between the visited cities, the individual with the lowest distance value is considered to be the best solution. After that, the individuals are sorted according to their fitness values, where the best solution is ranked first and the individual with the largest total distance (fitness) value last. 3.1.4. Backward mapping method In order to continue with the continuous flow of the standard DE, the discrete solutions resulted from the repairing method are mapped back to their equivalent continuous values by using the backward mapping method. This proposed method reverses the work of the ROV mapping method, by generating a random continuous vector in ascending order, where the first value of the vector is the smallest value and the last is the largest, for each discrete solution. This generated vector is sorted to have the same order as the corresponding discrete solution and so is the equivalent continuous solution of the discrete one. The backward mapping method is illustrated in Algorithm 6, and an example is given in Table 3. In this example, the first and second values of the final continuous vector are marked as bold and italic respectively, to better illustrate

15

the processes of the method. Once again to reduce the computational time, the backward mapping of each individual is executed in individual’s fitness calculation step. Then, both discrete and continuous populations are sorted according to the fitness values of the individuals. Having the continuous population, consistent with the discrete population, enables the natural DE’s processes (such as mutation and crossover) to continue. Table 3: Backward mapping method example (the first and second values of the final continuous vector are marked as bold and italic respectively)

Index j

1

2

3

4

5

6

Random continuous

0.45

0.91

1.29

1.54

1.96

2.35

Discrete solution

4

1

2

6

3

5

Ordered continuous based on discrete solution

1.54

0.45

0.91

2.35

1.29

1.96

Algorithm 6 Pseudo-code of the backward mapping method 1: ConV ect(1)← rand[0, 1] 2: for j from 2 to N do 3: ConV ect(j)← ConV ect(j − 1) + rand[0, 1] 4: end for 5: Sort ConV ect accorrding to the indices of the corresponding discrete solution

3.2. Ensemble of mutation strategies In DE, a mutation operator is used to increase the diversity of the population and thus enable DE to explore promising areas of the search space. Many research works, which use the ensemble of mutation can be found in the literature [59, 18]. However, we could not find works that used this technique for solving TSPs. In the proposed DE algorithm, three mutation strategies have been tested for solving TSP instances with different dimensions. The three mutations are as follows: (1) DE/rand/1, which uses the current vector and the difference between two other random vectors to generate the mutant vectors; (2) DE/best/2, where the mutant vector of any generation can be produced by incorporating the best target vector of that generation and the differences between other four vectors (instead of two in equation 3), randomly selected from the population, as shown in equation 8 ; and (3) Trigonometric mutation [24], which uses the fact that DE is performed usually on the basis of three individuals in the current population and that an individual is disturbed with a scaled vector differential from the other two individuals to produce a mutated vector/individual. Consequently, a hyper-geometric triangle is formed in the search space where the three selected individuals (i.e., triangle head points) exist as the vertices. Then, the selected-individuals’ objective function values together with this triangle can be

16

used in the mutation operation. Applying the objective function value information biases the process mutation operation towards the points (the vertices of the hyper-geometric triangle) with the best/lowest objective function values, as in equation 9. The equations of the three mutation strategies are given as follows: (1) DE/rand/1 → − − − − v i,g+1 = → x r1 ,g + F × (→ x r2 ,g − → x r3 ,g ), r1 6= r2 6= r3

(7)

(2) DE/best/2 → − − − − − − v i,g+1 = → x best,g + F × (→ x r1 ,g − → x r2 ,g ) + F × (→ x r3 ,g − → x r4 ,g )

(8)

(3) Trigonometric mutation → − − − − v i,g+1 = (→ x r1 ,g + → x r2 ,g + → x r3 ,g )/3 → − → − +(p2 − p1 )( x r1 ,g − x r2 g ) − − +(p3 − p2 )(→ x r2 ,g − → x r3 ,g ) → − → − +(p1 − p3 )( x r3 ,g − x r1 ,g ) where − p1 = |f (→ x r1 ,g )|/´ p → p2 = |f (− x r2 ,g )|/´ p − p3 = |f (→ x r3 ,g )|/´ p

(9)

and − − − p´ = |f (→ x r1 ,g )| + |f (→ x r2 ,g )| + |f (→ x r3 ,g )| → − → − → − x r1 ,g , x r2 ,g , x r3 ,g are three randomly chosen vectors, where r1 6= r2 6= r3 , − − are not equal to either each other or the target vector (→ x i,g ), and f (→ x r,g ) is the → − fitness value of individual x in generation g. Trigonometric mutation is applied to generate a new mutant vector by using its original vector (in the current generation) as a center point of a hyper-geometric triangle. A modification, which is summed three weighted vector differentials and uses their fitness values to bias the new vector towards a point with better objective value, is then applied to this vector. The results showed that the mutation strategies based on the best solution, i.e. DE/best/2, has a fast convergence speed and performs well on large-scale problems. However, it is more likely to get stuck at local optimum when solving small problems, as the exploitation of the best solution-based strategy has an effect on population diversity [7]. On the contrary, the mutation operators with randomly selected vectors, i.e. DE/rand/1, could enhance the DE exploration

17

process, as a high degree of variability affects the mutation and hence population diversity. Thus, this random-based mutation performed well on small-scale problems, where there are small chances of finding new solutions. Trigonometric mutation enables DE to get a better trade-off between the convergence rate and robustness and thereby obtain an acceptable solution. This mutation showed a balanced performance and was best for the medium-scaled instances. Algorithm 7 Pseudo-code of the ensemble of mutation strategies 1: M axIter ← M axF it/P S 2: iter ← 1 3: CycIter ← 0 4: M utationP ar ← 0 5: while iter<=M axIter do 6: iter ← iter + 1 7: if M utationP ar equal 0 // explore P S using the three mutations then 8: for t from 1 to 3 do 9: P Smt ← round(P S/3) 10: for s from 1 to P Smt do 11: Select non taken random solution from P S 12: end for 13: Apply mt mutation to P Smt 14: end for 15: else 16: Apply mutation m with the highest success rate to the whole P S 17: end if 18: if mod(CycIter,CS) equal 0 and CycIter < 2 × CS //Calculate the success rates of the three mutations then 19: CycIter ← CycIter + 1 20: for y from 1 to 3 do 21: Calculate the success rate (Succratem ) for every my mutation using equation 10 22: end for 23: M utationP ar ← my with the highest success rate value 24: else if CycIter > 3 × CS then 25: M utationP ar ← 0 26: CycIter ← 0 27: end if 28: end while In order to get the most benefits from the DE mutation step, an ensemble of the three mutation strategies is proposed, to efficiently explore the candidate solutions in the population. The pseudo-code for the implementation of the ensemble of mutations, is provided in Algorithm 7. In this technique, during the first generation, the population is divided into three sub-populations. Each sub-population is explored by one of the three mentioned mutations. At the end of each generation, the success rates (Succrate ) of the three mutations are 18

calculated, such that, Succratem,g =

sim,g P Sm,g

∀m = 1, 2, or 3

(10)

where sim,g is the total number of individuals successfully improved by the m mutation strategy at the current generation (g). P Sm,g is the sub-population size assigned to m mutation in g generation. After a predefined number of generations (CS), called cycle, the performances of the three mutations are evaluated by calculating their average success rate per cycle. For the next cycle, individuals use the mutation with the highest average Succrate , and the two worst individuals of that mutation are replaced by the best individual of each of the other mutations. At the end of next cycle, the three mutations are re-run in parallel again, and the Succrate of each is re-calculated for the new cycle. The trade-off between the three mutation strategies in every cycle is shown in Figure 7. Also, the process of the proposed ensemble mutation is shown in the discrete framework of DE in Figure 3.

Figure 7: The trade-off between mutation strategies in every cycle

3.3. Crossover In DE, trial vectors are generated by combining the target and mutant vectors according to a predefined possibility (CR). As the resultant individuals from the mutation are still in continuous variables, the exponential crossover is used to further explore the search space. 3.4. Mapping method and fitness re-evaluation In order to measure the fitness values of the individuals in the new population, the individuals’ variables need to be converted from their current continuous format to discrete format (i.e., cities). To do that, the best-matched value mapping method (BMV) [5, 6] is extended, to map continuous variables to discrete, not binary, variables. Then, the fitness values of the new discrete individuals are calculated. Finally, the individuals are sorted according to their fitness, where the solution with the best fitness is located first. 19

3.5. Selection A pairwise selection between two corresponding individuals (one from the new population and its old version in the previous population) is conducted, according to the greedy selection operator mentioned in equation 5. At the same time, an archive is created to save the improvement rate of each solution in the new population against its corresponding version in the old population. mutation strategy achieved it, and hence calculate the Succrate of each mutation at the end of the current cycle. In the next cycle, the mutation strategy with the highest Succrate value will be selected to perturb/explore the population. 3.6. Local searches In this subsection, after improving the exploration capability of the proposed DE by applying the ensemble of mutation strategies in subsection 3.2, local searches are used to enhance the exploitation capability of DE. In NDDE, two local searches are implemented and applied in two positions, as follows: 1) 3-opt local search is applied to the best individual, immediately after the new population is created by the selection operator. In its process, the algorithm first does 2-opt and falls back to 3-opt when the 2-opt opportunities are exhausted. This can be detected if no improvements happen after 5 consecutive trials. If the algorithm reaches a local optimum after applying 3-opt, then NDDE restarts the local search from a random route. 2) The double bridge technique is applied, if more than 100 consecutive generations provide the same best fitness value (i.e. NDDE is stuck in a locally optimum solution). In this technique, list of individuals, which are sorted by fitness, is processed in order. If any individual is found, that has the same fitness of a previously processed individual, then the double bridge technique is applied to it. The pseudo-code of the local search strategies is provided in Algorithm 8. In the proposed DE, the process of the double bridge local search is as follows: two consecutive individuals are selected and the Euclidean distance between two corresponding cities/points from the two individuals are calculated. Seeking to guarantee an efficient distribution of individuals over the search space, if the distance between the two solutions is smaller than a certain threshold (α), then one of the two solutions is perturbed in such a way that makes it different than the other solution. In the modified double bridge technique, the solution is divided into 4 parts/sub-groups, then these sub-groups are reordered by grouping each sub-tour with its nearest neighborhood sub-group, according to the Euclidean distances between their centroid points. In order to maintain diversity in the early iterations, α starts with a value very close to zero (similarity is increased when α is close to zero). As the iterations progress, the value of α is gradually moved away from zero, using equation 11, to allow NDDE to converge to the (near-)optimal solution. α = αinit +

cf e M axF it

(11)

where M axF it is the maximum number of fitness evaluations, cf e is the current fitness evaluations, rand is a random number in the range [0, 1], and αinit is the 20

Algorithm 8 Pseudo-code of the local search strategies 1: M axIter ← M axF it/P S 2: iter ← 1 3: N oImprovements ← 0 4: while iter<=M axIter and N oImprovements ≤ 100 do 5: DP op ← the discrete population after the DE’s selection and ranking processes 6: BestSched ← first (best) individual in DP op 7: N ewBestSched ←apply the 3-Opt local search to the BestSched solution 8: f (N ewBestSched) ←measure the fitness value of the new solution 9: if f (N ewBestSched) < f (BestSched) then 10: Accept the new solution and its fitness as the best solution and fitness 11: BestSched ← N ewBestSched 12: DP op ← BestSched 13: f (BestSched) ← f (N ewBestSched) 14: end if 15: if iter > 1 and f (BestSched)iter == f (BestSched)iter−1 then 16: N oImprovements ← N oImprovements + 1 17: else 18: N oImprovements ← 0 19: end if 20: if N oImprovements ≥ 10 then 21: Apply the double-bridge local search to all the solution in DP op with the same fitness values 22: end if 23: end while

21

initial value of the threshold. The pseudo-code of the modified double bridge local search is given in Algorithm 9. Algorithm 9 Double bridge technique 1: for i from 1 to P S − 1 do 2: pre_sol← P op(i) 3: cur_sol← P op(i + 1) 4: Euc_Dis←Norm difference (current_sol - previous_sol) 5: if Euc_Dis ≤ α then then 6: Divide cur_sol into four segments/sub-tours 7: Calculate the centroid for each sub-tour 8: Reconnect sub-tours with the minimum distance between their centroid points 9: Repeat until the tour construction is completed 10: end if 11: end for The two original local searches are explained in subsections 2.3.1 and 2.3.2. 3.7. Termination conditions In this study, in order to save computational time, two termination conditions are set and applied simultaneously. 1) Maximum allowed number of objective function evaluations (cf e), as the algorithm is terminated if cf e > M axF it , where M axF it is the maximum fitness evaluations. 2) The performance’s progress of the algorithm, as the algorithm is terminated if no improvement occurs in 100 consecutive iterations. 4. Experimental Analyses and Results This section starts with parametric analyses to decide the best parameters values to run the algorithm with. This is followed by a summary of the results obtained by applying the proposed NDDE. After that, several comparisons are conducted for evaluating the performance of the proposed NDDE, as follows: 1) self-comparison, where NDDE is compared with 5 other versions of NDDE to assess the significance of the different components, 2) comparison with DE-based algorithms, and finally 3) comparison with other state-of-the-art algorithms. In order to judge the performance of the proposed DE algorithm, computational experiments were carried out for solving 45 instances of TSPs with different numbers of cities. The instances were selected from the well-known TSP library (TSPLIB), which is described in [53]. The proposed DE was coded in MATLAB R2017b, and was executed on a PC with an i7 Processor and 16 GB memory. In our algorithm, each instance was independently executed 30 times with (500×N ) maximum fitness evaluations (M axF it ) [58]. In order to assess the quality of the algorithms, average error (AE%) value was used. AE% presents the average error between the result obtained by an algorithm and the 22

optimal/best-known solution of a problem. AE% can be calculated by equation 12.   S X Best − Opt Sol Sol s s AE% =  × 100 /S (12) OptSols s=1 where S is the total number of instances to be solved, BestSol is the best solution achieved by an algorithm for s instance and OptSol the pre-known optimal/best known solution of instance s. 4.1. Parameters tuning and settings Three sets of experiments were designed to analyze the effects of the three main parameters in DE algorithm, which are i) crossover rate (CR); ii) mutation factor (F ); and iii) population size (P S), on the performance of the proposed DE. Moreover, another two experiments have been set to analyze iv) the number of individuals to be repaired by the k-means clustering technique as a probability of P S in the initial population (ClusP rob), and the number of iterations in each cycle (CS) in the mutation step. For all experiments, the proposed DE was executed to solve 6 TSPs with a different number of cities, for 30 runs with up to 5000 generations in each. In these experiments, the selection of parameters was done in a sequential manner, in which the best parameter found in an experiment was fixed in the subsequent experiments and the other parameters, which not analyzed yet, are initialized with arbitrarily chosen initial values. For particular, while analyzing CR, the other parameters F , P S, ClusP rob and CS are set to the initial values of 1.5, 50, 50% and 15, respectively. 4.1.1. The influence of the DE parameters: crossover (CR), mutation factor (F ) and population size (P S) In the first set of experiments, the effect of the CR, F and P S operators were studied by running the proposed DE using different values of CR = [0.1, 0.3, 0.5, 0.7, 0.9]; F = [0.5, 1, 1.5, 2], Adapt[0.8-0.6], and Adapt[2-0.5] which were linearly decreased from 0.8 to 0.6 and from 2 to 0.5 respectively, according to equation 13; and P S=[25, 50, 100, 200].   M inF − M axF F = × cf e + M axF (13) M axF it where M inF and M axF are the lower and upper values of F . To provide an indication about the effect of the CR parameter, the AE% values that are resulted from solving 6 TSPs by using different values of CR is shown in Figure 8. The figure shows the AE% values as bars and the standard error with 95% confidence interval as small horizontal bars. From the figure, it can be noticed that CR=0.5 had the lowest AE% values for all tested problems, with significant difference in problem “pr76”. So, applying the exponential crossover to 50% of the individuals is used from here on.

23

Figure 8: Average error (AE%) of 6 TSPs with different values of crossover rate (CR)

Another variants of NDDE were run while varying F and P S as described above. Table 4 gives the mean AE% values of 6 TSPs solved using different values of F and P S. The computational results given in Table 4 shows that F , with value 1, is more suitable for solving problems with small size. However, F with an adaptive value, which decreased from 0.8 to 0.6 performed well with larger problems. The average values over the 6 TSPs suggested that F =Adapt[0.8-0.6] is the best value for the mutation parameter. Moreover, from the results, it can be noticed that increasing the value of P S to more than 50 had a negative influence on the performance of the algorithm. At the same time, a smaller population size, such as P S=25 reduces diversity, and hence degrades the performance of the algorithm due to becoming stuck in a locally optimum solution. Therefore, the AE% results show that the variant with P S=50 was the best for most of the solved TSPs. 4.1.2. The influence of the probability of individuals’ repairing (ClusP rob) and the number of iterations in each cycle (CS) parameters In order to obtain the best performance of the applied repairing method (k-means clustering), experiments, which are applying the k-means clustering to (0%, 10%, 30%, 50%, 70%, 90%, 100%) of the individuals in the initial population were designed and executed. Based on Table 5, there is no significant difference between AE% of the different values of ClusP rob. Thus, Friedman rank test was used to rank the different variants of NDDE based on their qualities. According to the Friedman rank, it was found that applying the repairing method for 10% of individuals (ClusP rob =10%) is the best. The results demonstrated the importance of using the repairing method over the initial population, as applying it for 10% of individuals can achieve better mean AE% value than

24

Table 4: AE% values of 6 TSPs using different values of mutation rate (F ) and population size (P S) F

TSPs

PS

0.5

1

1.5

2

Adapt[0.8-0.6]

Adapt[2-0.5]

25

50

100

200

eil51

2.12

2.07

2.14

2.06

2.13

2.12

2.36

1.96

2.13

2.64

berlin52

2.50

2.15

2.22

2.45

2.16

2.46

2.41

1.98

2.38

3.57

st70

2.57

2.43

2.39

2.66

2.32

2.37

2.52

2.45

2.43

3.71

pr76

2.23

2.05

1.99

1.90

1.90

2.01

2.92

1.99

2.29

2.79

kroA100

2.75

2.66

2.60

2.69

2.51

2.52

3.22

2.58

3.52

4.87

lin318

6.46

6.38

6.38

6.41

6.30

6.41

6.93

6.36

7.30

8.71

Avg.

3.11

2.96

2.95

3.03

2.88

2.98

3.39

2.89

3.34

4.38

Table 5: Friedman test for different values of the probability of individuals’ repairing (ClusP rob)

ClusP rob

Mean AE%

STD

Friedman mean

Rank

0% 10% 30% 50% 70% 90% 100%

4.26 4.06 4.13 4.17 4.32 4.20 4.44

1.08 1.16 1.00 1.11 1.05 1.07 1.06

4.33 2.83 3 3 3.5 5.5 5.83

5 1 2 3 4 6 7

applying it to 0% of individuals. However, it can be noticed that increasing the number of individuals to be repaired can reduce the efficiency of the algorithm, as the diversity between the individuals within the initial population is decreased. For determining the best number of iterations for each cycle, the algorithm was run using cycle (CS) equals to {5, 10, 15, 20, 25} iterations, as shown in Figure 9. From the figure, it can be concluded that there is no certain number of iterations, which is good for all problems, as the standard error bars didn’t show any significant difference between all CS’s values. Therefore, the number of iterations per cycle (CS) is designed to be selected adaptively by the algorithm itself. So, for the first 100 iterations, the algorithm was initially run with all possible CS values, with each value used for an equally sized sub-population. At the end of the initial iterations, the success rates of each value of CS are calculated and based on it, the CS value with the highest success rate is selected to continue solving the current problem through the coming iterations. 4.1.3. Parameters settings Based on the above experiments, the best parameters setup for the proposed NDDE are found to be as shown in Table 6. 25

Figure 9: AE% for 7 traveling salesman problems (TSPs) with various scales using different numbers of iterations in each cycle (CS)

Table 6: Best parameters settings of NDDE General parameters DE’s basic parameters NDDE’s parameters

M axF it Runs CR F PS ClusP rob CS

500×N 30 0.5 Adapt[0.8-0.6] 50 10% Selected adaptively based on the performance of NDDE for each TSP

4.2. Results and comparisons In this subsection, the results obtained from solving 45 different instances of TSPs with dimensions between 51 to 15112 number of cities are examined. According to the results, the proposed NDDE was able to achieve the optimal solutions for 80% (36 out of 45 problems) of the used data-set. Moreover, several self-comparisons between different versions of DE, in order to show the effect of each component on the overall performance of the proposed DE algorithm, are conducted. Eventually, the performance of the proposed NDDE is compared with those from state-of-the-art algorithms. To maintain a fair comparison, some earlier versions of the proposed DE, in addition to the final one, will be included to be compared with those which have the same components (i.e., proposed DE without local search V.S. standard DE without local search and so on). 4.2.1. Statistical analyses In order to statistically study the difference between the algorithms in this paper, statistical significance testing was performed using the Wilcoxon signed rank test [17]. This test can be used to judge the difference between paired 26

scores when the population cannot be assumed to be normally distributed, as an alternative to the paired-samples t-test. As a null hypothesis, it was assumed that there was no significant difference between the mean values of two samples for the number of the instances of TSPs with a 95% confidence level. 4.2.2. Self-comparisons Several versions of DE, starting from the standard DE with ranked-order value (ROV) mapping method to the final proposed version of DE, which incorporates a repairing method, an ensemble of mutation strategies, and local searches, are described in this sub-section. Six versions of DE were designed and are illustrated in Table 7. Table 7: Description of DE versions

Version

Description

1

Standard DE with ROV mapping

2

Standard DE with the proposed mapping method, Best-matched value (BMV)

3

DE with the adaptive ensemble of mutation operators and BMV mapping

4

Version3 + Modified k-means clustering as a repairing method

5

Version4 + Double-bridge local search

6

Version4 + 3-Opt local search (Final Version)

The significant difference between the 6 versions, based on the mean and best error values (AE%) from the optimal solution, are summarized in Table 8. The results presented in the table demonstrate the efficiency of each component in enhancing the performance of the standard DE algorithm. This is because each added component significantly improved the performance of the algorithm, by getting lower best and mean AE% values and solving more problems in comparison with its previous version. The statistical test, based on the best AE%, did not show a significant difference between versions 2 and 3 of DE, and this is probably because the ensemble of mutations focuses on enhancing the exploration rather than the exploitation capability of the proposed DE. However, version 3 significantly enhanced the mean AE% and solved more problems than version 2. 4.2.3. Comparison with DE-based algorithms In this subsection, the proposed NDDE will be compared with different DE-based algorithms with and without applying 3-Opt as a local search. For providing a fair comparison, two versions of NDDE were run. One version used double bridge local search and the other used 3-Opt. DE algorithms without 3-Opt local search: the performance of the proposed NDDE is compared with two DE algorithms from the literature: 1) hybrid differential evolution algorithm for traveling salesman problem (HDE) [57] and 2) an improved differential evolution algorithm for TSP (Improved DE) [40],

27

Table 8: Wilcoxon signed ranks test based on the mean and best error values with a 95% confidence level Mean AE% Best AE% DE Versions Better Equal Worse P Value Better Equal Worse P Value Version2 - Version1

23

0

7

0.001

26

0

4

0.000

Version3 - Version2

21

0

9

0.010

16

0

14

0.441

Version4 - Version3

30

0

0

0.000

30

0

0

0.000

Version5 - Version4

30

0

0

0.000

29

1

0

0.000

Version6 - Version5

30

0

0

0.000

29

1

0

0.000

which incorporates hill-climbing algorithm and Liuhai crossover, respectively as local searches. In order to make fair comparison, the double bridge technique was applied as a local search in the proposed DE. The best and mean distances of 5 small size TSPs using NDDE and the two DE algorithms are shown in Table 9. The results show that DE can obtain much better best and mean distances than the other two algorithms and demonstrate its capability for solving TSPs with small size. Table 9: Comparison of the proposed NDDE with other DE algorithms using the best and mean distances of five TSPs

Prob.

NDDE

Opt.

HDE

Improved DE

Best

Mean

Best

Mean

Best

eil51

426

428

432.5

439

444

429

berlin52

7542

7542

7679

7544

7816

-

st70

675

676

688

684

691

-

eil76

538

543

555.27

558

568

545

pr76

108159

108917

110481.2

109491

110539

-

DE algorithms with 3-Opt local search: the proposed DE is also compared with another DE algorithm called “A Real Adjacency Matrix-Coded Differential Evolution Algorithm for Traveling Salesman Problems” (RAMDE), which besides the real adjacency matrix-coding mechanism, also uses 3-Opt local search for further improving the results of TSPs. From the differences between the best and average error percentages (Avg. Error(%) − B. Error(%)) in Table 10, it can be noticed that RAMDE has more consistent performance than the proposed DE. However, RAMDE seems to be stuck in local optima and could not achieve any optimal solution for any problems. On the other hand, the proposed NDDE with the ensemble of mutations can explore the search space more efficiently than RAMDE and could reach the optimal solutions for 83% (5/6) of solved TSPs.

28

Table 10: Best and average distances and errors from the optimal solution for 6 large-scale TSPs, where bold means lower distances and better errors

Instance

Best known

lin318

42029

pcb442

d493

d657

u724

fl1400

50778

35002

48912

41910

20127

Algorithm

Best

B. Error (%)

Average

Avg. Error (%)

RAMDE

42513

1.15

42513

1.15

NDDE

42029

0

42280

0.6

RAMDE

51396

1.22

51532.9

1.49

NDDE

50778

0

51142.3

0.72

RAMDE

35642

1.83

35717.8

2.05

NDDE

35002

0

35189.6

0.54

RAMDE

49953

2.13

50304.2

2.84

NDDE

48912

0

49376.6

0.95

RAMDE

42885

2.33

43025.4

2.66

NDDE

41910

0

42376.6

0.6

RAMDE

20455

1.63

20455

1.63

NDDE

20164

0.18

20306.8

0.89

4.2.4. Comparison with other state-of-the-art algorithms In this subsection, a comparison between NDDE and other top meta-heuristic algorithms, namely: a set-based particle swarm optimization (S-CLPSO) [14], and ant colony system (ACS) [21], is conducted and the results are shown in Table 11. In the table, the best and mean distances, average error from the mean, and the computational time in millisecond for 6 TSPs with large size are shown. The results demonstrated the superiority of the proposed NDDE for solving large instances in terms of average error from the best known solutions, and showed that NDDE can save on average 84.6% and 73.5% of computational time, which are required by ACS and S-CLPSO respectively, for solving the same problems. Further comparisons between the proposed NDDE and three versions of SCLPSO, which adapted the 3-Opt local search, were conducted and the results are shown in Table 12. The three versions were described in [14] as they are all using 3-Opt local search but in different positions. So, in the versions “gbest”, “iteration best”, and “all positions”, the 3-Opt local search were performed on the global best position, iteration best position and on all positions, respectively. The results in Table 12 show that NDDE could get better best and mean AE% values than gbest version by 97.56% and 71.56%, respectively. Also, NDDE was able to find best AE% values equal to the other two versions (except for Fl1400) in much lower computational time. Although the “iteration best” and “all positions” versions could achieve better mean AE% values, the “iteration best” consumed computational time 37.68% more than the time required for NDDE to solve the same set of problems. It is also expected that applying the 3-Opt to both global and iteration positions (i.e. the “All Position” version) requires much more computational time. 29

Table 11: Comparison between NDDE and other top-algorithms, where bold means lower distances and better errors

Instance

Best known

lin318

42029

pcb442

Algorithm

Best

Mean

Avg. Error

Time (ms)

NDDE

42029

42280

0.60%

12580

ACS

48739

50690.38

20.61%

383877

S-CLPSO

42719

43518.4

3.54%

314635

NDDE

50778

51142.3

0.72%

16600

ACS

61879

66013.6

30.00%

-

S-CLPSO

-

51577.7

1.57%

-

NDDE

35002

35189.6

0.54%

14740

ACS

39050

40271.4

15.05%

156588

S-CLPSO

-

37028.83

5.79%

78751

NDDE

48912

49376.6

0.95%

100590

ACS

58405

59851.86

22.37%

299591

S-CLPSO

-

51799.53

5.90%

167444

NDDE

41910

42376.6

0.60%

31320

ACS

58405

59704.8

42.46%

342648

S-CLPSO

-

43373.3

3.49%

188186

NDDE

20164

20306.8

0.89%

470980

ACS

29376

30468

51.38%

2201702

S-CLPSO

22141

22025.2

9.43%

1419804

50778

d493

35002

d657

48912

u724

41910

fl1400

20127

Table 12: Comparison between NDDE and other versions of S-CLPSO with 3-Opt local search Instance

Best known

gbest (AE%)

Iteration (AE%)

All Positions (AE%)

NDDE (AE%)

Best

Mean

Best

Mean

Time (ms)

Best

Mean

Best

Mean

Time (ms)

D198

15780

0.13

1.48

0

0

9141

0

0

0

0.42

4800

Lin318

42029

0.76

1.98

0

0.24

24558

0

0

0

0.6

12580

D493

35002

1.04

3.1

0

0.21

63823

0

0

0

0.54

14740

D657

48912

2.51

3.71

0

0.09

117928

0

0.01

0

0.95

100590

U724

41910

1.69

2.53

0

0.1

133549

0

0.01

0

0.6

31320

D1291

50801

0.75

1.88

0

0.11

501289

0

0.01

0

1.1

98110

Fl1400

20127

1.1

4.02

0

0.05

757030

0

0.02

0.18

0.89

470980

Fl1577

22249

0.93

2.57

0.02

0.03

693490

0.02

0.02

0.02

0.95

700820

30

Figure 10: Compare the best, average AE% and computational time between NDDE and other 11 GA-based algorithms

Moreover, comparison between the NDDE and some genetic- and evolutionarybased algorithms, namely 1) Genetic Algorithm (GA) [27]; 2) Island based distributed GA (IDGA) [2]; 3) Evolutionary simulated annealing (ESA) [62] with the results have been reported from [47] is given in Table 13. Also, the comparison includes some other discrete algorithms: improved bat algorithm (IBA) [47]; Discrete firefly algorithm (DFA) and a discrete imperialist competitive algorithm (DICA), which were implemented by [47]. Table 14 provides the results from this comparison. For demonstrating the performance of NDDE, large TSPs instances with cities from 5,934 up to 15,112, have been solved and compared with two recent algorithms, namely: Discrete pigeon-inspired optimization algorithm with Metropolis acceptance criterion for large-scale traveling salesman problem (DPIO) [67] and Hybrid discrete artificial bee colony algorithm with threshold acceptance criterion for traveling salesman problem (HDABC) [65]. The results, which are given in Table 15, show the efficiency of the proposed NDDE for solving large-sized TSPs, as it could achieve better AE% values than both DPIO and HDABC. Noting that, the variations in the computational time back to the used programming language in the implementation, as DPIO used C++ and HDABC was applied by Java, and both of them are faster than MATLAB. More comparisons with several variants of GA-based algorithms proposed in [20, 3] are finally shown in Figure 10. Lastly, the Friedman rank test is used for ordering the compared algorithms according to their performance in solving TSPs. The mean AE%, average time in seconds, Friedman mean rank values, and the order of each algorithm, are shown in Table 16. The results obtained from the comparisons demonstrate the outstanding performance of the proposed DE, in comparison with the performance of the other discrete algorithms, which are mainly used to solve discrete problems, such as TSPs. The results also show that NDDE achieved higher quality solutions than other discrete algorithms for all problems in much faster computational time than others, by on average 87.98%. This is thanks to the proposed repairing method, which enables DE to start with near-optimal initial 31

32

629

58,537

73,682

48,191

Eil101

Pr144

Pr152

Pr299

Avg.

21,294

KroD100

538

Eil76

21,282

675

St70

20,749

7,542

Berlin52

KroC100

426

Eil51

KroA100

Opt.

Instance

Algorithm

48,447.50

74,005.20

58,550.60

633.1

21,356.00

20,807.90

21,327.10

540.4

677.6

7,559.10

426.7

Avg.

48,191

73,682

58,537

629

21,294

20,749

21,282

538

675

7,542

426

Best

225.80

811.12

1,022.05

43.01

12.97

196.06

186.26

142.62

7.59

5.8

54.07

2.21

S. dev.

NDDE

2.93

7.3

3

3.5

3.2

2.4

2.2

2.8

2.1

2.3

1.7

1.7

Time (s)

50,532.30

74,969.50

58,807.30

658.4

21,726.50

21,170.40

21,481.70

553.7

682.1

7,542

431.6

Avg.

49,242

74,172

58,574

650

21,500

20,749

21,282

546

675

7,542

426

Best

195.15

915.8

498.9

220.9

4.4

156.9

188.7

150.1

4.2

3.9

0

2.9

S. dev.

28.04

158.7

39.5

33.9

16.3

15.9

15.4

14

5.8

4.5

2.3

2.1

Time (s)

ESA - evolutionary simulated annealing

50,817.10

75,658.30

60,591.40

673.8

22,184.60

21,510.40

49,659

74,520

58,599

655

21,492

20,861

21,350

545

675

7,542

427

Best

553.69

1585.7

910.8

2342.8

12.5

405

390.2

420.8

9.8

5.7

0

7.3

S. dev.

24.37

147.6

33.4

32.8

10.6

9.7

10.2

9.9

5.6

4.2

2.4

1.7

Time (s)

GA - Genetic Algorithm

21,812.40

565.4

709.8

7,542

440.8

Avg.

50,513.30

75,126.70

58,893.00

660.7

21,696.90

21,298.70

21,731.80

557.7

690.2

7,542

434.4

Avg.

49,572

74,249

58,581

650

21,582

20,830

21,345

545

675

7,542

426

Best

394.99

1257.9

1005.7

1012.4

7.5

408.9

290.7

340.7

6.8

9.8

0

4.5

S. dev.

24.80

149.94

32

32.5

11.7

12.1

11.2

10.7

5.1

4.1

2.3

1.2

Time (s)

IDGA - the Island based Distributed Genetic Algorithm

Table 13: Comparison between NDDE and other GA-based algorithms

33

426

7542

675

538

21,282

20,749

21,294

629

58,537

73,682

48,191

Eil51

Berlin52

St70

Eil76

KroA100

KroC100

KroD100

Eil101

Pr144

Pr152

Pr299

Avg.

Opt.

Instances

Algorithms

48,447.50

74,005.20

58,550.60

633.1

21,356

20,807.90

21,327.10

540.4

677.6

7,559.10

426.7

Avg.

48,191

73,682

58,537

629

21,294

20,749

21,282

538

675

7,542

426

Best

225.80

811.12

1,022.05

43.01

12.97

196.06

186.26

142.62

7.59

5.8

54.07

2.21

S. dev.

NDDE

2.93

7.3

3

3.5

3.2

2.4

2.2

2.8

2.1

2.3

1.7

1.7

Time (s)

49,674.10

74,676.90

58,876.20

646.4

21,593.40

21,050.00

21,445.30

548.1

679.1

7,542

428.1

Avg.

48,310

73,921

58,537

634

21,294

20,749

21,282

539

675

7,542

426

Best

214.37

1200.1

426.5

295.6

4.9

141.6

164.7

116.5

3.8

2.8

0

1.6

S. dev.

24.43

147.2

31

30.3

13.1

11.7

12

10.6

5.1

3.9

2.1

1.7

Time (s)

IBA - Improved Bat Algorithm

49,839.70

74,934.30

58,993.30

659

21,683.80

21,096.30

21,483.60

556.8

685.3

7,542

430.8

Avg.

48,579

74,033

58,546

643

21,408

20,756

21,282

543

675

7,542

426

Best

234.93

1305.4

483.7

300.1

8.1

163.7

148.3

163.7

4.9

4

0

2.3

S. dev.

24.94

149.1

32.1

30.9

13.3

12.4

12.8

10.3

5.3

4.3

2.2

1.6

Time (s)

DFA - Discrete Firefly Algorithm

Table 14: Comparison between NDDE and other discrete algorithms

49,880.30

74,886.70

59,070.90

663.8

21,666.80

21,103.90

21,500.30

557.6

684.7

7,542

432.3

Avg.

48,600

74,052

58,563

644

21,399

20,756

21,282

544

675

7,542

426

Best

253.75

1413.7

513.9

323

9.6

174

161.1

183.4

5.8

3.7

0

3.1

S. dev.

24.89

150.3

32

30.7

12

12.6

11.7

10.8

5.3

4.1

2.5

1.8

Time (s)

DICA - Discrete Imperialist Competitive Algorithm

Table 15: AE% and computational time in seconds of NDDE, DPIO and HDABC

Instance

Optimal

NDDE

DPIO

HDABC

AE%

Time (s)

AE%

Time (s)

AE%

Time (s)

Rl5934

556045

0.00

74.26

1.041

68.7

1.79

63.74

Pla7397

23260728

0.02

132.28

1.441

113.6

1.47

108.34

Usa13509

19982859

0.02

419.67

1.168

318.3

1.57

434.04

Brd14051

469388

0.01

436.81

1.051

347.7

1.67

452.87

D15112

1573084

0.01

651.72

0.984

522.1

-

-

0.012

342.948

1.137

274.08

1.625

264.75

Avg.

Table 16: Friedman ranks of NDDE and other 6 meta-heuristic algorithms

Algorithms

AE%

Average Time(s)

Mean Rank

Order

NDDE

0.33

2.93

1.55

1

ESA

1.94

28.04

3.95

4

GA

3.72

24.37

6.68

7

IDGA

2.37

24.8

5.32

6

IBA

1.28

24.43

2.23

2

DFA

1.88

24.94

3.86

3

DICA

2.00

24.89

4.41

5

34

solutions. This allows it to achieve optimal solutions quickly after a few numbers of iterations. It is also due to the modified mapping method, which maps continuous variables of solutions to discrete solutions, and at the same time, directs the produced discrete solutions towards optimality. The significant saving in computational time that was achieved by the proposed discrete version of DE, demonstrates the efficiency of the introduced components, in improving the performance and the convergence speed of DE algorithm. On the other hand, it can be noticed that the average standard error of NDDE is higher than both the ESA and IBA algorithms in Tables 13 and 14 respectively. This is because of the ensemble of mutations in NDDE, that allows it to explore and define new solutions in the search space. However, as a result of this mutation, it leads to higher variation in the population of the candidate solutions, and hence more diversity in the population. As GA has been considered one of the most powerful evolutionary algorithms for solving discrete problems, Figure 10 compares the best and average AE%, and the computational time (in millisecond) of the proposed NDDE against 11 discrete algorithms, which were proposed based on GA. The figure also confirms the efficiency of the proposed NDDE, as it achieved lower best and average error values than all of the GA-based algorithms. It is also did this in lower time, except for DISP. However, the average best error of DISP is 5.8% and its average mean error is 7.9%, against 0% and 0.9% for NDDE, and the average time difference is only 1 second, calculated as 6.5 sec. (NDDE) minus 5.5 sec. (DISP). This indicates that DISP is unlikely to achieve better results than NDDE, if it runs for the same time. 5. Conclusion and Future Work Several meta-heuristic algorithms have been introduced for solving TSPs. However, some shortcomings, such as slow convergence speed, computational time and their ability to solve large-sized discrete problems can still be identified. Also, the fact that no single mutation operator has proven to be the best or that can perform consistently well over a wide range of TSPs, makes the existence of an algorithm, which can solve these problems a necessary requirement. Consequently, a novel discrete DE (NDDE) algorithm, which incorporates useful components was proposed. NDDE uses ROV and the previously proposed mapping method (BMV) for mapping continuous variables to discrete ones. Also, the well-known k-means clustering mechanism was modified and applied as a repairing method for obtaining near-optimal solutions within the initial stage. A backward mapping method was also proposed to maintain the continuous flow of the original DE in the evolving operators. After that, an ensemble of mutation strategies has been designed for increasing the exploration capability of the DE. In order to measure the fitness values of the individuals in the newly generated population, the BMV mapping method was adopted. Finally, two local searches were applied as follows: 3-opt was applied to the best solution in each generation, and the double bridge to all individuals with the same fitness value for maintaining diversity within the population. 35

In order to assess the performance of the proposed NDDE, 45 instances of TSPs taken from the well-known TSP library have been solved. The results showed how the proposed NDDE could converge quickly towards the optimal solutions for most of problems in the solved data set. From the results, the following conclusions can be drawn: 1) DE can be used for effectively solving discrete problems with integer variables, and this can enable the use of research work previously proposed for DE in the continuous domain, to be used for solving complex discrete problems. 2) Using an ensemble of DE mutations can enhance the performance of DE, and at the same time result higher standard deviation of the AE% values, which means more variations in the obtained solutions. This emphasizes that the ensemble mutation can explore the search space more effectively than a single mutation, and hence discover more solutions in a wider search space. In future work, the proposed DE will be tested for solving more complex discrete problems, such as traveling thief problems (TTPs). Moreover, a version of DE, which can be used for solving generic combinatorial optimization problems, will be designed and applied for solving some well-known optimization problems with a discrete domain. [1] Ahmed, Z. H., 2010. Genetic algorithm for the traveling salesman problem using sequential constructive crossover operator. International Journal of Biometrics & Bioinformatics (IJBB) 3 (6), 96–105. [2] Alba, E., Troya, J. M., 1999. A survey of parallel distributed genetic algorithms. Complexity 4 (4), 31–52. [3] Albayrak, M., Allahverdi, N., 2011. Development a new mutation operator to solve the traveling salesman problem by aid of genetic algorithms. Expert Systems with Applications 38 (3), 1313–1320. [4] Ali, I. M., Elsayed, S. M., Ray, T., Sarker, R. A., 2015. Memetic algorithm for solving resource constrained project scheduling problems. In: Evolutionary Computation (CEC), 2015 IEEE Congress on. IEEE, pp. 2761–2767. [5] Ali, I. M., Essam, D., Kasmarik, K., 2018. An efficient differential evolution algorithm for solving 0–1 knapsack problems. In: 2018 IEEE Congress on Evolutionary Computation (CEC). IEEE, pp. 1–8. [6] Ali, I. M., Essam, D., Kasmarik, K., 2019. A novel differential evolution mapping technique for generic combinatorial optimization problems. Applied Soft Computing 80, 297–309. [7] Arabas, J., Bartnik, L., Opara, K., 2011. Dmea- an algorithm that combines differential mutation with the fitness proportionate selection. In: Differential Evolution (SDE), 2011 IEEE Symposium on. IEEE, pp. 1–8. [8] Awad, N. H., Ali, M. Z., Mallipeddi, R., Suganthan, P. N., 2019. An efficient differential evolution algorithm for stochastic opf based active– reactive power dispatch problem considering renewable generators. Applied Soft Computing 76, 445–458. 36

[9] Baraldi, P., Bonfanti, G., Zio, E., 2018. Differential evolution-based multiobjective optimization for the definition of a health indicator for fault diagnostics and prognostics. Mechanical Systems and Signal Processing 102, 382–400. [10] Bean, J. C., 1994. Genetic algorithms and random keys for sequencing and optimization. ORSA journal on computing 6 (2), 154–160. [11] Cao, B., Zhao, J., Yang, P., Lv, Z., Liu, X., Min, G., 2018. 3-d multiobjective deployment of an industrial wireless sensor network for maritime applications utilizing a distributed parallel algorithm. IEEE Transactions on Industrial Informatics 14 (12), 5487–5495. [12] Cekmez, U., Ozsiginan, M., Sahingoz, O. K., 2013. Adapting the ga approach to solve traveling salesman problems on cuda architecture. In: Computational Intelligence and Informatics (CINTI), 2013 IEEE 14th International Symposium on. IEEE, pp. 423–428. [13] Chakravarthi, M., Chandramohan, B., 2019. Estimation of sampling time offsets in an n-channel time-interleaved adc network using differential evolution algorithm and correction using fractional delay filters. In: Machine Intelligence and Signal Analysis. Springer, pp. 267–278. [14] Chen, W.-N., Zhang, J., Chung, H. S., Zhong, W.-L., Wu, W.-G., Shi, Y.H., 2010. A novel set-based particle swarm optimization method for discrete optimization problems. IEEE Transactions on evolutionary computation 14 (2), 278–300. [15] Cheng, J., Zhang, G., Li, Z., Li, Y., 2012. Multi-objective ant colony optimization based on decomposition for bi-objective traveling salesman problems. Soft Computing 16 (4), 597–614. [16] Choong, S. S., Wong, L.-P., Lim, C. P., 2018. An artificial bee colony algorithm with a modified choice function for the traveling salesman problem. Swarm and Evolutionary Computation. [17] Corder, G. W., Foreman, D. I., 2009. Comparing two related samples: The wilcoxon signed ranks test. Nonparametric Statistics for Non-Statisticians. John Wiley & Sons, Inc, 38–56. [18] Cubukcuoglu, C., Ekici, B., Tasgetiren, M. F., Sariyildiz, S., 2019. Optimus: Self-adaptive differential evolution with ensemble of mutation strategies for grasshopper algorithmic modeling. Algorithms 12 (7), 141. [19] Das, S., Suganthan, P. N., 2011. Differential evolution: A survey of the state-of-the-art. IEEE transactions on evolutionary computation 15 (1), 4–31.

37

[20] Deng, Y., Liu, Y., Zhou, D., 2015. An improved genetic algorithm with initial population strategy for symmetric tsp. Mathematical Problems in Engineering 2015, 6 pages. [21] Dorigo, M., Gambardella, L. M., 1997. Ant colony system: a cooperative learning approach to the traveling salesman problem. IEEE Transactions on evolutionary computation 1 (1), 53–66. [22] Du, W., Leung, S. Y. S., Tang, Y., Vasilakos, A. V., 2016. Differential evolution with event-triggered impulsive control. IEEE transactions on cybernetics 47 (1), 244–257. [23] Ezugwu, A. E.-S., Adewumi, A. O., 2017. Discrete symbiotic organisms search algorithm for travelling salesman problem. Expert Systems with Applications 87, 70–78. [24] Fan, H.-Y., Lampinen, J., 2003. A trigonometric mutation operation to differential evolution. Journal of global optimization 27 (1), 105–129. [25] Fleming, P. J., Purshouse, R. C., 2002. Evolutionary algorithms in control systems engineering: a survey. Control engineering practice 10 (11), 1223– 1241. [26] Geng, X., Chen, Z., Yang, W., Shi, D., Zhao, K., 2011. Solving the traveling salesman problem based on an adaptive simulated annealing algorithm with greedy search. Applied Soft Computing 11 (4), 3680–3689. [27] Goldberg, D. E., Holland, J. H., 1988. Genetic algorithms and machine learning. Machine learning 3 (2), 95–99. [28] Gülcü, Ş., Mahi, M., Baykan, Ö. K., Kodaz, H., 2018. A parallel cooperative hybrid method based on ant colony optimization and 3-opt algorithm for solving traveling salesman problem. Soft Computing 22 (5), 1669–1685. [29] Gupta, I. K., Shakil, S., Shakil, S., 2019. A hybrid ga-pso algorithm to solve traveling salesman problem. In: Computational Intelligence: Theories, Applications and Future Directions-Volume I. Springer, pp. 453–462. [30] Han, Y.-Y., Gong, D., Sun, X., 2015. A discrete artificial bee colony algorithm incorporating differential evolution for the flow-shop scheduling problem with blocking. Engineering Optimization 47 (7), 927–946. [31] Hancer, E., Xue, B., Zhang, M., 2018. Differential evolution for filter feature selection based on information theory and feature ranking. KnowledgeBased Systems 140, 103–119. [32] Jünger, M., Reinelt, G., Rinaldi, G., 1995. The traveling salesman problem. Handbooks in operations research and management science 7, 225–330.

38

[33] Khan, I., Pal, S., Maiti, M. K., 2018. A modified particle swarm optimization algorithm for solving traveling salesman problem with imprecise cost matrix. In: 2018 4th International Conference on Recent Advances in Information Technology (RAIT). IEEE, pp. 1–8. [34] Laporte, G., 1992. The traveling salesman problem: An overview of exact and approximate algorithms. European Journal of Operational Research 59 (2), 231–247. [35] Lin, S., Kernighan, B. W., 1973. An effective heuristic algorithm for the traveling-salesman problem. Operations research 21 (2), 498–516. [36] Liu, B., Wang, L., Jin, Y.-H., 2007. An effective pso-based memetic algorithm for flow shop scheduling. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 37 (1), 18–27. [37] Lo, K. M., Yi, W. Y., Wong, P.-K., Leung, K.-S., Leung, Y., Mak, S.-T., 2018. A genetic algorithm with new local operators for multiple traveling salesman problems. Int. J. Comput. Intell. Syst. 11 (1), 692–705. [38] Mahi, M., Baykan, Ö. K., Kodaz, H., 2015. A new hybrid method based on particle swarm optimization, ant colony optimization and 3-opt algorithms for traveling salesman problem. Applied Soft Computing 30, 484–490. [39] Makuchowski, M., 2018. Effective algorithm of simulated annealing for the symmetric traveling salesman problem. In: International Conference on Dependability and Complex Systems. Springer, pp. 348–359. [40] Mi, M., Huifeng, X., Ming, Z., Yu, G., 2010. An improved differential evolution algorithm for tsp problem. In: Intelligent Computation Technology and Automation (ICICTA), 2010 International Conference on. Vol. 1. IEEE, pp. 544–547. [41] Mohamed, A. W., 2018. A novel differential evolution algorithm for solving constrained engineering optimization problems. Journal of Intelligent Manufacturing 29 (3), 659–692. [42] Myszkowski, P. B., Olech, Ł. P., Laszczyk, M., Skowroński, M. E., 2018. Hybrid differential evolution and greedy algorithm (degr) for solving multiskill resource-constrained project scheduling problem. Applied Soft Computing 62, 1–14. [43] Nahas, N., Nourelfath, M., 2018. Joint optimization of maintenance, buffers and machines in manufacturing lines. Engineering Optimization 50 (1), 37– 54. [44] Nguyen, T. T., Vu Quynh, N., Duong, M. Q., Van Dai, L., et al., 2018. Modified differential evolution algorithm: A novel approach to optimize the operation of hydrothermal power systems while considering the different constraints and valve point loading effects. Energies 11 (3), 540. 39

[45] Nouiri, M., Bekrar, A., Jemai, A., Niar, S., Ammari, A. C., 2018. An effective and distributed particle swarm optimization algorithm for flexible job-shop scheduling problem. Journal of Intelligent Manufacturing 29 (3), 603–615. [46] Opara, K. R., Arabas, J., 2019. Differential evolution: A survey of theoretical analyses. Swarm and evolutionary computation 44, 546–558. [47] Osaba, E., Yang, X.-S., Diaz, F., Lopez-Garcia, P., Carballedo, R., 2016. An improved discrete bat algorithm for symmetric and asymmetric traveling salesman problems. Engineering Applications of Artificial Intelligence 48, 59–71. [48] Papadimitriou, C. H., 1977. The euclidean travelling salesman problem is np-complete. Theoretical computer science 4 (3), 237–244. [49] Podder, T., Bhattachya, D., Chakraborty, S., 2019. Adaptive differential evolution with intersect mutation and repaired crossover rate. International Journal of Computational Intelligence & IoT 2 (1). [50] Price, K., Storn, R. M., Lampinen, J. A., 2006. Differential evolution: a practical approach to global optimization. Springer Science & Business Media. [51] Reddy, S. S., 2019. Optimal power flow using hybrid differential evolution and harmony search algorithm. International Journal of Machine Learning and Cybernetics 10 (5), 1077–1091. [52] Rego, C., Gamboa, D., Glover, F., Osterman, C., 2011. Traveling salesman problem heuristics: Leading methods, implementations and latest advances. European Journal of Operational Research 211 (3), 427–441. [53] Reinelt, G., 1991. Tsplib-a traveling salesman problem library. ORSA journal on computing 3 (4), 376–384. [54] Richter, D., Goldengorin, B., Jager, G., Molitor, P., 2007. Improving the efficiency of helsgaun lin-kernighan heuristic for the symmetric tsp. In: Workshop on Combinatorial and Algorithmic Aspects of Networking. Springer, pp. 99–111. [55] Storn, R., Price, K., 1997. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization 11 (4), 341–359. [56] Tanabe, R., Fukunaga, A., 2014. Reevaluating exponential crossover in differential evolution. In: International Conference on Parallel Problem Solving from Nature. Springer, pp. 201–210. [57] Wang, X., Xu, G., 2011. Hybrid differential evolution algorithm for traveling salesman problem. Procedia Engineering 15, 2716–2720. 40

[58] Wei, H., Hao, Z., Huang, H., Li, G., Chen, Q., 2016. A real adjacency matrix-coded differential evolution algorithm for traveling salesman problems. In: Bio-Inspired Computing-Theories and Applications. Springer, pp. 135–140. [59] Wu, G., Mallipeddi, R., Suganthan, P. N., 2019. Ensemble strategies for population-based optimization algorithms–a survey. Swarm and evolutionary computation 44, 695–711. [60] Wu, Y., Liu, X., Guan, W., Chen, B., Chen, X., Xie, C., 2018. High-speed 3d indoor localization system based on visible light communication using differential evolution algorithm. Optics Communications 424, 177–189. [61] Xu, X., Rong, H., Trovati, M., Liptrott, M., Bessis, N., 2018. Cs-pso: chaotic particle swarm optimization algorithm for solving combinatorial optimization problems. Soft Computing 22 (3), 783–795. [62] Yip, P. P., Pao, Y.-H., 1995. Combinatorial optimization with use of guided evolutionary simulated annealing. IEEE Transactions on neural networks 6 (2), 290–295. [63] Yu, X., Chen, W.-N., Gu, T., Zhang, H., Yuan, H., Kwong, S., Zhang, J., 2018. Set-based discrete particle swarm optimization based on decomposition for permutation-based multiobjective combinatorial optimization problems. IEEE transactions on cybernetics 48 (7), 2139–2153. [64] Zhang, C., Yang, H., Li, J., July 2018. Hybrid discrete differential evolution algorithm for motor train-sets scheduling. In: 2018 37th Chinese Control Conference (CCC). pp. 3245–3248. [65] Zhong, Y., Lin, J., Wang, L., Zhang, H., 2017. Hybrid discrete artificial bee colony algorithm with threshold acceptance criterion for traveling salesman problem. Information Sciences 421, 70–84. [66] Zhong, Y., Lin, J., Wang, L., Zhang, H., 2018. Discrete comprehensive learning particle swarm optimization algorithm with metropolis acceptance criterion for traveling salesman problem. [67] Zhong, Y., Wang, L., Lin, M., Zhang, H., 2019. Discrete pigeon-inspired optimization algorithm with metropolis acceptance criterion for large-scale traveling salesman problem. Swarm and Evolutionary Computation 48, 134–144. [68] Zhou, B.-h., Hu, L.-m., Zhong, Z.-y., 2018. A hybrid differential evolution algorithm with estimation of distribution algorithm for reentrant hybrid flow shop scheduling problem. Neural Computing and Applications 30 (1), 193–209. [69] Zhou, Y.-Z., Yi, W.-C., Gao, L., Li, X.-Y., 2017. Adaptive differential evolution with sorting crossover rate for continuous optimization problems. IEEE transactions on cybernetics 47 (9), 2742–2753. 41